35int main(
int argc,
char *argv[]) {
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"
52 "-snes_converged_reason \n"
61 "-eigen_pos_flag 0 \n";
64 if (!
static_cast<bool>(ifstream(
param_file))) {
65 std::ofstream file(
param_file.c_str(), std::ios::ate);
76 moab::Core mb_instance;
77 moab::Interface &moab = mb_instance;
81 PetscBool flg_file_out;
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;
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;
110 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Elastic Config",
"none");
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);
117 CHKERR PetscOptionsInt(
"-order",
"approximation order of spatial positions",
118 "", 1, &
order, PETSC_NULL);
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);
129 CHKERR PetscOptionsInt(
"-step_num",
"number of steps",
"", nb_sub_steps,
130 &nb_sub_steps, PETSC_NULL);
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);
137 CHKERR PetscOptionsReal(
"-cn_value",
"default regularisation cn value",
"",
138 1., &cn_value, PETSC_NULL);
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,
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);
158 CHKERR PetscOptionsReal(
"-scale_factor",
"scale factor",
"", scale_factor,
159 &scale_factor, PETSC_NULL);
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,
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);
180 ierr = PetscOptionsEnd();
184 if (flg_file != PETSC_TRUE) {
185 SETERRQ(PETSC_COMM_SELF, 1,
"*** ERROR -file_name (MESH FILE NEEDED)");
188 if (is_partitioned) {
190 "Partitioned mesh is not supported");
204 std::vector<BitRefLevel> bit_levels;
207 CHKERR bit_ref_manager->setBitRefLevelByDim(0, 3, bit_levels.back());
209 auto add_prism_interface = [&](std::vector<BitRefLevel> &bit_levels) {
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());
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(
229 CHKERR prism_interface->getSides(cubit_meshset, bit_levels.back(),
232 bit_levels.push_back(
BitRefLevel().set(bit_levels.size()));
234 CHKERR prism_interface->splitSides(ref_level_meshset,
235 bit_levels.back(), cubit_meshset,
238 CHKERR moab.delete_entities(&ref_level_meshset, 1);
243 CHKERR bit_ref_manager->updateMeshsetByEntitiesChildren(
244 cubit_meshset, bit_levels.back(), cubit_meshset, MBVERTEX,
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);
254 CHKERR bit_ref_manager->shiftRightBitRef(1);
255 bit_levels.pop_back();
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) {
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);
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);
287 Range simple_contact_prisms, simple_master_tris, simple_slave_tris;
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);
295 Range mortar_contact_prisms, mortar_master_tris, mortar_slave_tris;
298 if (it->getName().compare(0, 13,
"MORTAR_MASTER") == 0) {
300 it->meshset, MBTRI, mortar_master_tris,
true);
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);
313 boost::shared_ptr<ContactSearchKdTree::ContactCommonData_multiIndex>
314 contact_commondata_multi_index;
316 contact_commondata_multi_index =
317 boost::shared_ptr<ContactSearchKdTree::ContactCommonData_multiIndex>(
323 mortar_master_tris, mortar_slave_tris, contact_commondata_multi_index,
324 mortar_contact_prisms);
327 CHKERR moab.create_meshset(MESHSET_SET, meshset_slave_master_prisms);
330 moab.add_entities(meshset_slave_master_prisms, mortar_contact_prisms);
333 meshset_slave_master_prisms, 3, bit_levels.back());
335 Range tris_from_prism;
337 CHKERR m_field.
get_moab().get_adjacencies(mortar_contact_prisms, 2,
true,
339 moab::Interface::UNION);
341 tris_from_prism = tris_from_prism.subset_by_type(MBTRI);
342 mortar_slave_tris = intersect(tris_from_prism, mortar_slave_tris);
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);
354 if (!eigen_pos_flag) {
375 if (!eigen_pos_flag) {
389 auto set_contact_order = [&](
Range &contact_prisms,
int order_contact,
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);
411 ho_ents = unite(ho_ents, ents);
412 ho_ents = unite(ho_ents, tets);
421 if (!ignore_contact && order_contact >
order) {
422 CHKERR set_contact_order(all_contact_prisms, order_contact, nb_ho_levels);
440 CHKERR moab.get_connectivity(all_slave_tris, slave_verts,
true);
442 slave_verts,
"LAGMULT");
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);
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();
456 CHKERR HookeElement::addElasticElement(m_field, block_sets_ptr,
"ELASTIC",
458 "MESH_NODE_POSITIONS",
false);
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](
476 t_thermal_strain(
i,
j) = 0.0;
481 temp = init_temp - 1.0;
482 t_thermal_strain(
i,
j) =
483 thermal_expansion_coef * (
temp - init_temp) *
t_kd(
i,
j);
487 t_thermal_strain(
i,
j) = -thermal_expansion_coef *
t_kd(
i,
j);
490 t_thermal_strain(2, 2) = thermal_expansion_coef;
493 t_thermal_strain(
i,
j) = thermal_expansion_coef *
t_kd(
i,
j);
499 return t_thermal_strain;
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,
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));
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,
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);
540 Skinner skinner(&moab);
542 CHKERR skinner.find_skin(0, all_tets,
false, skin_tris);
553 "MESH_NODE_POSITIONS");
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");
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");
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");
577 auto make_mortar_contact_common_data = [&]() {
578 return boost::make_shared<MortarContactProblem::CommonDataMortarContact>(
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;
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;
596 auto get_mortar_master_traction_rhs = [&](
auto contact_problem,
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;
607 auto get_mortar_master_traction_lhs = [&](
auto contact_problem,
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;
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;
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;
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;
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);
655 mortar_contact_problem->addMortarContactElement(
656 "MORTAR_CONTACT_ELEM",
"SPATIAL_POSITION",
"LAGMULT",
657 mortar_contact_prisms);
659 mortar_contact_problem->addMortarContactElement(
660 "MORTAR_CONTACT_ELEM",
"SPATIAL_POSITION",
"LAGMULT",
661 mortar_contact_prisms,
"MESH_NODE_POSITIONS", eigen_pos_flag,
664 auto make_simple_contact_element = [&]() {
665 return boost::make_shared<SimpleContactProblem::SimpleContactElement>(
669 auto make_convective_simple_master_element = [&]() {
670 return boost::make_shared<
672 m_field,
"SPATIAL_POSITION",
"MESH_NODE_POSITIONS");
675 auto make_convective_simple_slave_element = [&]() {
676 return boost::make_shared<
678 m_field,
"SPATIAL_POSITION",
"MESH_NODE_POSITIONS");
681 auto make_simple_contact_common_data = [&]() {
682 return boost::make_shared<SimpleContactProblem::CommonDataSimpleContact>(
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;
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;
700 auto get_simple_master_traction_rhs = [&](
auto contact_problem,
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;
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;
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;
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);
737 simple_contact_problem->addPostProcContactElement(
738 "CONTACT_POST_PROC",
"SPATIAL_POSITION",
"LAGMULT",
739 "MESH_NODE_POSITIONS", all_slave_tris);
745 "MESH_NODE_POSITIONS");
760 DMType dm_name =
"DMMOFEM";
767 CHKERR DMSetType(dm, dm_name);
771 CHKERR DMSetFromOptions(dm);
793 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
794 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
797 CHKERR VecGhostUpdateBegin(
F, INSERT_VALUES, SCATTER_FORWARD);
798 CHKERR VecGhostUpdateEnd(
F, INSERT_VALUES, SCATTER_FORWARD);
800 CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
801 CHKERR MatZeroEntries(Aij);
803 fe_elastic_rhs_ptr->snes_f =
F;
804 fe_elastic_lhs_ptr->snes_B = Aij;
807 boost::shared_ptr<FEMethod> dirichlet_bc_ptr =
809 m_field,
"SPATIAL_POSITION", Aij,
D,
F));
812 dirichlet_bc_ptr->snes_x =
D;
815 boost::ptr_map<std::string, NeumannForcesSurface> neumann_forces;
817 m_field, neumann_forces, NULL,
"SPATIAL_POSITION");
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);
829 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr(
831 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr(
835 m_field, fe_spring_lhs_ptr, fe_spring_rhs_ptr,
"SPATIAL_POSITION",
836 "MESH_NODE_POSITIONS");
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);
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);
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);
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);
898 dirichlet_bc_ptr.get());
900 boost::shared_ptr<FEMethod> 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);
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);
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);
940 get_contact_lhs(simple_contact_problem,
941 make_simple_contact_element,
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);
958 char testing_options[] =
"-ksp_type fgmres "
960 "-pc_factor_mat_solver_type mumps "
961 "-snes_type newtonls "
962 "-snes_linesearch_type basic "
966 CHKERR PetscOptionsInsertString(PETSC_NULL, testing_options);
970 CHKERR SNESSetDM(snes, dm);
974 CHKERR SNESSetDM(snes, dm);
978 CHKERR SNESSetFromOptions(snes);
989 struct OpGetFieldGradientValuesOnSkin
990 :
public FaceElementForcesAndSourcesCore::UserDataOperator {
992 const std::string feVolName;
993 boost::shared_ptr<VolSideFe> sideOpFe;
995 OpGetFieldGradientValuesOnSkin(
const std::string
field_name,
996 const std::string vol_fe_name,
997 boost::shared_ptr<VolSideFe> side_fe)
998 : FaceElementForcesAndSourcesCore::UserDataOperator(
1000 feVolName(vol_fe_name), sideOpFe(side_fe) {}
1005 if (type != MBVERTEX)
1007 CHKERR loopSideVolumes(feVolName, *sideOpFe);
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));
1022 new OpGetFieldGradientValuesOnSkin(
"SPATIAL_POSITION",
"ELASTIC",
1026 "SPATIAL_POSITION", data_hooke_element_at_pts->spatPosMat));
1029 "MESH_NODE_POSITIONS", data_hooke_element_at_pts->meshNodePosMat));
1034 "SPATIAL_POSITION", data_hooke_element_at_pts,
1038 for (
int ss = 0; ss != nb_sub_steps; ++ss) {
1039 if (!ignore_pressure) {
1043 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Ignoring pressure...\n");
1045 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Load scale: %6.4e\n",
1048 CHKERR SNESSolve(snes, PETSC_NULL,
D);
1050 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
1051 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
1058 CHKERR HookeElement::calculateEnergy(dm, block_sets_ptr,
"SPATIAL_POSITION",
1059 "MESH_NODE_POSITIONS",
false,
false,
1062 CHKERR VecSum(v_energy, &energy);
1064 MOFEM_LOG_C(
"WORLD", Sev::inform,
"Elastic energy: %8.8e", energy);
1068 moab::Interface &moab_proc = mb_post;
1070 auto common_data_mortar_contact = make_mortar_contact_common_data();
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();
1077 fe_post_proc_mortar_contact = make_mortar_contact_element();
1080 std::array<double, 2> nb_gauss_pts;
1081 std::array<double, 2> contact_area;
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})
1093 CHKERR VecRestoreArrayRead(vec, &array);
1097 mb_post.delete_mesh();
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,
1105 CHKERR VecZeroEntries(common_data_mortar_contact->gaussPtsStateVec);
1106 CHKERR VecZeroEntries(common_data_mortar_contact->contactAreaVec);
1109 fe_post_proc_mortar_contact);
1111 CHKERR get_contact_data(common_data_mortar_contact->gaussPtsStateVec,
1113 CHKERR get_contact_data(common_data_mortar_contact->contactAreaVec,
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]);
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.);
1129 auto common_data_simple_contact = make_simple_contact_common_data();
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();
1136 fe_post_proc_simple_contact = make_simple_contact_element();
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,
1145 CHKERR VecZeroEntries(common_data_simple_contact->gaussPtsStateVec);
1146 CHKERR VecZeroEntries(common_data_simple_contact->contactAreaVec);
1149 fe_post_proc_simple_contact);
1151 CHKERR get_contact_data(common_data_simple_contact->gaussPtsStateVec,
1153 CHKERR get_contact_data(common_data_simple_contact->contactAreaVec,
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.);
1168 if (!ignore_contact) {
1170 std::ostringstream strm;
1171 strm <<
"out_contact_integ_pts"
1174 MOFEM_LOG_C(
"WORLD", Sev::inform,
"Write file %s\n",
1177 "PARALLEL=WRITE_PART");
1180 boost::shared_ptr<PostProcFaceOnRefinedMesh> post_proc_contact_ptr(
1183 CHKERR post_proc_contact_ptr->generateReferenceElementMesh();
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");
1190 if (!all_slave_tris.empty()) {
1192 post_proc_contact_ptr);
1194 std::ostringstream stm;
1195 stm <<
"out_contact_surface"
1198 MOFEM_LOG_C(
"WORLD", Sev::inform,
"Write file %s",
1200 CHKERR post_proc_contact_ptr->postProcMesh.write_file(
1208 PetscPrintf(PETSC_COMM_WORLD,
"Loop post proc on the skin\n");
1212 stm <<
"out_skin.h5m";
1214 MOFEM_LOG_C(
"WORLD", Sev::inform,
"Write file %s",
1222 dm,
"ELASTIC", fe_elastic_rhs_ptr, 0, n_parts);
1224 Range faces, tris, quads, tris_edges, quads_edges, ents_to_delete;
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);
1235 ents_to_delete.merge(all_contact_prisms);
1236 ents_to_delete.merge(quads);
1237 ents_to_delete.merge(subtract(quads_edges, tris_edges));
1239 CHKERR moab.delete_entities(ents_to_delete);
1241 MOFEM_LOG_C(
"WORLD", Sev::inform,
"Write file %s", output_mesh_name);
1242 CHKERR moab.write_file(output_mesh_name,
"MOAB");
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,
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;
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.};
1266 internal_stress_ref = {0., 0., -100., 0., 0., 0., 0., 0., 0.};
1269 nb_gauss_pts_ref = {216, 492};
1270 contact_area_ref = {0.125, 0.25};
1273 SETERRQ1(PETSC_COMM_SELF,
MOFEM_NOT_FOUND,
"Test number %d not found",
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 "
1289 (
int)nb_gauss_pts_ref[0], (
int)nb_gauss_pts[0]);
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 "
1295 (
int)nb_gauss_pts_ref[1], (
int)nb_gauss_pts[1]);
1297 if (std::abs(contact_area_ref[0] - contact_area[0]) >
eps) {
1299 "Wrong active contact area: should be %g "
1301 contact_area_ref[0], contact_area[0]);
1303 if (std::abs(contact_area_ref[1] - contact_area[1]) >
eps) {
1305 "Wrong potential contact area: should be %g "
1307 contact_area_ref[1], contact_area[1]);
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 "
1317 i, internal_stress_ref[
i], internal_stress[
i]);
1319 if (std::abs(actual_stress[
i] - actual_stress_ref[
i]) >
eps) {
1321 "Wrong component %d of actual stress: should be %g "
1323 i, actual_stress_ref[
i], actual_stress[
i]);
const std::string default_options
#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.
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.
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)
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