v0.14.0
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 (auto &data : dataOnElement) {
599 if (data) {
600 for (auto &dat : data->dataOnEntities) {
601 for (auto &ent_dat : dat) {
602 ent_dat.getEntDataBitRefLevel().clear();
603 }
604 }
605 }
606 }
607
608 auto &field_ents = getDataFieldEnts();
609 for (auto it : field_ents) {
610 if (auto e = it.lock()) {
611 const FieldSpace space = e->getSpace();
612 if (space > NOFIELD) {
613 const EntityType type = e->getEntType();
614 const signed char side = e->getSideNumberPtr()->side_number;
615 if (side >= 0) {
616 if (auto &data = dataOnElement[space]) {
617 if (type == MBVERTEX) {
618 auto &dat = data->dataOnEntities[type][0];
619 dat.getEntDataBitRefLevel().resize(getNumberOfNodes(), false);
620 dat.getEntDataBitRefLevel()[side] = e->getBitRefLevel();
621 } else {
622 auto &dat = data->dataOnEntities[type][side];
623 dat.getEntDataBitRefLevel().resize(1, false);
624 dat.getEntDataBitRefLevel()[0] = e->getBitRefLevel();
625 }
626 }
627 }
628 } else {
629 if (auto &data = dataOnElement[NOFIELD]) {
630 auto &dat = data->dataOnEntities[MBENTITYSET][0];
631 dat.getEntDataBitRefLevel().resize(1, false);
632 dat.getEntDataBitRefLevel()[0] = e->getBitRefLevel();
633 }
634 }
635 }
636 }
637
639};
640
643 const int bit_number) const {
644
645 auto get_nodes_field_data = [&](VectorDouble &nodes_data,
646 VectorFieldEntities &field_entities,
647 VectorDofs &nodes_dofs, FieldSpace &space,
649 VectorInt &bb_node_order) {
651
652 nodes_data.resize(0, false);
653 nodes_dofs.resize(0, false);
654 field_entities.resize(0, false);
655
656 auto field_it =
657 fieldsPtr->get<BitFieldId_mi_tag>().find(BitFEId().set(bit_number - 1));
658 if (field_it != fieldsPtr->get<BitFieldId_mi_tag>().end()) {
659
660#ifndef NDEBUG
661 if ((*field_it)->getBitNumber() != bit_number)
662 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong bit number");
663#endif
664 const int nb_dofs_on_vert = (*field_it)->getNbOfCoeffs();
665 space = (*field_it)->getSpace();
666 base = (*field_it)->getApproxBase();
667
668 auto &field_ents = getDataFieldEnts();
669 const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
670 bit_number, get_id_for_min_type<MBVERTEX>());
671 auto lo = std::lower_bound(field_ents.begin(), field_ents.end(), lo_uid,
672 cmp_uid_lo);
673 if (lo != field_ents.end()) {
674 const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
675 bit_number, get_id_for_max_type<MBVERTEX>());
676 auto hi = std::upper_bound(lo, field_ents.end(), hi_uid, cmp_uid_hi);
677 if (lo != hi) {
678
679 int nb_dofs = 0;
680 for (auto it = lo; it != hi; ++it) {
681 if (auto e = it->lock()) {
682 if (auto cache = e->entityCacheDataDofs.lock()) {
683 if (cache->loHi[0] != cache->loHi[1]) {
684 nb_dofs += std::distance(cache->loHi[0], cache->loHi[1]);
685 break;
686 }
687 }
688 }
689 }
690
691 if (nb_dofs) {
692
693 const int num_nodes = getNumberOfNodes();
694 bb_node_order.resize(num_nodes, false);
695 bb_node_order.clear();
696 const int max_nb_dofs = nb_dofs_on_vert * num_nodes;
697 nodes_data.resize(max_nb_dofs, false);
698 nodes_dofs.resize(max_nb_dofs, false);
699 field_entities.resize(num_nodes, false);
700 std::fill(nodes_data.begin(), nodes_data.end(), 0);
701 std::fill(nodes_dofs.begin(), nodes_dofs.end(), nullptr);
702 std::fill(field_entities.begin(), field_entities.end(), nullptr);
703
704 std::vector<boost::weak_ptr<FieldEntity>> brother_ents_vec;
705
706 for (auto it = lo; it != hi; ++it) {
707 if (auto e = it->lock()) {
708 const auto &sn = e->getSideNumberPtr();
709 // Some field entities on skeleton can have negative side
710 // number
711 if (const auto side_number = sn->side_number;
712 side_number >= 0) {
713 const int brother_side_number = sn->brother_side_number;
714
715 field_entities[side_number] = e.get();
716 if (brother_side_number != -1) {
717 brother_ents_vec.emplace_back(e);
718 field_entities[side_number] = field_entities[side_number];
719 }
720
721 bb_node_order[side_number] = e->getMaxOrder();
722 int pos = side_number * nb_dofs_on_vert;
723 auto ent_filed_data_vec = e->getEntFieldData();
724 if (auto cache = e->entityCacheDataDofs.lock()) {
725 for (auto dit = cache->loHi[0]; dit != cache->loHi[1];
726 ++dit) {
727 const auto dof_idx = (*dit)->getEntDofIdx();
728 nodes_data[pos + dof_idx] = ent_filed_data_vec[dof_idx];
729 nodes_dofs[pos + dof_idx] =
730 reinterpret_cast<FEDofEntity *>((*dit).get());
731 }
732 }
733 }
734 }
735 }
736
737 for (auto &it : brother_ents_vec) {
738 if (const auto e = it.lock()) {
739 const auto &sn = e->getSideNumberPtr();
740 const int side_number = sn->side_number;
741 const int brother_side_number = sn->brother_side_number;
742 bb_node_order[brother_side_number] = bb_node_order[side_number];
743 int pos = side_number * nb_dofs_on_vert;
744 int brother_pos = brother_side_number * nb_dofs_on_vert;
745 for (int ii = 0; ii != nb_dofs_on_vert; ++ii) {
746 nodes_data[brother_pos] = nodes_data[pos];
747 nodes_dofs[brother_pos] = nodes_dofs[pos];
748 ++pos;
749 ++brother_pos;
750 }
751 }
752 }
753 }
754 }
755 }
756 }
757
759 };
760
761 return get_nodes_field_data(
762 data.dataOnEntities[MBVERTEX][0].getFieldData(),
763 data.dataOnEntities[MBVERTEX][0].getFieldEntities(),
764 data.dataOnEntities[MBVERTEX][0].getFieldDofs(),
765 data.dataOnEntities[MBVERTEX][0].getSpace(),
766 data.dataOnEntities[MBVERTEX][0].getBase(),
767 data.dataOnEntities[MBVERTEX][0].getBBNodeOrder());
768}
769
771 EntitiesFieldData &data, const int bit_number, const EntityType type_lo,
772 const EntityType type_hi) const {
774 for (EntityType t = type_lo; t != type_hi; ++t) {
775 for (auto &dat : data.dataOnEntities[t]) {
776 dat.getOrder() = 0;
777 dat.getBase() = NOBASE;
778 dat.getSpace() = NOSPACE;
779 dat.getFieldData().resize(0, false);
780 dat.getFieldDofs().resize(0, false);
781 dat.getFieldEntities().resize(0, false);
782 }
783 }
784
785 auto &field_ents = getDataFieldEnts();
786 const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
787 bit_number, get_id_for_min_type(type_lo));
788 auto lo = std::lower_bound(field_ents.begin(), field_ents.end(), lo_uid,
789 cmp_uid_lo);
790 if (lo != field_ents.end()) {
791 const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
792 bit_number, get_id_for_max_type(type_hi));
793 auto hi = std::upper_bound(lo, field_ents.end(), hi_uid, cmp_uid_hi);
794 if (lo != hi) {
795
796 std::vector<boost::weak_ptr<FieldEntity>> brother_ents_vec;
797
798 for (auto it = lo; it != hi; ++it)
799 if (auto e = it->lock()) {
800 auto side_ptr = e->getSideNumberPtr();
801 if (const auto side = side_ptr->side_number; side >= 0) {
802 const EntityType type = e->getEntType();
803 auto &dat = data.dataOnEntities[type][side];
804 auto &ent_field = dat.getFieldEntities();
805 auto &ent_field_dofs = dat.getFieldDofs();
806 auto &ent_field_data = dat.getFieldData();
807
808 const int brother_side = side_ptr->brother_side_number;
809 if (brother_side != -1)
810 brother_ents_vec.emplace_back(e);
811
812 dat.getBase() = e->getApproxBase();
813 dat.getSpace() = e->getSpace();
814 const int ent_order = e->getMaxOrder();
815 dat.getOrder() =
816 dat.getOrder() > ent_order ? dat.getOrder() : ent_order;
817
818 auto ent_data = e->getEntFieldData();
819 ent_field_data.resize(ent_data.size(), false);
820 noalias(ent_field_data) = ent_data;
821 ent_field_dofs.resize(ent_data.size(), false);
822 std::fill(ent_field_dofs.begin(), ent_field_dofs.end(), nullptr);
823 ent_field.resize(1, false);
824 ent_field[0] = e.get();
825 if (auto cache = e->entityCacheDataDofs.lock()) {
826 for (auto dit = cache->loHi[0]; dit != cache->loHi[1]; ++dit) {
827 ent_field_dofs[(*dit)->getEntDofIdx()] =
828 reinterpret_cast<FEDofEntity *>((*dit).get());
829 }
830 }
831 }
832 }
833
834 for (auto &it : brother_ents_vec) {
835 if (const auto e = it.lock()) {
836 const EntityType type = e->getEntType();
837 const int side = e->getSideNumberPtr()->side_number;
838 const int brother_side = e->getSideNumberPtr()->brother_side_number;
839 auto &dat = data.dataOnEntities[type][side];
840 auto &dat_brother = data.dataOnEntities[type][brother_side];
841 dat_brother.getBase() = dat.getBase();
842 dat_brother.getSpace() = dat.getSpace();
843 dat_brother.getOrder() = dat.getOrder();
844 dat_brother.getFieldData() = dat.getFieldData();
845 dat_brother.getFieldDofs() = dat.getFieldDofs();
846 dat_brother.getFieldEntities() = dat.getFieldEntities();
847 }
848 }
849 }
850 }
851
853}
854
856 const int bit_number, VectorDouble &ent_field_data,
857 VectorDofs &ent_field_dofs, VectorFieldEntities &ent_field) const {
858
860
861 ent_field_data.resize(0, false);
862 ent_field_dofs.resize(0, false);
863 ent_field.resize(0, false);
864
865 auto &field_ents = getDataFieldEnts();
866 const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
867 bit_number, get_id_for_min_type<MBVERTEX>());
868 auto lo = std::lower_bound(field_ents.begin(), field_ents.end(), lo_uid,
869 cmp_uid_lo);
870 if (lo != field_ents.end()) {
871
872 ent_field.resize(field_ents.size(), false);
873 std::fill(ent_field.begin(), ent_field.end(), nullptr);
874
875 const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
876 bit_number, get_id_for_max_type<MBENTITYSET>());
877 auto hi = std::upper_bound(lo, field_ents.end(), hi_uid, cmp_uid_hi);
878 if (lo != hi) {
879
880 int side = 0;
881 for (auto it = lo; it != hi; ++it, ++side)
882 if (auto e = it->lock()) {
883
884 const auto size = e->getNbDofsOnEnt();
885 ent_field_data.resize(size, false);
886 ent_field_dofs.resize(size, false);
887 ent_field[side] = e.get();
888 noalias(ent_field_data) = e->getEntFieldData();
889
890 if (auto cache = e->entityCacheDataDofs.lock()) {
891 for (auto dit = cache->loHi[0]; dit != cache->loHi[1]; ++dit) {
892 ent_field_dofs[(*dit)->getEntDofIdx()] =
893 reinterpret_cast<FEDofEntity *>((*dit).get());
894 }
895 }
896 }
897 }
898 }
899
901}
902
905 const int bit_number) const {
907 if (data.dataOnEntities[MBENTITYSET].size() == 0)
908 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
909 "No space to insert data");
910
912 bit_number, data.dataOnEntities[MBENTITYSET][0].getFieldData(),
913 data.dataOnEntities[MBENTITYSET][0].getFieldDofs(),
914 data.dataOnEntities[MBENTITYSET][0].getFieldEntities());
916}
917
918// ** Face **
919
923 auto &face_nodes = data.facesNodes;
924 auto &face_nodes_order = data.facesNodesOrder;
925 auto &side_table = const_cast<SideNumber_multiIndex &>(
926 numeredEntFiniteElementPtr->getSideNumberTable());
927 const auto ent = numeredEntFiniteElementPtr->getEnt();
928 const auto type = numeredEntFiniteElementPtr->getEntType();
929 const auto nb_faces = CN::NumSubEntities(type, 2);
930 const EntityHandle *conn_ele;
931 int num_nodes_ele;
932 CHKERR mField.get_moab().get_connectivity(ent, conn_ele, num_nodes_ele, true);
933 auto side_ptr_it = side_table.get<1>().lower_bound(
934 boost::make_tuple(CN::TypeDimensionMap[2].first, 0));
935 auto hi_side_ptr_it = side_table.get<1>().upper_bound(
936 boost::make_tuple(CN::TypeDimensionMap[2].second, 100));
937
938 for (; side_ptr_it != hi_side_ptr_it; ++side_ptr_it) {
939 const auto side = (*side_ptr_it)->side_number;
940 const auto sense = (*side_ptr_it)->sense;
941 const auto offset = (*side_ptr_it)->offset;
942
943 EntityType face_type;
944 int nb_nodes_face;
945 auto face_indices =
946 CN::SubEntityVertexIndices(type, 2, side, face_type, nb_nodes_face);
947 face_nodes.resize(nb_faces, nb_nodes_face);
948 face_nodes_order.resize(nb_faces, nb_nodes_face);
949
950 if (sense == 1)
951 for (int n = 0; n != nb_nodes_face; ++n)
952 face_nodes_order(side, n) = (n + offset) % nb_nodes_face;
953 else
954 for (int n = 0; n != nb_nodes_face; ++n)
955 face_nodes_order(side, n) =
956 (nb_nodes_face - (n - offset) % nb_nodes_face) % nb_nodes_face;
957
958 for (int n = 0; n != nb_nodes_face; ++n)
959 face_nodes(side, n) = face_indices[face_nodes_order(side, n)];
960
961#ifndef NDEBUG
962 const auto face_ent = (*side_ptr_it)->ent;
963 auto check = [&]() {
965 const EntityHandle *conn_face;
966 // int nb_nodes_face;
967 CHKERR mField.get_moab().get_connectivity(face_ent, conn_face,
968 nb_nodes_face, true);
969 face_nodes.resize(nb_faces, nb_nodes_face);
970 for (int nn = 0; nn != nb_nodes_face; ++nn) {
971 if (face_nodes(side, nn) !=
972 std::distance(
973 conn_ele,
974 std::find(conn_ele, &conn_ele[num_nodes_ele], conn_face[nn]))) {
975 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
976 "Wrong face numeration");
977 }
978 }
980 };
981 CHKERR check();
982#endif
983 }
984
986}
987
988// ** Space and Base **
989
991 EntitiesFieldData &data) const {
993
994 if (nInTheLoop == 0) {
995 data.sPace.reset();
996 data.bAse.reset();
997 for (EntityType t = MBVERTEX; t != MBMAXTYPE; ++t) {
998 data.spacesOnEntities[t].reset();
999 data.basesOnEntities[t].reset();
1000 }
1001 for (int s = 0; s != LASTSPACE; ++s) {
1002 data.basesOnSpaces[s].reset();
1003 }
1004 }
1005
1006 auto fe_type = numeredEntFiniteElementPtr->getEntType();
1007
1008 if (getDataFieldEntsPtr())
1009 for (auto e : getDataFieldEnts()) {
1010 if (auto ptr = e.lock()) {
1011 // get data from entity
1012 const EntityType type = ptr->getEntType();
1013 if (DefaultElementAdjacency::getDefTypeMap(fe_type, type)) {
1014 const FieldSpace space = ptr->getSpace();
1015 const FieldApproximationBase approx = ptr->getApproxBase();
1016 // set data
1017 data.sPace.set(space);
1018 data.bAse.set(approx);
1019 data.spacesOnEntities[type].set(space);
1020 data.basesOnEntities[type].set(approx);
1021 data.basesOnSpaces[space].set(approx);
1022 }
1023 }
1024 }
1025 else
1026 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1027 "data fields ents not allocated on element");
1028
1030}
1031
1033 const FieldApproximationBase b) {
1035 if (dataOnElement[H1]->bAse.test(b)) {
1036 switch (static_cast<FieldApproximationBase>(b)) {
1037 case NOBASE:
1038 break;
1040 // BERNSTEIN_BEZIER_BASE is not hierarchical base
1041 break;
1046 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1047 "Functions genrating approximation base not defined");
1048
1049 for (int space = H1; space != LASTSPACE; ++space) {
1050 if (dataOnElement[H1]->sPace.test(space) &&
1051 dataOnElement[H1]->bAse.test(b) &&
1052 dataOnElement[H1]->basesOnSpaces[space].test(b)) {
1053 CHKERR getElementPolynomialBase()->getValue(
1054 gaussPts,
1055 boost::make_shared<EntPolynomialBaseCtx>(
1056 *dataOnElement[space], static_cast<FieldSpace>(space),
1057 static_cast<FieldApproximationBase>(b), NOBASE));
1058 }
1059 }
1060 break;
1061 case USER_BASE:
1062 if (!getUserPolynomialBase())
1063 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1064 "Functions generating user approximation base not defined");
1065
1066 for (int space = H1; space != LASTSPACE; ++space)
1067 if (dataOnElement[H1]->sPace.test(space) &&
1068 dataOnElement[H1]->bAse.test(b) &&
1069 dataOnElement[H1]->basesOnSpaces[space].test(b)) {
1070 CHKERR getUserPolynomialBase()->getValue(
1071 gaussPts,
1072 boost::make_shared<EntPolynomialBaseCtx>(
1073 *dataOnElement[space], static_cast<FieldSpace>(space),
1074 static_cast<FieldApproximationBase>(b), NOBASE));
1075 }
1076 break;
1077 default:
1079 "Base <%s> not yet implemented",
1081 }
1082 }
1084}
1085
1088 /// Use the some node base. Node base is usually used for construction other
1089 /// bases.
1090 for (int space = HCURL; space != LASTSPACE; ++space) {
1091 dataOnElement[space]->dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE) =
1092 dataOnElement[H1]->dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE);
1093 dataOnElement[space]->dataOnEntities[MBVERTEX][0].getDiffNSharedPtr(
1094 NOBASE) =
1095 dataOnElement[H1]->dataOnEntities[MBVERTEX][0].getDiffNSharedPtr(
1096 NOBASE);
1097 }
1098 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
1100 static_cast<FieldApproximationBase>(b));
1101 }
1103}
1104
1108
1109 const auto ele_type = numeredEntFiniteElementPtr->getEntType();
1110
1111 auto get_nodal_base_data = [&](EntitiesFieldData &data, auto field_ptr) {
1113 auto &space = data.dataOnEntities[MBVERTEX][0].getSpace();
1114 auto &base = data.dataOnEntities[MBVERTEX][0].getBase();
1115 auto &bb_node_order = data.dataOnEntities[MBVERTEX][0].getBBNodeOrder();
1116
1117 auto &field_ents = getDataFieldEnts();
1118 auto bit_number = field_ptr->getBitNumber();
1119 const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
1120 bit_number, get_id_for_min_type<MBVERTEX>());
1121 auto lo = std::lower_bound(field_ents.begin(), field_ents.end(), lo_uid,
1122 cmp_uid_lo);
1123 if (lo != field_ents.end()) {
1124 const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
1125 bit_number, get_id_for_max_type<MBVERTEX>());
1126 auto hi = std::upper_bound(lo, field_ents.end(), hi_uid, cmp_uid_hi);
1127 if (lo != hi) {
1128
1129 for (auto it = lo; it != hi; ++it)
1130 if (auto first_e = it->lock()) {
1131 space = first_e->getSpace();
1132 base = first_e->getApproxBase();
1133 const int num_nodes = getNumberOfNodes();
1134 bb_node_order.resize(num_nodes, false);
1135 bb_node_order.clear();
1136
1137 std::vector<boost::weak_ptr<FieldEntity>> brother_ents_vec;
1138
1139 for (; it != hi; ++it) {
1140 if (auto e = it->lock()) {
1141 const auto &sn = e->getSideNumberPtr();
1142 const int side_number = sn->side_number;
1143 const int brother_side_number = sn->brother_side_number;
1144 if (brother_side_number != -1)
1145 brother_ents_vec.emplace_back(e);
1146 bb_node_order[side_number] = e->getMaxOrder();
1147 }
1148 }
1149
1150 for (auto &it : brother_ents_vec) {
1151 if (const auto e = it.lock()) {
1152 const auto &sn = e->getSideNumberPtr();
1153 const int side_number = sn->side_number;
1154 const int brother_side_number = sn->brother_side_number;
1155 bb_node_order[brother_side_number] = bb_node_order[side_number];
1156 }
1157 }
1158
1159 break;
1160 }
1161 }
1162 }
1163
1165 };
1166
1167 auto get_entity_base_data = [&](EntitiesFieldData &data, auto field_ptr,
1168 const EntityType type_lo,
1169 const EntityType type_hi) {
1171 for (EntityType t = MBEDGE; t != MBPOLYHEDRON; ++t) {
1172 for (auto &dat : data.dataOnEntities[t]) {
1173 dat.getOrder() = 0;
1174 dat.getBase() = NOBASE;
1175 dat.getSpace() = NOSPACE;
1176 dat.getFieldData().resize(0, false);
1177 dat.getFieldDofs().resize(0, false);
1178 }
1179 }
1180
1181 auto &field_ents = getDataFieldEnts();
1182 auto bit_number = field_ptr->getBitNumber();
1183 const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
1184 bit_number, get_id_for_min_type(type_lo));
1185 auto lo = std::lower_bound(field_ents.begin(), field_ents.end(), lo_uid,
1186 cmp_uid_lo);
1187 if (lo != field_ents.end()) {
1188 const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
1189 bit_number, get_id_for_max_type(type_hi));
1190 auto hi = std::upper_bound(lo, field_ents.end(), hi_uid, cmp_uid_hi);
1191
1192 std::vector<boost::weak_ptr<FieldEntity>> brother_ents_vec;
1193 for (; lo != hi; ++lo) {
1194 if (auto e = lo->lock()) {
1195 if (auto cache = e->entityCacheDataDofs.lock()) {
1196 if (cache->loHi[0] != cache->loHi[1]) {
1197 if (const auto side = e->getSideNumberPtr()->side_number;
1198 side >= 0) {
1199 const EntityType type = e->getEntType();
1200 auto &dat = data.dataOnEntities[type][side];
1201 const int brother_side =
1202 e->getSideNumberPtr()->brother_side_number;
1203 if (brother_side != -1)
1204 brother_ents_vec.emplace_back(e);
1205 dat.getBase() = e->getApproxBase();
1206 dat.getSpace() = e->getSpace();
1207 const auto ent_order = e->getMaxOrder();
1208 dat.getOrder() =
1209 dat.getOrder() > ent_order ? dat.getOrder() : ent_order;
1210 }
1211 }
1212 }
1213 }
1214 }
1215
1216 for (auto &ent_ptr : brother_ents_vec) {
1217 if (auto e = ent_ptr.lock()) {
1218 const EntityType type = e->getEntType();
1219 const int side = e->getSideNumberPtr()->side_number;
1220 const int brother_side = e->getSideNumberPtr()->brother_side_number;
1221 auto &dat = data.dataOnEntities[type][side];
1222 auto &dat_brother = data.dataOnEntities[type][brother_side];
1223 dat_brother.getBase() = dat.getBase();
1224 dat_brother.getSpace() = dat.getSpace();
1225 dat_brother.getOrder() = dat.getOrder();
1226 }
1227 }
1228 }
1230 };
1231
1232 for (auto &e : getDataFieldEnts()) {
1233 if (auto ent_data_ptr = e.lock()) {
1234 if (ent_data_ptr->getApproxBase() == AINSWORTH_BERNSTEIN_BEZIER_BASE) {
1235 auto space = ent_data_ptr->getSpace();
1236 for (EntityType t = MBVERTEX; t != MBPOLYHEDRON; ++t) {
1237 for (auto &dat : (*dataOnElement[space]).dataOnEntities[t]) {
1238 for (auto &ptr : dat.getBBAlphaIndicesByOrderArray())
1239 ptr.reset();
1240 for (auto &ptr : dat.getBBNByOrderArray())
1241 ptr.reset();
1242 for (auto &ptr : dat.getBBDiffNByOrderArray())
1243 ptr.reset();
1244 }
1245 }
1246 }
1247 }
1248 }
1249
1250 auto check_space = [&](const auto space) {
1251 switch (space) {
1252 case H1:
1253 for (auto t = MBVERTEX; t <= ele_type; ++t) {
1254 if (dataOnElement[H1]->spacesOnEntities[t].test(H1))
1255 return true;
1256 }
1257 return false;
1258 case HCURL:
1259 for (auto t = MBEDGE; t <= ele_type; ++t) {
1260 if (dataOnElement[HCURL]->spacesOnEntities[t].test(HCURL))
1261 return true;
1262 }
1263 return false;
1264 case HDIV:
1265 for (auto t = MBTRI; t <= ele_type; ++t) {
1266 if (dataOnElement[HDIV]->spacesOnEntities[t].test(HDIV))
1267 return true;
1268 }
1269 return false;
1270 case L2:
1271 return dataOnElement[L2]->spacesOnEntities[ele_type].test(L2);
1272 break;
1273 default:
1274 THROW_MESSAGE("Not implemented");
1275 }
1276 };
1277
1278 std::set<string> fields_list;
1279 for (auto &e : getDataFieldEnts()) {
1280 if (auto ent_data_ptr = e.lock()) {
1281 if (ent_data_ptr->getApproxBase() == AINSWORTH_BERNSTEIN_BEZIER_BASE) {
1282 auto field_name = ent_data_ptr->getName();
1283 if (fields_list.find(field_name) == fields_list.end()) {
1284 auto field_ptr = ent_data_ptr->getFieldRawPtr();
1285 auto space = ent_data_ptr->getSpace();
1286 CHKERR get_nodal_base_data(*dataOnElement[space], field_ptr);
1287 CHKERR get_entity_base_data(*dataOnElement[space], field_ptr, MBEDGE,
1288 ele_type);
1289 if (check_space(space)) {
1290 CHKERR getElementPolynomialBase()->getValue(
1291 gaussPts, boost::make_shared<EntPolynomialBaseCtx>(
1292 *dataOnElement[space], field_name,
1293 static_cast<FieldSpace>(space),
1295 fields_list.insert(field_name);
1296 }
1297 }
1298 }
1299 }
1300 }
1301
1303};
1304
1307
1308 // Data on elements for proper spaces
1309 for (int space = H1; space != LASTSPACE; ++space) {
1310 dataOnElement[space]->setElementType(type);
1311 derivedDataOnElement[space]->setElementType(type);
1312 }
1313
1315}
1316
1317#define FUNCTION_NAME_WITH_OP_NAME(OP) \
1318 std::ostringstream ss; \
1319 ss << "(Calling user data operator " \
1320 << boost::typeindex::type_id_runtime(OP).pretty_name() << " rowField " \
1321 << (OP).rowFieldName << " colField " << (OP).colFieldName << ") "
1322
1323#define CATCH_OP_ERRORS(OP) \
1324 catch (MoFEMExceptionInitial const &ex) { \
1325 FUNCTION_NAME_WITH_OP_NAME(OP) << PETSC_FUNCTION_NAME; \
1326 return PetscError(PETSC_COMM_SELF, ex.lINE, ss.str().c_str(), __FILE__, \
1327 ex.errorCode, PETSC_ERROR_INITIAL, ex.what()); \
1328 } \
1329 catch (MoFEMExceptionRepeat const &ex) { \
1330 FUNCTION_NAME_WITH_OP_NAME(OP) << PETSC_FUNCTION_NAME; \
1331 return PetscError(PETSC_COMM_SELF, ex.lINE, ss.str().c_str(), __FILE__, \
1332 ex.errorCode, PETSC_ERROR_REPEAT, " "); \
1333 } \
1334 catch (MoFEMException const &ex) { \
1335 FUNCTION_NAME_WITH_OP_NAME(OP) << ex.errorMessage; \
1336 SETERRQ(PETSC_COMM_SELF, ex.errorCode, ss.str().c_str()); \
1337 } \
1338 catch (std::exception const &ex) { \
1339 std::string message("Error: " + std::string(ex.what()) + " at " + \
1340 boost::lexical_cast<std::string>(__LINE__) + " : " + \
1341 std::string(__FILE__) + " in " + \
1342 std::string(PETSC_FUNCTION_NAME)); \
1343 FUNCTION_NAME_WITH_OP_NAME(OP) << message; \
1344 SETERRQ(PETSC_COMM_SELF, MOFEM_STD_EXCEPTION_THROW, ss.str().c_str()); \
1345 }
1346
1349
1350 using UDO = UserDataOperator;
1351
1352 std::array<std::string, 2> field_name;
1353 std::array<const Field *, 3> field_struture;
1354 std::array<int, 2>
1355 field_id; // note the this is field bit number (nor field bit)
1356 std::array<FieldSpace, 2> space;
1357 std::array<FieldApproximationBase, 2> base;
1358
1359 constexpr std::array<UDO::OpType, 2> types = {UDO::OPROW, UDO::OPCOL};
1360 std::array<int, 2> last_eval_field_id = {0, 0};
1361
1362 std::array<boost::shared_ptr<EntitiesFieldData>, 2> op_data;
1363
1364 auto swap_bases = [&](auto &op) {
1366 for (size_t ss = 0; ss != 2; ++ss) {
1367 if (op.getOpType() & types[ss] || op.getOpType() & UDO::OPROWCOL) {
1368 switch (base[ss]) {
1370 CHKERR op_data[ss]->baseSwap(field_name[ss],
1372 default:
1373 break;
1374 }
1375 }
1376 }
1377
1379 };
1380
1381 const EntityType type = numeredEntFiniteElementPtr->getEntType();
1382
1383 // evaluate entity data only, no field specific data provided or known
1384 auto evaluate_op_space = [&](auto &op) {
1386
1387 // reseat all data which all field dependent
1388 dataOnElement[op.sPace]->resetFieldDependentData();
1389 std::fill(last_eval_field_id.begin(), last_eval_field_id.end(), 0);
1390
1391 switch (op.sPace) {
1392 case NOSPACE:
1393 try {
1394 CHKERR op.doWork(
1395 0, MBENTITYSET,
1396 dataOnElement[NOSPACE]->dataOnEntities[MBENTITYSET][0]);
1397 }
1398 CATCH_OP_ERRORS(op);
1399 break;
1400 case NOFIELD:
1401 case H1:
1402 case HCURL:
1403 case HDIV:
1404 case L2:
1405 try {
1406
1407 for (EntityType t = MBVERTEX; t != MBMAXTYPE; ++t) {
1408 for (auto &e : dataOnElement[op.sPace]->dataOnEntities[t]) {
1409 e.getSpace() = op.sPace;
1410 }
1411 }
1412
1413 CHKERR op.opRhs(*dataOnElement[op.sPace], false);
1414 }
1415 CATCH_OP_ERRORS(op);
1416 break;
1417 default:
1418 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1419 "Not implemented for this space", op.sPace);
1420 }
1421
1423 };
1424
1425 // set entity data
1426 auto set_op_entityties_data = [&](auto ss, auto &op) {
1428
1429#ifndef NDEBUG
1430 if ((op.getNumeredEntFiniteElementPtr()->getBitFieldIdData() &
1432 .none()) {
1433 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1434 "no data field < %s > on finite element < %s >",
1435 field_name[ss].c_str(), getFEName().c_str());
1436 }
1437#endif
1438
1439 op_data[ss] =
1440 !ss ? dataOnElement[space[ss]] : derivedDataOnElement[space[ss]];
1441
1442 for (auto &data : op_data[ss]->dataOnEntities[MBENTITYSET]) {
1443 CHKERR data.resetFieldDependentData();
1444 }
1445
1446 auto get_data_for_nodes = [&]() {
1448 if (!ss)
1449 CHKERR getRowNodesIndices(*op_data[ss], field_id[ss]);
1450 else
1451 CHKERR getColNodesIndices(*op_data[ss], field_id[ss]);
1452 CHKERR getNodesFieldData(*op_data[ss], field_id[ss]);
1454 };
1455
1456 // get data on entities but not nodes
1457 auto get_data_for_entities = [&]() {
1459 CHKERR getEntityFieldData(*op_data[ss], field_id[ss], MBEDGE);
1460 if (!ss)
1461 CHKERR getEntityRowIndices(*op_data[ss], field_id[ss], MBEDGE);
1462 else
1463 CHKERR getEntityColIndices(*op_data[ss], field_id[ss], MBEDGE);
1465 };
1466
1467 auto get_data_for_meshset = [&]() {
1469 if (!ss) {
1470 CHKERR getNoFieldRowIndices(*op_data[ss], field_id[ss]);
1471 } else {
1472 CHKERR getNoFieldColIndices(*op_data[ss], field_id[ss]);
1473 }
1474 CHKERR getNoFieldFieldData(*op_data[ss], field_id[ss]);
1476 };
1477
1478 switch (space[ss]) {
1479 case H1:
1480 CHKERR get_data_for_nodes();
1481 case HCURL:
1482 case HDIV:
1483 CHKERR get_data_for_entities();
1484 break;
1485 case L2:
1486 switch (type) {
1487 case MBVERTEX:
1488 CHKERR get_data_for_nodes();
1489 break;
1490 default:
1491 CHKERR get_data_for_entities();
1492 }
1493 break;
1494 case NOFIELD:
1495 // if (!getNinTheLoop()) {
1496 // NOFIELD data are the same for each element, can be
1497 // retrieved only once
1498 CHKERR get_data_for_meshset();
1499 // }
1500 break;
1501 default:
1502 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1503 "not implemented for this space < %s >",
1504 FieldSpaceNames[space[ss]]);
1505 }
1507 };
1508
1509 // evalate operators with field data
1510 auto evaluate_op_for_fields = [&](auto &op) {
1512
1513 if (op.getOpType() & UDO::OPROW) {
1514 try {
1515 CHKERR op.opRhs(*op_data[0], false);
1516 }
1517 CATCH_OP_ERRORS(op);
1518 }
1519
1520 if (op.getOpType() & UDO::OPCOL) {
1521 try {
1522 CHKERR op.opRhs(*op_data[1], false);
1523 }
1524 CATCH_OP_ERRORS(op);
1525 }
1526
1527 if (op.getOpType() & UDO::OPROWCOL) {
1528 try {
1529 CHKERR op.opLhs(*op_data[0], *op_data[1]);
1530 }
1531 CATCH_OP_ERRORS(op);
1532 }
1534 };
1535
1536 // Collect bit ref level on all entities, fields and spaces
1538
1539 auto oit = opPtrVector.begin();
1540 auto hi_oit = opPtrVector.end();
1541
1542 // interate over all operators on element
1543 for (; oit != hi_oit; oit++) {
1544
1545 try {
1546
1547 CHKERR oit->setPtrFE(this);
1548
1549 if ((oit->opType & UDO::OPSPACE) == UDO::OPSPACE) {
1550
1551 // run operator for space but specific field
1552 CHKERR evaluate_op_space(*oit);
1553
1554 } else if (
1555
1556 (oit->opType & (UDO::OPROW | UDO::OPCOL | UDO::OPROWCOL)) ==
1557 oit->opType
1558
1559 ) {
1560
1561 field_name[0] = oit->rowFieldName;
1562 field_name[1] = oit->colFieldName;
1563
1564 for (size_t ss = 0; ss != 2; ss++) {
1565
1566#ifndef NDEBUG
1567 if (field_name[ss].empty()) {
1568 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1569 "Not set Field name in operator %d (0-row, 1-column) in "
1570 "operator %s",
1571 ss,
1572 (boost::typeindex::type_id_runtime(*oit).pretty_name())
1573 .c_str());
1574 }
1575#endif
1576
1577 field_struture[ss] = mField.get_field_structure(field_name[ss]);
1578 field_id[ss] = field_struture[ss]->getBitNumber();
1579 space[ss] = field_struture[ss]->getSpace();
1580 base[ss] = field_struture[ss]->getApproxBase();
1581 }
1582
1583 // not that if field name do not change between operators, entity field
1584 // data are nor rebuild
1585 for (size_t ss = 0; ss != 2; ss++) {
1586
1587 if (oit->getOpType() & types[ss] ||
1588 oit->getOpType() & UDO::OPROWCOL) {
1589 if (last_eval_field_id[ss] != field_id[ss]) {
1590 CHKERR set_op_entityties_data(ss, *oit);
1591 last_eval_field_id[ss] = field_id[ss];
1592 }
1593 }
1594 }
1595
1596 CHKERR swap_bases(*oit);
1597
1598 // run operator for given field or row, column or both
1599 CHKERR evaluate_op_for_fields(*oit);
1600
1601 CHKERR swap_bases(*oit);
1602
1603 } else {
1604
1605 for (int i = 0; i != 5; ++i)
1606 if (oit->opType & (1 << i))
1607 MOFEM_LOG("SELF", Sev::error) << UDO::OpTypeNames[i];
1608 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1609 "Impossible operator type");
1610 }
1611 }
1612 CATCH_OP_ERRORS(*oit);
1613 }
1615}
1616
1617const char *const ForcesAndSourcesCore::UserDataOperator::OpTypeNames[] = {
1618 "OPROW", " OPCOL", "OPROWCOL", "OPSPACE", "OPLAST"};
1619
1620MoFEMErrorCode ForcesAndSourcesCore::UserDataOperator::getProblemRowIndices(
1621 const std::string field_name, const EntityType type, const int side,
1622 VectorInt &indices) const {
1624 if (ptrFE == NULL)
1625 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1626
1627 switch (type) {
1628 case MBVERTEX:
1630 break;
1631 default:
1632 CHKERR ptrFE->getProblemTypeRowIndices(field_name, type, side, indices);
1633 }
1635}
1636
1637MoFEMErrorCode ForcesAndSourcesCore::UserDataOperator::getProblemColIndices(
1638 const std::string field_name, const EntityType type, const int side,
1639 VectorInt &indices) const {
1641 if (ptrFE == NULL)
1642 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1643
1644 switch (type) {
1645 case MBVERTEX:
1646 CHKERR ptrFE->getProblemNodesColIndices(field_name, indices);
1647 break;
1648 default:
1649 CHKERR ptrFE->getProblemTypeColIndices(field_name, type, side, indices);
1650 }
1652}
1653
1657 sidePtrFE = const_cast<ForcesAndSourcesCore *>(side_fe_ptr);
1659}
1660
1662 const ForcesAndSourcesCore *refine_fe_ptr) {
1664 refinePtrFE = const_cast<ForcesAndSourcesCore *>(refine_fe_ptr);
1666}
1667
1668MoFEMErrorCode ForcesAndSourcesCore::UserDataOperator::loopSide(
1669 const string &fe_name, ForcesAndSourcesCore *side_fe, const size_t side_dim,
1670 const EntityHandle ent_for_side, const int verb,
1671 const LogManager::SeverityLevel sev, AdjCache *adj_cache) {
1673
1674 const auto *problem_ptr = getFEMethod()->problemPtr;
1675 auto &numered_fe =
1676 problem_ptr->numeredFiniteElementsPtr->get<Unique_mi_tag>();
1677
1678 auto fe_miit = ptrFE->mField.get_finite_elements()
1680 .find(fe_name);
1681 if (fe_miit != ptrFE->mField.get_finite_elements()
1683 .end()) {
1684
1685 const auto ent = ent_for_side ? ent_for_side : getFEEntityHandle();
1686
1687 side_fe->feName = fe_name;
1688
1689 CHKERR side_fe->setSideFEPtr(ptrFE);
1690 CHKERR side_fe->copyBasicMethod(*getFEMethod());
1691 CHKERR side_fe->copyPetscData(*getFEMethod());
1692 CHKERR side_fe->copyKsp(*getFEMethod());
1693 CHKERR side_fe->copySnes(*getFEMethod());
1694 CHKERR side_fe->copyTs(*getFEMethod());
1695
1696 side_fe->cacheWeakPtr = getFEMethod()->cacheWeakPtr;
1697
1698 CHKERR side_fe->preProcess();
1699
1700 std::vector<boost::weak_ptr<NumeredEntFiniteElement>> fe_vec;
1701 auto get_numered_fe_ptr = [&](auto &fe_uid, Range &&adjacent_ents)
1702 -> std::vector<boost::weak_ptr<NumeredEntFiniteElement>> & {
1703 fe_vec.reserve(adjacent_ents.size());
1704 for (auto fe_ent : adjacent_ents) {
1705 auto miit = numered_fe.find(
1707 if (miit != numered_fe.end()) {
1708 fe_vec.emplace_back(*miit);
1709 }
1710 }
1711 return fe_vec;
1712 };
1713
1714 auto get_bit_entity_adjacency = [&]() {
1715 Range adjacent_ents;
1716 CHKERR ptrFE->mField.getInterface<BitRefManager>()->getAdjacenciesAny(
1717 ent, side_dim, adjacent_ents);
1718 return adjacent_ents;
1719 };
1720
1721 auto get_adj = [&](auto &fe_uid)
1722 -> std::vector<boost::weak_ptr<NumeredEntFiniteElement>> & {
1723 if (adj_cache) {
1724 try {
1725 return (*adj_cache).at(ent);
1726 } catch (const std::out_of_range &) {
1727 return (*adj_cache)[ent] =
1728 get_numered_fe_ptr(fe_uid, get_bit_entity_adjacency());
1729 }
1730 } else {
1731 return get_numered_fe_ptr(fe_uid, get_bit_entity_adjacency());
1732 }
1733 };
1734
1735 auto adj = get_adj((*fe_miit)->getFEUId());
1736
1737 int nn = 0;
1738 side_fe->loopSize = adj.size();
1739 for (auto fe_weak_ptr : adj) {
1740 if (auto fe_ptr = fe_weak_ptr.lock()) {
1741 if (verb >= VERBOSE)
1742 MOFEM_LOG("SELF", sev) << "Side finite element "
1743 << "(" << nn << "): " << *fe_ptr;
1744 side_fe->nInTheLoop = nn;
1745 side_fe->numeredEntFiniteElementPtr = fe_ptr;
1746 CHKERR (*side_fe)();
1747 }
1748 ++nn;
1749 }
1750
1751 CHKERR side_fe->postProcess();
1752 }
1753
1755}
1756
1757MoFEMErrorCode ForcesAndSourcesCore::UserDataOperator::loopThis(
1758 const string &fe_name, ForcesAndSourcesCore *this_fe, const int verb,
1759 const LogManager::SeverityLevel sev) {
1761
1762 if (verb >= VERBOSE)
1763 MOFEM_LOG("SELF", sev) << "This finite element: "
1764 << *getNumeredEntFiniteElementPtr();
1765
1766 this_fe->feName = fe_name;
1767
1768 CHKERR this_fe->setRefineFEPtr(ptrFE);
1769 CHKERR this_fe->copyBasicMethod(*getFEMethod());
1770 CHKERR this_fe->copyPetscData(*getFEMethod());
1771 CHKERR this_fe->copyKsp(*getFEMethod());
1772 CHKERR this_fe->copySnes(*getFEMethod());
1773 CHKERR this_fe->copyTs(*getFEMethod());
1774
1775 this_fe->cacheWeakPtr = getFEMethod()->cacheWeakPtr;
1776
1777 CHKERR this_fe->preProcess();
1778
1779 this_fe->nInTheLoop = getNinTheLoop();
1780 this_fe->numeredEntFiniteElementPtr = getNumeredEntFiniteElementPtr();
1781
1782 CHKERR (*this_fe)();
1783
1784 CHKERR this_fe->postProcess();
1785
1787}
1788
1789MoFEMErrorCode ForcesAndSourcesCore::UserDataOperator::loopParent(
1790 const string &fe_name, ForcesAndSourcesCore *parent_fe, const int verb,
1791 const LogManager::SeverityLevel sev) {
1793
1794 auto &fes =
1795 ptrFE->mField.get_finite_elements()->get<FiniteElement_name_mi_tag>();
1796 auto fe_miit = fes.find(fe_name);
1797 if (fe_miit != fes.end()) {
1798
1799 const auto *problem_ptr = getFEMethod()->problemPtr;
1800 auto &numered_fe =
1801 problem_ptr->numeredFiniteElementsPtr->get<Unique_mi_tag>();
1802
1803
1804 parent_fe->feName = fe_name;
1805 CHKERR parent_fe->setRefineFEPtr(ptrFE);
1806 CHKERR parent_fe->copyBasicMethod(*getFEMethod());
1807 CHKERR parent_fe->copyPetscData(*getFEMethod());
1808 CHKERR parent_fe->copyKsp(*getFEMethod());
1809 CHKERR parent_fe->copySnes(*getFEMethod());
1810 CHKERR parent_fe->copyTs(*getFEMethod());
1811
1812 parent_fe->cacheWeakPtr = getFEMethod()->cacheWeakPtr;
1813
1814 const auto parent_ent = getNumeredEntFiniteElementPtr()->getParentEnt();
1815 auto miit = numered_fe.find(EntFiniteElement::getLocalUniqueIdCalculate(
1816 parent_ent, (*fe_miit)->getFEUId()));
1817 if (miit != numered_fe.end()) {
1818 if (verb >= VERBOSE)
1819 MOFEM_LOG("SELF", sev) << "Parent finite element: " << **miit;
1820 parent_fe->loopSize = 1;
1821 parent_fe->nInTheLoop = 0;
1822 parent_fe->numeredEntFiniteElementPtr = *miit;
1823 CHKERR parent_fe->preProcess();
1824 CHKERR (*parent_fe)();
1825 CHKERR parent_fe->postProcess();
1826 } else {
1827 if (verb >= VERBOSE)
1828 MOFEM_LOG("SELF", sev) << "Parent finite element: no parent";
1829 parent_fe->loopSize = 0;
1830 parent_fe->nInTheLoop = 0;
1831 CHKERR parent_fe->preProcess();
1832 CHKERR parent_fe->postProcess();
1833 }
1834 }
1835
1837}
1838
1839MoFEMErrorCode ForcesAndSourcesCore::UserDataOperator::loopChildren(
1840 const string &fe_name, ForcesAndSourcesCore *child_fe, const int verb,
1841 const LogManager::SeverityLevel sev) {
1843
1844 auto fe_miit = ptrFE->mField.get_finite_elements()
1846 .find(fe_name);
1847 if (fe_miit != ptrFE->mField.get_finite_elements()
1849 .end()) {
1850
1851 const auto *problem_ptr = getFEMethod()->problemPtr;
1852 auto &ref_ents = *getPtrFE()->mField.get_ref_ents();
1853 auto &numered_fe =
1854 problem_ptr->numeredFiniteElementsPtr->get<Unique_mi_tag>();
1855
1856 const auto parent_ent = getNumeredEntFiniteElementPtr()->getEnt();
1857 const auto parent_type = getNumeredEntFiniteElementPtr()->getEntType();
1858 auto range =
1859 ref_ents.get<Composite_ParentEnt_And_EntType_mi_tag>().equal_range(
1860 boost::make_tuple(parent_type, parent_ent));
1861
1862 if (auto size = std::distance(range.first, range.second)) {
1863
1864 std::vector<EntityHandle> childs_vec;
1865 childs_vec.reserve(size);
1866 for (; range.first != range.second; ++range.first)
1867 childs_vec.emplace_back((*range.first)->getEnt());
1868
1869 Range childs;
1870
1871 if ((childs_vec.back() - childs_vec.front() + 1) == size)
1872 childs = Range(childs_vec.front(), childs_vec.back());
1873 else
1874 childs.insert_list(childs_vec.begin(), childs_vec.end());
1875
1876 child_fe->feName = fe_name;
1877
1878 CHKERR child_fe->setRefineFEPtr(ptrFE);
1879 CHKERR child_fe->copyBasicMethod(*getFEMethod());
1880 CHKERR child_fe->copyPetscData(*getFEMethod());
1881 CHKERR child_fe->copyKsp(*getFEMethod());
1882 CHKERR child_fe->copySnes(*getFEMethod());
1883 CHKERR child_fe->copyTs(*getFEMethod());
1884
1885 child_fe->cacheWeakPtr = getFEMethod()->cacheWeakPtr;
1886
1887 CHKERR child_fe->preProcess();
1888
1889 int nn = 0;
1890 child_fe->loopSize = size;
1891
1892 for (auto p = childs.pair_begin(); p != childs.pair_end(); ++p) {
1893
1894 auto miit =
1895 numered_fe.lower_bound(EntFiniteElement::getLocalUniqueIdCalculate(
1896 p->first, (*fe_miit)->getFEUId()));
1897 auto hi_miit =
1898 numered_fe.upper_bound(EntFiniteElement::getLocalUniqueIdCalculate(
1899 p->second, (*fe_miit)->getFEUId()));
1900
1901 for (; miit != hi_miit; ++miit) {
1902
1903 if (verb >= VERBOSE)
1904 MOFEM_LOG("SELF", sev) << "Child finite element "
1905 << "(" << nn << "): " << **miit;
1906
1907 child_fe->nInTheLoop = nn++;
1908 child_fe->numeredEntFiniteElementPtr = *miit;
1909 CHKERR (*child_fe)();
1910 }
1911 }
1912 }
1913
1914 CHKERR child_fe->postProcess();
1915 }
1916
1918}
1919
1920int ForcesAndSourcesCore::getRule(int order_row, int order_col,
1921 int order_data) {
1922 return getRuleHook ? getRuleHook(order_row, order_col, order_data)
1923 : getRule(order_data);
1924}
1925
1927 int order_data) {
1928 return setRuleHook ? setRuleHook(this, order_row, order_col, order_data)
1929 : setGaussPts(order_data);
1930}
1931
1933
1934/** \deprecated setGaussPts(int row_order, int col_order, int data order);
1935 */
1938 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Sorry, not implemented");
1940}
1941
1942ForcesAndSourcesCore::UserDataOperator::UserDataOperator(const FieldSpace space,
1943 const char type,
1944 const bool symm)
1945 : DataOperator(symm), opType(type), sPace(space), ptrFE(nullptr) {}
1946
1947ForcesAndSourcesCore::UserDataOperator::UserDataOperator(
1948 const std::string field_name, const char type, const bool symm)
1949 : DataOperator(symm), opType(type), rowFieldName(field_name),
1950 colFieldName(field_name), sPace(LASTSPACE), ptrFE(nullptr) {}
1951
1952ForcesAndSourcesCore::UserDataOperator::UserDataOperator(
1953 const std::string row_field_name, const std::string col_field_name,
1954 const char type, const bool symm)
1955 : DataOperator(symm), opType(type), rowFieldName(row_field_name),
1956 colFieldName(col_field_name), sPace(LASTSPACE), ptrFE(nullptr) {}
1957
1960 if (preProcessHook) {
1961 ierr = preProcessHook();
1962 CHKERRG(ierr);
1963 }
1965}
1968 if (operatorHook) {
1969 ierr = operatorHook();
1970 CHKERRG(ierr);
1971 } else {
1972#ifndef NDEBUG
1973 MOFEM_LOG("SELF", Sev::warning)
1974 << "No method operator() overloaded on element entity on finite "
1975 "element <"
1976 << boost::typeindex::type_id_runtime(*this).pretty_name() << ">";
1977#endif
1978 }
1980}
1983 if (postProcessHook) {
1985 CHKERRG(ierr);
1986 }
1988}
1989
1991ForcesAndSourcesCore::UserDataOperator::setPtrFE(ForcesAndSourcesCore *ptr) {
1993 if (!(ptrFE = dynamic_cast<ForcesAndSourcesCore *>(ptr)))
1994 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1995 "User operator and finite element do not work together");
1997}
1998
1999} // 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.
boost::weak_ptr< CacheTuple > cacheWeakPtr
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...
std::map< EntityHandle, std::vector< boost::weak_ptr< NumeredEntFiniteElement > > > AdjCache
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