v0.15.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
11
13 int order;
15 mutable MatrixDouble N;
17 };
18
20
21 int order;
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;
33 };
34
35 using BaseCacheMI = boost::multi_index_container<
37 boost::multi_index::indexed_by<
38
39 boost::multi_index::hashed_unique<
40
41 composite_key<
42
44 member<BaseCacheItem, int, &BaseCacheItem::order>,
45 member<BaseCacheItem, int, &BaseCacheItem::nb_gauss_pts>>>>
46
47 >;
48
49 using HDivBaseFaceCacheMI = boost::multi_index_container<
51 boost::multi_index::indexed_by<
52
53 boost::multi_index::hashed_unique<
54
55 composite_key<
56
58 member<HDivBaseCacheItem, int, &HDivBaseCacheItem::order>,
59 member<HDivBaseCacheItem, int,
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>
69 static std::array<std::map<const void *, BaseCacheMI>, LASTBASE>
71 static std::array<std::map<const void *, HDivBaseFaceCacheMI>, LASTBASE>
73 static std::array<std::map<const void *, BaseCacheMI>, LASTBASE>
75};
76
77std::array<std::map<const void *, TetBaseCache::BaseCacheMI>, LASTBASE>
79std::array<std::map<const void *, TetBaseCache::BaseCacheMI>, LASTBASE>
81std::array<std::map<const void *, TetBaseCache::HDivBaseFaceCacheMI>, LASTBASE>
83std::array<std::map<const void *, TetBaseCache::BaseCacheMI>, LASTBASE>
85
87TetPolynomialBase::query_interface(boost::typeindex::type_index type_index,
88 UnknownInterface **iface) const {
89
91 *iface = const_cast<TetPolynomialBase *>(this);
93}
94
95TetPolynomialBase::TetPolynomialBase(const void *ptr) : vPtr(ptr) {}
96
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) {
109 }
110 }
111}
112
115
116 switch (cTx->bAse) {
120 break;
123 break;
124 default:
125 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
126 }
127
129}
130
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) =
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 }
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 }
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);
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
240}
241
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",
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);
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);
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 }
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
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);
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
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
523}
524
527
528 switch (cTx->bAse) {
533 break;
536 break;
537 default:
538 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
539 }
540
542}
543
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) =
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)
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;
583 }
584 }
585
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
602}
603
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",
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) +
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
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)};
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 };
746 face_n.data(), tet_face_ptr.data());
747 if (order > 3)
749 order,
750 &tet_alpha(
751 4 + 6 * NBEDGE_H1(order) + 4 * NBFACETRI_H1(order), 0));
752 }
753 }
754
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
779}
780
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) {
798
799 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
800 double *diffL, const int dim) =
802
803 int nb_gauss_pts = pts.size2();
804
805 // face shape functions
806
807 double *phi_f_e[4][3];
808 double *phi_f[4];
809 double *diff_phi_f_e[4][3];
810 double *diff_phi_f[4];
811
812 N_face_edge.resize(4, 3, false);
813 N_face_bubble.resize(4, false);
814 diffN_face_edge.resize(4, 3, false);
815 diffN_face_bubble.resize(4, false);
816
817 for (int ff = 0; ff != 4; ++ff) {
818 const auto face_edge_dofs = NBFACETRI_AINSWORTH_EDGE_HDIV(
819 broken_nbfacetri_edge_hdiv(faces_order[ff]));
820 // three edges on face
821 for (int ee = 0; ee < 3; ee++) {
822 N_face_edge(ff, ee).resize(nb_gauss_pts, 3 * face_edge_dofs, false);
823 diffN_face_edge(ff, ee).resize(nb_gauss_pts, 9 * face_edge_dofs, false);
824 phi_f_e[ff][ee] = &*N_face_edge(ff, ee).data().begin();
825 diff_phi_f_e[ff][ee] = &*diffN_face_edge(ff, ee).data().begin();
826 }
827 auto face_bubble_dofs = NBFACETRI_AINSWORTH_FACE_HDIV(
828 broken_nbfacetri_face_hdiv(faces_order[ff]));
829 N_face_bubble[ff].resize(nb_gauss_pts, 3 * face_bubble_dofs, false);
830 diffN_face_bubble[ff].resize(nb_gauss_pts, 9 * face_bubble_dofs, false);
831 phi_f[ff] = &*(N_face_bubble[ff].data().begin());
832 diff_phi_f[ff] = &*(diffN_face_bubble[ff].data().begin());
833 }
834
835 constexpr int nb_nodes_on_tet = 4;
836
837 for (int ff = 0; ff < 4; ff++) {
839 &faces_nodes[3 * ff], broken_nbfacetri_edge_hdiv(faces_order[ff]),
840 &*shape_functions.data().begin(), &*diff_shape_functions.data().begin(),
841 phi_f_e[ff], diff_phi_f_e[ff], nb_gauss_pts, nb_nodes_on_tet,
842 base_polynomials);
843 }
844
845 for (int ff = 0; ff < 4; ff++) {
847 &faces_nodes[3 * ff], broken_nbfacetri_face_hdiv(faces_order[ff]),
848 &*shape_functions.data().begin(), &*diff_shape_functions.data().begin(),
849 phi_f[ff], diff_phi_f[ff], nb_gauss_pts, nb_nodes_on_tet,
850 base_polynomials);
851 }
852
853 // volume shape functions
854
855 double *phi_v_e[6];
856 double *phi_v_f[4];
857 double *phi_v;
858 double *diff_phi_v_e[6];
859 double *diff_phi_v_f[4];
860 double *diff_phi_v;
861
862 const auto volume_edge_dofs = NBVOLUMETET_AINSWORTH_EDGE_HDIV(
863 broken_nbvolumetet_edge_hdiv(volume_order));
864 N_volume_edge.resize(6, false);
865 diffN_volume_edge.resize(6, false);
866 for (int ee = 0; ee != 6; ++ee) {
867 N_volume_edge[ee].resize(nb_gauss_pts, 3 * volume_edge_dofs, false);
868 diffN_volume_edge[ee].resize(nb_gauss_pts, 9 * volume_edge_dofs, false);
869 phi_v_e[ee] = &*(N_volume_edge[ee].data().begin());
870 diff_phi_v_e[ee] = &*(diffN_volume_edge[ee].data().begin());
871 }
872 if (volume_edge_dofs)
874 broken_nbvolumetet_edge_hdiv(volume_order),
875 &*shape_functions.data().begin(), &*diff_shape_functions.data().begin(),
876 phi_v_e, diff_phi_v_e, nb_gauss_pts, base_polynomials);
877
878 const auto volume_face_dofs = NBVOLUMETET_AINSWORTH_FACE_HDIV(
879 broken_nbvolumetet_face_hdiv(volume_order));
880 N_volume_face.resize(4, false);
881 diffN_volume_face.resize(4, false);
882 for (int ff = 0; ff != 4; ++ff) {
883 N_volume_face[ff].resize(nb_gauss_pts, 3 * volume_face_dofs, false);
884 diffN_volume_face[ff].resize(nb_gauss_pts, 9 * volume_face_dofs, false);
885 phi_v_f[ff] = &*(N_volume_face[ff].data().begin());
886 diff_phi_v_f[ff] = &*(diffN_volume_face[ff].data().begin());
887 }
888 if (volume_face_dofs)
890 broken_nbvolumetet_face_hdiv(volume_order),
891 &*shape_functions.data().begin(), &*diff_shape_functions.data().begin(),
892 phi_v_f, diff_phi_v_f, nb_gauss_pts, base_polynomials);
893
894 auto volume_bubble_dofs = NBVOLUMETET_AINSWORTH_VOLUME_HDIV(
895 broken_nbvolumetet_volume_hdiv(volume_order));
896 N_volume_bubble.resize(nb_gauss_pts, 3 * volume_bubble_dofs, false);
897 diffN_volume_bubble.resize(nb_gauss_pts, 9 * volume_bubble_dofs, false);
898 phi_v = &*(N_volume_bubble.data().begin());
899 diff_phi_v = &*(diffN_volume_bubble.data().begin());
900 if (volume_bubble_dofs)
902 broken_nbvolumetet_volume_hdiv(volume_order),
903 &*shape_functions.data().begin(), &*diff_shape_functions.data().begin(),
904 phi_v, diff_phi_v, nb_gauss_pts, base_polynomials);
905
907}
908
911
912 std::array<int, 4> faces_order;
913 std::array<int, 4 * 3> faces_nodes;
914
916 EntitiesFieldData &data = cTx->dAta;
917
918 int volume_order = data.dataOnEntities[MBTET][0].getOrder();
919 std::copy(data.facesNodes.data().begin(), data.facesNodes.data().end(),
920 faces_nodes.begin());
921 for (int ff = 0; ff != 4; ff++) {
922 faces_order[ff] = cTx->dAta.dataOnEntities[MBTRI][ff].getOrder();
923 }
924
926 pts, data.dataOnEntities[MBVERTEX][0].getN(base),
927 data.dataOnEntities[MBVERTEX][0].getDiffN(base), volume_order,
928 faces_order, faces_nodes,
929
935
936 );
937
938 // Set shape functions into data structure Shape functions hast to be put
939 // in arrays in order which guarantee hierarchical series of degrees of
940 // freedom, i.e. in other words dofs form sub-entities has to be group
941 // by order.
942
943 FTENSOR_INDEX(3, i);
944 FTENSOR_INDEX(3, j);
945
946 int nb_gauss_pts = pts.size2();
947
948 // faces
949 if (data.dataOnEntities[MBTRI].size() != 4) {
950 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
951 }
952
953 // face-face
955 getFTensor1FromPtr<3>(&*(N_face_bubble[0].data().begin())),
956 getFTensor1FromPtr<3>(&*(N_face_bubble[1].data().begin())),
957 getFTensor1FromPtr<3>(&*(N_face_bubble[2].data().begin())),
958 getFTensor1FromPtr<3>(&*(N_face_bubble[3].data().begin()))};
959 FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_f_f[] = {
960 getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_bubble[0].data().begin())),
961 getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_bubble[1].data().begin())),
962 getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_bubble[2].data().begin())),
963 getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_bubble[3].data().begin()))};
964 // face-edge
966 getFTensor1FromPtr<3>(&*(N_face_edge(0, 0).data().begin())),
967 getFTensor1FromPtr<3>(&*(N_face_edge(0, 1).data().begin())),
968 getFTensor1FromPtr<3>(&*(N_face_edge(0, 2).data().begin())),
969 getFTensor1FromPtr<3>(&*(N_face_edge(1, 0).data().begin())),
970 getFTensor1FromPtr<3>(&*(N_face_edge(1, 1).data().begin())),
971 getFTensor1FromPtr<3>(&*(N_face_edge(1, 2).data().begin())),
972 getFTensor1FromPtr<3>(&*(N_face_edge(2, 0).data().begin())),
973 getFTensor1FromPtr<3>(&*(N_face_edge(2, 1).data().begin())),
974 getFTensor1FromPtr<3>(&*(N_face_edge(2, 2).data().begin())),
975 getFTensor1FromPtr<3>(&*(N_face_edge(3, 0).data().begin())),
976 getFTensor1FromPtr<3>(&*(N_face_edge(3, 1).data().begin())),
977 getFTensor1FromPtr<3>(&*(N_face_edge(3, 2).data().begin()))};
978 FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_f_e[] = {
979 getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(0, 0).data().begin())),
980 getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(0, 1).data().begin())),
981 getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(0, 2).data().begin())),
982 getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(1, 0).data().begin())),
983 getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(1, 1).data().begin())),
984 getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(1, 2).data().begin())),
985 getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(2, 0).data().begin())),
986 getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(2, 1).data().begin())),
987 getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(2, 2).data().begin())),
988 getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(3, 0).data().begin())),
989 getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(3, 1).data().begin())),
990 getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(3, 2).data().begin()))};
991
992 for (int ff = 0; ff != 4; ff++) {
993 int face_order = faces_order[ff];
994 auto face_dofs =
999 data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts,
1000 3 * face_dofs, false);
1001 data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
1002 9 * face_dofs, false);
1003 if (face_dofs) {
1004 double *base_ptr =
1005 &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
1006 double *diff_base_ptr =
1007 &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
1008 auto t_base = getFTensor1FromPtr<3>(base_ptr);
1009 auto t_diff_base = getFTensor2HVecFromPtr<3, 3>(diff_base_ptr);
1010
1011 auto max_face_order =
1012 std::max(face_order,
1014 max_face_order =
1015 std::max(max_face_order,
1017
1018 for (int gg = 0; gg != nb_gauss_pts; gg++) {
1019 for (int oo = 0; oo != max_face_order; oo++) {
1020
1021 // face-edge
1023 for (int dd = NBFACETRI_AINSWORTH_EDGE_HDIV(oo);
1024 dd != NBFACETRI_AINSWORTH_EDGE_HDIV(oo + 1); dd++) {
1025 for (int ee = 0; ee != 3; ++ee) {
1026 t_base(i) = t_base_f_e[ff * 3 + ee](i);
1027 ++t_base;
1028 ++t_base_f_e[ff * 3 + ee];
1029 }
1030 for (int ee = 0; ee != 3; ++ee) {
1031 t_diff_base(i, j) = t_diff_base_f_e[ff * 3 + ee](i, j);
1032 ++t_diff_base;
1033 ++t_diff_base_f_e[ff * 3 + ee];
1034 }
1035 }
1036
1037 // face-face
1039 for (int dd = NBFACETRI_AINSWORTH_FACE_HDIV(oo);
1040 dd != NBFACETRI_AINSWORTH_FACE_HDIV(oo + 1); dd++) {
1041 t_base(i) = t_base_f_f[ff](i);
1042 ++t_base;
1043 ++t_base_f_f[ff];
1044 t_diff_base(i, j) = t_diff_base_f_f[ff](i, j);
1045 ++t_diff_base;
1046 ++t_diff_base_f_f[ff];
1047 }
1048 }
1049 }
1050 }
1051 }
1052
1053 // volume
1054 int volume_dofs =
1061 data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 3 * volume_dofs,
1062 false);
1063 data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
1064 9 * volume_dofs, false);
1065 if (volume_dofs) {
1066 double *base_ptr =
1067 &*data.dataOnEntities[MBTET][0].getN(base).data().begin();
1068 double *diff_base_ptr =
1069 &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin();
1070 auto t_base = getFTensor1FromPtr<3>(base_ptr);
1071 auto t_diff_base = getFTensor2HVecFromPtr<3, 3>(diff_base_ptr);
1072
1073 // volume-edge
1075 getFTensor1FromPtr<3>(&*N_volume_edge[0].data().begin()),
1076 getFTensor1FromPtr<3>(&*N_volume_edge[1].data().begin()),
1077 getFTensor1FromPtr<3>(&*N_volume_edge[2].data().begin()),
1078 getFTensor1FromPtr<3>(&*N_volume_edge[3].data().begin()),
1079 getFTensor1FromPtr<3>(&*N_volume_edge[4].data().begin()),
1080 getFTensor1FromPtr<3>(&*N_volume_edge[5].data().begin())};
1081 FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_v_e[] = {
1087 getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_edge[5].data().begin())};
1088
1089 // volume-faces
1091 getFTensor1FromPtr<3>(&*N_volume_face[0].data().begin()),
1092 getFTensor1FromPtr<3>(&*N_volume_face[1].data().begin()),
1093 getFTensor1FromPtr<3>(&*N_volume_face[2].data().begin()),
1094 getFTensor1FromPtr<3>(&*N_volume_face[3].data().begin())};
1095 FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_v_f[] = {
1099 getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_face[3].data().begin())};
1100
1101 // volume-bubble
1102 base_ptr = &*(N_volume_bubble.data().begin());
1103 diff_base_ptr = &*(diffN_volume_bubble.data().begin());
1104 auto t_base_v = getFTensor1FromPtr<3>(base_ptr);
1105 auto t_diff_base_v = getFTensor2HVecFromPtr<3, 3>(diff_base_ptr);
1106
1107 auto max_volume_order = std::max(
1108 volume_order,
1110 max_volume_order = std::max(
1111 max_volume_order,
1113 max_volume_order = std::max(
1114 max_volume_order,
1116
1117 for (int gg = 0; gg != nb_gauss_pts; gg++) {
1118 for (int oo = 0; oo < max_volume_order; oo++) {
1119
1120 // volume-edge
1121 if (oo <
1123 for (int dd = NBVOLUMETET_AINSWORTH_EDGE_HDIV(oo);
1124 dd != NBVOLUMETET_AINSWORTH_EDGE_HDIV(oo + 1); dd++) {
1125 for (int ee = 0; ee < 6; ee++) {
1126 t_base(i) = t_base_v_e[ee](i);
1127 ++t_base;
1128 ++t_base_v_e[ee];
1129 t_diff_base(i, j) = t_diff_base_v_e[ee](i, j);
1130 ++t_diff_base;
1131 ++t_diff_base_v_e[ee];
1132 }
1133 }
1134
1135 // volume-face
1136 if (oo <
1138 for (int dd = NBVOLUMETET_AINSWORTH_FACE_HDIV(oo);
1139 dd < NBVOLUMETET_AINSWORTH_FACE_HDIV(oo + 1); dd++) {
1140 for (int ff = 0; ff < 4; ff++) {
1141 t_base(i) = t_base_v_f[ff](i);
1142 ++t_base;
1143 ++t_base_v_f[ff];
1144 t_diff_base(i, j) = t_diff_base_v_f[ff](i, j);
1145 ++t_diff_base;
1146 ++t_diff_base_v_f[ff];
1147 }
1148 }
1149
1150 // volume-bubble
1151 if (oo <
1153 for (int dd = NBVOLUMETET_AINSWORTH_VOLUME_HDIV(oo);
1154 dd < NBVOLUMETET_AINSWORTH_VOLUME_HDIV(oo + 1); dd++) {
1155 t_base(i) = t_base_v(i);
1156 ++t_base;
1157 ++t_base_v;
1158 t_diff_base(i, j) = t_diff_base_v(i, j);
1159 ++t_diff_base;
1160 ++t_diff_base_v;
1161 }
1162 }
1163 }
1164 }
1165
1167}
1168
1172
1173 // Set shape functions into data structure Shape functions has to be put
1174 // in arrays in order which guarantee hierarchical series of degrees of
1175 // freedom, i.e. in other words dofs form sub-entities has to be group
1176 // by order.
1177
1179 EntitiesFieldData &data = cTx->dAta;
1180
1181 FTENSOR_INDEX(3, i);
1182 FTENSOR_INDEX(3, j);
1183
1184 int volume_order = data.dataOnEntities[MBTET][0].getOrder();
1185 int nb_gauss_pts = pts.size2();
1186 int nb_dofs_face =
1191 int nb_dofs_volume =
1198
1199 int nb_dofs = 4 * nb_dofs_face + nb_dofs_volume;
1200 data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 3 * nb_dofs,
1201 false);
1202 data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts, 9 * nb_dofs,
1203 false);
1204 if (nb_dofs == 0)
1206
1207 auto get_interior_cache = [this](auto base) -> TetBaseCache::BaseCacheMI * {
1208 if (vPtr) {
1209 auto it = TetBaseCache::hdivBrokenBaseInterior[base].find(vPtr);
1210 if (it != TetBaseCache::hdivBrokenBaseInterior[base].end()) {
1211 return &it->second;
1212 }
1213 }
1214 return nullptr;
1215 };
1216
1217 auto interior_cache_ptr = get_interior_cache(base);
1218
1219 if (interior_cache_ptr) {
1220 auto it =
1221 interior_cache_ptr->find(boost::make_tuple(volume_order, nb_gauss_pts));
1222 if (it != interior_cache_ptr->end()) {
1223 noalias(data.dataOnEntities[MBTET][0].getN(base)) = it->N;
1224 noalias(data.dataOnEntities[MBTET][0].getDiffN(base)) = it->diffN;
1226 }
1227 }
1228
1229 std::array<int, 4 * 3> faces_nodes = {0, 1, 3, 1, 2, 3, 0, 3, 2, 0, 2, 1};
1230 std::array<int, 4> faces_order{volume_order, volume_order, volume_order,
1231 volume_order};
1233 pts, data.dataOnEntities[MBVERTEX][0].getN(base),
1234 data.dataOnEntities[MBVERTEX][0].getDiffN(base), volume_order,
1235 faces_order, faces_nodes,
1236
1242
1243 );
1244
1245 auto *base_ptr = &*data.dataOnEntities[MBTET][0].getN(base).data().begin();
1246 auto *diff_base_ptr =
1247 &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin();
1248 auto t_base = getFTensor1FromPtr<3>(base_ptr);
1249 auto t_diff_base = getFTensor2HVecFromPtr<3, 3>(diff_base_ptr);
1250
1251 // face-edge
1253 getFTensor1FromPtr<3>(&*(N_face_edge(0, 0).data().begin())),
1254 getFTensor1FromPtr<3>(&*(N_face_edge(0, 1).data().begin())),
1255 getFTensor1FromPtr<3>(&*(N_face_edge(0, 2).data().begin())),
1256 getFTensor1FromPtr<3>(&*(N_face_edge(1, 0).data().begin())),
1257 getFTensor1FromPtr<3>(&*(N_face_edge(1, 1).data().begin())),
1258 getFTensor1FromPtr<3>(&*(N_face_edge(1, 2).data().begin())),
1259 getFTensor1FromPtr<3>(&*(N_face_edge(2, 0).data().begin())),
1260 getFTensor1FromPtr<3>(&*(N_face_edge(2, 1).data().begin())),
1261 getFTensor1FromPtr<3>(&*(N_face_edge(2, 2).data().begin())),
1262 getFTensor1FromPtr<3>(&*(N_face_edge(3, 0).data().begin())),
1263 getFTensor1FromPtr<3>(&*(N_face_edge(3, 1).data().begin())),
1264 getFTensor1FromPtr<3>(&*(N_face_edge(3, 2).data().begin()))};
1265 FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_f_e[] = {
1266 getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(0, 0).data().begin())),
1267 getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(0, 1).data().begin())),
1268 getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(0, 2).data().begin())),
1269 getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(1, 0).data().begin())),
1270 getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(1, 1).data().begin())),
1271 getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(1, 2).data().begin())),
1272 getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(2, 0).data().begin())),
1273 getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(2, 1).data().begin())),
1274 getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(2, 2).data().begin())),
1275 getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(3, 0).data().begin())),
1276 getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(3, 1).data().begin())),
1277 getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(3, 2).data().begin()))};
1278
1279 // face-face
1281 getFTensor1FromPtr<3>(&*(N_face_bubble[0].data().begin())),
1282 getFTensor1FromPtr<3>(&*(N_face_bubble[1].data().begin())),
1283 getFTensor1FromPtr<3>(&*(N_face_bubble[2].data().begin())),
1284 getFTensor1FromPtr<3>(&*(N_face_bubble[3].data().begin()))};
1285 FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_f_f[] = {
1286 getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_bubble[0].data().begin())),
1287 getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_bubble[1].data().begin())),
1288 getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_bubble[2].data().begin())),
1289 getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_bubble[3].data().begin()))};
1290
1291 // volume-edge
1293 getFTensor1FromPtr<3>(&*N_volume_edge[0].data().begin()),
1294 getFTensor1FromPtr<3>(&*N_volume_edge[1].data().begin()),
1295 getFTensor1FromPtr<3>(&*N_volume_edge[2].data().begin()),
1296 getFTensor1FromPtr<3>(&*N_volume_edge[3].data().begin()),
1297 getFTensor1FromPtr<3>(&*N_volume_edge[4].data().begin()),
1298 getFTensor1FromPtr<3>(&*N_volume_edge[5].data().begin())};
1299 FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_v_e[] = {
1305 getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_edge[5].data().begin())};
1306
1307 // volume-faces
1309 getFTensor1FromPtr<3>(&*N_volume_face[0].data().begin()),
1310 getFTensor1FromPtr<3>(&*N_volume_face[1].data().begin()),
1311 getFTensor1FromPtr<3>(&*N_volume_face[2].data().begin()),
1312 getFTensor1FromPtr<3>(&*N_volume_face[3].data().begin())};
1313 FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_v_f[] = {
1317 getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_face[3].data().begin())};
1318
1319 // volume-bubble
1320 auto *base_vol_ptr = &*(N_volume_bubble.data().begin());
1321 auto *diff_base_vol_ptr = &*(diffN_volume_bubble.data().begin());
1322 auto t_base_v = getFTensor1FromPtr<3>(base_vol_ptr);
1323 auto t_diff_base_v = getFTensor2HVecFromPtr<3, 3>(diff_base_vol_ptr);
1324
1325 int count_dofs = 0;
1326 int count_dofs_face = 0;
1327 int count_dofs_volume = 0;
1328
1329 auto max_volume_order =
1330 std::max(volume_order,
1332 max_volume_order =
1333 std::max(max_volume_order,
1335 max_volume_order =
1336 std::max(max_volume_order,
1338 max_volume_order =
1339 std::max(max_volume_order,
1341 max_volume_order = std::max(
1342 max_volume_order,
1344
1345 for (int gg = 0; gg != nb_gauss_pts; gg++) {
1346 for (int oo = 0; oo < max_volume_order; oo++) {
1347
1348 // faces-edge (((P) > 0) ? (P) : 0)
1350 for (int dd = NBFACETRI_AINSWORTH_EDGE_HDIV(oo);
1351 dd != NBFACETRI_AINSWORTH_EDGE_HDIV(oo + 1); dd++) {
1352 for (auto ff = 0; ff != 4; ++ff) {
1353 for (int ee = 0; ee != 3; ++ee) {
1354 t_base(i) = t_base_f_e[ff * 3 + ee](i);
1355 ++t_base;
1356 ++t_base_f_e[ff * 3 + ee];
1357 ++count_dofs;
1358 ++count_dofs_face;
1359 }
1360 for (int ee = 0; ee != 3; ++ee) {
1361 t_diff_base(i, j) = t_diff_base_f_e[ff * 3 + ee](i, j);
1362 ++t_diff_base;
1363 ++t_diff_base_f_e[ff * 3 + ee];
1364 }
1365 }
1366 }
1367
1368 // face-face (P - 1) * (P - 2) / 2
1370 for (int dd = NBFACETRI_AINSWORTH_FACE_HDIV(oo);
1371 dd != NBFACETRI_AINSWORTH_FACE_HDIV(oo + 1); dd++) {
1372 for (auto ff = 0; ff != 4; ++ff) {
1373 t_base(i) = t_base_f_f[ff](i);
1374 ++t_base;
1375 ++t_base_f_f[ff];
1376 t_diff_base(i, j) = t_diff_base_f_f[ff](i, j);
1377 ++t_diff_base;
1378 ++t_diff_base_f_f[ff];
1379 ++count_dofs;
1380 ++count_dofs_face;
1381 }
1382 }
1383
1384 // volume-edge (P - 1)
1386 for (int dd = NBVOLUMETET_AINSWORTH_EDGE_HDIV(oo);
1387 dd != NBVOLUMETET_AINSWORTH_EDGE_HDIV(oo + 1); dd++) {
1388 for (int ee = 0; ee < 6; ++ee) {
1389 t_base(i) = t_base_v_e[ee](i);
1390 ++t_base;
1391 ++t_base_v_e[ee];
1392 t_diff_base(i, j) = t_diff_base_v_e[ee](i, j);
1393 ++t_diff_base;
1394 ++t_diff_base_v_e[ee];
1395 ++count_dofs;
1396 ++count_dofs_volume;
1397 }
1398 }
1399
1400 // volume-face (P - 1) * (P - 2)
1402 for (int dd = NBVOLUMETET_AINSWORTH_FACE_HDIV(oo);
1403 dd != NBVOLUMETET_AINSWORTH_FACE_HDIV(oo + 1); dd++) {
1404 for (int ff = 0; ff < 4; ff++) {
1405 t_base(i) = t_base_v_f[ff](i);
1406 ++t_base;
1407 ++t_base_v_f[ff];
1408 t_diff_base(i, j) = t_diff_base_v_f[ff](i, j);
1409 ++t_diff_base;
1410 ++t_diff_base_v_f[ff];
1411 ++count_dofs;
1412 ++count_dofs_volume;
1413 }
1414 }
1415
1416 // volume-bubble (P - 3) * (P - 2) * (P - 1) / 2
1417 if (oo <
1419 for (int dd = NBVOLUMETET_AINSWORTH_VOLUME_HDIV(oo);
1420 dd != NBVOLUMETET_AINSWORTH_VOLUME_HDIV(oo + 1); dd++) {
1421 t_base(i) = t_base_v(i);
1422 ++t_base;
1423 ++t_base_v;
1424 t_diff_base(i, j) = t_diff_base_v(i, j);
1425 ++t_diff_base;
1426 ++t_diff_base_v;
1427 ++count_dofs;
1428 ++count_dofs_volume;
1429 }
1430 }
1431 }
1432
1433#ifndef NDEBUG
1434 if (nb_dofs != count_dofs / nb_gauss_pts) {
1435 MOFEM_LOG_CHANNEL("SELF");
1436 MOFEM_LOG("SELF", Sev::error) << "Nb dofs face: " << 4 * nb_dofs_face
1437 << " -> " << count_dofs_face / nb_gauss_pts;
1438 MOFEM_LOG("SELF", Sev::error) << "Nb dofs volume: " << nb_dofs_volume
1439 << " -> " << count_dofs_volume / nb_gauss_pts;
1440 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1441 "Number of dofs %d is different than expected %d",
1442 count_dofs / nb_gauss_pts, nb_dofs);
1443 }
1444#endif // NDEBUG
1445
1446 if (interior_cache_ptr) {
1447 auto p = interior_cache_ptr->emplace(
1448 TetBaseCache::BaseCacheItem{volume_order, nb_gauss_pts});
1449 p.first->N = data.dataOnEntities[MBTET][0].getN(base);
1450 p.first->diffN = data.dataOnEntities[MBTET][0].getDiffN(base);
1451 }
1452
1454}
1455
1458
1459 EntitiesFieldData &data = cTx->dAta;
1460 const FieldApproximationBase base = cTx->bAse;
1461 if (base != DEMKOWICZ_JACOBI_BASE) {
1462 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1463 "This should be used only with DEMKOWICZ_JACOBI_BASE "
1464 "but base is %s",
1466 }
1467 int nb_gauss_pts = pts.size2();
1468
1469 int volume_order = data.dataOnEntities[MBTET][0].getOrder();
1470
1471 int p_f[4];
1472 double *phi_f[4];
1473 double *diff_phi_f[4];
1474
1475 auto get_face_cache_ptr = [this]() -> TetBaseCache::HDivBaseFaceCacheMI * {
1476 if (vPtr) {
1479 return &it->second;
1480 }
1481 }
1482 return nullptr;
1483 };
1484
1485 auto face_cache_ptr = get_face_cache_ptr();
1486
1487 // Calculate base function on tet faces
1488 for (int ff = 0; ff != 4; ff++) {
1489 int face_order = data.dataOnEntities[MBTRI][ff].getOrder();
1490 int order = volume_order > face_order ? volume_order : face_order;
1491 data.dataOnEntities[MBTRI][ff].getN(base).resize(
1492 nb_gauss_pts, 3 * NBFACETRI_DEMKOWICZ_HDIV(order), false);
1493 data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(
1494 nb_gauss_pts, 9 * NBFACETRI_DEMKOWICZ_HDIV(order), false);
1496 continue;
1497
1498 if (face_cache_ptr) {
1499 auto it = face_cache_ptr->find(boost::make_tuple(
1500
1501 face_order, nb_gauss_pts,
1502
1503 data.facesNodes(ff, 0), data.facesNodes(ff, 1), data.facesNodes(ff, 2)
1504
1505 ));
1506 if (it != face_cache_ptr->end()) {
1507 noalias(data.dataOnEntities[MBTRI][ff].getN(base)) = it->N;
1508 noalias(data.dataOnEntities[MBTRI][ff].getDiffN(base)) = it->diffN;
1509 continue;
1510 }
1511 }
1512
1513 p_f[ff] = order;
1514 phi_f[ff] = &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
1515 diff_phi_f[ff] =
1516 &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
1517
1519 &data.facesNodes(ff, 0), order,
1520 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1521 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1522 phi_f[ff], diff_phi_f[ff], nb_gauss_pts, 4);
1523 if (face_cache_ptr) {
1524 auto p = face_cache_ptr->emplace(TetBaseCache::HDivBaseCacheItem{
1525 face_order, nb_gauss_pts, data.facesNodes(ff, 0),
1526 data.facesNodes(ff, 1), data.facesNodes(ff, 2)});
1527 p.first->N = data.dataOnEntities[MBTRI][ff].getN(base);
1528 p.first->diffN = data.dataOnEntities[MBTRI][ff].getDiffN(base);
1529 }
1530 }
1531
1532 auto get_interior_cache = [this]() -> TetBaseCache::BaseCacheMI * {
1533 if (vPtr) {
1534 auto it =
1537 return &it->second;
1538 }
1539 }
1540 return nullptr;
1541 };
1542
1543 auto interior_cache_ptr = get_interior_cache();
1544
1545 // Calculate base functions in tet interior
1546 if (NBVOLUMETET_DEMKOWICZ_HDIV(volume_order) > 0) {
1547 data.dataOnEntities[MBTET][0].getN(base).resize(
1548 nb_gauss_pts, 3 * NBVOLUMETET_DEMKOWICZ_HDIV(volume_order), false);
1549 data.dataOnEntities[MBTET][0].getDiffN(base).resize(
1550 nb_gauss_pts, 9 * NBVOLUMETET_DEMKOWICZ_HDIV(volume_order), false);
1551
1552 for (int v = 0; v != 1; ++v) {
1553 if (interior_cache_ptr) {
1554 auto it = interior_cache_ptr->find(
1555 boost::make_tuple(volume_order, nb_gauss_pts));
1556 if (it != interior_cache_ptr->end()) {
1557 noalias(data.dataOnEntities[MBTET][0].getN(base)) = it->N;
1558 noalias(data.dataOnEntities[MBTET][0].getDiffN(base)) = it->diffN;
1559 continue;
1560 }
1561 }
1562
1563 double *phi_v = &*data.dataOnEntities[MBTET][0].getN(base).data().begin();
1564 double *diff_phi_v =
1565 &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin();
1566
1568 volume_order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
1569 &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), p_f, phi_f,
1570 diff_phi_f, phi_v, diff_phi_v, nb_gauss_pts);
1571 if (interior_cache_ptr) {
1572 auto p = interior_cache_ptr->emplace(
1573 TetBaseCache::BaseCacheItem{volume_order, nb_gauss_pts});
1574 p.first->N = data.dataOnEntities[MBTET][0].getN(base);
1575 p.first->diffN = data.dataOnEntities[MBTET][0].getDiffN(base);
1576 }
1577 }
1578 }
1579
1580 // Set size of face base correctly
1581 for (int ff = 0; ff != 4; ff++) {
1582 int face_order = data.dataOnEntities[MBTRI][ff].getOrder();
1583 data.dataOnEntities[MBTRI][ff].getN(base).resize(
1584 nb_gauss_pts, 3 * NBFACETRI_DEMKOWICZ_HDIV(face_order), true);
1585 data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(
1586 nb_gauss_pts, 9 * NBFACETRI_DEMKOWICZ_HDIV(face_order), true);
1587 }
1588
1590}
1591
1595
1596 EntitiesFieldData &data = cTx->dAta;
1597 const FieldApproximationBase base = cTx->bAse;
1598 if (base != DEMKOWICZ_JACOBI_BASE) {
1599 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1600 "This should be used only with DEMKOWICZ_JACOBI_BASE "
1601 "but base is %s",
1603 }
1604 int nb_gauss_pts = pts.size2();
1605
1606 int volume_order = data.dataOnEntities[MBTET][0].getOrder();
1607 int nb_dofs_face = NBFACETRI_DEMKOWICZ_HDIV(volume_order);
1608 int nb_dofs_volume = NBVOLUMETET_DEMKOWICZ_HDIV(volume_order);
1609 int nb_dofs = 4 * nb_dofs_face + nb_dofs_volume;
1610 data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 3 * nb_dofs,
1611 false);
1612 data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts, 9 * nb_dofs,
1613 false);
1614 if (nb_dofs == 0)
1616
1617 auto get_interior_cache = [this]() -> TetBaseCache::BaseCacheMI * {
1618 if (vPtr) {
1619 auto it =
1621 vPtr);
1622 if (it !=
1624 return &it->second;
1625 }
1626 }
1627 return nullptr;
1628 };
1629
1630 auto interior_cache_ptr = get_interior_cache();
1631
1632 if (interior_cache_ptr) {
1633 auto it =
1634 interior_cache_ptr->find(boost::make_tuple(volume_order, nb_gauss_pts));
1635 if (it != interior_cache_ptr->end()) {
1636 noalias(data.dataOnEntities[MBTET][0].getN(base)) = it->N;
1637 noalias(data.dataOnEntities[MBTET][0].getDiffN(base)) = it->diffN;
1639 }
1640 }
1641
1642 std::array<MatrixDouble, 4> face_base_fun{
1643 MatrixDouble(nb_gauss_pts, 3 * nb_dofs_face),
1644 MatrixDouble(nb_gauss_pts, 3 * nb_dofs_face),
1645 MatrixDouble(nb_gauss_pts, 3 * nb_dofs_face),
1646 MatrixDouble(nb_gauss_pts, 3 * nb_dofs_face)};
1647 std::array<MatrixDouble, 4> face_diff_base{
1648 MatrixDouble(nb_gauss_pts, 9 * nb_dofs_face),
1649 MatrixDouble(nb_gauss_pts, 9 * nb_dofs_face),
1650 MatrixDouble(nb_gauss_pts, 9 * nb_dofs_face),
1651 MatrixDouble(nb_gauss_pts, 9 * nb_dofs_face)};
1652
1653 int faces_nodes[4][3] = {{0, 1, 3}, {1, 2, 3}, {0, 3, 2}, {0, 2, 1}};
1654
1655 std::array<int, 4> p_f{volume_order, volume_order, volume_order,
1656 volume_order};
1657 std::array<double *, 4> phi_f{
1658 &*face_base_fun[0].data().begin(), &*face_base_fun[1].data().begin(),
1659 &*face_base_fun[2].data().begin(), &*face_base_fun[3].data().begin()};
1660 std::array<double *, 4> diff_phi_f{
1661 &*face_diff_base[0].data().begin(), &*face_diff_base[1].data().begin(),
1662 &*face_diff_base[2].data().begin(), &*face_diff_base[3].data().begin()};
1663
1664 // Calculate base function on tet faces
1665 for (int ff = 0; ff != 4; ff++) {
1667 // &data.facesNodes(ff, 0)
1668 faces_nodes[ff], p_f[ff],
1669 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1670 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1671 phi_f[ff], diff_phi_f[ff], nb_gauss_pts, 4);
1672 }
1673
1674 MatrixDouble vol_bases(nb_gauss_pts, 3 * nb_dofs_volume);
1675 MatrixDouble vol_diff_bases(nb_gauss_pts, 9 * nb_dofs_volume);
1676 auto *phi_v = &*vol_bases.data().begin();
1677 auto *diff_phi_v = &*vol_diff_bases.data().begin();
1679 volume_order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
1680 &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), p_f.data(),
1681 phi_f.data(), diff_phi_f.data(), phi_v, diff_phi_v, nb_gauss_pts);
1682
1683 // faces
1685 getFTensor1FromPtr<3>(phi_f[0]), getFTensor1FromPtr<3>(phi_f[1]),
1686 getFTensor1FromPtr<3>(phi_f[2]), getFTensor1FromPtr<3>(phi_f[3])};
1687 FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_v_f[] = {
1688 getFTensor2HVecFromPtr<3, 3>(diff_phi_f[0]),
1689 getFTensor2HVecFromPtr<3, 3>(diff_phi_f[1]),
1690 getFTensor2HVecFromPtr<3, 3>(diff_phi_f[2]),
1691 getFTensor2HVecFromPtr<3, 3>(diff_phi_f[3])};
1692
1693 // volumes
1695 getFTensor1FromPtr<3>(&*vol_bases.data().begin());
1697 getFTensor2HVecFromPtr<3, 3>(&*vol_diff_bases.data().begin());
1698
1699 auto t_base = getFTensor1FromPtr<3>(
1700 &*data.dataOnEntities[MBTET][0].getN(base).data().begin());
1701 auto t_diff_base = getFTensor2HVecFromPtr<3, 3>(
1702 &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin());
1703
1704 FTENSOR_INDEX(3, i);
1705 FTENSOR_INDEX(3, j);
1706
1707 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
1708 for (int oo = 0; oo < volume_order; oo++) {
1709 // face
1710 for (auto dd = NBFACETRI_DEMKOWICZ_HDIV(oo);
1711 dd != NBFACETRI_DEMKOWICZ_HDIV(oo + 1); ++dd) {
1712 for (auto ff = 0; ff != 4; ++ff) {
1713 t_base(i) = t_base_v_f[ff](i);
1714 ++t_base;
1715 ++t_base_v_f[ff];
1716 t_diff_base(i, j) = t_diff_base_v_f[ff](i, j);
1717 ++t_diff_base;
1718 ++t_diff_base_v_f[ff];
1719 }
1720 }
1721 // volume
1722 for (auto dd = NBVOLUMETET_DEMKOWICZ_HDIV(oo);
1723 dd != NBVOLUMETET_DEMKOWICZ_HDIV(oo + 1); ++dd) {
1724 t_base(i) = t_base_v(i);
1725 ++t_base;
1726 ++t_base_v;
1727 t_diff_base(i, j) = t_diff_base_v(i, j);
1728 ++t_diff_base;
1729 ++t_diff_base_v;
1730 }
1731 }
1732 }
1733
1734 if (interior_cache_ptr) {
1735 auto p = interior_cache_ptr->emplace(
1736 TetBaseCache::BaseCacheItem{volume_order, nb_gauss_pts});
1737 p.first->N = data.dataOnEntities[MBTET][0].getN(base);
1738 p.first->diffN = data.dataOnEntities[MBTET][0].getDiffN(base);
1739 }
1740
1742}
1743
1746
1747 switch (cTx->spaceContinuity) {
1748 case CONTINUOUS:
1749 switch (cTx->bAse) {
1755 default:
1756 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
1757 }
1758 break;
1759 case DISCONTINUOUS:
1760 switch (cTx->bAse) {
1766 default:
1767 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
1768 }
1769 break;
1770
1771 default:
1772 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown continuity");
1773 }
1774
1776}
1777
1781
1782 EntitiesFieldData &data = cTx->dAta;
1783 const FieldApproximationBase base = cTx->bAse;
1784 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
1785 double *diffL, const int dim) =
1787
1788 int nb_gauss_pts = pts.size2();
1789
1790 // edges
1791 if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
1792 int sense[6], order[6];
1793 if (data.dataOnEntities[MBEDGE].size() != 6) {
1794 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1795 }
1796 double *hcurl_edge_n[6], *diff_hcurl_edge_n[6];
1797 for (int ee = 0; ee != 6; ee++) {
1798 if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
1799 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1800 "data inconsistency");
1801 }
1802 sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
1803 order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
1804 int nb_dofs =
1805 NBEDGE_AINSWORTH_HCURL(data.dataOnEntities[MBEDGE][ee].getOrder());
1806 data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts,
1807 3 * nb_dofs, false);
1808 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
1809 9 * nb_dofs, false);
1810 hcurl_edge_n[ee] =
1811 &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
1812 diff_hcurl_edge_n[ee] =
1813 &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
1814 }
1816 sense, order,
1817 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1818 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1819 hcurl_edge_n, diff_hcurl_edge_n, nb_gauss_pts, base_polynomials);
1820 } else {
1821 for (int ee = 0; ee != 6; ee++) {
1822 data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0, false);
1823 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
1824 false);
1825 }
1826 }
1827
1828 // triangles
1829 if (data.spacesOnEntities[MBTRI].test(HCURL)) {
1830 int order[4];
1831 // faces
1832 if (data.dataOnEntities[MBTRI].size() != 4) {
1833 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1834 }
1835 double *hcurl_base_n[4], *diff_hcurl_base_n[4];
1836 for (int ff = 0; ff != 4; ff++) {
1837 if (data.dataOnEntities[MBTRI][ff].getSense() == 0) {
1838 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1839 "data inconsistency");
1840 }
1841 order[ff] = data.dataOnEntities[MBTRI][ff].getOrder();
1842 int nb_dofs = NBFACETRI_AINSWORTH_HCURL(order[ff]);
1843 data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts,
1844 3 * nb_dofs, false);
1845 data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
1846 9 * nb_dofs, false);
1847 hcurl_base_n[ff] =
1848 &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
1849 diff_hcurl_base_n[ff] =
1850 &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
1851 }
1852 if (data.facesNodes.size1() != 4) {
1853 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1854 }
1855 if (data.facesNodes.size2() != 3) {
1856 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1857 }
1859 &*data.facesNodes.data().begin(), order,
1860 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1861 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1862 hcurl_base_n, diff_hcurl_base_n, nb_gauss_pts, base_polynomials);
1863 } else {
1864 for (int ff = 0; ff != 4; ff++) {
1865 data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, 0, false);
1866 data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts, 0,
1867 false);
1868 }
1869 }
1870
1871 if (data.spacesOnEntities[MBTET].test(HCURL)) {
1872
1873 // volume
1874 int order = data.dataOnEntities[MBTET][0].getOrder();
1875 int nb_vol_dofs = NBVOLUMETET_AINSWORTH_HCURL(order);
1876 data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts,
1877 3 * nb_vol_dofs, false);
1878 data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
1879 9 * nb_vol_dofs, false);
1881 data.dataOnEntities[MBTET][0].getOrder(),
1882 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1883 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1884 &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
1885 &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
1886 nb_gauss_pts, base_polynomials);
1887
1888 } else {
1889 data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 0, false);
1890 data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
1891 }
1892
1894}
1895
1899
1900 EntitiesFieldData &data = cTx->dAta;
1901 const FieldApproximationBase base = cTx->bAse;
1902 if (base != DEMKOWICZ_JACOBI_BASE) {
1903 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1904 "This should be used only with DEMKOWICZ_JACOBI_BASE "
1905 "but base is %s",
1907 }
1908
1909 int nb_gauss_pts = pts.size2();
1910
1911 // edges
1912 if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
1913 int sense[6], order[6];
1914 if (data.dataOnEntities[MBEDGE].size() != 6) {
1915 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1916 "wrong size of data structure, expected space for six edges "
1917 "but is %ld",
1918 data.dataOnEntities[MBEDGE].size());
1919 }
1920 double *hcurl_edge_n[6], *diff_hcurl_edge_n[6];
1921 for (int ee = 0; ee != 6; ee++) {
1922 if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
1923 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1924 "orintation of edges is not set");
1925 }
1926 sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
1927 order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
1928 int nb_dofs =
1929 NBEDGE_DEMKOWICZ_HCURL(data.dataOnEntities[MBEDGE][ee].getOrder());
1930 data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts,
1931 3 * nb_dofs, false);
1932 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
1933 9 * nb_dofs, false);
1934 hcurl_edge_n[ee] =
1935 &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
1936 diff_hcurl_edge_n[ee] =
1937 &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
1938 }
1940 sense, order,
1941 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1942 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1943 hcurl_edge_n, diff_hcurl_edge_n, nb_gauss_pts);
1944 } else {
1945 // No DOFs on edges, resize base function matrices, indicating that no
1946 // dofs on them.
1947 for (int ee = 0; ee != 6; ee++) {
1948 data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0, false);
1949 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
1950 false);
1951 }
1952 }
1953
1954 // triangles
1955 if (data.spacesOnEntities[MBTRI].test(HCURL)) {
1956 int order[4];
1957 // faces
1958 if (data.dataOnEntities[MBTRI].size() != 4) {
1959 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1960 "data structure for storing face h-curl base have wrong size "
1961 "should be four but is %ld",
1962 data.dataOnEntities[MBTRI].size());
1963 }
1964 double *hcurl_base_n[4], *diff_hcurl_base_n[4];
1965 for (int ff = 0; ff != 4; ff++) {
1966 if (data.dataOnEntities[MBTRI][ff].getSense() == 0) {
1967 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1968 "orintation of face is not set");
1969 }
1970 order[ff] = data.dataOnEntities[MBTRI][ff].getOrder();
1971 int nb_dofs = NBFACETRI_DEMKOWICZ_HCURL(order[ff]);
1972 data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts,
1973 3 * nb_dofs, false);
1974 data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
1975 9 * nb_dofs, false);
1976 hcurl_base_n[ff] =
1977 &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
1978 diff_hcurl_base_n[ff] =
1979 &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
1980 }
1981 if (data.facesNodes.size1() != 4) {
1982 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1983 "data inconsistency, should be four faces");
1984 }
1985 if (data.facesNodes.size2() != 3) {
1986 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1987 "data inconsistency, should be three nodes on face");
1988 }
1990 &*data.facesNodes.data().begin(), order,
1991 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1992 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1993 hcurl_base_n, diff_hcurl_base_n, nb_gauss_pts);
1994 } else {
1995 // No DOFs on faces, resize base function matrices, indicating that no
1996 // dofs on them.
1997 for (int ff = 0; ff != 4; ff++) {
1998 data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, 0, false);
1999 data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts, 0,
2000 false);
2001 }
2002 }
2003
2004 if (data.spacesOnEntities[MBTET].test(HCURL)) {
2005 // volume
2006 int order = data.dataOnEntities[MBTET][0].getOrder();
2007 int nb_vol_dofs = NBVOLUMETET_DEMKOWICZ_HCURL(order);
2008 data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts,
2009 3 * nb_vol_dofs, false);
2010 data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
2011 9 * nb_vol_dofs, false);
2013 data.dataOnEntities[MBTET][0].getOrder(),
2014 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
2015 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
2016 &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
2017 &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
2018 nb_gauss_pts);
2019 } else {
2020 // No DOFs on faces, resize base function matrices, indicating that no
2021 // dofs on them.
2022 data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 0, false);
2023 data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
2024 }
2025
2027}
2028
2031
2032 switch (cTx->bAse) {
2036 break;
2039 break;
2040 default:
2041 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
2042 }
2043
2045}
2046
2049 boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
2051
2053
2054 int nb_gauss_pts = pts.size2();
2055 if (!nb_gauss_pts)
2057
2058 if (pts.size1() < 3)
2059 SETERRQ(
2060 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2061 "Wrong dimension of pts, should be at least 3 rows with coordinates");
2062
2064 const FieldApproximationBase base = cTx->bAse;
2065 EntitiesFieldData &data = cTx->dAta;
2066 if (cTx->copyNodeBase == LASTBASE) {
2067 data.dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 4,
2068 false);
2070 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
2071 &pts(0, 0), &pts(1, 0), &pts(2, 0), nb_gauss_pts);
2072 } else {
2073 data.dataOnEntities[MBVERTEX][0].getNSharedPtr(base) =
2074 data.dataOnEntities[MBVERTEX][0].getNSharedPtr(cTx->copyNodeBase);
2075 }
2076 if (data.dataOnEntities[MBVERTEX][0].getN(base).size1() !=
2077 (unsigned int)nb_gauss_pts) {
2078 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2079 "Base functions or nodes has wrong number of integration points "
2080 "for base %s",
2082 }
2083 data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(4, 3, false);
2084 std::copy(Tools::diffShapeFunMBTET.begin(), Tools::diffShapeFunMBTET.end(),
2085 data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin());
2086 }
2087
2088 switch (cTx->spaceContinuity) {
2089 case CONTINUOUS:
2090
2091 switch (cTx->sPace) {
2092 case H1:
2093 CHKERR getValueH1(pts);
2094 break;
2095 case HDIV:
2096 CHKERR getValueHdiv(pts);
2097 break;
2098 case HCURL:
2099 CHKERR getValueHcurl(pts);
2100 break;
2101 case L2:
2102 CHKERR getValueL2(pts);
2103 break;
2104 default:
2105 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown space %s",
2107 }
2108 break;
2109
2110 case DISCONTINUOUS:
2111
2112 switch (cTx->sPace) {
2113 case HDIV:
2114 CHKERR getValueHdiv(pts);
2115 break;
2116 default:
2117 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown space %s",
2119 }
2120 break;
2121
2122 default:
2123 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown continuity");
2124 }
2125
2127}
2128
2130
2131 const FieldSpace space, const FieldContinuity continuity,
2132 const FieldApproximationBase base, DofsSideMap &dofs_side_map
2133
2134) {
2136
2137 switch (continuity) {
2138 case DISCONTINUOUS:
2139
2140 switch (space) {
2141 case HDIV:
2142 CHKERR setDofsSideMapHdiv(continuity, base, dofs_side_map);
2143 break;
2144 default:
2145 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown space %s",
2146 FieldSpaceNames[space]);
2147 }
2148 break;
2149
2150 default:
2151 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
2152 "Unknown (or not implemented) continuity");
2153 }
2154
2156}
2157
2160 const FieldApproximationBase base,
2161 DofsSideMap &dofs_side_map) {
2163
2164 // That has to be consistent with implementation of getValueHdiv for
2165 // particular base functions.
2166
2167 auto set_ainsworth = [&dofs_side_map]() {
2169
2170 dofs_side_map.clear();
2171
2172 int dof = 0;
2173 for (int oo = 0; oo < Field::maxBrokenDofsOrder; oo++) {
2174
2175 // faces-edge (((P) > 0) ? (P) : 0)
2176 for (int dd = NBFACETRI_AINSWORTH_EDGE_HDIV(oo);
2177 dd != NBFACETRI_AINSWORTH_EDGE_HDIV(oo + 1); dd++) {
2178 for (auto ff = 0; ff != 4; ++ff) {
2179 for (int ee = 0; ee != 3; ++ee) {
2180 dofs_side_map.insert(DofsSideMapData{MBTRI, ff, dof});
2181 ++dof;
2182 }
2183 }
2184 }
2185
2186 // face-face (P - 1) * (P - 2) / 2
2187 for (int dd = NBFACETRI_AINSWORTH_FACE_HDIV(oo);
2188 dd != NBFACETRI_AINSWORTH_FACE_HDIV(oo + 1); dd++) {
2189 for (auto ff = 0; ff != 4; ++ff) {
2190 dofs_side_map.insert(DofsSideMapData{MBTRI, ff, dof});
2191 ++dof;
2192 }
2193 }
2194
2195 // volume-edge (P - 1)
2196 for (int dd = NBVOLUMETET_AINSWORTH_EDGE_HDIV(oo);
2197 dd != NBVOLUMETET_AINSWORTH_EDGE_HDIV(oo + 1); dd++) {
2198 for (int ee = 0; ee < 6; ee++) {
2199 dofs_side_map.insert(DofsSideMapData{MBTET, 0, dof});
2200 ++dof;
2201 }
2202 }
2203 // volume-face (P - 1) * (P - 2)
2204 for (int dd = NBVOLUMETET_AINSWORTH_FACE_HDIV(oo);
2205 dd < NBVOLUMETET_AINSWORTH_FACE_HDIV(oo + 1); dd++) {
2206 for (int ff = 0; ff < 4; ff++) {
2207 dofs_side_map.insert(DofsSideMapData{MBTET, 0, dof});
2208 ++dof;
2209 }
2210 }
2211 // volume-bubble (P - 3) * (P - 2) * (P - 1) / 2
2212 for (int dd = NBVOLUMETET_AINSWORTH_VOLUME_HDIV(oo);
2213 dd < NBVOLUMETET_AINSWORTH_VOLUME_HDIV(oo + 1); dd++) {
2214 dofs_side_map.insert(DofsSideMapData{MBTET, 0, dof});
2215 ++dof;
2216 }
2217 }
2218
2220 };
2221
2222 auto set_demkowicz = [&dofs_side_map]() {
2224
2225 dofs_side_map.clear();
2226
2227 int dof = 0;
2228 for (int oo = 0; oo < Field::maxBrokenDofsOrder; oo++) {
2229
2230 // face
2231 for (auto dd = NBFACETRI_DEMKOWICZ_HDIV(oo);
2232 dd != NBFACETRI_DEMKOWICZ_HDIV(oo + 1); ++dd) {
2233 for (auto ff = 0; ff != 4; ++ff) {
2234 dofs_side_map.insert(DofsSideMapData{MBTRI, ff, dof});
2235 ++dof;
2236 }
2237 }
2238 // volume
2239 for (auto dd = NBVOLUMETET_DEMKOWICZ_HDIV(oo);
2240 dd != NBVOLUMETET_DEMKOWICZ_HDIV(oo + 1); ++dd) {
2241 dofs_side_map.insert(DofsSideMapData{MBTET, 0, dof});
2242 ++dof;
2243 }
2244 }
2245
2247 };
2248
2249 switch (continuity) {
2250 case DISCONTINUOUS:
2251 switch (base) {
2254 MoFEMFunctionReturnHot(set_ainsworth());
2256 MoFEMFunctionReturnHot(set_demkowicz());
2257 default:
2258 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
2259 }
2260 break;
2261
2262 default:
2263 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
2264 "Unknown (or not implemented) continuity");
2265 }
2266
2268}
2269
2270template <typename T>
2271auto tetCacheSwitch(const void *ptr, T &cache, std::string cache_name) {
2272 auto it = cache.find(ptr);
2273 if (it != cache.end()) {
2274 MOFEM_LOG_CHANNEL("WORLD");
2275 MOFEM_TAG_AND_LOG("WORLD", Sev::noisy, "TetPolynomialBase")
2276 << "Cache off " << cache_name << ": " << it->second.size();
2277 cache.erase(it);
2278 return false;
2279 } else {
2280 MOFEM_LOG_CHANNEL("WORLD");
2281 MOFEM_TAG_AND_LOG("WORLD", Sev::noisy, "TetPolynomialBase")
2282 << "Cache on " << cache_name;
2283 cache[ptr];
2284 return true;
2285 }
2286}
2287
2288template <>
2290 void *ptr) {
2292 std::string("hDivBaseFace") +
2294}
2295
2296template <>
2298 FieldApproximationBase base, void *ptr) {
2300 std::string("hdivBaseInterior") +
2302}
2303
2304template <>
2306 FieldApproximationBase base, void *ptr) {
2308 std::string("hdivBrokenBaseInterior") +
2310}
2311
2312template <>
2314 std::vector<void *> v) {
2315 for (auto fe_ptr : v) {
2318 }
2321 }
2324 }
2325 }
2326}
2327
2328template <>
2330 std::vector<void *> v) {
2331 for (auto fe_ptr : v) {
2334 }
2337 }
2340 }
2341 }
2342}
2343
2344template <>
2346 for (auto b = 0; b != LASTBASE; ++b) {
2348 }
2349}
2350
2351template <>
2353 for (auto b = 0; b != LASTBASE; ++b) {
2355 }
2356}
2357
2358template <>
2360 void *ptr) {
2362 std::string("hdivBaseInterior") +
2364}
2365
2366template <>
2368 std::vector<void *> v) {
2369 for (auto fe_ptr : v) {
2372 }
2373 }
2374}
2375
2376template <>
2378 std::vector<void *> v) {
2379 for (auto fe_ptr : v) {
2382 }
2383 }
2384}
2385
2386template <>
2388 for (auto b = 0; b != LASTBASE; ++b) {
2390 }
2391}
2392
2393template <>
2395 for (auto b = 0; b != LASTBASE; ++b) {
2397 }
2398}
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
#define FTENSOR_INDEX(DIM, I)
auto tetCacheSwitch(const void *ptr, T &cache, std::string cache_name)
constexpr double a
FieldApproximationBase
approximation base
Definition definitions.h:58
@ LASTBASE
Definition definitions.h:69
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
Definition definitions.h:60
@ AINSWORTH_LOBATTO_BASE
Definition definitions.h:62
@ NOBASE
Definition definitions.h:59
@ DEMKOWICZ_JACOBI_BASE
Definition definitions.h:66
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition definitions.h:64
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
FieldSpace
approximation spaces
Definition definitions.h:82
@ L2
field with C-1 continuity
Definition definitions.h:88
@ H1
continuous field
Definition definitions.h:85
@ HCURL
field with continuous tangents
Definition definitions.h:86
@ HDIV
field with continuous normal traction
Definition definitions.h:87
FieldContinuity
Field continuity.
Definition definitions.h:99
@ CONTINUOUS
Regular field.
@ DISCONTINUOUS
Broken continuity (No effect on L2 space)
static const char *const FieldSpaceNames[]
Definition definitions.h:92
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
static const char *const ApproximationBaseNames[]
Definition definitions.h:72
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
constexpr int order
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
#define NBVOLUMETET_AINSWORTH_EDGE_HDIV(P)
#define NBEDGE_DEMKOWICZ_HCURL(P)
#define NBVOLUMETET_H1(P)
Number of base functions on tetrahedron for H1 space.
#define NBVOLUMETET_AINSWORTH_HCURL(P)
#define NBFACETRI_AINSWORTH_HCURL(P)
#define NBVOLUMETET_DEMKOWICZ_HDIV(P)
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
#define NBVOLUMETET_AINSWORTH_FACE_HDIV(P)
#define NBEDGE_H1(P)
Number of base function on edge for H1 space.
#define NBFACETRI_DEMKOWICZ_HDIV(P)
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
#define NBVOLUMETET_DEMKOWICZ_HCURL(P)
#define NBFACETRI_DEMKOWICZ_HCURL(P)
#define NBFACETRI_AINSWORTH_FACE_HDIV(P)
#define NBVOLUMETET_AINSWORTH_VOLUME_HDIV(P)
#define NBFACETRI_AINSWORTH_EDGE_HDIV(P)
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
#define NBEDGE_AINSWORTH_HCURL(P)
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
#define NBFACETRI_H1(P)
Number of base function on triangle for H1 space.
#define NBVOLUMETET_L2(P)
Number of base functions on tetrahedron for L2 space.
FTensor::Index< 'i', SPACE_DIM > i
static double lambda
const double v
phase velocity of light in medium (cm/ns)
const double n
refractive index of diffusive medium
FTensor::Index< 'j', 3 > j
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasMatrix< int > MatrixInt
Definition Types.hpp:76
UBlasMatrix< double > MatrixDouble
Definition Types.hpp:77
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
FTensor::Tensor2< FTensor::PackPtr< double *, 9 >, 3, 3 > getFTensor2HVecFromPtr< 3, 3 >(double *ptr)
MoFEMErrorCode Hdiv_Demkowicz_Face_MBTET_ON_FACE(int *faces_nodes, int p, double *N, double *diffN, double *phi_f, double *diff_phi_f, int gdim, int nb)
Definition Hdiv.cpp:634
MoFEMErrorCode Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET_ON_FACE(int *faces_nodes, int p, double *N, double *diffN, double *phi_f_e[3], double *diff_phi_f_e[3], int gdim, int nb, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Hdiv base functions, Edge-based face functions by Ainsworth nme:nme847.
Definition Hdiv.cpp:47
MoFEMErrorCode Hdiv_Demkowicz_Interior_MBTET(int p, double *N, double *diffN, int p_face[], double *phi_f[4], double *diff_phi_f[4], double *phi_v, double *diff_phi_v, int gdim)
Definition Hdiv.cpp:780
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
MoFEMErrorCode Hdiv_Ainsworth_FaceBasedVolumeShapeFunctions_MBTET(int p, double *N, double *diffN, double *phi_v_f[], double *diff_phi_v_f[], int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Definition Hdiv.cpp:401
MoFEMErrorCode Hdiv_Ainsworth_FaceBubbleShapeFunctions_ON_FACE(int *faces_nodes, int p, double *N, double *diffN, double *phi_f, double *diff_phi_f, int gdim, int nb, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Face bubble functions by Ainsworth nme:nme847.
Definition Hdiv.cpp:174
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
FTensor::Tensor1< FTensor::PackPtr< double *, S >, DIM > getFTensor1FromPtr(double *ptr)
Make Tensor1 from pointer.
MoFEMErrorCode Hcurl_Demkowicz_FaceBaseFunctions_MBTET(int *faces_nodes, int *p, double *n, double *diff_n, double *phi[], double *diff_phi[], int nb_integration_pts)
Face base interior function.
Definition Hcurl.cpp:2395
MoFEMErrorCode Hdiv_Ainsworth_EdgeBasedVolumeShapeFunctions_MBTET(int p, double *N, double *diffN, double *phi_v_e[6], double *diff_phi_v_e[6], int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Hdiv base function, Edge-based interior (volume) functions by Ainsworth nme:nme847.
Definition Hdiv.cpp:307
MoFEMErrorCode Hcurl_Demkowicz_VolumeBaseFunctions_MBTET(int p, double *n, double *diff_n, double *phi, double *diff_phi, int nb_integration_pts)
Volume base interior function.
Definition Hcurl.cpp:2468
MoFEMErrorCode Hcurl_Demkowicz_EdgeBaseFunctions_MBTET(int *sense, int *p, double *n, double *diff_n, double *phi[], double *diff_phi[], int nb_integration_pts)
Edge based H-curl base functions on tetrahedral.
Definition Hcurl.cpp:2068
MoFEMErrorCode Hcurl_Ainsworth_FaceFunctions_MBTET(int *face_nodes, int *p, double *N, double *diffN, double *phi_f[4], double *diff_phi_f[4], int nb_integration_pts, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Face H-curl functions.
Definition Hcurl.cpp:1052
MoFEMErrorCode Hdiv_Ainsworth_VolumeBubbleShapeFunctions_MBTET(int p, double *N, double *diffN, double *phi_v, double *diff_phi_v, int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Interior bubble functions by Ainsworth nme:nme847.
Definition Hdiv.cpp:507
constexpr auto field_name
constexpr double g
static boost::function< int(int)> broken_nbvolumetet_edge_hdiv
Definition Hdiv.hpp:27
static boost::function< int(int)> broken_nbvolumetet_face_hdiv
Definition Hdiv.hpp:28
static boost::function< int(int)> broken_nbfacetri_face_hdiv
Definition Hdiv.hpp:26
static boost::function< int(int)> broken_nbvolumetet_volume_hdiv
Definition Hdiv.hpp:29
static boost::function< int(int)> broken_nbfacetri_edge_hdiv
Definition Hdiv.hpp:25
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.
static MoFEMErrorCode generateIndicesTriTet(const int N[], int *alpha[])
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)
static MoFEMErrorCode generateIndicesEdgeTet(const int N[], int *alpha[])
static MoFEMErrorCode generateIndicesVertexTet(const int N, int *alpha)
static MoFEMErrorCode generateIndicesTetTet(const int N, int *alpha)
Class used to pass element data to calculate base functions on tet,triangle,edge.
const FieldApproximationBase copyNodeBase
PetscErrorCode(* basePolynomialsType0)(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
const FieldContinuity spaceContinuity
const FieldApproximationBase bAse
data structure for finite element entity
MatrixInt facesNodes
nodes on finite element faces
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
static constexpr int maxBrokenDofsOrder
max number of broken dofs
Calculate base functions on tetrahedral.
ublas::matrix< MatrixDouble > diffN_face_edge
EntPolynomialBaseCtx * cTx
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)
static void switchCacheBaseOff(FieldApproximationBase base, std::vector< void * > v)
MoFEMErrorCode getValueH1AinsworthBase(MatrixDouble &pts)
MoFEMErrorCode getValueHdivAinsworthBrokenBase(MatrixDouble &pts)
ublas::vector< MatrixDouble > diffN_volume_face
MoFEMErrorCode getValueL2BernsteinBezierBase(MatrixDouble &pts)
MoFEMErrorCode getValueHcurlDemkowiczBase(MatrixDouble &pts)
static bool switchCacheBaseFace(FieldApproximationBase base, void *ptr)
MoFEMErrorCode getValueHcurl(MatrixDouble &pts)
Get base functions for Hcurl space.
MoFEMErrorCode getValueH1BernsteinBezierBase(MatrixDouble &pts)
static MoFEMErrorCode setDofsSideMap(const FieldSpace space, const FieldContinuity continuity, const FieldApproximationBase base, DofsSideMap &)
Set map of dof to side number.
MoFEMErrorCode getValueH1(MatrixDouble &pts)
Get base functions for H1 space.
ublas::vector< MatrixDouble > diffN_face_bubble
ublas::vector< MatrixDouble > N_volume_edge
ublas::vector< MatrixDouble > N_volume_face
static void switchCacheBaseOn(FieldApproximationBase base, std::vector< void * > v)
MoFEMErrorCode getValueHdiv(MatrixDouble &pts)
Get base functions for Hdiv space.
MoFEMErrorCode getValueHdivAinsworthBase(MatrixDouble &pts)
static bool switchCacheBrokenBaseInterior(FieldApproximationBase base, void *ptr)
MoFEMErrorCode getValueL2AinsworthBase(MatrixDouble &pts)
ublas::vector< MatrixDouble > N_face_bubble
ublas::matrix< MatrixDouble > N_face_edge
MoFEMErrorCode getValueL2(MatrixDouble &pts)
Get base functions for L2 space.
MoFEMErrorCode getValueHdivDemkowiczBase(MatrixDouble &pts)
TetPolynomialBase(const void *ptr=nullptr)
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
static bool switchCacheBaseInterior(FieldApproximationBase base, void *ptr)
MoFEMErrorCode getValueHdivDemkowiczBrokenBase(MatrixDouble &pts)
ublas::vector< MatrixDouble > diffN_volume_edge
MoFEMErrorCode getValueHcurlAinsworthBase(MatrixDouble &pts)
static MoFEMErrorCode setDofsSideMapHdiv(const FieldContinuity continuity, const FieldApproximationBase base, DofsSideMap &dofs_side_map)
Set the Dofs Side Map Hdiv object.
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:747
static constexpr std::array< double, 12 > diffShapeFunMBTET
Definition Tools.hpp:271
base class for all interface classes
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
static std::array< std::map< const void *, BaseCacheMI >, LASTBASE > hdivBaseInterior
static std::array< std::map< const void *, HDivBaseFaceCacheMI >, LASTBASE > hDivBaseFace
static std::array< std::map< const void *, BaseCacheMI >, LASTBASE > l2BaseInterior
static std::array< std::map< const void *, BaseCacheMI >, LASTBASE > hdivBrokenBaseInterior
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
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