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

[NEW FIX 4] Adapted TMDMN to StochDMN class

parent 9e364366
No related branches found
No related tags found
1 merge request!430Added 1 modified class for tempGrad derivatives in Network Interactions,...
...@@ -5017,16 +5017,17 @@ StochTMDMNDG3DMaterialLaw::StochTMDMNDG3DMaterialLaw(const int num, const double ...@@ -5017,16 +5017,17 @@ StochTMDMNDG3DMaterialLaw::StochTMDMNDG3DMaterialLaw(const int num, const double
Msg::Error("This combination of flags doesn't exist. Redefine in the constructor."); Msg::Error("This combination of flags doesn't exist. Redefine in the constructor.");
} }
if(_ReadTree){ if(_ReadTree){
StochDMNDG3DMaterialLaw::read_TreeData(ParaFile); StochTMDMNDG3DMaterialLaw::read_TreeData(ParaFile);
} }
else{ else{
StochTMDMNDG3DMaterialLaw::read_parameter(ParaFile); StochTMDMNDG3DMaterialLaw::read_parameter(ParaFile);
}
   
StochTMDMNDG3DMaterialLaw::fill_Matrices(); StochTMDMNDG3DMaterialLaw::fill_Matrices();
} }
   
StochTMDMNDG3DMaterialLaw::StochTMDMNDG3DMaterialLaw(const StochTMDMNDG3DMaterialLaw &src): StochDMNDG3DMaterialLaw(src), StochTMDMNDG3DMaterialLaw::StochTMDMNDG3DMaterialLaw(const StochTMDMNDG3DMaterialLaw &source): StochDMNDG3DMaterialLaw(source),
_col_Di(src._col_Di), _row_Di(src._row_Di), _flag_isothermal(src._flag_isothermal), _flag_microTempFixed(src._flag_microTempFixed), _Para_Alpha(src._Para_Alpha){} _col_Di(source._col_Di), _row_Di(source._row_Di), _flag_isothermal(source._flag_isothermal), _flag_microTempFixed(source._flag_microTempFixed), _Para_Alpha(source._Para_Alpha){};
   
StochTMDMNDG3DMaterialLaw::~StochTMDMNDG3DMaterialLaw(){}; StochTMDMNDG3DMaterialLaw::~StochTMDMNDG3DMaterialLaw(){};
   
...@@ -5206,7 +5207,7 @@ void StochTMDMNDG3DMaterialLaw::fill_Matrices(){ ...@@ -5206,7 +5207,7 @@ void StochTMDMNDG3DMaterialLaw::fill_Matrices(){
} }
} }
else{ // Composite nodes - all other levels else{ // Composite nodes - all other levels
StochDMNDG3DMaterialLaw::fill_W_WI(pos,l,W_WI); StochDMNDG3DMaterialLaw::fill_W_WI(pos,W_WI);
_VA[pos0] = W_WI[0][0]; _VA[pos0] = W_WI[0][0];
_VA[pos1] = W_WI[0][1]; _VA[pos1] = W_WI[0][1];
if(L_ch0 < _TotalLevel) _VI[pos0]=W_WI[1][0]; if(L_ch0 < _TotalLevel) _VI[pos0]=W_WI[1][0];
...@@ -5276,55 +5277,43 @@ void StochTMDMNDG3DMaterialLaw::read_parameter(const char *ParaFile){ ...@@ -5276,55 +5277,43 @@ void StochTMDMNDG3DMaterialLaw::read_parameter(const char *ParaFile){
_NInterface = _NLeaf-1; _NInterface = _NLeaf-1;
_N_Node = _NLeaf + _NInterface; _N_Node = _NLeaf + _NInterface;
   
// START HERE int row = _col_Di*_NLeaf;
int row = _col_Di*_NRoot;
int col = _row_Di*_NInterface; int col = _row_Di*_NInterface;
resizeFlag = _Nv.resize(row, col, true); resizeFlag = _Nv.resize(row, col, true);
   
resizeFlag = _I.resize(row, _col_Di, true); resizeFlag = _I.resize(row, _col_Di, true);
   
row = _col_Di*(_N_Node-1); row = _col_Di*(_N_Node-1);
col = _col_Di*_NRoot; col = _col_Di*_NLeaf;
resizeFlag = _V.resize(row, col, true); resizeFlag = _V.resize(row, col, true);
row = _row_Di*_NInterface; row = _row_Di*_NInterface;
resizeFlag = _NTV.resize(row, col, true); resizeFlag = _NTV.resize(row, col, true);
resizeFlag = _NT.resize(_row_Di*_NInterface, _col_Di*(_N_Node-1), true); resizeFlag = _NT.resize(row, _col_Di*(_N_Node-1), true);
okf = fscanf(Para, "%lf\n", &_Vf); okf = fscanf(Para, "%lf\n", &_Vf);
   
_NPara_Norm = _NInterface; _NPara_Norm = _NInterface;
resizeFlag = _C.resize(_col_Di*_NRoot,_col_Di*_NRoot,true); resizeFlag = _C.resize(col,col,true);
resizeFlag = Jacobian.resize(_row_Di*_NInterface, _row_Di*_NInterface,true); resizeFlag = Jacobian.resize(row,row,true);
_ParaAngle.resize(_NInterface,2,true); resizeFlag = _ParaAngle.resize(_NInterface,2,true);
double ParaN[2]; // W_alpha and w_beta -> Eq. 37 in Ling's draft double ParaN[2]; // W_alpha and w_beta -> Eq. 37 in Ling's draft
double Pi(3.14159265359); double Pi(3.14159265359);
SPoint3 Para_Norm; SPoint3 Para_Norm;
SPoint2 Para_Wt; SPoint2 Para_Wt;
   
if(_Dim==2){
ParaN[1] = 0.0;
sin1 = sin(Pi*ParaN[1]);
cos1 = cos(Pi*ParaN[1]);
for(int i=0;i<_NInterface;i++){ for(int i=0;i<_NInterface;i++){
if(_Dim==2){
okf = fscanf(Para, "%lf\n", &ParaN[0]); okf = fscanf(Para, "%lf\n", &ParaN[0]);
sin0 = sin(2.0*Pi*ParaN[0]); ParaN[1] = 0.0;
cos0 = cos(2.0*Pi*ParaN[0]);
Para_Norm[0] = cos1*sin0;
Para_Norm[1] = cos1*cos0;
Para_Norm[2] = sin1;
_Para_Norm.push_back(Para_Norm);
_ParaAngle(i,0) = ParaN[0];
}
} }
else if(_Dim==3){ else if(_Dim==3){
for(int i=0;i<_NInterface;i++){
okf = fscanf(Para, "%lf %lf\n", &ParaN[0], &ParaN[1]); okf = fscanf(Para, "%lf %lf\n", &ParaN[0], &ParaN[1]);
}
sin0 = sin(2.0*Pi*ParaN[0]); sin0 = sin(2.0*Pi*ParaN[0]);
cos0 = cos(2.0*Pi*ParaN[0]); cos0 = cos(2.0*Pi*ParaN[0]);
sin1 = sin(Pi*ParaN[1]); sin1 = sin(Pi*ParaN[1]);
...@@ -5336,7 +5325,6 @@ void StochTMDMNDG3DMaterialLaw::read_parameter(const char *ParaFile){ ...@@ -5336,7 +5325,6 @@ void StochTMDMNDG3DMaterialLaw::read_parameter(const char *ParaFile){
_ParaAngle(i,0) = ParaN[0]; _ParaAngle(i,0) = ParaN[0];
_ParaAngle(i,1) = ParaN[1]; _ParaAngle(i,1) = ParaN[1];
} }
}
if(_porous){ if(_porous){
_NPara_Wt = _NInterface; _NPara_Wt = _NInterface;
...@@ -5349,49 +5337,178 @@ void StochTMDMNDG3DMaterialLaw::read_parameter(const char *ParaFile){ ...@@ -5349,49 +5337,178 @@ void StochTMDMNDG3DMaterialLaw::read_parameter(const char *ParaFile){
_Para_Wt.push_back(Para_Wt); _Para_Wt.push_back(Para_Wt);
} }
} }
Msg::Info("Finish reading parameter file, Model = %d, level = %d, Dim = %d, NRoot = %d",Model, _TotalLevel, _Dim, _NRoot); // creat data as for tree data
int n = 0;
int m = -10;
for(int l=0; l<=_TotalLevel; l++){
std::vector<int> ln;
_LevelNode.push_back(ln);
while(n < pow(2,l+1)-1){
_LevelNode[l].push_back(n);
_NodeLevel.push_back(l);
if(n==0){
_Parent.push_back(0);
}
else{
_Parent.push_back(int((n-1)/2));
}
if(l==_TotalLevel){
_Mat_Id.push_back((n+1)%2);
}
else{
std::vector<int> ch;
ch.push_back(2*n+1);
ch.push_back(2*n+2);
_Child.push_back(ch);
_Mat_Id.push_back(m);
}
n +=1;
}
}
Msg::Info("Finish reading parameter filefor Complete Binary Tree, level = %d, Dim = %d, NLeaf = %d", _TotalLevel, _Dim, _NLeaf);
} }
// end of read_parameter // end of read_parameter
   
void StochTMDMNDG3DMaterialLaw::read_TreeData(const char *TreeFile){
std::ifstream file(TreeFile);
if(!file.is_open()) {
Msg::Error( "Error: Could not open the tree file." );
return ;
}
int level_value, node_value, parent_value;
double sin0, cos0, sin1,cos1;
double ParaN[2];
int pn = 0;
double Pi(3.14159265359);
SPoint3 Para_Norm;
SPoint2 Para_Wt;
bool resizeFlag;
std::string line;
int DataId = 0;
while(std::getline(file, line)) {
if(line.empty()) continue;
std::istringstream stream(line);
std::string value;
std::vector<std::string> row;
while (std::getline(stream, value, ',')) {
row.push_back(value);
}
if (std::isalpha(row[0][0])){
DataId +=1;
if(DataId == 4){
_NLeaf = static_cast<int>(_LevelNode[_TotalLevel].size());
_N_Node = static_cast<int>(_Parent.size());
_NInterface = _N_Node - _NLeaf;
resizeFlag = _ParaAngle.resize(_NInterface,2,true);
}
}
else{
if(DataId == 1){ // read the tree structure------------------
_TotalLevel = std::stoi(row[0]);
for(int i=0; i<=_TotalLevel; i++){
std::vector<int> LN;
_LevelNode.push_back(LN);
}
}
else if(DataId == 2){
level_value = std::stoi(row[0]);
node_value = std::stoi(row[1]);
parent_value = std::stoi(row[2]);
_Mat_Id.push_back(std::stoi(row[3]));
_LevelNode[level_value].push_back(node_value);
_NodeLevel.push_back(level_value);
_Parent.push_back(parent_value);
if(level_value < _TotalLevel){
std::vector<int> ch;
_Child.push_back(ch);
}
if(parent_value != node_value) _Child[parent_value].push_back(node_value);
}
else if(DataId == 3){ // read the DMN Parameter -----------------------
_Dim = std::stoi(row[0]);
_Vf = std::stod(row[1]);
}
else if(DataId == 4){
if(_Dim==2){
ParaN[0] = std::stod(row[0]);
ParaN[1] = 0.0;
}
else if(_Dim==3){
ParaN[0] = std::stod(row[0]);
ParaN[1] = std::stod(row[1]);
}
sin0 = sin(2.0*Pi*ParaN[0]);
cos0 = cos(2.0*Pi*ParaN[0]);
sin1 = sin(Pi*ParaN[1]);
cos1 = cos(Pi*ParaN[1]);
Para_Norm[0] = cos1*sin0;
Para_Norm[1] = cos1*cos0;
Para_Norm[2] = sin1;
_Para_Norm.push_back(Para_Norm);
_ParaAngle(pn,0) = ParaN[0];
_ParaAngle(pn,1) = ParaN[1];
pn +=1;
}
else if(DataId == 5){
Para_Wt[0] = std::stod(row[0]);
Para_Wt[1] = std::stod(row[1]);
_Para_Wt.push_back(Para_Wt);
}
}
}
int row = _col_Di*_NLeaf;
int col = _row_Di*_NInterface;
resizeFlag = _Nv.resize(row, col, true);
resizeFlag = _I.resize(row, _col_Di, true);
row = _col_Di*(_N_Node-1);
col = _col_Di*_NLeaf;
resizeFlag = _V.resize(row, col, true);
row = _row_Di*_NInterface;
resizeFlag = _NT.resize(row, _col_Di*(_N_Node-1), true);
resizeFlag = _NTV.resize(row, col, true);
resizeFlag = _C.resize(col, col,true);
resizeFlag = Jacobian.resize(row, row,true);
_NPara_Norm = _NInterface;
_NPara_Wt = static_cast<int>(_Para_Wt.size());
Msg::Info("Finish reading Tree Data and parameter file, level = %d, Dim = %d, NLeaf = %d", _TotalLevel, _Dim, _NLeaf);
}
void StochTMDMNDG3DMaterialLaw::reset_Parameter(std::vector<std::vector<double>>& Para_Norm, std::vector<std::vector<double>>& Para_Wt, std::vector<double>& Para_Alpha, const double Vf){ void StochTMDMNDG3DMaterialLaw::reset_Parameter(std::vector<std::vector<double>>& Para_Norm, std::vector<std::vector<double>>& Para_Wt, std::vector<double>& Para_Alpha, const double Vf){
double sin0, cos0, sin1,cos1; double sin0, cos0, sin1,cos1;
double Pi(3.14159265359); double Pi(3.14159265359);
for(int i=0;i<_NInterface;i++){
if(_Dim==2){ if(_Dim==2){
sin1 = 0.0; sin1 = 0.0;
cos1 = 1.0; cos1 = 1.0;
for(int i=0;i<_NInterface;i++){
sin0 = sin(2.0*Pi*Para_Norm[i][0]);
cos0 = cos(2.0*Pi*Para_Norm[i][0]);
_Para_Norm[i][0] = cos1*sin0;
_Para_Norm[i][1] = cos1*cos0;
_Para_Norm[i][2] = sin1;
_ParaAngle(i,0) = Para_Norm[i][0]; _ParaAngle(i,0) = Para_Norm[i][0];
} }
}
else if(_Dim==3){ else if(_Dim==3){
for(int i=0;i<_NInterface;i++){
sin0 = sin(2.0*Pi*Para_Norm[i][0]);
cos0 = cos(2.0*Pi*Para_Norm[i][0]);
sin1 = sin(Pi*Para_Norm[i][1]); sin1 = sin(Pi*Para_Norm[i][1]);
cos1 = cos(Pi*Para_Norm[i][1]); cos1 = cos(Pi*Para_Norm[i][1]);
_Para_Norm[i][0] = cos1*sin0;
_Para_Norm[i][1] = cos1*cos0;
_Para_Norm[i][2] = sin1;
_ParaAngle(i,0) = Para_Norm[i][0]; _ParaAngle(i,0) = Para_Norm[i][0];
_ParaAngle(i,1) = Para_Norm[i][1]; _ParaAngle(i,1) = Para_Norm[i][1];
} }
sin0 = sin(2.0*Pi*Para_Norm[i][0]);
cos0 = cos(2.0*Pi*Para_Norm[i][0]);
_Para_Norm[i][0] = cos1*sin0;
_Para_Norm[i][1] = cos1*cos0;
_Para_Norm[i][2] = sin1;
} }
if(_porous){
_NPara_Wt = _NInterface;
}
else{
_NPara_Wt = int(pow(2,(_TotalLevel-1)))-1;
}
for(int i=0;i<_NPara_Wt;i++){ for(int i=0;i<_NPara_Wt;i++){
_Para_Wt[i][0] = Para_Wt[i][0]; _Para_Wt[i][0] = Para_Wt[i][0];
_Para_Wt[i][1] = Para_Wt[i][1]; _Para_Wt[i][1] = Para_Wt[i][1];
if(_NodeLevel[_Child[i][1]]==_TotalLevel) _Para_Wt[i][1] = 1.0;
} }
for(int i=0;i<Para_Alpha.size();i++){ for(int i=0;i<Para_Alpha.size();i++){
_Para_Alpha[i] = _Para_Alpha[i]; _Para_Alpha[i] = _Para_Alpha[i];
...@@ -5419,7 +5536,17 @@ void StochTMDMNDG3DMaterialLaw::reset_Parameter(const char* Para){ ...@@ -5419,7 +5536,17 @@ void StochTMDMNDG3DMaterialLaw::reset_Parameter(const char* Para){
_Para_Norm.clear(); _Para_Norm.clear();
_Para_Wt.clear(); _Para_Wt.clear();
_Para_Alpha.clear(); _Para_Alpha.clear();
_LevelNode.clear();
_NodeLevel.clear();
_Mat_Id.clear();
_Parent.clear();
_Child.clear();
if(_ReadTree){
StochTMDMNDG3DMaterialLaw::read_TreeData(Para);
}
else{
StochTMDMNDG3DMaterialLaw::read_parameter(Para); StochTMDMNDG3DMaterialLaw::read_parameter(Para);
}
StochTMDMNDG3DMaterialLaw::fill_Matrices(); StochTMDMNDG3DMaterialLaw::fill_Matrices();
} }
   
...@@ -5466,21 +5593,20 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, ...@@ -5466,21 +5593,20 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
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; fullVector<double> dwdtvec, dmechSrcdTvec, dwdFvec, dmechSrcdFvec;
dwdtvec.resize(_NRoot); dmechSrcdTvec.resize(_NRoot); dwdtvec.resize(_NLeaf); dmechSrcdTvec.resize(_NLeaf);
dwdFvec.resize(9*_NRoot); dmechSrcdFvec.resize(9*_NRoot); dwdFvec.resize(9*_NLeaf); dmechSrcdFvec.resize(9*_NLeaf);
   
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(_row_Di*_NInterface); static fullVector<double> Res_node(_row_Di*_NInterface);
static fullVector<double> ab_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*_NLeaf);
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);
static fullMatrix<double> tmp(_col_Di*_NRoot, _col_Di); static fullMatrix<double> tmp(_col_Di*_NLeaf, _col_Di);
static fullMatrix<double> dRdF(_row_Di*_NInterface, _col_Di); static fullMatrix<double> dRdF(_row_Di*_NInterface, _col_Di);
static fullMatrix<double> dadF(_row_Di*_NInterface, _col_Di); static fullMatrix<double> dadF(_row_Di*_NInterface, _col_Di);
Res_node_t.setAll(0.0);
Res_node.setAll(0.0); Res_node.setAll(0.0);
ab_step.setAll(0.0); ab_step.setAll(0.0);
_C.setAll(0.0); _C.setAll(0.0);
...@@ -5491,17 +5617,15 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, ...@@ -5491,17 +5617,15 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
tmp.setAll(0.0); tmp.setAll(0.0);
int ite = 0; int ite = 0;
double r = 1.0; double r = 1.0;
bool last = false; int pos_start = _NInterface;
int pos_start = _N_Node-_NRoot;
int max_iter= 15; int max_iter= 15;
bool stiff_loc = true; bool stiff_loc = true;
   
   
while( (r>_tol and ite < max_iter) or last){ while(r>_tol and ite < max_iter){
// if((1+ite)%5 == 0) stiff_loc=false;
ite +=1; ite +=1;
_Nv.mult(ab_cur, FHvec); _Nv.mult(ab_cur, FHvec);
for(int i=0; i<_NRoot; i++){ for(int i=0; i<_NLeaf; i++){
if (_VA[pos_start+i] > 1.e-6){ if (_VA[pos_start+i] > 1.e-6){
ThermoMechanicsDG3DIPVariableBase* ipv_i = dynamic_cast<ThermoMechanicsDG3DIPVariableBase*>(ipvcur->getIPv(i)); ThermoMechanicsDG3DIPVariableBase* ipv_i = dynamic_cast<ThermoMechanicsDG3DIPVariableBase*>(ipvcur->getIPv(i));
const ThermoMechanicsDG3DIPVariableBase* ipvprev_i = dynamic_cast<const ThermoMechanicsDG3DIPVariableBase*>(ipv_prev->getIPv(i)); const ThermoMechanicsDG3DIPVariableBase* ipvprev_i = dynamic_cast<const ThermoMechanicsDG3DIPVariableBase*>(ipv_prev->getIPv(i));
...@@ -5541,15 +5665,7 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, ...@@ -5541,15 +5665,7 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
Msg::Error("This combination of flags doesn't exist. Redefine in the constructor."); Msg::Error("This combination of flags doesn't exist. Redefine in the constructor.");
} }
   
int mat_IP = 0; int mat_IP = _Mat_Id[_NInterface + i];
if(_Modified_DMN){
if(i%3 == 1) mat_IP = 1; // CHECK THIS FLAG
}
else{
mat_IP = i%2;
}
if(_porous) mat_IP=0;
_mapLaw[mat_IP]->stress(ipv_i, ipvprev_i, stiff_loc, checkfrac, dTangent); _mapLaw[mat_IP]->stress(ipv_i, ipvprev_i, stiff_loc, checkfrac, dTangent);
const STensor3& P_i = ipv_i->getConstRefToFirstPiolaKirchhoffStress(); const STensor3& P_i = ipv_i->getConstRefToFirstPiolaKirchhoffStress();
const SVector3& Q_i = ipv_i->getConstRefToThermalFlux(); const SVector3& Q_i = ipv_i->getConstRefToThermalFlux();
...@@ -5653,57 +5769,18 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, ...@@ -5653,57 +5769,18 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
_NTV.mult(PQvec, Res_node); _NTV.mult(PQvec, Res_node);
r = Res_node.norm()/_NInterface; r = Res_node.norm()/_NInterface;
if(r<_tol){
if(last or stiff_loc){
_NTV.mult(_C, NTVC); _NTV.mult(_C, NTVC);
NTVC.mult(_Nv, Jacobian); NTVC.mult(_Nv, Jacobian);
Jacobian.invertInPlace(); Jacobian.invertInPlace();
break;
}
else{
last = true;
stiff_loc = true;
}
}
else{
if(stiff_loc){
_NTV.mult(_C, NTVC);
NTVC.mult(_Nv, Jacobian);
Jacobian.invertInPlace();
}
else{
// Updating 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*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) -= ab_step(i)*Res_node_t(j)*Rho;
}
}
Modify_Jacobian.mult(Jacobian,Jacobian_inv);
Modify_Jacobian.transposeInPlace();
Jacobian_inv.mult(Modify_Jacobian,Jacobian);
for(int i =0; i<3*_NInterface; i++){
for(int j =0; j<3*_NInterface; j++){
Jacobian(i,j) += ab_step(i)*ab_step(j)*Rho;
}
}
stiff_loc = true;
}
Res_node_t = Res_node;
Jacobian.mult(Res_node, ab_step); Jacobian.mult(Res_node, ab_step);
ab_step.scale(-1.0); ab_step.scale(-1.0);
ab_cur.axpy(ab_step, 1.0); ab_cur.axpy(ab_step, 1.0);
} }
}
Jacobian.scale(-1.0); Jacobian.scale(-1.0);
int pos0 = _Child[0][0];
int pos1 = _Child[0][1];
   
if(r>_tol){ if(r>_tol){
Msg::Error("TMDMN residual (= %lf) didn't converge to tolenerce after %d iteration",r,ite); Msg::Error("TMDMN residual (= %lf) didn't converge to tolenerce after %d iteration",r,ite);
...@@ -5712,14 +5789,14 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, ...@@ -5712,14 +5789,14 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
else{ else{
for(int m=0;m<_col_Di;m++){ for(int m=0;m<_col_Di;m++){
PQ_hom[m] = 0.0; PQ_hom[m] = 0.0;
for(int i=0;i<_NRoot;i++){ for(int i=0;i<_NLeaf;i++){
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 PQ_hom[m] += (_VA[pos0]*_V((pos0-1)*_col_Di+m, _col_Di*i+m) + _VA[pos1]*_V((pos1-1)*_col_Di+m, _col_Di*i+m))*PQvec(_col_Di*i+m); // Stress average
} }
} }
double temp = 0.; double temp = 0.;
for(int i=0;i<_NRoot;i++){ for(int i=0;i<_NLeaf;i++){
temp = (_VA[1]*_V(0, _col_Di*i) + _VA[2]*_V(_col_Di, _col_Di*i)); temp = (_VA[pos0]*_V(0, _col_Di*i) + _VA[pos1]*_V(_col_Di, _col_Di*i));
Cp_hom += temp*Cpvec(i); // Cp average Cp_hom += temp*Cpvec(i); // Cp average
thermSrc_hom += temp*thermSrcvec(i); // w average thermSrc_hom += temp*thermSrcvec(i); // w average
mechSource_hom += temp*mechSrcvec(i); // mechSrc average mechSource_hom += temp*mechSrcvec(i); // mechSrc average
...@@ -5754,26 +5831,26 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, ...@@ -5754,26 +5831,26 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
} }
   
if(stiff){ if(stiff){
// All derivatives must be complete
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); // Check this structure of _I. tmp = dFidF tmp.add(_I);
   
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<_NLeaf;i++){
for(int j=0;j< _col_Di*_NRoot;j++){ for(int j=0;j< _col_Di*_NLeaf;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); C_hom[m][n] += (_VA[pos0]*_V((pos0-1)*_col_Di+m,m+_col_Di*i) + _VA[pos1]*_V((pos1-1)*_col_Di+m,m+_col_Di*i))*_C(m+_col_Di*i,j)*tmp(j,n);
} }
} }
} }
} }
   
double temp(0.), dTidT(1.); // = Cp_hom/(temp*Cpvec(i)) double temp(0.), dTidT(1.); // = Cp_hom/(temp*Cpvec(i))
for(int i=0;i<_NRoot;i++){ for(int i=0;i<_NLeaf;i++){
temp = (_VA[1]*_V(0, _col_Di*i) + _VA[2]*_V(_col_Di, _col_Di*i)); temp = (_VA[pos0]*_V(0, _col_Di*i) + _VA[pos1]*_V(_col_Di, _col_Di*i));
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++)
......
...@@ -784,6 +784,7 @@ class StochTMDMNDG3DMaterialLaw : public StochDMNDG3DMaterialLaw{ ...@@ -784,6 +784,7 @@ class StochTMDMNDG3DMaterialLaw : public StochDMNDG3DMaterialLaw{
std::vector<double> _Para_Alpha; // vector of coefficients std::vector<double> _Para_Alpha; // vector of coefficients
void fill_Matrices(); void fill_Matrices();
void read_parameter(const char *ParaFile); void read_parameter(const char *ParaFile);
void read_TreeData(const char *TreeFile);
void fill_NLocal(const int pos, std::vector<std::vector<double>> N) const; void fill_NLocal(const int pos, std::vector<std::vector<double>> N) const;
#endif //SWIG #endif //SWIG
public: public:
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment