v0.11.1
continuity_check_on_contact_prism_side_ele.cpp
Go to the documentation of this file.
1 /** \file continuity_check_on_contact_prism_side_ele.cpp
2  * \example continuity_check_on_contact_prism_side_ele.cpp
3  *
4  * \brief Testing integration on volume side on contact element
5  *
6  * Checking continuity of H1 ( later hdiv and hcurl spaces on faces), and
7  * testing methods for integration on volume side on contact element.
8  *
9  */
10 
11 /* This file is part of MoFEM.
12  * MoFEM is free software: you can redistribute it and/or modify it under
13  * the terms of the GNU Lesser General Public License as published by the
14  * Free Software Foundation, either version 3 of the License, or (at your
15  * option) any later version.
16  *
17  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
18  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
19  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
20  * License for more details.
21  *
22  * You should have received a copy of the GNU Lesser General Public
23  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
24 
25 #include <MoFEM.hpp>
26 using namespace MoFEM;
27 
28 static char help[] = "...\n\n";
29 
30 int main(int argc, char *argv[]) {
31 
32  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
33 
34  try {
35 
36  moab::Core mb_instance;
37  moab::Interface &moab = mb_instance;
38  int rank;
39  MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
40 
41  PetscBool flg = PETSC_TRUE;
42  char mesh_file_name[255];
43 #if PETSC_VERSION_GE(3, 6, 4)
44  CHKERR PetscOptionsGetString(PETSC_NULL, "", "-my_file", mesh_file_name,
45  255, &flg);
46 #else
47  CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-my_file",
48  mesh_file_name, 255, &flg);
49 #endif
50  if (flg != PETSC_TRUE) {
51  SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
52  }
53 
54  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
55  if (pcomm == NULL)
56  pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
57 
58  const char *option;
59  option = "";
60  CHKERR moab.load_file(mesh_file_name, 0, option);
61 
62  // Create MoFEM database
63  MoFEM::Core core(moab);
64  // Access to database through interface
65  MoFEM::Interface &m_field = core;
66 
67  enum side_contact { MASTERSIDE, SLAVESIDE };
68 
69  // Fields
70  auto get_base = []() -> FieldApproximationBase {
71  enum bases { AINSWORTH, DEMKOWICZ, LASTBASEOP };
72  const char *list_bases[] = {"ainsworth", "demkowicz"};
73  PetscBool flg;
74  PetscInt choice_base_value = AINSWORTH;
75  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
76  LASTBASEOP, &choice_base_value, &flg);
77  if (flg == PETSC_TRUE) {
79  if (choice_base_value == AINSWORTH)
81  else if (choice_base_value == DEMKOWICZ)
82  base = DEMKOWICZ_JACOBI_BASE;
83  return base;
84  }
85  return LASTBASE;
86  };
87 
88  auto get_space = []() -> FieldSpace {
89  enum spaces { H1, HDIV, HCURL, LASTSPACEOP };
90  const char *list_spaces[] = {"h1", "hdiv", "hcurl"};
91  PetscBool flg;
92  PetscInt choice_space_value = HDIV;
93  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-space", list_spaces,
94  LASTSPACEOP, &choice_space_value, &flg);
95  if (flg == PETSC_TRUE) {
97  if (choice_space_value == HDIV)
98  space = FieldSpace::HDIV;
99  else if (choice_space_value == HCURL)
100  space = FieldSpace::HCURL;
101  else if (choice_space_value == H1)
102  space = FieldSpace::H1;
103  return space;
104  }
105  return LASTSPACE;
106  };
107 
108  CHKERR m_field.add_field("F2", get_space(), get_base(), 1);
109 
110  auto add_prism_interface = [&](Range &tets, Range &prisms,
111  Range &master_tris, Range &slave_tris,
112  EntityHandle &meshset_tets,
113  EntityHandle &meshset_prisms,
114  std::vector<BitRefLevel> &bit_levels) {
116  PrismInterface *interface;
117  CHKERR m_field.getInterface(interface);
118 
119  int ll = 1;
120  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, SIDESET, cit)) {
121  if (cit->getName().compare(0, 11, "INT_CONTACT") == 0) {
122  CHKERR PetscPrintf(PETSC_COMM_WORLD, "Insert %s (id: %d)\n",
123  cit->getName().c_str(), cit->getMeshsetId());
124  EntityHandle cubit_meshset = cit->getMeshset();
125  Range tris;
126  CHKERR moab.get_entities_by_type(cubit_meshset, MBTRI, tris, true);
127  master_tris.merge(tris);
128 
129  {
130  // get tet entities from back bit_level
131  EntityHandle ref_level_meshset = 0;
132  CHKERR moab.create_meshset(MESHSET_SET, ref_level_meshset);
134  ->getEntitiesByTypeAndRefLevel(bit_levels.back(),
135  BitRefLevel().set(), MBTET,
136  ref_level_meshset);
138  ->getEntitiesByTypeAndRefLevel(bit_levels.back(),
139  BitRefLevel().set(), MBPRISM,
140  ref_level_meshset);
141  Range ref_level_tets;
142  CHKERR moab.get_entities_by_handle(ref_level_meshset,
143  ref_level_tets, true);
144  // get faces and tets to split
145  CHKERR interface->getSides(cubit_meshset, bit_levels.back(), true,
146  0);
147  // set new bit level
148  bit_levels.push_back(BitRefLevel().set(ll++));
149  // split faces and tets
150  CHKERR interface->splitSides(ref_level_meshset, bit_levels.back(),
151  cubit_meshset, true, true, 0);
152  // clean meshsets
153  CHKERR moab.delete_entities(&ref_level_meshset, 1);
154  }
155  // update cubit meshsets
156  for (_IT_CUBITMESHSETS_FOR_LOOP_(m_field, ciit)) {
157  EntityHandle cubit_meshset = ciit->meshset;
159  ->updateMeshsetByEntitiesChildren(
160  cubit_meshset, bit_levels.back(), cubit_meshset, MBVERTEX,
161  true);
163  ->updateMeshsetByEntitiesChildren(cubit_meshset,
164  bit_levels.back(),
165  cubit_meshset, MBEDGE, true);
167  ->updateMeshsetByEntitiesChildren(cubit_meshset,
168  bit_levels.back(),
169  cubit_meshset, MBTRI, true);
171  ->updateMeshsetByEntitiesChildren(cubit_meshset,
172  bit_levels.back(),
173  cubit_meshset, MBTET, true);
174  }
175  }
176  }
177 
178  for (unsigned int ll = 0; ll != bit_levels.size() - 1; ll++) {
179  CHKERR m_field.delete_ents_by_bit_ref(bit_levels[ll], bit_levels[ll],
180  true);
181  }
182  CHKERR m_field.getInterface<BitRefManager>()->shiftRightBitRef(
183  bit_levels.size() - 1);
184 
185  CHKERR moab.create_meshset(MESHSET_SET, meshset_tets);
186  CHKERR moab.create_meshset(MESHSET_SET, meshset_prisms);
187 
189  ->getEntitiesByTypeAndRefLevel(bit_levels[0], BitRefLevel().set(),
190  MBTET, meshset_tets);
191  CHKERR moab.get_entities_by_handle(meshset_tets, tets, true);
192 
194  ->getEntitiesByTypeAndRefLevel(bit_levels[0], BitRefLevel().set(),
195  MBPRISM, meshset_prisms);
196  CHKERR moab.get_entities_by_handle(meshset_prisms, prisms);
197 
198  Range tris;
199  CHKERR moab.get_adjacencies(prisms, 2, false, tris,
200  moab::Interface::UNION);
201  tris = tris.subset_by_type(MBTRI);
202  slave_tris = subtract(tris, master_tris);
203 
205  };
206 
207  Range all_tets, contact_prisms, master_tris, slave_tris;
208  EntityHandle meshset_tets, meshset_prisms;
209  std::vector<BitRefLevel> bit_levels;
210 
211  bit_levels.push_back(BitRefLevel().set(0));
212  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
213  0, 3, bit_levels[0]);
214 
215  CHKERR add_prism_interface(all_tets, contact_prisms, master_tris,
216  slave_tris, meshset_tets, meshset_prisms,
217  bit_levels);
218 
219  // meshset consisting all entities in mesh
220  EntityHandle root_set = moab.get_root_set();
221  // add entities to field
222  CHKERR m_field.add_ents_to_field_by_type(root_set, MBTET, "F2");
223  CHKERR m_field.add_ents_to_field_by_type(root_set, MBPRISM, "F2");
224 
225  // set app. order
226  int order = 2;
227  CHKERR m_field.set_field_order(root_set, MBTET, "F2", order);
228  CHKERR m_field.set_field_order(root_set, MBTRI, "F2", order);
229  CHKERR m_field.set_field_order(root_set, MBEDGE, "F2", order);
230 
231  if (get_space() == FieldSpace::H1) {
232  CHKERR m_field.set_field_order(root_set, MBVERTEX, "F2", 1);
233  }
234 
235  CHKERR m_field.build_fields();
236 
237  // add elements
238  CHKERR m_field.add_finite_element("V1");
239  CHKERR m_field.add_finite_element("C2");
240  CHKERR m_field.modify_finite_element_add_field_row("V1", "F2");
241  CHKERR m_field.modify_finite_element_add_field_col("V1", "F2");
242  CHKERR m_field.modify_finite_element_add_field_data("V1", "F2");
243  CHKERR m_field.modify_finite_element_add_field_row("C2", "F2");
244  CHKERR m_field.modify_finite_element_add_field_col("C2", "F2");
245  CHKERR m_field.modify_finite_element_add_field_data("C2", "F2");
246 
247  CHKERR m_field.add_ents_to_finite_element_by_type(root_set, MBTET, "V1");
248 
249  Range prism;
250  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
251  bit_levels[0], BitRefLevel().set(), MBPRISM, prism);
252  CHKERR m_field.add_ents_to_finite_element_by_type(root_set, MBPRISM, "C2");
253 
254  CHKERR m_field.build_finite_elements();
255  CHKERR m_field.build_adjacencies(bit_levels[0]);
256 
257  // Problems
258  CHKERR m_field.add_problem("P1");
259 
260  // set refinement level for problem
261  CHKERR m_field.modify_problem_ref_level_add_bit("P1", bit_levels[0]);
262  CHKERR m_field.modify_problem_add_finite_element("P1", "V1");
263  CHKERR m_field.modify_problem_add_finite_element("P1", "C2");
264 
265  // build problems
266  ProblemsManager *prb_mng_ptr;
267  CHKERR m_field.getInterface(prb_mng_ptr);
268  CHKERR prb_mng_ptr->buildProblem("P1", true);
269 
270  struct CommonData {
271  MatrixDouble dotNormalFace;
272  MatrixDouble dotNormalEleLeft;
273  MatrixDouble dotNormalEleRight;
274  MatrixDouble shapeFunH1Values;
275  MatrixDouble shapeFunH1VolSide;
276  };
277 
278  struct OnContactSideMaster
280 
282  struct OpVolSide : public VolumeElementForcesAndSourcesCoreOnContactPrismSide::
284 
285  CommonData &elemData;
286  OpVolSide(CommonData &elem_data)
288  "F2", UserDataOperator::OPROW),
289  elemData(elem_data) {}
290  MoFEMErrorCode doWork(int side, EntityType type,
293 
294  if (type == MBTRI && side == getFaceSideNumber()) {
295 
296  MatrixDouble diff =
297  getCoordsAtGaussPts() - getMasterCoordsAtGaussPts();
298  constexpr double eps = 1e-12;
299  if (norm_inf(diff) > eps)
300  SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
301  "coordinates at integration pts are different");
302 
303  const size_t nb_dofs = data.getN().size2() / 3;
304  const size_t nb_integration_pts = data.getN().size1();
305 
306  if (data.getSpace() == H1) {
307  MatrixDouble *ptr_elem_data = nullptr;
308  ptr_elem_data = &elemData.shapeFunH1VolSide;
309 
310  MatrixDouble &elem_data = *ptr_elem_data;
311  elem_data.resize(nb_integration_pts, nb_dofs, false);
312  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
313  FTensor::Tensor0<double *> t_base(&data.getN()(gg, 0));
314  for (size_t bb = 0; bb != nb_dofs; ++bb) {
315  elem_data(gg, bb) = t_base;
316  ++t_base;
317  }
318  }
319  }
320  }
322  }
323  };
324 
325  CommonData &elemData;
326  OnContactSideMaster(MoFEM::Interface &m_field, CommonData &elem_data)
328  "F2", UserDataOperator::OPROW,
330  FACEMASTER),
331  volSideFe(m_field), elemData(elem_data) {
332  volSideFe.getOpPtrVector().push_back(
333  new OnContactSideMaster::OpVolSide(elemData));
334  }
335 
336  MoFEMErrorCode doWork(int side, EntityType type,
338 
340  if (type == MBTRI && side == 0) {
341  const size_t nb_dofs = data.getN().size2() / 3;
342  const size_t nb_integration_pts = data.getN().size1();
343  if (data.getSpace() == H1) {
344  elemData.shapeFunH1Values.resize(nb_integration_pts, nb_dofs,
345  false);
346  elemData.shapeFunH1VolSide.resize(nb_integration_pts, nb_dofs,
347  false);
348 
349  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
350  FTensor::Tensor0<double *> t_base(&data.getN()(gg, 0));
351  for (size_t bb = 0; bb != nb_dofs; ++bb) {
352  elemData.shapeFunH1Values(gg, bb) = t_base;
353  ++t_base;
354  }
355  }
356  }
357  std::string side_fe_name = "V1";
358  const EntityHandle tri_master = getSideEntity(3, MBTRI);
359  CHKERR loopSideVolumes(side_fe_name, volSideFe, 3, tri_master);
360 
361  auto check_continuity_of_base = [&](auto &vol_dot_data) {
363 
364  if (vol_dot_data.size1() != elemData.dotNormalFace.size1())
365  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
366  "Inconsistent number of integration points");
367 
368  if (vol_dot_data.size2() != elemData.dotNormalFace.size2())
369  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
370  "Inconsistent number of base functions");
371  constexpr double eps = 1e-12;
372  for (size_t gg = 0; gg != vol_dot_data.size1(); ++gg)
373  for (size_t bb = 0; bb != vol_dot_data.size2(); ++bb) {
374  const double error = std::abs(vol_dot_data(gg, bb) -
375  elemData.dotNormalFace(gg, bb));
376  if (error > eps)
377  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
378  "Inconsistency %3.4e != %3.4e", vol_dot_data(gg, bb),
379  elemData.dotNormalFace(gg, bb));
380  }
382  };
383 
384  auto check_continuity_of_h1_base = [&](auto &vol_data) {
386 
387  if (vol_data.size1() != elemData.shapeFunH1Values.size1())
388  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
389  "Inconsistent number of integration points");
390 
391  if (vol_data.size2() != elemData.shapeFunH1Values.size2())
392  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
393  "Inconsistent number of base functions");
394  constexpr double eps = 1e-12;
395  for (size_t gg = 0; gg != vol_data.size1(); ++gg)
396  for (size_t bb = 0; bb != vol_data.size2(); ++bb) {
397  const double error = std::abs(
398  vol_data(gg, bb) - elemData.shapeFunH1Values(gg, bb));
399  if (error > eps)
400  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
401  "Inconsistency %3.4e != %3.4e", vol_data(gg, bb),
402  elemData.shapeFunH1Values(gg, bb));
403  }
405  };
406 
407  if (elemData.dotNormalEleLeft.size2() != 0)
408  CHKERR check_continuity_of_base(elemData.dotNormalEleLeft);
409  else if (elemData.dotNormalEleRight.size2() != 0)
410  CHKERR check_continuity_of_base(elemData.dotNormalEleRight);
411  else if (elemData.shapeFunH1VolSide.size2() != 0)
412  CHKERR check_continuity_of_h1_base(elemData.shapeFunH1VolSide);
413  }
415  }
416  };
417 
418  struct OnContactSideSlave
420 
422  struct OpVolSide : public VolumeElementForcesAndSourcesCoreOnContactPrismSide::
424 
425  CommonData &elemData;
426  OpVolSide(CommonData &elem_data)
428  "F2", UserDataOperator::OPROW),
429  elemData(elem_data) {}
430  MoFEMErrorCode doWork(int side, EntityType type,
433 
434  if (type == MBTRI && side == getFaceSideNumber()) {
435 
436  MatrixDouble diff =
437  getCoordsAtGaussPts() - getMasterCoordsAtGaussPts();
438  constexpr double eps = 1e-12;
439  if (norm_inf(diff) > eps)
440  SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
441  "coordinates at integration pts are different");
442 
443  const size_t nb_dofs = data.getN().size2() / 3;
444  const size_t nb_integration_pts = data.getN().size1();
445 
446  if (data.getSpace() == H1) {
447  MatrixDouble *ptr_elem_data = nullptr;
448  ptr_elem_data = &elemData.shapeFunH1VolSide;
449 
450  MatrixDouble &elem_data = *ptr_elem_data;
451  elem_data.resize(nb_integration_pts, nb_dofs, false);
452  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
453  FTensor::Tensor0<double *> t_base(&data.getN()(gg, 0));
454  for (size_t bb = 0; bb != nb_dofs; ++bb) {
455  elem_data(gg, bb) = t_base;
456  ++t_base;
457  }
458  }
459  }
460  }
462  }
463  };
464 
465  CommonData &elemData;
466  OnContactSideSlave(MoFEM::Interface &m_field, CommonData &elem_data)
468  "F2", UserDataOperator::OPROW,
470  FACESLAVE),
471  volSideFe(m_field), elemData(elem_data) {
472  volSideFe.getOpPtrVector().push_back(
473  new OnContactSideSlave::OpVolSide(elemData));
474  }
475 
476  MoFEMErrorCode doWork(int side, EntityType type,
478 
480  if (type == MBTRI && side == 0) {
481  const size_t nb_dofs = data.getN().size2() / 3;
482  const size_t nb_integration_pts = data.getN().size1();
483  if (data.getSpace() == H1) {
484  elemData.shapeFunH1Values.resize(nb_integration_pts, nb_dofs,
485  false);
486  elemData.shapeFunH1VolSide.resize(nb_integration_pts, nb_dofs,
487  false);
488 
489  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
490  FTensor::Tensor0<double *> t_base(&data.getN()(gg, 0));
491  for (size_t bb = 0; bb != nb_dofs; ++bb) {
492  elemData.shapeFunH1Values(gg, bb) = t_base;
493  ++t_base;
494  }
495  }
496  }
497 
498  std::string side_fe_name = "V1";
499  const EntityHandle tri_slave = getSideEntity(4, MBTRI);
500  CHKERR loopSideVolumes(side_fe_name, volSideFe, 3, tri_slave);
501 
502  auto check_continuity_of_base = [&](auto &vol_dot_data) {
504 
505  if (vol_dot_data.size1() != elemData.dotNormalFace.size1())
506  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
507  "Inconsistent number of integration points");
508 
509  if (vol_dot_data.size2() != elemData.dotNormalFace.size2())
510  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
511  "Inconsistent number of base functions");
512  constexpr double eps = 1e-12;
513  for (size_t gg = 0; gg != vol_dot_data.size1(); ++gg)
514  for (size_t bb = 0; bb != vol_dot_data.size2(); ++bb) {
515  const double error = std::abs(vol_dot_data(gg, bb) -
516  elemData.dotNormalFace(gg, bb));
517  if (error > eps)
518  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
519  "Inconsistency %3.4e != %3.4e", vol_dot_data(gg, bb),
520  elemData.dotNormalFace(gg, bb));
521  }
523  };
524 
525  auto check_continuity_of_h1_base = [&](auto &vol_data) {
527 
528  if (vol_data.size1() != elemData.shapeFunH1Values.size1())
529  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
530  "Inconsistent number of integration points");
531 
532  if (vol_data.size2() != elemData.shapeFunH1Values.size2())
533  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
534  "Inconsistent number of base functions");
535  constexpr double eps = 1e-12;
536  for (size_t gg = 0; gg != vol_data.size1(); ++gg)
537  for (size_t bb = 0; bb != vol_data.size2(); ++bb) {
538  const double error = std::abs(
539  vol_data(gg, bb) - elemData.shapeFunH1Values(gg, bb));
540  if (error > eps)
541  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
542  "Inconsistency %3.4e != %3.4e", vol_data(gg, bb),
543  elemData.shapeFunH1Values(gg, bb));
544  }
546  };
547 
548  if (elemData.dotNormalEleLeft.size2() != 0)
549  CHKERR check_continuity_of_base(elemData.dotNormalEleLeft);
550  else if (elemData.dotNormalEleRight.size2() != 0)
551  CHKERR check_continuity_of_base(elemData.dotNormalEleRight);
552  else if (elemData.shapeFunH1VolSide.size2() != 0)
553  CHKERR check_continuity_of_h1_base(elemData.shapeFunH1VolSide);
554  }
556  }
557  };
558 
559  CommonData elem_data;
560 
561  ContactPrismElementForcesAndSourcesCore contact_prism_fe_master(m_field);
562  ContactPrismElementForcesAndSourcesCore contact_prism_fe_slave(m_field);
563 
564  // OnContactSideMaster
565  contact_prism_fe_master.getOpPtrVector().push_back(
566  new OnContactSideMaster(m_field, elem_data));
567  // OnContactSideSlave
568  contact_prism_fe_slave.getOpPtrVector().push_back(
569  new OnContactSideSlave(m_field, elem_data));
570 
571  CHKERR m_field.loop_finite_elements("P1", "C2", contact_prism_fe_master);
572  CHKERR m_field.loop_finite_elements("P1", "C2", contact_prism_fe_slave);
573  }
574  CATCH_ERRORS;
575 
577 
578  return 0;
579 }
int main(int argc, char *argv[])
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:441
FieldApproximationBase
approximation base
Definition: definitions.h:150
@ LASTBASE
Definition: definitions.h:161
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:158
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
FieldSpace
approximation spaces
Definition: definitions.h:174
@ LASTSPACE
FieldSpace in [ 0, LASTSPACE )
Definition: definitions.h:181
@ H1
continuous field
Definition: definitions.h:177
@ HCURL
field with continuous tangents
Definition: definitions.h:178
@ HDIV
field with continuous normal traction
Definition: definitions.h:179
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:292
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
@ SIDESET
Definition: definitions.h:216
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:132
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
#define CHKERR
Inline error check.
Definition: definitions.h:604
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
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_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 add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
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 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 add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets in a moFEM field.
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
MoFEMErrorCode buildProblem(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures
virtual MoFEMErrorCode modify_problem_add_finite_element(const std::string &name_problem, const std::string &fe_name)=0
add finite element to problem, this add entities assigned to finite element to a particular problem
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
char mesh_file_name[255]
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
CoreTmp< 0 > Core
Definition: Core.hpp:1140
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1953
Managing BitRefLevels.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MoFEMErrorCode delete_ents_by_bit_ref(const BitRefLevel bit, const BitRefLevel mask, const bool remove_parent=false, int verb=DEFAULT_VERBOSITY)=0
delete entities form mofem and moab 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.
Core (interface) class.
Definition: Core.hpp:77
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:60
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:100
Data on single entity (This is passed as argument to DataOperator::doWork)
FieldSpace & getSpace()
Get field space.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
Deprecated interface functions.
DEPRECATED MoFEMErrorCode loop_finite_elements(const Problem *problem_ptr, const std::string &fe_name, FEMethod &method, int lower_rank, int upper_rank, MoFEMTypes bh, int verb=DEFAULT_VERBOSITY)
Data operator to do calculations at integration points.
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Create interface from given surface and insert flat prisms in-between.
MoFEMErrorCode getSides(const int msId, const CubitBCType cubit_bc_type, const BitRefLevel mesh_bit_level, const bool recursive, int verb=QUIET)
Store tetrahedra from each side of the interface separately in two child meshsets of the parent meshs...
MoFEMErrorCode splitSides(const EntityHandle meshset, const BitRefLevel &bit, const int msId, const CubitBCType cubit_bc_type, const bool add_interface_entities, const bool recursive=false, int verb=QUIET)
Split nodes and other entities of tetrahedra on both sides of the interface and insert flat prisms in...
Problem manager is used to build and partition problems.
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
VolumeElementForcesAndSourcesCoreOnContactPrismSideBase::UserDataOperator UserDataOperator