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