23 using namespace MoFEM;
28 static char help[] =
"...\n\n";
34 std::map<int, NonlinearElasticElement::BlockData> &
setOfBlocks;
47 std::map<int, NonlinearElasticElement::BlockData> &set_of_blocks,
50 :
FEMethod(), mField(m_field), postProc(m_field),
51 setOfBlocks(set_of_blocks), feElasticEnergy(fe_elastic_energy),
52 feKineticEnergy(fe_kinetic_energy), iNit(false) {
59 "_TsStep_", 1, MB_TYPE_INTEGER, th_step,
60 MB_TAG_CREAT | MB_TAG_EXCL | MB_TAG_MESH, &def_t_val);
61 if (
rval == MB_ALREADY_ALLOCATED) {
63 (
const void **)&step));
68 (
const void **)&step));
71 PetscBool flg = PETSC_TRUE;
73 "-my_output_prt", &pRT, &flg),
74 "Can not get option");
75 if (flg != PETSC_TRUE) {
93 std::map<int, NonlinearElasticElement::BlockData>::iterator sit =
95 for (; sit != setOfBlocks.end(); sit++) {
104 if ((*step) % pRT == 0) {
106 std::ostringstream sss;
107 sss <<
"out_values_" << (*step) <<
".h5m";
118 double E = feElasticEnergy.eNergy;
119 double T = feKineticEnergy.eNergy;
121 "DYNAMIC", Sev::inform,
122 "%d Time %3.2e Elastic energy %3.2e Kinetic Energy %3.2e Total %3.2e\n",
123 ts_step, ts_t,
E, T,
E + T);
147 double def_t_val = 0;
153 "_TsTime_", 1, MB_TYPE_DOUBLE, th_time,
154 MB_TAG_CREAT | MB_TAG_EXCL | MB_TAG_MESH, &def_t_val);
155 if (
rval == MB_ALREADY_ALLOCATED) {
156 rval = m_field.
get_moab().tag_get_by_ptr(th_time, &root_meshset, 1,
157 (
const void **)&time);
159 ierr = TSSetTime(ts, *time);
160 CHKERRABORT(PETSC_COMM_WORLD,
ierr);
162 rval = m_field.
get_moab().tag_set_data(th_time, &root_meshset, 1,
165 rval = m_field.
get_moab().tag_get_by_ptr(th_time, &root_meshset, 1,
166 (
const void **)&time);
171 "_TsStep_", 1, MB_TYPE_INTEGER, th_step,
172 MB_TAG_CREAT | MB_TAG_EXCL | MB_TAG_MESH, &def_t_val);
173 if (
rval == MB_ALREADY_ALLOCATED) {
174 rval = m_field.
get_moab().tag_get_by_ptr(th_step, &root_meshset, 1,
175 (
const void **)&step);
178 rval = m_field.
get_moab().tag_set_data(th_step, &root_meshset, 1,
181 rval = m_field.
get_moab().tag_get_by_ptr(th_step, &root_meshset, 1,
182 (
const void **)&step);
186 PetscBool flg = PETSC_TRUE;
189 CHKERRABORT(PETSC_COMM_WORLD,
ierr);
190 if (flg != PETSC_TRUE) {
224 #include <TimeForceScale.hpp>
226 int main(
int argc,
char *argv[]) {
228 const string default_options =
"-ksp_type fgmres \n"
230 "-pc_factor_mat_solver_type mumps \n"
231 "-mat_mumps_icntl_20 0 \n"
235 "-snes_type newtonls \n"
236 "-snes_linesearch_type basic \n"
237 "-snes_max_it 100 \n"
243 string param_file =
"param_file.petsc";
244 if (!
static_cast<bool>(ifstream(param_file))) {
245 std::ofstream file(param_file.c_str(), std::ios::ate);
246 if (file.is_open()) {
247 file << default_options;
255 auto core_log = logging::core::get();
266 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
267 auto moab_comm_wrap =
268 boost::make_shared<WrapMPIComm>(PETSC_COMM_WORLD,
false);
270 pcomm =
new ParallelComm(&moab, moab_comm_wrap->get_comm());
274 PetscBool is_partitioned = PETSC_FALSE;
275 PetscBool linear = PETSC_TRUE;
276 PetscInt disp_order = 1;
277 PetscInt vel_order = 1;
278 PetscBool is_solve_at_time_zero = PETSC_FALSE;
280 auto read_command_line_parameters = [&]() {
282 PetscBool flg = PETSC_TRUE;
285 if (flg != PETSC_TRUE)
286 SETERRQ(PETSC_COMM_SELF, 1,
"Error -my_file (mesh file needed)");
291 &is_partitioned, &flg);
296 enum bases { LEGENDRE, LOBATTO, BERNSTEIN_BEZIER, LASBASETOP };
297 const char *list_bases[] = {
"legendre",
"lobatto",
"bernstein_bezier"};
298 PetscInt choice_base_value = BERNSTEIN_BEZIER;
300 LASBASETOP, &choice_base_value, PETSC_NULL);
301 if (choice_base_value == LEGENDRE)
303 else if (choice_base_value == LOBATTO)
305 else if (choice_base_value == BERNSTEIN_BEZIER)
310 if (flg != PETSC_TRUE)
315 if (flg != PETSC_TRUE)
316 vel_order = disp_order;
319 "-my_solve_at_time_zero",
320 &is_solve_at_time_zero, &flg);
325 auto read_mesh = [&]() {
327 if (is_partitioned == PETSC_TRUE) {
330 option =
"PARALLEL=BCAST_DELETE;"
331 "PARALLEL_RESOLVE_SHARED_ENTS;"
332 "PARTITION=PARALLEL_PARTITION;";
342 CHKERR read_command_line_parameters();
352 CHKERR moab.create_meshset(MESHSET_SET, meshset_level0);
367 bool check_if_spatial_field_exist = m_field.
check_field(
"DISPLACEMENT");
392 "DISPLACEMENT", PETSC_NULL,
false,
true);
444 false,
false,
false);
459 string name =
"-my_accelerogram";
460 char time_file_name[255];
463 time_file_name, 255, &flg);
464 if (flg == PETSC_TRUE) {
486 "DAMPER",
"MESH_NODE_POSITIONS");
488 std::map<int, KelvinVoigtDamper::BlockMaterialData>::iterator
bit =
491 bit->second.lInear = linear;
504 elastic.getLoopFeEnergy(),
505 inertia.getLoopFeEnergy());
520 if (!check_if_spatial_field_exist) {
522 "MESH_NODE_POSITIONS");
538 "FLUID_PRESSURE_FE");
543 if (is_partitioned) {
562 "FLUID_PRESSURE_FE");
572 if (is_partitioned) {
591 CHKERR TSCreate(PETSC_COMM_WORLD, &ts);
592 CHKERR TSSetType(ts, TSBEULER);
598 ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>(
"Kuu",
600 CHKERR MatDuplicate(shellAij_ctx->
K, MAT_DO_NOT_COPY_VALUES,
604 D,
"DYNAMICS",
COL, shellAij_ctx->
u,
"Kuu",
COL,
607 D,
"DYNAMICS",
"VELOCITY",
COL, shellAij_ctx->
v,
"Kuu",
"DISPLACEMENT",
615 problem_ptr->
getNbDofsRow(), (
void *)shellAij_ctx, &shell_Aij);
616 CHKERR MatShellSetOperation(shell_Aij, MATOP_MULT,
618 CHKERR MatShellSetOperation(
619 shell_Aij, MATOP_ZERO_ENTRIES,
622 ConvectiveMassElement::ShellMatrixElement shell_matrix_element(m_field);
624 m_field,
"DISPLACEMENT", shellAij_ctx->barK, PETSC_NULL, PETSC_NULL);
627 shell_matrix_element.problemName =
"Kuu";
628 shell_matrix_element.shellMatCtx = shellAij_ctx;
629 shell_matrix_element.DirichletBcPtr = &shell_dirichlet_bc;
630 shell_matrix_element.loopK.push_back(
631 ConvectiveMassElement::ShellMatrixElement::PairNameFEMethodPtr(
632 "ELASTIC", &elastic.getLoopFeLhs()));
634 shell_matrix_element.loopK.push_back(
635 ConvectiveMassElement::ShellMatrixElement::PairNameFEMethodPtr(
636 "ELASTIC", &damper.
feLhs));
640 "MESH_NODE_POSITIONS", linear);
644 shell_matrix_element.loopM.push_back(
645 ConvectiveMassElement::ShellMatrixElement::PairNameFEMethodPtr(
646 "ELASTIC", &inertia.getLoopFeMassLhs()));
649 shell_matrix_element.loopAuxM.push_back(
650 ConvectiveMassElement::ShellMatrixElement::PairNameFEMethodPtr(
651 "ELASTIC", &inertia.getLoopFeMassAuxLhs()));
658 boost::ptr_map<std::string, NeumannForcesSurface> surface_forces;
660 string fe_name_str =
"FORCE_FE";
663 surface_forces.at(fe_name_str).getLoopFe(),
false,
667 CHKERR surface_forces.at(fe_name_str)
668 .addForce(
"DISPLACEMENT", PETSC_NULL, it->getMeshsetId(),
true);
669 surface_forces.at(fe_name_str)
674 boost::ptr_map<std::string, NeumannForcesSurface> surface_pressure;
676 string fe_name_str =
"PRESSURE_FE";
679 surface_pressure.at(fe_name_str).getLoopFe(),
false,
683 CHKERR surface_pressure.at(fe_name_str)
684 .addPressure(
"DISPLACEMENT", PETSC_NULL, it->getMeshsetId(),
true);
685 surface_pressure.at(fe_name_str)
691 boost::ptr_map<std::string, EdgeForce> edge_forces;
693 string fe_name_str =
"FORCE_FE";
694 edge_forces.insert(fe_name_str,
new EdgeForce(m_field));
697 CHKERR edge_forces.at(fe_name_str)
698 .addForce(
"DISPLACEMENT", PETSC_NULL, it->getMeshsetId(),
true);
699 edge_forces.at(fe_name_str).methodsOp.push_back(
new TimeForceScale());
704 boost::ptr_map<std::string, NodalForce> nodal_forces;
706 string fe_name_str =
"FORCE_FE";
707 nodal_forces.insert(fe_name_str,
new NodalForce(m_field));
710 CHKERR nodal_forces.at(fe_name_str)
711 .addForce(
"DISPLACEMENT",
F, it->getMeshsetId(),
true);
712 nodal_forces.at(fe_name_str).methodsOp.push_back(
new TimeForceScale());
718 m_field, ts,
"VELOCITY",
"DISPLACEMENT");
731 auto add_static_rhs = [&](
auto &loops_to_do_Rhs) {
733 loops_to_do_Rhs.push_back(
735 for (
auto fit = surface_forces.begin(); fit != surface_forces.end();
737 loops_to_do_Rhs.push_back(
740 for (
auto fit = surface_pressure.begin(); fit != surface_pressure.end();
742 loops_to_do_Rhs.push_back(
745 for (
auto fit = edge_forces.begin(); fit != edge_forces.end(); fit++) {
746 loops_to_do_Rhs.push_back(
749 for (
auto fit = nodal_forces.begin(); fit != nodal_forces.end(); fit++) {
750 loops_to_do_Rhs.push_back(
754 "FLUID_PRESSURE_FE", &fluid_pressure_fe.
getLoopFe()));
758 CHKERR add_static_rhs(loops_to_do_Rhs);
761 loops_to_do_Rhs.push_back(
778 loopsMonitor.push_back(
780 loopsMonitor.push_back(
789 CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
791 CHKERR TSSetFromOptions(ts);
794 CHKERR TSGetSNES(ts, &snes);
797 CHKERR SNESGetKSP(snes, &ksp);
798 CHKERR KSPSetFromOptions(ksp);
800 CHKERR KSPGetPC(ksp, &pc);
801 CHKERR PCSetType(pc, PCSHELL);
803 CHKERR PCShellSetContext(pc, (
void *)&pc_shell_ctx);
809 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
810 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
812 "DYNAMICS",
COL,
D, INSERT_VALUES, SCATTER_REVERSE);
815 if (is_solve_at_time_zero) {
817 Mat Aij = shellAij_ctx->K;
825 "Kuu",
COL,
D, INSERT_VALUES, SCATTER_FORWARD);
826 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
827 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
829 SnesCtx snes_ctx(m_field,
"Kuu");
832 CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
833 CHKERR SNESSetApplicationContext(snes, &snes_ctx);
836 CHKERR SNESSetFromOptions(snes);
845 CHKERR add_static_rhs(loops_to_do_Rhs);
851 loops_to_do_Mat.push_back(
861 "Kuu",
COL,
D, INSERT_VALUES, SCATTER_FORWARD);
863 CHKERR SNESSolve(snes, PETSC_NULL,
D);
865 CHKERR SNESGetIterationNumber(snes, &its);
866 MOFEM_LOG_C(
"DYNAMIC", Sev::inform,
"number of Newton iterations = %d\n",
871 "Kuu",
COL,
D, INSERT_VALUES, SCATTER_REVERSE);
875 CHKERR SNESDestroy(&snes);
878 if (is_solve_at_time_zero) {
880 "DYNAMICS",
COL,
D, INSERT_VALUES, SCATTER_FORWARD);
881 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
882 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
886 #if PETSC_VERSION_GE(3, 7, 0)
887 CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
890 CHKERR TSGetTime(ts, &ftime);
892 PetscInt steps, snesfails, rejects, nonlinits, linits;
893 CHKERR TSGetTimeStepNumber(ts, &steps);
894 CHKERR TSGetSNESFailures(ts, &snesfails);
895 CHKERR TSGetStepRejections(ts, &rejects);
896 CHKERR TSGetSNESIterations(ts, &nonlinits);
897 CHKERR TSGetKSPIterations(ts, &linits);
899 "steps %d (%d rejected, %D SNES fails), ftime %g, nonlinits "
901 steps, rejects, snesfails, ftime, nonlinits, linits);
906 CHKERR MatDestroy(&shellAij_ctx->K);
907 CHKERR MatDestroy(&shellAij_ctx->M);
908 CHKERR VecScatterDestroy(&shellAij_ctx->scatterU);
909 CHKERR VecScatterDestroy(&shellAij_ctx->scatterV);
910 CHKERR MatDestroy(&shell_Aij);