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

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 MeshTopologyTests::testEntityConstraints()
{
  bool success = true;

  // make two simple meshes
  MeshTopologyPtr mesh2D = makeRectMesh(0.0, 0.0, 2.0, 1.0,
                                        2, 1);
  MeshTopologyPtr mesh3D = makeHexMesh(0.0, 0.0, 0.0, 2.0, 4.0, 3.0,
                                       2, 2, 1);

  unsigned vertexDim = 0;
  unsigned edgeDim = 1;
  unsigned faceDim = 2;

  // first, check that unconstrained edges and faces are unconstrained

  set< unsigned > boundaryEdges;
  set< unsigned > internalEdges;

  for (unsigned cellIndex=0; cellIndex<mesh2D->cellCount(); cellIndex++)
  {
    CellPtr cell = mesh2D->getCell(cellIndex);
    unsigned sideCount = cell->getSideCount();

    for (unsigned sideOrdinal=0; sideOrdinal<sideCount; sideOrdinal++)
    {
      unsigned edgeIndex = cell->entityIndex(edgeDim, sideOrdinal);
      unsigned numCells = mesh2D->getActiveCellCount(edgeDim,edgeIndex);
      if (numCells == 1)   // boundary edge
      {
        boundaryEdges.insert(edgeIndex);
      }
      else if (numCells == 2)
      {
        internalEdges.insert(edgeIndex);
      }
      else
      {
        success = false;
        cout << "testEntityConstraints: In initial 2D mesh, edge " << edgeIndex << " has active cell count of " << numCells << ".\n";
      }
    }
  }
  if (internalEdges.size() != 1)
  {
    success = false;
    cout << "testEntityConstraints: In initial 2D mesh, there are " << internalEdges.size() << " internal edges (expected 1).\n";
  }
  for (set<unsigned>::iterator edgeIt=internalEdges.begin(); edgeIt != internalEdges.end(); edgeIt++)
  {
    unsigned edgeIndex = *edgeIt;
    unsigned constrainingEntityIndex = mesh2D->getConstrainingEntity(edgeDim,edgeIndex).first;
    if (constrainingEntityIndex != edgeIndex)
    {
      success = false;
      cout << "testEntityConstraints: In initial 2D mesh, internal edge is constrained by a different edge.\n";
    }
  }

  set<unsigned> boundaryFaces;
  set<unsigned> internalFaces;
  map<unsigned, vector<unsigned> > faceToEdges;
  for (unsigned cellIndex=0; cellIndex<mesh3D->cellCount(); cellIndex++)
  {
    CellPtr cell = mesh3D->getCell(cellIndex);
    unsigned sideCount = cell->getSideCount();

    for (unsigned sideOrdinal=0; sideOrdinal<sideCount; sideOrdinal++)
    {
      unsigned faceIndex = cell->entityIndex(faceDim, sideOrdinal);
      unsigned numCells = mesh3D->getActiveCellCount(faceDim,faceIndex);
      if (numCells == 1)   // boundary face
      {
        boundaryFaces.insert(faceIndex);
      }
      else if (numCells == 2)
      {
        internalFaces.insert(faceIndex);
      }
      else
      {
        success = false;
        cout << "testEntityConstraints: In initial 3D mesh, face " << faceIndex << " has active cell count of " << numCells << ".\n";
      }

      if (faceToEdges.find(faceIndex) == faceToEdges.end())
      {
        CellTopoPtr faceTopo = cell->topology()->getSubcell(faceDim, sideOrdinal);
        unsigned numEdges = faceTopo->getSubcellCount(edgeDim);
        vector<unsigned> edgeIndices(numEdges);
        for (unsigned edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
        {
          edgeIndices[edgeOrdinal] = mesh3D->getFaceEdgeIndex(faceIndex, edgeOrdinal);
        }
      }
    }
  }

  if (internalFaces.size() != 4)
  {
//.........这里部分代码省略.........