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