96 {
97
98 const string default_options = "-ksp_type gmres \n"
99 "-pc_type lu \n"
100 "-pc_factor_mat_solver_type mumps \n"
101 "-mat_mumps_icntl_20 0 \n"
102 "-ksp_monitor \n"
103 "-snes_type newtonls \n"
104 "-snes_linesearch_type basic \n"
105 "-snes_atol 1e-8 \n"
106 "-snes_rtol 1e-8 \n"
107 "-snes_monitor \n"
108 "-ts_monitor \n"
109 "-ts_type beuler \n";
110
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;
116 file.close();
117 }
118 }
119
121 auto core_log = logging::core::get();
122 core_log->add_sink(
126
127 core_log->add_sink(
131
132 try {
133
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;
142
143
144 enum bases { LEGENDRE, LOBATTO, BERNSTEIN_BEZIER, JACOBI, LASBASETOP };
145 const char *list_bases[] = {"legendre", "lobatto", "bernstein_bezier",
146 "jacobi"};
147 PetscInt choice_base_value = LOBATTO;
148
149
150 PetscOptionsBegin(PETSC_COMM_WORLD, "", "Elastic Config", "none");
151 CHKERR PetscOptionsString(
"-my_file",
"mesh file name",
"",
"mesh.h5m",
153
154 CHKERR PetscOptionsInt(
"-my_order",
"default approximation order",
"",
156
157 CHKERR PetscOptionsEList(
"-base",
"approximation base",
"", list_bases,
158 LASBASETOP, list_bases[choice_base_value],
159 &choice_base_value, PETSC_NULLPTR);
160
161 CHKERR PetscOptionsInt(
"-is_atom_test",
"ctest number",
"", test_nb,
162 &test_nb, PETSC_NULLPTR);
163
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);
168
169 CHKERR PetscOptionsString(
"-my_block_config",
"elastic configure file name",
170 "", "block_conf.in", block_config_file, 255,
171 &flg_block_config);
172
174 "-my_is_calculating_frequency", "set if frequency will be calculated",
175 "", is_calculating_frequency, &is_calculating_frequency, PETSC_NULLPTR);
176
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);
180
181 PetscOptionsEnd();
182
183
184 if (flg_file != PETSC_TRUE) {
185 SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
186 }
187
188
189 moab::Core mb_instance;
190 moab::Interface &moab = mb_instance;
191
192
193
194
195
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);
199 if (pcomm == NULL)
200 pcomm = new ParallelComm(&moab, moab_comm_world);
201
202
203 if (is_partitioned == PETSC_TRUE) {
204
205
206 const char *option;
207 option = "PARALLEL=READ_PART;"
208 "PARALLEL_RESOLVE_SHARED_ENTS;"
209 "PARTITION=PARALLEL_PARTITION;";
211 } else {
212
213
214
215
216 const char *option;
217 option = "";
219 }
220
221
224
225
231
232 bool mesh_has_tets = false;
233 bool mesh_has_prisms = false;
234 int nb_tets = 0;
235 int nb_hexs = 0;
236 int nb_prisms = 0;
237
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);
241
242 mesh_has_tets = (nb_tets + nb_hexs) > 0;
243 mesh_has_prisms = nb_prisms > 0;
244
245
246
248 bit_level0.set(0);
250 0, 3, bit_level0);
251
252
253
254
255
257 if (
bit->getName().compare(0, 3,
"ROD") == 0) {
259 0, 1, bit_level0);
260 }
261 }
262
263
265 switch (choice_base_value) {
266 case LEGENDRE:
268 break;
269 case LOBATTO:
271 break;
272 case BERNSTEIN_BEZIER:
274 break;
275 case JACOBI:
277 break;
278 default:
280 };
283
284
287
288
289
290
291
294
295
297 CHKERR m_field.
get_moab().get_entities_by_type(0, MBEDGE, all_edges,
true);
298
299
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);
307 }
308 }
309
311 "DISPLACEMENT");
312
313
315
316
317 Range edges_to_set_order;
318 edges_to_set_order = subtract(all_edges, edges_in_simple_rod);
319
320
321
322
330
333 else
335
336
337
338 auto setting_second_order_geometry = [&m_field]() {
340
344 moab::Interface::UNION);
345
346
347
348
349
350
351
352
353
354
355
358
360 };
361 CHKERR setting_second_order_geometry();
362
363
364
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;
375 it)) {
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));
383
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));
390
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));
397
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)
404 ->default_value(0));
405 }
406 po::parsed_options parsed =
407 parse_config_file(ini_file, config_file_options, true);
408 store(parsed, vm);
409 po::notify(vm);
411 it)) {
412 if (block_data[it->getMeshsetId()].oRder == -1)
413 continue;
414 if (block_data[it->getMeshsetId()].oRder ==
order)
415 continue;
416 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Set block %d order to %d",
417 it->getMeshsetId(),
418 block_data[it->getMeshsetId()].oRder);
420 CHKERR moab.get_entities_by_handle(it->getMeshset(), block_ents,
421 true);
422 Range ents_to_set_order;
423 CHKERR moab.get_adjacencies(block_ents, 3,
false,
424 ents_to_set_order,
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,
428 ents_to_set_order,
429 moab::Interface::UNION);
430 CHKERR moab.get_adjacencies(block_ents, 1,
false,
431 ents_to_set_order,
432 moab::Interface::UNION);
434 ents_to_set_order);
435
437 ents_to_set_order, "DISPLACEMENT",
438 block_data[it->getMeshsetId()].oRder);
439 }
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",
447 vit->c_str());
448 }
449 }
450
451
452
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",
464 bd.PoissonRatio);
465 }
466
468 };
469
470
471
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);
476
477 boost::shared_ptr<std::map<int, MassBlockData>> mass_block_sets_ptr =
478 boost::make_shared<std::map<int, MassBlockData>>();
480
481 auto fe_lhs_ptr =
482 boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
483 auto fe_rhs_ptr =
484 boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
485 fe_lhs_ptr->getRuleHook =
VolRule();
486 fe_rhs_ptr->getRuleHook =
VolRule();
487
489 false);
491 false);
492
493 boost::shared_ptr<ForcesAndSourcesCore> prism_fe_lhs_ptr(
495 boost::shared_ptr<ForcesAndSourcesCore> prism_fe_rhs_ptr(
497
498 CHKERR HookeElement::addElasticElement(m_field, block_sets_ptr,
"ELASTIC",
499 "DISPLACEMENT",
500 "MESH_NODE_POSITIONS", false);
501
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);
515 } else {
516 proc_skin = skin_faces;
517 }
520 "DISPLACEMENT");
522 "DISPLACEMENT");
524 "DISPLACEMENT");
526
527
528
529
531 "TEMP");
532 }
534 "POST_PROC_SKIN", "MESH_NODE_POSITIONS");
536 "POST_PROC_SKIN");
538 };
539 CHKERR add_skin_element_for_post_processing();
540
541 auto data_at_pts = boost::make_shared<HookeElement::DataAtIntegrationPts>();
542 if (mesh_has_tets) {
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);
546 }
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);
551 }
552
553 if (test_nb == 4) {
554
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;
564 };
565
566 fe_rhs_ptr->getOpPtrVector().push_back(
567 new HookeElement::OpAnalyticalInternalStrain_dx<0>(
568 "DISPLACEMENT", data_at_pts, thermal_strain));
569 }
570
571 boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_mass_ptr(
573
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));
580 }
581 }
582
583
584
586 "MESH_NODE_POSITIONS");
587
588
589
590 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr(
592 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr(
594
596 fe_spring_rhs_ptr, "DISPLACEMENT",
597 "MESH_NODE_POSITIONS");
598
599
600
602 "MESH_NODE_POSITIONS");
603
604
605
606
607
608
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));
611
612
614 m_field, fe_simple_rod_lhs_ptr, fe_simple_rod_rhs_ptr, "DISPLACEMENT",
615 "MESH_NODE_POSITIONS");
616
617
618
621 "DISPLACEMENT");
623 "DISPLACEMENT");
625 "DISPLACEMENT");
627 "MESH_NODE_POSITIONS");
628
632 CHKERR m_field.
get_moab().get_entities_by_dimension(it->meshset, 3, tets,
633 true);
635 }
637
638
639
643
644
645
647
648
649 fluid_pressure_fe.addNeumannFluidPressureBCElements("DISPLACEMENT");
650
651
653
655 bool add_temp_field = false;
657 if (block_data[it->getMeshsetId()].initTemp != 0) {
658 add_temp_field = true;
659 break;
660 }
661 }
662 if (add_temp_field) {
665
668 }
669 }
671 CHKERR thermal_stress_elem.addThermalStressElement(
672 "ELASTIC", "DISPLACEMENT", "TEMP");
673 }
674
675
677
678
679
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,
695 "TEMP");
696 }
697 }
698 }
699
700
701
703
705
706
708
709
710 auto dm =
createDM(PETSC_COMM_WORLD,
"MOFEM");
712 CHKERR DMSetFromOptions(dm);
714
724
725
726
733 CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
734
735
737 if (is_calculating_frequency == PETSC_TRUE) {
739 CHKERR MatSetOption(Mij, MAT_SPD, PETSC_TRUE);
740
741 }
742
743
744 fe_spring_lhs_ptr->ksp_B = Aij;
745 fe_spring_rhs_ptr->ksp_f =
F;
746
747
748 fe_simple_rod_lhs_ptr->ksp_B = Aij;
749 fe_simple_rod_rhs_ptr->ksp_f =
F;
750
751
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);
760
761
762
763
764
765
766
767
768
769
770
771 auto dirichlet_bc_ptr = boost::make_shared<DirichletDisplacementBc>(
772 m_field,
"DISPLACEMENT", Aij, D0,
F);
773
774
775
776 dirichlet_bc_ptr->snes_ctx = FEMethod::CTX_SNESNONE;
777 dirichlet_bc_ptr->ts_ctx = FEMethod::CTX_TSNONE;
778
779
780 CHKERR VecZeroEntries(D0);
781 CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
782 CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
784
785
786
788
789
790 CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
791 CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
793
794
795
796
797
798
799
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";
806
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";
813
814
817
818
821
822 if (is_calculating_frequency == PETSC_TRUE) {
823
824 fe_mass_ptr->snes_B = Mij;
825 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Calculate mass matrix ...";
827 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
828 }
829
830
831
832
833
834
835 boost::ptr_map<std::string, NeumannForcesSurface> neumann_forces;
838
839 {
840 boost::ptr_map<std::string, NeumannForcesSurface>::iterator mit =
841 neumann_forces.begin();
842 for (; mit != neumann_forces.end(); mit++) {
844 &mit->second->getLoopFe());
845 }
846 }
847
848 boost::ptr_map<std::string, NodalForce> nodal_forces;
850 "DISPLACEMENT");
851
852 {
853 boost::ptr_map<std::string, NodalForce>::iterator fit =
854 nodal_forces.begin();
855 for (; fit != nodal_forces.end(); fit++) {
857 &fit->second->getLoopFe());
858 }
859 }
860
861 boost::ptr_map<std::string, EdgeForce> edge_forces;
863 "DISPLACEMENT");
864 {
865 auto fit = edge_forces.begin();
866 for (; fit != edge_forces.end(); fit++) {
867 auto &fe = fit->second->getLoopFe();
869 }
870 }
871
875 CHKERR body_forces_methods.addBlock(
"DISPLACEMENT",
F,
876 it->getMeshsetId());
877 }
879 &body_forces_methods.getLoopFe());
880
882 false, false);
883 CHKERR fluid_pressure_fe.setNeumannFluidPressureFiniteElementOperators(
884 "DISPLACEMENT",
F,
false,
true);
885
887 &fluid_pressure_fe.getLoopFe());
888
890
891
892 PetscViewerPushFormat(
893 PETSC_VIEWER_STDOUT_SELF,
894 PETSC_VIEWER_ASCII_MATLAB);
895
896
897
898
899
900
901 if (is_calculating_frequency == PETSC_TRUE) {
902 CHKERR MatAssemblyBegin(Mij, MAT_FINAL_ASSEMBLY);
903 CHKERR MatAssemblyEnd(Mij, MAT_FINAL_ASSEMBLY);
904 }
905
906
907
908
909 CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
910 CHKERR VecGhostUpdateBegin(
F, ADD_VALUES, SCATTER_REVERSE);
911 CHKERR VecGhostUpdateEnd(
F, ADD_VALUES, SCATTER_REVERSE);
915
916
917 auto solver =
createKSP(PETSC_COMM_WORLD);
918 CHKERR KSPSetDM(solver, dm);
919 CHKERR KSPSetFromOptions(solver);
920 CHKERR KSPSetOperators(solver, Aij, Aij);
921
922
923 {
924
925 PetscBool same = PETSC_FALSE;
926 PC pc;
927 CHKERR KSPGetPC(solver, &pc);
928 PetscObjectTypeCompare((PetscObject)pc, PCMG, &same);
929 if (same) {
932 CHKERR PCSetFromOptions(pc);
933 } else {
934
935 CHKERR KSPSetDMActive(solver, PETSC_FALSE);
936 }
937 }
938 CHKERR KSPSetInitialGuessKnoll(solver, PETSC_FALSE);
939 CHKERR KSPSetInitialGuessNonzero(solver, PETSC_TRUE);
940
942
943 auto set_post_proc_skin = [&](auto &post_proc_skin) {
946 false);
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");
957 }
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,
963 true));
965 };
966
967 auto set_post_proc_tets = [&](auto &post_proc) {
969
970 CHKERR post_proc.generateReferenceElementMesh();
972 false);
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");
979 }
980
982 m_field, post_proc.postProcMesh, post_proc.mapGaussPts,
983 "DISPLACEMENT", post_proc.commonData, block_sets_ptr.get()));
985 };
986
987 auto set_post_proc_edge = [&](auto &post_proc_edge) {
989 CHKERR post_proc_edge.generateReferenceElementMesh();
990 CHKERR post_proc_edge.addFieldValuesPostProc(
"DISPLACEMENT");
992 };
993
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()));
1009 };
1010
1015
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);
1020
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,
1027 &field_eval_flag);
1028
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>();
1032
1033 if (field_eval_flag) {
1035 ->getData<VolumeElementForcesAndSourcesCore>();
1037 field_eval_data, "ELASTIC");
1038
1039 field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
1040 auto no_rule = [](int, int, int) { return -1; };
1041
1042 auto field_eval_fe_ptr = field_eval_data->feMethodPtr.lock();
1043 field_eval_fe_ptr->getRuleHook = no_rule;
1044
1046 field_eval_fe_ptr->getOpPtrVector().push_back(
1048 }
1049 field_eval_fe_ptr->getOpPtrVector().push_back(
1051 field_eval_fe_ptr->getOpPtrVector().push_back(
1053 }
1054
1055
1057
1058
1059 Vec F_thermal;
1060 CHKERR VecDuplicate(
F, &F_thermal);
1061
1062
1063
1064 CHKERR thermal_stress_elem.setThermalStressRhsOperators(
1065 "DISPLACEMENT", "TEMP", F_thermal);
1066
1069
1070
1071
1072
1073
1074
1075
1076
1078
1079
1080
1082 sit)) {
1083 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Process step %d",
1084 sit->get_step_number());
1085
1087 sit->get_step_number());
1088
1089 CHKERR VecZeroEntries(F_thermal);
1090 CHKERR VecGhostUpdateBegin(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1091 CHKERR VecGhostUpdateEnd(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1092
1093
1094
1095
1096
1097
1099 dm, "ELASTIC", &thermal_stress_elem.getLoopThermalStressRhs());
1100
1101
1102 CHKERR VecAssemblyBegin(F_thermal);
1103 CHKERR VecAssemblyEnd(F_thermal);
1104
1105 CHKERR VecGhostUpdateBegin(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1106 CHKERR VecGhostUpdateEnd(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1107
1108
1109 PetscReal nrm_F;
1110 CHKERR VecNorm(
F, NORM_2, &nrm_F);
1111
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",
1116 nrm_F_thermal);
1117
1118 CHKERR VecScale(F_thermal, -1);
1119
1120 CHKERR VecAXPY(F_thermal, 1,
F);
1121
1122
1123 dirichlet_bc_ptr->snes_x =
D;
1124 dirichlet_bc_ptr->snes_f = F_thermal;
1126
1127
1128 CHKERR KSPSolve(solver, F_thermal,
D);
1129
1131 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
1132 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
1133
1134
1136
1137
1139
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;
1150 }
1151 if (vector_field_ptr->size1()) {
1155 << "Eval point DISPLACEMENT magnitude: "
1156 << sqrt(t_disp(
i) * t_disp(
i));
1157 }
1158 if (tensor_field_ptr->size1()) {
1162 <<
"Eval point DISPLACEMENT_GRAD trace: " << t_disp_grad(
i,
i);
1163 }
1164
1166 }
1167
1168
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";
1174 if (!test_nb)
1175 CHKERR post_proc.writeFile(o1.str().c_str());
1176 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done ...";
1177 }
1178
1179 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Write output file skin ...";
1181 &post_proc_skin);
1182 std::ostringstream o1_skin;
1183 o1_skin << "out_skin_" << sit->step_number << ".h5m";
1184 if (!test_nb)
1185 CHKERR post_proc_skin.writeFile(o1_skin.str().c_str());
1186 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done ...";
1187 }
1188 } else {
1189
1190
1191 CHKERR VecZeroEntries(F_thermal);
1192 CHKERR VecGhostUpdateBegin(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1193 CHKERR VecGhostUpdateEnd(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1194
1195
1197 dm, "ELASTIC", &thermal_stress_elem.getLoopThermalStressRhs());
1198
1199
1200 CHKERR VecAssemblyBegin(F_thermal);
1201 CHKERR VecAssemblyEnd(F_thermal);
1202
1203
1204 CHKERR VecGhostUpdateBegin(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1205 CHKERR VecGhostUpdateEnd(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1206
1207
1208 PetscReal nrm_F;
1209 CHKERR VecNorm(
F, NORM_2, &nrm_F);
1210
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);
1214
1215 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"norm2 F_thermal = %6.4e",
1216 nrm_F_thermal);
1217
1218
1219 CHKERR VecScale(F_thermal, -1);
1220 CHKERR VecAXPY(F_thermal, 1,
F);
1221
1222
1223 dirichlet_bc_ptr->snes_x =
D;
1224 dirichlet_bc_ptr->snes_f = F_thermal;
1226
1227
1228 CHKERR KSPSolve(solver, F_thermal,
D);
1230
1231
1232 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
1233 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
1235
1236
1237 if (is_post_proc_volume == PETSC_TRUE) {
1238 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Write output file ...";
1240
1241 if (!test_nb)
1242 CHKERR post_proc.writeFile(
"out.h5m");
1243 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
1244 }
1245
1246 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Write output file skin ...";
1248 if (!test_nb)
1249 CHKERR post_proc_skin.writeFile(
"out_skin.h5m");
1250 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
1251 }
1252
1253
1254 CHKERR VecDestroy(&F_thermal);
1255 } else {
1256
1257
1258
1260
1261
1262
1263
1264
1266
1267 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
1268 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
1269
1271
1272 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Post-process start ...";
1273 if (is_post_proc_volume == PETSC_TRUE) {
1275 }
1279 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
1280
1281
1282 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Write output file ...";
1283 if (mesh_has_tets) {
1284 if (is_post_proc_volume == PETSC_TRUE) {
1285 if (!test_nb)
1286 CHKERR post_proc.writeFile(
"out.h5m");
1287 }
1288 if (!test_nb)
1289 CHKERR post_proc_skin.writeFile(
"out_skin.h5m");
1290 }
1291 if (mesh_has_prisms) {
1292 if (!test_nb)
1293 CHKERR prism_post_proc.writeFile(
"prism_out.h5m");
1294 }
1295 if (!edges_in_simple_rod.empty())
1296 if (!test_nb)
1297 CHKERR post_proc_edge.writeFile(
"out_edge.h5m");
1298 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
1299 }
1300
1301 if (is_calculating_frequency == PETSC_TRUE) {
1302
1303 Vec u1;
1304 VecDuplicate(
D, &u1);
1306 double mode_mass;
1307 CHKERR VecDot(u1,
D, &mode_mass);
1308 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Mode mass %6.4e\n", mode_mass);
1309
1310 Vec v1;
1311 VecDuplicate(
D, &v1);
1313
1314 double mode_stiffness;
1315 CHKERR VecDot(v1,
D, &mode_stiffness);
1316 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Mode stiffness %6.4e\n",
1317 mode_stiffness);
1318
1319 double frequency;
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);
1323 }
1324
1325
1326 auto calculate_strain_energy = [&]() {
1328
1330 CHKERR HookeElement::calculateEnergy(dm, block_sets_ptr,
"DISPLACEMENT",
1331 "MESH_NODE_POSITIONS", false, true,
1332 v_energy);
1333
1334
1335 double energy;
1336 CHKERR VecSum(v_energy, &energy);
1337 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Elastic energy %6.4e", energy);
1338
1339 switch (test_nb) {
1340 case 1:
1341 if (fabs(energy - 17.129) > 1e-3)
1343 "atom test diverged!");
1344 break;
1345 case 2:
1346 if (fabs(energy - 5.6475e-03) > 1e-4)
1348 "atom test diverged!");
1349 break;
1350 case 3:
1351 if (fabs(energy - 7.4679e-03) > 1e-4)
1353 "atom test diverged!");
1354 break;
1355 case 4:
1356 if (fabs(energy - 2.4992e+00) > 1e-3)
1358 "atom test diverged!");
1359 break;
1360
1361 case 8: {
1362 double min;
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);
1368 } break;
1369 case 9: {
1370 if (fabs(energy - 4.7416e-04) > 1e-8)
1372 "atom test diverged!");
1373 }
1374 default:
1375 break;
1376 }
1377
1379 };
1380 CHKERR calculate_strain_energy();
1381
1382 MPI_Comm_free(&moab_comm_world);
1383 }
1385
1387
1388 return 0;
1389}
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
#define MOFEM_LOG_C(channel, severity, format,...)
Kronecker Delta class symmetric.
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
#define CATCH_ERRORS
Catch errors.
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
@ 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 ...
@ 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.
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)
MoFEMErrorCode DMRegister_MGViaApproxOrders(const char sname[])
Register DM for Multi-Grid via approximation orders.
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.
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_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
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
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
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
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.
boost::shared_ptr< PCMGSetUpViaApproxOrdersCtx > createPCMGSetUpViaApproxOrdersCtx(DM dm, Mat A, bool use_shell_mat)
createPCMGSetUpViaApproxOrdersCtx
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
FTensor::Tensor1< FTensor::PackPtr< T *, S >, Tensor_Dim > getFTensor1FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 1 (vector) form data matrix.
FTensor::Tensor2< FTensor::PackPtr< double *, 1 >, Tensor_Dim1, Tensor_Dim2 > getFTensor2FromMat(MatrixDouble &data)
Get tensor rank 2 (matrix) form data matrix.
MoFEMErrorCode PCMGSetUpViaApproxOrders(PC pc, boost::shared_ptr< PCMGSetUpViaApproxOrdersCtx > ctx, int verb)
Function build MG structure.
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 setBlocks()
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.
virtual MPI_Comm & get_comm() const =0
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.
Log manager is used to build and partition problems.
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.
Get value at integration points for scalar field.
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Get values at integration pts for tensor field 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 reference to pointer of interface.
Volume finite element base.
Operator post-procesing stresses for Hook isotropic material.
Implentation of thermal stress element.
Set integration rule to volume elements.