v0.10.0
FieldCore.cpp
Go to the documentation of this file.
1 /** \file FieldCore.cpp
2  * \brief Core interface methods for managing fields.
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 FieldCoreFunctionBegin \
22  MoFEMFunctionBegin; \
23  MOFEM_LOG_CHANNEL("WORLD"); \
24  MOFEM_LOG_CHANNEL("SYNC"); \
25  MOFEM_LOG_FUNCTION(); \
26  MOFEM_LOG_TAG("WORLD", "FieldCore"); \
27  MOFEM_LOG_TAG("SYNC", "FieldCore");
28 
29 namespace MoFEM {
30 
31 BitFieldId Core::getBitFieldId(const std::string &name) const {
32  auto &set = fIelds.get<FieldName_mi_tag>();
33  auto miit = set.find(name);
34  if (miit == set.end()) {
35  THROW_MESSAGE("field < " + name +
36  " > not in database (top tip: check spelling)");
37  }
38  return (*miit)->getId();
39 }
40 
41 FieldBitNumber Core::get_field_bit_number(const std::string name) const {
42  auto &set = fIelds.get<FieldName_mi_tag>();
43  auto miit = set.find(name);
44  if (miit == set.end())
45  THROW_MESSAGE("field not in database (top tip: check spelling)");
46  return (*miit)->getBitNumber();
47 }
48 
50  auto &set = fIelds.get<BitFieldId_mi_tag>();
51  auto miit = set.find(id);
52  if (miit == set.end())
53  THROW_MESSAGE("field not in database (top tip: check spelling)");
54  return (*miit)->meshSet;
55 }
56 
57 EntityHandle Core::get_field_meshset(const std::string name) const {
58  return get_field_meshset(getBitFieldId(name));
59 }
60 
61 bool Core::check_field(const std::string &name) const {
62  auto &set = fIelds.get<FieldName_mi_tag>();
63  auto miit = set.find(name);
64  if (miit == set.end())
65  return false;
66  return true;
67 }
68 
69 const Field *Core::get_field_structure(const std::string &name) {
70  auto &set = fIelds.get<FieldName_mi_tag>();
71  auto miit = set.find(name);
72  if (miit == set.end()) {
73  throw MoFEMException(
75  std::string("field < " + name +
76  " > not in database (top tip: check spelling)")
77  .c_str());
78  }
79  return miit->get();
80 }
81 
83  int dim,
84  Range &ents) const {
85 
87  EntityHandle meshset = get_field_meshset(name);
88  CHKERR get_moab().get_entities_by_dimension(meshset, dim, ents, true);
90 }
91 
93  EntityType type,
94  Range &ents) const {
96  EntityHandle meshset = get_field_meshset(name);
97  CHKERR get_moab().get_entities_by_type(meshset, type, ents, true);
99 }
100 
102  Range &ents) const {
104  EntityHandle meshset = get_field_meshset(name);
105  CHKERR get_moab().get_entities_by_handle(meshset, ents, true);
107 }
108 
109 MoFEMErrorCode Core::addField(const std::string &name, const FieldSpace space,
110  const FieldApproximationBase base,
111  const FieldCoefficientsNumber nb_of_coefficients,
112  const TagType tag_type, const enum MoFEMTypes bh,
113  int verb) {
114  MOFEM_LOG_CHANNEL("WORLD");
115  MOFEM_LOG_TAG("WORLD", "FieldCore");
116  MOFEM_LOG_CHANNEL("SYNC");
117  MOFEM_LOG_TAG("SYNC", "FieldCore");
120 
121  if (verb == -1)
122  verb = verbose;
123  *buildMoFEM = 0;
124  auto fit = fIelds.get<FieldName_mi_tag>().find(name);
125  if (fit != fIelds.get<FieldName_mi_tag>().end()) {
126  if (bh == MF_EXCL)
127  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
128  "field is <%s> in database", name.c_str());
129 
130  } else {
131 
132  EntityHandle meshset;
133  CHKERR get_moab().create_meshset(MESHSET_SET, meshset);
134 
135  // Add field mesh set to partion meshset. In case of no elements
136  // on processor part, when mesh file is read, finite element meshset is
137  // prevented from deletion by moab reader.
138  auto add_meshset_to_partition = [&](auto meshset) {
140  const void *tag_vals[] = {&rAnk};
141  ParallelComm *pcomm = ParallelComm::get_pcomm(
142  &get_moab(), get_basic_entity_data_ptr()->pcommID);
143  Tag part_tag = pcomm->part_tag();
144  Range tagged_sets;
145  CHKERR get_moab().get_entities_by_type_and_tag(0, MBENTITYSET, &part_tag,
146  tag_vals, 1, tagged_sets,
147  moab::Interface::UNION);
148  for (auto s : tagged_sets)
149  CHKERR get_moab().add_entities(s, &meshset, 1);
151  };
152  CHKERR add_meshset_to_partition(meshset);
153 
154  // id
155  BitFieldId id = getFieldShift();
156 
157  auto create_tags = [&]() {
159  CHKERR
160  get_moab().tag_set_data(th_FieldId, &meshset, 1, &id);
161  // space
162  CHKERR get_moab().tag_set_data(th_FieldSpace, &meshset, 1, &space);
163  // base
164  CHKERR get_moab().tag_set_data(th_FieldBase, &meshset, 1, &base);
165 
166  // name
167  void const *tag_data[] = {name.c_str()};
168  int tag_sizes[1];
169  tag_sizes[0] = name.size();
170  CHKERR get_moab().tag_set_by_ptr(th_FieldName, &meshset, 1, tag_data,
171  tag_sizes);
172  // name data prefix
173  const std::string name_data_prefix("_App_Data");
174  void const *tag_prefix_data[] = {name_data_prefix.c_str()};
175  int tag_prefix_sizes[1];
176  tag_prefix_sizes[0] = name_data_prefix.size();
177  CHKERR get_moab().tag_set_by_ptr(th_FieldName_DataNamePrefix, &meshset, 1,
178  tag_prefix_data, tag_prefix_sizes);
179  Tag th_app_order, th_field_data, th_field_data_vert, th_rank;
180  // data
181  std::string tag_data_name = name_data_prefix + name;
182  const int def_len = 0;
183  CHKERR get_moab().tag_get_handle(
184  tag_data_name.c_str(), def_len, MB_TYPE_DOUBLE, th_field_data,
185  MB_TAG_CREAT | MB_TAG_VARLEN | MB_TAG_SPARSE, NULL);
186  std::string tag_data_name_verts = name_data_prefix + name + "V";
187  VectorDouble def_vert_data(nb_of_coefficients);
188  def_vert_data.clear();
189  CHKERR get_moab().tag_get_handle(
190  tag_data_name_verts.c_str(), nb_of_coefficients, MB_TYPE_DOUBLE,
191  th_field_data_vert, MB_TAG_CREAT | tag_type, &*def_vert_data.begin());
192  // order
193  ApproximationOrder def_ApproximationOrder = -1;
194  const std::string Tag_ApproximationOrder_name = "_App_Order_" + name;
195  CHKERR get_moab().tag_get_handle(
196  Tag_ApproximationOrder_name.c_str(), 1, MB_TYPE_INTEGER, th_app_order,
197  MB_TAG_CREAT | tag_type, &def_ApproximationOrder);
198  // rank
199  int def_rank = 1;
200  const std::string tag_rank_name = "_Field_Rank_" + name;
201  CHKERR get_moab().tag_get_handle(tag_rank_name.c_str(), 1,
202  MB_TYPE_INTEGER, th_rank,
203  MB_TAG_CREAT | MB_TAG_SPARSE, &def_rank);
204  CHKERR get_moab().tag_set_data(th_rank, &meshset, 1, &nb_of_coefficients);
205 
207  };
208 
209  auto create_undefined_cs = [&](auto &undefined_cs_ptr) {
211  const int unidefined_cs_dim[] = {-1, 0, 0, 0};
212  CHKERR getInterface<CoordSystemsManager>()->addCoordinateSystem(
213  unidefined_cs_dim, "UNDEFINED", MF_ZERO);
214  CHKERR getInterface<CoordSystemsManager>()->getCoordSysPtr(
215  "UNDEFINED", undefined_cs_ptr);
217  };
218 
219  auto add_field_meshset_to_cs = [&](auto &undefined_cs_ptr) {
221  int cs_name_size[1];
222  cs_name_size[0] = undefined_cs_ptr->getName().size();
223  void const *cs_name[] = {&*undefined_cs_ptr->getNameRef().begin()};
224  CHKERR get_moab().tag_set_by_ptr(
225  getInterface<CoordSystemsManager>()->get_th_CoordSysName(), &meshset,
226  1, cs_name, cs_name_size);
227  EntityHandle coord_sys_meshset = undefined_cs_ptr->getMeshset();
228  CHKERR get_moab().add_entities(coord_sys_meshset, &meshset, 1);
230  };
231 
232  CHKERR create_tags();
233  boost::shared_ptr<CoordSys> undefined_cs_ptr;
234  CHKERR create_undefined_cs(undefined_cs_ptr);
235  CHKERR add_field_meshset_to_cs(undefined_cs_ptr);
236 
237  auto p = fIelds.insert(
238  boost::make_shared<Field>(moab, meshset, undefined_cs_ptr));
239  if (verb > QUIET) {
240  MOFEM_LOG("WORLD", Sev::inform) << "Add field " << **p.first;
241  MOFEM_LOG("WORLD", Sev::noisy)
242  << "Field " << (*p.first)->getName() << " core value < "
243  << this->getValue() << " > field value ) "
244  << (*p.first)->getBitNumber() << " )";
245  }
246 
247  if (!p.second)
248  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND,
249  "field not inserted %s (top tip, it could be already "
250  "there)",
251  Field(moab, meshset, undefined_cs_ptr).getName().c_str());
252  }
253 
255 }
256 
257 MoFEMErrorCode Core::add_field(const std::string &name, const FieldSpace space,
258  const FieldApproximationBase base,
259  const FieldCoefficientsNumber nb_of_coefficients,
260  const TagType tag_type, const enum MoFEMTypes bh,
261  int verb) {
262  return this->addField(name, space, base, nb_of_coefficients, tag_type, bh,
263  verb);
264 }
265 
266 MoFEMErrorCode Core::addEntsToFieldByDim(const Range &ents, const int dim,
267  const std::string &name, int verb) {
268 
269  *buildMoFEM = 0;
270  EntityHandle idm = no_handle;
271  if (verb == -1)
272  verb = verbose;
274  idm = get_field_meshset(name);
275  FieldSpace space;
276  CHKERR get_moab().tag_get_data(th_FieldSpace, &idm, 1, &space);
277  std::vector<int> nb_ents_on_dim(3, 0);
278  switch (space) {
279  case L2:
280  CHKERR get_moab().add_entities(idm, ents);
281  if (verb >= VERY_VERBOSE) {
282  std::ostringstream ss;
283  ss << "add entities to field " << name;
284  ss << " nb. add ents " << ents.size();
285  ss << std::endl;
286  PetscSynchronizedPrintf(cOmm, ss.str().c_str());
287  }
288  break;
289  case H1:
290  CHKERR get_moab().add_entities(idm, ents);
291  for (int dd = 0; dd != dim; ++dd) {
292  Range adj_ents;
293  CHKERR get_moab().get_adjacencies(ents, dd, false, adj_ents,
294  moab::Interface::UNION);
295  if (dd == 0) {
296  Range topo_nodes;
297  CHKERR get_moab().get_connectivity(ents, topo_nodes, true);
298  Range mid_nodes;
299  CHKERR get_moab().get_connectivity(ents, mid_nodes, false);
300  mid_nodes = subtract(mid_nodes, topo_nodes);
301  adj_ents = subtract(adj_ents, mid_nodes);
302  }
303  CHKERR get_moab().add_entities(idm, adj_ents);
304  nb_ents_on_dim[dd] = adj_ents.size();
305  }
306  break;
307  case HCURL:
308  CHKERR get_moab().add_entities(idm, ents);
309  for (int dd = 1; dd != dim; ++dd) {
310  Range adj_ents;
311  CHKERR get_moab().get_adjacencies(ents, dd, false, adj_ents,
312  moab::Interface::UNION);
313  CHKERR get_moab().add_entities(idm, adj_ents);
314  nb_ents_on_dim[dd] = adj_ents.size();
315  }
316  break;
317  case HDIV:
318  CHKERR get_moab().add_entities(idm, ents);
319  if (dim > 2) {
320  Range adj_ents;
321  CHKERR get_moab().get_adjacencies(ents, 2, false, adj_ents,
322  moab::Interface::UNION);
323  CHKERR get_moab().add_entities(idm, adj_ents);
324  nb_ents_on_dim[2] = adj_ents.size();
325  }
326  break;
327  default:
328  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
329  "sorry, unknown space added to entity");
330  }
331  if (verb >= VERY_VERBOSE) {
332  std::ostringstream ss;
333  ss << "add entities to field " << name;
334  ss << " nb. add ents " << ents.size();
335  ss << " nb. add faces " << nb_ents_on_dim[2];
336  ss << " nb. add edges " << nb_ents_on_dim[1];
337  ss << " nb. add nodes " << nb_ents_on_dim[0];
338  ss << std::endl;
339  PetscSynchronizedPrintf(cOmm, ss.str().c_str());
340  PetscSynchronizedFlush(cOmm, PETSC_STDOUT);
341  }
343 }
344 
345 MoFEMErrorCode Core::add_ents_to_field_by_dim(const Range &ents, const int dim,
346  const std::string &name,
347  int verb) {
348  Range ents_dim = ents.subset_by_dimension(dim);
349  return addEntsToFieldByDim(ents_dim, dim, name, verb);
350 }
351 
353  const EntityType type,
354  const std::string &name,
355  int verb) {
357  Range ents_type = ents.subset_by_type(type);
358  if (!ents_type.empty()) {
359  const int dim = get_moab().dimension_from_handle(ents_type[0]);
360  CHKERR addEntsToFieldByDim(ents_type, dim, name, verb);
361  }
363 }
364 
366  const int dim,
367  const std::string &name,
368  const bool recursive, int verb) {
370  Range ents;
371  CHKERR get_moab().get_entities_by_dimension(meshset, dim, ents, recursive);
372  CHKERR addEntsToFieldByDim(ents, dim, name, verb);
374 }
375 
377  const EntityType type,
378  const std::string &name,
379  const bool recursive, int verb) {
381  Range ents;
382  CHKERR get_moab().get_entities_by_type(meshset, type, ents, recursive);
383  if (!ents.empty()) {
384  const int dim = get_moab().dimension_from_handle(ents[0]);
385  CHKERR addEntsToFieldByDim(ents, dim, name, verb);
386  }
388 }
389 
391  const double coords[],
392  int size, int verb) {
394 
395  if (verb == DEFAULT_VERBOSITY)
396  verb = verbose;
397 
398  Range verts;
399 
400  auto create_vertices = [&]() {
402 
403  vector<double *> arrays_coord;
404  EntityHandle startv = 0;
405  ReadUtilIface *iface;
406  CHKERR get_moab().query_interface(iface);
407  CHKERR iface->get_node_coords(3, size, 0, startv, arrays_coord);
408  verts.insert(startv, startv + size - 1);
409  for (int n = 0; n != size; ++n)
410  for (auto d : {0, 1, 2})
411  arrays_coord[d][n] = coords[3 * n + d];
412 
414  };
415 
416  auto add_verts_to_field = [&]() {
418  EntityHandle field_meshset = get_field_meshset(name);
419  CHKERR get_moab().add_entities(field_meshset, verts);
421  };
422 
423  CHKERR create_vertices();
424  CHKERR add_verts_to_field();
425 
427 }
428 
429 MoFEMErrorCode Core::setFieldOrderImpl(boost::shared_ptr<Field> field_ptr,
430  const Range &ents,
432  int verb) {
434 
435  if (verb == DEFAULT_VERBOSITY)
436  verb = verbose;
437  *buildMoFEM = 0;
438 
439  MOFEM_LOG("WORLD", Sev::noisy) << "Test field " << *field_ptr;
440 
441  const auto field_meshset = field_ptr->getMeshset();
442  const auto bit_number = field_ptr->getBitNumber();
443 
444  // intersection with field meshset
445  Range ents_of_id_meshset;
446  CHKERR get_moab().get_entities_by_handle(field_meshset, ents_of_id_meshset,
447  false);
448  Range field_ents = intersect(ents, ents_of_id_meshset);
449  if (verb > QUIET)
450  MOFEM_LOG_C("SYNC", Sev::noisy,
451  "change nb. of ents for order in the field <%s> %d",
452  field_ptr->getName().c_str(), field_ents.size());
453 
454  // ent view by field id (in set all MoabEnts has the same FieldId)
455  auto eiit = entsFields.get<Unique_mi_tag>().lower_bound(
456  FieldEntity::getLoBitNumberUId(field_ptr->getBitNumber()));
457  FieldEntity_multiIndex_ent_view ents_id_view;
458  if (eiit != entsFields.get<Unique_mi_tag>().end()) {
459  auto hi_eiit = entsFields.get<Unique_mi_tag>().upper_bound(
460  FieldEntity::getHiBitNumberUId(field_ptr->getBitNumber()));
461  std::copy(eiit, hi_eiit, std::back_inserter(ents_id_view));
462  }
463 
464  if (verb > QUIET)
465  MOFEM_LOG_C("SYNC", Sev::noisy,
466  "current nb. of ents in the multi index field <%s> %d",
467  field_ptr->getName().c_str(), ents_id_view.size());
468 
469  // loop over ents
470  int nb_ents_set_order_up = 0;
471  int nb_ents_set_order_down = 0;
472  int nb_ents_set_order_new = 0;
473 
474  FieldEntity_change_order modify_order_no_size_change(order, false);
475  FieldEntity_change_order modify_order_size_change(order, true);
476 
477  for (auto pit = field_ents.const_pair_begin();
478  pit != field_ents.const_pair_end(); pit++) {
479  EntityHandle first = pit->first;
480  EntityHandle second = pit->second;
481 
482  // Sanity check
483  switch (field_ptr->getSpace()) {
484  case H1:
485  break;
486  case HCURL:
487  if (get_moab().type_from_handle(first) == MBVERTEX)
488  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
489  "Hcurl space on vertices makes no sense");
490 
491  break;
492  case HDIV:
493  if (get_moab().type_from_handle(first) == MBVERTEX)
494  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
495  "Hdiv space on vertices makes no sense");
496 
497  if (get_moab().type_from_handle(first) == MBEDGE)
498  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
499  "Hdiv space on edges makes no sense");
500 
501  break;
502  default:
503  break;
504  }
505 
506  // Entity is in database, change order only if needed
507  Range ents_in_database;
508  auto vit = ents_id_view.get<1>().lower_bound(first);
509  auto hi_vit = ents_id_view.get<1>().upper_bound(second);
510  if (order >= 0) {
511  for (; vit != hi_vit; ++vit) {
512  ents_in_database.insert(vit->get()->getEnt());
513  // entity is in database and order is changed or reset
514  const ApproximationOrder old_approximation_order =
515  (*vit)->getMaxOrder();
516 
517  if (old_approximation_order != order) {
518 
519  FieldEntity_multiIndex::iterator miit =
520  entsFields.get<Unique_mi_tag>().find((*vit)->getLocalUniqueId());
521 
522  if ((*miit)->getMaxOrder() < order)
523  nb_ents_set_order_up++;
524  if ((*miit)->getMaxOrder() > order)
525  nb_ents_set_order_down++;
526 
527  // set dofs inactive if order is reduced, and set new order to entity
528  // if order is increased (note that dofs are not build if order is
529  // increased)
530 
531  bool can_change_size = true;
532  auto dit = dofsField.get<Unique_mi_tag>().lower_bound(
534  (*miit)->getEnt()));
535  if (dit != dofsField.get<Unique_mi_tag>().end()) {
536  auto hi_dit = dofsField.get<Unique_mi_tag>().upper_bound(
538  (*miit)->getEnt()));
539 
540  if (dit != hi_dit)
541  can_change_size = false;
542  for (; dit != hi_dit; dit++) {
543  if ((*dit)->getDofOrder() > order) {
544  bool success = dofsField.modify(dofsField.project<0>(dit),
545  DofEntity_active_change(false));
546  if (!success)
547  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
548  "modification unsuccessful");
549  }
550  }
551  }
552 
553  bool success =
554  entsFields.modify(entsFields.project<0>(miit),
555  can_change_size ? modify_order_size_change
556  : modify_order_no_size_change);
557  if (!success)
558  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
559  "modification unsuccessful");
560  }
561  }
562  }
563 
564  Range new_ents = subtract(Range(first, second), ents_in_database);
565  for (Range::const_pair_iterator pit = new_ents.const_pair_begin();
566  pit != new_ents.const_pair_end(); ++pit) {
567  EntityHandle first = pit->first;
568  EntityHandle second = pit->second;
569  const EntityType ent_type = get_moab().type_from_handle(first);
570  auto get_nb_dofs_on_order = [&](const int order) {
571  return order >= 0 ? (field_ptr->getFieldOrderTable()[ent_type])(order)
572  : 0;
573  };
574  const int nb_dofs_on_order = get_nb_dofs_on_order(order);
575  if (nb_dofs_on_order || order == -1) {
576 
577  const int field_rank = field_ptr->getNbOfCoeffs();
578  const int nb_dofs = nb_dofs_on_order * field_rank;
579 
580  // Entity is not in database and order is changed or reset
581  auto miit_ref_ent =
582  refinedEntities.get<Ent_mi_tag>().lower_bound(first);
583 
584  auto create_tags_for_max_order = [&](const Range &ents) {
586  if (order >= 0) {
587  std::vector<ApproximationOrder> o_vec(ents.size(), order);
588  CHKERR get_moab().tag_set_data(field_ptr->th_AppOrder, ents,
589  &*o_vec.begin());
590  }
592  };
593 
594  auto create_tags_for_data = [&](const Range &ents) {
596  if (order >= 0) {
597 
598  if (nb_dofs > 0) {
599  if (ent_type == MBVERTEX) {
600  std::vector<FieldData> d_vec(nb_dofs * ents.size(), 0);
601  CHKERR get_moab().tag_set_data(field_ptr->th_FieldDataVerts,
602  ents, &*d_vec.begin());
603  } else {
604  std::vector<int> tag_size(ents.size(), nb_dofs);
605  std::vector<FieldData> d_vec(nb_dofs, 0);
606  std::vector<void const *> d_vec_ptr(ents.size(),
607  &*d_vec.begin());
608  CHKERR get_moab().tag_set_by_ptr(field_ptr->th_FieldData, ents,
609  &*d_vec_ptr.begin(),
610  &*tag_size.begin());
611  }
612  }
613  }
614 
616  };
617 
618  auto get_ents_in_ref_ent = [&](auto miit_ref_ent) {
619  auto hi = refinedEntities.get<Ent_mi_tag>().upper_bound(second);
620  Range in;
621  for (; miit_ref_ent != hi; ++miit_ref_ent)
622  in.insert(miit_ref_ent->get()->getEnt());
623  return in;
624  };
625 
626  auto get_ents_max_order = [&](const Range &ents) {
627  boost::shared_ptr<std::vector<const void *>> vec(
628  new std::vector<const void *>());
629  vec->resize(ents.size());
630  CHKERR get_moab().tag_get_by_ptr(field_ptr->th_AppOrder, ents,
631  &*vec->begin());
632  return vec;
633  };
634 
635  auto get_ents_field_data_vector_adaptor =
636  [&](const Range &ents,
637  boost::shared_ptr<std::vector<const void *>> &ents_max_orders) {
638  // create shared pointer and reserve memory
639  boost::shared_ptr<std::vector<double *>> vec(
640  new std::vector<double *>());
641  vec->reserve(ents.size());
642 
643  if (order >= 0 && nb_dofs == 0) {
644  // set empty vector adaptor
645  for (int i = 0; i != ents.size(); ++i)
646  vec->emplace_back(nullptr);
647  } else {
648  moab::ErrorCode rval;
649  std::vector<int> tag_size(ents.size());
650  std::vector<const void *> d_vec_ptr(ents.size());
651 
652  // get tags data
653  if (ent_type == MBVERTEX)
654  rval = get_moab().tag_get_by_ptr(field_ptr->th_FieldDataVerts,
655  ents, &*d_vec_ptr.begin(),
656  &*tag_size.begin());
657  else
658  rval = get_moab().tag_get_by_ptr(field_ptr->th_FieldData,
659  ents, &*d_vec_ptr.begin(),
660  &*tag_size.begin());
661 
662  auto cast = [](auto p) {
663  return const_cast<FieldData *const>(
664  static_cast<const FieldData *>(p));
665  };
666 
667  // some of entities has tag not set or zero dofs on entity
668  if (rval == MB_SUCCESS) {
669  // all is ok, all entities has tag set
670  auto tit = d_vec_ptr.begin();
671  auto oit = ents_max_orders->begin();
672  for (auto sit = tag_size.begin(); sit != tag_size.end();
673  ++sit, ++tit, ++oit)
674  vec->emplace_back(cast(*tit));
675 
676  } else {
677  // set empty vector adaptor
678  for (int i = 0; i != ents.size(); ++i)
679  vec->emplace_back(nullptr);
680 
681  // check order on all entities, and if for that order non zero
682  // dofs are expected get pointer to tag data and reset vector
683  // adaptor
684  auto oit = ents_max_orders->begin();
685  auto dit = vec->begin();
686  for (auto eit = ents.begin(); eit != ents.end();
687  ++eit, ++oit, ++dit) {
688 
689  const int ent_order =
690  *static_cast<const ApproximationOrder *>(*oit);
691  const int ent_nb_dofs = get_nb_dofs_on_order(ent_order);
692 
693  if (ent_nb_dofs) {
694  int tag_size;
695  const void *ret_val;
696  if (ent_type == MBVERTEX) {
697  rval = get_moab().tag_get_by_ptr(
698  field_ptr->th_FieldDataVerts, &*eit, 1, &ret_val,
699  &tag_size);
700  } else {
701  rval = get_moab().tag_get_by_ptr(
702  field_ptr->th_FieldData, &*eit, 1, &ret_val,
703  &tag_size);
704  if (rval != MB_SUCCESS) {
705 
706  const int set_tag_size[] = {ent_nb_dofs * field_rank};
707  std::array<FieldData, MAX_DOFS_ON_ENTITY> set_d_vec;
708  std::fill(set_d_vec.begin(),
709  &set_d_vec[set_tag_size[0]], 0);
710  const void *set_d_vec_ptr[] = {set_d_vec.data()};
711  CHKERR get_moab().tag_set_by_ptr(
712  field_ptr->th_FieldData, &*eit, 1, set_d_vec_ptr,
713  set_tag_size);
714  rval = get_moab().tag_get_by_ptr(
715  field_ptr->th_FieldData, &*eit, 1, &ret_val,
716  &tag_size);
717 
718  if (rval != MB_SUCCESS) {
719  MOFEM_LOG_ATTRIBUTES("SELF",
722  MOFEM_LOG("SELF", Sev::error)
723  << "Error is triggered in MOAB, field tag data "
724  "for same reason can not be for accessed.";
725  MOFEM_LOG("SELF", Sev::error)
726  << "Set order: " << order;
727  MOFEM_LOG("SELF", Sev::error)
728  << "Nb. dofs on entity for given order: "
729  << set_tag_size[0];
730  MOFEM_LOG("SELF", Sev::error)
731  << "Entity type: "
732  << moab::CN::EntityTypeName(ent_type);
733  MOFEM_LOG("SELF", Sev::error)
734  << "Field: " << *field_ptr;
735  CHKERR rval;
736  }
737  }
738  }
739  const_cast<FieldData *&>(*dit) = cast(ret_val);
740  }
741  }
742  }
743  }
744  return vec;
745  };
746 
747  auto ents_in_ref_ent = get_ents_in_ref_ent(miit_ref_ent);
748 
749  CHKERR create_tags_for_max_order(ents_in_ref_ent);
750  CHKERR create_tags_for_data(ents_in_ref_ent);
751  auto ents_max_order = get_ents_max_order(ents_in_ref_ent);
752  auto ent_field_data =
753  get_ents_field_data_vector_adaptor(ents_in_ref_ent, ents_max_order);
754 
755  // reserve memory for field dofs
756  auto ents_array = boost::make_shared<std::vector<FieldEntity>>();
757  // Add sequence to field data structure. Note that entities are
758  // allocated once into vector. This vector is passed into sequence as a
759  // weak_ptr. Vector is destroyed at the point last entity inside that
760  // vector is destroyed.
761  ents_array->reserve(second - first + 1);
762  auto vit_max_order = ents_max_order->begin();
763  auto vit_field_data = ent_field_data->begin();
764  for (auto ent : ents_in_ref_ent) {
765  ents_array->emplace_back(
766  field_ptr, *miit_ref_ent,
767  boost::shared_ptr<double *const>(ent_field_data,
768  &*vit_field_data),
769  boost::shared_ptr<const int>(
770  ents_max_order, static_cast<const int *>(*vit_max_order)));
771  ++vit_max_order;
772  ++vit_field_data;
773  ++miit_ref_ent;
774  }
775  if (!ents_array->empty())
776  if ((*ents_array)[0].getFieldRawPtr() != field_ptr.get())
777  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
778  "Get field ent poiter and field pointer do not match for "
779  "field %s",
780  field_ptr->getName().c_str());
781  nb_ents_set_order_new += ents_array->size();
782 
783  // Check if any of entities in the range has bit level but is not added
784  // to database. That generate data inconsistency and error.
785  if (ents_in_ref_ent.size() < (second - first + 1)) {
786  Range ents_not_in_database =
787  subtract(Range(first, second), ents_in_ref_ent);
788  std::vector<const void *> vec_bits(ents_not_in_database.size());
789  CHKERR get_moab().tag_get_by_ptr(
790  get_basic_entity_data_ptr()->th_RefBitLevel, ents_not_in_database,
791  &*vec_bits.begin());
792  auto cast = [](auto p) {
793  return static_cast<const BitRefLevel *>(p);
794  };
795  for (auto v : vec_bits)
796  if (cast(v)->any())
797  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
798  "Try to add entities which are not seeded or added to "
799  "database");
800  }
801 
802  // Add entities to database
803  auto hint = entsFields.end();
804  for (auto &v : *ents_array)
805  hint = entsFields.emplace_hint(hint, ents_array, &v);
806  }
807  }
808  }
809 
810  if (verb > QUIET) {
811  MOFEM_LOG_C("SYNC", Sev::noisy,
812  "nb. of entities in field <%s> for which order was "
813  "increased %d (order %d)",
814  field_ptr->getName().c_str(), nb_ents_set_order_up, order);
815  MOFEM_LOG_C("SYNC", Sev::noisy,
816  "nb. of entities in field <%s> for which order was "
817  "reduced %d (order %d)",
818  field_ptr->getName().c_str(), nb_ents_set_order_down, order);
819  MOFEM_LOG_C(
820  "SYNC", Sev::noisy,
821  "nb. of entities in field <%s> for which order set %d (order %d)",
822  field_ptr->getName().c_str(), nb_ents_set_order_new, order);
823  MOFEM_LOG_SYNCHRONISE(cOmm);
824  }
825 
826  if (verb > QUIET) {
827  auto eiit = entsFields.get<Unique_mi_tag>().lower_bound(
828  FieldEntity::getLoBitNumberUId(field_ptr->getBitNumber()));
829  auto hi_eiit = entsFields.get<Unique_mi_tag>().upper_bound(
830  FieldEntity::getHiBitNumberUId(field_ptr->getBitNumber()));
831  MOFEM_LOG_C("SYNC", Sev::noisy,
832  "nb. of ents in the multi index field <%s> %d",
833  field_ptr->getName().c_str(), std::distance(eiit, hi_eiit));
834  }
835 
837 }
838 
839 MoFEMErrorCode Core::set_field_order(const Range &ents, const BitFieldId id,
840  const ApproximationOrder order, int verb) {
841  MOFEM_LOG_CHANNEL("WORLD");
842  MOFEM_LOG_TAG("WORLD", "FieldCore");
845 
846  // check field & meshset
847  auto field_it = fIelds.get<BitFieldId_mi_tag>().find(id);
848  if (field_it == fIelds.get<BitFieldId_mi_tag>().end())
849  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "no filed found");
850 
851  MOFEM_LOG("WORLD", Sev::noisy)
852  << "Field " << (*field_it)->getName() << " core value < "
853  << this->getValue() << " > field value ( " << (*field_it)->getBitNumber()
854  << " )";
855 
856  CHKERR this->setFieldOrderImpl(*field_it, ents, order, verb);
857 
859 }
860 
862  const EntityType type, const BitFieldId id,
863  const ApproximationOrder order, int verb) {
865  if (verb == -1)
866  verb = verbose;
867  *buildMoFEM = 0;
868  Range ents;
869  CHKERR get_moab().get_entities_by_type(meshset, type, ents);
870  if (verb > VERBOSE) {
871  PetscSynchronizedPrintf(cOmm, "nb. of ents for order change %d\n",
872  ents.size());
873  }
874  CHKERR this->set_field_order(ents, id, order, verb);
875  if (verb > VERBOSE) {
876  PetscSynchronizedFlush(cOmm, PETSC_STDOUT);
877  }
879 }
881  const EntityType type,
882  const std::string &name,
883  const ApproximationOrder order, int verb) {
885  if (verb == -1)
886  verb = verbose;
887  *buildMoFEM = 0;
888  CHKERR this->set_field_order(meshset, type, getBitFieldId(name), order, verb);
890 }
891 MoFEMErrorCode Core::set_field_order(const Range &ents, const std::string &name,
892  const ApproximationOrder order, int verb) {
894  if (verb == -1)
895  verb = verbose;
896  *buildMoFEM = 0;
897  CHKERR this->set_field_order(ents, getBitFieldId(name), order, verb);
899 }
901  const BitRefLevel &bit, const BitRefLevel &mask, const EntityType type,
902  const BitFieldId id, const ApproximationOrder order, int verb) {
904  if (verb == -1)
905  verb = verbose;
906  *buildMoFEM = 0;
907  Range ents;
909  ents, verb);
910  CHKERR this->set_field_order(ents, id, order, verb);
912 }
913 
915  const BitRefLevel &bit, const BitRefLevel &mask, const EntityType type,
916  const std::string &name, const ApproximationOrder order, int verb) {
918  if (verb == -1)
919  verb = verbose;
920  *buildMoFEM = 0;
921  Range ents;
923  ents, verb);
924  CHKERR this->set_field_order(ents, getBitFieldId(name), order, verb);
926 }
927 
929 Core::buildFieldForNoFieldImpl(boost::shared_ptr<Field> field_ptr,
930  std::map<EntityType, int> &dof_counter,
931  int verb) {
933 
934  const auto bit_number = field_ptr->getBitNumber();
935 
936  // ents in the field meshset
937  Range ents;
938  CHKERR get_moab().get_entities_by_handle(field_ptr->getMeshset(), ents,
939  false);
940  if (verb > VERBOSE)
941  MOFEM_LOG_C("SYNC", Sev::noisy, "Ents in field %s meshset %d\n",
942  field_ptr->getName().c_str(), ents.size());
943 
944  // ent view by field id (in set all MoabEnts has the same FieldId)
945  auto eiit = entsFields.get<Unique_mi_tag>().lower_bound(
946  FieldEntity::getLoBitNumberUId(field_ptr->getBitNumber()));
947  FieldEntity_multiIndex_ent_view ents_id_view;
948  if (eiit != entsFields.get<Unique_mi_tag>().end()) {
949  auto hi_eiit = entsFields.get<Unique_mi_tag>().upper_bound(
950  FieldEntity::getHiBitNumberUId(field_ptr->getBitNumber()));
951  std::copy(eiit, hi_eiit, std::back_inserter(ents_id_view));
952  }
953 
954  boost::shared_ptr<const int> zero_order(new const int(0));
955 
956  for (auto ent : ents) {
957  // search if field meshset is in database
958  auto ref_ent_it = refinedEntities.get<Ent_mi_tag>().find(ent);
959  if (ref_ent_it == refinedEntities.get<Ent_mi_tag>().end())
960  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
961  "Entity is not in MoFEM database, entities in field meshset need "
962  "to be seeded (i.e. bit ref level add to them)");
963 
964  auto add_dofs = [&](auto field_eit) {
966  // create dofs on this entity (nb. of dofs is equal to rank)
967  for (FieldCoefficientsNumber rank = 0; rank < field_ptr->getNbOfCoeffs();
968  rank++) {
969  // insert dof
970  auto p = dofsField.insert(
971  boost::make_shared<DofEntity>(field_eit, 0, rank, rank));
972  if (p.second) {
973  dof_counter[MBENTITYSET]++; // Count entities in the meshset
974  } else
975  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
976  "Dof expected to be created");
977  }
979  };
980 
981  // create database entity
982  auto field_ent_it = ents_id_view.get<1>().find(ent);
983  if (field_ent_it == ents_id_view.get<1>().end()) {
984 
985  auto p = entsFields.insert(
986 
987  boost::make_shared<FieldEntity>(
988  field_ptr, *ref_ent_it,
990  *ref_ent_it),
991  boost::shared_ptr<const int>(zero_order, zero_order.get()))
992 
993  );
994 
995  if ((*p.first)->getFieldRawPtr() != field_ptr.get())
996  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
997  "Get field ent poiter and field pointer do not match for "
998  "field %s",
999  field_ptr->getName().c_str());
1000 
1001  if (!p.second)
1002  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1003  "Entity should be created");
1004 
1005  CHKERR add_dofs(*(p.first));
1006 
1007  } else {
1008 
1009  // If there are DOFs in that range is more pragmatic to remove them
1010  // rather than to find sub-ranges or make them inactive
1011  auto dit = dofsField.get<Unique_mi_tag>().lower_bound(
1012  FieldEntity::getLoLocalEntityBitNumber(bit_number, ent));
1013  auto hi_dit = dofsField.get<Unique_mi_tag>().upper_bound(
1014  FieldEntity::getHiLocalEntityBitNumber(bit_number, ent));
1015  dofsField.get<Unique_mi_tag>().erase(dit, hi_dit);
1016  CHKERR add_dofs(*field_ent_it);
1017  }
1018  }
1019 
1020  if (verb > VERBOSE) {
1021  auto lo_dof = dofsField.get<Unique_mi_tag>().lower_bound(
1022  FieldEntity::getLoBitNumberUId(field_ptr->getBitNumber()));
1023  auto hi_dof = dofsField.get<Unique_mi_tag>().upper_bound(
1024  FieldEntity::getHiBitNumberUId(field_ptr->getBitNumber()));
1025  for (; lo_dof != hi_dof; lo_dof++)
1026  MOFEM_LOG("SYNC", Sev::noisy) << **lo_dof;
1027  MOFEM_LOG_SYNCHRONISE(cOmm);
1028  }
1029 
1031 }
1032 
1035  std::map<EntityType, int> &dof_counter, int verb) {
1037 
1038  if (verb == -1)
1039  verb = verbose;
1040 
1041  // find fields
1042  auto field_it = fIelds.get<BitFieldId_mi_tag>().find(id);
1043  if (field_it == fIelds.get<BitFieldId_mi_tag>().end())
1044  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "Field not found");
1045 
1046  if (verb > QUIET)
1047  MOFEM_LOG("WORLD", Sev::noisy)
1048  << "Field " << (*field_it)->getName() << " core value < "
1049  << this->getValue() << " > field value () "
1050  << (*field_it)->getBitNumber() << " )";
1051 
1052  CHKERR this->buildFieldForNoFieldImpl(*field_it, dof_counter, verb);
1053 
1055 }
1056 
1058  const BitFieldId id, std::map<EntityType, int> &dof_counter,
1059  std::map<EntityType, int> &inactive_dof_counter, int verb) {
1061  if (verb == -1)
1062  verb = verbose;
1063 
1064  // Find field
1065  auto &set_id = fIelds.get<BitFieldId_mi_tag>();
1066  auto field_it = set_id.find(id);
1067  if (field_it == set_id.end()) {
1068  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "Field not found");
1069  }
1070  const int bit_number = field_it->get()->getBitNumber();
1071  const int rank = field_it->get()->getNbOfCoeffs();
1072  const boost::string_ref &field_name = field_it->get()->getNameRef();
1073 
1074  // Ents in the field meshset
1075  Range ents_of_id_meshset;
1076  CHKERR get_moab().get_entities_by_handle((*field_it)->meshSet,
1077  ents_of_id_meshset, false);
1078  if (verb > VERY_NOISY) {
1079  PetscSynchronizedPrintf(PETSC_COMM_SELF, "Ents in field %s meshset %d\n",
1080  (*field_it)->getName().c_str(),
1081  ents_of_id_meshset.size());
1082  }
1083 
1084  for (auto p_eit = ents_of_id_meshset.pair_begin();
1085  p_eit != ents_of_id_meshset.pair_end(); ++p_eit) {
1086 
1087  const EntityHandle first = p_eit->first;
1088  const EntityHandle second = p_eit->second;
1089  const auto lo_uid =
1090  FieldEntity::getLoLocalEntityBitNumber(bit_number, first);
1091  const auto hi_uid =
1092  FieldEntity::getHiLocalEntityBitNumber(bit_number, second);
1093 
1094  auto feit = entsFields.get<Unique_mi_tag>().lower_bound(lo_uid);
1095  if (feit == entsFields.get<Unique_mi_tag>().end())
1096  continue;
1097  auto hi_feit = entsFields.get<Unique_mi_tag>().upper_bound(hi_uid);
1098 
1099  // If there are DOFs in that range is more pragmatic to remove them
1100  // rather than to find sub-ranges or make them inactive
1101  auto dit = dofsField.get<Unique_mi_tag>().lower_bound(lo_uid);
1102  auto hi_dit = dofsField.get<Unique_mi_tag>().upper_bound(hi_uid);
1103  dofsField.get<Unique_mi_tag>().erase(dit, hi_dit);
1104 
1105  // Add vertices DOFs by bulk
1106  boost::shared_ptr<std::vector<DofEntity>> dofs_array =
1107  boost::make_shared<std::vector<DofEntity>>(std::vector<DofEntity>());
1108  // Add Sequence of DOFs to sequence container as weak_ptr
1109  int nb_dofs_on_ents = 0;
1110  for (auto tmp_feit = feit; tmp_feit != hi_feit; ++tmp_feit) {
1111  nb_dofs_on_ents += rank * tmp_feit->get()->getOrderNbDofs(
1112  tmp_feit->get()->getMaxOrder());
1113  }
1114 
1115  // Reserve memory
1116  dofs_array->reserve(nb_dofs_on_ents);
1117 
1118  // Create DOFs
1119  for (; feit != hi_feit; ++feit) {
1120  // Create dofs instances and shared pointers
1121  int DD = 0;
1122  // Loop orders (loop until max entity order is set)
1123  for (int oo = 0; oo <= feit->get()->getMaxOrder(); ++oo) {
1124  // Loop nb. dofs at order oo
1125  for (int dd = 0; dd < feit->get()->getOrderNbDofsDiff(oo); ++dd) {
1126  // Loop rank
1127  for (int rr = 0; rr < rank; ++rr, ++DD) {
1128  dofs_array->emplace_back(*feit, oo, rr, DD);
1129  ++dof_counter[feit->get()->getEntType()];
1130  }
1131  }
1132  }
1133  if (DD > feit->get()->getNbDofsOnEnt()) {
1134  std::ostringstream ss;
1135  ss << "rank " << rAnk << " ";
1136  ss << **feit << std::endl;
1137  SETERRQ3(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1138  "Expected number of DOFs on entity not equal to number added "
1139  "to database (DD = %d != %d = "
1140  "feit->get()->getNbDofsOnEnt())\n"
1141  "%s",
1142  DD, feit->get()->getNbDofsOnEnt(), ss.str().c_str());
1143  }
1144  }
1145 
1146  // Insert into Multi-Index container
1147  int dofs_field_size0 = dofsField.size();
1148  auto hint = dofsField.end();
1149  for (auto &v : *dofs_array)
1150  hint = dofsField.emplace_hint(hint, dofs_array, &v);
1151 
1152  // Add Sequence of DOFs to sequence container as weak_ptr
1153  field_it->get()->getDofSequenceContainer().push_back(dofs_array);
1154 
1155  // Check data consistency
1156  if (PetscUnlikely(static_cast<int>(dofs_array.use_count()) !=
1157  static_cast<int>(dofs_array->size() + 1))) {
1158  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1159  "Wrong use count %d != %d", dofs_array.use_count(),
1160  dofs_array->size() + 1);
1161  }
1162  if (dofs_field_size0 + dofs_array->size() != dofsField.size()) {
1163  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1164  "Wrong number of inserted DOFs %d != %d", dofs_array->size(),
1165  dofsField.size() - dofs_field_size0);
1166  }
1167  }
1169 }
1170 
1171 MoFEMErrorCode Core::buildField(const boost::shared_ptr<Field> &field,
1172  int verb) {
1174  if (verb == -1)
1175  verb = verbose;
1176  if (verb > QUIET)
1177  MOFEM_LOG("SYNC", Sev::verbose) << "Build field " << field->getName();
1178 
1179  std::map<EntityType, int> dof_counter;
1180  std::map<EntityType, int> inactive_dof_counter;
1181 
1182  // Need to rebuild order table since number of dofs on each order when
1183  // field was created.
1184  if (field->getApproxBase() == USER_BASE)
1185  CHKERR field->rebuildDofsOrderMap();
1186 
1187  switch (field->getSpace()) {
1188  case NOFIELD:
1189  CHKERR this->buildFieldForNoField(field->getId(), dof_counter, verb);
1190  break;
1191  case L2:
1192  case H1:
1193  case HCURL:
1194  case HDIV:
1195  CHKERR this->buildFieldForL2H1HcurlHdiv(field->getId(), dof_counter,
1196  inactive_dof_counter, verb);
1197  break;
1198  default:
1199  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
1200  }
1201 
1202  if (verb > QUIET) {
1203  int nb_added_dofs = 0;
1204  int nb_inactive_added_dofs = 0;
1205  for (auto const &it : dof_counter) {
1206  switch (it.first) {
1207  case MBVERTEX:
1208  MOFEM_LOG("SYNC", Sev::verbose)
1209  << "Nb. of dofs (vertices) " << it.second << " (inactive "
1210  << inactive_dof_counter[it.first] << ")";
1211  break;
1212  case MBEDGE:
1213  MOFEM_LOG("SYNC", Sev::verbose)
1214  << "Nb. of dofs (edge) " << it.second << " (inactive "
1215  << inactive_dof_counter[it.first] << ")";
1216  break;
1217  case MBTRI:
1218  MOFEM_LOG("SYNC", Sev::verbose)
1219  << "Nb. of dofs (triangles) " << it.second << " (inactive "
1220  << inactive_dof_counter[it.first] << ")";
1221  break;
1222  case MBQUAD:
1223  MOFEM_LOG("SYNC", Sev::verbose)
1224  << "Nb. of dofs (quads) " << it.second << " (inactive "
1225  << inactive_dof_counter[it.first] << ")";
1226  break;
1227  case MBTET:
1228  MOFEM_LOG("SYNC", Sev::verbose)
1229  << "Nb. of dofs (tetrahedra) " << it.second << " (inactive "
1230  << inactive_dof_counter[it.first] << ")";
1231  break;
1232  case MBPRISM:
1233  MOFEM_LOG("SYNC", Sev::verbose)
1234  << "Nb. of dofs (prisms) " << it.second << " (inactive "
1235  << inactive_dof_counter[it.first] << ")";
1236  break;
1237  case MBENTITYSET:
1238  MOFEM_LOG("SYNC", Sev::verbose)
1239  << "Nb. of dofs (meshsets) " << it.second << " (inactive "
1240  << inactive_dof_counter[it.first] << ")";
1241  break;
1242  default:
1243  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
1244  }
1245  nb_added_dofs += it.second;
1246  nb_inactive_added_dofs += inactive_dof_counter[it.first];
1247  }
1248  if (verb > QUIET) {
1249  MOFEM_LOG("SYNC", Sev::verbose)
1250  << "Nb. added dofs " << nb_added_dofs << " (number of inactive dofs "
1251  << nb_inactive_added_dofs << " )";
1252  }
1253  }
1255 }
1256 
1257 MoFEMErrorCode Core::build_field(const std::string field_name, int verb) {
1259  auto field_it = fIelds.get<FieldName_mi_tag>().find(field_name);
1260  if (field_it == fIelds.get<FieldName_mi_tag>().end())
1261  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "Field < %s > not found",
1262  field_name.c_str());
1263 
1264  CHKERR this->buildField(*field_it, verb);
1265  if (verb > QUIET)
1266  MOFEM_LOG_SYNCHRONISE(cOmm);
1268 }
1269 
1272  if (verb == -1)
1273  verb = verbose;
1274 
1275  for (auto field : fIelds.get<BitFieldId_mi_tag>())
1276  CHKERR this->buildField(field, verb);
1277 
1278  *buildMoFEM = 1 << 0;
1279  if (verb > QUIET) {
1280  MOFEM_LOG("SYNC", Sev::inform) << "Number of dofs " << dofsField.size();
1281  MOFEM_LOG_SYNCHRONISE(cOmm);
1282  }
1283 
1285 }
1286 
1288 Core::list_dofs_by_field_name(const std::string &field_name) const {
1290  auto dit = dofsField.get<Unique_mi_tag>().lower_bound(
1291  FieldEntity::getLoBitNumberUId(get_field_bit_number((field_name))));
1292  auto hi_dit = dofsField.get<Unique_mi_tag>().upper_bound(
1293  FieldEntity::getHiBitNumberUId(get_field_bit_number(field_name)));
1294  MOFEM_LOG("SYNC", Sev::inform) << "List DOFs:";
1295  for (; dit != hi_dit; dit++)
1296  MOFEM_LOG("SYNC", Sev::inform) << *dit;
1297 
1298  MOFEM_LOG_SYNCHRONISE(cOmm);
1300 }
1301 
1304  MOFEM_LOG("SYNC", Sev::inform) << "List Fields:";
1305  for (auto &miit : fIelds.get<BitFieldId_mi_tag>())
1306  MOFEM_LOG("SYNC", Sev::inform) << *miit;
1307 
1308  MOFEM_LOG_SYNCHRONISE(cOmm);
1310 }
1311 
1313 Core::get_problem_finite_elements_entities(const std::string &problem_name,
1314  const std::string &fe_name,
1315  const EntityHandle meshset) {
1317  auto &prb = pRoblems.get<Problem_mi_tag>();
1318  auto p_miit = prb.find(problem_name);
1319  if (p_miit == prb.end())
1320  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1321  "No such problem like < %s >", problem_name.c_str());
1322  auto miit = p_miit->numeredFiniteElementsPtr->get<FiniteElement_name_mi_tag>()
1323  .lower_bound(fe_name);
1324  auto hi_miit = p_miit->numeredFiniteElementsPtr->get<FiniteElement_name_mi_tag>()
1325  .upper_bound(fe_name);
1326  for (; miit != hi_miit; miit++) {
1327  EntityHandle ent = (*miit)->getEnt();
1328  CHKERR get_moab().add_entities(meshset, &ent, 1);
1329  const int part = (*miit)->getPart();
1330  CHKERR get_moab().tag_set_data(th_Part, &ent, 1, &part);
1331  }
1333 }
1334 
1335 FieldEntityByUId::iterator
1336 Core::get_ent_field_by_name_begin(const std::string &field_name) const {
1337  return entsFields.get<Unique_mi_tag>().lower_bound(
1338  FieldEntity::getLoBitNumberUId(get_field_bit_number(field_name)));
1339 }
1340 FieldEntityByUId::iterator
1341 Core::get_ent_field_by_name_end(const std::string &field_name) const {
1342  return entsFields.get<Unique_mi_tag>().upper_bound(
1343  FieldEntity::getHiBitNumberUId(get_field_bit_number(field_name)));
1344 }
1345 DofEntityByUId::iterator
1346 Core::get_dofs_by_name_begin(const std::string &field_name) const {
1347  return dofsField.get<Unique_mi_tag>().lower_bound(
1348  FieldEntity::getLoBitNumberUId(get_field_bit_number(field_name)));
1349 }
1350 DofEntityByUId::iterator
1351 Core::get_dofs_by_name_end(const std::string &field_name) const {
1352  return dofsField.get<Unique_mi_tag>().upper_bound(
1353  FieldEntity::getHiBitNumberUId(get_field_bit_number(field_name)));
1354 }
1355 DofEntityByUId::iterator
1356 Core::get_dofs_by_name_and_ent_begin(const std::string &field_name,
1357  const EntityHandle ent) const {
1358  return dofsField.get<Unique_mi_tag>().lower_bound(
1359  FieldEntity::getLoLocalEntityBitNumber(get_field_bit_number(field_name),
1360  ent));
1361 }
1362 DofEntityByUId::iterator
1363 Core::get_dofs_by_name_and_ent_end(const std::string &field_name,
1364  const EntityHandle ent) const {
1365  return dofsField.get<Unique_mi_tag>().upper_bound(
1366  FieldEntity::getHiLocalEntityBitNumber(get_field_bit_number(field_name),
1367  ent));
1368 }
1369 DofEntityByUId::iterator
1370 Core::get_dofs_by_name_and_type_begin(const std::string &field_name,
1371  const EntityType type) const {
1372  return dofsField.get<Unique_mi_tag>().lower_bound(
1373  FieldEntity::getLoLocalEntityBitNumber(get_field_bit_number(field_name),
1375 }
1376 DofEntityByUId::iterator
1377 Core::get_dofs_by_name_and_type_end(const std::string &field_name,
1378  const EntityType type) const {
1379  return dofsField.get<Unique_mi_tag>().upper_bound(
1380  FieldEntity::getHiLocalEntityBitNumber(get_field_bit_number(field_name),
1382 }
1384 Core::check_number_of_ents_in_ents_field(const std::string &name) const {
1386  auto it = fIelds.get<FieldName_mi_tag>().find(name);
1387  if (it == fIelds.get<FieldName_mi_tag>().end()) {
1388  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1389  "field not found < %s >", name.c_str());
1390  }
1391  EntityHandle meshset = (*it)->getMeshset();
1392  int num_entities;
1393  CHKERR get_moab().get_number_entities_by_handle(meshset, num_entities);
1394 
1395  auto count_field_ents = [&]() {
1396  auto bit_number = (*it)->getBitNumber();
1397  auto low_eit = entsFields.get<Unique_mi_tag>().lower_bound(
1398  FieldEntity::getLoBitNumberUId(bit_number));
1399  auto hi_eit = entsFields.get<Unique_mi_tag>().upper_bound(
1400  FieldEntity::getHiBitNumberUId(bit_number));
1401  return std::distance(low_eit, hi_eit);
1402  };
1403 
1404  if (count_field_ents() > (unsigned int)num_entities) {
1405  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1406  "not equal number of entities in meshset and field multiindex "
1407  "< %s >",
1408  name.c_str());
1409  }
1411 }
1414  for (auto &it : fIelds.get<FieldName_mi_tag>()) {
1415  if (it->getSpace() == NOFIELD)
1416  continue; // FIXME: should be treated properly, not test is just
1417  // skipped for this NOFIELD space
1418  EntityHandle meshset = it->getMeshset();
1419  int num_entities;
1420  CHKERR get_moab().get_number_entities_by_handle(meshset, num_entities);
1421 
1422  auto count_field_ents = [&]() {
1423  auto bit_number = it->getBitNumber();
1424  auto low_eit = entsFields.get<Unique_mi_tag>().lower_bound(
1425  FieldEntity::getLoBitNumberUId(bit_number));
1426  auto hi_eit = entsFields.get<Unique_mi_tag>().upper_bound(
1427  FieldEntity::getHiBitNumberUId(bit_number));
1428  return std::distance(low_eit, hi_eit);
1429  };
1430 
1431  if (count_field_ents() > (unsigned int)num_entities) {
1432  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1433  "not equal number of entities in meshset and field "
1434  "multiindex < %s >",
1435  it->getName().c_str());
1436  }
1437  }
1439 }
1440 
1441 } // namespace MoFEM
FieldEntityByUId::iterator get_ent_field_by_name_end(const std::string &field_name) const
get begin iterator of filed dofs of given name (instead you can use IT_GET_ENT_FIELD_BY_NAME_FOR_LOOP...
Definition: FieldCore.cpp:1341
EntityHandle get_id_for_max_type()
MoFEMErrorCode addField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_coefficients, const TagType tag_type, const enum MoFEMTypes bh, int verb)
Template for add_field.
Definition: FieldCore.cpp:109
MoFEMErrorCode list_fields() const
list entities in the field
Definition: FieldCore.cpp:1302
EntityHandle get_id_for_min_type()
multi_index_container< boost::shared_ptr< FieldEntity >, indexed_by< sequenced<>, ordered_non_unique< tag< Ent_mi_tag >, const_mem_fun< FieldEntity::interface_type_RefEntity, EntityHandle, &FieldEntity::getEnt > > > > FieldEntity_multiIndex_ent_view
DofEntityByUId::iterator get_dofs_by_name_and_ent_end(const std::string &field_name, const EntityHandle ent) const
get begin iterator of filed dofs of given name and ent (instead you can use IT_GET_DOFS_FIELD_BY_NAME...
Definition: FieldCore.cpp:1363
MoFEMErrorCode buildFieldForNoFieldImpl(boost::shared_ptr< Field > field_ptr, std::map< EntityType, int > &dof_counter, int verb)
Definition: FieldCore.cpp:929
field with continuous normal traction
Definition: definitions.h:179
static UId getHiLocalEntityBitNumber(const char bit_number, const EntityHandle ent)
EntityHandle get_field_meshset(const BitFieldId id) const
Definition: FieldCore.cpp:49
user implemented approximation base
Definition: definitions.h:160
FieldBitNumber get_field_bit_number(const std::string name) const
get field bit number
Definition: FieldCore.cpp:41
MoFEMErrorCode set_field_order_by_entity_type_and_bit_ref(const BitRefLevel &bit, const BitRefLevel &mask, const EntityType type, const BitFieldId id, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)
Definition: FieldCore.cpp:900
Provide data structure for (tensor) field approximation.The Field is intended to provide support for ...
BitFieldId getBitFieldId(const std::string &name) const
Definition: FieldCore.cpp:31
MoFEMErrorCode getEntitiesByTypeAndRefLevel(const BitRefLevel bit, const BitRefLevel mask, const EntityType type, const EntityHandle meshset, int verb=0) const
add all ents from ref level given by bit to meshset
DofEntityByUId::iterator get_dofs_by_name_and_type_begin(const std::string &field_name, const EntityType type) const
get begin iterator of filed dofs of given name and ent type (instead you can use IT_GET_DOFS_FIELD_BY...
Definition: FieldCore.cpp:1370
int ApproximationOrder
Approximation on the entity.
Definition: Types.hpp:37
MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)
Add entities to field meshset.
Definition: FieldCore.cpp:352
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
scalar or vector of scalars describe (no true field)
Definition: definitions.h:176
#define MOFEM_LOG_TAG(channel, tag)
Tag channelTag channel tag is set until MOFEM_LOG_CHANNEL is called, then new tag can be set.
Definition: LogManager.hpp:334
MoFEMErrorCode set_field_order(const Range &ents, const BitFieldId id, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)
Definition: FieldCore.cpp:839
#define FieldCoreFunctionBegin
Definition: FieldCore.cpp:21
static Index< 'p', 3 > p
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:303
static UId getLoLocalEntityBitNumber(const char bit_number, const EntityHandle ent)
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:306
#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
MoFEMErrorCode get_problem_finite_elements_entities(const std::string &name, const std::string &fe_name, const EntityHandle meshset)
add finite elements to the meshset
Definition: FieldCore.cpp:1313
MoFEMErrorCode buildFieldForL2H1HcurlHdiv(const BitFieldId id, std::map< EntityType, int > &dof_counter, std::map< EntityType, int > &inactive_dof_counter, int verb=DEFAULT_VERBOSITY)
Definition: FieldCore.cpp:1057
MoFEMErrorCode get_field_entities_by_handle(const std::string name, Range &ents) const
get entities in the field by handle
Definition: FieldCore.cpp:101
#define THROW_MESSAGE(a)
Throw MoFEM exception.
Definition: definitions.h:628
MoFEMErrorCode check_number_of_ents_in_ents_field() const
check data consistency in entitiesPtr
Definition: FieldCore.cpp:1412
MoFEMErrorCode setFieldOrderImpl(boost::shared_ptr< Field > field_ptr, const Range &ents, const ApproximationOrder order, int verb)
Definition: FieldCore.cpp:429
MoFEMErrorCode build_field(const std::string field_name, int verb=DEFAULT_VERBOSITY)
build field by name
Definition: FieldCore.cpp:1257
bool check_field(const std::string &name) const
check if field is in database
Definition: FieldCore.cpp:61
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
char FieldBitNumber
Field bit number.
Definition: Types.hpp:39
MoFEMErrorCode add_ents_to_field_by_dim(const Range &ents, const int dim, const std::string &name, int verb=DEFAULT_VERBOSITY)
Add entities to field meshset.
Definition: FieldCore.cpp:345
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
DofEntityByUId::iterator get_dofs_by_name_and_type_end(const std::string &field_name, const EntityType ent) const
get begin iterator of filed dofs of given name end ent type(instead you can use IT_GET_DOFS_FIELD_BY_...
Definition: FieldCore.cpp:1377
MoFEMErrorCode create_vertices_and_add_to_field(const std::string name, const double coords[], int size, int verb=DEFAULT_VERBOSITY)
Create a vertices and add to field object.
Definition: FieldCore.cpp:390
static Index< 'n', 3 > n
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:85
const int dim
MoFEMErrorCode buildField(const boost::shared_ptr< Field > &field, int verb=DEFAULT_VERBOSITY)
Definition: FieldCore.cpp:1171
FieldEntityByUId::iterator get_ent_field_by_name_begin(const std::string &field_name) const
get begin iterator of filed ents of given name (instead you can use IT_GET_ENT_FIELD_BY_NAME_FOR_LOOP...
Definition: FieldCore.cpp:1336
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
constexpr int order
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
static Index< 'i', 3 > i
FieldApproximationBase
approximation base
Definition: definitions.h:150
MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)
Definition: FieldCore.cpp:1270
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
MoFEMErrorCode get_field_entities_by_type(const std::string name, EntityType type, Range &ents) const
get entities in the field by type
Definition: FieldCore.cpp:92
Managing BitRefLevels.
DofEntityByUId::iterator get_dofs_by_name_begin(const std::string &field_name) const
get begin iterator of filed dofs of given name (instead you can use IT_GET_DOFS_FIELD_BY_NAME_FOR_LOO...
Definition: FieldCore.cpp:1346
field with continuous tangents
Definition: definitions.h:178
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:279
MoFEMErrorCode addEntsToFieldByDim(const Range &ents, const int dim, const std::string &name, int verb=DEFAULT_VERBOSITY)
Definition: FieldCore.cpp:266
FieldSpace
approximation spaces
Definition: definitions.h:174
MoFEMErrorCode get_field_entities_by_dimension(const std::string name, int dim, Range &ents) const
get entities in the field by dimension
Definition: FieldCore.cpp:82
std::bitset< BITFIELDID_SIZE > BitFieldId
Field Id.
Definition: Types.hpp:53
#define CHKERR
Inline error check.
Definition: definitions.h:604
static boost::shared_ptr< FieldData *const > makeSharedFieldDataAdaptorPtr(const boost::shared_ptr< Field > &field_ptr, const boost::shared_ptr< RefEntity > &ref_ents_ptr)
Return shared pointer to entity field data vector adaptor.
const int getValue() const
Definition: Core.hpp:33
DofEntityByUId::iterator get_dofs_by_name_and_ent_begin(const std::string &field_name, const EntityHandle ent) const
get begin iterator of filed dofs of given name and ent(instead you can use IT_GET_DOFS_FIELD_BY_NAME_...
Definition: FieldCore.cpp:1356
const Field * get_field_structure(const std::string &name)
get field structure
Definition: FieldCore.cpp:69
structure to change FieldEntity order
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
MultiIndex Tag for field name.
MoFEMErrorCode list_dofs_by_field_name(const std::string &name) const
Definition: FieldCore.cpp:1288
static UId getLoBitNumberUId(const FieldBitNumber bit_number)
#define MOFEM_LOG_ATTRIBUTES(channel, bit)
Add attributes to channel.
Definition: LogManager.hpp:291
int FieldCoefficientsNumber
Number of field coefficients.
Definition: Types.hpp:38
DofEntityByUId::iterator get_dofs_by_name_end(const std::string &field_name) const
get begin iterator of filed dofs of given name (instead you can use IT_GET_DOFS_FIELD_BY_NAME_FOR_LOO...
Definition: FieldCore.cpp:1351
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
continuous field
Definition: definitions.h:177
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74
#define MOFEM_LOG_FUNCTION()
Set scopeMacro for function scope markup. The scope name is constructed with help of compiler and con...
Definition: LogManager.hpp:320
MoFEMErrorCode buildFieldForNoField(const BitFieldId id, std::map< EntityType, int > &dof_counter, int verb=DEFAULT_VERBOSITY)
Definition: FieldCore.cpp:1034
MoFEMTypes
Those types control how functions respond on arguments, f.e. error handling.
Definition: definitions.h:189
static UId getHiBitNumberUId(const FieldBitNumber bit_number)
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:340
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)
Add filed.
Definition: FieldCore.cpp:257
field with C-1 continuity
Definition: definitions.h:180