Actual source code: dmmbfem.cxx

petsc-3.13.6 2020-09-29
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: /*@C
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
111: .vb
112:     1-------2        1----3----2
113:       EDGE2             EDGE3
114: .ve

116:   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: @*/
131: PetscErrorCode Compute_Lagrange_Basis_1D_Internal ( const PetscInt nverts, const PetscReal *coords, const PetscInt npts, const PetscReal *quad, PetscReal *phypts,
132:     PetscReal *jxw, PetscReal *phi, PetscReal *dphidx,
133:     PetscReal *jacobian, PetscReal *ijacobian, PetscReal *volume)
134: {
135:   int i, j;
136:   PetscErrorCode  ierr;

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

240:   Notes:

242:   Example Physical Element (QUAD4)
243: .vb
244:     4------3        s
245:     |      |        |
246:     |      |        |
247:     |      |        |
248:     1------2        0-------r
249: .ve

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

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

264:   Level: advanced

266: @*/
267: PetscErrorCode Compute_Lagrange_Basis_2D_Internal ( const PetscInt nverts, const PetscReal *coords, const PetscInt npts, const PetscReal *quad, PetscReal *phypts,
268:     PetscReal *jxw, PetscReal *phi, PetscReal *dphidx, PetscReal *dphidy,
269:     PetscReal *jacobian, PetscReal *ijacobian, PetscReal *volume)
270: {
271:   PetscInt i, j, k;
272:   PetscErrorCode   ierr;

278:   PetscArrayzero(phi, npts);
279:   if (phypts) {
280:     PetscArrayzero(phypts, npts * 3);
281:   }
282:   if (dphidx) { /* Reset arrays. */
283:     PetscArrayzero(dphidx, npts * nverts);
284:     PetscArrayzero(dphidy, npts * nverts);
285:   }
286:   if (nverts == 4) { /* Linear Quadrangle */

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

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

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

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

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

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

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

336:     PetscArrayzero(jacobian, 4);
337:     PetscArrayzero(ijacobian, 4);

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

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

345:     /* invert the jacobian */
346:     DMatrix_Invert_2x2_Internal(jacobian, ijacobian, volume);
347:     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);

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

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

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

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

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

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

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

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


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

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

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

415:   Notes:

417:   Example Physical Element (HEX8)
418: .vb
419:       8------7
420:      /|     /|        t  s
421:     5------6 |        | /
422:     | |    | |        |/
423:     | 4----|-3        0-------r
424:     |/     |/
425:     1------2
426: .ve

428:   Input Parameter:
429: +  PetscInt  nverts -          the number of element vertices
430: .  PetscReal coords[3*nverts] - the physical coordinates of the vertices (in canonical numbering)
431: .  PetscInt  npts -            the number of evaluation points (quadrature points)
432: -  PetscReal quad[3*npts] -    the evaluation points (quadrature points) in the reference space

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

442:   Level: advanced

444: @*/
445: PetscErrorCode Compute_Lagrange_Basis_3D_Internal ( const PetscInt nverts, const PetscReal *coords, const PetscInt npts, const PetscReal *quad, PetscReal *phypts,
446:     PetscReal *jxw, PetscReal *phi, PetscReal *dphidx, PetscReal *dphidy, PetscReal *dphidz,
447:     PetscReal *jacobian, PetscReal *ijacobian, PetscReal *volume)
448: {
449:   PetscInt i, j, k;

456:   /* Reset arrays. */
457:   PetscArrayzero(phi, npts);
458:   if (phypts) {
459:     PetscArrayzero(phypts, npts * 3);
460:   }
461:   if (dphidx) {
462:     PetscArrayzero(dphidx, npts * nverts);
463:     PetscArrayzero(dphidy, npts * nverts);
464:     PetscArrayzero(dphidz, npts * nverts);
465:   }

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

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

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

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

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

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

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

535:       /* invert the jacobian */
536:       DMatrix_Invert_3x3_Internal(jacobian, ijacobian, volume);
537:       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);

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

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

553:     PetscArrayzero(jacobian, 9);
554:     PetscArrayzero(ijacobian, 9);

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

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

563:     /* invert the jacobian */
564:     DMatrix_Invert_3x3_Internal(jacobian, ijacobian, volume);
565:     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);

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

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

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

619:       phi[offset + 0] = 1.0 - r - s - t;
620:       phi[offset + 1] = r;
621:       phi[offset + 2] = s;
622:       phi[offset + 3] = t;

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

633:       /* Now set the derivatives */
634:       if (dphidx) {
635:         dphidx[0 + offset] = Dx[0];
636:         dphidx[1 + offset] = Dx[1];
637:         dphidx[2 + offset] = Dx[2];
638:         dphidx[3 + offset] = Dx[3];

640:         dphidy[0 + offset] = Dy[0];
641:         dphidy[1 + offset] = Dy[1];
642:         dphidy[2 + offset] = Dy[2];
643:         dphidy[3 + offset] = Dy[3];

645:         dphidz[0 + offset] = Dz[0];
646:         dphidz[1 + offset] = Dz[1];
647:         dphidz[2 + offset] = Dz[2];
648:         dphidz[3 + offset] = Dz[3];
649:       }

651:     } /* Tetrahedra -- ends */
652:   }
653:   else
654:   {
655:     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);
656:   }
657: #if 0
658:   /* verify if the computed basis functions are consistent */
659:   for ( j = 0; j < npts; j++ ) {
660:     PetscScalar phisum = 0, dphixsum = 0, dphiysum = 0, dphizsum = 0;
661:     const int offset = j * nverts;
662:     for ( i = 0; i < nverts; i++ ) {
663:       phisum += phi[i + offset];
664:       if (dphidx) dphixsum += dphidx[i + offset];
665:       if (dphidy) dphiysum += dphidy[i + offset];
666:       if (dphidz) dphizsum += dphidz[i + offset];
667:       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]);
668:     }
669:     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);
670:   }
671: #endif
672:   return(0);
673: }



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

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

683:   Input Parameter:
684: +  PetscInt  nverts,           the number of element vertices
685: .  PetscReal coords[3*nverts], the physical coordinates of the vertices (in canonical numbering)
686: .  PetscInt  npts,             the number of evaluation points (quadrature points)
687: -  PetscReal quad[3*npts],     the evaluation points (quadrature points) in the reference space

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

695:   Level: advanced

697: @*/
698: PetscErrorCode DMMoabFEMComputeBasis ( const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscQuadrature quadrature,
699:                                        PetscReal *phypts, PetscReal *jacobian_quadrature_weight_product,
700:                                        PetscReal *fe_basis, PetscReal **fe_basis_derivatives)
701: {
702:   PetscErrorCode  ierr;
703:   PetscInt        npoints,idim;
704:   bool            compute_der;
705:   const PetscReal *quadpts, *quadwts;
706:   PetscReal       jacobian[9], ijacobian[9], volume;

712:   compute_der = (fe_basis_derivatives != NULL);

714:   /* Get the quadrature points and weights for the given quadrature rule */
715:   PetscQuadratureGetData(quadrature, &idim, NULL, &npoints, &quadpts, &quadwts);
716:   if (idim != dim) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Dimension mismatch: provided (%D) vs quadrature (%D)\n",idim,dim);
717:   if (jacobian_quadrature_weight_product) {
718:     PetscArraycpy(jacobian_quadrature_weight_product, quadwts, npoints);
719:   }

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



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

755:   Input Parameter:

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

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

763:   Level: advanced

765: @*/
766: PetscErrorCode DMMoabFEMCreateQuadratureDefault ( const PetscInt dim, const PetscInt nverts, PetscQuadrature *quadrature )
767: {
768:   PetscReal *w, *x;
769:   PetscInt nc=1;
770:   PetscErrorCode  ierr;

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

827:         w[0] = w[1] = w[2] = w[3] = 1.0 / 24.0;
828:         order = 4;
829:       }
830:       else if (npoints == 10) { /*  KEAST3, K3  N=10, O=10 */
831:         const PetscReal x_10[30] = { 0.5684305841968444, 0.1438564719343852, 0.1438564719343852,
832:                                      0.1438564719343852, 0.1438564719343852, 0.1438564719343852,
833:                                      0.1438564719343852, 0.1438564719343852, 0.5684305841968444,
834:                                      0.1438564719343852, 0.5684305841968444, 0.1438564719343852,
835:                                      0.0000000000000000, 0.5000000000000000, 0.5000000000000000,
836:                                      0.5000000000000000, 0.0000000000000000, 0.5000000000000000,
837:                                      0.5000000000000000, 0.5000000000000000, 0.0000000000000000,
838:                                      0.5000000000000000, 0.0000000000000000, 0.0000000000000000,
839:                                      0.0000000000000000, 0.5000000000000000, 0.0000000000000000,
840:                                      0.0000000000000000, 0.0000000000000000, 0.5000000000000000
841:                                    };
842:         PetscArraycpy(x, x_10, 30);

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

866: /* Compute Jacobians */
867: PetscErrorCode ComputeJacobian_Internal ( const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal *quad, PetscReal *phypts,
868:   PetscReal *jacobian, PetscReal *ijacobian, PetscReal* dvolume)
869: {
870:   PetscInt i;
871:   PetscReal volume=1.0;

878:   PetscArrayzero(jacobian, dim * dim);
879:   if (ijacobian) {
880:     PetscArrayzero(ijacobian, dim * dim);
881:   }
882:   if (phypts) {
883:     PetscArrayzero(phypts, /*npts=1 * */ 3);
884:   }

886:   if (dim == 1) {

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

892:       for (i = 0; i < nverts; ++i) {
893:         const PetscReal* vertices = coordinates + i * 3;
894:         jacobian[0] += dNi_dxi[i] * vertices[0];
895:       }
896:     }
897:     else if (nverts == 3) { /* Quadratic Edge */
898:       const PetscReal dNi_dxi[3]  = { 4 * r - 3.0, 4 * ( 1.0 - 2.0 * r ), 4.0 * r - 1.0};

900:       for (i = 0; i < nverts; ++i) {
901:         const PetscReal* vertices = coordinates + i * 3;
902:         jacobian[0] += dNi_dxi[i] * vertices[0];
903:       }
904:     }
905:     else {
906:       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);
907:     }

909:     if (ijacobian) {
910:       /* invert the jacobian */
911:       ijacobian[0] = 1.0 / jacobian[0];
912:     }

914:   }
915:   else if (dim == 2) {

917:     if (nverts == 4) { /* Linear Quadrangle */
918:       const PetscReal& r = quad[0];
919:       const PetscReal& s = quad[1];

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

924:       for (i = 0; i < nverts; ++i) {
925:         const PetscReal* vertices = coordinates + i * 3;
926:         jacobian[0] += dNi_dxi[i]  * vertices[0];
927:         jacobian[2] += dNi_dxi[i]  * vertices[1];
928:         jacobian[1] += dNi_deta[i] * vertices[0];
929:         jacobian[3] += dNi_deta[i] * vertices[1];
930:       }
931:     }
932:     else if (nverts == 3) { /* Linear triangle */
933:       const PetscReal x2 = coordinates[2 * 3 + 0], y2 = coordinates[2 * 3 + 1];

935:       /* Jacobian is constant */
936:       jacobian[0] = (coordinates[0 * 3 + 0] - x2); jacobian[1] = (coordinates[1 * 3 + 0] - x2);
937:       jacobian[2] = (coordinates[0 * 3 + 1] - y2); jacobian[3] = (coordinates[1 * 3 + 1] - y2);
938:     }
939:     else {
940:       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);
941:     }

943:     /* invert the jacobian */
944:     if (ijacobian) {
945:       DMatrix_Invert_2x2_Internal(jacobian, ijacobian, &volume);
946:     }

948:   }
949:   else {

951:     if (nverts == 8) { /* Linear Hexahedra */
952:       const PetscReal& r = quad[0];
953:       const PetscReal& s = quad[1];
954:       const PetscReal& t = quad[2];
955:       const PetscReal dNi_dxi[8]  = {- ( 1.0 - s ) * ( 1.0 - t ),
956:                                        ( 1.0 - s ) * ( 1.0 - t ),
957:                                        (       s ) * ( 1.0 - t ),
958:                                      - (       s ) * ( 1.0 - t ),
959:                                      - ( 1.0 - s ) * (       t ),
960:                                        ( 1.0 - s ) * (       t ),
961:                                        (       s ) * (       t ),
962:                                      - (       s ) * (       t )
963:                                     };

965:       const PetscReal dNi_deta[8]  = { - ( 1.0 - r ) * ( 1.0 - t ),
966:                                        - (       r ) * ( 1.0 - t ),
967:                                          (       r ) * ( 1.0 - t ),
968:                                          ( 1.0 - r ) * ( 1.0 - t ),
969:                                        - ( 1.0 - r ) * (       t ),
970:                                        - (       r ) * (       t ),
971:                                          (       r ) * (       t ),
972:                                          ( 1.0 - r ) * (       t )
973:                                       };

975:       const PetscReal dNi_dzeta[8]  = { - ( 1.0 - r ) * ( 1.0 - s ),
976:                                         - (       r ) * ( 1.0 - s ),
977:                                         - (       r ) * (       s ),
978:                                         - ( 1.0 - r ) * (       s ),
979:                                           ( 1.0 - r ) * ( 1.0 - s ),
980:                                           (       r ) * ( 1.0 - s ),
981:                                           (       r ) * (       s ),
982:                                           ( 1.0 - r ) * (       s )
983:                                      };

985:       for (i = 0; i < nverts; ++i) {
986:         const PetscReal* vertex = coordinates + i * 3;
987:         jacobian[0] += dNi_dxi[i]   * vertex[0];
988:         jacobian[3] += dNi_dxi[i]   * vertex[1];
989:         jacobian[6] += dNi_dxi[i]   * vertex[2];
990:         jacobian[1] += dNi_deta[i]  * vertex[0];
991:         jacobian[4] += dNi_deta[i]  * vertex[1];
992:         jacobian[7] += dNi_deta[i]  * vertex[2];
993:         jacobian[2] += dNi_dzeta[i] * vertex[0];
994:         jacobian[5] += dNi_dzeta[i] * vertex[1];
995:         jacobian[8] += dNi_dzeta[i] * vertex[2];
996:       }

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

1002:       /* compute the jacobian */
1003:       jacobian[0] = coordinates[1 * 3 + 0] - x0;  jacobian[1] = coordinates[2 * 3 + 0] - x0; jacobian[2] = coordinates[3 * 3 + 0] - x0;
1004:       jacobian[3] = coordinates[1 * 3 + 1] - y0;  jacobian[4] = coordinates[2 * 3 + 1] - y0; jacobian[5] = coordinates[3 * 3 + 1] - y0;
1005:       jacobian[6] = coordinates[1 * 3 + 2] - z0;  jacobian[7] = coordinates[2 * 3 + 2] - z0; jacobian[8] = coordinates[3 * 3 + 2] - z0;
1006:     } /* Tetrahedra -- ends */
1007:     else
1008:     {
1009:       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);
1010:     }

1012:     if (ijacobian) {
1013:       /* invert the jacobian */
1014:       DMatrix_Invert_3x3_Internal(jacobian, ijacobian, &volume);
1015:     }

1017:   }
1018:   if ( volume < 1e-12 ) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Element has zero volume: %g. Degenerate element or invalid connectivity\n", volume);
1019:   if (dvolume) *dvolume = volume;
1020:   return(0);
1021: }


1024: PetscErrorCode FEMComputeBasis_JandF ( const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal *quadrature, PetscReal *phypts,
1025:                                        PetscReal *phibasis, PetscReal *jacobian, PetscReal *ijacobian, PetscReal* volume  )
1026: {
1027:   PetscErrorCode  ierr;

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



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

1056:   Input Parameter:
1057: +  PetscInt  dim -         the element dimension (1=EDGE, 2=QUAD/TRI, 3=HEX/TET)
1058: .  PetscInt nverts -       the number of vertices in the physical element
1059: .  PetscReal coordinates - the coordinates of vertices in the physical element
1060: -  PetscReal[3] xphy -     the coordinates of physical point for which natural coordinates (in reference frame) are sought

1062:   Output Parameter:
1063: +  PetscReal[3] natparam - the natural coordinates (in reference frame) corresponding to xphy
1064: -  PetscReal[nverts] phi - the basis functions evaluated at the natural coordinates (natparam)

1066:   Level: advanced

1068: @*/
1069: PetscErrorCode DMMoabPToRMapping( const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal* xphy, PetscReal* natparam, PetscReal* phi)
1070: {
1071:   /* Perform inverse evaluation for the mapping with use of Newton Raphson iteration */
1072:   const PetscReal tol = 1.0e-10;
1073:   const PetscInt max_iterations = 10;
1074:   const PetscReal error_tol_sqr = tol*tol;
1075:   PetscReal phibasis[8], jacobian[9], ijacobian[9], volume;
1076:   PetscReal phypts[3] = {0.0, 0.0, 0.0};
1077:   PetscReal delta[3] = {0.0, 0.0, 0.0};
1078:   PetscErrorCode  ierr;
1079:   PetscInt iters=0;
1080:   PetscReal error=1.0;


1087:   PetscArrayzero(jacobian, dim * dim);
1088:   PetscArrayzero(ijacobian, dim * dim);
1089:   PetscArrayzero(phibasis, nverts);

1091:   /* zero initial guess */
1092:   natparam[0] = natparam[1] = natparam[2] = 0.0;

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

1097:   error = 0.0;
1098:   switch(dim) {
1099:     case 3:
1100:       delta[2] = phypts[2] - xphy[2];
1101:       error += (delta[2]*delta[2]);
1102:     case 2:
1103:       delta[1] = phypts[1] - xphy[1];
1104:       error += (delta[1]*delta[1]);
1105:     case 1:
1106:       delta[0] = phypts[0] - xphy[0];
1107:       error += (delta[0]*delta[0]);
1108:       break;
1109:   }

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

1115:   while (error > error_tol_sqr) {

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

1120:     /* Compute natparam -= J.inverse() * delta */
1121:     switch(dim) {
1122:       case 1:
1123:         natparam[0] -= ijacobian[0] * delta[0];
1124:         break;
1125:       case 2:
1126:         natparam[0] -= ijacobian[0] * delta[0] + ijacobian[1] * delta[1];
1127:         natparam[1] -= ijacobian[2] * delta[0] + ijacobian[3] * delta[1];
1128:         break;
1129:       case 3:
1130:         natparam[0] -= ijacobian[0] * delta[0] + ijacobian[1] * delta[1] + ijacobian[2] * delta[2];
1131:         natparam[1] -= ijacobian[3] * delta[0] + ijacobian[4] * delta[1] + ijacobian[5] * delta[2];
1132:         natparam[2] -= ijacobian[6] * delta[0] + ijacobian[7] * delta[1] + ijacobian[8] * delta[2];
1133:         break;
1134:     }

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

1139:     error = 0.0;
1140:     switch(dim) {
1141:       case 3:
1142:         delta[2] = phypts[2] - xphy[2];
1143:         error += (delta[2]*delta[2]);
1144:       case 2:
1145:         delta[1] = phypts[1] - xphy[1];
1146:         error += (delta[1]*delta[1]);
1147:       case 1:
1148:         delta[0] = phypts[0] - xphy[0];
1149:         error += (delta[0]*delta[0]);
1150:         break;
1151:     }
1152: #if 0
1153:     PetscInfo5(NULL, "Newton [%D] : Current point in reference space : (%g, %g, %g) and error : %g\n", iters, natparam[0], natparam[1], natparam[2], error);
1154: #endif
1155:   }
1156:   if (phi) {
1157:     PetscArraycpy(phi, phibasis, nverts);
1158:   }
1159:   return(0);
1160: }