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

void DofOrdering::copyLikeCoefficients( FieldContainer<double> &newValues, Teuchos::RCP<DofOrdering> oldDofOrdering,
                          const FieldContainer<double> &oldValues ) {
  // copy the coefficients for the bases that agree between the two DofOrderings
  // requires that "like" bases are actually pointers to the same memory location
  TEUCHOS_TEST_FOR_EXCEPTION( newValues.rank() != 1, std::invalid_argument, "newValues.rank() != 1");
  TEUCHOS_TEST_FOR_EXCEPTION( newValues.size() != totalDofs(), std::invalid_argument, "newValues.size() != totalDofs()");
  TEUCHOS_TEST_FOR_EXCEPTION( oldValues.rank() != 1, std::invalid_argument, "oldValues.rank() != 1");
  TEUCHOS_TEST_FOR_EXCEPTION( oldValues.size() != oldDofOrdering->totalDofs(), std::invalid_argument, "oldValues.size() != oldDofOrdering->totalDofs()");
  
  newValues.initialize(0.0);
  
  for (set<int>::iterator varIDIt = varIDs.begin(); varIDIt != varIDs.end(); varIDIt++) {
    int varID = *varIDIt;
    int numSides = getNumSidesForVarID(varID);
    if ( numSides == oldDofOrdering->getNumSidesForVarID(varID) ) {
      for (int sideIndex=0; sideIndex < numSides; sideIndex++) {
        BasisPtr basis = getBasis(varID,sideIndex);
        if (basis.get() == oldDofOrdering->getBasis(varID,sideIndex).get() ) {
          // bases alike: copy coefficients
          int cardinality = basis->getCardinality();
          for (int dofOrdinal=0; dofOrdinal < cardinality; dofOrdinal++) {
            int dofIndex = getDofIndex(varID,dofOrdinal,sideIndex);
            newValues(dofIndex) = oldValues( oldDofOrdering->getDofIndex(varID,dofOrdinal,sideIndex) );
          }
        }
      }
    }
  }
}

int DofOrdering::getDofIndex(int varID, int basisDofOrdinal, int sideIndex, int subSideIndex) { 
  TEUCHOS_TEST_FOR_EXCEPTION( ( _indexNeedsToBeRebuilt ),
                     std::invalid_argument,
                     "getDofIndex called when _indexNeedsToBeRebuilt = true.  Call rebuildIndex() first.");
  if (subSideIndex >= 0) {
    // then we've got a MultiBasis, and the basisDofOrdinal we have is *relative* to the subbasis
    BasisPtr basis = getBasis(varID,sideIndex);
    if ( ! BasisFactory::basisFactory()->isMultiBasis(basis) ) {
      TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "subSideIndex >= 0 for non-MultiBasis...");
    }
    MultiBasis<>* multiBasis = (MultiBasis<>*) basis.get();
    //cout << "(basisDofOrdinal, subSideIndex) : (" << basisDofOrdinal << ", " << subSideIndex << ") --> ";
    basisDofOrdinal = multiBasis->relativeToAbsoluteDofOrdinal(basisDofOrdinal,subSideIndex);
    //cout << basisDofOrdinal << endl;
  }
  
  pair<int,int> key = make_pair(varID, sideIndex);
  map< pair<int,int>, vector<int> >::iterator entryIt = indices.find(key);
  if ( entryIt != indices.end() ) {
    int dofIndex = ((*entryIt).second)[basisDofOrdinal];
    if ((dofIndex < 0) || (dofIndex >= _nextIndex)) {
      cout << "dofIndex out of bounds.\n";
      TEUCHOS_TEST_FOR_EXCEPTION( (dofIndex < 0) || (dofIndex >= _nextIndex), std::invalid_argument, "dofIndex out of bounds.");
    }
    return dofIndex;
  } else {
    cout << "No entry found for dofIndex\n";
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "No entry found for DofIndex.");
    return -1;
  }
}

int DofOrdering::getTotalBasisCardinality() { // total number of *distinct* basis functions
  int totalBasisCardinality = 0;
  set< ::Camellia::Basis<>* > basesCounted;
  for (map< pair<int, int>, BasisPtr >::iterator basisIt = bases.begin(); basisIt != bases.end(); basisIt++) {
    BasisPtr basis = basisIt->second;
    if (basesCounted.find(basis.get())==basesCounted.end()) {
      basesCounted.insert(basis.get());
      totalBasisCardinality += basis->getCardinality();
    }
  }
  return totalBasisCardinality;
}

VectorBasisPtr basisForTransformation(ElementTypePtr cellType)
{
  // int polyOrder = std::max(cellType->trialOrderPtr->maxBasisDegree(), cellType->testOrderPtr->maxBasisDegree());
  int polyOrder = cellType->trialOrderPtr->maxBasisDegree();

  CellTopoPtr cellTopo = cellType->cellTopoPtr;
  if (cellTopo->getTensorialDegree() > 0)
  {
    // for now, we assume that this means space-time.  (At some point, we may support tensor-product spatial cell topologies for
    // things like fast quadrature support, and this would need revisiting then.)
    // we also assume that the curvilinearity is purely spatial (and in fact, just 2D for now).
    cellTopo = CellTopology::cellTopology(cellTopo->getShardsTopology(), cellTopo->getTensorialDegree() - 1);
  }

  BasisPtr basis = BasisFactory::basisFactory()->getBasis(polyOrder, cellTopo, Camellia::FUNCTION_SPACE_VECTOR_HGRAD);
  VectorBasisPtr vectorBasis = Teuchos::rcp( (VectorizedBasis<> *)basis.get(),false); // dynamic cast would be better
  return vectorBasis;
}

bool VectorizedBasisTestSuite::testVectorizedBasis()
{
  bool success = true;

  string myName = "testVectorizedBasis";

  shards::CellTopology quad_4(shards::getCellTopologyData<shards::Quadrilateral<4> >() );

  int polyOrder = 3, numPoints = 5, spaceDim = 2;

  BasisPtr hgradBasis = BasisFactory::basisFactory()->getBasis(polyOrder, quad_4.getKey(), Camellia::FUNCTION_SPACE_HGRAD);

  // first test: make a single-component vector basis.  This should agree in every entry with the basis itself, but its field container will have one higher rank...
  VectorizedBasis<> oneComp(hgradBasis, 1);

  FieldContainer<double> linePoints(numPoints, spaceDim);
  for (int i=0; i<numPoints; i++)
  {
    for (int j=0; j<spaceDim; j++)
    {
      linePoints(i,j) = ((double)(i + j)) / (numPoints + spaceDim);
    }
  }

  FieldContainer<double> compValues(hgradBasis->getCardinality(),numPoints);
  hgradBasis->getValues(compValues, linePoints, Intrepid::OPERATOR_VALUE);

  FieldContainer<double> values(hgradBasis->getCardinality(),linePoints.dimension(0),1); // one component
  oneComp.getValues(values, linePoints, Intrepid::OPERATOR_VALUE);

  for (int i=0; i<compValues.size(); i++)
  {
    double diff = abs(values[i]-compValues[i]);
    if (diff != 0.0)
    {
      success = false;
      cout << myName << ": one-component vector basis doesn't produce same values as component basis." << endl;
      cout << "difference: " << diff << " in enumerated value " << i << endl;
      cout << "values:\n" << values;
      cout << "compValues:\n" << compValues;
      return success;
    }
  }

  vector< BasisPtr > twoComps;
  twoComps.push_back( Teuchos::rcp( new VectorizedBasis<>(hgradBasis, 2) ) );
  twoComps.push_back( BasisFactory::basisFactory()->getBasis( polyOrder,
                      quad_4.getKey(), Camellia::FUNCTION_SPACE_VECTOR_HGRAD) );


  vector< BasisPtr >::iterator twoCompIt;
  for (twoCompIt = twoComps.begin(); twoCompIt != twoComps.end(); twoCompIt++)
  {
    BasisPtr twoComp = *twoCompIt;

    int componentCardinality = hgradBasis->getCardinality();

    if (twoComp->getCardinality() != 2 * hgradBasis->getCardinality() )
    {
      success = false;
      cout << myName << ": two-component vector basis cardinality != one-component cardinality * 2." << endl;
      cout << "twoComp->getCardinality(): " << twoComp->getCardinality() << endl;
      cout << "oneComp->getCardinality(): " << oneComp.getCardinality() << endl;
    }

    values.resize(twoComp->getCardinality(),linePoints.dimension(0),2); // two components
    twoComp->getValues(values, linePoints, Intrepid::OPERATOR_VALUE);
    for (int basisIndex=0; basisIndex<twoComp->getCardinality(); basisIndex++)
    {
      for (int k=0; k<numPoints; k++)
      {
        double xValueExpected = (basisIndex < componentCardinality) ? compValues(basisIndex,k) : 0;
        double xValueActual = values(basisIndex,k,0);
        double yValueExpected = (basisIndex >= componentCardinality) ? compValues(basisIndex - componentCardinality,k) : 0;
        double yValueActual = values(basisIndex,k,1);
        if ( ( abs(xValueActual - xValueExpected) != 0) || ( abs(yValueActual - yValueExpected) != 0) )
        {
          success = false;
          cout << myName << ": expected differs from actual\n";
          cout << "component\n" << compValues;
          cout << "vector values:\n" << values;
          return success;
        }
      }
    }

    // test the mapping from oneComp dofOrdinal to twoComp:
    VectorizedBasis<>* twoCompAsVectorBasis = (VectorizedBasis<>  *) twoComp.get();

    for (int compDofOrdinal=0; compDofOrdinal<oneComp.getCardinality(); compDofOrdinal++)
    {
      int dofOrdinal_0 = twoCompAsVectorBasis->getDofOrdinalFromComponentDofOrdinal(compDofOrdinal, 0);
      int dofOrdinal_1 = twoCompAsVectorBasis->getDofOrdinalFromComponentDofOrdinal(compDofOrdinal, 1);
      // we expect the lists to be stacked (this is implicit in the test above)
      // dofOrdinal_0 we expect to be == compDofOrdinal
      // dofOrdinal_1 we expect to be == compDofOrdinal + oneComp.getCardinality()
      if (dofOrdinal_0 != compDofOrdinal)
      {
        success = false;
        cout << "getDofOrdinalFromComponentDofOrdinal() not returning expected value in first component.\n";
//.........这里部分代码省略.........

void ParametricSurface::basisWeightsForEdgeInterpolant(FieldContainer<double> &edgeInterpolationCoefficients,
    VectorBasisPtr basis,
    MeshPtr mesh, int cellID)
{
  vector< ParametricCurvePtr > curves = mesh->parametricEdgesForCell(cellID);
  Teuchos::RCP<TransfiniteInterpolatingSurface> exactSurface = Teuchos::rcp( new TransfiniteInterpolatingSurface(curves) );
  exactSurface->setNeglectVertices(false);

  int basisDegree = basis->getDegree();
  shards::CellTopology line_2(shards::getCellTopologyData<shards::Line<2> >() );
  BasisPtr basis1D = BasisFactory::basisFactory()->getBasis(basisDegree, line_2.getKey(),
                     Camellia::FUNCTION_SPACE_HGRAD);

  BasisPtr compBasis = basis->getComponentBasis();
  int numComponents = basis->getNumComponents();
  if (numComponents != 2)
  {
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Only 2D surfaces supported right now");
  }

  edgeInterpolationCoefficients.resize(basis->getCardinality());

  set<int> edgeNodeFieldIndices = BasisFactory::basisFactory()->sideFieldIndices(basis,true); // true: include vertex dofs

  FieldContainer<double> dofCoords(compBasis->getCardinality(),2);
  IntrepidBasisWrapper< double, Intrepid::FieldContainer<double> >* intrepidBasisWrapper = dynamic_cast< IntrepidBasisWrapper< double, Intrepid::FieldContainer<double> >* >(compBasis.get());
  if (!intrepidBasisWrapper)
  {
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "compBasis does not appear to be an instance of IntrepidBasisWrapper");
  }
  Basis_HGRAD_QUAD_Cn_FEM<double, Intrepid::FieldContainer<double> >* intrepidBasis = dynamic_cast< Basis_HGRAD_QUAD_Cn_FEM<double, Intrepid::FieldContainer<double> >* >(intrepidBasisWrapper->intrepidBasis().get());
  if (!intrepidBasis)
  {
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "IntrepidBasisWrapper does not appear to wrap Basis_HGRAD_QUAD_Cn_FEM.");
  }
  intrepidBasis->getDofCoords(dofCoords);

  int edgeDim = 1;
  int vertexDim = 0;

  // set vertex dofs:
  for (int vertexIndex=0; vertexIndex<curves.size(); vertexIndex++)
  {
    double x = exactSurface->vertices()[vertexIndex].first;
    double y = exactSurface->vertices()[vertexIndex].second;
    int compDofOrdinal = compBasis->getDofOrdinal(vertexDim, vertexIndex, 0);
    int basisDofOrdinal_x = basis->getDofOrdinalFromComponentDofOrdinal(compDofOrdinal, 0);
    int basisDofOrdinal_y = basis->getDofOrdinalFromComponentDofOrdinal(compDofOrdinal, 1);
    edgeInterpolationCoefficients[basisDofOrdinal_x] = x;
    edgeInterpolationCoefficients[basisDofOrdinal_y] = y;
  }

  for (int edgeIndex=0; edgeIndex<curves.size(); edgeIndex++)
  {
    bool edgeDofsFlipped = edgeIndex >= 2; // because Intrepid's ordering of dofs on the quad is not CCW but tensor-product, we need to flip for the opposite edges
    // (what makes things worse is that the vertex/edge numbering *is* CCW)
    if (curves.size() != 4)
    {
      cout << "WARNING: have not worked out the rule for flipping or not flipping edge dofs for anything but quads.\n";
    }
    double edgeLength = curves[edgeIndex]->linearLength();

    //    cout << "edgeIndex " << edgeIndex << endl;
    for (int comp=0; comp<numComponents; comp++)
    {
      FieldContainer<double> basisCoefficients_comp;
      bool useH1ForEdgeInterpolant = true; // an experiment
      curves[edgeIndex]->projectionBasedInterpolant(basisCoefficients_comp, basis1D, comp, edgeLength, useH1ForEdgeInterpolant);
      //      cout << "for edge " << edgeIndex << " and comp " << comp << ", projection-based interpolant dofs:\n";
      //      cout << basisCoefficients_comp;
      ////      cout << "basis dof coords:\n" << dofCoords;
      //      int basisDofOrdinal = basis->getDofOrdinalFromComponentDofOrdinal(v0_dofOrdinal_comp, comp);
      //      edgeInterpolationCoefficients[basisDofOrdinal] = basisCoefficients_comp[v0_dofOrdinal_1D];

      if (compBasis->getDegree() >= 2)   // then there are some "middle" nodes on the edge
      {
        // get the first dofOrdinal for the edge, so we can check the number of edge basis functions
        int firstEdgeDofOrdinal = compBasis->getDofOrdinal(edgeDim, edgeIndex, 0);

        //        cout << "first edge dofOrdinal: " << firstEdgeDofOrdinal << endl;

        int numEdgeDofs = compBasis->getDofTag(firstEdgeDofOrdinal)[3];
        if (numEdgeDofs != basis1D->getCardinality() - 2)
        {
          TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "numEdgeDofs does not match 1D basis cardinality");
        }
        for (int edgeDofOrdinal=0; edgeDofOrdinal<numEdgeDofs; edgeDofOrdinal++)
        {
          // determine the index into basisCoefficients_comp:
          int edgeDofOrdinalIn1DBasis = edgeDofsFlipped ? numEdgeDofs - 1 - edgeDofOrdinal : edgeDofOrdinal;
          int dofOrdinal1D = basis1D->getDofOrdinal(edgeDim, 0, edgeDofOrdinalIn1DBasis);
          // determine the ordinal of the edge dof in the component basis:
          int compDofOrdinal = compBasis->getDofOrdinal(edgeDim, edgeIndex, edgeDofOrdinal);
          // now, determine its ordinal in the vector basis
          int basisDofOrdinal = basis->getDofOrdinalFromComponentDofOrdinal(compDofOrdinal, comp);

          //          cout << "edge dof ordinal " << edgeDofOrdinal << " has basis weight " << basisCoefficients_comp[dofOrdinal1D] << " for component " << comp << endl;
          //          cout << "node on cell is at (" << dofCoords(compDofOrdinal,0) << ", " << dofCoords(compDofOrdinal,1) << ")\n";
          //          cout << "mapping to basisDofOrdinal " << basisDofOrdinal << endl;

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