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

bool ParametricCurveTests::testGradientWrapper()
{
  bool success = true;

  // create an artificial function whose gradient is "interesting" and known
  FunctionPtr t1 = Function::xn(1);
  FunctionPtr t2 = Function::yn(1);
  FunctionPtr xt = t1 + t1 * t2;
  FunctionPtr yt = t2 + 2 * t1 * t2;
  FunctionPtr xt_dt1 = 1 + t2;
  FunctionPtr xt_dt2 = t1;
  FunctionPtr yt_dt1 = 2 * t2;
  FunctionPtr yt_dt2 = 1 + 2 * t1;

  FunctionPtr ft = Function::vectorize(xt, yt);
  FunctionPtr ft_dt1 = Function::vectorize(xt_dt1, yt_dt1);
  FunctionPtr ft_dt2 = Function::vectorize(xt_dt2, yt_dt2);

  FunctionPtr ft_gradt = Function::vectorize(ft_dt1, ft_dt2);

  // first test: confirm that on a parametric quad, the wrapped function agrees with the naked one
  int cubatureDegree = 5;
  BasisCachePtr parametricQuadCache = BasisCache::parametricQuadCache(cubatureDegree);
  FunctionPtr fx_gradx = ParametricCurve::parametricGradientWrapper(ft_gradt, true);

  double tol = 1e-14;
  if (! ft_gradt->equals(fx_gradx, parametricQuadCache))
  {
    success = false;
    cout << "on a parametric quad, the wrapped gradient doesn't agree with the naked one";
    reportFunctionValueDifferences(ft_gradt, fx_gradx, parametricQuadCache, tol);
  }
  if (! ft_gradt->equals(ft->grad(), parametricQuadCache))
  {
    success = false;
    cout << "on a parametric quad, manual gradient disagrees with automatic one (error in test construction, likely).";
    reportFunctionValueDifferences(ft_gradt, ft->grad(), parametricQuadCache, tol);
  }

  // on the quad domain defined by (0,0), (1,0), (2,3), (0,1),
  // some algebra shows that for x and y as functions of the parametric
  // coordinates, we have
  // x = t1 +     t1 * t2
  // y = t2 + 2 * t1 * t2
  // which gives the result that our original function f(t1,t2) = (t1 + t1 * t2, t2 + 2 * t1 * t2) =  (x, y)

  FunctionPtr x = Function::xn(1); // understood in physical space
  FunctionPtr y = Function::yn(1);
  FunctionPtr f1_xy = x;
  FunctionPtr f2_xy = y;
  FunctionPtr f_xy = Function::vectorize(f1_xy, f2_xy);

  // set up the quad domain
  FieldContainer<double> physicalCellNodes(1,4,2); // (C,P,D)
  physicalCellNodes(0,0,0) = 0;
  physicalCellNodes(0,0,1) = 0;

  physicalCellNodes(0,1,0) = 1;
  physicalCellNodes(0,1,1) = 0;

  physicalCellNodes(0,2,0) = 2;
  physicalCellNodes(0,2,1) = 3;

  physicalCellNodes(0,3,0) = 0;
  physicalCellNodes(0,3,1) = 1;

  // physical space BasisCache:
  shards::CellTopology quad_4(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
  BasisCachePtr basisCache = Teuchos::rcp( new BasisCache(physicalCellNodes, quad_4, cubatureDegree));

  // as a preliminary test, check that the Jacobian values and inverse values agree with our expectations
  // we expect the Jacobian to be:
  //        [ 1 + t2      t1 ]
  // 1/2 *  [                ]
  //        [ 2 * t2  1 + t2 ]
  // where (t1,t2) are parametric coordinates and the 1/2 comes from the transformation from reference
  // to parametric space
  int numCells = 1;
  int numPoints = basisCache->getRefCellPoints().dimension(0);
  int spaceDim = 2;
  FieldContainer<double> jacobianExpected(numCells,numPoints,spaceDim,spaceDim);
  FieldContainer<double> jacobianInvExpected(numCells,numPoints,spaceDim,spaceDim);
  // also check that the function we've chosen has the expected values
  // by first computing its gradient in parametric space and then dividing by 2 to account
  // for the transformation from reference to parametric space
  FieldContainer<double> fgrad_based_jacobian(numCells,numPoints,spaceDim,spaceDim);
  ft_gradt->values(fgrad_based_jacobian, parametricQuadCache);
  FieldContainer<double> parametricPoints = basisCache->computeParametricPoints();
  for (int ptIndex=0; ptIndex<numPoints; ptIndex++)
  {
    double t1 = parametricPoints(0,ptIndex,0);
    double t2 = parametricPoints(0,ptIndex,1);
    jacobianExpected(0,ptIndex,0,0) = 0.5 * (1 + t2);
    jacobianExpected(0,ptIndex,0,1) = 0.5 * (t1);
    jacobianExpected(0,ptIndex,1,0) = 0.5 * (2 * t2);
    jacobianExpected(0,ptIndex,1,1) = 0.5 * (1 + 2 * t1);
    jacobianInvExpected(0,ptIndex,0,0) = (2.0 / (1 + 2 * t1 + t2) ) * (1 + 2 * t1);
    jacobianInvExpected(0,ptIndex,0,1) = (2.0 / (1 + 2 * t1 + t2) ) * (- t1);
    jacobianInvExpected(0,ptIndex,1,0) = (2.0 / (1 + 2 * t1 + t2) ) * (- 2 * t2);
    jacobianInvExpected(0,ptIndex,1,1) = (2.0 / (1 + 2 * t1 + t2) ) * (1 + t2);
//.........这里部分代码省略.........

bool ParametricCurveTests::testTransfiniteInterpolant()
{
  bool success = true;

  // to begin, just let's do a simple quad mesh
  double width=4, height=3;
  double x0 = 0, y0 = 0;
  double x1 = width, y1 = 0;
  double x2 = width, y2 = height;
  double x3 = 0, y3 = height;

  vector< ParametricCurvePtr > edges(4);
  edges[0] = ParametricCurve::line(x0, y0, x1, y1);
  edges[1] = ParametricCurve::line(x1, y1, x2, y2);
  edges[2] = ParametricCurve::line(x2, y2, x3, y3);
  edges[3] = ParametricCurve::line(x3, y3, x0, y0);

  ParametricSurfacePtr transfiniteInterpolant = ParametricSurface::transfiniteInterpolant(edges);

  double x0_actual, y0_actual, x2_actual, y2_actual;
  transfiniteInterpolant->value(0, 0, x0_actual, y0_actual);
  transfiniteInterpolant->value(1, 1, x2_actual, y2_actual);

  double tol=1e-14;
  if ((abs(x0_actual-x0) > tol) || (abs(y0_actual-y0) > tol))
  {
    success = false;
    cout << "transfinite interpolant doesn't interpolate (x0,y0).\n";
  }
  if ((abs(x2_actual-x2) > tol) || (abs(y2_actual-y2) > tol))
  {
    success = false;
    cout << "transfinite interpolant doesn't interpolate (x2,y2).\n";
  }

  // the transfinite interpolant should be just (4t1, 3t2)
  FunctionPtr t1 = Function::xn(1);
  FunctionPtr t2 = Function::yn(1);
  FunctionPtr xPart = 4 * t1;
  FunctionPtr yPart = 3 * t2;
  FunctionPtr expected_tfi = Function::vectorize(xPart, yPart);

  int cubatureDegree = 4;
  BasisCachePtr parametricCache = BasisCache::parametricQuadCache(cubatureDegree);

  // a couple quick sanity checks:
  if (! expected_tfi->equals(expected_tfi, parametricCache))
  {
    success = false;
    cout << "ERROR in Function::equals(): vector Function does not equal itself.\n";
  }
  if (! transfiniteInterpolant->equals(transfiniteInterpolant, parametricCache))
  {
    success = false;
    cout << "Weird error: transfiniteInterpolant does not equal itself.\n";
  }

  // check that the transfiniteInterpolant's value() method matches values()
  {
    int numCells = 1;
    int numPoints = parametricCache->getRefCellPoints().dimension(0);
    int spaceDim = 2;
    FieldContainer<double> values(numCells,numPoints,spaceDim);
    transfiniteInterpolant->values(values, parametricCache);
    int cellIndex = 0;
    for (int ptIndex=0; ptIndex<numPoints; ptIndex++)
    {
      double t1 = parametricCache->getPhysicalCubaturePoints()(cellIndex,ptIndex,0);
      double t2 = parametricCache->getPhysicalCubaturePoints()(cellIndex,ptIndex,1);
      double x, y;
      transfiniteInterpolant->value(t1, t2, x, y);
      double x_expected = values(cellIndex,ptIndex,0);
      double y_expected = values(cellIndex,ptIndex,1);
      if (abs(x-x_expected) > tol)
      {
        cout << "transfinite interpolant value() does not match values()\n";
        success = false;
      }
      if (abs(y-y_expected) > tol)
      {
        cout << "transfinite interpolant value() does not match values()\n";
        success = false;
      }
    }
  }

  if (! expected_tfi->equals(transfiniteInterpolant, parametricCache, tol))
  {
    cout << "transfinite interpolant doesn't match expected.\n";
    success = false;
    int numCells = 1;
    int numPoints = parametricCache->getRefCellPoints().dimension(0);
    int spaceDim = 2;
    FieldContainer<double> values(numCells,numPoints,spaceDim);
    FieldContainer<double> expected_values(numCells,numPoints,spaceDim);
    expected_tfi->values(expected_values, parametricCache);
    transfiniteInterpolant->values(values, parametricCache);
    reportFunctionValueDifferences(parametricCache->getPhysicalCubaturePoints(), expected_values,
                                   values, tol);
  }
//.........这里部分代码省略.........

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