41 const string default_options =
"-ksp_type fgmres \n"
43 "-pc_factor_mat_solver_type mumps \n"
44 "-mat_mumps_icntl_13 1 \n"
45 "-mat_mumps_icntl_14 800 \n"
46 "-mat_mumps_icntl_20 0 \n"
47 "-mat_mumps_icntl_24 1 \n"
48 "-snes_type newtonls \n"
49 "-snes_linesearch_type basic \n"
50 "-snes_divergence_tolerance 0 \n"
56 "-snes_converged_reason \n"
58 "-my_order_lambda 1 \n"
59 "-my_order_contact 2 \n"
60 "-my_ho_levels_num 1 \n"
65 "-my_eigen_pos_flag 0 \n";
67 string param_file =
"param_file.petsc";
68 if (!
static_cast<bool>(ifstream(param_file))) {
69 std::ofstream file(param_file.c_str(), std::ios::ate);
71 file << default_options;
85 PetscBool flg_file_out;
88 char output_mesh_name[255];
90 PetscInt order_contact = 1;
91 PetscInt nb_ho_levels = 0;
92 PetscInt order_lambda = 1;
93 PetscReal r_value = 1.;
94 PetscReal cn_value = -1;
95 PetscInt nb_sub_steps = 1;
96 PetscBool is_partitioned = PETSC_FALSE;
97 PetscBool is_newton_cotes = PETSC_FALSE;
98 PetscInt test_num = 0;
99 PetscBool convect_pts = PETSC_FALSE;
100 PetscBool print_contact_state = PETSC_FALSE;
101 PetscBool alm_flag = PETSC_FALSE;
102 PetscBool eigen_pos_flag = PETSC_FALSE;
104 PetscReal thermal_expansion_coef = 1e-5;
106 PetscReal scale_factor = 1.0;
107 PetscBool analytical_input = PETSC_TRUE;
108 char stress_tag_name[255];
109 PetscBool flg_tag_name;
110 PetscBool save_mean_stress = PETSC_TRUE;
111 PetscBool ignore_contact = PETSC_FALSE;
112 PetscBool ignore_pressure = PETSC_FALSE;
114 PetscBool deform_flat_flag = PETSC_FALSE;
115 PetscReal flat_shift = 1.0;
116 PetscReal mesh_height = 1.0;
118 PetscBool wave_surf_flag = PETSC_FALSE;
119 PetscInt wave_dim = 2;
120 PetscReal wave_length = 1.0;
121 PetscReal wave_ampl = 0.01;
123 PetscBool delete_prisms = PETSC_FALSE;
125 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Elastic Config",
"none");
127 CHKERR PetscOptionsString(
"-my_file",
"mesh file name",
"",
"mesh.h5m",
129 CHKERR PetscOptionsString(
"-my_output_mesh_file",
"output mesh file name",
130 "",
"mesh.h5m", output_mesh_name, 255,
133 CHKERR PetscOptionsInt(
"-my_order",
134 "approximation order of spatial positions",
"", 1,
138 "approximation order of spatial positions in contact interface",
"", 1,
139 &order_contact, PETSC_NULL);
140 CHKERR PetscOptionsInt(
"-my_ho_levels_num",
"number of higher order levels",
141 "", 0, &nb_ho_levels, PETSC_NULL);
142 CHKERR PetscOptionsInt(
"-my_order_lambda",
143 "approximation order of Lagrange multipliers",
"", 1,
144 &order_lambda, PETSC_NULL);
146 CHKERR PetscOptionsInt(
"-my_step_num",
"number of steps",
"", nb_sub_steps,
147 &nb_sub_steps, PETSC_NULL);
149 CHKERR PetscOptionsBool(
"-my_is_partitioned",
150 "set if mesh is partitioned (this result that each "
151 "process keeps only part of the mes",
152 "", PETSC_FALSE, &is_partitioned, PETSC_NULL);
154 CHKERR PetscOptionsReal(
"-my_cn_value",
"default regularisation cn value",
155 "", 1., &cn_value, PETSC_NULL);
157 CHKERR PetscOptionsBool(
"-my_is_newton_cotes",
158 "set if Newton-Cotes quadrature rules are used",
"",
159 PETSC_FALSE, &is_newton_cotes, PETSC_NULL);
160 CHKERR PetscOptionsInt(
"-my_test_num",
"test number",
"", 0, &test_num,
162 CHKERR PetscOptionsBool(
"-my_convect",
"set to convect integration pts",
"",
163 PETSC_FALSE, &convect_pts, PETSC_NULL);
164 CHKERR PetscOptionsBool(
"-my_print_contact_state",
165 "output number of active gp at every iteration",
"",
166 PETSC_FALSE, &print_contact_state, PETSC_NULL);
167 CHKERR PetscOptionsBool(
"-my_alm_flag",
168 "if set use ALM, if not use C-function",
"",
169 PETSC_FALSE, &alm_flag, PETSC_NULL);
170 CHKERR PetscOptionsBool(
"-my_eigen_pos_flag",
171 "if set use eigen spatial positions are taken into "
172 "account for predeformed configuration",
173 "", PETSC_FALSE, &eigen_pos_flag, PETSC_NULL);
175 CHKERR PetscOptionsReal(
"-my_scale_factor",
"scale factor",
"",
176 scale_factor, &scale_factor, PETSC_NULL);
178 CHKERR PetscOptionsBool(
"-my_ignore_contact",
"if set true, ignore contact",
179 "", PETSC_FALSE, &ignore_contact, PETSC_NULL);
180 CHKERR PetscOptionsBool(
"-my_ignore_pressure",
181 "if set true, ignore pressure",
"", PETSC_FALSE,
182 &ignore_pressure, PETSC_NULL);
183 CHKERR PetscOptionsBool(
"-my_analytical_input",
184 "if set true, use analytical strain",
"",
185 PETSC_TRUE, &analytical_input, PETSC_NULL);
186 CHKERR PetscOptionsBool(
"-my_save_mean_stress",
187 "if set true, save mean stress",
"", PETSC_TRUE,
188 &save_mean_stress, PETSC_NULL);
189 CHKERR PetscOptionsString(
190 "-my_stress_tag_name",
"stress tag name file name",
"",
191 "INTERNAL_STRESS", stress_tag_name, 255, &flg_tag_name);
193 "-my_thermal_expansion_coef",
"thermal expansion coef ",
"",
194 thermal_expansion_coef, &thermal_expansion_coef, PETSC_NULL);
198 CHKERR PetscOptionsReal(
"-my_mesh_height",
199 "vertical dimension of the mesh ",
"", mesh_height,
200 &mesh_height, PETSC_NULL);
201 CHKERR PetscOptionsBool(
"-my_deform_flat",
"if set true, deform flat",
"",
202 PETSC_FALSE, &deform_flat_flag, PETSC_NULL);
203 CHKERR PetscOptionsReal(
"-my_flat_shift",
"flat shift ",
"", flat_shift,
204 &flat_shift, PETSC_NULL);
206 CHKERR PetscOptionsBool(
"-my_wave_surf",
207 "if set true, make one of the surfaces wavy",
"",
208 PETSC_FALSE, &wave_surf_flag, PETSC_NULL);
209 CHKERR PetscOptionsInt(
"-my_wave_dim",
"dimension (2 or 3)",
"", wave_dim,
210 &wave_dim, PETSC_NULL);
211 CHKERR PetscOptionsReal(
"-my_wave_length",
"profile wavelength",
"",
212 wave_length, &wave_length, PETSC_NULL);
213 CHKERR PetscOptionsReal(
"-my_wave_ampl",
"profile amplitude",
"", wave_ampl,
214 &wave_ampl, PETSC_NULL);
216 CHKERR PetscOptionsBool(
"-my_delete_prisms",
"if set true, delete prisms",
217 "", PETSC_FALSE, &delete_prisms, PETSC_NULL);
219 ierr = PetscOptionsEnd();
223 if (flg_file != PETSC_TRUE) {
224 SETERRQ(PETSC_COMM_SELF, 1,
"*** ERROR -my_file (MESH FILE NEEDED)");
227 if (is_partitioned == PETSC_TRUE) {
229 "Partitioned mesh is not supported");
244 std::vector<BitRefLevel> bit_levels;
247 CHKERR bit_ref_manager->setBitRefLevelByDim(0, 3, bit_levels.back());
249 auto add_prism_interface = [&](std::vector<BitRefLevel> &bit_levels) {
254 if (cit->getName().compare(0, 11,
"INT_CONTACT") == 0) {
255 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Insert %s (id: %d)\n",
256 cit->getName().c_str(), cit->getMeshsetId());
261 CHKERR moab.create_meshset(MESHSET_SET, ref_level_meshset);
262 CHKERR bit_ref_manager->getEntitiesByTypeAndRefLevel(
263 bit_levels.back(),
BitRefLevel().set(), MBTET, ref_level_meshset);
264 CHKERR bit_ref_manager->getEntitiesByTypeAndRefLevel(
269 CHKERR prism_interface->getSides(cubit_meshset, bit_levels.back(),
272 bit_levels.push_back(
BitRefLevel().set(bit_levels.size()));
274 CHKERR prism_interface->splitSides(ref_level_meshset,
275 bit_levels.back(), cubit_meshset,
278 CHKERR moab.delete_entities(&ref_level_meshset, 1);
283 CHKERR bit_ref_manager->updateMeshsetByEntitiesChildren(
284 cubit_meshset, bit_levels.back(), cubit_meshset, MBVERTEX,
286 CHKERR bit_ref_manager->updateMeshsetByEntitiesChildren(
287 cubit_meshset, bit_levels.back(), cubit_meshset, MBEDGE,
true);
288 CHKERR bit_ref_manager->updateMeshsetByEntitiesChildren(
289 cubit_meshset, bit_levels.back(), cubit_meshset, MBTRI,
true);
290 CHKERR bit_ref_manager->updateMeshsetByEntitiesChildren(
291 cubit_meshset, bit_levels.back(), cubit_meshset, MBTET,
true);
294 CHKERR bit_ref_manager->shiftRightBitRef(1);
295 bit_levels.pop_back();
302 auto find_contact_prisms = [&](std::vector<BitRefLevel> &bit_levels,
308 CHKERR moab.create_meshset(MESHSET_SET, meshset_prisms);
309 CHKERR bit_ref_manager->getEntitiesByTypeAndRefLevel(
310 bit_levels.back(),
BitRefLevel().set(), MBPRISM, meshset_prisms);
311 CHKERR moab.get_entities_by_handle(meshset_prisms, contact_prisms);
312 CHKERR moab.delete_entities(&meshset_prisms, 1);
315 for (Range::iterator pit = contact_prisms.begin();
316 pit != contact_prisms.end(); pit++) {
317 CHKERR moab.side_element(*pit, 2, 3, tri);
318 master_tris.insert(tri);
319 CHKERR moab.side_element(*pit, 2, 4, tri);
320 slave_tris.insert(tri);
326 Range contact_prisms, master_tris, slave_tris;
328 if (!ignore_contact) {
329 if (analytical_input) {
330 CHKERR add_prism_interface(bit_levels);
332 CHKERR find_contact_prisms(bit_levels, contact_prisms, master_tris,
336 auto deform_flat_surface = [&](
int block_id,
double shift,
double height) {
338 Range all_tets, all_nodes;
340 if (
bit->getName().compare(0, 11,
"MAT_ELASTIC") == 0) {
341 const int id =
bit->getMeshsetId();
343 if (
id == block_id) {
345 bit->getMeshset(), 3, tets,
true);
346 all_tets.merge(tets);
352 for (Range::iterator nit = all_nodes.begin(); nit != all_nodes.end();
354 CHKERR moab.get_coords(&*nit, 1, coords);
355 double x = coords[0];
356 double y = coords[1];
357 double z = coords[2];
359 CHKERR moab.set_coords(&*nit, 1, coords);
364 auto make_wavy_surface = [&](
int block_id,
int dim,
double lambda,
365 double delta,
double height) {
367 Range all_tets, all_nodes;
369 if (
bit->getName().compare(0, 11,
"MAT_ELASTIC") == 0) {
370 const int id =
bit->getMeshsetId();
372 if (
id == block_id) {
374 bit->getMeshset(), 3, tets,
true);
375 all_tets.merge(tets);
381 for (Range::iterator nit = all_nodes.begin(); nit != all_nodes.end();
383 CHKERR moab.get_coords(&*nit, 1, coords);
384 double x = coords[0];
385 double y = coords[1];
386 double z = coords[2];
387 double coef = (height + z) / height;
390 coords[2] -= coef *
delta * (1. - cos(2. * M_PI * x /
lambda));
395 (1. - cos(2. * M_PI * x /
lambda) * cos(2. * M_PI * y /
lambda));
399 "Wrong dimension = %d", dim);
402 CHKERR moab.set_coords(&*nit, 1, coords);
407 if (deform_flat_flag && analytical_input) {
408 CHKERR deform_flat_surface(1, flat_shift, mesh_height);
409 CHKERR deform_flat_surface(2, -flat_shift, mesh_height);
412 if (wave_surf_flag && analytical_input) {
413 CHKERR make_wavy_surface(1, wave_dim, wave_length, wave_ampl,
424 if (!eigen_pos_flag) {
445 if (!eigen_pos_flag) {
458 auto set_contact_order = [&](
Range &contact_prisms,
int order_contact,
461 Range contact_tris, contact_edges;
462 CHKERR moab.get_adjacencies(contact_prisms, 2,
false, contact_tris,
463 moab::Interface::UNION);
464 contact_tris = contact_tris.subset_by_type(MBTRI);
465 CHKERR moab.get_adjacencies(contact_tris, 1,
false, contact_edges,
466 moab::Interface::UNION);
468 ho_ents.merge(contact_tris);
469 ho_ents.merge(contact_edges);
470 for (
int ll = 0; ll < nb_ho_levels; ll++) {
471 Range ents, verts, tets;
472 CHKERR moab.get_connectivity(ho_ents, verts,
true);
473 CHKERR moab.get_adjacencies(verts, 3,
false, tets,
474 moab::Interface::UNION);
475 tets = tets.subset_by_type(MBTET);
476 for (
auto d : {1, 2}) {
477 CHKERR moab.get_adjacencies(tets,
d,
false, ents,
478 moab::Interface::UNION);
480 ho_ents = unite(ho_ents, ents);
481 ho_ents = unite(ho_ents, tets);
490 if (!ignore_contact && order_contact >
order) {
491 CHKERR set_contact_order(contact_prisms, order_contact, nb_ho_levels);
509 CHKERR moab.get_connectivity(slave_tris, slave_verts,
true);
511 slave_verts,
"LAGMULT");
514 boost::shared_ptr<std::map<int, BlockData>> block_sets_ptr =
515 boost::make_shared<std::map<int, BlockData>>();
518 boost::shared_ptr<ForcesAndSourcesCore> fe_elastic_lhs_ptr(
520 boost::shared_ptr<ForcesAndSourcesCore> fe_elastic_rhs_ptr(
522 fe_elastic_lhs_ptr->getRuleHook =
VolRule();
523 fe_elastic_rhs_ptr->getRuleHook =
VolRule();
527 "MESH_NODE_POSITIONS",
false);
529 auto data_hooke_element_at_pts =
530 boost::make_shared<HookeInternalStressElement::DataAtIntegrationPts>();
532 block_sets_ptr,
"SPATIAL_POSITION",
533 "MESH_NODE_POSITIONS",
false,
false,
534 MBTET, data_hooke_element_at_pts);
535 auto thermal_strain =
536 [&thermal_expansion_coef, &
init_temp, &test_num](
544 t_thermal_strain(
i,
j) = 0.0;
550 t_thermal_strain(
i,
j) =
555 t_thermal_strain(
i,
j) = -thermal_expansion_coef *
t_kd(
i,
j);
558 t_thermal_strain(2, 2) = thermal_expansion_coef;
561 t_thermal_strain(
i,
j) = thermal_expansion_coef *
t_kd(
i,
j);
566 return t_thermal_strain;
569 if (analytical_input) {
570 fe_elastic_rhs_ptr->getOpPtrVector().push_back(
571 new HookeElement::OpAnalyticalInternalStrain_dx<0>(
572 "SPATIAL_POSITION", data_hooke_element_at_pts, thermal_strain));
573 fe_elastic_rhs_ptr->getOpPtrVector().push_back(
575 "SPATIAL_POSITION",
"SPATIAL_POSITION", data_hooke_element_at_pts,
578 fe_elastic_rhs_ptr->getOpPtrVector().push_back(
580 "SPATIAL_POSITION",
"SPATIAL_POSITION", data_hooke_element_at_pts,
581 moab, stress_tag_name));
582 fe_elastic_rhs_ptr->getOpPtrVector().push_back(
584 "SPATIAL_POSITION", data_hooke_element_at_pts));
587 fe_elastic_rhs_ptr->getOpPtrVector().push_back(
589 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat));
590 fe_elastic_rhs_ptr->getOpPtrVector().push_back(
592 "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
593 fe_elastic_rhs_ptr->getOpPtrVector().push_back(
595 "SPATIAL_POSITION",
"SPATIAL_POSITION", data_hooke_element_at_pts,
596 *block_sets_ptr.get(), moab, scale_factor, save_mean_stress,
false,
601 if (
bit->getName().compare(0, 11,
"MAT_ELASTIC") == 0) {
603 CHKERR moab.get_entities_by_dimension(
bit->getMeshset(), 3, tets,
true);
604 all_tets.merge(tets);
607 Skinner skinner(&moab);
609 CHKERR skinner.find_skin(0, all_tets,
false, skin_tris);
620 "MESH_NODE_POSITIONS");
624 auto make_contact_element = [&]() {
625 return boost::make_shared<SimpleContactProblem::SimpleContactElement>(
629 auto make_convective_master_element = [&]() {
630 return boost::make_shared<
632 m_field,
"SPATIAL_POSITION",
"MESH_NODE_POSITIONS");
635 auto make_convective_slave_element = [&]() {
636 return boost::make_shared<
638 m_field,
"SPATIAL_POSITION",
"MESH_NODE_POSITIONS");
641 auto make_contact_common_data = [&]() {
642 return boost::make_shared<SimpleContactProblem::CommonDataSimpleContact>(
646 auto get_contact_rhs = [&](
auto contact_problem,
auto make_element,
647 bool is_alm =
false) {
648 auto fe_rhs_simple_contact = make_element();
649 auto common_data_simple_contact = make_contact_common_data();
650 if (print_contact_state) {
651 fe_rhs_simple_contact->contactStateVec =
652 common_data_simple_contact->gaussPtsStateVec;
654 contact_problem->setContactOperatorsRhs(
655 fe_rhs_simple_contact, common_data_simple_contact,
"SPATIAL_POSITION",
656 "LAGMULT", is_alm, eigen_pos_flag,
"EIGEN_POSITIONS",
true);
657 return fe_rhs_simple_contact;
660 auto get_master_traction_rhs = [&](
auto contact_problem,
auto make_element,
661 bool is_alm =
false) {
662 auto fe_rhs_simple_contact = make_element();
663 auto common_data_simple_contact = make_contact_common_data();
664 contact_problem->setMasterForceOperatorsRhs(
665 fe_rhs_simple_contact, common_data_simple_contact,
"SPATIAL_POSITION",
666 "LAGMULT", is_alm, eigen_pos_flag,
"EIGEN_POSITIONS",
true);
667 return fe_rhs_simple_contact;
670 auto get_master_traction_lhs = [&](
auto contact_problem,
auto make_element,
671 bool is_alm =
false) {
672 auto fe_lhs_simple_contact = make_element();
673 auto common_data_simple_contact = make_contact_common_data();
674 contact_problem->setMasterForceOperatorsLhs(
675 fe_lhs_simple_contact, common_data_simple_contact,
"SPATIAL_POSITION",
676 "LAGMULT", is_alm, eigen_pos_flag,
"EIGEN_POSITIONS",
true);
677 return fe_lhs_simple_contact;
680 auto get_contact_lhs = [&](
auto contact_problem,
auto make_element,
681 bool is_alm =
false) {
682 auto fe_lhs_simple_contact = make_element();
683 auto common_data_simple_contact = make_contact_common_data();
684 contact_problem->setContactOperatorsLhs(
685 fe_lhs_simple_contact, common_data_simple_contact,
"SPATIAL_POSITION",
686 "LAGMULT", is_alm, eigen_pos_flag,
"EIGEN_POSITIONS",
true);
687 return fe_lhs_simple_contact;
690 auto cn_value_ptr = boost::make_shared<double>(cn_value);
691 auto contact_problem = boost::make_shared<SimpleContactProblem>(
692 m_field, cn_value_ptr, is_newton_cotes);
696 contact_problem->addContactElement(
"CONTACT_ELEM",
"SPATIAL_POSITION",
697 "LAGMULT", contact_prisms);
699 contact_problem->addContactElement(
"CONTACT_ELEM",
"SPATIAL_POSITION",
700 "LAGMULT", contact_prisms,
701 eigen_pos_flag,
"EIGEN_POSITIONS");
703 contact_problem->addPostProcContactElement(
704 "CONTACT_POST_PROC",
"SPATIAL_POSITION",
"LAGMULT",
705 "MESH_NODE_POSITIONS", slave_tris);
711 "MESH_NODE_POSITIONS");
726 DMType dm_name =
"DMMOFEM";
733 CHKERR DMSetType(dm, dm_name);
737 CHKERR DMSetFromOptions(dm);
758 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
759 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
762 CHKERR VecGhostUpdateBegin(
F, INSERT_VALUES, SCATTER_FORWARD);
763 CHKERR VecGhostUpdateEnd(
F, INSERT_VALUES, SCATTER_FORWARD);
765 CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
766 CHKERR MatZeroEntries(Aij);
768 fe_elastic_rhs_ptr->snes_f =
F;
769 fe_elastic_lhs_ptr->snes_B = Aij;
772 boost::shared_ptr<FEMethod> dirichlet_bc_ptr =
774 m_field,
"SPATIAL_POSITION", Aij,
D,
F));
776 dirichlet_bc_ptr->snes_ctx = SnesMethod::CTX_SNESNONE;
777 dirichlet_bc_ptr->snes_x =
D;
780 boost::ptr_map<std::string, NeumannForcesSurface> neumann_forces;
782 m_field, neumann_forces, NULL,
"SPATIAL_POSITION");
784 boost::ptr_map<std::string, NeumannForcesSurface>::iterator mit =
785 neumann_forces.begin();
786 for (; mit != neumann_forces.end(); mit++) {
789 &mit->second->getLoopFe(), NULL, NULL);
794 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr(
796 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr(
800 m_field, fe_spring_lhs_ptr, fe_spring_rhs_ptr,
"SPATIAL_POSITION",
801 "MESH_NODE_POSITIONS");
804 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
805 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
808 dirichlet_bc_ptr.get(), NULL);
809 if (convect_pts == PETSC_TRUE) {
812 get_contact_rhs(contact_problem, make_convective_master_element),
813 PETSC_NULL, PETSC_NULL);
816 get_master_traction_rhs(contact_problem,
817 make_convective_slave_element),
818 PETSC_NULL, PETSC_NULL);
822 get_contact_rhs(contact_problem, make_contact_element, alm_flag),
823 PETSC_NULL, PETSC_NULL);
826 get_master_traction_rhs(contact_problem, make_contact_element,
828 PETSC_NULL, PETSC_NULL);
837 dirichlet_bc_ptr.get());
839 boost::shared_ptr<FEMethod> fe_null;
842 if (convect_pts == PETSC_TRUE) {
845 get_contact_lhs(contact_problem, make_convective_master_element),
846 PETSC_NULL, PETSC_NULL);
849 get_master_traction_lhs(contact_problem,
850 make_convective_slave_element),
851 PETSC_NULL, PETSC_NULL);
855 get_contact_lhs(contact_problem, make_contact_element, alm_flag),
856 PETSC_NULL, PETSC_NULL);
859 get_master_traction_lhs(contact_problem, make_contact_element,
861 PETSC_NULL, PETSC_NULL);
871 char testing_options[] =
"-ksp_type fgmres "
873 "-pc_factor_mat_solver_type mumps "
874 "-snes_type newtonls "
875 "-snes_linesearch_type basic "
879 CHKERR PetscOptionsInsertString(PETSC_NULL, testing_options);
883 CHKERR SNESSetDM(snes, dm);
887 CHKERR SNESSetDM(snes, dm);
891 CHKERR SNESSetFromOptions(snes);
896 CHKERR post_proc_skin.generateReferenceElementMesh();
898 CHKERR post_proc_skin.addFieldValuesPostProc(
"SPATIAL_POSITION");
899 CHKERR post_proc_skin.addFieldValuesPostProc(
"MESH_NODE_POSITIONS");
900 CHKERR post_proc_skin.addFieldValuesPostProc(
"EIGEN_POSITIONS");
902 struct OpGetFieldGradientValuesOnSkin
905 const std::string feVolName;
906 boost::shared_ptr<VolSideFe> sideOpFe;
908 OpGetFieldGradientValuesOnSkin(
const std::string
field_name,
909 const std::string vol_fe_name,
910 boost::shared_ptr<VolSideFe> side_fe)
913 feVolName(vol_fe_name), sideOpFe(side_fe) {}
918 if (
type != MBVERTEX)
925 boost::shared_ptr<VolSideFe> my_vol_side_fe_ptr =
926 boost::make_shared<VolSideFe>(m_field);
927 my_vol_side_fe_ptr->getOpPtrVector().push_back(
929 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat));
930 my_vol_side_fe_ptr->getOpPtrVector().push_back(
932 "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
934 post_proc_skin.getOpPtrVector().push_back(
935 new OpGetFieldGradientValuesOnSkin(
"SPATIAL_POSITION",
"ELASTIC",
936 my_vol_side_fe_ptr));
937 post_proc_skin.getOpPtrVector().push_back(
939 "SPATIAL_POSITION", data_hooke_element_at_pts->spatPosMat));
940 post_proc_skin.getOpPtrVector().push_back(
942 "MESH_NODE_POSITIONS", data_hooke_element_at_pts->meshNodePosMat));
944 post_proc_skin.getOpPtrVector().push_back(
947 "SPATIAL_POSITION", data_hooke_element_at_pts,
948 *block_sets_ptr.get(), post_proc_skin.postProcMesh,
949 post_proc_skin.mapGaussPts,
false,
false));
951 for (
int ss = 0; ss != nb_sub_steps; ++ss) {
952 if (!ignore_pressure) {
956 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Ignoring pressure...\n");
958 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Load scale: %6.4e\n",
961 CHKERR SNESSolve(snes, PETSC_NULL,
D);
963 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
964 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
972 "MESH_NODE_POSITIONS",
false,
false,
974 const double *eng_ptr;
975 CHKERR VecGetArrayRead(v_energy, &eng_ptr);
977 PetscPrintf(PETSC_COMM_WORLD,
"Elastic energy: %8.8e\n", *eng_ptr);
980 PetscPrintf(PETSC_COMM_WORLD,
"Loop post proc on the skin\n");
984 stm <<
"out_skin.h5m";
986 PetscPrintf(PETSC_COMM_WORLD,
"Write file %s\n",
out_file_name.c_str());
987 CHKERR post_proc_skin.writeFile(stm.str());
994 auto common_data_simple_contact = make_contact_common_data();
996 boost::shared_ptr<SimpleContactProblem::SimpleContactElement>
997 fe_post_proc_simple_contact;
998 if (convect_pts == PETSC_TRUE) {
999 fe_post_proc_simple_contact = make_convective_master_element();
1001 fe_post_proc_simple_contact = make_contact_element();
1004 std::array<double, 2> nb_gauss_pts;
1005 std::array<double, 2> contact_area;
1007 if (!ignore_contact) {
1008 contact_problem->setContactOperatorsForPostProc(
1009 fe_post_proc_simple_contact, common_data_simple_contact, m_field,
1010 "SPATIAL_POSITION",
"LAGMULT", mb_post, alm_flag, eigen_pos_flag,
1011 "EIGEN_POSITIONS",
true);
1013 mb_post.delete_mesh();
1015 CHKERR VecZeroEntries(common_data_simple_contact->gaussPtsStateVec);
1016 CHKERR VecZeroEntries(common_data_simple_contact->contactAreaVec);
1019 fe_post_proc_simple_contact);
1021 auto get_contact_data = [&](
auto vec, std::array<double, 2> &data) {
1023 CHKERR VecAssemblyBegin(vec);
1024 CHKERR VecAssemblyEnd(vec);
1025 const double *array;
1026 CHKERR VecGetArrayRead(vec, &array);
1028 for (
int i : {0, 1})
1031 CHKERR VecRestoreArrayRead(vec, &array);
1035 CHKERR get_contact_data(common_data_simple_contact->gaussPtsStateVec,
1037 CHKERR get_contact_data(common_data_simple_contact->contactAreaVec,
1041 PetscPrintf(PETSC_COMM_SELF,
"Active gauss pts: %d out of %d\n",
1042 (
int)nb_gauss_pts[0], (
int)nb_gauss_pts[1]);
1044 PetscPrintf(PETSC_COMM_SELF,
1045 "Active contact area: %8.8f out of %8.8f (%8.8f% %)\n",
1046 contact_area[0], contact_area[1],
1047 contact_area[0] / contact_area[1] * 100.);
1051 std::ostringstream strm;
1052 strm <<
"out_contact_integ_pts"
1055 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Write file %s\n",
1058 "PARALLEL=WRITE_PART");
1061 boost::shared_ptr<PostProcFaceOnRefinedMesh> post_proc_contact_ptr(
1064 CHKERR post_proc_contact_ptr->generateReferenceElementMesh();
1067 CHKERR post_proc_contact_ptr->addFieldValuesPostProc(
"LAGMULT");
1068 CHKERR post_proc_contact_ptr->addFieldValuesPostProc(
"SPATIAL_POSITION");
1069 CHKERR post_proc_contact_ptr->addFieldValuesPostProc(
"MESH_NODE_POSITIONS");
1071 if (!ignore_contact) {
1073 post_proc_contact_ptr);
1075 std::ostringstream stm;
1076 stm <<
"out_contact"
1079 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Write file %s\n",
1081 CHKERR post_proc_contact_ptr->postProcMesh.write_file(
1091 dm,
"ELASTIC", fe_elastic_rhs_ptr, 0, n_parts);
1093 if (delete_prisms) {
1094 Range faces, tris, quads, tris_edges, quads_edges, ents_to_delete;
1096 CHKERR moab.get_adjacencies(contact_prisms, 2,
true, faces,
1097 moab::Interface::UNION);
1098 tris = faces.subset_by_type(MBTRI);
1099 quads = faces.subset_by_type(MBQUAD);
1100 CHKERR moab.get_adjacencies(tris, 1,
true, tris_edges,
1101 moab::Interface::UNION);
1102 CHKERR moab.get_adjacencies(quads, 1,
true, quads_edges,
1103 moab::Interface::UNION);
1105 ents_to_delete.merge(contact_prisms);
1106 ents_to_delete.merge(quads);
1107 ents_to_delete.merge(subtract(quads_edges, tris_edges));
1109 CHKERR moab.delete_entities(ents_to_delete);
1112 PetscPrintf(PETSC_COMM_WORLD,
"Write file %s\n", output_mesh_name);
1113 CHKERR moab.write_file(output_mesh_name,
"MOAB");
1116 auto get_tag_handle = [&](
auto name,
auto size) {
1118 std::vector<double> def_vals(size, 0.0);
1119 CHKERR moab.tag_get_handle(name, size, MB_TYPE_DOUBLE,
th,
1120 MB_TAG_CREAT | MB_TAG_SPARSE,
1127 CHKERR moab.get_entities_by_dimension(0, 3, tets);
1129 std::array<double, 9> internal_stress, actual_stress;
1130 std::array<double, 9> internal_stress_ref, actual_stress_ref;
1131 std::array<double, 2> nb_gauss_pts_ref, contact_area_ref;
1134 internal_stress_ref = {5., 5., 5., 0., 0., 0., 0., 0., 0.};
1135 actual_stress_ref = {0., 0., 1., 0., 0., 0., 0., 0., 0.};
1138 internal_stress_ref = {5., 5., 5., 0., 0., 0., 0., 0., 0.};
1139 actual_stress_ref = {0., 5. / 3., 5. / 3., 0., 0., 0., 0., 0., 0.};
1142 actual_stress_ref = {0., 0., -100., 0., 0., 0., 0., 0., 0.};
1143 if (strcmp(stress_tag_name,
"INTERNAL_STRESS") == 0)
1144 internal_stress_ref = {0., 0., -200., 0., 0., 0., 0., 0., 0.};
1146 internal_stress_ref = {0., 0., -100., 0., 0., 0., 0., 0., 0.};
1149 nb_gauss_pts_ref = {96, 192};
1150 contact_area_ref = {0.125, 0.25};
1153 SETERRQ1(PETSC_COMM_SELF,
MOFEM_NOT_FOUND,
"Test number %d not found",
1157 auto th_internal_stress = get_tag_handle(
"MED_INTERNAL_STRESS", 9);
1158 auto th_actual_stress = get_tag_handle(
"MED_ACTUAL_STRESS", 9);
1159 CHKERR moab.tag_get_data(th_internal_stress, &tet, 1,
1160 internal_stress.data());
1161 CHKERR moab.tag_get_data(th_actual_stress, &tet, 1,
1162 actual_stress.data());
1163 const double eps = 1e-10;
1164 if (test_num == 4) {
1165 if (std::abs(nb_gauss_pts_ref[0] - nb_gauss_pts[0]) >
eps) {
1167 "Wrong number of active contact gauss pts: should be %d "
1169 (
int)nb_gauss_pts_ref[0], (
int)nb_gauss_pts[0]);
1171 if (std::abs(nb_gauss_pts_ref[1] - nb_gauss_pts[1]) >
eps) {
1173 "Wrong total number of contact gauss pts: should be %d "
1175 (
int)nb_gauss_pts_ref[1], (
int)nb_gauss_pts[1]);
1177 if (std::abs(contact_area_ref[0] - contact_area[0]) >
eps) {
1179 "Wrong active contact area: should be %g "
1181 contact_area_ref[0], contact_area[0]);
1183 if (std::abs(contact_area_ref[1] - contact_area[1]) >
eps) {
1185 "Wrong potential contact area: should be %g "
1187 contact_area_ref[1], contact_area[1]);
1190 if (save_mean_stress) {
1191 for (
int i = 0;
i < 9;
i++) {
1192 if (std::abs(internal_stress[
i] - internal_stress_ref[
i]) >
eps) {
1194 "Wrong component %d of internal stress: should be %g "
1196 i, internal_stress_ref[
i], internal_stress[
i]);
1198 if (std::abs(actual_stress[
i] - actual_stress_ref[
i]) >
eps) {
1200 "Wrong component %d of actual stress: should be %g "
1202 i, actual_stress_ref[
i], actual_stress[
i]);