v0.9.1
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 namespace MoFEM {
20 
24  *fe_ptr = &finiteElements;
26 }
27 
29  const EntFiniteElement_multiIndex **fe_ent_ptr) const {
31  *fe_ent_ptr = &entsFiniteElements;
33 }
34 
35 bool Core::check_finite_element(const std::string &name) const {
37  FeSetByName;
38  const FeSetByName &set = finiteElements.get<FiniteElement_name_mi_tag>();
39  FeSetByName::iterator miit = set.find(name);
40  if (miit == set.end())
41  return false;
42  return true;
43 }
44 
45 MoFEMErrorCode Core::add_finite_element(const std::string &fe_name,
46  enum MoFEMTypes bh, int verb) {
48  *buildMoFEM &= 1 << 0;
49  if (verb == -1) {
50  verb = verbose;
51  }
53  FiniteElements_by_name;
54  FiniteElements_by_name &finite_element_name_set =
56  FiniteElements_by_name::iterator it_fe =
57  finite_element_name_set.find(fe_name);
58  if (bh == MF_EXCL) {
59  if (it_fe != finite_element_name_set.end()) {
60  SETERRQ1(cOmm, MOFEM_NOT_FOUND, "this < %s > is there", fe_name.c_str());
61  }
62  } else {
63  if (it_fe != finite_element_name_set.end())
65  }
66  EntityHandle meshset;
67  CHKERR get_moab().create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER, meshset);
68 
69  // id
70  BitFEId id = getFEShift();
71  CHKERR get_moab().tag_set_data(th_FEId, &meshset, 1, &id);
72 
73  // Add finite element meshset to partion meshset. In case of no elements
74  // on processor part, when mesh file is read, finite element meshset is
75  // prevented from deletion by moab reader.
76  auto add_meshset_to_partition = [&](auto meshset) {
78  const void *tag_vals[] = {&rAnk};
79  ParallelComm *pcomm = ParallelComm::get_pcomm(
80  &get_moab(), get_basic_entity_data_ptr()->pcommID);
81  Tag part_tag = pcomm->part_tag();
82  Range tagged_sets;
83  CHKERR get_moab().get_entities_by_type_and_tag(0, MBENTITYSET, &part_tag,
84  tag_vals, 1, tagged_sets,
85  moab::Interface::UNION);
86  for (auto s : tagged_sets)
87  CHKERR get_moab().add_entities(s, &meshset, 1);
89  };
90  CHKERR add_meshset_to_partition(meshset);
91 
92  // id name
93  void const *tag_data[] = {fe_name.c_str()};
94  int tag_sizes[1];
95  tag_sizes[0] = fe_name.size();
96  CHKERR get_moab().tag_set_by_ptr(th_FEName, &meshset, 1, tag_data, tag_sizes);
97 
98  // add FiniteElement
99  std::pair<FiniteElement_multiIndex::iterator, bool> p = finiteElements.insert(
100  boost::shared_ptr<FiniteElement>(new FiniteElement(moab, meshset)));
101  if (!p.second)
102  SETERRQ(cOmm, MOFEM_OPERATION_UNSUCCESSFUL, "FiniteElement not inserted");
103  if (verb > 0) {
104  std::ostringstream ss;
105  ss << "add finite element: " << fe_name << std::endl;
106  PetscPrintf(cOmm, ss.str().c_str());
107  }
109 }
110 
113  const EntityType type,
114  ElementAdjacencyFunct function) {
116  *buildMoFEM &= 1 << 0;
118  FiniteElements_by_name;
119  FiniteElements_by_name &finite_element_name_set =
121  FiniteElements_by_name::iterator it_fe =
122  finite_element_name_set.find(fe_name);
123  if (it_fe == finite_element_name_set.end())
124  SETERRQ(cOmm, MOFEM_NOT_FOUND,
125  "This finite element is not defined (advise: check spelling)");
126  boost::shared_ptr<FiniteElement> fe;
127  fe = *it_fe;
128  fe->elementAdjacencyTable[type] = function;
130 }
131 
134  const std::string &name_data) {
136  *buildMoFEM &= 1 << 0;
138  FiniteElements_by_name;
139  FiniteElements_by_name &finite_element_name_set =
141  FiniteElements_by_name::iterator it_fe =
142  finite_element_name_set.find(fe_name);
143  if (it_fe == finite_element_name_set.end())
144  SETERRQ(cOmm, MOFEM_NOT_FOUND,
145  "This finite element is not defined (advise: check spelling)");
146  bool success = finite_element_name_set.modify(
147  it_fe, FiniteElement_change_bit_add(getBitFieldId(name_data)));
148  if (!success)
149  SETERRQ(cOmm, MOFEM_OPERATION_UNSUCCESSFUL, "modification unsuccessful");
151 }
152 
154 Core::modify_finite_element_add_field_row(const std::string &fe_name,
155  const std::string &name_row) {
157  *buildMoFEM &= 1 << 0;
159  FiniteElements_by_name;
160  FiniteElements_by_name &finite_element_name_set =
162  FiniteElements_by_name::iterator it_fe =
163  finite_element_name_set.find(fe_name);
164  if (it_fe == finite_element_name_set.end())
165  SETERRQ1(cOmm, MOFEM_NOT_FOUND, "this < %s > is not there",
166  fe_name.c_str());
167  bool success = finite_element_name_set.modify(
169  if (!success)
170  SETERRQ(cOmm, MOFEM_OPERATION_UNSUCCESSFUL, "modification unsuccessful");
172 }
173 
175 Core::modify_finite_element_add_field_col(const std::string &fe_name,
176  const std::string &name_col) {
178  *buildMoFEM &= 1 << 0;
180  FiniteElements_by_name;
181  FiniteElements_by_name &finite_element_name_set =
183  FiniteElements_by_name::iterator it_fe =
184  finite_element_name_set.find(fe_name);
185  if (it_fe == finite_element_name_set.end())
186  SETERRQ(cOmm, MOFEM_OPERATION_UNSUCCESSFUL, "this FiniteElement is there");
187  bool success = finite_element_name_set.modify(
189  if (!success)
190  SETERRQ(cOmm, MOFEM_OPERATION_UNSUCCESSFUL, "modification unsuccessful");
192 }
193 
196  const std::string &name_data) {
198  *buildMoFEM &= 1 << 0;
199  auto &finite_element_name_set =
201  auto it_fe = finite_element_name_set.find(fe_name);
202  if (it_fe == finite_element_name_set.end())
203  SETERRQ(cOmm, MOFEM_NOT_FOUND, "this FiniteElement is there");
204  bool success = finite_element_name_set.modify(
205  it_fe, FiniteElement_change_bit_off(getBitFieldId(name_data)));
206  if (!success)
207  SETERRQ(cOmm, MOFEM_OPERATION_UNSUCCESSFUL, "modification unsuccessful");
209 }
210 
212 Core::modify_finite_element_off_field_row(const std::string &fe_name,
213  const std::string &name_row) {
215  *buildMoFEM &= 1 << 0;
216  auto &finite_element_name_set =
218  auto it_fe = finite_element_name_set.find(fe_name);
219  if (it_fe == finite_element_name_set.end())
220  SETERRQ1(cOmm, MOFEM_NOT_FOUND, "this < %s > is not there",
221  fe_name.c_str());
222  bool success = finite_element_name_set.modify(
224  if (!success)
225  SETERRQ(cOmm, MOFEM_OPERATION_UNSUCCESSFUL, "modification unsuccessful");
227 }
228 
230 Core::modify_finite_element_off_field_col(const std::string &fe_name,
231  const std::string &name_col) {
233  *buildMoFEM &= 1 << 0;
234  auto &finite_element_name_set =
236  auto it_fe = finite_element_name_set.find(fe_name);
237  if (it_fe == finite_element_name_set.end())
238  SETERRQ(cOmm, MOFEM_NOT_FOUND, "this FiniteElement is there");
239  bool success = finite_element_name_set.modify(
241  if (!success)
242  SETERRQ(cOmm, MOFEM_OPERATION_UNSUCCESSFUL, "modification unsuccessful");
244 }
245 
246 BitFEId Core::getBitFEId(const std::string &name) const {
248  FiniteElements_by_name;
249  const FiniteElements_by_name &set =
251  FiniteElements_by_name::iterator miit = set.find(name);
252  if (miit == set.end())
254  ("finite element < " + name + " > not found (top tip: check spelling)")
255  .c_str());
256  return (*miit)->getId();
257 }
258 
259 std::string Core::getBitFEIdName(const BitFEId id) const {
261  finiteElements_by_id;
262  const finiteElements_by_id &set = finiteElements.get<BitFEId_mi_tag>();
263  finiteElements_by_id::iterator miit = set.find(id);
264  assert(miit != set.end());
265  return (*miit)->getName();
266 }
267 
270  finiteElements_by_id;
271  const finiteElements_by_id &set = finiteElements.get<BitFEId_mi_tag>();
272  finiteElements_by_id::iterator miit = set.find(id);
273  if (miit == set.end())
274  THROW_MESSAGE("finite element not found");
275  return (*miit)->meshset;
276 }
277 
278 EntityHandle Core::get_finite_element_meshset(const std::string &name) const {
280 }
281 
283 Core::get_finite_element_entities_by_dimension(const std::string name, int dim,
284  Range &ents) const {
285 
287 
289  CHKERR get_moab().get_entities_by_dimension(meshset, dim, ents, true);
291 }
292 
294  EntityType type,
295  Range &ents) const {
296 
298 
300  CHKERR get_moab().get_entities_by_type(meshset, type, ents, true);
301 
303 }
304 
307  Range &ents) const {
308 
310 
312  CHKERR get_moab().get_entities_by_handle(meshset, ents, true);
313 
315 }
316 
320  finiteElements_by_id;
321  const finiteElements_by_id &BitFEId_set =
323  finiteElements_by_id::iterator miit = BitFEId_set.begin();
324  for (; miit != BitFEId_set.end(); miit++) {
325  std::ostringstream ss;
326  ss << *miit << std::endl;
327  PetscSynchronizedPrintf(cOmm, ss.str().c_str());
328  }
329  PetscSynchronizedFlush(cOmm, PETSC_STDOUT);
331 }
332 
334  const EntityHandle meshset, const EntityType type, const std::string &name,
335  const bool recursive) {
336  *buildMoFEM &= 1 << 0;
337  EntityHandle idm = no_handle;
339 
341  Range ents;
342  CHKERR get_moab().get_entities_by_type(meshset, type, ents, recursive);
343  CHKERR getInterface<BitRefManager>()->setElementsBitRefLevel(ents);
344  CHKERR get_moab().add_entities(idm, ents);
345 
347 }
348 
351  const int dim, const std::string &name,
352  const bool recursive) {
353  EntityHandle idm = no_handle;
354  *buildMoFEM &= 1 << 0;
357  Range ents;
358  CHKERR get_moab().get_entities_by_dimension(meshset, dim, ents, recursive);
359  CHKERR getInterface<BitRefManager>()->setElementsBitRefLevel(ents);
360  CHKERR get_moab().add_entities(idm, ents);
362 }
363 
365  const Range &ents, const EntityType type, const std::string &name) {
366  EntityHandle idm = no_handle;
367  *buildMoFEM &= 1 << 0;
370  CHKERR getInterface<BitRefManager>()->setElementsBitRefLevel(
371  ents.subset_by_type(type));
372  CHKERR get_moab().add_entities(idm, ents.subset_by_type(type));
374 } // namespace MoFEM
375 
377 Core::add_ents_to_finite_element_by_dim(const Range &ents, const int dim,
378  const std::string &name) {
379  EntityHandle idm = no_handle;
380  *buildMoFEM &= 1 << 0;
383  CHKERR getInterface<BitRefManager>()->setElementsBitRefLevel(
384  ents.subset_by_dimension(dim));
385  CHKERR get_moab().add_entities(idm, ents.subset_by_dimension(dim));
387 }
388 
390  const EntityHandle meshset, const std::string &name, const bool recursive) {
392  CHKERR add_ents_to_finite_element_by_type(meshset, MBEDGE, name, recursive);
394 }
397  const std::string &name) {
399  CHKERR add_ents_to_finite_element_by_type(edges, MBEDGE, name);
401 }
404  const std::string &name) {
406  CHKERR add_ents_to_finite_element_by_type(vert, MBVERTEX, name);
408 }
411  const std::string &name) {
413  CHKERR add_ents_to_finite_element_by_type(tris, MBTRI, name);
415 }
417  const EntityHandle meshset, const std::string &name, const bool recursive) {
419  CHKERR add_ents_to_finite_element_by_type(meshset, MBTRI, name, recursive);
421 }
424  const std::string &name) {
426  CHKERR add_ents_to_finite_element_by_type(tets, MBTET, name);
428 }
430  const EntityHandle meshset, const std::string &name, const bool recursive) {
432  CHKERR add_ents_to_finite_element_by_type(meshset, MBTET, name, recursive);
434 }
437  const std::string &name) {
439  CHKERR add_ents_to_finite_element_by_type(prims, MBPRISM, name);
441 }
443  const EntityHandle meshset, const std::string &name, const bool recursive) {
445  CHKERR add_ents_to_finite_element_by_type(meshset, MBPRISM, name, recursive);
447 }
448 
451  const std::string &name,
452  EntityType type, int verb) {
455  type, verb);
456 
458 }
459 
461  const BitRefLevel &bit, const BitRefLevel &mask, const std::string &name,
462  EntityType type, int verb) {
464  CHKERR add_ents_to_finite_element_by_bit_ref(bit, mask, name, type, verb);
465 
467 }
468 
470  const BitRefLevel &bit, const BitRefLevel &mask, const std::string &name,
471  EntityType type, int verb) {
473 
474  if (verb == -1)
475  verb = verbose;
476  *buildMoFEM &= 1 << 0;
477  const BitFEId id = getBitFEId(name);
479  typedef RefElement_multiIndex::index<EntType_mi_tag>::type refMoabFE_by_type;
480  refMoabFE_by_type &ref_MoFEMFiniteElement =
482  refMoabFE_by_type::iterator miit = ref_MoFEMFiniteElement.lower_bound(type);
483  refMoabFE_by_type::iterator hi_miit =
484  ref_MoFEMFiniteElement.upper_bound(type);
485  if (verb > 1) {
486  PetscSynchronizedPrintf(cOmm, "nb. ref elements in database %d\n",
487  std::distance(miit, hi_miit));
488  }
489  int nb_add_FEs = 0;
490  for (; miit != hi_miit; miit++) {
491  BitRefLevel bit2 = miit->get()->getBitRefLevel();
492  if ((bit2 & mask) != bit2)
493  continue;
494  if ((bit2 & bit).any()) {
495  EntityHandle ent = miit->get()->getRefEnt();
496  CHKERR get_moab().add_entities(idm, &ent, 1);
497  nb_add_FEs++;
498  }
499  }
500  if (verb > 0) {
501  std::ostringstream ss;
502  ss << "Add Nb. FEs " << nb_add_FEs << " form BitRef " << bit << std::endl;
503  PetscSynchronizedPrintf(cOmm, "%s", ss.str().c_str());
504  PetscSynchronizedFlush(cOmm, PETSC_STDOUT);
505  }
506 
508 }
509 
511  const EntityHandle meshset, const std::string &name, const bool recursive) {
513  *buildMoFEM &= 1 << 0;
514  const BitFEId id = getBitFEId(name);
516  if (recursive == false) {
517  CHKERR get_moab().add_entities(idm, &meshset, 1);
518  } else {
519  Range meshsets;
520  CHKERR get_moab().get_entities_by_type(meshset, MBENTITYSET, meshsets,
521  false);
522  CHKERR get_moab().add_entities(idm, meshsets);
523  }
525 }
526 
527 template <int I> struct BuildFiniteElements {
528 
529  template <typename T1, typename T2>
530  static inline void addToData(T1 &range_dit, T2 &fe_vec) {
531  static_assert(I == DATA, "t should be set to DATA");
532 
533  for (auto dit = range_dit.first; dit != range_dit.second; ++dit) {
534  const EntityHandle dof_ent = dit->get()->getEnt();
535  // Fill array
536  for (auto fe_it : fe_vec) {
537  if (auto fe_ptr = fe_it.lock()) {
538  // Add FEDofEntity, first create dofs, one by one, note that memory
539  // is already reserved. Then create shared pointers and finally add
540  // th_FEName to element multi-index
541  // There are data dofs on this element
542  auto &side_number_ptr = fe_ptr->getSideNumberPtr(dof_ent);
543  fe_ptr->getDofsSequence().lock()->emplace_back(side_number_ptr, *dit);
544  }
545  }
546  }
547  }
548 
549  template <typename T> static inline void emplaceHint(T &fe_vec) {
550  static_assert(I == DATA, "t should be set to DATA");
551 
552  // Add to data in FE
553  for (auto fe_it : fe_vec) {
554  if (auto fe_ptr = fe_it.lock()) {
555  // It is a but unsafe, since if one mess up something in
556  // buildFiniteElements, weak_ptr will not give pointer
557  auto data_dofs_array_vec = fe_ptr->getDofsSequence().lock();
558  // Create shared pointers vector
559  auto hint = fe_ptr->data_dofs->end();
560  for (auto &vit : *data_dofs_array_vec)
561  hint =
562  fe_ptr->data_dofs->emplace_hint(hint, data_dofs_array_vec, &vit);
563  }
564  }
565  }
566 };
567 
569 Core::buildFiniteElements(const boost::shared_ptr<FiniteElement> &fe,
570  const Range *ents_ptr, int verb) {
572  if (verb == DEFAULT_VERBOSITY)
573  verb = verbose;
574 
575  auto &fields_by_id = fIelds.get<BitFieldId_mi_tag>();
576 
577  // Get id of mofem fields for row, col and data
578  enum IntLoop { ROW = 0, COL, DATA, LAST };
579  std::array<BitFieldId, LAST> fe_fields = {fe.get()->getBitFieldIdRow(),
580  fe.get()->getBitFieldIdCol(),
581  fe.get()->getBitFieldIdData()};
582 
583  // Get finite element meshset
584  EntityHandle meshset = get_finite_element_meshset(fe.get()->getId());
585 
586  // Get entities from finite element meshset // if meshset
587  Range fe_ents;
588  CHKERR get_moab().get_entities_by_handle(meshset, fe_ents, false);
589 
590  if (ents_ptr)
591  fe_ents = intersect(fe_ents, *ents_ptr);
592 
593  // Map entity uid to pointers
594  typedef std::vector<boost::weak_ptr<EntFiniteElement>> VecOfWeakFEPtrs;
595  typedef std::map<const UId *, VecOfWeakFEPtrs> MapEntUIdAndVecOfWeakFEPtrs;
596  MapEntUIdAndVecOfWeakFEPtrs ent_uid_and_fe_vec;
597  std::map<EntityHandle, boost::shared_ptr<std::vector<FEDofEntity>>>
598  data_dofs_array;
599  VecOfWeakFEPtrs processed_fes;
600  processed_fes.reserve(fe_ents.size());
601 
602  int last_row_field_ents_view_size = 0;
603  int last_col_field_ents_view_size = 0;
604 
605  // View of field entities on element
606  FieldEntity_vector_view data_field_ents_view;
607 
608  // Loop meshset finite element ents and add finite elements
609  for (Range::const_pair_iterator peit = fe_ents.const_pair_begin();
610  peit != fe_ents.const_pair_end(); peit++) {
611 
612  EntityHandle first = peit->first;
613  EntityHandle second = peit->second;
614 
615  // Find range of ref entities that is sequence
616  // note: iterator is a wrapper
617  // check if is in refinedFiniteElements database
618  auto ref_fe_miit =
619  refinedFiniteElements.get<Ent_mi_tag>().lower_bound(first);
620  if (ref_fe_miit == refinedFiniteElements.get<Ent_mi_tag>().end()) {
621  std::ostringstream ss;
622  ss << "refinedFiniteElements not in database ent = " << first;
623  ss << " type " << get_moab().type_from_handle(first);
624  ss << " " << *fe;
625  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, ss.str().c_str());
626  }
627  auto hi_ref_fe_miit =
628  refinedFiniteElements.get<Ent_mi_tag>().upper_bound(second);
629 
630  EntFiniteElement_multiIndex::iterator hint_p = entsFiniteElements.end();
631  for (; ref_fe_miit != hi_ref_fe_miit; ref_fe_miit++) {
632 
633  // Add finite element to database
634  hint_p = entsFiniteElements.emplace_hint(
635  hint_p, boost::make_shared<EntFiniteElement>(*ref_fe_miit, fe));
636  processed_fes.emplace_back(*hint_p);
637  auto fe_raw_ptr = hint_p->get();
638 
639  // Allocate space for etities view
640  data_field_ents_view.clear();
641  fe_raw_ptr->row_field_ents_view->reserve(last_row_field_ents_view_size);
642  // Create shared pointer for entities view
643  if (fe_fields[ROW] == fe_fields[COL]) {
644  fe_raw_ptr->col_field_ents_view = fe_raw_ptr->row_field_ents_view;
645  } else {
646  // row and columns are diffent
647  if (fe_raw_ptr->col_field_ents_view == fe_raw_ptr->row_field_ents_view)
648  fe_raw_ptr->col_field_ents_view =
649  boost::make_shared<FieldEntity_vector_view>();
650  fe_raw_ptr->col_field_ents_view->reserve(last_col_field_ents_view_size);
651  }
652 
653  int nb_dofs_on_data = 0;
654 
655  // Iterate over all field and check which one is on the element
656  for (unsigned int ii = 0; ii != BitFieldId().size(); ++ii) {
657 
658  // Common field id for ROW, COL and DATA
659  BitFieldId id_common = 0;
660  // Check if the field (ii) is added to finite element
661  for (int ss = 0; ss < LAST; ss++) {
662  id_common |= fe_fields[ss] & BitFieldId().set(ii);
663  }
664  if (id_common.none())
665  continue;
666 
667  // Find in database data associated with the field (ii)
668  const BitFieldId field_id = BitFieldId().set(ii);
669  auto miit = fields_by_id.find(field_id);
670  if (miit == fields_by_id.end()) {
671  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
672  "Data inconsistency");
673  }
674 
675  // Loop over adjacencies of element and find field entities on those
676  // adjacencies, that create hash map map_uid_fe which is used later
677  const std::string field_name = miit->get()->getName();
678  const bool add_to_data = (field_id & fe_fields[DATA]).any();
679  const bool add_to_row = (field_id & fe_fields[ROW]).any();
680  const bool add_to_col = (field_id & fe_fields[COL]).any();
681 
682  // Entities adjacent to entities
683  Range adj_ents;
684 
685  // Resolve entities on element, those entities are used to build tag
686  // with dof uids on finite element tag
687  CHKERR fe_raw_ptr->getElementAdjacency(*miit, adj_ents);
688 
689  for (Range::pair_iterator p_eit = adj_ents.pair_begin();
690  p_eit != adj_ents.pair_end(); ++p_eit) {
691 
692  const EntityHandle first = p_eit->first;
693  const EntityHandle second = p_eit->second;
694 
695  typedef FieldEntity_multiIndex::index<
696  Composite_Name_And_Ent_mi_tag>::type FieldEntityByComposite;
697  auto &field_ents_by_name_and_ent =
699  FieldEntityByComposite::iterator meit;
700 
701  // If one entity in the pair search for one, otherwise search for
702  // range
703  if (first == second)
704  meit = field_ents_by_name_and_ent.find(
705  boost::make_tuple(field_name, first));
706  else
707  meit = field_ents_by_name_and_ent.lower_bound(
708  boost::make_tuple(field_name, first));
709 
710  if (meit != field_ents_by_name_and_ent.end()) {
711 
712  decltype(meit) hi_meit;
713 
714  if (first == second) {
715  hi_meit = meit;
716  ++hi_meit;
717  } else
718  hi_meit = field_ents_by_name_and_ent.upper_bound(
719  boost::make_tuple(field_name, second));
720 
721  // Add to view and create list of finite elements with this dof UId
722  for (; meit != hi_meit; ++meit) {
723  // Add entity to map with key entity uids pointers and data
724  // finite elements weak ptrs. I using pointers to uids instead
725  // uids because this is faster.
726  const UId *uid_ptr = &(meit->get()->getGlobalUniqueId());
727  auto &fe_vec = ent_uid_and_fe_vec[uid_ptr];
728  // get number of dofs on entities to pre-allocate memory for
729  // element
730  const int nb_dofs_on_ent = meit->get()->getNbDofsOnEnt();
731  if (add_to_data) {
732  nb_dofs_on_data += nb_dofs_on_ent;
733  data_field_ents_view.emplace_back(*meit);
734  }
735  if (add_to_row) {
736  fe_raw_ptr->row_field_ents_view->emplace_back(*meit);
737  }
738  if (add_to_col) {
739  if (fe_raw_ptr->col_field_ents_view !=
740  fe_raw_ptr->row_field_ents_view)
741  fe_raw_ptr->col_field_ents_view->emplace_back(*meit);
742  }
743  // add finite element to processed list
744  fe_vec.emplace_back(*hint_p);
745  }
746  }
747  }
748  }
749 
750  // Sort field ents by uid
751  auto uid_comp = [](const auto &a, const auto &b) {
752  return a.lock()->getGlobalUniqueId() < b.lock()->getGlobalUniqueId();
753  };
754 
755  // Sort all views
756 
757  // Data
758  sort(data_field_ents_view.begin(), data_field_ents_view.end(), uid_comp);
759  for (auto e : data_field_ents_view)
760  fe_raw_ptr->data_field_ents_view->emplace_back(e);
761 
762  // Row
763  sort(fe_raw_ptr->row_field_ents_view->begin(),
764  fe_raw_ptr->row_field_ents_view->end(), uid_comp);
765  last_row_field_ents_view_size = fe_raw_ptr->row_field_ents_view->size();
766 
767  // Column
768  if (fe_raw_ptr->col_field_ents_view != fe_raw_ptr->row_field_ents_view) {
769  sort(fe_raw_ptr->col_field_ents_view->begin(),
770  fe_raw_ptr->col_field_ents_view->end(), uid_comp);
771  last_col_field_ents_view_size = fe_raw_ptr->col_field_ents_view->size();
772  }
773 
774  // Clear finite element data structures
775  fe_raw_ptr->data_dofs->clear();
776 
777  // Reserve memory for data FE Dofs
778  auto data_dofs_array_vec = boost::make_shared<std::vector<FEDofEntity>>();
779  data_dofs_array[fe_raw_ptr->getEnt()] = data_dofs_array_vec;
780  data_dofs_array_vec->reserve(nb_dofs_on_data);
781 
782  fe_raw_ptr->getDofsSequence() = data_dofs_array_vec;
783  }
784  }
785 
786  auto &dofs_by_ent_uid = dofsField.get<Unique_Ent_mi_tag>();
787 
788  // Loop over hash map, which has all entities on given elemnts
789  boost::shared_ptr<SideNumber> side_number_ptr;
790  for (auto &mit : ent_uid_and_fe_vec) {
791  auto range_dit = dofs_by_ent_uid.equal_range(*mit.first);
792  if (range_dit.first != range_dit.second) {
793  const BitFieldId field_id = range_dit.first->get()->getId();
794  if ((field_id & fe_fields[DATA]).any())
795  BuildFiniteElements<DATA>::addToData(range_dit, mit.second);
796  }
797  }
798 
800 
802 }
803 
806  if (verb == DEFAULT_VERBOSITY)
807  verb = verbose;
808 
809  // loop Finite Elements
810  for (auto &fe : finiteElements) {
811  if (verb >= VERBOSE)
812  PetscPrintf(cOmm, "Build Finite Elements %s\n", fe->getName().c_str());
813  CHKERR buildFiniteElements(fe, NULL, verb);
814  }
815 
816  if (verb >= VERBOSE) {
817  PetscSynchronizedPrintf(cOmm, "Nb. FEs %u\n", entsFiniteElements.size());
818  PetscSynchronizedFlush(cOmm, PETSC_STDOUT);
819  auto &finite_elements_by_id = entsFiniteElements.get<BitFEId_mi_tag>();
820  for (auto &fe : finiteElements) {
821  auto miit = finite_elements_by_id.lower_bound(fe->getId());
822  auto hi_miit = finite_elements_by_id.upper_bound(fe->getId());
823  int count = std::distance(miit, hi_miit);
824  std::ostringstream ss;
825  ss << *fe << " Nb. FEs " << count << std::endl;
826  PetscSynchronizedPrintf(cOmm, ss.str().c_str());
827  PetscSynchronizedFlush(cOmm, PETSC_STDOUT);
828  }
829  }
830 
831  *buildMoFEM |= 1 << 1;
833 }
834 
837  SETERRQ(cOmm, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
839 }
840 
842  const Range *const ents_ptr,
843  int verb) {
845  if (verb == -1)
846  verb = verbose;
847 
848  auto fe_miit = finiteElements.get<FiniteElement_name_mi_tag>().find(fe_name);
849  if (fe_miit == finiteElements.get<FiniteElement_name_mi_tag>().end())
850  SETERRQ1(cOmm, MOFEM_NOT_FOUND, "Finite element <%s> not found",
851  fe_name.c_str());
852 
853  if (verb >= VERBOSE)
854  PetscPrintf(cOmm, "Build Finite Elements %s\n", fe_name.c_str());
855  CHKERR buildFiniteElements(*fe_miit, ents_ptr, verb);
856 
857  if (verb >= VERBOSE) {
858  auto &finite_elements_by_id = entsFiniteElements.get<BitFEId_mi_tag>();
859  auto miit = finite_elements_by_id.lower_bound((*fe_miit)->getId());
860  auto hi_miit = finite_elements_by_id.upper_bound((*fe_miit)->getId());
861  int count = std::distance(miit, hi_miit);
862  std::ostringstream ss;
863  ss << *(*fe_miit) << " Nb. FEs " << count << std::endl;
864  PetscSynchronizedPrintf(cOmm, ss.str().c_str());
865  PetscSynchronizedFlush(cOmm, PETSC_STDOUT);
866  }
867 
868  *buildMoFEM |= 1 << 1;
870 }
871 
872 MoFEMErrorCode Core::build_adjacencies(const Range &ents, int verb) {
874  if (verb == DEFAULT_VERBOSITY)
875  verb = verbose;
876  if (!((*buildMoFEM) & BUILD_FIELD))
877  SETERRQ(cOmm, MOFEM_NOT_FOUND, "field not build");
878  if (!((*buildMoFEM) & BUILD_FE))
879  SETERRQ(cOmm, MOFEM_NOT_FOUND, "fe not build");
880  for (Range::const_pair_iterator peit = ents.pair_begin();
881  peit != ents.pair_end(); ++peit) {
882  EntFiniteElement_multiIndex::index<Ent_mi_tag>::type::iterator fit, hi_fit;
883  fit = entsFiniteElements.get<Ent_mi_tag>().lower_bound(peit->first);
884  hi_fit = entsFiniteElements.get<Ent_mi_tag>().upper_bound(peit->second);
885  for (; fit != hi_fit; ++fit) {
886  if ((*fit)->getBitFieldIdRow().none() &&
887  (*fit)->getBitFieldIdCol().none() &&
888  (*fit)->getBitFieldIdData().none())
889  continue;
890  int by = BYROW;
891  if ((*fit)->getBitFieldIdRow() != (*fit)->getBitFieldIdCol())
892  by |= BYCOL;
893  if ((*fit)->getBitFieldIdRow() != (*fit)->getBitFieldIdData())
894  by |= BYDATA;
896  auto hint = entFEAdjacencies.end();
897  for (auto e : *(*fit)->row_field_ents_view) {
898  hint = entFEAdjacencies.emplace_hint(hint, e.lock(), *fit);
899  bool success = entFEAdjacencies.modify(hint, modify_row);
900  if (!success)
901  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
902  "modification unsuccessful");
903  }
904  if ((*fit)->getBitFieldIdRow() != (*fit)->getBitFieldIdCol()) {
905  int by = BYCOL;
906  if ((*fit)->getBitFieldIdCol() != (*fit)->getBitFieldIdData())
907  by |= BYDATA;
909  auto hint = entFEAdjacencies.end();
910  for (auto e : *(*fit)->col_field_ents_view) {
911  hint = entFEAdjacencies.emplace_hint(hint, e.lock(), *fit);
912  bool success = entFEAdjacencies.modify(hint, modify_col);
913  if (!success)
914  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
915  "modification unsuccessful");
916  }
917  }
918  if ((*fit)->getBitFieldIdRow() != (*fit)->getBitFieldIdData() ||
919  (*fit)->getBitFieldIdCol() != (*fit)->getBitFieldIdData()) {
921  BYDATA);
922  auto hint = entFEAdjacencies.end();
923  for (auto &e : *(*fit)->data_field_ents_view) {
924  hint = entFEAdjacencies.emplace_hint(hint, e, *fit);
925  bool success = entFEAdjacencies.modify(hint, modify_data);
926  if (!success)
927  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
928  "modification unsuccessful");
929  }
930  }
931  }
932  }
933  if (verb >= VERY_NOISY) {
935  }
936  if (verb >= VERBOSE) {
937  PetscSynchronizedPrintf(cOmm, "Nb. entFEAdjacencies %u\n",
938  entFEAdjacencies.size());
939  PetscSynchronizedFlush(cOmm, PETSC_STDOUT);
940  }
941  *buildMoFEM |= 1 << 2;
943 }
944 
946  const BitRefLevel &mask, int verb) {
948  if (verb == -1)
949  verb = verbose;
950  Range ents;
951  CHKERR BitRefManager(*this).getEntitiesByRefLevel(bit, mask, ents);
952 
953  CHKERR build_adjacencies(ents, verb);
954 
956 }
959  if (verb == -1)
960  verb = verbose;
961  CHKERR build_adjacencies(bit, BitRefLevel().set(), verb);
962 
964 }
965 
966 EntFiniteElementByName::iterator
967 Core::get_fe_by_name_begin(const std::string &fe_name) const {
968  return entsFiniteElements.get<FiniteElement_name_mi_tag>().lower_bound(
969  fe_name);
970 }
971 EntFiniteElementByName::iterator
972 Core::get_fe_by_name_end(const std::string &fe_name) const {
973  return entsFiniteElements.get<FiniteElement_name_mi_tag>().upper_bound(
974  fe_name);
975 }
976 
978  const std::string &name) const {
980  FiniteElement_multiIndex::index<FiniteElement_name_mi_tag>::type::iterator it;
981  it = finiteElements.get<FiniteElement_name_mi_tag>().find(name);
982  if (it == finiteElements.get<FiniteElement_name_mi_tag>().end()) {
983  SETERRQ1(cOmm, 1, "finite element not found < %s >", name.c_str());
984  }
985  EntityHandle meshset = (*it)->getMeshset();
986 
987  int num_entities;
988  CHKERR get_moab().get_number_entities_by_handle(meshset, num_entities);
989 
991  (*it)->getName().c_str()) != (unsigned int)num_entities) {
992  SETERRQ1(cOmm, 1,
993  "not equal number of entities in meshset and finite elements "
994  "multiindex < %s >",
995  (*it)->getName().c_str());
996  }
998 }
999 
1002  FiniteElement_multiIndex::index<FiniteElement_name_mi_tag>::type::iterator it;
1003  it = finiteElements.get<FiniteElement_name_mi_tag>().begin();
1004  for (; it != finiteElements.get<FiniteElement_name_mi_tag>().end(); it++) {
1005  EntityHandle meshset = (*it)->getMeshset();
1006 
1007  int num_entities;
1008  CHKERR get_moab().get_number_entities_by_handle(meshset, num_entities);
1009 
1011  (*it)->getName().c_str()) != (unsigned int)num_entities) {
1012  SETERRQ1(cOmm, 1,
1013  "not equal number of entities in meshset and finite elements "
1014  "multiindex < %s >",
1015  (*it)->getName().c_str());
1016  }
1017  }
1019 }
1020 } // 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:133
EntityHandle get_finite_element_meshset(const BitFEId id) const
Definition: FECore.cpp:268
MoFEMErrorCode check_number_of_ents_in_ents_finite_element() const
check data consistency in entsFiniteElements
Definition: FECore.cpp:1000
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:510
MPI_Comm cOmm
MoFEM communicator.
Definition: Core.hpp:846
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:350
FieldEntity_multiIndex entsFields
entities on fields
Definition: Core.hpp:248
multi_index_container< boost::shared_ptr< FiniteElement >, indexed_by< hashed_unique< tag< FiniteElement_Meshset_mi_tag >, member< FiniteElement, EntityHandle, &FiniteElement::meshset > >, hashed_unique< tag< BitFEId_mi_tag >, const_mem_fun< FiniteElement, BitFEId, &FiniteElement::getId >, HashBit< BitFEId >, EqBit< BitFEId > >, ordered_unique< tag< FiniteElement_name_mi_tag >, const_mem_fun< FiniteElement, boost::string_ref, &FiniteElement::getNameRef > > > > FiniteElement_multiIndex
MultiIndex for entities for FiniteElement.
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:403
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:396
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:506
std::vector< boost::weak_ptr< FieldEntity > > FieldEntity_vector_view
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
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:283
MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)
build adjacencies
Definition: FECore.cpp:872
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:35
MoFEMErrorCode list_finite_elements() const
list finite elements in database
Definition: FECore.cpp:317
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:804
#define THROW_MESSAGE(a)
Throw MoFEM exception.
Definition: definitions.h:625
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:513
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:469
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:195
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:293
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:410
BitFEId getBitFEId(const std::string &name) const
Get field Id.
Definition: FECore.cpp:246
std::bitset< BITFEID_SIZE > BitFEId
Finite element Id.
Definition: Types.hpp:54
std::string getBitFEIdName(const BitFEId id) const
Get field name.
Definition: FECore.cpp:259
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:972
int rAnk
MOFEM communicator rank.
Definition: Core.hpp:850
static void emplaceHint(T &fe_vec)
Definition: FECore.cpp:549
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:230
MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)
add finite element
Definition: FECore.cpp:45
BitFEId getFEShift()
Return unique finite element Id.
Definition: Core.cpp:238
bool check_finite_element(const std::string &name) const
Check if finite element is in database.
Definition: FECore.cpp:35
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:212
int * buildMoFEM
keeps flags/semaphores for different stages
Definition: Core.hpp:896
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
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
std::bitset< BITFIELDID_SIZE > BitFieldId
Field Id.
Definition: Types.hpp:53
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:175
#define CHKERR
Inline error check.
Definition: definitions.h:601
multi_index_container< boost::shared_ptr< EntFiniteElement >, indexed_by< ordered_unique< tag< Unique_mi_tag >, member< EntFiniteElement, UId, &EntFiniteElement::globalUId > >, ordered_non_unique< tag< Ent_mi_tag >, const_mem_fun< EntFiniteElement, EntityHandle, &EntFiniteElement::getEnt > >, ordered_non_unique< tag< FiniteElement_name_mi_tag >, const_mem_fun< EntFiniteElement::interface_type_FiniteElement, boost::string_ref, &EntFiniteElement::getNameRef > >, ordered_non_unique< tag< BitFEId_mi_tag >, const_mem_fun< EntFiniteElement::interface_type_FiniteElement, BitFEId, &EntFiniteElement::getId >, LtBit< BitFEId > >, ordered_non_unique< tag< EntType_mi_tag >, const_mem_fun< EntFiniteElement::interface_type_RefEntity, EntityType, &EntFiniteElement::getEntType > >, ordered_non_unique< tag< Composite_Name_And_Ent_mi_tag >, composite_key< EntFiniteElement, const_mem_fun< EntFiniteElement::interface_type_FiniteElement, boost::string_ref, &EntFiniteElement::getNameRef >, const_mem_fun< EntFiniteElement, EntityHandle, &EntFiniteElement::getEnt > > > > > EntFiniteElement_multiIndex
MultiIndex container for EntFiniteElement.
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:333
Field_multiIndex fIelds
fields
Definition: Core.hpp:247
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
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:21
static void addToData(T1 &range_dit, T2 &fe_vec)
Definition: FECore.cpp:530
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:154
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:450
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:112
int verbose
Verbosity level.
Definition: Core.hpp:885
boost::function< MoFEMErrorCode(Interface &moab, const Field &field, const EntFiniteElement &fe, Range &adjacency)> ElementAdjacencyFunct
user adjacency function
Tag th_FEId
Definition: Core.hpp:206
MoFEMErrorCode get_finite_elements(const FiniteElement_multiIndex **fe_ptr) const
Get finite elements multi-index.
Definition: FECore.cpp:22
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412
MoFEMErrorCode buildFiniteElements(const boost::shared_ptr< FiniteElement > &fe, const Range *ents_ptr=NULL, int verb=DEFAULT_VERBOSITY)
Definition: FECore.cpp:569
MoFEMErrorCode get_ents_finite_elements(const EntFiniteElement_multiIndex **fe_ent_ptr) const
Get entities finite elements multi-index.
Definition: FECore.cpp:28
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
MoFEMErrorCode list_adjacencies() const
list adjacencies
Definition: FieldCore.cpp:1183
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:306
MoFEMTypes
Those types control how functions respond on arguments, f.e. error handling.
Definition: definitions.h:188
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:967
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:423