55using namespace boost::numeric;
57static char help[] =
"-my_block_config set block data\n"
58 "-my_order approximation order\n"
59 "-my_is_partitioned set if mesh is partitioned\n"
86 FatPrismElementForcesAndSourcesCore;
96int main(
int argc,
char *argv[]) {
100 "-pc_factor_mat_solver_type mumps \n"
101 "-mat_mumps_icntl_20 0 \n"
103 "-snes_type newtonls \n"
104 "-snes_linesearch_type basic \n"
109 "-ts_type beuler \n";
112 if (!
static_cast<bool>(ifstream(
param_file))) {
113 std::ofstream file(
param_file.c_str(), std::ios::ate);
114 if (file.is_open()) {
121 auto core_log = logging::core::get();
129 PetscBool flg_block_config, flg_file;
131 char block_config_file[255];
132 PetscInt test_nb = 0;
134 PetscBool is_partitioned = PETSC_FALSE;
135 PetscBool is_calculating_frequency = PETSC_FALSE;
136 PetscBool is_post_proc_volume = PETSC_TRUE;
139 enum bases { LEGENDRE, LOBATTO, BERNSTEIN_BEZIER, JACOBI, LASBASETOP };
140 const char *list_bases[] = {
"legendre",
"lobatto",
"bernstein_bezier",
142 PetscInt choice_base_value = LOBATTO;
145 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Elastic Config",
"none");
147 CHKERR PetscOptionsString(
"-my_file",
"mesh file name",
"",
"mesh.h5m",
150 CHKERR PetscOptionsInt(
"-my_order",
"default approximation order",
"",
153 CHKERR PetscOptionsEList(
"-base",
"approximation base",
"", list_bases,
154 LASBASETOP, list_bases[choice_base_value],
155 &choice_base_value, PETSC_NULL);
157 CHKERR PetscOptionsInt(
"-is_atom_test",
"ctest number",
"", test_nb,
158 &test_nb, PETSC_NULL);
160 CHKERR PetscOptionsBool(
"-my_is_partitioned",
161 "set if mesh is partitioned (this result that each "
162 "process keeps only one part of the mesh)",
163 "", is_partitioned, &is_partitioned, PETSC_NULL);
165 CHKERR PetscOptionsString(
"-my_block_config",
"elastic configure file name",
166 "",
"block_conf.in", block_config_file, 255,
170 "-my_is_calculating_frequency",
"set if frequency will be calculated",
171 "", is_calculating_frequency, &is_calculating_frequency, PETSC_NULL);
173 CHKERR PetscOptionsBool(
"-my_is_post_proc_volume",
174 "if true post proc volume",
"", is_post_proc_volume,
175 &is_post_proc_volume, PETSC_NULL);
177 ierr = PetscOptionsEnd();
181 if (flg_file != PETSC_TRUE) {
182 SETERRQ(PETSC_COMM_SELF, 1,
"*** ERROR -my_file (MESH FILE NEEDED)");
186 moab::Core mb_instance;
187 moab::Interface &moab = mb_instance;
193 MPI_Comm moab_comm_world;
194 MPI_Comm_dup(PETSC_COMM_WORLD, &moab_comm_world);
195 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
197 pcomm =
new ParallelComm(&moab, moab_comm_world);
200 if (is_partitioned == PETSC_TRUE) {
204 option =
"PARALLEL=READ_PART;"
205 "PARALLEL_RESOLVE_SHARED_ENTS;"
206 "PARTITION=PARALLEL_PARTITION;";
229 bool mesh_has_tets =
false;
230 bool mesh_has_prisms =
false;
235 CHKERR moab.get_number_entities_by_type(0, MBTET, nb_tets,
true);
236 CHKERR moab.get_number_entities_by_type(0, MBHEX, nb_hexs,
true);
237 CHKERR moab.get_number_entities_by_type(0, MBPRISM, nb_prisms,
true);
239 mesh_has_tets = (nb_tets + nb_hexs) > 0;
240 mesh_has_prisms = nb_prisms > 0;
254 if (
bit->getName().compare(0, 3,
"ROD") == 0) {
262 switch (choice_base_value) {
269 case BERNSTEIN_BEZIER:
294 CHKERR m_field.
get_moab().get_entities_by_type(0, MBEDGE, all_edges,
true);
297 Range edges_in_simple_rod;
299 if (
bit->getName().compare(0, 3,
"ROD") == 0) {
302 MBEDGE, edges,
true);
303 edges_in_simple_rod.merge(edges);
314 Range edges_to_set_order;
315 edges_to_set_order = subtract(all_edges, edges_in_simple_rod);
335 auto setting_second_order_geometry = [&m_field]() {
341 moab::Interface::UNION);
358 CHKERR setting_second_order_geometry();
362 std::map<int, BlockOptionData> block_data;
363 auto setting_blocks_data_and_order_from_config_file =
364 [&m_field, &moab, &block_data, flg_block_config, block_config_file,
365 order](boost::shared_ptr<std::map<int, BlockData>> &block_sets_ptr) {
367 if (flg_block_config) {
368 ifstream ini_file(block_config_file);
369 po::variables_map vm;
370 po::options_description config_file_options;
373 std::ostringstream str_order;
374 str_order <<
"block_" << it->getMeshsetId()
375 <<
".displacement_order";
376 config_file_options.add_options()(
377 str_order.str().c_str(),
378 po::value<int>(&block_data[it->getMeshsetId()].oRder)
379 ->default_value(
order));
381 std::ostringstream str_cond;
382 str_cond <<
"block_" << it->getMeshsetId() <<
".young_modulus";
383 config_file_options.add_options()(
384 str_cond.str().c_str(),
385 po::value<double>(&block_data[it->getMeshsetId()].yOung)
386 ->default_value(-1));
388 std::ostringstream str_capa;
389 str_capa <<
"block_" << it->getMeshsetId() <<
".poisson_ratio";
390 config_file_options.add_options()(
391 str_capa.str().c_str(),
392 po::value<double>(&block_data[it->getMeshsetId()].pOisson)
393 ->default_value(-2));
395 std::ostringstream str_init_temp;
396 str_init_temp <<
"block_" << it->getMeshsetId()
397 <<
".initial_temperature";
398 config_file_options.add_options()(
399 str_init_temp.str().c_str(),
400 po::value<double>(&block_data[it->getMeshsetId()].initTemp)
403 po::parsed_options parsed =
404 parse_config_file(ini_file, config_file_options,
true);
409 if (block_data[it->getMeshsetId()].oRder == -1)
411 if (block_data[it->getMeshsetId()].oRder ==
order)
413 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Set block %d order to %d",
415 block_data[it->getMeshsetId()].oRder);
417 CHKERR moab.get_entities_by_handle(it->getMeshset(), block_ents,
419 Range ents_to_set_order;
420 CHKERR moab.get_adjacencies(block_ents, 3,
false,
422 moab::Interface::UNION);
423 ents_to_set_order = ents_to_set_order.subset_by_dimension(3);
424 CHKERR moab.get_adjacencies(block_ents, 2,
false,
426 moab::Interface::UNION);
427 CHKERR moab.get_adjacencies(block_ents, 1,
false,
429 moab::Interface::UNION);
434 ents_to_set_order,
"DISPLACEMENT",
435 block_data[it->getMeshsetId()].oRder);
437 std::vector<std::string> additional_parameters;
438 additional_parameters =
439 collect_unrecognized(parsed.options, po::include_positional);
440 for (std::vector<std::string>::iterator vit =
441 additional_parameters.begin();
442 vit != additional_parameters.end(); vit++) {
443 MOFEM_LOG_C(
"ELASTIC", Sev::warning,
"Unrecognized option %s",
452 const int id = it->getMeshsetId();
453 auto &bd = (*block_sets_ptr)[id];
454 if (block_data[
id].yOung > 0)
455 bd.E = block_data[id].yOung;
456 if (block_data[
id].pOisson >= -1)
457 bd.PoissonRatio = block_data[id].pOisson;
458 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Block %d",
id);
459 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"\tYoung modulus %3.4g", bd.E);
460 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"\tPoisson ratio %3.4g",
469 boost::shared_ptr<std::map<int, HookeElement::BlockData>> block_sets_ptr =
470 boost::make_shared<std::map<int, HookeElement::BlockData>>();
471 CHKERR HookeElement::setBlocks(m_field, block_sets_ptr);
472 CHKERR setting_blocks_data_and_order_from_config_file(block_sets_ptr);
474 boost::shared_ptr<std::map<int, MassBlockData>> mass_block_sets_ptr =
475 boost::make_shared<std::map<int, MassBlockData>>();
479 boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
481 boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
482 fe_lhs_ptr->getRuleHook =
VolRule();
483 fe_rhs_ptr->getRuleHook =
VolRule();
490 boost::shared_ptr<ForcesAndSourcesCore> prism_fe_lhs_ptr(
492 boost::shared_ptr<ForcesAndSourcesCore> prism_fe_rhs_ptr(
495 CHKERR HookeElement::addElasticElement(m_field, block_sets_ptr,
"ELASTIC",
497 "MESH_NODE_POSITIONS",
false);
499 auto add_skin_element_for_post_processing = [&]() {
501 Range elastic_element_ents;
503 "ELASTIC", 3, elastic_element_ents);
506 CHKERR skin.find_skin(0, elastic_element_ents,
false, skin_faces);
508 if (is_partitioned) {
509 CHKERR pcomm->filter_pstatus(skin_faces,
510 PSTATUS_SHARED | PSTATUS_MULTISHARED,
511 PSTATUS_NOT, -1, &proc_skin);
513 proc_skin = skin_faces;
523 "POST_PROC_SKIN",
"MESH_NODE_POSITIONS");
528 CHKERR add_skin_element_for_post_processing();
530 auto data_at_pts = boost::make_shared<HookeElement::DataAtIntegrationPts>();
532 CHKERR HookeElement::setOperators(fe_lhs_ptr, fe_rhs_ptr, block_sets_ptr,
533 "DISPLACEMENT",
"MESH_NODE_POSITIONS",
534 false,
true, MBTET, data_at_pts);
536 if (mesh_has_prisms) {
537 CHKERR HookeElement::setOperators(
538 prism_fe_lhs_ptr, prism_fe_rhs_ptr, block_sets_ptr,
"DISPLACEMENT",
539 "MESH_NODE_POSITIONS",
false,
true, MBPRISM, data_at_pts);
544 auto thermal_strain =
547 constexpr double alpha = 1;
551 t_thermal_strain(
i,
j) = alpha * t_coords(2) *
t_kd(
i,
j);
552 return t_thermal_strain;
555 fe_rhs_ptr->getOpPtrVector().push_back(
556 new HookeElement::OpAnalyticalInternalStrain_dx<0>(
557 "DISPLACEMENT", data_at_pts, thermal_strain));
560 boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_mass_ptr(
563 for (
auto &sit : *block_sets_ptr) {
564 for (
auto &mit : *mass_block_sets_ptr) {
565 fe_mass_ptr->getOpPtrVector().push_back(
566 new HookeElement::OpCalculateMassMatrix(
"DISPLACEMENT",
567 "DISPLACEMENT", sit.second,
568 mit.second, data_at_pts));
575 "MESH_NODE_POSITIONS");
579 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr(
581 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr(
585 fe_spring_rhs_ptr,
"DISPLACEMENT",
586 "MESH_NODE_POSITIONS");
591 "MESH_NODE_POSITIONS");
598 boost::shared_ptr<EdgeEle> fe_simple_rod_lhs_ptr(
new EdgeEle(m_field));
599 boost::shared_ptr<EdgeEle> fe_simple_rod_rhs_ptr(
new EdgeEle(m_field));
603 m_field, fe_simple_rod_lhs_ptr, fe_simple_rod_rhs_ptr,
"DISPLACEMENT",
604 "MESH_NODE_POSITIONS");
616 "MESH_NODE_POSITIONS");
621 CHKERR m_field.
get_moab().get_entities_by_dimension(it->meshset, 3, tets,
644 bool add_temp_field =
false;
646 if (block_data[it->getMeshsetId()].initTemp != 0) {
647 add_temp_field =
true;
651 if (add_temp_field) {
661 "ELASTIC",
"DISPLACEMENT",
"TEMP");
670 "MESH_NODE_POSITIONS");
674 if (block_data[it->getMeshsetId()].initTemp != 0) {
676 "Set block %d temperature to %3.2g\n", it->getMeshsetId(),
677 block_data[it->getMeshsetId()].initTemp);
679 CHKERR moab.get_entities_by_handle(it->meshset, block_ents,
true);
681 CHKERR moab.get_connectivity(block_ents, vertices,
true);
683 block_data[it->getMeshsetId()].initTemp, MBVERTEX, vertices,
699 auto dm =
createDM(PETSC_COMM_WORLD,
"MOFEM");
701 CHKERR DMSetFromOptions(dm);
722 CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
726 if (is_calculating_frequency == PETSC_TRUE) {
728 CHKERR MatSetOption(Mij, MAT_SPD, PETSC_TRUE);
733 fe_spring_lhs_ptr->ksp_B = Aij;
734 fe_spring_rhs_ptr->ksp_f =
F;
737 fe_simple_rod_lhs_ptr->ksp_B = Aij;
738 fe_simple_rod_rhs_ptr->ksp_f =
F;
742 CHKERR VecGhostUpdateBegin(
F, INSERT_VALUES, SCATTER_FORWARD);
743 CHKERR VecGhostUpdateEnd(
F, INSERT_VALUES, SCATTER_FORWARD);
745 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
746 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
748 CHKERR MatZeroEntries(Aij);
760 auto dirichlet_bc_ptr = boost::make_shared<DirichletDisplacementBc>(
761 m_field,
"DISPLACEMENT", Aij, D0,
F);
765 dirichlet_bc_ptr->snes_ctx = FEMethod::CTX_SNESNONE;
766 dirichlet_bc_ptr->ts_ctx = FEMethod::CTX_TSNONE;
769 CHKERR VecZeroEntries(D0);
770 CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
771 CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
779 CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
780 CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
789 fe_rhs_ptr->snes_f =
F;
790 prism_fe_rhs_ptr->snes_f =
F;
791 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Assemble external force vector ...";
794 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
796 fe_lhs_ptr->snes_B = Aij;
797 prism_fe_lhs_ptr->snes_B = Aij;
798 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Calculate stiffness matrix ...";
801 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
811 if (is_calculating_frequency == PETSC_TRUE) {
813 fe_mass_ptr->snes_B = Mij;
814 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Calculate mass matrix ...";
816 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
824 boost::ptr_map<std::string, NeumannForcesSurface> neumann_forces;
829 boost::ptr_map<std::string, NeumannForcesSurface>::iterator mit =
830 neumann_forces.begin();
831 for (; mit != neumann_forces.end(); mit++) {
833 &mit->second->getLoopFe());
837 boost::ptr_map<std::string, NodalForce> nodal_forces;
842 boost::ptr_map<std::string, NodalForce>::iterator fit =
843 nodal_forces.begin();
844 for (; fit != nodal_forces.end(); fit++) {
846 &fit->second->getLoopFe());
850 boost::ptr_map<std::string, EdgeForce> edge_forces;
854 auto fit = edge_forces.begin();
855 for (; fit != edge_forces.end(); fit++) {
856 auto &fe = fit->second->getLoopFe();
873 "DISPLACEMENT",
F,
false,
true);
881 PetscViewerPushFormat(
882 PETSC_VIEWER_STDOUT_SELF,
883 PETSC_VIEWER_ASCII_MATLAB);
890 if (is_calculating_frequency == PETSC_TRUE) {
891 CHKERR MatAssemblyBegin(Mij, MAT_FINAL_ASSEMBLY);
892 CHKERR MatAssemblyEnd(Mij, MAT_FINAL_ASSEMBLY);
898 CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
899 CHKERR VecGhostUpdateBegin(
F, ADD_VALUES, SCATTER_REVERSE);
900 CHKERR VecGhostUpdateEnd(
F, ADD_VALUES, SCATTER_REVERSE);
906 auto solver =
createKSP(PETSC_COMM_WORLD);
907 CHKERR KSPSetDM(solver, dm);
908 CHKERR KSPSetFromOptions(solver);
909 CHKERR KSPSetOperators(solver, Aij, Aij);
914 PetscBool same = PETSC_FALSE;
916 CHKERR KSPGetPC(solver, &pc);
917 PetscObjectTypeCompare((PetscObject)pc, PCMG, &same);
921 CHKERR PCSetFromOptions(pc);
924 CHKERR KSPSetDMActive(solver, PETSC_FALSE);
927 CHKERR KSPSetInitialGuessKnoll(solver, PETSC_FALSE);
928 CHKERR KSPSetInitialGuessNonzero(solver, PETSC_TRUE);
932 auto set_post_proc_skin = [&](
auto &post_proc_skin) {
936 CHKERR post_proc_skin.generateReferenceElementMesh();
937 CHKERR post_proc_skin.addFieldValuesPostProc(
"DISPLACEMENT");
938 CHKERR post_proc_skin.addFieldValuesPostProc(
"MESH_NODE_POSITIONS");
939 CHKERR post_proc_skin.addFieldValuesGradientPostProcOnSkin(
940 "DISPLACEMENT",
"ELASTIC", data_at_pts->hMat,
true);
941 CHKERR post_proc_skin.addFieldValuesGradientPostProcOnSkin(
943 "MESH_NODE_POSITIONS",
"ELASTIC", data_at_pts->HMat,
false);
944 post_proc_skin.getOpPtrVector().push_back(
945 new HookeElement::OpPostProcHookeElement<
947 "DISPLACEMENT", data_at_pts, *block_sets_ptr,
948 post_proc_skin.postProcMesh, post_proc_skin.mapGaussPts,
true,
953 auto set_post_proc_tets = [&](
auto &post_proc) {
956 CHKERR post_proc.generateReferenceElementMesh();
959 CHKERR post_proc.addFieldValuesPostProc(
"DISPLACEMENT");
960 CHKERR post_proc.addFieldValuesPostProc(
"MESH_NODE_POSITIONS");
961 CHKERR post_proc.addFieldValuesGradientPostProc(
"DISPLACEMENT");
964 m_field, post_proc.postProcMesh, post_proc.mapGaussPts,
965 "DISPLACEMENT", post_proc.commonData, block_sets_ptr.get()));
969 auto set_post_proc_edge = [&](
auto &post_proc_edge) {
971 CHKERR post_proc_edge.generateReferenceElementMesh();
972 CHKERR post_proc_edge.addFieldValuesPostProc(
"DISPLACEMENT");
976 auto set_post_proc_prisms = [&](
auto &prism_post_proc) {
978 CHKERR prism_post_proc.generateReferenceElementMesh();
979 boost::shared_ptr<MatrixDouble> inv_jac_ptr(
new MatrixDouble);
980 prism_post_proc.getOpPtrVector().push_back(
982 prism_post_proc.getOpPtrVector().push_back(
984 CHKERR prism_post_proc.addFieldValuesPostProc(
"DISPLACEMENT");
985 CHKERR prism_post_proc.addFieldValuesPostProc(
"MESH_NODE_POSITIONS");
986 CHKERR prism_post_proc.addFieldValuesGradientPostProc(
"DISPLACEMENT");
988 m_field, prism_post_proc.postProcMesh, prism_post_proc.mapGaussPts,
989 "DISPLACEMENT", prism_post_proc.commonData, block_sets_ptr.get()));
998 CHKERR set_post_proc_skin(post_proc_skin);
999 CHKERR set_post_proc_tets(post_proc);
1000 CHKERR set_post_proc_prisms(prism_post_proc);
1001 CHKERR set_post_proc_edge(post_proc_edge);
1008 CHKERR VecDuplicate(
F, &F_thermal);
1013 "DISPLACEMENT",
"TEMP", F_thermal);
1031 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Process step %d",
1032 sit->get_step_number());
1035 sit->get_step_number());
1037 CHKERR VecZeroEntries(F_thermal);
1038 CHKERR VecGhostUpdateBegin(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1039 CHKERR VecGhostUpdateEnd(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1050 CHKERR VecAssemblyBegin(F_thermal);
1051 CHKERR VecAssemblyEnd(F_thermal);
1053 CHKERR VecGhostUpdateBegin(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1054 CHKERR VecGhostUpdateEnd(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1058 CHKERR VecNorm(
F, NORM_2, &nrm_F);
1060 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"norm2 F = %6.4e", nrm_F);
1061 PetscReal nrm_F_thermal;
1062 CHKERR VecNorm(F_thermal, NORM_2, &nrm_F_thermal);
1063 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"norm2 F_thermal = %6.4e",
1066 CHKERR VecScale(F_thermal, -1);
1068 CHKERR VecAXPY(F_thermal, 1,
F);
1071 dirichlet_bc_ptr->snes_x =
D;
1072 dirichlet_bc_ptr->snes_f = F_thermal;
1076 CHKERR KSPSolve(solver, F_thermal,
D);
1079 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
1080 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
1090 if (is_post_proc_volume == PETSC_TRUE) {
1091 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Write output file ...";
1093 std::ostringstream o1;
1094 o1 <<
"out_" << sit->step_number <<
".h5m";
1097 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done ...";
1100 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Write output file skin ...";
1103 std::ostringstream o1_skin;
1104 o1_skin <<
"out_skin" << sit->step_number <<
".h5m";
1107 MOFEM_LOG(
"POST_PROC_SKIN", Sev::inform) <<
"done ...";
1112 CHKERR VecZeroEntries(F_thermal);
1113 CHKERR VecGhostUpdateBegin(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1114 CHKERR VecGhostUpdateEnd(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1121 CHKERR VecAssemblyBegin(F_thermal);
1122 CHKERR VecAssemblyEnd(F_thermal);
1125 CHKERR VecGhostUpdateBegin(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1126 CHKERR VecGhostUpdateEnd(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1130 CHKERR VecNorm(
F, NORM_2, &nrm_F);
1132 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"norm2 F = %6.4e", nrm_F);
1133 PetscReal nrm_F_thermal;
1134 CHKERR VecNorm(F_thermal, NORM_2, &nrm_F_thermal);
1136 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"norm2 F_thermal = %6.4e",
1140 CHKERR VecScale(F_thermal, -1);
1141 CHKERR VecAXPY(F_thermal, 1,
F);
1144 dirichlet_bc_ptr->snes_x =
D;
1145 dirichlet_bc_ptr->snes_f = F_thermal;
1149 CHKERR KSPSolve(solver, F_thermal,
D);
1153 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
1154 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
1158 if (is_post_proc_volume == PETSC_TRUE) {
1159 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Write output file ...";
1164 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
1167 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Write output file skin ...";
1171 MOFEM_LOG(
"POST_PROC_SKIN", Sev::inform) <<
"done";
1175 CHKERR VecDestroy(&F_thermal);
1188 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
1189 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
1193 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Post-process start ...";
1194 if (is_post_proc_volume == PETSC_TRUE) {
1200 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
1203 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Write output file ...";
1204 if (mesh_has_tets) {
1205 if (is_post_proc_volume == PETSC_TRUE) {
1212 if (mesh_has_prisms) {
1216 if (!edges_in_simple_rod.empty())
1219 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
1222 if (is_calculating_frequency == PETSC_TRUE) {
1225 VecDuplicate(
D, &u1);
1228 CHKERR VecDot(u1,
D, &mode_mass);
1229 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Mode mass %6.4e\n", mode_mass);
1232 VecDuplicate(
D, &v1);
1235 double mode_stiffness;
1236 CHKERR VecDot(v1,
D, &mode_stiffness);
1237 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Mode stiffness %6.4e\n",
1241 double pi = 3.14159265359;
1242 frequency = std::sqrt(mode_stiffness / mode_mass) / (2 * pi);
1243 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Frequency %6.4e", frequency);
1247 auto calculate_strain_energy = [&]() {
1251 CHKERR HookeElement::calculateEnergy(dm, block_sets_ptr,
"DISPLACEMENT",
1252 "MESH_NODE_POSITIONS",
false,
true,
1257 CHKERR VecSum(v_energy, &energy);
1258 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Elastic energy %6.4e", energy);
1262 if (fabs(energy - 17.129) > 1e-3)
1264 "atom test diverged!");
1267 if (fabs(energy - 5.6475e-03) > 1e-4)
1269 "atom test diverged!");
1272 if (fabs(energy - 7.4679e-03) > 1e-4)
1274 "atom test diverged!");
1277 if (fabs(energy - 2.4992e+00) > 1e-3)
1279 "atom test diverged!");
1284 CHKERR VecMin(
D, PETSC_NULL, &min);
1285 constexpr double expected_val = 0.10001;
1286 if (fabs(min + expected_val) > 1e-10)
1288 "atom test diverged! %3.4e != %3.4e", min, expected_val);
1291 if (fabs(energy - 4.7416e-04) > 1e-8)
1293 "atom test diverged!");
1301 CHKERR calculate_strain_energy();
1303 MPI_Comm_free(&moab_comm_world);
const std::string default_options
Implementation of linear elastic material.
#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 DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
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 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 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_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
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 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
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"
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#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
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
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
auto createKSP(MPI_Comm comm)
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
SmartPetscObj< Mat > matDuplicate(Mat mat, MatDuplicateOption op)
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
MoFEMErrorCode addBlock(const std::string field_name, Vec F, int ms_id)
data for calculation inertia forces
MoFEMErrorCode setBlocks()
MyTriangleFE & getLoopFe()
MoFEMErrorCode setNeumannFluidPressureFiniteElementOperators(string field_name, Vec F, bool allow_negative_pressure=true, bool ho_geometry=false)
MoFEMErrorCode addNeumannFluidPressureBCElements(const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
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.
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
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.
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Volume finite element base.
data for calculation heat conductivity and heat capacity elements
Set data structures of MG pre-conditioner via approximation orders.
Operator post-procesing stresses for Hook isotropic material.
int getRuleThroughThickness(int order)
int getRuleTrianglesOnly(int order)
Implentation of thermal stress element.
MyVolumeFE & getLoopThermalStressRhs()
MoFEMErrorCode addThermalStressElement(const std::string fe_name, const std::string field_name, const std::string thermal_field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
MoFEMErrorCode setThermalStressRhsOperators(string field_name, string thermal_field_name, Vec &F, int verb=0)
int operator()(int, int, int order) const