v0.15.0
Loading...
Searching...
No Matches
Hdiv.cpp
Go to the documentation of this file.
1/** \file Hdiv.cpp
2
3 \brief Implementation of H-curl base
4
5 Based on Hierarchic Finite Element Bases on Unstructured Tetrahedral
6 Meshes, by Mark Ainsworth and Joe Coyle
7 Shape functions for MBTRI/MBTET and HCurl space
8
9*/
10
11using namespace MoFEM;
12
13boost::function<int(int)> AinsworthOrderHooks::broken_nbfacetri_edge_hdiv =
14 [](int p) { return p; };
15boost::function<int(int)> AinsworthOrderHooks::broken_nbfacetri_face_hdiv =
16 [](int p) { return p; };
17boost::function<int(int)> AinsworthOrderHooks::broken_nbvolumetet_edge_hdiv =
18 [](int p) { return p; };
19boost::function<int(int)> AinsworthOrderHooks::broken_nbvolumetet_face_hdiv =
20 [](int p) { return p; };
21boost::function<int(int)> AinsworthOrderHooks::broken_nbvolumetet_volume_hdiv =
22 [](int p) { return p; };
23
24MoFEMErrorCode MoFEM::Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET(
25 int *faces_nodes, int *p, double *N, double *diffN, double *phi_f_e[4][3],
26 double *diff_phi_f_e[4][3], int gdim,
27 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
28 double *L, double *diffL,
29 const int dim)) {
30
31 MoFEMFunctionBeginHot;
32#ifndef NDEBUG
33 if (!diff_phi_f_e)
34 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
35 "expected to return derivatives");
36#endif
37
38 for (int ff = 0; ff < 4; ff++) {
39 CHKERR Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET_ON_FACE(
40 &faces_nodes[3 * ff], p[ff], N, diffN, phi_f_e[ff], diff_phi_f_e[ff],
41 gdim, 4, base_polynomials);
42 }
43
44 MoFEMFunctionReturnHot(0);
45}
46
47MoFEMErrorCode MoFEM::Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET_ON_FACE(
48 int *faces_nodes, int p, double *N, double *diffN, double *phi_f_e[3],
49 double *diff_phi_f_e[3], int gdim, int nb,
50 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
51 double *L, double *diffL,
52 const int dim)) {
53
54 constexpr int face_edges_nodes[3][2] = {{0, 1}, {1, 2}, {2, 0}};
55 constexpr int face_opposite_edges_node[] = {2, 0, 1};
56 FTENSOR_INDEX(3, i);
57 FTENSOR_INDEX(3, j);
58 FTENSOR_INDEX(3, k);
59
60 MoFEMFunctionBeginHot;
61 if (p < 1)
62 MoFEMFunctionReturnHot(0);
63
64 FTensor::Tensor1<double, 3> t_edge_cross[3];
65 FTensor::Tensor1<double, 3> t_node_diff_ksi[4];
66 FTensor::Tensor1<double, 3> t_diff_ksi[3];
67 if (diffN) {
68 t_node_diff_ksi[0] =
69 FTensor::Tensor1<double, 3>(diffN[0], diffN[1], diffN[2]);
70 t_node_diff_ksi[1] =
71 FTensor::Tensor1<double, 3>(diffN[3], diffN[4], diffN[5]);
72 t_node_diff_ksi[2] =
73 FTensor::Tensor1<double, 3>(diffN[6], diffN[7], diffN[8]);
74 t_node_diff_ksi[3] =
75 FTensor::Tensor1<double, 3>(diffN[9], diffN[10], diffN[11]);
76 for (int ee = 0; ee < 3; ee++) {
77 const int n0 = faces_nodes[face_edges_nodes[ee][0]];
78 const int n1 = faces_nodes[face_edges_nodes[ee][1]];
79 t_diff_ksi[ee](i) = t_node_diff_ksi[n1](i) - t_node_diff_ksi[n0](i);
80 t_edge_cross[ee](i) = levi_civita(i, j, k) * t_node_diff_ksi[n0](j) *
81 t_node_diff_ksi[n1](k);
82 }
83 } else {
84 for (int ee = 0; ee < 3; ee++) {
85 t_edge_cross[ee](0) = 1;
86 t_edge_cross[ee](1) = 0;
87 t_edge_cross[ee](2) = 0;
88 }
89 }
90 double psi_l[p + 1];
91
92 for (int ee = 0; ee != 3; ee++) {
93 const int i0 = faces_nodes[face_edges_nodes[ee][0]];
94 const int i1 = faces_nodes[face_edges_nodes[ee][1]];
95 const int iO = faces_nodes[face_opposite_edges_node[ee]];
96 auto t_psi_f_e = getFTensor1FromPtr<3>(phi_f_e[ee]);
97 for (int ii = 0; ii != nb * gdim; ii += nb) {
98 const double n0 = N[ii + i0];
99 const double n1 = N[ii + i1];
100 const double lambda = N[ii + iO];
101 const double ksi = n1 - n0;
102 base_polynomials(p, ksi, NULL, psi_l, NULL, 3);
103 auto t_psi_l = FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(psi_l);
104 int jj = 0;
105 for (int l = 0; l <= p - 1; l++) {
106 t_psi_f_e(i) = (lambda * t_psi_l) * t_edge_cross[ee](i);
107 ++t_psi_f_e;
108 ++t_psi_l;
109 ++jj;
110 }
111#ifndef NDEBUG
112 if (jj != NBFACETRI_AINSWORTH_EDGE_HDIV(p)) {
113 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
114 "wrong order %d != %d", jj, NBFACETRI_AINSWORTH_EDGE_HDIV(p));
115 }
116#endif
117 }
118 }
119
120 if (diff_phi_f_e) {
121 double diff_psi_l[3 * (p + 1)];
122 for (int ee = 0; ee != 3; ee++) {
123 const int i0 = faces_nodes[face_edges_nodes[ee][0]];
124 const int i1 = faces_nodes[face_edges_nodes[ee][1]];
125 const int iO = faces_nodes[face_opposite_edges_node[ee]];
126 auto t_diff_phi_f_e = getFTensor2HVecFromPtr<3, 3>(diff_phi_f_e[ee]);
127 for (int ii = 0; ii != nb * gdim; ii += nb) {
128 const double n0 = N[ii + i0];
129 const double n1 = N[ii + i1];
130 const double lambda = N[ii + iO];
131 const double ksi = n1 - n0;
132 base_polynomials(p, ksi, &t_diff_ksi[ee](0), psi_l, diff_psi_l, 3);
133 auto t_psi_l = FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(psi_l);
134 auto t_diff_psi_l = FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3>{
135 &diff_psi_l[0], &diff_psi_l[p + 1], &diff_psi_l[2 * p + 2]};
136 for (int l = 0; l <= p - 1; l++) {
137 t_diff_phi_f_e(i, j) =
138 (t_node_diff_ksi[iO](j) * t_psi_l + lambda * t_diff_psi_l(j)) *
139 t_edge_cross[ee](i);
140 ++t_diff_psi_l;
141 ++t_diff_phi_f_e;
142 ++t_psi_l;
143 }
144 }
145 }
146 }
147
148 MoFEMFunctionReturnHot(0);
149}
150
151MoFEMErrorCode MoFEM::Hdiv_Ainsworth_FaceBubbleShapeFunctions(
152 int *faces_nodes, int *p, double *N, double *diffN, double *phi_f[],
153 double *diff_phi_f[], int gdim,
154 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
155 double *L, double *diffL,
156 const int dim)) {
157
158 MoFEMFunctionBeginHot;
159
160#ifndef NDEBUG
161 if (!diff_phi_f)
162 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
163 "expected to return derivatives");
164#endif
165
166 for (int ff = 0; ff < 4; ff++) {
167 CHKERR Hdiv_Ainsworth_FaceBubbleShapeFunctions_ON_FACE(
168 &faces_nodes[3 * ff], p[ff], N, diffN, phi_f[ff], diff_phi_f[ff], gdim,
169 4, base_polynomials);
170 }
171 MoFEMFunctionReturnHot(0);
172}
173
174MoFEMErrorCode MoFEM::Hdiv_Ainsworth_FaceBubbleShapeFunctions_ON_FACE(
175 int *face_nodes, int p, double *N, double *diffN, double *phi_f,
176 double *diff_phi_f, int gdim, int nb,
177 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
178 double *L, double *diffL,
179 const int dim)) {
180 FTENSOR_INDEX(3, i);
181 FTENSOR_INDEX(3, j);
182 FTENSOR_INDEX(3, k);
183
184 MoFEMFunctionBeginHot;
185 if (p < 3)
186 MoFEMFunctionReturnHot(0);
187
188 const int vert_i = face_nodes[1];
189 const int vert_j = face_nodes[2];
190 const int i0 = face_nodes[0];
191 FTensor::Tensor1<double, 3> t_cross;
192 FTensor::Tensor1<double, 3> t_node_diff_ksi[4];
193 FTensor::Tensor1<double, 3> t_diff_ksi0i;
194 FTensor::Tensor1<double, 3> t_diff_ksi0j;
195
196 if (diffN) {
197 t_node_diff_ksi[0] =
198 FTensor::Tensor1<double, 3>(diffN[0], diffN[1], diffN[2]);
199 t_node_diff_ksi[1] =
200 FTensor::Tensor1<double, 3>(diffN[3], diffN[4], diffN[5]);
201 t_node_diff_ksi[2] =
202 FTensor::Tensor1<double, 3>(diffN[6], diffN[7], diffN[8]);
203 t_node_diff_ksi[3] =
204 FTensor::Tensor1<double, 3>(diffN[9], diffN[10], diffN[11]);
205 t_diff_ksi0i(i) = t_node_diff_ksi[vert_i](i) - t_node_diff_ksi[i0](i);
206 t_diff_ksi0j(i) = t_node_diff_ksi[vert_j](i) - t_node_diff_ksi[i0](i);
207 t_cross(i) = levi_civita(i, j, k) * t_node_diff_ksi[vert_i](j) *
208 t_node_diff_ksi[vert_j](k);
209 } else {
210 t_cross(0) = 1;
211 t_cross(1) = 0;
212 t_cross(2) = 0;
213 }
214
215 double psi_l[p + 1], diff_psi_l[3 * (p + 1)];
216 double psi_m[p + 1], diff_psi_m[3 * (p + 1)];
217 auto t_psi_f = getFTensor1FromPtr<3>(phi_f);
218
219 for (int ii = 0; ii < gdim; ii++) {
220
221 const int node_shift = ii * nb;
222 const double ni = N[node_shift + vert_i];
223 const double nj = N[node_shift + vert_j];
224 const double n0 = N[node_shift + i0];
225 const double ksi0i = ni - n0;
226 const double ksi0j = nj - n0;
227 double beta_0ij = n0 * ni * nj;
228 base_polynomials(p, ksi0i, NULL, psi_l, NULL, 3);
229 base_polynomials(p, ksi0j, NULL, psi_m, NULL, 3);
230
231 int jj = 0;
232 int oo = 0;
233 for (; oo <= p - 3; oo++) {
234 auto t_psi_l = FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(psi_l);
235 for (int l = 0; l <= oo; l++) {
236 int m = oo - l;
237 if (m >= 0) {
238 t_psi_f(i) = (beta_0ij * t_psi_l * psi_m[m]) * t_cross(i);
239 ++t_psi_f;
240 }
241 ++t_psi_l;
242 ++jj;
243 }
244 }
245#ifndef NDEBUG
246 if (jj != NBFACETRI_AINSWORTH_FACE_HDIV(p)) {
247 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
248 "wrong order %d != %d", jj, NBFACETRI_AINSWORTH_FACE_HDIV(p));
249 }
250#endif
251 }
252
253 if (diff_phi_f) {
254 auto t_diff_phi_f = getFTensor2HVecFromPtr<3, 3>(diff_phi_f);
255
256 for (int ii = 0; ii < gdim; ii++) {
257
258 const int node_shift = ii * nb;
259 const double ni = N[node_shift + vert_i];
260 const double nj = N[node_shift + vert_j];
261 const double n0 = N[node_shift + i0];
262 const double ksi0i = ni - n0;
263 const double ksi0j = nj - n0;
264 double beta_0ij = n0 * ni * nj;
265 FTensor::Tensor1<double, 3> t_diff_beta_0ij;
266 t_diff_beta_0ij(i) = (ni * nj) * t_node_diff_ksi[i0](i) +
267 (n0 * nj) * t_node_diff_ksi[vert_i](i) +
268 (n0 * ni) * t_node_diff_ksi[vert_j](i);
269 base_polynomials(p, ksi0i, &t_diff_ksi0i(0), psi_l, diff_psi_l, 3);
270 base_polynomials(p, ksi0j, &t_diff_ksi0j(0), psi_m, diff_psi_m, 3);
271
272 int jj = 0;
273 int oo = 0;
274 for (; oo <= p - 3; oo++) {
275 auto t_psi_l = FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(psi_l);
276 auto t_diff_psi_l = FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3>{
277 &diff_psi_l[0], &diff_psi_l[p + 1], &diff_psi_l[2 * p + 2]};
278 for (int l = 0; l <= oo; l++) {
279 int m = oo - l;
280 if (m >= 0) {
281 auto t_diff_psi_m =
282 FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3>{
283 &diff_psi_m[m], &diff_psi_m[p + 1 + m],
284 &diff_psi_m[2 * p + 2 + m]};
285 t_diff_phi_f(i, j) = ((t_psi_l * psi_m[m]) * t_diff_beta_0ij(j) +
286 (beta_0ij * psi_m[m]) * t_diff_psi_l(j) +
287 (beta_0ij * t_psi_l) * t_diff_psi_m(j)) *
288 t_cross(i);
289 ++t_diff_phi_f;
290 }
291 ++t_psi_l;
292 ++t_diff_psi_l;
293 ++jj;
294 }
295 }
296#ifndef NDEBUG
297 if (jj != NBFACETRI_AINSWORTH_FACE_HDIV(p)) {
298 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
299 "wrong order %d != %d", jj, NBFACETRI_AINSWORTH_FACE_HDIV(p));
300 }
301#endif
302 }
303 }
304 MoFEMFunctionReturnHot(0);
305}
306
307MoFEMErrorCode MoFEM::Hdiv_Ainsworth_EdgeBasedVolumeShapeFunctions_MBTET(
308 int p, double *N, double *diffN, double *phi_v_e[6],
309 double *diff_phi_v_e[6], int gdim,
310 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
311 double *L, double *diffL,
312 const int dim)) {
313
314 constexpr int edges_nodes[6][2] = {{0, 1}, {1, 2}, {2, 0},
315 {0, 3}, {1, 3}, {2, 3}};
316
317 MoFEMFunctionBeginHot;
318 if (p < 2)
319 MoFEMFunctionReturnHot(0);
320 if (diffN == NULL) {
321 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
322 }
323
324 FTensor::Tensor1<double, 3> t_coords[4] = {
325 FTensor::Tensor1<double, 3>(0., 0., 0.),
326 FTensor::Tensor1<double, 3>(1., 0., 0.),
327 FTensor::Tensor1<double, 3>(0., 1., 0.),
328 FTensor::Tensor1<double, 3>(0., 0., 1.)};
329 FTensor::Tensor1<double *, 3> t_node_diff_ksi[4] = {
330 FTensor::Tensor1<double *, 3>(&diffN[0], &diffN[1], &diffN[2]),
331 FTensor::Tensor1<double *, 3>(&diffN[3], &diffN[4], &diffN[5]),
332 FTensor::Tensor1<double *, 3>(&diffN[6], &diffN[7], &diffN[8]),
333 FTensor::Tensor1<double *, 3>(&diffN[9], &diffN[10], &diffN[11])};
334
335 FTENSOR_INDEX(3, i);
336 FTENSOR_INDEX(3, j);
337
338 FTensor::Tensor1<double, 3> t_tou_e;
339 FTensor::Tensor1<double, 3> t_diff_ksi0i;
340 FTensor::Tensor1<double, 3> t_diff_beta_e;
341
342 double psi_l[p + 1];
343 double diff_psi_l[3 * (p + 1)];
344
345 for (int ee = 0; ee != 6; ee++) {
346 t_tou_e(i) =
347 t_coords[edges_nodes[ee][1]](i) - t_coords[edges_nodes[ee][0]](i);
348 t_diff_ksi0i(i) = t_node_diff_ksi[edges_nodes[ee][1]](i) -
349 t_node_diff_ksi[edges_nodes[ee][0]](i);
350 auto t_phi_v_e = getFTensor1FromPtr<3>(phi_v_e[ee]);
351 for (int ii = 0; ii != gdim; ii++) {
352 const int node_shift = ii * 4;
353 const double ni = N[node_shift + edges_nodes[ee][1]];
354 const double n0 = N[node_shift + edges_nodes[ee][0]];
355 const double beta_e = ni * n0;
356 const double ksi0i = ni - n0;
357 base_polynomials(p, ksi0i, NULL, psi_l, NULL, 3);
358 auto t_psi_l = FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(psi_l);
359 for (int l = 0; l <= p - 2; l++) {
360 t_phi_v_e(i) = (beta_e * t_psi_l) * t_tou_e(i);
361 ++t_phi_v_e;
362 ++t_psi_l;
363 }
364 }
365 }
366
367 if (diff_phi_v_e) {
368 for (int ee = 0; ee != 6; ee++) {
369 t_tou_e(i) =
370 t_coords[edges_nodes[ee][1]](i) - t_coords[edges_nodes[ee][0]](i);
371 t_diff_ksi0i(i) = t_node_diff_ksi[edges_nodes[ee][1]](i) -
372 t_node_diff_ksi[edges_nodes[ee][0]](i);
373 auto t_diff_phi_v_e = getFTensor2HVecFromPtr<3, 3>(diff_phi_v_e[ee]);
374 for (int ii = 0; ii != gdim; ii++) {
375 const int node_shift = ii * 4;
376 const double ni = N[node_shift + edges_nodes[ee][1]];
377 const double n0 = N[node_shift + edges_nodes[ee][0]];
378 const double beta_e = ni * n0;
379 const double ksi0i = ni - n0;
380 t_diff_beta_e(i) = ni * t_node_diff_ksi[edges_nodes[ee][0]](i) +
381 t_node_diff_ksi[edges_nodes[ee][1]](i) * n0;
382 base_polynomials(p, ksi0i, &t_diff_ksi0i(0), psi_l, diff_psi_l, 3);
383 auto t_psi_l = FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(psi_l);
384 auto t_diff_psi_l = FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3>(
385 diff_psi_l, &diff_psi_l[p + 1], &diff_psi_l[2 * p + 2]);
386 for (int l = 0; l <= p - 2; l++) {
387 t_diff_phi_v_e(i, j) =
388 (t_diff_beta_e(j) * t_psi_l + beta_e * t_diff_psi_l(j)) *
389 t_tou_e(i);
390 ++t_diff_phi_v_e;
391 ++t_diff_psi_l;
392 ++t_psi_l;
393 }
394 }
395 }
396 }
397
398 MoFEMFunctionReturnHot(0);
399}
400
401MoFEMErrorCode MoFEM::Hdiv_Ainsworth_FaceBasedVolumeShapeFunctions_MBTET(
402 int p, double *N, double *diffN, double *phi_v_f[], double *diff_phi_v_f[],
403 int gdim,
404 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
405 double *L, double *diffL,
406 const int dim)) {
407
408 constexpr int faces_nodes[4][3] = {
409 {0, 1, 3}, {1, 2, 3}, {0, 3, 2}, {0, 2, 1}};
410
411 MoFEMFunctionBeginHot;
412 if (p < 3)
413 MoFEMFunctionReturnHot(0);
414
415 FTensor::Tensor1<double, 3> t_coords[4] = {
416 FTensor::Tensor1<double, 3>(0., 0., 0.),
417 FTensor::Tensor1<double, 3>(1., 0., 0.),
418 FTensor::Tensor1<double, 3>(0., 1., 0.),
419 FTensor::Tensor1<double, 3>(0., 0., 1.)};
420
421 FTensor::Tensor1<double *, 3> t_node_diff_ksi[4] = {
422 FTensor::Tensor1<double *, 3>(&diffN[0], &diffN[1], &diffN[2]),
423 FTensor::Tensor1<double *, 3>(&diffN[3], &diffN[4], &diffN[5]),
424 FTensor::Tensor1<double *, 3>(&diffN[6], &diffN[7], &diffN[8]),
425 FTensor::Tensor1<double *, 3>(&diffN[9], &diffN[10], &diffN[11])};
426
427 FTENSOR_INDEX(3, i);
428 FTENSOR_INDEX(3, j);
429
430 FTensor::Tensor1<double, 3> t_tau0i[4], t_tau0j[4];
431 FTensor::Tensor1<double, 3> t_diff_ksi0i[4], t_diff_ksi0j[4];
432 for (int ff = 0; ff != 4; ff++) {
433 const int v0 = faces_nodes[ff][0];
434 const int vi = faces_nodes[ff][1];
435 const int vj = faces_nodes[ff][2];
436 t_tau0i[ff](i) = t_coords[vi](i) - t_coords[v0](i);
437 t_tau0j[ff](i) = t_coords[vj](i) - t_coords[v0](i);
438 t_diff_ksi0i[ff](i) = t_node_diff_ksi[vi](i) - t_node_diff_ksi[v0](i);
439 t_diff_ksi0j[ff](i) = t_node_diff_ksi[vj](i) - t_node_diff_ksi[v0](i);
440 }
441
442 double psi_l[p + 1], psi_m[p + 1];
443 double diff_psi_l[3 * (p + 1)], diff_psi_m[3 * (p + 1)];
444 for (int ff = 0; ff != 4; ff++) {
445 const int v0 = faces_nodes[ff][0];
446 const int vi = faces_nodes[ff][1];
447 const int vj = faces_nodes[ff][2];
448 auto t_phi_v_f = getFTensor1FromPtr<3>(phi_v_f[ff]);
449 auto t_diff_phi_v_f = getFTensor2HVecFromPtr<3, 3>(diff_phi_v_f[ff]);
450 for (int ii = 0; ii < gdim; ii++) {
451 const int node_shift = 4 * ii;
452 const double n0 = N[node_shift + v0];
453 const double ni = N[node_shift + vi];
454 const double nj = N[node_shift + vj];
455 const double beta_f = n0 * ni * nj;
456 FTensor::Tensor1<double, 3> t_diff_beta_f;
457 t_diff_beta_f(i) = (ni * nj) * t_node_diff_ksi[v0](i) +
458 (n0 * nj) * t_node_diff_ksi[vi](i) +
459 (n0 * ni) * t_node_diff_ksi[vj](i);
460 const double ksi0i = ni - n0;
461 const double ksi0j = nj - n0;
462 base_polynomials(p, ksi0i, &t_diff_ksi0i[ff](0), psi_l, diff_psi_l, 3);
463 base_polynomials(p, ksi0j, &t_diff_ksi0j[ff](0), psi_m, diff_psi_m, 3);
464 FTensor::Tensor1<double, 3> t_diff_a;
465 int jj = 0;
466 for (int oo = 0; oo <= p - 3; oo++) {
467 auto t_psi_l = FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(psi_l);
468 auto t_diff_psi_l = FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3>(
469 diff_psi_l, &diff_psi_l[p + 1], &diff_psi_l[2 * p + 2]);
470 for (int l = 0; l <= oo; l++) {
471 int m = oo - l;
472 if (m >= 0) {
473 auto t_diff_psi_m = FTensor::Tensor1<double, 3>{
474 diff_psi_m[m], diff_psi_m[p + 1 + m],
475 diff_psi_m[2 * p + 2 + m]};
476 const double a = beta_f * t_psi_l * psi_m[m];
477 t_phi_v_f(i) = a * t_tau0i[ff](i);
478 ++t_phi_v_f;
479 ++jj;
480 t_phi_v_f(i) = a * t_tau0j[ff](i);
481 ++t_phi_v_f;
482 ++jj;
483
484 t_diff_a(j) = (t_psi_l * psi_m[m]) * t_diff_beta_f(j) +
485 (beta_f * psi_m[m]) * t_diff_psi_l(j) +
486 (beta_f * t_psi_l) * t_diff_psi_m(j);
487 t_diff_phi_v_f(i, j) = t_diff_a(j) * t_tau0i[ff](i);
488 ++t_diff_phi_v_f;
489 t_diff_phi_v_f(i, j) = t_diff_a(j) * t_tau0j[ff](i);
490 ++t_diff_phi_v_f;
491
492 ++t_psi_l;
493 ++t_diff_psi_l;
494 }
495 }
496 }
497 if (jj != NBVOLUMETET_AINSWORTH_FACE_HDIV(p)) {
498 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
499 "wrong order %d != %d", jj,
500 NBVOLUMETET_AINSWORTH_FACE_HDIV(p));
501 }
502 }
503 }
504 MoFEMFunctionReturnHot(0);
505}
506
507MoFEMErrorCode MoFEM::Hdiv_Ainsworth_VolumeBubbleShapeFunctions_MBTET(
508 int p, double *N, double *diffN, double *phi_v, double *diff_phi_v,
509 int gdim,
510 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
511 double *L, double *diffL,
512 const int dim)) {
513
514 MoFEMFunctionBeginHot;
515 if (p < 4)
516 MoFEMFunctionReturnHot(0);
517
518 FTensor::Tensor1<double *, 3> t_node_diff_ksi[4] = {
519 FTensor::Tensor1<double *, 3>(&diffN[0], &diffN[1], &diffN[2]),
520 FTensor::Tensor1<double *, 3>(&diffN[3], &diffN[4], &diffN[5]),
521 FTensor::Tensor1<double *, 3>(&diffN[6], &diffN[7], &diffN[8]),
522 FTensor::Tensor1<double *, 3>(&diffN[9], &diffN[10], &diffN[11])};
523
524 FTensor::Index<'i', 3> i;
525 FTensor::Index<'j', 3> j;
526 FTensor::Number<0> N0;
527 FTensor::Number<1> N1;
528 FTensor::Number<2> N2;
529
530 FTensor::Tensor1<double, 3> t_diff_ksi0i;
531 FTensor::Tensor1<double, 3> t_diff_ksi0j;
532 FTensor::Tensor1<double, 3> t_diff_ksi0k;
533
534 t_diff_ksi0i(i) = t_node_diff_ksi[1](i) - t_node_diff_ksi[0](i);
535 t_diff_ksi0j(i) = t_node_diff_ksi[2](i) - t_node_diff_ksi[0](i);
536 t_diff_ksi0k(i) = t_node_diff_ksi[3](i) - t_node_diff_ksi[0](i);
537
538 double psi_l[p + 1];
539 double diff_psi_l[3 * (p + 1)];
540 double psi_m[p + 1];
541 double diff_psi_m[3 * (p + 1)];
542 double psi_n[p + 1];
543 double diff_psi_n[3 * (p + 1)];
544
545 auto t_phi_v = getFTensor1FromPtr<3>(phi_v);
546 auto t_diff_phi_v = getFTensor2HVecFromPtr<3, 3>(diff_phi_v);
547
548 FTensor::Tensor1<double, 3> t_diff_beta_v;
549 for (int ii = 0; ii < gdim; ii++) {
550 const int node_shift = ii * 4;
551 const double n0 = N[node_shift + 0];
552 const double ni = N[node_shift + 1];
553 const double nj = N[node_shift + 2];
554 const double nk = N[node_shift + 3];
555 const double ksi0i = ni - n0;
556 const double ksi0j = nj - n0;
557 const double ksi0k = nk - n0;
558 const double beta_v = n0 * ni * nj * nk;
559 t_diff_beta_v(i) = (ni * nj * nk) * t_node_diff_ksi[0](i) +
560 (n0 * nj * nk) * t_node_diff_ksi[1](i) +
561 (n0 * ni * nk) * t_node_diff_ksi[2](i) +
562 (n0 * ni * nj) * t_node_diff_ksi[3](i);
563 base_polynomials(p, ksi0i, &t_diff_ksi0i(0), psi_l, diff_psi_l, 3);
564 base_polynomials(p, ksi0j, &t_diff_ksi0j(0), psi_m, diff_psi_m, 3);
565 base_polynomials(p, ksi0k, &t_diff_ksi0k(0), psi_n, diff_psi_n, 3);
566
567 FTensor::Tensor1<double, 3> t_diff_a;
568
569 int jj = 0;
570 for (int oo = 0; oo <= p - 4; oo++) {
571 auto t_psi_l = FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(psi_l);
572 auto t_diff_psi_l = FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3>(
573 diff_psi_l, &diff_psi_l[p + 1], &diff_psi_l[2 * p + 2]);
574 for (int l = 0; l <= oo; l++) {
575 auto t_psi_m = FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(psi_m);
576 auto t_diff_psi_m = FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3>(
577 diff_psi_m, &diff_psi_m[p + 1], &diff_psi_m[2 * p + 2]);
578 for (int m = 0; (l + m) <= oo; m++) {
579 int n = oo - l - m;
580 if (n >= 0) {
581 FTensor::Tensor1<double, 3> t_diff_psi_n(diff_psi_n[n],
582 diff_psi_n[p + 1 + n],
583 diff_psi_n[2 * p + 2 + n]);
584 const double a = beta_v * t_psi_l * t_psi_m * psi_n[n];
585 t_phi_v(0) = a;
586 t_phi_v(1) = 0;
587 t_phi_v(2) = 0;
588 ++t_phi_v;
589 t_phi_v(0) = 0;
590 t_phi_v(1) = a;
591 t_phi_v(2) = 0;
592 ++t_phi_v;
593 t_phi_v(0) = 0;
594 t_phi_v(1) = 0;
595 t_phi_v(2) = a;
596 ++t_phi_v;
597 t_diff_a(j) = (t_psi_l * t_psi_m * psi_n[n]) * t_diff_beta_v(j) +
598 (beta_v * t_psi_m * psi_n[n]) * t_diff_psi_l(j) +
599 (beta_v * t_psi_l * psi_n[n]) * t_diff_psi_m(j) +
600 (beta_v * t_psi_l * t_psi_m) * t_diff_psi_n(j);
601 t_diff_phi_v(N0, j) = t_diff_a(j);
602 t_diff_phi_v(N1, j) = 0;
603 t_diff_phi_v(N2, j) = 0;
604 ++t_diff_phi_v;
605 t_diff_phi_v(N0, j) = 0;
606 t_diff_phi_v(N1, j) = t_diff_a(j);
607 t_diff_phi_v(N2, j) = 0;
608 ++t_diff_phi_v;
609 t_diff_phi_v(N0, j) = 0;
610 t_diff_phi_v(N1, j) = 0;
611 t_diff_phi_v(N2, j) = t_diff_a(j);
612 ++t_diff_phi_v;
613 ++jj;
614 }
615 ++t_psi_m;
616 ++t_diff_psi_m;
617 }
618 ++t_psi_l;
619 ++t_diff_psi_l;
620 }
621 }
622
623 if (3 * jj != NBVOLUMETET_AINSWORTH_VOLUME_HDIV(p)) {
624 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
625 "wrong order %d != %d", jj,
626 NBVOLUMETET_AINSWORTH_VOLUME_HDIV(p));
627 }
628 }
629
630 MoFEMFunctionReturnHot(0);
631}
632
633MoFEMErrorCode
634MoFEM::Hdiv_Demkowicz_Face_MBTET_ON_FACE(int *faces_nodes, int p, double *N,
635 double *diffN, double *phi_f,
636 double *diff_phi_f, int gdim, int nb) {
637 const int face_edges_nodes[3][2] = {{0, 1}, {1, 2}, {2, 0}};
638 const int face_opposite_edges_node[] = {2, 0, 1};
639
640 MoFEMFunctionBeginHot;
641
642 FTensor::Index<'i', 3> i;
643 FTensor::Index<'j', 3> j;
644
645 FTensor::Tensor1<double, 3> t_cross[3];
646 FTensor::Tensor2<double, 3, 3> t_diff_cross(0., 0., 0., 0., 0., 0., 0., 0.,
647 0.);
648 FTensor::Tensor1<double, 3> t_node_diff_ksi[4];
649 FTensor::Tensor1<double, 3> t_node_diff_sum_n0_n1;
650
651 const int i0 = faces_nodes[0];
652 const int i1 = faces_nodes[1];
653 const int i2 = faces_nodes[2];
654 const int o[] = {faces_nodes[face_opposite_edges_node[0]],
655 faces_nodes[face_opposite_edges_node[1]],
656 faces_nodes[face_opposite_edges_node[2]]};
657
658 FTensor::Tensor1<double, 3> t_diff_n0_p_n1;
659 FTensor::Tensor1<double, 3> t_diff_n0_p_n1_p_n2;
660
661 if (diff_phi_f) {
662 t_node_diff_ksi[0] =
663 FTensor::Tensor1<double, 3>(diffN[0], diffN[1], diffN[2]);
664 t_node_diff_ksi[1] =
665 FTensor::Tensor1<double, 3>(diffN[3], diffN[4], diffN[5]);
666 t_node_diff_ksi[2] =
667 FTensor::Tensor1<double, 3>(diffN[6], diffN[7], diffN[8]);
668 t_node_diff_ksi[3] =
669 FTensor::Tensor1<double, 3>(diffN[9], diffN[10], diffN[11]);
670 t_diff_cross(i, j) = 0;
671 for (int ee = 0; ee != 3; ee++) {
672 int ei0 = faces_nodes[face_edges_nodes[ee][0]];
673 int ei1 = faces_nodes[face_edges_nodes[ee][1]];
674 t_cross[ee](0) = t_node_diff_ksi[ei0](1) * t_node_diff_ksi[ei1](2) -
675 t_node_diff_ksi[ei0](2) * t_node_diff_ksi[ei1](1);
676 t_cross[ee](1) = t_node_diff_ksi[ei0](2) * t_node_diff_ksi[ei1](0) -
677 t_node_diff_ksi[ei0](0) * t_node_diff_ksi[ei1](2);
678 t_cross[ee](2) = t_node_diff_ksi[ei0](0) * t_node_diff_ksi[ei1](1) -
679 t_node_diff_ksi[ei0](1) * t_node_diff_ksi[ei1](0);
680 FTensor::Tensor1<double, 3> t_diff_o(
681 diffN[3 * o[ee] + 0], diffN[3 * o[ee] + 1], diffN[3 * o[ee] + 2]);
682 t_diff_cross(i, j) += t_cross[ee](i) * t_diff_o(j);
683 // cerr << t_cross[ee](0) << " " << t_cross[ee](1) << " " <<
684 // t_cross[ee](2) << endl;
685 }
686 // cerr << endl << endl;
687 t_diff_n0_p_n1(i) = t_node_diff_ksi[i0](i) + t_node_diff_ksi[i1](i);
688 t_diff_n0_p_n1_p_n2(i) = t_diff_n0_p_n1(i) + t_node_diff_ksi[i2](i);
689 } else {
690 for (int ee = 0; ee != 3; ee++) {
691 t_cross[ee](0) = 1;
692 t_cross[ee](1) = 0;
693 t_cross[ee](2) = 0;
694 }
695 }
696
697 FTensor::Tensor1<double *, 3> t_phi(&phi_f[HVEC0], &phi_f[HVEC1],
698 &phi_f[HVEC2], 3);
699 boost::shared_ptr<FTensor::Tensor2<double *, 3, 3>> t_diff_phi_ptr;
700 if (diff_phi_f) {
701 t_diff_phi_ptr = boost::shared_ptr<FTensor::Tensor2<double *, 3, 3>>(
702 new FTensor::Tensor2<double *, 3, 3>(
703 &diff_phi_f[HVEC0_0], &diff_phi_f[HVEC0_1], &diff_phi_f[HVEC0_2],
704 &diff_phi_f[HVEC1_0], &diff_phi_f[HVEC1_1], &diff_phi_f[HVEC1_2],
705 &diff_phi_f[HVEC2_0], &diff_phi_f[HVEC2_1], &diff_phi_f[HVEC2_2],
706 9));
707 }
708
709 double fi[p + 1], diff_fi[3 * p + 3];
710 double fj[p + 1], diff_fj[3 * p + 3];
711 double tmp_fj[p + 1], tmp_diff_fj[3 * p + 3];
712 for (int ii = 0; ii != gdim; ii++) {
713 const int shift = ii * nb;
714 double n0 = N[shift + i0];
715 double n1 = N[shift + i1];
716 double n2 = N[shift + i2];
717 double *diff_n1 = (diff_phi_f) ? &t_node_diff_ksi[i1](0) : NULL;
718 double *diff_n0_p_n1 = (diff_phi_f) ? &t_diff_n0_p_n1(0) : NULL;
719 ierr = Jacobi_polynomials(p, 0, n1, n0 + n1, diff_n1, diff_n0_p_n1, fi,
720 diff_phi_f ? diff_fi : NULL, 3);
721 CHKERRG(ierr);
722 for (int pp = 0; pp <= p; pp++) {
723 double *diff_n2 = (diff_phi_f) ? &t_node_diff_ksi[i2](0) : NULL;
724 double *diff_n0_p_n1_p_n2 = (diff_phi_f) ? &t_diff_n0_p_n1_p_n2(0) : NULL;
725 ierr = Jacobi_polynomials(pp, 2 * pp + 1, n2, n0 + n1 + n2, diff_n2,
726 diff_n0_p_n1_p_n2, tmp_fj,
727 diff_phi_f ? tmp_diff_fj : NULL, 3);
728 CHKERRG(ierr);
729 fj[pp] = tmp_fj[pp];
730 if (diff_phi_f) {
731 diff_fj[0 * (p + 1) + pp] = tmp_diff_fj[0 * (pp + 1) + pp];
732 diff_fj[1 * (p + 1) + pp] = tmp_diff_fj[1 * (pp + 1) + pp];
733 diff_fj[2 * (p + 1) + pp] = tmp_diff_fj[2 * (pp + 1) + pp];
734 }
735 }
736 double no0 = N[shift + o[0]];
737 double no1 = N[shift + o[1]];
738 double no2 = N[shift + o[2]];
739 FTensor::Tensor1<double, 3> base0;
740 base0(i) = no0 * t_cross[0](i) + no1 * t_cross[1](i) + no2 * t_cross[2](i);
741 int jj = 0;
742 for (int oo = 0; oo < p; oo++) {
743 FTensor::Tensor0<double *> t_fi(fi);
744 FTensor::Tensor1<double *, 3> t_diff_fi(&diff_fi[0], &diff_fi[p + 1],
745 &diff_fi[2 * p + 2], 1);
746 for (int ll = 0; ll <= oo; ll++) {
747 const int mm = oo - ll;
748 if (mm >= 0) {
749 const double a = t_fi * fj[mm];
750 // cerr << ll << " " << mm << " " << a << endl;
751 t_phi(i) = a * base0(i);
752 if (diff_phi_f) {
753 FTensor::Tensor1<double *, 3> t_diff_fj(
754 &diff_fj[0 + mm], &diff_fj[p + 1 + mm],
755 &diff_fj[2 * p + 2 + mm], 1);
756 (*t_diff_phi_ptr)(i, j) =
757 a * t_diff_cross(i, j) +
758 (t_diff_fi(j) * fj[mm] + t_fi * t_diff_fj(j)) * base0(i);
759 ++(*t_diff_phi_ptr);
760 ++t_diff_fi;
761 }
762 ++t_fi;
763 ++t_phi;
764 ++jj;
765 }
766 }
767 }
768 if (jj != NBFACETRI_DEMKOWICZ_HDIV(p)) {
769 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
770 "wrong number of base functions "
771 "jj!=NBFACETRI_DEMKOWICZ_HDIV(p) "
772 "%d!=%d",
773 jj, NBFACETRI_DEMKOWICZ_HDIV(p));
774 }
775 }
776
777 MoFEMFunctionReturnHot(0);
778}
779
780MoFEMErrorCode MoFEM::Hdiv_Demkowicz_Interior_MBTET(
781 int p, double *N, double *diffN, int p_f[], double *phi_f[4],
782 double *diff_phi_f[4], double *phi_v, double *diff_phi_v, int gdim) {
783
784 const int opposite_face_node[4] = {2, 0, 1, 3};
785 // list of zero node faces
786 const int znf[] = {0, 2, 3};
787 MoFEMFunctionBegin;
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 t_node_diff_ksi[0] =
794 FTensor::Tensor1<double, 3>(diffN[0], diffN[1], diffN[2]);
795 t_node_diff_ksi[1] =
796 FTensor::Tensor1<double, 3>(diffN[3], diffN[4], diffN[5]);
797 t_node_diff_ksi[2] =
798 FTensor::Tensor1<double, 3>(diffN[6], diffN[7], diffN[8]);
799 t_node_diff_ksi[3] =
800 FTensor::Tensor1<double, 3>(diffN[9], diffN[10], diffN[11]);
801 FTensor::Tensor1<double, 3> t_m_node_diff_ksi[4];
802 for (int ff = 0; ff != 4; ++ff) {
803 t_m_node_diff_ksi[ff](i) = -t_node_diff_ksi[ff](i);
804 }
805
806 FTensor::Tensor1<double *, 3> t_phi_v(&phi_v[HVEC0], &phi_v[HVEC1],
807 &phi_v[HVEC2], 3);
808 FTensor::Tensor2<double *, 3, 3> t_diff_phi_v(
809 &diff_phi_v[HVEC0_0], &diff_phi_v[HVEC0_1], &diff_phi_v[HVEC0_2],
810 &diff_phi_v[HVEC1_0], &diff_phi_v[HVEC1_1], &diff_phi_v[HVEC1_2],
811 &diff_phi_v[HVEC2_0], &diff_phi_v[HVEC2_1], &diff_phi_v[HVEC2_2], 9);
812
813 MatrixDouble fk(3, p + 1), diff_fk(3, 3 * p + 3);
814
815 for (int ii = 0; ii != gdim; ii++) {
816 const int shift = 4 * ii;
817
818 for (int ff = 0; ff != 3; ff++) {
819 const int fff = znf[ff];
820 const int iO = opposite_face_node[fff];
821 const double nO = N[shift + iO];
822 for (int pp = 1; pp <= p; pp++) {
823 CHKERR IntegratedJacobi_polynomials(
824 pp, 2 * pp + 2, nO, 1 - nO, &t_node_diff_ksi[iO](0),
825 &t_m_node_diff_ksi[iO](0), &fk(ff, 0), &diff_fk(ff, 0), 3);
826 }
827 }
828
829 int jj = 0;
830 for (int oo = 2; oo <= p; oo++) {
831 for (int k = 1; k != oo; k++) {
832 int OO = oo - k;
833 if (OO >= 0) {
834 int s = NBFACETRI_DEMKOWICZ_HDIV(OO - 1);
835 // Note that we do faces 0,2,3, skipping 1. All the faces which have
836 // zero node in it.
837 int nb_dofs = NBFACETRI_DEMKOWICZ_HDIV(p_f[znf[0]]);
838 int sp[] = {ii * 3 * nb_dofs + 3 * s, ii * 3 * nb_dofs + 3 * s,
839 ii * 3 * nb_dofs + 3 * s};
840 FTensor::Tensor1<double *, 3> t_phi_f[] = {
841 FTensor::Tensor1<double *, 3>(&phi_f[znf[0]][sp[0] + HVEC0],
842 &phi_f[znf[0]][sp[0] + HVEC1],
843 &phi_f[znf[0]][sp[0] + HVEC2], 3),
844 FTensor::Tensor1<double *, 3>(&phi_f[znf[1]][sp[1] + HVEC0],
845 &phi_f[znf[1]][sp[1] + HVEC1],
846 &phi_f[znf[1]][sp[1] + HVEC2], 3),
847 FTensor::Tensor1<double *, 3>(&phi_f[znf[2]][sp[2] + HVEC0],
848 &phi_f[znf[2]][sp[2] + HVEC1],
849 &phi_f[znf[2]][sp[2] + HVEC2], 3)};
850 int sdp[] = {ii * 9 * nb_dofs + 9 * s, ii * 9 * nb_dofs + 9 * s,
851 ii * 9 * nb_dofs + 9 * s};
852 FTensor::Tensor2<double *, 3, 3> t_diff_phi_f[] = {
853 FTensor::Tensor2<double *, 3, 3>(
854 &diff_phi_f[znf[0]][sdp[0] + HVEC0_0],
855 &diff_phi_f[znf[0]][sdp[0] + HVEC0_1],
856 &diff_phi_f[znf[0]][sdp[0] + HVEC0_2],
857 &diff_phi_f[znf[0]][sdp[0] + HVEC1_0],
858 &diff_phi_f[znf[0]][sdp[0] + HVEC1_1],
859 &diff_phi_f[znf[0]][sdp[0] + HVEC1_2],
860 &diff_phi_f[znf[0]][sdp[0] + HVEC2_0],
861 &diff_phi_f[znf[0]][sdp[0] + HVEC2_1],
862 &diff_phi_f[znf[0]][sdp[0] + HVEC2_2], 9),
863 FTensor::Tensor2<double *, 3, 3>(
864 &diff_phi_f[znf[1]][sdp[1] + HVEC0_0],
865 &diff_phi_f[znf[1]][sdp[1] + HVEC0_1],
866 &diff_phi_f[znf[1]][sdp[1] + HVEC0_2],
867 &diff_phi_f[znf[1]][sdp[1] + HVEC1_0],
868 &diff_phi_f[znf[1]][sdp[1] + HVEC1_1],
869 &diff_phi_f[znf[1]][sdp[1] + HVEC1_2],
870 &diff_phi_f[znf[1]][sdp[1] + HVEC2_0],
871 &diff_phi_f[znf[1]][sdp[1] + HVEC2_1],
872 &diff_phi_f[znf[1]][sdp[1] + HVEC2_2], 9),
873 FTensor::Tensor2<double *, 3, 3>(
874 &diff_phi_f[znf[2]][sdp[2] + HVEC0_0],
875 &diff_phi_f[znf[2]][sdp[2] + HVEC0_1],
876 &diff_phi_f[znf[2]][sdp[2] + HVEC0_2],
877 &diff_phi_f[znf[2]][sdp[2] + HVEC1_0],
878 &diff_phi_f[znf[2]][sdp[2] + HVEC1_1],
879 &diff_phi_f[znf[2]][sdp[2] + HVEC1_2],
880 &diff_phi_f[znf[2]][sdp[2] + HVEC2_0],
881 &diff_phi_f[znf[2]][sdp[2] + HVEC2_1],
882 &diff_phi_f[znf[2]][sdp[2] + HVEC2_2], 9)};
883 for (int ij = s; ij != NBFACETRI_DEMKOWICZ_HDIV(OO); ij++) {
884 for (int ff = 0; ff != 3; ff++) {
885 FTensor::Tensor1<double, 3> t_diff_fk(diff_fk(ff, 0 * p + k - 1),
886 diff_fk(ff, 1 * p + k - 1),
887 diff_fk(ff, 2 * p + k - 1));
888 t_phi_v(i) = fk(ff, k - 1) * t_phi_f[ff](i);
889 t_diff_phi_v(i, j) = t_diff_fk(j) * t_phi_f[ff](i) +
890 fk(ff, k - 1) * t_diff_phi_f[ff](i, j);
891 ++t_phi_v;
892 ++t_diff_phi_v;
893 ++t_phi_f[ff];
894 ++t_diff_phi_f[ff];
895 ++jj;
896 }
897 }
898 }
899 }
900 }
901 if (jj != NBVOLUMETET_DEMKOWICZ_HDIV(p)) {
902 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
903 "wrong number of base functions "
904 "jj!=NBVOLUMETET_DEMKOWICZ_HDIV(p) "
905 "%d!=%d",
906 jj, NBVOLUMETET_DEMKOWICZ_HDIV(p));
907 }
908 }
909
910 MoFEMFunctionReturn(0);
911}
FTENSOR_INDEX
#define FTENSOR_INDEX(DIM, I)
Definition Templates.hpp:2049
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
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
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition definitions.h:499
HVEC0
@ HVEC0
Definition definitions.h:199
HVEC1
@ HVEC1
Definition definitions.h:199
HVEC2
@ HVEC2
Definition definitions.h:199
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
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
HVEC1_2
@ HVEC1_2
Definition definitions.h:212
HVEC2_2
@ HVEC2_2
Definition definitions.h:213
HVEC2_0
@ HVEC2_0
Definition definitions.h:207
HVEC0_2
@ HVEC0_2
Definition definitions.h:211
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
NBVOLUMETET_DEMKOWICZ_HDIV
#define NBVOLUMETET_DEMKOWICZ_HDIV(P)
Definition h1_hdiv_hcurl_l2.h:140
NBVOLUMETET_AINSWORTH_FACE_HDIV
#define NBVOLUMETET_AINSWORTH_FACE_HDIV(P)
Definition h1_hdiv_hcurl_l2.h:134
NBFACETRI_DEMKOWICZ_HDIV
#define NBFACETRI_DEMKOWICZ_HDIV(P)
Definition h1_hdiv_hcurl_l2.h:139
NBFACETRI_AINSWORTH_FACE_HDIV
#define NBFACETRI_AINSWORTH_FACE_HDIV(P)
Definition h1_hdiv_hcurl_l2.h:131
NBVOLUMETET_AINSWORTH_VOLUME_HDIV
#define NBVOLUMETET_AINSWORTH_VOLUME_HDIV(P)
Definition h1_hdiv_hcurl_l2.h:135
NBFACETRI_AINSWORTH_EDGE_HDIV
#define NBFACETRI_AINSWORTH_EDGE_HDIV(P)
Definition h1_hdiv_hcurl_l2.h:130
i
FTensor::Index< 'i', SPACE_DIM > i
Definition hcurl_divergence_operator_2d.cpp:27
lambda
static double lambda
Definition incompressible_elasticity.cpp:181
n
const double n
refractive index of diffusive medium
Definition initial_diffusion.cpp:38
l
FTensor::Index< 'l', 3 > l
Definition matrix_function.cpp:21
j
FTensor::Index< 'j', 3 > j
Definition matrix_function.cpp:19
k
FTensor::Index< 'k', 3 > k
Definition matrix_function.cpp:20
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition Exceptions.hpp:76
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition Exceptions.hpp:56
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::Hdiv_Ainsworth_FaceBubbleShapeFunctions
MoFEMErrorCode Hdiv_Ainsworth_FaceBubbleShapeFunctions(int *faces_nodes, int *p, double *N, double *diffN, double *phi_f[], double *diff_phi_f[], int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Face bubble functions by Ainsworth nme:nme847.
Definition Hdiv.cpp:151
MoFEM::Hdiv_Demkowicz_Face_MBTET_ON_FACE
MoFEMErrorCode Hdiv_Demkowicz_Face_MBTET_ON_FACE(int *faces_nodes, int p, double *N, double *diffN, double *phi_f, double *diff_phi_f, int gdim, int nb)
Definition Hdiv.cpp:634
MoFEM::Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET
MoFEMErrorCode Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET(int *faces_nodes, int *p, double *N, double *diffN, double *phi_f_e[4][3], double *diff_phi_f_e[4][3], int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Hdiv base functions, Edge-based face functions by Ainsworth nme:nme847.
Definition Hdiv.cpp:24
MoFEM::Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET_ON_FACE
MoFEMErrorCode Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET_ON_FACE(int *faces_nodes, int p, double *N, double *diffN, double *phi_f_e[3], double *diff_phi_f_e[3], int gdim, int nb, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Hdiv base functions, Edge-based face functions by Ainsworth nme:nme847.
Definition Hdiv.cpp:47
MoFEM::Hdiv_Demkowicz_Interior_MBTET
MoFEMErrorCode Hdiv_Demkowicz_Interior_MBTET(int p, double *N, double *diffN, int p_face[], double *phi_f[4], double *diff_phi_f[4], double *phi_v, double *diff_phi_v, int gdim)
Definition Hdiv.cpp:780
MoFEM::Hdiv_Ainsworth_FaceBasedVolumeShapeFunctions_MBTET
MoFEMErrorCode Hdiv_Ainsworth_FaceBasedVolumeShapeFunctions_MBTET(int p, double *N, double *diffN, double *phi_v_f[], double *diff_phi_v_f[], int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Definition Hdiv.cpp:401
MoFEM::Hdiv_Ainsworth_FaceBubbleShapeFunctions_ON_FACE
MoFEMErrorCode Hdiv_Ainsworth_FaceBubbleShapeFunctions_ON_FACE(int *faces_nodes, int p, double *N, double *diffN, double *phi_f, double *diff_phi_f, int gdim, int nb, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Face bubble functions by Ainsworth nme:nme847.
Definition Hdiv.cpp:174
MoFEM::getFTensor1FromPtr
FTensor::Tensor1< FTensor::PackPtr< double *, S >, DIM > getFTensor1FromPtr(double *ptr)
Make Tensor1 from pointer.
Definition Templates.hpp:838
MoFEM::Hdiv_Ainsworth_EdgeBasedVolumeShapeFunctions_MBTET
MoFEMErrorCode Hdiv_Ainsworth_EdgeBasedVolumeShapeFunctions_MBTET(int p, double *N, double *diffN, double *phi_v_e[6], double *diff_phi_v_e[6], int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Hdiv base function, Edge-based interior (volume) functions by Ainsworth nme:nme847.
Definition Hdiv.cpp:307
MoFEM::Hdiv_Ainsworth_VolumeBubbleShapeFunctions_MBTET
MoFEMErrorCode Hdiv_Ainsworth_VolumeBubbleShapeFunctions_MBTET(int p, double *N, double *diffN, double *phi_v, double *diff_phi_v, int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Interior bubble functions by Ainsworth nme:nme847.
Definition Hdiv.cpp:507
m
FTensor::Index< 'm', 3 > m
Definition shallow_wave.cpp:80
N
const int N
Definition speed_test.cpp:3
MoFEM::AinsworthOrderHooks::broken_nbvolumetet_edge_hdiv
static boost::function< int(int)> broken_nbvolumetet_edge_hdiv
Definition Hdiv.hpp:27
MoFEM::AinsworthOrderHooks::broken_nbvolumetet_face_hdiv
static boost::function< int(int)> broken_nbvolumetet_face_hdiv
Definition Hdiv.hpp:28
MoFEM::AinsworthOrderHooks::broken_nbfacetri_face_hdiv
static boost::function< int(int)> broken_nbfacetri_face_hdiv
Definition Hdiv.hpp:26
MoFEM::AinsworthOrderHooks::broken_nbvolumetet_volume_hdiv
static boost::function< int(int)> broken_nbvolumetet_volume_hdiv
Definition Hdiv.hpp:29
MoFEM::AinsworthOrderHooks::broken_nbfacetri_edge_hdiv
static boost::function< int(int)> broken_nbfacetri_edge_hdiv
Definition Hdiv.hpp:25
Generated by Doxygen 1.12.0 and hosted at University of Glasgow