C++ BasisPtr::functionSpace方法代码示例

FCPtr BasisEvaluation::getTransformedValues(BasisPtr basis, Camellia::EOperator op,
                                            const FieldContainer<double> &referencePoints,
                                            int numCells,
                                            BasisCache* basisCache)
{
  Camellia::EFunctionSpace fs = basis->functionSpace();
  int spaceDim = referencePoints.dimension(1);
  int componentOfInterest;
  Intrepid::EOperator relatedOp;
  relatedOp = relatedOperator(op, fs, spaceDim, componentOfInterest);
  
  FCPtr referenceValues = getValues(basis,(Camellia::EOperator) relatedOp, referencePoints);
  return getTransformedValuesWithBasisValues(basis,op,referenceValues,numCells,basisCache);
}

void Projector::projectFunctionOntoBasisInterpolating(FieldContainer<double> &basisCoefficients, FunctionPtr fxn,
                                                      BasisPtr basis, BasisCachePtr domainBasisCache) {
  basisCoefficients.initialize(0);
  CellTopoPtr domainTopo = basis->domainTopology();
  unsigned domainDim = domainTopo->getDimension();
  
  IPPtr ip;
  
  bool traceVar = domainBasisCache->isSideCache();
  
  pair<IPPtr, VarPtr> ipVarPair = IP::standardInnerProductForFunctionSpace(basis->functionSpace(), traceVar, domainDim);
  ip = ipVarPair.first;
  VarPtr v = ipVarPair.second;
  
  IPPtr ip_l2 = Teuchos::rcp( new IP );
  ip_l2->addTerm(v);
  
  // for now, make all projections use L^2... (having some issues with gradients and cell Jacobians--I think we need the restriction of the cell Jacobian to the subcell, e.g., and it's not clear how to do that...)
  ip = ip_l2;
  
  FieldContainer<double> referenceDomainNodes(domainTopo->getVertexCount(),domainDim);
  CamelliaCellTools::refCellNodesForTopology(referenceDomainNodes, domainTopo);
  
  int basisCardinality = basis->getCardinality();
  
  set<int> allDofs;
  for (int i=0; i<basisCardinality; i++) {
    allDofs.insert(i);
  }
  
  for (int d=0; d<=domainDim; d++) {
    FunctionPtr projectionThusFar = NewBasisSumFunction::basisSumFunction(basis, basisCoefficients);
    FunctionPtr fxnToApproximate = fxn - projectionThusFar;
    int subcellCount = domainTopo->getSubcellCount(d);
    for (int subcord=0; subcord<subcellCount; subcord++) {
      set<int> subcellDofOrdinals = basis->dofOrdinalsForSubcell(d, subcord);
      if (subcellDofOrdinals.size() > 0) {
        FieldContainer<double> refCellPoints;
        FieldContainer<double> cubatureWeightsSubcell; // allows us to integrate over the fine subcell even when domain is higher-dimensioned
        if (d == 0) {
          refCellPoints.resize(1,domainDim);
          for (int d1=0; d1<domainDim; d1++) {
            refCellPoints(0,d1) = referenceDomainNodes(subcord,d1);
          }
          cubatureWeightsSubcell.resize(1);
          cubatureWeightsSubcell(0) = 1.0;
        } else {
          CellTopoPtr subcellTopo = domainTopo->getSubcell(d, subcord);
//          Teuchos::RCP<Cubature<double> > subcellCubature = cubFactory.create(subcellTopo, domainBasisCache->cubatureDegree());
          BasisCachePtr subcellCache = Teuchos::rcp( new BasisCache(subcellTopo, domainBasisCache->cubatureDegree(), false) );
          int numPoints = subcellCache->getRefCellPoints().dimension(0);
          refCellPoints.resize(numPoints,domainDim);
          cubatureWeightsSubcell = subcellCache->getCubatureWeights();
          
          if (d == domainDim) {
            refCellPoints = subcellCache->getRefCellPoints();
          } else {
            CamelliaCellTools::mapToReferenceSubcell(refCellPoints, subcellCache->getRefCellPoints(), d,
                                                     subcord, domainTopo);
          }
        }
        domainBasisCache->setRefCellPoints(refCellPoints, cubatureWeightsSubcell);
        IPPtr ipForProjection = (d==0) ? ip_l2 : ip; // just use values at vertices (ignore derivatives)
        set<int> dofsToSkip = allDofs;
        for (set<int>::iterator dofOrdinalIt=subcellDofOrdinals.begin(); dofOrdinalIt != subcellDofOrdinals.end(); dofOrdinalIt++) {
          dofsToSkip.erase(*dofOrdinalIt);
        }
        FieldContainer<double> newBasisCoefficients;
        projectFunctionOntoBasis(newBasisCoefficients, fxnToApproximate, basis, domainBasisCache, ipForProjection, v, dofsToSkip);
        for (int cellOrdinal=0; cellOrdinal<newBasisCoefficients.dimension(0); cellOrdinal++) {
          for (set<int>::iterator dofOrdinalIt=subcellDofOrdinals.begin(); dofOrdinalIt != subcellDofOrdinals.end(); dofOrdinalIt++) {
            basisCoefficients(cellOrdinal,*dofOrdinalIt) = newBasisCoefficients(cellOrdinal,*dofOrdinalIt);
          }
        }
      }
    }
  }
}

FCPtr BasisEvaluation::getValues(BasisPtr basis, Camellia::EOperator op,
                                 const FieldContainer<double> &refPoints)
{
  int numPoints = refPoints.dimension(0);
  int spaceDim = refPoints.dimension(1);  // points dimensions are (numPoints, spaceDim)
  int basisCardinality = basis->getCardinality();
  int componentOfInterest = -1;
  Camellia::EFunctionSpace fs = basis->functionSpace();
  Intrepid::EOperator relatedOp = relatedOperator(op, fs, spaceDim, componentOfInterest);

  int opCardinality = spaceDim; // for now, we assume basis values are in the same spaceDim as points (e.g. vector 1D has just 1 component)
  
  bool relatedOpIsDkOperator = (relatedOp >= OPERATOR_D1) && (relatedOp <= OPERATOR_D10);
  if (relatedOpIsDkOperator)
  {
    opCardinality = Intrepid::getDkCardinality(relatedOp, spaceDim);
  }
  
  if ((Camellia::EOperator)relatedOp != op)
  {
    // we can assume relatedResults has dimensions (numPoints,basisCardinality,spaceDimOut)
    FCPtr relatedResults = Teuchos::rcp(new FieldContainer<double>(basisCardinality,numPoints,opCardinality));
    basis->getValues(*(relatedResults.get()), refPoints, relatedOp);
    FCPtr result = getComponentOfInterest(relatedResults,op,fs,componentOfInterest);
    if ( result.get() == 0 )
    {
      result = relatedResults;
    }
    return result;
  }
  // if we get here, we should have a standard Intrepid operator, in which case we should
  // be able to: size a FieldContainer appropriately, and then call basis->getValues

  // But let's do just check that we have a standard Intrepid operator
  if ( (op >= Camellia::OP_X) || (op <  Camellia::OP_VALUE) )
  {
    TEUCHOS_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Unknown operator.");
  }

  // result dimensions should be either (numPoints,basisCardinality) or (numPoints,basisCardinality,spaceDimOut);
  Teuchos::Array<int> dimensions;
  dimensions.push_back(basisCardinality);
  dimensions.push_back(numPoints);
  int basisRank = basis->rangeRank();
  if (((basisRank == 0) && (op ==  Camellia::OP_VALUE)))
  {
    // scalar; leave as is
  }
  else if ((basisRank == 0) && relatedOpIsDkOperator)
  {
    dimensions.push_back(opCardinality);
  }
  else if ( ( ( basisRank == 1) && (op ==  Camellia::OP_VALUE) ) )
  {
    dimensions.push_back(opCardinality);
  }
  else if (
    ( ( basisRank == 0) && (op == Camellia::OP_GRAD) )
    ||
    ( ( basisRank == 0) && (op == Camellia::OP_CURL) ) )
  {
    dimensions.push_back(spaceDim);
  }
  else if ( (basis->rangeRank() == 1) && (op == Camellia::OP_GRAD) )
  {
    // grad of vector: a tensor
    dimensions.push_back(spaceDim);
    dimensions.push_back(opCardinality);
  }
  FCPtr result = Teuchos::rcp(new FieldContainer<double>(dimensions));
  basis->getValues(*(result.get()), refPoints, (Intrepid::EOperator)op);
  return result;
}

FCPtr BasisEvaluation::getTransformedValuesWithBasisValues(BasisPtr basis, Camellia::EOperator op,
                                                           constFCPtr referenceValues, int numCells,
                                                           BasisCache* basisCache)
{
  typedef FunctionSpaceTools fst;
  //  int numCells = cellJacobian.dimension(0);
  
  int spaceDim = basisCache->getSpaceDim(); // changed 3-21-16
  
  int componentOfInterest;
  Camellia::EFunctionSpace fs = basis->functionSpace();
  Intrepid::EOperator relatedOp = relatedOperator(op,fs,spaceDim, componentOfInterest);
  Teuchos::Array<int> dimensions;
  referenceValues->dimensions(dimensions);
  dimensions.insert(dimensions.begin(), numCells);
  Teuchos::RCP<FieldContainer<double> > transformedValues = Teuchos::rcp(new FieldContainer<double>(dimensions));
  bool vectorizedBasis = functionSpaceIsVectorized(fs);
  if (vectorizedBasis && (op ==  Camellia::OP_VALUE))
  {
    TEUCHOS_TEST_FOR_EXCEPTION( vectorizedBasis && (op ==  Camellia::OP_VALUE),
                               std::invalid_argument, "Vector HGRAD/HVOL with OP_VALUE not supported by getTransformedValuesWithBasisValues.  Please use getTransformedVectorValuesWithComponentBasisValues instead.");
  }
  switch(relatedOp)
  {
    case(Intrepid::OPERATOR_VALUE):
      switch(fs)
    {
      case Camellia::FUNCTION_SPACE_REAL_SCALAR:
        //          cout << "Reference values for FUNCTION_SPACE_REAL_SCALAR: " << *referenceValues;
      case Camellia::FUNCTION_SPACE_HGRAD:
      case Camellia::FUNCTION_SPACE_HGRAD_DISC:
        fst::HGRADtransformVALUE<double>(*transformedValues,*referenceValues);
        break;
      case Camellia::FUNCTION_SPACE_HCURL:
      case Camellia::FUNCTION_SPACE_HCURL_DISC:
        fst::HCURLtransformVALUE<double>(*transformedValues,basisCache->getJacobianInv(),*referenceValues);
        break;
      case Camellia::FUNCTION_SPACE_HDIV:
      case Camellia::FUNCTION_SPACE_HDIV_DISC:
      case Camellia::FUNCTION_SPACE_HDIV_FREE:
        fst::HDIVtransformVALUE<double>(*transformedValues,basisCache->getJacobian(),basisCache->getJacobianDet(),*referenceValues);
        break;
      case Camellia::FUNCTION_SPACE_HVOL:
      case Camellia::FUNCTION_SPACE_HVOL_SPACE_HGRAD_TIME:
      case Camellia::FUNCTION_SPACE_HGRAD_SPACE_HVOL_TIME:
        //        {
        //          static bool haveWarned = false;
        //          if (!haveWarned) {
        //            cout << "WARNING: for the moment, switching to the standard HVOLtransformVALUE method.\n";
        //            haveWarned = true;
        //          }
        //        }
        //          fst::HVOLtransformVALUE<double>(*transformedValues, cellJacobianDet, *referenceValues);
        // for the moment, use the fact that we know the HVOL basis is always an HGRAD basis:
        // (I think using the below amounts to solving for the HVOL variables scaled by Jacobian)
        fst::HGRADtransformVALUE<double>(*transformedValues,*referenceValues);
        break;
      default:
        TEUCHOS_TEST_FOR_EXCEPTION(true,std::invalid_argument, "unhandled transformation");
        break;
    }
      break;
    case(Intrepid::OPERATOR_GRAD):
    case(Intrepid::OPERATOR_D1):
      switch(fs)
    {
      case Camellia::FUNCTION_SPACE_HVOL:
      case Camellia::FUNCTION_SPACE_HGRAD:
      case Camellia::FUNCTION_SPACE_HGRAD_DISC:
        fst::HGRADtransformGRAD<double>(*transformedValues,basisCache->getJacobianInv(),*referenceValues);
        break;
      case Camellia::FUNCTION_SPACE_VECTOR_HVOL:
      case Camellia::FUNCTION_SPACE_VECTOR_HGRAD:
      case Camellia::FUNCTION_SPACE_VECTOR_HGRAD_DISC:
        // referenceValues has dimensions (F,P,D1,D2).  D1 is our component dimension, and D2 is the one that came from the gradient.
        // HGRADtransformGRAD expects (F,P,D) for input, and (C,F,P,D) for output.
        // If we split referenceValues into (F,P,D1=0,D2) and (F,P,D1=1,D2), then we can transform each of those, and then interleave the results…
      {
        // block off so we can create new stuff inside the switch case
        Teuchos::Array<int> dimensions;
        referenceValues->dimensions(dimensions);
        int numFields = dimensions[0];
        int numPoints = dimensions[1];
        int D1 = dimensions[dimensions.size()-2];
        int D2 = dimensions[dimensions.size()-1];
        dimensions[dimensions.size()-2] = D2; // put D2 in the D1 spot
        dimensions.pop_back(); // get rid of original D2
        FieldContainer<double> refValuesSlice(dimensions);
        dimensions.insert(dimensions.begin(),numCells);
        FieldContainer<double> transformedValuesSlice(dimensions);
        
        //          int numEntriesPerSlice = refValuesSlice.size();
        //          int numEntriesPerTransformedSlice = transformedValuesSlice.size();
        
        for (int compIndex1=0; compIndex1<D1; compIndex1++)
        {
          // could speed the following along by doing the enumeration arithmetic in place...
          for (int fieldIndex=0; fieldIndex<numFields; fieldIndex++)
          {
            for (int ptIndex=0; ptIndex<numPoints; ptIndex++)
//.........这里部分代码省略.........