33 {
34
35 const string default_options = "-ksp_type fgmres \n"
36 "-pc_type lu \n"
37 "-pc_factor_mat_solver_type mumps \n"
38 "-mat_mumps_icntl_13 1 \n"
39 "-mat_mumps_icntl_14 800 \n"
40 "-mat_mumps_icntl_20 0 \n"
41 "-mat_mumps_icntl_24 1 \n"
42 "-snes_type newtonls \n"
43 "-snes_linesearch_type basic \n"
44 "-snes_divergence_tolerance 0 \n"
45 "-snes_max_it 50 \n"
46 "-snes_atol 1e-8 \n"
47 "-snes_rtol 1e-10 \n"
48 "-snes_monitor \n"
49 "-ksp_monitor \n"
50 "-snes_converged_reason \n"
51 "-my_order 1 \n"
52 "-my_order_lambda 1 \n"
53 "-my_order_contact 2 \n"
54 "-my_ho_levels_num 1 \n"
55 "-my_step_num 1 \n"
56 "-my_cn_value 1. \n"
57 "-my_r_value 1. \n"
58 "-my_alm_flag 0 \n";
59
60 string param_file = "param_file.petsc";
61 if (!static_cast<bool>(ifstream(param_file))) {
62 std::ofstream file(param_file.c_str(), std::ios::ate);
63 if (file.is_open()) {
64 file << default_options;
65 file.close();
66 }
67 }
68
71 FOUR_SEASONS = 2,
75 PLANE_AXI = 6,
76 ARC_THREE_SURF = 7,
82 };
83
84
86
87
88 moab::Core mb_instance;
89 moab::Interface &moab = mb_instance;
90
91 try {
92 PetscBool flg_file;
93
96 PetscInt order_contact = 1;
97 PetscInt nb_ho_levels = 0;
98 PetscInt order_lambda = 1;
99 PetscReal r_value = 1.;
100 PetscReal cn_value = -1;
101 PetscInt nb_sub_steps = 1;
102 PetscBool is_partitioned = PETSC_FALSE;
103 PetscBool is_newton_cotes = PETSC_FALSE;
104 PetscInt test_num = 0;
105 PetscBool convect_pts = PETSC_FALSE;
106 PetscBool print_contact_state = PETSC_FALSE;
107 PetscBool alm_flag = PETSC_FALSE;
108 PetscBool wave_surf_flag = PETSC_FALSE;
109 PetscInt wave_dim = 2;
110 PetscInt wave_surf_block_id = 1;
111 PetscReal wave_length = 1.0;
112 PetscReal wave_ampl = 0.01;
113 PetscReal mesh_height = 1.0;
114
115 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Elastic Config",
"none");
116
117 CHKERR PetscOptionsString(
"-my_file",
"mesh file name",
"",
"mesh.h5m",
119
120 CHKERR PetscOptionsInt(
"-my_order",
121 "approximation order of spatial positions", "", 1,
124 "-my_order_contact",
125 "approximation order of spatial positions in contact interface", "", 1,
126 &order_contact, PETSC_NULL);
127 CHKERR PetscOptionsInt(
"-my_ho_levels_num",
"number of higher order levels",
128 "", 0, &nb_ho_levels, PETSC_NULL);
129 CHKERR PetscOptionsInt(
"-my_order_lambda",
130 "approximation order of Lagrange multipliers", "", 1,
131 &order_lambda, PETSC_NULL);
132
133 CHKERR PetscOptionsInt(
"-my_step_num",
"number of steps",
"", nb_sub_steps,
134 &nb_sub_steps, PETSC_NULL);
135
136 CHKERR PetscOptionsBool(
"-my_is_partitioned",
137 "set if mesh is partitioned (this result that each "
138 "process keeps only part of the mes",
139 "", PETSC_FALSE, &is_partitioned, PETSC_NULL);
140
141 CHKERR PetscOptionsReal(
"-my_cn_value",
"default regularisation cn value",
142 "", 1., &cn_value, PETSC_NULL);
143
144 CHKERR PetscOptionsBool(
"-my_is_newton_cotes",
145 "set if Newton-Cotes quadrature rules are used", "",
146 PETSC_FALSE, &is_newton_cotes, PETSC_NULL);
147 CHKERR PetscOptionsInt(
"-my_test_num",
"test number",
"", 0, &test_num,
148 PETSC_NULL);
149 CHKERR PetscOptionsBool(
"-my_convect",
"set to convect integration pts",
"",
150 PETSC_FALSE, &convect_pts, PETSC_NULL);
151 CHKERR PetscOptionsBool(
"-my_print_contact_state",
152 "output number of active gp at every iteration", "",
153 PETSC_FALSE, &print_contact_state, PETSC_NULL);
154 CHKERR PetscOptionsBool(
"-my_alm_flag",
155 "if set use ALM, if not use C-function", "",
156 PETSC_FALSE, &alm_flag, PETSC_NULL);
157
158 CHKERR PetscOptionsBool(
"-my_wave_surf",
159 "if set true, make one of the surfaces wavy", "",
160 PETSC_FALSE, &wave_surf_flag, PETSC_NULL);
161 CHKERR PetscOptionsInt(
"-my_wave_surf_block_id",
162 "make wavy surface of the block with this id", "",
163 wave_surf_block_id, &wave_surf_block_id, PETSC_NULL);
164 CHKERR PetscOptionsInt(
"-my_wave_dim",
"dimension (2 or 3)",
"", wave_dim,
165 &wave_dim, PETSC_NULL);
166 CHKERR PetscOptionsReal(
"-my_wave_length",
"profile wavelength",
"",
167 wave_length, &wave_length, PETSC_NULL);
168 CHKERR PetscOptionsReal(
"-my_wave_ampl",
"profile amplitude",
"", wave_ampl,
169 &wave_ampl, PETSC_NULL);
170 CHKERR PetscOptionsReal(
"-my_mesh_height",
171 "vertical dimension of the mesh ", "", mesh_height,
172 &mesh_height, PETSC_NULL);
173
174 ierr = PetscOptionsEnd();
176
177
178 if (flg_file != PETSC_TRUE) {
179 SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
180 }
181
182 if (is_partitioned == PETSC_TRUE) {
184 "Partitioned mesh is not supported");
185 }
186
187 const char *option;
188 option = "";
190
191
194
199
201
202 auto add_prism_interface = [&](
Range &contact_prisms,
Range &master_tris,
204 std::vector<BitRefLevel> &bit_levels) {
208
210 if (cit->getName().compare(0, 11, "INT_CONTACT") == 0) {
211 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Insert %s (id: %d)\n",
212 cit->getName().c_str(), cit->getMeshsetId());
214
215
217 CHKERR moab.create_meshset(MESHSET_SET, ref_level_meshset);
219 ->getEntitiesByTypeAndRefLevel(bit_levels.back(),
221 ref_level_meshset);
223 ->getEntitiesByTypeAndRefLevel(bit_levels.back(),
225 ref_level_meshset);
226
227
228 CHKERR interface->
getSides(cubit_meshset, bit_levels.back(),
true, 0);
229
230 bit_levels.push_back(
BitRefLevel().set(bit_levels.size()));
231
233 cubit_meshset, true, true, 0);
234
235 CHKERR moab.delete_entities(&ref_level_meshset, 1);
236
237
241 ->updateMeshsetByEntitiesChildren(
242 cubit_meshset, bit_levels.back(), cubit_meshset, MBVERTEX,
243 true);
245 ->updateMeshsetByEntitiesChildren(cubit_meshset,
246 bit_levels.back(),
247 cubit_meshset, MBEDGE, true);
249 ->updateMeshsetByEntitiesChildren(cubit_meshset,
250 bit_levels.back(),
251 cubit_meshset, MBTRI, true);
253 ->updateMeshsetByEntitiesChildren(cubit_meshset,
254 bit_levels.back(),
255 cubit_meshset, MBTET, true);
256 }
257
259 bit_levels.pop_back();
260 }
261 }
262
264 CHKERR moab.create_meshset(MESHSET_SET, meshset_prisms);
266 ->getEntitiesByTypeAndRefLevel(bit_levels.back(),
BitRefLevel().set(),
267 MBPRISM, meshset_prisms);
268 CHKERR moab.get_entities_by_handle(meshset_prisms, contact_prisms);
269 CHKERR moab.delete_entities(&meshset_prisms, 1);
270
272 for (Range::iterator pit = contact_prisms.begin();
273 pit != contact_prisms.end(); pit++) {
274 CHKERR moab.side_element(*pit, 2, 3, tri);
275 master_tris.insert(tri);
276 CHKERR moab.side_element(*pit, 2, 4, tri);
277 slave_tris.insert(tri);
278 }
279
281 };
282
283 auto make_wavy_surface = [&](
int block_id,
int dim,
double lambda,
284 double delta,
double height) {
286 Range all_tets, all_nodes;
288 if (
bit->getName().compare(0, 11,
"MAT_ELASTIC") == 0) {
289 const int id =
bit->getMeshsetId();
291 if (id == block_id) {
293 bit->getMeshset(), 3, tets,
true);
294 all_tets.merge(tets);
295 }
296 }
297 }
299 double coords[3];
300 for (Range::iterator nit = all_nodes.begin(); nit != all_nodes.end();
301 nit++) {
302 CHKERR moab.get_coords(&*nit, 1, coords);
303 double x = coords[0];
304 double y = coords[1];
305 double z = coords[2];
306 double coef = (height + z) / height;
307 switch (dim) {
308 case 2:
309 coords[2] -= coef *
delta * (1. - cos(2. * M_PI * x /
lambda));
310 break;
311 case 3:
312 coords[2] -=
314 (1. - cos(2. * M_PI * x /
lambda) * cos(2. * M_PI * y /
lambda));
315 break;
316 default:
318 "Wrong dimension = %d", dim);
319 }
320
321 CHKERR moab.set_coords(&*nit, 1, coords);
322 }
324 };
325
326 auto set_contact_order = [&](
Range &contact_prisms,
int order_contact,
327 int nb_ho_levels) {
329 Range contact_tris, contact_edges;
330 CHKERR moab.get_adjacencies(contact_prisms, 2,
false, contact_tris,
331 moab::Interface::UNION);
332 contact_tris = contact_tris.subset_by_type(MBTRI);
333 CHKERR moab.get_adjacencies(contact_tris, 1,
false, contact_edges,
334 moab::Interface::UNION);
336 ho_ents.merge(contact_tris);
337 ho_ents.merge(contact_edges);
338 for (int ll = 0; ll < nb_ho_levels; ll++) {
339 Range ents, verts, tets;
340 CHKERR moab.get_connectivity(ho_ents, verts,
true);
341 CHKERR moab.get_adjacencies(verts, 3,
false, tets,
342 moab::Interface::UNION);
343 tets = tets.subset_by_type(MBTET);
344 for (auto d : {1, 2}) {
345 CHKERR moab.get_adjacencies(tets, d,
false, ents,
346 moab::Interface::UNION);
347 }
348 ho_ents = unite(ho_ents, ents);
349 ho_ents = unite(ho_ents, tets);
350 }
351
353 order_contact);
354
356 };
357
358 Range contact_prisms, master_tris, slave_tris;
359 std::vector<BitRefLevel> bit_levels;
360
363 0, 3, bit_levels.back());
364
365 CHKERR add_prism_interface(contact_prisms, master_tris, slave_tris,
366 bit_levels);
367
368 if (wave_surf_flag) {
369 CHKERR make_wavy_surface(wave_surf_block_id, wave_dim, wave_length,
370 wave_ampl, mesh_height);
371 }
372
375
378
381
387
388
394
399
400 if (order_contact >
order) {
401 CHKERR set_contact_order(contact_prisms, order_contact, nb_ho_levels);
402 }
403
404
406
407
408 {
411 }
412
413 {
416 }
417
418
419 boost::shared_ptr<Hooke<adouble>> hooke_adouble_ptr(
new Hooke<adouble>());
420 boost::shared_ptr<Hooke<double>> hooke_double_ptr(
new Hooke<double>());
422 CHKERR elastic.setBlocks(hooke_double_ptr, hooke_adouble_ptr);
423 CHKERR elastic.addElement(
"ELASTIC",
"SPATIAL_POSITION");
424
425 CHKERR elastic.setOperators(
"SPATIAL_POSITION",
"MESH_NODE_POSITIONS",
426 false, false);
427
428 auto make_contact_element = [&]() {
429 return boost::make_shared<SimpleContactProblem::SimpleContactElement>(
430 m_field);
431 };
432
433 auto make_convective_master_element = [&]() {
434 return boost::make_shared<
436 m_field, "SPATIAL_POSITION", "MESH_NODE_POSITIONS");
437 };
438
439 auto make_convective_slave_element = [&]() {
440 return boost::make_shared<
442 m_field, "SPATIAL_POSITION", "MESH_NODE_POSITIONS");
443 };
444
445 auto make_contact_common_data = [&]() {
446 return boost::make_shared<SimpleContactProblem::CommonDataSimpleContact>(
447 m_field);
448 };
449
450 auto get_contact_rhs = [&](auto contact_problem, auto make_element,
451 bool is_alm = false) {
452 auto fe_rhs_simple_contact = make_element();
453 auto common_data_simple_contact = make_contact_common_data();
454 if (print_contact_state) {
455 fe_rhs_simple_contact->contactStateVec =
456 common_data_simple_contact->gaussPtsStateVec;
457 }
458 contact_problem->setContactOperatorsRhs(
459 fe_rhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
460 "LAGMULT", is_alm);
461 return fe_rhs_simple_contact;
462 };
463
464 auto get_master_traction_rhs = [&](auto contact_problem, auto make_element,
465 bool is_alm = false) {
466 auto fe_rhs_simple_contact = make_element();
467 auto common_data_simple_contact = make_contact_common_data();
468 contact_problem->setMasterForceOperatorsRhs(
469 fe_rhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
470 "LAGMULT", is_alm);
471 return fe_rhs_simple_contact;
472 };
473
474 auto get_master_traction_lhs = [&](auto contact_problem, auto make_element,
475 bool is_alm = false) {
476 auto fe_lhs_simple_contact = make_element();
477 auto common_data_simple_contact = make_contact_common_data();
478 contact_problem->setMasterForceOperatorsLhs(
479 fe_lhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
480 "LAGMULT", is_alm);
481 return fe_lhs_simple_contact;
482 };
483
484 auto get_contact_lhs = [&](auto contact_problem, auto make_element,
485 bool is_alm = false) {
486 auto fe_lhs_simple_contact = make_element();
487 auto common_data_simple_contact = make_contact_common_data();
488 contact_problem->setContactOperatorsLhs(
489 fe_lhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
490 "LAGMULT", is_alm);
491 return fe_lhs_simple_contact;
492 };
493
494 auto cn_value_ptr = boost::make_shared<double>(cn_value);
495 auto contact_problem = boost::make_shared<SimpleContactProblem>(
496 m_field, cn_value_ptr, is_newton_cotes);
497
498
499 contact_problem->addContactElement("CONTACT_ELEM", "SPATIAL_POSITION",
500 "LAGMULT", contact_prisms);
501
502 contact_problem->addPostProcContactElement(
503 "CONTACT_POST_PROC", "SPATIAL_POSITION", "LAGMULT",
504 "MESH_NODE_POSITIONS", slave_tris);
505
507
508
510 "MESH_NODE_POSITIONS");
511
512
514
515
517
518
520
521
523 bit_levels.back());
524
525 DMType dm_name = "DMMOFEM";
527
530
531
532 CHKERR DMSetType(dm, dm_name);
533
534
536 CHKERR DMSetFromOptions(dm);
538
544
546
547
550
551
553
556 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
557 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
558
560 CHKERR VecGhostUpdateBegin(
F, INSERT_VALUES, SCATTER_FORWARD);
561 CHKERR VecGhostUpdateEnd(
F, INSERT_VALUES, SCATTER_FORWARD);
562
563 CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
564 CHKERR MatZeroEntries(Aij);
565
566
567 boost::shared_ptr<FEMethod> dirichlet_bc_ptr =
569 m_field,
"SPATIAL_POSITION", Aij,
D,
F));
570
572 dirichlet_bc_ptr->snes_x =
D;
573
574
575 boost::ptr_map<std::string, NeumannForcesSurface> neumann_forces;
577 m_field, neumann_forces, NULL, "SPATIAL_POSITION");
578
579 boost::ptr_map<std::string, NeumannForcesSurface>::iterator mit =
580 neumann_forces.begin();
581 for (; mit != neumann_forces.end(); mit++) {
584 &mit->second->getLoopFe(), NULL, NULL);
585 }
586
587
588
589 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr(
591 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr(
593
595 m_field, fe_spring_lhs_ptr, fe_spring_rhs_ptr, "SPATIAL_POSITION",
596 "MESH_NODE_POSITIONS");
597
599 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
600 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
603 dirichlet_bc_ptr.get(), NULL);
604 if (convect_pts == PETSC_TRUE) {
606 dm, "CONTACT_ELEM",
607 get_contact_rhs(contact_problem, make_convective_master_element),
608 PETSC_NULL, PETSC_NULL);
610 dm, "CONTACT_ELEM",
611 get_master_traction_rhs(contact_problem,
612 make_convective_slave_element),
613 PETSC_NULL, PETSC_NULL);
614 } else {
616 dm, "CONTACT_ELEM",
617 get_contact_rhs(contact_problem, make_contact_element, alm_flag),
618 PETSC_NULL, PETSC_NULL);
620 dm, "CONTACT_ELEM",
621 get_master_traction_rhs(contact_problem, make_contact_element,
622 alm_flag),
623 PETSC_NULL, PETSC_NULL);
624 }
625
627 PETSC_NULL, PETSC_NULL);
629 PETSC_NULL);
631 dirichlet_bc_ptr.get());
632
633 boost::shared_ptr<FEMethod> fe_null;
635 fe_null);
636 if (convect_pts == PETSC_TRUE) {
638 dm, "CONTACT_ELEM",
639 get_contact_lhs(contact_problem, make_convective_master_element),
640 PETSC_NULL, PETSC_NULL);
642 dm, "CONTACT_ELEM",
643 get_master_traction_lhs(contact_problem,
644 make_convective_slave_element),
645 PETSC_NULL, PETSC_NULL);
646 } else {
648 dm, "CONTACT_ELEM",
649 get_contact_lhs(contact_problem, make_contact_element, alm_flag),
650 PETSC_NULL, PETSC_NULL);
652 dm, "CONTACT_ELEM",
653 get_master_traction_lhs(contact_problem, make_contact_element,
654 alm_flag),
655 PETSC_NULL, PETSC_NULL);
656 }
658 PETSC_NULL, PETSC_NULL);
660 PETSC_NULL);
662 dirichlet_bc_ptr);
663
664 if (test_num) {
665 char testing_options[] = "-ksp_type fgmres "
666 "-pc_type lu "
667 "-pc_factor_mat_solver_type mumps "
668 "-snes_type newtonls "
669 "-snes_linesearch_type basic "
670 "-snes_max_it 10 "
671 "-snes_atol 1e-8 "
672 "-snes_rtol 1e-8 ";
673 CHKERR PetscOptionsInsertString(PETSC_NULL, testing_options);
674 }
675
677 CHKERR SNESSetDM(snes, dm);
678 SNESConvergedReason snes_reason;
680
681 {
682 CHKERR SNESSetDM(snes, dm);
686 CHKERR SNESSetFromOptions(snes);
687 }
688
690
691 CHKERR post_proc.generateReferenceElementMesh();
692 CHKERR post_proc.addFieldValuesPostProc(
"SPATIAL_POSITION");
693 CHKERR post_proc.addFieldValuesPostProc(
"MESH_NODE_POSITIONS");
694 CHKERR post_proc.addFieldValuesGradientPostProc(
"SPATIAL_POSITION");
695
696 std::map<int, NonlinearElasticElement::BlockData>::iterator sit =
697 elastic.setOfBlocks.begin();
698 for (; sit != elastic.setOfBlocks.end(); sit++) {
700 post_proc.postProcMesh, post_proc.mapGaussPts, "SPATIAL_POSITION",
701 sit->second, post_proc.commonData));
702 }
703
704 for (int ss = 0; ss != nb_sub_steps; ++ss) {
706 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Load scale: %6.4e\n",
708
709 CHKERR SNESSolve(snes, PETSC_NULL,
D);
710
711 CHKERR SNESGetConvergedReason(snes, &snes_reason);
712
713 int its;
714 CHKERR SNESGetIterationNumber(snes, &its);
715 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Number of Newton iterations = %D\n",
716 its);
717
718 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
719 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
720 }
721
722
723 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
724 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
726
727 PetscPrintf(PETSC_COMM_WORLD, "Loop post proc\n");
729
731 elastic.getLoopFeEnergy().eNergy = 0;
732 PetscPrintf(PETSC_COMM_WORLD, "Loop energy\n");
734
735 PetscPrintf(PETSC_COMM_WORLD, "Elastic energy %9.9f\n",
736 elastic.getLoopFeEnergy().eNergy);
737
738 {
740 std::ostringstream stm;
741 stm << "out"
742 << ".h5m";
745 PetscPrintf(PETSC_COMM_WORLD,
"Write file %s\n",
out_file_name.c_str());
747 "PARALLEL=WRITE_PART");
748 }
749
750
751 moab::Core mb_post;
752 moab::Interface &moab_proc = mb_post;
753
754 auto common_data_simple_contact = make_contact_common_data();
755
756 boost::shared_ptr<SimpleContactProblem::SimpleContactElement>
757 fe_post_proc_simple_contact;
758 if (convect_pts == PETSC_TRUE) {
759 fe_post_proc_simple_contact = make_convective_master_element();
760 } else {
761 fe_post_proc_simple_contact = make_contact_element();
762 }
763
764 contact_problem->setContactOperatorsForPostProc(
765 fe_post_proc_simple_contact, common_data_simple_contact, m_field,
766 "SPATIAL_POSITION", "LAGMULT", mb_post, alm_flag);
767
768 mb_post.delete_mesh();
769
770 CHKERR VecZeroEntries(common_data_simple_contact->gaussPtsStateVec);
771 CHKERR VecZeroEntries(common_data_simple_contact->contactAreaVec);
772
774 fe_post_proc_simple_contact);
775
776 std::array<double, 2> nb_gauss_pts;
777 std::array<double, 2> contact_area;
778
779 auto get_contact_data = [&](auto vec, std::array<double, 2> &data) {
781 CHKERR VecAssemblyBegin(vec);
782 CHKERR VecAssemblyEnd(vec);
783 const double *array;
784 CHKERR VecGetArrayRead(vec, &array);
788 }
789 CHKERR VecRestoreArrayRead(vec, &array);
791 };
792
793 CHKERR get_contact_data(common_data_simple_contact->gaussPtsStateVec,
794 nb_gauss_pts);
795 CHKERR get_contact_data(common_data_simple_contact->contactAreaVec,
796 contact_area);
797
799 PetscPrintf(PETSC_COMM_SELF, "Active gauss pts: %d out of %d\n",
800 (int)nb_gauss_pts[0], (int)nb_gauss_pts[1]);
801
802 PetscPrintf(PETSC_COMM_SELF, "Active contact area: %9.9f out of %9.9f\n",
803 contact_area[0], contact_area[1]);
804 }
805
806 if (test_num) {
807 double expected_energy, expected_contact_area;
808 int expected_nb_gauss_pts;
809 constexpr double eps = 1e-6;
810 switch (test_num) {
812 expected_energy = 3.0e-04;
813 expected_contact_area = 3.0;
814 expected_nb_gauss_pts = 576;
815 break;
816 case FOUR_SEASONS:
817 expected_energy = 1.2e-01;
818 expected_contact_area = 106.799036701;
819 expected_nb_gauss_pts = 672;
820 break;
822 expected_energy = 3.0e-04;
823 expected_contact_area = 1.75;
824 expected_nb_gauss_pts = 336;
825 break;
827 expected_energy = 3.125e-04;
828 expected_contact_area = 0.25;
829 expected_nb_gauss_pts = 84;
830 break;
832 expected_energy = 0.000096432;
833 expected_contact_area = 0.25;
834 expected_nb_gauss_pts = 336;
835 break;
836 case PLANE_AXI:
837 expected_energy = 0.000438889;
838 expected_contact_area = 0.784409608;
839 expected_nb_gauss_pts = 300;
840 break;
841 case ARC_THREE_SURF:
842 expected_energy = 0.002573411;
843 expected_contact_area = 2.831455633;
844 expected_nb_gauss_pts = 228;
845 break;
847 expected_energy = 0.000733553;
848 expected_contact_area = 3.0;
849 expected_nb_gauss_pts = 144;
850 break;
852 expected_energy = 0.000733621;
853 expected_contact_area = 3.0;
854 expected_nb_gauss_pts = 144;
855 break;
857 expected_energy = 0.008537863;
858 expected_contact_area = 0.125;
859 expected_nb_gauss_pts = 384;
860 break;
862 expected_energy = 0.008538894;
863 expected_contact_area = 0.125;
864 expected_nb_gauss_pts = 384;
865 break;
866 default:
868 "Unknown test number %d", test_num);
869 }
870 if (std::abs(elastic.getLoopFeEnergy().eNergy - expected_energy) >
eps) {
872 "Wrong energy %6.4e != %6.4e", expected_energy,
873 elastic.getLoopFeEnergy().eNergy);
874 }
876 if ((int)nb_gauss_pts[0] != expected_nb_gauss_pts) {
878 "Wrong number of active gauss pts %d != %d",
879 expected_nb_gauss_pts, (int)nb_gauss_pts[0]);
880 }
881 if (std::abs(contact_area[0] - expected_contact_area) >
eps) {
883 "Wrong active contact area %6.4e != %6.4e",
884 expected_contact_area, contact_area[0]);
885 }
886 }
887 }
888
889 {
891 std::ostringstream strm;
892 strm << "out_contact_integ_pts"
893 << ".h5m";
895 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Write file %s\n",
898 "PARALLEL=WRITE_PART");
899 }
900
901 boost::shared_ptr<PostProcFaceOnRefinedMesh> post_proc_contact_ptr(
903
905 false);
906 CHKERR post_proc_contact_ptr->generateReferenceElementMesh();
907 CHKERR post_proc_contact_ptr->addFieldValuesPostProc(
"LAGMULT");
908 CHKERR post_proc_contact_ptr->addFieldValuesPostProc(
"SPATIAL_POSITION");
909 CHKERR post_proc_contact_ptr->addFieldValuesPostProc(
"MESH_NODE_POSITIONS");
910
912 post_proc_contact_ptr);
913
914 {
916 std::ostringstream stm;
917 stm << "out_contact"
918 << ".h5m";
920 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Write file %s\n",
922 CHKERR post_proc_contact_ptr->postProcMesh.write_file(
924 }
925 }
927
928
930
931 return 0;
932}
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_ATOM_TEST_INVALID
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
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 DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
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.
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
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
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)
static constexpr double delta
Set Dirichlet boundary conditions on spatial displacements.
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
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 MPI_Comm & get_comm() const =0
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.
Interface for managing meshsets containing materials and boundary conditions.
MoFEMErrorCode printForceSet() const
print meshsets with force boundary conditions data structure
MoFEMErrorCode printMaterialsSet() const
print meshsets with material data structure set on it
MoFEMErrorCode printDisplacementSet() const
print meshsets with displacement boundary conditions data structure
Create interface from given surface and insert flat prisms in-between.
MoFEMErrorCode getSides(const int msId, const CubitBCType cubit_bc_type, const BitRefLevel mesh_bit_level, const bool recursive, int verb=QUIET)
Store tetrahedra from each side of the interface separately in two child meshsets of the parent meshs...
MoFEMErrorCode splitSides(const EntityHandle meshset, const BitRefLevel &bit, const int msId, const CubitBCType cubit_bc_type, const bool add_interface_entities, const bool recursive=false, int verb=QUIET)
Split nodes and other entities of tetrahedra on both sides of the interface and insert flat prisms in...
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 reference to pointer of interface.
structure grouping operators and data used for calculation of nonlinear elastic element