C++ IPPtr::evaluate方法代码示例

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

  // v terms:
  confusionBF->addTerm( sigma1, v->dx() );
  confusionBF->addTerm( sigma2, v->dy() );
  confusionBF->addTerm( -u, beta * v->grad() );
  confusionBF->addTerm( beta_n_u_minus_sigma_n, v);

  ////////////////////   DEFINE INNER PRODUCT(S)   ///////////////////////

  // robust test norm
  IPPtr ip = Teuchos::rcp(new IP);

  // choose the mesh-independent norm even though it may have boundary layers
  ip->addTerm(v->grad());
  ip->addTerm(v);
  ip->addTerm(tau);
  ip->addTerm(tau->div());

  ////////////////////   SPECIFY RHS AND HELPFUL FUNCTIONS   ///////////////////////

  FunctionPtr n = Function::normal();
  vector<double> e1,e2;
  e1.push_back(1.0);
  e1.push_back(0.0);
  e2.push_back(0.0);
  e2.push_back(1.0);
  FunctionPtr one = Function::constant(1.0);

  FunctionPtr zero = Function::constant(0.0);
  RHSPtr rhs = RHS::rhs();
  FunctionPtr f = one; // if this is set to zero instead, we pass the test (a clue?)
  rhs->addTerm( f * v );

  ////////////////////   CREATE BCs   ///////////////////////
  BCPtr bc = BC::bc();
  SpatialFilterPtr squareBoundary = Teuchos::rcp( new SquareBoundary );

  bc->addDirichlet(uhat, squareBoundary, one);

  ////////////////////   BUILD MESH   ///////////////////////
  // define nodes for mesh
  int order = 2;
  int H1Order = order+1;
  int pToAdd = 2;

  // create a pointer to a new mesh:
  Teuchos::RCP<Mesh> mesh = MeshUtilities::buildUnitQuadMesh(nCells,confusionBF, H1Order, H1Order+pToAdd);

  ////////////////////   SOLVE & REFINE   ///////////////////////

  Teuchos::RCP<Solution> solution;
  solution = Teuchos::rcp( new Solution(mesh, bc, rhs, ip) );
  solution->solve(false);
  double energyError = solution->energyErrorTotal();

  LinearTermPtr residual = rhs->linearTermCopy();
  residual->addTerm(-confusionBF->testFunctional(solution),true);

//  FunctionPtr uh = Function::solution(uhat,solution);
//  FunctionPtr fn = Function::solution(beta_n_u_minus_sigma_n,solution);
//  FunctionPtr uF = Function::solution(u,solution);
//  FunctionPtr sigma = e1*Function::solution(sigma1,solution)+e2*Function::solution(sigma2,solution);
//  residual->addTerm(- (fn*v - uh*tau->dot_normal()));
//  residual->addTerm(- (uF*(tau->div() - beta*v->grad()) + sigma*((1/eps)*tau + v->grad())));
//  residual->addTerm(-(fn*v - uF*beta*v->grad() + sigma*v->grad())); // just v portion
//  residual->addTerm(uh*tau->dot_normal() - uF*tau->div() - sigma*((1/eps)*tau)); // just tau portion

  Teuchos::RCP<RieszRep> rieszResidual = Teuchos::rcp(new RieszRep(mesh, ip, residual));
  rieszResidual->computeRieszRep();
  double energyErrorLT = rieszResidual->getNorm();

  int cubEnrich = 0;
  bool testVsTest = true;
  FunctionPtr e_v = RieszRep::repFunction(v,rieszResidual);
  FunctionPtr e_tau = RieszRep::repFunction(tau,rieszResidual);
  // experiment by Nate: manually specify the error (this appears to produce identical results, as it should)
//  FunctionPtr err = e_v * e_v + e_tau * e_tau + e_v->grad() * e_v->grad() + e_tau->div() * e_tau->div();
  map<int,FunctionPtr> errFxns;
  errFxns[v->ID()] = e_v;
  errFxns[tau->ID()] = e_tau;
  LinearTermPtr ipAtErrFxns = ip->evaluate(errFxns);
  FunctionPtr err = ip->evaluate(errFxns)->evaluate(errFxns);
  double energyErrorIntegrated = sqrt(err->integrate(mesh,cubEnrich,testVsTest));

  // check that energy error computed thru Solution and through rieszRep are the same
  bool success1 = abs(energyError-energyErrorLT)<tol;
  // checks that matrix-computed and integrated errors are the same
  bool success2 = abs(energyErrorLT-energyErrorIntegrated)<tol;
  success = success1==true && success2==true;
  if (!success)
  {
    if (rank==0)
      cout << "Failed testLTResidual; energy error = " << energyError << ", while linearTerm error is computed to be " << energyErrorLT << ", and when computing through integration of the Riesz rep function, error = " << energyErrorIntegrated << endl;
  }
  //  VTKExporter exporter(solution, mesh, varFactory);
  //  exporter.exportSolution("testLTRes");
  //  cout << endl;

  return success;
}

// tests residual computation on simple convection
bool ScratchPadTests::testLTResidualSimple()
{
  double tol = 1e-11;
  int rank = Teuchos::GlobalMPISession::getRank();

  bool success = true;

  int nCells = 2;

  ////////////////////   DECLARE VARIABLES   ///////////////////////

  // define test variables
  VarFactoryPtr varFactory = VarFactory::varFactory();
  VarPtr v = varFactory->testVar("v", HGRAD);

  // define trial variables
  VarPtr beta_n_u = varFactory->fluxVar("\\widehat{\\beta \\cdot n u - \\sigma_{n}}");
  VarPtr u = varFactory->fieldVar("u");

  vector<double> beta;
  beta.push_back(1.0);
  beta.push_back(1.0);

  ////////////////////   DEFINE BILINEAR FORM   ///////////////////////

  BFPtr confusionBF = Teuchos::rcp( new BF(varFactory) );
  // v terms:
  confusionBF->addTerm( -u, beta * v->grad() );
  confusionBF->addTerm( beta_n_u, v);

  ////////////////////   DEFINE INNER PRODUCT(S)   ///////////////////////

  // robust test norm
  IPPtr ip = Teuchos::rcp(new IP);

  // choose the mesh-independent norm even though it may have BLs
  ip->addTerm(v->grad());
  ip->addTerm(v);

  ////////////////////   SPECIFY RHS AND HELPFUL FUNCTIONS   ///////////////////////

  FunctionPtr n = Function::normal();
  vector<double> e1,e2;
  e1.push_back(1.0);
  e1.push_back(0.0);
  e2.push_back(0.0);
  e2.push_back(1.0);
  FunctionPtr one = Function::constant(1.0);

  FunctionPtr zero = Function::constant(0.0);
  RHSPtr rhs = RHS::rhs();
  FunctionPtr f = one;
  rhs->addTerm( f * v );

  ////////////////////   CREATE BCs   ///////////////////////
  BCPtr bc = BC::bc();
  SpatialFilterPtr boundary = Teuchos::rcp( new InflowSquareBoundary );
  FunctionPtr u_in = Teuchos::rcp(new Uinflow);
  bc->addDirichlet(beta_n_u, boundary, beta*n*u_in);

  ////////////////////   BUILD MESH   ///////////////////////
  // define nodes for mesh
  int order = 2;
  int H1Order = order+1;
  int pToAdd = 2;

  // create a pointer to a new mesh:
  Teuchos::RCP<Mesh> mesh = MeshUtilities::buildUnitQuadMesh(nCells,confusionBF, H1Order, H1Order+pToAdd);

  ////////////////////   SOLVE & REFINE   ///////////////////////

  int cubEnrich = 0;

  Teuchos::RCP<Solution> solution;
  solution = Teuchos::rcp( new Solution(mesh, bc, rhs, ip) );
  solution->solve(false);
  double energyError = solution->energyErrorTotal();

  LinearTermPtr residual = rhs->linearTermCopy();
  residual->addTerm(-confusionBF->testFunctional(solution),true);

  Teuchos::RCP<RieszRep> rieszResidual = Teuchos::rcp(new RieszRep(mesh, ip, residual));
  rieszResidual->computeRieszRep(cubEnrich);
  double energyErrorLT = rieszResidual->getNorm();

  bool testVsTest = true;
  FunctionPtr e_v = RieszRep::repFunction(v,rieszResidual);
  map<int,FunctionPtr> errFxns;
  errFxns[v->ID()] = e_v;
  FunctionPtr err = (ip->evaluate(errFxns,false))->evaluate(errFxns,false); // don't need boundary terms unless they're in IP
  double energyErrorIntegrated = sqrt(err->integrate(mesh,cubEnrich,testVsTest));
  // check that energy error computed thru Solution and through rieszRep are the same
  success = abs(energyError-energyErrorLT) < tol;
  if (success==false)
  {
    if (rank==0)
      cout << "Failed testLTResidualSimple; energy error = " << energyError << ", while linearTerm error is computed to be " << energyErrorLT << endl;
    return success;
  }
//.........这里部分代码省略.........