QuadCnGradient.h

Go to the documentation of this file.

// Gradients of the arbitrary-order Bernstein-Bezier basis on the quadrilateral.
//
// Each 1-d Bernstein polynomial and its derivative are computed along each
// parametric axis (see EdgeCnBasis.h and EdgeCnGradient.h, including the
// dx/dr = 1/2 factor the change of variables introduces), then combined by the
// product rule.
WORKSPACE(RealT, rtmp, order[0] + 1);
WORKSPACE(RealT, drtmp, order[0] + 1);
WORKSPACE(RealT, stmp, order[1] + 1);
WORKSPACE(RealT, dstmp, order[1] + 1);
 
{
  RealT xx = 0.5 * (rr + 1.);
  RealT yy = 1. - xx;
  WORKSPACE(RealT, xpow, order[0] + 1);
  WORKSPACE(RealT, ypow, order[0] + 1);
  xpow[0] = 1.;
  for (int kk = 1; kk <= order[0]; ++kk)
  {
    xpow[kk] = xpow[kk - 1] * xx;
  }
  ypow[order[0]] = 1.;
  for (int kk = order[0] - 1; kk >= 0; --kk)
  {
    ypow[kk] = ypow[kk + 1] * yy;
  }
  RealT coeff = 1.;
  for (int kk = 0; kk <= order[0]; ++kk)
  {
    rtmp[kk] = coeff * xpow[kk] * ypow[kk];
    RealT ddx = 0.;
    if (kk > 0)
    {
      ddx += RealT(kk) * xpow[kk - 1] * ypow[kk];
    }
    if (kk < order[0])
    {
      ddx -= RealT(order[0] - kk) * xpow[kk] * ypow[kk + 1];
    }
    drtmp[kk] = 0.5 * coeff * ddx;
    coeff = coeff * RealT(order[0] - kk) / RealT(kk + 1);
  }
}
{
  RealT xx = 0.5 * (ss + 1.);
  RealT yy = 1. - xx;
  WORKSPACE(RealT, xpow, order[1] + 1);
  WORKSPACE(RealT, ypow, order[1] + 1);
  xpow[0] = 1.;
  for (int kk = 1; kk <= order[1]; ++kk)
  {
    xpow[kk] = xpow[kk - 1] * xx;
  }
  ypow[order[1]] = 1.;
  for (int kk = order[1] - 1; kk >= 0; --kk)
  {
    ypow[kk] = ypow[kk + 1] * yy;
  }
  RealT coeff = 1.;
  for (int kk = 0; kk <= order[1]; ++kk)
  {
    stmp[kk] = coeff * xpow[kk] * ypow[kk];
    RealT ddx = 0.;
    if (kk > 0)
    {
      ddx += RealT(kk) * xpow[kk - 1] * ypow[kk];
    }
    if (kk < order[1])
    {
      ddx -= RealT(order[1] - kk) * xpow[kk] * ypow[kk + 1];
    }
    dstmp[kk] = 0.5 * coeff * ddx;
    coeff = coeff * RealT(order[1] - kk) / RealT(kk + 1);
  }
}
 
// Compute the tensor product of basis derivatives. Like the fixed-order HGrad
// gradients, store a full 3-tuple per basis function and zero the axes the cell
// does not parameterize.
int idx = 0;
for (int pp = 0; pp <= order[1]; ++pp)
{
  for (int oo = 0; oo <= order[0]; ++oo)
  {
    basisGradient[idx++] = drtmp[oo] * stmp[pp]; // d/dr
    basisGradient[idx++] = rtmp[oo] * dstmp[pp]; // d/ds
    basisGradient[idx++] = 0.;
  }
}