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

void ParametricSurface::values(FieldContainer<double> &values, BasisCachePtr basisCache)
{
  FieldContainer<double> parametricPoints = basisCache->computeParametricPoints();
  int numCells = parametricPoints.dimension(0);
  int numPoints = parametricPoints.dimension(1);

  for (int cellIndex=0; cellIndex<numCells; cellIndex++)
  {
    double x, y;
    for (int ptIndex=0; ptIndex<numPoints; ptIndex++)
    {
      double t1, t2;
      t1 = parametricPoints(cellIndex,ptIndex,0);
      t2 = parametricPoints(cellIndex,ptIndex,1);
      this->value(t1, t2, x, y);
      values(cellIndex,ptIndex,0) = x;
      values(cellIndex,ptIndex,1) = y;
    }
  }
}

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