Actual source code: dmmbfem.cxx
2: #include <petscconf.h>
3: #include <petscdt.h>
4: #include <petsc/private/dmmbimpl.h>
6: /* Utility functions */
7: static inline PetscReal DMatrix_Determinant_2x2_Internal (const PetscReal inmat[2 * 2])
8: {
9: return inmat[0] * inmat[3] - inmat[1] * inmat[2];
10: }
12: static inline PetscErrorCode DMatrix_Invert_2x2_Internal(const PetscReal *inmat, PetscReal *outmat, PetscReal *determinant)
13: {
14: PetscReal det = DMatrix_Determinant_2x2_Internal(inmat);
15: if (outmat) {
16: outmat[0] = inmat[3] / det;
17: outmat[1] = -inmat[1] / det;
18: outmat[2] = -inmat[2] / det;
19: outmat[3] = inmat[0] / det;
20: }
21: if (determinant) *determinant = det;
22: return 0;
23: }
25: static inline PetscReal DMatrix_Determinant_3x3_Internal(const PetscReal inmat[3 * 3])
26: {
27: return inmat[0] * (inmat[8] * inmat[4] - inmat[7] * inmat[5])
28: - inmat[3] * (inmat[8] * inmat[1] - inmat[7] * inmat[2])
29: + inmat[6] * (inmat[5] * inmat[1] - inmat[4] * inmat[2]);
30: }
32: static inline PetscErrorCode DMatrix_Invert_3x3_Internal (const PetscReal *inmat, PetscReal *outmat, PetscScalar *determinant)
33: {
34: PetscReal det = DMatrix_Determinant_3x3_Internal(inmat);
35: if (outmat) {
36: outmat[0] = (inmat[8] * inmat[4] - inmat[7] * inmat[5]) / det;
37: outmat[1] = -(inmat[8] * inmat[1] - inmat[7] * inmat[2]) / det;
38: outmat[2] = (inmat[5] * inmat[1] - inmat[4] * inmat[2]) / det;
39: outmat[3] = -(inmat[8] * inmat[3] - inmat[6] * inmat[5]) / det;
40: outmat[4] = (inmat[8] * inmat[0] - inmat[6] * inmat[2]) / det;
41: outmat[5] = -(inmat[5] * inmat[0] - inmat[3] * inmat[2]) / det;
42: outmat[6] = (inmat[7] * inmat[3] - inmat[6] * inmat[4]) / det;
43: outmat[7] = -(inmat[7] * inmat[0] - inmat[6] * inmat[1]) / det;
44: outmat[8] = (inmat[4] * inmat[0] - inmat[3] * inmat[1]) / det;
45: }
46: if (determinant) *determinant = det;
47: return 0;
48: }
50: inline PetscReal DMatrix_Determinant_4x4_Internal(PetscReal inmat[4 * 4])
51: {
52: return
53: inmat[0 + 0 * 4] * (
54: inmat[1 + 1 * 4] * (inmat[2 + 2 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 2 * 4])
55: - inmat[1 + 2 * 4] * (inmat[2 + 1 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 1 * 4])
56: + inmat[1 + 3 * 4] * (inmat[2 + 1 * 4] * inmat[3 + 2 * 4] - inmat[2 + 2 * 4] * inmat[3 + 1 * 4]))
57: - inmat[0 + 1 * 4] * (
58: inmat[1 + 0 * 4] * (inmat[2 + 2 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 2 * 4])
59: - inmat[1 + 2 * 4] * (inmat[2 + 0 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 0 * 4])
60: + inmat[1 + 3 * 4] * (inmat[2 + 0 * 4] * inmat[3 + 2 * 4] - inmat[2 + 2 * 4] * inmat[3 + 0 * 4]))
61: + inmat[0 + 2 * 4] * (
62: inmat[1 + 0 * 4] * (inmat[2 + 1 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 1 * 4])
63: - inmat[1 + 1 * 4] * (inmat[2 + 0 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 0 * 4])
64: + inmat[1 + 3 * 4] * (inmat[2 + 0 * 4] * inmat[3 + 1 * 4] - inmat[2 + 1 * 4] * inmat[3 + 0 * 4]))
65: - inmat[0 + 3 * 4] * (
66: inmat[1 + 0 * 4] * (inmat[2 + 1 * 4] * inmat[3 + 2 * 4] - inmat[2 + 2 * 4] * inmat[3 + 1 * 4])
67: - inmat[1 + 1 * 4] * (inmat[2 + 0 * 4] * inmat[3 + 2 * 4] - inmat[2 + 2 * 4] * inmat[3 + 0 * 4])
68: + inmat[1 + 2 * 4] * (inmat[2 + 0 * 4] * inmat[3 + 1 * 4] - inmat[2 + 1 * 4] * inmat[3 + 0 * 4]));
69: }
71: inline PetscErrorCode DMatrix_Invert_4x4_Internal (PetscReal *inmat, PetscReal *outmat, PetscScalar *determinant)
72: {
73: PetscReal det = DMatrix_Determinant_4x4_Internal(inmat);
74: if (outmat) {
75: outmat[0] = (inmat[5] * inmat[10] * inmat[15] + inmat[6] * inmat[11] * inmat[13] + inmat[7] * inmat[9] * inmat[14] - inmat[5] * inmat[11] * inmat[14] - inmat[6] * inmat[9] * inmat[15] - inmat[7] * inmat[10] * inmat[13]) / det;
76: outmat[1] = (inmat[1] * inmat[11] * inmat[14] + inmat[2] * inmat[9] * inmat[15] + inmat[3] * inmat[10] * inmat[13] - inmat[1] * inmat[10] * inmat[15] - inmat[2] * inmat[11] * inmat[13] - inmat[3] * inmat[9] * inmat[14]) / det;
77: outmat[2] = (inmat[1] * inmat[6] * inmat[15] + inmat[2] * inmat[7] * inmat[13] + inmat[3] * inmat[5] * inmat[14] - inmat[1] * inmat[7] * inmat[14] - inmat[2] * inmat[5] * inmat[15] - inmat[3] * inmat[6] * inmat[13]) / det;
78: outmat[3] = (inmat[1] * inmat[7] * inmat[10] + inmat[2] * inmat[5] * inmat[11] + inmat[3] * inmat[6] * inmat[9] - inmat[1] * inmat[6] * inmat[11] - inmat[2] * inmat[7] * inmat[9] - inmat[3] * inmat[5] * inmat[10]) / det;
79: outmat[4] = (inmat[4] * inmat[11] * inmat[14] + inmat[6] * inmat[8] * inmat[15] + inmat[7] * inmat[10] * inmat[12] - inmat[4] * inmat[10] * inmat[15] - inmat[6] * inmat[11] * inmat[12] - inmat[7] * inmat[8] * inmat[14]) / det;
80: outmat[5] = (inmat[0] * inmat[10] * inmat[15] + inmat[2] * inmat[11] * inmat[12] + inmat[3] * inmat[8] * inmat[14] - inmat[0] * inmat[11] * inmat[14] - inmat[2] * inmat[8] * inmat[15] - inmat[3] * inmat[10] * inmat[12]) / det;
81: outmat[6] = (inmat[0] * inmat[7] * inmat[14] + inmat[2] * inmat[4] * inmat[15] + inmat[3] * inmat[6] * inmat[12] - inmat[0] * inmat[6] * inmat[15] - inmat[2] * inmat[7] * inmat[12] - inmat[3] * inmat[4] * inmat[14]) / det;
82: outmat[7] = (inmat[0] * inmat[6] * inmat[11] + inmat[2] * inmat[7] * inmat[8] + inmat[3] * inmat[4] * inmat[10] - inmat[0] * inmat[7] * inmat[10] - inmat[2] * inmat[4] * inmat[11] - inmat[3] * inmat[6] * inmat[8]) / det;
83: outmat[8] = (inmat[4] * inmat[9] * inmat[15] + inmat[5] * inmat[11] * inmat[12] + inmat[7] * inmat[8] * inmat[13] - inmat[4] * inmat[11] * inmat[13] - inmat[5] * inmat[8] * inmat[15] - inmat[7] * inmat[9] * inmat[12]) / det;
84: outmat[9] = (inmat[0] * inmat[11] * inmat[13] + inmat[1] * inmat[8] * inmat[15] + inmat[3] * inmat[9] * inmat[12] - inmat[0] * inmat[9] * inmat[15] - inmat[1] * inmat[11] * inmat[12] - inmat[3] * inmat[8] * inmat[13]) / det;
85: outmat[10] = (inmat[0] * inmat[5] * inmat[15] + inmat[1] * inmat[7] * inmat[12] + inmat[3] * inmat[4] * inmat[13] - inmat[0] * inmat[7] * inmat[13] - inmat[1] * inmat[4] * inmat[15] - inmat[3] * inmat[5] * inmat[12]) / det;
86: outmat[11] = (inmat[0] * inmat[7] * inmat[9] + inmat[1] * inmat[4] * inmat[11] + inmat[3] * inmat[5] * inmat[8] - inmat[0] * inmat[5] * inmat[11] - inmat[1] * inmat[7] * inmat[8] - inmat[3] * inmat[4] * inmat[9]) / det;
87: outmat[12] = (inmat[4] * inmat[10] * inmat[13] + inmat[5] * inmat[8] * inmat[14] + inmat[6] * inmat[9] * inmat[12] - inmat[4] * inmat[9] * inmat[14] - inmat[5] * inmat[10] * inmat[12] - inmat[6] * inmat[8] * inmat[13]) / det;
88: outmat[13] = (inmat[0] * inmat[9] * inmat[14] + inmat[1] * inmat[10] * inmat[12] + inmat[2] * inmat[8] * inmat[13] - inmat[0] * inmat[10] * inmat[13] - inmat[1] * inmat[8] * inmat[14] - inmat[2] * inmat[9] * inmat[12]) / det;
89: outmat[14] = (inmat[0] * inmat[6] * inmat[13] + inmat[1] * inmat[4] * inmat[14] + inmat[2] * inmat[5] * inmat[12] - inmat[0] * inmat[5] * inmat[14] - inmat[1] * inmat[6] * inmat[12] - inmat[2] * inmat[4] * inmat[13]) / det;
90: outmat[15] = (inmat[0] * inmat[5] * inmat[10] + inmat[1] * inmat[6] * inmat[8] + inmat[2] * inmat[4] * inmat[9] - inmat[0] * inmat[6] * inmat[9] - inmat[1] * inmat[4] * inmat[10] - inmat[2] * inmat[5] * inmat[8]) / det;
91: }
92: if (determinant) *determinant = det;
93: return 0;
94: }
96: /*@C
97: Compute_Lagrange_Basis_1D_Internal - Evaluate bases and derivatives at quadrature points for a EDGE2 or EDGE3 element.
99: The routine is given the coordinates of the vertices of a linear or quadratic edge element.
101: The routine evaluates the basis functions associated with each quadrature point provided,
102: and their derivatives with respect to X.
104: Notes:
106: Example Physical Element
107: .vb
108: 1-------2 1----3----2
109: EDGE2 EDGE3
110: .ve
112: Input Parameters:
113: + PetscInt nverts - the number of element vertices
114: . PetscReal coords[3*nverts] - the physical coordinates of the vertices (in canonical numbering)
115: . PetscInt npts - the number of evaluation points (quadrature points)
116: - PetscReal quad[3*npts] - the evaluation points (quadrature points) in the reference space
118: Output Parameters:
119: + PetscReal phypts[3*npts] - the evaluation points (quadrature points) transformed to the physical space
120: . PetscReal jxw[npts] - the jacobian determinant * quadrature weight necessary for assembling discrete contributions
121: . PetscReal phi[npts] - the bases evaluated at the specified quadrature points
122: . PetscReal dphidx[npts] - the derivative of the bases wrt X-direction evaluated at the specified quadrature points
123: . jacobian -
124: . ijacobian -
125: - volume
127: Level: advanced
129: @*/
130: PetscErrorCode Compute_Lagrange_Basis_1D_Internal(const PetscInt nverts, const PetscReal *coords, const PetscInt npts, const PetscReal *quad, PetscReal *phypts,
131: PetscReal *jxw, PetscReal *phi, PetscReal *dphidx,
132: PetscReal *jacobian, PetscReal *ijacobian, PetscReal *volume)
133: {
134: int i, j;
135: PetscErrorCode ierr;
141: if (phypts) {
142: PetscArrayzero(phypts, npts * 3);
143: }
144: if (dphidx) { /* Reset arrays. */
145: PetscArrayzero(dphidx, npts * nverts);
146: }
147: if (nverts == 2) { /* Linear Edge */
149: for (j = 0; j < npts; j++) {
150: const PetscInt offset = j * nverts;
151: const PetscReal r = quad[j];
153: phi[0 + offset] = ( 1.0 - r);
154: phi[1 + offset] = ( r);
156: const PetscReal dNi_dxi[2] = { -1.0, 1.0 };
158: jacobian[0] = ijacobian[0] = volume[0] = 0.0;
159: for (i = 0; i < nverts; ++i) {
160: const PetscReal* vertices = coords + i * 3;
161: jacobian[0] += dNi_dxi[i] * vertices[0];
162: if (phypts) {
163: phypts[3 * j + 0] += phi[i + offset] * vertices[0];
164: }
165: }
167: /* invert the jacobian */
168: *volume = jacobian[0];
169: ijacobian[0] = 1.0 / jacobian[0];
170: jxw[j] *= *volume;
172: /* Divide by element jacobian. */
173: for (i = 0; i < nverts; i++) {
174: if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[0];
175: }
176: }
177: } else if (nverts == 3) { /* Quadratic Edge */
179: for (j = 0; j < npts; j++) {
180: const PetscInt offset = j * nverts;
181: const PetscReal r = quad[j];
183: phi[0 + offset] = 1.0 + r * ( 2.0 * r - 3.0);
184: phi[1 + offset] = 4.0 * r * ( 1.0 - r);
185: phi[2 + offset] = r * ( 2.0 * r - 1.0);
187: const PetscReal dNi_dxi[3] = { 4 * r - 3.0, 4 * ( 1.0 - 2.0 * r), 4.0 * r - 1.0};
189: jacobian[0] = ijacobian[0] = volume[0] = 0.0;
190: for (i = 0; i < nverts; ++i) {
191: const PetscReal* vertices = coords + i * 3;
192: jacobian[0] += dNi_dxi[i] * vertices[0];
193: if (phypts) {
194: phypts[3 * j + 0] += phi[i + offset] * vertices[0];
195: }
196: }
198: /* invert the jacobian */
199: *volume = jacobian[0];
200: ijacobian[0] = 1.0 / jacobian[0];
201: if (jxw) jxw[j] *= *volume;
203: /* Divide by element jacobian. */
204: for (i = 0; i < nverts; i++) {
205: if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[0];
206: }
207: }
208: } 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);
209: return(0);
210: }
212: /*@C
213: Compute_Lagrange_Basis_2D_Internal - Evaluate bases and derivatives at quadrature points for a QUAD4 or TRI3 element.
215: The routine is given the coordinates of the vertices of a quadrangle or triangle.
217: The routine evaluates the basis functions associated with each quadrature point provided,
218: and their derivatives with respect to X and Y.
220: Notes:
222: Example Physical Element (QUAD4)
223: .vb
224: 4------3 s
225: | | |
226: | | |
227: | | |
228: 1------2 0-------r
229: .ve
231: Input Parameters:
232: + PetscInt nverts - the number of element vertices
233: . PetscReal coords[3*nverts] - the physical coordinates of the vertices (in canonical numbering)
234: . PetscInt npts - the number of evaluation points (quadrature points)
235: - PetscReal quad[3*npts] - the evaluation points (quadrature points) in the reference space
237: Output Parameters:
238: + PetscReal phypts[3*npts] - the evaluation points (quadrature points) transformed to the physical space
239: . PetscReal jxw[npts] - the jacobian determinant * quadrature weight necessary for assembling discrete contributions
240: . PetscReal phi[npts] - the bases evaluated at the specified quadrature points
241: . PetscReal dphidx[npts] - the derivative of the bases wrt X-direction evaluated at the specified quadrature points
242: . PetscReal dphidy[npts] - the derivative of the bases wrt Y-direction evaluated at the specified quadrature points
243: . jacobian -
244: . ijacobian -
245: - volume
247: Level: advanced
249: @*/
250: PetscErrorCode Compute_Lagrange_Basis_2D_Internal(const PetscInt nverts, const PetscReal *coords, const PetscInt npts, const PetscReal *quad, PetscReal *phypts,
251: PetscReal *jxw, PetscReal *phi, PetscReal *dphidx, PetscReal *dphidy,
252: PetscReal *jacobian, PetscReal *ijacobian, PetscReal *volume)
253: {
254: PetscInt i, j, k;
261: PetscArrayzero(phi, npts);
262: if (phypts) {
263: PetscArrayzero(phypts, npts * 3);
264: }
265: if (dphidx) { /* Reset arrays. */
266: PetscArrayzero(dphidx, npts * nverts);
267: PetscArrayzero(dphidy, npts * nverts);
268: }
269: if (nverts == 4) { /* Linear Quadrangle */
271: for (j = 0; j < npts; j++)
272: {
273: const PetscInt offset = j * nverts;
274: const PetscReal r = quad[0 + j * 2];
275: const PetscReal s = quad[1 + j * 2];
277: phi[0 + offset] = ( 1.0 - r) * ( 1.0 - s);
278: phi[1 + offset] = r * ( 1.0 - s);
279: phi[2 + offset] = r * s;
280: phi[3 + offset] = ( 1.0 - r) * s;
282: const PetscReal dNi_dxi[4] = { -1.0 + s, 1.0 - s, s, -s };
283: const PetscReal dNi_deta[4] = { -1.0 + r, -r, r, 1.0 - r };
285: PetscArrayzero(jacobian, 4);
286: PetscArrayzero(ijacobian, 4);
287: for (i = 0; i < nverts; ++i) {
288: const PetscReal* vertices = coords + i * 3;
289: jacobian[0] += dNi_dxi[i] * vertices[0];
290: jacobian[2] += dNi_dxi[i] * vertices[1];
291: jacobian[1] += dNi_deta[i] * vertices[0];
292: jacobian[3] += dNi_deta[i] * vertices[1];
293: if (phypts) {
294: phypts[3 * j + 0] += phi[i + offset] * vertices[0];
295: phypts[3 * j + 1] += phi[i + offset] * vertices[1];
296: phypts[3 * j + 2] += phi[i + offset] * vertices[2];
297: }
298: }
300: /* invert the jacobian */
301: DMatrix_Invert_2x2_Internal(jacobian, ijacobian, volume);
302: 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);
304: if (jxw) jxw[j] *= *volume;
306: /* Let us compute the bases derivatives by scaling with inverse jacobian. */
307: if (dphidx) {
308: for (i = 0; i < nverts; i++) {
309: for (k = 0; k < 2; ++k) {
310: if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[k * 2 + 0];
311: if (dphidy) dphidy[i + offset] += dNi_deta[i] * ijacobian[k * 2 + 1];
312: } /* for k=[0..2) */
313: } /* for i=[0..nverts) */
314: } /* if (dphidx) */
315: }
316: } else if (nverts == 3) { /* Linear triangle */
317: const PetscReal x2 = coords[2 * 3 + 0], y2 = coords[2 * 3 + 1];
319: PetscArrayzero(jacobian, 4);
320: PetscArrayzero(ijacobian, 4);
322: /* Jacobian is constant */
323: jacobian[0] = (coords[0 * 3 + 0] - x2); jacobian[1] = (coords[1 * 3 + 0] - x2);
324: jacobian[2] = (coords[0 * 3 + 1] - y2); jacobian[3] = (coords[1 * 3 + 1] - y2);
326: /* invert the jacobian */
327: DMatrix_Invert_2x2_Internal(jacobian, ijacobian, volume);
328: 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);
330: const PetscReal Dx[3] = { ijacobian[0], ijacobian[2], - ijacobian[0] - ijacobian[2] };
331: const PetscReal Dy[3] = { ijacobian[1], ijacobian[3], - ijacobian[1] - ijacobian[3] };
333: for (j = 0; j < npts; j++) {
334: const PetscInt offset = j * nverts;
335: const PetscReal r = quad[0 + j * 2];
336: const PetscReal s = quad[1 + j * 2];
338: if (jxw) jxw[j] *= 0.5;
340: /* const PetscReal Ni[3] = { r, s, 1.0 - r - s }; */
341: const PetscReal phipts_x = coords[2 * 3 + 0] + jacobian[0] * r + jacobian[1] * s;
342: const PetscReal phipts_y = coords[2 * 3 + 1] + jacobian[2] * r + jacobian[3] * s;
343: if (phypts) {
344: phypts[offset + 0] = phipts_x;
345: phypts[offset + 1] = phipts_y;
346: }
348: /* \phi_0 = (b.y - c.y) x + (b.x - c.x) y + c.x b.y - b.x c.y */
349: phi[0 + offset] = ( ijacobian[0] * (phipts_x - x2) + ijacobian[1] * (phipts_y - y2));
350: /* \phi_1 = (c.y - a.y) x + (a.x - c.x) y + c.x a.y - a.x c.y */
351: phi[1 + offset] = ( ijacobian[2] * (phipts_x - x2) + ijacobian[3] * (phipts_y - y2));
352: phi[2 + offset] = 1.0 - phi[0 + offset] - phi[1 + offset];
354: if (dphidx) {
355: dphidx[0 + offset] = Dx[0];
356: dphidx[1 + offset] = Dx[1];
357: dphidx[2 + offset] = Dx[2];
358: }
360: if (dphidy) {
361: dphidy[0 + offset] = Dy[0];
362: dphidy[1 + offset] = Dy[1];
363: dphidy[2 + offset] = Dy[2];
364: }
366: }
367: } 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);
368: return(0);
369: }
371: /*@C
372: Compute_Lagrange_Basis_3D_Internal - Evaluate bases and derivatives at quadrature points for a HEX8 or TET4 element.
374: The routine is given the coordinates of the vertices of a hexahedra or tetrahedra.
376: The routine evaluates the basis functions associated with each quadrature point provided,
377: and their derivatives with respect to X, Y and Z.
379: Notes:
381: Example Physical Element (HEX8)
382: .vb
383: 8------7
384: /| /| t s
385: 5------6 | | /
386: | | | | |/
387: | 4----|-3 0-------r
388: |/ |/
389: 1------2
390: .ve
392: Input Parameters:
393: + PetscInt nverts - the number of element vertices
394: . PetscReal coords[3*nverts] - the physical coordinates of the vertices (in canonical numbering)
395: . PetscInt npts - the number of evaluation points (quadrature points)
396: - PetscReal quad[3*npts] - the evaluation points (quadrature points) in the reference space
398: Output Parameters:
399: + PetscReal phypts[3*npts] - the evaluation points (quadrature points) transformed to the physical space
400: . PetscReal jxw[npts] - the jacobian determinant * quadrature weight necessary for assembling discrete contributions
401: . PetscReal phi[npts] - the bases evaluated at the specified quadrature points
402: . PetscReal dphidx[npts] - the derivative of the bases wrt X-direction evaluated at the specified quadrature points
403: . PetscReal dphidy[npts] - the derivative of the bases wrt Y-direction evaluated at the specified quadrature points
404: . PetscReal dphidz[npts] - the derivative of the bases wrt Z-direction evaluated at the specified quadrature points
405: . jacobian -
406: . ijacobian -
407: - volume
409: Level: advanced
411: @*/
412: PetscErrorCode Compute_Lagrange_Basis_3D_Internal(const PetscInt nverts, const PetscReal *coords, const PetscInt npts, const PetscReal *quad, PetscReal *phypts,
413: PetscReal *jxw, PetscReal *phi, PetscReal *dphidx, PetscReal *dphidy, PetscReal *dphidz,
414: PetscReal *jacobian, PetscReal *ijacobian, PetscReal *volume)
415: {
416: PetscInt i, j, k;
424: PetscArrayzero(phi, npts);
425: if (phypts) {
426: PetscArrayzero(phypts, npts * 3);
427: }
428: if (dphidx) {
429: PetscArrayzero(dphidx, npts * nverts);
430: PetscArrayzero(dphidy, npts * nverts);
431: PetscArrayzero(dphidz, npts * nverts);
432: }
434: if (nverts == 8) { /* Linear Hexahedra */
436: for (j = 0; j < npts; j++) {
437: const PetscInt offset = j * nverts;
438: const PetscReal& r = quad[j * 3 + 0];
439: const PetscReal& s = quad[j * 3 + 1];
440: const PetscReal& t = quad[j * 3 + 2];
442: phi[offset + 0] = ( 1.0 - r) * ( 1.0 - s) * ( 1.0 - t); /* 0,0,0 */
443: phi[offset + 1] = ( r) * ( 1.0 - s) * ( 1.0 - t); /* 1,0,0 */
444: phi[offset + 2] = ( r) * ( s) * ( 1.0 - t); /* 1,1,0 */
445: phi[offset + 3] = ( 1.0 - r) * ( s) * ( 1.0 - t); /* 0,1,0 */
446: phi[offset + 4] = ( 1.0 - r) * ( 1.0 - s) * ( t); /* 0,0,1 */
447: phi[offset + 5] = ( r) * ( 1.0 - s) * ( t); /* 1,0,1 */
448: phi[offset + 6] = ( r) * ( s) * ( t); /* 1,1,1 */
449: phi[offset + 7] = ( 1.0 - r) * ( s) * ( t); /* 0,1,1 */
451: const PetscReal dNi_dxi[8] = {- ( 1.0 - s) * ( 1.0 - t),
452: ( 1.0 - s) * ( 1.0 - t),
453: ( s) * ( 1.0 - t),
454: - ( s) * ( 1.0 - t),
455: - ( 1.0 - s) * ( t),
456: ( 1.0 - s) * ( t),
457: ( s) * ( t),
458: - ( s) * ( t)
459: };
461: const PetscReal dNi_deta[8] = { - ( 1.0 - r) * ( 1.0 - t),
462: - ( r) * ( 1.0 - t),
463: ( r) * ( 1.0 - t),
464: ( 1.0 - r) * ( 1.0 - t),
465: - ( 1.0 - r) * ( t),
466: - ( r) * ( t),
467: ( r) * ( t),
468: ( 1.0 - r) * ( t)
469: };
471: const PetscReal dNi_dzeta[8] = { - ( 1.0 - r) * ( 1.0 - s),
472: - ( r) * ( 1.0 - s),
473: - ( r) * ( s),
474: - ( 1.0 - r) * ( s),
475: ( 1.0 - r) * ( 1.0 - s),
476: ( r) * ( 1.0 - s),
477: ( r) * ( s),
478: ( 1.0 - r) * ( s)
479: };
481: PetscArrayzero(jacobian, 9);
482: PetscArrayzero(ijacobian, 9);
483: for (i = 0; i < nverts; ++i) {
484: const PetscReal* vertex = coords + i * 3;
485: jacobian[0] += dNi_dxi[i] * vertex[0];
486: jacobian[3] += dNi_dxi[i] * vertex[1];
487: jacobian[6] += dNi_dxi[i] * vertex[2];
488: jacobian[1] += dNi_deta[i] * vertex[0];
489: jacobian[4] += dNi_deta[i] * vertex[1];
490: jacobian[7] += dNi_deta[i] * vertex[2];
491: jacobian[2] += dNi_dzeta[i] * vertex[0];
492: jacobian[5] += dNi_dzeta[i] * vertex[1];
493: jacobian[8] += dNi_dzeta[i] * vertex[2];
494: if (phypts) {
495: phypts[3 * j + 0] += phi[i + offset] * vertex[0];
496: phypts[3 * j + 1] += phi[i + offset] * vertex[1];
497: phypts[3 * j + 2] += phi[i + offset] * vertex[2];
498: }
499: }
501: /* invert the jacobian */
502: DMatrix_Invert_3x3_Internal(jacobian, ijacobian, volume);
503: 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);
505: if (jxw) jxw[j] *= (*volume);
507: /* Divide by element jacobian. */
508: for (i = 0; i < nverts; ++i) {
509: for (k = 0; k < 3; ++k) {
510: if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[0 * 3 + k];
511: if (dphidy) dphidy[i + offset] += dNi_deta[i] * ijacobian[1 * 3 + k];
512: if (dphidz) dphidz[i + offset] += dNi_dzeta[i] * ijacobian[2 * 3 + k];
513: }
514: }
515: }
516: } else if (nverts == 4) { /* Linear Tetrahedra */
517: PetscReal Dx[4]={0,0,0,0}, Dy[4]={0,0,0,0}, Dz[4]={0,0,0,0};
518: const PetscReal x0 = coords[/*0 * 3 +*/ 0], y0 = coords[/*0 * 3 +*/ 1], z0 = coords[/*0 * 3 +*/ 2];
520: PetscArrayzero(jacobian, 9);
521: PetscArrayzero(ijacobian, 9);
523: /* compute the jacobian */
524: jacobian[0] = coords[1 * 3 + 0] - x0; jacobian[1] = coords[2 * 3 + 0] - x0; jacobian[2] = coords[3 * 3 + 0] - x0;
525: jacobian[3] = coords[1 * 3 + 1] - y0; jacobian[4] = coords[2 * 3 + 1] - y0; jacobian[5] = coords[3 * 3 + 1] - y0;
526: jacobian[6] = coords[1 * 3 + 2] - z0; jacobian[7] = coords[2 * 3 + 2] - z0; jacobian[8] = coords[3 * 3 + 2] - z0;
528: /* invert the jacobian */
529: DMatrix_Invert_3x3_Internal(jacobian, ijacobian, volume);
530: 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);
532: if (dphidx) {
533: Dx[0] = ( coords[1 + 2 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3])
534: - coords[1 + 1 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3])
535: - coords[1 + 3 * 3] * ( coords[2 + 1 * 3] - coords[2 + 2 * 3])) / *volume;
536: Dx[1] = - ( coords[1 + 2 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3])
537: - coords[1 + 0 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3])
538: - coords[1 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 2 * 3])) / *volume;
539: Dx[2] = ( coords[1 + 1 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3])
540: - coords[1 + 0 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3])
541: - coords[1 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 1 * 3])) / *volume;
542: Dx[3] = - ( Dx[0] + Dx[1] + Dx[2]);
543: Dy[0] = ( coords[0 + 1 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3])
544: - coords[0 + 2 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3])
545: + coords[0 + 3 * 3] * ( coords[2 + 1 * 3] - coords[2 + 2 * 3])) / *volume;
546: Dy[1] = - ( coords[0 + 0 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3])
547: - coords[0 + 2 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3])
548: + coords[0 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 2 * 3])) / *volume;
549: Dy[2] = ( coords[0 + 0 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3])
550: - coords[0 + 1 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3])
551: + coords[0 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 1 * 3])) / *volume;
552: Dy[3] = - ( Dy[0] + Dy[1] + Dy[2]);
553: Dz[0] = ( coords[0 + 1 * 3] * (coords[1 + 3 * 3] - coords[1 + 2 * 3])
554: - coords[0 + 2 * 3] * (coords[1 + 3 * 3] - coords[1 + 1 * 3])
555: + coords[0 + 3 * 3] * (coords[1 + 2 * 3] - coords[1 + 1 * 3])) / *volume;
556: Dz[1] = - ( coords[0 + 0 * 3] * (coords[1 + 3 * 3] - coords[1 + 2 * 3])
557: + coords[0 + 2 * 3] * (coords[1 + 0 * 3] - coords[1 + 3 * 3])
558: - coords[0 + 3 * 3] * (coords[1 + 0 * 3] - coords[1 + 2 * 3])) / *volume;
559: Dz[2] = ( coords[0 + 0 * 3] * (coords[1 + 3 * 3] - coords[1 + 1 * 3])
560: + coords[0 + 1 * 3] * (coords[1 + 0 * 3] - coords[1 + 3 * 3])
561: - coords[0 + 3 * 3] * (coords[1 + 0 * 3] - coords[1 + 1 * 3])) / *volume;
562: Dz[3] = - ( Dz[0] + Dz[1] + Dz[2]);
563: }
565: for (j = 0; j < npts; j++) {
566: const PetscInt offset = j * nverts;
567: const PetscReal& r = quad[j * 3 + 0];
568: const PetscReal& s = quad[j * 3 + 1];
569: const PetscReal& t = quad[j * 3 + 2];
571: if (jxw) jxw[j] *= *volume;
573: phi[offset + 0] = 1.0 - r - s - t;
574: phi[offset + 1] = r;
575: phi[offset + 2] = s;
576: phi[offset + 3] = t;
578: if (phypts) {
579: for (i = 0; i < nverts; ++i) {
580: const PetscScalar* vertices = coords + i * 3;
581: phypts[3 * j + 0] += phi[i + offset] * vertices[0];
582: phypts[3 * j + 1] += phi[i + offset] * vertices[1];
583: phypts[3 * j + 2] += phi[i + offset] * vertices[2];
584: }
585: }
587: /* Now set the derivatives */
588: if (dphidx) {
589: dphidx[0 + offset] = Dx[0];
590: dphidx[1 + offset] = Dx[1];
591: dphidx[2 + offset] = Dx[2];
592: dphidx[3 + offset] = Dx[3];
594: dphidy[0 + offset] = Dy[0];
595: dphidy[1 + offset] = Dy[1];
596: dphidy[2 + offset] = Dy[2];
597: dphidy[3 + offset] = Dy[3];
599: dphidz[0 + offset] = Dz[0];
600: dphidz[1 + offset] = Dz[1];
601: dphidz[2 + offset] = Dz[2];
602: dphidz[3 + offset] = Dz[3];
603: }
605: } /* Tetrahedra -- ends */
606: } 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);
607: return(0);
608: }
610: /*@C
611: DMMoabFEMComputeBasis - Evaluate bases and derivatives at quadrature points for a linear EDGE/QUAD/TRI/HEX/TET element.
613: The routine takes the coordinates of the vertices of an element and computes basis functions associated with
614: each quadrature point provided, and their derivatives with respect to X, Y and Z as appropriate.
616: Input Parameters:
617: + PetscInt nverts, the number of element vertices
618: . PetscReal coords[3*nverts], the physical coordinates of the vertices (in canonical numbering)
619: . PetscInt npts, the number of evaluation points (quadrature points)
620: - PetscReal quad[3*npts], the evaluation points (quadrature points) in the reference space
622: Output Parameters:
623: + PetscReal phypts[3*npts], the evaluation points (quadrature points) transformed to the physical space
624: . PetscReal jxw[npts], the jacobian determinant * quadrature weight necessary for assembling discrete contributions
625: . PetscReal fe_basis[npts], the bases values evaluated at the specified quadrature points
626: - 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
628: Level: advanced
630: @*/
631: PetscErrorCode DMMoabFEMComputeBasis(const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscQuadrature quadrature,
632: PetscReal *phypts, PetscReal *jacobian_quadrature_weight_product,
633: PetscReal *fe_basis, PetscReal **fe_basis_derivatives)
634: {
635: PetscErrorCode ierr;
636: PetscInt npoints,idim;
637: bool compute_der;
638: const PetscReal *quadpts, *quadwts;
639: PetscReal jacobian[9], ijacobian[9], volume;
645: compute_der = (fe_basis_derivatives != NULL);
647: /* Get the quadrature points and weights for the given quadrature rule */
648: PetscQuadratureGetData(quadrature, &idim, NULL, &npoints, &quadpts, &quadwts);
649: if (idim != dim) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Dimension mismatch: provided (%D) vs quadrature (%D)\n",idim,dim);
650: if (jacobian_quadrature_weight_product) {
651: PetscArraycpy(jacobian_quadrature_weight_product, quadwts, npoints);
652: }
654: switch (dim) {
655: case 1:
656: Compute_Lagrange_Basis_1D_Internal(nverts, coordinates, npoints, quadpts, phypts,
657: jacobian_quadrature_weight_product, fe_basis,
658: (compute_der ? fe_basis_derivatives[0] : NULL),
659: jacobian, ijacobian, &volume);
660: break;
661: case 2:
662: Compute_Lagrange_Basis_2D_Internal(nverts, coordinates, npoints, quadpts, phypts,
663: jacobian_quadrature_weight_product, fe_basis,
664: (compute_der ? fe_basis_derivatives[0] : NULL),
665: (compute_der ? fe_basis_derivatives[1] : NULL),
666: jacobian, ijacobian, &volume);
667: break;
668: case 3:
669: Compute_Lagrange_Basis_3D_Internal(nverts, coordinates, npoints, quadpts, phypts,
670: jacobian_quadrature_weight_product, fe_basis,
671: (compute_der ? fe_basis_derivatives[0] : NULL),
672: (compute_der ? fe_basis_derivatives[1] : NULL),
673: (compute_der ? fe_basis_derivatives[2] : NULL),
674: jacobian, ijacobian, &volume);
675: break;
676: default:
677: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %D", dim);
678: }
679: return(0);
680: }
682: /*@C
683: DMMoabFEMCreateQuadratureDefault - Create default quadrature rules for integration over an element with a given
684: dimension and polynomial order (deciphered from number of element vertices).
686: Input Parameters:
688: + PetscInt dim - the element dimension (1=EDGE, 2=QUAD/TRI, 3=HEX/TET)
689: - PetscInt nverts - the number of vertices in the physical element
691: Output Parameter:
692: . PetscQuadrature quadrature - the quadrature object with default settings to integrate polynomials defined over the element
694: Level: advanced
696: @*/
697: PetscErrorCode DMMoabFEMCreateQuadratureDefault(const PetscInt dim, const PetscInt nverts, PetscQuadrature *quadrature)
698: {
699: PetscReal *w, *x;
700: PetscInt nc=1;
701: PetscErrorCode ierr;
704: /* Create an appropriate quadrature rule to sample basis */
705: switch (dim)
706: {
707: case 1:
708: /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */
709: PetscDTStroudConicalQuadrature(1, nc, nverts, 0, 1.0, quadrature);
710: break;
711: case 2:
712: /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */
713: if (nverts == 3) {
714: const PetscInt order = 2;
715: const PetscInt npoints = (order == 2 ? 3 : 6);
716: PetscMalloc2(npoints * 2, &x, npoints, &w);
717: if (npoints == 3) {
718: x[0] = x[1] = x[2] = x[5] = 1.0 / 6.0;
719: x[3] = x[4] = 2.0 / 3.0;
720: w[0] = w[1] = w[2] = 1.0 / 3.0;
721: } else if (npoints == 6) {
722: x[0] = x[1] = x[2] = 0.44594849091597;
723: x[3] = x[4] = 0.10810301816807;
724: x[5] = 0.44594849091597;
725: x[6] = x[7] = x[8] = 0.09157621350977;
726: x[9] = x[10] = 0.81684757298046;
727: x[11] = 0.09157621350977;
728: w[0] = w[1] = w[2] = 0.22338158967801;
729: w[3] = w[4] = w[5] = 0.10995174365532;
730: } else SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Triangle quadrature rules for points 3 and 6 supported; npoints : %D", npoints);
731: PetscQuadratureCreate(PETSC_COMM_SELF, quadrature);
732: PetscQuadratureSetOrder(*quadrature, order);
733: PetscQuadratureSetData(*quadrature, dim, nc, npoints, x, w);
734: /* PetscDTStroudConicalQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature); */
735: } else {
736: PetscDTGaussTensorQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature);
737: }
738: break;
739: case 3:
740: /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */
741: if (nverts == 4) {
742: PetscInt order;
743: const PetscInt npoints = 4; // Choose between 4 and 10
744: PetscMalloc2(npoints * 3, &x, npoints, &w);
745: if (npoints == 4) { /* KEAST1, K1, N=4, O=4 */
746: const PetscReal x_4[12] = { 0.5854101966249685, 0.1381966011250105, 0.1381966011250105,
747: 0.1381966011250105, 0.1381966011250105, 0.1381966011250105,
748: 0.1381966011250105, 0.1381966011250105, 0.5854101966249685,
749: 0.1381966011250105, 0.5854101966249685, 0.1381966011250105
750: };
751: PetscArraycpy(x, x_4, 12);
753: w[0] = w[1] = w[2] = w[3] = 1.0 / 24.0;
754: order = 4;
755: } else if (npoints == 10) { /* KEAST3, K3 N=10, O=10 */
756: const PetscReal x_10[30] = { 0.5684305841968444, 0.1438564719343852, 0.1438564719343852,
757: 0.1438564719343852, 0.1438564719343852, 0.1438564719343852,
758: 0.1438564719343852, 0.1438564719343852, 0.5684305841968444,
759: 0.1438564719343852, 0.5684305841968444, 0.1438564719343852,
760: 0.0000000000000000, 0.5000000000000000, 0.5000000000000000,
761: 0.5000000000000000, 0.0000000000000000, 0.5000000000000000,
762: 0.5000000000000000, 0.5000000000000000, 0.0000000000000000,
763: 0.5000000000000000, 0.0000000000000000, 0.0000000000000000,
764: 0.0000000000000000, 0.5000000000000000, 0.0000000000000000,
765: 0.0000000000000000, 0.0000000000000000, 0.5000000000000000
766: };
767: PetscArraycpy(x, x_10, 30);
769: w[0] = w[1] = w[2] = w[3] = 0.2177650698804054;
770: w[4] = w[5] = w[6] = w[7] = w[8] = w[9] = 0.0214899534130631;
771: order = 10;
772: } else SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Tetrahedral quadrature rules for points 4 and 10 supported; npoints : %D", npoints);
773: PetscQuadratureCreate(PETSC_COMM_SELF, quadrature);
774: PetscQuadratureSetOrder(*quadrature, order);
775: PetscQuadratureSetData(*quadrature, dim, nc, npoints, x, w);
776: /* PetscDTStroudConicalQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature); */
777: } else {
778: PetscDTGaussTensorQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature);
779: }
780: break;
781: default:
782: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %D", dim);
783: }
784: return(0);
785: }
787: /* Compute Jacobians */
788: PetscErrorCode ComputeJacobian_Internal (const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal *quad, PetscReal *phypts,
789: PetscReal *jacobian, PetscReal *ijacobian, PetscReal* dvolume)
790: {
791: PetscInt i;
792: PetscReal volume=1.0;
799: PetscArrayzero(jacobian, dim * dim);
800: if (ijacobian) {
801: PetscArrayzero(ijacobian, dim * dim);
802: }
803: if (phypts) {
804: PetscArrayzero(phypts, /*npts=1 * */ 3);
805: }
807: if (dim == 1) {
808: const PetscReal& r = quad[0];
809: if (nverts == 2) { /* Linear Edge */
810: const PetscReal dNi_dxi[2] = { -1.0, 1.0 };
812: for (i = 0; i < nverts; ++i) {
813: const PetscReal* vertices = coordinates + i * 3;
814: jacobian[0] += dNi_dxi[i] * vertices[0];
815: }
816: } else if (nverts == 3) { /* Quadratic Edge */
817: const PetscReal dNi_dxi[3] = { 4 * r - 3.0, 4 * ( 1.0 - 2.0 * r), 4.0 * r - 1.0};
819: for (i = 0; i < nverts; ++i) {
820: const PetscReal* vertices = coordinates + i * 3;
821: jacobian[0] += dNi_dxi[i] * vertices[0];
822: }
823: } else {
824: 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);
825: }
827: if (ijacobian) {
828: /* invert the jacobian */
829: ijacobian[0] = 1.0 / jacobian[0];
830: }
831: } else if (dim == 2) {
833: if (nverts == 4) { /* Linear Quadrangle */
834: const PetscReal& r = quad[0];
835: const PetscReal& s = quad[1];
837: const PetscReal dNi_dxi[4] = { -1.0 + s, 1.0 - s, s, -s };
838: const PetscReal dNi_deta[4] = { -1.0 + r, -r, r, 1.0 - r };
840: for (i = 0; i < nverts; ++i) {
841: const PetscReal* vertices = coordinates + i * 3;
842: jacobian[0] += dNi_dxi[i] * vertices[0];
843: jacobian[2] += dNi_dxi[i] * vertices[1];
844: jacobian[1] += dNi_deta[i] * vertices[0];
845: jacobian[3] += dNi_deta[i] * vertices[1];
846: }
847: } else if (nverts == 3) { /* Linear triangle */
848: const PetscReal x2 = coordinates[2 * 3 + 0], y2 = coordinates[2 * 3 + 1];
850: /* Jacobian is constant */
851: jacobian[0] = (coordinates[0 * 3 + 0] - x2); jacobian[1] = (coordinates[1 * 3 + 0] - x2);
852: jacobian[2] = (coordinates[0 * 3 + 1] - y2); jacobian[3] = (coordinates[1 * 3 + 1] - y2);
853: } 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);
855: /* invert the jacobian */
856: if (ijacobian) {
857: DMatrix_Invert_2x2_Internal(jacobian, ijacobian, &volume);
858: }
859: } else {
861: if (nverts == 8) { /* Linear Hexahedra */
862: const PetscReal &r = quad[0];
863: const PetscReal &s = quad[1];
864: const PetscReal &t = quad[2];
865: const PetscReal dNi_dxi[8] = {- ( 1.0 - s) * ( 1.0 - t),
866: ( 1.0 - s) * ( 1.0 - t),
867: ( s) * ( 1.0 - t),
868: - ( s) * ( 1.0 - t),
869: - ( 1.0 - s) * ( t),
870: ( 1.0 - s) * ( t),
871: ( s) * ( t),
872: - ( s) * ( t)
873: };
875: const PetscReal dNi_deta[8] = { - ( 1.0 - r) * ( 1.0 - t),
876: - ( r) * ( 1.0 - t),
877: ( r) * ( 1.0 - t),
878: ( 1.0 - r) * ( 1.0 - t),
879: - ( 1.0 - r) * ( t),
880: - ( r) * ( t),
881: ( r) * ( t),
882: ( 1.0 - r) * ( t)
883: };
885: const PetscReal dNi_dzeta[8] = { - ( 1.0 - r) * ( 1.0 - s),
886: - ( r) * ( 1.0 - s),
887: - ( r) * ( s),
888: - ( 1.0 - r) * ( s),
889: ( 1.0 - r) * ( 1.0 - s),
890: ( r) * ( 1.0 - s),
891: ( r) * ( s),
892: ( 1.0 - r) * ( s)
893: };
895: for (i = 0; i < nverts; ++i) {
896: const PetscReal* vertex = coordinates + i * 3;
897: jacobian[0] += dNi_dxi[i] * vertex[0];
898: jacobian[3] += dNi_dxi[i] * vertex[1];
899: jacobian[6] += dNi_dxi[i] * vertex[2];
900: jacobian[1] += dNi_deta[i] * vertex[0];
901: jacobian[4] += dNi_deta[i] * vertex[1];
902: jacobian[7] += dNi_deta[i] * vertex[2];
903: jacobian[2] += dNi_dzeta[i] * vertex[0];
904: jacobian[5] += dNi_dzeta[i] * vertex[1];
905: jacobian[8] += dNi_dzeta[i] * vertex[2];
906: }
907: } else if (nverts == 4) { /* Linear Tetrahedra */
908: const PetscReal x0 = coordinates[/*0 * 3 +*/ 0], y0 = coordinates[/*0 * 3 +*/ 1], z0 = coordinates[/*0 * 3 +*/ 2];
910: /* compute the jacobian */
911: jacobian[0] = coordinates[1 * 3 + 0] - x0; jacobian[1] = coordinates[2 * 3 + 0] - x0; jacobian[2] = coordinates[3 * 3 + 0] - x0;
912: jacobian[3] = coordinates[1 * 3 + 1] - y0; jacobian[4] = coordinates[2 * 3 + 1] - y0; jacobian[5] = coordinates[3 * 3 + 1] - y0;
913: jacobian[6] = coordinates[1 * 3 + 2] - z0; jacobian[7] = coordinates[2 * 3 + 2] - z0; jacobian[8] = coordinates[3 * 3 + 2] - z0;
914: } 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);
916: if (ijacobian) {
917: /* invert the jacobian */
918: DMatrix_Invert_3x3_Internal(jacobian, ijacobian, &volume);
919: }
921: }
922: if (volume < 1e-12) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Element has zero volume: %g. Degenerate element or invalid connectivity\n", volume);
923: if (dvolume) *dvolume = volume;
924: return(0);
925: }
927: PetscErrorCode FEMComputeBasis_JandF(const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal *quadrature, PetscReal *phypts,
928: PetscReal *phibasis, PetscReal *jacobian, PetscReal *ijacobian, PetscReal* volume)
929: {
933: switch (dim) {
934: case 1:
935: Compute_Lagrange_Basis_1D_Internal(nverts, coordinates, 1, quadrature, phypts,
936: NULL, phibasis, NULL, jacobian, ijacobian, volume);
937: break;
938: case 2:
939: Compute_Lagrange_Basis_2D_Internal(nverts, coordinates, 1, quadrature, phypts,
940: NULL, phibasis, NULL, NULL, jacobian, ijacobian, volume);
941: break;
942: case 3:
943: Compute_Lagrange_Basis_3D_Internal(nverts, coordinates, 1, quadrature, phypts,
944: NULL, phibasis, NULL, NULL, NULL, jacobian, ijacobian, volume);
945: break;
946: default:
947: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %D", dim);
948: }
949: return(0);
950: }
952: /*@C
953: DMMoabPToRMapping - Compute the mapping from the physical coordinate system for a given element to the
954: canonical reference element. In addition to finding the inverse mapping evaluation through Newton iteration,
955: the basis function at the parametric point is also evaluated optionally.
957: Input Parameters:
958: + PetscInt dim - the element dimension (1=EDGE, 2=QUAD/TRI, 3=HEX/TET)
959: . PetscInt nverts - the number of vertices in the physical element
960: . PetscReal coordinates - the coordinates of vertices in the physical element
961: - PetscReal[3] xphy - the coordinates of physical point for which natural coordinates (in reference frame) are sought
963: Output Parameters:
964: + PetscReal[3] natparam - the natural coordinates (in reference frame) corresponding to xphy
965: - PetscReal[nverts] phi - the basis functions evaluated at the natural coordinates (natparam)
967: Level: advanced
969: @*/
970: PetscErrorCode DMMoabPToRMapping(const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal* xphy, PetscReal* natparam, PetscReal* phi)
971: {
972: /* Perform inverse evaluation for the mapping with use of Newton Raphson iteration */
973: const PetscReal tol = 1.0e-10;
974: const PetscInt max_iterations = 10;
975: const PetscReal error_tol_sqr = tol*tol;
976: PetscReal phibasis[8], jacobian[9], ijacobian[9], volume;
977: PetscReal phypts[3] = {0.0, 0.0, 0.0};
978: PetscReal delta[3] = {0.0, 0.0, 0.0};
979: PetscErrorCode ierr;
980: PetscInt iters=0;
981: PetscReal error=1.0;
988: PetscArrayzero(jacobian, dim * dim);
989: PetscArrayzero(ijacobian, dim * dim);
990: PetscArrayzero(phibasis, nverts);
992: /* zero initial guess */
993: natparam[0] = natparam[1] = natparam[2] = 0.0;
995: /* Compute delta = evaluate( xi) - x */
996: FEMComputeBasis_JandF(dim, nverts, coordinates, natparam, phypts, phibasis, jacobian, ijacobian, &volume);
998: error = 0.0;
999: switch(dim) {
1000: case 3:
1001: delta[2] = phypts[2] - xphy[2];
1002: error += (delta[2]*delta[2]);
1003: case 2:
1004: delta[1] = phypts[1] - xphy[1];
1005: error += (delta[1]*delta[1]);
1006: case 1:
1007: delta[0] = phypts[0] - xphy[0];
1008: error += (delta[0]*delta[0]);
1009: break;
1010: }
1012: while (error > error_tol_sqr) {
1014: if (++iters > max_iterations)
1015: 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]);
1017: /* Compute natparam -= J.inverse() * delta */
1018: switch(dim) {
1019: case 1:
1020: natparam[0] -= ijacobian[0] * delta[0];
1021: break;
1022: case 2:
1023: natparam[0] -= ijacobian[0] * delta[0] + ijacobian[1] * delta[1];
1024: natparam[1] -= ijacobian[2] * delta[0] + ijacobian[3] * delta[1];
1025: break;
1026: case 3:
1027: natparam[0] -= ijacobian[0] * delta[0] + ijacobian[1] * delta[1] + ijacobian[2] * delta[2];
1028: natparam[1] -= ijacobian[3] * delta[0] + ijacobian[4] * delta[1] + ijacobian[5] * delta[2];
1029: natparam[2] -= ijacobian[6] * delta[0] + ijacobian[7] * delta[1] + ijacobian[8] * delta[2];
1030: break;
1031: }
1033: /* Compute delta = evaluate( xi) - x */
1034: FEMComputeBasis_JandF(dim, nverts, coordinates, natparam, phypts, phibasis, jacobian, ijacobian, &volume);
1036: error = 0.0;
1037: switch(dim) {
1038: case 3:
1039: delta[2] = phypts[2] - xphy[2];
1040: error += (delta[2]*delta[2]);
1041: case 2:
1042: delta[1] = phypts[1] - xphy[1];
1043: error += (delta[1]*delta[1]);
1044: case 1:
1045: delta[0] = phypts[0] - xphy[0];
1046: error += (delta[0]*delta[0]);
1047: break;
1048: }
1049: }
1050: if (phi) {
1051: PetscArraycpy(phi, phibasis, nverts);
1052: }
1053: return(0);
1054: }