From 673ba8804cc922e2becbc284e8c4a642a5384fb7 Mon Sep 17 00:00:00 2001
From: FLE_Knight <ujwalkishore.jinaga@uliege.be>
Date: Thu, 29 May 2025 15:25:39 +0200
Subject: [PATCH] [MINOR PATCH and NEW FEATURE] Fixed a bug in StochDMN with
regards to rotation matrix implementation. Added a new feature for user input
euler angles
---
NonLinearSolver/materialLaw/mlawGenericTM.cpp | 2 +-
dG3D/src/dG3DMaterialLaw.cpp | 15 ++++++++-------
dG3D/src/dG3DMaterialLaw.h | 8 +++++---
3 files changed, 14 insertions(+), 11 deletions(-)
diff --git a/NonLinearSolver/materialLaw/mlawGenericTM.cpp b/NonLinearSolver/materialLaw/mlawGenericTM.cpp
index 34284fd4d..37e849bdd 100644
--- a/NonLinearSolver/materialLaw/mlawGenericTM.cpp
+++ b/NonLinearSolver/materialLaw/mlawGenericTM.cpp
@@ -28,7 +28,7 @@ mlawGenericTM::mlawGenericTM(const int num,
_cp = new constantScalarFunction(cp);
- // Euler angles to Rotation Matrix -> ZXZ convention
+ // Euler angles to Extrinsic Rotation Matrix -> ZXZ convention
STensor3 R; //3x3 rotation matrix
double c1,c2,c3,s1,s2,s3;
double s1c2, c1c2;
diff --git a/dG3D/src/dG3DMaterialLaw.cpp b/dG3D/src/dG3DMaterialLaw.cpp
index 6b0ed3fbd..7a3edfaab 100644
--- a/dG3D/src/dG3DMaterialLaw.cpp
+++ b/dG3D/src/dG3DMaterialLaw.cpp
@@ -3634,7 +3634,7 @@ double torchANNBasedDG3DMaterialLaw::soundSpeed() const
StochDMNDG3DMaterialLaw::StochDMNDG3DMaterialLaw(const int num, const double rho, const double E,const double nu, const char *ParaFile, const bool ReadTree, const bool porous,const double tol, const bool flag_isothermal):
- dG3DMaterialLaw(num,rho,false), _ReadTree(ReadTree), _porous(porous),_tol(tol), _flag_isothermal(flag_isothermal), _flag_rotation(false), _type_of_R(0){
+ dG3DMaterialLaw(num,rho,false), _ReadTree(ReadTree), _porous(porous),_tol(tol), _flag_isothermal(flag_isothermal), _flag_rotation(false), _type_of_R(0), _alpha(0.), _beta(0.), _gamma(0.){
if(_flag_isothermal){
if(_ReadTree){
@@ -3654,7 +3654,7 @@ StochDMNDG3DMaterialLaw::StochDMNDG3DMaterialLaw(const StochDMNDG3DMaterialLaw &
_NInterface(src._NInterface), _N_Node(src._N_Node), _NPara_Wt(src._NPara_Wt), _NPara_Norm(src._NPara_Norm), _Dim(src._Dim), _Vf(src._Vf), _tol(src._tol),
_Nv(src._Nv), _NTV(src._NTV), _V(src._V), _Para_Wt(src._Para_Wt), _Para_Norm(src._Para_Norm), _VA(src._VA), _VI(src._VI), _mapLaw(src._mapLaw),
_flag_isothermal(src._flag_isothermal), _flag_rotation(src._flag_rotation), _NPara_EulerAngles(src._NPara_EulerAngles), _Para_EulerAngles(src._Para_EulerAngles),
- _type_of_R(src._type_of_R), _R_user_input(src._R_user_input){}
+ _type_of_R(src._type_of_R), _R_user_input(src._R_user_input), _alpha(src._alpha), _beta(src._beta), _gamma(src._gamma){}
StochDMNDG3DMaterialLaw::~StochDMNDG3DMaterialLaw(){};
@@ -4347,6 +4347,9 @@ void StochDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, c
else if (_type_of_R == 1){
R = _R_user_input;
}
+ else if (_type_of_R == 2){
+ StochDMNDG3DMaterialLaw::eulerZXZToRotationMatrix(_alpha,_beta,_gamma,R); // user input euler angles
+ }
static STensor3 RT; STensorOperation::zero(RT);
STensorOperation::transposeSTensor3(R,RT);
@@ -4462,11 +4465,6 @@ void StochDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, c
}
}
- if (_flag_rotation){
- const STensor3& R = ipvcur->getConstReftoRotationMatrix();
- STensorOperation::multSTensor3InPlace2nd(R,P);
- }
-
if(stiff){
NTVC.mult(_I, dRdF);
Jacobian.mult(dRdF, dadF); // Eq. 68 in Ling's draft
@@ -5782,6 +5780,9 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
else if (_type_of_R == 1){
R = _R_user_input;
}
+ else if (_type_of_R == 2){
+ StochDMNDG3DMaterialLaw::eulerZXZToRotationMatrix(_alpha,_beta,_gamma,R); // user input euler angles
+ }
static STensor3 RT; STensorOperation::zero(RT);
STensorOperation::transposeSTensor3(R,RT);
diff --git a/dG3D/src/dG3DMaterialLaw.h b/dG3D/src/dG3DMaterialLaw.h
index 5a1bb7636..1ca18fe93 100644
--- a/dG3D/src/dG3DMaterialLaw.h
+++ b/dG3D/src/dG3DMaterialLaw.h
@@ -711,7 +711,7 @@ class StochDMNDG3DMaterialLaw : public dG3DMaterialLaw{
std::vector<SPoint2> _Para_Wt; // vector of weights
std::vector<SPoint3> _Para_Norm; // vector of normals
- std::vector<double> _Para_EulerAngles; // vector of rotation parameters
+ std::vector<double> _Para_EulerAngles; // vector of rotation parameters - traineable, not user defined
fullMatrix<double> _ParaAngle;
std::vector<dG3DMaterialLaw*> _mapLaw;
@@ -725,7 +725,8 @@ class StochDMNDG3DMaterialLaw : public dG3DMaterialLaw{
fullMatrix<double> _C;
fullMatrix<double> Jacobian;
- STensor3 _R_user_input;
+ double _alpha, _beta, _gamma; // user defined euler angles
+ STensor3 _R_user_input; // user defined rotation matrix
int _type_of_R; // type of rotation matrix
bool _flag_isothermal, _flag_rotation;
@@ -783,7 +784,7 @@ class StochDMNDG3DMaterialLaw : public dG3DMaterialLaw{
virtual void reset_Parameter(const char* Para);
virtual void reset_Parameter(std::vector<double>& Para_EulerAngles);
virtual void set_flag_DMN_response_rotation(const bool flag_rotation, const int type_of_R){_flag_rotation = flag_rotation; _type_of_R = type_of_R;};
- // virtual void setEulerAngles_DMN_response_rotation(const double alpha, const double beta, const double gamma){_alpha = alpha; _beta = beta; _gamma = gamma;};
+ virtual void set_euler_angles_user_input(const double alpha, const double beta, const double gamma){_alpha = alpha; _beta = beta; _gamma = gamma;};
};
// FLE
@@ -836,6 +837,7 @@ class StochTMDMNDG3DMaterialLaw : public StochDMNDG3DMaterialLaw{
virtual void reset_Parameter(const char* Para);
virtual void set_flag_DMN_response_rotation(const bool flag_rotation, const int type_of_R){_flag_rotation = flag_rotation; _type_of_R = type_of_R;};
virtual void writeMatrixToCSV(const fullMatrix<double>& matrix, const double& nrows, const double& ncols, const std::string& filename);
+ virtual void set_euler_angles_user_input(const double alpha, const double beta, const double gamma){_alpha = alpha; _beta = beta; _gamma = gamma;};
};
// FLE
--
GitLab