111 "-pc_factor_mat_solver_type mumps \n"
112 "-mat_mumps_icntl_20 0 \n"
114 "-snes_type newtonls \n"
115 "-snes_linesearch_type basic \n"
120 "-ts_type beuler \n";
123 if (!
static_cast<bool>(ifstream(
param_file))) {
124 std::ofstream file(
param_file.c_str(), std::ios::ate);
125 if (file.is_open()) {
132 auto core_log = logging::core::get();
134 LogManager::createSink(LogManager::getStrmWorld(),
"ELASTIC"));
135 LogManager::setLog(
"ELASTIC");
140 PetscBool flg_block_config, flg_file;
142 char block_config_file[255];
143 PetscInt test_nb = 0;
145 PetscBool is_partitioned = PETSC_FALSE;
146 PetscBool is_calculating_frequency = PETSC_FALSE;
147 PetscBool is_post_proc_volume = PETSC_TRUE;
150 enum bases { LEGENDRE, LOBATTO, BERNSTEIN_BEZIER, JACOBI, LASBASETOP };
151 const char *list_bases[] = {
"legendre",
"lobatto",
"bernstein_bezier",
153 PetscInt choice_base_value = LOBATTO;
156 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Elastic Config",
"none");
158 CHKERR PetscOptionsString(
"-my_file",
"mesh file name",
"",
"mesh.h5m",
161 CHKERR PetscOptionsInt(
"-my_order",
"default approximation order",
"",
164 CHKERR PetscOptionsEList(
"-base",
"approximation base",
"", list_bases,
165 LASBASETOP, list_bases[choice_base_value],
166 &choice_base_value, PETSC_NULL);
168 CHKERR PetscOptionsInt(
"-is_atom_test",
"ctest number",
"", test_nb,
169 &test_nb, PETSC_NULL);
171 CHKERR PetscOptionsBool(
"-my_is_partitioned",
172 "set if mesh is partitioned (this result that each "
173 "process keeps only one part of the mesh)",
174 "", is_partitioned, &is_partitioned, PETSC_NULL);
176 CHKERR PetscOptionsString(
"-my_block_config",
"elastic configure file name",
177 "",
"block_conf.in", block_config_file, 255,
181 "-my_is_calculating_frequency",
"set if frequency will be calculated",
182 "", is_calculating_frequency, &is_calculating_frequency, PETSC_NULL);
184 CHKERR PetscOptionsBool(
"-my_is_post_proc_volume",
185 "if true post proc volume",
"", is_post_proc_volume,
186 &is_post_proc_volume, PETSC_NULL);
188 ierr = PetscOptionsEnd();
192 if (flg_file != PETSC_TRUE) {
193 SETERRQ(PETSC_COMM_SELF, 1,
"*** ERROR -my_file (MESH FILE NEEDED)");
204 MPI_Comm moab_comm_world;
205 MPI_Comm_dup(PETSC_COMM_WORLD, &moab_comm_world);
206 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
208 pcomm =
new ParallelComm(&moab, moab_comm_world);
211 if (is_partitioned == PETSC_TRUE) {
215 option =
"PARALLEL=READ_PART;"
216 "PARALLEL_RESOLVE_SHARED_ENTS;"
217 "PARTITION=PARALLEL_PARTITION;";
240 bool mesh_has_tets =
false;
241 bool mesh_has_prisms =
false;
246 CHKERR moab.get_number_entities_by_type(0, MBTET, nb_tets,
true);
247 CHKERR moab.get_number_entities_by_type(0, MBHEX, nb_hexs,
true);
248 CHKERR moab.get_number_entities_by_type(0, MBPRISM, nb_prisms,
true);
250 mesh_has_tets = (nb_tets + nb_hexs) > 0;
251 mesh_has_prisms = nb_prisms > 0;
265 if (
bit->getName().compare(0, 3,
"ROD") == 0) {
273 switch (choice_base_value) {
280 case BERNSTEIN_BEZIER:
305 CHKERR m_field.
get_moab().get_entities_by_type(0, MBEDGE, all_edges,
true);
308 Range edges_in_simple_rod;
310 if (
bit->getName().compare(0, 3,
"ROD") == 0) {
313 MBEDGE, edges,
true);
314 edges_in_simple_rod.merge(edges);
325 Range edges_to_set_order;
326 edges_to_set_order = subtract(all_edges, edges_in_simple_rod);
346 auto setting_second_order_geometry = [&m_field]() {
352 moab::Interface::UNION);
369 CHKERR setting_second_order_geometry();
373 std::map<int, BlockOptionData> block_data;
374 auto setting_blocks_data_and_order_from_config_file =
375 [&m_field, &moab, &block_data, flg_block_config, block_config_file,
376 order](boost::shared_ptr<std::map<int, BlockData>> &block_sets_ptr) {
378 if (flg_block_config) {
379 ifstream ini_file(block_config_file);
380 po::variables_map vm;
381 po::options_description config_file_options;
384 std::ostringstream str_order;
385 str_order <<
"block_" << it->getMeshsetId()
386 <<
".displacement_order";
387 config_file_options.add_options()(
388 str_order.str().c_str(),
389 po::value<int>(&block_data[it->getMeshsetId()].oRder)
390 ->default_value(
order));
392 std::ostringstream str_cond;
393 str_cond <<
"block_" << it->getMeshsetId() <<
".young_modulus";
394 config_file_options.add_options()(
395 str_cond.str().c_str(),
396 po::value<double>(&block_data[it->getMeshsetId()].yOung)
397 ->default_value(-1));
399 std::ostringstream str_capa;
400 str_capa <<
"block_" << it->getMeshsetId() <<
".poisson_ratio";
401 config_file_options.add_options()(
402 str_capa.str().c_str(),
403 po::value<double>(&block_data[it->getMeshsetId()].pOisson)
404 ->default_value(-2));
406 std::ostringstream str_init_temp;
407 str_init_temp <<
"block_" << it->getMeshsetId()
408 <<
".initial_temperature";
409 config_file_options.add_options()(
410 str_init_temp.str().c_str(),
411 po::value<double>(&block_data[it->getMeshsetId()].initTemp)
414 po::parsed_options parsed =
415 parse_config_file(ini_file, config_file_options,
true);
420 if (block_data[it->getMeshsetId()].oRder == -1)
422 if (block_data[it->getMeshsetId()].oRder ==
order)
424 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Set block %d order to %d",
426 block_data[it->getMeshsetId()].oRder);
428 CHKERR moab.get_entities_by_handle(it->getMeshset(), block_ents,
430 Range ents_to_set_order;
431 CHKERR moab.get_adjacencies(block_ents, 3,
false,
433 moab::Interface::UNION);
434 ents_to_set_order = ents_to_set_order.subset_by_dimension(3);
435 CHKERR moab.get_adjacencies(block_ents, 2,
false,
437 moab::Interface::UNION);
438 CHKERR moab.get_adjacencies(block_ents, 1,
false,
440 moab::Interface::UNION);
445 ents_to_set_order,
"DISPLACEMENT",
446 block_data[it->getMeshsetId()].oRder);
448 std::vector<std::string> additional_parameters;
449 additional_parameters =
450 collect_unrecognized(parsed.options, po::include_positional);
451 for (std::vector<std::string>::iterator vit =
452 additional_parameters.begin();
453 vit != additional_parameters.end(); vit++) {
454 MOFEM_LOG_C(
"ELASTIC", Sev::warning,
"Unrecognized option %s",
463 const int id = it->getMeshsetId();
464 auto &bd = (*block_sets_ptr)[id];
465 if (block_data[
id].yOung > 0)
466 bd.E = block_data[id].yOung;
467 if (block_data[
id].pOisson >= -1)
468 bd.PoissonRatio = block_data[id].pOisson;
469 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Block %d",
id);
470 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"\tYoung modulus %3.4g", bd.E);
471 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"\tPoisson ratio %3.4g",
480 boost::shared_ptr<std::map<int, HookeElement::BlockData>> block_sets_ptr =
481 boost::make_shared<std::map<int, HookeElement::BlockData>>();
483 CHKERR setting_blocks_data_and_order_from_config_file(block_sets_ptr);
485 boost::shared_ptr<std::map<int, MassBlockData>> mass_block_sets_ptr =
486 boost::make_shared<std::map<int, MassBlockData>>();
490 boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
492 boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
493 fe_lhs_ptr->getRuleHook =
VolRule();
494 fe_rhs_ptr->getRuleHook =
VolRule();
501 boost::shared_ptr<ForcesAndSourcesCore> prism_fe_lhs_ptr(
503 boost::shared_ptr<ForcesAndSourcesCore> prism_fe_rhs_ptr(
508 "MESH_NODE_POSITIONS",
false);
510 auto add_skin_element_for_post_processing = [&]() {
512 Range elastic_element_ents;
514 "ELASTIC", 3, elastic_element_ents);
517 CHKERR skin.find_skin(0, elastic_element_ents,
false, skin_faces);
519 if (is_partitioned) {
520 CHKERR pcomm->filter_pstatus(skin_faces,
521 PSTATUS_SHARED | PSTATUS_MULTISHARED,
522 PSTATUS_NOT, -1, &proc_skin);
524 proc_skin = skin_faces;
534 "POST_PROC_SKIN",
"MESH_NODE_POSITIONS");
539 CHKERR add_skin_element_for_post_processing();
541 auto data_at_pts = boost::make_shared<HookeElement::DataAtIntegrationPts>();
544 "DISPLACEMENT",
"MESH_NODE_POSITIONS",
545 false,
true, MBTET, data_at_pts);
547 if (mesh_has_prisms) {
549 prism_fe_lhs_ptr, prism_fe_rhs_ptr, block_sets_ptr,
"DISPLACEMENT",
550 "MESH_NODE_POSITIONS",
false,
true, MBPRISM, data_at_pts);
555 auto thermal_strain =
558 constexpr
double alpha = 1;
563 return t_thermal_strain;
566 fe_rhs_ptr->getOpPtrVector().push_back(
567 new HookeElement::OpAnalyticalInternalStrain_dx<0>(
568 "DISPLACEMENT", data_at_pts, thermal_strain));
571 boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_mass_ptr(
574 for (
auto &sit : *block_sets_ptr) {
575 for (
auto &mit : *mass_block_sets_ptr) {
576 fe_mass_ptr->getOpPtrVector().push_back(
577 new HookeElement::OpCalculateMassMatrix(
"DISPLACEMENT",
578 "DISPLACEMENT", sit.second,
579 mit.second, data_at_pts));
586 "MESH_NODE_POSITIONS");
590 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr(
592 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr(
596 fe_spring_rhs_ptr,
"DISPLACEMENT",
597 "MESH_NODE_POSITIONS");
602 "MESH_NODE_POSITIONS");
609 boost::shared_ptr<EdgeEle> fe_simple_rod_lhs_ptr(
new EdgeEle(m_field));
610 boost::shared_ptr<EdgeEle> fe_simple_rod_rhs_ptr(
new EdgeEle(m_field));
614 m_field, fe_simple_rod_lhs_ptr, fe_simple_rod_rhs_ptr,
"DISPLACEMENT",
615 "MESH_NODE_POSITIONS");
627 "MESH_NODE_POSITIONS");
632 CHKERR m_field.
get_moab().get_entities_by_dimension(it->meshset, 3, tets,
649 fluid_pressure_fe.addNeumannFluidPressureBCElements(
"DISPLACEMENT");
655 bool add_temp_field =
false;
657 if (block_data[it->getMeshsetId()].initTemp != 0) {
658 add_temp_field =
true;
662 if (add_temp_field) {
671 CHKERR thermal_stress_elem.addThermalStressElement(
672 "ELASTIC",
"DISPLACEMENT",
"TEMP");
681 "MESH_NODE_POSITIONS");
685 if (block_data[it->getMeshsetId()].initTemp != 0) {
687 "Set block %d temperature to %3.2g\n", it->getMeshsetId(),
688 block_data[it->getMeshsetId()].initTemp);
690 CHKERR moab.get_entities_by_handle(it->meshset, block_ents,
true);
692 CHKERR moab.get_connectivity(block_ents, vertices,
true);
694 block_data[it->getMeshsetId()].initTemp, MBVERTEX, vertices,
712 CHKERR DMSetFromOptions(dm);
733 CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
737 if (is_calculating_frequency == PETSC_TRUE) {
739 CHKERR MatSetOption(Mij, MAT_SPD, PETSC_TRUE);
744 fe_spring_lhs_ptr->ksp_B = Aij;
745 fe_spring_rhs_ptr->ksp_f = F;
748 fe_simple_rod_lhs_ptr->ksp_B = Aij;
749 fe_simple_rod_rhs_ptr->ksp_f = F;
753 CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
754 CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
756 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
757 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
759 CHKERR MatZeroEntries(Aij);
771 auto dirichlet_bc_ptr = boost::make_shared<DirichletDisplacementBc>(
772 m_field,
"DISPLACEMENT", Aij, D0, F);
776 dirichlet_bc_ptr->snes_ctx = FEMethod::CTX_SNESNONE;
777 dirichlet_bc_ptr->ts_ctx = FEMethod::CTX_TSNONE;
780 CHKERR VecZeroEntries(D0);
781 CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
782 CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
790 CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
791 CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
800 fe_rhs_ptr->snes_f = F;
801 prism_fe_rhs_ptr->snes_f = F;
802 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Assemble external force vector ...";
805 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
807 fe_lhs_ptr->snes_B = Aij;
808 prism_fe_lhs_ptr->snes_B = Aij;
809 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Calculate stiffness matrix ...";
812 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
822 if (is_calculating_frequency == PETSC_TRUE) {
824 fe_mass_ptr->snes_B = Mij;
825 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Calculate mass matrix ...";
827 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
835 boost::ptr_map<std::string, NeumannForcesSurface> neumann_forces;
840 boost::ptr_map<std::string, NeumannForcesSurface>::iterator mit =
841 neumann_forces.begin();
842 for (; mit != neumann_forces.end(); mit++) {
844 &mit->second->getLoopFe());
848 boost::ptr_map<std::string, NodalForce> nodal_forces;
853 boost::ptr_map<std::string, NodalForce>::iterator fit =
854 nodal_forces.begin();
855 for (; fit != nodal_forces.end(); fit++) {
857 &fit->second->getLoopFe());
861 boost::ptr_map<std::string, EdgeForce> edge_forces;
865 auto fit = edge_forces.begin();
866 for (; fit != edge_forces.end(); fit++) {
867 auto &fe = fit->second->getLoopFe();
875 CHKERR body_forces_methods.addBlock(
"DISPLACEMENT", F,
879 &body_forces_methods.getLoopFe());
883 CHKERR fluid_pressure_fe.setNeumannFluidPressureFiniteElementOperators(
884 "DISPLACEMENT", F,
false,
true);
887 &fluid_pressure_fe.getLoopFe());
892 PetscViewerPushFormat(
893 PETSC_VIEWER_STDOUT_SELF,
894 PETSC_VIEWER_ASCII_MATLAB);
901 if (is_calculating_frequency == PETSC_TRUE) {
902 CHKERR MatAssemblyBegin(Mij, MAT_FINAL_ASSEMBLY);
903 CHKERR MatAssemblyEnd(Mij, MAT_FINAL_ASSEMBLY);
909 CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
910 CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
911 CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
912 CHKERR VecAssemblyBegin(F);
917 auto solver =
createKSP(PETSC_COMM_WORLD);
918 CHKERR KSPSetDM(solver, dm);
919 CHKERR KSPSetFromOptions(solver);
920 CHKERR KSPSetOperators(solver, Aij, Aij);
925 PetscBool same = PETSC_FALSE;
927 CHKERR KSPGetPC(solver, &pc);
928 PetscObjectTypeCompare((PetscObject)pc, PCMG, &same);
932 CHKERR PCSetFromOptions(pc);
935 CHKERR KSPSetDMActive(solver, PETSC_FALSE);
938 CHKERR KSPSetInitialGuessKnoll(solver, PETSC_FALSE);
939 CHKERR KSPSetInitialGuessNonzero(solver, PETSC_TRUE);
943 auto set_post_proc_skin = [&](
auto &post_proc_skin) {
947 CHKERR post_proc_skin.generateReferenceElementMesh();
948 CHKERR post_proc_skin.addFieldValuesPostProc(
"DISPLACEMENT");
949 CHKERR post_proc_skin.addFieldValuesPostProc(
"MESH_NODE_POSITIONS");
950 CHKERR post_proc_skin.addFieldValuesGradientPostProcOnSkin(
951 "DISPLACEMENT",
"ELASTIC", data_at_pts->hMat,
true);
952 CHKERR post_proc_skin.addFieldValuesGradientPostProcOnSkin(
954 "MESH_NODE_POSITIONS",
"ELASTIC", data_at_pts->HMat,
false);
955 post_proc_skin.getOpPtrVector().push_back(
956 new HookeElement::OpPostProcHookeElement<
958 "DISPLACEMENT", data_at_pts, *block_sets_ptr,
959 post_proc_skin.postProcMesh, post_proc_skin.mapGaussPts,
true,
964 auto set_post_proc_tets = [&](
auto &post_proc) {
967 CHKERR post_proc.generateReferenceElementMesh();
970 CHKERR post_proc.addFieldValuesPostProc(
"DISPLACEMENT");
971 CHKERR post_proc.addFieldValuesPostProc(
"MESH_NODE_POSITIONS");
972 CHKERR post_proc.addFieldValuesGradientPostProc(
"DISPLACEMENT");
975 m_field, post_proc.postProcMesh, post_proc.mapGaussPts,
976 "DISPLACEMENT", post_proc.commonData, block_sets_ptr.get()));
980 auto set_post_proc_edge = [&](
auto &post_proc_edge) {
982 CHKERR post_proc_edge.generateReferenceElementMesh();
983 CHKERR post_proc_edge.addFieldValuesPostProc(
"DISPLACEMENT");
987 auto set_post_proc_prisms = [&](
auto &prism_post_proc) {
989 CHKERR prism_post_proc.generateReferenceElementMesh();
990 boost::shared_ptr<MatrixDouble> inv_jac_ptr(
new MatrixDouble);
991 prism_post_proc.getOpPtrVector().push_back(
993 prism_post_proc.getOpPtrVector().push_back(
995 CHKERR prism_post_proc.addFieldValuesPostProc(
"DISPLACEMENT");
996 CHKERR prism_post_proc.addFieldValuesPostProc(
"MESH_NODE_POSITIONS");
997 CHKERR prism_post_proc.addFieldValuesGradientPostProc(
"DISPLACEMENT");
999 m_field, prism_post_proc.postProcMesh, prism_post_proc.mapGaussPts,
1000 "DISPLACEMENT", prism_post_proc.commonData, block_sets_ptr.get()));
1009 CHKERR set_post_proc_skin(post_proc_skin);
1010 CHKERR set_post_proc_tets(post_proc);
1011 CHKERR set_post_proc_prisms(prism_post_proc);
1012 CHKERR set_post_proc_edge(post_proc_edge);
1019 CHKERR VecDuplicate(F, &F_thermal);
1023 CHKERR thermal_stress_elem.setThermalStressRhsOperators(
1024 "DISPLACEMENT",
"TEMP", F_thermal);
1042 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Process step %d",
1043 sit->get_step_number());
1046 sit->get_step_number());
1048 CHKERR VecZeroEntries(F_thermal);
1049 CHKERR VecGhostUpdateBegin(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1050 CHKERR VecGhostUpdateEnd(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1058 dm,
"ELASTIC", &thermal_stress_elem.getLoopThermalStressRhs());
1061 CHKERR VecAssemblyBegin(F_thermal);
1062 CHKERR VecAssemblyEnd(F_thermal);
1064 CHKERR VecGhostUpdateBegin(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1065 CHKERR VecGhostUpdateEnd(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1069 CHKERR VecNorm(F, NORM_2, &nrm_F);
1071 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"norm2 F = %6.4e", nrm_F);
1072 PetscReal nrm_F_thermal;
1073 CHKERR VecNorm(F_thermal, NORM_2, &nrm_F_thermal);
1074 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"norm2 F_thermal = %6.4e",
1077 CHKERR VecScale(F_thermal, -1);
1079 CHKERR VecAXPY(F_thermal, 1, F);
1082 dirichlet_bc_ptr->snes_x =
D;
1083 dirichlet_bc_ptr->snes_f = F_thermal;
1087 CHKERR KSPSolve(solver, F_thermal,
D);
1090 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
1091 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
1101 if (is_post_proc_volume == PETSC_TRUE) {
1102 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Write output file ...";
1104 std::ostringstream o1;
1105 o1 <<
"out_" << sit->step_number <<
".h5m";
1107 CHKERR post_proc.writeFile(o1.str().c_str());
1108 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done ...";
1111 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Write output file skin ...";
1114 std::ostringstream o1_skin;
1115 o1_skin <<
"out_skin" << sit->step_number <<
".h5m";
1117 CHKERR post_proc_skin.writeFile(o1_skin.str().c_str());
1118 MOFEM_LOG(
"POST_PROC_SKIN", Sev::inform) <<
"done ...";
1123 CHKERR VecZeroEntries(F_thermal);
1124 CHKERR VecGhostUpdateBegin(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1125 CHKERR VecGhostUpdateEnd(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1129 dm,
"ELASTIC", &thermal_stress_elem.getLoopThermalStressRhs());
1132 CHKERR VecAssemblyBegin(F_thermal);
1133 CHKERR VecAssemblyEnd(F_thermal);
1136 CHKERR VecGhostUpdateBegin(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1137 CHKERR VecGhostUpdateEnd(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1141 CHKERR VecNorm(F, NORM_2, &nrm_F);
1143 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"norm2 F = %6.4e", nrm_F);
1144 PetscReal nrm_F_thermal;
1145 CHKERR VecNorm(F_thermal, NORM_2, &nrm_F_thermal);
1147 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"norm2 F_thermal = %6.4e",
1151 CHKERR VecScale(F_thermal, -1);
1152 CHKERR VecAXPY(F_thermal, 1, F);
1155 dirichlet_bc_ptr->snes_x =
D;
1156 dirichlet_bc_ptr->snes_f = F_thermal;
1160 CHKERR KSPSolve(solver, F_thermal,
D);
1164 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
1165 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
1169 if (is_post_proc_volume == PETSC_TRUE) {
1170 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Write output file ...";
1174 CHKERR post_proc.writeFile(
"out.h5m");
1175 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
1178 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Write output file skin ...";
1181 CHKERR post_proc_skin.writeFile(
"out_skin.h5m");
1182 MOFEM_LOG(
"POST_PROC_SKIN", Sev::inform) <<
"done";
1186 CHKERR VecDestroy(&F_thermal);
1191 CHKERR KSPSolve(solver, F,
D);
1199 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
1200 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
1204 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Post-process start ...";
1205 if (is_post_proc_volume == PETSC_TRUE) {
1211 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
1214 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Write output file ...";
1215 if (mesh_has_tets) {
1216 if (is_post_proc_volume == PETSC_TRUE) {
1218 CHKERR post_proc.writeFile(
"out.h5m");
1221 CHKERR post_proc_skin.writeFile(
"out_skin.h5m");
1223 if (mesh_has_prisms) {
1225 CHKERR prism_post_proc.writeFile(
"prism_out.h5m");
1227 if (!edges_in_simple_rod.empty())
1229 CHKERR post_proc_edge.writeFile(
"out_edge.h5m");
1230 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
1233 if (is_calculating_frequency == PETSC_TRUE) {
1236 VecDuplicate(
D, &u1);
1239 CHKERR VecDot(u1,
D, &mode_mass);
1240 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Mode mass %6.4e\n", mode_mass);
1243 VecDuplicate(
D, &v1);
1246 double mode_stiffness;
1247 CHKERR VecDot(v1,
D, &mode_stiffness);
1248 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Mode stiffness %6.4e\n",
1252 double pi = 3.14159265359;
1253 frequency = std::sqrt(mode_stiffness / mode_mass) / (2 *
pi);
1254 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Frequency %6.4e", frequency);
1258 auto calculate_strain_energy = [&]() {
1263 "MESH_NODE_POSITIONS",
false,
true,
1268 CHKERR VecSum(v_energy, &energy);
1269 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Elastic energy %6.4e", energy);
1273 if (fabs(energy - 17.129) > 1e-3)
1275 "atom test diverged!");
1278 if (fabs(energy - 5.6475e-03) > 1e-4)
1280 "atom test diverged!");
1283 if (fabs(energy - 7.4679e-03) > 1e-4)
1285 "atom test diverged!");
1288 if (fabs(energy - 2.4992e+00) > 1e-3)
1290 "atom test diverged!");
1295 CHKERR VecMin(
D, PETSC_NULL, &min);
1296 constexpr
double expected_val = 0.10001;
1297 if (fabs(min + expected_val) > 1e-10)
1299 "atom test diverged! %3.4e != %3.4e", min, expected_val);
1302 if (fabs(energy - 4.7416e-04) > 1e-8)
1304 "atom test diverged!");
1312 CHKERR calculate_strain_energy();
1314 MPI_Comm_free(&moab_comm_world);
const std::string default_options
static MoFEMErrorCode addElasticElement(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData >> &block_sets_ptr, const std::string element_name, const std::string x_field, const std::string X_field, const bool ale)
static MoFEMErrorCode setOperators(boost::shared_ptr< ForcesAndSourcesCore > fe_lhs_ptr, boost::shared_ptr< ForcesAndSourcesCore > fe_rhs_ptr, boost::shared_ptr< map< int, BlockData >> block_sets_ptr, const std::string x_field, const std::string X_field, const bool ale, const bool field_disp, const EntityType type=MBTET, boost::shared_ptr< DataAtIntegrationPts > data_at_pts=nullptr)
static MoFEMErrorCode setBlocks(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData >> &block_sets_ptr)
static MoFEMErrorCode calculateEnergy(DM dm, boost::shared_ptr< map< int, BlockData >> block_sets_ptr, const std::string x_field, const std::string X_field, const bool ale, const bool field_disp, SmartPetscObj< Vec > &v_energy_ptr)
#define MOFEM_LOG_C(channel, severity, format,...)
MoFEMErrorCode PCMGSetUpViaApproxOrders(PC pc, PCMGSetUpViaApproxOrdersCtx *ctx, int verb)
Function build MG structure.
Kronecker Delta class symmetric.
#define CATCH_ERRORS
Catch errors.
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
@ BODYFORCESSET
block name is "BODY_FORCES"
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
@ 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.
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
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 DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
MoFEMErrorCode DMRegister_MGViaApproxOrders(const char sname[])
Register DM for Multi-Grid via approximation orders.
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
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 add_ents_to_finite_element_by_dim(const EntityHandle entities, const int dim, const std::string &name, const bool recursive=true)=0
add entities to finite element
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_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
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 build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode add_ents_to_field_by_dim(const Range &ents, const int dim, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
virtual MoFEMErrorCode get_finite_element_entities_by_dimension(const std::string name, int dim, Range &ents) const =0
get entities in the finite element by dimension
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
virtual bool check_field(const std::string &name) const =0
check if field is in database
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
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_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet 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 load_series_data(const std::string &serie_name, const int step_number)
virtual bool check_series(const std::string &name) const
check if series is in database
#define _IT_SERIES_STEPS_BY_NAME_FOR_LOOP_(RECORDER, NAME, IT)
loop over recorded series step
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
const FTensor::Tensor2< T, Dim, Dim > Vec
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
UBlasMatrix< double > MatrixDouble
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
auto createSmartDM
Creates smart DM object.
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
SmartPetscObj< Mat > smartMatDuplicate(Mat &mat, MatDuplicateOption op)
DeprecatedCoreInterface Interface
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
const double D
diffusivity
MoFEMErrorCode setBlocks()
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
virtual moab::Interface & get_moab()=0
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
Interface for managing meshsets containing materials and boundary conditions.
MoFEMErrorCode printForceSet() const
print meshsets with force boundary conditions data structure
MoFEMErrorCode printMaterialsSet() const
print meshsets with material data structure set on it
MoFEMErrorCode printDisplacementSet() const
print meshsets with displacement boundary conditions data structure
Calculate inverse of jacobian for face element.
Transform local reference derivatives of shape functions to global derivatives.
Projection of edge entities with one mid-node on hierarchical basis.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Set data structures of MG pre-conditioner via approximation orders.
Operator post-procesing stresses for Hook isotropic material.
Implentation of thermal stress element.