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