52 using namespace MoFEM;
55 using namespace boost::numeric;
57 static 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"
87 int getRuleTrianglesOnly(
int order);
88 int getRuleThroughThickness(
int order);
96 int 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 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Elastic Config",
"none");
152 CHKERR PetscOptionsString(
"-my_file",
"mesh file name",
"",
"mesh.h5m",
155 CHKERR PetscOptionsInt(
"-my_order",
"default approximation order",
"",
158 CHKERR PetscOptionsEList(
"-base",
"approximation base",
"", list_bases,
159 LASBASETOP, list_bases[choice_base_value],
160 &choice_base_value, PETSC_NULL);
162 CHKERR PetscOptionsInt(
"-is_atom_test",
"ctest number",
"", test_nb,
163 &test_nb, PETSC_NULL);
165 CHKERR PetscOptionsBool(
"-my_is_partitioned",
166 "set if mesh is partitioned (this result that each "
167 "process keeps only one part of the mesh)",
168 "", is_partitioned, &is_partitioned, PETSC_NULL);
170 CHKERR PetscOptionsString(
"-my_block_config",
"elastic configure file name",
171 "",
"block_conf.in", block_config_file, 255,
175 "-my_is_calculating_frequency",
"set if frequency will be calculated",
176 "", is_calculating_frequency, &is_calculating_frequency, PETSC_NULL);
178 CHKERR PetscOptionsBool(
"-my_is_post_proc_volume",
179 "if true post proc volume",
"", is_post_proc_volume,
180 &is_post_proc_volume, PETSC_NULL);
182 ierr = PetscOptionsEnd();
186 if (flg_file != PETSC_TRUE) {
187 SETERRQ(PETSC_COMM_SELF, 1,
"*** ERROR -my_file (MESH FILE NEEDED)");
198 MPI_Comm moab_comm_world;
199 MPI_Comm_dup(PETSC_COMM_WORLD, &moab_comm_world);
200 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
202 pcomm =
new ParallelComm(&moab, moab_comm_world);
205 if (is_partitioned == PETSC_TRUE) {
209 option =
"PARALLEL=READ_PART;"
210 "PARALLEL_RESOLVE_SHARED_ENTS;"
211 "PARTITION=PARALLEL_PARTITION;";
234 bool mesh_has_tets =
false;
235 bool mesh_has_prisms =
false;
240 CHKERR moab.get_number_entities_by_type(0, MBTET, nb_tets,
true);
241 CHKERR moab.get_number_entities_by_type(0, MBHEX, nb_hexs,
true);
242 CHKERR moab.get_number_entities_by_type(0, MBPRISM, nb_prisms,
true);
244 mesh_has_tets = (nb_tets + nb_hexs) > 0;
245 mesh_has_prisms = nb_prisms > 0;
259 if (
bit->getName().compare(0, 3,
"ROD") == 0) {
267 switch (choice_base_value) {
274 case BERNSTEIN_BEZIER:
299 CHKERR m_field.
get_moab().get_entities_by_type(0, MBEDGE, all_edges,
true);
302 Range edges_in_simple_rod;
304 if (
bit->getName().compare(0, 3,
"ROD") == 0) {
307 MBEDGE, edges,
true);
308 edges_in_simple_rod.merge(edges);
319 Range edges_to_set_order;
320 edges_to_set_order = subtract(all_edges, edges_in_simple_rod);
340 auto setting_second_order_geometry = [&m_field]() {
346 moab::Interface::UNION);
363 CHKERR setting_second_order_geometry();
367 std::map<int, BlockOptionData> block_data;
368 auto setting_blocks_data_and_order_from_config_file =
369 [&m_field, &moab, &block_data, flg_block_config, block_config_file,
370 order](boost::shared_ptr<std::map<int, BlockData>> &block_sets_ptr) {
372 if (flg_block_config) {
373 ifstream ini_file(block_config_file);
374 po::variables_map vm;
375 po::options_description config_file_options;
378 std::ostringstream str_order;
379 str_order <<
"block_" << it->getMeshsetId()
380 <<
".displacement_order";
381 config_file_options.add_options()(
382 str_order.str().c_str(),
383 po::value<int>(&block_data[it->getMeshsetId()].oRder)
384 ->default_value(
order));
386 std::ostringstream str_cond;
387 str_cond <<
"block_" << it->getMeshsetId() <<
".young_modulus";
388 config_file_options.add_options()(
389 str_cond.str().c_str(),
390 po::value<double>(&block_data[it->getMeshsetId()].yOung)
391 ->default_value(-1));
393 std::ostringstream str_capa;
394 str_capa <<
"block_" << it->getMeshsetId() <<
".poisson_ratio";
395 config_file_options.add_options()(
396 str_capa.str().c_str(),
397 po::value<double>(&block_data[it->getMeshsetId()].pOisson)
398 ->default_value(-2));
400 std::ostringstream str_init_temp;
401 str_init_temp <<
"block_" << it->getMeshsetId()
402 <<
".initial_temperature";
403 config_file_options.add_options()(
404 str_init_temp.str().c_str(),
405 po::value<double>(&block_data[it->getMeshsetId()].initTemp)
408 po::parsed_options parsed =
409 parse_config_file(ini_file, config_file_options,
true);
414 if (block_data[it->getMeshsetId()].oRder == -1)
416 if (block_data[it->getMeshsetId()].oRder ==
order)
418 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Set block %d order to %d",
420 block_data[it->getMeshsetId()].oRder);
422 CHKERR moab.get_entities_by_handle(it->getMeshset(), block_ents,
424 Range ents_to_set_order;
425 CHKERR moab.get_adjacencies(block_ents, 3,
false,
427 moab::Interface::UNION);
428 ents_to_set_order = ents_to_set_order.subset_by_dimension(3);
429 CHKERR moab.get_adjacencies(block_ents, 2,
false,
431 moab::Interface::UNION);
432 CHKERR moab.get_adjacencies(block_ents, 1,
false,
434 moab::Interface::UNION);
439 ents_to_set_order,
"DISPLACEMENT",
440 block_data[it->getMeshsetId()].oRder);
442 std::vector<std::string> additional_parameters;
443 additional_parameters =
444 collect_unrecognized(parsed.options, po::include_positional);
445 for (std::vector<std::string>::iterator vit =
446 additional_parameters.begin();
447 vit != additional_parameters.end(); vit++) {
448 MOFEM_LOG_C(
"ELASTIC", Sev::warning,
"Unrecognized option %s",
457 const int id = it->getMeshsetId();
458 auto &bd = (*block_sets_ptr)[id];
459 if (block_data[
id].yOung > 0)
460 bd.E = block_data[id].yOung;
461 if (block_data[
id].pOisson >= -1)
462 bd.PoissonRatio = block_data[id].pOisson;
463 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Block %d",
id);
464 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"\tYoung modulus %3.4g", bd.E);
465 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"\tPoisson ratio %3.4g",
474 boost::shared_ptr<std::map<int, HookeElement::BlockData>> block_sets_ptr =
475 boost::make_shared<std::map<int, HookeElement::BlockData>>();
477 CHKERR setting_blocks_data_and_order_from_config_file(block_sets_ptr);
479 boost::shared_ptr<std::map<int, MassBlockData>> mass_block_sets_ptr =
480 boost::make_shared<std::map<int, MassBlockData>>();
484 boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
486 boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
487 fe_lhs_ptr->getRuleHook =
VolRule();
488 fe_rhs_ptr->getRuleHook =
VolRule();
495 boost::shared_ptr<ForcesAndSourcesCore> prism_fe_lhs_ptr(
497 boost::shared_ptr<ForcesAndSourcesCore> prism_fe_rhs_ptr(
502 "MESH_NODE_POSITIONS",
false);
504 auto add_skin_element_for_post_processing = [&]() {
506 Range elastic_element_ents;
508 "ELASTIC", 3, elastic_element_ents);
511 CHKERR skin.find_skin(0, elastic_element_ents,
false, skin_faces);
513 if (is_partitioned) {
514 CHKERR pcomm->filter_pstatus(skin_faces,
515 PSTATUS_SHARED | PSTATUS_MULTISHARED,
516 PSTATUS_NOT, -1, &proc_skin);
518 proc_skin = skin_faces;
536 "POST_PROC_SKIN",
"MESH_NODE_POSITIONS");
541 CHKERR add_skin_element_for_post_processing();
543 auto data_at_pts = boost::make_shared<HookeElement::DataAtIntegrationPts>();
546 "DISPLACEMENT",
"MESH_NODE_POSITIONS",
547 false,
true, MBTET, data_at_pts);
549 if (mesh_has_prisms) {
551 prism_fe_lhs_ptr, prism_fe_rhs_ptr, block_sets_ptr,
"DISPLACEMENT",
552 "MESH_NODE_POSITIONS",
false,
true, MBPRISM, data_at_pts);
557 auto thermal_strain =
560 constexpr
double alpha = 1;
564 t_thermal_strain(
i,
j) = alpha * t_coords(2) *
t_kd(
i,
j);
565 return t_thermal_strain;
568 fe_rhs_ptr->getOpPtrVector().push_back(
569 new HookeElement::OpAnalyticalInternalStrain_dx<0>(
570 "DISPLACEMENT", data_at_pts, thermal_strain));
573 boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_mass_ptr(
576 for (
auto &sit : *block_sets_ptr) {
577 for (
auto &mit : *mass_block_sets_ptr) {
578 fe_mass_ptr->getOpPtrVector().push_back(
579 new HookeElement::OpCalculateMassMatrix(
"DISPLACEMENT",
580 "DISPLACEMENT", sit.second,
581 mit.second, data_at_pts));
588 "MESH_NODE_POSITIONS");
592 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr(
594 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr(
598 fe_spring_rhs_ptr,
"DISPLACEMENT",
599 "MESH_NODE_POSITIONS");
604 "MESH_NODE_POSITIONS");
611 boost::shared_ptr<EdgeEle> fe_simple_rod_lhs_ptr(
new EdgeEle(m_field));
612 boost::shared_ptr<EdgeEle> fe_simple_rod_rhs_ptr(
new EdgeEle(m_field));
616 m_field, fe_simple_rod_lhs_ptr, fe_simple_rod_rhs_ptr,
"DISPLACEMENT",
617 "MESH_NODE_POSITIONS");
629 "MESH_NODE_POSITIONS");
634 CHKERR m_field.
get_moab().get_entities_by_dimension(it->meshset, 3, tets,
657 bool add_temp_field =
false;
659 if (block_data[it->getMeshsetId()].initTemp != 0) {
660 add_temp_field =
true;
664 if (add_temp_field) {
674 "ELASTIC",
"DISPLACEMENT",
"TEMP");
683 "MESH_NODE_POSITIONS");
687 if (block_data[it->getMeshsetId()].initTemp != 0) {
689 "Set block %d temperature to %3.2g\n", it->getMeshsetId(),
690 block_data[it->getMeshsetId()].initTemp);
692 CHKERR moab.get_entities_by_handle(it->meshset, block_ents,
true);
694 CHKERR moab.get_connectivity(block_ents, vertices,
true);
696 block_data[it->getMeshsetId()].initTemp, MBVERTEX, vertices,
712 auto dm =
createDM(PETSC_COMM_WORLD,
"MOFEM");
714 CHKERR DMSetFromOptions(dm);
735 CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
739 if (is_calculating_frequency == PETSC_TRUE) {
741 CHKERR MatSetOption(Mij, MAT_SPD, PETSC_TRUE);
746 fe_spring_lhs_ptr->ksp_B = Aij;
747 fe_spring_rhs_ptr->ksp_f =
F;
750 fe_simple_rod_lhs_ptr->ksp_B = Aij;
751 fe_simple_rod_rhs_ptr->ksp_f =
F;
755 CHKERR VecGhostUpdateBegin(
F, INSERT_VALUES, SCATTER_FORWARD);
756 CHKERR VecGhostUpdateEnd(
F, INSERT_VALUES, SCATTER_FORWARD);
758 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
759 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
761 CHKERR MatZeroEntries(Aij);
773 auto dirichlet_bc_ptr = boost::make_shared<DirichletDisplacementBc>(
774 m_field,
"DISPLACEMENT", Aij, D0,
F);
782 CHKERR VecZeroEntries(D0);
783 CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
784 CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
792 CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
793 CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
802 fe_rhs_ptr->snes_f =
F;
803 prism_fe_rhs_ptr->snes_f =
F;
804 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Assemble external force vector ...";
807 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
809 fe_lhs_ptr->snes_B = Aij;
810 prism_fe_lhs_ptr->snes_B = Aij;
811 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Calculate stiffness matrix ...";
814 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
824 if (is_calculating_frequency == PETSC_TRUE) {
826 fe_mass_ptr->snes_B = Mij;
827 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Calculate mass matrix ...";
829 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
837 boost::ptr_map<std::string, NeumannForcesSurface> neumann_forces;
842 boost::ptr_map<std::string, NeumannForcesSurface>::iterator mit =
843 neumann_forces.begin();
844 for (; mit != neumann_forces.end(); mit++) {
846 &mit->second->getLoopFe());
850 boost::ptr_map<std::string, NodalForce> nodal_forces;
855 boost::ptr_map<std::string, NodalForce>::iterator fit =
856 nodal_forces.begin();
857 for (; fit != nodal_forces.end(); fit++) {
859 &fit->second->getLoopFe());
863 boost::ptr_map<std::string, EdgeForce> edge_forces;
867 auto fit = edge_forces.begin();
868 for (; fit != edge_forces.end(); fit++) {
869 auto &fe = fit->second->getLoopFe();
886 "DISPLACEMENT",
F,
false,
true);
894 PetscViewerPushFormat(
895 PETSC_VIEWER_STDOUT_SELF,
896 PETSC_VIEWER_ASCII_MATLAB);
903 if (is_calculating_frequency == PETSC_TRUE) {
904 CHKERR MatAssemblyBegin(Mij, MAT_FINAL_ASSEMBLY);
905 CHKERR MatAssemblyEnd(Mij, MAT_FINAL_ASSEMBLY);
911 CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
912 CHKERR VecGhostUpdateBegin(
F, ADD_VALUES, SCATTER_REVERSE);
913 CHKERR VecGhostUpdateEnd(
F, ADD_VALUES, SCATTER_REVERSE);
919 auto solver =
createKSP(PETSC_COMM_WORLD);
920 CHKERR KSPSetDM(solver, dm);
921 CHKERR KSPSetFromOptions(solver);
922 CHKERR KSPSetOperators(solver, Aij, Aij);
927 PetscBool same = PETSC_FALSE;
929 CHKERR KSPGetPC(solver, &pc);
930 PetscObjectTypeCompare((PetscObject)pc, PCMG, &same);
934 CHKERR PCSetFromOptions(pc);
937 CHKERR KSPSetDMActive(solver, PETSC_FALSE);
940 CHKERR KSPSetInitialGuessKnoll(solver, PETSC_FALSE);
941 CHKERR KSPSetInitialGuessNonzero(solver, PETSC_TRUE);
945 auto set_post_proc_skin = [&](
auto &post_proc_skin) {
949 CHKERR post_proc_skin.generateReferenceElementMesh();
950 CHKERR post_proc_skin.addFieldValuesPostProc(
"DISPLACEMENT");
951 CHKERR post_proc_skin.addFieldValuesPostProc(
"MESH_NODE_POSITIONS");
952 CHKERR post_proc_skin.addFieldValuesGradientPostProcOnSkin(
953 "DISPLACEMENT",
"ELASTIC", data_at_pts->hMat,
true);
954 CHKERR post_proc_skin.addFieldValuesGradientPostProcOnSkin(
955 "MESH_NODE_POSITIONS",
"ELASTIC", data_at_pts->HMat,
false);
957 CHKERR post_proc_skin.addFieldValuesPostProc(
"TEMP");
958 CHKERR post_proc_skin.addFieldValuesGradientPostProc(
"TEMP");
960 post_proc_skin.getOpPtrVector().push_back(
961 new HookeElement::OpPostProcHookeElement<
963 "DISPLACEMENT", data_at_pts, *block_sets_ptr,
964 post_proc_skin.postProcMesh, post_proc_skin.mapGaussPts,
true,
969 auto set_post_proc_tets = [&](
auto &post_proc) {
972 CHKERR post_proc.generateReferenceElementMesh();
975 CHKERR post_proc.addFieldValuesPostProc(
"DISPLACEMENT");
976 CHKERR post_proc.addFieldValuesPostProc(
"MESH_NODE_POSITIONS");
977 CHKERR post_proc.addFieldValuesGradientPostProc(
"DISPLACEMENT");
979 CHKERR post_proc.addFieldValuesPostProc(
"TEMP");
980 CHKERR post_proc.addFieldValuesGradientPostProc(
"TEMP");
984 m_field, post_proc.postProcMesh, post_proc.mapGaussPts,
985 "DISPLACEMENT", post_proc.commonData, block_sets_ptr.get()));
989 auto set_post_proc_edge = [&](
auto &post_proc_edge) {
991 CHKERR post_proc_edge.generateReferenceElementMesh();
992 CHKERR post_proc_edge.addFieldValuesPostProc(
"DISPLACEMENT");
996 auto set_post_proc_prisms = [&](
auto &prism_post_proc) {
998 CHKERR prism_post_proc.generateReferenceElementMesh();
999 boost::shared_ptr<MatrixDouble> inv_jac_ptr(
new MatrixDouble);
1000 prism_post_proc.getOpPtrVector().push_back(
1002 prism_post_proc.getOpPtrVector().push_back(
1004 CHKERR prism_post_proc.addFieldValuesPostProc(
"DISPLACEMENT");
1005 CHKERR prism_post_proc.addFieldValuesPostProc(
"MESH_NODE_POSITIONS");
1006 CHKERR prism_post_proc.addFieldValuesGradientPostProc(
"DISPLACEMENT");
1008 m_field, prism_post_proc.postProcMesh, prism_post_proc.mapGaussPts,
1009 "DISPLACEMENT", prism_post_proc.commonData, block_sets_ptr.get()));
1018 CHKERR set_post_proc_skin(post_proc_skin);
1019 CHKERR set_post_proc_tets(post_proc);
1020 CHKERR set_post_proc_prisms(prism_post_proc);
1021 CHKERR set_post_proc_edge(post_proc_edge);
1023 PetscBool field_eval_flag = PETSC_FALSE;
1024 std::array<double, 3> field_eval_coords;
1025 boost::shared_ptr<FieldEvaluatorInterface::SetPtsData> field_eval_data;
1026 PetscInt coords_dim = 3;
1028 field_eval_coords.data(), &coords_dim,
1031 auto scalar_field_ptr = boost::make_shared<VectorDouble>();
1032 auto vector_field_ptr = boost::make_shared<MatrixDouble>();
1033 auto tensor_field_ptr = boost::make_shared<MatrixDouble>();
1035 if (field_eval_flag) {
1037 ->getData<VolumeElementForcesAndSourcesCore>();
1039 field_eval_data,
"ELASTIC");
1041 field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
1042 auto no_rule = [](
int,
int,
int) {
return -1; };
1044 auto field_eval_fe_ptr = field_eval_data->feMethodPtr.lock();
1045 field_eval_fe_ptr->getRuleHook = no_rule;
1048 field_eval_fe_ptr->getOpPtrVector().push_back(
1051 field_eval_fe_ptr->getOpPtrVector().push_back(
1053 field_eval_fe_ptr->getOpPtrVector().push_back(
1062 CHKERR VecDuplicate(
F, &F_thermal);
1067 "DISPLACEMENT",
"TEMP", F_thermal);
1085 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Process step %d",
1086 sit->get_step_number());
1089 sit->get_step_number());
1091 CHKERR VecZeroEntries(F_thermal);
1092 CHKERR VecGhostUpdateBegin(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1093 CHKERR VecGhostUpdateEnd(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1104 CHKERR VecAssemblyBegin(F_thermal);
1105 CHKERR VecAssemblyEnd(F_thermal);
1107 CHKERR VecGhostUpdateBegin(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1108 CHKERR VecGhostUpdateEnd(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1112 CHKERR VecNorm(
F, NORM_2, &nrm_F);
1114 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"norm2 F = %6.4e", nrm_F);
1115 PetscReal nrm_F_thermal;
1116 CHKERR VecNorm(F_thermal, NORM_2, &nrm_F_thermal);
1117 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"norm2 F_thermal = %6.4e",
1120 CHKERR VecScale(F_thermal, -1);
1122 CHKERR VecAXPY(F_thermal, 1,
F);
1125 dirichlet_bc_ptr->snes_x =
D;
1126 dirichlet_bc_ptr->snes_f = F_thermal;
1130 CHKERR KSPSolve(solver, F_thermal,
D);
1133 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
1134 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
1142 if (field_eval_flag) {
1144 ->evalFEAtThePoint3D(
1145 field_eval_coords.data(), 1e-12,
"ELASTIC_PROB",
"ELASTIC",
1148 if (scalar_field_ptr->size()) {
1151 <<
"Eval point TEMP: " << t_temp;
1153 if (vector_field_ptr->size1()) {
1155 auto t_disp = getFTensor1FromMat<3>(*vector_field_ptr);
1157 <<
"Eval point DISPLACEMENT magnitude: "
1158 << sqrt(t_disp(
i) * t_disp(
i));
1160 if (tensor_field_ptr->size1()) {
1162 auto t_disp_grad = getFTensor2FromMat<3, 3>(*tensor_field_ptr);
1164 <<
"Eval point DISPLACEMENT_GRAD trace: " << t_disp_grad(
i,
i);
1171 if (is_post_proc_volume == PETSC_TRUE) {
1172 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Write output file ...";
1174 std::ostringstream o1;
1175 o1 <<
"out_" << sit->step_number <<
".h5m";
1178 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done ...";
1181 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Write output file skin ...";
1184 std::ostringstream o1_skin;
1185 o1_skin <<
"out_skin_" << sit->step_number <<
".h5m";
1188 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done ...";
1193 CHKERR VecZeroEntries(F_thermal);
1194 CHKERR VecGhostUpdateBegin(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1195 CHKERR VecGhostUpdateEnd(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1202 CHKERR VecAssemblyBegin(F_thermal);
1203 CHKERR VecAssemblyEnd(F_thermal);
1206 CHKERR VecGhostUpdateBegin(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1207 CHKERR VecGhostUpdateEnd(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1211 CHKERR VecNorm(
F, NORM_2, &nrm_F);
1213 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"norm2 F = %6.4e", nrm_F);
1214 PetscReal nrm_F_thermal;
1215 CHKERR VecNorm(F_thermal, NORM_2, &nrm_F_thermal);
1217 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"norm2 F_thermal = %6.4e",
1221 CHKERR VecScale(F_thermal, -1);
1222 CHKERR VecAXPY(F_thermal, 1,
F);
1225 dirichlet_bc_ptr->snes_x =
D;
1226 dirichlet_bc_ptr->snes_f = F_thermal;
1230 CHKERR KSPSolve(solver, F_thermal,
D);
1234 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
1235 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
1239 if (is_post_proc_volume == PETSC_TRUE) {
1240 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Write output file ...";
1245 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
1248 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Write output file skin ...";
1252 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
1256 CHKERR VecDestroy(&F_thermal);
1269 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
1270 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
1274 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Post-process start ...";
1275 if (is_post_proc_volume == PETSC_TRUE) {
1281 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
1284 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Write output file ...";
1285 if (mesh_has_tets) {
1286 if (is_post_proc_volume == PETSC_TRUE) {
1293 if (mesh_has_prisms) {
1297 if (!edges_in_simple_rod.empty())
1300 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
1303 if (is_calculating_frequency == PETSC_TRUE) {
1306 VecDuplicate(
D, &u1);
1309 CHKERR VecDot(u1,
D, &mode_mass);
1310 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Mode mass %6.4e\n", mode_mass);
1313 VecDuplicate(
D, &v1);
1316 double mode_stiffness;
1317 CHKERR VecDot(v1,
D, &mode_stiffness);
1318 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Mode stiffness %6.4e\n",
1322 double pi = 3.14159265359;
1323 frequency = std::sqrt(mode_stiffness / mode_mass) / (2 *
pi);
1324 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Frequency %6.4e", frequency);
1328 auto calculate_strain_energy = [&]() {
1333 "MESH_NODE_POSITIONS",
false,
true,
1338 CHKERR VecSum(v_energy, &energy);
1339 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Elastic energy %6.4e", energy);
1343 if (fabs(energy - 17.129) > 1e-3)
1345 "atom test diverged!");
1348 if (fabs(energy - 5.6475e-03) > 1e-4)
1350 "atom test diverged!");
1353 if (fabs(energy - 7.4679e-03) > 1e-4)
1355 "atom test diverged!");
1358 if (fabs(energy - 2.4992e+00) > 1e-3)
1360 "atom test diverged!");
1365 CHKERR VecMin(
D, PETSC_NULL, &min);
1366 constexpr
double expected_val = 0.10001;
1367 if (fabs(min + expected_val) > 1e-10)
1369 "atom test diverged! %3.4e != %3.4e", min, expected_val);
1372 if (fabs(energy - 4.7416e-04) > 1e-8)
1374 "atom test diverged!");
1382 CHKERR calculate_strain_energy();
1384 MPI_Comm_free(&moab_comm_world);