C++ CellPtr::getNeighborInfo方法代码示例

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

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

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

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

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

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

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

      unsigned numTestPoints = coarseCellRefPoints.dimension(0);

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

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

      BasisCachePtr fineSideBasisCache = fineCellBasisCache->getSideBasisCache(sideOrdinal);

      fineCellBasisCache->setRefCellPoints(fineCellRefPoints);

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

      coarseCellBasisCache->setRefCellPoints(coarseCellRefPoints);

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

      // sanity check: do physical points match?

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

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

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

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

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

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

bool MeshTestUtility::determineRefTestPointsForNeighbors(MeshTopologyViewPtr meshTopo, CellPtr fineCell, unsigned int sideOrdinal,
    FieldContainer<double> &fineSideRefPoints, FieldContainer<double> &fineCellRefPoints,
    FieldContainer<double> &coarseSideRefPoints, FieldContainer<double> &coarseCellRefPoints)
{
  unsigned spaceDim = meshTopo->getDimension();
  unsigned sideDim = spaceDim - 1;

  if (spaceDim == 1)
  {
    fineSideRefPoints.resize(0,0);
    coarseSideRefPoints.resize(0,0);
    fineCellRefPoints.resize(1,1);
    coarseCellRefPoints.resize(1,1);

    FieldContainer<double> lineRefNodes(2,1);
    CellTopoPtr line = CellTopology::line();

    CamelliaCellTools::refCellNodesForTopology(lineRefNodes, line);

    fineCellRefPoints[0] = lineRefNodes[sideOrdinal];
    unsigned neighborSideOrdinal = fineCell->getNeighborInfo(sideOrdinal, meshTopo).second;
    if (neighborSideOrdinal != -1)
    {
      coarseCellRefPoints[0] = lineRefNodes[neighborSideOrdinal];
      return true;
    }
    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);
//.........这里部分代码省略.........