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();
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)");
191 moab::Core mb_instance;
192 moab::Interface &moab = mb_instance;
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>>();
476 CHKERR HookeElement::setBlocks(m_field, block_sets_ptr);
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(
500 CHKERR HookeElement::addElasticElement(m_field, block_sets_ptr,
"ELASTIC",
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>();
545 CHKERR HookeElement::setOperators(fe_lhs_ptr, fe_rhs_ptr, block_sets_ptr,
546 "DISPLACEMENT",
"MESH_NODE_POSITIONS",
547 false,
true, MBTET, data_at_pts);
549 if (mesh_has_prisms) {
550 CHKERR HookeElement::setOperators(
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);
778 dirichlet_bc_ptr->snes_ctx = FEMethod::CTX_SNESNONE;
779 dirichlet_bc_ptr->ts_ctx = FEMethod::CTX_TSNONE;
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 = [&]() {
1332 CHKERR HookeElement::calculateEnergy(dm, block_sets_ptr,
"DISPLACEMENT",
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);
const std::string default_options
Implementation of linear elastic material.
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
#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.
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data 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 MPI_Comm & get_comm() const =0
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
virtual int get_comm_rank() const =0
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
Field evaluator interface.
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.
static boost::shared_ptr< std::ostream > getStrmSync()
Get the strm sync 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.
Get value at integration points for scalar field.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
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