v0.8.23
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 #ifndef GENERATE_VTK_WITH_CURL_BASE
12 
13 using namespace MoFEM;
14 
16  int *sense, int *p, double *N, double *diffN, double *edge_n[],
17  double *diff_edge_n[], int nb_integration_pts,
18  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
19  double *L, double *diffL,
20  const int dim)) {
22 
23  const int edges_nodes[6][2] = {{0, 1}, {1, 2}, {2, 0},
24  {0, 3}, {1, 3}, {2, 3}};
25  int P[6];
26  for (int ee = 0; ee != 6; ee++)
27  P[ee] = NBEDGE_AINSWORTH_HCURL(p[ee]);
28 
31 
32  FTensor::Tensor1<double *, 3> t_node_diff_ksi[4] = {
33  FTensor::Tensor1<double *, 3>(&diffN[0], &diffN[1], &diffN[2]),
34  FTensor::Tensor1<double *, 3>(&diffN[3], &diffN[4], &diffN[5]),
35  FTensor::Tensor1<double *, 3>(&diffN[6], &diffN[7], &diffN[8]),
36  FTensor::Tensor1<double *, 3>(&diffN[9], &diffN[10], &diffN[11])};
37  double edge_diff_ksi[6][3];
38  FTensor::Tensor1<double *, 3> t_edge_diff_ksi[6] = {
39  FTensor::Tensor1<double *, 3>(&edge_diff_ksi[0][0], &edge_diff_ksi[0][1],
40  &edge_diff_ksi[0][2]),
41  FTensor::Tensor1<double *, 3>(&edge_diff_ksi[1][0], &edge_diff_ksi[1][1],
42  &edge_diff_ksi[1][2]),
43  FTensor::Tensor1<double *, 3>(&edge_diff_ksi[2][0], &edge_diff_ksi[2][1],
44  &edge_diff_ksi[2][2]),
45  FTensor::Tensor1<double *, 3>(&edge_diff_ksi[3][0], &edge_diff_ksi[3][1],
46  &edge_diff_ksi[3][2]),
47  FTensor::Tensor1<double *, 3>(&edge_diff_ksi[4][0], &edge_diff_ksi[4][1],
48  &edge_diff_ksi[4][2]),
49  FTensor::Tensor1<double *, 3>(&edge_diff_ksi[5][0], &edge_diff_ksi[5][1],
50  &edge_diff_ksi[5][2])};
51  for (int ee = 0; ee != 6; ee++) {
52  t_edge_diff_ksi[ee](i) = (t_node_diff_ksi[edges_nodes[ee][1]](i) -
53  t_node_diff_ksi[edges_nodes[ee][0]](i)) *
54  sense[ee];
55  }
56 
57  FTensor::Tensor1<double *, 3> t_edge_n[6] = {
58  FTensor::Tensor1<double *, 3>(&edge_n[0][0], &edge_n[0][1], &edge_n[0][2],
59  3),
60  FTensor::Tensor1<double *, 3>(&edge_n[1][0], &edge_n[1][1], &edge_n[1][2],
61  3),
62  FTensor::Tensor1<double *, 3>(&edge_n[2][0], &edge_n[2][1], &edge_n[2][2],
63  3),
64  FTensor::Tensor1<double *, 3>(&edge_n[3][0], &edge_n[3][1], &edge_n[3][2],
65  3),
66  FTensor::Tensor1<double *, 3>(&edge_n[4][0], &edge_n[4][1], &edge_n[4][2],
67  3),
68  FTensor::Tensor1<double *, 3>(&edge_n[5][0], &edge_n[5][1], &edge_n[5][2],
69  3)};
70  FTensor::Tensor2<double *, 3, 3> t_diff_edge_n[6] = {
72  &diff_edge_n[0][0], &diff_edge_n[0][3], &diff_edge_n[0][6],
73  &diff_edge_n[0][1], &diff_edge_n[0][4], &diff_edge_n[0][7],
74  &diff_edge_n[0][2], &diff_edge_n[0][5], &diff_edge_n[0][8], 9),
76  &diff_edge_n[1][0], &diff_edge_n[1][3], &diff_edge_n[1][6],
77  &diff_edge_n[1][1], &diff_edge_n[1][4], &diff_edge_n[1][7],
78  &diff_edge_n[1][2], &diff_edge_n[1][5], &diff_edge_n[1][8], 9),
80  &diff_edge_n[2][0], &diff_edge_n[2][3], &diff_edge_n[2][6],
81  &diff_edge_n[2][1], &diff_edge_n[2][4], &diff_edge_n[2][7],
82  &diff_edge_n[2][2], &diff_edge_n[2][5], &diff_edge_n[2][8], 9),
84  &diff_edge_n[3][0], &diff_edge_n[3][3], &diff_edge_n[3][6],
85  &diff_edge_n[3][1], &diff_edge_n[3][4], &diff_edge_n[3][7],
86  &diff_edge_n[3][2], &diff_edge_n[3][5], &diff_edge_n[3][8], 9),
88  &diff_edge_n[4][0], &diff_edge_n[4][3], &diff_edge_n[4][6],
89  &diff_edge_n[4][1], &diff_edge_n[4][4], &diff_edge_n[4][7],
90  &diff_edge_n[4][2], &diff_edge_n[4][5], &diff_edge_n[4][8], 9),
92  &diff_edge_n[5][0], &diff_edge_n[5][3], &diff_edge_n[5][6],
93  &diff_edge_n[5][1], &diff_edge_n[5][4], &diff_edge_n[5][7],
94  &diff_edge_n[5][2], &diff_edge_n[5][5], &diff_edge_n[5][8], 9)};
95  FTensor::Tensor1<double, 3> t_psi_e_0, t_psi_e_1;
96  FTensor::Tensor2<double, 3, 3> t_diff_psi_e_0, t_diff_psi_e_1;
97 
98  for (int ii = 0; ii != nb_integration_pts; ii++) {
99 
100  const int node_shift = ii * 4;
101  for (int ee = 0; ee != 6; ee++) {
102 
103  if (P[ee] == 0)
104  continue;
105 
106  t_psi_e_0(i) = (N[node_shift + edges_nodes[ee][1]] *
107  t_node_diff_ksi[edges_nodes[ee][0]](i) -
108  N[node_shift + edges_nodes[ee][0]] *
109  t_node_diff_ksi[edges_nodes[ee][1]](i)) *
110  sense[ee];
111  t_diff_psi_e_0(i, j) = (t_node_diff_ksi[edges_nodes[ee][1]](j) *
112  t_node_diff_ksi[edges_nodes[ee][0]](i) -
113  t_node_diff_ksi[edges_nodes[ee][0]](j) *
114  t_node_diff_ksi[edges_nodes[ee][1]](i)) *
115  sense[ee];
116 
117  t_psi_e_1(i) = N[node_shift + edges_nodes[ee][1]] *
118  t_node_diff_ksi[edges_nodes[ee][0]](i) +
119  N[node_shift + edges_nodes[ee][0]] *
120  t_node_diff_ksi[edges_nodes[ee][1]](i);
121  t_diff_psi_e_1(i, j) = t_node_diff_ksi[edges_nodes[ee][1]](j) *
122  t_node_diff_ksi[edges_nodes[ee][0]](i) +
123  t_node_diff_ksi[edges_nodes[ee][0]](j) *
124  t_node_diff_ksi[edges_nodes[ee][1]](i);
125 
126  (t_edge_n[ee])(i) = t_psi_e_0(i);
127  ++(t_edge_n[ee]);
128  (t_edge_n[ee])(i) = t_psi_e_1(i);
129  ++(t_edge_n[ee]);
130 
131  (t_diff_edge_n[ee])(i, j) = t_diff_psi_e_0(i, j);
132  ++(t_diff_edge_n[ee]);
133  (t_diff_edge_n[ee])(i, j) = t_diff_psi_e_1(i, j);
134  ++(t_diff_edge_n[ee]);
135 
136  if (p[ee] > 1) {
137 
138  const double ksi_0i = (N[node_shift + edges_nodes[ee][1]] -
139  N[node_shift + edges_nodes[ee][0]]) *
140  sense[ee];
141  double psi_l[p[ee] + 1], diff_psi_l[3 * p[ee] + 3];
142  ierr = base_polynomials(p[ee], ksi_0i, &edge_diff_ksi[ee][0], psi_l,
143  diff_psi_l, 3);
144  CHKERRG(ierr);
145  FTensor::Tensor1<double *, 3> t_diff_psi_l(
146  &diff_psi_l[0], &diff_psi_l[p[ee] + 1], &diff_psi_l[2 * p[ee] + 2],
147  1);
148 
149  for (int ll = 2; ll != P[ee]; ll++) {
150 
151  const double a = (double)(2 * ll + 1) / (double)(ll + 1);
152  const double b = (double)(ll) / (double)(ll + 1);
153 
154  (t_edge_n[ee])(i) = a * psi_l[ll - 1] * t_psi_e_1(i) -
155  b * psi_l[ll - 2] * t_psi_e_0(i);
156  ++(t_edge_n[ee]);
157 
158  (t_diff_edge_n[ee])(i, j) =
159  -b * (t_diff_psi_l(j) * t_psi_e_0(i) +
160  psi_l[ll - 2] * t_diff_psi_e_0(i, j));
161  ++t_diff_psi_l;
162  (t_diff_edge_n[ee])(i, j) +=
163  a * (t_diff_psi_l(j) * t_psi_e_1(i) +
164  psi_l[ll - 1] * t_diff_psi_e_1(i, j));
165  ++(t_diff_edge_n[ee]);
166  }
167  }
168  }
169  }
170 
172 }
173 
175  int sense, int p, double *N, double *diffN, double *edge_n,
176  double *diff_edge_n, int nb_integration_pts,
177  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
178  double *L, double *diffL,
179  const int dim)) {
181 
182  if (NBEDGE_AINSWORTH_HCURL(p) == 0)
184  if (diff_edge_n != NULL) {
185  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
186  "Calculation of derivatives not implemented");
187  }
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 
198  FTensor::Tensor1<double *, 3> t_edge_n(&edge_n[0], &edge_n[1], &edge_n[2], 3);
199  FTensor::Tensor1<double, 3> t_psi_e_0, t_psi_e_1;
200 
201  for (int ii = 0; ii != nb_integration_pts; ii++) {
202 
203  const int node_shift = ii * 2;
204 
205  t_psi_e_0(i) = (N[node_shift + 1] * t_node_diff_ksi[0](i) -
206  N[node_shift + 0] * t_node_diff_ksi[1](i)) *
207  sense;
208  t_psi_e_1(i) = N[node_shift + 1] * t_node_diff_ksi[0](i) +
209  N[node_shift + 0] * t_node_diff_ksi[1](i);
210  t_edge_n(i) = t_psi_e_0(i);
211  ++t_edge_n;
212  t_edge_n(i) = t_psi_e_1(i);
213  ++t_edge_n;
214 
215  if (p > 1) {
216  const double ksi_0i = (N[node_shift + 1] - N[node_shift + 0]) * sense;
217  double psi_l[p + 1];
218  ierr = base_polynomials(p, ksi_0i, NULL, psi_l, NULL, 3);
219  CHKERRG(ierr);
220  for (int ll = 2; ll != NBEDGE_AINSWORTH_HCURL(p); ll++) {
221  const double a = (double)(2 * ll + 1) / (double)(ll + 1);
222  const double b = (double)(ll) / (double)(ll + 1);
223  t_edge_n(i) =
224  a * psi_l[ll - 1] * t_psi_e_1(i) - b * psi_l[ll - 2] * t_psi_e_0(i);
225  ++t_edge_n;
226  }
227  }
228  }
229 
231 }
232 
234  int *sense, int *p, double *N, double *diffN, double *edge_n[],
235  double *diff_edge_n[], int nb_integration_pts,
236  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
237  double *L, double *diffL,
238  const int dim)) {
239 
241 
242  // TODO This is not by atom tests properly
243 
244  const int edges_nodes[3][2] = {{0, 1}, {1, 2}, {2, 0}};
245  int P[3];
246  for (int ee = 0; ee < 3; ee++)
247  P[ee] = NBEDGE_AINSWORTH_HCURL(p[ee]);
248 
251 
252  FTensor::Tensor1<double, 3> t_node_diff_ksi[3] = {
253  FTensor::Tensor1<double, 3>(diffN[0], diffN[1], 0.),
254  FTensor::Tensor1<double, 3>(diffN[2], diffN[3], 0.),
255  FTensor::Tensor1<double, 3>(diffN[4], diffN[5], 0.),
256  };
257  FTensor::Tensor1<double, 2> t_2d_diff_ksi[3] = {
258  FTensor::Tensor1<double, 2>(diffN[0], diffN[1]),
259  FTensor::Tensor1<double, 2>(diffN[2], diffN[3]),
260  FTensor::Tensor1<double, 2>(diffN[4], diffN[5])};
261 
262  FTensor::Tensor1<double *, 3> t_edge_n[3] = {
263  FTensor::Tensor1<double *, 3>(&edge_n[0][0], &edge_n[0][1], &edge_n[0][2],
264  3),
265  FTensor::Tensor1<double *, 3>(&edge_n[1][0], &edge_n[1][1], &edge_n[1][2],
266  3),
267  FTensor::Tensor1<double *, 3>(&edge_n[2][0], &edge_n[2][1], &edge_n[2][2],
268  3)};
269  FTensor::Tensor2<FTensor::PackPtr<double *, 6>, 3, 2> t_diff_edge_n[3] = {
271  &diff_edge_n[0][HVEC0_0], &diff_edge_n[0][HVEC0_1],
272  &diff_edge_n[0][HVEC1_0], &diff_edge_n[0][HVEC1_1],
273  &diff_edge_n[0][HVEC2_0], &diff_edge_n[0][HVEC2_1]),
275  &diff_edge_n[1][HVEC0_0], &diff_edge_n[1][HVEC0_1],
276  &diff_edge_n[1][HVEC1_0], &diff_edge_n[1][HVEC1_1],
277  &diff_edge_n[1][HVEC2_0], &diff_edge_n[1][HVEC2_1]),
279  &diff_edge_n[2][HVEC0_0], &diff_edge_n[2][HVEC0_1],
280  &diff_edge_n[2][HVEC1_0], &diff_edge_n[2][HVEC1_1],
281  &diff_edge_n[2][HVEC2_0], &diff_edge_n[2][HVEC2_1])};
282 
283  FTensor::Tensor1<double, 3> t_psi_e_0, t_psi_e_1;
284  FTensor::Tensor2<double, 3, 2> t_diff_psi_e_0, t_diff_psi_e_1;
285 
286  for (int ee = 0; ee != 3; ee++) {
287 
288  if (P[ee] == 0)
289  continue;
290  const int node0 = edges_nodes[ee][0];
291  const int node1 = edges_nodes[ee][1];
292  const int sense_edge = sense[ee];
293 
294  t_diff_psi_e_0(i, j) =
295  (t_node_diff_ksi[node0](i) * t_2d_diff_ksi[node1](j) -
296  t_node_diff_ksi[node1](i) * t_2d_diff_ksi[node0](j)) *
297  sense_edge;
298  t_diff_psi_e_1(i, j) = t_node_diff_ksi[node0](i) * t_2d_diff_ksi[node1](j) +
299  t_node_diff_ksi[node1](i) * t_2d_diff_ksi[node0](j);
300 
301  for (int ii = 0; ii != nb_integration_pts; ii++) {
302 
303  const int node_shift = ii * 3;
304  const double n0 = N[node_shift + node0];
305  const double n1 = N[node_shift + node1];
306 
307  t_psi_e_0(i) =
308  (n1 * t_node_diff_ksi[node0](i) - n0 * t_node_diff_ksi[node1](i)) *
309  sense_edge;
310  t_psi_e_1(i) =
311  n1 * t_node_diff_ksi[node0](i) + n0 * t_node_diff_ksi[node1](i);
312 
313  (t_edge_n[ee])(i) = t_psi_e_0(i);
314  (t_diff_edge_n[ee])(i, j) = t_diff_psi_e_0(i, j);
315  ++(t_edge_n[ee]);
316  ++(t_diff_edge_n[ee]);
317  (t_edge_n[ee])(i) = t_psi_e_1(i);
318  (t_diff_edge_n[ee])(i, j) = t_diff_psi_e_1(i, j);
319  ++(t_edge_n[ee]);
320  ++(t_diff_edge_n[ee]);
321 
322  if (p[ee] > 1) {
323  const double ksi_0i = (n1 - n0) * sense_edge;
324  double diff_ksi_0i[] = {
325  ((t_2d_diff_ksi[node1])(0) - (t_2d_diff_ksi[node0])(0)) *
326  sense_edge,
327  ((t_2d_diff_ksi[node1])(1) - (t_2d_diff_ksi[node0])(1)) *
328  sense_edge};
329 
330  double psi_l[p[ee] + 1], diff_psi_l[2 * p[ee] + 2];
331  ierr =
332  base_polynomials(p[ee], ksi_0i, diff_ksi_0i, psi_l, diff_psi_l, 2);
333  CHKERRG(ierr);
334  FTensor::Tensor1<double *, 2> t_diff_psi_ll_m1(
335  &diff_psi_l[0 + 2 - 1], &diff_psi_l[p[ee] + 1 + 2 - 1], 1);
336  FTensor::Tensor1<double *, 2> t_diff_psi_ll_m2(
337  &diff_psi_l[0 + 2 - 2], &diff_psi_l[p[ee] + 1 + 2 - 2], 1);
338  for (int ll = 2; ll != P[ee]; ll++) {
339  const double a = (double)(2 * ll + 1) / (double)(ll + 1);
340  const double b = (double)(ll) / (double)(ll + 1);
341  (t_edge_n[ee])(i) = a * psi_l[ll - 1] * t_psi_e_1(i) -
342  b * psi_l[ll - 2] * t_psi_e_0(i);
343  (t_diff_edge_n[ee])(i, j) = a * t_psi_e_1(i) * t_diff_psi_ll_m1(j) +
344  a * psi_l[ll - 1] * t_diff_psi_e_1(i, j) -
345  b * t_psi_e_0(i) * t_diff_psi_ll_m2(j) -
346  b * psi_l[ll - 2] * t_diff_psi_e_0(i, j);
347  ++(t_edge_n[ee]);
348  ++(t_diff_edge_n[ee]);
349  ++t_diff_psi_ll_m1;
350  ++t_diff_psi_ll_m2;
351  }
352  }
353  }
354  }
355 
357 }
358 
360  int *faces_nodes, int *p, double *N, double *diffN, double *phi_f_e[4][3],
361  double *diff_phi_f_e[4][3], int nb_integration_pts,
362  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
363  double *L, double *diffL,
364  const int dim)) {
365 
367  const int edges[3][2] = {{0, 1}, {1, 2}, {2, 0}};
368 
371 
372  FTensor::Tensor1<double *, 3> t_node_diff_ksi[4] = {
373  FTensor::Tensor1<double *, 3>(&diffN[0], &diffN[1], &diffN[2]),
374  FTensor::Tensor1<double *, 3>(&diffN[3], &diffN[4], &diffN[5]),
375  FTensor::Tensor1<double *, 3>(&diffN[6], &diffN[7], &diffN[8]),
376  FTensor::Tensor1<double *, 3>(&diffN[9], &diffN[10], &diffN[11])};
377  FTensor::Tensor1<double, 3> t_edge_diff_ksi;
378  FTensor::Tensor1<double, 3> t_diff_beta_e;
379 
380  for (int ff = 0; ff != 4; ff++) {
381 
382  const int o_nodes[3] = {faces_nodes[3 * ff + 2], faces_nodes[3 * ff + 0],
383  faces_nodes[3 * ff + 1]};
384  FTensor::Tensor1<double *, 3> t_opposite_node_diff[3] = {
385  FTensor::Tensor1<double *, 3>(&diffN[3 * o_nodes[0] + 0],
386  &diffN[3 * o_nodes[0] + 1],
387  &diffN[3 * o_nodes[0] + 2]),
388  FTensor::Tensor1<double *, 3>(&diffN[3 * o_nodes[1] + 0],
389  &diffN[3 * o_nodes[1] + 1],
390  &diffN[3 * o_nodes[1] + 2]),
391  FTensor::Tensor1<double *, 3>(&diffN[3 * o_nodes[2] + 0],
392  &diffN[3 * o_nodes[2] + 1],
393  &diffN[3 * o_nodes[2] + 2])};
394  double psi_l[p[ff] + 1], diff_psi_l[3 * p[ff] + 3];
395 
396  const int nb_base_fun_on_face = NBFACETRI_AINSWORTH_EDGE_HCURL(p[ff]);
397  // cerr << nb_base_fun_on_face << " " << p[ff] << endl;
398  if (nb_base_fun_on_face == 0)
399  continue;
400 
401  for (int ee = 0; ee != 3; ee++) {
402 
403  FTensor::Tensor1<double *, 3> t_face_edge_base(
404  &phi_f_e[ff][ee][0], &phi_f_e[ff][ee][1], &phi_f_e[ff][ee][2], 3);
405  FTensor::Tensor2<double *, 3, 3> t_diff_face_edge_base(
406  &diff_phi_f_e[ff][ee][0], &diff_phi_f_e[ff][ee][3],
407  &diff_phi_f_e[ff][ee][6], &diff_phi_f_e[ff][ee][1],
408  &diff_phi_f_e[ff][ee][4], &diff_phi_f_e[ff][ee][7],
409  &diff_phi_f_e[ff][ee][2], &diff_phi_f_e[ff][ee][5],
410  &diff_phi_f_e[ff][ee][8], 9);
411  const int en[] = {faces_nodes[3 * ff + edges[ee][0]],
412  faces_nodes[3 * ff + edges[ee][1]]};
413  t_edge_diff_ksi(i) =
414  t_node_diff_ksi[en[1]](i) - t_node_diff_ksi[en[0]](i);
415 
416  for (int ii = 0; ii != nb_integration_pts; ii++) {
417 
418  const int node_shift = ii * 4;
419  const double n[] = {N[node_shift + faces_nodes[3 * ff + edges[ee][0]]],
420  N[node_shift + faces_nodes[3 * ff + edges[ee][1]]]};
421  const double ksi_0i = n[1] - n[0];
422  ierr = base_polynomials(p[ff], ksi_0i, &t_edge_diff_ksi(0), psi_l,
423  diff_psi_l, 3);
424  CHKERRG(ierr);
425  FTensor::Tensor1<double *, 3> t_diff_psi_l(
426  &diff_psi_l[0], &diff_psi_l[p[ff] + 1], &diff_psi_l[2 * p[ff] + 2],
427  1);
428 
429  const double beta_e = n[0] * n[1];
430  t_diff_beta_e(j) =
431  t_node_diff_ksi[en[0]](j) * n[1] + n[0] * t_node_diff_ksi[en[1]](j);
432 
433  for (int ll = 0; ll != nb_base_fun_on_face; ll++) {
434  // if(ll == nb_base_fun_on_face-1) cerr << psi_l[ll] << endl;
435 
436  t_face_edge_base(i) =
437  beta_e * psi_l[ll] * t_opposite_node_diff[ee](i);
438  ++t_face_edge_base;
439 
440  t_diff_face_edge_base(i, j) =
441  (beta_e * t_diff_psi_l(j) + t_diff_beta_e(j) * psi_l[ll]) *
442  t_opposite_node_diff[ee](i);
443 
444  ++t_diff_face_edge_base;
445  ++t_diff_psi_l;
446  }
447  }
448  }
449  }
450 
452 }
453 
455  int *faces_nodes, int p, double *N, double *diffN, double *phi_f_e[3],
456  double *diff_phi_f_e[3], int nb_integration_pts,
457  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
458  double *L, double *diffL,
459  const int dim)) {
460 
462 
463  const int nb_base_fun_on_face = NBFACETRI_AINSWORTH_EDGE_HCURL(p);
464  if (nb_base_fun_on_face == 0)
466 
467  const int edges[3][2] = {{0, 1}, {1, 2}, {2, 0}};
468 
471 
472  const int o_nodes[3] = {2, 0, 1};
473  FTensor::Tensor2<double, 3, 3> t_opposite_node_diff(
474  diffN[2 * o_nodes[0] + 0], diffN[2 * o_nodes[0] + 1], 0.,
475  diffN[2 * o_nodes[1] + 0], diffN[2 * o_nodes[1] + 1], 0.,
476  diffN[2 * o_nodes[2] + 0], diffN[2 * o_nodes[2] + 1], 0.);
477  double psi_l[p + 1];
478  double diff_psi_l[2 * p + 2];
479 
480  FTensor::Tensor1<double *, 3> t_face_edge_base[3] = {
481  FTensor::Tensor1<double *, 3>(&phi_f_e[0][HVEC0], &phi_f_e[0][HVEC1],
482  &phi_f_e[0][HVEC2], 3),
483  FTensor::Tensor1<double *, 3>(&phi_f_e[1][HVEC0], &phi_f_e[1][HVEC1],
484  &phi_f_e[1][HVEC2], 3),
485  FTensor::Tensor1<double *, 3>(&phi_f_e[2][HVEC0], &phi_f_e[2][HVEC1],
486  &phi_f_e[2][HVEC2], 3),
487  };
488  FTensor::Tensor2<FTensor::PackPtr<double *, 6>, 3, 2> t_diff_face_edge_base[3] = {
490  &diff_phi_f_e[0][HVEC0_0], &diff_phi_f_e[0][HVEC0_1],
491  &diff_phi_f_e[0][HVEC1_0], &diff_phi_f_e[0][HVEC1_1],
492  &diff_phi_f_e[0][HVEC2_0], &diff_phi_f_e[0][HVEC2_1]),
494  &diff_phi_f_e[1][HVEC0_0], &diff_phi_f_e[1][HVEC0_1],
495  &diff_phi_f_e[1][HVEC1_0], &diff_phi_f_e[1][HVEC1_1],
496  &diff_phi_f_e[1][HVEC2_0], &diff_phi_f_e[1][HVEC2_1]),
498  &diff_phi_f_e[2][HVEC0_0], &diff_phi_f_e[2][HVEC0_1],
499  &diff_phi_f_e[2][HVEC1_0], &diff_phi_f_e[2][HVEC1_1],
500  &diff_phi_f_e[2][HVEC2_0], &diff_phi_f_e[2][HVEC2_1])};
501 
502  for (int ee = 0; ee != 3; ee++) {
503 
504  const int node0 = faces_nodes[edges[ee][0]];
505  const int node1 = faces_nodes[edges[ee][1]];
506  double diff_ksi0i[] = {diffN[2 * node1 + 0] - diffN[2 * node0 + 0],
507  diffN[2 * node1 + 1] - diffN[2 * node0 + 1]};
508 
509  for (int ii = 0; ii != nb_integration_pts; ii++) {
510 
511  const int node_shift = ii * 3;
512  const double n0 = N[node_shift + node0];
513  const double n1 = N[node_shift + node1];
514  const double ksi_0i = n1 - n0;
515  ierr = base_polynomials(p, ksi_0i, diff_ksi0i, psi_l, diff_psi_l, 2);
516  CHKERRG(ierr);
517  const double beta_e = n0 * n1;
518  FTensor::Tensor1<double, 2> t_diff_beta_e(
519  diffN[2 * node0 + 0] * n1 + n0 * diffN[2 * node1 + 0],
520  diffN[2 * node0 + 1] * n1 + n0 * diffN[2 * node1 + 1]);
521  FTensor::Tensor1<double *, 2> t_diff_psi_l(&diff_psi_l[0],
522  &diff_psi_l[p + 1], 1);
523 
524  for (int ll = 0; ll != nb_base_fun_on_face; ll++) {
525  t_face_edge_base[ee](i) =
526  beta_e * psi_l[ll] * t_opposite_node_diff(ee, i);
527  t_diff_face_edge_base[ee](i, j) =
528  beta_e * t_opposite_node_diff(ee, i) * t_diff_psi_l(j) +
529  psi_l[ll] * t_opposite_node_diff(ee, i) * t_diff_beta_e(j);
530  ++t_face_edge_base[ee];
531  ++t_diff_face_edge_base[ee];
532  ++t_diff_psi_l;
533  }
534  }
535  }
536 
538 }
539 
541  int *faces_nodes, int *p, double *N, double *diffN, double *phi_f[4],
542  double *diff_phi_f[4], int nb_integration_pts,
543  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
544  double *L, double *diffL,
545  const int dim)) {
546 
548 
551 
552  // double coords[] = { 0,0,0, 1,0,0, 0,1,0, 0,0,1 };
553  // FTensor::Tensor1<double*,3> t_coords[4] = {
554  // FTensor::Tensor1<double*,3>(&coords[0],&coords[ 1],&coords[ 2]),
555  // FTensor::Tensor1<double*,3>(&coords[3],&coords[ 4],&coords[ 5]),
556  // FTensor::Tensor1<double*,3>(&coords[6],&coords[ 7],&coords[ 8]),
557  // FTensor::Tensor1<double*,3>(&coords[9],&coords[10],&coords[11])
558  // };
559  FTensor::Tensor1<double *, 3> t_node_diff_ksi[4] = {
560  FTensor::Tensor1<double *, 3>(&diffN[0], &diffN[1], &diffN[2]),
561  FTensor::Tensor1<double *, 3>(&diffN[3], &diffN[4], &diffN[5]),
562  FTensor::Tensor1<double *, 3>(&diffN[6], &diffN[7], &diffN[8]),
563  FTensor::Tensor1<double *, 3>(&diffN[9], &diffN[10], &diffN[11])};
564  FTensor::Tensor1<double, 3> t_diff_ksi0i, t_diff_ksi0j;
565  FTensor::Tensor1<double, 3> diff_beta_0ij;
566 
569 
570  for (int ff = 0; ff != 4; ff++) {
571 
572  if (NBFACETRI_AINSWORTH_FACE_HCURL(p[ff]) == 0)
573  continue;
574 
575  int n0 = faces_nodes[3 * ff + 0];
576  int n1 = faces_nodes[3 * ff + 1];
577  int n2 = faces_nodes[3 * ff + 2];
578 
579  // tou_0i(i) = t_coords[n1](i)-t_coords[n0](i);
580  // tou_0j(i) = t_coords[n2](i)-t_coords[n0](i);
581  tou_0i(i) = t_node_diff_ksi[n1](i) - t_node_diff_ksi[n0](i);
582  tou_0j(i) = t_node_diff_ksi[n2](i) - t_node_diff_ksi[n0](i);
583 
584  t_diff_ksi0i(i) = t_node_diff_ksi[n1](i) - t_node_diff_ksi[n0](i);
585  t_diff_ksi0j(i) = t_node_diff_ksi[n2](i) - t_node_diff_ksi[n0](i);
586 
587  double psi_l_0i[p[ff] + 1], diff_psi_l_0i[3 * p[ff] + 3];
588  double psi_l_0j[p[ff] + 1], diff_psi_l_0j[3 * p[ff] + 3];
589 
590  FTensor::Tensor1<double *, 3> t_phi_f(&phi_f[ff][0], &phi_f[ff][1],
591  &phi_f[ff][2], 3);
593  &diff_phi_f[ff][0], &diff_phi_f[ff][3], &diff_phi_f[ff][6],
594  &diff_phi_f[ff][1], &diff_phi_f[ff][4], &diff_phi_f[ff][7],
595  &diff_phi_f[ff][2], &diff_phi_f[ff][5], &diff_phi_f[ff][8], 9);
597 
598  for (int ii = 0; ii != nb_integration_pts; ii++) {
599 
600  const int node_shift = ii * 4;
601  const double beta_0ij =
602  N[node_shift + n0] * N[node_shift + n1] * N[node_shift + n2];
603  diff_beta_0ij(i) =
604  t_node_diff_ksi[n0](i) * N[node_shift + n1] * N[node_shift + n2] +
605  N[node_shift + n0] * t_node_diff_ksi[n1](i) * N[node_shift + n2] +
606  N[node_shift + n0] * N[node_shift + n1] * t_node_diff_ksi[n2](i);
607 
608  const double ksi_0i = N[node_shift + n1] - N[node_shift + n0];
609  ierr = base_polynomials(p[ff], ksi_0i, &t_diff_ksi0i(0), psi_l_0i,
610  diff_psi_l_0i, 3);
611  CHKERRG(ierr);
612  const double ksi_0j = N[node_shift + n2] - N[node_shift + n0];
613  ierr = base_polynomials(p[ff], ksi_0j, &t_diff_ksi0j(0), psi_l_0j,
614  diff_psi_l_0j, 3);
615  CHKERRG(ierr);
616 
617  int cc = 0;
618  for (int oo = 0; oo <= (p[ff] - 3); oo++) {
619  FTensor::Tensor1<double *, 3> t_diff_psi_l_0i(
620  &diff_psi_l_0i[0], &diff_psi_l_0i[p[ff] + 1],
621  &diff_psi_l_0i[2 * p[ff] + 2], 1);
622  for (int pp0 = 0; pp0 <= oo; pp0++) {
623  const int pp1 = oo - pp0;
624  if (pp1 >= 0) {
625  FTensor::Tensor1<double *, 3> t_diff_psi_l_0j(
626  &diff_psi_l_0j[pp1], &diff_psi_l_0j[p[ff] + 1 + pp1],
627  &diff_psi_l_0j[2 * p[ff] + 2 + pp1], 1);
628  // base functions
629  const double a = beta_0ij * psi_l_0i[pp0] * psi_l_0j[pp1];
630  t_phi_f(i) = a * tou_0i(i);
631  ++t_phi_f;
632  ++cc;
633  t_phi_f(i) = a * tou_0j(i);
634  ++t_phi_f;
635  ++cc;
636  // diff base functions
637  t_b(j) = diff_beta_0ij(j) * psi_l_0i[pp0] * psi_l_0j[pp1] +
638  beta_0ij * t_diff_psi_l_0i(j) * psi_l_0j[pp1] +
639  beta_0ij * psi_l_0i[pp0] * t_diff_psi_l_0j(j);
640  t_diff_phi_f(i, j) = t_b(j) * tou_0i(i);
641  ++t_diff_phi_f;
642  t_diff_phi_f(i, j) = t_b(j) * tou_0j(i);
643  ++t_diff_phi_f;
644  ++t_diff_psi_l_0i;
645  }
646  }
647  }
648  const int nb_base_fun_on_face = NBFACETRI_AINSWORTH_FACE_HCURL(p[ff]);
649  if (cc != nb_base_fun_on_face) {
650  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
651  "Wrong number of base functions %d != %d", cc,
652  nb_base_fun_on_face);
653  }
654  }
655  }
657 }
658 
660  int *faces_nodes, int p, double *N, double *diffN, double *phi_f,
661  double *diff_phi_f, int nb_integration_pts,
662  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
663  double *L, double *diffL,
664  const int dim)) {
665 
667 
668  double zero = 0;
669  FTensor::Tensor2<double *, 3, 3> t_node_diff_ksi(&diffN[0], &diffN[1], &zero,
670  &diffN[2], &diffN[3], &zero,
671  &diffN[4], &diffN[5], &zero);
672 
675 
676  if (NBFACETRI_AINSWORTH_FACE_HCURL(p) == 0)
678 
682 
683  const int node0 = faces_nodes[0];
684  const int node1 = faces_nodes[1];
685  const int node2 = faces_nodes[2];
686 
689 
690  tou_0i(i) = t_node_diff_ksi(N1, i) - t_node_diff_ksi(N0, i);
691  ;
692  tou_0j(i) = t_node_diff_ksi(N2, i) - t_node_diff_ksi(N0, i);
693  ;
694 
695  double psi_l_0i[p + 1], psi_l_0j[p + 1];
696  double diff_psi_l_0i[2 * p + 2], diff_psi_l_0j[2 * p + 2];
697  FTensor::Tensor1<double *, 3> t_phi_f(&phi_f[0], &phi_f[1], &phi_f[2], 3);
699  &diff_phi_f[HVEC0_0], &diff_phi_f[HVEC0_1], &diff_phi_f[HVEC1_0],
700  &diff_phi_f[HVEC1_1], &diff_phi_f[HVEC2_0], &diff_phi_f[HVEC2_1]);
701 
702  double diff_ksi_0i[] = {t_node_diff_ksi(node1, 0) - t_node_diff_ksi(node0, 0),
703  t_node_diff_ksi(node1, 1) -
704  t_node_diff_ksi(node0, 1)};
705  double diff_ksi_0j[] = {t_node_diff_ksi(node2, 0) - t_node_diff_ksi(node0, 0),
706  t_node_diff_ksi(node2, 1) -
707  t_node_diff_ksi(node0, 1)};
708 
709  for (int ii = 0; ii != nb_integration_pts; ii++) {
710 
711  const int node_shift = ii * 3;
712  const double n0 = N[node_shift + node0];
713  const double n1 = N[node_shift + node1];
714  const double n2 = N[node_shift + node2];
715 
716  const double beta_0ij = n0 * n1 * n2;
717  FTensor::Tensor1<double, 2> diff_beta_0ij(
718  t_node_diff_ksi(node0, 0) * n1 * n2 +
719  n0 * t_node_diff_ksi(node1, 0) * n2 +
720  n0 * n1 * t_node_diff_ksi(node2, 0),
721  t_node_diff_ksi(node0, 1) * n1 * n2 +
722  n0 * t_node_diff_ksi(node1, 1) * n2 +
723  n0 * n1 * t_node_diff_ksi(node2, 1));
724 
725  const double ksi_0i = N[node_shift + node1] - N[node_shift + node0];
726  ierr = base_polynomials(p, ksi_0i, diff_ksi_0i, psi_l_0i, diff_psi_l_0i, 2);
727  CHKERRG(ierr);
728  const double ksi_0j = N[node_shift + node2] - N[node_shift + node0];
729  ierr = base_polynomials(p, ksi_0j, diff_ksi_0j, psi_l_0j, diff_psi_l_0j, 2);
730  CHKERRG(ierr);
731 
732  int cc = 0;
734  for (int oo = 0; oo <= (p - 3); oo++) {
735  for (int pp0 = 0; pp0 <= oo; pp0++) {
736  const int pp1 = oo - pp0;
737  if (pp1 >= 0) {
738  FTensor::Tensor1<double, 2> t_diff_psi_l_0i(
739  diff_psi_l_0i[0 + pp0], diff_psi_l_0i[p + 1 + pp0]);
740  FTensor::Tensor1<double, 2> t_diff_psi_l_0j(
741  diff_psi_l_0j[0 + pp1], diff_psi_l_0j[p + 1 + pp1]);
742  const double a = beta_0ij * psi_l_0i[pp0] * psi_l_0j[pp1];
743  t_diff_a(j) = diff_beta_0ij(j) * psi_l_0i[pp0] * psi_l_0j[pp1] +
744  beta_0ij * psi_l_0i[pp0] * t_diff_psi_l_0j(j) +
745  beta_0ij * psi_l_0j[pp1] * t_diff_psi_l_0i(j);
746 
747  t_phi_f(i) = a * tou_0i(i);
748  t_diff_phi_f(i, j) = tou_0i(i) * t_diff_a(j);
749  ++t_phi_f;
750  ++t_diff_phi_f;
751  ++cc;
752  t_phi_f(i) = a * tou_0j(i);
753  t_diff_phi_f(i, j) = tou_0j(i) * t_diff_a(j);
754  ++t_phi_f;
755  ++t_diff_phi_f;
756  ++cc;
757  }
758  }
759  }
760 
761  const int nb_base_fun_on_face = NBFACETRI_AINSWORTH_FACE_HCURL(p);
762  if (cc != nb_base_fun_on_face) {
763  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
764  "Wrong number of base functions %d != %d", cc,
765  nb_base_fun_on_face);
766  }
767  }
768 
770 }
771 
773  int *faces_nodes, int p, double *N, double *diffN, double *phi_v,
774  double *diff_phi_v, int nb_integration_pts,
775  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
776  double *L, double *diffL,
777  const int dim)) {
778 
780 
783 
784  const int face_opposite_nodes[] = {2, 0, 1, 3};
785 
788 
789  FTensor::Tensor1<double *, 3> t_node_diff_ksi[4] = {
790  FTensor::Tensor1<double *, 3>(&diffN[0], &diffN[1], &diffN[2]),
791  FTensor::Tensor1<double *, 3>(&diffN[3], &diffN[4], &diffN[5]),
792  FTensor::Tensor1<double *, 3>(&diffN[6], &diffN[7], &diffN[8]),
793  FTensor::Tensor1<double *, 3>(&diffN[9], &diffN[10], &diffN[11])};
794  FTensor::Tensor1<double, 3> t_diff_ksi0i, t_diff_ksi0j;
795 
796  MatrixDouble m_psi_l_0i(4, p + 1);
797  MatrixDouble m_psi_l_0j(4, p + 1);
798  MatrixDouble m_diff_psi_l_0i(4, 3 * p + 3);
799  MatrixDouble m_diff_psi_l_0j(4, 3 * p + 3);
800 
801  double *psi_l_0i[] = {&m_psi_l_0i(0, 0), &m_psi_l_0i(1, 0), &m_psi_l_0i(2, 0),
802  &m_psi_l_0i(3, 0)};
803  double *psi_l_0j[] = {&m_psi_l_0j(0, 0), &m_psi_l_0j(1, 0), &m_psi_l_0j(2, 0),
804  &m_psi_l_0j(3, 0)};
805  double *diff_psi_l_0i[] = {&m_diff_psi_l_0i(0, 0), &m_diff_psi_l_0i(1, 0),
806  &m_diff_psi_l_0i(2, 0), &m_diff_psi_l_0i(3, 0)};
807  double *diff_psi_l_0j[] = {&m_diff_psi_l_0j(0, 0), &m_diff_psi_l_0j(1, 0),
808  &m_diff_psi_l_0j(2, 0), &m_diff_psi_l_0j(3, 0)};
809  double beta_f[4];
810 
811  FTensor::Tensor1<double, 3> t_diff_beta_f[4];
812 
813  FTensor::Tensor1<double *, 3> t_phi_v(&phi_v[0], &phi_v[1], &phi_v[2], 3);
815  &diff_phi_v[0], &diff_phi_v[3], &diff_phi_v[6], &diff_phi_v[1],
816  &diff_phi_v[4], &diff_phi_v[7], &diff_phi_v[2], &diff_phi_v[5],
817  &diff_phi_v[8], 9);
818 
819  for (int ii = 0; ii != nb_integration_pts; ii++) {
820 
821  for (int ff = 0; ff != 4; ff++) {
822 
823  t_diff_ksi0i(i) = t_node_diff_ksi[faces_nodes[3 * ff + 1]](i) -
824  t_node_diff_ksi[faces_nodes[3 * ff + 0]](i);
825  t_diff_ksi0j(i) = t_node_diff_ksi[faces_nodes[3 * ff + 2]](i) -
826  t_node_diff_ksi[faces_nodes[3 * ff + 0]](i);
827 
828  const int node_shift = ii * 4;
829 
830  beta_f[ff] = N[node_shift + faces_nodes[3 * ff + 0]] *
831  N[node_shift + faces_nodes[3 * ff + 1]] *
832  N[node_shift + faces_nodes[3 * ff + 2]];
833 
834  t_diff_beta_f[ff](j) = t_node_diff_ksi[faces_nodes[3 * ff + 0]](j) *
835  N[node_shift + faces_nodes[3 * ff + 1]] *
836  N[node_shift + faces_nodes[3 * ff + 2]] +
837  N[node_shift + faces_nodes[3 * ff + 0]] *
838  t_node_diff_ksi[faces_nodes[3 * ff + 1]](j) *
839  N[node_shift + faces_nodes[3 * ff + 2]] +
840  N[node_shift + faces_nodes[3 * ff + 0]] *
841  N[node_shift + faces_nodes[3 * ff + 1]] *
842  t_node_diff_ksi[faces_nodes[3 * ff + 2]](j);
843 
844  const double ksi_0i = N[node_shift + faces_nodes[3 * ff + 1]] -
845  N[node_shift + faces_nodes[3 * ff + 0]];
846  ierr = base_polynomials(p, ksi_0i, &t_diff_ksi0i(0), psi_l_0i[ff],
847  diff_psi_l_0i[ff], 3);
848  CHKERRG(ierr);
849 
850  const double ksi_0j = N[node_shift + faces_nodes[3 * ff + 2]] -
851  N[node_shift + faces_nodes[3 * ff + 0]];
852  ierr = base_polynomials(p, ksi_0j, &t_diff_ksi0j(0), psi_l_0j[ff],
853  diff_psi_l_0j[ff], 3);
854  CHKERRG(ierr);
855  }
856 
857  int cc = 0;
858  for (int oo = 0; oo <= (p - 3); oo++) {
859  FTensor::Tensor1<double *, 3> t_diff_psi_l_0i[] = {
860  FTensor::Tensor1<double *, 3>(&diff_psi_l_0i[0][0],
861  &diff_psi_l_0i[0][p + 1],
862  &diff_psi_l_0i[0][2 * p + 2], 1),
863  FTensor::Tensor1<double *, 3>(&diff_psi_l_0i[1][0],
864  &diff_psi_l_0i[1][p + 1],
865  &diff_psi_l_0i[1][2 * p + 2], 1),
866  FTensor::Tensor1<double *, 3>(&diff_psi_l_0i[2][0],
867  &diff_psi_l_0i[2][p + 1],
868  &diff_psi_l_0i[2][2 * p + 2], 1),
869  FTensor::Tensor1<double *, 3>(&diff_psi_l_0i[3][0],
870  &diff_psi_l_0i[3][p + 1],
871  &diff_psi_l_0i[3][2 * p + 2], 1),
872  };
873  for (int pp0 = 0; pp0 <= oo; pp0++) {
874  const int pp1 = oo - pp0;
875  if (pp1 >= 0) {
876  for (int ff = 0; ff != 4; ff++) {
877  FTensor::Tensor1<double *, 3> t_diff_psi_l_0j(
878  &m_diff_psi_l_0j(ff, pp1), &m_diff_psi_l_0j(ff, p + 1 + pp1),
879  &m_diff_psi_l_0j(ff, 2 * p + 2 + pp1), 1);
880  const double t = psi_l_0i[ff][pp0] * psi_l_0j[ff][pp1];
881  const double a = beta_f[ff] * t;
882  t_phi_v(i) = a * t_node_diff_ksi[face_opposite_nodes[ff]](i);
883  ++t_phi_v;
884  ++cc;
885  t_diff_phi_v(i, j) =
886  (t_diff_beta_f[ff](j) * t +
887  beta_f[ff] * t_diff_psi_l_0i[ff](j) * psi_l_0j[ff][pp1] +
888  beta_f[ff] * psi_l_0i[ff][pp0] * t_diff_psi_l_0j(j)) *
889  t_node_diff_ksi[face_opposite_nodes[ff]](i);
890  ++t_diff_phi_v;
891  ++t_diff_psi_l_0i[ff];
892  }
893  }
894  }
895  }
896 
897  const int nb_base_fun_on_face = NBVOLUMETET_AINSWORTH_FACE_HCURL(p);
898  if (cc != nb_base_fun_on_face) {
899  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
900  "Wrong number of base functions %d != %d", cc,
901  nb_base_fun_on_face);
902  }
903  }
904 
906 }
907 
909  int p, double *N, double *diffN, double *phi_v, double *diff_phi_v,
910  int nb_integration_pts,
911  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
912  double *L, double *diffL,
913  const int dim)) {
914 
916 
919 
925 
926  FTensor::Tensor1<double *, 3> t_node_diff_ksi[4] = {
927  FTensor::Tensor1<double *, 3>(&diffN[0], &diffN[1], &diffN[2]),
928  FTensor::Tensor1<double *, 3>(&diffN[3], &diffN[4], &diffN[5]),
929  FTensor::Tensor1<double *, 3>(&diffN[6], &diffN[7], &diffN[8]),
930  FTensor::Tensor1<double *, 3>(&diffN[9], &diffN[10], &diffN[11])};
931 
932  double diff_ksi0i[3], diff_ksi0j[3], diff_ksi0k[3];
933  FTensor::Tensor1<double *, 3> t_diff_ksi0i(diff_ksi0i, &diff_ksi0i[1],
934  &diff_ksi0i[2]);
935  FTensor::Tensor1<double *, 3> t_diff_ksi0j(diff_ksi0j, &diff_ksi0j[1],
936  &diff_ksi0j[2]);
937  FTensor::Tensor1<double *, 3> t_diff_ksi0k(diff_ksi0k, &diff_ksi0k[1],
938  &diff_ksi0k[2]);
939  t_diff_ksi0i(i) = t_node_diff_ksi[1](i) - t_node_diff_ksi[0](i);
940  t_diff_ksi0j(i) = t_node_diff_ksi[2](i) - t_node_diff_ksi[0](i);
941  t_diff_ksi0k(i) = t_node_diff_ksi[3](i) - t_node_diff_ksi[0](i);
942 
943  std::vector<double> v_psi_l_0i(p + 1), v_diff_psi_l_0i(3 * p + 3);
944  std::vector<double> v_psi_l_0j(p + 1), v_diff_psi_l_0j(3 * p + 3);
945  std::vector<double> v_psi_l_0k(p + 1), v_diff_psi_l_0k(3 * p + 3);
946  double *psi_l_0i = &*v_psi_l_0i.begin();
947  double *diff_psi_l_0i = &*v_diff_psi_l_0i.begin();
948  double *psi_l_0j = &*v_psi_l_0j.begin();
949  double *diff_psi_l_0j = &*v_diff_psi_l_0j.begin();
950  double *psi_l_0k = &*v_psi_l_0k.begin();
951  double *diff_psi_l_0k = &*v_diff_psi_l_0k.begin();
952 
953  FTensor::Tensor1<double *, 3> t_phi_v(&phi_v[0], &phi_v[1], &phi_v[2], 3);
955  &diff_phi_v[0], &diff_phi_v[3], &diff_phi_v[6], &diff_phi_v[1],
956  &diff_phi_v[4], &diff_phi_v[7], &diff_phi_v[2], &diff_phi_v[5],
957  &diff_phi_v[8], 9);
959 
960  for (int ii = 0; ii != nb_integration_pts; ii++) {
961 
962  const int node_shift = ii * 4;
963  const int n0 = node_shift + 0;
964  const int n1 = node_shift + 1;
965  const int n2 = node_shift + 2;
966  const int n3 = node_shift + 3;
967 
968  const double beta_v = N[n0] * N[n1] * N[n2] * N[n3];
969 
970  const double ksi_0i = N[n1] - N[n0];
971  ierr = base_polynomials(p, ksi_0i, diff_ksi0i, psi_l_0i, diff_psi_l_0i, 3);
972  CHKERRG(ierr);
973  const double ksi_0j = N[n2] - N[n0];
974  ierr = base_polynomials(p, ksi_0j, diff_ksi0j, psi_l_0j, diff_psi_l_0j, 3);
975  CHKERRG(ierr);
976  const double ksi_0k = N[n3] - N[n0];
977  ierr = base_polynomials(p, ksi_0k, diff_ksi0k, psi_l_0k, diff_psi_l_0k, 3);
978  CHKERRG(ierr);
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  CHKERRG(ierr);
1090 
1091  VectorDouble base_face_bubble_functions[4];
1092  VectorDouble diff_base_face_bubble_functions[4];
1093  double *phi_f_f[4];
1094  double *diff_phi_f_f[4];
1095  for (int ff = 0; ff != 4; ff++) {
1096  int nb_dofs = NBFACETRI_AINSWORTH_FACE_HCURL(p[ff]);
1097  if (nb_dofs == 0) {
1098  phi_f_f[ff] = NULL;
1099  diff_phi_f_f[ff] = NULL;
1100  } else {
1101  base_face_bubble_functions[ff].resize(3 * nb_dofs * nb_integration_pts,
1102  false);
1103  diff_base_face_bubble_functions[ff].resize(
1104  9 * nb_dofs * nb_integration_pts, false);
1105  phi_f_f[ff] = &*base_face_bubble_functions[ff].data().begin();
1106  diff_phi_f_f[ff] = &*diff_base_face_bubble_functions[ff].data().begin();
1107  }
1108  }
1110  face_nodes, p, N, diffN, phi_f_f, diff_phi_f_f, nb_integration_pts,
1111  base_polynomials);
1112  CHKERRG(ierr);
1113 
1116 
1117  for (int ff = 0; ff != 4; ff++) {
1118 
1119  if (NBFACETRI_AINSWORTH_EDGE_HCURL(p[ff]) == 0)
1120  continue;
1121 
1122  FTensor::Tensor1<double *, 3> t_face_edge_base[] = {
1123  FTensor::Tensor1<double *, 3>(&phi_f_e[ff][0][0], &phi_f_e[ff][0][1],
1124  &phi_f_e[ff][0][2], 3),
1125  FTensor::Tensor1<double *, 3>(&phi_f_e[ff][1][0], &phi_f_e[ff][1][1],
1126  &phi_f_e[ff][1][2], 3),
1127  FTensor::Tensor1<double *, 3>(&phi_f_e[ff][2][0], &phi_f_e[ff][2][1],
1128  &phi_f_e[ff][2][2], 3)};
1129  FTensor::Tensor2<double *, 3, 3> t_diff_face_edge_base[] = {
1131  &diff_phi_f_e[ff][0][0], &diff_phi_f_e[ff][0][3],
1132  &diff_phi_f_e[ff][0][6], &diff_phi_f_e[ff][0][1],
1133  &diff_phi_f_e[ff][0][4], &diff_phi_f_e[ff][0][7],
1134  &diff_phi_f_e[ff][0][2], &diff_phi_f_e[ff][0][5],
1135  &diff_phi_f_e[ff][0][8], 9),
1137  &diff_phi_f_e[ff][1][0], &diff_phi_f_e[ff][1][3],
1138  &diff_phi_f_e[ff][1][6], &diff_phi_f_e[ff][1][1],
1139  &diff_phi_f_e[ff][1][4], &diff_phi_f_e[ff][1][7],
1140  &diff_phi_f_e[ff][1][2], &diff_phi_f_e[ff][1][5],
1141  &diff_phi_f_e[ff][1][8], 9),
1143  &diff_phi_f_e[ff][2][0], &diff_phi_f_e[ff][2][3],
1144  &diff_phi_f_e[ff][2][6], &diff_phi_f_e[ff][2][1],
1145  &diff_phi_f_e[ff][2][4], &diff_phi_f_e[ff][2][7],
1146  &diff_phi_f_e[ff][2][2], &diff_phi_f_e[ff][2][5],
1147  &diff_phi_f_e[ff][2][8], 9)};
1148 
1149  FTensor::Tensor1<double *, 3> t_face_base(&phi_f[ff][0], &phi_f[ff][1],
1150  &phi_f[ff][2], 3);
1151  FTensor::Tensor2<double *, 3, 3> t_diff_face_base(
1152  &diff_phi_f[ff][0], &diff_phi_f[ff][3], &diff_phi_f[ff][6],
1153  &diff_phi_f[ff][1], &diff_phi_f[ff][4], &diff_phi_f[ff][7],
1154  &diff_phi_f[ff][2], &diff_phi_f[ff][5], &diff_phi_f[ff][8], 9);
1155 
1156  if (NBFACETRI_AINSWORTH_FACE_HCURL(p[ff]) > 0) {
1157  FTensor::Tensor1<double *, 3> t_face_face_base(
1158  &phi_f_f[ff][0], &phi_f_f[ff][1], &phi_f_f[ff][2], 3);
1159  FTensor::Tensor2<double *, 3, 3> t_diff_face_face_base(
1160  &diff_phi_f_f[ff][0], &diff_phi_f_f[ff][3], &diff_phi_f_f[ff][6],
1161  &diff_phi_f_f[ff][1], &diff_phi_f_f[ff][4], &diff_phi_f_f[ff][7],
1162  &diff_phi_f_f[ff][2], &diff_phi_f_f[ff][5], &diff_phi_f_f[ff][8],
1163  9);
1164  for (int ii = 0; ii != nb_integration_pts; ii++) {
1165  int cc = 0;
1166  for (int oo = 0; oo <= p[ff]; oo++) {
1167  // Face-edge base
1168  if (oo > 1) {
1169  for (int ee = 0; ee != 3; ee++) {
1170  for (int ll = NBFACETRI_AINSWORTH_EDGE_HCURL(oo - 1);
1171  ll != NBFACETRI_AINSWORTH_EDGE_HCURL(oo); ll++) {
1172  t_face_base(i) = t_face_edge_base[ee](i);
1173  ++cc;
1174  ++t_face_base;
1175  ++t_face_edge_base[ee];
1176  t_diff_face_base(i, j) = t_diff_face_edge_base[ee](i, j);
1177  ++t_diff_face_base;
1178  ++t_diff_face_edge_base[ee];
1179  // cerr << oo << " " << ll << " " << cc << " " <<
1180  // NBFACETRI_AINSWORTH_EDGE_HCURL(oo) << endl;
1181  }
1182  }
1183  }
1184  // Face-face base
1185  if (oo > 2) {
1186  for (int ll = NBFACETRI_AINSWORTH_FACE_HCURL(oo - 1);
1187  ll != NBFACETRI_AINSWORTH_FACE_HCURL(oo); ll++) {
1188  t_face_base(i) = t_face_face_base(i);
1189  ++cc;
1190  ++t_face_base;
1191  ++t_face_face_base;
1192  t_diff_face_base(i, j) = t_diff_face_face_base(i, j);
1193  ++t_diff_face_base;
1194  ++t_diff_face_face_base;
1195  }
1196  }
1197  }
1198  // check consistency
1199  const int nb_base_fun_on_face = NBFACETRI_AINSWORTH_HCURL(p[ff]);
1200  if (cc != nb_base_fun_on_face) {
1201  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1202  "Wrong number of base functions %d != %d", cc,
1203  nb_base_fun_on_face);
1204  }
1205  }
1206  } else {
1207  for (int ii = 0; ii != nb_integration_pts; ii++) {
1208  int cc = 0;
1209  for (int oo = 0; oo <= p[ff]; oo++) {
1210  // Face-edge base
1211  if (oo > 1) {
1212  for (int ee = 0; ee != 3; ee++) {
1213  for (int ll = NBFACETRI_AINSWORTH_EDGE_HCURL(oo - 1);
1214  ll != NBFACETRI_AINSWORTH_EDGE_HCURL(oo); ll++) {
1215  t_face_base(i) = t_face_edge_base[ee](i);
1216  ++cc;
1217  ++t_face_base;
1218  ++t_face_edge_base[ee];
1219  t_diff_face_base(i, j) = t_diff_face_edge_base[ee](i, j);
1220  ++t_diff_face_base;
1221  ++t_diff_face_edge_base[ee];
1222  // cerr << oo << " " << ll << " " << cc << " " <<
1223  // NBFACETRI_AINSWORTH_EDGE_HCURL(oo) << endl;
1224  }
1225  }
1226  }
1227  }
1228  // check consistency
1229  const int nb_base_fun_on_face = NBFACETRI_AINSWORTH_HCURL(p[ff]);
1230  if (cc != nb_base_fun_on_face) {
1231  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1232  "Wrong number of base functions %d != %d", cc,
1233  nb_base_fun_on_face);
1234  }
1235  }
1236  }
1237  }
1238 
1239  } catch (MoFEMException const &e) {
1240  SETERRQ(PETSC_COMM_SELF, e.errorCode, e.errorMessage);
1241  } catch (std::exception &ex) {
1242  std::ostringstream ss;
1243  ss << "thorw in method: " << ex.what() << " at line " << __LINE__
1244  << " in file " << __FILE__;
1245  SETERRQ(PETSC_COMM_SELF, MOFEM_STD_EXCEPTION_THROW, ss.str().c_str());
1246  }
1247 
1249 }
1250 
1252  int *faces_nodes, int p, double *N, double *diffN, double *phi_f,
1253  double *diff_phi_f, int nb_integration_pts,
1254  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
1255  double *L, double *diffL,
1256  const int dim)) {
1257 
1259 
1260  if (NBFACETRI_AINSWORTH_EDGE_HCURL(p) == 0)
1262 
1263  MatrixDouble base_face_edge_functions, diff_base_face_edge_functions;
1264  double *phi_f_e[3];
1265  double *diff_phi_f_e[3];
1266  base_face_edge_functions.resize(3, 3 * NBFACETRI_AINSWORTH_EDGE_HCURL(p) *
1267  nb_integration_pts);
1268  diff_base_face_edge_functions.resize(
1269  3, 2 * 3 * NBFACETRI_AINSWORTH_EDGE_HCURL(p) * nb_integration_pts);
1270  // base_face_edge_functions.clear();
1271  for (int ee = 0; ee != 3; ee++) {
1272  phi_f_e[ee] = &base_face_edge_functions(ee, 0);
1273  diff_phi_f_e[ee] = &diff_base_face_edge_functions(ee, 0);
1274  }
1276  faces_nodes, p, N, diffN, phi_f_e, diff_phi_f_e, nb_integration_pts,
1277  base_polynomials);
1278  CHKERRG(ierr);
1279 
1280  VectorDouble base_face_bubble_functions;
1281  VectorDouble diff_base_face_bubble_functions;
1282  double *phi_f_f, *diff_phi_f_f;
1283  base_face_bubble_functions.resize(3 * NBFACETRI_AINSWORTH_FACE_HCURL(p) *
1284  nb_integration_pts);
1285  diff_base_face_bubble_functions.resize(
1286  2 * 3 * NBFACETRI_AINSWORTH_FACE_HCURL(p) * nb_integration_pts);
1287  phi_f_f = &*base_face_bubble_functions.data().begin();
1288  diff_phi_f_f = &*diff_base_face_bubble_functions.data().begin();
1290  faces_nodes, p, N, diffN, phi_f_f, diff_phi_f_f, nb_integration_pts,
1291  base_polynomials);
1292  CHKERRG(ierr);
1293  // cerr << diff_base_face_bubble_functions << endl;
1294 
1297 
1298  FTensor::Tensor1<double *, 3> t_face_edge_base[] = {
1299  FTensor::Tensor1<double *, 3>(&phi_f_e[0][HVEC0], &phi_f_e[0][HVEC1],
1300  &phi_f_e[0][HVEC2], 3),
1301  FTensor::Tensor1<double *, 3>(&phi_f_e[1][HVEC0], &phi_f_e[1][HVEC1],
1302  &phi_f_e[1][HVEC2], 3),
1303  FTensor::Tensor1<double *, 3>(&phi_f_e[2][HVEC0], &phi_f_e[2][HVEC1],
1304  &phi_f_e[2][HVEC2], 3)};
1305  FTensor::Tensor2<FTensor::PackPtr<double *, 6>, 3, 2> t_diff_face_edge_base[] = {
1307  &diff_phi_f_e[0][HVEC0_0], &diff_phi_f_e[0][HVEC0_1],
1308  &diff_phi_f_e[0][HVEC1_0], &diff_phi_f_e[0][HVEC1_1],
1309  &diff_phi_f_e[0][HVEC2_0], &diff_phi_f_e[0][HVEC2_1]),
1311  &diff_phi_f_e[1][HVEC0_0], &diff_phi_f_e[1][HVEC0_1],
1312  &diff_phi_f_e[1][HVEC1_0], &diff_phi_f_e[1][HVEC1_1],
1313  &diff_phi_f_e[1][HVEC2_0], &diff_phi_f_e[1][HVEC2_1]),
1315  &diff_phi_f_e[2][HVEC0_0], &diff_phi_f_e[2][HVEC0_1],
1316  &diff_phi_f_e[2][HVEC1_0], &diff_phi_f_e[2][HVEC1_1],
1317  &diff_phi_f_e[2][HVEC2_0], &diff_phi_f_e[2][HVEC2_1])};
1318 
1319  FTensor::Tensor1<double *, 3> t_face_base(&phi_f[0], &phi_f[1], &phi_f[2], 3);
1320  FTensor::Tensor2<FTensor::PackPtr<double *, 6>, 3, 2> t_diff_face_base(
1321  &diff_phi_f[HVEC0_0], &diff_phi_f[HVEC0_1], &diff_phi_f[HVEC1_0],
1322  &diff_phi_f[HVEC1_1], &diff_phi_f[HVEC2_0], &diff_phi_f[HVEC2_1]);
1323 
1324  if (NBFACETRI_AINSWORTH_FACE_HCURL(p) > 0) {
1325  FTensor::Tensor1<double *, 3> t_face_face_base(
1326  &phi_f_f[HVEC0], &phi_f_f[HVEC1], &phi_f_f[HVEC2], 3);
1327  FTensor::Tensor2<FTensor::PackPtr<double *, 6>, 3, 2> t_diff_face_face_base(
1328  &diff_phi_f_f[HVEC0_0], &diff_phi_f_f[HVEC0_1],
1329  &diff_phi_f_f[HVEC1_0], &diff_phi_f_f[HVEC1_1],
1330  &diff_phi_f_f[HVEC2_0], &diff_phi_f_f[HVEC2_1]);
1331  for (int ii = 0; ii != nb_integration_pts; ii++) {
1332  int cc = 0;
1333  for (int oo = 0; oo <= p; oo++) {
1334  // Face-edge base
1335  if (oo > 1) {
1336  for (int ee = 0; ee != 3; ee++) {
1337  for (int ll = NBFACETRI_AINSWORTH_EDGE_HCURL(oo - 1);
1338  ll != NBFACETRI_AINSWORTH_EDGE_HCURL(oo); ll++) {
1339  t_face_base(i) = t_face_edge_base[ee](i);
1340  t_diff_face_base(i, j) = t_diff_face_edge_base[ee](i, j);
1341  ++cc;
1342  ++t_face_base;
1343  ++t_face_edge_base[ee];
1344  ++t_diff_face_base;
1345  ++t_diff_face_edge_base[ee];
1346  }
1347  }
1348  }
1349  // Face-face base
1350  if (oo > 2) {
1351  for (int ll = NBFACETRI_AINSWORTH_FACE_HCURL(oo - 1);
1352  ll != NBFACETRI_AINSWORTH_FACE_HCURL(oo); ll++) {
1353  t_face_base(i) = t_face_face_base(i);
1354  t_diff_face_base(i, j) = t_diff_face_face_base(i, j);
1355  ++cc;
1356  ++t_face_base;
1357  ++t_face_face_base;
1358  ++t_diff_face_base;
1359  ++t_diff_face_face_base;
1360  }
1361  }
1362  }
1363  // check consistency
1364  const int nb_base_fun_on_face = NBFACETRI_AINSWORTH_HCURL(p);
1365  if (cc != nb_base_fun_on_face) {
1366  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1367  "Wrong number of base functions %d != %d", cc,
1368  nb_base_fun_on_face);
1369  }
1370  }
1371  } else {
1372  for (int ii = 0; ii != nb_integration_pts; ii++) {
1373  int cc = 0;
1374  for (int oo = 0; oo <= p; oo++) {
1375  // Face-edge base
1376  if (oo > 1) {
1377  for (int ee = 0; ee != 3; ee++) {
1378  for (int ll = NBFACETRI_AINSWORTH_EDGE_HCURL(oo - 1);
1379  ll != NBFACETRI_AINSWORTH_EDGE_HCURL(oo); ll++) {
1380  t_face_base(i) = t_face_edge_base[ee](i);
1381  t_diff_face_base(i, j) = t_diff_face_edge_base[ee](i, j);
1382  ++cc;
1383  ++t_face_base;
1384  ++t_face_edge_base[ee];
1385  ++t_diff_face_base;
1386  ++t_diff_face_edge_base[ee];
1387  // cerr << oo << " " << ll << " " << cc << " " <<
1388  // NBFACETRI_AINSWORTH_EDGE_HCURL(oo) << endl;
1389  }
1390  }
1391  }
1392  }
1393  // check consistency
1394  const int nb_base_fun_on_face = NBFACETRI_AINSWORTH_HCURL(p);
1395  if (cc != nb_base_fun_on_face) {
1396  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1397  "Wrong number of base functions %d != %d", cc,
1398  nb_base_fun_on_face);
1399  }
1400  }
1401  }
1402 
1404 }
1405 
1407  int p, double *N, double *diffN, double *phi_v, double *diff_phi_v,
1408  int nb_integration_pts,
1409  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
1410  double *L, double *diffL,
1411  const int dim)) {
1412 
1414 
1415  VectorDouble base_face_inetrior_functions(
1416  3 * NBVOLUMETET_AINSWORTH_FACE_HCURL(p) * nb_integration_pts);
1417  VectorDouble diff_base_face_inetrior_functions(
1418  9 * NBVOLUMETET_AINSWORTH_FACE_HCURL(p) * nb_integration_pts);
1419  // base_face_inetrior_functions.clear();
1420  // diff_base_face_inetrior_functions.clear();
1421  double *phi_v_f = &*base_face_inetrior_functions.data().begin();
1422  double *diff_phi_v_f = &*diff_base_face_inetrior_functions.data().begin();
1423  int faces_nodes[] = {0, 1, 3, 1, 2, 3, 0, 2, 3, 0, 1, 2};
1425  faces_nodes, p, N, diffN, phi_v_f, diff_phi_v_f, nb_integration_pts,
1426  base_polynomials);
1427  CHKERRG(ierr);
1428 
1429  VectorDouble base_interior_functions(3 * NBVOLUMETET_AINSWORTH_TET_HCURL(p) *
1430  nb_integration_pts);
1431  VectorDouble diff_base_interior_functions(
1432  9 * NBVOLUMETET_AINSWORTH_TET_HCURL(p) * nb_integration_pts);
1433  // base_interior_functions.clear();
1434  // diff_base_interior_functions.clear();
1435  double *phi_v_v = &*base_interior_functions.data().begin();
1436  double *diff_phi_v_v = &*diff_base_interior_functions.data().begin();
1438  p, N, diffN, phi_v_v, diff_phi_v_v, nb_integration_pts, base_polynomials);
1439  CHKERRG(ierr);
1440 
1443 
1444  FTensor::Tensor1<double *, 3> t_face_interior(&phi_v_f[0], &phi_v_f[1],
1445  &phi_v_f[2], 3);
1446  FTensor::Tensor2<double *, 3, 3> t_diff_face_interior(
1447  &diff_phi_v_f[0], &diff_phi_v_f[3], &diff_phi_v_f[6], &diff_phi_v_f[1],
1448  &diff_phi_v_f[4], &diff_phi_v_f[7], &diff_phi_v_f[2], &diff_phi_v_f[5],
1449  &diff_phi_v_f[8], 9);
1450 
1451  FTensor::Tensor1<double *, 3> t_phi_v(&phi_v[0], &phi_v[1], &phi_v[2], 3);
1452  FTensor::Tensor2<double *, 3, 3> t_diff_phi_v(
1453  &diff_phi_v[0], &diff_phi_v[3], &diff_phi_v[6], &diff_phi_v[1],
1454  &diff_phi_v[4], &diff_phi_v[7], &diff_phi_v[2], &diff_phi_v[5],
1455  &diff_phi_v[8], 9);
1456 
1457  if (NBVOLUMETET_AINSWORTH_TET_HCURL(p) > 0) {
1458  FTensor::Tensor1<double *, 3> t_volume_interior(&phi_v_v[0], &phi_v_v[1],
1459  &phi_v_v[2], 3);
1460  FTensor::Tensor2<double *, 3, 3> t_diff_volume_interior(
1461  &diff_phi_v_v[0], &diff_phi_v_v[3], &diff_phi_v_v[6], &diff_phi_v_v[1],
1462  &diff_phi_v_v[4], &diff_phi_v_v[7], &diff_phi_v_v[2], &diff_phi_v_v[5],
1463  &diff_phi_v_v[8], 9);
1464  for (int ii = 0; ii != nb_integration_pts; ii++) {
1465  int cc = 0;
1466  for (int oo = 0; oo <= p; oo++) {
1467  for (int ll = NBVOLUMETET_AINSWORTH_FACE_HCURL(oo - 1);
1468  ll != NBVOLUMETET_AINSWORTH_FACE_HCURL(oo); ll++) {
1469  t_phi_v(i) = t_face_interior(i);
1470  ++t_phi_v;
1471  ++t_face_interior;
1472  ++cc;
1473  t_diff_phi_v(i, j) = t_diff_face_interior(i, j);
1474  ++t_diff_phi_v;
1475  ++t_diff_face_interior;
1476  }
1477  for (int ll = NBVOLUMETET_AINSWORTH_TET_HCURL(oo - 1);
1478  ll != NBVOLUMETET_AINSWORTH_TET_HCURL(oo); ll++) {
1479  t_phi_v(i) = t_volume_interior(i);
1480  ++t_phi_v;
1481  ++t_volume_interior;
1482  ++cc;
1483  t_diff_phi_v(i, j) = t_diff_volume_interior(i, j);
1484  ++t_diff_phi_v;
1485  ++t_diff_volume_interior;
1486  }
1487  }
1488  // check consistency
1489  const int nb_base_fun_on_face = NBVOLUMETET_AINSWORTH_HCURL(p);
1490  if (cc != nb_base_fun_on_face) {
1491  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1492  "Wrong number of base functions %d != %d", cc,
1493  nb_base_fun_on_face);
1494  }
1495  }
1496  } else {
1497  for (int ii = 0; ii != nb_integration_pts; ii++) {
1498  int cc = 0;
1499  for (int oo = 0; oo <= p; oo++) {
1500  for (int ll = NBVOLUMETET_AINSWORTH_FACE_HCURL(oo - 1);
1501  ll != NBVOLUMETET_AINSWORTH_FACE_HCURL(oo); ll++) {
1502  t_phi_v(i) = t_face_interior(i);
1503  ++t_phi_v;
1504  ++t_face_interior;
1505  ++cc;
1506  t_diff_phi_v(i, j) = t_diff_face_interior(i, j);
1507  ++t_diff_phi_v;
1508  ++t_diff_face_interior;
1509  }
1510  }
1511  // check consistency
1512  const int nb_base_fun_on_face = NBVOLUMETET_AINSWORTH_HCURL(p);
1513  if (cc != nb_base_fun_on_face) {
1514  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1515  "Wrong number of base functions %d != %d", cc,
1516  nb_base_fun_on_face);
1517  }
1518  }
1519  }
1520 
1522 }
1523 
1524 #endif // Not GENERATE_VTK_WITH_CURL_BASE
1525 
1526 #ifdef GENERATE_VTK_WITH_CURL_BASE
1527 
1528 #include <MoFEM.hpp>
1529 #include <Hcurl.hpp>
1530 using namespace MoFEM;
1531 using namespace boost::numeric;
1532 
1533 MoFEMErrorCode VTK_Ainsworth_Hcurl_MBTET(const string file_name) {
1535 
1536  ErrorCode rval;
1537 
1538  double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
1539 
1540  moab::Core core_ref;
1541  moab::Interface &moab_ref = core_ref;
1542 
1543  EntityHandle nodes[4];
1544  for (int nn = 0; nn < 4; nn++) {
1545  rval = moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
1546  CHKERRG(rval);
1547  }
1548  EntityHandle tet;
1549  rval = moab_ref.create_element(MBTET, nodes, 4, tet);
1550  CHKERRG(rval);
1551 
1552  MoFEM::Core m_core_ref(moab_ref, PETSC_COMM_SELF, -2);
1553  MoFEM::Interface &m_field_ref = m_core_ref;
1554 
1555  ierr = m_field_ref.getInterface<BitRefManager>()->setBitRefLevelByDim(
1556  0, 3, BitRefLevel().set(0));
1557  CHKERRG(ierr);
1558 
1559  const int max_level = 4;
1560  for (int ll = 0; ll != max_level; ll++) {
1561  Range edges;
1562  ierr =
1563  m_field_ref.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1564  BitRefLevel().set(ll), BitRefLevel().set(), MBEDGE, edges);
1565  CHKERRG(ierr);
1566  Range tets;
1567  ierr =
1568  m_field_ref.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1569  BitRefLevel().set(ll), BitRefLevel(ll).set(), MBTET, tets);
1570  CHKERRG(ierr);
1571  // refine mesh
1572  MeshRefinement *m_ref;
1573  ierr = m_field_ref.getInterface(m_ref);
1574  CHKERRG(ierr);
1576  BitRefLevel().set(ll + 1));
1577  CHKERRG(ierr);
1578  ierr = m_ref->refine_TET(tets, BitRefLevel().set(ll + 1));
1579  CHKERRG(ierr);
1580  }
1581 
1582  Range tets;
1583  ierr =
1584  m_field_ref.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1585  BitRefLevel().set(max_level), BitRefLevel().set(max_level), MBTET,
1586  tets);
1587  CHKERRG(ierr);
1588 
1589  // Use 10 node tets to print base
1590  if (1) {
1591 
1592  // Range edges;
1593  // rval = moab_ref.get_adjacencies(tets,1,true,edges); CHKERRG(rval);
1594  EntityHandle meshset;
1595  rval = moab_ref.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER, meshset);
1596  CHKERRG(rval);
1597  rval = moab_ref.add_entities(meshset, tets);
1598  CHKERRG(rval);
1599  rval = moab_ref.convert_entities(meshset, true, false, false);
1600  CHKERRG(rval);
1601  rval = moab_ref.delete_entities(&meshset, 1);
1602  CHKERRG(rval);
1603  }
1604 
1605  Range elem_nodes;
1606  rval = moab_ref.get_connectivity(tets, elem_nodes, false);
1607  CHKERRG(rval);
1608 
1609  const int nb_gauss_pts = elem_nodes.size();
1610  MatrixDouble gauss_pts(nb_gauss_pts, 4);
1611  gauss_pts.clear();
1612  Range::iterator nit = elem_nodes.begin();
1613  for (int gg = 0; nit != elem_nodes.end(); nit++, gg++) {
1614  rval = moab_ref.get_coords(&*nit, 1, &gauss_pts(gg, 0));
1615  CHKERRG(rval);
1616  }
1617  gauss_pts = trans(gauss_pts);
1618 
1619  MatrixDouble shape_fun;
1620  shape_fun.resize(nb_gauss_pts, 4);
1621  ierr = ShapeMBTET(&*shape_fun.data().begin(), &gauss_pts(0, 0),
1622  &gauss_pts(1, 0), &gauss_pts(2, 0), nb_gauss_pts);
1623  CHKERRG(ierr);
1624 
1625  double diff_shape_fun[12];
1626  ierr = ShapeDiffMBTET(diff_shape_fun);
1627  CHKERRG(ierr);
1628 
1629  // int edge_sense[6] = { 1,1,1, 1,1,1 };
1630  const int order = 5;
1631  // int edge_order[6] = { order,order,order, order,order,order };
1632  double def_val[] = {0, 0, 0, 0, 0, 0};
1633  int faces_order[] = {order, order, order, order};
1634  int faces_nodes[] = {0, 1, 3, 1, 2, 3, 0, 2, 3, 0, 1, 2};
1635 
1636  // cout << "NBEDGE_AINSWORTH_HCURL " << NBEDGE_AINSWORTH_HCURL(order) <<
1637  // endl; MatrixDouble
1638  // base_edge_functions(6,3*nb_gauss_pts*NBEDGE_AINSWORTH_HCURL(order));
1639  // double* edge_n[] = {
1640  // &base_edge_functions(0,0),
1641  // &base_edge_functions(1,0),
1642  // &base_edge_functions(2,0),
1643  // &base_edge_functions(3,0),
1644  // &base_edge_functions(4,0),
1645  // &base_edge_functions(5,0)
1646  // };
1647  //
1648  // MatrixDouble
1649  // diff_base_edge_functions(6,9*nb_gauss_pts*NBEDGE_AINSWORTH_HCURL(order));
1650  // double* diff_edge_n[] = {
1651  // &diff_base_edge_functions(0,0),
1652  // &diff_base_edge_functions(1,0),
1653  // &diff_base_edge_functions(2,0),
1654  // &diff_base_edge_functions(3,0),
1655  // &diff_base_edge_functions(4,0),
1656  // &diff_base_edge_functions(5,0)
1657  // };
1658  //
1659  // ierr = Hcurl_Ainsworth_EdgeBaseFunctions_MBTET(
1660  // edge_sense,
1661  // edge_order,
1662  // &*shape_fun.data().begin(),
1663  // diff_shape_fun,
1664  // edge_n,
1665  // diff_edge_n,
1666  // nb_gauss_pts,
1667  // Legendre_polynomials
1668  // );
1669  //
1670  //
1671  // for(int ee = 0;ee!=6;ee++) {
1672  // for(int ll = 0;ll!=NBEDGE_AINSWORTH_HCURL(order);ll++) {
1673  // std::ostringstream ss;
1674  // ss << "curl_edge_" << ee << "_" << ll;
1675  // Tag th;
1676  // rval = moab_ref.tag_get_handle(
1677  // ss.str().c_str(),3,MB_TYPE_DOUBLE,th,MB_TAG_CREAT|MB_TAG_SPARSE,def_val
1678  // ); CHKERRG(rval);
1679  // std::ostringstream ss_grad;
1680  // ss_grad << "grad_curl_edge_" << ee << "_" << ll;
1681  // Tag th_grad;
1682  // rval = moab_ref.tag_get_handle(
1683  // ss_grad.str().c_str(),9,MB_TYPE_DOUBLE,th_grad,MB_TAG_CREAT|MB_TAG_SPARSE,def_val
1684  // ); CHKERRG(rval);
1685  //
1686  // int gg = 0;
1687  // for(Range::iterator nit =
1688  // elem_nodes.begin();nit!=elem_nodes.end();nit++,gg++) {
1689  // rval = moab_ref.tag_set_data(
1690  // th,&*nit,1,&(edge_n[ee][gg*3*NBEDGE_AINSWORTH_HCURL(order)+ll*3])
1691  // ); CHKERRG(rval);
1692  // int sh = gg*9*NBEDGE_AINSWORTH_HCURL(order)+ll*9;
1693  // double grad[9] = {
1694  // diff_edge_n[ee][sh+0],diff_edge_n[ee][sh+3],diff_edge_n[ee][sh+6],
1695  // diff_edge_n[ee][sh+1],diff_edge_n[ee][sh+4],diff_edge_n[ee][sh+7],
1696  // diff_edge_n[ee][sh+2],diff_edge_n[ee][sh+5],diff_edge_n[ee][sh+8]
1697  // };
1698  // rval = moab_ref.tag_set_data(th_grad,&*nit,1,grad);
1699  // CHKERRG(rval);
1700  // }
1701  // }
1702  // }
1703 
1704  // cout << "NBFACETRI_AINSWORTH_EDGE_HCURL(order) " <<
1705  // NBFACETRI_AINSWORTH_EDGE_HCURL(order) << endl; MatrixDouble
1706  // base_face_edge_functions(
1707  // 4,3*3*NBFACETRI_AINSWORTH_EDGE_HCURL(order)*nb_gauss_pts
1708  // );
1709  // MatrixDouble diff_base_face_edge_functions(
1710  // 4,3*9*NBFACETRI_AINSWORTH_EDGE_HCURL(order)*nb_gauss_pts
1711  // );
1712  // double *phi_f_e[4][3];
1713  // double *diff_phi_f_e[4][3];
1714  // for(int ff = 0;ff!=4;ff++) {
1715  // for(int ee = 0;ee!=3;ee++) {
1716  // phi_f_e[ff][ee] =
1717  // &base_face_edge_functions(ff,ee*3*NBFACETRI_AINSWORTH_EDGE_HCURL(order)*nb_gauss_pts);
1718  // diff_phi_f_e[ff][ee] =
1719  // &diff_base_face_edge_functions(ff,ee*9*NBFACETRI_AINSWORTH_EDGE_HCURL(order)*nb_gauss_pts);
1720  // }
1721  // }
1722  //
1723  // ierr = Hcurl_Ainsworth_EdgeBasedFaceFunctions_MBTET(
1724  // faces_nodes,
1725  // faces_order,
1726  // &*shape_fun.data().begin(),
1727  // diff_shape_fun,
1728  // phi_f_e,
1729  // diff_phi_f_e,
1730  // nb_gauss_pts,
1731  // Legendre_polynomials
1732  // ); CHKERRG(ierr);
1733  //
1734  // for(int ff = 0;ff!=4;ff++) {
1735  // for(int ee = 0;ee!=3;ee++) {
1736  // for(int ll = 0;ll!=NBFACETRI_AINSWORTH_EDGE_HCURL(order);ll++) {
1737  // std::ostringstream ss;
1738  // ss << "curl_face_edge_" << ff << "_" << ee << "_" << ll;
1739  // Tag th;
1740  // rval = moab_ref.tag_get_handle(
1741  // ss.str().c_str(),3,MB_TYPE_DOUBLE,th,MB_TAG_CREAT|MB_TAG_SPARSE,def_val
1742  // ); CHKERRG(rval);
1743  //
1744  // std::ostringstream ss_grad;
1745  // ss_grad << "grad_curl_face_edge_" << ff << "_" << ee << "_" << ll;
1746  // Tag th_grad;
1747  // rval = moab_ref.tag_get_handle(
1748  // ss_grad.str().c_str(),9,MB_TYPE_DOUBLE,th_grad,MB_TAG_CREAT|MB_TAG_SPARSE,def_val
1749  // ); CHKERRG(rval);
1750  //
1751  // int gg = 0;
1752  // for(Range::iterator nit =
1753  // elem_nodes.begin();nit!=elem_nodes.end();nit++,gg++) {
1754  //
1755  // int idx =
1756  // 3*NBFACETRI_AINSWORTH_EDGE_HCURL(order)*gg+ll*3;
1757  // if(idx >= base_face_edge_functions.size2()) {
1758  // cerr << ff << " " << ee << " " << ll << " " << gg << endl;
1759  // }
1760  //
1761  // rval = moab_ref.tag_set_data(th,&*nit,1,&(phi_f_e[ff][ee][idx]));
1762  // CHKERRG(rval);
1763  //
1764  // int sh = gg*9*NBFACETRI_AINSWORTH_EDGE_HCURL(order)+ll*9;
1765  // double grad[9] = {
1766  // diff_phi_f_e[ff][ee][sh+0],diff_phi_f_e[ff][ee][sh+3],diff_phi_f_e[ff][ee][sh+6],
1767  // diff_phi_f_e[ff][ee][sh+1],diff_phi_f_e[ff][ee][sh+4],diff_phi_f_e[ff][ee][sh+7],
1768  // diff_phi_f_e[ff][ee][sh+2],diff_phi_f_e[ff][ee][sh+5],diff_phi_f_e[ff][ee][sh+8]
1769  // };
1770  // rval = moab_ref.tag_set_data(th_grad,&*nit,1,grad);
1771  // CHKERRG(rval);
1772  //
1773  // }
1774  // }
1775  // }
1776  // }
1777  //
1778  cout << "NBFACETRI_AINSWORTH_FACE_HCURL "
1780  MatrixDouble base_face_bubble_functions(
1781  4, 3 * NBFACETRI_AINSWORTH_FACE_HCURL(order) * nb_gauss_pts);
1782  MatrixDouble diff_base_face_bubble_functions(
1783  4, 9 * NBFACETRI_AINSWORTH_FACE_HCURL(order) * nb_gauss_pts);
1784  double *phi_f[4];
1785  double *diff_phi_f[4];
1786  for (int ff = 0; ff != 4; ff++) {
1787  phi_f[ff] = &base_face_bubble_functions(ff, 0);
1788  diff_phi_f[ff] = &diff_base_face_bubble_functions(ff, 0);
1789  }
1790 
1792  faces_nodes, faces_order, &*shape_fun.data().begin(), diff_shape_fun,
1793  phi_f, diff_phi_f, nb_gauss_pts, Legendre_polynomials);
1794  CHKERRG(ierr);
1795 
1796  for (int ff = 0; ff != 4; ff++) {
1797  for (int ll = 0; ll != NBFACETRI_AINSWORTH_FACE_HCURL(order); ll++) {
1798  std::ostringstream ss;
1799  ss << "curl_face_bubble_" << ff << "_" << ll;
1800  Tag th;
1801  rval = moab_ref.tag_get_handle(ss.str().c_str(), 3, MB_TYPE_DOUBLE, th,
1802  MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
1803  CHKERRG(rval);
1804  std::ostringstream grad_ss;
1805  grad_ss << "grad_curl_face_bubble_" << ff << "_" << ll;
1806  Tag th_grad;
1807  rval = moab_ref.tag_get_handle(grad_ss.str().c_str(), 9, MB_TYPE_DOUBLE,
1808  th_grad, MB_TAG_CREAT | MB_TAG_SPARSE,
1809  def_val);
1810  CHKERRG(rval);
1811 
1812  int gg = 0;
1813  for (Range::iterator nit = elem_nodes.begin(); nit != elem_nodes.end();
1814  nit++, gg++) {
1815  int idx = 3 * NBFACETRI_AINSWORTH_FACE_HCURL(order) * gg + ll * 3;
1816  rval = moab_ref.tag_set_data(th, &*nit, 1, &(phi_f[ff][idx]));
1817  CHKERRG(rval);
1818  int sh = gg * 9 * NBFACETRI_AINSWORTH_FACE_HCURL(order) + ll * 9;
1819  double grad[9] = {diff_phi_f[ff][sh + 0], diff_phi_f[ff][sh + 3],
1820  diff_phi_f[ff][sh + 6], diff_phi_f[ff][sh + 1],
1821  diff_phi_f[ff][sh + 4], diff_phi_f[ff][sh + 7],
1822  diff_phi_f[ff][sh + 2], diff_phi_f[ff][sh + 5],
1823  diff_phi_f[ff][sh + 8]};
1824  rval = moab_ref.tag_set_data(th_grad, &*nit, 1, grad);
1825  CHKERRG(rval);
1826  }
1827  }
1828  }
1829 
1830  // cout << "NBVOLUMETET_AINSWORTH_FACE_HCURL " <<
1831  // NBVOLUMETET_AINSWORTH_FACE_HCURL(order) << endl; VectorDouble
1832  // base_face_inetrior_functions(3*NBVOLUMETET_AINSWORTH_FACE_HCURL(order)*nb_gauss_pts);
1833  // VectorDouble
1834  // diff_base_face_inetrior_functions(9*NBVOLUMETET_AINSWORTH_FACE_HCURL(order)*nb_gauss_pts);
1835  // double *phi_v_f = &base_face_inetrior_functions[0];
1836  // double *diff_phi_v_f = &diff_base_face_inetrior_functions[0];
1837  // ierr = Hcurl_Ainsworth_FaceInteriorFunctions_MBTET(
1838  // faces_nodes,
1839  // order,
1840  // &*shape_fun.data().begin(),
1841  // diff_shape_fun,
1842  // phi_v_f,
1843  // diff_phi_v_f,
1844  // nb_gauss_pts,
1845  // Legendre_polynomials
1846  // ); CHKERRG(ierr);
1847  // for(int ll = 0;ll!=NBVOLUMETET_AINSWORTH_FACE_HCURL(order);ll++) {
1848  //
1849  // std::ostringstream ss;
1850  // ss << "curl_face_interior_" << ll;
1851  // Tag th;
1852  // rval = moab_ref.tag_get_handle(
1853  // ss.str().c_str(),3,MB_TYPE_DOUBLE,th,MB_TAG_CREAT|MB_TAG_SPARSE,def_val
1854  // ); CHKERRG(rval);
1855  //
1856  // std::ostringstream ss_grad;
1857  // ss_grad << "grad_curl_face_interior_" << ll;
1858  // Tag th_grad;
1859  // rval = moab_ref.tag_get_handle(
1860  // ss_grad.str().c_str(),9,MB_TYPE_DOUBLE,th_grad,MB_TAG_CREAT|MB_TAG_SPARSE,def_val
1861  // ); CHKERRG(rval);
1862  //
1863  //
1864  // int gg = 0;
1865  // for(Range::iterator nit =
1866  // elem_nodes.begin();nit!=elem_nodes.end();nit++,gg++) {
1867  // int idx = 3*NBVOLUMETET_AINSWORTH_FACE_HCURL(order)*gg+ll*3;
1868  // rval = moab_ref.tag_set_data(th,&*nit,1,&(phi_v_f[idx]));
1869  // CHKERRG(rval); int sh =
1870  // gg*9*NBVOLUMETET_AINSWORTH_FACE_HCURL(order)+ll*9; double grad[9] = {
1871  // diff_phi_v_f[sh+0],diff_phi_v_f[sh+3],diff_phi_v_f[sh+6],
1872  // diff_phi_v_f[sh+1],diff_phi_v_f[sh+4],diff_phi_v_f[sh+7],
1873  // diff_phi_v_f[sh+2],diff_phi_v_f[sh+5],diff_phi_v_f[sh+8]
1874  // };
1875  // rval = moab_ref.tag_set_data(th_grad,&*nit,1,grad); CHKERRG(rval);
1876  // }
1877  // }
1878 
1879  // cout << "NBVOLUMETET_AINSWORTH_TET_HCURL " <<
1880  // NBVOLUMETET_AINSWORTH_TET_HCURL(order) << endl; VectorDouble
1881  // base_interior_functions(3*NBVOLUMETET_AINSWORTH_TET_HCURL(order)*nb_gauss_pts);
1882  // VectorDouble
1883  // diff_base_interior_functions(9*NBVOLUMETET_AINSWORTH_TET_HCURL(order)*nb_gauss_pts);
1884  // double *phi_v = &base_interior_functions[0];
1885  // double *diff_phi_v = &diff_base_interior_functions[0];
1886  // ierr = Hcurl_Ainsworth_VolumeInteriorFunctions_MBTET(
1887  // order,
1888  // &*shape_fun.data().begin(),
1889  // diff_shape_fun,
1890  // phi_v,
1891  // diff_phi_v,
1892  // nb_gauss_pts,
1893  // Legendre_polynomials
1894  // ); CHKERRG(ierr);
1895  // for(int ll = 0;ll!=NBVOLUMETET_AINSWORTH_TET_HCURL(order);ll++) {
1896  //
1897  // std::ostringstream ss;
1898  // ss << "curl_interior_" << ll;
1899  // Tag th;
1900  // rval = moab_ref.tag_get_handle(
1901  // ss.str().c_str(),3,MB_TYPE_DOUBLE,th,MB_TAG_CREAT|MB_TAG_SPARSE,def_val
1902  // ); CHKERRG(rval);
1903  //
1904  // std::ostringstream ss_gard;
1905  // ss_gard << "grad_curl_interior_" << ll;
1906  // Tag th_grad;
1907  // rval = moab_ref.tag_get_handle(
1908  // ss_gard.str().c_str(),9,MB_TYPE_DOUBLE,th_grad,MB_TAG_CREAT|MB_TAG_SPARSE,def_val
1909  // ); CHKERRG(rval);
1910  //
1911  // int gg = 0;
1912  // for(Range::iterator nit =
1913  // elem_nodes.begin();nit!=elem_nodes.end();nit++,gg++) {
1914  // int idx = 3*NBVOLUMETET_AINSWORTH_TET_HCURL(order)*gg+ll*3;
1915  // rval = moab_ref.tag_set_data(th,&*nit,1,&(phi_v[idx]));
1916  // CHKERRG(rval); int sh =
1917  // gg*9*NBVOLUMETET_AINSWORTH_TET_HCURL(order)+ll*9; double grad[9] = {
1918  // diff_phi_v[sh+0],diff_phi_v[sh+3],diff_phi_v[sh+6],
1919  // diff_phi_v[sh+1],diff_phi_v[sh+4],diff_phi_v[sh+7],
1920  // diff_phi_v[sh+2],diff_phi_v[sh+5],diff_phi_v[sh+8]
1921  // };
1922  // rval = moab_ref.tag_set_data(th_grad,&*nit,1,grad); CHKERRG(rval);
1923  // }
1924  // }
1925 
1926  // cout << "NBFACETRI_AINSWORTH_HCURL(order) " <<
1927  // NBFACETRI_AINSWORTH_HCURL(order) << endl; MatrixDouble base_face_functions(
1928  // 4,3*NBFACETRI_AINSWORTH_HCURL(order)*nb_gauss_pts
1929  // );
1930  // MatrixDouble diff_base_face_functions(
1931  // 4,9*NBFACETRI_AINSWORTH_HCURL(order)*nb_gauss_pts
1932  // );
1933  // for(int ff=0;ff!=4;ff++) {
1934  // phi_f[ff] = &base_face_functions(ff,0);
1935  // diff_phi_f[ff] = &diff_base_face_functions(ff,0);
1936  // }
1937  // ierr = Hcurl_Ainsworth_FaceFunctions_MBTET(
1938  // faces_nodes,
1939  // faces_order,
1940  // &*shape_fun.data().begin(),
1941  // diff_shape_fun,
1942  // phi_f,
1943  // diff_phi_f,
1944  // nb_gauss_pts,
1945  // Legendre_polynomials
1946  // ); CHKERRG(ierr);
1947  // for(int ff = 0;ff!=4;ff++) {
1948  // for(int ll = 0;ll!=NBFACETRI_AINSWORTH_HCURL(order);ll++) {
1949  // std::ostringstream ss;
1950  // ss << "curl_face_" << ff << "_" << ll;
1951  // Tag th;
1952  // rval = moab_ref.tag_get_handle(
1953  // ss.str().c_str(),3,MB_TYPE_DOUBLE,th,MB_TAG_CREAT|MB_TAG_SPARSE,def_val
1954  // ); CHKERRG(rval);
1955  //
1956  // int gg = 0;
1957  // for(Range::iterator nit =
1958  // elem_nodes.begin();nit!=elem_nodes.end();nit++,gg++) {
1959  // int idx = 3*NBFACETRI_AINSWORTH_HCURL (order)*gg+ll*3;
1960  // rval = moab_ref.tag_set_data(th,&*nit,1,&(phi_f[ff][idx]));
1961  // CHKERRG(rval);
1962  // }
1963  // }
1964  // }
1965  //
1966  // cout << "NBVOLUMETET_AINSWORTH_TET_HCURL(order) " <<
1967  // NBVOLUMETET_AINSWORTH_HCURL(order) << endl; VectorDouble
1968  // base_volume_functions(3*NBVOLUMETET_AINSWORTH_HCURL(order)*nb_gauss_pts);
1969  // VectorDouble
1970  // diff_base_volume_functions(9*NBVOLUMETET_AINSWORTH_HCURL(order)*nb_gauss_pts); phi_v
1971  // = &base_volume_functions[0]; diff_phi_v = &diff_base_volume_functions[0];
1972  // ierr = MoFEM::Hcurl_Ainsworth_VolumeFunctions_MBTET(
1973  // order,
1974  // &*shape_fun.data().begin(),
1975  // diff_shape_fun,
1976  // phi_v,
1977  // diff_phi_v,
1978  // nb_gauss_pts,
1979  // Legendre_polynomials
1980  // ); CHKERRG(ierr);
1981  // for(int ll = 0;ll!=NBVOLUMETET_AINSWORTH_HCURL(order);ll++) {
1982  // std::ostringstream ss;
1983  // ss << "curl_volume_" << ll;
1984  // Tag th;
1985  // rval = moab_ref.tag_get_handle(
1986  // ss.str().c_str(),3,MB_TYPE_DOUBLE,th,MB_TAG_CREAT|MB_TAG_SPARSE,def_val
1987  // ); CHKERRG(rval);
1988  //
1989  // int gg = 0;
1990  // for(Range::iterator nit =
1991  // elem_nodes.begin();nit!=elem_nodes.end();nit++,gg++) {
1992  // int idx = 3*NBVOLUMETET_AINSWORTH_HCURL(order)*gg+ll*3;
1993  // rval = moab_ref.tag_set_data(th,&*nit,1,&(phi_v[idx]));
1994  // CHKERRG(rval);
1995  // }
1996  // }
1997 
1998  EntityHandle meshset;
1999  rval = moab_ref.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER, meshset);
2000  CHKERRG(rval);
2001  rval = moab_ref.add_entities(meshset, tets);
2002  CHKERRG(rval);
2003  rval = moab_ref.write_file(file_name.c_str(), "VTK", "", &meshset, 1);
2004  CHKERRG(rval);
2005 
2007 }
2008 
2009 #endif // GENERATE_VTK_WITH_CURL_BASE
2010 
2011 #ifndef GENERATE_VTK_WITH_CURL_BASE
2012 
2014  MoFEMErrorCode operator()(int shift, int p, int nb_integration_pts,
2015  int n0_idx, int n1_idx, double n[],
2016  FTensor::Tensor1<double, 3> t_grad_n[],
2018  FTensor::Tensor2<double *, 3, 3> &t_diff_phi) {
2019 
2022 
2024 
2025  FTensor::Tensor1<double, 3> &t_grad_n0 = t_grad_n[n0_idx];
2026  FTensor::Tensor1<double, 3> &t_grad_n1 = t_grad_n[n1_idx];
2027  FTensor::Tensor1<double, 3> t_grad_n0_p_n1;
2028  t_grad_n0_p_n1(i) = t_grad_n0(i) + t_grad_n1(i);
2030  FTensor::Tensor2<double, 3, 3> t_diff_phi_0;
2031 
2032  VectorDouble fi(p + 1);
2033  VectorDouble diff_fi(3 * p + 3);
2034 
2035  for (int gg = 0; gg != nb_integration_pts; ++gg) {
2036  const int shift_n = shift * gg;
2037  const double n0 = n[shift_n + n0_idx];
2038  const double n1 = n[shift_n + n1_idx];
2039 
2040  t_phi(i) = n0 * t_grad_n1(i) - n1 * t_grad_n0(i);
2041  t_diff_phi(i, j) =
2042  t_grad_n0(j) * t_grad_n1(i) - t_grad_n1(j) * t_grad_n0(i);
2043 
2044  t_phi_0(i) = t_phi(i);
2045  t_diff_phi_0(i, j) = t_diff_phi(i, j);
2046 
2047  ++t_phi;
2048  ++t_diff_phi;
2049 
2050  if (p > 1) {
2051 
2052  CHKERR Jacobi_polynomials(p, 0, n1, n0 + n1, &t_grad_n1(0),
2053  &t_grad_n0_p_n1(0), &*fi.data().begin(),
2054  &*diff_fi.data().begin(), 3);
2055 
2056  FTensor::Tensor1<double *, 3> t_diff_fi(&diff_fi[0 * (p + 1) + 1],
2057  &diff_fi[1 * (p + 1) + 1],
2058  &diff_fi[2 * (p + 1) + 1], 1);
2059 
2060  FTensor::Tensor1<double, 3> t_diff_b;
2061  for (int oo = 1; oo <= p - 1; ++oo) {
2062  const double b = pow(n0 + n1, oo);
2063  t_diff_b(i) =
2064  oo * pow(n0 + n1, oo - 1) * (t_grad_n0(i) + t_grad_n1(i));
2065  t_phi(i) = b * fi[oo] * t_phi_0(i);
2066  t_diff_phi(i, j) = (b * fi[oo]) * t_diff_phi_0(i, j) +
2067  (b * t_diff_fi(j)) * t_phi_0(i) +
2068  t_diff_b(j) * fi[oo] * t_phi_0(i);
2069 
2070  ++t_diff_fi;
2071  ++t_phi;
2072  ++t_diff_phi;
2073  }
2074  }
2075  }
2076 
2078  }
2079 
2081  int shift, int p, int nb_integration_pts, int n0_idx, int n1_idx,
2082  double n[], FTensor::Tensor1<double, 3> t_grad_n[],
2084  FTensor::Tensor2<FTensor::PackPtr<double *, 6>, 3, 2> &t_diff_phi) {
2085 
2088 
2090 
2091  FTensor::Tensor1<double, 3> &t_grad_n0 = t_grad_n[n0_idx];
2092  FTensor::Tensor1<double, 3> &t_grad_n1 = t_grad_n[n1_idx];
2093  FTensor::Tensor1<double, 3> t_grad_n0_p_n1;
2094  t_grad_n0_p_n1(i) = t_grad_n0(i) + t_grad_n1(i);
2096  FTensor::Tensor2<double, 3, 2> t_diff_phi_0;
2097 
2098  FTensor::Tensor1<double, 2> t_grad_n0_j(t_grad_n0(0), t_grad_n0(1));
2099  FTensor::Tensor1<double, 2> t_grad_n1_j(t_grad_n1(0), t_grad_n1(1));
2100 
2101  VectorDouble fi(p + 1);
2102  MatrixDouble diff_fi(2, p + 1);
2103 
2104  for (int gg = 0; gg != nb_integration_pts; ++gg) {
2105  const int shift_n = shift * gg;
2106  const double n0 = n[shift_n + n0_idx];
2107  const double n1 = n[shift_n + n1_idx];
2108 
2109  t_phi(i) = n0 * t_grad_n1(i) - n1 * t_grad_n0(i);
2110  t_diff_phi(i, j) =
2111  t_grad_n1(i) * t_grad_n0_j(j) - t_grad_n0(i) * t_grad_n1_j(j);
2112 
2113  t_phi_0(i) = t_phi(i);
2114  t_diff_phi_0(i, j) = t_diff_phi(i, j);
2115 
2116  ++t_phi;
2117  ++t_diff_phi;
2118 
2119  if (p > 1) {
2120 
2121  CHKERR Jacobi_polynomials(p, 0, n1, n0 + n1, &t_grad_n1(0),
2122  &t_grad_n0_p_n1(0), &*fi.data().begin(),
2123  &*diff_fi.data().begin(), 2);
2124 
2125  FTensor::Tensor1<double *, 2> t_diff_fi(&diff_fi(0, 1), &diff_fi(1, 1));
2126 
2127  // cerr << diff_fi << endl;
2128 
2129  FTensor::Tensor1<double, 2> t_diff_b;
2130  for (int oo = 1; oo <= p - 1; ++oo) {
2131  const double b = pow(n0 + n1, oo);
2132  t_diff_b(j) =
2133  oo * pow(n0 + n1, oo - 1) * (t_grad_n0_j(j) + t_grad_n1_j(j));
2134  t_phi(i) = b * fi[oo] * t_phi_0(i);
2135  t_diff_phi(i, j) = (b * fi[oo]) * t_diff_phi_0(i, j) +
2136  (b * t_phi_0(i) * t_diff_fi(j)) +
2137  fi[oo] * t_phi_0(i) * t_diff_b(j);
2138 
2139  ++t_diff_fi;
2140  ++t_phi;
2141  ++t_diff_phi;
2142  }
2143  }
2144  }
2145 
2147  }
2148 
2149  MoFEMErrorCode operator()(int shift, int p, int nb_integration_pts,
2150  int n0_idx, int n1_idx, double n[],
2151  FTensor::Tensor1<double, 3> t_grad_n[],
2153 
2155 
2157 
2158  FTensor::Tensor1<double, 3> &t_grad_n0 = t_grad_n[n0_idx];
2159  FTensor::Tensor1<double, 3> &t_grad_n1 = t_grad_n[n1_idx];
2161 
2162  double fi[p + 1];
2163 
2164  for (int gg = 0; gg != nb_integration_pts; ++gg) {
2165  const int shift_n = shift * gg;
2166  const double n0 = n[shift_n + n0_idx];
2167  const double n1 = n[shift_n + n1_idx];
2168 
2169  t_phi(i) = n0 * t_grad_n1(i) - n1 * t_grad_n0(i);
2170  t_phi_0(i) = t_phi(i);
2171 
2172  ++t_phi;
2173 
2174  if (p > 1) {
2175 
2176  CHKERR Jacobi_polynomials(p, 0, n1, n0 + n1, NULL, NULL, fi, NULL, 1);
2177 
2178  FTensor::Tensor1<double, 2> t_diff_b;
2179  for (int oo = 1; oo <= p - 1; ++oo) {
2180  const double b = pow(n0 + n1, oo);
2181  t_phi(i) = b * fi[oo] * t_phi_0(i);
2182  ++t_phi;
2183  }
2184  }
2185  }
2186 
2188  }
2189 };
2190 
2192  int *sense, int *p, double *n, double *diff_n, double *phi[],
2193  double *diff_phi[], int nb_integration_pts) {
2194 
2195  const int e_nodes[6][2] = {{0, 1}, {1, 2}, {2, 0}, {0, 3}, {1, 3}, {2, 3}};
2196 
2198 
2199  FTensor::Tensor1<double, 3> t_grad_n[4];
2200  for (int nn = 0; nn != 4; ++nn) {
2201  t_grad_n[nn] = FTensor::Tensor1<double, 3>(
2202  diff_n[3 * nn + 0], diff_n[3 * nn + 1], diff_n[3 * nn + 2]);
2203  // cerr << t_grad_n[nn](0) << " " << t_grad_n[nn](1) << " " << t_grad_n[nn](2)
2204  // << endl;
2205  };
2206 
2207 
2208  HcurlEdgeBase h_curl_base_on_edge;
2209 
2210  for (int ee = 0; ee != 6; ++ee) {
2211 
2212  FTensor::Tensor1<double *, 3> t_phi(&phi[ee][HVEC0], &phi[ee][HVEC1],
2213  &phi[ee][HVEC2], 3);
2214 
2216  &diff_phi[ee][HVEC0_0], &diff_phi[ee][HVEC0_1],
2217  &diff_phi[ee][HVEC0_2], &diff_phi[ee][HVEC1_0],
2218  &diff_phi[ee][HVEC1_1], &diff_phi[ee][HVEC1_2],
2219  &diff_phi[ee][HVEC2_0], &diff_phi[ee][HVEC2_1],
2220  &diff_phi[ee][HVEC2_2], 9);
2221 
2222  int n0_idx = e_nodes[ee][0];
2223  int n1_idx = e_nodes[ee][1];
2224  if(sense[ee]==-1) {
2225  int n_tmp = n0_idx;
2226  n0_idx = n1_idx;
2227  n1_idx = n_tmp;
2228  }
2229 
2230  CHKERR h_curl_base_on_edge(4, p[ee], nb_integration_pts, n0_idx, n1_idx, n,
2231  t_grad_n, t_phi, t_diff_phi);
2232  }
2233 
2235 }
2236 
2238  int *sense, int *p, double *n, double *diff_n, double *phi[],
2239  double *diff_phi[], int nb_integration_pts) {
2240 
2241  const int e_nodes[3][2] = {{0, 1}, {1, 2}, {2, 0}};
2242 
2244 
2245  FTensor::Tensor1<double, 3> t_grad_n[3];
2246  for (int nn = 0; nn != 3; ++nn) {
2247  t_grad_n[nn] =
2248  FTensor::Tensor1<double, 3>(diff_n[2 * nn + 0], diff_n[2 * nn + 1], 0.);
2249  };
2250 
2251  HcurlEdgeBase h_curl_base_on_edge;
2252 
2253  for (int ee = 0; ee != 3; ++ee) {
2254 
2255  FTensor::Tensor1<double *, 3> t_phi(&phi[ee][HVEC0], &phi[ee][HVEC1],
2256  &phi[ee][HVEC2], 3);
2257 
2259  &diff_phi[ee][HVEC0_0], &diff_phi[ee][HVEC0_1],
2260  &diff_phi[ee][HVEC1_0], &diff_phi[ee][HVEC1_1],
2261  &diff_phi[ee][HVEC2_0], &diff_phi[ee][HVEC2_1]);
2262 
2263  int n0_idx = e_nodes[ee][0];
2264  int n1_idx = e_nodes[ee][1];
2265  if(sense[ee]==-1) {
2266  int n_tmp = n0_idx;
2267  n0_idx = n1_idx;
2268  n1_idx = n_tmp;
2269  }
2270 
2271  CHKERR h_curl_base_on_edge(3, p[ee], nb_integration_pts, n0_idx, n1_idx, n,
2272  t_grad_n, t_phi, t_diff_phi);
2273  }
2274 
2276 }
2277 
2279  int sense, int p, double *n, double *diff_n, double *phi,
2280  double *diff_phi, int nb_integration_pts) {
2282 
2283  FTensor::Tensor1<double, 3> t_grad_n[2];
2284  for (int nn = 0; nn != 2; ++nn) {
2285  t_grad_n[nn] = FTensor::Tensor1<double, 3>(diff_n[nn], 0., 0.);
2286  };
2287 
2288  HcurlEdgeBase h_curl_base_on_edge;
2289 
2290  FTensor::Tensor1<double *, 3> t_phi(&phi[HVEC0], &phi[HVEC1],
2291  &phi[HVEC2], 3);
2292 
2293  if (diff_phi != NULL) {
2294  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
2295  "Not implemented direvatives for edge for Hcurl Demkowicz base");
2296  }
2297 
2298  int n0_idx = 0;
2299  int n1_idx = 1;
2300  if (sense == -1) {
2301  int n_tmp = n0_idx;
2302  n0_idx = n1_idx;
2303  n1_idx = n_tmp;
2304  }
2305 
2306  CHKERR h_curl_base_on_edge(1, p, nb_integration_pts, n0_idx, n1_idx, n,
2307  t_grad_n, t_phi);
2308 
2310 }
2311 
2314  VectorDouble f0PhiII,diffF0PhiII;
2315  VectorDouble f1PhiII,diffF1PhiII;
2316  VectorDouble iFiF0,diffIFiF0;
2317  VectorDouble iFiF1,diffIFiF1;
2318 
2320 
2321  MoFEMErrorCode operator()(int shift, int p, int nb_integration_pts,
2322  int n0f0_idx, int n1f0_idx, int n2f0_idx,
2323  double n[], FTensor::Tensor1<double, 3> t_grad_n[],
2325  FTensor::Tensor2<double *, 3, 3> &t_diff_phi) {
2326 
2328 
2330  f0PhiII.resize(3 * NBEDGE_DEMKOWICZ_HCURL(p) * nb_integration_pts,false);
2331  diffF0PhiII.resize(9 * NBEDGE_DEMKOWICZ_HCURL(p) * nb_integration_pts,false);
2332 
2333  // edge base for family I
2334  double *f0_phi_ii = &*f0PhiII.data().begin();
2335  double *diff_f0_phi_ii = &*diffF0PhiII.data().begin();
2336  FTensor::Tensor1<double *, 3> t_f0_phi_ii(
2337  &f0_phi_ii[HVEC0], &f0_phi_ii[HVEC1], &f0_phi_ii[HVEC2], 3);
2338  FTensor::Tensor2<double *, 3, 3> t_diff_f0_phi_ii(
2339  &diff_f0_phi_ii[HVEC0_0], &diff_f0_phi_ii[HVEC0_1],
2340  &diff_f0_phi_ii[HVEC0_2], &diff_f0_phi_ii[HVEC1_0],
2341  &diff_f0_phi_ii[HVEC1_1], &diff_f0_phi_ii[HVEC1_2],
2342  &diff_f0_phi_ii[HVEC2_0], &diff_f0_phi_ii[HVEC2_1],
2343  &diff_f0_phi_ii[HVEC2_2], 9);
2344  CHKERR hCurlBaseOnEdge(4, p-1, nb_integration_pts, n0f0_idx, n1f0_idx, n,
2345  t_grad_n, t_f0_phi_ii, t_diff_f0_phi_ii);
2346 
2347  FTensor::Tensor1<double, 3> &t_grad_n0f0 = t_grad_n[n0f0_idx];
2348  FTensor::Tensor1<double, 3> &t_grad_n1f0 = t_grad_n[n1f0_idx];
2349  FTensor::Tensor1<double, 3> &t_grad_n2f0 = t_grad_n[n2f0_idx];
2350  FTensor::Tensor1<double, 3> t_grad_n0f0_p_n1f0;
2351  t_grad_n0f0_p_n1f0(i) = t_grad_n0f0(i) + t_grad_n1f0(i) + t_grad_n2f0(i);
2352 
2353  iFiF0.resize(p + 1, false);
2354  diffIFiF0.resize(3 * p + 3, false);
2355  double *ifif0 = &*iFiF0.data().begin();
2356  double *diff_ifif0 = &*diffIFiF0.data().begin();
2357 
2358  for (int gg = 0; gg != nb_integration_pts; ++gg) {
2359 
2360  const int shift_n = shift * gg;
2361  const double n0f0 = n[shift_n + n0f0_idx];
2362  const double n1f0 = n[shift_n + n1f0_idx];
2363  const double n2f0 = n[shift_n + n2f0_idx];
2364 
2365  int phi_shift = 3 * NBEDGE_DEMKOWICZ_HCURL(p - 1) * gg;
2366  int diff_phi_shift = 9 * NBEDGE_DEMKOWICZ_HCURL(p - 1) * gg;
2367 
2368  for (int oo = 2; oo <= p; ++oo) {
2369 
2370  FTensor::Tensor1<double *, 3> t_f0_phi_ii(
2371  &f0_phi_ii[phi_shift + HVEC0], &f0_phi_ii[phi_shift + HVEC1],
2372  &f0_phi_ii[phi_shift + HVEC2], 3);
2373  FTensor::Tensor2<double *, 3, 3> t_diff_f0_phi_ii(
2374  &diff_f0_phi_ii[diff_phi_shift + HVEC0_0],
2375  &diff_f0_phi_ii[diff_phi_shift + HVEC0_1],
2376  &diff_f0_phi_ii[diff_phi_shift + HVEC0_2],
2377  &diff_f0_phi_ii[diff_phi_shift + HVEC1_0],
2378  &diff_f0_phi_ii[diff_phi_shift + HVEC1_1],
2379  &diff_f0_phi_ii[diff_phi_shift + HVEC1_2],
2380  &diff_f0_phi_ii[diff_phi_shift + HVEC2_0],
2381  &diff_f0_phi_ii[diff_phi_shift + HVEC2_1],
2382  &diff_f0_phi_ii[diff_phi_shift + HVEC2_2], 9);
2383 
2384  for (int ii = 0; ii <= oo - 2; ii++) {
2385 
2386  int jj = oo - 2 - ii;
2387 
2388  // family I
2390  jj + 1, 2 * ii + 1, n2f0, n0f0 + n1f0 + n2f0, &t_grad_n2f0(0),
2391  &t_grad_n0f0_p_n1f0(0), ifif0, diff_ifif0, 3);
2392  FTensor::Tensor1<double, 3> t_diff_ifif0(
2393  diff_ifif0[0 + jj], diff_ifif0[(jj + 1) + jj],
2394  diff_ifif0[2 * (jj + 1) + jj]);
2395  t_phi(i) = ifif0[jj] * t_f0_phi_ii(i);
2396  t_diff_phi(i, j) = ifif0[jj] * t_diff_f0_phi_ii(i, j) +
2397  t_diff_ifif0(j) * t_f0_phi_ii(i);
2398  ++t_phi;
2399  ++t_diff_phi;
2400  ++t_f0_phi_ii;
2401  ++t_diff_f0_phi_ii;
2402 
2403  }
2404  }
2405  }
2406 
2408  }
2409 
2410  MoFEMErrorCode operator()(int shift, int p, int nb_integration_pts,
2411  int n0f0_idx, int n1f0_idx, int n2f0_idx,
2412  int n0f1_idx, int n1f1_idx, int n2f1_idx,
2413  double n[], FTensor::Tensor1<double, 3> t_grad_n[],
2415  FTensor::Tensor2<double *, 3, 3> &t_diff_phi) {
2416 
2418 
2420 
2421  f0PhiII.resize(3 * NBEDGE_DEMKOWICZ_HCURL(p) * nb_integration_pts,false);
2422  diffF0PhiII.resize(9 * NBEDGE_DEMKOWICZ_HCURL(p) * nb_integration_pts,false);
2423  f1PhiII.resize(3 * NBEDGE_DEMKOWICZ_HCURL(p) * nb_integration_pts,false);
2424  diffF1PhiII.resize(9 * NBEDGE_DEMKOWICZ_HCURL(p) * nb_integration_pts,false);
2425 
2426  // edge base for family I
2427  double *f0_phi_ii = &*f0PhiII.data().begin();
2428  double *diff_f0_phi_ii = &*diffF0PhiII.data().begin();
2429  FTensor::Tensor1<double *, 3> t_f0_phi_ii(
2430  &f0_phi_ii[HVEC0], &f0_phi_ii[HVEC1], &f0_phi_ii[HVEC2], 3);
2431  FTensor::Tensor2<double *, 3, 3> t_diff_f0_phi_ii(
2432  &diff_f0_phi_ii[HVEC0_0], &diff_f0_phi_ii[HVEC0_1],
2433  &diff_f0_phi_ii[HVEC0_2], &diff_f0_phi_ii[HVEC1_0],
2434  &diff_f0_phi_ii[HVEC1_1], &diff_f0_phi_ii[HVEC1_2],
2435  &diff_f0_phi_ii[HVEC2_0], &diff_f0_phi_ii[HVEC2_1],
2436  &diff_f0_phi_ii[HVEC2_2], 9);
2437  CHKERR hCurlBaseOnEdge(4, p-1, nb_integration_pts, n0f0_idx, n1f0_idx, n,
2438  t_grad_n, t_f0_phi_ii, t_diff_f0_phi_ii);
2439 
2440  // edge base for family II
2441  double *f1_phi_ii = &*f1PhiII.data().begin();
2442  double *diff_f1_phi_ii = &*diffF1PhiII.data().begin();
2443  FTensor::Tensor1<double *, 3> t_f1_phi_ii(
2444  &f1_phi_ii[HVEC0], &f1_phi_ii[HVEC1], &f1_phi_ii[HVEC2], 3);
2445  FTensor::Tensor2<double *, 3, 3> t_diff_f1_phi_ii(
2446  &diff_f1_phi_ii[HVEC0_0], &diff_f1_phi_ii[HVEC0_1],
2447  &diff_f1_phi_ii[HVEC0_2], &diff_f1_phi_ii[HVEC1_0],
2448  &diff_f1_phi_ii[HVEC1_1], &diff_f1_phi_ii[HVEC1_2],
2449  &diff_f1_phi_ii[HVEC2_0], &diff_f1_phi_ii[HVEC2_1],
2450  &diff_f1_phi_ii[HVEC2_2], 9);
2451  CHKERR hCurlBaseOnEdge(4, p-1, nb_integration_pts, n0f1_idx, n1f1_idx, n,
2452  t_grad_n, t_f1_phi_ii, t_diff_f1_phi_ii);
2453 
2454  FTensor::Tensor1<double, 3> &t_grad_n0f0 = t_grad_n[n0f0_idx];
2455  FTensor::Tensor1<double, 3> &t_grad_n1f0 = t_grad_n[n1f0_idx];
2456  FTensor::Tensor1<double, 3> &t_grad_n2f0 = t_grad_n[n2f0_idx];
2457  FTensor::Tensor1<double, 3> t_grad_n0f0_p_n1f0;
2458  t_grad_n0f0_p_n1f0(i) = t_grad_n0f0(i) + t_grad_n1f0(i);
2459 
2460  FTensor::Tensor1<double, 3> &t_grad_n0f1 = t_grad_n[n0f1_idx];
2461  FTensor::Tensor1<double, 3> &t_grad_n1f1 = t_grad_n[n1f1_idx];
2462  FTensor::Tensor1<double, 3> &t_grad_n2f1 = t_grad_n[n2f1_idx];
2463  FTensor::Tensor1<double, 3> t_grad_n0f1_p_n1f1;
2464  t_grad_n0f1_p_n1f1(i) = t_grad_n0f1(i) + t_grad_n1f1(i);
2465 
2466  iFiF0.resize(p + 1, false);
2467  diffIFiF0.resize(3 * p + 3, false);
2468  double *ifif0 = &*iFiF0.data().begin();
2469  double *diff_ifif0 = &*diffIFiF0.data().begin();
2470  iFiF1.resize(p + 1, false);
2471  diffIFiF1.resize(3 * p + 3, false);
2472  double *ifif1 = &*iFiF1.data().begin();
2473  double *diff_ifif1 = &*diffIFiF1.data().begin();
2474 
2475  for (int gg = 0; gg != nb_integration_pts; ++gg) {
2476 
2477  const int shift_n = shift * gg;
2478  const double n0f0 = n[shift_n + n0f0_idx];
2479  const double n1f0 = n[shift_n + n1f0_idx];
2480  const double n2f0 = n[shift_n + n2f0_idx];
2481  const double n0f1 = n[shift_n + n0f1_idx];
2482  const double n1f1 = n[shift_n + n1f1_idx];
2483  const double n2f1 = n[shift_n + n2f1_idx];
2484 
2485  int phi_shift = 3 * NBEDGE_DEMKOWICZ_HCURL(p - 1) * gg;
2486  int diff_phi_shift = 9 * NBEDGE_DEMKOWICZ_HCURL(p - 1) * gg;
2487 
2488  int kk = 0;
2489  for (int oo = 2; oo <= p; ++oo) {
2490 
2491  FTensor::Tensor1<double *, 3> t_f0_phi_ii(
2492  &f0_phi_ii[phi_shift + HVEC0], &f0_phi_ii[phi_shift + HVEC1],
2493  &f0_phi_ii[phi_shift + HVEC2], 3);
2494  FTensor::Tensor2<double *, 3, 3> t_diff_f0_phi_ii(
2495  &diff_f0_phi_ii[diff_phi_shift + HVEC0_0],
2496  &diff_f0_phi_ii[diff_phi_shift + HVEC0_1],
2497  &diff_f0_phi_ii[diff_phi_shift + HVEC0_2],
2498  &diff_f0_phi_ii[diff_phi_shift + HVEC1_0],
2499  &diff_f0_phi_ii[diff_phi_shift + HVEC1_1],
2500  &diff_f0_phi_ii[diff_phi_shift + HVEC1_2],
2501  &diff_f0_phi_ii[diff_phi_shift + HVEC2_0],
2502  &diff_f0_phi_ii[diff_phi_shift + HVEC2_1],
2503  &diff_f0_phi_ii[diff_phi_shift + HVEC2_2], 9);
2504 
2505  FTensor::Tensor1<double *, 3> t_f1_phi_ii(
2506  &f1_phi_ii[phi_shift + HVEC0],
2507  &f1_phi_ii[phi_shift + HVEC1],
2508  &f1_phi_ii[phi_shift + HVEC2], 3);
2509  FTensor::Tensor2<double *, 3, 3> t_diff_f1_phi_ii(
2510  &diff_f1_phi_ii[diff_phi_shift + HVEC0_0],
2511  &diff_f1_phi_ii[diff_phi_shift + HVEC0_1],
2512  &diff_f1_phi_ii[diff_phi_shift + HVEC0_2],
2513  &diff_f1_phi_ii[diff_phi_shift + HVEC1_0],
2514  &diff_f1_phi_ii[diff_phi_shift + HVEC1_1],
2515  &diff_f1_phi_ii[diff_phi_shift + HVEC1_2],
2516  &diff_f1_phi_ii[diff_phi_shift + HVEC2_0],
2517  &diff_f1_phi_ii[diff_phi_shift + HVEC2_1],
2518  &diff_f1_phi_ii[diff_phi_shift + HVEC2_2], 9);
2519 
2520  for (int ii = 0; ii <= oo - 2; ii++) {
2521 
2522  int jj = oo - 2 - ii;
2523 
2524  // family I
2526  jj + 1, 2 * ii + 1, n2f0, n0f0 + n1f0, &t_grad_n2f0(0),
2527  &t_grad_n0f0_p_n1f0(0), ifif0, diff_ifif0, 3);
2528  FTensor::Tensor1<double, 3> t_diff_ifif0(
2529  diff_ifif0[0 + jj], diff_ifif0[(jj + 1) + jj],
2530  diff_ifif0[2 * (jj + 1) + jj]);
2531 
2532  t_phi(i) = ifif0[jj] * t_f0_phi_ii(i);
2533  t_diff_phi(i, j) = ifif0[jj] * t_diff_f0_phi_ii(i, j) +
2534  t_diff_ifif0(j) * t_f0_phi_ii(i);
2535 
2536  ++t_phi;
2537  ++t_diff_phi;
2538  ++t_f0_phi_ii;
2539  ++t_diff_f0_phi_ii;
2540  ++kk;
2541 
2542  // family II
2544  jj + 1, 2 * ii + 1, n2f1, n0f1 + n1f1, &t_grad_n2f1(0),
2545  &t_grad_n0f1_p_n1f1(0), ifif1, diff_ifif1, 3);
2546  FTensor::Tensor1<double, 3> t_diff_ifif1(
2547  diff_ifif1[0 + jj], diff_ifif1[(jj + 1) + jj],
2548  diff_ifif1[2 * (jj + 1) + jj]);
2549  t_phi(i) = ifif1[jj] * t_f1_phi_ii(i);
2550  t_diff_phi(i, j) = ifif1[jj] * t_diff_f1_phi_ii(i, j) +
2551  t_diff_ifif1(j) * t_f1_phi_ii(i);
2552  ++t_phi;
2553  ++t_diff_phi;
2554  ++t_f1_phi_ii;
2555  ++t_diff_f1_phi_ii;
2556  ++kk;
2557  }
2558  }
2559  if(kk != NBFACETRI_DEMKOWICZ_HCURL(p)) {
2560  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2561  "Wrong number of base functions");
2562  }
2563 
2564  }
2566  }
2567 
2569  int shift, int p, int nb_integration_pts, int n0f0_idx, int n1f0_idx,
2570  int n2f0_idx, int n0f1_idx, int n1f1_idx, int n2f1_idx, double n[],
2571  FTensor::Tensor1<double, 3> t_grad_n[],
2573  FTensor::Tensor2<FTensor::PackPtr<double *, 6>, 3, 2> &t_diff_phi) {
2575 
2577 
2578  f0PhiII.resize(3 * NBEDGE_DEMKOWICZ_HCURL(p) * nb_integration_pts, false);
2579  diffF0PhiII.resize(6 * NBEDGE_DEMKOWICZ_HCURL(p) * nb_integration_pts,
2580  false);
2581  f1PhiII.resize(3 * NBEDGE_DEMKOWICZ_HCURL(p) * nb_integration_pts, false);
2582  diffF1PhiII.resize(6 * NBEDGE_DEMKOWICZ_HCURL(p) * nb_integration_pts,
2583  false);
2584 
2585  // edge base for family I
2586  double *f0_phi_ii = &*f0PhiII.data().begin();
2587  double *diff_f0_phi_ii = &*diffF0PhiII.data().begin();
2588  FTensor::Tensor1<double *, 3> t_f0_phi_ii(
2589  &f0_phi_ii[HVEC0], &f0_phi_ii[HVEC1], &f0_phi_ii[HVEC2], 3);
2590  FTensor::Tensor2<FTensor::PackPtr<double *, 6>, 3, 2> t_diff_f0_phi_ii(
2591  &diff_f0_phi_ii[HVEC0_0], &diff_f0_phi_ii[HVEC0_1],
2592  &diff_f0_phi_ii[HVEC1_0], &diff_f0_phi_ii[HVEC1_1],
2593  &diff_f0_phi_ii[HVEC2_0], &diff_f0_phi_ii[HVEC2_1]);
2594  CHKERR hCurlBaseOnEdge(3, p - 1, nb_integration_pts, n0f0_idx, n1f0_idx, n,
2595  t_grad_n, t_f0_phi_ii, t_diff_f0_phi_ii);
2596 
2597  // edge base for family II
2598  double *f1_phi_ii = &*f1PhiII.data().begin();
2599  double *diff_f1_phi_ii = &*diffF1PhiII.data().begin();
2600  FTensor::Tensor1<double *, 3> t_f1_phi_ii(
2601  &f1_phi_ii[HVEC0], &f1_phi_ii[HVEC1], &f1_phi_ii[HVEC2], 3);
2602  FTensor::Tensor2<FTensor::PackPtr<double *, 6>, 3, 2> t_diff_f1_phi_ii(
2603  &diff_f1_phi_ii[HVEC0_0], &diff_f1_phi_ii[HVEC0_1],
2604  &diff_f1_phi_ii[HVEC1_0], &diff_f1_phi_ii[HVEC1_1],
2605  &diff_f1_phi_ii[HVEC2_0], &diff_f1_phi_ii[HVEC2_1]);
2606  CHKERR hCurlBaseOnEdge(3, p - 1, nb_integration_pts, n0f1_idx, n1f1_idx, n,
2607  t_grad_n, t_f1_phi_ii, t_diff_f1_phi_ii);
2608 
2609  FTensor::Tensor1<double, 3> &t_grad_n0f0 = t_grad_n[n0f0_idx];
2610  FTensor::Tensor1<double, 3> &t_grad_n1f0 = t_grad_n[n1f0_idx];
2611  FTensor::Tensor1<double, 3> &t_grad_n2f0 = t_grad_n[n2f0_idx];
2612  FTensor::Tensor1<double, 3> t_grad_n0f0_p_n1f0;
2613  t_grad_n0f0_p_n1f0(i) = t_grad_n0f0(i) + t_grad_n1f0(i);
2614 
2615  FTensor::Tensor1<double, 3> &t_grad_n0f1 = t_grad_n[n0f1_idx];
2616  FTensor::Tensor1<double, 3> &t_grad_n1f1 = t_grad_n[n1f1_idx];
2617  FTensor::Tensor1<double, 3> &t_grad_n2f1 = t_grad_n[n2f1_idx];
2618  FTensor::Tensor1<double, 3> t_grad_n0f1_p_n1f1;
2619  t_grad_n0f1_p_n1f1(i) = t_grad_n0f1(i) + t_grad_n1f1(i);
2620 
2621  iFiF0.resize((p+1), false);
2622  diffIFiF0.resize(2 * p + 2, false);
2623  double *ifif0 = &*iFiF0.data().begin();
2624  double *diff_ifif0 = &*diffIFiF0.data().begin();
2625  iFiF1.resize(p + 1, false);
2626  diffIFiF1.resize(2 * p + 2, false);
2627  double *ifif1 = &*iFiF1.data().begin();
2628  double *diff_ifif1 = &*diffIFiF1.data().begin();
2629 
2630  FTensor::Tensor1<double, 2> t_grad_n2f0_j(t_grad_n2f0(0), t_grad_n2f0(1));
2631 
2632  for (int gg = 0; gg != nb_integration_pts; ++gg) {
2633 
2634  const int shift_n = shift * gg;
2635  const double n0f0 = n[shift_n + n0f0_idx];
2636  const double n1f0 = n[shift_n + n1f0_idx];
2637  const double n2f0 = n[shift_n + n2f0_idx];
2638  const double n0f1 = n[shift_n + n0f1_idx];
2639  const double n1f1 = n[shift_n + n1f1_idx];
2640  const double n2f1 = n[shift_n + n2f1_idx];
2641 
2642  int phi_shift = 3 * NBEDGE_DEMKOWICZ_HCURL(p - 1) * gg;
2643  int diff_phi_shift = 6 * NBEDGE_DEMKOWICZ_HCURL(p - 1) * gg;
2644 
2645  int kk = 0;
2646  for (int oo = 2; oo <= p; ++oo) {
2647 
2648  FTensor::Tensor1<double *, 3> t_f0_phi_ii(
2649  &f0_phi_ii[phi_shift + HVEC0], &f0_phi_ii[phi_shift + HVEC1],
2650  &f0_phi_ii[phi_shift + HVEC2], 3);
2651  FTensor::Tensor2<FTensor::PackPtr<double *, 6>, 3, 2> t_diff_f0_phi_ii(
2652  &diff_f0_phi_ii[diff_phi_shift + HVEC0_0],
2653  &diff_f0_phi_ii[diff_phi_shift + HVEC0_1],
2654  &diff_f0_phi_ii[diff_phi_shift + HVEC1_0],
2655  &diff_f0_phi_ii[diff_phi_shift + HVEC1_1],
2656  &diff_f0_phi_ii[diff_phi_shift + HVEC2_0],
2657  &diff_f0_phi_ii[diff_phi_shift + HVEC2_1]);
2658 
2659  FTensor::Tensor1<double *, 3> t_f1_phi_ii(
2660  &f1_phi_ii[phi_shift + HVEC0], &f1_phi_ii[phi_shift + HVEC1],
2661  &f1_phi_ii[phi_shift + HVEC2], 3);
2662  FTensor::Tensor2<FTensor::PackPtr<double *, 6>, 3, 2> t_diff_f1_phi_ii(
2663  &diff_f1_phi_ii[diff_phi_shift + HVEC0_0],
2664  &diff_f1_phi_ii[diff_phi_shift + HVEC0_1],
2665  &diff_f1_phi_ii[diff_phi_shift + HVEC1_0],
2666  &diff_f1_phi_ii[diff_phi_shift + HVEC1_1],
2667  &diff_f1_phi_ii[diff_phi_shift + HVEC2_0],
2668  &diff_f1_phi_ii[diff_phi_shift + HVEC2_1]);
2669 
2670  for (int ii = 0; ii <= oo - 2; ii++) {
2671 
2672  int jj = oo - 2 - ii;
2673 
2674  // family I
2676  jj + 1, 2 * ii + 1, n2f0, n0f0 + n1f0, &t_grad_n2f0(0),
2677  &t_grad_n0f0_p_n1f0(0), ifif0, diff_ifif0, 2);
2678  FTensor::Tensor1<double, 2> t_diff_ifif0(
2679  diff_ifif0[0 + jj], diff_ifif0[(jj + 1) + jj]);
2680  t_phi(i) = ifif0[jj] * t_f0_phi_ii(i);
2681  t_diff_phi(i, j) = ifif0[jj] * t_diff_f0_phi_ii(i, j) +
2682  t_f0_phi_ii(i) * t_diff_ifif0(j);
2683 
2684  ++t_phi;
2685  ++t_diff_phi;
2686  ++t_f0_phi_ii;
2687  ++t_diff_f0_phi_ii;
2688  ++kk;
2689 
2690  // family II
2691  CHKERR IntegratedJacobi_polynomials(jj + 1, 2 * ii + 1, n2f1, n0f1+n1f1,
2692  &t_grad_n2f1(0), &t_grad_n0f1_p_n1f1(0),
2693  ifif1, diff_ifif1, 2);
2694  FTensor::Tensor1<double, 2> t_diff_ifif1(diff_ifif1[0 + jj],
2695  diff_ifif1[(jj + 1) + jj]);
2696  t_phi(i) = ifif1[jj] * t_f1_phi_ii(i);
2697  t_diff_phi(i, j) = ifif1[jj] * t_diff_f1_phi_ii(i, j) +
2698  t_f1_phi_ii(i) * t_diff_ifif1(j);
2699  ++t_phi;
2700  ++t_diff_phi;
2701  ++t_f1_phi_ii;
2702  ++t_diff_f1_phi_ii;
2703  ++kk;
2704  }
2705  }
2706  if(kk != NBFACETRI_DEMKOWICZ_HCURL(p)) {
2707  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2708  "Wrong number of base functions");
2709  }
2710  }
2712  }
2713 
2714 };
2715 
2718  double *n, double *diff_n,
2719  double *phi[],
2720  double *diff_phi[],
2721  int nb_integration_pts) {
2723 
2724  FTensor::Tensor1<double, 3> t_grad_n[4];
2725  for (int nn = 0; nn != 4; ++nn) {
2726  t_grad_n[nn] = FTensor::Tensor1<double, 3>(
2727  diff_n[3 * nn + 0], diff_n[3 * nn + 1], diff_n[3 * nn + 2]);
2728  };
2729 
2730  HcurlFaceBase h_curl_face_base;
2731 
2732  for (int ff = 0; ff != 4; ++ff) {
2733 
2734  if (p[ff] < 2)
2735  continue;
2736 
2737  FTensor::Tensor1<double *, 3> t_phi(&phi[ff][HVEC0], &phi[ff][HVEC1],
2738  &phi[ff][HVEC2], 3);
2739 
2741  &diff_phi[ff][HVEC0_0], &diff_phi[ff][HVEC0_1],
2742  &diff_phi[ff][HVEC0_2], &diff_phi[ff][HVEC1_0],
2743  &diff_phi[ff][HVEC1_1], &diff_phi[ff][HVEC1_2],
2744  &diff_phi[ff][HVEC2_0], &diff_phi[ff][HVEC2_1],
2745  &diff_phi[ff][HVEC2_2], 9);
2746 
2747  // f0, f1 - family I and family II
2748  const int n0f0_idx = faces_nodes[3 * ff + 0];
2749  const int n1f0_idx = faces_nodes[3 * ff + 1];
2750  const int n2f0_idx = faces_nodes[3 * ff + 2];
2751  // family II
2752  const int n0f1_idx = faces_nodes[3 * ff + 1];
2753  const int n1f1_idx = faces_nodes[3 * ff + 2];
2754  const int n2f1_idx = faces_nodes[3 * ff + 0];
2755 
2756  CHKERR h_curl_face_base(4, p[ff], nb_integration_pts, n0f0_idx, n1f0_idx,
2757  n2f0_idx, n0f1_idx, n1f1_idx, n2f1_idx, n, t_grad_n,
2758  t_phi, t_diff_phi);
2759  }
2760 
2762 }
2763 
2766  double *n, double *diff_n,
2767  double *phi,
2768  double *diff_phi,
2769  int nb_integration_pts) {
2771 
2772  FTensor::Tensor1<double, 3> t_grad_n[3];
2773  for (int nn = 0; nn != 3; ++nn) {
2774  t_grad_n[nn] = FTensor::Tensor1<double, 3>(
2775  diff_n[2 * nn + 0], diff_n[2 * nn + 1], 0.);
2776  };
2777 
2778  HcurlFaceBase h_curl_face_base;
2779 
2780  if (p < 2)
2782 
2783  FTensor::Tensor1<double *, 3> t_phi(&phi[HVEC0], &phi[HVEC1],
2784  &phi[HVEC2], 3);
2785 
2787  &diff_phi[HVEC0_0], &diff_phi[HVEC0_1], &diff_phi[HVEC1_0],
2788  &diff_phi[HVEC1_1], &diff_phi[HVEC2_0], &diff_phi[HVEC2_1]);
2789 
2790  // f0, f1 - family I and family II
2791  const int n0f0_idx = faces_nodes[0];
2792  const int n1f0_idx = faces_nodes[1];
2793  const int n2f0_idx = faces_nodes[2];
2794  // family II
2795  const int n0f1_idx = faces_nodes[1];
2796  const int n1f1_idx = faces_nodes[2];
2797  const int n2f1_idx = faces_nodes[0];
2798 
2799  CHKERR h_curl_face_base(3, p, nb_integration_pts, n0f0_idx, n1f0_idx,
2800  n2f0_idx, n0f1_idx, n1f1_idx, n2f1_idx, n, t_grad_n,
2801  t_phi, t_diff_phi);
2802 
2804 }
2805 
2807  int p, double *n, double *diff_n, double *phi, double *diff_phi,
2808  int nb_integration_pts) {
2809 
2810  int family[3][4] = {{0, 1, 2, 3}, {1, 2, 3, 0}, {2, 3, 0, 1}};
2813 
2815 
2816  if (p < 3)
2818 
2819  FTensor::Tensor1<double *, 3> t_phi(&phi[HVEC0], &phi[HVEC1], &phi[HVEC2],
2820  3);
2821 
2823  &diff_phi[HVEC0_0], &diff_phi[HVEC0_1], &diff_phi[HVEC0_2],
2824  &diff_phi[HVEC1_0], &diff_phi[HVEC1_1], &diff_phi[HVEC1_2],
2825  &diff_phi[HVEC2_0], &diff_phi[HVEC2_1], &diff_phi[HVEC2_2], 9);
2826 
2827  FTensor::Tensor1<double, 3> t_grad_n[4];
2828  for (int nn = 0; nn != 4; ++nn) {
2829  t_grad_n[nn] = FTensor::Tensor1<double, 3>(
2830  diff_n[3 * nn + 0], diff_n[3 * nn + 1], diff_n[3 * nn + 2]);
2831  };
2832 
2833  int nb_face_functions = (NBFACETRI_DEMKOWICZ_HCURL(p-1)) / 2;
2834  MatrixDouble phi_ij(3, 3 * nb_face_functions * nb_integration_pts);
2835  MatrixDouble diff_phi_ij(3, 9 * nb_face_functions * nb_integration_pts);
2836  MatrixDouble fi_k(3, p + 1);
2837  MatrixDouble diff_fi_k(3, 3 * p + 3);
2838  HcurlFaceBase h_curl_face_base;
2839 
2840  // calate face base for each family
2841  for (int ff = 0; ff != 3; ++ff) {
2842  double *phi_ij_ptr = &phi_ij(ff,0);
2843  double *diff_phi_ij_ptr = &diff_phi_ij(ff,0);
2844 
2846  &phi_ij_ptr[HVEC0], &phi_ij_ptr[HVEC1], &phi_ij_ptr[HVEC2], 3);
2847 
2848  FTensor::Tensor2<double *, 3, 3> t_diff_phi_ij(
2849  &diff_phi_ij_ptr[HVEC0_0], &diff_phi_ij_ptr[HVEC0_1],
2850  &diff_phi_ij_ptr[HVEC0_2], &diff_phi_ij_ptr[HVEC1_0],
2851  &diff_phi_ij_ptr[HVEC1_1], &diff_phi_ij_ptr[HVEC1_2],
2852  &diff_phi_ij_ptr[HVEC2_0], &diff_phi_ij_ptr[HVEC2_1],
2853  &diff_phi_ij_ptr[HVEC2_2], 9);
2854 
2855  const int n0_idx = family[ff][0];
2856  const int n1_idx = family[ff][1];
2857  const int n2_idx = family[ff][2];
2858 
2859  CHKERR h_curl_face_base(4, p - 1, nb_integration_pts, n0_idx, n1_idx,
2860  n2_idx, n, t_grad_n, t_phi_ij, t_diff_phi_ij);
2861  }
2862 
2863  FTensor::Tensor1<double, 3> &t_grad_n3f0 = t_grad_n[family[0][3]];
2864  FTensor::Tensor1<double, 3> &t_grad_n3f1 = t_grad_n[family[1][3]];
2865  FTensor::Tensor1<double, 3> &t_grad_n3f2 = t_grad_n[family[2][3]];
2866 
2867  FTensor::Tensor1<double, 3> t_sum_f0;
2868  t_sum_f0(i) = -t_grad_n3f0(i);
2869  FTensor::Tensor1<double, 3> t_sum_f1;
2870  t_sum_f1(i) = -t_grad_n3f1(i);
2871  FTensor::Tensor1<double, 3> t_sum_f2;
2872  t_sum_f2(i) = -t_grad_n3f2(i);
2873 
2874  for (int gg = 0; gg != nb_integration_pts; ++gg) {
2875 
2876  int shift_n = 4 * gg;
2877 
2878  double n3f0 = n[shift_n + family[0][3]];
2879  double n3f1 = n[shift_n + family[1][3]];
2880  double n3f2 = n[shift_n + family[2][3]];
2881 
2882  int kk = 0;
2883  for (int oo = 3; oo <= p; ++oo) {
2884 
2885  int phi_shift = 3 * nb_face_functions * gg;
2886  int diff_phi_shift = 9 * nb_face_functions * gg;
2887 
2888  FTensor::Tensor1<double *, 3> t_phi_face_f0(
2889  &phi_ij(0, phi_shift + HVEC0), &phi_ij(0, phi_shift + HVEC1),
2890  &phi_ij(0, phi_shift + HVEC2), 3);
2891  FTensor::Tensor2<double *, 3, 3> t_diff_phi_face_f0(
2892  &diff_phi_ij(0, diff_phi_shift + HVEC0_0),
2893  &diff_phi_ij(0, diff_phi_shift + HVEC0_1),
2894  &diff_phi_ij(0, diff_phi_shift + HVEC0_2),
2895  &diff_phi_ij(0, diff_phi_shift + HVEC1_0),
2896  &diff_phi_ij(0, diff_phi_shift + HVEC1_1),
2897  &diff_phi_ij(0, diff_phi_shift + HVEC1_2),
2898  &diff_phi_ij(0, diff_phi_shift + HVEC2_0),
2899  &diff_phi_ij(0, diff_phi_shift + HVEC2_1),
2900  &diff_phi_ij(0, diff_phi_shift + HVEC2_2), 9);
2901 
2902  FTensor::Tensor1<double *, 3> t_phi_face_f1(
2903  &phi_ij(1, phi_shift + HVEC0), &phi_ij(1, phi_shift + HVEC1),
2904  &phi_ij(1, phi_shift + HVEC2), 3);
2905  FTensor::Tensor2<double *, 3, 3> t_diff_phi_face_f1(
2906  &diff_phi_ij(1, diff_phi_shift + HVEC0_0),
2907  &diff_phi_ij(1, diff_phi_shift + HVEC0_1),
2908  &diff_phi_ij(1, diff_phi_shift + HVEC0_2),
2909  &diff_phi_ij(1, diff_phi_shift + HVEC1_0),
2910  &diff_phi_ij(1, diff_phi_shift + HVEC1_1),
2911  &diff_phi_ij(1, diff_phi_shift + HVEC1_2),
2912  &diff_phi_ij(1, diff_phi_shift + HVEC2_0),
2913  &diff_phi_ij(1, diff_phi_shift + HVEC2_1),
2914  &diff_phi_ij(1, diff_phi_shift + HVEC2_2), 9);
2915 
2916  FTensor::Tensor1<double *, 3> t_phi_face_f2(
2917  &phi_ij(2, phi_shift + HVEC0), &phi_ij(2, phi_shift + HVEC1),
2918  &phi_ij(2, phi_shift + HVEC2), 3);
2919  FTensor::Tensor2<double *, 3, 3> t_diff_phi_face_f2(
2920  &diff_phi_ij(2, diff_phi_shift + HVEC0_0),
2921  &diff_phi_ij(2, diff_phi_shift + HVEC0_1),
2922  &diff_phi_ij(2, diff_phi_shift + HVEC0_2),
2923  &diff_phi_ij(2, diff_phi_shift + HVEC1_0),
2924  &diff_phi_ij(2, diff_phi_shift + HVEC1_1),
2925  &diff_phi_ij(2, diff_phi_shift + HVEC1_2),
2926  &diff_phi_ij(2, diff_phi_shift + HVEC2_0),
2927  &diff_phi_ij(2, diff_phi_shift + HVEC2_1),
2928  &diff_phi_ij(2, diff_phi_shift + HVEC2_2), 9);
2929 
2930  int ij = 0;
2931  for (int oo_ij = 2; oo_ij != oo; ++oo_ij) {
2932  int k = oo-oo_ij;
2933 
2935  k, 2 * oo_ij, n3f0, 1 - n3f0, &t_grad_n3f0(0),
2936  &t_sum_f0(0), &fi_k(0, 0), &diff_fi_k(0, 0), 3);
2938  k, 2 * oo_ij, n3f1, 1 - n3f1, &t_grad_n3f1(0),
2939  &t_sum_f1(0), &fi_k(1, 0), &diff_fi_k(1, 0), 3);
2941  k, 2 * oo_ij, n3f2, 1 - n3f2, &t_grad_n3f2(0),
2942  &t_sum_f2(0), &fi_k(2, 0), &diff_fi_k(2, 0), 3);
2943 
2944  FTensor::Tensor1<double, 3> t_diff_fi_k_f0(diff_fi_k(0, 0 + k - 1),
2945  diff_fi_k(0, k + k - 1),
2946  diff_fi_k(0, 2 * k + k - 1));
2947  FTensor::Tensor1<double, 3> t_diff_fi_k_f1(diff_fi_k(1, 0 + k - 1),
2948  diff_fi_k(1, k + k - 1),
2949  diff_fi_k(1, 2 * k + k - 1));
2950  FTensor::Tensor1<double, 3> t_diff_fi_k_f2(diff_fi_k(2, 0 + k - 1),
2951  diff_fi_k(2, k + k - 1),
2952  diff_fi_k(2, 2 * k + k - 1));
2953 
2954  for (; ij != NBFACETRI_DEMKOWICZ_HCURL(oo_ij) / 2; ++ij) {
2955  t_phi(i) = fi_k(0, k - 1) * t_phi_face_f0(i);
2956  t_diff_phi(i, j) = t_diff_fi_k_f0(j) * t_phi_face_f0(i) +
2957  fi_k(0, k - 1) * t_diff_phi_face_f0(i, j);
2958  ++t_phi;
2959  ++t_diff_phi;
2960  ++t_phi_face_f0;
2961  ++t_diff_phi_face_f0;
2962  ++kk;
2963 
2964  t_phi(i) = fi_k(1, k - 1) * t_phi_face_f1(i);
2965  t_diff_phi(i, j) = t_diff_fi_k_f1(j) * t_phi_face_f1(i) +
2966  fi_k(1, k - 1) * t_diff_phi_face_f1(i, j);
2967  ++t_phi;
2968  ++t_diff_phi;
2969  ++t_phi_face_f1;
2970  ++t_diff_phi_face_f1;
2971  ++kk;
2972 
2973  t_phi(i) = fi_k(2, k - 1) * t_phi_face_f2(i);
2974  t_diff_phi(i, j) = t_diff_fi_k_f2(j) * t_phi_face_f2(i) +
2975  fi_k(2, k - 1) * t_diff_phi_face_f2(i, j);
2976  ++t_phi;
2977  ++t_diff_phi;
2978  ++t_phi_face_f2;
2979  ++t_diff_phi_face_f2;
2980  ++kk;
2981  }
2982  }
2983  }
2984  if (kk != NBVOLUMETET_DEMKOWICZ_HCURL(p)) {
2985  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2986  "Wrong number of base functions");
2987  }
2988  }
2989 
2991 }
2992 
2993 #endif
2994 
2995 #ifdef GENERATE_VTK_WITH_CURL_BASE
2996 
2997 MoFEMErrorCode VTK_Demkowicz_Hcurl_MBTET(const string file_name) {
2999 
3000  double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
3001 
3002  moab::Core core_ref;
3003  moab::Interface &moab_ref = core_ref;
3004 
3005  EntityHandle nodes[4];
3006  for (int nn = 0; nn < 4; nn++) {
3007  CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
3008  }
3009  EntityHandle tet;
3010  CHKERR moab_ref.create_element(MBTET, nodes, 4, tet);
3011 
3012  MoFEM::Core m_core_ref(moab_ref, PETSC_COMM_SELF, -2);
3013  MoFEM::Interface &m_field_ref = m_core_ref;
3014 
3015  CHKERR m_field_ref.getInterface<BitRefManager>()->setBitRefLevelByDim(
3016  0, 3, BitRefLevel().set(0));
3017 
3018  const int max_level = 3;
3019  for (int ll = 0; ll != max_level; ll++) {
3020  Range edges;
3021  CHKERR m_field_ref.getInterface<BitRefManager>()
3022  ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
3023  BitRefLevel().set(), MBEDGE, edges);
3024  Range tets;
3025  CHKERR m_field_ref.getInterface<BitRefManager>()
3026  ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
3027  BitRefLevel(ll).set(), MBTET, tets);
3028  // refine mesh
3029  MeshRefinement *m_ref;
3030  CHKERR m_field_ref.getInterface(m_ref);
3031  CHKERR m_ref->add_vertices_in_the_middel_of_edges(edges,
3032  BitRefLevel().set(ll + 1));
3033  CHKERR m_ref->refine_TET(tets, BitRefLevel().set(ll + 1));
3034  }
3035 
3036  Range tets;
3037  CHKERR m_field_ref.getInterface<BitRefManager>()
3038  ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(max_level),
3039  BitRefLevel().set(max_level), MBTET, tets);
3040 
3041  // Use 10 node tets to print base
3042  if (1) {
3043  EntityHandle meshset;
3044  CHKERR moab_ref.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER, meshset);
3045  CHKERR moab_ref.add_entities(meshset, tets);
3046  CHKERR moab_ref.convert_entities(meshset, true, false, false);
3047  CHKERR moab_ref.delete_entities(&meshset, 1);
3048  }
3049 
3050  Range elem_nodes;
3051  rval = moab_ref.get_connectivity(tets, elem_nodes, false);
3052  CHKERRG(rval);
3053 
3054  const int nb_gauss_pts = elem_nodes.size();
3055  MatrixDouble gauss_pts(nb_gauss_pts, 4);
3056  gauss_pts.clear();
3057  Range::iterator nit = elem_nodes.begin();
3058  for (int gg = 0; nit != elem_nodes.end(); nit++, gg++) {
3059  rval = moab_ref.get_coords(&*nit, 1, &gauss_pts(gg, 0));
3060  CHKERRG(rval);
3061  }
3062  gauss_pts = trans(gauss_pts);
3063 
3064  MatrixDouble shape_fun;
3065  shape_fun.resize(nb_gauss_pts, 4);
3066  CHKERR ShapeMBTET(&*shape_fun.data().begin(), &gauss_pts(0, 0),
3067  &gauss_pts(1, 0), &gauss_pts(2, 0), nb_gauss_pts);
3068 
3069  double diff_shape_fun[12];
3070  CHKERR ShapeDiffMBTET(diff_shape_fun);
3071 
3072  int edge_sense[6] = {1, 1, 1, 1, 1, 1};
3073  const int order = 4;
3074  int edge_order[6] = {order, order, order, order, order, order};
3075 
3076  MatrixDouble edge_phi(6, 3 * NBEDGE_DEMKOWICZ_HCURL(order) * nb_gauss_pts);
3077  MatrixDouble edge_diff_phi(6,
3078  9 * NBEDGE_DEMKOWICZ_HCURL(order) * nb_gauss_pts);
3079 
3080  edge_phi.clear();
3081  edge_diff_phi.clear();
3082 
3083  double *edge_phi_ptr[] = {&edge_phi(0, 0), &edge_phi(1, 0), &edge_phi(2, 0),
3084  &edge_phi(3, 0), &edge_phi(4, 0), &edge_phi(5, 0)};
3085  double *edge_diff_phi_ptr[] = {&edge_diff_phi(0, 0), &edge_diff_phi(1, 0),
3086  &edge_diff_phi(2, 0), &edge_diff_phi(3, 0),
3087  &edge_diff_phi(4, 0), &edge_diff_phi(5, 0)};
3088 
3090  edge_sense, edge_order, &*shape_fun.data().begin(), diff_shape_fun,
3091  edge_phi_ptr, edge_diff_phi_ptr, nb_gauss_pts);
3092 
3093  // cerr << edge_phi << endl;
3094 
3095  for (int ee = 0; ee != 6; ++ee) {
3096  for (int ll = 0; ll != NBEDGE_DEMKOWICZ_HCURL(order); ++ll) {
3097  double def_val[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
3098  std::string tag_name = "E" + boost::lexical_cast<std::string>(ee) + "_" +
3099  boost::lexical_cast<std::string>(ll);
3100  Tag th;
3101  CHKERR moab_ref.tag_get_handle(
3102  tag_name.c_str(),
3103  3, MB_TYPE_DOUBLE, th, MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
3104 
3105  int gg = 0;
3106  for (Range::iterator nit = elem_nodes.begin(); nit != elem_nodes.end();
3107  nit++, gg++) {
3108  int idx = 3 * NBEDGE_DEMKOWICZ_HCURL(order) * gg + 3 * ll;
3109  CHKERR moab_ref.tag_set_data(th, &*nit, 1, &edge_phi(ee,idx));
3110  }
3111  }
3112  }
3113 
3114  int faces_order[] = {order, order, order, order};
3115  int faces_nodes[] = {0, 1, 3, 1, 2, 3, 0, 2, 3, 0, 1, 2};
3116 
3117  MatrixDouble face_phi(4, 3 * NBFACETRI_DEMKOWICZ_HCURL(order) * nb_gauss_pts);
3118  MatrixDouble face_diff_phi(4, 9 * NBFACETRI_DEMKOWICZ_HCURL(order) *
3119  nb_gauss_pts);
3120  face_phi.clear();
3121  face_diff_phi.clear();
3122 
3123  double *face_phi_ptr[] = {&face_phi(0, 0), &face_phi(1, 0), &face_phi(2, 0),
3124  &face_phi(3, 0)};
3125  double *face_diff_phi_ptr[] = {&face_diff_phi(0, 0), &face_diff_phi(1, 0),
3126  &face_diff_phi(2, 0), &face_diff_phi(3, 0)};
3127 
3129  faces_nodes, faces_order, &*shape_fun.data().begin(), diff_shape_fun,
3130  face_phi_ptr, face_diff_phi_ptr, nb_gauss_pts);
3131 
3132  for (int ff = 0; ff != 4; ++ff) {
3133  for (int ll = 0; ll != NBFACETRI_DEMKOWICZ_HCURL(order); ++ll) {
3134  double def_val[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
3135  std::string tag_name = "F" + boost::lexical_cast<std::string>(ff) + "_" +
3136  boost::lexical_cast<std::string>(ll);
3137  Tag th;
3138  CHKERR moab_ref.tag_get_handle(tag_name.c_str(), 3, MB_TYPE_DOUBLE, th,
3139  MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
3140 
3141  int gg = 0;
3142  for (Range::iterator nit = elem_nodes.begin(); nit != elem_nodes.end();
3143  nit++, gg++) {
3144  int idx = 3 * NBFACETRI_DEMKOWICZ_HCURL(order) * gg + 3 * ll;
3145  CHKERR moab_ref.tag_set_data(th, &*nit, 1, &face_phi(ff, idx));
3146  }
3147  }
3148  }
3149 
3150  MatrixDouble vol_phi(nb_gauss_pts, 3 * NBVOLUMETET_DEMKOWICZ_HCURL(order));
3151  MatrixDouble diff_vol_phi(nb_gauss_pts,
3153 
3155  order, &*shape_fun.data().begin(), diff_shape_fun, &vol_phi(0,0),
3156  &diff_vol_phi(0,0), nb_gauss_pts);
3157 
3158  for (int ll = 0; ll != NBVOLUMETET_DEMKOWICZ_HCURL(order); ++ll) {
3159  double def_val[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
3160  std::string tag_name = "V_" + boost::lexical_cast<std::string>(ll);
3161  Tag th;
3162  CHKERR moab_ref.tag_get_handle(tag_name.c_str(), 3, MB_TYPE_DOUBLE, th,
3163  MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
3164 
3165  int gg = 0;
3166  for (Range::iterator nit = elem_nodes.begin(); nit != elem_nodes.end();
3167  nit++, gg++) {
3168  int idx = 3 * ll;
3169  CHKERR moab_ref.tag_set_data(th, &*nit, 1, &vol_phi(gg, idx));
3170  }
3171  }
3172 
3173  EntityHandle meshset;
3174  CHKERR moab_ref.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER, meshset);
3175  CHKERR moab_ref.add_entities(meshset, tets);
3176  CHKERR moab_ref.write_file(file_name.c_str(), "VTK", "", &meshset, 1);
3177 
3179 }
3180 
3181 static char help[] = "...\n\n";
3182 
3183 int main(int argc, char *argv[]) {
3184 
3185  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
3186 
3187  try {
3188  // CHKERR VTK_Ainsworth_Hcurl_MBTET("out_curl_vtk_ainsworth_base_on_tet.vtk");
3189  CHKERR VTK_Demkowicz_Hcurl_MBTET("out_curl_vtk_demkowicz_base_on_tet.vtk");
3190  }
3191  CATCH_ERRORS;
3192 
3194 
3195  return 0;
3196 }
3197 
3198 #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.
virtual MoFEMErrorCode refine_TET(const EntityHandle meshset, const BitRefLevel &bit, const bool respect_interface=false, int verb=QUIET, Range *ref_edges=NULL)
refine TET in the meshset
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:1251
#define NBVOLUMETET_AINSWORTH_FACE_HCURL(P)
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
virtual MoFEMErrorCode add_vertices_in_the_middel_of_edges(const EntityHandle meshset, const BitRefLevel &bit, const bool recursive=false, int verb=QUIET, EntityHandle start_v=0)
make vertices in the middle of edges in meshset and add them to refinement levels defined by bit
VectorDouble f1PhiII
Definition: Hcurl.cpp:2315
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:500
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:233
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:2806
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:77
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:2765
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:476
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:543
MoFEMErrorCode operator()(int shift, 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< double *, 3 > &t_phi, FTensor::Tensor2< double *, 3, 3 > &t_diff_phi)
Definition: Hcurl.cpp:2321
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:454
MoFEMErrorCode operator()(int shift, int p, int nb_integration_pts, int n0_idx, int n1_idx, double n[], FTensor::Tensor1< double, 3 > t_grad_n[], FTensor::Tensor1< double *, 3 > &t_phi)
Definition: Hcurl.cpp:2149
Core (interface) class.
Definition: Core.hpp:50
MoFEMErrorCode operator()(int shift, int p, int nb_integration_pts, int n0_idx, int n1_idx, double n[], FTensor::Tensor1< double, 3 > t_grad_n[], FTensor::Tensor1< double *, 3 > &t_phi, FTensor::Tensor2< FTensor::PackPtr< double *, 6 >, 3, 2 > &t_diff_phi)
Definition: Hcurl.cpp:2080
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:507
HcurlEdgeBase hCurlBaseOnEdge
Definition: Hcurl.cpp:2313
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:908
T data[Tensor_Dim]
FTensor::Index< 'i', 3 > i
Definition: Hcurl.cpp:2319
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
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:84
MoFEMErrorCode operator()(int shift, int p, int nb_integration_pts, int n0_idx, int n1_idx, double n[], FTensor::Tensor1< double, 3 > t_grad_n[], FTensor::Tensor1< double *, 3 > &t_phi, FTensor::Tensor2< double *, 3, 3 > &t_diff_phi)
Definition: Hcurl.cpp:2014
MoFEMErrorCode operator()(int shift, 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< double *, 3 > &t_phi, FTensor::Tensor2< FTensor::PackPtr< double *, 6 >, 3, 2 > &t_diff_phi)
Definition: Hcurl.cpp:2568
VectorDouble f0PhiII
Definition: Hcurl.cpp:2314
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
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:2717
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:359
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:2191
#define NBEDGE_AINSWORTH_HCURL(P)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
Managing BitRefLevels.
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:2278
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:659
Mesh refinement interface.
#define CHKERR
Inline error check.
Definition: definitions.h:595
#define NBFACETRI_DEMKOWICZ_HCURL(P)
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:1406
#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:15
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
VectorDouble iFiF1
Definition: Hcurl.cpp:2317
#define NBVOLUMETET_AINSWORTH_HCURL(P)
#define NBFACETRI_AINSWORTH_HCURL(P)
#define NBFACETRI_AINSWORTH_EDGE_HCURL(P)
VectorDouble iFiF0
Definition: Hcurl.cpp:2316
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:432
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:772
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:76
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.
const int N
Definition: speed_test.cpp:3
MoFEMErrorCode operator()(int shift, 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< double *, 3 > &t_phi, FTensor::Tensor2< double *, 3, 3 > &t_diff_phi)
Definition: Hcurl.cpp:2410
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:540
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:2237
#define NBVOLUMETET_DEMKOWICZ_HCURL(P)
int main(int argc, char *argv[])
const int order
approximation order
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:174