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