Actual source code: dmmbfem.cxx

      1: #include <petscconf.h>
      2: #include <petscdt.h>
      3: #include <petsc/private/dmmbimpl.h>

      5: /* Utility functions */
      6: static inline PetscReal DMatrix_Determinant_2x2_Internal(const PetscReal inmat[2 * 2])
      7: {
      8:   return inmat[0] * inmat[3] - inmat[1] * inmat[2];
      9: }

     11: static inline PetscErrorCode DMatrix_Invert_2x2_Internal(const PetscReal *inmat, PetscReal *outmat, PetscReal *determinant)
     12: {
     13:   PetscReal det = DMatrix_Determinant_2x2_Internal(inmat);
     14:   if (outmat) {
     15:     outmat[0] = inmat[3] / det;
     16:     outmat[1] = -inmat[1] / det;
     17:     outmat[2] = -inmat[2] / det;
     18:     outmat[3] = inmat[0] / det;
     19:   }
     20:   if (determinant) *determinant = det;
     21:   return PETSC_SUCCESS;
     22: }

     24: static inline PetscReal DMatrix_Determinant_3x3_Internal(const PetscReal inmat[3 * 3])
     25: {
     26:   return inmat[0] * (inmat[8] * inmat[4] - inmat[7] * inmat[5]) - inmat[3] * (inmat[8] * inmat[1] - inmat[7] * inmat[2]) + inmat[6] * (inmat[5] * inmat[1] - inmat[4] * inmat[2]);
     27: }

     29: static inline PetscErrorCode DMatrix_Invert_3x3_Internal(const PetscReal *inmat, PetscReal *outmat, PetscReal *determinant)
     30: {
     31:   PetscReal det = DMatrix_Determinant_3x3_Internal(inmat);
     32:   if (outmat) {
     33:     outmat[0] = (inmat[8] * inmat[4] - inmat[7] * inmat[5]) / det;
     34:     outmat[1] = -(inmat[8] * inmat[1] - inmat[7] * inmat[2]) / det;
     35:     outmat[2] = (inmat[5] * inmat[1] - inmat[4] * inmat[2]) / det;
     36:     outmat[3] = -(inmat[8] * inmat[3] - inmat[6] * inmat[5]) / det;
     37:     outmat[4] = (inmat[8] * inmat[0] - inmat[6] * inmat[2]) / det;
     38:     outmat[5] = -(inmat[5] * inmat[0] - inmat[3] * inmat[2]) / det;
     39:     outmat[6] = (inmat[7] * inmat[3] - inmat[6] * inmat[4]) / det;
     40:     outmat[7] = -(inmat[7] * inmat[0] - inmat[6] * inmat[1]) / det;
     41:     outmat[8] = (inmat[4] * inmat[0] - inmat[3] * inmat[1]) / det;
     42:   }
     43:   if (determinant) *determinant = det;
     44:   return PETSC_SUCCESS;
     45: }

     47: static inline PetscReal DMatrix_Determinant_4x4_Internal(PetscReal inmat[4 * 4])
     48: {
     49:   return inmat[0 + 0 * 4] * (inmat[1 + 1 * 4] * (inmat[2 + 2 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 2 * 4]) - inmat[1 + 2 * 4] * (inmat[2 + 1 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 1 * 4]) + inmat[1 + 3 * 4] * (inmat[2 + 1 * 4] * inmat[3 + 2 * 4] - inmat[2 + 2 * 4] * inmat[3 + 1 * 4])) - inmat[0 + 1 * 4] * (inmat[1 + 0 * 4] * (inmat[2 + 2 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 2 * 4]) - inmat[1 + 2 * 4] * (inmat[2 + 0 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 0 * 4]) + inmat[1 + 3 * 4] * (inmat[2 + 0 * 4] * inmat[3 + 2 * 4] - inmat[2 + 2 * 4] * inmat[3 + 0 * 4])) + inmat[0 + 2 * 4] * (inmat[1 + 0 * 4] * (inmat[2 + 1 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 1 * 4]) - inmat[1 + 1 * 4] * (inmat[2 + 0 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 0 * 4]) + inmat[1 + 3 * 4] * (inmat[2 + 0 * 4] * inmat[3 + 1 * 4] - inmat[2 + 1 * 4] * inmat[3 + 0 * 4])) - inmat[0 + 3 * 4] * (inmat[1 + 0 * 4] * (inmat[2 + 1 * 4] * inmat[3 + 2 * 4] - inmat[2 + 2 * 4] * inmat[3 + 1 * 4]) - inmat[1 + 1 * 4] * (inmat[2 + 0 * 4] * inmat[3 + 2 * 4] - inmat[2 + 2 * 4] * inmat[3 + 0 * 4]) + inmat[1 + 2 * 4] * (inmat[2 + 0 * 4] * inmat[3 + 1 * 4] - inmat[2 + 1 * 4] * inmat[3 + 0 * 4]));
     50: }

     52: static inline PETSC_UNUSED PetscErrorCode DMatrix_Invert_4x4_Internal(PetscReal *inmat, PetscReal *outmat, PetscScalar *determinant)
     53: {
     54:   PetscReal det = DMatrix_Determinant_4x4_Internal(inmat);
     55:   if (outmat) {
     56:     outmat[0]  = (inmat[5] * inmat[10] * inmat[15] + inmat[6] * inmat[11] * inmat[13] + inmat[7] * inmat[9] * inmat[14] - inmat[5] * inmat[11] * inmat[14] - inmat[6] * inmat[9] * inmat[15] - inmat[7] * inmat[10] * inmat[13]) / det;
     57:     outmat[1]  = (inmat[1] * inmat[11] * inmat[14] + inmat[2] * inmat[9] * inmat[15] + inmat[3] * inmat[10] * inmat[13] - inmat[1] * inmat[10] * inmat[15] - inmat[2] * inmat[11] * inmat[13] - inmat[3] * inmat[9] * inmat[14]) / det;
     58:     outmat[2]  = (inmat[1] * inmat[6] * inmat[15] + inmat[2] * inmat[7] * inmat[13] + inmat[3] * inmat[5] * inmat[14] - inmat[1] * inmat[7] * inmat[14] - inmat[2] * inmat[5] * inmat[15] - inmat[3] * inmat[6] * inmat[13]) / det;
     59:     outmat[3]  = (inmat[1] * inmat[7] * inmat[10] + inmat[2] * inmat[5] * inmat[11] + inmat[3] * inmat[6] * inmat[9] - inmat[1] * inmat[6] * inmat[11] - inmat[2] * inmat[7] * inmat[9] - inmat[3] * inmat[5] * inmat[10]) / det;
     60:     outmat[4]  = (inmat[4] * inmat[11] * inmat[14] + inmat[6] * inmat[8] * inmat[15] + inmat[7] * inmat[10] * inmat[12] - inmat[4] * inmat[10] * inmat[15] - inmat[6] * inmat[11] * inmat[12] - inmat[7] * inmat[8] * inmat[14]) / det;
     61:     outmat[5]  = (inmat[0] * inmat[10] * inmat[15] + inmat[2] * inmat[11] * inmat[12] + inmat[3] * inmat[8] * inmat[14] - inmat[0] * inmat[11] * inmat[14] - inmat[2] * inmat[8] * inmat[15] - inmat[3] * inmat[10] * inmat[12]) / det;
     62:     outmat[6]  = (inmat[0] * inmat[7] * inmat[14] + inmat[2] * inmat[4] * inmat[15] + inmat[3] * inmat[6] * inmat[12] - inmat[0] * inmat[6] * inmat[15] - inmat[2] * inmat[7] * inmat[12] - inmat[3] * inmat[4] * inmat[14]) / det;
     63:     outmat[7]  = (inmat[0] * inmat[6] * inmat[11] + inmat[2] * inmat[7] * inmat[8] + inmat[3] * inmat[4] * inmat[10] - inmat[0] * inmat[7] * inmat[10] - inmat[2] * inmat[4] * inmat[11] - inmat[3] * inmat[6] * inmat[8]) / det;
     64:     outmat[8]  = (inmat[4] * inmat[9] * inmat[15] + inmat[5] * inmat[11] * inmat[12] + inmat[7] * inmat[8] * inmat[13] - inmat[4] * inmat[11] * inmat[13] - inmat[5] * inmat[8] * inmat[15] - inmat[7] * inmat[9] * inmat[12]) / det;
     65:     outmat[9]  = (inmat[0] * inmat[11] * inmat[13] + inmat[1] * inmat[8] * inmat[15] + inmat[3] * inmat[9] * inmat[12] - inmat[0] * inmat[9] * inmat[15] - inmat[1] * inmat[11] * inmat[12] - inmat[3] * inmat[8] * inmat[13]) / det;
     66:     outmat[10] = (inmat[0] * inmat[5] * inmat[15] + inmat[1] * inmat[7] * inmat[12] + inmat[3] * inmat[4] * inmat[13] - inmat[0] * inmat[7] * inmat[13] - inmat[1] * inmat[4] * inmat[15] - inmat[3] * inmat[5] * inmat[12]) / det;
     67:     outmat[11] = (inmat[0] * inmat[7] * inmat[9] + inmat[1] * inmat[4] * inmat[11] + inmat[3] * inmat[5] * inmat[8] - inmat[0] * inmat[5] * inmat[11] - inmat[1] * inmat[7] * inmat[8] - inmat[3] * inmat[4] * inmat[9]) / det;
     68:     outmat[12] = (inmat[4] * inmat[10] * inmat[13] + inmat[5] * inmat[8] * inmat[14] + inmat[6] * inmat[9] * inmat[12] - inmat[4] * inmat[9] * inmat[14] - inmat[5] * inmat[10] * inmat[12] - inmat[6] * inmat[8] * inmat[13]) / det;
     69:     outmat[13] = (inmat[0] * inmat[9] * inmat[14] + inmat[1] * inmat[10] * inmat[12] + inmat[2] * inmat[8] * inmat[13] - inmat[0] * inmat[10] * inmat[13] - inmat[1] * inmat[8] * inmat[14] - inmat[2] * inmat[9] * inmat[12]) / det;
     70:     outmat[14] = (inmat[0] * inmat[6] * inmat[13] + inmat[1] * inmat[4] * inmat[14] + inmat[2] * inmat[5] * inmat[12] - inmat[0] * inmat[5] * inmat[14] - inmat[1] * inmat[6] * inmat[12] - inmat[2] * inmat[4] * inmat[13]) / det;
     71:     outmat[15] = (inmat[0] * inmat[5] * inmat[10] + inmat[1] * inmat[6] * inmat[8] + inmat[2] * inmat[4] * inmat[9] - inmat[0] * inmat[6] * inmat[9] - inmat[1] * inmat[4] * inmat[10] - inmat[2] * inmat[5] * inmat[8]) / det;
     72:   }
     73:   if (determinant) *determinant = det;
     74:   return PETSC_SUCCESS;
     75: }

     77: /*
     78:   Compute_Lagrange_Basis_1D_Internal - Evaluate bases and derivatives at quadrature points for a EDGE2 or EDGE3 element.

     80:   The routine is given the coordinates of the vertices of a linear or quadratic edge element.

     82:   The routine evaluates the basis functions associated with each quadrature point provided,
     83:   and their derivatives with respect to X.

     85:   Notes:

     87:   Example Physical Element
     88: .vb
     89:     1-------2        1----3----2
     90:       EDGE2             EDGE3
     91: .ve

     93:   Input Parameters:
     94: +  PetscInt  nverts -          the number of element vertices
     95: .  PetscReal coords[3*nverts] - the physical coordinates of the vertices (in canonical numbering)
     96: .  PetscInt  npts -            the number of evaluation points (quadrature points)
     97: -  PetscReal quad[3*npts] -    the evaluation points (quadrature points) in the reference space

     99:   Output Parameters:
    100: +  PetscReal phypts[3*npts] -  the evaluation points (quadrature points) transformed to the physical space
    101: .  PetscReal jxw[npts] -       the jacobian determinant * quadrature weight necessary for assembling discrete contributions
    102: .  PetscReal phi[npts] -       the bases evaluated at the specified quadrature points
    103: .  PetscReal dphidx[npts] -    the derivative of the bases wrt X-direction evaluated at the specified quadrature points
    104: . jacobian  - jacobian
    105: . ijacobian - ijacobian
    106: - volume    - volume

    108:   Level: advanced
    109: */
    110: static PetscErrorCode Compute_Lagrange_Basis_1D_Internal(const PetscInt nverts, const PetscReal *coords, const PetscInt npts, const PetscReal *quad, PetscReal *phypts, PetscReal *jxw, PetscReal *phi, PetscReal *dphidx, PetscReal *jacobian, PetscReal *ijacobian, PetscReal *volume)
    111: {
    112:   int i, j;

    114:   PetscFunctionBegin;
    115:   PetscAssertPointer(jacobian, 9);
    116:   PetscAssertPointer(ijacobian, 10);
    117:   PetscAssertPointer(volume, 11);
    118:   if (phypts) PetscCall(PetscArrayzero(phypts, npts * 3));
    119:   if (dphidx) { /* Reset arrays. */
    120:     PetscCall(PetscArrayzero(dphidx, npts * nverts));
    121:   }
    122:   if (nverts == 2) { /* Linear Edge */

    124:     for (j = 0; j < npts; j++) {
    125:       const PetscInt  offset = j * nverts;
    126:       const PetscReal r      = quad[j];

    128:       phi[0 + offset] = (1.0 - r);
    129:       phi[1 + offset] = (r);

    131:       const PetscReal dNi_dxi[2] = {-1.0, 1.0};

    133:       jacobian[0] = ijacobian[0] = volume[0] = 0.0;
    134:       for (i = 0; i < nverts; ++i) {
    135:         const PetscReal *vertices = coords + i * 3;
    136:         jacobian[0] += dNi_dxi[i] * vertices[0];
    137:         if (phypts) phypts[3 * j + 0] += phi[i + offset] * vertices[0];
    138:       }

    140:       /* invert the jacobian */
    141:       *volume      = jacobian[0];
    142:       ijacobian[0] = 1.0 / jacobian[0];
    143:       jxw[j] *= *volume;

    145:       /*  Divide by element jacobian. */
    146:       for (i = 0; i < nverts; i++) {
    147:         if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[0];
    148:       }
    149:     }
    150:   } else if (nverts == 3) { /* Quadratic Edge */

    152:     for (j = 0; j < npts; j++) {
    153:       const PetscInt  offset = j * nverts;
    154:       const PetscReal r      = quad[j];

    156:       phi[0 + offset] = 1.0 + r * (2.0 * r - 3.0);
    157:       phi[1 + offset] = 4.0 * r * (1.0 - r);
    158:       phi[2 + offset] = r * (2.0 * r - 1.0);

    160:       const PetscReal dNi_dxi[3] = {4 * r - 3.0, 4 * (1.0 - 2.0 * r), 4.0 * r - 1.0};

    162:       jacobian[0] = ijacobian[0] = volume[0] = 0.0;
    163:       for (i = 0; i < nverts; ++i) {
    164:         const PetscReal *vertices = coords + i * 3;
    165:         jacobian[0] += dNi_dxi[i] * vertices[0];
    166:         if (phypts) phypts[3 * j + 0] += phi[i + offset] * vertices[0];
    167:       }

    169:       /* invert the jacobian */
    170:       *volume      = jacobian[0];
    171:       ijacobian[0] = 1.0 / jacobian[0];
    172:       if (jxw) jxw[j] *= *volume;

    174:       /*  Divide by element jacobian. */
    175:       for (i = 0; i < nverts; i++) {
    176:         if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[0];
    177:       }
    178:     }
    179:   } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "The number of entity vertices are invalid. Currently only support EDGE2 and EDGE3 basis evaluations in 1-D : %" PetscInt_FMT, nverts);
    180:   PetscFunctionReturn(PETSC_SUCCESS);
    181: }

    183: /*
    184:   Compute_Lagrange_Basis_2D_Internal - Evaluate bases and derivatives at quadrature points for a QUAD4 or TRI3 element.

    186:   The routine is given the coordinates of the vertices of a quadrangle or triangle.

    188:   The routine evaluates the basis functions associated with each quadrature point provided,
    189:   and their derivatives with respect to X and Y.

    191:   Notes:

    193:   Example Physical Element (QUAD4)
    194: .vb
    195:     4------3        s
    196:     |      |        |
    197:     |      |        |
    198:     |      |        |
    199:     1------2        0-------r
    200: .ve

    202:   Input Parameters:
    203: +  PetscInt  nverts -          the number of element vertices
    204: .  PetscReal coords[3*nverts] - the physical coordinates of the vertices (in canonical numbering)
    205: .  PetscInt  npts -            the number of evaluation points (quadrature points)
    206: -  PetscReal quad[3*npts] -    the evaluation points (quadrature points) in the reference space

    208:   Output Parameters:
    209: +  PetscReal phypts[3*npts] -  the evaluation points (quadrature points) transformed to the physical space
    210: .  PetscReal jxw[npts] -       the jacobian determinant * quadrature weight necessary for assembling discrete contributions
    211: .  PetscReal phi[npts] -       the bases evaluated at the specified quadrature points
    212: .  PetscReal dphidx[npts] -    the derivative of the bases wrt X-direction evaluated at the specified quadrature points
    213: .  PetscReal dphidy[npts] -    the derivative of the bases wrt Y-direction evaluated at the specified quadrature points
    214: . jacobian  - jacobian
    215: . ijacobian - ijacobian
    216: - volume    - volume

    218:   Level: advanced
    219: */
    220: static PetscErrorCode Compute_Lagrange_Basis_2D_Internal(const PetscInt nverts, const PetscReal *coords, const PetscInt npts, const PetscReal *quad, PetscReal *phypts, PetscReal *jxw, PetscReal *phi, PetscReal *dphidx, PetscReal *dphidy, PetscReal *jacobian, PetscReal *ijacobian, PetscReal *volume)
    221: {
    222:   PetscInt i, j, k;

    224:   PetscFunctionBegin;
    225:   PetscAssertPointer(jacobian, 10);
    226:   PetscAssertPointer(ijacobian, 11);
    227:   PetscAssertPointer(volume, 12);
    228:   PetscCall(PetscArrayzero(phi, npts));
    229:   if (phypts) PetscCall(PetscArrayzero(phypts, npts * 3));
    230:   if (dphidx) { /* Reset arrays. */
    231:     PetscCall(PetscArrayzero(dphidx, npts * nverts));
    232:     PetscCall(PetscArrayzero(dphidy, npts * nverts));
    233:   }
    234:   if (nverts == 4) { /* Linear Quadrangle */

    236:     for (j = 0; j < npts; j++) {
    237:       const PetscInt  offset = j * nverts;
    238:       const PetscReal r      = quad[0 + j * 2];
    239:       const PetscReal s      = quad[1 + j * 2];

    241:       phi[0 + offset] = (1.0 - r) * (1.0 - s);
    242:       phi[1 + offset] = r * (1.0 - s);
    243:       phi[2 + offset] = r * s;
    244:       phi[3 + offset] = (1.0 - r) * s;

    246:       const PetscReal dNi_dxi[4]  = {-1.0 + s, 1.0 - s, s, -s};
    247:       const PetscReal dNi_deta[4] = {-1.0 + r, -r, r, 1.0 - r};

    249:       PetscCall(PetscArrayzero(jacobian, 4));
    250:       PetscCall(PetscArrayzero(ijacobian, 4));
    251:       for (i = 0; i < nverts; ++i) {
    252:         const PetscReal *vertices = coords + i * 3;
    253:         jacobian[0] += dNi_dxi[i] * vertices[0];
    254:         jacobian[2] += dNi_dxi[i] * vertices[1];
    255:         jacobian[1] += dNi_deta[i] * vertices[0];
    256:         jacobian[3] += dNi_deta[i] * vertices[1];
    257:         if (phypts) {
    258:           phypts[3 * j + 0] += phi[i + offset] * vertices[0];
    259:           phypts[3 * j + 1] += phi[i + offset] * vertices[1];
    260:           phypts[3 * j + 2] += phi[i + offset] * vertices[2];
    261:         }
    262:       }

    264:       /* invert the jacobian */
    265:       PetscCall(DMatrix_Invert_2x2_Internal(jacobian, ijacobian, volume));
    266:       PetscCheck(*volume >= 1e-12, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Quadrangular element has zero volume: %g. Degenerate element or invalid connectivity", *volume);

    268:       if (jxw) jxw[j] *= *volume;

    270:       /*  Let us compute the bases derivatives by scaling with inverse jacobian. */
    271:       if (dphidx) {
    272:         for (i = 0; i < nverts; i++) {
    273:           for (k = 0; k < 2; ++k) {
    274:             if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[k * 2 + 0];
    275:             if (dphidy) dphidy[i + offset] += dNi_deta[i] * ijacobian[k * 2 + 1];
    276:           } /* for k=[0..2) */
    277:         } /* for i=[0..nverts) */
    278:       } /* if (dphidx) */
    279:     }
    280:   } else if (nverts == 3) { /* Linear triangle */
    281:     const PetscReal x2 = coords[2 * 3 + 0], y2 = coords[2 * 3 + 1];

    283:     PetscCall(PetscArrayzero(jacobian, 4));
    284:     PetscCall(PetscArrayzero(ijacobian, 4));

    286:     /* Jacobian is constant */
    287:     jacobian[0] = (coords[0 * 3 + 0] - x2);
    288:     jacobian[1] = (coords[1 * 3 + 0] - x2);
    289:     jacobian[2] = (coords[0 * 3 + 1] - y2);
    290:     jacobian[3] = (coords[1 * 3 + 1] - y2);

    292:     /* invert the jacobian */
    293:     PetscCall(DMatrix_Invert_2x2_Internal(jacobian, ijacobian, volume));
    294:     PetscCheck(*volume >= 1e-12, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Triangular element has zero volume: %g. Degenerate element or invalid connectivity", (double)*volume);

    296:     const PetscReal Dx[3] = {ijacobian[0], ijacobian[2], -ijacobian[0] - ijacobian[2]};
    297:     const PetscReal Dy[3] = {ijacobian[1], ijacobian[3], -ijacobian[1] - ijacobian[3]};

    299:     for (j = 0; j < npts; j++) {
    300:       const PetscInt  offset = j * nverts;
    301:       const PetscReal r      = quad[0 + j * 2];
    302:       const PetscReal s      = quad[1 + j * 2];

    304:       if (jxw) jxw[j] *= 0.5;

    306:       /* const PetscReal Ni[3]  = { r, s, 1.0 - r - s }; */
    307:       const PetscReal phipts_x = coords[2 * 3 + 0] + jacobian[0] * r + jacobian[1] * s;
    308:       const PetscReal phipts_y = coords[2 * 3 + 1] + jacobian[2] * r + jacobian[3] * s;
    309:       if (phypts) {
    310:         phypts[offset + 0] = phipts_x;
    311:         phypts[offset + 1] = phipts_y;
    312:       }

    314:       /* \phi_0 = (b.y - c.y) x + (b.x - c.x) y + c.x b.y - b.x c.y */
    315:       phi[0 + offset] = (ijacobian[0] * (phipts_x - x2) + ijacobian[1] * (phipts_y - y2));
    316:       /* \phi_1 = (c.y - a.y) x + (a.x - c.x) y + c.x a.y - a.x c.y */
    317:       phi[1 + offset] = (ijacobian[2] * (phipts_x - x2) + ijacobian[3] * (phipts_y - y2));
    318:       phi[2 + offset] = 1.0 - phi[0 + offset] - phi[1 + offset];

    320:       if (dphidx) {
    321:         dphidx[0 + offset] = Dx[0];
    322:         dphidx[1 + offset] = Dx[1];
    323:         dphidx[2 + offset] = Dx[2];
    324:       }

    326:       if (dphidy) {
    327:         dphidy[0 + offset] = Dy[0];
    328:         dphidy[1 + offset] = Dy[1];
    329:         dphidy[2 + offset] = Dy[2];
    330:       }
    331:     }
    332:   } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "The number of entity vertices are invalid. Currently only support QUAD4 and TRI3 basis evaluations in 2-D : %" PetscInt_FMT, nverts);
    333:   PetscFunctionReturn(PETSC_SUCCESS);
    334: }

    336: /*
    337:   Compute_Lagrange_Basis_3D_Internal - Evaluate bases and derivatives at quadrature points for a HEX8 or TET4 element.

    339:   The routine is given the coordinates of the vertices of a hexahedra or tetrahedra.

    341:   The routine evaluates the basis functions associated with each quadrature point provided,
    342:   and their derivatives with respect to X, Y and Z.

    344:   Notes:

    346:   Example Physical Element (HEX8)
    347: .vb
    348:       8------7
    349:      /|     /|        t  s
    350:     5------6 |        | /
    351:     | |    | |        |/
    352:     | 4----|-3        0-------r
    353:     |/     |/
    354:     1------2
    355: .ve

    357:   Input Parameters:
    358: +  PetscInt  nverts -          the number of element vertices
    359: .  PetscReal coords[3*nverts] - the physical coordinates of the vertices (in canonical numbering)
    360: .  PetscInt  npts -            the number of evaluation points (quadrature points)
    361: -  PetscReal quad[3*npts] -    the evaluation points (quadrature points) in the reference space

    363:   Output Parameters:
    364: +  PetscReal phypts[3*npts] -  the evaluation points (quadrature points) transformed to the physical space
    365: .  PetscReal jxw[npts] -       the jacobian determinant * quadrature weight necessary for assembling discrete contributions
    366: .  PetscReal phi[npts] -       the bases evaluated at the specified quadrature points
    367: .  PetscReal dphidx[npts] -    the derivative of the bases wrt X-direction evaluated at the specified quadrature points
    368: .  PetscReal dphidy[npts] -    the derivative of the bases wrt Y-direction evaluated at the specified quadrature points
    369: .  PetscReal dphidz[npts] -    the derivative of the bases wrt Z-direction evaluated at the specified quadrature points
    370: . jacobian  - jacobian
    371: . ijacobian - ijacobian
    372: - volume    - volume

    374:   Level: advanced
    375: */
    376: static PetscErrorCode Compute_Lagrange_Basis_3D_Internal(const PetscInt nverts, const PetscReal *coords, const PetscInt npts, const PetscReal *quad, PetscReal *phypts, PetscReal *jxw, PetscReal *phi, PetscReal *dphidx, PetscReal *dphidy, PetscReal *dphidz, PetscReal *jacobian, PetscReal *ijacobian, PetscReal *volume)
    377: {
    378:   PetscInt i, j, k;

    380:   PetscFunctionBegin;
    381:   PetscAssertPointer(jacobian, 11);
    382:   PetscAssertPointer(ijacobian, 12);
    383:   PetscAssertPointer(volume, 13);

    385:   PetscCall(PetscArrayzero(phi, npts));
    386:   if (phypts) PetscCall(PetscArrayzero(phypts, npts * 3));
    387:   if (dphidx) {
    388:     PetscCall(PetscArrayzero(dphidx, npts * nverts));
    389:     PetscCall(PetscArrayzero(dphidy, npts * nverts));
    390:     PetscCall(PetscArrayzero(dphidz, npts * nverts));
    391:   }

    393:   if (nverts == 8) { /* Linear Hexahedra */

    395:     for (j = 0; j < npts; j++) {
    396:       const PetscInt   offset = j * nverts;
    397:       const PetscReal &r      = quad[j * 3 + 0];
    398:       const PetscReal &s      = quad[j * 3 + 1];
    399:       const PetscReal &t      = quad[j * 3 + 2];

    401:       phi[offset + 0] = (1.0 - r) * (1.0 - s) * (1.0 - t); /* 0,0,0 */
    402:       phi[offset + 1] = (r) * (1.0 - s) * (1.0 - t);       /* 1,0,0 */
    403:       phi[offset + 2] = (r) * (s) * (1.0 - t);             /* 1,1,0 */
    404:       phi[offset + 3] = (1.0 - r) * (s) * (1.0 - t);       /* 0,1,0 */
    405:       phi[offset + 4] = (1.0 - r) * (1.0 - s) * (t);       /* 0,0,1 */
    406:       phi[offset + 5] = (r) * (1.0 - s) * (t);             /* 1,0,1 */
    407:       phi[offset + 6] = (r) * (s) * (t);                   /* 1,1,1 */
    408:       phi[offset + 7] = (1.0 - r) * (s) * (t);             /* 0,1,1 */

    410:       const PetscReal dNi_dxi[8] = {-(1.0 - s) * (1.0 - t), (1.0 - s) * (1.0 - t), (s) * (1.0 - t), -(s) * (1.0 - t), -(1.0 - s) * (t), (1.0 - s) * (t), (s) * (t), -(s) * (t)};

    412:       const PetscReal dNi_deta[8] = {-(1.0 - r) * (1.0 - t), -(r) * (1.0 - t), (r) * (1.0 - t), (1.0 - r) * (1.0 - t), -(1.0 - r) * (t), -(r) * (t), (r) * (t), (1.0 - r) * (t)};

    414:       const PetscReal dNi_dzeta[8] = {-(1.0 - r) * (1.0 - s), -(r) * (1.0 - s), -(r) * (s), -(1.0 - r) * (s), (1.0 - r) * (1.0 - s), (r) * (1.0 - s), (r) * (s), (1.0 - r) * (s)};

    416:       PetscCall(PetscArrayzero(jacobian, 9));
    417:       PetscCall(PetscArrayzero(ijacobian, 9));
    418:       for (i = 0; i < nverts; ++i) {
    419:         const PetscReal *vertex = coords + i * 3;
    420:         jacobian[0] += dNi_dxi[i] * vertex[0];
    421:         jacobian[3] += dNi_dxi[i] * vertex[1];
    422:         jacobian[6] += dNi_dxi[i] * vertex[2];
    423:         jacobian[1] += dNi_deta[i] * vertex[0];
    424:         jacobian[4] += dNi_deta[i] * vertex[1];
    425:         jacobian[7] += dNi_deta[i] * vertex[2];
    426:         jacobian[2] += dNi_dzeta[i] * vertex[0];
    427:         jacobian[5] += dNi_dzeta[i] * vertex[1];
    428:         jacobian[8] += dNi_dzeta[i] * vertex[2];
    429:         if (phypts) {
    430:           phypts[3 * j + 0] += phi[i + offset] * vertex[0];
    431:           phypts[3 * j + 1] += phi[i + offset] * vertex[1];
    432:           phypts[3 * j + 2] += phi[i + offset] * vertex[2];
    433:         }
    434:       }

    436:       /* invert the jacobian */
    437:       PetscCall(DMatrix_Invert_3x3_Internal(jacobian, ijacobian, volume));
    438:       PetscCheck(*volume >= 1e-8, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Hexahedral element has zero volume: %g. Degenerate element or invalid connectivity", *volume);

    440:       if (jxw) jxw[j] *= (*volume);

    442:       /*  Divide by element jacobian. */
    443:       for (i = 0; i < nverts; ++i) {
    444:         for (k = 0; k < 3; ++k) {
    445:           if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[0 * 3 + k];
    446:           if (dphidy) dphidy[i + offset] += dNi_deta[i] * ijacobian[1 * 3 + k];
    447:           if (dphidz) dphidz[i + offset] += dNi_dzeta[i] * ijacobian[2 * 3 + k];
    448:         }
    449:       }
    450:     }
    451:   } else if (nverts == 4) { /* Linear Tetrahedra */
    452:     PetscReal       Dx[4] = {0, 0, 0, 0}, Dy[4] = {0, 0, 0, 0}, Dz[4] = {0, 0, 0, 0};
    453:     const PetscReal x0 = coords[/*0 * 3 +*/ 0], y0 = coords[/*0 * 3 +*/ 1], z0 = coords[/*0 * 3 +*/ 2];

    455:     PetscCall(PetscArrayzero(jacobian, 9));
    456:     PetscCall(PetscArrayzero(ijacobian, 9));

    458:     /* compute the jacobian */
    459:     jacobian[0] = coords[1 * 3 + 0] - x0;
    460:     jacobian[1] = coords[2 * 3 + 0] - x0;
    461:     jacobian[2] = coords[3 * 3 + 0] - x0;
    462:     jacobian[3] = coords[1 * 3 + 1] - y0;
    463:     jacobian[4] = coords[2 * 3 + 1] - y0;
    464:     jacobian[5] = coords[3 * 3 + 1] - y0;
    465:     jacobian[6] = coords[1 * 3 + 2] - z0;
    466:     jacobian[7] = coords[2 * 3 + 2] - z0;
    467:     jacobian[8] = coords[3 * 3 + 2] - z0;

    469:     /* invert the jacobian */
    470:     PetscCall(DMatrix_Invert_3x3_Internal(jacobian, ijacobian, volume));
    471:     PetscCheck(*volume >= 1e-8, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Tetrahedral element has zero volume: %g. Degenerate element or invalid connectivity", (double)*volume);

    473:     if (dphidx) {
    474:       Dx[0] = (coords[1 + 2 * 3] * (coords[2 + 1 * 3] - coords[2 + 3 * 3]) - coords[1 + 1 * 3] * (coords[2 + 2 * 3] - coords[2 + 3 * 3]) - coords[1 + 3 * 3] * (coords[2 + 1 * 3] - coords[2 + 2 * 3])) / *volume;
    475:       Dx[1] = -(coords[1 + 2 * 3] * (coords[2 + 0 * 3] - coords[2 + 3 * 3]) - coords[1 + 0 * 3] * (coords[2 + 2 * 3] - coords[2 + 3 * 3]) - coords[1 + 3 * 3] * (coords[2 + 0 * 3] - coords[2 + 2 * 3])) / *volume;
    476:       Dx[2] = (coords[1 + 1 * 3] * (coords[2 + 0 * 3] - coords[2 + 3 * 3]) - coords[1 + 0 * 3] * (coords[2 + 1 * 3] - coords[2 + 3 * 3]) - coords[1 + 3 * 3] * (coords[2 + 0 * 3] - coords[2 + 1 * 3])) / *volume;
    477:       Dx[3] = -(Dx[0] + Dx[1] + Dx[2]);
    478:       Dy[0] = (coords[0 + 1 * 3] * (coords[2 + 2 * 3] - coords[2 + 3 * 3]) - coords[0 + 2 * 3] * (coords[2 + 1 * 3] - coords[2 + 3 * 3]) + coords[0 + 3 * 3] * (coords[2 + 1 * 3] - coords[2 + 2 * 3])) / *volume;
    479:       Dy[1] = -(coords[0 + 0 * 3] * (coords[2 + 2 * 3] - coords[2 + 3 * 3]) - coords[0 + 2 * 3] * (coords[2 + 0 * 3] - coords[2 + 3 * 3]) + coords[0 + 3 * 3] * (coords[2 + 0 * 3] - coords[2 + 2 * 3])) / *volume;
    480:       Dy[2] = (coords[0 + 0 * 3] * (coords[2 + 1 * 3] - coords[2 + 3 * 3]) - coords[0 + 1 * 3] * (coords[2 + 0 * 3] - coords[2 + 3 * 3]) + coords[0 + 3 * 3] * (coords[2 + 0 * 3] - coords[2 + 1 * 3])) / *volume;
    481:       Dy[3] = -(Dy[0] + Dy[1] + Dy[2]);
    482:       Dz[0] = (coords[0 + 1 * 3] * (coords[1 + 3 * 3] - coords[1 + 2 * 3]) - coords[0 + 2 * 3] * (coords[1 + 3 * 3] - coords[1 + 1 * 3]) + coords[0 + 3 * 3] * (coords[1 + 2 * 3] - coords[1 + 1 * 3])) / *volume;
    483:       Dz[1] = -(coords[0 + 0 * 3] * (coords[1 + 3 * 3] - coords[1 + 2 * 3]) + coords[0 + 2 * 3] * (coords[1 + 0 * 3] - coords[1 + 3 * 3]) - coords[0 + 3 * 3] * (coords[1 + 0 * 3] - coords[1 + 2 * 3])) / *volume;
    484:       Dz[2] = (coords[0 + 0 * 3] * (coords[1 + 3 * 3] - coords[1 + 1 * 3]) + coords[0 + 1 * 3] * (coords[1 + 0 * 3] - coords[1 + 3 * 3]) - coords[0 + 3 * 3] * (coords[1 + 0 * 3] - coords[1 + 1 * 3])) / *volume;
    485:       Dz[3] = -(Dz[0] + Dz[1] + Dz[2]);
    486:     }

    488:     for (j = 0; j < npts; j++) {
    489:       const PetscInt   offset = j * nverts;
    490:       const PetscReal &r      = quad[j * 3 + 0];
    491:       const PetscReal &s      = quad[j * 3 + 1];
    492:       const PetscReal &t      = quad[j * 3 + 2];

    494:       if (jxw) jxw[j] *= *volume;

    496:       phi[offset + 0] = 1.0 - r - s - t;
    497:       phi[offset + 1] = r;
    498:       phi[offset + 2] = s;
    499:       phi[offset + 3] = t;

    501:       if (phypts) {
    502:         for (i = 0; i < nverts; ++i) {
    503:           const PetscReal *vertices = coords + i * 3;
    504:           phypts[3 * j + 0] += phi[i + offset] * vertices[0];
    505:           phypts[3 * j + 1] += phi[i + offset] * vertices[1];
    506:           phypts[3 * j + 2] += phi[i + offset] * vertices[2];
    507:         }
    508:       }

    510:       /* Now set the derivatives */
    511:       if (dphidx) {
    512:         dphidx[0 + offset] = Dx[0];
    513:         dphidx[1 + offset] = Dx[1];
    514:         dphidx[2 + offset] = Dx[2];
    515:         dphidx[3 + offset] = Dx[3];

    517:         dphidy[0 + offset] = Dy[0];
    518:         dphidy[1 + offset] = Dy[1];
    519:         dphidy[2 + offset] = Dy[2];
    520:         dphidy[3 + offset] = Dy[3];

    522:         dphidz[0 + offset] = Dz[0];
    523:         dphidz[1 + offset] = Dz[1];
    524:         dphidz[2 + offset] = Dz[2];
    525:         dphidz[3 + offset] = Dz[3];
    526:       }

    528:     } /* Tetrahedra -- ends */
    529:   } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "The number of entity vertices are invalid. Currently only support HEX8 and TET4 basis evaluations in 3-D : %" PetscInt_FMT, nverts);
    530:   PetscFunctionReturn(PETSC_SUCCESS);
    531: }

    533: /*@C
    534:   DMMoabFEMComputeBasis - Evaluate bases and derivatives at quadrature points for a linear EDGE/QUAD/TRI/HEX/TET element.

    536:   The routine takes the coordinates of the vertices of an element and computes basis functions associated with
    537:   each quadrature point provided, and their derivatives with respect to X, Y and Z as appropriate.

    539:   Input Parameters:
    540: + dim         - the dimension
    541: . nverts      - the number of element vertices
    542: . coordinates - the physical coordinates of the vertices (in canonical numbering)
    543: - quadrature  - the evaluation points (quadrature points) in the reference space

    545:   Output Parameters:
    546: + phypts                             - the evaluation points (quadrature points) transformed to the physical space
    547: . jacobian_quadrature_weight_product - the jacobian determinant * quadrature weight necessary for assembling discrete contributions
    548: . fe_basis                           - the bases values evaluated at the specified quadrature points
    549: - fe_basis_derivatives               - the derivative of the bases wrt (X,Y,Z)-directions (depending on the dimension) evaluated at the specified quadrature points

    551:   Level: advanced

    553: .seealso: `DMMoabCreate()`
    554: @*/
    555: PetscErrorCode DMMoabFEMComputeBasis(const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscQuadrature quadrature, PetscReal *phypts, PetscReal *jacobian_quadrature_weight_product, PetscReal *fe_basis, PetscReal **fe_basis_derivatives)
    556: {
    557:   PetscInt         npoints, idim;
    558:   bool             compute_der;
    559:   const PetscReal *quadpts, *quadwts;
    560:   PetscReal        jacobian[9], ijacobian[9], volume;

    562:   PetscFunctionBegin;
    563:   PetscAssertPointer(coordinates, 3);
    565:   PetscAssertPointer(fe_basis, 7);
    566:   compute_der = (fe_basis_derivatives != NULL);

    568:   /* Get the quadrature points and weights for the given quadrature rule */
    569:   PetscCall(PetscQuadratureGetData(quadrature, &idim, NULL, &npoints, &quadpts, &quadwts));
    570:   PetscCheck(idim == dim, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Dimension mismatch: provided (%" PetscInt_FMT ") vs quadrature (%" PetscInt_FMT ")", idim, dim);
    571:   if (jacobian_quadrature_weight_product) PetscCall(PetscArraycpy(jacobian_quadrature_weight_product, quadwts, npoints));

    573:   switch (dim) {
    574:   case 1:
    575:     PetscCall(Compute_Lagrange_Basis_1D_Internal(nverts, coordinates, npoints, quadpts, phypts, jacobian_quadrature_weight_product, fe_basis, compute_der ? fe_basis_derivatives[0] : NULL, jacobian, ijacobian, &volume));
    576:     break;
    577:   case 2:
    578:     PetscCall(Compute_Lagrange_Basis_2D_Internal(nverts, coordinates, npoints, quadpts, phypts, jacobian_quadrature_weight_product, fe_basis, compute_der ? fe_basis_derivatives[0] : NULL, compute_der ? fe_basis_derivatives[1] : NULL, jacobian, ijacobian, &volume));
    579:     break;
    580:   case 3:
    581:     PetscCall(Compute_Lagrange_Basis_3D_Internal(nverts, coordinates, npoints, quadpts, phypts, jacobian_quadrature_weight_product, fe_basis, compute_der ? fe_basis_derivatives[0] : NULL, compute_der ? fe_basis_derivatives[1] : NULL, compute_der ? fe_basis_derivatives[2] : NULL, jacobian, ijacobian, &volume));
    582:     break;
    583:   default:
    584:     SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %" PetscInt_FMT, dim);
    585:   }
    586:   PetscFunctionReturn(PETSC_SUCCESS);
    587: }

    589: /*@C
    590:   DMMoabFEMCreateQuadratureDefault - Create default quadrature rules for integration over an element with a given
    591:   dimension and polynomial order (deciphered from number of element vertices).

    593:   Input Parameters:
    594: + dim    - the element dimension (1=EDGE, 2=QUAD/TRI, 3=HEX/TET)
    595: - nverts - the number of vertices in the physical element

    597:   Output Parameter:
    598: . quadrature - the quadrature object with default settings to integrate polynomials defined over the element

    600:   Level: advanced

    602: .seealso: `DMMoabCreate()`
    603: @*/
    604: PetscErrorCode DMMoabFEMCreateQuadratureDefault(const PetscInt dim, const PetscInt nverts, PetscQuadrature *quadrature)
    605: {
    606:   PetscReal *w, *x;
    607:   PetscInt   nc = 1;

    609:   PetscFunctionBegin;
    610:   /* Create an appropriate quadrature rule to sample basis */
    611:   switch (dim) {
    612:   case 1:
    613:     /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */
    614:     PetscCall(PetscDTStroudConicalQuadrature(1, nc, nverts, 0, 1.0, quadrature));
    615:     break;
    616:   case 2:
    617:     /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */
    618:     if (nverts == 3) {
    619:       const PetscInt order   = 2;
    620:       const PetscInt npoints = (order == 2 ? 3 : 6);
    621:       PetscCall(PetscMalloc2(npoints * 2, &x, npoints, &w));
    622:       if (npoints == 3) {
    623:         x[0] = x[1] = x[2] = x[5] = 1.0 / 6.0;
    624:         x[3] = x[4] = 2.0 / 3.0;
    625:         w[0] = w[1] = w[2] = 1.0 / 3.0;
    626:       } else if (npoints == 6) {
    627:         x[0] = x[1] = x[2] = 0.44594849091597;
    628:         x[3] = x[4] = 0.10810301816807;
    629:         x[5]        = 0.44594849091597;
    630:         x[6] = x[7] = x[8] = 0.09157621350977;
    631:         x[9] = x[10] = 0.81684757298046;
    632:         x[11]        = 0.09157621350977;
    633:         w[0] = w[1] = w[2] = 0.22338158967801;
    634:         w[3] = w[4] = w[5] = 0.10995174365532;
    635:       } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Triangle quadrature rules for points 3 and 6 supported; npoints : %" PetscInt_FMT, npoints);
    636:       PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, quadrature));
    637:       PetscCall(PetscQuadratureSetOrder(*quadrature, order));
    638:       PetscCall(PetscQuadratureSetData(*quadrature, dim, nc, npoints, x, w));
    639:       /* PetscCall(PetscDTStroudConicalQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature)); */
    640:     } else PetscCall(PetscDTGaussTensorQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature));
    641:     break;
    642:   case 3:
    643:     /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */
    644:     if (nverts == 4) {
    645:       PetscInt       order;
    646:       const PetscInt npoints = 4; // Choose between 4 and 10
    647:       PetscCall(PetscMalloc2(npoints * 3, &x, npoints, &w));
    648:       if (npoints == 4) { /*  KEAST1, K1,  N=4, O=4 */
    649:         const PetscReal x_4[12] = {0.5854101966249685, 0.1381966011250105, 0.1381966011250105, 0.1381966011250105, 0.1381966011250105, 0.1381966011250105,
    650:                                    0.1381966011250105, 0.1381966011250105, 0.5854101966249685, 0.1381966011250105, 0.5854101966249685, 0.1381966011250105};
    651:         PetscCall(PetscArraycpy(x, x_4, 12));

    653:         w[0] = w[1] = w[2] = w[3] = 1.0 / 24.0;
    654:         order                     = 4;
    655:       } else if (npoints == 10) { /*  KEAST3, K3  N=10, O=10 */
    656:         const PetscReal x_10[30] = {0.5684305841968444, 0.1438564719343852, 0.1438564719343852, 0.1438564719343852, 0.1438564719343852, 0.1438564719343852, 0.1438564719343852, 0.1438564719343852, 0.5684305841968444, 0.1438564719343852,
    657:                                     0.5684305841968444, 0.1438564719343852, 0.0000000000000000, 0.5000000000000000, 0.5000000000000000, 0.5000000000000000, 0.0000000000000000, 0.5000000000000000, 0.5000000000000000, 0.5000000000000000,
    658:                                     0.0000000000000000, 0.5000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.5000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.5000000000000000};
    659:         PetscCall(PetscArraycpy(x, x_10, 30));

    661:         w[0] = w[1] = w[2] = w[3] = 0.2177650698804054;
    662:         w[4] = w[5] = w[6] = w[7] = w[8] = w[9] = 0.0214899534130631;
    663:         order                                   = 10;
    664:       } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Tetrahedral quadrature rules for points 4 and 10 supported; npoints : %" PetscInt_FMT, npoints);
    665:       PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, quadrature));
    666:       PetscCall(PetscQuadratureSetOrder(*quadrature, order));
    667:       PetscCall(PetscQuadratureSetData(*quadrature, dim, nc, npoints, x, w));
    668:       /* PetscCall(PetscDTStroudConicalQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature)); */
    669:     } else PetscCall(PetscDTGaussTensorQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature));
    670:     break;
    671:   default:
    672:     SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %" PetscInt_FMT, dim);
    673:   }
    674:   PetscFunctionReturn(PETSC_SUCCESS);
    675: }

    677: /* Compute Jacobians */
    678: static PetscErrorCode FEMComputeBasis_JandF(const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal *quadrature, PetscReal *phypts, PetscReal *phibasis, PetscReal *jacobian, PetscReal *ijacobian, PetscReal *volume)
    679: {
    680:   PetscFunctionBegin;
    681:   switch (dim) {
    682:   case 1:
    683:     PetscCall(Compute_Lagrange_Basis_1D_Internal(nverts, coordinates, 1, quadrature, phypts, NULL, phibasis, NULL, jacobian, ijacobian, volume));
    684:     break;
    685:   case 2:
    686:     PetscCall(Compute_Lagrange_Basis_2D_Internal(nverts, coordinates, 1, quadrature, phypts, NULL, phibasis, NULL, NULL, jacobian, ijacobian, volume));
    687:     break;
    688:   case 3:
    689:     PetscCall(Compute_Lagrange_Basis_3D_Internal(nverts, coordinates, 1, quadrature, phypts, NULL, phibasis, NULL, NULL, NULL, jacobian, ijacobian, volume));
    690:     break;
    691:   default:
    692:     SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %" PetscInt_FMT, dim);
    693:   }
    694:   PetscFunctionReturn(PETSC_SUCCESS);
    695: }

    697: /*@C
    698:   DMMoabPToRMapping - Compute the mapping from the physical coordinate system for a given element to the
    699:   canonical reference element.

    701:   Input Parameters:
    702: + dim         - the element dimension (1=EDGE, 2=QUAD/TRI, 3=HEX/TET)
    703: . nverts      - the number of vertices in the physical element
    704: . coordinates - the coordinates of vertices in the physical element
    705: - xphy        - the coordinates of physical point for which natural coordinates (in reference frame) are sought

    707:   Output Parameters:
    708: + natparam - the natural coordinates (in reference frame) corresponding to xphy
    709: - phi      - the basis functions evaluated at the natural coordinates (natparam)

    711:   Level: advanced

    713:   Notes:
    714:   In addition to finding the inverse mapping evaluation through Newton iteration, the basis
    715:   function at the parametric point is also evaluated optionally.

    717: .seealso: `DMMoabCreate()`
    718: @*/
    719: PetscErrorCode DMMoabPToRMapping(const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal *xphy, PetscReal *natparam, PetscReal *phi)
    720: {
    721:   /* Perform inverse evaluation for the mapping with use of Newton Raphson iteration */
    722:   const PetscReal tol            = 1.0e-10;
    723:   const PetscInt  max_iterations = 10;
    724:   const PetscReal error_tol_sqr  = tol * tol;
    725:   PetscReal       phibasis[8], jacobian[9], ijacobian[9], volume;
    726:   PetscReal       phypts[3] = {0.0, 0.0, 0.0};
    727:   PetscReal       delta[3]  = {0.0, 0.0, 0.0};
    728:   PetscInt        iters     = 0;
    729:   PetscReal       error     = 1.0;

    731:   PetscFunctionBegin;
    732:   PetscAssertPointer(coordinates, 3);
    733:   PetscAssertPointer(xphy, 4);
    734:   PetscAssertPointer(natparam, 5);

    736:   PetscCall(PetscArrayzero(jacobian, dim * dim));
    737:   PetscCall(PetscArrayzero(ijacobian, dim * dim));
    738:   PetscCall(PetscArrayzero(phibasis, nverts));

    740:   /* zero initial guess */
    741:   natparam[0] = natparam[1] = natparam[2] = 0.0;

    743:   /* Compute delta = evaluate( xi) - x */
    744:   PetscCall(FEMComputeBasis_JandF(dim, nverts, coordinates, natparam, phypts, phibasis, jacobian, ijacobian, &volume));

    746:   error = 0.0;
    747:   switch (dim) {
    748:   case 3:
    749:     delta[2] = phypts[2] - xphy[2];
    750:     error += (delta[2] * delta[2]);
    751:   case 2:
    752:     delta[1] = phypts[1] - xphy[1];
    753:     error += (delta[1] * delta[1]);
    754:   case 1:
    755:     delta[0] = phypts[0] - xphy[0];
    756:     error += (delta[0] * delta[0]);
    757:     break;
    758:   }

    760:   while (error > error_tol_sqr) {
    761:     PetscCheck(++iters <= max_iterations, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Maximum Newton iterations (10) reached. Current point in reference space : (%g, %g, %g)", natparam[0], natparam[1], natparam[2]);

    763:     /* Compute natparam -= J.inverse() * delta */
    764:     switch (dim) {
    765:     case 1:
    766:       natparam[0] -= ijacobian[0] * delta[0];
    767:       break;
    768:     case 2:
    769:       natparam[0] -= ijacobian[0] * delta[0] + ijacobian[1] * delta[1];
    770:       natparam[1] -= ijacobian[2] * delta[0] + ijacobian[3] * delta[1];
    771:       break;
    772:     case 3:
    773:       natparam[0] -= ijacobian[0] * delta[0] + ijacobian[1] * delta[1] + ijacobian[2] * delta[2];
    774:       natparam[1] -= ijacobian[3] * delta[0] + ijacobian[4] * delta[1] + ijacobian[5] * delta[2];
    775:       natparam[2] -= ijacobian[6] * delta[0] + ijacobian[7] * delta[1] + ijacobian[8] * delta[2];
    776:       break;
    777:     }

    779:     /* Compute delta = evaluate( xi) - x */
    780:     PetscCall(FEMComputeBasis_JandF(dim, nverts, coordinates, natparam, phypts, phibasis, jacobian, ijacobian, &volume));

    782:     error = 0.0;
    783:     switch (dim) {
    784:     case 3:
    785:       delta[2] = phypts[2] - xphy[2];
    786:       error += (delta[2] * delta[2]);
    787:     case 2:
    788:       delta[1] = phypts[1] - xphy[1];
    789:       error += (delta[1] * delta[1]);
    790:     case 1:
    791:       delta[0] = phypts[0] - xphy[0];
    792:       error += (delta[0] * delta[0]);
    793:       break;
    794:     }
    795:   }
    796:   if (phi) PetscCall(PetscArraycpy(phi, phibasis, nverts));
    797:   PetscFunctionReturn(PETSC_SUCCESS);
    798: }