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

bool checkVertexOrdinalsQuad(BasisPtr basis, vector<int> &vertexOrdinals)
{
  // check that the given indices are exactly the vertex basis functions:
  // a) these are (1,0) or (0,1) at the corresponding vertices
  // b) others are (0,0) at the vertices

  int numVertices = 4;

  FieldContainer<double> refCellPoints(numVertices,2); // vertices, in order
  refCellPoints(0,0) = -1;
  refCellPoints(0,1) = -1;
  refCellPoints(1,0) =  1;
  refCellPoints(1,1) = -1;
  refCellPoints(2,0) =  1;
  refCellPoints(2,1) =  1;
  refCellPoints(3,0) = -1;
  refCellPoints(3,1) =  1;

  // assume scalar basis for now -- we'll throw an exception if not...
  FieldContainer<double> values(basis->getCardinality(), numVertices); // F, P
  basis->getValues(values, refCellPoints, OPERATOR_VALUE);

  double tol = 1e-14;
  for (int vertexIndex=0; vertexIndex<numVertices; vertexIndex++)
  {
    int vertexOrdinal = vertexOrdinals[vertexIndex];
    for (int fieldIndex=0; fieldIndex<basis->getCardinality(); fieldIndex++)
    {
      double value = values(fieldIndex,vertexIndex);
      if (fieldIndex==vertexOrdinal)
      {
        // expect non-zero
        if (value < tol)
        {
          return false;
        }
      }
      else
      {
        // expect zero
        if (value > tol)
        {
          return false;
        }
      }
    }
  }
  return true;
}

bool testBasisClassifications(BasisPtr basis) {
  bool success = true;
  
  CellTopoPtr cellTopo = basis->domainTopology();
  
  int numVertices = cellTopo->getVertexCount();
  int numEdges = cellTopo->getEdgeCount();
  
  int degree = basis->getDegree();
  
  // TODO: finish this
  vector<int> vertexOrdinals;
  for (int vertexIndex=0; vertexIndex < numVertices; vertexIndex++) {
    set<int> dofOrdinals = basis->dofOrdinalsForVertex(vertexIndex);
    if (dofOrdinals.size() == 0) TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "No dofOrdinal for vertex...");
    if (dofOrdinals.size() > 1) TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "More than one dofOrdinal per vertex...");
    vertexOrdinals.push_back(*(dofOrdinals.begin()));
  }
//  
//  if (! checkVertexOrdinalsQuad(basis, vertexOrdinals) ) {
//    success = false;
//    cout << "vertex dof ordinals don't match expected\n";
//    cout << "ordinals: ";
//    for (int vertexIndex=0; vertexIndex < numVertices; vertexIndex++) {
//      cout << vertexOrdinals[vertexIndex] << " ";
//    }
//    cout << endl;
//  }
  
  // get the points in reference space for each vertex
  FieldContainer<double> points;
  if (numVertices == 2) { // line
    points.resize(2,1);
    points(0,0) = -1;
    points(1,0) = 1;
  } else if (numVertices == 3) { // triangle
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "triangles not yet supported");
  } else if (numVertices == 4) { // quad
    points.resize(4,2);
    points(0,0) = -1;
    points(0,1) = -1;
    points(1,0) =  1;
    points(1,1) = -1;
    points(2,0) =  1;
    points(2,1) =  1;
    points(3,0) = -1;
    points(3,1) =  1;
  } else {
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "unsupported topology");
  }
  
  FieldContainer<double> vertexValues;
  if (basis->rangeRank() > 0) {
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "rank > 0 bases not yet supported");
  } else {
    vertexValues.resize(basis->getCardinality(),numVertices);
  }
  
  basis->getValues(vertexValues, points, Intrepid::OPERATOR_VALUE);
  
  // test that the points are correctly classified
  for (int fieldIndex=0; fieldIndex<basis->getCardinality(); fieldIndex++) {
    for (int ptIndex=0; ptIndex<numVertices; ptIndex++) {
      int dofOrdinalForPoint = vertexOrdinals[ptIndex];
      bool expectZero = (dofOrdinalForPoint != fieldIndex);
      if (expectZero) {
        if (vertexValues(fieldIndex,ptIndex) != 0) {
          success = false;
          cout << "Expected 0 for fieldIndex " << fieldIndex << " and ptIndex " << ptIndex;
          cout << ", but got " << vertexValues(fieldIndex,ptIndex) << endl;
        }
      } else {
        if (vertexValues(fieldIndex,ptIndex) == 0) {
          cout << "Expected nonzero for fieldIndex " << fieldIndex << " and ptIndex " << ptIndex << endl;
          success = false;
        }
      }
    }
  }
  
  if (!success) {
    cout << "Failed testBasisClassifications; vertexValues:\n" << vertexValues;
  }
  
  return success;
}

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

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