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

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

  shards::CellTopology quad_4(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
  int numVertices = quad_4.getVertexCount();
  int numComponents = quad_4.getDimension();
  int vertexDim = 0;

  for (int polyOrder = 1; polyOrder<10; polyOrder++)
  {
    BasisPtr hgradBasis =  BasisFactory::basisFactory()->getConformingBasis(polyOrder,
                           quad_4.getKey(),
                           Camellia::FUNCTION_SPACE_HGRAD);
    BasisPtr vectorHGradBasis = BasisFactory::basisFactory()->getConformingBasis( polyOrder,
                                quad_4.getKey(),
                                Camellia::FUNCTION_SPACE_VECTOR_HGRAD);
    vector<int> vertexNodeFieldIndices;
    for (int vertexIndex=0; vertexIndex<numVertices; vertexIndex++)
    {
      for (int comp=0; comp<numComponents; comp++)
      {
        int vertexNodeFieldIndex = vectorHGradBasis->getDofOrdinal(vertexDim, vertexIndex, comp);
        vertexNodeFieldIndices.push_back(vertexNodeFieldIndex);
//        cout << "vertexNodeFieldIndex for vertex index " << vertexIndex << ", comp " << comp;
//        cout << " = " << vertexNodeFieldIndex << endl;
      }
    }
    if (!checkVertexNodalIndicesQuad(vectorHGradBasis, vertexNodeFieldIndices) )
    {
      success = false;
      cout << "testVectorizedBasisTags: Vertex tags for vectorized HGRAD basis";
      cout << " of order " << polyOrder << " are incorrect.\n";
    }
  }

  return success;
}

//.........这里部分代码省略.........
  shards::CellTopology line_2(shards::getCellTopologyData<shards::Line<2> >() );
  BasisPtr quadraticBasis = BasisFactory::basisFactory()->getBasis(2, line_2.getKey(), Camellia::FUNCTION_SPACE_HGRAD);

  // figure out what the weights for the quadratic "middle node" basis function should be:
  double expected_H1_weight_x, expected_H1_weight_y;
  double expected_L2_weight_x, expected_L2_weight_y;

  FunctionPtr middleBasis;
  {
    FunctionPtr t = Function::xn(1);
    middleBasis = 4 * t * (1-t);
  }

  double middleBasisL2_squared = (middleBasis*middleBasis)->integrate(basisCache);
  double middleBasisH1_squared = ( middleBasis->dx() * middleBasis->dx() )->integrate(basisCache) + middleBasisL2_squared;

  ParametricCurvePtr circularArcBubble = ParametricCurve::bubble(circularArc);

  FunctionPtr bubble_x = circularArcBubble->x();
  FunctionPtr bubble_y = circularArcBubble->y();

  double x_against_middle_L2 = (bubble_x * middleBasis)->integrate(basisCache);
  double x_against_middle_H1 = (bubble_x->dx() * middleBasis->dx())->integrate(basisCache) + x_against_middle_L2;

  double y_against_middle_L2 = (bubble_y * middleBasis)->integrate(basisCache);
  double y_against_middle_H1 = (bubble_y->dx() * middleBasis->dx())->integrate(basisCache) + y_against_middle_L2;

  expected_L2_weight_x = x_against_middle_L2 / middleBasisL2_squared;
  expected_H1_weight_x = x_against_middle_H1 / middleBasisH1_squared;

  expected_L2_weight_y = y_against_middle_L2 / middleBasisL2_squared;
  expected_H1_weight_y = y_against_middle_H1 / middleBasisH1_squared;

  int middleBasisOrdinal = quadraticBasis->getDofOrdinal(1,0,0);

  FieldContainer<double> basisCoefficients_x, basisCoefficients_y;
  bool useH1 = false; // just trying to diagnose whether the issue is in derivatives or values (most likely derivatives)
  double lengthScale = 1.0;
  circularArcBubble->projectionBasedInterpolant(basisCoefficients_x, quadraticBasis, 0, lengthScale, useH1);
  circularArcBubble->projectionBasedInterpolant(basisCoefficients_y, quadraticBasis, 1, lengthScale, useH1);

  double weightError_x = abs(expected_L2_weight_x-basisCoefficients_x[middleBasisOrdinal]);
  double weightError_y = abs(expected_L2_weight_y-basisCoefficients_y[middleBasisOrdinal]);

  if (weightError_x > tol)
  {
    success = false;
    cout << "testCircularArc(): L2 projection doesn't match expected basis weight in x.\n";
    cout << "expected " << expected_L2_weight_x << ", was " << basisCoefficients_x[middleBasisOrdinal] << endl;
  }
  if (weightError_y > tol)
  {
    success = false;
    cout << "testCircularArc(): L2 projection doesn't match expected basis weight in y.\n";
    cout << "expected " << expected_L2_weight_y << ", was " << basisCoefficients_y[middleBasisOrdinal] << endl;
  }

  useH1 = true;
  circularArcBubble->projectionBasedInterpolant(basisCoefficients_x, quadraticBasis, 0, lengthScale, useH1);
  circularArcBubble->projectionBasedInterpolant(basisCoefficients_y, quadraticBasis, 1, lengthScale, useH1);

  weightError_x = abs(expected_H1_weight_x-basisCoefficients_x[middleBasisOrdinal]);
  weightError_y = abs(expected_H1_weight_y-basisCoefficients_y[middleBasisOrdinal]);

  if (weightError_x > tol)
  {

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;

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