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

// computes riesz representation over a single element - map is from int (testID) to FieldContainer of values (sized cellIndex, numPoints)
void RieszRep::computeRepresentationValues(FieldContainer<double> &values, int testID, IntrepidExtendedTypes::EOperatorExtended op, BasisCachePtr basisCache){

  if (_repsNotComputed){
    cout << "Computing riesz rep dofs" << endl;
    computeRieszRep();
  }

  int spaceDim = _mesh->getTopology()->getSpaceDim();
  int numCells = values.dimension(0);
  int numPoints = values.dimension(1);
  vector<GlobalIndexType> cellIDs = basisCache->cellIDs();

  // all elems coming in should be of same type
  ElementPtr elem = _mesh->getElement(cellIDs[0]);
  ElementTypePtr elemTypePtr = elem->elementType();   
  DofOrderingPtr testOrderingPtr = elemTypePtr->testOrderPtr;
  CellTopoPtrLegacy cellTopoPtr = elemTypePtr->cellTopoPtr;
  int numTestDofsForVarID = testOrderingPtr->getBasisCardinality(testID, 0);
  BasisPtr testBasis = testOrderingPtr->getBasis(testID);
  
  bool testBasisIsVolumeBasis = (spaceDim == testBasis->domainTopology()->getDimension());  
  bool useCubPointsSideRefCell = testBasisIsVolumeBasis && basisCache->isSideCache();
  
  Teuchos::RCP< const FieldContainer<double> > transformedBasisValues = basisCache->getTransformedValues(testBasis,op,useCubPointsSideRefCell);
  
  int rank = values.rank() - 2; // if values are shaped as (C,P), scalar...
  if (rank > 1) {
    cout << "ranks greater than 1 not presently supported...\n";
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "ranks greater than 1 not presently supported...");
  }
  
//  Camellia::print("cellIDs",cellIDs);
  
  values.initialize(0.0);
  for (int cellIndex = 0;cellIndex<numCells;cellIndex++){
    int cellID = cellIDs[cellIndex];
    for (int j = 0;j<numTestDofsForVarID;j++) {
      int dofIndex = testOrderingPtr->getDofIndex(testID, j);
      for (int i = 0;i<numPoints;i++) {
        if (rank==0) {
          double basisValue = (*transformedBasisValues)(cellIndex,j,i);
          values(cellIndex,i) += basisValue*_rieszRepDofsGlobal[cellID](dofIndex);
        } else {
          for (int d = 0; d<spaceDim; d++) {
            double basisValue = (*transformedBasisValues)(cellIndex,j,i,d);
            values(cellIndex,i,d) += basisValue*_rieszRepDofsGlobal[cellID](dofIndex);
          }
        }
      }
    }
  }
//  TestSuite::serializeOutput("rep values", values);
}

void PenaltyMethodFilter::filter(FieldContainer<double> &localStiffnessMatrix, FieldContainer<double> &localRHSVector,
                                 BasisCachePtr basisCache, Teuchos::RCP<Mesh> mesh, Teuchos::RCP<BC> bc){ 
  
  // assumption: filter gets elements of all the same type  
  TEUCHOS_TEST_FOR_EXCEPTION(basisCache->cellIDs().size()==0,std::invalid_argument,"no cell IDs given to filter");
  
  ElementTypePtr elemTypePtr = mesh->getElement(basisCache->cellIDs()[0])->elementType();
  int numCells = localStiffnessMatrix.dimension(0);
  
  DofOrderingPtr trialOrderPtr = elemTypePtr->trialOrderPtr;
  
  unsigned numSides = CamelliaCellTools::getSideCount( *elemTypePtr->cellTopoPtr );
  // only allows for L2 inner products at the moment. 
  IntrepidExtendedTypes::EOperatorExtended trialOperator =  IntrepidExtendedTypes::OP_VALUE;
	
  // loop over sides first 
  for (unsigned int sideIndex = 0; sideIndex<numSides; sideIndex++){
    
    // GET INTEGRATION INFO - get cubature points and side normals to send to Constraints (Cell,Point, spaceDim)
    FieldContainer<double> sideCubPoints = basisCache->getPhysicalCubaturePointsForSide(sideIndex);
    FieldContainer<double> sideNormals = basisCache->getSideUnitNormals(sideIndex);        
    
    int numPts = sideCubPoints.dimension(1);
    
    // GET CONSTRAINT INFO
    vector<map<int, FieldContainer<double> > > constrCoeffsVector;
    vector<FieldContainer<double> > constraintValuesVector;
    vector<FieldContainer<bool> > imposeHereVector;
    _constraints->getConstraints(sideCubPoints,sideNormals,constrCoeffsVector,constraintValuesVector);
    
    //loop thru constraints
    int i = 0;
    for (vector<map<int,FieldContainer<double> > >::iterator constrIt = constrCoeffsVector.begin();
         constrIt !=constrCoeffsVector.end(); constrIt++) {
      map<int,FieldContainer<double> > constrCoeffs = *constrIt;
      FieldContainer<double> constrValues = constraintValuesVector[i];
      i++;
      
      double penaltyParameter = 1e7; // (single_precision)^(-1) - perhaps have this computed relative to terms in the matrix?
      
      for (map<int,FieldContainer<double> >::iterator constrTestIDIt = constrCoeffs.begin();
           constrTestIDIt !=constrCoeffs.end(); constrTestIDIt++) {
        pair<int,FieldContainer<double> > constrTestPair = *constrTestIDIt;
        int testTrialID = constrTestPair.first;
        
        // get basis to integrate for testing fxns
        BasisPtr testTrialBasis = trialOrderPtr->getBasis(testTrialID,sideIndex);
        FieldContainer<double> testTrialValuesTransformedWeighted = *(basisCache->getTransformedWeightedValues(testTrialBasis,trialOperator,
                                                                                                              sideIndex,false));
        // make copies b/c we can't fudge with return values from basisCache (const) - dimensions (Cell,Field - basis ordinal, Point)
        FieldContainer<double> testTrialValuesWeightedCopy = testTrialValuesTransformedWeighted;
        
        int numDofs2 = trialOrderPtr->getBasisCardinality(testTrialID,sideIndex); 
        for (int cellIndex=0; cellIndex<numCells; cellIndex++){
          for (int dofIndex=0; dofIndex<numDofs2; dofIndex++){
            for (int ptIndex=0; ptIndex<numPts; ptIndex++){
              testTrialValuesWeightedCopy(cellIndex, dofIndex, ptIndex) *= constrTestPair.second(cellIndex, ptIndex); // scale by constraint coeff
            }	   
          }
        }
        
        // loop thru pairs of trialIDs and constr coeffs
        for (map<int,FieldContainer<double> >::iterator constrIDIt = constrCoeffs.begin();
             constrIDIt !=constrCoeffs.end(); constrIDIt++) {
          pair<int,FieldContainer<double> > constrPair = *constrIDIt;
          int trialID = constrPair.first;
          
          // get basis to integrate
          BasisPtr trialBasis1 = trialOrderPtr->getBasis(trialID,sideIndex);
          // for trial: the value lives on the side, so we don't use the volume coords either:
          FieldContainer<double> trialValuesTransformed = *(basisCache->getTransformedValues(trialBasis1,trialOperator,sideIndex,false));
          // make copies b/c we can't fudge with return values from basisCache (const) - dimensions (Cell,Field - basis ordinal, Point)
          FieldContainer<double> trialValuesCopy = trialValuesTransformed;
          
          // transform trial values
          int numDofs1 = trialOrderPtr->getBasisCardinality(trialID,sideIndex); 
          for (int dofIndex=0; dofIndex<numDofs1; dofIndex++){
            for (int cellIndex=0; cellIndex<numCells; cellIndex++){
              for (int ptIndex=0; ptIndex<numPts; ptIndex++){
                trialValuesCopy(cellIndex, dofIndex, ptIndex) *= constrPair.second(cellIndex, ptIndex); // scale by constraint coeff
              }	   
            }
          }
          
          /////////////////////////////////////////////////////////////////////////////////////
          
          
          // integrate the transformed values, add them to the relevant trial/testTrialID dof combos
          FieldContainer<double> unweightedPenaltyMatrix(numCells,numDofs1,numDofs2);
          FunctionSpaceTools::integrate<double>(unweightedPenaltyMatrix,trialValuesCopy,testTrialValuesWeightedCopy,COMP_BLAS);
          
          for (int cellIndex=0; cellIndex<numCells; cellIndex++){
            for (int testDofIndex=0; testDofIndex<numDofs2; testDofIndex++){		
              int localTestDof = trialOrderPtr->getDofIndex(testTrialID, testDofIndex, sideIndex);
              for (int trialDofIndex=0; trialDofIndex<numDofs1; trialDofIndex++){		
                int localTrialDof = trialOrderPtr->getDofIndex(trialID, trialDofIndex, sideIndex);
                localStiffnessMatrix(cellIndex,localTrialDof,localTestDof) 
                                 += penaltyParameter*unweightedPenaltyMatrix(cellIndex,trialDofIndex,testDofIndex);
              }
            }
//.........这里部分代码省略.........

  void values(FieldContainer<double> &values, BasisCachePtr basisCache)
  {
    // sets values(_cellIndex,P,D)
    TEUCHOS_TEST_FOR_EXCEPTION(_cellIndex == -1, std::invalid_argument, "must call setCellIndex before calling values!");

//    cout << "_basisCoefficients:\n" << _basisCoefficients;

    BasisCachePtr spaceTimeBasisCache;
    if (basisCache->cellTopologyIsSpaceTime())
    {
      // then we require that the basisCache provided be a space-time basis cache
      SpaceTimeBasisCache* spaceTimeCache = dynamic_cast<SpaceTimeBasisCache*>(basisCache.get());
      TEUCHOS_TEST_FOR_EXCEPTION(!spaceTimeCache, std::invalid_argument, "space-time requires a SpaceTimeBasisCache");
      spaceTimeBasisCache = basisCache;
      basisCache = spaceTimeCache->getSpatialBasisCache();
    }

    int numDofs = _basis->getCardinality();
    int spaceDim = basisCache->getSpaceDim();

    bool basisIsVolumeBasis = (spaceDim == _basis->domainTopology()->getDimension());
    bool useCubPointsSideRefCell = basisIsVolumeBasis && basisCache->isSideCache();

    int numPoints = values.dimension(1);

    // check if we're taking a temporal derivative
    int component;
    Intrepid::EOperator relatedOp = BasisEvaluation::relatedOperator(_op, _basis->functionSpace(), spaceDim, component);
    if ((relatedOp == Intrepid::OPERATOR_GRAD) && (component==spaceDim)) {
      // then we are taking the temporal part of the Jacobian of the reference to curvilinear-reference space
      // based on our assumptions that curvilinearity is just in the spatial direction (and is orthogonally extruded in the
      // temporal direction), this is always the identity.
      for (int ptIndex=0; ptIndex<numPoints; ptIndex++)
      {
        for (int d=0; d<values.dimension(2); d++)
        {
          if (d < spaceDim)
            values(_cellIndex,ptIndex,d) = 0.0;
          else
            values(_cellIndex,ptIndex,d) = 1.0;
        }
      }
      return;
    }
    constFCPtr transformedValues = basisCache->getTransformedValues(_basis, _op, useCubPointsSideRefCell);

    // transformedValues has dimensions (C,F,P,[D,D])
    // therefore, the rank of the sum is transformedValues->rank() - 3
    int rank = transformedValues->rank() - 3;
    TEUCHOS_TEST_FOR_EXCEPTION(rank != values.rank()-2, std::invalid_argument, "values rank is incorrect.");


    int spaceTimeSideOrdinal = (spaceTimeBasisCache != Teuchos::null) ? spaceTimeBasisCache->getSideIndex() : -1;
    // I'm pretty sure much of this treatment of the time dimension could be simplified by taking advantage of SpaceTimeBasisCache::getTemporalBasisCache()...
    double t0 = -1, t1 = -1;
    if ((spaceTimeSideOrdinal != -1) && (!spaceTimeBasisCache->cellTopology()->sideIsSpatial(spaceTimeSideOrdinal)))
    {
      unsigned sideTime0 = spaceTimeBasisCache->cellTopology()->getTemporalSideOrdinal(0);
      unsigned sideTime1 = spaceTimeBasisCache->cellTopology()->getTemporalSideOrdinal(1);
      // get first node of each of the time-orthogonal sides, and use that to determine t0 and t1:
      unsigned spaceTimeNodeTime0 = spaceTimeBasisCache->cellTopology()->getNodeMap(spaceDim, sideTime0, 0);
      unsigned spaceTimeNodeTime1 = spaceTimeBasisCache->cellTopology()->getNodeMap(spaceDim, sideTime1, 0);
      t0 = spaceTimeBasisCache->getPhysicalCellNodes()(_cellIndex,spaceTimeNodeTime0,spaceDim);
      t1 = spaceTimeBasisCache->getPhysicalCellNodes()(_cellIndex,spaceTimeNodeTime1,spaceDim);
    }

    // initialize the values we're responsible for setting
    if (_op == OP_VALUE)
    {
      for (int ptIndex=0; ptIndex<numPoints; ptIndex++)
      {
        for (int d=0; d<values.dimension(2); d++)
        {
          if (d < spaceDim)
            values(_cellIndex,ptIndex,d) = 0.0;
          else if ((spaceTimeBasisCache != Teuchos::null) && (spaceTimeSideOrdinal == -1))
            values(_cellIndex,ptIndex,spaceDim) = spaceTimeBasisCache->getPhysicalCubaturePoints()(_cellIndex,ptIndex,spaceDim);
          else if ((spaceTimeBasisCache != Teuchos::null) && (spaceTimeSideOrdinal != -1))
          {
            if (spaceTimeBasisCache->cellTopology()->sideIsSpatial(spaceTimeSideOrdinal))
            {
              values(_cellIndex,ptIndex,spaceDim) = spaceTimeBasisCache->getPhysicalCubaturePoints()(_cellIndex,ptIndex,spaceDim-1);
            }
            else
            {
              double temporalPoint;
              unsigned temporalNode = spaceTimeBasisCache->cellTopology()->getTemporalComponentSideOrdinal(spaceTimeSideOrdinal);
              if (temporalNode==0)
                temporalPoint = t0;
              else
                temporalPoint = t1;
              values(_cellIndex,ptIndex,spaceDim) = temporalPoint;
            }
          }
        }
      }
    }
    else if ((_op == OP_DX) || (_op == OP_DY) || (_op == OP_DZ))
    {
      for (int ptIndex=0; ptIndex<numPoints; ptIndex++)
//.........这里部分代码省略.........