v0.9.2
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 /* MoFEM is free software: you can redistribute it and/or modify it under
6  * the terms of the GNU Lesser General Public License as published by the
7  * Free Software Foundation, either version 3 of the License, or (at your
8  * option) any later version.
9  *
10  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
11  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
13  * License for more details.
14  *
15  * You should have received a copy of the GNU Lesser General Public
16  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>
17  */
18 
19 #define FECoreFunctionBegin \
20  MoFEMFunctionBegin; \
21  MOFEM_LOG_CHANNEL("WORLD"); \
22  MOFEM_LOG_CHANNEL("SYNC"); \
23  MOFEM_LOG_FUNCTION(); \
24  MOFEM_LOG_TAG("WORLD", "FECore"); \
25  MOFEM_LOG_TAG("SYNC", "FECore");
26 
27 namespace MoFEM {
28 
29 bool Core::check_finite_element(const std::string &name) const {
31  FeSetByName;
32  const FeSetByName &set = finiteElements.get<FiniteElement_name_mi_tag>();
33  FeSetByName::iterator miit = set.find(name);
34  if (miit == set.end())
35  return false;
36  return true;
37 }
38 
39 MoFEMErrorCode Core::add_finite_element(const std::string &fe_name,
40  enum MoFEMTypes bh, int verb) {
42  *buildMoFEM &= 1 << 0;
43  if (verb == -1) {
44  verb = verbose;
45  }
47  FiniteElements_by_name;
48  FiniteElements_by_name &finite_element_name_set =
50  FiniteElements_by_name::iterator it_fe =
51  finite_element_name_set.find(fe_name);
52  if (bh == MF_EXCL) {
53  if (it_fe != finite_element_name_set.end()) {
54  SETERRQ1(cOmm, MOFEM_NOT_FOUND, "this < %s > is there", fe_name.c_str());
55  }
56  } else {
57  if (it_fe != finite_element_name_set.end())
59  }
60  EntityHandle meshset;
61  CHKERR get_moab().create_meshset(MESHSET_SET, meshset);
62 
63  // id
64  BitFEId id = getFEShift();
65  CHKERR get_moab().tag_set_data(th_FEId, &meshset, 1, &id);
66 
67  // Add finite element meshset to partion meshset. In case of no elements
68  // on processor part, when mesh file is read, finite element meshset is
69  // prevented from deletion by moab reader.
70  auto add_meshset_to_partition = [&](auto meshset) {
72  const void *tag_vals[] = {&rAnk};
73  ParallelComm *pcomm = ParallelComm::get_pcomm(
74  &get_moab(), get_basic_entity_data_ptr()->pcommID);
75  Tag part_tag = pcomm->part_tag();
76  Range tagged_sets;
77  CHKERR get_moab().get_entities_by_type_and_tag(0, MBENTITYSET, &part_tag,
78  tag_vals, 1, tagged_sets,
79  moab::Interface::UNION);
80  for (auto s : tagged_sets)
81  CHKERR get_moab().add_entities(s, &meshset, 1);
83  };
84  CHKERR add_meshset_to_partition(meshset);
85 
86  // id name
87  void const *tag_data[] = {fe_name.c_str()};
88  int tag_sizes[1];
89  tag_sizes[0] = fe_name.size();
90  CHKERR get_moab().tag_set_by_ptr(th_FEName, &meshset, 1, tag_data, tag_sizes);
91 
92  // add FiniteElement
93  auto p = finiteElements.insert(
94  boost::shared_ptr<FiniteElement>(new FiniteElement(moab, meshset)));
95  if (!p.second)
96  SETERRQ(cOmm, MOFEM_OPERATION_UNSUCCESSFUL, "FiniteElement not inserted");
97 
98  if (verb > QUIET)
99  MOFEM_LOG("WORLD", Sev::inform) << "Add finite element " << fe_name;
100 
102 }
103 
106  const EntityType type,
107  ElementAdjacencyFunct function) {
109  *buildMoFEM &= 1 << 0;
111  FiniteElements_by_name;
112  FiniteElements_by_name &finite_element_name_set =
114  FiniteElements_by_name::iterator it_fe =
115  finite_element_name_set.find(fe_name);
116  if (it_fe == finite_element_name_set.end())
117  SETERRQ(cOmm, MOFEM_NOT_FOUND,
118  "This finite element is not defined (advise: check spelling)");
119  boost::shared_ptr<FiniteElement> fe;
120  fe = *it_fe;
121  fe->elementAdjacencyTable[type] = function;
123 }
124 
127  const std::string &name_data) {
129  *buildMoFEM &= 1 << 0;
131  FiniteElements_by_name;
132  FiniteElements_by_name &finite_element_name_set =
134  FiniteElements_by_name::iterator it_fe =
135  finite_element_name_set.find(fe_name);
136  if (it_fe == finite_element_name_set.end())
137  SETERRQ(cOmm, MOFEM_NOT_FOUND,
138  "This finite element is not defined (advise: check spelling)");
139  bool success = finite_element_name_set.modify(
140  it_fe, FiniteElement_change_bit_add(getBitFieldId(name_data)));
141  if (!success)
142  SETERRQ(cOmm, MOFEM_OPERATION_UNSUCCESSFUL, "modification unsuccessful");
144 }
145 
147 Core::modify_finite_element_add_field_row(const std::string &fe_name,
148  const std::string &name_row) {
150  *buildMoFEM &= 1 << 0;
152  FiniteElements_by_name;
153  FiniteElements_by_name &finite_element_name_set =
155  FiniteElements_by_name::iterator it_fe =
156  finite_element_name_set.find(fe_name);
157  if (it_fe == finite_element_name_set.end())
158  SETERRQ1(cOmm, MOFEM_NOT_FOUND, "this < %s > is not there",
159  fe_name.c_str());
160  bool success = finite_element_name_set.modify(
162  if (!success)
163  SETERRQ(cOmm, MOFEM_OPERATION_UNSUCCESSFUL, "modification unsuccessful");
165 }
166 
168 Core::modify_finite_element_add_field_col(const std::string &fe_name,
169  const std::string &name_col) {
171  *buildMoFEM &= 1 << 0;
173  FiniteElements_by_name;
174  FiniteElements_by_name &finite_element_name_set =
176  FiniteElements_by_name::iterator it_fe =
177  finite_element_name_set.find(fe_name);
178  if (it_fe == finite_element_name_set.end())
179  SETERRQ(cOmm, MOFEM_OPERATION_UNSUCCESSFUL, "this FiniteElement is there");
180  bool success = finite_element_name_set.modify(
182  if (!success)
183  SETERRQ(cOmm, MOFEM_OPERATION_UNSUCCESSFUL, "modification unsuccessful");
185 }
186 
189  const std::string &name_data) {
191  *buildMoFEM &= 1 << 0;
192  auto &finite_element_name_set =
194  auto it_fe = finite_element_name_set.find(fe_name);
195  if (it_fe == finite_element_name_set.end())
196  SETERRQ(cOmm, MOFEM_NOT_FOUND, "this FiniteElement is there");
197  bool success = finite_element_name_set.modify(
198  it_fe, FiniteElement_change_bit_off(getBitFieldId(name_data)));
199  if (!success)
200  SETERRQ(cOmm, MOFEM_OPERATION_UNSUCCESSFUL, "modification unsuccessful");
202 }
203 
205 Core::modify_finite_element_off_field_row(const std::string &fe_name,
206  const std::string &name_row) {
208  *buildMoFEM &= 1 << 0;
209  auto &finite_element_name_set =
211  auto it_fe = finite_element_name_set.find(fe_name);
212  if (it_fe == finite_element_name_set.end())
213  SETERRQ1(cOmm, MOFEM_NOT_FOUND, "this < %s > is not there",
214  fe_name.c_str());
215  bool success = finite_element_name_set.modify(
217  if (!success)
218  SETERRQ(cOmm, MOFEM_OPERATION_UNSUCCESSFUL, "modification unsuccessful");
220 }
221 
223 Core::modify_finite_element_off_field_col(const std::string &fe_name,
224  const std::string &name_col) {
226  *buildMoFEM &= 1 << 0;
227  auto &finite_element_name_set =
229  auto it_fe = finite_element_name_set.find(fe_name);
230  if (it_fe == finite_element_name_set.end())
231  SETERRQ(cOmm, MOFEM_NOT_FOUND, "this FiniteElement is there");
232  bool success = finite_element_name_set.modify(
234  if (!success)
235  SETERRQ(cOmm, MOFEM_OPERATION_UNSUCCESSFUL, "modification unsuccessful");
237 }
238 
239 BitFEId Core::getBitFEId(const std::string &name) const {
241  FiniteElements_by_name;
242  const FiniteElements_by_name &set =
244  FiniteElements_by_name::iterator miit = set.find(name);
245  if (miit == set.end())
247  ("finite element < " + name + " > not found (top tip: check spelling)")
248  .c_str());
249  return (*miit)->getId();
250 }
251 
252 std::string Core::getBitFEIdName(const BitFEId id) const {
254  finiteElements_by_id;
255  const finiteElements_by_id &set = finiteElements.get<BitFEId_mi_tag>();
256  finiteElements_by_id::iterator miit = set.find(id);
257  assert(miit != set.end());
258  return (*miit)->getName();
259 }
260 
263  finiteElements_by_id;
264  const finiteElements_by_id &set = finiteElements.get<BitFEId_mi_tag>();
265  finiteElements_by_id::iterator miit = set.find(id);
266  if (miit == set.end())
267  THROW_MESSAGE("finite element not found");
268  return (*miit)->meshset;
269 }
270 
271 EntityHandle Core::get_finite_element_meshset(const std::string &name) const {
273 }
274 
277  Range &ents) const {
278 
280 
282  CHKERR get_moab().get_entities_by_dimension(meshset, dim, ents, true);
284 }
285 
287  EntityType type,
288  Range &ents) const {
289 
291 
293  CHKERR get_moab().get_entities_by_type(meshset, type, ents, true);
294 
296 }
297 
300  Range &ents) const {
301 
303 
305  CHKERR get_moab().get_entities_by_handle(meshset, ents, true);
306 
308 }
309 
313  finiteElements_by_id;
314  const finiteElements_by_id &BitFEId_set =
316  finiteElements_by_id::iterator miit = BitFEId_set.begin();
317  for (; miit != BitFEId_set.end(); miit++) {
318  std::ostringstream ss;
319  ss << *miit << std::endl;
320  PetscSynchronizedPrintf(cOmm, ss.str().c_str());
321  }
322  PetscSynchronizedFlush(cOmm, PETSC_STDOUT);
324 }
325 
327  const EntityHandle meshset, const EntityType type, const std::string &name,
328  const bool recursive) {
329  *buildMoFEM &= 1 << 0;
330  EntityHandle idm = no_handle;
332 
334  Range ents;
335  CHKERR get_moab().get_entities_by_type(meshset, type, ents, recursive);
336  CHKERR getInterface<BitRefManager>()->setElementsBitRefLevel(ents);
337  CHKERR get_moab().add_entities(idm, ents);
338 
340 }
341 
344  const int dim, const std::string &name,
345  const bool recursive) {
346  EntityHandle idm = no_handle;
347  *buildMoFEM &= 1 << 0;
350  Range ents;
351  CHKERR get_moab().get_entities_by_dimension(meshset, dim, ents, recursive);
352  CHKERR getInterface<BitRefManager>()->setElementsBitRefLevel(ents);
353  CHKERR get_moab().add_entities(idm, ents);
355 }
356 
358  const Range &ents, const EntityType type, const std::string &name) {
359  EntityHandle idm = no_handle;
360  *buildMoFEM &= 1 << 0;
363  CHKERR getInterface<BitRefManager>()->setElementsBitRefLevel(
364  ents.subset_by_type(type));
365  CHKERR get_moab().add_entities(idm, ents.subset_by_type(type));
367 } // namespace MoFEM
368 
370 Core::add_ents_to_finite_element_by_dim(const Range &ents, const int dim,
371  const std::string &name) {
372  EntityHandle idm = no_handle;
373  *buildMoFEM &= 1 << 0;
376  CHKERR getInterface<BitRefManager>()->setElementsBitRefLevel(
377  ents.subset_by_dimension(dim));
378  CHKERR get_moab().add_entities(idm, ents.subset_by_dimension(dim));
380 }
381 
383  const EntityHandle meshset, const std::string &name, const bool recursive) {
385  CHKERR add_ents_to_finite_element_by_type(meshset, MBEDGE, name, recursive);
387 }
390  const std::string &name) {
392  CHKERR add_ents_to_finite_element_by_type(edges, MBEDGE, name);
394 }
397  const std::string &name) {
399  CHKERR add_ents_to_finite_element_by_type(vert, MBVERTEX, name);
401 }
404  const std::string &name) {
406  CHKERR add_ents_to_finite_element_by_type(tris, MBTRI, name);
408 }
410  const EntityHandle meshset, const std::string &name, const bool recursive) {
412  CHKERR add_ents_to_finite_element_by_type(meshset, MBTRI, name, recursive);
414 }
417  const std::string &name) {
419  CHKERR add_ents_to_finite_element_by_type(tets, MBTET, name);
421 }
423  const EntityHandle meshset, const std::string &name, const bool recursive) {
425  CHKERR add_ents_to_finite_element_by_type(meshset, MBTET, name, recursive);
427 }
430  const std::string &name) {
432  CHKERR add_ents_to_finite_element_by_type(prims, MBPRISM, name);
434 }
436  const EntityHandle meshset, const std::string &name, const bool recursive) {
438  CHKERR add_ents_to_finite_element_by_type(meshset, MBPRISM, name, recursive);
440 }
441 
444  const std::string &name,
445  EntityType type, int verb) {
448  type, verb);
449 
451 }
452 
454  const BitRefLevel &bit, const BitRefLevel &mask, const std::string &name,
455  EntityType type, int verb) {
457  CHKERR add_ents_to_finite_element_by_bit_ref(bit, mask, name, type, verb);
458 
460 }
461 
463  const BitRefLevel &bit, const BitRefLevel &mask, const std::string &name,
464  EntityType type, int verb) {
466 
467  if (verb == -1)
468  verb = verbose;
469  *buildMoFEM &= 1 << 0;
470  const BitFEId id = getBitFEId(name);
472 
473  auto &ref_MoFEMFiniteElement = refinedFiniteElements.get<EntType_mi_tag>();
474  auto miit = ref_MoFEMFiniteElement.lower_bound(type);
475  auto hi_miit = ref_MoFEMFiniteElement.upper_bound(type);
476 
477  int nb_add_fes = 0;
478  for (; miit != hi_miit; miit++) {
479  BitRefLevel bit2 = miit->get()->getBitRefLevel();
480  if ((bit2 & mask) != bit2)
481  continue;
482  if ((bit2 & bit).any()) {
483  EntityHandle ent = miit->get()->getRefEnt();
484  CHKERR get_moab().add_entities(idm, &ent, 1);
485  nb_add_fes++;
486  }
487  }
488 
489  MOFEM_LOG("SYNC", Sev::inform)
490  << "Finite element " << name << " added. Nb. of elements added "
491  << nb_add_fes << " out of " << std::distance(miit, hi_miit);
492 
494 
496 }
497 
499  const EntityHandle meshset, const std::string &name, const bool recursive) {
501  *buildMoFEM &= 1 << 0;
502  const BitFEId id = getBitFEId(name);
504  if (recursive == false) {
505  CHKERR get_moab().add_entities(idm, &meshset, 1);
506  } else {
507  Range meshsets;
508  CHKERR get_moab().get_entities_by_type(meshset, MBENTITYSET, meshsets,
509  false);
510  CHKERR get_moab().add_entities(idm, meshsets);
511  }
513 }
514 
515 template <int I> struct BuildFiniteElements {
516 
517  template <typename T1, typename T2>
518  static inline void addToData(T1 &range_dit, T2 &fe_vec) {
519  static_assert(I == DATA, "t should be set to DATA");
520 
521  for (auto dit = range_dit.first; dit != range_dit.second; ++dit) {
522  const EntityHandle dof_ent = dit->get()->getEnt();
523  // Fill array
524  for (auto fe_it : fe_vec) {
525  if (auto fe_ptr = fe_it.lock()) {
526  // Add FEDofEntity, first create dofs, one by one, note that memory
527  // is already reserved. Then create shared pointers and finally add
528  // th_FEName to element multi-index
529  // There are data dofs on this element
530  auto &side_number_ptr = fe_ptr->getSideNumberPtr(dof_ent);
531  fe_ptr->getDofsSequence().lock()->emplace_back(side_number_ptr, *dit);
532  }
533  }
534  }
535  }
536 
537  template <typename T> static inline void emplaceHint(T &fe_vec) {
538  static_assert(I == DATA, "t should be set to DATA");
539 
540  // Add to data in FE
541  for (auto fe_it : fe_vec) {
542  if (auto fe_ptr = fe_it.lock()) {
543  // It is a but unsafe, since if one mess up something in
544  // buildFiniteElements, weak_ptr will not give pointer
545  auto data_dofs_array_vec = fe_ptr->getDofsSequence().lock();
546  // Create shared pointers vector
547  auto hint = fe_ptr->data_dofs->end();
548  for (auto &vit : *data_dofs_array_vec)
549  hint =
550  fe_ptr->data_dofs->emplace_hint(hint, data_dofs_array_vec, &vit);
551  }
552  }
553  }
554 };
555 
557 Core::buildFiniteElements(const boost::shared_ptr<FiniteElement> &fe,
558  const Range *ents_ptr, int verb) {
560  if (verb == DEFAULT_VERBOSITY)
561  verb = verbose;
562 
563  if (verb > QUIET)
564  MOFEM_LOG("SYNC", Sev::verbose)
565  << "Build Finite Elements " << fe->getName();
566 
567  auto &fields_by_id = fIelds.get<BitFieldId_mi_tag>();
568 
569  // Get id of mofem fields for row, col and data
570  enum IntLoop { ROW = 0, COL, DATA, LAST };
571  std::array<BitFieldId, LAST> fe_fields = {fe.get()->getBitFieldIdRow(),
572  fe.get()->getBitFieldIdCol(),
573  fe.get()->getBitFieldIdData()};
574 
575  // Get finite element meshset
576  EntityHandle meshset = get_finite_element_meshset(fe.get()->getId());
577 
578  // Get entities from finite element meshset // if meshset
579  Range fe_ents;
580  CHKERR get_moab().get_entities_by_handle(meshset, fe_ents, false);
581 
582  if (ents_ptr)
583  fe_ents = intersect(fe_ents, *ents_ptr);
584 
585  // Map entity uid to pointers
586  typedef std::vector<boost::weak_ptr<EntFiniteElement>> VecOfWeakFEPtrs;
587  typedef std::map<const UId *, VecOfWeakFEPtrs> MapEntUIdAndVecOfWeakFEPtrs;
588  MapEntUIdAndVecOfWeakFEPtrs ent_uid_and_fe_vec;
589  std::map<EntityHandle, boost::shared_ptr<std::vector<FEDofEntity>>>
590  data_dofs_array;
591  VecOfWeakFEPtrs processed_fes;
592  processed_fes.reserve(fe_ents.size());
593 
594  int last_row_field_ents_view_size = 0;
595  int last_col_field_ents_view_size = 0;
596 
597  // View of field entities on element
598  FieldEntity_vector_view data_field_ents_view;
599 
600  // Loop meshset finite element ents and add finite elements
601  for (Range::const_pair_iterator peit = fe_ents.const_pair_begin();
602  peit != fe_ents.const_pair_end(); peit++) {
603 
604  EntityHandle first = peit->first;
605  EntityHandle second = peit->second;
606 
607  // Find range of ref entities that is sequence
608  // note: iterator is a wrapper
609  // check if is in refinedFiniteElements database
610  auto ref_fe_miit =
611  refinedFiniteElements.get<Ent_mi_tag>().lower_bound(first);
612  if (ref_fe_miit == refinedFiniteElements.get<Ent_mi_tag>().end()) {
613  std::ostringstream ss;
614  ss << "refinedFiniteElements not in database ent = " << first;
615  ss << " type " << get_moab().type_from_handle(first);
616  ss << " " << *fe;
617  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, ss.str().c_str());
618  }
619  auto hi_ref_fe_miit =
620  refinedFiniteElements.get<Ent_mi_tag>().upper_bound(second);
621 
622  EntFiniteElement_multiIndex::iterator hint_p = entsFiniteElements.end();
623  for (; ref_fe_miit != hi_ref_fe_miit; ref_fe_miit++) {
624 
625  // Add finite element to database
626  hint_p = entsFiniteElements.emplace_hint(
627  hint_p, boost::make_shared<EntFiniteElement>(*ref_fe_miit, fe));
628  processed_fes.emplace_back(*hint_p);
629  auto fe_raw_ptr = hint_p->get();
630 
631  // Allocate space for etities view
632  data_field_ents_view.clear();
633  fe_raw_ptr->row_field_ents_view->reserve(last_row_field_ents_view_size);
634  // Create shared pointer for entities view
635  if (fe_fields[ROW] == fe_fields[COL]) {
636  fe_raw_ptr->col_field_ents_view = fe_raw_ptr->row_field_ents_view;
637  } else {
638  // row and columns are diffent
639  if (fe_raw_ptr->col_field_ents_view == fe_raw_ptr->row_field_ents_view)
640  fe_raw_ptr->col_field_ents_view =
641  boost::make_shared<FieldEntity_vector_view>();
642  fe_raw_ptr->col_field_ents_view->reserve(last_col_field_ents_view_size);
643  }
644 
645  int nb_dofs_on_data = 0;
646 
647  // Iterate over all field and check which one is on the element
648  for (unsigned int ii = 0; ii != BitFieldId().size(); ++ii) {
649 
650  // Common field id for ROW, COL and DATA
651  BitFieldId id_common = 0;
652  // Check if the field (ii) is added to finite element
653  for (int ss = 0; ss < LAST; ss++) {
654  id_common |= fe_fields[ss] & BitFieldId().set(ii);
655  }
656  if (id_common.none())
657  continue;
658 
659  // Find in database data associated with the field (ii)
660  const BitFieldId field_id = BitFieldId().set(ii);
661  auto miit = fields_by_id.find(field_id);
662  if (miit == fields_by_id.end()) {
663  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
664  "Data inconsistency");
665  }
666 
667  // Loop over adjacencies of element and find field entities on those
668  // adjacencies, that create hash map map_uid_fe which is used later
669  const std::string field_name = miit->get()->getName();
670  const bool add_to_data = (field_id & fe_fields[DATA]).any();
671  const bool add_to_row = (field_id & fe_fields[ROW]).any();
672  const bool add_to_col = (field_id & fe_fields[COL]).any();
673 
674  // Entities adjacent to entities
675  Range adj_ents;
676 
677  // Resolve entities on element, those entities are used to build tag
678  // with dof uids on finite element tag
679  CHKERR fe_raw_ptr->getElementAdjacency(*miit, adj_ents);
680 
681  for (Range::pair_iterator p_eit = adj_ents.pair_begin();
682  p_eit != adj_ents.pair_end(); ++p_eit) {
683 
684  const EntityHandle first = p_eit->first;
685  const EntityHandle second = p_eit->second;
686 
687  typedef FieldEntity_multiIndex::index<
688  Composite_Name_And_Ent_mi_tag>::type FieldEntityByComposite;
689  auto &field_ents_by_name_and_ent =
691  FieldEntityByComposite::iterator meit;
692 
693  // If one entity in the pair search for one, otherwise search for
694  // range
695  if (first == second)
696  meit = field_ents_by_name_and_ent.find(
697  boost::make_tuple(field_name, first));
698  else
699  meit = field_ents_by_name_and_ent.lower_bound(
700  boost::make_tuple(field_name, first));
701 
702  if (meit != field_ents_by_name_and_ent.end()) {
703 
704  decltype(meit) hi_meit;
705 
706  if (first == second) {
707  hi_meit = meit;
708  ++hi_meit;
709  } else
710  hi_meit = field_ents_by_name_and_ent.upper_bound(
711  boost::make_tuple(field_name, second));
712 
713  // Add to view and create list of finite elements with this dof UId
714  for (; meit != hi_meit; ++meit) {
715  // Add entity to map with key entity uids pointers and data
716  // finite elements weak ptrs. I using pointers to uids instead
717  // uids because this is faster.
718  const UId *uid_ptr = &(meit->get()->getGlobalUniqueId());
719  auto &fe_vec = ent_uid_and_fe_vec[uid_ptr];
720  // get number of dofs on entities to pre-allocate memory for
721  // element
722  const int nb_dofs_on_ent = meit->get()->getNbDofsOnEnt();
723  if (add_to_data) {
724  nb_dofs_on_data += nb_dofs_on_ent;
725  data_field_ents_view.emplace_back(*meit);
726  }
727  if (add_to_row) {
728  fe_raw_ptr->row_field_ents_view->emplace_back(*meit);
729  }
730  if (add_to_col) {
731  if (fe_raw_ptr->col_field_ents_view !=
732  fe_raw_ptr->row_field_ents_view)
733  fe_raw_ptr->col_field_ents_view->emplace_back(*meit);
734  }
735  // add finite element to processed list
736  fe_vec.emplace_back(*hint_p);
737  }
738  }
739  }
740  }
741 
742  // Sort field ents by uid
743  auto uid_comp = [](const auto &a, const auto &b) {
744  return a.lock()->getGlobalUniqueId() < b.lock()->getGlobalUniqueId();
745  };
746 
747  // Sort all views
748 
749  // Data
750  sort(data_field_ents_view.begin(), data_field_ents_view.end(), uid_comp);
751  for (auto e : data_field_ents_view)
752  fe_raw_ptr->data_field_ents_view->emplace_back(e);
753 
754  // Row
755  sort(fe_raw_ptr->row_field_ents_view->begin(),
756  fe_raw_ptr->row_field_ents_view->end(), uid_comp);
757  last_row_field_ents_view_size = fe_raw_ptr->row_field_ents_view->size();
758 
759  // Column
760  if (fe_raw_ptr->col_field_ents_view != fe_raw_ptr->row_field_ents_view) {
761  sort(fe_raw_ptr->col_field_ents_view->begin(),
762  fe_raw_ptr->col_field_ents_view->end(), uid_comp);
763  last_col_field_ents_view_size = fe_raw_ptr->col_field_ents_view->size();
764  }
765 
766  // Clear finite element data structures
767  fe_raw_ptr->data_dofs->clear();
768 
769  // Reserve memory for data FE Dofs
770  auto data_dofs_array_vec = boost::make_shared<std::vector<FEDofEntity>>();
771  data_dofs_array[fe_raw_ptr->getEnt()] = data_dofs_array_vec;
772  data_dofs_array_vec->reserve(nb_dofs_on_data);
773 
774  fe_raw_ptr->getDofsSequence() = data_dofs_array_vec;
775  }
776  }
777 
778  auto &dofs_by_ent_uid = dofsField.get<Unique_Ent_mi_tag>();
779 
780  // Loop over hash map, which has all entities on given elemnts
781  boost::shared_ptr<SideNumber> side_number_ptr;
782  for (auto &mit : ent_uid_and_fe_vec) {
783  auto range_dit = dofs_by_ent_uid.equal_range(*mit.first);
784  if (range_dit.first != range_dit.second) {
785  const BitFieldId field_id = range_dit.first->get()->getId();
786  if ((field_id & fe_fields[DATA]).any())
787  BuildFiniteElements<DATA>::addToData(range_dit, mit.second);
788  }
789  }
790 
792 
794 }
795 
798 
799  if (verb == DEFAULT_VERBOSITY)
800  verb = verbose;
801 
802  // loop Finite Elements
803  for (auto &fe : finiteElements)
804  CHKERR buildFiniteElements(fe, NULL, verb);
805 
806  if (verb > QUIET) {
807 
808  auto &fe_ents = entsFiniteElements.get<BitFEId_mi_tag>();
809  for (auto &fe : finiteElements) {
810  auto miit = fe_ents.lower_bound(fe->getId());
811  auto hi_miit = fe_ents.upper_bound(fe->getId());
812  const auto count = std::distance(miit, hi_miit);
813  MOFEM_LOG("SYNC", Sev::inform)
814  << "Finite element " << fe->getName()
815  << " added. Nb. of elements added " << count;
816  MOFEM_LOG("SYNC", Sev::noisy) << *fe;
817 
818  auto slg = MoFEM::LogManager::getLog("SYNC");
819  for (auto &field : fIelds) {
820  auto rec = slg.open_record(keywords::severity = Sev::verbose);
821  if (rec) {
822  logging::record_ostream strm(rec);
823  strm << "Field " << field->getName() << " on finite element: ";
824  if ((field->getId() & fe->getBitFieldIdRow()).any())
825  strm << "row ";
826  if ((field->getId() & fe->getBitFieldIdCol()).any())
827  strm << "columns ";
828  if ((field->getId() & fe->getBitFieldIdData()).any())
829  strm << "data";
830  strm.flush();
831  slg.push_record(boost::move(rec));
832  }
833  }
834  }
835 
837  }
838 
839  *buildMoFEM |= 1 << 1;
841 }
842 
845  SETERRQ(cOmm, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
847 }
848 
850  const Range *const ents_ptr,
851  int verb) {
853  if (verb == -1)
854  verb = verbose;
855 
856  auto fe_miit = finiteElements.get<FiniteElement_name_mi_tag>().find(fe_name);
857  if (fe_miit == finiteElements.get<FiniteElement_name_mi_tag>().end())
858  SETERRQ1(cOmm, MOFEM_NOT_FOUND, "Finite element <%s> not found",
859  fe_name.c_str());
860 
861  CHKERR buildFiniteElements(*fe_miit, ents_ptr, verb);
862 
863  if (verb >= VERBOSE) {
864  auto &fe_ents = entsFiniteElements.get<BitFEId_mi_tag>();
865  auto miit = fe_ents.lower_bound((*fe_miit)->getId());
866  auto hi_miit = fe_ents.upper_bound((*fe_miit)->getId());
867  const auto count = std::distance(miit, hi_miit);
868  MOFEM_LOG("SYNC", Sev::inform) << "Finite element " << fe_name
869  << " added. Nb. of elements added " << count;
871  }
872 
873  *buildMoFEM |= 1 << 1;
875 }
876 
877 MoFEMErrorCode Core::build_adjacencies(const Range &ents, int verb) {
879  if (verb == DEFAULT_VERBOSITY)
880  verb = verbose;
881 
882  if (!((*buildMoFEM) & BUILD_FIELD))
883  SETERRQ(cOmm, MOFEM_NOT_FOUND, "field not build");
884  if (!((*buildMoFEM) & BUILD_FE))
885  SETERRQ(cOmm, MOFEM_NOT_FOUND, "fe not build");
886  for (Range::const_pair_iterator peit = ents.pair_begin();
887  peit != ents.pair_end(); ++peit) {
888  EntFiniteElement_multiIndex::index<Ent_mi_tag>::type::iterator fit, hi_fit;
889  fit = entsFiniteElements.get<Ent_mi_tag>().lower_bound(peit->first);
890  hi_fit = entsFiniteElements.get<Ent_mi_tag>().upper_bound(peit->second);
891  for (; fit != hi_fit; ++fit) {
892  if ((*fit)->getBitFieldIdRow().none() &&
893  (*fit)->getBitFieldIdCol().none() &&
894  (*fit)->getBitFieldIdData().none())
895  continue;
896  int by = BYROW;
897  if ((*fit)->getBitFieldIdRow() != (*fit)->getBitFieldIdCol())
898  by |= BYCOL;
899  if ((*fit)->getBitFieldIdRow() != (*fit)->getBitFieldIdData())
900  by |= BYDATA;
902  auto hint = entFEAdjacencies.end();
903  for (auto e : *(*fit)->row_field_ents_view) {
904  hint = entFEAdjacencies.emplace_hint(hint, e.lock(), *fit);
905  bool success = entFEAdjacencies.modify(hint, modify_row);
906  if (!success)
907  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
908  "modification unsuccessful");
909  }
910  if ((*fit)->getBitFieldIdRow() != (*fit)->getBitFieldIdCol()) {
911  int by = BYCOL;
912  if ((*fit)->getBitFieldIdCol() != (*fit)->getBitFieldIdData())
913  by |= BYDATA;
915  auto hint = entFEAdjacencies.end();
916  for (auto e : *(*fit)->col_field_ents_view) {
917  hint = entFEAdjacencies.emplace_hint(hint, e.lock(), *fit);
918  bool success = entFEAdjacencies.modify(hint, modify_col);
919  if (!success)
920  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
921  "modification unsuccessful");
922  }
923  }
924  if ((*fit)->getBitFieldIdRow() != (*fit)->getBitFieldIdData() ||
925  (*fit)->getBitFieldIdCol() != (*fit)->getBitFieldIdData()) {
927  BYDATA);
928  auto hint = entFEAdjacencies.end();
929  for (auto &e : *(*fit)->data_field_ents_view) {
930  hint = entFEAdjacencies.emplace_hint(hint, e, *fit);
931  bool success = entFEAdjacencies.modify(hint, modify_data);
932  if (!success)
933  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
934  "modification unsuccessful");
935  }
936  }
937  }
938  }
939 
940  if (verb >= VERBOSE) {
941  MOFEM_LOG("WORLD", Sev::inform)
942  << "Number of adjacencies " << entFEAdjacencies.size();
944  }
945 
946  *buildMoFEM |= 1 << 2;
948 }
949 
951  const BitRefLevel &mask, int verb) {
953  if (verb == -1)
954  verb = verbose;
955  Range ents;
956  CHKERR BitRefManager(*this).getEntitiesByRefLevel(bit, mask, ents);
957 
958  CHKERR build_adjacencies(ents, verb);
959 
961 }
964  if (verb == -1)
965  verb = verbose;
966  CHKERR build_adjacencies(bit, BitRefLevel().set(), verb);
967 
969 }
970 
971 EntFiniteElementByName::iterator
972 Core::get_fe_by_name_begin(const std::string &fe_name) const {
973  return entsFiniteElements.get<FiniteElement_name_mi_tag>().lower_bound(
974  fe_name);
975 }
976 EntFiniteElementByName::iterator
977 Core::get_fe_by_name_end(const std::string &fe_name) const {
978  return entsFiniteElements.get<FiniteElement_name_mi_tag>().upper_bound(
979  fe_name);
980 }
981 
983  const std::string &name) const {
985  FiniteElement_multiIndex::index<FiniteElement_name_mi_tag>::type::iterator it;
986  it = finiteElements.get<FiniteElement_name_mi_tag>().find(name);
987  if (it == finiteElements.get<FiniteElement_name_mi_tag>().end()) {
988  SETERRQ1(cOmm, 1, "finite element not found < %s >", name.c_str());
989  }
990  EntityHandle meshset = (*it)->getMeshset();
991 
992  int num_entities;
993  CHKERR get_moab().get_number_entities_by_handle(meshset, num_entities);
994 
996  (*it)->getName().c_str()) != (unsigned int)num_entities) {
997  SETERRQ1(cOmm, 1,
998  "not equal number of entities in meshset and finite elements "
999  "multiindex < %s >",
1000  (*it)->getName().c_str());
1001  }
1003 }
1004 
1007  FiniteElement_multiIndex::index<FiniteElement_name_mi_tag>::type::iterator it;
1008  it = finiteElements.get<FiniteElement_name_mi_tag>().begin();
1009  for (; it != finiteElements.get<FiniteElement_name_mi_tag>().end(); it++) {
1010  EntityHandle meshset = (*it)->getMeshset();
1011 
1012  int num_entities;
1013  CHKERR get_moab().get_number_entities_by_handle(meshset, num_entities);
1014 
1016  (*it)->getName().c_str()) != (unsigned int)num_entities) {
1017  SETERRQ1(cOmm, 1,
1018  "not equal number of entities in meshset and finite elements "
1019  "multiindex < %s >",
1020  (*it)->getName().c_str());
1021  }
1022  }
1024 }
1025 } // namespace MoFEM
DofEntity_multiIndex dofsField
dofs on fields
Definition: Core.hpp:249
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:126
EntityHandle get_finite_element_meshset(const BitFEId id) const
Definition: FECore.cpp:261
MoFEMErrorCode check_number_of_ents_in_ents_finite_element() const
check data consistency in entsFiniteElements
Definition: FECore.cpp:1005
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:498
MPI_Comm cOmm
MoFEM communicator.
Definition: Core.hpp:856
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:343
FieldEntity_multiIndex entsFields
entities on fields
Definition: Core.hpp:248
DEPRECATED MoFEMErrorCode add_ents_to_finite_element_by_VERTICEs(const Range &vert, const std::string &name)
add VERTICES entities from range to finite element database given by name
Definition: FECore.cpp:396
DEPRECATED MoFEMErrorCode add_ents_to_finite_element_by_EDGEs(const Range &vert, const std::string &name)
add EDGES entities from range to finite element database given by name
Definition: FECore.cpp:389
moab::Interface & get_moab()
Definition: Core.hpp:266
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:507
std::vector< boost::weak_ptr< FieldEntity > > FieldEntity_vector_view
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:303
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
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:276
MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)
build adjacencies
Definition: FECore.cpp:877
const EntityHandle no_handle
No entity handle is indicated by zero handle, i.e. root meshset.
Definition: Common.hpp:27
FTensor::Index< 'I', 3 > I
Definition: PlasticOps.hpp:70
MoFEMErrorCode list_finite_elements() const
list finite elements in database
Definition: FECore.cpp:310
MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)
Build finite elementsBuild finite element data structures. Have to be run before problem and adjacenc...
Definition: FECore.cpp:796
#define THROW_MESSAGE(a)
Throw MoFEM exception.
Definition: definitions.h:626
Tag th_FEName
Definition: Core.hpp:206
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:514
std::reference_wrapper< moab::Interface > moab
moab database
Definition: Core.hpp:265
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:462
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:188
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
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:286
DEPRECATED MoFEMErrorCode add_ents_to_finite_element_by_PRISMs(const Range &prims, const BitFEId id)
DEPRECATED MoFEMErrorCode add_ents_to_finite_element_by_TRIs(const Range &tris, const std::string &name)
add TRI entities from range to finite element database given by name
Definition: FECore.cpp:403
BitFEId getBitFEId(const std::string &name) const
Get field Id.
Definition: FECore.cpp:239
std::bitset< BITFEID_SIZE > BitFEId
Finite element Id.
Definition: Types.hpp:53
const int dim
std::string getBitFEIdName(const BitFEId id) const
Get field name.
Definition: FECore.cpp:252
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:977
int rAnk
MOFEM communicator rank.
Definition: Core.hpp:860
static void emplaceHint(T &fe_vec)
Definition: FECore.cpp:537
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:223
MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)
add finite element
Definition: FECore.cpp:39
BitFEId getFEShift()
Return unique finite element Id.
Definition: Core.cpp:260
bool check_finite_element(const std::string &name) const
Check if finite element is in database.
Definition: FECore.cpp:29
Finite element definition.
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:205
int * buildMoFEM
keeps flags/semaphores for different stages
Definition: Core.hpp:906
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
FieldEntityEntFiniteElementAdjacencyMap_multiIndex entFEAdjacencies
adjacencies of elements to dofs
Definition: Core.hpp:255
Managing BitRefLevels.
RefElement_multiIndex refinedFiniteElements
refined elements
Definition: Core.hpp:245
FiniteElement_multiIndex finiteElements
finite elements
Definition: Core.hpp:251
#define FECoreFunctionBegin
Definition: FECore.cpp:19
static LoggerType & getLog(const std::string channel)
Get logger by channel.
Definition: LogManager.cpp:369
std::bitset< BITFIELDID_SIZE > BitFieldId
Field Id.
Definition: Types.hpp:52
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:168
#define CHKERR
Inline error check.
Definition: definitions.h:602
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:326
Field_multiIndex fIelds
fields
Definition: Core.hpp:247
#define MOFEM_LOG_SYNCHORMISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:340
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:50
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
BitFieldId getBitFieldId(const std::string &name) const
Definition: FieldCore.cpp:28
static void addToData(T1 &range_dit, T2 &fe_vec)
Definition: FECore.cpp:518
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:147
DEPRECATED MoFEMErrorCode add_ents_to_finite_element_EntType_by_bit_ref(const BitRefLevel &bit, const std::string &name, EntityType type, int verb=DEFAULT_VERBOSITY)
add TET elements from given refinement level to finite element database given by name
Definition: FECore.cpp:443
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:105
int verbose
Verbosity level.
Definition: Core.hpp:895
boost::function< MoFEMErrorCode(Interface &moab, const Field &field, const EntFiniteElement &fe, Range &adjacency)> ElementAdjacencyFunct
user adjacency function
Tag th_FEId
Definition: Core.hpp:206
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413
MoFEMErrorCode buildFiniteElements(const boost::shared_ptr< FiniteElement > &fe, const Range *ents_ptr=NULL, int verb=DEFAULT_VERBOSITY)
Definition: FECore.cpp:557
EntFiniteElement_multiIndex entsFiniteElements
finite element entities
Definition: Core.hpp:252
boost::shared_ptr< BasicEntityData > & get_basic_entity_data_ptr()
Get pointer to basic entity data.
Definition: Core.hpp:234
uint128_t UId
Unique Id.
Definition: Types.hpp:41
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:299
MoFEMTypes
Those types control how functions respond on arguments, f.e. error handling.
Definition: definitions.h:189
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:972
DEPRECATED MoFEMErrorCode add_ents_to_finite_element_by_TETs(const Range &tets, const std::string &name)
add TET entities from range to finite element database given by name
Definition: FECore.cpp:416