Actual source code: dmmbfem.cxx
petsc-3.9.4 2018-09-11
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: if (!inmat) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_POINTER, "Invalid input matrix specified for 2x2 inversion.");
15: PetscReal det = DMatrix_Determinant_2x2_Internal(inmat);
16: if (outmat) {
17: outmat[0] = inmat[3] / det;
18: outmat[1] = -inmat[1] / det;
19: outmat[2] = -inmat[2] / det;
20: outmat[3] = inmat[0] / det;
21: }
22: if (determinant) *determinant = det;
23: return(0);
24: }
26: static inline PetscReal DMatrix_Determinant_3x3_Internal ( const PetscReal inmat[3 * 3] )
27: {
28: return inmat[0] * (inmat[8] * inmat[4] - inmat[7] * inmat[5])
29: - inmat[3] * (inmat[8] * inmat[1] - inmat[7] * inmat[2])
30: + inmat[6] * (inmat[5] * inmat[1] - inmat[4] * inmat[2]);
31: }
33: static inline PetscErrorCode DMatrix_Invert_3x3_Internal (const PetscReal *inmat, PetscReal *outmat, PetscScalar *determinant)
34: {
35: if (!inmat) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_POINTER, "Invalid input matrix specified for 3x3 inversion.");
36: PetscReal det = DMatrix_Determinant_3x3_Internal(inmat);
37: if (outmat) {
38: outmat[0] = (inmat[8] * inmat[4] - inmat[7] * inmat[5]) / det;
39: outmat[1] = -(inmat[8] * inmat[1] - inmat[7] * inmat[2]) / det;
40: outmat[2] = (inmat[5] * inmat[1] - inmat[4] * inmat[2]) / det;
41: outmat[3] = -(inmat[8] * inmat[3] - inmat[6] * inmat[5]) / det;
42: outmat[4] = (inmat[8] * inmat[0] - inmat[6] * inmat[2]) / det;
43: outmat[5] = -(inmat[5] * inmat[0] - inmat[3] * inmat[2]) / det;
44: outmat[6] = (inmat[7] * inmat[3] - inmat[6] * inmat[4]) / det;
45: outmat[7] = -(inmat[7] * inmat[0] - inmat[6] * inmat[1]) / det;
46: outmat[8] = (inmat[4] * inmat[0] - inmat[3] * inmat[1]) / det;
47: }
48: if (determinant) *determinant = det;
49: return(0);
50: }
52: inline PetscReal DMatrix_Determinant_4x4_Internal ( PetscReal inmat[4 * 4] )
53: {
54: return
55: inmat[0 + 0 * 4] * (
56: inmat[1 + 1 * 4] * ( inmat[2 + 2 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 2 * 4] )
57: - inmat[1 + 2 * 4] * ( inmat[2 + 1 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 1 * 4] )
58: + inmat[1 + 3 * 4] * ( inmat[2 + 1 * 4] * inmat[3 + 2 * 4] - inmat[2 + 2 * 4] * inmat[3 + 1 * 4] ) )
59: - inmat[0 + 1 * 4] * (
60: inmat[1 + 0 * 4] * ( inmat[2 + 2 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 2 * 4] )
61: - inmat[1 + 2 * 4] * ( inmat[2 + 0 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 0 * 4] )
62: + inmat[1 + 3 * 4] * ( inmat[2 + 0 * 4] * inmat[3 + 2 * 4] - inmat[2 + 2 * 4] * inmat[3 + 0 * 4] ) )
63: + inmat[0 + 2 * 4] * (
64: inmat[1 + 0 * 4] * ( inmat[2 + 1 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 1 * 4] )
65: - inmat[1 + 1 * 4] * ( inmat[2 + 0 * 4] * inmat[3 + 3 * 4] - inmat[2 + 3 * 4] * inmat[3 + 0 * 4] )
66: + inmat[1 + 3 * 4] * ( inmat[2 + 0 * 4] * inmat[3 + 1 * 4] - inmat[2 + 1 * 4] * inmat[3 + 0 * 4] ) )
67: - inmat[0 + 3 * 4] * (
68: inmat[1 + 0 * 4] * ( inmat[2 + 1 * 4] * inmat[3 + 2 * 4] - inmat[2 + 2 * 4] * inmat[3 + 1 * 4] )
69: - inmat[1 + 1 * 4] * ( inmat[2 + 0 * 4] * inmat[3 + 2 * 4] - inmat[2 + 2 * 4] * inmat[3 + 0 * 4] )
70: + inmat[1 + 2 * 4] * ( inmat[2 + 0 * 4] * inmat[3 + 1 * 4] - inmat[2 + 1 * 4] * inmat[3 + 0 * 4] ) );
71: }
73: inline PetscErrorCode DMatrix_Invert_4x4_Internal (PetscReal *inmat, PetscReal *outmat, PetscScalar *determinant)
74: {
75: if (!inmat) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_POINTER, "Invalid input matrix specified for 4x4 inversion.");
76: PetscReal det = DMatrix_Determinant_4x4_Internal(inmat);
77: if (outmat) {
78: 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;
79: 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;
80: 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;
81: 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;
82: 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;
83: 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;
84: 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;
85: 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;
86: 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;
87: 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;
88: 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;
89: 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;
90: 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;
91: 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;
92: 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;
93: 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;
94: }
95: if (determinant) *determinant = det;
96: return(0);
97: }
100: /*@
101: Compute_Lagrange_Basis_1D_Internal - Evaluate bases and derivatives at quadrature points for a EDGE2 or EDGE3 element.
103: The routine is given the coordinates of the vertices of a linear or quadratic edge element.
105: The routine evaluates the basis functions associated with each quadrature point provided,
106: and their derivatives with respect to X.
108: Notes:
110: Example Physical Element
112: 1-------2 1----3----2
113: EDGE2 EDGE3
115: Input Parameter:
117: . PetscInt nverts, the number of element vertices
118: . PetscReal coords[3*nverts], the physical coordinates of the vertices (in canonical numbering)
119: . PetscInt npts, the number of evaluation points (quadrature points)
120: . PetscReal quad[3*npts], the evaluation points (quadrature points) in the reference space
122: Output Parameter:
123: . PetscReal phypts[3*npts], the evaluation points (quadrature points) transformed to the physical space
124: . PetscReal jxw[npts], the jacobian determinant * quadrature weight necessary for assembling discrete contributions
125: . PetscReal phi[npts], the bases evaluated at the specified quadrature points
126: . PetscReal dphidx[npts], the derivative of the bases wrt X-direction evaluated at the specified quadrature points
128: Level: advanced
130: .keywords: DMMoab, FEM, 1-D
131: @*/
132: PetscErrorCode Compute_Lagrange_Basis_1D_Internal ( const PetscInt nverts, const PetscReal *coords, const PetscInt npts, const PetscReal *quad, PetscReal *phypts,
133: PetscReal *jxw, PetscReal *phi, PetscReal *dphidx,
134: PetscReal *jacobian, PetscReal *ijacobian, PetscReal *volume)
135: {
136: int i, j;
137: PetscErrorCode ierr;
143: if (phypts) {
144: PetscMemzero(phypts, npts * 3 * sizeof(PetscReal));
145: }
146: if (dphidx) { /* Reset arrays. */
147: PetscMemzero(dphidx, npts * nverts * sizeof(PetscReal));
148: }
149: if (nverts == 2) { /* Linear Edge */
151: for (j = 0; j < npts; j++)
152: {
153: const PetscInt offset = j * nverts;
154: const PetscReal r = quad[j];
156: phi[0 + offset] = ( 1.0 - r );
157: phi[1 + offset] = ( r );
159: const PetscReal dNi_dxi[2] = { -1.0, 1.0 };
161: jacobian[0] = ijacobian[0] = volume[0] = 0.0;
162: for (i = 0; i < nverts; ++i) {
163: const PetscReal* vertices = coords + i * 3;
164: jacobian[0] += dNi_dxi[i] * vertices[0];
165: if (phypts) {
166: phypts[3 * j + 0] += phi[i + offset] * vertices[0];
167: }
168: }
170: /* invert the jacobian */
171: *volume = jacobian[0];
172: ijacobian[0] = 1.0 / jacobian[0];
173: jxw[j] *= *volume;
175: /* Divide by element jacobian. */
176: for ( i = 0; i < nverts; i++ ) {
177: if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[0];
178: }
180: }
181: }
182: else if (nverts == 3) { /* Quadratic Edge */
184: for (j = 0; j < npts; j++)
185: {
186: const PetscInt offset = j * nverts;
187: const PetscReal r = quad[j];
189: phi[0 + offset] = 1.0 + r * ( 2.0 * r - 3.0 );
190: phi[1 + offset] = 4.0 * r * ( 1.0 - r );
191: phi[2 + offset] = r * ( 2.0 * r - 1.0 );
193: const PetscReal dNi_dxi[3] = { 4 * r - 3.0, 4 * ( 1.0 - 2.0 * r ), 4.0 * r - 1.0};
195: jacobian[0] = ijacobian[0] = volume[0] = 0.0;
196: for (i = 0; i < nverts; ++i) {
197: const PetscReal* vertices = coords + i * 3;
198: jacobian[0] += dNi_dxi[i] * vertices[0];
199: if (phypts) {
200: phypts[3 * j + 0] += phi[i + offset] * vertices[0];
201: }
202: }
204: /* invert the jacobian */
205: *volume = jacobian[0];
206: ijacobian[0] = 1.0 / jacobian[0];
207: if (jxw) jxw[j] *= *volume;
209: /* Divide by element jacobian. */
210: for ( i = 0; i < nverts; i++ ) {
211: if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[0];
212: }
213: }
214: }
215: else {
216: 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);
217: }
218: #if 0
219: /* verify if the computed basis functions are consistent */
220: for ( j = 0; j < npts; j++ ) {
221: PetscScalar phisum = 0, dphixsum = 0;
222: for ( i = 0; i < nverts; i++ ) {
223: phisum += phi[i + j * nverts];
224: if (dphidx) dphixsum += dphidx[i + j * nverts];
225: }
226: PetscPrintf(PETSC_COMM_WORLD, "Sum of basis at quadrature point %D = %g, %g, %g\n", j, phisum, dphixsum);
227: }
228: #endif
229: return(0);
230: }
233: /*@
234: Compute_Lagrange_Basis_2D_Internal - Evaluate bases and derivatives at quadrature points for a QUAD4 or TRI3 element.
236: The routine is given the coordinates of the vertices of a quadrangle or triangle.
238: The routine evaluates the basis functions associated with each quadrature point provided,
239: and their derivatives with respect to X and Y.
241: Notes:
243: Example Physical Element (QUAD4)
245: 4------3 s
246: | | |
247: | | |
248: | | |
249: 1------2 0-------r
251: Input Parameter:
253: . PetscInt nverts, the number of element vertices
254: . PetscReal coords[3*nverts], the physical coordinates of the vertices (in canonical numbering)
255: . PetscInt npts, the number of evaluation points (quadrature points)
256: . PetscReal quad[3*npts], the evaluation points (quadrature points) in the reference space
258: Output Parameter:
259: . PetscReal phypts[3*npts], the evaluation points (quadrature points) transformed to the physical space
260: . PetscReal jxw[npts], the jacobian determinant * quadrature weight necessary for assembling discrete contributions
261: . PetscReal phi[npts], the bases evaluated at the specified quadrature points
262: . PetscReal dphidx[npts], the derivative of the bases wrt X-direction evaluated at the specified quadrature points
263: . PetscReal dphidy[npts], the derivative of the bases wrt Y-direction evaluated at the specified quadrature points
265: Level: advanced
267: .keywords: DMMoab, FEM, 2-D
268: @*/
269: PetscErrorCode Compute_Lagrange_Basis_2D_Internal ( const PetscInt nverts, const PetscReal *coords, const PetscInt npts, const PetscReal *quad, PetscReal *phypts,
270: PetscReal *jxw, PetscReal *phi, PetscReal *dphidx, PetscReal *dphidy,
271: PetscReal *jacobian, PetscReal *ijacobian, PetscReal *volume)
272: {
273: PetscInt i, j, k;
274: PetscErrorCode ierr;
280: PetscMemzero(phi, npts * sizeof(PetscReal));
281: if (phypts) {
282: PetscMemzero(phypts, npts * 3 * sizeof(PetscReal));
283: }
284: if (dphidx) { /* Reset arrays. */
285: PetscMemzero(dphidx, npts * nverts * sizeof(PetscReal));
286: PetscMemzero(dphidy, npts * nverts * sizeof(PetscReal));
287: }
288: if (nverts == 4) { /* Linear Quadrangle */
290: for (j = 0; j < npts; j++)
291: {
292: const PetscInt offset = j * nverts;
293: const PetscReal r = quad[0 + j * 2];
294: const PetscReal s = quad[1 + j * 2];
296: phi[0 + offset] = ( 1.0 - r ) * ( 1.0 - s );
297: phi[1 + offset] = r * ( 1.0 - s );
298: phi[2 + offset] = r * s;
299: phi[3 + offset] = ( 1.0 - r ) * s;
301: const PetscReal dNi_dxi[4] = { -1.0 + s, 1.0 - s, s, -s };
302: const PetscReal dNi_deta[4] = { -1.0 + r, -r, r, 1.0 - r };
304: PetscMemzero(jacobian, 4 * sizeof(PetscReal));
305: PetscMemzero(ijacobian, 4 * sizeof(PetscReal));
306: for (i = 0; i < nverts; ++i) {
307: const PetscReal* vertices = coords + i * 3;
308: jacobian[0] += dNi_dxi[i] * vertices[0];
309: jacobian[2] += dNi_dxi[i] * vertices[1];
310: jacobian[1] += dNi_deta[i] * vertices[0];
311: jacobian[3] += dNi_deta[i] * vertices[1];
312: if (phypts) {
313: phypts[3 * j + 0] += phi[i + offset] * vertices[0];
314: phypts[3 * j + 1] += phi[i + offset] * vertices[1];
315: phypts[3 * j + 2] += phi[i + offset] * vertices[2];
316: }
317: }
319: /* invert the jacobian */
320: DMatrix_Invert_2x2_Internal(jacobian, ijacobian, volume);
321: 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);
323: if (jxw) jxw[j] *= *volume;
325: /* Let us compute the bases derivatives by scaling with inverse jacobian. */
326: if (dphidx) {
327: for ( i = 0; i < nverts; i++ ) {
328: for ( k = 0; k < 2; ++k) {
329: if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[k * 2 + 0];
330: if (dphidy) dphidy[i + offset] += dNi_deta[i] * ijacobian[k * 2 + 1];
331: } /* for k=[0..2) */
332: } /* for i=[0..nverts) */
333: } /* if (dphidx) */
334: }
335: }
336: else if (nverts == 3) { /* Linear triangle */
338: PetscMemzero(jacobian, 4 * sizeof(PetscReal));
339: PetscMemzero(ijacobian, 4 * sizeof(PetscReal));
341: const PetscReal x2 = coords[2 * 3 + 0], y2 = coords[2 * 3 + 1];
343: /* Jacobian is constant */
344: jacobian[0] = (coords[0 * 3 + 0] - x2); jacobian[1] = (coords[1 * 3 + 0] - x2);
345: jacobian[2] = (coords[0 * 3 + 1] - y2); jacobian[3] = (coords[1 * 3 + 1] - y2);
347: /* invert the jacobian */
348: DMatrix_Invert_2x2_Internal(jacobian, ijacobian, volume);
349: if ( *volume < 1e-12 ) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Triangular element has zero volume: %g. Degenerate element or invalid connectivity\n", *volume);
351: const PetscReal Dx[3] = { ijacobian[0], ijacobian[2], - ijacobian[0] - ijacobian[2] };
352: const PetscReal Dy[3] = { ijacobian[1], ijacobian[3], - ijacobian[1] - ijacobian[3] };
354: for (j = 0; j < npts; j++)
355: {
356: const PetscInt offset = j * nverts;
357: const PetscReal r = quad[0 + j * 2];
358: const PetscReal s = quad[1 + j * 2];
360: if (jxw) jxw[j] *= 0.5;
362: /* const PetscReal Ni[3] = { r, s, 1.0 - r - s }; */
363: const PetscReal phipts_x = coords[2 * 3 + 0] + jacobian[0] * r + jacobian[1] * s;
364: const PetscReal phipts_y = coords[2 * 3 + 1] + jacobian[2] * r + jacobian[3] * s;
365: if (phypts) {
366: phypts[offset + 0] = phipts_x;
367: phypts[offset + 1] = phipts_y;
368: }
370: /* \phi_0 = (b.y � c.y) x + (b.x � c.x) y + c.x b.y � b.x c.y */
371: phi[0 + offset] = ( ijacobian[0] * (phipts_x - x2) + ijacobian[1] * (phipts_y - y2) );
372: /* \phi_1 = (c.y � a.y) x + (a.x � c.x) y + c.x a.y � a.x c.y */
373: phi[1 + offset] = ( ijacobian[2] * (phipts_x - x2) + ijacobian[3] * (phipts_y - y2) );
374: phi[2 + offset] = 1.0 - phi[0 + offset] - phi[1 + offset];
376: if (dphidx) {
377: dphidx[0 + offset] = Dx[0];
378: dphidx[1 + offset] = Dx[1];
379: dphidx[2 + offset] = Dx[2];
380: }
382: if (dphidy) {
383: dphidy[0 + offset] = Dy[0];
384: dphidy[1 + offset] = Dy[1];
385: dphidy[2 + offset] = Dy[2];
386: }
388: }
389: }
390: else {
391: 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);
392: }
393: #if 0
394: /* verify if the computed basis functions are consistent */
395: for ( j = 0; j < npts; j++ ) {
396: PetscScalar phisum = 0, dphixsum = 0, dphiysum = 0;
397: for ( i = 0; i < nverts; i++ ) {
398: phisum += phi[i + j * nverts];
399: if (dphidx) dphixsum += dphidx[i + j * nverts];
400: if (dphidy) dphiysum += dphidy[i + j * nverts];
401: }
402: PetscPrintf(PETSC_COMM_WORLD, "Sum of basis at quadrature point %D = %g, %g, %g\n", j, phisum, dphixsum, dphiysum);
403: }
404: #endif
405: return(0);
406: }
409: /*@
410: Compute_Lagrange_Basis_3D_Internal - Evaluate bases and derivatives at quadrature points for a HEX8 or TET4 element.
412: The routine is given the coordinates of the vertices of a hexahedra or tetrahedra.
414: The routine evaluates the basis functions associated with each quadrature point provided,
415: and their derivatives with respect to X, Y and Z.
417: Notes:
419: Example Physical Element (HEX8)
421: 8------7
422: /| /| t s
423: 5------6 | | /
424: | | | | |/
425: | 4----|-3 0-------r
426: |/ |/
427: 1------2
429: Input Parameter:
431: . PetscInt nverts, the number of element vertices
432: . PetscReal coords[3*nverts], the physical coordinates of the vertices (in canonical numbering)
433: . PetscInt npts, the number of evaluation points (quadrature points)
434: . PetscReal quad[3*npts], the evaluation points (quadrature points) in the reference space
436: Output Parameter:
437: . PetscReal phypts[3*npts], the evaluation points (quadrature points) transformed to the physical space
438: . PetscReal jxw[npts], the jacobian determinant * quadrature weight necessary for assembling discrete contributions
439: . PetscReal phi[npts], the bases evaluated at the specified quadrature points
440: . PetscReal dphidx[npts], the derivative of the bases wrt X-direction evaluated at the specified quadrature points
441: . PetscReal dphidy[npts], the derivative of the bases wrt Y-direction evaluated at the specified quadrature points
442: . PetscReal dphidz[npts], the derivative of the bases wrt Z-direction evaluated at the specified quadrature points
444: Level: advanced
446: .keywords: DMMoab, FEM, 3-D
447: @*/
448: PetscErrorCode Compute_Lagrange_Basis_3D_Internal ( const PetscInt nverts, const PetscReal *coords, const PetscInt npts, const PetscReal *quad, PetscReal *phypts,
449: PetscReal *jxw, PetscReal *phi, PetscReal *dphidx, PetscReal *dphidy, PetscReal *dphidz,
450: PetscReal *jacobian, PetscReal *ijacobian, PetscReal *volume)
451: {
452: PetscInt i, j, k;
459: /* Reset arrays. */
460: PetscMemzero(phi, npts * sizeof(PetscReal));
461: if (phypts) {
462: PetscMemzero(phypts, npts * 3 * sizeof(PetscReal));
463: }
464: if (dphidx) {
465: PetscMemzero(dphidx, npts * nverts * sizeof(PetscReal));
466: PetscMemzero(dphidy, npts * nverts * sizeof(PetscReal));
467: PetscMemzero(dphidz, npts * nverts * sizeof(PetscReal));
468: }
470: if (nverts == 8) { /* Linear Hexahedra */
472: for (j = 0; j < npts; j++)
473: {
474: const PetscInt offset = j * nverts;
475: const PetscReal& r = quad[j * 3 + 0];
476: const PetscReal& s = quad[j * 3 + 1];
477: const PetscReal& t = quad[j * 3 + 2];
479: phi[offset + 0] = ( 1.0 - r ) * ( 1.0 - s ) * ( 1.0 - t ); /* 0,0,0 */
480: phi[offset + 1] = ( r ) * ( 1.0 - s ) * ( 1.0 - t ); /* 1,0,0 */
481: phi[offset + 2] = ( r ) * ( s ) * ( 1.0 - t ); /* 1,1,0 */
482: phi[offset + 3] = ( 1.0 - r ) * ( s ) * ( 1.0 - t ); /* 0,1,0 */
483: phi[offset + 4] = ( 1.0 - r ) * ( 1.0 - s ) * ( t ); /* 0,0,1 */
484: phi[offset + 5] = ( r ) * ( 1.0 - s ) * ( t ); /* 1,0,1 */
485: phi[offset + 6] = ( r ) * ( s ) * ( t ); /* 1,1,1 */
486: phi[offset + 7] = ( 1.0 - r ) * ( s ) * ( t ); /* 0,1,1 */
488: const PetscReal dNi_dxi[8] = {- ( 1.0 - s ) * ( 1.0 - t ),
489: ( 1.0 - s ) * ( 1.0 - t ),
490: ( s ) * ( 1.0 - t ),
491: - ( s ) * ( 1.0 - t ),
492: - ( 1.0 - s ) * ( t ),
493: ( 1.0 - s ) * ( t ),
494: ( s ) * ( t ),
495: - ( s ) * ( t )
496: };
498: const PetscReal dNi_deta[8] = { - ( 1.0 - r ) * ( 1.0 - t ),
499: - ( r ) * ( 1.0 - t ),
500: ( r ) * ( 1.0 - t ),
501: ( 1.0 - r ) * ( 1.0 - t ),
502: - ( 1.0 - r ) * ( t ),
503: - ( r ) * ( t ),
504: ( r ) * ( t ),
505: ( 1.0 - r ) * ( t )
506: };
508: const PetscReal dNi_dzeta[8] = { - ( 1.0 - r ) * ( 1.0 - s ),
509: - ( r ) * ( 1.0 - s ),
510: - ( r ) * ( s ),
511: - ( 1.0 - r ) * ( s ),
512: ( 1.0 - r ) * ( 1.0 - s ),
513: ( r ) * ( 1.0 - s ),
514: ( r ) * ( s ),
515: ( 1.0 - r ) * ( s )
516: };
518: PetscMemzero(jacobian, 9 * sizeof(PetscReal));
519: PetscMemzero(ijacobian, 9 * sizeof(PetscReal));
520: for (i = 0; i < nverts; ++i) {
521: const PetscReal* vertex = coords + i * 3;
522: jacobian[0] += dNi_dxi[i] * vertex[0];
523: jacobian[3] += dNi_dxi[i] * vertex[1];
524: jacobian[6] += dNi_dxi[i] * vertex[2];
525: jacobian[1] += dNi_deta[i] * vertex[0];
526: jacobian[4] += dNi_deta[i] * vertex[1];
527: jacobian[7] += dNi_deta[i] * vertex[2];
528: jacobian[2] += dNi_dzeta[i] * vertex[0];
529: jacobian[5] += dNi_dzeta[i] * vertex[1];
530: jacobian[8] += dNi_dzeta[i] * vertex[2];
531: if (phypts) {
532: phypts[3 * j + 0] += phi[i + offset] * vertex[0];
533: phypts[3 * j + 1] += phi[i + offset] * vertex[1];
534: phypts[3 * j + 2] += phi[i + offset] * vertex[2];
535: }
536: }
538: /* invert the jacobian */
539: DMatrix_Invert_3x3_Internal(jacobian, ijacobian, volume);
540: 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);
542: if (jxw) jxw[j] *= (*volume);
544: /* Divide by element jacobian. */
545: for ( i = 0; i < nverts; ++i ) {
546: for ( k = 0; k < 3; ++k) {
547: if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[0 * 3 + k];
548: if (dphidy) dphidy[i + offset] += dNi_deta[i] * ijacobian[1 * 3 + k];
549: if (dphidz) dphidz[i + offset] += dNi_dzeta[i] * ijacobian[2 * 3 + k];
550: }
551: }
552: }
553: }
554: else if (nverts == 4) { /* Linear Tetrahedra */
556: PetscMemzero(jacobian, 9 * sizeof(PetscReal));
557: PetscMemzero(ijacobian, 9 * sizeof(PetscReal));
559: const PetscReal x0 = coords[/*0 * 3 +*/ 0], y0 = coords[/*0 * 3 +*/ 1], z0 = coords[/*0 * 3 +*/ 2];
561: /* compute the jacobian */
562: jacobian[0] = coords[1 * 3 + 0] - x0; jacobian[1] = coords[2 * 3 + 0] - x0; jacobian[2] = coords[3 * 3 + 0] - x0;
563: jacobian[3] = coords[1 * 3 + 1] - y0; jacobian[4] = coords[2 * 3 + 1] - y0; jacobian[5] = coords[3 * 3 + 1] - y0;
564: jacobian[6] = coords[1 * 3 + 2] - z0; jacobian[7] = coords[2 * 3 + 2] - z0; jacobian[8] = coords[3 * 3 + 2] - z0;
566: /* invert the jacobian */
567: DMatrix_Invert_3x3_Internal(jacobian, ijacobian, volume);
568: if ( *volume < 1e-8 ) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Tetrahedral element has zero volume: %g. Degenerate element or invalid connectivity\n", *volume);
570: PetscReal Dx[4]={0,0,0,0}, Dy[4]={0,0,0,0}, Dz[4]={0,0,0,0};
571: if (dphidx) {
572: Dx[0] = ( coords[1 + 2 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3] )
573: - coords[1 + 1 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3] )
574: - coords[1 + 3 * 3] * ( coords[2 + 1 * 3] - coords[2 + 2 * 3] )
575: ) / *volume;
576: Dx[1] = - ( coords[1 + 2 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3] )
577: - coords[1 + 0 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3] )
578: - coords[1 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 2 * 3] )
579: ) / *volume;
580: Dx[2] = ( coords[1 + 1 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3] )
581: - coords[1 + 0 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3] )
582: - coords[1 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 1 * 3] )
583: ) / *volume;
584: Dx[3] = - ( Dx[0] + Dx[1] + Dx[2] );
585: Dy[0] = ( coords[0 + 1 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3] )
586: - coords[0 + 2 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3] )
587: + coords[0 + 3 * 3] * ( coords[2 + 1 * 3] - coords[2 + 2 * 3] )
588: ) / *volume;
589: Dy[1] = - ( coords[0 + 0 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3] )
590: - coords[0 + 2 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3] )
591: + coords[0 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 2 * 3] )
592: ) / *volume;
593: Dy[2] = ( coords[0 + 0 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3] )
594: - coords[0 + 1 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3] )
595: + coords[0 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 1 * 3] )
596: ) / *volume;
597: Dy[3] = - ( Dy[0] + Dy[1] + Dy[2] );
598: Dz[0] = ( coords[0 + 1 * 3] * (coords[1 + 3 * 3] - coords[1 + 2 * 3] )
599: - coords[0 + 2 * 3] * (coords[1 + 3 * 3] - coords[1 + 1 * 3] )
600: + coords[0 + 3 * 3] * (coords[1 + 2 * 3] - coords[1 + 1 * 3] )
601: ) / *volume;
602: Dz[1] = - ( coords[0 + 0 * 3] * (coords[1 + 3 * 3] - coords[1 + 2 * 3] )
603: + coords[0 + 2 * 3] * (coords[1 + 0 * 3] - coords[1 + 3 * 3] )
604: - coords[0 + 3 * 3] * (coords[1 + 0 * 3] - coords[1 + 2 * 3] )
605: ) / *volume;
606: Dz[2] = ( coords[0 + 0 * 3] * (coords[1 + 3 * 3] - coords[1 + 1 * 3] )
607: + coords[0 + 1 * 3] * (coords[1 + 0 * 3] - coords[1 + 3 * 3] )
608: - coords[0 + 3 * 3] * (coords[1 + 0 * 3] - coords[1 + 1 * 3] )
609: ) / *volume;
610: Dz[3] = - ( Dz[0] + Dz[1] + Dz[2] );
611: }
613: for ( j = 0; j < npts; j++ )
614: {
615: const PetscInt offset = j * nverts;
616: const PetscReal& r = quad[j * 3 + 0];
617: const PetscReal& s = quad[j * 3 + 1];
618: const PetscReal& t = quad[j * 3 + 2];
620: if (jxw) jxw[j] *= *volume;
622: phi[offset + 0] = 1.0 - r - s - t;
623: phi[offset + 1] = r;
624: phi[offset + 2] = s;
625: phi[offset + 3] = t;
627: if (phypts) {
628: for (i = 0; i < nverts; ++i) {
629: const PetscScalar* vertices = coords + i * 3;
630: phypts[3 * j + 0] += phi[i + offset] * vertices[0];
631: phypts[3 * j + 1] += phi[i + offset] * vertices[1];
632: phypts[3 * j + 2] += phi[i + offset] * vertices[2];
633: }
634: }
636: /* Now set the derivatives */
637: if (dphidx) {
638: dphidx[0 + offset] = Dx[0];
639: dphidx[1 + offset] = Dx[1];
640: dphidx[2 + offset] = Dx[2];
641: dphidx[3 + offset] = Dx[3];
642:
643: dphidy[0 + offset] = Dy[0];
644: dphidy[1 + offset] = Dy[1];
645: dphidy[2 + offset] = Dy[2];
646: dphidy[3 + offset] = Dy[3];
647:
648: dphidz[0 + offset] = Dz[0];
649: dphidz[1 + offset] = Dz[1];
650: dphidz[2 + offset] = Dz[2];
651: dphidz[3 + offset] = Dz[3];
652: }
654: } /* Tetrahedra -- ends */
655: }
656: else
657: {
658: 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);
659: }
660: #if 0
661: /* verify if the computed basis functions are consistent */
662: for ( j = 0; j < npts; j++ ) {
663: PetscScalar phisum = 0, dphixsum = 0, dphiysum = 0, dphizsum = 0;
664: const int offset = j * nverts;
665: for ( i = 0; i < nverts; i++ ) {
666: phisum += phi[i + offset];
667: if (dphidx) dphixsum += dphidx[i + offset];
668: if (dphidy) dphiysum += dphidy[i + offset];
669: if (dphidz) dphizsum += dphidz[i + offset];
670: if (dphidx) PetscPrintf(PETSC_COMM_WORLD, "\t Values [%d]: [JxW] [phi, dphidx, dphidy, dphidz] = %g, %g, %g, %g, %g\n", j, jxw[j], phi[i + offset], dphidx[i + offset], dphidy[i + offset], dphidz[i + offset]);
671: }
672: if (dphidx) PetscPrintf(PETSC_COMM_WORLD, "Sum of basis at quadrature point %D (%g, %g, %g) = %g, %g, %g, %g\n", j, quad[3 * j + 0], quad[3 * j + 1], quad[3 * j + 2], phisum, dphixsum, dphiysum, dphizsum);
673: }
674: #endif
675: return(0);
676: }
680: /*@
681: DMMoabFEMComputeBasis - Evaluate bases and derivatives at quadrature points for a linear EDGE/QUAD/TRI/HEX/TET element.
683: The routine takes the coordinates of the vertices of an element and computes basis functions associated with
684: each quadrature point provided, and their derivatives with respect to X, Y and Z as appropriate.
686: Input Parameter:
688: . PetscInt nverts, the number of element vertices
689: . PetscReal coords[3*nverts], the physical coordinates of the vertices (in canonical numbering)
690: . PetscInt npts, the number of evaluation points (quadrature points)
691: . PetscReal quad[3*npts], the evaluation points (quadrature points) in the reference space
693: Output Parameter:
694: . PetscReal phypts[3*npts], the evaluation points (quadrature points) transformed to the physical space
695: . PetscReal jxw[npts], the jacobian determinant * quadrature weight necessary for assembling discrete contributions
696: . PetscReal fe_basis[npts], the bases values evaluated at the specified quadrature points
697: . 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
699: Level: advanced
701: .keywords: DMMoab, FEM, 3-D
702: @*/
703: PetscErrorCode DMMoabFEMComputeBasis ( const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscQuadrature quadrature,
704: PetscReal *phypts, PetscReal *jacobian_quadrature_weight_product,
705: PetscReal *fe_basis, PetscReal **fe_basis_derivatives)
706: {
707: PetscErrorCode ierr;
708: PetscInt npoints,idim;
709: bool compute_der;
710: const PetscReal *quadpts, *quadwts;
711: PetscReal jacobian[9], ijacobian[9], volume;
717: compute_der = (fe_basis_derivatives != NULL);
719: /* Get the quadrature points and weights for the given quadrature rule */
720: PetscQuadratureGetData(quadrature, &idim, NULL, &npoints, &quadpts, &quadwts);
721: if (idim != dim) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Dimension mismatch: provided (%D) vs quadrature (%D)\n",idim,dim);
722: if (jacobian_quadrature_weight_product) {
723: PetscMemcpy(jacobian_quadrature_weight_product, quadwts, npoints * sizeof(PetscReal));
724: }
726: switch (dim) {
727: case 1:
728: Compute_Lagrange_Basis_1D_Internal(nverts, coordinates, npoints, quadpts, phypts,
729: jacobian_quadrature_weight_product, fe_basis,
730: (compute_der ? fe_basis_derivatives[0] : NULL),
731: jacobian, ijacobian, &volume);
732: break;
733: case 2:
734: Compute_Lagrange_Basis_2D_Internal(nverts, coordinates, npoints, quadpts, phypts,
735: jacobian_quadrature_weight_product, fe_basis,
736: (compute_der ? fe_basis_derivatives[0] : NULL),
737: (compute_der ? fe_basis_derivatives[1] : NULL),
738: jacobian, ijacobian, &volume);
739: break;
740: case 3:
741: Compute_Lagrange_Basis_3D_Internal(nverts, coordinates, npoints, quadpts, phypts,
742: jacobian_quadrature_weight_product, fe_basis,
743: (compute_der ? fe_basis_derivatives[0] : NULL),
744: (compute_der ? fe_basis_derivatives[1] : NULL),
745: (compute_der ? fe_basis_derivatives[2] : NULL),
746: jacobian, ijacobian, &volume);
747: break;
748: default:
749: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %D", dim);
750: }
751: return(0);
752: }
756: /*@
757: DMMoabFEMCreateQuadratureDefault - Create default quadrature rules for integration over an element with a given
758: dimension and polynomial order (deciphered from number of element vertices).
760: Input Parameter:
762: . PetscInt dim, the element dimension (1=EDGE, 2=QUAD/TRI, 3=HEX/TET)
763: . PetscInt nverts, the number of vertices in the physical element
765: Output Parameter:
766: . PetscQuadrature quadrature, the quadrature object with default settings to integrate polynomials defined over the element
768: Level: advanced
770: .keywords: DMMoab, Quadrature, PetscDT
771: @*/
772: PetscErrorCode DMMoabFEMCreateQuadratureDefault ( const PetscInt dim, const PetscInt nverts, PetscQuadrature *quadrature )
773: {
774: PetscReal *w, *x;
775: PetscInt nc=1;
776: PetscErrorCode ierr;
779: /* Create an appropriate quadrature rule to sample basis */
780: switch (dim)
781: {
782: case 1:
783: /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */
784: PetscDTGaussJacobiQuadrature(1, nc, nverts, 0, 1.0, quadrature);
785: break;
786: case 2:
787: /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */
788: if (nverts == 3) {
789: const PetscInt order = 2;
790: const PetscInt npoints = (order == 2 ? 3 : 6);
791: PetscMalloc2(npoints * 2, &x, npoints, &w);
792: if (npoints == 3) {
793: x[0] = x[1] = x[2] = x[5] = 1.0 / 6.0;
794: x[3] = x[4] = 2.0 / 3.0;
795: w[0] = w[1] = w[2] = 1.0 / 3.0;
796: }
797: else if (npoints == 6) {
798: x[0] = x[1] = x[2] = 0.44594849091597;
799: x[3] = x[4] = 0.10810301816807;
800: x[5] = 0.44594849091597;
801: x[6] = x[7] = x[8] = 0.09157621350977;
802: x[9] = x[10] = 0.81684757298046;
803: x[11] = 0.09157621350977;
804: w[0] = w[1] = w[2] = 0.22338158967801;
805: w[3] = w[4] = w[5] = 0.10995174365532;
806: }
807: else {
808: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Triangle quadrature rules for points 3 and 6 supported; npoints : %D", npoints);
809: }
810: PetscQuadratureCreate(PETSC_COMM_SELF, quadrature);
811: PetscQuadratureSetOrder(*quadrature, order);
812: PetscQuadratureSetData(*quadrature, dim, nc, npoints, x, w);
813: /* PetscDTGaussJacobiQuadrature(dim, nverts, 0.0, 1.0, quadrature); */
814: }
815: else {
816: PetscDTGaussTensorQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature);
817: }
818: break;
819: case 3:
820: /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */
821: if (nverts == 4) {
822: PetscInt order;
823: const PetscInt npoints = 4; // Choose between 4 and 10
824: PetscMalloc2(npoints * 3, &x, npoints, &w);
825: if (npoints == 4) { /* KEAST1, K1, N=4, O=4 */
826: const PetscReal x_4[12] = { 0.5854101966249685, 0.1381966011250105, 0.1381966011250105,
827: 0.1381966011250105, 0.1381966011250105, 0.1381966011250105,
828: 0.1381966011250105, 0.1381966011250105, 0.5854101966249685,
829: 0.1381966011250105, 0.5854101966249685, 0.1381966011250105
830: };
831: PetscMemcpy(x, x_4, 12 * sizeof(PetscReal));
833: w[0] = w[1] = w[2] = w[3] = 1.0 / 24.0;
834: order = 4;
835: }
836: else if (npoints == 10) { /* KEAST3, K3 N=10, O=10 */
837: const PetscReal x_10[30] = { 0.5684305841968444, 0.1438564719343852, 0.1438564719343852,
838: 0.1438564719343852, 0.1438564719343852, 0.1438564719343852,
839: 0.1438564719343852, 0.1438564719343852, 0.5684305841968444,
840: 0.1438564719343852, 0.5684305841968444, 0.1438564719343852,
841: 0.0000000000000000, 0.5000000000000000, 0.5000000000000000,
842: 0.5000000000000000, 0.0000000000000000, 0.5000000000000000,
843: 0.5000000000000000, 0.5000000000000000, 0.0000000000000000,
844: 0.5000000000000000, 0.0000000000000000, 0.0000000000000000,
845: 0.0000000000000000, 0.5000000000000000, 0.0000000000000000,
846: 0.0000000000000000, 0.0000000000000000, 0.5000000000000000
847: };
848: PetscMemcpy(x, x_10, 30 * sizeof(PetscReal));
850: w[0] = w[1] = w[2] = w[3] = 0.2177650698804054;
851: w[4] = w[5] = w[6] = w[7] = w[8] = w[9] = 0.0214899534130631;
852: order = 10;
853: }
854: else {
855: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Tetrahedral quadrature rules for points 4 and 10 supported; npoints : %D", npoints);
856: }
857: PetscQuadratureCreate(PETSC_COMM_SELF, quadrature);
858: PetscQuadratureSetOrder(*quadrature, order);
859: PetscQuadratureSetData(*quadrature, dim, nc, npoints, x, w);
860: /* PetscDTGaussJacobiQuadrature(dim, nverts, 0.0, 1.0, quadrature); */
861: }
862: else {
863: PetscDTGaussTensorQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature);
864: }
865: break;
866: default:
867: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %D", dim);
868: }
869: return(0);
870: }
872: /* Compute Jacobians */
873: PetscErrorCode ComputeJacobian_Internal ( const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal *quad, PetscReal *phypts,
874: PetscReal *jacobian, PetscReal *ijacobian, PetscReal* dvolume)
875: {
876: PetscInt i;
877: PetscReal volume=1.0;
884: PetscMemzero(jacobian, dim * dim * sizeof(PetscReal));
885: if (ijacobian) {
886: PetscMemzero(ijacobian, dim * dim * sizeof(PetscReal));
887: }
888: if (phypts) {
889: PetscMemzero(phypts, /*npts=1 * */ 3 * sizeof(PetscReal));
890: }
892: if (dim == 1) {
894: const PetscReal& r = quad[0];
895: if (nverts == 2) { /* Linear Edge */
896: const PetscReal dNi_dxi[2] = { -1.0, 1.0 };
898: for (i = 0; i < nverts; ++i) {
899: const PetscReal* vertices = coordinates + i * 3;
900: jacobian[0] += dNi_dxi[i] * vertices[0];
901: }
902: }
903: else if (nverts == 3) { /* Quadratic Edge */
904: const PetscReal dNi_dxi[3] = { 4 * r - 3.0, 4 * ( 1.0 - 2.0 * r ), 4.0 * r - 1.0};
906: for (i = 0; i < nverts; ++i) {
907: const PetscReal* vertices = coordinates + i * 3;
908: jacobian[0] += dNi_dxi[i] * vertices[0];
909: }
910: }
911: else {
912: 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);
913: }
915: if (ijacobian) {
916: /* invert the jacobian */
917: ijacobian[0] = 1.0 / jacobian[0];
918: }
920: }
921: else if (dim == 2) {
923: if (nverts == 4) { /* Linear Quadrangle */
924: const PetscReal& r = quad[0];
925: const PetscReal& s = quad[1];
927: const PetscReal dNi_dxi[4] = { -1.0 + s, 1.0 - s, s, -s };
928: const PetscReal dNi_deta[4] = { -1.0 + r, -r, r, 1.0 - r };
930: for (i = 0; i < nverts; ++i) {
931: const PetscReal* vertices = coordinates + i * 3;
932: jacobian[0] += dNi_dxi[i] * vertices[0];
933: jacobian[2] += dNi_dxi[i] * vertices[1];
934: jacobian[1] += dNi_deta[i] * vertices[0];
935: jacobian[3] += dNi_deta[i] * vertices[1];
936: }
937: }
938: else if (nverts == 3) { /* Linear triangle */
939: const PetscReal x2 = coordinates[2 * 3 + 0], y2 = coordinates[2 * 3 + 1];
941: /* Jacobian is constant */
942: jacobian[0] = (coordinates[0 * 3 + 0] - x2); jacobian[1] = (coordinates[1 * 3 + 0] - x2);
943: jacobian[2] = (coordinates[0 * 3 + 1] - y2); jacobian[3] = (coordinates[1 * 3 + 1] - y2);
944: }
945: else {
946: 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);
947: }
949: /* invert the jacobian */
950: if (ijacobian) {
951: DMatrix_Invert_2x2_Internal(jacobian, ijacobian, &volume);
952: }
954: }
955: else {
957: if (nverts == 8) { /* Linear Hexahedra */
958: const PetscReal& r = quad[0];
959: const PetscReal& s = quad[1];
960: const PetscReal& t = quad[2];
961: const PetscReal dNi_dxi[8] = {- ( 1.0 - s ) * ( 1.0 - t ),
962: ( 1.0 - s ) * ( 1.0 - t ),
963: ( s ) * ( 1.0 - t ),
964: - ( s ) * ( 1.0 - t ),
965: - ( 1.0 - s ) * ( t ),
966: ( 1.0 - s ) * ( t ),
967: ( s ) * ( t ),
968: - ( s ) * ( t )
969: };
971: const PetscReal dNi_deta[8] = { - ( 1.0 - r ) * ( 1.0 - t ),
972: - ( r ) * ( 1.0 - t ),
973: ( r ) * ( 1.0 - t ),
974: ( 1.0 - r ) * ( 1.0 - t ),
975: - ( 1.0 - r ) * ( t ),
976: - ( r ) * ( t ),
977: ( r ) * ( t ),
978: ( 1.0 - r ) * ( t )
979: };
981: const PetscReal dNi_dzeta[8] = { - ( 1.0 - r ) * ( 1.0 - s ),
982: - ( r ) * ( 1.0 - s ),
983: - ( r ) * ( s ),
984: - ( 1.0 - r ) * ( s ),
985: ( 1.0 - r ) * ( 1.0 - s ),
986: ( r ) * ( 1.0 - s ),
987: ( r ) * ( s ),
988: ( 1.0 - r ) * ( s )
989: };
991: for (i = 0; i < nverts; ++i) {
992: const PetscReal* vertex = coordinates + i * 3;
993: jacobian[0] += dNi_dxi[i] * vertex[0];
994: jacobian[3] += dNi_dxi[i] * vertex[1];
995: jacobian[6] += dNi_dxi[i] * vertex[2];
996: jacobian[1] += dNi_deta[i] * vertex[0];
997: jacobian[4] += dNi_deta[i] * vertex[1];
998: jacobian[7] += dNi_deta[i] * vertex[2];
999: jacobian[2] += dNi_dzeta[i] * vertex[0];
1000: jacobian[5] += dNi_dzeta[i] * vertex[1];
1001: jacobian[8] += dNi_dzeta[i] * vertex[2];
1002: }
1004: }
1005: else if (nverts == 4) { /* Linear Tetrahedra */
1006: const PetscReal x0 = coordinates[/*0 * 3 +*/ 0], y0 = coordinates[/*0 * 3 +*/ 1], z0 = coordinates[/*0 * 3 +*/ 2];
1008: /* compute the jacobian */
1009: jacobian[0] = coordinates[1 * 3 + 0] - x0; jacobian[1] = coordinates[2 * 3 + 0] - x0; jacobian[2] = coordinates[3 * 3 + 0] - x0;
1010: jacobian[3] = coordinates[1 * 3 + 1] - y0; jacobian[4] = coordinates[2 * 3 + 1] - y0; jacobian[5] = coordinates[3 * 3 + 1] - y0;
1011: jacobian[6] = coordinates[1 * 3 + 2] - z0; jacobian[7] = coordinates[2 * 3 + 2] - z0; jacobian[8] = coordinates[3 * 3 + 2] - z0;
1012: } /* Tetrahedra -- ends */
1013: else
1014: {
1015: 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);
1016: }
1018: if (ijacobian) {
1019: /* invert the jacobian */
1020: DMatrix_Invert_3x3_Internal(jacobian, ijacobian, &volume);
1021: }
1023: }
1024: if ( volume < 1e-12 ) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Element has zero volume: %g. Degenerate element or invalid connectivity\n", volume);
1025: if (dvolume) *dvolume = volume;
1026: return(0);
1027: }
1030: PetscErrorCode FEMComputeBasis_JandF ( const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal *quadrature, PetscReal *phypts,
1031: PetscReal *phibasis, PetscReal *jacobian, PetscReal *ijacobian, PetscReal* volume )
1032: {
1033: PetscErrorCode ierr;
1036: switch (dim) {
1037: case 1:
1038: Compute_Lagrange_Basis_1D_Internal(nverts, coordinates, 1, quadrature, phypts,
1039: NULL, phibasis, NULL, jacobian, ijacobian, volume);
1040: break;
1041: case 2:
1042: Compute_Lagrange_Basis_2D_Internal(nverts, coordinates, 1, quadrature, phypts,
1043: NULL, phibasis, NULL, NULL, jacobian, ijacobian, volume);
1044: break;
1045: case 3:
1046: Compute_Lagrange_Basis_3D_Internal(nverts, coordinates, 1, quadrature, phypts,
1047: NULL, phibasis, NULL, NULL, NULL, jacobian, ijacobian, volume);
1048: break;
1049: default:
1050: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %D", dim);
1051: }
1052: return(0);
1053: }
1057: /*@
1058: DMMoabPToRMapping - Compute the mapping from the physical coordinate system for a given element to the
1059: canonical reference element. In addition to finding the inverse mapping evaluation through Newton iteration,
1060: the basis function at the parametric point is also evaluated optionally.
1062: Input Parameter:
1064: . PetscInt dim, the element dimension (1=EDGE, 2=QUAD/TRI, 3=HEX/TET)
1065: . PetscInt nverts, the number of vertices in the physical element
1066: . PetscReal coordinates, the coordinates of vertices in the physical element
1067: . PetscReal[3] xphy, the coordinates of physical point for which natural coordinates (in reference frame) are sought
1069: Output Parameter:
1070: . PetscReal[3] natparam, the natural coordinates (in reference frame) corresponding to xphy
1071: . PetscReal[nverts] phi, the basis functions evaluated at the natural coordinates (natparam)
1073: Level: advanced
1075: .keywords: DMMoab, Mapping, FEM
1076: @*/
1077: PetscErrorCode DMMoabPToRMapping( const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal* xphy, PetscReal* natparam, PetscReal* phi)
1078: {
1079: /* Perform inverse evaluation for the mapping with use of Newton Raphson iteration */
1080: const PetscReal tol = 1.0e-10;
1081: const PetscInt max_iterations = 10;
1082: const PetscReal error_tol_sqr = tol*tol;
1083: PetscReal phibasis[8], jacobian[9], ijacobian[9], volume;
1084: PetscReal phypts[3] = {0.0, 0.0, 0.0};
1085: PetscReal delta[3] = {0.0, 0.0, 0.0};
1086: PetscErrorCode ierr;
1087: PetscInt iters=0;
1088: PetscReal error=1.0;
1095: PetscMemzero(jacobian, dim * dim * sizeof(PetscReal));
1096: PetscMemzero(ijacobian, dim * dim * sizeof(PetscReal));
1097: PetscMemzero(phibasis, nverts * sizeof(PetscReal));
1099: /* zero initial guess */
1100: natparam[0] = natparam[1] = natparam[2] = 0.0;
1102: /* Compute delta = evaluate( xi ) - x */
1103: FEMComputeBasis_JandF ( dim, nverts, coordinates, natparam, phypts, phibasis, jacobian, ijacobian, &volume );
1105: error = 0.0;
1106: switch(dim) {
1107: case 3:
1108: delta[2] = phypts[2] - xphy[2];
1109: error += (delta[2]*delta[2]);
1110: case 2:
1111: delta[1] = phypts[1] - xphy[1];
1112: error += (delta[1]*delta[1]);
1113: case 1:
1114: delta[0] = phypts[0] - xphy[0];
1115: error += (delta[0]*delta[0]);
1116: break;
1117: }
1119: #if 0
1120: PetscInfo4(NULL, "Current point in physical space : (%g, %g, %g) and error : %g\n", xphy[0], xphy[1], xphy[2], error);
1121: #endif
1123: while (error > error_tol_sqr) {
1125: if(++iters > max_iterations)
1126: 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]);
1128: /* Compute natparam -= J.inverse() * delta */
1129: switch(dim) {
1130: case 1:
1131: natparam[0] -= ijacobian[0] * delta[0];
1132: break;
1133: case 2:
1134: natparam[0] -= ijacobian[0] * delta[0] + ijacobian[1] * delta[1];
1135: natparam[1] -= ijacobian[2] * delta[0] + ijacobian[3] * delta[1];
1136: break;
1137: case 3:
1138: natparam[0] -= ijacobian[0] * delta[0] + ijacobian[1] * delta[1] + ijacobian[2] * delta[2];
1139: natparam[1] -= ijacobian[3] * delta[0] + ijacobian[4] * delta[1] + ijacobian[5] * delta[2];
1140: natparam[2] -= ijacobian[6] * delta[0] + ijacobian[7] * delta[1] + ijacobian[8] * delta[2];
1141: break;
1142: }
1144: /* Compute delta = evaluate( xi ) - x */
1145: FEMComputeBasis_JandF ( dim, nverts, coordinates, natparam, phypts, phibasis, jacobian, ijacobian, &volume );
1147: error = 0.0;
1148: switch(dim) {
1149: case 3:
1150: delta[2] = phypts[2] - xphy[2];
1151: error += (delta[2]*delta[2]);
1152: case 2:
1153: delta[1] = phypts[1] - xphy[1];
1154: error += (delta[1]*delta[1]);
1155: case 1:
1156: delta[0] = phypts[0] - xphy[0];
1157: error += (delta[0]*delta[0]);
1158: break;
1159: }
1160: #if 0
1161: PetscInfo5(NULL, "Newton [%D] : Current point in reference space : (%g, %g, %g) and error : %g\n", iters, natparam[0], natparam[1], natparam[2], error);
1162: #endif
1163: }
1164: if (phi) {
1165: PetscMemcpy(phi, phibasis, nverts * sizeof(PetscReal));
1166: }
1167: return(0);
1168: }