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

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

    fgrad_based_jacobian(0,ptIndex,0,0) /= 2.0;
    fgrad_based_jacobian(0,ptIndex,0,1) /= 2.0;
    fgrad_based_jacobian(0,ptIndex,1,0) /= 2.0;
    fgrad_based_jacobian(0,ptIndex,1,1) /= 2.0;
  }
  FieldContainer<double> jacobian = basisCache->getJacobian();
  FieldContainer<double> jacobianInv = basisCache->getJacobianInv();
  double maxDiff = 0;
  if (! fcsAgree(jacobianExpected, jacobian, tol, maxDiff))
  {
    success = false;
    cout << "Jacobian expected does not match actual.\n";
    reportFunctionValueDifferences(parametricPoints, jacobian, jacobianExpected, tol);
  }
  if (! fcsAgree(jacobianInvExpected, jacobianInv, tol, maxDiff))
  {
    success = false;
    cout << "Jacobian inverse expected does not match actual.\n";
    reportFunctionValueDifferences(parametricPoints, jacobianInv, jacobianInvExpected, tol);
  }
  if (! fcsAgree(fgrad_based_jacobian, jacobianExpected, tol, maxDiff))
  {
    success = false;
    cout << "Jacobian from fgrad does not agree with the transformation jacobian (problem with test?).\n";
    reportFunctionValueDifferences(parametricPoints, fgrad_based_jacobian, jacobianExpected, tol);
  }

  // test that the gradient values agree
  if (! fx_gradx->equals(f_xy->grad(), basisCache))
  {
    success = false;
    cout << "wrapped gradient does not agree with analytically transformed function.\n";
    reportFunctionValueDifferences(fx_gradx, f_xy->grad(), basisCache, tol);
  }

  // finally, although this isn't really the right place for this, it is convenient here
  // to test the TFI for the "mesh" we were concerned with above.
  int H1Order = 5;
  BFPtr bf = VGPStokesFormulation(1.0).bf();
  physicalCellNodes.resize(4,2);
  int horizontalElements = 1, verticalElements = 1;
  MeshPtr mesh = MeshFactory::quadMesh(bf, H1Order, physicalCellNodes);

  int cellID = 0;
  vector< ParametricCurvePtr > edges = mesh->parametricEdgesForCell(cellID);
  ParametricSurfacePtr tfi = ParametricSurface::transfiniteInterpolant(edges);

  double v2[2];
  edges[2]->value(0, v2[0], v2[1]);
//  cout << "v2 = (" << v2[0] << ", " << v2[1] << ")\n";

  if ( ! tfi->equals(f_xy, basisCache) )
  {
    success = false;
    cout << "TFI does not agree with analytically constructed transformation function.\n";
    reportFunctionValueDifferences(tfi, f_xy, basisCache, tol);
  }
  if ( ! tfi->grad()->equals(f_xy->grad(), basisCache) )
  {
    success = false;
    cout << "TFI does not agree with analytically constructed transformation function.\n";
    reportFunctionValueDifferences(tfi->grad(), f_xy->grad(), basisCache, tol);
  }

  return success;
}