v0.13.1
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
7
8#ifdef __cplusplus
9extern "C" {
10#endif
11#include <cblas.h>
12#include <lapack_wrap.h>
13// #include <gm_rule.h>
14#ifdef __cplusplus
15}
16#endif
17
18namespace MoFEM {
19
20static auto cmp_uid_lo(const boost::weak_ptr<FieldEntity> &a, const UId &b) {
21 if (auto a_ptr = a.lock()) {
22 if (a_ptr->getLocalUniqueId() < b)
23 return true;
24 else
25 return false;
26 } else {
27 return false;
28 }
29}
30
31static auto cmp_uid_hi(const UId &b, const boost::weak_ptr<FieldEntity> &a) {
32 if (auto a_ptr = a.lock()) {
33 if (b < a_ptr->getLocalUniqueId())
34 return true;
35 else
36 return false;
37 } else {
38 return true;
39 }
40}
41
43 :
44
45 mField(m_field), getRuleHook(0), setRuleHook(0),
46 dataOnElement{
47
48 boost::make_shared<EntitiesFieldData>(MBENTITYSET), // NOSPACE,
49 boost::make_shared<EntitiesFieldData>(MBENTITYSET), // NOFIELD
50 boost::make_shared<EntitiesFieldData>(MBENTITYSET), // H1
51 boost::make_shared<EntitiesFieldData>(MBENTITYSET), // HCURL
52 boost::make_shared<EntitiesFieldData>(MBENTITYSET), // HDIV
53 boost::make_shared<EntitiesFieldData>(MBENTITYSET) // L2
54
55 },
56 derivedDataOnElement{
57
58 nullptr,
59 boost::make_shared<DerivedEntitiesFieldData>(
60 dataOnElement[NOFIELD]), // NOFIELD
61 boost::make_shared<DerivedEntitiesFieldData>(dataOnElement[H1]), // H1
62 boost::make_shared<DerivedEntitiesFieldData>(
63 dataOnElement[HCURL]), // HCURL
64 boost::make_shared<DerivedEntitiesFieldData>(
65 dataOnElement[HDIV]), // HDIV
66 boost::make_shared<DerivedEntitiesFieldData>(dataOnElement[L2]) // L2
67
68 },
69 dataNoField(*dataOnElement[NOFIELD].get()),
70 dataH1(*dataOnElement[H1].get()), dataHcurl(*dataOnElement[HCURL].get()),
71 dataHdiv(*dataOnElement[HDIV].get()), dataL2(*dataOnElement[L2].get()),
72 lastEvaluatedElementEntityType(MBMAXTYPE), sidePtrFE(nullptr),
73 refinePtrFE(nullptr) {
74
75 dataOnElement[NOSPACE]->dataOnEntities[MBENTITYSET].push_back(
77
78 dataOnElement[NOFIELD]->dataOnEntities[MBENTITYSET].push_back(
80 derivedDataOnElement[NOFIELD]->dataOnEntities[MBENTITYSET].push_back(
82}
83
84// ** Sense **
85
87 const EntityType type,
88 boost::ptr_vector<EntitiesFieldData::EntData> &data) const {
90
91 auto &side_table = numeredEntFiniteElementPtr->getSideNumberTable().get<0>();
92 auto sit = side_table.lower_bound(get_id_for_min_type(type));
93 if (sit != side_table.end()) {
94 auto hi_sit = side_table.upper_bound(get_id_for_max_type(type));
95 for (; sit != hi_sit; ++sit) {
96 if (const auto side_number = (*sit)->side_number; side_number >= 0) {
97 const int brother_side_number = (*sit)->brother_side_number;
98 const int sense = (*sit)->sense;
99
100 data[side_number].getSense() = sense;
101 if (brother_side_number != -1)
102 data[brother_side_number].getSense() = sense;
103 }
104 }
105 }
107}
108
109// ** Order **
110
111template <typename ENTMULTIINDEX>
112static inline int getMaxOrder(const ENTMULTIINDEX &multi_index) {
113 int max_order = 0;
114 for (auto ent_field_weak_ptr : multi_index)
115 if (auto e = ent_field_weak_ptr.lock()) {
116 const int order = e->getMaxOrder();
117 max_order = (max_order < order) ? order : max_order;
118 }
119 return max_order;
120}
121
123 int max_order = 0;
124 for (auto e : getDataFieldEnts()) {
125 if (auto ptr = e.lock()) {
126 const int order = ptr->getMaxOrder();
127 max_order = (max_order < order) ? order : max_order;
128 }
129 }
130 return max_order;
131}
132
135}
136
139}
140
142 const EntityType type, const FieldSpace space,
143 boost::ptr_vector<EntitiesFieldData::EntData> &data) const {
145
146 auto set_order = [&]() {
148 auto &side_table = numeredEntFiniteElementPtr->getSideNumberTable();
149
150 for (unsigned int s = 0; s != data.size(); ++s)
151 data[s].getOrder() = 0;
152
153 const FieldEntity_vector_view &data_field_ent = getDataFieldEnts();
154
155 for (auto r = fieldsPtr->get<BitFieldId_space_mi_tag>().equal_range(space);
156 r.first != r.second; ++r.first) {
157
158 const auto field_bit_number = (*r.first)->getBitNumber();
159 const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
160 field_bit_number, get_id_for_min_type(type));
161 auto lo = std::lower_bound(data_field_ent.begin(), data_field_ent.end(),
162 lo_uid, cmp_uid_lo);
163 if (lo != data_field_ent.end()) {
164 const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
165 field_bit_number, get_id_for_max_type(type));
166 auto hi =
167 std::upper_bound(lo, data_field_ent.end(), hi_uid, cmp_uid_hi);
168 for (; lo != hi; ++lo) {
169
170 if (auto ptr = lo->lock()) {
171
172 auto &e = *ptr;
173 auto sit = side_table.find(e.getEnt());
174 if (sit != side_table.end()) {
175 auto &side = *sit;
176 if (const auto side_number = side->side_number; side_number >= 0) {
177 ApproximationOrder ent_order = e.getMaxOrder();
178 auto &dat = data[side_number];
179 dat.getOrder() =
180 dat.getOrder() > ent_order ? dat.getOrder() : ent_order;
181 }
182 } else
183 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
184 "Entity on side of the element not found");
185 }
186 }
187 }
188 }
189
191 };
192
193 auto set_order_on_brother = [&]() {
195 auto &side_table =
196 numeredEntFiniteElementPtr->getSideNumberTable().get<0>();
197 auto sit = side_table.lower_bound(get_id_for_min_type(type));
198 if (sit != side_table.end()) {
199 auto hi_sit = side_table.upper_bound(get_id_for_max_type(type));
200 for (; sit != hi_sit; ++sit) {
201 const int brother_side_number = (*sit)->brother_side_number;
202 if (brother_side_number != -1) {
203 const int side_number = (*sit)->side_number;
204 data[brother_side_number].getOrder() = data[side_number].getOrder();
205 }
206 }
207 }
209 };
210
211 CHKERR set_order();
212 CHKERR set_order_on_brother();
213
215}
216
217// ** Indices **
218
219template <typename EXTRACTOR>
221 const int bit_number, FieldEntity_vector_view &ents_field,
222 VectorInt &nodes_indices, VectorInt &local_nodes_indices,
223 EXTRACTOR &&extractor) const {
225
226 auto field_it = fieldsPtr->get<BitFieldId_mi_tag>().find(
227 BitFieldId().set(bit_number - 1));
228 if (field_it != fieldsPtr->get<BitFieldId_mi_tag>().end()) {
229
230#ifndef NDEBUG
231 if ((*field_it)->getBitNumber() != bit_number)
232 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong bit number");
233#endif
234 const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
235 bit_number, get_id_for_min_type<MBVERTEX>());
236 auto lo = std::lower_bound(ents_field.begin(), ents_field.end(), lo_uid,
237 cmp_uid_lo);
238 if (lo != ents_field.end()) {
239 const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
240 bit_number, get_id_for_max_type<MBVERTEX>());
241 auto hi = std::upper_bound(lo, ents_field.end(), hi_uid, cmp_uid_hi);
242
243 const int num_nodes = getNumberOfNodes();
244 const int nb_dofs_on_vert = (*field_it)->getNbOfCoeffs();
245 const int max_nb_dofs = nb_dofs_on_vert * num_nodes;
246
247 int nb_dofs = 0;
248 for (auto it = lo; it != hi; ++it) {
249 if (auto e = it->lock()) {
250 if (auto cache = extractor(e).lock()) {
251 if (cache->loHi[0] != cache->loHi[1]) {
252 nb_dofs += std::distance(cache->loHi[0], cache->loHi[1]);
253 break;
254 }
255 }
256 }
257 }
258
259 if (nb_dofs) {
260 nodes_indices.resize(max_nb_dofs, false);
261 local_nodes_indices.resize(max_nb_dofs, false);
262 } else {
263 nodes_indices.resize(0, false);
264 local_nodes_indices.resize(0, false);
265 }
266
267 if (nb_dofs != max_nb_dofs) {
268 std::fill(nodes_indices.begin(), nodes_indices.end(), -1);
269 std::fill(local_nodes_indices.begin(), local_nodes_indices.end(), -1);
270 }
271
272 for (auto it = lo; it != hi; ++it) {
273 if (auto e = it->lock()) {
274 auto side_ptr = e->getSideNumberPtr();
275 if (const auto side_number = side_ptr->side_number; 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
598
599 // for (int s = H1; s != LASTSPACE; ++s) {
600 // dataOnElement[s]->dataOnEntities[MBENTITYSET].resize(0);
601 // }
602
603 for (auto &data : dataOnElement) {
604 if (data) {
605 for (auto &dat : data->dataOnEntities) {
606 for (auto &ent_dat : dat) {
607 ent_dat.getEntDataBitRefLevel().reset();
608 }
609 }
610 }
611 }
612
613 auto &field_ents = getDataFieldEnts();
614 for (auto it : field_ents) {
615 if (auto e = it.lock()) {
616 const FieldSpace space = e->getSpace();
617 if (space > NOFIELD) {
618 const EntityType type = e->getEntType();
619 const signed char side =
620 type == MBVERTEX ? 0 : e->getSideNumberPtr()->side_number;
621 if (side >= 0) {
622 if (auto &data = dataOnElement[space]) {
623 auto &dat = data->dataOnEntities[type][side];
624 dat.getEntDataBitRefLevel() |= e->getBitRefLevel();
625 }
626 }
627 } else {
628 if (auto &data = dataOnElement[NOFIELD]) {
629 auto &dat = data->dataOnEntities[MBENTITYSET][0];
630 dat.getEntDataBitRefLevel() |= e->getBitRefLevel();
631 }
632 }
633 }
634 }
635
637};
638
641 const int bit_number) const {
642
643 auto get_nodes_field_data = [&](VectorDouble &nodes_data,
644 VectorFieldEntities &field_entities,
645 VectorDofs &nodes_dofs, FieldSpace &space,
647 VectorInt &bb_node_order) {
649
650 nodes_data.resize(0, false);
651 nodes_dofs.resize(0, false);
652 field_entities.resize(0, false);
653
654 auto field_it =
655 fieldsPtr->get<BitFieldId_mi_tag>().find(BitFEId().set(bit_number - 1));
656 if (field_it != fieldsPtr->get<BitFieldId_mi_tag>().end()) {
657
658#ifndef NDEBUG
659 if ((*field_it)->getBitNumber() != bit_number)
660 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong bit number");
661#endif
662 const int nb_dofs_on_vert = (*field_it)->getNbOfCoeffs();
663 space = (*field_it)->getSpace();
664 base = (*field_it)->getApproxBase();
665
666 auto &field_ents = getDataFieldEnts();
667 const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
668 bit_number, get_id_for_min_type<MBVERTEX>());
669 auto lo = std::lower_bound(field_ents.begin(), field_ents.end(), lo_uid,
670 cmp_uid_lo);
671 if (lo != field_ents.end()) {
672 const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
673 bit_number, get_id_for_max_type<MBVERTEX>());
674 auto hi = std::upper_bound(lo, field_ents.end(), hi_uid, cmp_uid_hi);
675 if (lo != hi) {
676
677 int nb_dofs = 0;
678 for (auto it = lo; it != hi; ++it) {
679 if (auto e = it->lock()) {
680 nb_dofs += e->getNbDofsOnEnt();
681 }
682 }
683
684 if (nb_dofs) {
685
686 const int num_nodes = getNumberOfNodes();
687 bb_node_order.resize(num_nodes, false);
688 bb_node_order.clear();
689 const int max_nb_dofs = nb_dofs_on_vert * num_nodes;
690 nodes_data.resize(max_nb_dofs, false);
691 nodes_dofs.resize(max_nb_dofs, false);
692 field_entities.resize(num_nodes, false);
693 std::fill(nodes_data.begin(), nodes_data.end(), 0);
694 std::fill(nodes_dofs.begin(), nodes_dofs.end(), nullptr);
695 std::fill(field_entities.begin(), field_entities.end(), nullptr);
696
697 std::vector<boost::weak_ptr<FieldEntity>> brother_ents_vec;
698
699 for (auto it = lo; it != hi; ++it) {
700 if (auto e = it->lock()) {
701 const auto &sn = e->getSideNumberPtr();
702 // Some field entities on skeleton can have negative side
703 // number
704 if (const auto side_number = sn->side_number; 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<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 = problem_ptr->numeredFiniteElementsPtr
1685
1686 int nn = 0;
1687 side_fe->loopSize = adjacent_ents.size();
1688 for (auto fe_ent : adjacent_ents) {
1689 auto miit = numered_fe.find(boost::make_tuple(fe_name, fe_ent));
1690 if (miit != numered_fe.end()) {
1691
1692 if (verb >= VERBOSE)
1693 MOFEM_LOG("SELF", sev) << "Side finite element "
1694 << "(" << nn << "): " << **miit;
1695
1696 side_fe->nInTheLoop = nn++;
1697 side_fe->numeredEntFiniteElementPtr = *miit;
1698 CHKERR (*side_fe)();
1699 }
1700 }
1701
1702 CHKERR side_fe->postProcess();
1703
1705}
1706
1707MoFEMErrorCode ForcesAndSourcesCore::UserDataOperator::loopThis(
1708 const string &fe_name, ForcesAndSourcesCore *this_fe, const int verb,
1709 const LogManager::SeverityLevel sev) {
1711
1712 if (verb >= VERBOSE)
1713 MOFEM_LOG("SELF", sev) << "This finite element: "
1714 << *getNumeredEntFiniteElementPtr();
1715
1716 const auto *problem_ptr = getFEMethod()->problemPtr;
1717 this_fe->feName = fe_name;
1718
1719 CHKERR this_fe->setRefineFEPtr(ptrFE);
1720 CHKERR this_fe->copyBasicMethod(*getFEMethod());
1721 CHKERR this_fe->copyPetscData(*getFEMethod());
1722 CHKERR this_fe->copyKsp(*getFEMethod());
1723 CHKERR this_fe->copySnes(*getFEMethod());
1724 CHKERR this_fe->copyTs(*getFEMethod());
1725
1726 CHKERR this_fe->preProcess();
1727
1728 this_fe->nInTheLoop = getNinTheLoop();
1729 this_fe->numeredEntFiniteElementPtr = getNumeredEntFiniteElementPtr();
1730
1731 CHKERR (*this_fe)();
1732
1733 CHKERR this_fe->postProcess();
1734
1736}
1737
1738MoFEMErrorCode ForcesAndSourcesCore::UserDataOperator::loopParent(
1739 const string &fe_name, ForcesAndSourcesCore *parent_fe, const int verb,
1740 const LogManager::SeverityLevel sev) {
1742
1743 const auto *problem_ptr = getFEMethod()->problemPtr;
1744 auto &numered_fe = problem_ptr->numeredFiniteElementsPtr
1746
1747 const auto parent_ent = getNumeredEntFiniteElementPtr()->getParentEnt();
1748 auto miit = numered_fe.find(boost::make_tuple(fe_name, parent_ent));
1749 if (miit != numered_fe.end()) {
1750
1751 if (verb >= VERBOSE)
1752 MOFEM_LOG("SELF", sev) << "Parent finite element: " << **miit;
1753
1754 parent_fe->feName = fe_name;
1755
1756 CHKERR parent_fe->setRefineFEPtr(ptrFE);
1757 CHKERR parent_fe->copyBasicMethod(*getFEMethod());
1758 CHKERR parent_fe->copyPetscData(*getFEMethod());
1759 CHKERR parent_fe->copyKsp(*getFEMethod());
1760 CHKERR parent_fe->copySnes(*getFEMethod());
1761 CHKERR parent_fe->copyTs(*getFEMethod());
1762
1763 CHKERR parent_fe->preProcess();
1764
1765 parent_fe->loopSize = 1;
1766 parent_fe->nInTheLoop = 0;
1767 parent_fe->numeredEntFiniteElementPtr = *miit;
1768 CHKERR (*parent_fe)();
1769
1770 CHKERR parent_fe->postProcess();
1771 }
1772
1774}
1775
1776MoFEMErrorCode ForcesAndSourcesCore::UserDataOperator::loopChildren(
1777 const string &fe_name, ForcesAndSourcesCore *child_fe, const int verb,
1778 const LogManager::SeverityLevel sev) {
1780
1781 const auto *problem_ptr = getFEMethod()->problemPtr;
1782 auto &ref_ents = *getPtrFE()->mField.get_ref_ents();
1783 auto &numered_fe = problem_ptr->numeredFiniteElementsPtr
1785
1786 const auto parent_ent = getNumeredEntFiniteElementPtr()->getEnt();
1787 const auto parent_type = getNumeredEntFiniteElementPtr()->getEntType();
1788 auto range =
1789 ref_ents.get<Composite_ParentEnt_And_EntType_mi_tag>().equal_range(
1790 boost::make_tuple(parent_type, parent_ent));
1791
1792 if (auto size = std::distance(range.first, range.second)) {
1793
1794 std::vector<EntityHandle> childs_vec;
1795 childs_vec.reserve(size);
1796 for (; range.first != range.second; ++range.first)
1797 childs_vec.emplace_back((*range.first)->getEnt());
1798
1799 Range childs;
1800
1801 if ((childs_vec.back() - childs_vec.front() + 1) == size)
1802 childs = Range(childs_vec.front(), childs_vec.back());
1803 else
1804 childs.insert_list(childs_vec.begin(), childs_vec.end());
1805
1806 child_fe->feName = fe_name;
1807
1808 CHKERR child_fe->setRefineFEPtr(ptrFE);
1809 CHKERR child_fe->copyBasicMethod(*getFEMethod());
1810 CHKERR child_fe->copyPetscData(*getFEMethod());
1811 CHKERR child_fe->copyKsp(*getFEMethod());
1812 CHKERR child_fe->copySnes(*getFEMethod());
1813 CHKERR child_fe->copyTs(*getFEMethod());
1814
1815 CHKERR child_fe->preProcess();
1816
1817 int nn = 0;
1818 child_fe->loopSize = size;
1819
1820 for (auto p = childs.pair_begin(); p != childs.pair_end(); ++p) {
1821
1822 auto miit = numered_fe.lower_bound(boost::make_tuple(fe_name, p->first));
1823 auto hi_miit =
1824 numered_fe.upper_bound(boost::make_tuple(fe_name, p->second));
1825
1826 for (; miit != hi_miit; ++miit) {
1827
1828 if (verb >= VERBOSE)
1829 MOFEM_LOG("SELF", sev) << "Child finite element "
1830 << "(" << nn << "): " << **miit;
1831
1832 child_fe->nInTheLoop = nn++;
1833 child_fe->numeredEntFiniteElementPtr = *miit;
1834 CHKERR (*child_fe)();
1835 }
1836 }
1837
1838 CHKERR child_fe->postProcess();
1839 };
1840
1842}
1843
1844int ForcesAndSourcesCore::getRule(int order_row, int order_col,
1845 int order_data) {
1846 return getRuleHook ? getRuleHook(order_row, order_col, order_data)
1847 : getRule(order_data);
1848}
1849
1851 int order_data) {
1852 return setRuleHook ? setRuleHook(this, order_row, order_col, order_data)
1853 : setGaussPts(order_data);
1854}
1855
1857
1858/** \deprecated setGaussPts(int row_order, int col_order, int data order);
1859 */
1862 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Sorry, not implemented");
1864}
1865
1866ForcesAndSourcesCore::UserDataOperator::UserDataOperator(const FieldSpace space,
1867 const char type,
1868 const bool symm)
1869 : DataOperator(symm), opType(type), sPace(space), ptrFE(nullptr) {}
1870
1871ForcesAndSourcesCore::UserDataOperator::UserDataOperator(
1872 const std::string field_name, const char type, const bool symm)
1873 : DataOperator(symm), opType(type), rowFieldName(field_name),
1874 colFieldName(field_name), sPace(LASTSPACE), ptrFE(nullptr) {}
1875
1876ForcesAndSourcesCore::UserDataOperator::UserDataOperator(
1877 const std::string row_field_name, const std::string col_field_name,
1878 const char type, const bool symm)
1879 : DataOperator(symm), opType(type), rowFieldName(row_field_name),
1880 colFieldName(col_field_name), sPace(LASTSPACE), ptrFE(nullptr) {}
1881
1884 if (preProcessHook) {
1885 ierr = preProcessHook();
1886 CHKERRG(ierr);
1887 }
1889}
1892 if (operatorHook) {
1893 ierr = operatorHook();
1894 CHKERRG(ierr);
1895 } else {
1896#ifndef NDEBUG
1897 MOFEM_LOG("SELF", Sev::warning)
1898 << "No method operator() overloaded on element entity on finite "
1899 "element <"
1900 << boost::typeindex::type_id_runtime(*this).pretty_name() << ">";
1901#endif
1902 }
1904}
1907 if (postProcessHook) {
1909 CHKERRG(ierr);
1910 }
1912}
1913
1915ForcesAndSourcesCore::UserDataOperator::setPtrFE(ForcesAndSourcesCore *ptr) {
1917 if (!(ptrFE = dynamic_cast<ForcesAndSourcesCore *>(ptr)))
1918 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1919 "User operator and finite element do not work together");
1921}
1922
1923} // 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 Field * get_field_structure(const std::string &name)=0
get field structure
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
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
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
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.
boost::ptr_vector< UserDataOperator > opPtrVector
Vector of finite element users data operators.
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