25 MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
27 PetscBool flg = PETSC_TRUE;
29 #if PETSC_VERSION_GE(3, 6, 4)
36 if (flg != PETSC_TRUE) {
37 SETERRQ(PETSC_COMM_SELF, 1,
"*** ERROR -my_file (MESH FILE NEEDED)");
49 enum side_contact { MASTERSIDE, SLAVESIDE };
53 enum bases { AINSWORTH, DEMKOWICZ, LASTBASEOP };
54 const char *list_bases[] = {
"ainsworth",
"demkowicz"};
56 PetscInt choice_base_value = AINSWORTH;
58 LASTBASEOP, &choice_base_value, &flg);
59 if (flg == PETSC_TRUE) {
61 if (choice_base_value == AINSWORTH)
63 else if (choice_base_value == DEMKOWICZ)
72 const char *list_spaces[] = {
"h1",
"hdiv",
"hcurl"};
74 PetscInt choice_space_value =
HDIV;
76 LASTSPACEOP, &choice_space_value, &flg);
77 if (flg == PETSC_TRUE) {
79 if (choice_space_value ==
HDIV)
81 else if (choice_space_value ==
HCURL)
83 else if (choice_space_value ==
H1)
92 auto add_prism_interface = [&](
Range &tets,
Range &prisms,
96 std::vector<BitRefLevel> &bit_levels) {
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());
108 CHKERR moab.get_entities_by_type(cubit_meshset, MBTRI, tris,
true);
109 master_tris.merge(tris);
114 CHKERR moab.create_meshset(MESHSET_SET, ref_level_meshset);
116 ->getEntitiesByTypeAndRefLevel(bit_levels.back(),
120 ->getEntitiesByTypeAndRefLevel(bit_levels.back(),
123 Range ref_level_tets;
124 CHKERR moab.get_entities_by_handle(ref_level_meshset,
125 ref_level_tets,
true);
127 CHKERR interface->
getSides(cubit_meshset, bit_levels.back(),
true,
133 cubit_meshset,
true,
true, 0);
135 CHKERR moab.delete_entities(&ref_level_meshset, 1);
141 ->updateMeshsetByEntitiesChildren(
142 cubit_meshset, bit_levels.back(), cubit_meshset, MBMAXTYPE,
148 for (
unsigned int ll = 0; ll != bit_levels.size() - 1; ll++) {
153 bit_levels.size() - 1);
155 CHKERR moab.create_meshset(MESHSET_SET, meshset_tets);
156 CHKERR moab.create_meshset(MESHSET_SET, meshset_prisms);
159 ->getEntitiesByTypeAndRefLevel(bit_levels[0],
BitRefLevel().set(),
160 MBTET, meshset_tets);
161 CHKERR moab.get_entities_by_handle(meshset_tets, tets,
true);
164 ->getEntitiesByTypeAndRefLevel(bit_levels[0],
BitRefLevel().set(),
165 MBPRISM, meshset_prisms);
166 CHKERR moab.get_entities_by_handle(meshset_prisms, prisms);
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);
177 Range all_tets, contact_prisms, master_tris, slave_tris;
179 std::vector<BitRefLevel> bit_levels;
183 0, 3, bit_levels[0]);
185 CHKERR add_prism_interface(all_tets, contact_prisms, master_tris,
186 slave_tris, meshset_tets, meshset_prisms,
221 bit_levels[0],
BitRefLevel().set(), MBPRISM, prism);
248 struct OnContactSideMaster
259 elemData(elem_data) {}
264 if (
type == MBTRI && side == getFaceSideNumber()) {
267 getCoordsAtGaussPts() - getMasterCoordsAtGaussPts();
268 constexpr
double eps = 1e-12;
269 if (norm_inf(diff) >
eps)
271 "coordinates at integration pts are different");
273 const size_t nb_dofs = data.
getN().size2() / 3;
274 const size_t nb_integration_pts = data.
getN().size1();
278 ptr_elem_data = &elemData.shapeFunH1VolSide;
281 elem_data.resize(nb_integration_pts, nb_dofs,
false);
282 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
284 for (
size_t bb = 0; bb != nb_dofs; ++bb) {
285 elem_data(gg, bb) = t_base;
301 volSideFe(m_field), elemData(elem_data) {
303 new OnContactSideMaster::OpVolSide(elemData));
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();
314 elemData.shapeFunH1Values.resize(nb_integration_pts, nb_dofs,
316 elemData.shapeFunH1VolSide.resize(nb_integration_pts, nb_dofs,
319 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
321 for (
size_t bb = 0; bb != nb_dofs; ++bb) {
322 elemData.shapeFunH1Values(gg, bb) = t_base;
327 std::string side_fe_name =
"V1";
331 auto check_continuity_of_base = [&](
auto &vol_dot_data) {
334 if (vol_dot_data.size1() != elemData.dotNormalFace.size1())
336 "Inconsistent number of integration points");
338 if (vol_dot_data.size2() != elemData.dotNormalFace.size2())
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));
348 "Inconsistency %3.4e != %3.4e", vol_dot_data(gg, bb),
349 elemData.dotNormalFace(gg, bb));
354 auto check_continuity_of_h1_base = [&](
auto &vol_data) {
357 if (vol_data.size1() != elemData.shapeFunH1Values.size1())
359 "Inconsistent number of integration points");
361 if (vol_data.size2() != elemData.shapeFunH1Values.size2())
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));
371 "Inconsistency %3.4e != %3.4e", vol_data(gg, bb),
372 elemData.shapeFunH1Values(gg, bb));
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);
388 struct OnContactSideSlave
399 elemData(elem_data) {}
404 if (
type == MBTRI && side == getFaceSideNumber()) {
407 getCoordsAtGaussPts() - getMasterCoordsAtGaussPts();
408 constexpr
double eps = 1e-12;
409 if (norm_inf(diff) >
eps)
411 "coordinates at integration pts are different");
413 const size_t nb_dofs = data.
getN().size2() / 3;
414 const size_t nb_integration_pts = data.
getN().size1();
418 ptr_elem_data = &elemData.shapeFunH1VolSide;
421 elem_data.resize(nb_integration_pts, nb_dofs,
false);
422 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
424 for (
size_t bb = 0; bb != nb_dofs; ++bb) {
425 elem_data(gg, bb) = t_base;
441 volSideFe(m_field), elemData(elem_data) {
443 new OnContactSideSlave::OpVolSide(elemData));
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();
454 elemData.shapeFunH1Values.resize(nb_integration_pts, nb_dofs,
456 elemData.shapeFunH1VolSide.resize(nb_integration_pts, nb_dofs,
459 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
461 for (
size_t bb = 0; bb != nb_dofs; ++bb) {
462 elemData.shapeFunH1Values(gg, bb) = t_base;
468 std::string side_fe_name =
"V1";
472 auto check_continuity_of_base = [&](
auto &vol_dot_data) {
475 if (vol_dot_data.size1() != elemData.dotNormalFace.size1())
477 "Inconsistent number of integration points");
479 if (vol_dot_data.size2() != elemData.dotNormalFace.size2())
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));
489 "Inconsistency %3.4e != %3.4e", vol_dot_data(gg, bb),
490 elemData.dotNormalFace(gg, bb));
495 auto check_continuity_of_h1_base = [&](
auto &vol_data) {
498 if (vol_data.size1() != elemData.shapeFunH1Values.size1())
500 "Inconsistent number of integration points");
502 if (vol_data.size2() != elemData.shapeFunH1Values.size2())
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));
512 "Inconsistency %3.4e != %3.4e", vol_data(gg, bb),
513 elemData.shapeFunH1Values(gg, bb));
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);
535 contact_prism_fe_master.getOpPtrVector().push_back(
536 new OnContactSideMaster(m_field, elem_data));
538 contact_prism_fe_slave.getOpPtrVector().push_back(
539 new OnContactSideSlave(m_field, elem_data));