29 {
30
31
33
34
35 moab::Core mb_instance;
36 moab::Interface &moab = mb_instance;
37
38 try {
39 PetscBool flg_file;
42 PetscInt order_lambda = 1;
43 PetscReal r_value = 1.;
44 PetscReal cn_value = 1.;
45 PetscBool is_newton_cotes = PETSC_FALSE;
46 PetscBool test_jacobian = PETSC_FALSE;
47 PetscBool convect_pts = PETSC_FALSE;
48 PetscBool test_ale = PETSC_FALSE;
49 PetscBool alm_flag = PETSC_FALSE;
50 PetscBool eigen_pos_flag = PETSC_FALSE;
51 PetscBool use_reference_coordinates = PETSC_TRUE;
52
53 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Elastic Config",
"none");
54
56 PETSC_NULL);
57
58 CHKERR PetscOptionsString(
"-my_file",
"mesh file name",
"",
"mesh.h5m",
60
61 CHKERR PetscOptionsInt(
"-my_order",
62 "approximation order for spatial positions", "", 1,
65 "-my_order_lambda",
66 "default approximation order of Lagrange multipliers", "", 1,
67 &order_lambda, PETSC_NULL);
68
69 CHKERR PetscOptionsReal(
"-my_cn_value",
"default regularisation cn value",
70 "", 1., &cn_value, PETSC_NULL);
71
72 CHKERR PetscOptionsBool(
"-my_is_newton_cotes",
73 "set if Newton-Cotes integration rules are used",
74 "", PETSC_FALSE, &is_newton_cotes, PETSC_NULL);
75 CHKERR PetscOptionsBool(
"-my_convect",
"set to convect integration pts",
"",
76 PETSC_FALSE, &convect_pts, PETSC_NULL);
77 CHKERR PetscOptionsBool(
"-my_alm_flag",
"set to convect integration pts",
78 "", PETSC_FALSE, &alm_flag, PETSC_NULL);
79
80 CHKERR PetscOptionsBool(
"-my_eigen_pos_flag",
81 "if set use eigen spatial positions are taken into "
82 "account for predeformed configuration",
83 "", PETSC_FALSE, &eigen_pos_flag, PETSC_NULL);
84
86 PETSC_NULL);
87
89 &use_reference_coordinates, PETSC_NULL);
90
91 ierr = PetscOptionsEnd();
93
94
95 if (flg_file != PETSC_TRUE) {
96 SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
97 }
98
99
100 const char *option;
101 option = "";
103
104
107
108 auto add_prism_interface = [&](
Range &contact_prisms,
Range &master_tris,
110 std::vector<BitRefLevel> &bit_levels) {
114
116 if (cit->getName().compare(0, 11, "INT_CONTACT") == 0) {
117 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Insert %s (id: %d)\n",
118 cit->getName().c_str(), cit->getMeshsetId());
120
121
123 CHKERR moab.create_meshset(MESHSET_SET, ref_level_meshset);
125 ->getEntitiesByTypeAndRefLevel(bit_levels.back(),
127 ref_level_meshset);
129 ->getEntitiesByTypeAndRefLevel(bit_levels.back(),
131 ref_level_meshset);
132
133
134 CHKERR interface->
getSides(cubit_meshset, bit_levels.back(),
true, 0);
135
136 bit_levels.push_back(
BitRefLevel().set(bit_levels.size()));
137
139 cubit_meshset, true, true, 0);
140
141 CHKERR moab.delete_entities(&ref_level_meshset, 1);
142
144 bit_levels.pop_back();
145 }
146 }
147
149 CHKERR moab.create_meshset(MESHSET_SET, meshset_prisms);
151 ->getEntitiesByTypeAndRefLevel(bit_levels.back(),
BitRefLevel().set(),
152 MBPRISM, meshset_prisms);
153 CHKERR moab.get_entities_by_handle(meshset_prisms, contact_prisms);
154 CHKERR moab.delete_entities(&meshset_prisms, 1);
155
157 for (Range::iterator pit = contact_prisms.begin();
158 pit != contact_prisms.end(); pit++) {
159 CHKERR moab.side_element(*pit, 2, 3, tri);
160 master_tris.insert(tri);
161 CHKERR moab.side_element(*pit, 2, 4, tri);
162 slave_tris.insert(tri);
163 }
164
166 };
167
168 Range contact_prisms, master_tris, slave_tris;
169 std::vector<BitRefLevel> bit_levels;
170
173 0, 3, bit_levels.back());
174
175 CHKERR add_prism_interface(contact_prisms, master_tris, slave_tris,
176 bit_levels);
177
180
181
182
188
191
196
199
205
206 if (eigen_pos_flag) {
214 }
215
216
218
219 PetscRandom rctx;
220 PetscRandomCreate(PETSC_COMM_WORLD, &rctx);
221
222 auto set_coord = [&](
VectorAdaptor &&field_data,
double *x,
double *y,
223 double *z) {
225 double value;
227 PetscRandomGetValue(rctx, &value);
228 field_data[0] = (*x) + (value - 0.5) *
scale;
229 PetscRandomGetValue(rctx, &value);
230 field_data[1] = (*y) + (value - 0.5) *
scale;
231 PetscRandomGetValue(rctx, &value);
232 field_data[2] = (*z) + (value - 0.5) *
scale;
234 };
235
236 auto set_pressure = [&](
VectorAdaptor &&field_data,
double *x,
double *y,
237 double *z) {
239 double value;
241 PetscRandomGetValueReal(rctx, &value);
242 field_data[0] = value *
scale;
244 };
245
247 "SPATIAL_POSITION");
249 "LAGMULT");
250
251 if (eigen_pos_flag) {
253 set_coord, "SPATIAL_POSITION");
254 }
255
256 if (test_ale == PETSC_TRUE) {
258 set_coord, "MESH_NODE_POSITIONS");
259 } else {
260
261 {
264 }
265 }
266
267 PetscRandomDestroy(&rctx);
268
269 auto cn_value_ptr = boost::make_shared<double>(cn_value);
270 auto contact_problem = boost::make_shared<SimpleContactProblem>(
271 m_field, cn_value_ptr, is_newton_cotes);
272
273 auto make_contact_element = [&]() {
274 return boost::make_shared<SimpleContactProblem::SimpleContactElement>(
275 m_field);
276 };
277
278 auto make_convective_master_element = [&]() {
279 return boost::make_shared<
281 m_field, "SPATIAL_POSITION", "MESH_NODE_POSITIONS");
282 };
283
284 auto make_convective_slave_element = [&]() {
285 return boost::make_shared<
287 m_field, "SPATIAL_POSITION", "MESH_NODE_POSITIONS");
288 };
289
290 auto make_contact_common_data = [&]() {
291 return boost::make_shared<SimpleContactProblem::CommonDataSimpleContact>(
292 m_field);
293 };
294
295 auto get_contact_rhs = [&](auto contact_problem, auto make_element,
296 bool is_alm = false) {
297 auto fe_rhs_simple_contact = make_element();
298 auto common_data_simple_contact = make_contact_common_data();
299 contact_problem->setContactOperatorsRhs(
300 fe_rhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
301 "LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS",
302 use_reference_coordinates);
303 return fe_rhs_simple_contact;
304 };
305
306 auto get_master_contact_lhs = [&](auto contact_problem, auto make_element,
307 bool is_alm = false) {
308 auto fe_lhs_simple_contact = make_element();
309 auto common_data_simple_contact = make_contact_common_data();
310 contact_problem->setContactOperatorsLhs(
311 fe_lhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
312 "LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS",
313 use_reference_coordinates);
314 return fe_lhs_simple_contact;
315 };
316
317 auto get_master_traction_rhs = [&](auto contact_problem, auto make_element,
318 bool alm_flag = false) {
319 auto fe_rhs_simple_contact = make_element();
320 auto common_data_simple_contact = make_contact_common_data();
321 contact_problem->setMasterForceOperatorsRhs(
322 fe_rhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
323 "LAGMULT", alm_flag, eigen_pos_flag, "EIGEN_POSITIONS",
324 use_reference_coordinates);
325 return fe_rhs_simple_contact;
326 };
327
328 auto get_master_traction_lhs = [&](auto contact_problem, auto make_element,
329 bool alm_flag = false) {
330 auto fe_lhs_simple_contact = make_element();
331 auto common_data_simple_contact = make_contact_common_data();
332 contact_problem->setMasterForceOperatorsLhs(
333 fe_lhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
334 "LAGMULT", alm_flag, eigen_pos_flag, "EIGEN_POSITIONS",
335 use_reference_coordinates);
336 return fe_lhs_simple_contact;
337 };
338
339 auto get_contact_material_rhs = [&](auto contact_problem, auto make_element,
341 auto fe_rhs_simple_contact_ale_material = make_element();
342 auto common_data_simple_contact = make_contact_common_data();
343 common_data_simple_contact->forcesOnlyOnEntitiesRow.clear();
344 common_data_simple_contact->forcesOnlyOnEntitiesRow = ale_nodes;
345 contact_problem->setContactOperatorsRhsALEMaterial(
346 fe_rhs_simple_contact_ale_material, common_data_simple_contact,
347 "SPATIAL_POSITION", "MESH_NODE_POSITIONS", "LAGMULT", "MATERIAL");
348 return fe_rhs_simple_contact_ale_material;
349 };
350
351 auto get_simple_contact_ale_lhs = [&](auto contact_problem,
352 auto make_element) {
353 auto fe_lhs_simple_contact_ale = make_element();
354 auto common_data_simple_contact = make_contact_common_data();
355 contact_problem->setContactOperatorsLhsALE(
356 fe_lhs_simple_contact_ale, common_data_simple_contact,
357 "SPATIAL_POSITION", "MESH_NODE_POSITIONS", "LAGMULT", eigen_pos_flag,
358 "EIGEN_POSITIONS");
359 return fe_lhs_simple_contact_ale;
360 };
361
362 auto get_simple_contact_ale_material_lhs =
363 [&](
auto contact_problem,
auto make_element,
Range &ale_nodes) {
364 auto fe_lhs_simple_contact_material_ale = make_element();
365 auto common_data_simple_contact = make_contact_common_data();
366 common_data_simple_contact->forcesOnlyOnEntitiesRow.clear();
367 common_data_simple_contact->forcesOnlyOnEntitiesRow = ale_nodes;
368 contact_problem->setContactOperatorsLhsALEMaterial(
369 fe_lhs_simple_contact_material_ale, common_data_simple_contact,
370 "SPATIAL_POSITION", "MESH_NODE_POSITIONS", "LAGMULT", "MATERIAL");
371 return fe_lhs_simple_contact_material_ale;
372 };
373
374
375 if (!eigen_pos_flag)
376 contact_problem->addContactElement("CONTACT_ELEM", "SPATIAL_POSITION",
377 "LAGMULT", contact_prisms);
378 else
379 contact_problem->addContactElement("CONTACT_ELEM", "SPATIAL_POSITION",
380 "LAGMULT", contact_prisms,
381 eigen_pos_flag, "EIGEN_POSITIONS");
382
384 if (test_ale == PETSC_TRUE) {
385 if (!eigen_pos_flag)
386 contact_problem->addContactElementALE(
387 "ALE_CONTACT_ELEM", "SPATIAL_POSITION", "MESH_NODE_POSITIONS",
388 "LAGMULT", contact_prisms);
389 else
390 contact_problem->addContactElementALE(
391 "ALE_CONTACT_ELEM", "SPATIAL_POSITION", "MESH_NODE_POSITIONS",
392 "LAGMULT", contact_prisms, eigen_pos_flag, "EIGEN_POSITIONS");
393
395 CHKERR moab.get_adjacencies(contact_prisms, 2,
false, faces,
396 moab::Interface::UNION);
397 Range tris = faces.subset_by_type(MBTRI);
398
399 CHKERR moab.get_adjacencies(tris, 3,
false, all_tets,
400 moab::Interface::UNION);
401
402
405 "SPATIAL_POSITION");
407 "SPATIAL_POSITION");
409 "MESH_NODE_POSITIONS");
411 "MESH_NODE_POSITIONS");
413 "SPATIAL_POSITION");
415 "MATERIAL", "MESH_NODE_POSITIONS");
417 "MATERIAL");
419 }
420
421
423
424
426
427
429
430
432 bit_levels.back());
433
434 DMType dm_name = "DMMOFEM";
436
437
440 CHKERR DMSetType(dm, dm_name);
441
442
444 CHKERR DMSetFromOptions(dm);
446
448
449 if (test_ale == PETSC_TRUE) {
452 }
453
455
456
459
460
462
464 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
465 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
466
468 CHKERR VecGhostUpdateBegin(
F, INSERT_VALUES, SCATTER_FORWARD);
469 CHKERR VecGhostUpdateEnd(
F, INSERT_VALUES, SCATTER_FORWARD);
470
471 CHKERR MatSetOption(A, MAT_SPD, PETSC_TRUE);
473
475
476 if (convect_pts == PETSC_TRUE) {
478 dm, "CONTACT_ELEM",
479 get_contact_rhs(contact_problem, make_convective_master_element),
480 PETSC_NULL, PETSC_NULL);
482 dm, "CONTACT_ELEM",
483 get_master_contact_lhs(contact_problem,
484 make_convective_master_element),
485 NULL, NULL);
487 dm, "CONTACT_ELEM",
488 get_master_traction_rhs(contact_problem,
489 make_convective_slave_element),
490 PETSC_NULL, PETSC_NULL);
492 dm, "CONTACT_ELEM",
493 get_master_traction_lhs(contact_problem,
494 make_convective_slave_element),
495 NULL, NULL);
496 } else {
498 dm, "CONTACT_ELEM",
499 get_contact_rhs(contact_problem, make_contact_element, alm_flag),
500 PETSC_NULL, PETSC_NULL);
502 dm, "CONTACT_ELEM",
503 get_master_traction_rhs(contact_problem, make_contact_element,
504 alm_flag),
505 PETSC_NULL, PETSC_NULL);
507 get_master_contact_lhs(contact_problem,
508 make_contact_element,
509 alm_flag),
510 PETSC_NULL, PETSC_NULL);
512 dm, "CONTACT_ELEM",
513 get_master_traction_lhs(contact_problem, make_contact_element,
514 alm_flag),
515 PETSC_NULL, PETSC_NULL);
516 }
517
518 if (test_ale == PETSC_TRUE) {
520 CHKERR moab.get_connectivity(all_tets, nodes,
false);
521
523 dm, "ALE_CONTACT_ELEM",
524 get_contact_material_rhs(contact_problem, make_contact_element,
525 nodes),
526 PETSC_NULL, PETSC_NULL);
527
529 dm, "ALE_CONTACT_ELEM",
530 get_simple_contact_ale_lhs(contact_problem, make_contact_element),
531 NULL, NULL);
532
534 dm, "ALE_CONTACT_ELEM",
535 get_simple_contact_ale_material_lhs(contact_problem,
536 make_contact_element, nodes),
537 NULL, NULL);
538 }
539
540 if (test_jacobian == PETSC_TRUE) {
541 char testing_options[] =
542 "-snes_test_jacobian -snes_test_jacobian_display "
543 "-snes_no_convergence_test -snes_atol 0 -snes_rtol 0 "
544 "-snes_max_it "
545 "1 ";
546 CHKERR PetscOptionsInsertString(NULL, testing_options);
547 } else {
548 char testing_options[] = "-snes_no_convergence_test -snes_atol 0 "
549 "-snes_rtol 0 "
550 "-snes_max_it 1 ";
551 CHKERR PetscOptionsInsertString(NULL, testing_options);
552 }
553
555 SNESConvergedReason snes_reason;
557
558
559 {
563 CHKERR SNESSetFromOptions(snes);
564 }
565
566 CHKERR SNESSolve(snes, PETSC_NULL,
D);
567
568 if (test_jacobian == PETSC_FALSE) {
569 double nrm_A0;
570 CHKERR MatNorm(A, NORM_INFINITY, &nrm_A0);
571
572 char testing_options_fd[] = "-snes_fd";
573 CHKERR PetscOptionsInsertString(NULL, testing_options_fd);
574
577 CHKERR SNESSetFromOptions(snes);
578
579 CHKERR SNESSolve(snes, NULL,
D);
580 CHKERR MatAXPY(A, -1, fdA, SUBSET_NONZERO_PATTERN);
581
582 double nrm_A;
583 CHKERR MatNorm(A, NORM_INFINITY, &nrm_A);
584 PetscPrintf(PETSC_COMM_WORLD, "Matrix norms %3.4e %3.4e\n", nrm_A,
585 nrm_A / nrm_A0);
586 nrm_A /= nrm_A0;
587
588 constexpr double tol = 1e-6;
591 "Difference between hand-calculated tangent matrix and finite "
592 "difference matrix is too big");
593 }
594 }
595 }
597
598
600
601 return 0;
602}
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#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.
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
PetscErrorCode DMMoFEMSNESSetFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES residual evaluation function
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.
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMMoFEMSNESSetJacobian(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES Jacobian evaluation function
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
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_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
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 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_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
#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.
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
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
VectorShallowArrayAdaptor< double > VectorAdaptor
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
DEPRECATED auto smartCreateDMMatrix(DM dm)
auto createSNES(MPI_Comm comm)
PetscErrorCode SnesMat(SNES snes, Vec x, Mat A, Mat B, void *ctx)
This is MoFEM implementation for the left hand side (tangent matrix) evaluation in SNES solver.
DEPRECATED SmartPetscObj< Vec > smartVectorDuplicate(Vec vec)
DEPRECATED auto createSmartDM(MPI_Comm comm, const std::string dm_type_name)
PetscErrorCode SnesRhs(SNES snes, Vec x, Vec f, void *ctx)
This is MoFEM implementation for the right hand side (residual vector) evaluation in SNES solver.
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
DEPRECATED auto smartCreateDMVector(DM dm)
DEPRECATED SmartPetscObj< Mat > smartMatDuplicate(Mat mat, MatDuplicateOption op)
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MPI_Comm & get_comm() const =0
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.
Deprecated interface functions.
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...
Projection of edge entities with one mid-node on hierarchical basis.
intrusive_ptr for managing petsc objects
Interface for nonlinear (SNES) solver.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.