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

 void values(FieldContainer<double> &values, BasisCachePtr sliceBasisCache) {
   vector<GlobalIndexType> sliceCellIDs = sliceBasisCache->cellIDs();
   
   Teuchos::Array<int> dim;
   values.dimensions(dim);
   dim[0] = 1; // one cell
   Teuchos::Array<int> offset(dim.size());
   
   for (int cellOrdinal = 0; cellOrdinal < sliceCellIDs.size(); cellOrdinal++) {
     offset[0] = cellOrdinal;
     int enumeration = values.getEnumeration(offset);
     FieldContainer<double>valuesForCell(dim,&values[enumeration]);
     GlobalIndexType sliceCellID = sliceCellIDs[cellOrdinal];
     int numPoints = sliceBasisCache->getPhysicalCubaturePoints().dimension(1);
     int spaceDim = sliceBasisCache->getPhysicalCubaturePoints().dimension(2);
     FieldContainer<double> spaceTimePhysicalPoints(1,numPoints,spaceDim+1);
     for (int ptOrdinal=0; ptOrdinal<numPoints; ptOrdinal++) {
       for (int d=0; d<spaceDim; d++) {
         spaceTimePhysicalPoints(0,ptOrdinal,d) = sliceBasisCache->getPhysicalCubaturePoints()(cellOrdinal,ptOrdinal,d);
       }
       spaceTimePhysicalPoints(0,ptOrdinal,spaceDim) = _t;
     }
     
     GlobalIndexType cellID = _cellIDMap[sliceCellID];
     BasisCachePtr spaceTimeBasisCache = BasisCache::basisCacheForCell(_spaceTimeMesh, cellID);
     
     FieldContainer<double> spaceTimeRefPoints(1,numPoints,spaceDim+1);
     CamelliaCellTools::mapToReferenceFrame(spaceTimeRefPoints, spaceTimePhysicalPoints, _spaceTimeMesh, cellID);
     spaceTimeRefPoints.resize(numPoints,spaceDim+1);
     spaceTimeBasisCache->setRefCellPoints(spaceTimeRefPoints);
     _spaceTimeFunction->values(valuesForCell, spaceTimeBasisCache);
   }
 }

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

double basisSumAtParametricPoint(FieldContainer<double> &basisCoefficients, BasisPtr basis, double tValue, BasisCachePtr parametricBasisCache)
{
  int numPoints = 1;
  int spaceDim = 1;
  FieldContainer<double> parametricPoints(numPoints,spaceDim);
  parametricPoints[0] = tValue;
  FieldContainer<double> refCellPoints = parametricBasisCache->getRefCellPointsForPhysicalPoints(parametricPoints);
  parametricBasisCache->setRefCellPoints(refCellPoints);
  Teuchos::RCP< const FieldContainer<double> > basisValues = parametricBasisCache->getValues(basis, OP_VALUE);
  double sum = 0;
  int ptIndex = 0; // one point
  for (int fieldIndex=0; fieldIndex<basisCoefficients.size(); fieldIndex++)
  {
    sum += (*basisValues)(fieldIndex,ptIndex) * basisCoefficients(fieldIndex);
  }
  return sum;
}

void Projector::projectFunctionOntoBasisInterpolating(FieldContainer<double> &basisCoefficients, FunctionPtr fxn,
                                                      BasisPtr basis, BasisCachePtr domainBasisCache) {
  basisCoefficients.initialize(0);
  CellTopoPtr domainTopo = basis->domainTopology();
  unsigned domainDim = domainTopo->getDimension();
  
  IPPtr ip;
  
  bool traceVar = domainBasisCache->isSideCache();
  
  pair<IPPtr, VarPtr> ipVarPair = IP::standardInnerProductForFunctionSpace(basis->functionSpace(), traceVar, domainDim);
  ip = ipVarPair.first;
  VarPtr v = ipVarPair.second;
  
  IPPtr ip_l2 = Teuchos::rcp( new IP );
  ip_l2->addTerm(v);
  
  // for now, make all projections use L^2... (having some issues with gradients and cell Jacobians--I think we need the restriction of the cell Jacobian to the subcell, e.g., and it's not clear how to do that...)
  ip = ip_l2;
  
  FieldContainer<double> referenceDomainNodes(domainTopo->getVertexCount(),domainDim);
  CamelliaCellTools::refCellNodesForTopology(referenceDomainNodes, domainTopo);
  
  int basisCardinality = basis->getCardinality();
  
  set<int> allDofs;
  for (int i=0; i<basisCardinality; i++) {
    allDofs.insert(i);
  }
  
  for (int d=0; d<=domainDim; d++) {
    FunctionPtr projectionThusFar = NewBasisSumFunction::basisSumFunction(basis, basisCoefficients);
    FunctionPtr fxnToApproximate = fxn - projectionThusFar;
    int subcellCount = domainTopo->getSubcellCount(d);
    for (int subcord=0; subcord<subcellCount; subcord++) {
      set<int> subcellDofOrdinals = basis->dofOrdinalsForSubcell(d, subcord);
      if (subcellDofOrdinals.size() > 0) {
        FieldContainer<double> refCellPoints;
        FieldContainer<double> cubatureWeightsSubcell; // allows us to integrate over the fine subcell even when domain is higher-dimensioned
        if (d == 0) {
          refCellPoints.resize(1,domainDim);
          for (int d1=0; d1<domainDim; d1++) {
            refCellPoints(0,d1) = referenceDomainNodes(subcord,d1);
          }
          cubatureWeightsSubcell.resize(1);
          cubatureWeightsSubcell(0) = 1.0;
        } else {
          CellTopoPtr subcellTopo = domainTopo->getSubcell(d, subcord);
//          Teuchos::RCP<Cubature<double> > subcellCubature = cubFactory.create(subcellTopo, domainBasisCache->cubatureDegree());
          BasisCachePtr subcellCache = Teuchos::rcp( new BasisCache(subcellTopo, domainBasisCache->cubatureDegree(), false) );
          int numPoints = subcellCache->getRefCellPoints().dimension(0);
          refCellPoints.resize(numPoints,domainDim);
          cubatureWeightsSubcell = subcellCache->getCubatureWeights();
          
          if (d == domainDim) {
            refCellPoints = subcellCache->getRefCellPoints();
          } else {
            CamelliaCellTools::mapToReferenceSubcell(refCellPoints, subcellCache->getRefCellPoints(), d,
                                                     subcord, domainTopo);
          }
        }
        domainBasisCache->setRefCellPoints(refCellPoints, cubatureWeightsSubcell);
        IPPtr ipForProjection = (d==0) ? ip_l2 : ip; // just use values at vertices (ignore derivatives)
        set<int> dofsToSkip = allDofs;
        for (set<int>::iterator dofOrdinalIt=subcellDofOrdinals.begin(); dofOrdinalIt != subcellDofOrdinals.end(); dofOrdinalIt++) {
          dofsToSkip.erase(*dofOrdinalIt);
        }
        FieldContainer<double> newBasisCoefficients;
        projectFunctionOntoBasis(newBasisCoefficients, fxnToApproximate, basis, domainBasisCache, ipForProjection, v, dofsToSkip);
        for (int cellOrdinal=0; cellOrdinal<newBasisCoefficients.dimension(0); cellOrdinal++) {
          for (set<int>::iterator dofOrdinalIt=subcellDofOrdinals.begin(); dofOrdinalIt != subcellDofOrdinals.end(); dofOrdinalIt++) {
            basisCoefficients(cellOrdinal,*dofOrdinalIt) = newBasisCoefficients(cellOrdinal,*dofOrdinalIt);
          }
        }
      }
    }
  }
}

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

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

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

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

  int numCubPoints = fineSideTopoCub->getNumPoints();

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

  fineSideTopoCub->getCubature(cubPoints, cubWeights);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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