[Question] Does anyone know logic of iso-parametric interpolation functions for different unstructured type of cells?

[kuldip-wagh]

We can see the r,s,t based interolation code for unstructured cell type types. Does anyone know logic/algorithm explanation behind it? The code looks something like this:

//----------------------------------------------------------------------------
// Compute iso-parametric interpolation functions
//
static inline void pyramidInterpolationFunctions(float pcoords[3], float sf[5])
{
  float rm, sm, tm;

  rm = 1.f - pcoords[0];
  sm = 1.f - pcoords[1];
  tm = 1.f - pcoords[2];

  sf[0] = rm * sm * tm;
  sf[1] = pcoords[0] * sm * tm;
  sf[2] = pcoords[0] * pcoords[1] * tm;
  sf[3] = rm * pcoords[1] * tm;
  sf[4] = pcoords[2];
}

//----------------------------------------------------------------------------
static inline void pyramidInterpolationDerivs(float pcoords[3],
                                              float derivs[15])
{
  // r-derivatives
  derivs[0] = -(pcoords[1] - 1.f) * (pcoords[2] - 1.f);
  derivs[1] = (pcoords[1] - 1.f) * (pcoords[2] - 1.f);
  derivs[2] = pcoords[1] - pcoords[1] * pcoords[2];
  derivs[3] = pcoords[1] * (pcoords[2] - 1.f);
  derivs[4] = 0.f;

  // s-derivatives
  derivs[5] = -(pcoords[0] - 1.f) * (pcoords[2] - 1.f);
  derivs[6] = pcoords[0] * (pcoords[2] - 1.f);
  derivs[7] = pcoords[0] - pcoords[0] * pcoords[2];
  derivs[8] = (pcoords[0] - 1.f) * (pcoords[2] - 1.f);
  derivs[9] = 0.f;

  // t-derivatives
  derivs[10] = -(pcoords[0] - 1.f) * (pcoords[1] - 1.f);
  derivs[11] = pcoords[0] * (pcoords[1] - 1.f);
  derivs[12] = -pcoords[0] * pcoords[1];
  derivs[13] = (pcoords[0] - 1.f) * pcoords[1];
  derivs[14] = 1.f;
}

static const uniform float PYRAMID_DIVERGED               = 1.e6;
static const uniform int PYRAMID_MAX_ITERATION            = 10;
static const uniform float PYRAMID_CONVERGED              = 1.e-04;
static const uniform float PYRAMID_OUTSIDE_CELL_TOLERANCE = 1.e-06;

static bool intersectAndSamplePyramid(const void *uniform userData,
                                      uniform uint64 id,
                                      uniform bool assumeInside,
                                      float &result,
                                      vec3f samplePos)
{
  const VKLUnstructuredVolume *uniform self =
      (const VKLUnstructuredVolume *uniform)userData;

  float pcoords[3] = {0.5, 0.5, 0.5};
  float derivs[15];
  float weights[5];

  // Get cell offset in index buffer
  const uniform uint64 cOffset             = getCellOffset(self, id);
  const uniform float determinantTolerance = self->iterativeTolerance[id];

  // Enter iteration loop
  bool converged = false;
  for (uniform int iteration = 0;
       !converged && (iteration < PYRAMID_MAX_ITERATION);
       iteration++) {
    unmasked
    {
      // Calculate element interpolation functions and derivatives
      pyramidInterpolationFunctions(pcoords, weights);
      pyramidInterpolationDerivs(pcoords, derivs);

      // Calculate newton functions
      vec3f fcol = make_vec3f(0.f, 0.f, 0.f);
      vec3f rcol = make_vec3f(0.f, 0.f, 0.f);
      vec3f scol = make_vec3f(0.f, 0.f, 0.f);
      vec3f tcol = make_vec3f(0.f, 0.f, 0.f);
      for (uniform int i = 0; i < 5; i++) {
        const uniform vec3f pt =
            get_vec3f(self->vertex, getVertexId(self, cOffset + i));
        fcol = fcol + pt * weights[i];
        rcol = rcol + pt * derivs[i];
        scol = scol + pt * derivs[i + 5];
        tcol = tcol + pt * derivs[i + 10];
      }

      fcol = fcol - samplePos;

      // Compute determinants and generate improvements
      const float d = det(make_LinearSpace3f(rcol, scol, tcol));
    }

    if (absf(d) < determinantTolerance) {
      return false;
    }

    const float d0 = det(make_LinearSpace3f(fcol, scol, tcol)) / d;
    const float d1 = det(make_LinearSpace3f(rcol, fcol, tcol)) / d;
    const float d2 = det(make_LinearSpace3f(rcol, scol, fcol)) / d;

    pcoords[0] = pcoords[0] - d0;
    pcoords[1] = pcoords[1] - d1;
    pcoords[2] = pcoords[2] - d2;

    // Convergence/divergence test - if neither, repeat
    if ((absf(d0) < PYRAMID_CONVERGED) & (absf(d1) < PYRAMID_CONVERGED) &
        (absf(d2) < PYRAMID_CONVERGED)) {
      converged = true;
    } else if ((absf(pcoords[0]) > PYRAMID_DIVERGED) |
               (absf(pcoords[1]) > PYRAMID_DIVERGED) |
               (absf(pcoords[2]) > PYRAMID_DIVERGED)) {
      return false;
    }
  }

  if (!converged) {
    return false;
  }

  const uniform float lowerlimit = 0.0 - PYRAMID_OUTSIDE_CELL_TOLERANCE;
  const uniform float upperlimit = 1.0 + PYRAMID_OUTSIDE_CELL_TOLERANCE;
  if (assumeInside || (pcoords[0] >= lowerlimit && pcoords[0] <= upperlimit &&
                       pcoords[1] >= lowerlimit && pcoords[1] <= upperlimit &&
                       pcoords[2] >= lowerlimit && pcoords[2] <= upperlimit)) {
    // Evaluation
    if (isValid(self->cellValue)) {
      result = get_float(self->cellValue, id);
    } else {
      float val = 0.f;
      for (uniform int i = 0; i < 5; i++) {
        val += weights[i] *
               get_float(self->vertexValue, getVertexId(self, cOffset + i));
      }
      result = val;
    }

    return true;
  }

  return false;