v0.14.0
Loading...
Searching...
No Matches
Classes | Functions | Variables
simple_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 to volume elements. More...
 

Functions

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

Variables

static char help [] = "\n"
 

Function Documentation

◆ main()

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

Post proc on the skin

Definition at line 39 of file simple_contact_thermal.cpp.

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

Variable Documentation

◆ help

char help[] = "\n"
static

Definition at line 30 of file simple_contact_thermal.cpp.

Generated by Doxygen 1.12.0 and hosted at University of Glasgow