96 {
97
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
112 if (!
static_cast<bool>(ifstream(
param_file))) {
113 std::ofstream file(
param_file.c_str(), std::ios::ate);
114 if (file.is_open()) {
116 file.close();
117 }
118 }
119
121 auto core_log = logging::core::get();
122 core_log->add_sink(
126
127 try {
128
129 PetscBool flg_block_config, flg_file;
131 char block_config_file[255];
132 PetscInt test_nb = 0;
134 PetscBool is_partitioned = PETSC_FALSE;
135 PetscBool is_calculating_frequency = PETSC_FALSE;
136 PetscBool is_post_proc_volume = PETSC_TRUE;
137
138
139 enum bases { LEGENDRE, LOBATTO, BERNSTEIN_BEZIER, JACOBI, LASBASETOP };
140 const char *list_bases[] = {"legendre", "lobatto", "bernstein_bezier",
141 "jacobi"};
142 PetscInt choice_base_value = LOBATTO;
143
144
145 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Elastic Config",
"none");
147 CHKERR PetscOptionsString(
"-my_file",
"mesh file name",
"",
"mesh.h5m",
149
150 CHKERR PetscOptionsInt(
"-my_order",
"default approximation order",
"",
152
153 CHKERR PetscOptionsEList(
"-base",
"approximation base",
"", list_bases,
154 LASBASETOP, list_bases[choice_base_value],
155 &choice_base_value, PETSC_NULL);
156
157 CHKERR PetscOptionsInt(
"-is_atom_test",
"ctest number",
"", test_nb,
158 &test_nb, PETSC_NULL);
159
160 CHKERR PetscOptionsBool(
"-my_is_partitioned",
161 "set if mesh is partitioned (this result that each "
162 "process keeps only one part of the mesh)",
163 "", is_partitioned, &is_partitioned, PETSC_NULL);
164
165 CHKERR PetscOptionsString(
"-my_block_config",
"elastic configure file name",
166 "", "block_conf.in", block_config_file, 255,
167 &flg_block_config);
168
170 "-my_is_calculating_frequency", "set if frequency will be calculated",
171 "", is_calculating_frequency, &is_calculating_frequency, PETSC_NULL);
172
173 CHKERR PetscOptionsBool(
"-my_is_post_proc_volume",
174 "if true post proc volume", "", is_post_proc_volume,
175 &is_post_proc_volume, PETSC_NULL);
176
177 ierr = PetscOptionsEnd();
179
180
181 if (flg_file != PETSC_TRUE) {
182 SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
183 }
184
185
186 moab::Core mb_instance;
187 moab::Interface &moab = mb_instance;
188
189
190
191
192
193 MPI_Comm moab_comm_world;
194 MPI_Comm_dup(PETSC_COMM_WORLD, &moab_comm_world);
195 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
196 if (pcomm == NULL)
197 pcomm = new ParallelComm(&moab, moab_comm_world);
198
199
200 if (is_partitioned == PETSC_TRUE) {
201
202
203 const char *option;
204 option = "PARALLEL=READ_PART;"
205 "PARALLEL_RESOLVE_SHARED_ENTS;"
206 "PARTITION=PARALLEL_PARTITION;";
208 } else {
209
210
211
212
213 const char *option;
214 option = "";
216 }
217
218
221
222
228
229 bool mesh_has_tets = false;
230 bool mesh_has_prisms = false;
231 int nb_tets = 0;
232 int nb_hexs = 0;
233 int nb_prisms = 0;
234
235 CHKERR moab.get_number_entities_by_type(0, MBTET, nb_tets,
true);
236 CHKERR moab.get_number_entities_by_type(0, MBHEX, nb_hexs,
true);
237 CHKERR moab.get_number_entities_by_type(0, MBPRISM, nb_prisms,
true);
238
239 mesh_has_tets = (nb_tets + nb_hexs) > 0;
240 mesh_has_prisms = nb_prisms > 0;
241
242
243
245 bit_level0.set(0);
247 0, 3, bit_level0);
248
249
250
251
252
254 if (
bit->getName().compare(0, 3,
"ROD") == 0) {
256 0, 1, bit_level0);
257 }
258 }
259
260
262 switch (choice_base_value) {
263 case LEGENDRE:
265 break;
266 case LOBATTO:
268 break;
269 case BERNSTEIN_BEZIER:
271 break;
272 case JACOBI:
274 break;
275 default:
277 };
280
281
284
285
286
287
288
291
292
294 CHKERR m_field.
get_moab().get_entities_by_type(0, MBEDGE, all_edges,
true);
295
296
297 Range edges_in_simple_rod;
299 if (
bit->getName().compare(0, 3,
"ROD") == 0) {
302 MBEDGE, edges, true);
303 edges_in_simple_rod.merge(edges);
304 }
305 }
306
308 "DISPLACEMENT");
309
310
312
313
314 Range edges_to_set_order;
315 edges_to_set_order = subtract(all_edges, edges_in_simple_rod);
316
317
318
319
327
330 else
332
333
334
335 auto setting_second_order_geometry = [&m_field]() {
337
341 moab::Interface::UNION);
342
343
344
345
346
347
348
349
350
351
352
355
357 };
358 CHKERR setting_second_order_geometry();
359
360
361
362 std::map<int, BlockOptionData> block_data;
363 auto setting_blocks_data_and_order_from_config_file =
364 [&m_field, &moab, &block_data, flg_block_config, block_config_file,
365 order](boost::shared_ptr<std::map<int, BlockData>> &block_sets_ptr) {
367 if (flg_block_config) {
368 ifstream ini_file(block_config_file);
369 po::variables_map vm;
370 po::options_description config_file_options;
372 it)) {
373 std::ostringstream str_order;
374 str_order << "block_" << it->getMeshsetId()
375 << ".displacement_order";
376 config_file_options.add_options()(
377 str_order.str().c_str(),
378 po::value<int>(&block_data[it->getMeshsetId()].oRder)
379 ->default_value(
order));
380
381 std::ostringstream str_cond;
382 str_cond << "block_" << it->getMeshsetId() << ".young_modulus";
383 config_file_options.add_options()(
384 str_cond.str().c_str(),
385 po::value<double>(&block_data[it->getMeshsetId()].yOung)
386 ->default_value(-1));
387
388 std::ostringstream str_capa;
389 str_capa << "block_" << it->getMeshsetId() << ".poisson_ratio";
390 config_file_options.add_options()(
391 str_capa.str().c_str(),
392 po::value<double>(&block_data[it->getMeshsetId()].pOisson)
393 ->default_value(-2));
394
395 std::ostringstream str_init_temp;
396 str_init_temp << "block_" << it->getMeshsetId()
397 << ".initial_temperature";
398 config_file_options.add_options()(
399 str_init_temp.str().c_str(),
400 po::value<double>(&block_data[it->getMeshsetId()].initTemp)
401 ->default_value(0));
402 }
403 po::parsed_options parsed =
404 parse_config_file(ini_file, config_file_options, true);
405 store(parsed, vm);
406 po::notify(vm);
408 it)) {
409 if (block_data[it->getMeshsetId()].oRder == -1)
410 continue;
411 if (block_data[it->getMeshsetId()].oRder ==
order)
412 continue;
413 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Set block %d order to %d",
414 it->getMeshsetId(),
415 block_data[it->getMeshsetId()].oRder);
417 CHKERR moab.get_entities_by_handle(it->getMeshset(), block_ents,
418 true);
419 Range ents_to_set_order;
420 CHKERR moab.get_adjacencies(block_ents, 3,
false,
421 ents_to_set_order,
422 moab::Interface::UNION);
423 ents_to_set_order = ents_to_set_order.subset_by_dimension(3);
424 CHKERR moab.get_adjacencies(block_ents, 2,
false,
425 ents_to_set_order,
426 moab::Interface::UNION);
427 CHKERR moab.get_adjacencies(block_ents, 1,
false,
428 ents_to_set_order,
429 moab::Interface::UNION);
431 ents_to_set_order);
432
434 ents_to_set_order, "DISPLACEMENT",
435 block_data[it->getMeshsetId()].oRder);
436 }
437 std::vector<std::string> additional_parameters;
438 additional_parameters =
439 collect_unrecognized(parsed.options, po::include_positional);
440 for (std::vector<std::string>::iterator vit =
441 additional_parameters.begin();
442 vit != additional_parameters.end(); vit++) {
443 MOFEM_LOG_C(
"ELASTIC", Sev::warning,
"Unrecognized option %s",
444 vit->c_str());
445 }
446 }
447
448
449
452 const int id = it->getMeshsetId();
453 auto &bd = (*block_sets_ptr)[id];
454 if (block_data[id].yOung > 0)
455 bd.E = block_data[id].yOung;
456 if (block_data[id].pOisson >= -1)
457 bd.PoissonRatio = block_data[id].pOisson;
458 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Block %d",
id);
459 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"\tYoung modulus %3.4g", bd.E);
460 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"\tPoisson ratio %3.4g",
461 bd.PoissonRatio);
462 }
463
465 };
466
467
468
469 boost::shared_ptr<std::map<int, HookeElement::BlockData>> block_sets_ptr =
470 boost::make_shared<std::map<int, HookeElement::BlockData>>();
471 CHKERR HookeElement::setBlocks(m_field, block_sets_ptr);
472 CHKERR setting_blocks_data_and_order_from_config_file(block_sets_ptr);
473
474 boost::shared_ptr<std::map<int, MassBlockData>> mass_block_sets_ptr =
475 boost::make_shared<std::map<int, MassBlockData>>();
477
478 auto fe_lhs_ptr =
479 boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
480 auto fe_rhs_ptr =
481 boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
482 fe_lhs_ptr->getRuleHook =
VolRule();
483 fe_rhs_ptr->getRuleHook =
VolRule();
484
486 false);
488 false);
489
490 boost::shared_ptr<ForcesAndSourcesCore> prism_fe_lhs_ptr(
492 boost::shared_ptr<ForcesAndSourcesCore> prism_fe_rhs_ptr(
494
495 CHKERR HookeElement::addElasticElement(m_field, block_sets_ptr,
"ELASTIC",
496 "DISPLACEMENT",
497 "MESH_NODE_POSITIONS", false);
498
499 auto add_skin_element_for_post_processing = [&]() {
501 Range elastic_element_ents;
503 "ELASTIC", 3, elastic_element_ents);
506 CHKERR skin.find_skin(0, elastic_element_ents,
false, skin_faces);
508 if (is_partitioned) {
509 CHKERR pcomm->filter_pstatus(skin_faces,
510 PSTATUS_SHARED | PSTATUS_MULTISHARED,
511 PSTATUS_NOT, -1, &proc_skin);
512 } else {
513 proc_skin = skin_faces;
514 }
517 "DISPLACEMENT");
519 "DISPLACEMENT");
521 "DISPLACEMENT");
523 "POST_PROC_SKIN", "MESH_NODE_POSITIONS");
525 "POST_PROC_SKIN");
527 };
528 CHKERR add_skin_element_for_post_processing();
529
530 auto data_at_pts = boost::make_shared<HookeElement::DataAtIntegrationPts>();
531 if (mesh_has_tets) {
532 CHKERR HookeElement::setOperators(fe_lhs_ptr, fe_rhs_ptr, block_sets_ptr,
533 "DISPLACEMENT", "MESH_NODE_POSITIONS",
534 false, true, MBTET, data_at_pts);
535 }
536 if (mesh_has_prisms) {
537 CHKERR HookeElement::setOperators(
538 prism_fe_lhs_ptr, prism_fe_rhs_ptr, block_sets_ptr, "DISPLACEMENT",
539 "MESH_NODE_POSITIONS", false, true, MBPRISM, data_at_pts);
540 }
541
542 if (test_nb == 4) {
543
544 auto thermal_strain =
547 constexpr double alpha = 1;
552 return t_thermal_strain;
553 };
554
555 fe_rhs_ptr->getOpPtrVector().push_back(
556 new HookeElement::OpAnalyticalInternalStrain_dx<0>(
557 "DISPLACEMENT", data_at_pts, thermal_strain));
558 }
559
560 boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_mass_ptr(
562
563 for (auto &sit : *block_sets_ptr) {
564 for (auto &mit : *mass_block_sets_ptr) {
565 fe_mass_ptr->getOpPtrVector().push_back(
566 new HookeElement::OpCalculateMassMatrix("DISPLACEMENT",
567 "DISPLACEMENT", sit.second,
568 mit.second, data_at_pts));
569 }
570 }
571
572
573
575 "MESH_NODE_POSITIONS");
576
577
578
579 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr(
581 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr(
583
585 fe_spring_rhs_ptr, "DISPLACEMENT",
586 "MESH_NODE_POSITIONS");
587
588
589
591 "MESH_NODE_POSITIONS");
592
593
594
595
596
597
598 boost::shared_ptr<EdgeEle> fe_simple_rod_lhs_ptr(
new EdgeEle(m_field));
599 boost::shared_ptr<EdgeEle> fe_simple_rod_rhs_ptr(
new EdgeEle(m_field));
600
601
603 m_field, fe_simple_rod_lhs_ptr, fe_simple_rod_rhs_ptr, "DISPLACEMENT",
604 "MESH_NODE_POSITIONS");
605
606
607
610 "DISPLACEMENT");
612 "DISPLACEMENT");
614 "DISPLACEMENT");
616 "MESH_NODE_POSITIONS");
617
621 CHKERR m_field.
get_moab().get_entities_by_dimension(it->meshset, 3, tets,
622 true);
624 }
626
627
628
632
633
634
636
637
638 fluid_pressure_fe.addNeumannFluidPressureBCElements("DISPLACEMENT");
639
640
642
644 bool add_temp_field = false;
646 if (block_data[it->getMeshsetId()].initTemp != 0) {
647 add_temp_field = true;
648 break;
649 }
650 }
651 if (add_temp_field) {
654
657 }
658 }
660 CHKERR thermal_stress_elem.addThermalStressElement(
661 "ELASTIC", "DISPLACEMENT", "TEMP");
662 }
663
664
666
667
668
670 "MESH_NODE_POSITIONS");
674 if (block_data[it->getMeshsetId()].initTemp != 0) {
676 "Set block %d temperature to %3.2g\n", it->getMeshsetId(),
677 block_data[it->getMeshsetId()].initTemp);
679 CHKERR moab.get_entities_by_handle(it->meshset, block_ents,
true);
681 CHKERR moab.get_connectivity(block_ents, vertices,
true);
683 block_data[it->getMeshsetId()].initTemp, MBVERTEX, vertices,
684 "TEMP");
685 }
686 }
687 }
688
689
690
692
694
695
697
698
701 CHKERR DMSetFromOptions(dm);
703
713
714
715
722 CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
723
724
726 if (is_calculating_frequency == PETSC_TRUE) {
728 CHKERR MatSetOption(Mij, MAT_SPD, PETSC_TRUE);
729
730 }
731
732
733 fe_spring_lhs_ptr->ksp_B = Aij;
734 fe_spring_rhs_ptr->ksp_f = F;
735
736
737 fe_simple_rod_lhs_ptr->ksp_B = Aij;
738 fe_simple_rod_rhs_ptr->ksp_f = F;
739
740
742 CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
743 CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
745 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
746 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
748 CHKERR MatZeroEntries(Aij);
749
750
751
752
753
754
755
756
757
758
759
760 auto dirichlet_bc_ptr = boost::make_shared<DirichletDisplacementBc>(
761 m_field, "DISPLACEMENT", Aij, D0, F);
762
763
764
765 dirichlet_bc_ptr->snes_ctx = FEMethod::CTX_SNESNONE;
766 dirichlet_bc_ptr->ts_ctx = FEMethod::CTX_TSNONE;
767
768
769 CHKERR VecZeroEntries(D0);
770 CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
771 CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
773
774
775
777
778
779 CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
780 CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
782
783
784
785
786
787
788
789 fe_rhs_ptr->snes_f = F;
790 prism_fe_rhs_ptr->snes_f = F;
791 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Assemble external force vector ...";
794 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
795
796 fe_lhs_ptr->snes_B = Aij;
797 prism_fe_lhs_ptr->snes_B = Aij;
798 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Calculate stiffness matrix ...";
801 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
802
803
806
807
810
811 if (is_calculating_frequency == PETSC_TRUE) {
812
813 fe_mass_ptr->snes_B = Mij;
814 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Calculate mass matrix ...";
816 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
817 }
818
819
820
821
822
823
824 boost::ptr_map<std::string, NeumannForcesSurface> neumann_forces;
826 F, "DISPLACEMENT");
827
828 {
829 boost::ptr_map<std::string, NeumannForcesSurface>::iterator mit =
830 neumann_forces.begin();
831 for (; mit != neumann_forces.end(); mit++) {
833 &mit->second->getLoopFe());
834 }
835 }
836
837 boost::ptr_map<std::string, NodalForce> nodal_forces;
839 "DISPLACEMENT");
840
841 {
842 boost::ptr_map<std::string, NodalForce>::iterator fit =
843 nodal_forces.begin();
844 for (; fit != nodal_forces.end(); fit++) {
846 &fit->second->getLoopFe());
847 }
848 }
849
850 boost::ptr_map<std::string, EdgeForce> edge_forces;
852 "DISPLACEMENT");
853 {
854 auto fit = edge_forces.begin();
855 for (; fit != edge_forces.end(); fit++) {
856 auto &fe = fit->second->getLoopFe();
858 }
859 }
860
864 CHKERR body_forces_methods.addBlock(
"DISPLACEMENT", F,
865 it->getMeshsetId());
866 }
868 &body_forces_methods.getLoopFe());
869
871 false, false);
872 CHKERR fluid_pressure_fe.setNeumannFluidPressureFiniteElementOperators(
873 "DISPLACEMENT", F, false, true);
874
876 &fluid_pressure_fe.getLoopFe());
877
879
880
881 PetscViewerPushFormat(
882 PETSC_VIEWER_STDOUT_SELF,
883 PETSC_VIEWER_ASCII_MATLAB);
884
885
886
887
888
889
890 if (is_calculating_frequency == PETSC_TRUE) {
891 CHKERR MatAssemblyBegin(Mij, MAT_FINAL_ASSEMBLY);
892 CHKERR MatAssemblyEnd(Mij, MAT_FINAL_ASSEMBLY);
893 }
894
895
896
897
898 CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
899 CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
900 CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
901 CHKERR VecAssemblyBegin(F);
904
905
906 auto solver =
createKSP(PETSC_COMM_WORLD);
907 CHKERR KSPSetDM(solver, dm);
908 CHKERR KSPSetFromOptions(solver);
909 CHKERR KSPSetOperators(solver, Aij, Aij);
910
911
912 {
913
914 PetscBool same = PETSC_FALSE;
915 PC pc;
916 CHKERR KSPGetPC(solver, &pc);
917 PetscObjectTypeCompare((PetscObject)pc, PCMG, &same);
918 if (same) {
921 CHKERR PCSetFromOptions(pc);
922 } else {
923
924 CHKERR KSPSetDMActive(solver, PETSC_FALSE);
925 }
926 }
927 CHKERR KSPSetInitialGuessKnoll(solver, PETSC_FALSE);
928 CHKERR KSPSetInitialGuessNonzero(solver, PETSC_TRUE);
929
931
932 auto set_post_proc_skin = [&](auto &post_proc_skin) {
935 false);
936 CHKERR post_proc_skin.generateReferenceElementMesh();
937 CHKERR post_proc_skin.addFieldValuesPostProc(
"DISPLACEMENT");
938 CHKERR post_proc_skin.addFieldValuesPostProc(
"MESH_NODE_POSITIONS");
939 CHKERR post_proc_skin.addFieldValuesGradientPostProcOnSkin(
940 "DISPLACEMENT", "ELASTIC", data_at_pts->hMat, true);
941 CHKERR post_proc_skin.addFieldValuesGradientPostProcOnSkin(
942
943 "MESH_NODE_POSITIONS", "ELASTIC", data_at_pts->HMat, false);
944 post_proc_skin.getOpPtrVector().push_back(
945 new HookeElement::OpPostProcHookeElement<
947 "DISPLACEMENT", data_at_pts, *block_sets_ptr,
948 post_proc_skin.postProcMesh, post_proc_skin.mapGaussPts, true,
949 true));
951 };
952
953 auto set_post_proc_tets = [&](auto &post_proc) {
955
956 CHKERR post_proc.generateReferenceElementMesh();
958 false);
959 CHKERR post_proc.addFieldValuesPostProc(
"DISPLACEMENT");
960 CHKERR post_proc.addFieldValuesPostProc(
"MESH_NODE_POSITIONS");
961 CHKERR post_proc.addFieldValuesGradientPostProc(
"DISPLACEMENT");
962
964 m_field, post_proc.postProcMesh, post_proc.mapGaussPts,
965 "DISPLACEMENT", post_proc.commonData, block_sets_ptr.get()));
967 };
968
969 auto set_post_proc_edge = [&](auto &post_proc_edge) {
971 CHKERR post_proc_edge.generateReferenceElementMesh();
972 CHKERR post_proc_edge.addFieldValuesPostProc(
"DISPLACEMENT");
974 };
975
976 auto set_post_proc_prisms = [&](auto &prism_post_proc) {
978 CHKERR prism_post_proc.generateReferenceElementMesh();
979 boost::shared_ptr<MatrixDouble> inv_jac_ptr(
new MatrixDouble);
980 prism_post_proc.getOpPtrVector().push_back(
982 prism_post_proc.getOpPtrVector().push_back(
984 CHKERR prism_post_proc.addFieldValuesPostProc(
"DISPLACEMENT");
985 CHKERR prism_post_proc.addFieldValuesPostProc(
"MESH_NODE_POSITIONS");
986 CHKERR prism_post_proc.addFieldValuesGradientPostProc(
"DISPLACEMENT");
988 m_field, prism_post_proc.postProcMesh, prism_post_proc.mapGaussPts,
989 "DISPLACEMENT", prism_post_proc.commonData, block_sets_ptr.get()));
991 };
992
997
998 CHKERR set_post_proc_skin(post_proc_skin);
999 CHKERR set_post_proc_tets(post_proc);
1000 CHKERR set_post_proc_prisms(prism_post_proc);
1001 CHKERR set_post_proc_edge(post_proc_edge);
1002
1003
1005
1006
1008 CHKERR VecDuplicate(F, &F_thermal);
1009
1010
1011
1012 CHKERR thermal_stress_elem.setThermalStressRhsOperators(
1013 "DISPLACEMENT", "TEMP", F_thermal);
1014
1017
1018
1019
1020
1021
1022
1023
1024
1026
1027
1028
1030 sit)) {
1031 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Process step %d",
1032 sit->get_step_number());
1033
1035 sit->get_step_number());
1036
1037 CHKERR VecZeroEntries(F_thermal);
1038 CHKERR VecGhostUpdateBegin(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1039 CHKERR VecGhostUpdateEnd(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1040
1041
1042
1043
1044
1045
1047 dm, "ELASTIC", &thermal_stress_elem.getLoopThermalStressRhs());
1048
1049
1050 CHKERR VecAssemblyBegin(F_thermal);
1051 CHKERR VecAssemblyEnd(F_thermal);
1052
1053 CHKERR VecGhostUpdateBegin(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1054 CHKERR VecGhostUpdateEnd(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1055
1056
1057 PetscReal nrm_F;
1058 CHKERR VecNorm(F, NORM_2, &nrm_F);
1059
1060 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"norm2 F = %6.4e", nrm_F);
1061 PetscReal nrm_F_thermal;
1062 CHKERR VecNorm(F_thermal, NORM_2, &nrm_F_thermal);
1063 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"norm2 F_thermal = %6.4e",
1064 nrm_F_thermal);
1065
1066 CHKERR VecScale(F_thermal, -1);
1067
1068 CHKERR VecAXPY(F_thermal, 1, F);
1069
1070
1071 dirichlet_bc_ptr->snes_x =
D;
1072 dirichlet_bc_ptr->snes_f = F_thermal;
1074
1075
1076 CHKERR KSPSolve(solver, F_thermal,
D);
1077
1079 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
1080 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
1081
1082
1084
1085
1087
1088
1089
1090 if (is_post_proc_volume == PETSC_TRUE) {
1091 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Write output file ...";
1093 std::ostringstream o1;
1094 o1 << "out_" << sit->step_number << ".h5m";
1095 if (!test_nb)
1096 CHKERR post_proc.writeFile(o1.str().c_str());
1097 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done ...";
1098 }
1099
1100 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Write output file skin ...";
1102 &post_proc_skin);
1103 std::ostringstream o1_skin;
1104 o1_skin << "out_skin" << sit->step_number << ".h5m";
1105 if (!test_nb)
1106 CHKERR post_proc_skin.writeFile(o1_skin.str().c_str());
1107 MOFEM_LOG(
"POST_PROC_SKIN", Sev::inform) <<
"done ...";
1108 }
1109 } else {
1110
1111
1112 CHKERR VecZeroEntries(F_thermal);
1113 CHKERR VecGhostUpdateBegin(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1114 CHKERR VecGhostUpdateEnd(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1115
1116
1118 dm, "ELASTIC", &thermal_stress_elem.getLoopThermalStressRhs());
1119
1120
1121 CHKERR VecAssemblyBegin(F_thermal);
1122 CHKERR VecAssemblyEnd(F_thermal);
1123
1124
1125 CHKERR VecGhostUpdateBegin(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1126 CHKERR VecGhostUpdateEnd(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1127
1128
1129 PetscReal nrm_F;
1130 CHKERR VecNorm(F, NORM_2, &nrm_F);
1131
1132 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"norm2 F = %6.4e", nrm_F);
1133 PetscReal nrm_F_thermal;
1134 CHKERR VecNorm(F_thermal, NORM_2, &nrm_F_thermal);
1135
1136 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"norm2 F_thermal = %6.4e",
1137 nrm_F_thermal);
1138
1139
1140 CHKERR VecScale(F_thermal, -1);
1141 CHKERR VecAXPY(F_thermal, 1, F);
1142
1143
1144 dirichlet_bc_ptr->snes_x =
D;
1145 dirichlet_bc_ptr->snes_f = F_thermal;
1147
1148
1149 CHKERR KSPSolve(solver, F_thermal,
D);
1151
1152
1153 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
1154 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
1156
1157
1158 if (is_post_proc_volume == PETSC_TRUE) {
1159 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Write output file ...";
1161
1162 if (!test_nb)
1163 CHKERR post_proc.writeFile(
"out.h5m");
1164 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
1165 }
1166
1167 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Write output file skin ...";
1169 if (!test_nb)
1170 CHKERR post_proc_skin.writeFile(
"out_skin.h5m");
1171 MOFEM_LOG(
"POST_PROC_SKIN", Sev::inform) <<
"done";
1172 }
1173
1174
1175 CHKERR VecDestroy(&F_thermal);
1176 } else {
1177
1178
1179
1180 CHKERR KSPSolve(solver, F,
D);
1181
1182
1183
1184
1185
1187
1188 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
1189 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
1190
1192
1193 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Post-process start ...";
1194 if (is_post_proc_volume == PETSC_TRUE) {
1196 }
1200 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
1201
1202
1203 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Write output file ...";
1204 if (mesh_has_tets) {
1205 if (is_post_proc_volume == PETSC_TRUE) {
1206 if (!test_nb)
1207 CHKERR post_proc.writeFile(
"out.h5m");
1208 }
1209 if (!test_nb)
1210 CHKERR post_proc_skin.writeFile(
"out_skin.h5m");
1211 }
1212 if (mesh_has_prisms) {
1213 if (!test_nb)
1214 CHKERR prism_post_proc.writeFile(
"prism_out.h5m");
1215 }
1216 if (!edges_in_simple_rod.empty())
1217 if (!test_nb)
1218 CHKERR post_proc_edge.writeFile(
"out_edge.h5m");
1219 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
1220 }
1221
1222 if (is_calculating_frequency == PETSC_TRUE) {
1223
1225 VecDuplicate(
D, &u1);
1227 double mode_mass;
1228 CHKERR VecDot(u1,
D, &mode_mass);
1229 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Mode mass %6.4e\n", mode_mass);
1230
1232 VecDuplicate(
D, &v1);
1234
1235 double mode_stiffness;
1236 CHKERR VecDot(v1,
D, &mode_stiffness);
1237 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Mode stiffness %6.4e\n",
1238 mode_stiffness);
1239
1240 double frequency;
1241 double pi = 3.14159265359;
1242 frequency = std::sqrt(mode_stiffness / mode_mass) / (2 * pi);
1243 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Frequency %6.4e", frequency);
1244 }
1245
1246
1247 auto calculate_strain_energy = [&]() {
1249
1251 CHKERR HookeElement::calculateEnergy(dm, block_sets_ptr,
"DISPLACEMENT",
1252 "MESH_NODE_POSITIONS", false, true,
1253 v_energy);
1254
1255
1256 double energy;
1257 CHKERR VecSum(v_energy, &energy);
1258 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Elastic energy %6.4e", energy);
1259
1260 switch (test_nb) {
1261 case 1:
1262 if (fabs(energy - 17.129) > 1e-3)
1264 "atom test diverged!");
1265 break;
1266 case 2:
1267 if (fabs(energy - 5.6475e-03) > 1e-4)
1269 "atom test diverged!");
1270 break;
1271 case 3:
1272 if (fabs(energy - 7.4679e-03) > 1e-4)
1274 "atom test diverged!");
1275 break;
1276 case 4:
1277 if (fabs(energy - 2.4992e+00) > 1e-3)
1279 "atom test diverged!");
1280 break;
1281
1282 case 8: {
1283 double min;
1284 CHKERR VecMin(
D, PETSC_NULL, &min);
1285 constexpr double expected_val = 0.10001;
1286 if (fabs(min + expected_val) > 1e-10)
1288 "atom test diverged! %3.4e != %3.4e", min, expected_val);
1289 } break;
1290 case 9: {
1291 if (fabs(energy - 4.7416e-04) > 1e-8)
1293 "atom test diverged!");
1294 }
1295 default:
1296 break;
1297 }
1298
1300 };
1301 CHKERR calculate_strain_energy();
1302
1303 MPI_Comm_free(&moab_comm_world);
1304 }
1306
1308
1309 return 0;
1310}
const std::string default_options
#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
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
const FTensor::Tensor2< T, Dim, Dim > Vec
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
SmartPetscObj< Mat > smartMatDuplicate(Mat &mat, MatDuplicateOption op)
auto createKSP(MPI_Comm comm)
auto createSmartDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
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.
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Interface for managing meshsets containing materials and boundary conditions.
MoFEMErrorCode printForceSet() const
print meshsets with force boundary conditions data structure
MoFEMErrorCode printMaterialsSet() const
print meshsets with material data structure set on it
MoFEMErrorCode printDisplacementSet() const
print meshsets with displacement boundary conditions data structure
Calculate inverse of jacobian for face element.
Transform local reference derivatives of shape functions to global derivatives.
Projection of edge entities with one mid-node on hierarchical basis.
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Volume finite element base.
Set data structures of MG pre-conditioner via approximation orders.
Operator post-procesing stresses for Hook isotropic material.
Implentation of thermal stress element.