[GenericCrackPhaseField] - Bug Fix

NonLinearSolver/materialLaw/mlawGenericCrackPhaseField.cpp

@@ -139,7 +139,7 @@ bool mlawGenericCrackPhaseField::withEnergyDissipation() const {
     if(_matlaw == NULL) {
         Msg::Error("withEnergyDissipation: Mechanic law is null");
     }
-    return _matlaw->withEnergyDissipation();
+    return true;
 }
 
 void mlawGenericCrackPhaseField::ElasticStiffness(STensor43 *elasticStiffness) const {
@@ -286,15 +286,30 @@ void mlawGenericCrackPhaseField::constitutive(
     // Compare with effenergy of the last step and only update if its greater than it
     const double hh0 = qPF0.getConstRefToeffEner();
     const STensor3 PMax0 = qPF0.getConstRefToPMax();
-    if( hh <= hh0){
+    
+    if( hh <= 0.0){
+        if (hh0 >= 0.0){
             hh = hh0;
             PMax = PMax0;
         }
-    // also update the max values of effEner and PMax if they have increased 
         else{
+            hh = 0.0;
+            STensor3 P_zero;
+            STensorOperation::zero(P_zero);
+            PMax = P_zero;
+        }
+    }
+
+    if( hh <= hh0){
+        hh = hh0;
+        PMax = PMax0;
+        
+    }
+    //update the max values of effEner and PMax
+    
     qPF1.seteffEner(hh);
     qPF1.setPMax(PMax);
-    }
+
 
     // Compute total P
     STensorOperation::scale(Pmeca, degrade, P);
@@ -307,19 +322,24 @@ void mlawGenericCrackPhaseField::constitutive(
     
     double _lc = qPF1.getRefToCharacteristicLength()[0];
 
+    _lc = sqrt(_lc);
+
     // Update the local damage variable
     // qPF1.setlocalDamageVariable(((-2.0 * _lc) / _gc) * (1.0 - nonlocalVariable1) * hh);
-    qPF1.setlocalDamageVariable(((2.0 * _lc) / _gc) * (1.0 - nonlocalVariable1) * hh);
+
+    int test_factor = 1.0;
+
+    qPF1.setlocalDamageVariable(((test_factor * 2.0 * _lc) / _gc) * (1.0 - nonlocalVariable1) * hh);
 
     if(getNbNonLocalVariables() == 1){
 
-        double factor = ((2.0 * _lc) / _gc) * (1.0 - nonlocalVariable1);
+        double factor = ((test_factor * 2.0 * _lc) / _gc) * (1.0 - nonlocalVariable1);
         STensorOperation::scale(PMax, factor, dLocalPlasticStrainDStrain[0]);
 
-        double factor1 = 2.0 * (1.0 - nonlocalVariable1);
+        double factor1 = -2.0 * (1.0 - nonlocalVariable1);
         STensorOperation::scale(Pmeca, factor1, dStressDNonLocalPlasticStrain[0]);
 
-        dLocalPlasticStrainDNonLocalPlasticStrain(0,0) = ((-2.0 * _lc) / _gc) * hh;
+        dLocalPlasticStrainDNonLocalPlasticStrain(0,0) = ((test_factor * -2.0 * _lc) / _gc) * hh;
     }
 
     //TODO: 
@@ -328,13 +348,14 @@ void mlawGenericCrackPhaseField::constitutive(
     // q1->getRefToDefoEnergy() = qMecan->defoEnergy() * degrade * DamageEnergy // TODO Fix 
     // q1->getRefToplastic() = ... compute from mecha
     
-    double defoEnergy = qMecan.defoEnergy() * degrade; 
+    double defoEnergy = test_factor * qMecan.defoEnergy() * degrade; 
     qPF1.setdefoEnergy(defoEnergy);
 
-    double plasticEnergy = qMecan.plasticEnergy();
+    double plasticEnergy = test_factor * qMecan.plasticEnergy();
     qPF1.setplasticEnergy(plasticEnergy);
     
-    double damageEnergy = (1.0 - degrade) * hh;
+    // double damageEnergy = (1.0 - degrade) * hh;
+    double damageEnergy = test_factor * (1.0 - degrade) * hh;
     qPF1.setdamageEnergy(damageEnergy);
     
     // Path following only - irreversible energy
@@ -363,7 +384,7 @@ void mlawGenericCrackPhaseField::constitutive(
 
         double& DIrrevEnergyDNonlocalVar = qPF1.getRefToDIrreversibleEnergyDNonLocalVariable();
 
-        DIrrevEnergyDF = Pmeca;
+        DIrrevEnergyDF = test_factor * Pmeca;
 
         if (this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::DEFO_ENERGY)
         {
@@ -381,15 +402,15 @@ void mlawGenericCrackPhaseField::constitutive(
             //     // DIrrevEnergyDNonlocalVar =  -effEner*q1->getDDamageDp();
             //     DIrrevEnergyDNonlocalVar =  qMecan.defoEnergy() * -2 * (1 - nonlocalVariable1);
             // }
-            DIrrevEnergyDNonlocalVar =  qMecan.defoEnergy() * -2 * (1 - nonlocalVariable1);
+            DIrrevEnergyDNonlocalVar =  test_factor * qMecan.defoEnergy() * -2 * (1 - nonlocalVariable1);
         }
 
         else if ((this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::DISSIPATION_ENERGY) or (this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::DAMAGE_ENERGY))
         {
 
             // DIrrevEnergyDF *= (D - q0->getDamage());
-            DIrrevEnergyDF = (1.0 - degrade) * PMax;
-    
+            DIrrevEnergyDF = test_factor * (1.0 - degrade) * PMax;
+            //DIrrevEnergyDF *= (1.0 - degrade);
             // if (qPF1->getNonLocalToLocal()){
 
             //     DIrrevEnergyDF.daxpy(dLocalEffectiveStrainsDStrains,q1->getDDamageDp()*effEner);
@@ -401,7 +422,7 @@ void mlawGenericCrackPhaseField::constitutive(
             //     DIrrevEnergyDNonlocalVar =  hh * 2.0 * (1.0 - nonlocalVariable1);
 
             // }
-            DIrrevEnergyDNonlocalVar =  hh * 2.0 * (1.0 - nonlocalVariable1);
+            DIrrevEnergyDNonlocalVar = test_factor * hh * 2.0 * (1.0 - nonlocalVariable1);
 
         }
 

dG3D/benchmarks/CrackPhaseField/pf.py

 #coding-Utf-8-*-
 from gmshpy import*
-# from dG3DpyDebug import*
-from dG3Dpy import*
+from dG3DpyDebug import*
+#from dG3Dpy import*
 
 
 
@@ -31,9 +31,7 @@ rho = 2450.0	# kg/m3 :	density
 E = 210.0e9		# Pa : 		Young modulus
 nu = 0.3 		# - : 		poisson coefficient
 
-cl = (0.015*1.0e-3)*(0.015*1.0e-3)#*1.10*1.10	# m2: non-local parameter (== NL-length^2) 
-
-
+cl = (0.0375*1.0e-3) * (0.0375*1.0e-3) # m2: non-local parameter (== NL-length^2) 
 
 
 # # material parameters (2) - cohesive law
@@ -68,7 +66,7 @@ law1.setUseBarF(False)
 
 # FracLaw1 = FractureByCohesive3DLaw(FracLawNum1,BulkLawNum1,InterfaceLawNum1)
 
-mylaw1 = GenericCrackPhaseFieldDG3DMaterialLaw(lawnum3, rho, Cl_Law1, 2700)
+mylaw1 = GenericCrackPhaseFieldDG3DMaterialLaw(lawnum3, rho, Cl_Law1, 2700.0)
 mylaw1.setMechanicalMaterialLaw(law1)
 
 # Solver parameters
@@ -81,7 +79,7 @@ if Solving_mode == 1:
 	ftime =1.0   		# Final time
 	tol=1.e-5  		# Relative tolerance for NR scheme
 	tolAbs = 1.e-6		# Absolute tolerance for NR scheme
-	nstepArch= 4	# Number of step between 2 archiving
+	nstepArch= 50	# Number of step between 2 archiving
 	nstepArchEnergy = 1	# Number of step between 2 energy computation and archiving
 	nstepArchForce = 1	# Number of step between 2 force archiving
 	MaxIter = 200		# Maximum number of iterations
@@ -164,7 +162,6 @@ if Solving_mode == 1 :
 	# solver archiving
 	mysolver.stepBetweenArchiving(nstepArch) 	# archiving frequency for IPField
 	# mysolver.energyComputation(nstepArchEnergy)	# archiving frequency for energy
-
 elif Solving_mode == 2 :	
 	mysolver.Scheme(soltype) 			# solver scheme
 	mysolver.Solver(sol) 				# solver choice
@@ -254,9 +251,16 @@ elif LoadingMode == 1:
 	mysolver.displacementBC("Edge",102,1,0.)
 	mysolver.displacementBC("Edge",104,1,0.)
 
+elif LoadingMode == 2:
+	tot_disp = 0.006 * 1.e-3
 
+	mysolver.displacementBC("Edge",101,0,0.)		# face x = 0
+	mysolver.displacementBC("Edge",101,1,0.)		# face x = 0
 	
+	mysolver.displacementBC("Face",100,2,0.)
 
+	mysolver.displacementBC("Edge",103,1, tot_disp/ftime)
+	mysolver.displacementBC("Edge",103,0, 0.)