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

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

#ifdef HAVE_MPI
  Teuchos::GlobalMPISession mpiSession(&argc, &argv,0);

  Epetra_MpiComm Comm(MPI_COMM_WORLD);
#else
  Epetra_SerialComm Comm;
#endif

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

  Comm.Barrier(); // set breakpoint here to allow debugger attachment to other MPI processes than the one you automatically attached to.

  Teuchos::CommandLineProcessor cmdp(false,true); // false: don't throw exceptions; true: do return errors for unrecognized options

  // problem parameters:
  double mu = 0.1;
  double permCoef = 1e4;
  int numRefs = 0;
  int k = 2, delta_k = 2;
  string norm = "Graph";
  cmdp.setOption("polyOrder",&k,"polynomial order for field variable u");
  cmdp.setOption("delta_k", &delta_k, "test space polynomial order enrichment");
  cmdp.setOption("numRefs",&numRefs,"number of refinements");
  cmdp.setOption("norm", &norm, "norm");
  cmdp.setOption("mu", &mu, "mu");
  cmdp.setOption("permCoef", &permCoef, "Permeability coefficient");

  if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL)
  {
#ifdef HAVE_MPI
    MPI_Finalize();
#endif
    return -1;
  }

  FunctionPtr zero = TFunction<double>::zero();
  FunctionPtr one = TFunction<double>::constant(1);
  FunctionPtr sin2pix = Teuchos::rcp( new Sin_ax(2*pi) );
  FunctionPtr cos2pix = Teuchos::rcp( new Cos_ax(2*pi) );
  FunctionPtr sin2piy = Teuchos::rcp( new Sin_ay(2*pi) );
  FunctionPtr cos2piy = Teuchos::rcp( new Cos_ay(2*pi) );
  FunctionPtr u1_exact = sin2pix*cos2piy;
  FunctionPtr u2_exact = -cos2pix*sin2piy;
  FunctionPtr x2 = TFunction<double>::xn(2);
  FunctionPtr y2 = TFunction<double>::yn(2);
  FunctionPtr p_exact = x2*y2 - 1./9;
  FunctionPtr permInv = permCoef*(sin2pix + 1.1);

  VarFactoryPtr vf = VarFactory::varFactory();
  //fields:
  VarPtr sigma1 = vf->fieldVar("sigma1", VECTOR_L2);
  VarPtr sigma2 = vf->fieldVar("sigma2", VECTOR_L2);
  VarPtr u1 = vf->fieldVar("u1", L2);
  VarPtr u2 = vf->fieldVar("u2", L2);
  VarPtr p = vf->fieldVar("p", L2);

  // traces:
  VarPtr u1hat = vf->traceVar("u1hat");
  VarPtr u2hat = vf->traceVar("u2hat");
  VarPtr t1c = vf->fluxVar("t1c");
  VarPtr t2c = vf->fluxVar("t2c");

  // test:
  VarPtr v1 = vf->testVar("v1", HGRAD);
  VarPtr v2 = vf->testVar("v2", HGRAD);
  VarPtr tau1 = vf->testVar("tau1", HDIV);
  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);
//.........这里部分代码省略.........

PoissonFormulation::PoissonFormulation(int spaceDim, bool useConformingTraces, PoissonFormulationChoice formulationChoice)
{
  _spaceDim = spaceDim;

  if (formulationChoice == ULTRAWEAK)
  {
    Space tauSpace = (spaceDim > 1) ? HDIV : HGRAD;
    Space phi_hat_space = useConformingTraces ? HGRAD : L2;
    Space psiSpace = (spaceDim > 1) ? VECTOR_L2 : L2;

    // fields
    VarPtr phi;
    VarPtr psi;

    // traces
    VarPtr phi_hat, psi_n_hat;

    // tests
    VarPtr q;
    VarPtr tau;

    VarFactoryPtr vf = VarFactory::varFactory();
    phi = vf->fieldVar(S_PHI);
    psi = vf->fieldVar(S_PSI, psiSpace);

    TFunctionPtr<double> parity = TFunction<double>::sideParity();
    
    if (spaceDim > 1)
      phi_hat = vf->traceVar(S_PHI_HAT, phi, phi_hat_space);
    else
      phi_hat = vf->fluxVar(S_PHI_HAT, phi * (Function::normal_1D() * parity), phi_hat_space); // for spaceDim==1, the "normal" component is in the flux-ness of phi_hat (it's a plus or minus 1)

    TFunctionPtr<double> n = TFunction<double>::normal();

    if (spaceDim > 1)
      psi_n_hat = vf->fluxVar(S_PSI_N_HAT, psi * (n * parity));
    else
      psi_n_hat = vf->fluxVar(S_PSI_N_HAT, psi * (Function::normal_1D() * parity));

    q = vf->testVar(S_Q, HGRAD);
    tau = vf->testVar(S_TAU, tauSpace);

    _poissonBF = Teuchos::rcp( new BF(vf) );

    if (spaceDim==1)
    {
      // for spaceDim==1, the "normal" component is in the flux-ness of phi_hat (it's a plus or minus 1)
      _poissonBF->addTerm(phi, tau->dx());
      _poissonBF->addTerm(psi, tau);
      _poissonBF->addTerm(-phi_hat, tau);

      _poissonBF->addTerm(-psi, q->dx());
      _poissonBF->addTerm(psi_n_hat, q);
    }
    else
    {
      _poissonBF->addTerm(phi, tau->div());
      _poissonBF->addTerm(psi, tau);
      _poissonBF->addTerm(-phi_hat, tau->dot_normal());

      _poissonBF->addTerm(-psi, q->grad());
      _poissonBF->addTerm(psi_n_hat, q);
    }
  }
  else if ((formulationChoice == PRIMAL) || (formulationChoice == CONTINUOUS_GALERKIN))
  {
    // field
    VarPtr phi;
    
    // flux
    VarPtr psi_n_hat;
    
    // tests
    VarPtr q;
    
    VarFactoryPtr vf = VarFactory::varFactory();
    phi = vf->fieldVar(S_PHI, HGRAD);
    
    TFunctionPtr<double> parity = TFunction<double>::sideParity();
    TFunctionPtr<double> n = TFunction<double>::normal();
    
    if (formulationChoice == PRIMAL)
    {
      if (spaceDim > 1)
        psi_n_hat = vf->fluxVar(S_PSI_N_HAT, phi->grad() * (n * parity));
      else
        psi_n_hat = vf->fluxVar(S_PSI_N_HAT, phi->dx() * (Function::normal_1D() * parity));
    }
    q = vf->testVar(S_Q, HGRAD);
    
    _poissonBF = BF::bf(vf);
    _poissonBF->addTerm(-phi->grad(), q->grad());

    if (formulationChoice == CONTINUOUS_GALERKIN)
    {
      FunctionPtr boundaryIndicator = Function::meshBoundaryCharacteristic();
      _poissonBF->addTerm(phi->grad() * n, boundaryIndicator * q);
    }
    else // primal
    {
//.........这里部分代码省略.........

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

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

//  double Re = 10.0;  // Cockburn Kanschat Stokes
//  double Re = 40.0; // Evans Hughes Navier-Stokes
//  double Re = 1000.0;

  string formulationTypeStr = "vgp";

  FunctionPtr u1_exact, u2_exact, p_exact;

  int numCellsFineMesh = 20; // for computing a zero-mean pressure
  int H1OrderFineMesh = 5;

//  VGPNavierStokesProblem(double Re, Intrepid::FieldContainer<double> &quadPoints, int horizontalCells,
//                         int verticalCells, int H1Order, int pToAdd,
//                         TFunctionPtr<double> u1_0, TFunctionPtr<double> u2_0, TFunctionPtr<double> f1, TFunctionPtr<double> f2,
//                         bool enrichVelocity = false, bool enhanceFluxes = false)
  
  FunctionPtr zero = Function::zero();
  VGPNavierStokesProblem zeroProblem = VGPNavierStokesProblem(Re, quadPointsKovasznay,
                                       numCellsFineMesh, numCellsFineMesh,
                                       H1OrderFineMesh, pToAdd,
                                       zero, zero, zero, useCompliantNorm || useStokesCompliantNorm);

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


  VarPtr v1_vgp = varFactory->testVar(VGP_V1_S, HGRAD);
  VarPtr v2_vgp = varFactory->testVar(VGP_V2_S, HGRAD);

//  if (rank==0) {
//    cout << "bilinear form with zero background flow:\n";
//    zeroProblem.bf()->printTrialTestInteractions();
//  }

  VGPStokesFormulation stokesForm(1/Re);

  NavierStokesFormulation::setKovasznay(Re, zeroProblem.mesh(), u1_exact, u2_exact, p_exact);

//  if (rank==0) cout << "Stokes bilinearForm: " << stokesForm.bf()->displayString() << endl;
  
  map< string, string > convergenceDataForMATLAB; // key: field file name

  for (int polyOrder = minPolyOrder; polyOrder <= maxPolyOrder; polyOrder++)
  {
    int H1Order = polyOrder + 1;

    int numCells1D = pow(2.0,minLogElements);

    if (rank==0)

void FunctionTests::setup()
{
  ////////////////////   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_const;
  beta_const.push_back(2.0);
  beta_const.pu