Actual source code: dmmbfem.cxx
1: #include <petscconf.h>
2: #include <petscdt.h>
3: #include <petsc/private/dmmbimpl.h>
5: /* Utility functions */
6: static inline PetscReal DMatrix_Determinant_2x2_Internal(const PetscReal inmat[2 * 2])
7: {
8: return inmat[0] * inmat[3] - inmat[1] * inmat[2];
9: }
11: static inline PetscErrorCode DMatrix_Invert_2x2_Internal(const PetscReal *inmat, PetscReal *outmat, PetscReal *determinant)
12: {
13: PetscReal det = DMatrix_Determinant_2x2_Internal(inmat);
14: if (outmat) {
15: outmat[0] = inmat[3] / det;
16: outmat[1] = -inmat[1] / det;
17: outmat[2] = -inmat[2] / det;
18: outmat[3] = inmat[0] / det;
19: }
20: if (determinant) *determinant = det;
21: return PETSC_SUCCESS;
22: }
24: static inline PetscReal DMatrix_Determinant_3x3_Internal(const PetscReal inmat[3 * 3])
25: {
26: 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]);
27: }
29: static inline PetscErrorCode DMatrix_Invert_3x3_Internal(const PetscReal *inmat, PetscReal *outmat, PetscReal *determinant)
30: {
31: PetscReal det = DMatrix_Determinant_3x3_Internal(inmat);
32: if (outmat) {
33: outmat[0] = (inmat[8] * inmat[4] - inmat[7] * inmat[5]) / det;
34: outmat[1] = -(inmat[8] * inmat[1] - inmat[7] * inmat[2]) / det;
35: outmat[2] = (inmat[5] * inmat[1] - inmat[4] * inmat[2]) / det;
36: outmat[3] = -(inmat[8] * inmat[3] - inmat[6] * inmat[5]) / det;
37: outmat[4] = (inmat[8] * inmat[0] - inmat[6] * inmat[2]) / det;
38: outmat[5] = -(inmat[5] * inmat[0] - inmat[3] * inmat[2]) / det;
39: outmat[6] = (inmat[7] * inmat[3] - inmat[6] * inmat[4]) / det;
40: outmat[7] = -(inmat[7] * inmat[0] - inmat[6] * inmat[1]) / det;
41: outmat[8] = (inmat[4] * inmat[0] - inmat[3] * inmat[1]) / det;
42: }
43: if (determinant) *determinant = det;
44: return PETSC_SUCCESS;
45: }
47: static inline PetscReal DMatrix_Determinant_4x4_Internal(PetscReal inmat[4 * 4])
48: {
49: 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]));
50: }
52: static inline PetscErrorCode DMatrix_Invert_4x4_Internal(PetscReal *inmat, PetscReal *outmat, PetscScalar *determinant)
53: {
54: PetscReal det = DMatrix_Determinant_4x4_Internal(inmat);
55: if (outmat) {
56: 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;
57: 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;
58: 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;
59: 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;
60: 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;
61: 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;
62: 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;
63: 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;
64: 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;
65: 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;
66: 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;
67: 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;
68: 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;
69: 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;
70: 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;
71: 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;
72: }
73: if (determinant) *determinant = det;
74: return PETSC_SUCCESS;
75: }
77: /*
78: Compute_Lagrange_Basis_1D_Internal - Evaluate bases and derivatives at quadrature points for a EDGE2 or EDGE3 element.
80: The routine is given the coordinates of the vertices of a linear or quadratic edge element.
82: The routine evaluates the basis functions associated with each quadrature point provided,
83: and their derivatives with respect to X.
85: Notes:
87: Example Physical Element
88: .vb
89: 1-------2 1----3----2
90: EDGE2 EDGE3
91: .ve
93: Input Parameters:
94: + PetscInt nverts - the number of element vertices
95: . PetscReal coords[3*nverts] - the physical coordinates of the vertices (in canonical numbering)
96: . PetscInt npts - the number of evaluation points (quadrature points)
97: - PetscReal quad[3*npts] - the evaluation points (quadrature points) in the reference space
99: Output Parameters:
100: + PetscReal phypts[3*npts] - the evaluation points (quadrature points) transformed to the physical space
101: . PetscReal jxw[npts] - the jacobian determinant * quadrature weight necessary for assembling discrete contributions
102: . PetscReal phi[npts] - the bases evaluated at the specified quadrature points
103: . PetscReal dphidx[npts] - the derivative of the bases wrt X-direction evaluated at the specified quadrature points
104: . jacobian - jacobian
105: . ijacobian - ijacobian
106: - volume - volume
108: Level: advanced
109: */
110: static 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)
111: {
112: int i, j;
114: PetscFunctionBegin;
115: PetscAssertPointer(jacobian, 9);
116: PetscAssertPointer(ijacobian, 10);
117: PetscAssertPointer(volume, 11);
118: if (phypts) PetscCall(PetscArrayzero(phypts, npts * 3));
119: if (dphidx) { /* Reset arrays. */
120: PetscCall(PetscArrayzero(dphidx, npts * nverts));
121: }
122: if (nverts == 2) { /* Linear Edge */
124: for (j = 0; j < npts; j++) {
125: const PetscInt offset = j * nverts;
126: const PetscReal r = quad[j];
128: phi[0 + offset] = (1.0 - r);
129: phi[1 + offset] = (r);
131: const PetscReal dNi_dxi[2] = {-1.0, 1.0};
133: jacobian[0] = ijacobian[0] = volume[0] = 0.0;
134: for (i = 0; i < nverts; ++i) {
135: const PetscReal *vertices = coords + i * 3;
136: jacobian[0] += dNi_dxi[i] * vertices[0];
137: if (phypts) phypts[3 * j + 0] += phi[i + offset] * vertices[0];
138: }
140: /* invert the jacobian */
141: *volume = jacobian[0];
142: ijacobian[0] = 1.0 / jacobian[0];
143: jxw[j] *= *volume;
145: /* Divide by element jacobian. */
146: for (i = 0; i < nverts; i++) {
147: if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[0];
148: }
149: }
150: } else if (nverts == 3) { /* Quadratic Edge */
152: for (j = 0; j < npts; j++) {
153: const PetscInt offset = j * nverts;
154: const PetscReal r = quad[j];
156: phi[0 + offset] = 1.0 + r * (2.0 * r - 3.0);
157: phi[1 + offset] = 4.0 * r * (1.0 - r);
158: phi[2 + offset] = r * (2.0 * r - 1.0);
160: const PetscReal dNi_dxi[3] = {4 * r - 3.0, 4 * (1.0 - 2.0 * r), 4.0 * r - 1.0};
162: jacobian[0] = ijacobian[0] = volume[0] = 0.0;
163: for (i = 0; i < nverts; ++i) {
164: const PetscReal *vertices = coords + i * 3;
165: jacobian[0] += dNi_dxi[i] * vertices[0];
166: if (phypts) phypts[3 * j + 0] += phi[i + offset] * vertices[0];
167: }
169: /* invert the jacobian */
170: *volume = jacobian[0];
171: ijacobian[0] = 1.0 / jacobian[0];
172: if (jxw) jxw[j] *= *volume;
174: /* Divide by element jacobian. */
175: for (i = 0; i < nverts; i++) {
176: if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[0];
177: }
178: }
179: } 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);
180: PetscFunctionReturn(PETSC_SUCCESS);
181: }
183: /*
184: Compute_Lagrange_Basis_2D_Internal - Evaluate bases and derivatives at quadrature points for a QUAD4 or TRI3 element.
186: The routine is given the coordinates of the vertices of a quadrangle or triangle.
188: The routine evaluates the basis functions associated with each quadrature point provided,
189: and their derivatives with respect to X and Y.
191: Notes:
193: Example Physical Element (QUAD4)
194: .vb
195: 4------3 s
196: | | |
197: | | |
198: | | |
199: 1------2 0-------r
200: .ve
202: Input Parameters:
203: + PetscInt nverts - the number of element vertices
204: . PetscReal coords[3*nverts] - the physical coordinates of the vertices (in canonical numbering)
205: . PetscInt npts - the number of evaluation points (quadrature points)
206: - PetscReal quad[3*npts] - the evaluation points (quadrature points) in the reference space
208: Output Parameters:
209: + PetscReal phypts[3*npts] - the evaluation points (quadrature points) transformed to the physical space
210: . PetscReal jxw[npts] - the jacobian determinant * quadrature weight necessary for assembling discrete contributions
211: . PetscReal phi[npts] - the bases evaluated at the specified quadrature points
212: . PetscReal dphidx[npts] - the derivative of the bases wrt X-direction evaluated at the specified quadrature points
213: . PetscReal dphidy[npts] - the derivative of the bases wrt Y-direction evaluated at the specified quadrature points
214: . jacobian - jacobian
215: . ijacobian - ijacobian
216: - volume - volume
218: Level: advanced
219: */
220: static 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)
221: {
222: PetscInt i, j, k;
224: PetscFunctionBegin;
225: PetscAssertPointer(jacobian, 10);
226: PetscAssertPointer(ijacobian, 11);
227: PetscAssertPointer(volume, 12);
228: PetscCall(PetscArrayzero(phi, npts));
229: if (phypts) PetscCall(PetscArrayzero(phypts, npts * 3));
230: if (dphidx) { /* Reset arrays. */
231: PetscCall(PetscArrayzero(dphidx, npts * nverts));
232: PetscCall(PetscArrayzero(dphidy, npts * nverts));
233: }
234: if (nverts == 4) { /* Linear Quadrangle */
236: for (j = 0; j < npts; j++) {
237: const PetscInt offset = j * nverts;
238: const PetscReal r = quad[0 + j * 2];
239: const PetscReal s = quad[1 + j * 2];
241: phi[0 + offset] = (1.0 - r) * (1.0 - s);
242: phi[1 + offset] = r * (1.0 - s);
243: phi[2 + offset] = r * s;
244: phi[3 + offset] = (1.0 - r) * s;
246: const PetscReal dNi_dxi[4] = {-1.0 + s, 1.0 - s, s, -s};
247: const PetscReal dNi_deta[4] = {-1.0 + r, -r, r, 1.0 - r};
249: PetscCall(PetscArrayzero(jacobian, 4));
250: PetscCall(PetscArrayzero(ijacobian, 4));
251: for (i = 0; i < nverts; ++i) {
252: const PetscReal *vertices = coords + i * 3;
253: jacobian[0] += dNi_dxi[i] * vertices[0];
254: jacobian[2] += dNi_dxi[i] * vertices[1];
255: jacobian[1] += dNi_deta[i] * vertices[0];
256: jacobian[3] += dNi_deta[i] * vertices[1];
257: if (phypts) {
258: phypts[3 * j + 0] += phi[i + offset] * vertices[0];
259: phypts[3 * j + 1] += phi[i + offset] * vertices[1];
260: phypts[3 * j + 2] += phi[i + offset] * vertices[2];
261: }
262: }
264: /* invert the jacobian */
265: PetscCall(DMatrix_Invert_2x2_Internal(jacobian, ijacobian, volume));
266: PetscCheck(*volume >= 1e-12, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Quadrangular element has zero volume: %g. Degenerate element or invalid connectivity", *volume);
268: if (jxw) jxw[j] *= *volume;
270: /* Let us compute the bases derivatives by scaling with inverse jacobian. */
271: if (dphidx) {
272: for (i = 0; i < nverts; i++) {
273: for (k = 0; k < 2; ++k) {
274: if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[k * 2 + 0];
275: if (dphidy) dphidy[i + offset] += dNi_deta[i] * ijacobian[k * 2 + 1];
276: } /* for k=[0..2) */
277: } /* for i=[0..nverts) */
278: } /* if (dphidx) */
279: }
280: } else if (nverts == 3) { /* Linear triangle */
281: const PetscReal x2 = coords[2 * 3 + 0], y2 = coords[2 * 3 + 1];
283: PetscCall(PetscArrayzero(jacobian, 4));
284: PetscCall(PetscArrayzero(ijacobian, 4));
286: /* Jacobian is constant */
287: jacobian[0] = (coords[0 * 3 + 0] - x2);
288: jacobian[1] = (coords[1 * 3 + 0] - x2);
289: jacobian[2] = (coords[0 * 3 + 1] - y2);
290: jacobian[3] = (coords[1 * 3 + 1] - y2);
292: /* invert the jacobian */
293: PetscCall(DMatrix_Invert_2x2_Internal(jacobian, ijacobian, volume));
294: PetscCheck(*volume >= 1e-12, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Triangular element has zero volume: %g. Degenerate element or invalid connectivity", (double)*volume);
296: const PetscReal Dx[3] = {ijacobian[0], ijacobian[2], -ijacobian[0] - ijacobian[2]};
297: const PetscReal Dy[3] = {ijacobian[1], ijacobian[3], -ijacobian[1] - ijacobian[3]};
299: for (j = 0; j < npts; j++) {
300: const PetscInt offset = j * nverts;
301: const PetscReal r = quad[0 + j * 2];
302: const PetscReal s = quad[1 + j * 2];
304: if (jxw) jxw[j] *= 0.5;
306: /* const PetscReal Ni[3] = { r, s, 1.0 - r - s }; */
307: const PetscReal phipts_x = coords[2 * 3 + 0] + jacobian[0] * r + jacobian[1] * s;
308: const PetscReal phipts_y = coords[2 * 3 + 1] + jacobian[2] * r + jacobian[3] * s;
309: if (phypts) {
310: phypts[offset + 0] = phipts_x;
311: phypts[offset + 1] = phipts_y;
312: }
314: /* \phi_0 = (b.y - c.y) x + (b.x - c.x) y + c.x b.y - b.x c.y */
315: phi[0 + offset] = (ijacobian[0] * (phipts_x - x2) + ijacobian[1] * (phipts_y - y2));
316: /* \phi_1 = (c.y - a.y) x + (a.x - c.x) y + c.x a.y - a.x c.y */
317: phi[1 + offset] = (ijacobian[2] * (phipts_x - x2) + ijacobian[3] * (phipts_y - y2));
318: phi[2 + offset] = 1.0 - phi[0 + offset] - phi[1 + offset];
320: if (dphidx) {
321: dphidx[0 + offset] = Dx[0];
322: dphidx[1 + offset] = Dx[1];
323: dphidx[2 + offset] = Dx[2];
324: }
326: if (dphidy) {
327: dphidy[0 + offset] = Dy[0];
328: dphidy[1 + offset] = Dy[1];
329: dphidy[2 + offset] = Dy[2];
330: }
331: }
332: } 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);
333: PetscFunctionReturn(PETSC_SUCCESS);
334: }
336: /*
337: Compute_Lagrange_Basis_3D_Internal - Evaluate bases and derivatives at quadrature points for a HEX8 or TET4 element.
339: The routine is given the coordinates of the vertices of a hexahedra or tetrahedra.
341: The routine evaluates the basis functions associated with each quadrature point provided,
342: and their derivatives with respect to X, Y and Z.
344: Notes:
346: Example Physical Element (HEX8)
347: .vb
348: 8------7
349: /| /| t s
350: 5------6 | | /
351: | | | | |/
352: | 4----|-3 0-------r
353: |/ |/
354: 1------2
355: .ve
357: Input Parameters:
358: + PetscInt nverts - the number of element vertices
359: . PetscReal coords[3*nverts] - the physical coordinates of the vertices (in canonical numbering)
360: . PetscInt npts - the number of evaluation points (quadrature points)
361: - PetscReal quad[3*npts] - the evaluation points (quadrature points) in the reference space
363: Output Parameters:
364: + PetscReal phypts[3*npts] - the evaluation points (quadrature points) transformed to the physical space
365: . PetscReal jxw[npts] - the jacobian determinant * quadrature weight necessary for assembling discrete contributions
366: . PetscReal phi[npts] - the bases evaluated at the specified quadrature points
367: . PetscReal dphidx[npts] - the derivative of the bases wrt X-direction evaluated at the specified quadrature points
368: . PetscReal dphidy[npts] - the derivative of the bases wrt Y-direction evaluated at the specified quadrature points
369: . PetscReal dphidz[npts] - the derivative of the bases wrt Z-direction evaluated at the specified quadrature points
370: . jacobian - jacobian
371: . ijacobian - ijacobian
372: - volume - volume
374: Level: advanced
375: */
376: static 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)
377: {
378: PetscInt i, j, k;
380: PetscFunctionBegin;
381: PetscAssertPointer(jacobian, 11);
382: PetscAssertPointer(ijacobian, 12);
383: PetscAssertPointer(volume, 13);
385: PetscCall(PetscArrayzero(phi, npts));
386: if (phypts) PetscCall(PetscArrayzero(phypts, npts * 3));
387: if (dphidx) {
388: PetscCall(PetscArrayzero(dphidx, npts * nverts));
389: PetscCall(PetscArrayzero(dphidy, npts * nverts));
390: PetscCall(PetscArrayzero(dphidz, npts * nverts));
391: }
393: if (nverts == 8) { /* Linear Hexahedra */
395: for (j = 0; j < npts; j++) {
396: const PetscInt offset = j * nverts;
397: const PetscReal &r = quad[j * 3 + 0];
398: const PetscReal &s = quad[j * 3 + 1];
399: const PetscReal &t = quad[j * 3 + 2];
401: phi[offset + 0] = (1.0 - r) * (1.0 - s) * (1.0 - t); /* 0,0,0 */
402: phi[offset + 1] = (r) * (1.0 - s) * (1.0 - t); /* 1,0,0 */
403: phi[offset + 2] = (r) * (s) * (1.0 - t); /* 1,1,0 */
404: phi[offset + 3] = (1.0 - r) * (s) * (1.0 - t); /* 0,1,0 */
405: phi[offset + 4] = (1.0 - r) * (1.0 - s) * (t); /* 0,0,1 */
406: phi[offset + 5] = (r) * (1.0 - s) * (t); /* 1,0,1 */
407: phi[offset + 6] = (r) * (s) * (t); /* 1,1,1 */
408: phi[offset + 7] = (1.0 - r) * (s) * (t); /* 0,1,1 */
410: 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)};
412: 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)};
414: 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)};
416: PetscCall(PetscArrayzero(jacobian, 9));
417: PetscCall(PetscArrayzero(ijacobian, 9));
418: for (i = 0; i < nverts; ++i) {
419: const PetscReal *vertex = coords + i * 3;
420: jacobian[0] += dNi_dxi[i] * vertex[0];
421: jacobian[3] += dNi_dxi[i] * vertex[1];
422: jacobian[6] += dNi_dxi[i] * vertex[2];
423: jacobian[1] += dNi_deta[i] * vertex[0];
424: jacobian[4] += dNi_deta[i] * vertex[1];
425: jacobian[7] += dNi_deta[i] * vertex[2];
426: jacobian[2] += dNi_dzeta[i] * vertex[0];
427: jacobian[5] += dNi_dzeta[i] * vertex[1];
428: jacobian[8] += dNi_dzeta[i] * vertex[2];
429: if (phypts) {
430: phypts[3 * j + 0] += phi[i + offset] * vertex[0];
431: phypts[3 * j + 1] += phi[i + offset] * vertex[1];
432: phypts[3 * j + 2] += phi[i + offset] * vertex[2];
433: }
434: }
436: /* invert the jacobian */
437: PetscCall(DMatrix_Invert_3x3_Internal(jacobian, ijacobian, volume));
438: PetscCheck(*volume >= 1e-8, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Hexahedral element has zero volume: %g. Degenerate element or invalid connectivity", *volume);
440: if (jxw) jxw[j] *= (*volume);
442: /* Divide by element jacobian. */
443: for (i = 0; i < nverts; ++i) {
444: for (k = 0; k < 3; ++k) {
445: if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[0 * 3 + k];
446: if (dphidy) dphidy[i + offset] += dNi_deta[i] * ijacobian[1 * 3 + k];
447: if (dphidz) dphidz[i + offset] += dNi_dzeta[i] * ijacobian[2 * 3 + k];
448: }
449: }
450: }
451: } else if (nverts == 4) { /* Linear Tetrahedra */
452: PetscReal Dx[4] = {0, 0, 0, 0}, Dy[4] = {0, 0, 0, 0}, Dz[4] = {0, 0, 0, 0};
453: const PetscReal x0 = coords[/*0 * 3 +*/ 0], y0 = coords[/*0 * 3 +*/ 1], z0 = coords[/*0 * 3 +*/ 2];
455: PetscCall(PetscArrayzero(jacobian, 9));
456: PetscCall(PetscArrayzero(ijacobian, 9));
458: /* compute the jacobian */
459: jacobian[0] = coords[1 * 3 + 0] - x0;
460: jacobian[1] = coords[2 * 3 + 0] - x0;
461: jacobian[2] = coords[3 * 3 + 0] - x0;
462: jacobian[3] = coords[1 * 3 + 1] - y0;
463: jacobian[4] = coords[2 * 3 + 1] - y0;
464: jacobian[5] = coords[3 * 3 + 1] - y0;
465: jacobian[6] = coords[1 * 3 + 2] - z0;
466: jacobian[7] = coords[2 * 3 + 2] - z0;
467: jacobian[8] = coords[3 * 3 + 2] - z0;
469: /* invert the jacobian */
470: PetscCall(DMatrix_Invert_3x3_Internal(jacobian, ijacobian, volume));
471: PetscCheck(*volume >= 1e-8, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Tetrahedral element has zero volume: %g. Degenerate element or invalid connectivity", (double)*volume);
473: if (dphidx) {
474: 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;
475: 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;
476: 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;
477: Dx[3] = -(Dx[0] + Dx[1] + Dx[2]);
478: 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;
479: 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;
480: 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;
481: Dy[3] = -(Dy[0] + Dy[1] + Dy[2]);
482: 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;
483: 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;
484: 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;
485: Dz[3] = -(Dz[0] + Dz[1] + Dz[2]);
486: }
488: for (j = 0; j < npts; j++) {
489: const PetscInt offset = j * nverts;
490: const PetscReal &r = quad[j * 3 + 0];
491: const PetscReal &s = quad[j * 3 + 1];
492: const PetscReal &t = quad[j * 3 + 2];
494: if (jxw) jxw[j] *= *volume;
496: phi[offset + 0] = 1.0 - r - s - t;
497: phi[offset + 1] = r;
498: phi[offset + 2] = s;
499: phi[offset + 3] = t;
501: if (phypts) {
502: for (i = 0; i < nverts; ++i) {
503: const PetscReal *vertices = coords + i * 3;
504: phypts[3 * j + 0] += phi[i + offset] * vertices[0];
505: phypts[3 * j + 1] += phi[i + offset] * vertices[1];
506: phypts[3 * j + 2] += phi[i + offset] * vertices[2];
507: }
508: }
510: /* Now set the derivatives */
511: if (dphidx) {
512: dphidx[0 + offset] = Dx[0];
513: dphidx[1 + offset] = Dx[1];
514: dphidx[2 + offset] = Dx[2];
515: dphidx[3 + offset] = Dx[3];
517: dphidy[0 + offset] = Dy[0];
518: dphidy[1 + offset] = Dy[1];
519: dphidy[2 + offset] = Dy[2];
520: dphidy[3 + offset] = Dy[3];
522: dphidz[0 + offset] = Dz[0];
523: dphidz[1 + offset] = Dz[1];
524: dphidz[2 + offset] = Dz[2];
525: dphidz[3 + offset] = Dz[3];
526: }
528: } /* Tetrahedra -- ends */
529: } 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);
530: PetscFunctionReturn(PETSC_SUCCESS);
531: }
533: /*@C
534: DMMoabFEMComputeBasis - Evaluate bases and derivatives at quadrature points for a linear EDGE/QUAD/TRI/HEX/TET element.
536: The routine takes the coordinates of the vertices of an element and computes basis functions associated with
537: each quadrature point provided, and their derivatives with respect to X, Y and Z as appropriate.
539: Input Parameters:
540: + dim - the dimension
541: . nverts - the number of element vertices
542: . coordinates - the physical coordinates of the vertices (in canonical numbering)
543: - quadrature - the evaluation points (quadrature points) in the reference space
545: Output Parameters:
546: + phypts - the evaluation points (quadrature points) transformed to the physical space
547: . jacobian_quadrature_weight_product - the jacobian determinant * quadrature weight necessary for assembling discrete contributions
548: . fe_basis - the bases values evaluated at the specified quadrature points
549: - fe_basis_derivatives - the derivative of the bases wrt (X,Y,Z)-directions (depending on the dimension) evaluated at the specified quadrature points
551: Level: advanced
553: .seealso: `DMMoabCreate()`
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;
562: PetscFunctionBegin;
563: PetscAssertPointer(coordinates, 3);
565: PetscAssertPointer(fe_basis, 7);
566: compute_der = (fe_basis_derivatives != NULL);
568: /* Get the quadrature points and weights for the given quadrature rule */
569: PetscCall(PetscQuadratureGetData(quadrature, &idim, NULL, &npoints, &quadpts, &quadwts));
570: PetscCheck(idim == dim, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Dimension mismatch: provided (%" PetscInt_FMT ") vs quadrature (%" PetscInt_FMT ")", idim, dim);
571: if (jacobian_quadrature_weight_product) PetscCall(PetscArraycpy(jacobian_quadrature_weight_product, quadwts, npoints));
573: switch (dim) {
574: case 1:
575: 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));
576: break;
577: case 2:
578: 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));
579: break;
580: case 3:
581: 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));
582: break;
583: default:
584: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %" PetscInt_FMT, dim);
585: }
586: PetscFunctionReturn(PETSC_SUCCESS);
587: }
589: /*@C
590: DMMoabFEMCreateQuadratureDefault - Create default quadrature rules for integration over an element with a given
591: dimension and polynomial order (deciphered from number of element vertices).
593: Input Parameters:
594: + dim - the element dimension (1=EDGE, 2=QUAD/TRI, 3=HEX/TET)
595: - nverts - the number of vertices in the physical element
597: Output Parameter:
598: . quadrature - the quadrature object with default settings to integrate polynomials defined over the element
600: Level: advanced
602: .seealso: `DMMoabCreate()`
603: @*/
604: PetscErrorCode DMMoabFEMCreateQuadratureDefault(const PetscInt dim, const PetscInt nverts, PetscQuadrature *quadrature)
605: {
606: PetscReal *w, *x;
607: PetscInt nc = 1;
609: PetscFunctionBegin;
610: /* Create an appropriate quadrature rule to sample basis */
611: switch (dim) {
612: case 1:
613: /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */
614: PetscCall(PetscDTStroudConicalQuadrature(1, nc, nverts, 0, 1.0, quadrature));
615: break;
616: case 2:
617: /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */
618: if (nverts == 3) {
619: const PetscInt order = 2;
620: const PetscInt npoints = (order == 2 ? 3 : 6);
621: PetscCall(PetscMalloc2(npoints * 2, &x, npoints, &w));
622: if (npoints == 3) {
623: x[0] = x[1] = x[2] = x[5] = 1.0 / 6.0;
624: x[3] = x[4] = 2.0 / 3.0;
625: w[0] = w[1] = w[2] = 1.0 / 3.0;
626: } else if (npoints == 6) {
627: x[0] = x[1] = x[2] = 0.44594849091597;
628: x[3] = x[4] = 0.10810301816807;
629: x[5] = 0.44594849091597;
630: x[6] = x[7] = x[8] = 0.09157621350977;
631: x[9] = x[10] = 0.81684757298046;
632: x[11] = 0.09157621350977;
633: w[0] = w[1] = w[2] = 0.22338158967801;
634: w[3] = w[4] = w[5] = 0.10995174365532;
635: } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Triangle quadrature rules for points 3 and 6 supported; npoints : %" PetscInt_FMT, npoints);
636: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, quadrature));
637: PetscCall(PetscQuadratureSetOrder(*quadrature, order));
638: PetscCall(PetscQuadratureSetData(*quadrature, dim, nc, npoints, x, w));
639: /* PetscCall(PetscDTStroudConicalQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature)); */
640: } else PetscCall(PetscDTGaussTensorQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature));
641: break;
642: case 3:
643: /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */
644: if (nverts == 4) {
645: PetscInt order;
646: const PetscInt npoints = 4; // Choose between 4 and 10
647: PetscCall(PetscMalloc2(npoints * 3, &x, npoints, &w));
648: if (npoints == 4) { /* KEAST1, K1, N=4, O=4 */
649: const PetscReal x_4[12] = {0.5854101966249685, 0.1381966011250105, 0.1381966011250105, 0.1381966011250105, 0.1381966011250105, 0.1381966011250105,
650: 0.1381966011250105, 0.1381966011250105, 0.5854101966249685, 0.1381966011250105, 0.5854101966249685, 0.1381966011250105};
651: PetscCall(PetscArraycpy(x, x_4, 12));
653: w[0] = w[1] = w[2] = w[3] = 1.0 / 24.0;
654: order = 4;
655: } else if (npoints == 10) { /* KEAST3, K3 N=10, O=10 */
656: 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,
657: 0.5684305841968444, 0.1438564719343852, 0.0000000000000000, 0.5000000000000000, 0.5000000000000000, 0.5000000000000000, 0.0000000000000000, 0.5000000000000000, 0.5000000000000000, 0.5000000000000000,
658: 0.0000000000000000, 0.5000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.5000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.5000000000000000};
659: PetscCall(PetscArraycpy(x, x_10, 30));
661: w[0] = w[1] = w[2] = w[3] = 0.2177650698804054;
662: w[4] = w[5] = w[6] = w[7] = w[8] = w[9] = 0.0214899534130631;
663: order = 10;
664: } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Tetrahedral quadrature rules for points 4 and 10 supported; npoints : %" PetscInt_FMT, npoints);
665: PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, quadrature));
666: PetscCall(PetscQuadratureSetOrder(*quadrature, order));
667: PetscCall(PetscQuadratureSetData(*quadrature, dim, nc, npoints, x, w));
668: /* PetscCall(PetscDTStroudConicalQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature)); */
669: } else PetscCall(PetscDTGaussTensorQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature));
670: break;
671: default:
672: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %" PetscInt_FMT, dim);
673: }
674: PetscFunctionReturn(PETSC_SUCCESS);
675: }
677: /* Compute Jacobians */
678: static 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)
679: {
680: PetscFunctionBegin;
681: switch (dim) {
682: case 1:
683: PetscCall(Compute_Lagrange_Basis_1D_Internal(nverts, coordinates, 1, quadrature, phypts, NULL, phibasis, NULL, jacobian, ijacobian, volume));
684: break;
685: case 2:
686: PetscCall(Compute_Lagrange_Basis_2D_Internal(nverts, coordinates, 1, quadrature, phypts, NULL, phibasis, NULL, NULL, jacobian, ijacobian, volume));
687: break;
688: case 3:
689: PetscCall(Compute_Lagrange_Basis_3D_Internal(nverts, coordinates, 1, quadrature, phypts, NULL, phibasis, NULL, NULL, NULL, jacobian, ijacobian, volume));
690: break;
691: default:
692: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %" PetscInt_FMT, dim);
693: }
694: PetscFunctionReturn(PETSC_SUCCESS);
695: }
697: /*@C
698: DMMoabPToRMapping - Compute the mapping from the physical coordinate system for a given element to the
699: canonical reference element.
701: Input Parameters:
702: + dim - the element dimension (1=EDGE, 2=QUAD/TRI, 3=HEX/TET)
703: . nverts - the number of vertices in the physical element
704: . coordinates - the coordinates of vertices in the physical element
705: - xphy - the coordinates of physical point for which natural coordinates (in reference frame) are sought
707: Output Parameters:
708: + natparam - the natural coordinates (in reference frame) corresponding to xphy
709: - phi - the basis functions evaluated at the natural coordinates (natparam)
711: Level: advanced
713: Notes:
714: In addition to finding the inverse mapping evaluation through Newton iteration, the basis
715: function at the parametric point is also evaluated optionally.
717: .seealso: `DMMoabCreate()`
718: @*/
719: PetscErrorCode DMMoabPToRMapping(const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal *xphy, PetscReal *natparam, PetscReal *phi)
720: {
721: /* Perform inverse evaluation for the mapping with use of Newton Raphson iteration */
722: const PetscReal tol = 1.0e-10;
723: const PetscInt max_iterations = 10;
724: const PetscReal error_tol_sqr = tol * tol;
725: PetscReal phibasis[8], jacobian[9], ijacobian[9], volume;
726: PetscReal phypts[3] = {0.0, 0.0, 0.0};
727: PetscReal delta[3] = {0.0, 0.0, 0.0};
728: PetscInt iters = 0;
729: PetscReal error = 1.0;
731: PetscFunctionBegin;
732: PetscAssertPointer(coordinates, 3);
733: PetscAssertPointer(xphy, 4);
734: PetscAssertPointer(natparam, 5);
736: PetscCall(PetscArrayzero(jacobian, dim * dim));
737: PetscCall(PetscArrayzero(ijacobian, dim * dim));
738: PetscCall(PetscArrayzero(phibasis, nverts));
740: /* zero initial guess */
741: natparam[0] = natparam[1] = natparam[2] = 0.0;
743: /* Compute delta = evaluate( xi) - x */
744: PetscCall(FEMComputeBasis_JandF(dim, nverts, coordinates, natparam, phypts, phibasis, jacobian, ijacobian, &volume));
746: error = 0.0;
747: switch (dim) {
748: case 3:
749: delta[2] = phypts[2] - xphy[2];
750: error += (delta[2] * delta[2]);
751: case 2:
752: delta[1] = phypts[1] - xphy[1];
753: error += (delta[1] * delta[1]);
754: case 1:
755: delta[0] = phypts[0] - xphy[0];
756: error += (delta[0] * delta[0]);
757: break;
758: }
760: while (error > error_tol_sqr) {
761: 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]);
763: /* Compute natparam -= J.inverse() * delta */
764: switch (dim) {
765: case 1:
766: natparam[0] -= ijacobian[0] * delta[0];
767: break;
768: case 2:
769: natparam[0] -= ijacobian[0] * delta[0] + ijacobian[1] * delta[1];
770: natparam[1] -= ijacobian[2] * delta[0] + ijacobian[3] * delta[1];
771: break;
772: case 3:
773: natparam[0] -= ijacobian[0] * delta[0] + ijacobian[1] * delta[1] + ijacobian[2] * delta[2];
774: natparam[1] -= ijacobian[3] * delta[0] + ijacobian[4] * delta[1] + ijacobian[5] * delta[2];
775: natparam[2] -= ijacobian[6] * delta[0] + ijacobian[7] * delta[1] + ijacobian[8] * delta[2];
776: break;
777: }
779: /* Compute delta = evaluate( xi) - x */
780: PetscCall(FEMComputeBasis_JandF(dim, nverts, coordinates, natparam, phypts, phibasis, jacobian, ijacobian, &volume));
782: error = 0.0;
783: switch (dim) {
784: case 3:
785: delta[2] = phypts[2] - xphy[2];
786: error += (delta[2] * delta[2]);
787: case 2:
788: delta[1] = phypts[1] - xphy[1];
789: error += (delta[1] * delta[1]);
790: case 1:
791: delta[0] = phypts[0] - xphy[0];
792: error += (delta[0] * delta[0]);
793: break;
794: }
795: }
796: if (phi) PetscCall(PetscArraycpy(phi, phibasis, nverts));
797: PetscFunctionReturn(PETSC_SUCCESS);
798: }