v0.10.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  }
49  FiniteElements_by_name;
50  FiniteElements_by_name &finite_element_name_set =
51  finiteElements.get<FiniteElement_name_mi_tag>();
52  FiniteElements_by_name::iterator it_fe =
53  finite_element_name_set.find(fe_name);
54  if (bh == MF_EXCL) {
55  if (it_fe != finite_element_name_set.end()) {
56  SETERRQ1(cOmm, MOFEM_NOT_FOUND, "this < %s > is there", fe_name.c_str());
57  }
58  } else {
59  if (it_fe != finite_element_name_set.end())
61  }
62  EntityHandle meshset;
63  CHKERR get_moab().create_meshset(MESHSET_SET, meshset);
64 
65  // id
66  BitFEId id = getFEShift();
67  CHKERR get_moab().tag_set_data(th_FEId, &meshset, 1, &id);
68 
69  // Add finite element meshset to partion meshset. In case of no elements
70  // on processor part, when mesh file is read, finite element meshset is
71  // prevented from deletion by moab reader.
72  auto add_meshset_to_partition = [&](auto meshset) {
74  const void *tag_vals[] = {&rAnk};
75  ParallelComm *pcomm = ParallelComm::get_pcomm(
76  &get_moab(), get_basic_entity_data_ptr()->pcommID);
77  Tag part_tag = pcomm->part_tag();
78  Range tagged_sets;
79  CHKERR get_moab().get_entities_by_type_and_tag(0, MBENTITYSET, &part_tag,
80  tag_vals, 1, tagged_sets,
81  moab::Interface::UNION);
82  for (auto s : tagged_sets)
83  CHKERR get_moab().add_entities(s, &meshset, 1);
85  };
86  CHKERR add_meshset_to_partition(meshset);
87 
88  // id name
89  void const *tag_data[] = {fe_name.c_str()};
90  int tag_sizes[1];
91  tag_sizes[0] = fe_name.size();
92  CHKERR get_moab().tag_set_by_ptr(th_FEName, &meshset, 1, tag_data, tag_sizes);
93 
94  // add FiniteElement
95  auto p = finiteElements.insert(
96  boost::shared_ptr<FiniteElement>(new FiniteElement(moab, meshset)));
97  if (!p.second)
98  SETERRQ(cOmm, MOFEM_OPERATION_UNSUCCESSFUL, "FiniteElement not inserted");
99 
100  if (verb > QUIET)
101  MOFEM_LOG("WORLD", Sev::inform) << "Add finite element " << fe_name;
102 
104 }
105 
108  const EntityType type,
109  ElementAdjacencyFunct function) {
111  *buildMoFEM &= 1 << 0;
113  FiniteElements_by_name;
114  FiniteElements_by_name &finite_element_name_set =
115  finiteElements.get<FiniteElement_name_mi_tag>();
116  FiniteElements_by_name::iterator it_fe =
117  finite_element_name_set.find(fe_name);
118  if (it_fe == finite_element_name_set.end())
119  SETERRQ(cOmm, MOFEM_NOT_FOUND,
120  "This finite element is not defined (advise: check spelling)");
121  boost::shared_ptr<FiniteElement> fe;
122  fe = *it_fe;
123  fe->elementAdjacencyTable[type] = function;
125 }
126 
129  const std::string &name_data) {
131  *buildMoFEM &= 1 << 0;
133  FiniteElements_by_name;
134  FiniteElements_by_name &finite_element_name_set =
135  finiteElements.get<FiniteElement_name_mi_tag>();
136  FiniteElements_by_name::iterator it_fe =
137  finite_element_name_set.find(fe_name);
138  if (it_fe == finite_element_name_set.end())
139  SETERRQ(cOmm, MOFEM_NOT_FOUND,
140  "This finite element is not defined (advise: check spelling)");
141  bool success = finite_element_name_set.modify(
142  it_fe, FiniteElement_change_bit_add(getBitFieldId(name_data)));
143  if (!success)
144  SETERRQ(cOmm, MOFEM_OPERATION_UNSUCCESSFUL, "modification unsuccessful");
146 }
147 
149 Core::modify_finite_element_add_field_row(const std::string &fe_name,
150  const std::string &name_row) {
152  *buildMoFEM &= 1 << 0;
154  FiniteElements_by_name;
155  FiniteElements_by_name &finite_element_name_set =
156  finiteElements.get<FiniteElement_name_mi_tag>();
157  FiniteElements_by_name::iterator it_fe =
158  finite_element_name_set.find(fe_name);
159  if (it_fe == finite_element_name_set.end())
160  SETERRQ1(cOmm, MOFEM_NOT_FOUND, "this < %s > is not there",
161  fe_name.c_str());
162  bool success = finite_element_name_set.modify(
163  it_fe, FiniteElement_row_change_bit_add(getBitFieldId(name_row)));
164  if (!success)
165  SETERRQ(cOmm, MOFEM_OPERATION_UNSUCCESSFUL, "modification unsuccessful");
167 }
168 
170 Core::modify_finite_element_add_field_col(const std::string &fe_name,
171  const std::string &name_col) {
173  *buildMoFEM &= 1 << 0;
175  FiniteElements_by_name;
176  FiniteElements_by_name &finite_element_name_set =
177  finiteElements.get<FiniteElement_name_mi_tag>();
178  FiniteElements_by_name::iterator it_fe =
179  finite_element_name_set.find(fe_name);
180  if (it_fe == finite_element_name_set.end())
181  SETERRQ(cOmm, MOFEM_OPERATION_UNSUCCESSFUL, "this FiniteElement is there");
182  bool success = finite_element_name_set.modify(
183  it_fe, FiniteElement_col_change_bit_add(getBitFieldId(name_col)));
184  if (!success)
185  SETERRQ(cOmm, MOFEM_OPERATION_UNSUCCESSFUL, "modification unsuccessful");
187 }
188 
191  const std::string &name_data) {
193  *buildMoFEM &= 1 << 0;
194  auto &finite_element_name_set =
195  finiteElements.get<FiniteElement_name_mi_tag>();
196  auto it_fe = finite_element_name_set.find(fe_name);
197  if (it_fe == finite_element_name_set.end())
198  SETERRQ(cOmm, MOFEM_NOT_FOUND, "this FiniteElement is there");
199  bool success = finite_element_name_set.modify(
200  it_fe, FiniteElement_change_bit_off(getBitFieldId(name_data)));
201  if (!success)
202  SETERRQ(cOmm, MOFEM_OPERATION_UNSUCCESSFUL, "modification unsuccessful");
204 }
205 
207 Core::modify_finite_element_off_field_row(const std::string &fe_name,
208  const std::string &name_row) {
210  *buildMoFEM &= 1 << 0;
211  auto &finite_element_name_set =
212  finiteElements.get<FiniteElement_name_mi_tag>();
213  auto it_fe = finite_element_name_set.find(fe_name);
214  if (it_fe == finite_element_name_set.end())
215  SETERRQ1(cOmm, MOFEM_NOT_FOUND, "this < %s > is not there",
216  fe_name.c_str());
217  bool success = finite_element_name_set.modify(
218  it_fe, FiniteElement_row_change_bit_off(getBitFieldId(name_row)));
219  if (!success)
220  SETERRQ(cOmm, MOFEM_OPERATION_UNSUCCESSFUL, "modification unsuccessful");
222 }
223 
225 Core::modify_finite_element_off_field_col(const std::string &fe_name,
226  const std::string &name_col) {
228  *buildMoFEM &= 1 << 0;
229  auto &finite_element_name_set =
230  finiteElements.get<FiniteElement_name_mi_tag>();
231  auto it_fe = finite_element_name_set.find(fe_name);
232  if (it_fe == finite_element_name_set.end())
233  SETERRQ(cOmm, MOFEM_NOT_FOUND, "this FiniteElement is there");
234  bool success = finite_element_name_set.modify(
235  it_fe, FiniteElement_col_change_bit_off(getBitFieldId(name_col)));
236  if (!success)
237  SETERRQ(cOmm, MOFEM_OPERATION_UNSUCCESSFUL, "modification unsuccessful");
239 }
240 
241 BitFEId Core::getBitFEId(const std::string &name) const {
242  auto &fe_by_id = finiteElements.get<FiniteElement_name_mi_tag>();
243  auto miit = fe_by_id.find(name);
244  if (miit == fe_by_id.end())
246  ("finite element < " + name + " > not found (top tip: check spelling)")
247  .c_str());
248  return (*miit)->getId();
249 }
250 
251 std::string Core::getBitFEIdName(const BitFEId id) const {
252  auto &fe_by_id = finiteElements.get<BitFEId_mi_tag>();
253  auto miit = fe_by_id.find(id);
254  if (miit == fe_by_id.end())
255  THROW_MESSAGE("finite element not found");
256  return (*miit)->getName();
257 }
258 
260  auto &fe_by_id = finiteElements.get<BitFEId_mi_tag>();
261  auto miit = fe_by_id.find(id);
262  if (miit == fe_by_id.end())
263  THROW_MESSAGE("finite element not found");
264  return (*miit)->meshset;
265 }
266 
267 EntityHandle Core::get_finite_element_meshset(const std::string &name) const {
268  return get_finite_element_meshset(getBitFEId(name));
269 }
270 
273  Range &ents) const {
274 
276 
277  EntityHandle meshset = get_finite_element_meshset(name);
278  CHKERR get_moab().get_entities_by_dimension(meshset, dim, ents, true);
280 }
281 
283  EntityType type,
284  Range &ents) const {
285 
287 
288  EntityHandle meshset = get_finite_element_meshset(name);
289  CHKERR get_moab().get_entities_by_type(meshset, type, ents, true);
290 
292 }
293 
296  Range &ents) const {
297 
299 
300  EntityHandle meshset = get_finite_element_meshset(name);
301  CHKERR get_moab().get_entities_by_handle(meshset, ents, true);
302 
304 }
305 
308  for (auto &fe : finiteElements.get<FiniteElement_name_mi_tag>())
309  MOFEM_LOG("SYNC", Sev::inform) << fe;
310 
311  MOFEM_LOG_SYNCHRONISE(cOmm);
313 }
314 
316  const EntityHandle meshset, const EntityType type, const std::string &name,
317  const bool recursive) {
318  *buildMoFEM &= 1 << 0;
319  EntityHandle idm = no_handle;
321 
322  idm = get_finite_element_meshset(getBitFEId(name));
323  Range ents;
324  CHKERR get_moab().get_entities_by_type(meshset, type, ents, recursive);
325  CHKERR getInterface<BitRefManager>()->setElementsBitRefLevel(ents);
326  CHKERR get_moab().add_entities(idm, ents);
327 
329 }
330 
333  const int dim, const std::string &name,
334  const bool recursive) {
335  EntityHandle idm = no_handle;
336  *buildMoFEM &= 1 << 0;
338  idm = get_finite_element_meshset(getBitFEId(name));
339  Range ents;
340  CHKERR get_moab().get_entities_by_dimension(meshset, dim, ents, recursive);
341  CHKERR getInterface<BitRefManager>()->setElementsBitRefLevel(ents);
342  CHKERR get_moab().add_entities(idm, ents);
344 }
345 
347  const Range &ents, const EntityType type, const std::string &name) {
348  EntityHandle idm = no_handle;
349  *buildMoFEM &= 1 << 0;
351  idm = get_finite_element_meshset(getBitFEId(name));
352  CHKERR getInterface<BitRefManager>()->setElementsBitRefLevel(
353  ents.subset_by_type(type));
354  CHKERR get_moab().add_entities(idm, ents.subset_by_type(type));
356 } // namespace MoFEM
357 
359 Core::add_ents_to_finite_element_by_dim(const Range &ents, const int dim,
360  const std::string &name) {
361  EntityHandle idm = no_handle;
362  *buildMoFEM &= 1 << 0;
364  idm = get_finite_element_meshset(getBitFEId(name));
365  CHKERR getInterface<BitRefManager>()->setElementsBitRefLevel(
366  ents.subset_by_dimension(dim));
367  CHKERR get_moab().add_entities(idm, ents.subset_by_dimension(dim));
369 }
370 
372  const EntityHandle meshset, const std::string &name, const bool recursive) {
374  CHKERR add_ents_to_finite_element_by_type(meshset, MBEDGE, name, recursive);
376 }
379  const std::string &name) {
381  CHKERR add_ents_to_finite_element_by_type(edges, MBEDGE, name);
383 }
386  const std::string &name) {
388  CHKERR add_ents_to_finite_element_by_type(vert, MBVERTEX, name);
390 }
393  const std::string &name) {
395  CHKERR add_ents_to_finite_element_by_type(tris, MBTRI, name);
397 }
399  const EntityHandle meshset, const std::string &name, const bool recursive) {
401  CHKERR add_ents_to_finite_element_by_type(meshset, MBTRI, name, recursive);
403 }
406  const std::string &name) {
408  CHKERR add_ents_to_finite_element_by_type(tets, MBTET, name);
410 }
412  const EntityHandle meshset, const std::string &name, const bool recursive) {
414  CHKERR add_ents_to_finite_element_by_type(meshset, MBTET, name, recursive);
416 }
419  const std::string &name) {
421  CHKERR add_ents_to_finite_element_by_type(prims, MBPRISM, name);
423 }
425  const EntityHandle meshset, const std::string &name, const bool recursive) {
427  CHKERR add_ents_to_finite_element_by_type(meshset, MBPRISM, name, recursive);
429 }
430 
433  const std::string &name,
434  EntityType type, int verb) {
436  CHKERR add_ents_to_finite_element_by_bit_ref(bit, BitRefLevel().set(), name,
437  type, verb);
438 
440 }
441 
443  const BitRefLevel &bit, const BitRefLevel &mask, const std::string &name,
444  EntityType type, int verb) {
446  CHKERR add_ents_to_finite_element_by_bit_ref(bit, mask, name, type, verb);
447 
449 }
450 
452  const BitRefLevel &bit, const BitRefLevel &mask, const std::string &name,
453  EntityType type, int verb) {
455 
456  if (verb == -1)
457  verb = verbose;
458  *buildMoFEM &= 1 << 0;
459  const BitFEId id = getBitFEId(name);
460  const EntityHandle idm = get_finite_element_meshset(id);
461 
462  auto &ref_MoFEMFiniteElement = refinedFiniteElements.get<Ent_mi_tag>();
463  auto miit = ref_MoFEMFiniteElement.lower_bound(get_id_for_min_type(type));
464  auto hi_miit = ref_MoFEMFiniteElement.upper_bound(get_id_for_max_type(type));
465 
466  int nb_add_fes = 0;
467  for (; miit != hi_miit; miit++) {
468  BitRefLevel bit2 = miit->get()->getBitRefLevel();
469  if ((bit2 & mask) != bit2)
470  continue;
471  if ((bit2 & bit).any()) {
472  EntityHandle ent = miit->get()->getEnt();
473  CHKERR get_moab().add_entities(idm, &ent, 1);
474  nb_add_fes++;
475  }
476  }
477 
478  MOFEM_LOG("SYNC", Sev::inform)
479  << "Finite element " << name << " added. Nb. of elements added "
480  << nb_add_fes << " out of " << std::distance(miit, hi_miit);
481 
483 
485 }
486 
488  const EntityHandle meshset, const std::string &name, const bool recursive) {
490  *buildMoFEM &= 1 << 0;
491  const BitFEId id = getBitFEId(name);
492  const EntityHandle idm = get_finite_element_meshset(id);
493  if (recursive == false) {
494  CHKERR get_moab().add_entities(idm, &meshset, 1);
495  } else {
496  Range meshsets;
497  CHKERR get_moab().get_entities_by_type(meshset, MBENTITYSET, meshsets,
498  false);
499  CHKERR get_moab().add_entities(idm, meshsets);
500  }
502 }
503 
505 Core::buildFiniteElements(const boost::shared_ptr<FiniteElement> &fe,
506  const Range *ents_ptr, int verb) {
508  if (verb == DEFAULT_VERBOSITY)
509  verb = verbose;
510 
511  if (verb > QUIET)
512  MOFEM_LOG("SYNC", Sev::verbose)
513  << "Build Finite Elements " << fe->getName();
514 
515  auto &fields_by_id = fIelds.get<BitFieldId_mi_tag>();
516 
517  // Get id of mofem fields for row, col and data
518  enum IntLoop { ROW = 0, COL, DATA, LAST };
519  std::array<BitFieldId, LAST> fe_fields = {fe.get()->getBitFieldIdRow(),
520  fe.get()->getBitFieldIdCol(),
521  fe.get()->getBitFieldIdData()};
522 
523  // Get finite element meshset
524  EntityHandle meshset = get_finite_element_meshset(fe.get()->getId());
525 
526  // Get entities from finite element meshset // if meshset
527  Range fe_ents;
528  CHKERR get_moab().get_entities_by_handle(meshset, fe_ents, false);
529 
530  if (ents_ptr)
531  fe_ents = intersect(fe_ents, *ents_ptr);
532 
533  // Map entity uid to pointers
534  typedef std::vector<boost::weak_ptr<EntFiniteElement>> VecOfWeakFEPtrs;
535  typedef std::map<const UId *, VecOfWeakFEPtrs> MapEntUIdAndVecOfWeakFEPtrs;
536  MapEntUIdAndVecOfWeakFEPtrs ent_uid_and_fe_vec;
537  VecOfWeakFEPtrs processed_fes;
538  processed_fes.reserve(fe_ents.size());
539 
540  int last_data_field_ents_view_size = 0;
541  int last_row_field_ents_view_size = 0;
542  int last_col_field_ents_view_size = 0;
543 
544  // Loop meshset finite element ents and add finite elements
545  for (Range::const_pair_iterator peit = fe_ents.const_pair_begin();
546  peit != fe_ents.const_pair_end(); peit++) {
547 
548  const auto first = peit->first;
549  const auto second = peit->second;
550 
551  // Find range of ref entities that is sequence
552  // note: iterator is a wrapper
553  // check if is in refinedFiniteElements database
554  auto ref_fe_miit =
555  refinedFiniteElements.get<Ent_mi_tag>().lower_bound(first);
556  if (ref_fe_miit == refinedFiniteElements.get<Ent_mi_tag>().end()) {
557  std::ostringstream ss;
558  ss << "refinedFiniteElements not in database ent = " << first;
559  ss << " type " << get_moab().type_from_handle(first);
560  ss << " " << *fe;
561  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, ss.str().c_str());
562  }
563  auto hi_ref_fe_miit =
564  refinedFiniteElements.get<Ent_mi_tag>().upper_bound(second);
565 
566  EntFiniteElement_multiIndex::iterator hint_p = entsFiniteElements.end();
567  for (; ref_fe_miit != hi_ref_fe_miit; ref_fe_miit++) {
568 
569  // Add finite element to database
570  hint_p = entsFiniteElements.emplace_hint(
571  hint_p, boost::make_shared<EntFiniteElement>(*ref_fe_miit, fe));
572  processed_fes.emplace_back(*hint_p);
573  auto fe_raw_ptr = hint_p->get();
574 
575  // Allocate space for etities view
576  bool row_as_data = false, col_as_row = false;
577  if (fe_fields[DATA] == fe_fields[ROW])
578  row_as_data = true;
579  if (fe_fields[ROW] == fe_fields[COL])
580  col_as_row = true;
581 
582  fe_raw_ptr->getDataFieldEntsPtr()->reserve(
583  last_data_field_ents_view_size);
584 
585  if (row_as_data) {
586  fe_raw_ptr->getRowFieldEntsPtr() = fe_raw_ptr->getDataFieldEntsPtr();
587  } else {
588  // row and col are diffent
589  if (fe_raw_ptr->getRowFieldEntsPtr() ==
590  fe_raw_ptr->getDataFieldEntsPtr())
591  fe_raw_ptr->getRowFieldEntsPtr() =
592  boost::make_shared<FieldEntity_vector_view>();
593  fe_raw_ptr->getRowFieldEntsPtr()->reserve(
594  last_row_field_ents_view_size);
595  }
596 
597  if (row_as_data && col_as_row) {
598  fe_raw_ptr->getColFieldEntsPtr() = fe_raw_ptr->getDataFieldEntsPtr();
599  } else if (col_as_row) {
600  fe_raw_ptr->getColFieldEntsPtr() = fe_raw_ptr->getRowFieldEntsPtr();
601  } else {
602  if (
603 
604  fe_raw_ptr->getColFieldEntsPtr() ==
605  fe_raw_ptr->getRowFieldEntsPtr() ||
606  fe_raw_ptr->getColFieldEntsPtr() ==
607  fe_raw_ptr->getDataFieldEntsPtr()
608 
609  )
610  fe_raw_ptr->getColFieldEntsPtr() =
611  boost::make_shared<FieldEntity_vector_view>();
612  fe_raw_ptr->getColFieldEntsPtr()->reserve(
613  last_col_field_ents_view_size);
614  }
615 
616  // Iterate over all field and check which one is on the element
617  for (unsigned int ii = 0; ii != BitFieldId().size(); ++ii) {
618 
619  // Common field id for ROW, COL and DATA
620  BitFieldId id_common = 0;
621  // Check if the field (ii) is added to finite element
622  for (int ss = 0; ss < LAST; ss++) {
623  id_common |= fe_fields[ss] & BitFieldId().set(ii);
624  }
625  if (id_common.none())
626  continue;
627 
628  // Find in database data associated with the field (ii)
629  const BitFieldId field_id = BitFieldId().set(ii);
630  auto miit = fields_by_id.find(field_id);
631  if (miit == fields_by_id.end())
632  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Field not found");
633  auto field_bit_number = (*miit)->getBitNumber();
634 
635  // Loop over adjacencies of element and find field entities on those
636  // adjacencies, that create hash map map_uid_fe which is used later
637  const std::string field_name = miit->get()->getName();
638  const bool add_to_data = (field_id & fe_fields[DATA]).any();
639  const bool add_to_row = (field_id & fe_fields[ROW]).any();
640  const bool add_to_col = (field_id & fe_fields[COL]).any();
641 
642  // Entities adjacent to entities
643  Range adj_ents;
644 
645  // Resolve entities on element, those entities are used to build tag
646  // with dof uids on finite element tag
647  CHKERR fe_raw_ptr->getElementAdjacency(*miit, adj_ents);
648 
649  for (Range::pair_iterator p_eit = adj_ents.pair_begin();
650  p_eit != adj_ents.pair_end(); ++p_eit) {
651 
652  const EntityHandle first = p_eit->first;
653  const EntityHandle second = p_eit->second;
654 
656  FieldEntityByComposite;
657  auto &field_ents_by_name_and_ent = entsFields.get<Unique_mi_tag>();
658  FieldEntityByComposite::iterator meit;
659 
660  const auto lo_uid =
661  FieldEntity::getLocalUniqueIdCalculate(field_bit_number, first);
662 
663  // If one entity in the pair search for one, otherwise search for
664  // range
665  if (first == second)
666  meit = field_ents_by_name_and_ent.find(lo_uid);
667  else
668  meit = field_ents_by_name_and_ent.lower_bound(lo_uid);
669 
670  if (meit != field_ents_by_name_and_ent.end()) {
671 
672  decltype(meit) hi_meit;
673 
674  if (first == second) {
675  hi_meit = meit;
676  ++hi_meit;
677  } else
678  hi_meit = field_ents_by_name_and_ent.upper_bound(
680  second));
681 
682  // Add to view and create list of finite elements with this dof UId
683  for (; meit != hi_meit; ++meit) {
684  // Add entity to map with key entity uids pointers and data
685  // finite elements weak ptrs. I using pointers to uids instead
686  // uids because this is faster.
687  const UId *uid_ptr = &(meit->get()->getLocalUniqueId());
688  auto &fe_vec = ent_uid_and_fe_vec[uid_ptr];
689  if (add_to_data) {
690  fe_raw_ptr->getDataFieldEntsPtr()->emplace_back(*meit);
691  }
692  if (add_to_row && !row_as_data) {
693  fe_raw_ptr->getRowFieldEntsPtr()->emplace_back(*meit);
694  }
695  if (add_to_col && !col_as_row) {
696  fe_raw_ptr->getColFieldEntsPtr()->emplace_back(*meit);
697  }
698 
699  // add finite element to processed list
700  fe_vec.emplace_back(*hint_p);
701  }
702  }
703  }
704  }
705 
706  // Sort field ents by uid
707  auto uid_comp = [](const auto &a, const auto &b) {
708  return a.lock()->getLocalUniqueId() < b.lock()->getLocalUniqueId();
709  };
710 
711  // Sort all views
712 
713  // Data
714  sort(fe_raw_ptr->getDataFieldEntsPtr()->begin(),
715  fe_raw_ptr->getDataFieldEntsPtr()->end(), uid_comp);
716  last_data_field_ents_view_size =
717  fe_raw_ptr->getDataFieldEntsPtr()->size();
718 
719  // Row
720  if (!row_as_data) {
721  sort(fe_raw_ptr->getRowFieldEntsPtr()->begin(),
722  fe_raw_ptr->getRowFieldEntsPtr()->end(), uid_comp);
723  last_row_field_ents_view_size =
724  fe_raw_ptr->getRowFieldEntsPtr()->size();
725  }
726 
727  // Column
728  if (!col_as_row) {
729  sort(fe_raw_ptr->getColFieldEntsPtr()->begin(),
730  fe_raw_ptr->getColFieldEntsPtr()->end(), uid_comp);
731  last_col_field_ents_view_size =
732  fe_raw_ptr->getColFieldEntsPtr()->size();
733  }
734 
735  }
736  }
737 
738  auto &dofs_by_ent_uid = dofsField.get<Unique_mi_tag>();
739 
741 }
742 
745 
746  if (verb == DEFAULT_VERBOSITY)
747  verb = verbose;
748 
749  // loop Finite Elements
750  for (auto &fe : finiteElements)
751  CHKERR buildFiniteElements(fe, NULL, verb);
752 
753  if (verb > QUIET) {
754 
755  auto &fe_ents = entsFiniteElements.get<FiniteElement_name_mi_tag>();
756  for (auto &fe : finiteElements) {
757  auto miit = fe_ents.lower_bound(fe->getName());
758  auto hi_miit = fe_ents.upper_bound(fe->getName());
759  const auto count = std::distance(miit, hi_miit);
760  MOFEM_LOG("SYNC", Sev::inform)
761  << "Finite element " << fe->getName()
762  << " added. Nb. of elements added " << count;
763  MOFEM_LOG("SYNC", Sev::noisy) << *fe;
764 
765  auto slg = MoFEM::LogManager::getLog("SYNC");
766  for (auto &field : fIelds) {
767  auto rec = slg.open_record(keywords::severity = Sev::verbose);
768  if (rec) {
769  logging::record_ostream strm(rec);
770  strm << "Field " << field->getName() << " on finite element: ";
771  if ((field->getId() & fe->getBitFieldIdRow()).any())
772  strm << "row ";
773  if ((field->getId() & fe->getBitFieldIdCol()).any())
774  strm << "columns ";
775  if ((field->getId() & fe->getBitFieldIdData()).any())
776  strm << "data";
777  strm.flush();
778  slg.push_record(boost::move(rec));
779  }
780  }
781  }
782 
783  MOFEM_LOG_SYNCHRONISE(cOmm);
784  }
785 
786  *buildMoFEM |= 1 << 1;
788 }
789 
792  SETERRQ(cOmm, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
794 }
795 
796 MoFEMErrorCode Core::build_finite_elements(const string fe_name,
797  const Range *const ents_ptr,
798  int verb) {
800  if (verb == -1)
801  verb = verbose;
802 
803  auto fe_miit = finiteElements.get<FiniteElement_name_mi_tag>().find(fe_name);
804  if (fe_miit == finiteElements.get<FiniteElement_name_mi_tag>().end())
805  SETERRQ1(cOmm, MOFEM_NOT_FOUND, "Finite element <%s> not found",
806  fe_name.c_str());
807 
808  CHKERR buildFiniteElements(*fe_miit, ents_ptr, verb);
809 
810  if (verb >= VERBOSE) {
811  auto &fe_ents = entsFiniteElements.get<FiniteElement_name_mi_tag>();
812  auto miit = fe_ents.lower_bound((*fe_miit)->getName());
813  auto hi_miit = fe_ents.upper_bound((*fe_miit)->getName());
814  const auto count = std::distance(miit, hi_miit);
815  MOFEM_LOG("SYNC", Sev::inform) << "Finite element " << fe_name
816  << " added. Nb. of elements added " << count;
817  MOFEM_LOG_SYNCHRONISE(cOmm);
818  }
819 
820  *buildMoFEM |= 1 << 1;
822 }
823 
824 MoFEMErrorCode Core::build_adjacencies(const Range &ents, int verb) {
826  if (verb == DEFAULT_VERBOSITY)
827  verb = verbose;
828 
829  if (!((*buildMoFEM) & BUILD_FIELD))
830  SETERRQ(cOmm, MOFEM_NOT_FOUND, "field not build");
831  if (!((*buildMoFEM) & BUILD_FE))
832  SETERRQ(cOmm, MOFEM_NOT_FOUND, "fe not build");
833  for (Range::const_pair_iterator peit = ents.pair_begin();
834  peit != ents.pair_end(); ++peit) {
835  EntFiniteElement_multiIndex::index<Ent_mi_tag>::type::iterator fit, hi_fit;
836  fit = entsFiniteElements.get<Ent_mi_tag>().lower_bound(peit->first);
837  hi_fit = entsFiniteElements.get<Ent_mi_tag>().upper_bound(peit->second);
838  for (; fit != hi_fit; ++fit) {
839  if ((*fit)->getBitFieldIdRow().none() &&
840  (*fit)->getBitFieldIdCol().none() &&
841  (*fit)->getBitFieldIdData().none())
842  continue;
843  int by = BYROW;
844  if ((*fit)->getBitFieldIdRow() != (*fit)->getBitFieldIdCol())
845  by |= BYCOL;
846  if ((*fit)->getBitFieldIdRow() != (*fit)->getBitFieldIdData())
847  by |= BYDATA;
849  auto hint = entFEAdjacencies.end();
850  for (auto e : *(*fit)->getRowFieldEntsPtr()) {
851  hint = entFEAdjacencies.emplace_hint(hint, e.lock(), *fit);
852  bool success = entFEAdjacencies.modify(hint, modify_row);
853  if (!success)
854  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
855  "modification unsuccessful");
856  }
857  if ((*fit)->getBitFieldIdRow() != (*fit)->getBitFieldIdCol()) {
858  int by = BYCOL;
859  if ((*fit)->getBitFieldIdCol() != (*fit)->getBitFieldIdData())
860  by |= BYDATA;
862  auto hint = entFEAdjacencies.end();
863  for (auto e : *(*fit)->getColFieldEntsPtr()) {
864  hint = entFEAdjacencies.emplace_hint(hint, e.lock(), *fit);
865  bool success = entFEAdjacencies.modify(hint, modify_col);
866  if (!success)
867  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
868  "modification unsuccessful");
869  }
870  }
871  if ((*fit)->getBitFieldIdRow() != (*fit)->getBitFieldIdData() ||
872  (*fit)->getBitFieldIdCol() != (*fit)->getBitFieldIdData()) {
874  BYDATA);
875  auto hint = entFEAdjacencies.end();
876  for (auto &e : (*fit)->getDataFieldEnts()) {
877  hint = entFEAdjacencies.emplace_hint(hint, e.lock(), *fit);
878  bool success = entFEAdjacencies.modify(hint, modify_data);
879  if (!success)
880  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
881  "modification unsuccessful");
882  }
883  }
884  }
885  }
886 
887  if (verb >= VERBOSE) {
888  MOFEM_LOG("WORLD", Sev::inform)
889  << "Number of adjacencies " << entFEAdjacencies.size();
891  }
892 
893  *buildMoFEM |= 1 << 2;
895 }
896 
898  const BitRefLevel &mask, int verb) {
900  if (verb == -1)
901  verb = verbose;
902  Range ents;
903  CHKERR BitRefManager(*this).getEntitiesByRefLevel(bit, mask, ents);
904 
905  CHKERR build_adjacencies(ents, verb);
906 
908 }
911  if (verb == -1)
912  verb = verbose;
913  CHKERR build_adjacencies(bit, BitRefLevel().set(), verb);
914 
916 }
917 
918 EntFiniteElementByName::iterator
919 Core::get_fe_by_name_begin(const std::string &fe_name) const {
920  return entsFiniteElements.get<FiniteElement_name_mi_tag>().lower_bound(
921  fe_name);
922 }
923 EntFiniteElementByName::iterator
924 Core::get_fe_by_name_end(const std::string &fe_name) const {
925  return entsFiniteElements.get<FiniteElement_name_mi_tag>().upper_bound(
926  fe_name);
927 }
928 
930  const std::string &name) const {
932  FiniteElement_multiIndex::index<FiniteElement_name_mi_tag>::type::iterator it;
933  it = finiteElements.get<FiniteElement_name_mi_tag>().find(name);
934  if (it == finiteElements.get<FiniteElement_name_mi_tag>().end()) {
935  SETERRQ1(cOmm, 1, "finite element not found < %s >", name.c_str());
936  }
937  EntityHandle meshset = (*it)->getMeshset();
938 
939  int num_entities;
940  CHKERR get_moab().get_number_entities_by_handle(meshset, num_entities);
941 
942  if (entsFiniteElements.get<FiniteElement_name_mi_tag>().count(
943  (*it)->getName().c_str()) != (unsigned int)num_entities) {
944  SETERRQ1(cOmm, 1,
945  "not equal number of entities in meshset and finite elements "
946  "multiindex < %s >",
947  (*it)->getName().c_str());
948  }
950 }
951 
954  FiniteElement_multiIndex::index<FiniteElement_name_mi_tag>::type::iterator it;
955  it = finiteElements.get<FiniteElement_name_mi_tag>().begin();
956  for (; it != finiteElements.get<FiniteElement_name_mi_tag>().end(); it++) {
957  EntityHandle meshset = (*it)->getMeshset();
958 
959  int num_entities;
960  CHKERR get_moab().get_number_entities_by_handle(meshset, num_entities);
961 
962  if (entsFiniteElements.get<FiniteElement_name_mi_tag>().count(
963  (*it)->getName().c_str()) != (unsigned int)num_entities) {
964  SETERRQ1(cOmm, 1,
965  "not equal number of entities in meshset and finite elements "
966  "multiindex < %s >",
967  (*it)->getName().c_str());
968  }
969  }
971 }
972 } // namespace MoFEM
EntityHandle get_id_for_max_type()
EntityHandle get_id_for_min_type()
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:282
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:743
UId getLocalUniqueIdCalculate()
Get the Local Unique Id Calculate object.
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:149
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
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:128
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:385
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
static Index< 'p', 3 > p
#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:485
const EntityHandle no_handle
No entity handle is indicated by zero handle, i.e. root meshset.
Definition: Common.hpp:27
#define THROW_MESSAGE(a)
Throw MoFEM exception.
Definition: definitions.h:628
MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)
build adjacencies
Definition: FECore.cpp:824
MoFEMErrorCode list_finite_elements() const
list finite elements in database
Definition: FECore.cpp:306
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
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:924
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:295
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
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:451
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:378
bool check_finite_element(const std::string &name) const
Check if finite element is in database.
Definition: FECore.cpp:31
std::bitset< BITFEID_SIZE > BitFEId
Finite element Id.
Definition: Types.hpp:54
const int dim
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:225
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:919
Finite element definition.
BitFEId getBitFEId(const std::string &name) const
Get field Id.
Definition: FECore.cpp:241
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
Managing BitRefLevels.
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:432
#define FECoreFunctionBegin
Definition: FECore.cpp:21
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:487
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:332
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:405
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:53
#define CHKERR
Inline error check.
Definition: definitions.h:604
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:392
EntityHandle get_finite_element_meshset(const BitFEId id) const
Definition: FECore.cpp:259
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
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:107
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
DEPRECATED MoFEMErrorCode add_ents_to_finite_element_by_PRISMs(const Range &prims, const BitFEId id)
std::string getBitFEIdName(const BitFEId id) const
Get field name.
Definition: FECore.cpp:251
MoFEMErrorCode check_number_of_ents_in_ents_finite_element() const
check data consistency in entsFiniteElements
Definition: FECore.cpp:952
boost::function< MoFEMErrorCode(Interface &moab, const Field &field, const EntFiniteElement &fe, Range &adjacency)> ElementAdjacencyFunct
user adjacency function
MoFEMErrorCode modify_finite_element_off_field_row(const std::string &fe_name, const std::string &name_row)
unset field row which finite element use
Definition: FECore.cpp:207
MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string &name_col)
set field col which finite element use
Definition: FECore.cpp:170
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
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:315
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:190
MoFEMErrorCode buildFiniteElements(const boost::shared_ptr< FiniteElement > &fe, const Range *ents_ptr=NULL, int verb=DEFAULT_VERBOSITY)
Definition: FECore.cpp:505
uint128_t UId
Unique Id.
Definition: Types.hpp:42
MoFEMTypes
Those types control how functions respond on arguments, f.e. error handling.
Definition: definitions.h:189
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:272
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:340