v0.13.2
Loading...
Searching...
No Matches
ForcesAndSourcesCore.cpp
Go to the documentation of this file.
1/** \file ForcesAndSourcesCore.cpp
2
3\brief Implementation of Elements on Entities for Forces and Sources
4*/
5
6#ifdef __cplusplus
7extern "C" {
8#endif
9#include <cblas.h>
10#include <lapack_wrap.h>
11// #include <gm_rule.h>
12#ifdef __cplusplus
13}
14#endif
15
16namespace MoFEM {
17
18static auto cmp_uid_lo(const boost::weak_ptr<FieldEntity> &a, const UId &b) {
19 if (auto a_ptr = a.lock()) {
20 if (a_ptr->getLocalUniqueId() < b)
21 return true;
22 else
23 return false;
24 } else {
25 return false;
26 }
27}
28
29static auto cmp_uid_hi(const UId &b, const boost::weak_ptr<FieldEntity> &a) {
30 if (auto a_ptr = a.lock()) {
31 if (b < a_ptr->getLocalUniqueId())
32 return true;
33 else
34 return false;
35 } else {
36 return true;
37 }
38}
39
41 :
42
43 mField(m_field), getRuleHook(0), setRuleHook(0),
44 dataOnElement{
45
46 boost::make_shared<EntitiesFieldData>(MBENTITYSET), // NOSPACE,
47 boost::make_shared<EntitiesFieldData>(MBENTITYSET), // NOFIELD
48 boost::make_shared<EntitiesFieldData>(MBENTITYSET), // H1
49 boost::make_shared<EntitiesFieldData>(MBENTITYSET), // HCURL
50 boost::make_shared<EntitiesFieldData>(MBENTITYSET), // HDIV
51 boost::make_shared<EntitiesFieldData>(MBENTITYSET) // L2
52
53 },
54 derivedDataOnElement{
55
56 nullptr,
57 boost::make_shared<DerivedEntitiesFieldData>(
58 dataOnElement[NOFIELD]), // NOFIELD
59 boost::make_shared<DerivedEntitiesFieldData>(dataOnElement[H1]), // H1
60 boost::make_shared<DerivedEntitiesFieldData>(
61 dataOnElement[HCURL]), // HCURL
62 boost::make_shared<DerivedEntitiesFieldData>(
63 dataOnElement[HDIV]), // HDIV
64 boost::make_shared<DerivedEntitiesFieldData>(dataOnElement[L2]) // L2
65
66 },
67 dataNoField(*dataOnElement[NOFIELD].get()),
68 dataH1(*dataOnElement[H1].get()), dataHcurl(*dataOnElement[HCURL].get()),
69 dataHdiv(*dataOnElement[HDIV].get()), dataL2(*dataOnElement[L2].get()),
70 lastEvaluatedElementEntityType(MBMAXTYPE), sidePtrFE(nullptr),
71 refinePtrFE(nullptr) {
72
73 dataOnElement[NOSPACE]->dataOnEntities[MBENTITYSET].push_back(
75
76 dataOnElement[NOFIELD]->dataOnEntities[MBENTITYSET].push_back(
78 derivedDataOnElement[NOFIELD]->dataOnEntities[MBENTITYSET].push_back(
80}
81
82// ** Sense **
83
85 const EntityType type,
86 boost::ptr_vector<EntitiesFieldData::EntData> &data) const {
88
89 auto &side_table = numeredEntFiniteElementPtr->getSideNumberTable().get<0>();
90 auto sit = side_table.lower_bound(get_id_for_min_type(type));
91 if (sit != side_table.end()) {
92 auto hi_sit = side_table.upper_bound(get_id_for_max_type(type));
93 for (; sit != hi_sit; ++sit) {
94 if (const auto side_number = (*sit)->side_number; side_number >= 0) {
95 const int brother_side_number = (*sit)->brother_side_number;
96 const int sense = (*sit)->sense;
97
98 data[side_number].getSense() = sense;
99 if (brother_side_number != -1)
100 data[brother_side_number].getSense() = sense;
101 }
102 }
103 }
105}
106
107// ** Order **
108
109template <typename ENTMULTIINDEX>
110static inline int getMaxOrder(const ENTMULTIINDEX &multi_index) {
111 int max_order = 0;
112 for (auto ent_field_weak_ptr : multi_index)
113 if (auto e = ent_field_weak_ptr.lock()) {
114 const int order = e->getMaxOrder();
115 max_order = (max_order < order) ? order : max_order;
116 }
117 return max_order;
118}
119
121 int max_order = 0;
122 for (auto e : getDataFieldEnts()) {
123 if (auto ptr = e.lock()) {
124 const int order = ptr->getMaxOrder();
125 max_order = (max_order < order) ? order : max_order;
126 }
127 }
128 return max_order;
129}
130
133}
134
137}
138
140 const EntityType type, const FieldSpace space,
141 boost::ptr_vector<EntitiesFieldData::EntData> &data) const {
143
144 auto set_order = [&]() {
146 auto &side_table = numeredEntFiniteElementPtr->getSideNumberTable();
147
148 for (unsigned int s = 0; s != data.size(); ++s)
149 data[s].getOrder() = 0;
150
151 const FieldEntity_vector_view &data_field_ent = getDataFieldEnts();
152
153 for (auto r = fieldsPtr->get<BitFieldId_space_mi_tag>().equal_range(space);
154 r.first != r.second; ++r.first) {
155
156 const auto field_bit_number = (*r.first)->getBitNumber();
157 const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
158 field_bit_number, get_id_for_min_type(type));
159 auto lo = std::lower_bound(data_field_ent.begin(), data_field_ent.end(),
160 lo_uid, cmp_uid_lo);
161 if (lo != data_field_ent.end()) {
162 const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
163 field_bit_number, get_id_for_max_type(type));
164 auto hi =
165 std::upper_bound(lo, data_field_ent.end(), hi_uid, cmp_uid_hi);
166 for (; lo != hi; ++lo) {
167
168 if (auto ptr = lo->lock()) {
169
170 auto &e = *ptr;
171 auto sit = side_table.find(e.getEnt());
172 if (sit != side_table.end()) {
173 auto &side = *sit;
174 if (const auto side_number = side->side_number;
175 side_number >= 0) {
176 ApproximationOrder ent_order = e.getMaxOrder();
177 auto &dat = data[side_number];
178 dat.getOrder() =
179 dat.getOrder() > ent_order ? dat.getOrder() : ent_order;
180 }
181 } else
182 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
183 "Entity on side of the element not found");
184 }
185 }
186 }
187 }
188
190 };
191
192 auto set_order_on_brother = [&]() {
194 auto &side_table =
195 numeredEntFiniteElementPtr->getSideNumberTable().get<0>();
196 auto sit = side_table.lower_bound(get_id_for_min_type(type));
197 if (sit != side_table.end()) {
198 auto hi_sit = side_table.upper_bound(get_id_for_max_type(type));
199 for (; sit != hi_sit; ++sit) {
200 const int brother_side_number = (*sit)->brother_side_number;
201 if (brother_side_number != -1) {
202 const int side_number = (*sit)->side_number;
203 data[brother_side_number].getOrder() = data[side_number].getOrder();
204 }
205 }
206 }
208 };
209
210 CHKERR set_order();
211 CHKERR set_order_on_brother();
212
214}
215
216// ** Indices **
217
218template <typename EXTRACTOR>
220 const int bit_number, FieldEntity_vector_view &ents_field,
221 VectorInt &nodes_indices, VectorInt &local_nodes_indices,
222 EXTRACTOR &&extractor) const {
224
225 auto field_it = fieldsPtr->get<BitFieldId_mi_tag>().find(
226 BitFieldId().set(bit_number - 1));
227 if (field_it != fieldsPtr->get<BitFieldId_mi_tag>().end()) {
228
229#ifndef NDEBUG
230 if ((*field_it)->getBitNumber() != bit_number)
231 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong bit number");
232#endif
233 const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
234 bit_number, get_id_for_min_type<MBVERTEX>());
235 auto lo = std::lower_bound(ents_field.begin(), ents_field.end(), lo_uid,
236 cmp_uid_lo);
237 if (lo != ents_field.end()) {
238 const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
239 bit_number, get_id_for_max_type<MBVERTEX>());
240 auto hi = std::upper_bound(lo, ents_field.end(), hi_uid, cmp_uid_hi);
241
242 const int num_nodes = getNumberOfNodes();
243 const int nb_dofs_on_vert = (*field_it)->getNbOfCoeffs();
244 const int max_nb_dofs = nb_dofs_on_vert * num_nodes;
245
246 int nb_dofs = 0;
247 for (auto it = lo; it != hi; ++it) {
248 if (auto e = it->lock()) {
249 if (auto cache = extractor(e).lock()) {
250 if (cache->loHi[0] != cache->loHi[1]) {
251 nb_dofs += std::distance(cache->loHi[0], cache->loHi[1]);
252 break;
253 }
254 }
255 }
256 }
257
258 if (nb_dofs) {
259 nodes_indices.resize(max_nb_dofs, false);
260 local_nodes_indices.resize(max_nb_dofs, false);
261 } else {
262 nodes_indices.resize(0, false);
263 local_nodes_indices.resize(0, false);
264 }
265
266 if (nb_dofs != max_nb_dofs) {
267 std::fill(nodes_indices.begin(), nodes_indices.end(), -1);
268 std::fill(local_nodes_indices.begin(), local_nodes_indices.end(), -1);
269 }
270
271 for (auto it = lo; it != hi; ++it) {
272 if (auto e = it->lock()) {
273 auto side_ptr = e->getSideNumberPtr();
274 if (const auto side_number = side_ptr->side_number;
275 side_number >= 0) {
276 const auto brother_side_number = side_ptr->brother_side_number;
277 if (auto cache = extractor(e).lock()) {
278 for (auto dit = cache->loHi[0]; dit != cache->loHi[1]; ++dit) {
279 auto &dof = **dit;
280 const int idx = dof.getPetscGlobalDofIdx();
281 const int local_idx = dof.getPetscLocalDofIdx();
282 const int pos =
283 side_number * nb_dofs_on_vert + dof.getDofCoeffIdx();
284 nodes_indices[pos] = idx;
285 local_nodes_indices[pos] = local_idx;
286 if (brother_side_number != -1) {
287 const int elem_idx = brother_side_number * nb_dofs_on_vert +
288 (*dit)->getDofCoeffIdx();
289 nodes_indices[elem_idx] = idx;
290 local_nodes_indices[elem_idx] = local_idx;
291 }
292 }
293 }
294 }
295 }
296 }
297 } else {
298 nodes_indices.resize(0, false);
299 local_nodes_indices.resize(0, false);
300 }
301
302 } else {
303 nodes_indices.resize(0, false);
304 local_nodes_indices.resize(0, false);
305 }
306
308}
309
312 const int bit_number) const {
313
314 struct Extractor {
315 boost::weak_ptr<EntityCacheNumeredDofs>
316 operator()(boost::shared_ptr<FieldEntity> &e) {
317 return e->entityCacheRowDofs;
318 }
319 };
320
321 return getNodesIndices(bit_number, getRowFieldEnts(),
322 data.dataOnEntities[MBVERTEX][0].getIndices(),
323 data.dataOnEntities[MBVERTEX][0].getLocalIndices(),
324 Extractor());
325}
326
329 const int bit_number) const {
330
331 struct Extractor {
332 boost::weak_ptr<EntityCacheNumeredDofs>
333 operator()(boost::shared_ptr<FieldEntity> &e) {
334 return e->entityCacheColDofs;
335 }
336 };
337
338 return getNodesIndices(bit_number, getColFieldEnts(),
339 data.dataOnEntities[MBVERTEX][0].getIndices(),
340 data.dataOnEntities[MBVERTEX][0].getLocalIndices(),
341 Extractor());
342}
343
344template <typename EXTRACTOR>
346 EntitiesFieldData &data, const int bit_number,
347 FieldEntity_vector_view &ents_field, const EntityType type_lo,
348 const EntityType type_hi, EXTRACTOR &&extractor) const {
350
351 for (EntityType t = type_lo; t != type_hi; ++t) {
352 for (auto &dat : data.dataOnEntities[t]) {
353 dat.getIndices().resize(0, false);
354 dat.getLocalIndices().resize(0, false);
355 }
356 }
357
358 const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
359 bit_number, get_id_for_min_type(type_lo));
360 auto lo = std::lower_bound(ents_field.begin(), ents_field.end(), lo_uid,
361 cmp_uid_lo);
362 if (lo != ents_field.end()) {
363 const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
364 bit_number, get_id_for_max_type(type_hi));
365 auto hi = std::upper_bound(lo, ents_field.end(), hi_uid, cmp_uid_hi);
366 if (lo != hi) {
367
368 std::vector<boost::weak_ptr<FieldEntity>> brother_ents_vec;
369
370 for (auto it = lo; it != hi; ++it)
371 if (auto e = it->lock()) {
372
373 const EntityType type = e->getEntType();
374 auto side_ptr = e->getSideNumberPtr();
375 if (const auto side = side_ptr->side_number; side >= 0) {
376 const auto nb_dofs_on_ent = e->getNbDofsOnEnt();
377 const auto brother_side = side_ptr->brother_side_number;
378 auto &dat = data.dataOnEntities[type][side];
379 auto &ent_field_indices = dat.getIndices();
380 auto &ent_field_local_indices = dat.getLocalIndices();
381
382 ent_field_indices.resize(nb_dofs_on_ent, false);
383 ent_field_local_indices.resize(nb_dofs_on_ent, false);
384 std::fill(ent_field_indices.data().begin(),
385 ent_field_indices.data().end(), -1);
386 std::fill(ent_field_local_indices.data().begin(),
387 ent_field_local_indices.data().end(), -1);
388
389 if (auto cache = extractor(e).lock()) {
390 for (auto dit = cache->loHi[0]; dit != cache->loHi[1]; ++dit) {
391 const int idx = (*dit)->getEntDofIdx();
392 ent_field_indices[idx] = (*dit)->getPetscGlobalDofIdx();
393 ent_field_local_indices[idx] = (*dit)->getPetscLocalDofIdx();
394 }
395 }
396
397 if (brother_side != -1) {
398 auto &dat_brother = data.dataOnEntities[type][brother_side];
399 dat_brother.getIndices().resize(nb_dofs_on_ent, false);
400 dat_brother.getLocalIndices().resize(nb_dofs_on_ent, false);
401 noalias(dat_brother.getIndices()) = dat.getIndices();
402 noalias(dat_brother.getLocalIndices()) = dat.getLocalIndices();
403 }
404 }
405 }
406 }
407 }
408
410}
411
413 EntitiesFieldData &data, const int bit_number, const EntityType type_lo,
414 const EntityType type_hi) const {
415
416 struct Extractor {
417 boost::weak_ptr<EntityCacheNumeredDofs>
418 operator()(boost::shared_ptr<FieldEntity> &e) {
419 return e->entityCacheRowDofs;
420 }
421 };
422
423 return getEntityIndices(data, bit_number, getRowFieldEnts(), type_lo, type_hi,
424 Extractor());
425}
426
428 EntitiesFieldData &data, const int bit_number, const EntityType type_lo,
429 const EntityType type_hi) const {
430
431 struct Extractor {
432 boost::weak_ptr<EntityCacheNumeredDofs>
433 operator()(boost::shared_ptr<FieldEntity> &e) {
434 return e->entityCacheColDofs;
435 }
436 };
437
438 return getEntityIndices(data, bit_number, getColFieldEnts(), type_lo, type_hi,
439 Extractor());
440}
441
443 const int bit_number, boost::shared_ptr<FENumeredDofEntity_multiIndex> dofs,
444 VectorInt &indices) const {
446 auto dit = dofs->get<Unique_mi_tag>().lower_bound(
448 auto hi_dit = dofs->get<Unique_mi_tag>().upper_bound(
450 indices.resize(std::distance(dit, hi_dit));
451 for (; dit != hi_dit; dit++) {
452 int idx = (*dit)->getPetscGlobalDofIdx();
453 indices[(*dit)->getDofCoeffIdx()] = idx;
454 }
456}
457
460 const int bit_number) const {
462#ifndef NDEBUG
463 if (data.dataOnEntities[MBENTITYSET].size() == 0) {
464 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
465 "data.dataOnEntities[MBENTITYSET] is empty");
466 }
467#endif
469 data.dataOnEntities[MBENTITYSET][0].getIndices());
471}
472
475 const int bit_number) const {
477#ifndef NDEBUG
478 if (data.dataOnEntities[MBENTITYSET].size() == 0) {
479 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
480 "data.dataOnEntities[MBENTITYSET] is empty");
481 }
482#endif
484 data.dataOnEntities[MBENTITYSET][0].getIndices());
486}
487
488// ** Indices from problem **
489
491 const std::string &field_name, const NumeredDofEntity_multiIndex &dofs,
492 VectorInt &nodes_indices) const {
494
495 const Field *field_struture = mField.get_field_structure(field_name);
496 if (field_struture->getSpace() == H1) {
497
498 const int num_nodes = getNumberOfNodes();
499 nodes_indices.resize(field_struture->getNbOfCoeffs() * num_nodes, false);
500 std::fill(nodes_indices.begin(), nodes_indices.end(), -1);
501
502 auto &side_table = const_cast<SideNumber_multiIndex &>(
503 numeredEntFiniteElementPtr->getSideNumberTable());
504 auto siit = side_table.get<1>().lower_bound(boost::make_tuple(MBVERTEX, 0));
505 auto hi_siit =
506 side_table.get<1>().lower_bound(boost::make_tuple(MBVERTEX, 10000));
507
508 int nn = 0;
509 for (; siit != hi_siit; siit++, nn++) {
510 if (siit->get()->side_number >= 0) {
511 auto bit_number = mField.get_field_bit_number(field_name);
512 const EntityHandle ent = siit->get()->ent;
513 auto dit = dofs.get<Unique_mi_tag>().lower_bound(
515 auto hi_dit = dofs.get<Unique_mi_tag>().upper_bound(
517 for (; dit != hi_dit; dit++) {
518 nodes_indices[siit->get()->side_number * (*dit)->getNbOfCoeffs() +
519 (*dit)->getDofCoeffIdx()] =
520 (*dit)->getPetscGlobalDofIdx();
521 }
522 }
523 }
524 } else {
525 nodes_indices.resize(0, false);
526 }
527
529}
530
532 const std::string &field_name, const NumeredDofEntity_multiIndex &dofs,
533 EntityType type, int side_number, VectorInt &indices) const {
535
536 indices.resize(0);
537
538 auto &side_table = const_cast<SideNumber_multiIndex &>(
539 numeredEntFiniteElementPtr->getSideNumberTable());
540 auto siit =
541 side_table.get<1>().lower_bound(boost::make_tuple(type, side_number));
542 auto hi_siit =
543 side_table.get<1>().upper_bound(boost::make_tuple(type, side_number));
544
545 for (; siit != hi_siit; siit++) {
546 if (siit->get()->side_number >= 0) {
547
548 const EntityHandle ent = siit->get()->ent;
549 auto bit_number = mField.get_field_bit_number(field_name);
550 auto dit = dofs.get<Unique_mi_tag>().lower_bound(
552 auto hi_dit = dofs.get<Unique_mi_tag>().upper_bound(
554 indices.resize(std::distance(dit, hi_dit));
555 for (; dit != hi_dit; dit++) {
556
557 indices[(*dit)->getEntDofIdx()] = (*dit)->getPetscGlobalDofIdx();
558 }
559 }
560 }
561
563}
564
566 const std::string &field_name, VectorInt &nodes_indices) const {
568 nodes_indices);
569}
570
573 EntityType type, int side_number,
574 VectorInt &indices) const {
576 type, side_number, indices);
577}
578
580 const std::string &field_name, VectorInt &nodes_indices) const {
582 nodes_indices);
583}
584
587 EntityType type, int side_number,
588 VectorInt &indices) const {
590 type, side_number, indices);
591}
592
593// ** Data **
594
597
598 // for (int s = H1; s != LASTSPACE; ++s) {
599 // dataOnElement[s]->dataOnEntities[MBENTITYSET].resize(0);
600 // }
601
602 for (auto &data : dataOnElement) {
603 if (data) {
604 for (auto &dat : data->dataOnEntities) {
605 for (auto &ent_dat : dat) {
606 ent_dat.getEntDataBitRefLevel().reset();
607 }
608 }
609 }
610 }
611
612 auto &field_ents = getDataFieldEnts();
613 for (auto it : field_ents) {
614 if (auto e = it.lock()) {
615 const FieldSpace space = e->getSpace();
616 if (space > NOFIELD) {
617 const EntityType type = e->getEntType();
618 const signed char side =
619 type == MBVERTEX ? 0 : e->getSideNumberPtr()->side_number;
620 if (side >= 0) {
621 if (auto &data = dataOnElement[space]) {
622 auto &dat = data->dataOnEntities[type][side];
623 dat.getEntDataBitRefLevel() |= e->getBitRefLevel();
624 }
625 }
626 } else {
627 if (auto &data = dataOnElement[NOFIELD]) {
628 auto &dat = data->dataOnEntities[MBENTITYSET][0];
629 dat.getEntDataBitRefLevel() |= e->getBitRefLevel();
630 }
631 }
632 }
633 }
634
636};
637
640 const int bit_number) const {
641
642 auto get_nodes_field_data = [&](VectorDouble &nodes_data,
643 VectorFieldEntities &field_entities,
644 VectorDofs &nodes_dofs, FieldSpace &space,
646 VectorInt &bb_node_order) {
648
649 nodes_data.resize(0, false);
650 nodes_dofs.resize(0, false);
651 field_entities.resize(0, false);
652
653 auto field_it =
654 fieldsPtr->get<BitFieldId_mi_tag>().find(BitFEId().set(bit_number - 1));
655 if (field_it != fieldsPtr->get<BitFieldId_mi_tag>().end()) {
656
657#ifndef NDEBUG
658 if ((*field_it)->getBitNumber() != bit_number)
659 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong bit number");
660#endif
661 const int nb_dofs_on_vert = (*field_it)->getNbOfCoeffs();
662 space = (*field_it)->getSpace();
663 base = (*field_it)->getApproxBase();
664
665 auto &field_ents = getDataFieldEnts();
666 const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
667 bit_number, get_id_for_min_type<MBVERTEX>());
668 auto lo = std::lower_bound(field_ents.begin(), field_ents.end(), lo_uid,
669 cmp_uid_lo);
670 if (lo != field_ents.end()) {
671 const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
672 bit_number, get_id_for_max_type<MBVERTEX>());
673 auto hi = std::upper_bound(lo, field_ents.end(), hi_uid, cmp_uid_hi);
674 if (lo != hi) {
675
676 int nb_dofs = 0;
677 for (auto it = lo; it != hi; ++it) {
678 if (auto e = it->lock()) {
679 nb_dofs += e->getNbDofsOnEnt();
680 }
681 }
682
683 if (nb_dofs) {
684
685 const int num_nodes = getNumberOfNodes();
686 bb_node_order.resize(num_nodes, false);
687 bb_node_order.clear();
688 const int max_nb_dofs = nb_dofs_on_vert * num_nodes;
689 nodes_data.resize(max_nb_dofs, false);
690 nodes_dofs.resize(max_nb_dofs, false);
691 field_entities.resize(num_nodes, false);
692 std::fill(nodes_data.begin(), nodes_data.end(), 0);
693 std::fill(nodes_dofs.begin(), nodes_dofs.end(), nullptr);
694 std::fill(field_entities.begin(), field_entities.end(), nullptr);
695
696 std::vector<boost::weak_ptr<FieldEntity>> brother_ents_vec;
697
698 for (auto it = lo; it != hi; ++it) {
699 if (auto e = it->lock()) {
700 const auto &sn = e->getSideNumberPtr();
701 // Some field entities on skeleton can have negative side
702 // number
703 if (const auto side_number = sn->side_number;
704 side_number >= 0) {
705 const int brother_side_number = sn->brother_side_number;
706
707 field_entities[side_number] = e.get();
708 if (brother_side_number != -1) {
709 brother_ents_vec.emplace_back(e);
710 field_entities[side_number] = field_entities[side_number];
711 }
712
713 bb_node_order[side_number] = e->getMaxOrder();
714 int pos = side_number * nb_dofs_on_vert;
715 auto ent_filed_data_vec = e->getEntFieldData();
716 if (auto cache = e->entityCacheDataDofs.lock()) {
717 for (auto dit = cache->loHi[0]; dit != cache->loHi[1];
718 ++dit) {
719 const auto dof_idx = (*dit)->getEntDofIdx();
720 nodes_data[pos + dof_idx] = ent_filed_data_vec[dof_idx];
721 nodes_dofs[pos + dof_idx] =
722 reinterpret_cast<FEDofEntity *>((*dit).get());
723 }
724 }
725 }
726 }
727 }
728
729 for (auto &it : brother_ents_vec) {
730 if (const auto e = it.lock()) {
731 const auto &sn = e->getSideNumberPtr();
732 const int side_number = sn->side_number;
733 const int brother_side_number = sn->brother_side_number;
734 bb_node_order[brother_side_number] = bb_node_order[side_number];
735 int pos = side_number * nb_dofs_on_vert;
736 int brother_pos = brother_side_number * nb_dofs_on_vert;
737 for (int ii = 0; ii != nb_dofs_on_vert; ++ii) {
738 nodes_data[brother_pos] = nodes_data[pos];
739 nodes_dofs[brother_pos] = nodes_dofs[pos];
740 ++pos;
741 ++brother_pos;
742 }
743 }
744 }
745 }
746 }
747 }
748 }
749
751 };
752
753 return get_nodes_field_data(
754 data.dataOnEntities[MBVERTEX][0].getFieldData(),
755 data.dataOnEntities[MBVERTEX][0].getFieldEntities(),
756 data.dataOnEntities[MBVERTEX][0].getFieldDofs(),
757 data.dataOnEntities[MBVERTEX][0].getSpace(),
758 data.dataOnEntities[MBVERTEX][0].getBase(),
759 data.dataOnEntities[MBVERTEX][0].getBBNodeOrder());
760}
761
763 EntitiesFieldData &data, const int bit_number, const EntityType type_lo,
764 const EntityType type_hi) const {
766 for (EntityType t = type_lo; t != type_hi; ++t) {
767 for (auto &dat : data.dataOnEntities[t]) {
768 dat.getOrder() = 0;
769 dat.getBase() = NOBASE;
770 dat.getSpace() = NOSPACE;
771 dat.getFieldData().resize(0, false);
772 dat.getFieldDofs().resize(0, false);
773 dat.getFieldEntities().resize(0, false);
774 }
775 }
776
777 auto &field_ents = getDataFieldEnts();
778 const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
779 bit_number, get_id_for_min_type(type_lo));
780 auto lo = std::lower_bound(field_ents.begin(), field_ents.end(), lo_uid,
781 cmp_uid_lo);
782 if (lo != field_ents.end()) {
783 const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
784 bit_number, get_id_for_max_type(type_hi));
785 auto hi = std::upper_bound(lo, field_ents.end(), hi_uid, cmp_uid_hi);
786 if (lo != hi) {
787
788 std::vector<boost::weak_ptr<FieldEntity>> brother_ents_vec;
789
790 for (auto it = lo; it != hi; ++it)
791 if (auto e = it->lock()) {
792 auto side_ptr = e->getSideNumberPtr();
793 if (const auto side = side_ptr->side_number; side >= 0) {
794 const EntityType type = e->getEntType();
795 auto &dat = data.dataOnEntities[type][side];
796 auto &ent_field = dat.getFieldEntities();
797 auto &ent_field_dofs = dat.getFieldDofs();
798 auto &ent_field_data = dat.getFieldData();
799
800 const int brother_side = side_ptr->brother_side_number;
801 if (brother_side != -1)
802 brother_ents_vec.emplace_back(e);
803
804 dat.getBase() = e->getApproxBase();
805 dat.getSpace() = e->getSpace();
806 const int ent_order = e->getMaxOrder();
807 dat.getOrder() =
808 dat.getOrder() > ent_order ? dat.getOrder() : ent_order;
809
810 auto ent_data = e->getEntFieldData();
811 ent_field_data.resize(ent_data.size(), false);
812 noalias(ent_field_data) = ent_data;
813 ent_field_dofs.resize(ent_data.size(), false);
814 std::fill(ent_field_dofs.begin(), ent_field_dofs.end(), nullptr);
815 ent_field.resize(1, false);
816 ent_field[0] = e.get();
817 if (auto cache = e->entityCacheDataDofs.lock()) {
818 for (auto dit = cache->loHi[0]; dit != cache->loHi[1]; ++dit) {
819 ent_field_dofs[(*dit)->getEntDofIdx()] =
820 reinterpret_cast<FEDofEntity *>((*dit).get());
821 }
822 }
823 }
824 }
825
826 for (auto &it : brother_ents_vec) {
827 if (const auto e = it.lock()) {
828 const EntityType type = e->getEntType();
829 const int side = e->getSideNumberPtr()->side_number;
830 const int brother_side = e->getSideNumberPtr()->brother_side_number;
831 auto &dat = data.dataOnEntities[type][side];
832 auto &dat_brother = data.dataOnEntities[type][brother_side];
833 dat_brother.getBase() = dat.getBase();
834 dat_brother.getSpace() = dat.getSpace();
835 dat_brother.getOrder() = dat.getOrder();
836 dat_brother.getFieldData() = dat.getFieldData();
837 dat_brother.getFieldDofs() = dat.getFieldDofs();
838 dat_brother.getFieldEntities() = dat.getFieldEntities();
839 }
840 }
841 }
842 }
843
845}
846
848 const int bit_number, VectorDouble &ent_field_data,
849 VectorDofs &ent_field_dofs, VectorFieldEntities &ent_field) const {
850
852
853 ent_field_data.resize(0, false);
854 ent_field_dofs.resize(0, false);
855 ent_field.resize(0, false);
856
857 auto &field_ents = getDataFieldEnts();
858 const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
859 bit_number, get_id_for_min_type<MBVERTEX>());
860 auto lo = std::lower_bound(field_ents.begin(), field_ents.end(), lo_uid,
861 cmp_uid_lo);
862 if (lo != field_ents.end()) {
863
864 ent_field.resize(field_ents.size(), false);
865 std::fill(ent_field.begin(), ent_field.end(), nullptr);
866
867 const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
868 bit_number, get_id_for_max_type<MBENTITYSET>());
869 auto hi = std::upper_bound(lo, field_ents.end(), hi_uid, cmp_uid_hi);
870 if (lo != hi) {
871
872 int side = 0;
873 for (auto it = lo; it != hi; ++it, ++side)
874 if (auto e = it->lock()) {
875
876 const auto size = e->getNbDofsOnEnt();
877 ent_field_data.resize(size, false);
878 ent_field_dofs.resize(size, false);
879 ent_field[side] = e.get();
880 noalias(ent_field_data) = e->getEntFieldData();
881
882 if (auto cache = e->entityCacheDataDofs.lock()) {
883 for (auto dit = cache->loHi[0]; dit != cache->loHi[1]; ++dit) {
884 ent_field_dofs[(*dit)->getEntDofIdx()] =
885 reinterpret_cast<FEDofEntity *>((*dit).get());
886 }
887 }
888 }
889 }
890 }
891
893}
894
897 const int bit_number) const {
899 if (data.dataOnEntities[MBENTITYSET].size() == 0)
900 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
901 "No space to insert data");
902
904 bit_number, data.dataOnEntities[MBENTITYSET][0].getFieldData(),
905 data.dataOnEntities[MBENTITYSET][0].getFieldDofs(),
906 data.dataOnEntities[MBENTITYSET][0].getFieldEntities());
908}
909
910// ** Face **
911
915 auto &face_nodes = data.facesNodes;
916 auto &face_nodes_order = data.facesNodesOrder;
917 auto &side_table = const_cast<SideNumber_multiIndex &>(
918 numeredEntFiniteElementPtr->getSideNumberTable());
919 const auto ent = numeredEntFiniteElementPtr->getEnt();
920 const auto type = numeredEntFiniteElementPtr->getEntType();
921 const auto nb_faces = CN::NumSubEntities(type, 2);
922 const EntityHandle *conn_ele;
923 int num_nodes_ele;
924 CHKERR mField.get_moab().get_connectivity(ent, conn_ele, num_nodes_ele, true);
925 auto side_ptr_it = side_table.get<1>().lower_bound(
926 boost::make_tuple(CN::TypeDimensionMap[2].first, 0));
927 auto hi_side_ptr_it = side_table.get<1>().upper_bound(
928 boost::make_tuple(CN::TypeDimensionMap[2].second, 100));
929
930 for (; side_ptr_it != hi_side_ptr_it; ++side_ptr_it) {
931 const auto side = (*side_ptr_it)->side_number;
932 const auto sense = (*side_ptr_it)->sense;
933 const auto offset = (*side_ptr_it)->offset;
934 const auto face_ent = (*side_ptr_it)->ent;
935
936 EntityType face_type;
937 int nb_nodes_face;
938 auto face_indices =
939 CN::SubEntityVertexIndices(type, 2, side, face_type, nb_nodes_face);
940 face_nodes.resize(nb_faces, nb_nodes_face);
941 face_nodes_order.resize(nb_faces, nb_nodes_face);
942
943 if (sense == 1)
944 for (int n = 0; n != nb_nodes_face; ++n)
945 face_nodes_order(side, n) = (n + offset) % nb_nodes_face;
946 else
947 for (int n = 0; n != nb_nodes_face; ++n)
948 face_nodes_order(side, n) =
949 (nb_nodes_face - (n - offset) % nb_nodes_face) % nb_nodes_face;
950
951 for (int n = 0; n != nb_nodes_face; ++n)
952 face_nodes(side, n) = face_indices[face_nodes_order(side, n)];
953
954#ifndef NDEBUG
955 auto check = [&]() {
957 const EntityHandle *conn_face;
958 // int nb_nodes_face;
959 CHKERR mField.get_moab().get_connectivity(face_ent, conn_face,
960 nb_nodes_face, true);
961 face_nodes.resize(nb_faces, nb_nodes_face);
962 for (int nn = 0; nn != nb_nodes_face; ++nn) {
963 if (face_nodes(side, nn) !=
964 std::distance(
965 conn_ele,
966 std::find(conn_ele, &conn_ele[num_nodes_ele], conn_face[nn]))) {
967 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
968 "Wrong face numeration");
969 }
970 }
972 };
973 CHKERR check();
974#endif
975 }
976
978}
979
980// ** Space and Base **
981
983 EntitiesFieldData &data) const {
985
986 if (nInTheLoop == 0) {
987 data.sPace.reset();
988 data.bAse.reset();
989 for (EntityType t = MBVERTEX; t != MBMAXTYPE; ++t) {
990 data.spacesOnEntities[t].reset();
991 data.basesOnEntities[t].reset();
992 }
993 for (int s = 0; s != LASTSPACE; ++s) {
994 data.basesOnSpaces[s].reset();
995 }
996 }
997
998 auto fe_type = numeredEntFiniteElementPtr->getEntType();
999
1000 if (getDataFieldEntsPtr())
1001 for (auto e : getDataFieldEnts()) {
1002 if (auto ptr = e.lock()) {
1003 // get data from entity
1004 const EntityType type = ptr->getEntType();
1005 if (DefaultElementAdjacency::getDefTypeMap(fe_type, type)) {
1006 const FieldSpace space = ptr->getSpace();
1007 const FieldApproximationBase approx = ptr->getApproxBase();
1008 // set data
1009 data.sPace.set(space);
1010 data.bAse.set(approx);
1011 data.spacesOnEntities[type].set(space);
1012 data.basesOnEntities[type].set(approx);
1013 data.basesOnSpaces[space].set(approx);
1014 }
1015 }
1016 }
1017 else
1018 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1019 "data fields ents not allocated on element");
1020
1022}
1023
1025 const FieldApproximationBase b) {
1027 if (dataOnElement[H1]->bAse.test(b)) {
1028 switch (static_cast<FieldApproximationBase>(b)) {
1029 case NOBASE:
1030 break;
1032 // BERNSTEIN_BEZIER_BASE is not hierarchical base
1033 break;
1038 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1039 "Functions genrating approximation base not defined");
1040
1041 for (int space = H1; space != LASTSPACE; ++space) {
1042 if (dataOnElement[H1]->sPace.test(space) &&
1043 dataOnElement[H1]->bAse.test(b) &&
1044 dataOnElement[H1]->basesOnSpaces[space].test(b)) {
1045 CHKERR getElementPolynomialBase()->getValue(
1046 gaussPts,
1047 boost::make_shared<EntPolynomialBaseCtx>(
1048 *dataOnElement[space], static_cast<FieldSpace>(space),
1049 static_cast<FieldApproximationBase>(b), NOBASE));
1050 }
1051 }
1052 break;
1053 case USER_BASE:
1054 if (!getUserPolynomialBase())
1055 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1056 "Functions generating user approximation base not defined");
1057
1058 for (int space = H1; space != LASTSPACE; ++space)
1059 if (dataOnElement[H1]->sPace.test(space) &&
1060 dataOnElement[H1]->bAse.test(b) &&
1061 dataOnElement[H1]->basesOnSpaces[space].test(b)) {
1062 CHKERR getUserPolynomialBase()->getValue(
1063 gaussPts,
1064 boost::make_shared<EntPolynomialBaseCtx>(
1065 *dataOnElement[space], static_cast<FieldSpace>(space),
1066 static_cast<FieldApproximationBase>(b), NOBASE));
1067 }
1068 break;
1069 default:
1071 "Base <%s> not yet implemented",
1073 }
1074 }
1076}
1077
1080 /// Use the some node base. Node base is usually used for construction other
1081 /// bases.
1082 for (int space = HCURL; space != LASTSPACE; ++space) {
1083 dataOnElement[space]->dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE) =
1084 dataOnElement[H1]->dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE);
1085 dataOnElement[space]->dataOnEntities[MBVERTEX][0].getDiffNSharedPtr(
1086 NOBASE) =
1087 dataOnElement[H1]->dataOnEntities[MBVERTEX][0].getDiffNSharedPtr(
1088 NOBASE);
1089 }
1090 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
1092 static_cast<FieldApproximationBase>(b));
1093 }
1095}
1096
1100
1101 const auto ele_type = numeredEntFiniteElementPtr->getEntType();
1102
1103 auto get_nodal_base_data = [&](EntitiesFieldData &data, auto field_ptr) {
1105 auto &space = data.dataOnEntities[MBVERTEX][0].getSpace();
1106 auto &base = data.dataOnEntities[MBVERTEX][0].getBase();
1107 auto &bb_node_order = data.dataOnEntities[MBVERTEX][0].getBBNodeOrder();
1108
1109 auto &field_ents = getDataFieldEnts();
1110 auto bit_number = field_ptr->getBitNumber();
1111 const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
1112 bit_number, get_id_for_min_type<MBVERTEX>());
1113 auto lo = std::lower_bound(field_ents.begin(), field_ents.end(), lo_uid,
1114 cmp_uid_lo);
1115 if (lo != field_ents.end()) {
1116 const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
1117 bit_number, get_id_for_max_type<MBVERTEX>());
1118 auto hi = std::upper_bound(lo, field_ents.end(), hi_uid, cmp_uid_hi);
1119 if (lo != hi) {
1120
1121 for (auto it = lo; it != hi; ++it)
1122 if (auto first_e = it->lock()) {
1123 space = first_e->getSpace();
1124 base = first_e->getApproxBase();
1125 const int num_nodes = getNumberOfNodes();
1126 bb_node_order.resize(num_nodes, false);
1127 bb_node_order.clear();
1128 const int nb_dof_idx = first_e->getNbOfCoeffs();
1129
1130 std::vector<boost::weak_ptr<FieldEntity>> brother_ents_vec;
1131
1132 for (; it != hi; ++it) {
1133 if (auto e = it->lock()) {
1134 const auto &sn = e->getSideNumberPtr();
1135 const int side_number = sn->side_number;
1136 const int brother_side_number = sn->brother_side_number;
1137 if (brother_side_number != -1)
1138 brother_ents_vec.emplace_back(e);
1139 bb_node_order[side_number] = e->getMaxOrder();
1140 }
1141 }
1142
1143 for (auto &it : brother_ents_vec) {
1144 if (const auto e = it.lock()) {
1145 const auto &sn = e->getSideNumberPtr();
1146 const int side_number = sn->side_number;
1147 const int brother_side_number = sn->brother_side_number;
1148 bb_node_order[brother_side_number] = bb_node_order[side_number];
1149 }
1150 }
1151
1152 break;
1153 }
1154 }
1155 }
1156
1158 };
1159
1160 auto get_entity_base_data = [&](EntitiesFieldData &data, auto field_ptr,
1161 const EntityType type_lo,
1162 const EntityType type_hi) {
1164 for (EntityType t = MBEDGE; t != MBPOLYHEDRON; ++t) {
1165 for (auto &dat : data.dataOnEntities[t]) {
1166 dat.getOrder() = 0;
1167 dat.getBase() = NOBASE;
1168 dat.getSpace() = NOSPACE;
1169 dat.getFieldData().resize(0, false);
1170 dat.getFieldDofs().resize(0, false);
1171 }
1172 }
1173
1174 auto &field_ents = getDataFieldEnts();
1175 auto bit_number = field_ptr->getBitNumber();
1176 const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
1177 bit_number, get_id_for_min_type(type_lo));
1178 auto lo = std::lower_bound(field_ents.begin(), field_ents.end(), lo_uid,
1179 cmp_uid_lo);
1180 if (lo != field_ents.end()) {
1181 const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
1182 bit_number, get_id_for_max_type(type_hi));
1183 auto hi = std::upper_bound(lo, field_ents.end(), hi_uid, cmp_uid_hi);
1184
1185 std::vector<boost::weak_ptr<FieldEntity>> brother_ents_vec;
1186 for (; lo != hi; ++lo) {
1187 if (auto e = lo->lock()) {
1188 if (auto cache = e->entityCacheDataDofs.lock()) {
1189 if (cache->loHi[0] != cache->loHi[1]) {
1190 if (const auto side = e->getSideNumberPtr()->side_number;
1191 side >= 0) {
1192 const EntityType type = e->getEntType();
1193 auto &dat = data.dataOnEntities[type][side];
1194 const int brother_side =
1195 e->getSideNumberPtr()->brother_side_number;
1196 if (brother_side != -1)
1197 brother_ents_vec.emplace_back(e);
1198 dat.getBase() = e->getApproxBase();
1199 dat.getSpace() = e->getSpace();
1200 const auto ent_order = e->getMaxOrder();
1201 dat.getOrder() =
1202 dat.getOrder() > ent_order ? dat.getOrder() : ent_order;
1203 }
1204 }
1205 }
1206 }
1207 }
1208
1209 for (auto &ent_ptr : brother_ents_vec) {
1210 if (auto e = ent_ptr.lock()) {
1211 const EntityType type = e->getEntType();
1212 const int side = e->getSideNumberPtr()->side_number;
1213 const int brother_side = e->getSideNumberPtr()->brother_side_number;
1214 auto &dat = data.dataOnEntities[type][side];
1215 auto &dat_brother = data.dataOnEntities[type][brother_side];
1216 dat_brother.getBase() = dat.getBase();
1217 dat_brother.getSpace() = dat.getSpace();
1218 dat_brother.getOrder() = dat.getOrder();
1219 }
1220 }
1221 }
1223 };
1224
1225 for (auto &e : getDataFieldEnts()) {
1226 if (auto ent_data_ptr = e.lock()) {
1227 if (ent_data_ptr->getApproxBase() == AINSWORTH_BERNSTEIN_BEZIER_BASE) {
1228 auto space = ent_data_ptr->getSpace();
1229 for (EntityType t = MBVERTEX; t != MBPOLYHEDRON; ++t) {
1230 for (auto &dat : (*dataOnElement[space]).dataOnEntities[t]) {
1231 for (auto &ptr : dat.getBBAlphaIndicesByOrderArray())
1232 ptr.reset();
1233 for (auto &ptr : dat.getBBNByOrderArray())
1234 ptr.reset();
1235 for (auto &ptr : dat.getBBDiffNByOrderArray())
1236 ptr.reset();
1237 }
1238 }
1239 }
1240 }
1241 }
1242
1243 auto check_space = [&](const auto space) {
1244 switch (space) {
1245 case H1:
1246 for (auto t = MBVERTEX; t <= ele_type; ++t) {
1247 if (dataOnElement[H1]->spacesOnEntities[t].test(H1))
1248 return true;
1249 }
1250 return false;
1251 case HCURL:
1252 for (auto t = MBEDGE; t <= ele_type; ++t) {
1253 if (dataOnElement[HCURL]->spacesOnEntities[t].test(HCURL))
1254 return true;
1255 }
1256 return false;
1257 case HDIV:
1258 for (auto t = MBTRI; t <= ele_type; ++t) {
1259 if (dataOnElement[HDIV]->spacesOnEntities[t].test(HDIV))
1260 return true;
1261 }
1262 return false;
1263 case L2:
1264 return dataOnElement[L2]->spacesOnEntities[ele_type].test(L2);
1265 break;
1266 default:
1267 THROW_MESSAGE("Not implemented");
1268 }
1269 };
1270
1271 std::set<string> fields_list;
1272 for (auto &e : getDataFieldEnts()) {
1273 if (auto ent_data_ptr = e.lock()) {
1274 if (ent_data_ptr->getApproxBase() == AINSWORTH_BERNSTEIN_BEZIER_BASE) {
1275 auto field_name = ent_data_ptr->getName();
1276 if (fields_list.find(field_name) == fields_list.end()) {
1277 auto field_ptr = ent_data_ptr->getFieldRawPtr();
1278 auto space = ent_data_ptr->getSpace();
1279 CHKERR get_nodal_base_data(*dataOnElement[space], field_ptr);
1280 CHKERR get_entity_base_data(*dataOnElement[space], field_ptr, MBEDGE,
1281 ele_type);
1282 if (check_space(space)) {
1283 CHKERR getElementPolynomialBase()->getValue(
1284 gaussPts, boost::make_shared<EntPolynomialBaseCtx>(
1285 *dataOnElement[space], field_name,
1286 static_cast<FieldSpace>(space),
1288 fields_list.insert(field_name);
1289 }
1290 }
1291 }
1292 }
1293 }
1294
1296};
1297
1300
1301 // Data on elements for proper spaces
1302 for (int space = H1; space != LASTSPACE; ++space) {
1303 dataOnElement[space]->setElementType(type);
1304 derivedDataOnElement[space]->setElementType(type);
1305 }
1306
1308}
1309
1310#define FUNCTION_NAME_WITH_OP_NAME(OP) \
1311 std::ostringstream ss; \
1312 ss << "(Calling user data operator " \
1313 << boost::typeindex::type_id_runtime(OP).pretty_name() << " rowField " \
1314 << (OP).rowFieldName << " colField " << (OP).colFieldName << ") "
1315
1316#define CATCH_OP_ERRORS(OP) \
1317 catch (MoFEMExceptionInitial const &ex) { \
1318 FUNCTION_NAME_WITH_OP_NAME(OP) << PETSC_FUNCTION_NAME; \
1319 return PetscError(PETSC_COMM_SELF, ex.lINE, ss.str().c_str(), __FILE__, \
1320 ex.errorCode, PETSC_ERROR_INITIAL, ex.what()); \
1321 } \
1322 catch (MoFEMExceptionRepeat const &ex) { \
1323 FUNCTION_NAME_WITH_OP_NAME(OP) << PETSC_FUNCTION_NAME; \
1324 return PetscError(PETSC_COMM_SELF, ex.lINE, ss.str().c_str(), __FILE__, \
1325 ex.errorCode, PETSC_ERROR_REPEAT, " "); \
1326 } \
1327 catch (MoFEMException const &ex) { \
1328 FUNCTION_NAME_WITH_OP_NAME(OP) << ex.errorMessage; \
1329 SETERRQ(PETSC_COMM_SELF, ex.errorCode, ss.str().c_str()); \
1330 } \
1331 catch (std::exception const &ex) { \
1332 std::string message("Error: " + std::string(ex.what()) + " at " + \
1333 boost::lexical_cast<std::string>(__LINE__) + " : " + \
1334 std::string(__FILE__) + " in " + \
1335 std::string(PETSC_FUNCTION_NAME)); \
1336 FUNCTION_NAME_WITH_OP_NAME(OP) << message; \
1337 SETERRQ(PETSC_COMM_SELF, MOFEM_STD_EXCEPTION_THROW, ss.str().c_str()); \
1338 }
1339
1342
1343 using UDO = UserDataOperator;
1344
1345 std::array<std::string, 2> field_name;
1346 std::array<const Field *, 3> field_struture;
1347 std::array<int, 2>
1348 field_id; // note the this is field bit number (nor field bit)
1349 std::array<FieldSpace, 2> space;
1350 std::array<FieldApproximationBase, 2> base;
1351
1352 constexpr std::array<UDO::OpType, 2> types = {UDO::OPROW, UDO::OPCOL};
1353 std::array<int, 2> last_eval_field_id = {0, 0};
1354
1355 std::array<boost::shared_ptr<EntitiesFieldData>, 2> op_data;
1356
1357 auto swap_bases = [&](auto &op) {
1359 for (size_t ss = 0; ss != 2; ++ss) {
1360 if (op.getOpType() & types[ss] || op.getOpType() & UDO::OPROWCOL) {
1361 switch (base[ss]) {
1363 CHKERR op_data[ss]->baseSwap(field_name[ss],
1365 default:
1366 break;
1367 }
1368 }
1369 }
1370
1372 };
1373
1374 const EntityType type = numeredEntFiniteElementPtr->getEntType();
1375
1376 // evaluate entity data only, no field specific data provided or known
1377 auto evaluate_op_space = [&](auto &op) {
1379
1380 // reseat all data which all field dependent
1381 dataOnElement[op.sPace]->resetFieldDependentData();
1382 std::fill(last_eval_field_id.begin(), last_eval_field_id.end(), 0);
1383
1384 switch (op.sPace) {
1385 case NOSPACE:
1386 try {
1387 CHKERR op.doWork(
1388 0, MBENTITYSET,
1389 dataOnElement[NOSPACE]->dataOnEntities[MBENTITYSET][0]);
1390 }
1391 CATCH_OP_ERRORS(op);
1392 break;
1393 case NOFIELD:
1394 case H1:
1395 case HCURL:
1396 case HDIV:
1397 case L2:
1398 try {
1399
1400 for (EntityType t = MBVERTEX; t != MBMAXTYPE; ++t) {
1401 for (auto &e : dataOnElement[op.sPace]->dataOnEntities[t]) {
1402 e.getSpace() = op.sPace;
1403 }
1404 }
1405
1406 CHKERR op.opRhs(*dataOnElement[op.sPace], false);
1407 }
1408 CATCH_OP_ERRORS(op);
1409 break;
1410 default:
1411 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1412 "Not implemented for this space", op.sPace);
1413 }
1414
1416 };
1417
1418 // set entity data
1419 auto set_op_entityties_data = [&](auto ss, auto &op) {
1421
1422#ifndef NDEBUG
1423 if ((op.getNumeredEntFiniteElementPtr()->getBitFieldIdData() &
1425 .none()) {
1426 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1427 "no data field < %s > on finite element < %s >",
1428 field_name[ss].c_str(), getFEName().c_str());
1429 }
1430#endif
1431
1432 op_data[ss] =
1433 !ss ? dataOnElement[space[ss]] : derivedDataOnElement[space[ss]];
1434
1435 for (auto &data : op_data[ss]->dataOnEntities[MBENTITYSET]) {
1436 CHKERR data.resetFieldDependentData();
1437 }
1438
1439 auto get_data_for_nodes = [&]() {
1441 if (!ss)
1442 CHKERR getRowNodesIndices(*op_data[ss], field_id[ss]);
1443 else
1444 CHKERR getColNodesIndices(*op_data[ss], field_id[ss]);
1445 CHKERR getNodesFieldData(*op_data[ss], field_id[ss]);
1447 };
1448
1449 // get data on entities but not nodes
1450 auto get_data_for_entities = [&]() {
1452 CHKERR getEntityFieldData(*op_data[ss], field_id[ss], MBEDGE);
1453 if (!ss)
1454 CHKERR getEntityRowIndices(*op_data[ss], field_id[ss], MBEDGE);
1455 else
1456 CHKERR getEntityColIndices(*op_data[ss], field_id[ss], MBEDGE);
1458 };
1459
1460 auto get_data_for_meshset = [&]() {
1462 if (!ss) {
1463 CHKERR getNoFieldRowIndices(*op_data[ss], field_id[ss]);
1464 } else {
1465 CHKERR getNoFieldColIndices(*op_data[ss], field_id[ss]);
1466 }
1467 CHKERR getNoFieldFieldData(*op_data[ss], field_id[ss]);
1469 };
1470
1471 switch (space[ss]) {
1472 case H1:
1473 CHKERR get_data_for_nodes();
1474 case HCURL:
1475 case HDIV:
1476 CHKERR get_data_for_entities();
1477 break;
1478 case L2:
1479 switch (type) {
1480 case MBVERTEX:
1481 CHKERR get_data_for_nodes();
1482 break;
1483 default:
1484 CHKERR get_data_for_entities();
1485 }
1486 break;
1487 case NOFIELD:
1488 // if (!getNinTheLoop()) {
1489 // NOFIELD data are the same for each element, can be
1490 // retrieved only once
1491 CHKERR get_data_for_meshset();
1492 // }
1493 break;
1494 default:
1495 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1496 "not implemented for this space < %s >",
1497 FieldSpaceNames[space[ss]]);
1498 }
1500 };
1501
1502 // evalate operators with field data
1503 auto evaluate_op_for_fields = [&](auto &op) {
1505
1506 if (op.getOpType() & UDO::OPROW) {
1507 try {
1508 CHKERR op.opRhs(*op_data[0], false);
1509 }
1510 CATCH_OP_ERRORS(op);
1511 }
1512
1513 if (op.getOpType() & UDO::OPCOL) {
1514 try {
1515 CHKERR op.opRhs(*op_data[1], false);
1516 }
1517 CATCH_OP_ERRORS(op);
1518 }
1519
1520 if (op.getOpType() & UDO::OPROWCOL) {
1521 try {
1522 CHKERR op.opLhs(*op_data[0], *op_data[1]);
1523 }
1524 CATCH_OP_ERRORS(op);
1525 }
1527 };
1528
1529 // Collect bit ref level on all entities, fields and spaces
1531
1532 auto oit = opPtrVector.begin();
1533 auto hi_oit = opPtrVector.end();
1534
1535 // interate over all operators on element
1536 for (; oit != hi_oit; oit++) {
1537
1538 try {
1539
1540 CHKERR oit->setPtrFE(this);
1541
1542 if ((oit->opType & UDO::OPSPACE) == UDO::OPSPACE) {
1543
1544 // run operator for space but specific field
1545 CHKERR evaluate_op_space(*oit);
1546
1547 } else if (
1548
1549 (oit->opType & (UDO::OPROW | UDO::OPCOL | UDO::OPROWCOL)) ==
1550 oit->opType
1551
1552 ) {
1553
1554 field_name[0] = oit->rowFieldName;
1555 field_name[1] = oit->colFieldName;
1556
1557 for (size_t ss = 0; ss != 2; ss++) {
1558
1559#ifndef NDEBUG
1560 if (field_name[ss].empty()) {
1561 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1562 "Not set Field name in operator %d (0-row, 1-column) in "
1563 "operator %s",
1564 ss,
1565 (boost::typeindex::type_id_runtime(*oit).pretty_name())
1566 .c_str());
1567 }
1568#endif
1569
1570 field_struture[ss] = mField.get_field_structure(field_name[ss]);
1571 field_id[ss] = field_struture[ss]->getBitNumber();
1572 space[ss] = field_struture[ss]->getSpace();
1573 base[ss] = field_struture[ss]->getApproxBase();
1574 }
1575
1576 // not that if field name do not change between operators, entity field
1577 // data are nor rebuild
1578 for (size_t ss = 0; ss != 2; ss++) {
1579
1580 if (oit->getOpType() & types[ss] ||
1581 oit->getOpType() & UDO::OPROWCOL) {
1582 if (last_eval_field_id[ss] != field_id[ss]) {
1583 CHKERR set_op_entityties_data(ss, *oit);
1584 last_eval_field_id[ss] = field_id[ss];
1585 }
1586 }
1587 }
1588
1589 CHKERR swap_bases(*oit);
1590
1591 // run operator for given field or row, column or both
1592 CHKERR evaluate_op_for_fields(*oit);
1593
1594 CHKERR swap_bases(*oit);
1595
1596 } else {
1597
1598 for (int i = 0; i != 5; ++i)
1599 if (oit->opType & (1 << i))
1600 MOFEM_LOG("SELF", Sev::error) << UDO::OpTypeNames[i];
1601 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1602 "Impossible operator type");
1603 }
1604 }
1605 CATCH_OP_ERRORS(*oit);
1606 }
1608}
1609
1610const char *const ForcesAndSourcesCore::UserDataOperator::OpTypeNames[] = {
1611 "OPROW", " OPCOL", "OPROWCOL", "OPSPACE", "OPLAST"};
1612
1613MoFEMErrorCode ForcesAndSourcesCore::UserDataOperator::getProblemRowIndices(
1614 const std::string field_name, const EntityType type, const int side,
1615 VectorInt &indices) const {
1617 if (ptrFE == NULL)
1618 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1619
1620 switch (type) {
1621 case MBVERTEX:
1623 break;
1624 default:
1625 CHKERR ptrFE->getProblemTypeRowIndices(field_name, type, side, indices);
1626 }
1628}
1629
1630MoFEMErrorCode ForcesAndSourcesCore::UserDataOperator::getProblemColIndices(
1631 const std::string field_name, const EntityType type, const int side,
1632 VectorInt &indices) const {
1634 if (ptrFE == NULL)
1635 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1636
1637 switch (type) {
1638 case MBVERTEX:
1639 CHKERR ptrFE->getProblemNodesColIndices(field_name, indices);
1640 break;
1641 default:
1642 CHKERR ptrFE->getProblemTypeColIndices(field_name, type, side, indices);
1643 }
1645}
1646
1650 sidePtrFE = const_cast<ForcesAndSourcesCore *>(side_fe_ptr);
1652}
1653
1655 const ForcesAndSourcesCore *refine_fe_ptr) {
1657 refinePtrFE = const_cast<ForcesAndSourcesCore *>(refine_fe_ptr);
1659}
1660
1661MoFEMErrorCode ForcesAndSourcesCore::UserDataOperator::loopSide(
1662 const string &fe_name, ForcesAndSourcesCore *side_fe, const size_t side_dim,
1663 const EntityHandle ent_for_side, const int verb,
1664 const LogManager::SeverityLevel sev) {
1666 const auto ent = ent_for_side ? ent_for_side : getFEEntityHandle();
1667 const auto *problem_ptr = getFEMethod()->problemPtr;
1668
1669 side_fe->feName = fe_name;
1670
1671 CHKERR side_fe->setSideFEPtr(ptrFE);
1672 CHKERR side_fe->copyBasicMethod(*getFEMethod());
1673 CHKERR side_fe->copyPetscData(*getFEMethod());
1674 CHKERR side_fe->copyKsp(*getFEMethod());
1675 CHKERR side_fe->copySnes(*getFEMethod());
1676 CHKERR side_fe->copyTs(*getFEMethod());
1677
1678 CHKERR side_fe->preProcess();
1679
1680 Range adjacent_ents;
1681 CHKERR ptrFE->mField.getInterface<BitRefManager>()->getAdjacenciesAny(
1682 ent, side_dim, adjacent_ents);
1683 auto &numered_fe =
1684 problem_ptr->numeredFiniteElementsPtr->get<Unique_mi_tag>();
1685
1686 auto fe_miit = ptrFE->mField.get_finite_elements()
1688 .find(fe_name);
1689 if (fe_miit != ptrFE->mField.get_finite_elements()
1691 .end()) {
1692 int nn = 0;
1693 side_fe->loopSize = adjacent_ents.size();
1694 for (auto fe_ent : adjacent_ents) {
1695 auto miit = numered_fe.find(EntFiniteElement::getLocalUniqueIdCalculate(
1696 fe_ent, (*fe_miit)->getFEUId()));
1697 if (miit != numered_fe.end()) {
1698 if (verb >= VERBOSE)
1699 MOFEM_LOG("SELF", sev) << "Side finite element "
1700 << "(" << nn << "): " << **miit;
1701 side_fe->nInTheLoop = nn++;
1702 side_fe->numeredEntFiniteElementPtr = *miit;
1703 CHKERR (*side_fe)();
1704 }
1705 }
1706 }
1707
1708 CHKERR side_fe->postProcess();
1709
1711}
1712
1713MoFEMErrorCode ForcesAndSourcesCore::UserDataOperator::loopThis(
1714 const string &fe_name, ForcesAndSourcesCore *this_fe, const int verb,
1715 const LogManager::SeverityLevel sev) {
1717
1718 if (verb >= VERBOSE)
1719 MOFEM_LOG("SELF", sev) << "This finite element: "
1720 << *getNumeredEntFiniteElementPtr();
1721
1722 const auto *problem_ptr = getFEMethod()->problemPtr;
1723 this_fe->feName = fe_name;
1724
1725 CHKERR this_fe->setRefineFEPtr(ptrFE);
1726 CHKERR this_fe->copyBasicMethod(*getFEMethod());
1727 CHKERR this_fe->copyPetscData(*getFEMethod());
1728 CHKERR this_fe->copyKsp(*getFEMethod());
1729 CHKERR this_fe->copySnes(*getFEMethod());
1730 CHKERR this_fe->copyTs(*getFEMethod());
1731
1732 CHKERR this_fe->preProcess();
1733
1734 this_fe->nInTheLoop = getNinTheLoop();
1735 this_fe->numeredEntFiniteElementPtr = getNumeredEntFiniteElementPtr();
1736
1737 CHKERR (*this_fe)();
1738
1739 CHKERR this_fe->postProcess();
1740
1742}
1743
1744MoFEMErrorCode ForcesAndSourcesCore::UserDataOperator::loopParent(
1745 const string &fe_name, ForcesAndSourcesCore *parent_fe, const int verb,
1746 const LogManager::SeverityLevel sev) {
1748
1749 const auto *problem_ptr = getFEMethod()->problemPtr;
1750 auto &numered_fe =
1751 problem_ptr->numeredFiniteElementsPtr->get<Unique_mi_tag>();
1752
1753 parent_fe->feName = fe_name;
1754 CHKERR parent_fe->setRefineFEPtr(ptrFE);
1755 CHKERR parent_fe->copyBasicMethod(*getFEMethod());
1756 CHKERR parent_fe->copyPetscData(*getFEMethod());
1757 CHKERR parent_fe->copyKsp(*getFEMethod());
1758 CHKERR parent_fe->copySnes(*getFEMethod());
1759 CHKERR parent_fe->copyTs(*getFEMethod());
1760
1761 CHKERR parent_fe->preProcess();
1762
1763 auto fe_miit = ptrFE->mField.get_finite_elements()
1765 .find(fe_name);
1766 if (fe_miit != ptrFE->mField.get_finite_elements()
1768 .end()) {
1769 const auto parent_ent = getNumeredEntFiniteElementPtr()->getParentEnt();
1770 auto miit = numered_fe.find(EntFiniteElement::getLocalUniqueIdCalculate(
1771 parent_ent, (*fe_miit)->getFEUId()));
1772 if (miit != numered_fe.end()) {
1773 if (verb >= VERBOSE)
1774 MOFEM_LOG("SELF", sev) << "Parent finite element: " << **miit;
1775 parent_fe->loopSize = 1;
1776 parent_fe->nInTheLoop = 0;
1777 parent_fe->numeredEntFiniteElementPtr = *miit;
1778 CHKERR (*parent_fe)();
1779 }
1780 }
1781
1782 CHKERR parent_fe->postProcess();
1783
1785}
1786
1787MoFEMErrorCode ForcesAndSourcesCore::UserDataOperator::loopChildren(
1788 const string &fe_name, ForcesAndSourcesCore *child_fe, const int verb,
1789 const LogManager::SeverityLevel sev) {
1791
1792 const auto *problem_ptr = getFEMethod()->problemPtr;
1793 auto &ref_ents = *getPtrFE()->mField.get_ref_ents();
1794 auto &numered_fe =
1795 problem_ptr->numeredFiniteElementsPtr->get<Unique_mi_tag>();
1796
1797 const auto parent_ent = getNumeredEntFiniteElementPtr()->getEnt();
1798 const auto parent_type = getNumeredEntFiniteElementPtr()->getEntType();
1799 auto range =
1800 ref_ents.get<Composite_ParentEnt_And_EntType_mi_tag>().equal_range(
1801 boost::make_tuple(parent_type, parent_ent));
1802
1803 if (auto size = std::distance(range.first, range.second)) {
1804
1805 std::vector<EntityHandle> childs_vec;
1806 childs_vec.reserve(size);
1807 for (; range.first != range.second; ++range.first)
1808 childs_vec.emplace_back((*range.first)->getEnt());
1809
1810 Range childs;
1811
1812 if ((childs_vec.back() - childs_vec.front() + 1) == size)
1813 childs = Range(childs_vec.front(), childs_vec.back());
1814 else
1815 childs.insert_list(childs_vec.begin(), childs_vec.end());
1816
1817 child_fe->feName = fe_name;
1818
1819 CHKERR child_fe->setRefineFEPtr(ptrFE);
1820 CHKERR child_fe->copyBasicMethod(*getFEMethod());
1821 CHKERR child_fe->copyPetscData(*getFEMethod());
1822 CHKERR child_fe->copyKsp(*getFEMethod());
1823 CHKERR child_fe->copySnes(*getFEMethod());
1824 CHKERR child_fe->copyTs(*getFEMethod());
1825
1826 CHKERR child_fe->preProcess();
1827
1828 int nn = 0;
1829 child_fe->loopSize = size;
1830
1831 auto fe_miit = ptrFE->mField.get_finite_elements()
1833 .find(fe_name);
1834 if (fe_miit != ptrFE->mField.get_finite_elements()
1836 .end()) {
1837
1838 for (auto p = childs.pair_begin(); p != childs.pair_end(); ++p) {
1839
1840 auto miit =
1841 numered_fe.lower_bound(EntFiniteElement::getLocalUniqueIdCalculate(
1842 p->first, (*fe_miit)->getFEUId()));
1843 auto hi_miit =
1844 numered_fe.upper_bound(EntFiniteElement::getLocalUniqueIdCalculate(
1845 p->second, (*fe_miit)->getFEUId()));
1846
1847 for (; miit != hi_miit; ++miit) {
1848
1849 if (verb >= VERBOSE)
1850 MOFEM_LOG("SELF", sev) << "Child finite element "
1851 << "(" << nn << "): " << **miit;
1852
1853 child_fe->nInTheLoop = nn++;
1854 child_fe->numeredEntFiniteElementPtr = *miit;
1855 CHKERR (*child_fe)();
1856 }
1857 }
1858 }
1859
1860 CHKERR child_fe->postProcess();
1861 };
1862
1864}
1865
1866int ForcesAndSourcesCore::getRule(int order_row, int order_col,
1867 int order_data) {
1868 return getRuleHook ? getRuleHook(order_row, order_col, order_data)
1869 : getRule(order_data);
1870}
1871
1873 int order_data) {
1874 return setRuleHook ? setRuleHook(this, order_row, order_col, order_data)
1875 : setGaussPts(order_data);
1876}
1877
1879
1880/** \deprecated setGaussPts(int row_order, int col_order, int data order);
1881 */
1884 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Sorry, not implemented");
1886}
1887
1888ForcesAndSourcesCore::UserDataOperator::UserDataOperator(const FieldSpace space,
1889 const char type,
1890 const bool symm)
1891 : DataOperator(symm), opType(type), sPace(space), ptrFE(nullptr) {}
1892
1893ForcesAndSourcesCore::UserDataOperator::UserDataOperator(
1894 const std::string field_name, const char type, const bool symm)
1895 : DataOperator(symm), opType(type), rowFieldName(field_name),
1896 colFieldName(field_name), sPace(LASTSPACE), ptrFE(nullptr) {}
1897
1898ForcesAndSourcesCore::UserDataOperator::UserDataOperator(
1899 const std::string row_field_name, const std::string col_field_name,
1900 const char type, const bool symm)
1901 : DataOperator(symm), opType(type), rowFieldName(row_field_name),
1902 colFieldName(col_field_name), sPace(LASTSPACE), ptrFE(nullptr) {}
1903
1906 if (preProcessHook) {
1907 ierr = preProcessHook();
1908 CHKERRG(ierr);
1909 }
1911}
1914 if (operatorHook) {
1915 ierr = operatorHook();
1916 CHKERRG(ierr);
1917 } else {
1918#ifndef NDEBUG
1919 MOFEM_LOG("SELF", Sev::warning)
1920 << "No method operator() overloaded on element entity on finite "
1921 "element <"
1922 << boost::typeindex::type_id_runtime(*this).pretty_name() << ">";
1923#endif
1924 }
1926}
1929 if (postProcessHook) {
1931 CHKERRG(ierr);
1932 }
1934}
1935
1937ForcesAndSourcesCore::UserDataOperator::setPtrFE(ForcesAndSourcesCore *ptr) {
1939 if (!(ptrFE = dynamic_cast<ForcesAndSourcesCore *>(ptr)))
1940 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1941 "User operator and finite element do not work together");
1943}
1944
1945} // namespace MoFEM
static Index< 'p', 3 > p
#define CATCH_OP_ERRORS(OP)
constexpr double a
@ VERBOSE
Definition: definitions.h:209
FieldApproximationBase
approximation base
Definition: definitions.h:58
@ LASTBASE
Definition: definitions.h:69
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:62
@ USER_BASE
user implemented approximation base
Definition: definitions.h:68
@ 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()
Definition: definitions.h:447
FieldSpace
approximation spaces
Definition: definitions.h:82
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ LASTSPACE
FieldSpace in [ 0, LASTSPACE )
Definition: definitions.h:89
@ NOFIELD
scalar or vector of scalars describe (no true field)
Definition: definitions.h:84
@ H1
continuous field
Definition: definitions.h:85
@ NOSPACE
Definition: definitions.h:83
@ HCURL
field with continuous tangents
Definition: definitions.h:86
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
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 ...
Definition: definitions.h:346
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
@ 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()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:561
FTensor::Index< 'n', SPACE_DIM > n
multi_index_container< boost::shared_ptr< NumeredDofEntity >, indexed_by< ordered_unique< tag< Unique_mi_tag >, const_mem_fun< NumeredDofEntity::interface_type_DofEntity, UId, &NumeredDofEntity::getLocalUniqueId > >, ordered_non_unique< tag< Part_mi_tag >, member< NumeredDofEntity, unsigned int, &NumeredDofEntity::pArt > >, ordered_non_unique< tag< Idx_mi_tag >, member< NumeredDofEntity, DofIdx, &NumeredDofEntity::dofIdx > >, ordered_non_unique< tag< PetscGlobalIdx_mi_tag >, member< NumeredDofEntity, DofIdx, &NumeredDofEntity::petscGloablDofIdx > >, ordered_non_unique< tag< PetscLocalIdx_mi_tag >, member< NumeredDofEntity, DofIdx, &NumeredDofEntity::petscLocalDofIdx > >, ordered_non_unique< tag< Ent_mi_tag >, const_mem_fun< NumeredDofEntity::interface_type_DofEntity, EntityHandle, &NumeredDofEntity::getEnt > > > > NumeredDofEntity_multiIndex
MultiIndex container keeps NumeredDofEntity.
multi_index_container< boost::shared_ptr< SideNumber >, indexed_by< ordered_unique< member< SideNumber, EntityHandle, &SideNumber::ent > >, ordered_non_unique< composite_key< SideNumber, const_mem_fun< SideNumber, EntityType, &SideNumber::getEntType >, member< SideNumber, signed char, &SideNumber::side_number > > > > > SideNumber_multiIndex
SideNumber_multiIndex for SideNumber.
virtual const Field * get_field_structure(const std::string &name, enum MoFEMTypes bh=MF_EXIST) const =0
get field structure
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
SeverityLevel
Severity levels.
Definition: LogManager.hpp:33
FTensor::Index< 'i', SPACE_DIM > i
BlockParamData * cache
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
std::bitset< BITFEID_SIZE > BitFEId
Finite element Id.
Definition: Types.hpp:43
int ApproximationOrder
Approximation on the entity.
Definition: Types.hpp:26
std::bitset< BITFIELDID_SIZE > BitFieldId
Field Id.
Definition: Types.hpp:42
uint128_t UId
Unique Id.
Definition: Types.hpp:31
UBlasVector< int > VectorInt
Definition: Types.hpp:67
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
static auto cmp_uid_lo(const boost::weak_ptr< FieldEntity > &a, const UId &b)
EntityHandle get_id_for_max_type()
static auto cmp_uid_hi(const UId &b, const boost::weak_ptr< FieldEntity > &a)
EntityHandle get_id_for_min_type()
ublas::vector< FieldEntity *, FieldEntAllocator > VectorFieldEntities
std::vector< boost::weak_ptr< FieldEntity > > FieldEntity_vector_view
Field_multiIndex::index< FieldName_mi_tag >::type::iterator field_it
static int getMaxOrder(const ENTMULTIINDEX &multi_index)
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
constexpr double t
plate stiffness
Definition: plate.cpp:59
constexpr auto field_name
int loopSize
local number oe methods to process
MoFEMErrorCode copyBasicMethod(const BasicMethod &basic)
Copy data from other base method to this base method.
boost::function< MoFEMErrorCode()> preProcessHook
Hook function for pre-processing.
const Field_multiIndex * fieldsPtr
raw pointer to fields container
int nInTheLoop
number currently of processed method
const Problem * problemPtr
raw pointer to problem
boost::function< MoFEMErrorCode()> operatorHook
Hook function for operator.
int getNinTheLoop() const
get number of evaluated element in the loop
boost::function< MoFEMErrorCode()> postProcessHook
Hook function for post-processing.
Managing BitRefLevels.
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
virtual moab::Interface & get_moab()=0
virtual BitFieldId get_field_id(const std::string &name) const =0
Get field Id.
virtual MPI_Comm & get_comm() const =0
base operator to do operations at Gauss Pt. level
static bool getDefTypeMap(const EntityType fe_type, const EntityType ent_type)
Deprecated interface functions.
this class derive data form other data structure
UId getLocalUniqueIdCalculate() const
Generate UId for finite element entity.
Data on single entity (This is passed as argument to DataOperator::doWork)
data structure for finite element entity
std::bitset< LASTSPACE > sPace
spaces on element
MatrixInt facesNodes
nodes on finite element faces
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
std::bitset< LASTBASE > bAse
bases on element
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
std::array< std::bitset< LASTBASE >, LASTSPACE > basesOnSpaces
base on spaces
MatrixInt facesNodesOrder
order of face nodes on element
std::array< std::bitset< LASTBASE >, MBMAXTYPE > basesOnEntities
bases on entity types
keeps information about indexed dofs for the finite element
std::string feName
Name of finite element.
auto & getRowFieldEnts() const
auto getFEName() const
get finite element name
const FieldEntity_vector_view & getDataFieldEnts() const
boost::shared_ptr< FieldEntity_vector_view > & getDataFieldEntsPtr() const
auto getColDofsPtr() const
EntityHandle getFEEntityHandle() const
auto getNumberOfNodes() const
auto getRowDofsPtr() const
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
auto & getColFieldEnts() const
static UId getHiBitNumberUId(const FieldBitNumber bit_number)
UId getLocalUniqueIdCalculate()
Get the Local Unique Id Calculate object.
static UId getLoBitNumberUId(const FieldBitNumber bit_number)
static UId getLoLocalEntityBitNumber(const char bit_number, const EntityHandle ent)
static UId getHiLocalEntityBitNumber(const char bit_number, const EntityHandle ent)
Provide data structure for (tensor) field approximation.
FieldCoefficientsNumber getNbOfCoeffs() const
Get number of field coefficients.
FieldSpace getSpace() const
Get field approximation space.
FieldBitNumber getBitNumber() const
Get number of set bit in Field ID. Each field has uid, get getBitNumber get number of bit set for giv...
structure to get information form mofem into EntitiesFieldData
MoFEMErrorCode loopOverOperators()
Iterate user data operators.
MoFEMErrorCode getProblemNodesRowIndices(const std::string &field_name, VectorInt &nodes_indices) const
MoFEMErrorCode getEntityIndices(EntitiesFieldData &data, const int bit_number, FieldEntity_vector_view &ents_field, const EntityType type_lo, const EntityType type_hi, EXTRACTOR &&extractor) const
MoFEMErrorCode getEntityRowIndices(EntitiesFieldData &data, const int bit_number, const EntityType type_lo=MBVERTEX, const EntityType type_hi=MBPOLYHEDRON) const
MoFEMErrorCode getRowNodesIndices(EntitiesFieldData &data, const int bit_number) const
get row node indices from FENumeredDofEntity_multiIndex
int getMaxRowOrder() const
Get max order of approximation for field in rows.
auto & getUserPolynomialBase()
Get the User Polynomial Base object.
MoFEMErrorCode getSpacesAndBaseOnEntities(EntitiesFieldData &data) const
Get field approximation space and base on entities.
MoFEMErrorCode getProblemTypeColIndices(const std::string &field_name, EntityType type, int side_number, VectorInt &indices) const
MoFEMErrorCode getNoFieldIndices(const int bit_number, boost::shared_ptr< FENumeredDofEntity_multiIndex > dofs, VectorInt &nodes_indices) const
get NoField indices
virtual MoFEMErrorCode operator()()
function is run for every finite element
MoFEMErrorCode getFaceNodes(EntitiesFieldData &data) const
Get nodes on faces.
const std::array< boost::shared_ptr< EntitiesFieldData >, LASTSPACE > derivedDataOnElement
Entity data on element entity columns fields.
MoFEMErrorCode getNoFieldRowIndices(EntitiesFieldData &data, const int bit_number) const
get col NoField indices
virtual MoFEMErrorCode setGaussPts(int order_row, int order_col, int order_data)
set user specific integration rule
MoFEMErrorCode getNodesIndices(const int bit_number, FieldEntity_vector_view &ents_field, VectorInt &nodes_indices, VectorInt &local_nodes_indices, EXTRACTOR &&extractor) const
get node indices
MoFEMErrorCode getEntityFieldData(EntitiesFieldData &data, const int bit_number, const EntityType type_lo=MBVERTEX, const EntityType type_hi=MBPOLYHEDRON) const
MoFEMErrorCode setRefineFEPtr(const ForcesAndSourcesCore *refine_fe_ptr)
Set the pointer to face element refined.
auto & getElementPolynomialBase()
Get the Entity Polynomial Base object.
ForcesAndSourcesCore * sidePtrFE
Element to integrate on the sides.
ForcesAndSourcesCore(Interface &m_field)
virtual int getRule(int order_row, int order_col, int order_data)
another variant of getRule
boost::ptr_deque< UserDataOperator > opPtrVector
Vector of finite element users data operators.
MoFEMErrorCode getEntitySense(const EntityType type, boost::ptr_vector< EntitiesFieldData::EntData > &data) const
get sense (orientation) of entity
const std::array< boost::shared_ptr< EntitiesFieldData >, LASTSPACE > dataOnElement
Entity data on element entity rows fields.
MoFEMErrorCode getProblemNodesColIndices(const std::string &field_name, VectorInt &nodes_indices) const
MoFEMErrorCode calHierarchicalBaseFunctionsOnElement()
Calculate base functions.
MatrixDouble gaussPts
Matrix of integration points.
MoFEMErrorCode calBernsteinBezierBaseFunctionsOnElement()
Calculate Bernstein-Bezier base.
virtual MoFEMErrorCode preProcess()
function is run at the beginning of loop
MoFEMErrorCode getColNodesIndices(EntitiesFieldData &data, const int bit_number) const
get col node indices from FENumeredDofEntity_multiIndex
MoFEMErrorCode setSideFEPtr(const ForcesAndSourcesCore *side_fe_ptr)
Set the pointer to face element on the side.
int getMaxColOrder() const
Get max order of approximation for field in columns.
MoFEMErrorCode getProblemNodesIndices(const std::string &field_name, const NumeredDofEntity_multiIndex &dofs, VectorInt &nodes_indices) const
get indices of nodal indices which are declared for problem but not this particular element
GaussHookFun setRuleHook
Set function to calculate integration rule.
MoFEMErrorCode getNodesFieldData(EntitiesFieldData &data, const int bit_number) const
Get data on nodes.
MoFEMErrorCode getEntityColIndices(EntitiesFieldData &data, const int bit_number, const EntityType type_lo=MBVERTEX, const EntityType type_hi=MBPOLYHEDRON) const
MoFEMErrorCode getEntityDataOrder(const EntityType type, const FieldSpace space, boost::ptr_vector< EntitiesFieldData::EntData > &data) const
Get the entity data order.
ForcesAndSourcesCore * refinePtrFE
Element to integrate parent or child.
MoFEMErrorCode getNoFieldColIndices(EntitiesFieldData &data, const int bit_number) const
get col NoField indices
int getMaxDataOrder() const
Get max order of approximation for data fields.
MoFEMErrorCode getNoFieldFieldData(const int bit_number, VectorDouble &ent_field_data, VectorDofs &ent_field_dofs, VectorFieldEntities &ent_field) const
Get field data on nodes.
MoFEMErrorCode getProblemTypeRowIndices(const std::string &field_name, EntityType type, int side_number, VectorInt &indices) const
MoFEMErrorCode getProblemTypeIndices(const std::string &field_name, const NumeredDofEntity_multiIndex &dofs, EntityType type, int side_number, VectorInt &indices) const
get indices by type (generic function) which are declared for problem but not this particular element
RuleHookFun getRuleHook
Hook to get rule.
MoFEMErrorCode createDataOnElement(EntityType type)
Create a entity data on element object.
virtual MoFEMErrorCode postProcess()
function is run at the end of loop
MoFEMErrorCode copyKsp(const KspMethod &ksp)
copy data form another method
Definition: LoopMethods.cpp:55
MoFEMErrorCode copyPetscData(const PetscData &petsc_data)
Definition: LoopMethods.cpp:31
boost::shared_ptr< NumeredDofEntity_multiIndex > numeredRowDofsPtr
store DOFs on rows for this problem
boost::shared_ptr< NumeredDofEntity_multiIndex > numeredColDofsPtr
store DOFs on columns for this problem
MoFEMErrorCode copySnes(const SnesMethod &snes)
Copy snes data.
Definition: LoopMethods.cpp:75
MoFEMErrorCode copyTs(const TSMethod &ts)
Copy TS solver data.
Definition: LoopMethods.cpp:96
boost::shared_ptr< SideNumber > getSideNumberPtr() const