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]) - inmat[3] * (inmat[8] * inmat[1] - inmat[7] * inmat[2]) + inmat[6] * (inmat[5] * inmat[1] - inmat[4] * inmat[2]);
28: }
30: static inline PetscErrorCode DMatrix_Invert_3x3_Internal(const PetscReal *inmat, PetscReal *outmat, PetscScalar *determinant)
31: {
32: PetscReal det = DMatrix_Determinant_3x3_Internal(inmat);
33: if (outmat) {
34: outmat[0] = (inmat[8] * inmat[4] - inmat[7] * inmat[5]) / det;
35: outmat[1] = -(inmat[8] * inmat[1] - inmat[7] * inmat[2]) / det;
36: outmat[2] = (inmat[5] * inmat[1] - inmat[4] * inmat[2]) / det;
37: outmat[3] = -(inmat[8] * inmat[3] - inmat[6] * inmat[5]) / det;
38: outmat[4] = (inmat[8] * inmat[0] - inmat[6] * inmat[2]) / det;
39: outmat[5] = -(inmat[5] * inmat[0] - inmat[3] * inmat[2]) / det;
40: outmat[6] = (inmat[7] * inmat[3] - inmat[6] * inmat[4]) / det;
41: outmat[7] = -(inmat[7] * inmat[0] - inmat[6] * inmat[1]) / det;
42: outmat[8] = (inmat[4] * inmat[0] - inmat[3] * inmat[1]) / det;
43: }
44: if (determinant) *determinant = det;
45: return 0;
46: }
48: inline PetscReal DMatrix_Determinant_4x4_Internal(PetscReal inmat[4 * 4])
49: {
50: return inmat[0 + 0 * 4] * (inmat[1 + 1 * 4] * (inmat[2 + 2 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 2 * 4]) - inmat[1 + 2 * 4] * (inmat[2 + 1 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 1 * 4]) + inmat[1 + 3 * 4] * (inmat[2 + 1 * 4] * inmat[3 + 2 * 4] - inmat[2 + 2 * 4] * inmat[3 + 1 * 4])) - inmat[0 + 1 * 4] * (inmat[1 + 0 * 4] * (inmat[2 + 2 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 2 * 4]) - inmat[1 + 2 * 4] * (inmat[2 + 0 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 0 * 4]) + inmat[1 + 3 * 4] * (inmat[2 + 0 * 4] * inmat[3 + 2 * 4] - inmat[2 + 2 * 4] * inmat[3 + 0 * 4])) + inmat[0 + 2 * 4] * (inmat[1 + 0 * 4] * (inmat[2 + 1 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 1 * 4]) - inmat[1 + 1 * 4] * (inmat[2 + 0 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 0 * 4]) + inmat[1 + 3 * 4] * (inmat[2 + 0 * 4] * inmat[3 + 1 * 4] - inmat[2 + 1 * 4] * inmat[3 + 0 * 4])) - inmat[0 + 3 * 4] * (inmat[1 + 0 * 4] * (inmat[2 + 1 * 4] * inmat[3 + 2 * 4] - inmat[2 + 2 * 4] * inmat[3 + 1 * 4]) - inmat[1 + 1 * 4] * (inmat[2 + 0 * 4] * inmat[3 + 2 * 4] - inmat[2 + 2 * 4] * inmat[3 + 0 * 4]) + inmat[1 + 2 * 4] * (inmat[2 + 0 * 4] * inmat[3 + 1 * 4] - inmat[2 + 1 * 4] * inmat[3 + 0 * 4]));
51: }
53: inline PetscErrorCode DMatrix_Invert_4x4_Internal(PetscReal *inmat, PetscReal *outmat, PetscScalar *determinant)
54: {
55: PetscReal det = DMatrix_Determinant_4x4_Internal(inmat);
56: if (outmat) {
57: outmat[0] = (inmat[5] * inmat[10] * inmat[15] + inmat[6] * inmat[11] * inmat[13] + inmat[7] * inmat[9] * inmat[14] - inmat[5] * inmat[11] * inmat[14] - inmat[6] * inmat[9] * inmat[15] - inmat[7] * inmat[10] * inmat[13]) / det;
58: outmat[1] = (inmat[1] * inmat[11] * inmat[14] + inmat[2] * inmat[9] * inmat[15] + inmat[3] * inmat[10] * inmat[13] - inmat[1] * inmat[10] * inmat[15] - inmat[2] * inmat[11] * inmat[13] - inmat[3] * inmat[9] * inmat[14]) / det;
59: outmat[2] = (inmat[1] * inmat[6] * inmat[15] + inmat[2] * inmat[7] * inmat[13] + inmat[3] * inmat[5] * inmat[14] - inmat[1] * inmat[7] * inmat[14] - inmat[2] * inmat[5] * inmat[15] - inmat[3] * inmat[6] * inmat[13]) / det;
60: outmat[3] = (inmat[1] * inmat[7] * inmat[10] + inmat[2] * inmat[5] * inmat[11] + inmat[3] * inmat[6] * inmat[9] - inmat[1] * inmat[6] * inmat[11] - inmat[2] * inmat[7] * inmat[9] - inmat[3] * inmat[5] * inmat[10]) / det;
61: outmat[4] = (inmat[4] * inmat[11] * inmat[14] + inmat[6] * inmat[8] * inmat[15] + inmat[7] * inmat[10] * inmat[12] - inmat[4] * inmat[10] * inmat[15] - inmat[6] * inmat[11] * inmat[12] - inmat[7] * inmat[8] * inmat[14]) / det;
62: outmat[5] = (inmat[0] * inmat[10] * inmat[15] + inmat[2] * inmat[11] * inmat[12] + inmat[3] * inmat[8] * inmat[14] - inmat[0] * inmat[11] * inmat[14] - inmat[2] * inmat[8] * inmat[15] - inmat[3] * inmat[10] * inmat[12]) / det;
63: outmat[6] = (inmat[0] * inmat[7] * inmat[14] + inmat[2] * inmat[4] * inmat[15] + inmat[3] * inmat[6] * inmat[12] - inmat[0] * inmat[6] * inmat[15] - inmat[2] * inmat[7] * inmat[12] - inmat[3] * inmat[4] * inmat[14]) / det;
64: outmat[7] = (inmat[0] * inmat[6] * inmat[11] + inmat[2] * inmat[7] * inmat[8] + inmat[3] * inmat[4] * inmat[10] - inmat[0] * inmat[7] * inmat[10] - inmat[2] * inmat[4] * inmat[11] - inmat[3] * inmat[6] * inmat[8]) / det;
65: outmat[8] = (inmat[4] * inmat[9] * inmat[15] + inmat[5] * inmat[11] * inmat[12] + inmat[7] * inmat[8] * inmat[13] - inmat[4] * inmat[11] * inmat[13] - inmat[5] * inmat[8] * inmat[15] - inmat[7] * inmat[9] * inmat[12]) / det;
66: outmat[9] = (inmat[0] * inmat[11] * inmat[13] + inmat[1] * inmat[8] * inmat[15] + inmat[3] * inmat[9] * inmat[12] - inmat[0] * inmat[9] * inmat[15] - inmat[1] * inmat[11] * inmat[12] - inmat[3] * inmat[8] * inmat[13]) / det;
67: outmat[10] = (inmat[0] * inmat[5] * inmat[15] + inmat[1] * inmat[7] * inmat[12] + inmat[3] * inmat[4] * inmat[13] - inmat[0] * inmat[7] * inmat[13] - inmat[1] * inmat[4] * inmat[15] - inmat[3] * inmat[5] * inmat[12]) / det;
68: outmat[11] = (inmat[0] * inmat[7] * inmat[9] + inmat[1] * inmat[4] * inmat[11] + inmat[3] * inmat[5] * inmat[8] - inmat[0] * inmat[5] * inmat[11] - inmat[1] * inmat[7] * inmat[8] - inmat[3] * inmat[4] * inmat[9]) / det;
69: outmat[12] = (inmat[4] * inmat[10] * inmat[13] + inmat[5] * inmat[8] * inmat[14] + inmat[6] * inmat[9] * inmat[12] - inmat[4] * inmat[9] * inmat[14] - inmat[5] * inmat[10] * inmat[12] - inmat[6] * inmat[8] * inmat[13]) / det;
70: outmat[13] = (inmat[0] * inmat[9] * inmat[14] + inmat[1] * inmat[10] * inmat[12] + inmat[2] * inmat[8] * inmat[13] - inmat[0] * inmat[10] * inmat[13] - inmat[1] * inmat[8] * inmat[14] - inmat[2] * inmat[9] * inmat[12]) / det;
71: outmat[14] = (inmat[0] * inmat[6] * inmat[13] + inmat[1] * inmat[4] * inmat[14] + inmat[2] * inmat[5] * inmat[12] - inmat[0] * inmat[5] * inmat[14] - inmat[1] * inmat[6] * inmat[12] - inmat[2] * inmat[4] * inmat[13]) / det;
72: outmat[15] = (inmat[0] * inmat[5] * inmat[10] + inmat[1] * inmat[6] * inmat[8] + inmat[2] * inmat[4] * inmat[9] - inmat[0] * inmat[6] * inmat[9] - inmat[1] * inmat[4] * inmat[10] - inmat[2] * inmat[5] * inmat[8]) / det;
73: }
74: if (determinant) *determinant = det;
75: return 0;
76: }
78: /*@C
79: Compute_Lagrange_Basis_1D_Internal - Evaluate bases and derivatives at quadrature points for a EDGE2 or EDGE3 element.
81: The routine is given the coordinates of the vertices of a linear or quadratic edge element.
83: The routine evaluates the basis functions associated with each quadrature point provided,
84: and their derivatives with respect to X.
86: Notes:
88: Example Physical Element
89: .vb
90: 1-------2 1----3----2
91: EDGE2 EDGE3
92: .ve
94: Input Parameters:
95: + PetscInt nverts - the number of element vertices
96: . PetscReal coords[3*nverts] - the physical coordinates of the vertices (in canonical numbering)
97: . PetscInt npts - the number of evaluation points (quadrature points)
98: - PetscReal quad[3*npts] - the evaluation points (quadrature points) in the reference space
100: Output Parameters:
101: + PetscReal phypts[3*npts] - the evaluation points (quadrature points) transformed to the physical space
102: . PetscReal jxw[npts] - the jacobian determinant * quadrature weight necessary for assembling discrete contributions
103: . PetscReal phi[npts] - the bases evaluated at the specified quadrature points
104: . PetscReal dphidx[npts] - the derivative of the bases wrt X-direction evaluated at the specified quadrature points
105: . jacobian - jacobian
106: . ijacobian - ijacobian
107: - volume - volume
109: Level: advanced
111: @*/
112: PetscErrorCode Compute_Lagrange_Basis_1D_Internal(const PetscInt nverts, const PetscReal *coords, const PetscInt npts, const PetscReal *quad, PetscReal *phypts, PetscReal *jxw, PetscReal *phi, PetscReal *dphidx, PetscReal *jacobian, PetscReal *ijacobian, PetscReal *volume)
113: {
114: int i, j;
119: if (phypts) PetscArrayzero(phypts, npts * 3);
120: if (dphidx) { /* Reset arrays. */
121: PetscArrayzero(dphidx, npts * nverts);
122: }
123: if (nverts == 2) { /* Linear Edge */
125: for (j = 0; j < npts; j++) {
126: const PetscInt offset = j * nverts;
127: const PetscReal r = quad[j];
129: phi[0 + offset] = (1.0 - r);
130: phi[1 + offset] = (r);
132: const PetscReal dNi_dxi[2] = {-1.0, 1.0};
134: jacobian[0] = ijacobian[0] = volume[0] = 0.0;
135: for (i = 0; i < nverts; ++i) {
136: const PetscReal *vertices = coords + i * 3;
137: jacobian[0] += dNi_dxi[i] * vertices[0];
138: if (phypts) phypts[3 * j + 0] += phi[i + offset] * vertices[0];
139: }
141: /* invert the jacobian */
142: *volume = jacobian[0];
143: ijacobian[0] = 1.0 / jacobian[0];
144: jxw[j] *= *volume;
146: /* Divide by element jacobian. */
147: for (i = 0; i < nverts; i++) {
148: if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[0];
149: }
150: }
151: } else if (nverts == 3) { /* Quadratic Edge */
153: for (j = 0; j < npts; j++) {
154: const PetscInt offset = j * nverts;
155: const PetscReal r = quad[j];
157: phi[0 + offset] = 1.0 + r * (2.0 * r - 3.0);
158: phi[1 + offset] = 4.0 * r * (1.0 - r);
159: phi[2 + offset] = r * (2.0 * r - 1.0);
161: const PetscReal dNi_dxi[3] = {4 * r - 3.0, 4 * (1.0 - 2.0 * r), 4.0 * r - 1.0};
163: jacobian[0] = ijacobian[0] = volume[0] = 0.0;
164: for (i = 0; i < nverts; ++i) {
165: const PetscReal *vertices = coords + i * 3;
166: jacobian[0] += dNi_dxi[i] * vertices[0];
167: if (phypts) phypts[3 * j + 0] += phi[i + offset] * vertices[0];
168: }
170: /* invert the jacobian */
171: *volume = jacobian[0];
172: ijacobian[0] = 1.0 / jacobian[0];
173: if (jxw) jxw[j] *= *volume;
175: /* Divide by element jacobian. */
176: for (i = 0; i < nverts; i++) {
177: if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[0];
178: }
179: }
180: } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "The number of entity vertices are invalid. Currently only support EDGE2 and EDGE3 basis evaluations in 1-D : %" PetscInt_FMT, nverts);
181: return 0;
182: }
184: /*@C
185: Compute_Lagrange_Basis_2D_Internal - Evaluate bases and derivatives at quadrature points for a QUAD4 or TRI3 element.
187: The routine is given the coordinates of the vertices of a quadrangle or triangle.
189: The routine evaluates the basis functions associated with each quadrature point provided,
190: and their derivatives with respect to X and Y.
192: Notes:
194: Example Physical Element (QUAD4)
195: .vb
196: 4------3 s
197: | | |
198: | | |
199: | | |
200: 1------2 0-------r
201: .ve
203: Input Parameters:
204: + PetscInt nverts - the number of element vertices
205: . PetscReal coords[3*nverts] - the physical coordinates of the vertices (in canonical numbering)
206: . PetscInt npts - the number of evaluation points (quadrature points)
207: - PetscReal quad[3*npts] - the evaluation points (quadrature points) in the reference space
209: Output Parameters:
210: + PetscReal phypts[3*npts] - the evaluation points (quadrature points) transformed to the physical space
211: . PetscReal jxw[npts] - the jacobian determinant * quadrature weight necessary for assembling discrete contributions
212: . PetscReal phi[npts] - the bases evaluated at the specified quadrature points
213: . PetscReal dphidx[npts] - the derivative of the bases wrt X-direction evaluated at the specified quadrature points
214: . PetscReal dphidy[npts] - the derivative of the bases wrt Y-direction evaluated at the specified quadrature points
215: . jacobian - jacobian
216: . ijacobian - ijacobian
217: - volume - volume
219: Level: advanced
221: @*/
222: PetscErrorCode Compute_Lagrange_Basis_2D_Internal(const PetscInt nverts, const PetscReal *coords, const PetscInt npts, const PetscReal *quad, PetscReal *phypts, PetscReal *jxw, PetscReal *phi, PetscReal *dphidx, PetscReal *dphidy, PetscReal *jacobian, PetscReal *ijacobian, PetscReal *volume)
223: {
224: PetscInt i, j, k;
229: PetscArrayzero(phi, npts);
230: if (phypts) PetscArrayzero(phypts, npts * 3);
231: if (dphidx) { /* Reset arrays. */
232: PetscArrayzero(dphidx, npts * nverts);
233: PetscArrayzero(dphidy, npts * nverts);
234: }
235: if (nverts == 4) { /* Linear Quadrangle */
237: for (j = 0; j < npts; j++) {
238: const PetscInt offset = j * nverts;
239: const PetscReal r = quad[0 + j * 2];
240: const PetscReal s = quad[1 + j * 2];
242: phi[0 + offset] = (1.0 - r) * (1.0 - s);
243: phi[1 + offset] = r * (1.0 - s);
244: phi[2 + offset] = r * s;
245: phi[3 + offset] = (1.0 - r) * s;
247: const PetscReal dNi_dxi[4] = {-1.0 + s, 1.0 - s, s, -s};
248: const PetscReal dNi_deta[4] = {-1.0 + r, -r, r, 1.0 - r};
250: PetscArrayzero(jacobian, 4);
251: PetscArrayzero(ijacobian, 4);
252: for (i = 0; i < nverts; ++i) {
253: const PetscReal *vertices = coords + i * 3;
254: jacobian[0] += dNi_dxi[i] * vertices[0];
255: jacobian[2] += dNi_dxi[i] * vertices[1];
256: jacobian[1] += dNi_deta[i] * vertices[0];
257: jacobian[3] += dNi_deta[i] * vertices[1];
258: if (phypts) {
259: phypts[3 * j + 0] += phi[i + offset] * vertices[0];
260: phypts[3 * j + 1] += phi[i + offset] * vertices[1];
261: phypts[3 * j + 2] += phi[i + offset] * vertices[2];
262: }
263: }
265: /* invert the jacobian */
266: DMatrix_Invert_2x2_Internal(jacobian, ijacobian, volume);
269: if (jxw) jxw[j] *= *volume;
271: /* Let us compute the bases derivatives by scaling with inverse jacobian. */
272: if (dphidx) {
273: for (i = 0; i < nverts; i++) {
274: for (k = 0; k < 2; ++k) {
275: if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[k * 2 + 0];
276: if (dphidy) dphidy[i + offset] += dNi_deta[i] * ijacobian[k * 2 + 1];
277: } /* for k=[0..2) */
278: } /* for i=[0..nverts) */
279: } /* if (dphidx) */
280: }
281: } else if (nverts == 3) { /* Linear triangle */
282: const PetscReal x2 = coords[2 * 3 + 0], y2 = coords[2 * 3 + 1];
284: PetscArrayzero(jacobian, 4);
285: PetscArrayzero(ijacobian, 4);
287: /* Jacobian is constant */
288: jacobian[0] = (coords[0 * 3 + 0] - x2);
289: jacobian[1] = (coords[1 * 3 + 0] - x2);
290: jacobian[2] = (coords[0 * 3 + 1] - y2);
291: jacobian[3] = (coords[1 * 3 + 1] - y2);
293: /* invert the jacobian */
294: DMatrix_Invert_2x2_Internal(jacobian, ijacobian, volume);
297: const PetscReal Dx[3] = {ijacobian[0], ijacobian[2], -ijacobian[0] - ijacobian[2]};
298: const PetscReal Dy[3] = {ijacobian[1], ijacobian[3], -ijacobian[1] - ijacobian[3]};
300: for (j = 0; j < npts; j++) {
301: const PetscInt offset = j * nverts;
302: const PetscReal r = quad[0 + j * 2];
303: const PetscReal s = quad[1 + j * 2];
305: if (jxw) jxw[j] *= 0.5;
307: /* const PetscReal Ni[3] = { r, s, 1.0 - r - s }; */
308: const PetscReal phipts_x = coords[2 * 3 + 0] + jacobian[0] * r + jacobian[1] * s;
309: const PetscReal phipts_y = coords[2 * 3 + 1] + jacobian[2] * r + jacobian[3] * s;
310: if (phypts) {
311: phypts[offset + 0] = phipts_x;
312: phypts[offset + 1] = phipts_y;
313: }
315: /* \phi_0 = (b.y - c.y) x + (b.x - c.x) y + c.x b.y - b.x c.y */
316: phi[0 + offset] = (ijacobian[0] * (phipts_x - x2) + ijacobian[1] * (phipts_y - y2));
317: /* \phi_1 = (c.y - a.y) x + (a.x - c.x) y + c.x a.y - a.x c.y */
318: phi[1 + offset] = (ijacobian[2] * (phipts_x - x2) + ijacobian[3] * (phipts_y - y2));
319: phi[2 + offset] = 1.0 - phi[0 + offset] - phi[1 + offset];
321: if (dphidx) {
322: dphidx[0 + offset] = Dx[0];
323: dphidx[1 + offset] = Dx[1];
324: dphidx[2 + offset] = Dx[2];
325: }
327: if (dphidy) {
328: dphidy[0 + offset] = Dy[0];
329: dphidy[1 + offset] = Dy[1];
330: dphidy[2 + offset] = Dy[2];
331: }
332: }
333: } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "The number of entity vertices are invalid. Currently only support QUAD4 and TRI3 basis evaluations in 2-D : %" PetscInt_FMT, nverts);
334: return 0;
335: }
337: /*@C
338: Compute_Lagrange_Basis_3D_Internal - Evaluate bases and derivatives at quadrature points for a HEX8 or TET4 element.
340: The routine is given the coordinates of the vertices of a hexahedra or tetrahedra.
342: The routine evaluates the basis functions associated with each quadrature point provided,
343: and their derivatives with respect to X, Y and Z.
345: Notes:
347: Example Physical Element (HEX8)
348: .vb
349: 8------7
350: /| /| t s
351: 5------6 | | /
352: | | | | |/
353: | 4----|-3 0-------r
354: |/ |/
355: 1------2
356: .ve
358: Input Parameters:
359: + PetscInt nverts - the number of element vertices
360: . PetscReal coords[3*nverts] - the physical coordinates of the vertices (in canonical numbering)
361: . PetscInt npts - the number of evaluation points (quadrature points)
362: - PetscReal quad[3*npts] - the evaluation points (quadrature points) in the reference space
364: Output Parameters:
365: + PetscReal phypts[3*npts] - the evaluation points (quadrature points) transformed to the physical space
366: . PetscReal jxw[npts] - the jacobian determinant * quadrature weight necessary for assembling discrete contributions
367: . PetscReal phi[npts] - the bases evaluated at the specified quadrature points
368: . PetscReal dphidx[npts] - the derivative of the bases wrt X-direction evaluated at the specified quadrature points
369: . PetscReal dphidy[npts] - the derivative of the bases wrt Y-direction evaluated at the specified quadrature points
370: . PetscReal dphidz[npts] - the derivative of the bases wrt Z-direction evaluated at the specified quadrature points
371: . jacobian - jacobian
372: . ijacobian - ijacobian
373: - volume - volume
375: Level: advanced
377: @*/
378: PetscErrorCode Compute_Lagrange_Basis_3D_Internal(const PetscInt nverts, const PetscReal *coords, const PetscInt npts, const PetscReal *quad, PetscReal *phypts, PetscReal *jxw, PetscReal *phi, PetscReal *dphidx, PetscReal *dphidy, PetscReal *dphidz, PetscReal *jacobian, PetscReal *ijacobian, PetscReal *volume)
379: {
380: PetscInt i, j, k;
386: PetscArrayzero(phi, npts);
387: if (phypts) PetscArrayzero(phypts, npts * 3);
388: if (dphidx) {
389: PetscArrayzero(dphidx, npts * nverts);
390: PetscArrayzero(dphidy, npts * nverts);
391: PetscArrayzero(dphidz, npts * nverts);
392: }
394: if (nverts == 8) { /* Linear Hexahedra */
396: for (j = 0; j < npts; j++) {
397: const PetscInt offset = j * nverts;
398: const PetscReal &r = quad[j * 3 + 0];
399: const PetscReal &s = quad[j * 3 + 1];
400: const PetscReal &t = quad[j * 3 + 2];
402: phi[offset + 0] = (1.0 - r) * (1.0 - s) * (1.0 - t); /* 0,0,0 */
403: phi[offset + 1] = (r) * (1.0 - s) * (1.0 - t); /* 1,0,0 */
404: phi[offset + 2] = (r) * (s) * (1.0 - t); /* 1,1,0 */
405: phi[offset + 3] = (1.0 - r) * (s) * (1.0 - t); /* 0,1,0 */
406: phi[offset + 4] = (1.0 - r) * (1.0 - s) * (t); /* 0,0,1 */
407: phi[offset + 5] = (r) * (1.0 - s) * (t); /* 1,0,1 */
408: phi[offset + 6] = (r) * (s) * (t); /* 1,1,1 */
409: phi[offset + 7] = (1.0 - r) * (s) * (t); /* 0,1,1 */
411: const PetscReal dNi_dxi[8] = {-(1.0 - s) * (1.0 - t), (1.0 - s) * (1.0 - t), (s) * (1.0 - t), -(s) * (1.0 - t), -(1.0 - s) * (t), (1.0 - s) * (t), (s) * (t), -(s) * (t)};
413: const PetscReal dNi_deta[8] = {-(1.0 - r) * (1.0 - t), -(r) * (1.0 - t), (r) * (1.0 - t), (1.0 - r) * (1.0 - t), -(1.0 - r) * (t), -(r) * (t), (r) * (t), (1.0 - r) * (t)};
415: const PetscReal dNi_dzeta[8] = {-(1.0 - r) * (1.0 - s), -(r) * (1.0 - s), -(r) * (s), -(1.0 - r) * (s), (1.0 - r) * (1.0 - s), (r) * (1.0 - s), (r) * (s), (1.0 - r) * (s)};
417: PetscArrayzero(jacobian, 9);
418: PetscArrayzero(ijacobian, 9);
419: for (i = 0; i < nverts; ++i) {
420: const PetscReal *vertex = coords + i * 3;
421: jacobian[0] += dNi_dxi[i] * vertex[0];
422: jacobian[3] += dNi_dxi[i] * vertex[1];
423: jacobian[6] += dNi_dxi[i] * vertex[2];
424: jacobian[1] += dNi_deta[i] * vertex[0];
425: jacobian[4] += dNi_deta[i] * vertex[1];
426: jacobian[7] += dNi_deta[i] * vertex[2];
427: jacobian[2] += dNi_dzeta[i] * vertex[0];
428: jacobian[5] += dNi_dzeta[i] * vertex[1];
429: jacobian[8] += dNi_dzeta[i] * vertex[2];
430: if (phypts) {
431: phypts[3 * j + 0] += phi[i + offset] * vertex[0];
432: phypts[3 * j + 1] += phi[i + offset] * vertex[1];
433: phypts[3 * j + 2] += phi[i + offset] * vertex[2];
434: }
435: }
437: /* invert the jacobian */
438: DMatrix_Invert_3x3_Internal(jacobian, ijacobian, volume);
441: if (jxw) jxw[j] *= (*volume);
443: /* Divide by element jacobian. */
444: for (i = 0; i < nverts; ++i) {
445: for (k = 0; k < 3; ++k) {
446: if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[0 * 3 + k];
447: if (dphidy) dphidy[i + offset] += dNi_deta[i] * ijacobian[1 * 3 + k];
448: if (dphidz) dphidz[i + offset] += dNi_dzeta[i] * ijacobian[2 * 3 + k];
449: }
450: }
451: }
452: } else if (nverts == 4) { /* Linear Tetrahedra */
453: PetscReal Dx[4] = {0, 0, 0, 0}, Dy[4] = {0, 0, 0, 0}, Dz[4] = {0, 0, 0, 0};
454: const PetscReal x0 = coords[/*0 * 3 +*/ 0], y0 = coords[/*0 * 3 +*/ 1], z0 = coords[/*0 * 3 +*/ 2];
456: PetscArrayzero(jacobian, 9);
457: PetscArrayzero(ijacobian, 9);
459: /* compute the jacobian */
460: jacobian[0] = coords[1 * 3 + 0] - x0;
461: jacobian[1] = coords[2 * 3 + 0] - x0;
462: jacobian[2] = coords[3 * 3 + 0] - x0;
463: jacobian[3] = coords[1 * 3 + 1] - y0;
464: jacobian[4] = coords[2 * 3 + 1] - y0;
465: jacobian[5] = coords[3 * 3 + 1] - y0;
466: jacobian[6] = coords[1 * 3 + 2] - z0;
467: jacobian[7] = coords[2 * 3 + 2] - z0;
468: jacobian[8] = coords[3 * 3 + 2] - z0;
470: /* invert the jacobian */
471: DMatrix_Invert_3x3_Internal(jacobian, ijacobian, volume);
474: if (dphidx) {
475: Dx[0] = (coords[1 + 2 * 3] * (coords[2 + 1 * 3] - coords[2 + 3 * 3]) - coords[1 + 1 * 3] * (coords[2 + 2 * 3] - coords[2 + 3 * 3]) - coords[1 + 3 * 3] * (coords[2 + 1 * 3] - coords[2 + 2 * 3])) / *volume;
476: Dx[1] = -(coords[1 + 2 * 3] * (coords[2 + 0 * 3] - coords[2 + 3 * 3]) - coords[1 + 0 * 3] * (coords[2 + 2 * 3] - coords[2 + 3 * 3]) - coords[1 + 3 * 3] * (coords[2 + 0 * 3] - coords[2 + 2 * 3])) / *volume;
477: Dx[2] = (coords[1 + 1 * 3] * (coords[2 + 0 * 3] - coords[2 + 3 * 3]) - coords[1 + 0 * 3] * (coords[2 + 1 * 3] - coords[2 + 3 * 3]) - coords[1 + 3 * 3] * (coords[2 + 0 * 3] - coords[2 + 1 * 3])) / *volume;
478: Dx[3] = -(Dx[0] + Dx[1] + Dx[2]);
479: Dy[0] = (coords[0 + 1 * 3] * (coords[2 + 2 * 3] - coords[2 + 3 * 3]) - coords[0 + 2 * 3] * (coords[2 + 1 * 3] - coords[2 + 3 * 3]) + coords[0 + 3 * 3] * (coords[2 + 1 * 3] - coords[2 + 2 * 3])) / *volume;
480: Dy[1] = -(coords[0 + 0 * 3] * (coords[2 + 2 * 3] - coords[2 + 3 * 3]) - coords[0 + 2 * 3] * (coords[2 + 0 * 3] - coords[2 + 3 * 3]) + coords[0 + 3 * 3] * (coords[2 + 0 * 3] - coords[2 + 2 * 3])) / *volume;
481: Dy[2] = (coords[0 + 0 * 3] * (coords[2 + 1 * 3] - coords[2 + 3 * 3]) - coords[0 + 1 * 3] * (coords[2 + 0 * 3] - coords[2 + 3 * 3]) + coords[0 + 3 * 3] * (coords[2 + 0 * 3] - coords[2 + 1 * 3])) / *volume;
482: Dy[3] = -(Dy[0] + Dy[1] + Dy[2]);
483: Dz[0] = (coords[0 + 1 * 3] * (coords[1 + 3 * 3] - coords[1 + 2 * 3]) - coords[0 + 2 * 3] * (coords[1 + 3 * 3] - coords[1 + 1 * 3]) + coords[0 + 3 * 3] * (coords[1 + 2 * 3] - coords[1 + 1 * 3])) / *volume;
484: Dz[1] = -(coords[0 + 0 * 3] * (coords[1 + 3 * 3] - coords[1 + 2 * 3]) + coords[0 + 2 * 3] * (coords[1 + 0 * 3] - coords[1 + 3 * 3]) - coords[0 + 3 * 3] * (coords[1 + 0 * 3] - coords[1 + 2 * 3])) / *volume;
485: Dz[2] = (coords[0 + 0 * 3] * (coords[1 + 3 * 3] - coords[1 + 1 * 3]) + coords[0 + 1 * 3] * (coords[1 + 0 * 3] - coords[1 + 3 * 3]) - coords[0 + 3 * 3] * (coords[1 + 0 * 3] - coords[1 + 1 * 3])) / *volume;
486: Dz[3] = -(Dz[0] + Dz[1] + Dz[2]);
487: }
489: for (j = 0; j < npts; j++) {
490: const PetscInt offset = j * nverts;
491: const PetscReal &r = quad[j * 3 + 0];
492: const PetscReal &s = quad[j * 3 + 1];
493: const PetscReal &t = quad[j * 3 + 2];
495: if (jxw) jxw[j] *= *volume;
497: phi[offset + 0] = 1.0 - r - s - t;
498: phi[offset + 1] = r;
499: phi[offset + 2] = s;
500: phi[offset + 3] = t;
502: if (phypts) {
503: for (i = 0; i < nverts; ++i) {
504: const PetscScalar *vertices = coords + i * 3;
505: phypts[3 * j + 0] += phi[i + offset] * vertices[0];
506: phypts[3 * j + 1] += phi[i + offset] * vertices[1];
507: phypts[3 * j + 2] += phi[i + offset] * vertices[2];
508: }
509: }
511: /* Now set the derivatives */
512: if (dphidx) {
513: dphidx[0 + offset] = Dx[0];
514: dphidx[1 + offset] = Dx[1];
515: dphidx[2 + offset] = Dx[2];
516: dphidx[3 + offset] = Dx[3];
518: dphidy[0 + offset] = Dy[0];
519: dphidy[1 + offset] = Dy[1];
520: dphidy[2 + offset] = Dy[2];
521: dphidy[3 + offset] = Dy[3];
523: dphidz[0 + offset] = Dz[0];
524: dphidz[1 + offset] = Dz[1];
525: dphidz[2 + offset] = Dz[2];
526: dphidz[3 + offset] = Dz[3];
527: }
529: } /* Tetrahedra -- ends */
530: } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "The number of entity vertices are invalid. Currently only support HEX8 and TET4 basis evaluations in 3-D : %" PetscInt_FMT, nverts);
531: return 0;
532: }
534: /*@C
535: DMMoabFEMComputeBasis - Evaluate bases and derivatives at quadrature points for a linear EDGE/QUAD/TRI/HEX/TET element.
537: The routine takes the coordinates of the vertices of an element and computes basis functions associated with
538: each quadrature point provided, and their derivatives with respect to X, Y and Z as appropriate.
540: Input Parameters:
541: + PetscInt nverts - the number of element vertices
542: . PetscReal coords[3*nverts] - the physical coordinates of the vertices (in canonical numbering)
543: . PetscInt npts - the number of evaluation points (quadrature points)
544: - PetscReal quad[3*npts] - the evaluation points (quadrature points) in the reference space
546: Output Parameters:
547: + PetscReal phypts[3*npts] - the evaluation points (quadrature points) transformed to the physical space
548: . PetscReal jxw[npts] - the jacobian determinant * quadrature weight necessary for assembling discrete contributions
549: . PetscReal fe_basis[npts] - the bases values evaluated at the specified quadrature points
550: - 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
552: Level: advanced
554: @*/
555: PetscErrorCode DMMoabFEMComputeBasis(const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscQuadrature quadrature, PetscReal *phypts, PetscReal *jacobian_quadrature_weight_product, PetscReal *fe_basis, PetscReal **fe_basis_derivatives)
556: {
557: PetscInt npoints, idim;
558: bool compute_der;
559: const PetscReal *quadpts, *quadwts;
560: PetscReal jacobian[9], ijacobian[9], volume;
565: compute_der = (fe_basis_derivatives != NULL);
567: /* Get the quadrature points and weights for the given quadrature rule */
568: PetscQuadratureGetData(quadrature, &idim, NULL, &npoints, &quadpts, &quadwts);
570: if (jacobian_quadrature_weight_product) PetscArraycpy(jacobian_quadrature_weight_product, quadwts, npoints);
572: switch (dim) {
573: case 1:
574: Compute_Lagrange_Basis_1D_Internal(nverts, coordinates, npoints, quadpts, phypts, jacobian_quadrature_weight_product, fe_basis, (compute_der ? fe_basis_derivatives[0] : NULL), jacobian, ijacobian, &volume);
575: break;
576: case 2:
577: Compute_Lagrange_Basis_2D_Internal(nverts, coordinates, npoints, quadpts, phypts, jacobian_quadrature_weight_product, fe_basis, (compute_der ? fe_basis_derivatives[0] : NULL), (compute_der ? fe_basis_derivatives[1] : NULL), jacobian, ijacobian, &volume);
578: break;
579: case 3:
580: Compute_Lagrange_Basis_3D_Internal(nverts, coordinates, npoints, quadpts, phypts, jacobian_quadrature_weight_product, fe_basis, (compute_der ? fe_basis_derivatives[0] : NULL), (compute_der ? fe_basis_derivatives[1] : NULL), (compute_der ? fe_basis_derivatives[2] : NULL), jacobian, ijacobian, &volume);
581: break;
582: default:
583: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %" PetscInt_FMT, dim);
584: }
585: return 0;
586: }
588: /*@C
589: DMMoabFEMCreateQuadratureDefault - Create default quadrature rules for integration over an element with a given
590: dimension and polynomial order (deciphered from number of element vertices).
592: Input Parameters:
594: + PetscInt dim - the element dimension (1=EDGE, 2=QUAD/TRI, 3=HEX/TET)
595: - PetscInt nverts - the number of vertices in the physical element
597: Output Parameter:
598: . PetscQuadrature quadrature - the quadrature object with default settings to integrate polynomials defined over the element
600: Level: advanced
602: @*/
603: PetscErrorCode DMMoabFEMCreateQuadratureDefault(const PetscInt dim, const PetscInt nverts, PetscQuadrature *quadrature)
604: {
605: PetscReal *w, *x;
606: PetscInt nc = 1;
608: /* Create an appropriate quadrature rule to sample basis */
609: switch (dim) {
610: case 1:
611: /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */
612: PetscDTStroudConicalQuadrature(1, nc, nverts, 0, 1.0, quadrature);
613: break;
614: case 2:
615: /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */
616: if (nverts == 3) {
617: const PetscInt order = 2;
618: const PetscInt npoints = (order == 2 ? 3 : 6);
619: PetscMalloc2(npoints * 2, &x, npoints, &w);
620: if (npoints == 3) {
621: x[0] = x[1] = x[2] = x[5] = 1.0 / 6.0;
622: x[3] = x[4] = 2.0 / 3.0;
623: w[0] = w[1] = w[2] = 1.0 / 3.0;
624: } else if (npoints == 6) {
625: x[0] = x[1] = x[2] = 0.44594849091597;
626: x[3] = x[4] = 0.10810301816807;
627: x[5] = 0.44594849091597;
628: x[6] = x[7] = x[8] = 0.09157621350977;
629: x[9] = x[10] = 0.81684757298046;
630: x[11] = 0.09157621350977;
631: w[0] = w[1] = w[2] = 0.22338158967801;
632: w[3] = w[4] = w[5] = 0.10995174365532;
633: } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Triangle quadrature rules for points 3 and 6 supported; npoints : %" PetscInt_FMT, npoints);
634: PetscQuadratureCreate(PETSC_COMM_SELF, quadrature);
635: PetscQuadratureSetOrder(*quadrature, order);
636: PetscQuadratureSetData(*quadrature, dim, nc, npoints, x, w);
637: /* PetscDTStroudConicalQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature); */
638: } else PetscDTGaussTensorQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature);
639: break;
640: case 3:
641: /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */
642: if (nverts == 4) {
643: PetscInt order;
644: const PetscInt npoints = 4; // Choose between 4 and 10
645: PetscMalloc2(npoints * 3, &x, npoints, &w);
646: if (npoints == 4) { /* KEAST1, K1, N=4, O=4 */
647: const PetscReal x_4[12] = {0.5854101966249685, 0.1381966011250105, 0.1381966011250105, 0.1381966011250105, 0.1381966011250105, 0.1381966011250105,
648: 0.1381966011250105, 0.1381966011250105, 0.5854101966249685, 0.1381966011250105, 0.5854101966249685, 0.1381966011250105};
649: PetscArraycpy(x, x_4, 12);
651: w[0] = w[1] = w[2] = w[3] = 1.0 / 24.0;
652: order = 4;
653: } else if (npoints == 10) { /* KEAST3, K3 N=10, O=10 */
654: const PetscReal x_10[30] = {0.5684305841968444, 0.1438564719343852, 0.1438564719343852, 0.1438564719343852, 0.1438564719343852, 0.1438564719343852, 0.1438564719343852, 0.1438564719343852, 0.5684305841968444, 0.1438564719343852,
655: 0.5684305841968444, 0.1438564719343852, 0.0000000000000000, 0.5000000000000000, 0.5000000000000000, 0.5000000000000000, 0.0000000000000000, 0.5000000000000000, 0.5000000000000000, 0.5000000000000000,
656: 0.0000000000000000, 0.5000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.5000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.5000000000000000};
657: PetscArraycpy(x, x_10, 30);
659: w[0] = w[1] = w[2] = w[3] = 0.2177650698804054;
660: w[4] = w[5] = w[6] = w[7] = w[8] = w[9] = 0.0214899534130631;
661: order = 10;
662: } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Tetrahedral quadrature rules for points 4 and 10 supported; npoints : %" PetscInt_FMT, npoints);
663: PetscQuadratureCreate(PETSC_COMM_SELF, quadrature);
664: PetscQuadratureSetOrder(*quadrature, order);
665: PetscQuadratureSetData(*quadrature, dim, nc, npoints, x, w);
666: /* PetscDTStroudConicalQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature); */
667: } else PetscDTGaussTensorQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature);
668: break;
669: default:
670: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %" PetscInt_FMT, dim);
671: }
672: return 0;
673: }
675: /* Compute Jacobians */
676: PetscErrorCode ComputeJacobian_Internal(const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal *quad, PetscReal *phypts, PetscReal *jacobian, PetscReal *ijacobian, PetscReal *dvolume)
677: {
678: PetscInt i;
679: PetscReal volume = 1.0;
684: PetscArrayzero(jacobian, dim * dim);
685: if (ijacobian) PetscArrayzero(ijacobian, dim * dim);
686: if (phypts) PetscArrayzero(phypts, /*npts=1 * */ 3);
688: if (dim == 1) {
689: const PetscReal &r = quad[0];
690: if (nverts == 2) { /* Linear Edge */
691: const PetscReal dNi_dxi[2] = {-1.0, 1.0};
693: for (i = 0; i < nverts; ++i) {
694: const PetscReal *vertices = coordinates + i * 3;
695: jacobian[0] += dNi_dxi[i] * vertices[0];
696: }
697: } else if (nverts == 3) { /* Quadratic Edge */
698: const PetscReal dNi_dxi[3] = {4 * r - 3.0, 4 * (1.0 - 2.0 * r), 4.0 * r - 1.0};
700: for (i = 0; i < nverts; ++i) {
701: const PetscReal *vertices = coordinates + i * 3;
702: jacobian[0] += dNi_dxi[i] * vertices[0];
703: }
704: } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "The number of 1-D entity vertices are invalid. Currently only support EDGE2 and EDGE3 basis evaluations in 1-D : %" PetscInt_FMT, nverts);
706: if (ijacobian) {
707: /* invert the jacobian */
708: ijacobian[0] = 1.0 / jacobian[0];
709: }
710: } else if (dim == 2) {
711: if (nverts == 4) { /* Linear Quadrangle */
712: const PetscReal &r = quad[0];
713: const PetscReal &s = quad[1];
715: const PetscReal dNi_dxi[4] = {-1.0 + s, 1.0 - s, s, -s};
716: const PetscReal dNi_deta[4] = {-1.0 + r, -r, r, 1.0 - r};
718: for (i = 0; i < nverts; ++i) {
719: const PetscReal *vertices = coordinates + i * 3;
720: jacobian[0] += dNi_dxi[i] * vertices[0];
721: jacobian[2] += dNi_dxi[i] * vertices[1];
722: jacobian[1] += dNi_deta[i] * vertices[0];
723: jacobian[3] += dNi_deta[i] * vertices[1];
724: }
725: } else if (nverts == 3) { /* Linear triangle */
726: const PetscReal x2 = coordinates[2 * 3 + 0], y2 = coordinates[2 * 3 + 1];
728: /* Jacobian is constant */
729: jacobian[0] = (coordinates[0 * 3 + 0] - x2);
730: jacobian[1] = (coordinates[1 * 3 + 0] - x2);
731: jacobian[2] = (coordinates[0 * 3 + 1] - y2);
732: jacobian[3] = (coordinates[1 * 3 + 1] - y2);
733: } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "The number of 2-D entity vertices are invalid. Currently only support QUAD4 and TRI3 basis evaluations in 2-D : %" PetscInt_FMT, nverts);
735: /* invert the jacobian */
736: if (ijacobian) DMatrix_Invert_2x2_Internal(jacobian, ijacobian, &volume);
737: } else {
738: if (nverts == 8) { /* Linear Hexahedra */
739: const PetscReal &r = quad[0];
740: const PetscReal &s = quad[1];
741: const PetscReal &t = quad[2];
742: const PetscReal dNi_dxi[8] = {-(1.0 - s) * (1.0 - t), (1.0 - s) * (1.0 - t), (s) * (1.0 - t), -(s) * (1.0 - t), -(1.0 - s) * (t), (1.0 - s) * (t), (s) * (t), -(s) * (t)};
744: const PetscReal dNi_deta[8] = {-(1.0 - r) * (1.0 - t), -(r) * (1.0 - t), (r) * (1.0 - t), (1.0 - r) * (1.0 - t), -(1.0 - r) * (t), -(r) * (t), (r) * (t), (1.0 - r) * (t)};
746: const PetscReal dNi_dzeta[8] = {-(1.0 - r) * (1.0 - s), -(r) * (1.0 - s), -(r) * (s), -(1.0 - r) * (s), (1.0 - r) * (1.0 - s), (r) * (1.0 - s), (r) * (s), (1.0 - r) * (s)};
748: for (i = 0; i < nverts; ++i) {
749: const PetscReal *vertex = coordinates + i * 3;
750: jacobian[0] += dNi_dxi[i] * vertex[0];
751: jacobian[3] += dNi_dxi[i] * vertex[1];
752: jacobian[6] += dNi_dxi[i] * vertex[2];
753: jacobian[1] += dNi_deta[i] * vertex[0];
754: jacobian[4] += dNi_deta[i] * vertex[1];
755: jacobian[7] += dNi_deta[i] * vertex[2];
756: jacobian[2] += dNi_dzeta[i] * vertex[0];
757: jacobian[5] += dNi_dzeta[i] * vertex[1];
758: jacobian[8] += dNi_dzeta[i] * vertex[2];
759: }
760: } else if (nverts == 4) { /* Linear Tetrahedra */
761: const PetscReal x0 = coordinates[/*0 * 3 +*/ 0], y0 = coordinates[/*0 * 3 +*/ 1], z0 = coordinates[/*0 * 3 +*/ 2];
763: /* compute the jacobian */
764: jacobian[0] = coordinates[1 * 3 + 0] - x0;
765: jacobian[1] = coordinates[2 * 3 + 0] - x0;
766: jacobian[2] = coordinates[3 * 3 + 0] - x0;
767: jacobian[3] = coordinates[1 * 3 + 1] - y0;
768: jacobian[4] = coordinates[2 * 3 + 1] - y0;
769: jacobian[5] = coordinates[3 * 3 + 1] - y0;
770: jacobian[6] = coordinates[1 * 3 + 2] - z0;
771: jacobian[7] = coordinates[2 * 3 + 2] - z0;
772: jacobian[8] = coordinates[3 * 3 + 2] - z0;
773: } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "The number of 3-D entity vertices are invalid. Currently only support HEX8 and TET4 basis evaluations in 3-D : %" PetscInt_FMT, nverts);
775: if (ijacobian) {
776: /* invert the jacobian */
777: DMatrix_Invert_3x3_Internal(jacobian, ijacobian, &volume);
778: }
779: }
781: if (dvolume) *dvolume = volume;
782: return 0;
783: }
785: PetscErrorCode FEMComputeBasis_JandF(const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal *quadrature, PetscReal *phypts, PetscReal *phibasis, PetscReal *jacobian, PetscReal *ijacobian, PetscReal *volume)
786: {
787: switch (dim) {
788: case 1:
789: Compute_Lagrange_Basis_1D_Internal(nverts, coordinates, 1, quadrature, phypts, NULL, phibasis, NULL, jacobian, ijacobian, volume);
790: break;
791: case 2:
792: Compute_Lagrange_Basis_2D_Internal(nverts, coordinates, 1, quadrature, phypts, NULL, phibasis, NULL, NULL, jacobian, ijacobian, volume);
793: break;
794: case 3:
795: Compute_Lagrange_Basis_3D_Internal(nverts, coordinates, 1, quadrature, phypts, NULL, phibasis, NULL, NULL, NULL, jacobian, ijacobian, volume);
796: break;
797: default:
798: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %" PetscInt_FMT, dim);
799: }
800: return 0;
801: }
803: /*@C
804: DMMoabPToRMapping - Compute the mapping from the physical coordinate system for a given element to the
805: canonical reference element. In addition to finding the inverse mapping evaluation through Newton iteration,
806: the basis function at the parametric point is also evaluated optionally.
808: Input Parameters:
809: + PetscInt dim - the element dimension (1=EDGE, 2=QUAD/TRI, 3=HEX/TET)
810: . PetscInt nverts - the number of vertices in the physical element
811: . PetscReal coordinates - the coordinates of vertices in the physical element
812: - PetscReal[3] xphy - the coordinates of physical point for which natural coordinates (in reference frame) are sought
814: Output Parameters:
815: + PetscReal[3] natparam - the natural coordinates (in reference frame) corresponding to xphy
816: - PetscReal[nverts] phi - the basis functions evaluated at the natural coordinates (natparam)
818: Level: advanced
820: @*/
821: PetscErrorCode DMMoabPToRMapping(const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal *xphy, PetscReal *natparam, PetscReal *phi)
822: {
823: /* Perform inverse evaluation for the mapping with use of Newton Raphson iteration */
824: const PetscReal tol = 1.0e-10;
825: const PetscInt max_iterations = 10;
826: const PetscReal error_tol_sqr = tol * tol;
827: PetscReal phibasis[8], jacobian[9], ijacobian[9], volume;
828: PetscReal phypts[3] = {0.0, 0.0, 0.0};
829: PetscReal delta[3] = {0.0, 0.0, 0.0};
830: PetscInt iters = 0;
831: PetscReal error = 1.0;
837: PetscArrayzero(jacobian, dim * dim);
838: PetscArrayzero(ijacobian, dim * dim);
839: PetscArrayzero(phibasis, nverts);
841: /* zero initial guess */
842: natparam[0] = natparam[1] = natparam[2] = 0.0;
844: /* Compute delta = evaluate( xi) - x */
845: FEMComputeBasis_JandF(dim, nverts, coordinates, natparam, phypts, phibasis, jacobian, ijacobian, &volume);
847: error = 0.0;
848: switch (dim) {
849: case 3:
850: delta[2] = phypts[2] - xphy[2];
851: error += (delta[2] * delta[2]);
852: case 2:
853: delta[1] = phypts[1] - xphy[1];
854: error += (delta[1] * delta[1]);
855: case 1:
856: delta[0] = phypts[0] - xphy[0];
857: error += (delta[0] * delta[0]);
858: break;
859: }
861: while (error > error_tol_sqr) {
864: /* Compute natparam -= J.inverse() * delta */
865: switch (dim) {
866: case 1:
867: natparam[0] -= ijacobian[0] * delta[0];
868: break;
869: case 2:
870: natparam[0] -= ijacobian[0] * delta[0] + ijacobian[1] * delta[1];
871: natparam[1] -= ijacobian[2] * delta[0] + ijacobian[3] * delta[1];
872: break;
873: case 3:
874: natparam[0] -= ijacobian[0] * delta[0] + ijacobian[1] * delta[1] + ijacobian[2] * delta[2];
875: natparam[1] -= ijacobian[3] * delta[0] + ijacobian[4] * delta[1] + ijacobian[5] * delta[2];
876: natparam[2] -= ijacobian[6] * delta[0] + ijacobian[7] * delta[1] + ijacobian[8] * delta[2];
877: break;
878: }
880: /* Compute delta = evaluate( xi) - x */
881: FEMComputeBasis_JandF(dim, nverts, coordinates, natparam, phypts, phibasis, jacobian, ijacobian, &volume);
883: error = 0.0;
884: switch (dim) {
885: case 3:
886: delta[2] = phypts[2] - xphy[2];
887: error += (delta[2] * delta[2]);
888: case 2:
889: delta[1] = phypts[1] - xphy[1];
890: error += (delta[1] * delta[1]);
891: case 1:
892: delta[0] = phypts[0] - xphy[0];
893: error += (delta[0] * delta[0]);
894: break;
895: }
896: }
897: if (phi) PetscArraycpy(phi, phibasis, nverts);
898: return 0;
899: }