C++ SolutionPtr::writeToVTK方法代码示例

//.........这里部分代码省略.........
    //      cout << "fakeRHSIntegrals:\n" << fakeRHSIntegrals;
    for (int i=0; i<elems.size(); i++)
    {
      int cellID = cellIDs[i];
      // pick out the ones for testOne:
      massFluxIntegral[cellID] = fakeRHSIntegrals(i,testOneIndex);
    }
    // find the largest:
    for (int i=0; i<elems.size(); i++)
    {
      int cellID = cellIDs[i];
      maxMassFluxIntegral = max(abs(massFluxIntegral[cellID]), maxMassFluxIntegral);
    }
    for (int i=0; i<elems.size(); i++)
    {
      int cellID = cellIDs[i];
      maxCellMeasure = max(maxCellMeasure,cellMeasures(i));
      minCellMeasure = min(minCellMeasure,cellMeasures(i));
      maxMassFluxIntegral = max(abs(massFluxIntegral[cellID]), maxMassFluxIntegral);
      totalMassFlux += massFluxIntegral[cellID];
      totalAbsMassFlux += abs( massFluxIntegral[cellID] );
    }
  }
  if (rank==0)
  {
    cout << "largest mass flux: " << maxMassFluxIntegral << endl;
    cout << "total mass flux: " << totalMassFlux << endl;
    cout << "sum of mass flux absolute value: " << totalAbsMassFlux << endl;
    cout << "largest h: " << sqrt(maxCellMeasure) << endl;
    cout << "smallest h: " << sqrt(minCellMeasure) << endl;
    cout << "ratio of largest / smallest h: " << sqrt(maxCellMeasure) / sqrt(minCellMeasure) << endl;
  }
  if (rank == 0)
  {
    cout << "phi ID: " << phi->ID() << endl;
    cout << "psi1 ID: " << psi_1->ID() << endl;
    cout << "psi2 ID: " << psi_2->ID() << endl;

    cout << "streamMesh has " << streamMesh->numActiveElements() << " elements.\n";
    cout << "solving for approximate stream function...\n";
  }

  streamSolution->solve(useMumps);
  energyErrorTotal = streamSolution->energyErrorTotal();
  if (rank == 0)
  {
    cout << "...solved.\n";
    cout << "Stream mesh has energy error: " << energyErrorTotal << endl;
  }

  if (rank==0)
  {
    solution->writeToVTK("nsCavitySoln.vtk");
    if (! meshHasTriangles )
    {
      massFlux->writeBoundaryValuesToMATLABFile(solution->mesh(), "massFlux.dat");
      u_mag->writeValuesToMATLABFile(solution->mesh(), "u_mag.m");
      u_div->writeValuesToMATLABFile(solution->mesh(), "u_div.m");
      solution->writeFieldsToFile(u1->ID(), "u1.m");
      solution->writeFluxesToFile(u1hat->ID(), "u1_hat.dat");
      solution->writeFieldsToFile(u2->ID(), "u2.m");
      solution->writeFluxesToFile(u2hat->ID(), "u2_hat.dat");
      solution->writeFieldsToFile(p->ID(), "p.m");
      streamSolution->writeFieldsToFile(phi->ID(), "phi.m");

      streamSolution->writeFluxesToFile(phi_hat->ID(), "phi_hat.dat");
      streamSolution->writeFieldsToFile(psi_1->ID(), "psi1.m");
      streamSolution->writeFieldsToFile(psi_2->ID(), "psi2.m");
      vorticity->writeValuesToMATLABFile(streamMesh, "vorticity.m");

      FunctionPtr ten = Teuchos::rcp( new ConstantScalarFunction(10) );
      ten->writeBoundaryValuesToMATLABFile(solution->mesh(), "skeleton.dat");
      cout << "wrote files: u_mag.m, u_div.m, u1.m, u1_hat.dat, u2.m, u2_hat.dat, p.m, phi.m, vorticity.m.\n";
    }
    else
    {
      solution->writeToFile(u1->ID(), "u1.dat");
      solution->writeToFile(u2->ID(), "u2.dat");
      solution->writeToFile(u2->ID(), "p.dat");
      cout << "wrote files: u1.dat, u2.dat, p.dat\n";
    }

    FieldContainer<double> points = pointGrid(0, 1, 0, 1, 100);
    FieldContainer<double> pointData = solutionData(points, streamSolution, phi);
    GnuPlotUtil::writeXYPoints("phi_patch_navierStokes_cavity.dat", pointData);
    set<double> patchContourLevels = diagonalContourLevels(pointData,1);
    vector<string> patchDataPath;
    patchDataPath.push_back("phi_patch_navierStokes_cavity.dat");
    GnuPlotUtil::writeContourPlotScript(patchContourLevels, patchDataPath, "lidCavityNavierStokes.p");

    GnuPlotUtil::writeExactMeshSkeleton("lid_navierStokes_continuation_adaptive", mesh, 2);

    writePatchValues(0, 1, 0, 1, streamSolution, phi, "phi_patch.m");
    writePatchValues(0, .1, 0, .1, streamSolution, phi, "phi_patch_detail.m");
    writePatchValues(0, .01, 0, .01, streamSolution, phi, "phi_patch_minute_detail.m");
    writePatchValues(0, .001, 0, .001, streamSolution, phi, "phi_patch_minute_minute_detail.m");
  }

  return 0;
}

//.........这里部分代码省略.........

  ////////////////////////////////////////////////////////////////////
  // HESSIAN BIT + CHECKS ON GRADIENT + HESSIAN
  ////////////////////////////////////////////////////////////////////

  VarFactory hessianVars = varFactory.getBubnovFactory(VarFactory::BUBNOV_TRIAL);
  VarPtr du = hessianVars.test(u->ID());
  BFPtr hessianBF = Teuchos::rcp( new BF(hessianVars) ); // initialize bilinear form

  FunctionPtr du_current  = Teuchos::rcp( new PreviousSolutionFunction(solution, u) );

  FunctionPtr fnhat = Teuchos::rcp(new PreviousSolutionFunction(solution,fn));
  LinearTermPtr residual = Teuchos::rcp(new LinearTerm);// residual
  residual->addTerm(fnhat*v,true);
  residual->addTerm( - (e1 * (u_prev_squared_div2) + e2 * (u_prev)) * v->grad(),true);

  LinearTermPtr Bdu = Teuchos::rcp(new LinearTerm);// residual
  Bdu->addTerm( - du_current*(beta*v->grad()));

  Teuchos::RCP<RieszRep> riesz = Teuchos::rcp(new RieszRep(mesh, ip, residual));
  Teuchos::RCP<RieszRep> duRiesz = Teuchos::rcp(new RieszRep(mesh, ip, Bdu));
  riesz->computeRieszRep();
  FunctionPtr e_v = Teuchos::rcp(new RepFunction(v,riesz));
  e_v->writeValuesToMATLABFile(mesh, "e_v.m");
  FunctionPtr posErrPart = Teuchos::rcp(new PositivePart(e_v->dx()));
  hessianBF->addTerm(e_v->dx()*u,du);
  //  hessianBF->addTerm(posErrPart*u,du);
  Teuchos::RCP<HessianFilter> hessianFilter = Teuchos::rcp(new HessianFilter(hessianBF));

  if (useHessian)
  {
    solution->setWriteMatrixToFile(true,"hessianStiffness.dat");
  }
  else
  {
    solution->setWriteMatrixToFile(true,"stiffness.dat");
  }

  Teuchos::RCP< LineSearchStep > LS_Step = Teuchos::rcp(new LineSearchStep(riesz));

  double NL_residual = 9e99;
  for (int i = 0; i<numSteps; i++)
  {
    solution->solve(false); // do one solve to initialize things...
    double stepLength = 1.0;
    stepLength = LS_Step->stepSize(backgroundFlow,solution, NL_residual);
    if (useHessian)
    {
      solution->setFilter(hessianFilter);
    }
    backgroundFlow->addSolution(solution,stepLength);
    NL_residual = LS_Step->getNLResidual();
    if (rank==0)
    {
      cout << "NL residual after adding = " << NL_residual << " with step size " << stepLength << endl;
    }

    int numGlobalDofs = mesh->numGlobalDofs();
    double fd_gradient;
    for (int dofIndex = 0; dofIndex<numGlobalDofs; dofIndex++)
    {

      TestingUtilities::initializeSolnCoeffs(solnPerturbation);
      TestingUtilities::setSolnCoeffForGlobalDofIndex(solnPerturbation,1.0,dofIndex);

      fd_gradient = FiniteDifferenceUtilities::finiteDifferenceGradient(mesh, riesz, backgroundFlow, dofIndex);

      // CHECK GRADIENT
      LinearTermPtr b_u =  bf->testFunctional(solnPerturbation);
      map<int,FunctionPtr> NL_err_rep_map;

      NL_err_rep_map[v->ID()] = Teuchos::rcp(new RepFunction(v,riesz));
      FunctionPtr gradient = b_u->evaluate(NL_err_rep_map, TestingUtilities::isFluxOrTraceDof(mesh,dofIndex)); // use boundary part only if flux or trace
      double grad;
      if (TestingUtilities::isFluxOrTraceDof(mesh,dofIndex))
      {
        grad = gradient->integralOfJump(mesh,10);
      }
      else
      {
        grad = gradient->integrate(mesh,10);
      }
      double fdgrad = fd_gradient;
      double diff = grad-fdgrad;
      if (abs(diff)>1e-6 && i>0)
      {
        cout << "Found difference of " << diff << ", " << " with fd val = " << fdgrad << " and gradient = " << grad << " in dof " << dofIndex << ", isTraceDof = " << TestingUtilities::isFluxOrTraceDof(mesh,dofIndex) << endl;
      }
    }
  }

  if (rank==0)
  {
    backgroundFlow->writeToVTK("BurgersTest.vtu",min(H1Order+1,4));
    solution->writeFluxesToFile(fn->ID(),"fn.dat");
    cout << "wrote solution files" << endl;
  }

  return 0;
}

//.........这里部分代码省略.........
    mesh->registerSolution(solution);

    ////////////////////////////////////////////////////////////////////
    // WARNING: UNFINISHED HESSIAN BIT
    ////////////////////////////////////////////////////////////////////
    VarFactory hessianVars = varFactory.getBubnovFactory(VarFactory::BUBNOV_TRIAL);
    VarPtr du = hessianVars.test(u->ID());
    BFPtr hessianBF = Teuchos::rcp( new BF(hessianVars) ); // initialize bilinear form
    //  FunctionPtr e_v = Function::constant(1.0); // dummy error rep function for now - should do nothing

    FunctionPtr u_current  = Teuchos::rcp( new PreviousSolutionFunction(solution, u) );

    FunctionPtr sig1_prev = Teuchos::rcp( new PreviousSolutionFunction(solution, sigma1) );
    FunctionPtr sig2_prev = Teuchos::rcp( new PreviousSolutionFunction(solution, sigma2) );
    FunctionPtr sig_prev = (e1*sig1_prev + e2*sig2_prev);
    FunctionPtr fnhat = Teuchos::rcp(new PreviousSolutionFunction(solution,beta_n_u_minus_sigma_hat));
    FunctionPtr uhat_prev = Teuchos::rcp(new PreviousSolutionFunction(solution,uhat));
    LinearTermPtr residual = Teuchos::rcp(new LinearTerm);// residual
    residual->addTerm(fnhat*v - (e1 * (u_prev_squared_div2 - sig1_prev) + e2 * (u_prev - sig2_prev)) * v->grad());
    residual->addTerm((1/epsilon)*sig_prev * tau + u_prev * tau->div() - uhat_prev*tau->dot_normal());

    LinearTermPtr Bdu = Teuchos::rcp(new LinearTerm);// residual
    Bdu->addTerm( u_current*tau->div() - u_current*(beta*v->grad()));

    Teuchos::RCP<RieszRep> riesz = Teuchos::rcp(new RieszRep(mesh, ip, residual));
    Teuchos::RCP<RieszRep> duRiesz = Teuchos::rcp(new RieszRep(mesh, ip, Bdu));
    riesz->computeRieszRep();
    FunctionPtr e_v = Teuchos::rcp(new RepFunction(v,riesz));
    e_v->writeValuesToMATLABFile(mesh, "e_v.m");
    FunctionPtr posErrPart = Teuchos::rcp(new PositivePart(e_v->dx()));
    hessianBF->addTerm(e_v->dx()*u,du);
    //  hessianBF->addTerm(posErrPart*u,du);
    Teuchos::RCP<HessianFilter> hessianFilter = Teuchos::rcp(new HessianFilter(hessianBF));

    if (useHessian)
    {
        solution->setWriteMatrixToFile(true,"hessianStiffness.dat");
    }
    else
    {
        solution->setWriteMatrixToFile(true,"stiffness.dat");
    }

    Teuchos::RCP< LineSearchStep > LS_Step = Teuchos::rcp(new LineSearchStep(riesz));
    ofstream out;
    out.open("Burgers.txt");
    double NL_residual = 9e99;
    for (int i = 0; i<numSteps; i++)
    {
        solution->solve(false); // do one solve to initialize things...
        double stepLength = 1.0;
        stepLength = LS_Step->stepSize(backgroundFlow,solution, NL_residual);
        if (useHessian)
        {
            solution->setFilter(hessianFilter);
        }
        backgroundFlow->addSolution(solution,stepLength);
        NL_residual = LS_Step->getNLResidual();
        if (rank==0)
        {
            cout << "NL residual after adding = " << NL_residual << " with step size " << stepLength << endl;
            out << NL_residual << endl; // saves initial NL error
        }
    }
    out.close();


    ////////////////////////////////////////////////////////////////////
    // DEFINE REFINEMENT STRATEGY
    ////////////////////////////////////////////////////////////////////
    Teuchos::RCP<RefinementStrategy> refinementStrategy;
    refinementStrategy = Teuchos::rcp(new RefinementStrategy(solution,energyThreshold));

    int numRefs = 0;

    Teuchos::RCP<NonlinearStepSize> stepSize = Teuchos::rcp(new NonlinearStepSize(nonlinearStepSize));
    Teuchos::RCP<NonlinearSolveStrategy> solveStrategy;
    solveStrategy = Teuchos::rcp( new NonlinearSolveStrategy(backgroundFlow, solution, stepSize,
                                  nonlinearRelativeEnergyTolerance));

    ////////////////////////////////////////////////////////////////////
    // SOLVE
    ////////////////////////////////////////////////////////////////////

    for (int refIndex=0; refIndex<numRefs; refIndex++)
    {
        solveStrategy->solve(rank==0);       // print to console on rank 0
        refinementStrategy->refine(rank==0); // print to console on rank 0
    }
    //  solveStrategy->solve(rank==0);

    if (rank==0)
    {
        backgroundFlow->writeToVTK("Burgers.vtu",min(H1Order+1,4));
        solution->writeFluxesToFile(uhat->ID(), "burgers.dat");
        cout << "wrote solution files" << endl;
    }

    return 0;
}

//.........这里部分代码省略.........
  }

  ////////////////////////////////////////////////////////////////////
  // DEFINE REFINEMENT STRATEGY
  ////////////////////////////////////////////////////////////////////
  Teuchos::RCP<RefinementStrategy> refinementStrategy;
  refinementStrategy = Teuchos::rcp(new RefinementStrategy(solution,energyThreshold));

  ////////////////////////////////////////////////////////////////////
  // SOLVE
  ////////////////////////////////////////////////////////////////////

  for (int refIndex=0; refIndex<=numRefs; refIndex++)
  {
    double L2Update = 1e7;
    int iterCount = 0;
    while (L2Update > nonlinearRelativeEnergyTolerance && iterCount < maxNewtonIterations)
    {
      solution->solve();
      L2Update = solution->L2NormOfSolutionGlobal(u->ID());
      cout << "L2 Norm of Update = " << L2Update << endl;
      // backgroundFlow->clear();
      backgroundFlow->addSolution(solution, newtonStepSize);
      iterCount++;
    }
    cout << endl;

    // check conservation
    VarPtr testOne = varFactory.testVar("1", CONSTANT_SCALAR);
    // Create a fake bilinear form for the testing
    BFPtr fakeBF = Teuchos::rcp( new BF(varFactory) );
    // Define our mass flux
    FunctionPtr massFlux = Teuchos::rcp( new PreviousSolutionFunction(solution, fhat) );
    LinearTermPtr massFluxTerm = massFlux * testOne;

    Teuchos::RCP<shards::CellTopology> quadTopoPtr = Teuchos::rcp(new shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >() ));
    DofOrderingFactory dofOrderingFactory(fakeBF);
    int fakeTestOrder = H1Order;
    DofOrderingPtr testOrdering = dofOrderingFactory.testOrdering(fakeTestOrder, *quadTopoPtr);

    int testOneIndex = testOrdering->getDofIndex(testOne->ID(),0);
    vector< ElementTypePtr > elemTypes = mesh->elementTypes(); // global element types
    map<int, double> massFluxIntegral; // cellID -> integral
    double maxMassFluxIntegral = 0.0;
    double totalMassFlux = 0.0;
    double totalAbsMassFlux = 0.0;
    for (vector< ElementTypePtr >::iterator elemTypeIt = elemTypes.begin(); elemTypeIt != elemTypes.end(); elemTypeIt++)
    {
      ElementTypePtr elemType = *elemTypeIt;
      vector< ElementPtr > elems = mesh->elementsOfTypeGlobal(elemType);
      vector<int> cellIDs;
      for (int i=0; i<elems.size(); i++)
      {
        cellIDs.push_back(elems[i]->cellID());
      }
      FieldContainer<double> physicalCellNodes = mesh->physicalCellNodesGlobal(elemType);
      BasisCachePtr basisCache = Teuchos::rcp( new BasisCache(elemType,mesh) );
      basisCache->setPhysicalCellNodes(physicalCellNodes,cellIDs,true); // true: create side caches
      FieldContainer<double> cellMeasures = basisCache->getCellMeasures();
      FieldContainer<double> fakeRHSIntegrals(elems.size(),testOrdering->totalDofs());
      massFluxTerm->integrate(fakeRHSIntegrals,testOrdering,basisCache,true); // true: force side evaluation
      for (int i=0; i<elems.size(); i++)
      {
        int cellID = cellIDs[i];
        // pick out the ones for testOne:
        massFluxIntegral[cellID] = fakeRHSIntegrals(i,testOneIndex);
      }
      // find the largest:
      for (int i=0; i<elems.size(); i++)
      {
        int cellID = cellIDs[i];
        maxMassFluxIntegral = max(abs(massFluxIntegral[cellID]), maxMassFluxIntegral);
      }
      for (int i=0; i<elems.size(); i++)
      {
        int cellID = cellIDs[i];
        maxMassFluxIntegral = max(abs(massFluxIntegral[cellID]), maxMassFluxIntegral);
        totalMassFlux += massFluxIntegral[cellID];
        totalAbsMassFlux += abs( massFluxIntegral[cellID] );
      }
    }
    if (rank==0)
    {
      cout << endl;
      cout << "largest mass flux: " << maxMassFluxIntegral << endl;
      cout << "total mass flux: " << totalMassFlux << endl;
      cout << "sum of mass flux absolute value: " << totalAbsMassFlux << endl;
      cout << endl;

      stringstream outfile;
      outfile << "burgers_" << refIndex;
      backgroundFlow->writeToVTK(outfile.str(), 5);
    }

    if (refIndex < numRefs)
      refinementStrategy->refine(rank==0); // print to console on rank 0
  }

  return 0;
}