96int main(
int argc,
char *argv[]) {
98 const string default_options =
"-ksp_type gmres \n"
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";
111 string param_file =
"param_file.petsc";
112 if (!
static_cast<bool>(ifstream(param_file))) {
113 std::ofstream file(param_file.c_str(), std::ios::ate);
114 if (file.is_open()) {
115 file << default_options;
121 auto core_log = logging::core::get();
134 PetscBool flg_block_config, flg_file;
136 char block_config_file[255];
137 PetscInt test_nb = 0;
139 PetscBool is_partitioned = PETSC_FALSE;
140 PetscBool is_calculating_frequency = PETSC_FALSE;
141 PetscBool is_post_proc_volume = PETSC_TRUE;
144 enum bases { LEGENDRE, LOBATTO, BERNSTEIN_BEZIER, JACOBI, LASBASETOP };
145 const char *list_bases[] = {
"legendre",
"lobatto",
"bernstein_bezier",
147 PetscInt choice_base_value = LOBATTO;
150 PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Elastic Config",
"none");
151 CHKERR PetscOptionsString(
"-my_file",
"mesh file name",
"",
"mesh.h5m",
154 CHKERR PetscOptionsInt(
"-my_order",
"default approximation order",
"",
157 CHKERR PetscOptionsEList(
"-base",
"approximation base",
"", list_bases,
158 LASBASETOP, list_bases[choice_base_value],
159 &choice_base_value, PETSC_NULLPTR);
161 CHKERR PetscOptionsInt(
"-is_atom_test",
"ctest number",
"", test_nb,
162 &test_nb, PETSC_NULLPTR);
164 CHKERR PetscOptionsBool(
"-my_is_partitioned",
165 "set if mesh is partitioned (this result that each "
166 "process keeps only one part of the mesh)",
167 "", is_partitioned, &is_partitioned, PETSC_NULLPTR);
169 CHKERR PetscOptionsString(
"-my_block_config",
"elastic configure file name",
170 "",
"block_conf.in", block_config_file, 255,
174 "-my_is_calculating_frequency",
"set if frequency will be calculated",
175 "", is_calculating_frequency, &is_calculating_frequency, PETSC_NULLPTR);
177 CHKERR PetscOptionsBool(
"-my_is_post_proc_volume",
178 "if true post proc volume",
"", is_post_proc_volume,
179 &is_post_proc_volume, PETSC_NULLPTR);
184 if (flg_file != PETSC_TRUE) {
185 SETERRQ(PETSC_COMM_SELF, 1,
"*** ERROR -my_file (MESH FILE NEEDED)");
189 moab::Core mb_instance;
190 moab::Interface &moab = mb_instance;
196 MPI_Comm moab_comm_world;
197 MPI_Comm_dup(PETSC_COMM_WORLD, &moab_comm_world);
198 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
200 pcomm =
new ParallelComm(&moab, moab_comm_world);
203 if (is_partitioned == PETSC_TRUE) {
207 option =
"PARALLEL=READ_PART;"
208 "PARALLEL_RESOLVE_SHARED_ENTS;"
209 "PARTITION=PARALLEL_PARTITION;";
232 bool mesh_has_tets =
false;
233 bool mesh_has_prisms =
false;
238 CHKERR moab.get_number_entities_by_type(0, MBTET, nb_tets,
true);
239 CHKERR moab.get_number_entities_by_type(0, MBHEX, nb_hexs,
true);
240 CHKERR moab.get_number_entities_by_type(0, MBPRISM, nb_prisms,
true);
242 mesh_has_tets = (nb_tets + nb_hexs) > 0;
243 mesh_has_prisms = nb_prisms > 0;
257 if (
bit->getName().compare(0, 3,
"ROD") == 0) {
265 switch (choice_base_value) {
272 case BERNSTEIN_BEZIER:
297 CHKERR m_field.
get_moab().get_entities_by_type(0, MBEDGE, all_edges,
true);
300 Range edges_in_simple_rod;
302 if (
bit->getName().compare(0, 3,
"ROD") == 0) {
305 MBEDGE, edges,
true);
306 edges_in_simple_rod.merge(edges);
317 Range edges_to_set_order;
318 edges_to_set_order = subtract(all_edges, edges_in_simple_rod);
338 auto setting_second_order_geometry = [&m_field]() {
344 moab::Interface::UNION);
361 CHKERR setting_second_order_geometry();
365 std::map<int, BlockOptionData> block_data;
366 auto setting_blocks_data_and_order_from_config_file =
367 [&m_field, &moab, &block_data, flg_block_config, block_config_file,
368 order](boost::shared_ptr<std::map<int, BlockData>> &block_sets_ptr) {
370 if (flg_block_config) {
371 ifstream ini_file(block_config_file);
372 po::variables_map vm;
373 po::options_description config_file_options;
376 std::ostringstream str_order;
377 str_order <<
"block_" << it->getMeshsetId()
378 <<
".displacement_order";
379 config_file_options.add_options()(
380 str_order.str().c_str(),
381 po::value<int>(&block_data[it->getMeshsetId()].oRder)
382 ->default_value(
order));
384 std::ostringstream str_cond;
385 str_cond <<
"block_" << it->getMeshsetId() <<
".young_modulus";
386 config_file_options.add_options()(
387 str_cond.str().c_str(),
388 po::value<double>(&block_data[it->getMeshsetId()].yOung)
389 ->default_value(-1));
391 std::ostringstream str_capa;
392 str_capa <<
"block_" << it->getMeshsetId() <<
".poisson_ratio";
393 config_file_options.add_options()(
394 str_capa.str().c_str(),
395 po::value<double>(&block_data[it->getMeshsetId()].pOisson)
396 ->default_value(-2));
398 std::ostringstream str_init_temp;
399 str_init_temp <<
"block_" << it->getMeshsetId()
400 <<
".initial_temperature";
401 config_file_options.add_options()(
402 str_init_temp.str().c_str(),
403 po::value<double>(&block_data[it->getMeshsetId()].initTemp)
406 po::parsed_options parsed =
407 parse_config_file(ini_file, config_file_options,
true);
412 if (block_data[it->getMeshsetId()].oRder == -1)
414 if (block_data[it->getMeshsetId()].oRder ==
order)
416 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Set block %d order to %d",
418 block_data[it->getMeshsetId()].oRder);
420 CHKERR moab.get_entities_by_handle(it->getMeshset(), block_ents,
422 Range ents_to_set_order;
423 CHKERR moab.get_adjacencies(block_ents, 3,
false,
425 moab::Interface::UNION);
426 ents_to_set_order = ents_to_set_order.subset_by_dimension(3);
427 CHKERR moab.get_adjacencies(block_ents, 2,
false,
429 moab::Interface::UNION);
430 CHKERR moab.get_adjacencies(block_ents, 1,
false,
432 moab::Interface::UNION);
437 ents_to_set_order,
"DISPLACEMENT",
438 block_data[it->getMeshsetId()].oRder);
440 std::vector<std::string> additional_parameters;
441 additional_parameters =
442 collect_unrecognized(parsed.options, po::include_positional);
443 for (std::vector<std::string>::iterator vit =
444 additional_parameters.begin();
445 vit != additional_parameters.end(); vit++) {
446 MOFEM_LOG_C(
"ELASTIC", Sev::warning,
"Unrecognized option %s",
455 const int id = it->getMeshsetId();
456 auto &bd = (*block_sets_ptr)[id];
457 if (block_data[
id].yOung > 0)
458 bd.E = block_data[id].yOung;
459 if (block_data[
id].pOisson >= -1)
460 bd.PoissonRatio = block_data[id].pOisson;
461 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Block %d",
id);
462 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"\tYoung modulus %3.4g", bd.E);
463 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"\tPoisson ratio %3.4g",
472 boost::shared_ptr<std::map<int, HookeElement::BlockData>> block_sets_ptr =
473 boost::make_shared<std::map<int, HookeElement::BlockData>>();
474 CHKERR HookeElement::setBlocks(m_field, block_sets_ptr);
475 CHKERR setting_blocks_data_and_order_from_config_file(block_sets_ptr);
477 boost::shared_ptr<std::map<int, MassBlockData>> mass_block_sets_ptr =
478 boost::make_shared<std::map<int, MassBlockData>>();
482 boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
484 boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
485 fe_lhs_ptr->getRuleHook =
VolRule();
486 fe_rhs_ptr->getRuleHook =
VolRule();
493 boost::shared_ptr<ForcesAndSourcesCore> prism_fe_lhs_ptr(
495 boost::shared_ptr<ForcesAndSourcesCore> prism_fe_rhs_ptr(
498 CHKERR HookeElement::addElasticElement(m_field, block_sets_ptr,
"ELASTIC",
500 "MESH_NODE_POSITIONS",
false);
502 auto add_skin_element_for_post_processing = [&]() {
504 Range elastic_element_ents;
506 "ELASTIC", 3, elastic_element_ents);
509 CHKERR skin.find_skin(0, elastic_element_ents,
false, skin_faces);
511 if (is_partitioned) {
512 CHKERR pcomm->filter_pstatus(skin_faces,
513 PSTATUS_SHARED | PSTATUS_MULTISHARED,
514 PSTATUS_NOT, -1, &proc_skin);
516 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>();
543 CHKERR HookeElement::setOperators(fe_lhs_ptr, fe_rhs_ptr, block_sets_ptr,
544 "DISPLACEMENT",
"MESH_NODE_POSITIONS",
545 false,
true, MBTET, data_at_pts);
547 if (mesh_has_prisms) {
548 CHKERR HookeElement::setOperators(
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;
562 t_thermal_strain(
i,
j) = alpha * t_coords(2) *
t_kd(
i,
j);
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,
655 bool add_temp_field =
false;
657 if (block_data[it->getMeshsetId()].initTemp != 0) {
658 add_temp_field =
true;
662 if (add_temp_field) {
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,
710 auto dm =
createDM(PETSC_COMM_WORLD,
"MOFEM");
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();
884 "DISPLACEMENT",
F,
false,
true);
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);
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(
953 "MESH_NODE_POSITIONS",
"ELASTIC", data_at_pts->HMat,
false);
955 CHKERR post_proc_skin.addFieldValuesPostProc(
"TEMP");
956 CHKERR post_proc_skin.addFieldValuesGradientPostProc(
"TEMP");
958 post_proc_skin.getOpPtrVector().push_back(
959 new HookeElement::OpPostProcHookeElement<
961 "DISPLACEMENT", data_at_pts, *block_sets_ptr,
962 post_proc_skin.postProcMesh, post_proc_skin.mapGaussPts,
true,
967 auto set_post_proc_tets = [&](
auto &post_proc) {
970 CHKERR post_proc.generateReferenceElementMesh();
973 CHKERR post_proc.addFieldValuesPostProc(
"DISPLACEMENT");
974 CHKERR post_proc.addFieldValuesPostProc(
"MESH_NODE_POSITIONS");
975 CHKERR post_proc.addFieldValuesGradientPostProc(
"DISPLACEMENT");
977 CHKERR post_proc.addFieldValuesPostProc(
"TEMP");
978 CHKERR post_proc.addFieldValuesGradientPostProc(
"TEMP");
982 m_field, post_proc.postProcMesh, post_proc.mapGaussPts,
983 "DISPLACEMENT", post_proc.commonData, block_sets_ptr.get()));
987 auto set_post_proc_edge = [&](
auto &post_proc_edge) {
989 CHKERR post_proc_edge.generateReferenceElementMesh();
990 CHKERR post_proc_edge.addFieldValuesPostProc(
"DISPLACEMENT");
994 auto set_post_proc_prisms = [&](
auto &prism_post_proc) {
996 CHKERR prism_post_proc.generateReferenceElementMesh();
997 boost::shared_ptr<MatrixDouble> inv_jac_ptr(
new MatrixDouble);
998 prism_post_proc.getOpPtrVector().push_back(
1000 prism_post_proc.getOpPtrVector().push_back(
1002 CHKERR prism_post_proc.addFieldValuesPostProc(
"DISPLACEMENT");
1003 CHKERR prism_post_proc.addFieldValuesPostProc(
"MESH_NODE_POSITIONS");
1004 CHKERR prism_post_proc.addFieldValuesGradientPostProc(
"DISPLACEMENT");
1006 m_field, prism_post_proc.postProcMesh, prism_post_proc.mapGaussPts,
1007 "DISPLACEMENT", prism_post_proc.commonData, block_sets_ptr.get()));
1016 CHKERR set_post_proc_skin(post_proc_skin);
1017 CHKERR set_post_proc_tets(post_proc);
1018 CHKERR set_post_proc_prisms(prism_post_proc);
1019 CHKERR set_post_proc_edge(post_proc_edge);
1021 PetscBool field_eval_flag = PETSC_FALSE;
1022 std::array<double, 3> field_eval_coords;
1023 boost::shared_ptr<FieldEvaluatorInterface::SetPtsData> field_eval_data;
1024 PetscInt coords_dim = 3;
1026 field_eval_coords.data(), &coords_dim,
1029 auto scalar_field_ptr = boost::make_shared<VectorDouble>();
1030 auto vector_field_ptr = boost::make_shared<MatrixDouble>();
1031 auto tensor_field_ptr = boost::make_shared<MatrixDouble>();
1033 if (field_eval_flag) {
1035 ->getData<VolumeElementForcesAndSourcesCore>();
1037 field_eval_data,
"ELASTIC");
1039 field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
1040 auto no_rule = [](int, int, int) {
return -1; };
1042 auto field_eval_fe_ptr = field_eval_data->feMethodPtr.lock();
1043 field_eval_fe_ptr->getRuleHook = no_rule;
1046 field_eval_fe_ptr->getOpPtrVector().push_back(
1049 field_eval_fe_ptr->getOpPtrVector().push_back(
1051 field_eval_fe_ptr->getOpPtrVector().push_back(
1060 CHKERR VecDuplicate(
F, &F_thermal);
1065 "DISPLACEMENT",
"TEMP", F_thermal);
1083 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Process step %d",
1084 sit->get_step_number());
1087 sit->get_step_number());
1089 CHKERR VecZeroEntries(F_thermal);
1090 CHKERR VecGhostUpdateBegin(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1091 CHKERR VecGhostUpdateEnd(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1102 CHKERR VecAssemblyBegin(F_thermal);
1103 CHKERR VecAssemblyEnd(F_thermal);
1105 CHKERR VecGhostUpdateBegin(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1106 CHKERR VecGhostUpdateEnd(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1110 CHKERR VecNorm(
F, NORM_2, &nrm_F);
1112 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"norm2 F = %6.4e", nrm_F);
1113 PetscReal nrm_F_thermal;
1114 CHKERR VecNorm(F_thermal, NORM_2, &nrm_F_thermal);
1115 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"norm2 F_thermal = %6.4e",
1118 CHKERR VecScale(F_thermal, -1);
1120 CHKERR VecAXPY(F_thermal, 1,
F);
1123 dirichlet_bc_ptr->snes_x =
D;
1124 dirichlet_bc_ptr->snes_f = F_thermal;
1128 CHKERR KSPSolve(solver, F_thermal,
D);
1131 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
1132 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
1140 if (field_eval_flag) {
1142 ->evalFEAtThePoint<3>(
1143 field_eval_coords.data(), 1e-12,
"ELASTIC_PROB",
"ELASTIC",
1146 if (scalar_field_ptr->size()) {
1149 <<
"Eval point TEMP: " << t_temp;
1151 if (vector_field_ptr->size1()) {
1155 <<
"Eval point DISPLACEMENT magnitude: "
1156 << sqrt(t_disp(
i) * t_disp(
i));
1158 if (tensor_field_ptr->size1()) {
1162 <<
"Eval point DISPLACEMENT_GRAD trace: " << t_disp_grad(
i,
i);
1169 if (is_post_proc_volume == PETSC_TRUE) {
1170 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Write output file ...";
1172 std::ostringstream o1;
1173 o1 <<
"out_" << sit->step_number <<
".h5m";
1176 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done ...";
1179 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Write output file skin ...";
1182 std::ostringstream o1_skin;
1183 o1_skin <<
"out_skin_" << sit->step_number <<
".h5m";
1186 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done ...";
1191 CHKERR VecZeroEntries(F_thermal);
1192 CHKERR VecGhostUpdateBegin(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1193 CHKERR VecGhostUpdateEnd(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1200 CHKERR VecAssemblyBegin(F_thermal);
1201 CHKERR VecAssemblyEnd(F_thermal);
1204 CHKERR VecGhostUpdateBegin(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1205 CHKERR VecGhostUpdateEnd(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1209 CHKERR VecNorm(
F, NORM_2, &nrm_F);
1211 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"norm2 F = %6.4e", nrm_F);
1212 PetscReal nrm_F_thermal;
1213 CHKERR VecNorm(F_thermal, NORM_2, &nrm_F_thermal);
1215 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"norm2 F_thermal = %6.4e",
1219 CHKERR VecScale(F_thermal, -1);
1220 CHKERR VecAXPY(F_thermal, 1,
F);
1223 dirichlet_bc_ptr->snes_x =
D;
1224 dirichlet_bc_ptr->snes_f = F_thermal;
1228 CHKERR KSPSolve(solver, F_thermal,
D);
1232 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
1233 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
1237 if (is_post_proc_volume == PETSC_TRUE) {
1238 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Write output file ...";
1243 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
1246 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Write output file skin ...";
1250 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
1254 CHKERR VecDestroy(&F_thermal);
1267 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
1268 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
1272 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Post-process start ...";
1273 if (is_post_proc_volume == PETSC_TRUE) {
1279 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
1282 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Write output file ...";
1283 if (mesh_has_tets) {
1284 if (is_post_proc_volume == PETSC_TRUE) {
1291 if (mesh_has_prisms) {
1295 if (!edges_in_simple_rod.empty())
1298 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
1301 if (is_calculating_frequency == PETSC_TRUE) {
1304 VecDuplicate(
D, &u1);
1307 CHKERR VecDot(u1,
D, &mode_mass);
1308 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Mode mass %6.4e\n", mode_mass);
1311 VecDuplicate(
D, &v1);
1314 double mode_stiffness;
1315 CHKERR VecDot(v1,
D, &mode_stiffness);
1316 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Mode stiffness %6.4e\n",
1320 double pi = 3.14159265359;
1321 frequency = std::sqrt(mode_stiffness / mode_mass) / (2 * pi);
1322 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Frequency %6.4e", frequency);
1326 auto calculate_strain_energy = [&]() {
1330 CHKERR HookeElement::calculateEnergy(dm, block_sets_ptr,
"DISPLACEMENT",
1331 "MESH_NODE_POSITIONS",
false,
true,
1336 CHKERR VecSum(v_energy, &energy);
1337 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Elastic energy %6.4e", energy);
1341 if (fabs(energy - 17.129) > 1e-3)
1343 "atom test diverged!");
1346 if (fabs(energy - 5.6475e-03) > 1e-4)
1348 "atom test diverged!");
1351 if (fabs(energy - 7.4679e-03) > 1e-4)
1353 "atom test diverged!");
1356 if (fabs(energy - 2.4992e+00) > 1e-3)
1358 "atom test diverged!");
1363 CHKERR VecMin(
D, PETSC_NULLPTR, &min);
1364 constexpr double expected_val = 0.10001;
1365 if (fabs(min + expected_val) > 1e-10)
1367 "atom test diverged! %3.4e != %3.4e", min, expected_val);
1370 if (fabs(energy - 4.7416e-04) > 1e-8)
1372 "atom test diverged!");
1380 CHKERR calculate_strain_energy();
1382 MPI_Comm_free(&moab_comm_world);