Actual source code: dmmbfem.cxx


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

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

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

 25: static inline PetscReal DMatrix_Determinant_3x3_Internal(const PetscReal inmat[3 * 3])
 26: {
 27:   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]);
 28: }

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

 48: inline PetscReal DMatrix_Determinant_4x4_Internal(PetscReal inmat[4 * 4])
 49: {
 50:   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]));
 51: }

 53: inline PetscErrorCode DMatrix_Invert_4x4_Internal(PetscReal *inmat, PetscReal *outmat, PetscScalar *determinant)
 54: {
 55:   PetscReal det = DMatrix_Determinant_4x4_Internal(inmat);
 56:   if (outmat) {
 57:     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;
 58:     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;
 59:     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;
 60:     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;
 61:     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;
 62:     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;
 63:     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;
 64:     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;
 65:     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;
 66:     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;
 67:     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;
 68:     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;
 69:     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;
 70:     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;
 71:     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;
 72:     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;
 73:   }
 74:   if (determinant) *determinant = det;
 75:   return PETSC_SUCCESS;
 76: }

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

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

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

 86:   Notes:

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

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

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

109:   Level: advanced

111: @*/
112: 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)
113: {
114:   int i, j;

116:   PetscFunctionBegin;
120:   if (phypts) PetscCall(PetscArrayzero(phypts, npts * 3));
121:   if (dphidx) { /* Reset arrays. */
122:     PetscCall(PetscArrayzero(dphidx, npts * nverts));
123:   }
124:   if (nverts == 2) { /* Linear Edge */

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

130:       phi[0 + offset] = (1.0 - r);
131:       phi[1 + offset] = (r);

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

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

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

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

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

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

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

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

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

176:       /*  Divide by element jacobian. */
177:       for (i = 0; i < nverts; i++) {
178:         if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[0];
179:       }
180:     }
181:   } 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);
182:   PetscFunctionReturn(PETSC_SUCCESS);
183: }

185: /*@C
186:   Compute_Lagrange_Basis_2D_Internal - Evaluate bases and derivatives at quadrature points for a QUAD4 or TRI3 element.

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

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

193:   Notes:

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

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

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

220:   Level: advanced

222: @*/
223: 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)
224: {
225:   PetscInt i, j, k;

227:   PetscFunctionBegin;
231:   PetscCall(PetscArrayzero(phi, npts));
232:   if (phypts) PetscCall(PetscArrayzero(phypts, npts * 3));
233:   if (dphidx) { /* Reset arrays. */
234:     PetscCall(PetscArrayzero(dphidx, npts * nverts));
235:     PetscCall(PetscArrayzero(dphidy, npts * nverts));
236:   }
237:   if (nverts == 4) { /* Linear Quadrangle */

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

244:       phi[0 + offset] = (1.0 - r) * (1.0 - s);
245:       phi[1 + offset] = r * (1.0 - s);
246:       phi[2 + offset] = r * s;
247:       phi[3 + offset] = (1.0 - r) * s;

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

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

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

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

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

286:     PetscCall(PetscArrayzero(jacobian, 4));
287:     PetscCall(PetscArrayzero(ijacobian, 4));

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

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

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

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

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

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

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

323:       if (dphidx) {
324:         dphidx[0 + offset] = Dx[0];
325:         dphidx[1 + offset] = Dx[1];
326:         dphidx[2 + offset] = Dx[2];
327:       }

329:       if (dphidy) {
330:         dphidy[0 + offset] = Dy[0];
331:         dphidy[1 + offset] = Dy[1];
332:         dphidy[2 + offset] = Dy[2];
333:       }
334:     }
335:   } 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);
336:   PetscFunctionReturn(PETSC_SUCCESS);
337: }

339: /*@C
340:   Compute_Lagrange_Basis_3D_Internal - Evaluate bases and derivatives at quadrature points for a HEX8 or TET4 element.

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

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

347:   Notes:

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

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

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

377:   Level: advanced

379: @*/
380: 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)
381: {
382:   PetscInt i, j, k;

384:   PetscFunctionBegin;

389:   PetscCall(PetscArrayzero(phi, npts));
390:   if (phypts) PetscCall(PetscArrayzero(phypts, npts * 3));
391:   if (dphidx) {
392:     PetscCall(PetscArrayzero(dphidx, npts * nverts));
393:     PetscCall(PetscArrayzero(dphidy, npts * nverts));
394:     PetscCall(PetscArrayzero(dphidz, npts * nverts));
395:   }

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

399:     for (j = 0; j < npts; j++) {
400:       const PetscInt   offset = j * nverts;
401:       const PetscReal &r      = quad[j * 3 + 0];
402:       const PetscReal &s      = quad[j * 3 + 1];
403:       const PetscReal &t      = quad[j * 3 + 2];

405:       phi[offset + 0] = (1.0 - r) * (1.0 - s) * (1.0 - t); /* 0,0,0 */
406:       phi[offset + 1] = (r) * (1.0 - s) * (1.0 - t);       /* 1,0,0 */
407:       phi[offset + 2] = (r) * (s) * (1.0 - t);             /* 1,1,0 */
408:       phi[offset + 3] = (1.0 - r) * (s) * (1.0 - t);       /* 0,1,0 */
409:       phi[offset + 4] = (1.0 - r) * (1.0 - s) * (t);       /* 0,0,1 */
410:       phi[offset + 5] = (r) * (1.0 - s) * (t);             /* 1,0,1 */
411:       phi[offset + 6] = (r) * (s) * (t);                   /* 1,1,1 */
412:       phi[offset + 7] = (1.0 - r) * (s) * (t);             /* 0,1,1 */

414:       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)};

416:       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)};

418:       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)};

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

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

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

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

459:     PetscCall(PetscArrayzero(jacobian, 9));
460:     PetscCall(PetscArrayzero(ijacobian, 9));

462:     /* compute the jacobian */
463:     jacobian[0] = coords[1 * 3 + 0] - x0;
464:     jacobian[1] = coords[2 * 3 + 0] - x0;
465:     jacobian[2] = coords[3 * 3 + 0] - x0;
466:     jacobian[3] = coords[1 * 3 + 1] - y0;
467:     jacobian[4] = coords[2 * 3 + 1] - y0;
468:     jacobian[5] = coords[3 * 3 + 1] - y0;
469:     jacobian[6] = coords[1 * 3 + 2] - z0;
470:     jacobian[7] = coords[2 * 3 + 2] - z0;
471:     jacobian[8] = coords[3 * 3 + 2] - z0;

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

477:     if (dphidx) {
478:       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;
479:       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;
480:       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;
481:       Dx[3] = -(Dx[0] + Dx[1] + Dx[2]);
482:       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;
483:       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;
484:       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;
485:       Dy[3] = -(Dy[0] + Dy[1] + Dy[2]);
486:       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;
487:       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;
488:       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;
489:       Dz[3] = -(Dz[0] + Dz[1] + Dz[2]);
490:     }

492:     for (j = 0; j < npts; j++) {
493:       const PetscInt   offset = j * nverts;
494:       const PetscReal &r      = quad[j * 3 + 0];
495:       const PetscReal &s      = quad[j * 3 + 1];
496:       const PetscReal &t      = quad[j * 3 + 2];

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

500:       phi[offset + 0] = 1.0 - r - s - t;
501:       phi[offset + 1] = r;
502:       phi[offset + 2] = s;
503:       phi[offset + 3] = t;

505:       if (phypts) {
506:         for (i = 0; i < nverts; ++i) {
507:           const PetscScalar *vertices = coords + i * 3;
508:           phypts[3 * j + 0] += phi[i + offset] * vertices[0];
509:           phypts[3 * j + 1] += phi[i + offset] * vertices[1];
510:           phypts[3 * j + 2] += phi[i + offset] * vertices[2];
511:         }
512:       }

514:       /* Now set the derivatives */
515:       if (dphidx) {
516:         dphidx[0 + offset] = Dx[0];
517:         dphidx[1 + offset] = Dx[1];
518:         dphidx[2 + offset] = Dx[2];
519:         dphidx[3 + offset] = Dx[3];

521:         dphidy[0 + offset] = Dy[0];
522:         dphidy[1 + offset] = Dy[1];
523:         dphidy[2 + offset] = Dy[2];
524:         dphidy[3 + offset] = Dy[3];

526:         dphidz[0 + offset] = Dz[0];
527:         dphidz[1 + offset] = Dz[1];
528:         dphidz[2 + offset] = Dz[2];
529:         dphidz[3 + offset] = Dz[3];
530:       }

532:     } /* Tetrahedra -- ends */
533:   } 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);
534:   PetscFunctionReturn(PETSC_SUCCESS);
535: }

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

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

543:   Input Parameters:
544: +  PetscInt  nverts -           the number of element vertices
545: .  PetscReal coords[3*nverts] - the physical coordinates of the vertices (in canonical numbering)
546: .  PetscInt  npts -             the number of evaluation points (quadrature points)
547: -  PetscReal quad[3*npts] -     the evaluation points (quadrature points) in the reference space

549:   Output Parameters:
550: +  PetscReal phypts[3*npts] -   the evaluation points (quadrature points) transformed to the physical space
551: .  PetscReal jxw[npts] -        the jacobian determinant * quadrature weight necessary for assembling discrete contributions
552: .  PetscReal fe_basis[npts] -   the bases values evaluated at the specified quadrature points
553: -  PetscReal fe_basis_derivatives[dim][npts] - the derivative of the bases wrt (X,Y,Z)-directions (depending on the dimension) evaluated at the specified quadrature points

555:   Level: advanced

557: @*/
558: 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)
559: {
560:   PetscInt         npoints, idim;
561:   bool             compute_der;
562:   const PetscReal *quadpts, *quadwts;
563:   PetscReal        jacobian[9], ijacobian[9], volume;

565:   PetscFunctionBegin;
569:   compute_der = (fe_basis_derivatives != NULL);

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

576:   switch (dim) {
577:   case 1:
578:     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));
579:     break;
580:   case 2:
581:     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));
582:     break;
583:   case 3:
584:     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));
585:     break;
586:   default:
587:     SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %" PetscInt_FMT, dim);
588:   }
589:   PetscFunctionReturn(PETSC_SUCCESS);
590: }

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

596:   Input Parameters:

598: +  PetscInt  dim   -   the element dimension (1=EDGE, 2=QUAD/TRI, 3=HEX/TET)
599: -  PetscInt nverts -   the number of vertices in the physical element

601:   Output Parameter:
602: .  PetscQuadrature quadrature -  the quadrature object with default settings to integrate polynomials defined over the element

604:   Level: advanced

606: @*/
607: PetscErrorCode DMMoabFEMCreateQuadratureDefault(const PetscInt dim, const PetscInt nverts, PetscQuadrature *quadrature)
608: {
609:   PetscReal *w, *x;
610:   PetscInt   nc = 1;

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

656:         w[0] = w[1] = w[2] = w[3] = 1.0 / 24.0;
657:         order                     = 4;
658:       } else if (npoints == 10) { /*  KEAST3, K3  N=10, O=10 */
659:         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,
660:                                     0.5684305841968444, 0.1438564719343852, 0.0000000000000000, 0.5000000000000000, 0.5000000000000000, 0.5000000000000000, 0.0000000000000000, 0.5000000000000000, 0.5000000000000000, 0.5000000000000000,
661:                                     0.0000000000000000, 0.5000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.5000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.5000000000000000};
662:         PetscCall(PetscArraycpy(x, x_10, 30));

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

680: /* Compute Jacobians */
681: PetscErrorCode ComputeJacobian_Internal(const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal *quad, PetscReal *phypts, PetscReal *jacobian, PetscReal *ijacobian, PetscReal *dvolume)
682: {
683:   PetscInt  i;
684:   PetscReal volume = 1.0;

686:   PetscFunctionBegin;
690:   PetscCall(PetscArrayzero(jacobian, dim * dim));
691:   if (ijacobian) PetscCall(PetscArrayzero(ijacobian, dim * dim));
692:   if (phypts) PetscCall(PetscArrayzero(phypts, /*npts=1 * */ 3));

694:   if (dim == 1) {
695:     const PetscReal &r = quad[0];
696:     if (nverts == 2) { /* Linear Edge */
697:       const PetscReal dNi_dxi[2] = {-1.0, 1.0};

699:       for (i = 0; i < nverts; ++i) {
700:         const PetscReal *vertices = coordinates + i * 3;
701:         jacobian[0] += dNi_dxi[i] * vertices[0];
702:       }
703:     } else if (nverts == 3) { /* Quadratic Edge */
704:       const PetscReal dNi_dxi[3] = {4 * r - 3.0, 4 * (1.0 - 2.0 * r), 4.0 * r - 1.0};

706:       for (i = 0; i < nverts; ++i) {
707:         const PetscReal *vertices = coordinates + i * 3;
708:         jacobian[0] += dNi_dxi[i] * vertices[0];
709:       }
710:     } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "The number of 1-D entity vertices are invalid. Currently only support EDGE2 and EDGE3 basis evaluations in 1-D : %" PetscInt_FMT, nverts);

712:     if (ijacobian) {
713:       /* invert the jacobian */
714:       ijacobian[0] = 1.0 / jacobian[0];
715:     }
716:   } else if (dim == 2) {
717:     if (nverts == 4) { /* Linear Quadrangle */
718:       const PetscReal &r = quad[0];
719:       const PetscReal &s = quad[1];

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

724:       for (i = 0; i < nverts; ++i) {
725:         const PetscReal *vertices = coordinates + i * 3;
726:         jacobian[0] += dNi_dxi[i] * vertices[0];
727:         jacobian[2] += dNi_dxi[i] * vertices[1];
728:         jacobian[1] += dNi_deta[i] * vertices[0];
729:         jacobian[3] += dNi_deta[i] * vertices[1];
730:       }
731:     } else if (nverts == 3) { /* Linear triangle */
732:       const PetscReal x2 = coordinates[2 * 3 + 0], y2 = coordinates[2 * 3 + 1];

734:       /* Jacobian is constant */
735:       jacobian[0] = (coordinates[0 * 3 + 0] - x2);
736:       jacobian[1] = (coordinates[1 * 3 + 0] - x2);
737:       jacobian[2] = (coordinates[0 * 3 + 1] - y2);
738:       jacobian[3] = (coordinates[1 * 3 + 1] - y2);
739:     } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "The number of 2-D entity vertices are invalid. Currently only support QUAD4 and TRI3 basis evaluations in 2-D : %" PetscInt_FMT, nverts);

741:     /* invert the jacobian */
742:     if (ijacobian) PetscCall(DMatrix_Invert_2x2_Internal(jacobian, ijacobian, &volume));
743:   } else {
744:     if (nverts == 8) { /* Linear Hexahedra */
745:       const PetscReal &r          = quad[0];
746:       const PetscReal &s          = quad[1];
747:       const PetscReal &t          = quad[2];
748:       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)};

750:       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)};

752:       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)};

754:       for (i = 0; i < nverts; ++i) {
755:         const PetscReal *vertex = coordinates + i * 3;
756:         jacobian[0] += dNi_dxi[i] * vertex[0];
757:         jacobian[3] += dNi_dxi[i] * vertex[1];
758:         jacobian[6] += dNi_dxi[i] * vertex[2];
759:         jacobian[1] += dNi_deta[i] * vertex[0];
760:         jacobian[4] += dNi_deta[i] * vertex[1];
761:         jacobian[7] += dNi_deta[i] * vertex[2];
762:         jacobian[2] += dNi_dzeta[i] * vertex[0];
763:         jacobian[5] += dNi_dzeta[i] * vertex[1];
764:         jacobian[8] += dNi_dzeta[i] * vertex[2];
765:       }
766:     } else if (nverts == 4) { /* Linear Tetrahedra */
767:       const PetscReal x0 = coordinates[/*0 * 3 +*/ 0], y0 = coordinates[/*0 * 3 +*/ 1], z0 = coordinates[/*0 * 3 +*/ 2];

769:       /* compute the jacobian */
770:       jacobian[0] = coordinates[1 * 3 + 0] - x0;
771:       jacobian[1] = coordinates[2 * 3 + 0] - x0;
772:       jacobian[2] = coordinates[3 * 3 + 0] - x0;
773:       jacobian[3] = coordinates[1 * 3 + 1] - y0;
774:       jacobian[4] = coordinates[2 * 3 + 1] - y0;
775:       jacobian[5] = coordinates[3 * 3 + 1] - y0;
776:       jacobian[6] = coordinates[1 * 3 + 2] - z0;
777:       jacobian[7] = coordinates[2 * 3 + 2] - z0;
778:       jacobian[8] = coordinates[3 * 3 + 2] - z0;
779:     } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "The number of 3-D entity vertices are invalid. Currently only support HEX8 and TET4 basis evaluations in 3-D : %" PetscInt_FMT, nverts);

781:     if (ijacobian) {
782:       /* invert the jacobian */
783:       PetscCall(DMatrix_Invert_3x3_Internal(jacobian, ijacobian, &volume));
784:     }
785:   }
786:   PetscCheck(volume >= 1e-12, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Element has zero volume: %g. Degenerate element or invalid connectivity", volume);
787:   if (dvolume) *dvolume = volume;
788:   PetscFunctionReturn(PETSC_SUCCESS);
789: }

791: 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)
792: {
793:   PetscFunctionBegin;
794:   switch (dim) {
795:   case 1:
796:     PetscCall(Compute_Lagrange_Basis_1D_Internal(nverts, coordinates, 1, quadrature, phypts, NULL, phibasis, NULL, jacobian, ijacobian, volume));
797:     break;
798:   case 2:
799:     PetscCall(Compute_Lagrange_Basis_2D_Internal(nverts, coordinates, 1, quadrature, phypts, NULL, phibasis, NULL, NULL, jacobian, ijacobian, volume));
800:     break;
801:   case 3:
802:     PetscCall(Compute_Lagrange_Basis_3D_Internal(nverts, coordinates, 1, quadrature, phypts, NULL, phibasis, NULL, NULL, NULL, jacobian, ijacobian, volume));
803:     break;
804:   default:
805:     SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %" PetscInt_FMT, dim);
806:   }
807:   PetscFunctionReturn(PETSC_SUCCESS);
808: }

810: /*@C
811:   DMMoabPToRMapping - Compute the mapping from the physical coordinate system for a given element to the
812:   canonical reference element. In addition to finding the inverse mapping evaluation through Newton iteration,
813:   the basis function at the parametric point is also evaluated optionally.

815:   Input Parameters:
816: +  PetscInt  dim -         the element dimension (1=EDGE, 2=QUAD/TRI, 3=HEX/TET)
817: .  PetscInt nverts -       the number of vertices in the physical element
818: .  PetscReal coordinates - the coordinates of vertices in the physical element
819: -  PetscReal[3] xphy -     the coordinates of physical point for which natural coordinates (in reference frame) are sought

821:   Output Parameters:
822: +  PetscReal[3] natparam - the natural coordinates (in reference frame) corresponding to xphy
823: -  PetscReal[nverts] phi - the basis functions evaluated at the natural coordinates (natparam)

825:   Level: advanced

827: @*/
828: PetscErrorCode DMMoabPToRMapping(const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal *xphy, PetscReal *natparam, PetscReal *phi)
829: {
830:   /* Perform inverse evaluation for the mapping with use of Newton Raphson iteration */
831:   const PetscReal tol            = 1.0e-10;
832:   const PetscInt  max_iterations = 10;
833:   const PetscReal error_tol_sqr  = tol * tol;
834:   PetscReal       phibasis[8], jacobian[9], ijacobian[9], volume;
835:   PetscReal       phypts[3] = {0.0, 0.0, 0.0};
836:   PetscReal       delta[3]  = {0.0, 0.0, 0.0};
837:   PetscInt        iters     = 0;
838:   PetscReal       error     = 1.0;

840:   PetscFunctionBegin;

845:   PetscCall(PetscArrayzero(jacobian, dim * dim));
846:   PetscCall(PetscArrayzero(ijacobian, dim * dim));
847:   PetscCall(PetscArrayzero(phibasis, nverts));

849:   /* zero initial guess */
850:   natparam[0] = natparam[1] = natparam[2] = 0.0;

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

855:   error = 0.0;
856:   switch (dim) {
857:   case 3:
858:     delta[2] = phypts[2] - xphy[2];
859:     error += (delta[2] * delta[2]);
860:   case 2:
861:     delta[1] = phypts[1] - xphy[1];
862:     error += (delta[1] * delta[1]);
863:   case 1:
864:     delta[0] = phypts[0] - xphy[0];
865:     error += (delta[0] * delta[0]);
866:     break;
867:   }

869:   while (error > error_tol_sqr) {
870:     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]);

872:     /* Compute natparam -= J.inverse() * delta */
873:     switch (dim) {
874:     case 1:
875:       natparam[0] -= ijacobian[0] * delta[0];
876:       break;
877:     case 2:
878:       natparam[0] -= ijacobian[0] * delta[0] + ijacobian[1] * delta[1];
879:       natparam[1] -= ijacobian[2] * delta[0] + ijacobian[3] * delta[1];
880:       break;
881:     case 3:
882:       natparam[0] -= ijacobian[0] * delta[0] + ijacobian[1] * delta[1] + ijacobian[2] * delta[2];
883:       natparam[1] -= ijacobian[3] * delta[0] + ijacobian[4] * delta[1] + ijacobian[5] * delta[2];
884:       natparam[2] -= ijacobian[6] * delta[0] + ijacobian[7] * delta[1] + ijacobian[8] * delta[2];
885:       break;
886:     }

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

891:     error = 0.0;
892:     switch (dim) {
893:     case 3:
894:       delta[2] = phypts[2] - xphy[2];
895:       error += (delta[2] * delta[2]);
896:     case 2:
897:       delta[1] = phypts[1] - xphy[1];
898:       error += (delta[1] * delta[1]);
899:     case 1:
900:       delta[0] = phypts[0] - xphy[0];
901:       error += (delta[0] * delta[0]);
902:       break;
903:     }
904:   }
905:   if (phi) PetscCall(PetscArraycpy(phi, phibasis, nverts));
906:   PetscFunctionReturn(PETSC_SUCCESS);
907: }