C++ BasisCachePtr::setPhysicalCellNodes方法代码示例

bool ScratchPadTests::testConstantFunctionProduct()
{
  bool success = true;
  // set up basisCache (even though it won't really be used here)
  BasisCachePtr basisCache = Teuchos::rcp( new BasisCache( _elemType, _spectralConfusionMesh ) );
  vector<GlobalIndexType> cellIDs;
  int cellID = 0;
  cellIDs.push_back(cellID);
  basisCache->setPhysicalCellNodes( _spectralConfusionMesh->physicalCellNodesForCell(cellID),
                                    cellIDs, true );

  int numCells = _basisCache->getPhysicalCubaturePoints().dimension(0);
  int numPoints = _testPoints.dimension(0);
  FunctionPtr three = Function::constant(3.0);
  FunctionPtr two = Function::constant(2.0);

  FieldContainer<double> values(numCells,numPoints);
  two->values(values,basisCache);
  three->scalarMultiplyBasisValues( values, basisCache );

  FieldContainer<double> expectedValues(numCells,numPoints);
  expectedValues.initialize( 3.0 * 2.0 );

  double tol = 1e-15;
  double maxDiff = 0.0;
  if ( ! fcsAgree(expectedValues, values, tol, maxDiff) )
  {
    success = false;
    cout << "Expected product differs from actual; maxDiff: " << maxDiff << endl;
  }
  return success;
}

bool FunctionTests::testAdaptiveIntegrate()
{
  bool success = true;

  // we must create our own basisCache here because _basisCache
  // has had its ref cell points set, which basically means it's
  // opted out of having any help with integration.
  BasisCachePtr basisCache = Teuchos::rcp( new BasisCache( _elemType, _spectralConfusionMesh ) );
  vector<GlobalIndexType> cellIDs;
  cellIDs.push_back(0);
  basisCache->setPhysicalCellNodes( _spectralConfusionMesh->physicalCellNodesForCell(0), cellIDs, true );

  double eps = .1; //
  FunctionPtr boundaryLayerFunction = Teuchos::rcp( new BoundaryLayerFunction(eps) );
  int numCells = basisCache->cellIDs().size();
  FieldContainer<double> integrals(numCells);
  double quadtol = 1e-2;
  double computedIntegral = boundaryLayerFunction->integrate(_spectralConfusionMesh,quadtol);
  double trueIntegral = (eps*(exp(1/eps) - exp(-1/eps)))*2.0;
  double diff = trueIntegral-computedIntegral;
  double relativeError = abs(diff)/abs(trueIntegral); // relative error

  double tol = 1e-2;
  if (relativeError > tol)
  {
    success = false;
    cout << "failing testAdaptiveIntegrate() with computed integral " << computedIntegral << " and true integral " << trueIntegral << endl;
  }
  return success;
}

void PreviousSolutionFunction::values(FieldContainer<double> &values, BasisCachePtr basisCache) {
  int rank = Teuchos::GlobalMPISession::getRank();
  if (_overrideMeshCheck) {
    _solnExpression->evaluate(values, _soln, basisCache);
    return;
  }
  if (!basisCache.get()) cout << "basisCache is nil!\n";
  if (!_soln.get()) cout << "_soln is nil!\n";
  // values are stored in (C,P,D) order
  if (basisCache->mesh().get() == _soln->mesh().get()) {
    _solnExpression->evaluate(values, _soln, basisCache);
  } else {
    static bool warningIssued = false;
    if (!warningIssued) {
      if (rank==0)
        cout << "NOTE: In PreviousSolutionFunction, basisCache's mesh doesn't match solution's.  If this is not what you intended, it would be a good idea to make sure that the mesh is passed in on BasisCache construction; the evaluation will be a lot slower without it...\n";
      warningIssued = true;
    }
    // get the physicalPoints, and make a basisCache for each...
    FieldContainer<double> physicalPoints = basisCache->getPhysicalCubaturePoints();
    FieldContainer<double> value(1,1); // assumes scalar-valued solution function.
    int numCells = physicalPoints.dimension(0);
    int numPoints = physicalPoints.dimension(1);
    int spaceDim = physicalPoints.dimension(2);
    physicalPoints.resize(numCells*numPoints,spaceDim);
    vector< ElementPtr > elements = _soln->mesh()->elementsForPoints(physicalPoints, false); // false: don't make elements null just because they're off-rank.
    FieldContainer<double> point(1,spaceDim);
    FieldContainer<double> refPoint(1,spaceDim);
    int combinedIndex = 0;
    vector<GlobalIndexType> cellID;
    cellID.push_back(-1);
    BasisCachePtr basisCacheOnePoint;
    for (int cellIndex=0; cellIndex<numCells; cellIndex++) {
      for (int ptIndex=0; ptIndex<numPoints; ptIndex++, combinedIndex++) {
        if (elements[combinedIndex].get()==NULL) continue; // no element found for point; skip it…
        ElementTypePtr elemType = elements[combinedIndex]->elementType();
        for (int d=0; d<spaceDim; d++) {
          point(0,d) = physicalPoints(combinedIndex,d);
        }
        if (elements[combinedIndex]->cellID() != cellID[0]) {
          cellID[0] = elements[combinedIndex]->cellID();
          basisCacheOnePoint = Teuchos::rcp( new BasisCache(elemType, _soln->mesh()) );
          basisCacheOnePoint->setPhysicalCellNodes(_soln->mesh()->physicalCellNodesForCell(cellID[0]),cellID,false); // false: don't createSideCacheToo
        }
        // compute the refPoint:
        typedef CellTools<double>  CellTools;
        int whichCell = 0;
        CellTools::mapToReferenceFrame(refPoint,point,_soln->mesh()->physicalCellNodesForCell(cellID[0]),
                                       *(elemType->cellTopoPtr),whichCell);
        basisCacheOnePoint->setRefCellPoints(refPoint);
        //          cout << "refCellPoints:\n " << refPoint;
        //          cout << "physicalCubaturePoints:\n " << basisCacheOnePoint->getPhysicalCubaturePoints();
        _solnExpression->evaluate(value, _soln, basisCacheOnePoint);
        //          cout << "value at point (" << point(0,0) << ", " << point(0,1) << ") = " << value(0,0) << endl;
        values(cellIndex,ptIndex) = value(0,0);
      }
    }
  }
}

bool RHSTests::testIntegrateAgainstStandardBasis()
{
    bool success = true;
    double tol = 1e-14;

    Teuchos::RCP<ElementType> elemType = _mesh->getElement(0)->elementType();

    int rank = _mesh->Comm()->MyPID();
    vector< Teuchos::RCP< Element > > elemsInPartitionOfType = _mesh->elementsOfType(rank, elemType);
    FieldContainer<double> physicalCellNodes = _mesh->physicalCellNodes(elemType);

    int numCells = elemsInPartitionOfType.size();
    int numTestDofs = elemType->testOrderPtr->totalDofs();

// set up diagonal testWeights matrices so we can reuse the existing computeRHS, and compare results…
    FieldContainer<double> testWeights(numCells,numTestDofs,numTestDofs);
    for (int cellIndex=0; cellIndex<numCells; cellIndex++)
    {
        for (int i=0; i<numTestDofs; i++)
        {
            testWeights(cellIndex,i,i) = 1.0;
        }
    }

    FieldContainer<double> rhsExpected(numCells,numTestDofs);
    FieldContainer<double> rhsActual(numCells,numTestDofs);

    // determine cellIDs
    vector<GlobalIndexType> cellIDs = _mesh->globalDofAssignment()->cellIDsOfElementType(rank, elemType);

    if (numCells > 0)
    {
        // prepare basisCache and cellIDs
        BasisCachePtr basisCache = Teuchos::rcp(new BasisCache(elemType,_mesh));
        bool createSideCacheToo = true;
        basisCache->setPhysicalCellNodes(physicalCellNodes,cellIDs,createSideCacheToo);
        _rhs->integrateAgainstStandardBasis(rhsActual, elemType->testOrderPtr, basisCache);
        _rhs->integrateAgainstOptimalTests(rhsExpected, testWeights, elemType->testOrderPtr, basisCache);
    }

    double maxDiff = 0.0;

    if ( ! fcsAgree(rhsExpected,rhsActual,tol,maxDiff) )
    {
        success = false;
        cout << "Failed testIntegrateAgainstStandardBasis: maxDiff = " << maxDiff << endl;
    }

    // check success across MPI nodes
    return allSuccess(success);
}

bool RHSTests::testTrivialRHS()
{

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

    int rank = _mesh->Comm()->MyPID();

    Teuchos::RCP<ElementType> elemType = _mesh->getElement(0)->elementType();

    vector< Teuchos::RCP< Element > > elemsInPartitionOfType = _mesh->elementsOfType(rank, elemType);
    FieldContainer<double> physicalCellNodes = _mesh->physicalCellNodes(elemType);

    int numCells = elemsInPartitionOfType.size();
    int numTestDofs = elemType->testOrderPtr->totalDofs();

    // determine cellIDs
    vector<GlobalIndexType> cellIDs = _mesh->globalDofAssignment()->cellIDsOfElementType(rank, elemType);

    if (numCells > 0)
    {
        // prepare basisCache and cellIDs
        BasisCachePtr basisCache = Teuchos::rcp(new BasisCache(elemType,_mesh));
        bool createSideCacheToo = true;
        basisCache->setPhysicalCellNodes(physicalCellNodes,cellIDs,createSideCacheToo);

        FieldContainer<double> rhsExpected(numCells,numTestDofs);

        FunctionPtr zero = Function::constant(0.0);
        RHSPtr rhs = RHS::rhs();
        FunctionPtr f = zero;
        rhs->addTerm( f * _v ); // obviously, with f = 0 adding this term is not necessary!
        rhs->integrateAgainstStandardBasis(rhsExpected, elemType->testOrderPtr, basisCache);

        for (int i = 0; i<numCells; i++)
        {
            for (int j = 0; j<numTestDofs; j++)
            {
                if (abs(rhsExpected(i,j))>tol)
                {
                    success = false;
                }
            }
        }
    }

    return allSuccess(success);
}

bool RHSTests::testRHSEasy()
{
    bool success = true;
    double tol = 1e-14;

    int rank = _mesh->Comm()->MyPID();

    int numProcs = Teuchos::GlobalMPISession::getNProc();

    Teuchos::RCP<ElementType> elemType = _mesh->getElement(0)->elementType();

    vector< Teuchos::RCP< Element > > elemsInPartitionOfType = _mesh->elementsOfType(rank, elemType);
    FieldContainer<double> physicalCellNodes = _mesh->physicalCellNodes(elemType);

    int numCells = elemsInPartitionOfType.size();
    int numTestDofs = elemType->testOrderPtr->totalDofs();

    FieldContainer<double> rhsExpected(numCells,numTestDofs);
    FieldContainer<double> rhsActual(numCells,numTestDofs);

    // determine cellIDs
    vector<GlobalIndexType> cellIDs = _mesh->globalDofAssignment()->cellIDsOfElementType(rank, elemType);

    // prepare basisCache and cellIDs
    if (numCells > 0)
    {
        BasisCachePtr basisCache = Teuchos::rcp(new BasisCache(elemType,_mesh));
        bool createSideCacheToo = true;
        basisCache->setPhysicalCellNodes(physicalCellNodes,cellIDs,createSideCacheToo);

        _rhs->integrateAgainstStandardBasis(rhsActual, elemType->testOrderPtr, basisCache);
        _rhsEasy->integrateAgainstStandardBasis(rhsExpected, elemType->testOrderPtr, basisCache);
    }

    double maxDiff = 0.0;

    if ( ! fcsAgree(rhsExpected,rhsActual,tol,maxDiff) )
    {
        success = false;
        cout << "Failed testRHSEasy: maxDiff = " << maxDiff << endl;
        cout << "Expected values:\n" << rhsExpected;
        cout << "Actual values:\n" << rhsActual;
        cout << "Test dof ordering:\n" << *(elemType->testOrderPtr);
    }

    return allSuccess(success);
}

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

//.........这里部分代码省略.........
  if (rank==0)
    cout << "Integral of pressure: " << p_avg << endl;

  // integrate massFlux over each element (a test):
  // fake a new bilinear form so we can integrate against 1
  VarPtr testOne = varFactory.testVar("1",CONSTANT_SCALAR);
  BFPtr fakeBF = Teuchos::rcp( new BF(varFactory) );
  LinearTermPtr massFluxTerm = massFlux * testOne;

  CellTopoPtrLegacy quadTopoPtr = Teuchos::rcp(new shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >() ));
  DofOrderingFactory dofOrderingFactory(fakeBF);
  int fakeTestOrder = H1Order;
  DofOrderingPtr testOrdering = dofOrderingFactory.testOrdering(fakeTestOrder, *quadTopoPtr);

  int testOneIndex = testOrdering->getDofIndex(testOne->ID(),0);
  vector< ElementTypePtr > elemTypes = mesh->elementTypes(); // global element types
  map<int, double> massFluxIntegral; // cellID -> integral
  double maxMassFluxIntegral = 0.0;
  double totalMassFlux = 0.0;
  double totalAbsMassFlux = 0.0;
  double maxCellMeasure = 0;
  double minCellMeasure = 1;
  for (vector< ElementTypePtr >::iterator elemTypeIt = elemTypes.begin(); elemTypeIt != elemTypes.end(); elemTypeIt++)
  {
    ElementTypePtr elemType = *elemTypeIt;
    vector< ElementPtr > elems = mesh->elementsOfTypeGlobal(elemType);
    vector<GlobalIndexType> cellIDs;
    for (int i=0; i<elems.size(); i++)
    {
      cellIDs.push_back(elems[i]->cellID());
    }
    FieldContainer<double> physicalCellNodes = mesh->physicalCellNodesGlobal(elemType);
    BasisCachePtr basisCache = Teuchos::rcp( new BasisCache(elemType,mesh,polyOrder) ); // enrich by trial space order
    basisCache->setPhysicalCellNodes(physicalCellNodes,cellIDs,true); // true: create side caches
    FieldContainer<double> cellMeasures = basisCache->getCellMeasures();
    FieldContainer<double> fakeRHSIntegrals(elems.size(),testOrdering->totalDofs());
    massFluxTerm->integrate(fakeRHSIntegrals,testOrdering,basisCache,true); // true: force side evaluation
    //      cout << "fakeRHSIntegrals:\n" << fakeRHSIntegrals;
    for (int i=0; i<elems.size(); i++)
    {
      int cellID = cellIDs[i];
      // pick out the ones for testOne:
      massFluxIntegral[cellID] = fakeRHSIntegrals(i,testOneIndex);
    }
    // find the largest:
    for (int i=0; i<elems.size(); i++)
    {
      int cellID = cellIDs[i];
      maxMassFluxIntegral = max(abs(massFluxIntegral[cellID]), maxMassFluxIntegral);
    }
    for (int i=0; i<elems.size(); i++)
    {
      int cellID = cellIDs[i];
      maxCellMeasure = max(maxCellMeasure,cellMeasures(i));
      minCellMeasure = min(minCellMeasure,cellMeasures(i));
      maxMassFluxIntegral = max(abs(massFluxIntegral[cellID]), maxMassFluxIntegral);
      totalMassFlux += massFluxIntegral[cellID];
      totalAbsMassFlux += abs( massFluxIntegral[cellID] );
    }
  }
  if (rank==0)
  {
    cout << "largest mass flux: " << maxMassFluxIntegral << endl;
    cout << "total mass flux: " << totalMassFlux << endl;
    cout << "sum of mass flux absolute value: " << totalAbsMassFlux << endl;
    cout << "largest h: " << sqrt(maxCellMeasure) << endl;

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

  // v:
  // (sigma, grad v)_K - (sigma_hat_n, v)_dK - (u, beta dot grad v) + (u_hat * n dot beta, v)_dK
  bf->addTerm( sigma1, v->dx() );
  bf->addTerm( sigma2, v->dy() );
  bf->addTerm( -u, beta * v->grad());
  bf->addTerm( beta_n_u_minus_sigma_hat, v);

  // ==================== SET INITIAL GUESS ==========================
  mesh->registerSolution(backgroundFlow);

  map<int, Teuchos::RCP<AbstractFunction> > functionMap;
  functionMap[BurgersBilinearForm::U] = Teuchos::rcp(new InitialGuess());
  functionMap[BurgersBilinearForm::SIGMA_1] = Teuchos::rcp(new ZeroFunction());
  functionMap[BurgersBilinearForm::SIGMA_2] = Teuchos::rcp(new ZeroFunction());

  backgroundFlow->projectOntoMesh(functionMap);

  // ==================== END SET INITIAL GUESS ==========================
  // compare stiffness matrices for first linear step:
  int trialOrder = 1;
  pToAdd = 0;
  int testOrder = trialOrder + pToAdd;
  CellTopoPtr quadTopoPtr = Teuchos::rcp(new shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >() ));
  DofOrderingFactory dofOrderingFactory(bf);
  DofOrderingPtr testOrdering = dofOrderingFactory.testOrdering(testOrder, *quadTopoPtr);
  DofOrderingPtr trialOrdering = dofOrderingFactory.trialOrdering(trialOrder, *quadTopoPtr);

  int numCells = 1;
  // just use testOrdering for both trial and test spaces (we only use to define BasisCache)
  ElementTypePtr elemType  = Teuchos::rcp( new ElementType(trialOrdering, testOrdering, quadTopoPtr) );
  BasisCachePtr basisCache = Teuchos::rcp( new BasisCache(elemType) );
  quadPoints.resize(1,quadPoints.dimension(0),quadPoints.dimension(1));
  basisCache->setPhysicalCellNodes(quadPoints,vector<int>(1),true); // true: do create side cache
  FieldContainer<double> cellSideParities(numCells,quadTopoPtr->getSideCount());
  cellSideParities.initialize(1.0); // not worried here about neighbors actually having opposite parity -- just want the two BF implementations to agree...
  FieldContainer<double> expectedValues(numCells, testOrdering->totalDofs(), trialOrdering->totalDofs() );
  FieldContainer<double> actualValues(numCells, testOrdering->totalDofs(), trialOrdering->totalDofs() );
  oldBurgersBF->stiffnessMatrix(expectedValues, elemType, cellSideParities, basisCache);
  bf->stiffnessMatrix(actualValues, elemType, cellSideParities, basisCache);

  // compare beta's as well
  FieldContainer<double> expectedBeta = oldBurgersBF->getBeta(basisCache);
  Teuchos::Array<int> dim;
  expectedBeta.dimensions(dim);
  FieldContainer<double> actualBeta(dim);
  beta->values(actualBeta,basisCache);

  double tol = 1e-14;
  double maxDiff;
  if (rank == 0)
  {
    if ( ! TestSuite::fcsAgree(expectedBeta,actualBeta,tol,maxDiff) )
    {
      cout << "Test failed: old Burgers beta differs from new; maxDiff " << maxDiff << ".\n";
      cout << "Old beta: \n" << expectedBeta;
      cout << "New beta: \n" << actualBeta;
    }
    else
    {
      cout << "Old and new Burgers beta agree!!\n";
    }

    if ( ! TestSuite::fcsAgree(expectedValues,actualValues,tol,maxDiff) )
    {
      cout << "Test failed: old Burgers stiffness differs from new; maxDiff " << maxDiff << ".\n";

void ExactSolution::L2NormOfError(FieldContainer<double> &errorSquaredPerCell, Solution &solution, ElementTypePtr elemTypePtr, int trialID, int sideIndex, int cubDegree, double solutionLift) {
//  BasisCache(ElementTypePtr elemType, Teuchos::RCP<Mesh> mesh = Teuchos::rcp( (Mesh*) NULL ), bool testVsTest=false, int cubatureDegreeEnrichment = 0)

  DofOrdering dofOrdering = *(elemTypePtr->trialOrderPtr.get());
  BasisPtr basis = dofOrdering.getBasis(trialID,sideIndex);
  
  bool boundaryIntegral = solution.mesh()->bilinearForm()->isFluxOrTrace(trialID);
  
  BasisCachePtr basisCache;
  if (cubDegree <= 0) { // then take the default cub. degree
    basisCache = Teuchos::rcp( new BasisCache( elemTypePtr, solution.mesh() ) );
  } else {
    // we could eliminate the logic below if we just added BasisCache::setCubatureDegree()...
    // (the logic below is just to make the enriched cubature match the requested cubature degree...)
    int maxTrialDegree;
    if (!boundaryIntegral) {
      maxTrialDegree = elemTypePtr->trialOrderPtr->maxBasisDegreeForVolume();
    } else {
      maxTrialDegree = elemTypePtr->trialOrderPtr->maxBasisDegree(); // generally, this will be the trace degree
    }
    int maxTestDegree = elemTypePtr->testOrderPtr->maxBasisDegree();
    int cubDegreeEnrichment = max(cubDegree - (maxTrialDegree + maxTestDegree), 0);
    basisCache = Teuchos::rcp( new BasisCache( elemTypePtr, solution.mesh(), false, cubDegreeEnrichment) );
  }
  
  // much of this code is the same as what's in the volume integration in computeStiffness...
  FieldContainer<double> physicalCellNodes = solution.mesh()->physicalCellNodes(elemTypePtr);
  vector<GlobalIndexType> cellIDs = solution.mesh()->cellIDsOfType(elemTypePtr);
  basisCache->setPhysicalCellNodes(physicalCellNodes, cellIDs, true);
  
  if (boundaryIntegral) {
    basisCache = basisCache->getSideBasisCache(sideIndex);
  }
  
  FieldContainer<double> weightedMeasure = basisCache->getWeightedMeasures();
  FieldContainer<double> weightedErrorSquared;
  
  int numCells = basisCache->getPhysicalCubaturePoints().dimension(0);
  int numCubPoints = basisCache->getPhysicalCubaturePoints().dimension(1);
  int spaceDim = basisCache->getPhysicalCubaturePoints().dimension(2);
  
  Teuchos::Array<int> dimensions;
  dimensions.push_back(numCells);
  dimensions.push_back(numCubPoints);
  
  int basisRank = BasisFactory::basisFactory()->getBasisRank(basis);
  if (basisRank==1) {
    dimensions.push_back(spaceDim);
  }
  
  FieldContainer<double> computedValues(dimensions);
  FieldContainer<double> exactValues(dimensions);
  
  if (solutionLift != 0.0) {
    int size = computedValues.size();
    for (int i=0; i<size; i++) {
      computedValues[i] += solutionLift;
    }
  }
  
  solution.solutionValues(computedValues, trialID, basisCache);
  this->solutionValues(exactValues, trialID, basisCache);
  
//  cout << "ExactSolution: exact values:\n" << exactValues;
//  cout << "ExactSolution: computed values:\n" << computedValues;
  
  FieldContainer<double> errorSquared(numCells,numCubPoints);
  
  squaredDifference(errorSquared,computedValues,exactValues);
  
  weightedErrorSquared.resize(numCells,numCubPoints);
  for (int cellIndex=0; cellIndex<numCells; cellIndex++) {
    for (int ptIndex=0; ptIndex<numCubPoints; ptIndex++) {
      // following two lines for viewing in the debugger:
      double weight = weightedMeasure(cellIndex,ptIndex);
      double errorSquaredVal = errorSquared(cellIndex,ptIndex);
      weightedErrorSquared(cellIndex,ptIndex) = errorSquared(cellIndex,ptIndex) * weightedMeasure(cellIndex,ptIndex);
    }
  }
  
  // compute the integral
  errorSquaredPerCell.initialize(0.0);
  int numPoints = weightedErrorSquared.dimension(1);
  for (int cellIndex=0; cellIndex<numCells; cellIndex++) {
    for (int ptIndex=0; ptIndex<numPoints; ptIndex++) {
      errorSquaredPerCell(cellIndex) += weightedErrorSquared(cellIndex,ptIndex);
    }
  }
}

bool LobattoBasisTests::testSimpleStiffnessMatrix() {
  bool success = true;
  
  int rank = Teuchos::GlobalMPISession::getRank();
  
  VarFactory varFactory;
  VarPtr u = varFactory.fieldVar("u");
  VarPtr un = varFactory.fluxVar("un_hat");
  VarPtr v = varFactory.testVar("v", HGRAD);
  
  BFPtr bf = Teuchos::rcp( new BF(varFactory) );
  vector<double> beta;
  beta.push_back(1);
  beta.push_back(1);
  bf->addTerm(beta * u, v->grad());
  bf->addTerm(un, v);
  
  DofOrderingPtr trialSpace = Teuchos::rcp( new DofOrdering );
  DofOrderingPtr testSpace = Teuchos::rcp( new DofOrdering );
  
  const int numSides = 4;
  const int spaceDim = 2;
  
  int fieldOrder = 3;
  int testOrder = fieldOrder+2;
  BasisPtr fieldBasis = Camellia::intrepidQuadHGRAD(fieldOrder);
  BasisPtr fluxBasis = Camellia::intrepidLineHGRAD(fieldOrder);
  trialSpace->addEntry(u->ID(), fieldBasis, fieldBasis->rangeRank());
  for (int i=0; i<numSides; i++) {
    trialSpace->addEntry(un->ID(), fluxBasis, fluxBasis->rangeRank(), i);
  }
  
  BasisPtr testBasis = Camellia::lobattoQuadHGRAD(testOrder+1,false); // +1 because it lives in HGRAD
  testSpace->addEntry(v->ID(), testBasis, testBasis->rangeRank());
  
  int numTrialDofs = trialSpace->totalDofs();
  int numTestDofs = testSpace->totalDofs();
  int numCells = 1;
  
  FieldContainer<double> cellNodes(numCells,numSides,spaceDim);
  cellNodes(0,0,0) = 0;
  cellNodes(0,0,1) = 0;
  cellNodes(0,1,0) = 1;
  cellNodes(0,1,1) = 0;
  cellNodes(0,2,0) = 1;
  cellNodes(0,2,1) = 1;
  cellNodes(0,3,0) = 0;
  cellNodes(0,3,1) = 1;
  
  FieldContainer<double> stiffness(numCells,numTestDofs,numTrialDofs);
  
  FieldContainer<double> cellSideParities(numCells,numSides);
  cellSideParities.initialize(1.0);
  
  Teuchos::RCP<shards::CellTopology> quad_4 = Teuchos::rcp( new shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >() ) );
  Teuchos::RCP<ElementType> elemType = Teuchos::rcp( new ElementType(trialSpace, testSpace, quad_4));
  
  BasisCachePtr basisCache = Teuchos::rcp( new BasisCache(elemType) );
  vector<GlobalIndexType> cellIDs;
  cellIDs.push_back(0);
  basisCache->setPhysicalCellNodes(cellNodes, cellIDs, true);
  bf->stiffnessMatrix(stiffness, elemType, cellSideParities, basisCache);

  // TODO: finish this test
  
//  cout << stiffness;
  if (rank==0)
    cout << "Warning: testSimpleStiffnessMatrix() unfinished.\n";
  
  return success;
}

//.........这里部分代码省略.........
    else
    {
      return false;
    }
  }
  pair<GlobalIndexType, unsigned> neighborInfo = fineCell->getNeighborInfo(sideOrdinal, meshTopo);
  if (neighborInfo.first == -1)
  {
    // boundary
    return false;
  }
  CellPtr neighborCell = meshTopo->getCell(neighborInfo.first);
  if (neighborCell->isParent(meshTopo))
  {
    return false; // fineCell isn't the finer of the two...
  }

  CellTopoPtr fineSideTopo = fineCell->topology()->getSubcell(sideDim, sideOrdinal);

  CubatureFactory cubFactory;
  int cubDegree = 4; // fairly arbitrary choice: enough to get a decent number of test points...
  Teuchos::RCP<Cubature<double> > fineSideTopoCub = cubFactory.create(fineSideTopo, cubDegree);

  int numCubPoints = fineSideTopoCub->getNumPoints();

  FieldContainer<double> cubPoints(numCubPoints, sideDim);
  FieldContainer<double> cubWeights(numCubPoints); // we neglect these...

  fineSideTopoCub->getCubature(cubPoints, cubWeights);

  FieldContainer<double> sideRefCellNodes(fineSideTopo->getNodeCount(),sideDim);
  CamelliaCellTools::refCellNodesForTopology(sideRefCellNodes, fineSideTopo);

  int numTestPoints = numCubPoints + fineSideTopo->getNodeCount();

  FieldContainer<double> testPoints(numTestPoints, sideDim);
  for (int ptOrdinal=0; ptOrdinal<testPoints.dimension(0); ptOrdinal++)
  {
    if (ptOrdinal<fineSideTopo->getNodeCount())
    {
      for (int d=0; d<sideDim; d++)
      {
        testPoints(ptOrdinal,d) = sideRefCellNodes(ptOrdinal,d);
      }
    }
    else
    {
      for (int d=0; d<sideDim; d++)
      {
        testPoints(ptOrdinal,d) = cubPoints(ptOrdinal-fineSideTopo->getNodeCount(),d);
      }
    }
  }

  fineSideRefPoints = testPoints;
  fineCellRefPoints.resize(numTestPoints, spaceDim);

  CamelliaCellTools::mapToReferenceSubcell(fineCellRefPoints, testPoints, sideDim, sideOrdinal, fineCell->topology());

  CellTopoPtr coarseSideTopo = neighborCell->topology()->getSubcell(sideDim, neighborInfo.second);

  unsigned fineSideAncestorPermutation = fineCell->ancestralPermutationForSubcell(sideDim, sideOrdinal, meshTopo);
  unsigned coarseSidePermutation = neighborCell->subcellPermutation(sideDim, neighborInfo.second);

  unsigned coarseSideAncestorPermutationInverse = CamelliaCellTools::permutationInverse(coarseSideTopo, coarseSidePermutation);

  unsigned composedPermutation = CamelliaCellTools::permutationComposition(coarseSideTopo, coarseSideAncestorPermutationInverse, fineSideAncestorPermutation); // goes from coarse ordering to fine.

  RefinementBranch fineRefBranch = fineCell->refinementBranchForSide(sideOrdinal, meshTopo);

  FieldContainer<double> fineSideNodes(fineSideTopo->getNodeCount(), sideDim);  // relative to the ancestral cell whose neighbor is compatible
  if (fineRefBranch.size() == 0)
  {
    CamelliaCellTools::refCellNodesForTopology(fineSideNodes, coarseSideTopo, composedPermutation);
  }
  else
  {
    FieldContainer<double> ancestralSideNodes(coarseSideTopo->getNodeCount(), sideDim);
    CamelliaCellTools::refCellNodesForTopology(ancestralSideNodes, coarseSideTopo, composedPermutation);

    RefinementBranch fineSideRefBranch = RefinementPattern::sideRefinementBranch(fineRefBranch, sideOrdinal);
    fineSideNodes = RefinementPattern::descendantNodes(fineSideRefBranch, ancestralSideNodes);
  }

  BasisCachePtr sideTopoCache = Teuchos::rcp( new BasisCache(fineSideTopo, 1, false) );
  sideTopoCache->setRefCellPoints(testPoints);

  // add cell dimension
  fineSideNodes.resize(1,fineSideNodes.dimension(0), fineSideNodes.dimension(1));
  sideTopoCache->setPhysicalCellNodes(fineSideNodes, vector<GlobalIndexType>(), false);
  coarseSideRefPoints = sideTopoCache->getPhysicalCubaturePoints();

  // strip off the cell dimension:
  coarseSideRefPoints.resize(coarseSideRefPoints.dimension(1),coarseSideRefPoints.dimension(2));

  coarseCellRefPoints.resize(numTestPoints,spaceDim);
  CamelliaCellTools::mapToReferenceSubcell(coarseCellRefPoints, coarseSideRefPoints, sideDim, neighborInfo.second, neighborCell->topology());

  return true; // containers filled....
}

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

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

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

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

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

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

bool LinearTermTests::testRieszInversionAsProjection()
{
  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 = 2;
  int pToAdd = 2;

  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 = 2;
  int horizontalCells = nCells, verticalCells = nCells;
  // create a new mesh:
  MeshPtr 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 = myMesh->cellIDsOfTypeGlobal(elemType);
  bool createSideCacheToo = true;

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

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

  FunctionPtr x = Function::xn(1);
  FunctionPtr y = Function::yn(1);
  FunctionPtr testFxn1 = x;
  FunctionPtr testFxn2 = y;
  FunctionPtr fxnToProject = x * y + 1.0;

  integrand->addTerm(fxnToProject * v);

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

  Teuchos::RCP<RieszRep> riesz = Teuchos::rcp(new RieszRep(myMesh, ip_L2, integrand));
  riesz->computeRieszRep();

  FunctionPtr rieszFxn = RieszRep::repFunction(v,riesz);
  int numCells = basisCache->getPhysicalCubaturePoints().dimension(0);
  int numPts = basisCache->getPhysicalCubaturePoints().dimension(1);

  FieldContainer<double> valProject( numCells, numPts );
  FieldContainer<double> valExpected( numCells, numPts );

  rieszFxn->values(valProject,basisCache);
  fxnToProject->values(valExpected,basisCache);

//  int rank = Teuchos::GlobalMPISession::getRank();
//  if (rank==0) cout << "physicalCubaturePoints:\n" << basisCache->getPhysicalCubaturePoints();

  double maxDiff;
  double tol = 1e-12;
  success = TestSuite::fcsAgree(valProject,valExpected,tol,maxDiff);
  if (success==false)
//.........这里部分代码省略.........

//.........这里部分代码省略.........
        if (transient)
          outfile << "TransientConfusion_" << refIndex << "_" << timestepCount;
        else
          outfile << "TransientConfusion_" << refIndex;
        solution->writeToVTK(outfile.str(), 5);
      }

      //////////////////////////////////////////////////////////////////////////
      // Check conservation by testing against one
      //////////////////////////////////////////////////////////////////////////
      VarPtr testOne = varFactory.testVar("1", CONSTANT_SCALAR);
      // Create a fake bilinear form for the testing
      BFPtr fakeBF = Teuchos::rcp( new BF(varFactory) );
      // Define our mass flux
      FunctionPtr flux_current_time = Teuchos::rcp( new PreviousSolutionFunction(solution, beta_n_u_minus_sigma_n) );
      FunctionPtr delta_u = Teuchos::rcp( new PreviousSolutionFunction(flowResidual, u) );
      LinearTermPtr surfaceFlux = -1.0 * flux_current_time * testOne;
      LinearTermPtr volumeChange = invDt * delta_u * testOne;
      LinearTermPtr massFluxTerm;
      if (transient)
      {
        massFluxTerm = volumeChange;
        // massFluxTerm->addTerm(surfaceFlux);
      }
      else
      {
        massFluxTerm = surfaceFlux;
      }
      // cout << "surface case = " << surfaceFlux->summands()[0].first->boundaryValueOnly() << " volume case = " << volumeChange->summands()[0].first->boundaryValueOnly() << endl;

      // FunctionPtr massFlux= Teuchos::rcp( new PreviousSolutionFunction(solution, beta_n_u_minus_sigma_n) );
      // LinearTermPtr massFluxTerm = massFlux * testOne;

      Teuchos::RCP<shards::CellTopology> quadTopoPtr = Teuchos::rcp(new shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >() ));
      DofOrderingFactory dofOrderingFactory(fakeBF);
      int fakeTestOrder = H1Order;
      DofOrderingPtr testOrdering = dofOrderingFactory.testOrdering(fakeTestOrder, *quadTopoPtr);

      int testOneIndex = testOrdering->getDofIndex(testOne->ID(),0);
      vector< ElementTypePtr > elemTypes = mesh->elementTypes(); // global element types
      map<int, double> massFluxIntegral; // cellID -> integral
      double maxMassFluxIntegral = 0.0;
      double totalMassFlux = 0.0;
      double totalAbsMassFlux = 0.0;
      for (vector< ElementTypePtr >::iterator elemTypeIt = elemTypes.begin(); elemTypeIt != elemTypes.end(); elemTypeIt++) 
      {
        ElementTypePtr elemType = *elemTypeIt;
        vector< ElementPtr > elems = mesh->elementsOfTypeGlobal(elemType);
        vector<int> cellIDs;
        for (int i=0; i<elems.size(); i++) {
          cellIDs.push_back(elems[i]->cellID());
        }
        FieldContainer<double> physicalCellNodes = mesh->physicalCellNodesGlobal(elemType);
        BasisCachePtr basisCache = Teuchos::rcp( new BasisCache(elemType,mesh) );
        basisCache->setPhysicalCellNodes(physicalCellNodes,cellIDs,true); // true: create side caches
        FieldContainer<double> cellMeasures = basisCache->getCellMeasures();
        FieldContainer<double> fakeRHSIntegrals(elems.size(),testOrdering->totalDofs());
        massFluxTerm->integrate(fakeRHSIntegrals,testOrdering,basisCache,true); // true: force side evaluation
        for (int i=0; i<elems.size(); i++) {
          int cellID = cellIDs[i];
          // pick out the ones for testOne:
          massFluxIntegral[cellID] = fakeRHSIntegrals(i,testOneIndex);
        }
        // find the largest:
        for (int i=0; i<elems.size(); i++) {
          int cellID = cellIDs[i];
          maxMassFluxIntegral = max(abs(massFluxIntegral[cellID]), maxMassFluxIntegral);
        }
        for (int i=0; i<elems.size(); i++) {
          int cellID = cellIDs[i];
          maxMassFluxIntegral = max(abs(massFluxIntegral[cellID]), maxMassFluxIntegral);
          totalMassFlux += massFluxIntegral[cellID];
          totalAbsMassFlux += abs( massFluxIntegral[cellID] );
        }
      }

      // Print results from processor with rank 0
      if (rank == 0)
      {
        cout << "largest mass flux: " << maxMassFluxIntegral << endl;
        cout << "total mass flux: " << totalMassFlux << endl;
        cout << "sum of mass flux absolute value: " << totalAbsMassFlux << endl;
      }

      prevTimeFlow->setSolution(solution); // reset previous time solution to current time sol
      timestepCount++;
    }

    if (refIndex < numRefs){
      if (rank==0){
        cout << "Performing refinement number " << refIndex << endl;
      }     
      refinementStrategy->refine(rank==0);    
      // RESET solution every refinement - make sure discretization error doesn't creep in
      // prevTimeFlow->projectOntoMesh(functionMap);
    }
  }

  return 0;
}