C++ DofOrderingPtr类代码示例

map<GlobalIndexType,FieldContainer<double> > RieszRep::integrateRHS() {
  // NVR: changed this to only return integrated values for rank-local cells.

  map<GlobalIndexType,FieldContainer<double> > cellRHS;
  set<GlobalIndexType> cellIDs = _mesh->cellIDsInPartition();
  for (set<GlobalIndexType>::iterator cellIDIt=cellIDs.begin(); cellIDIt !=cellIDs.end(); cellIDIt++){
    GlobalIndexType cellID = *cellIDIt;
    ElementTypePtr elemTypePtr = _mesh->getElementType(cellID);
    DofOrderingPtr testOrderingPtr = elemTypePtr->testOrderPtr;
    int numTestDofs = testOrderingPtr->totalDofs();

    int cubEnrich = 0; // set to zero for release
    BasisCachePtr basisCache = BasisCache::basisCacheForCell(_mesh,cellID,true,cubEnrich);

    FieldContainer<double> rhsValues(1,numTestDofs);
    _rhs->integrate(rhsValues, testOrderingPtr, basisCache);

    FieldContainer<double> rhsVals(numTestDofs);
    for (int i = 0;i<numTestDofs;i++){
      rhsVals(i) = rhsValues(0,i);
    }
    cellRHS[cellID] = rhsVals;
  }
  return cellRHS;
}

double MeshUtilities::computeMaxLocalConditionNumber(Teuchos::RCP< DPGInnerProduct > ip, MeshPtr mesh, bool jacobiScaling, string sparseFileToWriteTo) {
  int rank = Teuchos::GlobalMPISession::getRank();
  int numProcs = Teuchos::GlobalMPISession::getNProc();
  vector< ElementPtr > elements = mesh->elementsInPartition(rank);

//  cout << "Checking condition numbers on rank " << rank << endl;
  
  FieldContainer<double> maxConditionNumberIPMatrix;
  int maxCellID = -1;
  double myMaxConditionNumber = -1;
  for (vector< ElementPtr >::iterator elemIt = elements.begin(); elemIt != elements.end(); elemIt++) {
    int cellID = (*elemIt)->cellID();
    bool testVsTest = true;
    BasisCachePtr cellBasisCache = BasisCache::basisCacheForCell(mesh, cellID, testVsTest);
    DofOrderingPtr testSpace = (*elemIt)->elementType()->testOrderPtr;
    int testDofs = testSpace->totalDofs();
    int numCells = 1;
    FieldContainer<double> innerProductMatrix(numCells,testDofs,testDofs);
    ip->computeInnerProductMatrix(innerProductMatrix, testSpace, cellBasisCache);
    // reshape:
    innerProductMatrix.resize(testDofs,testDofs);
    if (jacobiScaling)
      SerialDenseMatrixUtility::jacobiScaleMatrix(innerProductMatrix);
//    double conditionNumber = SerialDenseMatrixUtility::estimate1NormConditionNumber(innerProductMatrix);
    double conditionNumber = SerialDenseMatrixUtility::estimate2NormConditionNumber(innerProductMatrix);
    if (conditionNumber > myMaxConditionNumber) {
      myMaxConditionNumber = conditionNumber;
      maxConditionNumberIPMatrix = innerProductMatrix;
      maxCellID = cellID;
    } else if (maxConditionNumberIPMatrix.size()==0) {
      // could be that the estimation failed--we still want a copy of the matrix written to file.
      maxConditionNumberIPMatrix = innerProductMatrix;
    }
  }
//  cout << "Determined condition numbers on rank " << rank << endl;
  FieldContainer<double> maxConditionNumbers(numProcs);
  maxConditionNumbers[rank] = myMaxConditionNumber;
  MPIWrapper::entryWiseSum(maxConditionNumbers);
  
  double maxConditionNumber = maxConditionNumbers[0];
  int maxConditionNumberOwner = 0; // the MPI node with the max condition number
  for (int i=1; i<numProcs; i++) {
    if (maxConditionNumber < maxConditionNumbers[i]) {
      maxConditionNumber = maxConditionNumbers[i];
      maxConditionNumberOwner = i;
    }
  }
  
  if (rank==maxConditionNumberOwner) { // owner is responsible for writing to file
//    cout << "max condition number is on rank " << rank << endl;
    if (sparseFileToWriteTo.length() > 0) {
      if (maxConditionNumberIPMatrix.size() > 0) {
        DataIO::writeMatrixToSparseDataFile(maxConditionNumberIPMatrix, sparseFileToWriteTo);
      }
    }
  }
//  cout << "max condition number occurs in cellID " << maxCellID << endl;
  return maxConditionNumber;
}

void BilinearFormUtility<Scalar>::computeStiffnessMatrixForCell(FieldContainer<Scalar> &stiffness, Teuchos::RCP<Mesh> mesh, int cellID)
{
    DofOrderingPtr trialOrder = mesh->getElementType(cellID)->trialOrderPtr;
    DofOrderingPtr testOrder  = mesh->getElementType(cellID)->testOrderPtr;
    CellTopoPtr     cellTopo  = mesh->getElementType(cellID)->cellTopoPtr;
    FieldContainer<double> physicalCellNodes = mesh->physicalCellNodesForCell(cellID);
    FieldContainer<double> cellSideParities  = mesh->cellSideParitiesForCell(cellID);
    int numCells = 1;
    stiffness.resize(numCells,testOrder->totalDofs(),trialOrder->totalDofs());
    computeStiffnessMatrix(stiffness,mesh->bilinearForm(),trialOrder,testOrder,cellTopo,physicalCellNodes,cellSideParities);
}

// computes riesz representation over a single element - map is from int (testID) to FieldContainer of values (sized cellIndex, numPoints)
void RieszRep::computeRepresentationValues(FieldContainer<double> &values, int testID, IntrepidExtendedTypes::EOperatorExtended op, BasisCachePtr basisCache){

  if (_repsNotComputed){
    cout << "Computing riesz rep dofs" << endl;
    computeRieszRep();
  }

  int spaceDim = _mesh->getTopology()->getSpaceDim();
  int numCells = values.dimension(0);
  int numPoints = values.dimension(1);
  vector<GlobalIndexType> cellIDs = basisCache->cellIDs();

  // all elems coming in should be of same type
  ElementPtr elem = _mesh->getElement(cellIDs[0]);
  ElementTypePtr elemTypePtr = elem->elementType();   
  DofOrderingPtr testOrderingPtr = elemTypePtr->testOrderPtr;
  CellTopoPtrLegacy cellTopoPtr = elemTypePtr->cellTopoPtr;
  int numTestDofsForVarID = testOrderingPtr->getBasisCardinality(testID, 0);
  BasisPtr testBasis = testOrderingPtr->getBasis(testID);
  
  bool testBasisIsVolumeBasis = (spaceDim == testBasis->domainTopology()->getDimension());  
  bool useCubPointsSideRefCell = testBasisIsVolumeBasis && basisCache->isSideCache();
  
  Teuchos::RCP< const FieldContainer<double> > transformedBasisValues = basisCache->getTransformedValues(testBasis,op,useCubPointsSideRefCell);
  
  int rank = values.rank() - 2; // if values are shaped as (C,P), scalar...
  if (rank > 1) {
    cout << "ranks greater than 1 not presently supported...\n";
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "ranks greater than 1 not presently supported...");
  }
  
//  Camellia::print("cellIDs",cellIDs);
  
  values.initialize(0.0);
  for (int cellIndex = 0;cellIndex<numCells;cellIndex++){
    int cellID = cellIDs[cellIndex];
    for (int j = 0;j<numTestDofsForVarID;j++) {
      int dofIndex = testOrderingPtr->getDofIndex(testID, j);
      for (int i = 0;i<numPoints;i++) {
        if (rank==0) {
          double basisValue = (*transformedBasisValues)(cellIndex,j,i);
          values(cellIndex,i) += basisValue*_rieszRepDofsGlobal[cellID](dofIndex);
        } else {
          for (int d = 0; d<spaceDim; d++) {
            double basisValue = (*transformedBasisValues)(cellIndex,j,i,d);
            values(cellIndex,i,d) += basisValue*_rieszRepDofsGlobal[cellID](dofIndex);
          }
        }
      }
    }
  }
//  TestSuite::serializeOutput("rep values", values);
}

bool checkLTSumConsistency(LinearTermPtr a, LinearTermPtr b, DofOrderingPtr dofOrdering, BasisCachePtr basisCache)
{
  double tol = 1e-14;

  int numCells = basisCache->cellIDs().size();
  int numDofs = dofOrdering->totalDofs();
  bool forceBoundaryTerm = false;
  FieldContainer<double> aValues(numCells,numDofs), bValues(numCells,numDofs), sumValues(numCells,numDofs);
  a->integrate(aValues,dofOrdering,basisCache,forceBoundaryTerm);
  b->integrate(bValues,dofOrdering,basisCache,forceBoundaryTerm);
  (a+b)->integrate(sumValues, dofOrdering, basisCache, forceBoundaryTerm);

  int size = aValues.size();

  for (int i=0; i<size; i++)
  {
    double expectedValue = aValues[i] + bValues[i];
    double diff = abs( expectedValue - sumValues[i] );
    if (diff > tol)
    {
      return false;
    }
  }
  return true;
}

map< int, vector<DofInfo> > constructGlobalDofToLocalDofInfoMap(MeshPtr mesh)
{
  // go through the mesh as a whole, and collect info for each dof
  map< int, vector<DofInfo> > infoMap;
  DofInfo info;
  set<GlobalIndexType> activeCellIDs = mesh->getActiveCellIDsGlobal();
  for (set<GlobalIndexType>::iterator cellIt = activeCellIDs.begin(); cellIt != activeCellIDs.end(); cellIt++)
  {
    GlobalIndexType cellID = *cellIt;
    info.cellID = cellID;
    ElementPtr element = mesh->getElement(cellID);
    DofOrderingPtr trialOrder = element->elementType()->trialOrderPtr;
    set<int> trialIDs = trialOrder->getVarIDs();
    info.totalDofs = trialOrder->totalDofs();
    for (set<int>::iterator trialIt=trialIDs.begin(); trialIt != trialIDs.end(); trialIt++)
    {
      info.trialID = *trialIt;
      const vector<int>* sidesForVar = &trialOrder->getSidesForVarID(info.trialID);
      for (vector<int>::const_iterator sideIt = sidesForVar->begin(); sideIt != sidesForVar->end(); sideIt++)
      {
        int sideIndex = *sideIt;
        info.sideIndex = sideIndex;
        info.basisCardinality = trialOrder->getBasisCardinality(info.trialID, info.sideIndex);
        for (int basisOrdinal=0; basisOrdinal < info.basisCardinality; basisOrdinal++)
        {
          info.basisOrdinal = basisOrdinal;
          info.localDofIndex = trialOrder->getDofIndex(info.trialID, info.basisOrdinal, info.sideIndex);
          pair<int, int> localDofIndexKey = make_pair(info.cellID, info.localDofIndex);
          int globalDofIndex = mesh->getLocalToGlobalMap().find(localDofIndexKey)->second;
//          cout << "(" << info.cellID << "," << info.localDofIndex << ") --> " << globalDofIndex << endl;
          infoMap[globalDofIndex].push_back(info);
        }
      }
    }
  }
  return infoMap;
}

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

  vector< FunctionPtr > vectorFxns;

  double xValue = 3, yValue = 4;
  FunctionPtr simpleVector = Function::vectorize(Function::constant(xValue), Function::constant(yValue));
  vectorFxns.push_back(simpleVector);
  FunctionPtr x = Function::xn(1);
  FunctionPtr y = Function::yn(1);
  vectorFxns.push_back( Function::vectorize(x*x, x*y) );

  VGPStokesFormulation vgpStokes = VGPStokesFormulation(1.0);
  BFPtr bf = vgpStokes.bf();

  int h1Order = 1;
  MeshPtr mesh = MeshFactory::quadMesh(bf, h1Order);

  int cellID=0; // the only cell
  BasisCachePtr basisCache = BasisCache::basisCacheForCell(mesh, cellID);

  for (int i=0; i<vectorFxns.size(); i++)
  {
    FunctionPtr vectorFxn_i = vectorFxns[i];
    for (int j=0; j<vectorFxns.size(); j++)
    {
      FunctionPtr vectorFxn_j = vectorFxns[j];
      FunctionPtr dotProduct = vectorFxn_i * vectorFxn_j;
      FunctionPtr expectedDotProduct = vectorFxn_i->x() * vectorFxn_j->x() + vectorFxn_i->y() * vectorFxn_j->y();
      if (! expectedDotProduct->equals(dotProduct, basisCache))
      {
        cout << "testValuesDottedWithTensor() failed: expected dot product does not match dotProduct.\n";
        success = false;
        double tol = 1e-14;
        reportFunctionValueDifferences(dotProduct, expectedDotProduct, basisCache, tol);
      }
    }
  }

  // now, let's try the same thing, but for a LinearTerm dot product
  VarFactoryPtr vf = VarFactory::varFactory();
  VarPtr v = vf->testVar("v", HGRAD);

  DofOrderingPtr dofOrdering = Teuchos::rcp( new DofOrdering(CellTopology::quad()) );
  shards::CellTopology quad_4(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
  BasisPtr basis = BasisFactory::basisFactory()->getBasis(h1Order, quad_4.getKey(), Camellia::FUNCTION_SPACE_HGRAD);
  dofOrdering->addEntry(v->ID(), basis, v->rank());

  int numCells = 1;
  int numFields = basis->getCardinality();

  for (int i=0; i<vectorFxns.size(); i++)
  {
    FunctionPtr f_i = vectorFxns[i];
    LinearTermPtr lt_i = f_i * v;
    LinearTermPtr lt_i_x = f_i->x() * v;
    LinearTermPtr lt_i_y = f_i->y() * v;
    for (int j=0; j<vectorFxns.size(); j++)
    {
      FunctionPtr f_j = vectorFxns[j];
      LinearTermPtr lt_j = f_j * v;
      LinearTermPtr lt_j_x = f_j->x() * v;
      LinearTermPtr lt_j_y = f_j->y() * v;
      FieldContainer<double> values(numCells,numFields,numFields);
      lt_i->integrate(values, dofOrdering, lt_j, dofOrdering, basisCache);
      FieldContainer<double> values_expected(numCells,numFields,numFields);
      lt_i_x->integrate(values_expected,dofOrdering,lt_j_x,dofOrdering,basisCache);
      lt_i_y->integrate(values_expected,dofOrdering,lt_j_y,dofOrdering,basisCache);
      double tol = 1e-14;
      double maxDiff = 0;
      if (!fcsAgree(values, values_expected, tol, maxDiff))
      {
        cout << "FunctionTests::testValuesDottedWithTensor: ";
        cout << "dot product and sum of the products of scalar components differ by maxDiff " << maxDiff;
        cout << " in LinearTerm::integrate().\n";
        success = false;
      }
    }
  }

//  // finally, let's try the same sort of thing, but now with a vector-valued basis
//  BasisPtr vectorBasisTemp = BasisFactory::basisFactory()->getBasis(h1Order, quad_4.getKey(), Camellia::FUNCTION_SPACE_VECTOR_HGRAD);
//  VectorBasisPtr vectorBasis = Teuchos::rcp( (VectorizedBasis<double, FieldContainer<double> > *)vectorBasisTemp.get(),false);
//
//  BasisPtr compBasis = vectorBasis->getComponentBasis();
//
//  // create a new v, and a new dofOrdering
//  VarPtr v_vector = vf->testVar("v_vector", VECTOR_HGRAD);
//  dofOrdering = Teuchos::rcp( new DofOrdering );
//  dofOrdering->addEntry(v_vector->ID(), vectorBasis, v_vector->rank());
//
//  DofOrderingPtr dofOrderingComp = Teuchos::rcp( new DofOrdering );
//  dofOrderingComp->addEntry(v->ID(), compBasis, v->rank());
//

  return success;
}

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

  ////////////////////////////////////////////////////////////////////
  // DEFINE REFINEMENT STRATEGY
  ////////////////////////////////////////////////////////////////////
  Teuchos::RCP<RefinementStrategy> refinementStrategy;
  refinementStrategy = Teuchos::rcp(new RefinementStrategy(solution,energyThreshold));

  ////////////////////////////////////////////////////////////////////
  // SOLVE
  ////////////////////////////////////////////////////////////////////

  for (int refIndex=0; refIndex<=numRefs; refIndex++)
  {
    double L2Update = 1e7;
    int iterCount = 0;
    while (L2Update > nonlinearRelativeEnergyTolerance && iterCount < maxNewtonIterations)
    {
      solution->solve();
      L2Update = solution->L2NormOfSolutionGlobal(u->ID());
      cout << "L2 Norm of Update = " << L2Update << endl;
      // backgroundFlow->clear();
      backgroundFlow->addSolution(solution, newtonStepSize);
      iterCount++;
    }
    cout << endl;

    // check conservation
    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 massFlux = Teuchos::rcp( new PreviousSolutionFunction(solution, fhat) );
    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] );
      }
    }
    if (rank==0)
    {
      cout << endl;
      cout << "largest mass flux: " << maxMassFluxIntegral << endl;
      cout << "total mass flux: " << totalMassFlux << endl;
      cout << "sum of mass flux absolute value: " << totalAbsMassFlux << endl;
      cout << endl;

      stringstream outfile;
      outfile << "burgers_" << refIndex;
      backgroundFlow->writeToVTK(outfile.str(), 5);
    }

    if (refIndex < numRefs)
      refinementStrategy->refine(rank==0); // print to console on rank 0
  }

  return 0;
}

void PenaltyMethodFilter::filter(FieldContainer<double> &localStiffnessMatrix, FieldContainer<double> &localRHSVector,
                                 BasisCachePtr basisCache, Teuchos::RCP<Mesh> mesh, Teuchos::RCP<BC> bc){ 
  
  // assumption: filter gets elements of all the same type  
  TEUCHOS_TEST_FOR_EXCEPTION(basisCache->cellIDs().size()==0,std::invalid_argument,"no cell IDs given to filter");
  
  ElementTypePtr elemTypePtr = mesh->getElement(basisCache->cellIDs()[0])->elementType();
  int numCells = localStiffnessMatrix.dimension(0);
  
  DofOrderingPtr trialOrderPtr = elemTypePtr->trialOrderPtr;
  
  unsigned numSides = CamelliaCellTools::getSideCount( *elemTypePtr->cellTopoPtr );
  // only allows for L2 inner products at the moment. 
  IntrepidExtendedTypes::EOperatorExtended trialOperator =  IntrepidExtendedTypes::OP_VALUE;
	
  // loop over sides first 
  for (unsigned int sideIndex = 0; sideIndex<numSides; sideIndex++){
    
    // GET INTEGRATION INFO - get cubature points and side normals to send to Constraints (Cell,Point, spaceDim)
    FieldContainer<double> sideCubPoints = basisCache->getPhysicalCubaturePointsForSide(sideIndex);
    FieldContainer<double> sideNormals = basisCache->getSideUnitNormals(sideIndex);        
    
    int numPts = sideCubPoints.dimension(1);
    
    // GET CONSTRAINT INFO
    vector<map<int, FieldContainer<double> > > constrCoeffsVector;
    vector<FieldContainer<double> > constraintValuesVector;
    vector<FieldContainer<bool> > imposeHereVector;
    _constraints->getConstraints(sideCubPoints,sideNormals,constrCoeffsVector,constraintValuesVector);
    
    //loop thru constraints
    int i = 0;
    for (vector<map<int,FieldContainer<double> > >::iterator constrIt = constrCoeffsVector.begin();
         constrIt !=constrCoeffsVector.end(); constrIt++) {
      map<int,FieldContainer<double> > constrCoeffs = *constrIt;
      FieldContainer<double> constrValues = constraintValuesVector[i];
      i++;
      
      double penaltyParameter = 1e7; // (single_precision)^(-1) - perhaps have this computed relative to terms in the matrix?
      
      for (map<int,FieldContainer<double> >::iterator constrTestIDIt = constrCoeffs.begin();
           constrTestIDIt !=constrCoeffs.end(); constrTestIDIt++) {
        pair<int,FieldContainer<double> > constrTestPair = *constrTestIDIt;
        int testTrialID = constrTestPair.first;
        
        // get basis to integrate for testing fxns
        BasisPtr testTrialBasis = trialOrderPtr->getBasis(testTrialID,sideIndex);
        FieldContainer<double> testTrialValuesTransformedWeighted = *(basisCache->getTransformedWeightedValues(testTrialBasis,trialOperator,
                                                                                                              sideIndex,false));
        // make copies b/c we can't fudge with return values from basisCache (const) - dimensions (Cell,Field - basis ordinal, Point)
        FieldContainer<double> testTrialValuesWeightedCopy = testTrialValuesTransformedWeighted;
        
        int numDofs2 = trialOrderPtr->getBasisCardinality(testTrialID,sideIndex); 
        for (int cellIndex=0; cellIndex<numCells; cellIndex++){
          for (int dofIndex=0; dofIndex<numDofs2; dofIndex++){
            for (int ptIndex=0; ptIndex<numPts; ptIndex++){
              testTrialValuesWeightedCopy(cellIndex, dofIndex, ptIndex) *= constrTestPair.second(cellIndex, ptIndex); // scale by constraint coeff
            }	   
          }
        }
        
        // loop thru pairs of trialIDs and constr coeffs
        for (map<int,FieldContainer<double> >::iterator constrIDIt = constrCoeffs.begin();
             constrIDIt !=constrCoeffs.end(); constrIDIt++) {
          pair<int,FieldContainer<double> > constrPair = *constrIDIt;
          int trialID = constrPair.first;
          
          // get basis to integrate
          BasisPtr trialBasis1 = trialOrderPtr->getBasis(trialID,sideIndex);
          // for trial: the value lives on the side, so we don't use the volume coords either:
          FieldContainer<double> trialValuesTransformed = *(basisCache->getTransformedValues(trialBasis1,trialOperator,sideIndex,false));
          // make copies b/c we can't fudge with return values from basisCache (const) - dimensions (Cell,Field - basis ordinal, Point)
          FieldContainer<double> trialValuesCopy = trialValuesTransformed;
          
          // transform trial values
          int numDofs1 = trialOrderPtr->getBasisCardinality(trialID,sideIndex); 
          for (int dofIndex=0; dofIndex<numDofs1; dofIndex++){
            for (int cellIndex=0; cellIndex<numCells; cellIndex++){
              for (int ptIndex=0; ptIndex<numPts; ptIndex++){
                trialValuesCopy(cellIndex, dofIndex, ptIndex) *= constrPair.second(cellIndex, ptIndex); // scale by constraint coeff
              }	   
            }
          }
          
          /////////////////////////////////////////////////////////////////////////////////////
          
          
          // integrate the transformed values, add them to the relevant trial/testTrialID dof combos
          FieldContainer<double> unweightedPenaltyMatrix(numCells,numDofs1,numDofs2);
          FunctionSpaceTools::integrate<double>(unweightedPenaltyMatrix,trialValuesCopy,testTrialValuesWeightedCopy,COMP_BLAS);
          
          for (int cellIndex=0; cellIndex<numCells; cellIndex++){
            for (int testDofIndex=0; testDofIndex<numDofs2; testDofIndex++){		
              int localTestDof = trialOrderPtr->getDofIndex(testTrialID, testDofIndex, sideIndex);
              for (int trialDofIndex=0; trialDofIndex<numDofs1; trialDofIndex++){		
                int localTrialDof = trialOrderPtr->getDofIndex(trialID, trialDofIndex, sideIndex);
                localStiffnessMatrix(cellIndex,localTrialDof,localTestDof) 
                                 += penaltyParameter*unweightedPenaltyMatrix(cellIndex,trialDofIndex,testDofIndex);
              }
            }
//.........这里部分代码省略.........

int main(int argc, char *argv[])
{
#ifdef HAVE_MPI
  Teuchos::GlobalMPISession mpiSession(&argc, &argv,0);
  int rank=mpiSession.getRank();
  int numProcs=mpiSession.getNProc();
#else
  int rank = 0;
  int numProcs = 1;
#endif
  int polyOrder = 3;
  int pToAdd = 2; // for tests

  // define our manufactured solution or problem bilinear form:
  double epsilon = 1e-2;
  bool useTriangles = false;

  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 = polyOrder + 1;
  int horizontalCells = 1, verticalCells = 1;

  double energyThreshold = 0.2; // for mesh refinements
  double nonlinearStepSize = 0.5;
  double nonlinearRelativeEnergyTolerance = 0.015; // used to determine convergence of the nonlinear solution

  ////////////////////////////////////////////////////////////////////
  // SET UP PROBLEM
  ////////////////////////////////////////////////////////////////////

  Teuchos::RCP<BurgersBilinearForm> oldBurgersBF = Teuchos::rcp(new BurgersBilinearForm(epsilon));

  // new-style bilinear form definition
  VarFactory varFactory;
  VarPtr uhat = varFactory.traceVar("\\widehat{u}");
  VarPtr beta_n_u_minus_sigma_hat = varFactory.fluxVar("\\widehat{\\beta_n u - \\sigma_n}");
  VarPtr u = varFactory.fieldVar("u");
  VarPtr sigma1 = varFactory.fieldVar("\\sigma_1");
  VarPtr sigma2 = varFactory.fieldVar("\\sigma_2");

  VarPtr tau = varFactory.testVar("\\tau",HDIV);
  VarPtr v = varFactory.testVar("v",HGRAD);
  BFPtr bf = Teuchos::rcp( new BF(varFactory) );

  // create a pointer to a new mesh:
  Teuchos::RCP<Mesh> mesh = Mesh::buildQuadMesh(quadPoints, horizontalCells, verticalCells, bf, H1Order, H1Order+pToAdd, useTriangles);
  mesh->setPartitionPolicy(Teuchos::rcp(new ZoltanMeshPartitionPolicy("HSFC")));

  Teuchos::RCP<Solution> backgroundFlow = Teuchos::rcp(new Solution(mesh, Teuchos::rcp((BC*)NULL) , Teuchos::rcp((RHS*)NULL), Teuchos::rcp((DPGInnerProduct*)NULL))); // create null solution
  oldBurgersBF->setBackgroundFlow(backgroundFlow);

  // tau parts:
  // 1/eps (sigma, tau)_K + (u, div tau)_K - (u_hat, tau_n)_dK
  bf->addTerm(sigma1 / epsilon, tau->x());
  bf->addTerm(sigma2 / epsilon, tau->y());
  bf->addTerm(u, tau->div());
  bf->addTerm( - uhat, tau->dot_normal() );

  vector<double> e1(2); // (1,0)
  e1[0] = 1;
  vector<double> e2(2); // (0,1)
  e2[1] = 1;

  FunctionPtr u_prev = Teuchos::rcp( new PreviousSolutionFunction(backgroundFlow, u) );
  FunctionPtr beta = e1 * u_prev + Teuchos::rcp( new ConstantVectorFunction( e2 ) );

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

void Projector::projectFunctionOntoBasis(FieldContainer<double> &basisCoefficients, FunctionPtr fxn, 
                                         BasisPtr basis, BasisCachePtr basisCache, IPPtr ip, VarPtr v,
                                         set<int> fieldIndicesToSkip) {
  CellTopoPtr cellTopo = basis->domainTopology();
  DofOrderingPtr dofOrderPtr = Teuchos::rcp(new DofOrdering());
  
  if (! fxn.get()) {
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "fxn cannot be null!");
  }
  
  int cardinality = basis->getCardinality();
  int numCells = basisCache->getPhysicalCubaturePoints().dimension(0);
  int numDofs = cardinality - fieldIndicesToSkip.size();
  if (numDofs==0) {
    // we're skipping all the fields, so just initialize basisCoefficients to 0 and return
    basisCoefficients.resize(numCells,cardinality);
    basisCoefficients.initialize(0);
    return;
  }
  
  FieldContainer<double> gramMatrix(numCells,cardinality,cardinality);
  FieldContainer<double> ipVector(numCells,cardinality);

  // fake a DofOrdering
  DofOrderingPtr dofOrdering = Teuchos::rcp( new DofOrdering );
  if (! basisCache->isSideCache()) {
    dofOrdering->addEntry(v->ID(), basis, v->rank());
  } else {
    dofOrdering->addEntry(v->ID(), basis, v->rank(), basisCache->getSideIndex());
  }
  
  ip->computeInnerProductMatrix(gramMatrix, dofOrdering, basisCache);
  ip->computeInnerProductVector(ipVector, v, fxn, dofOrdering, basisCache);
  
//  cout << "physical points for projection:\n" << basisCache->getPhysicalCubaturePoints();
//  cout << "gramMatrix:\n" << gramMatrix;
//  cout << "ipVector:\n" << ipVector;
  
  map<int,int> oldToNewIndices;
  if (fieldIndicesToSkip.size() > 0) {
    // the code to do with fieldIndicesToSkip might not be terribly efficient...
    // (but it's not likely to be called too frequently)
    int i_indices_skipped = 0;
    for (int i=0; i<cardinality; i++) {
      int new_index;
      if (fieldIndicesToSkip.find(i) != fieldIndicesToSkip.end()) {
        i_indices_skipped++;
        new_index = -1;
      } else {
        new_index = i - i_indices_skipped;
      }
      oldToNewIndices[i] = new_index;
    }
    
    FieldContainer<double> gramMatrixFiltered(numCells,numDofs,numDofs);
    FieldContainer<double> ipVectorFiltered(numCells,numDofs);
    // now filter out the values that we're to skip
    
    for (int cellIndex=0; cellIndex<numCells; cellIndex++) {
      for (int i=0; i<cardinality; i++) {
        int i_filtered = oldToNewIndices[i];
        if (i_filtered == -1) {
          continue;
        }
        ipVectorFiltered(cellIndex,i_filtered) = ipVector(cellIndex,i);
        
        for (int j=0; j<cardinality; j++) {
          int j_filtered = oldToNewIndices[j];
          if (j_filtered == -1) {
            continue;
          }
          gramMatrixFiltered(cellIndex,i_filtered,j_filtered) = gramMatrix(cellIndex,i,j);
        }
      }
    }
//    cout << "gramMatrixFiltered:\n" << gramMatrixFiltered;
//    cout << "ipVectorFiltered:\n" << ipVectorFiltered;
    gramMatrix = gramMatrixFiltered;
    ipVector = ipVectorFiltered;
  }
  
  for (int cellIndex=0; cellIndex<numCells; cellIndex++){
    
    // TODO: rewrite to take advantage of SerialDenseWrapper...
    Epetra_SerialDenseSolver solver;
    
    Epetra_SerialDenseMatrix A(Copy,
                               &gramMatrix(cellIndex,0,0),
                               gramMatrix.dimension(2), 
                               gramMatrix.dimension(2),  
                               gramMatrix.dimension(1)); // stride -- fc stores in row-major order (a.o.t. SDM)
    
    Epetra_SerialDenseVector b(Copy,
                               &ipVector(cellIndex,0),
                               ipVector.dimension(1));
    
    Epetra_SerialDenseVector x(gramMatrix.dimension(1));
    
    solver.SetMatrix(A);
    int info = solver.SetVectors(x,b);
//.........这里部分代码省略.........

void RieszRep::distributeDofs(){
  int myRank = Teuchos::GlobalMPISession::getRank();
  int numRanks = Teuchos::GlobalMPISession::getNProc();
#ifdef HAVE_MPI
  Epetra_MpiComm Comm(MPI_COMM_WORLD);
  //cout << "rank: " << rank << " of " << numProcs << endl;
#else
  Epetra_SerialComm Comm;
#endif  

  // the code below could stand to be reworked; I'm pretty sure this is not the best way to distribute the data, and it would also be best to get rid of the iteration over the global set of active elements.  But a similar point could be made about this method as a whole: do we really need to distribute all the dofs to every rank?  It may be best to eliminate this method altogether.
  
  vector<GlobalIndexType> cellIDsByPartitionOrdering;
  for (int rank=0; rank<numRanks; rank++) {
    set<GlobalIndexType> cellIDsForRank = _mesh->globalDofAssignment()->cellsInPartition(rank);
    cellIDsByPartitionOrdering.insert(cellIDsByPartitionOrdering.end(), cellIDsForRank.begin(), cellIDsForRank.end());
  }
  // determine inverse map:
  map<GlobalIndexType,int> ordinalForCellID;
  for (int ordinal=0; ordinal<cellIDsByPartitionOrdering.size(); ordinal++) {
    GlobalIndexType cellID = cellIDsByPartitionOrdering[ordinal];
    ordinalForCellID[cellID] = ordinal;
//    cout << "ordinalForCellID[" << cellID << "] = " << ordinal << endl;
  }
  
  for (int cellOrdinal=0; cellOrdinal<cellIDsByPartitionOrdering.size(); cellOrdinal++) {
    GlobalIndexType cellID = cellIDsByPartitionOrdering[cellOrdinal];
    ElementTypePtr elemTypePtr = _mesh->getElementType(cellID);
    DofOrderingPtr testOrderingPtr = elemTypePtr->testOrderPtr;
    int numDofs = testOrderingPtr->totalDofs();
    
    int cellIDPartition = _mesh->partitionForCellID(cellID);
    bool isInPartition = (cellIDPartition == myRank);

    int numMyDofs;
    FieldContainer<double> dofs(numDofs);
    if (isInPartition){  // if in partition
      numMyDofs = numDofs;
      dofs = _rieszRepDofs[cellID];
    } else{
      numMyDofs = 0;
    }
    
    Epetra_Map dofMap(numDofs,numMyDofs,0,Comm); 
    Epetra_Vector distributedRieszDofs(dofMap);
    if (isInPartition) {
      for (int i = 0;i<numMyDofs;i++) { // shouldn't activate on off-proc partitions
        distributedRieszDofs.ReplaceGlobalValues(1,&dofs(i),&i);
      }     
    }
    Epetra_Map importMap(numDofs,numDofs,0,Comm); // every proc should own their own copy of the dofs
    Epetra_Import testDofImporter(importMap, dofMap); 
    Epetra_Vector globalRieszDofs(importMap);
    globalRieszDofs.Import(distributedRieszDofs, testDofImporter, Insert);  
    if (!isInPartition){
      for (int i = 0;i<numDofs;i++){
        dofs(i) = globalRieszDofs[i];
      }      
    }
    _rieszRepDofsGlobal[cellID] = dofs;
//    { // debugging
//      ostringstream cellIDlabel;
//      cellIDlabel << "cell " << cellID << " _rieszRepDofsGlobal, after global import";
//      TestSuite::serializeOutput(cellIDlabel.str(), _rieszRepDofsGlobal[cellID]);
//    }
  }
  
  // distribute norms as well
  GlobalIndexType numElems = _mesh->numActiveElements();
  set<GlobalIndexType> rankLocalCellIDs = _mesh->cellIDsInPartition();
  IndexType numMyElems = rankLocalCellIDs.size();
  GlobalIndexType myElems[numMyElems];
  // build cell index
  GlobalIndexType myCellOrdinal = 0;

  double rankLocalRieszNorms[numMyElems];
  
  for (set<GlobalIndexType>::iterator cellIDIt = rankLocalCellIDs.begin(); cellIDIt != rankLocalCellIDs.end(); cellIDIt++) {
    GlobalIndexType cellID = *cellIDIt;
    myElems[myCellOrdinal] = ordinalForCellID[cellID];
    rankLocalRieszNorms[myCellOrdinal] = _rieszRepNormSquared[cellID];
    myCellOrdinal++;
  }
  Epetra_Map normMap((GlobalIndexTypeToCast)numElems,(int)numMyElems,(GlobalIndexTypeToCast *)myElems,(GlobalIndexTypeToCast)0,Comm);

  Epetra_Vector distributedRieszNorms(normMap);
  int err = distributedRieszNorms.ReplaceGlobalValues(numMyElems,rankLocalRieszNorms,(GlobalIndexTypeToCast *)myElems);
  if (err != 0) {
    cout << "RieszRep::distributeDofs(): on rank" << myRank << ", ReplaceGlobalValues returned error code " << err << endl;
  }

  Epetra_Map normImportMap((GlobalIndexTypeToCast)numElems,(GlobalIndexTypeToCast)numElems,0,Comm);
  Epetra_Import normImporter(normImportMap,normMap); 
  Epetra_Vector globalNorms(normImportMap);
  globalNorms.Import(distributedRieszNorms, normImporter, Add);  // add should be OK (everything should be zeros)

  for (int cellOrdinal=0; cellOrdinal<cellIDsByPartitionOrdering.size(); cellOrdinal++) {
    GlobalIndexType cellID = cellIDsByPartitionOrdering[cellOrdinal];
    _rieszRepNormSquaredGlobal[cellID] = globalNorms[cellOrdinal];
//    if (myRank==0) cout << "_rieszRepNormSquaredGlobal[" << cellID << "] = " << globalNorms[cellOrdinal] << endl;
//.........这里部分代码省略.........

bool MeshTestUtility::neighborBasesAgreeOnSides(Teuchos::RCP<Mesh> mesh, Epetra_MultiVector &globalSolutionCoefficients)
{
  bool success = true;
  MeshTopologyViewPtr meshTopo = mesh->getTopology();
  int spaceDim = meshTopo->getDimension();

  set<IndexType> activeCellIndices = meshTopo->getActiveCellIndices();
  for (set<IndexType>::iterator cellIt=activeCellIndices.begin(); cellIt != activeCellIndices.end(); cellIt++)
  {
    IndexType cellIndex = *cellIt;
    CellPtr cell = meshTopo->getCell(cellIndex);

    BasisCachePtr fineCellBasisCache = BasisCache::basisCacheForCell(mesh, cellIndex);
    ElementTypePtr fineElemType = mesh->getElementType(cellIndex);
    DofOrderingPtr fineElemTrialOrder = fineElemType->trialOrderPtr;

    FieldContainer<double> fineSolutionCoefficients(fineElemTrialOrder->totalDofs());
    mesh->globalDofAssignment()->interpretGlobalCoefficients(cellIndex, fineSolutionCoefficients, globalSolutionCoefficients);
//    if ((cellIndex==0) || (cellIndex==2)) {
//      cout << "MeshTestUtility: local coefficients for cell " << cellIndex << ":\n" << fineSolutionCoefficients;
//    }

    unsigned sideCount = cell->getSideCount();
    for (unsigned sideOrdinal=0; sideOrdinal<sideCount; sideOrdinal++)
    {
      FieldContainer<double> fineSideRefPoints, fineCellRefPoints, coarseSideRefPoints, coarseCellRefPoints;
      bool hasCoarserNeighbor = determineRefTestPointsForNeighbors(meshTopo, cell, sideOrdinal, fineSideRefPoints, fineCellRefPoints, coarseSideRefPoints, coarseCellRefPoints);
      if (!hasCoarserNeighbor) continue;

      pair<GlobalIndexType, unsigned> neighborInfo = cell->getNeighborInfo(sideOrdinal, meshTopo);

      CellPtr neighborCell = meshTopo->getCell(neighborInfo.first);

      unsigned numTestPoints = coarseCellRefPoints.dimension(0);

      //        cout << "testing neighbor agreement between cell " << cellIndex << " and " << neighborCell->cellIndex() << endl;

      // if we get here, the cell has a neighbor on this side, and is at least as fine as that neighbor.

      BasisCachePtr fineSideBasisCache = fineCellBasisCache->getSideBasisCache(sideOrdinal);

      fineCellBasisCache->setRefCellPoints(fineCellRefPoints);

      BasisCachePtr coarseCellBasisCache = BasisCache::basisCacheForCell(mesh, neighborInfo.first);
      BasisCachePtr coarseSideBasisCache = coarseCellBasisCache->getSideBasisCache(neighborInfo.second);

      coarseCellBasisCache->setRefCellPoints(coarseCellRefPoints);

      bool hasSidePoints = (fineSideRefPoints.size() > 0);
      if (hasSidePoints)
      {
        fineSideBasisCache->setRefCellPoints(fineSideRefPoints);
        coarseSideBasisCache->setRefCellPoints(coarseSideRefPoints);
      }

      // sanity check: do physical points match?

      FieldContainer<double> fineCellPhysicalCubaturePoints = fineCellBasisCache->getPhysicalCubaturePoints();
      FieldContainer<double> coarseCellPhysicalCubaturePoints = coarseCellBasisCache->getPhysicalCubaturePoints();

      double tol = 1e-14;
      double maxDiff = 0;
      if (! fcsAgree(coarseCellPhysicalCubaturePoints, fineCellPhysicalCubaturePoints, tol, maxDiff) )
      {
        cout << "ERROR: MeshTestUtility::neighborBasesAgreeOnSides internal error: fine and coarse cell cubature points do not agree.\n";
        success = false;
        continue;
      }

      if (hasSidePoints)
      {
        FieldContainer<double> fineSidePhysicalCubaturePoints = fineSideBasisCache->getPhysicalCubaturePoints();
        FieldContainer<double> coarseSidePhysicalCubaturePoints = coarseSideBasisCache->getPhysicalCubaturePoints();

        double tol = 1e-14;
        double maxDiff = 0;
        if (! fcsAgree(fineSidePhysicalCubaturePoints, fineCellPhysicalCubaturePoints, tol, maxDiff) )
        {
          cout << "ERROR: MeshTestUtility::neighborBasesAgreeOnSides internal error: fine side and cell cubature points do not agree.\n";
          success = false;
          continue;
        }
        if (! fcsAgree(coarseSidePhysicalCubaturePoints, coarseCellPhysicalCubaturePoints, tol, maxDiff) )
        {
          cout << "ERROR: MeshTestUtility::neighborBasesAgreeOnSides internal error: coarse side and cell cubature points do not agree.\n";
          success = false;
          continue;
        }
        if (! fcsAgree(coarseSidePhysicalCubaturePoints, fineSidePhysicalCubaturePoints, tol, maxDiff) )
        {
          cout << "ERROR: MeshTestUtility::neighborBasesAgreeOnSides internal error: fine and coarse side cubature points do not agree.\n";
          success = false;
          continue;
        }
      }

      ElementTypePtr coarseElementType = mesh->getElementType(neighborInfo.first);
      DofOrderingPtr coarseElemTrialOrder = coarseElementType->trialOrderPtr;

      FieldContainer<double> coarseSolutionCoefficients(coarseElemTrialOrder->totalDofs());
//.........这里部分代码省略.........

bool MeshTestUtility::checkLocalGlobalConsistency(MeshPtr mesh, double tol)
{
  bool success = true;

  set<GlobalIndexType> cellIDs = mesh->getActiveCellIDs();

  GlobalDofAssignmentPtr gda = mesh->globalDofAssignment();

  // TODO: make this only check the locally-owned cells (right now does the whole mesh on every rank)

  int numGlobalDofs = gda->globalDofCount();
  FieldContainer<double> globalCoefficients(numGlobalDofs);
  for (int i=0; i<numGlobalDofs; i++)
  {
    globalCoefficients(i) = 2*i + 1; // arbitrary cofficients
  }
  FieldContainer<double> globalCoefficientsExpected = globalCoefficients;
  FieldContainer<double> globalCoefficientsActual(numGlobalDofs);
  FieldContainer<double> localCoefficients;

  Epetra_SerialComm Comm;
  Epetra_BlockMap map(numGlobalDofs, 1, 0, Comm);
  Epetra_Vector globalCoefficientsVector(map);
  for (int i=0; i<numGlobalDofs; i++)
  {
    globalCoefficientsVector[i] = globalCoefficients(i);
  }

  cellIDs = mesh->getActiveCellIDs();
  for (set<GlobalIndexType>::iterator cellIDIt = cellIDs.begin(); cellIDIt != cellIDs.end(); cellIDIt++)
  {
    GlobalIndexType cellID = *cellIDIt;
    gda->interpretGlobalCoefficients(cellID, localCoefficients, globalCoefficientsVector);
    FieldContainer<GlobalIndexType> globalDofIndices;
    FieldContainer<double> globalCoefficientsForCell;

    DofOrderingPtr trialOrder = mesh->getElementType(cellID)->trialOrderPtr;
    set<int> varIDs = trialOrder->getVarIDs();

    for (set<int>::iterator varIt=varIDs.begin(); varIt != varIDs.end(); varIt++)
    {
      int varID = *varIt;
      const vector<int>* sidesForVar = &trialOrder->getSidesForVarID(varID);
      for (vector<int>::const_iterator sideIt = sidesForVar->begin(); sideIt != sidesForVar->end(); sideIt++)
      {
        int sideOrdinal = *sideIt;

        BasisPtr basis;
        if (sidesForVar->size() == 1)
        {
          basis = trialOrder->getBasis(varID);
        }
        else
        {
          basis = trialOrder->getBasis(varID,sideOrdinal);
        }

        FieldContainer<double> basisCoefficients(basis->getCardinality());

        for (int dofOrdinal=0; dofOrdinal<basis->getCardinality(); dofOrdinal++)
        {
          int localDofIndex;
          if (sidesForVar->size() == 1)
          {
            localDofIndex = trialOrder->getDofIndex(varID, dofOrdinal);
          }
          else
          {
            localDofIndex = trialOrder->getDofIndex(varID, dofOrdinal, sideOrdinal);
          }

          basisCoefficients(dofOrdinal) = localCoefficients(localDofIndex);
        }

        gda->interpretLocalBasisCoefficients(cellID, varID, sideOrdinal,
                                             basisCoefficients, globalCoefficientsForCell, globalDofIndices);

        //        if ( (cellID==2) && (sideOrdinal==1) && (varID==0) ) {
        //          cout << "DEBUGGING: (cellID==2) && (sideOrdinal==1).\n";
        //          cout << "globalDofIndices:\n" << globalDofIndices;
        //          cout << "globalCoefficientsForCell:\n" << globalCoefficientsForCell;
        //          cout << "basisCoefficients:\n" << basisCoefficients;
        //        }

        for (int dofOrdinal=0; dofOrdinal < globalDofIndices.size(); dofOrdinal++)
        {
          GlobalIndexType dofIndex = globalDofIndices(dofOrdinal);
          globalCoefficientsActual(dofIndex) = globalCoefficientsForCell(dofOrdinal);

          double diff = abs(globalCoefficientsForCell(dofOrdinal) - globalCoefficientsExpected(dofIndex)) / globalCoefficientsExpected(dofIndex);
          if (diff > tol)
          {
            cout << "In mapping for cell " << cellID << " and var " << varID << " on side " << sideOrdinal << ", ";
            cout << "expected coefficient for global dof index " << dofIndex << " to be " << globalCoefficientsExpected(dofIndex);
            cout << ", but was " << globalCoefficientsForCell(dofOrdinal);
            cout << " (relative diff = " << diff << "; tol = " << tol << ")\n";
            success = false;
          }
        }
      }
//.........这里部分代码省略.........

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