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

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 ExactSolution::solutionValues(FieldContainer<double> &values, int trialID, BasisCachePtr basisCache) {
  if (_exactFunctions.find(trialID) != _exactFunctions.end() ) {
    _exactFunctions[trialID]->values(values,basisCache);
    return;
  }
  
  // TODO: change ExactSolution::solutionValues (below) to take a *const* points FieldContainer, to avoid this copy:
  FieldContainer<double> points = basisCache->getPhysicalCubaturePoints();
  if (basisCache->getSideIndex() >= 0) {
    FieldContainer<double> unitNormals = basisCache->getSideNormals();
    this->solutionValues(values,trialID,points,unitNormals);
  } else {
    this->solutionValues(values,trialID,points);
  }
}

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