C++ RCP::L2NormOfSolutionGlobal方法代码示例

//.........这里部分代码省略.........
  cellFace3.push_back(make_pair(mesh->getElement(9 ), 3));
  cellFace3.push_back(make_pair(mesh->getElement(8 ), 3));
  cellFace3.push_back(make_pair(mesh->getElement(19), 3));
  cellFace3.push_back(make_pair(mesh->getElement(18), 3));
  cellFace4.push_back(make_pair(mesh->getElement(9 ), 1));
  cellFace4.push_back(make_pair(mesh->getElement(8 ), 1));
  cellFace4.push_back(make_pair(mesh->getElement(19), 1));
  cellFace4.push_back(make_pair(mesh->getElement(18), 1));

  // // for loading refinement history
  // if (replayFile.length() > 0) {
  //   RefinementHistory refHistory;
  //   replayFile = replayFile;
  //   refHistory.loadFromFile(replayFile);
  //   refHistory.playback(mesh);
  //   int numElems = mesh->numActiveElements();
  //   if (commRank==0){
  //     double minSideLength = meshInfo.getMinCellSideLength() ;
  //     cout << "after replay, num elems = " << numElems << " and min side length = " << minSideLength << endl;
  //   }
  // }

  for (int i = 0; i < uniformRefinements; i++)
    refinementStrategy.hRefineUniformly(mesh);

  double nonlinearRelativeEnergyTolerance = 1e-5; // used to determine convergence of the nonlinear solution
  for (int refIndex=0; refIndex<=numRefs; refIndex++)
  {
    double L2Update = 1e10;
    int iterCount = 0;
    while (L2Update > nonlinearRelativeEnergyTolerance && iterCount < maxNewtonIterations)
    {
      solution->solve(false);
      double u1L2Update = solution->L2NormOfSolutionGlobal(u1->ID());
      double u2L2Update = solution->L2NormOfSolutionGlobal(u2->ID());
      L2Update = sqrt(u1L2Update*u1L2Update + u2L2Update*u2L2Update);
      double energy_error = solution->energyErrorTotal();

      // Check local conservation
      if (commRank == 0)
      {
        FunctionPtr n = Function::normal();
        FunctionPtr u1_prev = Function::solution(u1hat, solution);
        FunctionPtr u2_prev = Function::solution(u2hat, solution);
        FunctionPtr flux = u1_prev*n->x() + u2_prev*n->y();
        Teuchos::Tuple<double, 3> fluxImbalances = checkConservation(flux, zero, mesh);
        // cout << "Mass flux: Largest Local = " << fluxImbalances[0]
        //   << ", Global = " << fluxImbalances[1] << ", Sum Abs = " << fluxImbalances[2] << endl;

        errOut << mesh->numGlobalDofs() << " " << energy_error << " "
               << fluxImbalances[0] << " " << fluxImbalances[1] << " " << fluxImbalances[2] << endl;

        double massFlux0 = computeFluxOverElementSides(u1_prev, mesh, cellFace0);
        double massFlux1 = computeFluxOverElementSides(u1_prev, mesh, cellFace1);
        double massFlux2 = computeFluxOverElementSides(u1_prev, mesh, cellFace2);
        double massFlux3 = computeFluxOverElementSides(u1_prev, mesh, cellFace3);
        double massFlux4 = computeFluxOverElementSides(u1_prev, mesh, cellFace4);
        fluxOut << massFlux0 << " " << massFlux1 << " " << massFlux2 << " " << massFlux3 << " " << massFlux4 << " " << endl;
        cout << "Total mass flux = " << massFlux0 << " " << massFlux1 << " " << massFlux2 << " " << massFlux3 << " " << massFlux4 << " " << endl;

        // if (saveFile.length() > 0) {
        //   std::ostringstream oss;
        //   oss << string(saveFile) << refIndex ;
        //   cout << "on refinement " << refIndex << " saving mesh file to " << oss.str() << endl;
        //   refHistory->saveToFile(oss.str());
        // }

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

  ////////////////////   DEFINE INNER PRODUCT(S)   ///////////////////////
  IPPtr ip = Teuchos::rcp(new IP);
  if (norm == 0)
  {
    ip = bf->graphNorm();
  }
  else if (norm == 1)
  {
    // ip = bf->l2Norm();
  }

  ////////////////////   CREATE BCs   ///////////////////////
  Teuchos::RCP<BCEasy> bc = Teuchos::rcp( new BCEasy );
  // Teuchos::RCP<PenaltyConstraints> pc = Teuchos::rcp( new PenaltyConstraints );
  SpatialFilterPtr left = Teuchos::rcp( new ConstantXBoundary(-0.5) );
  SpatialFilterPtr right = Teuchos::rcp( new ConstantXBoundary(1) );
  SpatialFilterPtr top = Teuchos::rcp( new ConstantYBoundary(-0.5) );
  SpatialFilterPtr bottom = Teuchos::rcp( new ConstantYBoundary(1.5) );
  bc->addDirichlet(u1hat, left, u1Exact);
  bc->addDirichlet(u2hat, left, u2Exact);
  bc->addDirichlet(u1hat, right, u1Exact);
  bc->addDirichlet(u2hat, right, u2Exact);
  bc->addDirichlet(u1hat, top, u1Exact);
  bc->addDirichlet(u2hat, top, u2Exact);
  bc->addDirichlet(u1hat, bottom, u1Exact);
  bc->addDirichlet(u2hat, bottom, u2Exact);

  // zero mean constraint on pressure
  bc->addZeroMeanConstraint(p);

  // pc->addConstraint(u1hat*u2hat-t1hat == zero, top);
  // pc->addConstraint(u2hat*u2hat-t2hat == zero, top);

  Teuchos::RCP<Solution> solution = Teuchos::rcp( new Solution(mesh, bc, rhs, ip) );
  // solution->setFilter(pc);

  // if (enforceLocalConservation) {
  //   solution->lagrangeConstraints()->addConstraint(u1hat->times_normal_x() + u2hat->times_normal_y() == zero);
  // }

  // ==================== Register Solutions ==========================
  mesh->registerSolution(solution);
  mesh->registerSolution(backgroundFlow);

  // Teuchos::RCP< RefinementHistory > refHistory = Teuchos::rcp( new RefinementHistory );
  // mesh->registerObserver(refHistory);

  ////////////////////   SOLVE & REFINE   ///////////////////////
  double energyThreshold = 0.2; // for mesh refinements
  RefinementStrategy refinementStrategy( solution, energyThreshold );
  VTKExporter exporter(backgroundFlow, mesh, varFactory);
  stringstream outfile;
  outfile << "kovasznay" << "_" << 0;
  exporter.exportSolution(outfile.str());

  double nonlinearRelativeEnergyTolerance = 1e-5; // used to determine convergence of the nonlinear solution
  for (int refIndex=0; refIndex<=numRefs; refIndex++)
  {
    double L2Update = 1e10;
    int iterCount = 0;
    while (L2Update > nonlinearRelativeEnergyTolerance && iterCount < maxNewtonIterations)
    {
      solution->solve(false);
      double u1L2Update = solution->L2NormOfSolutionGlobal(u1->ID());
      double u2L2Update = solution->L2NormOfSolutionGlobal(u2->ID());
      L2Update = sqrt(u1L2Update*u1L2Update + u2L2Update*u2L2Update);

      // Check local conservation
      if (commRank == 0)
      {
        cout << "L2 Norm of Update = " << L2Update << endl;

        // if (saveFile.length() > 0) {
        //   std::ostringstream oss;
        //   oss << string(saveFile) << refIndex ;
        //   cout << "on refinement " << refIndex << " saving mesh file to " << oss.str() << endl;
        //   refHistory->saveToFile(oss.str());
        // }
      }

      // line search algorithm
      double alpha = 1.0;
      backgroundFlow->addSolution(solution, alpha);
      iterCount++;
    }

    if (commRank == 0)
    {
      stringstream outfile;
      outfile << "kovasznay" << "_" << refIndex+1;
      exporter.exportSolution(outfile.str());
    }

    if (refIndex < numRefs)
      refinementStrategy.refine(commRank==0); // print to console on commRank 0
  }

  return 0;
}

//.........这里部分代码省略.........
  inflowBC->addDirichlet(fhat, bottomBoundary, -u0 );
  inflowBC->addDirichlet(fhat, leftBoundary, leftVal );
  inflowBC->addDirichlet(fhat, rightBoundary, rightVal );

  ////////////////////////////////////////////////////////////////////
  // CREATE SOLUTION OBJECT
  ////////////////////////////////////////////////////////////////////
  Teuchos::RCP<Solution> solution = Teuchos::rcp(new Solution(mesh, inflowBC, rhs, ip));
  mesh->registerSolution(solution);

  if (enforceLocalConservation)
  {
    FunctionPtr zero = Teuchos::rcp( new ConstantScalarFunction(0.0) );
    solution->lagrangeConstraints()->addConstraint(fhat == zero);
  }

  ////////////////////////////////////////////////////////////////////
  // 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++)