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