v0.14.0
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 
15 static PetscErrorCode ierr;
16 
17 PetscErrorCode 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 }
191 PetscErrorCode 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 }
274 PetscErrorCode 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 }
373 PetscErrorCode 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 }
475 PetscErrorCode 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 }
556 PetscErrorCode 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 }
574 PetscErrorCode 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 }
592 PetscErrorCode 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 }
609 PetscErrorCode 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 }
620 PetscErrorCode 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 }
631 PetscErrorCode 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 }
642 PetscErrorCode 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 }
790 PetscErrorCode 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 
959 PetscErrorCode 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 
1091 PetscErrorCode 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 }
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
NBEDGE_H1
#define NBEDGE_H1(P)
Numer of base function on edge for H1 space.
Definition: h1_hdiv_hcurl_l2.h:55
sdf_hertz.d
float d
Definition: sdf_hertz.py:5
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_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
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
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_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
NBVOLUMETET_H1
#define NBVOLUMETET_H1(P)
Number of base functions on tetrahedron for H1 space.
Definition: h1_hdiv_hcurl_l2.h:75
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_FaceShapeDiffMBTETinvJ
PetscErrorCode H1_FaceShapeDiffMBTETinvJ(int *base_p, int *p, double *face_diffN[], double *invJac, double *face_diffNinvJac[], int GDIM)
Definition: h1.c:574
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
EshelbianPlasticity::P
@ P
Definition: EshelbianContact.cpp:193
H1_FaceGradientOfDeformation_hierarchical
PetscErrorCode H1_FaceGradientOfDeformation_hierarchical(int p, double *diffN, double *dofs, double *F)
Definition: h1.c:620
MoFEM::L
VectorDouble L
Definition: Projection10NodeCoordsOnField.cpp:124
NBFACEQUAD_H1
#define NBFACEQUAD_H1(P)
Number of base functions on quad for H1 space.
Definition: h1_hdiv_hcurl_l2.h:65
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
definitions.h
useful compiler directives and definitions
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
N
const int N
Definition: speed_test.cpp:3
ierr
static PetscErrorCode ierr
Definition: h1.c:15
H1_VolumeShapeDiffMBTETinvJ
PetscErrorCode H1_VolumeShapeDiffMBTETinvJ(int base_p, int p, double *volume_diffN, double *invJac, double *volume_diffNinvJac, int GDIM)
Definition: h1.c:592
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
FTensor::dd
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
NBFACETRI_H1
#define NBFACETRI_H1(P)
Number of base function on triangle for H1 space.
Definition: h1_hdiv_hcurl_l2.h:60
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
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
h1_hdiv_hcurl_l2.h
Functions to approximate hierarchical spaces.
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
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_EdgeShapeDiffMBTETinvJ
PetscErrorCode H1_EdgeShapeDiffMBTETinvJ(int *base_p, int *p, double *edge_diffN[], double *invJac, double *edge_diffNinvJac[], int GDIM)
Definition: h1.c:556
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
NBVOLUMEPRISM_H1
#define NBVOLUMEPRISM_H1(P)
Number of base functions on prism for H1 space.
Definition: h1_hdiv_hcurl_l2.h:80
invJac
MatrixDouble invJac
Definition: HookeElement.hpp:683
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
F
@ F
Definition: free_surface.cpp:394
Generated by Doxygen 1.8.17 and hosted at University of Glasgow