C++ MeshPtr::bilinearForm方法代码示例

MeshPtr MeshTools::timeSliceMesh(MeshPtr spaceTimeMesh, double t,
                                 map<GlobalIndexType, GlobalIndexType> &sliceCellIDToSpaceTimeCellID, int H1OrderForSlice) {
  MeshTopologyPtr meshTopo = spaceTimeMesh->getTopology();
  set<IndexType> cellIDsToCheck = meshTopo->getRootCellIndices();
  set<IndexType> activeCellIDsForTime;
  
  set<IndexType> allActiveCellIDs = meshTopo->getActiveCellIndices();
  
  int spaceDim = meshTopo->getSpaceDim() - 1;  // # of true spatial dimensions
  
  MeshTopologyPtr sliceTopo = Teuchos::rcp( new MeshTopology(spaceDim) );
  set<IndexType> rootCellIDs = meshTopo->getRootCellIndices();
  for (set<IndexType>::iterator rootCellIt = rootCellIDs.begin(); rootCellIt != rootCellIDs.end(); rootCellIt++) {
    IndexType rootCellID = *rootCellIt;
    FieldContainer<double> physicalNodes = spaceTimeMesh->physicalCellNodesForCell(rootCellID);
    if (cellMatches(physicalNodes, t)) { // cell and some subset of its descendents should be included in slice mesh
      vector< vector< double > > sliceNodes = timeSliceForCell(physicalNodes, t);
      CellTopoPtrLegacy cellTopo = getBottomTopology(meshTopo, rootCellID);
      CellPtr sliceCell = sliceTopo->addCell(cellTopo, sliceNodes);
      sliceCellIDToSpaceTimeCellID[sliceCell->cellIndex()] = rootCellID;
    }
  }
  
  MeshPtr sliceMesh = Teuchos::rcp( new Mesh(sliceTopo, spaceTimeMesh->bilinearForm(), H1OrderForSlice, spaceDim) );
  
  // process refinements.  For now, we assume isotropic refinements, which means that each refinement in spacetime induces a refinement in the spatial slice
  set<IndexType> sliceCellIDsToCheckForRefinement = sliceTopo->getActiveCellIndices();
  while (sliceCellIDsToCheckForRefinement.size() > 0) {
    set<IndexType>::iterator cellIt = sliceCellIDsToCheckForRefinement.begin();
    IndexType sliceCellID = *cellIt;
    sliceCellIDsToCheckForRefinement.erase(cellIt);
    
    CellPtr sliceCell = sliceTopo->getCell(sliceCellID);
    CellPtr spaceTimeCell = meshTopo->getCell(sliceCellIDToSpaceTimeCellID[sliceCellID]);
    if (spaceTimeCell->isParent()) {
      set<GlobalIndexType> cellsToRefine;
      cellsToRefine.insert(sliceCellID);
      sliceMesh->hRefine(cellsToRefine, RefinementPattern::regularRefinementPattern(sliceCell->topology()->getKey()));
      vector<IndexType> spaceTimeChildren = spaceTimeCell->getChildIndices();
      for (int childOrdinal=0; childOrdinal<spaceTimeChildren.size(); childOrdinal++) {
        IndexType childID = spaceTimeChildren[childOrdinal];
        FieldContainer<double> childNodes = meshTopo->physicalCellNodesForCell(childID);
        if (cellMatches(childNodes, t)) {
          vector< vector<double> > childSlice = timeSliceForCell(childNodes, t);
          CellPtr childSliceCell = sliceTopo->findCellWithVertices(childSlice);
          sliceCellIDToSpaceTimeCellID[childSliceCell->cellIndex()] = childID;
          sliceCellIDsToCheckForRefinement.insert(childSliceCell->cellIndex());
        }
      }
    }
  }
  
  return sliceMesh;
}

MeshPtr MeshTools::timeSliceMesh(MeshPtr spaceTimeMesh, double t,
                                 map<GlobalIndexType, GlobalIndexType> &sliceCellIDToSpaceTimeCellID, int H1OrderForSlice)
{
  MeshTopology* meshTopo = dynamic_cast<MeshTopology*>(spaceTimeMesh->getTopology().get());
  TEUCHOS_TEST_FOR_EXCEPTION(!meshTopo, std::invalid_argument, "timeSliceMesh() called with spaceTimeMesh that appears to be pure MeshTopologyView.  This is not supported.");
  set<IndexType> cellIDsToCheck = meshTopo->getRootCellIndices();
  set<IndexType> activeCellIDsForTime;

  set<IndexType> allActiveCellIDs = meshTopo->getActiveCellIndices();

  int spaceDim = meshTopo->getDimension() - 1;  // # of true spatial dimensions

  MeshTopologyPtr sliceTopo = Teuchos::rcp( new MeshTopology(spaceDim) );
  set<IndexType> rootCellIDs = meshTopo->getRootCellIndices();
  for (set<IndexType>::iterator rootCellIt = rootCellIDs.begin(); rootCellIt != rootCellIDs.end(); rootCellIt++)
  {
    IndexType rootCellID = *rootCellIt;
    FieldContainer<double> physicalNodes = spaceTimeMesh->physicalCellNodesForCell(rootCellID);
    if (cellMatches(physicalNodes, t))   // cell and some subset of its descendents should be included in slice mesh
    {
      vector< vector< double > > sliceNodes = timeSliceForCell(physicalNodes, t);
      CellTopoPtr cellTopo = getBottomTopology(meshTopo, rootCellID);
      GlobalIndexType newCellID = sliceTopo->cellCount();
      CellPtr sliceCell = sliceTopo->addCell(newCellID, cellTopo, sliceNodes); // for consistency, this is only valid if run on every MPI rank.
      sliceCellIDToSpaceTimeCellID[sliceCell->cellIndex()] = rootCellID;
    }
  }

  MeshPtr sliceMesh = Teuchos::rcp( new Mesh(sliceTopo, spaceTimeMesh->bilinearForm(), H1OrderForSlice, spaceDim) );

  // process refinements.  For now, we assume isotropic refinements, which means that each refinement in spacetime induces a refinement in the spatial slice
  set<IndexType> sliceCellIDsToCheckForRefinement = sliceTopo->getActiveCellIndices();
  while (sliceCellIDsToCheckForRefinement.size() > 0)
  {
    set<IndexType>::iterator cellIt = sliceCellIDsToCheckForRefinement.begin();
    IndexType sliceCellID = *cellIt;
    sliceCellIDsToCheckForRefinement.erase(cellIt);

    CellPtr sliceCell = sliceTopo->getCell(sliceCellID);
    CellPtr spaceTimeCell = meshTopo->getCell(sliceCellIDToSpaceTimeCellID[sliceCellID]);
    if (spaceTimeCell->isParent(spaceTimeMesh->getTopology()))
    {
      set<GlobalIndexType> cellsToRefine;
      cellsToRefine.insert(sliceCellID);
      sliceMesh->hRefine(cellsToRefine, RefinementPattern::regularRefinementPattern(sliceCell->topology()));
      vector<IndexType> spaceTimeChildren = spaceTimeCell->getChildIndices(spaceTimeMesh->getTopology());
      for (int childOrdinal=0; childOrdinal<spaceTimeChildren.size(); childOrdinal++)
      {
        IndexType childID = spaceTimeChildren[childOrdinal];
        FieldContainer<double> childNodes = meshTopo->physicalCellNodesForCell(childID);
        if (cellMatches(childNodes, t))
        {
          vector< vector<double> > childSlice = timeSliceForCell(childNodes, t);
          CellPtr childSliceCell = sliceTopo->findCellWithVertices(childSlice);
          sliceCellIDToSpaceTimeCellID[childSliceCell->cellIndex()] = childID;
          sliceCellIDsToCheckForRefinement.insert(childSliceCell->cellIndex());
        }
      }
    }
  }

  MeshPartitionPolicyPtr partitionPolicy = MeshPartitionPolicy::inducedPartitionPolicy(sliceMesh, spaceTimeMesh, sliceCellIDToSpaceTimeCellID);

  sliceMesh->setPartitionPolicy(partitionPolicy);

  return sliceMesh;
}

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