Actual source code: dmmbfem.cxx
petsc-3.14.6 2021-03-30
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 Parameter:
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 Parameter:
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
124: Level: advanced
126: @*/
127: PetscErrorCode Compute_Lagrange_Basis_1D_Internal(const PetscInt nverts, const PetscReal *coords, const PetscInt npts, const PetscReal *quad, PetscReal *phypts,
128: PetscReal *jxw, PetscReal *phi, PetscReal *dphidx,
129: PetscReal *jacobian, PetscReal *ijacobian, PetscReal *volume)
130: {
131: int i, j;
132: PetscErrorCode ierr;
138: if (phypts) {
139: PetscArrayzero(phypts, npts * 3);
140: }
141: if (dphidx) { /* Reset arrays. */
142: PetscArrayzero(dphidx, npts * nverts);
143: }
144: if (nverts == 2) { /* Linear Edge */
146: for (j = 0; j < npts; j++) {
147: const PetscInt offset = j * nverts;
148: const PetscReal r = quad[j];
150: phi[0 + offset] = ( 1.0 - r);
151: phi[1 + offset] = ( r);
153: const PetscReal dNi_dxi[2] = { -1.0, 1.0 };
155: jacobian[0] = ijacobian[0] = volume[0] = 0.0;
156: for (i = 0; i < nverts; ++i) {
157: const PetscReal* vertices = coords + i * 3;
158: jacobian[0] += dNi_dxi[i] * vertices[0];
159: if (phypts) {
160: phypts[3 * j + 0] += phi[i + offset] * vertices[0];
161: }
162: }
164: /* invert the jacobian */
165: *volume = jacobian[0];
166: ijacobian[0] = 1.0 / jacobian[0];
167: jxw[j] *= *volume;
169: /* Divide by element jacobian. */
170: for (i = 0; i < nverts; i++) {
171: if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[0];
172: }
173: }
174: } else if (nverts == 3) { /* Quadratic Edge */
176: for (j = 0; j < npts; j++) {
177: const PetscInt offset = j * nverts;
178: const PetscReal r = quad[j];
180: phi[0 + offset] = 1.0 + r * ( 2.0 * r - 3.0);
181: phi[1 + offset] = 4.0 * r * ( 1.0 - r);
182: phi[2 + offset] = r * ( 2.0 * r - 1.0);
184: const PetscReal dNi_dxi[3] = { 4 * r - 3.0, 4 * ( 1.0 - 2.0 * r), 4.0 * r - 1.0};
186: jacobian[0] = ijacobian[0] = volume[0] = 0.0;
187: for (i = 0; i < nverts; ++i) {
188: const PetscReal* vertices = coords + i * 3;
189: jacobian[0] += dNi_dxi[i] * vertices[0];
190: if (phypts) {
191: phypts[3 * j + 0] += phi[i + offset] * vertices[0];
192: }
193: }
195: /* invert the jacobian */
196: *volume = jacobian[0];
197: ijacobian[0] = 1.0 / jacobian[0];
198: if (jxw) jxw[j] *= *volume;
200: /* Divide by element jacobian. */
201: for (i = 0; i < nverts; i++) {
202: if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[0];
203: }
204: }
205: } else SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "The number of entity vertices are invalid. Currently only support EDGE2 and EDGE3 basis evaluations in 1-D : %D", nverts);
206: return(0);
207: }
209: /*@C
210: Compute_Lagrange_Basis_2D_Internal - Evaluate bases and derivatives at quadrature points for a QUAD4 or TRI3 element.
212: The routine is given the coordinates of the vertices of a quadrangle or triangle.
214: The routine evaluates the basis functions associated with each quadrature point provided,
215: and their derivatives with respect to X and Y.
217: Notes:
219: Example Physical Element (QUAD4)
220: .vb
221: 4------3 s
222: | | |
223: | | |
224: | | |
225: 1------2 0-------r
226: .ve
228: Input Parameter:
229: + PetscInt nverts - the number of element vertices
230: . PetscReal coords[3*nverts] - the physical coordinates of the vertices (in canonical numbering)
231: . PetscInt npts - the number of evaluation points (quadrature points)
232: - PetscReal quad[3*npts] - the evaluation points (quadrature points) in the reference space
234: Output Parameter:
235: + PetscReal phypts[3*npts] - the evaluation points (quadrature points) transformed to the physical space
236: . PetscReal jxw[npts] - the jacobian determinant * quadrature weight necessary for assembling discrete contributions
237: . PetscReal phi[npts] - the bases evaluated at the specified quadrature points
238: . PetscReal dphidx[npts] - the derivative of the bases wrt X-direction evaluated at the specified quadrature points
239: - PetscReal dphidy[npts] - the derivative of the bases wrt Y-direction evaluated at the specified quadrature points
241: Level: advanced
243: @*/
244: PetscErrorCode Compute_Lagrange_Basis_2D_Internal(const PetscInt nverts, const PetscReal *coords, const PetscInt npts, const PetscReal *quad, PetscReal *phypts,
245: PetscReal *jxw, PetscReal *phi, PetscReal *dphidx, PetscReal *dphidy,
246: PetscReal *jacobian, PetscReal *ijacobian, PetscReal *volume)
247: {
248: PetscInt i, j, k;
255: PetscArrayzero(phi, npts);
256: if (phypts) {
257: PetscArrayzero(phypts, npts * 3);
258: }
259: if (dphidx) { /* Reset arrays. */
260: PetscArrayzero(dphidx, npts * nverts);
261: PetscArrayzero(dphidy, npts * nverts);
262: }
263: if (nverts == 4) { /* Linear Quadrangle */
265: for (j = 0; j < npts; j++)
266: {
267: const PetscInt offset = j * nverts;
268: const PetscReal r = quad[0 + j * 2];
269: const PetscReal s = quad[1 + j * 2];
271: phi[0 + offset] = ( 1.0 - r) * ( 1.0 - s);
272: phi[1 + offset] = r * ( 1.0 - s);
273: phi[2 + offset] = r * s;
274: phi[3 + offset] = ( 1.0 - r) * s;
276: const PetscReal dNi_dxi[4] = { -1.0 + s, 1.0 - s, s, -s };
277: const PetscReal dNi_deta[4] = { -1.0 + r, -r, r, 1.0 - r };
279: PetscArrayzero(jacobian, 4);
280: PetscArrayzero(ijacobian, 4);
281: for (i = 0; i < nverts; ++i) {
282: const PetscReal* vertices = coords + i * 3;
283: jacobian[0] += dNi_dxi[i] * vertices[0];
284: jacobian[2] += dNi_dxi[i] * vertices[1];
285: jacobian[1] += dNi_deta[i] * vertices[0];
286: jacobian[3] += dNi_deta[i] * vertices[1];
287: if (phypts) {
288: phypts[3 * j + 0] += phi[i + offset] * vertices[0];
289: phypts[3 * j + 1] += phi[i + offset] * vertices[1];
290: phypts[3 * j + 2] += phi[i + offset] * vertices[2];
291: }
292: }
294: /* invert the jacobian */
295: DMatrix_Invert_2x2_Internal(jacobian, ijacobian, volume);
296: 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);
298: if (jxw) jxw[j] *= *volume;
300: /* Let us compute the bases derivatives by scaling with inverse jacobian. */
301: if (dphidx) {
302: for (i = 0; i < nverts; i++) {
303: for (k = 0; k < 2; ++k) {
304: if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[k * 2 + 0];
305: if (dphidy) dphidy[i + offset] += dNi_deta[i] * ijacobian[k * 2 + 1];
306: } /* for k=[0..2) */
307: } /* for i=[0..nverts) */
308: } /* if (dphidx) */
309: }
310: } else if (nverts == 3) { /* Linear triangle */
311: const PetscReal x2 = coords[2 * 3 + 0], y2 = coords[2 * 3 + 1];
313: PetscArrayzero(jacobian, 4);
314: PetscArrayzero(ijacobian, 4);
316: /* Jacobian is constant */
317: jacobian[0] = (coords[0 * 3 + 0] - x2); jacobian[1] = (coords[1 * 3 + 0] - x2);
318: jacobian[2] = (coords[0 * 3 + 1] - y2); jacobian[3] = (coords[1 * 3 + 1] - y2);
320: /* invert the jacobian */
321: DMatrix_Invert_2x2_Internal(jacobian, ijacobian, volume);
322: if (*volume < 1e-12) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Triangular element has zero volume: %g. Degenerate element or invalid connectivity\n", (double)*volume);
324: const PetscReal Dx[3] = { ijacobian[0], ijacobian[2], - ijacobian[0] - ijacobian[2] };
325: const PetscReal Dy[3] = { ijacobian[1], ijacobian[3], - ijacobian[1] - ijacobian[3] };
327: for (j = 0; j < npts; j++) {
328: const PetscInt offset = j * nverts;
329: const PetscReal r = quad[0 + j * 2];
330: const PetscReal s = quad[1 + j * 2];
332: if (jxw) jxw[j] *= 0.5;
334: /* const PetscReal Ni[3] = { r, s, 1.0 - r - s }; */
335: const PetscReal phipts_x = coords[2 * 3 + 0] + jacobian[0] * r + jacobian[1] * s;
336: const PetscReal phipts_y = coords[2 * 3 + 1] + jacobian[2] * r + jacobian[3] * s;
337: if (phypts) {
338: phypts[offset + 0] = phipts_x;
339: phypts[offset + 1] = phipts_y;
340: }
342: /* \phi_0 = (b.y - c.y) x + (b.x - c.x) y + c.x b.y - b.x c.y */
343: phi[0 + offset] = ( ijacobian[0] * (phipts_x - x2) + ijacobian[1] * (phipts_y - y2));
344: /* \phi_1 = (c.y - a.y) x + (a.x - c.x) y + c.x a.y - a.x c.y */
345: phi[1 + offset] = ( ijacobian[2] * (phipts_x - x2) + ijacobian[3] * (phipts_y - y2));
346: phi[2 + offset] = 1.0 - phi[0 + offset] - phi[1 + offset];
348: if (dphidx) {
349: dphidx[0 + offset] = Dx[0];
350: dphidx[1 + offset] = Dx[1];
351: dphidx[2 + offset] = Dx[2];
352: }
354: if (dphidy) {
355: dphidy[0 + offset] = Dy[0];
356: dphidy[1 + offset] = Dy[1];
357: dphidy[2 + offset] = Dy[2];
358: }
360: }
361: } else SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "The number of entity vertices are invalid. Currently only support QUAD4 and TRI3 basis evaluations in 2-D : %D", nverts);
362: return(0);
363: }
365: /*@C
366: Compute_Lagrange_Basis_3D_Internal - Evaluate bases and derivatives at quadrature points for a HEX8 or TET4 element.
368: The routine is given the coordinates of the vertices of a hexahedra or tetrahedra.
370: The routine evaluates the basis functions associated with each quadrature point provided,
371: and their derivatives with respect to X, Y and Z.
373: Notes:
375: Example Physical Element (HEX8)
376: .vb
377: 8------7
378: /| /| t s
379: 5------6 | | /
380: | | | | |/
381: | 4----|-3 0-------r
382: |/ |/
383: 1------2
384: .ve
386: Input Parameter:
387: + PetscInt nverts - the number of element vertices
388: . PetscReal coords[3*nverts] - the physical coordinates of the vertices (in canonical numbering)
389: . PetscInt npts - the number of evaluation points (quadrature points)
390: - PetscReal quad[3*npts] - the evaluation points (quadrature points) in the reference space
392: Output Parameter:
393: + PetscReal phypts[3*npts] - the evaluation points (quadrature points) transformed to the physical space
394: . PetscReal jxw[npts] - the jacobian determinant * quadrature weight necessary for assembling discrete contributions
395: . PetscReal phi[npts] - the bases evaluated at the specified quadrature points
396: . PetscReal dphidx[npts] - the derivative of the bases wrt X-direction evaluated at the specified quadrature points
397: . PetscReal dphidy[npts] - the derivative of the bases wrt Y-direction evaluated at the specified quadrature points
398: - PetscReal dphidz[npts] - the derivative of the bases wrt Z-direction evaluated at the specified quadrature points
400: Level: advanced
402: @*/
403: PetscErrorCode Compute_Lagrange_Basis_3D_Internal(const PetscInt nverts, const PetscReal *coords, const PetscInt npts, const PetscReal *quad, PetscReal *phypts,
404: PetscReal *jxw, PetscReal *phi, PetscReal *dphidx, PetscReal *dphidy, PetscReal *dphidz,
405: PetscReal *jacobian, PetscReal *ijacobian, PetscReal *volume)
406: {
407: PetscInt i, j, k;
415: PetscArrayzero(phi, npts);
416: if (phypts) {
417: PetscArrayzero(phypts, npts * 3);
418: }
419: if (dphidx) {
420: PetscArrayzero(dphidx, npts * nverts);
421: PetscArrayzero(dphidy, npts * nverts);
422: PetscArrayzero(dphidz, npts * nverts);
423: }
425: if (nverts == 8) { /* Linear Hexahedra */
427: for (j = 0; j < npts; j++) {
428: const PetscInt offset = j * nverts;
429: const PetscReal& r = quad[j * 3 + 0];
430: const PetscReal& s = quad[j * 3 + 1];
431: const PetscReal& t = quad[j * 3 + 2];
433: phi[offset + 0] = ( 1.0 - r) * ( 1.0 - s) * ( 1.0 - t); /* 0,0,0 */
434: phi[offset + 1] = ( r) * ( 1.0 - s) * ( 1.0 - t); /* 1,0,0 */
435: phi[offset + 2] = ( r) * ( s) * ( 1.0 - t); /* 1,1,0 */
436: phi[offset + 3] = ( 1.0 - r) * ( s) * ( 1.0 - t); /* 0,1,0 */
437: phi[offset + 4] = ( 1.0 - r) * ( 1.0 - s) * ( t); /* 0,0,1 */
438: phi[offset + 5] = ( r) * ( 1.0 - s) * ( t); /* 1,0,1 */
439: phi[offset + 6] = ( r) * ( s) * ( t); /* 1,1,1 */
440: phi[offset + 7] = ( 1.0 - r) * ( s) * ( t); /* 0,1,1 */
442: const PetscReal dNi_dxi[8] = {- ( 1.0 - s) * ( 1.0 - t),
443: ( 1.0 - s) * ( 1.0 - t),
444: ( s) * ( 1.0 - t),
445: - ( s) * ( 1.0 - t),
446: - ( 1.0 - s) * ( t),
447: ( 1.0 - s) * ( t),
448: ( s) * ( t),
449: - ( s) * ( t)
450: };
452: const PetscReal dNi_deta[8] = { - ( 1.0 - r) * ( 1.0 - t),
453: - ( r) * ( 1.0 - t),
454: ( r) * ( 1.0 - t),
455: ( 1.0 - r) * ( 1.0 - t),
456: - ( 1.0 - r) * ( t),
457: - ( r) * ( t),
458: ( r) * ( t),
459: ( 1.0 - r) * ( t)
460: };
462: const PetscReal dNi_dzeta[8] = { - ( 1.0 - r) * ( 1.0 - s),
463: - ( r) * ( 1.0 - s),
464: - ( r) * ( s),
465: - ( 1.0 - r) * ( s),
466: ( 1.0 - r) * ( 1.0 - s),
467: ( r) * ( 1.0 - s),
468: ( r) * ( s),
469: ( 1.0 - r) * ( s)
470: };
472: PetscArrayzero(jacobian, 9);
473: PetscArrayzero(ijacobian, 9);
474: for (i = 0; i < nverts; ++i) {
475: const PetscReal* vertex = coords + i * 3;
476: jacobian[0] += dNi_dxi[i] * vertex[0];
477: jacobian[3] += dNi_dxi[i] * vertex[1];
478: jacobian[6] += dNi_dxi[i] * vertex[2];
479: jacobian[1] += dNi_deta[i] * vertex[0];
480: jacobian[4] += dNi_deta[i] * vertex[1];
481: jacobian[7] += dNi_deta[i] * vertex[2];
482: jacobian[2] += dNi_dzeta[i] * vertex[0];
483: jacobian[5] += dNi_dzeta[i] * vertex[1];
484: jacobian[8] += dNi_dzeta[i] * vertex[2];
485: if (phypts) {
486: phypts[3 * j + 0] += phi[i + offset] * vertex[0];
487: phypts[3 * j + 1] += phi[i + offset] * vertex[1];
488: phypts[3 * j + 2] += phi[i + offset] * vertex[2];
489: }
490: }
492: /* invert the jacobian */
493: DMatrix_Invert_3x3_Internal(jacobian, ijacobian, volume);
494: 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);
496: if (jxw) jxw[j] *= (*volume);
498: /* Divide by element jacobian. */
499: for (i = 0; i < nverts; ++i) {
500: for (k = 0; k < 3; ++k) {
501: if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[0 * 3 + k];
502: if (dphidy) dphidy[i + offset] += dNi_deta[i] * ijacobian[1 * 3 + k];
503: if (dphidz) dphidz[i + offset] += dNi_dzeta[i] * ijacobian[2 * 3 + k];
504: }
505: }
506: }
507: } else if (nverts == 4) { /* Linear Tetrahedra */
508: PetscReal Dx[4]={0,0,0,0}, Dy[4]={0,0,0,0}, Dz[4]={0,0,0,0};
509: const PetscReal x0 = coords[/*0 * 3 +*/ 0], y0 = coords[/*0 * 3 +*/ 1], z0 = coords[/*0 * 3 +*/ 2];
511: PetscArrayzero(jacobian, 9);
512: PetscArrayzero(ijacobian, 9);
514: /* compute the jacobian */
515: jacobian[0] = coords[1 * 3 + 0] - x0; jacobian[1] = coords[2 * 3 + 0] - x0; jacobian[2] = coords[3 * 3 + 0] - x0;
516: jacobian[3] = coords[1 * 3 + 1] - y0; jacobian[4] = coords[2 * 3 + 1] - y0; jacobian[5] = coords[3 * 3 + 1] - y0;
517: jacobian[6] = coords[1 * 3 + 2] - z0; jacobian[7] = coords[2 * 3 + 2] - z0; jacobian[8] = coords[3 * 3 + 2] - z0;
519: /* invert the jacobian */
520: DMatrix_Invert_3x3_Internal(jacobian, ijacobian, volume);
521: if (*volume < 1e-8) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Tetrahedral element has zero volume: %g. Degenerate element or invalid connectivity\n", (double)*volume);
523: if (dphidx) {
524: Dx[0] = ( coords[1 + 2 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3])
525: - coords[1 + 1 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3])
526: - coords[1 + 3 * 3] * ( coords[2 + 1 * 3] - coords[2 + 2 * 3])) / *volume;
527: Dx[1] = - ( coords[1 + 2 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3])
528: - coords[1 + 0 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3])
529: - coords[1 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 2 * 3])) / *volume;
530: Dx[2] = ( coords[1 + 1 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3])
531: - coords[1 + 0 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3])
532: - coords[1 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 1 * 3])) / *volume;
533: Dx[3] = - ( Dx[0] + Dx[1] + Dx[2]);
534: Dy[0] = ( coords[0 + 1 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3])
535: - coords[0 + 2 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3])
536: + coords[0 + 3 * 3] * ( coords[2 + 1 * 3] - coords[2 + 2 * 3])) / *volume;
537: Dy[1] = - ( coords[0 + 0 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3])
538: - coords[0 + 2 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3])
539: + coords[0 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 2 * 3])) / *volume;
540: Dy[2] = ( coords[0 + 0 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3])
541: - coords[0 + 1 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3])
542: + coords[0 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 1 * 3])) / *volume;
543: Dy[3] = - ( Dy[0] + Dy[1] + Dy[2]);
544: Dz[0] = ( coords[0 + 1 * 3] * (coords[1 + 3 * 3] - coords[1 + 2 * 3])
545: - coords[0 + 2 * 3] * (coords[1 + 3 * 3] - coords[1 + 1 * 3])
546: + coords[0 + 3 * 3] * (coords[1 + 2 * 3] - coords[1 + 1 * 3])) / *volume;
547: Dz[1] = - ( coords[0 + 0 * 3] * (coords[1 + 3 * 3] - coords[1 + 2 * 3])
548: + coords[0 + 2 * 3] * (coords[1 + 0 * 3] - coords[1 + 3 * 3])
549: - coords[0 + 3 * 3] * (coords[1 + 0 * 3] - coords[1 + 2 * 3])) / *volume;
550: Dz[2] = ( coords[0 + 0 * 3] * (coords[1 + 3 * 3] - coords[1 + 1 * 3])
551: + coords[0 + 1 * 3] * (coords[1 + 0 * 3] - coords[1 + 3 * 3])
552: - coords[0 + 3 * 3] * (coords[1 + 0 * 3] - coords[1 + 1 * 3])) / *volume;
553: Dz[3] = - ( Dz[0] + Dz[1] + Dz[2]);
554: }
556: for (j = 0; j < npts; j++) {
557: const PetscInt offset = j * nverts;
558: const PetscReal& r = quad[j * 3 + 0];
559: const PetscReal& s = quad[j * 3 + 1];
560: const PetscReal& t = quad[j * 3 + 2];
562: if (jxw) jxw[j] *= *volume;
564: phi[offset + 0] = 1.0 - r - s - t;
565: phi[offset + 1] = r;
566: phi[offset + 2] = s;
567: phi[offset + 3] = t;
569: if (phypts) {
570: for (i = 0; i < nverts; ++i) {
571: const PetscScalar* vertices = coords + i * 3;
572: phypts[3 * j + 0] += phi[i + offset] * vertices[0];
573: phypts[3 * j + 1] += phi[i + offset] * vertices[1];
574: phypts[3 * j + 2] += phi[i + offset] * vertices[2];
575: }
576: }
578: /* Now set the derivatives */
579: if (dphidx) {
580: dphidx[0 + offset] = Dx[0];
581: dphidx[1 + offset] = Dx[1];
582: dphidx[2 + offset] = Dx[2];
583: dphidx[3 + offset] = Dx[3];
585: dphidy[0 + offset] = Dy[0];
586: dphidy[1 + offset] = Dy[1];
587: dphidy[2 + offset] = Dy[2];
588: dphidy[3 + offset] = Dy[3];
590: dphidz[0 + offset] = Dz[0];
591: dphidz[1 + offset] = Dz[1];
592: dphidz[2 + offset] = Dz[2];
593: dphidz[3 + offset] = Dz[3];
594: }
596: } /* Tetrahedra -- ends */
597: } else SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "The number of entity vertices are invalid. Currently only support HEX8 and TET4 basis evaluations in 3-D : %D", nverts);
598: return(0);
599: }
601: /*@C
602: DMMoabFEMComputeBasis - Evaluate bases and derivatives at quadrature points for a linear EDGE/QUAD/TRI/HEX/TET element.
604: The routine takes the coordinates of the vertices of an element and computes basis functions associated with
605: each quadrature point provided, and their derivatives with respect to X, Y and Z as appropriate.
607: Input Parameter:
608: + PetscInt nverts, the number of element vertices
609: . PetscReal coords[3*nverts], the physical coordinates of the vertices (in canonical numbering)
610: . PetscInt npts, the number of evaluation points (quadrature points)
611: - PetscReal quad[3*npts], the evaluation points (quadrature points) in the reference space
613: Output Parameter:
614: + PetscReal phypts[3*npts], the evaluation points (quadrature points) transformed to the physical space
615: . PetscReal jxw[npts], the jacobian determinant * quadrature weight necessary for assembling discrete contributions
616: . PetscReal fe_basis[npts], the bases values evaluated at the specified quadrature points
617: - 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
619: Level: advanced
621: @*/
622: PetscErrorCode DMMoabFEMComputeBasis(const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscQuadrature quadrature,
623: PetscReal *phypts, PetscReal *jacobian_quadrature_weight_product,
624: PetscReal *fe_basis, PetscReal **fe_basis_derivatives)
625: {
626: PetscErrorCode ierr;
627: PetscInt npoints,idim;
628: bool compute_der;
629: const PetscReal *quadpts, *quadwts;
630: PetscReal jacobian[9], ijacobian[9], volume;
636: compute_der = (fe_basis_derivatives != NULL);
638: /* Get the quadrature points and weights for the given quadrature rule */
639: PetscQuadratureGetData(quadrature, &idim, NULL, &npoints, &quadpts, &quadwts);
640: if (idim != dim) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Dimension mismatch: provided (%D) vs quadrature (%D)\n",idim,dim);
641: if (jacobian_quadrature_weight_product) {
642: PetscArraycpy(jacobian_quadrature_weight_product, quadwts, npoints);
643: }
645: switch (dim) {
646: case 1:
647: Compute_Lagrange_Basis_1D_Internal(nverts, coordinates, npoints, quadpts, phypts,
648: jacobian_quadrature_weight_product, fe_basis,
649: (compute_der ? fe_basis_derivatives[0] : NULL),
650: jacobian, ijacobian, &volume);
651: break;
652: case 2:
653: Compute_Lagrange_Basis_2D_Internal(nverts, coordinates, npoints, quadpts, phypts,
654: jacobian_quadrature_weight_product, fe_basis,
655: (compute_der ? fe_basis_derivatives[0] : NULL),
656: (compute_der ? fe_basis_derivatives[1] : NULL),
657: jacobian, ijacobian, &volume);
658: break;
659: case 3:
660: Compute_Lagrange_Basis_3D_Internal(nverts, coordinates, npoints, quadpts, phypts,
661: jacobian_quadrature_weight_product, fe_basis,
662: (compute_der ? fe_basis_derivatives[0] : NULL),
663: (compute_der ? fe_basis_derivatives[1] : NULL),
664: (compute_der ? fe_basis_derivatives[2] : NULL),
665: jacobian, ijacobian, &volume);
666: break;
667: default:
668: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %D", dim);
669: }
670: return(0);
671: }
673: /*@C
674: DMMoabFEMCreateQuadratureDefault - Create default quadrature rules for integration over an element with a given
675: dimension and polynomial order (deciphered from number of element vertices).
677: Input Parameter:
679: + PetscInt dim - the element dimension (1=EDGE, 2=QUAD/TRI, 3=HEX/TET)
680: - PetscInt nverts - the number of vertices in the physical element
682: Output Parameter:
683: . PetscQuadrature quadrature - the quadrature object with default settings to integrate polynomials defined over the element
685: Level: advanced
687: @*/
688: PetscErrorCode DMMoabFEMCreateQuadratureDefault(const PetscInt dim, const PetscInt nverts, PetscQuadrature *quadrature)
689: {
690: PetscReal *w, *x;
691: PetscInt nc=1;
692: PetscErrorCode ierr;
695: /* Create an appropriate quadrature rule to sample basis */
696: switch (dim)
697: {
698: case 1:
699: /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */
700: PetscDTStroudConicalQuadrature(1, nc, nverts, 0, 1.0, quadrature);
701: break;
702: case 2:
703: /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */
704: if (nverts == 3) {
705: const PetscInt order = 2;
706: const PetscInt npoints = (order == 2 ? 3 : 6);
707: PetscMalloc2(npoints * 2, &x, npoints, &w);
708: if (npoints == 3) {
709: x[0] = x[1] = x[2] = x[5] = 1.0 / 6.0;
710: x[3] = x[4] = 2.0 / 3.0;
711: w[0] = w[1] = w[2] = 1.0 / 3.0;
712: } else if (npoints == 6) {
713: x[0] = x[1] = x[2] = 0.44594849091597;
714: x[3] = x[4] = 0.10810301816807;
715: x[5] = 0.44594849091597;
716: x[6] = x[7] = x[8] = 0.09157621350977;
717: x[9] = x[10] = 0.81684757298046;
718: x[11] = 0.09157621350977;
719: w[0] = w[1] = w[2] = 0.22338158967801;
720: w[3] = w[4] = w[5] = 0.10995174365532;
721: } else SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Triangle quadrature rules for points 3 and 6 supported; npoints : %D", npoints);
722: PetscQuadratureCreate(PETSC_COMM_SELF, quadrature);
723: PetscQuadratureSetOrder(*quadrature, order);
724: PetscQuadratureSetData(*quadrature, dim, nc, npoints, x, w);
725: /* PetscDTStroudConicalQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature); */
726: } else {
727: PetscDTGaussTensorQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature);
728: }
729: break;
730: case 3:
731: /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */
732: if (nverts == 4) {
733: PetscInt order;
734: const PetscInt npoints = 4; // Choose between 4 and 10
735: PetscMalloc2(npoints * 3, &x, npoints, &w);
736: if (npoints == 4) { /* KEAST1, K1, N=4, O=4 */
737: const PetscReal x_4[12] = { 0.5854101966249685, 0.1381966011250105, 0.1381966011250105,
738: 0.1381966011250105, 0.1381966011250105, 0.1381966011250105,
739: 0.1381966011250105, 0.1381966011250105, 0.5854101966249685,
740: 0.1381966011250105, 0.5854101966249685, 0.1381966011250105
741: };
742: PetscArraycpy(x, x_4, 12);
744: w[0] = w[1] = w[2] = w[3] = 1.0 / 24.0;
745: order = 4;
746: } else if (npoints == 10) { /* KEAST3, K3 N=10, O=10 */
747: const PetscReal x_10[30] = { 0.5684305841968444, 0.1438564719343852, 0.1438564719343852,
748: 0.1438564719343852, 0.1438564719343852, 0.1438564719343852,
749: 0.1438564719343852, 0.1438564719343852, 0.5684305841968444,
750: 0.1438564719343852, 0.5684305841968444, 0.1438564719343852,
751: 0.0000000000000000, 0.5000000000000000, 0.5000000000000000,
752: 0.5000000000000000, 0.0000000000000000, 0.5000000000000000,
753: 0.5000000000000000, 0.5000000000000000, 0.0000000000000000,
754: 0.5000000000000000, 0.0000000000000000, 0.0000000000000000,
755: 0.0000000000000000, 0.5000000000000000, 0.0000000000000000,
756: 0.0000000000000000, 0.0000000000000000, 0.5000000000000000
757: };
758: PetscArraycpy(x, x_10, 30);
760: w[0] = w[1] = w[2] = w[3] = 0.2177650698804054;
761: w[4] = w[5] = w[6] = w[7] = w[8] = w[9] = 0.0214899534130631;
762: order = 10;
763: } else SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Tetrahedral quadrature rules for points 4 and 10 supported; npoints : %D", npoints);
764: PetscQuadratureCreate(PETSC_COMM_SELF, quadrature);
765: PetscQuadratureSetOrder(*quadrature, order);
766: PetscQuadratureSetData(*quadrature, dim, nc, npoints, x, w);
767: /* PetscDTStroudConicalQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature); */
768: } else {
769: PetscDTGaussTensorQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature);
770: }
771: break;
772: default:
773: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %D", dim);
774: }
775: return(0);
776: }
778: /* Compute Jacobians */
779: PetscErrorCode ComputeJacobian_Internal (const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal *quad, PetscReal *phypts,
780: PetscReal *jacobian, PetscReal *ijacobian, PetscReal* dvolume)
781: {
782: PetscInt i;
783: PetscReal volume=1.0;
790: PetscArrayzero(jacobian, dim * dim);
791: if (ijacobian) {
792: PetscArrayzero(ijacobian, dim * dim);
793: }
794: if (phypts) {
795: PetscArrayzero(phypts, /*npts=1 * */ 3);
796: }
798: if (dim == 1) {
799: const PetscReal& r = quad[0];
800: if (nverts == 2) { /* Linear Edge */
801: const PetscReal dNi_dxi[2] = { -1.0, 1.0 };
803: for (i = 0; i < nverts; ++i) {
804: const PetscReal* vertices = coordinates + i * 3;
805: jacobian[0] += dNi_dxi[i] * vertices[0];
806: }
807: } else if (nverts == 3) { /* Quadratic Edge */
808: const PetscReal dNi_dxi[3] = { 4 * r - 3.0, 4 * ( 1.0 - 2.0 * r), 4.0 * r - 1.0};
810: for (i = 0; i < nverts; ++i) {
811: const PetscReal* vertices = coordinates + i * 3;
812: jacobian[0] += dNi_dxi[i] * vertices[0];
813: }
814: } else {
815: 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);
816: }
818: if (ijacobian) {
819: /* invert the jacobian */
820: ijacobian[0] = 1.0 / jacobian[0];
821: }
822: } else if (dim == 2) {
824: if (nverts == 4) { /* Linear Quadrangle */
825: const PetscReal& r = quad[0];
826: const PetscReal& s = quad[1];
828: const PetscReal dNi_dxi[4] = { -1.0 + s, 1.0 - s, s, -s };
829: const PetscReal dNi_deta[4] = { -1.0 + r, -r, r, 1.0 - r };
831: for (i = 0; i < nverts; ++i) {
832: const PetscReal* vertices = coordinates + i * 3;
833: jacobian[0] += dNi_dxi[i] * vertices[0];
834: jacobian[2] += dNi_dxi[i] * vertices[1];
835: jacobian[1] += dNi_deta[i] * vertices[0];
836: jacobian[3] += dNi_deta[i] * vertices[1];
837: }
838: } else if (nverts == 3) { /* Linear triangle */
839: const PetscReal x2 = coordinates[2 * 3 + 0], y2 = coordinates[2 * 3 + 1];
841: /* Jacobian is constant */
842: jacobian[0] = (coordinates[0 * 3 + 0] - x2); jacobian[1] = (coordinates[1 * 3 + 0] - x2);
843: jacobian[2] = (coordinates[0 * 3 + 1] - y2); jacobian[3] = (coordinates[1 * 3 + 1] - y2);
844: } else SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "The number of 2-D entity vertices are invalid. Currently only support QUAD4 and TRI3 basis evaluations in 2-D : %D", nverts);
846: /* invert the jacobian */
847: if (ijacobian) {
848: DMatrix_Invert_2x2_Internal(jacobian, ijacobian, &volume);
849: }
850: } else {
852: if (nverts == 8) { /* Linear Hexahedra */
853: const PetscReal &r = quad[0];
854: const PetscReal &s = quad[1];
855: const PetscReal &t = quad[2];
856: const PetscReal dNi_dxi[8] = {- ( 1.0 - s) * ( 1.0 - t),
857: ( 1.0 - s) * ( 1.0 - t),
858: ( s) * ( 1.0 - t),
859: - ( s) * ( 1.0 - t),
860: - ( 1.0 - s) * ( t),
861: ( 1.0 - s) * ( t),
862: ( s) * ( t),
863: - ( s) * ( t)
864: };
866: const PetscReal dNi_deta[8] = { - ( 1.0 - r) * ( 1.0 - t),
867: - ( r) * ( 1.0 - t),
868: ( r) * ( 1.0 - t),
869: ( 1.0 - r) * ( 1.0 - t),
870: - ( 1.0 - r) * ( t),
871: - ( r) * ( t),
872: ( r) * ( t),
873: ( 1.0 - r) * ( t)
874: };
876: const PetscReal dNi_dzeta[8] = { - ( 1.0 - r) * ( 1.0 - s),
877: - ( r) * ( 1.0 - s),
878: - ( r) * ( s),
879: - ( 1.0 - r) * ( s),
880: ( 1.0 - r) * ( 1.0 - s),
881: ( r) * ( 1.0 - s),
882: ( r) * ( s),
883: ( 1.0 - r) * ( s)
884: };
886: for (i = 0; i < nverts; ++i) {
887: const PetscReal* vertex = coordinates + i * 3;
888: jacobian[0] += dNi_dxi[i] * vertex[0];
889: jacobian[3] += dNi_dxi[i] * vertex[1];
890: jacobian[6] += dNi_dxi[i] * vertex[2];
891: jacobian[1] += dNi_deta[i] * vertex[0];
892: jacobian[4] += dNi_deta[i] * vertex[1];
893: jacobian[7] += dNi_deta[i] * vertex[2];
894: jacobian[2] += dNi_dzeta[i] * vertex[0];
895: jacobian[5] += dNi_dzeta[i] * vertex[1];
896: jacobian[8] += dNi_dzeta[i] * vertex[2];
897: }
898: } else if (nverts == 4) { /* Linear Tetrahedra */
899: const PetscReal x0 = coordinates[/*0 * 3 +*/ 0], y0 = coordinates[/*0 * 3 +*/ 1], z0 = coordinates[/*0 * 3 +*/ 2];
901: /* compute the jacobian */
902: jacobian[0] = coordinates[1 * 3 + 0] - x0; jacobian[1] = coordinates[2 * 3 + 0] - x0; jacobian[2] = coordinates[3 * 3 + 0] - x0;
903: jacobian[3] = coordinates[1 * 3 + 1] - y0; jacobian[4] = coordinates[2 * 3 + 1] - y0; jacobian[5] = coordinates[3 * 3 + 1] - y0;
904: jacobian[6] = coordinates[1 * 3 + 2] - z0; jacobian[7] = coordinates[2 * 3 + 2] - z0; jacobian[8] = coordinates[3 * 3 + 2] - z0;
905: } else SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "The number of 3-D entity vertices are invalid. Currently only support HEX8 and TET4 basis evaluations in 3-D : %D", nverts);
907: if (ijacobian) {
908: /* invert the jacobian */
909: DMatrix_Invert_3x3_Internal(jacobian, ijacobian, &volume);
910: }
912: }
913: if (volume < 1e-12) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Element has zero volume: %g. Degenerate element or invalid connectivity\n", volume);
914: if (dvolume) *dvolume = volume;
915: return(0);
916: }
919: PetscErrorCode FEMComputeBasis_JandF(const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal *quadrature, PetscReal *phypts,
920: PetscReal *phibasis, PetscReal *jacobian, PetscReal *ijacobian, PetscReal* volume)
921: {
925: switch (dim) {
926: case 1:
927: Compute_Lagrange_Basis_1D_Internal(nverts, coordinates, 1, quadrature, phypts,
928: NULL, phibasis, NULL, jacobian, ijacobian, volume);
929: break;
930: case 2:
931: Compute_Lagrange_Basis_2D_Internal(nverts, coordinates, 1, quadrature, phypts,
932: NULL, phibasis, NULL, NULL, jacobian, ijacobian, volume);
933: break;
934: case 3:
935: Compute_Lagrange_Basis_3D_Internal(nverts, coordinates, 1, quadrature, phypts,
936: NULL, phibasis, NULL, NULL, NULL, jacobian, ijacobian, volume);
937: break;
938: default:
939: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %D", dim);
940: }
941: return(0);
942: }
944: /*@C
945: DMMoabPToRMapping - Compute the mapping from the physical coordinate system for a given element to the
946: canonical reference element. In addition to finding the inverse mapping evaluation through Newton iteration,
947: the basis function at the parametric point is also evaluated optionally.
949: Input Parameter:
950: + PetscInt dim - the element dimension (1=EDGE, 2=QUAD/TRI, 3=HEX/TET)
951: . PetscInt nverts - the number of vertices in the physical element
952: . PetscReal coordinates - the coordinates of vertices in the physical element
953: - PetscReal[3] xphy - the coordinates of physical point for which natural coordinates (in reference frame) are sought
955: Output Parameter:
956: + PetscReal[3] natparam - the natural coordinates (in reference frame) corresponding to xphy
957: - PetscReal[nverts] phi - the basis functions evaluated at the natural coordinates (natparam)
959: Level: advanced
961: @*/
962: PetscErrorCode DMMoabPToRMapping(const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal* xphy, PetscReal* natparam, PetscReal* phi)
963: {
964: /* Perform inverse evaluation for the mapping with use of Newton Raphson iteration */
965: const PetscReal tol = 1.0e-10;
966: const PetscInt max_iterations = 10;
967: const PetscReal error_tol_sqr = tol*tol;
968: PetscReal phibasis[8], jacobian[9], ijacobian[9], volume;
969: PetscReal phypts[3] = {0.0, 0.0, 0.0};
970: PetscReal delta[3] = {0.0, 0.0, 0.0};
971: PetscErrorCode ierr;
972: PetscInt iters=0;
973: PetscReal error=1.0;
980: PetscArrayzero(jacobian, dim * dim);
981: PetscArrayzero(ijacobian, dim * dim);
982: PetscArrayzero(phibasis, nverts);
984: /* zero initial guess */
985: natparam[0] = natparam[1] = natparam[2] = 0.0;
987: /* Compute delta = evaluate( xi) - x */
988: FEMComputeBasis_JandF(dim, nverts, coordinates, natparam, phypts, phibasis, jacobian, ijacobian, &volume);
990: error = 0.0;
991: switch(dim) {
992: case 3:
993: delta[2] = phypts[2] - xphy[2];
994: error += (delta[2]*delta[2]);
995: case 2:
996: delta[1] = phypts[1] - xphy[1];
997: error += (delta[1]*delta[1]);
998: case 1:
999: delta[0] = phypts[0] - xphy[0];
1000: error += (delta[0]*delta[0]);
1001: break;
1002: }
1004: while (error > error_tol_sqr) {
1006: if (++iters > max_iterations)
1007: 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]);
1009: /* Compute natparam -= J.inverse() * delta */
1010: switch(dim) {
1011: case 1:
1012: natparam[0] -= ijacobian[0] * delta[0];
1013: break;
1014: case 2:
1015: natparam[0] -= ijacobian[0] * delta[0] + ijacobian[1] * delta[1];
1016: natparam[1] -= ijacobian[2] * delta[0] + ijacobian[3] * delta[1];
1017: break;
1018: case 3:
1019: natparam[0] -= ijacobian[0] * delta[0] + ijacobian[1] * delta[1] + ijacobian[2] * delta[2];
1020: natparam[1] -= ijacobian[3] * delta[0] + ijacobian[4] * delta[1] + ijacobian[5] * delta[2];
1021: natparam[2] -= ijacobian[6] * delta[0] + ijacobian[7] * delta[1] + ijacobian[8] * delta[2];
1022: break;
1023: }
1025: /* Compute delta = evaluate( xi) - x */
1026: FEMComputeBasis_JandF(dim, nverts, coordinates, natparam, phypts, phibasis, jacobian, ijacobian, &volume);
1028: error = 0.0;
1029: switch(dim) {
1030: case 3:
1031: delta[2] = phypts[2] - xphy[2];
1032: error += (delta[2]*delta[2]);
1033: case 2:
1034: delta[1] = phypts[1] - xphy[1];
1035: error += (delta[1]*delta[1]);
1036: case 1:
1037: delta[0] = phypts[0] - xphy[0];
1038: error += (delta[0]*delta[0]);
1039: break;
1040: }
1041: }
1042: if (phi) {
1043: PetscArraycpy(phi, phibasis, nverts);
1044: }
1045: return(0);
1046: }