v0.13.1
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
6#include <MoFEM.hpp>
7
8#define FECoreFunctionBegin \
9 MoFEMFunctionBegin; \
10 MOFEM_LOG_CHANNEL("WORLD"); \
11 MOFEM_LOG_CHANNEL("SYNC"); \
12 MOFEM_LOG_FUNCTION(); \
13 MOFEM_LOG_TAG("WORLD", "FECore"); \
14 MOFEM_LOG_TAG("SYNC", "FECore");
15
16namespace MoFEM {
17
18bool Core::check_finite_element(const std::string &name) const {
19 typedef FiniteElement_multiIndex::index<FiniteElement_name_mi_tag>::type
20 FeSetByName;
21 const FeSetByName &set = finiteElements.get<FiniteElement_name_mi_tag>();
22 FeSetByName::iterator miit = set.find(name);
23 if (miit == set.end())
24 return false;
25 return true;
26}
27
28MoFEMErrorCode Core::add_finite_element(const std::string &fe_name,
29 enum MoFEMTypes bh, int verb) {
31 *buildMoFEM &= 1 << 0;
32 if (verb == -1) {
33 verb = verbose;
34 }
35
36 // Add finite element meshset to partion meshset. In case of no elements
37 // on processor part, when mesh file is read, finite element meshset is
38 // prevented from deletion by moab reader.
39 auto add_meshset_to_partition = [&](auto meshset) {
41 const void *tag_vals[] = {&rAnk};
42 ParallelComm *pcomm = ParallelComm::get_pcomm(
43 &get_moab(), get_basic_entity_data_ptr()->pcommID);
44 Tag part_tag = pcomm->part_tag();
45 Range tagged_sets;
46 CHKERR get_moab().get_entities_by_type_and_tag(0, MBENTITYSET, &part_tag,
47 tag_vals, 1, tagged_sets,
48 moab::Interface::UNION);
49 for (auto s : tagged_sets)
50 CHKERR get_moab().add_entities(s, &meshset, 1);
52 };
53
54 auto &finite_element_name_set =
55 finiteElements.get<FiniteElement_name_mi_tag>();
56 auto it_fe = finite_element_name_set.find(fe_name);
57
58 if (bh == MF_EXCL) {
59 if (it_fe != finite_element_name_set.end()) {
60 SETERRQ1(mofemComm, MOFEM_NOT_FOUND, "this < %s > is there",
61 fe_name.c_str());
62 }
63
64 } else {
65 if (it_fe != finite_element_name_set.end())
67 }
68 EntityHandle meshset;
69 CHKERR get_moab().create_meshset(MESHSET_SET, meshset);
70 CHKERR add_meshset_to_partition(meshset);
71
72 // id
73 int fe_shift = 0;
74 for (; finiteElements.get<BitFEId_mi_tag>().find(BitFEId().set(fe_shift)) !=
75 finiteElements.get<BitFEId_mi_tag>().end();
76 ++fe_shift) {
77 }
78
79 auto id = BitFEId().set(fe_shift);
80 CHKERR get_moab().tag_set_data(th_FEId, &meshset, 1, &id);
81
82 // id name
83 void const *tag_data[] = {fe_name.c_str()};
84 int tag_sizes[1];
85 tag_sizes[0] = fe_name.size();
86 CHKERR get_moab().tag_set_by_ptr(th_FEName, &meshset, 1, tag_data, tag_sizes);
87
88 // add FiniteElement
89 auto p = finiteElements.insert(
90 boost::shared_ptr<FiniteElement>(new FiniteElement(moab, meshset)));
91 if (!p.second)
92 SETERRQ(mofemComm, MOFEM_OPERATION_UNSUCCESSFUL,
93 "FiniteElement not inserted");
94
95 if (verb > QUIET)
96 MOFEM_LOG("WORLD", Sev::inform) << "Add finite element " << fe_name;
97
99}
100
103 const EntityType type,
104 ElementAdjacencyFunct function) {
106 *buildMoFEM &= 1 << 0;
107 typedef FiniteElement_multiIndex::index<FiniteElement_name_mi_tag>::type
108 FiniteElements_by_name;
109 FiniteElements_by_name &finite_element_name_set =
110 finiteElements.get<FiniteElement_name_mi_tag>();
111 FiniteElements_by_name::iterator it_fe =
112 finite_element_name_set.find(fe_name);
113 if (it_fe == finite_element_name_set.end())
114 SETERRQ(mofemComm, MOFEM_NOT_FOUND,
115 "This finite element is not defined (advise: check spelling)");
116 boost::shared_ptr<FiniteElement> fe;
117 fe = *it_fe;
118 fe->elementAdjacencyTable[type] = function;
120}
121
124 const std::string &name_data) {
126 *buildMoFEM &= 1 << 0;
127 typedef FiniteElement_multiIndex::index<FiniteElement_name_mi_tag>::type
128 FiniteElements_by_name;
129 FiniteElements_by_name &finite_element_name_set =
130 finiteElements.get<FiniteElement_name_mi_tag>();
131 FiniteElements_by_name::iterator it_fe =
132 finite_element_name_set.find(fe_name);
133 if (it_fe == finite_element_name_set.end())
134 SETERRQ(mofemComm, MOFEM_NOT_FOUND,
135 "This finite element is not defined (advise: check spelling)");
136 bool success = finite_element_name_set.modify(
137 it_fe, FiniteElement_change_bit_add(get_field_id(name_data)));
138 if (!success)
139 SETERRQ(mofemComm, MOFEM_OPERATION_UNSUCCESSFUL,
140 "modification unsuccessful");
142}
143
146 const std::string &name_row) {
148 *buildMoFEM &= 1 << 0;
149 typedef FiniteElement_multiIndex::index<FiniteElement_name_mi_tag>::type
150 FiniteElements_by_name;
151 FiniteElements_by_name &finite_element_name_set =
152 finiteElements.get<FiniteElement_name_mi_tag>();
153 FiniteElements_by_name::iterator it_fe =
154 finite_element_name_set.find(fe_name);
155 if (it_fe == finite_element_name_set.end())
156 SETERRQ1(mofemComm, MOFEM_NOT_FOUND, "this < %s > is not there",
157 fe_name.c_str());
158 bool success = finite_element_name_set.modify(
159 it_fe, FiniteElement_row_change_bit_add(get_field_id(name_row)));
160 if (!success)
161 SETERRQ(mofemComm, MOFEM_OPERATION_UNSUCCESSFUL,
162 "modification unsuccessful");
164}
165
168 const std::string &name_col) {
170 *buildMoFEM &= 1 << 0;
171 typedef FiniteElement_multiIndex::index<FiniteElement_name_mi_tag>::type
172 FiniteElements_by_name;
173 FiniteElements_by_name &finite_element_name_set =
174 finiteElements.get<FiniteElement_name_mi_tag>();
175 FiniteElements_by_name::iterator it_fe =
176 finite_element_name_set.find(fe_name);
177 if (it_fe == finite_element_name_set.end())
178 SETERRQ(mofemComm, MOFEM_OPERATION_UNSUCCESSFUL,
179 "this FiniteElement is there");
180 bool success = finite_element_name_set.modify(
181 it_fe, FiniteElement_col_change_bit_add(get_field_id(name_col)));
182 if (!success)
183 SETERRQ(mofemComm, MOFEM_OPERATION_UNSUCCESSFUL,
184 "modification unsuccessful");
186}
187
190 const std::string &name_data) {
192 *buildMoFEM &= 1 << 0;
193 auto &finite_element_name_set =
194 finiteElements.get<FiniteElement_name_mi_tag>();
195 auto it_fe = finite_element_name_set.find(fe_name);
196 if (it_fe == finite_element_name_set.end())
197 SETERRQ(mofemComm, MOFEM_NOT_FOUND, "this FiniteElement is there");
198 bool success = finite_element_name_set.modify(
199 it_fe, FiniteElement_change_bit_off(get_field_id(name_data)));
200 if (!success)
201 SETERRQ(mofemComm, MOFEM_OPERATION_UNSUCCESSFUL,
202 "modification unsuccessful");
204}
205
208 const std::string &name_row) {
210 *buildMoFEM &= 1 << 0;
211 auto &finite_element_name_set =
212 finiteElements.get<FiniteElement_name_mi_tag>();
213 auto it_fe = finite_element_name_set.find(fe_name);
214 if (it_fe == finite_element_name_set.end())
215 SETERRQ1(mofemComm, MOFEM_NOT_FOUND, "this < %s > is not there",
216 fe_name.c_str());
217 bool success = finite_element_name_set.modify(
218 it_fe, FiniteElement_row_change_bit_off(get_field_id(name_row)));
219 if (!success)
220 SETERRQ(mofemComm, MOFEM_OPERATION_UNSUCCESSFUL,
221 "modification unsuccessful");
223}
224
227 const std::string &name_col) {
229 *buildMoFEM &= 1 << 0;
230 auto &finite_element_name_set =
231 finiteElements.get<FiniteElement_name_mi_tag>();
232 auto it_fe = finite_element_name_set.find(fe_name);
233 if (it_fe == finite_element_name_set.end())
234 SETERRQ(mofemComm, MOFEM_NOT_FOUND, "this FiniteElement is there");
235 bool success = finite_element_name_set.modify(
236 it_fe, FiniteElement_col_change_bit_off(get_field_id(name_col)));
237 if (!success)
238 SETERRQ(mofemComm, MOFEM_OPERATION_UNSUCCESSFUL,
239 "modification unsuccessful");
241}
242
243BitFEId Core::getBitFEId(const std::string &name) const {
244 auto &fe_by_id = finiteElements.get<FiniteElement_name_mi_tag>();
245 auto miit = fe_by_id.find(name);
246 if (miit == fe_by_id.end())
248 ("finite element < " + name + " > not found (top tip: check spelling)")
249 .c_str());
250 return (*miit)->getId();
251}
252
253std::string Core::getBitFEIdName(const BitFEId id) const {
254 auto &fe_by_id = finiteElements.get<BitFEId_mi_tag>();
255 auto miit = fe_by_id.find(id);
256 if (miit == fe_by_id.end())
257 THROW_MESSAGE("finite element not found");
258 return (*miit)->getName();
259}
260
262 auto &fe_by_id = finiteElements.get<BitFEId_mi_tag>();
263 auto miit = fe_by_id.find(id);
264 if (miit == fe_by_id.end())
265 THROW_MESSAGE("finite element not found");
266 return (*miit)->meshset;
267}
268
269EntityHandle Core::get_finite_element_meshset(const std::string &name) const {
270 return get_finite_element_meshset(getBitFEId(name));
271}
272
275 Range &ents) const {
276
278
279 EntityHandle meshset = get_finite_element_meshset(name);
280 CHKERR get_moab().get_entities_by_dimension(meshset, dim, ents, true);
282}
283
285 EntityType type,
286 Range &ents) const {
287
289
290 EntityHandle meshset = get_finite_element_meshset(name);
291 CHKERR get_moab().get_entities_by_type(meshset, type, ents, true);
292
294}
295
298 Range &ents) const {
299
301
302 EntityHandle meshset = get_finite_element_meshset(name);
303 CHKERR get_moab().get_entities_by_handle(meshset, ents, true);
304
306}
307
310 for (auto &fe : finiteElements.get<FiniteElement_name_mi_tag>())
311 MOFEM_LOG("SYNC", Sev::inform) << fe;
312
313 MOFEM_LOG_SYNCHRONISE(mofemComm);
315}
316
318 const EntityHandle meshset, const EntityType type, const std::string &name,
319 const bool recursive) {
320 *buildMoFEM &= 1 << 0;
323
324 idm = get_finite_element_meshset(getBitFEId(name));
325 Range ents;
326 CHKERR get_moab().get_entities_by_type(meshset, type, ents, recursive);
327 CHKERR getInterface<BitRefManager>()->setElementsBitRefLevel(ents);
328 CHKERR get_moab().add_entities(idm, ents);
329
331}
332
335 const int dim, const std::string &name,
336 const bool recursive) {
338 *buildMoFEM &= 1 << 0;
340 idm = get_finite_element_meshset(getBitFEId(name));
341 Range ents;
342 CHKERR get_moab().get_entities_by_dimension(meshset, dim, ents, recursive);
343 CHKERR getInterface<BitRefManager>()->setElementsBitRefLevel(ents);
344 CHKERR get_moab().add_entities(idm, ents);
346}
347
349 const Range &ents, const EntityType type, const std::string &name) {
351 *buildMoFEM &= 1 << 0;
353 idm = get_finite_element_meshset(getBitFEId(name));
354 CHKERR getInterface<BitRefManager>()->setElementsBitRefLevel(
355 ents.subset_by_type(type));
356 CHKERR get_moab().add_entities(idm, ents.subset_by_type(type));
358} // namespace MoFEM
359
362 const std::string &name) {
364 *buildMoFEM &= 1 << 0;
366 idm = get_finite_element_meshset(getBitFEId(name));
367 CHKERR getInterface<BitRefManager>()->setElementsBitRefLevel(
368 ents.subset_by_dimension(dim));
369 CHKERR get_moab().add_entities(idm, ents.subset_by_dimension(dim));
371}
372
375 const std::string &name,
376 EntityType type, int verb) {
378 CHKERR add_ents_to_finite_element_by_bit_ref(bit, BitRefLevel().set(), name,
379 type, verb);
380
382}
383
385 const BitRefLevel &bit, const BitRefLevel &mask, const std::string &name,
386 EntityType type, int verb) {
388 CHKERR add_ents_to_finite_element_by_bit_ref(bit, mask, name, type, verb);
389
391}
392
394 const BitRefLevel &bit, const BitRefLevel &mask, const std::string &name,
395 EntityType type, int verb) {
397
398 if (verb == -1)
399 verb = verbose;
400 *buildMoFEM &= 1 << 0;
401 const BitFEId id = getBitFEId(name);
402 const EntityHandle idm = get_finite_element_meshset(id);
403
404 auto &ref_MoFEMFiniteElement = refinedFiniteElements.get<Ent_mi_tag>();
405 auto miit = ref_MoFEMFiniteElement.lower_bound(get_id_for_min_type(type));
406 auto hi_miit = ref_MoFEMFiniteElement.upper_bound(get_id_for_max_type(type));
407
408 int nb_add_fes = 0;
409 for (; miit != hi_miit; miit++) {
410 const auto &bit2 = miit->get()->getBitRefLevel();
411 if ((bit2 & mask) != bit2)
412 continue;
413 if ((bit2 & bit).any()) {
414 EntityHandle ent = miit->get()->getEnt();
415 CHKERR get_moab().add_entities(idm, &ent, 1);
416 nb_add_fes++;
417 }
418 }
419
420 MOFEM_LOG("SYNC", Sev::inform)
421 << "Finite element " << name << " added. Nb. of elements added "
422 << nb_add_fes << " out of " << std::distance(miit, hi_miit);
423
424 MOFEM_LOG_SYNCHRONISE(mofemComm)
425
427}
428
430 const EntityHandle meshset, const std::string &name, const bool recursive) {
432 *buildMoFEM &= 1 << 0;
433 const BitFEId id = getBitFEId(name);
434 const EntityHandle idm = get_finite_element_meshset(id);
435 if (recursive == false) {
436 CHKERR get_moab().add_entities(idm, &meshset, 1);
437 } else {
438 Range meshsets;
439 CHKERR get_moab().get_entities_by_type(meshset, MBENTITYSET, meshsets,
440 false);
441 CHKERR get_moab().add_entities(idm, meshsets);
442 }
444}
445
447Core::buildFiniteElements(const boost::shared_ptr<FiniteElement> &fe,
448 const Range *ents_ptr, int verb) {
450 if (verb == DEFAULT_VERBOSITY)
451 verb = verbose;
452
453 if (verb > QUIET)
454 MOFEM_LOG("SYNC", Sev::verbose)
455 << "Build Finite Elements " << fe->getName();
456
457 auto &fields_by_id = fIelds.get<BitFieldId_mi_tag>();
458
459 // Get id of mofem fields for row, col and data
460 enum IntLoop { ROW = 0, COL, DATA, LAST };
461 std::array<BitFieldId, LAST> fe_fields = {fe.get()->getBitFieldIdRow(),
462 fe.get()->getBitFieldIdCol(),
463 fe.get()->getBitFieldIdData()};
464
465 // Get finite element meshset
466 EntityHandle meshset = get_finite_element_meshset(fe.get()->getId());
467
468 // Get entities from finite element meshset // if meshset
469 Range fe_ents;
470 CHKERR get_moab().get_entities_by_handle(meshset, fe_ents, false);
471
472 if (ents_ptr)
473 fe_ents = intersect(fe_ents, *ents_ptr);
474
475 // Map entity uid to pointers
476 typedef std::vector<boost::weak_ptr<EntFiniteElement>> VecOfWeakFEPtrs;
477 typedef std::map<const UId *, VecOfWeakFEPtrs> MapEntUIdAndVecOfWeakFEPtrs;
478 MapEntUIdAndVecOfWeakFEPtrs ent_uid_and_fe_vec;
479 VecOfWeakFEPtrs processed_fes;
480 processed_fes.reserve(fe_ents.size());
481
482 int last_data_field_ents_view_size = 0;
483 int last_row_field_ents_view_size = 0;
484 int last_col_field_ents_view_size = 0;
485
486 // Entities adjacent to entities
487 std::vector<EntityHandle> adj_ents;
488
489 // Loop meshset finite element ents and add finite elements
490 for (Range::const_pair_iterator peit = fe_ents.const_pair_begin();
491 peit != fe_ents.const_pair_end(); peit++) {
492
493 const auto first = peit->first;
494 const auto second = peit->second;
495
496 // Find range of ref entities that is sequence
497 // note: iterator is a wrapper
498 // check if is in refinedFiniteElements database
499 auto ref_fe_miit =
500 refinedFiniteElements.get<Ent_mi_tag>().lower_bound(first);
501 if (ref_fe_miit == refinedFiniteElements.get<Ent_mi_tag>().end()) {
502 std::ostringstream ss;
503 ss << "refinedFiniteElements not in database ent = " << first << " type "
504 << type_from_handle(first) << " " << *fe;
505 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, ss.str().c_str());
506 }
507 auto hi_ref_fe_miit =
508 refinedFiniteElements.get<Ent_mi_tag>().upper_bound(second);
509
510 EntFiniteElement_multiIndex::iterator hint_p = entsFiniteElements.end();
511 for (; ref_fe_miit != hi_ref_fe_miit; ref_fe_miit++) {
512
513 // Add finite element to database
514 hint_p = entsFiniteElements.emplace_hint(
515 hint_p, boost::make_shared<EntFiniteElement>(*ref_fe_miit, fe));
516 processed_fes.emplace_back(*hint_p);
517 auto fe_raw_ptr = hint_p->get();
518
519 // Allocate space for etities view
520 bool row_as_data = false, col_as_row = false;
521 if (fe_fields[DATA] == fe_fields[ROW])
522 row_as_data = true;
523 if (fe_fields[ROW] == fe_fields[COL])
524 col_as_row = true;
525
526 fe_raw_ptr->getDataFieldEntsPtr()->reserve(
527 last_data_field_ents_view_size);
528
529 if (row_as_data) {
530 fe_raw_ptr->getRowFieldEntsPtr() = fe_raw_ptr->getDataFieldEntsPtr();
531 } else {
532 // row and col are diffent
533 if (fe_raw_ptr->getRowFieldEntsPtr() ==
534 fe_raw_ptr->getDataFieldEntsPtr())
535 fe_raw_ptr->getRowFieldEntsPtr() =
536 boost::make_shared<FieldEntity_vector_view>();
537 fe_raw_ptr->getRowFieldEntsPtr()->reserve(
538 last_row_field_ents_view_size);
539 }
540
541 if (row_as_data && col_as_row) {
542 fe_raw_ptr->getColFieldEntsPtr() = fe_raw_ptr->getDataFieldEntsPtr();
543 } else if (col_as_row) {
544 fe_raw_ptr->getColFieldEntsPtr() = fe_raw_ptr->getRowFieldEntsPtr();
545 } else {
546 if (
547
548 fe_raw_ptr->getColFieldEntsPtr() ==
549 fe_raw_ptr->getRowFieldEntsPtr() ||
550 fe_raw_ptr->getColFieldEntsPtr() ==
551 fe_raw_ptr->getDataFieldEntsPtr()
552
553 )
554 fe_raw_ptr->getColFieldEntsPtr() =
555 boost::make_shared<FieldEntity_vector_view>();
556 fe_raw_ptr->getColFieldEntsPtr()->reserve(
557 last_col_field_ents_view_size);
558 }
559
560 // Iterate over all field and check which one is on the element
561 for (unsigned int ii = 0; ii != BitFieldId().size(); ++ii) {
562
563 // Common field id for ROW, COL and DATA
564 BitFieldId id_common = 0;
565 // Check if the field (ii) is added to finite element
566 for (int ss = 0; ss < LAST; ss++) {
567 id_common |= fe_fields[ss] & BitFieldId().set(ii);
568 }
569 if (id_common.none())
570 continue;
571
572 // Find in database data associated with the field (ii)
573 const BitFieldId field_id = BitFieldId().set(ii);
574 auto miit = fields_by_id.find(field_id);
575 if (miit == fields_by_id.end())
576 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Field not found");
577 auto field_bit_number = (*miit)->getBitNumber();
578
579 // Loop over adjacencies of element and find field entities on those
580 // adjacencies, that create hash map map_uid_fe which is used later
581 const std::string field_name = miit->get()->getName();
582 const bool add_to_data = (field_id & fe_fields[DATA]).any();
583 const bool add_to_row = (field_id & fe_fields[ROW]).any();
584 const bool add_to_col = (field_id & fe_fields[COL]).any();
585
586 // Resolve entities on element, those entities are used to build tag
587 // with dof uids on finite element tag
588 adj_ents.clear();
589 CHKERR fe_raw_ptr->getElementAdjacency(*miit, adj_ents);
590
591 for(auto ent : adj_ents) {
592
593 auto dof_it = entsFields.get<Unique_mi_tag>().find(
594 FieldEntity::getLocalUniqueIdCalculate(field_bit_number, ent));
595 if(dof_it!=entsFields.get<Unique_mi_tag>().end()) {
596 // Add entity to map with key entity uids pointers and data
597 // finite elements weak ptrs. I using pointers to uids instead
598 // uids because this is faster.
599 const UId *uid_ptr = &(dof_it->get()->getLocalUniqueId());
600 auto &fe_vec = ent_uid_and_fe_vec[uid_ptr];
601 if (add_to_data) {
602 fe_raw_ptr->getDataFieldEntsPtr()->emplace_back(*dof_it);
603 }
604 if (add_to_row && !row_as_data) {
605 fe_raw_ptr->getRowFieldEntsPtr()->emplace_back(*dof_it);
606 }
607 if (add_to_col && !col_as_row) {
608 fe_raw_ptr->getColFieldEntsPtr()->emplace_back(*dof_it);
609 }
610
611 // add finite element to processed list
612 fe_vec.emplace_back(*hint_p);
613 }
614
615 }
616
617
618 }
619
620 // Sort field ents by uid
621 auto uid_comp = [](const auto &a, const auto &b) {
622 return a.lock()->getLocalUniqueId() < b.lock()->getLocalUniqueId();
623 };
624
625 // Sort all views
626
627 // Data
628 sort(fe_raw_ptr->getDataFieldEntsPtr()->begin(),
629 fe_raw_ptr->getDataFieldEntsPtr()->end(), uid_comp);
630 last_data_field_ents_view_size =
631 fe_raw_ptr->getDataFieldEntsPtr()->size();
632
633 // Row
634 if (!row_as_data) {
635 sort(fe_raw_ptr->getRowFieldEntsPtr()->begin(),
636 fe_raw_ptr->getRowFieldEntsPtr()->end(), uid_comp);
637 last_row_field_ents_view_size =
638 fe_raw_ptr->getRowFieldEntsPtr()->size();
639 }
640
641 // Column
642 if (!col_as_row) {
643 sort(fe_raw_ptr->getColFieldEntsPtr()->begin(),
644 fe_raw_ptr->getColFieldEntsPtr()->end(), uid_comp);
645 last_col_field_ents_view_size =
646 fe_raw_ptr->getColFieldEntsPtr()->size();
647 }
648 }
649 }
650
651 auto &dofs_by_ent_uid = dofsField.get<Unique_mi_tag>();
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<FiniteElement_name_mi_tag>();
669 for (auto &fe : finiteElements) {
670 auto miit = fe_ents.lower_bound(fe->getName());
671 auto hi_miit = fe_ents.upper_bound(fe->getName());
672 const auto count = std::distance(miit, hi_miit);
673 MOFEM_LOG("SYNC", Sev::inform)
674 << "Finite element " << fe->getName()
675 << " added. Nb. of elements added " << count;
676 MOFEM_LOG("SYNC", Sev::noisy) << *fe;
677
678 auto slg = MoFEM::LogManager::getLog("SYNC");
679 for (auto &field : fIelds) {
680 auto rec = slg.open_record(keywords::severity = Sev::verbose);
681 if (rec) {
682 logging::record_ostream strm(rec);
683 strm << "Field " << field->getName() << " on finite element: ";
684 if ((field->getId() & fe->getBitFieldIdRow()).any())
685 strm << "row ";
686 if ((field->getId() & fe->getBitFieldIdCol()).any())
687 strm << "columns ";
688 if ((field->getId() & fe->getBitFieldIdData()).any())
689 strm << "data";
690 strm.flush();
691 slg.push_record(boost::move(rec));
692 }
693 }
694 }
695
696 MOFEM_LOG_SYNCHRONISE(mofemComm);
697 }
698
699 *buildMoFEM |= 1 << 1;
701}
702
705 SETERRQ(mofemComm, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
707}
708
710 const Range *const ents_ptr,
711 int verb) {
713 if (verb == -1)
714 verb = verbose;
715
716 auto fe_miit = finiteElements.get<FiniteElement_name_mi_tag>().find(fe_name);
717 if (fe_miit == finiteElements.get<FiniteElement_name_mi_tag>().end())
718 SETERRQ1(mofemComm, MOFEM_NOT_FOUND, "Finite element <%s> not found",
719 fe_name.c_str());
720
721 CHKERR buildFiniteElements(*fe_miit, ents_ptr, verb);
722
723 if (verb >= VERBOSE) {
724 auto &fe_ents = entsFiniteElements.get<FiniteElement_name_mi_tag>();
725 auto miit = fe_ents.lower_bound((*fe_miit)->getName());
726 auto hi_miit = fe_ents.upper_bound((*fe_miit)->getName());
727 const auto count = std::distance(miit, hi_miit);
728 MOFEM_LOG("SYNC", Sev::inform) << "Finite element " << fe_name
729 << " added. Nb. of elements added " << count;
730 MOFEM_LOG_SYNCHRONISE(mofemComm);
731 }
732
733 *buildMoFEM |= 1 << 1;
735}
736
739 if (verb == DEFAULT_VERBOSITY)
740 verb = verbose;
741
742 if (!((*buildMoFEM) & BUILD_FIELD))
743 SETERRQ(mofemComm, MOFEM_NOT_FOUND, "field not build");
744 if (!((*buildMoFEM) & BUILD_FE))
745 SETERRQ(mofemComm, MOFEM_NOT_FOUND, "fe not build");
746 for (auto peit = ents.pair_begin(); peit != ents.pair_end(); ++peit) {
747 auto fit = entsFiniteElements.get<Ent_mi_tag>().lower_bound(peit->first);
748 auto hi_fit = entsFiniteElements.get<Ent_mi_tag>().upper_bound(peit->second);
749 for (; fit != hi_fit; ++fit) {
750 if ((*fit)->getBitFieldIdRow().none() &&
751 (*fit)->getBitFieldIdCol().none() &&
752 (*fit)->getBitFieldIdData().none())
753 continue;
754 int by = BYROW;
755 if ((*fit)->getBitFieldIdRow() != (*fit)->getBitFieldIdCol())
756 by |= BYCOL;
757 if ((*fit)->getBitFieldIdRow() != (*fit)->getBitFieldIdData())
758 by |= BYDATA;
760 auto hint = entFEAdjacencies.end();
761 for (auto e : *(*fit)->getRowFieldEntsPtr()) {
762 hint = entFEAdjacencies.emplace_hint(hint, e.lock(), *fit);
763 bool success = entFEAdjacencies.modify(hint, modify_row);
764 if (!success)
765 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
766 "modification unsuccessful");
767 }
768 if ((*fit)->getBitFieldIdRow() != (*fit)->getBitFieldIdCol()) {
769 int by = BYCOL;
770 if ((*fit)->getBitFieldIdCol() != (*fit)->getBitFieldIdData())
771 by |= BYDATA;
773 auto hint = entFEAdjacencies.end();
774 for (auto e : *(*fit)->getColFieldEntsPtr()) {
775 hint = entFEAdjacencies.emplace_hint(hint, e.lock(), *fit);
776 bool success = entFEAdjacencies.modify(hint, modify_col);
777 if (!success)
778 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
779 "modification unsuccessful");
780 }
781 }
782 if ((*fit)->getBitFieldIdRow() != (*fit)->getBitFieldIdData() ||
783 (*fit)->getBitFieldIdCol() != (*fit)->getBitFieldIdData()) {
785 BYDATA);
786 auto hint = entFEAdjacencies.end();
787 for (auto &e : (*fit)->getDataFieldEnts()) {
788 hint = entFEAdjacencies.emplace_hint(hint, e.lock(), *fit);
789 bool success = entFEAdjacencies.modify(hint, modify_data);
790 if (!success)
791 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
792 "modification unsuccessful");
793 }
794 }
795 }
796 }
797
798 if (verb >= VERBOSE) {
799 MOFEM_LOG("WORLD", Sev::inform)
800 << "Number of adjacencies " << entFEAdjacencies.size();
801 MOFEM_LOG_SYNCHRONISE(mofemComm)
802 }
803
804 *buildMoFEM |= 1 << 2;
806}
807
809 const BitRefLevel &mask, int verb) {
811 if (verb == -1)
812 verb = verbose;
813 Range ents;
814 CHKERR BitRefManager(*this).getEntitiesByRefLevel(bit, mask, ents);
815
816 CHKERR build_adjacencies(ents, verb);
817
819}
822 if (verb == -1)
823 verb = verbose;
824 CHKERR build_adjacencies(bit, BitRefLevel().set(), verb);
825
827}
828
829EntFiniteElementByName::iterator
830Core::get_fe_by_name_begin(const std::string &fe_name) const {
831 return entsFiniteElements.get<FiniteElement_name_mi_tag>().lower_bound(
832 fe_name);
833}
834EntFiniteElementByName::iterator
835Core::get_fe_by_name_end(const std::string &fe_name) const {
836 return entsFiniteElements.get<FiniteElement_name_mi_tag>().upper_bound(
837 fe_name);
838}
839
841 const std::string &name) const {
843 FiniteElement_multiIndex::index<FiniteElement_name_mi_tag>::type::iterator it;
844 it = finiteElements.get<FiniteElement_name_mi_tag>().find(name);
845 if (it == finiteElements.get<FiniteElement_name_mi_tag>().end()) {
846 SETERRQ1(mofemComm, 1, "finite element not found < %s >", name.c_str());
847 }
848 EntityHandle meshset = (*it)->getMeshset();
849
850 int num_entities;
851 CHKERR get_moab().get_number_entities_by_handle(meshset, num_entities);
852
853 if (entsFiniteElements.get<FiniteElement_name_mi_tag>().count(
854 (*it)->getName().c_str()) != (unsigned int)num_entities) {
855 SETERRQ1(mofemComm, 1,
856 "not equal number of entities in meshset and finite elements "
857 "multiindex < %s >",
858 (*it)->getName().c_str());
859 }
861}
862
865 FiniteElement_multiIndex::index<FiniteElement_name_mi_tag>::type::iterator it;
866 it = finiteElements.get<FiniteElement_name_mi_tag>().begin();
867 for (; it != finiteElements.get<FiniteElement_name_mi_tag>().end(); it++) {
868 EntityHandle meshset = (*it)->getMeshset();
869
870 int num_entities;
871 CHKERR get_moab().get_number_entities_by_handle(meshset, num_entities);
872
873 if (entsFiniteElements.get<FiniteElement_name_mi_tag>().count(
874 (*it)->getName().c_str()) != (unsigned int)num_entities) {
875 SETERRQ1(mofemComm, 1,
876 "not equal number of entities in meshset and finite elements "
877 "multiindex < %s >",
878 (*it)->getName().c_str());
879 }
880 }
882}
883} // namespace MoFEM
static Index< 'p', 3 > p
#define FECoreFunctionBegin
Definition: FECore.cpp:8
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:338
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_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:301
static LoggerType & getLog(const std::string channel)
Get logger by channel.
Definition: LogManager.cpp:370
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
uint128_t UId
Unique Id.
Definition: Types.hpp:31
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:1592
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:274
MoFEMErrorCode check_number_of_ents_in_ents_finite_element() const
check data consistency in entsFiniteElements
Definition: FECore.cpp:863
EntFiniteElementByName::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:835
EntityHandle get_finite_element_meshset(const BitFEId id) const
Definition: FECore.cpp:261
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:226
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:189
MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)
add finite element
Definition: FECore.cpp:28
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:429
std::string getBitFEIdName(const BitFEId id) const
Get field name.
Definition: FECore.cpp:253
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:334
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:297
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:167
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:317
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:123
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:393
bool check_finite_element(const std::string &name) const
Check if finite element is in database.
Definition: FECore.cpp:18
MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)
Build finite elements.
Definition: FECore.cpp:656
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:374
MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)
build adjacencies
Definition: FECore.cpp:737
BitFEId getBitFEId(const std::string &name) const
Get field Id.
Definition: FECore.cpp:243
EntFiniteElementByName::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:830
MoFEMErrorCode buildFiniteElements(const boost::shared_ptr< FiniteElement > &fe, const Range *ents_ptr=NULL, int verb=DEFAULT_VERBOSITY)
Definition: FECore.cpp:447
MoFEMErrorCode list_finite_elements() const
list finite elements in database
Definition: FECore.cpp:308
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:102
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:284
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:207
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:145
UId getLocalUniqueIdCalculate()
Get the Local Unique Id Calculate object.
Finite element definition.