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

[NEW FIX] Averaged specific heat fixed in TMDMN

parent 0e11b824
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
+ 64
− 54
View file @ b4610983
......@@ -3884,12 +3884,12 @@ void StochDMNDG3DMaterialLaw::fill_NLocal(const int pos, double N[][3])const{
 
void StochDMNDG3DMaterialLaw::fill_W_WI(const int pos, const int level, double W[][2]){
 
// code to fill a matrix related to weights and interfaces
// code to fill a matrix related to weights and interfaces -> For composite nodes?
 
double VI = _VI[pos];
double VI = _VI[pos]; // V0 = 1 - VI
double Wtmp0 = _Para_Wt[pos][0];
double Wtmp1 = _Para_Wt[pos][1];
double WA = Wtmp0*(1.0-VI)+ Wtmp1*VI;
double WA = Wtmp0*(1.0-VI)+ Wtmp1*VI; // Eq. 39 in Ling's draft
 
for(int i=0; i<2; i++){
for(int j=0; j<2; j++){
......@@ -3898,8 +3898,8 @@ void StochDMNDG3DMaterialLaw::fill_W_WI(const int pos, const int level, double W
}
W[0][0]= WA;
W[0][1]= 1.0-WA;
W[1][0]= (Wtmp1*VI)/WA;
W[1][1]= ((1.0-Wtmp1)*VI)/(1.0-WA);
W[1][0]= (Wtmp1*VI)/WA; // Eq. 40 in Ling's draft, v_1A
W[1][1]= ((1.0-Wtmp1)*VI)/(1.0-WA); // Eq. 40 in Ling's draft, v_1B
}
 
 
......@@ -4011,7 +4011,7 @@ void StochDMNDG3DMaterialLaw::fill_Matrices(){
}
}
 
// information of voulme fraction to stress on node
// information of volume fraction to stress on node
Row_start = col*(_N_Node-1-_NRoot);
for(int nc=0; nc<col*_NRoot; nc++){
_V.set(Row_start+nc,nc,1.0);
......@@ -4303,13 +4303,13 @@ void StochDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, c
}
}
 
_NTV.mult(Pvec, Res_node);
_NTV.mult(Pvec, Res_node); // Eq 60 in Ling's draft, Res_node = residual at the node
r = Res_node.norm()/_NInterface;
 
if(r<_tol){
if(last or stiff_loc){
_NTV.mult(_C, NTVC);
NTVC.mult(_Nv, Jacobian);
NTVC.mult(_Nv, Jacobian); // Eq 65 in Ling's draft
Jacobian.invertInPlace();
break;
}
......@@ -4326,7 +4326,7 @@ void StochDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, c
Jacobian.invertInPlace();
}
else{
// What is this? HOw does the following update the jacobian?
// What is this? How does the following update the jacobian?
_NTV.mult(Pvec, Res_node);
Res_node_t.scale(-1.0);
Res_node_t.axpy(Res_node, 1.0);
......@@ -4379,7 +4379,7 @@ void StochDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, c
 
if(stiff){
NTVC.mult(_I, dRdF);
Jacobian.mult(dRdF, dadF);
Jacobian.mult(dRdF, dadF); // Eq. 68 in Ling's draft
_Nv.mult(dadF, tmp);
tmp.add(_I);
 
......@@ -4388,7 +4388,7 @@ void StochDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, c
C_hom[m][n] = 0.0;
for(int i=0;i<_NRoot;i++){
for(int j=0;j< 9*_NRoot;j++){
C_hom[m][n] += (_VA[1]*_V(m,m+9*i) + _VA[2]*_V(9+m,m+9*i))*_C(m+9*i,j)*tmp(j,n);
C_hom[m][n] += (_VA[1]*_V(m,m+9*i) + _VA[2]*_V(9+m,m+9*i))*_C(m+9*i,j)*tmp(j,n); // Eq. 70 in Ling's draft
}
}
}
......@@ -4408,7 +4408,7 @@ void StochDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, c
ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
}
// end of stress
// Compute derivatives repect to DMN parameters at trainning satge
// Compute derivatives respect to DMN parameters at training satge
void StochDMNDG3DMaterialLaw::getdVAdP(std::vector< std::vector <SPoint2> >& dVAdPv) const{
int pos0, pos1;
static std::vector< std::vector <SPoint2> > dVIdPv;
......@@ -5070,13 +5070,14 @@ StochTMDMNDG3DMaterialLaw::StochTMDMNDG3DMaterialLaw(const int num, const double
const bool porous, const bool flag_isothermal, const bool flag_microTempFixed, const double tol):
StochDMNDG3DMaterialLaw(num,rho,E,nu,ParaFile,porous,tol),_flag_isothermal(flag_isothermal),_flag_microTempFixed(flag_microTempFixed){
// Meanings of _col_Di and _row_Di are interchanged in this class.
if (_flag_isothermal && _flag_microTempFixed){
_col_Di = 9; _row_Di = 3;
}
else if(!_flag_isothermal && _flag_microTempFixed){
_col_Di = 12; _row_Di = 4;
}
else if(!_flag_isothermal && !_flag_microTempFixed){
else if(!_flag_isothermal && !_flag_microTempFixed){ // Deprecated option - Impossible to implement without length scale.
_col_Di = 13; _row_Di = 5;
}
else{
......@@ -5277,7 +5278,7 @@ void StochTMDMNDG3DMaterialLaw::fill_Matrices(){
double tmp;
 
// set the size of rows and cols in Di matrix
int col = _col_Di; // = 9;
int col = _col_Di; // = 9, 12;
 
for(int i=0; i <_N_Node; i++){
_VA.push_back(1.0); // Vector for position of nodes for Stress
......@@ -5289,7 +5290,7 @@ void StochTMDMNDG3DMaterialLaw::fill_Matrices(){
if(pos == int(pow(2,l+1)-1)) l +=1; // if its the first position (interface) of every level, l+=1
fill_NLocal(pos,Norm); // size of Norm depends on the flags set inside this function
pos_kids(pos, l, pos0, pos1);
Row_start = pos*_row_Di; // 3;
Row_start = pos*_row_Di; // 3, 4;
Col_start0 = (pos0-1)*col;
Col_start1 = (pos1-1)*col;
 
......@@ -5303,7 +5304,7 @@ void StochTMDMNDG3DMaterialLaw::fill_Matrices(){
}
}
 
if(l ==_TotalLevel-1){
if(l ==_TotalLevel-1){ // Last level
if(_porous){
_VA[pos0] = (1.0-_VI[pos])*_Para_Wt[pos][0];
_VA[pos1] = (1.0-_VI[pos])*(1.0-_Para_Wt[pos][0]);
......@@ -5313,13 +5314,13 @@ void StochTMDMNDG3DMaterialLaw::fill_Matrices(){
_VA[pos1] = _VI[pos];
}
}
else if(_Modified_DMN and l ==_TotalLevel-2){
else if(_Modified_DMN and l ==_TotalLevel-2){ // Composite nodes - 2nd last level
StochDMNDG3DMaterialLaw::fill_W_WI(pos,l,W_WI);
_VA[pos0] = W_WI[0][0];
_VA[pos1] = W_WI[0][1];
_VI.push_back(W_WI[1][0]);
_VI.push_back(W_WI[1][0]); // Collect v_1A for all composite nodes
}
else{
else{ // Composite nodes - all other levels
StochDMNDG3DMaterialLaw::fill_W_WI(pos,l,W_WI);
_VA[pos0] = W_WI[0][0];
_VA[pos1] = W_WI[0][1];
......@@ -5327,7 +5328,7 @@ void StochTMDMNDG3DMaterialLaw::fill_Matrices(){
_VI.push_back(W_WI[1][1]);
}
//information of homogenous strain to local strain
Col_start = pos*_row_Di; // *3
Col_start = pos*_row_Di; // *3, 4 // Note that the meanings of _row_Di and _col_Di have been swapped
if(!_Modified_DMN or l< _TotalLevel-2){
k = int(pow(2,l));
node = _NRoot/k/2;
......@@ -5375,7 +5376,7 @@ void StochTMDMNDG3DMaterialLaw::fill_Matrices(){
}
}
 
// information of voulme fraction to stress on node
// information of volume fraction to stress on node
Row_start = col*(_N_Node-1-_NRoot);
for(int nc=0; nc<col*_NRoot; nc++){
_V.set(Row_start+nc,nc,1.0);
......@@ -5449,7 +5450,7 @@ void StochTMDMNDG3DMaterialLaw::read_parameter(const char *ParaFile){
_ParaAngle.resize(_NInterface,2,true);
double ParaN[2];
double ParaN[2]; // W_alpha and w_beta -> Eq. 37 in Ling's draft
double Pi(3.14159265359);
SPoint3 Para_Norm;
SPoint2 Para_Wt;
......@@ -5601,8 +5602,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 P_hom[_col_Di]; //static
// std::vector<double> P_hom;
double PQ_hom[_col_Di]; //static
double Cp_hom(0.), thermSrc_hom(0.), mechSource_hom(0.);
// std::vector<double> PQ_hom;
 
fullVector<double>& FHvec = ipvcur->getRefToCombinedGradientVect();
fullVector<double>& PQvec = ipvcur->getRefToCombinedStressVect();
......@@ -5624,7 +5626,7 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
Res_node_t.setAll(0.0);
Res_node.setAll(0.0);
a_step.setAll(0.0);
ab_step.setAll(0.0);
_C.setAll(0.0);
NTVC.setAll(0.0);
Jacobian.setAll(0.0);
......@@ -5750,7 +5752,6 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
const double& dmechSourcedt_i = ipv_i->getRefTodMechanicalSourcedField()(0,0);
const STensor3& dmechSourcedf_i = ipv_i->getRefTodMechanicalSourcedF()[0];
// CHECK with python AND DEBUG the tangents
if (_flag_isothermal && _flag_microTempFixed){
for(int m=0; m<3; m++)
for(int n=0; n<3; n++)
......@@ -5762,7 +5763,7 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
for(int m=0; m<3; m++)
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+m) = dqdgradT_i(m,n);
_C(i*_col_Di + 9+m, i*_col_Di + 9+n) = dqdgradT_i(m,n);
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
......@@ -5813,7 +5814,6 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
stiff_loc = true;
}
}
else{
if(stiff_loc){
_NTV.mult(_C, NTVC);
......@@ -5821,17 +5821,17 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
Jacobian.invertInPlace();
}
else{
// What is this? HOw does the following update the jacobian?
// What is this? How does the following update the jacobian?
_NTV.mult(PQvec, Res_node);
Res_node_t.scale(-1.0);
Res_node_t.axpy(Res_node, 1.0);
 
double Rho = 1.0/(Res_node_t*a_step);
double Rho = 1.0/(Res_node_t*ab_step);
Modify_Jacobian.setAll(0.0);
for(int i =0; i<3*_NInterface; i++){
Modify_Jacobian(i,i) = 1.0;
for(int j =0; j<3*_NInterface; j++){
Modify_Jacobian(i,j) -= a_step(i)*Res_node_t(j)*Rho;
Modify_Jacobian(i,j) -= ab_step(i)*Res_node_t(j)*Rho;
}
}
Modify_Jacobian.mult(Jacobian,Jacobian_inv);
......@@ -5840,7 +5840,7 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
 
for(int i =0; i<3*_NInterface; i++){
for(int j =0; j<3*_NInterface; j++){
Jacobian(i,j) += a_step(i)*a_step(j)*Rho;
Jacobian(i,j) += ab_step(i)*ab_step(j)*Rho;
}
}
stiff_loc = true;
......@@ -5848,56 +5848,66 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
}
 
Res_node_t = Res_node;
Jacobian.mult(Res_node, a_step);
a_step.scale(-1.0);
ab_cur.axpy(a_step, 1.0);
Jacobian.mult(Res_node, ab_step);
ab_step.scale(-1.0);
ab_cur.axpy(ab_step, 1.0);
}
}
Jacobian.scale(-1.0);
 
if(r>_tol){
Msg::Error("DMN residual (= %lf) did n't converge to tolenerce after %d iteration",r,ite);
Msg::Error("TMDMN residual (= %lf) didn't converge to tolenerce after %d iteration",r,ite);
return;
}
else{
for(int m=0;m<_col_Di;m++){
P_hom[m] = 0.0;
PQ_hom[m] = 0.0;
for(int i=0;i<_NRoot;i++){
P_hom[m] += (_VA[1]*_V(m, _col_Di*i+m) + _VA[2]*_V(_col_Di+m, _col_Di*i+m)) *PQvec(_col_Di*i+m); // CHECK
PQ_hom[m] += (_VA[1]*_V(m, _col_Di*i+m) + _VA[2]*_V(_col_Di+m, _col_Di*i+m)) *PQvec(_col_Di*i+m); // Stress average
}
}
Cp_hom, thermSrc_hom, mechSource_hom = 0.;
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
thermSrc_hom += (_VA[1]*_V(0, _col_Di*i) + _VA[2]*_V(_col_Di, _col_Di*i)) *thermSrcvec(i); // w average
mechSource_hom += (_VA[1]*_V(0, _col_Di*i) + _VA[2]*_V(_col_Di, _col_Di*i)) *mechSrcvec(i); // mechSrc average
}
if (_flag_isothermal && _flag_microTempFixed){
// Cp =
Cp = Cp_hom; w = 0.; mechSource = 0.;
for(int i=0;i<3;i++){
Q(i) = 0.;
for(int j=0;j<3;j++){
P(i,j) = P_hom[i+_row_Di*j]; // CHECK
P(i,j) = PQ_hom[i+_row_Di*j]; // CHECK
}
}
}
else if(!_flag_isothermal && _flag_microTempFixed){
Cp = Cp_hom; w = thermSrc_hom; mechSource = mechSource_hom;
for(int i=0;i<3;i++){
// Q(i) = P_hom[(i+1)*_row_Di-1];
Q(i) = P_hom[9+i];
// Q(i) = PQ_hom[(i+1)*_row_Di-1];
Q(i) = PQ_hom[9+i];
for(int j=0;j<3;j++){
// P(i,j) = P_hom[i+_row_Di*j]; // CHECK
P(i,j) = P_hom[i+3*j]; // CHECK
// P(i,j) = PQ_hom[i+_row_Di*j];
P(i,j) = PQ_hom[i+3*j];
}
}
}
else if(!_flag_isothermal && !_flag_microTempFixed){
Cp = P_hom[_col_Di-1];
Cp = PQ_hom[_col_Di-1]; w = thermSrc_hom; mechSource = mechSource_hom;
for(int i=0;i<3;i++){
// Q(i) = P_hom[(i+1)*(_row_Di-1)-1];
Q(i) = P_hom[9+i];
// Q(i) = PQ_hom[(i+1)*(_row_Di-1)-1];
Q(i) = PQ_hom[9+i];
for(int j=0;j<3;j++){
// P(i,j) = P_hom[i+(_row_Di-1)*j]; // CHECK
P(i,j) = P_hom[i+3*j]; // CHECK
// P(i,j) = PQ_hom[i+(_row_Di-1)*j];
P(i,j) = PQ_hom[i+3*j];
}
}
}
 
if(stiff){
// All derivatives must be complete
NTVC.mult(_I, dRdF);
Jacobian.mult(dRdF, dadF);
_Nv.mult(dadF, tmp);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment