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 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 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Elastic Config",
"none");
152 CHKERR PetscOptionsString(
"-my_file",
"mesh file name",
"",
"mesh.h5m",
154
155 CHKERR PetscOptionsInt(
"-my_order",
"default approximation order",
"",
157
158 CHKERR PetscOptionsEList(
"-base",
"approximation base",
"", list_bases,
159 LASBASETOP, list_bases[choice_base_value],
160 &choice_base_value, PETSC_NULL);
161
162 CHKERR PetscOptionsInt(
"-is_atom_test",
"ctest number",
"", test_nb,
163 &test_nb, PETSC_NULL);
164
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);
169
170 CHKERR PetscOptionsString(
"-my_block_config",
"elastic configure file name",
171 "", "block_conf.in", block_config_file, 255,
172 &flg_block_config);
173
175 "-my_is_calculating_frequency", "set if frequency will be calculated",
176 "", is_calculating_frequency, &is_calculating_frequency, PETSC_NULL);
177
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);
181
182 ierr = PetscOptionsEnd();
184
185
186 if (flg_file != PETSC_TRUE) {
187 SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
188 }
189
190
191 moab::Core mb_instance;
192 moab::Interface &moab = mb_instance;
193
194
195
196
197
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);
201 if (pcomm == NULL)
202 pcomm = new ParallelComm(&moab, moab_comm_world);
203
204
205 if (is_partitioned == PETSC_TRUE) {
206
207
208 const char *option;
209 option = "PARALLEL=READ_PART;"
210 "PARALLEL_RESOLVE_SHARED_ENTS;"
211 "PARTITION=PARALLEL_PARTITION;";
213 } else {
214
215
216
217
218 const char *option;
219 option = "";
221 }
222
223
226
227
233
234 bool mesh_has_tets = false;
235 bool mesh_has_prisms = false;
236 int nb_tets = 0;
237 int nb_hexs = 0;
238 int nb_prisms = 0;
239
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);
243
244 mesh_has_tets = (nb_tets + nb_hexs) > 0;
245 mesh_has_prisms = nb_prisms > 0;
246
247
248
250 bit_level0.set(0);
252 0, 3, bit_level0);
253
254
255
256
257
259 if (
bit->getName().compare(0, 3,
"ROD") == 0) {
261 0, 1, bit_level0);
262 }
263 }
264
265
267 switch (choice_base_value) {
268 case LEGENDRE:
270 break;
271 case LOBATTO:
273 break;
274 case BERNSTEIN_BEZIER:
276 break;
277 case JACOBI:
279 break;
280 default:
282 };
285
286
289
290
291
292
293
296
297
299 CHKERR m_field.
get_moab().get_entities_by_type(0, MBEDGE, all_edges,
true);
300
301
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);
309 }
310 }
311
313 "DISPLACEMENT");
314
315
317
318
319 Range edges_to_set_order;
320 edges_to_set_order = subtract(all_edges, edges_in_simple_rod);
321
322
323
324
332
335 else
337
338
339
340 auto setting_second_order_geometry = [&m_field]() {
342
346 moab::Interface::UNION);
347
348
349
350
351
352
353
354
355
356
357
360
362 };
363 CHKERR setting_second_order_geometry();
364
365
366
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;
377 it)) {
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));
385
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));
392
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));
399
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)
406 ->default_value(0));
407 }
408 po::parsed_options parsed =
409 parse_config_file(ini_file, config_file_options, true);
410 store(parsed, vm);
411 po::notify(vm);
413 it)) {
414 if (block_data[it->getMeshsetId()].oRder == -1)
415 continue;
416 if (block_data[it->getMeshsetId()].oRder ==
order)
417 continue;
418 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Set block %d order to %d",
419 it->getMeshsetId(),
420 block_data[it->getMeshsetId()].oRder);
422 CHKERR moab.get_entities_by_handle(it->getMeshset(), block_ents,
423 true);
424 Range ents_to_set_order;
425 CHKERR moab.get_adjacencies(block_ents, 3,
false,
426 ents_to_set_order,
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,
430 ents_to_set_order,
431 moab::Interface::UNION);
432 CHKERR moab.get_adjacencies(block_ents, 1,
false,
433 ents_to_set_order,
434 moab::Interface::UNION);
436 ents_to_set_order);
437
439 ents_to_set_order, "DISPLACEMENT",
440 block_data[it->getMeshsetId()].oRder);
441 }
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",
449 vit->c_str());
450 }
451 }
452
453
454
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",
466 bd.PoissonRatio);
467 }
468
470 };
471
472
473
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);
478
479 boost::shared_ptr<std::map<int, MassBlockData>> mass_block_sets_ptr =
480 boost::make_shared<std::map<int, MassBlockData>>();
482
483 auto fe_lhs_ptr =
484 boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
485 auto fe_rhs_ptr =
486 boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
487 fe_lhs_ptr->getRuleHook =
VolRule();
488 fe_rhs_ptr->getRuleHook =
VolRule();
489
491 false);
493 false);
494
495 boost::shared_ptr<ForcesAndSourcesCore> prism_fe_lhs_ptr(
497 boost::shared_ptr<ForcesAndSourcesCore> prism_fe_rhs_ptr(
499
500 CHKERR HookeElement::addElasticElement(m_field, block_sets_ptr,
"ELASTIC",
501 "DISPLACEMENT",
502 "MESH_NODE_POSITIONS", false);
503
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);
517 } else {
518 proc_skin = skin_faces;
519 }
522 "DISPLACEMENT");
524 "DISPLACEMENT");
526 "DISPLACEMENT");
528
529
530
531
533 "TEMP");
534 }
536 "POST_PROC_SKIN", "MESH_NODE_POSITIONS");
538 "POST_PROC_SKIN");
540 };
541 CHKERR add_skin_element_for_post_processing();
542
543 auto data_at_pts = boost::make_shared<HookeElement::DataAtIntegrationPts>();
544 if (mesh_has_tets) {
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);
548 }
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);
553 }
554
555 if (test_nb == 4) {
556
557 auto thermal_strain =
560 constexpr double alpha = 1;
565 return t_thermal_strain;
566 };
567
568 fe_rhs_ptr->getOpPtrVector().push_back(
569 new HookeElement::OpAnalyticalInternalStrain_dx<0>(
570 "DISPLACEMENT", data_at_pts, thermal_strain));
571 }
572
573 boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_mass_ptr(
575
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));
582 }
583 }
584
585
586
588 "MESH_NODE_POSITIONS");
589
590
591
592 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr(
594 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr(
596
598 fe_spring_rhs_ptr, "DISPLACEMENT",
599 "MESH_NODE_POSITIONS");
600
601
602
604 "MESH_NODE_POSITIONS");
605
606
607
608
609
610
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));
613
614
616 m_field, fe_simple_rod_lhs_ptr, fe_simple_rod_rhs_ptr, "DISPLACEMENT",
617 "MESH_NODE_POSITIONS");
618
619
620
623 "DISPLACEMENT");
625 "DISPLACEMENT");
627 "DISPLACEMENT");
629 "MESH_NODE_POSITIONS");
630
634 CHKERR m_field.
get_moab().get_entities_by_dimension(it->meshset, 3, tets,
635 true);
637 }
639
640
641
645
646
647
649
650
651 fluid_pressure_fe.addNeumannFluidPressureBCElements("DISPLACEMENT");
652
653
655
657 bool add_temp_field = false;
659 if (block_data[it->getMeshsetId()].initTemp != 0) {
660 add_temp_field = true;
661 break;
662 }
663 }
664 if (add_temp_field) {
667
670 }
671 }
673 CHKERR thermal_stress_elem.addThermalStressElement(
674 "ELASTIC", "DISPLACEMENT", "TEMP");
675 }
676
677
679
680
681
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,
697 "TEMP");
698 }
699 }
700 }
701
702
703
705
707
708
710
711
712 auto dm =
createDM(PETSC_COMM_WORLD,
"MOFEM");
714 CHKERR DMSetFromOptions(dm);
716
726
727
728
735 CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
736
737
739 if (is_calculating_frequency == PETSC_TRUE) {
741 CHKERR MatSetOption(Mij, MAT_SPD, PETSC_TRUE);
742
743 }
744
745
746 fe_spring_lhs_ptr->ksp_B = Aij;
747 fe_spring_rhs_ptr->ksp_f =
F;
748
749
750 fe_simple_rod_lhs_ptr->ksp_B = Aij;
751 fe_simple_rod_rhs_ptr->ksp_f =
F;
752
753
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);
762
763
764
765
766
767
768
769
770
771
772
773 auto dirichlet_bc_ptr = boost::make_shared<DirichletDisplacementBc>(
774 m_field,
"DISPLACEMENT", Aij, D0,
F);
775
776
777
778 dirichlet_bc_ptr->snes_ctx = FEMethod::CTX_SNESNONE;
779 dirichlet_bc_ptr->ts_ctx = FEMethod::CTX_TSNONE;
780
781
782 CHKERR VecZeroEntries(D0);
783 CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
784 CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
786
787
788
790
791
792 CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
793 CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
795
796
797
798
799
800
801
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";
808
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";
815
816
819
820
823
824 if (is_calculating_frequency == PETSC_TRUE) {
825
826 fe_mass_ptr->snes_B = Mij;
827 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Calculate mass matrix ...";
829 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
830 }
831
832
833
834
835
836
837 boost::ptr_map<std::string, NeumannForcesSurface> neumann_forces;
840
841 {
842 boost::ptr_map<std::string, NeumannForcesSurface>::iterator mit =
843 neumann_forces.begin();
844 for (; mit != neumann_forces.end(); mit++) {
846 &mit->second->getLoopFe());
847 }
848 }
849
850 boost::ptr_map<std::string, NodalForce> nodal_forces;
852 "DISPLACEMENT");
853
854 {
855 boost::ptr_map<std::string, NodalForce>::iterator fit =
856 nodal_forces.begin();
857 for (; fit != nodal_forces.end(); fit++) {
859 &fit->second->getLoopFe());
860 }
861 }
862
863 boost::ptr_map<std::string, EdgeForce> edge_forces;
865 "DISPLACEMENT");
866 {
867 auto fit = edge_forces.begin();
868 for (; fit != edge_forces.end(); fit++) {
869 auto &fe = fit->second->getLoopFe();
871 }
872 }
873
877 CHKERR body_forces_methods.addBlock(
"DISPLACEMENT",
F,
878 it->getMeshsetId());
879 }
881 &body_forces_methods.getLoopFe());
882
884 false, false);
885 CHKERR fluid_pressure_fe.setNeumannFluidPressureFiniteElementOperators(
886 "DISPLACEMENT",
F,
false,
true);
887
889 &fluid_pressure_fe.getLoopFe());
890
892
893
894 PetscViewerPushFormat(
895 PETSC_VIEWER_STDOUT_SELF,
896 PETSC_VIEWER_ASCII_MATLAB);
897
898
899
900
901
902
903 if (is_calculating_frequency == PETSC_TRUE) {
904 CHKERR MatAssemblyBegin(Mij, MAT_FINAL_ASSEMBLY);
905 CHKERR MatAssemblyEnd(Mij, MAT_FINAL_ASSEMBLY);
906 }
907
908
909
910
911 CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
912 CHKERR VecGhostUpdateBegin(
F, ADD_VALUES, SCATTER_REVERSE);
913 CHKERR VecGhostUpdateEnd(
F, ADD_VALUES, SCATTER_REVERSE);
917
918
919 auto solver =
createKSP(PETSC_COMM_WORLD);
920 CHKERR KSPSetDM(solver, dm);
921 CHKERR KSPSetFromOptions(solver);
922 CHKERR KSPSetOperators(solver, Aij, Aij);
923
924
925 {
926
927 PetscBool same = PETSC_FALSE;
928 PC pc;
929 CHKERR KSPGetPC(solver, &pc);
930 PetscObjectTypeCompare((PetscObject)pc, PCMG, &same);
931 if (same) {
934 CHKERR PCSetFromOptions(pc);
935 } else {
936
937 CHKERR KSPSetDMActive(solver, PETSC_FALSE);
938 }
939 }
940 CHKERR KSPSetInitialGuessKnoll(solver, PETSC_FALSE);
941 CHKERR KSPSetInitialGuessNonzero(solver, PETSC_TRUE);
942
944
945 auto set_post_proc_skin = [&](auto &post_proc_skin) {
948 false);
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");
959 }
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,
965 true));
967 };
968
969 auto set_post_proc_tets = [&](auto &post_proc) {
971
972 CHKERR post_proc.generateReferenceElementMesh();
974 false);
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");
981 }
982
984 m_field, post_proc.postProcMesh, post_proc.mapGaussPts,
985 "DISPLACEMENT", post_proc.commonData, block_sets_ptr.get()));
987 };
988
989 auto set_post_proc_edge = [&](auto &post_proc_edge) {
991 CHKERR post_proc_edge.generateReferenceElementMesh();
992 CHKERR post_proc_edge.addFieldValuesPostProc(
"DISPLACEMENT");
994 };
995
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()));
1011 };
1012
1017
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);
1022
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,
1029 &field_eval_flag);
1030
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>();
1034
1035 if (field_eval_flag) {
1037 ->getData<VolumeElementForcesAndSourcesCore>();
1039 field_eval_data, "ELASTIC");
1040
1041 field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
1042 auto no_rule = [](
int,
int,
int) {
return -1; };
1043
1044 auto field_eval_fe_ptr = field_eval_data->feMethodPtr.lock();
1045 field_eval_fe_ptr->getRuleHook = no_rule;
1046
1048 field_eval_fe_ptr->getOpPtrVector().push_back(
1050 }
1051 field_eval_fe_ptr->getOpPtrVector().push_back(
1053 field_eval_fe_ptr->getOpPtrVector().push_back(
1055 }
1056
1057
1059
1060
1062 CHKERR VecDuplicate(
F, &F_thermal);
1063
1064
1065
1066 CHKERR thermal_stress_elem.setThermalStressRhsOperators(
1067 "DISPLACEMENT", "TEMP", F_thermal);
1068
1071
1072
1073
1074
1075
1076
1077
1078
1080
1081
1082
1084 sit)) {
1085 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Process step %d",
1086 sit->get_step_number());
1087
1089 sit->get_step_number());
1090
1091 CHKERR VecZeroEntries(F_thermal);
1092 CHKERR VecGhostUpdateBegin(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1093 CHKERR VecGhostUpdateEnd(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1094
1095
1096
1097
1098
1099
1101 dm, "ELASTIC", &thermal_stress_elem.getLoopThermalStressRhs());
1102
1103
1104 CHKERR VecAssemblyBegin(F_thermal);
1105 CHKERR VecAssemblyEnd(F_thermal);
1106
1107 CHKERR VecGhostUpdateBegin(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1108 CHKERR VecGhostUpdateEnd(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1109
1110
1111 PetscReal nrm_F;
1112 CHKERR VecNorm(
F, NORM_2, &nrm_F);
1113
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",
1118 nrm_F_thermal);
1119
1120 CHKERR VecScale(F_thermal, -1);
1121
1122 CHKERR VecAXPY(F_thermal, 1,
F);
1123
1124
1125 dirichlet_bc_ptr->snes_x =
D;
1126 dirichlet_bc_ptr->snes_f = F_thermal;
1128
1129
1130 CHKERR KSPSolve(solver, F_thermal,
D);
1131
1133 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
1134 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
1135
1136
1138
1139
1141
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;
1152 }
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));
1159 }
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);
1165 }
1166
1168 }
1169
1170
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";
1176 if (!test_nb)
1177 CHKERR post_proc.writeFile(o1.str().c_str());
1178 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done ...";
1179 }
1180
1181 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Write output file skin ...";
1183 &post_proc_skin);
1184 std::ostringstream o1_skin;
1185 o1_skin << "out_skin_" << sit->step_number << ".h5m";
1186 if (!test_nb)
1187 CHKERR post_proc_skin.writeFile(o1_skin.str().c_str());
1188 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done ...";
1189 }
1190 } else {
1191
1192
1193 CHKERR VecZeroEntries(F_thermal);
1194 CHKERR VecGhostUpdateBegin(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1195 CHKERR VecGhostUpdateEnd(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
1196
1197
1199 dm, "ELASTIC", &thermal_stress_elem.getLoopThermalStressRhs());
1200
1201
1202 CHKERR VecAssemblyBegin(F_thermal);
1203 CHKERR VecAssemblyEnd(F_thermal);
1204
1205
1206 CHKERR VecGhostUpdateBegin(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1207 CHKERR VecGhostUpdateEnd(F_thermal, ADD_VALUES, SCATTER_REVERSE);
1208
1209
1210 PetscReal nrm_F;
1211 CHKERR VecNorm(
F, NORM_2, &nrm_F);
1212
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);
1216
1217 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"norm2 F_thermal = %6.4e",
1218 nrm_F_thermal);
1219
1220
1221 CHKERR VecScale(F_thermal, -1);
1222 CHKERR VecAXPY(F_thermal, 1,
F);
1223
1224
1225 dirichlet_bc_ptr->snes_x =
D;
1226 dirichlet_bc_ptr->snes_f = F_thermal;
1228
1229
1230 CHKERR KSPSolve(solver, F_thermal,
D);
1232
1233
1234 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
1235 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
1237
1238
1239 if (is_post_proc_volume == PETSC_TRUE) {
1240 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Write output file ...";
1242
1243 if (!test_nb)
1244 CHKERR post_proc.writeFile(
"out.h5m");
1245 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
1246 }
1247
1248 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Write output file skin ...";
1250 if (!test_nb)
1251 CHKERR post_proc_skin.writeFile(
"out_skin.h5m");
1252 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
1253 }
1254
1255
1256 CHKERR VecDestroy(&F_thermal);
1257 } else {
1258
1259
1260
1262
1263
1264
1265
1266
1268
1269 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
1270 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
1271
1273
1274 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Post-process start ...";
1275 if (is_post_proc_volume == PETSC_TRUE) {
1277 }
1281 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
1282
1283
1284 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Write output file ...";
1285 if (mesh_has_tets) {
1286 if (is_post_proc_volume == PETSC_TRUE) {
1287 if (!test_nb)
1288 CHKERR post_proc.writeFile(
"out.h5m");
1289 }
1290 if (!test_nb)
1291 CHKERR post_proc_skin.writeFile(
"out_skin.h5m");
1292 }
1293 if (mesh_has_prisms) {
1294 if (!test_nb)
1295 CHKERR prism_post_proc.writeFile(
"prism_out.h5m");
1296 }
1297 if (!edges_in_simple_rod.empty())
1298 if (!test_nb)
1299 CHKERR post_proc_edge.writeFile(
"out_edge.h5m");
1300 MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done";
1301 }
1302
1303 if (is_calculating_frequency == PETSC_TRUE) {
1304
1306 VecDuplicate(
D, &u1);
1308 double mode_mass;
1309 CHKERR VecDot(u1,
D, &mode_mass);
1310 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Mode mass %6.4e\n", mode_mass);
1311
1313 VecDuplicate(
D, &v1);
1315
1316 double mode_stiffness;
1317 CHKERR VecDot(v1,
D, &mode_stiffness);
1318 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Mode stiffness %6.4e\n",
1319 mode_stiffness);
1320
1321 double frequency;
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);
1325 }
1326
1327
1328 auto calculate_strain_energy = [&]() {
1330
1332 CHKERR HookeElement::calculateEnergy(dm, block_sets_ptr,
"DISPLACEMENT",
1333 "MESH_NODE_POSITIONS", false, true,
1334 v_energy);
1335
1336
1337 double energy;
1338 CHKERR VecSum(v_energy, &energy);
1339 MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Elastic energy %6.4e", energy);
1340
1341 switch (test_nb) {
1342 case 1:
1343 if (fabs(energy - 17.129) > 1e-3)
1345 "atom test diverged!");
1346 break;
1347 case 2:
1348 if (fabs(energy - 5.6475e-03) > 1e-4)
1350 "atom test diverged!");
1351 break;
1352 case 3:
1353 if (fabs(energy - 7.4679e-03) > 1e-4)
1355 "atom test diverged!");
1356 break;
1357 case 4:
1358 if (fabs(energy - 2.4992e+00) > 1e-3)
1360 "atom test diverged!");
1361 break;
1362
1363 case 8: {
1364 double min;
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);
1370 } break;
1371 case 9: {
1372 if (fabs(energy - 4.7416e-04) > 1e-8)
1374 "atom test diverged!");
1375 }
1376 default:
1377 break;
1378 }
1379
1381 };
1382 CHKERR calculate_strain_energy();
1383
1384 MPI_Comm_free(&moab_comm_world);
1385 }
1387
1389
1390 return 0;
1391}
const std::string default_options
#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
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.
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 setBlocks()
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.
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 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.
Set data structures of MG pre-conditioner via approximation orders.
Operator post-procesing stresses for Hook isotropic material.
Implentation of thermal stress element.