C++ BCPtr::addDirichlet方法代码示例

int main(int argc, char *argv[])
{
  Teuchos::GlobalMPISession mpiSession(&argc, &argv, NULL); // initialize MPI
  
  Teuchos::CommandLineProcessor cmdp(false,true); // false: don't throw exceptions; true: do return errors for unrecognized options
  
  int numElements = 3;
  double xLeft = 0.0, xRight = 1.0;
  int polyOrder = 1, delta_k = 1;
  
  cmdp.setOption("numElements", &numElements );
  cmdp.setOption("polyOrder", &polyOrder );
  cmdp.setOption("delta_k", &delta_k );
  cmdp.setOption("xLeft", &xLeft );
  cmdp.setOption("xRight", &xRight );
  
  if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL)
  {
#ifdef HAVE_MPI
    MPI_Finalize();
#endif
    return -1;
  }
  int spaceDim = 1;
  bool conformingTraces = true; // conformingTraces argument has no effect in 1D
  PoissonFormulation poissonForm(spaceDim, conformingTraces);
  
  MeshPtr mesh = MeshFactory::intervalMesh(poissonForm.bf(), xLeft, xRight, numElements, polyOrder + 1, delta_k); // 1D equispaced
  
  RHSPtr rhs = RHS::rhs(); // zero RHS
  IPPtr ip = poissonForm.bf()->graphNorm();
  BCPtr bc = BC::bc();
  bc->addDirichlet(poissonForm.phi_hat(), SpatialFilter::allSpace(), Function::zero());
  
  SolutionPtr solution = Solution::solution(poissonForm.bf(), mesh, bc, rhs, ip);
  solution->solve();
  
  GDAMinimumRule* minRule = dynamic_cast<GDAMinimumRule*>(mesh->globalDofAssignment().get());
//  minRule->printGlobalDofInfo();
  
  Teuchos::RCP<Epetra_CrsMatrix> A = solution->getStiffnessMatrix();
  EpetraExt::RowMatrixToMatrixMarketFile("A.dat",*A, NULL, NULL, false);
  
  HDF5Exporter exporter(mesh);
  
  return 0;
}

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

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

  // quasi-optimal norm
  IPPtr qoptIP = Teuchos::rcp(new IP);
  
  if (!useCompliantGraphNorm) {
    qoptIP->addTerm( tau / eps + v->grad() );
    qoptIP->addTerm( beta_const * v->grad() - tau->div() );
    
    qoptIP->addTerm( v );
  } else {
    FunctionPtr h = Teuchos::rcp( new hFunction );
    
    // here, we're aiming at optimality in 1/h^2 |u|^2 + 1/eps^2 |sigma|^2
    
    qoptIP->addTerm( tau + eps * v->grad() );
    qoptIP->addTerm( h * beta_const * v->grad() - tau->div() );
    qoptIP->addTerm(v);
    qoptIP->addTerm(tau);
  }
  
  ////////////////////   SPECIFY RHS   ///////////////////////
  RHSPtr rhs = RHS::rhs();
  FunctionPtr f = Teuchos::rcp( new ConstantScalarFunction(0.0) );
  rhs->addTerm( f * v ); // obviously, with f = 0 adding this term is not necessary!

  ////////////////////   CREATE BCs   ///////////////////////
  BCPtr bc = BC::bc();
  SpatialFilterPtr inflowBoundary = Teuchos::rcp( new InflowSquareBoundary );
  SpatialFilterPtr outflowBoundary = Teuchos::rcp( new OutflowSquareBoundary );
  FunctionPtr u0 = Teuchos::rcp( new U0 );
  bc->addDirichlet(uhat, outflowBoundary, u0);

  bc->addDirichlet(uhat, inflowBoundary, u0);
  
//  Teuchos::RCP<PenaltyConstraints> pc = Teuchos::rcp(new PenaltyConstraints);
//  pc->addConstraint(uhat==u0,inflowBoundary);

  ////////////////////   BUILD MESH   ///////////////////////
  // create a new mesh on a single-cell, unit square domain
  Teuchos::RCP<Mesh> mesh = MeshFactory::quadMeshMinRule(confusionBF, H1Order, delta_k);
  
  ////////////////////   SOLVE & REFINE   ///////////////////////
  Teuchos::RCP<Solution> solution = Teuchos::rcp( new Solution(mesh, bc, rhs, qoptIP) );
//  solution->setFilter(pc);
  
  double energyThreshold = 0.2; // for mesh refinements
  
  bool useRieszRepBasedRefStrategy = true;
  
  if (rank==0) {
    if (useRieszRepBasedRefStrategy) {
      cout << "using RieszRep-based refinement strategy.\n";
    } else {
      cout << "using solution-based refinement strategy.\n";
    }
  }
  Teuchos::RCP<RefinementStrategy> refinementStrategy;
  if (!useRieszRepBasedRefStrategy) {
    refinementStrategy = Teuchos::rcp( new RefinementStrategy( solution, energyThreshold ) );
  } else {
    LinearTermPtr residual = confusionBF->testFunctional(solution) - rhs->linearTerm();
    refinementStrategy = Teuchos::rcp( new RefinementStrategy( mesh, residual, qoptIP, energyThreshold ) );
  }

bool VectorizedBasisTestSuite::testPoisson()
{
  bool success = true;

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

  // define trial variables
  VarPtr uhat = varFactory->traceVar("\\widehat{u}");
  VarPtr sigma_n = varFactory->fluxVar("\\widehat{\\sigma_{n}}");
  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( -sigma_n, v);

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

  ////////////////////   SPECIFY RHS   ///////////////////////
  RHSPtr rhs = RHS::rhs();
  FunctionPtr f = Function::constant(1.0);
  rhs->addTerm( f * v );

  ////////////////////   CREATE BCs   ///////////////////////
  BCPtr bc = BC::bc();
  SpatialFilterPtr boundary = SpatialFilter::allSpace();
  FunctionPtr zero = Function::zero();
  bc->addDirichlet(uhat, boundary, zero);

  ////////////////////   BUILD MESH   ///////////////////////
  int H1Order = 3, pToAdd = 2;
  // define nodes for mesh
  FieldContainer<double> meshBoundary(4,2);

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

  int horizontalCells = 1, verticalCells = 1;

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

  ////////////////////   SOLVE & REFINE   ///////////////////////
  Teuchos::RCP<Solution> solution = Teuchos::rcp( new Solution(mesh, bc, rhs, ip) );
  double energyThreshold = 0.2; // for mesh refinements
  RefinementStrategy refinementStrategy( solution, energyThreshold );
#ifdef USE_VTK
  VTKExporter exporter(solution, mesh, varFactory);
#endif

  for (int refIndex=0; refIndex<=4; refIndex++)
  {
    solution->solve(false);
#ifdef USE_VTK
    // output commented out because it's not properly part of the test.
//    stringstream outfile;
//    outfile << "test_" << refIndex;
//    exporter.exportSolution(outfile.str());
#endif

    if (refIndex < 4)
      refinementStrategy.refine(false); // don't print to console
  }
  return success;
}

//.........这里部分代码省略.........
    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;
    inflowBC->addDirichlet(beta_n_u_minus_sigma_hat,inflowBoundary,
                           ( e1 * u0_squared_div_2 + e2 * u0) * n );

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

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

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

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

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

    Teuchos::RCP<RieszRep> riesz = Teuchos::rcp(new RieszRep(mesh, ip, residual));
    Teuchos::RCP<RieszRep> duRiesz = Teuchos::rcp(new RieszRep(mesh, ip, Bdu));

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

    ////////////////////   CREATE BCs   ///////////////////////
    // Teuchos::RCP<BC> bc = Teuchos::rcp( new BCEasy );
    BCPtr bc = BC::bc();
    FunctionPtr zero = Function::zero();
    SpatialFilterPtr entireBoundary = Teuchos::rcp( new EntireBoundary );
    bc->addDirichlet(uhat, entireBoundary, zero);

    ////////////////////   SOLVE & REFINE   ///////////////////////
    Teuchos::RCP<Solution> solution = Teuchos::rcp( new Solution(mesh, bc, rhs, ip) );
    solution->solve(false);
    RefinementStrategy refinementStrategy( solution, 0.2);

    // Output solution
    FunctionPtr uSoln = Function::solution(u, solution);
    FunctionPtr sigmaSoln = Function::solution(sigma, solution);
    FunctionPtr uhatSoln = Function::solution(uhat, solution);
    FunctionPtr fhatSoln = Function::solution(fhat, solution);
    {
      XDMFExporter exporter(meshTopology, "Poisson", false);
      exporter.exportFunction(uSoln, "u", 0, 4);
      exporter.exportFunction(uSoln, "u", 1, 5);
      exporter.exportFunction(uhatSoln, "uhat", 0, 4);
      exporter.exportFunction(uhatSoln, "uhat", 1, 5);
      // exporter.exportFunction(fhatSoln, "fhat", 0, 4);
      // exporter.exportFunction(fhatSoln, "fhat", 1, 5);
    }
    {
      XDMFExporter exporter(meshTopology, "PoissonSolution", false);
      exporter.exportSolution(solution, mesh, varFactory, 0, 2, cellIDToSubdivision(mesh, 10));
      refinementStrategy.refine(true);
      solution->solve(false);
      exporter.exportSolution(solution, mesh, varFactory, 1, 2, cellIDToSubdivision(mesh, 10));
    }
    // exporter.exportFunction(sigmaSoln, "Poisson-s", "sigma", 0, 5);
    // exporter.exportFunction(uhatSoln, "Poisson-uhat", "uhat", 1, 6);
  }

  {

//.........这里部分代码省略.........
  if (useMumpsIfAvailable) solver = Teuchos::rcp( new MumpsSolver );
#endif

//  double errorPercentage = 0.5; // for mesh refinements: ask to refine elements that account for 80% of the error in each step
//  Teuchos::RCP<RefinementStrategy> refinementStrategy;
//  refinementStrategy = Teuchos::rcp( new ErrorPercentageRefinementStrategy( soln, errorPercentage ));

  if (maxRefinements != 0)
  {
    cout << "Warning: maxRefinements is not 0, but the slice exporter implicitly assumes there won't be any refinements.\n";
  }

  MeshPtr mesh;

  MeshPtr prevMesh;
  SolutionPtr prevSoln;

  mesh = MeshFactory::rectilinearMesh(bf, dimensions, elementCounts, H1Order, delta_k, origin);

  if (rank==0) cout << "Initial mesh has " << mesh->getTopology()->activeCellCount() << " active (leaf) cells " << "and " << mesh->globalDofCount() << " degrees of freedom.\n";

  FunctionPtr sideParity = Function::sideParity();

  int lastFrameOutputted = -1;

  SolutionPtr soln;

  IPPtr ip;
  ip = bf->graphNorm();

  FunctionPtr u0 = Teuchos::rcp( new Cone_U0(0.0, 0.25, 0.1, 1.0, false) );

  BCPtr bc = BC::bc();
  bc->addDirichlet(qHat, inflowFilter, Function::zero()); // zero BCs enforced at the inflow boundary.
  bc->addDirichlet(qHat, SpatialFilter::matchingZ(t0), u0);

  MeshPtr initialMesh = mesh;

  int startingSlabNumber;
  if (previousSolutionTimeSlabNumber != -1)
  {
    startingSlabNumber = previousSolutionTimeSlabNumber + 1;

    if (rank==0) cout << "Loading mesh from " << previousMeshFile << endl;

    prevMesh = MeshFactory::loadFromHDF5(bf, previousMeshFile);
    prevSoln = Solution::solution(mesh, bc, RHS::rhs(), ip); // include BC and IP objects for sake of condensed dof interpreter setup...
    prevSoln->setUseCondensedSolve(useCondensedSolve);

    if (rank==0) cout << "Loading solution from " << previousSolutionFile << endl;
    prevSoln->loadFromHDF5(previousSolutionFile);

    double tn = (previousSolutionTimeSlabNumber+1) * timeLengthPerSlab;
    origin[2] = tn;
    mesh = MeshFactory::rectilinearMesh(bf, dimensions, elementCounts, H1Order, delta_k, origin);

    FunctionPtr q_prev = Function::solution(qHat, prevSoln);
    FunctionPtr q_transfer = Teuchos::rcp( new MeshTransferFunction(-q_prev, prevMesh, mesh, tn) ); // negate because the normals go in opposite directions

    bc = BC::bc();
    bc->addDirichlet(qHat, inflowFilter, Function::zero()); // zero BCs enforced at the inflow boundary.
    bc->addDirichlet(qHat, SpatialFilter::matchingZ(tn), q_transfer);

    double t_slab_final = (previousSolutionTimeSlabNumber+1) * timeLengthPerSlab;
    int frameOrdinal = 0;

void TransientTests::SetUp()
{
  ////////////////////   DECLARE VARIABLES   ///////////////////////
  // define test variables
  VarPtr v = varFactory.testVar("v", HGRAD);
  
  // define trial variables
  beta_n_u_hat = varFactory.fluxVar("\\widehat{\\beta \\cdot n }");
  u = varFactory.fieldVar("u");

  vector<double> beta;
  beta.push_back(1.0);
  beta.push_back(0.0);
  
  ////////////////////   BUILD MESH   ///////////////////////
  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 = 4, verticalCells = 4;
  
  // create a pointer to a new mesh:
  mesh = MeshFactory::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);
  prevTimeFlow = Teuchos::rcp(new Solution(mesh, nullBC, nullRHS, nullIP) );  
  flowResidual = Teuchos::rcp(new Solution(mesh, nullBC, nullRHS, nullIP) );  

  FunctionPtr u_prev_time = Teuchos::rcp( new PreviousSolutionFunction(prevTimeFlow, u) );
  
  ////////////////////   DEFINE BILINEAR FORM   ///////////////////////
  RHSPtr rhs = RHS::rhs();
  FunctionPtr invDt = Teuchos::rcp(new ScalarParamFunction(1.0/dt));    
  
  // v terms:
  bf->addTerm( beta * u, - v->grad() );
  bf->addTerm( beta_n_u_hat, v);

  // transient terms
  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();

  ////////////////////   CREATE BCs   ///////////////////////
  BCPtr bc = BC::bc();

  SpatialFilterPtr lBoundary = Teuchos::rcp( new LeftBoundary );
  FunctionPtr u1 = Teuchos::rcp( new InletBC );
  bc->addDirichlet(beta_n_u_hat, lBoundary, -u1);

  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 ConstantScalarFunction(u_free) );

  // prevTimeFlow->projectOntoMesh(functionMap);
  
}

int main(int argc, char *argv[])
{
  Teuchos::GlobalMPISession mpiSession(&argc, &argv, 0);
  
  int spaceDim = 2;
  int meshWidth = 2;
  bool conformingTraces = true;
  int H1Order = 2, delta_k = 3;
  double domainWidth = 1.0e-3;
  bool diagScaling = false;
  double h = domainWidth / meshWidth;
  double weight = h / 4.0; // ratio of area of square with sidelength h to its perimeter
  
  double sigma_weight = 1.0; // h / 4.0; // sigma = sigma_weight * u->grad()

  Space uHatSpace = conformingTraces ? HGRAD : L2;
  
  VarFactoryPtr vf = VarFactory::varFactory();
  
  // fields
  VarPtr u = vf->fieldVar("u");
  VarPtr sigma = vf->fieldVar("sigma", VECTOR_L2);
  
  // traces
  VarPtr u_hat = vf->traceVar("u_hat", uHatSpace);
  VarPtr sigma_n = vf->fluxVar("sigma_n");
  
  // tests
  VarPtr v = vf->testVar("v", HGRAD);
  VarPtr tau = vf->testVar("tau", HDIV);
  
  BFPtr bf = BF::bf(vf);
  
// standard BF:
//  bf->addTerm(sigma, v->grad());
//  bf->addTerm(sigma_n, v);
//  
//  bf->addTerm(sigma, tau);
//  bf->addTerm(u, tau->div());
//  bf->addTerm(-u_hat, tau->dot_normal());
  
  // weighted BF:
  bf->addTerm(sigma, v->grad());
  bf->addTerm(weight * sigma_n, v);
  
  bf->addTerm(sigma, tau);
  bf->addTerm(sigma_weight * u, tau->div());
  bf->addTerm(- sigma_weight * weight * u_hat, tau->dot_normal());
  
  IPPtr ip = IP::ip();
// standard IP:
  ip->addTerm(tau + v->grad());
  ip->addTerm(tau->div());
  ip->addTerm(v);
  ip->addTerm(tau);
  
  // weighted IP:
//  ip->addTerm(tau + v->grad());
//  ip->addTerm(sigma_weight * tau->div());
//  ip->addTerm(max(sigma_weight,1e-3) * v);
//  ip->addTerm(sigma_weight * weight * tau);
  
  BCPtr bc = BC::bc();
  bc->addDirichlet(u_hat, SpatialFilter::allSpace(), Function::zero());
  
  RHSPtr rhs = RHS::rhs();
  rhs->addTerm(1.0 * sigma_weight * v);
  
  vector<double> dimensions(spaceDim,domainWidth);
  vector<int> elementCounts(spaceDim,meshWidth);
  
  MeshPtr mesh = MeshFactory::rectilinearMesh(bf, dimensions, elementCounts, H1Order, delta_k);
  
  SolutionPtr soln = Solution::solution(mesh, bc, rhs, ip);
  
  soln->setUseCondensedSolve(true);
  soln->initializeLHSVector();
  soln->initializeStiffnessAndLoad();
  soln->populateStiffnessAndLoad();
  
  Teuchos::RCP<Epetra_RowMatrix> stiffness = soln->getStiffnessMatrix();
  
  double condNumber = conditionNumberLAPACK(*stiffness, diagScaling);
  
  cout << "condest (1-norm): " << condNumber << endl;
  
  return 0;
}

// tests to make sure the energy error calculated thru direct integration works for vector valued test functions too
bool ScratchPadTests::testLTResidual()
{
  double tol = 1e-11;
  int rank = Teuchos::GlobalMPISession::getRank();

  bool success = true;

  int nCells = 2;
  double eps = .1;

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

  // define test variables
  VarFactoryPtr varFactory = 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);

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

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

int main(int argc, char *argv[])
{
  Teuchos::GlobalMPISession mpiSession(&argc, &argv, NULL); // initialize MPI
  
  Teuchos::CommandLineProcessor cmdp(false,true); // false: don't throw exceptions; true: do return errors for unrecognized options
  
  int numElements = 3;
  vector<vector<double>> domainDim(3,vector<double>{0.0,1.0}); // first index: spaceDim; second: 0/1 for x0, x1, etc.
  int polyOrder = 2, delta_k = 1;
  int spaceDim = 2;
  
  cmdp.setOption("numElements", &numElements );
  cmdp.setOption("polyOrder", &polyOrder );
  cmdp.setOption("delta_k", &delta_k );
  cmdp.setOption("x0", &domainDim[0][0] );
  cmdp.setOption("x1", &domainDim[0][1] );
  cmdp.setOption("y0", &domainDim[1][0] );
  cmdp.setOption("y1", &domainDim[1][1] );
  cmdp.setOption("z0", &domainDim[2][0] );
  cmdp.setOption("z1", &domainDim[2][1] );
  cmdp.setOption("spaceDim", &spaceDim);
  
  if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL)
  {
#ifdef HAVE_MPI
    MPI_Finalize();
#endif
    return -1;
  }
  
  vector<double> x0(spaceDim);
  vector<double> domainSize(spaceDim);
  vector<int> elementCounts(spaceDim);
  for (int d=0; d<spaceDim; d++)
  {
    x0[d] = domainDim[d][0];
    domainSize[d] = domainDim[d][1] - x0[d];
    elementCounts[d] = numElements;
  }
  
  bool conformingTraces = true; // no difference for primal/continuous formulations
  PoissonFormulation formCG(spaceDim, conformingTraces, PoissonFormulation::CONTINUOUS_GALERKIN);
  VarPtr q = formCG.q();
  VarPtr phi = formCG.phi();
  BFPtr bf = formCG.bf();
  
  MeshPtr bubnovMesh = MeshFactory::rectilinearMesh(bf, domainSize, elementCounts, polyOrder, 0, x0);

  // Right now, hanging nodes don't work with continuous field variables
  // there is a GDAMinimumRule test demonstrating the failure, SolvePoisson2DContinuousGalerkinHangingNode.
  // make a mesh with hanging nodes (when spaceDim > 1)
//  {
//    set<GlobalIndexType> cellsToRefine = {0};
//    bubnovMesh->hRefine(cellsToRefine);
//  }

  RHSPtr rhs = RHS::rhs();
  rhs->addTerm(1.0 * q); // unit forcing
  
  IPPtr ip = Teuchos::null; // will give Bubnov-Galerkin
  BCPtr bc = BC::bc();
  bc->addDirichlet(phi, SpatialFilter::allSpace(), Function::zero());
  
  SolutionPtr solution = Solution::solution(bf, bubnovMesh, bc, rhs, ip);
  solution->solve();
  
  HDF5Exporter exporter(bubnovMesh, "PoissonContinuousGalerkin");
  exporter.exportSolution(solution);

  /**** Sort-of-primal experiment ****/
  // an experiment: try doing "primal" DPG with IBP to the boundary
//  ip = IP::ip();
//  ip->addTerm(q->grad());
//  ip->addTerm(q);
//
//  solution = Solution::solution(bf, bubnovMesh, bc, rhs, ip);
//  solution->solve();
//  
//  HDF5Exporter primalNoFluxExporter(bubnovMesh, "PoissonPrimalNoFlux");
//  primalNoFluxExporter.exportSolution(solution);
  
  //*** Primal Formulation ***//
  PoissonFormulation form(spaceDim, conformingTraces, PoissonFormulation::PRIMAL);
  q = form.q();
  phi = form.phi();
  bf = form.bf();
  
  bc = BC::bc();
  bc->addDirichlet(phi, SpatialFilter::allSpace(), Function::zero());
  
  rhs = RHS::rhs();
  rhs->addTerm(1.0 * q); // unit forcing
  
  MeshPtr primalMesh = MeshFactory::rectilinearMesh(bf, domainSize, elementCounts, polyOrder, delta_k, x0);
  
  ip = IP::ip();
  ip->addTerm(q->grad());
  ip->addTerm(q);

  // Right now, hanging nodes don't work with continuous field variables
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    ip->addTerm( tau22 );
  }
  else if (norm == 2)
  {
    // ip->addTerm( 0.5/sqrt(nu)*tau11-0.5/nu*tau22 );
    // ip->addTerm( 1./sqrt(nu)*tau12 );
    // ip->addTerm( 1./sqrt(nu)*tau12 );
    // ip->addTerm( 0.5/sqrt(nu)*tau22-0.5/nu*tau11 );
    ip->addTerm( tau11 );
    ip->addTerm( tau12 );
    ip->addTerm( tau12 );
    ip->addTerm( tau22 );

    ip->addTerm( 2*tau11->dx() + 2*tau12->dy() - 2*u1_prev*v1->dx() - u2_prev*v1->dy() - u2_prev*v2->dx() );
    ip->addTerm( 2*tau12->dx() + 2*tau22->dy() - 2*u2_prev*v2->dy() - u1_prev*v1->dy() - u1_prev*v2->dx() );

    ip->addTerm( 2*u1_prev*v1->dx() + u2_prev*v1->dy() + u2_prev*v2->dx() );
    ip->addTerm( 2*u2_prev*v2->dy() + u1_prev*v1->dy() + u1_prev*v2->dx() );

    ip->addTerm( sqrt(nu)*v1->grad() );
    ip->addTerm( sqrt(nu)*v2->grad() );

    ip->addTerm( v1 );
    ip->addTerm( v2 );
  }

  ////////////////////   CREATE BCs   ///////////////////////
  BCPtr bc = BC::bc();
  // 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);
  // bc->addDirichlet(u1hat, left, zero);
  // bc->addDirichlet(u2hat, left, zero);
  // bc->addDirichlet(u1hat, right, zero);
  // bc->addDirichlet(u2hat, right, zero);
  // bc->addDirichlet(u1hat, top, zero);
  // bc->addDirichlet(u2hat, top, zero);
  // bc->addDirichlet(u1hat, bottom, zero);
  // bc->addDirichlet(u2hat, bottom, zero);

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

//.........这里部分代码省略.........
    c = Function::vectorize(y-0.5, 0.5-x);
  }
//  FunctionPtr c = Function::vectorize(y, x);
  FunctionPtr n = Function::normal();

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

  bf->addTerm(u / dt, v);
  bf->addTerm(- theta * u, c * v->grad());
//  bf->addTerm(theta * u_hat, (c * n) * v);
  bf->addTerm(qHat, v);

  double width = 2.0, height = 2.0;
  int horizontalCells = numCells, verticalCells = numCells;
  double x0 = -0.5;
  double y0 = -0.5;

  if (usePeriodicBCs)
  {
    x0 = 0.0;
    y0 = 0.0;
    width = 1.0;
    height = 1.0;
  }

  BCPtr bc = BC::bc();

  SpatialFilterPtr inflowFilter  = Teuchos::rcp( new InflowFilterForClockwisePlanarRotation (x0,x0+width,y0,y0+height,0.5,0.5));

  vector< PeriodicBCPtr > periodicBCs;
  if (! usePeriodicBCs)
  {
    //  bc->addDirichlet(u_hat, SpatialFilter::allSpace(), Function::zero());
    bc->addDirichlet(qHat, inflowFilter, Function::zero()); // zero BCs enforced at the inflow boundary.
  }
  else
  {
    periodicBCs.push_back(PeriodicBC::xIdentification(x0, x0+width));
    periodicBCs.push_back(PeriodicBC::yIdentification(y0, y0+height));
  }

  MeshPtr mesh = MeshFactory::quadMeshMinRule(bf, H1Order, delta_k, width, height,
                 horizontalCells, verticalCells, false, x0, y0, periodicBCs);

  FunctionPtr u0 = Teuchos::rcp( new Cone_U0(0.0, 0.25, 0.1, 1.0, usePeriodicBCs) );

  RHSPtr initialRHS = RHS::rhs();
  initialRHS->addTerm(u0 / dt * v);
  initialRHS->addTerm((1-theta) * u0 * c * v->grad());

  IPPtr ip;
//  ip = Teuchos::rcp( new IP );
//  ip->addTerm(v);
//  ip->addTerm(c * v->grad());
  ip = bf->graphNorm();

  // create two Solution objects; we'll switch between these for time steps
  SolutionPtr soln0 = Solution::solution(mesh, bc, initialRHS, ip);
  soln0->setCubatureEnrichmentDegree(5);
  FunctionPtr u_soln0 = Function::solution(u, soln0);
  FunctionPtr qHat_soln0 = Function::solution(qHat, soln0);

  RHSPtr rhs1 = RHS::rhs();
  rhs1->addTerm(u_soln0 / dt * v);
  rhs1->addTerm((1-theta) * u_soln0 * c * v->grad());

//.........这里部分代码省略.........
  VarPtr tau2 = vf->testVar("tau2", HDIV);
  VarPtr q = vf->testVar("q", HGRAD);

  BFPtr bf = Teuchos::rcp( new BF(vf) );

  bf->addTerm(1./mu*sigma1, tau1);
  bf->addTerm(1./mu*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(sigma1, v1->grad());
  bf->addTerm(sigma2, v2->grad());
  bf->addTerm(-p, v1->dx());
  bf->addTerm(-p, v2->dy());
  bf->addTerm(t1c, v1);
  bf->addTerm(t2c, v2);
  bf->addTerm(mu*permInv*u1, v1);
  bf->addTerm(mu*permInv*u2, v2);

  bf->addTerm(-u1, q->dx());
  bf->addTerm(-u2, q->dy());
  bf->addTerm(u1hat, q->times_normal_x());
  bf->addTerm(u2hat, q->times_normal_y());

  RHSPtr rhs = RHS::rhs();

  BCPtr bc = BC::bc();

  SpatialFilterPtr y_equals_one = SpatialFilter::matchingY(1.0);
  SpatialFilterPtr y_equals_zero = SpatialFilter::matchingY(0);
  SpatialFilterPtr x_equals_one = SpatialFilter::matchingX(1.0);
  SpatialFilterPtr x_equals_zero = SpatialFilter::matchingX(0.0);
  bc->addDirichlet(u1hat, y_equals_zero, u1_exact);
  bc->addDirichlet(u2hat, y_equals_zero, u2_exact);
  bc->addDirichlet(u1hat, x_equals_zero, u1_exact);
  bc->addDirichlet(u2hat, x_equals_zero, u2_exact);
  bc->addDirichlet(u1hat, y_equals_one, u1_exact);
  bc->addDirichlet(u2hat, y_equals_one, u2_exact);
  bc->addDirichlet(u1hat, x_equals_one, u1_exact);
  bc->addDirichlet(u2hat, x_equals_one, u2_exact);
  bc->addZeroMeanConstraint(p);

  MeshPtr mesh = MeshFactory::quadMesh(bf, k+1, delta_k, 1, 1, 4, 4);

  map<string, IPPtr> brinkmanIPs;
  brinkmanIPs["Graph"] = bf->graphNorm();

  brinkmanIPs["Decoupled"] = Teuchos::rcp(new IP);
  brinkmanIPs["Decoupled"]->addTerm(tau1);
  brinkmanIPs["Decoupled"]->addTerm(tau2);
  brinkmanIPs["Decoupled"]->addTerm(tau1->div());
  brinkmanIPs["Decoupled"]->addTerm(tau2->div());
  brinkmanIPs["Decoupled"]->addTerm(permInv*v1);
  brinkmanIPs["Decoupled"]->addTerm(permInv*v2);
  brinkmanIPs["Decoupled"]->addTerm(v1->grad());
  brinkmanIPs["Decoupled"]->addTerm(v2->grad());
  brinkmanIPs["Decoupled"]->addTerm(q);
  brinkmanIPs["Decoupled"]->addTerm(q->grad());

  // brinkmanIPs["CoupledRobust"] = Teuchos::rcp(new IP);
  // brinkmanIPs["CoupledRobust"]->addTerm(tau->div()-beta*v->grad());
  // brinkmanIPs["CoupledRobust"]->addTerm(Function<double>::min(one/Function<double>::h(),Function<double>::constant(1./sqrt(epsilon)))*tau);
  // brinkmanIPs["CoupledRobust"]->addTerm(sqrt(epsilon)*v->grad());
  // brinkmanIPs["CoupledRobust"]->addTerm(beta*v->grad());
  // brinkmanIPs["CoupledRobust"]->addTerm(Function<double>::min(sqrt(epsilon)*one/Function<double>::h(),one)*v);

  IPPtr ip = brinkmanIPs[norm];

  SolutionPtr soln = TSolution<double>::solution(mesh, bc, rhs, ip);

  double threshold = 0.20;
  RefinementStrategy refStrategy(soln, threshold);

  ostringstream refName;
  refName << "brinkman";
  HDF5Exporter exporter(mesh,refName.str());

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

    double energyError = soln->energyErrorTotal();
    if (commRank == 0)
    {
      // if (refIndex > 0)
      // refStrategy.printRefinementStatistics(refIndex-1);
      cout << "Refinement:\t " << refIndex << " \tElements:\t " << mesh->numActiveElements()
           << " \tDOFs:\t " << mesh->numGlobalDofs() << " \tEnergy Error:\t " << energyError << endl;
    }

    exporter.exportSolution(soln, refIndex);

    if (refIndex != numRefs)
      refStrategy.refine();
  }

  return 0;
}

bool ScratchPadTests::testResidualMemoryError()
{

  int rank = Teuchos::GlobalMPISession::getRank();

  double tol = 1e-11;
  bool success = true;

  int nCells = 2;
  double eps = 1e-2;

  ////////////////////   DECLARE VARIABLES   ///////////////////////
  // define test variables
  VarFactoryPtr varFactory = 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);

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

  // robust test norm
  IPPtr robIP = Teuchos::rcp(new IP);
  robIP->addTerm(tau);
  robIP->addTerm(tau->div());
  robIP->addTerm(v->grad());
  robIP->addTerm(v);

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

  FunctionPtr zero = Function::constant(0.0);
  FunctionPtr one = Function::constant(1.0);
  RHSPtr rhs = RHS::rhs();
  FunctionPtr f = zero;
  //  FunctionPtr f = one;
  rhs->addTerm( f * v ); // obviously, with f = 0 adding this term is not necessary!

  ////////////////////   CREATE BCs   ///////////////////////
  BCPtr bc = BC::bc();
  SpatialFilterPtr inflowBoundary = Teuchos::rcp( new LRInflowSquareBoundary );
  SpatialFilterPtr outflowBoundary = Teuchos::rcp( new LROutflowSquareBoundary);

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

  bc->addDirichlet(beta_n_u_minus_sigma_n, inflowBoundary, beta*n*one);
  bc->addDirichlet(uhat, outflowBoundary, zero);

  ////////////////////   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);
  //  mesh->setPartitionPolicy(Teuchos::rcp(new ZoltanMeshPartitionPolicy("HSFC")));

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

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

  LinearTermPtr residual = rhs->linearTermCopy();
  residual->addTerm(-confusionBF->testFunctional(solution));
  RieszRepPtr rieszResidual = Teuchos::rcp(new RieszRep(mesh, robIP, residual));
  rieszResidual->computeRieszRep();
  FunctionPtr e_v = RieszRep::repFunction(v,rieszResidual);
//.........这里部分代码省略.........