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

[NEW FIX 3] Compatibility with StochDMN

parent ae75fc05
Branches
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
+ 72
− 154
View file @ 9e364366
......@@ -4999,9 +4999,9 @@ void StochDMNDG3DMaterialLaw::setTime(const double t,const double dtime){
}
 
// FLE
StochTMDMNDG3DMaterialLaw::StochTMDMNDG3DMaterialLaw(const int num, const double rho, const double E,const double nu, const char *ParaFile,
StochTMDMNDG3DMaterialLaw::StochTMDMNDG3DMaterialLaw(const int num, const double rho, const double E,const double nu, const char *ParaFile, const bool ReadTree,
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){
StochDMNDG3DMaterialLaw(num,rho,E,nu,ParaFile,ReadTree,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){
......@@ -5016,8 +5016,12 @@ StochTMDMNDG3DMaterialLaw::StochTMDMNDG3DMaterialLaw(const int num, const double
else{
Msg::Error("This combination of flags doesn't exist. Redefine in the constructor.");
}
if(_ReadTree){
StochDMNDG3DMaterialLaw::read_TreeData(ParaFile);
}
else{
StochTMDMNDG3DMaterialLaw::read_parameter(ParaFile);
StochTMDMNDG3DMaterialLaw::fill_Matrices();
}
 
......@@ -5031,29 +5035,22 @@ void StochTMDMNDG3DMaterialLaw::createIPState(IPStateBase* &ips, bool hasBodyFor
bool inter=true;
const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
if(iele==NULL) inter=false;
IPVariable* ipvi = new StochTMDMNDG3DIPVariable(_NRoot, _NInterface, _col_Di, _row_Di, hasBodyForce,inter);
IPVariable* ipv1 = new StochTMDMNDG3DIPVariable(_NRoot, _NInterface, _col_Di, _row_Di, hasBodyForce,inter);
IPVariable* ipv2 = new StochTMDMNDG3DIPVariable(_NRoot, _NInterface, _col_Di, _row_Di, hasBodyForce,inter);
IPVariable* ipvi = new StochTMDMNDG3DIPVariable(_NLeaf, _NInterface, _col_Di, _row_Di, hasBodyForce,inter);
IPVariable* ipv1 = new StochTMDMNDG3DIPVariable(_NLeaf, _NInterface, _col_Di, _row_Di, hasBodyForce,inter);
IPVariable* ipv2 = new StochTMDMNDG3DIPVariable(_NLeaf, _NInterface, _col_Di, _row_Di, hasBodyForce,inter);
 
StochTMDMNDG3DIPVariable* ipvi_root = static_cast<StochTMDMNDG3DIPVariable*>(ipvi);
StochTMDMNDG3DIPVariable* ipv1_root = static_cast<StochTMDMNDG3DIPVariable*>(ipv1);
StochTMDMNDG3DIPVariable* ipv2_root = static_cast<StochTMDMNDG3DIPVariable*>(ipv2);
for (int i=0; i< _NRoot; i++){
StochTMDMNDG3DIPVariable* ipvi_Leaf = static_cast<StochTMDMNDG3DIPVariable*>(ipvi);
StochTMDMNDG3DIPVariable* ipv1_Leaf = static_cast<StochTMDMNDG3DIPVariable*>(ipv1);
StochTMDMNDG3DIPVariable* ipv2_Leaf = static_cast<StochTMDMNDG3DIPVariable*>(ipv2);
for (int i=0; i< _NLeaf; i++){
int mat_IP = 0;
IPStateBase* ips_i = NULL;
if(_Modified_DMN){
if(i%3 == 1) mat_IP = 1;
}
else{
mat_IP = i%2;
}
if(_porous) mat_IP=0;
_mapLaw[mat_IP]->createIPState(ips_i, hasBodyForce, state_, ele, nbFF_, GP, gpt);
std::vector<IPVariable*> ip_all;
ips_i->getAllIPVariable(ip_all);
ipvi_root->addIPv(i, ip_all[0]);
ipv1_root->addIPv(i, ip_all[1]);
ipv2_root->addIPv(i, ip_all[2]);
ipvi_Leaf->addIPv(i, ip_all[0]);
ipv1_Leaf->addIPv(i, ip_all[1]);
ipv2_Leaf->addIPv(i, ip_all[2]);
}
if(ips != NULL) delete ips;
ips = new IP3State(state_,ipvi,ipv1,ipv2);
......@@ -5066,20 +5063,14 @@ void StochTMDMNDG3DMaterialLaw::createIPVariable(IPVariable* &ipv, bool hasBodyF
if(iele == NULL){
inter=false;
}
ipv = new StochTMDMNDG3DIPVariable(_NRoot, _NInterface, _col_Di, _row_Di, hasBodyForce,inter);
StochTMDMNDG3DIPVariable* ipv_root = static_cast<StochTMDMNDG3DIPVariable*>(ipv);
for (int i=0; i< _NRoot; i++){
ipv = new StochTMDMNDG3DIPVariable(_NLeaf, _NInterface, _col_Di, _row_Di, hasBodyForce,inter);
StochTMDMNDG3DIPVariable* ipv_Leaf = static_cast<StochTMDMNDG3DIPVariable*>(ipv);
int mat_IP;
for (int i=0; i< _NLeaf; i++){
IPVariable* ipv_i = NULL;
int mat_IP =0;
if(_Modified_DMN){
if(i%3 == 1) mat_IP = 1;
}
else{
mat_IP = i%2;
}
if(_porous) mat_IP=0;
mat_IP = _Mat_Id[_NInterface + i];
_mapLaw[mat_IP]->createIPVariable(ipv_i, hasBodyForce, ele, nbFF_, GP, gpt);
ipv_root->addIPv(i, ipv_i);
ipv_Leaf->addIPv(i, ipv_i);
}
};
 
......@@ -5089,16 +5080,10 @@ void StochTMDMNDG3DMaterialLaw::initialIPVariable(IPVariable* ipv, bool stiff){
if (ipv_all == NULL){
Msg::Error("StochTMDMNDG3DIPVariable must be used in StochTMDMNDG3DMaterialLaw::initialIPVariable");
}
for (int i=0; i< _NRoot; i++){
int mat_IP;
for (int i=0; i< _NLeaf; i++){
IPVariable* ipv_i = ipv_all->getIPv(i);
int mat_IP =0;
if(_Modified_DMN){
if(i%3 == 1) mat_IP = 1;
}
else{
mat_IP = i%2;
}
if(_porous) mat_IP=0;
mat_IP = _Mat_Id[_NInterface + i];
_mapLaw[mat_IP]->initialIPVariable(ipv_i, stiff);
}
}
......@@ -5109,17 +5094,11 @@ void StochTMDMNDG3DMaterialLaw::checkInternalState(IPVariable* ipv, const IPVari
if (ipv_all == NULL){
Msg::Error("StochTMDMNDG3DIPVariable must be used in StochTMDMNDG3DMaterialLaw::checkInternalState");
}
for (int i=0; i< _NRoot; i++){
int mat_IP;
for (int i=0; i< _NLeaf; i++){
IPVariable* ipv_i = ipv_all->getIPv(i);
const IPVariable* ipvprev_i = ipvprev_all->getIPv(i);
int mat_IP = 0;
if(_Modified_DMN){
if(i%3 == 1) mat_IP = 1;
}
else{
mat_IP = i%2;
}
if(_porous) mat_IP=0;
mat_IP = _Mat_Id[_NInterface + i];
_mapLaw[mat_IP]->checkInternalState(ipv_i, ipvprev_i);
}
}
......@@ -5151,9 +5130,6 @@ void StochTMDMNDG3DMaterialLaw::fill_NLocal(const int pos, std::vector<std::vect
N.resize(_col_Di);
for(auto& row : N)
row.resize(_row_Di,0.);
// N[0][0] = n0; N[1][1] = n0; N[2][2] = n0; N[3][3] = n0;
// N[4][0] = n1; N[5][1] = n1; N[6][2] = n1; N[7][3] = n1;
// N[8][0] = n2; N[9][1] = n2; N[10][2] = n2; N[11][3] = n2;
N[0][0] = n0; N[1][1] = n0; N[2][2] = n0;
N[3][0] = n1; N[4][1] = n1; N[5][2] = n1;
N[6][0] = n2; N[7][1] = n2; N[8][2] = n2;
......@@ -5163,11 +5139,6 @@ void StochTMDMNDG3DMaterialLaw::fill_NLocal(const int pos, std::vector<std::vect
N.resize(_col_Di);
for(auto& row : N)
row.resize(_row_Di,0.);
// N[0][0] = n0; N[1][1] = n0; N[2][2] = n0; N[3][3] = n0;
// N[4][0] = n1; N[5][1] = n1; N[6][2] = n1; N[7][3] = n1;
// N[8][0] = n2; N[9][1] = n2; N[10][2] = n2; N[11][3] = n2;
// N[12][4] = 1.;
N[0][0] = n0; N[1][1] = n0; N[2][2] = n0;
N[3][0] = n1; N[4][1] = n1; N[5][2] = n1;
N[6][0] = n2; N[7][1] = n2; N[8][2] = n2;
......@@ -5179,24 +5150,8 @@ void StochTMDMNDG3DMaterialLaw::fill_NLocal(const int pos, std::vector<std::vect
}
}
 
void StochTMDMNDG3DMaterialLaw::pos_kids(const int pos, const int level, int& pos0, int& pos1) const {
// code to calculate positions of child nodes using given "pos" and "level"
if(!_Modified_DMN or level<_TotalLevel-2 ){
pos0 = 2*pos+1;
pos1 = 2*pos+2;
}
else if(level ==_TotalLevel-2){
pos0 = pos + pow(2, level);
int k = pos - (pow(2, level)-1);
pos1 = int(3*pow(2, level)-1) + k*3+2; // CHECK
}
else if(level==_TotalLevel-1){
int k = pos - (pow(2, level)-1);
pos0 = int(3*pow(2, level-1)-1) +k*3; // CHECK
pos1 = pos0+1;
}
void StochTMDMNDG3DMaterialLaw::initLaws(const std::map<int,materialLaw*> &maplaw){
// To define this function, StochDMNDG3DMaterialLaw::TwoPhasesInteract has to be redefined in this class.
}
 
// fill the matrices which keep the topological information
......@@ -5204,25 +5159,28 @@ void StochTMDMNDG3DMaterialLaw::fill_Matrices(){
 
std::vector<std::vector<double>> Norm; // double Norm[9][3];
double W_WI[2][2];
int pos0=0;
int pos1=0;
int k, node;
int pos0, pos1;
int L_ch0, L_ch1;
int node_start0, node_end0, node_start1, node_end1;
int Col_start, Col_start0, Col_start1, Row_start, Row_start0, Row_start1;
double tmp;
 
// set the size of rows and cols in Di matrix
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
}
// fill matrix for equilibrium of interface
_VI.push_back(_Vf);
int l = -1;
_VA.resize(_N_Node,0.0); // Vector for position of nodes for Stress
_VA[0] = 1.0;
_VI.resize(_NInterface, 0.0); // fill matrix for equilibrium of interface
_VI[0] = _Vf;
int l;
for(int pos=0; pos <_NInterface; pos++){
if(pos == int(pow(2,l+1)-1)) l +=1; // if its the first position (interface) of every level, l+=1
l = _NodeLevel[pos];
fill_NLocal(pos,Norm); // size of Norm depends on the flags set inside this function
pos_kids(pos, l, pos0, pos1);
pos0 = _Child[pos][0];
pos1 = _Child[pos][1];
L_ch0 = _NodeLevel[pos0];
L_ch1 = _NodeLevel[pos1];
Row_start = pos*_row_Di; // 3, 4;
Col_start0 = (pos0-1)*col;
Col_start1 = (pos1-1)*col;
......@@ -5247,62 +5205,30 @@ void StochTMDMNDG3DMaterialLaw::fill_Matrices(){
_VA[pos1] = _VI[pos];
}
}
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]); // Collect v_1A for all composite nodes
}
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];
_VI.push_back(W_WI[1][0]);
_VI.push_back(W_WI[1][1]);
if(L_ch0 < _TotalLevel) _VI[pos0]=W_WI[1][0];
if(L_ch1 < _TotalLevel) _VI[pos1]=W_WI[1][1];
// Msg::Info("Node= %d,%d,%d, VI = %f,%f,%f",pos,pos0,pos1,_VI[pos],W_WI[1][0],W_WI[1][1]);
}
//information of homogenous strain to local strain
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;
int i = pos-int(pow(2,l))+1;
for(int r=0; r<node; r++){
Row_start0 = (i*2*node+r)*col;
Row_start1 = (i*2*node+r+node)*col;
for(int nr=0; nr<col;nr++){
for(int nc=0; nc<_row_Di;nc++){
_Nv.set(Row_start0+nr, Col_start+nc,Norm[nr][nc]/_VA[pos0]);
_Nv.set(Row_start1+nr, Col_start+nc,-S_ba*Norm[nr][nc]/_VA[pos1]);
}
}
}
}
else if(l == _TotalLevel-2){
node = _row_Di; // 3;
int i = pos-int(pow(2,l))+1;
for(int r=0; r<2; r++){
Row_start0 = (i*node+r)*col;
StochDMNDG3DMaterialLaw::getLeafOfInterface(pos0, node_start0, node_end0);
StochDMNDG3DMaterialLaw::getLeafOfInterface(pos1, node_start1, node_end1);
for(int r=node_start0; r<=node_end0; r++){
Row_start0 = (r-_NInterface)*col;
for(int nr=0; nr<col;nr++){
for(int nc=0; nc<_row_Di;nc++){
_Nv.set(Row_start0+nr, Col_start+nc, Norm[nr][nc]/_VA[pos0]);
}
}
}
Row_start1 = (i*node+2)*col;
for(int nr=0; nr<col;nr++){
for(int nc=0; nc<_row_Di; nc++){
_Nv.set(Row_start1+nr, Col_start+nc,-S_ba*Norm[nr][nc]/_VA[pos1]);
}
}
}
else if(l == _TotalLevel-1){
node = _row_Di; // 3;
int i = pos-int(pow(2,l))+1;
Row_start0 = (i*node)*col;
Row_start1 = (i*node+1)*col;
for(int r=node_start1; r<=node_end1; r++){
Row_start1 = (r-_NInterface)*col;
for(int nr=0; nr<col;nr++){
for(int nc=0; nc<_row_Di;nc++){
_Nv.set(Row_start0+nr, Col_start+nc,Norm[nr][nc]/_VA[pos0]);
_Nv.set(Row_start1+nr, Col_start+nc,-S_ba*Norm[nr][nc]/_VA[pos1]);
}
}
......@@ -5310,20 +5236,19 @@ void StochTMDMNDG3DMaterialLaw::fill_Matrices(){
}
 
// information of volume fraction to stress on node
Row_start = col*(_N_Node-1-_NRoot);
for(int nc=0; nc<col*_NRoot; nc++){
Row_start = col*(_N_Node-1-_NLeaf);
for(int nc=0; nc<col*_NLeaf; nc++){
_V.set(Row_start+nc,nc,1.0);
}
 
l = _TotalLevel-1;
for(int pos = _NInterface-1; pos > 0; pos--){
if(pos == int(pow(2,l)-2)) l -=1;
pos_kids(pos, l, pos0, pos1);
pos0 = _Child[pos][0];
pos1 = _Child[pos][1];
Row_start = (pos-1)*col;
Row_start0 = (pos0-1)*col;
Row_start1 = (pos1-1)*col;
for(int nr=0; nr<col;nr++){
for(int nc=0; nc<col*_NRoot ; nc++){
for(int nc=0; nc<col*_NLeaf ; nc++){
tmp = _VA[pos0]*_V(Row_start0+nr,nc) + _VA[pos1]*_V(Row_start1+nr,nc);
_V.set(Row_start+nr,nc, tmp);
}
......@@ -5331,7 +5256,7 @@ void StochTMDMNDG3DMaterialLaw::fill_Matrices(){
}
_NT.mult(_V, _NTV);
 
for(int i=0; i<col*_NRoot;i++){
for(int i=0; i<col*_NLeaf;i++){
_I(i, i%col) = 1.0;
}
}
......@@ -5341,24 +5266,17 @@ void StochTMDMNDG3DMaterialLaw::fill_Matrices(){
void StochTMDMNDG3DMaterialLaw::read_parameter(const char *ParaFile){
 
double sin0, cos0, sin1, cos1;
int Model;
FILE *Para = fopen(ParaFile, "r");
if ( Para != NULL ){
int okf = fscanf(Para, "%d\n", &Model);
_Modified_DMN = static_cast<bool>(Model);
okf = fscanf(Para, "%d %d\n", &_Dim, &_TotalLevel);
int okf = fscanf(Para, "%d %d\n", &_Dim, &_TotalLevel);
 
bool resizeFlag;
if(_Modified_DMN){
_NRoot = int(pow(2,(_TotalLevel-2))*3); // CHECK THIS FLAG
}
else{
_NRoot = int(pow(2,_TotalLevel));
}
_NInterface = _NRoot-1;
_N_Node = _NRoot + _NInterface;
_NLeaf = int(pow(2,_TotalLevel));
_NInterface = _NLeaf-1;
_N_Node = _NLeaf + _NInterface;
 
// START HERE
int row = _col_Di*_NRoot;
int col = _row_Di*_NInterface;
resizeFlag = _Nv.resize(row, col, true);
......
dG3D/src/dG3DMaterialLaw.h
+ 6
− 6
View file @ 9e364366
......@@ -709,8 +709,8 @@ class StochDMNDG3DMaterialLaw : public dG3DMaterialLaw{
std::vector<std::vector<int>> _Child; // Child[Node][0] = left child ; Child[Node][1] = right child of Node
std::vector<SPoint2> _Para_Wt;
std::vector<SPoint3> _Para_Norm;
std::vector<SPoint2> _Para_Wt; // vector of weights
std::vector<SPoint3> _Para_Norm; // vector of normals
fullMatrix<double> _ParaAngle;
std::vector<dG3DMaterialLaw*> _mapLaw;
fullMatrix<double> _Nv;
......@@ -785,11 +785,10 @@ class StochTMDMNDG3DMaterialLaw : public StochDMNDG3DMaterialLaw{
void fill_Matrices();
void read_parameter(const char *ParaFile);
void fill_NLocal(const int pos, std::vector<std::vector<double>> N) const;
void pos_kids(const int pos, const int level, int& pos0, int& pos1) const;
#endif //SWIG
public:
StochTMDMNDG3DMaterialLaw(const int num, const double rho, const double E,const double nu, const char *ParaFile,
const bool porous, const bool flag_isothermal, const bool flag_microTempFixed, const double tol=1e-6);
StochTMDMNDG3DMaterialLaw(const int num, const double rho, const double E,const double nu, const char *ParaFile, const bool ReadTree=false,
const bool porous = false, const bool flag_isothermal = false, const bool flag_microTempFixed = true, const double tol=1e-6);
#ifndef SWIG
StochTMDMNDG3DMaterialLaw(const StochTMDMNDG3DMaterialLaw &source);
virtual ~StochTMDMNDG3DMaterialLaw();
......@@ -800,6 +799,7 @@ class StochTMDMNDG3DMaterialLaw : public StochDMNDG3DMaterialLaw{
virtual void createIPState(IPStateBase* &ips, bool hasBodyForce, const bool* state_=NULL,const MElement *ele=NULL, const int nbFF_=0, const IntPt *GP=NULL, const int gpt = 0) const;
virtual void createIPVariable(IPVariable* &ipv, bool hasBodyForce, const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const;
virtual void initialIPVariable(IPVariable* ipv, bool stiff);
virtual void initLaws(const std::map<int,materialLaw*> &maplaw);
virtual double scaleFactor() const {return 1.;};
virtual double soundSpeed() const {return 0.;};
virtual materialLaw* clone() const {return new StochTMDMNDG3DMaterialLaw(*this);};
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment