v0.9.0
Hcurl.cpp
Go to the documentation of this file.
1 /** \file Hcurl.cpp
2 
3  \brief Implementation of H-curl base
4 
5  Hierarchic Finite Element Bases on Unstructured Tetrahedral
6  Meshes, by Mark Ainsworth and Joe Coyle and by Demkowicz
7  Shape functions for MBTRI/MBTET and HCurl space
8 
9 */
10 
11 // #undef GENERATE_VTK_WITH_CURL_BASE
12 #ifndef GENERATE_VTK_WITH_CURL_BASE
13 
14 using namespace MoFEM;
15 
17  int *sense, int *p, double *N, double *diffN, double *edge_n[],
18  double *diff_edge_n[], int nb_integration_pts,
19  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
20  double *L, double *diffL,
21  const int dim)) {
23 
24  const int edges_nodes[6][2] = {{0, 1}, {1, 2}, {2, 0},
25  {0, 3}, {1, 3}, {2, 3}};
26  int P[6];
27  for (int ee = 0; ee != 6; ee++)
28  P[ee] = NBEDGE_AINSWORTH_HCURL(p[ee]);
29 
32 
33  FTensor::Tensor1<double *, 3> t_node_diff_ksi[4] = {
34  FTensor::Tensor1<double *, 3>(&diffN[0], &diffN[1], &diffN[2]),
35  FTensor::Tensor1<double *, 3>(&diffN[3], &diffN[4], &diffN[5]),
36  FTensor::Tensor1<double *, 3>(&diffN[6], &diffN[7], &diffN[8]),
37  FTensor::Tensor1<double *, 3>(&diffN[9], &diffN[10], &diffN[11])};
38  double edge_diff_ksi[6][3];
39  FTensor::Tensor1<double *, 3> t_edge_diff_ksi[6] = {
40  FTensor::Tensor1<double *, 3>(&edge_diff_ksi[0][0], &edge_diff_ksi[0][1],
41  &edge_diff_ksi[0][2]),
42  FTensor::Tensor1<double *, 3>(&edge_diff_ksi[1][0], &edge_diff_ksi[1][1],
43  &edge_diff_ksi[1][2]),
44  FTensor::Tensor1<double *, 3>(&edge_diff_ksi[2][0], &edge_diff_ksi[2][1],
45  &edge_diff_ksi[2][2]),
46  FTensor::Tensor1<double *, 3>(&edge_diff_ksi[3][0], &edge_diff_ksi[3][1],
47  &edge_diff_ksi[3][2]),
48  FTensor::Tensor1<double *, 3>(&edge_diff_ksi[4][0], &edge_diff_ksi[4][1],
49  &edge_diff_ksi[4][2]),
50  FTensor::Tensor1<double *, 3>(&edge_diff_ksi[5][0], &edge_diff_ksi[5][1],
51  &edge_diff_ksi[5][2])};
52  for (int ee = 0; ee != 6; ee++) {
53  t_edge_diff_ksi[ee](i) = (t_node_diff_ksi[edges_nodes[ee][1]](i) -
54  t_node_diff_ksi[edges_nodes[ee][0]](i)) *
55  sense[ee];
56  }
57 
58  FTensor::Tensor1<double *, 3> t_edge_n[6] = {
59  FTensor::Tensor1<double *, 3>(&edge_n[0][0], &edge_n[0][1], &edge_n[0][2],
60  3),
61  FTensor::Tensor1<double *, 3>(&edge_n[1][0], &edge_n[1][1], &edge_n[1][2],
62  3),
63  FTensor::Tensor1<double *, 3>(&edge_n[2][0], &edge_n[2][1], &edge_n[2][2],
64  3),
65  FTensor::Tensor1<double *, 3>(&edge_n[3][0], &edge_n[3][1], &edge_n[3][2],
66  3),
67  FTensor::Tensor1<double *, 3>(&edge_n[4][0], &edge_n[4][1], &edge_n[4][2],
68  3),
69  FTensor::Tensor1<double *, 3>(&edge_n[5][0], &edge_n[5][1], &edge_n[5][2],
70  3)};
71  FTensor::Tensor2<double *, 3, 3> t_diff_edge_n[6] = {
73  &diff_edge_n[0][0], &diff_edge_n[0][3], &diff_edge_n[0][6],
74  &diff_edge_n[0][1], &diff_edge_n[0][4], &diff_edge_n[0][7],
75  &diff_edge_n[0][2], &diff_edge_n[0][5], &diff_edge_n[0][8], 9),
77  &diff_edge_n[1][0], &diff_edge_n[1][3], &diff_edge_n[1][6],
78  &diff_edge_n[1][1], &diff_edge_n[1][4], &diff_edge_n[1][7],
79  &diff_edge_n[1][2], &diff_edge_n[1][5], &diff_edge_n[1][8], 9),
81  &diff_edge_n[2][0], &diff_edge_n[2][3], &diff_edge_n[2][6],
82  &diff_edge_n[2][1], &diff_edge_n[2][4], &diff_edge_n[2][7],
83  &diff_edge_n[2][2], &diff_edge_n[2][5], &diff_edge_n[2][8], 9),
85  &diff_edge_n[3][0], &diff_edge_n[3][3], &diff_edge_n[3][6],
86  &diff_edge_n[3][1], &diff_edge_n[3][4], &diff_edge_n[3][7],
87  &diff_edge_n[3][2], &diff_edge_n[3][5], &diff_edge_n[3][8], 9),
89  &diff_edge_n[4][0], &diff_edge_n[4][3], &diff_edge_n[4][6],
90  &diff_edge_n[4][1], &diff_edge_n[4][4], &diff_edge_n[4][7],
91  &diff_edge_n[4][2], &diff_edge_n[4][5], &diff_edge_n[4][8], 9),
93  &diff_edge_n[5][0], &diff_edge_n[5][3], &diff_edge_n[5][6],
94  &diff_edge_n[5][1], &diff_edge_n[5][4], &diff_edge_n[5][7],
95  &diff_edge_n[5][2], &diff_edge_n[5][5], &diff_edge_n[5][8], 9)};
96  FTensor::Tensor1<double, 3> t_psi_e_0, t_psi_e_1;
97  FTensor::Tensor2<double, 3, 3> t_diff_psi_e_0, t_diff_psi_e_1;
98 
99  for (int ii = 0; ii != nb_integration_pts; ii++) {
100 
101  const int node_shift = ii * 4;
102  for (int ee = 0; ee != 6; ee++) {
103 
104  if (P[ee] == 0)
105  continue;
106 
107  t_psi_e_0(i) = (N[node_shift + edges_nodes[ee][1]] *
108  t_node_diff_ksi[edges_nodes[ee][0]](i) -
109  N[node_shift + edges_nodes[ee][0]] *
110  t_node_diff_ksi[edges_nodes[ee][1]](i)) *
111  sense[ee];
112  t_diff_psi_e_0(i, j) = (t_node_diff_ksi[edges_nodes[ee][1]](j) *
113  t_node_diff_ksi[edges_nodes[ee][0]](i) -
114  t_node_diff_ksi[edges_nodes[ee][0]](j) *
115  t_node_diff_ksi[edges_nodes[ee][1]](i)) *
116  sense[ee];
117 
118  t_psi_e_1(i) = N[node_shift + edges_nodes[ee][1]] *
119  t_node_diff_ksi[edges_nodes[ee][0]](i) +
120  N[node_shift + edges_nodes[ee][0]] *
121  t_node_diff_ksi[edges_nodes[ee][1]](i);
122  t_diff_psi_e_1(i, j) = t_node_diff_ksi[edges_nodes[ee][1]](j) *
123  t_node_diff_ksi[edges_nodes[ee][0]](i) +
124  t_node_diff_ksi[edges_nodes[ee][0]](j) *
125  t_node_diff_ksi[edges_nodes[ee][1]](i);
126 
127  (t_edge_n[ee])(i) = t_psi_e_0(i);
128  ++(t_edge_n[ee]);
129  (t_edge_n[ee])(i) = t_psi_e_1(i);
130  ++(t_edge_n[ee]);
131 
132  (t_diff_edge_n[ee])(i, j) = t_diff_psi_e_0(i, j);
133  ++(t_diff_edge_n[ee]);
134  (t_diff_edge_n[ee])(i, j) = t_diff_psi_e_1(i, j);
135  ++(t_diff_edge_n[ee]);
136 
137  if (p[ee] > 1) {
138 
139  const double ksi_0i = (N[node_shift + edges_nodes[ee][1]] -
140  N[node_shift + edges_nodes[ee][0]]) *
141  sense[ee];
142  double psi_l[p[ee] + 1], diff_psi_l[3 * p[ee] + 3];
143  CHKERR base_polynomials(p[ee], ksi_0i, &edge_diff_ksi[ee][0], psi_l,
144  diff_psi_l, 3);
145 
146  FTensor::Tensor1<double *, 3> t_diff_psi_l(
147  &diff_psi_l[0], &diff_psi_l[p[ee] + 1], &diff_psi_l[2 * p[ee] + 2],
148  1);
149 
150  for (int ll = 2; ll != P[ee]; ll++) {
151 
152  const double a = (double)(2 * ll + 1) / (double)(ll + 1);
153  const double b = (double)(ll) / (double)(ll + 1);
154 
155  (t_edge_n[ee])(i) = a * psi_l[ll - 1] * t_psi_e_1(i) -
156  b * psi_l[ll - 2] * t_psi_e_0(i);
157  ++(t_edge_n[ee]);
158 
159  (t_diff_edge_n[ee])(i, j) =
160  -b * (t_diff_psi_l(j) * t_psi_e_0(i) +
161  psi_l[ll - 2] * t_diff_psi_e_0(i, j));
162  ++t_diff_psi_l;
163  (t_diff_edge_n[ee])(i, j) +=
164  a * (t_diff_psi_l(j) * t_psi_e_1(i) +
165  psi_l[ll - 1] * t_diff_psi_e_1(i, j));
166  ++(t_diff_edge_n[ee]);
167  }
168  }
169  }
170  }
171 
173 }
174 
176  int sense, int p, double *N, double *diffN, double *edge_n,
177  double *diff_edge_n, int nb_integration_pts,
178  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
179  double *L, double *diffL,
180  const int dim)) {
182 
183  if (NBEDGE_AINSWORTH_HCURL(p) == 0)
185  if (diff_edge_n != NULL)
186  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
187  "Calculation of derivatives not implemented");
188 
190  FTensor::Tensor1<double, 3> t_node_diff_ksi[2];
191  t_node_diff_ksi[0](0) = diffN[0];
192  t_node_diff_ksi[0](1) = 0;
193  t_node_diff_ksi[0](2) = 0;
194  t_node_diff_ksi[1](0) = diffN[1];
195  t_node_diff_ksi[1](1) = 0;
196  t_node_diff_ksi[1](2) = 0;
197 
199  &edge_n[0], &edge_n[1], &edge_n[2]);
200  FTensor::Tensor1<double, 3> t_psi_e_0, t_psi_e_1;
201 
202  for (int ii = 0; ii != nb_integration_pts; ii++) {
203 
204  const int node_shift = ii * 2;
205 
206  t_psi_e_0(i) = (N[node_shift + 1] * t_node_diff_ksi[0](i) -
207  N[node_shift + 0] * t_node_diff_ksi[1](i)) *
208  sense;
209  t_psi_e_1(i) = N[node_shift + 1] * t_node_diff_ksi[0](i) +
210  N[node_shift + 0] * t_node_diff_ksi[1](i);
211 
212  t_edge_n(i) = t_psi_e_0(i);
213  ++t_edge_n;
214 
215  t_edge_n(i) = t_psi_e_1(i);
216  ++t_edge_n;
217 
218  if (p > 1) {
219 
220  const double ksi_0i = (N[node_shift + 1] - N[node_shift + 0]) * sense;
221  double psi_l[p + 1];
222  CHKERR base_polynomials(p, ksi_0i, NULL, psi_l, NULL, 3);
223 
224  for (int ll = 2; ll != NBEDGE_AINSWORTH_HCURL(p); ll++) {
225  const double a = (double)(2 * ll + 1) / (double)(ll + 1);
226  const double b = (double)(ll) / (double)(ll + 1);
227  t_edge_n(i) =
228  a * psi_l[ll - 1] * t_psi_e_1(i) - b * psi_l[ll - 2] * t_psi_e_0(i);
229  ++t_edge_n;
230  }
231  }
232  }
233 
235 }
236 
238  int *sense, int *p, double *N, double *diffN, double *edge_n[],
239  double *diff_edge_n[], int nb_integration_pts,
240  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
241  double *L, double *diffL,
242  const int dim)) {
243 
245 
246  // TODO This is not by atom tests properly
247 
248  const int edges_nodes[3][2] = {{0, 1}, {1, 2}, {2, 0}};
249  int P[3];
250  for (int ee = 0; ee < 3; ee++)
251  P[ee] = NBEDGE_AINSWORTH_HCURL(p[ee]);
252 
255 
256  FTensor::Tensor1<double, 3> t_node_diff_ksi[3] = {
257  FTensor::Tensor1<double, 3>(diffN[0], diffN[1], 0.),
258  FTensor::Tensor1<double, 3>(diffN[2], diffN[3], 0.),
259  FTensor::Tensor1<double, 3>(diffN[4], diffN[5], 0.),
260  };
261  FTensor::Tensor1<double, 2> t_2d_diff_ksi[3] = {
262  FTensor::Tensor1<double, 2>(diffN[0], diffN[1]),
263  FTensor::Tensor1<double, 2>(diffN[2], diffN[3]),
264  FTensor::Tensor1<double, 2>(diffN[4], diffN[5])};
265 
266  FTensor::Tensor1<double *, 3> t_edge_n[3] = {
267  FTensor::Tensor1<double *, 3>(&edge_n[0][0], &edge_n[0][1], &edge_n[0][2],
268  3),
269  FTensor::Tensor1<double *, 3>(&edge_n[1][0], &edge_n[1][1], &edge_n[1][2],
270  3),
271  FTensor::Tensor1<double *, 3>(&edge_n[2][0], &edge_n[2][1], &edge_n[2][2],
272  3)};
273  FTensor::Tensor2<FTensor::PackPtr<double *, 6>, 3, 2> t_diff_edge_n[3] = {
275  &diff_edge_n[0][HVEC0_0], &diff_edge_n[0][HVEC0_1],
276  &diff_edge_n[0][HVEC1_0], &diff_edge_n[0][HVEC1_1],
277  &diff_edge_n[0][HVEC2_0], &diff_edge_n[0][HVEC2_1]),
279  &diff_edge_n[1][HVEC0_0], &diff_edge_n[1][HVEC0_1],
280  &diff_edge_n[1][HVEC1_0], &diff_edge_n[1][HVEC1_1],
281  &diff_edge_n[1][HVEC2_0], &diff_edge_n[1][HVEC2_1]),
283  &diff_edge_n[2][HVEC0_0], &diff_edge_n[2][HVEC0_1],
284  &diff_edge_n[2][HVEC1_0], &diff_edge_n[2][HVEC1_1],
285  &diff_edge_n[2][HVEC2_0], &diff_edge_n[2][HVEC2_1])};
286 
287  FTensor::Tensor1<double, 3> t_psi_e_0, t_psi_e_1;
288  FTensor::Tensor2<double, 3, 2> t_diff_psi_e_0, t_diff_psi_e_1;
289 
290  for (int ee = 0; ee != 3; ee++) {
291 
292  if (P[ee] == 0)
293  continue;
294  const int node0 = edges_nodes[ee][0];
295  const int node1 = edges_nodes[ee][1];
296  const int sense_edge = sense[ee];
297 
298  t_diff_psi_e_0(i, j) =
299  (t_node_diff_ksi[node0](i) * t_2d_diff_ksi[node1](j) -
300  t_node_diff_ksi[node1](i) * t_2d_diff_ksi[node0](j)) *
301  sense_edge;
302  t_diff_psi_e_1(i, j) = t_node_diff_ksi[node0](i) * t_2d_diff_ksi[node1](j) +
303  t_node_diff_ksi[node1](i) * t_2d_diff_ksi[node0](j);
304 
305  for (int ii = 0; ii != nb_integration_pts; ii++) {
306 
307  const int node_shift = ii * 3;
308  const double n0 = N[node_shift + node0];
309  const double n1 = N[node_shift + node1];
310 
311  t_psi_e_0(i) =
312  (n1 * t_node_diff_ksi[node0](i) - n0 * t_node_diff_ksi[node1](i)) *
313  sense_edge;
314  t_psi_e_1(i) =
315  n1 * t_node_diff_ksi[node0](i) + n0 * t_node_diff_ksi[node1](i);
316 
317  (t_edge_n[ee])(i) = t_psi_e_0(i);
318  (t_diff_edge_n[ee])(i, j) = t_diff_psi_e_0(i, j);
319  ++(t_edge_n[ee]);
320  ++(t_diff_edge_n[ee]);
321  (t_edge_n[ee])(i) = t_psi_e_1(i);
322  (t_diff_edge_n[ee])(i, j) = t_diff_psi_e_1(i, j);
323  ++(t_edge_n[ee]);
324  ++(t_diff_edge_n[ee]);
325 
326  if (p[ee] > 1) {
327  const double ksi_0i = (n1 - n0) * sense_edge;
328  double diff_ksi_0i[] = {
329  ((t_2d_diff_ksi[node1])(0) - (t_2d_diff_ksi[node0])(0)) *
330  sense_edge,
331  ((t_2d_diff_ksi[node1])(1) - (t_2d_diff_ksi[node0])(1)) *
332  sense_edge};
333 
334  double psi_l[p[ee] + 1], diff_psi_l[2 * p[ee] + 2];
335  CHKERR
336  base_polynomials(p[ee], ksi_0i, diff_ksi_0i, psi_l, diff_psi_l, 2);
337 
338  FTensor::Tensor1<double *, 2> t_diff_psi_ll_m1(
339  &diff_psi_l[0 + 2 - 1], &diff_psi_l[p[ee] + 1 + 2 - 1], 1);
340  FTensor::Tensor1<double *, 2> t_diff_psi_ll_m2(
341  &diff_psi_l[0 + 2 - 2], &diff_psi_l[p[ee] + 1 + 2 - 2], 1);
342  for (int ll = 2; ll != P[ee]; ll++) {
343  const double a = (double)(2 * ll + 1) / (double)(ll + 1);
344  const double b = (double)(ll) / (double)(ll + 1);
345  (t_edge_n[ee])(i) = a * psi_l[ll - 1] * t_psi_e_1(i) -
346  b * psi_l[ll - 2] * t_psi_e_0(i);
347  (t_diff_edge_n[ee])(i, j) = a * t_psi_e_1(i) * t_diff_psi_ll_m1(j) +
348  a * psi_l[ll - 1] * t_diff_psi_e_1(i, j) -
349  b * t_psi_e_0(i) * t_diff_psi_ll_m2(j) -
350  b * psi_l[ll - 2] * t_diff_psi_e_0(i, j);
351  ++(t_edge_n[ee]);
352  ++(t_diff_edge_n[ee]);
353  ++t_diff_psi_ll_m1;
354  ++t_diff_psi_ll_m2;
355  }
356  }
357  }
358  }
359 
361 }
362 
364  int *faces_nodes, int *p, double *N, double *diffN, double *phi_f_e[4][3],
365  double *diff_phi_f_e[4][3], int nb_integration_pts,
366  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
367  double *L, double *diffL,
368  const int dim)) {
369 
371  const int edges[3][2] = {{0, 1}, {1, 2}, {2, 0}};
372 
375 
376  FTensor::Tensor1<double *, 3> t_node_diff_ksi[4] = {
377  FTensor::Tensor1<double *, 3>(&diffN[0], &diffN[1], &diffN[2]),
378  FTensor::Tensor1<double *, 3>(&diffN[3], &diffN[4], &diffN[5]),
379  FTensor::Tensor1<double *, 3>(&diffN[6], &diffN[7], &diffN[8]),
380  FTensor::Tensor1<double *, 3>(&diffN[9], &diffN[10], &diffN[11])};
381  FTensor::Tensor1<double, 3> t_edge_diff_ksi;
382  FTensor::Tensor1<double, 3> t_diff_beta_e;
383 
384  for (int ff = 0; ff != 4; ff++) {
385 
386  const int o_nodes[3] = {faces_nodes[3 * ff + 2], faces_nodes[3 * ff + 0],
387  faces_nodes[3 * ff + 1]};
388  FTensor::Tensor1<double *, 3> t_opposite_node_diff[3] = {
389  FTensor::Tensor1<double *, 3>(&diffN[3 * o_nodes[0] + 0],
390  &diffN[3 * o_nodes[0] + 1],
391  &diffN[3 * o_nodes[0] + 2]),
392  FTensor::Tensor1<double *, 3>(&diffN[3 * o_nodes[1] + 0],
393  &diffN[3 * o_nodes[1] + 1],
394  &diffN[3 * o_nodes[1] + 2]),
395  FTensor::Tensor1<double *, 3>(&diffN[3 * o_nodes[2] + 0],
396  &diffN[3 * o_nodes[2] + 1],
397  &diffN[3 * o_nodes[2] + 2])};
398  double psi_l[p[ff] + 1], diff_psi_l[3 * p[ff] + 3];
399 
400  const int nb_base_fun_on_face = NBFACETRI_AINSWORTH_EDGE_HCURL(p[ff]);
401  // cerr << nb_base_fun_on_face << " " << p[ff] << endl;
402  if (nb_base_fun_on_face == 0)
403  continue;
404 
405  for (int ee = 0; ee != 3; ee++) {
406 
407  FTensor::Tensor1<double *, 3> t_face_edge_base(
408  &phi_f_e[ff][ee][0], &phi_f_e[ff][ee][1], &phi_f_e[ff][ee][2], 3);
409  FTensor::Tensor2<double *, 3, 3> t_diff_face_edge_base(
410  &diff_phi_f_e[ff][ee][0], &diff_phi_f_e[ff][ee][3],
411  &diff_phi_f_e[ff][ee][6], &diff_phi_f_e[ff][ee][1],
412  &diff_phi_f_e[ff][ee][4], &diff_phi_f_e[ff][ee][7],
413  &diff_phi_f_e[ff][ee][2], &diff_phi_f_e[ff][ee][5],
414  &diff_phi_f_e[ff][ee][8], 9);
415  const int en[] = {faces_nodes[3 * ff + edges[ee][0]],
416  faces_nodes[3 * ff + edges[ee][1]]};
417  t_edge_diff_ksi(i) =
418  t_node_diff_ksi[en[1]](i) - t_node_diff_ksi[en[0]](i);
419 
420  for (int ii = 0; ii != nb_integration_pts; ii++) {
421 
422  const int node_shift = ii * 4;
423  const double n[] = {N[node_shift + faces_nodes[3 * ff + edges[ee][0]]],
424  N[node_shift + faces_nodes[3 * ff + edges[ee][1]]]};
425  const double ksi_0i = n[1] - n[0];
426  CHKERR base_polynomials(p[ff], ksi_0i, &t_edge_diff_ksi(0), psi_l,
427  diff_psi_l, 3);
428 
429  FTensor::Tensor1<double *, 3> t_diff_psi_l(
430  &diff_psi_l[0], &diff_psi_l[p[ff] + 1], &diff_psi_l[2 * p[ff] + 2],
431  1);
432 
433  const double beta_e = n[0] * n[1];
434  t_diff_beta_e(j) =
435  t_node_diff_ksi[en[0]](j) * n[1] + n[0] * t_node_diff_ksi[en[1]](j);
436 
437  for (int ll = 0; ll != nb_base_fun_on_face; ll++) {
438  // if(ll == nb_base_fun_on_face-1) cerr << psi_l[ll] << endl;
439 
440  t_face_edge_base(i) =
441  beta_e * psi_l[ll] * t_opposite_node_diff[ee](i);
442  ++t_face_edge_base;
443 
444  t_diff_face_edge_base(i, j) =
445  (beta_e * t_diff_psi_l(j) + t_diff_beta_e(j) * psi_l[ll]) *
446  t_opposite_node_diff[ee](i);
447 
448  ++t_diff_face_edge_base;
449  ++t_diff_psi_l;
450  }
451  }
452  }
453  }
454 
456 }
457 
459  int *faces_nodes, int p, double *N, double *diffN, double *phi_f_e[3],
460  double *diff_phi_f_e[3], int nb_integration_pts,
461  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
462  double *L, double *diffL,
463  const int dim)) {
464 
466 
467  const int nb_base_fun_on_face = NBFACETRI_AINSWORTH_EDGE_HCURL(p);
468  if (nb_base_fun_on_face == 0)
470 
471  const int edges[3][2] = {{0, 1}, {1, 2}, {2, 0}};
472 
475 
476  const int o_nodes[3] = {2, 0, 1};
477  FTensor::Tensor2<double, 3, 3> t_opposite_node_diff(
478  diffN[2 * o_nodes[0] + 0], diffN[2 * o_nodes[0] + 1], 0.,
479  diffN[2 * o_nodes[1] + 0], diffN[2 * o_nodes[1] + 1], 0.,
480  diffN[2 * o_nodes[2] + 0], diffN[2 * o_nodes[2] + 1], 0.);
481  double psi_l[p + 1];
482  double diff_psi_l[2 * p + 2];
483 
484  FTensor::Tensor1<double *, 3> t_face_edge_base[3] = {
485  FTensor::Tensor1<double *, 3>(&phi_f_e[0][HVEC0], &phi_f_e[0][HVEC1],
486  &phi_f_e[0][HVEC2], 3),
487  FTensor::Tensor1<double *, 3>(&phi_f_e[1][HVEC0], &phi_f_e[1][HVEC1],
488  &phi_f_e[1][HVEC2], 3),
489  FTensor::Tensor1<double *, 3>(&phi_f_e[2][HVEC0], &phi_f_e[2][HVEC1],
490  &phi_f_e[2][HVEC2], 3),
491  };
493  t_diff_face_edge_base[3] = {
495  &diff_phi_f_e[0][HVEC0_0], &diff_phi_f_e[0][HVEC0_1],
496  &diff_phi_f_e[0][HVEC1_0], &diff_phi_f_e[0][HVEC1_1],
497  &diff_phi_f_e[0][HVEC2_0], &diff_phi_f_e[0][HVEC2_1]),
499  &diff_phi_f_e[1][HVEC0_0], &diff_phi_f_e[1][HVEC0_1],
500  &diff_phi_f_e[1][HVEC1_0], &diff_phi_f_e[1][HVEC1_1],
501  &diff_phi_f_e[1][HVEC2_0], &diff_phi_f_e[1][HVEC2_1]),
503  &diff_phi_f_e[2][HVEC0_0], &diff_phi_f_e[2][HVEC0_1],
504  &diff_phi_f_e[2][HVEC1_0], &diff_phi_f_e[2][HVEC1_1],
505  &diff_phi_f_e[2][HVEC2_0], &diff_phi_f_e[2][HVEC2_1])};
506 
507  for (int ee = 0; ee != 3; ee++) {
508 
509  const int node0 = faces_nodes[edges[ee][0]];
510  const int node1 = faces_nodes[edges[ee][1]];
511  double diff_ksi0i[] = {diffN[2 * node1 + 0] - diffN[2 * node0 + 0],
512  diffN[2 * node1 + 1] - diffN[2 * node0 + 1]};
513 
514  for (int ii = 0; ii != nb_integration_pts; ii++) {
515 
516  const int node_shift = ii * 3;
517  const double n0 = N[node_shift + node0];
518  const double n1 = N[node_shift + node1];
519  const double ksi_0i = n1 - n0;
520  CHKERR base_polynomials(p, ksi_0i, diff_ksi0i, psi_l, diff_psi_l, 2);
521 
522  const double beta_e = n0 * n1;
523  FTensor::Tensor1<double, 2> t_diff_beta_e(
524  diffN[2 * node0 + 0] * n1 + n0 * diffN[2 * node1 + 0],
525  diffN[2 * node0 + 1] * n1 + n0 * diffN[2 * node1 + 1]);
526  FTensor::Tensor1<double *, 2> t_diff_psi_l(&diff_psi_l[0],
527  &diff_psi_l[p + 1], 1);
528 
529  for (int ll = 0; ll != nb_base_fun_on_face; ll++) {
530  t_face_edge_base[ee](i) =
531  beta_e * psi_l[ll] * t_opposite_node_diff(ee, i);
532  t_diff_face_edge_base[ee](i, j) =
533  beta_e * t_opposite_node_diff(ee, i) * t_diff_psi_l(j) +
534  psi_l[ll] * t_opposite_node_diff(ee, i) * t_diff_beta_e(j);
535  ++t_face_edge_base[ee];
536  ++t_diff_face_edge_base[ee];
537  ++t_diff_psi_l;
538  }
539  }
540  }
541 
543 }
544 
546  int *faces_nodes, int *p, double *N, double *diffN, double *phi_f[4],
547  double *diff_phi_f[4], int nb_integration_pts,
548  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
549  double *L, double *diffL,
550  const int dim)) {
551 
553 
556 
557  // double coords[] = { 0,0,0, 1,0,0, 0,1,0, 0,0,1 };
558  // FTensor::Tensor1<double*,3> t_coords[4] = {
559  // FTensor::Tensor1<double*,3>(&coords[0],&coords[ 1],&coords[ 2]),
560  // FTensor::Tensor1<double*,3>(&coords[3],&coords[ 4],&coords[ 5]),
561  // FTensor::Tensor1<double*,3>(&coords[6],&coords[ 7],&coords[ 8]),
562  // FTensor::Tensor1<double*,3>(&coords[9],&coords[10],&coords[11])
563  // };
564  FTensor::Tensor1<double *, 3> t_node_diff_ksi[4] = {
565  FTensor::Tensor1<double *, 3>(&diffN[0], &diffN[1], &diffN[2]),
566  FTensor::Tensor1<double *, 3>(&diffN[3], &diffN[4], &diffN[5]),
567  FTensor::Tensor1<double *, 3>(&diffN[6], &diffN[7], &diffN[8]),
568  FTensor::Tensor1<double *, 3>(&diffN[9], &diffN[10], &diffN[11])};
569  FTensor::Tensor1<double, 3> t_diff_ksi0i, t_diff_ksi0j;
570  FTensor::Tensor1<double, 3> diff_beta_0ij;
571 
574 
575  for (int ff = 0; ff != 4; ff++) {
576 
577  if (NBFACETRI_AINSWORTH_FACE_HCURL(p[ff]) == 0)
578  continue;
579 
580  int n0 = faces_nodes[3 * ff + 0];
581  int n1 = faces_nodes[3 * ff + 1];
582  int n2 = faces_nodes[3 * ff + 2];
583 
584  // tou_0i(i) = t_coords[n1](i)-t_coords[n0](i);
585  // tou_0j(i) = t_coords[n2](i)-t_coords[n0](i);
586  tou_0i(i) = t_node_diff_ksi[n1](i) - t_node_diff_ksi[n0](i);
587  tou_0j(i) = t_node_diff_ksi[n2](i) - t_node_diff_ksi[n0](i);
588 
589  t_diff_ksi0i(i) = t_node_diff_ksi[n1](i) - t_node_diff_ksi[n0](i);
590  t_diff_ksi0j(i) = t_node_diff_ksi[n2](i) - t_node_diff_ksi[n0](i);
591 
592  double psi_l_0i[p[ff] + 1], diff_psi_l_0i[3 * p[ff] + 3];
593  double psi_l_0j[p[ff] + 1], diff_psi_l_0j[3 * p[ff] + 3];
594 
595  FTensor::Tensor1<double *, 3> t_phi_f(&phi_f[ff][0], &phi_f[ff][1],
596  &phi_f[ff][2], 3);
598  &diff_phi_f[ff][0], &diff_phi_f[ff][3], &diff_phi_f[ff][6],
599  &diff_phi_f[ff][1], &diff_phi_f[ff][4], &diff_phi_f[ff][7],
600  &diff_phi_f[ff][2], &diff_phi_f[ff][5], &diff_phi_f[ff][8], 9);
602 
603  for (int ii = 0; ii != nb_integration_pts; ii++) {
604 
605  const int node_shift = ii * 4;
606  const double beta_0ij =
607  N[node_shift + n0] * N[node_shift + n1] * N[node_shift + n2];
608  diff_beta_0ij(i) =
609  t_node_diff_ksi[n0](i) * N[node_shift + n1] * N[node_shift + n2] +
610  N[node_shift + n0] * t_node_diff_ksi[n1](i) * N[node_shift + n2] +
611  N[node_shift + n0] * N[node_shift + n1] * t_node_diff_ksi[n2](i);
612 
613  const double ksi_0i = N[node_shift + n1] - N[node_shift + n0];
614  CHKERR base_polynomials(p[ff], ksi_0i, &t_diff_ksi0i(0), psi_l_0i,
615  diff_psi_l_0i, 3);
616 
617  const double ksi_0j = N[node_shift + n2] - N[node_shift + n0];
618  CHKERR base_polynomials(p[ff], ksi_0j, &t_diff_ksi0j(0), psi_l_0j,
619  diff_psi_l_0j, 3);
620 
621  int cc = 0;
622  for (int oo = 0; oo <= (p[ff] - 3); oo++) {
623  FTensor::Tensor1<double *, 3> t_diff_psi_l_0i(
624  &diff_psi_l_0i[0], &diff_psi_l_0i[p[ff] + 1],
625  &diff_psi_l_0i[2 * p[ff] + 2], 1);
626  for (int pp0 = 0; pp0 <= oo; pp0++) {
627  const int pp1 = oo - pp0;
628  if (pp1 >= 0) {
629  FTensor::Tensor1<double *, 3> t_diff_psi_l_0j(
630  &diff_psi_l_0j[pp1], &diff_psi_l_0j[p[ff] + 1 + pp1],
631  &diff_psi_l_0j[2 * p[ff] + 2 + pp1], 1);
632  // base functions
633  const double a = beta_0ij * psi_l_0i[pp0] * psi_l_0j[pp1];
634  t_phi_f(i) = a * tou_0i(i);
635  ++t_phi_f;
636  ++cc;
637  t_phi_f(i) = a * tou_0j(i);
638  ++t_phi_f;
639  ++cc;
640  // diff base functions
641  t_b(j) = diff_beta_0ij(j) * psi_l_0i[pp0] * psi_l_0j[pp1] +
642  beta_0ij * t_diff_psi_l_0i(j) * psi_l_0j[pp1] +
643  beta_0ij * psi_l_0i[pp0] * t_diff_psi_l_0j(j);
644  t_diff_phi_f(i, j) = t_b(j) * tou_0i(i);
645  ++t_diff_phi_f;
646  t_diff_phi_f(i, j) = t_b(j) * tou_0j(i);
647  ++t_diff_phi_f;
648  ++t_diff_psi_l_0i;
649  }
650  }
651  }
652  const int nb_base_fun_on_face = NBFACETRI_AINSWORTH_FACE_HCURL(p[ff]);
653  if (cc != nb_base_fun_on_face) {
654  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
655  "Wrong number of base functions %d != %d", cc,
656  nb_base_fun_on_face);
657  }
658  }
659  }
661 }
662 
664  int *faces_nodes, int p, double *N, double *diffN, double *phi_f,
665  double *diff_phi_f, int nb_integration_pts,
666  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
667  double *L, double *diffL,
668  const int dim)) {
669 
671 
672  double zero = 0;
673  FTensor::Tensor2<double *, 3, 3> t_node_diff_ksi(&diffN[0], &diffN[1], &zero,
674  &diffN[2], &diffN[3], &zero,
675  &diffN[4], &diffN[5], &zero);
676 
679 
680  if (NBFACETRI_AINSWORTH_FACE_HCURL(p) == 0)
682 
686 
687  const int node0 = faces_nodes[0];
688  const int node1 = faces_nodes[1];
689  const int node2 = faces_nodes[2];
690 
693 
694  tou_0i(i) = t_node_diff_ksi(N1, i) - t_node_diff_ksi(N0, i);
695  ;
696  tou_0j(i) = t_node_diff_ksi(N2, i) - t_node_diff_ksi(N0, i);
697  ;
698 
699  double psi_l_0i[p + 1], psi_l_0j[p + 1];
700  double diff_psi_l_0i[2 * p + 2], diff_psi_l_0j[2 * p + 2];
701  FTensor::Tensor1<double *, 3> t_phi_f(&phi_f[0], &phi_f[1], &phi_f[2], 3);
703  &diff_phi_f[HVEC0_0], &diff_phi_f[HVEC0_1], &diff_phi_f[HVEC1_0],
704  &diff_phi_f[HVEC1_1], &diff_phi_f[HVEC2_0], &diff_phi_f[HVEC2_1]);
705 
706  double diff_ksi_0i[] = {t_node_diff_ksi(node1, 0) - t_node_diff_ksi(node0, 0),
707  t_node_diff_ksi(node1, 1) -
708  t_node_diff_ksi(node0, 1)};
709  double diff_ksi_0j[] = {t_node_diff_ksi(node2, 0) - t_node_diff_ksi(node0, 0),
710  t_node_diff_ksi(node2, 1) -
711  t_node_diff_ksi(node0, 1)};
712 
713  for (int ii = 0; ii != nb_integration_pts; ii++) {
714 
715  const int node_shift = ii * 3;
716  const double n0 = N[node_shift + node0];
717  const double n1 = N[node_shift + node1];
718  const double n2 = N[node_shift + node2];
719 
720  const double beta_0ij = n0 * n1 * n2;
721  FTensor::Tensor1<double, 2> diff_beta_0ij(
722  t_node_diff_ksi(node0, 0) * n1 * n2 +
723  n0 * t_node_diff_ksi(node1, 0) * n2 +
724  n0 * n1 * t_node_diff_ksi(node2, 0),
725  t_node_diff_ksi(node0, 1) * n1 * n2 +
726  n0 * t_node_diff_ksi(node1, 1) * n2 +
727  n0 * n1 * t_node_diff_ksi(node2, 1));
728 
729  const double ksi_0i = N[node_shift + node1] - N[node_shift + node0];
730  CHKERR base_polynomials(p, ksi_0i, diff_ksi_0i, psi_l_0i, diff_psi_l_0i, 2);
731 
732  const double ksi_0j = N[node_shift + node2] - N[node_shift + node0];
733  CHKERR base_polynomials(p, ksi_0j, diff_ksi_0j, psi_l_0j, diff_psi_l_0j, 2);
734 
735  int cc = 0;
737  for (int oo = 0; oo <= (p - 3); oo++) {
738  for (int pp0 = 0; pp0 <= oo; pp0++) {
739  const int pp1 = oo - pp0;
740  if (pp1 >= 0) {
741  FTensor::Tensor1<double, 2> t_diff_psi_l_0i(
742  diff_psi_l_0i[0 + pp0], diff_psi_l_0i[p + 1 + pp0]);
743  FTensor::Tensor1<double, 2> t_diff_psi_l_0j(
744  diff_psi_l_0j[0 + pp1], diff_psi_l_0j[p + 1 + pp1]);
745  const double a = beta_0ij * psi_l_0i[pp0] * psi_l_0j[pp1];
746  t_diff_a(j) = diff_beta_0ij(j) * psi_l_0i[pp0] * psi_l_0j[pp1] +
747  beta_0ij * psi_l_0i[pp0] * t_diff_psi_l_0j(j) +
748  beta_0ij * psi_l_0j[pp1] * t_diff_psi_l_0i(j);
749 
750  t_phi_f(i) = a * tou_0i(i);
751  t_diff_phi_f(i, j) = tou_0i(i) * t_diff_a(j);
752  ++t_phi_f;
753  ++t_diff_phi_f;
754  ++cc;
755  t_phi_f(i) = a * tou_0j(i);
756  t_diff_phi_f(i, j) = tou_0j(i) * t_diff_a(j);
757  ++t_phi_f;
758  ++t_diff_phi_f;
759  ++cc;
760  }
761  }
762  }
763 
764  const int nb_base_fun_on_face = NBFACETRI_AINSWORTH_FACE_HCURL(p);
765  if (cc != nb_base_fun_on_face) {
766  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
767  "Wrong number of base functions %d != %d", cc,
768  nb_base_fun_on_face);
769  }
770  }
771 
773 }
774 
776  int *faces_nodes, int p, double *N, double *diffN, double *phi_v,
777  double *diff_phi_v, int nb_integration_pts,
778  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
779  double *L, double *diffL,
780  const int dim)) {
781 
783 
786 
787  const int face_opposite_nodes[] = {2, 0, 1, 3};
788 
791 
792  FTensor::Tensor1<double *, 3> t_node_diff_ksi[4] = {
793  FTensor::Tensor1<double *, 3>(&diffN[0], &diffN[1], &diffN[2]),
794  FTensor::Tensor1<double *, 3>(&diffN[3], &diffN[4], &diffN[5]),
795  FTensor::Tensor1<double *, 3>(&diffN[6], &diffN[7], &diffN[8]),
796  FTensor::Tensor1<double *, 3>(&diffN[9], &diffN[10], &diffN[11])};
797  FTensor::Tensor1<double, 3> t_diff_ksi0i, t_diff_ksi0j;
798 
799  MatrixDouble m_psi_l_0i(4, p + 1);
800  MatrixDouble m_psi_l_0j(4, p + 1);
801  MatrixDouble m_diff_psi_l_0i(4, 3 * p + 3);
802  MatrixDouble m_diff_psi_l_0j(4, 3 * p + 3);
803 
804  double *psi_l_0i[] = {&m_psi_l_0i(0, 0), &m_psi_l_0i(1, 0), &m_psi_l_0i(2, 0),
805  &m_psi_l_0i(3, 0)};
806  double *psi_l_0j[] = {&m_psi_l_0j(0, 0), &m_psi_l_0j(1, 0), &m_psi_l_0j(2, 0),
807  &m_psi_l_0j(3, 0)};
808  double *diff_psi_l_0i[] = {&m_diff_psi_l_0i(0, 0), &m_diff_psi_l_0i(1, 0),
809  &m_diff_psi_l_0i(2, 0), &m_diff_psi_l_0i(3, 0)};
810  double *diff_psi_l_0j[] = {&m_diff_psi_l_0j(0, 0), &m_diff_psi_l_0j(1, 0),
811  &m_diff_psi_l_0j(2, 0), &m_diff_psi_l_0j(3, 0)};
812  double beta_f[4];
813 
814  FTensor::Tensor1<double, 3> t_diff_beta_f[4];
815 
816  FTensor::Tensor1<double *, 3> t_phi_v(&phi_v[0], &phi_v[1], &phi_v[2], 3);
818  &diff_phi_v[0], &diff_phi_v[3], &diff_phi_v[6], &diff_phi_v[1],
819  &diff_phi_v[4], &diff_phi_v[7], &diff_phi_v[2], &diff_phi_v[5],
820  &diff_phi_v[8], 9);
821 
822  for (int ii = 0; ii != nb_integration_pts; ii++) {
823 
824  for (int ff = 0; ff != 4; ff++) {
825 
826  t_diff_ksi0i(i) = t_node_diff_ksi[faces_nodes[3 * ff + 1]](i) -
827  t_node_diff_ksi[faces_nodes[3 * ff + 0]](i);
828  t_diff_ksi0j(i) = t_node_diff_ksi[faces_nodes[3 * ff + 2]](i) -
829  t_node_diff_ksi[faces_nodes[3 * ff + 0]](i);
830 
831  const int node_shift = ii * 4;
832 
833  beta_f[ff] = N[node_shift + faces_nodes[3 * ff + 0]] *
834  N[node_shift + faces_nodes[3 * ff + 1]] *
835  N[node_shift + faces_nodes[3 * ff + 2]];
836 
837  t_diff_beta_f[ff](j) = t_node_diff_ksi[faces_nodes[3 * ff + 0]](j) *
838  N[node_shift + faces_nodes[3 * ff + 1]] *
839  N[node_shift + faces_nodes[3 * ff + 2]] +
840  N[node_shift + faces_nodes[3 * ff + 0]] *
841  t_node_diff_ksi[faces_nodes[3 * ff + 1]](j) *
842  N[node_shift + faces_nodes[3 * ff + 2]] +
843  N[node_shift + faces_nodes[3 * ff + 0]] *
844  N[node_shift + faces_nodes[3 * ff + 1]] *
845  t_node_diff_ksi[faces_nodes[3 * ff + 2]](j);
846 
847  const double ksi_0i = N[node_shift + faces_nodes[3 * ff + 1]] -
848  N[node_shift + faces_nodes[3 * ff + 0]];
849  CHKERR base_polynomials(p, ksi_0i, &t_diff_ksi0i(0), psi_l_0i[ff],
850  diff_psi_l_0i[ff], 3);
851 
852  const double ksi_0j = N[node_shift + faces_nodes[3 * ff + 2]] -
853  N[node_shift + faces_nodes[3 * ff + 0]];
854  CHKERR base_polynomials(p, ksi_0j, &t_diff_ksi0j(0), psi_l_0j[ff],
855  diff_psi_l_0j[ff], 3);
856  }
857 
858  int cc = 0;
859  for (int oo = 0; oo <= (p - 3); oo++) {
860  FTensor::Tensor1<double *, 3> t_diff_psi_l_0i[] = {
861  FTensor::Tensor1<double *, 3>(&diff_psi_l_0i[0][0],
862  &diff_psi_l_0i[0][p + 1],
863  &diff_psi_l_0i[0][2 * p + 2], 1),
864  FTensor::Tensor1<double *, 3>(&diff_psi_l_0i[1][0],
865  &diff_psi_l_0i[1][p + 1],
866  &diff_psi_l_0i[1][2 * p + 2], 1),
867  FTensor::Tensor1<double *, 3>(&diff_psi_l_0i[2][0],
868  &diff_psi_l_0i[2][p + 1],
869  &diff_psi_l_0i[2][2 * p + 2], 1),
870  FTensor::Tensor1<double *, 3>(&diff_psi_l_0i[3][0],
871  &diff_psi_l_0i[3][p + 1],
872  &diff_psi_l_0i[3][2 * p + 2], 1),
873  };
874  for (int pp0 = 0; pp0 <= oo; pp0++) {
875  const int pp1 = oo - pp0;
876  if (pp1 >= 0) {
877  for (int ff = 0; ff != 4; ff++) {
878  FTensor::Tensor1<double *, 3> t_diff_psi_l_0j(
879  &m_diff_psi_l_0j(ff, pp1), &m_diff_psi_l_0j(ff, p + 1 + pp1),
880  &m_diff_psi_l_0j(ff, 2 * p + 2 + pp1), 1);
881  const double t = psi_l_0i[ff][pp0] * psi_l_0j[ff][pp1];
882  const double a = beta_f[ff] * t;
883  t_phi_v(i) = a * t_node_diff_ksi[face_opposite_nodes[ff]](i);
884  ++t_phi_v;
885  ++cc;
886  t_diff_phi_v(i, j) =
887  (t_diff_beta_f[ff](j) * t +
888  beta_f[ff] * t_diff_psi_l_0i[ff](j) * psi_l_0j[ff][pp1] +
889  beta_f[ff] * psi_l_0i[ff][pp0] * t_diff_psi_l_0j(j)) *
890  t_node_diff_ksi[face_opposite_nodes[ff]](i);
891  ++t_diff_phi_v;
892  ++t_diff_psi_l_0i[ff];
893  }
894  }
895  }
896  }
897 
898  const int nb_base_fun_on_face = NBVOLUMETET_AINSWORTH_FACE_HCURL(p);
899  if (cc != nb_base_fun_on_face) {
900  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
901  "Wrong number of base functions %d != %d", cc,
902  nb_base_fun_on_face);
903  }
904  }
905 
907 }
908 
910  int p, double *N, double *diffN, double *phi_v, double *diff_phi_v,
911  int nb_integration_pts,
912  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
913  double *L, double *diffL,
914  const int dim)) {
915 
917 
920 
926 
927  FTensor::Tensor1<double *, 3> t_node_diff_ksi[4] = {
928  FTensor::Tensor1<double *, 3>(&diffN[0], &diffN[1], &diffN[2]),
929  FTensor::Tensor1<double *, 3>(&diffN[3], &diffN[4], &diffN[5]),
930  FTensor::Tensor1<double *, 3>(&diffN[6], &diffN[7], &diffN[8]),
931  FTensor::Tensor1<double *, 3>(&diffN[9], &diffN[10], &diffN[11])};
932 
933  double diff_ksi0i[3], diff_ksi0j[3], diff_ksi0k[3];
934  FTensor::Tensor1<double *, 3> t_diff_ksi0i(diff_ksi0i, &diff_ksi0i[1],
935  &diff_ksi0i[2]);
936  FTensor::Tensor1<double *, 3> t_diff_ksi0j(diff_ksi0j, &diff_ksi0j[1],
937  &diff_ksi0j[2]);
938  FTensor::Tensor1<double *, 3> t_diff_ksi0k(diff_ksi0k, &diff_ksi0k[1],
939  &diff_ksi0k[2]);
940  t_diff_ksi0i(i) = t_node_diff_ksi[1](i) - t_node_diff_ksi[0](i);
941  t_diff_ksi0j(i) = t_node_diff_ksi[2](i) - t_node_diff_ksi[0](i);
942  t_diff_ksi0k(i) = t_node_diff_ksi[3](i) - t_node_diff_ksi[0](i);
943 
944  std::vector<double> v_psi_l_0i(p + 1), v_diff_psi_l_0i(3 * p + 3);
945  std::vector<double> v_psi_l_0j(p + 1), v_diff_psi_l_0j(3 * p + 3);
946  std::vector<double> v_psi_l_0k(p + 1), v_diff_psi_l_0k(3 * p + 3);
947  double *psi_l_0i = &*v_psi_l_0i.begin();
948  double *diff_psi_l_0i = &*v_diff_psi_l_0i.begin();
949  double *psi_l_0j = &*v_psi_l_0j.begin();
950  double *diff_psi_l_0j = &*v_diff_psi_l_0j.begin();
951  double *psi_l_0k = &*v_psi_l_0k.begin();
952  double *diff_psi_l_0k = &*v_diff_psi_l_0k.begin();
953 
954  FTensor::Tensor1<double *, 3> t_phi_v(&phi_v[0], &phi_v[1], &phi_v[2], 3);
956  &diff_phi_v[0], &diff_phi_v[3], &diff_phi_v[6], &diff_phi_v[1],
957  &diff_phi_v[4], &diff_phi_v[7], &diff_phi_v[2], &diff_phi_v[5],
958  &diff_phi_v[8], 9);
960 
961  for (int ii = 0; ii != nb_integration_pts; ii++) {
962 
963  const int node_shift = ii * 4;
964  const int n0 = node_shift + 0;
965  const int n1 = node_shift + 1;
966  const int n2 = node_shift + 2;
967  const int n3 = node_shift + 3;
968 
969  const double beta_v = N[n0] * N[n1] * N[n2] * N[n3];
970 
971  const double ksi_0i = N[n1] - N[n0];
972  CHKERR base_polynomials(p, ksi_0i, diff_ksi0i, psi_l_0i, diff_psi_l_0i, 3);
973 
974  const double ksi_0j = N[n2] - N[n0];
975  CHKERR base_polynomials(p, ksi_0j, diff_ksi0j, psi_l_0j, diff_psi_l_0j, 3);
976 
977  const double ksi_0k = N[n3] - N[n0];
978  CHKERR base_polynomials(p, ksi_0k, diff_ksi0k, psi_l_0k, diff_psi_l_0k, 3);
979 
980  FTensor::Tensor1<double, 3> t_diff_beta_v;
981  t_diff_beta_v(j) = t_node_diff_ksi[0](j) * N[n1] * N[n2] * N[n3] +
982  N[n0] * t_node_diff_ksi[1](j) * N[n2] * N[n3] +
983  N[n0] * N[n1] * t_node_diff_ksi[2](j) * N[n3] +
984  N[n0] * N[n1] * N[n2] * t_node_diff_ksi[3](j);
985 
986  int cc = 0;
987  for (int oo = 0; oo <= (p - 4); oo++) {
988  FTensor::Tensor1<double *, 3> t_diff_psi_l_0i(
989  &diff_psi_l_0i[0], &diff_psi_l_0i[p + 1], &diff_psi_l_0i[2 * p + 2],
990  1);
991  for (int pp0 = 0; pp0 <= oo; pp0++) {
992  FTensor::Tensor1<double *, 3> t_diff_psi_l_0j(
993  &diff_psi_l_0j[0], &diff_psi_l_0j[p + 1], &diff_psi_l_0j[2 * p + 2],
994  1);
995  for (int pp1 = 0; (pp0 + pp1) <= oo; pp1++) {
996  const int pp2 = oo - pp0 - pp1;
997  if (pp2 >= 0) {
998  FTensor::Tensor1<double *, 3> t_diff_psi_l_0k(
999  &diff_psi_l_0k[0 + pp2], &diff_psi_l_0k[p + 1 + pp2],
1000  &diff_psi_l_0k[2 * p + 2 + pp2], 1);
1001  const double t = psi_l_0i[pp0] * psi_l_0j[pp1] * psi_l_0k[pp2];
1002  const double a = beta_v * t;
1003  t_phi_v(0) = a;
1004  t_phi_v(1) = 0;
1005  t_phi_v(2) = 0;
1006  ++t_phi_v;
1007  ++cc;
1008  t_phi_v(0) = 0;
1009  t_phi_v(1) = a;
1010  t_phi_v(2) = 0;
1011  ++t_phi_v;
1012  ++cc;
1013  t_phi_v(0) = 0;
1014  t_phi_v(1) = 0;
1015  t_phi_v(2) = a;
1016  ++t_phi_v;
1017  ++cc;
1018  t_b(j) =
1019  t_diff_beta_v(j) * t +
1020  beta_v * (t_diff_psi_l_0i(j) * psi_l_0j[pp1] * psi_l_0k[pp2] +
1021  psi_l_0i[pp0] * t_diff_psi_l_0j(j) * psi_l_0k[pp2] +
1022  psi_l_0i[pp0] * psi_l_0j[pp1] * t_diff_psi_l_0k(j));
1023  t_diff_phi_v(N0, j) = t_b(j);
1024  t_diff_phi_v(N1, j) = 0;
1025  t_diff_phi_v(N2, j) = 0;
1026  ++t_diff_phi_v;
1027  t_diff_phi_v(N0, j) = 0;
1028  t_diff_phi_v(N1, j) = t_b(j);
1029  t_diff_phi_v(N2, j) = 0;
1030  ++t_diff_phi_v;
1031  t_diff_phi_v(N0, j) = 0;
1032  t_diff_phi_v(N1, j) = 0;
1033  t_diff_phi_v(N2, j) = t_b(j);
1034  ++t_diff_phi_v;
1035  }
1036  ++t_diff_psi_l_0j;
1037  }
1038  ++t_diff_psi_l_0i;
1039  }
1040  }
1041 
1042  const int nb_base_fun_on_face = NBVOLUMETET_AINSWORTH_TET_HCURL(p);
1043  if (cc != nb_base_fun_on_face) {
1044  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1045  "Wrong number of base functions %d != %d", cc,
1046  nb_base_fun_on_face);
1047  }
1048  }
1050 }
1051 
1053  int *face_nodes, int *p, double *N, double *diffN, double *phi_f[4],
1054  double *diff_phi_f[4], int nb_integration_pts,
1055  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
1056  double *L, double *diffL,
1057  const int dim)) {
1058 
1060 
1061  try {
1062 
1063  MatrixDouble base_face_edge_functions[4];
1064  MatrixDouble diff_base_face_edge_functions[4];
1065  double *phi_f_e[4][3];
1066  double *diff_phi_f_e[4][3];
1067  for (int ff = 0; ff != 4; ff++) {
1068  if (NBFACETRI_AINSWORTH_EDGE_HCURL(p[ff]) == 0) {
1069  for (int ee = 0; ee != 3; ee++) {
1070  phi_f_e[ff][ee] = NULL;
1071  diff_phi_f_e[ff][ee] = NULL;
1072  }
1073  } else {
1074  base_face_edge_functions[ff].resize(
1075  3, 3 * NBFACETRI_AINSWORTH_EDGE_HCURL(p[ff]) * nb_integration_pts);
1076  diff_base_face_edge_functions[ff].resize(
1077  3, 9 * NBFACETRI_AINSWORTH_EDGE_HCURL(p[ff]) * nb_integration_pts);
1078  // base_face_edge_functions[ff].clear();
1079  // diff_base_face_edge_functions[ff].clear();
1080  for (int ee = 0; ee != 3; ee++) {
1081  phi_f_e[ff][ee] = &base_face_edge_functions[ff](ee, 0);
1082  diff_phi_f_e[ff][ee] = &diff_base_face_edge_functions[ff](ee, 0);
1083  }
1084  }
1085  }
1087  face_nodes, p, N, diffN, phi_f_e, diff_phi_f_e, nb_integration_pts,
1088  base_polynomials);
1089 
1090  VectorDouble base_face_bubble_functions[4];
1091  VectorDouble diff_base_face_bubble_functions[4];
1092  double *phi_f_f[4];
1093  double *diff_phi_f_f[4];
1094  for (int ff = 0; ff != 4; ff++) {
1095  int nb_dofs = NBFACETRI_AINSWORTH_FACE_HCURL(p[ff]);
1096  if (nb_dofs == 0) {
1097  phi_f_f[ff] = NULL;
1098  diff_phi_f_f[ff] = NULL;
1099  } else {
1100  base_face_bubble_functions[ff].resize(3 * nb_dofs * nb_integration_pts,
1101  false);
1102  diff_base_face_bubble_functions[ff].resize(
1103  9 * nb_dofs * nb_integration_pts, false);
1104  phi_f_f[ff] = &*base_face_bubble_functions[ff].data().begin();
1105  diff_phi_f_f[ff] = &*diff_base_face_bubble_functions[ff].data().begin();
1106  }
1107  }
1109  face_nodes, p, N, diffN, phi_f_f, diff_phi_f_f, nb_integration_pts,
1110  base_polynomials);
1111 
1114 
1115  for (int ff = 0; ff != 4; ff++) {
1116 
1117  if (NBFACETRI_AINSWORTH_EDGE_HCURL(p[ff]) == 0)
1118  continue;
1119 
1120  FTensor::Tensor1<double *, 3> t_face_edge_base[] = {
1121  FTensor::Tensor1<double *, 3>(&phi_f_e[ff][0][0], &phi_f_e[ff][0][1],
1122  &phi_f_e[ff][0][2], 3),
1123  FTensor::Tensor1<double *, 3>(&phi_f_e[ff][1][0], &phi_f_e[ff][1][1],
1124  &phi_f_e[ff][1][2], 3),
1125  FTensor::Tensor1<double *, 3>(&phi_f_e[ff][2][0], &phi_f_e[ff][2][1],
1126  &phi_f_e[ff][2][2], 3)};
1127  FTensor::Tensor2<double *, 3, 3> t_diff_face_edge_base[] = {
1129  &diff_phi_f_e[ff][0][0], &diff_phi_f_e[ff][0][3],
1130  &diff_phi_f_e[ff][0][6], &diff_phi_f_e[ff][0][1],
1131  &diff_phi_f_e[ff][0][4], &diff_phi_f_e[ff][0][7],
1132  &diff_phi_f_e[ff][0][2], &diff_phi_f_e[ff][0][5],
1133  &diff_phi_f_e[ff][0][8], 9),
1135  &diff_phi_f_e[ff][1][0], &diff_phi_f_e[ff][1][3],
1136  &diff_phi_f_e[ff][1][6], &diff_phi_f_e[ff][1][1],
1137  &diff_phi_f_e[ff][1][4], &diff_phi_f_e[ff][1][7],
1138  &diff_phi_f_e[ff][1][2], &diff_phi_f_e[ff][1][5],
1139  &diff_phi_f_e[ff][1][8], 9),
1141  &diff_phi_f_e[ff][2][0], &diff_phi_f_e[ff][2][3],
1142  &diff_phi_f_e[ff][2][6], &diff_phi_f_e[ff][2][1],
1143  &diff_phi_f_e[ff][2][4], &diff_phi_f_e[ff][2][7],
1144  &diff_phi_f_e[ff][2][2], &diff_phi_f_e[ff][2][5],
1145  &diff_phi_f_e[ff][2][8], 9)};
1146 
1147  FTensor::Tensor1<double *, 3> t_face_base(&phi_f[ff][0], &phi_f[ff][1],
1148  &phi_f[ff][2], 3);
1149  FTensor::Tensor2<double *, 3, 3> t_diff_face_base(
1150  &diff_phi_f[ff][0], &diff_phi_f[ff][3], &diff_phi_f[ff][6],
1151  &diff_phi_f[ff][1], &diff_phi_f[ff][4], &diff_phi_f[ff][7],
1152  &diff_phi_f[ff][2], &diff_phi_f[ff][5], &diff_phi_f[ff][8], 9);
1153 
1154  if (NBFACETRI_AINSWORTH_FACE_HCURL(p[ff]) > 0) {
1155  FTensor::Tensor1<double *, 3> t_face_face_base(
1156  &phi_f_f[ff][0], &phi_f_f[ff][1], &phi_f_f[ff][2], 3);
1157  FTensor::Tensor2<double *, 3, 3> t_diff_face_face_base(
1158  &diff_phi_f_f[ff][0], &diff_phi_f_f[ff][3], &diff_phi_f_f[ff][6],
1159  &diff_phi_f_f[ff][1], &diff_phi_f_f[ff][4], &diff_phi_f_f[ff][7],
1160  &diff_phi_f_f[ff][2], &diff_phi_f_f[ff][5], &diff_phi_f_f[ff][8],
1161  9);
1162  for (int ii = 0; ii != nb_integration_pts; ii++) {
1163  int cc = 0;
1164  for (int oo = 0; oo <= p[ff]; oo++) {
1165  // Face-edge base
1166  if (oo > 1) {
1167  for (int ee = 0; ee != 3; ee++) {
1168  for (int ll = NBFACETRI_AINSWORTH_EDGE_HCURL(oo - 1);
1169  ll != NBFACETRI_AINSWORTH_EDGE_HCURL(oo); ll++) {
1170  t_face_base(i) = t_face_edge_base[ee](i);
1171  ++cc;
1172  ++t_face_base;
1173  ++t_face_edge_base[ee];
1174  t_diff_face_base(i, j) = t_diff_face_edge_base[ee](i, j);
1175  ++t_diff_face_base;
1176  ++t_diff_face_edge_base[ee];
1177  // cerr << oo << " " << ll << " " << cc << " " <<
1178  // NBFACETRI_AINSWORTH_EDGE_HCURL(oo) << endl;
1179  }
1180  }
1181  }
1182  // Face-face base
1183  if (oo > 2) {
1184  for (int ll = NBFACETRI_AINSWORTH_FACE_HCURL(oo - 1);
1185  ll != NBFACETRI_AINSWORTH_FACE_HCURL(oo); ll++) {
1186  t_face_base(i) = t_face_face_base(i);
1187  ++cc;
1188  ++t_face_base;
1189  ++t_face_face_base;
1190  t_diff_face_base(i, j) = t_diff_face_face_base(i, j);
1191  ++t_diff_face_base;
1192  ++t_diff_face_face_base;
1193  }
1194  }
1195  }
1196  // check consistency
1197  const int nb_base_fun_on_face = NBFACETRI_AINSWORTH_HCURL(p[ff]);
1198  if (cc != nb_base_fun_on_face) {
1199  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1200  "Wrong number of base functions %d != %d", cc,
1201  nb_base_fun_on_face);
1202  }
1203  }
1204  } else {
1205  for (int ii = 0; ii != nb_integration_pts; ii++) {
1206  int cc = 0;
1207  for (int oo = 0; oo <= p[ff]; oo++) {
1208  // Face-edge base
1209  if (oo > 1) {
1210  for (int ee = 0; ee != 3; ee++) {
1211  for (int ll = NBFACETRI_AINSWORTH_EDGE_HCURL(oo - 1);
1212  ll != NBFACETRI_AINSWORTH_EDGE_HCURL(oo); ll++) {
1213  t_face_base(i) = t_face_edge_base[ee](i);
1214  ++cc;
1215  ++t_face_base;
1216  ++t_face_edge_base[ee];
1217  t_diff_face_base(i, j) = t_diff_face_edge_base[ee](i, j);
1218  ++t_diff_face_base;
1219  ++t_diff_face_edge_base[ee];
1220  // cerr << oo << " " << ll << " " << cc << " " <<
1221  // NBFACETRI_AINSWORTH_EDGE_HCURL(oo) << endl;
1222  }
1223  }
1224  }
1225  }
1226  // check consistency
1227  const int nb_base_fun_on_face = NBFACETRI_AINSWORTH_HCURL(p[ff]);
1228  if (cc != nb_base_fun_on_face) {
1229  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1230  "Wrong number of base functions %d != %d", cc,
1231  nb_base_fun_on_face);
1232  }
1233  }
1234  }
1235  }
1236 
1237  } catch (MoFEMException const &e) {
1238  SETERRQ(PETSC_COMM_SELF, e.errorCode, e.errorMessage);
1239  } catch (std::exception &ex) {
1240  std::ostringstream ss;
1241  ss << "thorw in method: " << ex.what() << " at line " << __LINE__
1242  << " in file " << __FILE__;
1243  SETERRQ(PETSC_COMM_SELF, MOFEM_STD_EXCEPTION_THROW, ss.str().c_str());
1244  }
1245 
1247 }
1248 
1250  int *faces_nodes, int p, double *N, double *diffN, double *phi_f,
1251  double *diff_phi_f, int nb_integration_pts,
1252  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
1253  double *L, double *diffL,
1254  const int dim)) {
1255 
1257 
1258  if (NBFACETRI_AINSWORTH_EDGE_HCURL(p) == 0)
1260 
1261  MatrixDouble base_face_edge_functions, diff_base_face_edge_functions;
1262  double *phi_f_e[3];
1263  double *diff_phi_f_e[3];
1264  base_face_edge_functions.resize(3, 3 * NBFACETRI_AINSWORTH_EDGE_HCURL(p) *
1265  nb_integration_pts);
1266  diff_base_face_edge_functions.resize(
1267  3, 2 * 3 * NBFACETRI_AINSWORTH_EDGE_HCURL(p) * nb_integration_pts);
1268  // base_face_edge_functions.clear();
1269  for (int ee = 0; ee != 3; ee++) {
1270  phi_f_e[ee] = &base_face_edge_functions(ee, 0);
1271  diff_phi_f_e[ee] = &diff_base_face_edge_functions(ee, 0);
1272  }
1274  faces_nodes, p, N, diffN, phi_f_e, diff_phi_f_e, nb_integration_pts,
1275  base_polynomials);
1276 
1277  VectorDouble base_face_bubble_functions;
1278  VectorDouble diff_base_face_bubble_functions;
1279  double *phi_f_f, *diff_phi_f_f;
1280  base_face_bubble_functions.resize(3 * NBFACETRI_AINSWORTH_FACE_HCURL(p) *
1281  nb_integration_pts);
1282  diff_base_face_bubble_functions.resize(
1283  2 * 3 * NBFACETRI_AINSWORTH_FACE_HCURL(p) * nb_integration_pts);
1284  phi_f_f = &*base_face_bubble_functions.data().begin();
1285  diff_phi_f_f = &*diff_base_face_bubble_functions.data().begin();
1287  faces_nodes, p, N, diffN, phi_f_f, diff_phi_f_f, nb_integration_pts,
1288  base_polynomials);
1289 
1290  // cerr << diff_base_face_bubble_functions << endl;
1291 
1294 
1295  FTensor::Tensor1<double *, 3> t_face_edge_base[] = {
1296  FTensor::Tensor1<double *, 3>(&phi_f_e[0][HVEC0], &phi_f_e[0][HVEC1],
1297  &phi_f_e[0][HVEC2], 3),
1298  FTensor::Tensor1<double *, 3>(&phi_f_e[1][HVEC0], &phi_f_e[1][HVEC1],
1299  &phi_f_e[1][HVEC2], 3),
1300  FTensor::Tensor1<double *, 3>(&phi_f_e[2][HVEC0], &phi_f_e[2][HVEC1],
1301  &phi_f_e[2][HVEC2], 3)};
1303  t_diff_face_edge_base[] = {
1305  &diff_phi_f_e[0][HVEC0_0], &diff_phi_f_e[0][HVEC0_1],
1306  &diff_phi_f_e[0][HVEC1_0], &diff_phi_f_e[0][HVEC1_1],
1307  &diff_phi_f_e[0][HVEC2_0], &diff_phi_f_e[0][HVEC2_1]),
1309  &diff_phi_f_e[1][HVEC0_0], &diff_phi_f_e[1][HVEC0_1],
1310  &diff_phi_f_e[1][HVEC1_0], &diff_phi_f_e[1][HVEC1_1],
1311  &diff_phi_f_e[1][HVEC2_0], &diff_phi_f_e[1][HVEC2_1]),
1313  &diff_phi_f_e[2][HVEC0_0], &diff_phi_f_e[2][HVEC0_1],
1314  &diff_phi_f_e[2][HVEC1_0], &diff_phi_f_e[2][HVEC1_1],
1315  &diff_phi_f_e[2][HVEC2_0], &diff_phi_f_e[2][HVEC2_1])};
1316 
1317  FTensor::Tensor1<double *, 3> t_face_base(&phi_f[0], &phi_f[1], &phi_f[2], 3);
1318  FTensor::Tensor2<FTensor::PackPtr<double *, 6>, 3, 2> t_diff_face_base(
1319  &diff_phi_f[HVEC0_0], &diff_phi_f[HVEC0_1], &diff_phi_f[HVEC1_0],
1320  &diff_phi_f[HVEC1_1], &diff_phi_f[HVEC2_0], &diff_phi_f[HVEC2_1]);
1321 
1322  if (NBFACETRI_AINSWORTH_FACE_HCURL(p) > 0) {
1323  FTensor::Tensor1<double *, 3> t_face_face_base(
1324  &phi_f_f[HVEC0], &phi_f_f[HVEC1], &phi_f_f[HVEC2], 3);
1325  FTensor::Tensor2<FTensor::PackPtr<double *, 6>, 3, 2> t_diff_face_face_base(
1326  &diff_phi_f_f[HVEC0_0], &diff_phi_f_f[HVEC0_1], &diff_phi_f_f[HVEC1_0],
1327  &diff_phi_f_f[HVEC1_1], &diff_phi_f_f[HVEC2_0], &diff_phi_f_f[HVEC2_1]);
1328  for (int ii = 0; ii != nb_integration_pts; ii++) {
1329  int cc = 0;
1330  for (int oo = 0; oo <= p; oo++) {
1331  // Face-edge base
1332  if (oo > 1) {
1333  for (int ee = 0; ee != 3; ee++) {
1334  for (int ll = NBFACETRI_AINSWORTH_EDGE_HCURL(oo - 1);
1335  ll != NBFACETRI_AINSWORTH_EDGE_HCURL(oo); ll++) {
1336  t_face_base(i) = t_face_edge_base[ee](i);
1337  t_diff_face_base(i, j) = t_diff_face_edge_base[ee](i, j);
1338  ++cc;
1339  ++t_face_base;
1340  ++t_face_edge_base[ee];
1341  ++t_diff_face_base;
1342  ++t_diff_face_edge_base[ee];
1343  }
1344  }
1345  }
1346  // Face-face base
1347  if (oo > 2) {
1348  for (int ll = NBFACETRI_AINSWORTH_FACE_HCURL(oo - 1);
1349  ll != NBFACETRI_AINSWORTH_FACE_HCURL(oo); ll++) {
1350  t_face_base(i) = t_face_face_base(i);
1351  t_diff_face_base(i, j) = t_diff_face_face_base(i, j);
1352  ++cc;
1353  ++t_face_base;
1354  ++t_face_face_base;
1355  ++t_diff_face_base;
1356  ++t_diff_face_face_base;
1357  }
1358  }
1359  }
1360  // check consistency
1361  const int nb_base_fun_on_face = NBFACETRI_AINSWORTH_HCURL(p);
1362  if (cc != nb_base_fun_on_face) {
1363  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1364  "Wrong number of base functions %d != %d", cc,
1365  nb_base_fun_on_face);
1366  }
1367  }
1368  } else {
1369  for (int ii = 0; ii != nb_integration_pts; ii++) {
1370  int cc = 0;
1371  for (int oo = 0; oo <= p; oo++) {
1372  // Face-edge base
1373  if (oo > 1) {
1374  for (int ee = 0; ee != 3; ee++) {
1375  for (int ll = NBFACETRI_AINSWORTH_EDGE_HCURL(oo - 1);
1376  ll != NBFACETRI_AINSWORTH_EDGE_HCURL(oo); ll++) {
1377  t_face_base(i) = t_face_edge_base[ee](i);
1378  t_diff_face_base(i, j) = t_diff_face_edge_base[ee](i, j);
1379  ++cc;
1380  ++t_face_base;
1381  ++t_face_edge_base[ee];
1382  ++t_diff_face_base;
1383  ++t_diff_face_edge_base[ee];
1384  // cerr << oo << " " << ll << " " << cc << " " <<
1385  // NBFACETRI_AINSWORTH_EDGE_HCURL(oo) << endl;
1386  }
1387  }
1388  }
1389  }
1390  // check consistency
1391  const int nb_base_fun_on_face = NBFACETRI_AINSWORTH_HCURL(p);
1392  if (cc != nb_base_fun_on_face) {
1393  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1394  "Wrong number of base functions %d != %d", cc,
1395  nb_base_fun_on_face);
1396  }
1397  }
1398  }
1399 
1401 }
1402 
1404  int p, double *N, double *diffN, double *phi_v, double *diff_phi_v,
1405  int nb_integration_pts,
1406  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
1407  double *L, double *diffL,
1408  const int dim)) {
1409 
1411 
1412  VectorDouble base_face_inetrior_functions(
1413  3 * NBVOLUMETET_AINSWORTH_FACE_HCURL(p) * nb_integration_pts);
1414  VectorDouble diff_base_face_inetrior_functions(
1415  9 * NBVOLUMETET_AINSWORTH_FACE_HCURL(p) * nb_integration_pts);
1416  // base_face_inetrior_functions.clear();
1417  // diff_base_face_inetrior_functions.clear();
1418  double *phi_v_f = &*base_face_inetrior_functions.data().begin();
1419  double *diff_phi_v_f = &*diff_base_face_inetrior_functions.data().begin();
1420  int faces_nodes[] = {0, 1, 3, 1, 2, 3, 0, 2, 3, 0, 1, 2};
1422  faces_nodes, p, N, diffN, phi_v_f, diff_phi_v_f, nb_integration_pts,
1423  base_polynomials);
1424 
1425  VectorDouble base_interior_functions(3 * NBVOLUMETET_AINSWORTH_TET_HCURL(p) *
1426  nb_integration_pts);
1427  VectorDouble diff_base_interior_functions(
1428  9 * NBVOLUMETET_AINSWORTH_TET_HCURL(p) * nb_integration_pts);
1429  // base_interior_functions.clear();
1430  // diff_base_interior_functions.clear();
1431  double *phi_v_v = &*base_interior_functions.data().begin();
1432  double *diff_phi_v_v = &*diff_base_interior_functions.data().begin();
1434  p, N, diffN, phi_v_v, diff_phi_v_v, nb_integration_pts, base_polynomials);
1435 
1438 
1439  FTensor::Tensor1<double *, 3> t_face_interior(&phi_v_f[0], &phi_v_f[1],
1440  &phi_v_f[2], 3);
1441  FTensor::Tensor2<double *, 3, 3> t_diff_face_interior(
1442  &diff_phi_v_f[0], &diff_phi_v_f[3], &diff_phi_v_f[6], &diff_phi_v_f[1],
1443  &diff_phi_v_f[4], &diff_phi_v_f[7], &diff_phi_v_f[2], &diff_phi_v_f[5],
1444  &diff_phi_v_f[8], 9);
1445 
1446  FTensor::Tensor1<double *, 3> t_phi_v(&phi_v[0], &phi_v[1], &phi_v[2], 3);
1447  FTensor::Tensor2<double *, 3, 3> t_diff_phi_v(
1448  &diff_phi_v[0], &diff_phi_v[3], &diff_phi_v[6], &diff_phi_v[1],
1449  &diff_phi_v[4], &diff_phi_v[7], &diff_phi_v[2], &diff_phi_v[5],
1450  &diff_phi_v[8], 9);
1451 
1452  if (NBVOLUMETET_AINSWORTH_TET_HCURL(p) > 0) {
1453  FTensor::Tensor1<double *, 3> t_volume_interior(&phi_v_v[0], &phi_v_v[1],
1454  &phi_v_v[2], 3);
1455  FTensor::Tensor2<double *, 3, 3> t_diff_volume_interior(
1456  &diff_phi_v_v[0], &diff_phi_v_v[3], &diff_phi_v_v[6], &diff_phi_v_v[1],
1457  &diff_phi_v_v[4], &diff_phi_v_v[7], &diff_phi_v_v[2], &diff_phi_v_v[5],
1458  &diff_phi_v_v[8], 9);
1459  for (int ii = 0; ii != nb_integration_pts; ii++) {
1460  int cc = 0;
1461  for (int oo = 0; oo <= p; oo++) {
1462  for (int ll = NBVOLUMETET_AINSWORTH_FACE_HCURL(oo - 1);
1463  ll != NBVOLUMETET_AINSWORTH_FACE_HCURL(oo); ll++) {
1464  t_phi_v(i) = t_face_interior(i);
1465  ++t_phi_v;
1466  ++t_face_interior;
1467  ++cc;
1468  t_diff_phi_v(i, j) = t_diff_face_interior(i, j);
1469  ++t_diff_phi_v;
1470  ++t_diff_face_interior;
1471  }
1472  for (int ll = NBVOLUMETET_AINSWORTH_TET_HCURL(oo - 1);
1473  ll != NBVOLUMETET_AINSWORTH_TET_HCURL(oo); ll++) {
1474  t_phi_v(i) = t_volume_interior(i);
1475  ++t_phi_v;
1476  ++t_volume_interior;
1477  ++cc;
1478  t_diff_phi_v(i, j) = t_diff_volume_interior(i, j);
1479  ++t_diff_phi_v;
1480  ++t_diff_volume_interior;
1481  }
1482  }
1483  // check consistency
1484  const int nb_base_fun_on_face = NBVOLUMETET_AINSWORTH_HCURL(p);
1485  if (cc != nb_base_fun_on_face) {
1486  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1487  "Wrong number of base functions %d != %d", cc,
1488  nb_base_fun_on_face);
1489  }
1490  }
1491  } else {
1492  for (int ii = 0; ii != nb_integration_pts; ii++) {
1493  int cc = 0;
1494  for (int oo = 0; oo <= p; oo++) {
1495  for (int ll = NBVOLUMETET_AINSWORTH_FACE_HCURL(oo - 1);
1496  ll != NBVOLUMETET_AINSWORTH_FACE_HCURL(oo); ll++) {
1497  t_phi_v(i) = t_face_interior(i);
1498  ++t_phi_v;
1499  ++t_face_interior;
1500  ++cc;
1501  t_diff_phi_v(i, j) = t_diff_face_interior(i, j);
1502  ++t_diff_phi_v;
1503  ++t_diff_face_interior;
1504  }
1505  }
1506  // check consistency
1507  const int nb_base_fun_on_face = NBVOLUMETET_AINSWORTH_HCURL(p);
1508  if (cc != nb_base_fun_on_face) {
1509  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1510  "Wrong number of base functions %d != %d", cc,
1511  nb_base_fun_on_face);
1512  }
1513  }
1514  }
1515 
1517 }
1518 
1519 #endif // Not GENERATE_VTK_WITH_CURL_BASE
1520 
1521 #ifdef GENERATE_VTK_WITH_CURL_BASE
1522 
1523 #include <MoFEM.hpp>
1524 #include <Hcurl.hpp>
1525 using namespace MoFEM;
1526 using namespace boost::numeric;
1527 
1528 MoFEMErrorCode VTK_Ainsworth_Hcurl_MBTET(const string file_name) {
1530 
1531  double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
1532 
1533  moab::Core core_ref;
1534  moab::Interface &moab_ref = core_ref;
1535 
1536  EntityHandle nodes[4];
1537  for (int nn = 0; nn < 4; nn++) {
1538  CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
1539  }
1540  EntityHandle tet;
1541  CHKERR moab_ref.create_element(MBTET, nodes, 4, tet);
1542 
1543  MoFEM::Core m_core_ref(moab_ref, PETSC_COMM_SELF, -2);
1544  MoFEM::Interface &m_field_ref = m_core_ref;
1545 
1546  CHKERR m_field_ref.getInterface<BitRefManager>()->setBitRefLevelByDim(
1547  0, 3, BitRefLevel().set(0));
1548 
1549  const int max_level = 4;
1550  for (int ll = 0; ll != max_level; ll++) {
1551  Range edges;
1552  CHKERR m_field_ref.getInterface<BitRefManager>()
1553  ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
1554  BitRefLevel().set(), MBEDGE, edges);
1555  Range tets;
1556  CHKERR m_field_ref.getInterface<BitRefManager>()
1557  ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
1558  BitRefLevel(ll).set(), MBTET, tets);
1559  // refine mesh
1560  MeshRefinement *m_ref;
1561  CHKERR m_field_ref.getInterface(m_ref);
1562  CHKERR m_ref->add_vertices_in_the_middle_of_edges(
1563  edges, BitRefLevel().set(ll + 1));
1564  CHKERR m_ref->refine_TET(tets, BitRefLevel().set(ll + 1));
1565  }
1566 
1567  Range tets;
1568  CHKERR m_field_ref.getInterface<BitRefManager>()
1569  ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(max_level),
1570  BitRefLevel().set(max_level), MBTET, tets);
1571 
1572  // Use 10 node tets to print base
1573  if (1) {
1574 
1575  EntityHandle meshset;
1576  CHKERR moab_ref.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER, meshset);
1577  CHKERR moab_ref.add_entities(meshset, tets);
1578  CHKERR moab_ref.convert_entities(meshset, true, false, false);
1579  CHKERR moab_ref.delete_entities(&meshset, 1);
1580  }
1581 
1582  Range elem_nodes;
1583  CHKERR moab_ref.get_connectivity(tets, elem_nodes, false);
1584 
1585  const int nb_gauss_pts = elem_nodes.size();
1586  MatrixDouble gauss_pts(nb_gauss_pts, 4);
1587  gauss_pts.clear();
1588  Range::iterator nit = elem_nodes.begin();
1589  for (int gg = 0; nit != elem_nodes.end(); nit++, gg++) {
1590  CHKERR moab_ref.get_coords(&*nit, 1, &gauss_pts(gg, 0));
1591  }
1592  gauss_pts = trans(gauss_pts);
1593 
1594  MatrixDouble shape_fun;
1595  shape_fun.resize(nb_gauss_pts, 4);
1596  CHKERR ShapeMBTET(&*shape_fun.data().begin(), &gauss_pts(0, 0),
1597  &gauss_pts(1, 0), &gauss_pts(2, 0), nb_gauss_pts);
1598 
1599  double diff_shape_fun[12];
1600  CHKERR ShapeDiffMBTET(diff_shape_fun);
1601 
1602  // int edge_sense[6] = { 1,1,1, 1,1,1 };
1603  const int order = 5;
1604  // int edge_order[6] = { order,order,order, order,order,order };
1605  double def_val[] = {0, 0, 0, 0, 0, 0};
1606  int faces_order[] = {order, order, order, order};
1607  int faces_nodes[] = {0, 1, 3, 1, 2, 3, 0, 2, 3, 0, 1, 2};
1608 
1609  // cout << "NBEDGE_AINSWORTH_HCURL " << NBEDGE_AINSWORTH_HCURL(order) <<
1610  // endl; MatrixDouble
1611  // base_edge_functions(6,3*nb_gauss_pts*NBEDGE_AINSWORTH_HCURL(order));
1612  // double* edge_n[] = {
1613  // &base_edge_functions(0,0),
1614  // &base_edge_functions(1,0),
1615  // &base_edge_functions(2,0),
1616  // &base_edge_functions(3,0),
1617  // &base_edge_functions(4,0),
1618  // &base_edge_functions(5,0)
1619  // };
1620  //
1621  // MatrixDouble
1622  // diff_base_edge_functions(6,9*nb_gauss_pts*NBEDGE_AINSWORTH_HCURL(order));
1623  // double* diff_edge_n[] = {
1624  // &diff_base_edge_functions(0,0),
1625  // &diff_base_edge_functions(1,0),
1626  // &diff_base_edge_functions(2,0),
1627  // &diff_base_edge_functions(3,0),
1628  // &diff_base_edge_functions(4,0),
1629  // &diff_base_edge_functions(5,0)
1630  // };
1631  //
1632  // CHKERR Hcurl_Ainsworth_EdgeBaseFunctions_MBTET(
1633  // edge_sense,
1634  // edge_order,
1635  // &*shape_fun.data().begin(),
1636  // diff_shape_fun,
1637  // edge_n,
1638  // diff_edge_n,
1639  // nb_gauss_pts,
1640  // Legendre_polynomials
1641  // );
1642  //
1643  //
1644  // for(int ee = 0;ee!=6;ee++) {
1645  // for(int ll = 0;ll!=NBEDGE_AINSWORTH_HCURL(order);ll++) {
1646  // std::ostringstream ss;
1647  // ss << "curl_edge_" << ee << "_" << ll;
1648  // Tag th;
1649  // CHKERR moab_ref.tag_get_handle(
1650  // ss.str().c_str(),3,MB_TYPE_DOUBLE,th,MB_TAG_CREAT|MB_TAG_SPARSE,def_val
1651  // );
1652  // std::ostringstream ss_grad;
1653  // ss_grad << "grad_curl_edge_" << ee << "_" << ll;
1654  // Tag th_grad;
1655  // CHKERR moab_ref.tag_get_handle(
1656  // ss_grad.str().c_str(),9,MB_TYPE_DOUBLE,th_grad,MB_TAG_CREAT|MB_TAG_SPARSE,def_val
1657  // );
1658  //
1659  // int gg = 0;
1660  // for(Range::iterator nit =
1661  // elem_nodes.begin();nit!=elem_nodes.end();nit++,gg++) {
1662  // CHKERR moab_ref.tag_set_data(
1663  // th,&*nit,1,&(edge_n[ee][gg*3*NBEDGE_AINSWORTH_HCURL(order)+ll*3])
1664  // );
1665  // int sh = gg*9*NBEDGE_AINSWORTH_HCURL(order)+ll*9;
1666  // double grad[9] = {
1667  // diff_edge_n[ee][sh+0],diff_edge_n[ee][sh+3],diff_edge_n[ee][sh+6],
1668  // diff_edge_n[ee][sh+1],diff_edge_n[ee][sh+4],diff_edge_n[ee][sh+7],
1669  // diff_edge_n[ee][sh+2],diff_edge_n[ee][sh+5],diff_edge_n[ee][sh+8]
1670  // };
1671  // CHKERR moab_ref.tag_set_data(th_grad,&*nit,1,grad);
1672  //
1673  // }
1674  // }
1675  // }
1676 
1677  // cout << "NBFACETRI_AINSWORTH_EDGE_HCURL(order) " <<
1678  // NBFACETRI_AINSWORTH_EDGE_HCURL(order) << endl; MatrixDouble
1679  // base_face_edge_functions(
1680  // 4,3*3*NBFACETRI_AINSWORTH_EDGE_HCURL(order)*nb_gauss_pts
1681  // );
1682  // MatrixDouble diff_base_face_edge_functions(
1683  // 4,3*9*NBFACETRI_AINSWORTH_EDGE_HCURL(order)*nb_gauss_pts
1684  // );
1685  // double *phi_f_e[4][3];
1686  // double *diff_phi_f_e[4][3];
1687  // for(int ff = 0;ff!=4;ff++) {
1688  // for(int ee = 0;ee!=3;ee++) {
1689  // phi_f_e[ff][ee] =
1690  // &base_face_edge_functions(ff,ee*3*NBFACETRI_AINSWORTH_EDGE_HCURL(order)*nb_gauss_pts);
1691  // diff_phi_f_e[ff][ee] =
1692  // &diff_base_face_edge_functions(ff,ee*9*NBFACETRI_AINSWORTH_EDGE_HCURL(order)*nb_gauss_pts);
1693  // }
1694  // }
1695  //
1696  // CHKERR Hcurl_Ainsworth_EdgeBasedFaceFunctions_MBTET(
1697  // faces_nodes,
1698  // faces_order,
1699  // &*shape_fun.data().begin(),
1700  // diff_shape_fun,
1701  // phi_f_e,
1702  // diff_phi_f_e,
1703  // nb_gauss_pts,
1704  // Legendre_polynomials
1705  // );
1706  //
1707  // for(int ff = 0;ff!=4;ff++) {
1708  // for(int ee = 0;ee!=3;ee++) {
1709  // for(int ll = 0;ll!=NBFACETRI_AINSWORTH_EDGE_HCURL(order);ll++) {
1710  // std::ostringstream ss;
1711  // ss << "curl_face_edge_" << ff << "_" << ee << "_" << ll;
1712  // Tag th;
1713  // CHKERR moab_ref.tag_get_handle(
1714  // ss.str().c_str(),3,MB_TYPE_DOUBLE,th,MB_TAG_CREAT|MB_TAG_SPARSE,def_val
1715  // );
1716  //
1717  // std::ostringstream ss_grad;
1718  // ss_grad << "grad_curl_face_edge_" << ff << "_" << ee << "_" << ll;
1719  // Tag th_grad;
1720  // CHKERR moab_ref.tag_get_handle(
1721  // ss_grad.str().c_str(),9,MB_TYPE_DOUBLE,th_grad,MB_TAG_CREAT|MB_TAG_SPARSE,def_val
1722  // );
1723  //
1724  // int gg = 0;
1725  // for(Range::iterator nit =
1726  // elem_nodes.begin();nit!=elem_nodes.end();nit++,gg++) {
1727  //
1728  // int idx =
1729  // 3*NBFACETRI_AINSWORTH_EDGE_HCURL(order)*gg+ll*3;
1730  // if(idx >= base_face_edge_functions.size2()) {
1731  // cerr << ff << " " << ee << " " << ll << " " << gg << endl;
1732  // }
1733  //
1734  // CHKERR moab_ref.tag_set_data(th,&*nit,1,&(phi_f_e[ff][ee][idx]));
1735  //
1736  //
1737  // int sh = gg*9*NBFACETRI_AINSWORTH_EDGE_HCURL(order)+ll*9;
1738  // double grad[9] = {
1739  // diff_phi_f_e[ff][ee][sh+0],diff_phi_f_e[ff][ee][sh+3],diff_phi_f_e[ff][ee][sh+6],
1740  // diff_phi_f_e[ff][ee][sh+1],diff_phi_f_e[ff][ee][sh+4],diff_phi_f_e[ff][ee][sh+7],
1741  // diff_phi_f_e[ff][ee][sh+2],diff_phi_f_e[ff][ee][sh+5],diff_phi_f_e[ff][ee][sh+8]
1742  // };
1743  // CHKERR moab_ref.tag_set_data(th_grad,&*nit,1,grad);
1744  //
1745  //
1746  // }
1747  // }
1748  // }
1749  // }
1750  //
1751  cout << "NBFACETRI_AINSWORTH_FACE_HCURL "
1752  << NBFACETRI_AINSWORTH_FACE_HCURL(order) << endl;
1753  MatrixDouble base_face_bubble_functions(
1754  4, 3 * NBFACETRI_AINSWORTH_FACE_HCURL(order) * nb_gauss_pts);
1755  MatrixDouble diff_base_face_bubble_functions(
1756  4, 9 * NBFACETRI_AINSWORTH_FACE_HCURL(order) * nb_gauss_pts);
1757  double *phi_f[4];
1758  double *diff_phi_f[4];
1759  for (int ff = 0; ff != 4; ff++) {
1760  phi_f[ff] = &base_face_bubble_functions(ff, 0);
1761  diff_phi_f[ff] = &diff_base_face_bubble_functions(ff, 0);
1762  }
1763 
1765  faces_nodes, faces_order, &*shape_fun.data().begin(), diff_shape_fun,
1766  phi_f, diff_phi_f, nb_gauss_pts, Legendre_polynomials);
1767 
1768  for (int ff = 0; ff != 4; ff++) {
1769  for (int ll = 0; ll != NBFACETRI_AINSWORTH_FACE_HCURL(order); ll++) {
1770  std::ostringstream ss;
1771  ss << "curl_face_bubble_" << ff << "_" << ll;
1772  Tag th;
1773  CHKERR moab_ref.tag_get_handle(ss.str().c_str(), 3, MB_TYPE_DOUBLE, th,
1774  MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
1775  std::ostringstream grad_ss;
1776  grad_ss << "grad_curl_face_bubble_" << ff << "_" << ll;
1777  Tag th_grad;
1778  CHKERR moab_ref.tag_get_handle(grad_ss.str().c_str(), 9, MB_TYPE_DOUBLE,
1779  th_grad, MB_TAG_CREAT | MB_TAG_SPARSE,
1780  def_val);
1781 
1782  int gg = 0;
1783  for (Range::iterator nit = elem_nodes.begin(); nit != elem_nodes.end();
1784  nit++, gg++) {
1785  int idx = 3 * NBFACETRI_AINSWORTH_FACE_HCURL(order) * gg + ll * 3;
1786  CHKERR moab_ref.tag_set_data(th, &*nit, 1, &(phi_f[ff][idx]));
1787  int sh = gg * 9 * NBFACETRI_AINSWORTH_FACE_HCURL(order) + ll * 9;
1788  double grad[9] = {diff_phi_f[ff][sh + 0], diff_phi_f[ff][sh + 3],
1789  diff_phi_f[ff][sh + 6], diff_phi_f[ff][sh + 1],
1790  diff_phi_f[ff][sh + 4], diff_phi_f[ff][sh + 7],
1791  diff_phi_f[ff][sh + 2], diff_phi_f[ff][sh + 5],
1792  diff_phi_f[ff][sh + 8]};
1793  CHKERR moab_ref.tag_set_data(th_grad, &*nit, 1, grad);
1794  }
1795  }
1796  }
1797 
1798  // cout << "NBVOLUMETET_AINSWORTH_FACE_HCURL " <<
1799  // NBVOLUMETET_AINSWORTH_FACE_HCURL(order) << endl; VectorDouble
1800  // base_face_inetrior_functions(3*NBVOLUMETET_AINSWORTH_FACE_HCURL(order)*nb_gauss_pts);
1801  // VectorDouble
1802  // diff_base_face_inetrior_functions(9*NBVOLUMETET_AINSWORTH_FACE_HCURL(order)*nb_gauss_pts);
1803  // double *phi_v_f = &base_face_inetrior_functions[0];
1804  // double *diff_phi_v_f = &diff_base_face_inetrior_functions[0];
1805  // CHKERR Hcurl_Ainsworth_FaceInteriorFunctions_MBTET(
1806  // faces_nodes,
1807  // order,
1808  // &*shape_fun.data().begin(),
1809  // diff_shape_fun,
1810  // phi_v_f,
1811  // diff_phi_v_f,
1812  // nb_gauss_pts,
1813  // Legendre_polynomials
1814  // );
1815  // for(int ll = 0;ll!=NBVOLUMETET_AINSWORTH_FACE_HCURL(order);ll++) {
1816  //
1817  // std::ostringstream ss;
1818  // ss << "curl_face_interior_" << ll;
1819  // Tag th;
1820  // CHKERR moab_ref.tag_get_handle(
1821  // ss.str().c_str(),3,MB_TYPE_DOUBLE,th,MB_TAG_CREAT|MB_TAG_SPARSE,def_val
1822  // );
1823  //
1824  // std::ostringstream ss_grad;
1825  // ss_grad << "grad_curl_face_interior_" << ll;
1826  // Tag th_grad;
1827  // CHKERR moab_ref.tag_get_handle(
1828  // ss_grad.str().c_str(),9,MB_TYPE_DOUBLE,th_grad,MB_TAG_CREAT|MB_TAG_SPARSE,def_val
1829  // );
1830  //
1831  // int gg = 0;
1832  // for(Range::iterator nit =
1833  // elem_nodes.begin();nit!=elem_nodes.end();nit++,gg++) {
1834  // int idx = 3*NBVOLUMETET_AINSWORTH_FACE_HCURL(order)*gg+ll*3;
1835  // CHKERR moab_ref.tag_set_data(th,&*nit,1,&(phi_v_f[idx]));
1836  // int sh =
1837  // gg*9*NBVOLUMETET_AINSWORTH_FACE_HCURL(order)+ll*9; double grad[9] = {
1838  // diff_phi_v_f[sh+0],diff_phi_v_f[sh+3],diff_phi_v_f[sh+6],
1839  // diff_phi_v_f[sh+1],diff_phi_v_f[sh+4],diff_phi_v_f[sh+7],
1840  // diff_phi_v_f[sh+2],diff_phi_v_f[sh+5],diff_phi_v_f[sh+8]
1841  // };
1842  // CHKERR moab_ref.tag_set_data(th_grad,&*nit,1,grad);
1843  // }
1844  // }
1845 
1846  // cout << "NBVOLUMETET_AINSWORTH_TET_HCURL " <<
1847  // NBVOLUMETET_AINSWORTH_TET_HCURL(order) << endl; VectorDouble
1848  // base_interior_functions(3*NBVOLUMETET_AINSWORTH_TET_HCURL(order)*nb_gauss_pts);
1849  // VectorDouble
1850  // diff_base_interior_functions(9*NBVOLUMETET_AINSWORTH_TET_HCURL(order)*nb_gauss_pts);
1851  // double *phi_v = &base_interior_functions[0];
1852  // double *diff_phi_v = &diff_base_interior_functions[0];
1853  // CHKERR Hcurl_Ainsworth_VolumeInteriorFunctions_MBTET(
1854  // order,
1855  // &*shape_fun.data().begin(),
1856  // diff_shape_fun,
1857  // phi_v,
1858  // diff_phi_v,
1859  // nb_gauss_pts,
1860  // Legendre_polynomials
1861  // );
1862  // for(int ll = 0;ll!=NBVOLUMETET_AINSWORTH_TET_HCURL(order);ll++) {
1863  //
1864  // std::ostringstream ss;
1865  // ss << "curl_interior_" << ll;
1866  // Tag th;
1867  // CHKERR moab_ref.tag_get_handle(
1868  // ss.str().c_str(),3,MB_TYPE_DOUBLE,th,MB_TAG_CREAT|MB_TAG_SPARSE,def_val
1869  // );
1870  //
1871  // std::ostringstream ss_gard;
1872  // ss_gard << "grad_curl_interior_" << ll;
1873  // Tag th_grad;
1874  // CHKERR moab_ref.tag_get_handle(
1875  // ss_gard.str().c_str(),9,MB_TYPE_DOUBLE,th_grad,MB_TAG_CREAT|MB_TAG_SPARSE,def_val
1876  // );
1877  //
1878  // int gg = 0;
1879  // for(Range::iterator nit =
1880  // elem_nodes.begin();nit!=elem_nodes.end();nit++,gg++) {
1881  // int idx = 3*NBVOLUMETET_AINSWORTH_TET_HCURL(order)*gg+ll*3;
1882  // CHKERR moab_ref.tag_set_data(th,&*nit,1,&(phi_v[idx]));
1883  // int sh =
1884  // gg*9*NBVOLUMETET_AINSWORTH_TET_HCURL(order)+ll*9; double grad[9] = {
1885  // diff_phi_v[sh+0],diff_phi_v[sh+3],diff_phi_v[sh+6],
1886  // diff_phi_v[sh+1],diff_phi_v[sh+4],diff_phi_v[sh+7],
1887  // diff_phi_v[sh+2],diff_phi_v[sh+5],diff_phi_v[sh+8]
1888  // };
1889  // CHKERR moab_ref.tag_set_data(th_grad,&*nit,1,grad);
1890  // }
1891  // }
1892 
1893  // cout << "NBFACETRI_AINSWORTH_HCURL(order) " <<
1894  // NBFACETRI_AINSWORTH_HCURL(order) << endl; MatrixDouble base_face_functions(
1895  // 4,3*NBFACETRI_AINSWORTH_HCURL(order)*nb_gauss_pts
1896  // );
1897  // MatrixDouble diff_base_face_functions(
1898  // 4,9*NBFACETRI_AINSWORTH_HCURL(order)*nb_gauss_pts
1899  // );
1900  // for(int ff=0;ff!=4;ff++) {
1901  // phi_f[ff] = &base_face_functions(ff,0);
1902  // diff_phi_f[ff] = &diff_base_face_functions(ff,0);
1903  // }
1904  // CHKERR Hcurl_Ainsworth_FaceFunctions_MBTET(
1905  // faces_nodes,
1906  // faces_order,
1907  // &*shape_fun.data().begin(),
1908  // diff_shape_fun,
1909  // phi_f,
1910  // diff_phi_f,
1911  // nb_gauss_pts,
1912  // Legendre_polynomials
1913  // );
1914  // for(int ff = 0;ff!=4;ff++) {
1915  // for(int ll = 0;ll!=NBFACETRI_AINSWORTH_HCURL(order);ll++) {
1916  // std::ostringstream ss;
1917  // ss << "curl_face_" << ff << "_" << ll;
1918  // Tag th;
1919  // CHKERR moab_ref.tag_get_handle(
1920  // ss.str().c_str(),3,MB_TYPE_DOUBLE,th,MB_TAG_CREAT|MB_TAG_SPARSE,def_val
1921  // );
1922  //
1923  // int gg = 0;
1924  // for(Range::iterator nit =
1925  // elem_nodes.begin();nit!=elem_nodes.end();nit++,gg++) {
1926  // int idx = 3*NBFACETRI_AINSWORTH_HCURL (order)*gg+ll*3;
1927  // CHKERR moab_ref.tag_set_data(th,&*nit,1,&(phi_f[ff][idx]));
1928  //
1929  // }
1930  // }
1931  // }
1932  //
1933  // cout << "NBVOLUMETET_AINSWORTH_TET_HCURL(order) " <<
1934  // NBVOLUMETET_AINSWORTH_HCURL(order) << endl; VectorDouble
1935  // base_volume_functions(3*NBVOLUMETET_AINSWORTH_HCURL(order)*nb_gauss_pts);
1936  // VectorDouble
1937  // diff_base_volume_functions(9*NBVOLUMETET_AINSWORTH_HCURL(order)*nb_gauss_pts);
1938  // phi_v = &base_volume_functions[0]; diff_phi_v =
1939  // &diff_base_volume_functions[0]; CHKERR
1940  // MoFEM::Hcurl_Ainsworth_VolumeFunctions_MBTET(
1941  // order,
1942  // &*shape_fun.data().begin(),
1943  // diff_shape_fun,
1944  // phi_v,
1945  // diff_phi_v,
1946  // nb_gauss_pts,
1947  // Legendre_polynomials
1948  // );
1949  // for(int ll = 0;ll!=NBVOLUMETET_AINSWORTH_HCURL(order);ll++) {
1950  // std::ostringstream ss;
1951  // ss << "curl_volume_" << ll;
1952  // Tag th;
1953  // CHKERR moab_ref.tag_get_handle(
1954  // ss.str().c_str(),3,MB_TYPE_DOUBLE,th,MB_TAG_CREAT|MB_TAG_SPARSE,def_val
1955  // );
1956  //
1957  // int gg = 0;
1958  // for(Range::iterator nit =
1959  // elem_nodes.begin();nit!=elem_nodes.end();nit++,gg++) {
1960  // int idx = 3*NBVOLUMETET_AINSWORTH_HCURL(order)*gg+ll*3;
1961  // CHKERR moab_ref.tag_set_data(th,&*nit,1,&(phi_v[idx]));
1962  //
1963  // }
1964  // }
1965 
1966  EntityHandle meshset;
1967  CHKERR moab_ref.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER, meshset);
1968  CHKERR moab_ref.add_entities(meshset, tets);
1969  CHKERR moab_ref.write_file(file_name.c_str(), "VTK", "", &meshset, 1);
1970 
1972 }
1973 
1974 #endif // GENERATE_VTK_WITH_CURL_BASE
1975 
1976 #ifndef GENERATE_VTK_WITH_CURL_BASE
1977 
1979 
1981 
1985 
1988 
1989  template <int DIM, bool CALCULATE_DIRVATIVES>
1991  calculate(int p, int nb_integration_pts, int n0_idx, int n1_idx, double n[],
1992  FTensor::Tensor1<double, 3> t_grad_n[],
1995  *t_diff_phi_ptr) {
1996 
1998 
2000 
2001  FTensor::Tensor1<double, 3> &t_grad_n0 = t_grad_n[n0_idx];
2002  FTensor::Tensor1<double, 3> &t_grad_n1 = t_grad_n[n1_idx];
2003  tGradN0pN1(i) = t_grad_n0(i) + t_grad_n1(i);
2004 
2005  fI.resize(p + 1);
2006  diffFi.resize(3, p + 1);
2007  diffFi.clear();
2008 
2010 
2011  for (int gg = 0; gg != nb_integration_pts; ++gg) {
2012 
2013  const int shift_n = (DIM + 1) * gg;
2014  const double n0 = n[shift_n + n0_idx];
2015  const double n1 = n[shift_n + n1_idx];
2016 
2017  tPhi0(i) = n0 * t_grad_n1(i) - n1 * t_grad_n0(i);
2018  t_phi(i) = tPhi0(i);
2019 
2020  ++t_phi;
2021 
2022  if (CALCULATE_DIRVATIVES) {
2023 
2024  t_diff_phi0(i, j) =
2025  t_grad_n0(j) * t_grad_n1(i) - t_grad_n1(j) * t_grad_n0(i);
2026  (*t_diff_phi_ptr)(i, j) = t_diff_phi0(i, j);
2027  ++(*t_diff_phi_ptr);
2028  }
2029 
2030  if (p > 1) {
2031 
2032  if (CALCULATE_DIRVATIVES)
2033  CHKERR Jacobi_polynomials(p, 0, n1, n0 + n1, &t_grad_n1(0),
2034  &tGradN0pN1(0), &*fI.data().begin(),
2035  &*diffFi.data().begin(), DIM);
2036  else
2037  CHKERR Jacobi_polynomials(p, 0, n1, n0 + n1, nullptr, nullptr,
2038  &*fI.data().begin(), nullptr, DIM);
2039 
2041  &diffFi(0, 1), &diffFi(1, 1), &diffFi(2, 1));
2042 
2043  for (int oo = 1; oo <= p - 1; ++oo) {
2044 
2045  const double b = pow(n0 + n1, oo);
2046  t_phi(i) = b * fI[oo] * tPhi0(i);
2047 
2048  if (CALCULATE_DIRVATIVES) {
2049 
2050  tDiffb(i) =
2051  oo * pow(n0 + n1, oo - 1) * (t_grad_n0(i) + t_grad_n1(i));
2052  (*t_diff_phi_ptr)(i, j) = (b * fI[oo]) * t_diff_phi0(i, j) +
2053  (b * t_diff_fi(j)) * tPhi0(i) +
2054  tDiffb(j) * fI[oo] * tPhi0(i);
2055  ++t_diff_fi;
2056  ++(*t_diff_phi_ptr);
2057  }
2058 
2059  ++t_phi;
2060  }
2061  }
2062  }
2063 
2065  }
2066 };
2067 
2069  int *sense, int *p, double *n, double *diff_n, double *phi[],
2070  double *diff_phi[], int nb_integration_pts) {
2071 
2072  constexpr int e_nodes[6][2] = {{0, 1}, {1, 2}, {2, 0},
2073  {0, 3}, {1, 3}, {2, 3}};
2074 
2076 
2077  FTensor::Tensor1<double, 3> t_grad_n[4];
2078  for (int nn = 0; nn != 4; ++nn)
2079  t_grad_n[nn] = FTensor::Tensor1<double, 3>(
2080  diff_n[3 * nn + 0], diff_n[3 * nn + 1], diff_n[3 * nn + 2]);
2081 
2082  HcurlEdgeBase h_curl_base_on_edge;
2083 
2084  for (int ee = 0; ee != 6; ++ee) {
2085 
2086  auto t_phi = getFTensor1FromPtr<3>(phi[ee]);
2087  auto t_diff_phi = getFTensor2FromPtr<3, 3>(diff_phi[ee]);
2088 
2089  int n0_idx = e_nodes[ee][0];
2090  int n1_idx = e_nodes[ee][1];
2091  if (sense[ee] == -1) {
2092  int n_tmp = n0_idx;
2093  n0_idx = n1_idx;
2094  n1_idx = n_tmp;
2095  }
2096 
2097  CHKERR h_curl_base_on_edge.calculate<3, true>(p[ee], nb_integration_pts,
2098  n0_idx, n1_idx, n, t_grad_n,
2099  t_phi, &t_diff_phi);
2100  }
2101 
2103 }
2104 
2106  int *sense, int *p, double *n, double *diff_n, double *phi[],
2107  double *diff_phi[], int nb_integration_pts) {
2108 
2109  constexpr int e_nodes[3][2] = {{0, 1}, {1, 2}, {2, 0}};
2110 
2112 
2113  FTensor::Tensor1<double, 3> t_grad_n[3];
2114  for (int nn = 0; nn != 3; ++nn)
2115  t_grad_n[nn] =
2116  FTensor::Tensor1<double, 3>(diff_n[2 * nn + 0], diff_n[2 * nn + 1], 0.);
2117 
2118  HcurlEdgeBase h_curl_base_on_edge;
2119 
2120  for (int ee = 0; ee != 3; ++ee) {
2121 
2122  auto t_phi = getFTensor1FromPtr<3>(phi[ee]);
2123  auto t_diff_phi = getFTensor2FromPtr<3, 2>(diff_phi[ee]);
2124 
2125  int n0_idx = e_nodes[ee][0];
2126  int n1_idx = e_nodes[ee][1];
2127  if (sense[ee] == -1) {
2128  int n_tmp = n0_idx;
2129  n0_idx = n1_idx;
2130  n1_idx = n_tmp;
2131  }
2132 
2133  CHKERR h_curl_base_on_edge.calculate<2, true>(p[ee], nb_integration_pts,
2134  n0_idx, n1_idx, n, t_grad_n,
2135  t_phi, &t_diff_phi);
2136  }
2137 
2139 }
2140 
2142  int sense, int p, double *n, double *diff_n, double *phi, double *diff_phi,
2143  int nb_integration_pts) {
2145 
2146  FTensor::Tensor1<double, 3> t_grad_n[2];
2147  for (int nn = 0; nn != 2; ++nn)
2148  t_grad_n[nn] = FTensor::Tensor1<double, 3>(diff_n[nn], 0., 0.);
2149 
2150  HcurlEdgeBase h_curl_base_on_edge;
2151 
2153  &phi[HVEC0], &phi[HVEC1], &phi[HVEC2]);
2154 
2155  if (diff_phi != NULL)
2156  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
2157  "Not implemented direvatives for edge for Hcurl Demkowicz base");
2158 
2159  int n0_idx = 0;
2160  int n1_idx = 1;
2161  if (sense == -1) {
2162  int n_tmp = n0_idx;
2163  n0_idx = n1_idx;
2164  n1_idx = n_tmp;
2165  }
2166 
2167  CHKERR h_curl_base_on_edge.calculate<1, false>(
2168  p, nb_integration_pts, n0_idx, n1_idx, n, t_grad_n, t_phi, nullptr);
2169 
2171 }
2172 
2175  VectorDouble f0PhiII, diffF0PhiII;
2176  VectorDouble f1PhiII, diffF1PhiII;
2177  VectorDouble iFiF0, diffIFiF0;
2178  VectorDouble iFiF1, diffIFiF1;
2179 
2181 
2182  template <int DIM>
2184  int p, int nb_integration_pts, int n0f0_idx, int n1f0_idx, int n2f0_idx,
2185  double n[], FTensor::Tensor1<double, 3> t_grad_n[],
2188  &t_diff_phi) {
2189 
2191 
2193  f0PhiII.resize(3 * NBEDGE_DEMKOWICZ_HCURL(p) * nb_integration_pts, false);
2194  diffF0PhiII.resize(3 * DIM * NBEDGE_DEMKOWICZ_HCURL(p) * nb_integration_pts,
2195  false);
2196 
2197  // edge base for family I
2198  double *f0_phi_ii = &*f0PhiII.data().begin();
2199  double *diff_f0_phi_ii = &*diffF0PhiII.data().begin();
2200  auto t_f0_phi_ii = getFTensor1FromPtr<3>(f0_phi_ii);
2201  auto t_diff_f0_phi_ii = getFTensor2FromPtr<3, DIM>(diff_f0_phi_ii);
2202 
2203  CHKERR hCurlBaseOnEdge.calculate<DIM, true>(p - 1, nb_integration_pts,
2204  n0f0_idx, n1f0_idx, n, t_grad_n,
2205  t_f0_phi_ii, &t_diff_f0_phi_ii);
2206 
2207  FTensor::Tensor1<double, 3> &t_grad_n0f0 = t_grad_n[n0f0_idx];
2208  FTensor::Tensor1<double, 3> &t_grad_n1f0 = t_grad_n[n1f0_idx];
2209  FTensor::Tensor1<double, 3> &t_grad_n2f0 = t_grad_n[n2f0_idx];
2210  FTensor::Tensor1<double, 3> t_grad_n0f0_p_n1f0;
2211  t_grad_n0f0_p_n1f0(i) = t_grad_n0f0(i) + t_grad_n1f0(i) + t_grad_n2f0(i);
2212 
2213  iFiF0.resize(p + 1, false);
2214  diffIFiF0.resize(3 * p + 3, false);
2215  diffIFiF0.clear();
2216 
2217  double *ifif0 = &*iFiF0.data().begin();
2218  double *diff_ifif0 = &*diffIFiF0.data().begin();
2219 
2220  for (int gg = 0; gg != nb_integration_pts; ++gg) {
2221 
2222  const int shift_n = (DIM + 1) * gg;
2223  const double n0f0 = n[shift_n + n0f0_idx];
2224  const double n1f0 = n[shift_n + n1f0_idx];
2225  const double n2f0 = n[shift_n + n2f0_idx];
2226 
2227  int phi_shift = 3 * NBEDGE_DEMKOWICZ_HCURL(p - 1) * gg;
2228  int diff_phi_shift = (3 * DIM) * NBEDGE_DEMKOWICZ_HCURL(p - 1) * gg;
2229 
2230  for (int oo = 2; oo <= p; ++oo) {
2231 
2232  auto t_f0_phi_ii = getFTensor1FromPtr<3>(&f0_phi_ii[phi_shift]);
2233  auto t_diff_f0_phi_ii =
2234  getFTensor2FromPtr<3, DIM>(&diff_f0_phi_ii[diff_phi_shift]);
2235 
2236  for (int ii = 0; ii <= oo - 2; ii++) {
2237 
2238  int jj = oo - 2 - ii;
2239 
2240  // family I
2242  jj + 1, 2 * ii + 1, n2f0, n0f0 + n1f0 + n2f0, &t_grad_n2f0(0),
2243  &t_grad_n0f0_p_n1f0(0), ifif0, diff_ifif0, DIM);
2244  FTensor::Tensor1<double, 3> t_diff_ifif0(
2245  diff_ifif0[0 + jj], diff_ifif0[(jj + 1) + jj],
2246  diff_ifif0[2 * (jj + 1) + jj]);
2247  t_phi(i) = ifif0[jj] * t_f0_phi_ii(i);
2248  t_diff_phi(i, j) = ifif0[jj] * t_diff_f0_phi_ii(i, j) +
2249  t_diff_ifif0(j) * t_f0_phi_ii(i);
2250  ++t_phi;
2251  ++t_diff_phi;
2252  ++t_f0_phi_ii;
2253  ++t_diff_f0_phi_ii;
2254  }
2255  }
2256  }
2257 
2259  }
2260 
2261  template <int DIM>
2263  calculateTwoFamily(int p, int nb_integration_pts, int n0f0_idx, int n1f0_idx,
2264  int n2f0_idx, int n0f1_idx, int n1f1_idx, int n2f1_idx,
2265  double n[], FTensor::Tensor1<double, 3> t_grad_n[],
2268  DIM> &t_diff_phi) {
2269 
2271 
2273 
2274  f0PhiII.resize(3 * NBEDGE_DEMKOWICZ_HCURL(p) * nb_integration_pts, false);
2275  diffF0PhiII.resize(3 * DIM * NBEDGE_DEMKOWICZ_HCURL(p) * nb_integration_pts,
2276  false);
2277  f1PhiII.resize(3 * NBEDGE_DEMKOWICZ_HCURL(p) * nb_integration_pts, false);
2278  diffF1PhiII.resize(3 * DIM * NBEDGE_DEMKOWICZ_HCURL(p) * nb_integration_pts,
2279  false);
2280 
2281  // edge base for family I
2282  double *f0_phi_ii = &*f0PhiII.data().begin();
2283  double *diff_f0_phi_ii = &*diffF0PhiII.data().begin();
2284  auto t_f0_phi_ii = getFTensor1FromPtr<3>(f0_phi_ii);
2285  auto t_diff_f0_phi_ii = getFTensor2FromPtr<3, DIM>(diff_f0_phi_ii);
2286  CHKERR hCurlBaseOnEdge.calculate<DIM, true>(p - 1, nb_integration_pts,
2287  n0f0_idx, n1f0_idx, n, t_grad_n,
2288  t_f0_phi_ii, &t_diff_f0_phi_ii);
2289 
2290  // edge base for family II
2291  double *f1_phi_ii = &*f1PhiII.data().begin();
2292  double *diff_f1_phi_ii = &*diffF1PhiII.data().begin();
2293  auto t_f1_phi_ii = getFTensor1FromPtr<3>(f1_phi_ii);
2294  auto t_diff_f1_phi_ii = getFTensor2FromPtr<3, DIM>(diff_f1_phi_ii);
2295  CHKERR hCurlBaseOnEdge.calculate<DIM, true>(p - 1, nb_integration_pts,
2296  n0f1_idx, n1f1_idx, n, t_grad_n,
2297  t_f1_phi_ii, &t_diff_f1_phi_ii);
2298 
2299  FTensor::Tensor1<double, 3> &t_grad_n0f0 = t_grad_n[n0f0_idx];
2300  FTensor::Tensor1<double, 3> &t_grad_n1f0 = t_grad_n[n1f0_idx];
2301  FTensor::Tensor1<double, 3> &t_grad_n2f0 = t_grad_n[n2f0_idx];
2302  FTensor::Tensor1<double, 3> t_grad_n0f0_p_n1f0;
2303  t_grad_n0f0_p_n1f0(i) = t_grad_n0f0(i) + t_grad_n1f0(i);
2304 
2305  FTensor::Tensor1<double, 3> &t_grad_n0f1 = t_grad_n[n0f1_idx];
2306  FTensor::Tensor1<double, 3> &t_grad_n1f1 = t_grad_n[n1f1_idx];
2307  FTensor::Tensor1<double, 3> &t_grad_n2f1 = t_grad_n[n2f1_idx];
2308  FTensor::Tensor1<double, 3> t_grad_n0f1_p_n1f1;
2309  t_grad_n0f1_p_n1f1(i) = t_grad_n0f1(i) + t_grad_n1f1(i);
2310 
2311  iFiF0.resize(p + 1, false);
2312  diffIFiF0.resize(3 * p + 3, false);
2313  diffIFiF0.clear();
2314  double *ifif0 = &*iFiF0.data().begin();
2315  double *diff_ifif0 = &*diffIFiF0.data().begin();
2316  iFiF1.resize(p + 1, false);
2317  diffIFiF1.resize(3 * p + 3, false);
2318  diffIFiF1.clear();
2319  double *ifif1 = &*iFiF1.data().begin();
2320  double *diff_ifif1 = &*diffIFiF1.data().begin();
2321 
2322  for (int gg = 0; gg != nb_integration_pts; ++gg) {
2323 
2324  const int shift_n = (DIM + 1) * gg;
2325  const double n0f0 = n[shift_n + n0f0_idx];
2326  const double n1f0 = n[shift_n + n1f0_idx];
2327  const double n2f0 = n[shift_n + n2f0_idx];
2328  const double n0f1 = n[shift_n + n0f1_idx];
2329  const double n1f1 = n[shift_n + n1f1_idx];
2330  const double n2f1 = n[shift_n + n2f1_idx];
2331 
2332  int phi_shift = 3 * NBEDGE_DEMKOWICZ_HCURL(p - 1) * gg;
2333  int diff_phi_shift = 3 * DIM * NBEDGE_DEMKOWICZ_HCURL(p - 1) * gg;
2334 
2335  int kk = 0;
2336  for (int oo = 2; oo <= p; ++oo) {
2337 
2338  auto t_f0_phi_ii = getFTensor1FromPtr<3>(&f0_phi_ii[phi_shift]);
2339  auto t_diff_f0_phi_ii =
2340  getFTensor2FromPtr<3, DIM>(&diff_f0_phi_ii[diff_phi_shift]);
2341  auto t_f1_phi_ii = getFTensor1FromPtr<3>(&f1_phi_ii[phi_shift]);
2342  auto t_diff_f1_phi_ii =
2343  getFTensor2FromPtr<3, DIM>(&diff_f1_phi_ii[diff_phi_shift]);
2344 
2345  for (int ii = 0; ii <= oo - 2; ii++) {
2346 
2347  int jj = oo - 2 - ii;
2348 
2349  // family I
2351  jj + 1, 2 * ii + 1, n2f0, n0f0 + n1f0, &t_grad_n2f0(0),
2352  &t_grad_n0f0_p_n1f0(0), ifif0, diff_ifif0, 3);
2353  FTensor::Tensor1<double, 3> t_diff_ifif0(
2354  diff_ifif0[0 + jj], diff_ifif0[(jj + 1) + jj],
2355  diff_ifif0[2 * (jj + 1) + jj]);
2356 
2357  t_phi(i) = ifif0[jj] * t_f0_phi_ii(i);
2358  t_diff_phi(i, j) = ifif0[jj] * t_diff_f0_phi_ii(i, j) +
2359  t_diff_ifif0(j) * t_f0_phi_ii(i);
2360 
2361  ++t_phi;
2362  ++t_diff_phi;
2363  ++t_f0_phi_ii;
2364  ++t_diff_f0_phi_ii;
2365  ++kk;
2366 
2367  // family II
2369  jj + 1, 2 * ii + 1, n2f1, n0f1 + n1f1, &t_grad_n2f1(0),
2370  &t_grad_n0f1_p_n1f1(0), ifif1, diff_ifif1, 3);
2371  FTensor::Tensor1<double, 3> t_diff_ifif1(
2372  diff_ifif1[0 + jj], diff_ifif1[(jj + 1) + jj],
2373  diff_ifif1[2 * (jj + 1) + jj]);
2374  t_phi(i) = ifif1[jj] * t_f1_phi_ii(i);
2375  t_diff_phi(i, j) = ifif1[jj] * t_diff_f1_phi_ii(i, j) +
2376  t_diff_ifif1(j) * t_f1_phi_ii(i);
2377  ++t_phi;
2378  ++t_diff_phi;
2379  ++t_f1_phi_ii;
2380  ++t_diff_f1_phi_ii;
2381  ++kk;
2382  }
2383  }
2384  if (kk != NBFACETRI_DEMKOWICZ_HCURL(p))
2385  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2386  "Wrong number of base functions");
2387  }
2389  }
2390 };
2391 
2393  int *faces_nodes, int *p, double *n, double *diff_n, double *phi[],
2394  double *diff_phi[], int nb_integration_pts) {
2396 
2397  FTensor::Tensor1<double, 3> t_grad_n[4];
2398  for (int nn = 0; nn != 4; ++nn) {
2399  t_grad_n[nn] = FTensor::Tensor1<double, 3>(
2400  diff_n[3 * nn + 0], diff_n[3 * nn + 1], diff_n[3 * nn + 2]);
2401  };
2402 
2403  HcurlFaceBase h_curl_face_base;
2404 
2405  for (int ff = 0; ff != 4; ++ff) {
2406 
2407  if (p[ff] > 1) {
2408 
2409  auto t_phi = getFTensor1FromPtr<3>(phi[ff]);
2410  auto t_diff_phi = getFTensor2FromPtr<3, 3>(diff_phi[ff]);
2411 
2412  // f0, f1 - family I and family II
2413  const int n0f0_idx = faces_nodes[3 * ff + 0];
2414  const int n1f0_idx = faces_nodes[3 * ff + 1];
2415  const int n2f0_idx = faces_nodes[3 * ff + 2];
2416  // family II
2417  const int n0f1_idx = faces_nodes[3 * ff + 1];
2418  const int n1f1_idx = faces_nodes[3 * ff + 2];
2419  const int n2f1_idx = faces_nodes[3 * ff + 0];
2420 
2421  CHKERR h_curl_face_base.calculateTwoFamily<3>(
2422  p[ff], nb_integration_pts, n0f0_idx, n1f0_idx, n2f0_idx, n0f1_idx,
2423  n1f1_idx, n2f1_idx, n, t_grad_n, t_phi, t_diff_phi);
2424  }
2425  }
2426 
2428 }
2429 
2431  int *faces_nodes, int p, double *n, double *diff_n, double *phi,
2432  double *diff_phi, int nb_integration_pts) {
2434 
2435  FTensor::Tensor1<double, 3> t_grad_n[3];
2436  for (int nn = 0; nn != 3; ++nn) {
2437  t_grad_n[nn] =
2438  FTensor::Tensor1<double, 3>(diff_n[2 * nn + 0], diff_n[2 * nn + 1], 0.);
2439  };
2440 
2441  HcurlFaceBase h_curl_face_base;
2442 
2443  if (p > 1) {
2444 
2445  auto t_phi = getFTensor1FromPtr<3>(phi);
2446  auto t_diff_phi = getFTensor2FromPtr<3, 2>(diff_phi);
2447 
2448  // f0, f1 - family I and family II
2449  const int n0f0_idx = faces_nodes[0];
2450  const int n1f0_idx = faces_nodes[1];
2451  const int n2f0_idx = faces_nodes[2];
2452  // family II
2453  const int n0f1_idx = faces_nodes[1];
2454  const int n1f1_idx = faces_nodes[2];
2455  const int n2f1_idx = faces_nodes[0];
2456 
2457  CHKERR h_curl_face_base.calculateTwoFamily<2>(
2458  p, nb_integration_pts, n0f0_idx, n1f0_idx, n2f0_idx, n0f1_idx, n1f1_idx,
2459  n2f1_idx, n, t_grad_n, t_phi, t_diff_phi);
2460  }
2461 
2463 }
2464 
2466  int p, double *n, double *diff_n, double *phi, double *diff_phi,
2467  int nb_integration_pts) {
2468 
2469  constexpr int family[3][4] = {{0, 1, 2, 3}, {1, 2, 3, 0}, {2, 3, 0, 1}};
2472 
2474 
2475  if (p > 2) {
2476 
2477  auto t_phi = getFTensor1FromPtr<3>(phi);
2478  auto t_diff_phi = getFTensor2FromPtr<3, 3>(diff_phi);
2479 
2480  FTensor::Tensor1<double, 3> t_grad_n[4];
2481  for (int nn = 0; nn != 4; ++nn) {
2482  t_grad_n[nn] = FTensor::Tensor1<double, 3>(
2483  diff_n[3 * nn + 0], diff_n[3 * nn + 1], diff_n[3 * nn + 2]);
2484  };
2485 
2486  int nb_face_functions = (NBFACETRI_DEMKOWICZ_HCURL(p - 1)) / 2;
2487  MatrixDouble phi_ij(3, 3 * nb_face_functions * nb_integration_pts);
2488  MatrixDouble diff_phi_ij(3, 9 * nb_face_functions * nb_integration_pts);
2489  MatrixDouble fi_k(3, p + 1);
2490  MatrixDouble diff_fi_k(3, 3 * p + 3);
2491  HcurlFaceBase h_curl_face_base;
2492 
2493  // calate face base for each family
2494  for (int ff = 0; ff != 3; ++ff) {
2495  double *phi_ij_ptr = &phi_ij(ff, 0);
2496  double *diff_phi_ij_ptr = &diff_phi_ij(ff, 0);
2497 
2498  auto t_phi_ij = getFTensor1FromPtr<3>(phi_ij_ptr);
2499  auto t_diff_phi_ij = getFTensor2FromPtr<3, 3>(diff_phi_ij_ptr);
2500 
2501  const int n0_idx = family[ff][0];
2502  const int n1_idx = family[ff][1];
2503  const int n2_idx = family[ff][2];
2504 
2505  CHKERR h_curl_face_base.calculateOneFamily<3>(
2506  p - 1, nb_integration_pts, n0_idx, n1_idx, n2_idx, n, t_grad_n,
2507  t_phi_ij, t_diff_phi_ij);
2508  }
2509 
2510  FTensor::Tensor1<double, 3> &t_grad_n3f0 = t_grad_n[family[0][3]];
2511  FTensor::Tensor1<double, 3> &t_grad_n3f1 = t_grad_n[family[1][3]];
2512  FTensor::Tensor1<double, 3> &t_grad_n3f2 = t_grad_n[family[2][3]];
2513 
2514  FTensor::Tensor1<double, 3> t_sum_f0;
2515  t_sum_f0(i) = -t_grad_n3f0(i);
2516  FTensor::Tensor1<double, 3> t_sum_f1;
2517  t_sum_f1(i) = -t_grad_n3f1(i);
2518  FTensor::Tensor1<double, 3> t_sum_f2;
2519  t_sum_f2(i) = -t_grad_n3f2(i);
2520 
2521  for (int gg = 0; gg != nb_integration_pts; ++gg) {
2522 
2523  int shift_n = 4 * gg;
2524 
2525  double n3f0 = n[shift_n + family[0][3]];
2526  double n3f1 = n[shift_n + family[1][3]];
2527  double n3f2 = n[shift_n + family[2][3]];
2528 
2529  int kk = 0;
2530  for (int oo = 3; oo <= p; ++oo) {
2531 
2532  int phi_shift = 3 * nb_face_functions * gg;
2533  int diff_phi_shift = 9 * nb_face_functions * gg;
2534 
2535  auto t_phi_face_f0 = getFTensor1FromPtr<3>(&phi_ij(0, phi_shift));
2536  auto t_diff_phi_face_f0 =
2537  getFTensor2FromPtr<3, 3>(&diff_phi_ij(0, diff_phi_shift));
2538  auto t_phi_face_f1 = getFTensor1FromPtr<3>(&phi_ij(1, phi_shift));
2539  auto t_diff_phi_face_f1 =
2540  getFTensor2FromPtr<3, 3>(&diff_phi_ij(1, diff_phi_shift));
2541  auto t_phi_face_f2 = getFTensor1FromPtr<3>(&phi_ij(2, phi_shift));
2542  auto t_diff_phi_face_f2 =
2543  getFTensor2FromPtr<3, 3>(&diff_phi_ij(2, diff_phi_shift));
2544 
2545  int ij = 0;
2546  for (int oo_ij = 2; oo_ij != oo; ++oo_ij) {
2547  int k = oo - oo_ij;
2548 
2549  CHKERR IntegratedJacobi_polynomials(k, 2 * oo_ij, n3f0, 1 - n3f0,
2550  &t_grad_n3f0(0), &t_sum_f0(0),
2551  &fi_k(0, 0), &diff_fi_k(0, 0), 3);
2552  CHKERR IntegratedJacobi_polynomials(k, 2 * oo_ij, n3f1, 1 - n3f1,
2553  &t_grad_n3f1(0), &t_sum_f1(0),
2554  &fi_k(1, 0), &diff_fi_k(1, 0), 3);
2555  CHKERR IntegratedJacobi_polynomials(k, 2 * oo_ij, n3f2, 1 - n3f2,
2556  &t_grad_n3f2(0), &t_sum_f2(0),
2557  &fi_k(2, 0), &diff_fi_k(2, 0), 3);
2558 
2559  FTensor::Tensor1<double, 3> t_diff_fi_k_f0(
2560  diff_fi_k(0, 0 + k - 1), diff_fi_k(0, k + k - 1),
2561  diff_fi_k(0, 2 * k + k - 1));
2562  FTensor::Tensor1<double, 3> t_diff_fi_k_f1(
2563  diff_fi_k(1, 0 + k - 1), diff_fi_k(1, k + k - 1),
2564  diff_fi_k(1, 2 * k + k - 1));
2565  FTensor::Tensor1<double, 3> t_diff_fi_k_f2(
2566  diff_fi_k(2, 0 + k - 1), diff_fi_k(2, k + k - 1),
2567  diff_fi_k(2, 2 * k + k - 1));
2568 
2569  for (; ij != NBFACETRI_DEMKOWICZ_HCURL(oo_ij) / 2; ++ij) {
2570  t_phi(i) = fi_k(0, k - 1) * t_phi_face_f0(i);
2571  t_diff_phi(i, j) = t_diff_fi_k_f0(j) * t_phi_face_f0(i) +
2572  fi_k(0, k - 1) * t_diff_phi_face_f0(i, j);
2573  ++t_phi;
2574  ++t_diff_phi;
2575  ++t_phi_face_f0;
2576  ++t_diff_phi_face_f0;
2577  ++kk;
2578 
2579  t_phi(i) = fi_k(1, k - 1) * t_phi_face_f1(i);
2580  t_diff_phi(i, j) = t_diff_fi_k_f1(j) * t_phi_face_f1(i) +
2581  fi_k(1, k - 1) * t_diff_phi_face_f1(i, j);
2582  ++t_phi;
2583  ++t_diff_phi;
2584  ++t_phi_face_f1;
2585  ++t_diff_phi_face_f1;
2586  ++kk;
2587 
2588  t_phi(i) = fi_k(2, k - 1) * t_phi_face_f2(i);
2589  t_diff_phi(i, j) = t_diff_fi_k_f2(j) * t_phi_face_f2(i) +
2590  fi_k(2, k - 1) * t_diff_phi_face_f2(i, j);
2591  ++t_phi;
2592  ++t_diff_phi;
2593  ++t_phi_face_f2;
2594  ++t_diff_phi_face_f2;
2595  ++kk;
2596  }
2597  }
2598  }
2599  if (kk != NBVOLUMETET_DEMKOWICZ_HCURL(p))
2600  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2601  "Wrong number of base functions");
2602  }
2603  }
2604 
2606 }
2607 
2608 #endif
2609 
2610 #ifdef GENERATE_VTK_WITH_CURL_BASE
2611 
2612 MoFEMErrorCode VTK_Demkowicz_Hcurl_MBTET(const string file_name) {
2614 
2615  double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
2616 
2617  moab::Core core_ref;
2618  moab::Interface &moab_ref = core_ref;
2619 
2620  EntityHandle nodes[4];
2621  for (int nn = 0; nn < 4; nn++) {
2622  CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
2623  }
2624  EntityHandle tet;
2625  CHKERR moab_ref.create_element(MBTET, nodes, 4, tet);
2626 
2627  MoFEM::Core m_core_ref(moab_ref, PETSC_COMM_SELF, -2);
2628  MoFEM::Interface &m_field_ref = m_core_ref;
2629 
2630  CHKERR m_field_ref.getInterface<BitRefManager>()->setBitRefLevelByDim(
2631  0, 3, BitRefLevel().set(0));
2632 
2633  const int max_level = 3;
2634  for (int ll = 0; ll != max_level; ll++) {
2635  Range edges;
2636  CHKERR m_field_ref.getInterface<BitRefManager>()
2637  ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
2638  BitRefLevel().set(), MBEDGE, edges);
2639  Range tets;
2640  CHKERR m_field_ref.getInterface<BitRefManager>()
2641  ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
2642  BitRefLevel(ll).set(), MBTET, tets);
2643  // refine mesh
2644  MeshRefinement *m_ref;
2645  CHKERR m_field_ref.getInterface(m_ref);
2646  CHKERR m_ref->add_vertices_in_the_middle_of_edges(
2647  edges, BitRefLevel().set(ll + 1));
2648  CHKERR m_ref->refine_TET(tets, BitRefLevel().set(ll + 1));
2649  }
2650 
2651  Range tets;
2652  CHKERR m_field_ref.getInterface<BitRefManager>()
2653  ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(max_level),
2654  BitRefLevel().set(max_level), MBTET, tets);
2655 
2656  // Use 10 node tets to print base
2657  if (1) {
2658  EntityHandle meshset;
2659  CHKERR moab_ref.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER, meshset);
2660  CHKERR moab_ref.add_entities(meshset, tets);
2661  CHKERR moab_ref.convert_entities(meshset, true, false, false);
2662  CHKERR moab_ref.delete_entities(&meshset, 1);
2663  }
2664 
2665  Range elem_nodes;
2666  CHKERR moab_ref.get_connectivity(tets, elem_nodes, false);
2667 
2668  const int nb_gauss_pts = elem_nodes.size();
2669  MatrixDouble gauss_pts(nb_gauss_pts, 4);
2670  gauss_pts.clear();
2671  Range::iterator nit = elem_nodes.begin();
2672  for (int gg = 0; nit != elem_nodes.end(); nit++, gg++) {
2673  CHKERR moab_ref.get_coords(&*nit, 1, &gauss_pts(gg, 0));
2674  }
2675  gauss_pts = trans(gauss_pts);
2676 
2677  MatrixDouble shape_fun;
2678  shape_fun.resize(nb_gauss_pts, 4);
2679  CHKERR ShapeMBTET(&*shape_fun.data().begin(), &gauss_pts(0, 0),
2680  &gauss_pts(1, 0), &gauss_pts(2, 0), nb_gauss_pts);
2681 
2682  double diff_shape_fun[12];
2683  CHKERR ShapeDiffMBTET(diff_shape_fun);
2684 
2685  int edge_sense[6] = {1, 1, 1, 1, 1, 1};
2686  const int order = 4;
2687  int edge_order[6] = {order, order, order, order, order, order};
2688 
2689  MatrixDouble edge_phi(6, 3 * NBEDGE_DEMKOWICZ_HCURL(order) * nb_gauss_pts);
2690  MatrixDouble edge_diff_phi(6,
2691  9 * NBEDGE_DEMKOWICZ_HCURL(order) * nb_gauss_pts);
2692 
2693  edge_phi.clear();
2694  edge_diff_phi.clear();
2695 
2696  double *edge_phi_ptr[] = {&edge_phi(0, 0), &edge_phi(1, 0), &edge_phi(2, 0),
2697  &edge_phi(3, 0), &edge_phi(4, 0), &edge_phi(5, 0)};
2698  double *edge_diff_phi_ptr[] = {&edge_diff_phi(0, 0), &edge_diff_phi(1, 0),
2699  &edge_diff_phi(2, 0), &edge_diff_phi(3, 0),
2700  &edge_diff_phi(4, 0), &edge_diff_phi(5, 0)};
2701 
2703  edge_sense, edge_order, &*shape_fun.data().begin(), diff_shape_fun,
2704  edge_phi_ptr, edge_diff_phi_ptr, nb_gauss_pts);
2705 
2706  // cerr << edge_phi << endl;
2707 
2708  for (int ee = 0; ee != 6; ++ee) {
2709  for (int ll = 0; ll != NBEDGE_DEMKOWICZ_HCURL(order); ++ll) {
2710  double def_val[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
2711  std::string tag_name = "E" + boost::lexical_cast<std::string>(ee) + "_" +
2712  boost::lexical_cast<std::string>(ll);
2713  Tag th;
2714  CHKERR moab_ref.tag_get_handle(tag_name.c_str(), 3, MB_TYPE_DOUBLE, th,
2715  MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
2716 
2717  int gg = 0;
2718  for (Range::iterator nit = elem_nodes.begin(); nit != elem_nodes.end();
2719  nit++, gg++) {
2720  int idx = 3 * NBEDGE_DEMKOWICZ_HCURL(order) * gg + 3 * ll;
2721  CHKERR moab_ref.tag_set_data(th, &*nit, 1, &edge_phi(ee, idx));
2722  }
2723  }
2724  }
2725 
2726  int faces_order[] = {order, order, order, order};
2727  int faces_nodes[] = {0, 1, 3, 1, 2, 3, 0, 2, 3, 0, 1, 2};
2728 
2729  MatrixDouble face_phi(4, 3 * NBFACETRI_DEMKOWICZ_HCURL(order) * nb_gauss_pts);
2730  MatrixDouble face_diff_phi(4, 9 * NBFACETRI_DEMKOWICZ_HCURL(order) *
2731  nb_gauss_pts);
2732  face_phi.clear();
2733  face_diff_phi.clear();
2734 
2735  double *face_phi_ptr[] = {&face_phi(0, 0), &face_phi(1, 0), &face_phi(2, 0),
2736  &face_phi(3, 0)};
2737  double *face_diff_phi_ptr[] = {&face_diff_phi(0, 0), &face_diff_phi(1, 0),
2738  &face_diff_phi(2, 0), &face_diff_phi(3, 0)};
2739 
2741  faces_nodes, faces_order, &*shape_fun.data().begin(), diff_shape_fun,
2742  face_phi_ptr, face_diff_phi_ptr, nb_gauss_pts);
2743 
2744  for (int ff = 0; ff != 4; ++ff) {
2745  for (int ll = 0; ll != NBFACETRI_DEMKOWICZ_HCURL(order); ++ll) {
2746  double def_val[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
2747  std::string tag_name = "F" + boost::lexical_cast<std::string>(ff) + "_" +
2748  boost::lexical_cast<std::string>(ll);
2749  Tag th;
2750  CHKERR moab_ref.tag_get_handle(tag_name.c_str(), 3, MB_TYPE_DOUBLE, th,
2751  MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
2752 
2753  int gg = 0;
2754  for (Range::iterator nit = elem_nodes.begin(); nit != elem_nodes.end();
2755  nit++, gg++) {
2756  int idx = 3 * NBFACETRI_DEMKOWICZ_HCURL(order) * gg + 3 * ll;
2757  CHKERR moab_ref.tag_set_data(th, &*nit, 1, &face_phi(ff, idx));
2758  }
2759  }
2760  }
2761 
2762  MatrixDouble vol_phi(nb_gauss_pts, 3 * NBVOLUMETET_DEMKOWICZ_HCURL(order));
2763  MatrixDouble diff_vol_phi(nb_gauss_pts,
2764  9 * NBVOLUMETET_DEMKOWICZ_HCURL(order));
2765 
2767  order, &*shape_fun.data().begin(), diff_shape_fun, &vol_phi(0, 0),
2768  &diff_vol_phi(0, 0), nb_gauss_pts);
2769 
2770  for (int ll = 0; ll != NBVOLUMETET_DEMKOWICZ_HCURL(order); ++ll) {
2771  double def_val[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
2772  std::string tag_name = "V_" + boost::lexical_cast<std::string>(ll);
2773  Tag th;
2774  CHKERR moab_ref.tag_get_handle(tag_name.c_str(), 3, MB_TYPE_DOUBLE, th,
2775  MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
2776 
2777  int gg = 0;
2778  for (Range::iterator nit = elem_nodes.begin(); nit != elem_nodes.end();
2779  nit++, gg++) {
2780  int idx = 3 * ll;
2781  CHKERR moab_ref.tag_set_data(th, &*nit, 1, &vol_phi(gg, idx));
2782  }
2783  }
2784 
2785  EntityHandle meshset;
2786  CHKERR moab_ref.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER, meshset);
2787  CHKERR moab_ref.add_entities(meshset, tets);
2788  CHKERR moab_ref.write_file(file_name.c_str(), "VTK", "", &meshset, 1);
2789 
2791 }
2792 
2793 static char help[] = "...\n\n";
2794 
2795 int main(int argc, char *argv[]) {
2796 
2797  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
2798 
2799  try {
2800  // CHKERR
2801  // VTK_Ainsworth_Hcurl_MBTET("out_curl_vtk_ainsworth_base_on_tet.vtk");
2802  CHKERR VTK_Demkowicz_Hcurl_MBTET("out_curl_vtk_demkowicz_base_on_tet.vtk");
2803  }
2804  CATCH_ERRORS;
2805 
2807 
2808  return 0;
2809 }
2810 
2811 #endif // GENERATE_VTK_WITH_CURL_BASE
PetscErrorCode Legendre_polynomials(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
Calculate Legendre approximation basis.
MoFEMErrorCode Hcurl_Ainsworth_FaceFunctions_MBTET_ON_FACE(int *faces_nodes, int p, double *N, double *diffN, double *phi_f, double *diff_phi_f, int nb_integration_pts, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Face H-curl functions.
Definition: Hcurl.cpp:1249
#define NBVOLUMETET_AINSWORTH_FACE_HCURL(P)
FTensor::Tensor1< double, 3 > tGradN0pN1
Definition: Hcurl.cpp:1982
PetscErrorCode Jacobi_polynomials(int p, double alpha, double x, double t, double *diff_x, double *diff_t, double *L, double *diffL, const int dim)
Calculate Jacobi approximation basis.
PetscErrorCode ShapeMBTET(double *N, const double *G_X, const double *G_Y, const double *G_Z, int DIM)
calculate shape functions
Definition: fem_tools.c:318
VectorDouble f1PhiII
Definition: Hcurl.cpp:2176
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
MoFEMErrorCode Hcurl_Ainsworth_EdgeBaseFunctions_MBTET_ON_FACE(int *sense, int *p, double *N, double *diffN, double *edgeN[], double *diff_edgeN[], int nb_integration_pts, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Edge based H-curl base functions on face.
Definition: Hcurl.cpp:237
MoFEMErrorCode Hcurl_Demkowicz_VolumeBaseFunctions_MBTET(int p, double *n, double *diff_n, double *phi, double *diff_phi, int nb_integration_pts)
Volume base interior function.
Definition: Hcurl.cpp:2465
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
MoFEMErrorCode Hcurl_Demkowicz_FaceBaseFunctions_MBTRI(int *faces_nodes, int p, double *n, double *diff_n, double *phi, double *diff_phi, int nb_integration_pts)
Face base interior function.
Definition: Hcurl.cpp:2430
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
MoFEMErrorCode calculateOneFamily(int p, int nb_integration_pts, int n0f0_idx, int n1f0_idx, int n2f0_idx, double n[], FTensor::Tensor1< double, 3 > t_grad_n[], FTensor::Tensor1< FTensor::PackPtr< double *, 3 >, 3 > &t_phi, FTensor::Tensor2< FTensor::PackPtr< double *, DIM *3 >, 3, DIM > &t_diff_phi)
Definition: Hcurl.cpp:2183
MoFEMErrorCode Hcurl_Ainsworth_EdgeBasedFaceFunctions_MBTET_ON_FACE(int *faces_nodes, int p, double *N, double *diffN, double *phi_f_e[3], double *diff_phi_f_e[3], int nb_integration_pts, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Face edge base functions of Hcurl space.
Definition: Hcurl.cpp:458
Core (interface) class.
Definition: Core.hpp:50
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
HcurlEdgeBase hCurlBaseOnEdge
Definition: Hcurl.cpp:2174
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
MoFEMErrorCode Hcurl_Ainsworth_VolumeInteriorFunctions_MBTET(int p, double *N, double *diffN, double *phi_v, double *diff_phi_v, int nb_integration_pts, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Volume interior function.
Definition: Hcurl.cpp:909
T data[Tensor_Dim]
FTensor::Index< 'i', 3 > i
Definition: Hcurl.cpp:2180
MoFEMErrorCode calculate(int p, int nb_integration_pts, int n0_idx, int n1_idx, double n[], FTensor::Tensor1< double, 3 > t_grad_n[], FTensor::Tensor1< FTensor::PackPtr< double *, 3 >, 3 > &t_phi, FTensor::Tensor2< FTensor::PackPtr< double *, 3 *DIM >, 3, DIM > *t_diff_phi_ptr)
Definition: Hcurl.cpp:1991
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:51
MoFEMErrorCode calculateTwoFamily(int p, int nb_integration_pts, int n0f0_idx, int n1f0_idx, int n2f0_idx, int n0f1_idx, int n1f1_idx, int n2f1_idx, double n[], FTensor::Tensor1< double, 3 > t_grad_n[], FTensor::Tensor1< FTensor::PackPtr< double *, 3 >, 3 > &t_phi, FTensor::Tensor2< FTensor::PackPtr< double *, 3 *DIM >, 3, DIM > &t_diff_phi)
Definition: Hcurl.cpp:2263
FTensor::Tensor1< FTensor::PackPtr< double *, 3 >, 3 > getFTensor1FromPtr< 3 >(double *ptr)
Definition: Templates.hpp:336
VectorDouble f0PhiII
Definition: Hcurl.cpp:2175
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
FTensor::Tensor2< FTensor::PackPtr< double *, 9 >, 3, 3 > getFTensor2FromPtr< 3, 3 >(double *ptr)
Definition: Templates.hpp:369
MoFEMErrorCode Hcurl_Demkowicz_FaceBaseFunctions_MBTET(int *faces_nodes, int *p, double *n, double *diff_n, double *phi[], double *diff_phi[], int nb_integration_pts)
Face base interior function.
Definition: Hcurl.cpp:2392
PetscErrorCode ShapeDiffMBTET(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:331
MoFEMErrorCode Hcurl_Ainsworth_EdgeBasedFaceFunctions_MBTET(int *faces_nodes, int *p, double *N, double *diffN, double *phi_f_e[4][3], double *diff_phi_f_e[4][3], int nb_integration_pts, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Face edge base functions of Hcurl space on tetrahedral.
Definition: Hcurl.cpp:363
MoFEMErrorCode Hcurl_Demkowicz_EdgeBaseFunctions_MBTET(int *sense, int *p, double *n, double *diff_n, double *phi[], double *diff_phi[], int nb_integration_pts)
Edge based H-curl base functions on tetrahedral.
Definition: Hcurl.cpp:2068
#define NBEDGE_AINSWORTH_HCURL(P)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
Managing BitRefLevels.
FTensor::Tensor1< double, 3 > tPhi0
Definition: Hcurl.cpp:1983
MoFEMErrorCode Hcurl_Ainsworth_FaceFunctions_MBTET(int *face_nodes, int *p, double *N, double *diffN, double *phi_f[4], double *diff_phi_f[4], int nb_integration_pts, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Face H-curl functions.
Definition: Hcurl.cpp:1052
MoFEMErrorCode Hcurl_Demkowicz_EdgeBaseFunctions_MBEDGE(int sense, int p, double *n, double *diff_n, double *phi, double *diff_phi, int nb_integration_pts)
Edge based H-curl base functions on edge.
Definition: Hcurl.cpp:2141
FTensor::Index< 'i', 3 > i
Definition: Hcurl.cpp:1980
MoFEMErrorCode Hcurl_Ainsworth_BubbleFaceFunctions_MBTET_ON_FACE(int *faces_nodes, int p, double *N, double *diffN, double *phi_f, double *diff_phi_f, int nb_integration_pts, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Face edge base functions of Hcurl space on face.
Definition: Hcurl.cpp:663
Mesh refinement interface.
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define NBFACETRI_DEMKOWICZ_HCURL(P)
FTensor::Tensor2< FTensor::PackPtr< double *, 6 >, 3, 2 > getFTensor2FromPtr< 3, 2 >(double *ptr)
Definition: Templates.hpp:358
MatrixDouble diffFi
Definition: Hcurl.cpp:1987
MoFEMErrorCode Hcurl_Ainsworth_VolumeFunctions_MBTET(int p, double *N, double *diffN, double *phi_v, double *diff_phi_v, int nb_integration_pts, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
H-curl volume base functions.
Definition: Hcurl.cpp:1403
#define NBEDGE_DEMKOWICZ_HCURL(P)
static char help[]
MoFEMErrorCode Hcurl_Ainsworth_EdgeBaseFunctions_MBTET(int *sense, int *p, double *N, double *diffN, double *edgeN[], double *diff_edgeN[], int nb_integration_pts, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Edge based H-curl base functions on tetrahedral.
Definition: Hcurl.cpp:16
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
FTensor::Tensor1< double, 3 > tDiffb
Definition: Hcurl.cpp:1984
VectorDouble iFiF1
Definition: Hcurl.cpp:2178
#define NBVOLUMETET_AINSWORTH_HCURL(P)
#define NBFACETRI_AINSWORTH_HCURL(P)
#define NBFACETRI_AINSWORTH_EDGE_HCURL(P)
VectorDouble iFiF0
Definition: Hcurl.cpp:2177
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:433
MoFEMErrorCode Hcurl_Ainsworth_FaceInteriorFunctions_MBTET(int *faces_nodes, int p, double *N, double *diffN, double *phi_v, double *diff_phi_v, int nb_integration_pts, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Face base interior function.
Definition: Hcurl.cpp:775
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:72
PetscErrorCode IntegratedJacobi_polynomials(int p, double alpha, double x, double t, double *diff_x, double *diff_t, double *L, double *diffL, const int dim)
Calculate integrated Jacobi approximation basis.
VectorDouble fI
Definition: Hcurl.cpp:1986
const int N
Definition: speed_test.cpp:3
MoFEMErrorCode Hcurl_Ainsworth_BubbleFaceFunctions_MBTET(int *faces_nodes, int *p, double *N, double *diffN, double *phi_f[4], double *diff_phi_f[4], int nb_integration_pts, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Face edge base functions of Hcurl space on face on tetrahedral.
Definition: Hcurl.cpp:545
Implementation of H-curl base function.
#define NBVOLUMETET_AINSWORTH_TET_HCURL(P)
MoFEMErrorCode Hcurl_Demkowicz_EdgeBaseFunctions_MBTRI(int *sense, int *p, double *n, double *diff_n, double *phi[], double *diff_phi[], int nb_integration_pts)
Edge based H-curl base functions on teriangle.
Definition: Hcurl.cpp:2105
#define NBVOLUMETET_DEMKOWICZ_HCURL(P)
int main(int argc, char *argv[])
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:61
#define NBFACETRI_AINSWORTH_FACE_HCURL(P)
MoFEMErrorCode Hcurl_Ainsworth_EdgeBaseFunctions_MBTET_ON_EDGE(int sense, int p, double *N, double *diffN, double *edgeN, double *diff_edgeN, int nb_integration_pts, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Edge based H-curl base functions on edge.
Definition: Hcurl.cpp:175