Merge branch 'vdg-cm3'

NonLinearSolver/internalPoints/ipCohesive.h

@@ -113,7 +113,7 @@ class IPCohesive : public IPVariable{
 		virtual const double& getConstRefToReversibleEnergy() const {return _reversibleEnergy;};
 		
 		virtual double& getRefToIrreversibleEnergy() {return _irreversibleEnergy;};
-		virtual const double& getConstRefToIrreversibleEnergy() const {return _irreversibleEnergy;};
+		virtual double irreversibleEnergy() const {return _irreversibleEnergy;};
 		
 		virtual double& getRefToEffectivePotential() {return _effectivePotential;};
 		virtual const double& getConstRefToEffectivePotential() const {return _effectivePotential;};

NonLinearSolver/materialLaw/CMakeLists.txt

@@ -59,6 +59,7 @@ set(SRC
 	cohesiveDamageLaw.cpp
   mlawHyperelasticDamage.cpp
 	elasticPotential.cpp
+  cohesiveElasticPotential.cpp
   mlawAnisotropicStoch.cpp
   #  src/op_eshelby.cpp
    # Headers without cpp change this ??

NonLinearSolver/materialLaw/cohesiveDamageLaw.h

@@ -44,7 +44,7 @@ class BilinearCohesiveDamageLaw : public cohesiveDamageLaw{
 															 const double D0, double& D, 
 															 const bool stiff, double& dDdp) const;
 		virtual cohesiveDamageLaw::cohesiveDamageLawName getType() {return cohesiveDamageLaw::Bilinear;};
-		virtual cohesiveDamageLaw* clone() const {new BilinearCohesiveDamageLaw(*this);};
+		virtual cohesiveDamageLaw* clone() const {return new BilinearCohesiveDamageLaw(*this);};
 		#endif // SWIG
 };
 
@@ -61,7 +61,7 @@ class ExponentialCohesiveDamageLaw  : public cohesiveDamageLaw{
 															 const double D0, double& D, 
 															 const bool stiff, double& dDdp) const;
 		virtual cohesiveDamageLaw::cohesiveDamageLawName getType() {return cohesiveDamageLaw::Exponential;};
-		virtual cohesiveDamageLaw* clone() const {new ExponentialCohesiveDamageLaw(*this);};
+		virtual cohesiveDamageLaw* clone() const {return new ExponentialCohesiveDamageLaw(*this);};
 		#endif // SWIG
 };
 

NonLinearSolver/materialLaw/cohesiveElasticPotential.cpp

+//
+// C++ Interface: material law
+//
+// Author:  <Van Dung Nguyen>, (C) 2017
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+
+#include "cohesiveElasticPotential.h"
+
+linearCohesiveElasticPotential::linearCohesiveElasticPotential(const int num): cohesiveElasticPotential(num){};
+
+double linearCohesiveElasticPotential::get(const SVector3& ujump, const SVector3& normal, const double Kn, const double Kt) const{
+  double deltan = dot(ujump,normal);
+  double Psi = 0.5*Kn*deltan*deltan;
+  
+  SVector3 deltat(ujump);
+  for (int i=0; i< 3; i++){
+    deltat(i) -= deltan*normal(i);
+  }
+  
+  Psi += 0.5*Kt*dot(deltat,deltat);
+  
+  return Psi;
+};
+double linearCohesiveElasticPotential::getPositiveValue(const SVector3& ujump, const SVector3& normal, const double Kn, const double Kt) const{
+  double deltan = dot(ujump,normal);
+  double PoPlus = 0.;
+  if (deltan >=0.){
+    PoPlus += 0.5*Kn*deltan*deltan;
+  }
+  SVector3 deltat(ujump);
+  for (int i=0; i< 3; i++){
+    deltat(i) -= deltan*normal(i);
+  }
+  PoPlus += 0.5*Kt*dot(deltat,deltat);
+  return PoPlus;
+};
+double linearCohesiveElasticPotential::getNegativeValue(const SVector3& ujump, const SVector3& normal, const double Kn, const double Kt) const{
+  double deltan = dot(ujump,normal);
+  double PoNeg = 0.;
+  if (deltan < 0.){
+    PoNeg += 0.5*Kn*deltan*deltan;
+  }
+  return PoNeg;
+};
+
+void linearCohesiveElasticPotential::constitutivePositive(const SVector3& ujump, const SVector3& normal, const double Kn, const double Kt, SVector3& T, const bool stiff, STensor3* L) const{
+  double deltan = dot(ujump,normal);
+  SVector3 deltat(ujump);
+  for (int i=0; i< 3; i++){
+    deltat(i) -= deltan*normal(i);
+  }
+  T *= 0.;
+  
+  for (int i=0; i<3; i++){
+    if (deltan >=0.){
+      T(i) += Kn*deltan*normal(i);
+    }
+    T(i) += Kt*deltat(i);
+  }
+  
+  if (stiff){
+    (*L) *= 0.;
+    static STensor3 I(1.);
+    for (int i=0; i<3; i++){
+      for (int j=0; j<3; j++){
+        if (deltan >=0.){
+          (*L)(i,j) += Kn*normal(i)*normal(j);
+        }
+        (*L)(i,j) += Kt*(I(i,j) -normal(i)*normal(j)) ;
+      }
+    }
+  }
+};
+void linearCohesiveElasticPotential::constitutiveNegative(const SVector3& ujump, const SVector3& normal, const double Kn, const double Kt, SVector3& T, const bool stiff, STensor3* L) const{
+    double deltan = dot(ujump,normal);
+  T *= 0.;
+  if (deltan <0.){
+    for (int i=0; i<3; i++){
+      T(i) = Kn*deltan*normal(i);
+    }
+  }
+  
+  if (stiff){
+    (*L) *= 0.;
+    if (deltan <0.){
+      for (int i=0; i<3; i++){
+        for (int j=0; j<3; j++){
+          (*L)(i,j) += Kn*normal(i)*normal(j);
+        }
+      }
+    }
+  }
+};
+
+void linearCohesiveElasticPotential::getEffectiveJump(const SVector3& ujump, const SVector3& normal, const double Kn, const double Kt, double& effJump, const bool stiff, SVector3* deffJumpdUjump) const{
+  double betaSq = Kt/Kn;
+  double deltan = dot(ujump,normal);
+  SVector3 deltat(ujump);
+  for (int i=0; i< 3; i++){
+    deltat(i) -= deltan*normal(i);
+  }
+  
+  effJump = 0.;
+  if (deltan >=0.){
+    effJump = sqrt(deltan*deltan+ betaSq*dot(deltat,deltat));
+  }
+  else{
+    effJump = sqrt(betaSq*dot(deltat,deltat));
+  }
+  
+  if (stiff){
+    (*deffJumpdUjump)*= 0.;
+   
+    for (int i=0; i<3; i++){
+      (*deffJumpdUjump)(i) += betaSq*deltat(i)/effJump;
+      if (deltan >=0.){
+        (*deffJumpdUjump)(i) += deltan*normal(i)/effJump;
+      }
+    }
+  }
+};
+    
\ No newline at end of file

NonLinearSolver/materialLaw/cohesiveElasticPotential.h

+//
+// C++ Interface: material law
+//
+// Author:  <Van Dung Nguyen>, (C) 2017
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+
+#ifndef COHESIVEELASTICPOTENTIAL_H_
+#define COHESIVEELASTICPOTENTIAL_H_
+
+#ifndef SWIG
+#include "STensor3.h"
+#endif //SWIG
+
+class cohesiveElasticPotential{
+  #ifndef SWIG
+  protected:
+    int _num;
+  
+  public:
+    cohesiveElasticPotential(const int num): _num(num){}
+    cohesiveElasticPotential(const cohesiveElasticPotential& src): _num(src._num){}
+    virtual ~cohesiveElasticPotential(){}
+    
+    virtual double get(const SVector3& ujump, const SVector3& normal, const double Kn, const double Kt) const = 0;
+    virtual double getPositiveValue(const SVector3& ujump, const SVector3& normal, const double Kn, const double Kt) const = 0;
+    virtual double getNegativeValue(const SVector3& ujump, const SVector3& normal, const double Kn, const double Kt) const = 0;
+    virtual void constitutivePositive(const SVector3& ujump, const SVector3& normal, const double Kn, const double Kt, SVector3& T, const bool sitff, STensor3* L) const = 0;
+    virtual void constitutiveNegative(const SVector3& ujump, const SVector3& normal, const double Kn, const double Kt, SVector3& T, const bool sitff, STensor3* L) const = 0;
+    virtual void getEffectiveJump(const SVector3& ujump, const SVector3& normal, const double Kn, const double Kt, double& effJump, 
+                                    const bool stiff, SVector3* deffJumpdUjump) const = 0;
+    
+    virtual cohesiveElasticPotential* clone() const = 0;
+  #endif //SWIG
+};
+
+class linearCohesiveElasticPotential : public cohesiveElasticPotential{
+  protected:
+  
+  public:
+    linearCohesiveElasticPotential(const int num);
+    #ifndef SWIG
+    linearCohesiveElasticPotential(const linearCohesiveElasticPotential& src):cohesiveElasticPotential(src) {}
+    virtual ~linearCohesiveElasticPotential(){}
+    
+    virtual double get(const SVector3& ujump, const SVector3& normal, const double Kn, const double Kt) const;
+    virtual double getPositiveValue(const SVector3& ujump, const SVector3& normal, const double Kn, const double Kt) const;
+    virtual double getNegativeValue(const SVector3& ujump, const SVector3& normal, const double Kn, const double Kt) const;
+    virtual void constitutivePositive(const SVector3& ujump, const SVector3& normal, const double Kn, const double Kt, SVector3& T, const bool sitff, STensor3* L) const;
+    virtual void constitutiveNegative(const SVector3& ujump, const SVector3& normal, const double Kn, const double Kt, SVector3& T, const bool sitff, STensor3* L) const;
+    virtual void getEffectiveJump(const SVector3& ujump, const SVector3& normal, const double Kn, const double Kt, double& effJump, 
+                                  const bool stiff, SVector3* deffJumpdUjump) const;
+    
+    virtual cohesiveElasticPotential* clone() const {return new linearCohesiveElasticPotential(*this);};
+    #endif //SWIG
+};
+
+
+#endif //COHESIVEELASTICPOTENTIAL_H_

NonLinearSolver/materialLaw/mlawCohesive.cpp

@@ -18,6 +18,147 @@ CohesiveLawBase& CohesiveLawBase::operator=(const materialLaw &source)
   return *this;
 }
 
+generalElasticDamageCohesiveLaw::generalElasticDamageCohesiveLaw(const int num, const cohesiveElasticPotential& elP, const cohesiveDamageLaw& dam):CohesiveLawBase(num,true){
+  _elasticPotential  =elP.clone();
+  _damageLaw = dam.clone();
+};
+generalElasticDamageCohesiveLaw::generalElasticDamageCohesiveLaw(const generalElasticDamageCohesiveLaw& src): CohesiveLawBase(src){
+  _elasticPotential = NULL;
+  if (src._elasticPotential != NULL){
+    _elasticPotential = src._elasticPotential->clone();
+  }
+  
+  _damageLaw = NULL;
+  if (src._damageLaw != NULL){
+    _damageLaw = src._damageLaw->clone();
+  }
+}
+
+generalElasticDamageCohesiveLaw& generalElasticDamageCohesiveLaw::operator=(const materialLaw &source){
+  CohesiveLawBase::operator=(source);
+	const generalElasticDamageCohesiveLaw* psrc = dynamic_cast<const generalElasticDamageCohesiveLaw*>(&source);
+	if (psrc !=NULL){
+    if (_elasticPotential) {delete _elasticPotential; _elasticPotential = NULL;}
+    if (psrc->_elasticPotential != NULL){
+      _elasticPotential = psrc->_elasticPotential->clone();
+    }
+		if (_damageLaw) {delete _damageLaw; _damageLaw = NULL;};
+		if (psrc->_damageLaw != NULL){
+			_damageLaw = psrc->_damageLaw->clone();
+		}
+	}
+  return *this;
+};
+
+generalElasticDamageCohesiveLaw::~generalElasticDamageCohesiveLaw(){
+  if (_elasticPotential) delete _elasticPotential;
+  if (_damageLaw) delete _damageLaw;
+};
+
+void generalElasticDamageCohesiveLaw::constitutive(const SVector3& jump, const SVector3& normDir,
+                              const IPCohesive* q0, IPCohesive* q, SVector3& T, const bool stiff, STensor3& DTDjump) const{
+
+  // take the same internal state as previous
+  q->operator=(*dynamic_cast<const IPVariable*>(q0));
+  
+  // estimation of current kinematic variables	
+	// normal part
+	double deltan = dot(jump,normDir);
+	if (deltan>=0){
+		q->getRefToTensionFlag() = true;
+	}
+	else{
+		q->getRefToTensionFlag() = false;
+	}
+	
+  const double& beta = q->getConstRefToMixedModeParameter();
+  const double& Kp = q->getConstRefToInitialPenality();
+  
+  double betaSq = beta*beta;
+  double Kn = Kp;
+  double Kt = betaSq*Kp;
+  
+  
+  double energyPos = _elasticPotential->getPositiveValue(jump,normDir,Kn,Kt);
+  double energyNeg = _elasticPotential->getNegativeValue(jump,normDir,Kn,Kt);
+  
+  SVector3 TPos(0.), TNeg(0.);
+  STensor3 dTPosDjump(0.), dTNegDjump(0.);
+  
+
+  _elasticPotential->constitutivePositive(jump,normDir,Kn,Kt,TPos,stiff,&dTPosDjump);
+  _elasticPotential->constitutiveNegative(jump,normDir,Kn,Kt,TNeg,stiff,&dTNegDjump);
+  
+  double delta(0.);
+  SVector3 DdeltaDjump(0.);
+  _elasticPotential->getEffectiveJump(jump,normDir,Kn,Kt,delta,stiff,&DdeltaDjump);
+  
+  // compute maximal effective opening
+	double& deltamax = q->getRefToMaximalEffectiveJump();
+	const double& deltamaxPrev = q0->getConstRefToMaximalEffectiveJump();
+	deltamax = std::max(deltamaxPrev,delta);
+  
+  // compute damage
+	const double& delta0 = q->getConstRefToEffectiveJumpDamageOnset();
+  const double& deltac = q->getConstRefToCriticalEffectiveJump();
+  const double& sigmac = q->getConstRefToCriticalEffectiveStress();
+	
+	double& D = q->getRefToCohesiveDamage();
+	const double& Dprev = q0->getConstRefToCohesiveDamage();
+	double& dDDdelta = q->getRefToDCohesiveDamageDEffectiveJump();
+	
+	D = Dprev;
+	dDDdelta = 0.;
+    
+	_damageLaw->computeDamage(delta0,deltac,deltamaxPrev,deltamax,Dprev,D, stiff, dDDdelta);
+
+	  // parameters need to estimate for current step
+  double& sigmamax = q->getRefToMaximalEffectiveStress();
+	sigmamax = (1-D)*Kp*deltamax;
+	
+	if (D > 0.9999999999){
+		q->Delete(true);
+	}
+	else{
+		q->Delete(false);
+	}
+  
+  T = TPos;
+  T *= (1.-D);
+  T += TNeg;
+  	
+  double& fracEnergy  = q->getRefToFractureEnergy();
+  fracEnergy +=  0.5*deltac*(q0->getConstRefToMaximalEffectiveStress() - sigmamax);  
+	
+  double revEnerg = q->getRefToReversibleEnergy();
+  revEnerg = (1.-D)*energyPos+energyNeg;
+  
+  double & irreEnerg  = q->getRefToIrreversibleEnergy();
+  irreEnerg = q0->irreversibleEnergy();
+	irreEnerg +=  energyPos*(D-Dprev);
+
+  if (stiff){
+    DTDjump = dTPosDjump;
+    DTDjump*= (1.-D);
+    DTDjump += dTNegDjump;
+    for(int i =0; i<3; i++){
+      for(int j=0; j<3; j++){
+        DTDjump(i,j) -= TPos(i)*dDDdelta*DdeltaDjump(j);
+      }
+    }
+    
+    
+		SVector3& dirreEnergDjump = q->getRefToDIrreversibleEnergyDJump();
+		dirreEnergDjump *= 0.;
+
+    for (int i=0; i<3; i++){
+      dirreEnergDjump(i) += (TPos(i)*(D-Dprev) + energyPos*dDDdelta*DdeltaDjump(i));
+    };
+  }
+  
+};
+
+
 BiLinearCohesiveLaw& BiLinearCohesiveLaw::operator=(const materialLaw &source)
 {
   CohesiveLawBase::operator=(source);

NonLinearSolver/materialLaw/mlawCohesive.h

@@ -18,6 +18,7 @@
 #include "STensor33.h"
 #include "mlaw.h"
 #include "cohesiveDamageLaw.h"
+#include "cohesiveElasticPotential.h"
 
 class CohesiveLawBase : public materialLaw{
   protected:
@@ -49,6 +50,28 @@ class CohesiveLawBase : public materialLaw{
     virtual materialLaw* clone() const = 0;
 };
 
+class generalElasticDamageCohesiveLaw : public CohesiveLawBase{
+  protected:
+    cohesiveElasticPotential* _elasticPotential;
+    cohesiveDamageLaw* _damageLaw;
+    
+  public:
+    generalElasticDamageCohesiveLaw(const int num, const cohesiveElasticPotential& elP, const cohesiveDamageLaw& damageLaw);
+    generalElasticDamageCohesiveLaw(const generalElasticDamageCohesiveLaw& src);
+    virtual generalElasticDamageCohesiveLaw& operator=(const materialLaw &source);
+    virtual ~generalElasticDamageCohesiveLaw();
+    
+    virtual void createIPState(IPStateBase* &ips,const bool* state_=NULL,const MElement *ele=NULL, const int nbFF_=0, const IntPt *GP=NULL, const int gpt = 0) const{
+      Msg::Fatal("generalElasticDamageCohesiveLaw::createIPState has not been implemented");
+    };
+
+    virtual void constitutive(const SVector3& jump1, const SVector3& normal1,
+                              const IPCohesive* q0, IPCohesive* q1, SVector3& T, const bool stiff, STensor3& DTDjump) const;
+
+    virtual materialLaw* clone() const {return new generalElasticDamageCohesiveLaw(*this);};
+    virtual bool withEnergyDissipation() const {return true;};
+};
+
 
 // bilinear cohesive law
 class BiLinearCohesiveLaw : public CohesiveLawBase{

NonLinearSolver/nlmechsolpy.i

@@ -34,6 +34,8 @@
   #include "geometryRotation.h"
   #include "nonLinearMicroBC.h"
   #include "elasticPotential.h"
+  #include "cohesiveElasticPotential.h"
+  #include "cohesiveDamageLaw.h"
 %}
 
 %init%{
@@ -74,6 +76,8 @@
 %nodefaultctor GeometryRotation;
 %nodefaultctor nonLinearMicroBC;
 %nodefaultctor elasticPotential;
+%nodefaultctor cohesiveElasticPotential;
+%nodefaultctor cohesiveDamageLaw;
 %nodefaultctor levelSet;
 
 %include "NonLinearSolverConfig.h"
@@ -112,3 +116,5 @@
 %include "geometryRotation.h"
 %include "nonLinearMicroBC.h"
 %include "elasticPotential.h"
+%include "cohesiveElasticPotential.h"
+%include "cohesiveDamageLaw.h"

dG3D/benchmarks/CMakeLists.txt

@@ -121,3 +121,4 @@ add_subdirectory(elasticFE2)
 add_subdirectory(multiscaleCohesiveTest2D_fullDG_microPathFollowing)
 add_subdirectory(multiscaleCohesiveTest2D_mixedMode_compression)
 add_subdirectory(compositeHyperDamage)
+add_subdirectory(mixedmodeDelaminationExpoLaw)

dG3D/benchmarks/mixedmodeDelaminationExpoLaw/CMakeLists.txt

+# test file
+
+set(PYFILE delam.py)
+
+set(FILES2DELETE 
+  disp*.msh  
+  energy.csv
+  force*.csv
+  NodalDisplacement*.csv
+	stress*.msh
+)
+
+add_cm3python_test(${PYFILE} "${FILES2DELETE}")

dG3D/benchmarks/mixedmodeDelaminationExpoLaw/delam.geo

+mm = 1e-3;
+L = 2*mm;
+H = 0.5*mm;
+t = 0.1*mm;
+
+lsca = t;
+
+Point(1) = {0,0,0,lsca};
+Point(2) = {H,0,0,lsca};
+Line(1) = {1, 2};
+Extrude {0, L, 0} {
+  Line{1}; Layers{25}; 
+}
+Extrude {0, H, 0} {
+  Line{2}; Layers{5};
+}
+Extrude {0, 0, -1.5*t} {
+  Surface{5}; Layers{1};
+}
+Extrude {0, 0, t} {
+  Surface{5, 9}; Layers{1};
+}
+Physical Volume(51) = {2, 3};
+Physical Volume(52) = {1};
+Physical Surface(53) = {40, 18};
+Physical Surface(56) = {70};
+Physical Point(58) = {36};

dG3D/benchmarks/mixedmodeDelaminationExpoLaw/delam.py

+#coding-Utf-8-*-
+from gmshpy import *
+from dG3Dpy import*
+
+#script to launch beam problem with a python script
+
+# material law
+lawnumEP = 2 
+rhoEP   = 1000.
+youngEP = 32e9
+nuEP    = 0.3
+sy0EP   = 25.e16
+hEP     = 7.1019
+
+lawcnumInt = 3
+Gc   = 600.
+GcII = 1200.
+sigmac   = 20.e6
+sigmacII = 20.e6
+alphaGc = 1.
+fsmin = 1.
+fsmax = 1.
+Kp = 1e15
+
+# geometry
+meshfile="delam.msh" # name of mesh file
+
+# solver
+sol = 2  # Gmm=0 (default) Taucs=1 PETsc=2
+soltype = 1 # StaticLinear=0 (default) StaticNonLinear=1
+nstep = 200   # number of step (used only if soltype=1)
+ftime =1 # Final time (used only if soltype=1)
+tol=1.e-4   # relative tolerance for NR scheme (used only if soltype=1)
+nstepArch=10 # Number of step between 2 archiving (used only if soltype=1)
+fullDg =0 #O = CG, 1 = DG
+space1 = 0 # function space (Lagrange=0)
+beta1  = 100
+
+#  compute solution and BC (given directly to the solver
+# creation of law
+lawEP = J2LinearDG3DMaterialLaw(lawnumEP,rhoEP,youngEP,nuEP,sy0EP,hEP)
+
+elasPotential = linearCohesiveElasticPotential(1)
+beta = 5e-4
+alp = beta-1
+damageLaw = ExponentialCohesiveDamageLaw(1,alp,beta)
+
+lawInter  = DelaminationLinearCohesive3DLaw(lawcnumInt,elasPotential,damageLaw,Gc,GcII,sigmac,sigmacII,alphaGc,fsmin,fsmax,Kp)
+
+
+# creation of ElasticField
+
+myfieldEP1 = dG3DDomain(10,51,space1,lawnumEP,fullDg)
+myfieldEP1.stabilityParameters(beta1)
+
+myfieldEP2 = dG3DDomain(11,52,space1,lawnumEP,fullDg)
+myfieldEP2.stabilityParameters(beta1)
+
+myinterfield = interDomainBetween3D(12,myfieldEP1,myfieldEP2,lawcnumInt)
+myinterfield.stabilityParameters(beta1)
+
+# creation of Solver
+mysolver = nonLinearMechSolver(1000)
+mysolver.loadModel(meshfile)
+mysolver.addDomain(myfieldEP1)
+mysolver.addDomain(myfieldEP2)
+mysolver.addMaterialLaw(lawEP)
+
+mysolver.addDomain(myinterfield)
+mysolver.addMaterialLaw(lawInter)
+
+
+mysolver.Scheme(soltype)
+mysolver.Solver(sol)
+mysolver.snlData(nstep,ftime,tol)
+mysolver.stepBetweenArchiving(nstepArch)
+
+
+#mysolver.explicitSpectralRadius(ftime,0.5,0.)
+#mysolver.dynamicRelaxation(0.1, ftime, 1.e-3,2)
+#mysolver.explicitTimeStepEvaluation(nstep)
+#mysolver.iterativeProcedure(0)
+
+
+# BC
+#mysolver.displacementBC("Volume",51,1,0.)
+#mysolver.displacementBC("Volume",52,1,0.)
+
+
+mysolver.displacementBC("Face",53,0,0.)
+mysolver.displacementBC("Face",53,1,0.)
+mysolver.displacementBC("Face",53,2,0.)
+
+mysolver.displacementBC("Face",56,2,1e-3)
+
+
+
+mysolver.internalPointBuildView("svm",IPField.SVM, 1, 1)
+mysolver.internalPointBuildView("sig_xx",IPField.SIG_XX, 1, 1)
+mysolver.internalPointBuildView("sig_yy",IPField.SIG_YY, 1, 1)
+mysolver.internalPointBuildView("sig_zz",IPField.SIG_ZZ, 1, 1)
+mysolver.internalPointBuildView("sig_xy",IPField.SIG_XY, 1, 1)
+mysolver.internalPointBuildView("sig_yz",IPField.SIG_YZ, 1, 1)
+mysolver.internalPointBuildView("sig_xz",IPField.SIG_XZ, 1, 1)
+mysolver.internalPointBuildView("epl",IPField.PLASTICSTRAIN, 1, 1)
+
+mysolver.archivingForceOnPhysicalGroup("Face", 53, 2, 1)
+mysolver.archivingForceOnPhysicalGroup("Face", 56, 2, 1)
+
+mysolver.archivingNodeDisplacement(58,2,1)
+
+mysolver.solve()
+
+check = TestCheck()
+check.equal(1.124978e+00,mysolver.getArchivedForceOnPhysicalGroup("Face", 56, 2),2.e-4)
+
+

dG3D/src/dG3DCohesiveIPVariable.h

@@ -232,7 +232,7 @@ class BaseCohesive3DIPVariable : public Cohesive3DIPVariable{
   virtual STensor3& getRefToBulkDeformationGradientAtFailureOnset() {Msg::Error("BaseCohesive3DIPVariable::getRefToBulkDeformationGradientAtFailureOnset is not defined");};
 	
 	
-	virtual double irreversibleEnergy() const {return _qBilinear->getConstRefToIrreversibleEnergy();};
+	virtual double irreversibleEnergy() const {return _qBilinear->irreversibleEnergy();};
 	virtual double& getRefToIrreversibleEnergy() {return _qBilinear->getRefToIrreversibleEnergy();};
 	
 	virtual SVector3& getRefToDIrreversibleEnergyDJump() {return _qBilinear->getRefToDIrreversibleEnergyDJump();};

dG3D/src/dG3DCohesiveMaterialLaw.cpp

@@ -25,7 +25,7 @@ Cohesive3DLaw::Cohesive3DLaw(const Cohesive3DLaw &source) : materialLaw(source)
 
 BaseCohesive3DLaw::BaseCohesive3DLaw(const int num, const double fractureStrengthFactorMin,
                                                const double fractureStrengthFactorMax, const double Kp) :
-									Cohesive3DLaw(num,true),
+									Cohesive3DLaw(num,true), _cohesiveLaw(NULL),
 									_fscmin(fractureStrengthFactorMin),
 									_fscmax(fractureStrengthFactorMax),_Kp(Kp)
 {
@@ -111,6 +111,17 @@ LinearCohesive3DLaw::LinearCohesive3DLaw(const int num, const double Gc, const d
   _cohesiveLaw = new BiLinearCohesiveLaw(num,true);
 }
 
+LinearCohesive3DLaw::LinearCohesive3DLaw(const int num, const cohesiveElasticPotential& poten, const cohesiveDamageLaw& dam, const double Gc, const double sigmac,
+                                               const double beta,const double mu,
+                                               const double fractureStrengthFactorMin,
+                                               const double fractureStrengthFactorMax, const double Kp) :
+									BaseCohesive3DLaw(num,fractureStrengthFactorMin, fractureStrengthFactorMax, Kp),
+									_Gc(Gc), _sigmac(sigmac),_beta(beta), _mu(mu)
+{
+  _cohesiveLaw = new generalElasticDamageCohesiveLaw(num,poten,dam);
+}
+
+
 LinearCohesive3DLaw::LinearCohesive3DLaw(const LinearCohesive3DLaw &source) : BaseCohesive3DLaw(source),
       _Gc(source._Gc), _sigmac(source._sigmac),
       _beta(source._beta), _mu(source._mu)
@@ -531,6 +542,17 @@ NonLocalDamageLinearCohesive3DLaw::NonLocalDamageLinearCohesive3DLaw(const int n
 
 }
 
+NonLocalDamageLinearCohesive3DLaw::NonLocalDamageLinearCohesive3DLaw(const int num, const cohesiveElasticPotential& poten, const cohesiveDamageLaw& dam, const double Gc, const double sigmac, const double Dc,
+                                               const double beta,
+                                               const double fractureStrengthFactorMin,
+                                               const double fractureStrengthFactorMax, const double Kp,
+																							 const bool damageCheck) :
+								LinearCohesive3DLaw(num,poten,dam,Gc,sigmac,beta,0.,fractureStrengthFactorMin,fractureStrengthFactorMax,Kp),
+                                                                                         _Dc(Dc),_damageBased(damageCheck)
+{
+
+}
+
 NonLocalDamageLinearCohesive3DLaw::NonLocalDamageLinearCohesive3DLaw(const NonLocalDamageLinearCohesive3DLaw &source) : LinearCohesive3DLaw(source),
                                                                                         _Dc(source._Dc), _damageBased(source._damageBased)
 {
@@ -716,6 +738,19 @@ DelaminationLinearCohesive3DLaw::DelaminationLinearCohesive3DLaw(const int num,
     alphaGc=_alphaGc;
 }
 
+DelaminationLinearCohesive3DLaw::DelaminationLinearCohesive3DLaw(const int num, const cohesiveElasticPotential& poten, const cohesiveDamageLaw& dam, const double GcI, const double GcII, const double sigmacI,const double sigmacII, const double _alphaGc,
+              const double fractureStrengthFactorMin, const double fractureStrengthFactorMax,
+              const double Kp) :
+	                                          LinearCohesive3DLaw(num,poten,dam,GcI,sigmacI,0,-1.,fractureStrengthFactorMin,fractureStrengthFactorMax,Kp)
+{
+    if(sigmacI <= 0. or sigmacII <= 0. or GcI <= 0. or GcII <= 0. or _alphaGc <= 0.)
+      Msg::Error("DelaminationLinearCohesive3DLaw: wrong inputs");
+    mixedModeSigmaRatio=sigmacII/sigmacI;
+    mixedModeGcRatio=GcII/GcI;
+    alphaGc=_alphaGc;
+}
+
+
 DelaminationLinearCohesive3DLaw::DelaminationLinearCohesive3DLaw(const DelaminationLinearCohesive3DLaw &source) : LinearCohesive3DLaw(source)
 {
     mixedModeSigmaRatio=source.mixedModeSigmaRatio;

dG3D/src/dG3DCohesiveMaterialLaw.h

@@ -109,6 +109,10 @@ class LinearCohesive3DLaw : public BaseCohesive3DLaw
   LinearCohesive3DLaw(const int num, const double Gc, const double sigmac, const double beta, const double mu,
                          const double fractureStrengthFactorMin = 1., const double fractureStrengthFactorMax = 1.,
                          const double Kp= 1e6);
+                         
+  LinearCohesive3DLaw(const int num, const cohesiveElasticPotential& poten, const cohesiveDamageLaw& dam, const double Gc, const double sigmac, const double beta, const double mu,
+                         const double fractureStrengthFactorMin = 1., const double fractureStrengthFactorMax = 1.,
+                         const double Kp= 1e6);
 #ifndef SWIG
   LinearCohesive3DLaw(const LinearCohesive3DLaw &source);
   virtual ~LinearCohesive3DLaw();
@@ -232,6 +236,11 @@ class NonLocalDamageLinearCohesive3DLaw : public LinearCohesive3DLaw
   NonLocalDamageLinearCohesive3DLaw(const int num, const double Gc, const double sigmac, const double Dc, const double beta,
                          const double fractureStrengthFactorMin = 1., const double fractureStrengthFactorMax = 1., const double Kp = 1e6,
 												 const bool damageCheck = false);
+                         
+  NonLocalDamageLinearCohesive3DLaw(const int num, const cohesiveElasticPotential& poten, const cohesiveDamageLaw& dam, const double Gc, const double sigmac, const double Dc, const double beta,
+                         const double fractureStrengthFactorMin = 1., const double fractureStrengthFactorMax = 1., const double Kp = 1e6,
+												 const bool damageCheck = false);
+                         
 #ifndef SWIG
   NonLocalDamageLinearCohesive3DLaw(const NonLocalDamageLinearCohesive3DLaw &source);
   virtual ~NonLocalDamageLinearCohesive3DLaw(){}
@@ -262,6 +271,9 @@ class DelaminationLinearCohesive3DLaw : public LinearCohesive3DLaw
     DelaminationLinearCohesive3DLaw(const int num, const double GcI, const double GcII, const double sigmacI,const double sigmacII, const double alphaGc,
                   const double fractureStrengthFactorMin = 1., const double fractureStrengthFactorMax = 1.,
                   const double Kp = 1e6);
+    DelaminationLinearCohesive3DLaw(const int num, const cohesiveElasticPotential& poten, const cohesiveDamageLaw& dam, const double GcI, const double GcII, const double sigmacI,const double sigmacII, const double alphaGc,
+                  const double fractureStrengthFactorMin = 1., const double fractureStrengthFactorMax = 1.,
+                  const double Kp = 1e6);
 #ifndef SWIG
     DelaminationLinearCohesive3DLaw(const DelaminationLinearCohesive3DLaw &source);
     virtual ~DelaminationLinearCohesive3DLaw(){}