Actual source code: dmmbfem.cxx

petsc-3.10.5 2019-03-28
Report Typos and Errors

  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:   if (!inmat) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_POINTER, "Invalid input matrix specified for 2x2 inversion.");
 15:   PetscReal det = DMatrix_Determinant_2x2_Internal(inmat);
 16:   if (outmat) {
 17:     outmat[0] = inmat[3] / det;
 18:     outmat[1] = -inmat[1] / det;
 19:     outmat[2] = -inmat[2] / det;
 20:     outmat[3] = inmat[0] / det;
 21:   }
 22:   if (determinant) *determinant = det;
 23:   return(0);
 24: }

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

 33: static inline PetscErrorCode DMatrix_Invert_3x3_Internal (const PetscReal *inmat, PetscReal *outmat, PetscScalar *determinant)
 34: {
 35:   if (!inmat) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_POINTER, "Invalid input matrix specified for 3x3 inversion.");
 36:   PetscReal det = DMatrix_Determinant_3x3_Internal(inmat);
 37:   if (outmat) {
 38:     outmat[0] = (inmat[8] * inmat[4] - inmat[7] * inmat[5]) / det;
 39:     outmat[1] = -(inmat[8] * inmat[1] - inmat[7] * inmat[2]) / det;
 40:     outmat[2] = (inmat[5] * inmat[1] - inmat[4] * inmat[2]) / det;
 41:     outmat[3] = -(inmat[8] * inmat[3] - inmat[6] * inmat[5]) / det;
 42:     outmat[4] = (inmat[8] * inmat[0] - inmat[6] * inmat[2]) / det;
 43:     outmat[5] = -(inmat[5] * inmat[0] - inmat[3] * inmat[2]) / det;
 44:     outmat[6] = (inmat[7] * inmat[3] - inmat[6] * inmat[4]) / det;
 45:     outmat[7] = -(inmat[7] * inmat[0] - inmat[6] * inmat[1]) / det;
 46:     outmat[8] = (inmat[4] * inmat[0] - inmat[3] * inmat[1]) / det;
 47:   }
 48:   if (determinant) *determinant = det;
 49:   return(0);
 50: }

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

 73: inline PetscErrorCode DMatrix_Invert_4x4_Internal (PetscReal *inmat, PetscReal *outmat, PetscScalar *determinant)
 74: {
 75:   if (!inmat) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_POINTER, "Invalid input matrix specified for 4x4 inversion.");
 76:   PetscReal det = DMatrix_Determinant_4x4_Internal(inmat);
 77:   if (outmat) {
 78:     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;
 79:     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;
 80:     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;
 81:     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;
 82:     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;
 83:     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;
 84:     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;
 85:     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;
 86:     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;
 87:     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;
 88:     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;
 89:     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;
 90:     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;
 91:     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;
 92:     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;
 93:     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;
 94:   }
 95:   if (determinant) *determinant = det;
 96:   return(0);
 97: }


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

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

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

108:   Notes:

110:   Example Physical Element

112:     1-------2        1----3----2
113:       EDGE2             EDGE3

115:   Input Parameter:

117: .  PetscInt  nverts,           the number of element vertices
118: .  PetscReal coords[3*nverts], the physical coordinates of the vertices (in canonical numbering)
119: .  PetscInt  npts,             the number of evaluation points (quadrature points)
120: .  PetscReal quad[3*npts],     the evaluation points (quadrature points) in the reference space

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

128:   Level: advanced

130: .keywords: DMMoab, FEM, 1-D
131: @*/
132: PetscErrorCode Compute_Lagrange_Basis_1D_Internal ( const PetscInt nverts, const PetscReal *coords, const PetscInt npts, const PetscReal *quad, PetscReal *phypts,
133:     PetscReal *jxw, PetscReal *phi, PetscReal *dphidx,
134:     PetscReal *jacobian, PetscReal *ijacobian, PetscReal *volume)
135: {
136:   int i, j;
137:   PetscErrorCode  ierr;

143:   if (phypts) {
144:     PetscMemzero(phypts, npts * 3 * sizeof(PetscReal));
145:   }
146:   if (dphidx) { /* Reset arrays. */
147:     PetscMemzero(dphidx, npts * nverts * sizeof(PetscReal));
148:   }
149:   if (nverts == 2) { /* Linear Edge */

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

156:       phi[0 + offset] = ( 1.0 - r );
157:       phi[1 + offset] = (       r );

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

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

170:       /* invert the jacobian */
171:       *volume = jacobian[0];
172:       ijacobian[0] = 1.0 / jacobian[0];
173:       jxw[j] *= *volume;

175:       /*  Divide by element jacobian. */
176:       for ( i = 0; i < nverts; i++ ) {
177:         if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[0];
178:       }

180:     }
181:   }
182:   else if (nverts == 3) { /* Quadratic Edge */

184:     for (j = 0; j < npts; j++)
185:     {
186:       const PetscInt offset = j * nverts;
187:       const PetscReal r = quad[j];

189:       phi[0 + offset] = 1.0 + r * ( 2.0 * r - 3.0 );
190:       phi[1 + offset] = 4.0 * r * ( 1.0 - r );
191:       phi[2 + offset] = r * ( 2.0 * r - 1.0 );

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

195:       jacobian[0] = ijacobian[0] = volume[0] = 0.0;
196:       for (i = 0; i < nverts; ++i) {
197:         const PetscReal* vertices = coords + i * 3;
198:         jacobian[0] += dNi_dxi[i] * vertices[0];
199:         if (phypts) {
200:           phypts[3 * j + 0] += phi[i + offset] * vertices[0];
201:         }
202:       }

204:       /* invert the jacobian */
205:       *volume = jacobian[0];
206:       ijacobian[0] = 1.0 / jacobian[0];
207:       if (jxw) jxw[j] *= *volume;

209:       /*  Divide by element jacobian. */
210:       for ( i = 0; i < nverts; i++ ) {
211:         if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[0];
212:       }
213:     }
214:   }
215:   else {
216:     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);
217:   }
218: #if 0
219:   /* verify if the computed basis functions are consistent */
220:   for ( j = 0; j < npts; j++ ) {
221:     PetscScalar phisum = 0, dphixsum = 0;
222:     for ( i = 0; i < nverts; i++ ) {
223:       phisum += phi[i + j * nverts];
224:       if (dphidx) dphixsum += dphidx[i + j * nverts];
225:     }
226:     PetscPrintf(PETSC_COMM_WORLD, "Sum of basis at quadrature point %D = %g, %g, %g\n", j, phisum, dphixsum);
227:   }
228: #endif
229:   return(0);
230: }


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

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

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

241:   Notes:

243:   Example Physical Element (QUAD4)

245:     4------3        s
246:     |      |        |
247:     |      |        |
248:     |      |        |
249:     1------2        0-------r

251:   Input Parameter:

253: .  PetscInt  nverts,           the number of element vertices
254: .  PetscReal coords[3*nverts], the physical coordinates of the vertices (in canonical numbering)
255: .  PetscInt  npts,             the number of evaluation points (quadrature points)
256: .  PetscReal quad[3*npts],     the evaluation points (quadrature points) in the reference space

258:   Output Parameter:
259: .  PetscReal phypts[3*npts],   the evaluation points (quadrature points) transformed to the physical space
260: .  PetscReal jxw[npts],        the jacobian determinant * quadrature weight necessary for assembling discrete contributions
261: .  PetscReal phi[npts],        the bases evaluated at the specified quadrature points
262: .  PetscReal dphidx[npts],     the derivative of the bases wrt X-direction evaluated at the specified quadrature points
263: .  PetscReal dphidy[npts],     the derivative of the bases wrt Y-direction evaluated at the specified quadrature points

265:   Level: advanced

267: .keywords: DMMoab, FEM, 2-D
268: @*/
269: PetscErrorCode Compute_Lagrange_Basis_2D_Internal ( const PetscInt nverts, const PetscReal *coords, const PetscInt npts, const PetscReal *quad, PetscReal *phypts,
270:     PetscReal *jxw, PetscReal *phi, PetscReal *dphidx, PetscReal *dphidy,
271:     PetscReal *jacobian, PetscReal *ijacobian, PetscReal *volume)
272: {
273:   PetscInt i, j, k;
274:   PetscErrorCode   ierr;

280:   PetscMemzero(phi, npts * sizeof(PetscReal));
281:   if (phypts) {
282:     PetscMemzero(phypts, npts * 3 * sizeof(PetscReal));
283:   }
284:   if (dphidx) { /* Reset arrays. */
285:     PetscMemzero(dphidx, npts * nverts * sizeof(PetscReal));
286:     PetscMemzero(dphidy, npts * nverts * sizeof(PetscReal));
287:   }
288:   if (nverts == 4) { /* Linear Quadrangle */

290:     for (j = 0; j < npts; j++)
291:     {
292:       const PetscInt offset = j * nverts;
293:       const PetscReal r = quad[0 + j * 2];
294:       const PetscReal s = quad[1 + j * 2];

296:       phi[0 + offset] = ( 1.0 - r ) * ( 1.0 - s );
297:       phi[1 + offset] =         r   * ( 1.0 - s );
298:       phi[2 + offset] =         r   *         s;
299:       phi[3 + offset] = ( 1.0 - r ) *         s;

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

304:       PetscMemzero(jacobian, 4 * sizeof(PetscReal));
305:       PetscMemzero(ijacobian, 4 * sizeof(PetscReal));
306:       for (i = 0; i < nverts; ++i) {
307:         const PetscReal* vertices = coords + i * 3;
308:         jacobian[0] += dNi_dxi[i] * vertices[0];
309:         jacobian[2] += dNi_dxi[i] * vertices[1];
310:         jacobian[1] += dNi_deta[i] * vertices[0];
311:         jacobian[3] += dNi_deta[i] * vertices[1];
312:         if (phypts) {
313:           phypts[3 * j + 0] += phi[i + offset] * vertices[0];
314:           phypts[3 * j + 1] += phi[i + offset] * vertices[1];
315:           phypts[3 * j + 2] += phi[i + offset] * vertices[2];
316:         }
317:       }

319:       /* invert the jacobian */
320:       DMatrix_Invert_2x2_Internal(jacobian, ijacobian, volume);
321:       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);

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

325:       /*  Let us compute the bases derivatives by scaling with inverse jacobian. */
326:       if (dphidx) {
327:         for ( i = 0; i < nverts; i++ ) {
328:           for ( k = 0; k < 2; ++k) {
329:             if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[k * 2 + 0];
330:             if (dphidy) dphidy[i + offset] += dNi_deta[i] * ijacobian[k * 2 + 1];
331:           } /* for k=[0..2) */
332:         } /* for i=[0..nverts) */
333:       } /* if (dphidx) */
334:     }
335:   }
336:   else if (nverts == 3) { /* Linear triangle */

338:     PetscMemzero(jacobian, 4 * sizeof(PetscReal));
339:     PetscMemzero(ijacobian, 4 * sizeof(PetscReal));

341:     const PetscReal x2 = coords[2 * 3 + 0], y2 = coords[2 * 3 + 1];

343:     /* Jacobian is constant */
344:     jacobian[0] = (coords[0 * 3 + 0] - x2); jacobian[1] = (coords[1 * 3 + 0] - x2);
345:     jacobian[2] = (coords[0 * 3 + 1] - y2); jacobian[3] = (coords[1 * 3 + 1] - y2);

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

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

354:     for (j = 0; j < npts; j++)
355:     {
356:       const PetscInt offset = j * nverts;
357:       const PetscReal r = quad[0 + j * 2];
358:       const PetscReal s = quad[1 + j * 2];

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

362:       /* const PetscReal Ni[3]  = { r, s, 1.0 - r - s }; */
363:       const PetscReal phipts_x = coords[2 * 3 + 0] + jacobian[0] * r + jacobian[1] * s;
364:       const PetscReal phipts_y = coords[2 * 3 + 1] + jacobian[2] * r + jacobian[3] * s;
365:       if (phypts) {
366:         phypts[offset + 0] = phipts_x;
367:         phypts[offset + 1] = phipts_y;
368:       }

370:       /* \phi_0 = (b.y � c.y) x + (b.x � c.x) y + c.x b.y � b.x c.y */
371:       phi[0 + offset] = (  ijacobian[0] * (phipts_x - x2) + ijacobian[1] * (phipts_y - y2) );
372:       /* \phi_1 = (c.y � a.y) x + (a.x � c.x) y + c.x a.y � a.x c.y */
373:       phi[1 + offset] = (  ijacobian[2] * (phipts_x - x2) + ijacobian[3] * (phipts_y - y2) );
374:       phi[2 + offset] = 1.0 - phi[0 + offset] - phi[1 + offset];

376:       if (dphidx) {
377:         dphidx[0 + offset] = Dx[0];
378:         dphidx[1 + offset] = Dx[1];
379:         dphidx[2 + offset] = Dx[2];
380:       }

382:       if (dphidy) {
383:         dphidy[0 + offset] = Dy[0];
384:         dphidy[1 + offset] = Dy[1];
385:         dphidy[2 + offset] = Dy[2];
386:       }

388:     }
389:   }
390:   else {
391:     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);
392:   }
393: #if 0
394:   /* verify if the computed basis functions are consistent */
395:   for ( j = 0; j < npts; j++ ) {
396:     PetscScalar phisum = 0, dphixsum = 0, dphiysum = 0;
397:     for ( i = 0; i < nverts; i++ ) {
398:       phisum += phi[i + j * nverts];
399:       if (dphidx) dphixsum += dphidx[i + j * nverts];
400:       if (dphidy) dphiysum += dphidy[i + j * nverts];
401:     }
402:     PetscPrintf(PETSC_COMM_WORLD, "Sum of basis at quadrature point %D = %g, %g, %g\n", j, phisum, dphixsum, dphiysum);
403:   }
404: #endif
405:   return(0);
406: }


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

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

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

417:   Notes:

419:   Example Physical Element (HEX8)

421:       8------7
422:      /|     /|        t  s
423:     5------6 |        | /
424:     | |    | |        |/
425:     | 4----|-3        0-------r
426:     |/     |/
427:     1------2

429:   Input Parameter:

431: .  PetscInt  nverts,           the number of element vertices
432: .  PetscReal coords[3*nverts], the physical coordinates of the vertices (in canonical numbering)
433: .  PetscInt  npts,             the number of evaluation points (quadrature points)
434: .  PetscReal quad[3*npts],     the evaluation points (quadrature points) in the reference space

436:   Output Parameter:
437: .  PetscReal phypts[3*npts],   the evaluation points (quadrature points) transformed to the physical space
438: .  PetscReal jxw[npts],        the jacobian determinant * quadrature weight necessary for assembling discrete contributions
439: .  PetscReal phi[npts],        the bases evaluated at the specified quadrature points
440: .  PetscReal dphidx[npts],     the derivative of the bases wrt X-direction evaluated at the specified quadrature points
441: .  PetscReal dphidy[npts],     the derivative of the bases wrt Y-direction evaluated at the specified quadrature points
442: .  PetscReal dphidz[npts],     the derivative of the bases wrt Z-direction evaluated at the specified quadrature points

444:   Level: advanced

446: .keywords: DMMoab, FEM, 3-D
447: @*/
448: PetscErrorCode Compute_Lagrange_Basis_3D_Internal ( const PetscInt nverts, const PetscReal *coords, const PetscInt npts, const PetscReal *quad, PetscReal *phypts,
449:     PetscReal *jxw, PetscReal *phi, PetscReal *dphidx, PetscReal *dphidy, PetscReal *dphidz,
450:     PetscReal *jacobian, PetscReal *ijacobian, PetscReal *volume)
451: {
452:   PetscInt i, j, k;

459:   /* Reset arrays. */
460:   PetscMemzero(phi, npts * sizeof(PetscReal));
461:   if (phypts) {
462:     PetscMemzero(phypts, npts * 3 * sizeof(PetscReal));
463:   }
464:   if (dphidx) {
465:     PetscMemzero(dphidx, npts * nverts * sizeof(PetscReal));
466:     PetscMemzero(dphidy, npts * nverts * sizeof(PetscReal));
467:     PetscMemzero(dphidz, npts * nverts * sizeof(PetscReal));
468:   }

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

472:     for (j = 0; j < npts; j++)
473:     {
474:       const PetscInt offset = j * nverts;
475:       const PetscReal& r = quad[j * 3 + 0];
476:       const PetscReal& s = quad[j * 3 + 1];
477:       const PetscReal& t = quad[j * 3 + 2];

479:       phi[offset + 0] = ( 1.0 - r ) * ( 1.0 - s ) * ( 1.0 - t ); /* 0,0,0 */
480:       phi[offset + 1] = (       r ) * ( 1.0 - s ) * ( 1.0 - t ); /* 1,0,0 */
481:       phi[offset + 2] = (       r ) * (       s ) * ( 1.0 - t ); /* 1,1,0 */
482:       phi[offset + 3] = ( 1.0 - r ) * (       s ) * ( 1.0 - t ); /* 0,1,0 */
483:       phi[offset + 4] = ( 1.0 - r ) * ( 1.0 - s ) * (       t ); /* 0,0,1 */
484:       phi[offset + 5] = (       r ) * ( 1.0 - s ) * (       t ); /* 1,0,1 */
485:       phi[offset + 6] = (       r ) * (       s ) * (       t ); /* 1,1,1 */
486:       phi[offset + 7] = ( 1.0 - r ) * (       s ) * (       t ); /* 0,1,1 */

488:       const PetscReal dNi_dxi[8]  = {- ( 1.0 - s ) * ( 1.0 - t ),
489:                                        ( 1.0 - s ) * ( 1.0 - t ),
490:                                        (       s ) * ( 1.0 - t ),
491:                                      - (       s ) * ( 1.0 - t ),
492:                                      - ( 1.0 - s ) * (       t ),
493:                                        ( 1.0 - s ) * (       t ),
494:                                        (       s ) * (       t ),
495:                                      - (       s ) * (       t )
496:                                     };

498:       const PetscReal dNi_deta[8]  = { - ( 1.0 - r ) * ( 1.0 - t ),
499:                                        - (       r ) * ( 1.0 - t ),
500:                                          (       r ) * ( 1.0 - t ),
501:                                          ( 1.0 - r ) * ( 1.0 - t ),
502:                                        - ( 1.0 - r ) * (       t ),
503:                                        - (       r ) * (       t ),
504:                                          (       r ) * (       t ),
505:                                          ( 1.0 - r ) * (       t )
506:                                       };

508:       const PetscReal dNi_dzeta[8]  = { - ( 1.0 - r ) * ( 1.0 - s ),
509:                                         - (       r ) * ( 1.0 - s ),
510:                                         - (       r ) * (       s ),
511:                                         - ( 1.0 - r ) * (       s ),
512:                                           ( 1.0 - r ) * ( 1.0 - s ),
513:                                           (       r ) * ( 1.0 - s ),
514:                                           (       r ) * (       s ),
515:                                           ( 1.0 - r ) * (       s )
516:                                      };

518:       PetscMemzero(jacobian, 9 * sizeof(PetscReal));
519:       PetscMemzero(ijacobian, 9 * sizeof(PetscReal));
520:       for (i = 0; i < nverts; ++i) {
521:         const PetscReal* vertex = coords + i * 3;
522:         jacobian[0] += dNi_dxi[i]   * vertex[0];
523:         jacobian[3] += dNi_dxi[i]   * vertex[1];
524:         jacobian[6] += dNi_dxi[i]   * vertex[2];
525:         jacobian[1] += dNi_deta[i]  * vertex[0];
526:         jacobian[4] += dNi_deta[i]  * vertex[1];
527:         jacobian[7] += dNi_deta[i]  * vertex[2];
528:         jacobian[2] += dNi_dzeta[i] * vertex[0];
529:         jacobian[5] += dNi_dzeta[i] * vertex[1];
530:         jacobian[8] += dNi_dzeta[i] * vertex[2];
531:         if (phypts) {
532:           phypts[3 * j + 0] += phi[i + offset] * vertex[0];
533:           phypts[3 * j + 1] += phi[i + offset] * vertex[1];
534:           phypts[3 * j + 2] += phi[i + offset] * vertex[2];
535:         }
536:       }

538:       /* invert the jacobian */
539:       DMatrix_Invert_3x3_Internal(jacobian, ijacobian, volume);
540:       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);

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

544:       /*  Divide by element jacobian. */
545:       for ( i = 0; i < nverts; ++i ) {
546:         for ( k = 0; k < 3; ++k) {
547:           if (dphidx) dphidx[i + offset] += dNi_dxi[i]   * ijacobian[0 * 3 + k];
548:           if (dphidy) dphidy[i + offset] += dNi_deta[i]  * ijacobian[1 * 3 + k];
549:           if (dphidz) dphidz[i + offset] += dNi_dzeta[i] * ijacobian[2 * 3 + k];
550:         }
551:       }
552:     }
553:   }
554:   else if (nverts == 4) { /* Linear Tetrahedra */

556:     PetscMemzero(jacobian, 9 * sizeof(PetscReal));
557:     PetscMemzero(ijacobian, 9 * sizeof(PetscReal));

559:     const PetscReal x0 = coords[/*0 * 3 +*/ 0], y0 = coords[/*0 * 3 +*/ 1], z0 = coords[/*0 * 3 +*/ 2];

561:     /* compute the jacobian */
562:     jacobian[0] = coords[1 * 3 + 0] - x0;  jacobian[1] = coords[2 * 3 + 0] - x0; jacobian[2] = coords[3 * 3 + 0] - x0;
563:     jacobian[3] = coords[1 * 3 + 1] - y0;  jacobian[4] = coords[2 * 3 + 1] - y0; jacobian[5] = coords[3 * 3 + 1] - y0;
564:     jacobian[6] = coords[1 * 3 + 2] - z0;  jacobian[7] = coords[2 * 3 + 2] - z0; jacobian[8] = coords[3 * 3 + 2] - z0;

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

570:     PetscReal Dx[4]={0,0,0,0}, Dy[4]={0,0,0,0}, Dz[4]={0,0,0,0};
571:     if (dphidx) {
572:       Dx[0] =   ( coords[1 + 2 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3] )
573:                  - coords[1 + 1 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3] )
574:                  - coords[1 + 3 * 3] * ( coords[2 + 1 * 3] - coords[2 + 2 * 3] )
575:                 ) / *volume;
576:       Dx[1] = - ( coords[1 + 2 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3] )
577:                  - coords[1 + 0 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3] )
578:                  - coords[1 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 2 * 3] )
579:                 ) / *volume;
580:       Dx[2] =   ( coords[1 + 1 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3] )
581:                  - coords[1 + 0 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3] )
582:                  - coords[1 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 1 * 3] )
583:                 ) / *volume;
584:       Dx[3] =  - ( Dx[0] + Dx[1] + Dx[2] );
585:       Dy[0] =   ( coords[0 + 1 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3] )
586:                  - coords[0 + 2 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3] )
587:                  + coords[0 + 3 * 3] * ( coords[2 + 1 * 3] - coords[2 + 2 * 3] )
588:                 ) / *volume;
589:       Dy[1] = - ( coords[0 + 0 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3] )
590:                  - coords[0 + 2 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3] )
591:                  + coords[0 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 2 * 3] )
592:                 ) / *volume;
593:       Dy[2] =   ( coords[0 + 0 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3] )
594:                  - coords[0 + 1 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3] )
595:                  + coords[0 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 1 * 3] )
596:                 ) / *volume;
597:       Dy[3] =  - ( Dy[0] + Dy[1] + Dy[2] );
598:       Dz[0] =   ( coords[0 + 1 * 3] * (coords[1 + 3 * 3] - coords[1 + 2 * 3] )
599:                  - coords[0 + 2 * 3] * (coords[1 + 3 * 3] - coords[1 + 1 * 3] )
600:                  + coords[0 + 3 * 3] * (coords[1 + 2 * 3] - coords[1 + 1 * 3] )
601:                 ) / *volume;
602:       Dz[1] = - ( coords[0 + 0 * 3] * (coords[1 + 3 * 3] - coords[1 + 2 * 3] )
603:                   + coords[0 + 2 * 3] * (coords[1 + 0 * 3] - coords[1 + 3 * 3] )
604:                   - coords[0 + 3 * 3] * (coords[1 + 0 * 3] - coords[1 + 2 * 3] )
605:                 ) / *volume;
606:       Dz[2] =   ( coords[0 + 0 * 3] * (coords[1 + 3 * 3] - coords[1 + 1 * 3] )
607:                  + coords[0 + 1 * 3] * (coords[1 + 0 * 3] - coords[1 + 3 * 3] )
608:                  - coords[0 + 3 * 3] * (coords[1 + 0 * 3] - coords[1 + 1 * 3] )
609:                 ) / *volume;
610:       Dz[3] =  - ( Dz[0] + Dz[1] + Dz[2] );
611:     }

613:     for ( j = 0; j < npts; j++ )
614:     {
615:       const PetscInt offset = j * nverts;
616:       const PetscReal& r = quad[j * 3 + 0];
617:       const PetscReal& s = quad[j * 3 + 1];
618:       const PetscReal& t = quad[j * 3 + 2];

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

622:       phi[offset + 0] = 1.0 - r - s - t;
623:       phi[offset + 1] = r;
624:       phi[offset + 2] = s;
625:       phi[offset + 3] = t;

627:       if (phypts) {
628:         for (i = 0; i < nverts; ++i) {
629:           const PetscScalar* vertices = coords + i * 3;
630:           phypts[3 * j + 0] += phi[i + offset] * vertices[0];
631:           phypts[3 * j + 1] += phi[i + offset] * vertices[1];
632:           phypts[3 * j + 2] += phi[i + offset] * vertices[2];
633:         }
634:       }

636:       /* Now set the derivatives */
637:       if (dphidx) {
638:         dphidx[0 + offset] = Dx[0];
639:         dphidx[1 + offset] = Dx[1];
640:         dphidx[2 + offset] = Dx[2];
641:         dphidx[3 + offset] = Dx[3];
642: 
643:         dphidy[0 + offset] = Dy[0];
644:         dphidy[1 + offset] = Dy[1];
645:         dphidy[2 + offset] = Dy[2];
646:         dphidy[3 + offset] = Dy[3];
647: 
648:         dphidz[0 + offset] = Dz[0];
649:         dphidz[1 + offset] = Dz[1];
650:         dphidz[2 + offset] = Dz[2];
651:         dphidz[3 + offset] = Dz[3];
652:       }

654:     } /* Tetrahedra -- ends */
655:   }
656:   else
657:   {
658:     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);
659:   }
660: #if 0
661:   /* verify if the computed basis functions are consistent */
662:   for ( j = 0; j < npts; j++ ) {
663:     PetscScalar phisum = 0, dphixsum = 0, dphiysum = 0, dphizsum = 0;
664:     const int offset = j * nverts;
665:     for ( i = 0; i < nverts; i++ ) {
666:       phisum += phi[i + offset];
667:       if (dphidx) dphixsum += dphidx[i + offset];
668:       if (dphidy) dphiysum += dphidy[i + offset];
669:       if (dphidz) dphizsum += dphidz[i + offset];
670:       if (dphidx) PetscPrintf(PETSC_COMM_WORLD, "\t Values [%d]: [JxW] [phi, dphidx, dphidy, dphidz] = %g, %g, %g, %g, %g\n", j, jxw[j], phi[i + offset], dphidx[i + offset], dphidy[i + offset], dphidz[i + offset]);
671:     }
672:     if (dphidx) PetscPrintf(PETSC_COMM_WORLD, "Sum of basis at quadrature point %D (%g, %g, %g) = %g, %g, %g, %g\n", j, quad[3 * j + 0], quad[3 * j + 1], quad[3 * j + 2], phisum, dphixsum, dphiysum, dphizsum);
673:   }
674: #endif
675:   return(0);
676: }



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

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

686:   Input Parameter:

688: .  PetscInt  nverts,           the number of element vertices
689: .  PetscReal coords[3*nverts], the physical coordinates of the vertices (in canonical numbering)
690: .  PetscInt  npts,             the number of evaluation points (quadrature points)
691: .  PetscReal quad[3*npts],     the evaluation points (quadrature points) in the reference space

693:   Output Parameter:
694: .  PetscReal phypts[3*npts],   the evaluation points (quadrature points) transformed to the physical space
695: .  PetscReal jxw[npts],        the jacobian determinant * quadrature weight necessary for assembling discrete contributions
696: .  PetscReal fe_basis[npts],        the bases values evaluated at the specified quadrature points
697: .  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

699:   Level: advanced

701: .keywords: DMMoab, FEM, 3-D
702: @*/
703: PetscErrorCode DMMoabFEMComputeBasis ( const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscQuadrature quadrature, 
704:                                        PetscReal *phypts, PetscReal *jacobian_quadrature_weight_product, 
705:                                        PetscReal *fe_basis, PetscReal **fe_basis_derivatives)
706: {
707:   PetscErrorCode  ierr;
708:   PetscInt        npoints,idim;
709:   bool            compute_der;
710:   const PetscReal *quadpts, *quadwts;
711:   PetscReal       jacobian[9], ijacobian[9], volume;

717:   compute_der = (fe_basis_derivatives != NULL);

719:   /* Get the quadrature points and weights for the given quadrature rule */
720:   PetscQuadratureGetData(quadrature, &idim, NULL, &npoints, &quadpts, &quadwts);
721:   if (idim != dim) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Dimension mismatch: provided (%D) vs quadrature (%D)\n",idim,dim);
722:   if (jacobian_quadrature_weight_product) {
723:     PetscMemcpy(jacobian_quadrature_weight_product, quadwts, npoints * sizeof(PetscReal));
724:   }

726:   switch (dim) {
727:   case 1:
728:     Compute_Lagrange_Basis_1D_Internal(nverts, coordinates, npoints, quadpts, phypts,
729:            jacobian_quadrature_weight_product, fe_basis,
730:            (compute_der ? fe_basis_derivatives[0] : NULL),
731:            jacobian, ijacobian, &volume);
732:     break;
733:   case 2:
734:     Compute_Lagrange_Basis_2D_Internal(nverts, coordinates, npoints, quadpts, phypts,
735:            jacobian_quadrature_weight_product, fe_basis,
736:            (compute_der ? fe_basis_derivatives[0] : NULL),
737:            (compute_der ? fe_basis_derivatives[1] : NULL),
738:            jacobian, ijacobian, &volume);
739:     break;
740:   case 3:
741:     Compute_Lagrange_Basis_3D_Internal(nverts, coordinates, npoints, quadpts, phypts,
742:            jacobian_quadrature_weight_product, fe_basis,
743:            (compute_der ? fe_basis_derivatives[0] : NULL),
744:            (compute_der ? fe_basis_derivatives[1] : NULL),
745:            (compute_der ? fe_basis_derivatives[2] : NULL),
746:            jacobian, ijacobian, &volume);
747:     break;
748:   default:
749:     SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %D", dim);
750:   }
751:   return(0);
752: }



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

760:   Input Parameter:

762: .  PetscInt  dim,           the element dimension (1=EDGE, 2=QUAD/TRI, 3=HEX/TET)
763: .  PetscInt nverts,      the number of vertices in the physical element

765:   Output Parameter:
766: .  PetscQuadrature quadrature,  the quadrature object with default settings to integrate polynomials defined over the element

768:   Level: advanced

770: .keywords: DMMoab, Quadrature, PetscDT
771: @*/
772: PetscErrorCode DMMoabFEMCreateQuadratureDefault ( const PetscInt dim, const PetscInt nverts, PetscQuadrature *quadrature )
773: {
774:   PetscReal *w, *x;
775:   PetscInt nc=1;
776:   PetscErrorCode  ierr;

779:   /* Create an appropriate quadrature rule to sample basis */
780:   switch (dim)
781:   {
782:   case 1:
783:     /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */
784:     PetscDTGaussJacobiQuadrature(1, nc, nverts, 0, 1.0, quadrature);
785:     break;
786:   case 2:
787:     /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */
788:     if (nverts == 3) {
789:       const PetscInt order = 2;
790:       const PetscInt npoints = (order == 2 ? 3 : 6);
791:       PetscMalloc2(npoints * 2, &x, npoints, &w);
792:       if (npoints == 3) {
793:         x[0] = x[1] = x[2] = x[5] = 1.0 / 6.0;
794:         x[3] = x[4] = 2.0 / 3.0;
795:         w[0] = w[1] = w[2] = 1.0 / 3.0;
796:       }
797:       else if (npoints == 6) {
798:         x[0] = x[1] = x[2] = 0.44594849091597;
799:         x[3] = x[4] = 0.10810301816807;
800:         x[5] = 0.44594849091597;
801:         x[6] = x[7] = x[8] = 0.09157621350977;
802:         x[9] = x[10] = 0.81684757298046;
803:         x[11] = 0.09157621350977;
804:         w[0] = w[1] = w[2] = 0.22338158967801;
805:         w[3] = w[4] = w[5] = 0.10995174365532;
806:       }
807:       else {
808:         SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Triangle quadrature rules for points 3 and 6 supported; npoints : %D", npoints);
809:       }
810:       PetscQuadratureCreate(PETSC_COMM_SELF, quadrature);
811:       PetscQuadratureSetOrder(*quadrature, order);
812:       PetscQuadratureSetData(*quadrature, dim, nc, npoints, x, w);
813:       /* PetscDTGaussJacobiQuadrature(dim, nverts, 0.0, 1.0, quadrature); */
814:     }
815:     else {
816:       PetscDTGaussTensorQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature);
817:     }
818:     break;
819:   case 3:
820:     /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */
821:     if (nverts == 4) {
822:       PetscInt order;
823:       const PetscInt npoints = 4; // Choose between 4 and 10
824:       PetscMalloc2(npoints * 3, &x, npoints, &w);
825:       if (npoints == 4) { /*  KEAST1, K1,  N=4, O=4 */
826:         const PetscReal x_4[12] = { 0.5854101966249685, 0.1381966011250105, 0.1381966011250105,
827:                                     0.1381966011250105, 0.1381966011250105, 0.1381966011250105,
828:                                     0.1381966011250105, 0.1381966011250105, 0.5854101966249685,
829:                                     0.1381966011250105, 0.5854101966249685, 0.1381966011250105
830:                                   };
831:         PetscMemcpy(x, x_4, 12 * sizeof(PetscReal));

833:         w[0] = w[1] = w[2] = w[3] = 1.0 / 24.0;
834:         order = 4;
835:       }
836:       else if (npoints == 10) { /*  KEAST3, K3  N=10, O=10 */
837:         const PetscReal x_10[30] = { 0.5684305841968444, 0.1438564719343852, 0.1438564719343852,
838:                                      0.1438564719343852, 0.1438564719343852, 0.1438564719343852,
839:                                      0.1438564719343852, 0.1438564719343852, 0.5684305841968444,
840:                                      0.1438564719343852, 0.5684305841968444, 0.1438564719343852,
841:                                      0.0000000000000000, 0.5000000000000000, 0.5000000000000000,
842:                                      0.5000000000000000, 0.0000000000000000, 0.5000000000000000,
843:                                      0.5000000000000000, 0.5000000000000000, 0.0000000000000000,
844:                                      0.5000000000000000, 0.0000000000000000, 0.0000000000000000,
845:                                      0.0000000000000000, 0.5000000000000000, 0.0000000000000000,
846:                                      0.0000000000000000, 0.0000000000000000, 0.5000000000000000
847:                                    };
848:         PetscMemcpy(x, x_10, 30 * sizeof(PetscReal));

850:         w[0] = w[1] = w[2] = w[3] = 0.2177650698804054;
851:         w[4] = w[5] = w[6] = w[7] = w[8] = w[9] = 0.0214899534130631;
852:         order = 10;
853:       }
854:       else {
855:         SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Tetrahedral quadrature rules for points 4 and 10 supported; npoints : %D", npoints);
856:       }
857:       PetscQuadratureCreate(PETSC_COMM_SELF, quadrature);
858:       PetscQuadratureSetOrder(*quadrature, order);
859:       PetscQuadratureSetData(*quadrature, dim, nc, npoints, x, w);
860:       /* PetscDTGaussJacobiQuadrature(dim, nverts, 0.0, 1.0, quadrature); */
861:     }
862:     else {
863:       PetscDTGaussTensorQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature);
864:     }
865:     break;
866:   default:
867:     SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %D", dim);
868:   }
869:   return(0);
870: }

872: /* Compute Jacobians */
873: PetscErrorCode ComputeJacobian_Internal ( const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal *quad, PetscReal *phypts,
874:   PetscReal *jacobian, PetscReal *ijacobian, PetscReal* dvolume)
875: {
876:   PetscInt i;
877:   PetscReal volume=1.0;

884:   PetscMemzero(jacobian, dim * dim * sizeof(PetscReal));
885:   if (ijacobian) {
886:     PetscMemzero(ijacobian, dim * dim * sizeof(PetscReal));
887:   }
888:   if (phypts) {
889:     PetscMemzero(phypts, /*npts=1 * */ 3 * sizeof(PetscReal));
890:   }

892:   if (dim == 1) {

894:     const PetscReal& r = quad[0];
895:     if (nverts == 2) { /* Linear Edge */
896:       const PetscReal dNi_dxi[2]  = { -1.0, 1.0 };

898:       for (i = 0; i < nverts; ++i) {
899:         const PetscReal* vertices = coordinates + i * 3;
900:         jacobian[0] += dNi_dxi[i] * vertices[0];
901:       }
902:     }
903:     else if (nverts == 3) { /* Quadratic Edge */
904:       const PetscReal dNi_dxi[3]  = { 4 * r - 3.0, 4 * ( 1.0 - 2.0 * r ), 4.0 * r - 1.0};

906:       for (i = 0; i < nverts; ++i) {
907:         const PetscReal* vertices = coordinates + i * 3;
908:         jacobian[0] += dNi_dxi[i] * vertices[0];
909:       }
910:     }
911:     else {
912:       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);
913:     }

915:     if (ijacobian) {
916:       /* invert the jacobian */
917:       ijacobian[0] = 1.0 / jacobian[0];
918:     }

920:   }
921:   else if (dim == 2) {

923:     if (nverts == 4) { /* Linear Quadrangle */
924:       const PetscReal& r = quad[0];
925:       const PetscReal& s = quad[1];

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

930:       for (i = 0; i < nverts; ++i) {
931:         const PetscReal* vertices = coordinates + i * 3;
932:         jacobian[0] += dNi_dxi[i]  * vertices[0];
933:         jacobian[2] += dNi_dxi[i]  * vertices[1];
934:         jacobian[1] += dNi_deta[i] * vertices[0];
935:         jacobian[3] += dNi_deta[i] * vertices[1];
936:       }
937:     }
938:     else if (nverts == 3) { /* Linear triangle */
939:       const PetscReal x2 = coordinates[2 * 3 + 0], y2 = coordinates[2 * 3 + 1];

941:       /* Jacobian is constant */
942:       jacobian[0] = (coordinates[0 * 3 + 0] - x2); jacobian[1] = (coordinates[1 * 3 + 0] - x2);
943:       jacobian[2] = (coordinates[0 * 3 + 1] - y2); jacobian[3] = (coordinates[1 * 3 + 1] - y2);
944:     }
945:     else {
946:       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);
947:     }

949:     /* invert the jacobian */
950:     if (ijacobian) {
951:       DMatrix_Invert_2x2_Internal(jacobian, ijacobian, &volume);
952:     }

954:   }
955:   else {

957:     if (nverts == 8) { /* Linear Hexahedra */
958:       const PetscReal& r = quad[0];
959:       const PetscReal& s = quad[1];
960:       const PetscReal& t = quad[2];
961:       const PetscReal dNi_dxi[8]  = {- ( 1.0 - s ) * ( 1.0 - t ),
962:                                        ( 1.0 - s ) * ( 1.0 - t ),
963:                                        (       s ) * ( 1.0 - t ),
964:                                      - (       s ) * ( 1.0 - t ),
965:                                      - ( 1.0 - s ) * (       t ),
966:                                        ( 1.0 - s ) * (       t ),
967:                                        (       s ) * (       t ),
968:                                      - (       s ) * (       t )
969:                                     };

971:       const PetscReal dNi_deta[8]  = { - ( 1.0 - r ) * ( 1.0 - t ),
972:                                        - (       r ) * ( 1.0 - t ),
973:                                          (       r ) * ( 1.0 - t ),
974:                                          ( 1.0 - r ) * ( 1.0 - t ),
975:                                        - ( 1.0 - r ) * (       t ),
976:                                        - (       r ) * (       t ),
977:                                          (       r ) * (       t ),
978:                                          ( 1.0 - r ) * (       t )
979:                                       };

981:       const PetscReal dNi_dzeta[8]  = { - ( 1.0 - r ) * ( 1.0 - s ),
982:                                         - (       r ) * ( 1.0 - s ),
983:                                         - (       r ) * (       s ),
984:                                         - ( 1.0 - r ) * (       s ),
985:                                           ( 1.0 - r ) * ( 1.0 - s ),
986:                                           (       r ) * ( 1.0 - s ),
987:                                           (       r ) * (       s ),
988:                                           ( 1.0 - r ) * (       s )
989:                                      };

991:       for (i = 0; i < nverts; ++i) {
992:         const PetscReal* vertex = coordinates + i * 3;
993:         jacobian[0] += dNi_dxi[i]   * vertex[0];
994:         jacobian[3] += dNi_dxi[i]   * vertex[1];
995:         jacobian[6] += dNi_dxi[i]   * vertex[2];
996:         jacobian[1] += dNi_deta[i]  * vertex[0];
997:         jacobian[4] += dNi_deta[i]  * vertex[1];
998:         jacobian[7] += dNi_deta[i]  * vertex[2];
999:         jacobian[2] += dNi_dzeta[i] * vertex[0];
1000:         jacobian[5] += dNi_dzeta[i] * vertex[1];
1001:         jacobian[8] += dNi_dzeta[i] * vertex[2];
1002:       }

1004:     }
1005:     else if (nverts == 4) { /* Linear Tetrahedra */
1006:       const PetscReal x0 = coordinates[/*0 * 3 +*/ 0], y0 = coordinates[/*0 * 3 +*/ 1], z0 = coordinates[/*0 * 3 +*/ 2];

1008:       /* compute the jacobian */
1009:       jacobian[0] = coordinates[1 * 3 + 0] - x0;  jacobian[1] = coordinates[2 * 3 + 0] - x0; jacobian[2] = coordinates[3 * 3 + 0] - x0;
1010:       jacobian[3] = coordinates[1 * 3 + 1] - y0;  jacobian[4] = coordinates[2 * 3 + 1] - y0; jacobian[5] = coordinates[3 * 3 + 1] - y0;
1011:       jacobian[6] = coordinates[1 * 3 + 2] - z0;  jacobian[7] = coordinates[2 * 3 + 2] - z0; jacobian[8] = coordinates[3 * 3 + 2] - z0;
1012:     } /* Tetrahedra -- ends */
1013:     else
1014:     {
1015:       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);
1016:     }

1018:     if (ijacobian) {
1019:       /* invert the jacobian */
1020:       DMatrix_Invert_3x3_Internal(jacobian, ijacobian, &volume);
1021:     }

1023:   }
1024:   if ( volume < 1e-12 ) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Element has zero volume: %g. Degenerate element or invalid connectivity\n", volume);
1025:   if (dvolume) *dvolume = volume;
1026:   return(0);
1027: }


1030: PetscErrorCode FEMComputeBasis_JandF ( const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal *quadrature, PetscReal *phypts,
1031:                                        PetscReal *phibasis, PetscReal *jacobian, PetscReal *ijacobian, PetscReal* volume  )
1032: {
1033:   PetscErrorCode  ierr;

1036:   switch (dim) {
1037:     case 1:
1038:       Compute_Lagrange_Basis_1D_Internal(nverts, coordinates, 1, quadrature, phypts,
1039:             NULL, phibasis, NULL, jacobian, ijacobian, volume);
1040:       break;
1041:     case 2:
1042:       Compute_Lagrange_Basis_2D_Internal(nverts, coordinates, 1, quadrature, phypts,
1043:             NULL, phibasis, NULL, NULL, jacobian, ijacobian, volume);
1044:       break;
1045:     case 3:
1046:       Compute_Lagrange_Basis_3D_Internal(nverts, coordinates, 1, quadrature, phypts,
1047:             NULL, phibasis, NULL, NULL, NULL, jacobian, ijacobian, volume);
1048:       break;
1049:     default:
1050:       SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %D", dim);
1051:   }
1052:   return(0);
1053: }



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

1062:   Input Parameter:

1064: .  PetscInt  dim,          the element dimension (1=EDGE, 2=QUAD/TRI, 3=HEX/TET)
1065: .  PetscInt nverts,        the number of vertices in the physical element
1066: .  PetscReal coordinates,  the coordinates of vertices in the physical element
1067: .  PetscReal[3] xphy,      the coordinates of physical point for which natural coordinates (in reference frame) are sought

1069:   Output Parameter:
1070: .  PetscReal[3] natparam,  the natural coordinates (in reference frame) corresponding to xphy
1071: .  PetscReal[nverts] phi,  the basis functions evaluated at the natural coordinates (natparam)

1073:   Level: advanced

1075: .keywords: DMMoab, Mapping, FEM
1076: @*/
1077: PetscErrorCode DMMoabPToRMapping( const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal* xphy, PetscReal* natparam, PetscReal* phi)
1078: {
1079:   /* Perform inverse evaluation for the mapping with use of Newton Raphson iteration */
1080:   const PetscReal tol = 1.0e-10;
1081:   const PetscInt max_iterations = 10;
1082:   const PetscReal error_tol_sqr = tol*tol;
1083:   PetscReal phibasis[8], jacobian[9], ijacobian[9], volume;
1084:   PetscReal phypts[3] = {0.0, 0.0, 0.0};
1085:   PetscReal delta[3] = {0.0, 0.0, 0.0};
1086:   PetscErrorCode  ierr;
1087:   PetscInt iters=0;
1088:   PetscReal error=1.0;


1095:   PetscMemzero(jacobian, dim * dim * sizeof(PetscReal));
1096:   PetscMemzero(ijacobian, dim * dim * sizeof(PetscReal));
1097:   PetscMemzero(phibasis, nverts * sizeof(PetscReal));

1099:   /* zero initial guess */
1100:   natparam[0] = natparam[1] = natparam[2] = 0.0;

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

1105:   error = 0.0;
1106:   switch(dim) {
1107:     case 3:
1108:       delta[2] = phypts[2] - xphy[2];
1109:       error += (delta[2]*delta[2]);
1110:     case 2:
1111:       delta[1] = phypts[1] - xphy[1];
1112:       error += (delta[1]*delta[1]);
1113:     case 1:
1114:       delta[0] = phypts[0] - xphy[0];
1115:       error += (delta[0]*delta[0]);
1116:       break;
1117:   }

1119: #if 0
1120:   PetscInfo4(NULL, "Current point in physical space : (%g, %g, %g) and error : %g\n", xphy[0], xphy[1], xphy[2], error);
1121: #endif

1123:   while (error > error_tol_sqr) {

1125:     if(++iters > max_iterations)
1126:       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]);

1128:     /* Compute natparam -= J.inverse() * delta */
1129:     switch(dim) {
1130:       case 1:
1131:         natparam[0] -= ijacobian[0] * delta[0];
1132:         break;
1133:       case 2:
1134:         natparam[0] -= ijacobian[0] * delta[0] + ijacobian[1] * delta[1];
1135:         natparam[1] -= ijacobian[2] * delta[0] + ijacobian[3] * delta[1];
1136:         break;
1137:       case 3:
1138:         natparam[0] -= ijacobian[0] * delta[0] + ijacobian[1] * delta[1] + ijacobian[2] * delta[2];
1139:         natparam[1] -= ijacobian[3] * delta[0] + ijacobian[4] * delta[1] + ijacobian[5] * delta[2];
1140:         natparam[2] -= ijacobian[6] * delta[0] + ijacobian[7] * delta[1] + ijacobian[8] * delta[2];
1141:         break;
1142:     }

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

1147:     error = 0.0;
1148:     switch(dim) {
1149:       case 3:
1150:         delta[2] = phypts[2] - xphy[2];
1151:         error += (delta[2]*delta[2]);
1152:       case 2:
1153:         delta[1] = phypts[1] - xphy[1];
1154:         error += (delta[1]*delta[1]);
1155:       case 1:
1156:         delta[0] = phypts[0] - xphy[0];
1157:         error += (delta[0]*delta[0]);
1158:         break;
1159:     }
1160: #if 0
1161:     PetscInfo5(NULL, "Newton [%D] : Current point in reference space : (%g, %g, %g) and error : %g\n", iters, natparam[0], natparam[1], natparam[2], error);
1162: #endif
1163:   }
1164:   if (phi) {
1165:     PetscMemcpy(phi, phibasis, nverts * sizeof(PetscReal));
1166:   }
1167:   return(0);
1168: }