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

void LagrangeConstraints::getCoefficients(FieldContainer<double> &lhs, FieldContainer<double> &rhs,
    int elemConstraintIndex, DofOrderingPtr trialOrdering,
    BasisCachePtr basisCache)
{
  LinearTermPtr lt = _constraints[elemConstraintIndex].linearTerm();
  TFunctionPtr<double> f = _constraints[elemConstraintIndex].f();
  lt->integrate(lhs, trialOrdering, basisCache);
  bool onBoundary = f->boundaryValueOnly();
  if ( !onBoundary )
  {
    f->integrate(rhs, basisCache);
  }
  else
  {
    int numSides = basisCache->cellTopology()->getSideCount();
    rhs.initialize(0);
    for (int sideIndex=0; sideIndex<numSides; sideIndex++)
    {
      f->integrate(rhs, basisCache->getSideBasisCache(sideIndex), true); // true: sumInto
    }
  }
}

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

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

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

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