From a1491eec6a25a0e74ffe2e91955fdd60e16781ca Mon Sep 17 00:00:00 2001
From: Mohib <smmustafa@uliege.be>
Date: Thu, 1 May 2025 01:35:15 +0200
Subject: [PATCH] [Minor - torchSCULatticeLaw] - HouseKeeping

---
 dG3D/src/dG3DMaterialLaw.cpp | 125 ++++++++++-------------------------
 dG3D/src/dG3DMaterialLaw.h   |  69 +++++++++----------
 2 files changed, 68 insertions(+), 126 deletions(-)

diff --git a/dG3D/src/dG3DMaterialLaw.cpp b/dG3D/src/dG3DMaterialLaw.cpp
index 48855bd72..9ca07d50d 100644
--- a/dG3D/src/dG3DMaterialLaw.cpp
+++ b/dG3D/src/dG3DMaterialLaw.cpp
@@ -3072,11 +3072,11 @@ void torchANNBasedDG3DMaterialLaw::RNNstress_stiff(const fullMatrix<double>& E0,
 
 #if defined(HAVE_TORCH)
 // subroutine for computing stress and stiffness using extra inputs (other than Strain and internal variables) @Mohib TODO: Rename
-//appended arguments with E0 and Extra0 to feature enable SCU @Mohib
-void torchANNBasedDG3DMaterialLaw::RNNstressGeo_stiff(const fullMatrix<double>& E0, const fullMatrix<double>& E1, const fullMatrix<double>& Extra0, const fullMatrix<double>& Extra1, const torch::Tensor& h0, torch::Tensor& h1,
-                                      fullMatrix<double>& S, const bool stiff, fullMatrix<double>& DSDE)
-{
-
+// appended arguments with E0 and Extra0 to feature enable SCU @Mohib
+void torchANNBasedDG3DMaterialLaw::RNNstressGeo_stiff(const fullMatrix<double>& E0, const fullMatrix<double>& E1, 
+                                                      const fullMatrix<double>& Extra0, const fullMatrix<double>& Extra1, 
+                                                      const torch::Tensor& h0, torch::Tensor& h1, 
+                                                      fullMatrix<double>& S, const bool stiff, fullMatrix<double>& DSDE){
    // Containers for Strain and extra inputs @ Mohib
    vector<float> E_vec(_numberOfInput - _numberOfExtraInput);
    // Appended E0 to feature enable SCU @Mohib
@@ -3096,7 +3096,6 @@ void torchANNBasedDG3DMaterialLaw::RNNstressGeo_stiff(const fullMatrix<double>&
    // Container for  normalized combined input exteries(Extra + Input) @ Mohib
    torch::Tensor Combine0_norm;
 
-   // @Mohib TODO: Remove b4 commit
    torch::Tensor E_norm;
    // Appended E0 and dE_norm to feature enable SCU
    torch::Tensor E0_norm, dE_norm;
@@ -3118,41 +3117,18 @@ void torchANNBasedDG3DMaterialLaw::RNNstressGeo_stiff(const fullMatrix<double>&
       Normalize_strain(E0, E0_vec);
       // Normalize extra inputs at previous time step and store in Extra0_vec. Necessary for feature enabling SCU @Mohib
       Normalize_geo(Extra0, Extra0_vec);
-
-// Changes in SCU with Radius correction dont require this combine vector anymore TODO: Cleanup b4 commit
-//      // Populate the container "Combine_vec0" with the normalized quantities in the correct order.
-//      // Order is : Radius + E0 for SCU. Last entry in Extra0_vec contains time and isnt needed here @Mohib.
-//      Combine0_vec.insert(Combine0_vec.end(), Extra0_vec.begin(), Extra0_vec.end() - 1);
-//      Combine0_vec.insert(Combine0_vec.end(), E0_vec.begin(), E0_vec.end());
-
    }
 
    if(_NeedExtraNorm and _DoubleInput){
-
-// Changes in SCU with Radius correction dont require this combine vector anymore TODO: Cleanup b4 commit
-//      // Populate the container "Combine_vec" with the normalized quantities in the correct order.
-//      // Order is : Radius + E for SCU . Last entry in Extra_vec contains time and isnt needed here @Mohib.
-//      Combine_vec.insert(Combine_vec.end(), Extra_vec.begin(), Extra_vec.end() - 1);
-//      Combine_vec.insert(Combine_vec.end(), E_vec.begin(), E_vec.end());
-
       if(stiff){
          E0_norm = torch::from_blob(E0_vec.data(), {1,1, _numberOfInput - _numberOfExtraInput}, torch::requires_grad());
-
-// Changes in SCU with Radius correction dont require this combine vector anymore TODO: Cleanup b4 commit
-//         Combine0_norm = torch::from_blob(Combine0_vec.data(), {1,1, _numberOfInput - 1}, torch::requires_grad());
       }
       else{
          E0_norm = torch::from_blob(E0_vec.data(), {1,1, _numberOfInput - _numberOfExtraInput}, torch::requires_grad(false));
-// Changes in SCU with Radius correction dont require this combine vector anymore TODO: Cleanup b4 commit
-//         Combine0_norm = torch::from_blob(Combine0_vec.data(), {1,1, _numberOfInput - 1}, torch::requires_grad(false));
       }
-
-// Changes in SCU with Radius correction dont require pushing this combine vector into the input jit anymore TODO: Cleanup b4 commit
-//      inputs.push_back(Combine0_norm);
       inputs.push_back(E0_norm);
    }
    else{
-
       // Populate the container "Combine_vec" with the normalized quantities in the correct order.
       // Order is : Radius, Size, TSTEP + E for GRU @Mohib.
       Combine_vec.insert(Combine_vec.end(), Extra_vec.begin(), Extra_vec.end());
@@ -3168,10 +3144,7 @@ void torchANNBasedDG3DMaterialLaw::RNNstressGeo_stiff(const fullMatrix<double>&
       }
       inputs.push_back(Combine_norm);
    }
-
-  //  // @Mohib TODO: fix 3 and 6 cases
-  //  E_norm = torch::from_blob(E_vec.data(), {1,1, _numberOfInput - _numOfExtra}, torch::requires_grad());
-   
+ 
    // added norm of E1 and delta E_vec to feature enable SCU @Mohib
    vector<float> norm_E1(1);
    // delta DE_vec doesnt include time and radius so decrementing by num of extea variables @Mohib.
@@ -3180,9 +3153,7 @@ void torchANNBasedDG3DMaterialLaw::RNNstressGeo_stiff(const fullMatrix<double>&
   
    if(_kinematicInput == torchANNBasedDG3DMaterialLaw::EGLInc or _NeedExtraNorm){     
        if(_kinematicInput == torchANNBasedDG3DMaterialLaw::EGL and _NeedExtraNorm){           
-          for (int i = 0; i <_numberOfInput - _numberOfExtraInput; i++) {
-            // Changes in SCU with Radius correction dont require this combine vector anymore TODO: Cleanup b4 commit
-            //DE_vec[i] = Combine_vec[i]- Combine0_vec[i];
+          for (int i = 0; i <_numberOfInput - _numberOfExtraInput; i++){
           	DE_vec[i] = E_vec[i]-E0_vec[i];
           }
           if(_DoubleInput){
@@ -3194,32 +3165,16 @@ void torchANNBasedDG3DMaterialLaw::RNNstressGeo_stiff(const fullMatrix<double>&
             }  
           inputs.push_back(dE_norm); 
           }   
-//          norm_E1[0] = sqrt(std::inner_product(DE_vec.begin(), DE_vec.end(), DE_vec.begin(), 0.0));
+          // norm_E1[0] = sqrt(std::inner_product(DE_vec.begin(), DE_vec.end(), DE_vec.begin(), 0.0));
           norm_E1[0] = 0;
           for (int i = 0; i <_numberOfInput - _numberOfExtraInput; i++ ) {
-			norm_E1[0] += DE_vec[i] * DE_vec[i];
+			      norm_E1[0] += DE_vec[i] * DE_vec[i];
           }
           norm_E1[0] = sqrt(norm_E1[0]);
-       } 
-       else{
-
-// Changes in SCU with Radius correction dont require this combine vector anymore TODO: Cleanup b4 commit
-//          fullMatrix<double> CombineE1(1, _numberOfInput - 1);
-//
-//
-//          for (int i = 0; i < _numberOfExtraInput - 1; i++)
-//          {
-//            CombineE1(0, i) = Extra1(0, i);
-//          }
-//
-//          for (int i = 0; i < _numberOfInput - _numberOfExtraInput; i++)
-//          {
-//            CombineE1(0, _numberOfExtraInput - 1 + i) = E1(0, i);
-//          }
-//			norm_E1[0] = CombineE1.norm();
-          norm_E1[0] = E1.norm();
-
-       }       
+      } 
+      else{
+        norm_E1[0] = E1.norm();
+      }   
             
       if(stiff){
         norm = torch::from_blob(norm_E1.data(), {1,1,1}, torch::requires_grad());
@@ -3235,21 +3190,18 @@ void torchANNBasedDG3DMaterialLaw::RNNstressGeo_stiff(const fullMatrix<double>&
    torch::Tensor dt_norm; 
 
    if(_kinematicInput == torchANNBasedDG3DMaterialLaw::EGLInc or _NeedExtraNorm){     
-       if(_kinematicInput == torchANNBasedDG3DMaterialLaw::EGL and _NeedExtraNorm){   
-              // For SCU dt needs to be seperately normalized. @Mohib
-              Normalize_dt(Extra1, Extra_vec);
-              Normalize_dt(Extra0, Extra0_vec);
-
-              dt_vec[0] = Extra_vec[_numberOfExtraInput - 1] - Extra0_vec[_numberOfExtraInput - 1];
-
-
-              if(stiff){
-               dt_norm = torch::from_blob(dt_vec.data(), {1,1,1}, torch::requires_grad());
-              }
-              else{
-               dt_norm = torch::from_blob(dt_vec.data(), {1,1,1}, torch::requires_grad(false));
-            }  
-          inputs.push_back(dt_norm); 
+       if(_kinematicInput == torchANNBasedDG3DMaterialLaw::EGL and _NeedExtraNorm){
+        // For SCU dt needs to be seperately normalized. @Mohib
+        Normalize_dt(Extra1, Extra_vec);
+        Normalize_dt(Extra0, Extra0_vec);
+        dt_vec[0] = Extra_vec[_numberOfExtraInput - 1] - Extra0_vec[_numberOfExtraInput - 1];
+        if(stiff){
+          dt_norm = torch::from_blob(dt_vec.data(), {1,1,1}, torch::requires_grad());
+        }
+        else{
+          dt_norm = torch::from_blob(dt_vec.data(), {1,1,1}, torch::requires_grad(false));
+        }  
+        inputs.push_back(dt_norm); 
        }
    }
 
@@ -3259,21 +3211,16 @@ void torchANNBasedDG3DMaterialLaw::RNNstressGeo_stiff(const fullMatrix<double>&
 
    if(_kinematicInput == torchANNBasedDG3DMaterialLaw::EGLInc or _NeedExtraNorm){
        if(_kinematicInput == torchANNBasedDG3DMaterialLaw::EGL and _NeedExtraNorm){
-//              // For SCU dt needs to be seperately normalized. @Mohib
-//              Normalize_dt(Extra1, Extra_vec);
-//              Normalize_dt(Extra0, Extra0_vec);
-
-              r_vec[0] = Extra_vec[0];
-
-
-              if(stiff){
-               r_norm = torch::from_blob(r_vec.data(), {1,1,1}, torch::requires_grad());
-              }
-              else{
-               r_norm = torch::from_blob(r_vec.data(), {1,1,1}, torch::requires_grad(false));
-            }
-          inputs.push_back(r_norm);
-       }
+        // Extra_vec already contains normalized value
+        r_vec[0] = Extra_vec[0];
+        if(stiff){
+          r_norm = torch::from_blob(r_vec.data(), {1,1,1}, torch::requires_grad());
+        }
+        else{
+          r_norm = torch::from_blob(r_vec.data(), {1,1,1}, torch::requires_grad(false));
+        }
+        inputs.push_back(r_norm);
+      }
    }
 
    inputs.push_back(h0);
@@ -3564,7 +3511,7 @@ void torchANNBasedDG3DMaterialLaw::Normalize_geo(const fullMatrix<double>& Extra
 {
     if(_NeedExtraNorm and _DoubleInput){
       int j = 0;
-      // For SCU we only need to normalize radius dt will be done seperatly. @Mohib
+      // For SCU we only need to normalize radius, dt will be done seperatly. @Mohib
       for(int i=0; i<_numberOfExtraInput - 1; i++){
         if(_normParams_ExtraInput[j + 1] != 0.0){
             Geo_norm[i] = (Extra1(0, i) - _normParams_ExtraInput[j]) / _normParams_ExtraInput[j + 1];
diff --git a/dG3D/src/dG3DMaterialLaw.h b/dG3D/src/dG3DMaterialLaw.h
index a059ea83f..f4aabb343 100644
--- a/dG3D/src/dG3DMaterialLaw.h
+++ b/dG3D/src/dG3DMaterialLaw.h
@@ -602,45 +602,40 @@ class torchANNBasedDG3DMaterialLaw : public dG3DMaterialLaw{
                 const double SZZstd, const double SZXmean, const double SZXstd, bool pert=false, double tol = 1e-5);
 
 		void setKinematicInput(const int i);
-		void setNeedExtraNorm (const bool NeedExtraNorm){
-		    _NeedExtraNorm = NeedExtraNorm;
-		    if(_NeedExtraNorm) Msg::Info("Extra Norm is applied");
-                };
-
-                void setDoubleInput (const bool DoubleInput){
-		    _DoubleInput = DoubleInput;
-		    if(_DoubleInput) Msg::Info("Double Input is applied");
-                };
-
-
-        void setInitialHValue(double initialValue_h) {_initialValue_h=initialValue_h;};
-        // Setter for no of Extra Input Parameters @Mohib
-        void setNumExtraInput(const int value);
-        // Setter for storing time position (arg) in the extra array @Mohib
-        void setTimeArg(const int value);
-        // Setter for init values of Extra Inputs @Mohib
-        void setInitialExtraInput(const double value);
-        // Setter for storing norm params of Extra Inputs @Mohib
-        void setNormExtraInp(const double mean, const double std);
-
-
-
-        void setFileExtraInputsGp(const std::string fname)
-        {
-            _fileExtraInputsGp=fname;
-            _extraVariablesAtGaussPoints.fill(fname);
-            //_extraVariablesAtGaussPoints.print();
-            //const std::vector<double> &rf=_extraVariablesAtGaussPoints.getFieldAtClosestLocation(SVector3(25,28.873,28.873));
-            //double val=_extraVariablesAtGaussPoints.returnFieldAtClosestLocation(0,SVector3(25,28.873,28.873));
-            //Msg::Info(" Field 0 at GP is %f",val);
-        }
+		
+    void setNeedExtraNorm (const bool NeedExtraNorm){
+      _NeedExtraNorm = NeedExtraNorm;
+		  if(_NeedExtraNorm) Msg::Info("Extra Norm is applied");
+    };
+    void setDoubleInput (const bool DoubleInput){
+      _DoubleInput = DoubleInput;
+		  if(_DoubleInput) Msg::Info("Double Input is applied");
+    };
 
-        void fillExtraInputAtGp(int nb1, int nb2)
-        {
-           _extraVariablesAtGaussPoints.fillExtraInputAtGp(nb1,nb2);
-        }
+    void setInitialHValue(double initialValue_h) {_initialValue_h=initialValue_h;};
+    // Setter for no of Extra Input Parameters @Mohib
+    void setNumExtraInput(const int value);
+    // Setter for storing time position (arg) in the extra array @Mohib
+    void setTimeArg(const int value);
+    // Setter for init values of Extra Inputs @Mohib
+    void setInitialExtraInput(const double value);
+    // Setter for storing norm params of Extra Inputs @Mohib
+    void setNormExtraInp(const double mean, const double std);
+
+    void setFileExtraInputsGp(const std::string fname){
+      _fileExtraInputsGp=fname;
+      _extraVariablesAtGaussPoints.fill(fname);
+      //_extraVariablesAtGaussPoints.print();
+      //const std::vector<double> &rf=_extraVariablesAtGaussPoints.getFieldAtClosestLocation(SVector3(25,28.873,28.873));
+      //double val=_extraVariablesAtGaussPoints.returnFieldAtClosestLocation(0,SVector3(25,28.873,28.873));
+      //Msg::Info(" Field 0 at GP is %f",val);
+    }
+    void fillExtraInputAtGp(int nb1, int nb2){
+      _extraVariablesAtGaussPoints.fillExtraInputAtGp(nb1,nb2);
+    }
     #ifndef SWIG
-		torchANNBasedDG3DMaterialLaw(const torchANNBasedDG3DMaterialLaw& src);
+		
+    torchANNBasedDG3DMaterialLaw(const torchANNBasedDG3DMaterialLaw& src);
 		virtual ~torchANNBasedDG3DMaterialLaw();
 		virtual materialLaw::matname getType() const {return materialLaw::torchANN;}
     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;
-- 
GitLab