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

  void values(FieldContainer<double> &values, BasisCachePtr basisCache)
  {
    int numCells = values.dimension(0);
    int numPoints = values.dimension(1);

    MeshPtr mesh = basisCache->mesh();
    vector<int> cellIDs = basisCache->cellIDs();

    double tol=1e-14;
    for (int cellIndex=0; cellIndex<numCells; cellIndex++)
    {
      double h = 1.0;
      if (_spatialCoord==0)
      {
        h = mesh->getCellXSize(cellIDs[cellIndex]);
      }
      else if (_spatialCoord==1)
      {
        h = mesh->getCellYSize(cellIDs[cellIndex]);
      }
      for (int ptIndex=0; ptIndex<numPoints; ptIndex++)
      {
        values(cellIndex,ptIndex) = sqrt(h);
      }
    }
  }

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);
      }
    }
  }
}

void SideParityFunction::values(Intrepid::FieldContainer<double> &values, BasisCachePtr sideBasisCache)
{
  this->CHECK_VALUES_RANK(values);
  int numCells = values.dimension(0);
  int numPoints = values.dimension(1);
  int sideIndex = sideBasisCache->getSideIndex();
  if (sideIndex == -1)
  {
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "non-sideBasisCache passed into SideParityFunction");
  }
  if (sideBasisCache->getCellSideParities().size() > 0)
  {
    // then we'll use this, and won't require that mesh and cellIDs are set
    if (sideBasisCache->getCellSideParities().dimension(0) != numCells)
    {
      TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "sideBasisCache->getCellSideParities() is non-empty, but the cell dimension doesn't match that of the values FieldContainer.");
    }

    for (int cellOrdinal=0; cellOrdinal<numCells; cellOrdinal++)
    {
      int parity = sideBasisCache->getCellSideParities()(cellOrdinal,sideIndex);
      for (int ptOrdinal=0; ptOrdinal<numPoints; ptOrdinal++)
      {
        values(cellOrdinal,ptOrdinal) = parity;
      }
    }
  }
  else
  {
    vector<GlobalIndexType> cellIDs = sideBasisCache->cellIDs();
    if (cellIDs.size() != numCells)
    {
      TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "cellIDs.size() != numCells");
    }
    Teuchos::RCP<Mesh> mesh = sideBasisCache->mesh();
    if (! mesh.get())
    {
      TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "mesh unset in BasisCache.");
    }
    for (int cellIndex=0; cellIndex<numCells; cellIndex++)
    {
      int parity = mesh->cellSideParitiesForCell(cellIDs[cellIndex])(0,sideIndex);
      for (int ptIndex=0; ptIndex<numPoints; ptIndex++)
      {
        values(cellIndex,ptIndex) = parity;
      }
    }
  }
}

  void values(FieldContainer<double> &values, BasisCachePtr basisCache){
    MeshPtr mesh = basisCache->mesh();
    vector<int> cellIDs = basisCache->cellIDs();
    int numPoints = values.dimension(1);
    for (int i = 0;i<cellIDs.size();i++){
      double h = mesh->getCellXSize(cellIDs[i]); // default to x-direction
      if (_xyDim == 1){
	h = mesh->getCellYSize(cellIDs[i]);
      }
      double hPower = pow(h,_power);      
      for (int j = 0;j<numPoints;j++){
	values(i,j) = hPower;
      }
    }
  }

    virtual void localStiffnessMatrixAndRHS(FieldContainer<double> &localStiffness, FieldContainer<double> &rhsVector,
                                            IPPtr ip, BasisCachePtr ipBasisCache,
                                            RHSPtr rhs, BasisCachePtr basisCache) {
        double testMatrixAssemblyTime = 0, testMatrixInversionTime = 0, localStiffnessDeterminationFromTestsTime = 0;

#ifdef HAVE_MPI
        Epetra_MpiComm Comm(MPI_COMM_WORLD);
        //cout << "rank: " << rank << " of " << numProcs << endl;
#else
        Epetra_SerialComm Comm;
#endif

        Epetra_Time timer(Comm);

        // localStiffness should have dim. (numCells, numTrialFields, numTrialFields)
        MeshPtr mesh = basisCache->mesh();
        if (mesh.get() == NULL) {
            cout << "localStiffnessMatrix requires BasisCache to have mesh set.\n";
            TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "localStiffnessMatrix requires BasisCache to have mesh set.");
        }
        const vector<GlobalIndexType>* cellIDs = &basisCache->cellIDs();
        int numCells = cellIDs->size();
        if (numCells != localStiffness.dimension(0)) {
            cout << "localStiffnessMatrix requires basisCache->cellIDs() to have the same # of cells as the first dimension of localStiffness\n";
            TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "localStiffnessMatrix requires basisCache->cellIDs() to have the same # of cells as the first dimension of localStiffness");
        }

        ElementTypePtr elemType = mesh->getElementType((*cellIDs)[0]); // we assume all cells provided are of the same type
        DofOrderingPtr trialOrder = elemType->trialOrderPtr;
        DofOrderingPtr fieldOrder = mesh->getDofOrderingFactory().getFieldOrdering(trialOrder);
        DofOrderingPtr traceOrder = mesh->getDofOrderingFactory().getTraceOrdering(trialOrder);

        map<int, int> stiffnessIndexForTraceIndex;
        map<int, int> stiffnessIndexForFieldIndex;
        set<int> varIDs = trialOrder->getVarIDs();
        for (set<int>::iterator varIt = varIDs.begin(); varIt != varIDs.end(); varIt++) {
            int varID = *varIt;
            const vector<int>* sidesForVar = &trialOrder->getSidesForVarID(varID);
            bool isTrace = (sidesForVar->size() > 1);
            for (vector<int>::const_iterator sideIt = sidesForVar->begin(); sideIt != sidesForVar->end(); sideIt++) {
                int sideOrdinal = *sideIt;
                vector<int> dofIndices = trialOrder->getDofIndices(varID,sideOrdinal);
                if (isTrace) {
                    vector<int> traceDofIndices = traceOrder->getDofIndices(varID,sideOrdinal);
                    for (int i=0; i<traceDofIndices.size(); i++) {
                        stiffnessIndexForTraceIndex[traceDofIndices[i]] = dofIndices[i];
                    }
                } else {
                    vector<int> fieldDofIndices = fieldOrder->getDofIndices(varID);
                    for (int i=0; i<fieldDofIndices.size(); i++) {
                        stiffnessIndexForFieldIndex[fieldDofIndices[i]] = dofIndices[i];
                    }
                }
            }
        }

        int numTrialDofs = trialOrder->totalDofs();
        if ((numTrialDofs != localStiffness.dimension(1)) || (numTrialDofs != localStiffness.dimension(2))) {
            cout << "localStiffness should have dimensions (C,numTrialFields,numTrialFields).\n";
            TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "localStiffness should have dimensions (C,numTrialFields,numTrialFields).");
        }

        map<int,int> traceTestMap, fieldTestMap;
        int numEquations = _virtualTerms.getFieldTestVars().size();
        for (int eqn=0; eqn<numEquations; eqn++) {
            VarPtr testVar = _virtualTerms.getFieldTestVars()[eqn];
            VarPtr traceVar = _virtualTerms.getAssociatedTrace(testVar);
            VarPtr fieldVar = _virtualTerms.getAssociatedField(testVar);
            traceTestMap[traceVar->ID()] = testVar->ID();
            fieldTestMap[fieldVar->ID()] = testVar->ID();
        }

        int maxDegreeField = fieldOrder->maxBasisDegree();
        int testDegreeInterior = maxDegreeField + _virtualTerms.getTestEnrichment();
        int testDegreeTrace = testDegreeInterior + 2;

        cout << "ERROR in Virtual: getRelabeledDofOrdering() is commented out in DofOrderingFactory.  Need to rewrite for the new caching scheme.\n";
        DofOrderingPtr testOrderInterior; // = mesh->getDofOrderingFactory().getRelabeledDofOrdering(fieldOrder, fieldTestMap);
        testOrderInterior = mesh->getDofOrderingFactory().setBasisDegree(testOrderInterior, testDegreeInterior, false);
        DofOrderingPtr testOrderTrace = mesh->getDofOrderingFactory().setBasisDegree(testOrderInterior, testDegreeTrace, true); // this has a bunch of extra dofs (interior guys)

        map<int, int> remappedTraceIndices; // go from the index that includes the interior dofs to one that doesn't
        set<int> testIDs = testOrderTrace->getVarIDs();
        int testTraceIndex = 0;
        for (set<int>::iterator testIDIt = testIDs.begin(); testIDIt != testIDs.end(); testIDIt++) {
            int testID = *testIDIt;
            BasisPtr basis = testOrderTrace->getBasis(testID);
            vector<int> interiorDofs = basis->dofOrdinalsForInterior();
            for (int basisOrdinal=0; basisOrdinal<basis->getCardinality(); basisOrdinal++) {
                if (std::find(interiorDofs.begin(),interiorDofs.end(),basisOrdinal) == interiorDofs.end()) {
                    int dofIndex = testOrderTrace->getDofIndex(testID, basisOrdinal);
                    remappedTraceIndices[dofIndex] = testTraceIndex;
                    testTraceIndex++;
                }
            }
        }

//    DofOrderingPtr testOrderTrace = mesh->getDofOrderingFactory().getRelabeledDofOrdering(traceOrder, traceTestMap);
//    testOrderTrace = mesh->getDofOrderingFactory().setBasisDegree(testOrderTrace, testDegreeTrace);

//.........这里部分代码省略.........