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