v0.13.2
Loading...
Searching...
No Matches
h1.c
Go to the documentation of this file.
1/** \file h1.c
2
3 Based on Hierarchic Finite Element Bases on Unstructured Tetrahedral
4 Meshes, by Mark Ainsworth and Joe Coyle
5 Shape functions for MBTRI and H1 approximation
6
7*/
8
9#include <petscsys.h>
10#include <definitions.h>
11#include <cblas.h>
12
13#include <h1_hdiv_hcurl_l2.h>
14
15static PetscErrorCode ierr;
16
17PetscErrorCode H1_EdgeShapeFunctions_MBTRI(
18 int *sense, int *p, double *N, double *diffN, double *edgeN[3],
19 double *diff_edgeN[3], int GDIM,
20 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
21 double *L, double *diffL,
22 const int dim)) {
23 MoFEMFunctionBeginHot;
24
25 double *edgeN01 = NULL, *edgeN12 = NULL, *edgeN20 = NULL;
26 if (edgeN != NULL) {
27 edgeN01 = edgeN[0];
28 edgeN12 = edgeN[1];
29 edgeN20 = edgeN[2];
30 }
31 double *diff_edgeN01 = NULL, *diff_edgeN12 = NULL, *diff_edgeN20 = NULL;
32 if (diff_edgeN != NULL) {
33 diff_edgeN01 = diff_edgeN[0];
34 diff_edgeN12 = diff_edgeN[1];
35 diff_edgeN20 = diff_edgeN[2];
36 }
37 int P[3];
38 int ee = 0;
39 for (; ee < 3; ee++) {
40 P[ee] = NBEDGE_H1(p[ee]);
41 }
42 int dd = 0;
43 double diff_ksi01[2], diff_ksi12[2], diff_ksi20[2];
44 if (diff_edgeN != NULL) {
45 for (; dd < 2; dd++) {
46 diff_ksi01[dd] = (diffN[1 * 2 + dd] - diffN[0 * 2 + dd]) * sense[0];
47 diff_ksi12[dd] = (diffN[2 * 2 + dd] - diffN[1 * 2 + dd]) * sense[1];
48 diff_ksi20[dd] = (diffN[0 * 2 + dd] - diffN[2 * 2 + dd]) * sense[2];
49 }
50 }
51 int ii = 0;
52 for (; ii < GDIM; ii++) {
53 int node_shift = ii * 3;
54 double ksi01 = (N[node_shift + 1] - N[node_shift + 0]) * sense[0];
55 double ksi12 = (N[node_shift + 2] - N[node_shift + 1]) * sense[1];
56 double ksi20 = (N[node_shift + 0] - N[node_shift + 2]) * sense[2];
57 double L01[p[0] + 1], L12[p[1] + 1], L20[p[2] + 1];
58 double diffL01[2 * (p[0] + 1)], diffL12[2 * (p[1] + 1)],
59 diffL20[2 * (p[2] + 1)];
60 ierr = base_polynomials(p[0], ksi01, diff_ksi01, L01, diffL01, 2);
61 CHKERRQ(ierr);
62 ierr = base_polynomials(p[1], ksi12, diff_ksi12, L12, diffL12, 2);
63 CHKERRQ(ierr);
64 ierr = base_polynomials(p[2], ksi20, diff_ksi20, L20, diffL20, 2);
65 CHKERRQ(ierr);
66 int shift;
67 if (edgeN != NULL) {
68 // edge 01
69 shift = ii * (P[0]);
70 {
71 double *edge_n_ptr = &edgeN01[shift];
72 double *l_ptr = L01;
73 double scalar = N[node_shift + 0] * N[node_shift + 1];
74 int size = P[0];
75 for (size_t jj = 0; jj != size; ++jj, ++l_ptr) {
76 *edge_n_ptr = (*l_ptr) * scalar;
77 ++edge_n_ptr;
78 }
79 }
80
81 // edge 12
82 shift = ii * (P[1]);
83 {
84 double *edge_n_ptr = &edgeN12[shift];
85 double *l_ptr = L12;
86 double scalar = N[node_shift + 1] * N[node_shift + 2];
87 int size = P[1];
88 for (size_t jj = 0; jj != size; ++jj, ++l_ptr) {
89 *edge_n_ptr = (*l_ptr) * scalar;
90 ++edge_n_ptr;
91 }
92 }
93
94 // edge 20
95 shift = ii * (P[2]);
96 {
97 double *edge_n_ptr = &edgeN20[shift];
98 double *l_ptr = L20;
99 double scalar = N[node_shift + 2] * N[node_shift + 0];
100 int size = P[2];
101 for (size_t jj = 0; jj != size; ++jj, ++l_ptr) {
102 *edge_n_ptr = (*l_ptr) * scalar;
103 ++edge_n_ptr;
104 }
105 }
106 }
107 if (diff_edgeN != NULL) {
108 if (P[0] > 0) {
109 // edge01
110 shift = ii * (P[0]);
111 {
112 double *diff_edge_n_ptr = &diff_edgeN01[2 * shift];
113 double *diff_l_x = &diffL01[0 * (p[0] + 1)];
114 double *diff_l_y = &diffL01[1 * (p[0] + 1)];
115 double scalar_x = diffN[2 * 0 + 0] * N[node_shift + 1] +
116 N[node_shift + 0] * diffN[2 * 1 + 0];
117 double scalar_y = diffN[2 * 0 + 1] * N[node_shift + 1] +
118 N[node_shift + 0] * diffN[2 * 1 + 1];
119 double v = N[node_shift + 0] * N[node_shift + 1];
120 double *l_ptr = L01;
121
122 int size = P[0];
123 for (size_t jj = 0; jj != size;
124 ++jj, ++diff_l_x, ++diff_l_y, ++l_ptr) {
125
126 *diff_edge_n_ptr = v * (*diff_l_x) + scalar_x * (*l_ptr);
127 ++diff_edge_n_ptr;
128 *diff_edge_n_ptr = v * (*diff_l_y) + scalar_y * (*l_ptr);
129 ++diff_edge_n_ptr;
130 }
131 }
132
133 }
134 if (P[1] > 0) {
135 // edge12
136 shift = ii * (P[1]);
137 {
138 double *diff_edge_n_ptr = &diff_edgeN12[2 * shift];
139 double *diff_l_x = &diffL12[0 * (p[1] + 1)];
140 double *diff_l_y = &diffL12[1 * (p[1] + 1)];
141 double scalar_x = diffN[2 * 1 + 0] * N[node_shift + 2] +
142 N[node_shift + 1] * diffN[2 * 2 + 0];
143 double scalar_y = diffN[2 * 1 + 1] * N[node_shift + 2] +
144 N[node_shift + 1] * diffN[2 * 2 + 1];
145 double v = N[node_shift + 1] * N[node_shift + 2];
146 double *l_ptr = L12;
147
148 int size = P[1];
149 for (size_t jj = 0; jj != size;
150 ++jj, ++diff_l_x, ++diff_l_y, ++l_ptr) {
151
152 *diff_edge_n_ptr = v * (*diff_l_x) + scalar_x * (*l_ptr);
153 ++diff_edge_n_ptr;
154 *diff_edge_n_ptr = v * (*diff_l_y) + scalar_y * (*l_ptr);
155 ++diff_edge_n_ptr;
156 }
157 }
158
159 }
160 if (P[2] > 0) {
161 // edge20
162 shift = ii * (P[2]);
163
164 {
165 double *diff_edge_n_ptr = &diff_edgeN20[2 * shift];
166 double *diff_l_x = &diffL20[0 * (p[2] + 1)];
167 double *diff_l_y = &diffL20[1 * (p[2] + 1)];
168 double scalar_x = diffN[2 * 2 + 0] * N[node_shift + 0] +
169 N[node_shift + 2] * diffN[2 * 0 + 0];
170 double scalar_y = diffN[2 * 2 + 1] * N[node_shift + 0] +
171 N[node_shift + 2] * diffN[2 * 0 + 1];
172 double v = N[node_shift + 2] * N[node_shift + 0];
173 double *l_ptr = L20;
174
175 int size = P[2];
176 for (size_t jj = 0; jj != size;
177 ++jj, ++diff_l_x, ++diff_l_y, ++l_ptr) {
178
179 *diff_edge_n_ptr = v * (*diff_l_x) + scalar_x * (*l_ptr);
180 ++diff_edge_n_ptr;
181 *diff_edge_n_ptr = v * (*diff_l_y) + scalar_y * (*l_ptr);
182 ++diff_edge_n_ptr;
183 }
184 }
185
186 }
187 }
188 }
189 MoFEMFunctionReturnHot(0);
190}
191PetscErrorCode H1_FaceShapeFunctions_MBTRI(
192 const int *face_nodes, int p, double *N, double *diffN, double *faceN,
193 double *diff_faceN, int GDIM,
194 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
195 double *L, double *diffL,
196 const int dim)) {
197 MoFEMFunctionBeginHot;
198
199 int P = NBFACETRI_H1(p);
200 if (P == 0)
201 MoFEMFunctionReturnHot(0);
202 double diff_ksiL0[2], diff_ksiL1[2];
203 double *diff_ksi_faces[] = {diff_ksiL0, diff_ksiL1};
204 int dd = 0;
205 if (diff_faceN != NULL) {
206 for (; dd < 2; dd++) {
207 diff_ksi_faces[0][dd] =
208 diffN[face_nodes[1] * 2 + dd] - diffN[face_nodes[0] * 2 + dd];
209 diff_ksi_faces[1][dd] =
210 diffN[face_nodes[2] * 2 + dd] - diffN[face_nodes[0] * 2 + dd];
211 }
212 }
213 int ii = 0;
214 for (; ii < GDIM; ii++) {
215 int node_shift = ii * 3;
216 double ksi_faces[2];
217 ksi_faces[0] =
218 N[node_shift + face_nodes[1]] - N[node_shift + face_nodes[0]];
219 ksi_faces[1] =
220 N[node_shift + face_nodes[2]] - N[node_shift + face_nodes[0]];
221 double L0[p + 1], L1[p + 1];
222 double diffL0[2 * (p + 1)], diffL1[2 * (p + 1)];
223 ierr = base_polynomials(p, ksi_faces[0], diff_ksi_faces[0], L0, diffL0, 2);
224 CHKERRQ(ierr);
225 ierr = base_polynomials(p, ksi_faces[1], diff_ksi_faces[1], L1, diffL1, 2);
226 CHKERRQ(ierr);
227 double v = N[node_shift + face_nodes[0]] * N[node_shift + face_nodes[1]] *
228 N[node_shift + face_nodes[2]];
229 double v2[2] = {0, 0};
230 if (diff_faceN != NULL) {
231 double n1n2 =
232 N[node_shift + face_nodes[1]] * N[node_shift + face_nodes[2]];
233 double n0n2 =
234 N[node_shift + face_nodes[0]] * N[node_shift + face_nodes[2]];
235 double n0n1 =
236 N[node_shift + face_nodes[0]] * N[node_shift + face_nodes[1]];
237 dd = 0;
238 for (; dd < 2; dd++) {
239 v2[dd] = diffN[face_nodes[0] * 2 + dd] * n1n2 +
240 diffN[face_nodes[1] * 2 + dd] * n0n2 +
241 diffN[face_nodes[2] * 2 + dd] * n0n1;
242 }
243 }
244 int shift = ii * P;
245 int jj = 0;
246 int oo = 0;
247 for (; oo <= (p - 3); oo++) {
248 int pp0 = 0;
249 for (; pp0 <= oo; pp0++) {
250 int pp1 = oo - pp0;
251 if (pp1 >= 0) {
252 if (faceN != NULL) {
253 faceN[shift + jj] = v * L0[pp0] * L1[pp1];
254 }
255 if (diff_faceN != NULL) {
256 dd = 0;
257 for (; dd < 2; dd++) {
258 double *val = &diff_faceN[2 * shift + 2 * jj + dd];
259 *val = (L0[pp0] * diffL1[dd * (p + 1) + pp1] +
260 diffL0[dd * (p + 1) + pp0] * L1[pp1]) *
261 v;
262 *val += L0[pp0] * L1[pp1] * v2[dd];
263 }
264 }
265 jj++;
266 }
267 }
268 }
269 if (jj != P)
270 SETERRQ2(PETSC_COMM_SELF, 1, "wrong order %d != %d", jj, P);
271 }
272 MoFEMFunctionReturnHot(0);
273}
274PetscErrorCode H1_EdgeShapeFunctions_MBTET(
275 int *sense, int *p, double *N, double *diffN, double *edgeN[],
276 double *diff_edgeN[], int GDIM,
277 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
278 double *L, double *diffL,
279 const int dim)) {
280 MoFEMFunctionBeginHot;
281
282 int P[6];
283 int ee = 0;
284 for (; ee < 6; ee++)
285 P[ee] = NBEDGE_H1(p[ee]);
286 int edges_nodes[2 * 6] = {0, 1, 1, 2, 2, 0, 0, 3, 1, 3, 2, 3};
287 double diff_ksi01[3], diff_ksi12[3], diff_ksi20[3];
288 double diff_ksi03[3], diff_ksi13[3], diff_ksi23[3];
289 double *edges_diff_ksi[6] = {diff_ksi01, diff_ksi12, diff_ksi20,
290 diff_ksi03, diff_ksi13, diff_ksi23};
291 for (ee = 0; ee < 6; ee++) {
292 int dd = 0;
293 for (; dd < 3; dd++) {
294 edges_diff_ksi[ee][dd] = (diffN[edges_nodes[2 * ee + 1] * 3 + dd] -
295 diffN[edges_nodes[2 * ee + 0] * 3 + dd]) *
296 sense[ee];
297 }
298 }
299 int ii = 0;
300 for (; ii < GDIM; ii++) {
301 int node_shift = ii * 4;
302 double edges_ksi[6];
303 for (ee = 0; ee < 6; ee++) {
304 edges_ksi[ee] = (N[node_shift + edges_nodes[2 * ee + 1]] -
305 N[node_shift + edges_nodes[2 * ee + 0]]) *
306 sense[ee];
307 }
308 for (ee = 0; ee < 6; ee++) {
309 if (P[ee] == 0)
310 continue;
311 double L[p[ee] + 1], diffL[3 * (p[ee] + 1)];
312 ierr = base_polynomials(p[ee], edges_ksi[ee], edges_diff_ksi[ee], L,
313 diffL, 3);
314 CHKERRQ(ierr);
315 double v = N[node_shift + edges_nodes[2 * ee + 0]] *
316 N[node_shift + edges_nodes[2 * ee + 1]];
317 if (edgeN != NULL)
318 if (edgeN[ee] != NULL) {
319 int shift = ii * P[ee];
320 {
321 double *edge_n_ptr = &edgeN[ee][shift];
322 double *l_ptr = L;
323 double scalar = N[node_shift + edges_nodes[2 * ee + 0]] *
324 N[node_shift + edges_nodes[2 * ee + 1]];
325 int size = P[ee];
326 for (size_t jj = 0; jj != size; ++jj, ++l_ptr) {
327 *edge_n_ptr = (*l_ptr) * scalar;
328 ++edge_n_ptr;
329 }
330 }
331 }
332 if (diff_edgeN != NULL)
333 if (diff_edgeN[ee] != NULL) {
334 int shift = ii * P[ee];
335 {
336 double *diff_edge_n_ptr = &diff_edgeN[ee][3 * shift];
337 double *diff_l_x = &diffL[0 * (p[ee] + 1)];
338 double *diff_l_y = &diffL[1 * (p[ee] + 1)];
339 double *diff_l_z = &diffL[2 * (p[ee] + 1)];
340 double *l_ptr = L;
341 double scalar_x = diffN[3 * edges_nodes[2 * ee + 0] + 0] *
342 N[node_shift + edges_nodes[2 * ee + 1]] +
343 N[node_shift + edges_nodes[2 * ee + 0]] *
344 diffN[3 * edges_nodes[2 * ee + 1] + 0];
345 double scalar_y = diffN[3 * edges_nodes[2 * ee + 0] + 1] *
346 N[node_shift + edges_nodes[2 * ee + 1]] +
347 N[node_shift + edges_nodes[2 * ee + 0]] *
348 diffN[3 * edges_nodes[2 * ee + 1] + 1];
349 double scalar_z = diffN[3 * edges_nodes[2 * ee + 0] + 2] *
350 N[node_shift + edges_nodes[2 * ee + 1]] +
351 N[node_shift + edges_nodes[2 * ee + 0]] *
352 diffN[3 * edges_nodes[2 * ee + 1] + 2];
353
354 int size = P[ee];
355 for (size_t jj = 0; jj != size;
356 ++jj, ++diff_l_x, ++diff_l_y, ++diff_l_z, ++l_ptr) {
357
358 *diff_edge_n_ptr = v * (*diff_l_x) + scalar_x * (*l_ptr);
359 ++diff_edge_n_ptr;
360 *diff_edge_n_ptr = v * (*diff_l_y) + scalar_y * (*l_ptr);
361 ++diff_edge_n_ptr;
362 *diff_edge_n_ptr = v * (*diff_l_z) + scalar_z * (*l_ptr);
363 ++diff_edge_n_ptr;
364
365 }
366 }
367
368 }
369 }
370 }
371 MoFEMFunctionReturnHot(0);
372}
373PetscErrorCode H1_FaceShapeFunctions_MBTET(
374 int *faces_nodes, int *p, double *N, double *diffN, double *faceN[],
375 double *diff_faceN[], int GDIM,
376 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
377 double *L, double *diffL,
378 const int dim)) {
379 MoFEMFunctionBeginHot;
380
381 int P[4];
382 int ff = 0;
383 for (; ff < 4; ff++)
384 P[ff] = NBFACETRI_H1(p[ff]);
385 double diff_ksiL0F0[3], diff_ksiL1F0[3];
386 double diff_ksiL0F1[3], diff_ksiL1F1[3];
387 double diff_ksiL0F2[3], diff_ksiL1F2[3];
388 double diff_ksiL0F3[3], diff_ksiL1F3[3];
389 double *diff_ksi_faces[] = {diff_ksiL0F0, diff_ksiL1F0, diff_ksiL0F1,
390 diff_ksiL1F1, diff_ksiL0F2, diff_ksiL1F2,
391 diff_ksiL0F3, diff_ksiL1F3};
392 int dd = 0;
393 for (; dd < 3; dd++) {
394 for (ff = 0; ff < 4; ff++) {
395 diff_ksi_faces[ff * 2 + 0][dd] =
396 (diffN[faces_nodes[3 * ff + 1] * 3 + dd] -
397 diffN[faces_nodes[3 * ff + 0] * 3 + dd]);
398 diff_ksi_faces[ff * 2 + 1][dd] =
399 (diffN[faces_nodes[3 * ff + 2] * 3 + dd] -
400 diffN[faces_nodes[3 * ff + 0] * 3 + dd]);
401 }
402 }
403 int ii = 0;
404 for (; ii < GDIM; ii++) {
405 int node_shift = ii * 4;
406 double ksi_faces[8];
407 for (ff = 0; ff < 4; ff++) {
408 ksi_faces[2 * ff + 0] = N[node_shift + faces_nodes[3 * ff + 1]] -
409 N[node_shift + faces_nodes[3 * ff + 0]];
410 ksi_faces[2 * ff + 1] = N[node_shift + faces_nodes[3 * ff + 2]] -
411 N[node_shift + faces_nodes[3 * ff + 0]];
412 }
413 int shift;
414 for (ff = 0; ff < 4; ff++) {
415 if (P[ff] == 0)
416 continue;
417 double L0[p[ff] + 1], L1[p[ff] + 1];
418 double diffL0[3 * (p[ff] + 1)], diffL1[3 * (p[ff] + 1)];
419 ierr = base_polynomials(p[ff], ksi_faces[ff * 2 + 0],
420 diff_ksi_faces[ff * 2 + 0], L0, diffL0, 3);
421 CHKERRQ(ierr);
422 ierr = base_polynomials(p[ff], ksi_faces[ff * 2 + 1],
423 diff_ksi_faces[ff * 2 + 1], L1, diffL1, 3);
424 CHKERRQ(ierr);
425 double v = N[node_shift + faces_nodes[3 * ff + 0]] *
426 N[node_shift + faces_nodes[3 * ff + 1]] *
427 N[node_shift + faces_nodes[3 * ff + 2]];
428 double v2[3] = {0, 0, 0};
429 dd = 0;
430 double n1n2 = N[node_shift + faces_nodes[3 * ff + 1]] *
431 N[node_shift + faces_nodes[3 * ff + 2]];
432 double n0n2 = N[node_shift + faces_nodes[3 * ff + 0]] *
433 N[node_shift + faces_nodes[3 * ff + 2]];
434 double n0n1 = N[node_shift + faces_nodes[3 * ff + 0]] *
435 N[node_shift + faces_nodes[3 * ff + 1]];
436 for (; dd < 3; dd++) {
437 v2[dd] = diffN[3 * faces_nodes[3 * ff + 0] + dd] * n1n2 +
438 diffN[3 * faces_nodes[3 * ff + 1] + dd] * n0n2 +
439 diffN[3 * faces_nodes[3 * ff + 2] + dd] * n0n1;
440 }
441 shift = ii * P[ff];
442 int jj = 0;
443 int oo = 0;
444 for (; oo <= (p[ff] - 3); oo++) {
445 int pp0 = 0;
446 for (; pp0 <= oo; pp0++) {
447 int pp1 = oo - pp0;
448 if (pp1 >= 0) {
449 if (faceN != NULL)
450 if (faceN[ff] != NULL) {
451 faceN[ff][shift + jj] = v * L0[pp0] * L1[pp1];
452 }
453 if (diff_faceN != NULL)
454 if (diff_faceN[ff] != NULL) {
455 dd = 0;
456 double L0L1 = L0[pp0] * L1[pp1];
457 for (; dd < 3; dd++) {
458 diff_faceN[ff][3 * shift + 3 * jj + dd] =
459 (L0[pp0] * diffL1[dd * (p[ff] + 1) + pp1] +
460 diffL0[dd * (p[ff] + 1) + pp0] * L1[pp1]) *
461 v;
462 diff_faceN[ff][3 * shift + 3 * jj + dd] += L0L1 * v2[dd];
463 }
464 }
465 jj++;
466 }
467 }
468 }
469 if (jj != P[ff])
470 SETERRQ2(PETSC_COMM_SELF, 1, "wrong order %d != %d", jj, P[ff]);
471 }
472 }
473 MoFEMFunctionReturnHot(0);
474}
475PetscErrorCode H1_VolumeShapeFunctions_MBTET(
476 int p, double *N, double *diffN, double *volumeN, double *diff_volumeN,
477 int GDIM,
478 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
479 double *L, double *diffL,
480 const int dim)) {
481 MoFEMFunctionBeginHot;
482
483 int P = NBVOLUMETET_H1(p);
484 if (P == 0)
485 MoFEMFunctionReturnHot(0);
486 double diff_ksiL0[3], diff_ksiL1[3], diff_ksiL2[3];
487 int dd = 0;
488 for (; dd < 3; dd++) {
489 diff_ksiL0[dd] = (diffN[1 * 3 + dd] - diffN[0 * 3 + dd]);
490 diff_ksiL1[dd] = (diffN[2 * 3 + dd] - diffN[0 * 3 + dd]);
491 diff_ksiL2[dd] = (diffN[3 * 3 + dd] - diffN[0 * 3 + dd]);
492 }
493 int ii = 0;
494 for (; ii < GDIM; ii++) {
495 int node_shift = ii * 4;
496 double ksiL0 = N[node_shift + 1] - N[node_shift + 0];
497 double ksiL1 = N[node_shift + 2] - N[node_shift + 0];
498 double ksiL2 = N[node_shift + 3] - N[node_shift + 0];
499 double L0[p + 1], L1[p + 1], L2[p + 1];
500 double diffL0[3 * (p + 1)], diffL1[3 * (p + 1)], diffL2[3 * (p + 1)];
501 ierr = base_polynomials(p, ksiL0, diff_ksiL0, L0, diffL0, 3);
502 CHKERRQ(ierr);
503 ierr = base_polynomials(p, ksiL1, diff_ksiL1, L1, diffL1, 3);
504 CHKERRQ(ierr);
505 ierr = base_polynomials(p, ksiL2, diff_ksiL2, L2, diffL2, 3);
506 CHKERRQ(ierr);
507 double v = N[node_shift + 0] * N[node_shift + 1] * N[node_shift + 2] *
508 N[node_shift + 3];
509 double v2[3] = {0, 0, 0};
510 dd = 0;
511 for (; dd < 3; dd++) {
512 v2[dd] = diffN[3 * 0 + dd] * N[node_shift + 1] * N[node_shift + 2] *
513 N[node_shift + 3] +
514 N[node_shift + 0] * diffN[3 * 1 + dd] * N[node_shift + 2] *
515 N[node_shift + 3] +
516 N[node_shift + 0] * N[node_shift + 1] * diffN[3 * 2 + dd] *
517 N[node_shift + 3] +
518 N[node_shift + 0] * N[node_shift + 1] * N[node_shift + 2] *
519 diffN[3 * 3 + dd];
520 }
521 int shift = ii * P;
522 int jj = 0;
523 int oo = 0;
524 for (; oo <= (p - 4); oo++) {
525 int pp0 = 0;
526 for (; pp0 <= oo; pp0++) {
527 int pp1 = 0;
528 for (; (pp0 + pp1) <= oo; pp1++) {
529 int pp2 = oo - pp0 - pp1;
530 if (pp2 >= 0) {
531 if (volumeN != NULL) {
532 volumeN[shift + jj] = L0[pp0] * L1[pp1] * L2[pp2] * v;
533 }
534 if (diff_volumeN != NULL) {
535 dd = 0;
536 for (; dd < 3; dd++) {
537 diff_volumeN[3 * shift + 3 * jj + dd] =
538 (diffL0[dd * (p + 1) + pp0] * L1[pp1] * L2[pp2] +
539 L0[pp0] * diffL1[dd * (p + 1) + pp1] * L2[pp2] +
540 L0[pp0] * L1[pp1] * diffL2[dd * (p + 1) + pp2]) *
541 v;
542 diff_volumeN[3 * shift + 3 * jj + dd] +=
543 L0[pp0] * L1[pp1] * L2[pp2] * v2[dd];
544 }
545 }
546 jj++;
547 }
548 }
549 }
550 }
551 if (jj != P)
552 SETERRQ1(PETSC_COMM_SELF, 1, "wrong order %d", jj);
553 }
554 MoFEMFunctionReturnHot(0);
555}
556PetscErrorCode H1_EdgeShapeDiffMBTETinvJ(int *base_p, int *p,
557 double *edge_diffN[], double *invJac,
558 double *edge_diffNinvJac[], int GDIM) {
559 MoFEMFunctionBeginHot;
560 int ee = 0, ii, gg;
561 for (; ee < 6; ee++) {
562 for (ii = 0; ii < NBEDGE_H1(p[ee]); ii++) {
563 for (gg = 0; gg < GDIM; gg++) {
564 int shift1 = NBEDGE_H1(base_p[ee]) * gg;
565 int shift2 = NBEDGE_H1(p[ee]) * gg;
566 cblas_dgemv(CblasRowMajor, CblasTrans, 3, 3, 1., invJac, 3,
567 &(edge_diffN[ee])[3 * shift1 + 3 * ii], 1, 0.,
568 &(edge_diffNinvJac[ee])[3 * shift2 + 3 * ii], 1);
569 }
570 }
571 }
572 MoFEMFunctionReturnHot(0);
573}
574PetscErrorCode H1_FaceShapeDiffMBTETinvJ(int *base_p, int *p,
575 double *face_diffN[], double *invJac,
576 double *face_diffNinvJac[], int GDIM) {
577 MoFEMFunctionBeginHot;
578 int ff = 0, ii, gg;
579 for (; ff < 4; ff++) {
580 for (ii = 0; ii < NBFACETRI_H1(p[ff]); ii++) {
581 for (gg = 0; gg < GDIM; gg++) {
582 int shift1 = NBFACETRI_H1(base_p[ff]) * gg;
583 int shift2 = NBFACETRI_H1(p[ff]) * gg;
584 cblas_dgemv(CblasRowMajor, CblasTrans, 3, 3, 1., invJac, 3,
585 &(face_diffN[ff])[3 * shift1 + 3 * ii], 1, 0.,
586 &(face_diffNinvJac[ff])[3 * shift2 + 3 * ii], 1);
587 }
588 }
589 }
590 MoFEMFunctionReturnHot(0);
591}
592PetscErrorCode H1_VolumeShapeDiffMBTETinvJ(int base_p, int p,
593 double *volume_diffN, double *invJac,
594 double *volume_diffNinvJac,
595 int GDIM) {
596 MoFEMFunctionBeginHot;
597 int ii, gg;
598 for (ii = 0; ii < NBVOLUMETET_H1(p); ii++) {
599 for (gg = 0; gg < GDIM; gg++) {
600 int shift1 = NBVOLUMETET_H1(base_p) * gg;
601 int shift2 = NBVOLUMETET_H1(p) * gg;
602 cblas_dgemv(CblasRowMajor, CblasTrans, 3, 3, 1., invJac, 3,
603 &(volume_diffN)[3 * shift1 + 3 * ii], 1, 0.,
604 &(volume_diffNinvJac)[3 * shift2 + 3 * ii], 1);
605 }
606 }
607 MoFEMFunctionReturnHot(0);
608}
609PetscErrorCode H1_EdgeGradientOfDeformation_hierarchical(int p, double *diffN,
610 double *dofs,
611 double *F) {
612 MoFEMFunctionBeginHot;
613 int col, row = 0;
614 for (; row < 3; row++)
615 for (col = 0; col < 3; col++)
616 F[3 * row + col] =
617 cblas_ddot(NBEDGE_H1(p), &diffN[col], 3, &dofs[row], 3);
618 MoFEMFunctionReturnHot(0);
619}
620PetscErrorCode H1_FaceGradientOfDeformation_hierarchical(int p, double *diffN,
621 double *dofs,
622 double *F) {
623 MoFEMFunctionBeginHot;
624 int col, row = 0;
625 for (; row < 3; row++)
626 for (col = 0; col < 3; col++)
627 F[3 * row + col] =
628 cblas_ddot(NBFACETRI_H1(p), &diffN[col], 3, &dofs[row], 3);
629 MoFEMFunctionReturnHot(0);
630}
631PetscErrorCode H1_VolumeGradientOfDeformation_hierarchical(int p, double *diffN,
632 double *dofs,
633 double *F) {
634 MoFEMFunctionBeginHot;
635 int col, row = 0;
636 for (; row < 3; row++)
637 for (col = 0; col < 3; col++)
638 F[3 * row + col] =
639 cblas_ddot(NBVOLUMETET_H1(p), &diffN[col], 3, &dofs[row], 3);
640 MoFEMFunctionReturnHot(0);
641}
642PetscErrorCode H1_QuadShapeFunctions_MBPRISM(
643 int *faces_nodes, int *p, double *N, double *diffN, double *faceN[],
644 double *diff_faceN[], int GDIM,
645 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
646 double *L, double *diffL,
647 const int dim)) {
648 MoFEMFunctionBeginHot;
649 // TODO: return separately components of the tensor product between two edges
650
651 int P[3];
652 int ff = 0;
653 for (; ff < 3; ff++)
654 P[ff] = NBFACEQUAD_H1(p[ff]);
655
656 double ksi_faces[6];
657 double diff_ksiL0F0[3], diff_ksiL3F0[3];
658 double diff_ksiL0F1[3], diff_ksiL3F1[3];
659 double diff_ksiL0F2[3], diff_ksiL3F2[3];
660 double *diff_ksi_faces[] = {diff_ksiL0F0, diff_ksiL3F0, diff_ksiL0F1,
661 diff_ksiL3F1, diff_ksiL0F2, diff_ksiL3F2};
662 int ii = 0;
663 for (; ii < GDIM; ii++) {
664 int node_shift = ii * 6;
665 int node_diff_shift = 3 * node_shift;
666 int ff = 0;
667 for (; ff < 3; ff++) {
668 if (P[ff] == 0)
669 continue;
670 int n0 = faces_nodes[4 * ff + 0];
671 int n1 = faces_nodes[4 * ff + 1];
672 int n2 = faces_nodes[4 * ff + 2];
673 int n3 = faces_nodes[4 * ff + 3];
674 int e0 = 2 * ff + 0;
675 int e1 = 2 * ff + 1;
676 double ksi0 = N[node_shift + n0] + N[node_shift + n3];
677 double ksi1 = N[node_shift + n1] + N[node_shift + n2];
678 double eta0 = N[node_shift + n0] + N[node_shift + n1];
679 double eta1 = N[node_shift + n2] + N[node_shift + n3];
680 ksi_faces[e0] = ksi1 - ksi0;
681 ksi_faces[e1] = eta1 - eta0;
682
683 int dd = 0;
684 for (; dd < 3; dd++) {
685 double diff_ksi0 = diffN[node_diff_shift + 3 * n0 + dd] +
686 diffN[node_diff_shift + 3 * n3 + dd];
687 double diff_ksi1 = diffN[node_diff_shift + 3 * n1 + dd] +
688 diffN[node_diff_shift + 3 * n2 + dd];
689 double diff_eta0 = diffN[node_diff_shift + 3 * n0 + dd] +
690 diffN[node_diff_shift + 3 * n1 + dd];
691 double diff_eta1 = diffN[node_diff_shift + 3 * n2 + dd] +
692 diffN[node_diff_shift + 3 * n3 + dd];
693 diff_ksi_faces[e0][dd] = diff_ksi1 - diff_ksi0;
694 diff_ksi_faces[e1][dd] = diff_eta1 - diff_eta0;
695 }
696 double L0[p[ff] + 1], L1[p[ff] + 1];
697 double diffL0[3 * (p[ff] + 1)], diffL1[3 * (p[ff] + 1)];
698 ierr = base_polynomials(p[ff], ksi_faces[e0], diff_ksi_faces[e0], L0,
699 diffL0, 3);
700 CHKERRQ(ierr);
701 ierr = base_polynomials(p[ff], ksi_faces[e1], diff_ksi_faces[e1], L1,
702 diffL1, 3);
703 CHKERRQ(ierr);
704
705 double v = N[node_shift + n0] * N[node_shift + n2] +
706 N[node_shift + n1] * N[node_shift + n3];
707
708 double diff_v[3];
709 dd = 0;
710 for (; dd < 3; ++dd)
711 diff_v[dd] =
712
713 diffN[node_diff_shift + 3 * n0 + dd] * N[node_shift + n2] +
714 N[node_shift + n0] * diffN[node_diff_shift + 3 * n2 + dd] +
715
716 diffN[node_diff_shift + 3 * n1 + dd] * N[node_shift + n3] +
717 N[node_shift + n1] * diffN[node_diff_shift + 3 * n3 + dd];
718
719 int shift;
720 shift = ii * P[ff];
721
722 if (faceN != NULL) {
723 if (faceN[ff] != NULL) {
724 int jj = 0;
725 int oo = 0;
726 for (; oo <= p[ff] - 2; ++oo) {
727 int pp0 = 0;
728 for (; pp0 < oo; ++pp0) {
729 int pp1 = oo;
730 faceN[ff][shift + jj] = L0[pp0] * L1[pp1] * v;
731 ++jj;
732 faceN[ff][shift + jj] = L0[pp1] * L1[pp0] * v;
733 ++jj;
734 }
735 faceN[ff][shift + jj] = L0[oo] * L1[oo] * v;
736 ++jj;
737 }
738 if (jj != P[ff])
739 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
740 "Inconsistent implementation (bug in the code) %d != %d",
741 jj, P);
742 }
743 }
744
745 if (diff_faceN != NULL) {
746 if (diff_faceN[ff] != NULL) {
747 int jj = 0;
748 int oo = 0;
749 for (; oo <= p[ff] - 2; ++oo) {
750 int pp0 = 0;
751 for (; pp0 < oo; ++pp0) {
752 int pp1 = oo;
753 int dd;
754 for (dd = 0; dd < 3; dd++) {
755 diff_faceN[ff][3 * shift + 3 * jj + dd] =
756 (L0[pp0] * diffL1[dd * (p[ff] + 1) + pp1] +
757 diffL0[dd * (p[ff] + 1) + pp0] * L1[pp1]) *
758 v +
759 L0[pp0] * L1[pp1] * diff_v[dd];
760 }
761 ++jj;
762 for (dd = 0; dd < 3; dd++) {
763 diff_faceN[ff][3 * shift + 3 * jj + dd] =
764 (L0[pp1] * diffL1[dd * (p[ff] + 1) + pp0] +
765 diffL0[dd * (p[ff] + 1) + pp1] * L1[pp0]) *
766 v +
767 L0[pp1] * L1[pp0] * diff_v[dd];
768 }
769 ++jj;
770 }
771 for (dd = 0; dd < 3; dd++) {
772 diff_faceN[ff][3 * shift + 3 * jj + dd] =
773 (L0[oo] * diffL1[dd * (p[ff] + 1) + oo] +
774 diffL0[dd * (p[ff] + 1) + oo] * L1[oo]) *
775 v +
776 L0[oo] * L1[oo] * diff_v[dd];
777 }
778 ++jj;
779 }
780 if (jj != P[ff])
781 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
782 "Inconsistent implementation (bug in the code) %d != %d",
783 jj, P);
784 }
785 }
786 }
787 }
788 MoFEMFunctionReturnHot(0);
789}
790PetscErrorCode H1_VolumeShapeFunctions_MBPRISM(
791 int p, double *N, double *diffN, double *volumeN, double *diff_volumeN,
792 int GDIM,
793 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
794 double *L, double *diffL,
795 const int dim)) {
796 MoFEMFunctionBeginHot;
797
798 int P = NBVOLUMEPRISM_H1(p);
799 if (P == 0)
800 MoFEMFunctionReturnHot(0);
801 double diff_ksiL0[3], diff_ksiL1[3], diff_ksiL2[3];
802 int ii = 0;
803 for (; ii < GDIM; ii++) {
804 int node_shift = ii * 6;
805 int node_diff_shift = ii * 18;
806
807 double n0 = N[node_shift + 0];
808 double n1 = N[node_shift + 1];
809 double n2 = N[node_shift + 2];
810 double n3 = N[node_shift + 3];
811 double n4 = N[node_shift + 4];
812 double n5 = N[node_shift + 5];
813
814 double ksiL0 = n1 + n4 - n0 - n3;
815 double ksiL1 = n2 + n5 - n0 - n3;
816 double ksiL2 = (n3 + n4 + n5) - (n0 + n1 + n2);
817
818 int dd = 0;
819 for (; dd < 3; dd++) {
820 double diff_n0 = diffN[node_diff_shift + 3 * 0 + dd];
821 double diff_n1 = diffN[node_diff_shift + 3 * 1 + dd];
822 double diff_n2 = diffN[node_diff_shift + 3 * 2 + dd];
823 double diff_n3 = diffN[node_diff_shift + 3 * 3 + dd];
824 double diff_n4 = diffN[node_diff_shift + 3 * 4 + dd];
825 double diff_n5 = diffN[node_diff_shift + 3 * 5 + dd];
826 diff_ksiL0[dd] = (diff_n1 + diff_n4) - (diff_n0 + diff_n3);
827 diff_ksiL1[dd] = (diff_n2 + diff_n5) - (diff_n0 + diff_n3);
828 diff_ksiL2[dd] =
829 (diff_n3 + diff_n4 + diff_n5) - (diff_n0 + diff_n1 + diff_n2);
830 }
831
832 double L0[p + 1], L1[p + 1], L2[p + 1];
833 double diffL0[3 * (p + 1)], diffL1[3 * (p + 1)], diffL2[3 * (p + 1)];
834 ierr = base_polynomials(p, ksiL0, diff_ksiL0, L0, diffL0, 3);
835 CHKERRQ(ierr);
836 ierr = base_polynomials(p, ksiL1, diff_ksiL1, L1, diffL1, 3);
837 CHKERRQ(ierr);
838 ierr = base_polynomials(p, ksiL2, diff_ksiL2, L2, diffL2, 3);
839 CHKERRQ(ierr);
840
841 double v_tri = (n0 + n3) * (n1 + n4) * (n2 + n5);
842 double v_edge = (n0 + n1 + n2) * (n3 + n4 + n5);
843 double v = v_tri * v_edge;
844
845 double diff_v_tri[3];
846 double diff_v_edge[3];
847 double diff_v[3];
848 dd = 0;
849 for (; dd < 3; dd++) {
850 double diff_n0 = diffN[node_diff_shift + 3 * 0 + dd];
851 double diff_n1 = diffN[node_diff_shift + 3 * 1 + dd];
852 double diff_n2 = diffN[node_diff_shift + 3 * 2 + dd];
853 double diff_n3 = diffN[node_diff_shift + 3 * 3 + dd];
854 double diff_n4 = diffN[node_diff_shift + 3 * 4 + dd];
855 double diff_n5 = diffN[node_diff_shift + 3 * 5 + dd];
856 diff_v_tri[dd] = (diff_n0 + diff_n3) * (n1 + n4) * (n2 + n5) +
857 (diff_n1 + diff_n4) * (n0 + n3) * (n2 + n5) +
858 (diff_n2 + diff_n5) * (n0 + n3) * (n1 + n4);
859 diff_v_edge[dd] = (diff_n0 + diff_n1 + diff_n2) * (n3 + n4 + n5) +
860 (diff_n3 + diff_n4 + diff_n5) * (n0 + n1 + n2);
861 diff_v[dd] = diff_v_tri[dd] * v_edge + v_tri * diff_v_edge[dd];
862 }
863
864 int shift = ii * P;
865
866 if (volumeN != NULL) {
867 int jj = 0;
868 int oo, pp0, pp1, zz;
869 for (oo = 0; oo <= p - 3; ++oo) {
870 for (pp0 = 0; pp0 < oo; ++pp0) {
871 for (pp1 = 0; pp1 < oo; ++pp1) {
872 const int perm[3][3] = {
873 {oo, pp0, pp1}, {pp0, oo, pp1}, {pp0, pp1, oo}};
874 for (zz = 0; zz != 3; ++zz) {
875 volumeN[shift + jj] =
876 L0[perm[zz][0]] * L1[perm[zz][1]] * L2[perm[zz][2]] * v;
877 ++jj;
878 }
879 }
880 const int perm[3][3] = {{pp0, oo, oo}, {oo, pp0, oo}, {oo, oo, pp0}};
881 for (zz = 0; zz != 3; ++zz) {
882 volumeN[shift + jj] =
883 L0[perm[zz][0]] * L1[perm[zz][1]] * L2[perm[zz][2]] * v;
884 ++jj;
885 }
886 }
887 volumeN[shift + jj] = L0[oo] * L1[oo] * L2[oo] * v;
888 ++jj;
889 }
890 if (jj != P)
891 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
892 "Inconsistent implementation (bug in the code) %d != %d", jj,
893 P);
894 }
895
896 if (diff_volumeN != NULL) {
897 int jj = 0;
898 int oo, pp0, pp1, zz, dd;
899 for (oo = 0; oo <= p - 3; ++oo) {
900 for (pp0 = 0; pp0 < oo; ++pp0) {
901 for (pp1 = 0; pp1 < oo; ++pp1) {
902 const int perm[3][3] = {
903 {oo, pp0, pp1}, {pp0, oo, pp1}, {pp0, pp1, oo}};
904 for (zz = 0; zz != 3; ++zz) {
905 const int i = perm[zz][0];
906 const int j = perm[zz][1];
907 const int k = perm[zz][2];
908 for (dd = 0; dd != 3; ++dd) {
909 diff_volumeN[3 * shift + 3 * jj + dd] =
910 (diffL0[dd * (p + 1) + i] * L1[j] * L2[k] +
911 L0[i] * diffL1[dd * (p + 1) + j] * L2[k] +
912 L0[i] * L1[j] * diffL2[dd * (p + 1) + k]) *
913 v +
914 L0[i] * L1[j] * L2[k] * diff_v[dd];
915 }
916 ++jj;
917 }
918 }
919 const int perm[3][3] = {{pp0, oo, oo}, {oo, pp0, oo}, {oo, oo, pp0}};
920 for (zz = 0; zz != 3; ++zz) {
921 const int i = perm[zz][0];
922 const int j = perm[zz][1];
923 const int k = perm[zz][2];
924 for (dd = 0; dd != 3; ++dd) {
925 diff_volumeN[3 * shift + 3 * jj + dd] =
926 (diffL0[dd * (p + 1) + i] * L1[j] * L2[k] +
927 L0[i] * diffL1[dd * (p + 1) + j] * L2[k] +
928 L0[i] * L1[j] * diffL2[dd * (p + 1) + k]) *
929 v +
930 L0[i] * L1[j] * L2[k] * diff_v[dd];
931 }
932 ++jj;
933 }
934 }
935
936 const int i = oo;
937 const int j = oo;
938 const int k = oo;
939 for (dd = 0; dd != 3; ++dd) {
940 diff_volumeN[3 * shift + 3 * jj + dd] =
941 (diffL0[dd * (p + 1) + i] * L1[j] * L2[k] +
942 L0[i] * diffL1[dd * (p + 1) + j] * L2[k] +
943 L0[i] * L1[j] * diffL2[dd * (p + 1) + k]) *
944 v +
945 L0[i] * L1[j] * L2[k] * diff_v[dd];
946 }
947
948 ++jj;
949 }
950 if (jj != P)
951 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
952 "Inconsistent implementation (bug in the code) %d != %d", jj,
953 P);
954 }
955 }
956 MoFEMFunctionReturnHot(0);
957}
958
959PetscErrorCode H1_QuadShapeFunctions_MBQUAD(
960 int *faces_nodes, int p, double *N, double *diffN, double *faceN,
961 double *diff_faceN, int GDIM,
962 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
963 double *L, double *diffL,
964 const int dim)) {
965 MoFEMFunctionBeginHot;
966 // TODO: return separately components of the tensor product between two edges
967
968 int P = NBFACEQUAD_H1(p);
969 if (P == 0)
970 MoFEMFunctionReturnHot(0);
971 double diff_ksiL0F0[2], diff_ksiL1F0[2];
972 double *diff_ksi_faces[] = {diff_ksiL0F0, diff_ksiL1F0};
973 int ii = 0;
974 for (; ii < GDIM; ii++) {
975 int node_shift = ii * 4;
976 int node_diff_shift = 2 * node_shift;
977
978 int n0 = faces_nodes[0];
979 int n1 = faces_nodes[1];
980 int n2 = faces_nodes[2];
981 int n3 = faces_nodes[3];
982 double ksi0 = N[node_shift + n0] + N[node_shift + n3];
983 double ksi1 = N[node_shift + n1] + N[node_shift + n2];
984 double eta0 = N[node_shift + n0] + N[node_shift + n1];
985 double eta1 = N[node_shift + n2] + N[node_shift + n3];
986 double ksi_faces = ksi1 - ksi0;
987 double eta_faces = eta1 - eta0;
988 double diff_ksi_faces[2];
989 double diff_eta_faces[2];
990 int dd = 0;
991 for (; dd < 2; dd++) {
992 double diff_ksi0 = diffN[node_diff_shift + 2 * n0 + dd] +
993 diffN[node_diff_shift + 2 * n3 + dd];
994 double diff_ksi1 = diffN[node_diff_shift + 2 * n1 + dd] +
995 diffN[node_diff_shift + 2 * n2 + dd];
996 double diff_eta0 = diffN[node_diff_shift + 2 * n0 + dd] +
997 diffN[node_diff_shift + 2 * n1 + dd];
998 double diff_eta1 = diffN[node_diff_shift + 2 * n2 + dd] +
999 diffN[node_diff_shift + 2 * n3 + dd];
1000
1001 diff_ksi_faces[dd] = diff_ksi1 - diff_ksi0;
1002 diff_eta_faces[dd] = diff_eta1 - diff_eta0;
1003 }
1004 double L0[p + 1], L1[p + 1];
1005 double diffL0[2 * (p + 1)], diffL1[2 * (p + 1)];
1006 ierr = base_polynomials(p, ksi_faces, diff_ksi_faces, L0, diffL0, 2);
1007 CHKERRQ(ierr);
1008 ierr = base_polynomials(p, eta_faces, diff_eta_faces, L1, diffL1, 2);
1009 CHKERRQ(ierr);
1010
1011 double v = N[node_shift + n0] * N[node_shift + n2] +
1012 N[node_shift + n1] * N[node_shift + n3];
1013
1014 double diff_v[2];
1015 dd = 0;
1016 for (; dd < 2; ++dd)
1017 diff_v[dd] =
1018
1019 diffN[node_diff_shift + 2 * n0 + dd] * N[node_shift + n2] +
1020 N[node_shift + n0] * diffN[node_diff_shift + 2 * n2 + dd] +
1021
1022 diffN[node_diff_shift + 2 * n1 + dd] * N[node_shift + n3] +
1023 N[node_shift + n1] * diffN[node_diff_shift + 2 * n3 + dd];
1024
1025 const int shift = ii * P;
1026
1027 if (faceN != NULL) {
1028 int jj = 0;
1029 int oo = 0;
1030 for (; oo <= p - 2; ++oo) {
1031 int pp0 = 0;
1032 for (; pp0 < oo; ++pp0) {
1033 int pp1 = oo;
1034 faceN[shift + jj] = L0[pp0] * L1[pp1] * v;
1035 ++jj;
1036 faceN[shift + jj] = L0[pp1] * L1[pp0] * v;
1037 ++jj;
1038 }
1039 faceN[shift + jj] = L0[oo] * L1[oo] * v;
1040 ++jj;
1041 }
1042 if (jj != P)
1043 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1044 "Inconsistent implementation (bug in the code) %d != %d", jj,
1045 P);
1046 }
1047
1048 if (diff_faceN != NULL) {
1049 int jj = 0;
1050 int oo = 0;
1051 for (; oo <= p - 2; ++oo) {
1052 int pp0 = 0;
1053 for (; pp0 < oo; ++pp0) {
1054 int pp1 = oo;
1055 int dd;
1056 for (dd = 0; dd < 2; dd++) {
1057 diff_faceN[2 * shift + 2 * jj + dd] =
1058 (L0[pp0] * diffL1[dd * (p + 1) + pp1] +
1059 diffL0[dd * (p + 1) + pp0] * L1[pp1]) *
1060 v +
1061 L0[pp0] * L1[pp1] * diff_v[dd];
1062 }
1063 ++jj;
1064 for (dd = 0; dd < 2; dd++) {
1065 diff_faceN[2 * shift + 2 * jj + dd] =
1066 (L0[pp1] * diffL1[dd * (p + 1) + pp0] +
1067 diffL0[dd * (p + 1) + pp1] * L1[pp0]) *
1068 v +
1069 L0[pp1] * L1[pp0] * diff_v[dd];
1070 }
1071 ++jj;
1072 }
1073 for (dd = 0; dd < 2; dd++) {
1074 diff_faceN[2 * shift + 2 * jj + dd] =
1075 (L0[oo] * diffL1[dd * (p + 1) + oo] +
1076 diffL0[dd * (p + 1) + oo] * L1[oo]) *
1077 v +
1078 L0[oo] * L1[oo] * diff_v[dd];
1079 }
1080 ++jj;
1081 }
1082 if (jj != P)
1083 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1084 "Inconsistent implementation (bug in the code) %d != %d", jj,
1085 P);
1086 }
1087 }
1088 MoFEMFunctionReturnHot(0);
1089}
1090
1091PetscErrorCode H1_EdgeShapeFunctions_MBQUAD(
1092 int *sense, int *p, double *N, double *diffN, double *edgeN[4],
1093 double *diff_edgeN[4], int GDIM,
1094 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
1095 double *L, double *diffL,
1096 const int dim)) {
1097 MoFEMFunctionBeginHot;
1098
1099 double *edgeN01 = NULL, *edgeN12 = NULL, *edgeN23 = NULL, *edgeN30 = NULL;
1100 if (edgeN != NULL) {
1101 edgeN01 = edgeN[0];
1102 edgeN12 = edgeN[1];
1103 edgeN23 = edgeN[2];
1104 edgeN30 = edgeN[3];
1105 }
1106 double *diff_edgeN01 = NULL, *diff_edgeN12 = NULL, *diff_edgeN23 = NULL,
1107 *diff_edgeN30 = NULL;
1108 if (diff_edgeN != NULL) {
1109 diff_edgeN01 = diff_edgeN[0];
1110 diff_edgeN12 = diff_edgeN[1];
1111 diff_edgeN23 = diff_edgeN[2];
1112 diff_edgeN30 = diff_edgeN[3];
1113 }
1114 int P[4];
1115 int ee = 0;
1116 for (; ee < 4; ee++)
1117 P[ee] = NBEDGE_H1(p[ee]);
1118
1119 int n0 = 0;
1120 int n1 = 1;
1121 int n2 = 2;
1122 int n3 = 3;
1123
1124 int ii = 0;
1125 for (; ii < GDIM; ii++) {
1126 int node_shift = ii * 4;
1127 int node_diff_shift = 2 * node_shift;
1128
1129 double shape0 = N[node_shift + n0];
1130 double shape1 = N[node_shift + n1];
1131 double shape2 = N[node_shift + n2];
1132 double shape3 = N[node_shift + n3];
1133
1134 double ksi01 = (shape1 + shape2 - shape0 - shape3) * sense[n0];
1135 double ksi12 = (shape2 + shape3 - shape1 - shape0) * sense[n1];
1136 double ksi23 = (shape3 + shape0 - shape2 - shape1) * sense[n2];
1137 double ksi30 = (shape0 + shape1 - shape3 - shape2) * sense[n3];
1138
1139 double extrude_zeta01 = shape0 + shape1;
1140 double extrude_ksi12 = shape1 + shape2;
1141 double extrude_zeta23 = shape2 + shape3;
1142 double extrude_ksi30 = shape0 + shape3;
1143
1144 double bubble_ksi = extrude_ksi12 * extrude_ksi30;
1145 double bubble_zeta = extrude_zeta01 * extrude_zeta23;
1146
1147 double diff_ksi01[2], diff_ksi12[2], diff_ksi23[2], diff_ksi30[2];
1148 double diff_extrude_zeta01[2];
1149 double diff_extrude_ksi12[2];
1150 double diff_extrude_zeta23[2];
1151 double diff_extrude_ksi30[2];
1152 double diff_bubble_ksi[2];
1153 double diff_bubble_zeta[2];
1154
1155 int d = 0;
1156 for (; d < 2; d++) {
1157 double diff_shape0 = diffN[node_diff_shift + 2 * n0 + d];
1158 double diff_shape1 = diffN[node_diff_shift + 2 * n1 + d];
1159 double diff_shape2 = diffN[node_diff_shift + 2 * n2 + d];
1160 double diff_shape3 = diffN[node_diff_shift + 2 * n3 + d];
1161 diff_ksi01[d] =
1162 (diff_shape1 + diff_shape2 - diff_shape0 - diff_shape3) * sense[n0];
1163 diff_ksi12[d] =
1164 (diff_shape2 + diff_shape3 - diff_shape1 - diff_shape0) * sense[n1];
1165 diff_ksi23[d] =
1166 (diff_shape3 + diff_shape0 - diff_shape2 - diff_shape1) * sense[n2];
1167 diff_ksi30[d] =
1168 (diff_shape0 + diff_shape1 - diff_shape3 - diff_shape2) * sense[n3];
1169 diff_extrude_zeta01[d] = diff_shape0 + diff_shape1;
1170 diff_extrude_ksi12[d] = diff_shape1 + diff_shape2;
1171 diff_extrude_zeta23[d] = diff_shape2 + diff_shape3;
1172 diff_extrude_ksi30[d] = diff_shape0 + diff_shape3;
1173 diff_bubble_ksi[d] = diff_extrude_ksi12[d] * extrude_ksi30 +
1174 extrude_ksi12 * diff_extrude_ksi30[d];
1175 diff_bubble_zeta[d] = diff_extrude_zeta01[d] * extrude_zeta23 +
1176 extrude_zeta01 * diff_extrude_zeta23[d];
1177 }
1178
1179 double L01[p[0] + 1], L12[p[1] + 1], L23[p[2] + 1], L30[p[3] + 1];
1180 double diffL01[2 * (p[0] + 1)], diffL12[2 * (p[1] + 1)],
1181 diffL23[2 * (p[2] + 1)], diffL30[2 * (p[3] + 1)];
1182 ierr = base_polynomials(p[0], ksi01, diff_ksi01, L01, diffL01, 2);
1183 CHKERRQ(ierr);
1184 ierr = base_polynomials(p[1], ksi12, diff_ksi12, L12, diffL12, 2);
1185 CHKERRQ(ierr);
1186 ierr = base_polynomials(p[2], ksi23, diff_ksi23, L23, diffL23, 2);
1187 CHKERRQ(ierr);
1188 ierr = base_polynomials(p[3], ksi30, diff_ksi30, L30, diffL30, 2);
1189 CHKERRQ(ierr);
1190
1191 int shift;
1192 if (edgeN != NULL) {
1193 // edge01
1194 shift = ii * (P[0]);
1195 cblas_dcopy(P[0], L01, 1, &edgeN01[shift], 1);
1196 cblas_dscal(P[0], bubble_ksi * extrude_zeta01, &edgeN01[shift], 1);
1197 // edge12
1198 shift = ii * (P[1]);
1199 cblas_dcopy(P[1], L12, 1, &edgeN12[shift], 1);
1200 cblas_dscal(P[1], bubble_zeta * extrude_ksi12, &edgeN12[shift], 1);
1201 // edge23
1202 shift = ii * (P[2]);
1203 cblas_dcopy(P[2], L23, 1, &edgeN23[shift], 1);
1204 cblas_dscal(P[2], bubble_ksi * extrude_zeta23, &edgeN23[shift], 1);
1205 // edge30
1206 shift = ii * (P[3]);
1207 cblas_dcopy(P[3], L30, 1, &edgeN30[shift], 1);
1208 cblas_dscal(P[3], bubble_zeta * extrude_ksi30, &edgeN30[shift], 1);
1209 }
1210 if (diff_edgeN != NULL) {
1211 if (P[0] > 0) {
1212 // edge01
1213 shift = ii * (P[0]);
1214 bzero(&diff_edgeN01[2 * shift], sizeof(double) * 2 * (P[0]));
1215 int d = 0;
1216 for (; d != 2; ++d) {
1217 cblas_daxpy(P[0], bubble_ksi * extrude_zeta01,
1218 &diffL01[d * (p[0] + 1)], 1, &diff_edgeN01[2 * shift + d],
1219 2);
1220 cblas_daxpy(P[0],
1221 diff_bubble_ksi[d] * extrude_zeta01 +
1222 bubble_ksi * diff_extrude_zeta01[d],
1223 L01, 1, &diff_edgeN01[2 * shift + d], 2);
1224 }
1225 }
1226 if (P[1] > 0) {
1227 // edge12
1228 shift = ii * (P[1]);
1229 bzero(&diff_edgeN12[2 * shift], sizeof(double) * 2 * (P[1]));
1230 int d = 0;
1231 for (; d != 2; ++d) {
1232 cblas_daxpy(P[1], bubble_zeta * extrude_ksi12,
1233 &diffL12[d * (p[1] + 1)], 1, &diff_edgeN12[2 * shift + d],
1234 2);
1235 cblas_daxpy(P[1],
1236 diff_bubble_zeta[d] * extrude_ksi12 +
1237 bubble_zeta * diff_extrude_ksi12[d],
1238 L12, 1, &diff_edgeN12[2 * shift + d], 2);
1239 }
1240 }
1241 if (P[2] > 0) {
1242 // edge23
1243 shift = ii * (P[2]);
1244 bzero(&diff_edgeN23[2 * shift], sizeof(double) * 2 * (P[2]));
1245 int d = 0;
1246 for (; d != 2; ++d) {
1247 cblas_daxpy(P[2], bubble_ksi * extrude_zeta23,
1248 &diffL23[d * (p[2] + 1)], 1, &diff_edgeN23[2 * shift + d],
1249 2);
1250 cblas_daxpy(P[2],
1251 diff_bubble_ksi[d] * extrude_zeta23 +
1252 bubble_ksi * diff_extrude_zeta23[d],
1253 L23, 1, &diff_edgeN23[2 * shift + d], 2);
1254 }
1255 }
1256 if (P[3] > 0) {
1257 // edge30
1258 shift = ii * (P[3]);
1259 bzero(&diff_edgeN30[2 * shift], sizeof(double) * 2 * (P[3]));
1260 int d = 0;
1261 for (; d != 2; ++d) {
1262 cblas_daxpy(P[3], bubble_zeta * extrude_ksi30,
1263 &diffL30[d * (p[3] + 1)], 1, &diff_edgeN30[2 * shift + d],
1264 2);
1265 cblas_daxpy(P[3],
1266 diff_bubble_zeta[d] * extrude_ksi30 +
1267 bubble_zeta * diff_extrude_ksi30[d],
1268 L30, 1, &diff_edgeN30[2 * shift + d], 2);
1269 }
1270 }
1271 }
1272 }
1273 MoFEMFunctionReturnHot(0);
1274}
L
static Index< 'L', 3 > L
Definition: BasicFeTools.hpp:95
p
static Index< 'p', 3 > p
Definition: BasicFeTools.hpp:90
invJac
MatrixDouble invJac
Definition: HookeElement.hpp:683
definitions.h
useful compiler directives and definitions
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
dim
const int dim
Definition: elec_phys_2D.cpp:12
F
@ F
Definition: free_surface.cpp:401
H1_QuadShapeFunctions_MBPRISM
PetscErrorCode H1_QuadShapeFunctions_MBPRISM(int *faces_nodes, int *p, double *N, double *diffN, double *faceN[], double *diff_faceN[], int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Definition: h1.c:642
H1_VolumeShapeDiffMBTETinvJ
PetscErrorCode H1_VolumeShapeDiffMBTETinvJ(int base_p, int p, double *volume_diffN, double *invJac, double *volume_diffNinvJac, int GDIM)
Definition: h1.c:592
H1_EdgeShapeFunctions_MBTET
PetscErrorCode H1_EdgeShapeFunctions_MBTET(int *sense, int *p, double *N, double *diffN, double *edgeN[], double *diff_edgeN[], int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Definition: h1.c:274
H1_FaceGradientOfDeformation_hierarchical
PetscErrorCode H1_FaceGradientOfDeformation_hierarchical(int p, double *diffN, double *dofs, double *F)
Definition: h1.c:620
ierr
static PetscErrorCode ierr
Definition: h1.c:15
H1_QuadShapeFunctions_MBQUAD
PetscErrorCode H1_QuadShapeFunctions_MBQUAD(int *faces_nodes, int p, double *N, double *diffN, double *faceN, double *diff_faceN, int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Definition: h1.c:959
H1_VolumeGradientOfDeformation_hierarchical
PetscErrorCode H1_VolumeGradientOfDeformation_hierarchical(int p, double *diffN, double *dofs, double *F)
Definition: h1.c:631
H1_FaceShapeDiffMBTETinvJ
PetscErrorCode H1_FaceShapeDiffMBTETinvJ(int *base_p, int *p, double *face_diffN[], double *invJac, double *face_diffNinvJac[], int GDIM)
Definition: h1.c:574
H1_EdgeShapeDiffMBTETinvJ
PetscErrorCode H1_EdgeShapeDiffMBTETinvJ(int *base_p, int *p, double *edge_diffN[], double *invJac, double *edge_diffNinvJac[], int GDIM)
Definition: h1.c:556
H1_EdgeShapeFunctions_MBQUAD
PetscErrorCode H1_EdgeShapeFunctions_MBQUAD(int *sense, int *p, double *N, double *diffN, double *edgeN[4], double *diff_edgeN[4], int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Definition: h1.c:1091
H1_EdgeShapeFunctions_MBTRI
PetscErrorCode H1_EdgeShapeFunctions_MBTRI(int *sense, int *p, double *N, double *diffN, double *edgeN[3], double *diff_edgeN[3], int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
H1_EdgeShapeFunctions_MBTRI.
Definition: h1.c:17
H1_FaceShapeFunctions_MBTRI
PetscErrorCode H1_FaceShapeFunctions_MBTRI(const int *face_nodes, int p, double *N, double *diffN, double *faceN, double *diff_faceN, int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Definition: h1.c:191
H1_EdgeGradientOfDeformation_hierarchical
PetscErrorCode H1_EdgeGradientOfDeformation_hierarchical(int p, double *diffN, double *dofs, double *F)
Definition: h1.c:609
H1_VolumeShapeFunctions_MBTET
PetscErrorCode H1_VolumeShapeFunctions_MBTET(int p, double *N, double *diffN, double *volumeN, double *diff_volumeN, int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Definition: h1.c:475
H1_FaceShapeFunctions_MBTET
PetscErrorCode H1_FaceShapeFunctions_MBTET(int *faces_nodes, int *p, double *N, double *diffN, double *faceN[], double *diff_faceN[], int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Definition: h1.c:373
H1_VolumeShapeFunctions_MBPRISM
PetscErrorCode H1_VolumeShapeFunctions_MBPRISM(int p, double *N, double *diffN, double *volumeN, double *diff_volumeN, int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Definition: h1.c:790
h1_hdiv_hcurl_l2.h
Functions to approximate hierarchical spaces.
NBVOLUMETET_H1
#define NBVOLUMETET_H1(P)
Number of base functions on tetrahedron for H1 space.
Definition: h1_hdiv_hcurl_l2.h:75
NBFACEQUAD_H1
#define NBFACEQUAD_H1(P)
Number of base functions on quad for H1 space.
Definition: h1_hdiv_hcurl_l2.h:65
NBVOLUMEPRISM_H1
#define NBVOLUMEPRISM_H1(P)
Number of base functions on prism for H1 space.
Definition: h1_hdiv_hcurl_l2.h:80
NBEDGE_H1
#define NBEDGE_H1(P)
Numer of base function on edge for H1 space.
Definition: h1_hdiv_hcurl_l2.h:55
NBFACETRI_H1
#define NBFACETRI_H1(P)
Number of base function on triangle for H1 space.
Definition: h1_hdiv_hcurl_l2.h:60
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
N
const int N
Definition: speed_test.cpp:3
Generated by Doxygen 1.9.5 and hosted at University of Glasgow