Method InterpolateFunctions

Common/DataModel/vtkLagrangeTriangle.cxx:52–113

------------------------------------------------------------------------------

//------------------------------------------------------------------------------
void vtkLagrangeTriangle::InterpolateFunctions(const double pcoords[3], double* weights)
{
  // Adapted from P. Silvester, "High-Order Polynomial Triangular Finite
  // Elements for Potential Problems". Int. J. Engng Sci. Vol. 7, pp. 849-861.
  // Pergamon Press, 1969. The generic method is valid for all orders, but we
  // unroll the first two orders to reduce computational cost.

  double tau[3] = { pcoords[0], pcoords[1], 1. - pcoords[0] - pcoords[1] };

  vtkIdType n = this->GetOrder();

  if (n == 1)
  {
    // for the linear case, we simply return the parametric coordinates, rotated
    // into the parametric frame (e.g. barycentric tau_2 = parametric x).
    weights[0] = tau[2];
    weights[1] = tau[0];
    weights[2] = tau[1];
  }
  else if (n == 2)
  {
#ifdef SEVEN_POINT_TRIANGLE
    if (this->GetPoints()->GetNumberOfPoints() == 7)
    {
      double rs = tau[0] * tau[1];
      double rt = tau[0] * tau[2];
      double st = tau[1] * tau[2];
      double rst = rs * tau[2];
      weights[0] = tau[2] + 3.0 * rst - 2.0 * rt - 2.0 * st;
      weights[1] = tau[0] + 3.0 * rst - 2.0 * rt - 2.0 * rs;
      weights[2] = tau[1] + 3.0 * rst - 2.0 * rs - 2.0 * st;
      weights[3] = 4.0 * rt - 12.0 * rst;
      weights[4] = 4.0 * rs - 12.0 * rst;
      weights[5] = 4.0 * st - 12.0 * rst;
      weights[6] = 27.0 * rst;
      return;
    }
#endif
    weights[0] = tau[2] * (2.0 * tau[2] - 1.0);
    weights[1] = tau[0] * (2.0 * tau[0] - 1.0);
    weights[2] = tau[1] * (2.0 * tau[1] - 1.0);
    weights[3] = 4.0 * tau[0] * tau[2];
    weights[4] = 4.0 * tau[0] * tau[1];
    weights[5] = 4.0 * tau[1] * tau[2];
  }
  else
  {
    vtkIdType nPoints = this->GetPoints()->GetNumberOfPoints();

    for (vtkIdType idx = 0; idx < nPoints; idx++)
    {
      weights[idx] = 1.;
      vtkIdType lambda[3];
      this->ToBarycentricIndex(idx, lambda);

      for (vtkIdType dim = 0; dim < 3; dim++)
      {
        weights[idx] *= Eta(n, lambda[dim], tau[dim]);