Actual source code: dmmbfem.cxx
2: #include <petscconf.h>
3: #include <petscdt.h>
4: #include <petsc/private/dmmbimpl.h>
6: /* Utility functions */
7: static inline PetscReal DMatrix_Determinant_2x2_Internal (const PetscReal inmat[2 * 2])
8: {
9: return inmat[0] * inmat[3] - inmat[1] * inmat[2];
10: }
12: static inline PetscErrorCode DMatrix_Invert_2x2_Internal(const PetscReal *inmat, PetscReal *outmat, PetscReal *determinant)
13: {
14: PetscReal det = DMatrix_Determinant_2x2_Internal(inmat);
15: if (outmat) {
16: outmat[0] = inmat[3] / det;
17: outmat[1] = -inmat[1] / det;
18: outmat[2] = -inmat[2] / det;
19: outmat[3] = inmat[0] / det;
20: }
21: if (determinant) *determinant = det;
22: return 0;
23: }
25: static inline PetscReal DMatrix_Determinant_3x3_Internal(const PetscReal inmat[3 * 3])
26: {
27: return inmat[0] * (inmat[8] * inmat[4] - inmat[7] * inmat[5])
28: - inmat[3] * (inmat[8] * inmat[1] - inmat[7] * inmat[2])
29: + inmat[6] * (inmat[5] * inmat[1] - inmat[4] * inmat[2]);
30: }
32: static inline PetscErrorCode DMatrix_Invert_3x3_Internal (const PetscReal *inmat, PetscReal *outmat, PetscScalar *determinant)
33: {
34: PetscReal det = DMatrix_Determinant_3x3_Internal(inmat);
35: if (outmat) {
36: outmat[0] = (inmat[8] * inmat[4] - inmat[7] * inmat[5]) / det;
37: outmat[1] = -(inmat[8] * inmat[1] - inmat[7] * inmat[2]) / det;
38: outmat[2] = (inmat[5] * inmat[1] - inmat[4] * inmat[2]) / det;
39: outmat[3] = -(inmat[8] * inmat[3] - inmat[6] * inmat[5]) / det;
40: outmat[4] = (inmat[8] * inmat[0] - inmat[6] * inmat[2]) / det;
41: outmat[5] = -(inmat[5] * inmat[0] - inmat[3] * inmat[2]) / det;
42: outmat[6] = (inmat[7] * inmat[3] - inmat[6] * inmat[4]) / det;
43: outmat[7] = -(inmat[7] * inmat[0] - inmat[6] * inmat[1]) / det;
44: outmat[8] = (inmat[4] * inmat[0] - inmat[3] * inmat[1]) / det;
45: }
46: if (determinant) *determinant = det;
47: return 0;
48: }
50: inline PetscReal DMatrix_Determinant_4x4_Internal(PetscReal inmat[4 * 4])
51: {
52: return
53: inmat[0 + 0 * 4] * (
54: inmat[1 + 1 * 4] * (inmat[2 + 2 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 2 * 4])
55: - inmat[1 + 2 * 4] * (inmat[2 + 1 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 1 * 4])
56: + inmat[1 + 3 * 4] * (inmat[2 + 1 * 4] * inmat[3 + 2 * 4] - inmat[2 + 2 * 4] * inmat[3 + 1 * 4]))
57: - inmat[0 + 1 * 4] * (
58: inmat[1 + 0 * 4] * (inmat[2 + 2 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 2 * 4])
59: - inmat[1 + 2 * 4] * (inmat[2 + 0 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 0 * 4])
60: + inmat[1 + 3 * 4] * (inmat[2 + 0 * 4] * inmat[3 + 2 * 4] - inmat[2 + 2 * 4] * inmat[3 + 0 * 4]))
61: + inmat[0 + 2 * 4] * (
62: inmat[1 + 0 * 4] * (inmat[2 + 1 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 1 * 4])
63: - inmat[1 + 1 * 4] * (inmat[2 + 0 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 0 * 4])
64: + inmat[1 + 3 * 4] * (inmat[2 + 0 * 4] * inmat[3 + 1 * 4] - inmat[2 + 1 * 4] * inmat[3 + 0 * 4]))
65: - inmat[0 + 3 * 4] * (
66: inmat[1 + 0 * 4] * (inmat[2 + 1 * 4] * inmat[3 + 2 * 4] - inmat[2 + 2 * 4] * inmat[3 + 1 * 4])
67: - inmat[1 + 1 * 4] * (inmat[2 + 0 * 4] * inmat[3 + 2 * 4] - inmat[2 + 2 * 4] * inmat[3 + 0 * 4])
68: + inmat[1 + 2 * 4] * (inmat[2 + 0 * 4] * inmat[3 + 1 * 4] - inmat[2 + 1 * 4] * inmat[3 + 0 * 4]));
69: }
71: inline PetscErrorCode DMatrix_Invert_4x4_Internal (PetscReal *inmat, PetscReal *outmat, PetscScalar *determinant)
72: {
73: PetscReal det = DMatrix_Determinant_4x4_Internal(inmat);
74: if (outmat) {
75: outmat[0] = (inmat[5] * inmat[10] * inmat[15] + inmat[6] * inmat[11] * inmat[13] + inmat[7] * inmat[9] * inmat[14] - inmat[5] * inmat[11] * inmat[14] - inmat[6] * inmat[9] * inmat[15] - inmat[7] * inmat[10] * inmat[13]) / det;
76: outmat[1] = (inmat[1] * inmat[11] * inmat[14] + inmat[2] * inmat[9] * inmat[15] + inmat[3] * inmat[10] * inmat[13] - inmat[1] * inmat[10] * inmat[15] - inmat[2] * inmat[11] * inmat[13] - inmat[3] * inmat[9] * inmat[14]) / det;
77: outmat[2] = (inmat[1] * inmat[6] * inmat[15] + inmat[2] * inmat[7] * inmat[13] + inmat[3] * inmat[5] * inmat[14] - inmat[1] * inmat[7] * inmat[14] - inmat[2] * inmat[5] * inmat[15] - inmat[3] * inmat[6] * inmat[13]) / det;
78: outmat[3] = (inmat[1] * inmat[7] * inmat[10] + inmat[2] * inmat[5] * inmat[11] + inmat[3] * inmat[6] * inmat[9] - inmat[1] * inmat[6] * inmat[11] - inmat[2] * inmat[7] * inmat[9] - inmat[3] * inmat[5] * inmat[10]) / det;
79: outmat[4] = (inmat[4] * inmat[11] * inmat[14] + inmat[6] * inmat[8] * inmat[15] + inmat[7] * inmat[10] * inmat[12] - inmat[4] * inmat[10] * inmat[15] - inmat[6] * inmat[11] * inmat[12] - inmat[7] * inmat[8] * inmat[14]) / det;
80: outmat[5] = (inmat[0] * inmat[10] * inmat[15] + inmat[2] * inmat[11] * inmat[12] + inmat[3] * inmat[8] * inmat[14] - inmat[0] * inmat[11] * inmat[14] - inmat[2] * inmat[8] * inmat[15] - inmat[3] * inmat[10] * inmat[12]) / det;
81: outmat[6] = (inmat[0] * inmat[7] * inmat[14] + inmat[2] * inmat[4] * inmat[15] + inmat[3] * inmat[6] * inmat[12] - inmat[0] * inmat[6] * inmat[15] - inmat[2] * inmat[7] * inmat[12] - inmat[3] * inmat[4] * inmat[14]) / det;
82: outmat[7] = (inmat[0] * inmat[6] * inmat[11] + inmat[2] * inmat[7] * inmat[8] + inmat[3] * inmat[4] * inmat[10] - inmat[0] * inmat[7] * inmat[10] - inmat[2] * inmat[4] * inmat[11] - inmat[3] * inmat[6] * inmat[8]) / det;
83: outmat[8] = (inmat[4] * inmat[9] * inmat[15] + inmat[5] * inmat[11] * inmat[12] + inmat[7] * inmat[8] * inmat[13] - inmat[4] * inmat[11] * inmat[13] - inmat[5] * inmat[8] * inmat[15] - inmat[7] * inmat[9] * inmat[12]) / det;
84: outmat[9] = (inmat[0] * inmat[11] * inmat[13] + inmat[1] * inmat[8] * inmat[15] + inmat[3] * inmat[9] * inmat[12] - inmat[0] * inmat[9] * inmat[15] - inmat[1] * inmat[11] * inmat[12] - inmat[3] * inmat[8] * inmat[13]) / det;
85: outmat[10] = (inmat[0] * inmat[5] * inmat[15] + inmat[1] * inmat[7] * inmat[12] + inmat[3] * inmat[4] * inmat[13] - inmat[0] * inmat[7] * inmat[13] - inmat[1] * inmat[4] * inmat[15] - inmat[3] * inmat[5] * inmat[12]) / det;
86: outmat[11] = (inmat[0] * inmat[7] * inmat[9] + inmat[1] * inmat[4] * inmat[11] + inmat[3] * inmat[5] * inmat[8] - inmat[0] * inmat[5] * inmat[11] - inmat[1] * inmat[7] * inmat[8] - inmat[3] * inmat[4] * inmat[9]) / det;
87: outmat[12] = (inmat[4] * inmat[10] * inmat[13] + inmat[5] * inmat[8] * inmat[14] + inmat[6] * inmat[9] * inmat[12] - inmat[4] * inmat[9] * inmat[14] - inmat[5] * inmat[10] * inmat[12] - inmat[6] * inmat[8] * inmat[13]) / det;
88: outmat[13] = (inmat[0] * inmat[9] * inmat[14] + inmat[1] * inmat[10] * inmat[12] + inmat[2] * inmat[8] * inmat[13] - inmat[0] * inmat[10] * inmat[13] - inmat[1] * inmat[8] * inmat[14] - inmat[2] * inmat[9] * inmat[12]) / det;
89: outmat[14] = (inmat[0] * inmat[6] * inmat[13] + inmat[1] * inmat[4] * inmat[14] + inmat[2] * inmat[5] * inmat[12] - inmat[0] * inmat[5] * inmat[14] - inmat[1] * inmat[6] * inmat[12] - inmat[2] * inmat[4] * inmat[13]) / det;
90: outmat[15] = (inmat[0] * inmat[5] * inmat[10] + inmat[1] * inmat[6] * inmat[8] + inmat[2] * inmat[4] * inmat[9] - inmat[0] * inmat[6] * inmat[9] - inmat[1] * inmat[4] * inmat[10] - inmat[2] * inmat[5] * inmat[8]) / det;
91: }
92: if (determinant) *determinant = det;
93: return 0;
94: }
96: /*@C
97: Compute_Lagrange_Basis_1D_Internal - Evaluate bases and derivatives at quadrature points for a EDGE2 or EDGE3 element.
99: The routine is given the coordinates of the vertices of a linear or quadratic edge element.
101: The routine evaluates the basis functions associated with each quadrature point provided,
102: and their derivatives with respect to X.
104: Notes:
106: Example Physical Element
107: .vb
108: 1-------2 1----3----2
109: EDGE2 EDGE3
110: .ve
112: Input Parameters:
113: + PetscInt nverts - the number of element vertices
114: . PetscReal coords[3*nverts] - the physical coordinates of the vertices (in canonical numbering)
115: . PetscInt npts - the number of evaluation points (quadrature points)
116: - PetscReal quad[3*npts] - the evaluation points (quadrature points) in the reference space
118: Output Parameters:
119: + PetscReal phypts[3*npts] - the evaluation points (quadrature points) transformed to the physical space
120: . PetscReal jxw[npts] - the jacobian determinant * quadrature weight necessary for assembling discrete contributions
121: . PetscReal phi[npts] - the bases evaluated at the specified quadrature points
122: . PetscReal dphidx[npts] - the derivative of the bases wrt X-direction evaluated at the specified quadrature points
123: . jacobian - jacobian
124: . ijacobian - ijacobian
125: - volume - volume
127: Level: advanced
129: @*/
130: PetscErrorCode Compute_Lagrange_Basis_1D_Internal(const PetscInt nverts, const PetscReal *coords, const PetscInt npts, const PetscReal *quad, PetscReal *phypts,
131: PetscReal *jxw, PetscReal *phi, PetscReal *dphidx,
132: PetscReal *jacobian, PetscReal *ijacobian, PetscReal *volume)
133: {
134: int i, j;
139: if (phypts) {
140: PetscArrayzero(phypts, npts * 3);
141: }
142: if (dphidx) { /* Reset arrays. */
143: PetscArrayzero(dphidx, npts * nverts);
144: }
145: if (nverts == 2) { /* Linear Edge */
147: for (j = 0; j < npts; j++) {
148: const PetscInt offset = j * nverts;
149: const PetscReal r = quad[j];
151: phi[0 + offset] = ( 1.0 - r);
152: phi[1 + offset] = ( r);
154: const PetscReal dNi_dxi[2] = { -1.0, 1.0 };
156: jacobian[0] = ijacobian[0] = volume[0] = 0.0;
157: for (i = 0; i < nverts; ++i) {
158: const PetscReal* vertices = coords + i * 3;
159: jacobian[0] += dNi_dxi[i] * vertices[0];
160: if (phypts) {
161: phypts[3 * j + 0] += phi[i + offset] * vertices[0];
162: }
163: }
165: /* invert the jacobian */
166: *volume = jacobian[0];
167: ijacobian[0] = 1.0 / jacobian[0];
168: jxw[j] *= *volume;
170: /* Divide by element jacobian. */
171: for (i = 0; i < nverts; i++) {
172: if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[0];
173: }
174: }
175: } else if (nverts == 3) { /* Quadratic Edge */
177: for (j = 0; j < npts; j++) {
178: const PetscInt offset = j * nverts;
179: const PetscReal r = quad[j];
181: phi[0 + offset] = 1.0 + r * ( 2.0 * r - 3.0);
182: phi[1 + offset] = 4.0 * r * ( 1.0 - r);
183: phi[2 + offset] = r * ( 2.0 * r - 1.0);
185: const PetscReal dNi_dxi[3] = { 4 * r - 3.0, 4 * ( 1.0 - 2.0 * r), 4.0 * r - 1.0};
187: jacobian[0] = ijacobian[0] = volume[0] = 0.0;
188: for (i = 0; i < nverts; ++i) {
189: const PetscReal* vertices = coords + i * 3;
190: jacobian[0] += dNi_dxi[i] * vertices[0];
191: if (phypts) {
192: phypts[3 * j + 0] += phi[i + offset] * vertices[0];
193: }
194: }
196: /* invert the jacobian */
197: *volume = jacobian[0];
198: ijacobian[0] = 1.0 / jacobian[0];
199: if (jxw) jxw[j] *= *volume;
201: /* Divide by element jacobian. */
202: for (i = 0; i < nverts; i++) {
203: if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[0];
204: }
205: }
206: } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "The number of entity vertices are invalid. Currently only support EDGE2 and EDGE3 basis evaluations in 1-D : %D", nverts);
207: return 0;
208: }
210: /*@C
211: Compute_Lagrange_Basis_2D_Internal - Evaluate bases and derivatives at quadrature points for a QUAD4 or TRI3 element.
213: The routine is given the coordinates of the vertices of a quadrangle or triangle.
215: The routine evaluates the basis functions associated with each quadrature point provided,
216: and their derivatives with respect to X and Y.
218: Notes:
220: Example Physical Element (QUAD4)
221: .vb
222: 4------3 s
223: | | |
224: | | |
225: | | |
226: 1------2 0-------r
227: .ve
229: Input Parameters:
230: + PetscInt nverts - the number of element vertices
231: . PetscReal coords[3*nverts] - the physical coordinates of the vertices (in canonical numbering)
232: . PetscInt npts - the number of evaluation points (quadrature points)
233: - PetscReal quad[3*npts] - the evaluation points (quadrature points) in the reference space
235: Output Parameters:
236: + PetscReal phypts[3*npts] - the evaluation points (quadrature points) transformed to the physical space
237: . PetscReal jxw[npts] - the jacobian determinant * quadrature weight necessary for assembling discrete contributions
238: . PetscReal phi[npts] - the bases evaluated at the specified quadrature points
239: . PetscReal dphidx[npts] - the derivative of the bases wrt X-direction evaluated at the specified quadrature points
240: . PetscReal dphidy[npts] - the derivative of the bases wrt Y-direction evaluated at the specified quadrature points
241: . jacobian - jacobian
242: . ijacobian - ijacobian
243: - volume - volume
245: Level: advanced
247: @*/
248: PetscErrorCode Compute_Lagrange_Basis_2D_Internal(const PetscInt nverts, const PetscReal *coords, const PetscInt npts, const PetscReal *quad, PetscReal *phypts,
249: PetscReal *jxw, PetscReal *phi, PetscReal *dphidx, PetscReal *dphidy,
250: PetscReal *jacobian, PetscReal *ijacobian, PetscReal *volume)
251: {
252: PetscInt i, j, k;
257: PetscArrayzero(phi, npts);
258: if (phypts) {
259: PetscArrayzero(phypts, npts * 3);
260: }
261: if (dphidx) { /* Reset arrays. */
262: PetscArrayzero(dphidx, npts * nverts);
263: PetscArrayzero(dphidy, npts * nverts);
264: }
265: if (nverts == 4) { /* Linear Quadrangle */
267: for (j = 0; j < npts; j++)
268: {
269: const PetscInt offset = j * nverts;
270: const PetscReal r = quad[0 + j * 2];
271: const PetscReal s = quad[1 + j * 2];
273: phi[0 + offset] = ( 1.0 - r) * ( 1.0 - s);
274: phi[1 + offset] = r * ( 1.0 - s);
275: phi[2 + offset] = r * s;
276: phi[3 + offset] = ( 1.0 - r) * s;
278: const PetscReal dNi_dxi[4] = { -1.0 + s, 1.0 - s, s, -s };
279: const PetscReal dNi_deta[4] = { -1.0 + r, -r, r, 1.0 - r };
281: PetscArrayzero(jacobian, 4);
282: PetscArrayzero(ijacobian, 4);
283: for (i = 0; i < nverts; ++i) {
284: const PetscReal* vertices = coords + i * 3;
285: jacobian[0] += dNi_dxi[i] * vertices[0];
286: jacobian[2] += dNi_dxi[i] * vertices[1];
287: jacobian[1] += dNi_deta[i] * vertices[0];
288: jacobian[3] += dNi_deta[i] * vertices[1];
289: if (phypts) {
290: phypts[3 * j + 0] += phi[i + offset] * vertices[0];
291: phypts[3 * j + 1] += phi[i + offset] * vertices[1];
292: phypts[3 * j + 2] += phi[i + offset] * vertices[2];
293: }
294: }
296: /* invert the jacobian */
297: DMatrix_Invert_2x2_Internal(jacobian, ijacobian, volume);
300: if (jxw) jxw[j] *= *volume;
302: /* Let us compute the bases derivatives by scaling with inverse jacobian. */
303: if (dphidx) {
304: for (i = 0; i < nverts; i++) {
305: for (k = 0; k < 2; ++k) {
306: if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[k * 2 + 0];
307: if (dphidy) dphidy[i + offset] += dNi_deta[i] * ijacobian[k * 2 + 1];
308: } /* for k=[0..2) */
309: } /* for i=[0..nverts) */
310: } /* if (dphidx) */
311: }
312: } else if (nverts == 3) { /* Linear triangle */
313: const PetscReal x2 = coords[2 * 3 + 0], y2 = coords[2 * 3 + 1];
315: PetscArrayzero(jacobian, 4);
316: PetscArrayzero(ijacobian, 4);
318: /* Jacobian is constant */
319: jacobian[0] = (coords[0 * 3 + 0] - x2); jacobian[1] = (coords[1 * 3 + 0] - x2);
320: jacobian[2] = (coords[0 * 3 + 1] - y2); jacobian[3] = (coords[1 * 3 + 1] - y2);
322: /* invert the jacobian */
323: DMatrix_Invert_2x2_Internal(jacobian, ijacobian, volume);
326: const PetscReal Dx[3] = { ijacobian[0], ijacobian[2], - ijacobian[0] - ijacobian[2] };
327: const PetscReal Dy[3] = { ijacobian[1], ijacobian[3], - ijacobian[1] - ijacobian[3] };
329: for (j = 0; j < npts; j++) {
330: const PetscInt offset = j * nverts;
331: const PetscReal r = quad[0 + j * 2];
332: const PetscReal s = quad[1 + j * 2];
334: if (jxw) jxw[j] *= 0.5;
336: /* const PetscReal Ni[3] = { r, s, 1.0 - r - s }; */
337: const PetscReal phipts_x = coords[2 * 3 + 0] + jacobian[0] * r + jacobian[1] * s;
338: const PetscReal phipts_y = coords[2 * 3 + 1] + jacobian[2] * r + jacobian[3] * s;
339: if (phypts) {
340: phypts[offset + 0] = phipts_x;
341: phypts[offset + 1] = phipts_y;
342: }
344: /* \phi_0 = (b.y - c.y) x + (b.x - c.x) y + c.x b.y - b.x c.y */
345: phi[0 + offset] = ( ijacobian[0] * (phipts_x - x2) + ijacobian[1] * (phipts_y - y2));
346: /* \phi_1 = (c.y - a.y) x + (a.x - c.x) y + c.x a.y - a.x c.y */
347: phi[1 + offset] = ( ijacobian[2] * (phipts_x - x2) + ijacobian[3] * (phipts_y - y2));
348: phi[2 + offset] = 1.0 - phi[0 + offset] - phi[1 + offset];
350: if (dphidx) {
351: dphidx[0 + offset] = Dx[0];
352: dphidx[1 + offset] = Dx[1];
353: dphidx[2 + offset] = Dx[2];
354: }
356: if (dphidy) {
357: dphidy[0 + offset] = Dy[0];
358: dphidy[1 + offset] = Dy[1];
359: dphidy[2 + offset] = Dy[2];
360: }
362: }
363: } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "The number of entity vertices are invalid. Currently only support QUAD4 and TRI3 basis evaluations in 2-D : %D", nverts);
364: return 0;
365: }
367: /*@C
368: Compute_Lagrange_Basis_3D_Internal - Evaluate bases and derivatives at quadrature points for a HEX8 or TET4 element.
370: The routine is given the coordinates of the vertices of a hexahedra or tetrahedra.
372: The routine evaluates the basis functions associated with each quadrature point provided,
373: and their derivatives with respect to X, Y and Z.
375: Notes:
377: Example Physical Element (HEX8)
378: .vb
379: 8------7
380: /| /| t s
381: 5------6 | | /
382: | | | | |/
383: | 4----|-3 0-------r
384: |/ |/
385: 1------2
386: .ve
388: Input Parameters:
389: + PetscInt nverts - the number of element vertices
390: . PetscReal coords[3*nverts] - the physical coordinates of the vertices (in canonical numbering)
391: . PetscInt npts - the number of evaluation points (quadrature points)
392: - PetscReal quad[3*npts] - the evaluation points (quadrature points) in the reference space
394: Output Parameters:
395: + PetscReal phypts[3*npts] - the evaluation points (quadrature points) transformed to the physical space
396: . PetscReal jxw[npts] - the jacobian determinant * quadrature weight necessary for assembling discrete contributions
397: . PetscReal phi[npts] - the bases evaluated at the specified quadrature points
398: . PetscReal dphidx[npts] - the derivative of the bases wrt X-direction evaluated at the specified quadrature points
399: . PetscReal dphidy[npts] - the derivative of the bases wrt Y-direction evaluated at the specified quadrature points
400: . PetscReal dphidz[npts] - the derivative of the bases wrt Z-direction evaluated at the specified quadrature points
401: . jacobian - jacobian
402: . ijacobian - ijacobian
403: - volume - volume
405: Level: advanced
407: @*/
408: PetscErrorCode Compute_Lagrange_Basis_3D_Internal(const PetscInt nverts, const PetscReal *coords, const PetscInt npts, const PetscReal *quad, PetscReal *phypts,
409: PetscReal *jxw, PetscReal *phi, PetscReal *dphidx, PetscReal *dphidy, PetscReal *dphidz,
410: PetscReal *jacobian, PetscReal *ijacobian, PetscReal *volume)
411: {
412: PetscInt i, j, k;
418: PetscArrayzero(phi, npts);
419: if (phypts) {
420: PetscArrayzero(phypts, npts * 3);
421: }
422: if (dphidx) {
423: PetscArrayzero(dphidx, npts * nverts);
424: PetscArrayzero(dphidy, npts * nverts);
425: PetscArrayzero(dphidz, npts * nverts);
426: }
428: if (nverts == 8) { /* Linear Hexahedra */
430: for (j = 0; j < npts; j++) {
431: const PetscInt offset = j * nverts;
432: const PetscReal& r = quad[j * 3 + 0];
433: const PetscReal& s = quad[j * 3 + 1];
434: const PetscReal& t = quad[j * 3 + 2];
436: phi[offset + 0] = ( 1.0 - r) * ( 1.0 - s) * ( 1.0 - t); /* 0,0,0 */
437: phi[offset + 1] = ( r) * ( 1.0 - s) * ( 1.0 - t); /* 1,0,0 */
438: phi[offset + 2] = ( r) * ( s) * ( 1.0 - t); /* 1,1,0 */
439: phi[offset + 3] = ( 1.0 - r) * ( s) * ( 1.0 - t); /* 0,1,0 */
440: phi[offset + 4] = ( 1.0 - r) * ( 1.0 - s) * ( t); /* 0,0,1 */
441: phi[offset + 5] = ( r) * ( 1.0 - s) * ( t); /* 1,0,1 */
442: phi[offset + 6] = ( r) * ( s) * ( t); /* 1,1,1 */
443: phi[offset + 7] = ( 1.0 - r) * ( s) * ( t); /* 0,1,1 */
445: const PetscReal dNi_dxi[8] = {- ( 1.0 - s) * ( 1.0 - t),
446: ( 1.0 - s) * ( 1.0 - t),
447: ( s) * ( 1.0 - t),
448: - ( s) * ( 1.0 - t),
449: - ( 1.0 - s) * ( t),
450: ( 1.0 - s) * ( t),
451: ( s) * ( t),
452: - ( s) * ( t)
453: };
455: const PetscReal dNi_deta[8] = { - ( 1.0 - r) * ( 1.0 - t),
456: - ( r) * ( 1.0 - t),
457: ( r) * ( 1.0 - t),
458: ( 1.0 - r) * ( 1.0 - t),
459: - ( 1.0 - r) * ( t),
460: - ( r) * ( t),
461: ( r) * ( t),
462: ( 1.0 - r) * ( t)
463: };
465: const PetscReal dNi_dzeta[8] = { - ( 1.0 - r) * ( 1.0 - s),
466: - ( r) * ( 1.0 - s),
467: - ( r) * ( s),
468: - ( 1.0 - r) * ( s),
469: ( 1.0 - r) * ( 1.0 - s),
470: ( r) * ( 1.0 - s),
471: ( r) * ( s),
472: ( 1.0 - r) * ( s)
473: };
475: PetscArrayzero(jacobian, 9);
476: PetscArrayzero(ijacobian, 9);
477: for (i = 0; i < nverts; ++i) {
478: const PetscReal* vertex = coords + i * 3;
479: jacobian[0] += dNi_dxi[i] * vertex[0];
480: jacobian[3] += dNi_dxi[i] * vertex[1];
481: jacobian[6] += dNi_dxi[i] * vertex[2];
482: jacobian[1] += dNi_deta[i] * vertex[0];
483: jacobian[4] += dNi_deta[i] * vertex[1];
484: jacobian[7] += dNi_deta[i] * vertex[2];
485: jacobian[2] += dNi_dzeta[i] * vertex[0];
486: jacobian[5] += dNi_dzeta[i] * vertex[1];
487: jacobian[8] += dNi_dzeta[i] * vertex[2];
488: if (phypts) {
489: phypts[3 * j + 0] += phi[i + offset] * vertex[0];
490: phypts[3 * j + 1] += phi[i + offset] * vertex[1];
491: phypts[3 * j + 2] += phi[i + offset] * vertex[2];
492: }
493: }
495: /* invert the jacobian */
496: DMatrix_Invert_3x3_Internal(jacobian, ijacobian, volume);
499: if (jxw) jxw[j] *= (*volume);
501: /* Divide by element jacobian. */
502: for (i = 0; i < nverts; ++i) {
503: for (k = 0; k < 3; ++k) {
504: if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[0 * 3 + k];
505: if (dphidy) dphidy[i + offset] += dNi_deta[i] * ijacobian[1 * 3 + k];
506: if (dphidz) dphidz[i + offset] += dNi_dzeta[i] * ijacobian[2 * 3 + k];
507: }
508: }
509: }
510: } else if (nverts == 4) { /* Linear Tetrahedra */
511: PetscReal Dx[4]={0,0,0,0}, Dy[4]={0,0,0,0}, Dz[4]={0,0,0,0};
512: const PetscReal x0 = coords[/*0 * 3 +*/ 0], y0 = coords[/*0 * 3 +*/ 1], z0 = coords[/*0 * 3 +*/ 2];
514: PetscArrayzero(jacobian, 9);
515: PetscArrayzero(ijacobian, 9);
517: /* compute the jacobian */
518: jacobian[0] = coords[1 * 3 + 0] - x0; jacobian[1] = coords[2 * 3 + 0] - x0; jacobian[2] = coords[3 * 3 + 0] - x0;
519: jacobian[3] = coords[1 * 3 + 1] - y0; jacobian[4] = coords[2 * 3 + 1] - y0; jacobian[5] = coords[3 * 3 + 1] - y0;
520: jacobian[6] = coords[1 * 3 + 2] - z0; jacobian[7] = coords[2 * 3 + 2] - z0; jacobian[8] = coords[3 * 3 + 2] - z0;
522: /* invert the jacobian */
523: DMatrix_Invert_3x3_Internal(jacobian, ijacobian, volume);
526: if (dphidx) {
527: Dx[0] = ( coords[1 + 2 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3])
528: - coords[1 + 1 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3])
529: - coords[1 + 3 * 3] * ( coords[2 + 1 * 3] - coords[2 + 2 * 3])) / *volume;
530: Dx[1] = - ( coords[1 + 2 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3])
531: - coords[1 + 0 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3])
532: - coords[1 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 2 * 3])) / *volume;
533: Dx[2] = ( coords[1 + 1 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3])
534: - coords[1 + 0 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3])
535: - coords[1 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 1 * 3])) / *volume;
536: Dx[3] = - ( Dx[0] + Dx[1] + Dx[2]);
537: Dy[0] = ( coords[0 + 1 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3])
538: - coords[0 + 2 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3])
539: + coords[0 + 3 * 3] * ( coords[2 + 1 * 3] - coords[2 + 2 * 3])) / *volume;
540: Dy[1] = - ( coords[0 + 0 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3])
541: - coords[0 + 2 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3])
542: + coords[0 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 2 * 3])) / *volume;
543: Dy[2] = ( coords[0 + 0 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3])
544: - coords[0 + 1 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3])
545: + coords[0 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 1 * 3])) / *volume;
546: Dy[3] = - ( Dy[0] + Dy[1] + Dy[2]);
547: Dz[0] = ( coords[0 + 1 * 3] * (coords[1 + 3 * 3] - coords[1 + 2 * 3])
548: - coords[0 + 2 * 3] * (coords[1 + 3 * 3] - coords[1 + 1 * 3])
549: + coords[0 + 3 * 3] * (coords[1 + 2 * 3] - coords[1 + 1 * 3])) / *volume;
550: Dz[1] = - ( coords[0 + 0 * 3] * (coords[1 + 3 * 3] - coords[1 + 2 * 3])
551: + coords[0 + 2 * 3] * (coords[1 + 0 * 3] - coords[1 + 3 * 3])
552: - coords[0 + 3 * 3] * (coords[1 + 0 * 3] - coords[1 + 2 * 3])) / *volume;
553: Dz[2] = ( coords[0 + 0 * 3] * (coords[1 + 3 * 3] - coords[1 + 1 * 3])
554: + coords[0 + 1 * 3] * (coords[1 + 0 * 3] - coords[1 + 3 * 3])
555: - coords[0 + 3 * 3] * (coords[1 + 0 * 3] - coords[1 + 1 * 3])) / *volume;
556: Dz[3] = - ( Dz[0] + Dz[1] + Dz[2]);
557: }
559: for (j = 0; j < npts; j++) {
560: const PetscInt offset = j * nverts;
561: const PetscReal& r = quad[j * 3 + 0];
562: const PetscReal& s = quad[j * 3 + 1];
563: const PetscReal& t = quad[j * 3 + 2];
565: if (jxw) jxw[j] *= *volume;
567: phi[offset + 0] = 1.0 - r - s - t;
568: phi[offset + 1] = r;
569: phi[offset + 2] = s;
570: phi[offset + 3] = t;
572: if (phypts) {
573: for (i = 0; i < nverts; ++i) {
574: const PetscScalar* vertices = coords + i * 3;
575: phypts[3 * j + 0] += phi[i + offset] * vertices[0];
576: phypts[3 * j + 1] += phi[i + offset] * vertices[1];
577: phypts[3 * j + 2] += phi[i + offset] * vertices[2];
578: }
579: }
581: /* Now set the derivatives */
582: if (dphidx) {
583: dphidx[0 + offset] = Dx[0];
584: dphidx[1 + offset] = Dx[1];
585: dphidx[2 + offset] = Dx[2];
586: dphidx[3 + offset] = Dx[3];
588: dphidy[0 + offset] = Dy[0];
589: dphidy[1 + offset] = Dy[1];
590: dphidy[2 + offset] = Dy[2];
591: dphidy[3 + offset] = Dy[3];
593: dphidz[0 + offset] = Dz[0];
594: dphidz[1 + offset] = Dz[1];
595: dphidz[2 + offset] = Dz[2];
596: dphidz[3 + offset] = Dz[3];
597: }
599: } /* Tetrahedra -- ends */
600: } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "The number of entity vertices are invalid. Currently only support HEX8 and TET4 basis evaluations in 3-D : %D", nverts);
601: return 0;
602: }
604: /*@C
605: DMMoabFEMComputeBasis - Evaluate bases and derivatives at quadrature points for a linear EDGE/QUAD/TRI/HEX/TET element.
607: The routine takes the coordinates of the vertices of an element and computes basis functions associated with
608: each quadrature point provided, and their derivatives with respect to X, Y and Z as appropriate.
610: Input Parameters:
611: + PetscInt nverts - the number of element vertices
612: . PetscReal coords[3*nverts] - the physical coordinates of the vertices (in canonical numbering)
613: . PetscInt npts - the number of evaluation points (quadrature points)
614: - PetscReal quad[3*npts] - the evaluation points (quadrature points) in the reference space
616: Output Parameters:
617: + PetscReal phypts[3*npts] - the evaluation points (quadrature points) transformed to the physical space
618: . PetscReal jxw[npts] - the jacobian determinant * quadrature weight necessary for assembling discrete contributions
619: . PetscReal fe_basis[npts] - the bases values evaluated at the specified quadrature points
620: - 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
622: Level: advanced
624: @*/
625: PetscErrorCode DMMoabFEMComputeBasis(const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscQuadrature quadrature,
626: PetscReal *phypts, PetscReal *jacobian_quadrature_weight_product,
627: PetscReal *fe_basis, PetscReal **fe_basis_derivatives)
628: {
629: PetscInt npoints,idim;
630: bool compute_der;
631: const PetscReal *quadpts, *quadwts;
632: PetscReal jacobian[9], ijacobian[9], volume;
637: compute_der = (fe_basis_derivatives != NULL);
639: /* Get the quadrature points and weights for the given quadrature rule */
640: PetscQuadratureGetData(quadrature, &idim, NULL, &npoints, &quadpts, &quadwts);
642: if (jacobian_quadrature_weight_product) {
643: PetscArraycpy(jacobian_quadrature_weight_product, quadwts, npoints);
644: }
646: switch (dim) {
647: case 1:
648: PetscCall(Compute_Lagrange_Basis_1D_Internal(nverts, coordinates, npoints, quadpts, phypts,
649: jacobian_quadrature_weight_product, fe_basis,
650: (compute_der ? fe_basis_derivatives[0] : NULL),
651: jacobian, ijacobian, &volume));
652: break;
653: case 2:
654: PetscCall(Compute_Lagrange_Basis_2D_Internal(nverts, coordinates, npoints, quadpts, phypts,
655: jacobian_quadrature_weight_product, fe_basis,
656: (compute_der ? fe_basis_derivatives[0] : NULL),
657: (compute_der ? fe_basis_derivatives[1] : NULL),
658: jacobian, ijacobian, &volume));
659: break;
660: case 3:
661: PetscCall(Compute_Lagrange_Basis_3D_Internal(nverts, coordinates, npoints, quadpts, phypts,
662: jacobian_quadrature_weight_product, fe_basis,
663: (compute_der ? fe_basis_derivatives[0] : NULL),
664: (compute_der ? fe_basis_derivatives[1] : NULL),
665: (compute_der ? fe_basis_derivatives[2] : NULL),
666: jacobian, ijacobian, &volume));
667: break;
668: default:
669: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %D", dim);
670: }
671: return 0;
672: }
674: /*@C
675: DMMoabFEMCreateQuadratureDefault - Create default quadrature rules for integration over an element with a given
676: dimension and polynomial order (deciphered from number of element vertices).
678: Input Parameters:
680: + PetscInt dim - the element dimension (1=EDGE, 2=QUAD/TRI, 3=HEX/TET)
681: - PetscInt nverts - the number of vertices in the physical element
683: Output Parameter:
684: . PetscQuadrature quadrature - the quadrature object with default settings to integrate polynomials defined over the element
686: Level: advanced
688: @*/
689: PetscErrorCode DMMoabFEMCreateQuadratureDefault(const PetscInt dim, const PetscInt nverts, PetscQuadrature *quadrature)
690: {
691: PetscReal *w, *x;
692: PetscInt nc=1;
694: /* Create an appropriate quadrature rule to sample basis */
695: switch (dim)
696: {
697: case 1:
698: /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */
699: PetscDTStroudConicalQuadrature(1, nc, nverts, 0, 1.0, quadrature);
700: break;
701: case 2:
702: /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */
703: if (nverts == 3) {
704: const PetscInt order = 2;
705: const PetscInt npoints = (order == 2 ? 3 : 6);
706: PetscMalloc2(npoints * 2, &x, npoints, &w);
707: if (npoints == 3) {
708: x[0] = x[1] = x[2] = x[5] = 1.0 / 6.0;
709: x[3] = x[4] = 2.0 / 3.0;
710: w[0] = w[1] = w[2] = 1.0 / 3.0;
711: } else if (npoints == 6) {
712: x[0] = x[1] = x[2] = 0.44594849091597;
713: x[3] = x[4] = 0.10810301816807;
714: x[5] = 0.44594849091597;
715: x[6] = x[7] = x[8] = 0.09157621350977;
716: x[9] = x[10] = 0.81684757298046;
717: x[11] = 0.09157621350977;
718: w[0] = w[1] = w[2] = 0.22338158967801;
719: w[3] = w[4] = w[5] = 0.10995174365532;
720: } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Triangle quadrature rules for points 3 and 6 supported; npoints : %D", npoints);
721: PetscQuadratureCreate(PETSC_COMM_SELF, quadrature);
722: PetscQuadratureSetOrder(*quadrature, order);
723: PetscQuadratureSetData(*quadrature, dim, nc, npoints, x, w);
724: /* PetscDTStroudConicalQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature); */
725: } else {
726: PetscDTGaussTensorQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature);
727: }
728: break;
729: case 3:
730: /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */
731: if (nverts == 4) {
732: PetscInt order;
733: const PetscInt npoints = 4; // Choose between 4 and 10
734: PetscMalloc2(npoints * 3, &x, npoints, &w);
735: if (npoints == 4) { /* KEAST1, K1, N=4, O=4 */
736: const PetscReal x_4[12] = { 0.5854101966249685, 0.1381966011250105, 0.1381966011250105,
737: 0.1381966011250105, 0.1381966011250105, 0.1381966011250105,
738: 0.1381966011250105, 0.1381966011250105, 0.5854101966249685,
739: 0.1381966011250105, 0.5854101966249685, 0.1381966011250105
740: };
741: PetscArraycpy(x, x_4, 12);
743: w[0] = w[1] = w[2] = w[3] = 1.0 / 24.0;
744: order = 4;
745: } else if (npoints == 10) { /* KEAST3, K3 N=10, O=10 */
746: const PetscReal x_10[30] = { 0.5684305841968444, 0.1438564719343852, 0.1438564719343852,
747: 0.1438564719343852, 0.1438564719343852, 0.1438564719343852,
748: 0.1438564719343852, 0.1438564719343852, 0.5684305841968444,
749: 0.1438564719343852, 0.5684305841968444, 0.1438564719343852,
750: 0.0000000000000000, 0.5000000000000000, 0.5000000000000000,
751: 0.5000000000000000, 0.0000000000000000, 0.5000000000000000,
752: 0.5000000000000000, 0.5000000000000000, 0.0000000000000000,
753: 0.5000000000000000, 0.0000000000000000, 0.0000000000000000,
754: 0.0000000000000000, 0.5000000000000000, 0.0000000000000000,
755: 0.0000000000000000, 0.0000000000000000, 0.5000000000000000
756: };
757: PetscArraycpy(x, x_10, 30);
759: w[0] = w[1] = w[2] = w[3] = 0.2177650698804054;
760: w[4] = w[5] = w[6] = w[7] = w[8] = w[9] = 0.0214899534130631;
761: order = 10;
762: } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Tetrahedral quadrature rules for points 4 and 10 supported; npoints : %D", npoints);
763: PetscQuadratureCreate(PETSC_COMM_SELF, quadrature);
764: PetscQuadratureSetOrder(*quadrature, order);
765: PetscQuadratureSetData(*quadrature, dim, nc, npoints, x, w);
766: /* PetscDTStroudConicalQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature); */
767: } else {
768: PetscDTGaussTensorQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature);
769: }
770: break;
771: default:
772: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %D", dim);
773: }
774: return 0;
775: }
777: /* Compute Jacobians */
778: PetscErrorCode ComputeJacobian_Internal (const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal *quad, PetscReal *phypts,
779: PetscReal *jacobian, PetscReal *ijacobian, PetscReal* dvolume)
780: {
781: PetscInt i;
782: PetscReal volume=1.0;
787: PetscArrayzero(jacobian, dim * dim);
788: if (ijacobian) {
789: PetscArrayzero(ijacobian, dim * dim);
790: }
791: if (phypts) {
792: PetscArrayzero(phypts, /*npts=1 * */ 3);
793: }
795: if (dim == 1) {
796: const PetscReal& r = quad[0];
797: if (nverts == 2) { /* Linear Edge */
798: const PetscReal dNi_dxi[2] = { -1.0, 1.0 };
800: for (i = 0; i < nverts; ++i) {
801: const PetscReal* vertices = coordinates + i * 3;
802: jacobian[0] += dNi_dxi[i] * vertices[0];
803: }
804: } else if (nverts == 3) { /* Quadratic Edge */
805: const PetscReal dNi_dxi[3] = { 4 * r - 3.0, 4 * ( 1.0 - 2.0 * r), 4.0 * r - 1.0};
807: for (i = 0; i < nverts; ++i) {
808: const PetscReal* vertices = coordinates + i * 3;
809: jacobian[0] += dNi_dxi[i] * vertices[0];
810: }
811: } else {
812: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "The number of 1-D entity vertices are invalid. Currently only support EDGE2 and EDGE3 basis evaluations in 1-D : %D", nverts);
813: }
815: if (ijacobian) {
816: /* invert the jacobian */
817: ijacobian[0] = 1.0 / jacobian[0];
818: }
819: } else if (dim == 2) {
821: if (nverts == 4) { /* Linear Quadrangle */
822: const PetscReal& r = quad[0];
823: const PetscReal& s = quad[1];
825: const PetscReal dNi_dxi[4] = { -1.0 + s, 1.0 - s, s, -s };
826: const PetscReal dNi_deta[4] = { -1.0 + r, -r, r, 1.0 - r };
828: for (i = 0; i < nverts; ++i) {
829: const PetscReal* vertices = coordinates + i * 3;
830: jacobian[0] += dNi_dxi[i] * vertices[0];
831: jacobian[2] += dNi_dxi[i] * vertices[1];
832: jacobian[1] += dNi_deta[i] * vertices[0];
833: jacobian[3] += dNi_deta[i] * vertices[1];
834: }
835: } else if (nverts == 3) { /* Linear triangle */
836: const PetscReal x2 = coordinates[2 * 3 + 0], y2 = coordinates[2 * 3 + 1];
838: /* Jacobian is constant */
839: jacobian[0] = (coordinates[0 * 3 + 0] - x2); jacobian[1] = (coordinates[1 * 3 + 0] - x2);
840: jacobian[2] = (coordinates[0 * 3 + 1] - y2); jacobian[3] = (coordinates[1 * 3 + 1] - y2);
841: } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "The number of 2-D entity vertices are invalid. Currently only support QUAD4 and TRI3 basis evaluations in 2-D : %D", nverts);
843: /* invert the jacobian */
844: if (ijacobian) {
845: DMatrix_Invert_2x2_Internal(jacobian, ijacobian, &volume);
846: }
847: } else {
849: if (nverts == 8) { /* Linear Hexahedra */
850: const PetscReal &r = quad[0];
851: const PetscReal &s = quad[1];
852: const PetscReal &t = quad[2];
853: const PetscReal dNi_dxi[8] = {- ( 1.0 - s) * ( 1.0 - t),
854: ( 1.0 - s) * ( 1.0 - t),
855: ( s) * ( 1.0 - t),
856: - ( s) * ( 1.0 - t),
857: - ( 1.0 - s) * ( t),
858: ( 1.0 - s) * ( t),
859: ( s) * ( t),
860: - ( s) * ( t)
861: };
863: const PetscReal dNi_deta[8] = { - ( 1.0 - r) * ( 1.0 - t),
864: - ( r) * ( 1.0 - t),
865: ( r) * ( 1.0 - t),
866: ( 1.0 - r) * ( 1.0 - t),
867: - ( 1.0 - r) * ( t),
868: - ( r) * ( t),
869: ( r) * ( t),
870: ( 1.0 - r) * ( t)
871: };
873: const PetscReal dNi_dzeta[8] = { - ( 1.0 - r) * ( 1.0 - s),
874: - ( r) * ( 1.0 - s),
875: - ( r) * ( s),
876: - ( 1.0 - r) * ( s),
877: ( 1.0 - r) * ( 1.0 - s),
878: ( r) * ( 1.0 - s),
879: ( r) * ( s),
880: ( 1.0 - r) * ( s)
881: };
883: for (i = 0; i < nverts; ++i) {
884: const PetscReal* vertex = coordinates + i * 3;
885: jacobian[0] += dNi_dxi[i] * vertex[0];
886: jacobian[3] += dNi_dxi[i] * vertex[1];
887: jacobian[6] += dNi_dxi[i] * vertex[2];
888: jacobian[1] += dNi_deta[i] * vertex[0];
889: jacobian[4] += dNi_deta[i] * vertex[1];
890: jacobian[7] += dNi_deta[i] * vertex[2];
891: jacobian[2] += dNi_dzeta[i] * vertex[0];
892: jacobian[5] += dNi_dzeta[i] * vertex[1];
893: jacobian[8] += dNi_dzeta[i] * vertex[2];
894: }
895: } else if (nverts == 4) { /* Linear Tetrahedra */
896: const PetscReal x0 = coordinates[/*0 * 3 +*/ 0], y0 = coordinates[/*0 * 3 +*/ 1], z0 = coordinates[/*0 * 3 +*/ 2];
898: /* compute the jacobian */
899: jacobian[0] = coordinates[1 * 3 + 0] - x0; jacobian[1] = coordinates[2 * 3 + 0] - x0; jacobian[2] = coordinates[3 * 3 + 0] - x0;
900: jacobian[3] = coordinates[1 * 3 + 1] - y0; jacobian[4] = coordinates[2 * 3 + 1] - y0; jacobian[5] = coordinates[3 * 3 + 1] - y0;
901: jacobian[6] = coordinates[1 * 3 + 2] - z0; jacobian[7] = coordinates[2 * 3 + 2] - z0; jacobian[8] = coordinates[3 * 3 + 2] - z0;
902: } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "The number of 3-D entity vertices are invalid. Currently only support HEX8 and TET4 basis evaluations in 3-D : %D", nverts);
904: if (ijacobian) {
905: /* invert the jacobian */
906: DMatrix_Invert_3x3_Internal(jacobian, ijacobian, &volume);
907: }
909: }
911: if (dvolume) *dvolume = volume;
912: return 0;
913: }
915: PetscErrorCode FEMComputeBasis_JandF(const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal *quadrature, PetscReal *phypts,
916: PetscReal *phibasis, PetscReal *jacobian, PetscReal *ijacobian, PetscReal* volume)
917: {
918: switch (dim) {
919: case 1:
920: Compute_Lagrange_Basis_1D_Internal(nverts, coordinates, 1, quadrature, phypts, NULL, phibasis, NULL, jacobian, ijacobian, volume);
921: break;
922: case 2:
923: Compute_Lagrange_Basis_2D_Internal(nverts, coordinates, 1, quadrature, phypts, NULL, phibasis, NULL, NULL, jacobian, ijacobian, volume);
924: break;
925: case 3:
926: Compute_Lagrange_Basis_3D_Internal(nverts, coordinates, 1, quadrature, phypts, NULL, phibasis, NULL, NULL, NULL, jacobian, ijacobian, volume);
927: break;
928: default:
929: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %D", dim);
930: }
931: return 0;
932: }
934: /*@C
935: DMMoabPToRMapping - Compute the mapping from the physical coordinate system for a given element to the
936: canonical reference element. In addition to finding the inverse mapping evaluation through Newton iteration,
937: the basis function at the parametric point is also evaluated optionally.
939: Input Parameters:
940: + PetscInt dim - the element dimension (1=EDGE, 2=QUAD/TRI, 3=HEX/TET)
941: . PetscInt nverts - the number of vertices in the physical element
942: . PetscReal coordinates - the coordinates of vertices in the physical element
943: - PetscReal[3] xphy - the coordinates of physical point for which natural coordinates (in reference frame) are sought
945: Output Parameters:
946: + PetscReal[3] natparam - the natural coordinates (in reference frame) corresponding to xphy
947: - PetscReal[nverts] phi - the basis functions evaluated at the natural coordinates (natparam)
949: Level: advanced
951: @*/
952: PetscErrorCode DMMoabPToRMapping(const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal* xphy, PetscReal* natparam, PetscReal* phi)
953: {
954: /* Perform inverse evaluation for the mapping with use of Newton Raphson iteration */
955: const PetscReal tol = 1.0e-10;
956: const PetscInt max_iterations = 10;
957: const PetscReal error_tol_sqr = tol*tol;
958: PetscReal phibasis[8], jacobian[9], ijacobian[9], volume;
959: PetscReal phypts[3] = {0.0, 0.0, 0.0};
960: PetscReal delta[3] = {0.0, 0.0, 0.0};
961: PetscInt iters=0;
962: PetscReal error=1.0;
968: PetscArrayzero(jacobian, dim * dim);
969: PetscArrayzero(ijacobian, dim * dim);
970: PetscArrayzero(phibasis, nverts);
972: /* zero initial guess */
973: natparam[0] = natparam[1] = natparam[2] = 0.0;
975: /* Compute delta = evaluate( xi) - x */
976: FEMComputeBasis_JandF(dim, nverts, coordinates, natparam, phypts, phibasis, jacobian, ijacobian, &volume);
978: error = 0.0;
979: switch(dim) {
980: case 3:
981: delta[2] = phypts[2] - xphy[2];
982: error += (delta[2]*delta[2]);
983: case 2:
984: delta[1] = phypts[1] - xphy[1];
985: error += (delta[1]*delta[1]);
986: case 1:
987: delta[0] = phypts[0] - xphy[0];
988: error += (delta[0]*delta[0]);
989: break;
990: }
992: while (error > error_tol_sqr) {
994: if (++iters > max_iterations)
995: SETERRQ(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]);
997: /* Compute natparam -= J.inverse() * delta */
998: switch(dim) {
999: case 1:
1000: natparam[0] -= ijacobian[0] * delta[0];
1001: break;
1002: case 2:
1003: natparam[0] -= ijacobian[0] * delta[0] + ijacobian[1] * delta[1];
1004: natparam[1] -= ijacobian[2] * delta[0] + ijacobian[3] * delta[1];
1005: break;
1006: case 3:
1007: natparam[0] -= ijacobian[0] * delta[0] + ijacobian[1] * delta[1] + ijacobian[2] * delta[2];
1008: natparam[1] -= ijacobian[3] * delta[0] + ijacobian[4] * delta[1] + ijacobian[5] * delta[2];
1009: natparam[2] -= ijacobian[6] * delta[0] + ijacobian[7] * delta[1] + ijacobian[8] * delta[2];
1010: break;
1011: }
1013: /* Compute delta = evaluate( xi) - x */
1014: FEMComputeBasis_JandF(dim, nverts, coordinates, natparam, phypts, phibasis, jacobian, ijacobian, &volume);
1016: error = 0.0;
1017: switch(dim) {
1018: case 3:
1019: delta[2] = phypts[2] - xphy[2];
1020: error += (delta[2]*delta[2]);
1021: case 2:
1022: delta[1] = phypts[1] - xphy[1];
1023: error += (delta[1]*delta[1]);
1024: case 1:
1025: delta[0] = phypts[0] - xphy[0];
1026: error += (delta[0]*delta[0]);
1027: break;
1028: }
1029: }
1030: if (phi) {
1031: PetscArraycpy(phi, phibasis, nverts);
1032: }
1033: return 0;
1034: }