v0.16.0
Loading...
Searching...
No Matches
TetPolynomialBase.cpp
Go to the documentation of this file.
1/** \file TetPolynomialBase.cpp
2\brief Implementation of hierarchical bases on tetrahedral
3
4A l2, h1, h-div and h-curl spaces are implemented.
5
6*/
7
8using namespace MoFEM;
9
10struct TetBaseCache {
11
12 struct BaseCacheItem {
13 int order;
14 int nb_gauss_pts;
15 mutable MatrixDouble N;
16 mutable MatrixDouble diffN;
17 };
18
19 struct HDivBaseCacheItem {
20
21 int order;
22 int nb_gauss_pts;
23
24 // Number of permeations for tetrahedron
25 // That is P(3, 4) = 24
26
27 int n0;
28 int n1;
29 int n2;
30
31 mutable MatrixDouble N;
32 mutable MatrixDouble diffN;
33 };
34
35 using BaseCacheMI = boost::multi_index_container<
36 BaseCacheItem,
37 boost::multi_index::indexed_by<
38
39 boost::multi_index::hashed_unique<
40
41 composite_key<
42
43 BaseCacheItem,
44 member<BaseCacheItem, int, &BaseCacheItem::order>,
45 member<BaseCacheItem, int, &BaseCacheItem::nb_gauss_pts>>>>
46
47 >;
48
49 using HDivBaseFaceCacheMI = boost::multi_index_container<
50 HDivBaseCacheItem,
51 boost::multi_index::indexed_by<
52
53 boost::multi_index::hashed_unique<
54
55 composite_key<
56
57 HDivBaseCacheItem,
58 member<HDivBaseCacheItem, int, &HDivBaseCacheItem::order>,
59 member<HDivBaseCacheItem, int,
60 &HDivBaseCacheItem::nb_gauss_pts>,
61 member<HDivBaseCacheItem, int, &HDivBaseCacheItem::n0>,
62 member<HDivBaseCacheItem, int, &HDivBaseCacheItem::n1>,
63 member<HDivBaseCacheItem, int, &HDivBaseCacheItem::n2>>>>
64
65 >;
66
67 static std::array<std::map<const void *, BaseCacheMI>, LASTBASE>
68 hdivBaseInterior;
69 static std::array<std::map<const void *, BaseCacheMI>, LASTBASE>
70 hdivBrokenBaseInterior;
71 static std::array<std::map<const void *, HDivBaseFaceCacheMI>, LASTBASE>
72 hDivBaseFace;
73 static std::array<std::map<const void *, BaseCacheMI>, LASTBASE>
74 l2BaseInterior;
75};
76
77std::array<std::map<const void *, TetBaseCache::BaseCacheMI>, LASTBASE>
78 TetBaseCache::hdivBaseInterior;
79std::array<std::map<const void *, TetBaseCache::BaseCacheMI>, LASTBASE>
80 TetBaseCache::hdivBrokenBaseInterior;
81std::array<std::map<const void *, TetBaseCache::HDivBaseFaceCacheMI>, LASTBASE>
82 TetBaseCache::hDivBaseFace;
83std::array<std::map<const void *, TetBaseCache::BaseCacheMI>, LASTBASE>
84 TetBaseCache::l2BaseInterior;
85
86MoFEMErrorCode
87TetPolynomialBase::query_interface(boost::typeindex::type_index type_index,
88 UnknownInterface **iface) const {
89
90 MoFEMFunctionBeginHot;
91 *iface = const_cast<TetPolynomialBase *>(this);
92 MoFEMFunctionReturnHot(0);
93}
94
95TetPolynomialBase::TetPolynomialBase(const void *ptr) : vPtr(ptr) {}
96
97TetPolynomialBase::~TetPolynomialBase() {
98 if (vPtr) {
99
100 auto erase = [&](auto cache) {
101 if (cache.find(vPtr) != cache.end())
102 cache.erase(vPtr);
103 };
104
105 for (auto b = 0; b != LASTBASE; ++b) {
106 erase(TetBaseCache::hdivBaseInterior[b]);
107 erase(TetBaseCache::hdivBrokenBaseInterior[b]);
108 erase(TetBaseCache::hDivBaseFace[b]);
109 }
110 }
111}
112
113MoFEMErrorCode TetPolynomialBase::getValueH1(MatrixDouble &pts) {
114 MoFEMFunctionBegin;
115
116 switch (cTx->bAse) {
117 case AINSWORTH_LEGENDRE_BASE:
118 case AINSWORTH_LOBATTO_BASE:
119 CHKERR getValueH1AinsworthBase(pts);
120 break;
121 case AINSWORTH_BERNSTEIN_BEZIER_BASE:
122 CHKERR getValueH1BernsteinBezierBase(pts);
123 break;
124 default:
125 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
126 }
127
128 MoFEMFunctionReturn(0);
129}
130
131MoFEMErrorCode TetPolynomialBase::getValueH1AinsworthBase(MatrixDouble &pts) {
132 MoFEMFunctionBegin;
133
134 EntitiesFieldData &data = cTx->dAta;
135 const FieldApproximationBase base = cTx->bAse;
136 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
137 double *diffL, const int dim) =
138 cTx->basePolynomialsType0;
139
140 int nb_gauss_pts = pts.size2();
141
142 int sense[6], order[6];
143 if (data.spacesOnEntities[MBEDGE].test(H1)) {
144 // edges
145 if (data.dataOnEntities[MBEDGE].size() != 6) {
146 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
147 }
148 double *h1_edge_n[6], *diff_h1_egde_n[6];
149 for (int ee = 0; ee != 6; ++ee) {
150 if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
151 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
152 "data inconsistency");
153 }
154 sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
155 order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
156 int nb_dofs = NBEDGE_H1(data.dataOnEntities[MBEDGE][ee].getOrder());
157 data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, nb_dofs,
158 false);
159 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
160 3 * nb_dofs, false);
161 h1_edge_n[ee] =
162 &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
163 diff_h1_egde_n[ee] =
164 &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
165 }
166 CHKERR H1_EdgeShapeFunctions_MBTET(
167 sense, order,
168 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
169 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
170 h1_edge_n, diff_h1_egde_n, nb_gauss_pts, base_polynomials);
171 } else {
172 for (int ee = 0; ee != 6; ++ee) {
173 data.dataOnEntities[MBEDGE][ee].getN(base).resize(0, 0, false);
174 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(0, 0, false);
175 }
176 }
177
178 if (data.spacesOnEntities[MBTRI].test(H1)) {
179 // faces
180 if (data.dataOnEntities[MBTRI].size() != 4) {
181 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
182 }
183 double *h1_face_n[4], *diff_h1_face_n[4];
184 for (int ff = 0; ff != 4; ++ff) {
185 if (data.dataOnEntities[MBTRI][ff].getSense() == 0) {
186 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
187 "data inconsistency");
188 }
189 int nb_dofs = NBFACETRI_H1(data.dataOnEntities[MBTRI][ff].getOrder());
190 order[ff] = data.dataOnEntities[MBTRI][ff].getOrder();
191 data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, nb_dofs,
192 false);
193 data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
194 3 * nb_dofs, false);
195 h1_face_n[ff] =
196 &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
197 diff_h1_face_n[ff] =
198 &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
199 }
200 if (data.facesNodes.size1() != 4) {
201 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
202 }
203 if (data.facesNodes.size2() != 3) {
204 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
205 }
206 CHKERR H1_FaceShapeFunctions_MBTET(
207 &*data.facesNodes.data().begin(), order,
208 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
209 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
210 h1_face_n, diff_h1_face_n, nb_gauss_pts, base_polynomials);
211
212 } else {
213 for (int ff = 0; ff != 4; ++ff) {
214 data.dataOnEntities[MBTRI][ff].getN(base).resize(0, false);
215 data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(0, 0, false);
216 }
217 }
218
219 if (data.spacesOnEntities[MBTET].test(H1)) {
220 // volume
221 int order = data.dataOnEntities[MBTET][0].getOrder();
222 int nb_vol_dofs = NBVOLUMETET_H1(order);
223 data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, nb_vol_dofs,
224 false);
225 data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
226 3 * nb_vol_dofs, false);
227 CHKERR H1_VolumeShapeFunctions_MBTET(
228 data.dataOnEntities[MBTET][0].getOrder(),
229 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
230 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
231 &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
232 &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
233 nb_gauss_pts, base_polynomials);
234 } else {
235 data.dataOnEntities[MBTET][0].getN(base).resize(0, 0, false);
236 data.dataOnEntities[MBTET][0].getDiffN(base).resize(0, 0, false);
237 }
238
239 MoFEMFunctionReturn(0);
240}
241
242MoFEMErrorCode
243TetPolynomialBase::getValueH1BernsteinBezierBase(MatrixDouble &pts) {
244 MoFEMFunctionBegin;
245
246 EntitiesFieldData &data = cTx->dAta;
247 const std::string field_name = cTx->fieldName;
248 const int nb_gauss_pts = pts.size2();
249
250 if (data.dataOnEntities[MBVERTEX][0].getN(NOBASE).size1() !=
251 (unsigned int)nb_gauss_pts)
252 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
253 "Base functions or nodes has wrong number of integration points "
254 "for base %s",
255 ApproximationBaseNames[NOBASE]);
256 auto &lambda = data.dataOnEntities[MBVERTEX][0].getN(NOBASE);
257
258 auto get_alpha = [field_name](auto &data) -> MatrixInt & {
259 auto &ptr = data.getBBAlphaIndicesSharedPtr(field_name);
260 if (!ptr)
261 ptr.reset(new MatrixInt());
262 return *ptr;
263 };
264
265 auto get_base = [field_name](auto &data) -> MatrixDouble & {
266 auto &ptr = data.getBBNSharedPtr(field_name);
267 if (!ptr)
268 ptr.reset(new MatrixDouble());
269 return *ptr;
270 };
271
272 auto get_diff_base = [field_name](auto &data) -> MatrixDouble & {
273 auto &ptr = data.getBBDiffNSharedPtr(field_name);
274 if (!ptr)
275 ptr.reset(new MatrixDouble());
276 return *ptr;
277 };
278
279 auto get_alpha_by_name_ptr =
280 [](auto &data,
281 const std::string &field_name) -> boost::shared_ptr<MatrixInt> & {
282 return data.getBBAlphaIndicesSharedPtr(field_name);
283 };
284
285 auto get_base_by_name_ptr =
286 [](auto &data,
287 const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
288 return data.getBBNSharedPtr(field_name);
289 };
290
291 auto get_diff_base_by_name_ptr =
292 [](auto &data,
293 const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
294 return data.getBBDiffNSharedPtr(field_name);
295 };
296
297 auto get_alpha_by_order_ptr =
298 [](auto &data, const size_t o) -> boost::shared_ptr<MatrixInt> & {
299 return data.getBBAlphaIndicesByOrderSharedPtr(o);
300 };
301
302 auto get_base_by_order_ptr =
303 [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
304 return data.getBBNByOrderSharedPtr(o);
305 };
306
307 auto get_diff_base_by_order_ptr =
308 [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
309 return data.getBBDiffNByOrderSharedPtr(o);
310 };
311
312 auto &vert_ent_data = data.dataOnEntities[MBVERTEX][0];
313 auto &vertex_alpha = get_alpha(vert_ent_data);
314 vertex_alpha.resize(4, 4, false);
315 vertex_alpha.clear();
316 for (int n = 0; n != 4; ++n)
317 vertex_alpha(n, n) = data.dataOnEntities[MBVERTEX][0].getBBNodeOrder()[n];
318
319 auto &vert_get_n = get_base(vert_ent_data);
320 auto &vert_get_diff_n = get_diff_base(vert_ent_data);
321 vert_get_n.resize(nb_gauss_pts, 4, false);
322 vert_get_diff_n.resize(nb_gauss_pts, 12, false);
323 CHKERR BernsteinBezier::baseFunctionsTet(
324 1, lambda.size1(), vertex_alpha.size1(), &vertex_alpha(0, 0),
325 &lambda(0, 0), Tools::diffShapeFunMBTET.data(), &vert_get_n(0, 0),
326 &vert_get_diff_n(0, 0));
327 for (int n = 0; n != 4; ++n) {
328 const double f = boost::math::factorial<double>(
329 data.dataOnEntities[MBVERTEX][0].getBBNodeOrder()[n]);
330 for (int g = 0; g != nb_gauss_pts; ++g) {
331 vert_get_n(g, n) *= f;
332 for (int d = 0; d != 3; ++d)
333 vert_get_diff_n(g, 3 * n + d) *= f;
334 }
335 }
336
337 // edges
338 if (data.spacesOnEntities[MBEDGE].test(H1)) {
339 if (data.dataOnEntities[MBEDGE].size() != 6)
340 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
341 "Wrong size of ent data");
342
343 constexpr int edges_nodes[6][2] = {{0, 1}, {1, 2}, {2, 0},
344 {0, 3}, {1, 3}, {2, 3}};
345 for (int ee = 0; ee != 6; ++ee) {
346 auto &ent_data = data.dataOnEntities[MBEDGE][ee];
347 const int sense = ent_data.getSense();
348 if (sense == 0)
349 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
350 "Sense of the edge unknown");
351 const int order = ent_data.getOrder();
352 const int nb_dofs = NBEDGE_H1(order);
353
354 if (nb_dofs) {
355 if (get_alpha_by_order_ptr(ent_data, order)) {
356 get_alpha_by_name_ptr(ent_data, field_name) =
357 get_alpha_by_order_ptr(ent_data, order);
358 get_base_by_name_ptr(ent_data, field_name) =
359 get_base_by_order_ptr(ent_data, order);
360 get_diff_base_by_name_ptr(ent_data, field_name) =
361 get_diff_base_by_order_ptr(ent_data, order);
362 } else {
363 auto &get_n = get_base(ent_data);
364 auto &get_diff_n = get_diff_base(ent_data);
365 get_n.resize(nb_gauss_pts, nb_dofs, false);
366 get_diff_n.resize(nb_gauss_pts, 3 * nb_dofs, false);
367
368 auto &edge_alpha = get_alpha(data.dataOnEntities[MBEDGE][ee]);
369 edge_alpha.resize(nb_dofs, 4, false);
370 CHKERR BernsteinBezier::generateIndicesEdgeTet(ee, order,
371 &edge_alpha(0, 0));
372 if (sense == -1) {
373 for (int i = 0; i != edge_alpha.size1(); ++i) {
374 int a = edge_alpha(i, edges_nodes[ee][0]);
375 edge_alpha(i, edges_nodes[ee][0]) =
376 edge_alpha(i, edges_nodes[ee][1]);
377 edge_alpha(i, edges_nodes[ee][1]) = a;
378 }
379 }
380 CHKERR BernsteinBezier::baseFunctionsTet(
381 order, lambda.size1(), edge_alpha.size1(), &edge_alpha(0, 0),
382 &lambda(0, 0), Tools::diffShapeFunMBTET.data(), &get_n(0, 0),
383 &get_diff_n(0, 0));
384
385 get_alpha_by_order_ptr(ent_data, order) =
386 get_alpha_by_name_ptr(ent_data, field_name);
387 get_base_by_order_ptr(ent_data, order) =
388 get_base_by_name_ptr(ent_data, field_name);
389 get_diff_base_by_order_ptr(ent_data, order) =
390 get_diff_base_by_name_ptr(ent_data, field_name);
391 }
392 }
393 }
394 } else {
395 for (int ee = 0; ee != 6; ++ee) {
396 auto &ent_data = data.dataOnEntities[MBEDGE][ee];
397 ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
398 auto &get_n = get_base(ent_data);
399 auto &get_diff_n = get_diff_base(ent_data);
400 get_n.resize(nb_gauss_pts, 0, false);
401 get_diff_n.resize(nb_gauss_pts, 0, false);
402 }
403 }
404
405 // face
406 if (data.spacesOnEntities[MBTRI].test(H1)) {
407 if (data.dataOnEntities[MBTRI].size() != 4)
408 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
409 "Wrong size of ent data");
410 if (data.facesNodes.size1() != 4)
411 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
412 if (data.facesNodes.size2() != 3)
413 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
414
415 for (int ff = 0; ff != 4; ++ff) {
416 auto &ent_data = data.dataOnEntities[MBTRI][ff];
417 const int order = ent_data.getOrder();
418 const int nb_dofs = NBFACETRI_H1(order);
419
420 if (nb_dofs) {
421 if (get_alpha_by_order_ptr(ent_data, order)) {
422 get_alpha_by_name_ptr(ent_data, field_name) =
423 get_alpha_by_order_ptr(ent_data, order);
424 get_base_by_name_ptr(ent_data, field_name) =
425 get_base_by_order_ptr(ent_data, order);
426 get_diff_base_by_name_ptr(ent_data, field_name) =
427 get_diff_base_by_order_ptr(ent_data, order);
428 } else {
429
430 auto &get_n = get_base(ent_data);
431 auto &get_diff_n = get_diff_base(ent_data);
432 get_n.resize(nb_gauss_pts, nb_dofs, false);
433 get_diff_n.resize(nb_gauss_pts, 3 * nb_dofs, false);
434
435 auto &face_alpha = get_alpha(ent_data);
436 face_alpha.resize(nb_dofs, 4, false);
437
438 CHKERR BernsteinBezier::generateIndicesTriTet(ff, order,
439 &face_alpha(0, 0));
440 senseFaceAlpha.resize(face_alpha.size1(), face_alpha.size2(), false);
441 senseFaceAlpha.clear();
442 constexpr int tri_nodes[4][3] = {
443 {0, 1, 3}, {1, 2, 3}, {0, 2, 3}, {0, 1, 2}};
444 for (int d = 0; d != nb_dofs; ++d)
445 for (int n = 0; n != 3; ++n)
446 senseFaceAlpha(d, data.facesNodes(ff, n)) =
447 face_alpha(d, tri_nodes[ff][n]);
448 face_alpha.swap(senseFaceAlpha);
449 CHKERR BernsteinBezier::baseFunctionsTet(
450 order, lambda.size1(), face_alpha.size1(), &face_alpha(0, 0),
451 &lambda(0, 0), Tools::diffShapeFunMBTET.data(), &get_n(0, 0),
452 &get_diff_n(0, 0));
453
454 get_alpha_by_order_ptr(ent_data, order) =
455 get_alpha_by_name_ptr(ent_data, field_name);
456 get_base_by_order_ptr(ent_data, order) =
457 get_base_by_name_ptr(ent_data, field_name);
458 get_diff_base_by_order_ptr(ent_data, order) =
459 get_diff_base_by_name_ptr(ent_data, field_name);
460 }
461 }
462 }
463 } else {
464 for (int ff = 0; ff != 4; ++ff) {
465 auto &ent_data = data.dataOnEntities[MBTRI][ff];
466 ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
467 auto &get_n = get_base(ent_data);
468 auto &get_diff_n = get_diff_base(ent_data);
469 get_n.resize(nb_gauss_pts, 0, false);
470 get_diff_n.resize(nb_gauss_pts, 0, false);
471 }
472 }
473
474 if (data.spacesOnEntities[MBTET].test(H1)) {
475 if (data.dataOnEntities[MBTET].size() != 1)
476 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
477 "Wrong size ent of ent data");
478
479 auto &ent_data = data.dataOnEntities[MBTET][0];
480 const int order = ent_data.getOrder();
481 const int nb_dofs = NBVOLUMETET_H1(order);
482 if (get_alpha_by_order_ptr(ent_data, order)) {
483 get_alpha_by_name_ptr(ent_data, field_name) =
484 get_alpha_by_order_ptr(ent_data, order);
485 get_base_by_name_ptr(ent_data, field_name) =
486 get_base_by_order_ptr(ent_data, order);
487 get_diff_base_by_name_ptr(ent_data, field_name) =
488 get_diff_base_by_order_ptr(ent_data, order);
489 } else {
490
491 auto &get_n = get_base(ent_data);
492 auto &get_diff_n = get_diff_base(ent_data);
493 get_n.resize(nb_gauss_pts, nb_dofs, false);
494 get_diff_n.resize(nb_gauss_pts, 3 * nb_dofs, false);
495 if (nb_dofs) {
496 auto &tet_alpha = get_alpha(ent_data);
497 tet_alpha.resize(nb_dofs, 4, false);
498
499 CHKERR BernsteinBezier::generateIndicesTetTet(order, &tet_alpha(0, 0));
500 CHKERR BernsteinBezier::baseFunctionsTet(
501 order, lambda.size1(), tet_alpha.size1(), &tet_alpha(0, 0),
502 &lambda(0, 0), Tools::diffShapeFunMBTET.data(), &get_n(0, 0),
503 &get_diff_n(0, 0));
504
505 get_alpha_by_order_ptr(ent_data, order) =
506 get_alpha_by_name_ptr(ent_data, field_name);
507 get_base_by_order_ptr(ent_data, order) =
508 get_base_by_name_ptr(ent_data, field_name);
509 get_diff_base_by_order_ptr(ent_data, order) =
510 get_diff_base_by_name_ptr(ent_data, field_name);
511 }
512 }
513 } else {
514 auto &ent_data = data.dataOnEntities[MBTET][0];
515 ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
516 auto &get_n = get_base(ent_data);
517 auto &get_diff_n = get_diff_base(ent_data);
518 get_n.resize(nb_gauss_pts, 0, false);
519 get_diff_n.resize(nb_gauss_pts, 0, false);
520 }
521
522 MoFEMFunctionReturn(0);
523}
524
525MoFEMErrorCode TetPolynomialBase::getValueL2(MatrixDouble &pts) {
526 MoFEMFunctionBegin;
527
528 switch (cTx->bAse) {
529 case AINSWORTH_LEGENDRE_BASE:
530 case AINSWORTH_LOBATTO_BASE:
531 case DEMKOWICZ_JACOBI_BASE:
532 CHKERR getValueL2AinsworthBase(pts);
533 break;
534 case AINSWORTH_BERNSTEIN_BEZIER_BASE:
535 CHKERR getValueL2BernsteinBezierBase(pts);
536 break;
537 default:
538 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
539 }
540
541 MoFEMFunctionReturn(0);
542}
543
544MoFEMErrorCode TetPolynomialBase::getValueL2AinsworthBase(MatrixDouble &pts) {
545 MoFEMFunctionBegin;
546
547 EntitiesFieldData &data = cTx->dAta;
548 const FieldApproximationBase base = cTx->bAse;
549 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
550 double *diffL, const int dim) =
551 cTx->basePolynomialsType0;
552
553 int nb_gauss_pts = pts.size2();
554 int volume_order = data.dataOnEntities[MBTET][0].getOrder();
555 int nb_dofs = NBVOLUMETET_L2(volume_order);
556
557 data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, nb_dofs, false);
558 data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts, 3 * nb_dofs,
559 false);
560
561 if (!nb_dofs)
562 MoFEMFunctionReturnHot(0);
563
564 auto get_interior_cache = [this](auto base) -> TetBaseCache::BaseCacheMI * {
565 if (vPtr) {
566 auto it = TetBaseCache::l2BaseInterior[base].find(vPtr);
567 if (it != TetBaseCache::l2BaseInterior[base].end()) {
568 return &it->second;
569 }
570 }
571 return nullptr;
572 };
573
574 auto interior_cache_ptr = get_interior_cache(base);
575
576 if (interior_cache_ptr) {
577 auto it =
578 interior_cache_ptr->find(boost::make_tuple(volume_order, nb_gauss_pts));
579 if (it != interior_cache_ptr->end()) {
580 noalias(data.dataOnEntities[MBTET][0].getN(base)) = it->N;
581 noalias(data.dataOnEntities[MBTET][0].getDiffN(base)) = it->diffN;
582 MoFEMFunctionReturnHot(0);
583 }
584 }
585
586 CHKERR L2_Ainsworth_ShapeFunctions_MBTET(
587 data.dataOnEntities[MBTET][0].getOrder(),
588 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
589 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
590 &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
591 &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
592 nb_gauss_pts, base_polynomials);
593
594 if (interior_cache_ptr) {
595 auto p = interior_cache_ptr->emplace(
596 TetBaseCache::BaseCacheItem{volume_order, nb_gauss_pts});
597 p.first->N = data.dataOnEntities[MBTET][0].getN(base);
598 p.first->diffN = data.dataOnEntities[MBTET][0].getDiffN(base);
599 }
600
601 MoFEMFunctionReturn(0);
602}
603
604MoFEMErrorCode
605TetPolynomialBase::getValueL2BernsteinBezierBase(MatrixDouble &pts) {
606 MoFEMFunctionBegin;
607
608 EntitiesFieldData &data = cTx->dAta;
609 const std::string field_name = cTx->fieldName;
610 const int nb_gauss_pts = pts.size2();
611
612 if (data.dataOnEntities[MBVERTEX][0].getN(NOBASE).size1() !=
613 (unsigned int)nb_gauss_pts)
614 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
615 "Base functions or nodes has wrong number of integration points "
616 "for base %s",
617 ApproximationBaseNames[NOBASE]);
618 auto &lambda = data.dataOnEntities[MBVERTEX][0].getN(NOBASE);
619
620 auto get_alpha = [field_name](auto &data) -> MatrixInt & {
621 auto &ptr = data.getBBAlphaIndicesSharedPtr(field_name);
622 if (!ptr)
623 ptr.reset(new MatrixInt());
624 return *ptr;
625 };
626
627 auto get_base = [field_name](auto &data) -> MatrixDouble & {
628 auto &ptr = data.getBBNSharedPtr(field_name);
629 if (!ptr)
630 ptr.reset(new MatrixDouble());
631 return *ptr;
632 };
633
634 auto get_diff_base = [field_name](auto &data) -> MatrixDouble & {
635 auto &ptr = data.getBBDiffNSharedPtr(field_name);
636 if (!ptr)
637 ptr.reset(new MatrixDouble());
638 return *ptr;
639 };
640
641 auto get_alpha_by_name_ptr =
642 [](auto &data,
643 const std::string &field_name) -> boost::shared_ptr<MatrixInt> & {
644 return data.getBBAlphaIndicesSharedPtr(field_name);
645 };
646
647 auto get_base_by_name_ptr =
648 [](auto &data,
649 const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
650 return data.getBBNSharedPtr(field_name);
651 };
652
653 auto get_diff_base_by_name_ptr =
654 [](auto &data,
655 const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
656 return data.getBBDiffNSharedPtr(field_name);
657 };
658
659 auto get_alpha_by_order_ptr =
660 [](auto &data, const size_t o) -> boost::shared_ptr<MatrixInt> & {
661 return data.getBBAlphaIndicesByOrderSharedPtr(o);
662 };
663
664 auto get_base_by_order_ptr =
665 [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
666 return data.getBBNByOrderSharedPtr(o);
667 };
668
669 auto get_diff_base_by_order_ptr =
670 [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
671 return data.getBBDiffNByOrderSharedPtr(o);
672 };
673
674 if (data.spacesOnEntities[MBTET].test(L2)) {
675 if (data.dataOnEntities[MBTET].size() != 1)
676 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
677 "Wrong size ent of ent data");
678
679 auto &ent_data = data.dataOnEntities[MBTET][0];
680 const int order = ent_data.getOrder();
681 const int nb_dofs = NBVOLUMETET_L2(order);
682
683 if (get_alpha_by_order_ptr(ent_data, order)) {
684 get_alpha_by_name_ptr(ent_data, field_name) =
685 get_alpha_by_order_ptr(ent_data, order);
686 get_base_by_name_ptr(ent_data, field_name) =
687 get_base_by_order_ptr(ent_data, order);
688 get_diff_base_by_name_ptr(ent_data, field_name) =
689 get_diff_base_by_order_ptr(ent_data, order);
690 } else {
691
692 auto &get_n = get_base(ent_data);
693 auto &get_diff_n = get_diff_base(ent_data);
694 get_n.resize(nb_gauss_pts, nb_dofs, false);
695 get_diff_n.resize(nb_gauss_pts, 3 * nb_dofs, false);
696
697 if (nb_dofs) {
698
699 if (order == 0) {
700
701 if (nb_dofs != 1)
702 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
703 "Inconsistent number of DOFs");
704
705 auto &tri_alpha = get_alpha(ent_data);
706 tri_alpha.clear();
707 get_n(0, 0) = 1;
708 get_diff_n.clear();
709
710 } else {
711
712 if (nb_dofs != 4 + 6 * NBEDGE_H1(order) + 4 * NBFACETRI_H1(order) +
713 NBVOLUMETET_H1(order) ||
714 nb_dofs != NBVOLUMETET_L2(order))
715 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
716 "Inconsistent number of DOFs");
717
718 auto &tet_alpha = get_alpha(ent_data);
719 tet_alpha.resize(nb_dofs, 4, false);
720
721 CHKERR BernsteinBezier::generateIndicesVertexTet(order,
722 &tet_alpha(0, 0));
723 if (order > 1) {
724 std::array<int, 6> edge_n{order, order, order, order, order, order};
725 std::array<int *, 6> tet_edge_ptr{
726 &tet_alpha(4, 0),
727 &tet_alpha(4 + 1 * NBEDGE_H1(order), 0),
728 &tet_alpha(4 + 2 * NBEDGE_H1(order), 0),
729 &tet_alpha(4 + 3 * NBEDGE_H1(order), 0),
730 &tet_alpha(4 + 4 * NBEDGE_H1(order), 0),
731 &tet_alpha(4 + 5 * NBEDGE_H1(order), 0)};
732 CHKERR BernsteinBezier::generateIndicesEdgeTet(edge_n.data(),
733 tet_edge_ptr.data());
734 if (order > 2) {
735 std::array<int, 6> face_n{order, order, order, order};
736 std::array<int *, 6> tet_face_ptr{
737 &tet_alpha(4 + 6 * NBEDGE_H1(order), 0),
738 &tet_alpha(4 + 6 * NBEDGE_H1(order) + 1 * NBFACETRI_H1(order),
739 0),
740 &tet_alpha(4 + 6 * NBEDGE_H1(order) + 2 * NBFACETRI_H1(order),
741 0),
742 &tet_alpha(4 + 6 * NBEDGE_H1(order) + 3 * NBFACETRI_H1(order),
743 0),
744 };
745 CHKERR BernsteinBezier::generateIndicesTriTet(
746 face_n.data(), tet_face_ptr.data());
747 if (order > 3)
748 CHKERR BernsteinBezier::generateIndicesTetTet(
749 order,
750 &tet_alpha(
751 4 + 6 * NBEDGE_H1(order) + 4 * NBFACETRI_H1(order), 0));
752 }
753 }
754
755 CHKERR BernsteinBezier::baseFunctionsTet(
756 order, lambda.size1(), tet_alpha.size1(), &tet_alpha(0, 0),
757 &lambda(0, 0), Tools::diffShapeFunMBTET.data(), &get_n(0, 0),
758 &get_diff_n(0, 0));
759
760 get_alpha_by_order_ptr(ent_data, order) =
761 get_alpha_by_name_ptr(ent_data, field_name);
762 get_base_by_order_ptr(ent_data, order) =
763 get_base_by_name_ptr(ent_data, field_name);
764 get_diff_base_by_order_ptr(ent_data, order) =
765 get_diff_base_by_name_ptr(ent_data, field_name);
766 }
767 }
768 }
769 } else {
770 auto &ent_data = data.dataOnEntities[MBTET][0];
771 ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
772 auto &get_n = get_base(ent_data);
773 auto &get_diff_n = get_diff_base(ent_data);
774 get_n.resize(nb_gauss_pts, 0, false);
775 get_diff_n.resize(nb_gauss_pts, 0, false);
776 }
777
778 MoFEMFunctionReturn(0);
779}
780
781MoFEMErrorCode TetPolynomialBase::getValueHdivAinsworthBaseImpl(
782
783 MatrixDouble &pts,
784
785 MatrixDouble &shape_functions, MatrixDouble &diff_shape_functions,
786
787 int volume_order, std::array<int, 4> &faces_order,
788 std::array<int, 3 * 4> &faces_nodes,
789
790 boost::function<int(int)> broken_nbfacetri_edge_hdiv,
791 boost::function<int(int)> broken_nbfacetri_face_hdiv,
792 boost::function<int(int)> broken_nbvolumetet_edge_hdiv,
793 boost::function<int(int)> broken_nbvolumetet_face_hdiv,
794 boost::function<int(int)> broken_nbvolumetet_volume_hdiv
795
796) {
797 MoFEMFunctionBegin;
798
799 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
800 double *diffL, const int dim) =
801 cTx->basePolynomialsType0;
802
803 int nb_gauss_pts = pts.size2();
804 auto get_data_ptr = [](MatrixDouble &m) -> double * {
805 return (m.size1() && m.size2()) ? &*m.data().begin() : nullptr;
806 };
807
808 // face shape functions
809
810 double *phi_f_e[4][3];
811 double *phi_f[4];
812 double *diff_phi_f_e[4][3];
813 double *diff_phi_f[4];
814
815 N_face_edge.resize(4, 3, false);
816 N_face_bubble.resize(4, false);
817 diffN_face_edge.resize(4, 3, false);
818 diffN_face_bubble.resize(4, false);
819
820 for (int ff = 0; ff != 4; ++ff) {
821 const auto face_edge_dofs = NBFACETRI_AINSWORTH_EDGE_HDIV(
822 broken_nbfacetri_edge_hdiv(faces_order[ff]));
823 // three edges on face
824 for (int ee = 0; ee < 3; ee++) {
825 N_face_edge(ff, ee).resize(nb_gauss_pts, 3 * face_edge_dofs, false);
826 diffN_face_edge(ff, ee).resize(nb_gauss_pts, 9 * face_edge_dofs, false);
827 phi_f_e[ff][ee] = get_data_ptr(N_face_edge(ff, ee));
828 diff_phi_f_e[ff][ee] = get_data_ptr(diffN_face_edge(ff, ee));
829 }
830 auto face_bubble_dofs = NBFACETRI_AINSWORTH_FACE_HDIV(
831 broken_nbfacetri_face_hdiv(faces_order[ff]));
832 N_face_bubble[ff].resize(nb_gauss_pts, 3 * face_bubble_dofs, false);
833 diffN_face_bubble[ff].resize(nb_gauss_pts, 9 * face_bubble_dofs, false);
834 phi_f[ff] = get_data_ptr(N_face_bubble[ff]);
835 diff_phi_f[ff] = get_data_ptr(diffN_face_bubble[ff]);
836 }
837
838 constexpr int nb_nodes_on_tet = 4;
839
840 for (int ff = 0; ff < 4; ff++) {
841 CHKERR Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET_ON_FACE(
842 &faces_nodes[3 * ff], broken_nbfacetri_edge_hdiv(faces_order[ff]),
843 &*shape_functions.data().begin(), &*diff_shape_functions.data().begin(),
844 phi_f_e[ff], diff_phi_f_e[ff], nb_gauss_pts, nb_nodes_on_tet,
845 base_polynomials);
846 }
847
848 for (int ff = 0; ff < 4; ff++) {
849 CHKERR Hdiv_Ainsworth_FaceBubbleShapeFunctions_ON_FACE(
850 &faces_nodes[3 * ff], broken_nbfacetri_face_hdiv(faces_order[ff]),
851 &*shape_functions.data().begin(), &*diff_shape_functions.data().begin(),
852 phi_f[ff], diff_phi_f[ff], nb_gauss_pts, nb_nodes_on_tet,
853 base_polynomials);
854 }
855
856 // volume shape functions
857
858 double *phi_v_e[6];
859 double *phi_v_f[4];
860 double *phi_v;
861 double *diff_phi_v_e[6];
862 double *diff_phi_v_f[4];
863 double *diff_phi_v;
864
865 const auto volume_edge_dofs = NBVOLUMETET_AINSWORTH_EDGE_HDIV(
866 broken_nbvolumetet_edge_hdiv(volume_order));
867 N_volume_edge.resize(6, false);
868 diffN_volume_edge.resize(6, false);
869 for (int ee = 0; ee != 6; ++ee) {
870 N_volume_edge[ee].resize(nb_gauss_pts, 3 * volume_edge_dofs, false);
871 diffN_volume_edge[ee].resize(nb_gauss_pts, 9 * volume_edge_dofs, false);
872 phi_v_e[ee] = get_data_ptr(N_volume_edge[ee]);
873 diff_phi_v_e[ee] = get_data_ptr(diffN_volume_edge[ee]);
874 }
875 if (volume_edge_dofs)
876 CHKERR Hdiv_Ainsworth_EdgeBasedVolumeShapeFunctions_MBTET(
877 broken_nbvolumetet_edge_hdiv(volume_order),
878 &*shape_functions.data().begin(), &*diff_shape_functions.data().begin(),
879 phi_v_e, diff_phi_v_e, nb_gauss_pts, base_polynomials);
880
881 const auto volume_face_dofs = NBVOLUMETET_AINSWORTH_FACE_HDIV(
882 broken_nbvolumetet_face_hdiv(volume_order));
883 N_volume_face.resize(4, false);
884 diffN_volume_face.resize(4, false);
885 for (int ff = 0; ff != 4; ++ff) {
886 N_volume_face[ff].resize(nb_gauss_pts, 3 * volume_face_dofs, false);
887 diffN_volume_face[ff].resize(nb_gauss_pts, 9 * volume_face_dofs, false);
888 phi_v_f[ff] = get_data_ptr(N_volume_face[ff]);
889 diff_phi_v_f[ff] = get_data_ptr(diffN_volume_face[ff]);
890 }
891 if (volume_face_dofs)
892 CHKERR Hdiv_Ainsworth_FaceBasedVolumeShapeFunctions_MBTET(
893 broken_nbvolumetet_face_hdiv(volume_order),
894 &*shape_functions.data().begin(), &*diff_shape_functions.data().begin(),
895 phi_v_f, diff_phi_v_f, nb_gauss_pts, base_polynomials);
896
897 auto volume_bubble_dofs = NBVOLUMETET_AINSWORTH_VOLUME_HDIV(
898 broken_nbvolumetet_volume_hdiv(volume_order));
899 N_volume_bubble.resize(nb_gauss_pts, 3 * volume_bubble_dofs, false);
900 diffN_volume_bubble.resize(nb_gauss_pts, 9 * volume_bubble_dofs, false);
901 phi_v = get_data_ptr(N_volume_bubble);
902 diff_phi_v = get_data_ptr(diffN_volume_bubble);
903 if (volume_bubble_dofs)
904 CHKERR Hdiv_Ainsworth_VolumeBubbleShapeFunctions_MBTET(
905 broken_nbvolumetet_volume_hdiv(volume_order),
906 &*shape_functions.data().begin(), &*diff_shape_functions.data().begin(),
907 phi_v, diff_phi_v, nb_gauss_pts, base_polynomials);
908
909 MoFEMFunctionReturn(0);
910}
911
912MoFEMErrorCode TetPolynomialBase::getValueHdivAinsworthBase(MatrixDouble &pts) {
913 MoFEMFunctionBegin;
914
915 std::array<int, 4> faces_order;
916 std::array<int, 4 * 3> faces_nodes;
917
918 FieldApproximationBase base = cTx->bAse;
919 EntitiesFieldData &data = cTx->dAta;
920
921 int volume_order = data.dataOnEntities[MBTET][0].getOrder();
922 std::copy(data.facesNodes.data().begin(), data.facesNodes.data().end(),
923 faces_nodes.begin());
924 for (int ff = 0; ff != 4; ff++) {
925 faces_order[ff] = cTx->dAta.dataOnEntities[MBTRI][ff].getOrder();
926 }
927
928 CHKERR getValueHdivAinsworthBaseImpl(
929 pts, data.dataOnEntities[MBVERTEX][0].getN(base),
930 data.dataOnEntities[MBVERTEX][0].getDiffN(base), volume_order,
931 faces_order, faces_nodes,
932
933 AinsworthOrderHooks::broken_nbfacetri_edge_hdiv,
934 AinsworthOrderHooks::broken_nbfacetri_face_hdiv,
935 AinsworthOrderHooks::broken_nbvolumetet_edge_hdiv,
936 AinsworthOrderHooks::broken_nbvolumetet_face_hdiv,
937 AinsworthOrderHooks::broken_nbvolumetet_volume_hdiv
938
939 );
940
941 // Set shape functions into data structure Shape functions hast to be put
942 // in arrays in order which guarantee hierarchical series of degrees of
943 // freedom, i.e. in other words dofs form sub-entities has to be group
944 // by order.
945
946 FTENSOR_INDEX(3, i);
947 FTENSOR_INDEX(3, j);
948
949 int nb_gauss_pts = pts.size2();
950 auto get_data_ptr = [](MatrixDouble &m) -> double * {
951 return (m.size1() && m.size2()) ? &*m.data().begin() : nullptr;
952 };
953
954 // faces
955 if (data.dataOnEntities[MBTRI].size() != 4) {
956 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
957 }
958
959 // face-face
960 using Tensor1Ptr3 =
961 decltype(getFTensor1FromPtr<3>(static_cast<double *>(nullptr)));
962 Tensor1Ptr3 t_base_f_f[] = {
963 getFTensor1FromPtr<3>(get_data_ptr(N_face_bubble[0])),
964 getFTensor1FromPtr<3>(get_data_ptr(N_face_bubble[1])),
965 getFTensor1FromPtr<3>(get_data_ptr(N_face_bubble[2])),
966 getFTensor1FromPtr<3>(get_data_ptr(N_face_bubble[3]))};
967 FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_f_f[] = {
968 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_face_bubble[0])),
969 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_face_bubble[1])),
970 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_face_bubble[2])),
971 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_face_bubble[3]))};
972 // face-edge
973 Tensor1Ptr3 t_base_f_e[] = {
974 getFTensor1FromPtr<3>(get_data_ptr(N_face_edge(0, 0))),
975 getFTensor1FromPtr<3>(get_data_ptr(N_face_edge(0, 1))),
976 getFTensor1FromPtr<3>(get_data_ptr(N_face_edge(0, 2))),
977 getFTensor1FromPtr<3>(get_data_ptr(N_face_edge(1, 0))),
978 getFTensor1FromPtr<3>(get_data_ptr(N_face_edge(1, 1))),
979 getFTensor1FromPtr<3>(get_data_ptr(N_face_edge(1, 2))),
980 getFTensor1FromPtr<3>(get_data_ptr(N_face_edge(2, 0))),
981 getFTensor1FromPtr<3>(get_data_ptr(N_face_edge(2, 1))),
982 getFTensor1FromPtr<3>(get_data_ptr(N_face_edge(2, 2))),
983 getFTensor1FromPtr<3>(get_data_ptr(N_face_edge(3, 0))),
984 getFTensor1FromPtr<3>(get_data_ptr(N_face_edge(3, 1))),
985 getFTensor1FromPtr<3>(get_data_ptr(N_face_edge(3, 2)))};
986 FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_f_e[] = {
987 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_face_edge(0, 0))),
988 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_face_edge(0, 1))),
989 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_face_edge(0, 2))),
990 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_face_edge(1, 0))),
991 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_face_edge(1, 1))),
992 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_face_edge(1, 2))),
993 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_face_edge(2, 0))),
994 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_face_edge(2, 1))),
995 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_face_edge(2, 2))),
996 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_face_edge(3, 0))),
997 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_face_edge(3, 1))),
998 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_face_edge(3, 2)))};
999
1000 for (int ff = 0; ff != 4; ff++) {
1001 int face_order = faces_order[ff];
1002 auto face_dofs =
1003 3 * NBFACETRI_AINSWORTH_EDGE_HDIV(
1004 AinsworthOrderHooks::broken_nbfacetri_edge_hdiv(face_order)) +
1005 NBFACETRI_AINSWORTH_FACE_HDIV(
1006 AinsworthOrderHooks::broken_nbfacetri_face_hdiv(face_order));
1007 data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts,
1008 3 * face_dofs, false);
1009 data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
1010 9 * face_dofs, false);
1011 if (face_dofs) {
1012 double *base_ptr =
1013 &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
1014 double *diff_base_ptr =
1015 &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
1016 auto t_base = getFTensor1FromPtr<3>(base_ptr);
1017 auto t_diff_base = getFTensor2HVecFromPtr<3, 3>(diff_base_ptr);
1018
1019 auto max_face_order =
1020 std::max(face_order,
1021 AinsworthOrderHooks::broken_nbfacetri_edge_hdiv(face_order));
1022 max_face_order =
1023 std::max(max_face_order,
1024 AinsworthOrderHooks::broken_nbfacetri_face_hdiv(face_order));
1025
1026 for (int gg = 0; gg != nb_gauss_pts; gg++) {
1027 for (int oo = 0; oo != max_face_order; oo++) {
1028
1029 // face-edge
1030 if (oo < AinsworthOrderHooks::broken_nbfacetri_edge_hdiv(face_order))
1031 for (int dd = NBFACETRI_AINSWORTH_EDGE_HDIV(oo);
1032 dd != NBFACETRI_AINSWORTH_EDGE_HDIV(oo + 1); dd++) {
1033 for (int ee = 0; ee != 3; ++ee) {
1034 t_base(i) = t_base_f_e[ff * 3 + ee](i);
1035 ++t_base;
1036 ++t_base_f_e[ff * 3 + ee];
1037 }
1038 for (int ee = 0; ee != 3; ++ee) {
1039 t_diff_base(i, j) = t_diff_base_f_e[ff * 3 + ee](i, j);
1040 ++t_diff_base;
1041 ++t_diff_base_f_e[ff * 3 + ee];
1042 }
1043 }
1044
1045 // face-face
1046 if (oo < AinsworthOrderHooks::broken_nbfacetri_face_hdiv(face_order))
1047 for (int dd = NBFACETRI_AINSWORTH_FACE_HDIV(oo);
1048 dd != NBFACETRI_AINSWORTH_FACE_HDIV(oo + 1); dd++) {
1049 t_base(i) = t_base_f_f[ff](i);
1050 ++t_base;
1051 ++t_base_f_f[ff];
1052 t_diff_base(i, j) = t_diff_base_f_f[ff](i, j);
1053 ++t_diff_base;
1054 ++t_diff_base_f_f[ff];
1055 }
1056 }
1057 }
1058 }
1059 }
1060
1061 // volume
1062 int volume_dofs =
1063 6 * NBVOLUMETET_AINSWORTH_EDGE_HDIV(
1064 AinsworthOrderHooks::broken_nbvolumetet_edge_hdiv(volume_order)) +
1065 4 * NBVOLUMETET_AINSWORTH_FACE_HDIV(
1066 AinsworthOrderHooks::broken_nbvolumetet_face_hdiv(volume_order)) +
1067 NBVOLUMETET_AINSWORTH_VOLUME_HDIV(
1068 AinsworthOrderHooks::broken_nbvolumetet_volume_hdiv(volume_order));
1069 data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 3 * volume_dofs,
1070 false);
1071 data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
1072 9 * volume_dofs, false);
1073 if (volume_dofs) {
1074 double *base_ptr =
1075 &*data.dataOnEntities[MBTET][0].getN(base).data().begin();
1076 double *diff_base_ptr =
1077 &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin();
1078 auto t_base = getFTensor1FromPtr<3>(base_ptr);
1079 auto t_diff_base = getFTensor2HVecFromPtr<3, 3>(diff_base_ptr);
1080
1081 // volume-edge
1082 using Tensor1Ptr3 =
1083 decltype(getFTensor1FromPtr<3>(static_cast<double *>(nullptr)));
1084 Tensor1Ptr3 t_base_v_e[] = {
1085 getFTensor1FromPtr<3>(get_data_ptr(N_volume_edge[0])),
1086 getFTensor1FromPtr<3>(get_data_ptr(N_volume_edge[1])),
1087 getFTensor1FromPtr<3>(get_data_ptr(N_volume_edge[2])),
1088 getFTensor1FromPtr<3>(get_data_ptr(N_volume_edge[3])),
1089 getFTensor1FromPtr<3>(get_data_ptr(N_volume_edge[4])),
1090 getFTensor1FromPtr<3>(get_data_ptr(N_volume_edge[5]))};
1091 FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_v_e[] = {
1092 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_volume_edge[0])),
1093 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_volume_edge[1])),
1094 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_volume_edge[2])),
1095 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_volume_edge[3])),
1096 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_volume_edge[4])),
1097 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_volume_edge[5]))};
1098
1099 // volume-faces
1100 Tensor1Ptr3 t_base_v_f[] = {
1101 getFTensor1FromPtr<3>(get_data_ptr(N_volume_face[0])),
1102 getFTensor1FromPtr<3>(get_data_ptr(N_volume_face[1])),
1103 getFTensor1FromPtr<3>(get_data_ptr(N_volume_face[2])),
1104 getFTensor1FromPtr<3>(get_data_ptr(N_volume_face[3]))};
1105 FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_v_f[] = {
1106 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_volume_face[0])),
1107 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_volume_face[1])),
1108 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_volume_face[2])),
1109 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_volume_face[3]))};
1110
1111 // volume-bubble
1112 base_ptr = get_data_ptr(N_volume_bubble);
1113 diff_base_ptr = get_data_ptr(diffN_volume_bubble);
1114 auto t_base_v = getFTensor1FromPtr<3>(base_ptr);
1115 auto t_diff_base_v = getFTensor2HVecFromPtr<3, 3>(diff_base_ptr);
1116
1117 auto max_volume_order = std::max(
1118 volume_order,
1119 AinsworthOrderHooks::broken_nbvolumetet_edge_hdiv(volume_order));
1120 max_volume_order = std::max(
1121 max_volume_order,
1122 AinsworthOrderHooks::broken_nbvolumetet_face_hdiv(volume_order));
1123 max_volume_order = std::max(
1124 max_volume_order,
1125 AinsworthOrderHooks::broken_nbvolumetet_volume_hdiv(volume_order));
1126
1127 for (int gg = 0; gg != nb_gauss_pts; gg++) {
1128 for (int oo = 0; oo < max_volume_order; oo++) {
1129
1130 // volume-edge
1131 if (oo <
1132 AinsworthOrderHooks::broken_nbvolumetet_edge_hdiv(volume_order))
1133 for (int dd = NBVOLUMETET_AINSWORTH_EDGE_HDIV(oo);
1134 dd != NBVOLUMETET_AINSWORTH_EDGE_HDIV(oo + 1); dd++) {
1135 for (int ee = 0; ee < 6; ee++) {
1136 t_base(i) = t_base_v_e[ee](i);
1137 ++t_base;
1138 ++t_base_v_e[ee];
1139 t_diff_base(i, j) = t_diff_base_v_e[ee](i, j);
1140 ++t_diff_base;
1141 ++t_diff_base_v_e[ee];
1142 }
1143 }
1144
1145 // volume-face
1146 if (oo <
1147 AinsworthOrderHooks::broken_nbvolumetet_face_hdiv(volume_order))
1148 for (int dd = NBVOLUMETET_AINSWORTH_FACE_HDIV(oo);
1149 dd < NBVOLUMETET_AINSWORTH_FACE_HDIV(oo + 1); dd++) {
1150 for (int ff = 0; ff < 4; ff++) {
1151 t_base(i) = t_base_v_f[ff](i);
1152 ++t_base;
1153 ++t_base_v_f[ff];
1154 t_diff_base(i, j) = t_diff_base_v_f[ff](i, j);
1155 ++t_diff_base;
1156 ++t_diff_base_v_f[ff];
1157 }
1158 }
1159
1160 // volume-bubble
1161 if (oo <
1162 AinsworthOrderHooks::broken_nbvolumetet_volume_hdiv(volume_order))
1163 for (int dd = NBVOLUMETET_AINSWORTH_VOLUME_HDIV(oo);
1164 dd < NBVOLUMETET_AINSWORTH_VOLUME_HDIV(oo + 1); dd++) {
1165 t_base(i) = t_base_v(i);
1166 ++t_base;
1167 ++t_base_v;
1168 t_diff_base(i, j) = t_diff_base_v(i, j);
1169 ++t_diff_base;
1170 ++t_diff_base_v;
1171 }
1172 }
1173 }
1174 }
1175
1176 MoFEMFunctionReturn(0);
1177}
1178
1179MoFEMErrorCode
1180TetPolynomialBase::getValueHdivAinsworthBrokenBase(MatrixDouble &pts) {
1181 MoFEMFunctionBegin;
1182
1183 // Set shape functions into data structure Shape functions has to be put
1184 // in arrays in order which guarantee hierarchical series of degrees of
1185 // freedom, i.e. in other words dofs form sub-entities has to be group
1186 // by order.
1187
1188 FieldApproximationBase base = cTx->bAse;
1189 EntitiesFieldData &data = cTx->dAta;
1190
1191 FTENSOR_INDEX(3, i);
1192 FTENSOR_INDEX(3, j);
1193
1194 int volume_order = data.dataOnEntities[MBTET][0].getOrder();
1195 int nb_gauss_pts = pts.size2();
1196 auto get_data_ptr = [](MatrixDouble &m) -> double * {
1197 return (m.size1() && m.size2()) ? &*m.data().begin() : nullptr;
1198 };
1199 int nb_dofs_face =
1200 3 * NBFACETRI_AINSWORTH_EDGE_HDIV(
1201 AinsworthOrderHooks::broken_nbfacetri_edge_hdiv(volume_order)) +
1202 NBFACETRI_AINSWORTH_FACE_HDIV(
1203 AinsworthOrderHooks::broken_nbfacetri_face_hdiv(volume_order));
1204 int nb_dofs_volume =
1205 6 * NBVOLUMETET_AINSWORTH_EDGE_HDIV(
1206 AinsworthOrderHooks::broken_nbvolumetet_edge_hdiv(volume_order)) +
1207 4 * NBVOLUMETET_AINSWORTH_FACE_HDIV(
1208 AinsworthOrderHooks::broken_nbvolumetet_face_hdiv(volume_order)) +
1209 NBVOLUMETET_AINSWORTH_VOLUME_HDIV(
1210 AinsworthOrderHooks::broken_nbvolumetet_volume_hdiv(volume_order));
1211
1212 int nb_dofs = 4 * nb_dofs_face + nb_dofs_volume;
1213 data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 3 * nb_dofs,
1214 false);
1215 data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts, 9 * nb_dofs,
1216 false);
1217 if (nb_dofs == 0)
1218 MoFEMFunctionReturnHot(0);
1219
1220 auto get_interior_cache = [this](auto base) -> TetBaseCache::BaseCacheMI * {
1221 if (vPtr) {
1222 auto it = TetBaseCache::hdivBrokenBaseInterior[base].find(vPtr);
1223 if (it != TetBaseCache::hdivBrokenBaseInterior[base].end()) {
1224 return &it->second;
1225 }
1226 }
1227 return nullptr;
1228 };
1229
1230 auto interior_cache_ptr = get_interior_cache(base);
1231
1232 if (interior_cache_ptr) {
1233 auto it =
1234 interior_cache_ptr->find(boost::make_tuple(volume_order, nb_gauss_pts));
1235 if (it != interior_cache_ptr->end()) {
1236 noalias(data.dataOnEntities[MBTET][0].getN(base)) = it->N;
1237 noalias(data.dataOnEntities[MBTET][0].getDiffN(base)) = it->diffN;
1238 MoFEMFunctionReturnHot(0);
1239 }
1240 }
1241
1242 std::array<int, 4 * 3> faces_nodes = {0, 1, 3, 1, 2, 3, 0, 3, 2, 0, 2, 1};
1243 std::array<int, 4> faces_order{volume_order, volume_order, volume_order,
1244 volume_order};
1245 CHKERR getValueHdivAinsworthBaseImpl(
1246 pts, data.dataOnEntities[MBVERTEX][0].getN(base),
1247 data.dataOnEntities[MBVERTEX][0].getDiffN(base), volume_order,
1248 faces_order, faces_nodes,
1249
1250 AinsworthOrderHooks::broken_nbfacetri_edge_hdiv,
1251 AinsworthOrderHooks::broken_nbfacetri_face_hdiv,
1252 AinsworthOrderHooks::broken_nbvolumetet_edge_hdiv,
1253 AinsworthOrderHooks::broken_nbvolumetet_face_hdiv,
1254 AinsworthOrderHooks::broken_nbvolumetet_volume_hdiv
1255
1256 );
1257
1258 auto *base_ptr = &*data.dataOnEntities[MBTET][0].getN(base).data().begin();
1259 auto *diff_base_ptr =
1260 &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin();
1261 auto t_base = getFTensor1FromPtr<3>(base_ptr);
1262 auto t_diff_base = getFTensor2HVecFromPtr<3, 3>(diff_base_ptr);
1263
1264 // face-edge
1265 using Tensor1Ptr3 =
1266 decltype(getFTensor1FromPtr<3>(static_cast<double *>(nullptr)));
1267 Tensor1Ptr3 t_base_f_e[] = {
1268 getFTensor1FromPtr<3>(get_data_ptr(N_face_edge(0, 0))),
1269 getFTensor1FromPtr<3>(get_data_ptr(N_face_edge(0, 1))),
1270 getFTensor1FromPtr<3>(get_data_ptr(N_face_edge(0, 2))),
1271 getFTensor1FromPtr<3>(get_data_ptr(N_face_edge(1, 0))),
1272 getFTensor1FromPtr<3>(get_data_ptr(N_face_edge(1, 1))),
1273 getFTensor1FromPtr<3>(get_data_ptr(N_face_edge(1, 2))),
1274 getFTensor1FromPtr<3>(get_data_ptr(N_face_edge(2, 0))),
1275 getFTensor1FromPtr<3>(get_data_ptr(N_face_edge(2, 1))),
1276 getFTensor1FromPtr<3>(get_data_ptr(N_face_edge(2, 2))),
1277 getFTensor1FromPtr<3>(get_data_ptr(N_face_edge(3, 0))),
1278 getFTensor1FromPtr<3>(get_data_ptr(N_face_edge(3, 1))),
1279 getFTensor1FromPtr<3>(get_data_ptr(N_face_edge(3, 2)))};
1280 FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_f_e[] = {
1281 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_face_edge(0, 0))),
1282 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_face_edge(0, 1))),
1283 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_face_edge(0, 2))),
1284 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_face_edge(1, 0))),
1285 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_face_edge(1, 1))),
1286 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_face_edge(1, 2))),
1287 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_face_edge(2, 0))),
1288 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_face_edge(2, 1))),
1289 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_face_edge(2, 2))),
1290 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_face_edge(3, 0))),
1291 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_face_edge(3, 1))),
1292 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_face_edge(3, 2)))};
1293
1294 // face-face
1295 Tensor1Ptr3 t_base_f_f[] = {
1296 getFTensor1FromPtr<3>(get_data_ptr(N_face_bubble[0])),
1297 getFTensor1FromPtr<3>(get_data_ptr(N_face_bubble[1])),
1298 getFTensor1FromPtr<3>(get_data_ptr(N_face_bubble[2])),
1299 getFTensor1FromPtr<3>(get_data_ptr(N_face_bubble[3]))};
1300 FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_f_f[] = {
1301 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_face_bubble[0])),
1302 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_face_bubble[1])),
1303 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_face_bubble[2])),
1304 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_face_bubble[3]))};
1305
1306 // volume-edge
1307 Tensor1Ptr3 t_base_v_e[] = {
1308 getFTensor1FromPtr<3>(get_data_ptr(N_volume_edge[0])),
1309 getFTensor1FromPtr<3>(get_data_ptr(N_volume_edge[1])),
1310 getFTensor1FromPtr<3>(get_data_ptr(N_volume_edge[2])),
1311 getFTensor1FromPtr<3>(get_data_ptr(N_volume_edge[3])),
1312 getFTensor1FromPtr<3>(get_data_ptr(N_volume_edge[4])),
1313 getFTensor1FromPtr<3>(get_data_ptr(N_volume_edge[5]))};
1314 FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_v_e[] = {
1315 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_volume_edge[0])),
1316 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_volume_edge[1])),
1317 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_volume_edge[2])),
1318 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_volume_edge[3])),
1319 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_volume_edge[4])),
1320 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_volume_edge[5]))};
1321
1322 // volume-faces
1323 Tensor1Ptr3 t_base_v_f[] = {
1324 getFTensor1FromPtr<3>(get_data_ptr(N_volume_face[0])),
1325 getFTensor1FromPtr<3>(get_data_ptr(N_volume_face[1])),
1326 getFTensor1FromPtr<3>(get_data_ptr(N_volume_face[2])),
1327 getFTensor1FromPtr<3>(get_data_ptr(N_volume_face[3]))};
1328 FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_v_f[] = {
1329 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_volume_face[0])),
1330 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_volume_face[1])),
1331 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_volume_face[2])),
1332 getFTensor2HVecFromPtr<3, 3>(get_data_ptr(diffN_volume_face[3]))};
1333
1334 // volume-bubble
1335 auto *base_vol_ptr = get_data_ptr(N_volume_bubble);
1336 auto *diff_base_vol_ptr = get_data_ptr(diffN_volume_bubble);
1337 auto t_base_v = getFTensor1FromPtr<3>(base_vol_ptr);
1338 auto t_diff_base_v = getFTensor2HVecFromPtr<3, 3>(diff_base_vol_ptr);
1339
1340 int count_dofs = 0;
1341 int count_dofs_face = 0;
1342 int count_dofs_volume = 0;
1343
1344 auto max_volume_order =
1345 std::max(volume_order,
1346 AinsworthOrderHooks::broken_nbfacetri_edge_hdiv(volume_order));
1347 max_volume_order =
1348 std::max(max_volume_order,
1349 AinsworthOrderHooks::broken_nbfacetri_face_hdiv(volume_order));
1350 max_volume_order =
1351 std::max(max_volume_order,
1352 AinsworthOrderHooks::broken_nbvolumetet_edge_hdiv(volume_order));
1353 max_volume_order =
1354 std::max(max_volume_order,
1355 AinsworthOrderHooks::broken_nbvolumetet_face_hdiv(volume_order));
1356 max_volume_order = std::max(
1357 max_volume_order,
1358 AinsworthOrderHooks::broken_nbvolumetet_volume_hdiv(volume_order));
1359
1360 for (int gg = 0; gg != nb_gauss_pts; gg++) {
1361 for (int oo = 0; oo < max_volume_order; oo++) {
1362
1363 // faces-edge (((P) > 0) ? (P) : 0)
1364 if (oo < AinsworthOrderHooks::broken_nbfacetri_edge_hdiv(volume_order))
1365 for (int dd = NBFACETRI_AINSWORTH_EDGE_HDIV(oo);
1366 dd != NBFACETRI_AINSWORTH_EDGE_HDIV(oo + 1); dd++) {
1367 for (auto ff = 0; ff != 4; ++ff) {
1368 for (int ee = 0; ee != 3; ++ee) {
1369 t_base(i) = t_base_f_e[ff * 3 + ee](i);
1370 ++t_base;
1371 ++t_base_f_e[ff * 3 + ee];
1372 ++count_dofs;
1373 ++count_dofs_face;
1374 }
1375 for (int ee = 0; ee != 3; ++ee) {
1376 t_diff_base(i, j) = t_diff_base_f_e[ff * 3 + ee](i, j);
1377 ++t_diff_base;
1378 ++t_diff_base_f_e[ff * 3 + ee];
1379 }
1380 }
1381 }
1382
1383 // face-face (P - 1) * (P - 2) / 2
1384 if (oo < AinsworthOrderHooks::broken_nbfacetri_face_hdiv(volume_order))
1385 for (int dd = NBFACETRI_AINSWORTH_FACE_HDIV(oo);
1386 dd != NBFACETRI_AINSWORTH_FACE_HDIV(oo + 1); dd++) {
1387 for (auto ff = 0; ff != 4; ++ff) {
1388 t_base(i) = t_base_f_f[ff](i);
1389 ++t_base;
1390 ++t_base_f_f[ff];
1391 t_diff_base(i, j) = t_diff_base_f_f[ff](i, j);
1392 ++t_diff_base;
1393 ++t_diff_base_f_f[ff];
1394 ++count_dofs;
1395 ++count_dofs_face;
1396 }
1397 }
1398
1399 // volume-edge (P - 1)
1400 if (oo < AinsworthOrderHooks::broken_nbvolumetet_edge_hdiv(volume_order))
1401 for (int dd = NBVOLUMETET_AINSWORTH_EDGE_HDIV(oo);
1402 dd != NBVOLUMETET_AINSWORTH_EDGE_HDIV(oo + 1); dd++) {
1403 for (int ee = 0; ee < 6; ++ee) {
1404 t_base(i) = t_base_v_e[ee](i);
1405 ++t_base;
1406 ++t_base_v_e[ee];
1407 t_diff_base(i, j) = t_diff_base_v_e[ee](i, j);
1408 ++t_diff_base;
1409 ++t_diff_base_v_e[ee];
1410 ++count_dofs;
1411 ++count_dofs_volume;
1412 }
1413 }
1414
1415 // volume-face (P - 1) * (P - 2)
1416 if (oo < AinsworthOrderHooks::broken_nbvolumetet_face_hdiv(volume_order))
1417 for (int dd = NBVOLUMETET_AINSWORTH_FACE_HDIV(oo);
1418 dd != NBVOLUMETET_AINSWORTH_FACE_HDIV(oo + 1); dd++) {
1419 for (int ff = 0; ff < 4; ff++) {
1420 t_base(i) = t_base_v_f[ff](i);
1421 ++t_base;
1422 ++t_base_v_f[ff];
1423 t_diff_base(i, j) = t_diff_base_v_f[ff](i, j);
1424 ++t_diff_base;
1425 ++t_diff_base_v_f[ff];
1426 ++count_dofs;
1427 ++count_dofs_volume;
1428 }
1429 }
1430
1431 // volume-bubble (P - 3) * (P - 2) * (P - 1) / 2
1432 if (oo <
1433 AinsworthOrderHooks::broken_nbvolumetet_volume_hdiv(volume_order))
1434 for (int dd = NBVOLUMETET_AINSWORTH_VOLUME_HDIV(oo);
1435 dd != NBVOLUMETET_AINSWORTH_VOLUME_HDIV(oo + 1); dd++) {
1436 t_base(i) = t_base_v(i);
1437 ++t_base;
1438 ++t_base_v;
1439 t_diff_base(i, j) = t_diff_base_v(i, j);
1440 ++t_diff_base;
1441 ++t_diff_base_v;
1442 ++count_dofs;
1443 ++count_dofs_volume;
1444 }
1445 }
1446 }
1447
1448#ifndef NDEBUG
1449 if (nb_dofs != count_dofs / nb_gauss_pts) {
1450 MOFEM_LOG_CHANNEL("SELF");
1451 MOFEM_LOG("SELF", Sev::error) << "Nb dofs face: " << 4 * nb_dofs_face
1452 << " -> " << count_dofs_face / nb_gauss_pts;
1453 MOFEM_LOG("SELF", Sev::error) << "Nb dofs volume: " << nb_dofs_volume
1454 << " -> " << count_dofs_volume / nb_gauss_pts;
1455 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1456 "Number of dofs %d is different than expected %d",
1457 count_dofs / nb_gauss_pts, nb_dofs);
1458 }
1459#endif // NDEBUG
1460
1461 if (interior_cache_ptr) {
1462 auto p = interior_cache_ptr->emplace(
1463 TetBaseCache::BaseCacheItem{volume_order, nb_gauss_pts});
1464 p.first->N = data.dataOnEntities[MBTET][0].getN(base);
1465 p.first->diffN = data.dataOnEntities[MBTET][0].getDiffN(base);
1466 }
1467
1468 MoFEMFunctionReturn(0);
1469}
1470
1471MoFEMErrorCode TetPolynomialBase::getValueHdivDemkowiczBase(MatrixDouble &pts) {
1472 MoFEMFunctionBegin;
1473
1474 EntitiesFieldData &data = cTx->dAta;
1475 const FieldApproximationBase base = cTx->bAse;
1476 if (base != DEMKOWICZ_JACOBI_BASE) {
1477 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1478 "This should be used only with DEMKOWICZ_JACOBI_BASE "
1479 "but base is %s",
1480 ApproximationBaseNames[base]);
1481 }
1482 int nb_gauss_pts = pts.size2();
1483
1484 int volume_order = data.dataOnEntities[MBTET][0].getOrder();
1485
1486 int p_f[4];
1487 double *phi_f[4];
1488 double *diff_phi_f[4];
1489
1490 auto get_face_cache_ptr = [this]() -> TetBaseCache::HDivBaseFaceCacheMI * {
1491 if (vPtr) {
1492 auto it = TetBaseCache::hDivBaseFace[DEMKOWICZ_JACOBI_BASE].find(vPtr);
1493 if (it != TetBaseCache::hDivBaseFace[DEMKOWICZ_JACOBI_BASE].end()) {
1494 return &it->second;
1495 }
1496 }
1497 return nullptr;
1498 };
1499
1500 auto face_cache_ptr = get_face_cache_ptr();
1501
1502 // Calculate base function on tet faces
1503 for (int ff = 0; ff != 4; ff++) {
1504 int face_order = data.dataOnEntities[MBTRI][ff].getOrder();
1505 int order = volume_order > face_order ? volume_order : face_order;
1506 data.dataOnEntities[MBTRI][ff].getN(base).resize(
1507 nb_gauss_pts, 3 * NBFACETRI_DEMKOWICZ_HDIV(order), false);
1508 data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(
1509 nb_gauss_pts, 9 * NBFACETRI_DEMKOWICZ_HDIV(order), false);
1510 if (NBFACETRI_DEMKOWICZ_HDIV(order) == 0)
1511 continue;
1512
1513 if (face_cache_ptr) {
1514 auto it = face_cache_ptr->find(boost::make_tuple(
1515
1516 face_order, nb_gauss_pts,
1517
1518 data.facesNodes(ff, 0), data.facesNodes(ff, 1), data.facesNodes(ff, 2)
1519
1520 ));
1521 if (it != face_cache_ptr->end()) {
1522 noalias(data.dataOnEntities[MBTRI][ff].getN(base)) = it->N;
1523 noalias(data.dataOnEntities[MBTRI][ff].getDiffN(base)) = it->diffN;
1524 continue;
1525 }
1526 }
1527
1528 p_f[ff] = order;
1529 phi_f[ff] = &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
1530 diff_phi_f[ff] =
1531 &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
1532
1533 CHKERR Hdiv_Demkowicz_Face_MBTET_ON_FACE(
1534 &data.facesNodes(ff, 0), order,
1535 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1536 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1537 phi_f[ff], diff_phi_f[ff], nb_gauss_pts, 4);
1538 if (face_cache_ptr) {
1539 auto p = face_cache_ptr->emplace(TetBaseCache::HDivBaseCacheItem{
1540 face_order, nb_gauss_pts, data.facesNodes(ff, 0),
1541 data.facesNodes(ff, 1), data.facesNodes(ff, 2)});
1542 p.first->N = data.dataOnEntities[MBTRI][ff].getN(base);
1543 p.first->diffN = data.dataOnEntities[MBTRI][ff].getDiffN(base);
1544 }
1545 }
1546
1547 auto get_interior_cache = [this]() -> TetBaseCache::BaseCacheMI * {
1548 if (vPtr) {
1549 auto it =
1550 TetBaseCache::hdivBaseInterior[DEMKOWICZ_JACOBI_BASE].find(vPtr);
1551 if (it != TetBaseCache::hdivBaseInterior[DEMKOWICZ_JACOBI_BASE].end()) {
1552 return &it->second;
1553 }
1554 }
1555 return nullptr;
1556 };
1557
1558 auto interior_cache_ptr = get_interior_cache();
1559
1560 // Calculate base functions in tet interior
1561 if (NBVOLUMETET_DEMKOWICZ_HDIV(volume_order) > 0) {
1562 data.dataOnEntities[MBTET][0].getN(base).resize(
1563 nb_gauss_pts, 3 * NBVOLUMETET_DEMKOWICZ_HDIV(volume_order), false);
1564 data.dataOnEntities[MBTET][0].getDiffN(base).resize(
1565 nb_gauss_pts, 9 * NBVOLUMETET_DEMKOWICZ_HDIV(volume_order), false);
1566
1567 for (int v = 0; v != 1; ++v) {
1568 if (interior_cache_ptr) {
1569 auto it = interior_cache_ptr->find(
1570 boost::make_tuple(volume_order, nb_gauss_pts));
1571 if (it != interior_cache_ptr->end()) {
1572 noalias(data.dataOnEntities[MBTET][0].getN(base)) = it->N;
1573 noalias(data.dataOnEntities[MBTET][0].getDiffN(base)) = it->diffN;
1574 continue;
1575 }
1576 }
1577
1578 double *phi_v = &*data.dataOnEntities[MBTET][0].getN(base).data().begin();
1579 double *diff_phi_v =
1580 &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin();
1581
1582 CHKERR Hdiv_Demkowicz_Interior_MBTET(
1583 volume_order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
1584 &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), p_f, phi_f,
1585 diff_phi_f, phi_v, diff_phi_v, nb_gauss_pts);
1586 if (interior_cache_ptr) {
1587 auto p = interior_cache_ptr->emplace(
1588 TetBaseCache::BaseCacheItem{volume_order, nb_gauss_pts});
1589 p.first->N = data.dataOnEntities[MBTET][0].getN(base);
1590 p.first->diffN = data.dataOnEntities[MBTET][0].getDiffN(base);
1591 }
1592 }
1593 }
1594
1595 // Set size of face base correctly
1596 for (int ff = 0; ff != 4; ff++) {
1597 int face_order = data.dataOnEntities[MBTRI][ff].getOrder();
1598 data.dataOnEntities[MBTRI][ff].getN(base).resize(
1599 nb_gauss_pts, 3 * NBFACETRI_DEMKOWICZ_HDIV(face_order), true);
1600 data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(
1601 nb_gauss_pts, 9 * NBFACETRI_DEMKOWICZ_HDIV(face_order), true);
1602 }
1603
1604 MoFEMFunctionReturn(0);
1605}
1606
1607MoFEMErrorCode
1608TetPolynomialBase::getValueHdivDemkowiczBrokenBase(MatrixDouble &pts) {
1609 MoFEMFunctionBegin;
1610
1611 EntitiesFieldData &data = cTx->dAta;
1612 const FieldApproximationBase base = cTx->bAse;
1613 if (base != DEMKOWICZ_JACOBI_BASE) {
1614 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1615 "This should be used only with DEMKOWICZ_JACOBI_BASE "
1616 "but base is %s",
1617 ApproximationBaseNames[base]);
1618 }
1619 int nb_gauss_pts = pts.size2();
1620
1621 int volume_order = data.dataOnEntities[MBTET][0].getOrder();
1622 int nb_dofs_face = NBFACETRI_DEMKOWICZ_HDIV(volume_order);
1623 int nb_dofs_volume = NBVOLUMETET_DEMKOWICZ_HDIV(volume_order);
1624 int nb_dofs = 4 * nb_dofs_face + nb_dofs_volume;
1625 data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 3 * nb_dofs,
1626 false);
1627 data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts, 9 * nb_dofs,
1628 false);
1629 if (nb_dofs == 0)
1630 MoFEMFunctionReturnHot(0);
1631
1632 auto get_interior_cache = [this]() -> TetBaseCache::BaseCacheMI * {
1633 if (vPtr) {
1634 auto it =
1635 TetBaseCache::hdivBrokenBaseInterior[DEMKOWICZ_JACOBI_BASE].find(
1636 vPtr);
1637 if (it !=
1638 TetBaseCache::hdivBrokenBaseInterior[DEMKOWICZ_JACOBI_BASE].end()) {
1639 return &it->second;
1640 }
1641 }
1642 return nullptr;
1643 };
1644
1645 auto interior_cache_ptr = get_interior_cache();
1646
1647 if (interior_cache_ptr) {
1648 auto it =
1649 interior_cache_ptr->find(boost::make_tuple(volume_order, nb_gauss_pts));
1650 if (it != interior_cache_ptr->end()) {
1651 noalias(data.dataOnEntities[MBTET][0].getN(base)) = it->N;
1652 noalias(data.dataOnEntities[MBTET][0].getDiffN(base)) = it->diffN;
1653 MoFEMFunctionReturnHot(0);
1654 }
1655 }
1656
1657 std::array<MatrixDouble, 4> face_base_fun{
1658 MatrixDouble(nb_gauss_pts, 3 * nb_dofs_face),
1659 MatrixDouble(nb_gauss_pts, 3 * nb_dofs_face),
1660 MatrixDouble(nb_gauss_pts, 3 * nb_dofs_face),
1661 MatrixDouble(nb_gauss_pts, 3 * nb_dofs_face)};
1662 std::array<MatrixDouble, 4> face_diff_base{
1663 MatrixDouble(nb_gauss_pts, 9 * nb_dofs_face),
1664 MatrixDouble(nb_gauss_pts, 9 * nb_dofs_face),
1665 MatrixDouble(nb_gauss_pts, 9 * nb_dofs_face),
1666 MatrixDouble(nb_gauss_pts, 9 * nb_dofs_face)};
1667
1668 int faces_nodes[4][3] = {{0, 1, 3}, {1, 2, 3}, {0, 3, 2}, {0, 2, 1}};
1669
1670 std::array<int, 4> p_f{volume_order, volume_order, volume_order,
1671 volume_order};
1672 std::array<double *, 4> phi_f{
1673 &*face_base_fun[0].data().begin(), &*face_base_fun[1].data().begin(),
1674 &*face_base_fun[2].data().begin(), &*face_base_fun[3].data().begin()};
1675 std::array<double *, 4> diff_phi_f{
1676 &*face_diff_base[0].data().begin(), &*face_diff_base[1].data().begin(),
1677 &*face_diff_base[2].data().begin(), &*face_diff_base[3].data().begin()};
1678
1679 // Calculate base function on tet faces
1680 for (int ff = 0; ff != 4; ff++) {
1681 CHKERR Hdiv_Demkowicz_Face_MBTET_ON_FACE(
1682 // &data.facesNodes(ff, 0)
1683 faces_nodes[ff], p_f[ff],
1684 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1685 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1686 phi_f[ff], diff_phi_f[ff], nb_gauss_pts, 4);
1687 }
1688
1689 MatrixDouble vol_bases(nb_gauss_pts, 3 * nb_dofs_volume);
1690 MatrixDouble vol_diff_bases(nb_gauss_pts, 9 * nb_dofs_volume);
1691 auto *phi_v = &*vol_bases.data().begin();
1692 auto *diff_phi_v = &*vol_diff_bases.data().begin();
1693 CHKERR Hdiv_Demkowicz_Interior_MBTET(
1694 volume_order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
1695 &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), p_f.data(),
1696 phi_f.data(), diff_phi_f.data(), phi_v, diff_phi_v, nb_gauss_pts);
1697
1698 // faces
1699 using Tensor1Ptr3 =
1700 decltype(getFTensor1FromPtr<3>(static_cast<double *>(nullptr)));
1701 Tensor1Ptr3 t_base_v_f[] = {
1702 getFTensor1FromPtr<3>(phi_f[0]), getFTensor1FromPtr<3>(phi_f[1]),
1703 getFTensor1FromPtr<3>(phi_f[2]), getFTensor1FromPtr<3>(phi_f[3])};
1704 FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_v_f[] = {
1705 getFTensor2HVecFromPtr<3, 3>(diff_phi_f[0]),
1706 getFTensor2HVecFromPtr<3, 3>(diff_phi_f[1]),
1707 getFTensor2HVecFromPtr<3, 3>(diff_phi_f[2]),
1708 getFTensor2HVecFromPtr<3, 3>(diff_phi_f[3])};
1709
1710 // volumes
1711 auto t_base_v = getFTensor1FromPtr<3>(&*vol_bases.data().begin());
1712 FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_v =
1713 getFTensor2HVecFromPtr<3, 3>(&*vol_diff_bases.data().begin());
1714
1715 auto t_base = getFTensor1FromPtr<3>(
1716 &*data.dataOnEntities[MBTET][0].getN(base).data().begin());
1717 auto t_diff_base = getFTensor2HVecFromPtr<3, 3>(
1718 &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin());
1719
1720 FTENSOR_INDEX(3, i);
1721 FTENSOR_INDEX(3, j);
1722
1723 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
1724 for (int oo = 0; oo < volume_order; oo++) {
1725 // face
1726 for (auto dd = NBFACETRI_DEMKOWICZ_HDIV(oo);
1727 dd != NBFACETRI_DEMKOWICZ_HDIV(oo + 1); ++dd) {
1728 for (auto ff = 0; ff != 4; ++ff) {
1729 t_base(i) = t_base_v_f[ff](i);
1730 ++t_base;
1731 ++t_base_v_f[ff];
1732 t_diff_base(i, j) = t_diff_base_v_f[ff](i, j);
1733 ++t_diff_base;
1734 ++t_diff_base_v_f[ff];
1735 }
1736 }
1737 // volume
1738 for (auto dd = NBVOLUMETET_DEMKOWICZ_HDIV(oo);
1739 dd != NBVOLUMETET_DEMKOWICZ_HDIV(oo + 1); ++dd) {
1740 t_base(i) = t_base_v(i);
1741 ++t_base;
1742 ++t_base_v;
1743 t_diff_base(i, j) = t_diff_base_v(i, j);
1744 ++t_diff_base;
1745 ++t_diff_base_v;
1746 }
1747 }
1748 }
1749
1750 if (interior_cache_ptr) {
1751 auto p = interior_cache_ptr->emplace(
1752 TetBaseCache::BaseCacheItem{volume_order, nb_gauss_pts});
1753 p.first->N = data.dataOnEntities[MBTET][0].getN(base);
1754 p.first->diffN = data.dataOnEntities[MBTET][0].getDiffN(base);
1755 }
1756
1757 MoFEMFunctionReturn(0);
1758}
1759
1760MoFEMErrorCode TetPolynomialBase::getValueHdiv(MatrixDouble &pts) {
1761 MoFEMFunctionBeginHot;
1762
1763 switch (cTx->spaceContinuity) {
1764 case CONTINUOUS:
1765 switch (cTx->bAse) {
1766 case AINSWORTH_LEGENDRE_BASE:
1767 case AINSWORTH_LOBATTO_BASE:
1768 MoFEMFunctionReturnHot(getValueHdivAinsworthBase(pts));
1769 case DEMKOWICZ_JACOBI_BASE:
1770 MoFEMFunctionReturnHot(getValueHdivDemkowiczBase(pts));
1771 default:
1772 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
1773 }
1774 break;
1775 case DISCONTINUOUS:
1776 switch (cTx->bAse) {
1777 case AINSWORTH_LEGENDRE_BASE:
1778 case AINSWORTH_LOBATTO_BASE:
1779 MoFEMFunctionReturnHot(getValueHdivAinsworthBrokenBase(pts));
1780 case DEMKOWICZ_JACOBI_BASE:
1781 MoFEMFunctionReturnHot(getValueHdivDemkowiczBrokenBase(pts));
1782 default:
1783 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
1784 }
1785 break;
1786
1787 default:
1788 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown continuity");
1789 }
1790
1791 MoFEMFunctionReturnHot(0);
1792}
1793
1794MoFEMErrorCode
1795TetPolynomialBase::getValueHcurlAinsworthBase(MatrixDouble &pts) {
1796 MoFEMFunctionBegin;
1797
1798 EntitiesFieldData &data = cTx->dAta;
1799 const FieldApproximationBase base = cTx->bAse;
1800 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
1801 double *diffL, const int dim) =
1802 cTx->basePolynomialsType0;
1803
1804 int nb_gauss_pts = pts.size2();
1805
1806 // edges
1807 if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
1808 int sense[6], order[6];
1809 if (data.dataOnEntities[MBEDGE].size() != 6) {
1810 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1811 }
1812 double *hcurl_edge_n[6], *diff_hcurl_edge_n[6];
1813 for (int ee = 0; ee != 6; ee++) {
1814 if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
1815 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1816 "data inconsistency");
1817 }
1818 sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
1819 order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
1820 int nb_dofs =
1821 NBEDGE_AINSWORTH_HCURL(data.dataOnEntities[MBEDGE][ee].getOrder());
1822 data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts,
1823 3 * nb_dofs, false);
1824 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
1825 9 * nb_dofs, false);
1826 hcurl_edge_n[ee] =
1827 &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
1828 diff_hcurl_edge_n[ee] =
1829 &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
1830 }
1831 CHKERR Hcurl_Ainsworth_EdgeBaseFunctions_MBTET(
1832 sense, order,
1833 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1834 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1835 hcurl_edge_n, diff_hcurl_edge_n, nb_gauss_pts, base_polynomials);
1836 } else {
1837 for (int ee = 0; ee != 6; ee++) {
1838 data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0, false);
1839 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
1840 false);
1841 }
1842 }
1843
1844 // triangles
1845 if (data.spacesOnEntities[MBTRI].test(HCURL)) {
1846 int order[4];
1847 // faces
1848 if (data.dataOnEntities[MBTRI].size() != 4) {
1849 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1850 }
1851 double *hcurl_base_n[4], *diff_hcurl_base_n[4];
1852 for (int ff = 0; ff != 4; ff++) {
1853 if (data.dataOnEntities[MBTRI][ff].getSense() == 0) {
1854 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1855 "data inconsistency");
1856 }
1857 order[ff] = data.dataOnEntities[MBTRI][ff].getOrder();
1858 int nb_dofs = NBFACETRI_AINSWORTH_HCURL(order[ff]);
1859 data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts,
1860 3 * nb_dofs, false);
1861 data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
1862 9 * nb_dofs, false);
1863 hcurl_base_n[ff] =
1864 &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
1865 diff_hcurl_base_n[ff] =
1866 &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
1867 }
1868 if (data.facesNodes.size1() != 4) {
1869 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1870 }
1871 if (data.facesNodes.size2() != 3) {
1872 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1873 }
1874 CHKERR Hcurl_Ainsworth_FaceFunctions_MBTET(
1875 &*data.facesNodes.data().begin(), order,
1876 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1877 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1878 hcurl_base_n, diff_hcurl_base_n, nb_gauss_pts, base_polynomials);
1879 } else {
1880 for (int ff = 0; ff != 4; ff++) {
1881 data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, 0, false);
1882 data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts, 0,
1883 false);
1884 }
1885 }
1886
1887 if (data.spacesOnEntities[MBTET].test(HCURL)) {
1888
1889 // volume
1890 int order = data.dataOnEntities[MBTET][0].getOrder();
1891 int nb_vol_dofs = NBVOLUMETET_AINSWORTH_HCURL(order);
1892 data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts,
1893 3 * nb_vol_dofs, false);
1894 data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
1895 9 * nb_vol_dofs, false);
1896 CHKERR Hcurl_Ainsworth_VolumeFunctions_MBTET(
1897 data.dataOnEntities[MBTET][0].getOrder(),
1898 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1899 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1900 &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
1901 &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
1902 nb_gauss_pts, base_polynomials);
1903
1904 } else {
1905 data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 0, false);
1906 data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
1907 }
1908
1909 MoFEMFunctionReturn(0);
1910}
1911
1912MoFEMErrorCode
1913TetPolynomialBase::getValueHcurlDemkowiczBase(MatrixDouble &pts) {
1914 MoFEMFunctionBegin;
1915
1916 EntitiesFieldData &data = cTx->dAta;
1917 const FieldApproximationBase base = cTx->bAse;
1918 if (base != DEMKOWICZ_JACOBI_BASE) {
1919 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1920 "This should be used only with DEMKOWICZ_JACOBI_BASE "
1921 "but base is %s",
1922 ApproximationBaseNames[base]);
1923 }
1924
1925 int nb_gauss_pts = pts.size2();
1926
1927 // edges
1928 if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
1929 int sense[6], order[6];
1930 if (data.dataOnEntities[MBEDGE].size() != 6) {
1931 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1932 "wrong size of data structure, expected space for six edges "
1933 "but is %zu",
1934 data.dataOnEntities[MBEDGE].size());
1935 }
1936 double *hcurl_edge_n[6], *diff_hcurl_edge_n[6];
1937 for (int ee = 0; ee != 6; ee++) {
1938 if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
1939 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1940 "orintation of edges is not set");
1941 }
1942 sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
1943 order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
1944 int nb_dofs =
1945 NBEDGE_DEMKOWICZ_HCURL(data.dataOnEntities[MBEDGE][ee].getOrder());
1946 data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts,
1947 3 * nb_dofs, false);
1948 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
1949 9 * nb_dofs, false);
1950 hcurl_edge_n[ee] =
1951 &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
1952 diff_hcurl_edge_n[ee] =
1953 &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
1954 }
1955 CHKERR Hcurl_Demkowicz_EdgeBaseFunctions_MBTET(
1956 sense, order,
1957 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1958 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1959 hcurl_edge_n, diff_hcurl_edge_n, nb_gauss_pts);
1960 } else {
1961 // No DOFs on edges, resize base function matrices, indicating that no
1962 // dofs on them.
1963 for (int ee = 0; ee != 6; ee++) {
1964 data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0, false);
1965 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
1966 false);
1967 }
1968 }
1969
1970 // triangles
1971 if (data.spacesOnEntities[MBTRI].test(HCURL)) {
1972 int order[4];
1973 // faces
1974 if (data.dataOnEntities[MBTRI].size() != 4) {
1975 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1976 "data structure for storing face h-curl base have wrong size "
1977 "should be four but is %zu",
1978 data.dataOnEntities[MBTRI].size());
1979 }
1980 double *hcurl_base_n[4], *diff_hcurl_base_n[4];
1981 for (int ff = 0; ff != 4; ff++) {
1982 if (data.dataOnEntities[MBTRI][ff].getSense() == 0) {
1983 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1984 "orintation of face is not set");
1985 }
1986 order[ff] = data.dataOnEntities[MBTRI][ff].getOrder();
1987 int nb_dofs = NBFACETRI_DEMKOWICZ_HCURL(order[ff]);
1988 data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts,
1989 3 * nb_dofs, false);
1990 data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
1991 9 * nb_dofs, false);
1992 hcurl_base_n[ff] =
1993 &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
1994 diff_hcurl_base_n[ff] =
1995 &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
1996 }
1997 if (data.facesNodes.size1() != 4) {
1998 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1999 "data inconsistency, should be four faces");
2000 }
2001 if (data.facesNodes.size2() != 3) {
2002 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2003 "data inconsistency, should be three nodes on face");
2004 }
2005 CHKERR Hcurl_Demkowicz_FaceBaseFunctions_MBTET(
2006 &*data.facesNodes.data().begin(), order,
2007 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
2008 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
2009 hcurl_base_n, diff_hcurl_base_n, nb_gauss_pts);
2010 } else {
2011 // No DOFs on faces, resize base function matrices, indicating that no
2012 // dofs on them.
2013 for (int ff = 0; ff != 4; ff++) {
2014 data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, 0, false);
2015 data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts, 0,
2016 false);
2017 }
2018 }
2019
2020 if (data.spacesOnEntities[MBTET].test(HCURL)) {
2021 // volume
2022 int order = data.dataOnEntities[MBTET][0].getOrder();
2023 int nb_vol_dofs = NBVOLUMETET_DEMKOWICZ_HCURL(order);
2024 data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts,
2025 3 * nb_vol_dofs, false);
2026 data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
2027 9 * nb_vol_dofs, false);
2028 CHKERR Hcurl_Demkowicz_VolumeBaseFunctions_MBTET(
2029 data.dataOnEntities[MBTET][0].getOrder(),
2030 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
2031 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
2032 &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
2033 &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
2034 nb_gauss_pts);
2035 } else {
2036 // No DOFs on faces, resize base function matrices, indicating that no
2037 // dofs on them.
2038 data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 0, false);
2039 data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
2040 }
2041
2042 MoFEMFunctionReturn(0);
2043}
2044
2045MoFEMErrorCode TetPolynomialBase::getValueHcurl(MatrixDouble &pts) {
2046 MoFEMFunctionBegin;
2047
2048 switch (cTx->bAse) {
2049 case AINSWORTH_LEGENDRE_BASE:
2050 case AINSWORTH_LOBATTO_BASE:
2051 CHKERR getValueHcurlAinsworthBase(pts);
2052 break;
2053 case DEMKOWICZ_JACOBI_BASE:
2054 CHKERR getValueHcurlDemkowiczBase(pts);
2055 break;
2056 default:
2057 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
2058 }
2059
2060 MoFEMFunctionReturn(0);
2061}
2062
2063MoFEMErrorCode
2064TetPolynomialBase::getValue(MatrixDouble &pts,
2065 boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
2066 MoFEMFunctionBegin;
2067
2068 cTx = ctx_ptr->getInterface<EntPolynomialBaseCtx>();
2069
2070 int nb_gauss_pts = pts.size2();
2071 if (!nb_gauss_pts)
2072 MoFEMFunctionReturnHot(0);
2073
2074 if (pts.size1() < 3)
2075 SETERRQ(
2076 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2077 "Wrong dimension of pts, should be at least 3 rows with coordinates");
2078
2079 if (cTx->bAse != AINSWORTH_BERNSTEIN_BEZIER_BASE) {
2080 const FieldApproximationBase base = cTx->bAse;
2081 EntitiesFieldData &data = cTx->dAta;
2082 if (cTx->copyNodeBase == LASTBASE) {
2083 data.dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 4,
2084 false);
2085 CHKERR Tools::shapeFunMBTET(
2086 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
2087 &pts(0, 0), &pts(1, 0), &pts(2, 0), nb_gauss_pts);
2088 } else {
2089 data.dataOnEntities[MBVERTEX][0].getN(base) =
2090 data.dataOnEntities[MBVERTEX][0].getN(cTx->copyNodeBase);
2091 }
2092 if (data.dataOnEntities[MBVERTEX][0].getN(base).size1() !=
2093 (unsigned int)nb_gauss_pts) {
2094 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2095 "Base functions or nodes has wrong number of integration points "
2096 "for base %s",
2097 ApproximationBaseNames[base]);
2098 }
2099 data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(4, 3, false);
2100 std::copy(Tools::diffShapeFunMBTET.begin(), Tools::diffShapeFunMBTET.end(),
2101 data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin());
2102 }
2103
2104 switch (cTx->spaceContinuity) {
2105 case CONTINUOUS:
2106
2107 switch (cTx->sPace) {
2108 case H1:
2109 CHKERR getValueH1(pts);
2110 break;
2111 case HDIV:
2112 CHKERR getValueHdiv(pts);
2113 break;
2114 case HCURL:
2115 CHKERR getValueHcurl(pts);
2116 break;
2117 case L2:
2118 CHKERR getValueL2(pts);
2119 break;
2120 default:
2121 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown space %s",
2122 FieldSpaceNames[cTx->sPace]);
2123 }
2124 break;
2125
2126 case DISCONTINUOUS:
2127
2128 switch (cTx->sPace) {
2129 case HDIV:
2130 CHKERR getValueHdiv(pts);
2131 break;
2132 default:
2133 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown space %s",
2134 FieldSpaceNames[cTx->sPace]);
2135 }
2136 break;
2137
2138 default:
2139 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown continuity");
2140 }
2141
2142 MoFEMFunctionReturn(0);
2143}
2144
2145MoFEMErrorCode TetPolynomialBase::setDofsSideMap(
2146
2147 const FieldSpace space, const FieldContinuity continuity,
2148 const FieldApproximationBase base, DofsSideMap &dofs_side_map
2149
2150) {
2151 MoFEMFunctionBegin;
2152
2153 switch (continuity) {
2154 case DISCONTINUOUS:
2155
2156 switch (space) {
2157 case HDIV:
2158 CHKERR setDofsSideMapHdiv(continuity, base, dofs_side_map);
2159 break;
2160 default:
2161 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown space %s",
2162 FieldSpaceNames[space]);
2163 }
2164 break;
2165
2166 default:
2167 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
2168 "Unknown (or not implemented) continuity");
2169 }
2170
2171 MoFEMFunctionReturn(0);
2172}
2173
2174MoFEMErrorCode
2175TetPolynomialBase::setDofsSideMapHdiv(const FieldContinuity continuity,
2176 const FieldApproximationBase base,
2177 DofsSideMap &dofs_side_map) {
2178 MoFEMFunctionBegin;
2179
2180 // That has to be consistent with implementation of getValueHdiv for
2181 // particular base functions.
2182
2183 auto set_ainsworth = [&dofs_side_map]() {
2184 MoFEMFunctionBegin;
2185
2186 dofs_side_map.clear();
2187
2188 const auto map_order = Field::maxBrokenDofsOrder;
2189 const auto face_edge_order =
2190 AinsworthOrderHooks::broken_nbfacetri_edge_hdiv(map_order);
2191 const auto face_face_order =
2192 AinsworthOrderHooks::broken_nbfacetri_face_hdiv(map_order);
2193 const auto volume_edge_order =
2194 AinsworthOrderHooks::broken_nbvolumetet_edge_hdiv(map_order);
2195 const auto volume_face_order =
2196 AinsworthOrderHooks::broken_nbvolumetet_face_hdiv(map_order);
2197 const auto volume_bubble_order =
2198 AinsworthOrderHooks::broken_nbvolumetet_volume_hdiv(map_order);
2199
2200 auto max_order = std::max(map_order, face_edge_order);
2201 max_order = std::max(max_order, face_face_order);
2202 max_order = std::max(max_order, volume_edge_order);
2203 max_order = std::max(max_order, volume_face_order);
2204 max_order = std::max(max_order, volume_bubble_order);
2205
2206 int dof = 0;
2207 for (int oo = 0; oo < max_order; oo++) {
2208
2209 // faces-edge (((P) > 0) ? (P) : 0)
2210 if (oo < face_edge_order)
2211 for (int dd = NBFACETRI_AINSWORTH_EDGE_HDIV(oo);
2212 dd != NBFACETRI_AINSWORTH_EDGE_HDIV(oo + 1); dd++) {
2213 for (auto ff = 0; ff != 4; ++ff) {
2214 for (int ee = 0; ee != 3; ++ee) {
2215 dofs_side_map.insert(DofsSideMapData{MBTRI, ff, dof});
2216 ++dof;
2217 }
2218 }
2219 }
2220
2221 // face-face (P - 1) * (P - 2) / 2
2222 if (oo < face_face_order)
2223 for (int dd = NBFACETRI_AINSWORTH_FACE_HDIV(oo);
2224 dd != NBFACETRI_AINSWORTH_FACE_HDIV(oo + 1); dd++) {
2225 for (auto ff = 0; ff != 4; ++ff) {
2226 dofs_side_map.insert(DofsSideMapData{MBTRI, ff, dof});
2227 ++dof;
2228 }
2229 }
2230
2231 // volume-edge (P - 1)
2232 if (oo < volume_edge_order)
2233 for (int dd = NBVOLUMETET_AINSWORTH_EDGE_HDIV(oo);
2234 dd != NBVOLUMETET_AINSWORTH_EDGE_HDIV(oo + 1); dd++) {
2235 for (int ee = 0; ee < 6; ee++) {
2236 dofs_side_map.insert(DofsSideMapData{MBTET, 0, dof});
2237 ++dof;
2238 }
2239 }
2240
2241 // volume-face (P - 1) * (P - 2)
2242 if (oo < volume_face_order)
2243 for (int dd = NBVOLUMETET_AINSWORTH_FACE_HDIV(oo);
2244 dd < NBVOLUMETET_AINSWORTH_FACE_HDIV(oo + 1); dd++) {
2245 for (int ff = 0; ff < 4; ff++) {
2246 dofs_side_map.insert(DofsSideMapData{MBTET, 0, dof});
2247 ++dof;
2248 }
2249 }
2250
2251 // volume-bubble (P - 3) * (P - 2) * (P - 1) / 2
2252 if (oo < volume_bubble_order)
2253 for (int dd = NBVOLUMETET_AINSWORTH_VOLUME_HDIV(oo);
2254 dd < NBVOLUMETET_AINSWORTH_VOLUME_HDIV(oo + 1); dd++) {
2255 dofs_side_map.insert(DofsSideMapData{MBTET, 0, dof});
2256 ++dof;
2257 }
2258 }
2259
2260 MoFEMFunctionReturn(0);
2261 };
2262
2263 auto set_demkowicz = [&dofs_side_map]() {
2264 MoFEMFunctionBegin;
2265
2266 dofs_side_map.clear();
2267
2268 int dof = 0;
2269 for (int oo = 0; oo < Field::maxBrokenDofsOrder; oo++) {
2270
2271 // face
2272 for (auto dd = NBFACETRI_DEMKOWICZ_HDIV(oo);
2273 dd != NBFACETRI_DEMKOWICZ_HDIV(oo + 1); ++dd) {
2274 for (auto ff = 0; ff != 4; ++ff) {
2275 dofs_side_map.insert(DofsSideMapData{MBTRI, ff, dof});
2276 ++dof;
2277 }
2278 }
2279 // volume
2280 for (auto dd = NBVOLUMETET_DEMKOWICZ_HDIV(oo);
2281 dd != NBVOLUMETET_DEMKOWICZ_HDIV(oo + 1); ++dd) {
2282 dofs_side_map.insert(DofsSideMapData{MBTET, 0, dof});
2283 ++dof;
2284 }
2285 }
2286
2287 MoFEMFunctionReturn(0);
2288 };
2289
2290 switch (continuity) {
2291 case DISCONTINUOUS:
2292 switch (base) {
2293 case AINSWORTH_LEGENDRE_BASE:
2294 case AINSWORTH_LOBATTO_BASE:
2295 MoFEMFunctionReturnHot(set_ainsworth());
2296 case DEMKOWICZ_JACOBI_BASE:
2297 MoFEMFunctionReturnHot(set_demkowicz());
2298 default:
2299 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
2300 }
2301 break;
2302
2303 default:
2304 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
2305 "Unknown (or not implemented) continuity");
2306 }
2307
2308 MoFEMFunctionReturn(0);
2309}
2310
2311template <typename T>
2312auto tetCacheSwitch(const void *ptr, T &cache, std::string cache_name) {
2313 auto it = cache.find(ptr);
2314 if (it != cache.end()) {
2315 MOFEM_LOG_CHANNEL("WORLD");
2316 MOFEM_TAG_AND_LOG("WORLD", Sev::noisy, "TetPolynomialBase")
2317 << "Cache off " << cache_name << ": " << it->second.size();
2318 cache.erase(it);
2319 return false;
2320 } else {
2321 MOFEM_LOG_CHANNEL("WORLD");
2322 MOFEM_TAG_AND_LOG("WORLD", Sev::noisy, "TetPolynomialBase")
2323 << "Cache on " << cache_name;
2324 cache[ptr];
2325 return true;
2326 }
2327}
2328
2329template <>
2330bool TetPolynomialBase::switchCacheBaseFace<HDIV>(FieldApproximationBase base,
2331 void *ptr) {
2332 return tetCacheSwitch(ptr, TetBaseCache::hDivBaseFace[base],
2333 std::string("hDivBaseFace") +
2334 ApproximationBaseNames[base]);
2335}
2336
2337template <>
2338bool TetPolynomialBase::switchCacheBaseInterior<HDIV>(
2339 FieldApproximationBase base, void *ptr) {
2340 return tetCacheSwitch(ptr, TetBaseCache::hdivBaseInterior[base],
2341 std::string("hdivBaseInterior") +
2342 ApproximationBaseNames[base]);
2343}
2344
2345template <>
2346bool TetPolynomialBase::switchCacheBrokenBaseInterior<HDIV>(
2347 FieldApproximationBase base, void *ptr) {
2348 return tetCacheSwitch(ptr, TetBaseCache::hdivBrokenBaseInterior[base],
2349 std::string("hdivBrokenBaseInterior") +
2350 ApproximationBaseNames[base]);
2351}
2352
2353template <>
2354void TetPolynomialBase::switchCacheBaseOn<HDIV>(FieldApproximationBase base,
2355 std::vector<void *> v) {
2356 for (auto fe_ptr : v) {
2357 if (!TetPolynomialBase::switchCacheBaseFace<HDIV>(base, fe_ptr)) {
2358 TetPolynomialBase::switchCacheBaseFace<HDIV>(base, fe_ptr);
2359 }
2360 if (!TetPolynomialBase::switchCacheBaseInterior<HDIV>(base, fe_ptr)) {
2361 TetPolynomialBase::switchCacheBaseInterior<HDIV>(base, fe_ptr);
2362 }
2363 if (!TetPolynomialBase::switchCacheBrokenBaseInterior<HDIV>(base, fe_ptr)) {
2364 TetPolynomialBase::switchCacheBrokenBaseInterior<HDIV>(base, fe_ptr);
2365 }
2366 }
2367}
2368
2369template <>
2370void TetPolynomialBase::switchCacheBaseOff<HDIV>(FieldApproximationBase base,
2371 std::vector<void *> v) {
2372 for (auto fe_ptr : v) {
2373 if (TetPolynomialBase::switchCacheBaseFace<HDIV>(base, fe_ptr)) {
2374 TetPolynomialBase::switchCacheBaseFace<HDIV>(base, fe_ptr);
2375 }
2376 if (TetPolynomialBase::switchCacheBaseInterior<HDIV>(base, fe_ptr)) {
2377 TetPolynomialBase::switchCacheBaseInterior<HDIV>(base, fe_ptr);
2378 }
2379 if (TetPolynomialBase::switchCacheBrokenBaseInterior<HDIV>(base, fe_ptr)) {
2380 TetPolynomialBase::switchCacheBrokenBaseInterior<HDIV>(base, fe_ptr);
2381 }
2382 }
2383}
2384
2385template <>
2386void TetPolynomialBase::switchCacheBaseOn<HDIV>(std::vector<void *> v) {
2387 for (auto b = 0; b != LASTBASE; ++b) {
2388 switchCacheBaseOn<HDIV>(static_cast<FieldApproximationBase>(b), v);
2389 }
2390}
2391
2392template <>
2393void TetPolynomialBase::switchCacheBaseOff<HDIV>(std::vector<void *> v) {
2394 for (auto b = 0; b != LASTBASE; ++b) {
2395 switchCacheBaseOff<HDIV>(static_cast<FieldApproximationBase>(b), v);
2396 }
2397}
2398
2399template <>
2400bool TetPolynomialBase::switchCacheBaseInterior<L2>(FieldApproximationBase base,
2401 void *ptr) {
2402 return tetCacheSwitch(ptr, TetBaseCache::l2BaseInterior[base],
2403 std::string("hdivBaseInterior") +
2404 ApproximationBaseNames[base]);
2405}
2406
2407template <>
2408void TetPolynomialBase::switchCacheBaseOn<L2>(FieldApproximationBase base,
2409 std::vector<void *> v) {
2410 for (auto fe_ptr : v) {
2411 if (!TetPolynomialBase::switchCacheBaseInterior<L2>(base, fe_ptr)) {
2412 TetPolynomialBase::switchCacheBaseInterior<L2>(base, fe_ptr);
2413 }
2414 }
2415}
2416
2417template <>
2418void TetPolynomialBase::switchCacheBaseOff<L2>(FieldApproximationBase base,
2419 std::vector<void *> v) {
2420 for (auto fe_ptr : v) {
2421 if (TetPolynomialBase::switchCacheBaseInterior<L2>(base, fe_ptr)) {
2422 TetPolynomialBase::switchCacheBaseInterior<L2>(base, fe_ptr);
2423 }
2424 }
2425}
2426
2427template <>
2428void TetPolynomialBase::switchCacheBaseOn<L2>(std::vector<void *> v) {
2429 for (auto b = 0; b != LASTBASE; ++b) {
2430 switchCacheBaseOn<L2>(static_cast<FieldApproximationBase>(b), v);
2431 }
2432}
2433
2434template <>
2435void TetPolynomialBase::switchCacheBaseOff<L2>(std::vector<void *> v) {
2436 for (auto b = 0; b != LASTBASE; ++b) {
2437 switchCacheBaseOff<L2>(static_cast<FieldApproximationBase>(b), v);
2438 }
2439}
MOFEM_TAG_AND_LOG
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition LogManager.hpp:376
FTENSOR_INDEX
#define FTENSOR_INDEX(DIM, I)
Definition Templates.hpp:2561
tetCacheSwitch
auto tetCacheSwitch(const void *ptr, T &cache, std::string cache_name)
Definition TetPolynomialBase.cpp:2312
a
constexpr double a
Definition approx_sphere.cpp:30
FieldApproximationBase
FieldApproximationBase
approximation base
Definition definitions.h:58
LASTBASE
@ LASTBASE
Definition definitions.h:69
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
AINSWORTH_LOBATTO_BASE
@ AINSWORTH_LOBATTO_BASE
Definition definitions.h:62
NOBASE
@ NOBASE
Definition definitions.h:59
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition definitions.h:66
AINSWORTH_BERNSTEIN_BEZIER_BASE
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition definitions.h:64
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:463
FieldSpace
FieldSpace
approximation spaces
Definition definitions.h:82
L2
@ L2
field with C-1 continuity
Definition definitions.h:88
H1
@ H1
continuous field
Definition definitions.h:85
HCURL
@ HCURL
field with continuous tangents
Definition definitions.h:86
HDIV
@ HDIV
field with continuous normal traction
Definition definitions.h:87
FieldContinuity
FieldContinuity
Field continuity.
Definition definitions.h:99
CONTINUOUS
@ CONTINUOUS
Regular field.
Definition definitions.h:100
DISCONTINUOUS
@ DISCONTINUOUS
Broken continuity (No effect on L2 space)
Definition definitions.h:101
FieldSpaceNames
static const char *const FieldSpaceNames[]
Definition definitions.h:92
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
ApproximationBaseNames
static const char *const ApproximationBaseNames[]
Definition definitions.h:72
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:432
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:38
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition LogManager.hpp:316
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition LogManager.hpp:292
NBVOLUMETET_AINSWORTH_EDGE_HDIV
#define NBVOLUMETET_AINSWORTH_EDGE_HDIV(P)
Definition h1_hdiv_hcurl_l2.h:133
NBEDGE_DEMKOWICZ_HCURL
#define NBEDGE_DEMKOWICZ_HCURL(P)
Definition h1_hdiv_hcurl_l2.h:108
NBVOLUMETET_H1
#define NBVOLUMETET_H1(P)
Number of base functions on tetrahedron for H1 space.
Definition h1_hdiv_hcurl_l2.h:75
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_DEMKOWICZ_HDIV
#define NBVOLUMETET_DEMKOWICZ_HDIV(P)
Definition h1_hdiv_hcurl_l2.h:140
H1_EdgeShapeFunctions_MBTET
PetscErrorCode H1_EdgeShapeFunctions_MBTET(int *sense, int *p, double *N, double *diffN, double *edgeN[], double *diff_edgeN[], int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Definition h1.c:274
NBVOLUMETET_AINSWORTH_FACE_HDIV
#define NBVOLUMETET_AINSWORTH_FACE_HDIV(P)
Definition h1_hdiv_hcurl_l2.h:134
NBEDGE_H1
#define NBEDGE_H1(P)
Number of base function on edge for H1 space.
Definition h1_hdiv_hcurl_l2.h:55
NBFACETRI_DEMKOWICZ_HDIV
#define NBFACETRI_DEMKOWICZ_HDIV(P)
Definition h1_hdiv_hcurl_l2.h:139
L2_Ainsworth_ShapeFunctions_MBTET
PetscErrorCode L2_Ainsworth_ShapeFunctions_MBTET(int p, double *N, double *diffN, double *L2N, double *diff_L2N, int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Get base functions on tetrahedron for L2 space.
Definition l2.c:74
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_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
H1_VolumeShapeFunctions_MBTET
PetscErrorCode H1_VolumeShapeFunctions_MBTET(int p, double *N, double *diffN, double *volumeN, double *diff_volumeN, int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Definition h1.c:475
NBEDGE_AINSWORTH_HCURL
#define NBEDGE_AINSWORTH_HCURL(P)
Definition h1_hdiv_hcurl_l2.h:97
H1_FaceShapeFunctions_MBTET
PetscErrorCode H1_FaceShapeFunctions_MBTET(int *faces_nodes, int *p, double *N, double *diffN, double *faceN[], double *diff_faceN[], int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Definition h1.c:373
NBFACETRI_H1
#define NBFACETRI_H1(P)
Number of base function on triangle for H1 space.
Definition h1_hdiv_hcurl_l2.h:60
NBVOLUMETET_L2
#define NBVOLUMETET_L2(P)
Number of base functions on tetrahedron for L2 space.
Definition h1_hdiv_hcurl_l2.h:27
i
FTensor::Index< 'i', SPACE_DIM > i
Definition hcurl_divergence_operator_2d.cpp:27
lambda
static double lambda
Definition incompressible_elasticity.cpp:181
v
const double v
phase velocity of light in medium (cm/ns)
Definition initial_diffusion.cpp:40
n
const double n
refractive index of diffusive medium
Definition initial_diffusion.cpp:38
j
FTensor::Index< 'j', 3 > j
Definition matrix_function.cpp:21
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition Exceptions.hpp:56
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition Types.hpp:77
MoFEM::Types::MatrixInt
UBlasMatrix< int > MatrixInt
Definition Types.hpp:76
MoFEM
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
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:683
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 .
Definition Hdiv.cpp:97
MoFEM::getFTensor2HVecFromPtr< 3, 3 >
auto getFTensor2HVecFromPtr< 3, 3 >(double *ptr)
Definition Templates.hpp:641
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:829
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::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:450
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 .
Definition Hdiv.cpp:224
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_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:2402
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 .
Definition Hdiv.cpp:357
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:2475
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:2079
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::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 .
Definition Hdiv.cpp:556
field_name
constexpr auto field_name
Definition poisson_2d_homogeneous.cpp:13
g
constexpr double g
Definition shallow_wave.cpp:63
m
FTensor::Index< 'm', 3 > m
Definition shallow_wave.cpp:80
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
MoFEM::BaseFunction::DofsSideMap
multi_index_container< DofsSideMapData, indexed_by< ordered_non_unique< tag< TypeSide_mi_tag >, composite_key< DofsSideMapData, member< DofsSideMapData, EntityType, &DofsSideMapData::type >, member< DofsSideMapData, int, &DofsSideMapData::side > > >, ordered_unique< tag< EntDofIdx_mi_tag >, member< DofsSideMapData, int, &DofsSideMapData::dof > > > > DofsSideMap
Map entity stype and side to element/entity dof index.
Definition BaseFunction.hpp:71
MoFEM::BernsteinBezier::generateIndicesTriTet
static MoFEMErrorCode generateIndicesTriTet(const int N[], int *alpha[])
Definition BernsteinBezier.cpp:173
MoFEM::BernsteinBezier::baseFunctionsTet
static MoFEMErrorCode baseFunctionsTet(const int N, const int gdim, const int n_alpha, const int *alpha, const double *lambda, const double *grad_lambda, double *base, double *grad_base)
Definition BernsteinBezier.cpp:441
MoFEM::BernsteinBezier::generateIndicesEdgeTet
static MoFEMErrorCode generateIndicesEdgeTet(const int N[], int *alpha[])
Definition BernsteinBezier.cpp:138
MoFEM::BernsteinBezier::generateIndicesVertexTet
static MoFEMErrorCode generateIndicesVertexTet(const int N, int *alpha)
Definition BernsteinBezier.cpp:128
MoFEM::BernsteinBezier::generateIndicesTetTet
static MoFEMErrorCode generateIndicesTetTet(const int N, int *alpha)
Definition BernsteinBezier.cpp:201
MoFEM::EntPolynomialBaseCtx
Class used to pass element data to calculate base functions on tet,triangle,edge.
Definition EntPolynomialBaseCtx.hpp:22
MoFEM::EntPolynomialBaseCtx::copyNodeBase
const FieldApproximationBase copyNodeBase
Definition EntPolynomialBaseCtx.hpp:41
MoFEM::EntPolynomialBaseCtx::basePolynomialsType0
PetscErrorCode(* basePolynomialsType0)(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
Definition EntPolynomialBaseCtx.hpp:27
MoFEM::EntPolynomialBaseCtx::spaceContinuity
const FieldContinuity spaceContinuity
Definition EntPolynomialBaseCtx.hpp:38
MoFEM::EntPolynomialBaseCtx::bAse
const FieldApproximationBase bAse
Definition EntPolynomialBaseCtx.hpp:39
MoFEM::EntitiesFieldData
data structure for finite element entity
Definition EntitiesFieldData.hpp:38
MoFEM::EntitiesFieldData::facesNodes
MatrixInt facesNodes
nodes on finite element faces
Definition EntitiesFieldData.hpp:43
MoFEM::EntitiesFieldData::spacesOnEntities
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
Definition EntitiesFieldData.hpp:47
MoFEM::EntitiesFieldData::dataOnEntities
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
Definition EntitiesFieldData.hpp:55
MoFEM::Field::maxBrokenDofsOrder
static constexpr int maxBrokenDofsOrder
Maximum order for broken space DOFs.
Definition FieldMultiIndices.hpp:241
MoFEM::TetPolynomialBase
Calculate base functions on tetrahedral.
Definition TetPolynomialBase.hpp:17
MoFEM::TetPolynomialBase::diffN_face_edge
ublas::matrix< MatrixDouble > diffN_face_edge
Definition TetPolynomialBase.hpp:170
MoFEM::TetPolynomialBase::cTx
EntPolynomialBaseCtx * cTx
Definition TetPolynomialBase.hpp:71
MoFEM::TetPolynomialBase::getValueHdivAinsworthBaseImpl
MoFEMErrorCode getValueHdivAinsworthBaseImpl(MatrixDouble &pts, MatrixDouble &shape_functions, MatrixDouble &diff_shape_functions, int volume_order, std::array< int, 4 > &faces_order, std::array< int, 3 *4 > &faces_nodes, boost::function< int(int)> broken_nbfacetri_edge_hdiv, boost::function< int(int)> broken_nbfacetri_face_hdiv, boost::function< int(int)> broken_nbvolumetet_edge_hdiv, boost::function< int(int)> broken_nbvolumetet_face_hdiv, boost::function< int(int)> broken_nbvolumetet_volume_hdiv)
Definition TetPolynomialBase.cpp:781
MoFEM::TetPolynomialBase::getValueH1AinsworthBase
MoFEMErrorCode getValueH1AinsworthBase(MatrixDouble &pts)
Definition TetPolynomialBase.cpp:131
MoFEM::TetPolynomialBase::getValueHdivAinsworthBrokenBase
MoFEMErrorCode getValueHdivAinsworthBrokenBase(MatrixDouble &pts)
Definition TetPolynomialBase.cpp:1180
MoFEM::TetPolynomialBase::diffN_volume_face
ublas::vector< MatrixDouble > diffN_volume_face
Definition TetPolynomialBase.hpp:173
MoFEM::TetPolynomialBase::getValueL2BernsteinBezierBase
MoFEMErrorCode getValueL2BernsteinBezierBase(MatrixDouble &pts)
Definition TetPolynomialBase.cpp:605
MoFEM::TetPolynomialBase::getValueHcurlDemkowiczBase
MoFEMErrorCode getValueHcurlDemkowiczBase(MatrixDouble &pts)
Definition TetPolynomialBase.cpp:1913
MoFEM::TetPolynomialBase::getValueHcurl
MoFEMErrorCode getValueHcurl(MatrixDouble &pts)
Get base functions for Hcurl space.
Definition TetPolynomialBase.cpp:2045
MoFEM::TetPolynomialBase::getValueH1BernsteinBezierBase
MoFEMErrorCode getValueH1BernsteinBezierBase(MatrixDouble &pts)
Definition TetPolynomialBase.cpp:243
MoFEM::TetPolynomialBase::setDofsSideMap
static MoFEMErrorCode setDofsSideMap(const FieldSpace space, const FieldContinuity continuity, const FieldApproximationBase base, DofsSideMap &)
Set map of dof to side number.
Definition TetPolynomialBase.cpp:2145
MoFEM::TetPolynomialBase::getValueH1
MoFEMErrorCode getValueH1(MatrixDouble &pts)
Get base functions for H1 space.
Definition TetPolynomialBase.cpp:113
MoFEM::TetPolynomialBase::diffN_face_bubble
ublas::vector< MatrixDouble > diffN_face_bubble
Definition TetPolynomialBase.hpp:171
MoFEM::TetPolynomialBase::N_volume_edge
ublas::vector< MatrixDouble > N_volume_edge
Definition TetPolynomialBase.hpp:166
MoFEM::TetPolynomialBase::N_volume_face
ublas::vector< MatrixDouble > N_volume_face
Definition TetPolynomialBase.hpp:167
MoFEM::TetPolynomialBase::getValueHdiv
MoFEMErrorCode getValueHdiv(MatrixDouble &pts)
Get base functions for Hdiv space.
Definition TetPolynomialBase.cpp:1760
MoFEM::TetPolynomialBase::getValueHdivAinsworthBase
MoFEMErrorCode getValueHdivAinsworthBase(MatrixDouble &pts)
Definition TetPolynomialBase.cpp:912
MoFEM::TetPolynomialBase::getValueL2AinsworthBase
MoFEMErrorCode getValueL2AinsworthBase(MatrixDouble &pts)
Definition TetPolynomialBase.cpp:544
MoFEM::TetPolynomialBase::N_face_bubble
ublas::vector< MatrixDouble > N_face_bubble
Definition TetPolynomialBase.hpp:165
MoFEM::TetPolynomialBase::N_face_edge
ublas::matrix< MatrixDouble > N_face_edge
Definition TetPolynomialBase.hpp:164
MoFEM::TetPolynomialBase::getValueL2
MoFEMErrorCode getValueL2(MatrixDouble &pts)
Get base functions for L2 space.
Definition TetPolynomialBase.cpp:525
MoFEM::TetPolynomialBase::getValueHdivDemkowiczBase
MoFEMErrorCode getValueHdivDemkowiczBase(MatrixDouble &pts)
Definition TetPolynomialBase.cpp:1471
MoFEM::TetPolynomialBase::TetPolynomialBase
TetPolynomialBase(const void *ptr=nullptr)
Definition TetPolynomialBase.cpp:95
MoFEM::TetPolynomialBase::getValue
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Definition TetPolynomialBase.cpp:2064
MoFEM::TetPolynomialBase::getValueHdivDemkowiczBrokenBase
MoFEMErrorCode getValueHdivDemkowiczBrokenBase(MatrixDouble &pts)
Definition TetPolynomialBase.cpp:1608
MoFEM::TetPolynomialBase::diffN_volume_edge
ublas::vector< MatrixDouble > diffN_volume_edge
Definition TetPolynomialBase.hpp:172
MoFEM::TetPolynomialBase::getValueHcurlAinsworthBase
MoFEMErrorCode getValueHcurlAinsworthBase(MatrixDouble &pts)
Definition TetPolynomialBase.cpp:1795
MoFEM::TetPolynomialBase::setDofsSideMapHdiv
static MoFEMErrorCode setDofsSideMapHdiv(const FieldContinuity continuity, const FieldApproximationBase base, DofsSideMap &dofs_side_map)
Set the Dofs Side Map Hdiv object.
Definition TetPolynomialBase.cpp:2175
MoFEM::Tools::shapeFunMBTET
static MoFEMErrorCode shapeFunMBTET(double *shape, const double *ksi, const double *eta, const double *zeta, const double nb)
Calculate shape functions on tetrahedron.
Definition Tools.hpp:759
MoFEM::Tools::diffShapeFunMBTET
static constexpr std::array< double, 12 > diffShapeFunMBTET
Definition Tools.hpp:271
MoFEM::UnknownInterface
base class for all interface classes
Definition UnknownInterface.hpp:34
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition UnknownInterface.hpp:93
TetBaseCache::BaseCacheMI
boost::multi_index_container< BaseCacheItem, boost::multi_index::indexed_by< boost::multi_index::hashed_unique< composite_key< BaseCacheItem, member< BaseCacheItem, int, &BaseCacheItem::order >, member< BaseCacheItem, int, &BaseCacheItem::nb_gauss_pts > > > > > BaseCacheMI
Definition TetPolynomialBase.cpp:47
TetBaseCache::hdivBaseInterior
static std::array< std::map< const void *, BaseCacheMI >, LASTBASE > hdivBaseInterior
Definition TetPolynomialBase.cpp:68
TetBaseCache::HDivBaseFaceCacheMI
boost::multi_index_container< HDivBaseCacheItem, boost::multi_index::indexed_by< boost::multi_index::hashed_unique< composite_key< HDivBaseCacheItem, member< HDivBaseCacheItem, int, &HDivBaseCacheItem::order >, member< HDivBaseCacheItem, int, &HDivBaseCacheItem::nb_gauss_pts >, member< HDivBaseCacheItem, int, &HDivBaseCacheItem::n0 >, member< HDivBaseCacheItem, int, &HDivBaseCacheItem::n1 >, member< HDivBaseCacheItem, int, &HDivBaseCacheItem::n2 > > > > > HDivBaseFaceCacheMI
Definition TetPolynomialBase.cpp:65
TetBaseCache::hDivBaseFace
static std::array< std::map< const void *, HDivBaseFaceCacheMI >, LASTBASE > hDivBaseFace
Definition TetPolynomialBase.cpp:72
TetBaseCache::l2BaseInterior
static std::array< std::map< const void *, BaseCacheMI >, LASTBASE > l2BaseInterior
Definition TetPolynomialBase.cpp:74
TetBaseCache::hdivBrokenBaseInterior
static std::array< std::map< const void *, BaseCacheMI >, LASTBASE > hdivBrokenBaseInterior
Definition TetPolynomialBase.cpp:70
Generated by Doxygen 1.9.8 and hosted at University of Glasgow