159 const string default_options =
"-ksp_type fgmres \n"
161 "-pc_factor_mat_solver_type mumps\n"
162 "-mat_mumps_icntl_20 0\n"
167 "-snes_type newtonls \n"
168 "-snes_linesearch_type l2 \n"
169 "-snes_linesearch_monitor \n"
173 "-snes_converged_reason \n";
175 string param_file =
"param_file.petsc";
176 if (!
static_cast<bool>(ifstream(param_file))) {
177 std::ofstream file(param_file.c_str(), std::ios::ate);
178 if (file.is_open()) {
179 file << default_options;
190 MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
193 PetscBool flg = PETSC_TRUE;
197 if (flg != PETSC_TRUE) {
199 "*** ERROR -my_file (MESH FILE NEEDED)");
202 PetscScalar step_size_reduction;
204 &step_size_reduction, &flg);
205 if (flg != PETSC_TRUE) {
206 step_size_reduction = 1.;
212 if (flg != PETSC_TRUE) {
218 if (flg != PETSC_TRUE) {
225 if (flg != PETSC_TRUE) {
231 if (std::string(
mesh_file_name).find(
"restart") == std::string::npos) {
241 Tag th_step_size, th_step;
242 double def_step_size = 1;
243 rval = moab.tag_get_handle(
"_STEPSIZE", 1, MB_TYPE_DOUBLE, th_step_size,
244 MB_TAG_CREAT | MB_TAG_MESH, &def_step_size);
245 if (
rval == MB_ALREADY_ALLOCATED)
249 rval = moab.tag_get_handle(
"_STEP", 1, MB_TYPE_INTEGER, th_step,
250 MB_TAG_CREAT | MB_TAG_MESH, &def_step);
251 if (
rval == MB_ALREADY_ALLOCATED)
254 const void *tag_data_step_size[1];
256 CHKERR moab.tag_get_by_ptr(th_step_size, &root, 1, tag_data_step_size);
257 double &step_size = *(
double *)tag_data_step_size[0];
258 const void *tag_data_step[1];
259 CHKERR moab.tag_get_by_ptr(th_step, &root, 1, tag_data_step);
260 int &step = *(
int *)tag_data_step[0];
262 CHKERR PetscPrintf(PETSC_COMM_WORLD,
263 "Start step %D and step_size = %6.4e\n", step,
276 "_MY_REFINEMENT_LEVEL",
sizeof(
BitRefLevel), MB_TYPE_OPAQUE,
277 th_my_ref_level, MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_BYTES,
281 CHKERR m_field.
get_moab().tag_get_by_ptr(th_my_ref_level, &root_meshset, 1,
282 (
const void **)&ptr_bit_level0);
292 std::vector<BitRefLevel> bit_levels;
299 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Insert Interface %d\n",
300 cit->getMeshsetId());
305 CHKERR moab.create_meshset(MESHSET_SET, ref_level_meshset);
307 ->getEntitiesByTypeAndRefLevel(bit_levels.back(),
311 ->getEntitiesByTypeAndRefLevel(bit_levels.back(),
314 Range ref_level_tets;
315 CHKERR moab.get_entities_by_handle(ref_level_meshset, ref_level_tets,
318 CHKERR interface_ptr->
getSides(cubit_meshset, bit_levels.back(),
true,
324 cubit_meshset,
true,
true, 0);
326 CHKERR moab.delete_entities(&ref_level_meshset, 1);
332 ->updateMeshsetByEntitiesChildren(cubit_meshset,
334 cubit_meshset, MBMAXTYPE,
true);
338 bit_level0 = bit_levels.back();
339 problem_bit_level = bit_level0;
366 "ELASTIC",
"MESH_NODE_POSITIONS");
381 "INTERFACE",
"MESH_NODE_POSITIONS");
420 problem_bit_level,
BitRefLevel().set(),
"ELASTIC", MBTET);
422 problem_bit_level,
BitRefLevel().set(),
"INTERFACE", MBPRISM);
429 const double coords[] = {0, 0, 0};
431 Range range_no_field_vertex;
432 range_no_field_vertex.insert(no_field_vertex);
438 range_no_field_vertex);
443 CHKERR moab.create_meshset(MESHSET_SET, meshset_fe_arc_length);
444 CHKERR moab.add_entities(meshset_fe_arc_length, &no_field_vertex, 1);
450 meshset_fe_arc_length,
"ARC_LENGTH",
false);
502 CHKERR m_field.
get_moab().get_entities_by_type(it->meshset, MBTET, tets,
515 "MESH_NODE_POSITIONS");
545 Vec F, F_body_force,
D;
547 "ELASTIC_MECHANICS",
COL, &
F);
549 CHKERR VecDuplicate(
F, &F_body_force);
552 ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>(
"ELASTIC_MECHANICS",
559 boost::ptr_vector<CohesiveInterfaceElement::PhysicalEquation>
565 cout << std::endl << *it << std::endl;
568 string name = it->getName();
570 if (name.compare(0, 11,
"MAT_ELASTIC") == 0) {
572 CHKERR it->getAttributeDataStructure(mydata);
576 }
else if (name.compare(0, 10,
"MAT_INTERF") == 0) {
578 CHKERR it->getAttributeDataStructure(mydata);
581 interface_materials.push_back(
583 interface_materials.back().h = 1;
584 interface_materials.back().youngModulus = mydata.
data.alpha;
585 interface_materials.back().beta = mydata.data.beta;
586 interface_materials.back().ft = mydata.data.ft;
587 interface_materials.back().Gf = mydata.data.Gf;
591 rval = moab.get_entities_by_type(meshset, MBTRI, tris,
true);
594 rval = moab.get_adjacencies(tris, 3,
false, ents3d,
595 moab::Interface::UNION);
597 interface_materials.back().pRisms = ents3d.subset_by_type(MBPRISM);
602 boost::ptr_vector<CohesiveInterfaceElement::PhysicalEquation>::iterator
603 pit = interface_materials.begin();
604 for (; pit != interface_materials.end(); pit++) {
609 boost::shared_ptr<ArcLengthCtx> arc_ctx = boost::shared_ptr<ArcLengthCtx>(
611 boost::scoped_ptr<ArcLengthElement> my_arc_method_ptr(
613 ArcLengthSnesCtx snes_ctx(m_field,
"ELASTIC_MECHANICS", arc_ctx);
618 &my_dirichlet_bc.problemPtr);
619 CHKERR my_dirichlet_bc.iNitialize();
620 boost::shared_ptr<Hooke<adouble>> hooke_adouble_ptr(
new Hooke<adouble>);
621 boost::shared_ptr<Hooke<double>> hooke_double_ptr(
new Hooke<double>);
625 elastic.setOfBlocks[id].iD = id;
628 elastic.setOfBlocks[id].materialDoublePtr = hooke_double_ptr;
629 elastic.setOfBlocks[id].materialAdoublePtr = hooke_adouble_ptr;
632 elastic.setOfBlocks[
id].tEts,
true);
634 false,
false,
false);
636 false,
false,
false);
638 false,
false,
false);
639 CHKERR elastic.setOperators(
"DISPLACEMENT",
"MESH_NODE_POSITIONS",
false,
643 CHKERR cohesive_elements.addOps(
"DISPLACEMENT", interface_materials);
648 MatGetLocalSize(Aij, &
m, &
n);
649 boost::scoped_ptr<ArcLengthMatShell> mat_ctx(
652 CHKERR MatCreateShell(PETSC_COMM_WORLD,
m,
n,
M,
N, (
void *)mat_ctx.get(),
654 CHKERR MatShellSetOperation(ShellAij, MATOP_MULT,
661 CHKERR body_forces_methods.addBlock(
"DISPLACEMENT", F_body_force,
664 CHKERR VecZeroEntries(F_body_force);
665 CHKERR VecGhostUpdateBegin(F_body_force, INSERT_VALUES, SCATTER_FORWARD);
666 CHKERR VecGhostUpdateEnd(F_body_force, INSERT_VALUES, SCATTER_FORWARD);
668 body_forces_methods.getLoopFe());
669 CHKERR VecGhostUpdateBegin(F_body_force, ADD_VALUES, SCATTER_REVERSE);
670 CHKERR VecGhostUpdateEnd(F_body_force, ADD_VALUES, SCATTER_REVERSE);
671 CHKERR VecAssemblyBegin(F_body_force);
672 CHKERR VecAssemblyEnd(F_body_force);
675 boost::ptr_map<std::string, NeumannForcesSurface> neumann_forces;
676 string fe_name_str =
"FORCE_FE";
679 neumann_forces.at(fe_name_str).getLoopFe(),
false,
683 CHKERR neumann_forces.at(fe_name_str)
684 .addForce(
"DISPLACEMENT", arc_ctx->F_lambda, it->getMeshsetId());
686 fe_name_str =
"PRESSURE_FE";
689 neumann_forces.at(fe_name_str).getLoopFe(),
false,
693 CHKERR neumann_forces.at(fe_name_str)
694 .addPressure(
"DISPLACEMENT", arc_ctx->F_lambda, it->getMeshsetId());
697 boost::ptr_map<std::string, NodalForce> nodal_forces;
698 fe_name_str =
"FORCE_FE";
699 nodal_forces.insert(fe_name_str,
new NodalForce(m_field));
702 CHKERR nodal_forces.at(fe_name_str)
703 .addForce(
"DISPLACEMENT", arc_ctx->F_lambda, it->getMeshsetId());
707 CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
708 CHKERR SNESSetApplicationContext(snes, &snes_ctx);
710 CHKERR SNESSetJacobian(snes, ShellAij, Aij,
SnesMat, &snes_ctx);
711 CHKERR SNESSetFromOptions(snes);
714 CHKERR SNESGetKSP(snes, &ksp);
716 CHKERR KSPGetPC(ksp, &pc);
717 boost::scoped_ptr<PCArcLengthCtx> pc_ctx(
719 CHKERR PCSetType(pc, PCSHELL);
720 CHKERR PCShellSetContext(pc, pc_ctx.get());
726 snes_ctx.getComputeRhs();
727 snes_ctx.getPreProcComputeRhs().push_back(&my_dirichlet_bc);
728 snes_ctx.getPreProcComputeRhs().push_back(&pre_post_proc_fe);
730 "INTERFACE", &cohesive_elements.getFeRhs()));
731 loops_to_do_Rhs.push_back(
733 loops_to_do_Rhs.push_back(
735 snes_ctx.getPostProcComputeRhs().push_back(&pre_post_proc_fe);
736 snes_ctx.getPostProcComputeRhs().push_back(&my_dirichlet_bc);
740 snes_ctx.getSetOperators();
741 snes_ctx.getPreProcSetOperators().push_back(&my_dirichlet_bc);
743 "INTERFACE", &cohesive_elements.getFeLhs()));
744 loops_to_do_Mat.push_back(
746 loops_to_do_Mat.push_back(
748 snes_ctx.getPostProcSetOperators().push_back(&my_dirichlet_bc);
750 double gamma = 0.5, reduction = 1;
753 step_size = step_size_reduction;
755 reduction = step_size_reduction;
759 boost::ptr_map<std::string, NeumannForcesSurface>::iterator mit =
760 neumann_forces.begin();
761 CHKERR VecZeroEntries(arc_ctx->F_lambda);
762 CHKERR VecGhostUpdateBegin(arc_ctx->F_lambda, INSERT_VALUES,
764 CHKERR VecGhostUpdateEnd(arc_ctx->F_lambda, INSERT_VALUES, SCATTER_FORWARD);
765 for (; mit != neumann_forces.end(); mit++) {
767 mit->second->getLoopFe());
769 CHKERR VecGhostUpdateBegin(arc_ctx->F_lambda, ADD_VALUES, SCATTER_REVERSE);
770 CHKERR VecGhostUpdateEnd(arc_ctx->F_lambda, ADD_VALUES, SCATTER_REVERSE);
771 CHKERR VecAssemblyBegin(arc_ctx->F_lambda);
772 CHKERR VecAssemblyEnd(arc_ctx->F_lambda);
773 for (std::vector<int>::iterator vit = my_dirichlet_bc.dofsIndices.begin();
774 vit != my_dirichlet_bc.dofsIndices.end(); vit++) {
775 CHKERR VecSetValue(arc_ctx->F_lambda, *vit, 0, INSERT_VALUES);
777 CHKERR VecAssemblyBegin(arc_ctx->F_lambda);
778 CHKERR VecAssemblyEnd(arc_ctx->F_lambda);
780 CHKERR VecDot(arc_ctx->F_lambda, arc_ctx->F_lambda, &arc_ctx->F_lambda2);
781 PetscPrintf(PETSC_COMM_WORLD,
"\tFlambda2 = %6.4e\n", arc_ctx->F_lambda2);
785 "ELASTIC_MECHANICS",
COL,
D, INSERT_VALUES, SCATTER_FORWARD);
787 "ELASTIC_MECHANICS",
"DISPLACEMENT",
"X0_DISPLACEMENT",
COL,
788 arc_ctx->x0, INSERT_VALUES, SCATTER_FORWARD);
790 CHKERR VecNorm(arc_ctx->x0, NORM_2, &x0_nrm);
791 CHKERR PetscPrintf(PETSC_COMM_WORLD,
792 "\tRead x0_nrm = %6.4e dlambda = %6.4e\n", x0_nrm,
794 CHKERR arc_ctx->setAlphaBeta(1, 0);
797 CHKERR arc_ctx->setAlphaBeta(0, 1);
802 CHKERR post_proc.generateReferenceElementMesh();
803 CHKERR post_proc.addFieldValuesPostProc(
"DISPLACEMENT");
804 CHKERR post_proc.addFieldValuesGradientPostProc(
"DISPLACEMENT");
807 m_field, post_proc.postProcMesh, post_proc.mapGaussPts,
"DISPLACEMENT",
808 post_proc.commonData));
810 bool converged_state =
false;
811 for (; step < max_steps; step++) {
814 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Load Step %D step_size = %6.4e\n",
816 CHKERR arc_ctx->setS(step_size);
817 CHKERR arc_ctx->setAlphaBeta(0, 1);
818 CHKERR VecCopy(
D, arc_ctx->x0);
820 CHKERR my_arc_method_ptr->calculate_init_dlambda(&dlambda);
821 CHKERR my_arc_method_ptr->set_dlambda_to_x(
D, dlambda);
822 }
else if (step == 2) {
823 CHKERR arc_ctx->setAlphaBeta(1, 0);
824 CHKERR my_arc_method_ptr->calculate_dx_and_dlambda(
D);
825 CHKERR my_arc_method_ptr->calculate_lambda_int(step_size);
826 CHKERR arc_ctx->setS(step_size);
827 double dlambda = arc_ctx->dLambda;
829 CHKERR VecNorm(arc_ctx->dx, NORM_2, &dx_nrm);
830 CHKERR PetscPrintf(PETSC_COMM_WORLD,
831 "Load Step %D step_size = %6.4e dlambda0 = %6.4e "
832 "dx_nrm = %6.4e dx2 = %6.4e\n",
833 step, step_size, dlambda, dx_nrm, arc_ctx->dx2);
834 CHKERR VecCopy(
D, arc_ctx->x0);
835 CHKERR VecAXPY(
D, 1., arc_ctx->dx);
836 CHKERR my_arc_method_ptr->set_dlambda_to_x(
D, dlambda);
838 CHKERR my_arc_method_ptr->calculate_dx_and_dlambda(
D);
839 CHKERR my_arc_method_ptr->calculate_lambda_int(step_size);
842 step_size *= reduction;
843 CHKERR arc_ctx->setS(step_size);
844 double dlambda = reduction * arc_ctx->dLambda;
845 CHKERR VecScale(arc_ctx->dx, reduction);
847 CHKERR VecNorm(arc_ctx->dx, NORM_2, &dx_nrm);
848 CHKERR PetscPrintf(PETSC_COMM_WORLD,
849 "Load Step %D step_size = %6.4e dlambda0 = %6.4e "
850 "dx_nrm = %6.4e dx2 = %6.4e\n",
851 step, step_size, dlambda, dx_nrm, arc_ctx->dx2);
852 CHKERR VecCopy(
D, arc_ctx->x0);
853 CHKERR VecAXPY(
D, 1., arc_ctx->dx);
854 CHKERR my_arc_method_ptr->set_dlambda_to_x(
D, dlambda);
857 CHKERR SNESSolve(snes, PETSC_NULL,
D);
861 "ELASTIC_MECHANICS",
COL,
D, INSERT_VALUES, SCATTER_REVERSE);
863 cohesive_elements.getFeHistory(), 0,
866 CHKERR my_arc_method_ptr->remove_damaged_prisms_nodes();
869 CHKERR SNESGetIterationNumber(snes, &its);
870 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"number of Newton iterations = %D\n",
873 SNESConvergedReason reason;
874 CHKERR SNESGetConvergedReason(snes, &reason);
877 CHKERR arc_ctx->setAlphaBeta(1, 0);
879 converged_state =
false;
882 if (step > 1 && converged_state) {
883 reduction = pow((
double)its_d / (
double)(its + 1), gamma);
884 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"reduction step_size = %6.4e\n",
890 "ELASTIC_MECHANICS",
COL,
D, INSERT_VALUES, SCATTER_REVERSE);
892 "ELASTIC_MECHANICS",
"DISPLACEMENT",
"X0_DISPLACEMENT",
COL,
893 arc_ctx->x0, INSERT_VALUES, SCATTER_REVERSE);
894 converged_state =
true;
899 PetscFOpen(PETSC_COMM_SELF,
DATAFILENAME,
"a+", &datafile);
900 PetscFPrintf(PETSC_COMM_WORLD, datafile,
"%d %d ", reason, its);
902 CHKERR my_arc_method_ptr->postProcessLoadPath();
909 std::ostringstream ss;
910 ss <<
"out_values_" << step <<
".h5m";
911 CHKERR post_proc.writeFile(ss.str().c_str());
917 "ELASTIC_MECHANICS",
COL,
D, INSERT_VALUES, SCATTER_REVERSE);
922 CHKERR VecDestroy(&F_body_force);
924 CHKERR SNESDestroy(&snes);
925 CHKERR MatDestroy(&ShellAij);