C++ BasisPtr类代码示例

bool basisSumEqualsFunction(FieldContainer<double> &basisCoefficients, BasisPtr basis, FunctionPtr f)
{
  // tests on [0,1]
  FunctionPtr sumFunction = Teuchos::rcp( new BasisSumFunction(basis, basisCoefficients) );

  int cubatureDegree = basis->getDegree() * 2;
  BasisCachePtr basisCache = BasisCache::basisCache1D(0, 1, cubatureDegree);

  double tol = 1e-13;
  return sumFunction->equals(f, basisCache, tol);
}

bool checkVertexOrdinalsQuad(BasisPtr basis, vector<int> &vertexOrdinals) {
  // check that the given indices are exactly the vertex basis functions:
  // a) these are (1,0) or (0,1) at the corresponding vertices
  // b) others are (0,0) at the vertices
  
  int numVertices = 4;
  
  FieldContainer<double> refCellPoints(numVertices,2); // vertices, in order
  refCellPoints(0,0) = -1;
  refCellPoints(0,1) = -1;
  refCellPoints(1,0) =  1;
  refCellPoints(1,1) = -1;
  refCellPoints(2,0) =  1;
  refCellPoints(2,1) =  1;
  refCellPoints(3,0) = -1;
  refCellPoints(3,1) =  1;
  
  // assume scalar basis for now -- we'll throw an exception if not...
  FieldContainer<double> values(basis->getCardinality(), numVertices); // F, P
  basis->getValues(values, refCellPoints, OPERATOR_VALUE);
  
  double tol = 1e-14;
  for (int vertexIndex=0; vertexIndex<numVertices; vertexIndex++) {
    int vertexOrdinal = vertexOrdinals[vertexIndex];
    for (int fieldIndex=0; fieldIndex<basis->getCardinality(); fieldIndex++) {
      double value = values(fieldIndex,vertexIndex);
      if (fieldIndex==vertexOrdinal) {
        // expect non-zero
        if (value < tol) {
          return false;
        }
      } else {
        // expect zero
        if (value > tol) {
          return false;
        }
      }
    }
  }
  return true;
}

bool VectorizedBasisTestSuite::testVectorizedBasisTags()
{
  bool success = true;

  shards::CellTopology quad_4(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
  int numVertices = quad_4.getVertexCount();
  int numComponents = quad_4.getDimension();
  int vertexDim = 0;

  for (int polyOrder = 1; polyOrder<10; polyOrder++)
  {
    BasisPtr hgradBasis =  BasisFactory::basisFactory()->getConformingBasis(polyOrder,
                           quad_4.getKey(),
                           Camellia::FUNCTION_SPACE_HGRAD);
    BasisPtr vectorHGradBasis = BasisFactory::basisFactory()->getConformingBasis( polyOrder,
                                quad_4.getKey(),
                                Camellia::FUNCTION_SPACE_VECTOR_HGRAD);
    vector<int> vertexNodeFieldIndices;
    for (int vertexIndex=0; vertexIndex<numVertices; vertexIndex++)
    {
      for (int comp=0; comp<numComponents; comp++)
      {
        int vertexNodeFieldIndex = vectorHGradBasis->getDofOrdinal(vertexDim, vertexIndex, comp);
        vertexNodeFieldIndices.push_back(vertexNodeFieldIndex);
//        cout << "vertexNodeFieldIndex for vertex index " << vertexIndex << ", comp " << comp;
//        cout << " = " << vertexNodeFieldIndex << endl;
      }
    }
    if (!checkVertexNodalIndicesQuad(vectorHGradBasis, vertexNodeFieldIndices) )
    {
      success = false;
      cout << "testVectorizedBasisTags: Vertex tags for vectorized HGRAD basis";
      cout << " of order " << polyOrder << " are incorrect.\n";
    }
  }

  return success;
}

FCPtr BasisEvaluation::getTransformedValues(BasisPtr basis, Camellia::EOperator op,
                                            const FieldContainer<double> &referencePoints,
                                            int numCells,
                                            BasisCache* basisCache)
{
  Camellia::EFunctionSpace fs = basis->functionSpace();
  int spaceDim = referencePoints.dimension(1);
  int componentOfInterest;
  Intrepid::EOperator relatedOp;
  relatedOp = relatedOperator(op, fs, spaceDim, componentOfInterest);
  
  FCPtr referenceValues = getValues(basis,(Camellia::EOperator) relatedOp, referencePoints);
  return getTransformedValuesWithBasisValues(basis,op,referenceValues,numCells,basisCache);
}

bool basisSumInterpolatesCurveEndPoints(FieldContainer<double> &basisCoefficients_x, FieldContainer<double> &basisCoefficients_y,
                                        BasisPtr basis, ParametricCurvePtr curve)
{
  double curve_x0, curve_y0, curve_x1, curve_y1;
  curve->value(0, curve_x0, curve_y0);
  curve->value(1, curve_x1, curve_y1);
  BasisCachePtr basisCache = BasisCache::basisCache1D(0, 1, basis->getDegree()*2);
  double sum_x0 = basisSumAtParametricPoint(basisCoefficients_x, basis, 0, basisCache);
  double sum_x1 = basisSumAtParametricPoint(basisCoefficients_x, basis, 1, basisCache);
  double sum_y0 = basisSumAtParametricPoint(basisCoefficients_y, basis, 0, basisCache);
  double sum_y1 = basisSumAtParametricPoint(basisCoefficients_y, basis, 1, basisCache);
  double tol = 1e-13;
  double x0_err = abs(sum_x0-curve_x0);
  double y0_err = abs(sum_y0-curve_y0);
  double x1_err = abs(sum_x1-curve_x1);
  double y1_err = abs(sum_y1-curve_y1);
  double sum_err = x0_err + y0_err + x1_err + y1_err;
  return sum_err < tol;
}

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

void Projector::projectFunctionOntoBasis(FieldContainer<double> &basisCoefficients, Teuchos::RCP<AbstractFunction> fxn, BasisPtr basis,
                                         const FieldContainer<double> &physicalCellNodes) {

  CellTopoPtr cellTopo = basis->domainTopology();
  DofOrderingPtr dofOrderPtr = Teuchos::rcp(new DofOrdering());

  int basisRank = BasisFactory::basisFactory()->getBasisRank(basis);
  int ID = 0; // only one entry for this fake dofOrderPtr
  dofOrderPtr->addEntry(ID,basis,basisRank);
  int maxTrialDegree = dofOrderPtr->maxBasisDegree();

  // do not build side caches - no projections for sides supported at the moment
  if (cellTopo->getTensorialDegree() != 0) {
    cout << "Projector::projectFunctionOntoBasis() does not yet support tensorial degree > 0.\n";
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Projector::projectFunctionOntoBasis() does not yet support tensorial degree > 0.");
  }
  shards::CellTopology shardsTopo = cellTopo->getShardsTopology();
  BasisCache basisCache(physicalCellNodes, shardsTopo, *(dofOrderPtr), maxTrialDegree, false);
  // assume only L2 projections
  IntrepidExtendedTypes::EOperatorExtended op =  IntrepidExtendedTypes::OP_VALUE;

  // have information, build inner product matrix
  int numDofs = basis->getCardinality();
  FieldContainer<double> cubPoints = basisCache.getPhysicalCubaturePoints();    

  FieldContainer<double> basisValues = *(basisCache.getTransformedValues(basis, op));
  FieldContainer<double> testBasisValues = *(basisCache.getTransformedWeightedValues(basis, op));

  int numCells = physicalCellNodes.dimension(0);
  int numPts = cubPoints.dimension(1);
  FieldContainer<double> functionValues;
  fxn->getValues(functionValues, cubPoints);  

  FieldContainer<double> gramMatrix(numCells,numDofs,numDofs);
  FieldContainer<double> ipVector(numCells,numDofs);
  FunctionSpaceTools::integrate<double>(gramMatrix,basisValues,testBasisValues,COMP_BLAS);
  FunctionSpaceTools::integrate<double>(ipVector,functionValues,testBasisValues,COMP_BLAS); 
  
  basisCoefficients.resize(numCells,numDofs);
  for (int cellIndex=0; cellIndex<numCells; cellIndex++){

    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));
    
    /*
    cout << "matrix A = " << endl;
    for (int i=0;i<gramMatrix.dimension(2);i++){
      for (int j=0;j<gramMatrix.dimension(1);j++){
	cout << A(i,j) << " ";
      }
      cout << endl;
    }
    cout << endl;

    cout << "vector B = " << endl;
    for (int i=0;i<functionValues.dimension(1);i++){
      cout << b(i) << endl;
    }
    */

    solver.SetMatrix(A);
    int info = solver.SetVectors(x,b);
    if (info!=0){
      cout << "projectFunctionOntoBasis: failed to SetVectors with error " << info << endl;
    }

    bool equilibrated = false;
    if (solver.ShouldEquilibrate()){
      solver.EquilibrateMatrix();
      solver.EquilibrateRHS();      
      equilibrated = true;
    }   

    info = solver.Solve();
    if (info!=0){
      cout << "projectFunctionOntoBasis: failed to solve with error " << info << endl;
    }

    if (equilibrated) {
      int successLocal = solver.UnequilibrateLHS();
      if (successLocal != 0) {
        cout << "projection: unequilibration FAILED with error: " << successLocal << endl;
      }
    }

    for (int i=0;i<numDofs;i++){
      basisCoefficients(cellIndex,i) = x(i);
    }   
    
//.........这里部分代码省略.........

void TRieszRep<Scalar>::computeRepresentationValues(FieldContainer<Scalar> &values, int testID, Camellia::EOperator op, BasisCachePtr basisCache)
{

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

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

  // all elems coming in should be of same type
  ElementTypePtr elemTypePtr = _mesh->getElementType(cellIDs[0]);
  DofOrderingPtr testOrderingPtr = elemTypePtr->testOrderPtr;
  int numTestDofsForVarID = testOrderingPtr->getBasisCardinality(testID, VOLUME_INTERIOR_SIDE_ORDINAL);
  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 testBasisClassifications(BasisPtr basis) {
  bool success = true;
  
  CellTopoPtr cellTopo = basis->domainTopology();
  
  int numVertices = cellTopo->getVertexCount();
  int numEdges = cellTopo->getEdgeCount();
  
  int degree = basis->getDegree();
  
  // TODO: finish this
  vector<int> vertexOrdinals;
  for (int vertexIndex=0; vertexIndex < numVertices; vertexIndex++) {
    set<int> dofOrdinals = basis->dofOrdinalsForVertex(vertexIndex);
    if (dofOrdinals.size() == 0) TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "No dofOrdinal for vertex...");
    if (dofOrdinals.size() > 1) TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "More than one dofOrdinal per vertex...");
    vertexOrdinals.push_back(*(dofOrdinals.begin()));
  }
//  
//  if (! checkVertexOrdinalsQuad(basis, vertexOrdinals) ) {
//    success = false;
//    cout << "vertex dof ordinals don't match expected\n";
//    cout << "ordinals: ";
//    for (int vertexIndex=0; vertexIndex < numVertices; vertexIndex++) {
//      cout << vertexOrdinals[vertexIndex] << " ";
//    }
//    cout << endl;
//  }
  
  // get the points in reference space for each vertex
  FieldContainer<double> points;
  if (numVertices == 2) { // line
    points.resize(2,1);
    points(0,0) = -1;
    points(1,0) = 1;
  } else if (numVertices == 3) { // triangle
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "triangles not yet supported");
  } else if (numVertices == 4) { // quad
    points.resize(4,2);
    points(0,0) = -1;
    points(0,1) = -1;
    points(1,0) =  1;
    points(1,1) = -1;
    points(2,0) =  1;
    points(2,1) =  1;
    points(3,0) = -1;
    points(3,1) =  1;
  } else {
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "unsupported topology");
  }
  
  FieldContainer<double> vertexValues;
  if (basis->rangeRank() > 0) {
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "rank > 0 bases not yet supported");
  } else {
    vertexValues.resize(basis->getCardinality(),numVertices);
  }
  
  basis->getValues(vertexValues, points, Intrepid::OPERATOR_VALUE);
  
  // test that the points are correctly classified
  for (int fieldIndex=0; fieldIndex<basis->getCardinality(); fieldIndex++) {
    for (int ptIndex=0; ptIndex<numVertices; ptIndex++) {
      int dofOrdinalForPoint = vertexOrdinals[ptIndex];
      bool expectZero = (dofOrdinalForPoint != fieldIndex);
      if (expectZero) {
        if (vertexValues(fieldIndex,ptIndex) != 0) {
          success = false;
          cout << "Expected 0 for fieldIndex " << fieldIndex << " and ptIndex " << ptIndex;
          cout << ", but got " << vertexValues(fieldIndex,ptIndex) << endl;
        }
      } else {
        if (vertexValues(fieldIndex,ptIndex) == 0) {
          cout << "Expected nonzero for fieldIndex " << fieldIndex << " and ptIndex " << ptIndex << endl;
          success = false;
        }
      }
    }
  }
  
  if (!success) {
    cout << "Failed testBasisClassifications; vertexValues:\n" << vertexValues;
  }
  
  return success;
}

FCPtr BasisEvaluation::getValues(BasisPtr basis, Camellia::EOperator op,
                                 const FieldContainer<double> &refPoints)
{
  int numPoints = refPoints.dimension(0);
  int spaceDim = refPoints.dimension(1);  // points dimensions are (numPoints, spaceDim)
  int basisCardinality = basis->getCardinality();
  int componentOfInterest = -1;
  Camellia::EFunctionSpace fs = basis->functionSpace();
  Intrepid::EOperator relatedOp = relatedOperator(op, fs, spaceDim, componentOfInterest);

  int opCardinality = spaceDim; // for now, we assume basis values are in the same spaceDim as points (e.g. vector 1D has just 1 component)
  
  bool relatedOpIsDkOperator = (relatedOp >= OPERATOR_D1) && (relatedOp <= OPERATOR_D10);
  if (relatedOpIsDkOperator)
  {
    opCardinality = Intrepid::getDkCardinality(relatedOp, spaceDim);
  }
  
  if ((Camellia::EOperator)relatedOp != op)
  {
    // we can assume relatedResults has dimensions (numPoints,basisCardinality,spaceDimOut)
    FCPtr relatedResults = Teuchos::rcp(new FieldContainer<double>(basisCardinality,numPoints,opCardinality));
    basis->getValues(*(relatedResults.get()), refPoints, relatedOp);
    FCPtr result = getComponentOfInterest(relatedResults,op,fs,componentOfInterest);
    if ( result.get() == 0 )
    {
      result = relatedResults;
    }
    return result;
  }
  // if we get here, we should have a standard Intrepid operator, in which case we should
  // be able to: size a FieldContainer appropriately, and then call basis->getValues

  // But let's do just check that we have a standard Intrepid operator
  if ( (op >= Camellia::OP_X) || (op <  Camellia::OP_VALUE) )
  {
    TEUCHOS_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Unknown operator.");
  }

  // result dimensions should be either (numPoints,basisCardinality) or (numPoints,basisCardinality,spaceDimOut);
  Teuchos::Array<int> dimensions;
  dimensions.push_back(basisCardinality);
  dimensions.push_back(numPoints);
  int basisRank = basis->rangeRank();
  if (((basisRank == 0) && (op ==  Camellia::OP_VALUE)))
  {
    // scalar; leave as is
  }
  else if ((basisRank == 0) && relatedOpIsDkOperator)
  {
    dimensions.push_back(opCardinality);
  }
  else if ( ( ( basisRank == 1) && (op ==  Camellia::OP_VALUE) ) )
  {
    dimensions.push_back(opCardinality);
  }
  else if (
    ( ( basisRank == 0) && (op == Camellia::OP_GRAD) )
    ||
    ( ( basisRank == 0) && (op == Camellia::OP_CURL) ) )
  {
    dimensions.push_back(spaceDim);
  }
  else if ( (basis->rangeRank() == 1) && (op == Camellia::OP_GRAD) )
  {
    // grad of vector: a tensor
    dimensions.push_back(spaceDim);
    dimensions.push_back(opCardinality);
  }
  FCPtr result = Teuchos::rcp(new FieldContainer<double>(dimensions));
  basis->getValues(*(result.get()), refPoints, (Intrepid::EOperator)op);
  return result;
}

int DofOrdering::getBasisCardinality(int varID, int sideIndex) {
  BasisPtr basis = getBasis(varID,sideIndex);
  return basis->getCardinality();
//  return getBasis(varID,sideIndex)->getCardinality();
}

FCPtr BasisEvaluation::getTransformedValuesWithBasisValues(BasisPtr basis, Camellia::EOperator op, int spaceDim,
                                                           constFCPtr referenceValues, int numCells,
                                                           const FieldContainer<double> &cellJacobian,
                                                           const FieldContainer<double> &cellJacobianInv,
                                                           const FieldContainer<double> &cellJacobianDet)
{
  typedef FunctionSpaceTools fst;
//  int numCells = cellJacobian.dimension(0);
  
//  int spaceDim = basis->domainTopology()->getDimension(); // changed 2/18/15
//  // 6-16-14 NVR: getting the spaceDim from cellJacobian's dimensioning is the way we've historically done it.
//  // I think it might be better to do this using basis->domainTopology() generally, but for now we only make the
//  // switch in case the domain topology is a Node.
//  if (basis->domainTopology()->getDimension() == 0) {
//    spaceDim = 0;
//  } else {
//    spaceDim = cellJacobian.dimension(2);
//  }

  int componentOfInterest;
  Camellia::EFunctionSpace fs = basis->functionSpace();
  Intrepid::EOperator relatedOp = relatedOperator(op,fs,spaceDim, componentOfInterest);
  Teuchos::Array<int> dimensions;
  referenceValues->dimensions(dimensions);
  dimensions.insert(dimensions.begin(), numCells);
  Teuchos::RCP<FieldContainer<double> > transformedValues = Teuchos::rcp(new FieldContainer<double>(dimensions));
  bool vectorizedBasis = functionSpaceIsVectorized(fs);
  if (vectorizedBasis && (op ==  Camellia::OP_VALUE))
  {
    TEUCHOS_TEST_FOR_EXCEPTION( vectorizedBasis && (op ==  Camellia::OP_VALUE),
                                std::invalid_argument, "Vector HGRAD/HVOL with OP_VALUE not supported by getTransformedValuesWithBasisValues.  Please use getTransformedVectorValuesWithComponentBasisValues instead.");
  }
  switch(relatedOp)
  {
  case(Intrepid::OPERATOR_VALUE):
    switch(fs)
    {
    case Camellia::FUNCTION_SPACE_REAL_SCALAR:
    //          cout << "Reference values for FUNCTION_SPACE_REAL_SCALAR: " << *referenceValues;
    case Camellia::FUNCTION_SPACE_HGRAD:
    case Camellia::FUNCTION_SPACE_HGRAD_DISC:
      fst::HGRADtransformVALUE<double>(*transformedValues,*referenceValues);
      break;
    case Camellia::FUNCTION_SPACE_HCURL:
    case Camellia::FUNCTION_SPACE_HCURL_DISC:
      fst::HCURLtransformVALUE<double>(*transformedValues,cellJacobianInv,*referenceValues);
      break;
    case Camellia::FUNCTION_SPACE_HDIV:
    case Camellia::FUNCTION_SPACE_HDIV_DISC:
    case Camellia::FUNCTION_SPACE_HDIV_FREE:
      fst::HDIVtransformVALUE<double>(*transformedValues,cellJacobian,cellJacobianDet,*referenceValues);
        break;
      case Camellia::FUNCTION_SPACE_HVOL:
      case Camellia::FUNCTION_SPACE_HVOL_SPACE_HGRAD_TIME:
      case Camellia::FUNCTION_SPACE_HGRAD_SPACE_HVOL_TIME:
//        {
//          static bool haveWarned = false;
//          if (!haveWarned) {
//            cout << "WARNING: for the moment, switching to the standard HVOLtransformVALUE method.\n";
//            haveWarned = true;
//          }
//        }
//          fst::HVOLtransformVALUE<double>(*transformedValues, cellJacobianDet, *referenceValues);
      // for the moment, use the fact that we know the HVOL basis is always an HGRAD basis:
      // (I think using the below amounts to solving for the HVOL variables scaled by Jacobian)
      fst::HGRADtransformVALUE<double>(*transformedValues,*referenceValues);
      break;
    default:
      TEUCHOS_TEST_FOR_EXCEPTION(true,std::invalid_argument, "unhandled transformation");
      break;
    }
    break;
  case(Intrepid::OPERATOR_GRAD):
    switch(fs)
    {
    case Camellia::FUNCTION_SPACE_HVOL:
    case Camellia::FUNCTION_SPACE_HGRAD:
    case Camellia::FUNCTION_SPACE_HGRAD_DISC:
      fst::HGRADtransformGRAD<double>(*transformedValues,cellJacobianInv,*referenceValues);
      break;
    case Camellia::FUNCTION_SPACE_VECTOR_HVOL:
    case Camellia::FUNCTION_SPACE_VECTOR_HGRAD:
    case Camellia::FUNCTION_SPACE_VECTOR_HGRAD_DISC:
      // referenceValues has dimensions (F,P,D1,D2).  D1 is our component dimension, and D2 is the one that came from the gradient.
      // HGRADtransformGRAD expects (F,P,D) for input, and (C,F,P,D) for output.
      // If we split referenceValues into (F,P,D1=0,D2) and (F,P,D1=1,D2), then we can transform each of those, and then interleave the results…
    {
      // block off so we can create new stuff inside the switch case
      Teuchos::Array<int> dimensions;
      referenceValues->dimensions(dimensions);
      int numFields = dimensions[0];
      int numPoints = dimensions[1];
      int D1 = dimensions[dimensions.size()-2];
      int D2 = dimensions[dimensions.size()-1];
      dimensions[dimensions.size()-2] = D2; // put D2 in the D1 spot
      dimensions.pop_back(); // get rid of original D2
      FieldContainer<double> refValuesSlice(dimensions);
      dimensions.insert(dimensions.begin(),numCells);
      FieldContainer<double> transformedValuesSlice(dimensions);

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

bool ParametricCurveTests::testCircularArc()
{
  bool success = true;

  // the arc details are copied from CurvilinearMeshTests -- motivation is to diagnose test failure there with a more granular test here
  double radius = 1.0;
  double meshWidth = sqrt(2);

  ParametricCurvePtr circle = ParametricCurve::circle(radius, meshWidth / 2.0, meshWidth / 2.0);
  ParametricCurvePtr circularArc = ParametricCurve::subCurve(circle,  5.0/8.0, 7.0/8.0);

  BasisCachePtr basisCache = BasisCache::parametric1DCache(15); // overintegrate to be safe

  FunctionPtr cos_part = Teuchos::rcp( new Cos_ax(PI/2, 1.25*PI));
  FunctionPtr sin_part = Teuchos::rcp( new Sin_ax(PI/2, 1.25*PI));
  FunctionPtr x_t = meshWidth / 2 + cos_part;
  FunctionPtr y_t = meshWidth / 2 + sin_part;

  FunctionPtr dx_dt = (- PI / 2) * sin_part;
  FunctionPtr dy_dt = (PI / 2) * cos_part;

  double arcIntegral_x = circularArc->x()->integrate(basisCache);
  double x_t_integral = x_t->integrate(basisCache);

  // check that we have the right idea for all those functions:
  if (! x_t->equals(circularArc->x(),basisCache))
  {
    double x1_expected = Function::evaluate(x_t,1);
    double x1_actual;
    circularArc->xPart()->value(1,x1_actual);
    cout << "expected x(0) = " << x1_expected;
    cout << "; actual = " << x1_actual << endl;
    cout << "x part of circularArc doesn't match expected.\n";
    success = false;
  }
  if (! y_t->equals(circularArc->y(),basisCache))
  {
    double y1_actual;
    circularArc->yPart()->value(1,y1_actual);
    cout << "expected y(1) = " << Function::evaluate(y_t,1);
    cout << "; actual = " << y1_actual << endl;
    cout << "y part of circularArc doesn't match expected.\n";
    success = false;
  }

  if (! dx_dt->equals(circularArc->dt_parametric()->x(),basisCache))
  {
    cout << "dx/dt of circularArc doesn't match expected.\n";
    success = false;
  }
  if (! dy_dt->equals(circularArc->dt_parametric()->y(),basisCache))
  {
    cout << "dy/dt of circularArc doesn't match expected.\n";
    success = false;
  }

  // test exact curve at t=0.5

  double tol=1e-14;
  double t = 0.5;
  double x_expected = meshWidth / 2;
  double y_expected = meshWidth / 2 - radius;

  double x,y,xErr,yErr;
  // check value
  circularArc->value(t, x, y);
  xErr = abs(x-x_expected);
  yErr = abs(y-y_expected);
  if (xErr > tol)
  {
    cout << "exact arc x at t=0.5 is incorrect.\n";
    success = false;
  }
  if (yErr > tol)
  {
    cout << "exact arc y at t=0.5 is incorrect.\n";
    success = false;
  }

  // check derivatives
  // figuring out what the x derivative should be is a bit of work, I think,
  // but the y value is at a minimum, so its derivative should be zero
  y_expected = 0;
  circularArc->dt_parametric()->value(t, x, y);
  yErr = abs(y-y_expected);
  if (yErr > tol)
  {
    cout << "exact arc dy/dt at t=0.5 is nonzero.\n";
    success = false;
  }

  shards::CellTopology line_2(shards::getCellTopologyData<shards::Line<2> >() );
  BasisPtr quadraticBasis = BasisFactory::basisFactory()->getBasis(2, line_2.getKey(), Camellia::FUNCTION_SPACE_HGRAD);

  // figure out what the weights for the quadratic "middle node" basis function should be:
  double expected_H1_weight_x, expected_H1_weight_y;
  double expected_L2_weight_x, expected_L2_weight_y;

  FunctionPtr middleBasis;
  {
//.........这里部分代码省略.........

void Boundary::bcsToImpose(FieldContainer<GlobalIndexType> &globalIndices,
                           FieldContainer<Scalar> &globalValues, TBC<Scalar> &bc,
                           DofInterpreter* dofInterpreter)
{
  set< GlobalIndexType > rankLocalCells = _mesh->cellIDsInPartition();
  map< GlobalIndexType, double> bcGlobalIndicesAndValues;

  for (GlobalIndexType cellID : rankLocalCells)
  {
    bcsToImpose(bcGlobalIndicesAndValues, bc, cellID, dofInterpreter);
  }
  singletonBCsToImpose(bcGlobalIndicesAndValues, bc, dofInterpreter);
  
  // ****** New, tag-based BC imposition follows ******
  map< GlobalIndexType, double> bcTagGlobalIndicesAndValues;
  
  map< int, vector<pair<VarPtr, TFunctionPtr<Scalar>>>> tagBCs = bc.getDirichletTagBCs(); // keys are tags
  
  MeshTopology* meshTopo = dynamic_cast<MeshTopology*>(_mesh->getTopology().get());
  
  TEUCHOS_TEST_FOR_EXCEPTION(!meshTopo, std::invalid_argument, "pure MeshTopologyViews are not yet supported by new tag-based BC imposition");
  
  for (auto tagBC : tagBCs)
  {
    int tagID = tagBC.first;
    
    vector<EntitySetPtr> entitySets = meshTopo->getEntitySetsForTagID(DIRICHLET_SET_TAG_NAME, tagID);
    for (EntitySetPtr entitySet : entitySets)
    {
      // get rank-local cells that match the entity set:
      set<IndexType> matchingCellIDs = entitySet->cellIDsThatMatch(_mesh->getTopology(), rankLocalCells);
      for (IndexType cellID : matchingCellIDs)
      {
        ElementTypePtr elemType = _mesh->getElementType(cellID);
        BasisCachePtr basisCache = BasisCache::basisCacheForCell(_mesh, cellID);
        
        for (auto varFunctionPair : tagBC.second)
        {
          VarPtr var = varFunctionPair.first;
          FunctionPtr f = varFunctionPair.second;
          
          vector<int> sideOrdinals = elemType->trialOrderPtr->getSidesForVarID(var->ID());
          
          for (int sideOrdinal : sideOrdinals)
          {
            BasisPtr basis = elemType->trialOrderPtr->getBasis(var->ID(), sideOrdinal);
            bool isVolume = basis->domainTopology()->getDimension() == _mesh->getDimension();
            for (int d=0; d<_mesh->getDimension(); d++)
            {
              vector<unsigned> matchingSubcells;
              if (isVolume)
                matchingSubcells = entitySet->subcellOrdinals(_mesh->getTopology(), cellID, d);
              else
              {
                CellTopoPtr cellTopo = elemType->cellTopoPtr;
                int sideDim = cellTopo->getDimension() - 1;
                vector<unsigned> matchingSubcellsOnSide = entitySet->subcellOrdinalsOnSide(_mesh->getTopology(), cellID, sideOrdinal, d);
                for (unsigned sideSubcellOrdinal : matchingSubcellsOnSide)
                {
                  unsigned cellSubcellOrdinal = CamelliaCellTools::subcellOrdinalMap(cellTopo, sideDim, sideOrdinal, d, sideSubcellOrdinal);
                  matchingSubcells.push_back(cellSubcellOrdinal);
                }
              }
              
              if (matchingSubcells.size() == 0) continue; // nothing to impose
              
              /*
               What follows - projecting the function onto the basis on the whole domain - is more expensive than necessary,
               in the general case: we can do the projection on just the matching subcells, and if we had a way of taking the
               restriction of a basis to a subcell of the domain, then we could avoid computing the whole basis as well.
               
               But for now, this should work, and it's simple to implement.
               */
              BasisCachePtr basisCacheForImposition = isVolume ? basisCache : basisCache->getSideBasisCache(sideOrdinal);
              int numCells = 1;
              FieldContainer<double> basisCoefficients(numCells,basis->getCardinality());
              Projector<double>::projectFunctionOntoBasisInterpolating(basisCoefficients, f, basis, basisCacheForImposition);
              basisCoefficients.resize(basis->getCardinality());
              
              set<GlobalIndexType> matchingGlobalIndices;
              for (unsigned matchingSubcell : matchingSubcells)
              {
                set<GlobalIndexType> subcellGlobalIndices = dofInterpreter->globalDofIndicesForVarOnSubcell(var->ID(),cellID,d,matchingSubcell);
                matchingGlobalIndices.insert(subcellGlobalIndices.begin(),subcellGlobalIndices.end());
              }
              
              FieldContainer<double> globalData;
              FieldContainer<GlobalIndexType> globalDofIndices;
//              dofInterpreter->interpretLocalBasisCoefficients(cellID, var->ID(), sideOrdinal, basisCoefficientsToImpose, globalData, globalDofIndices);
              dofInterpreter->interpretLocalBasisCoefficients(cellID, var->ID(), sideOrdinal, basisCoefficients, globalData, globalDofIndices);
              for (int globalDofOrdinal=0; globalDofOrdinal<globalDofIndices.size(); globalDofOrdinal++)
              {
                GlobalIndexType globalDofIndex = globalDofIndices(globalDofOrdinal);
                if (matchingGlobalIndices.find(globalDofIndex) != matchingGlobalIndices.end())
                  bcTagGlobalIndicesAndValues[globalDofIndex] = globalData(globalDofOrdinal);
              }
            }
          }
        }
      }
//.........这里部分代码省略.........

bool FunctionTests::testBasisSumFunction()
{
  bool success = true;
  // on a single-element mesh, the BasisSumFunction should be identical to
  // the Solution with those coefficients

  // define a new mesh: more interesting if we're not on the ref cell
  int spaceDim = 2;
  FieldContainer<double> quadPoints(4,2);

  quadPoints(0,0) = 0.0; // x1
  quadPoints(0,1) = 0.0; // y1
  quadPoints(1,0) = 2.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 = 1, pToAdd = 0;
  int horizontalCells = 1, verticalCells = 1;

  // create a pointer to a new mesh:
  MeshPtr spectralConfusionMesh = MeshFactory::buildQuadMesh(quadPoints, horizontalCells, verticalCells,
                                  _confusionBF, H1Order, H1Order+pToAdd);

  BCPtr bc = BC::bc();
  SolutionPtr soln = Teuchos::rcp( new Solution(spectralConfusionMesh, bc) );

  soln->initializeLHSVector();

  int cellID = 0;
  double tol = 1e-16; // overly restrictive, just for now.

  DofOrderingPtr trialSpace = spectralConfusionMesh->getElementType(cellID)->trialOrderPtr;
  set<int> trialIDs = trialSpace->getVarIDs();

  BasisCachePtr volumeCache = BasisCache::basisCacheForCell(spectralConfusionMesh, cellID);

  for (set<int>::iterator trialIt=trialIDs.begin(); trialIt != trialIDs.end(); trialIt++)
  {
    int trialID = *trialIt;
    const vector<int>* sidesForVar = &trialSpace->getSidesForVarID(trialID);
    bool boundaryValued = sidesForVar->size() != 1;
    // note that for volume trialIDs, sideIndex = 0, and numSides = 1…
    for (vector<int>::const_iterator sideIt = sidesForVar->begin(); sideIt != sidesForVar->end(); sideIt++)
    {
      int sideIndex = *sideIt;

      BasisCachePtr sideCache = volumeCache->getSideBasisCache(sideIndex);
      BasisPtr basis = trialSpace->getBasis(trialID, sideIndex);
      int basisCardinality = basis->getCardinality();
      for (int basisOrdinal = 0; basisOrdinal<basisCardinality; basisOrdinal++)
      {
        FieldContainer<double> basisCoefficients(basisCardinality);
        basisCoefficients(basisOrdinal) = 1.0;
        soln->setSolnCoeffsForCellID(basisCoefficients, cellID, trialID, sideIndex);

        VarPtr v = Var::varForTrialID(trialID, spectralConfusionMesh->bilinearForm());
        FunctionPtr solnFxn = Function::solution(v, soln, false);
        FunctionPtr basisSumFxn = Teuchos::rcp( new BasisSumFunction(basis, basisCoefficients, Teuchos::rcp((BasisCache*)NULL), OP_VALUE, boundaryValued) );
        if (!boundaryValued)
        {
          double l2diff = (solnFxn - basisSumFxn)->l2norm(spectralConfusionMesh);
//          cout << "l2diff = " << l2diff << endl;
          if (l2diff > tol)
          {
            success = false;
            cout << "testBasisSumFunction: l2diff of " << l2diff << " exceeds tol of " << tol << endl;
            cout << "l2norm of basisSumFxn: " << basisSumFxn->l2norm(spectralConfusionMesh) << endl;
            cout << "l2norm of solnFxn: " << solnFxn->l2norm(spectralConfusionMesh) << endl;
          }
          l2diff = (solnFxn->dx() - basisSumFxn->dx())->l2norm(spectralConfusionMesh);
          //          cout << "l2diff = " << l2diff << endl;
          if (l2diff > tol)
          {
            success = false;
            cout << "testBasisSumFunction: l2diff of dx() " << l2diff << " exceeds tol of " << tol << endl;
            cout << "l2norm of basisSumFxn->dx(): " << basisSumFxn->dx()->l2norm(spectralConfusionMesh) << endl;
            cout << "l2norm of solnFxn->dx(): " << solnFxn->dx()->l2norm(spectralConfusionMesh) << endl;
          }

          // test that the restriction to a side works
          int numSides = volumeCache->cellTopology()->getSideCount();
          for (int i=0; i<numSides; i++)
          {
            BasisCachePtr mySideCache = volumeCache->getSideBasisCache(i);
            if (! solnFxn->equals(basisSumFxn, mySideCache, tol))
            {
              success = false;
              cout << "testBasisSumFunction: on side 0, l2diff of " << l2diff << " exceeds tol of " << tol << endl;
              reportFunctionValueDifferences(solnFxn, basisSumFxn, mySideCache, tol);
            }
            if (! solnFxn->grad(spaceDim)->equals(basisSumFxn->grad(spaceDim), mySideCache, tol))
            {
              success = false;
              cout << "testBasisSumFunction: on side 0, l2diff of dx() " << l2diff << " exceeds tol of " << tol << endl;
              reportFunctionValueDifferences(solnFxn->grad(spaceDim), basisSumFxn->grad(spaceDim), mySideCache, tol);
            }
          }
//.........这里部分代码省略.........