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;
53 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Elastic Config",
"none");
58 CHKERR PetscOptionsString(
"-my_file",
"mesh file name",
"",
"mesh.h5m",
61 CHKERR PetscOptionsInt(
"-my_order",
62 "approximation order for spatial positions",
"", 1,
66 "default approximation order of Lagrange multipliers",
"", 1,
67 &order_lambda, PETSC_NULL);
69 CHKERR PetscOptionsReal(
"-my_cn_value",
"default regularisation cn value",
70 "", 1., &cn_value, PETSC_NULL);
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);
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);
89 &use_reference_coordinates, PETSC_NULL);
91 ierr = PetscOptionsEnd();
95 if (flg_file != PETSC_TRUE) {
96 SETERRQ(PETSC_COMM_SELF, 1,
"*** ERROR -my_file (MESH FILE NEEDED)");
108 auto add_prism_interface = [&](
Range &contact_prisms,
Range &master_tris,
110 std::vector<BitRefLevel> &bit_levels) {
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());
123 CHKERR moab.create_meshset(MESHSET_SET, ref_level_meshset);
125 ->getEntitiesByTypeAndRefLevel(bit_levels.back(),
129 ->getEntitiesByTypeAndRefLevel(bit_levels.back(),
134 CHKERR interface->
getSides(cubit_meshset, bit_levels.back(),
true, 0);
136 bit_levels.push_back(
BitRefLevel().set(bit_levels.size()));
139 cubit_meshset,
true,
true, 0);
141 CHKERR moab.delete_entities(&ref_level_meshset, 1);
144 bit_levels.pop_back();
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);
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);
168 Range contact_prisms, master_tris, slave_tris;
169 std::vector<BitRefLevel> bit_levels;
173 0, 3, bit_levels.back());
175 CHKERR add_prism_interface(contact_prisms, master_tris, slave_tris,
206 if (eigen_pos_flag) {
220 PetscRandomCreate(PETSC_COMM_WORLD, &rctx);
222 auto set_coord = [&](
VectorAdaptor &&field_data,
double *x,
double *y,
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;
236 auto set_pressure = [&](
VectorAdaptor &&field_data,
double *x,
double *y,
241 PetscRandomGetValueReal(rctx, &value);
242 field_data[0] = value *
scale;
251 if (eigen_pos_flag) {
253 set_coord,
"SPATIAL_POSITION");
256 if (test_ale == PETSC_TRUE) {
258 set_coord,
"MESH_NODE_POSITIONS");
267 PetscRandomDestroy(&rctx);
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);
273 auto make_contact_element = [&]() {
274 return boost::make_shared<SimpleContactProblem::SimpleContactElement>(
278 auto make_convective_master_element = [&]() {
279 return boost::make_shared<
281 m_field,
"SPATIAL_POSITION",
"MESH_NODE_POSITIONS");
284 auto make_convective_slave_element = [&]() {
285 return boost::make_shared<
287 m_field,
"SPATIAL_POSITION",
"MESH_NODE_POSITIONS");
290 auto make_contact_common_data = [&]() {
291 return boost::make_shared<SimpleContactProblem::CommonDataSimpleContact>(
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;
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;
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;
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;
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;
351 auto get_simple_contact_ale_lhs = [&](
auto contact_problem,
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,
359 return fe_lhs_simple_contact_ale;
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;
376 contact_problem->addContactElement(
"CONTACT_ELEM",
"SPATIAL_POSITION",
377 "LAGMULT", contact_prisms);
379 contact_problem->addContactElement(
"CONTACT_ELEM",
"SPATIAL_POSITION",
380 "LAGMULT", contact_prisms,
381 eigen_pos_flag,
"EIGEN_POSITIONS");
384 if (test_ale == PETSC_TRUE) {
386 contact_problem->addContactElementALE(
387 "ALE_CONTACT_ELEM",
"SPATIAL_POSITION",
"MESH_NODE_POSITIONS",
388 "LAGMULT", contact_prisms);
390 contact_problem->addContactElementALE(
391 "ALE_CONTACT_ELEM",
"SPATIAL_POSITION",
"MESH_NODE_POSITIONS",
392 "LAGMULT", contact_prisms, eigen_pos_flag,
"EIGEN_POSITIONS");
395 CHKERR moab.get_adjacencies(contact_prisms, 2,
false, faces,
396 moab::Interface::UNION);
397 Range tris = faces.subset_by_type(MBTRI);
399 CHKERR moab.get_adjacencies(tris, 3,
false, all_tets,
400 moab::Interface::UNION);
409 "MESH_NODE_POSITIONS");
411 "MESH_NODE_POSITIONS");
415 "MATERIAL",
"MESH_NODE_POSITIONS");
434 DMType dm_name =
"DMMOFEM";
440 CHKERR DMSetType(dm, dm_name);
444 CHKERR DMSetFromOptions(dm);
449 if (test_ale == PETSC_TRUE) {
464 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
465 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
468 CHKERR VecGhostUpdateBegin(
F, INSERT_VALUES, SCATTER_FORWARD);
469 CHKERR VecGhostUpdateEnd(
F, INSERT_VALUES, SCATTER_FORWARD);
471 CHKERR MatSetOption(
A, MAT_SPD, PETSC_TRUE);
476 if (convect_pts == PETSC_TRUE) {
479 get_contact_rhs(contact_problem, make_convective_master_element),
480 PETSC_NULL, PETSC_NULL);
483 get_master_contact_lhs(contact_problem,
484 make_convective_master_element),
488 get_master_traction_rhs(contact_problem,
489 make_convective_slave_element),
490 PETSC_NULL, PETSC_NULL);
493 get_master_traction_lhs(contact_problem,
494 make_convective_slave_element),
499 get_contact_rhs(contact_problem, make_contact_element, alm_flag),
500 PETSC_NULL, PETSC_NULL);
503 get_master_traction_rhs(contact_problem, make_contact_element,
505 PETSC_NULL, PETSC_NULL);
507 get_master_contact_lhs(contact_problem,
508 make_contact_element,
510 PETSC_NULL, PETSC_NULL);
513 get_master_traction_lhs(contact_problem, make_contact_element,
515 PETSC_NULL, PETSC_NULL);
518 if (test_ale == PETSC_TRUE) {
520 CHKERR moab.get_connectivity(all_tets, nodes,
false);
523 dm,
"ALE_CONTACT_ELEM",
524 get_contact_material_rhs(contact_problem, make_contact_element,
526 PETSC_NULL, PETSC_NULL);
529 dm,
"ALE_CONTACT_ELEM",
530 get_simple_contact_ale_lhs(contact_problem, make_contact_element),
534 dm,
"ALE_CONTACT_ELEM",
535 get_simple_contact_ale_material_lhs(contact_problem,
536 make_contact_element, nodes),
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 "
546 CHKERR PetscOptionsInsertString(NULL, testing_options);
548 char testing_options[] =
"-snes_no_convergence_test -snes_atol 0 "
551 CHKERR PetscOptionsInsertString(NULL, testing_options);
555 SNESConvergedReason snes_reason;
563 CHKERR SNESSetFromOptions(snes);
566 CHKERR SNESSolve(snes, PETSC_NULL,
D);
568 if (test_jacobian == PETSC_FALSE) {
570 CHKERR MatNorm(
A, NORM_INFINITY, &nrm_A0);
572 char testing_options_fd[] =
"-snes_fd";
573 CHKERR PetscOptionsInsertString(NULL, testing_options_fd);
577 CHKERR SNESSetFromOptions(snes);
579 CHKERR SNESSolve(snes, NULL,
D);
580 CHKERR MatAXPY(
A, -1, fdA, SUBSET_NONZERO_PATTERN);
583 CHKERR MatNorm(
A, NORM_INFINITY, &nrm_A);
584 PetscPrintf(PETSC_COMM_WORLD,
"Matrix norms %3.4e %3.4e\n", nrm_A,
588 constexpr
double tol = 1e-6;
591 "Difference between hand-calculated tangent matrix and finite "
592 "difference matrix is too big");