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 0;
 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])
 28:            - inmat[3] * (inmat[8] * inmat[1] - inmat[7] * inmat[2])
 29:            + inmat[6] * (inmat[5] * inmat[1] - inmat[4] * inmat[2]);
 30: }

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

 50: inline PetscReal DMatrix_Determinant_4x4_Internal(PetscReal inmat[4 * 4])
 51: {
 52:   return
 53:     inmat[0 + 0 * 4] * (
 54:       inmat[1 + 1 * 4] * (inmat[2 + 2 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 2 * 4])
 55:       - inmat[1 + 2 * 4] * (inmat[2 + 1 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 1 * 4])
 56:       + inmat[1 + 3 * 4] * (inmat[2 + 1 * 4] * inmat[3 + 2 * 4] - inmat[2 + 2 * 4] * inmat[3 + 1 * 4]))
 57:     - inmat[0 + 1 * 4] * (
 58:       inmat[1 + 0 * 4] * (inmat[2 + 2 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 2 * 4])
 59:       - inmat[1 + 2 * 4] * (inmat[2 + 0 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 0 * 4])
 60:       + inmat[1 + 3 * 4] * (inmat[2 + 0 * 4] * inmat[3 + 2 * 4] - inmat[2 + 2 * 4] * inmat[3 + 0 * 4]))
 61:     + inmat[0 + 2 * 4] * (
 62:       inmat[1 + 0 * 4] * (inmat[2 + 1 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 1 * 4])
 63:       - inmat[1 + 1 * 4] * (inmat[2 + 0 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 0 * 4])
 64:       + inmat[1 + 3 * 4] * (inmat[2 + 0 * 4] * inmat[3 + 1 * 4] - inmat[2 + 1 * 4] * inmat[3 + 0 * 4]))
 65:     - inmat[0 + 3 * 4] * (
 66:       inmat[1 + 0 * 4] * (inmat[2 + 1 * 4] * inmat[3 + 2 * 4] - inmat[2 + 2 * 4] * inmat[3 + 1 * 4])
 67:       - inmat[1 + 1 * 4] * (inmat[2 + 0 * 4] * inmat[3 + 2 * 4] - inmat[2 + 2 * 4] * inmat[3 + 0 * 4])
 68:       + inmat[1 + 2 * 4] * (inmat[2 + 0 * 4] * inmat[3 + 1 * 4] - inmat[2 + 1 * 4] * inmat[3 + 0 * 4]));
 69: }

 71: inline PetscErrorCode DMatrix_Invert_4x4_Internal (PetscReal *inmat, PetscReal *outmat, PetscScalar *determinant)
 72: {
 73:   PetscReal det = DMatrix_Determinant_4x4_Internal(inmat);
 74:   if (outmat) {
 75:     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;
 76:     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;
 77:     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;
 78:     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;
 79:     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;
 80:     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;
 81:     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;
 82:     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;
 83:     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;
 84:     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;
 85:     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;
 86:     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;
 87:     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;
 88:     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;
 89:     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;
 90:     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;
 91:   }
 92:   if (determinant) *determinant = det;
 93:   return 0;
 94: }

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

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

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

104:   Notes:

106:   Example Physical Element
107: .vb
108:     1-------2        1----3----2
109:       EDGE2             EDGE3
110: .ve

112:   Input Parameters:
113: +  PetscInt  nverts -          the number of element vertices
114: .  PetscReal coords[3*nverts] - the physical coordinates of the vertices (in canonical numbering)
115: .  PetscInt  npts -            the number of evaluation points (quadrature points)
116: -  PetscReal quad[3*npts] -    the evaluation points (quadrature points) in the reference space

118:   Output Parameters:
119: +  PetscReal phypts[3*npts] -  the evaluation points (quadrature points) transformed to the physical space
120: .  PetscReal jxw[npts] -       the jacobian determinant * quadrature weight necessary for assembling discrete contributions
121: .  PetscReal phi[npts] -       the bases evaluated at the specified quadrature points
122: .  PetscReal dphidx[npts] -    the derivative of the bases wrt X-direction evaluated at the specified quadrature points
123: .  jacobian -
124: .  ijacobian -
125: -  volume

127:   Level: advanced

129: @*/
130: PetscErrorCode Compute_Lagrange_Basis_1D_Internal(const PetscInt nverts, const PetscReal *coords, const PetscInt npts, const PetscReal *quad, PetscReal *phypts,
131:                                                   PetscReal *jxw, PetscReal *phi, PetscReal *dphidx,
132:                                                   PetscReal *jacobian, PetscReal *ijacobian, PetscReal *volume)
133: {
134:   int             i, j;
135:   PetscErrorCode  ierr;

141:   if (phypts) {
142:     PetscArrayzero(phypts, npts * 3);
143:   }
144:   if (dphidx) { /* Reset arrays. */
145:     PetscArrayzero(dphidx, npts * nverts);
146:   }
147:   if (nverts == 2) { /* Linear Edge */

149:     for (j = 0; j < npts; j++) {
150:       const PetscInt offset = j * nverts;
151:       const PetscReal r = quad[j];

153:       phi[0 + offset] = ( 1.0 - r);
154:       phi[1 + offset] = (       r);

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

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

167:       /* invert the jacobian */
168:       *volume = jacobian[0];
169:       ijacobian[0] = 1.0 / jacobian[0];
170:       jxw[j] *= *volume;

172:       /*  Divide by element jacobian. */
173:       for (i = 0; i < nverts; i++) {
174:         if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[0];
175:       }
176:     }
177:   } else if (nverts == 3) { /* Quadratic Edge */

179:     for (j = 0; j < npts; j++) {
180:       const PetscInt offset = j * nverts;
181:       const PetscReal r = quad[j];

183:       phi[0 + offset] = 1.0 + r * ( 2.0 * r - 3.0);
184:       phi[1 + offset] = 4.0 * r * ( 1.0 - r);
185:       phi[2 + offset] = r * ( 2.0 * r - 1.0);

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

189:       jacobian[0] = ijacobian[0] = volume[0] = 0.0;
190:       for (i = 0; i < nverts; ++i) {
191:         const PetscReal* vertices = coords + i * 3;
192:         jacobian[0] += dNi_dxi[i] * vertices[0];
193:         if (phypts) {
194:           phypts[3 * j + 0] += phi[i + offset] * vertices[0];
195:         }
196:       }

198:       /* invert the jacobian */
199:       *volume = jacobian[0];
200:       ijacobian[0] = 1.0 / jacobian[0];
201:       if (jxw) jxw[j] *= *volume;

203:       /*  Divide by element jacobian. */
204:       for (i = 0; i < nverts; i++) {
205:         if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[0];
206:       }
207:     }
208:   } else SETERRQ1(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 : %D", nverts);
209:   return(0);
210: }

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

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

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

220:   Notes:

222:   Example Physical Element (QUAD4)
223: .vb
224:     4------3        s
225:     |      |        |
226:     |      |        |
227:     |      |        |
228:     1------2        0-------r
229: .ve

231:   Input Parameters:
232: +  PetscInt  nverts -          the number of element vertices
233: .  PetscReal coords[3*nverts] - the physical coordinates of the vertices (in canonical numbering)
234: .  PetscInt  npts -            the number of evaluation points (quadrature points)
235: -  PetscReal quad[3*npts] -    the evaluation points (quadrature points) in the reference space

237:   Output Parameters:
238: +  PetscReal phypts[3*npts] -  the evaluation points (quadrature points) transformed to the physical space
239: .  PetscReal jxw[npts] -       the jacobian determinant * quadrature weight necessary for assembling discrete contributions
240: .  PetscReal phi[npts] -       the bases evaluated at the specified quadrature points
241: .  PetscReal dphidx[npts] -    the derivative of the bases wrt X-direction evaluated at the specified quadrature points
242: .  PetscReal dphidy[npts] -    the derivative of the bases wrt Y-direction evaluated at the specified quadrature points
243: .  jacobian -
244: .  ijacobian -
245: -  volume

247:   Level: advanced

249: @*/
250: PetscErrorCode Compute_Lagrange_Basis_2D_Internal(const PetscInt nverts, const PetscReal *coords, const PetscInt npts, const PetscReal *quad, PetscReal *phypts,
251:                                                   PetscReal *jxw, PetscReal *phi, PetscReal *dphidx, PetscReal *dphidy,
252:                                                   PetscReal *jacobian, PetscReal *ijacobian, PetscReal *volume)
253: {
254:   PetscInt       i, j, k;

261:   PetscArrayzero(phi, npts);
262:   if (phypts) {
263:     PetscArrayzero(phypts, npts * 3);
264:   }
265:   if (dphidx) { /* Reset arrays. */
266:     PetscArrayzero(dphidx, npts * nverts);
267:     PetscArrayzero(dphidy, npts * nverts);
268:   }
269:   if (nverts == 4) { /* Linear Quadrangle */

271:     for (j = 0; j < npts; j++)
272:     {
273:       const PetscInt offset = j * nverts;
274:       const PetscReal r = quad[0 + j * 2];
275:       const PetscReal s = quad[1 + j * 2];

277:       phi[0 + offset] = ( 1.0 - r) * ( 1.0 - s);
278:       phi[1 + offset] =         r   * ( 1.0 - s);
279:       phi[2 + offset] =         r   *         s;
280:       phi[3 + offset] = ( 1.0 - r) *         s;

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

285:       PetscArrayzero(jacobian, 4);
286:       PetscArrayzero(ijacobian, 4);
287:       for (i = 0; i < nverts; ++i) {
288:         const PetscReal* vertices = coords + i * 3;
289:         jacobian[0] += dNi_dxi[i] * vertices[0];
290:         jacobian[2] += dNi_dxi[i] * vertices[1];
291:         jacobian[1] += dNi_deta[i] * vertices[0];
292:         jacobian[3] += dNi_deta[i] * vertices[1];
293:         if (phypts) {
294:           phypts[3 * j + 0] += phi[i + offset] * vertices[0];
295:           phypts[3 * j + 1] += phi[i + offset] * vertices[1];
296:           phypts[3 * j + 2] += phi[i + offset] * vertices[2];
297:         }
298:       }

300:       /* invert the jacobian */
301:       DMatrix_Invert_2x2_Internal(jacobian, ijacobian, volume);
302:       if (*volume < 1e-12) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Quadrangular element has zero volume: %g. Degenerate element or invalid connectivity\n", *volume);

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

306:       /*  Let us compute the bases derivatives by scaling with inverse jacobian. */
307:       if (dphidx) {
308:         for (i = 0; i < nverts; i++) {
309:           for (k = 0; k < 2; ++k) {
310:             if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[k * 2 + 0];
311:             if (dphidy) dphidy[i + offset] += dNi_deta[i] * ijacobian[k * 2 + 1];
312:           } /* for k=[0..2) */
313:         } /* for i=[0..nverts) */
314:       } /* if (dphidx) */
315:     }
316:   } else if (nverts == 3) { /* Linear triangle */
317:     const PetscReal x2 = coords[2 * 3 + 0], y2 = coords[2 * 3 + 1];

319:     PetscArrayzero(jacobian, 4);
320:     PetscArrayzero(ijacobian, 4);

322:     /* Jacobian is constant */
323:     jacobian[0] = (coords[0 * 3 + 0] - x2); jacobian[1] = (coords[1 * 3 + 0] - x2);
324:     jacobian[2] = (coords[0 * 3 + 1] - y2); jacobian[3] = (coords[1 * 3 + 1] - y2);

326:     /* invert the jacobian */
327:     DMatrix_Invert_2x2_Internal(jacobian, ijacobian, volume);
328:     if (*volume < 1e-12) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Triangular element has zero volume: %g. Degenerate element or invalid connectivity\n", (double)*volume);

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

333:     for (j = 0; j < npts; j++) {
334:       const PetscInt offset = j * nverts;
335:       const PetscReal r = quad[0 + j * 2];
336:       const PetscReal s = quad[1 + j * 2];

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

340:       /* const PetscReal Ni[3]  = { r, s, 1.0 - r - s }; */
341:       const PetscReal phipts_x = coords[2 * 3 + 0] + jacobian[0] * r + jacobian[1] * s;
342:       const PetscReal phipts_y = coords[2 * 3 + 1] + jacobian[2] * r + jacobian[3] * s;
343:       if (phypts) {
344:         phypts[offset + 0] = phipts_x;
345:         phypts[offset + 1] = phipts_y;
346:       }

348:       /* \phi_0 = (b.y - c.y) x + (b.x - c.x) y + c.x b.y - b.x c.y */
349:       phi[0 + offset] = (  ijacobian[0] * (phipts_x - x2) + ijacobian[1] * (phipts_y - y2));
350:       /* \phi_1 = (c.y - a.y) x + (a.x - c.x) y + c.x a.y - a.x c.y */
351:       phi[1 + offset] = (  ijacobian[2] * (phipts_x - x2) + ijacobian[3] * (phipts_y - y2));
352:       phi[2 + offset] = 1.0 - phi[0 + offset] - phi[1 + offset];

354:       if (dphidx) {
355:         dphidx[0 + offset] = Dx[0];
356:         dphidx[1 + offset] = Dx[1];
357:         dphidx[2 + offset] = Dx[2];
358:       }

360:       if (dphidy) {
361:         dphidy[0 + offset] = Dy[0];
362:         dphidy[1 + offset] = Dy[1];
363:         dphidy[2 + offset] = Dy[2];
364:       }

366:     }
367:   } else SETERRQ1(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 : %D", nverts);
368:   return(0);
369: }

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

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

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

379:   Notes:

381:   Example Physical Element (HEX8)
382: .vb
383:       8------7
384:      /|     /|        t  s
385:     5------6 |        | /
386:     | |    | |        |/
387:     | 4----|-3        0-------r
388:     |/     |/
389:     1------2
390: .ve

392:   Input Parameters:
393: +  PetscInt  nverts -          the number of element vertices
394: .  PetscReal coords[3*nverts] - the physical coordinates of the vertices (in canonical numbering)
395: .  PetscInt  npts -            the number of evaluation points (quadrature points)
396: -  PetscReal quad[3*npts] -    the evaluation points (quadrature points) in the reference space

398:   Output Parameters:
399: +  PetscReal phypts[3*npts] -  the evaluation points (quadrature points) transformed to the physical space
400: .  PetscReal jxw[npts] -       the jacobian determinant * quadrature weight necessary for assembling discrete contributions
401: .  PetscReal phi[npts] -       the bases evaluated at the specified quadrature points
402: .  PetscReal dphidx[npts] -    the derivative of the bases wrt X-direction evaluated at the specified quadrature points
403: .  PetscReal dphidy[npts] -    the derivative of the bases wrt Y-direction evaluated at the specified quadrature points
404: .  PetscReal dphidz[npts] -    the derivative of the bases wrt Z-direction evaluated at the specified quadrature points
405: .  jacobian -
406: .  ijacobian -
407: -  volume

409:   Level: advanced

411: @*/
412: PetscErrorCode Compute_Lagrange_Basis_3D_Internal(const PetscInt nverts, const PetscReal *coords, const PetscInt npts, const PetscReal *quad, PetscReal *phypts,
413:                                                   PetscReal *jxw, PetscReal *phi, PetscReal *dphidx, PetscReal *dphidy, PetscReal *dphidz,
414:                                                   PetscReal *jacobian, PetscReal *ijacobian, PetscReal *volume)
415: {
416:   PetscInt       i, j, k;


424:   PetscArrayzero(phi, npts);
425:   if (phypts) {
426:     PetscArrayzero(phypts, npts * 3);
427:   }
428:   if (dphidx) {
429:     PetscArrayzero(dphidx, npts * nverts);
430:     PetscArrayzero(dphidy, npts * nverts);
431:     PetscArrayzero(dphidz, npts * nverts);
432:   }

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

436:     for (j = 0; j < npts; j++) {
437:       const PetscInt offset = j * nverts;
438:       const PetscReal& r = quad[j * 3 + 0];
439:       const PetscReal& s = quad[j * 3 + 1];
440:       const PetscReal& t = quad[j * 3 + 2];

442:       phi[offset + 0] = ( 1.0 - r) * ( 1.0 - s) * ( 1.0 - t); /* 0,0,0 */
443:       phi[offset + 1] = (       r) * ( 1.0 - s) * ( 1.0 - t); /* 1,0,0 */
444:       phi[offset + 2] = (       r) * (       s) * ( 1.0 - t); /* 1,1,0 */
445:       phi[offset + 3] = ( 1.0 - r) * (       s) * ( 1.0 - t); /* 0,1,0 */
446:       phi[offset + 4] = ( 1.0 - r) * ( 1.0 - s) * (       t); /* 0,0,1 */
447:       phi[offset + 5] = (       r) * ( 1.0 - s) * (       t); /* 1,0,1 */
448:       phi[offset + 6] = (       r) * (       s) * (       t); /* 1,1,1 */
449:       phi[offset + 7] = ( 1.0 - r) * (       s) * (       t); /* 0,1,1 */

451:       const PetscReal dNi_dxi[8]  = {- ( 1.0 - s) * ( 1.0 - t),
452:                                        ( 1.0 - s) * ( 1.0 - t),
453:                                        (       s) * ( 1.0 - t),
454:                                      - (       s) * ( 1.0 - t),
455:                                      - ( 1.0 - s) * (       t),
456:                                        ( 1.0 - s) * (       t),
457:                                        (       s) * (       t),
458:                                      - (       s) * (       t)
459:                                     };

461:       const PetscReal dNi_deta[8]  = { - ( 1.0 - r) * ( 1.0 - t),
462:                                        - (       r) * ( 1.0 - t),
463:                                          (       r) * ( 1.0 - t),
464:                                          ( 1.0 - r) * ( 1.0 - t),
465:                                        - ( 1.0 - r) * (       t),
466:                                        - (       r) * (       t),
467:                                          (       r) * (       t),
468:                                          ( 1.0 - r) * (       t)
469:                                       };

471:       const PetscReal dNi_dzeta[8]  = { - ( 1.0 - r) * ( 1.0 - s),
472:                                         - (       r) * ( 1.0 - s),
473:                                         - (       r) * (       s),
474:                                         - ( 1.0 - r) * (       s),
475:                                           ( 1.0 - r) * ( 1.0 - s),
476:                                           (       r) * ( 1.0 - s),
477:                                           (       r) * (       s),
478:                                           ( 1.0 - r) * (       s)
479:                                      };

481:       PetscArrayzero(jacobian, 9);
482:       PetscArrayzero(ijacobian, 9);
483:       for (i = 0; i < nverts; ++i) {
484:         const PetscReal* vertex = coords + i * 3;
485:         jacobian[0] += dNi_dxi[i]   * vertex[0];
486:         jacobian[3] += dNi_dxi[i]   * vertex[1];
487:         jacobian[6] += dNi_dxi[i]   * vertex[2];
488:         jacobian[1] += dNi_deta[i]  * vertex[0];
489:         jacobian[4] += dNi_deta[i]  * vertex[1];
490:         jacobian[7] += dNi_deta[i]  * vertex[2];
491:         jacobian[2] += dNi_dzeta[i] * vertex[0];
492:         jacobian[5] += dNi_dzeta[i] * vertex[1];
493:         jacobian[8] += dNi_dzeta[i] * vertex[2];
494:         if (phypts) {
495:           phypts[3 * j + 0] += phi[i + offset] * vertex[0];
496:           phypts[3 * j + 1] += phi[i + offset] * vertex[1];
497:           phypts[3 * j + 2] += phi[i + offset] * vertex[2];
498:         }
499:       }

501:       /* invert the jacobian */
502:       DMatrix_Invert_3x3_Internal(jacobian, ijacobian, volume);
503:       if (*volume < 1e-8) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Hexahedral element has zero volume: %g. Degenerate element or invalid connectivity\n", *volume);

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

507:       /*  Divide by element jacobian. */
508:       for (i = 0; i < nverts; ++i) {
509:         for (k = 0; k < 3; ++k) {
510:           if (dphidx) dphidx[i + offset] += dNi_dxi[i]   * ijacobian[0 * 3 + k];
511:           if (dphidy) dphidy[i + offset] += dNi_deta[i]  * ijacobian[1 * 3 + k];
512:           if (dphidz) dphidz[i + offset] += dNi_dzeta[i] * ijacobian[2 * 3 + k];
513:         }
514:       }
515:     }
516:   } else if (nverts == 4) { /* Linear Tetrahedra */
517:     PetscReal       Dx[4]={0,0,0,0}, Dy[4]={0,0,0,0}, Dz[4]={0,0,0,0};
518:     const PetscReal x0 = coords[/*0 * 3 +*/ 0], y0 = coords[/*0 * 3 +*/ 1], z0 = coords[/*0 * 3 +*/ 2];

520:     PetscArrayzero(jacobian, 9);
521:     PetscArrayzero(ijacobian, 9);

523:     /* compute the jacobian */
524:     jacobian[0] = coords[1 * 3 + 0] - x0;  jacobian[1] = coords[2 * 3 + 0] - x0; jacobian[2] = coords[3 * 3 + 0] - x0;
525:     jacobian[3] = coords[1 * 3 + 1] - y0;  jacobian[4] = coords[2 * 3 + 1] - y0; jacobian[5] = coords[3 * 3 + 1] - y0;
526:     jacobian[6] = coords[1 * 3 + 2] - z0;  jacobian[7] = coords[2 * 3 + 2] - z0; jacobian[8] = coords[3 * 3 + 2] - z0;

528:     /* invert the jacobian */
529:     DMatrix_Invert_3x3_Internal(jacobian, ijacobian, volume);
530:     if (*volume < 1e-8) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Tetrahedral element has zero volume: %g. Degenerate element or invalid connectivity\n", (double)*volume);

532:     if (dphidx) {
533:       Dx[0] =   ( coords[1 + 2 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3])
534:                  - coords[1 + 1 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3])
535:                  - coords[1 + 3 * 3] * ( coords[2 + 1 * 3] - coords[2 + 2 * 3])) / *volume;
536:       Dx[1] = - ( coords[1 + 2 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3])
537:                  - coords[1 + 0 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3])
538:                  - coords[1 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 2 * 3])) / *volume;
539:       Dx[2] =   ( coords[1 + 1 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3])
540:                  - coords[1 + 0 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3])
541:                  - coords[1 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 1 * 3])) / *volume;
542:       Dx[3] =  - ( Dx[0] + Dx[1] + Dx[2]);
543:       Dy[0] =   ( coords[0 + 1 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3])
544:                  - coords[0 + 2 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3])
545:                  + coords[0 + 3 * 3] * ( coords[2 + 1 * 3] - coords[2 + 2 * 3])) / *volume;
546:       Dy[1] = - ( coords[0 + 0 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3])
547:                  - coords[0 + 2 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3])
548:                  + coords[0 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 2 * 3])) / *volume;
549:       Dy[2] =   ( coords[0 + 0 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3])
550:                  - coords[0 + 1 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3])
551:                  + coords[0 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 1 * 3])) / *volume;
552:       Dy[3] =  - ( Dy[0] + Dy[1] + Dy[2]);
553:       Dz[0] =   ( coords[0 + 1 * 3] * (coords[1 + 3 * 3] - coords[1 + 2 * 3])
554:                  - coords[0 + 2 * 3] * (coords[1 + 3 * 3] - coords[1 + 1 * 3])
555:                  + coords[0 + 3 * 3] * (coords[1 + 2 * 3] - coords[1 + 1 * 3])) / *volume;
556:       Dz[1] = - ( coords[0 + 0 * 3] * (coords[1 + 3 * 3] - coords[1 + 2 * 3])
557:                   + coords[0 + 2 * 3] * (coords[1 + 0 * 3] - coords[1 + 3 * 3])
558:                   - coords[0 + 3 * 3] * (coords[1 + 0 * 3] - coords[1 + 2 * 3])) / *volume;
559:       Dz[2] =   ( coords[0 + 0 * 3] * (coords[1 + 3 * 3] - coords[1 + 1 * 3])
560:                  + coords[0 + 1 * 3] * (coords[1 + 0 * 3] - coords[1 + 3 * 3])
561:                  - coords[0 + 3 * 3] * (coords[1 + 0 * 3] - coords[1 + 1 * 3])) / *volume;
562:       Dz[3] =  - ( Dz[0] + Dz[1] + Dz[2]);
563:     }

565:     for (j = 0; j < npts; j++) {
566:       const PetscInt offset = j * nverts;
567:       const PetscReal& r = quad[j * 3 + 0];
568:       const PetscReal& s = quad[j * 3 + 1];
569:       const PetscReal& t = quad[j * 3 + 2];

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

573:       phi[offset + 0] = 1.0 - r - s - t;
574:       phi[offset + 1] = r;
575:       phi[offset + 2] = s;
576:       phi[offset + 3] = t;

578:       if (phypts) {
579:         for (i = 0; i < nverts; ++i) {
580:           const PetscScalar* vertices = coords + i * 3;
581:           phypts[3 * j + 0] += phi[i + offset] * vertices[0];
582:           phypts[3 * j + 1] += phi[i + offset] * vertices[1];
583:           phypts[3 * j + 2] += phi[i + offset] * vertices[2];
584:         }
585:       }

587:       /* Now set the derivatives */
588:       if (dphidx) {
589:         dphidx[0 + offset] = Dx[0];
590:         dphidx[1 + offset] = Dx[1];
591:         dphidx[2 + offset] = Dx[2];
592:         dphidx[3 + offset] = Dx[3];

594:         dphidy[0 + offset] = Dy[0];
595:         dphidy[1 + offset] = Dy[1];
596:         dphidy[2 + offset] = Dy[2];
597:         dphidy[3 + offset] = Dy[3];

599:         dphidz[0 + offset] = Dz[0];
600:         dphidz[1 + offset] = Dz[1];
601:         dphidz[2 + offset] = Dz[2];
602:         dphidz[3 + offset] = Dz[3];
603:       }

605:     } /* Tetrahedra -- ends */
606:   } else SETERRQ1(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 : %D", nverts);
607:   return(0);
608: }

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

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

616:   Input Parameters:
617: +  PetscInt  nverts,           the number of element vertices
618: .  PetscReal coords[3*nverts], the physical coordinates of the vertices (in canonical numbering)
619: .  PetscInt  npts,             the number of evaluation points (quadrature points)
620: -  PetscReal quad[3*npts],     the evaluation points (quadrature points) in the reference space

622:   Output Parameters:
623: +  PetscReal phypts[3*npts],   the evaluation points (quadrature points) transformed to the physical space
624: .  PetscReal jxw[npts],        the jacobian determinant * quadrature weight necessary for assembling discrete contributions
625: .  PetscReal fe_basis[npts],        the bases values evaluated at the specified quadrature points
626: -  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

628:   Level: advanced

630: @*/
631: PetscErrorCode DMMoabFEMComputeBasis(const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscQuadrature quadrature,
632:                                      PetscReal *phypts, PetscReal *jacobian_quadrature_weight_product,
633:                                      PetscReal *fe_basis, PetscReal **fe_basis_derivatives)
634: {
635:   PetscErrorCode  ierr;
636:   PetscInt        npoints,idim;
637:   bool            compute_der;
638:   const PetscReal *quadpts, *quadwts;
639:   PetscReal       jacobian[9], ijacobian[9], volume;

645:   compute_der = (fe_basis_derivatives != NULL);

647:   /* Get the quadrature points and weights for the given quadrature rule */
648:   PetscQuadratureGetData(quadrature, &idim, NULL, &npoints, &quadpts, &quadwts);
649:   if (idim != dim) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Dimension mismatch: provided (%D) vs quadrature (%D)\n",idim,dim);
650:   if (jacobian_quadrature_weight_product) {
651:     PetscArraycpy(jacobian_quadrature_weight_product, quadwts, npoints);
652:   }

654:   switch (dim) {
655:   case 1:
656:     Compute_Lagrange_Basis_1D_Internal(nverts, coordinates, npoints, quadpts, phypts,
657:            jacobian_quadrature_weight_product, fe_basis,
658:            (compute_der ? fe_basis_derivatives[0] : NULL),
659:            jacobian, ijacobian, &volume);
660:     break;
661:   case 2:
662:     Compute_Lagrange_Basis_2D_Internal(nverts, coordinates, npoints, quadpts, phypts,
663:            jacobian_quadrature_weight_product, fe_basis,
664:            (compute_der ? fe_basis_derivatives[0] : NULL),
665:            (compute_der ? fe_basis_derivatives[1] : NULL),
666:            jacobian, ijacobian, &volume);
667:     break;
668:   case 3:
669:     Compute_Lagrange_Basis_3D_Internal(nverts, coordinates, npoints, quadpts, phypts,
670:            jacobian_quadrature_weight_product, fe_basis,
671:            (compute_der ? fe_basis_derivatives[0] : NULL),
672:            (compute_der ? fe_basis_derivatives[1] : NULL),
673:            (compute_der ? fe_basis_derivatives[2] : NULL),
674:            jacobian, ijacobian, &volume);
675:     break;
676:   default:
677:     SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %D", dim);
678:   }
679:   return(0);
680: }

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

686:   Input Parameters:

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

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

694:   Level: advanced

696: @*/
697: PetscErrorCode DMMoabFEMCreateQuadratureDefault(const PetscInt dim, const PetscInt nverts, PetscQuadrature *quadrature)
698: {
699:   PetscReal       *w, *x;
700:   PetscInt        nc=1;
701:   PetscErrorCode  ierr;

704:   /* Create an appropriate quadrature rule to sample basis */
705:   switch (dim)
706:   {
707:   case 1:
708:     /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */
709:     PetscDTStroudConicalQuadrature(1, nc, nverts, 0, 1.0, quadrature);
710:     break;
711:   case 2:
712:     /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */
713:     if (nverts == 3) {
714:       const PetscInt order = 2;
715:       const PetscInt npoints = (order == 2 ? 3 : 6);
716:       PetscMalloc2(npoints * 2, &x, npoints, &w);
717:       if (npoints == 3) {
718:         x[0] = x[1] = x[2] = x[5] = 1.0 / 6.0;
719:         x[3] = x[4] = 2.0 / 3.0;
720:         w[0] = w[1] = w[2] = 1.0 / 3.0;
721:       } else if (npoints == 6) {
722:         x[0] = x[1] = x[2] = 0.44594849091597;
723:         x[3] = x[4] = 0.10810301816807;
724:         x[5] = 0.44594849091597;
725:         x[6] = x[7] = x[8] = 0.09157621350977;
726:         x[9] = x[10] = 0.81684757298046;
727:         x[11] = 0.09157621350977;
728:         w[0] = w[1] = w[2] = 0.22338158967801;
729:         w[3] = w[4] = w[5] = 0.10995174365532;
730:       } else SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Triangle quadrature rules for points 3 and 6 supported; npoints : %D", npoints);
731:       PetscQuadratureCreate(PETSC_COMM_SELF, quadrature);
732:       PetscQuadratureSetOrder(*quadrature, order);
733:       PetscQuadratureSetData(*quadrature, dim, nc, npoints, x, w);
734:       /* PetscDTStroudConicalQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature); */
735:     } else {
736:       PetscDTGaussTensorQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature);
737:     }
738:     break;
739:   case 3:
740:     /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */
741:     if (nverts == 4) {
742:       PetscInt order;
743:       const PetscInt npoints = 4; // Choose between 4 and 10
744:       PetscMalloc2(npoints * 3, &x, npoints, &w);
745:       if (npoints == 4) { /*  KEAST1, K1,  N=4, O=4 */
746:         const PetscReal x_4[12] = { 0.5854101966249685, 0.1381966011250105, 0.1381966011250105,
747:                                     0.1381966011250105, 0.1381966011250105, 0.1381966011250105,
748:                                     0.1381966011250105, 0.1381966011250105, 0.5854101966249685,
749:                                     0.1381966011250105, 0.5854101966249685, 0.1381966011250105
750:                                   };
751:         PetscArraycpy(x, x_4, 12);

753:         w[0] = w[1] = w[2] = w[3] = 1.0 / 24.0;
754:         order = 4;
755:       } else if (npoints == 10) { /*  KEAST3, K3  N=10, O=10 */
756:         const PetscReal x_10[30] = { 0.5684305841968444, 0.1438564719343852, 0.1438564719343852,
757:                                      0.1438564719343852, 0.1438564719343852, 0.1438564719343852,
758:                                      0.1438564719343852, 0.1438564719343852, 0.5684305841968444,
759:                                      0.1438564719343852, 0.5684305841968444, 0.1438564719343852,
760:                                      0.0000000000000000, 0.5000000000000000, 0.5000000000000000,
761:                                      0.5000000000000000, 0.0000000000000000, 0.5000000000000000,
762:                                      0.5000000000000000, 0.5000000000000000, 0.0000000000000000,
763:                                      0.5000000000000000, 0.0000000000000000, 0.0000000000000000,
764:                                      0.0000000000000000, 0.5000000000000000, 0.0000000000000000,
765:                                      0.0000000000000000, 0.0000000000000000, 0.5000000000000000
766:                                    };
767:         PetscArraycpy(x, x_10, 30);

769:         w[0] = w[1] = w[2] = w[3] = 0.2177650698804054;
770:         w[4] = w[5] = w[6] = w[7] = w[8] = w[9] = 0.0214899534130631;
771:         order = 10;
772:       } else SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Tetrahedral quadrature rules for points 4 and 10 supported; npoints : %D", npoints);
773:       PetscQuadratureCreate(PETSC_COMM_SELF, quadrature);
774:       PetscQuadratureSetOrder(*quadrature, order);
775:       PetscQuadratureSetData(*quadrature, dim, nc, npoints, x, w);
776:       /* PetscDTStroudConicalQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature); */
777:     } else {
778:       PetscDTGaussTensorQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature);
779:     }
780:     break;
781:   default:
782:     SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %D", dim);
783:   }
784:   return(0);
785: }

787: /* Compute Jacobians */
788: PetscErrorCode ComputeJacobian_Internal (const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal *quad, PetscReal *phypts,
789:   PetscReal *jacobian, PetscReal *ijacobian, PetscReal* dvolume)
790: {
791:   PetscInt       i;
792:   PetscReal      volume=1.0;

799:   PetscArrayzero(jacobian, dim * dim);
800:   if (ijacobian) {
801:     PetscArrayzero(ijacobian, dim * dim);
802:   }
803:   if (phypts) {
804:     PetscArrayzero(phypts, /*npts=1 * */ 3);
805:   }

807:   if (dim == 1) {
808:     const PetscReal& r = quad[0];
809:     if (nverts == 2) { /* Linear Edge */
810:       const PetscReal dNi_dxi[2]  = { -1.0, 1.0 };

812:       for (i = 0; i < nverts; ++i) {
813:         const PetscReal* vertices = coordinates + i * 3;
814:         jacobian[0] += dNi_dxi[i] * vertices[0];
815:       }
816:     } else if (nverts == 3) { /* Quadratic Edge */
817:       const PetscReal dNi_dxi[3]  = { 4 * r - 3.0, 4 * ( 1.0 - 2.0 * r), 4.0 * r - 1.0};

819:       for (i = 0; i < nverts; ++i) {
820:         const PetscReal* vertices = coordinates + i * 3;
821:         jacobian[0] += dNi_dxi[i] * vertices[0];
822:       }
823:     } else {
824:       SETERRQ1(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 : %D", nverts);
825:     }

827:     if (ijacobian) {
828:       /* invert the jacobian */
829:       ijacobian[0] = 1.0 / jacobian[0];
830:     }
831:   } else if (dim == 2) {

833:     if (nverts == 4) { /* Linear Quadrangle */
834:       const PetscReal& r = quad[0];
835:       const PetscReal& s = quad[1];

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

840:       for (i = 0; i < nverts; ++i) {
841:         const PetscReal* vertices = coordinates + i * 3;
842:         jacobian[0] += dNi_dxi[i]  * vertices[0];
843:         jacobian[2] += dNi_dxi[i]  * vertices[1];
844:         jacobian[1] += dNi_deta[i] * vertices[0];
845:         jacobian[3] += dNi_deta[i] * vertices[1];
846:       }
847:     } else if (nverts == 3) { /* Linear triangle */
848:       const PetscReal x2 = coordinates[2 * 3 + 0], y2 = coordinates[2 * 3 + 1];

850:       /* Jacobian is constant */
851:       jacobian[0] = (coordinates[0 * 3 + 0] - x2); jacobian[1] = (coordinates[1 * 3 + 0] - x2);
852:       jacobian[2] = (coordinates[0 * 3 + 1] - y2); jacobian[3] = (coordinates[1 * 3 + 1] - y2);
853:     } else SETERRQ1(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 : %D", nverts);

855:     /* invert the jacobian */
856:     if (ijacobian) {
857:       DMatrix_Invert_2x2_Internal(jacobian, ijacobian, &volume);
858:     }
859:   } else {

861:     if (nverts == 8) { /* Linear Hexahedra */
862:       const PetscReal &r = quad[0];
863:       const PetscReal &s = quad[1];
864:       const PetscReal &t = quad[2];
865:       const PetscReal dNi_dxi[8]  = {- ( 1.0 - s) * ( 1.0 - t),
866:                                        ( 1.0 - s) * ( 1.0 - t),
867:                                        (       s) * ( 1.0 - t),
868:                                      - (       s) * ( 1.0 - t),
869:                                      - ( 1.0 - s) * (       t),
870:                                        ( 1.0 - s) * (       t),
871:                                        (       s) * (       t),
872:                                      - (       s) * (       t)
873:                                     };

875:       const PetscReal dNi_deta[8]  = { - ( 1.0 - r) * ( 1.0 - t),
876:                                        - (       r) * ( 1.0 - t),
877:                                          (       r) * ( 1.0 - t),
878:                                          ( 1.0 - r) * ( 1.0 - t),
879:                                        - ( 1.0 - r) * (       t),
880:                                        - (       r) * (       t),
881:                                          (       r) * (       t),
882:                                          ( 1.0 - r) * (       t)
883:                                       };

885:       const PetscReal dNi_dzeta[8]  = { - ( 1.0 - r) * ( 1.0 - s),
886:                                         - (       r) * ( 1.0 - s),
887:                                         - (       r) * (       s),
888:                                         - ( 1.0 - r) * (       s),
889:                                           ( 1.0 - r) * ( 1.0 - s),
890:                                           (       r) * ( 1.0 - s),
891:                                           (       r) * (       s),
892:                                           ( 1.0 - r) * (       s)
893:                                      };

895:       for (i = 0; i < nverts; ++i) {
896:         const PetscReal* vertex = coordinates + i * 3;
897:         jacobian[0] += dNi_dxi[i]   * vertex[0];
898:         jacobian[3] += dNi_dxi[i]   * vertex[1];
899:         jacobian[6] += dNi_dxi[i]   * vertex[2];
900:         jacobian[1] += dNi_deta[i]  * vertex[0];
901:         jacobian[4] += dNi_deta[i]  * vertex[1];
902:         jacobian[7] += dNi_deta[i]  * vertex[2];
903:         jacobian[2] += dNi_dzeta[i] * vertex[0];
904:         jacobian[5] += dNi_dzeta[i] * vertex[1];
905:         jacobian[8] += dNi_dzeta[i] * vertex[2];
906:       }
907:     } else if (nverts == 4) { /* Linear Tetrahedra */
908:       const PetscReal x0 = coordinates[/*0 * 3 +*/ 0], y0 = coordinates[/*0 * 3 +*/ 1], z0 = coordinates[/*0 * 3 +*/ 2];

910:       /* compute the jacobian */
911:       jacobian[0] = coordinates[1 * 3 + 0] - x0;  jacobian[1] = coordinates[2 * 3 + 0] - x0; jacobian[2] = coordinates[3 * 3 + 0] - x0;
912:       jacobian[3] = coordinates[1 * 3 + 1] - y0;  jacobian[4] = coordinates[2 * 3 + 1] - y0; jacobian[5] = coordinates[3 * 3 + 1] - y0;
913:       jacobian[6] = coordinates[1 * 3 + 2] - z0;  jacobian[7] = coordinates[2 * 3 + 2] - z0; jacobian[8] = coordinates[3 * 3 + 2] - z0;
914:     } else SETERRQ1(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 : %D", nverts);

916:     if (ijacobian) {
917:       /* invert the jacobian */
918:       DMatrix_Invert_3x3_Internal(jacobian, ijacobian, &volume);
919:     }

921:   }
922:   if (volume < 1e-12) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Element has zero volume: %g. Degenerate element or invalid connectivity\n", volume);
923:   if (dvolume) *dvolume = volume;
924:   return(0);
925: }

927: PetscErrorCode FEMComputeBasis_JandF(const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal *quadrature, PetscReal *phypts,
928:                                      PetscReal *phibasis, PetscReal *jacobian, PetscReal *ijacobian, PetscReal* volume)
929: {

933:   switch (dim) {
934:     case 1:
935:       Compute_Lagrange_Basis_1D_Internal(nverts, coordinates, 1, quadrature, phypts,
936:             NULL, phibasis, NULL, jacobian, ijacobian, volume);
937:       break;
938:     case 2:
939:       Compute_Lagrange_Basis_2D_Internal(nverts, coordinates, 1, quadrature, phypts,
940:             NULL, phibasis, NULL, NULL, jacobian, ijacobian, volume);
941:       break;
942:     case 3:
943:       Compute_Lagrange_Basis_3D_Internal(nverts, coordinates, 1, quadrature, phypts,
944:             NULL, phibasis, NULL, NULL, NULL, jacobian, ijacobian, volume);
945:       break;
946:     default:
947:       SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %D", dim);
948:   }
949:   return(0);
950: }

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

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

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

967:   Level: advanced

969: @*/
970: PetscErrorCode DMMoabPToRMapping(const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal* xphy, PetscReal* natparam, PetscReal* phi)
971: {
972:   /* Perform inverse evaluation for the mapping with use of Newton Raphson iteration */
973:   const PetscReal tol = 1.0e-10;
974:   const PetscInt  max_iterations = 10;
975:   const PetscReal error_tol_sqr = tol*tol;
976:   PetscReal       phibasis[8], jacobian[9], ijacobian[9], volume;
977:   PetscReal       phypts[3] = {0.0, 0.0, 0.0};
978:   PetscReal       delta[3] = {0.0, 0.0, 0.0};
979:   PetscErrorCode  ierr;
980:   PetscInt        iters=0;
981:   PetscReal       error=1.0;


988:   PetscArrayzero(jacobian, dim * dim);
989:   PetscArrayzero(ijacobian, dim * dim);
990:   PetscArrayzero(phibasis, nverts);

992:   /* zero initial guess */
993:   natparam[0] = natparam[1] = natparam[2] = 0.0;

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

998:   error = 0.0;
999:   switch(dim) {
1000:     case 3:
1001:       delta[2] = phypts[2] - xphy[2];
1002:       error += (delta[2]*delta[2]);
1003:     case 2:
1004:       delta[1] = phypts[1] - xphy[1];
1005:       error += (delta[1]*delta[1]);
1006:     case 1:
1007:       delta[0] = phypts[0] - xphy[0];
1008:       error += (delta[0]*delta[0]);
1009:       break;
1010:   }

1012:   while (error > error_tol_sqr) {

1014:     if (++iters > max_iterations)
1015:       SETERRQ3(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]);

1017:     /* Compute natparam -= J.inverse() * delta */
1018:     switch(dim) {
1019:       case 1:
1020:         natparam[0] -= ijacobian[0] * delta[0];
1021:         break;
1022:       case 2:
1023:         natparam[0] -= ijacobian[0] * delta[0] + ijacobian[1] * delta[1];
1024:         natparam[1] -= ijacobian[2] * delta[0] + ijacobian[3] * delta[1];
1025:         break;
1026:       case 3:
1027:         natparam[0] -= ijacobian[0] * delta[0] + ijacobian[1] * delta[1] + ijacobian[2] * delta[2];
1028:         natparam[1] -= ijacobian[3] * delta[0] + ijacobian[4] * delta[1] + ijacobian[5] * delta[2];
1029:         natparam[2] -= ijacobian[6] * delta[0] + ijacobian[7] * delta[1] + ijacobian[8] * delta[2];
1030:         break;
1031:     }

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

1036:     error = 0.0;
1037:     switch(dim) {
1038:       case 3:
1039:         delta[2] = phypts[2] - xphy[2];
1040:         error += (delta[2]*delta[2]);
1041:       case 2:
1042:         delta[1] = phypts[1] - xphy[1];
1043:         error += (delta[1]*delta[1]);
1044:       case 1:
1045:         delta[0] = phypts[0] - xphy[0];
1046:         error += (delta[0]*delta[0]);
1047:         break;
1048:     }
1049:   }
1050:   if (phi) {
1051:     PetscArraycpy(phi, phibasis, nverts);
1052:   }
1053:   return(0);
1054: }