v0.15.0
Loading...
Searching...
No Matches
FECore.cpp
Go to the documentation of this file.
1/** \file FECore.cpp
2 * \brief Core interface methods for managing deletions and insertion dofs
3 */
4
5#include <MoFEM.hpp>
6
7#define FECoreFunctionBegin \
8 MoFEMFunctionBegin; \
9 MOFEM_LOG_CHANNEL("WORLD"); \
10 MOFEM_LOG_CHANNEL("SYNC"); \
11 MOFEM_LOG_FUNCTION(); \
12 MOFEM_LOG_TAG("WORLD", "FECore"); \
13 MOFEM_LOG_TAG("SYNC", "FECore");
14
15namespace MoFEM {
16
17const FiniteElement *
18Core::get_finite_element_structure(const std::string &name,
19 enum MoFEMTypes bh) const {
20 auto miit = finiteElements.get<FiniteElement_name_mi_tag>().find(name);
21 if (miit == finiteElements.get<FiniteElement_name_mi_tag>().end()) {
22 if (bh == MF_EXIST) {
23 throw MoFEMException(
25 std::string("finite element < " + name +
26 " > not in database (top tip: check spelling)")
27 .c_str());
28 } else {
29 return nullptr;
30 }
31 }
32 return miit->get();
33}
34
35bool Core::check_finite_element(const std::string &name) const {
36 auto miit = finiteElements.get<FiniteElement_name_mi_tag>().find(name);
37 if (miit == finiteElements.get<FiniteElement_name_mi_tag>().end())
38 return false;
39 else
40 return true;
41}
42
43MoFEMErrorCode Core::add_finite_element(const std::string &fe_name,
44 enum MoFEMTypes bh, int verb) {
46 *buildMoFEM &= 1 << 0;
47 if (verb == -1) {
48 verb = verbose;
49 }
50
51 // Add finite element meshset to partition meshset. In case of no elements
52 // on processor part, when mesh file is read, finite element meshset is
53 // prevented from deletion by moab reader.
54 auto add_meshset_to_partition = [&](auto meshset) {
56 const void *tag_vals[] = {&rAnk};
57 ParallelComm *pcomm = ParallelComm::get_pcomm(
58 &get_moab(), get_basic_entity_data_ptr()->pcommID);
59 Tag part_tag = pcomm->part_tag();
60 Range tagged_sets;
61 CHKERR get_moab().get_entities_by_type_and_tag(0, MBENTITYSET, &part_tag,
62 tag_vals, 1, tagged_sets,
63 moab::Interface::UNION);
64 for (auto s : tagged_sets)
65 CHKERR get_moab().add_entities(s, &meshset, 1);
67 };
68
69 auto &finite_element_name_set =
70 finiteElements.get<FiniteElement_name_mi_tag>();
71 auto it_fe = finite_element_name_set.find(fe_name);
72
73 if (bh == MF_EXCL) {
74 if (it_fe != finite_element_name_set.end()) {
75 SETERRQ(mofemComm, MOFEM_NOT_FOUND, "this < %s > is there",
76 fe_name.c_str());
77 }
78
79 } else {
80 if (it_fe != finite_element_name_set.end())
82 }
83 EntityHandle meshset;
84 CHKERR get_moab().create_meshset(MESHSET_SET, meshset);
85 CHKERR add_meshset_to_partition(meshset);
86
87 // id
88 int fe_shift = 0;
89 for (; finiteElements.get<BitFEId_mi_tag>().find(BitFEId().set(fe_shift)) !=
90 finiteElements.get<BitFEId_mi_tag>().end();
91 ++fe_shift) {
92 }
93
94 auto id = BitFEId().set(fe_shift);
95 CHKERR get_moab().tag_set_data(th_FEId, &meshset, 1, &id);
96
97 // id name
98 void const *tag_data[] = {fe_name.c_str()};
99 int tag_sizes[1];
100 tag_sizes[0] = fe_name.size();
101 CHKERR get_moab().tag_set_by_ptr(th_FEName, &meshset, 1, tag_data, tag_sizes);
102
103 // add FiniteElement
104 auto p = finiteElements.insert(
105 boost::shared_ptr<FiniteElement>(new FiniteElement(moab, meshset)));
106 if (!p.second)
107 SETERRQ(mofemComm, MOFEM_OPERATION_UNSUCCESSFUL,
108 "FiniteElement not inserted");
109
110 if (verb > QUIET)
111 MOFEM_LOG("WORLD", Sev::inform) << "Add finite element " << fe_name;
112
114}
115
118 const EntityType type,
119 ElementAdjacencyFunct function) {
121 *buildMoFEM &= 1 << 0;
122 typedef FiniteElement_multiIndex::index<FiniteElement_name_mi_tag>::type
123 FiniteElements_by_name;
124 FiniteElements_by_name &finite_element_name_set =
125 finiteElements.get<FiniteElement_name_mi_tag>();
126 FiniteElements_by_name::iterator it_fe =
127 finite_element_name_set.find(fe_name);
128 if (it_fe == finite_element_name_set.end())
129 SETERRQ(mofemComm, MOFEM_NOT_FOUND,
130 "This finite element is not defined (advise: check spelling)");
131 boost::shared_ptr<FiniteElement> fe;
132 fe = *it_fe;
133 fe->elementAdjacencyTable[type] = function;
135}
136
139 const std::string name_data) {
141 *buildMoFEM &= 1 << 0;
142 typedef FiniteElement_multiIndex::index<FiniteElement_name_mi_tag>::type
143 FiniteElements_by_name;
144 FiniteElements_by_name &finite_element_name_set =
145 finiteElements.get<FiniteElement_name_mi_tag>();
146 FiniteElements_by_name::iterator it_fe =
147 finite_element_name_set.find(fe_name);
148 if (it_fe == finite_element_name_set.end())
149 SETERRQ(mofemComm, MOFEM_NOT_FOUND,
150 "This finite element is not defined (advise: check spelling)");
151 bool success = finite_element_name_set.modify(
152 it_fe, FiniteElement_change_bit_add(get_field_id(name_data)));
153 if (!success)
154 SETERRQ(mofemComm, MOFEM_OPERATION_UNSUCCESSFUL,
155 "modification unsuccessful");
157}
158
161 const std::string name_row) {
163 *buildMoFEM &= 1 << 0;
164 typedef FiniteElement_multiIndex::index<FiniteElement_name_mi_tag>::type
165 FiniteElements_by_name;
166 FiniteElements_by_name &finite_element_name_set =
167 finiteElements.get<FiniteElement_name_mi_tag>();
168 FiniteElements_by_name::iterator it_fe =
169 finite_element_name_set.find(fe_name);
170 if (it_fe == finite_element_name_set.end())
171 SETERRQ(mofemComm, MOFEM_NOT_FOUND, "this < %s > is not there",
172 fe_name.c_str());
173 bool success = finite_element_name_set.modify(
174 it_fe, FiniteElement_row_change_bit_add(get_field_id(name_row)));
175 if (!success)
176 SETERRQ(mofemComm, MOFEM_OPERATION_UNSUCCESSFUL,
177 "modification unsuccessful");
179}
180
183 const std::string name_col) {
185 *buildMoFEM &= 1 << 0;
186 auto &finite_element_name_set =
187 finiteElements.get<FiniteElement_name_mi_tag>();
188 auto it_fe = finite_element_name_set.find(fe_name);
189 if (it_fe == finite_element_name_set.end())
190 SETERRQ(mofemComm, MOFEM_OPERATION_UNSUCCESSFUL,
191 "this FiniteElement is there");
192 bool success = finite_element_name_set.modify(
193 it_fe, FiniteElement_col_change_bit_add(get_field_id(name_col)));
194 if (!success)
195 SETERRQ(mofemComm, MOFEM_OPERATION_UNSUCCESSFUL,
196 "modification unsuccessful");
198}
199
202 const std::string name_data) {
204 *buildMoFEM &= 1 << 0;
205 auto &finite_element_name_set =
206 finiteElements.get<FiniteElement_name_mi_tag>();
207 auto it_fe = finite_element_name_set.find(fe_name);
208 if (it_fe == finite_element_name_set.end())
209 SETERRQ(mofemComm, MOFEM_NOT_FOUND, "this FiniteElement is there");
210 bool success = finite_element_name_set.modify(
211 it_fe, FiniteElement_change_bit_off(get_field_id(name_data)));
212 if (!success)
213 SETERRQ(mofemComm, MOFEM_OPERATION_UNSUCCESSFUL,
214 "modification unsuccessful");
216}
217
220 const std::string name_row) {
222 *buildMoFEM &= 1 << 0;
223 auto &finite_element_name_set =
224 finiteElements.get<FiniteElement_name_mi_tag>();
225 auto it_fe = finite_element_name_set.find(fe_name);
226 if (it_fe == finite_element_name_set.end())
227 SETERRQ(mofemComm, MOFEM_NOT_FOUND, "this < %s > is not there",
228 fe_name.c_str());
229 bool success = finite_element_name_set.modify(
230 it_fe, FiniteElement_row_change_bit_off(get_field_id(name_row)));
231 if (!success)
232 SETERRQ(mofemComm, MOFEM_OPERATION_UNSUCCESSFUL,
233 "modification unsuccessful");
235}
236
239 const std::string name_col) {
241 *buildMoFEM &= 1 << 0;
242 auto &finite_element_name_set =
243 finiteElements.get<FiniteElement_name_mi_tag>();
244 auto it_fe = finite_element_name_set.find(fe_name);
245 if (it_fe == finite_element_name_set.end())
246 SETERRQ(mofemComm, MOFEM_NOT_FOUND, "this FiniteElement is there");
247 bool success = finite_element_name_set.modify(
248 it_fe, FiniteElement_col_change_bit_off(get_field_id(name_col)));
249 if (!success)
250 SETERRQ(mofemComm, MOFEM_OPERATION_UNSUCCESSFUL,
251 "modification unsuccessful");
253}
254
255BitFEId Core::getBitFEId(const std::string &fe_name) const {
256 auto &fe_by_id = finiteElements.get<FiniteElement_name_mi_tag>();
257 auto miit = fe_by_id.find(fe_name);
258 if (miit == fe_by_id.end())
260 ("finite element < " + fe_name + " > not found (top tip: check spelling)")
261 .c_str());
262 return (*miit)->getId();
263}
264
265std::string Core::getBitFEIdName(const BitFEId id) const {
266 auto &fe_by_id = finiteElements.get<BitFEId_mi_tag>();
267 auto miit = fe_by_id.find(id);
268 if (miit == fe_by_id.end())
269 THROW_MESSAGE("finite element not found");
270 return (*miit)->getName();
271}
272
274 auto &fe_by_id = finiteElements.get<BitFEId_mi_tag>();
275 auto miit = fe_by_id.find(id);
276 if (miit == fe_by_id.end())
277 THROW_MESSAGE("finite element not found");
278 return (*miit)->meshset;
279}
280
281EntityHandle Core::get_finite_element_meshset(const std::string name) const {
282 return get_finite_element_meshset(getBitFEId(name));
283}
284
287 Range &ents) const {
288
290
291 EntityHandle meshset = get_finite_element_meshset(name);
292 CHKERR get_moab().get_entities_by_dimension(meshset, dim, ents, true);
294}
295
297 EntityType type,
298 Range &ents) const {
299
301
302 EntityHandle meshset = get_finite_element_meshset(name);
303 CHKERR get_moab().get_entities_by_type(meshset, type, ents, true);
304
306}
307
310 Range &ents) const {
311
313
314 EntityHandle meshset = get_finite_element_meshset(name);
315 CHKERR get_moab().get_entities_by_handle(meshset, ents, true);
316
318}
319
322 for (auto &fe : finiteElements.get<FiniteElement_name_mi_tag>())
323 MOFEM_LOG("SYNC", Sev::inform) << fe;
324
325 MOFEM_LOG_SYNCHRONISE(mofemComm);
327}
328
330 const EntityHandle meshset, const EntityType type, const std::string name,
331 const bool recursive) {
332 *buildMoFEM &= 1 << 0;
335
336 idm = get_finite_element_meshset(getBitFEId(name));
337 Range ents;
338 CHKERR get_moab().get_entities_by_type(meshset, type, ents, recursive);
339 CHKERR getInterface<BitRefManager>()->setElementsBitRefLevel(ents);
340 CHKERR get_moab().add_entities(idm, ents);
341
343}
344
347 const int dim, const std::string name,
348 const bool recursive) {
350 *buildMoFEM &= 1 << 0;
352 idm = get_finite_element_meshset(getBitFEId(name));
353 Range ents;
354 CHKERR get_moab().get_entities_by_dimension(meshset, dim, ents, recursive);
355 CHKERR getInterface<BitRefManager>()->setElementsBitRefLevel(ents);
356 CHKERR get_moab().add_entities(idm, ents);
358}
359
361 const Range ents, const EntityType type, const std::string name) {
363 *buildMoFEM &= 1 << 0;
365 idm = get_finite_element_meshset(getBitFEId(name));
366 CHKERR getInterface<BitRefManager>()->setElementsBitRefLevel(
367 ents.subset_by_type(type));
368 CHKERR get_moab().add_entities(idm, ents.subset_by_type(type));
370} // namespace MoFEM
371
374 const std::string name) {
376 *buildMoFEM &= 1 << 0;
378 idm = get_finite_element_meshset(getBitFEId(name));
379 CHKERR getInterface<BitRefManager>()->setElementsBitRefLevel(
380 ents.subset_by_dimension(dim));
381 CHKERR get_moab().add_entities(idm, ents.subset_by_dimension(dim));
383}
384
387 const std::string &name,
388 EntityType type, int verb) {
390 CHKERR add_ents_to_finite_element_by_bit_ref(bit, BitRefLevel().set(), name,
391 type, verb);
392
394}
395
397 const BitRefLevel &bit, const BitRefLevel &mask, const std::string &name,
398 EntityType type, int verb) {
400 CHKERR add_ents_to_finite_element_by_bit_ref(bit, mask, name, type, verb);
401
403}
404
406 const BitRefLevel bit, const BitRefLevel mask, const std::string name,
407 EntityType type, int verb) {
409
410 if (verb == -1)
411 verb = verbose;
412 *buildMoFEM &= 1 << 0;
413 const BitFEId id = getBitFEId(name);
414 const EntityHandle idm = get_finite_element_meshset(id);
415
416 auto &ref_MoFEMFiniteElement = refinedFiniteElements.get<Ent_mi_tag>();
417 auto miit = ref_MoFEMFiniteElement.lower_bound(get_id_for_min_type(type));
418 auto hi_miit = ref_MoFEMFiniteElement.upper_bound(get_id_for_max_type(type));
419
420 int nb_add_fes = 0;
421 for (; miit != hi_miit; miit++) {
422 const auto &bit2 = miit->get()->getBitRefLevel();
423 if ((bit2 & mask) != bit2)
424 continue;
425 if ((bit2 & bit).any()) {
426 EntityHandle ent = miit->get()->getEnt();
427 CHKERR get_moab().add_entities(idm, &ent, 1);
428 nb_add_fes++;
429 }
430 }
431
432 MOFEM_LOG("SYNC", Sev::inform)
433 << "Finite element " << name << " added. Nb. of elements added "
434 << nb_add_fes << " out of " << std::distance(miit, hi_miit);
435
436 MOFEM_LOG_SYNCHRONISE(mofemComm)
437
439}
440
442 const EntityHandle meshset, const std::string &name, const bool recursive) {
444 *buildMoFEM &= 1 << 0;
445 const BitFEId id = getBitFEId(name);
446 const EntityHandle idm = get_finite_element_meshset(id);
447 if (recursive == false) {
448 CHKERR get_moab().add_entities(idm, &meshset, 1);
449 } else {
450 Range meshsets;
451 CHKERR get_moab().get_entities_by_type(meshset, MBENTITYSET, meshsets,
452 false);
453 CHKERR get_moab().add_entities(idm, meshsets);
454 }
456}
457
459Core::buildFiniteElements(const boost::shared_ptr<FiniteElement> &fe,
460 const Range *ents_ptr, int verb) {
462 if (verb == DEFAULT_VERBOSITY)
463 verb = verbose;
464
465 if (verb > QUIET)
466 MOFEM_LOG("SYNC", Sev::verbose)
467 << "Build Finite Elements " << fe->getName();
468
469 auto &fields_by_id = fIelds.get<BitFieldId_mi_tag>();
470
471 // Get id of mofem fields for row, col and data
472 enum IntLoop { ROW = 0, COL, DATA, LAST };
473 std::array<BitFieldId, LAST> fe_fields = {fe.get()->getBitFieldIdRow(),
474 fe.get()->getBitFieldIdCol(),
475 fe.get()->getBitFieldIdData()};
476
477 // Get finite element meshset
478 EntityHandle meshset = get_finite_element_meshset(fe.get()->getId());
479
480 // Get entities from finite element meshset // if meshset
481 Range fe_ents;
482 CHKERR get_moab().get_entities_by_handle(meshset, fe_ents, false);
483
484 if (ents_ptr)
485 fe_ents = intersect(fe_ents, *ents_ptr);
486
487 // Map entity uid to pointers
488 typedef std::vector<boost::weak_ptr<EntFiniteElement>> VecOfWeakFEPtrs;
489 VecOfWeakFEPtrs processed_fes;
490 processed_fes.reserve(fe_ents.size());
491
492 int last_data_field_ents_view_size = 0;
493 int last_row_field_ents_view_size = 0;
494 int last_col_field_ents_view_size = 0;
495
496 // Entities adjacent to entities
497 std::vector<EntityHandle> adj_ents;
498
499 // Loop meshset finite element ents and add finite elements
500 for (Range::const_pair_iterator peit = fe_ents.const_pair_begin();
501 peit != fe_ents.const_pair_end(); peit++) {
502
503 const auto first = peit->first;
504 const auto second = peit->second;
505
506 // Find range of ref entities that is sequence
507 // note: iterator is a wrapper
508 // check if is in refinedFiniteElements database
509 auto ref_fe_miit =
510 refinedFiniteElements.get<Ent_mi_tag>().lower_bound(first);
511 auto hi_ref_fe_miit =
512 refinedFiniteElements.get<Ent_mi_tag>().upper_bound(second);
513 if (std::distance(ref_fe_miit, hi_ref_fe_miit) != (second - first + 1)) {
514 MOFEM_LOG("SELF", Sev::noisy) << "Finite element " << fe->getName()
515 << " not defined on all entities in "
516 "the meshset, missing entities";
517 MOFEM_LOG("SELF", Sev::noisy)
518 << "Missing entities: " << first << " - " << second;
519 }
520
521 EntFiniteElement_multiIndex::iterator hint_p = entsFiniteElements.end();
522 for (; ref_fe_miit != hi_ref_fe_miit; ref_fe_miit++) {
523
524 // Add finite element to database
525 hint_p = entsFiniteElements.emplace_hint(
526 hint_p, boost::make_shared<EntFiniteElement>(*ref_fe_miit, fe));
527 processed_fes.emplace_back(*hint_p);
528 auto fe_raw_ptr = hint_p->get();
529
530 // Allocate space for entities view
531 bool row_as_data = false, col_as_row = false;
532 if (fe_fields[DATA] == fe_fields[ROW])
533 row_as_data = true;
534 if (fe_fields[ROW] == fe_fields[COL])
535 col_as_row = true;
536
537 fe_raw_ptr->getDataFieldEntsPtr()->reserve(
538 last_data_field_ents_view_size);
539
540 if (row_as_data) {
541 fe_raw_ptr->getRowFieldEntsPtr() = fe_raw_ptr->getDataFieldEntsPtr();
542 } else {
543 // row and col are different
544 if (fe_raw_ptr->getRowFieldEntsPtr() ==
545 fe_raw_ptr->getDataFieldEntsPtr())
546 fe_raw_ptr->getRowFieldEntsPtr() =
547 boost::make_shared<FieldEntity_vector_view>();
548 fe_raw_ptr->getRowFieldEntsPtr()->reserve(
549 last_row_field_ents_view_size);
550 }
551
552 if (row_as_data && col_as_row) {
553 fe_raw_ptr->getColFieldEntsPtr() = fe_raw_ptr->getDataFieldEntsPtr();
554 } else if (col_as_row) {
555 fe_raw_ptr->getColFieldEntsPtr() = fe_raw_ptr->getRowFieldEntsPtr();
556 } else {
557 if (
558
559 fe_raw_ptr->getColFieldEntsPtr() ==
560 fe_raw_ptr->getRowFieldEntsPtr() ||
561 fe_raw_ptr->getColFieldEntsPtr() ==
562 fe_raw_ptr->getDataFieldEntsPtr()
563
564 )
565 fe_raw_ptr->getColFieldEntsPtr() =
566 boost::make_shared<FieldEntity_vector_view>();
567 fe_raw_ptr->getColFieldEntsPtr()->reserve(
568 last_col_field_ents_view_size);
569 }
570
571 // Iterate over all field and check which one is on the element
572 for (unsigned int ii = 0; ii != BitFieldId().size(); ++ii) {
573
574 // Common field id for ROW, COL and DATA
575 BitFieldId id_common = 0;
576 // Check if the field (ii) is added to finite element
577 for (int ss = 0; ss < LAST; ss++) {
578 id_common |= fe_fields[ss] & BitFieldId().set(ii);
579 }
580 if (id_common.none())
581 continue;
582
583 // Find in database data associated with the field (ii)
584 const BitFieldId field_id = BitFieldId().set(ii);
585 auto miit = fields_by_id.find(field_id);
586 if (miit == fields_by_id.end())
587 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Field not found");
588 auto field_bit_number = (*miit)->getBitNumber();
589
590 // Loop over adjacencies of element and find field entities on those
591 // adjacencies, that create hash map map_uid_fe which is used later
592 const std::string field_name = miit->get()->getName();
593 const bool add_to_data = (field_id & fe_fields[DATA]).any();
594 const bool add_to_row = (field_id & fe_fields[ROW]).any();
595 const bool add_to_col = (field_id & fe_fields[COL]).any();
596
597 // Resolve entities on element, those entities are used to build tag
598 // with dof uids on finite element tag
599 adj_ents.clear();
600 CHKERR fe_raw_ptr->getElementAdjacency(*miit, adj_ents);
601
602 for (auto ent : adj_ents) {
603
604 auto dof_it = entsFields.get<Unique_mi_tag>().find(
605 FieldEntity::getLocalUniqueIdCalculate(field_bit_number, ent));
606 if (dof_it != entsFields.get<Unique_mi_tag>().end()) {
607
608 if (add_to_data) {
609 fe_raw_ptr->getDataFieldEntsPtr()->emplace_back(*dof_it);
610 }
611 if (add_to_row && !row_as_data) {
612 fe_raw_ptr->getRowFieldEntsPtr()->emplace_back(*dof_it);
613 }
614 if (add_to_col && !col_as_row) {
615 fe_raw_ptr->getColFieldEntsPtr()->emplace_back(*dof_it);
616 }
617
618 }
619 }
620 }
621
622 // Sort field ents by uid
623 auto uid_comp = [](const auto &a, const auto &b) {
624 return a.lock()->getLocalUniqueId() < b.lock()->getLocalUniqueId();
625 };
626
627 // Sort all views
628
629 // Data
630 sort(fe_raw_ptr->getDataFieldEntsPtr()->begin(),
631 fe_raw_ptr->getDataFieldEntsPtr()->end(), uid_comp);
632 last_data_field_ents_view_size =
633 fe_raw_ptr->getDataFieldEntsPtr()->size();
634
635 // Row
636 if (!row_as_data) {
637 sort(fe_raw_ptr->getRowFieldEntsPtr()->begin(),
638 fe_raw_ptr->getRowFieldEntsPtr()->end(), uid_comp);
639 last_row_field_ents_view_size =
640 fe_raw_ptr->getRowFieldEntsPtr()->size();
641 }
642
643 // Column
644 if (!col_as_row) {
645 sort(fe_raw_ptr->getColFieldEntsPtr()->begin(),
646 fe_raw_ptr->getColFieldEntsPtr()->end(), uid_comp);
647 last_col_field_ents_view_size =
648 fe_raw_ptr->getColFieldEntsPtr()->size();
649 }
650 }
651 }
652
654}
655
658
659 if (verb == DEFAULT_VERBOSITY)
660 verb = verbose;
661
662 // loop Finite Elements
663 for (auto &fe : finiteElements)
664 CHKERR buildFiniteElements(fe, NULL, verb);
665
666 if (verb > QUIET) {
667
668 auto &fe_ents = entsFiniteElements.get<Unique_mi_tag>();
669 for (auto &fe : finiteElements) {
670 auto miit = fe_ents.lower_bound(
672 auto hi_miit =
674 get_id_for_max_type<MBENTITYSET>(), fe->getFEUId()));
675 const auto count = std::distance(miit, hi_miit);
676 MOFEM_LOG("SYNC", Sev::inform)
677 << "Finite element " << fe->getName()
678 << " added. Nb. of elements added " << count;
679 MOFEM_LOG("SYNC", Sev::noisy) << *fe;
680
681 auto slg = MoFEM::LogManager::getLog("SYNC");
682 for (auto &field : fIelds) {
683 auto rec = slg.open_record(keywords::severity = Sev::verbose);
684 if (rec) {
685 logging::record_ostream strm(rec);
686 strm << "Field " << field->getName() << " on finite element: ";
687 if ((field->getId() & fe->getBitFieldIdRow()).any())
688 strm << "row ";
689 if ((field->getId() & fe->getBitFieldIdCol()).any())
690 strm << "columns ";
691 if ((field->getId() & fe->getBitFieldIdData()).any())
692 strm << "data";
693 strm.flush();
694 slg.push_record(boost::move(rec));
695 }
696 }
697 }
698
699 MOFEM_LOG_SYNCHRONISE(mofemComm);
700 }
701
702 *buildMoFEM |= 1 << 1;
704}
705
708 SETERRQ(mofemComm, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
710}
711
713 const Range *const ents_ptr,
714 int verb) {
716 if (verb == -1)
717 verb = verbose;
718
719 auto fe_miit = finiteElements.get<FiniteElement_name_mi_tag>().find(fe_name);
720 if (fe_miit == finiteElements.get<FiniteElement_name_mi_tag>().end())
721 SETERRQ(mofemComm, MOFEM_NOT_FOUND, "Finite element <%s> not found",
722 fe_name.c_str());
723
724 CHKERR buildFiniteElements(*fe_miit, ents_ptr, verb);
725
726 if (verb >= VERBOSE) {
727 auto &fe_ents = entsFiniteElements.get<Unique_mi_tag>();
728 auto miit = fe_ents.lower_bound(
729 EntFiniteElement::getLocalUniqueIdCalculate(0, (*fe_miit)->getFEUId()));
730 auto hi_miit =
732 get_id_for_max_type<MBENTITYSET>(), (*fe_miit)->getFEUId()));
733 const auto count = std::distance(miit, hi_miit);
734 MOFEM_LOG("SYNC", Sev::inform) << "Finite element " << fe_name
735 << " added. Nb. of elements added " << count;
736 MOFEM_LOG_SYNCHRONISE(mofemComm);
737 }
738
739 *buildMoFEM |= 1 << 1;
741}
742
745 if (verb == DEFAULT_VERBOSITY)
746 verb = verbose;
747
748 if (!((*buildMoFEM) & BUILD_FIELD))
749 SETERRQ(mofemComm, MOFEM_NOT_FOUND, "field not build");
750 if (!((*buildMoFEM) & BUILD_FE))
751 SETERRQ(mofemComm, MOFEM_NOT_FOUND, "fe not build");
752 for (auto peit = ents.pair_begin(); peit != ents.pair_end(); ++peit) {
753 auto fit = entsFiniteElements.get<Ent_mi_tag>().lower_bound(peit->first);
754 auto hi_fit =
755 entsFiniteElements.get<Ent_mi_tag>().upper_bound(peit->second);
756 for (; fit != hi_fit; ++fit) {
757 if ((*fit)->getBitFieldIdRow().none() &&
758 (*fit)->getBitFieldIdCol().none() &&
759 (*fit)->getBitFieldIdData().none())
760 continue;
761 int by = BYROW;
762 if ((*fit)->getBitFieldIdRow() != (*fit)->getBitFieldIdCol())
763 by |= BYCOL;
764 if ((*fit)->getBitFieldIdRow() != (*fit)->getBitFieldIdData())
765 by |= BYDATA;
767 auto hint = entFEAdjacencies.end();
768 for (auto e : *(*fit)->getRowFieldEntsPtr()) {
769 hint = entFEAdjacencies.emplace_hint(hint, e.lock(), *fit);
770 bool success = entFEAdjacencies.modify(hint, modify_row);
771 if (!success)
772 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
773 "modification unsuccessful");
774 }
775 if ((*fit)->getBitFieldIdRow() != (*fit)->getBitFieldIdCol()) {
776 int by = BYCOL;
777 if ((*fit)->getBitFieldIdCol() != (*fit)->getBitFieldIdData())
778 by |= BYDATA;
780 auto hint = entFEAdjacencies.end();
781 for (auto e : *(*fit)->getColFieldEntsPtr()) {
782 hint = entFEAdjacencies.emplace_hint(hint, e.lock(), *fit);
783 bool success = entFEAdjacencies.modify(hint, modify_col);
784 if (!success)
785 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
786 "modification unsuccessful");
787 }
788 }
789 if ((*fit)->getBitFieldIdRow() != (*fit)->getBitFieldIdData() ||
790 (*fit)->getBitFieldIdCol() != (*fit)->getBitFieldIdData()) {
792 BYDATA);
793 auto hint = entFEAdjacencies.end();
794 for (auto &e : (*fit)->getDataFieldEnts()) {
795 hint = entFEAdjacencies.emplace_hint(hint, e.lock(), *fit);
796 bool success = entFEAdjacencies.modify(hint, modify_data);
797 if (!success)
798 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
799 "modification unsuccessful");
800 }
801 }
802 }
803 }
804
805 if (verb >= VERBOSE) {
806 MOFEM_LOG("WORLD", Sev::inform)
807 << "Number of adjacencies " << entFEAdjacencies.size();
808 MOFEM_LOG_SYNCHRONISE(mofemComm)
809 }
810
811 *buildMoFEM |= 1 << 2;
813}
814
816 const BitRefLevel &mask, int verb) {
818 if (verb == -1)
819 verb = verbose;
820 Range ents;
821 CHKERR BitRefManager(*this).getEntitiesByRefLevel(bit, mask, ents);
822
823 CHKERR build_adjacencies(ents, verb);
824
826}
829 if (verb == -1)
830 verb = verbose;
831 CHKERR build_adjacencies(bit, BitRefLevel().set(), verb);
832
834}
835
836EntFiniteElement_multiIndex::index<Unique_mi_tag>::type::iterator
837Core::get_fe_by_name_begin(const std::string &fe_name) const {
838 auto miit = finiteElements.get<FiniteElement_name_mi_tag>().find(fe_name);
839 if (miit != finiteElements.get<FiniteElement_name_mi_tag>().end()) {
840 return entsFiniteElements.get<Unique_mi_tag>().lower_bound(
841 EntFiniteElement::getLocalUniqueIdCalculate(0, (*miit)->getFEUId()));
842 } else {
843 return entsFiniteElements.get<Unique_mi_tag>().end();
844 }
845}
846
847EntFiniteElement_multiIndex::index<Unique_mi_tag>::type::iterator
848Core::get_fe_by_name_end(const std::string &fe_name) const {
849 auto miit = finiteElements.get<FiniteElement_name_mi_tag>().find(fe_name);
850 if (miit != finiteElements.get<FiniteElement_name_mi_tag>().end()) {
851 return entsFiniteElements.get<Unique_mi_tag>().upper_bound(
853 get_id_for_max_type<MBENTITYSET>(), (*miit)->getFEUId()));
854 } else {
855 return entsFiniteElements.get<Unique_mi_tag>().end();
856 }
857}
858
860 const std::string &name) const {
862 FiniteElement_multiIndex::index<FiniteElement_name_mi_tag>::type::iterator it;
863 it = finiteElements.get<FiniteElement_name_mi_tag>().find(name);
864 if (it == finiteElements.get<FiniteElement_name_mi_tag>().end()) {
865 SETERRQ(mofemComm, 1, "finite element not found < %s >", name.c_str());
866 }
867 EntityHandle meshset = (*it)->getMeshset();
868
869 int num_entities;
870 CHKERR get_moab().get_number_entities_by_handle(meshset, num_entities);
871
872 auto counts_fes = [&]() {
873 return std::distance(get_fe_by_name_begin((*it)->getName()),
874 get_fe_by_name_end((*it)->getName()));
875 };
876
877 if (counts_fes() != static_cast<size_t>(num_entities)) {
878 SETERRQ(mofemComm, MOFEM_DATA_INCONSISTENCY,
879 "not equal number of entities in meshset and finite elements "
880 "multiindex < %s >",
881 (*it)->getName().c_str());
882 }
884}
885
888 FiniteElement_multiIndex::index<FiniteElement_name_mi_tag>::type::iterator it;
889 it = finiteElements.get<FiniteElement_name_mi_tag>().begin();
890 for (; it != finiteElements.get<FiniteElement_name_mi_tag>().end(); it++) {
891 EntityHandle meshset = (*it)->getMeshset();
892
893 int num_entities;
894 CHKERR get_moab().get_number_entities_by_handle(meshset, num_entities);
895
896 auto counts_fes = [&]() {
897 return std::distance(get_fe_by_name_begin((*it)->getName()),
898 get_fe_by_name_end((*it)->getName()));
899 };
900
901 if (counts_fes() != static_cast<size_t>(num_entities)) {
902 SETERRQ(mofemComm, MOFEM_DATA_INCONSISTENCY,
903 "not equal number of entities in meshset and finite elements "
904 "multiindex < %s >",
905 (*it)->getName().c_str());
906 }
907 }
909}
910
912Core::get_problem_finite_elements_entities(const std::string problem_name,
913 const std::string &fe_name,
914 const EntityHandle meshset) {
916 auto &prb = pRoblems.get<Problem_mi_tag>();
917 auto p_miit = prb.find(problem_name);
918 if (p_miit == prb.end())
919 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
920 "No such problem like < %s >", problem_name.c_str());
921
922 auto fe_miit = finiteElements.get<FiniteElement_name_mi_tag>().find(fe_name);
923 if (fe_miit != finiteElements.get<FiniteElement_name_mi_tag>().end()) {
924 auto miit =
925 p_miit->numeredFiniteElementsPtr->get<Unique_mi_tag>().lower_bound(
927 0, (*fe_miit)->getFEUId()));
928 auto hi_miit =
929 p_miit->numeredFiniteElementsPtr->get<Unique_mi_tag>().upper_bound(
931 get_id_for_max_type<MBENTITYSET>(), (*fe_miit)->getFEUId()));
932
933 if (miit != hi_miit) {
934 std::vector<EntityHandle> ents;
935 ents.reserve(std::distance(miit, hi_miit));
936 for (; miit != hi_miit; ++miit)
937 ents.push_back((*miit)->getEnt());
938 int part = (*miit)->getPart();
939 CHKERR get_moab().tag_clear_data(th_Part, &*ents.begin(), ents.size(),
940 &part);
941 CHKERR get_moab().add_entities(meshset, &*ents.begin(), ents.size());
942 }
943 }
944
946}
947
948} // namespace MoFEM
#define FECoreFunctionBegin
Definition FECore.cpp:7
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
constexpr double a
@ QUIET
@ DEFAULT_VERBOSITY
@ VERBOSE
@ COL
@ DATA
@ ROW
MoFEMTypes
Those types control how functions respond on arguments, f.e. error handling.
@ MF_EXIST
@ MF_EXCL
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_NOT_FOUND
Definition definitions.h:33
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition definitions.h:34
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ BYCOL
@ BYDATA
@ BYROW
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
boost::function< MoFEMErrorCode( Interface &moab, const Field &field, const EntFiniteElement &fe, std::vector< EntityHandle > &adjacency)> ElementAdjacencyFunct
user adjacency function
MoFEMErrorCode getEntitiesByRefLevel(const BitRefLevel bit, const BitRefLevel mask, const EntityHandle meshset, const int verb=QUIET) const
add all ents from ref level given by bit to meshset
#define MOFEM_LOG(channel, severity)
Log.
static LoggerType & getLog(const std::string channel)
Get logger by channel.
auto bit
set bit
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
std::bitset< BITFEID_SIZE > BitFEId
Finite element Id.
Definition Types.hpp:43
std::bitset< BITFIELDID_SIZE > BitFieldId
Field Id.
Definition Types.hpp:42
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition Types.hpp:40
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
EntityHandle get_id_for_max_type()
EntityHandle get_id_for_min_type()
const EntityHandle no_handle
No entity handle is indicated by zero handle, i.e. root meshset.
Definition Common.hpp:12
constexpr auto field_name
Managing BitRefLevels.
MoFEMErrorCode get_finite_element_entities_by_dimension(const std::string name, int dim, Range &ents) const
get entities in the finite element by dimension
Definition FECore.cpp:286
MoFEMErrorCode modify_finite_element_off_field_row(const std::string &fe_name, const std::string name_row)
unset field row which finite element use
Definition FECore.cpp:219
MoFEMErrorCode check_number_of_ents_in_ents_finite_element() const
check data consistency in entsFiniteElements
Definition FECore.cpp:886
MoFEMErrorCode get_problem_finite_elements_entities(const std::string name, const std::string &fe_name, const EntityHandle meshset)
add finite elements to the meshset
Definition FECore.cpp:912
MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_col)
set field col which finite element use
Definition FECore.cpp:182
EntityHandle get_finite_element_meshset(const BitFEId id) const
Definition FECore.cpp:273
EntFiniteElement_multiIndex::index< Unique_mi_tag >::type::iterator get_fe_by_name_end(const std::string &fe_name) const
get end iterator of finite elements of given name (instead you can use IT_GET_FES_BY_NAME_FOR_LOOP(MF...
Definition FECore.cpp:848
MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)
set finite element field data
Definition FECore.cpp:138
MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)
add finite element
Definition FECore.cpp:43
MoFEMErrorCode add_ents_to_finite_element_by_MESHSET(const EntityHandle meshset, const std::string &name, const bool recursive=false)
add MESHSET element to finite element database given by name
Definition FECore.cpp:441
EntFiniteElement_multiIndex::index< Unique_mi_tag >::type::iterator get_fe_by_name_begin(const std::string &fe_name) const
get begin iterator of finite elements of given name (instead you can use IT_GET_FES_BY_NAME_FOR_LOOP(...
Definition FECore.cpp:837
std::string getBitFEIdName(const BitFEId id) const
Get field name.
Definition FECore.cpp:265
MoFEMErrorCode get_finite_element_entities_by_handle(const std::string name, Range &ents) const
get entities in the finite element by handle
Definition FECore.cpp:309
MoFEMErrorCode modify_finite_element_off_field_col(const std::string &fe_name, const std::string name_col)
unset field col which finite element use
Definition FECore.cpp:238
MoFEMErrorCode modify_finite_element_off_field_data(const std::string &fe_name, const std::string name_filed)
unset finite element field data
Definition FECore.cpp:201
MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)
set field row which finite element use
Definition FECore.cpp:160
bool check_finite_element(const std::string &name) const
Check if finite element is in database.
Definition FECore.cpp:35
MoFEMErrorCode add_ents_to_finite_element_by_dim(const EntityHandle meshset, const int dim, const std::string name, const bool recursive=true)
add entities to finite element
Definition FECore.cpp:346
MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)
Build finite elements.
Definition FECore.cpp:656
BitFEId getBitFEId(const std::string &fe_name) const
Get field Id.
Definition FECore.cpp:255
DEPRECATED MoFEMErrorCode add_ents_to_finite_element_EntType_by_bit_ref(const BitRefLevel &bit, const std::string &name, EntityType type, int verb=DEFAULT_VERBOSITY)
Definition FECore.cpp:386
MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)
build adjacencies
Definition FECore.cpp:743
MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle meshset, const EntityType type, const std::string name, const bool recursive=true)
add entities to finite element
Definition FECore.cpp:329
MoFEMErrorCode add_ents_to_finite_element_by_bit_ref(const BitRefLevel bit, const BitRefLevel mask, const std::string name, EntityType type, int verb=DEFAULT_VERBOSITY)
add TET entities from given refinement level to finite element database given by name
Definition FECore.cpp:405
MoFEMErrorCode buildFiniteElements(const boost::shared_ptr< FiniteElement > &fe, const Range *ents_ptr=NULL, int verb=DEFAULT_VERBOSITY)
Definition FECore.cpp:459
const FiniteElement * get_finite_element_structure(const std::string &name, enum MoFEMTypes bh=MF_EXCL) const
get finite element struture
Definition FECore.cpp:18
MoFEMErrorCode list_finite_elements() const
list finite elements in database
Definition FECore.cpp:320
MoFEMErrorCode modify_finite_element_adjacency_table(const std::string &fe_name, const EntityType type, ElementAdjacencyFunct function)
modify finite element table, only for advanced user
Definition FECore.cpp:117
MoFEMErrorCode get_finite_element_entities_by_type(const std::string name, EntityType type, Range &ents) const
get entities in the finite element by type
Definition FECore.cpp:296
UId getLocalUniqueIdCalculate() const
Generate UId for finite element entity.
UId getLocalUniqueIdCalculate()
Get the Local Unique Id Calculate object.
Finite element definition.