From 73c230e6245b2b7a58fe4e91816d9d80a00bf4c9 Mon Sep 17 00:00:00 2001
From: jleclerc <julien.leclerc@ulg.ac.be>
Date: Fri, 18 Jan 2019 11:15:27 +0100
Subject: [PATCH 01/11] add some message in debug mode only

---
 NonLinearSolver/materialLaw/STensorOperations.cpp | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/NonLinearSolver/materialLaw/STensorOperations.cpp b/NonLinearSolver/materialLaw/STensorOperations.cpp
index 968fc48c1..3af562176 100644
--- a/NonLinearSolver/materialLaw/STensorOperations.cpp
+++ b/NonLinearSolver/materialLaw/STensorOperations.cpp
@@ -12,7 +12,9 @@ bool STensorOperation::polynomial(const STensor3& a, const int order, const full
   static const STensor43 I4(1.,1.);
 
   if (STensorOperation::isnan(a)){
+    #ifdef _DEBUG
     Msg::Error("STensorOperation::polynomial of a NAN tensor");
+    #endif // _DEBUG
     STensorOperation::zero(F);
     if (dF){
       (*dF) = I4;
@@ -127,7 +129,9 @@ void STensorOperation::logSTensor3(const STensor3& a, const int order, STensor3&
     
     bool ok =  STensorOperation::polynomial(ami,order,coeffs,loga,dloga,ddloga);
     if (!ok){
+      #ifdef _DEBUG
       Msg::Error("error in computing logSTensor3");
+      #endif // _DEBUG
     }
   }
   else if (order == -1)
@@ -235,7 +239,9 @@ void STensorOperation::expSTensor3(const STensor3& a, const int order, STensor3&
     
     bool ok = STensorOperation::polynomial(a,order,coeffs,expa,dexpa,ddexpa);
     if (!ok){
+      #ifdef _DEBUG
       Msg::Error("error in computing expSTensor3");
+      #endif // _DEBUG
     }
   }
   else if (order == -1)
-- 
GitLab


From 57b705b173eb35c4b12fe4e3a8ab26bebdd25ec6 Mon Sep 17 00:00:00 2001
From: jleclerc <julien.leclerc@ulg.ac.be>
Date: Fri, 18 Jan 2019 11:16:09 +0100
Subject: [PATCH 02/11] update checkFailed function and penalty term

---
 dG3D/src/dG3DCohesiveBandMaterialLaw.cpp | 48 +++++++++++++-----------
 1 file changed, 26 insertions(+), 22 deletions(-)

diff --git a/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp b/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp
index 9522593c1..db83cbf7b 100644
--- a/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp
+++ b/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp
@@ -1062,7 +1062,11 @@ void PorosityCohesiveBand3DLaw::checkInternalState(IPVariable* ipv, const IPVari
       // interfaceIP should be deleted after satisfying failure criterion
     if (ipvbulk->getIPNonLocalPorosity()->isFailed()){
       fipv->Delete(true);
-      Msg::Warning("ipv is deleted");
+      Msg::Info("ipv is deleted");
+    }
+    else if (ipvbulk->getIPNonLocalPorosity()->getYieldPorosity() < 1.e-5){
+      fipv->Delete(true);
+      Msg::Info("ipv is deleted due to compression");
     }
     else{
       fipv->Delete(false);
@@ -1092,17 +1096,17 @@ void PorosityCohesiveBand3DLaw::stress(IPVariable* ipv, const IPVariable* ipvpre
     const SVector3& ujump = fipv->getConstRefToJump();
     const SVector3& ujumpDG = fipv->getConstRefToDGJump();
     double realUJumpNormal = dot(ujump,curN) - dot(ujumpDG,curN);
-    if (realUJumpNormal < 0){
-      fipv->setTensionFlag(true);
-      if(!fipvprev->ifTension()) Msg::Info("Penalty term added");
+    if (realUJumpNormal < 0.){
+      fipv->setTensionFlag(false);
+      if(fipvprev->ifTension()) Msg::Info("Penalty term added");
       for (int i=0; i<3; i++){
-        //T(i) += _Kp*realUJumpNormal*curN(i);
-        T(i) -= _Kp*realUJumpNormal*realUJumpNormal*curN(i);
+        T(i) += _Kp*realUJumpNormal*curN(i);
+        //T(i) -= _Kp*realUJumpNormal*realUJumpNormal*curN(i);
       }
     }
     else{
-      fipv->setTensionFlag(false);
-      if(fipvprev->ifTension()) Msg::Info("Penalty term deleted");
+      fipv->setTensionFlag(true);
+      if(!fipvprev->ifTension()) Msg::Info("Penalty term deleted");
     }
 
 
@@ -1112,12 +1116,12 @@ void PorosityCohesiveBand3DLaw::stress(IPVariable* ipv, const IPVariable* ipvpre
       STensorOperation::zero(fipv->getRefToDInterfaceForceDjump());
       STensorOperation::zero(fipv->getRefToDInterfaceForceDjumpGradient());
 
-      if (realUJumpNormal < 0){
+      if (realUJumpNormal < 0.){
         STensor3& dTdjump = fipv->getRefToDInterfaceForceDjump();
         for (int i=0; i<3; i++){
           for (int j=0; j<3; j++){
-            //dTdjump(i,j) += _Kp*curN(i)*curN(j);
-            dTdjump(i,j) -= 2.*_Kp*realUJumpNormal*curN(i)*curN(j);
+            dTdjump(i,j) += _Kp*curN(i)*curN(j);
+            //dTdjump(i,j) -= 2.*_Kp*realUJumpNormal*curN(i)*curN(j);
           }
         }
       }
@@ -1162,17 +1166,17 @@ void PorosityCohesiveBand3DLaw::stress(IPVariable* ipv, const IPVariable* ipvpre
       const SVector3& ujumpDG = fipv->getConstRefToDGJump();
 
       double realUJumpNormal = dot(ujump,curN) - dot(ujumpDG,curN);
-      if (realUJumpNormal < 0){
-        fipv->setTensionFlag(true);
-        if(!fipvprev->ifTension()) Msg::Info("Add penalty term");
+      if (realUJumpNormal < 0.){
+        fipv->setTensionFlag(false);
+        if(fipvprev->ifTension()) Msg::Info("Add penalty term");
         for (int i=0; i<3; i++){
-          T(i) -= _Kp*realUJumpNormal*realUJumpNormal*curN(i);
-          //T(i) += _Kp*realUJumpNormal*curN(i);
+          T(i) += _Kp*realUJumpNormal*curN(i);
+          //T(i) -= _Kp*realUJumpNormal*realUJumpNormal*curN(i);
         }
       }
       else{
-        fipv->setTensionFlag(false);
-        if(fipvprev->ifTension()) Msg::Info("Delete penalty term");
+        fipv->setTensionFlag(true);
+        if(!fipvprev->ifTension()) Msg::Info("Delete penalty term");
       }
 
 
@@ -1214,12 +1218,12 @@ void PorosityCohesiveBand3DLaw::stress(IPVariable* ipv, const IPVariable* ipvpre
           }
         }
 
-        if (realUJumpNormal < 0){
+        if (realUJumpNormal < 0.){
           for (int i=0; i<3; i++){
             for (int j=0; j<3; j++){
-              dTdjump(i,j) -= 2.*_Kp*realUJumpNormal*curN(i)*curN(j);
-              //dTdjump(i,j) += _Kp*curN(i)*curN(j);
-            };
+              dTdjump(i,j) += _Kp*curN(i)*curN(j);
+              //dTdjump(i,j) -= 2.*_Kp*realUJumpNormal*curN(i)*curN(j);
+            }
           }
         };
 
-- 
GitLab


From 9a66de9af25f7d5156106e7c021f684623513b22 Mon Sep 17 00:00:00 2001
From: jleclerc <julien.leclerc@ulg.ac.be>
Date: Fri, 18 Jan 2019 11:27:07 +0100
Subject: [PATCH 03/11] add a warning message in transferInterfaceDataToBulk in
 CBMlaw

---
 dG3D/src/dG3DCohesiveBandMaterialLaw.cpp | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp b/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp
index db83cbf7b..5940c6996 100644
--- a/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp
+++ b/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp
@@ -1450,7 +1450,9 @@ void PorosityCohesiveBand3DLaw::transferInterfaceDataToBulk(IPVariable* ipv, con
         F(i,j) += _gradFactor*(Gradjump(i,j) - Gradjump0(i,j));
       }
     }
-
+    
+    double detF = STensorOperation::determinantSTensor3(F);
+    if (detF < 1.e-10){Msg::Error("PorosityCohesiveBand3DLaw:: Negative jacobian= %e", detF);};
 
     // Compute derivatives
     if (stiff)
-- 
GitLab


From 63f779d00d7400795b96e6dc1eedd7995652a51e Mon Sep 17 00:00:00 2001
From: jleclerc <julien.leclerc@ulg.ac.be>
Date: Wed, 23 Jan 2019 10:39:57 +0100
Subject: [PATCH 04/11] add something to delete problematic IPVs

---
 dG3D/src/dG3DCohesiveBandMaterialLaw.cpp | 12 ++++++++----
 1 file changed, 8 insertions(+), 4 deletions(-)

diff --git a/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp b/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp
index 5940c6996..a2ba69cfa 100644
--- a/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp
+++ b/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp
@@ -1064,14 +1064,18 @@ void PorosityCohesiveBand3DLaw::checkInternalState(IPVariable* ipv, const IPVari
       fipv->Delete(true);
       Msg::Info("ipv is deleted");
     }
-    else if (ipvbulk->getIPNonLocalPorosity()->getYieldPorosity() < 1.e-5){
+    else if (ipvbulk->getIPNonLocalPorosity()->getYieldPorosity() < 0.75*_mlawPorous->getInitialPorosity()){
+      fipv->Delete(true);
+      Msg::Error("ipv is deleted due to too much compression: fv= %e", ipvbulk->getIPNonLocalPorosity()->getYieldPorosity());
+    }
+    else if(ipvbulk->getJ() < 1.e-10){
+      Msg::Error("ipv is deleted due to negative jacobian: J= %e", ipvbulk->getJ());
       fipv->Delete(true);
-      Msg::Info("ipv is deleted due to compression");
     }
     else{
       fipv->Delete(false);
-    };
-  };
+    }
+  }
 };
 
 
-- 
GitLab


From 8b84181770771b71564de918d3081020d829407b Mon Sep 17 00:00:00 2001
From: jleclerc <julien.leclerc@ulg.ac.be>
Date: Wed, 23 Jan 2019 14:13:51 +0100
Subject: [PATCH 05/11] modify tol on cl

---
 NonLinearSolver/materialLaw/CLengthLaw.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/NonLinearSolver/materialLaw/CLengthLaw.cpp b/NonLinearSolver/materialLaw/CLengthLaw.cpp
index 00ad7a5d9..1fc3e56c6 100644
--- a/NonLinearSolver/materialLaw/CLengthLaw.cpp
+++ b/NonLinearSolver/materialLaw/CLengthLaw.cpp
@@ -72,7 +72,7 @@ CLengthLaw * ZeroCLengthLaw::clone() const
 IsotropicCLengthLaw::IsotropicCLengthLaw(const int num, const double _cl, const bool init) : CLengthLaw(num,init),
 											cl(_cl)
 {
- double tol=1.e-10;
+ double tol=1.e-12;
  if(cl < tol) cl=tol;
 }
 IsotropicCLengthLaw::IsotropicCLengthLaw(const IsotropicCLengthLaw &source) :
-- 
GitLab


From a3c7ca9899161d487c2949ee9b8baa4957bf9d42 Mon Sep 17 00:00:00 2001
From: jleclerc <julien.leclerc@ulg.ac.be>
Date: Tue, 29 Jan 2019 11:23:38 +0100
Subject: [PATCH 06/11] add a limited contribution of negative jump

---
 dG3D/src/dG3DCohesiveBandMaterialLaw.cpp | 87 +++++++++++++++++-------
 dG3D/src/dG3DCohesiveBandMaterialLaw.h   | 23 ++++---
 2 files changed, 79 insertions(+), 31 deletions(-)

diff --git a/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp b/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp
index a2ba69cfa..e4432258b 100644
--- a/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp
+++ b/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp
@@ -990,12 +990,12 @@ void NonlocalPorosityCohesiveBand3DLaw::checkInternalState(IPVariable* ipv, cons
 
 
 
-
-
+// -------------- PorosityCohesiveBand3DLaw
 PorosityCohesiveBand3DLaw::PorosityCohesiveBand3DLaw(const int num, const double h, const double Kp, const bool accountNormalBulk):
     Cohesive3DLaw(num,true),
     _h(h),_Kp(Kp),
-    _accountNormalBulk(accountNormalBulk), _coalescenceLaw(NULL), _mlawPorous(NULL)
+    _accountNormalBulk(accountNormalBulk), _enabledNegativeJumpN_limit(true),_maxPenetrationjumpN(-1.e-8),_enableDeletionOfBadConditionedIPVs(false),
+    _coalescenceLaw(NULL), _mlawPorous(NULL)
 {
 
 };
@@ -1003,7 +1003,9 @@ PorosityCohesiveBand3DLaw::PorosityCohesiveBand3DLaw(const int num, const double
 
 PorosityCohesiveBand3DLaw::PorosityCohesiveBand3DLaw(const PorosityCohesiveBand3DLaw& src):
     Cohesive3DLaw(src),
-    _h(src._h),_Kp(src._Kp), _accountNormalBulk(src._accountNormalBulk)
+    _h(src._h),_Kp(src._Kp), 
+    _accountNormalBulk(src._accountNormalBulk), _enabledNegativeJumpN_limit(src._enabledNegativeJumpN_limit),
+    _maxPenetrationjumpN(src._maxPenetrationjumpN),_enableDeletionOfBadConditionedIPVs(src._enableDeletionOfBadConditionedIPVs)
 {
   _coalescenceLaw = NULL;
   if (src._coalescenceLaw!=NULL){
@@ -1046,42 +1048,56 @@ void PorosityCohesiveBand3DLaw::createIPVariable(IPVariable* &ipv,const MElement
   ipv = new CohesiveBand3DIPVariable();
 }
 
-void PorosityCohesiveBand3DLaw::checkInternalState(IPVariable* ipv, const IPVariable* ipvp) const{
+void PorosityCohesiveBand3DLaw::checkInternalState(IPVariable* ipv, const IPVariable* ipvp) const
+{
+  // Determine if the ipv is deleted or not.
+  /// if deleted, nothing is computed, except strains and penalty terms 
+  /// otherwise, compute interface force and their derivatives from ipv bulk information
+  
+  // get needed ipv
   FractureCohesive3DIPVariable *fipv = static_cast < FractureCohesive3DIPVariable *> (ipv);
   const FractureCohesive3DIPVariable *fipvprev = static_cast < const FractureCohesive3DIPVariable *> (ipvp);
-  if (fipvprev->isDeleted())
-  {
+  
+  if (fipvprev->isDeleted()){
+    // If ipv is already deleted, it always stays deleted
     fipv->Delete(true);
   }
   else{
-    // If the element is not deleted, compute interface force and their derivatives for ipv bulk information
-      // get ipvbulk
+    // If the ipv is not deleted, check if a criterion is satisfied
+    
+    // get needed ipvbulk
     const	nonLocalPorosityDG3DIPVariable* ipvbulk = dynamic_cast<nonLocalPorosityDG3DIPVariable*>(fipv->getIPvBulk());
 
-    // Deletion criterion:
+    // 1st deletion criterion:
       // interfaceIP should be deleted after satisfying failure criterion
     if (ipvbulk->getIPNonLocalPorosity()->isFailed()){
       fipv->Delete(true);
       Msg::Info("ipv is deleted");
     }
-    else if (ipvbulk->getIPNonLocalPorosity()->getYieldPorosity() < 0.75*_mlawPorous->getInitialPorosity()){
-      fipv->Delete(true);
-      Msg::Error("ipv is deleted due to too much compression: fv= %e", ipvbulk->getIPNonLocalPorosity()->getYieldPorosity());
-    }
-    else if(ipvbulk->getJ() < 1.e-10){
-      Msg::Error("ipv is deleted due to negative jacobian: J= %e", ipvbulk->getJ());
-      fipv->Delete(true);
+    // 2nd deletion criterion:
+      // interfaceIP should be deleted if considred as bad conditionned
+    else if (_enableDeletionOfBadConditionedIPVs){
+      if (ipvbulk->getIPNonLocalPorosity()->getYieldPorosity() < 0.75*_mlawPorous->getInitialPorosity()){
+        fipv->Delete(true);
+        Msg::Error("ipv is deleted due to too much compression: fv= %e", ipvbulk->getIPNonLocalPorosity()->getYieldPorosity());
+      }
+      else if(ipvbulk->getJ() < 1.e-10){
+        Msg::Error("ipv is deleted due to negative jacobian: J= %e", ipvbulk->getJ());
+        fipv->Delete(true);
+      }
     }
     else{
+      // if none of criterion is satisfied : the ipv is not deleted
       fipv->Delete(false);
     }
+
   }
 };
 
 
 void PorosityCohesiveBand3DLaw::stress(IPVariable* ipv, const IPVariable* ipvprev, const bool stiff, const bool checkfrac)
 {
-  // Get ipvfrac
+  // get ipvfrac
   FractureCohesive3DIPVariable *fipv = static_cast < FractureCohesive3DIPVariable *> (ipv);
   const FractureCohesive3DIPVariable *fipvprev = static_cast < const FractureCohesive3DIPVariable *> (ipvprev);
 
@@ -1245,14 +1261,14 @@ void PorosityCohesiveBand3DLaw::stress(IPVariable* ipv, const IPVariable* ipvpre
 									dirrEnergDF(k,l) += _h*dirrEnergdFinf(i,j)*dFinterfacedF(i,j,k,l);
 								}
 							}
-
 						}
 					}
 				}
-
+      
 		  }
 		}
 	}
+  
 };
 
 
@@ -1391,19 +1407,21 @@ void PorosityCohesiveBand3DLaw::checkCohesiveInsertion(IPStateBase* ipsm, IPStat
   };
 };
 
+
+
 void PorosityCohesiveBand3DLaw::transferInterfaceDataToBulk(IPVariable* ipv, const IPVariable* ipvprev, const bool stiff) const
 {
   // Update kinematic data if the cohesive zone is active (== if broken), otherwise, nothing to do
   FractureCohesive3DIPVariable *fipv = static_cast < FractureCohesive3DIPVariable *> (ipv);
   const FractureCohesive3DIPVariable *fipvprev = static_cast < const FractureCohesive3DIPVariable *> (ipvprev);
 
-  if (fipv->isbroken())
-  {
+  if (fipv->isbroken()){
     CohesiveBand3DIPVariable* cohIpv = dynamic_cast<CohesiveBand3DIPVariable*>(fipv->getIPvFrac());
     const CohesiveBand3DIPVariable* cohIpvprev = dynamic_cast<const CohesiveBand3DIPVariable*>(fipvprev->getIPvFrac());
 
     // Update deformation gradient for bulk law
     STensor3& F = fipv->getRefToDeformationGradient();
+    
     // - start from previous
     F = fipvprev->getConstRefToDeformationGradient();
 
@@ -1444,6 +1462,22 @@ void PorosityCohesiveBand3DLaw::transferInterfaceDataToBulk(IPVariable* ipv, con
         F(i,j) +=  (ujump(i) - ujump0(i))/_h*refN(j);
       }
     }
+    
+    // - substract normal part if asked
+    static SVector3 curN;
+    curN = fipvprev->getConstRefToCurrentOutwardNormal();
+    curN.normalize();
+    double jumpN = dot(ujump,curN)- dot(ujump0,curN);
+    
+    if(_enabledNegativeJumpN_limit){
+      if(jumpN < _maxPenetrationjumpN){
+        for (int i=0; i<3; i++){
+          for (int j=0; j<3; j++){
+            F(i,j) -= (jumpN-_maxPenetrationjumpN)*curN(i)/_h*refN(j);
+          }
+        }
+      }
+    }
 
     // - add jump gradient contribution
     double _gradFactor = 0.5;
@@ -1474,6 +1508,13 @@ void PorosityCohesiveBand3DLaw::transferInterfaceDataToBulk(IPVariable* ipv, con
         for (int j=0; j<3; j++){
           for (int k=0; k<3; k++){
             dFinterfacedJump(i,j,k) += I(i,k)*refN(j)/_h;
+            
+            if(_enabledNegativeJumpN_limit){
+              if(jumpN < _maxPenetrationjumpN){
+                dFinterfacedJump(i,j,k) -= curN(k)*curN(i)*refN(j)/_h;
+              }
+            }
+
             for (int l=0; l<3; l++){
               dFinterfacedF(i,j,k,l) += I(i,k)*I(j,l);
               dFinterfacedgradJump(i,j,k,l) += _gradFactor*I(i,k)*I(j,l);
@@ -1484,8 +1525,8 @@ void PorosityCohesiveBand3DLaw::transferInterfaceDataToBulk(IPVariable* ipv, con
           }
         }
       }
-
     }
+
   }
 };
 
diff --git a/dG3D/src/dG3DCohesiveBandMaterialLaw.h b/dG3D/src/dG3DCohesiveBandMaterialLaw.h
index 2ac5c2d34..fee88b59d 100644
--- a/dG3D/src/dG3DCohesiveBandMaterialLaw.h
+++ b/dG3D/src/dG3DCohesiveBandMaterialLaw.h
@@ -81,18 +81,25 @@ class NonlocalPorosityCohesiveBand3DLaw : public CohesiveBand3DLaw{
 
 class PorosityCohesiveBand3DLaw : public Cohesive3DLaw{
     #ifndef SWIG
-  protected:
-    double _h;                      // cohesive band thickness
-    double _Kp;               // penalty parameter to avoid penetration
-    bool _accountNormalBulk; //false if the normal part of bulk defo is not accounted for in interfacical defo
-															// true otherwise
-    CoalescenceLaw* _coalescenceLaw; // will be deleted sooon
-    mlawNonLocalPorosity* _mlawPorous;
+protected:
+    // Cohesive band parameters
+    double _h;                        // cohesive band thickness
+    double _Kp;                       // penalty parameter to avoid penetration
+    
+    // Options and numerical parameters
+    bool _accountNormalBulk;          //false if the normal part of bulk defo is not accounted for in interfacical defo
+                                        // true otherwise (==all the components of the bulk F is used...)
+    bool _enabledNegativeJumpN_limit; // true if the normal jump is limited in the negative values
+    double _maxPenetrationjumpN;      // maximal (negative!) value of admissible penetration jump
+    bool _enableDeletionOfBadConditionedIPVs; // true if bad-conditionned ipvs are deleted for convergence purpose
     
+    // Porous law pointer (for mlaw-dependent case)
+    mlawNonLocalPorosity* _mlawPorous;
+    CoalescenceLaw* _coalescenceLaw;  // will be deleted sooon
     #endif
 
   public:
-    PorosityCohesiveBand3DLaw(const int num, const double h, const double Kp, const bool accountNormalBulk);
+    PorosityCohesiveBand3DLaw(const int num, const double h, const double Kp, const bool accountNormalBulk=true);
     void setCoalescenceLaw(const CoalescenceLaw& coalesLaw);
     void setPorousMaterialLaw(const NonLocalPorousDuctileDamageDG3DMaterialLaw& dG3DmlawPorous);
     #ifndef SWIG
-- 
GitLab


From 3f04a7a58df3690f854f8474572104d01d9c836f Mon Sep 17 00:00:00 2001
From: jleclerc <julien.leclerc@ulg.ac.be>
Date: Tue, 29 Jan 2019 11:34:04 +0100
Subject: [PATCH 07/11] add options settings

---
 dG3D/src/dG3DCohesiveBandMaterialLaw.cpp | 25 ++++++++++++++++++++++++
 dG3D/src/dG3DCohesiveBandMaterialLaw.h   |  2 ++
 2 files changed, 27 insertions(+)

diff --git a/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp b/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp
index e4432258b..7b9187129 100644
--- a/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp
+++ b/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp
@@ -1030,6 +1030,31 @@ void PorosityCohesiveBand3DLaw::setPorousMaterialLaw(const NonLocalPorousDuctile
   _coalescenceLaw = _mlawPorous->getCoalescenceLaw()->clone();
 };
 
+void PorosityCohesiveBand3DLaw::setMaximalAdmissiblePenetrationJump(bool activationFlag, double maxNegJump){
+  _enabledNegativeJumpN_limit = activationFlag;
+  _maxPenetrationjumpN = maxNegJump;
+  if (maxNegJump > 0.){
+    Msg::Fatal("PorosityCohesiveBand3DLaw::the admissible value shloud be negative instead of %e",maxNegJump);
+  }
+  
+  if (_enabledNegativeJumpN_limit){
+    Msg::Info("PorosityCohesiveBand3DLaw:: the normal penetration jump is limited to %e",maxNegJump);
+  }
+  else{
+     Msg::Info("PorosityCohesiveBand3DLaw:: the normal penetration jump is not limited");
+  }
+}
+
+void PorosityCohesiveBand3DLaw::setDeletionOfBadConditionnedIPV(bool activationFlag){
+  _enableDeletionOfBadConditionedIPVs = activationFlag;
+  if (_enableDeletionOfBadConditionedIPVs){
+    Msg::Info("PorosityCohesiveBand3DLaw:: the bad-conditionned ipvs will be deleted");
+  }
+  else{
+     Msg::Info("PorosityCohesiveBand3DLaw:: the bad-conditionned ipvs will NOT be deleted");
+  }
+}
+
 
 void PorosityCohesiveBand3DLaw::createIPState(IPStateBase* &ips,const bool* state_,
                                               const MElement *ele, const int nbFF_, const IntPt *GP,
diff --git a/dG3D/src/dG3DCohesiveBandMaterialLaw.h b/dG3D/src/dG3DCohesiveBandMaterialLaw.h
index fee88b59d..964853aa0 100644
--- a/dG3D/src/dG3DCohesiveBandMaterialLaw.h
+++ b/dG3D/src/dG3DCohesiveBandMaterialLaw.h
@@ -102,6 +102,8 @@ protected:
     PorosityCohesiveBand3DLaw(const int num, const double h, const double Kp, const bool accountNormalBulk=true);
     void setCoalescenceLaw(const CoalescenceLaw& coalesLaw);
     void setPorousMaterialLaw(const NonLocalPorousDuctileDamageDG3DMaterialLaw& dG3DmlawPorous);
+    void setMaximalAdmissiblePenetrationJump(bool activationFlag, double maxNegJump=-1.);
+    void setDeletionOfBadConditionnedIPV(bool activationFlag);
     #ifndef SWIG
     PorosityCohesiveBand3DLaw(const PorosityCohesiveBand3DLaw& src);
     virtual ~PorosityCohesiveBand3DLaw(){
-- 
GitLab


From 9589d0f64a8cca6e3a11c5532ab8033069fcfbff Mon Sep 17 00:00:00 2001
From: jleclerc <julien.leclerc@ulg.ac.be>
Date: Tue, 29 Jan 2019 11:35:06 +0100
Subject: [PATCH 08/11] adapt printed message

---
 dG3D/src/dG3DCohesiveBandMaterialLaw.cpp | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp b/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp
index 7b9187129..e6aaedd79 100644
--- a/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp
+++ b/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp
@@ -1514,9 +1514,11 @@ void PorosityCohesiveBand3DLaw::transferInterfaceDataToBulk(IPVariable* ipv, con
       }
     }
     
-    double detF = STensorOperation::determinantSTensor3(F);
-    if (detF < 1.e-10){Msg::Error("PorosityCohesiveBand3DLaw:: Negative jacobian= %e", detF);};
-
+    if(!fipvprev->isDeleted()){
+      double detF = STensorOperation::determinantSTensor3(F);
+      if (detF < 1.e-10){Msg::Error("PorosityCohesiveBand3DLaw:: Negative jacobian= %e", detF);};
+    }
+    
     // Compute derivatives
     if (stiff)
     {
-- 
GitLab


From 0c637c984dd805734cf32c638972dc2e80f6566a Mon Sep 17 00:00:00 2001
From: jleclerc <julien.leclerc@ulg.ac.be>
Date: Fri, 1 Feb 2019 15:54:25 +0100
Subject: [PATCH 09/11] update first version of predictor update

---
 NonLinearSolver/Domain/partDomain.h           |  3 +
 .../nlsolver/implicitHulbertChungPETSc.h      | 26 +++++-
 .../nlsolver/nonLinearMechSolver.cpp          | 93 +++++++++++++++----
 .../nlsolver/nonLinearMechSolver.h            |  2 +
 NonLinearSolver/nlsolver/nonLinearSystems.h   |  1 +
 dG3D/src/dG3DDomain.cpp                       | 44 +++++++++
 dG3D/src/dG3DDomain.h                         |  1 +
 7 files changed, 151 insertions(+), 19 deletions(-)

diff --git a/NonLinearSolver/Domain/partDomain.h b/NonLinearSolver/Domain/partDomain.h
index d347b1ae6..113e49771 100644
--- a/NonLinearSolver/Domain/partDomain.h
+++ b/NonLinearSolver/Domain/partDomain.h
@@ -198,6 +198,9 @@ public:
 
   virtual FunctionSpaceBase* getFunctionSpace() =0;
   virtual const FunctionSpaceBase* getFunctionSpace() const=0;
+  
+  virtual void getDofKeyswithZeroPredictor(std::vector<Dof>& DofVector) const{};
+  
 // some data of BC have to be set by domain
   virtual FunctionSpaceBase* getSpaceForBC(const nonLinearBoundaryCondition::type bc_type, const nonLinearBoundaryCondition::location bc_location,
                                            const int dof_comp,const groupOfElements *groupBC) const=0; // for dirichlet, neumann and initial
diff --git a/NonLinearSolver/nlsolver/implicitHulbertChungPETSc.h b/NonLinearSolver/nlsolver/implicitHulbertChungPETSc.h
index c9b8fd678..03c2126fe 100644
--- a/NonLinearSolver/nlsolver/implicitHulbertChungPETSc.h
+++ b/NonLinearSolver/nlsolver/implicitHulbertChungPETSc.h
@@ -85,10 +85,14 @@ class implicitHulbertChungPetsc : public nonLinearSystem<scalar>,
     bool _whichStep;
     double _alpham,_alphaf, _beta,_gamma; // parameters of scheme
     double _timeStep; // value of time step
-    int _nbRows; // To know the system size.
-
+    int _nbRows; // To know the system size
     implicitStepPetsc<scalar> *_currentStep, *_previousStep, *_initialStep, *_step1, *_step2;
     Mat _M; // mass matrix
+    
+    // Predictor management
+    int _nbZeroPredictorDofs; // number of dof for which the predictor is equal to the previous step
+    std::vector<int> _zeroPredictorDofKeys; // Keys of concerned dof
+    std::vector<scalar> _zeros; // Vectors of zeros
 
     int _numIteration;
     double OneMinusAlphaf; // 1 - alphaf
@@ -110,6 +114,7 @@ class implicitHulbertChungPetsc : public nonLinearSystem<scalar>,
     implicitHulbertChungPetsc(double alpham=0., double alphaf = 1., double beta=0., double gamma=0.5):
           _alphaf(alphaf), _alpham(alpham),_beta(beta),_gamma(gamma),
           linearSystemPETSc<scalar>(),_whichStep(true),_timeStep(0.), _nbRows(0),
+          _nbZeroPredictorDofs(0), _zeroPredictorDofKeys(), _zeros(),
           _numIteration(0),_flagb(false),_initialStep(NULL){
       OneMinusAlphaf = 1. - _alphaf;
       AlphafMinusOne = -1.* OneMinusAlphaf;
@@ -154,6 +159,16 @@ class implicitHulbertChungPetsc : public nonLinearSystem<scalar>,
       implicitStepPetsc<scalar>& dst = this->getState(destination);
       dst = src;
     }
+    
+    virtual void addZeroPredictorDofKeys(const std::vector<int>& dofKeys){
+      // Copy the vector dof keys and set _nbZeroPredictorDofs and _zeros accordingly
+      _zeroPredictorDofKeys.clear();
+      _zeroPredictorDofKeys = dofKeys;
+      _nbZeroPredictorDofs = _zeroPredictorDofKeys.size();
+      _zeros.clear();
+      _zeros.assign(_nbZeroPredictorDofs,0.);
+    }
+    
 
     void nextStep(){
       if(_whichStep){
@@ -463,6 +478,13 @@ class implicitHulbertChungPetsc : public nonLinearSystem<scalar>,
         _try(VecAssemblyBegin(_currentStep->_xddot));
         _try(VecAssemblyEnd(_currentStep->_xddot));
         _try(VecZeroEntries(this->_currentStep->_xddot));
+        
+        if(_nbZeroPredictorDofs > 0){
+          // Reset the predictor state to zero for some of the dofs
+          /* check if this is at the good place */
+          _try(VecSetValues(_currentStep->_xdot , _nbZeroPredictorDofs, &_zeroPredictorDofKeys[0], &_zeros[0], INSERT_VALUES));
+          _try(VecSetValues(_currentStep->_xddot, _nbZeroPredictorDofs, &_zeroPredictorDofKeys[0], &_zeros[0], INSERT_VALUES));
+        }
       }
       else{
         _try(MatAssemblyBegin(this->_a,MAT_FINAL_ASSEMBLY));
diff --git a/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp b/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
index 98c025c48..b4f3bed24 100644
--- a/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
+++ b/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
@@ -3695,6 +3695,12 @@ void nonLinearMechSolver::velocityBC(std::string onwhat, const int numphys, cons
   allDirichlet.back()._mycondition = nonLinearBoundaryCondition::velocity;
 }
 
+void nonLinearMechSolver::accelerationBC(std::string onwhat, const int numphys, const int comp, const double value)
+{
+  this->displacementBC(onwhat,numphys,comp,value);
+  allDirichlet.back()._mycondition = nonLinearBoundaryCondition::acceleration;
+}
+
 void nonLinearMechSolver::displacementBC(std::string onwhat, const int numphys, const int comp, const double value, elementFilter *filter)
 {
   simpleFunction<double> *fct = new simpleFunctionTime<double>(value);
@@ -5798,15 +5804,68 @@ void nonLinearMechSolver::initializeStaticScheme()
           AssembleMass(cdom->getMassTerm(),rcsp,pAssembler);
         }
       }
-    }
+    }    
   }
 
-
   std::string name = "A";
   linearSystem<double>* lsys =pAssembler->getLinearSystem(name);
-  nonLinearSystem<double>* nonsys = dynamic_cast<nonLinearSystem<double>*>(lsys);
+  nonLinearSystem<double>* nlsys = dynamic_cast<nonLinearSystem<double>*>(lsys);
   pathFollowingSystem<double>* pathSys = dynamic_cast<pathFollowingSystem<double>*>(lsys);
 
+  if (whatScheme == Implicit){ //now only for implicit and dgdomain
+  
+    // Collect all dofs number associated with components for which the predictor 
+    // (in the time integration scheme) needs to be zero (== equal to previous state)
+    std::vector<int> DofKeysContainer; // vector of all the considered Dof
+    std::vector<Dof> DofForOneDomain; // vector for one domain
+    
+    /// Loop on all domains
+    for(std::vector<partDomain*>::iterator it=domainVector.begin(); it!=domainVector.end(); ++it){
+      partDomain *dom = *it;
+      // take back the keys of dofs
+      dom->getDofKeyswithZeroPredictor(DofForOneDomain);
+      for(int i = 0; i < DofForOneDomain.size(); ++i){
+        // get the position of the dof and store it in the vector
+        int posDoF = pAssembler->getDofNumber(DofForOneDomain[i]);
+        if(posDoF >= 0){
+          DofKeysContainer.push_back(posDoF);
+        }
+      }
+    }
+    // transfer dof to non-linear system
+    nlsys->addZeroPredictorDofKeys(DofKeysContainer);
+    
+    
+      /*
+      /// Loop on all domains and elements
+      for(std::vector<partDomain*>::iterator it=domainVector.begin(); it!=domainVector.end(); ++it){
+        partDomain *dom = *it;
+        for(groupOfElements::elementContainer::iterator el_it=dom->g_cbegin(); el_it!= dom->g_cend(); ++el_it){
+          MElement* ele = *el_it;
+          /// get all the dof keys of an element
+          std::vector<Dof> R;
+          dom->getFunctionSpace()->getKeys(ele,R);
+          
+          /// Loop on the considered components
+          for(std::vector<int>::iterator comp_it=_noMassMechanicalDofComp.begin(); comp_it!=_noMassMechanicalDofComp.end(); ++comp_it){
+            /// create the filter for the desired components
+            FilterDof* compDofFilter = dom->createFilterDof(*comp_it);/// (maybe do it only once, when the domain is called and create a vector ?)
+            for(int i = 0; i < R.size(); ++i){
+              // if the actual considered dof is associated with the sought component, add its dof number in the container 
+              if(compDofFilter->operator()(R[i])){
+                int posDoF = pAssembler->getDofNumber(R[i]);
+                if(posDoF >= 0){
+                  DofKeysContainer.push_back(posDoF);
+                }
+              }
+            }
+          }
+        }
+      }
+      */
+
+  }
+
   if (_pathFollowing){
     this->computeLoadVector(_ufield);
   }
@@ -6595,7 +6654,7 @@ int nonLinearMechSolver::NewtonRaphson(dofManager<double> *pmanager,unknownField
       break;
     }
   }
-
+  //this->fixNodalDofs(); // reset values of converged state in accordance to BC
   return iter;
 }
 
@@ -12142,20 +12201,20 @@ void nonLinearMechSolver::readPBCDataFromFile(const std::string fn){
 
 
 void nonLinearMechSolver::saveStiffnessMatrixToFile(const int iter){
-linearSystem<double>* lsys=NULL;
-std::string name="A";
-lsys=pAssembler->getLinearSystem(name);
+  linearSystem<double>* lsys=NULL;
+  std::string name="A";
+  lsys=pAssembler->getLinearSystem(name);
 #if defined(HAVE_PETSC)
-linearSystemPETSc<double>* lpet=dynamic_cast<linearSystemPETSc<double>*>(lsys);
-if(lpet==NULL)return;
-else{
-  functionPETSc::MatToFile(lpet->getMatrix(),"data/Mat_Stiffness.txt");
-}
-implicitHulbertChungPetsc<double>* hcsys=dynamic_cast<implicitHulbertChungPetsc<double>*>(lsys);
-if(hcsys==NULL)return;
-else{
-  functionPETSc::MatToFile(hcsys->getMassMatrix(),"data/Mat_Mass.txt");
-}
+  linearSystemPETSc<double>* lpet=dynamic_cast<linearSystemPETSc<double>*>(lsys);
+  if(lpet==NULL)return;
+  else{
+    functionPETSc::MatToFile(lpet->getMatrix(),"data/Mat_Stiffness.txt");
+  }
+  implicitHulbertChungPetsc<double>* hcsys=dynamic_cast<implicitHulbertChungPetsc<double>*>(lsys);
+  if(hcsys==NULL)return;
+  else{
+    functionPETSc::MatToFile(hcsys->getMassMatrix(),"data/Mat_Mass.txt");
+  }
 #endif
 };
 
diff --git a/NonLinearSolver/nlsolver/nonLinearMechSolver.h b/NonLinearSolver/nlsolver/nonLinearMechSolver.h
index 1bc13d9d2..f2d8232f3 100644
--- a/NonLinearSolver/nlsolver/nonLinearMechSolver.h
+++ b/NonLinearSolver/nlsolver/nonLinearMechSolver.h
@@ -358,6 +358,7 @@ class nonLinearMechSolver
 
   selectiveUpdateBase* _selectiveObject;
 
+
 #ifndef SWIG
   void commonModel(); // common operation for loadModel and createModel
   void initMapRanks(const std::set<int>& parts);
@@ -601,6 +602,7 @@ class nonLinearMechSolver
   virtual void velocityBC(std::string onwhat, const int numphys, const int comp, const double value,elementFilter *filter);
   virtual void velocityBC(std::string onwhat, const int numphys, const int comp, simpleFunction<double> *fct);
   virtual void velocityBC(std::string onwhat, const int numphys, const int comp, elementFilter *filter, simpleFunction<double> *fct);
+  virtual void accelerationBC(std::string onwhat, const int numphys, const int comp, const double value);
   virtual void constraintBC(std::string onwhat, const int numphys, const int comp);
 
 	virtual void periodicBC(std::string onwhat, const int phys1, const int phys2, const int v1, const int v2, const int comp);
diff --git a/NonLinearSolver/nlsolver/nonLinearSystems.h b/NonLinearSolver/nlsolver/nonLinearSystems.h
index bf2e26675..6da50cfb5 100644
--- a/NonLinearSolver/nlsolver/nonLinearSystems.h
+++ b/NonLinearSolver/nlsolver/nonLinearSystems.h
@@ -80,6 +80,7 @@ class nonLinearSystem{
   virtual void getFromSolution(int row, scalar &val, const nonLinearBoundaryCondition::whichCondition wv) const=0;
   virtual void getFromMatrix(int row, int col, scalar &val,const typename nonLinearSystem<scalar>::whichMatrix wm) const=0;
   virtual void printAll() {};
+  virtual void addZeroPredictorDofKeys(const std::vector<int>& dofKeys){Msg::Fatal("nonLinearSystem::addZeroPredictorDofKeys should not be called");};
   #if defined(HAVE_MPI)
   virtual double* getArrayIndexPreviousRightHandSidePlusMPI(const int index)=0;
 //  virtual void setPreviousRightHandSidePlusMPI(const int row, const scalar &val)=0;
diff --git a/dG3D/src/dG3DDomain.cpp b/dG3D/src/dG3DDomain.cpp
index c90dcff51..bd96a00df 100644
--- a/dG3D/src/dG3DDomain.cpp
+++ b/dG3D/src/dG3DDomain.cpp
@@ -28,6 +28,7 @@
 #include "nonLocalDamageDG3DIPVariable.h"
 
 
+
 dG3DDomain::dG3DDomain (const int tag,const int phys, const int ws, const int lnum,
                                           const int fdg, const int dim, int nonLocalVar, int nonConstitutiveExtraDOFVar) : dgPartDomain(tag,phys,ws,fdg,nonLocalVar, nonConstitutiveExtraDOFVar),
                                                             _lnum(lnum),
@@ -456,6 +457,8 @@ QuadratureBase* dG3DDomain::getQuadratureRulesForNeumannBC(nonLinearNeumannBC &n
   return integ;
 }
 
+
+
 void dG3DDomain::createTerms(unknownField *uf,IPField*ip)
 {
   FunctionSpace<double>* dgspace = static_cast<FunctionSpace<double>*>(_space);
@@ -693,6 +696,47 @@ BilinearTermBase* dG3DDomain::createNeumannMatrixTerm(FunctionSpace<double> *spn
     return new BiNonLinearTermVoid();
   }
 }
+
+
+
+void dG3DDomain::getDofKeyswithZeroPredictor(std::vector<Dof>& DofVector) const
+{
+  // Send back the dof for which the predictor must be set to zero
+  DofVector.clear(); // clean it first
+  
+  // Create a vector of filters for desired dofs
+  /// /!\ now, it is only non-local dofs
+  std::vector<FilterDof*> FilterVec;
+  for(int i = 0; i < getNumNonLocalVariable(); ++i){
+    FilterVec.push_back(this->createFilterDof(i+3));
+  }
+
+  // Loop on all elements
+  std::vector<Dof> R;
+  const partDomain *dom = static_cast<const partDomain*>(this);
+  for(groupOfElements::elementContainer::const_iterator it=gi->begin(); it!=gi->end();++it){
+    // get dofs of an element
+    MElement *ele = *it;
+    R.clear();
+    _space->getKeys(ele,R);
+    
+    // Loop and all dof and apply all the filters 
+    int nbDof = _space->getNumKeys(ele);
+    for(int i = 0; i < nbDof; ++i){
+      for(std::vector<FilterDof*>::const_iterator it = FilterVec.cbegin(); it != FilterVec.cend(); ++it){
+        FilterDof* filter = *it;
+        if(filter->operator()(R[i])){
+          // if the dof belongs to the desired components, add it to the list and break the filter loop
+          DofVector.push_back(R[i]);
+          break;
+        }
+      }
+    }
+  }
+}
+  
+
+
 void dG3DDomain::computeAllIPStrain(AllIPState *aips,const unknownField *ufield,const IPStateBase::whichState ws, const bool virt){
   IntPt *GP;
   fullVector<double> dispm;
diff --git a/dG3D/src/dG3DDomain.h b/dG3D/src/dG3DDomain.h
index d0804496e..e872885d6 100644
--- a/dG3D/src/dG3DDomain.h
+++ b/dG3D/src/dG3DDomain.h
@@ -178,6 +178,7 @@ class dG3DDomain : public dgPartDomain{
   virtual int getMinusLawNum() const{return _lnum;}
   virtual int getPlusLawNum() const{return _lnum;}
   virtual void setGaussIntegrationRule();
+  virtual void getDofKeyswithZeroPredictor(std::vector<Dof>& DofVector) const; // get back the dof for which the predictor must be set to zero
   virtual void createTerms(unknownField *uf,IPField*ip);
   virtual void initializeTerms(unknownField *uf,IPField *ip);
   virtual LinearTermBase<double>* createNeumannTerm(FunctionSpace<double> *spneu,  const groupOfElements* g,
-- 
GitLab


From 21fdcd5603403c27e479eb6454c234ee1845e7a8 Mon Sep 17 00:00:00 2001
From: jleclerc <julien.leclerc@ulg.ac.be>
Date: Mon, 4 Feb 2019 11:44:14 +0100
Subject: [PATCH 10/11] correct mistake in parenthesis and indentation

---
 dG3D/src/dG3DTerms.cpp | 34 +++++++++++++++-------------------
 1 file changed, 15 insertions(+), 19 deletions(-)

diff --git a/dG3D/src/dG3DTerms.cpp b/dG3D/src/dG3DTerms.cpp
index a7e7c3b4d..14d7a30e1 100644
--- a/dG3D/src/dG3DTerms.cpp
+++ b/dG3D/src/dG3DTerms.cpp
@@ -1257,7 +1257,7 @@ void dG3DForceInter::get(MElement *ele, int npts, IntPt *GP, fullVector<double>
         if( ipvm->getNumberNonLocalVariable()>getNumNonLocalVariable() ) Msg::Fatal("Your minus material law uses more non local variable than your domain dG3DForceInter::get");
         if( ipvp->getNumberNonLocalVariable()>getNumNonLocalVariable() ) Msg::Fatal("Your plus material law uses more non local variable than your domain dG3DForceInter::get");
         if( ipvm->getNumberNonLocalVariable()!= ipvm->getNumberNonLocalVariable()  ) Msg::Fatal("Your minus and plus material laws do not have the same number of non local variables dG3DForceInter::get");
-        if( getNonLocalContinuity() && !broken)
+        if( getNonLocalContinuity())
         {
           // characteristic size
           const double nbetahs = getNonLocalStabilityParameter()/hs;
@@ -1361,7 +1361,7 @@ void dG3DForceInter::get(MElement *ele, int npts, IntPt *GP, fullVector<double>
             Msg::Fatal("Your plus and minus material laws do not use the same number of constitutive extra dof variables dG3DForceInter::get");
         for (int extraDOFField = 0; extraDOFField< ipvp->getNumConstitutiveExtraDofDiffusionVariable(); extraDOFField++){
          for (int extraDOFField2 = 0; extraDOFField2< ipvp->getNumConstitutiveExtraDofDiffusionVariable(); extraDOFField2++){
-          if( getConstitutiveExtraDofDiffusionContinuity() && !broken)
+          if( getConstitutiveExtraDofDiffusionContinuity())
           {
             const double nbetahs = getConstitutiveExtraDofDiffusionStabilityParameter()/hs;
             double fieldJump;
@@ -1494,19 +1494,19 @@ void dG3DForceInter::get(MElement *ele, int npts, IntPt *GP, fullVector<double>
             }
             //new
             if(extraDOFField==0){
-            double gamma=0.;
-            for(int j=0;j<nbFFm;j++)
-            {
-              m(j+indexField*nbFFm) += gamma*( Stiff_alphadialitationm_JumpN*Valsm[j+indexField*nbFFm])*(wJ/2.)*(-1.);
-            }
-            for(int j=0;j<nbFFp;j++)
-            {
-              m(j+indexField*nbFFp+nbdofm) += gamma*( Stiff_alphadialitationp_JumpN*Valsp[j+indexField*nbFFp])*(wJ/2.)*(-1.);
+              double gamma=0.;
+              for(int j=0;j<nbFFm;j++){
+                m(j+indexField*nbFFm) += gamma*( Stiff_alphadialitationm_JumpN*Valsm[j+indexField*nbFFm])*(wJ/2.)*(-1.);
+              }
+              for(int j=0;j<nbFFp;j++){
+                m(j+indexField*nbFFp+nbdofm) += gamma*( Stiff_alphadialitationp_JumpN*Valsp[j+indexField*nbFFp])*(wJ/2.)*(-1.);
+              }
             }
-           }
+          }
+         }
+        }
       }
-      else //broken case
-      {
+      else{ // broken case
         const FractureCohesive3DIPVariable *ipvmf = static_cast<const FractureCohesive3DIPVariable*>(ipvm);
         const FractureCohesive3DIPVariable *ipvpf = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
 
@@ -1529,14 +1529,10 @@ void dG3DForceInter::get(MElement *ele, int npts, IntPt *GP, fullVector<double>
             m(j+k*nbFFp+nbdofm) += interfaceForce[k]*(Valsp[j+0*nbFFp]*wJ);
           }
         }
-       }
       }
-     }
-     }
-     }
     }
-   //}
-    //m.print("Force inter");
+  }
+//m.print("Force inter");
 }
 
 void dG3DStiffnessInter::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &stiff) const
-- 
GitLab


From 046b20159f32d1fbda66eb543d71b3e7a901c0cd Mon Sep 17 00:00:00 2001
From: jleclerc <julien.leclerc@ulg.ac.be>
Date: Mon, 4 Feb 2019 13:41:43 +0100
Subject: [PATCH 11/11] clean code for the selection of code

---
 .../nlsolver/implicitHulbertChungPETSc.h      | 30 +++++++++++++++--
 .../nlsolver/nonLinearMechSolver.cpp          | 33 ++-----------------
 dG3D/src/dG3DDomain.cpp                       | 32 ++++++++++++++----
 3 files changed, 54 insertions(+), 41 deletions(-)

diff --git a/NonLinearSolver/nlsolver/implicitHulbertChungPETSc.h b/NonLinearSolver/nlsolver/implicitHulbertChungPETSc.h
index 03c2126fe..32a91d908 100644
--- a/NonLinearSolver/nlsolver/implicitHulbertChungPETSc.h
+++ b/NonLinearSolver/nlsolver/implicitHulbertChungPETSc.h
@@ -467,6 +467,29 @@ class implicitHulbertChungPetsc : public nonLinearSystem<scalar>,
           _try(VecAssemblyEnd(this->_previousStep->_xdot));
         }
         #endif // HAVE_MPI
+/*
+        // Version with predictor = 0
+        _try(VecCopy(this->_previousStep->_x,this->_currentStep->_x));
+        
+        _try(VecCopy(this->_previousStep->_xdot,this->_currentStep->_xdot));
+        if(_nbZeroPredictorDofs > 0){
+          _try(VecSetValues(_currentStep->_xdot, _nbZeroPredictorDofs, &_zeroPredictorDofKeys[0], &_zeros[0], INSERT_VALUES));
+        }
+        _try(VecAXPY(this->_currentStep->_x,_timeStep,this->_currentStep->_xdot));
+        
+        _try(VecAssemblyBegin(_currentStep->_xddot));
+        _try(VecAssemblyEnd(_currentStep->_xddot));
+        _try(VecCopy(this->_previousStep->_xddot,this->_currentStep->_xddot));
+                if(_nbZeroPredictorDofs > 0){
+          _try(VecSetValues(_currentStep->_xddot, _nbZeroPredictorDofs, &_zeroPredictorDofKeys[0], &_zeros[0], INSERT_VALUES));
+        }
+        _try(VecAXPY(this->_currentStep->_x,Deltat2halfMinusbeta,this->_currentStep->_xddot));
+        _try(VecAXPY(this->_currentStep->_xdot,DeltatOneMinusGamma,this->_currentStep->_xddot));
+        
+        _try(VecZeroEntries(this->_currentStep->_xddot));
+        
+*/
+
 
         _try(VecCopy(this->_previousStep->_x,this->_currentStep->_x));
         _try(VecAXPY(this->_currentStep->_x,_timeStep,this->_previousStep->_xdot));
@@ -478,13 +501,14 @@ class implicitHulbertChungPetsc : public nonLinearSystem<scalar>,
         _try(VecAssemblyBegin(_currentStep->_xddot));
         _try(VecAssemblyEnd(_currentStep->_xddot));
         _try(VecZeroEntries(this->_currentStep->_xddot));
-        
+
+/*
         if(_nbZeroPredictorDofs > 0){
-          // Reset the predictor state to zero for some of the dofs
-          /* check if this is at the good place */
+          // Reset the predictor state to zero for some of the dofs (xdot and xddot)
           _try(VecSetValues(_currentStep->_xdot , _nbZeroPredictorDofs, &_zeroPredictorDofKeys[0], &_zeros[0], INSERT_VALUES));
           _try(VecSetValues(_currentStep->_xddot, _nbZeroPredictorDofs, &_zeroPredictorDofKeys[0], &_zeros[0], INSERT_VALUES));
         }
+*/
       }
       else{
         _try(MatAssemblyBegin(this->_a,MAT_FINAL_ASSEMBLY));
diff --git a/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp b/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
index b4f3bed24..537e11774 100644
--- a/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
+++ b/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
@@ -5813,7 +5813,7 @@ void nonLinearMechSolver::initializeStaticScheme()
   pathFollowingSystem<double>* pathSys = dynamic_cast<pathFollowingSystem<double>*>(lsys);
 
   if (whatScheme == Implicit){ //now only for implicit and dgdomain
-  
+    Msg::Info("start dof checking");
     // Collect all dofs number associated with components for which the predictor 
     // (in the time integration scheme) needs to be zero (== equal to previous state)
     std::vector<int> DofKeysContainer; // vector of all the considered Dof
@@ -5834,36 +5834,7 @@ void nonLinearMechSolver::initializeStaticScheme()
     }
     // transfer dof to non-linear system
     nlsys->addZeroPredictorDofKeys(DofKeysContainer);
-    
-    
-      /*
-      /// Loop on all domains and elements
-      for(std::vector<partDomain*>::iterator it=domainVector.begin(); it!=domainVector.end(); ++it){
-        partDomain *dom = *it;
-        for(groupOfElements::elementContainer::iterator el_it=dom->g_cbegin(); el_it!= dom->g_cend(); ++el_it){
-          MElement* ele = *el_it;
-          /// get all the dof keys of an element
-          std::vector<Dof> R;
-          dom->getFunctionSpace()->getKeys(ele,R);
-          
-          /// Loop on the considered components
-          for(std::vector<int>::iterator comp_it=_noMassMechanicalDofComp.begin(); comp_it!=_noMassMechanicalDofComp.end(); ++comp_it){
-            /// create the filter for the desired components
-            FilterDof* compDofFilter = dom->createFilterDof(*comp_it);/// (maybe do it only once, when the domain is called and create a vector ?)
-            for(int i = 0; i < R.size(); ++i){
-              // if the actual considered dof is associated with the sought component, add its dof number in the container 
-              if(compDofFilter->operator()(R[i])){
-                int posDoF = pAssembler->getDofNumber(R[i]);
-                if(posDoF >= 0){
-                  DofKeysContainer.push_back(posDoF);
-                }
-              }
-            }
-          }
-        }
-      }
-      */
-
+    Msg::Info("Dof checking for predictor done: %i DoF are with a predictor = 0",DofKeysContainer.size());
   }
 
   if (_pathFollowing){
diff --git a/dG3D/src/dG3DDomain.cpp b/dG3D/src/dG3DDomain.cpp
index bd96a00df..7670448da 100644
--- a/dG3D/src/dG3DDomain.cpp
+++ b/dG3D/src/dG3DDomain.cpp
@@ -704,25 +704,25 @@ void dG3DDomain::getDofKeyswithZeroPredictor(std::vector<Dof>& DofVector) const
   // Send back the dof for which the predictor must be set to zero
   DofVector.clear(); // clean it first
   
+/*
   // Create a vector of filters for desired dofs
   /// /!\ now, it is only non-local dofs
   std::vector<FilterDof*> FilterVec;
   for(int i = 0; i < getNumNonLocalVariable(); ++i){
-    FilterVec.push_back(this->createFilterDof(i+3));
+    FilterVec.push_back(this->createFilterDof(getDim()+i));
   }
-
   // Loop on all elements
   std::vector<Dof> R;
   const partDomain *dom = static_cast<const partDomain*>(this);
-  for(groupOfElements::elementContainer::const_iterator it=gi->begin(); it!=gi->end();++it){
+  for(groupOfElements::elementContainer::const_iterator it=g->begin(); it!=g->end();++it){
     // get dofs of an element
     MElement *ele = *it;
     R.clear();
     _space->getKeys(ele,R);
     
-    // Loop and all dof and apply all the filters 
-    int nbDof = _space->getNumKeys(ele);
-    for(int i = 0; i < nbDof; ++i){
+    // Loop and all dof and apply all the filters
+    //int nbDof = _space->getNumKeys(ele);
+    for(int i = 0; i < R.size(); ++i){
       for(std::vector<FilterDof*>::const_iterator it = FilterVec.cbegin(); it != FilterVec.cend(); ++it){
         FilterDof* filter = *it;
         if(filter->operator()(R[i])){
@@ -733,8 +733,26 @@ void dG3DDomain::getDofKeyswithZeroPredictor(std::vector<Dof>& DofVector) const
       }
     }
   }
+*/
+  // Loop on all elements
+  std::vector<Dof> R;
+  const partDomain *dom = static_cast<const partDomain*>(this);
+  for(groupOfElements::elementContainer::const_iterator it=g->begin(); it!=g->end();++it){
+    // get dofs of an element
+    MElement *ele = *it;
+    R.clear();
+    _space->getKeys(ele,R);
+    
+    int nbFF = ele->getNumShapeFunctions();
+    
+    // Loop on all non-local var and shape functions
+    for(int nl = 0; nl < getNumNonLocalVariable(); ++nl){
+      for(int i = 0; i < nbFF; ++i){
+        DofVector.push_back(R[i+(3+nl)*nbFF]);
+      }
+    }
+  }
 }
-  
 
 
 void dG3DDomain::computeAllIPStrain(AllIPState *aips,const unknownField *ufield,const IPStateBase::whichState ws, const bool virt){
-- 
GitLab