16 {
17
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;
29#if PETSC_VERSION_GE(3, 6, 4)
31 255, &flg);
32#else
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 = "";
43
44
46
48
49 enum side_contact { MASTERSIDE, SLAVESIDE };
50
51
53 enum bases { AINSWORTH, DEMKOWICZ, LASTBASEOP };
54 const char *list_bases[] = {"ainsworth", "demkowicz"};
55 PetscBool flg;
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)
65 return base;
66 }
68 };
69
72 const char *list_spaces[] = {"h1", "hdiv", "hcurl"};
73 PetscBool flg;
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)
85 return space;
86 }
88 };
89
91
92 auto add_prism_interface = [&](
Range &tets,
Range &prisms,
96 std::vector<BitRefLevel> &bit_levels) {
100
101 int ll = 1;
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);
110
111 {
112
114 CHKERR moab.create_meshset(MESHSET_SET, ref_level_meshset);
116 ->getEntitiesByTypeAndRefLevel(bit_levels.back(),
118 ref_level_meshset);
120 ->getEntitiesByTypeAndRefLevel(bit_levels.back(),
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
127 CHKERR interface->
getSides(cubit_meshset, bit_levels.back(),
true,
128 0);
129
131
133 cubit_meshset, true, true, 0);
134
135 CHKERR moab.delete_entities(&ref_level_meshset, 1);
136 }
137
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++) {
150 true);
151 }
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
159 ->getEntitiesByTypeAndRefLevel(bit_levels[0],
BitRefLevel().set(),
160 MBTET, meshset_tets);
161 CHKERR moab.get_entities_by_handle(meshset_tets, tets,
true);
162
164 ->getEntitiesByTypeAndRefLevel(bit_levels[0],
BitRefLevel().set(),
165 MBPRISM, meshset_prisms);
166 CHKERR moab.get_entities_by_handle(meshset_prisms, prisms);
167
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
175 };
176
177 Range all_tets, contact_prisms, master_tris, slave_tris;
179 std::vector<BitRefLevel> bit_levels;
180
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
191
194
195
200
203 }
204
206
207
216
218
221 bit_levels[0],
BitRefLevel().set(), MBPRISM, prism);
223
226
227
229
230
234
235
239
246 };
247
248 struct OnContactSideMaster
250
254
259 elemData(elem_data) {}
263
264 if (type == MBTRI && side == getFaceSideNumber()) {
265
267 getCoordsAtGaussPts() - getMasterCoordsAtGaussPts();
268 constexpr double eps = 1e-12;
269 if (norm_inf(diff) >
eps)
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
278 ptr_elem_data = &elemData.shapeFunH1VolSide;
279
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;
286 ++t_base;
287 }
288 }
289 }
290 }
292 }
293 };
294
300 FACEMASTER),
301 volSideFe(m_field), elemData(elem_data) {
303 new OnContactSideMaster::OpVolSide(elemData));
304 }
305
308
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,
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) {
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";
330
331 auto check_continuity_of_base = [&](auto &vol_dot_data) {
333
334 if (vol_dot_data.size1() != elemData.dotNormalFace.size1())
336 "Inconsistent number of integration points");
337
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));
350 }
352 };
353
354 auto check_continuity_of_h1_base = [&](auto &vol_data) {
356
357 if (vol_data.size1() != elemData.shapeFunH1Values.size1())
359 "Inconsistent number of integration points");
360
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));
373 }
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 }
385 }
386 };
387
388 struct OnContactSideSlave
390
394
399 elemData(elem_data) {}
403
404 if (type == MBTRI && side == getFaceSideNumber()) {
405
407 getCoordsAtGaussPts() - getMasterCoordsAtGaussPts();
408 constexpr double eps = 1e-12;
409 if (norm_inf(diff) >
eps)
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
418 ptr_elem_data = &elemData.shapeFunH1VolSide;
419
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;
426 ++t_base;
427 }
428 }
429 }
430 }
432 }
433 };
434
440 FACESLAVE),
441 volSideFe(m_field), elemData(elem_data) {
443 new OnContactSideSlave::OpVolSide(elemData));
444 }
445
448
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,
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) {
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";
471
472 auto check_continuity_of_base = [&](auto &vol_dot_data) {
474
475 if (vol_dot_data.size1() != elemData.dotNormalFace.size1())
477 "Inconsistent number of integration points");
478
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));
491 }
493 };
494
495 auto check_continuity_of_h1_base = [&](auto &vol_data) {
497
498 if (vol_data.size1() != elemData.shapeFunH1Values.size1())
500 "Inconsistent number of integration points");
501
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));
514 }
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 }
526 }
527 };
528
530
533
534
535 contact_prism_fe_master.getOpPtrVector().push_back(
536 new OnContactSideMaster(m_field, elem_data));
537
538 contact_prism_fe_slave.getOpPtrVector().push_back(
539 new OnContactSideSlave(m_field, elem_data));
540
543 }
545
547
548 return 0;
549}
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define CATCH_ERRORS
Catch errors.
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
FieldSpace
approximation spaces
@ LASTSPACE
FieldSpace in [ 0, LASTSPACE )
@ HCURL
field with continuous tangents
@ HDIV
field with continuous normal traction
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_ATOM_TEST_INVALID
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
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_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
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_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
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.
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.
#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 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
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
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
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)
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
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
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.
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
virtual MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
FieldSpace & getSpace()
Get field space.
EntityHandle getSideEntity(const int side_number, const EntityType type)
Get the side entity.
boost::ptr_deque< 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(IFACE *&iface) const
Get interface reference to pointer of interface.