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

Go to the source code of this file.

Classes

struct  VolRule
 Set integration rule. More...
 

Typedefs

using BlockData = NonlinearElasticElement::BlockData
 
using VolSideFe = VolumeElementForcesAndSourcesCoreOnSide
 

Functions

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

Variables

static char help [] = "\n"
 

Typedef Documentation

◆ BlockData

using BlockData = NonlinearElasticElement::BlockData

Definition at line 28 of file mortar_contact_thermal.cpp.

◆ VolSideFe

using VolSideFe = VolumeElementForcesAndSourcesCoreOnSide

Definition at line 29 of file mortar_contact_thermal.cpp.

Function Documentation

◆ main()

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

Post proc on the skin

Definition at line 35 of file mortar_contact_thermal.cpp.

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

Variable Documentation

◆ help

char help[] = "\n"
static

Definition at line 26 of file mortar_contact_thermal.cpp.

Generated by Doxygen 1.9.5 and hosted at University of Glasgow