C++ VarPtr类代码示例

void PoissonExactSolution::setUseSinglePointBCForPHI(bool useSinglePointBCForPhi, IndexType vertexIndexForZeroValue)
{
  FunctionPtr phi_exact = phi();

  VarFactoryPtr vf = _bf->varFactory();
  VarPtr psi_hat_n = vf->fluxVar(PoissonBilinearForm::S_PSI_HAT_N);
  VarPtr q = vf->testVar(PoissonBilinearForm::S_Q, HGRAD);

  VarPtr phi = vf->fieldVar(PoissonBilinearForm::S_PHI);

  SpatialFilterPtr wholeBoundary = SpatialFilter::allSpace();

  FunctionPtr n = Function::normal();
  FunctionPtr psi_n_exact = phi_exact->grad() * n;

  _bc = BC::bc();
  _bc->addDirichlet(psi_hat_n, wholeBoundary, psi_n_exact);
  if (!useSinglePointBCForPhi)
  {
    _bc->addZeroMeanConstraint(phi);

  }
  else
  {
    std::vector<double> point = getPointForBCImposition();
    double value = Function::evaluate(phi_exact, point[0], point[1]);
//    cout << "PoissonExactSolution: imposing phi = " << value << " at (" << point[0] << ", " << point[1] << ")\n";
    _bc->addSpatialPointBC(phi->ID(), value, point);
  }
}

map<int, FunctionPtr > PreviousSolutionFunction::functionMap( vector< VarPtr > varPtrs, SolutionPtr soln) {
  map<int, FunctionPtr > functionMap;
  for (vector< VarPtr >::iterator varIt = varPtrs.begin(); varIt != varPtrs.end(); varIt++) {
    VarPtr var = *varIt;
    functionMap[var->ID()] = Teuchos::rcp( new PreviousSolutionFunction(soln, var));
  }
  return functionMap;
}

BFPtr TestBilinearFormDx::bf()
{
  VarFactoryPtr vf = VarFactory::varFactory();
  VarPtr u = vf->fieldVar("u",HGRAD);
  VarPtr v = vf->testVar("v",HGRAD);
  BFPtr bf = BF::bf(vf);
  bf->addTerm(u->dx(), v->dx());

  return bf;
}

bool LinearTermTests::testLinearTermEvaluation()
{
  bool success = true;
  double eps = .1;

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

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

  // define a couple LinearTerms
  LinearTermPtr vVecLT = Teuchos::rcp(new LinearTerm);
  LinearTermPtr tauVecLT = Teuchos::rcp(new LinearTerm);
  vVecLT->addTerm(sqrt(eps)*v->grad());
  tauVecLT->addTerm((1/sqrt(eps))*tau);

  //////////////////// evaluate LinearTerms /////////////////

  map<int,FunctionPtr> errRepMap;
  errRepMap[v->ID()] = one;
  errRepMap[tau->ID()] = one*e1+one*e2; // vector valued fxn (1,1)
  FunctionPtr errTau = tauVecLT->evaluate(errRepMap,false);
  FunctionPtr errV = vVecLT->evaluate(errRepMap,false);
  try
  {
    bool xTauZero = errTau->x()->isZero();
    bool yTauZero = errTau->y()->isZero();
    bool xVZero = errV->dx()->isZero();
    bool yVZero = errV->dy()->isZero();

  }
  catch (...)
  {
    cout << "testLinearTermEvaluation: Caught exception.\n";
    success = false;
  }
  /*
  FunctionPtr xErr = (errTau->x())*(errTau->x()) + (errV->dx())*(errV->dx());
  FunctionPtr yErr = (errTau->y())*(errTau->y()) + (errV->dy())*(errV->dy());
  double xErrVal = xErr->integrate(mesh,15,true);
  */

  // if we don't crash, return success
  return success;

}

void writePatchValues(double xMin, double xMax, double yMin, double yMax,
                      SolutionPtr solution, VarPtr u1, string filename,
                      int numPoints=100)
{
  FieldContainer<double> points = pointGrid(xMin,xMax,yMin,yMax,numPoints);
  FieldContainer<double> values(numPoints*numPoints);
  solution->solutionValues(values, u1->ID(), points);
  ofstream fout(filename.c_str());
  fout << setprecision(15);

  fout << "X = zeros(" << numPoints << ",1);\n";
  //    fout << "Y = zeros(numPoints);\n";
  fout << "U = zeros(" << numPoints << "," << numPoints << ");\n";
  for (int i=0; i<numPoints; i++)
  {
    fout << "X(" << i+1 << ")=" << points(i,0) << ";\n";
  }
  for (int i=0; i<numPoints; i++)
  {
    fout << "Y(" << i+1 << ")=" << points(i,1) << ";\n";
  }

  for (int i=0; i<numPoints; i++)
  {
    for (int j=0; j<numPoints; j++)
    {
      int pointIndex = i*numPoints + j;
      fout << "U("<<i+1<<","<<j+1<<")=" << values(pointIndex) << ";" << endl;
    }
  }
  fout.close();
}

PreviousSolutionFunction::PreviousSolutionFunction(SolutionPtr soln, VarPtr var, bool multiplyFluxesByCellParity) : Function(var->rank()) {
  _soln = soln;
  _solnExpression = 1.0 * var;
  _overrideMeshCheck = false;
  if ((var->varType() == FLUX) && multiplyFluxesByCellParity) {
    FunctionPtr parity = Teuchos::rcp( new SideParityFunction );
    _solnExpression = parity * var;
  }
}

double integralOverMesh(LinearTermPtr testTerm, VarPtr testVar, FunctionPtr fxnToSubstitute) {
  map<int, FunctionPtr > varAsFunction;
  varAsFunction[testVar->ID()] = fxnToSubstitute;
  
  FunctionPtr substituteOnBoundary = testTerm->evaluate(varAsFunction, true);
  FunctionPtr substituteOnInterior = testTerm->evaluate(varAsFunction, false);
  double integral = substituteOnBoundary->integrate(mesh);
  integral += substituteOnInterior->integrate(mesh);
  return integral;
}

 void VarItemDelegate::drawBar (VarPtr dt, QPainter * painter, const QStyleOptionViewItem & option, const QModelIndex & index ) const {
   (void)index;
   //  QItemDelegate::drawBackground(painter,option,index);
 
   if (dt!=0) {
     if (dt->getType() != VARTYPE_ID_BOOL && dt->hasValue() && dt->hasMinValue() && dt->hasMaxValue()) {
       painter->save();  
   
       QRectF rect=option.rect;
       double border_x=2.0;
       double width=rect.width()-(border_x*2.0);
       if (width < 0.0) width=0.0;
       rect.setX(rect.x()+border_x);
       rect.setWidth(width);
 
       double range=dt->getMaxValue() - dt->getMinValue();
       double value=dt->getValue()-dt->getMinValue();
       double pos=0.0;
       if (range!=0.0) pos=(value/range)*width;
       if (pos < 0.0) pos=0.0;
       if (pos > width) pos=width;
 
       QRectF bar=rect;
       double height=2;
       bar.setY(bar.y()+(bar.height() - (height + 1)));
       bar.setHeight(height);
       QRectF frame=bar;
       bar.setWidth(pos);
 
       QColor high_c;
       if (option.state & QStyle::State_Selected) {
         high_c=option.palette.color(QPalette::HighlightedText);
       } else {
         high_c=option.palette.color(QPalette::Highlight);
       }
       high_c.setAlpha(96);
       //painter->setPen(QPen(high_c,1));
       //painter->setBrush(Qt::NoBrush);
       painter->setPen(Qt::NoPen);
       painter->setBrush(QBrush(high_c));
       painter->drawRect(frame);
 
       high_c.setAlpha(192);
       QColor nocolor(255,255,255,0) ;
       painter->setPen(Qt::NoPen);
       painter->setBrush(QBrush(high_c));
       painter->drawRect(bar);
 
       painter->restore();
     }
   }
 }

PreviousSolutionFunction<Scalar>::PreviousSolutionFunction(TSolutionPtr<Scalar> soln, VarPtr var, bool multiplyFluxesByCellParity) : TFunction<Scalar>(var->rank())
{
  _soln = soln;
  _solnExpression = 1.0 * var;
  _overrideMeshCheck = false;
  if ((var->varType() == FLUX) && multiplyFluxesByCellParity)
  {
    TFunctionPtr<double> parity = TFunction<double>::sideParity();
    _solnExpression = parity * var;
  }
}

FieldContainer<double> solutionData(FieldContainer<double> &points, SolutionPtr solution, VarPtr u1) {
  int numPoints = points.dimension(0);
  FieldContainer<double> values(numPoints);
  solution->solutionValues(values, u1->ID(), points);
  
  FieldContainer<double> xyzData(numPoints, 3);
  for (int ptIndex=0; ptIndex<numPoints; ptIndex++) {
    xyzData(ptIndex,0) = points(ptIndex,0);
    xyzData(ptIndex,1) = points(ptIndex,1);
    xyzData(ptIndex,2) = values(ptIndex);
  }
  return xyzData;
}

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
  double epsilon = args.Input<double>("--epsilon", "diffusion parameter");
  int numRefs = args.Input<int>("--numRefs", "number of refinement steps");
  bool enforceLocalConservation = args.Input<bool>("--conserve", "enforce local conservation");
  int norm = args.Input<int>("--norm", "0 = graph\n    1 = robust\n    2 = modified robust");

  // Optional arguments (have defaults)
  bool zeroL2 = args.Input("--zeroL2", "take L2 term on v in robust norm to zero", false);
  args.Process();

  ////////////////////   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 sigma = varFactory.fieldVar("sigma", VECTOR_L2);

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

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

  // v terms:
  bf->addTerm( sigma, v->grad() );
  bf->addTerm( beta * u, - v->grad() );
  bf->addTerm( beta_n_u_minus_sigma_n, v);

  ////////////////////   DEFINE INNER PRODUCT(S)   ///////////////////////
  IPPtr ip = Teuchos::rcp(new IP);
  if (norm == 0)
  {
    ip = bf->graphNorm();
    FunctionPtr h2_scaling = Teuchos::rcp( new ZeroMeanScaling );
    ip->addZeroMeanTerm( h2_scaling*v );
  }
  // Robust norm
  else if (norm == 1)
  {
    // robust test norm
    FunctionPtr ip_scaling = Teuchos::rcp( new EpsilonScaling(epsilon) );
    FunctionPtr h2_scaling = Teuchos::rcp( new ZeroMeanScaling );
    if (!zeroL2)
      ip->addTerm( v );
    ip->addTerm( sqrt(epsilon) * v->grad() );
    // Weight these two terms for inflow
    ip->addTerm( beta * v->grad() );
    ip->addTerm( tau->div() );
    ip->addTerm( ip_scaling/sqrt(epsilon) * tau );
    if (zeroL2)
      ip->addZeroMeanTerm( h2_scaling*v );
  }
  // Modified robust norm
  else if (norm == 2)
  {
    // robust test norm
    FunctionPtr ip_scaling = Teuchos::rcp( new EpsilonScaling(epsilon) );
    FunctionPtr h2_scaling = Teuchos::rcp( new ZeroMeanScaling );
    // FunctionPtr ip_weight = Teuchos::rcp( new IPWeight() );
    if (!zeroL2)
      ip->addTerm( v );
    ip->addTerm( sqrt(epsilon) * v->grad() );
    ip->addTerm( beta * v->grad() );
    ip->addTerm( tau->div() - beta*v->grad() );
    ip->addTerm( ip_scaling/sqrt(epsilon) * tau );
    if (zeroL2)
      ip->addZeroMeanTerm( h2_scaling*v );
  }

  // // robust test norm
  // IPPtr robIP = Teuchos::rcp(new IP);
  // FunctionPtr ip_scaling = Teuchos::rcp( new EpsilonScaling(epsilon) );
  // if (!enforceLocalConservation)
  //   robIP->addTerm( ip_scaling * v );
  // robIP->addTerm( sqrt(epsilon) * v->grad() );
  // // Weight these two terms for inflow
  // FunctionPtr ip_weight = Teuchos::rcp( new IPWeight() );
  // robIP->addTerm( ip_weight * beta * v->grad() );
  // robIP->addTerm( ip_weight * tau->div() );
//.........这里部分代码省略.........

void MainWindow::notificationExample(VarPtr node) {
  printf("Value of node '%s' was changed to: '%s'\n",node->getName().c_str(),node->getString().c_str());
}

bool LinearTermTests::testIntegrateMixedBasis()
{
  bool success = true;

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

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

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

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

  // build CONSTANT SINGLE ELEMENT MESH
  int order = 0;
  int H1Order = order+1;
  int pToAdd = 1;
  int nCells = 2; // along a side

  // create a pointer to a new mesh:
  Teuchos::RCP<Mesh> mesh = MeshUtilities::buildUnitQuadMesh(nCells,convectionBF, H1Order, H1Order+pToAdd);
  ElementTypePtr elemType = mesh->getElement(0)->elementType();
  BasisCachePtr basisCache = Teuchos::rcp(new BasisCache(elemType, mesh));
  vector<GlobalIndexType> cellIDs;
  vector< ElementPtr > allElems = mesh->activeElements();
  vector< ElementPtr >::iterator elemIt;
  for (elemIt=allElems.begin(); elemIt!=allElems.end(); elemIt++)
  {
    cellIDs.push_back((*elemIt)->cellID());
  }
  bool createSideCacheToo = true;
  basisCache->setPhysicalCellNodes(mesh->physicalCellNodesGlobal(elemType), cellIDs, createSideCacheToo);

  int numTrialDofs = elemType->trialOrderPtr->totalDofs();
  int numCells = mesh->numActiveElements();
  double areaPerCell = 1.0 / numCells;
  FieldContainer<double> integrals(numCells,numTrialDofs);
  FieldContainer<double> expectedIntegrals(numCells,numTrialDofs);
  double sidelengthOfCell = 1.0 / nCells;
  DofOrderingPtr trialOrdering = elemType->trialOrderPtr;
  int dofForField = trialOrdering->getDofIndex(u->ID(), 0);
  vector<int> dofsForFlux;
  const vector<int>* sidesForFlux = &trialOrdering->getSidesForVarID(beta_n_u_hat->ID());
  for (vector<int>::const_iterator sideIt = sidesForFlux->begin(); sideIt != sidesForFlux->end(); sideIt++)
  {
    int sideIndex = *sideIt;
    dofsForFlux.push_back(trialOrdering->getDofIndex(beta_n_u_hat->ID(), 0, sideIndex));
  }
  for (int cellIndex = 0; cellIndex < numCells; cellIndex++)
  {
    expectedIntegrals(cellIndex, dofForField) = areaPerCell;
    for (vector<int>::iterator dofIt = dofsForFlux.begin(); dofIt != dofsForFlux.end(); dofIt++)
    {
      int fluxDofIndex = *dofIt;
      expectedIntegrals(cellIndex, fluxDofIndex) = sidelengthOfCell;
    }
  }

//  cout << "expectedIntegrals:\n" << expectedIntegrals;

  // setup: with constant degrees of freedom, we expect that the integral of int_dK (flux) + int_K (field) will be ones for each degree of freedom, assuming the basis functions for these constants field/flux variables are just C = 1.0.
  //
  //On a unit square, int_K (constant) = 1.0, and int_dK (u_i) = 1, for i = 0,...,3.

  LinearTermPtr lt = 1.0 * beta_n_u_hat;
  LinearTermPtr field =  1.0 * u;
  lt->addTerm(field,true);
  lt->integrate(integrals, elemType->trialOrderPtr, basisCache);

  double tol = 1e-12;
  double maxDiff;
  success = TestSuite::fcsAgree(integrals,expectedIntegrals,tol,maxDiff);
  if (success==false)
  {
    cout << "Failed testIntegrateMixedBasis with maxDiff = " << maxDiff << endl;
  }

  return success;
}

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)
  bool enforceLocalConservation = args.Input<bool>("--conserve", "enforce local conservation", false);
  double Re = args.Input("--Re", "Reynolds number", 40);
  double nu = 1./Re;
  double lambda = Re/2.-sqrt(Re*Re/4+4*pi*pi);
  int maxNewtonIterations = args.Input("--maxIterations", "maximum number of Newton iterations", 20);
  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();

  // if (commRank==0)
  // {
  //   cout << "saveFile is " << saveFile << endl;
  //   cout << "loadFile is " << replayFile << endl;
  // }

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

  ////////////////////   DECLARE VARIABLES   ///////////////////////
  // define test variables
  VarFactory varFactory;
  // VarPtr tau11 = varFactory.testVar("tau11", HGRAD);
  // VarPtr tau12 = varFactory.testVar("tau12", HGRAD);
  // VarPtr tau22 = varFactory.testVar("tau22", HGRAD);
  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 q = varFactory.testVar("q", HGRAD);

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

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

  // define nodes for mesh
  FieldContainer<double> meshBoundary(4,2);
  double xmin = -0.5;
  double xmax =  1.0;
  double ymin = -0.5;
  double ymax =  1.5;

  meshBoundary(0,0) =  xmin; // x1
  meshBoundary(0,1) =  ymin; // y1
  meshBoundary(1,0) =  xmax;
  meshBoundary(1,1) =  ymin;
  meshBoundary(2,0) =  xmax;
  meshBoundary(2,1) =  ymax;
  meshBoundary(3,0) =  xmin;
  meshBoundary(3,1) =  ymax;

  int horizontalCells = 6, verticalCells = 8;

  // create a pointer to a new mesh:
  Teuchos::RCP<Mesh> mesh = Mesh::buildQuadMesh(meshBoundary, horizontalCells, verticalCells,
                            bf, H1Order, 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;
//.........这里部分代码省略.........

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

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