35 {
36
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
64 if (!
static_cast<bool>(ifstream(
param_file))) {
65 std::ofstream file(
param_file.c_str(), std::ios::ate);
66 if (file.is_open()) {
68 file.close();
69 }
70 }
71
72
74
75
76 moab::Core mb_instance;
77 moab::Interface &moab = mb_instance;
78
79 try {
80 PetscBool flg_file;
81 PetscBool flg_file_out;
82
84 char output_mesh_name[255];
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",
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);
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);
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();
182
183
184 if (flg_file != PETSC_TRUE) {
185 SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -file_name (MESH FILE NEEDED)");
186 }
187
188 if (is_partitioned) {
190 "Partitioned mesh is not supported");
191 }
192
193 const char *option;
194 option = "";
199
200
203
204 std::vector<BitRefLevel> bit_levels;
207 CHKERR bit_ref_manager->setBitRefLevelByDim(0, 3, bit_levels.back());
208
209 auto add_prism_interface = [&](std::vector<BitRefLevel> &bit_levels) {
212
214 if (cit->getName().compare(0, 11, "INT_CONTACT") == 0) {
215 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Insert %s (id: %d)\n",
216 cit->getName().c_str(), cit->getMeshsetId());
218
219
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(
226 ref_level_meshset);
227
228
229 CHKERR prism_interface->getSides(cubit_meshset, bit_levels.back(),
230 true, 0);
231
232 bit_levels.push_back(
BitRefLevel().set(bit_levels.size()));
233
234 CHKERR prism_interface->splitSides(ref_level_meshset,
235 bit_levels.back(), cubit_meshset,
236 true, true, 0);
237
238 CHKERR moab.delete_entities(&ref_level_meshset, 1);
239
240
243 CHKERR bit_ref_manager->updateMeshsetByEntitiesChildren(
244 cubit_meshset, bit_levels.back(), cubit_meshset, MBVERTEX,
245 true);
246 CHKERR bit_ref_manager->updateMeshsetByEntitiesChildren(
247 cubit_meshset, bit_levels.back(), cubit_meshset, MBEDGE, true);
248 CHKERR bit_ref_manager->updateMeshsetByEntitiesChildren(
249 cubit_meshset, bit_levels.back(), cubit_meshset, MBTRI, true);
250 CHKERR bit_ref_manager->updateMeshsetByEntitiesChildren(
251 cubit_meshset, bit_levels.back(), cubit_meshset, MBTET, true);
252 }
253
254 CHKERR bit_ref_manager->shiftRightBitRef(1);
255 bit_levels.pop_back();
256 }
257 }
258
260 };
261
262 auto find_contact_prisms = [&](std::vector<BitRefLevel> &bit_levels,
263 Range &simple_contact_prisms,
264 Range &simple_master_tris,
265 Range &simple_slave_tris) {
267
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
276 for (Range::iterator pit = simple_contact_prisms.begin();
277 pit != simple_contact_prisms.end(); pit++) {
278 CHKERR moab.side_element(*pit, 2, 3, tri);
279 simple_master_tris.insert(tri);
280 CHKERR moab.side_element(*pit, 2, 4, tri);
281 simple_slave_tris.insert(tri);
282 }
283
285 };
286
287 Range simple_contact_prisms, simple_master_tris, simple_slave_tris;
288
289 if (!ignore_contact) {
290 CHKERR add_prism_interface(bit_levels);
291 CHKERR find_contact_prisms(bit_levels, simple_contact_prisms,
292 simple_master_tris, simple_slave_tris);
293 }
294
295 Range mortar_contact_prisms, mortar_master_tris, mortar_slave_tris;
296
298 if (it->getName().compare(0, 13, "MORTAR_MASTER") == 0) {
300 it->meshset, MBTRI, mortar_master_tris, true);
301 }
302 }
303
305 if (it->getName().compare(0, 12, "MORTAR_SLAVE") == 0) {
306 CHKERR m_field.
get_moab().get_entities_by_type(it->meshset, MBTRI,
307 mortar_slave_tris, true);
308 }
309 }
310
312
313 boost::shared_ptr<ContactSearchKdTree::ContactCommonData_multiIndex>
314 contact_commondata_multi_index;
315
316 contact_commondata_multi_index =
317 boost::shared_ptr<ContactSearchKdTree::ContactCommonData_multiIndex>(
319
320 CHKERR contact_search_kd_tree.buildTree(mortar_master_tris);
321
322 CHKERR contact_search_kd_tree.contactSearchAlgorithm(
323 mortar_master_tris, mortar_slave_tris, contact_commondata_multi_index,
324 mortar_contact_prisms);
325
327 CHKERR moab.create_meshset(MESHSET_SET, meshset_slave_master_prisms);
328
330 moab.add_entities(meshset_slave_master_prisms, mortar_contact_prisms);
331
333 meshset_slave_master_prisms, 3, bit_levels.back());
334
335 Range tris_from_prism;
336
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
350
353
354 if (!eigen_pos_flag) {
357 }
358
361
367
368
374
375 if (!eigen_pos_flag) {
381 }
382
384
388
389 auto set_contact_order = [&](
Range &contact_prisms,
int order_contact,
390 int nb_ho_levels) {
392 Range contact_tris, contact_edges;
393 CHKERR moab.get_adjacencies(contact_prisms, 2,
false, contact_tris,
394 moab::Interface::UNION);
395 contact_tris = contact_tris.subset_by_type(MBTRI);
396 CHKERR moab.get_adjacencies(contact_tris, 1,
false, contact_edges,
397 moab::Interface::UNION);
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
416 order_contact);
417
419 };
420
421 if (!ignore_contact && order_contact >
order) {
422 CHKERR set_contact_order(all_contact_prisms, order_contact, nb_ho_levels);
423 }
424
425
427
428
429 {
432 }
433
434 {
437 }
438
440 CHKERR moab.get_connectivity(all_slave_tris, slave_verts,
true);
442 slave_verts, "LAGMULT");
443
444
445 boost::shared_ptr<std::map<int, BlockData>> block_sets_ptr =
446 boost::make_shared<std::map<int, BlockData>>();
447 CHKERR HookeElement::setBlocks(m_field, block_sets_ptr);
448
449 boost::shared_ptr<ForcesAndSourcesCore> fe_elastic_lhs_ptr(
451 boost::shared_ptr<ForcesAndSourcesCore> fe_elastic_rhs_ptr(
453 fe_elastic_lhs_ptr->getRuleHook =
VolRule();
454 fe_elastic_rhs_ptr->getRuleHook =
VolRule();
455
456 CHKERR HookeElement::addElasticElement(m_field, block_sets_ptr,
"ELASTIC",
457 "SPATIAL_POSITION",
458 "MESH_NODE_POSITIONS", false);
459
460 auto data_hooke_element_at_pts =
461 boost::make_shared<HookeInternalStressElement::DataAtIntegrationPts>();
462 CHKERR HookeElement::setOperators(fe_elastic_lhs_ptr, fe_elastic_rhs_ptr,
463 block_sets_ptr, "SPATIAL_POSITION",
464 "MESH_NODE_POSITIONS", false, false,
465 MBTET, data_hooke_element_at_pts);
466 auto thermal_strain =
467 [&thermal_expansion_coef, &init_temp, &test_num](
472
475
476 t_thermal_strain(
i,
j) = 0.0;
477
478 switch (test_num) {
479 case 0:
480
481 temp = init_temp - 1.0;
482 t_thermal_strain(
i,
j) =
483 thermal_expansion_coef * (
temp - init_temp) *
t_kd(
i,
j);
484 break;
485 case 1:
486 case 2:
487 t_thermal_strain(
i,
j) = -thermal_expansion_coef *
t_kd(
i,
j);
488 break;
489 case 3:
490 t_thermal_strain(2, 2) = thermal_expansion_coef;
491 break;
492 case 4:
493 t_thermal_strain(
i,
j) = thermal_expansion_coef *
t_kd(
i,
j);
494 break;
495 default:
496 break;
497 }
498
499 return t_thermal_strain;
500 };
501
502 if (analytical_input) {
503 fe_elastic_rhs_ptr->getOpPtrVector().push_back(
504 new HookeElement::OpAnalyticalInternalStrain_dx<0>(
505 "SPATIAL_POSITION", data_hooke_element_at_pts, thermal_strain));
506 fe_elastic_rhs_ptr->getOpPtrVector().push_back(
508 "SPATIAL_POSITION", "SPATIAL_POSITION", data_hooke_element_at_pts,
509 thermal_strain));
510 } else {
511 fe_elastic_rhs_ptr->getOpPtrVector().push_back(
513 "SPATIAL_POSITION", "SPATIAL_POSITION", data_hooke_element_at_pts,
514 moab, stress_tag_name));
515 fe_elastic_rhs_ptr->getOpPtrVector().push_back(
517 "SPATIAL_POSITION", data_hooke_element_at_pts));
518 }
519
520 fe_elastic_rhs_ptr->getOpPtrVector().push_back(
522 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat));
523 fe_elastic_rhs_ptr->getOpPtrVector().push_back(
525 "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
526 fe_elastic_rhs_ptr->getOpPtrVector().push_back(
528 "SPATIAL_POSITION", "SPATIAL_POSITION", data_hooke_element_at_pts,
529 *block_sets_ptr.get(), moab, scale_factor, save_mean_stress, false,
530 false));
531
534 if (
bit->getName().compare(0, 11,
"MAT_ELASTIC") == 0) {
536 CHKERR moab.get_entities_by_dimension(
bit->getMeshset(), 3, tets,
true);
537 all_tets.merge(tets);
538 }
539 }
540 Skinner skinner(&moab);
542 CHKERR skinner.find_skin(0, all_tets,
false, skin_tris);
543
547 "SPATIAL_POSITION");
549 "SPATIAL_POSITION");
551 "SPATIAL_POSITION");
553 "MESH_NODE_POSITIONS");
555 "EIGEN_POSITIONS");
556
557 auto make_mortar_contact_element = [&]() {
558 return boost::make_shared<MortarContactProblem::MortarContactElement>(
559 m_field, contact_commondata_multi_index, "SPATIAL_POSITION",
560 "MESH_NODE_POSITIONS");
561 };
562
563 auto make_convective_mortar_master_element = [&]() {
564 return boost::make_shared<
566 m_field, contact_commondata_multi_index, "SPATIAL_POSITION",
567 "MESH_NODE_POSITIONS");
568 };
569
570 auto make_convective_mortar_slave_element = [&]() {
571 return boost::make_shared<
573 m_field, contact_commondata_multi_index, "SPATIAL_POSITION",
574 "MESH_NODE_POSITIONS");
575 };
576
577 auto make_mortar_contact_common_data = [&]() {
578 return boost::make_shared<MortarContactProblem::CommonDataMortarContact>(
579 m_field);
580 };
581
582 auto get_mortar_contact_rhs = [&](auto contact_problem, auto make_element,
583 bool is_alm = false) {
584 auto fe_rhs_mortar_contact = make_element();
585 auto common_data_mortar_contact = make_mortar_contact_common_data();
586 if (print_contact_state) {
587 fe_rhs_mortar_contact->contactStateVec =
588 common_data_mortar_contact->gaussPtsStateVec;
589 }
590 contact_problem->setContactOperatorsRhs(
591 fe_rhs_mortar_contact, common_data_mortar_contact, "SPATIAL_POSITION",
592 "LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS");
593 return fe_rhs_mortar_contact;
594 };
595
596 auto get_mortar_master_traction_rhs = [&](auto contact_problem,
597 auto make_element,
598 bool is_alm = false) {
599 auto fe_rhs_mortar_contact = make_element();
600 auto common_data_mortar_contact = make_mortar_contact_common_data();
601 contact_problem->setMasterForceOperatorsRhs(
602 fe_rhs_mortar_contact, common_data_mortar_contact, "SPATIAL_POSITION",
603 "LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS");
604 return fe_rhs_mortar_contact;
605 };
606
607 auto get_mortar_master_traction_lhs = [&](auto contact_problem,
608 auto make_element,
609 bool is_alm = false) {
610 auto fe_lhs_mortar_contact = make_element();
611 auto common_data_mortar_contact = make_mortar_contact_common_data();
612 contact_problem->setMasterForceOperatorsLhs(
613 fe_lhs_mortar_contact, common_data_mortar_contact, "SPATIAL_POSITION",
614 "LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS");
615 return fe_lhs_mortar_contact;
616 };
617
618 auto get_mortar_master_help_traction_lhs =
619 [&](auto contact_problem, auto make_element, bool is_alm = false) {
620 auto fe_lhs_mortar_contact = make_element();
621 auto common_data_mortar_contact = make_mortar_contact_common_data();
622 contact_problem->setMasterForceOperatorsLhs(
623 fe_lhs_mortar_contact, common_data_mortar_contact,
624 "SPATIAL_POSITION", "LAGMULT", is_alm);
625 return fe_lhs_mortar_contact;
626 };
627
628 auto get_mortar_contact_lhs = [&](auto contact_problem, auto make_element,
629 bool is_alm = false) {
630 auto fe_lhs_mortar_contact = make_element();
631 auto common_data_mortar_contact = make_mortar_contact_common_data();
632 contact_problem->setContactOperatorsLhs(
633 fe_lhs_mortar_contact, common_data_mortar_contact, "SPATIAL_POSITION",
634 "LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS");
635 return fe_lhs_mortar_contact;
636 };
637
638 auto get_mortar_contact_help_lhs =
639 [&](auto contact_problem, auto make_element, bool is_alm = false) {
640 auto fe_lhs_mortar_contact = make_element();
641 auto common_data_mortar_contact = make_mortar_contact_common_data();
642 contact_problem->setContactOperatorsLhs(
643 fe_lhs_mortar_contact, common_data_mortar_contact,
644 "SPATIAL_POSITION", "LAGMULT", is_alm);
645 return fe_lhs_mortar_contact;
646 };
647
648 auto cn_value_ptr = boost::make_shared<double>(cn_value);
649 auto mortar_contact_problem = boost::make_shared<MortarContactProblem>(
650 m_field, contact_commondata_multi_index, cn_value_ptr, is_newton_cotes);
651
652
653
654 if (!eigen_pos_flag)
655 mortar_contact_problem->addMortarContactElement(
656 "MORTAR_CONTACT_ELEM", "SPATIAL_POSITION", "LAGMULT",
657 mortar_contact_prisms);
658 else
659 mortar_contact_problem->addMortarContactElement(
660 "MORTAR_CONTACT_ELEM", "SPATIAL_POSITION", "LAGMULT",
661 mortar_contact_prisms, "MESH_NODE_POSITIONS", eigen_pos_flag,
662 "EIGEN_POSITIONS");
663
664 auto make_simple_contact_element = [&]() {
665 return boost::make_shared<SimpleContactProblem::SimpleContactElement>(
666 m_field);
667 };
668
669 auto make_convective_simple_master_element = [&]() {
670 return boost::make_shared<
672 m_field, "SPATIAL_POSITION", "MESH_NODE_POSITIONS");
673 };
674
675 auto make_convective_simple_slave_element = [&]() {
676 return boost::make_shared<
678 m_field, "SPATIAL_POSITION", "MESH_NODE_POSITIONS");
679 };
680
681 auto make_simple_contact_common_data = [&]() {
682 return boost::make_shared<SimpleContactProblem::CommonDataSimpleContact>(
683 m_field);
684 };
685
686 auto get_simple_contact_rhs = [&](auto contact_problem, auto make_element,
687 bool is_alm = false) {
688 auto fe_rhs_simple_contact = make_element();
689 auto common_data_simple_contact = make_simple_contact_common_data();
690 if (print_contact_state) {
691 fe_rhs_simple_contact->contactStateVec =
692 common_data_simple_contact->gaussPtsStateVec;
693 }
694 contact_problem->setContactOperatorsRhs(
695 fe_rhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
696 "LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS");
697 return fe_rhs_simple_contact;
698 };
699
700 auto get_simple_master_traction_rhs = [&](auto contact_problem,
701 auto make_element,
702 bool is_alm = false) {
703 auto fe_rhs_simple_contact = make_element();
704 auto common_data_simple_contact = make_simple_contact_common_data();
705 contact_problem->setMasterForceOperatorsRhs(
706 fe_rhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
707 "LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS");
708 return fe_rhs_simple_contact;
709 };
710
711 auto get_simple_master_traction_lhs =
712 [&](auto contact_problem, auto make_element, bool is_alm = false) {
713 auto fe_lhs_simple_contact = make_element();
714 auto common_data_simple_contact = make_simple_contact_common_data();
715 contact_problem->setMasterForceOperatorsLhs(
716 fe_lhs_simple_contact, common_data_simple_contact,
717 "SPATIAL_POSITION", "LAGMULT", is_alm);
718 return fe_lhs_simple_contact;
719 };
720
721 auto get_contact_lhs = [&](auto contact_problem, auto make_element,
722 bool is_alm = false) {
723 auto fe_lhs_simple_contact = make_element();
724 auto common_data_simple_contact = make_simple_contact_common_data();
725 contact_problem->setContactOperatorsLhs(
726 fe_lhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
727 "LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS");
728 return fe_lhs_simple_contact;
729 };
730
731 auto simple_contact_problem = boost::make_shared<SimpleContactProblem>(
732 m_field, cn_value_ptr, is_newton_cotes);
733 simple_contact_problem->addContactElement("SIMPLE_CONTACT_ELEM",
734 "SPATIAL_POSITION", "LAGMULT",
735 simple_contact_prisms);
736
737 simple_contact_problem->addPostProcContactElement(
738 "CONTACT_POST_PROC", "SPATIAL_POSITION", "LAGMULT",
739 "MESH_NODE_POSITIONS", all_slave_tris);
740
742
743
745 "MESH_NODE_POSITIONS");
746
747
749
750
752
753
755
756
758 bit_levels.back());
759
760 DMType dm_name = "DMMOFEM";
762
765
766
767 CHKERR DMSetType(dm, dm_name);
768
769
771 CHKERR DMSetFromOptions(dm);
773
781
783
784
787
788
790
793 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
794 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
795
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
807 boost::shared_ptr<FEMethod> dirichlet_bc_ptr =
809 m_field,
"SPATIAL_POSITION", Aij,
D,
F));
810
812 dirichlet_bc_ptr->snes_x =
D;
813
814
815 boost::ptr_map<std::string, NeumannForcesSurface> neumann_forces;
817 m_field, neumann_forces, NULL, "SPATIAL_POSITION");
818
819 boost::ptr_map<std::string, NeumannForcesSurface>::iterator mit =
820 neumann_forces.begin();
821 for (; mit != neumann_forces.end(); mit++) {
824 &mit->second->getLoopFe(), NULL, NULL);
825 }
826
827
828
829 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr(
831 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr(
833
835 m_field, fe_spring_lhs_ptr, fe_spring_rhs_ptr, "SPATIAL_POSITION",
836 "MESH_NODE_POSITIONS");
837
839 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
840 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
843 dirichlet_bc_ptr.get(), NULL);
844 if (convect_pts == PETSC_TRUE) {
846 dm, "MORTAR_CONTACT_ELEM",
847 get_mortar_contact_rhs(mortar_contact_problem,
848 make_convective_mortar_master_element),
849 PETSC_NULL, PETSC_NULL);
851 dm, "MORTAR_CONTACT_ELEM",
852 get_mortar_master_traction_rhs(mortar_contact_problem,
853 make_convective_mortar_slave_element),
854 PETSC_NULL, PETSC_NULL);
855 } else {
857 dm, "MORTAR_CONTACT_ELEM",
858 get_mortar_contact_rhs(mortar_contact_problem,
859 make_mortar_contact_element, alm_flag),
860 PETSC_NULL, PETSC_NULL);
862 dm, "MORTAR_CONTACT_ELEM",
863 get_mortar_master_traction_rhs(mortar_contact_problem,
864 make_mortar_contact_element, alm_flag),
865 PETSC_NULL, PETSC_NULL);
866 }
867
868 if (convect_pts == PETSC_TRUE) {
870 dm, "SIMPLE_CONTACT_ELEM",
871 get_simple_contact_rhs(simple_contact_problem,
872 make_convective_simple_master_element),
873 PETSC_NULL, PETSC_NULL);
875 dm, "SIMPLE_CONTACT_ELEM",
876 get_simple_master_traction_rhs(simple_contact_problem,
877 make_convective_simple_slave_element),
878 PETSC_NULL, PETSC_NULL);
879 } else {
881 dm, "SIMPLE_CONTACT_ELEM",
882 get_simple_contact_rhs(simple_contact_problem,
883 make_simple_contact_element, alm_flag),
884 PETSC_NULL, PETSC_NULL);
886 dm, "SIMPLE_CONTACT_ELEM",
887 get_simple_master_traction_rhs(simple_contact_problem,
888 make_simple_contact_element, alm_flag),
889 PETSC_NULL, PETSC_NULL);
890 }
891
893 PETSC_NULL);
894
896 PETSC_NULL);
898 dirichlet_bc_ptr.get());
899
900 boost::shared_ptr<FEMethod> fe_null;
902 fe_null);
903 if (convect_pts == PETSC_TRUE) {
905 dm, "MORTAR_CONTACT_ELEM",
906 get_mortar_contact_help_lhs(mortar_contact_problem,
907 make_convective_mortar_master_element),
908 PETSC_NULL, PETSC_NULL);
910 dm, "MORTAR_CONTACT_ELEM",
911 get_mortar_master_help_traction_lhs(
912 mortar_contact_problem, make_convective_mortar_slave_element),
913 PETSC_NULL, PETSC_NULL);
914 } else {
916 dm, "MORTAR_CONTACT_ELEM",
917 get_mortar_contact_lhs(mortar_contact_problem,
918 make_mortar_contact_element, alm_flag),
919 PETSC_NULL, PETSC_NULL);
921 dm, "MORTAR_CONTACT_ELEM",
922 get_mortar_master_traction_lhs(mortar_contact_problem,
923 make_mortar_contact_element, alm_flag),
924 PETSC_NULL, PETSC_NULL);
925 }
926
927 if (convect_pts == PETSC_TRUE) {
929 dm, "SIMPLE_CONTACT_ELEM",
930 get_contact_lhs(simple_contact_problem,
931 make_convective_simple_master_element),
932 PETSC_NULL, PETSC_NULL);
934 dm, "SIMPLE_CONTACT_ELEM",
935 get_simple_master_traction_lhs(simple_contact_problem,
936 make_convective_simple_slave_element),
937 PETSC_NULL, PETSC_NULL);
938 } else {
940 get_contact_lhs(simple_contact_problem,
941 make_simple_contact_element,
942 alm_flag),
943 PETSC_NULL, PETSC_NULL);
945 dm, "SIMPLE_CONTACT_ELEM",
946 get_simple_master_traction_lhs(simple_contact_problem,
947 make_simple_contact_element, alm_flag),
948 PETSC_NULL, PETSC_NULL);
949 }
951 PETSC_NULL);
953 PETSC_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
970 CHKERR SNESSetDM(snes, dm);
972
973 {
974 CHKERR SNESSetDM(snes, dm);
978 CHKERR SNESSetFromOptions(snes);
979 }
980
981
983 CHKERR post_proc_skin.generateReferenceElementMesh();
985 CHKERR post_proc_skin.addFieldValuesPostProc(
"SPATIAL_POSITION");
986 CHKERR post_proc_skin.addFieldValuesPostProc(
"MESH_NODE_POSITIONS");
987 CHKERR post_proc_skin.addFieldValuesPostProc(
"EIGEN_POSITIONS");
988
989 struct OpGetFieldGradientValuesOnSkin
990 : public FaceElementForcesAndSourcesCore::UserDataOperator {
991
992 const std::string feVolName;
993 boost::shared_ptr<VolSideFe> sideOpFe;
994
995 OpGetFieldGradientValuesOnSkin(
const std::string
field_name,
996 const std::string vol_fe_name,
997 boost::shared_ptr<VolSideFe> side_fe)
1000 feVolName(vol_fe_name), sideOpFe(side_fe) {}
1001
1005 if (type != MBVERTEX)
1007 CHKERR loopSideVolumes(feVolName, *sideOpFe);
1009 }
1010 };
1011
1012 boost::shared_ptr<VolSideFe> vol_side_fe_ptr =
1013 boost::make_shared<VolSideFe>(m_field);
1014 vol_side_fe_ptr->getOpPtrVector().push_back(
1016 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat));
1017 vol_side_fe_ptr->getOpPtrVector().push_back(
1019 "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
1020
1021 post_proc_skin.getOpPtrVector().push_back(
1022 new OpGetFieldGradientValuesOnSkin("SPATIAL_POSITION", "ELASTIC",
1023 vol_side_fe_ptr));
1024 post_proc_skin.getOpPtrVector().push_back(
1026 "SPATIAL_POSITION", data_hooke_element_at_pts->spatPosMat));
1027 post_proc_skin.getOpPtrVector().push_back(
1029 "MESH_NODE_POSITIONS", data_hooke_element_at_pts->meshNodePosMat));
1030
1031 post_proc_skin.getOpPtrVector().push_back(
1034 "SPATIAL_POSITION", data_hooke_element_at_pts,
1035 *block_sets_ptr.get(), post_proc_skin.postProcMesh,
1036 post_proc_skin.mapGaussPts, false, false));
1037
1038 for (int ss = 0; ss != nb_sub_steps; ++ss) {
1039 if (!ignore_pressure) {
1041 } else {
1043 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Ignoring pressure...\n");
1044 }
1045 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Load scale: %6.4e\n",
1047
1048 CHKERR SNESSolve(snes, PETSC_NULL,
D);
1049
1050 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
1051 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
1052 }
1053
1054
1056
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
1064 MOFEM_LOG_C(
"WORLD", Sev::inform,
"Elastic energy: %8.8e", energy);
1065
1066
1067 moab::Core mb_post;
1068 moab::Interface &moab_proc = mb_post;
1069
1070 auto common_data_mortar_contact = make_mortar_contact_common_data();
1071
1072 boost::shared_ptr<MortarContactProblem::MortarContactElement>
1073 fe_post_proc_mortar_contact;
1074 if (convect_pts == PETSC_TRUE) {
1075 fe_post_proc_mortar_contact = make_convective_mortar_master_element();
1076 } else {
1077 fe_post_proc_mortar_contact = make_mortar_contact_element();
1078 }
1079
1080 std::array<double, 2> nb_gauss_pts;
1081 std::array<double, 2> contact_area;
1082
1083 auto get_contact_data = [&](auto vec, std::array<double, 2> &data) {
1085 CHKERR VecAssemblyBegin(vec);
1086 CHKERR VecAssemblyEnd(vec);
1087 const double *array;
1088 CHKERR VecGetArrayRead(vec, &array);
1090 for (
int i : {0, 1})
1092 }
1093 CHKERR VecRestoreArrayRead(vec, &array);
1095 };
1096
1097 mb_post.delete_mesh();
1098
1099 if (!mortar_contact_prisms.empty()) {
1100 mortar_contact_problem->setContactOperatorsForPostProc(
1101 fe_post_proc_mortar_contact, common_data_mortar_contact, m_field,
1102 "SPATIAL_POSITION", "LAGMULT", mb_post, alm_flag, eigen_pos_flag,
1103 "EIGEN_POSITIONS");
1104
1105 CHKERR VecZeroEntries(common_data_mortar_contact->gaussPtsStateVec);
1106 CHKERR VecZeroEntries(common_data_mortar_contact->contactAreaVec);
1107
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
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
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
1158 "Active gauss pts (matching): %d out of %d",
1159 (int)nb_gauss_pts[0], (int)nb_gauss_pts[1]);
1161 "SELF", Sev::inform,
1162 "Active contact area (matching): %8.8f out of %8.8f (%8.8f% %)",
1163 contact_area[0], contact_area[1],
1164 contact_area[0] / contact_area[1] * 100.);
1165 }
1166 }
1167
1168 if (!ignore_contact) {
1170 std::ostringstream strm;
1171 strm << "out_contact_integ_pts"
1172 << ".h5m";
1174 MOFEM_LOG_C(
"WORLD", Sev::inform,
"Write file %s\n",
1177 "PARALLEL=WRITE_PART");
1178 }
1179
1180 boost::shared_ptr<PostProcFaceOnRefinedMesh> post_proc_contact_ptr(
1182
1183 CHKERR post_proc_contact_ptr->generateReferenceElementMesh();
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()) {
1192 post_proc_contact_ptr);
1194 std::ostringstream stm;
1195 stm << "out_contact_surface"
1196 << ".h5m";
1198 MOFEM_LOG_C(
"WORLD", Sev::inform,
"Write file %s",
1200 CHKERR post_proc_contact_ptr->postProcMesh.write_file(
1202 }
1203
1205 "EIGEN_POSITIONS");
1206
1207 {
1208 PetscPrintf(PETSC_COMM_WORLD, "Loop post proc on the skin\n");
1210 ostringstream stm;
1212 stm << "out_skin.h5m";
1214 MOFEM_LOG_C(
"WORLD", Sev::inform,
"Write file %s",
1216 CHKERR post_proc_skin.writeFile(stm.str());
1217 }
1218
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) {
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());
1251 };
1252
1253 if (test_num) {
1255 CHKERR moab.get_entities_by_dimension(0, 3, tets);
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) {
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) {
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) {
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) {
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) {
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) {
1321 "Wrong component %d of actual stress: should be %g "
1322 "but is %g",
1323 i, actual_stress_ref[
i], actual_stress[
i]);
1324 }
1325 }
1326 }
1327 }
1328 }
1329 }
1330 }
1332
1333
1335
1336 return 0;
1337}
const std::string default_options
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MOFEM_LOG_C(channel, severity, format,...)
Kronecker Delta class symmetric.
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_ATOM_TEST_INVALID
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
void temp(int x, int y=10)
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
PetscErrorCode DMMoFEMSNESSetFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES residual evaluation function
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.
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMMoFEMSNESSetJacobian(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES Jacobian evaluation function
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
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.
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
#define MOFEM_LOG(channel, severity)
Log.
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
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
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.
DEPRECATED auto smartCreateDMMatrix(DM dm)
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.
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.
DEPRECATED auto smartCreateDMVector(DM dm)
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
constexpr auto field_name
Set Dirichlet boundary conditions on spatial displacements.
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
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
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.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
static MoFEMErrorCode getFileVersion(moab::Interface &moab, Version &version)
Get database major version.
Volume finite element base.