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