C++ VarFactory类代码示例

bool LobattoBasisTests::testLobattoValues() {  
  bool success = true;
  
  FunctionPtr x = Function::xn(1);
  vector< FunctionPtr > lobattoFunctionsExpected;
  // Pavel Solin's first two Lobatto functions:
//  lobattoFunctionsExpected.push_back( (1 - x) / 2 );
//  lobattoFunctionsExpected.push_back( (1 + x) / 2 );
  // Demkowicz's:
  lobattoFunctionsExpected.push_back( Function::constant(1.0) );
  lobattoFunctionsExpected.push_back( x );
  
  lobattoFunctionsExpected.push_back( (x*x - 1) / 2);
  lobattoFunctionsExpected.push_back( (x*x - 1) * x / 2);
  lobattoFunctionsExpected.push_back( (x*x - 1) * (5 * x * x - 1) / 8);
  lobattoFunctionsExpected.push_back( (x*x - 1) * (7 * x * x - 3) * x / 8);
  
  vector< FunctionPtr > lobattoFunctions;
  
  bool conformingFalse = false;
  
  int n_max = 4;
  for (int n=0; n<=n_max; n++) {
    lobattoFunctions.push_back( Teuchos::rcp( new LobattoFunction<>(n,conformingFalse,false) ) );
  }
  
  VarFactory varFactory;
  varFactory.testVar("v", HGRAD);
  varFactory.fieldVar("u", L2);
  
  BFPtr bf = Teuchos::rcp( new BF(varFactory) );
  MeshPtr mesh = MeshFactory::quadMesh(bf, n_max);
  
  BasisCachePtr basisCache = BasisCache::basisCacheForCell(mesh, 0);
  
  double tol = 1e-12;
  for (int n=0; n<=n_max; n++) {
    if (! lobattoFunctions[n]->equals(lobattoFunctionsExpected[n], basisCache, tol) ) {
      cout << "Lobatto function " << n << " does not match expected.\n";
      success = false;
    }
  }
  
  return success;
}

bool LobattoBasisTests::testLobattoDerivativeValues() {
  bool success = true;
  FunctionPtr x = Function::xn(1);
  vector< FunctionPtr > legendreFunctions; // manually specified
  vector< FunctionPtr > lobattoDerivatives;
  
  FunctionPtr one = Function::constant(1);
  legendreFunctions.push_back(one);
  legendreFunctions.push_back(x);
  legendreFunctions.push_back(0.5 * (3 * x * x - 1) );
  legendreFunctions.push_back( (5 * x * x * x - 3 * x) / 2);
  legendreFunctions.push_back( (1.0 / 8.0) * (35 * x * x * x * x - 30 * x * x + 3) );
  legendreFunctions.push_back( (1.0 / 8.0) * (63 * x * x * x * x * x - 70 * x * x * x + 15 * x) );
  
  bool conformingFalse = false;
  
  int n_max = 4;
  for (int n=0; n<=n_max+1; n++) {
    lobattoDerivatives.push_back( Teuchos::rcp( new LobattoFunction<>(n,conformingFalse,true) ) );
  }
  
  VarFactory varFactory;
  varFactory.testVar("v", HGRAD);
  varFactory.fieldVar("u", L2);
  
  BFPtr bf = Teuchos::rcp( new BF(varFactory) );
  MeshPtr mesh = MeshFactory::quadMesh(bf, n_max);
  
  BasisCachePtr basisCache = BasisCache::basisCacheForCell(mesh, 0);
  
  double tol = 1e-8;
  for (int n=0; n<=n_max; n++) {
    if (! legendreFunctions[n]->equals(lobattoDerivatives[n+1], basisCache, tol) ) {
      cout << "Legendre function " << n << " != Lobatto function " << n + 1 << " derivative.\n";
      cout << "L_" << n << "(0.5) = " << Function::evaluate(legendreFunctions[n],0.5) << endl;
      cout << "l'_" << n+1 << "(0.5) = " << Function::evaluate(lobattoDerivatives[n+1],0.5) << endl;
      success = false;
    }
  }
  
  return success;
}

bool LobattoBasisTests::testLegendreValues() {  
  bool success = true;
  
  FunctionPtr x = Function::xn(1);
  vector< FunctionPtr > legendreFunctionsExpected;
  legendreFunctionsExpected.push_back( Function::constant(1.0) );
  legendreFunctionsExpected.push_back( x );
  legendreFunctionsExpected.push_back( (3 * x*x - 1) / 2);
  legendreFunctionsExpected.push_back( (5 * x*x*x - 3*x) / 2);
  legendreFunctionsExpected.push_back( (35 * x * x * x * x - 30 * x * x + 3) / 8);
  legendreFunctionsExpected.push_back( (63 * x * x * x * x * x - 70 * x * x * x + 15 * x) / 8);
  
  vector< FunctionPtr > legendreFunctions;
  
  int n_max = 4;
  for (int n=0; n<=n_max; n++) {
    legendreFunctions.push_back( Teuchos::rcp( new LegendreFunction(n) ) );
  }
  
  VarFactory varFactory;
  varFactory.testVar("v", HGRAD);
  varFactory.fieldVar("u", L2);
  
  BFPtr bf = Teuchos::rcp( new BF(varFactory) );
  MeshPtr mesh = MeshFactory::quadMesh(bf, n_max);
  
  BasisCachePtr basisCache = BasisCache::basisCacheForCell(mesh, 0);
  
  double tol = 1e-8; // relax to make sure that failure isn't just roundoff
  for (int n=0; n<=n_max; n++) {
    if (! legendreFunctions[n]->equals(legendreFunctionsExpected[n], basisCache, tol) ) {
      cout << "Legendre function " << n << " does not match expected.\n";
      success = false;
    }
  }
  
  return success;
}

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
  ////////////////////   DECLARE VARIABLES   ///////////////////////
  // define test variables
  VarFactory varFactory;
  VarPtr tau = varFactory.testVar("\\tau", HDIV);
  VarPtr v = varFactory.testVar("v", HGRAD);

  // define 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");

  vector<double> beta_const;
  double c = sqrt(1.25);
  beta_const.push_back(1.0/c);
  beta_const.push_back(.5/c);
//  FunctionPtr beta = Teuchos::rcp(new Beta());

  double eps = 1e-3;

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

  BFPtr confusionBF = Teuchos::rcp( new BF(varFactory) );
  // tau terms:
  confusionBF->addTerm(sigma1 / eps, tau->x());
  confusionBF->addTerm(sigma2 / eps, 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_const * u, - v->grad() );
  confusionBF->addTerm( beta_n_u_minus_sigma_n, v);

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

  // quasi-optimal norm
  IPPtr qoptIP = Teuchos::rcp(new IP);
  qoptIP->addTerm( v );
  qoptIP->addTerm( tau / eps + v->grad() );
  qoptIP->addTerm( beta_const * v->grad() - tau->div() );

  // robust test norm
  IPPtr robIP = Teuchos::rcp(new IP);
  FunctionPtr ip_scaling = Teuchos::rcp( new EpsilonScaling(eps) );
  if (enforceLocalConservation)
  {
    robIP->addZeroMeanTerm( v );
  }
  else
  {
    robIP->addTerm( ip_scaling * v );
  }
  robIP->addTerm( sqrt(eps) * v->grad() );

  bool useNewBC = false;
  FunctionPtr weight = Teuchos::rcp( new SqrtWeight(eps) );
  if (useNewBC)
  {
    robIP->addTerm( beta_const * v->grad() );
    robIP->addTerm( tau->div() );
    robIP->addTerm( ip_scaling/sqrt(eps) * tau );
  }
  else
  {
    robIP->addTerm( weight * beta_const * v->grad() );
    robIP->addTerm( weight * tau->div() );
    robIP->addTerm( weight * ip_scaling/sqrt(eps) * tau );
  }


  ////////////////////   SPECIFY RHS   ///////////////////////
  FunctionPtr zero = Teuchos::rcp( new ConstantScalarFunction(0.0) );
  Teuchos::RCP<RHSEasy> rhs = Teuchos::rcp( new RHSEasy );
  FunctionPtr f = zero;
  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 inflowBoundary = Teuchos::rcp( new InflowSquareBoundary );
  SpatialFilterPtr outflowBoundary = Teuchos::rcp( new OutflowSquareBoundary );

  FunctionPtr u0 = Teuchos::rcp( new U0 );
  FunctionPtr n = Teuchos::rcp( new UnitNormalFunction );
  bc->addDirichlet(uhat, outflowBoundary, zero);
  if (useNewBC)
  {
    bc->addDirichlet(beta_n_u_minus_sigma_n, inflowBoundary, beta_const*n*u0);
//.........这里部分代码省略.........

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

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

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   ///////////////////////
//.........这里部分代码省略.........

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 numRefs = args.Input<int>("--numRefs","num adaptive refinements",0);
  int numPreRefs = args.Input<int>("--numPreRefs","num preemptive adaptive refinements",0);
  int order = args.Input<int>("--order","order of approximation",2);
  double eps = args.Input<double>("--epsilon","diffusion parameter",1e-2);
  double energyThreshold = args.Input<double>("-energyThreshold","energy thresh for adaptivity", .5);
  double rampHeight = args.Input<double>("--rampHeight","ramp height at x = 2", 0.0);
  double ipSwitch = args.Input<double>("--ipSwitch","point at which to switch to graph norm", 0.0); // default to 0 to remain on robust norm
  bool useAnisotropy = args.Input<bool>("--useAnisotropy","aniso flag ", false);

  int H1Order = order+1; 
  int pToAdd = args.Input<int>("--pToAdd","test space enrichment", 2);

  FunctionPtr zero = Function::constant(0.0);
  FunctionPtr one = Function::constant(1.0);
  FunctionPtr n = Teuchos::rcp( new UnitNormalFunction );
  vector<double> e1,e2;
  e1.push_back(1.0);e1.push_back(0.0);
  e2.push_back(0.0);e2.push_back(1.0);

  ////////////////////   DECLARE VARIABLES   ///////////////////////
  // define test variables
  VarFactory varFactory; 
  VarPtr tau = varFactory.testVar("\\tau", HDIV);
  VarPtr v = varFactory.testVar("v", HGRAD);
  
  // define 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");

  vector<double> beta;
  beta.push_back(1.0);
  beta.push_back(0.0);
  
  ////////////////////   DEFINE BILINEAR FORM   ///////////////////////

  BFPtr confusionBF = Teuchos::rcp( new BF(varFactory) );
  // tau terms:
  confusionBF->addTerm(sigma1 / eps, tau->x());
  confusionBF->addTerm(sigma2 / eps, 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( -u, beta * v->grad() );
  confusionBF->addTerm( beta_n_u_minus_sigma_n, v);

  // first order term with magnitude alpha
  double alpha = 0.0;
  //  confusionBF->addTerm(alpha * u, v);

  ////////////////////   BUILD MESH   ///////////////////////


  // create a pointer to a new mesh:
  Teuchos::RCP<Mesh> mesh = MeshUtilities::buildUnitQuadMesh(nCells,confusionBF, H1Order, H1Order+pToAdd);
  mesh->setPartitionPolicy(Teuchos::rcp(new ZoltanMeshPartitionPolicy("HSFC")));  
  MeshInfo meshInfo(mesh); // gets info like cell measure, etc

  ////////////////////   DEFINE INNER PRODUCT(S)   ///////////////////////
  IPPtr ip = Teuchos::rcp(new IP);

  /*
   // robust test norm
  FunctionPtr C_h = Teuchos::rcp( new EpsilonScaling(eps) );  
  FunctionPtr invH = Teuchos::rcp(new InvHScaling);
  FunctionPtr invSqrtH = Teuchos::rcp(new InvSqrtHScaling);
  FunctionPtr sqrtH = Teuchos::rcp(new SqrtHScaling);
  FunctionPtr hSwitch = Teuchos::rcp(new HSwitch(ipSwitch,mesh));
  ip->addTerm(hSwitch*sqrt(eps) * v->grad() );
  ip->addTerm(hSwitch*beta * v->grad() );
  ip->addTerm(hSwitch*tau->div() );
  
  // graph norm
  ip->addTerm( (one-hSwitch)*((1.0/eps) * tau + v->grad()));
  ip->addTerm( (one-hSwitch)*(beta * v->grad() - tau->div()));

  // regularizing terms
  ip->addTerm(C_h/sqrt(eps) * tau );    
  ip->addTerm(invSqrtH*v);
  */

   // robust test norm
  IPPtr robIP = Teuchos::rcp(new IP);
//.........这里部分代码省略.........

int main(int argc, char *argv[])
{
  int rank = 0;
#ifdef HAVE_MPI
  // TODO: figure out the right thing to do here...
  // may want to modify argc and argv before we make the following call:
  Teuchos::GlobalMPISession mpiSession(&argc, &argv,0);
  rank=mpiSession.getRank();
#else
#endif
  bool useLineSearch = false;

  int pToAdd = 2; // for optimal test function approximation
  int pToAddForStreamFunction = 2;
  double nonlinearStepSize = 1.0;
  double dt = 0.5;
  double nonlinearRelativeEnergyTolerance = 0.015; // used to determine convergence of the nonlinear solution
  //  double nonlinearRelativeEnergyTolerance = 0.15; // used to determine convergence of the nonlinear solution
  double eps = 1.0/64.0; // width of ramp up to 1.0 for top BC;  eps == 0 ==> soln not in H1
  // epsilon above is chosen to match our initial 16x16 mesh, to avoid quadrature errors.
  //  double eps = 0.0; // John Evans's problem: not in H^1
  bool enforceLocalConservation = false;
  bool enforceOneIrregularity = true;
  bool reportPerCellErrors  = true;
  bool useMumps = true;

  int horizontalCells, verticalCells;

  int maxIters = 50; // for nonlinear steps

  vector<double> ReValues;

  // usage: polyOrder [numRefinements]
  // parse args:
  if (argc < 6)
  {
    cout << "Usage: NavierStokesCavityFlowContinuationFixedMesh fieldPolyOrder hCells vCells energyErrorGoal Re0 [Re1 ...]\n";
    return -1;
  }
  int polyOrder = atoi(argv[1]);
  horizontalCells = atoi(argv[2]);
  verticalCells = atoi(argv[3]);
  double energyErrorGoal = atof(argv[4]);
  for (int i=5; i<argc; i++)
  {
    ReValues.push_back(atof(argv[i]));
  }
  if (rank == 0)
  {
    cout << "L^2 order: " << polyOrder << endl;
    cout << "initial mesh size: " << horizontalCells << " x " << verticalCells << endl;
    cout << "energy error goal: " << energyErrorGoal << endl;
    cout << "Reynolds number values for continuation:\n";
    for (int i=0; i<ReValues.size(); i++)
    {
      cout << ReValues[i] << ", ";
    }
    cout << endl;
  }

  FieldContainer<double> quadPoints(4,2);

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

  // define meshes:
  int H1Order = polyOrder + 1;
  bool useTriangles = false;
  bool meshHasTriangles = useTriangles;

  double minL2Increment = 1e-8;

  // get variable definitions:
  VarFactory varFactory = VGPStokesFormulation::vgpVarFactory();
  u1 = varFactory.fieldVar(VGP_U1_S);
  u2 = varFactory.fieldVar(VGP_U2_S);
  sigma11 = varFactory.fieldVar(VGP_SIGMA11_S);
  sigma12 = varFactory.fieldVar(VGP_SIGMA12_S);
  sigma21 = varFactory.fieldVar(VGP_SIGMA21_S);
  sigma22 = varFactory.fieldVar(VGP_SIGMA22_S);
  p = varFactory.fieldVar(VGP_P_S);

  u1hat = varFactory.traceVar(VGP_U1HAT_S);
  u2hat = varFactory.traceVar(VGP_U2HAT_S);
  t1n = varFactory.fluxVar(VGP_T1HAT_S);
  t2n = varFactory.fluxVar(VGP_T2HAT_S);

  v1 = varFactory.testVar(VGP_V1_S, HGRAD);
  v2 = varFactory.testVar(VGP_V2_S, HGRAD);
  tau1 = varFactory.testVar(VGP_TAU1_S, HDIV);
  tau2 = varFactory.testVar(VGP_TAU2_S, HDIV);
  q = varFactory.testVar(VGP_Q_S, HGRAD);

  FunctionPtr u1_0 = Teuchos::rcp( new U1_0(eps) );
//.........这里部分代码省略.........

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");
  bool enforceLocalConservation = args.Input<bool>("--conserve", "enforce local conservation");
  bool steady = args.Input<bool>("--steady", "run steady rather than transient");

  // Optional arguments (have defaults)
  double dt = args.Input("--dt", "time step", 0.25);
  int numTimeSteps = args.Input("--nt", "number of time steps", 20);
  halfWidth = args.Input("--halfWidth", "half width of inlet profile", 1.0);
  args.Process();

  ////////////////////   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");

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

  ////////////////////   BUILD MESH   ///////////////////////
  BFPtr bf = Teuchos::rcp( new BF(varFactory) );
  // define nodes for mesh
  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 = 8, verticalCells = 8;

  // create a pointer to a new mesh:
  Teuchos::RCP<Mesh> mesh = Mesh::buildQuadMesh(meshBoundary, horizontalCells, verticalCells,
                            bf, 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:
  bf->addTerm( beta * u, - v->grad() );
  bf->addTerm( beta_n_u_hat, v);

  if (!steady)
  {
    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);
//.........这里部分代码省略.........

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
  ////////////////////   DECLARE VARIABLES   ///////////////////////
  // define test variables
  VarFactory varFactory;
  VarPtr tau = varFactory.testVar("\\tau", HDIV);
  VarPtr v = varFactory.testVar("v", HGRAD);

  // define 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");

  vector<double> beta_const;
  beta_const.push_back(1.0);
  beta_const.push_back(0.0);
//  FunctionPtr beta = Teuchos::rcp(new Beta());

  double eps = 1e-2;

  ////////////////////   DEFINE BILINEAR FORM   ///////////////////////
  BFPtr confusionBF = Teuchos::rcp( new BF(varFactory) );
  // tau terms:
  confusionBF->addTerm(sigma1 / eps, tau->x());
  confusionBF->addTerm(sigma2 / eps, 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_const * u, - v->grad() );
  confusionBF->addTerm( beta_n_u_minus_sigma_n, v);

  ////////////////////   DEFINE INNER PRODUCT(S)   ///////////////////////
  // mathematician's norm
  IPPtr mathIP = Teuchos::rcp(new IP());
  mathIP->addTerm(tau);
  mathIP->addTerm(tau->div());

  mathIP->addTerm(v);
  mathIP->addTerm(v->grad());

  // quasi-optimal norm
  IPPtr qoptIP = Teuchos::rcp(new IP);
  qoptIP->addTerm( v );
  qoptIP->addTerm( tau / eps + v->grad() );
  qoptIP->addTerm( beta_const * v->grad() - tau->div() );

  // robust test norm
  IPPtr robIP = Teuchos::rcp(new IP);
  FunctionPtr ip_scaling = Teuchos::rcp( new EpsilonScaling(eps) );
  if (enforceLocalConservation)
  {
    robIP->addZeroMeanTerm( v );
  }
  else
  {
    robIP->addTerm( ip_scaling * v );
  }

  robIP->addTerm( sqrt(eps) * v->grad() );
  robIP->addTerm( beta_const * v->grad() );
  robIP->addTerm( tau->div() );
  robIP->addTerm( ip_scaling/sqrt(eps) * tau );

  ////////////////////   SPECIFY RHS   ///////////////////////
  FunctionPtr zero = Teuchos::rcp( new ConstantScalarFunction(0.0) );
  Teuchos::RCP<RHSEasy> rhs = Teuchos::rcp( new RHSEasy );
  FunctionPtr f = zero;
  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 inflowBoundary = Teuchos::rcp( new InflowSquareBoundary );
  //  SpatialFilterPtr outflowBoundary = Teuchos::rcp( new OutflowSquareBoundary );

  SpatialFilterPtr inflowTop = Teuchos::rcp(new InflowLshapeTop);
  SpatialFilterPtr inflowBot = Teuchos::rcp(new InflowLshapeBottom);
  SpatialFilterPtr LshapeBot1 = Teuchos::rcp(new LshapeBottom1);
  SpatialFilterPtr LshapeBot2 = Teuchos::rcp(new LshapeBottom2);
  SpatialFilterPtr Top = Teuchos::rcp(new LshapeTop);
  SpatialFilterPtr Out = Teuchos::rcp(new LshapeOutflow);

  FunctionPtr u0 = Teuchos::rcp( new U0 );
  bc->addDirichlet(uhat, LshapeBot1, u0);
  bc->addDirichlet(uhat, LshapeBot2, u0);
  bc->addDirichlet(uhat, Top, u0);
  bc->addDirichlet(uhat, Out, 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 = 2;

    // define our manufactured solution or problem bilinear form:
    double epsilon = 1e-3;
    bool useTriangles = false;

    int pToAdd = 2;
    int nCells = 2;
    if ( argc > 1)
    {
        nCells = atoi(argv[1]);
        if (rank==0)
        {
            cout << "numCells = " << nCells << endl;
        }
    }
    int numSteps = 20;
    if ( argc > 2)
    {
        numSteps = atoi(argv[2]);
        if (rank==0)
        {
            cout << "num NR steps = " << numSteps << endl;
        }
    }
    int useHessian = 0; // defaults to "not use"
    if ( argc > 3)
    {
        useHessian = atoi(argv[3]);
        if (rank==0)
        {
            cout << "useHessian = " << useHessian << endl;
        }
    }

    int thresh = numSteps; // threshhold for when to apply linesearch/hessian
    if ( argc > 4)
    {
        thresh = atoi(argv[4]);
        if (rank==0)
        {
            cout << "thresh = " << thresh << endl;
        }
    }

    int H1Order = polyOrder + 1;

    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 uhat = varFactory.traceVar("\\widehat{u}");
    VarPtr beta_n_u_minus_sigma_hat = varFactory.fluxVar("\\widehat{\\beta_n u - \\sigma_n}");
    VarPtr u = varFactory.fieldVar("u");
    VarPtr sigma1 = varFactory.fieldVar("\\sigma_1");
    VarPtr sigma2 = varFactory.fieldVar("\\sigma_2");

    VarPtr tau = varFactory.testVar("\\tau",HDIV);
    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);
    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) );
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    vector<FunctionPtr> functions;
    functions.push_back(function);
    functions.push_back(vect);
    vector<string> functionNames;
    functionNames.push_back("function");
    functionNames.push_back("vect");
    vector<FunctionPtr> bdrfunctions;
    bdrfunctions.push_back(fbdr);
    bdrfunctions.push_back(fbdr);
    vector<string> bdrfunctionNames;
    bdrfunctionNames.push_back("bdr1");
    bdrfunctionNames.push_back("bdr2");

    map<int, int> cellIDToNum1DPts;
    cellIDToNum1DPts[1] = 4;

    {
      XDMFExporter exporter(meshTopology, "Grid2D", false);
      // exporter.exportFunction(function, "function2", 0, 10);
      // exporter.exportFunction(vect, "vect2", 1, 10, cellIDToNum1DPts);
      // exporter.exportFunction(fbdr, "boundary2", 0);
      exporter.exportFunction(functions, functionNames, 1, 10);
    }
    {
      XDMFExporter exporter(meshTopology, "BdrGrid2D", false);
      // exporter.exportFunction(function, "function2", 0, 10);
      // exporter.exportFunction(vect, "vect2", 1, 10, cellIDToNum1DPts);
      // exporter.exportFunction(fbdr, "boundary2", 0);
      exporter.exportFunction(bdrfunctions, bdrfunctionNames, 1, 10);
    }

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

    // define trial variables
    VarPtr uhat = varFactory.traceVar("uhat");
    VarPtr fhat = varFactory.fluxVar("fhat");
    VarPtr u = varFactory.fieldVar("u");
    VarPtr sigma = varFactory.fieldVar("sigma", VECTOR_L2);

    ////////////////////   DEFINE BILINEAR FORM   ///////////////////////
    BFPtr bf = Teuchos::rcp( new BF(varFactory) );
    // tau terms:
    bf->addTerm(sigma, tau);
    bf->addTerm(u, tau->div());
    bf->addTerm(-uhat, tau->dot_normal());

    // v terms:
    bf->addTerm( sigma, v->grad() );
    bf->addTerm( fhat, v);

    ////////////////////   BUILD MESH   ///////////////////////
    int H1Order = 4, pToAdd = 2;
    Teuchos::RCP<Mesh> mesh = Teuchos::rcp( new Mesh (meshTopology, bf, H1Order, pToAdd) );

    ////////////////////   DEFINE INNER PRODUCT(S)   ///////////////////////
    IPPtr ip = bf->graphNorm();

    ////////////////////   SPECIFY RHS   ///////////////////////
    RHSPtr rhs = RHS::rhs();
    // Teuchos::RCP<RHS> rhs = Teuchos::rcp( new RHS );
    FunctionPtr one = Function::constant(1.0);
    rhs->addTerm( one * v );

void FunctionTests::SetUp()
{
  ////////////////////   DECLARE VARIABLES   ///////////////////////
  // define test variables
  VarFactory varFactory; 
  VarPtr tau = varFactory.testVar("\\tau", HDIV);
  VarPtr v = varFactory.testVar("v", HGRAD);
  
  // define 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");
  
  vector<double> beta_const;
  beta_const.push_back(2.0);
  beta_const.push_back(1.0);
  
  double eps = 1e-2;
  
  // standard confusion bilinear form
  _confusionBF = Teuchos::rcp( new BF(varFactory) );
  // tau terms:
  _confusionBF->addTerm(sigma1 / eps, tau->x());
  _confusionBF->addTerm(sigma2 / eps, 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_const * u, - v->grad() );
  _confusionBF->addTerm( beta_n_u_minus_sigma_n, v);
    
  ////////////////////   BUILD MESH   ///////////////////////
  // define nodes for mesh
  FieldContainer<double> quadPoints(4,2);
  
  quadPoints(0,0) = -1.0; // x1
  quadPoints(0,1) = -1.0; // y1
  quadPoints(1,0) = 1.0;
  quadPoints(1,1) = -1.0;
  quadPoints(2,0) = 1.0;
  quadPoints(2,1) = 1.0;
  quadPoints(3,0) = -1.0;
  quadPoints(3,1) = 1.0;
  
  int H1Order = 1, pToAdd = 0;
  int horizontalCells = 1, verticalCells = 1;
  
  // create a pointer to a new mesh:
  _spectralConfusionMesh = MeshFactory::buildQuadMesh(quadPoints, horizontalCells, verticalCells,
                                               _confusionBF, H1Order, H1Order+pToAdd);
  
  // some 2D test points:
  // setup test points:
  static const int NUM_POINTS_1D = 10;
  double x[NUM_POINTS_1D] = {-1.0,-0.8,-0.6,-.4,-.2,0,0.2,0.4,0.6,0.8};
  double y[NUM_POINTS_1D] = {-0.8,-0.6,-.4,-.2,0,0.2,0.4,0.6,0.8,1.0};
  
  _testPoints = FieldContainer<double>(NUM_POINTS_1D*NUM_POINTS_1D,2);
  for (int i=0; i<NUM_POINTS_1D; i++) {
    for (int j=0; j<NUM_POINTS_1D; j++) {
      _testPoints(i*NUM_POINTS_1D + j, 0) = x[i];
      _testPoints(i*NUM_POINTS_1D + j, 1) = y[j];
    }
  }
  
  _elemType = _spectralConfusionMesh->getElement(0)->elementType();
  vector<GlobalIndexType> cellIDs;
  int cellID = 0;
  cellIDs.push_back(cellID);
  _basisCache = Teuchos::rcp( new BasisCache( _elemType, _spectralConfusionMesh ) );
  _basisCache->setRefCellPoints(_testPoints);
  
  _basisCache->setPhysicalCellNodes( _spectralConfusionMesh->physicalCellNodesForCell(cellID), cellIDs, true );
}

bool HConvergenceStudyTests::testBestApproximationErrorComputation() {
  bool success = true;

  bool enrichVelocity = false; // true would be for the "compliant" norm, which isn't working well yet
  
  int minLogElements = 0, maxLogElements = minLogElements;
  int numCells1D = pow(2.0,minLogElements);
  int H1Order = 1;
  int pToAdd = 2;
  
  double tol = 1e-16;
  double Re = 40.0;
  
  VarFactory varFactory = VGPStokesFormulation::vgpVarFactory();
  VarPtr u1_vgp = varFactory.fieldVar(VGP_U1_S);
  VarPtr u2_vgp = varFactory.fieldVar(VGP_U2_S);
  VarPtr sigma11_vgp = varFactory.fieldVar(VGP_SIGMA11_S);
  VarPtr sigma12_vgp = varFactory.fieldVar(VGP_SIGMA12_S);
  VarPtr sigma21_vgp = varFactory.fieldVar(VGP_SIGMA21_S);
  VarPtr sigma22_vgp = varFactory.fieldVar(VGP_SIGMA22_S);
  VarPtr p_vgp = varFactory.fieldVar(VGP_P_S);
  
  VGPStokesFormulation stokesForm(1/Re);
  
  int numCellsFineMesh = 20; // for computing a zero-mean pressure
  int H1OrderFineMesh = 5;
  
  // define Kovasznay domain:
  FieldContainer<double> quadPointsKovasznay(4,2);
  
  // Domain from Evans Hughes for Navier-Stokes:
  quadPointsKovasznay(0,0) =  0.0; // x1
  quadPointsKovasznay(0,1) = -0.5; // y1
  quadPointsKovasznay(1,0) =  1.0;
  quadPointsKovasznay(1,1) = -0.5;
  quadPointsKovasznay(2,0) =  1.0;
  quadPointsKovasznay(2,1) =  0.5;
  quadPointsKovasznay(3,0) =  0.0;
  quadPointsKovasznay(3,1) =  0.5;
  
  FunctionPtr zero = Function::zero();
  bool dontEnhanceFluxes = false;
  VGPNavierStokesProblem zeroProblem = VGPNavierStokesProblem(Re, quadPointsKovasznay,
                                                              numCellsFineMesh, numCellsFineMesh,
                                                              H1OrderFineMesh, pToAdd,
                                                              zero, zero, zero, enrichVelocity, dontEnhanceFluxes);
  
  FunctionPtr u1_exact, u2_exact, p_exact;
  NavierStokesFormulation::setKovasznay(Re, zeroProblem.mesh(), u1_exact, u2_exact, p_exact);
  
  
  VGPNavierStokesProblem problem = VGPNavierStokesProblem(Re,quadPointsKovasznay,
                                                          numCells1D,numCells1D,
                                                          H1Order, pToAdd,
                                                          u1_exact, u2_exact, p_exact, enrichVelocity, dontEnhanceFluxes);

  HConvergenceStudy study(problem.exactSolution(),
                          problem.mesh()->bilinearForm(),
                          problem.exactSolution()->rhs(),
                          problem.backgroundFlow()->bc(),
                          problem.bf()->graphNorm(),
                          minLogElements, maxLogElements,
                          H1Order, pToAdd, false, false, false);
  study.setReportRelativeErrors(false); // we want absolute errors

  Teuchos::RCP<Mesh> mesh = problem.mesh();
  
  int cubatureDegreeEnrichment = 10;
  
  int L2Order = H1Order - 1;
  int meshCubatureDegree = L2Order + H1Order + pToAdd;

  study.setCubatureDegreeForExact(cubatureDegreeEnrichment + meshCubatureDegree);
  
  FunctionPtr f = u1_exact;
  int trialID = u1_vgp->ID();
  {
    double fIntegral = f->integrate(mesh,cubatureDegreeEnrichment);
//    cout << "testBestApproximationErrorComputation: integral of f on whole mesh = " << fIntegral << endl;
    
    double l2ErrorOfAverage = (Function::constant(fIntegral) - f)->l2norm(mesh,cubatureDegreeEnrichment);
//    cout << "testBestApproximationErrorComputation: l2 error of fIntegral: " << l2ErrorOfAverage << endl;
    
    ElementTypePtr elemType = mesh->elementTypes()[0];
    vector<GlobalIndexType> cellIDs = mesh->cellIDsOfTypeGlobal(elemType);
    
    bool testVsTest = false;
    BasisCachePtr basisCache = Teuchos::rcp( new BasisCache(elemType, mesh, testVsTest, cubatureDegreeEnrichment) );
    basisCache->setPhysicalCellNodes(mesh->physicalCellNodesGlobal(elemType), cellIDs, false); // false: no side cache

    FieldContainer<double> projectionValues(cellIDs.size());
    f->integrate(projectionValues, basisCache);
    FieldContainer<double> cellMeasures = basisCache->getCellMeasures();
    
    for (int i=0; i<projectionValues.size(); i++) {
      projectionValues(i) /= cellMeasures(i);
    }
    
    // since we're not worried about the actual solution values at all, just use a single zero solution:
    vector< SolutionPtr > solutions;
//.........这里部分代码省略.........