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 PETSC_SUCCESS;
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 PETSC_SUCCESS;
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 PETSC_SUCCESS;
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;
116: PetscFunctionBegin;
120: if (phypts) PetscCall(PetscArrayzero(phypts, npts * 3));
121: if (dphidx) { /* Reset arrays. */
122: PetscCall(PetscArrayzero(dphidx, npts * nverts));
123: }
124: if (nverts == 2) { /* Linear Edge */
126: for (j = 0; j < npts; j++) {
127: const PetscInt offset = j * nverts;
128: const PetscReal r = quad[j];
130: phi[0 + offset] = (1.0 - r);
131: phi[1 + offset] = (r);
133: const PetscReal dNi_dxi[2] = {-1.0, 1.0};
135: jacobian[0] = ijacobian[0] = volume[0] = 0.0;
136: for (i = 0; i < nverts; ++i) {
137: const PetscReal *vertices = coords + i * 3;
138: jacobian[0] += dNi_dxi[i] * vertices[0];
139: if (phypts) phypts[3 * j + 0] += phi[i + offset] * vertices[0];
140: }
142: /* invert the jacobian */
143: *volume = jacobian[0];
144: ijacobian[0] = 1.0 / jacobian[0];
145: jxw[j] *= *volume;
147: /* Divide by element jacobian. */
148: for (i = 0; i < nverts; i++) {
149: if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[0];
150: }
151: }
152: } else if (nverts == 3) { /* Quadratic Edge */
154: for (j = 0; j < npts; j++) {
155: const PetscInt offset = j * nverts;
156: const PetscReal r = quad[j];
158: phi[0 + offset] = 1.0 + r * (2.0 * r - 3.0);
159: phi[1 + offset] = 4.0 * r * (1.0 - r);
160: phi[2 + offset] = r * (2.0 * r - 1.0);
162: const PetscReal dNi_dxi[3] = {4 * r - 3.0, 4 * (1.0 - 2.0 * r), 4.0 * r - 1.0};
164: jacobian[0] = ijacobian[0] = volume[0] = 0.0;
165: for (i = 0; i < nverts; ++i) {
166: const PetscReal *vertices = coords + i * 3;
167: jacobian[0] += dNi_dxi[i] * vertices[0];
168: if (phypts) phypts[3 * j + 0] += phi[i + offset] * vertices[0];
169: }
171: /* invert the jacobian */
172: *volume = jacobian[0];
173: ijacobian[0] = 1.0 / jacobian[0];
174: if (jxw) jxw[j] *= *volume;
176: /* Divide by element jacobian. */
177: for (i = 0; i < nverts; i++) {
178: if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[0];
179: }
180: }
181: } 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);
182: PetscFunctionReturn(PETSC_SUCCESS);
183: }
185: /*@C
186: Compute_Lagrange_Basis_2D_Internal - Evaluate bases and derivatives at quadrature points for a QUAD4 or TRI3 element.
188: The routine is given the coordinates of the vertices of a quadrangle or triangle.
190: The routine evaluates the basis functions associated with each quadrature point provided,
191: and their derivatives with respect to X and Y.
193: Notes:
195: Example Physical Element (QUAD4)
196: .vb
197: 4------3 s
198: | | |
199: | | |
200: | | |
201: 1------2 0-------r
202: .ve
204: Input Parameters:
205: + PetscInt nverts - the number of element vertices
206: . PetscReal coords[3*nverts] - the physical coordinates of the vertices (in canonical numbering)
207: . PetscInt npts - the number of evaluation points (quadrature points)
208: - PetscReal quad[3*npts] - the evaluation points (quadrature points) in the reference space
210: Output Parameters:
211: + PetscReal phypts[3*npts] - the evaluation points (quadrature points) transformed to the physical space
212: . PetscReal jxw[npts] - the jacobian determinant * quadrature weight necessary for assembling discrete contributions
213: . PetscReal phi[npts] - the bases evaluated at the specified quadrature points
214: . PetscReal dphidx[npts] - the derivative of the bases wrt X-direction evaluated at the specified quadrature points
215: . PetscReal dphidy[npts] - the derivative of the bases wrt Y-direction evaluated at the specified quadrature points
216: . jacobian - jacobian
217: . ijacobian - ijacobian
218: - volume - volume
220: Level: advanced
222: @*/
223: 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)
224: {
225: PetscInt i, j, k;
227: PetscFunctionBegin;
231: PetscCall(PetscArrayzero(phi, npts));
232: if (phypts) PetscCall(PetscArrayzero(phypts, npts * 3));
233: if (dphidx) { /* Reset arrays. */
234: PetscCall(PetscArrayzero(dphidx, npts * nverts));
235: PetscCall(PetscArrayzero(dphidy, npts * nverts));
236: }
237: if (nverts == 4) { /* Linear Quadrangle */
239: for (j = 0; j < npts; j++) {
240: const PetscInt offset = j * nverts;
241: const PetscReal r = quad[0 + j * 2];
242: const PetscReal s = quad[1 + j * 2];
244: phi[0 + offset] = (1.0 - r) * (1.0 - s);
245: phi[1 + offset] = r * (1.0 - s);
246: phi[2 + offset] = r * s;
247: phi[3 + offset] = (1.0 - r) * s;
249: const PetscReal dNi_dxi[4] = {-1.0 + s, 1.0 - s, s, -s};
250: const PetscReal dNi_deta[4] = {-1.0 + r, -r, r, 1.0 - r};
252: PetscCall(PetscArrayzero(jacobian, 4));
253: PetscCall(PetscArrayzero(ijacobian, 4));
254: for (i = 0; i < nverts; ++i) {
255: const PetscReal *vertices = coords + i * 3;
256: jacobian[0] += dNi_dxi[i] * vertices[0];
257: jacobian[2] += dNi_dxi[i] * vertices[1];
258: jacobian[1] += dNi_deta[i] * vertices[0];
259: jacobian[3] += dNi_deta[i] * vertices[1];
260: if (phypts) {
261: phypts[3 * j + 0] += phi[i + offset] * vertices[0];
262: phypts[3 * j + 1] += phi[i + offset] * vertices[1];
263: phypts[3 * j + 2] += phi[i + offset] * vertices[2];
264: }
265: }
267: /* invert the jacobian */
268: PetscCall(DMatrix_Invert_2x2_Internal(jacobian, ijacobian, volume));
269: PetscCheck(*volume >= 1e-12, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Quadrangular element has zero volume: %g. Degenerate element or invalid connectivity", *volume);
271: if (jxw) jxw[j] *= *volume;
273: /* Let us compute the bases derivatives by scaling with inverse jacobian. */
274: if (dphidx) {
275: for (i = 0; i < nverts; i++) {
276: for (k = 0; k < 2; ++k) {
277: if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[k * 2 + 0];
278: if (dphidy) dphidy[i + offset] += dNi_deta[i] * ijacobian[k * 2 + 1];
279: } /* for k=[0..2) */
280: } /* for i=[0..nverts) */
281: } /* if (dphidx) */
282: }
283: } else if (nverts == 3) { /* Linear triangle */
284: const PetscReal x2 = coords[2 * 3 + 0], y2 = coords[2 * 3 + 1];
286: PetscCall(PetscArrayzero(jacobian, 4));
287: PetscCall(PetscArrayzero(ijacobian, 4));
289: /* Jacobian is constant */
290: jacobian[0] = (coords[0 * 3 + 0] - x2);
291: jacobian[1] = (coords[1 * 3 + 0] - x2);
292: jacobian[2] = (coords[0 * 3 + 1] - y2);
293: jacobian[3] = (coords[1 * 3 + 1] - y2);
295: /* invert the jacobian */
296: PetscCall(DMatrix_Invert_2x2_Internal(jacobian, ijacobian, volume));
297: PetscCheck(*volume >= 1e-12, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Triangular element has zero volume: %g. Degenerate element or invalid connectivity", (double)*volume);
299: const PetscReal Dx[3] = {ijacobian[0], ijacobian[2], -ijacobian[0] - ijacobian[2]};
300: const PetscReal Dy[3] = {ijacobian[1], ijacobian[3], -ijacobian[1] - ijacobian[3]};
302: for (j = 0; j < npts; j++) {
303: const PetscInt offset = j * nverts;
304: const PetscReal r = quad[0 + j * 2];
305: const PetscReal s = quad[1 + j * 2];
307: if (jxw) jxw[j] *= 0.5;
309: /* const PetscReal Ni[3] = { r, s, 1.0 - r - s }; */
310: const PetscReal phipts_x = coords[2 * 3 + 0] + jacobian[0] * r + jacobian[1] * s;
311: const PetscReal phipts_y = coords[2 * 3 + 1] + jacobian[2] * r + jacobian[3] * s;
312: if (phypts) {
313: phypts[offset + 0] = phipts_x;
314: phypts[offset + 1] = phipts_y;
315: }
317: /* \phi_0 = (b.y - c.y) x + (b.x - c.x) y + c.x b.y - b.x c.y */
318: phi[0 + offset] = (ijacobian[0] * (phipts_x - x2) + ijacobian[1] * (phipts_y - y2));
319: /* \phi_1 = (c.y - a.y) x + (a.x - c.x) y + c.x a.y - a.x c.y */
320: phi[1 + offset] = (ijacobian[2] * (phipts_x - x2) + ijacobian[3] * (phipts_y - y2));
321: phi[2 + offset] = 1.0 - phi[0 + offset] - phi[1 + offset];
323: if (dphidx) {
324: dphidx[0 + offset] = Dx[0];
325: dphidx[1 + offset] = Dx[1];
326: dphidx[2 + offset] = Dx[2];
327: }
329: if (dphidy) {
330: dphidy[0 + offset] = Dy[0];
331: dphidy[1 + offset] = Dy[1];
332: dphidy[2 + offset] = Dy[2];
333: }
334: }
335: } 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);
336: PetscFunctionReturn(PETSC_SUCCESS);
337: }
339: /*@C
340: Compute_Lagrange_Basis_3D_Internal - Evaluate bases and derivatives at quadrature points for a HEX8 or TET4 element.
342: The routine is given the coordinates of the vertices of a hexahedra or tetrahedra.
344: The routine evaluates the basis functions associated with each quadrature point provided,
345: and their derivatives with respect to X, Y and Z.
347: Notes:
349: Example Physical Element (HEX8)
350: .vb
351: 8------7
352: /| /| t s
353: 5------6 | | /
354: | | | | |/
355: | 4----|-3 0-------r
356: |/ |/
357: 1------2
358: .ve
360: Input Parameters:
361: + PetscInt nverts - the number of element vertices
362: . PetscReal coords[3*nverts] - the physical coordinates of the vertices (in canonical numbering)
363: . PetscInt npts - the number of evaluation points (quadrature points)
364: - PetscReal quad[3*npts] - the evaluation points (quadrature points) in the reference space
366: Output Parameters:
367: + PetscReal phypts[3*npts] - the evaluation points (quadrature points) transformed to the physical space
368: . PetscReal jxw[npts] - the jacobian determinant * quadrature weight necessary for assembling discrete contributions
369: . PetscReal phi[npts] - the bases evaluated at the specified quadrature points
370: . PetscReal dphidx[npts] - the derivative of the bases wrt X-direction evaluated at the specified quadrature points
371: . PetscReal dphidy[npts] - the derivative of the bases wrt Y-direction evaluated at the specified quadrature points
372: . PetscReal dphidz[npts] - the derivative of the bases wrt Z-direction evaluated at the specified quadrature points
373: . jacobian - jacobian
374: . ijacobian - ijacobian
375: - volume - volume
377: Level: advanced
379: @*/
380: 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)
381: {
382: PetscInt i, j, k;
384: PetscFunctionBegin;
389: PetscCall(PetscArrayzero(phi, npts));
390: if (phypts) PetscCall(PetscArrayzero(phypts, npts * 3));
391: if (dphidx) {
392: PetscCall(PetscArrayzero(dphidx, npts * nverts));
393: PetscCall(PetscArrayzero(dphidy, npts * nverts));
394: PetscCall(PetscArrayzero(dphidz, npts * nverts));
395: }
397: if (nverts == 8) { /* Linear Hexahedra */
399: for (j = 0; j < npts; j++) {
400: const PetscInt offset = j * nverts;
401: const PetscReal &r = quad[j * 3 + 0];
402: const PetscReal &s = quad[j * 3 + 1];
403: const PetscReal &t = quad[j * 3 + 2];
405: phi[offset + 0] = (1.0 - r) * (1.0 - s) * (1.0 - t); /* 0,0,0 */
406: phi[offset + 1] = (r) * (1.0 - s) * (1.0 - t); /* 1,0,0 */
407: phi[offset + 2] = (r) * (s) * (1.0 - t); /* 1,1,0 */
408: phi[offset + 3] = (1.0 - r) * (s) * (1.0 - t); /* 0,1,0 */
409: phi[offset + 4] = (1.0 - r) * (1.0 - s) * (t); /* 0,0,1 */
410: phi[offset + 5] = (r) * (1.0 - s) * (t); /* 1,0,1 */
411: phi[offset + 6] = (r) * (s) * (t); /* 1,1,1 */
412: phi[offset + 7] = (1.0 - r) * (s) * (t); /* 0,1,1 */
414: 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)};
416: 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)};
418: 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)};
420: PetscCall(PetscArrayzero(jacobian, 9));
421: PetscCall(PetscArrayzero(ijacobian, 9));
422: for (i = 0; i < nverts; ++i) {
423: const PetscReal *vertex = coords + i * 3;
424: jacobian[0] += dNi_dxi[i] * vertex[0];
425: jacobian[3] += dNi_dxi[i] * vertex[1];
426: jacobian[6] += dNi_dxi[i] * vertex[2];
427: jacobian[1] += dNi_deta[i] * vertex[0];
428: jacobian[4] += dNi_deta[i] * vertex[1];
429: jacobian[7] += dNi_deta[i] * vertex[2];
430: jacobian[2] += dNi_dzeta[i] * vertex[0];
431: jacobian[5] += dNi_dzeta[i] * vertex[1];
432: jacobian[8] += dNi_dzeta[i] * vertex[2];
433: if (phypts) {
434: phypts[3 * j + 0] += phi[i + offset] * vertex[0];
435: phypts[3 * j + 1] += phi[i + offset] * vertex[1];
436: phypts[3 * j + 2] += phi[i + offset] * vertex[2];
437: }
438: }
440: /* invert the jacobian */
441: PetscCall(DMatrix_Invert_3x3_Internal(jacobian, ijacobian, volume));
442: PetscCheck(*volume >= 1e-8, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Hexahedral element has zero volume: %g. Degenerate element or invalid connectivity", *volume);
444: if (jxw) jxw[j] *= (*volume);
446: /* Divide by element jacobian. */
447: for (i = 0; i < nverts; ++i) {
448: for (k = 0; k < 3; ++k) {
449: if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[0 * 3 + k];
450: if (dphidy) dphidy[i + offset] += dNi_deta[i] * ijacobian[1 * 3 + k];
451: if (dphidz) dphidz[i + offset] += dNi_dzeta[i] * ijacobian[2 * 3 + k];
452: }
453: }
454: }
455: } else if (nverts == 4) { /* Linear Tetrahedra */
456: PetscReal Dx[4] = {0, 0, 0, 0}, Dy[4] = {0, 0, 0, 0}, Dz[4] = {0, 0, 0, 0};
457: const PetscReal x0 = coords[/*0 * 3 +*/ 0], y0 = coords[/*0 * 3 +*/ 1], z0 = coords[/*0 * 3 +*/ 2];
459: PetscCall(PetscArrayzero(jacobian, 9));
460: PetscCall(PetscArrayzero(ijacobian, 9));
462: /* compute the jacobian */
463: jacobian[0] = coords[1 * 3 + 0] - x0;
464: jacobian[1] = coords[2 * 3 + 0] - x0;
465: jacobian[2] = coords[3 * 3 + 0] - x0;
466: jacobian[3] = coords[1 * 3 + 1] - y0;
467: jacobian[4] = coords[2 * 3 + 1] - y0;
468: jacobian[5] = coords[3 * 3 + 1] - y0;
469: jacobian[6] = coords[1 * 3 + 2] - z0;
470: jacobian[7] = coords[2 * 3 + 2] - z0;
471: jacobian[8] = coords[3 * 3 + 2] - z0;
473: /* invert the jacobian */
474: PetscCall(DMatrix_Invert_3x3_Internal(jacobian, ijacobian, volume));
475: PetscCheck(*volume >= 1e-8, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Tetrahedral element has zero volume: %g. Degenerate element or invalid connectivity", (double)*volume);
477: if (dphidx) {
478: 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;
479: 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;
480: 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;
481: Dx[3] = -(Dx[0] + Dx[1] + Dx[2]);
482: 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;
483: 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;
484: 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;
485: Dy[3] = -(Dy[0] + Dy[1] + Dy[2]);
486: 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;
487: 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;
488: 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;
489: Dz[3] = -(Dz[0] + Dz[1] + Dz[2]);
490: }
492: for (j = 0; j < npts; j++) {
493: const PetscInt offset = j * nverts;
494: const PetscReal &r = quad[j * 3 + 0];
495: const PetscReal &s = quad[j * 3 + 1];
496: const PetscReal &t = quad[j * 3 + 2];
498: if (jxw) jxw[j] *= *volume;
500: phi[offset + 0] = 1.0 - r - s - t;
501: phi[offset + 1] = r;
502: phi[offset + 2] = s;
503: phi[offset + 3] = t;
505: if (phypts) {
506: for (i = 0; i < nverts; ++i) {
507: const PetscScalar *vertices = coords + i * 3;
508: phypts[3 * j + 0] += phi[i + offset] * vertices[0];
509: phypts[3 * j + 1] += phi[i + offset] * vertices[1];
510: phypts[3 * j + 2] += phi[i + offset] * vertices[2];
511: }
512: }
514: /* Now set the derivatives */
515: if (dphidx) {
516: dphidx[0 + offset] = Dx[0];
517: dphidx[1 + offset] = Dx[1];
518: dphidx[2 + offset] = Dx[2];
519: dphidx[3 + offset] = Dx[3];
521: dphidy[0 + offset] = Dy[0];
522: dphidy[1 + offset] = Dy[1];
523: dphidy[2 + offset] = Dy[2];
524: dphidy[3 + offset] = Dy[3];
526: dphidz[0 + offset] = Dz[0];
527: dphidz[1 + offset] = Dz[1];
528: dphidz[2 + offset] = Dz[2];
529: dphidz[3 + offset] = Dz[3];
530: }
532: } /* Tetrahedra -- ends */
533: } 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);
534: PetscFunctionReturn(PETSC_SUCCESS);
535: }
537: /*@C
538: DMMoabFEMComputeBasis - Evaluate bases and derivatives at quadrature points for a linear EDGE/QUAD/TRI/HEX/TET element.
540: The routine takes the coordinates of the vertices of an element and computes basis functions associated with
541: each quadrature point provided, and their derivatives with respect to X, Y and Z as appropriate.
543: Input Parameters:
544: + PetscInt nverts - the number of element vertices
545: . PetscReal coords[3*nverts] - the physical coordinates of the vertices (in canonical numbering)
546: . PetscInt npts - the number of evaluation points (quadrature points)
547: - PetscReal quad[3*npts] - the evaluation points (quadrature points) in the reference space
549: Output Parameters:
550: + PetscReal phypts[3*npts] - the evaluation points (quadrature points) transformed to the physical space
551: . PetscReal jxw[npts] - the jacobian determinant * quadrature weight necessary for assembling discrete contributions
552: . PetscReal fe_basis[npts] - the bases values evaluated at the specified quadrature points
553: - 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
555: Level: advanced
557: @*/
558: 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)
559: {
560: PetscInt npoints, idim;
561: bool compute_der;
562: const PetscReal *quadpts, *quadwts;
563: PetscReal jacobian[9], ijacobian[9], volume;
565: PetscFunctionBegin;
569: compute_der = (fe_basis_derivatives != NULL);
571: /* Get the quadrature points and weights for the given quadrature rule */
572: PetscCall(PetscQuadratureGetData(quadrature, &idim, NULL, &npoints, &quadpts, &quadwts));
573: PetscCheck(idim == dim, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Dimension mismatch: provided (%" PetscInt_FMT ") vs quadrature (%" PetscInt_FMT ")", idim, dim);
574: if (jacobian_quadrature_weight_product) PetscCall(PetscArraycpy(jacobian_quadrature_weight_product, quadwts, npoints));
576: switch (dim) {
577: case 1:
578: PetscCall(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));
579: break;
580: case 2:
581: PetscCall(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));
582: break;
583: case 3:
584: PetscCall(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));
585: break;
586: default:
587: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %" PetscInt_FMT, dim);
588: }
589: PetscFunctionReturn(PETSC_SUCCESS);
590: }
592: /*@C
593: DMMoabFEMCreateQuadratureDefault - Create default quadrature rules for integration over an element with a given
594: dimension and polynomial order (deciphered from number of element vertices).
596: Input Parameters:
598: + PetscInt dim - the element dimension (1=EDGE, 2=QUAD/TRI, 3=HEX/TET)
599: - PetscInt nverts - the number of vertices in the physical element
601: Output Parameter:
602: . PetscQuadrature quadrature - the quadrature object with default settings to integrate polynomials defined over the element
604: Level: advanced
606: @*/
607: PetscErrorCode DMMoabFEMCreateQuadratureDefault(const PetscInt dim, const PetscInt nverts, PetscQuadrature *quadrature)
608: {
609: PetscReal *w, *x;
610: PetscInt nc = 1;
612: PetscFunctionBegin;
613: /* Create an appropriate quadrature rule to sample basis */
614: switch (dim) {
615: case 1:
616: /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */
617: PetscCall(PetscDTStroudConicalQuadrature(1, nc, nverts, 0, 1.0, quadrature));
618: break;
619: case 2:
620: /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */
621: if (nverts == 3) {
622: const PetscInt order = 2;
623: const PetscInt npoints = (order == 2 ? 3 : 6);
624: PetscCall(PetscMalloc2(npoints * 2, &x, npoints, &w));
625: if (npoints == 3) {
626: x[0] = x[1] = x[2] = x[5] = 1.0 / 6.0;
627: x[3] = x[4] = 2.0 / 3.0;
628: w[0] = w[1] = w[2] = 1.0 / 3.0;
629: } else if (npoints == 6) {
630: x[0] = x[1] = x[2] = 0.44594849091597;
631: x[3] = x[4] = 0.10810301816807;
632: x[5] = 0.44594849091597;
633: x[6] = x[7] = x[8] = 0.09157621350977;
634: x[9] = x[10] = 0.81684757298046;
635: x[11] = 0.09157621350977;
636: w[0] = w[1] = w[2] = 0.22338158967801;
637: w[3] = w[4] = w[5] = 0.10995174365532;
638: } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Triangle quadrature rules for points 3 and 6 supported; npoints : %" PetscInt_FMT, npoints);
639: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, quadrature));
640: PetscCall(PetscQuadratureSetOrder(*quadrature, order));
641: PetscCall(PetscQuadratureSetData(*quadrature, dim, nc, npoints, x, w));
642: /* PetscCall(PetscDTStroudConicalQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature)); */
643: } else PetscCall(PetscDTGaussTensorQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature));
644: break;
645: case 3:
646: /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */
647: if (nverts == 4) {
648: PetscInt order;
649: const PetscInt npoints = 4; // Choose between 4 and 10
650: PetscCall(PetscMalloc2(npoints * 3, &x, npoints, &w));
651: if (npoints == 4) { /* KEAST1, K1, N=4, O=4 */
652: const PetscReal x_4[12] = {0.5854101966249685, 0.1381966011250105, 0.1381966011250105, 0.1381966011250105, 0.1381966011250105, 0.1381966011250105,
653: 0.1381966011250105, 0.1381966011250105, 0.5854101966249685, 0.1381966011250105, 0.5854101966249685, 0.1381966011250105};
654: PetscCall(PetscArraycpy(x, x_4, 12));
656: w[0] = w[1] = w[2] = w[3] = 1.0 / 24.0;
657: order = 4;
658: } else if (npoints == 10) { /* KEAST3, K3 N=10, O=10 */
659: 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,
660: 0.5684305841968444, 0.1438564719343852, 0.0000000000000000, 0.5000000000000000, 0.5000000000000000, 0.5000000000000000, 0.0000000000000000, 0.5000000000000000, 0.5000000000000000, 0.5000000000000000,
661: 0.0000000000000000, 0.5000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.5000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.5000000000000000};
662: PetscCall(PetscArraycpy(x, x_10, 30));
664: w[0] = w[1] = w[2] = w[3] = 0.2177650698804054;
665: w[4] = w[5] = w[6] = w[7] = w[8] = w[9] = 0.0214899534130631;
666: order = 10;
667: } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Tetrahedral quadrature rules for points 4 and 10 supported; npoints : %" PetscInt_FMT, npoints);
668: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, quadrature));
669: PetscCall(PetscQuadratureSetOrder(*quadrature, order));
670: PetscCall(PetscQuadratureSetData(*quadrature, dim, nc, npoints, x, w));
671: /* PetscCall(PetscDTStroudConicalQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature)); */
672: } else PetscCall(PetscDTGaussTensorQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature));
673: break;
674: default:
675: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %" PetscInt_FMT, dim);
676: }
677: PetscFunctionReturn(PETSC_SUCCESS);
678: }
680: /* Compute Jacobians */
681: PetscErrorCode ComputeJacobian_Internal(const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal *quad, PetscReal *phypts, PetscReal *jacobian, PetscReal *ijacobian, PetscReal *dvolume)
682: {
683: PetscInt i;
684: PetscReal volume = 1.0;
686: PetscFunctionBegin;
690: PetscCall(PetscArrayzero(jacobian, dim * dim));
691: if (ijacobian) PetscCall(PetscArrayzero(ijacobian, dim * dim));
692: if (phypts) PetscCall(PetscArrayzero(phypts, /*npts=1 * */ 3));
694: if (dim == 1) {
695: const PetscReal &r = quad[0];
696: if (nverts == 2) { /* Linear Edge */
697: const PetscReal dNi_dxi[2] = {-1.0, 1.0};
699: for (i = 0; i < nverts; ++i) {
700: const PetscReal *vertices = coordinates + i * 3;
701: jacobian[0] += dNi_dxi[i] * vertices[0];
702: }
703: } else if (nverts == 3) { /* Quadratic Edge */
704: const PetscReal dNi_dxi[3] = {4 * r - 3.0, 4 * (1.0 - 2.0 * r), 4.0 * r - 1.0};
706: for (i = 0; i < nverts; ++i) {
707: const PetscReal *vertices = coordinates + i * 3;
708: jacobian[0] += dNi_dxi[i] * vertices[0];
709: }
710: } 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);
712: if (ijacobian) {
713: /* invert the jacobian */
714: ijacobian[0] = 1.0 / jacobian[0];
715: }
716: } else if (dim == 2) {
717: if (nverts == 4) { /* Linear Quadrangle */
718: const PetscReal &r = quad[0];
719: const PetscReal &s = quad[1];
721: const PetscReal dNi_dxi[4] = {-1.0 + s, 1.0 - s, s, -s};
722: const PetscReal dNi_deta[4] = {-1.0 + r, -r, r, 1.0 - r};
724: for (i = 0; i < nverts; ++i) {
725: const PetscReal *vertices = coordinates + i * 3;
726: jacobian[0] += dNi_dxi[i] * vertices[0];
727: jacobian[2] += dNi_dxi[i] * vertices[1];
728: jacobian[1] += dNi_deta[i] * vertices[0];
729: jacobian[3] += dNi_deta[i] * vertices[1];
730: }
731: } else if (nverts == 3) { /* Linear triangle */
732: const PetscReal x2 = coordinates[2 * 3 + 0], y2 = coordinates[2 * 3 + 1];
734: /* Jacobian is constant */
735: jacobian[0] = (coordinates[0 * 3 + 0] - x2);
736: jacobian[1] = (coordinates[1 * 3 + 0] - x2);
737: jacobian[2] = (coordinates[0 * 3 + 1] - y2);
738: jacobian[3] = (coordinates[1 * 3 + 1] - y2);
739: } 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);
741: /* invert the jacobian */
742: if (ijacobian) PetscCall(DMatrix_Invert_2x2_Internal(jacobian, ijacobian, &volume));
743: } else {
744: if (nverts == 8) { /* Linear Hexahedra */
745: const PetscReal &r = quad[0];
746: const PetscReal &s = quad[1];
747: const PetscReal &t = quad[2];
748: 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)};
750: 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)};
752: 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)};
754: for (i = 0; i < nverts; ++i) {
755: const PetscReal *vertex = coordinates + i * 3;
756: jacobian[0] += dNi_dxi[i] * vertex[0];
757: jacobian[3] += dNi_dxi[i] * vertex[1];
758: jacobian[6] += dNi_dxi[i] * vertex[2];
759: jacobian[1] += dNi_deta[i] * vertex[0];
760: jacobian[4] += dNi_deta[i] * vertex[1];
761: jacobian[7] += dNi_deta[i] * vertex[2];
762: jacobian[2] += dNi_dzeta[i] * vertex[0];
763: jacobian[5] += dNi_dzeta[i] * vertex[1];
764: jacobian[8] += dNi_dzeta[i] * vertex[2];
765: }
766: } else if (nverts == 4) { /* Linear Tetrahedra */
767: const PetscReal x0 = coordinates[/*0 * 3 +*/ 0], y0 = coordinates[/*0 * 3 +*/ 1], z0 = coordinates[/*0 * 3 +*/ 2];
769: /* compute the jacobian */
770: jacobian[0] = coordinates[1 * 3 + 0] - x0;
771: jacobian[1] = coordinates[2 * 3 + 0] - x0;
772: jacobian[2] = coordinates[3 * 3 + 0] - x0;
773: jacobian[3] = coordinates[1 * 3 + 1] - y0;
774: jacobian[4] = coordinates[2 * 3 + 1] - y0;
775: jacobian[5] = coordinates[3 * 3 + 1] - y0;
776: jacobian[6] = coordinates[1 * 3 + 2] - z0;
777: jacobian[7] = coordinates[2 * 3 + 2] - z0;
778: jacobian[8] = coordinates[3 * 3 + 2] - z0;
779: } 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);
781: if (ijacobian) {
782: /* invert the jacobian */
783: PetscCall(DMatrix_Invert_3x3_Internal(jacobian, ijacobian, &volume));
784: }
785: }
786: PetscCheck(volume >= 1e-12, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Element has zero volume: %g. Degenerate element or invalid connectivity", volume);
787: if (dvolume) *dvolume = volume;
788: PetscFunctionReturn(PETSC_SUCCESS);
789: }
791: 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)
792: {
793: PetscFunctionBegin;
794: switch (dim) {
795: case 1:
796: PetscCall(Compute_Lagrange_Basis_1D_Internal(nverts, coordinates, 1, quadrature, phypts, NULL, phibasis, NULL, jacobian, ijacobian, volume));
797: break;
798: case 2:
799: PetscCall(Compute_Lagrange_Basis_2D_Internal(nverts, coordinates, 1, quadrature, phypts, NULL, phibasis, NULL, NULL, jacobian, ijacobian, volume));
800: break;
801: case 3:
802: PetscCall(Compute_Lagrange_Basis_3D_Internal(nverts, coordinates, 1, quadrature, phypts, NULL, phibasis, NULL, NULL, NULL, jacobian, ijacobian, volume));
803: break;
804: default:
805: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %" PetscInt_FMT, dim);
806: }
807: PetscFunctionReturn(PETSC_SUCCESS);
808: }
810: /*@C
811: DMMoabPToRMapping - Compute the mapping from the physical coordinate system for a given element to the
812: canonical reference element. In addition to finding the inverse mapping evaluation through Newton iteration,
813: the basis function at the parametric point is also evaluated optionally.
815: Input Parameters:
816: + PetscInt dim - the element dimension (1=EDGE, 2=QUAD/TRI, 3=HEX/TET)
817: . PetscInt nverts - the number of vertices in the physical element
818: . PetscReal coordinates - the coordinates of vertices in the physical element
819: - PetscReal[3] xphy - the coordinates of physical point for which natural coordinates (in reference frame) are sought
821: Output Parameters:
822: + PetscReal[3] natparam - the natural coordinates (in reference frame) corresponding to xphy
823: - PetscReal[nverts] phi - the basis functions evaluated at the natural coordinates (natparam)
825: Level: advanced
827: @*/
828: PetscErrorCode DMMoabPToRMapping(const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal *xphy, PetscReal *natparam, PetscReal *phi)
829: {
830: /* Perform inverse evaluation for the mapping with use of Newton Raphson iteration */
831: const PetscReal tol = 1.0e-10;
832: const PetscInt max_iterations = 10;
833: const PetscReal error_tol_sqr = tol * tol;
834: PetscReal phibasis[8], jacobian[9], ijacobian[9], volume;
835: PetscReal phypts[3] = {0.0, 0.0, 0.0};
836: PetscReal delta[3] = {0.0, 0.0, 0.0};
837: PetscInt iters = 0;
838: PetscReal error = 1.0;
840: PetscFunctionBegin;
845: PetscCall(PetscArrayzero(jacobian, dim * dim));
846: PetscCall(PetscArrayzero(ijacobian, dim * dim));
847: PetscCall(PetscArrayzero(phibasis, nverts));
849: /* zero initial guess */
850: natparam[0] = natparam[1] = natparam[2] = 0.0;
852: /* Compute delta = evaluate( xi) - x */
853: PetscCall(FEMComputeBasis_JandF(dim, nverts, coordinates, natparam, phypts, phibasis, jacobian, ijacobian, &volume));
855: error = 0.0;
856: switch (dim) {
857: case 3:
858: delta[2] = phypts[2] - xphy[2];
859: error += (delta[2] * delta[2]);
860: case 2:
861: delta[1] = phypts[1] - xphy[1];
862: error += (delta[1] * delta[1]);
863: case 1:
864: delta[0] = phypts[0] - xphy[0];
865: error += (delta[0] * delta[0]);
866: break;
867: }
869: while (error > error_tol_sqr) {
870: PetscCheck(++iters <= max_iterations, 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]);
872: /* Compute natparam -= J.inverse() * delta */
873: switch (dim) {
874: case 1:
875: natparam[0] -= ijacobian[0] * delta[0];
876: break;
877: case 2:
878: natparam[0] -= ijacobian[0] * delta[0] + ijacobian[1] * delta[1];
879: natparam[1] -= ijacobian[2] * delta[0] + ijacobian[3] * delta[1];
880: break;
881: case 3:
882: natparam[0] -= ijacobian[0] * delta[0] + ijacobian[1] * delta[1] + ijacobian[2] * delta[2];
883: natparam[1] -= ijacobian[3] * delta[0] + ijacobian[4] * delta[1] + ijacobian[5] * delta[2];
884: natparam[2] -= ijacobian[6] * delta[0] + ijacobian[7] * delta[1] + ijacobian[8] * delta[2];
885: break;
886: }
888: /* Compute delta = evaluate( xi) - x */
889: PetscCall(FEMComputeBasis_JandF(dim, nverts, coordinates, natparam, phypts, phibasis, jacobian, ijacobian, &volume));
891: error = 0.0;
892: switch (dim) {
893: case 3:
894: delta[2] = phypts[2] - xphy[2];
895: error += (delta[2] * delta[2]);
896: case 2:
897: delta[1] = phypts[1] - xphy[1];
898: error += (delta[1] * delta[1]);
899: case 1:
900: delta[0] = phypts[0] - xphy[0];
901: error += (delta[0] * delta[0]);
902: break;
903: }
904: }
905: if (phi) PetscCall(PetscArraycpy(phi, phibasis, nverts));
906: PetscFunctionReturn(PETSC_SUCCESS);
907: }