C++ VarFactoryPtr类代码示例

void RHSTests::setup()
{
    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;

    int H1Order = 3;
    int delta_p = 3; // for test functions
    int horizontalCells = 2;
    int verticalCells = 2;

    double eps = 1.0; // not really testing for sharp gradients right now--just want to see if things basically work
    double beta_x = 1.0;
    double beta_y = 1.0;

    BFPtr confusionBF = ConfusionBilinearForm::confusionBF(eps,beta_x,beta_y);
    Teuchos::RCP<ConfusionProblemLegacy> confusionProblem = Teuchos::rcp( new ConfusionProblemLegacy(confusionBF, beta_x, beta_y) );
    _rhs = confusionProblem;
    _mesh = MeshFactory::buildQuadMesh(quadPoints, horizontalCells, verticalCells, confusionBF, H1Order, H1Order+delta_p);
    _mesh->setUsePatchBasis(false);

    VarFactoryPtr varFactory = confusionBF->varFactory(); // Create test IDs that match the enum in ConfusionBilinearForm
    _tau = varFactory->testVar(ConfusionBilinearForm::S_TAU,HDIV);
    _v = varFactory->testVar(ConfusionBilinearForm::S_V,HGRAD);

    _rhsEasy = RHS::rhs();
    _rhsEasy->addTerm( _v );
}

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

VarPtr PoissonFormulation::tau()
{
  VarFactoryPtr vf = _poissonBF->varFactory();
  if (_spaceDim > 1)
    return vf->testVar(S_TAU, HDIV);
  else
    return vf->testVar(S_TAU, HGRAD);
}

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;

}

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

  VarFactoryPtr varFactory = 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) ) );
  }

  VarFactoryPtr varFactory = 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;
}

VarPtr PressurelessStokesFormulation::v(int i)
{
  if (i > _spaceDim)
  {
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "i must be less than or equal to _spaceDim");
  }
  VarFactoryPtr vf = _stokesBF->varFactory();
  switch (i)
  {
  case 1:
    return vf->testVar(S_V1, HGRAD);
  case 2:
    return vf->testVar(S_V2, HGRAD);
  case 3:
    return vf->testVar(S_V3, HGRAD);
  }
  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "unhandled i value");
}

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

  VarFactoryPtr varFactory = 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;
}

ConfusionManufacturedSolution::ConfusionManufacturedSolution(double epsilon, double beta_x, double beta_y)
{
  _epsilon = epsilon;
  _beta_x  = beta_x;
  _beta_y  = beta_y;

  // set the class variables from ExactSolution:
//  _bc = Teuchos::rcp(this,false);  // false: don't let the RCP own the memory
//  _rhs = Teuchos::rcp(this,false);

  BFPtr bf = ConfusionBilinearForm::confusionBF(epsilon,beta_x,beta_y);

  _bilinearForm = bf;

  VarFactoryPtr vf = bf->varFactory();

  VarPtr u = vf->fieldVar(ConfusionBilinearForm::S_U);
  VarPtr sigma1 = vf->fieldVar(ConfusionBilinearForm::S_SIGMA_1);
  VarPtr sigma2 = vf->fieldVar(ConfusionBilinearForm::S_SIGMA_2);

  _u_hat = vf->traceVar(ConfusionBilinearForm::S_U_HAT);
  _beta_n_u_minus_sigma_hat = vf->fluxVar(ConfusionBilinearForm::S_BETA_N_U_MINUS_SIGMA_HAT);

  _v = vf->testVar(ConfusionBilinearForm::S_V, HGRAD);

  FunctionPtr u_exact = this->u();
  FunctionPtr sigma_exact = epsilon * u_exact->grad();

  FunctionPtr u_exact_laplacian = u_exact->dx()->dx() + u_exact->dy()->dy();

  _rhs = RHS::rhs();
  FunctionPtr f = - _epsilon * u_exact_laplacian + _beta_x * u_exact->dx() + _beta_y * u_exact->dy();
  _rhs->addTerm( f * _v );

  _bc = BC::bc();
  _bc->addDirichlet(_u_hat, SpatialFilter::allSpace(), u_exact);

  FunctionPtr beta = Function::vectorize(Function::constant(_beta_x), Function::constant(_beta_y));
  FunctionPtr n = Function::normal();
  FunctionPtr one_skeleton = Function::meshSkeletonCharacteristic(); // allows restriction to skeleton
  FunctionPtr sigma_flux_exact = beta * ( n * u_exact - sigma_exact * one_skeleton);

  this->setSolutionFunction(u, u_exact);
  this->setSolutionFunction(sigma1, sigma_exact->x());
  this->setSolutionFunction(sigma2, sigma_exact->y());
  this->setSolutionFunction(_u_hat, u_exact);
  this->setSolutionFunction(_beta_n_u_minus_sigma_hat, sigma_flux_exact);
}

void LinearTermTests::setup()
{

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

  sine_x = Teuchos::rcp( new Sine_x );
  cos_y = Teuchos::rcp( new Cosine_y );

  VarFactoryPtr varFactory = VarFactory::varFactory();
  q1 = varFactory->testVar("q_1", HDIV);
  q2 = varFactory->testVar("q_2", HDIV);
  q3 = varFactory->testVar("q_3", HDIV);

  v1 = varFactory->testVar("v_1", HGRAD);
  v2 = varFactory->testVar("v_2", HGRAD);
  v3 = varFactory->testVar("v_3", HGRAD);

  u1 = varFactory->fieldVar("u_1", HGRAD);
  u2 = varFactory->fieldVar("u_2", HGRAD);
  u3 = varFactory->fieldVar("u_3", HGRAD);

  u1_hat = varFactory->traceVar("\\widehat{u}_1");
  u2_hat = varFactory->traceVar("\\widehat{u}_2");

  u3_hat_n = varFactory->fluxVar("\\widehat{u}_3n");

  bf = Teuchos::rcp(new BF(varFactory)); // made-up bf for Mesh + previous solution tests

  bf->addTerm(u1_hat, q1->dot_normal());
  bf->addTerm(u1, q1->x());
  bf->addTerm(u2, q1->y());

  bf->addTerm(u3_hat_n, v1);
  bf->addTerm(u3, v1);

//  DofOrderingFactory discreteSpaceFactory(bf);

  int polyOrder = 3, testToAdd = 2;
  Teuchos::RCP<shards::CellTopology> quadTopoPtr;
//  quadTopoPtr = Teuchos::rcp(new shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >() ));

  // 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 horizontalElements = 2, verticalElements = 2;

  mesh = MeshFactory::buildQuadMesh(quadPoints, horizontalElements, verticalElements, bf, polyOrder+1, polyOrder+1+testToAdd);

  ElementTypePtr elemType = mesh->getElement(0)->elementType();
  trialOrder = elemType->trialOrderPtr;
  testOrder = elemType->testOrderPtr;

  basisCache = Teuchos::rcp(new BasisCache(elemType, mesh));

  vector<GlobalIndexType> cellIDs;
  cellIDs.push_back(0);
  cellIDs.push_back(1);
  cellIDs.push_back(2);
  cellIDs.push_back(3);
  bool createSideCacheToo = true;

  basisCache->setPhysicalCellNodes(mesh->physicalCellNodesGlobal(elemType), cellIDs, createSideCacheToo);
}

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[])
{
  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;
}

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

// tests Riesz inversion by integration by parts
bool LinearTermTests::testRieszInversion()
{
  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 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);
  double eps = .01;

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

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

  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;

  int nCells = 1;
  int horizontalCells = nCells, verticalCells = nCells;
  // create a pointer to a new mesh:
  Teuchos::RCP<Mesh> myMesh = MeshFactory::buildQuadMesh(quadPoints, horizontalCells, verticalCells,
                              confusionBF, H1Order, H1Order+pToAdd);

  ElementTypePtr elemType = myMesh->getElement(0)->elementType();
  BasisCachePtr basisCache = Teuchos::rcp(new BasisCache(elemType, myMesh));

  vector<GlobalIndexType> cellIDs;
  vector<ElementPtr> elems = myMesh->activeElements();
  vector<ElementPtr>::iterator elemIt;
  for (elemIt=elems.begin(); elemIt!=elems.end(); elemIt++)
  {
    int cellID = (*elemIt)->cellID();
    cellIDs.push_back(cellID);
  }
  bool createSideCacheToo = true;

  basisCache->setPhysicalCellNodes(myMesh->physicalCellNodesGlobal(elemType), cellIDs, createSideCacheToo);

  LinearTermPtr integrand = Teuchos::rcp(new LinearTerm);// residual
  LinearTermPtr integrandIBP = Teuchos::rcp(new LinearTerm);// residual

  vector<double> e1(2); // (1,0)
  vector<double> e2(2); // (0,1)
  e1[0] = 1;
  e2[1] = 1;
  FunctionPtr n = Function::normal();
  FunctionPtr X = Function::xn(1);
  FunctionPtr Y = Function::yn(1);
  FunctionPtr testFxn1 = X;
  FunctionPtr testFxn2 = Y;
  FunctionPtr divTestFxn = testFxn1->dx() + testFxn2->dy();
  FunctionPtr vectorTest = testFxn1*e1 + testFxn2*e2;

  integrand->addTerm(divTestFxn*v);
  integrandIBP->addTerm(vectorTest*n*v - vectorTest*v->grad()); // boundary term

  IPPtr sobolevIP = Teuchos::rcp(new IP);
  sobolevIP->addTerm(v);
  sobolevIP->addTerm(tau);

  Teuchos::RCP<RieszRep> riesz = Teuchos::rcp(new RieszRep(myMesh, sobolevIP, integrand));
  //  riesz->setPrintOption(true);
  riesz->computeRieszRep();
//.........这里部分代码省略.........