Actual source code: dmmbfem.cxx
petsc-3.13.6 2020-09-29
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: /*@C
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
111: .vb
112: 1-------2 1----3----2
113: EDGE2 EDGE3
114: .ve
116: 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: @*/
131: PetscErrorCode Compute_Lagrange_Basis_1D_Internal ( const PetscInt nverts, const PetscReal *coords, const PetscInt npts, const PetscReal *quad, PetscReal *phypts,
132: PetscReal *jxw, PetscReal *phi, PetscReal *dphidx,
133: PetscReal *jacobian, PetscReal *ijacobian, PetscReal *volume)
134: {
135: int i, j;
136: PetscErrorCode ierr;
142: if (phypts) {
143: PetscArrayzero(phypts, npts * 3);
144: }
145: if (dphidx) { /* Reset arrays. */
146: PetscArrayzero(dphidx, npts * nverts);
147: }
148: if (nverts == 2) { /* Linear Edge */
150: for (j = 0; j < npts; j++)
151: {
152: const PetscInt offset = j * nverts;
153: const PetscReal r = quad[j];
155: phi[0 + offset] = ( 1.0 - r );
156: phi[1 + offset] = ( r );
158: const PetscReal dNi_dxi[2] = { -1.0, 1.0 };
160: jacobian[0] = ijacobian[0] = volume[0] = 0.0;
161: for (i = 0; i < nverts; ++i) {
162: const PetscReal* vertices = coords + i * 3;
163: jacobian[0] += dNi_dxi[i] * vertices[0];
164: if (phypts) {
165: phypts[3 * j + 0] += phi[i + offset] * vertices[0];
166: }
167: }
169: /* invert the jacobian */
170: *volume = jacobian[0];
171: ijacobian[0] = 1.0 / jacobian[0];
172: jxw[j] *= *volume;
174: /* Divide by element jacobian. */
175: for ( i = 0; i < nverts; i++ ) {
176: if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[0];
177: }
179: }
180: }
181: else if (nverts == 3) { /* Quadratic Edge */
183: for (j = 0; j < npts; j++)
184: {
185: const PetscInt offset = j * nverts;
186: const PetscReal r = quad[j];
188: phi[0 + offset] = 1.0 + r * ( 2.0 * r - 3.0 );
189: phi[1 + offset] = 4.0 * r * ( 1.0 - r );
190: phi[2 + offset] = r * ( 2.0 * r - 1.0 );
192: const PetscReal dNi_dxi[3] = { 4 * r - 3.0, 4 * ( 1.0 - 2.0 * r ), 4.0 * r - 1.0};
194: jacobian[0] = ijacobian[0] = volume[0] = 0.0;
195: for (i = 0; i < nverts; ++i) {
196: const PetscReal* vertices = coords + i * 3;
197: jacobian[0] += dNi_dxi[i] * vertices[0];
198: if (phypts) {
199: phypts[3 * j + 0] += phi[i + offset] * vertices[0];
200: }
201: }
203: /* invert the jacobian */
204: *volume = jacobian[0];
205: ijacobian[0] = 1.0 / jacobian[0];
206: if (jxw) jxw[j] *= *volume;
208: /* Divide by element jacobian. */
209: for ( i = 0; i < nverts; i++ ) {
210: if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[0];
211: }
212: }
213: }
214: else {
215: 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);
216: }
217: #if 0
218: /* verify if the computed basis functions are consistent */
219: for ( j = 0; j < npts; j++ ) {
220: PetscScalar phisum = 0, dphixsum = 0;
221: for ( i = 0; i < nverts; i++ ) {
222: phisum += phi[i + j * nverts];
223: if (dphidx) dphixsum += dphidx[i + j * nverts];
224: }
225: PetscPrintf(PETSC_COMM_WORLD, "Sum of basis at quadrature point %D = %g, %g, %g\n", j, phisum, dphixsum);
226: }
227: #endif
228: return(0);
229: }
232: /*@C
233: Compute_Lagrange_Basis_2D_Internal - Evaluate bases and derivatives at quadrature points for a QUAD4 or TRI3 element.
235: The routine is given the coordinates of the vertices of a quadrangle or triangle.
237: The routine evaluates the basis functions associated with each quadrature point provided,
238: and their derivatives with respect to X and Y.
240: Notes:
242: Example Physical Element (QUAD4)
243: .vb
244: 4------3 s
245: | | |
246: | | |
247: | | |
248: 1------2 0-------r
249: .ve
251: Input Parameter:
252: + PetscInt nverts - the number of element vertices
253: . PetscReal coords[3*nverts] - the physical coordinates of the vertices (in canonical numbering)
254: . PetscInt npts - the number of evaluation points (quadrature points)
255: - PetscReal quad[3*npts] - the evaluation points (quadrature points) in the reference space
257: Output Parameter:
258: + PetscReal phypts[3*npts] - the evaluation points (quadrature points) transformed to the physical space
259: . PetscReal jxw[npts] - the jacobian determinant * quadrature weight necessary for assembling discrete contributions
260: . PetscReal phi[npts] - the bases evaluated at the specified quadrature points
261: . PetscReal dphidx[npts] - the derivative of the bases wrt X-direction evaluated at the specified quadrature points
262: - PetscReal dphidy[npts] - the derivative of the bases wrt Y-direction evaluated at the specified quadrature points
264: Level: advanced
266: @*/
267: PetscErrorCode Compute_Lagrange_Basis_2D_Internal ( const PetscInt nverts, const PetscReal *coords, const PetscInt npts, const PetscReal *quad, PetscReal *phypts,
268: PetscReal *jxw, PetscReal *phi, PetscReal *dphidx, PetscReal *dphidy,
269: PetscReal *jacobian, PetscReal *ijacobian, PetscReal *volume)
270: {
271: PetscInt i, j, k;
272: PetscErrorCode ierr;
278: PetscArrayzero(phi, npts);
279: if (phypts) {
280: PetscArrayzero(phypts, npts * 3);
281: }
282: if (dphidx) { /* Reset arrays. */
283: PetscArrayzero(dphidx, npts * nverts);
284: PetscArrayzero(dphidy, npts * nverts);
285: }
286: if (nverts == 4) { /* Linear Quadrangle */
288: for (j = 0; j < npts; j++)
289: {
290: const PetscInt offset = j * nverts;
291: const PetscReal r = quad[0 + j * 2];
292: const PetscReal s = quad[1 + j * 2];
294: phi[0 + offset] = ( 1.0 - r ) * ( 1.0 - s );
295: phi[1 + offset] = r * ( 1.0 - s );
296: phi[2 + offset] = r * s;
297: phi[3 + offset] = ( 1.0 - r ) * s;
299: const PetscReal dNi_dxi[4] = { -1.0 + s, 1.0 - s, s, -s };
300: const PetscReal dNi_deta[4] = { -1.0 + r, -r, r, 1.0 - r };
302: PetscArrayzero(jacobian, 4);
303: PetscArrayzero(ijacobian, 4);
304: for (i = 0; i < nverts; ++i) {
305: const PetscReal* vertices = coords + i * 3;
306: jacobian[0] += dNi_dxi[i] * vertices[0];
307: jacobian[2] += dNi_dxi[i] * vertices[1];
308: jacobian[1] += dNi_deta[i] * vertices[0];
309: jacobian[3] += dNi_deta[i] * vertices[1];
310: if (phypts) {
311: phypts[3 * j + 0] += phi[i + offset] * vertices[0];
312: phypts[3 * j + 1] += phi[i + offset] * vertices[1];
313: phypts[3 * j + 2] += phi[i + offset] * vertices[2];
314: }
315: }
317: /* invert the jacobian */
318: DMatrix_Invert_2x2_Internal(jacobian, ijacobian, volume);
319: 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);
321: if (jxw) jxw[j] *= *volume;
323: /* Let us compute the bases derivatives by scaling with inverse jacobian. */
324: if (dphidx) {
325: for ( i = 0; i < nverts; i++ ) {
326: for ( k = 0; k < 2; ++k) {
327: if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[k * 2 + 0];
328: if (dphidy) dphidy[i + offset] += dNi_deta[i] * ijacobian[k * 2 + 1];
329: } /* for k=[0..2) */
330: } /* for i=[0..nverts) */
331: } /* if (dphidx) */
332: }
333: }
334: else if (nverts == 3) { /* Linear triangle */
336: PetscArrayzero(jacobian, 4);
337: PetscArrayzero(ijacobian, 4);
339: const PetscReal x2 = coords[2 * 3 + 0], y2 = coords[2 * 3 + 1];
341: /* Jacobian is constant */
342: jacobian[0] = (coords[0 * 3 + 0] - x2); jacobian[1] = (coords[1 * 3 + 0] - x2);
343: jacobian[2] = (coords[0 * 3 + 1] - y2); jacobian[3] = (coords[1 * 3 + 1] - y2);
345: /* invert the jacobian */
346: DMatrix_Invert_2x2_Internal(jacobian, ijacobian, volume);
347: 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);
349: const PetscReal Dx[3] = { ijacobian[0], ijacobian[2], - ijacobian[0] - ijacobian[2] };
350: const PetscReal Dy[3] = { ijacobian[1], ijacobian[3], - ijacobian[1] - ijacobian[3] };
352: for (j = 0; j < npts; j++)
353: {
354: const PetscInt offset = j * nverts;
355: const PetscReal r = quad[0 + j * 2];
356: const PetscReal s = quad[1 + j * 2];
358: if (jxw) jxw[j] *= 0.5;
360: /* const PetscReal Ni[3] = { r, s, 1.0 - r - s }; */
361: const PetscReal phipts_x = coords[2 * 3 + 0] + jacobian[0] * r + jacobian[1] * s;
362: const PetscReal phipts_y = coords[2 * 3 + 1] + jacobian[2] * r + jacobian[3] * s;
363: if (phypts) {
364: phypts[offset + 0] = phipts_x;
365: phypts[offset + 1] = phipts_y;
366: }
368: /* \phi_0 = (b.y - c.y) x + (b.x - c.x) y + c.x b.y - b.x c.y */
369: phi[0 + offset] = ( ijacobian[0] * (phipts_x - x2) + ijacobian[1] * (phipts_y - y2) );
370: /* \phi_1 = (c.y - a.y) x + (a.x - c.x) y + c.x a.y - a.x c.y */
371: phi[1 + offset] = ( ijacobian[2] * (phipts_x - x2) + ijacobian[3] * (phipts_y - y2) );
372: phi[2 + offset] = 1.0 - phi[0 + offset] - phi[1 + offset];
374: if (dphidx) {
375: dphidx[0 + offset] = Dx[0];
376: dphidx[1 + offset] = Dx[1];
377: dphidx[2 + offset] = Dx[2];
378: }
380: if (dphidy) {
381: dphidy[0 + offset] = Dy[0];
382: dphidy[1 + offset] = Dy[1];
383: dphidy[2 + offset] = Dy[2];
384: }
386: }
387: }
388: else {
389: 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);
390: }
391: #if 0
392: /* verify if the computed basis functions are consistent */
393: for ( j = 0; j < npts; j++ ) {
394: PetscScalar phisum = 0, dphixsum = 0, dphiysum = 0;
395: for ( i = 0; i < nverts; i++ ) {
396: phisum += phi[i + j * nverts];
397: if (dphidx) dphixsum += dphidx[i + j * nverts];
398: if (dphidy) dphiysum += dphidy[i + j * nverts];
399: }
400: PetscPrintf(PETSC_COMM_WORLD, "Sum of basis at quadrature point %D = %g, %g, %g\n", j, phisum, dphixsum, dphiysum);
401: }
402: #endif
403: return(0);
404: }
407: /*@C
408: Compute_Lagrange_Basis_3D_Internal - Evaluate bases and derivatives at quadrature points for a HEX8 or TET4 element.
410: The routine is given the coordinates of the vertices of a hexahedra or tetrahedra.
412: The routine evaluates the basis functions associated with each quadrature point provided,
413: and their derivatives with respect to X, Y and Z.
415: Notes:
417: Example Physical Element (HEX8)
418: .vb
419: 8------7
420: /| /| t s
421: 5------6 | | /
422: | | | | |/
423: | 4----|-3 0-------r
424: |/ |/
425: 1------2
426: .ve
428: Input Parameter:
429: + PetscInt nverts - the number of element vertices
430: . PetscReal coords[3*nverts] - the physical coordinates of the vertices (in canonical numbering)
431: . PetscInt npts - the number of evaluation points (quadrature points)
432: - PetscReal quad[3*npts] - the evaluation points (quadrature points) in the reference space
434: Output Parameter:
435: + PetscReal phypts[3*npts] - the evaluation points (quadrature points) transformed to the physical space
436: . PetscReal jxw[npts] - the jacobian determinant * quadrature weight necessary for assembling discrete contributions
437: . PetscReal phi[npts] - the bases evaluated at the specified quadrature points
438: . PetscReal dphidx[npts] - the derivative of the bases wrt X-direction evaluated at the specified quadrature points
439: . PetscReal dphidy[npts] - the derivative of the bases wrt Y-direction evaluated at the specified quadrature points
440: - PetscReal dphidz[npts] - the derivative of the bases wrt Z-direction evaluated at the specified quadrature points
442: Level: advanced
444: @*/
445: PetscErrorCode Compute_Lagrange_Basis_3D_Internal ( const PetscInt nverts, const PetscReal *coords, const PetscInt npts, const PetscReal *quad, PetscReal *phypts,
446: PetscReal *jxw, PetscReal *phi, PetscReal *dphidx, PetscReal *dphidy, PetscReal *dphidz,
447: PetscReal *jacobian, PetscReal *ijacobian, PetscReal *volume)
448: {
449: PetscInt i, j, k;
456: /* Reset arrays. */
457: PetscArrayzero(phi, npts);
458: if (phypts) {
459: PetscArrayzero(phypts, npts * 3);
460: }
461: if (dphidx) {
462: PetscArrayzero(dphidx, npts * nverts);
463: PetscArrayzero(dphidy, npts * nverts);
464: PetscArrayzero(dphidz, npts * nverts);
465: }
467: if (nverts == 8) { /* Linear Hexahedra */
469: for (j = 0; j < npts; j++)
470: {
471: const PetscInt offset = j * nverts;
472: const PetscReal& r = quad[j * 3 + 0];
473: const PetscReal& s = quad[j * 3 + 1];
474: const PetscReal& t = quad[j * 3 + 2];
476: phi[offset + 0] = ( 1.0 - r ) * ( 1.0 - s ) * ( 1.0 - t ); /* 0,0,0 */
477: phi[offset + 1] = ( r ) * ( 1.0 - s ) * ( 1.0 - t ); /* 1,0,0 */
478: phi[offset + 2] = ( r ) * ( s ) * ( 1.0 - t ); /* 1,1,0 */
479: phi[offset + 3] = ( 1.0 - r ) * ( s ) * ( 1.0 - t ); /* 0,1,0 */
480: phi[offset + 4] = ( 1.0 - r ) * ( 1.0 - s ) * ( t ); /* 0,0,1 */
481: phi[offset + 5] = ( r ) * ( 1.0 - s ) * ( t ); /* 1,0,1 */
482: phi[offset + 6] = ( r ) * ( s ) * ( t ); /* 1,1,1 */
483: phi[offset + 7] = ( 1.0 - r ) * ( s ) * ( t ); /* 0,1,1 */
485: const PetscReal dNi_dxi[8] = {- ( 1.0 - s ) * ( 1.0 - t ),
486: ( 1.0 - s ) * ( 1.0 - t ),
487: ( s ) * ( 1.0 - t ),
488: - ( s ) * ( 1.0 - t ),
489: - ( 1.0 - s ) * ( t ),
490: ( 1.0 - s ) * ( t ),
491: ( s ) * ( t ),
492: - ( s ) * ( t )
493: };
495: const PetscReal dNi_deta[8] = { - ( 1.0 - r ) * ( 1.0 - t ),
496: - ( r ) * ( 1.0 - t ),
497: ( r ) * ( 1.0 - t ),
498: ( 1.0 - r ) * ( 1.0 - t ),
499: - ( 1.0 - r ) * ( t ),
500: - ( r ) * ( t ),
501: ( r ) * ( t ),
502: ( 1.0 - r ) * ( t )
503: };
505: const PetscReal dNi_dzeta[8] = { - ( 1.0 - r ) * ( 1.0 - s ),
506: - ( r ) * ( 1.0 - s ),
507: - ( r ) * ( s ),
508: - ( 1.0 - r ) * ( s ),
509: ( 1.0 - r ) * ( 1.0 - s ),
510: ( r ) * ( 1.0 - s ),
511: ( r ) * ( s ),
512: ( 1.0 - r ) * ( s )
513: };
515: PetscArrayzero(jacobian, 9);
516: PetscArrayzero(ijacobian, 9);
517: for (i = 0; i < nverts; ++i) {
518: const PetscReal* vertex = coords + i * 3;
519: jacobian[0] += dNi_dxi[i] * vertex[0];
520: jacobian[3] += dNi_dxi[i] * vertex[1];
521: jacobian[6] += dNi_dxi[i] * vertex[2];
522: jacobian[1] += dNi_deta[i] * vertex[0];
523: jacobian[4] += dNi_deta[i] * vertex[1];
524: jacobian[7] += dNi_deta[i] * vertex[2];
525: jacobian[2] += dNi_dzeta[i] * vertex[0];
526: jacobian[5] += dNi_dzeta[i] * vertex[1];
527: jacobian[8] += dNi_dzeta[i] * vertex[2];
528: if (phypts) {
529: phypts[3 * j + 0] += phi[i + offset] * vertex[0];
530: phypts[3 * j + 1] += phi[i + offset] * vertex[1];
531: phypts[3 * j + 2] += phi[i + offset] * vertex[2];
532: }
533: }
535: /* invert the jacobian */
536: DMatrix_Invert_3x3_Internal(jacobian, ijacobian, volume);
537: 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);
539: if (jxw) jxw[j] *= (*volume);
541: /* Divide by element jacobian. */
542: for ( i = 0; i < nverts; ++i ) {
543: for ( k = 0; k < 3; ++k) {
544: if (dphidx) dphidx[i + offset] += dNi_dxi[i] * ijacobian[0 * 3 + k];
545: if (dphidy) dphidy[i + offset] += dNi_deta[i] * ijacobian[1 * 3 + k];
546: if (dphidz) dphidz[i + offset] += dNi_dzeta[i] * ijacobian[2 * 3 + k];
547: }
548: }
549: }
550: }
551: else if (nverts == 4) { /* Linear Tetrahedra */
553: PetscArrayzero(jacobian, 9);
554: PetscArrayzero(ijacobian, 9);
556: const PetscReal x0 = coords[/*0 * 3 +*/ 0], y0 = coords[/*0 * 3 +*/ 1], z0 = coords[/*0 * 3 +*/ 2];
558: /* compute the jacobian */
559: jacobian[0] = coords[1 * 3 + 0] - x0; jacobian[1] = coords[2 * 3 + 0] - x0; jacobian[2] = coords[3 * 3 + 0] - x0;
560: jacobian[3] = coords[1 * 3 + 1] - y0; jacobian[4] = coords[2 * 3 + 1] - y0; jacobian[5] = coords[3 * 3 + 1] - y0;
561: jacobian[6] = coords[1 * 3 + 2] - z0; jacobian[7] = coords[2 * 3 + 2] - z0; jacobian[8] = coords[3 * 3 + 2] - z0;
563: /* invert the jacobian */
564: DMatrix_Invert_3x3_Internal(jacobian, ijacobian, volume);
565: 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);
567: PetscReal Dx[4]={0,0,0,0}, Dy[4]={0,0,0,0}, Dz[4]={0,0,0,0};
568: if (dphidx) {
569: Dx[0] = ( coords[1 + 2 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3] )
570: - coords[1 + 1 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3] )
571: - coords[1 + 3 * 3] * ( coords[2 + 1 * 3] - coords[2 + 2 * 3] )
572: ) / *volume;
573: Dx[1] = - ( coords[1 + 2 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3] )
574: - coords[1 + 0 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3] )
575: - coords[1 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 2 * 3] )
576: ) / *volume;
577: Dx[2] = ( coords[1 + 1 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3] )
578: - coords[1 + 0 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3] )
579: - coords[1 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 1 * 3] )
580: ) / *volume;
581: Dx[3] = - ( Dx[0] + Dx[1] + Dx[2] );
582: Dy[0] = ( coords[0 + 1 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3] )
583: - coords[0 + 2 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3] )
584: + coords[0 + 3 * 3] * ( coords[2 + 1 * 3] - coords[2 + 2 * 3] )
585: ) / *volume;
586: Dy[1] = - ( coords[0 + 0 * 3] * ( coords[2 + 2 * 3] - coords[2 + 3 * 3] )
587: - coords[0 + 2 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3] )
588: + coords[0 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 2 * 3] )
589: ) / *volume;
590: Dy[2] = ( coords[0 + 0 * 3] * ( coords[2 + 1 * 3] - coords[2 + 3 * 3] )
591: - coords[0 + 1 * 3] * ( coords[2 + 0 * 3] - coords[2 + 3 * 3] )
592: + coords[0 + 3 * 3] * ( coords[2 + 0 * 3] - coords[2 + 1 * 3] )
593: ) / *volume;
594: Dy[3] = - ( Dy[0] + Dy[1] + Dy[2] );
595: Dz[0] = ( coords[0 + 1 * 3] * (coords[1 + 3 * 3] - coords[1 + 2 * 3] )
596: - coords[0 + 2 * 3] * (coords[1 + 3 * 3] - coords[1 + 1 * 3] )
597: + coords[0 + 3 * 3] * (coords[1 + 2 * 3] - coords[1 + 1 * 3] )
598: ) / *volume;
599: Dz[1] = - ( coords[0 + 0 * 3] * (coords[1 + 3 * 3] - coords[1 + 2 * 3] )
600: + coords[0 + 2 * 3] * (coords[1 + 0 * 3] - coords[1 + 3 * 3] )
601: - coords[0 + 3 * 3] * (coords[1 + 0 * 3] - coords[1 + 2 * 3] )
602: ) / *volume;
603: Dz[2] = ( coords[0 + 0 * 3] * (coords[1 + 3 * 3] - coords[1 + 1 * 3] )
604: + coords[0 + 1 * 3] * (coords[1 + 0 * 3] - coords[1 + 3 * 3] )
605: - coords[0 + 3 * 3] * (coords[1 + 0 * 3] - coords[1 + 1 * 3] )
606: ) / *volume;
607: Dz[3] = - ( Dz[0] + Dz[1] + Dz[2] );
608: }
610: for ( j = 0; j < npts; j++ )
611: {
612: const PetscInt offset = j * nverts;
613: const PetscReal& r = quad[j * 3 + 0];
614: const PetscReal& s = quad[j * 3 + 1];
615: const PetscReal& t = quad[j * 3 + 2];
617: if (jxw) jxw[j] *= *volume;
619: phi[offset + 0] = 1.0 - r - s - t;
620: phi[offset + 1] = r;
621: phi[offset + 2] = s;
622: phi[offset + 3] = t;
624: if (phypts) {
625: for (i = 0; i < nverts; ++i) {
626: const PetscScalar* vertices = coords + i * 3;
627: phypts[3 * j + 0] += phi[i + offset] * vertices[0];
628: phypts[3 * j + 1] += phi[i + offset] * vertices[1];
629: phypts[3 * j + 2] += phi[i + offset] * vertices[2];
630: }
631: }
633: /* Now set the derivatives */
634: if (dphidx) {
635: dphidx[0 + offset] = Dx[0];
636: dphidx[1 + offset] = Dx[1];
637: dphidx[2 + offset] = Dx[2];
638: dphidx[3 + offset] = Dx[3];
640: dphidy[0 + offset] = Dy[0];
641: dphidy[1 + offset] = Dy[1];
642: dphidy[2 + offset] = Dy[2];
643: dphidy[3 + offset] = Dy[3];
645: dphidz[0 + offset] = Dz[0];
646: dphidz[1 + offset] = Dz[1];
647: dphidz[2 + offset] = Dz[2];
648: dphidz[3 + offset] = Dz[3];
649: }
651: } /* Tetrahedra -- ends */
652: }
653: else
654: {
655: 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);
656: }
657: #if 0
658: /* verify if the computed basis functions are consistent */
659: for ( j = 0; j < npts; j++ ) {
660: PetscScalar phisum = 0, dphixsum = 0, dphiysum = 0, dphizsum = 0;
661: const int offset = j * nverts;
662: for ( i = 0; i < nverts; i++ ) {
663: phisum += phi[i + offset];
664: if (dphidx) dphixsum += dphidx[i + offset];
665: if (dphidy) dphiysum += dphidy[i + offset];
666: if (dphidz) dphizsum += dphidz[i + offset];
667: 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]);
668: }
669: 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);
670: }
671: #endif
672: return(0);
673: }
677: /*@C
678: DMMoabFEMComputeBasis - Evaluate bases and derivatives at quadrature points for a linear EDGE/QUAD/TRI/HEX/TET element.
680: The routine takes the coordinates of the vertices of an element and computes basis functions associated with
681: each quadrature point provided, and their derivatives with respect to X, Y and Z as appropriate.
683: Input Parameter:
684: + PetscInt nverts, the number of element vertices
685: . PetscReal coords[3*nverts], the physical coordinates of the vertices (in canonical numbering)
686: . PetscInt npts, the number of evaluation points (quadrature points)
687: - PetscReal quad[3*npts], the evaluation points (quadrature points) in the reference space
689: Output Parameter:
690: + PetscReal phypts[3*npts], the evaluation points (quadrature points) transformed to the physical space
691: . PetscReal jxw[npts], the jacobian determinant * quadrature weight necessary for assembling discrete contributions
692: . PetscReal fe_basis[npts], the bases values evaluated at the specified quadrature points
693: - 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
695: Level: advanced
697: @*/
698: PetscErrorCode DMMoabFEMComputeBasis ( const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscQuadrature quadrature,
699: PetscReal *phypts, PetscReal *jacobian_quadrature_weight_product,
700: PetscReal *fe_basis, PetscReal **fe_basis_derivatives)
701: {
702: PetscErrorCode ierr;
703: PetscInt npoints,idim;
704: bool compute_der;
705: const PetscReal *quadpts, *quadwts;
706: PetscReal jacobian[9], ijacobian[9], volume;
712: compute_der = (fe_basis_derivatives != NULL);
714: /* Get the quadrature points and weights for the given quadrature rule */
715: PetscQuadratureGetData(quadrature, &idim, NULL, &npoints, &quadpts, &quadwts);
716: if (idim != dim) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Dimension mismatch: provided (%D) vs quadrature (%D)\n",idim,dim);
717: if (jacobian_quadrature_weight_product) {
718: PetscArraycpy(jacobian_quadrature_weight_product, quadwts, npoints);
719: }
721: switch (dim) {
722: case 1:
723: Compute_Lagrange_Basis_1D_Internal(nverts, coordinates, npoints, quadpts, phypts,
724: jacobian_quadrature_weight_product, fe_basis,
725: (compute_der ? fe_basis_derivatives[0] : NULL),
726: jacobian, ijacobian, &volume);
727: break;
728: case 2:
729: Compute_Lagrange_Basis_2D_Internal(nverts, coordinates, npoints, quadpts, phypts,
730: jacobian_quadrature_weight_product, fe_basis,
731: (compute_der ? fe_basis_derivatives[0] : NULL),
732: (compute_der ? fe_basis_derivatives[1] : NULL),
733: jacobian, ijacobian, &volume);
734: break;
735: case 3:
736: Compute_Lagrange_Basis_3D_Internal(nverts, coordinates, npoints, quadpts, phypts,
737: jacobian_quadrature_weight_product, fe_basis,
738: (compute_der ? fe_basis_derivatives[0] : NULL),
739: (compute_der ? fe_basis_derivatives[1] : NULL),
740: (compute_der ? fe_basis_derivatives[2] : NULL),
741: jacobian, ijacobian, &volume);
742: break;
743: default:
744: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %D", dim);
745: }
746: return(0);
747: }
751: /*@C
752: DMMoabFEMCreateQuadratureDefault - Create default quadrature rules for integration over an element with a given
753: dimension and polynomial order (deciphered from number of element vertices).
755: Input Parameter:
757: + PetscInt dim - the element dimension (1=EDGE, 2=QUAD/TRI, 3=HEX/TET)
758: - PetscInt nverts - the number of vertices in the physical element
760: Output Parameter:
761: . PetscQuadrature quadrature - the quadrature object with default settings to integrate polynomials defined over the element
763: Level: advanced
765: @*/
766: PetscErrorCode DMMoabFEMCreateQuadratureDefault ( const PetscInt dim, const PetscInt nverts, PetscQuadrature *quadrature )
767: {
768: PetscReal *w, *x;
769: PetscInt nc=1;
770: PetscErrorCode ierr;
773: /* Create an appropriate quadrature rule to sample basis */
774: switch (dim)
775: {
776: case 1:
777: /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */
778: PetscDTStroudConicalQuadrature(1, nc, nverts, 0, 1.0, quadrature);
779: break;
780: case 2:
781: /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */
782: if (nverts == 3) {
783: const PetscInt order = 2;
784: const PetscInt npoints = (order == 2 ? 3 : 6);
785: PetscMalloc2(npoints * 2, &x, npoints, &w);
786: if (npoints == 3) {
787: x[0] = x[1] = x[2] = x[5] = 1.0 / 6.0;
788: x[3] = x[4] = 2.0 / 3.0;
789: w[0] = w[1] = w[2] = 1.0 / 3.0;
790: }
791: else if (npoints == 6) {
792: x[0] = x[1] = x[2] = 0.44594849091597;
793: x[3] = x[4] = 0.10810301816807;
794: x[5] = 0.44594849091597;
795: x[6] = x[7] = x[8] = 0.09157621350977;
796: x[9] = x[10] = 0.81684757298046;
797: x[11] = 0.09157621350977;
798: w[0] = w[1] = w[2] = 0.22338158967801;
799: w[3] = w[4] = w[5] = 0.10995174365532;
800: }
801: else {
802: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Triangle quadrature rules for points 3 and 6 supported; npoints : %D", npoints);
803: }
804: PetscQuadratureCreate(PETSC_COMM_SELF, quadrature);
805: PetscQuadratureSetOrder(*quadrature, order);
806: PetscQuadratureSetData(*quadrature, dim, nc, npoints, x, w);
807: /* PetscDTStroudConicalQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature); */
808: }
809: else {
810: PetscDTGaussTensorQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature);
811: }
812: break;
813: case 3:
814: /* Create Gauss quadrature rules with <order = nverts> in the span [-1, 1] */
815: if (nverts == 4) {
816: PetscInt order;
817: const PetscInt npoints = 4; // Choose between 4 and 10
818: PetscMalloc2(npoints * 3, &x, npoints, &w);
819: if (npoints == 4) { /* KEAST1, K1, N=4, O=4 */
820: const PetscReal x_4[12] = { 0.5854101966249685, 0.1381966011250105, 0.1381966011250105,
821: 0.1381966011250105, 0.1381966011250105, 0.1381966011250105,
822: 0.1381966011250105, 0.1381966011250105, 0.5854101966249685,
823: 0.1381966011250105, 0.5854101966249685, 0.1381966011250105
824: };
825: PetscArraycpy(x, x_4, 12);
827: w[0] = w[1] = w[2] = w[3] = 1.0 / 24.0;
828: order = 4;
829: }
830: else if (npoints == 10) { /* KEAST3, K3 N=10, O=10 */
831: const PetscReal x_10[30] = { 0.5684305841968444, 0.1438564719343852, 0.1438564719343852,
832: 0.1438564719343852, 0.1438564719343852, 0.1438564719343852,
833: 0.1438564719343852, 0.1438564719343852, 0.5684305841968444,
834: 0.1438564719343852, 0.5684305841968444, 0.1438564719343852,
835: 0.0000000000000000, 0.5000000000000000, 0.5000000000000000,
836: 0.5000000000000000, 0.0000000000000000, 0.5000000000000000,
837: 0.5000000000000000, 0.5000000000000000, 0.0000000000000000,
838: 0.5000000000000000, 0.0000000000000000, 0.0000000000000000,
839: 0.0000000000000000, 0.5000000000000000, 0.0000000000000000,
840: 0.0000000000000000, 0.0000000000000000, 0.5000000000000000
841: };
842: PetscArraycpy(x, x_10, 30);
844: w[0] = w[1] = w[2] = w[3] = 0.2177650698804054;
845: w[4] = w[5] = w[6] = w[7] = w[8] = w[9] = 0.0214899534130631;
846: order = 10;
847: }
848: else {
849: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Tetrahedral quadrature rules for points 4 and 10 supported; npoints : %D", npoints);
850: }
851: PetscQuadratureCreate(PETSC_COMM_SELF, quadrature);
852: PetscQuadratureSetOrder(*quadrature, order);
853: PetscQuadratureSetData(*quadrature, dim, nc, npoints, x, w);
854: /* PetscDTStroudConicalQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature); */
855: }
856: else {
857: PetscDTGaussTensorQuadrature(dim, nc, nverts, 0.0, 1.0, quadrature);
858: }
859: break;
860: default:
861: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %D", dim);
862: }
863: return(0);
864: }
866: /* Compute Jacobians */
867: PetscErrorCode ComputeJacobian_Internal ( const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal *quad, PetscReal *phypts,
868: PetscReal *jacobian, PetscReal *ijacobian, PetscReal* dvolume)
869: {
870: PetscInt i;
871: PetscReal volume=1.0;
878: PetscArrayzero(jacobian, dim * dim);
879: if (ijacobian) {
880: PetscArrayzero(ijacobian, dim * dim);
881: }
882: if (phypts) {
883: PetscArrayzero(phypts, /*npts=1 * */ 3);
884: }
886: if (dim == 1) {
888: const PetscReal& r = quad[0];
889: if (nverts == 2) { /* Linear Edge */
890: const PetscReal dNi_dxi[2] = { -1.0, 1.0 };
892: for (i = 0; i < nverts; ++i) {
893: const PetscReal* vertices = coordinates + i * 3;
894: jacobian[0] += dNi_dxi[i] * vertices[0];
895: }
896: }
897: else if (nverts == 3) { /* Quadratic Edge */
898: const PetscReal dNi_dxi[3] = { 4 * r - 3.0, 4 * ( 1.0 - 2.0 * r ), 4.0 * r - 1.0};
900: for (i = 0; i < nverts; ++i) {
901: const PetscReal* vertices = coordinates + i * 3;
902: jacobian[0] += dNi_dxi[i] * vertices[0];
903: }
904: }
905: else {
906: 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);
907: }
909: if (ijacobian) {
910: /* invert the jacobian */
911: ijacobian[0] = 1.0 / jacobian[0];
912: }
914: }
915: else if (dim == 2) {
917: if (nverts == 4) { /* Linear Quadrangle */
918: const PetscReal& r = quad[0];
919: const PetscReal& s = quad[1];
921: const PetscReal dNi_dxi[4] = { -1.0 + s, 1.0 - s, s, -s };
922: const PetscReal dNi_deta[4] = { -1.0 + r, -r, r, 1.0 - r };
924: for (i = 0; i < nverts; ++i) {
925: const PetscReal* vertices = coordinates + i * 3;
926: jacobian[0] += dNi_dxi[i] * vertices[0];
927: jacobian[2] += dNi_dxi[i] * vertices[1];
928: jacobian[1] += dNi_deta[i] * vertices[0];
929: jacobian[3] += dNi_deta[i] * vertices[1];
930: }
931: }
932: else if (nverts == 3) { /* Linear triangle */
933: const PetscReal x2 = coordinates[2 * 3 + 0], y2 = coordinates[2 * 3 + 1];
935: /* Jacobian is constant */
936: jacobian[0] = (coordinates[0 * 3 + 0] - x2); jacobian[1] = (coordinates[1 * 3 + 0] - x2);
937: jacobian[2] = (coordinates[0 * 3 + 1] - y2); jacobian[3] = (coordinates[1 * 3 + 1] - y2);
938: }
939: else {
940: 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);
941: }
943: /* invert the jacobian */
944: if (ijacobian) {
945: DMatrix_Invert_2x2_Internal(jacobian, ijacobian, &volume);
946: }
948: }
949: else {
951: if (nverts == 8) { /* Linear Hexahedra */
952: const PetscReal& r = quad[0];
953: const PetscReal& s = quad[1];
954: const PetscReal& t = quad[2];
955: const PetscReal dNi_dxi[8] = {- ( 1.0 - s ) * ( 1.0 - t ),
956: ( 1.0 - s ) * ( 1.0 - t ),
957: ( s ) * ( 1.0 - t ),
958: - ( s ) * ( 1.0 - t ),
959: - ( 1.0 - s ) * ( t ),
960: ( 1.0 - s ) * ( t ),
961: ( s ) * ( t ),
962: - ( s ) * ( t )
963: };
965: const PetscReal dNi_deta[8] = { - ( 1.0 - r ) * ( 1.0 - t ),
966: - ( r ) * ( 1.0 - t ),
967: ( r ) * ( 1.0 - t ),
968: ( 1.0 - r ) * ( 1.0 - t ),
969: - ( 1.0 - r ) * ( t ),
970: - ( r ) * ( t ),
971: ( r ) * ( t ),
972: ( 1.0 - r ) * ( t )
973: };
975: const PetscReal dNi_dzeta[8] = { - ( 1.0 - r ) * ( 1.0 - s ),
976: - ( r ) * ( 1.0 - s ),
977: - ( r ) * ( s ),
978: - ( 1.0 - r ) * ( s ),
979: ( 1.0 - r ) * ( 1.0 - s ),
980: ( r ) * ( 1.0 - s ),
981: ( r ) * ( s ),
982: ( 1.0 - r ) * ( s )
983: };
985: for (i = 0; i < nverts; ++i) {
986: const PetscReal* vertex = coordinates + i * 3;
987: jacobian[0] += dNi_dxi[i] * vertex[0];
988: jacobian[3] += dNi_dxi[i] * vertex[1];
989: jacobian[6] += dNi_dxi[i] * vertex[2];
990: jacobian[1] += dNi_deta[i] * vertex[0];
991: jacobian[4] += dNi_deta[i] * vertex[1];
992: jacobian[7] += dNi_deta[i] * vertex[2];
993: jacobian[2] += dNi_dzeta[i] * vertex[0];
994: jacobian[5] += dNi_dzeta[i] * vertex[1];
995: jacobian[8] += dNi_dzeta[i] * vertex[2];
996: }
998: }
999: else if (nverts == 4) { /* Linear Tetrahedra */
1000: const PetscReal x0 = coordinates[/*0 * 3 +*/ 0], y0 = coordinates[/*0 * 3 +*/ 1], z0 = coordinates[/*0 * 3 +*/ 2];
1002: /* compute the jacobian */
1003: jacobian[0] = coordinates[1 * 3 + 0] - x0; jacobian[1] = coordinates[2 * 3 + 0] - x0; jacobian[2] = coordinates[3 * 3 + 0] - x0;
1004: jacobian[3] = coordinates[1 * 3 + 1] - y0; jacobian[4] = coordinates[2 * 3 + 1] - y0; jacobian[5] = coordinates[3 * 3 + 1] - y0;
1005: jacobian[6] = coordinates[1 * 3 + 2] - z0; jacobian[7] = coordinates[2 * 3 + 2] - z0; jacobian[8] = coordinates[3 * 3 + 2] - z0;
1006: } /* Tetrahedra -- ends */
1007: else
1008: {
1009: 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);
1010: }
1012: if (ijacobian) {
1013: /* invert the jacobian */
1014: DMatrix_Invert_3x3_Internal(jacobian, ijacobian, &volume);
1015: }
1017: }
1018: if ( volume < 1e-12 ) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Element has zero volume: %g. Degenerate element or invalid connectivity\n", volume);
1019: if (dvolume) *dvolume = volume;
1020: return(0);
1021: }
1024: PetscErrorCode FEMComputeBasis_JandF ( const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal *quadrature, PetscReal *phypts,
1025: PetscReal *phibasis, PetscReal *jacobian, PetscReal *ijacobian, PetscReal* volume )
1026: {
1027: PetscErrorCode ierr;
1030: switch (dim) {
1031: case 1:
1032: Compute_Lagrange_Basis_1D_Internal(nverts, coordinates, 1, quadrature, phypts,
1033: NULL, phibasis, NULL, jacobian, ijacobian, volume);
1034: break;
1035: case 2:
1036: Compute_Lagrange_Basis_2D_Internal(nverts, coordinates, 1, quadrature, phypts,
1037: NULL, phibasis, NULL, NULL, jacobian, ijacobian, volume);
1038: break;
1039: case 3:
1040: Compute_Lagrange_Basis_3D_Internal(nverts, coordinates, 1, quadrature, phypts,
1041: NULL, phibasis, NULL, NULL, NULL, jacobian, ijacobian, volume);
1042: break;
1043: default:
1044: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension; should be in [1,3] : %D", dim);
1045: }
1046: return(0);
1047: }
1051: /*@C
1052: DMMoabPToRMapping - Compute the mapping from the physical coordinate system for a given element to the
1053: canonical reference element. In addition to finding the inverse mapping evaluation through Newton iteration,
1054: the basis function at the parametric point is also evaluated optionally.
1056: Input Parameter:
1057: + PetscInt dim - the element dimension (1=EDGE, 2=QUAD/TRI, 3=HEX/TET)
1058: . PetscInt nverts - the number of vertices in the physical element
1059: . PetscReal coordinates - the coordinates of vertices in the physical element
1060: - PetscReal[3] xphy - the coordinates of physical point for which natural coordinates (in reference frame) are sought
1062: Output Parameter:
1063: + PetscReal[3] natparam - the natural coordinates (in reference frame) corresponding to xphy
1064: - PetscReal[nverts] phi - the basis functions evaluated at the natural coordinates (natparam)
1066: Level: advanced
1068: @*/
1069: PetscErrorCode DMMoabPToRMapping( const PetscInt dim, const PetscInt nverts, const PetscReal *coordinates, const PetscReal* xphy, PetscReal* natparam, PetscReal* phi)
1070: {
1071: /* Perform inverse evaluation for the mapping with use of Newton Raphson iteration */
1072: const PetscReal tol = 1.0e-10;
1073: const PetscInt max_iterations = 10;
1074: const PetscReal error_tol_sqr = tol*tol;
1075: PetscReal phibasis[8], jacobian[9], ijacobian[9], volume;
1076: PetscReal phypts[3] = {0.0, 0.0, 0.0};
1077: PetscReal delta[3] = {0.0, 0.0, 0.0};
1078: PetscErrorCode ierr;
1079: PetscInt iters=0;
1080: PetscReal error=1.0;
1087: PetscArrayzero(jacobian, dim * dim);
1088: PetscArrayzero(ijacobian, dim * dim);
1089: PetscArrayzero(phibasis, nverts);
1091: /* zero initial guess */
1092: natparam[0] = natparam[1] = natparam[2] = 0.0;
1094: /* Compute delta = evaluate( xi ) - x */
1095: FEMComputeBasis_JandF ( dim, nverts, coordinates, natparam, phypts, phibasis, jacobian, ijacobian, &volume );
1097: error = 0.0;
1098: switch(dim) {
1099: case 3:
1100: delta[2] = phypts[2] - xphy[2];
1101: error += (delta[2]*delta[2]);
1102: case 2:
1103: delta[1] = phypts[1] - xphy[1];
1104: error += (delta[1]*delta[1]);
1105: case 1:
1106: delta[0] = phypts[0] - xphy[0];
1107: error += (delta[0]*delta[0]);
1108: break;
1109: }
1111: #if 0
1112: PetscInfo4(NULL, "Current point in physical space : (%g, %g, %g) and error : %g\n", xphy[0], xphy[1], xphy[2], error);
1113: #endif
1115: while (error > error_tol_sqr) {
1117: if(++iters > max_iterations)
1118: 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]);
1120: /* Compute natparam -= J.inverse() * delta */
1121: switch(dim) {
1122: case 1:
1123: natparam[0] -= ijacobian[0] * delta[0];
1124: break;
1125: case 2:
1126: natparam[0] -= ijacobian[0] * delta[0] + ijacobian[1] * delta[1];
1127: natparam[1] -= ijacobian[2] * delta[0] + ijacobian[3] * delta[1];
1128: break;
1129: case 3:
1130: natparam[0] -= ijacobian[0] * delta[0] + ijacobian[1] * delta[1] + ijacobian[2] * delta[2];
1131: natparam[1] -= ijacobian[3] * delta[0] + ijacobian[4] * delta[1] + ijacobian[5] * delta[2];
1132: natparam[2] -= ijacobian[6] * delta[0] + ijacobian[7] * delta[1] + ijacobian[8] * delta[2];
1133: break;
1134: }
1136: /* Compute delta = evaluate( xi ) - x */
1137: FEMComputeBasis_JandF ( dim, nverts, coordinates, natparam, phypts, phibasis, jacobian, ijacobian, &volume );
1139: error = 0.0;
1140: switch(dim) {
1141: case 3:
1142: delta[2] = phypts[2] - xphy[2];
1143: error += (delta[2]*delta[2]);
1144: case 2:
1145: delta[1] = phypts[1] - xphy[1];
1146: error += (delta[1]*delta[1]);
1147: case 1:
1148: delta[0] = phypts[0] - xphy[0];
1149: error += (delta[0]*delta[0]);
1150: break;
1151: }
1152: #if 0
1153: PetscInfo5(NULL, "Newton [%D] : Current point in reference space : (%g, %g, %g) and error : %g\n", iters, natparam[0], natparam[1], natparam[2], error);
1154: #endif
1155: }
1156: if (phi) {
1157: PetscArraycpy(phi, phibasis, nverts);
1158: }
1159: return(0);
1160: }