v0.14.0
Loading...
Searching...
No Matches
Classes | Typedefs | Functions | Variables
mortar_contact_thermal.cpp File Reference
#include <Mortar.hpp>
#include <HookeInternalStressElement.hpp>

Go to the source code of this file.

Classes

struct  VolRule
 Set integration rule. More...
 

Typedefs

using BlockData = NonlinearElasticElement::BlockData
 
using VolSideFe = VolumeElementForcesAndSourcesCoreOnSide
 

Functions

int main (int argc, char *argv[])
 

Variables

static char help [] = "\n"
 

Typedef Documentation

◆ BlockData

Definition at line 28 of file mortar_contact_thermal.cpp.

◆ VolSideFe

Definition at line 29 of file mortar_contact_thermal.cpp.

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)

Post proc on the skin

Definition at line 35 of file mortar_contact_thermal.cpp.

35 {
36
37 const string default_options = "-ksp_type fgmres \n"
38 "-pc_type lu \n"
39 "-pc_factor_mat_solver_type mumps \n"
40 "-mat_mumps_icntl_13 1 \n"
41 "-mat_mumps_icntl_14 800 \n"
42 "-mat_mumps_icntl_20 0 \n"
43 "-mat_mumps_icntl_24 1 \n"
44 "-snes_type newtonls \n"
45 "-snes_linesearch_type basic \n"
46 "-snes_divergence_tolerance 0 \n"
47 "-snes_max_it 50 \n"
48 "-snes_atol 1e-8 \n"
49 "-snes_rtol 1e-10 \n"
50 "-snes_monitor \n"
51 "-ksp_monitor \n"
52 "-snes_converged_reason \n"
53 "-order 1 \n"
54 "-order_lambda 1 \n"
55 "-order_contact 2 \n"
56 "-ho_levels_num 1 \n"
57 "-step_num 1 \n"
58 "-cn_value 1. \n"
59 "-r_value 1. \n"
60 "-alm_flag 0 \n"
61 "-eigen_pos_flag 0 \n";
62
63 string param_file = "param_file.petsc";
64 if (!static_cast<bool>(ifstream(param_file))) {
65 std::ofstream file(param_file.c_str(), std::ios::ate);
66 if (file.is_open()) {
67 file << default_options;
68 file.close();
69 }
70 }
71
72 // Initialize MoFEM
73 MoFEM::Core::Initialize(&argc, &argv, param_file.c_str(), help);
74
75 // Create mesh database
76 moab::Core mb_instance; // create database
77 moab::Interface &moab = mb_instance; // create interface to database
78
79 try {
80 PetscBool flg_file;
81 PetscBool flg_file_out;
82
83 char mesh_file_name[255];
84 char output_mesh_name[255];
85 PetscInt order = 1;
86 PetscInt order_contact = 1;
87 PetscInt nb_ho_levels = 0;
88 PetscInt order_lambda = 1;
89 PetscReal r_value = 1.;
90 PetscReal cn_value = -1;
91 PetscInt nb_sub_steps = 1;
92 PetscBool is_partitioned = PETSC_FALSE;
93 PetscBool is_newton_cotes = PETSC_FALSE;
94 PetscInt test_num = 0;
95 PetscBool convect_pts = PETSC_FALSE;
96 PetscBool print_contact_state = PETSC_FALSE;
97 PetscBool alm_flag = PETSC_FALSE;
98 PetscBool eigen_pos_flag = PETSC_FALSE;
99
100 PetscReal thermal_expansion_coef = 1e-5;
101 PetscReal init_temp = 250.0;
102 PetscReal scale_factor = 1.0;
103 PetscBool analytical_input = PETSC_TRUE;
104 char stress_tag_name[255];
105 PetscBool flg_tag_name;
106 PetscBool save_mean_stress = PETSC_TRUE;
107 PetscBool ignore_contact = PETSC_FALSE;
108 PetscBool ignore_pressure = PETSC_FALSE;
109
110 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Elastic Config", "none");
111
112 CHKERR PetscOptionsString("-file_name", "mesh file name", "", "mesh.h5m",
113 mesh_file_name, 255, &flg_file);
114 CHKERR PetscOptionsString("-output_mesh_file", "output mesh file name", "",
115 "mesh.h5m", output_mesh_name, 255, &flg_file_out);
116
117 CHKERR PetscOptionsInt("-order", "approximation order of spatial positions",
118 "", 1, &order, PETSC_NULL);
119 CHKERR PetscOptionsInt(
120 "-order_contact",
121 "approximation order of spatial positions in contact interface", "", 1,
122 &order_contact, PETSC_NULL);
123 CHKERR PetscOptionsInt("-ho_levels_num", "number of higher order levels",
124 "", 0, &nb_ho_levels, PETSC_NULL);
125 CHKERR PetscOptionsInt("-order_lambda",
126 "approximation order of Lagrange multipliers", "", 1,
127 &order_lambda, PETSC_NULL);
128
129 CHKERR PetscOptionsInt("-step_num", "number of steps", "", nb_sub_steps,
130 &nb_sub_steps, PETSC_NULL);
131
132 CHKERR PetscOptionsBool("-is_partitioned",
133 "set if mesh is partitioned (this result that each "
134 "process keeps only part of the mes",
135 "", PETSC_FALSE, &is_partitioned, PETSC_NULL);
136
137 CHKERR PetscOptionsReal("-cn_value", "default regularisation cn value", "",
138 1., &cn_value, PETSC_NULL);
139
140 CHKERR PetscOptionsBool("-is_newton_cotes",
141 "set if Newton-Cotes quadrature rules are used", "",
142 PETSC_FALSE, &is_newton_cotes, PETSC_NULL);
143 CHKERR PetscOptionsInt("-test_num", "test number", "", 0, &test_num,
144 PETSC_NULL);
145 CHKERR PetscOptionsBool("-convect", "set to convect integration pts", "",
146 PETSC_FALSE, &convect_pts, PETSC_NULL);
147 CHKERR PetscOptionsBool("-print_contact_state",
148 "output number of active gp at every iteration", "",
149 PETSC_FALSE, &print_contact_state, PETSC_NULL);
150 CHKERR PetscOptionsBool("-alm_flag",
151 "if set use ALM, if not use C-function", "",
152 PETSC_FALSE, &alm_flag, PETSC_NULL);
153 CHKERR PetscOptionsBool("-eigen_pos_flag",
154 "if set use eigen spatial positions are taken into "
155 "account for predeformed configuration",
156 "", PETSC_FALSE, &eigen_pos_flag, PETSC_NULL);
157
158 CHKERR PetscOptionsReal("-scale_factor", "scale factor", "", scale_factor,
159 &scale_factor, PETSC_NULL);
160
161 CHKERR PetscOptionsBool("-ignore_contact", "if set true, ignore contact",
162 "", PETSC_FALSE, &ignore_contact, PETSC_NULL);
163 CHKERR PetscOptionsBool("-ignore_pressure", "if set true, ignore pressure",
164 "", PETSC_FALSE, &ignore_pressure, PETSC_NULL);
165 CHKERR PetscOptionsBool("-analytical_input",
166 "if set true, use analytical strain", "",
167 PETSC_TRUE, &analytical_input, PETSC_NULL);
168 CHKERR PetscOptionsBool("-save_mean_stress",
169 "if set true, save mean stress", "", PETSC_TRUE,
170 &save_mean_stress, PETSC_NULL);
171 CHKERR PetscOptionsString("-stress_tag_name", "stress tag name file name",
172 "", "INTERNAL_STRESS", stress_tag_name, 255,
173 &flg_tag_name);
174 CHKERR PetscOptionsReal(
175 "-thermal_expansion_coef", "thermal expansion coef ", "",
176 thermal_expansion_coef, &thermal_expansion_coef, PETSC_NULL);
177 CHKERR PetscOptionsReal("-init_temp", "init_temp ", "", init_temp,
178 &init_temp, PETSC_NULL);
179
180 ierr = PetscOptionsEnd();
181 CHKERRQ(ierr);
182
183 // Check if mesh file was provided
184 if (flg_file != PETSC_TRUE) {
185 SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -file_name (MESH FILE NEEDED)");
186 }
187
188 if (is_partitioned) {
189 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
190 "Partitioned mesh is not supported");
191 }
192
193 const char *option;
194 option = "";
195 CHKERR moab.load_file(mesh_file_name, 0, option);
196 Version file_ver;
198 MOFEM_LOG("WORLD", Sev::inform) << "File version " << file_ver.strVersion();
199
200 // Create MoFEM database and link it to MoAB
201 MoFEM::Core core(moab);
202 MoFEM::Interface &m_field = core;
203
204 std::vector<BitRefLevel> bit_levels;
205 bit_levels.push_back(BitRefLevel().set(0));
206 auto bit_ref_manager = m_field.getInterface<BitRefManager>();
207 CHKERR bit_ref_manager->setBitRefLevelByDim(0, 3, bit_levels.back());
208
209 auto add_prism_interface = [&](std::vector<BitRefLevel> &bit_levels) {
211 auto prism_interface = m_field.getInterface<PrismInterface>();
212
214 if (cit->getName().compare(0, 11, "INT_CONTACT") == 0) {
215 CHKERR PetscPrintf(PETSC_COMM_WORLD, "Insert %s (id: %d)\n",
216 cit->getName().c_str(), cit->getMeshsetId());
217 EntityHandle cubit_meshset = cit->getMeshset();
218
219 // get tet entities from back bit_level
220 EntityHandle ref_level_meshset;
221 CHKERR moab.create_meshset(MESHSET_SET, ref_level_meshset);
222 CHKERR bit_ref_manager->getEntitiesByTypeAndRefLevel(
223 bit_levels.back(), BitRefLevel().set(), MBTET, ref_level_meshset);
224 CHKERR bit_ref_manager->getEntitiesByTypeAndRefLevel(
225 bit_levels.back(), BitRefLevel().set(), MBPRISM,
226 ref_level_meshset);
227
228 // get faces and tets to split
229 CHKERR prism_interface->getSides(cubit_meshset, bit_levels.back(),
230 true, 0);
231 // set new bit level
232 bit_levels.push_back(BitRefLevel().set(bit_levels.size()));
233 // split faces and tets
234 CHKERR prism_interface->splitSides(ref_level_meshset,
235 bit_levels.back(), cubit_meshset,
236 true, true, 0);
237 // clean meshsets
238 CHKERR moab.delete_entities(&ref_level_meshset, 1);
239
240 // update cubit meshsets
241 for (_IT_CUBITMESHSETS_FOR_LOOP_(m_field, ciit)) {
242 EntityHandle cubit_meshset = ciit->meshset;
243 CHKERR bit_ref_manager->updateMeshsetByEntitiesChildren(
244 cubit_meshset, bit_levels.back(), cubit_meshset, MBVERTEX,
245 true);
246 CHKERR bit_ref_manager->updateMeshsetByEntitiesChildren(
247 cubit_meshset, bit_levels.back(), cubit_meshset, MBEDGE, true);
248 CHKERR bit_ref_manager->updateMeshsetByEntitiesChildren(
249 cubit_meshset, bit_levels.back(), cubit_meshset, MBTRI, true);
250 CHKERR bit_ref_manager->updateMeshsetByEntitiesChildren(
251 cubit_meshset, bit_levels.back(), cubit_meshset, MBTET, true);
252 }
253
254 CHKERR bit_ref_manager->shiftRightBitRef(1);
255 bit_levels.pop_back();
256 }
257 }
258
260 };
261
262 auto find_contact_prisms = [&](std::vector<BitRefLevel> &bit_levels,
263 Range &simple_contact_prisms,
264 Range &simple_master_tris,
265 Range &simple_slave_tris) {
267
268 EntityHandle meshset_prisms;
269 CHKERR moab.create_meshset(MESHSET_SET, meshset_prisms);
270 CHKERR bit_ref_manager->getEntitiesByTypeAndRefLevel(
271 bit_levels.back(), BitRefLevel().set(), MBPRISM, meshset_prisms);
272 CHKERR moab.get_entities_by_handle(meshset_prisms, simple_contact_prisms);
273 CHKERR moab.delete_entities(&meshset_prisms, 1);
274
275 EntityHandle tri;
276 for (Range::iterator pit = simple_contact_prisms.begin();
277 pit != simple_contact_prisms.end(); pit++) {
278 CHKERR moab.side_element(*pit, 2, 3, tri);
279 simple_master_tris.insert(tri);
280 CHKERR moab.side_element(*pit, 2, 4, tri);
281 simple_slave_tris.insert(tri);
282 }
283
285 };
286
287 Range simple_contact_prisms, simple_master_tris, simple_slave_tris;
288
289 if (!ignore_contact) {
290 CHKERR add_prism_interface(bit_levels);
291 CHKERR find_contact_prisms(bit_levels, simple_contact_prisms,
292 simple_master_tris, simple_slave_tris);
293 }
294
295 Range mortar_contact_prisms, mortar_master_tris, mortar_slave_tris;
296
298 if (it->getName().compare(0, 13, "MORTAR_MASTER") == 0) {
299 CHKERR m_field.get_moab().get_entities_by_type(
300 it->meshset, MBTRI, mortar_master_tris, true);
301 }
302 }
303
305 if (it->getName().compare(0, 12, "MORTAR_SLAVE") == 0) {
306 CHKERR m_field.get_moab().get_entities_by_type(it->meshset, MBTRI,
307 mortar_slave_tris, true);
308 }
309 }
310
311 ContactSearchKdTree contact_search_kd_tree(m_field);
312
313 boost::shared_ptr<ContactSearchKdTree::ContactCommonData_multiIndex>
314 contact_commondata_multi_index;
315
316 contact_commondata_multi_index =
317 boost::shared_ptr<ContactSearchKdTree::ContactCommonData_multiIndex>(
319
320 CHKERR contact_search_kd_tree.buildTree(mortar_master_tris);
321
322 CHKERR contact_search_kd_tree.contactSearchAlgorithm(
323 mortar_master_tris, mortar_slave_tris, contact_commondata_multi_index,
324 mortar_contact_prisms);
325
326 EntityHandle meshset_slave_master_prisms;
327 CHKERR moab.create_meshset(MESHSET_SET, meshset_slave_master_prisms);
328
329 CHKERR
330 moab.add_entities(meshset_slave_master_prisms, mortar_contact_prisms);
331
332 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
333 meshset_slave_master_prisms, 3, bit_levels.back());
334
335 Range tris_from_prism;
336 // find actual masters and slave used
337 CHKERR m_field.get_moab().get_adjacencies(mortar_contact_prisms, 2, true,
338 tris_from_prism,
339 moab::Interface::UNION);
340
341 tris_from_prism = tris_from_prism.subset_by_type(MBTRI);
342 mortar_slave_tris = intersect(tris_from_prism, mortar_slave_tris);
343
344 Range all_contact_prisms, all_slave_tris;
345 all_contact_prisms = unite(simple_contact_prisms, mortar_contact_prisms);
346 all_slave_tris = unite(simple_slave_tris, mortar_slave_tris);
347
348 CHKERR m_field.add_field("SPATIAL_POSITION", H1, AINSWORTH_LEGENDRE_BASE, 3,
349 MB_TAG_SPARSE, MF_ZERO);
350
351 CHKERR m_field.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
352 3, MB_TAG_SPARSE, MF_ZERO);
353
354 if (!eigen_pos_flag) {
355 CHKERR m_field.add_field("EIGEN_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
356 3, MB_TAG_SPARSE, MF_ZERO);
357 }
358
359 CHKERR m_field.add_field("LAGMULT", H1, AINSWORTH_LEGENDRE_BASE, 1,
360 MB_TAG_SPARSE, MF_ZERO);
361
362 CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "MESH_NODE_POSITIONS");
363 CHKERR m_field.set_field_order(0, MBTET, "MESH_NODE_POSITIONS", 1);
364 CHKERR m_field.set_field_order(0, MBTRI, "MESH_NODE_POSITIONS", 1);
365 CHKERR m_field.set_field_order(0, MBEDGE, "MESH_NODE_POSITIONS", 1);
366 CHKERR m_field.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
367
368 // Declare problem add entities (by tets) to the field
369 CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "SPATIAL_POSITION");
370 CHKERR m_field.set_field_order(0, MBTET, "SPATIAL_POSITION", order);
371 CHKERR m_field.set_field_order(0, MBTRI, "SPATIAL_POSITION", order);
372 CHKERR m_field.set_field_order(0, MBEDGE, "SPATIAL_POSITION", order);
373 CHKERR m_field.set_field_order(0, MBVERTEX, "SPATIAL_POSITION", 1);
374
375 if (!eigen_pos_flag) {
376 CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "EIGEN_POSITIONS");
377 CHKERR m_field.set_field_order(0, MBTET, "EIGEN_POSITIONS", order);
378 CHKERR m_field.set_field_order(0, MBTRI, "EIGEN_POSITIONS", order);
379 CHKERR m_field.set_field_order(0, MBEDGE, "EIGEN_POSITIONS", order);
380 CHKERR m_field.set_field_order(0, MBVERTEX, "EIGEN_POSITIONS", 1);
381 }
382
383 CHKERR m_field.add_ents_to_field_by_type(all_slave_tris, MBTRI, "LAGMULT");
384
385 CHKERR m_field.set_field_order(0, MBTRI, "LAGMULT", order_lambda);
386 CHKERR m_field.set_field_order(0, MBEDGE, "LAGMULT", order_lambda);
387 CHKERR m_field.set_field_order(0, MBVERTEX, "LAGMULT", 1);
388
389 auto set_contact_order = [&](Range &contact_prisms, int order_contact,
390 int nb_ho_levels) {
392 Range contact_tris, contact_edges;
393 CHKERR moab.get_adjacencies(contact_prisms, 2, false, contact_tris,
394 moab::Interface::UNION);
395 contact_tris = contact_tris.subset_by_type(MBTRI);
396 CHKERR moab.get_adjacencies(contact_tris, 1, false, contact_edges,
397 moab::Interface::UNION);
398 Range ho_ents;
399 ho_ents.merge(contact_tris);
400 ho_ents.merge(contact_edges);
401 for (int ll = 0; ll < nb_ho_levels; ll++) {
402 Range ents, verts, tets;
403 CHKERR moab.get_connectivity(ho_ents, verts, true);
404 CHKERR moab.get_adjacencies(verts, 3, false, tets,
405 moab::Interface::UNION);
406 tets = tets.subset_by_type(MBTET);
407 for (auto d : {1, 2}) {
408 CHKERR moab.get_adjacencies(tets, d, false, ents,
409 moab::Interface::UNION);
410 }
411 ho_ents = unite(ho_ents, ents);
412 ho_ents = unite(ho_ents, tets);
413 }
414
415 CHKERR m_field.set_field_order(ho_ents, "SPATIAL_POSITION",
416 order_contact);
417
419 };
420
421 if (!ignore_contact && order_contact > order) {
422 CHKERR set_contact_order(all_contact_prisms, order_contact, nb_ho_levels);
423 }
424
425 // build field
426 CHKERR m_field.build_fields();
427
428 // Projection on "x" field
429 {
430 Projection10NodeCoordsOnField ent_method(m_field, "SPATIAL_POSITION");
431 CHKERR m_field.loop_dofs("SPATIAL_POSITION", ent_method);
432 }
433 // Projection on "X" field
434 {
435 Projection10NodeCoordsOnField ent_method(m_field, "MESH_NODE_POSITIONS");
436 CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method);
437 }
438
439 Range slave_verts;
440 CHKERR moab.get_connectivity(all_slave_tris, slave_verts, true);
441 CHKERR m_field.getInterface<FieldBlas>()->setField(0.0, MBVERTEX,
442 slave_verts, "LAGMULT");
443
444 // Add elastic element
445 boost::shared_ptr<std::map<int, BlockData>> block_sets_ptr =
446 boost::make_shared<std::map<int, BlockData>>();
447 CHKERR HookeElement::setBlocks(m_field, block_sets_ptr);
448
449 boost::shared_ptr<ForcesAndSourcesCore> fe_elastic_lhs_ptr(
451 boost::shared_ptr<ForcesAndSourcesCore> fe_elastic_rhs_ptr(
453 fe_elastic_lhs_ptr->getRuleHook = VolRule();
454 fe_elastic_rhs_ptr->getRuleHook = VolRule();
455
456 CHKERR HookeElement::addElasticElement(m_field, block_sets_ptr, "ELASTIC",
457 "SPATIAL_POSITION",
458 "MESH_NODE_POSITIONS", false);
459
460 auto data_hooke_element_at_pts =
461 boost::make_shared<HookeInternalStressElement::DataAtIntegrationPts>();
462 CHKERR HookeElement::setOperators(fe_elastic_lhs_ptr, fe_elastic_rhs_ptr,
463 block_sets_ptr, "SPATIAL_POSITION",
464 "MESH_NODE_POSITIONS", false, false,
465 MBTET, data_hooke_element_at_pts);
466 auto thermal_strain =
467 [&thermal_expansion_coef, &init_temp, &test_num](
472
474 double temp;
475
476 t_thermal_strain(i, j) = 0.0;
477
478 switch (test_num) {
479 case 0:
480 // Put here analytical formula which may depend on coordinates
481 temp = init_temp - 1.0;
482 t_thermal_strain(i, j) =
483 thermal_expansion_coef * (temp - init_temp) * t_kd(i, j);
484 break;
485 case 1:
486 case 2:
487 t_thermal_strain(i, j) = -thermal_expansion_coef * t_kd(i, j);
488 break;
489 case 3:
490 t_thermal_strain(2, 2) = thermal_expansion_coef;
491 break;
492 case 4:
493 t_thermal_strain(i, j) = thermal_expansion_coef * t_kd(i, j);
494 break;
495 default:
496 break;
497 }
498
499 return t_thermal_strain;
500 };
501
502 if (analytical_input) {
503 fe_elastic_rhs_ptr->getOpPtrVector().push_back(
504 new HookeElement::OpAnalyticalInternalStrain_dx<0>(
505 "SPATIAL_POSITION", data_hooke_element_at_pts, thermal_strain));
506 fe_elastic_rhs_ptr->getOpPtrVector().push_back(
508 "SPATIAL_POSITION", "SPATIAL_POSITION", data_hooke_element_at_pts,
509 thermal_strain));
510 } else {
511 fe_elastic_rhs_ptr->getOpPtrVector().push_back(
513 "SPATIAL_POSITION", "SPATIAL_POSITION", data_hooke_element_at_pts,
514 moab, stress_tag_name));
515 fe_elastic_rhs_ptr->getOpPtrVector().push_back(
517 "SPATIAL_POSITION", data_hooke_element_at_pts));
518 }
519
520 fe_elastic_rhs_ptr->getOpPtrVector().push_back(
522 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat));
523 fe_elastic_rhs_ptr->getOpPtrVector().push_back(
525 "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
526 fe_elastic_rhs_ptr->getOpPtrVector().push_back(
528 "SPATIAL_POSITION", "SPATIAL_POSITION", data_hooke_element_at_pts,
529 *block_sets_ptr.get(), moab, scale_factor, save_mean_stress, false,
530 false));
531
532 Range all_tets;
534 if (bit->getName().compare(0, 11, "MAT_ELASTIC") == 0) {
535 Range tets;
536 CHKERR moab.get_entities_by_dimension(bit->getMeshset(), 3, tets, true);
537 all_tets.merge(tets);
538 }
539 }
540 Skinner skinner(&moab);
541 Range skin_tris;
542 CHKERR skinner.find_skin(0, all_tets, false, skin_tris);
543
544 CHKERR m_field.add_finite_element("SKIN", MF_ZERO);
545 CHKERR m_field.add_ents_to_finite_element_by_type(skin_tris, MBTRI, "SKIN");
547 "SPATIAL_POSITION");
549 "SPATIAL_POSITION");
551 "SPATIAL_POSITION");
553 "MESH_NODE_POSITIONS");
555 "EIGEN_POSITIONS");
556
557 auto make_mortar_contact_element = [&]() {
558 return boost::make_shared<MortarContactProblem::MortarContactElement>(
559 m_field, contact_commondata_multi_index, "SPATIAL_POSITION",
560 "MESH_NODE_POSITIONS");
561 };
562
563 auto make_convective_mortar_master_element = [&]() {
564 return boost::make_shared<
566 m_field, contact_commondata_multi_index, "SPATIAL_POSITION",
567 "MESH_NODE_POSITIONS");
568 };
569
570 auto make_convective_mortar_slave_element = [&]() {
571 return boost::make_shared<
573 m_field, contact_commondata_multi_index, "SPATIAL_POSITION",
574 "MESH_NODE_POSITIONS");
575 };
576
577 auto make_mortar_contact_common_data = [&]() {
578 return boost::make_shared<MortarContactProblem::CommonDataMortarContact>(
579 m_field);
580 };
581
582 auto get_mortar_contact_rhs = [&](auto contact_problem, auto make_element,
583 bool is_alm = false) {
584 auto fe_rhs_mortar_contact = make_element();
585 auto common_data_mortar_contact = make_mortar_contact_common_data();
586 if (print_contact_state) {
587 fe_rhs_mortar_contact->contactStateVec =
588 common_data_mortar_contact->gaussPtsStateVec;
589 }
590 contact_problem->setContactOperatorsRhs(
591 fe_rhs_mortar_contact, common_data_mortar_contact, "SPATIAL_POSITION",
592 "LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS");
593 return fe_rhs_mortar_contact;
594 };
595
596 auto get_mortar_master_traction_rhs = [&](auto contact_problem,
597 auto make_element,
598 bool is_alm = false) {
599 auto fe_rhs_mortar_contact = make_element();
600 auto common_data_mortar_contact = make_mortar_contact_common_data();
601 contact_problem->setMasterForceOperatorsRhs(
602 fe_rhs_mortar_contact, common_data_mortar_contact, "SPATIAL_POSITION",
603 "LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS");
604 return fe_rhs_mortar_contact;
605 };
606
607 auto get_mortar_master_traction_lhs = [&](auto contact_problem,
608 auto make_element,
609 bool is_alm = false) {
610 auto fe_lhs_mortar_contact = make_element();
611 auto common_data_mortar_contact = make_mortar_contact_common_data();
612 contact_problem->setMasterForceOperatorsLhs(
613 fe_lhs_mortar_contact, common_data_mortar_contact, "SPATIAL_POSITION",
614 "LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS");
615 return fe_lhs_mortar_contact;
616 };
617
618 auto get_mortar_master_help_traction_lhs =
619 [&](auto contact_problem, auto make_element, bool is_alm = false) {
620 auto fe_lhs_mortar_contact = make_element();
621 auto common_data_mortar_contact = make_mortar_contact_common_data();
622 contact_problem->setMasterForceOperatorsLhs(
623 fe_lhs_mortar_contact, common_data_mortar_contact,
624 "SPATIAL_POSITION", "LAGMULT", is_alm);
625 return fe_lhs_mortar_contact;
626 };
627
628 auto get_mortar_contact_lhs = [&](auto contact_problem, auto make_element,
629 bool is_alm = false) {
630 auto fe_lhs_mortar_contact = make_element();
631 auto common_data_mortar_contact = make_mortar_contact_common_data();
632 contact_problem->setContactOperatorsLhs(
633 fe_lhs_mortar_contact, common_data_mortar_contact, "SPATIAL_POSITION",
634 "LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS");
635 return fe_lhs_mortar_contact;
636 };
637
638 auto get_mortar_contact_help_lhs =
639 [&](auto contact_problem, auto make_element, bool is_alm = false) {
640 auto fe_lhs_mortar_contact = make_element();
641 auto common_data_mortar_contact = make_mortar_contact_common_data();
642 contact_problem->setContactOperatorsLhs(
643 fe_lhs_mortar_contact, common_data_mortar_contact,
644 "SPATIAL_POSITION", "LAGMULT", is_alm);
645 return fe_lhs_mortar_contact;
646 };
647
648 auto cn_value_ptr = boost::make_shared<double>(cn_value);
649 auto mortar_contact_problem = boost::make_shared<MortarContactProblem>(
650 m_field, contact_commondata_multi_index, cn_value_ptr, is_newton_cotes);
651
652 // add fields to the global matrix by adding the element
653
654 if (!eigen_pos_flag)
655 mortar_contact_problem->addMortarContactElement(
656 "MORTAR_CONTACT_ELEM", "SPATIAL_POSITION", "LAGMULT",
657 mortar_contact_prisms);
658 else
659 mortar_contact_problem->addMortarContactElement(
660 "MORTAR_CONTACT_ELEM", "SPATIAL_POSITION", "LAGMULT",
661 mortar_contact_prisms, "MESH_NODE_POSITIONS", eigen_pos_flag,
662 "EIGEN_POSITIONS");
663
664 auto make_simple_contact_element = [&]() {
665 return boost::make_shared<SimpleContactProblem::SimpleContactElement>(
666 m_field);
667 };
668
669 auto make_convective_simple_master_element = [&]() {
670 return boost::make_shared<
672 m_field, "SPATIAL_POSITION", "MESH_NODE_POSITIONS");
673 };
674
675 auto make_convective_simple_slave_element = [&]() {
676 return boost::make_shared<
678 m_field, "SPATIAL_POSITION", "MESH_NODE_POSITIONS");
679 };
680
681 auto make_simple_contact_common_data = [&]() {
682 return boost::make_shared<SimpleContactProblem::CommonDataSimpleContact>(
683 m_field);
684 };
685
686 auto get_simple_contact_rhs = [&](auto contact_problem, auto make_element,
687 bool is_alm = false) {
688 auto fe_rhs_simple_contact = make_element();
689 auto common_data_simple_contact = make_simple_contact_common_data();
690 if (print_contact_state) {
691 fe_rhs_simple_contact->contactStateVec =
692 common_data_simple_contact->gaussPtsStateVec;
693 }
694 contact_problem->setContactOperatorsRhs(
695 fe_rhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
696 "LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS");
697 return fe_rhs_simple_contact;
698 };
699
700 auto get_simple_master_traction_rhs = [&](auto contact_problem,
701 auto make_element,
702 bool is_alm = false) {
703 auto fe_rhs_simple_contact = make_element();
704 auto common_data_simple_contact = make_simple_contact_common_data();
705 contact_problem->setMasterForceOperatorsRhs(
706 fe_rhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
707 "LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS");
708 return fe_rhs_simple_contact;
709 };
710
711 auto get_simple_master_traction_lhs =
712 [&](auto contact_problem, auto make_element, bool is_alm = false) {
713 auto fe_lhs_simple_contact = make_element();
714 auto common_data_simple_contact = make_simple_contact_common_data();
715 contact_problem->setMasterForceOperatorsLhs(
716 fe_lhs_simple_contact, common_data_simple_contact,
717 "SPATIAL_POSITION", "LAGMULT", is_alm);
718 return fe_lhs_simple_contact;
719 };
720
721 auto get_contact_lhs = [&](auto contact_problem, auto make_element,
722 bool is_alm = false) {
723 auto fe_lhs_simple_contact = make_element();
724 auto common_data_simple_contact = make_simple_contact_common_data();
725 contact_problem->setContactOperatorsLhs(
726 fe_lhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
727 "LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS");
728 return fe_lhs_simple_contact;
729 };
730
731 auto simple_contact_problem = boost::make_shared<SimpleContactProblem>(
732 m_field, cn_value_ptr, is_newton_cotes);
733 simple_contact_problem->addContactElement("SIMPLE_CONTACT_ELEM",
734 "SPATIAL_POSITION", "LAGMULT",
735 simple_contact_prisms);
736
737 simple_contact_problem->addPostProcContactElement(
738 "CONTACT_POST_PROC", "SPATIAL_POSITION", "LAGMULT",
739 "MESH_NODE_POSITIONS", all_slave_tris);
740
741 CHKERR MetaNeumannForces::addNeumannBCElements(m_field, "SPATIAL_POSITION");
742
743 // Add spring boundary condition applied on surfaces.
744 CHKERR MetaSpringBC::addSpringElements(m_field, "SPATIAL_POSITION",
745 "MESH_NODE_POSITIONS");
746
747 // build finite elemnts
749
750 // build adjacencies
751 CHKERR m_field.build_adjacencies(bit_levels.back());
752
753 // define problems
754 CHKERR m_field.add_problem("CONTACT_PROB", MF_ZERO);
755
756 // set refinement level for problem
757 CHKERR m_field.modify_problem_ref_level_add_bit("CONTACT_PROB",
758 bit_levels.back());
759
760 DMType dm_name = "DMMOFEM";
761 CHKERR DMRegister_MoFEM(dm_name);
762
764 dm = createSmartDM(m_field.get_comm(), dm_name);
765
766 // create dm instance
767 CHKERR DMSetType(dm, dm_name);
768
769 // set dm datastruture which created mofem datastructures
770 CHKERR DMMoFEMCreateMoFEM(dm, &m_field, "CONTACT_PROB", bit_levels.back());
771 CHKERR DMSetFromOptions(dm);
772 CHKERR DMMoFEMSetIsPartitioned(dm, is_partitioned);
773 // add elements to dm
774 CHKERR DMMoFEMAddElement(dm, "MORTAR_CONTACT_ELEM");
775 CHKERR DMMoFEMAddElement(dm, "SIMPLE_CONTACT_ELEM");
776 CHKERR DMMoFEMAddElement(dm, "ELASTIC");
777 CHKERR DMMoFEMAddElement(dm, "PRESSURE_FE");
778 CHKERR DMMoFEMAddElement(dm, "SPRING");
779 CHKERR DMMoFEMAddElement(dm, "CONTACT_POST_PROC");
780 CHKERR DMMoFEMAddElement(dm, "SKIN");
781
782 CHKERR DMSetUp(dm);
783
784 // Vector of DOFs and the RHS
785 auto D = smartCreateDMVector(dm);
786 auto F = smartVectorDuplicate(D);
787
788 // Stiffness matrix
789 auto Aij = smartCreateDMMatrix(dm);
790
791 CHKERR VecZeroEntries(D);
792 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_FORWARD);
793 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
794 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
795
796 CHKERR VecZeroEntries(F);
797 CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
798 CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
799
800 CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
801 CHKERR MatZeroEntries(Aij);
802
803 fe_elastic_rhs_ptr->snes_f = F;
804 fe_elastic_lhs_ptr->snes_B = Aij;
805
806 // Dirichlet BC
807 boost::shared_ptr<FEMethod> dirichlet_bc_ptr =
808 boost::shared_ptr<FEMethod>(new DirichletSpatialPositionsBc(
809 m_field, "SPATIAL_POSITION", Aij, D, F));
810
811 dirichlet_bc_ptr->snes_ctx = SnesMethod::CTX_SNESNONE;
812 dirichlet_bc_ptr->snes_x = D;
813
814 // Assemble pressure and traction forces
815 boost::ptr_map<std::string, NeumannForcesSurface> neumann_forces;
817 m_field, neumann_forces, NULL, "SPATIAL_POSITION");
818
819 boost::ptr_map<std::string, NeumannForcesSurface>::iterator mit =
820 neumann_forces.begin();
821 for (; mit != neumann_forces.end(); mit++) {
822 mit->second->methodsOp.push_back(new MortarContactProblem::LoadScale());
823 CHKERR DMMoFEMSNESSetFunction(dm, mit->first.c_str(),
824 &mit->second->getLoopFe(), NULL, NULL);
825 }
826
827 // Implementation of spring element
828 // Create new instances of face elements for springs
829 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr(
831 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr(
833
835 m_field, fe_spring_lhs_ptr, fe_spring_rhs_ptr, "SPATIAL_POSITION",
836 "MESH_NODE_POSITIONS");
837
838 CHKERR DMoFEMPreProcessFiniteElements(dm, dirichlet_bc_ptr.get());
839 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
840 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
841 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
843 dirichlet_bc_ptr.get(), NULL);
844 if (convect_pts == PETSC_TRUE) {
846 dm, "MORTAR_CONTACT_ELEM",
847 get_mortar_contact_rhs(mortar_contact_problem,
848 make_convective_mortar_master_element),
849 PETSC_NULL, PETSC_NULL);
851 dm, "MORTAR_CONTACT_ELEM",
852 get_mortar_master_traction_rhs(mortar_contact_problem,
853 make_convective_mortar_slave_element),
854 PETSC_NULL, PETSC_NULL);
855 } else {
857 dm, "MORTAR_CONTACT_ELEM",
858 get_mortar_contact_rhs(mortar_contact_problem,
859 make_mortar_contact_element, alm_flag),
860 PETSC_NULL, PETSC_NULL);
862 dm, "MORTAR_CONTACT_ELEM",
863 get_mortar_master_traction_rhs(mortar_contact_problem,
864 make_mortar_contact_element, alm_flag),
865 PETSC_NULL, PETSC_NULL);
866 }
867
868 if (convect_pts == PETSC_TRUE) {
870 dm, "SIMPLE_CONTACT_ELEM",
871 get_simple_contact_rhs(simple_contact_problem,
872 make_convective_simple_master_element),
873 PETSC_NULL, PETSC_NULL);
875 dm, "SIMPLE_CONTACT_ELEM",
876 get_simple_master_traction_rhs(simple_contact_problem,
877 make_convective_simple_slave_element),
878 PETSC_NULL, PETSC_NULL);
879 } else {
881 dm, "SIMPLE_CONTACT_ELEM",
882 get_simple_contact_rhs(simple_contact_problem,
883 make_simple_contact_element, alm_flag),
884 PETSC_NULL, PETSC_NULL);
886 dm, "SIMPLE_CONTACT_ELEM",
887 get_simple_master_traction_rhs(simple_contact_problem,
888 make_simple_contact_element, alm_flag),
889 PETSC_NULL, PETSC_NULL);
890 }
891
892 CHKERR DMMoFEMSNESSetFunction(dm, "ELASTIC", fe_elastic_rhs_ptr, PETSC_NULL,
893 PETSC_NULL);
894
895 CHKERR DMMoFEMSNESSetFunction(dm, "SPRING", fe_spring_rhs_ptr, PETSC_NULL,
896 PETSC_NULL);
897 CHKERR DMMoFEMSNESSetFunction(dm, DM_NO_ELEMENT, PETSC_NULL, PETSC_NULL,
898 dirichlet_bc_ptr.get());
899
900 boost::shared_ptr<FEMethod> fe_null;
901 CHKERR DMMoFEMSNESSetJacobian(dm, DM_NO_ELEMENT, fe_null, dirichlet_bc_ptr,
902 fe_null);
903 if (convect_pts == PETSC_TRUE) {
905 dm, "MORTAR_CONTACT_ELEM",
906 get_mortar_contact_help_lhs(mortar_contact_problem,
907 make_convective_mortar_master_element),
908 PETSC_NULL, PETSC_NULL);
910 dm, "MORTAR_CONTACT_ELEM",
911 get_mortar_master_help_traction_lhs(
912 mortar_contact_problem, make_convective_mortar_slave_element),
913 PETSC_NULL, PETSC_NULL);
914 } else {
916 dm, "MORTAR_CONTACT_ELEM",
917 get_mortar_contact_lhs(mortar_contact_problem,
918 make_mortar_contact_element, alm_flag),
919 PETSC_NULL, PETSC_NULL);
921 dm, "MORTAR_CONTACT_ELEM",
922 get_mortar_master_traction_lhs(mortar_contact_problem,
923 make_mortar_contact_element, alm_flag),
924 PETSC_NULL, PETSC_NULL);
925 }
926
927 if (convect_pts == PETSC_TRUE) {
929 dm, "SIMPLE_CONTACT_ELEM",
930 get_contact_lhs(simple_contact_problem,
931 make_convective_simple_master_element),
932 PETSC_NULL, PETSC_NULL);
934 dm, "SIMPLE_CONTACT_ELEM",
935 get_simple_master_traction_lhs(simple_contact_problem,
936 make_convective_simple_slave_element),
937 PETSC_NULL, PETSC_NULL);
938 } else {
939 CHKERR DMMoFEMSNESSetJacobian(dm, "SIMPLE_CONTACT_ELEM",
940 get_contact_lhs(simple_contact_problem,
941 make_simple_contact_element,
942 alm_flag),
943 PETSC_NULL, PETSC_NULL);
945 dm, "SIMPLE_CONTACT_ELEM",
946 get_simple_master_traction_lhs(simple_contact_problem,
947 make_simple_contact_element, alm_flag),
948 PETSC_NULL, PETSC_NULL);
949 }
950 CHKERR DMMoFEMSNESSetJacobian(dm, "ELASTIC", fe_elastic_lhs_ptr, PETSC_NULL,
951 PETSC_NULL);
952 CHKERR DMMoFEMSNESSetJacobian(dm, "SPRING", fe_spring_lhs_ptr, PETSC_NULL,
953 PETSC_NULL);
954 CHKERR DMMoFEMSNESSetJacobian(dm, DM_NO_ELEMENT, fe_null, fe_null,
955 dirichlet_bc_ptr);
956
957 if (test_num) {
958 char testing_options[] = "-ksp_type fgmres "
959 "-pc_type lu "
960 "-pc_factor_mat_solver_type mumps "
961 "-snes_type newtonls "
962 "-snes_linesearch_type basic "
963 "-snes_max_it 10 "
964 "-snes_atol 1e-8 "
965 "-snes_rtol 1e-8 ";
966 CHKERR PetscOptionsInsertString(PETSC_NULL, testing_options);
967 }
968
969 auto snes = MoFEM::createSNES(m_field.get_comm());
970 CHKERR SNESSetDM(snes, dm);
971 SnesCtx *snes_ctx;
972 // create snes nonlinear solver
973 {
974 CHKERR SNESSetDM(snes, dm);
975 CHKERR DMMoFEMGetSnesCtx(dm, &snes_ctx);
976 CHKERR SNESSetFunction(snes, F, SnesRhs, snes_ctx);
977 CHKERR SNESSetJacobian(snes, Aij, Aij, SnesMat, snes_ctx);
978 CHKERR SNESSetFromOptions(snes);
979 }
980
981 /// Post proc on the skin
982 PostProcFaceOnRefinedMesh post_proc_skin(m_field);
983 CHKERR post_proc_skin.generateReferenceElementMesh();
984 CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", post_proc_skin, false, false);
985 CHKERR post_proc_skin.addFieldValuesPostProc("SPATIAL_POSITION");
986 CHKERR post_proc_skin.addFieldValuesPostProc("MESH_NODE_POSITIONS");
987 CHKERR post_proc_skin.addFieldValuesPostProc("EIGEN_POSITIONS");
988
989 struct OpGetFieldGradientValuesOnSkin
990 : public FaceElementForcesAndSourcesCore::UserDataOperator {
991
992 const std::string feVolName;
993 boost::shared_ptr<VolSideFe> sideOpFe;
994
995 OpGetFieldGradientValuesOnSkin(const std::string field_name,
996 const std::string vol_fe_name,
997 boost::shared_ptr<VolSideFe> side_fe)
1000 feVolName(vol_fe_name), sideOpFe(side_fe) {}
1001
1002 MoFEMErrorCode doWork(int side, EntityType type,
1005 if (type != MBVERTEX)
1007 CHKERR loopSideVolumes(feVolName, *sideOpFe);
1009 }
1010 };
1011
1012 boost::shared_ptr<VolSideFe> vol_side_fe_ptr =
1013 boost::make_shared<VolSideFe>(m_field);
1014 vol_side_fe_ptr->getOpPtrVector().push_back(
1016 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat));
1017 vol_side_fe_ptr->getOpPtrVector().push_back(
1019 "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
1020
1021 post_proc_skin.getOpPtrVector().push_back(
1022 new OpGetFieldGradientValuesOnSkin("SPATIAL_POSITION", "ELASTIC",
1023 vol_side_fe_ptr));
1024 post_proc_skin.getOpPtrVector().push_back(
1026 "SPATIAL_POSITION", data_hooke_element_at_pts->spatPosMat));
1027 post_proc_skin.getOpPtrVector().push_back(
1029 "MESH_NODE_POSITIONS", data_hooke_element_at_pts->meshNodePosMat));
1030
1031 post_proc_skin.getOpPtrVector().push_back(
1034 "SPATIAL_POSITION", data_hooke_element_at_pts,
1035 *block_sets_ptr.get(), post_proc_skin.postProcMesh,
1036 post_proc_skin.mapGaussPts, false, false));
1037
1038 for (int ss = 0; ss != nb_sub_steps; ++ss) {
1039 if (!ignore_pressure) {
1040 MortarContactProblem::LoadScale::lAmbda = (ss + 1.0) / nb_sub_steps;
1041 } else {
1043 CHKERR PetscPrintf(PETSC_COMM_WORLD, "Ignoring pressure...\n");
1044 }
1045 CHKERR PetscPrintf(PETSC_COMM_WORLD, "Load scale: %6.4e\n",
1047
1048 CHKERR SNESSolve(snes, PETSC_NULL, D);
1049
1050 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
1051 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
1052 }
1053
1054 // save on mesh
1055 CHKERR DMoFEMMeshToGlobalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
1056
1057 SmartPetscObj<Vec> v_energy;
1058 CHKERR HookeElement::calculateEnergy(dm, block_sets_ptr, "SPATIAL_POSITION",
1059 "MESH_NODE_POSITIONS", false, false,
1060 v_energy);
1061 double energy;
1062 CHKERR VecSum(v_energy, &energy);
1063 // Print elastic energy
1064 MOFEM_LOG_C("WORLD", Sev::inform, "Elastic energy: %8.8e", energy);
1065
1066 // moab_instance
1067 moab::Core mb_post; // create database
1068 moab::Interface &moab_proc = mb_post; // create interface to database
1069
1070 auto common_data_mortar_contact = make_mortar_contact_common_data();
1071
1072 boost::shared_ptr<MortarContactProblem::MortarContactElement>
1073 fe_post_proc_mortar_contact;
1074 if (convect_pts == PETSC_TRUE) {
1075 fe_post_proc_mortar_contact = make_convective_mortar_master_element();
1076 } else {
1077 fe_post_proc_mortar_contact = make_mortar_contact_element();
1078 }
1079
1080 std::array<double, 2> nb_gauss_pts;
1081 std::array<double, 2> contact_area;
1082
1083 auto get_contact_data = [&](auto vec, std::array<double, 2> &data) {
1085 CHKERR VecAssemblyBegin(vec);
1086 CHKERR VecAssemblyEnd(vec);
1087 const double *array;
1088 CHKERR VecGetArrayRead(vec, &array);
1089 if (m_field.get_comm_rank() == 0) {
1090 for (int i : {0, 1})
1091 data[i] = array[i];
1092 }
1093 CHKERR VecRestoreArrayRead(vec, &array);
1095 };
1096
1097 mb_post.delete_mesh();
1098
1099 if (!mortar_contact_prisms.empty()) {
1100 mortar_contact_problem->setContactOperatorsForPostProc(
1101 fe_post_proc_mortar_contact, common_data_mortar_contact, m_field,
1102 "SPATIAL_POSITION", "LAGMULT", mb_post, alm_flag, eigen_pos_flag,
1103 "EIGEN_POSITIONS");
1104
1105 CHKERR VecZeroEntries(common_data_mortar_contact->gaussPtsStateVec);
1106 CHKERR VecZeroEntries(common_data_mortar_contact->contactAreaVec);
1107
1108 CHKERR DMoFEMLoopFiniteElements(dm, "MORTAR_CONTACT_ELEM",
1109 fe_post_proc_mortar_contact);
1110
1111 CHKERR get_contact_data(common_data_mortar_contact->gaussPtsStateVec,
1112 nb_gauss_pts);
1113 CHKERR get_contact_data(common_data_mortar_contact->contactAreaVec,
1114 contact_area);
1115
1116 if (m_field.get_comm_rank() == 0) {
1117 PetscPrintf(PETSC_COMM_SELF,
1118 "Active gauss pts (mortar): %d out of %d\n",
1119 (int)nb_gauss_pts[0], (int)nb_gauss_pts[1]);
1120
1121 PetscPrintf(
1122 PETSC_COMM_SELF,
1123 "Active contact area (mortar): %8.8f out of %8.8f (%8.8f% %)\n",
1124 contact_area[0], contact_area[1],
1125 contact_area[0] / contact_area[1] * 100.);
1126 }
1127 }
1128
1129 auto common_data_simple_contact = make_simple_contact_common_data();
1130
1131 boost::shared_ptr<SimpleContactProblem::SimpleContactElement>
1132 fe_post_proc_simple_contact;
1133 if (convect_pts == PETSC_TRUE) {
1134 fe_post_proc_simple_contact = make_convective_simple_master_element();
1135 } else {
1136 fe_post_proc_simple_contact = make_simple_contact_element();
1137 }
1138
1139 if (!simple_contact_prisms.empty()) {
1140 simple_contact_problem->setContactOperatorsForPostProc(
1141 fe_post_proc_simple_contact, common_data_simple_contact, m_field,
1142 "SPATIAL_POSITION", "LAGMULT", mb_post, alm_flag, eigen_pos_flag,
1143 "EIGEN_POSITIONS");
1144
1145 CHKERR VecZeroEntries(common_data_simple_contact->gaussPtsStateVec);
1146 CHKERR VecZeroEntries(common_data_simple_contact->contactAreaVec);
1147
1148 CHKERR DMoFEMLoopFiniteElements(dm, "SIMPLE_CONTACT_ELEM",
1149 fe_post_proc_simple_contact);
1150
1151 CHKERR get_contact_data(common_data_simple_contact->gaussPtsStateVec,
1152 nb_gauss_pts);
1153 CHKERR get_contact_data(common_data_simple_contact->contactAreaVec,
1154 contact_area);
1155
1156 if (m_field.get_comm_rank() == 0) {
1157 MOFEM_LOG_C("SELF", Sev::inform,
1158 "Active gauss pts (matching): %d out of %d",
1159 (int)nb_gauss_pts[0], (int)nb_gauss_pts[1]);
1161 "SELF", Sev::inform,
1162 "Active contact area (matching): %8.8f out of %8.8f (%8.8f% %)",
1163 contact_area[0], contact_area[1],
1164 contact_area[0] / contact_area[1] * 100.);
1165 }
1166 }
1167
1168 if (!ignore_contact) {
1169 string out_file_name;
1170 std::ostringstream strm;
1171 strm << "out_contact_integ_pts"
1172 << ".h5m";
1173 out_file_name = strm.str();
1174 MOFEM_LOG_C("WORLD", Sev::inform, "Write file %s\n",
1175 out_file_name.c_str());
1176 CHKERR mb_post.write_file(out_file_name.c_str(), "MOAB",
1177 "PARALLEL=WRITE_PART");
1178 }
1179
1180 boost::shared_ptr<PostProcFaceOnRefinedMesh> post_proc_contact_ptr(
1181 new PostProcFaceOnRefinedMesh(m_field));
1182
1183 CHKERR post_proc_contact_ptr->generateReferenceElementMesh();
1184 CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", *post_proc_contact_ptr, false,
1185 false);
1186 CHKERR post_proc_contact_ptr->addFieldValuesPostProc("LAGMULT");
1187 CHKERR post_proc_contact_ptr->addFieldValuesPostProc("SPATIAL_POSITION");
1188 CHKERR post_proc_contact_ptr->addFieldValuesPostProc("MESH_NODE_POSITIONS");
1189
1190 if (!all_slave_tris.empty()) {
1191 CHKERR DMoFEMLoopFiniteElements(dm, "CONTACT_POST_PROC",
1192 post_proc_contact_ptr);
1193 string out_file_name;
1194 std::ostringstream stm;
1195 stm << "out_contact_surface"
1196 << ".h5m";
1197 out_file_name = stm.str();
1198 MOFEM_LOG_C("WORLD", Sev::inform, "Write file %s",
1199 out_file_name.c_str());
1200 CHKERR post_proc_contact_ptr->postProcMesh.write_file(
1201 out_file_name.c_str(), "MOAB", "PARALLEL=WRITE_PART");
1202 }
1203
1204 CHKERR m_field.getInterface<FieldBlas>()->fieldCopy(1., "SPATIAL_POSITION",
1205 "EIGEN_POSITIONS");
1206
1207 {
1208 PetscPrintf(PETSC_COMM_WORLD, "Loop post proc on the skin\n");
1209 CHKERR DMoFEMLoopFiniteElements(dm, "SKIN", &post_proc_skin);
1210 ostringstream stm;
1211 string out_file_name;
1212 stm << "out_skin.h5m";
1213 out_file_name = stm.str();
1214 MOFEM_LOG_C("WORLD", Sev::inform, "Write file %s",
1215 out_file_name.c_str());
1216 CHKERR post_proc_skin.writeFile(stm.str());
1217 }
1218
1219 const int n_parts = m_field.get_comm_size();
1220 if (m_field.get_comm_rank() == 0) {
1222 dm, "ELASTIC", fe_elastic_rhs_ptr, 0, n_parts);
1223
1224 Range faces, tris, quads, tris_edges, quads_edges, ents_to_delete;
1225
1226 CHKERR moab.get_adjacencies(all_contact_prisms, 2, true, faces,
1227 moab::Interface::UNION);
1228 tris = faces.subset_by_type(MBTRI);
1229 quads = faces.subset_by_type(MBQUAD);
1230 CHKERR moab.get_adjacencies(tris, 1, true, tris_edges,
1231 moab::Interface::UNION);
1232 CHKERR moab.get_adjacencies(quads, 1, true, quads_edges,
1233 moab::Interface::UNION);
1234
1235 ents_to_delete.merge(all_contact_prisms);
1236 ents_to_delete.merge(quads);
1237 ents_to_delete.merge(subtract(quads_edges, tris_edges));
1238
1239 CHKERR moab.delete_entities(ents_to_delete);
1240
1241 MOFEM_LOG_C("WORLD", Sev::inform, "Write file %s", output_mesh_name);
1242 CHKERR moab.write_file(output_mesh_name, "MOAB");
1243
1244 auto get_tag_handle = [&](auto name, auto size) {
1245 Tag th;
1246 std::vector<double> def_vals(size, 0.0);
1247 CHKERR moab.tag_get_handle(name, size, MB_TYPE_DOUBLE, th,
1248 MB_TAG_CREAT | MB_TAG_SPARSE,
1249 def_vals.data());
1250 return th;
1251 };
1252
1253 if (test_num) {
1254 Range tets;
1255 CHKERR moab.get_entities_by_dimension(0, 3, tets);
1256 EntityHandle tet = tets.front();
1257 std::array<double, 9> internal_stress, actual_stress;
1258 std::array<double, 9> internal_stress_ref, actual_stress_ref;
1259 std::array<double, 2> nb_gauss_pts_ref, contact_area_ref;
1260 switch (test_num) {
1261 case 3:
1262 actual_stress_ref = {0., 0., -100., 0., 0., 0., 0., 0., 0.};
1263 if (strcmp(stress_tag_name, "INTERNAL_STRESS") == 0)
1264 internal_stress_ref = {0., 0., -200., 0., 0., 0., 0., 0., 0.};
1265 else
1266 internal_stress_ref = {0., 0., -100., 0., 0., 0., 0., 0., 0.};
1267 break;
1268 case 4:
1269 nb_gauss_pts_ref = {216, 492};
1270 contact_area_ref = {0.125, 0.25};
1271 break;
1272 default:
1273 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "Test number %d not found",
1274 test_num);
1275 }
1276
1277 auto th_internal_stress = get_tag_handle("MED_INTERNAL_STRESS", 9);
1278 auto th_actual_stress = get_tag_handle("MED_ACTUAL_STRESS", 9);
1279 CHKERR moab.tag_get_data(th_internal_stress, &tet, 1,
1280 internal_stress.data());
1281 CHKERR moab.tag_get_data(th_actual_stress, &tet, 1,
1282 actual_stress.data());
1283 if (test_num == 4) {
1284 const double eps = 1e-3;
1285 if (std::abs(nb_gauss_pts_ref[0] - nb_gauss_pts[0]) > eps) {
1286 SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1287 "Wrong number of active contact gauss pts: should be %d "
1288 "but is %d",
1289 (int)nb_gauss_pts_ref[0], (int)nb_gauss_pts[0]);
1290 }
1291 if (std::abs(nb_gauss_pts_ref[1] - nb_gauss_pts[1]) > eps) {
1292 SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1293 "Wrong total number of contact gauss pts: should be %d "
1294 "but is %d",
1295 (int)nb_gauss_pts_ref[1], (int)nb_gauss_pts[1]);
1296 }
1297 if (std::abs(contact_area_ref[0] - contact_area[0]) > eps) {
1298 SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1299 "Wrong active contact area: should be %g "
1300 "but is %g",
1301 contact_area_ref[0], contact_area[0]);
1302 }
1303 if (std::abs(contact_area_ref[1] - contact_area[1]) > eps) {
1304 SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1305 "Wrong potential contact area: should be %g "
1306 "but is %g",
1307 contact_area_ref[1], contact_area[1]);
1308 }
1309 } else {
1310 if (save_mean_stress) {
1311 const double eps = 1e-10;
1312 for (int i = 0; i < 9; i++) {
1313 if (std::abs(internal_stress[i] - internal_stress_ref[i]) > eps) {
1314 SETERRQ3(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1315 "Wrong component %d of internal stress: should be %g "
1316 "but is %g",
1317 i, internal_stress_ref[i], internal_stress[i]);
1318 }
1319 if (std::abs(actual_stress[i] - actual_stress_ref[i]) > eps) {
1320 SETERRQ3(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1321 "Wrong component %d of actual stress: should be %g "
1322 "but is %g",
1323 i, actual_stress_ref[i], actual_stress[i]);
1324 }
1325 }
1326 }
1327 }
1328 }
1329 }
1330 }
1332
1333 // finish work cleaning memory, getting statistics, etc
1335
1336 return 0;
1337}
#define DM_NO_ELEMENT
Definition: DMMoFEM.hpp:10
const std::string default_options
std::string param_file
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
static const double eps
Kronecker Delta class symmetric.
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ MF_ZERO
Definition: definitions.h:98
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ H1
continuous field
Definition: definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ BLOCKSET
Definition: definitions.h:148
@ MOFEM_NOT_FOUND
Definition: definitions.h:33
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
constexpr int order
@ F
constexpr auto t_kd
void temp(int x, int y=10)
Definition: simple.cpp:4
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1109
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:483
PetscErrorCode DMMoFEMSNESSetFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES residual evaluation function
Definition: DMMoFEM.cpp:704
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.
Definition: DMMoFEM.cpp:118
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:509
PetscErrorCode DMMoFEMSNESSetJacobian(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES Jacobian evaluation function
Definition: DMMoFEM.cpp:745
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:572
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
Definition: DMMoFEM.cpp:1080
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMoFEM.cpp:521
PetscErrorCode DMoFEMLoopFiniteElementsUpAndLowRank(DM dm, const char fe_name[], MoFEM::FEMethod *method, int low_rank, int up_rank, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:553
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:532
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.
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
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_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.
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
auto bit
set bit
FTensor::Index< 'i', SPACE_DIM > i
char out_file_name[255]
double D
FTensor::Index< 'j', 3 > j
char mesh_file_name[255]
static char help[]
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
DEPRECATED auto smartCreateDMMatrix(DM dm)
Definition: DMMoFEM.hpp:996
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.
Definition: SnesCtx.cpp:136
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.
Definition: SnesCtx.cpp:27
DEPRECATED auto smartCreateDMVector(DM dm)
Definition: DMMoFEM.hpp:1011
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
constexpr auto field_name
multi_index_container< boost::shared_ptr< ContactCommonData >, indexed_by< ordered_unique< tag< Prism_tag >, member< ContactCommonData, EntityHandle, &ContactCommonData::pRism > > > > ContactCommonData_multiIndex
Set Dirichlet boundary conditions on spatial displacements.
static MoFEMErrorCode addNeumannBCElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS", Range *intersect_ptr=NULL)
Declare finite element.
static MoFEMErrorCode setMomentumFluxOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, NeumannForcesSurface > &neumann_forces, Vec F, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Set operators to finite elements calculating right hand side vector.
static MoFEMErrorCode setSpringOperators(MoFEM::Interface &m_field, boost::shared_ptr< FaceElementForcesAndSourcesCore > fe_spring_lhs_ptr, boost::shared_ptr< FaceElementForcesAndSourcesCore > fe_spring_rhs_ptr, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS", double stiffness_scale=1.)
Implementation of spring element. Set operators to calculate LHS and RHS.
static MoFEMErrorCode addSpringElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Declare spring element.
Managing BitRefLevels.
virtual int get_comm_size() const =0
virtual moab::Interface & get_moab()=0
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.
virtual int get_comm_rank() const =0
Core (interface) class.
Definition: Core.hpp:82
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
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
Basic algebra on fields.
Definition: FieldBlas.hpp:21
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Create interface from given surface and insert flat prisms in-between.
Projection of edge entities with one mid-node on hierarchical basis.
intrusive_ptr for managing petsc objects
Interface for nonlinear (SNES) solver.
Definition: SnesCtx.hpp:13
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
static MoFEMErrorCode getFileVersion(moab::Interface &moab, Version &version)
Get database major version.
std::string strVersion()
Element used to integrate on slave surfaces. It convects integration points on slaves,...
Element used to integrate on master surfaces. It convects integration points on slaves,...
Postprocess on face.
Element used to integrate on slave surfaces. It convects integration points on slaves,...
Element used to integrate on master surfaces. It convects integration points on slaves,...
Set integration rule.

Variable Documentation

◆ help

char help[] = "\n"
static

Definition at line 26 of file mortar_contact_thermal.cpp.