diff --git a/FEM/rf_msp_new.cpp b/FEM/rf_msp_new.cpp index 260936851..34bb9121c 100644 --- a/FEM/rf_msp_new.cpp +++ b/FEM/rf_msp_new.cpp @@ -985,62 +985,104 @@ std::ios::pos_type CSolidProperties::Read(std::ifstream* msp_file) in_sd.clear(); } - if (line_string.find("$MICRO_STRUCTURE_PLAS") != - string::npos) // WX:09.2011 for anisotropic plasticity - { - if (line_string.find("NORETURNMAPPING") != string::npos) - Plasticity_type = - 44; // direct stress intergration anisotropic plas. - Plasticity_Bedding = true; - MicroStruTensor = new double[3]; - Bedding_Norm = new double[3]; - TransMicroStru = new Matrix(6, 6); - TransMicroStru_T = new Matrix(6, 6); - TransMicroStru_TInv = new Matrix(6, 6); - for (i = 0; i < 4; i++) - { - in_sd.str(GetLineFromFile1(msp_file)); - in_sd >> line_string; - if (line_string.find("MICRO_STRUCTURE_TENSOR") != string::npos) - { - in_sd >> MicroStruTensor[0] >> MicroStruTensor[1] >> - MicroStruTensor[2]; - } - else if (line_string.find("BEDDING_NORM") != string::npos) - { - in_sd >> Bedding_Norm[0] >> Bedding_Norm[1] >> - Bedding_Norm[2]; - } - /*else if(line_string.find("FRICTION_CURVE")!=string::npos) - { - in_sd>>bedding_fric_curve; - }*/ - else if (line_string.find("UNIAXI_COMP_CURVE") != string::npos) - { - in_sd >> bedding_uc_curve_order; - comp_para = new double[bedding_uc_curve_order + 1]; - for (int ii = 0; ii < bedding_uc_curve_order + 1; ii++) - in_sd >> comp_para[ii]; - } - else if (line_string.find("TENSION_CURVE") != string::npos) - { - in_sd >> bedding_tens_curve_order; - tens_para = new double[bedding_tens_curve_order + 1]; - for (int ii = 0; ii < bedding_tens_curve_order + 1; ii++) - in_sd >> tens_para[ii]; - } - in_sd.clear(); - } - CalTransMatrixMicroStru(TransMicroStru, Bedding_Norm); - TransMicroStru->GetTranspose(*TransMicroStru_T); - Cal_Inv_Matrix(6, TransMicroStru_T, TransMicroStru_TInv); - // TransMicroStru_T->Write(); - // TransMicroStru_TInv->Write(); - // TransMicroStru->Write(); - } - in_sd.clear(); - } - return position; + if (line_string.find("$MICRO_STRUCTURE_PLAS") != string::npos) // WX:09.2011 for anisotropic plasticity + { + if (line_string.find("NORETURNMAPPING") != string::npos) + Plasticity_type = 44; // direct stress intergration anisotropic plas. + Plasticity_Bedding = true; + MicroStruTensor = new double[3]; + Bedding_Norm = new double[3]; + TransMicroStru = new Matrix(6, 6); + TransMicroStru_T = new Matrix(6, 6); + TransMicroStru_TInv = new Matrix(6, 6); + for (i = 0; i < 4; i++) + { + in_sd.str(GetLineFromFile1(msp_file)); + in_sd >> line_string; + if (line_string.find("MICRO_STRUCTURE_TENSOR") != string::npos) + { + in_sd >> MicroStruTensor[0] >> MicroStruTensor[1] >> MicroStruTensor[2]; + } + else if (line_string.find("BEDDING_NORM") != string::npos) + { + in_sd >> Bedding_Norm[0] >> Bedding_Norm[1] >> Bedding_Norm[2]; + } + /*else if(line_string.find("FRICTION_CURVE")!=string::npos) + { + in_sd>>bedding_fric_curve; + }*/ + else if (line_string.find("UNIAXI_COMP_CURVE") != string::npos) + { + in_sd >> bedding_uc_curve_order; + comp_para = new double[bedding_uc_curve_order + 1]; + for (int ii = 0; ii < bedding_uc_curve_order + 1; ii++) + in_sd >> comp_para[ii]; + } + else if (line_string.find("TENSION_CURVE") != string::npos) + { + in_sd >> bedding_tens_curve_order; + tens_para = new double[bedding_tens_curve_order + 1]; + for (int ii = 0; ii < bedding_tens_curve_order + 1; ii++) + in_sd >> tens_para[ii]; + } + in_sd.clear(); + } + CalTransMatrixMicroStru(TransMicroStru, Bedding_Norm); + TransMicroStru->GetTranspose(*TransMicroStru_T); + Cal_Inv_Matrix(6, TransMicroStru_T, TransMicroStru_TInv); + // TransMicroStru_T->Write(); + // TransMicroStru_TInv->Write(); + // TransMicroStru->Write(); + } + in_sd.clear(); + if (line_string.find("$WEAKNESS_PLANE") != string::npos) // LU:04.2019 for slip on prescribed plane (only failure, not checked for general failure on the element) + { + Plasticity_type = 45; // direct stress intergration anisotropic plas. + Plasticity_Bedding = true; + MicroStruTensor = new double[3]; + Bedding_Norm = new double[3]; + TransMicroStru = new Matrix(6, 6); + TransMicroStru_T = new Matrix(6, 6); + TransMicroStru_TInv = new Matrix(6, 6); + for (i = 0; i < 4; i++) // get orientation of the joint + { + in_sd.str(GetLineFromFile1(msp_file)); + in_sd >> line_string; + if (line_string.find("MICRO_STRUCTURE_TENSOR") != string::npos) + { + in_sd >> MicroStruTensor[0] >> MicroStruTensor[1] >> MicroStruTensor[2]; + } + else if (line_string.find("WEAKPLANE_NORM") != string::npos) //same role as BEDDING_NORM previously defined + { + in_sd >> Bedding_Norm[0] >> Bedding_Norm[1] >> Bedding_Norm[2]; + } + in_sd.clear(); + } + // get properties of the joint + Size = 6; + /* + i parameter + 0 cohesion of joint + 1 friction angle of joint + 2 dilatance angle of joint + 3 tension strength of joint + 4 hardening curve for friction angle of joint + 5 hardening curve for cohesion of joint + */ + data_Plasticity_joint = new Matrix(Size_j); + for (i = 0; i < Size; i++) + { + in_sd.str(GetLineFromFile1(msp_file)); + in_sd >> (*data_Plasticity_joint)(i); + in_sd.clear(); + } + CalTransMatrixMicroStru(TransMicroStru, Bedding_Norm); + TransMicroStru->GetTranspose(*TransMicroStru_T); + Cal_Inv_Matrix(6, TransMicroStru_T, TransMicroStru_TInv); + } + in_sd.clear(); + } + return position; } //========================================================================== @@ -1223,6 +1265,8 @@ CSolidProperties::~CSolidProperties() delete data_Youngs; if (data_Plasticity) delete data_Plasticity; + if (data_Plasticity_joint) + delete data_Plasticity_joint; if (data_Capacity) delete data_Capacity; if (data_Conductivity) @@ -1239,6 +1283,7 @@ CSolidProperties::~CSolidProperties() data_Density = NULL; data_Youngs = NULL; data_Plasticity = NULL; + data_Plasticity_joint = NULL; data_Capacity = NULL; data_Conductivity = NULL; data_Creep = NULL; @@ -2764,6 +2809,24 @@ void CSolidProperties::CalculateCoefficent_MOHR( } } +void CSolidProperties::CalculateCoefficent_MOHRjoint( + double ep, double scalar_comp, double scalar_tens) // WX:11.2010, +// 09.2011 +{ + int valid = 1; + double thetaj = (*data_Plasticity_joint)(1) * PI / 180; + double phij = (*data_Plasticity_joint)(2) * PI / 180; + Y0j = (*data_Plasticity_joint)(0); + tensionj = (*data_Plasticity_joint)(3); + + if ((*data_Plasticity_joint)(5) > 0 && (*data_Plasticity_joint)(5) < 100) + thetaj = + GetCurveValue((int)(*data_Plasticity_joint)(5), 0, ep, &valid) * PI / 180; + if ((*data_Plasticity_joint)(4) > 0 && (*data_Plasticity_joint)(4) < 100) + Y0j = GetCurveValue((int)(*data_Plasticity)(4), 0, ep, &valid); + + +} void CSolidProperties::CalculateCoefficent_HOEKBROWN() // WX: 02.2011 { HoekB_a = (*data_Plasticity)(0); @@ -3684,6 +3747,804 @@ void CSolidProperties::TangentialDPwithTensionCorner(Matrix* Dep, double /*mm*/) // Dep->Write(); } +/******************************************************* +LU: Mohr Coulomb and Bedding plane = prescribed shear on one plane (anisotropy of parameters mantained for legacy) +routine is based on DirectStressIntegrationMOHR, first it checks and compute for rupture of the matrix +(DirectStressIntegrationMOHR with plasticity_beddings check removed) and then it checks for rupture on the +weakness plane (stress rotated and projected parallel and perpendicular to the weakplane normal vector) +*******************************************************/ +int CSolidProperties::StressIntegrationMOHR_Joint(const int GPiGPj, const ElementValue_DM* ele_val, double* TryStress, + const int Update, Matrix* Dep) +{ + int i, j, counter; + int yield = 0; + int Dim = 2; + int Size = std::min(std::size_t(6ul), ele_val->Stress->Rows()); + double TmpValue1, TmpValue2, dstrNorm = 0; + ; + //double LodeAngle, I1, J2, J3, sqrtJ2; + double shearsurf, tensionsurf, ep, dlamda = 0.0 /*, dlamda_0, dlamda_1*/; //, ddlamda, Jacob; + + // initialize all vectors + double dstrs[6] = { 0. }, TryStress_0[6] = { 0. }, TryStr_buff[6] = { 0. }, TmpStress0_m[6], TmpStress[6] = { 0. }, tmp_prin_str[6] = { 0. }, tmp_prin_dir[9] = { 0. }; + double dstrainP[6] = { 0. }, dstressP[6] = { 0. }; // TmpStress_m[6] = { 0. }, TmpStress0[6] = { 0. }, devStr[6] = { 0. }, + double prin_str[6] = { 0. }, prin_str0[6] = { 0. }, prin_dir[9] = { 0. }; + //double sqrt3 = sqrt(3.0); + //bool lode_out_range = false; + double TryStress_local[6] = { 0, }, dstrs_local[6] = { 0. }; + + *TmpDe = (0.); + *Inv_De = (0.); + *TMatrix = (0.); + *TmpMatrix = (0.); + *TmpMatrix2 = (0.); + + *TmpDe = *Dep; + + // ConstitutiveMatrix->resize(Size,Size); //in head already defined, and is used for later as global variable + + *ConstitutiveMatrix = (0.); + + ep = (*ele_val->pStrain)(GPiGPj); // get eff plas strain + double id_fs_n[6] = { 0, 0, 1, 0, 0, 0 }; // to remove dilation effect on stress {0, 0, 0, 0, 0, 0}; + double id_fs_s[6] = { 0, 0, 0, 0, 1, 1 }; + //double Kronecker[6] = { 1, 1, 1, 0, 0, 0 }; + + if (Size > 4) + Dim = 3; + + for (i = 0; i < Size; i++) + { + dstrs[i] = TryStress[i]; // d_stress + TryStress[i] = (*ele_val->Stress)(i, GPiGPj); // stress_0 + //devS[i] = TryStress[i] + dstrs[i]; + TryStr_buff[i] = TryStress[i] + dstrs[i]; + TryStress_0[i] = TryStr_buff[i]; + //TmpStress_m[i] = TryStr_buff[i]; + dstrNorm += dstrs[i] * dstrs[i]; + } + //if (Size == 4) + // TryStr_buff[4] = TryStr_buff[5] = 0.; + //if (dstrNorm < MKleinsteZahl) + // return 0; + + CalPrinStrDir(TryStr_buff, prin_str, prin_dir, Dim); + // evaluate Mohr_Coulomb on the matrix ; + CalTransMatrixA(prin_dir, TransMatrixA, Size); + TransMatrixA->GetTranspose(*TransMatrixA_T); + Cal_Inv_Matrix(Size, TransMatrixA, Inv_TransA); + Inv_TransA->GetTranspose(*Inv_TransA_T); + // Inv_TransA_T->multi(*Dep, *Inv_TransA, *PrinDe); //De in prin. (A^-T De A^-1) + *TmpDe = (0.); + Inv_TransA_T->multi(*Dep, *Inv_TransA, *TmpDe); // De in prin. coord. + for (i = 0; i < Size; i++) + { + TmpStress0_m[i] = TryStress[i]; + } + if (Size == 4) + TmpStress0_m[4] = TmpStress0_m[5] = 0.; + CalPrinStrDir(TmpStress0_m, tmp_prin_str, tmp_prin_dir, Dim); // prin. stresses t0 + + CalculateCoefficent_MOHR(ep, 0, 0); + + shearsurf = Ntheta * prin_str[0] - prin_str[2] - csn; + tensionsurf = prin_str[0] - tension; + + /* */ + //beginning of tension or shear rupture in the matrix + if (dstrNorm == 0) + { + shearsurf = -1; + tensionsurf = -1; + } + //std::cout << "Shearsurf " << shearsurf << " tensionsurf " << tensionsurf < 0) || (tensionsurf > 0)) + { + std::cout << " Stating plastic matrix calculation for element with GPiGPj " << GPiGPj << endl; + double P_12, P_31, P_41, P_63, P_64, P_52, P_85, P_74, P_78, P_98, P_45 /*, P_X7*/; + double t1, t2 /*, t1ra*/, t1r1 /*, t2ra, t2r2*/, t3r1 /*, t3r2*/; + double fkt1, fkt2; + double mm = 0.; + + double l1[6] = { 0. }, l2[6] = { 0. }, l1g[6] = { 0. }, l2g[6] = { 0. }, l1R[6] = { 0. }, l2R[6] = { 0. }, l3R[6] = { 0. }; + double rsp[6] = { 0. }, rtp[6] = { 0. }; + double sigA[6] = { 0. }, sig1R[6] = { 0. }, sig2R[6] = { 0. }, sigaR[6] = { 0. }; + double dFsdprin_s[6] = { 0. }, dFtdprin_s[6] = { 0. }, dGsdprin_s[6] = { 0. }, dGtdprin_s[6] = { 0. }; + double De_dGsdprin_s[6] = { 0. }, De_dGtdprin_s[6] = { 0. }; + double dStressP[6] = { 0. }, dStrainP[6] = { 0. }; + double cVec[6] = { 0 }; + + yield = 1; + fkt1 = 0.001; + fkt2 = 0.01; + sigA[0] = sigA[1] = sigA[2] = csn / (Ntheta - 1); + sig1R[0] = sig1R[1] = sig2R[0] = tension; + sig1R[2] = sig2R[1] = sig2R[2] = Ntheta * tension - csn; + sigaR[0] = sigaR[1] = sigaR[2] = tension; + + l1[0] = l1[1] = l1g[0] = l1g[1] = l2[0] = l2g[0] = 1.; + l1[2] = l2[1] = l2[2] = Ntheta; + l1g[2] = l2g[1] = l2g[2] = Nphi; + l1R[2] = l2R[1] = l2R[2] = l3R[1] = 1.; + + dFsdprin_s[0] = Ntheta; + dGsdprin_s[0] = Nphi; + dFsdprin_s[2] = dGsdprin_s[2] = -1.; + TmpDe->multi(dGsdprin_s, De_dGsdprin_s); + // PrinDe->multi(dGsdprin_s, De_dGsdprin_s);// + TmpValue1 = 0.; + for (i = 0; i < Size; i++) + TmpValue1 += dFsdprin_s[i] * De_dGsdprin_s[i]; + for (i = 0; i < Size; i++) + rsp[i] = De_dGsdprin_s[i] / TmpValue1; + + dFtdprin_s[0] = dGtdprin_s[0] = 1.; + TmpDe->multi(dGtdprin_s, De_dGtdprin_s); + // PrinDe->multi(dGtdprin_s, De_dGtdprin_s);// + TmpValue2 = 0.; + for (i = 0; i < Size; i++) + TmpValue2 += dFtdprin_s[i] * De_dGtdprin_s[i]; + for (i = 0; i < Size; i++) + rtp[i] = De_dGtdprin_s[i] / TmpValue2; + + double tmp_shearsurf = Ntheta * tmp_prin_str[0] - tmp_prin_str[2] - csn; + double tmp_tensionsurf = tmp_prin_str[0] - tension; + if (((tmp_tensionsurf) == 0 && (tmp_shearsurf) <= 0) || ((tmp_tensionsurf) <= 0 && (tmp_shearsurf) == 0)) + mm = 0.; + else if (prin_str[0] != tmp_prin_str[0]) + { + mm = (tension - tmp_prin_str[0]) / (prin_str[0] - tmp_prin_str[0]); + if (mm >= 0 && mm <= 1) + { + double tmp_prin_str_3 = tmp_prin_str[2] + mm * (prin_str[2] - tmp_prin_str[2]); + if (tmp_prin_str_3 < (Ntheta * tension - csn) || tmp_prin_str_3 > tension) + mm = (csn + tmp_prin_str[2] - Ntheta * tmp_prin_str[0]) + / (Ntheta * (prin_str[0] - tmp_prin_str[0]) - (prin_str[2] - tmp_prin_str[2])); + } + else + mm = (csn + tmp_prin_str[2] - Ntheta * tmp_prin_str[0]) + / (Ntheta * (prin_str[0] - tmp_prin_str[0]) - (prin_str[2] - tmp_prin_str[2])); + } + else + mm = (csn + tmp_prin_str[2] - Ntheta * tmp_prin_str[0]) + / (Ntheta * (prin_str[0] - tmp_prin_str[0]) - (prin_str[2] - tmp_prin_str[2])); + + Cal_Inv_Matrix(Size, TmpDe, Inv_De); + // Cal_Inv_Matrix(Size, PrinDe, Inv_De); + + t1 = CalVar_t(l1, l1g, Inv_De, prin_str, sig1R, Size); + t2 = CalVar_t(l2, l2g, Inv_De, prin_str, sig2R, Size); + t1r1 = CalVar_t(l1R, l1R, Inv_De, prin_str, sig1R, Size); + // t1ra = CalVar_t(l1R, l1R, Inv_De, prin_str, sigaR, Size); + // t2r2 = CalVar_t(l2R, l2R, Inv_De, prin_str, sig2R, Size); + // t2ra = CalVar_t(l2R, l2R, Inv_De, prin_str, sigaR, Size); + t3r1 = CalVar_t(l3R, l3R, Inv_De, prin_str, sig1R, Size); + // t3r2 = CalVar_t(l3R, l3R, Inv_De, prin_str, sig2R, Size); + + P_12 = CalVarP(rsp, l1, prin_str, sig1R); + P_31 = CalVarP(rsp, l2, prin_str, sig2R); + P_41 = CalVarP(rsp, l3R, prin_str, sig1R); + // P_63 = 0; + // P_63 = CalVarP(rsp, l2R, prin_str, sig2R); + P_63 = P_41; + P_45 = CalVarP(rsp, rtp, prin_str, sig1R); + // P_64 = CalVarP(rsp, rtp, prin_str, sig2R); + P_64 = CalVarP(rsp, rtp, prin_str, sig2R); + // P_52 = 0; + // P_52 = CalVarP(rsp, l1R, prin_str, sig1R); + P_52 = P_41; + P_74 = CalVarP(rtp, l3R, prin_str, sig1R); + // P_85 = 0; + // P_85 = CalVarP(rtp, l1, prin_str, sig1R); + P_85 = P_74; + P_78 = CalVarP(rtp, l1R, prin_str, sig1R); + // P_98 = 0; + P_98 = CalVarP(rtp, l2R, prin_str, sig2R); + // P_X7 = CalVarP(rtp, l2R, prin_str, sig2R); + + if (P_12 >= 0 && P_31 <= 0 && P_41 <= 0) // return to fmc + { + double tmpvalue = 0.; + // Matrix *tmpMatrix2 = new Matrix (Size,Size); + // Matrix *dGds_dFds = new Matrix (Size,Size); + *TmpMatrix2 = (0.); // WX:08.2011 + *dGds_dFds = (0.); // WX:08.2011 + for (i = 0; i < 3; i++) + { + tmpvalue += dFsdprin_s[i] * (prin_str[i] - sigA[i]); + } + for (i = 0; i < 3; i++) + { + prin_str0[i] = prin_str[i]; + prin_str[i] -= tmpvalue * rsp[i]; + dStressP[i] = prin_str[i] - prin_str0[i]; + } + Inv_De->multi(dStressP, dStrainP); + for (i = 0; i < 3; i++) + for (j = 0; j < 3; j++) + (*dGds_dFds)(i, j) = rsp[i] * dFsdprin_s[j]; //(D*dGds/(dFds*D*dGds))*dFds + dGds_dFds->multi(*TmpDe, *TmpMatrix2); + // dGds_dFds->multi(*PrinDe, *tmpMatrix2); + *TmpDe -= *TmpMatrix2; + //*PrinDe -= *tmpMatrix2; + // for(i=3; imulti(dStressP, dStrainP); // dstrainP = D-1 * dstressp + VecCrossProduct(dStrainP, l1g, cVec); // cVec = dstrainP X l1g + + CalDep_l(l1, l1g, Inv_De, Dep_l, 1.0); // Dep_l = l*lgT/(lT*D-1*lg) + CalDep_l(cVec, cVec, Inv_De, dDep_l, fkt2); // dDep_l = fkt2*cVec*cVecT/(cVecT*D-1*cVec) + //*TmpDe=(0.); + for (i = 0; i < 3; i++) + for (j = 0; j < 3; j++) + (*TmpDe)(i, j) = (*Dep_l)(i, j) + (*dDep_l)(i, j); + // for(i=3; i0 && t2<0 ) //return to l2 + else if (P_31 > 0 && P_63 < 0) // return to l2 + { + for (i = 0; i < 3; i++) + { + prin_str0[i] = prin_str[i]; + prin_str[i] = t2 * l2[i] + sig2R[i]; + dStressP[i] = prin_str[i] - prin_str0[i]; + } + Inv_De->multi(dStressP, dStrainP); + VecCrossProduct(dStrainP, l2g, cVec); + CalDep_l(l2, l2g, Inv_De, Dep_l, 1.0); + CalDep_l(cVec, cVec, Inv_De, dDep_l, fkt2); + //*TmpDe=(0.); + for (i = 0; i < 3; i++) + for (j = 0; j < 3; j++) + (*TmpDe)(i, j) = (*Dep_l)(i, j) + (*dDep_l)(i, j); + // for(i=3; i0 && P_74<0 && t3r2>0 && t3r1<0 ) //return to l3R + else if (P_41 > 0 && P_45 > 0 && P_64 < 0 && P_74 < 0) // return to l3R + { + for (i = 0; i < 3; i++) + { + prin_str0[i] = prin_str[i]; + prin_str[i] = t3r1 * l3R[i] + sig1R[i]; + dStressP[i] = prin_str[i] - prin_str0[i]; + } + Inv_De->multi(dStressP, dStrainP); + VecCrossProduct(dStrainP, l3R, cVec); + CalDep_l(l3R, l3R, Inv_De, Dep_l, 1.0); + CalDep_l(cVec, cVec, Inv_De, dDep_l, fkt2); + //*TmpDe=(0.); + for (i = 0; i < 3; i++) + for (j = 0; j < 3; j++) + (*TmpDe)(i, j) = (*Dep_l)(i, j) + (*dDep_l)(i, j); + // for(i=3; i= 0 && P_78 >= 0) // return to ft, && P_X7<=0 + { + double tmpvalue = 0.; + // Matrix *tmpMatrix2 = new Matrix (Size,Size); + // Matrix *dGds_dFds = new Matrix (Size,Size); + *TmpMatrix2 = (0.); // WX:08.2011 + *dGds_dFds = (0.); // WX:08.2011 + for (i = 0; i < 3; i++) + tmpvalue += dFtdprin_s[i] * (prin_str[i] - sigaR[i]); + for (i = 0; i < 3; i++) + { + prin_str0[i] = prin_str[i]; + prin_str[i] -= tmpvalue * rtp[i]; + dStressP[i] = prin_str[i] - prin_str0[i]; + } + Inv_De->multi(dStressP, dStrainP); + + for (i = 0; i < 3; i++) + for (j = 0; j < 3; j++) + (*dGds_dFds)(i, j) = rtp[i] * dFtdprin_s[j]; + dGds_dFds->multi(*TmpDe, *TmpMatrix2); + // dGds_dFds->multi(*PrinDe, *tmpMatrix2); + *TmpDe -= *TmpMatrix2; + //*PrinDe -= *tmpMatrix2; + for (i = 0; i < 3; i++) + for (j = 0; j < 3; j++) + if (abs((*TmpDe)(i, j)) < MKleinsteZahl) + (*TmpDe)(i, j) = 1.; //(*Dep)(i,j)*1e-9;//avoid "0" in Dep + // for(i=3; iWrite(); + //*TmpDe = *Dep; //to be improved + // delete tmpMatrix2; + // delete dGds_dFds; + } + + else if (P_85 >= 0 && P_78 <= 0 && P_98 <= 0) // return to l1R + { + for (i = 0; i < 3; i++) + { + prin_str0[i] = prin_str[i]; + prin_str[i] = t1r1 * l1R[i] + sig1R[i]; + dStressP[i] = prin_str[i] - prin_str0[i]; + } + Inv_De->multi(dStressP, dStrainP); + VecCrossProduct(dStrainP, l1R, cVec); + CalDep_l(l1R, l1R, Inv_De, Dep_l, 1.0); + CalDep_l(cVec, cVec, Inv_De, dDep_l, fkt2); + //*TmpDe=(0.); + for (i = 0; i < 3; i++) + for (j = 0; j < 3; j++) + (*TmpDe)(i, j) = (*Dep_l)(i, j) + (*dDep_l)(i, j); + // for(i=3; i=0 && t1ra>=0 ) //return to sigRa + else if (P_98 >= 0) // return to sigRa + { + double tmpvalue = 0.; + // Matrix *dGds_dFds = new Matrix (Size,Size); + *dGds_dFds = (0.); // WX:08.2011 + + for (i = 0; i < 3; i++) + { + prin_str0[i] = prin_str[i]; + prin_str[i] = sigaR[i]; + dStressP[i] = prin_str[i] - prin_str0[i]; + } + Inv_De->multi(dStressP, dStrainP); + for (i = 0; i < 3; i++) + for (j = 0; j < 3; j++) + (*dGds_dFds)(i, j) = dStressP[i] * dStressP[j]; + tmpvalue = 0.; + for (i = 0; i < 3; i++) + tmpvalue += dStrainP[i] * dStressP[i]; + *dGds_dFds /= tmpvalue; + *TmpDe -= *dGds_dFds; + //*PrinDe -= *dGds_dFds; + //*TmpDe *= fkt1; + for (i = 0; i < Size; i++) + for (j = 0; j < Size; j++) + (*TmpDe)(i, j) *= fkt1; + //*PrinDe *= fkt1; + + // for(i=3; i=0) //return to sig1R + else if (P_52 >= 0 && P_45 <= 0 && P_85 <= 0) // return to sig1R + { + double tmpvalue = 0.; + // Matrix *dGds_dFds = new Matrix (Size,Size); + *dGds_dFds = (0.); // WX:08.2011 + + for (i = 0; i < 3; i++) + { + prin_str0[i] = prin_str[i]; + prin_str[i] = sig1R[i]; + dStressP[i] = prin_str[i] - prin_str0[i]; + } + Inv_De->multi(dStressP, dStrainP); + for (i = 0; i < 3; i++) + for (j = 0; j < 3; j++) + (*dGds_dFds)(i, j) = dStressP[i] * dStressP[j]; + tmpvalue = 0.; + for (i = 0; i < 3; i++) + tmpvalue += dStrainP[i] * dStressP[i]; + *dGds_dFds /= tmpvalue; + *TmpDe -= *dGds_dFds; + //*PrinDe -= *dGds_dFds; + //*TmpDe *= fkt1; + for (i = 0; i < Size; i++) + for (j = 0; j < Size; j++) + (*TmpDe)(i, j) *= fkt1; + + } + + else if (P_63 >= 0 && P_64 >= 0) // return to sig2R + { + double tmpvalue = 0.; + // Matrix *dGds_dFds = new Matrix (Size,Size); + *dGds_dFds = (0.); // WX:08.2011 + + for (i = 0; i < 3; i++) + { + prin_str0[i] = prin_str[i]; + prin_str[i] = sig2R[i]; + dStressP[i] = prin_str[i] - prin_str0[i]; + } + Inv_De->multi(dStressP, dStrainP); + for (i = 0; i < 3; i++) + for (j = 0; j < 3; j++) + (*dGds_dFds)(i, j) = dStressP[i] * dStressP[j]; + tmpvalue = 0.; + for (i = 0; i < 3; i++) + tmpvalue += dStrainP[i] * dStressP[i]; + *dGds_dFds /= tmpvalue; + *TmpDe -= *dGds_dFds; + //*PrinDe -= *dGds_dFds; + //*TmpDe *= fkt1; + for (i = 0; i < Size; i++) + for (j = 0; j < Size; j++) + (*TmpDe)(i, j) *= fkt1; + + } + else + { + cout << "error" << prin_str[0] << "||" << prin_str[1] << "||" << prin_str[2] << endl; + cout << "error" + << "csn=" << csn << "||" + << "Ntheta=" << Ntheta << "||" + << "tension=" << tension << endl; + } + // update prin str to normal coordinate + for (i = 0; i < Size; i++) + TryStress[i] = 0.; + *ConstitutiveMatrix = (0.); + + TransMatrixA_T->multi(prin_str, TryStress); // updated stress + TransMatrixA_T->multi(*TmpDe, *TransMatrixA, *ConstitutiveMatrix); // updated Depc + + for (i = 3; i < Size; i++) + (*ConstitutiveMatrix)(i, i) = (*Dep)(i, i); + + for (i = 0; i < Size; i++) + for (j = 3; j < Size; j++) + { + if (i < j) + { + (*ConstitutiveMatrix)(i, j) = 0.; + (*ConstitutiveMatrix)(j, i) = 0.; + } + } + + // mm = 0;//WX:12.2011 corecte Dep + + for (i = 0; i < Size; i++) + for (j = 0; j < Size; j++) + (*ConstitutiveMatrix)(i, j) = mm * (*Dep)(i, j) + (1 - mm) * (*ConstitutiveMatrix)(i, j); + ep += sqrt(2.0 / 3.0 * (TensorMutiplication2(dStrainP, dStrainP, Dim))); // updated eff plas strain + // ConstitutiveMatrix->Write(); + // TransMatrixA->Write(); + // TransMatrixA_T->Write(); + // if (mmHM_load_Factor) + // ele_val->HM_load_Factor = mm;//WX subincrement. not sure + } // end check of Mohr-Copulomb rupture in the matrix + else + { // calculate shear or tension in the joint + double normdstr = 0.; + for (i = 0; i < Size; i++) + { + //dstrs[i] = TryStress[i]; // d_stress + TryStress[i] = (*ele_val->Stress)(i, GPiGPj); // stress_0 + TryStr_buff[i] = TryStress[i] + dstrs[i]; + //TryStress_0[i] = TryStr_buff[i]; + normdstr += fabs(dstrs[i]); + } + /*if (normdstr < MKleinsteZahl) + { + //std::cout << "normdstr smaller than MKleinsteZahl" <multi(TryStr_buff, TmpStress); + TransMicroStru_TInv->multi(dstrs, dstrs_local); + TransMicroStru_TInv->multi(TryStress, TryStress_local); + for (i = 0; i < Size; i++) + { + TryStr_buff[i] = TmpStress[i]; // all stress in local cordinate sys. + } + *TmpMatrix = (0.); + *TmpMatrix2 = (0.); + CalculateCoefficent_MOHRjoint(ep, 0, 0); + if (GPiGPj == 5) + { + std::cout << "Princ, 1 (vert) -5.50e+06 , 2 (along strike) -3.70e+06 , stress 3 (perp to strike) -2.80e+06 for elem " << (*ele_val->GetIndex(), ele_value_id) << endl; + std::cout << "OGS , stress_xx " << (*ele_val->Stress)(0, GPiGPj) << " , stress_yy " << (*ele_val->Stress)(1, GPiGPj) << " , stress_zz " << (*ele_val->Stress)(2, GPiGPj) << " , stress_xy " << (*ele_val->Stress)(3, GPiGPj) << " , stress_xz " << (*ele_val->Stress)(4, GPiGPj) << " , stress_yz " << (*ele_val->Stress)(5, GPiGPj) << endl; + std::cout << "stress local_11 " << TryStr_buff[0] << " , stress local_22 " << TryStr_buff[1] << " , stress local_33 " << TryStr_buff[2] << " SHEAR local_12 " << TryStr_buff[3] << " , local_13 " << TryStr_buff[4] << " , local_23 " << TryStr_buff[5] << endl; + shearsurf = sqrt((TryStr_buff[5] * TryStr_buff[5]) + (TryStr_buff[4] * TryStr_buff[4])) + TryStr_buff[2] * tan(thetaj) - Y0j; + } + tensionsurf = -1; //TryStr_buff[2]; + if (tensionsurf > 1e-5) + { + //std::cout << "Shearsurf on joint" << shearsurf << endl; + shearsurf = -1; //sqrt((TryStr_buff[5] * TryStr_buff[5]) + (TryStr_buff[4] * TryStr_buff[4])); + //std::cout << "New shearsurf on joint" << shearsurf << endl; + //id_fs_n[2]= 0.0; + //id_fs_s[5] =-1 ; + //id_fs_s[4] =-1 ; + } + double shearstress = sqrt((TryStr_buff[5] * TryStr_buff[5]) + (TryStr_buff[4] * TryStr_buff[4])); + if (fabs(shearstress < 150000 ) || ep > 0.0035) + { + shearsurf = -1; + } + //std::cout << "Shear stress 3,4,5 " << shearstress << " Normal Stress " << tensionsurf < 1e-5 ) // || tensionsurf > 1e-5) // start check of Mohr-Copulomb rupture in the joint + { + yield = 1; + //std::cout << " Stating plastic joint calculation for element with GPiGPj " << GPiGPj << endl; + //if (tensionsurf > 1e-5) + //{ + //double Ei_new = (*data_Youngs)(0); + //Ei_new = Ei*1e3; + //(*data_Youngs)(0) = Ei_new; + //return 0; + //} + //else + //{ + //std::cout << "Shearsurf on joint" << shearsurf << endl; + + + //std::cout << "First shearsurf, shear stresses " << TryStr_buff[4] << " , " << TryStr_buff[5] << " shear strength " << TryStr_buff[2] * tan((*data_Plasticity_joint)(1) * PI / 180) - Y0j << " normal comp " << TryStr_buff[2] << endl; + double shear_stress_module = 0.0; + double dsig_dlamda[6] = { 0. }, /*dsig_dlamda_k1[6]={0.}, dft_dsig[6] = {0.},*/ dfs_dsig[6] = { 0. }, + dgs_dsig[6] = { 0. }, shear_cos_direction[6] = { 0. }; + bool first_step; + double shearsurf_k1/*, tensionsurf_k1 , local_damp, tmp_pos = 1., tmp_neg=1.*/; // Newton downhill + double lamda_pos = 0., lamda_neg = 0., shearsurf_pos = 0, shearsurf_neg = 0; + + if (shearsurf > 1e-5) // LU check if it works or should be larger ... + { + dlamda = 0; //, dlamda_1=1, dlamda_0=0.; + counter = 0; + first_step = true; + // local_damp = 1.; + shearsurf_k1 = shearsurf; + for (i = 0; i < Size; i++) + { + dstrs_local[i] = TmpStress[i] - TryStress_local[i]; + TryStr_buff[i] = TmpStress[i]; + } + + + while (1) + { + shear_stress_module = sqrt((TmpStress[5] * TmpStress[5]) + (TmpStress[4] * TmpStress[4]));//+(TmpStress[3]*TmpStress[3])); + + for (i = 0; i < Size; i++) + shear_cos_direction[i] = id_fs_s[i] * TmpStress[i] / shear_stress_module; + + + if (first_step) + { + for (i = 0; i < Size; i++) + { + dfs_dsig[i] = id_fs_n[i] * tan(thetaj) + shear_cos_direction[i]; + dgs_dsig[i] = id_fs_n[i] * tan(phij) + shear_cos_direction[i];; + } + } + // + if (counter == 0) + { + //std::cout << " Counter is still 0, first_step " << first_step << " shearsurf_k1 is "<< shearsurf_k1 < 0) + { + lamda_pos = dlamda; + TmpValue1 = 0.; + if (first_step) + { + for (i = 0; i < Size; i++) + dsig_dlamda[i] = 0.; + Dep->multi(dgs_dsig, dsig_dlamda, -1); + lamda_pos = 0.; + shearsurf_pos = shearsurf; + if (fabs(shearsurf_k1) < 1e-3) + { + std::cout << " At first step, ssurf < 1e-3 stop " << endl; + break; + } + // first_step = false; + } + for (i = 0; i < Size; i++) + { + TmpValue1 += dfs_dsig[i] * (-1) * dsig_dlamda[i]; //(dF/dstr)T De dG/dstr + // TmpValue2 += dfs_dsig[i]*dstrs[i]; + } + if (TmpValue1 == 0) + { + dlamda = 0.; + //std::cout << " dlamda is 0 " <multi(dgs_dsig,dsig_dlamda,-1);//update De dG/dsig + //std::cout << " Yes " << endl; + break; + } + else if (shearsurf_k1 < 0) + { + for(i=0;imulti(dgs_dsig,dsig_dlamda,-1);//update De dG/dsig, reactivated 13.09.2019 + shearsurf_neg = shearsurf_k1; + lamda_neg = dlamda; + dlamda = lamda_pos - shearsurf_pos * (lamda_neg - lamda_pos) / (shearsurf_neg - shearsurf_pos); + //std::cout << "No, shearsurf_k1 negative, dlamda lamda_neg-lamda_pos" << (lamda_neg - lamda_pos) <multi(dgs_dsig,dsig_dlamda,-1);//update De dG/dsig, reactivated 13.09.2019 + lamda_pos = dlamda; + shearsurf_pos = shearsurf_k1; + // LU: test to avoid mohr failure crashes 09.2019 + + if (shearsurf_neg > -1.e-13) + dlamda = lamda_pos + shearsurf_pos * (lamda_neg - lamda_pos) / (shearsurf_pos); + else + dlamda = lamda_pos - shearsurf_pos * (lamda_neg - lamda_pos) / (shearsurf_neg - shearsurf_pos); + // LU: test to avoid mohr failure crashes 09.2019 + //dlamda = lamda_pos - shearsurf_pos * (lamda_neg - lamda_pos) / (shearsurf_neg - shearsurf_pos); + //std::cout << "No, shearsurf_k1 positive, dlamda lamda_neg-lamda_pos" << (lamda_neg - lamda_pos) <GetTranspose(*TmpMatrix2); for (i = 0; i < Size; i++) { TmpStrComp[i] = 0.; TmpStrTens[i] = 0.; if (tmp_prin_str[i] < MKleinsteZahl) TmpPrinStrComp[i] = tmp_prin_str[i]; else TmpPrinStrTens[i] = tmp_prin_str[i]; } TmpMatrix2->multi(TmpPrinStrComp, TmpStrComp); TmpMatrix2->multi(TmpPrinStrTens, TmpStrTens); AnisoParaComp = CalAnisoPara(TmpStrComp, MicroStruTensor); AnisoParaTens = CalAnisoPara(TmpStrTens, MicroStruTensor); //CalculateCoefficent_MOHR(ep, AnisoParaComp, AnisoParaTens); CalculateCoefficent_MOHR(ep, 0, 0); //ignore anisotropy LU 06.2019 + } + + */ + shear_stress_module = sqrt((TmpStress[5] * TmpStress[5]) + (TmpStress[4] * TmpStress[4]));//+(TmpStress[3]*TmpStress[3])); + for (i = 0; i < Size; i++) + shear_cos_direction[i] = id_fs_s[i] * TmpStress[i] / shear_stress_module; + + shearsurf_k1 = sqrt((TmpStress[5] * TmpStress[5]) + (TmpStress[4] * TmpStress[4])) + TmpStress[2] * tan(thetaj) - Y0j; + std::cout << "Iteration " << counter << " Ssurf_k1 on joint " << shearsurf_k1 << " Lambda pos " << lamda_pos << " ssurf pos " << shearsurf_pos << " Lambda neg " << lamda_neg << " ssurf neg " << shearsurf_neg << endl; + + // std::cout << "Lambda pos" << lamda_pos << "ssurf pos" << shearsurf_pos << endl; + // std::cout << "Lambda neg" << lamda_neg << "ssurf neg" << shearsurf_neg << endl; + if (tensionsurf > 1e-5) + { + shearsurf_k1 = -1;//sqrt((TmpStress[5] * TmpStress[5]) + (TmpStress[4] * TmpStress[4])); + std::cout << "But tensile surf is positive, New shearsurf_k1 on joint" << shearsurf_k1 << endl; + + if (fabs(shearsurf_k1) < 1e-6) + break; + } + //tensionsurf = -1.e5; + //std::cout << i <<"Freshly baked shearsurf_k1, shear stresses " << TmpStress[4] << " , " << TmpStress[5] <<" shear strength " << TmpStress[2] * tan((*data_Plasticity)(1) * PI / 180) - Y0 <<" normal comp" <Write(); + // TmpMatrix2->Write(); + // + for (i = 0; i < Size; i++) + for (j = 0; j < Size; j++) + (*TmpMatrix)(i, j) = dgs_dsig[i] * dfs_dsig[j]; + Dep->multi(*TmpMatrix, *Dep, *TmpMatrix2); // De (dG/dstr)T dF/dstr De + // + // TmpDe->Write(); + // TmpMatrix2->Write(); + // + for (i = 0; i < Size; i++) + { + for (j = 0; j < Size; j++) + (*TmpDe)(i, j) -= (*TmpMatrix2)(i, j) / TmpValue1; // Dep = De - De dG/ds dF/ds De / dF/ds de dG/ds + } + // TEST + // TransMicroStru->Write(); + // TransMicroStru_T->Write(); + // TransMicroStru_TInv->Write(); + // TmpDe->Write(); + // Dep->Write(); + // + } // end if shear + + for (i = 0; i < Size; i++) + TryStress[i] = 0; + TransMicroStru_T->multi(TmpStress, TryStress); + //std::cout << "Final shearsurf, shear stresses " << TmpStress[4] << " , " << TmpStress[5] << " shear strength " << TmpStress[2] * tan((*data_Plasticity_joint)(1) * PI / 180) - Y0j << " normal comp " << TmpStress[2] << endl; + for (i = 0; i < Size; i++) //LU output to check stress values 06.19 + { + + //std::cout << i <<"-component old stress " <<(*ele_val->Stress)(i, GPiGPj) << ",new stress " << TryStress[i] <multi(*TmpDe, *TransMicroStru, *ConstitutiveMatrix); + // TEST + // ConstitutiveMatrix->Write(); + // + Cal_Inv_Matrix(Size, Dep, Inv_De); + for (i = 0; i < Size; i++) + dstressP[i] = TryStress[i] - TryStress_0[i]; + Inv_De->multi(dstressP, dstrainP); + ep += sqrt(2.0 / 3.0 * (TensorMutiplication2(dstrainP, dstrainP, Dim))); + + } + //} + else // no rupture in the joint or in the matrix => update stress simply + { + for (i = 0; i < Size; i++) + TryStress[i] += dstrs[i]; + } + } + + if (Update) + { + (*ele_val->pStrain)(GPiGPj) = ep; + } + return yield; +} /******************************************************* WX: Mohr coulomb, directe stress integration, also for aniso.