C++ BasisCachePtr类代码示例

  bool matchesPoints(FieldContainer<bool> &pointsMatch, BasisCachePtr basisCache)
  {
    const FieldContainer<double> *points = &(basisCache->getPhysicalCubaturePoints());
    const FieldContainer<double> *normals = &(basisCache->getSideNormals());
    int numCells = (*points).dimension(0);
    int numPoints = (*points).dimension(1);

    double tol = 1e-3;
    bool somePointMatches = false;
    for (int cellIndex=0; cellIndex<numCells; cellIndex++)
    {
      for (int ptIndex=0; ptIndex<numPoints; ptIndex++)
      {
        double x = (*points)(cellIndex,ptIndex,0);
        double y = (*points)(cellIndex,ptIndex,1);
        double n1 = (*normals)(cellIndex,ptIndex,0);
        double n2 = (*normals)(cellIndex,ptIndex,1);
        double beta_n = _beta[0]*n1 + _beta[1]*n2 ;
        pointsMatch(cellIndex,ptIndex) = false;
        if (abs((x-.5)*(x-.5)+y*y) < 0.75+tol && beta_n < 0)
        {
          pointsMatch(cellIndex,ptIndex) = true;
          somePointMatches = true;
        }
      }
    }
    return somePointMatches;
  }

void FunctionTests::checkFunctionsAgree(FunctionPtr f1, FunctionPtr f2, BasisCachePtr basisCache) {
  ASSERT_EQ(f1->rank(), f2->rank())
    << "f1->rank() " << f1->rank() << " != f2->rank() " << f2->rank() << endl;

  int rank = f1->rank();
  int numCells = basisCache->getPhysicalCubaturePoints().dimension(0);
  int numPoints = basisCache->getPhysicalCubaturePoints().dimension(1);
  int spaceDim = basisCache->getPhysicalCubaturePoints().dimension(2);
  Teuchos::Array<int> dim;
  dim.append(numCells);
  dim.append(numPoints);
  for (int i=0; i<rank; i++) {
    dim.append(spaceDim);
  }
  FieldContainer<double> f1Values(dim);
  FieldContainer<double> f2Values(dim);
  f1->values(f1Values,basisCache);
  f2->values(f2Values,basisCache);
  
  double tol = 1e-14;
  double maxDiff;
  EXPECT_TRUE(fcsAgree(f1Values,f2Values,tol,maxDiff))
    << "Test failed: f1 and f2 disagree; maxDiff " << maxDiff << ".\n"
    << "f1Values: \n" << f1Values
    << "f2Values: \n" << f2Values;
}

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

  void values(FieldContainer<double> &values, BasisCachePtr basisCache)
  {
    int numCells = values.dimension(0);
    int numPoints = values.dimension(1);

    MeshPtr mesh = basisCache->mesh();
    vector<int> cellIDs = basisCache->cellIDs();

    double tol=1e-14;
    for (int cellIndex=0; cellIndex<numCells; cellIndex++)
    {
      double h = 1.0;
      if (_spatialCoord==0)
      {
        h = mesh->getCellXSize(cellIDs[cellIndex]);
      }
      else if (_spatialCoord==1)
      {
        h = mesh->getCellYSize(cellIDs[cellIndex]);
      }
      for (int ptIndex=0; ptIndex<numPoints; ptIndex++)
      {
        values(cellIndex,ptIndex) = sqrt(h);
      }
    }
  }

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

  // bool matchesPoint(double x, double y) {
  //   return true;
  // }
  bool matchesPoints(FieldContainer<bool> &pointsMatch, BasisCachePtr basisCache)
  {
    const FieldContainer<double> *points = &(basisCache->getPhysicalCubaturePoints());
    const FieldContainer<double> *normals = &(basisCache->getSideNormals());
    int numCells = (*points).dimension(0);
    int numPoints = (*points).dimension(1);

    FieldContainer<double> beta_pts(numCells,numPoints,2);
    _beta->values(beta_pts,basisCache);

    double tol=1e-14;
    bool somePointMatches = false;
    for (int cellIndex=0; cellIndex<numCells; cellIndex++)
    {
      for (int ptIndex=0; ptIndex<numPoints; ptIndex++)
      {
        double n1 = (*normals)(cellIndex,ptIndex,0);
        double n2 = (*normals)(cellIndex,ptIndex,1);
        double beta_n = beta_pts(cellIndex,ptIndex,0)*n1 + beta_pts(cellIndex,ptIndex,1)*n2 ;
        // cout << "beta = (" << beta_pts(cellIndex,ptIndex,0)<<", "<<
        //   beta_pts(cellIndex,ptIndex,1)<<"), n = ("<<n1<<", "<<n2<<")"<<endl;
        pointsMatch(cellIndex,ptIndex) = false;
        if (beta_n < 0)
        {
          pointsMatch(cellIndex,ptIndex) = true;
          somePointMatches = true;
        }
      }
    }
    return somePointMatches;
  }

void MeshTransformationFunction::values(FieldContainer<double> &values, BasisCachePtr basisCache)
{
  CHECK_VALUES_RANK(values);
  vector<GlobalIndexType> cellIDs = basisCache->cellIDs();
  // identity map is the right thing most of the time
  // we'll do something different only where necessary
  int spaceDim = values.dimension(2);
  TEUCHOS_TEST_FOR_EXCEPTION(basisCache->cellTopology()->getDimension() != spaceDim, std::invalid_argument, "cellTopology dimension does not match the shape of the values container");
  if (_op == OP_VALUE)
  {
    values = basisCache->getPhysicalCubaturePoints(); // identity
  }
  else if (_op == OP_DX)
  {
    // identity map is 1 in all the x slots, 0 in all others
    int mod_value = 0; // the x slots are the mod spaceDim = 0 slots;
    for (int i=0; i<values.size(); i++)
    {
      values[i] = (i%spaceDim == mod_value) ? 1.0 : 0.0;
    }
  }
  else if (_op == OP_DY)
  {
    // identity map is 1 in all the y slots, 0 in all others
    int mod_value = 1; // the y slots are the mod spaceDim = 1 slots;
    for (int i=0; i<values.size(); i++)
    {
      values[i] = (i%spaceDim == mod_value) ? 1.0 : 0.0;
    }
  }
  else if (_op == OP_DZ)
  {
    // identity map is 1 in all the z slots, 0 in all others
    int mod_value = 2; // the z slots are the mod spaceDim = 2 slots;
    for (int i=0; i<values.size(); i++)
    {
      values[i] = (i%spaceDim == mod_value) ? 1.0 : 0.0;
    }
  }
//  if (_op == OP_DX) {
//    cout << "values before cellTransformation:\n" << values;
//  }
  for (int cellIndex=0; cellIndex < cellIDs.size(); cellIndex++)
  {
    GlobalIndexType cellID = cellIDs[cellIndex];
    if (_cellTransforms.find(cellID) == _cellTransforms.end()) continue;
    TFunctionPtr<double> cellTransformation = _cellTransforms[cellID];
    ((CellTransformationFunction*)cellTransformation.get())->setCellIndex(cellIndex);
    cellTransformation->values(values, basisCache);
  }
//  if (_op == OP_DX) {
//    cout << "values after cellTransformation:\n" << values;
//  }
}

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

 void values(FieldContainer<double> &values, BasisCachePtr basisCache)
 {
   vector<GlobalIndexType> cellIDs = basisCache->cellIDs();
   int numPoints = values.dimension(1);
   FieldContainer<double> points = basisCache->getPhysicalCubaturePoints();
   for (int i = 0; i<cellIDs.size(); i++)
   {
     for (int j = 0; j<numPoints; j++)
     {
       values(i,j) = cellIDs[i];
     }
   }
 }

void UnitNormalFunction::values(Intrepid::FieldContainer<double> &values, BasisCachePtr basisCache)
{
    this->CHECK_VALUES_RANK(values);
    int numCells = values.dimension(0);
    int numPoints = values.dimension(1);
    int spaceDim = basisCache->getSpaceDim();
    if (_comp == -1)
    {
        // check the the "D" dimension of values is correct:
        if (_spaceTime)
        {
            TEUCHOS_TEST_FOR_EXCEPTION(values.dimension(2) != spaceDim+1, std::invalid_argument, "For space-time normals, values.dimension(2) should be spaceDim + 1.");
        }
        else
        {
            TEUCHOS_TEST_FOR_EXCEPTION(values.dimension(2) != spaceDim, std::invalid_argument, "For spatial normals, values.dimension(2) should be spaceDim.");
        }
    }
    const Intrepid::FieldContainer<double> *sideNormals = _spaceTime ? &(basisCache->getSideNormalsSpaceTime()) : &(basisCache->getSideNormals());

    int comp = _comp;
    if (comp == -2)
    {
        // want to select the temporal component, t()
        comp = spaceDim;
    }
    for (int cellIndex=0; cellIndex<numCells; cellIndex++)
    {
        for (int ptIndex=0; ptIndex<numPoints; ptIndex++)
        {
            if (comp == -1)
            {
                for (int d=0; d<spaceDim; d++)
                {
                    double nd = (*sideNormals)(cellIndex,ptIndex,d);
                    values(cellIndex,ptIndex,d) = nd;
                }
                if (_spaceTime)
                {
                    double nd = (*sideNormals)(cellIndex,ptIndex,spaceDim);
                    values(cellIndex,ptIndex,spaceDim) = nd;
                }
            }
            else
            {
                double ni = (*sideNormals)(cellIndex,ptIndex,comp);
                values(cellIndex,ptIndex) = ni;
            }
        }
    }
}

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

void ExactSolution::solutionValues(FieldContainer<double> &values, int trialID, BasisCachePtr basisCache) {
  if (_exactFunctions.find(trialID) != _exactFunctions.end() ) {
    _exactFunctions[trialID]->values(values,basisCache);
    return;
  }
  
  // TODO: change ExactSolution::solutionValues (below) to take a *const* points FieldContainer, to avoid this copy:
  FieldContainer<double> points = basisCache->getPhysicalCubaturePoints();
  if (basisCache->getSideIndex() >= 0) {
    FieldContainer<double> unitNormals = basisCache->getSideNormals();
    this->solutionValues(values,trialID,points,unitNormals);
  } else {
    this->solutionValues(values,trialID,points);
  }
}

  void values(FieldContainer<double> &values, BasisCachePtr basisCache){
    MeshPtr mesh = basisCache->mesh();
    vector<int> cellIDs = basisCache->cellIDs();
    int numPoints = values.dimension(1);
    for (int i = 0;i<cellIDs.size();i++){
      double h = mesh->getCellXSize(cellIDs[i]); // default to x-direction
      if (_xyDim == 1){
	h = mesh->getCellYSize(cellIDs[i]);
      }
      double hPower = pow(h,_power);      
      for (int j = 0;j<numPoints;j++){
	values(i,j) = hPower;
      }
    }
  }

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

    void values(FieldContainer<double> &values, BasisCachePtr basisCache) {
      int numCells = values.dimension(0);
      int numPoints = values.dimension(1);

      const FieldContainer<double> *points = &(basisCache->getPhysicalCubaturePoints());
      for (int cellIndex=0; cellIndex<numCells; cellIndex++) {
        double xCenter = 0;
        double yCenter = 0;
        int nPts = 0;
        for (int ptIndex=0; ptIndex<numPoints; ptIndex++) {
          double x = (*points)(cellIndex,ptIndex,0);
          double y = (*points)(cellIndex,ptIndex,1);
          xCenter += x;
          yCenter += y;
          nPts++;
        }
        xCenter /= nPts;
        yCenter /= nPts;
        for (int ptIndex=0; ptIndex<numPoints; ptIndex++) {
          double x = (*points)(cellIndex,ptIndex,0);
          double y = (*points)(cellIndex,ptIndex,1);
          if (abs(y) <= halfWidth && abs(yCenter) < halfWidth)
            values(cellIndex, ptIndex) = 1.0;
          else
            values(cellIndex, ptIndex) = 0.0;
        }
      }
    }