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

[NEW FIX 2] Fixes to TDMN tangents. Removal of older convention of gradient... 

[NEW FIX 2] Fixes to TDMN tangents. Removal of older convention of gradient vector stack. Bug patches to follow.
parent cba9325b
No related branches found
No related tags found
1 merge request!430Added 1 modified class for tempGrad derivatives in Network Interactions,...
Show whitespace changes
Inline Side-by-side
dG3D/src/dG3DMaterialLaw.cpp
+ 6
− 38
View file @ 440f8ac4
...@@ -5671,10 +5671,8 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, ...@@ -5671,10 +5671,8 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
else if(!_flag_isothermal && _flag_microTempFixed){ else if(!_flag_isothermal && _flag_microTempFixed){
Tloc = T; Tloc = T;
for(int m=0; m<3; m++){ for(int m=0; m<3; m++){
// Hloc(m) = H(m) + FHvec(i*_col_Di+ (m+1)*_row_Di-1);
Hloc(m) = H(m) + FHvec(i*_col_Di+ 9+m); Hloc(m) = H(m) + FHvec(i*_col_Di+ 9+m);
for(int n=0; n<3; n++){ for(int n=0; n<3; n++){
// Floc(m,n) = F(m,n) + FHvec(i*_col_Di + m + n*_row_Di);
Floc(m,n) = F(m,n) + FHvec(i*_col_Di + m + n*3); Floc(m,n) = F(m,n) + FHvec(i*_col_Di + m + n*3);
} }
} }
...@@ -5682,10 +5680,8 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, ...@@ -5682,10 +5680,8 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
else if(!_flag_isothermal && !_flag_microTempFixed){ else if(!_flag_isothermal && !_flag_microTempFixed){
Tloc = T + FHvec((i+1)*_col_Di-1); Tloc = T + FHvec((i+1)*_col_Di-1);
for(int m=0; m<3; m++){ for(int m=0; m<3; m++){
// Hloc(m) = H(m) + FHvec(i*_col_Di+ (m+1)*(_row_Di-1)-1);
Hloc(m) = H(m) + FHvec(i*_col_Di+ 9+m); Hloc(m) = H(m) + FHvec(i*_col_Di+ 9+m);
for(int n=0; n<3; n++){ for(int n=0; n<3; n++){
// Floc(m,n) = F(m,n) + FHvec(i*_col_Di + m + n*(_row_Di-1));
Floc(m,n) = F(m,n) + FHvec(i*_col_Di + m + n*3); Floc(m,n) = F(m,n) + FHvec(i*_col_Di + m + n*3);
} }
} }
...@@ -5721,10 +5717,8 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, ...@@ -5721,10 +5717,8 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
else if(!_flag_isothermal && _flag_microTempFixed){ else if(!_flag_isothermal && _flag_microTempFixed){
Cpvec(i) = Cp_i; thermSrcvec(i) = w_i; mechSrcvec(i) = mechSource_i; Cpvec(i) = Cp_i; thermSrcvec(i) = w_i; mechSrcvec(i) = mechSource_i;
for(int m=0; m<3; m++){ for(int m=0; m<3; m++){
// PQvec(i*_col_Di + (m+1)*_row_Di-1) = Q_i(m);
PQvec(i*_col_Di+ 9+m) = Q_i(m); PQvec(i*_col_Di+ 9+m) = Q_i(m);
for(int n=0; n<3; n++){ for(int n=0; n<3; n++){
// PQvec(i*_col_Di + m + n*_row_Di) = P_i(m,n);
PQvec(i*_col_Di + m + n*3) = P_i(m,n); PQvec(i*_col_Di + m + n*3) = P_i(m,n);
} }
} }
...@@ -5733,10 +5727,8 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, ...@@ -5733,10 +5727,8 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
thermSrcvec(i) = w_i; mechSrcvec(i) = mechSource_i; thermSrcvec(i) = w_i; mechSrcvec(i) = mechSource_i;
PQvec((i+1)*_col_Di-1) = Cp_i; PQvec((i+1)*_col_Di-1) = Cp_i;
for(int m=0; m<3; m++){ for(int m=0; m<3; m++){
// PQvec(i*_col_Di + (m+1)*(_row_Di-1)-1) = Q_i(m);
PQvec(i*_col_Di+ 9+m) = Q_i(m); PQvec(i*_col_Di+ 9+m) = Q_i(m);
for(int n=0; n<3; n++){ for(int n=0; n<3; n++){
// PQvec(i*_col_Di + m + n*(_row_Di-1)) = P_i(m,n);
PQvec(i*_col_Di + m + n*3) = P_i(m,n); PQvec(i*_col_Di + m + n*3) = P_i(m,n);
} }
} }
...@@ -5772,16 +5764,13 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, ...@@ -5772,16 +5764,13 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
dwdtvec(i) = dwdt_i; dmechSrcdTvec(i) = dmechSourcedt_i; 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 + 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); dwdFvec(i*9+m+_row_Di*n) = dwdf_i(m,n);
dmechSrcdFvec(i*9+m+_row_Di*n) = dmechSourcedf_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+_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
_C(i*_col_Di + 9+m, i*_col_Di+n+3*k) = dqdF_i(m,n,k); _C(i*_col_Di + 9+m, i*_col_Di+n+3*k) = dqdF_i(m,n,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+3*n, i*_col_Di+k+3*l) = L_i(m,n,k,l); _C(i*_col_Di+m+3*n, i*_col_Di+k+3*l) = L_i(m,n,k,l);
} }
} }
...@@ -5791,18 +5780,14 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, ...@@ -5791,18 +5780,14 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
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++){
// _C(i*_col_Di + (m+1)*(_row_Di-1)-1, i*_col_Di + (n+1)*(_row_Di-1)-1) = dqdgradT_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); dwdFvec(i*9+m+_row_Di*n) = dwdf_i(m,n);
dmechSrcdFvec(i*9+m+_row_Di*n) = dmechSourcedf_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+(_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
_C(i*_col_Di + 9+m, i*_col_Di+n+3*k) = dqdF_i(m,n,k); _C(i*_col_Di + 9+m, i*_col_Di+n+3*k) = dqdF_i(m,n,k);
for(int l=0; l<3; l++) for(int l=0; l<3; l++)
// _C(i*_col_Di+m+(_row_Di-1)*n, i*_col_Di+k+(_row_Di-1)*l) = L_i(m,n,k,l);
_C(i*_col_Di+m+3*n, i*_col_Di+k+3*l) = L_i(m,n,k,l); _C(i*_col_Di+m+3*n, i*_col_Di+k+3*l) = L_i(m,n,k,l);
} }
} }
...@@ -5836,7 +5821,7 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, ...@@ -5836,7 +5821,7 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
Jacobian.invertInPlace(); Jacobian.invertInPlace();
} }
else{ else{
// What is this? How does the following update the jacobian? // Updating Jacobian
_NTV.mult(PQvec, Res_node); _NTV.mult(PQvec, Res_node);
Res_node_t.scale(-1.0); Res_node_t.scale(-1.0);
Res_node_t.axpy(Res_node, 1.0); Res_node_t.axpy(Res_node, 1.0);
...@@ -5859,7 +5844,6 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, ...@@ -5859,7 +5844,6 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
} }
} }
stiff_loc = true; stiff_loc = true;
// What?
} }
   
Res_node_t = Res_node; Res_node_t = Res_node;
...@@ -5895,17 +5879,15 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, ...@@ -5895,17 +5879,15 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
for(int i=0;i<3;i++){ for(int i=0;i<3;i++){
Q(i) = 0.; Q(i) = 0.;
for(int j=0;j<3;j++){ for(int j=0;j<3;j++){
P(i,j) = PQ_hom[i+_row_Di*j]; // CHECK P(i,j) = PQ_hom[i+_row_Di*j];
} }
} }
} }
else if(!_flag_isothermal && _flag_microTempFixed){ else if(!_flag_isothermal && _flag_microTempFixed){
Cp = Cp_hom; w = thermSrc_hom; mechSource = mechSource_hom; Cp = Cp_hom; w = thermSrc_hom; mechSource = mechSource_hom;
for(int i=0;i<3;i++){ for(int i=0;i<3;i++){
// Q(i) = PQ_hom[(i+1)*_row_Di-1];
Q(i) = PQ_hom[9+i]; Q(i) = PQ_hom[9+i];
for(int j=0;j<3;j++){ for(int j=0;j<3;j++){
// P(i,j) = PQ_hom[i+_row_Di*j];
P(i,j) = PQ_hom[i+3*j]; P(i,j) = PQ_hom[i+3*j];
} }
} }
...@@ -5913,10 +5895,8 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, ...@@ -5913,10 +5895,8 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
else if(!_flag_isothermal && !_flag_microTempFixed){ else if(!_flag_isothermal && !_flag_microTempFixed){
Cp = PQ_hom[_col_Di-1]; w = thermSrc_hom; mechSource = mechSource_hom; Cp = PQ_hom[_col_Di-1]; w = thermSrc_hom; mechSource = mechSource_hom;
for(int i=0;i<3;i++){ for(int i=0;i<3;i++){
// Q(i) = PQ_hom[(i+1)*(_row_Di-1)-1];
Q(i) = PQ_hom[9+i]; Q(i) = PQ_hom[9+i];
for(int j=0;j<3;j++){ for(int j=0;j<3;j++){
// P(i,j) = PQ_hom[i+(_row_Di-1)*j];
P(i,j) = PQ_hom[i+3*j]; P(i,j) = PQ_hom[i+3*j];
} }
} }
...@@ -5940,15 +5920,9 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, ...@@ -5940,15 +5920,9 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
} }
} }
   
double temp(0.), dTidT(0.); double temp(0.), dTidT(1.); // = Cp_hom/(temp*Cpvec(i))
for(int i=0;i<_NRoot;i++){ for(int i=0;i<_NRoot;i++){
temp = (_VA[1]*_V(0, _col_Di*i) + _VA[2]*_V(_col_Di, _col_Di*i)); // Check temp = (_VA[1]*_V(0, _col_Di*i) + _VA[2]*_V(_col_Di, _col_Di*i));
if (abs(temp) > 0. && Cpvec(i) > 0.){
dTidT = Cp_hom/(temp*Cpvec(i));
}
else{
dTidT = 0.;
}
dwdt_hom += temp*dwdtvec(i)*dTidT; // dwdT average dwdt_hom += temp*dwdtvec(i)*dTidT; // dwdT average
dmechSourcedt_hom += temp*dmechSrcdTvec(i)*dTidT; // dmechSourcedt average dmechSourcedt_hom += temp*dmechSrcdTvec(i)*dTidT; // dmechSourcedt average
for(int m=0;m<9;m++) for(int m=0;m<9;m++)
...@@ -5969,16 +5943,13 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, ...@@ -5969,16 +5943,13 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
dwdt = dwdt_hom; dmechSourcedt = dmechSourcedt_hom; 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[9+m][9+n]; dqdgradT(m,n) = C_hom[9+m][9+n];
dwdf(m,n) = dwdf_hom[m+3*n]; dwdf(m,n) = dwdf_hom[m+3*n];
dmechSourcedf(m,n) = dmechSourcedf_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[9+m][n+3*k]; dqdF(m,n,k) = C_hom[9+m][n+3*k];
// dPdgradT(m,n,l) = C_hom[m+_row_Di*n][(k+1)*_row_Di-1]; // 0. This is for dPdgradT // dPdgradT(m,n,l) = C_hom[m+_row_Di*n][(k+1)*_row_Di-1]; // 0. This is for dPdgradT
for(int l=0; l<3; l++) for(int l=0; l<3; l++)
// L(m,n,k,l) = C_hom[m+_row_Di*n][k+_row_Di*l];
L(m,n,k,l) = C_hom[m+3*n][k+3*l]; L(m,n,k,l) = C_hom[m+3*n][k+3*l];
} }
} }
...@@ -5987,17 +5958,14 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, ...@@ -5987,17 +5958,14 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
for(int m=0; m<3; m++){ for(int m=0; m<3; m++){
dqdT(m) = C_hom[(m+1)*(_row_Di-1)-1][_col_Di-1]; dqdT(m) = C_hom[(m+1)*(_row_Di-1)-1][_col_Di-1];
for(int n=0; n<3; n++){ for(int n=0; n<3; n++){
// 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]; dwdf(m,n) = dwdf_hom[m+3*n];
dmechSourcedf(m,n) = dmechSourcedf_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[9+m][n+3*k]; dqdF(m,n,k) = C_hom[9+m][n+3*k];
// C_hom[m+(_row_Di-1)*n][(k+1)*(_row_Di-1)-1] = 0.; This is for dPdgradT // C_hom[m+(_row_Di-1)*n][(k+1)*(_row_Di-1)-1] = 0.; This is for dPdgradT
for(int l=0; l<3; l++) for(int l=0; l<3; l++)
// L(m,n,k,l) = C_hom[m+(_row_Di-1)*n][k+(_row_Di-1)*l];
L(m,n,k,l) = C_hom[m+3*n][k+3*l]; L(m,n,k,l) = C_hom[m+3*n][k+3*l];
} }
} }
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment