168 const string default_options =
"-ksp_type fgmres \n"
170 "-pc_factor_mat_solver_type mumps \n"
171 "-mat_mumps_icntl_20 0 \n"
175 "-snes_type newtonls \n"
176 "-snes_linesearch_type basic \n"
177 "-snes_max_it 100 \n"
183 string param_file =
"param_file.petsc";
184 if (!
static_cast<bool>(ifstream(param_file))) {
185 std::ofstream file(param_file.c_str(), std::ios::ate);
186 if (file.is_open()) {
187 file << default_options;
192 SlepcInitialize(&argc, &argv, param_file.c_str(),
help);
199 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
200 auto moab_comm_wrap =
201 boost::make_shared<WrapMPIComm>(PETSC_COMM_WORLD,
false);
203 pcomm =
new ParallelComm(&moab, moab_comm_wrap->get_comm());
205 PetscBool flg = PETSC_TRUE;
209 if (flg != PETSC_TRUE) {
210 SETERRQ(PETSC_COMM_SELF, 1,
"*** ERROR -my_file (MESH FILE NEEDED)");
215 PetscBool is_partitioned = PETSC_FALSE;
217 &is_partitioned, &flg);
219 if (is_partitioned == PETSC_TRUE) {
222 option =
"PARALLEL=BCAST_DELETE;PARALLEL_RESOLVE_SHARED_ENTS;PARTITION="
223 "PARALLEL_PARTITION;";
225 CHKERR pcomm->resolve_shared_ents(0, 3, 0);
226 CHKERR pcomm->resolve_shared_ents(0, 3, 1);
227 CHKERR pcomm->resolve_shared_ents(0, 3, 2);
237 Range CubitSIDESETs_meshsets;
239 SIDESET, CubitSIDESETs_meshsets);
245 CHKERR moab.create_meshset(MESHSET_SET, meshset_level0);
260 bool check_if_spatial_field_exist = m_field.
check_field(
"SPATIAL_POSITION");
273 boost::shared_ptr<Hooke<double>> mat_double =
274 boost::make_shared<Hooke<double>>();
275 boost::shared_ptr<MyMat<adouble>> mat_adouble =
276 boost::make_shared<MyMat<adouble>>();
279 CHKERR elastic.setBlocks(mat_double, mat_adouble);
280 CHKERR elastic.addElement(
"ELASTIC",
"SPATIAL_POSITION");
285 elastic.feRhs.getOpPtrVector().push_back(
287 "D0", elastic.commonData));
288 elastic.feLhs.getOpPtrVector().push_back(
290 "D0", elastic.commonData));
291 CHKERR elastic.setOperators(
"SPATIAL_POSITION");
307 if (flg != PETSC_TRUE) {
332 "MESH_NODE_POSITIONS");
338 CHKERR moab.get_entities_by_type(it->meshset, MBTRI, tris,
true);
345 CHKERR moab.get_entities_by_type(it->meshset, MBTRI, tris,
true);
357 if (!check_if_spatial_field_exist) {
359 "MESH_NODE_POSITIONS");
377 if (is_partitioned) {
392 "ELASTIC_MECHANICS",
ROW, &
F);
397 ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>(
"ELASTIC_MECHANICS",
403 neumann_forces.getLoopSpatialFe();
416 CHKERR VecGhostUpdateBegin(
F, INSERT_VALUES, SCATTER_FORWARD);
417 CHKERR VecGhostUpdateEnd(
F, INSERT_VALUES, SCATTER_FORWARD);
418 CHKERR MatZeroEntries(Aij);
421 "ELASTIC_MECHANICS",
COL,
D, INSERT_VALUES, SCATTER_FORWARD);
422 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
423 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
427 my_Dirichlet_bc.snes_ctx = SnesMethod::CTX_SNESSETFUNCTION;
428 my_Dirichlet_bc.snes_x =
D;
429 my_Dirichlet_bc.snes_f =
F;
432 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
433 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
435 "ELASTIC_MECHANICS",
COL,
D, INSERT_VALUES, SCATTER_REVERSE);
438 boost::ptr_map<std::string, NodalForce> nodal_forces;
439 string fe_name_str =
"FORCE_FE";
440 nodal_forces.insert(fe_name_str,
new NodalForce(m_field));
443 boost::ptr_map<std::string, NodalForce>::iterator fit =
444 nodal_forces.begin();
445 for (; fit != nodal_forces.end(); fit++) {
447 fit->second->getLoopFe());
450 neumann.
snes_ctx = SnesMethod::CTX_SNESSETFUNCTION;
455 elastic.getLoopFeRhs().snes_ctx = SnesMethod::CTX_SNESSETFUNCTION;
456 elastic.getLoopFeRhs().snes_x =
D;
457 elastic.getLoopFeRhs().snes_f =
F;
459 elastic.getLoopFeRhs());
466 my_Dirichlet_bc.snes_ctx = SnesMethod::CTX_SNESSETJACOBIAN;
467 my_Dirichlet_bc.snes_B = Aij;
476 elastic.getLoopFeLhs().snes_ctx = SnesMethod::CTX_SNESSETJACOBIAN;
477 elastic.getLoopFeLhs().snes_B = Aij;
479 elastic.getLoopFeLhs());
486 CHKERR VecGhostUpdateBegin(
F, ADD_VALUES, SCATTER_REVERSE);
487 CHKERR VecGhostUpdateEnd(
F, ADD_VALUES, SCATTER_REVERSE);
489 CHKERR MatAssemblyBegin(Aij, MAT_FINAL_ASSEMBLY);
490 CHKERR MatAssemblyEnd(Aij, MAT_FINAL_ASSEMBLY);
500 CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
501 CHKERR KSPSetOperators(solver, Aij, Aij);
502 CHKERR KSPSetFromOptions(solver);
507 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
508 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
511 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
512 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
515 "ELASTIC_MECHANICS",
"SPATIAL_POSITION",
"D0",
COL,
D, INSERT_VALUES,
519 CHKERR MatDuplicate(Aij, MAT_SHARE_NONZERO_PATTERN, &Bij);
521 CHKERR MatZeroEntries(Bij);
540 mat_adouble->doAotherwiseB =
false;
542 my_Dirichlet_bc.snes_ctx = SnesMethod::CTX_SNESSETJACOBIAN;
543 my_Dirichlet_bc.snes_B = Bij;
547 neumann.
snes_ctx = SnesMethod::CTX_SNESSETJACOBIAN;
549 PetscBool is_conservative = PETSC_FALSE;
551 &is_conservative, &flg);
552 if (is_conservative) {
557 elastic.getLoopFeLhs().snes_ctx = SnesMethod::CTX_SNESSETJACOBIAN;
558 elastic.getLoopFeLhs().snes_B = Bij;
560 elastic.getLoopFeLhs());
565 CHKERR MatSetOption(Bij, MAT_SPD, PETSC_TRUE);
566 CHKERR MatAssemblyBegin(Bij, MAT_FINAL_ASSEMBLY);
567 CHKERR MatAssemblyEnd(Bij, MAT_FINAL_ASSEMBLY);
579 PetscInt nev, maxit, its;
584 CHKERR EPSCreate(PETSC_COMM_WORLD, &
eps);
602 PetscPrintf(PETSC_COMM_WORLD,
" Number of iterations of the method: %D\n",
609 PetscPrintf(PETSC_COMM_WORLD,
" Solution method: %s\n",
type);
610 CHKERR EPSGetDimensions(
eps, &nev, NULL, NULL);
611 PetscPrintf(PETSC_COMM_WORLD,
" Number of requested eigenvalues: %D\n",
614 PetscPrintf(PETSC_COMM_WORLD,
" Stopping condition: tol=%.4g, maxit=%D\n",
619 CHKERR post_proc.generateReferenceElementMesh();
620 CHKERR post_proc.addFieldValuesGradientPostProc(
"SPATIAL_POSITION");
621 CHKERR post_proc.addFieldValuesPostProc(
"SPATIAL_POSITION");
622 CHKERR post_proc.addFieldValuesPostProc(
"MESH_NODE_POSITIONS");
623 CHKERR post_proc.addFieldValuesPostProc(
"EIGEN_VECTOR");
624 CHKERR post_proc.addFieldValuesGradientPostProc(
"EIGEN_VECTOR");
625 CHKERR post_proc.addFieldValuesPostProc(
"D0");
626 CHKERR post_proc.addFieldValuesGradientPostProc(
"D0");
629 CHKERR post_proc.writeFile(
"out.h5m");
631 PetscScalar eigr, eigi, nrm2r;
632 for (
int nn = 0; nn < nev; nn++) {
633 CHKERR EPSGetEigenpair(
eps, nn, &eigr, &eigi,
D, PETSC_NULL);
634 CHKERR VecNorm(
D, NORM_2, &nrm2r);
637 " ncov = %D eigr = %.4g eigi = %.4g (inv eigr = %.4g) nrm2r = %.4g\n",
638 nn, eigr, eigi, 1. / eigr, nrm2r);
639 std::ostringstream o1;
640 o1 <<
"eig_" << nn <<
".h5m";
642 "ELASTIC_MECHANICS",
"SPATIAL_POSITION",
"EIGEN_VECTOR",
COL,
D,
643 INSERT_VALUES, SCATTER_REVERSE);
646 CHKERR post_proc.writeFile(o1.str().c_str());
649 CHKERR KSPDestroy(&solver);