v0.9.1
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 namespace MoFEM {
20 
21 BitFieldId Core::getBitFieldId(const std::string &name) const {
22  auto &set = fIelds.get<FieldName_mi_tag>();
23  auto miit = set.find(name);
24  if (miit == set.end()) {
25  THROW_MESSAGE("field < " + name +
26  " > not in database (top tip: check spelling)");
27  }
28  return (*miit)->getId();
29 }
30 
32  auto &set = fIelds.get<BitFieldId_mi_tag>();
33  auto miit = set.find(id);
34  if (miit == set.end())
35  THROW_MESSAGE("field not in database (top tip: check spelling)");
36  return (*miit)->meshSet;
37 }
38 
39 EntityHandle Core::get_field_meshset(const std::string &name) const {
40  return get_field_meshset(getBitFieldId(name));
41 }
42 
43 bool Core::check_field(const std::string &name) const {
44  auto &set = fIelds.get<FieldName_mi_tag>();
45  auto miit = set.find(name);
46  if (miit == set.end())
47  return false;
48  return true;
49 }
50 
51 const Field *Core::get_field_structure(const std::string &name) {
52  auto &set = fIelds.get<FieldName_mi_tag>();
53  auto miit = set.find(name);
54  if (miit == set.end()) {
55  throw MoFEMException(
57  std::string("field < " + name +
58  " > not in database (top tip: check spelling)")
59  .c_str());
60  }
61  return miit->get();
62 }
63 
65  int dim,
66  Range &ents) const {
67 
69  EntityHandle meshset = get_field_meshset(name);
70  CHKERR get_moab().get_entities_by_dimension(meshset, dim, ents, true);
72 }
73 
75  EntityType type,
76  Range &ents) const {
78  EntityHandle meshset = get_field_meshset(name);
79  CHKERR get_moab().get_entities_by_type(meshset, type, ents, true);
81 }
82 
84  Range &ents) const {
86  EntityHandle meshset = get_field_meshset(name);
87  CHKERR get_moab().get_entities_by_handle(meshset, ents, true);
89 }
90 
91 MoFEMErrorCode Core::add_field(const std::string &name, const FieldSpace space,
92  const FieldApproximationBase base,
93  const FieldCoefficientsNumber nb_of_coefficients,
94  const TagType tag_type, const enum MoFEMTypes bh,
95  int verb) {
97  if (verb == -1)
98  verb = verbose;
99  *buildMoFEM = 0;
100  auto fit = fIelds.get<FieldName_mi_tag>().find(name);
101  if (fit != fIelds.get<FieldName_mi_tag>().end()) {
102 
103  if (bh == MF_EXCL)
104  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
105  "field is <%s> in database", name.c_str());
106 
107  } else {
108 
109  EntityHandle meshset;
110  CHKERR get_moab().create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER,
111  meshset);
112 
113  // Add field mesh set to partion meshset. In case of no elements
114  // on processor part, when mesh file is read, finite element meshset is
115  // prevented from deletion by moab reader.
116  auto add_meshset_to_partition = [&](auto meshset) {
118  const void *tag_vals[] = {&rAnk};
119  ParallelComm *pcomm = ParallelComm::get_pcomm(
120  &get_moab(), get_basic_entity_data_ptr()->pcommID);
121  Tag part_tag = pcomm->part_tag();
122  Range tagged_sets;
123  CHKERR get_moab().get_entities_by_type_and_tag(0, MBENTITYSET, &part_tag,
124  tag_vals, 1, tagged_sets,
125  moab::Interface::UNION);
126  for (auto s : tagged_sets)
127  CHKERR get_moab().add_entities(s, &meshset, 1);
129  };
130  CHKERR add_meshset_to_partition(meshset);
131 
132  // id
133  BitFieldId id = getFieldShift();
134 
135  auto create_tags = [&]() {
137  CHKERR
138  get_moab().tag_set_data(th_FieldId, &meshset, 1, &id);
139  // space
140  CHKERR get_moab().tag_set_data(th_FieldSpace, &meshset, 1, &space);
141  // base
142  CHKERR get_moab().tag_set_data(th_FieldBase, &meshset, 1, &base);
143 
144  // name
145  void const *tag_data[] = {name.c_str()};
146  int tag_sizes[1];
147  tag_sizes[0] = name.size();
148  CHKERR get_moab().tag_set_by_ptr(th_FieldName, &meshset, 1, tag_data,
149  tag_sizes);
150  // name data prefix
151  const std::string name_data_prefix("_App_Data");
152  void const *tag_prefix_data[] = {name_data_prefix.c_str()};
153  int tag_prefix_sizes[1];
154  tag_prefix_sizes[0] = name_data_prefix.size();
155  CHKERR get_moab().tag_set_by_ptr(th_FieldName_DataNamePrefix, &meshset, 1,
156  tag_prefix_data, tag_prefix_sizes);
157  Tag th_app_order, th_field_data, th_field_data_vert, th_rank;
158  // data
159  std::string tag_data_name = name_data_prefix + name;
160  const int def_len = 0;
161  CHKERR get_moab().tag_get_handle(
162  tag_data_name.c_str(), def_len, MB_TYPE_DOUBLE, th_field_data,
163  MB_TAG_CREAT | MB_TAG_VARLEN | MB_TAG_SPARSE, NULL);
164  std::string tag_data_name_verts = name_data_prefix + name + "V";
165  VectorDouble def_vert_data(nb_of_coefficients);
166  def_vert_data.clear();
167  CHKERR get_moab().tag_get_handle(
168  tag_data_name_verts.c_str(), nb_of_coefficients, MB_TYPE_DOUBLE,
169  th_field_data_vert, MB_TAG_CREAT | tag_type, &*def_vert_data.begin());
170  // order
171  ApproximationOrder def_ApproximationOrder = -1;
172  const std::string Tag_ApproximationOrder_name = "_App_Order_" + name;
173  CHKERR get_moab().tag_get_handle(
174  Tag_ApproximationOrder_name.c_str(), 1, MB_TYPE_INTEGER, th_app_order,
175  MB_TAG_CREAT | tag_type, &def_ApproximationOrder);
176  // rank
177  int def_rank = 1;
178  const std::string tag_rank_name = "_Field_Rank_" + name;
179  CHKERR get_moab().tag_get_handle(tag_rank_name.c_str(), 1,
180  MB_TYPE_INTEGER, th_rank,
181  MB_TAG_CREAT | MB_TAG_SPARSE, &def_rank);
182  CHKERR get_moab().tag_set_data(th_rank, &meshset, 1, &nb_of_coefficients);
183 
185  };
186 
187  auto create_undefined_cs = [&](auto &undefined_cs_ptr) {
189  const int unidefined_cs_dim[] = {-1, 0, 0, 0};
190  CHKERR getInterface<CoordSystemsManager>()->addCoordinateSystem(
191  unidefined_cs_dim, "UNDEFINED", MF_ZERO);
192  CHKERR getInterface<CoordSystemsManager>()->getCoordSysPtr(
193  "UNDEFINED", undefined_cs_ptr);
195  };
196 
197  auto add_field_meshset_to_cs = [&](auto &undefined_cs_ptr) {
199  int cs_name_size[1];
200  cs_name_size[0] = undefined_cs_ptr->getName().size();
201  void const *cs_name[] = {&*undefined_cs_ptr->getNameRef().begin()};
202  CHKERR get_moab().tag_set_by_ptr(
203  getInterface<CoordSystemsManager>()->get_th_CoordSysName(), &meshset,
204  1, cs_name, cs_name_size);
205  EntityHandle coord_sys_meshset = undefined_cs_ptr->getMeshset();
206  CHKERR get_moab().add_entities(coord_sys_meshset, &meshset, 1);
208  };
209 
210  CHKERR create_tags();
211  boost::shared_ptr<CoordSys> undefined_cs_ptr;
212  CHKERR create_undefined_cs(undefined_cs_ptr);
213  CHKERR add_field_meshset_to_cs(undefined_cs_ptr);
214 
215  auto p = fIelds.insert(
216  boost::make_shared<Field>(moab, meshset, undefined_cs_ptr));
217  if (bh == MF_EXCL) {
218  if (!p.second)
219  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND,
220  "field not inserted %s (top tip, it could be already there)",
221  Field(moab, meshset, undefined_cs_ptr).getName().c_str());
222  }
223 
224  if (verb > QUIET) {
225  std::ostringstream ss;
226  ss << "add: " << **p.first << std::endl;
227  PetscPrintf(cOmm, ss.str().c_str());
228  }
229  }
230 
232 }
233 
234 MoFEMErrorCode Core::addEntsToFieldByDim(const Range &ents, const int dim,
235  const std::string &name, int verb) {
236 
237  *buildMoFEM = 0;
238  EntityHandle idm = no_handle;
239  if (verb == -1)
240  verb = verbose;
242  idm = get_field_meshset(name);
243  FieldSpace space;
244  CHKERR get_moab().tag_get_data(th_FieldSpace, &idm, 1, &space);
245  std::vector<int> nb_ents_on_dim(3, 0);
246  switch (space) {
247  case L2:
248  CHKERR get_moab().add_entities(idm, ents);
249  if (verb >= VERY_VERBOSE) {
250  std::ostringstream ss;
251  ss << "add entities to field " << name;
252  ss << " nb. add ents " << ents.size();
253  ss << std::endl;
254  PetscSynchronizedPrintf(cOmm, ss.str().c_str());
255  }
256  break;
257  case H1:
258  CHKERR get_moab().add_entities(idm, ents);
259  for (int dd = 0; dd != dim; ++dd) {
260  Range adj_ents;
261  CHKERR get_moab().get_adjacencies(ents, dd, false, adj_ents,
262  moab::Interface::UNION);
263  if (dd == 0) {
264  Range topo_nodes;
265  CHKERR get_moab().get_connectivity(ents, topo_nodes, true);
266  Range mid_nodes;
267  CHKERR get_moab().get_connectivity(ents, mid_nodes, false);
268  mid_nodes = subtract(mid_nodes, topo_nodes);
269  adj_ents = subtract(adj_ents, mid_nodes);
270  }
271  CHKERR get_moab().add_entities(idm, adj_ents);
272  nb_ents_on_dim[dd] = adj_ents.size();
273  }
274  break;
275  case HCURL:
276  CHKERR get_moab().add_entities(idm, ents);
277  for (int dd = 1; dd != dim; ++dd) {
278  Range adj_ents;
279  CHKERR get_moab().get_adjacencies(ents, dd, false, adj_ents,
280  moab::Interface::UNION);
281  CHKERR get_moab().add_entities(idm, adj_ents);
282  nb_ents_on_dim[dd] = adj_ents.size();
283  }
284  break;
285  case HDIV:
286  CHKERR get_moab().add_entities(idm, ents);
287  if (dim > 2) {
288  Range adj_ents;
289  CHKERR get_moab().get_adjacencies(ents, 2, false, adj_ents,
290  moab::Interface::UNION);
291  CHKERR get_moab().add_entities(idm, adj_ents);
292  nb_ents_on_dim[2] = adj_ents.size();
293  }
294  break;
295  default:
296  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
297  "sorry, unknown space added to entity");
298  }
299  if (verb >= VERY_VERBOSE) {
300  std::ostringstream ss;
301  ss << "add entities to field " << name;
302  ss << " nb. add ents " << ents.size();
303  ss << " nb. add faces " << nb_ents_on_dim[2];
304  ss << " nb. add edges " << nb_ents_on_dim[1];
305  ss << " nb. add nodes " << nb_ents_on_dim[0];
306  ss << std::endl;
307  PetscSynchronizedPrintf(cOmm, ss.str().c_str());
308  PetscSynchronizedFlush(cOmm, PETSC_STDOUT);
309  }
311 }
312 
313 MoFEMErrorCode Core::add_ents_to_field_by_dim(const Range &ents, const int dim,
314  const std::string &name,
315  int verb) {
316  Range ents_dim = ents.subset_by_dimension(dim);
317  return addEntsToFieldByDim(ents_dim, dim, name, verb);
318 }
319 
321  const EntityType type,
322  const std::string &name,
323  int verb) {
325  Range ents_type = ents.subset_by_type(type);
326  if (!ents_type.empty()) {
327  const int dim = get_moab().dimension_from_handle(ents_type[0]);
328  CHKERR addEntsToFieldByDim(ents_type, dim, name, verb);
329  }
331 }
332 
334  const int dim,
335  const std::string &name,
336  const bool recursive, int verb) {
338  Range ents;
339  CHKERR get_moab().get_entities_by_dimension(meshset, dim, ents, recursive);
340  CHKERR addEntsToFieldByDim(ents, dim, name, verb);
342 }
343 
345  const EntityType type,
346  const std::string &name,
347  const bool recursive, int verb) {
349  Range ents;
350  CHKERR get_moab().get_entities_by_type(meshset, type, ents, recursive);
351  if (!ents.empty()) {
352  const int dim = get_moab().dimension_from_handle(ents[0]);
353  CHKERR addEntsToFieldByDim(ents, dim, name, verb);
354  }
356 }
357 
359  const double coords[],
360  int size, int verb) {
362 
363  if (verb == DEFAULT_VERBOSITY)
364  verb = verbose;
365 
366  Range verts;
367 
368  auto create_vertices = [&]() {
370 
371  vector<double *> arrays_coord;
372  EntityHandle startv = 0;
373  ReadUtilIface *iface;
374  CHKERR get_moab().query_interface(iface);
375  CHKERR iface->get_node_coords(3, size, 0, startv, arrays_coord);
376  verts.insert(startv, startv + size - 1);
377  for (int n = 0; n != size; ++n)
378  for (auto d : {0, 1, 2})
379  arrays_coord[d][n] = coords[3 * n + d];
380 
382  };
383 
384  auto add_verts_to_field = [&]() {
386  EntityHandle field_meshset = get_field_meshset(name);
387  CHKERR get_moab().add_entities(field_meshset, verts);
389  };
390 
391  CHKERR create_vertices();
392  CHKERR add_verts_to_field();
393 
395  }
396 
397  MoFEMErrorCode Core::set_field_order(const Range &ents, const BitFieldId id,
399  int verb) {
401  if (verb == DEFAULT_VERBOSITY)
402  verb = verbose;
403  *buildMoFEM = 0;
404 
405  // check field & meshset
406  auto miit = fIelds.get<BitFieldId_mi_tag>().find(id);
407  if (miit == fIelds.get<BitFieldId_mi_tag>().end())
408  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "no filed found");
409 
411  // intersection with field meshset
412  Range ents_of_id_meshset;
413  CHKERR get_moab().get_entities_by_handle(idm, ents_of_id_meshset, false);
414  Range field_ents = intersect(ents, ents_of_id_meshset);
415  if (verb > VERBOSE)
416  PetscSynchronizedPrintf(
417  cOmm, "nb. of ents for order change in the field <%s> %d\n",
418  miit->get()->getName().c_str(), field_ents.size());
419 
420  // ent view by field id (in set all MoabEnts has the same FieldId)
421  auto eiit = entsFields.get<FieldName_mi_tag>().lower_bound(
422  miit->get()->getNameRef());
423  FieldEntity_multiIndex_ent_view ents_id_view;
424  if (eiit != entsFields.get<FieldName_mi_tag>().end()) {
425  auto hi_eiit = entsFields.get<FieldName_mi_tag>().upper_bound(
426  miit->get()->getNameRef());
427  std::copy(eiit, hi_eiit, std::back_inserter(ents_id_view));
428  }
429  if (verb > VERBOSE)
430  PetscSynchronizedPrintf(
431  cOmm, "nb. of ents in the multi index field <%s> %d\n",
432  miit->get()->getName().c_str(), ents_id_view.size());
433 
434  // loop over ents
435  int nb_ents_set_order_up = 0;
436  int nb_ents_set_order_down = 0;
437  int nb_ents_set_order_new = 0;
438 
439  FieldEntity_change_order modify_order_no_size_change(order, false);
440  FieldEntity_change_order modify_order_size_change(order, true);
441 
442  for (auto pit = field_ents.const_pair_begin();
443  pit != field_ents.const_pair_end(); pit++) {
444  EntityHandle first = pit->first;
445  EntityHandle second = pit->second;
446 
447  // Sanity check
448  switch ((*miit)->getSpace()) {
449  case H1:
450  break;
451  case HCURL:
452  if (get_moab().type_from_handle(first) == MBVERTEX)
453  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
454  "Hcurl space on vertices makes no sense");
455 
456  break;
457  case HDIV:
458  if (get_moab().type_from_handle(first) == MBVERTEX)
459  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
460  "Hdiv space on vertices makes no sense");
461 
462  if (get_moab().type_from_handle(first) == MBEDGE)
463  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
464  "Hdiv space on edges makes no sense");
465 
466  break;
467  default:
468  break;
469  }
470 
471  // Entity is in database, change order only if needed
472  Range ents_in_database;
473  auto vit = ents_id_view.get<1>().lower_bound(first);
474  auto hi_vit = ents_id_view.get<1>().upper_bound(second);
475  if (order >= 0) {
476  for (; vit != hi_vit; ++vit) {
477  ents_in_database.insert(vit->get()->getEnt());
478  // entity is in database and order is changed or reset
479  const ApproximationOrder old_approximation_order =
480  (*vit)->getMaxOrder();
481 
482  if (old_approximation_order != order) {
483 
484  FieldEntity_multiIndex::iterator miit =
485  entsFields.get<Unique_mi_tag>().find(
486  (*vit)->getGlobalUniqueId());
487 
488  if ((*miit)->getMaxOrder() < order)
489  nb_ents_set_order_up++;
490  if ((*miit)->getMaxOrder() > order)
491  nb_ents_set_order_down++;
492 
493  // set dofs inactive if order is reduced, and set new order to entity
494  // if order is increased (note that dofs are not build if order is
495  // increased)
496 
497  bool can_change_size = true;
498  auto dit =
499  dofsField.get<Composite_Name_And_Ent_mi_tag>().lower_bound(
500  boost::make_tuple((*miit)->getNameRef(), (*miit)->getEnt()));
501  if (dit != dofsField.get<Composite_Name_And_Ent_mi_tag>().end()) {
502  auto hi_dit =
503  dofsField.get<Composite_Name_And_Ent_mi_tag>().upper_bound(
504  boost::make_tuple((*miit)->getNameRef(),
505  (*miit)->getEnt()));
506  if (dit != hi_dit)
507  can_change_size = false;
508  for (; dit != hi_dit; dit++) {
509  if ((*dit)->getDofOrder() > order) {
510  bool success = dofsField.modify(
511  dofsField.project<0>(dit), DofEntity_active_change(false));
512  if (!success)
513  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
514  "modification unsuccessful");
515  }
516  }
517  }
518 
519  bool success =
520  entsFields.modify(entsFields.project<0>(miit),
521  can_change_size ? modify_order_size_change
522  : modify_order_no_size_change);
523  if (!success)
524  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
525  "modification unsuccessful");
526  }
527  }
528  }
529 
530  Range new_ents = subtract(Range(first, second), ents_in_database);
531  for (Range::const_pair_iterator pit = new_ents.const_pair_begin();
532  pit != new_ents.const_pair_end(); ++pit) {
533  EntityHandle first = pit->first;
534  EntityHandle second = pit->second;
535  const EntityType ent_type = get_moab().type_from_handle(first);
536  auto get_nb_dofs_on_order = [&](const int order) {
537  return order >= 0 ? ((*miit)->getFieldOrderTable()[ent_type])(order)
538  : 0;
539  };
540  const int field_rank = (*miit)->getNbOfCoeffs();
541  const int nb_dofs_on_order = get_nb_dofs_on_order(order);
542  const int nb_dofs = nb_dofs_on_order * field_rank;
543 
544  // reserve memory for field dofs
545  boost::shared_ptr<std::vector<FieldEntity>> ents_array(
546  new std::vector<FieldEntity>());
547 
548  // Add sequence to field data structure. Note that entities are allocated
549  // once into vector. This vector is passed into sequence as a weak_ptr.
550  // Vector is destroyed at the point last entity inside that vector is
551  // destroyed.
552  miit->get()->getEntSequenceContainer().push_back(ents_array);
553  ents_array->reserve(second - first + 1);
554 
555  // Entity is not in database and order is changed or reset
556  auto miit_ref_ent = refinedEntities.get<Ent_mi_tag>().lower_bound(first);
557 
558  auto create_tags_for_max_order = [&](const Range &ents) {
560  if (order >= 0) {
561  std::vector<ApproximationOrder> o_vec(ents.size(), order);
562  CHKERR get_moab().tag_set_data((*miit)->th_AppOrder, ents,
563  &*o_vec.begin());
564  }
566  };
567 
568  auto create_tags_for_data = [&](const Range &ents) {
570  if (order >= 0) {
571 
572  if (nb_dofs > 0) {
573  if (ent_type == MBVERTEX) {
574  std::vector<FieldData> d_vec(nb_dofs * ents.size(), 0);
575  CHKERR get_moab().tag_set_data((*miit)->th_FieldDataVerts, ents,
576  &*d_vec.begin());
577  } else {
578  std::vector<int> tag_size(ents.size(), nb_dofs);
579  std::vector<FieldData> d_vec(nb_dofs, 0);
580  std::vector<void const *> d_vec_ptr(ents.size(),
581  &*d_vec.begin());
582  CHKERR get_moab().tag_set_by_ptr((*miit)->th_FieldData, ents,
583  &*d_vec_ptr.begin(),
584  &*tag_size.begin());
585  }
586  }
587  }
589  };
590 
591  auto get_ents_in_ref_ent = [&](auto miit_ref_ent) {
592  auto hi = refinedEntities.get<Ent_mi_tag>().upper_bound(second);
593  Range in;
594  for (; miit_ref_ent != hi; ++miit_ref_ent)
595  in.insert(miit_ref_ent->get()->getRefEnt());
596  return in;
597  };
598 
599  auto get_ents_max_order = [&](const Range &ents) {
600  boost::shared_ptr<std::vector<const void *>> vec(
601  new std::vector<const void *>());
602  vec->resize(ents.size());
603  CHKERR get_moab().tag_get_by_ptr((*miit)->th_AppOrder, ents,
604  &*vec->begin());
605  return vec;
606  };
607 
608  auto get_ents_field_data_vector_adaptor =
609  [&](const Range &ents,
610  boost::shared_ptr<std::vector<const void *>> &ents_max_orders) {
611  // create shared pointer and reserve memory
612  boost::shared_ptr<std::vector<double *>> vec(
613  new std::vector<double *>());
614  vec->reserve(ents.size());
615 
616  auto get_nb_dofs = [&](const auto order) {
617  return get_nb_dofs_on_order(order) * field_rank;
618  };
619 
620  if (order >= 0 && get_nb_dofs(order) == 0) {
621  // set empty vector adaptor
622  for (int i = 0; i != ents.size(); ++i)
623  vec->emplace_back(nullptr);
624  } else {
625  moab::ErrorCode rval;
626  std::vector<int> tag_size(ents.size());
627  std::vector<const void *> d_vec_ptr(ents.size());
628 
629  // get tags data
630  if (ent_type == MBVERTEX)
631  rval = get_moab().tag_get_by_ptr((*miit)->th_FieldDataVerts,
632  ents, &*d_vec_ptr.begin(),
633  &*tag_size.begin());
634  else
635  rval = get_moab().tag_get_by_ptr((*miit)->th_FieldData, ents,
636  &*d_vec_ptr.begin(),
637  &*tag_size.begin());
638 
639  auto cast = [](auto p) {
640  return const_cast<FieldData *const>(
641  static_cast<const FieldData *>(p));
642  };
643 
644  // some of entities has tag not set or zero dofs on entity
645  if (rval == MB_SUCCESS) {
646  // all is ok, all entities has tag set
647  auto tit = d_vec_ptr.begin();
648  auto oit = ents_max_orders->begin();
649  for (auto sit = tag_size.begin(); sit != tag_size.end();
650  ++sit, ++tit, ++oit)
651  vec->emplace_back(cast(*tit));
652 
653  } else {
654  // set empty vector adaptor
655  for (int i = 0; i != ents.size(); ++i)
656  vec->emplace_back(nullptr);
657  // check order on all entities, and if for that order non zero
658  // dofs are expected get pointer to tag data and reset vector
659  // adaptor
660  auto oit = ents_max_orders->begin();
661  auto dit = vec->begin();
662  for (auto eit = ents.begin(); eit != ents.end();
663  ++eit, ++oit, ++dit) {
664  if (get_nb_dofs(
665  *static_cast<const ApproximationOrder *>(*oit))) {
666  int tag_size;
667  const void *ret_val;
668  if (ent_type == MBVERTEX)
669  CHKERR get_moab().tag_get_by_ptr(
670  (*miit)->th_FieldDataVerts, &*eit, 1, &ret_val,
671  &tag_size);
672  else
673  CHKERR get_moab().tag_get_by_ptr((*miit)->th_FieldData,
674  &*eit, 1, &ret_val,
675  &tag_size);
676  const_cast<FieldData *&>(*dit) = cast(ret_val);
677  }
678  }
679  }
680  }
681  return vec;
682  };
683 
684  auto ents_in_ref_ent = get_ents_in_ref_ent(miit_ref_ent);
685 
686  CHKERR create_tags_for_max_order(ents_in_ref_ent);
687  CHKERR create_tags_for_data(ents_in_ref_ent);
688  auto ents_max_order = get_ents_max_order(ents_in_ref_ent);
689  auto ent_field_data =
690  get_ents_field_data_vector_adaptor(ents_in_ref_ent, ents_max_order);
691 
692  auto vit_max_order = ents_max_order->begin();
693  auto vit_field_data = ent_field_data->begin();
694  for (auto ent : ents_in_ref_ent) {
695  ents_array->emplace_back(
696  *miit, *miit_ref_ent,
697  boost::shared_ptr<double *const>(ent_field_data, &*vit_field_data),
698  boost::shared_ptr<const int>(
699  ents_max_order, static_cast<const int *>(*vit_max_order)));
700  ++miit_ref_ent;
701  ++vit_max_order;
702  ++vit_field_data;
703  }
704  nb_ents_set_order_new += ents_in_ref_ent.size();
705 
706  // Check if any of entities in the range has bit level but is not added
707  // to database. That generate data inconsistency and error.
708  if (ents_in_ref_ent.size() < (second - first + 1)) {
709  Range ents_not_in_database =
710  subtract(Range(first, second), ents_in_ref_ent);
711  std::vector<const void *> vec_bits(ents_not_in_database.size());
712  CHKERR get_moab().tag_get_by_ptr(
713  get_basic_entity_data_ptr()->th_RefBitLevel, ents_not_in_database,
714  &*vec_bits.begin());
715  auto cast = [](auto p) { return static_cast<const BitRefLevel *>(p); };
716  for (auto v : vec_bits)
717  if (cast(v)->any())
718  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
719  "Try to add entities which are not seeded or added to "
720  "database");
721  }
722 
723  // Add entities to database
724  auto hint = entsFields.end();
725  for (auto &v : *ents_array) {
726  hint = entsFields.emplace_hint(hint, ents_array, &v);
727  }
728  }
729  }
730 
731  if (verb >= VERY_VERBOSE) {
732  PetscSynchronizedPrintf(
733  cOmm,
734  "nb. of entities in field <%s> for which order was "
735  "increased %d (order %d)\n",
736  miit->get()->getName().c_str(), nb_ents_set_order_up, order);
737  PetscSynchronizedPrintf(
738  cOmm,
739  "nb. of entities in field <%s> for which order was "
740  "reduced %d (order %d)\n",
741  miit->get()->getName().c_str(), nb_ents_set_order_down, order);
742  PetscSynchronizedPrintf(
743  cOmm,
744  "nb. of entities in field <%s> for which order set %d (order %d)\n",
745  miit->get()->getName().c_str(), nb_ents_set_order_new, order);
746  PetscSynchronizedFlush(cOmm, PETSC_STDOUT);
747  }
748 
750 }
752  const EntityType type, const BitFieldId id,
753  const ApproximationOrder order, int verb) {
755  if (verb == -1)
756  verb = verbose;
757  *buildMoFEM = 0;
758  Range ents;
759  CHKERR get_moab().get_entities_by_type(meshset, type, ents);
760  if (verb > VERBOSE) {
761  PetscSynchronizedPrintf(cOmm, "nb. of ents for order change %d\n",
762  ents.size());
763  }
764  CHKERR set_field_order(ents, id, order, verb);
765  if (verb > VERBOSE) {
766  PetscSynchronizedFlush(cOmm, PETSC_STDOUT);
767  }
769 }
771  const EntityType type,
772  const std::string &name,
773  const ApproximationOrder order, int verb) {
775  if (verb == -1)
776  verb = verbose;
777  *buildMoFEM = 0;
778  CHKERR set_field_order(meshset, type, getBitFieldId(name), order, verb);
780 }
781 MoFEMErrorCode Core::set_field_order(const Range &ents, const std::string &name,
782  const ApproximationOrder order, int verb) {
784  if (verb == -1)
785  verb = verbose;
786  *buildMoFEM = 0;
787  CHKERR set_field_order(ents, getBitFieldId(name), order, verb);
789 }
791  const BitRefLevel &bit, const BitRefLevel &mask, const EntityType type,
792  const BitFieldId id, const ApproximationOrder order, int verb) {
794  if (verb == -1)
795  verb = verbose;
796  *buildMoFEM = 0;
797  Range ents;
799  ents, verb);
800  CHKERR set_field_order(ents, id, order, verb);
802 }
803 
805  const BitRefLevel &bit, const BitRefLevel &mask, const EntityType type,
806  const std::string &name, const ApproximationOrder order, int verb) {
808  if (verb == -1)
809  verb = verbose;
810  *buildMoFEM = 0;
811  Range ents;
813  ents, verb);
814  CHKERR set_field_order(ents, getBitFieldId(name), order, verb);
816 }
817 
820  std::map<EntityType, int> &dof_counter, int verb) {
822  if (verb == -1)
823  verb = verbose;
824  // field it
825  auto &set_id = fIelds.get<BitFieldId_mi_tag>();
826  // find fields
827  auto miit = set_id.find(id);
828  if (miit == set_id.end())
829  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "field not found");
830 
831  // ents in the field meshset
832  Range ents_of_id_meshset;
833  CHKERR get_moab().get_entities_by_handle((*miit)->meshSet, ents_of_id_meshset,
834  false);
835  if (verb > VERY_NOISY)
836  PetscSynchronizedPrintf(cOmm, "ents in field %s meshset %d\n",
837  (*miit)->getName().c_str(),
838  ents_of_id_meshset.size());
839 
840  // ent view by field id (in set all MoabEnts has the same FieldId)
841  auto eiit =
842  entsFields.get<FieldName_mi_tag>().lower_bound(miit->get()->getNameRef());
843  FieldEntity_multiIndex_ent_view ents_id_view;
844  if (eiit != entsFields.get<FieldName_mi_tag>().end()) {
845  auto hi_eiit = entsFields.get<FieldName_mi_tag>().upper_bound(
846  miit->get()->getNameRef());
847  std::copy(eiit, hi_eiit, std::back_inserter(ents_id_view));
848  }
849 
850  boost::shared_ptr<const int> zero_order(new const int(0));
851 
852  for (Range::iterator eit = ents_of_id_meshset.begin();
853  eit != ents_of_id_meshset.end(); eit++) {
854  // search if field meshset is in database
855  RefEntity_multiIndex::index<Ent_mi_tag>::type::iterator miit_ref_ent;
856  miit_ref_ent = refinedEntities.get<Ent_mi_tag>().find(*eit);
857  if (miit_ref_ent == refinedEntities.get<Ent_mi_tag>().end()) {
858  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
859  "Entity is not in MoFEM database, entities in field meshset need "
860  "to be seeded (i.e. bit ref level add to them)");
861  }
862 
863  auto add_dofs = [&](auto field_eit) {
865  // create dofs on this entity (nb. of dofs is equal to rank)
866  for (FieldCoefficientsNumber rank = 0; rank < (*miit)->getNbOfCoeffs();
867  rank++) {
868  std::pair<DofEntity_multiIndex::iterator, bool> d_miit;
869  // insert dof
870  d_miit = dofsField.insert(
871  boost::make_shared<DofEntity>(field_eit, 0, rank, rank, true));
872  if (d_miit.second) {
873  dof_counter[MBENTITYSET]++; // Count entities in the meshset
874  } else {
875  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
876  "Dof should be created");
877  }
878  }
880  };
881 
882  // create database entity
883  auto field_eit = ents_id_view.get<1>().find(*eit);
884  if (field_eit == ents_id_view.get<1>().end()) {
885 
886  std::pair<FieldEntity_multiIndex::iterator, bool> e_miit;
887  e_miit = entsFields.insert(boost::make_shared<FieldEntity>(
888  *miit, *miit_ref_ent,
889  FieldEntity::makeSharedFieldDataAdaptorPtr(*miit, *miit_ref_ent),
890  boost::shared_ptr<const int>(zero_order, zero_order.get())));
891 
892  if (!e_miit.second)
893  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
894  "Entity should be created");
895 
896  CHKERR add_dofs(*(e_miit.first));
897 
898 
899  } else {
900 
901  // If there are DOFs in that range is more pragmatic to remove them
902  // rather than to find sub-ranges or make them inactive
903  auto dit = dofsField.get<Composite_Name_And_Ent_mi_tag>().lower_bound(
904  boost::make_tuple((*miit)->getNameRef(), *eit));
905  auto hi_dit = dofsField.get<Composite_Name_And_Ent_mi_tag>().upper_bound(
906  boost::make_tuple((*miit)->getNameRef(), *eit));
907  dofsField.get<Composite_Name_And_Ent_mi_tag>().erase(dit, hi_dit);
908 
909  CHKERR add_dofs(*field_eit);
910  }
911  }
912 
913  if (verb > VERY_VERBOSE) {
914  auto lo_dof = dofsField.get<FieldName_mi_tag>().lower_bound(
915  miit->get()->getNameRef());
916  auto hi_dof = dofsField.get<FieldName_mi_tag>().upper_bound(
917  miit->get()->getNameRef());
918  for (; lo_dof != hi_dof; lo_dof++) {
919  std::ostringstream ss;
920  ss << *lo_dof << std::endl;
921  PetscSynchronizedPrintf(cOmm, ss.str().c_str());
922  }
923  PetscSynchronizedFlush(cOmm, PETSC_STDOUT);
924  }
925 
927 }
928 
930  const BitFieldId id, std::map<EntityType, int> &dof_counter,
931  std::map<EntityType, int> &inactive_dof_counter, int verb) {
933  if (verb == -1)
934  verb = verbose;
935 
936  // Find field
937  auto &set_id = fIelds.get<BitFieldId_mi_tag>();
938  auto field_it = set_id.find(id);
939  if (field_it == set_id.end()) {
940  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "Field not found");
941  }
942  const int rank = field_it->get()->getNbOfCoeffs();
943  const boost::string_ref &field_name = field_it->get()->getNameRef();
944 
945  // Ents in the field meshset
946  Range ents_of_id_meshset;
947  CHKERR get_moab().get_entities_by_handle((*field_it)->meshSet,
948  ents_of_id_meshset, false);
949  if (verb > VERY_NOISY) {
950  PetscSynchronizedPrintf(PETSC_COMM_SELF, "Ents in field %s meshset %d\n",
951  (*field_it)->getName().c_str(),
952  ents_of_id_meshset.size());
953  }
954 
955  for (auto p_eit = ents_of_id_meshset.pair_begin();
956  p_eit != ents_of_id_meshset.pair_end(); ++p_eit) {
957 
958  const EntityHandle first = p_eit->first;
959  const EntityHandle second = p_eit->second;
960 
961  auto feit = entsFields.get<Composite_Name_And_Ent_mi_tag>().lower_bound(
962  boost::make_tuple(field_name, first));
963  if (feit == entsFields.get<Composite_Name_And_Ent_mi_tag>().end())
964  continue;
965  auto hi_feit = entsFields.get<Composite_Name_And_Ent_mi_tag>().upper_bound(
966  boost::make_tuple(field_name, second));
967 
968  // If there are DOFs in that range is more pragmatic to remove them
969  // rather than to find sub-ranges or make them inactive
970  auto dit = dofsField.get<Composite_Name_And_Ent_mi_tag>().lower_bound(
971  boost::make_tuple(field_name, first));
972  auto hi_dit = dofsField.get<Composite_Name_And_Ent_mi_tag>().upper_bound(
973  boost::make_tuple(field_name, second));
974  dofsField.get<Composite_Name_And_Ent_mi_tag>().erase(dit, hi_dit);
975 
976  // Add vertices DOFs by bulk
977  boost::shared_ptr<std::vector<DofEntity>> dofs_array =
978  boost::make_shared<std::vector<DofEntity>>(std::vector<DofEntity>());
979  // Add Sequence of DOFs to sequence container as weak_ptr
980  int nb_dofs_on_ents = 0;
981  for (auto tmp_feit = feit; tmp_feit != hi_feit; ++tmp_feit) {
982  nb_dofs_on_ents += rank * tmp_feit->get()->getOrderNbDofs(
983  tmp_feit->get()->getMaxOrder());
984  }
985 
986  // Reserve memory
987  dofs_array->reserve(nb_dofs_on_ents);
988 
989  // Create DOFs
990  for (; feit != hi_feit; ++feit) {
991  // Create dofs instances and shared pointers
992  int DD = 0;
993  // Loop orders (loop until max entity order is set)
994  for (int oo = 0; oo <= feit->get()->getMaxOrder(); ++oo) {
995  // Loop nb. dofs at order oo
996  for (int dd = 0; dd < feit->get()->getOrderNbDofsDiff(oo); ++dd) {
997  // Loop rank
998  for (int rr = 0; rr < rank; ++rr, ++DD) {
999  dofs_array->emplace_back(*feit, oo, rr, DD, true);
1000  ++dof_counter[feit->get()->getEntType()];
1001  }
1002  }
1003  }
1004  if (DD > feit->get()->getNbDofsOnEnt()) {
1005  std::ostringstream ss;
1006  ss << "rank " << rAnk << " ";
1007  ss << **feit << std::endl;
1008  SETERRQ3(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1009  "Expected number of DOFs on entity not equal to number added "
1010  "to database (DD = %d != %d = "
1011  "feit->get()->getNbDofsOnEnt())\n"
1012  "%s",
1013  DD, feit->get()->getNbDofsOnEnt(), ss.str().c_str());
1014  }
1015  }
1016 
1017  // Insert into Multi-Index container
1018  int dofs_field_size0 = dofsField.size();
1019  auto hint = dofsField.end();
1020  for (auto &v : *dofs_array)
1021  hint = dofsField.emplace_hint(hint, dofs_array, &v);
1022 
1023  // Add Sequence of DOFs to sequence container as weak_ptr
1024  field_it->get()->getDofSequenceContainer().push_back(dofs_array);
1025 
1026  // Check data consistency
1027  if (PetscUnlikely(static_cast<int>(dofs_array.use_count()) !=
1028  static_cast<int>(dofs_array->size() + 1))) {
1029  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1030  "Wrong use count %d != %d", dofs_array.use_count(),
1031  dofs_array->size() + 1);
1032  }
1033  if (dofs_field_size0 + dofs_array->size() != dofsField.size()) {
1034  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1035  "Wrong number of inserted DOFs %d != %d", dofs_array->size(),
1036  dofsField.size() - dofs_field_size0);
1037  }
1038  }
1040 }
1041 
1042 MoFEMErrorCode Core::buildField(const boost::shared_ptr<Field> &field,
1043  int verb) {
1045  if (verb == -1)
1046  verb = verbose;
1047  if (verb > QUIET) {
1048  PetscSynchronizedPrintf(cOmm, "Build Field %s (rank %d)\n",
1049  field->getName().c_str(), rAnk);
1050  }
1051  std::map<EntityType, int> dof_counter;
1052  std::map<EntityType, int> inactive_dof_counter;
1053 
1054  // Need to rebuild order table since number of dofs on each order when
1055  // field was created.
1056  if (field->getApproxBase() == USER_BASE)
1057  CHKERR field->rebuildDofsOrderMap();
1058 
1059  switch (field->getSpace()) {
1060  case NOFIELD:
1061  CHKERR buildFieldForNoField(field->getId(), dof_counter, verb);
1062  break;
1063  case L2:
1064  case H1:
1065  case HCURL:
1066  case HDIV:
1067  CHKERR buildFieldForL2H1HcurlHdiv(field->getId(), dof_counter,
1068  inactive_dof_counter, verb);
1069  break;
1070  default:
1071  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
1072  }
1073 
1074  if (verb > QUIET) {
1075  int nb_added_dofs = 0;
1076  int nb_inactive_added_dofs = 0;
1077  for (auto const &it : dof_counter) {
1078  switch (it.first) {
1079  case MBVERTEX:
1080  PetscSynchronizedPrintf(cOmm,
1081  "nb added dofs (vertices) %d (inactive %d)\n",
1082  it.second, inactive_dof_counter[it.first]);
1083  break;
1084  case MBEDGE:
1085  PetscSynchronizedPrintf(cOmm,
1086  "nb added dofs (edges) %d (inactive %d)\n",
1087  it.second, inactive_dof_counter[it.first]);
1088  break;
1089  case MBTRI:
1090  PetscSynchronizedPrintf(cOmm,
1091  "nb added dofs (triangles) %d (inactive %d)\n",
1092  it.second, inactive_dof_counter[it.first]);
1093  break;
1094  case MBQUAD:
1095  PetscSynchronizedPrintf(cOmm,
1096  "nb added dofs (quads) %d (inactive %d)\n",
1097  it.second, inactive_dof_counter[it.first]);
1098  break;
1099  case MBTET:
1100  PetscSynchronizedPrintf(cOmm, "nb added dofs (tets) %d (inactive %d)\n",
1101  it.second, inactive_dof_counter[it.first]);
1102  break;
1103  case MBPRISM:
1104  PetscSynchronizedPrintf(cOmm,
1105  "nb added dofs (prisms) %d (inactive %d)\n",
1106  it.second, inactive_dof_counter[it.first]);
1107  break;
1108  case MBENTITYSET:
1109  PetscSynchronizedPrintf(cOmm,
1110  "nb added dofs (meshsets) %d (inactive %d)\n",
1111  it.second, inactive_dof_counter[it.first]);
1112  break;
1113  default:
1114  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
1115  }
1116  nb_added_dofs += it.second;
1117  nb_inactive_added_dofs += inactive_dof_counter[it.first];
1118  }
1119  if (verb > QUIET) {
1120  PetscSynchronizedPrintf(cOmm,
1121  "nb added dofs %d (number of inactive dofs %d)\n",
1122  nb_added_dofs, nb_inactive_added_dofs);
1123  }
1124  }
1126 }
1127 
1128 MoFEMErrorCode Core::build_field(const std::string field_name, int verb) {
1130  auto miit = fIelds.get<FieldName_mi_tag>().find(field_name);
1131  if (miit == fIelds.get<FieldName_mi_tag>().end()) {
1132  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "Field < %s > not found",
1133  field_name.c_str());
1134  }
1135  CHKERR buildField((*miit), verb);
1136  PetscSynchronizedFlush(cOmm, PETSC_STDOUT);
1138 }
1139 
1142  if (verb == -1)
1143  verb = verbose;
1144  auto &set_id = fIelds.get<BitFieldId_mi_tag>();
1145  for (auto miit = set_id.begin(); miit != set_id.end(); miit++) {
1146  CHKERR buildField((*miit), verb);
1147  }
1148  *buildMoFEM = 1 << 0;
1149  if (verb > QUIET) {
1150  PetscSynchronizedPrintf(cOmm, "Nb. dofs %u\n", dofsField.size());
1151  }
1152  PetscSynchronizedFlush(cOmm, PETSC_STDOUT);
1154  // return 0;
1155 }
1156 
1158 Core::list_dofs_by_field_name(const std::string &field_name) const {
1160  auto dit = dofsField.get<FieldName_mi_tag>().lower_bound(field_name);
1161  auto hi_dit = dofsField.get<FieldName_mi_tag>().upper_bound(field_name);
1162  for (; dit != hi_dit; dit++) {
1163  std::ostringstream ss;
1164  ss << "rank " << rAnk << " ";
1165  ss << *dit << std::endl;
1166  PetscSynchronizedPrintf(cOmm, ss.str().c_str());
1167  }
1168  PetscSynchronizedFlush(cOmm, PETSC_STDOUT);
1170 }
1173  auto &set_id = fIelds.get<BitFieldId_mi_tag>();
1174  for (auto &miit : set_id) {
1175  std::ostringstream ss;
1176  ss << *miit << std::endl;
1177  PetscSynchronizedPrintf(cOmm, ss.str().c_str());
1178  }
1179  PetscSynchronizedFlush(cOmm, PETSC_STDOUT);
1181 }
1182 
1185  for (auto &miit : entFEAdjacencies) {
1186  std::ostringstream ss;
1187  ss << miit << std::endl;
1188  PetscSynchronizedPrintf(cOmm, ss.str().c_str());
1189  }
1190  PetscSynchronizedFlush(cOmm, PETSC_STDOUT);
1192 }
1193 
1195 Core::get_problem_finite_elements_entities(const std::string &problem_name,
1196  const std::string &fe_name,
1197  const EntityHandle meshset) {
1199  auto &prb = pRoblems.get<Problem_mi_tag>();
1200  auto p_miit = prb.find(problem_name);
1201  if (p_miit == prb.end())
1202  SETERRQ1(PETSC_COMM_SELF, 1, "no such problem like < %s >",
1203  problem_name.c_str());
1204  auto miit = p_miit->numeredFiniteElements->get<FiniteElement_name_mi_tag>()
1205  .lower_bound(fe_name);
1206  auto hi_miit = p_miit->numeredFiniteElements->get<FiniteElement_name_mi_tag>()
1207  .upper_bound(fe_name);
1208  for (; miit != hi_miit; miit++) {
1209  EntityHandle ent = (*miit)->getEnt();
1210  CHKERR get_moab().add_entities(meshset, &ent, 1);
1211  int part = (*miit)->getPart();
1212  CHKERR get_moab().tag_set_data(th_Part, &ent, 1, &part);
1213  }
1215 }
1216 
1217 FieldEntityByFieldName::iterator
1218 Core::get_ent_field_by_name_begin(const std::string &field_name) const {
1219  return entsFields.get<FieldName_mi_tag>().lower_bound(field_name);
1220 }
1221 FieldEntityByFieldName::iterator
1222 Core::get_ent_field_by_name_end(const std::string &field_name) const {
1223  return entsFields.get<FieldName_mi_tag>().upper_bound(field_name);
1224 }
1225 DofEntityByFieldName::iterator
1226 Core::get_dofs_by_name_begin(const std::string &field_name) const {
1227  return dofsField.get<FieldName_mi_tag>().lower_bound(field_name);
1228 }
1229 DofEntityByFieldName::iterator
1230 Core::get_dofs_by_name_end(const std::string &field_name) const {
1231  return dofsField.get<FieldName_mi_tag>().upper_bound(field_name);
1232 }
1233 DofEntityByNameAndEnt::iterator
1234 Core::get_dofs_by_name_and_ent_begin(const std::string &field_name,
1235  const EntityHandle ent) const {
1236  return dofsField.get<Composite_Name_And_Ent_mi_tag>().lower_bound(
1237  boost::make_tuple(field_name, ent));
1238 }
1239 DofEntityByNameAndEnt::iterator
1240 Core::get_dofs_by_name_and_ent_end(const std::string &field_name,
1241  const EntityHandle ent) const {
1242  return dofsField.get<Composite_Name_And_Ent_mi_tag>().upper_bound(
1243  boost::make_tuple(field_name, ent));
1244 }
1245 DofEntityByNameAndType::iterator
1246 Core::get_dofs_by_name_and_type_begin(const std::string &field_name,
1247  const EntityType type) const {
1248  return dofsField.get<Composite_Name_And_Type_mi_tag>().lower_bound(
1249  boost::make_tuple(field_name, type));
1250 }
1251 DofEntityByNameAndType::iterator
1252 Core::get_dofs_by_name_and_type_end(const std::string &field_name,
1253  const EntityType type) const {
1254  return dofsField.get<Composite_Name_And_Type_mi_tag>().upper_bound(
1255  boost::make_tuple(field_name, type));
1256 }
1258 Core::check_number_of_ents_in_ents_field(const std::string &name) const {
1260  auto it = fIelds.get<FieldName_mi_tag>().find(name);
1261  if (it == fIelds.get<FieldName_mi_tag>().end()) {
1262  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1263  "field not found < %s >", name.c_str());
1264  }
1265  EntityHandle meshset = (*it)->getMeshset();
1266  int num_entities;
1267  CHKERR get_moab().get_number_entities_by_handle(meshset, num_entities);
1268  if (entsFields.get<FieldName_mi_tag>().count((*it)->getName()) >
1269  (unsigned int)num_entities) {
1270  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1271  "not equal number of entities in meshset and field multiindex "
1272  "< %s >",
1273  name.c_str());
1274  }
1276 }
1279  for (auto &it : fIelds.get<FieldName_mi_tag>()) {
1280  if (it->getSpace() == NOFIELD)
1281  continue; // FIXME: should be treated properly, not test is just
1282  // skipped for this NOFIELD space
1283  EntityHandle meshset = it->getMeshset();
1284  int num_entities;
1285  CHKERR get_moab().get_number_entities_by_handle(meshset, num_entities);
1286  if (entsFields.get<FieldName_mi_tag>().count(it->getName()) >
1287  (unsigned int)num_entities) {
1288  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1289  "not equal number of entities in meshset and field "
1290  "multiindex < %s >",
1291  it->getName().c_str());
1292  }
1293  }
1295 }
1296 
1297 } // namespace MoFEM
DofEntity_multiIndex dofsField
dofs on fields
Definition: Core.hpp:249
Tag th_FieldSpace
Definition: Core.hpp:204
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:1195
Tag th_FieldName_DataNamePrefix
Definition: Core.hpp:204
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:320
field with continuous normal traction
Definition: definitions.h:178
user implemented approximation base
Definition: definitions.h:159
DofEntityByNameAndEnt::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:1234
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:358
MPI_Comm cOmm
MoFEM communicator.
Definition: Core.hpp:846
RefEntity_multiIndex refinedEntities
refined entities
Definition: Core.hpp:244
Provide data structure for (tensor) field approximation.The Field is intended to provide support for ...
FieldEntity_multiIndex entsFields
entities on fields
Definition: Core.hpp:248
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
MoFEMErrorCode check_number_of_ents_in_ents_field() const
check data consistency in entitiesPtr
Definition: FieldCore.cpp:1277
int ApproximationOrder
Approximation on the entity.
Definition: Types.hpp:37
moab::Interface & get_moab()
Definition: Core.hpp:266
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:506
scalar or vector of scalars describe (no true field)
Definition: definitions.h:175
DofEntityByNameAndType::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:1246
DofEntityByFieldName::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:1230
static boost::shared_ptr< FieldData *const > makeSharedFieldDataAdaptorPtr(const boost::shared_ptr< Field > &field_ptr, const boost::shared_ptr< RefEntity > &ref_ent_ptr)
Return shared pointer to entity field data vector adaptor.
FTensor::Index< 'n', 2 > n
Definition: PlasticOps.hpp:68
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
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:313
const EntityHandle no_handle
No entity handle is indicated by zero handle, i.e. root meshset.
Definition: Common.hpp:27
#define THROW_MESSAGE(a)
Throw MoFEM exception.
Definition: definitions.h:625
Tag th_Part
Tag for partition number.
Definition: Core.hpp:200
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:513
MoFEMErrorCode list_dofs_by_field_name(const std::string &name) const
Definition: FieldCore.cpp:1158
std::reference_wrapper< moab::Interface > moab
moab database
Definition: Core.hpp:265
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:64
MoFEMErrorCode list_fields() const
list entities in the field
Definition: FieldCore.cpp:1171
const Field * get_field_structure(const std::string &name)
get field structure
Definition: FieldCore.cpp:51
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
Tag th_FieldId
Definition: Core.hpp:204
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:84
DofEntityByFieldName::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:1226
int rAnk
MOFEM communicator rank.
Definition: Core.hpp:850
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:91
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
FieldEntityByFieldName::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:1222
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
FieldApproximationBase
approximation base
Definition: definitions.h:149
MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)
Definition: FieldCore.cpp:1140
int * buildMoFEM
keeps flags/semaphores for different stages
Definition: Core.hpp:896
constexpr int order
Tag th_FieldBase
Definition: Core.hpp:204
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
FieldEntityEntFiniteElementAdjacencyMap_multiIndex entFEAdjacencies
adjacencies of elements to dofs
Definition: Core.hpp:255
MoFEMErrorCode build_field(const std::string field_name, int verb=DEFAULT_VERBOSITY)
build field by name
Definition: FieldCore.cpp:1128
Managing BitRefLevels.
Problem_multiIndex pRoblems
problems multi-index
Definition: Core.hpp:257
field with continuous tangents
Definition: definitions.h:177
multi_index_container< boost::shared_ptr< FieldEntity >, indexed_by< sequenced<>, ordered_non_unique< tag< Ent_mi_tag >, const_mem_fun< FieldEntity, EntityHandle, &FieldEntity::getEnt > > > > FieldEntity_multiIndex_ent_view
Tag th_FieldName
Definition: Core.hpp:204
MoFEMErrorCode get_field_entities_by_handle(const std::string name, Range &ents) const
get entities in the field by handle
Definition: FieldCore.cpp:83
FieldSpace
approximation spaces
Definition: definitions.h:173
std::bitset< BITFIELDID_SIZE > BitFieldId
Field Id.
Definition: Types.hpp:53
#define CHKERR
Inline error check.
Definition: definitions.h:601
Field_multiIndex fIelds
fields
Definition: Core.hpp:247
MoFEMErrorCode buildField(const boost::shared_ptr< Field > &field, int verb=DEFAULT_VERBOSITY)
Definition: FieldCore.cpp:1042
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
MoFEMErrorCode addEntsToFieldByDim(const Range &ents, const int dim, const std::string &name, int verb=DEFAULT_VERBOSITY)
Definition: FieldCore.cpp:234
MultiIndex Tag for field name.
MoFEMErrorCode set_field_order(const Range &ents, const BitFieldId id, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)
Definition: FieldCore.cpp:397
BitFieldId getBitFieldId(const std::string &name) const
Definition: FieldCore.cpp:21
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:929
int FieldCoefficientsNumber
Number of field coefficients.
Definition: Types.hpp:38
int verbose
Verbosity level.
Definition: Core.hpp:885
MoFEMErrorCode buildFieldForNoField(const BitFieldId id, std::map< EntityType, int > &dof_counter, int verb=DEFAULT_VERBOSITY)
Definition: FieldCore.cpp:819
FieldEntityByFieldName::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:1218
EntityHandle get_field_meshset(const BitFieldId id) const
Definition: FieldCore.cpp:31
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412
continuous field
Definition: definitions.h:176
Tag th_RefBitLevel
Definition: Core.hpp:201
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:72
boost::shared_ptr< BasicEntityData > & get_basic_entity_data_ptr()
Get pointer to basic entity data.
Definition: Core.hpp:234
MoFEMErrorCode list_adjacencies() const
list adjacencies
Definition: FieldCore.cpp:1183
DofEntityByNameAndType::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:1252
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:790
DofEntityByNameAndEnt::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:1240
bool check_field(const std::string &name) const
check if field is in database
Definition: FieldCore.cpp:43
MoFEMTypes
Those types control how functions respond on arguments, f.e. error handling.
Definition: definitions.h:188
BitFieldId getFieldShift()
Return unique field Id.
Definition: Core.cpp:226
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:74
FTensor::Index< 'i', 2 > i
[Common data]
Definition: ContactOps.hpp:26
field with C-1 continuity
Definition: definitions.h:179