v0.13.0
Simple.cpp
Go to the documentation of this file.
1 /** \file Simple.cpp
2  * \brief Implementation of simple interface
3  * \ingroup mofem_simple_interface
4  */
5 
6 /* MoFEM is free software: you can redistribute it and/or modify it under
7  * the terms of the GNU Lesser General Public License as published by the
8  * Free Software Foundation, either version 3 of the License, or (at your
9  * option) any later version.
10  *
11  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
12  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
14  * License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public
17  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>
18  */
19 
20 namespace MoFEM {
21 
22 MoFEMErrorCode Simple::query_interface(boost::typeindex::type_index type_index,
23  UnknownInterface **iface) const {
24  *iface = const_cast<Simple *>(this);
25  return 0;
26 }
27 
29  static_assert(DIM == 2 || DIM == 3, "not implemented");
30  return MOFEM_NOT_IMPLEMENTED;
31 }
32 
33 template <> MoFEMErrorCode Simple::setSkeletonAdjacency<2>() {
34  Interface &m_field = cOre;
36 
37  auto defaultSkeletonEdge =
38  [&](moab::Interface &moab, const Field &field, const EntFiniteElement &fe,
39  std::vector<EntityHandle> &adjacency) -> MoFEMErrorCode {
41 
42  CHKERR DefaultElementAdjacency::defaultEdge(moab, field, fe, adjacency);
43 
44  if (std::find(domainFields.begin(), domainFields.end(), field.getName()) !=
45  domainFields.end()) {
46 
47  const EntityHandle fe_ent = fe.getEnt();
48  std::vector<EntityHandle> bride_adjacency_edge;
49  CHKERR moab.get_adjacencies(&fe_ent, 1, 2, false, bride_adjacency_edge);
50 
51  switch (field.getSpace()) {
52  case H1:
53  CHKERR moab.get_connectivity(&*bride_adjacency_edge.begin(),
54  bride_adjacency_edge.size(), adjacency,
55  true);
56  case HCURL:
57  case HDIV:
58  CHKERR moab.get_adjacencies(&*bride_adjacency_edge.begin(),
59  bride_adjacency_edge.size(), 1, false,
60  adjacency, moab::Interface::UNION);
61  case L2:
62  adjacency.insert(adjacency.end(), bride_adjacency_edge.begin(),
63  bride_adjacency_edge.end());
64  break;
65  case NOFIELD:
66  break;
67  default:
68  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
69  "this field is not implemented for TRI finite element");
70  }
71 
72  std::sort(adjacency.begin(), adjacency.end());
73  auto it = std::unique(adjacency.begin(), adjacency.end());
74 
75  std::vector<EntityHandle> new_adjacency(
76  std::distance(adjacency.begin(), it));
77  std::copy(adjacency.begin(), it, new_adjacency.begin());
78 
79  for (auto e : new_adjacency) {
80  auto side_table = fe.getSideNumberTable();
81  if (side_table.find(e) == side_table.end())
82  const_cast<SideNumber_multiIndex &>(fe.getSideNumberTable())
83  .insert(
84  boost::shared_ptr<SideNumber>(new SideNumber(e, -1, 0, 0)));
85  }
86 
87  adjacency.swap(new_adjacency);
88  }
89 
91  };
92 
94  defaultSkeletonEdge);
95 
97 }
98 
99 template <> MoFEMErrorCode Simple::setSkeletonAdjacency<3>() {
100  Interface &m_field = cOre;
102 
103  auto defaultSkeletonEdge =
104  [&](moab::Interface &moab, const Field &field, const EntFiniteElement &fe,
105  std::vector<EntityHandle> &adjacency) -> MoFEMErrorCode {
107 
108  CHKERR DefaultElementAdjacency::defaultFace(moab, field, fe, adjacency);
109 
110  if (std::find(domainFields.begin(), domainFields.end(), field.getName()) !=
111  domainFields.end()) {
112 
113  const EntityHandle fe_ent = fe.getEnt();
114  std::vector<EntityHandle> bride_adjacency_edge;
115  CHKERR moab.get_adjacencies(&fe_ent, 1, 2, false, bride_adjacency_edge);
116 
117  switch (field.getSpace()) {
118  case H1:
119  CHKERR moab.get_connectivity(&*bride_adjacency_edge.begin(),
120  bride_adjacency_edge.size(), adjacency,
121  true);
122  case HCURL:
123  CHKERR moab.get_adjacencies(&*bride_adjacency_edge.begin(),
124  bride_adjacency_edge.size(), 1, false,
125  adjacency, moab::Interface::UNION);
126  case HDIV:
127  CHKERR moab.get_adjacencies(&*bride_adjacency_edge.begin(),
128  bride_adjacency_edge.size(), 2, false,
129  adjacency, moab::Interface::UNION);
130  case L2:
131  adjacency.insert(adjacency.end(), bride_adjacency_edge.begin(),
132  bride_adjacency_edge.end());
133  break;
134  case NOFIELD:
135  break;
136  default:
137  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
138  "this field is not implemented for TRI finite element");
139  }
140 
141  std::sort(adjacency.begin(), adjacency.end());
142  auto it = std::unique(adjacency.begin(), adjacency.end());
143 
144  std::vector<EntityHandle> new_adjacency(
145  std::distance(adjacency.begin(), it));
146  std::copy(adjacency.begin(), it, new_adjacency.begin());
147 
148  for (auto e : new_adjacency) {
149  auto side_table = fe.getSideNumberTable();
150  if (side_table.find(e) == side_table.end())
151  const_cast<SideNumber_multiIndex &>(fe.getSideNumberTable())
152  .insert(
153  boost::shared_ptr<SideNumber>(new SideNumber(e, -1, 0, 0)));
154  }
155 
156  adjacency.swap(new_adjacency);
157  }
158 
160  };
161 
163  defaultSkeletonEdge);
165  defaultSkeletonEdge);
166 
168 }
169 
172  switch (getDim()) {
173  case 1:
174  THROW_MESSAGE("Not implemented");
175  case 2:
176  return setSkeletonAdjacency<2>();
177  case 3:
178  return setSkeletonAdjacency<3>();
179  default:
180  THROW_MESSAGE("Not implemented");
181  }
183 }
184 
187  if (dim == -1)
188  dim = getDim();
189  switch (dim) {
190  case 2:
191  return setSkeletonAdjacency<2>();
192  case 3:
193  return setSkeletonAdjacency<3>();
194  default:
195  SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED, "Not implemented");
196  }
198 }
199 
200 Simple::Simple(const Core &core)
201  : cOre(const_cast<Core &>(core)), bitLevel(BitRefLevel().set(0)),
202  bitLevelMask(BitRefLevel().set()), meshSet(0), boundaryMeshset(0),
203  skeletonMeshset(0), nameOfProblem("SimpleProblem"), domainFE("dFE"),
204  boundaryFE("bFE"), skeletonFE("sFE"), dIm(-1), addSkeletonFE(false),
205  addBoundaryFE(false) {
206  PetscLogEventRegister("SimpleSetUp", 0, &MOFEM_EVENT_SimpleSetUP);
207  PetscLogEventRegister("SimpleLoadMesh", 0, &MOFEM_EVENT_SimpleLoadMesh);
208  PetscLogEventRegister("SimpleBuildFields", 0, &MOFEM_EVENT_SimpleBuildFields);
209  PetscLogEventRegister("SimpleBuildFiniteElements", 0,
211  PetscLogEventRegister("SimpleSetUp", 0, &MOFEM_EVENT_SimpleBuildProblem);
212  PetscLogEventRegister("SimpleKSPSolve", 0, &MOFEM_EVENT_SimpleKSPSolve);
213  strcpy(meshFileName, "mesh.h5m");
214 }
215 
217  PetscBool flg = PETSC_TRUE;
219  ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "", "Simple interface options",
220  "none");
221  CHKERRG(ierr);
222  ierr = PetscOptionsString("-file_name", "file name", "", "mesh.h5m",
223  meshFileName, 255, &flg);
224  CHKERRG(ierr);
225  ierr = PetscOptionsEnd();
226  CHKERRG(ierr);
228 }
229 
230 MoFEMErrorCode Simple::loadFile(const std::string options,
231  const std::string mesh_file_name) {
232  Interface &m_field = cOre;
234  PetscLogEventBegin(MOFEM_EVENT_SimpleLoadMesh, 0, 0, 0, 0);
235 
236  if (!mesh_file_name.empty())
237  strcpy(meshFileName, mesh_file_name.c_str());
238 
239  // This is a case of distributed mesh and algebra. In that case each processor
240  // keep only part of the problem.
241  CHKERR m_field.get_moab().load_file(meshFileName, 0, options.c_str());
242  CHKERR m_field.rebuild_database();
243 
244  // determine problem dimension
245  if (dIm == -1) {
246  int nb_ents_3d;
247  CHKERR m_field.get_moab().get_number_entities_by_dimension(
248  meshSet, 3, nb_ents_3d, true);
249  if (nb_ents_3d > 0) {
250  dIm = 3;
251  } else {
252  int nb_ents_2d;
253  CHKERR m_field.get_moab().get_number_entities_by_dimension(
254  meshSet, 2, nb_ents_2d, true);
255  if (nb_ents_2d > 0) {
256  dIm = 2;
257  } else {
258  dIm = 1;
259  }
260  }
261  }
262 
263  if (!boundaryMeshset)
265  if (!skeletonMeshset)
267  if (addSkeletonFE)
269 
270  Range ents;
271  CHKERR m_field.get_moab().get_entities_by_dimension(meshSet, dIm, ents, true);
272  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(ents, bitLevel,
273  false);
274 
275  PetscLogEventEnd(MOFEM_EVENT_SimpleLoadMesh, 0, 0, 0, 0);
277 }
278 
281  Interface &m_field = cOre;
282  if (m_field.get_comm_size() == 1)
284  else
285  CHKERR loadFile("PARALLEL=READ_PART;"
286  "PARALLEL_RESOLVE_SHARED_ENTS;"
287  "PARTITION=PARALLEL_PARTITION;",
288  meshFileName);
290 }
291 
293 Simple::addDomainField(const std::string &name, const FieldSpace space,
294  const FieldApproximationBase base,
295  const FieldCoefficientsNumber nb_of_cooficients,
296  const TagType tag_type, const enum MoFEMTypes bh,
297  int verb) {
298 
299  Interface &m_field = cOre;
301  CHKERR m_field.add_field(name, space, base, nb_of_cooficients, tag_type, bh,
302  verb);
303  if (space == NOFIELD)
304  noFieldFields.push_back(name);
305  else
306  domainFields.push_back(name);
308 }
309 
311 Simple::addBoundaryField(const std::string &name, const FieldSpace space,
312  const FieldApproximationBase base,
313  const FieldCoefficientsNumber nb_of_cooficients,
314  const TagType tag_type, const enum MoFEMTypes bh,
315  int verb) {
316  Interface &m_field = cOre;
318  CHKERR m_field.add_field(name, space, base, nb_of_cooficients, tag_type, bh,
319  verb);
320  boundaryFields.push_back(name);
321  if (space == NOFIELD)
322  SETERRQ(
323  PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
324  "NOFIELD space for boundary filed not implemented in Simple interface");
326 }
327 
329 Simple::addSkeletonField(const std::string &name, const FieldSpace space,
330  const FieldApproximationBase base,
331  const FieldCoefficientsNumber nb_of_cooficients,
332  const TagType tag_type, const enum MoFEMTypes bh,
333  int verb) {
334 
335  Interface &m_field = cOre;
337  CHKERR m_field.add_field(name, space, base, nb_of_cooficients, tag_type, bh,
338  verb);
339  skeletonFields.push_back(name);
340  if (space == NOFIELD)
341  SETERRQ(
342  PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
343  "NOFIELD space for boundary filed not implemented in Simple interface");
344 
346 }
347 
349 Simple::addDataField(const std::string &name, const FieldSpace space,
350  const FieldApproximationBase base,
351  const FieldCoefficientsNumber nb_of_cooficients,
352  const TagType tag_type, const enum MoFEMTypes bh,
353  int verb) {
354 
355  Interface &m_field = cOre;
357  CHKERR m_field.add_field(name, space, base, nb_of_cooficients, tag_type, bh,
358  verb);
359  if (space == NOFIELD)
360  noFieldDataFields.push_back(name);
361  else
362  dataFields.push_back(name);
364 }
365 
366 MoFEMErrorCode Simple::removeDomainField(const std::string &name) {
367  Interface &m_field = cOre;
369 
370  auto remove_field_from_list = [&](auto &vec) {
371  auto it = std::find(vec.begin(), vec.end(), name);
372  if (it != vec.end())
373  vec.erase(it);
374  };
375 
376  remove_field_from_list(noFieldFields);
377  remove_field_from_list(domainFields);
378 
380 }
381 
382 MoFEMErrorCode Simple::removeBoundaryField(const std::string &name) {
383  Interface &m_field = cOre;
385 
386  auto remove_field_from_list = [&](auto &vec) {
387  auto it = std::find(vec.begin(), vec.end(), name);
388  if (it != vec.end())
389  vec.erase(it);
390  };
391 
392  remove_field_from_list(boundaryFields);
393 
395 }
396 
397 MoFEMErrorCode Simple::removeSkeletonField(const std::string &name) {
398  Interface &m_field = cOre;
400 
401  auto remove_field_from_list = [&](auto &vec) {
402  auto it = std::find(vec.begin(), vec.end(), name);
403  if (it != vec.end())
404  vec.erase(it);
405  };
406 
407  remove_field_from_list(skeletonFields);
408 
410 }
411 
413  Interface &m_field = cOre;
415 
416  auto clear_rows_and_cols = [&](auto &fe_name) {
418  auto fe_ptr = m_field.get_finite_elements();
419  auto &fe_by_name = const_cast<FiniteElement_multiIndex *>(fe_ptr)
420  ->get<FiniteElement_name_mi_tag>();
421  auto it_fe = fe_by_name.find(fe_name);
422  if (it_fe != fe_by_name.end()) {
423 
424  if (!fe_by_name.modify(it_fe, FiniteElement_row_change_bit_reset()))
425  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
426  "modification unsuccessful");
427 
428  if (!fe_by_name.modify(it_fe, FiniteElement_col_change_bit_reset()))
429  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
430  "modification unsuccessful");
431  }
433  };
434  CHKERR clear_rows_and_cols(domainFE);
435  CHKERR clear_rows_and_cols(boundaryFE);
436  CHKERR clear_rows_and_cols(skeletonFE);
437 
438  // Define finite elements
440 
441  auto add_fields = [&](auto &fe_name, auto &fields) {
443  for (auto &field : fields) {
444  CHKERR m_field.modify_finite_element_add_field_row(fe_name, field);
445  CHKERR m_field.modify_finite_element_add_field_col(fe_name, field);
446  CHKERR m_field.modify_finite_element_add_field_data(fe_name, field);
447  }
449  };
450 
451  auto add_data_fields = [&](auto &fe_name, auto &fields) {
453  for (auto &field : fields)
454  CHKERR m_field.modify_finite_element_add_field_data(fe_name, field);
456  };
457 
458  CHKERR add_fields(domainFE, domainFields);
459  CHKERR add_data_fields(domainFE, dataFields);
460  CHKERR add_fields(domainFE, noFieldFields);
461  CHKERR add_data_fields(domainFE, noFieldDataFields);
462 
463  if (addBoundaryFE || !boundaryFields.empty()) {
465  CHKERR add_fields(boundaryFE, domainFields);
466  if (!boundaryFields.empty())
467  CHKERR add_fields(boundaryFE, boundaryFields);
468  CHKERR add_data_fields(boundaryFE, dataFields);
469  CHKERR add_data_fields(boundaryFE, noFieldDataFields);
470  CHKERR add_fields(boundaryFE, noFieldFields);
471  }
472  if (addSkeletonFE || !skeletonFields.empty()) {
474  CHKERR add_fields(skeletonFE, domainFields);
475  if (!skeletonFields.empty())
476  CHKERR add_fields(skeletonFE, skeletonFields);
477  CHKERR add_data_fields(skeletonFE, dataFields);
478  CHKERR add_data_fields(skeletonFE, noFieldDataFields);
479  CHKERR add_fields(skeletonFE, noFieldFields);
480  }
482 }
483 
484 MoFEMErrorCode Simple::defineProblem(const PetscBool is_partitioned) {
485  Interface &m_field = cOre;
487  // Create dm instance
488  dM = createSmartDM(m_field.get_comm(), "DMMOFEM");
489  // set dm data structure which created mofem data structures
490  CHKERR DMMoFEMCreateMoFEM(dM, &m_field, nameOfProblem.c_str(), bitLevel,
491  bitLevelMask);
492  CHKERR DMSetFromOptions(dM);
494  if (addBoundaryFE || !boundaryFields.empty()) {
496  }
497  if (addSkeletonFE || !skeletonFields.empty()) {
499  }
500  for (std::vector<std::string>::iterator fit = otherFEs.begin();
501  fit != otherFEs.end(); ++fit) {
502  CHKERR DMMoFEMAddElement(dM, fit->c_str());
503  }
504  CHKERR DMMoFEMSetIsPartitioned(dM, is_partitioned);
506 }
507 
508 MoFEMErrorCode Simple::setFieldOrder(const std::string field_name,
509  const int order, const Range *ents) {
511  fieldsOrder.insert(
512  {field_name, {order, ents == NULL ? Range() : Range(*ents)}});
514 }
515 
517  Interface &m_field = cOre;
519  PetscLogEventBegin(MOFEM_EVENT_SimpleBuildFields, 0, 0, 0, 0);
520 
521  auto comm_interface_ptr = m_field.getInterface<CommInterface>();
522  auto bit_ref_ptr = m_field.getInterface<BitRefManager>();
523 
524  // Add entities to the fields
525  auto add_ents_to_field = [&](auto meshset, auto dim, auto &fields) {
527  for (auto &field : fields) {
528  CHKERR m_field.add_ents_to_field_by_dim(meshset, dim, field);
529  CHKERR comm_interface_ptr->synchroniseFieldEntities(field, 0);
530  }
532  };
533 
534  auto make_no_field_ents = [&](auto &fields) {
536  for (auto &field : fields) {
537  std::array<double, 6> coords = {0, 0, 0, 0, 0, 0};
538  CHKERR m_field.create_vertices_and_add_to_field(field, coords.data(), 1);
539  CHKERR comm_interface_ptr->makeFieldEntitiesMultishared(field, 0);
540  CHKERR bit_ref_ptr->setFieldEntitiesBitRefLevel(field, bitLevel);
541  }
543  };
544 
545  CHKERR add_ents_to_field(meshSet, dIm, domainFields);
546  CHKERR add_ents_to_field(meshSet, dIm, dataFields);
547  CHKERR add_ents_to_field(boundaryMeshset, dIm - 1, boundaryFields);
548  CHKERR add_ents_to_field(skeletonMeshset, dIm - 1, skeletonFields);
549 
550  std::set<std::string> nofield_fields;
551  for (auto &f : noFieldFields)
552  nofield_fields.insert(f);
553  for (auto &f : noFieldDataFields)
554  nofield_fields.insert(f);
555 
556  CHKERR make_no_field_ents(nofield_fields);
557 
558  // Set order
559  auto set_order = [&](auto meshset, auto dim, auto &fields) {
561  for (auto &field : fields) {
562  if (fieldsOrder.find(field) == fieldsOrder.end()) {
563  SETERRQ1(PETSC_COMM_WORLD, MOFEM_INVALID_DATA,
564  "Order for field not set %s", field.c_str());
565  }
566  int dds = 0;
567  const Field *field_ptr = m_field.get_field_structure(field);
568  switch (field_ptr->getSpace()) {
569  case L2:
570  dds = dim;
571  break;
572  case HDIV:
573  dds = 2;
574  break;
575  case HCURL:
576  dds = 1;
577  break;
578  case H1:
579  dds = 1;
580  break;
581  default:
582  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
583  "Glasgow we have a problem");
584  }
585 
586  auto set_order = [&](auto field, auto &ents) {
588 
589  auto range = fieldsOrder.equal_range(field);
590  for (auto o = range.first; o != range.second; ++o) {
591  if (!o->second.second.empty())
592  ents = intersect(ents, o->second.second);
593  CHKERR m_field.set_field_order(ents, field, o->second.first);
594  }
595 
597  };
598 
599  if (field_ptr->getSpace() == H1) {
600  if (field_ptr->getApproxBase() == AINSWORTH_BERNSTEIN_BEZIER_BASE) {
601  Range ents;
602  CHKERR m_field.get_field_entities_by_dimension(field, 0, ents);
603  CHKERR set_order(field, ents);
604  } else {
605  CHKERR m_field.set_field_order(meshSet, MBVERTEX, field, 1);
606  }
607  }
608  for (int dd = dds; dd <= dim; dd++) {
609  Range ents;
610  CHKERR m_field.get_field_entities_by_dimension(field, dd, ents);
611  CHKERR set_order(field, ents);
612  }
613  }
615  };
616 
617  CHKERR set_order(meshSet, dIm, domainFields);
618  CHKERR set_order(meshSet, dIm, dataFields);
619  CHKERR set_order(meshSet, dIm - 1, boundaryFields);
620  CHKERR set_order(meshSet, dIm - 1, skeletonFields);
621 
622  // Build fields
623  CHKERR m_field.build_fields();
624  PetscLogEventEnd(MOFEM_EVENT_SimpleBuildFields, 0, 0, 0, 0);
626 }
627 
629  Interface &m_field = cOre;
631  PetscLogEventBegin(MOFEM_EVENT_SimpleBuildFiniteElements, 0, 0, 0, 0);
632  // Add finite elements
634  true);
636  if (addBoundaryFE || !boundaryFields.empty()) {
638  boundaryFE, true);
640  }
641  if (addSkeletonFE || !skeletonFields.empty()) {
643  skeletonFE, true);
645  }
646  for (std::vector<std::string>::iterator fit = otherFEs.begin();
647  fit != otherFEs.end(); ++fit) {
648  CHKERR m_field.build_finite_elements(*fit);
649  }
650  PetscLogEventEnd(MOFEM_EVENT_SimpleBuildFiniteElements, 0, 0, 0, 0);
652 }
653 
655  Interface &m_field = cOre;
657  PetscLogEventBegin(MOFEM_EVENT_SimpleBuildProblem, 0, 0, 0, 0);
659  // Set problem by the DOFs on the fields rather that by adding DOFs on the
660  // elements
661  m_field.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
663  m_field.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_FALSE;
664  PetscLogEventEnd(MOFEM_EVENT_SimpleBuildProblem, 0, 0, 0, 0);
666 }
667 
668 MoFEMErrorCode Simple::setUp(const PetscBool is_partitioned) {
669  Interface &m_field = cOre;
671 
672  PetscLogEventBegin(MOFEM_EVENT_SimpleSetUP, 0, 0, 0, 0);
673 
675  if (addSkeletonFE || !skeletonFields.empty())
677  CHKERR defineProblem(is_partitioned);
681 
682  PetscLogEventEnd(MOFEM_EVENT_SimpleSetUP, 0, 0, 0, 0);
684 }
685 
687  Interface &m_field = cOre;
689 
690  const Problem *problem_ptr;
691  CHKERR DMMoFEMGetProblemPtr(dM, &problem_ptr);
692  const auto problem_name = problem_ptr->getName();
693  CHKERR m_field.modify_problem_ref_level_set_bit(problem_name, bitLevel);
694  CHKERR m_field.modify_problem_mask_ref_level_set_bit(problem_name,
695  bitLevelMask);
696 
698  if (addSkeletonFE || !skeletonFields.empty())
702 
703  PetscLogEventBegin(MOFEM_EVENT_SimpleBuildProblem, 0, 0, 0, 0);
705  // Set problem by the DOFs on the fields rather that by adding DOFs on the
706  // elements
707  m_field.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
709  m_field.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_FALSE;
710  PetscLogEventEnd(MOFEM_EVENT_SimpleBuildProblem, 0, 0, 0, 0);
711 
713 }
714 
717  CHKERR PetscObjectReference(getPetscObject(dM.get()));
718  *dm = dM.get();
720 }
721 
722 /**
723  * @brief Delete dm
724  *
725  * @return MoFEMErrorCode
726  */
729  dM.reset();
731 }
732 
733 /**
734  * @brief Delete finite elements
735  *
736  * @return MoFEMErrorCode
737  */
739  Interface &m_field = cOre;
741  for (auto fe : {domainFE, boundaryFE, skeletonFE}) {
742  if (m_field.check_finite_element(fe)) {
743  CHKERR m_field.delete_finite_element(fe);
744  }
745  }
747 }
748 
750 Simple::addFieldToEmptyFieldBlocks(const std::string row_field,
751  const std::string col_field) const {
752  Interface &m_field = cOre;
755  getProblemName(), row_field, col_field);
757 }
758 
760  Interface &m_field = cOre;
762  ParallelComm *pcomm =
763  ParallelComm::get_pcomm(&m_field.get_moab(), MYPCOMM_INDEX);
764 
765  auto get_skin = [&](auto meshset) {
766  // filter not owned entities, those are not on boundary
767 
768  Range domain_ents;
769  CHKERR m_field.get_moab().get_entities_by_dimension(meshset, dIm,
770  domain_ents, true);
771  CHKERR pcomm->filter_pstatus(domain_ents,
772  PSTATUS_SHARED | PSTATUS_MULTISHARED,
773  PSTATUS_NOT, -1, nullptr);
774 
775  Skinner skin(&m_field.get_moab());
776  Range domain_skin;
777  CHKERR skin.find_skin(0, domain_ents, false, domain_skin);
778  CHKERR pcomm->filter_pstatus(domain_skin,
779  PSTATUS_SHARED | PSTATUS_MULTISHARED,
780  PSTATUS_NOT, -1, nullptr);
781  return domain_skin;
782  };
783 
784  auto create_boundary_meshset = [&](auto &&domain_skin) {
786  // create boundary meshset
787  if (boundaryMeshset != 0) {
789  }
790  CHKERR m_field.get_moab().create_meshset(MESHSET_SET, boundaryMeshset);
791  CHKERR m_field.get_moab().add_entities(boundaryMeshset, domain_skin);
792  for (int dd = 0; dd != dIm - 1; dd++) {
793  Range adj;
794  CHKERR m_field.get_moab().get_adjacencies(domain_skin, dd, false, adj,
795  moab::Interface::UNION);
796  CHKERR m_field.get_moab().add_entities(boundaryMeshset, adj);
797  }
799  };
800 
801  CHKERR create_boundary_meshset(get_skin(meshSet));
802 
804 }
805 
807  Interface &m_field = cOre;
809 
810  ParallelComm *pcomm =
811  ParallelComm::get_pcomm(&m_field.get_moab(), MYPCOMM_INDEX);
812 
813  auto create_skeleton_meshset = [&](auto meshset) {
815  // create boundary meshset
816  if (skeletonMeshset != 0) {
818  }
819  Range boundary_ents, skeleton_ents;
820  CHKERR m_field.get_moab().get_entities_by_dimension(boundaryMeshset,
821  dIm - 1, boundary_ents);
822  Range domain_ents;
823  CHKERR m_field.get_moab().get_entities_by_dimension(meshset, dIm,
824  domain_ents, true);
825  CHKERR m_field.get_moab().get_adjacencies(
826  domain_ents, dIm - 1, false, skeleton_ents, moab::Interface::UNION);
827  skeleton_ents = subtract(skeleton_ents, boundary_ents);
828  CHKERR pcomm->filter_pstatus(skeleton_ents, PSTATUS_NOT_OWNED, PSTATUS_NOT,
829  -1, nullptr);
830  CHKERR m_field.get_moab().create_meshset(MESHSET_SET, skeletonMeshset);
831  CHKERR m_field.get_moab().add_entities(skeletonMeshset, skeleton_ents);
833  };
834 
835  CHKERR create_skeleton_meshset(meshSet);
836 
838 }
839 
841  Interface &m_field = cOre;
843  MOFEM_LOG("WORLD", Sev::verbose) << "Exchange ghost cells";
844 
845  ParallelComm *pcomm =
846  ParallelComm::get_pcomm(&m_field.get_moab(), MYPCOMM_INDEX);
847  if (pcomm == NULL)
848  pcomm = new ParallelComm(&m_field.get_moab(), m_field.get_comm());
849 
850  Range verts;
851  CHKERR m_field.get_moab().get_entities_by_type(0, MBVERTEX, verts);
852  CHKERR pcomm->exchange_ghost_cells(getDim(), getDim() - 1, 1,
853  3 /**get all adjacent ghosted entities */,
854  true, false, meshSet ? &meshSet : nullptr);
855 
856  Range shared;
857  CHKERR m_field.get_moab().get_entities_by_dimension(0, dIm, shared);
858  for (auto d = dIm - 1; d >= 1; --d) {
859  CHKERR m_field.get_moab().get_adjacencies(shared, d, false, shared,
860  moab::Interface::UNION);
861  }
862  CHKERR pcomm->filter_pstatus(shared, PSTATUS_SHARED | PSTATUS_MULTISHARED,
863  PSTATUS_OR, -1, &shared);
864  Tag part_tag = pcomm->part_tag();
865  CHKERR pcomm->exchange_tags(part_tag, shared);
866 
868 }
869 
870 } // namespace MoFEM
static Index< 'o', 3 > o
MoFEMTypes
Those types control how functions respond on arguments, f.e. error handling.
Definition: definitions.h:110
@ MF_ZERO
Definition: definitions.h:111
FieldApproximationBase
approximation base
Definition: definitions.h:71
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:77
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
FieldSpace
approximation spaces
Definition: definitions.h:95
@ L2
field with C-1 continuity
Definition: definitions.h:101
@ NOFIELD
scalar or vector of scalars describe (no true field)
Definition: definitions.h:97
@ H1
continuous field
Definition: definitions.h:98
@ HCURL
field with continuous tangents
Definition: definitions.h:99
@ HDIV
field with continuous normal traction
Definition: definitions.h:100
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:496
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:47
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
@ MOFEM_INVALID_DATA
Definition: definitions.h:49
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:45
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:574
const int dim
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMMoFEM.cpp:1061
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition: DMMMoFEM.cpp:130
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:461
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMMoFEM.cpp:384
PetscErrorCode DMSetUp_MoFEM(DM dm)
Definition: DMMMoFEM.cpp:1203
multi_index_container< boost::shared_ptr< SideNumber >, indexed_by< ordered_unique< member< SideNumber, EntityHandle, &SideNumber::ent > >, ordered_non_unique< composite_key< SideNumber, const_mem_fun< SideNumber, EntityType, &SideNumber::getEntType >, member< SideNumber, char, &SideNumber::side_number > > > > > SideNumber_multiIndex
SideNumber_multiIndex for SideNumber.
multi_index_container< boost::shared_ptr< FiniteElement >, indexed_by< hashed_unique< tag< FiniteElement_Meshset_mi_tag >, member< FiniteElement, EntityHandle, &FiniteElement::meshset > >, hashed_unique< tag< BitFEId_mi_tag >, const_mem_fun< FiniteElement, BitFEId, &FiniteElement::getId >, HashBit< BitFEId >, EqBit< BitFEId > >, ordered_unique< tag< FiniteElement_name_mi_tag >, const_mem_fun< FiniteElement, boost::string_ref, &FiniteElement::getNameRef > > > > FiniteElement_multiIndex
MultiIndex for entities for FiniteElement.
virtual const FiniteElement_multiIndex * get_finite_elements() const =0
Get the finite elements object.
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string &name_row)=0
set field row which finite element use
virtual MoFEMErrorCode add_ents_to_finite_element_by_dim(const EntityHandle entities, const int dim, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode modify_finite_element_adjacency_table(const std::string &fe_name, const EntityType type, ElementAdjacencyFunct function)=0
modify finite element table, only for advanced user
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string &name_row)=0
set field col which finite element use
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode add_ents_to_field_by_dim(const Range &ents, const int dim, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode get_field_entities_by_dimension(const std::string name, int dim, Range &ents) const =0
get entities in the field by dimension
virtual Field * get_field_structure(const std::string &name)=0
get field structure
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:311
virtual MoFEMErrorCode modify_problem_mask_ref_level_set_bit(const std::string &name_problem, const BitRefLevel &bit)=0
set dof mask ref level for problem
virtual MoFEMErrorCode modify_problem_ref_level_set_bit(const std::string &name_problem, const BitRefLevel &bit)=0
set ref level for problem
char mesh_file_name[255]
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
int FieldCoefficientsNumber
Number of field coefficients.
Definition: Types.hpp:38
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
auto createSmartDM
Creates smart DM object.
PetscObject getPetscObject(T obj)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
Managing BitRefLevels.
Managing BitRefLevels.
virtual int get_comm_size() const =0
virtual bool check_finite_element(const std::string &name) const =0
Check if finite element is in database.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MoFEMErrorCode delete_finite_element(const std::string name, int verb=DEFAULT_VERBOSITY)=0
delete finite element from mofem database
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual MoFEMErrorCode rebuild_database(int verb=DEFAULT_VERBOSITY)=0
Clear database and initialize it once again.
virtual MoFEMErrorCode create_vertices_and_add_to_field(const std::string name, const double coords[], int size, int verb=DEFAULT_VERBOSITY)=0
Create a vertices and add to field object.
Core (interface) class.
Definition: Core.hpp:92
static MoFEMErrorCode defaultEdge(Interface &moab, const Field &field, const EntFiniteElement &fe, std::vector< EntityHandle > &adjacency)
static MoFEMErrorCode defaultFace(Interface &moab, const Field &field, const EntFiniteElement &fe, std::vector< EntityHandle > &adjacency)
Deprecated interface functions.
Finite element data for entity.
Provide data structure for (tensor) field approximation.
FieldApproximationBase getApproxBase() const
Get approximation base.
FieldSpace getSpace() const
Get field approximation space.
keeps basic data about problem
Problem manager is used to build and partition problems.
keeps information about side number for the finite element
Simple interface for fast problem set-up.
Definition: Simple.hpp:33
EntityHandle boundaryMeshset
meshset with boundary
Definition: Simple.hpp:444
MoFEMErrorCode buildProblem()
Build problem.
Definition: Simple.cpp:654
std::vector< std::string > dataFields
Data fields.
Definition: Simple.hpp:453
std::vector< std::string > domainFields
domain fields
Definition: Simple.hpp:450
MoFEMErrorCode removeDomainField(const std::string &name)
Remove field form domain.
Definition: Simple.cpp:366
MoFEMErrorCode defineFiniteElements()
Define finite elements.
Definition: Simple.cpp:412
MoFEM::Core & cOre
Definition: Simple.hpp:431
char meshFileName[255]
mesh file name
Definition: Simple.hpp:467
EntityHandle meshSet
domain meshset
Definition: Simple.hpp:443
int getDim() const
Get the problem dimension.
Definition: Simple.hpp:270
std::vector< std::string > otherFEs
Other finite elements.
Definition: Simple.hpp:465
MoFEMErrorCode setSkeletonAdjacency()
Definition: Simple.cpp:28
MoFEMErrorCode buildFiniteElements()
Build finite elements.
Definition: Simple.cpp:628
MoFEMErrorCode addDataField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
Definition: Simple.cpp:349
MoFEMErrorCode exchangeGhostCells()
Definition: Simple.cpp:840
std::multimap< std::string, std::pair< int, Range > > fieldsOrder
fields order
Definition: Simple.hpp:458
PetscLogEvent MOFEM_EVENT_SimpleBuildFields
Definition: Simple.hpp:438
PetscLogEvent MOFEM_EVENT_SimpleLoadMesh
Definition: Simple.hpp:437
MoFEMErrorCode addDomainField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
Definition: Simple.cpp:293
std::vector< std::string > noFieldFields
NOFIELD field name.
Definition: Simple.hpp:454
bool addSkeletonFE
Add skeleton FE.
Definition: Simple.hpp:447
bool addBoundaryFE
Add boundary FE.
Definition: Simple.hpp:448
MoFEMErrorCode addFieldToEmptyFieldBlocks(const std::string row_field, const std::string col_field) const
add empty block to problem
Definition: Simple.cpp:750
std::string skeletonFE
skeleton finite element
Definition: Simple.hpp:463
std::vector< std::string > noFieldDataFields
NOFIELD field name.
Definition: Simple.hpp:455
std::string boundaryFE
boundary finite element
Definition: Simple.hpp:462
std::string nameOfProblem
problem name
Definition: Simple.hpp:460
SmartPetscObj< DM > getDM()
Return smart DM object.
Definition: Simple.hpp:260
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:216
std::string domainFE
domain finite element
Definition: Simple.hpp:461
std::vector< std::string > boundaryFields
boundary fields
Definition: Simple.hpp:451
MoFEMErrorCode reSetUp()
Rebuild internal MoFEM data structures.
Definition: Simple.cpp:686
SmartPetscObj< DM > dM
Discrete manager (interface to PETSc/MoFEM functions)
Definition: Simple.hpp:471
MoFEMErrorCode buildFields()
Build fields.
Definition: Simple.cpp:516
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:230
BitRefLevel bitLevel
BitRefLevel of the probelm.
Definition: Simple.hpp:433
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:508
MoFEMErrorCode addBoundaryField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on boundary.
Definition: Simple.cpp:311
MoFEMErrorCode defineProblem(const PetscBool is_partitioned=PETSC_TRUE)
define problem
Definition: Simple.cpp:484
std::vector< std::string > skeletonFields
fields on the skeleton
Definition: Simple.hpp:452
MoFEMErrorCode addSkeletonField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on skeleton.
Definition: Simple.cpp:329
MoFEMErrorCode removeSkeletonField(const std::string &name)
Remove field form skeleton.
Definition: Simple.cpp:397
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
Definition: Simple.cpp:22
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:668
MoFEMErrorCode removeBoundaryField(const std::string &name)
Remove field form boundary.
Definition: Simple.cpp:382
const std::string getProblemName() const
Get the Problem Name.
Definition: Simple.hpp:340
int dIm
dimension of problem
Definition: Simple.hpp:468
MoFEMErrorCode deleteDM()
Delete dm.
Definition: Simple.cpp:727
PetscLogEvent MOFEM_EVENT_SimpleSetUP
Definition: Simple.hpp:436
PetscLogEvent MOFEM_EVENT_SimpleKSPSolve
Definition: Simple.hpp:441
MoFEMErrorCode createBoundaryMeshset()
Definition: Simple.cpp:759
PetscLogEvent MOFEM_EVENT_SimpleBuildFiniteElements
Definition: Simple.hpp:439
BitRefLevel bitLevelMask
BitRefLevel of the probelm.
Definition: Simple.hpp:434
PetscLogEvent MOFEM_EVENT_SimpleBuildProblem
Definition: Simple.hpp:440
MoFEMErrorCode deleteFiniteElements()
Delete finite elements.
Definition: Simple.cpp:738
Simple(const MoFEM::Core &core)
Definition: Simple.cpp:200
EntityHandle skeletonMeshset
skeleton meshset with boundary
Definition: Simple.hpp:445
MoFEMErrorCode createSkeletonMeshset()
Definition: Simple.cpp:806
base class for all interface classes
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.