v0.15.0
Loading...
Searching...
No Matches
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
14using namespace MoFEM;
15
16MoFEMErrorCode 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
175MoFEMErrorCode 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
237MoFEMErrorCode 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
363MoFEMErrorCode 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
458MoFEMErrorCode 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
545MoFEMErrorCode 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 SETERRQ(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
663MoFEMErrorCode 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 SETERRQ(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
775MoFEMErrorCode 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 SETERRQ(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
909MoFEMErrorCode 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 SETERRQ(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
1052MoFEMErrorCode 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 SETERRQ(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 SETERRQ(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, "%s", 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, "%s", ss.str().c_str());
1244 }
1245
1246 MoFEMFunctionReturnHot(0);
1247}
1248
1249MoFEMErrorCode 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 SETERRQ(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 SETERRQ(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
1403MoFEMErrorCode 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 SETERRQ(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 SETERRQ(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>
1525using namespace MoFEM;
1526using namespace boost::numeric;
1527
1528MoFEMErrorCode 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
1978struct 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
2068MoFEMErrorCode 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
2105MoFEMErrorCode 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
2144MoFEMErrorCode 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
2176struct 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
2395MoFEMErrorCode 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
2433MoFEMErrorCode 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
2468MoFEMErrorCode 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
2615MoFEMErrorCode 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
2796static char help[] = "...\n\n";
2797
2798int 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
Hcurl.hpp
Implementation of H-curl base function.
help
static char help[]
Definition activate_deactivate_dofs.cpp:13
main
int main()
Definition adol-c_atom.cpp:46
a
constexpr double a
Definition approx_sphere.cpp:30
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
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
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition definitions.h:385
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:463
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359
HVEC0
@ HVEC0
Definition definitions.h:199
HVEC1
@ HVEC1
Definition definitions.h:199
HVEC2
@ HVEC2
Definition definitions.h:199
MOFEM_STD_EXCEPTION_THROW
@ MOFEM_STD_EXCEPTION_THROW
Definition definitions.h:39
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:432
HVEC1_1
@ HVEC1_1
Definition definitions.h:209
HVEC0_1
@ HVEC0_1
Definition definitions.h:208
HVEC1_0
@ HVEC1_0
Definition definitions.h:206
HVEC2_1
@ HVEC2_1
Definition definitions.h:210
HVEC2_0
@ HVEC2_0
Definition definitions.h:207
HVEC0_0
@ HVEC0_0
Definition definitions.h:205
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:456
order
constexpr int order
Definition dg_projection.cpp:18
ShapeDiffMBTET
PetscErrorCode ShapeDiffMBTET(double *diffN)
calculate derivatives of shape functions
Definition fem_tools.c:319
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
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
NBEDGE_DEMKOWICZ_HCURL
#define NBEDGE_DEMKOWICZ_HCURL(P)
Definition h1_hdiv_hcurl_l2.h:108
NBVOLUMETET_AINSWORTH_TET_HCURL
#define NBVOLUMETET_AINSWORTH_TET_HCURL(P)
Definition h1_hdiv_hcurl_l2.h:103
NBFACETRI_AINSWORTH_FACE_HCURL
#define NBFACETRI_AINSWORTH_FACE_HCURL(P)
Definition h1_hdiv_hcurl_l2.h:99
NBVOLUMETET_AINSWORTH_HCURL
#define NBVOLUMETET_AINSWORTH_HCURL(P)
Definition h1_hdiv_hcurl_l2.h:105
NBFACETRI_AINSWORTH_HCURL
#define NBFACETRI_AINSWORTH_HCURL(P)
Definition h1_hdiv_hcurl_l2.h:100
NBVOLUMETET_AINSWORTH_FACE_HCURL
#define NBVOLUMETET_AINSWORTH_FACE_HCURL(P)
Definition h1_hdiv_hcurl_l2.h:101
NBVOLUMETET_DEMKOWICZ_HCURL
#define NBVOLUMETET_DEMKOWICZ_HCURL(P)
Definition h1_hdiv_hcurl_l2.h:110
NBFACETRI_DEMKOWICZ_HCURL
#define NBFACETRI_DEMKOWICZ_HCURL(P)
Definition h1_hdiv_hcurl_l2.h:109
NBFACETRI_AINSWORTH_EDGE_HCURL
#define NBFACETRI_AINSWORTH_EDGE_HCURL(P)
Definition h1_hdiv_hcurl_l2.h:98
NBEDGE_AINSWORTH_HCURL
#define NBEDGE_AINSWORTH_HCURL(P)
Definition h1_hdiv_hcurl_l2.h:97
i
FTensor::Index< 'i', SPACE_DIM > i
Definition hcurl_divergence_operator_2d.cpp:27
n
const double n
refractive index of diffusive medium
Definition initial_diffusion.cpp:38
j
FTensor::Index< 'j', 3 > j
Definition matrix_function.cpp:19
k
FTensor::Index< 'k', 3 > k
Definition matrix_function.cpp:20
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition Exceptions.hpp:56
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
MoFEM
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
MoFEM::getFTensor2HVecFromPtr< 3, 3 >
FTensor::Tensor2< FTensor::PackPtr< double *, 9 >, 3, 3 > getFTensor2HVecFromPtr< 3, 3 >(double *ptr)
Definition Templates.hpp:968
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
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::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_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::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
MoFEM::getFTensor2FromPtr
auto getFTensor2FromPtr(double *ptr)
Make Tensor2 from pointer.
Definition Templates.hpp:925
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
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
MoFEM::getFTensor2HVecFromPtr< 3, 2 >
FTensor::Tensor2< FTensor::PackPtr< double *, 6 >, 3, 2 > getFTensor2HVecFromPtr< 3, 2 >(double *ptr)
Definition Templates.hpp:957
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
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_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::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
MoFEM::getFTensor1FromPtr
FTensor::Tensor1< FTensor::PackPtr< double *, S >, DIM > getFTensor1FromPtr(double *ptr)
Make Tensor1 from pointer.
Definition Templates.hpp:838
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
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::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
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::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
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
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
t
constexpr double t
plate stiffness
Definition plate.cpp:58
phi
static double phi
Definition poisson_2d_dis_galerkin.cpp:30
N
const int N
Definition speed_test.cpp:3
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::tDiffb
FTensor::Tensor1< double, 3 > tDiffb
Definition Hcurl.cpp:1984
HcurlEdgeBase::i
FTensor::Index< 'i', 3 > i
Definition Hcurl.cpp:1980
HcurlEdgeBase::diffFi
MatrixDouble diffFi
Definition Hcurl.cpp:1987
HcurlEdgeBase::tGradN0pN1
FTensor::Tensor1< double, 3 > tGradN0pN1
Definition Hcurl.cpp:1982
HcurlEdgeBase::fI
VectorDouble fI
Definition Hcurl.cpp:1986
HcurlEdgeBase::tPhi0
FTensor::Tensor1< double, 3 > tPhi0
Definition Hcurl.cpp:1983
HcurlFaceBase::iFiF0
VectorDouble iFiF0
Definition Hcurl.cpp:2180
HcurlFaceBase::diffF0PhiII
VectorDouble diffF0PhiII
Definition Hcurl.cpp:2178
HcurlFaceBase::diffIFiF1
VectorDouble diffIFiF1
Definition Hcurl.cpp:2181
HcurlFaceBase::hCurlBaseOnEdge
HcurlEdgeBase hCurlBaseOnEdge
Definition Hcurl.cpp:2177
HcurlFaceBase::iFiF1
VectorDouble iFiF1
Definition Hcurl.cpp:2181
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
HcurlFaceBase::f1PhiII
VectorDouble f1PhiII
Definition Hcurl.cpp:2179
HcurlFaceBase::i
FTensor::Index< 'i', 3 > i
Definition Hcurl.cpp:2183
HcurlFaceBase::diffIFiF0
VectorDouble diffIFiF0
Definition Hcurl.cpp:2180
HcurlFaceBase::diffF1PhiII
VectorDouble diffF1PhiII
Definition Hcurl.cpp:2179
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
HcurlFaceBase::f0PhiII
VectorDouble f0PhiII
Definition Hcurl.cpp:2178
MoFEM::BitRefManager
Managing BitRefLevels.
Definition BitRefManager.hpp:21
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
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:118
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition DeprecatedCoreInterface.hpp:16
MoFEM::MeshRefinement
Mesh refinement interface.
Definition MeshRefinement.hpp:26
MoFEM::MeshRefinement::refineTets
MoFEMErrorCode refineTets(const EntityHandle meshset, const BitRefLevel &bit, int verb=QUIET, const bool debug=false)
refine TET in the meshset
Definition MeshRefinement.cpp:197
MoFEM::MeshRefinement::addVerticesInTheMiddleOfEdges
MoFEMErrorCode addVerticesInTheMiddleOfEdges(const EntityHandle meshset, const BitRefLevel &bit, const bool recursive=false, int verb=QUIET, EntityHandle start_v=0)
make vertices in the middle of edges in meshset and add them to refinement levels defined by bit
Definition MeshRefinement.cpp:42
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition UnknownInterface.hpp:93
Generated by Doxygen 1.12.0 and hosted at University of Glasgow