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

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

void Projector::projectFunctionOntoBasis(FieldContainer<double> &basisCoefficients, FunctionPtr fxn, 
                                         BasisPtr basis, BasisCachePtr basisCache, IPPtr ip, VarPtr v,
                                         set<int> fieldIndicesToSkip) {
  CellTopoPtr cellTopo = basis->domainTopology();
  DofOrderingPtr dofOrderPtr = Teuchos::rcp(new DofOrdering());
  
  if (! fxn.get()) {
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "fxn cannot be null!");
  }
  
  int cardinality = basis->getCardinality();
  int numCells = basisCache->getPhysicalCubaturePoints().dimension(0);
  int numDofs = cardinality - fieldIndicesToSkip.size();
  if (numDofs==0) {
    // we're skipping all the fields, so just initialize basisCoefficients to 0 and return
    basisCoefficients.resize(numCells,cardinality);
    basisCoefficients.initialize(0);
    return;
  }
  
  FieldContainer<double> gramMatrix(numCells,cardinality,cardinality);
  FieldContainer<double> ipVector(numCells,cardinality);

  // fake a DofOrdering
  DofOrderingPtr dofOrdering = Teuchos::rcp( new DofOrdering );
  if (! basisCache->isSideCache()) {
    dofOrdering->addEntry(v->ID(), basis, v->rank());
  } else {
    dofOrdering->addEntry(v->ID(), basis, v->rank(), basisCache->getSideIndex());
  }
  
  ip->computeInnerProductMatrix(gramMatrix, dofOrdering, basisCache);
  ip->computeInnerProductVector(ipVector, v, fxn, dofOrdering, basisCache);
  
//  cout << "physical points for projection:\n" << basisCache->getPhysicalCubaturePoints();
//  cout << "gramMatrix:\n" << gramMatrix;
//  cout << "ipVector:\n" << ipVector;
  
  map<int,int> oldToNewIndices;
  if (fieldIndicesToSkip.size() > 0) {
    // the code to do with fieldIndicesToSkip might not be terribly efficient...
    // (but it's not likely to be called too frequently)
    int i_indices_skipped = 0;
    for (int i=0; i<cardinality; i++) {
      int new_index;
      if (fieldIndicesToSkip.find(i) != fieldIndicesToSkip.end()) {
        i_indices_skipped++;
        new_index = -1;
      } else {
        new_index = i - i_indices_skipped;
      }
      oldToNewIndices[i] = new_index;
    }
    
    FieldContainer<double> gramMatrixFiltered(numCells,numDofs,numDofs);
    FieldContainer<double> ipVectorFiltered(numCells,numDofs);
    // now filter out the values that we're to skip
    
    for (int cellIndex=0; cellIndex<numCells; cellIndex++) {
      for (int i=0; i<cardinality; i++) {
        int i_filtered = oldToNewIndices[i];
        if (i_filtered == -1) {
          continue;
        }
        ipVectorFiltered(cellIndex,i_filtered) = ipVector(cellIndex,i);
        
        for (int j=0; j<cardinality; j++) {
          int j_filtered = oldToNewIndices[j];
          if (j_filtered == -1) {
            continue;
          }
          gramMatrixFiltered(cellIndex,i_filtered,j_filtered) = gramMatrix(cellIndex,i,j);
        }
      }
    }
//    cout << "gramMatrixFiltered:\n" << gramMatrixFiltered;
//    cout << "ipVectorFiltered:\n" << ipVectorFiltered;
    gramMatrix = gramMatrixFiltered;
    ipVector = ipVectorFiltered;
  }
  
  for (int cellIndex=0; cellIndex<numCells; cellIndex++){
    
    // TODO: rewrite to take advantage of SerialDenseWrapper...
    Epetra_SerialDenseSolver solver;
    
    Epetra_SerialDenseMatrix A(Copy,
                               &gramMatrix(cellIndex,0,0),
                               gramMatrix.dimension(2), 
                               gramMatrix.dimension(2),  
                               gramMatrix.dimension(1)); // stride -- fc stores in row-major order (a.o.t. SDM)
    
    Epetra_SerialDenseVector b(Copy,
                               &ipVector(cellIndex,0),
                               ipVector.dimension(1));
    
    Epetra_SerialDenseVector x(gramMatrix.dimension(1));
    
    solver.SetMatrix(A);
    int info = solver.SetVectors(x,b);
//.........这里部分代码省略.........

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

//.........这里部分代码省略.........
        if (printTimings) {
            cout << "numCells: " << numCells << endl;
            cout << "numTestDofs: " << numTestDofs << endl;
        }

        FieldContainer<double> rhsVectorTest(numCells,testOrderInterior->totalDofs()); // rhsVector is zero for the "trace" test dofs
        {
            // project the load f onto the space of interior test dofs.
            LinearTermPtr f = rhs->linearTerm();
            set<int> testIDs = f->varIDs();
            for (int eqn=0; eqn<numEquations; eqn++) {
                VarPtr v = _virtualTerms.getFieldTestVars()[eqn];

                if (testIDs.find(v->ID()) != testIDs.end()) {
                    BasisPtr testInteriorBasis = testOrderInterior->getBasis(v->ID());
                    FieldContainer<double> fValues(numCells,testInteriorBasis->getCardinality());
//          DofOrderingPtr oneVarOrderingTest = Teuchos::rcp(new DofOrdering(testInteriorBasis->domainTopology()));
                    DofOrderingPtr oneVarOrderingTest = Teuchos::rcp(new DofOrdering);
                    oneVarOrderingTest->addEntry(v->ID(), testInteriorBasis, testInteriorBasis->rangeRank());

                    LinearTermPtr f_v = Teuchos::rcp( new LinearTerm );
                    typedef pair< FunctionPtr, VarPtr > LinearSummand;
                    vector<LinearSummand> summands = f->summands();
                    for (int i=0; i<summands.size(); i++) {
                        FunctionPtr f = summands[i].first;
                        if (v->ID() == summands[i].second->ID()) {
                            f_v->addTerm(f * v);
                            f_v->integrate(fValues, oneVarOrderingTest, basisCache);
                        }
                    }

                    LinearTermPtr v_lt = 1.0 * v;
                    FieldContainer<double> l2(numCells,testInteriorBasis->getCardinality(),testInteriorBasis->getCardinality());
                    v_lt->integrate(l2,oneVarOrderingTest,v_lt,oneVarOrderingTest,basisCache,basisCache->isSideCache());

                    Teuchos::Array<int> testTestDim(2), testOneDim(2);
                    testTestDim[0] = testInteriorBasis->getCardinality();
                    testTestDim[1] = testInteriorBasis->getCardinality();
                    testOneDim[0] = testInteriorBasis->getCardinality();
                    testOneDim[1] = 1;
                    FieldContainer<double> projection(testOneDim);
                    for (int cellOrdinal=0; cellOrdinal<numCells; cellOrdinal++) {
                        FieldContainer<double> l2cell(testTestDim,&l2(cellOrdinal,0,0));
                        FieldContainer<double> f_cell(testOneDim,&fValues(cellOrdinal,0));

                        SerialDenseWrapper::solveSystemUsingQR(projection, l2cell, f_cell);

                        // rows in projection correspond to Ae_i, columns to the e_j.  I.e. projection coefficients for e_i are found in the ith row
                        for (int basisOrdinal_j=0; basisOrdinal_j<projection.dimension(0); basisOrdinal_j++) {
                            int testIndex = testOrderInterior->getDofIndex(v->ID(), basisOrdinal_j);
                            rhsVectorTest(cellOrdinal,testIndex) = projection(basisOrdinal_j,0);
                        }
                    }
                }
            }
        }

        // project strong operator applied to field terms, and use this to populate the top left portion of stiffness matrix:
        {
            FieldContainer<double> trialFieldTestInterior(numCells, fieldOrder->totalDofs(), testOrderInterior->totalDofs());
            for (int eqn=0; eqn<numEquations; eqn++) {
                LinearTermPtr Au = _virtualTerms.getFieldOperators()[eqn];
                VarPtr v = _virtualTerms.getFieldTestVars()[eqn];
                set<int> fieldIDs = Au->varIDs();
                for (set<int>::iterator fieldIt = fieldIDs.begin(); fieldIt != fieldIDs.end(); fieldIt++) {
                    int fieldID = *fieldIt;