v0.14.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 partion 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 SETERRQ1(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 SETERRQ1(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 typedef FiniteElement_multiIndex::index<FiniteElement_name_mi_tag>::type
187 FiniteElements_by_name;
188 FiniteElements_by_name &finite_element_name_set =
189 finiteElements.get<FiniteElement_name_mi_tag>();
190 FiniteElements_by_name::iterator it_fe =
191 finite_element_name_set.find(fe_name);
192 if (it_fe == finite_element_name_set.end())
193 SETERRQ(mofemComm, MOFEM_OPERATION_UNSUCCESSFUL,
194 "this FiniteElement is there");
195 bool success = finite_element_name_set.modify(
196 it_fe, FiniteElement_col_change_bit_add(get_field_id(name_col)));
197 if (!success)
198 SETERRQ(mofemComm, MOFEM_OPERATION_UNSUCCESSFUL,
199 "modification unsuccessful");
201}
202
205 const std::string name_data) {
207 *buildMoFEM &= 1 << 0;
208 auto &finite_element_name_set =
209 finiteElements.get<FiniteElement_name_mi_tag>();
210 auto it_fe = finite_element_name_set.find(fe_name);
211 if (it_fe == finite_element_name_set.end())
212 SETERRQ(mofemComm, MOFEM_NOT_FOUND, "this FiniteElement is there");
213 bool success = finite_element_name_set.modify(
214 it_fe, FiniteElement_change_bit_off(get_field_id(name_data)));
215 if (!success)
216 SETERRQ(mofemComm, MOFEM_OPERATION_UNSUCCESSFUL,
217 "modification unsuccessful");
219}
220
223 const std::string name_row) {
225 *buildMoFEM &= 1 << 0;
226 auto &finite_element_name_set =
227 finiteElements.get<FiniteElement_name_mi_tag>();
228 auto it_fe = finite_element_name_set.find(fe_name);
229 if (it_fe == finite_element_name_set.end())
230 SETERRQ1(mofemComm, MOFEM_NOT_FOUND, "this < %s > is not there",
231 fe_name.c_str());
232 bool success = finite_element_name_set.modify(
233 it_fe, FiniteElement_row_change_bit_off(get_field_id(name_row)));
234 if (!success)
235 SETERRQ(mofemComm, MOFEM_OPERATION_UNSUCCESSFUL,
236 "modification unsuccessful");
238}
239
242 const std::string name_col) {
244 *buildMoFEM &= 1 << 0;
245 auto &finite_element_name_set =
246 finiteElements.get<FiniteElement_name_mi_tag>();
247 auto it_fe = finite_element_name_set.find(fe_name);
248 if (it_fe == finite_element_name_set.end())
249 SETERRQ(mofemComm, MOFEM_NOT_FOUND, "this FiniteElement is there");
250 bool success = finite_element_name_set.modify(
251 it_fe, FiniteElement_col_change_bit_off(get_field_id(name_col)));
252 if (!success)
253 SETERRQ(mofemComm, MOFEM_OPERATION_UNSUCCESSFUL,
254 "modification unsuccessful");
256}
257
258BitFEId Core::getBitFEId(const std::string &fe_name) const {
259 auto &fe_by_id = finiteElements.get<FiniteElement_name_mi_tag>();
260 auto miit = fe_by_id.find(fe_name);
261 if (miit == fe_by_id.end())
263 ("finite element < " + fe_name + " > not found (top tip: check spelling)")
264 .c_str());
265 return (*miit)->getId();
266}
267
268std::string Core::getBitFEIdName(const BitFEId id) const {
269 auto &fe_by_id = finiteElements.get<BitFEId_mi_tag>();
270 auto miit = fe_by_id.find(id);
271 if (miit == fe_by_id.end())
272 THROW_MESSAGE("finite element not found");
273 return (*miit)->getName();
274}
275
277 auto &fe_by_id = finiteElements.get<BitFEId_mi_tag>();
278 auto miit = fe_by_id.find(id);
279 if (miit == fe_by_id.end())
280 THROW_MESSAGE("finite element not found");
281 return (*miit)->meshset;
282}
283
284EntityHandle Core::get_finite_element_meshset(const std::string name) const {
285 return get_finite_element_meshset(getBitFEId(name));
286}
287
290 Range &ents) const {
291
293
294 EntityHandle meshset = get_finite_element_meshset(name);
295 CHKERR get_moab().get_entities_by_dimension(meshset, dim, ents, true);
297}
298
300 EntityType type,
301 Range &ents) const {
302
304
305 EntityHandle meshset = get_finite_element_meshset(name);
306 CHKERR get_moab().get_entities_by_type(meshset, type, ents, true);
307
309}
310
313 Range &ents) const {
314
316
317 EntityHandle meshset = get_finite_element_meshset(name);
318 CHKERR get_moab().get_entities_by_handle(meshset, ents, true);
319
321}
322
325 for (auto &fe : finiteElements.get<FiniteElement_name_mi_tag>())
326 MOFEM_LOG("SYNC", Sev::inform) << fe;
327
328 MOFEM_LOG_SYNCHRONISE(mofemComm);
330}
331
333 const EntityHandle meshset, const EntityType type, const std::string &name,
334 const bool recursive) {
335 *buildMoFEM &= 1 << 0;
338
339 idm = get_finite_element_meshset(getBitFEId(name));
340 Range ents;
341 CHKERR get_moab().get_entities_by_type(meshset, type, ents, recursive);
342 CHKERR getInterface<BitRefManager>()->setElementsBitRefLevel(ents);
343 CHKERR get_moab().add_entities(idm, ents);
344
346}
347
350 const int dim, const std::string &name,
351 const bool recursive) {
353 *buildMoFEM &= 1 << 0;
355 idm = get_finite_element_meshset(getBitFEId(name));
356 Range ents;
357 CHKERR get_moab().get_entities_by_dimension(meshset, dim, ents, recursive);
358 CHKERR getInterface<BitRefManager>()->setElementsBitRefLevel(ents);
359 CHKERR get_moab().add_entities(idm, ents);
361}
362
364 const Range &ents, const EntityType type, const std::string &name) {
366 *buildMoFEM &= 1 << 0;
368 idm = get_finite_element_meshset(getBitFEId(name));
369 CHKERR getInterface<BitRefManager>()->setElementsBitRefLevel(
370 ents.subset_by_type(type));
371 CHKERR get_moab().add_entities(idm, ents.subset_by_type(type));
373} // namespace MoFEM
374
377 const std::string &name) {
379 *buildMoFEM &= 1 << 0;
381 idm = get_finite_element_meshset(getBitFEId(name));
382 CHKERR getInterface<BitRefManager>()->setElementsBitRefLevel(
383 ents.subset_by_dimension(dim));
384 CHKERR get_moab().add_entities(idm, ents.subset_by_dimension(dim));
386}
387
390 const std::string &name,
391 EntityType type, int verb) {
393 CHKERR add_ents_to_finite_element_by_bit_ref(bit, BitRefLevel().set(), name,
394 type, verb);
395
397}
398
400 const BitRefLevel &bit, const BitRefLevel &mask, const std::string &name,
401 EntityType type, int verb) {
403 CHKERR add_ents_to_finite_element_by_bit_ref(bit, mask, name, type, verb);
404
406}
407
409 const BitRefLevel &bit, const BitRefLevel &mask, const std::string &name,
410 EntityType type, int verb) {
412
413 if (verb == -1)
414 verb = verbose;
415 *buildMoFEM &= 1 << 0;
416 const BitFEId id = getBitFEId(name);
417 const EntityHandle idm = get_finite_element_meshset(id);
418
419 auto &ref_MoFEMFiniteElement = refinedFiniteElements.get<Ent_mi_tag>();
420 auto miit = ref_MoFEMFiniteElement.lower_bound(get_id_for_min_type(type));
421 auto hi_miit = ref_MoFEMFiniteElement.upper_bound(get_id_for_max_type(type));
422
423 int nb_add_fes = 0;
424 for (; miit != hi_miit; miit++) {
425 const auto &bit2 = miit->get()->getBitRefLevel();
426 if ((bit2 & mask) != bit2)
427 continue;
428 if ((bit2 & bit).any()) {
429 EntityHandle ent = miit->get()->getEnt();
430 CHKERR get_moab().add_entities(idm, &ent, 1);
431 nb_add_fes++;
432 }
433 }
434
435 MOFEM_LOG("SYNC", Sev::inform)
436 << "Finite element " << name << " added. Nb. of elements added "
437 << nb_add_fes << " out of " << std::distance(miit, hi_miit);
438
439 MOFEM_LOG_SYNCHRONISE(mofemComm)
440
442}
443
445 const EntityHandle meshset, const std::string &name, const bool recursive) {
447 *buildMoFEM &= 1 << 0;
448 const BitFEId id = getBitFEId(name);
449 const EntityHandle idm = get_finite_element_meshset(id);
450 if (recursive == false) {
451 CHKERR get_moab().add_entities(idm, &meshset, 1);
452 } else {
453 Range meshsets;
454 CHKERR get_moab().get_entities_by_type(meshset, MBENTITYSET, meshsets,
455 false);
456 CHKERR get_moab().add_entities(idm, meshsets);
457 }
459}
460
462Core::buildFiniteElements(const boost::shared_ptr<FiniteElement> &fe,
463 const Range *ents_ptr, int verb) {
465 if (verb == DEFAULT_VERBOSITY)
466 verb = verbose;
467
468 if (verb > QUIET)
469 MOFEM_LOG("SYNC", Sev::verbose)
470 << "Build Finite Elements " << fe->getName();
471
472 auto &fields_by_id = fIelds.get<BitFieldId_mi_tag>();
473
474 // Get id of mofem fields for row, col and data
475 enum IntLoop { ROW = 0, COL, DATA, LAST };
476 std::array<BitFieldId, LAST> fe_fields = {fe.get()->getBitFieldIdRow(),
477 fe.get()->getBitFieldIdCol(),
478 fe.get()->getBitFieldIdData()};
479
480 // Get finite element meshset
481 EntityHandle meshset = get_finite_element_meshset(fe.get()->getId());
482
483 // Get entities from finite element meshset // if meshset
484 Range fe_ents;
485 CHKERR get_moab().get_entities_by_handle(meshset, fe_ents, false);
486
487 if (ents_ptr)
488 fe_ents = intersect(fe_ents, *ents_ptr);
489
490 // Map entity uid to pointers
491 typedef std::vector<boost::weak_ptr<EntFiniteElement>> VecOfWeakFEPtrs;
492 VecOfWeakFEPtrs processed_fes;
493 processed_fes.reserve(fe_ents.size());
494
495 int last_data_field_ents_view_size = 0;
496 int last_row_field_ents_view_size = 0;
497 int last_col_field_ents_view_size = 0;
498
499 // Entities adjacent to entities
500 std::vector<EntityHandle> adj_ents;
501
502 // Loop meshset finite element ents and add finite elements
503 for (Range::const_pair_iterator peit = fe_ents.const_pair_begin();
504 peit != fe_ents.const_pair_end(); peit++) {
505
506 const auto first = peit->first;
507 const auto second = peit->second;
508
509 // Find range of ref entities that is sequence
510 // note: iterator is a wrapper
511 // check if is in refinedFiniteElements database
512 auto ref_fe_miit =
513 refinedFiniteElements.get<Ent_mi_tag>().lower_bound(first);
514 if (ref_fe_miit == refinedFiniteElements.get<Ent_mi_tag>().end()) {
515 std::ostringstream ss;
516 ss << "refinedFiniteElements not in database ent = " << first << " type "
517 << type_from_handle(first) << " " << *fe;
518 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, ss.str().c_str());
519 }
520 auto hi_ref_fe_miit =
521 refinedFiniteElements.get<Ent_mi_tag>().upper_bound(second);
522
523 EntFiniteElement_multiIndex::iterator hint_p = entsFiniteElements.end();
524 for (; ref_fe_miit != hi_ref_fe_miit; ref_fe_miit++) {
525
526 // Add finite element to database
527 hint_p = entsFiniteElements.emplace_hint(
528 hint_p, boost::make_shared<EntFiniteElement>(*ref_fe_miit, fe));
529 processed_fes.emplace_back(*hint_p);
530 auto fe_raw_ptr = hint_p->get();
531
532 // Allocate space for entities view
533 bool row_as_data = false, col_as_row = false;
534 if (fe_fields[DATA] == fe_fields[ROW])
535 row_as_data = true;
536 if (fe_fields[ROW] == fe_fields[COL])
537 col_as_row = true;
538
539 fe_raw_ptr->getDataFieldEntsPtr()->reserve(
540 last_data_field_ents_view_size);
541
542 if (row_as_data) {
543 fe_raw_ptr->getRowFieldEntsPtr() = fe_raw_ptr->getDataFieldEntsPtr();
544 } else {
545 // row and col are different
546 if (fe_raw_ptr->getRowFieldEntsPtr() ==
547 fe_raw_ptr->getDataFieldEntsPtr())
548 fe_raw_ptr->getRowFieldEntsPtr() =
549 boost::make_shared<FieldEntity_vector_view>();
550 fe_raw_ptr->getRowFieldEntsPtr()->reserve(
551 last_row_field_ents_view_size);
552 }
553
554 if (row_as_data && col_as_row) {
555 fe_raw_ptr->getColFieldEntsPtr() = fe_raw_ptr->getDataFieldEntsPtr();
556 } else if (col_as_row) {
557 fe_raw_ptr->getColFieldEntsPtr() = fe_raw_ptr->getRowFieldEntsPtr();
558 } else {
559 if (
560
561 fe_raw_ptr->getColFieldEntsPtr() ==
562 fe_raw_ptr->getRowFieldEntsPtr() ||
563 fe_raw_ptr->getColFieldEntsPtr() ==
564 fe_raw_ptr->getDataFieldEntsPtr()
565
566 )
567 fe_raw_ptr->getColFieldEntsPtr() =
568 boost::make_shared<FieldEntity_vector_view>();
569 fe_raw_ptr->getColFieldEntsPtr()->reserve(
570 last_col_field_ents_view_size);
571 }
572
573 // Iterate over all field and check which one is on the element
574 for (unsigned int ii = 0; ii != BitFieldId().size(); ++ii) {
575
576 // Common field id for ROW, COL and DATA
577 BitFieldId id_common = 0;
578 // Check if the field (ii) is added to finite element
579 for (int ss = 0; ss < LAST; ss++) {
580 id_common |= fe_fields[ss] & BitFieldId().set(ii);
581 }
582 if (id_common.none())
583 continue;
584
585 // Find in database data associated with the field (ii)
586 const BitFieldId field_id = BitFieldId().set(ii);
587 auto miit = fields_by_id.find(field_id);
588 if (miit == fields_by_id.end())
589 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Field not found");
590 auto field_bit_number = (*miit)->getBitNumber();
591
592 // Loop over adjacencies of element and find field entities on those
593 // adjacencies, that create hash map map_uid_fe which is used later
594 const std::string field_name = miit->get()->getName();
595 const bool add_to_data = (field_id & fe_fields[DATA]).any();
596 const bool add_to_row = (field_id & fe_fields[ROW]).any();
597 const bool add_to_col = (field_id & fe_fields[COL]).any();
598
599 // Resolve entities on element, those entities are used to build tag
600 // with dof uids on finite element tag
601 adj_ents.clear();
602 CHKERR fe_raw_ptr->getElementAdjacency(*miit, adj_ents);
603
604 for (auto ent : adj_ents) {
605
606 auto dof_it = entsFields.get<Unique_mi_tag>().find(
607 FieldEntity::getLocalUniqueIdCalculate(field_bit_number, ent));
608 if (dof_it != entsFields.get<Unique_mi_tag>().end()) {
609
610 if (add_to_data) {
611 fe_raw_ptr->getDataFieldEntsPtr()->emplace_back(*dof_it);
612 }
613 if (add_to_row && !row_as_data) {
614 fe_raw_ptr->getRowFieldEntsPtr()->emplace_back(*dof_it);
615 }
616 if (add_to_col && !col_as_row) {
617 fe_raw_ptr->getColFieldEntsPtr()->emplace_back(*dof_it);
618 }
619
620 }
621 }
622 }
623
624 // Sort field ents by uid
625 auto uid_comp = [](const auto &a, const auto &b) {
626 return a.lock()->getLocalUniqueId() < b.lock()->getLocalUniqueId();
627 };
628
629 // Sort all views
630
631 // Data
632 sort(fe_raw_ptr->getDataFieldEntsPtr()->begin(),
633 fe_raw_ptr->getDataFieldEntsPtr()->end(), uid_comp);
634 last_data_field_ents_view_size =
635 fe_raw_ptr->getDataFieldEntsPtr()->size();
636
637 // Row
638 if (!row_as_data) {
639 sort(fe_raw_ptr->getRowFieldEntsPtr()->begin(),
640 fe_raw_ptr->getRowFieldEntsPtr()->end(), uid_comp);
641 last_row_field_ents_view_size =
642 fe_raw_ptr->getRowFieldEntsPtr()->size();
643 }
644
645 // Column
646 if (!col_as_row) {
647 sort(fe_raw_ptr->getColFieldEntsPtr()->begin(),
648 fe_raw_ptr->getColFieldEntsPtr()->end(), uid_comp);
649 last_col_field_ents_view_size =
650 fe_raw_ptr->getColFieldEntsPtr()->size();
651 }
652 }
653 }
654
656}
657
660
661 if (verb == DEFAULT_VERBOSITY)
662 verb = verbose;
663
664 // loop Finite Elements
665 for (auto &fe : finiteElements)
666 CHKERR buildFiniteElements(fe, NULL, verb);
667
668 if (verb > QUIET) {
669
670 auto &fe_ents = entsFiniteElements.get<Unique_mi_tag>();
671 for (auto &fe : finiteElements) {
672 auto miit = fe_ents.lower_bound(
674 auto hi_miit =
676 get_id_for_max_type<MBENTITYSET>(), fe->getFEUId()));
677 const auto count = std::distance(miit, hi_miit);
678 MOFEM_LOG("SYNC", Sev::inform)
679 << "Finite element " << fe->getName()
680 << " added. Nb. of elements added " << count;
681 MOFEM_LOG("SYNC", Sev::noisy) << *fe;
682
683 auto slg = MoFEM::LogManager::getLog("SYNC");
684 for (auto &field : fIelds) {
685 auto rec = slg.open_record(keywords::severity = Sev::verbose);
686 if (rec) {
687 logging::record_ostream strm(rec);
688 strm << "Field " << field->getName() << " on finite element: ";
689 if ((field->getId() & fe->getBitFieldIdRow()).any())
690 strm << "row ";
691 if ((field->getId() & fe->getBitFieldIdCol()).any())
692 strm << "columns ";
693 if ((field->getId() & fe->getBitFieldIdData()).any())
694 strm << "data";
695 strm.flush();
696 slg.push_record(boost::move(rec));
697 }
698 }
699 }
700
701 MOFEM_LOG_SYNCHRONISE(mofemComm);
702 }
703
704 *buildMoFEM |= 1 << 1;
706}
707
710 SETERRQ(mofemComm, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
712}
713
715 const Range *const ents_ptr,
716 int verb) {
718 if (verb == -1)
719 verb = verbose;
720
721 auto fe_miit = finiteElements.get<FiniteElement_name_mi_tag>().find(fe_name);
722 if (fe_miit == finiteElements.get<FiniteElement_name_mi_tag>().end())
723 SETERRQ1(mofemComm, MOFEM_NOT_FOUND, "Finite element <%s> not found",
724 fe_name.c_str());
725
726 CHKERR buildFiniteElements(*fe_miit, ents_ptr, verb);
727
728 if (verb >= VERBOSE) {
729 auto &fe_ents = entsFiniteElements.get<Unique_mi_tag>();
730 auto miit = fe_ents.lower_bound(
731 EntFiniteElement::getLocalUniqueIdCalculate(0, (*fe_miit)->getFEUId()));
732 auto hi_miit =
734 get_id_for_max_type<MBENTITYSET>(), (*fe_miit)->getFEUId()));
735 const auto count = std::distance(miit, hi_miit);
736 MOFEM_LOG("SYNC", Sev::inform) << "Finite element " << fe_name
737 << " added. Nb. of elements added " << count;
738 MOFEM_LOG_SYNCHRONISE(mofemComm);
739 }
740
741 *buildMoFEM |= 1 << 1;
743}
744
747 if (verb == DEFAULT_VERBOSITY)
748 verb = verbose;
749
750 if (!((*buildMoFEM) & BUILD_FIELD))
751 SETERRQ(mofemComm, MOFEM_NOT_FOUND, "field not build");
752 if (!((*buildMoFEM) & BUILD_FE))
753 SETERRQ(mofemComm, MOFEM_NOT_FOUND, "fe not build");
754 for (auto peit = ents.pair_begin(); peit != ents.pair_end(); ++peit) {
755 auto fit = entsFiniteElements.get<Ent_mi_tag>().lower_bound(peit->first);
756 auto hi_fit =
757 entsFiniteElements.get<Ent_mi_tag>().upper_bound(peit->second);
758 for (; fit != hi_fit; ++fit) {
759 if ((*fit)->getBitFieldIdRow().none() &&
760 (*fit)->getBitFieldIdCol().none() &&
761 (*fit)->getBitFieldIdData().none())
762 continue;
763 int by = BYROW;
764 if ((*fit)->getBitFieldIdRow() != (*fit)->getBitFieldIdCol())
765 by |= BYCOL;
766 if ((*fit)->getBitFieldIdRow() != (*fit)->getBitFieldIdData())
767 by |= BYDATA;
769 auto hint = entFEAdjacencies.end();
770 for (auto e : *(*fit)->getRowFieldEntsPtr()) {
771 hint = entFEAdjacencies.emplace_hint(hint, e.lock(), *fit);
772 bool success = entFEAdjacencies.modify(hint, modify_row);
773 if (!success)
774 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
775 "modification unsuccessful");
776 }
777 if ((*fit)->getBitFieldIdRow() != (*fit)->getBitFieldIdCol()) {
778 int by = BYCOL;
779 if ((*fit)->getBitFieldIdCol() != (*fit)->getBitFieldIdData())
780 by |= BYDATA;
782 auto hint = entFEAdjacencies.end();
783 for (auto e : *(*fit)->getColFieldEntsPtr()) {
784 hint = entFEAdjacencies.emplace_hint(hint, e.lock(), *fit);
785 bool success = entFEAdjacencies.modify(hint, modify_col);
786 if (!success)
787 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
788 "modification unsuccessful");
789 }
790 }
791 if ((*fit)->getBitFieldIdRow() != (*fit)->getBitFieldIdData() ||
792 (*fit)->getBitFieldIdCol() != (*fit)->getBitFieldIdData()) {
794 BYDATA);
795 auto hint = entFEAdjacencies.end();
796 for (auto &e : (*fit)->getDataFieldEnts()) {
797 hint = entFEAdjacencies.emplace_hint(hint, e.lock(), *fit);
798 bool success = entFEAdjacencies.modify(hint, modify_data);
799 if (!success)
800 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
801 "modification unsuccessful");
802 }
803 }
804 }
805 }
806
807 if (verb >= VERBOSE) {
808 MOFEM_LOG("WORLD", Sev::inform)
809 << "Number of adjacencies " << entFEAdjacencies.size();
810 MOFEM_LOG_SYNCHRONISE(mofemComm)
811 }
812
813 *buildMoFEM |= 1 << 2;
815}
816
818 const BitRefLevel &mask, int verb) {
820 if (verb == -1)
821 verb = verbose;
822 Range ents;
823 CHKERR BitRefManager(*this).getEntitiesByRefLevel(bit, mask, ents);
824
825 CHKERR build_adjacencies(ents, verb);
826
828}
831 if (verb == -1)
832 verb = verbose;
833 CHKERR build_adjacencies(bit, BitRefLevel().set(), verb);
834
836}
837
838EntFiniteElement_multiIndex::index<Unique_mi_tag>::type::iterator
839Core::get_fe_by_name_begin(const std::string &fe_name) const {
840 auto miit = finiteElements.get<FiniteElement_name_mi_tag>().find(fe_name);
841 if (miit != finiteElements.get<FiniteElement_name_mi_tag>().end()) {
842 return entsFiniteElements.get<Unique_mi_tag>().lower_bound(
843 EntFiniteElement::getLocalUniqueIdCalculate(0, (*miit)->getFEUId()));
844 } else {
845 return entsFiniteElements.get<Unique_mi_tag>().end();
846 }
847}
848
849EntFiniteElement_multiIndex::index<Unique_mi_tag>::type::iterator
850Core::get_fe_by_name_end(const std::string &fe_name) const {
851 auto miit = finiteElements.get<FiniteElement_name_mi_tag>().find(fe_name);
852 if (miit != finiteElements.get<FiniteElement_name_mi_tag>().end()) {
853 return entsFiniteElements.get<Unique_mi_tag>().upper_bound(
855 get_id_for_max_type<MBENTITYSET>(), (*miit)->getFEUId()));
856 } else {
857 return entsFiniteElements.get<Unique_mi_tag>().end();
858 }
859}
860
862 const std::string &name) const {
864 FiniteElement_multiIndex::index<FiniteElement_name_mi_tag>::type::iterator it;
865 it = finiteElements.get<FiniteElement_name_mi_tag>().find(name);
866 if (it == finiteElements.get<FiniteElement_name_mi_tag>().end()) {
867 SETERRQ1(mofemComm, 1, "finite element not found < %s >", name.c_str());
868 }
869 EntityHandle meshset = (*it)->getMeshset();
870
871 int num_entities;
872 CHKERR get_moab().get_number_entities_by_handle(meshset, num_entities);
873
874 auto counts_fes = [&]() {
875 return std::distance(get_fe_by_name_begin((*it)->getName()),
876 get_fe_by_name_end((*it)->getName()));
877 };
878
879 if (counts_fes() != static_cast<size_t>(num_entities)) {
880 SETERRQ1(mofemComm, MOFEM_DATA_INCONSISTENCY,
881 "not equal number of entities in meshset and finite elements "
882 "multiindex < %s >",
883 (*it)->getName().c_str());
884 }
886}
887
890 FiniteElement_multiIndex::index<FiniteElement_name_mi_tag>::type::iterator it;
891 it = finiteElements.get<FiniteElement_name_mi_tag>().begin();
892 for (; it != finiteElements.get<FiniteElement_name_mi_tag>().end(); it++) {
893 EntityHandle meshset = (*it)->getMeshset();
894
895 int num_entities;
896 CHKERR get_moab().get_number_entities_by_handle(meshset, num_entities);
897
898 auto counts_fes = [&]() {
899 return std::distance(get_fe_by_name_begin((*it)->getName()),
900 get_fe_by_name_end((*it)->getName()));
901 };
902
903 if (counts_fes() != static_cast<size_t>(num_entities)) {
904 SETERRQ1(mofemComm, MOFEM_DATA_INCONSISTENCY,
905 "not equal number of entities in meshset and finite elements "
906 "multiindex < %s >",
907 (*it)->getName().c_str());
908 }
909 }
911}
912
914Core::get_problem_finite_elements_entities(const std::string problem_name,
915 const std::string &fe_name,
916 const EntityHandle meshset) {
918 auto &prb = pRoblems.get<Problem_mi_tag>();
919 auto p_miit = prb.find(problem_name);
920 if (p_miit == prb.end())
921 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
922 "No such problem like < %s >", problem_name.c_str());
923
924 auto fe_miit = finiteElements.get<FiniteElement_name_mi_tag>().find(fe_name);
925 if (fe_miit != finiteElements.get<FiniteElement_name_mi_tag>().end()) {
926 auto miit =
927 p_miit->numeredFiniteElementsPtr->get<Unique_mi_tag>().lower_bound(
929 0, (*fe_miit)->getFEUId()));
930 auto hi_miit =
931 p_miit->numeredFiniteElementsPtr->get<Unique_mi_tag>().upper_bound(
933 get_id_for_max_type<MBENTITYSET>(), (*fe_miit)->getFEUId()));
934
935 if (miit != hi_miit) {
936 std::vector<EntityHandle> ents;
937 ents.reserve(std::distance(miit, hi_miit));
938 for (; miit != hi_miit; ++miit)
939 ents.push_back((*miit)->getEnt());
940 int part = (*miit)->getPart();
941 CHKERR get_moab().tag_clear_data(th_Part, &*ents.begin(), ents.size(),
942 &part);
943 CHKERR get_moab().add_entities(meshset, &*ents.begin(), ents.size());
944 }
945 }
946
948}
949
950} // namespace MoFEM
static Index< 'p', 3 > p
#define FECoreFunctionBegin
Definition: FECore.cpp:7
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:345
constexpr double a
@ QUIET
Definition: definitions.h:208
@ DEFAULT_VERBOSITY
Definition: definitions.h:207
@ VERBOSE
Definition: definitions.h:209
@ COL
Definition: definitions.h:123
@ DATA
Definition: definitions.h:123
@ ROW
Definition: definitions.h:123
MoFEMTypes
Those types control how functions respond on arguments, f.e. error handling.
Definition: definitions.h:97
@ MF_EXIST
Definition: definitions.h:100
@ MF_EXCL
Definition: definitions.h:99
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ 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()
Definition: definitions.h:416
@ BYCOL
Definition: definitions.h:132
@ BYDATA
Definition: definitions.h:133
@ BYROW
Definition: definitions.h:131
#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
const int dim
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.
Definition: LogManager.hpp:308
static LoggerType & getLog(const std::string channel)
Get logger by channel.
Definition: LogManager.cpp:395
auto bit
set bit
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
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
auto type_from_handle(const EntityHandle h)
get type from entity handle
Definition: Templates.hpp:1869
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:289
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:222
MoFEMErrorCode check_number_of_ents_in_ents_finite_element() const
check data consistency in entsFiniteElements
Definition: FECore.cpp:888
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:914
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:276
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:850
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:444
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:839
std::string getBitFEIdName(const BitFEId id) const
Get field name.
Definition: FECore.cpp:268
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:349
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:312
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:332
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:241
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:204
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:408
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 build_finite_elements(int verb=DEFAULT_VERBOSITY)
Build finite elements.
Definition: FECore.cpp:658
BitFEId getBitFEId(const std::string &fe_name) const
Get field Id.
Definition: FECore.cpp:258
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:389
MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)
build adjacencies
Definition: FECore.cpp:745
MoFEMErrorCode buildFiniteElements(const boost::shared_ptr< FiniteElement > &fe, const Range *ents_ptr=NULL, int verb=DEFAULT_VERBOSITY)
Definition: FECore.cpp:462
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:323
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:299
UId getLocalUniqueIdCalculate() const
Generate UId for finite element entity.
UId getLocalUniqueIdCalculate()
Get the Local Unique Id Calculate object.
Finite element definition.