v0.14.0
Loading...
Searching...
No Matches
simple_contact_thermal.cpp
Go to the documentation of this file.
1/** \file simple_contact_thermal.cpp
2 * \example simple_contact_thermal.cpp
3 *
4 * Implementation of simple contact between surfaces with matching meshes
5 * taking into account internal stress resulting from the thermal expansion
6 *
7 **/
8
9/* MoFEM is free software: you can redistribute it and/or modify it under
10 * the terms of the GNU Lesser General Public License as published by the
11 * Free Software Foundation, either version 3 of the License, or (at your
12 * option) any later version.
13 *
14 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
15 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
16 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
17 * License for more details.
18 *
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>.
21 */
22
23#include <Mortar.hpp>
25
26using namespace MoFEM;
27using namespace std;
28
29
30static char help[] = "\n";
34
35struct VolRule {
36 int operator()(int, int, int order) const { return 2 * (order - 1); }
37};
38
39int main(int argc, char *argv[]) {
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;
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) {
251 auto prism_interface = m_field.getInterface<PrismInterface>();
252
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
300 };
301
302 auto find_contact_prisms = [&](std::vector<BitRefLevel> &bit_levels,
303 Range &contact_prisms, Range &master_tris,
304 Range &slave_tris) {
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
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) {
338 Range all_tets, all_nodes;
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 }
362 };
363
364 auto make_wavy_surface = [&](int block_id, int dim, double lambda,
365 double delta, double height) {
367 Range all_tets, all_nodes;
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 }
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) {
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
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(
520 boost::shared_ptr<ForcesAndSourcesCore> fe_elastic_rhs_ptr(
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](
539 FTensor::Index<'i', 3> i;
540 FTensor::Index<'j', 3> j;
541
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(
575 "SPATIAL_POSITION", "SPATIAL_POSITION", data_hooke_element_at_pts,
576 thermal_strain));
577 } else {
578 fe_elastic_rhs_ptr->getOpPtrVector().push_back(
580 "SPATIAL_POSITION", "SPATIAL_POSITION", data_hooke_element_at_pts,
581 moab, stress_tag_name));
582 fe_elastic_rhs_ptr->getOpPtrVector().push_back(
584 "SPATIAL_POSITION", data_hooke_element_at_pts));
585 }
586
587 fe_elastic_rhs_ptr->getOpPtrVector().push_back(
589 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat));
590 fe_elastic_rhs_ptr->getOpPtrVector().push_back(
592 "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
593 fe_elastic_rhs_ptr->getOpPtrVector().push_back(
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;
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");
614 "SPATIAL_POSITION");
616 "SPATIAL_POSITION");
618 "SPATIAL_POSITION");
620 "MESH_NODE_POSITIONS");
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<
632 m_field, "SPATIAL_POSITION", "MESH_NODE_POSITIONS");
633 };
634
635 auto make_convective_slave_element = [&]() {
636 return boost::make_shared<
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
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
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;
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(
796 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr(
798
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);
808 dirichlet_bc_ptr.get(), NULL);
809 if (convect_pts == PETSC_TRUE) {
811 dm, "CONTACT_ELEM",
812 get_contact_rhs(contact_problem, make_convective_master_element),
813 PETSC_NULL, PETSC_NULL);
815 dm, "CONTACT_ELEM",
816 get_master_traction_rhs(contact_problem,
817 make_convective_slave_element),
818 PETSC_NULL, PETSC_NULL);
819 } else {
821 dm, "CONTACT_ELEM",
822 get_contact_rhs(contact_problem, make_contact_element, alm_flag),
823 PETSC_NULL, PETSC_NULL);
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) {
844 dm, "CONTACT_ELEM",
845 get_contact_lhs(contact_problem, make_convective_master_element),
846 PETSC_NULL, PETSC_NULL);
848 dm, "CONTACT_ELEM",
849 get_master_traction_lhs(contact_problem,
850 make_convective_slave_element),
851 PETSC_NULL, PETSC_NULL);
852 } else {
854 dm, "CONTACT_ELEM",
855 get_contact_lhs(contact_problem, make_contact_element, alm_flag),
856 PETSC_NULL, PETSC_NULL);
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,
918 if (type != MBVERTEX)
920 CHKERR loopSideVolumes(feVolName, *sideOpFe);
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(
929 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat));
930 my_vol_side_fe_ptr->getOpPtrVector().push_back(
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(
939 "SPATIAL_POSITION", data_hooke_element_at_pts->spatPosMat));
940 post_proc_skin.getOpPtrVector().push_back(
942 "MESH_NODE_POSITIONS", data_hooke_element_at_pts->meshNodePosMat));
943
944 post_proc_skin.getOpPtrVector().push_back(
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 {
956 CHKERR PetscPrintf(PETSC_COMM_WORLD, "Ignoring pressure...\n");
957 }
958 CHKERR PetscPrintf(PETSC_COMM_WORLD, "Load scale: %6.4e\n",
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) {
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);
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) {
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 }
1211
1212 // finish work cleaning memory, getting statistics, etc
1214
1215 return 0;
1216}
#define DM_NO_ELEMENT
Definition DMMoFEM.hpp:10
const std::string default_options
std::string param_file
int main()
static const double eps
Kronecker Delta class symmetric.
#define CATCH_ERRORS
Catch errors.
@ 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()
@ 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 ...
@ BLOCKSET
@ MOFEM_NOT_FOUND
Definition definitions.h:33
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
constexpr int order
const int dim
@ 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
MoFEMErrorCode setBitRefLevelByDim(const EntityHandle meshset, const int dim, const BitRefLevel bit, int verb=QUIET) const
Set the Bit Ref Level By Dim object.
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.
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
MoFEMErrorCode writeFile(const std::string file_name, const char *file_type="MOAB", const char *file_options="PARALLEL=WRITE_PART")
wrote results in (MOAB) format, use "file_name.h5m"
#define MOFEM_LOG(channel, severity)
Log.
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
static double lambda
char out_file_name[255]
double D
FTensor::Index< 'j', 3 > j
char mesh_file_name[255]
VolumeElementForcesAndSourcesCoreOnSide VolSideFe
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition Types.hpp:40
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
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:139
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
static constexpr double delta
static char help[]
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
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
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()
data for calculation heat conductivity and heat capacity elements
Postprocess on face.
MoFEMErrorCode generateReferenceElementMesh()
std::vector< EntityHandle > mapGaussPts
moab::Interface & postProcMesh
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.
int operator()(int, int, int order) const