Skip to content
Snippets Groups Projects
Commit cba9325b authored by Ujwal Kishore Jinaga's avatar Ujwal Kishore Jinaga :clown:
Browse files

[NEW FIX 1] TMDMN tangents Fix 1. More to follow.

parent b4610983
No related branches found
No related tags found
1 merge request!430Added 1 modified class for tempGrad derivatives in Network Interactions,...
...@@ -5583,9 +5583,9 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, ...@@ -5583,9 +5583,9 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
   
STensor3& P = ipvcur->getRefToFirstPiolaKirchhoffStress(); STensorOperation::zero(P); STensor3& P = ipvcur->getRefToFirstPiolaKirchhoffStress(); STensorOperation::zero(P);
SVector3& Q = ipvcur->getRefToThermalFlux(); STensorOperation::zero(Q); SVector3& Q = ipvcur->getRefToThermalFlux(); STensorOperation::zero(Q);
double& w = ipvcur->getRefToThermalSource(); double& w = ipvcur->getRefToThermalSource(); w = 0.;
double& Cp = ipvcur->getRefToExtraDofFieldCapacityPerUnitField()(0); double& Cp = ipvcur->getRefToExtraDofFieldCapacityPerUnitField()(0);
double& mechSource = ipvcur->getRefToMechanicalSource()(0); double& mechSource = ipvcur->getRefToMechanicalSource()(0); mechSource = 0.;
   
STensor43& L = ipvcur->getRefToTangentModuli(); STensorOperation::zero(L); STensor43& L = ipvcur->getRefToTangentModuli(); STensorOperation::zero(L);
STensor3& dPdT = ipvcur->getRefTodPdT(); STensorOperation::zero(dPdT); STensor3& dPdT = ipvcur->getRefTodPdT(); STensorOperation::zero(dPdT);
...@@ -5594,7 +5594,7 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, ...@@ -5594,7 +5594,7 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
STensor33& dqdF = ipvcur->getRefTodThermalFluxdF(); STensorOperation::zero(dqdF); STensor33& dqdF = ipvcur->getRefTodThermalFluxdF(); STensorOperation::zero(dqdF);
double& dwdt = ipvcur->getRefTodThermalSourcedField(); double& dwdt = ipvcur->getRefTodThermalSourcedField();
STensor3& dwdf = ipvcur->getRefTodThermalSourcedF(); STensorOperation::zero(dwdf); STensor3& dwdf = ipvcur->getRefTodThermalSourcedF(); STensorOperation::zero(dwdf);
double& dmechSourcedt = ipvcur->getRefTodMechanicalSourcedField()(0,0); double& dmechSourcedt = ipvcur->getRefTodMechanicalSourcedField()(0,0); dmechSourcedt = 0.;
STensor3& dmechSourcedf = ipvcur->getRefTodMechanicalSourcedF()[0]; STensorOperation::zero(dmechSourcedf); STensor3& dmechSourcedf = ipvcur->getRefTodMechanicalSourcedF()[0]; STensorOperation::zero(dmechSourcedf);
double dCpdT(0.); double dCpdT(0.);
STensor3 dCpdF; STensorOperation::zero(dCpdF); STensor3 dCpdF; STensorOperation::zero(dCpdF);
...@@ -5603,7 +5603,9 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, ...@@ -5603,7 +5603,9 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
   
double C_hom[_col_Di][_col_Di]; // static // The tangent is combined with thermal tangents. double C_hom[_col_Di][_col_Di]; // static // The tangent is combined with thermal tangents.
double PQ_hom[_col_Di]; //static double PQ_hom[_col_Di]; //static
double Cp_hom(0.), thermSrc_hom(0.), mechSource_hom(0.); double dPQdT_hom[_col_Di]; // static
double Cp_hom(0.), thermSrc_hom(0.), mechSource_hom(0.), dwdt_hom(0.), dmechSourcedt_hom(0.);
static double dwdf_hom[9], dmechSourcedf_hom[9];
// std::vector<double> PQ_hom; // std::vector<double> PQ_hom;
   
fullVector<double>& FHvec = ipvcur->getRefToCombinedGradientVect(); fullVector<double>& FHvec = ipvcur->getRefToCombinedGradientVect();
...@@ -5612,11 +5614,14 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, ...@@ -5612,11 +5614,14 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
fullVector<double>& Cpvec = ipvcur->getRefToCpVect(); fullVector<double>& Cpvec = ipvcur->getRefToCpVect();
fullVector<double>& thermSrcvec = ipvcur->getRefToThermalSourceVect(); fullVector<double>& thermSrcvec = ipvcur->getRefToThermalSourceVect();
fullVector<double>& mechSrcvec = ipvcur->getRefToMechSourceVect(); fullVector<double>& mechSrcvec = ipvcur->getRefToMechSourceVect();
fullVector<double> dwdtvec, dmechSrcdTvec, dwdFvec, dmechSrcdFvec;
dwdtvec.resize(_NRoot); dmechSrcdTvec.resize(_NRoot);
dwdFvec.resize(9*_NRoot); dmechSrcdFvec.resize(9*_NRoot);
   
PQvec.setAll(0.0); Cpvec.setAll(0.0); thermSrcvec.setAll(0.0); mechSrcvec.setAll(0.0); PQvec.setAll(0.0); Cpvec.setAll(0.0); thermSrcvec.setAll(0.0); mechSrcvec.setAll(0.0);
static fullVector<double> Res_node_t(_row_Di*_NInterface); static fullVector<double> Res_node_t(_row_Di*_NInterface);
static fullVector<double> Res_node(_row_Di*_NInterface); static fullVector<double> Res_node(_row_Di*_NInterface);
static fullVector<double> a_step(_row_Di*_NInterface); static fullVector<double> ab_step(_row_Di*_NInterface);
static fullMatrix<double> NTVC(_row_Di*_NInterface, _col_Di*_NRoot); static fullMatrix<double> NTVC(_row_Di*_NInterface, _col_Di*_NRoot);
static fullMatrix<double> Jacobian_inv(_row_Di*_NInterface, _row_Di*_NInterface); static fullMatrix<double> Jacobian_inv(_row_Di*_NInterface, _row_Di*_NInterface);
static fullMatrix<double> Modify_Jacobian(_row_Di*_NInterface, _row_Di*_NInterface); static fullMatrix<double> Modify_Jacobian(_row_Di*_NInterface, _row_Di*_NInterface);
...@@ -5753,17 +5758,24 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, ...@@ -5753,17 +5758,24 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
const STensor3& dmechSourcedf_i = ipv_i->getRefTodMechanicalSourcedF()[0]; const STensor3& dmechSourcedf_i = ipv_i->getRefTodMechanicalSourcedF()[0];
if (_flag_isothermal && _flag_microTempFixed){ if (_flag_isothermal && _flag_microTempFixed){
dwdtvec(i) = 0.; dmechSrcdTvec(i) = 0.;
for(int m=0; m<3; m++) for(int m=0; m<3; m++)
for(int n=0; n<3; n++) for(int n=0; n<3; n++){
dwdFvec(i*_col_Di+m+_row_Di*n) = 0.;
dmechSrcdFvec(i*_col_Di+m+_row_Di*n) = 0.;
for(int k=0; k<3; k++) for(int k=0; k<3; k++)
for(int l=0; l<3; l++) for(int l=0; l<3; l++)
_C(i*_col_Di+m+_row_Di*n, i*_col_Di+k+_row_Di*l) = L_i(m,n,k,l); _C(i*_col_Di+m+_row_Di*n, i*_col_Di+k+_row_Di*l) = L_i(m,n,k,l);
} }
}
else if(!_flag_isothermal && _flag_microTempFixed){ else if(!_flag_isothermal && _flag_microTempFixed){
dwdtvec(i) = dwdt_i; dmechSrcdTvec(i) = dmechSourcedt_i;
for(int m=0; m<3; m++) for(int m=0; m<3; m++)
for(int n=0; n<3; n++){ for(int n=0; n<3; n++){
// _C(i*_col_Di + (m+1)*_row_Di-1, i*_col_Di + (n+1)*_row_Di-1) = dqdgradT_i(m,n); // _C(i*_col_Di + (m+1)*_row_Di-1, i*_col_Di + (n+1)*_row_Di-1) = dqdgradT_i(m,n);
_C(i*_col_Di + 9+m, i*_col_Di + 9+n) = dqdgradT_i(m,n); _C(i*_col_Di + 9+m, i*_col_Di + 9+n) = dqdgradT_i(m,n);
dwdFvec(i*9+m+_row_Di*n) = dwdf_i(m,n);
dmechSrcdFvec(i*9+m+_row_Di*n) = dmechSourcedf_i(m,n);
for(int k=0; k<3; k++){ for(int k=0; k<3; k++){
// _C(i*_col_Di + (m+1)*_row_Di-1, i*_col_Di+n+_row_Di*k) = dqdF_i(m,n,k); // _C(i*_col_Di + (m+1)*_row_Di-1, i*_col_Di+n+_row_Di*k) = dqdF_i(m,n,k);
// _C(i*_col_Di+m+_row_Di*n, i*_col_Di + (k+1)*_row_Di-1) = 0.; // This is for dPdgradT_i // _C(i*_col_Di+m+_row_Di*n, i*_col_Di + (k+1)*_row_Di-1) = 0.; // This is for dPdgradT_i
...@@ -5775,6 +5787,7 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, ...@@ -5775,6 +5787,7 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
} }
} }
else if(!_flag_isothermal && !_flag_microTempFixed){ else if(!_flag_isothermal && !_flag_microTempFixed){
dwdtvec(i) = dwdt_i; dmechSrcdTvec(i) = dmechSourcedt_i;
for(int m=0; m<3; m++){ for(int m=0; m<3; m++){
_C(i*_col_Di + (m+1)*(_row_Di-1)-1, (i+1)*_col_Di-1) = dqdT_i(m); _C(i*_col_Di + (m+1)*(_row_Di-1)-1, (i+1)*_col_Di-1) = dqdT_i(m);
for(int n=0; n<3; n++){ for(int n=0; n<3; n++){
...@@ -5782,6 +5795,8 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, ...@@ -5782,6 +5795,8 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
// _C(i*_col_Di + m + n*(_row_Di-1), (i+1)*_col_Di-1) = dPdT_i(m,n); // _C(i*_col_Di + m + n*(_row_Di-1), (i+1)*_col_Di-1) = dPdT_i(m,n);
_C(i*_col_Di + 9+m, i*_col_Di + 9+m) = dqdgradT_i(m,n); _C(i*_col_Di + 9+m, i*_col_Di + 9+m) = dqdgradT_i(m,n);
_C(i*_col_Di + m + 3*n, (i+1)*_col_Di-1) = dPdT_i(m,n); _C(i*_col_Di + m + 3*n, (i+1)*_col_Di-1) = dPdT_i(m,n);
dwdFvec(i*9+m+_row_Di*n) = dwdf_i(m,n);
dmechSrcdFvec(i*9+m+_row_Di*n) = dmechSourcedf_i(m,n);
for(int k=0; k<3; k++){ for(int k=0; k<3; k++){
// _C(i*_col_Di + (m+1)*(_row_Di-1)-1, i*_col_Di+n+(_row_Di-1)*k) = dqdF_i(m,n,k); // _C(i*_col_Di + (m+1)*(_row_Di-1)-1, i*_col_Di+n+(_row_Di-1)*k) = dqdF_i(m,n,k);
// _C(i*_col_Di+m+(_row_Di-1)*n, i*_col_Di + (k+1)*(_row_Di-1)-1) = 0.; // This is for dPdgradT_i // _C(i*_col_Di+m+(_row_Di-1)*n, i*_col_Di + (k+1)*(_row_Di-1)-1) = 0.; // This is for dPdgradT_i
...@@ -5867,11 +5882,12 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, ...@@ -5867,11 +5882,12 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
} }
} }
Cp_hom, thermSrc_hom, mechSource_hom = 0.; double temp = 0.;
for(int i=0;i<_NRoot;i++){ for(int i=0;i<_NRoot;i++){
Cp_hom += (_VA[1]*_V(0, _col_Di*i) + _VA[2]*_V(_col_Di, _col_Di*i)) *Cpvec(i); // Cp average temp = (_VA[1]*_V(0, _col_Di*i) + _VA[2]*_V(_col_Di, _col_Di*i));
thermSrc_hom += (_VA[1]*_V(0, _col_Di*i) + _VA[2]*_V(_col_Di, _col_Di*i)) *thermSrcvec(i); // w average Cp_hom += temp*Cpvec(i); // Cp average
mechSource_hom += (_VA[1]*_V(0, _col_Di*i) + _VA[2]*_V(_col_Di, _col_Di*i)) *mechSrcvec(i); // mechSrc average thermSrc_hom += temp*thermSrcvec(i); // w average
mechSource_hom += temp*mechSrcvec(i); // mechSrc average
} }
   
if (_flag_isothermal && _flag_microTempFixed){ if (_flag_isothermal && _flag_microTempFixed){
...@@ -5911,17 +5927,35 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, ...@@ -5911,17 +5927,35 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
NTVC.mult(_I, dRdF); NTVC.mult(_I, dRdF);
Jacobian.mult(dRdF, dadF); Jacobian.mult(dRdF, dadF);
_Nv.mult(dadF, tmp); _Nv.mult(dadF, tmp);
tmp.add(_I); tmp.add(_I); // Check this structure of _I. tmp = dFidF
   
for(int m=0;m<_col_Di;m++){ for(int m=0;m<_col_Di;m++){
for(int n=0;n<_col_Di;n++){ for(int n=0;n<_col_Di;n++){
C_hom[m][n] = 0.0; C_hom[m][n] = 0.0;
for(int i=0;i<_NRoot;i++){ for(int i=0;i<_NRoot;i++){
for(int j=0;j< _col_Di*_NRoot;j++){ for(int j=0;j< _col_Di*_NRoot;j++){
C_hom[m][n] += (_VA[1]*_V(m,m+_col_Di*i) + _VA[2]*_V(_col_Di+m,m+_col_Di*i))*_C(m+_col_Di*i,j)*tmp(j,n); // CHECK C_hom[m][n] += (_VA[1]*_V(m,m+_col_Di*i) + _VA[2]*_V(_col_Di+m,m+_col_Di*i))*_C(m+_col_Di*i,j)*tmp(j,n);
}
} }
} }
} }
double temp(0.), dTidT(0.);
for(int i=0;i<_NRoot;i++){
temp = (_VA[1]*_V(0, _col_Di*i) + _VA[2]*_V(_col_Di, _col_Di*i)); // Check
if (abs(temp) > 0. && Cpvec(i) > 0.){
dTidT = Cp_hom/(temp*Cpvec(i));
}
else{
dTidT = 0.;
}
dwdt_hom += temp*dwdtvec(i)*dTidT; // dwdT average
dmechSourcedt_hom += temp*dmechSrcdTvec(i)*dTidT; // dmechSourcedt average
for(int m=0;m<9;m++)
for(int n=0;n<9;n++){
dwdf_hom[n] += temp*dwdFvec(m+9*i)*tmp(9*i,n);
dmechSourcedf_hom[n] += temp*dmechSrcdTvec(m+9*i)*tmp(9*i,n);
}
} }
   
if (_flag_isothermal && _flag_microTempFixed){ if (_flag_isothermal && _flag_microTempFixed){
...@@ -5932,10 +5966,13 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, ...@@ -5932,10 +5966,13 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
L(m,n,k,l) = C_hom[m+_row_Di*n][k+_row_Di*l]; L(m,n,k,l) = C_hom[m+_row_Di*n][k+_row_Di*l];
} }
else if(!_flag_isothermal && _flag_microTempFixed){ else if(!_flag_isothermal && _flag_microTempFixed){
dwdt = dwdt_hom; dmechSourcedt = dmechSourcedt_hom;
for(int m=0; m<3; m++) for(int m=0; m<3; m++)
for(int n=0; n<3; n++){ for(int n=0; n<3; n++){
// dqdgradT(m,n) = C_hom[(m+1)*_row_Di-1][(n+1)*_row_Di-1]; // dqdgradT(m,n) = C_hom[(m+1)*_row_Di-1][(n+1)*_row_Di-1];
dqdgradT(m,n) = C_hom[9+m][9+n]; dqdgradT(m,n) = C_hom[9+m][9+n];
dwdf(m,n) = dwdf_hom[m+3*n];
dmechSourcedf(m,n) = dmechSourcedf_hom[m+3*n];
for(int k=0; k<3; k++){ for(int k=0; k<3; k++){
// dqdF(m,n,k) = C_hom[(m+1)*_row_Di-1][n+_row_Di*k]; // dqdF(m,n,k) = C_hom[(m+1)*_row_Di-1][n+_row_Di*k];
dqdF(m,n,k) = C_hom[9+m][n+3*k]; dqdF(m,n,k) = C_hom[9+m][n+3*k];
...@@ -5953,6 +5990,8 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, ...@@ -5953,6 +5990,8 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
// dqdgradT(m,n) = C_hom[(m+1)*(_row_Di-1)-1][(n+1)*(_row_Di-1)-1]; // dqdgradT(m,n) = C_hom[(m+1)*(_row_Di-1)-1][(n+1)*(_row_Di-1)-1];
dPdT(m,n) = C_hom[m + n*(_row_Di-1)][_col_Di-1]; dPdT(m,n) = C_hom[m + n*(_row_Di-1)][_col_Di-1];
dqdgradT(m,n) = C_hom[9+m][9+n]; dqdgradT(m,n) = C_hom[9+m][9+n];
dwdf(m,n) = dwdf_hom[m+3*n];
dmechSourcedf(m,n) = dmechSourcedf_hom[m+3*n];
for(int k=0; k<3; k++){ for(int k=0; k<3; k++){
// dqdF(m,n,k) = C_hom[(m+1)*(_row_Di-1)-1][n+(_row_Di-1)*k]; // dqdF(m,n,k) = C_hom[(m+1)*(_row_Di-1)-1][n+(_row_Di-1)*k];
dqdF(m,n,k) = C_hom[9+m][n+3*k]; dqdF(m,n,k) = C_hom[9+m][n+3*k];
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment