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

int main(int argc, char *argv[]) {

#ifdef HAVE_MPI
  Teuchos::GlobalMPISession mpiSession(&argc, &argv,0);
  choice::MpiArgs args( argc, argv );
#else
  choice::Args args( argc, argv );
#endif
  int rank = Teuchos::GlobalMPISession::getRank();
  int numProcs = Teuchos::GlobalMPISession::getNProc();

  int nCells = args.Input<int>("--nCells", "num cells",2);
  int numSteps = args.Input<int>("--numSteps", "num NR steps",20);

  int polyOrder = 0;
  
  // define our manufactured solution or problem bilinear form:
  bool useTriangles = false;
  
  int pToAdd = 1;

  args.Process();

  int H1Order = polyOrder + 1;
  
  ////////////////////////////////////////////////////////////////////
  // DEFINE VARIABLES 
  ////////////////////////////////////////////////////////////////////
  
  // new-style bilinear form definition
  VarFactory varFactory;
  VarPtr fn = varFactory.fluxVar("\\widehat{\\beta_n_u}");
  VarPtr u = varFactory.fieldVar("u");
  
  VarPtr v = varFactory.testVar("v",HGRAD);
  BFPtr bf = Teuchos::rcp( new BF(varFactory) ); // initialize bilinear form
  
  ////////////////////////////////////////////////////////////////////
  // CREATE MESH 
  ////////////////////////////////////////////////////////////////////
  
  // create a pointer to a new mesh:
  Teuchos::RCP<Mesh> mesh = MeshUtilities::buildUnitQuadMesh(nCells , bf, H1Order, H1Order+pToAdd);
  
  ////////////////////////////////////////////////////////////////////
  // INITIALIZE BACKGROUND FLOW FUNCTIONS
  ////////////////////////////////////////////////////////////////////
  BCPtr nullBC = Teuchos::rcp((BC*)NULL); RHSPtr nullRHS = Teuchos::rcp((RHS*)NULL); IPPtr nullIP = Teuchos::rcp((IP*)NULL);
  SolutionPtr backgroundFlow = Teuchos::rcp(new Solution(mesh, nullBC, nullRHS, nullIP) );
  SolutionPtr solnPerturbation = Teuchos::rcp(new Solution(mesh, nullBC, nullRHS, nullIP) );
  
  vector<double> e1(2),e2(2);
  e1[0] = 1; e2[1] = 1;
  
  FunctionPtr u_prev = Teuchos::rcp( new PreviousSolutionFunction(backgroundFlow, u) );
  FunctionPtr beta = e1 * u_prev + Teuchos::rcp( new ConstantVectorFunction( e2 ) );
  
  ////////////////////////////////////////////////////////////////////
  // DEFINE BILINEAR FORM
  ////////////////////////////////////////////////////////////////////
  
  // v:
  bf->addTerm( -u, beta * v->grad());
  bf->addTerm( fn, v);

  ////////////////////////////////////////////////////////////////////
  // DEFINE RHS
  ////////////////////////////////////////////////////////////////////

  Teuchos::RCP<RHSEasy> rhs = Teuchos::rcp( new RHSEasy );
  FunctionPtr u_prev_squared_div2 = 0.5 * u_prev * u_prev;  
  rhs->addTerm((e1 * u_prev_squared_div2 + e2 * u_prev) * v->grad());

  // ==================== SET INITIAL GUESS ==========================

  mesh->registerSolution(backgroundFlow);
  FunctionPtr zero = Function::constant(0.0);
  FunctionPtr u0 = Teuchos::rcp( new U0 );
  FunctionPtr n = Teuchos::rcp( new UnitNormalFunction );
  //  FunctionPtr parity = Teuchos::rcp(new SideParityFunction);

  FunctionPtr u0_squared_div_2 = 0.5 * u0 * u0;

  map<int, Teuchos::RCP<Function> > functionMap;
  functionMap[u->ID()] = u0;
  //  functionMap[fn->ID()] = -(e1 * u0_squared_div_2 + e2 * u0) * n * parity;
  backgroundFlow->projectOntoMesh(functionMap);

  // ==================== END SET INITIAL GUESS ==========================

  ////////////////////////////////////////////////////////////////////
  // DEFINE INNER PRODUCT
  ////////////////////////////////////////////////////////////////////

  IPPtr ip = Teuchos::rcp( new IP );
  ip->addTerm( v );
  ip->addTerm(v->grad());
  //  ip->addTerm( beta * v->grad() ); // omitting term to make IP non-dependent on u

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

bool LinearTermTests::testIntegration()
{
  // for now, we just check the consistency: for LinearTerm a = b + c, does a->integrate
  // give the same values as b->integrate + c->integrate ?
  bool success = true;

  //  VarPtr v1, v2, v3; // HGRAD members (test variables)
  //  VarPtr q1, q2, q3; // HDIV members (test variables)
  //  VarPtr u1, u2, u3; // L2 members (trial variables)
  //  VarPtr u1_hat, u2_hat; // trace variables
  //  VarPtr u3_hat_n; // flux variable
  //
  //  FunctionPtr sine_x;

  if ( ! checkLTSumConsistency(1 * v1, 1 * v2, testOrder, basisCache) )
  {
    cout << "(v1 + v2)->integrate not consistent with sum of summands integration.\n";
    success = false;
  }

  if ( ! checkLTSumConsistency(sine_x * v1, 1 * v2, testOrder, basisCache) )
  {
    cout << "(sine_x * v1 + v2)->integrate not consistent with sum of summands integration.\n";
    success = false;
  }

  if ( ! checkLTSumConsistency(1 * q1->div(), 1 * q2->x(), testOrder, basisCache) )
  {
    cout << "(q1->div() + q2->x())->integrate not consistent with sum of summands integration.\n";
    success = false;
  }

  if ( ! checkLTSumConsistency(1 * u1, 1 * u2, trialOrder, basisCache) )
  {
    cout << "(u1 + u2)->integrate not consistent with sum of summands integration.\n";
    success = false;
  }

  if ( ! checkLTSumConsistency(1 * u1, sine_x * u2, trialOrder, basisCache) )
  {
    cout << "(u1 + sine_x * u2)->integrate not consistent with sum of summands integration.\n";
    success = false;
  }

  // now, same thing, but with boundary-value-only functions in the mix:
  // this next is a fairly complex test; may want to add a more granular one above...
  IPPtr ip = Teuchos::rcp(new IP);
  RHSPtr rhs = RHS::rhs();
  BCPtr bc = BC::bc();
  SolutionPtr solution = Teuchos::rcp( new Solution(mesh,bc,rhs,ip) );
  // project some functions onto solution, so that something interesting is there:
  FunctionPtr u1_proj = sine_x;
  FunctionPtr u2_proj = cos_y;
  FunctionPtr u3_proj = u1_proj * u2_proj;
  map<int, FunctionPtr> solnToProject;
  solnToProject[u1->ID()] = u1_proj;
  solnToProject[u2->ID()] = u2_proj;
  solnToProject[u3->ID()] = u3_proj;
  solnToProject[u1_hat->ID()] = u1_proj;
  solnToProject[u2_hat->ID()] = u2_proj;
  // u3_hat_n isn't too much like a 'real' bilinear form, in that u3 itself is a scalar
  // this is just a test, so I'm not worried about it...
  solnToProject[u3_hat_n->ID()] = u3_proj;

  solution->projectOntoMesh(solnToProject);

  LinearTermPtr bfTestFunctional = bf->testFunctional(solution);

  // bf->addTerm(u1, q1->x());
  // bf->addTerm(u2, q1->y());
  // bf->addTerm(u3, v1);

//  bf->addTerm(u1_hat, q1->dot_normal());
//  bf->addTerm(u3_hat_n, v1);

  LinearTermPtr testFunctionalNoBoundaryValues = u1_proj * q1->x() + u2_proj * q1->y() + u3_proj * v1;

  FunctionPtr u1_hat_prev = Teuchos::rcp( new PreviousSolutionFunction<double>(solution, u1_hat) );
  FunctionPtr u2_hat_prev = Teuchos::rcp( new PreviousSolutionFunction<double>(solution, u2_hat) );
  FunctionPtr u3_hat_prev = Teuchos::rcp( new PreviousSolutionFunction<double>(solution, u3_hat_n) );
  LinearTermPtr testFunctionalBoundaryValues = u1_hat_prev * q1->dot_normal() + u3_hat_prev * v1;

  if ( ! checkLTSumConsistency(testFunctionalNoBoundaryValues, testFunctionalBoundaryValues,
                               testOrder, basisCache) )
  {
    cout << "bfTestFunctional->integrate not consistent with sum of summands integration.\n";
    success = false;
  }

  if ( ! checkLTSumConsistency(testFunctionalBoundaryValues, bfTestFunctional - testFunctionalBoundaryValues,
                               testOrder, basisCache) )
  {
    cout << "bfTestFunctional->integrate not consistent with sum of summands integration.\n";
    success = false;
  }

  if ( ! checkLTSumConsistency(testFunctionalNoBoundaryValues, bfTestFunctional - testFunctionalNoBoundaryValues,
                               testOrder, basisCache) )
  {
    cout << "bfTestFunctional->integrate not consistent with sum of summands integration.\n";
//.........这里部分代码省略.........

int main(int argc, char *argv[])
{
  // Process command line arguments
#ifdef HAVE_MPI
  Teuchos::GlobalMPISession mpiSession(&argc, &argv,0);
  int rank=mpiSession.getRank();
  int numProcs=mpiSession.getNProc();
#else
  int rank = 0;
  int numProcs = 1;
#endif

  ////////////////////   DECLARE VARIABLES   ///////////////////////
  // define test variables
  VarFactory varFactory;
  VarPtr v = varFactory.testVar("v", HGRAD);

  // define trial variables
  VarPtr beta_n_u_hat = varFactory.fluxVar("\\widehat{\\beta \\cdot n }");
  VarPtr u = varFactory.fieldVar("u");

  FunctionPtr beta = Teuchos::rcp(new Beta());

  ////////////////////   BUILD MESH   ///////////////////////
  BFPtr confusionBF = Teuchos::rcp( new BF(varFactory) );
  // define nodes for mesh
  FieldContainer<double> meshBoundary(4,2);

  meshBoundary(0,0) = -1.0; // x1
  meshBoundary(0,1) = -1.0; // y1
  meshBoundary(1,0) =  1.0;
  meshBoundary(1,1) = -1.0;
  meshBoundary(2,0) =  1.0;
  meshBoundary(2,1) =  1.0;
  meshBoundary(3,0) = -1.0;
  meshBoundary(3,1) =  1.0;

  int horizontalCells = 32, verticalCells = 32;

  // create a pointer to a new mesh:
  Teuchos::RCP<Mesh> mesh = Mesh::buildQuadMesh(meshBoundary, horizontalCells, verticalCells,
                            confusionBF, H1Order, H1Order+pToAdd);

  ////////////////////////////////////////////////////////////////////
  // INITIALIZE FLOW FUNCTIONS
  ////////////////////////////////////////////////////////////////////

  BCPtr nullBC = Teuchos::rcp((BC*)NULL);
  RHSPtr nullRHS = Teuchos::rcp((RHS*)NULL);
  IPPtr nullIP = Teuchos::rcp((IP*)NULL);
  SolutionPtr prevTimeFlow = Teuchos::rcp(new Solution(mesh, nullBC, nullRHS, nullIP) );
  SolutionPtr flowResidual = Teuchos::rcp(new Solution(mesh, nullBC, nullRHS, nullIP) );

  FunctionPtr u_prev_time = Teuchos::rcp( new PreviousSolutionFunction(prevTimeFlow, u) );

  ////////////////////   DEFINE BILINEAR FORM   ///////////////////////
  Teuchos::RCP<RHSEasy> rhs = Teuchos::rcp( new RHSEasy );
  FunctionPtr invDt = Teuchos::rcp(new ScalarParamFunction(1.0/dt));

  // v terms:
  confusionBF->addTerm( beta * u, - v->grad() );
  confusionBF->addTerm( beta_n_u_hat, v);

  confusionBF->addTerm( u, invDt*v );
  rhs->addTerm( u_prev_time * invDt * v );

  ////////////////////   SPECIFY RHS   ///////////////////////
  FunctionPtr f = Teuchos::rcp( new ConstantScalarFunction(0.0) );
  rhs->addTerm( f * v ); // obviously, with f = 0 adding this term is not necessary!

  ////////////////////   DEFINE INNER PRODUCT(S)   ///////////////////////
  // robust test norm
  IPPtr ip = confusionBF->graphNorm();
  // IPPtr ip = Teuchos::rcp(new IP);
  // ip->addTerm(v);
  // ip->addTerm(invDt*v - beta*v->grad());

  ////////////////////   CREATE BCs   ///////////////////////
  Teuchos::RCP<BCEasy> bc = Teuchos::rcp( new BCEasy );
  SpatialFilterPtr inflowBoundary = Teuchos::rcp( new InflowSquareBoundary(beta) );
  FunctionPtr u0 = Teuchos::rcp( new ConstantScalarFunction(0) );
  FunctionPtr n = Teuchos::rcp( new UnitNormalFunction );

  bc->addDirichlet(beta_n_u_hat, inflowBoundary, beta*n*u0);

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

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

  // ==================== SET INITIAL GUESS ==========================
  FunctionPtr u_init = Teuchos::rcp(new InitialCondition());
  map<int, Teuchos::RCP<Function> > functionMap;
  functionMap[u->ID()]      = u_init;

  prevTimeFlow->projectOntoMesh(functionMap);

  ////////////////////   SOLVE & REFINE   ///////////////////////
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  BCPtr nullBC = Teuchos::rcp((BC*)NULL);
  RHSPtr nullRHS = Teuchos::rcp((RHS*)NULL);
  IPPtr nullIP = Teuchos::rcp((IP*)NULL);
  SolutionPtr backgroundFlow = Teuchos::rcp(new Solution(mesh, nullBC, nullRHS, nullIP) );

  vector<double> e1(2); // (1,0)
  e1[0] = 1;
  vector<double> e2(2); // (0,1)
  e2[1] = 1;

  FunctionPtr u1_prev = Function::solution(u1, backgroundFlow);
  FunctionPtr u2_prev = Function::solution(u2, backgroundFlow);
  FunctionPtr sigma1_prev = Function::solution(sigma1, backgroundFlow);
  FunctionPtr sigma2_prev = Function::solution(sigma2, backgroundFlow);
  FunctionPtr p_prev = Function::solution(p, backgroundFlow);
  // FunctionPtr sigma11_prev = Function::solution(sigma11, backgroundFlow);
  // FunctionPtr sigma12_prev = Function::solution(sigma12, backgroundFlow);
  // FunctionPtr sigma22_prev = Function::solution(sigma22, backgroundFlow);

  FunctionPtr zero = Teuchos::rcp( new ConstantScalarFunction(0.0) );
  FunctionPtr one = Teuchos::rcp( new ConstantScalarFunction(1.0) );
  FunctionPtr u1Exact     = Teuchos::rcp( new ExactU1(lambda) );
  FunctionPtr u2Exact     = Teuchos::rcp( new ExactU2(lambda) );
  // FunctionPtr beta = e1 * u1_prev + e2 * u2_prev;

  // ==================== SET INITIAL GUESS ==========================
  map<int, Teuchos::RCP<Function> > functionMap;
  functionMap[u1->ID()] = u1Exact;
  functionMap[u2->ID()] = u2Exact;
  // functionMap[sigma1->ID()] = Function::vectorize(zero,zero);
  // functionMap[sigma2->ID()] = Function::vectorize(zero,zero);
  // functionMap[p->ID()] = zero;

  backgroundFlow->projectOntoMesh(functionMap);

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

  // // stress equation
  bf->addTerm( 1./nu*sigma1, tau1 );
  bf->addTerm( 1./nu*sigma2, tau2 );
  bf->addTerm( u1, tau1->div() );
  bf->addTerm( u2, tau2->div() );
  bf->addTerm( -u1hat, tau1->dot_normal() );
  bf->addTerm( -u2hat, tau2->dot_normal() );
  // bf->addTerm( 1./(2*nu)*sigma11, tau11 );
  // bf->addTerm( 1./(2*nu)*sigma12, tau12 );
  // bf->addTerm( 1./(2*nu)*sigma12, tau12 );
  // bf->addTerm( 1./(2*nu)*sigma22, tau22 );
  // bf->addTerm( u1, tau11->dx() );
  // bf->addTerm( u1, tau12->dy() );
  // bf->addTerm( u2, tau12->dx() );
  // bf->addTerm( u2, tau22->dy() );
  // bf->addTerm( -u1hat, tau11->times_normal_x() );
  // bf->addTerm( -u1hat, tau12->times_normal_y() );
  // bf->addTerm( -u2hat, tau12->times_normal_x() );
  // bf->addTerm( -u2hat, tau22->times_normal_y() );

  // momentum equation
  bf->addTerm( -2.*u1_prev*u1, v1->dx() );
  bf->addTerm( -u2_prev*u1, v1->dy() );
  bf->addTerm( -u1_prev*u2, v1->dy() );
  bf->addTerm( -u2_prev*u1, v2->dx() );
  bf->addTerm( -u1_prev*u2, v1->dy() );
  bf->addTerm( -2.*u2_prev*u2, v2->dy() );
  bf->addTerm( -p, v1->dx() );
  bf->addTerm( -p, v2->dy() );

//.........这里部分代码省略.........
  LinearTermPtr restLT = Teuchos::rcp(new LinearTerm);
  restLT->addTerm(alpha*v);
  restLT->addTerm(invSqrtH*v);
  restLT = restLT + beta * v->grad();
  restLT = restLT + tau->div();

  ////////////////////   SPECIFY RHS   ///////////////////////

  Teuchos::RCP<RHSEasy> rhs = Teuchos::rcp( new RHSEasy );
  FunctionPtr f = zero;
  //  f = one;
  rhs->addTerm( f * v ); // obviously, with f = 0 adding this term is not necessary!

  ////////////////////   CREATE BCs   ///////////////////////
  Teuchos::RCP<BCEasy> bc = Teuchos::rcp( new BCEasy );

  SpatialFilterPtr Inflow = Teuchos::rcp(new LeftInflow);
  SpatialFilterPtr wallBoundary = Teuchos::rcp(new WallBoundary);//MeshUtilities::rampBoundary(rampHeight);
  SpatialFilterPtr freeStream = Teuchos::rcp(new FreeStreamBoundary);

  bc->addDirichlet(uhat, wallBoundary, one);
  //  bc->addDirichlet(uhat, wallBoundary, Teuchos::rcp(new WallSmoothBC(eps)));
  bc->addDirichlet(beta_n_u_minus_sigma_n, Inflow, zero);
  bc->addDirichlet(beta_n_u_minus_sigma_n, freeStream, zero);

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

  Teuchos::RCP<Solution> solution;
  solution = Teuchos::rcp( new Solution(mesh, bc, rhs, ip) );
  BCPtr nullBC = Teuchos::rcp((BC*)NULL); RHSPtr nullRHS = Teuchos::rcp((RHS*)NULL); IPPtr nullIP = Teuchos::rcp((IP*)NULL);
  SolutionPtr backgroundFlow = Teuchos::rcp(new Solution(mesh, nullBC, nullRHS, nullIP) );  
  mesh->registerSolution(backgroundFlow); // to trigger issue with p-refinements
  map<int, Teuchos::RCP<Function> > functionMap; functionMap[u->ID()] = Function::constant(3.14);
  backgroundFlow->projectOntoMesh(functionMap);

  // lower p to p = 1 at SINGULARITY only
  vector<int> ids;
  /*
  for (int i = 0;i<mesh->numActiveElements();i++){
    bool cellIDset = false;
    int cellID = mesh->activeElements()[i]->cellID();
    int elemOrder = mesh->cellPolyOrder(cellID)-1;
    FieldContainer<double> vv(4,2); mesh->verticesForCell(vv, cellID);
    bool vertexOnWall = false; bool vertexAtSingularity = false;
    for (int j = 0;j<4;j++){
      if ((abs(vv(j,0)-.5) + abs(vv(j,1)))<1e-10){
	vertexAtSingularity = true;     
	cellIDset = true;
      }
    }	
    if (!vertexAtSingularity && elemOrder<2 && !cellIDset ){
      ids.push_back(cellID);
      cout << "celliD = " << cellID << endl;
    }
  }
  */
  ids.push_back(1);
  ids.push_back(3);
  mesh->pRefine(ids); // to put order = 1

  return 0;
  
  LinearTermPtr residual = rhs->linearTermCopy();
  residual->addTerm(-confusionBF->testFunctional(solution));  
  RieszRepPtr rieszResidual = Teuchos::rcp(new RieszRep(mesh, ip, residual));
  rieszResidual->computeRieszRep();

//.........这里部分代码省略.........
  
  qoptIP = Teuchos::rcp(new IP());
      
  if (useCompliantGraphNorm) {
    qoptIP->addTerm( mu * v1->dx() + tau1->x() ); // sigma11
    qoptIP->addTerm( mu * v1->dy() + tau1->y() ); // sigma12
    qoptIP->addTerm( mu * v2->dx() + tau2->x() ); // sigma21
    qoptIP->addTerm( mu * v2->dy() + tau2->y() ); // sigma22
    qoptIP->addTerm( mu * v1->dx() + mu * v2->dy() );   // pressure
    qoptIP->addTerm( h * tau1->div() - h * q->dx() );   // u1
    qoptIP->addTerm( h * tau2->div() - h * q->dy());    // u2
    
    qoptIP->addTerm( (mu / h) * v1 );
    qoptIP->addTerm( (mu / h) * v2 );
    qoptIP->addTerm( q );
    qoptIP->addTerm( tau1 );
    qoptIP->addTerm( tau2 );
  } else { // standard graph norm, then
    qoptIP = stokesBF->graphNorm();
  }

  ip = qoptIP;
  
  if (rank==0) 
    ip->printInteractions();
  
  // aim is just to answer one simple question:
  // have I figured out a trial-space preimage for optimal test function (q=1, tau=0, v=0)?
  
  SolutionPtr soln = Teuchos::rcp(new Solution(mesh));
  
  FunctionPtr x = Function::xn();
  FunctionPtr y = Function::yn();
  
  // u1 = u1_hat = x / 2
  FunctionPtr u1_exact = x / 2;
  
  // u2 = u2_hat = y / 2
  FunctionPtr u2_exact = y / 2;
  
  // sigma = 0.5 * I
  FunctionPtr sigma11_exact = Function::constant(0.5);
  FunctionPtr sigma22_exact = Function::constant(0.5);
  
  // tn_hat = 0.5 * n
  FunctionPtr n = Function::normal();
  FunctionPtr t1n_exact = n->x() / 2;
  FunctionPtr t2n_exact = n->y() / 2;
  
  map<int, FunctionPtr > exact_soln;
  exact_soln[u1->ID()] = u1_exact;
  exact_soln[u1hat->ID()] = u1_exact;
  exact_soln[u2->ID()] = u2_exact;
  exact_soln[u2hat->ID()] = u2_exact;
  exact_soln[sigma11->ID()] = sigma11_exact;
  exact_soln[sigma22->ID()] = sigma22_exact;
  exact_soln[t1n->ID()] = t1n_exact;
  exact_soln[t2n->ID()] = t2n_exact;
  
  exact_soln[p->ID()] = Function::zero();
  exact_soln[sigma12->ID()] = Function::zero();
  exact_soln[sigma21->ID()] = Function::zero();
  
  soln->projectOntoMesh(exact_soln);
  
  LinearTermPtr soln_functional = stokesBF->testFunctional(soln);
  
  RieszRepPtr rieszRep = Teuchos::rcp( new RieszRep(mesh, ip, soln_functional) );
  
  rieszRep->computeRieszRep();
  
  // get test functions:
  FunctionPtr q_fxn = Teuchos::rcp( new RepFunction(q, rieszRep) );
  FunctionPtr v1_fxn = Teuchos::rcp( new RepFunction(v1, rieszRep) );
  FunctionPtr v2_fxn = Teuchos::rcp( new RepFunction(v2, rieszRep) );
  FunctionPtr tau1_fxn = Teuchos::rcp( new RepFunction(tau1, rieszRep) );
  FunctionPtr tau2_fxn = Teuchos::rcp( new RepFunction(tau2, rieszRep) );
  
  cout << "L2 norm of (q-1) : " << (q_fxn - 1)->l2norm(mesh) << endl;
  cout << "L2 norm of (v1) : " << (v1_fxn)->l2norm(mesh) << endl;
  cout << "L2 norm of (v2) : " << (v2_fxn)->l2norm(mesh) << endl;
  cout << "L2 norm of (tau1) : " << (tau1_fxn)->l2norm(mesh) << endl;
  cout << "L2 norm of (tau2) : " << (tau2_fxn)->l2norm(mesh) << endl;
  
  VTKExporter exporter(soln, mesh, varFactory);
  exporter.exportSolution("conservationPreimage", H1Order*2);

  cout << "Checking that the soln_functional is what I expect:\n";
  
  FunctionPtr xyVector = Function::vectorize(x, y);
  
  cout << "With v1 = x, integral: " << integralOverMesh(soln_functional, v1, x) << endl;
  cout << "With v2 = y, integral: " << integralOverMesh(soln_functional, v2, y) << endl;
  cout << "With tau1=(x,y), integral: " << integralOverMesh(soln_functional, tau1, xyVector) << endl;
  cout << "With tau2=(x,y), integral: " << integralOverMesh(soln_functional, tau2, xyVector) << endl;
  cout << "With q   =x, integral: " << integralOverMesh(soln_functional, q, x) << endl;
  
  cout << "(Expect 0s all around, except for q, where we expect (1,x) == 0.5.)\n";
  return 0;
}

int main(int argc, char *argv[]) {
  // Process command line arguments
  if (argc > 1)
    numRefs = atof(argv[1]);
#ifdef HAVE_MPI
  Teuchos::GlobalMPISession mpiSession(&argc, &argv,0);
  int rank=mpiSession.getRank();
  int numProcs=mpiSession.getNProc();
#else
  int rank = 0;
  int numProcs = 1;
#endif

  FunctionPtr beta = Teuchos::rcp(new Beta());

  ////////////////////////////////////////////////////////////////////
  // DEFINE VARIABLES 
  ////////////////////////////////////////////////////////////////////
  // test variables
  VarFactory varFactory; 
  VarPtr tau = varFactory.testVar("\\tau", HDIV);
  VarPtr v = varFactory.testVar("v", HGRAD);

  // trial variables
  VarPtr uhat = varFactory.traceVar("\\widehat{u}");
  VarPtr beta_n_u_minus_sigma_n = varFactory.fluxVar("\\widehat{\\beta \\cdot n u - \\sigma_{n}}");
  VarPtr u = varFactory.fieldVar("u");
  VarPtr sigma1 = varFactory.fieldVar("\\sigma_1");
  VarPtr sigma2 = varFactory.fieldVar("\\sigma_2");

  ////////////////////////////////////////////////////////////////////
  // CREATE MESH 
  ////////////////////////////////////////////////////////////////////

  BFPtr confusionBF = Teuchos::rcp( new BF(varFactory) );

  FieldContainer<double> meshBoundary(4,2);

  meshBoundary(0,0) =  0.0; // x1
  meshBoundary(0,1) = -2.0; // y1
  meshBoundary(1,0) =  4.0;
  meshBoundary(1,1) = -2.0;
  meshBoundary(2,0) =  4.0;
  meshBoundary(2,1) =  2.0;
  meshBoundary(3,0) =  0.0;
  meshBoundary(3,1) =  2.0;

  int horizontalCells = 4, verticalCells = 4;

  // create a pointer to a new mesh:
  Teuchos::RCP<Mesh> mesh = Mesh::buildQuadMesh(meshBoundary, horizontalCells, verticalCells,
      confusionBF, H1Order, H1Order+pToAdd, false);

  ////////////////////////////////////////////////////////////////////
  // INITIALIZE BACKGROUND FLOW FUNCTIONS
  ////////////////////////////////////////////////////////////////////

  BCPtr nullBC = Teuchos::rcp((BC*)NULL);
  RHSPtr nullRHS = Teuchos::rcp((RHS*)NULL);
  IPPtr nullIP = Teuchos::rcp((IP*)NULL);
  SolutionPtr prevTimeFlow = Teuchos::rcp(new Solution(mesh, nullBC, nullRHS, nullIP) );  
  SolutionPtr flowResidual = Teuchos::rcp(new Solution(mesh, nullBC, nullRHS, nullIP) );  

  FunctionPtr u_prev_time = Teuchos::rcp( new PreviousSolutionFunction(prevTimeFlow, u) );

  // ==================== SET INITIAL GUESS ==========================
  double u_free = 0.0;
  double sigma1_free = 0.0;
  double sigma2_free = 0.0;
  map<int, Teuchos::RCP<Function> > functionMap;
  functionMap[u->ID()] = Teuchos::rcp( new ConstantScalarFunction(u_free) );
  functionMap[sigma1->ID()] = Teuchos::rcp( new ConstantScalarFunction(sigma1_free) );
  functionMap[sigma2->ID()] = Teuchos::rcp( new ConstantScalarFunction(sigma2_free) );

  prevTimeFlow->projectOntoMesh(functionMap);
  // ==================== END SET INITIAL GUESS ==========================

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

  // tau terms:
  confusionBF->addTerm(sigma1 / epsilon, tau->x());
  confusionBF->addTerm(sigma2 / epsilon, tau->y());
  confusionBF->addTerm(u, tau->div());
  confusionBF->addTerm(-uhat, tau->dot_normal());

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

  ////////////////////////////////////////////////////////////////////
  // TIMESTEPPING TERMS
  ////////////////////////////////////////////////////////////////////
  Teuchos::RCP<RHSEasy> rhs = Teuchos::rcp( new RHSEasy );

  double dt = 0.25;
  FunctionPtr invDt = Teuchos::rcp(new ScalarParamFunction(1.0/dt));    
//.........这里部分代码省略.........

int main(int argc, char *argv[])
{
#ifdef HAVE_MPI
  Teuchos::GlobalMPISession mpiSession(&argc, &argv,0);
  choice::MpiArgs args( argc, argv );
#else
  choice::Args args( argc, argv );
#endif
  int commRank = Teuchos::GlobalMPISession::getRank();
  int numProcs = Teuchos::GlobalMPISession::getNProc();

  // Required arguments
  int numRefs = args.Input<int>("--numRefs", "number of refinement steps");
  int norm = args.Input<int>("--norm", "0 = graph\n    1 = robust\n    2 = coupled robust");

  // Optional arguments (have defaults)
  int uniformRefinements = args.Input("--uniformRefinements", "number of uniform refinements", 0);
  bool enforceLocalConservation = args.Input<bool>("--conserve", "enforce local conservation", false);
  double radius = args.Input("--r", "cylinder radius", 0.6);
  int Re = args.Input("--Re", "Reynolds number", 1);
  int maxNewtonIterations = args.Input("--maxIterations", "maximum number of Newton iterations", 1);
  int polyOrder = args.Input("--polyOrder", "polynomial order for field variables", 2);
  int deltaP = args.Input("--deltaP", "how much to enrich test space", 2);
  // string saveFile = args.Input<string>("--meshSaveFile", "file to which to save refinement history", "");
  // string replayFile = args.Input<string>("--meshLoadFile", "file with refinement history to replay", "");
  args.Process();

  ////////////////////   PROBLEM DEFINITIONS   ///////////////////////
  int H1Order = polyOrder+1;

  ////////////////////   DECLARE VARIABLES   ///////////////////////
  // define test variables
  VarFactory varFactory;
  VarPtr tau1 = varFactory.testVar("tau1", HDIV);
  VarPtr tau2 = varFactory.testVar("tau2", HDIV);
  VarPtr v1 = varFactory.testVar("v1", HGRAD);
  VarPtr v2 = varFactory.testVar("v2", HGRAD);
  VarPtr vc = varFactory.testVar("vc", HGRAD);

  // define trial variables
  VarPtr u1 = varFactory.fieldVar("u1");
  VarPtr u2 = varFactory.fieldVar("u2");
  VarPtr p = varFactory.fieldVar("p");
  VarPtr u1hat = varFactory.traceVar("u1hat");
  VarPtr u2hat = varFactory.traceVar("u2hat");
  VarPtr t1hat = varFactory.fluxVar("t1hat");
  VarPtr t2hat = varFactory.fluxVar("t2hat");
  VarPtr sigma1 = varFactory.fieldVar("sigma1", VECTOR_L2);
  VarPtr sigma2 = varFactory.fieldVar("sigma2", VECTOR_L2);

  ////////////////////   BUILD MESH   ///////////////////////
  BFPtr bf = Teuchos::rcp( new BF(varFactory) );

  // create a pointer to a new mesh:
  Teuchos::RCP<Mesh> mesh = MeshFactory::shiftedHemkerMesh(-1, 3, 2, radius, bf, H1Order, deltaP);

  ////////////////////////////////////////////////////////////////////
  // INITIALIZE BACKGROUND FLOW FUNCTIONS
  ////////////////////////////////////////////////////////////////////

  BCPtr nullBC = Teuchos::rcp((BC*)NULL);
  RHSPtr nullRHS = Teuchos::rcp((RHS*)NULL);
  IPPtr nullIP = Teuchos::rcp((IP*)NULL);
  SolutionPtr backgroundFlow = Teuchos::rcp(new Solution(mesh, nullBC, nullRHS, nullIP) );

  vector<double> e1(2); // (1,0)
  e1[0] = 1;
  vector<double> e2(2); // (0,1)
  e2[1] = 1;

  FunctionPtr u1_prev = Function::solution(u1, backgroundFlow);
  FunctionPtr u2_prev = Function::solution(u2, backgroundFlow);
  FunctionPtr sigma1_prev = Function::solution(sigma1, backgroundFlow);
  FunctionPtr sigma2_prev = Function::solution(sigma2, backgroundFlow);

  FunctionPtr zero = Teuchos::rcp( new ConstantScalarFunction(0.0) );
  FunctionPtr one = Teuchos::rcp( new ConstantScalarFunction(1.0) );
  FunctionPtr beta = e1 * u1_prev + e2 * u2_prev;

  // ==================== SET INITIAL GUESS ==========================
  map<int, Teuchos::RCP<Function> > functionMap;
  functionMap[u1->ID()] = one;
  functionMap[u2->ID()] = zero;
  functionMap[sigma1->ID()] = Function::vectorize(zero,zero);
  functionMap[sigma2->ID()] = Function::vectorize(zero,zero);
  functionMap[p->ID()] = zero;

  backgroundFlow->projectOntoMesh(functionMap);

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

  // // stress equation
  bf->addTerm( sigma1, tau1 );
  bf->addTerm( sigma2, tau2 );
  bf->addTerm( u1, tau1->div() );
  bf->addTerm( u2, tau2->div() );
  bf->addTerm( -u1hat, tau1->dot_normal() );
  bf->addTerm( -u2hat, tau2->dot_normal() );

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

//.........这里部分代码省略.........
  {
    bf->addTerm( u, invDt*v );
    rhs->addTerm( u_prev_time * invDt * v );
  }

  ////////////////////   SPECIFY RHS   ///////////////////////
  FunctionPtr f = Teuchos::rcp( new ConstantScalarFunction(0.0) );
  rhs->addTerm( f * v ); // obviously, with f = 0 adding this term is not necessary!

  ////////////////////   DEFINE INNER PRODUCT(S)   ///////////////////////
  IPPtr ip = bf->graphNorm();
  // ip->addTerm(v);
  // ip->addTerm(beta*v->grad());

  ////////////////////   CREATE BCs   ///////////////////////
  Teuchos::RCP<BCEasy> bc = Teuchos::rcp( new BCEasy );
  SpatialFilterPtr lBoundary = Teuchos::rcp( new LeftBoundary );
  FunctionPtr u1 = Teuchos::rcp( new InletBC );
  bc->addDirichlet(beta_n_u_hat, lBoundary, -u1);

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

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

  // ==================== SET INITIAL GUESS ==========================
  double u_free = 0.0;
  map<int, Teuchos::RCP<Function> > functionMap;
  // functionMap[u->ID()]      = Teuchos::rcp( new ConInletBC
  functionMap[u->ID()]      = Teuchos::rcp( new InletBC );

  prevTimeFlow->projectOntoMesh(functionMap);

  ////////////////////   SOLVE & REFINE   ///////////////////////
  if (enforceLocalConservation)
  {
    if (steady)
    {
      FunctionPtr zero = Teuchos::rcp( new ConstantScalarFunction(0.0) );
      solution->lagrangeConstraints()->addConstraint(beta_n_u_hat == zero);
    }
    else
    {
      // FunctionPtr parity = Teuchos::rcp<Function>( new SideParityFunction );
      // LinearTermPtr conservedQuantity = Teuchos::rcp<LinearTerm>( new LinearTerm(parity, beta_n_u_minus_sigma_n) );
      LinearTermPtr conservedQuantity = Teuchos::rcp<LinearTerm>( new LinearTerm(1.0, beta_n_u_hat) );
      LinearTermPtr sourcePart = Teuchos::rcp<LinearTerm>( new LinearTerm(invDt, u) );
      conservedQuantity->addTerm(sourcePart, true);
      solution->lagrangeConstraints()->addConstraint(conservedQuantity == u_prev_time * invDt);
    }
  }

  double energyThreshold = 0.2; // for mesh refinements
  RefinementStrategy refinementStrategy( solution, energyThreshold );
  VTKExporter exporter(solution, mesh, varFactory);

  for (int refIndex=0; refIndex<=numRefs; refIndex++)
  {
    if (steady)
    {
      solution->solve(false);

      if (commRank == 0)
      {

//.........这里部分代码省略.........
    e2[1] = 1;

    FunctionPtr u_prev = Teuchos::rcp( new PreviousSolutionFunction(backgroundFlow, u) );
    FunctionPtr beta = e1 * u_prev + Teuchos::rcp( new ConstantVectorFunction( e2 ) );

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

    // tau parts:
    // 1/eps (sigma, tau)_K + (u, div tau)_K - (u_hat, tau_n)_dK
    bf->addTerm(sigma1 / epsilon, tau->x());
    bf->addTerm(sigma2 / epsilon, tau->y());
    bf->addTerm(u, tau->div());
    bf->addTerm( - uhat, tau->dot_normal() );

    // v:
    // (sigma, grad v)_K - (sigma_hat_n, v)_dK - (u, beta dot grad v) + (u_hat * n dot beta, v)_dK
    bf->addTerm( sigma1, v->dx() );
    bf->addTerm( sigma2, v->dy() );
    bf->addTerm( -u, beta * v->grad());
    bf->addTerm( beta_n_u_minus_sigma_hat, v);

    // ==================== SET INITIAL GUESS ==========================
    mesh->registerSolution(backgroundFlow);
    FunctionPtr zero = Teuchos::rcp( new ConstantScalarFunction(0.0) );
    FunctionPtr u0 = Teuchos::rcp( new U0 );

    map<int, Teuchos::RCP<Function> > functionMap;
    functionMap[u->ID()] = u0;
    functionMap[sigma1->ID()] = zero;
    functionMap[sigma2->ID()] = zero;

    backgroundFlow->projectOntoMesh(functionMap);
    // ==================== END SET INITIAL GUESS ==========================

    ////////////////////////////////////////////////////////////////////
    // DEFINE INNER PRODUCT
    ////////////////////////////////////////////////////////////////////
    // function to scale the squared guy by epsilon/h
    FunctionPtr epsilonOverHScaling = Teuchos::rcp( new EpsilonScaling(epsilon) );
    IPPtr ip = Teuchos::rcp( new IP );
    ip->addTerm( epsilonOverHScaling * (1.0/sqrt(epsilon))* tau);
    ip->addTerm( tau->div());
    //  ip->addTerm( epsilonOverHScaling * v );
    ip->addTerm( v );
    ip->addTerm( sqrt(epsilon) * v->grad() );
    ip->addTerm(v->grad());
    //  ip->addTerm( beta * v->grad() );

    ////////////////////////////////////////////////////////////////////
    // DEFINE RHS
    ////////////////////////////////////////////////////////////////////
    RHSPtr rhs = RHS::rhs();
    FunctionPtr u_prev_squared_div2 = 0.5 * u_prev * u_prev;

    rhs->addTerm((e1 * u_prev_squared_div2 + e2 * u_prev) * v->grad() - u_prev * tau->div());

    ////////////////////////////////////////////////////////////////////
    // DEFINE DIRICHLET BC
    ////////////////////////////////////////////////////////////////////
    FunctionPtr n = Teuchos::rcp( new UnitNormalFunction );
    SpatialFilterPtr outflowBoundary = Teuchos::rcp( new TopBoundary);
    SpatialFilterPtr inflowBoundary = Teuchos::rcp( new NegatedSpatialFilter(outflowBoundary) );
    BCPtr inflowBC = BC::bc();
    FunctionPtr u0_squared_div_2 = 0.5 * u0 * u0;

int main(int argc, char *argv[])
{
#ifdef HAVE_MPI
  Teuchos::GlobalMPISession mpiSession(&argc, &argv,0);
  int rank=mpiSession.getRank();
  int numProcs=mpiSession.getNProc();
#else
  int rank = 0;
  int numProcs = 1;
#endif
  int polyOrder = 3;
  int pToAdd = 2; // for tests

  // define our manufactured solution or problem bilinear form:
  bool useTriangles = false;

  FieldContainer<double> meshPoints(4,2);

  meshPoints(0,0) = 0.0; // x1
  meshPoints(0,1) = 0.0; // y1
  meshPoints(1,0) = 1.0;
  meshPoints(1,1) = 0.0;
  meshPoints(2,0) = 1.0;
  meshPoints(2,1) = 1.0;
  meshPoints(3,0) = 0.0;
  meshPoints(3,1) = 1.0;

  int H1Order = polyOrder + 1;
  int horizontalCells = 4, verticalCells = 4;

  double energyThreshold = 0.2; // for mesh refinements
  double nonlinearStepSize = 0.5;
  double nonlinearRelativeEnergyTolerance = 1e-8; // used to determine convergence of the nonlinear solution

  ////////////////////////////////////////////////////////////////////
  // DEFINE VARIABLES
  ////////////////////////////////////////////////////////////////////

  // new-style bilinear form definition
  VarFactory varFactory;
  VarPtr fhat = varFactory.fluxVar("\\widehat{f}");
  VarPtr u = varFactory.fieldVar("u");

  VarPtr v = varFactory.testVar("v",HGRAD);
  BFPtr bf = Teuchos::rcp( new BF(varFactory) ); // initialize bilinear form

  ////////////////////////////////////////////////////////////////////
  // CREATE MESH
  ////////////////////////////////////////////////////////////////////

  // create a pointer to a new mesh:
  Teuchos::RCP<Mesh> mesh = Mesh::buildQuadMesh(meshPoints, horizontalCells,
                            verticalCells, bf, H1Order,
                            H1Order+pToAdd, useTriangles);
  mesh->setPartitionPolicy(Teuchos::rcp(new ZoltanMeshPartitionPolicy("HSFC")));

  ////////////////////////////////////////////////////////////////////
  // INITIALIZE BACKGROUND FLOW FUNCTIONS
  ////////////////////////////////////////////////////////////////////
  BCPtr nullBC = Teuchos::rcp((BC*)NULL);
  RHSPtr nullRHS = Teuchos::rcp((RHS*)NULL);
  IPPtr nullIP = Teuchos::rcp((IP*)NULL);
  SolutionPtr backgroundFlow = Teuchos::rcp(new Solution(mesh, nullBC,
                               nullRHS, nullIP) );

  vector<double> e1(2); // (1,0)
  e1[0] = 1;
  vector<double> e2(2); // (0,1)
  e2[1] = 1;

  FunctionPtr u_prev = Teuchos::rcp( new PreviousSolutionFunction(backgroundFlow, u) );
  FunctionPtr beta = e1 * u_prev + Teuchos::rcp( new ConstantVectorFunction( e2 ) );

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

  // v:
  // (sigma, grad v)_K - (sigma_hat_n, v)_dK - (u, beta dot grad v) + (u_hat * n dot beta, v)_dK
  bf->addTerm( -u, beta * v->grad());
  bf->addTerm( fhat, v);

  // ==================== SET INITIAL GUESS ==========================
  mesh->registerSolution(backgroundFlow);
  FunctionPtr zero = Teuchos::rcp( new ConstantScalarFunction(0.0) );
  FunctionPtr u0 = Teuchos::rcp( new U0 );

  map<int, Teuchos::RCP<Function> > functionMap;
  functionMap[u->ID()] = u0;

  backgroundFlow->projectOntoMesh(functionMap);
  // ==================== END SET INITIAL GUESS ==========================

  ////////////////////////////////////////////////////////////////////
  // DEFINE INNER PRODUCT
  ////////////////////////////////////////////////////////////////////
  IPPtr ip = Teuchos::rcp( new IP );
  ip->addTerm( v );
  ip->addTerm( beta * v->grad() );

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