39int main(
int argc,
char *argv[]) {
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";
68 if (!
static_cast<bool>(ifstream(
param_file))) {
69 std::ofstream file(
param_file.c_str(), std::ios::ate);
80 moab::Core mb_instance;
81 moab::Interface &moab = mb_instance;
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;
105 PetscReal init_temp = 250.0;
106 PetscReal scale_factor = 1.0;
107 PetscBool analytical_input = PETSC_TRUE;
108 char stress_tag_name[255];
109 PetscBool flg_tag_name;
110 PetscBool save_mean_stress = PETSC_TRUE;
111 PetscBool ignore_contact = PETSC_FALSE;
112 PetscBool ignore_pressure = PETSC_FALSE;
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);
195 CHKERR PetscOptionsReal(
"-my_init_temp",
"init_temp ",
"", init_temp,
196 &init_temp, 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>>();
516 CHKERR HookeElement::setBlocks(m_field, block_sets_ptr);
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();
525 CHKERR HookeElement::addElasticElement(m_field, block_sets_ptr,
"ELASTIC",
527 "MESH_NODE_POSITIONS",
false);
529 auto data_hooke_element_at_pts =
530 boost::make_shared<HookeInternalStressElement::DataAtIntegrationPts>();
531 CHKERR HookeElement::setOperators(fe_elastic_lhs_ptr, fe_elastic_rhs_ptr,
532 block_sets_ptr,
"SPATIAL_POSITION",
533 "MESH_NODE_POSITIONS",
false,
false,
534 MBTET, data_hooke_element_at_pts);
535 auto thermal_strain =
536 [&thermal_expansion_coef, &init_temp, &test_num](
544 t_thermal_strain(
i,
j) = 0.0;
549 temp = init_temp - 1.0;
550 t_thermal_strain(
i,
j) =
551 thermal_expansion_coef * (
temp - init_temp) *
t_kd(
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));
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);
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)
920 CHKERR loopSideVolumes(feVolName, *sideOpFe);
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));
935 new OpGetFieldGradientValuesOnSkin(
"SPATIAL_POSITION",
"ELASTIC",
936 my_vol_side_fe_ptr));
939 "SPATIAL_POSITION", data_hooke_element_at_pts->spatPosMat));
942 "MESH_NODE_POSITIONS", data_hooke_element_at_pts->meshNodePosMat));
947 "SPATIAL_POSITION", data_hooke_element_at_pts,
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);
971 CHKERR HookeElement::calculateEnergy(dm, block_sets_ptr,
"SPATIAL_POSITION",
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());
992 moab::Interface &moab_proc = mb_post;
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]);
const std::string default_options
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
@ MOFEM_DATA_INCONSISTENCY
#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.
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
MoFEMErrorCode writeFile(const std::string file_name, const char *file_type="MOAB", const char *file_options="PARALLEL=WRITE_PART")
wrote results in (MOAB) format, use "file_name.h5m"
#define MOFEM_LOG(channel, severity)
Log.
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets in a moFEM field.
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
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.
implementation of Data Operators for Forces and Sources
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
static constexpr double delta
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)
default operator for TRI element
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Create interface from given surface and insert flat prisms in-between.
Projection of edge entities with one mid-node on hierarchical basis.
intrusive_ptr for managing petsc objects
Interface for nonlinear (SNES) solver.
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.
Base volume element used to integrate on skeleton.
data for calculation heat conductivity and heat capacity elements
MoFEMErrorCode generateReferenceElementMesh()
std::vector< EntityHandle > mapGaussPts
moab::Interface & postProcMesh
int operator()(int, int, int order) const