Skip to content
Snippets Groups Projects
Select Git revision
  • gmsh_4_12_2
  • master default protected
  • patches-4.14
  • steplayer
  • bl
  • pluginMeshQuality
  • fixBugsAmaury
  • hierarchical-basis
  • alphashapes
  • relaying
  • new_export_boris
  • oras_vs_osm
  • reassign_partitions
  • distributed_fwi
  • rename-classes
  • fix/fortran-api-example-t4
  • robust_partitions
  • reducing_files
  • fix_overlaps
  • 3115-issue-fix
  • 3023-Fillet2D-Update
  • gmsh_4_14_0
  • gmsh_4_13_1
  • gmsh_4_13_0
  • gmsh_4_12_1
  • gmsh_4_12_0
  • gmsh_4_11_1
  • gmsh_4_11_0
  • gmsh_4_10_5
  • gmsh_4_10_4
  • gmsh_4_10_3
  • gmsh_4_10_2
  • gmsh_4_10_1
  • gmsh_4_10_0
  • gmsh_4_9_5
  • gmsh_4_9_4
  • gmsh_4_9_3
  • gmsh_4_9_2
  • gmsh_4_9_1
  • gmsh_4_9_0
40 results

t16.geo

Blame
  • nonLocalDamageDG3DMaterialLaw.cpp 53.74 KiB
    //
    // C++ Interface: materialLaw
    //
    // Description: Class with definition of materialLaw for shell
    //
    //
    // Author:  <Gauthier BECKER>, (C) 2011
    //
    // Copyright: See COPYING file that comes with this distribution
    //
    //
    
    #include "nonLocalDamageDG3DMaterialLaw.h"
    #include "ipstate.h"
    #include "MInterfaceElement.h"
    #include "ipField.h"
    #include <iostream>
    #include <fstream>
    #include <iomanip>
    #include <vector>
    #include <string>
    #include "nonLocalDamageDG3DIPVariable.h"
    #include "dG3DCohesiveIPVariable.h"
    #include "FractureCohesiveDG3DIPVariable.h"
    
    virtualNonLocalDamageDG3DMaterialLaw::virtualNonLocalDamageDG3DMaterialLaw(const int num, const dG3DMaterialLaw& base, const int numNLVariable): 
      dG3DMaterialLaw(num,base.density(),true), numNonlocalVar(numNLVariable)
    {
      elasticStiffness  = (base.elasticStiffness);
      TotalelasticStiffness = (base.TotalelasticStiffness);
      _lawBase = dynamic_cast<dG3DMaterialLaw*>(base.clone());
      characteristicLength.resize(numNLVariable,STensor3(0.));
    };
    
    virtualNonLocalDamageDG3DMaterialLaw::virtualNonLocalDamageDG3DMaterialLaw(const virtualNonLocalDamageDG3DMaterialLaw &source):
      dG3DMaterialLaw(source), numNonlocalVar(source.numNonlocalVar),characteristicLength(source.characteristicLength){
      _lawBase = NULL;
      if (source._lawBase){
        _lawBase = dynamic_cast<dG3DMaterialLaw*>(source._lawBase->clone());
      }
    };
    
    virtualNonLocalDamageDG3DMaterialLaw::~virtualNonLocalDamageDG3DMaterialLaw(){
      if (_lawBase != NULL) delete _lawBase; _lawBase = NULL;
    };
    
    void virtualNonLocalDamageDG3DMaterialLaw::setNonLocalLengthScale(const int index, const CLengthLaw &cLLaw){
      IPCLength* ipCL = NULL;
      cLLaw.createIPVariable(ipCL);
      cLLaw.computeCL(0.,*ipCL);
      characteristicLength[index] = ipCL->getConstRefToCL();
      characteristicLength[index].print("characteristicLength");
      if (ipCL!= NULL) delete ipCL;
    };
    
    void virtualNonLocalDamageDG3DMaterialLaw::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
    {
      // check interface element
      IPStateBase* ipsBase= NULL;
      _lawBase->createIPState(ipsBase,state_,ele,nbFF_,GP,gpt);
      
      dG3DIPVariableBase* dgi = static_cast<dG3DIPVariableBase*>(ipsBase->getState(IPStateBase::initial));
      dG3DIPVariableBase* dg1 = static_cast<dG3DIPVariableBase*>(ipsBase->getState(IPStateBase::previous));
      dG3DIPVariableBase* dg2 = static_cast<dG3DIPVariableBase*>(ipsBase->getState(IPStateBase::current));
      
      
      IPVariable* ipvi = new  virtualNonLocalDamageDG3DIPVariable(numNonlocalVar,dgi,characteristicLength);
      IPVariable* ipv1 = new  virtualNonLocalDamageDG3DIPVariable(numNonlocalVar,dg1,characteristicLength);
      IPVariable* ipv2 = new  virtualNonLocalDamageDG3DIPVariable(numNonlocalVar,dg2,characteristicLength);
      
      if(ips != NULL) delete ips;
      ips = new IP3State(state_,ipvi,ipv1,ipv2);
      delete ipsBase;
    }
    
    void virtualNonLocalDamageDG3DMaterialLaw::createIPVariable(IPVariable* &ipv,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
    {
      IPVariable* ipvBase = NULL;
      _lawBase->createIPVariable(ipvBase,ele,nbFF_,GP,gpt);
      dG3DIPVariableBase* dgipv = static_cast<dG3DIPVariableBase*>(ipvBase);
      if(ipv !=NULL) delete ipv;
      ipv = new  virtualNonLocalDamageDG3DIPVariable(numNonlocalVar,dgipv,characteristicLength);
      delete ipvBase;
    }
    
    void virtualNonLocalDamageDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac)
    {
      /* get ipvariable */
      virtualNonLocalDamageDG3DIPVariable* ipvcur; 
      const virtualNonLocalDamageDG3DIPVariable* ipvprev; 
    
      FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
      if(ipvtmp !=NULL)
      {
    
        ipvcur = static_cast<virtualNonLocalDamageDG3DIPVariable*>(ipvtmp->getIPvBulk());
        const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
        ipvprev = static_cast<const virtualNonLocalDamageDG3DIPVariable*>(ipvtmp2->getIPvBulk());
       }
      else
      {
        ipvcur = static_cast<virtualNonLocalDamageDG3DIPVariable*>(ipv);
        ipvprev = static_cast<const virtualNonLocalDamageDG3DIPVariable*>(ipvp);
      }
      
      IPVariable* dgIPvCur = static_cast<IPVariable*>(ipvcur->getIPVBase());
      const IPVariable* dgIPvPrev = static_cast<const IPVariable*>(ipvprev->getIPVBase());
      
      _lawBase->stress(dgIPvCur,dgIPvPrev,stiff,checkfrac);
    };
    
    
    
    NonLocalDamageDG3DMaterialLaw::NonLocalDamageDG3DMaterialLaw(const int num, const double rho, const char *propName) :
                                                                  dG3DMaterialLaw(num,rho,true)
    {
      _nldlaw = new mlawNonLocalDamage(num,rho,propName);
      double nu = _nldlaw->poissonRatio();
      double mu = _nldlaw->shearModulus();
      double E = 2.*mu*(1.+nu);
      fillElasticStiffness(E, nu, elasticStiffness);
    }
    
    NonLocalDamageDG3DMaterialLaw::~NonLocalDamageDG3DMaterialLaw()
    {
      if (_nldlaw !=NULL) delete _nldlaw;
      _nldlaw = NULL;
    };
    
    
    NonLocalDamageDG3DMaterialLaw::NonLocalDamageDG3DMaterialLaw(const NonLocalDamageDG3DMaterialLaw &source) :
                                                                 dG3DMaterialLaw(source)
    {
    	_nldlaw = NULL;
    	if (source._nldlaw != NULL){
    		_nldlaw = dynamic_cast<materialLawWithDamage*>(source._nldlaw->clone());
    	}
    
    }
    
    void NonLocalDamageDG3DMaterialLaw::setTime(const double t,const double dtime){
      dG3DMaterialLaw::setTime(t,dtime);
      _nldlaw->setTime(t,dtime);
      return;
    }
    
    materialLaw::matname NonLocalDamageDG3DMaterialLaw::getType() const {return _nldlaw->getType();}
    
    bool NonLocalDamageDG3DMaterialLaw::withEnergyDissipation() const {return _nldlaw->withEnergyDissipation();};
    
    void NonLocalDamageDG3DMaterialLaw::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
    {
      // check interface element
      bool inter=true;
      const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
      if(iele==NULL) inter=false;
    
      int nsdv = _nldlaw->getNsdv();
      IPVariable* ipvi = new  nonLocalDamageDG3DIPVariable(nsdv,inter);
      IPVariable* ipv1 = new  nonLocalDamageDG3DIPVariable(nsdv,inter);
      IPVariable* ipv2 = new  nonLocalDamageDG3DIPVariable(nsdv,inter);
    
      if(ips != NULL) delete ips;
      ips = new IP3State(state_,ipvi,ipv1,ipv2);
      _nldlaw->createIPState((static_cast <nonLocalDamageDG3DIPVariable*> (ipvi))->getIPNonLocalDamage(),
                             (static_cast <nonLocalDamageDG3DIPVariable*> (ipv1))->getIPNonLocalDamage(),
                 			 (static_cast <nonLocalDamageDG3DIPVariable*> (ipv2))->getIPNonLocalDamage());
    
    }
    
    void NonLocalDamageDG3DMaterialLaw::createIPVariable(IPVariable* &ipv,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
    {
      if(ipv !=NULL) delete ipv;
      bool inter=true;
      const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
      if(iele == NULL) inter=false;
    
      ipv = new  nonLocalDamageDG3DIPVariable(_nldlaw->getNsdv(),inter);
    
      IPNonLocalDamage * ipvnl = static_cast <nonLocalDamageDG3DIPVariable*>(ipv)->getIPNonLocalDamage();
      _nldlaw->createIPVariable(ipvnl,ele,nbFF_,GP,gpt);
    }
    
    
    double NonLocalDamageDG3DMaterialLaw::soundSpeed() const{return _nldlaw->soundSpeed();}
    
    void NonLocalDamageDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac)
    {
      /* get ipvariable */
      nonLocalDamageDG3DIPVariable* ipvcur; //= static_cast < nonLocalDamageDG3DIPVariable *> (ipv);
      const nonLocalDamageDG3DIPVariable* ipvprev; //= static_cast <const  nonLocalDamageDG3DIPVariable *> (ipvp);
    
      FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
      if(ipvtmp !=NULL)
      {
    
        ipvcur = static_cast<nonLocalDamageDG3DIPVariable*>(ipvtmp->getIPvBulk());
        const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
        ipvprev = static_cast<const nonLocalDamageDG3DIPVariable*>(ipvtmp2->getIPvBulk());
       }
      else
      {
        ipvcur = static_cast<nonLocalDamageDG3DIPVariable*>(ipv);
        ipvprev = static_cast<const nonLocalDamageDG3DIPVariable*>(ipvp);
      }
    
      /* strain xyz */
      const STensor3& F1 = ipvcur->getConstRefToDeformationGradient();
      const STensor3& F0 = ipvprev->getConstRefToDeformationGradient();
    
      /* data for J2 law */
      IPNonLocalDamage* q1 = ipvcur->getIPNonLocalDamage();
      const IPNonLocalDamage* q0 = ipvprev->getIPNonLocalDamage();
    
      /* compute stress */
      _nldlaw->constitutive(F0,F1,ipvcur->getRefToFirstPiolaKirchhoffStress(),q0,q1,ipvcur->getRefToTangentModuli(),
                            ipvcur->getRefToDLocalVariableDStrain()[0],ipvcur->getRefToDStressDNonLocalVariable()[0],
                            ipvcur->getRefToDLocalVariableDNonLocalVariable()(0,0),stiff);
    
      ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
    
    }
    
    double NonLocalDamageDG3DMaterialLaw::scaleFactor() const{return _nldlaw->scaleFactor();};
    
    void NonLocalDamageDG3DMaterialLaw::setMacroSolver(const nonLinearMechSolver* sv){
    	dG3DMaterialLaw::setMacroSolver(sv);
    	if (_nldlaw != NULL){
    		_nldlaw->setMacroSolver(sv);
    	}
    };
    
    //**** NonLocalDamageDG3DMaterialLaw_Stochastic
    NonLocalDamageDG3DMaterialLaw_Stoch::NonLocalDamageDG3DMaterialLaw_Stoch(const int num, const double rho, const char *propName, 
                         const double Ori_x, const double Ori_y, const double Ori_z,const double Lx,const double Ly,const double Lz, 
                         const int dx, const int dy, const int dz, const char *RandProp, const double euler0,const double euler1,
                         const double euler2) : NonLocalDamageDG3DMaterialLaw(num,rho,propName){
      if(_nldlaw != NULL) delete _nldlaw;
      _nldlaw = new mlawNonLocalDamage_Stoch(num,rho,propName,Ori_x, Ori_y, Ori_z, Lx, Ly, Lz, dx, dy, dz, RandProp, 
                                            euler0,euler1,euler2);
      double nu = _nldlaw->poissonRatio();
      double mu = _nldlaw->shearModulus();
      double E = 2.*mu*(1.+nu);
      fillElasticStiffness(E, nu, elasticStiffness);
    }
    
    NonLocalDamageDG3DMaterialLaw_Stoch::NonLocalDamageDG3DMaterialLaw_Stoch(const int num, const double rho, const char *propName, const double Ori_x, const double Ori_y, const double Ori_z,const double Lx,const double Ly,const double Lz, const char *RandProp, const int intpl) : NonLocalDamageDG3DMaterialLaw(num,rho,propName){
      if(_nldlaw != NULL) delete _nldlaw;
      _nldlaw = new mlawNonLocalDamage_Stoch(num,rho,propName,Ori_x, Ori_y, Ori_z, Lx, Ly, Lz, RandProp,intpl);
      double nu = _nldlaw->poissonRatio();
      double mu = _nldlaw->shearModulus();
      double E = 2.*mu*(1.+nu);
      fillElasticStiffness(E, nu, elasticStiffness);
    }
    
    
    
    NonLocalDamageDG3DMaterialLaw_Stoch::~NonLocalDamageDG3DMaterialLaw_Stoch()
    {
    
     if (_nldlaw !=NULL) delete _nldlaw;
      _nldlaw = NULL;
    
    };
    
    
    NonLocalDamageDG3DMaterialLaw_Stoch::NonLocalDamageDG3DMaterialLaw_Stoch(const NonLocalDamageDG3DMaterialLaw_Stoch &source) :
                                                                 NonLocalDamageDG3DMaterialLaw(source) {}
    
    void NonLocalDamageDG3DMaterialLaw_Stoch::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
    {
      // check interface element
      bool inter=true;
      const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
      if(iele==NULL) inter=false;
      SPoint3 pt_Global;
      ele->pnt(GP[gpt].pt[0],GP[gpt].pt[1],GP[gpt].pt[2],pt_Global);
      SVector3 GaussP(pt_Global);
    
      int nsdv = _nldlaw->getNsdv();
      IPVariable* ipvi = new  nonLocalDamageDG3DIPVariable(nsdv,inter);
      IPVariable* ipv1 = new  nonLocalDamageDG3DIPVariable(nsdv,inter);
      IPVariable* ipv2 = new  nonLocalDamageDG3DIPVariable(nsdv,inter);
    
      if(ips != NULL) delete ips;
      ips = new IP3State(state_,ipvi,ipv1,ipv2);
     static_cast<mlawNonLocalDamage_Stoch*>(_nldlaw)->createIPState(GaussP, (static_cast <nonLocalDamageDG3DIPVariable*> (ipvi))->getIPNonLocalDamage(),
                             (static_cast <nonLocalDamageDG3DIPVariable*> (ipv1))->getIPNonLocalDamage(),
                 			 (static_cast <nonLocalDamageDG3DIPVariable*> (ipv2))->getIPNonLocalDamage());
    
    }
    
    void NonLocalDamageDG3DMaterialLaw_Stoch::createIPVariable(IPVariable* &ipv,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
    {
      if(ipv !=NULL) delete ipv;
      bool inter=true;
      const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
      if(iele == NULL) inter=false;
      SPoint3 pt_Global;
      ele->pnt(GP[gpt].pt[0],GP[gpt].pt[1],GP[gpt].pt[2],pt_Global);
      SVector3 GaussP(pt_Global);
    
      ipv = new  nonLocalDamageDG3DIPVariable(_nldlaw->getNsdv(),inter);
    
      IPNonLocalDamage * ipvnl = static_cast <nonLocalDamageDG3DIPVariable*>(ipv)->getIPNonLocalDamage();
      static_cast<mlawNonLocalDamage_Stoch*>(_nldlaw)->createIPVariable(GaussP, ipvnl,ele,nbFF_,GP,gpt);
    }
    // *********
    
    NonLocalDamageHyperelasticDG3DMaterialLaw::NonLocalDamageHyperelasticDG3DMaterialLaw(const int num, const double rho, const elasticPotential& elP, const CLengthLaw &cLLaw,
                                             const DamageLaw &damLaw): dG3DMaterialLaw(num,rho,true){
      _nlLaw = new mlawNonlocalDamageHyperelastic(num,rho,elP,cLLaw,damLaw);
    	
      double E = elP.getYoungModulus();
    	double nu = elP.getPoissonRatio();
    	fillElasticStiffness(E, nu, elasticStiffness); 
    };
    NonLocalDamageHyperelasticDG3DMaterialLaw::NonLocalDamageHyperelasticDG3DMaterialLaw(const NonLocalDamageHyperelasticDG3DMaterialLaw& src):
    dG3DMaterialLaw(src)
    {
      _nlLaw = NULL;
      if (src._nlLaw != NULL){
        _nlLaw = dynamic_cast<mlawNonlocalDamageHyperelastic*>(src._nlLaw->clone());
      }
    };
    void NonLocalDamageHyperelasticDG3DMaterialLaw::setTime(const double t,const double dtime)
    {
    	dG3DMaterialLaw::setTime(t,dtime);
    	_nlLaw->setTime(t,dtime);
    }
    
    void NonLocalDamageHyperelasticDG3DMaterialLaw::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
    {
    	bool inter=true;
    	const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
    	if(iele==NULL) inter=false;
    
    	IPVariable* ipvi = new  nonLocalDamageIsotropicElasticityDG3DIPVariable(*_nlLaw,inter);
    	IPVariable* ipv1 = new  nonLocalDamageIsotropicElasticityDG3DIPVariable(*_nlLaw,inter);
    	IPVariable* ipv2 = new  nonLocalDamageIsotropicElasticityDG3DIPVariable(*_nlLaw,inter);
    	if(ips != NULL) delete ips;
    	ips = new IP3State(state_,ipvi,ipv1,ipv2);
    }
    
    void NonLocalDamageHyperelasticDG3DMaterialLaw::createIPVariable(IPVariable* &ipv,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
    {
    	bool inter=true;
    	const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
    	if(iele == NULL) inter=false;
    
    	if(ipv !=NULL) delete ipv;
    	ipv = new  nonLocalDamageIsotropicElasticityDG3DIPVariable(*_nlLaw,inter);
    }
    
    
    void NonLocalDamageHyperelasticDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac)
    {
        /* get ipvariable */
        nonLocalDamageIsotropicElasticityDG3DIPVariable* ipvcur; //= static_cast < nonLocalDamageDG3DIPVariable *> (ipv);
        const nonLocalDamageIsotropicElasticityDG3DIPVariable* ipvprev; //= static_cast <const  nonLocalDamageDG3DIPVariable *> (ipvp);
    
        FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
        if(ipvtmp !=NULL)
    		{
    				ipvcur = static_cast<nonLocalDamageIsotropicElasticityDG3DIPVariable*>(ipvtmp->getIPvBulk());
    				const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
    				ipvprev = static_cast<const nonLocalDamageIsotropicElasticityDG3DIPVariable*>(ipvtmp2->getIPvBulk());
    		}
        else
    		{
    				ipvcur = static_cast<nonLocalDamageIsotropicElasticityDG3DIPVariable*>(ipv);
    				ipvprev = static_cast<const nonLocalDamageIsotropicElasticityDG3DIPVariable*>(ipvp);
    		}
    
        /* Deformation gradient F */
        const STensor3& F1 = ipvcur->getConstRefToDeformationGradient();
        const STensor3& F0 = ipvprev->getConstRefToDeformationGradient();
    
        /* data for non local isotropic elastic law */
        IPNonLocalDamageIsotropicElasticity* q1 = ipvcur->getIPNonLocalDamageIsotropicElasticity();
        const IPNonLocalDamageIsotropicElasticity* q0 = ipvprev->getIPNonLocalDamageIsotropicElasticity();
    
        /* compute stress */
        _nlLaw->constitutive(F0,F1,ipvcur->getRefToFirstPiolaKirchhoffStress(),q0,q1,ipvcur->getRefToTangentModuli(),
                                     ipvcur->getRefToDLocalVariableDStrain()[0],ipvcur->getRefToDStressDNonLocalVariable()[0],
                                     ipvcur->getRefToDLocalVariableDNonLocalVariable()(0,0),stiff);
    
        ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);   // Is it correct ????
    
    }
    
    NonLocalAnisotropicDamageHyperelasticDG3DMaterialLaw::NonLocalAnisotropicDamageHyperelasticDG3DMaterialLaw(const int num, const double rho, 
                                            const elasticPotential& elP, const CLengthLaw &cLLaw,
                                             const DamageLaw &damLaw): NonLocalDamageHyperelasticDG3DMaterialLaw(num,rho,elP,cLLaw,damLaw){
      if (_nlLaw != NULL) delete _nlLaw;
      _nlLaw = new mlawNonlocalAnisotropicDamageHyperelastic(num,rho,elP,cLLaw,damLaw);                                         
    };
    
    
    NonLocalDamageIsotropicElasticityDG3DMaterialLaw::NonLocalDamageIsotropicElasticityDG3DMaterialLaw(
        const int num,const double rho, double E,const double nu, const CLengthLaw &cLLaw,const DamageLaw &damLaw) :
        dG3DMaterialLaw(num,rho,true)
    {
        _nldIsotropicElasticitylaw = new mlawNonLocalDamageIsotropicElasticity(num,rho,E,nu,cLLaw,damLaw);
        fillElasticStiffness(E, nu, elasticStiffness); // Is it usefull ?
    }
    
    NonLocalDamageIsotropicElasticityDG3DMaterialLaw::~NonLocalDamageIsotropicElasticityDG3DMaterialLaw()
    {
    		if (_nldIsotropicElasticitylaw!=NULL){
    			delete _nldIsotropicElasticitylaw;
    			_nldIsotropicElasticitylaw = NULL;
    		}
    };
    
    NonLocalDamageIsotropicElasticityDG3DMaterialLaw::NonLocalDamageIsotropicElasticityDG3DMaterialLaw(const NonLocalDamageIsotropicElasticityDG3DMaterialLaw &source) :
        dG3DMaterialLaw(source)
    {
    	_nldIsotropicElasticitylaw = NULL;
    	if (source._nldIsotropicElasticitylaw != NULL){
        _nldIsotropicElasticitylaw = new mlawNonLocalDamageIsotropicElasticity(*(source._nldIsotropicElasticitylaw));
    	}
    }
    
    void NonLocalDamageIsotropicElasticityDG3DMaterialLaw::setNonLocalLimitingMethod(const int num){_nldIsotropicElasticitylaw->setNonLocalLimitingMethod(num);};
    
    void NonLocalDamageIsotropicElasticityDG3DMaterialLaw::setTime(const double t,const double dtime)
    {
        dG3DMaterialLaw::setTime(t,dtime);
        _nldIsotropicElasticitylaw->setTime(t,dtime);
    }
    
    materialLaw::matname NonLocalDamageIsotropicElasticityDG3DMaterialLaw::getType() const
    {
        return _nldIsotropicElasticitylaw->getType();
    }
    
    void NonLocalDamageIsotropicElasticityDG3DMaterialLaw::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
    {
        // check if interface element or not
        bool inter=true;
        const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
        if(iele==NULL) inter=false;
    
        IPVariable* ipvi = new  nonLocalDamageIsotropicElasticityDG3DIPVariable(*_nldIsotropicElasticitylaw,inter);
        IPVariable* ipv1 = new  nonLocalDamageIsotropicElasticityDG3DIPVariable(*_nldIsotropicElasticitylaw,inter);
        IPVariable* ipv2 = new  nonLocalDamageIsotropicElasticityDG3DIPVariable(*_nldIsotropicElasticitylaw,inter);
        if(ips != NULL) delete ips;
        ips = new IP3State(state_,ipvi,ipv1,ipv2);
    
        _nldIsotropicElasticitylaw->createIPState((static_cast <nonLocalDamageIsotropicElasticityDG3DIPVariable*> (ipvi))->getIPNonLocalDamageIsotropicElasticity(),
                                      (static_cast <nonLocalDamageIsotropicElasticityDG3DIPVariable*> (ipv1))->getIPNonLocalDamageIsotropicElasticity(),
                                      (static_cast <nonLocalDamageIsotropicElasticityDG3DIPVariable*> (ipv2))->getIPNonLocalDamageIsotropicElasticity());
    
    }
    
    void NonLocalDamageIsotropicElasticityDG3DMaterialLaw::createIPVariable(IPVariable* &ipv,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
    {
        // check if interface element or not
        bool inter=true;
        const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
        if(iele == NULL) inter=false;
    
        if(ipv !=NULL) delete ipv;
        ipv = new  nonLocalDamageIsotropicElasticityDG3DIPVariable(*_nldIsotropicElasticitylaw,inter);
        IPNonLocalDamageIsotropicElasticity * ipvnl = static_cast <nonLocalDamageIsotropicElasticityDG3DIPVariable*>(ipv)->getIPNonLocalDamageIsotropicElasticity();
        _nldIsotropicElasticitylaw->createIPVariable(ipvnl);
    }
    
    double NonLocalDamageIsotropicElasticityDG3DMaterialLaw::soundSpeed() const
    {
        return _nldIsotropicElasticitylaw->soundSpeed();
    }
    
    
    void NonLocalDamageIsotropicElasticityDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac)
    {
        /* get ipvariable */
        nonLocalDamageIsotropicElasticityDG3DIPVariable* ipvcur; //= static_cast < nonLocalDamageDG3DIPVariable *> (ipv);
        const nonLocalDamageIsotropicElasticityDG3DIPVariable* ipvprev; //= static_cast <const  nonLocalDamageDG3DIPVariable *> (ipvp);
    
        FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
        if(ipvtmp !=NULL)
            {
    
                ipvcur = static_cast<nonLocalDamageIsotropicElasticityDG3DIPVariable*>(ipvtmp->getIPvBulk());
                const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
                ipvprev = static_cast<const nonLocalDamageIsotropicElasticityDG3DIPVariable*>(ipvtmp2->getIPvBulk());
            }
        else
            {
                ipvcur = static_cast<nonLocalDamageIsotropicElasticityDG3DIPVariable*>(ipv);
                ipvprev = static_cast<const nonLocalDamageIsotropicElasticityDG3DIPVariable*>(ipvp);
            }
    
        /* Deformation gradient F */
        const STensor3& F1 = ipvcur->getConstRefToDeformationGradient();
        const STensor3& F0 = ipvprev->getConstRefToDeformationGradient();
    
        /* data for non local isotropic elastic law */
        IPNonLocalDamageIsotropicElasticity* q1 = ipvcur->getIPNonLocalDamageIsotropicElasticity();
        const IPNonLocalDamageIsotropicElasticity* q0 = ipvprev->getIPNonLocalDamageIsotropicElasticity();
    
        /* compute stress */
        _nldIsotropicElasticitylaw->constitutive(F0,F1,ipvcur->getRefToFirstPiolaKirchhoffStress(),q0,q1,ipvcur->getRefToTangentModuli(),
                                     ipvcur->getRefToDLocalVariableDStrain()[0],ipvcur->getRefToDStressDNonLocalVariable()[0],
                                     ipvcur->getRefToDLocalVariableDNonLocalVariable()(0,0),stiff);
    
        ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);   // Is it correct ????
    
    }
    
    //
    NonLocalDamageJ2HyperDG3DMaterialLaw::NonLocalDamageJ2HyperDG3DMaterialLaw(const int num,const double rho,
                       double E,const double nu,const J2IsotropicHardening &_j2IH,
                        const CLengthLaw &_cLLaw,const DamageLaw &_damLaw,const double tol,
                        const bool matrixBypert, const double tolPert) :
                                                                  dG3DMaterialLaw(num,rho,true)
    {
      _nldJ2Hyperlaw = new mlawNonLocalDamageJ2Hyper(num,E,nu,rho,_j2IH,_cLLaw,_damLaw,tol,matrixBypert,tolPert);
      fillElasticStiffness(E, nu, elasticStiffness);
    }
    
    NonLocalDamageJ2HyperDG3DMaterialLaw::~NonLocalDamageJ2HyperDG3DMaterialLaw() {
    	if (_nldJ2Hyperlaw!=NULL){
    		delete _nldJ2Hyperlaw;
    		_nldJ2Hyperlaw = NULL;
    	}
    };
    
    
    NonLocalDamageJ2HyperDG3DMaterialLaw::NonLocalDamageJ2HyperDG3DMaterialLaw(const NonLocalDamageJ2HyperDG3DMaterialLaw &source) :
                                                                 dG3DMaterialLaw(source)
    {
    	_nldJ2Hyperlaw = NULL;
    	if (source._nldJ2Hyperlaw != NULL){
    		_nldJ2Hyperlaw = new mlawNonLocalDamageJ2Hyper(*(source._nldJ2Hyperlaw));
    	}
    }
    
    void NonLocalDamageJ2HyperDG3DMaterialLaw::setStrainOrder(const int order){
      _nldJ2Hyperlaw->setStrainOrder(order);
    };
    
    void
    NonLocalDamageJ2HyperDG3DMaterialLaw::setTime(const double t,const double dtime){
      dG3DMaterialLaw::setTime(t,dtime);
      _nldJ2Hyperlaw->setTime(t,dtime);
    }
    
    materialLaw::matname
    NonLocalDamageJ2HyperDG3DMaterialLaw::getType() const {return _nldJ2Hyperlaw->getType();}
    
    void NonLocalDamageJ2HyperDG3DMaterialLaw::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
    {
      // check interface element
      bool inter=true;
      const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
      if(iele==NULL) inter=false;
      IPVariable* ipvi = new  nonLocalDamageJ2HyperDG3DIPVariable(*_nldJ2Hyperlaw,inter);
      IPVariable* ipv1 = new  nonLocalDamageJ2HyperDG3DIPVariable(*_nldJ2Hyperlaw,inter);
      IPVariable* ipv2 = new  nonLocalDamageJ2HyperDG3DIPVariable(*_nldJ2Hyperlaw,inter);
      if(ips != NULL) delete ips;
      ips = new IP3State(state_,ipvi,ipv1,ipv2);
      _nldJ2Hyperlaw->createIPState((static_cast <nonLocalDamageJ2HyperDG3DIPVariable*> (ipvi))->getIPNonLocalDamageJ2Hyper(),
                             (static_cast <nonLocalDamageJ2HyperDG3DIPVariable*> (ipv1))->getIPNonLocalDamageJ2Hyper(),
                 			 (static_cast <nonLocalDamageJ2HyperDG3DIPVariable*> (ipv2))->getIPNonLocalDamageJ2Hyper());
    
    }
    
    void NonLocalDamageJ2HyperDG3DMaterialLaw::createIPVariable(IPVariable* &ipv,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
    {
      if(ipv !=NULL) delete ipv;
      bool inter=true;
      const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
      if(iele == NULL){
        inter=false;
      }
      ipv = new  nonLocalDamageJ2HyperDG3DIPVariable(*_nldJ2Hyperlaw,inter);
      IPNonLocalDamageJ2Hyper * ipvnl = static_cast <nonLocalDamageJ2HyperDG3DIPVariable*>(ipv)->getIPNonLocalDamageJ2Hyper();
      _nldJ2Hyperlaw->createIPVariable(ipvnl);
    }
    
    double
    NonLocalDamageJ2HyperDG3DMaterialLaw::soundSpeed() const{return _nldJ2Hyperlaw->soundSpeed();} // or change value ??
    
    void NonLocalDamageJ2HyperDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac)
    {
      /* get ipvariable */
      nonLocalDamageJ2HyperDG3DIPVariable* ipvcur; //= static_cast < nonLocalDamageDG3DIPVariable *> (ipv);
      const nonLocalDamageJ2HyperDG3DIPVariable* ipvprev; //= static_cast <const  nonLocalDamageDG3DIPVariable *> (ipvp);
    
      FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
      if(ipvtmp !=NULL)
      {
    
        ipvcur = static_cast<nonLocalDamageJ2HyperDG3DIPVariable*>(ipvtmp->getIPvBulk());
        const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
        ipvprev = static_cast<const nonLocalDamageJ2HyperDG3DIPVariable*>(ipvtmp2->getIPvBulk());
       }
      else
      {
        ipvcur = static_cast<nonLocalDamageJ2HyperDG3DIPVariable*>(ipv);
        ipvprev = static_cast<const nonLocalDamageJ2HyperDG3DIPVariable*>(ipvp);
      }
    
      /* strain xyz */
      const STensor3& F1 = ipvcur->getConstRefToDeformationGradient();
      const STensor3& F0 = ipvprev->getConstRefToDeformationGradient();
    
      /* data for J2 law */
      IPNonLocalDamageJ2Hyper* q1 = ipvcur->getIPNonLocalDamageJ2Hyper();
      const IPNonLocalDamageJ2Hyper* q0 = ipvprev->getIPNonLocalDamageJ2Hyper();
      
      STensor43& elasticL = ipvcur->getRefToElasticTangentModuli();
    
      /* compute stress */
      _nldJ2Hyperlaw->constitutive(F0,F1,ipvcur->getRefToFirstPiolaKirchhoffStress(),q0,q1,ipvcur->getRefToTangentModuli(),
                            ipvcur->getRefToDLocalVariableDStrain()[0],ipvcur->getRefToDStressDNonLocalVariable()[0],
                            ipvcur->getRefToDLocalVariableDNonLocalVariable()(0,0),stiff,&elasticL);
    
      ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
    
    }
    
    
    //
    NonLocalDamageGursonDG3DMaterialLaw::NonLocalDamageGursonDG3DMaterialLaw(const int num,const double rho,
                        const double E,const double nu,
                        const double q1, const double q2, const double q3, const double fVinitial,
                        const J2IsotropicHardening &_j2IH,
                        const CLengthLaw &_cLLaw, const double tol,
                        const bool matrixByPert, const double pert) :
                                                                  dG3DMaterialLaw(num,rho,true), _nonLocalMethod(1)
    {
      _nldGursonlaw = new mlawNonLocalDamageGurson(num,E,nu,rho,q1, q2, q3,
                                                  fVinitial,_j2IH,_cLLaw,tol,matrixByPert,pert);
      fillElasticStiffness(E, nu, elasticStiffness);
    }
    
    NonLocalDamageGursonDG3DMaterialLaw::NonLocalDamageGursonDG3DMaterialLaw(const NonLocalDamageGursonDG3DMaterialLaw &source) :
                                                                 dG3DMaterialLaw(source), _nonLocalMethod(source._nonLocalMethod)
    {
    	_nldGursonlaw = NULL;
    	if (source._nldGursonlaw != NULL){
    		_nldGursonlaw = new mlawNonLocalDamageGurson(*(source._nldGursonlaw));
    	}
    }
    
    NonLocalDamageGursonDG3DMaterialLaw::~NonLocalDamageGursonDG3DMaterialLaw(){
    	if (_nldGursonlaw != NULL){
    		delete _nldGursonlaw;
    		_nldGursonlaw = NULL;
    	}
    }
    
    void NonLocalDamageGursonDG3DMaterialLaw::setOrderForLogExp(const int newOrder){_nldGursonlaw->setOrderForLogExp(newOrder);};
    void NonLocalDamageGursonDG3DMaterialLaw::setPredictorCorrectorParameters(const int InitialGuessMethod, const int maxite, const double theta,  const int sys)
      {_nldGursonlaw->setPredictorCorrectorParameters(InitialGuessMethod, maxite, theta, sys);};
    void NonLocalDamageGursonDG3DMaterialLaw::setViscosityLaw(const viscosityLaw& visco)
    {
      _nldGursonlaw->setViscosityLaw(visco);
    };
    void NonLocalDamageGursonDG3DMaterialLaw::setNucleationLaw(const NucleationLaw& added_GdnLaw)
    {
      _nldGursonlaw->setNucleationLaw(added_GdnLaw);
    };
    void NonLocalDamageGursonDG3DMaterialLaw::setScatterredInitialPorosity(const double fmin, const double fmax)
      {
        _nldGursonlaw->setScatterredInitialPorosity(fmin,fmax);
    };
    void NonLocalDamageGursonDG3DMaterialLaw::setCoalescenceLaw(const CoalescenceLaw& added_coalesLaw)
      {_nldGursonlaw->setCoalescenceLaw(added_coalesLaw);};
    
    void NonLocalDamageGursonDG3DMaterialLaw::setFailureTolerance(const double tol, const double elTol, const double localF){
      _nldGursonlaw->setFailureTolerance(tol,elTol,localF);
    };
    
    void NonLocalDamageGursonDG3DMaterialLaw::setSubStepping(const bool fl, const int maxNumStep){
    	_nldGursonlaw->setSubStepping(fl,maxNumStep);
    }
    
    void NonLocalDamageGursonDG3DMaterialLaw::setCorrectedRegularizedFunction(const scalarFunction& fct){
      _nldGursonlaw->setCorrectedRegularizedFunction(fct);
    };
    
    void NonLocalDamageGursonDG3DMaterialLaw::setLocalRegularizedFunction(const scalarFunction& fct){
      _nldGursonlaw->setLocalRegularizedFunction(fct);
    };
    
    void NonLocalDamageGursonDG3DMaterialLaw::setNonLocalMethod(const int i) {
      _nonLocalMethod = i;
      if (_nonLocalMethod == 1)
        Msg::Info("nonlocal porosity is used");
      else if (_nonLocalMethod == 2)
        Msg::Info("nonlocal equivalent plasticity is used");
      else if (_nonLocalMethod == 3)
        Msg::Info("nonlocal volumetric plasticity is used");
      else{
        Msg::Error("nonlocal method %d has not been defined in NonLocalDamageGursonDG3DMaterialLaw; use nonlocal porosity instead",_nonLocalMethod);
        _nonLocalMethod = 1;
      }
    };
    
    void NonLocalDamageGursonDG3DMaterialLaw::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
    {
      // check interface element
      bool inter=true;
      const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
      if(iele==NULL) inter=false;
      double fvInit = _nldGursonlaw->getInitialPorosity();
      IPVariable* ipvi = NULL;
      if (_nonLocalMethod == 1)
        ipvi = new  nonLocalPorosityDG3DIPVariable(*_nldGursonlaw,fvInit,inter);
      else if (_nonLocalMethod == 2)
        ipvi = new  nonLocalPorousEqPlasticDG3DIPVariable(*_nldGursonlaw,fvInit,inter);
      else if (_nonLocalMethod == 3)
        ipvi = new  nonLocalPorousVolumetricPlasticDG3DIPVariable(*_nldGursonlaw,fvInit,inter); 
      else
        Msg::Fatal("nonlocal method must be correctly defined in NonLocalDamageGursonDG3DMaterialLaw");
        
      IPVariable* ipv1= NULL;
      if (_nonLocalMethod == 1)
        ipv1 = new  nonLocalPorosityDG3DIPVariable(*_nldGursonlaw,fvInit,inter);
      else if (_nonLocalMethod == 2)
        ipv1 = new nonLocalPorousEqPlasticDG3DIPVariable(*_nldGursonlaw,fvInit,inter);
      else if (_nonLocalMethod == 3)
        ipv1 = new nonLocalPorousVolumetricPlasticDG3DIPVariable(*_nldGursonlaw,fvInit,inter);
      else
        Msg::Fatal("nonlocal method must be correctly defined in NonLocalDamageGursonDG3DMaterialLaw");
        
      IPVariable* ipv2 = NULL;
      if (_nonLocalMethod == 1)
        ipv2 = new  nonLocalPorosityDG3DIPVariable(*_nldGursonlaw,fvInit,inter);
      else if (_nonLocalMethod == 2)
        ipv2 = new  nonLocalPorousEqPlasticDG3DIPVariable(*_nldGursonlaw,fvInit,inter);
      else if (_nonLocalMethod == 3)
        ipv2 = new nonLocalPorousVolumetricPlasticDG3DIPVariable(*_nldGursonlaw,fvInit,inter);
      else
        Msg::Fatal("nonlocal method must be correctly defined in NonLocalDamageGursonDG3DMaterialLaw");
        
      if(ips != NULL) delete ips;
      ips = new IP3State(state_,ipvi,ipv1,ipv2);
      _nldGursonlaw->createIPState((static_cast <nonLocalPorosityDG3DIPVariable*> (ipvi))->getIPNonLocalPorosity(),
                             (static_cast <nonLocalPorosityDG3DIPVariable*> (ipv1))->getIPNonLocalPorosity(),
                 			 (static_cast <nonLocalPorosityDG3DIPVariable*> (ipv2))->getIPNonLocalPorosity());
    
    }
    
    void NonLocalDamageGursonDG3DMaterialLaw::createIPVariable(IPVariable* &ipv,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
    {
      if(ipv !=NULL) delete ipv;
      bool inter=true;
      const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
      if(iele == NULL){
        inter=false;
      }
      double fvInit = _nldGursonlaw->getInitialPorosity();
      if (_nonLocalMethod == 1)
        ipv = new  nonLocalPorosityDG3DIPVariable(*_nldGursonlaw,fvInit,inter);
      else if (_nonLocalMethod == 2)
        ipv = new  nonLocalPorousEqPlasticDG3DIPVariable(*_nldGursonlaw,fvInit,inter);
      else if (_nonLocalMethod == 3)
        ipv = new  nonLocalPorousVolumetricPlasticDG3DIPVariable(*_nldGursonlaw,fvInit,inter);
      else
        Msg::Fatal("nonlocal method must be correctly defined in NonLocalDamageGursonDG3DMaterialLaw");
        
      IPNonLocalPorosity * ipvnl = static_cast <nonLocalPorosityDG3DIPVariable*>(ipv)->getIPNonLocalPorosity();
      _nldGursonlaw->createIPVariable(ipvnl,ele,nbFF_,GP,gpt);
    }
    
    
    
    void NonLocalDamageGursonDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac)
    {
      /* get ipvariable */
      nonLocalPorosityDG3DIPVariable* ipvcur; //= static_cast < nonLocalDamageDG3DIPVariable *> (ipv);
      const nonLocalPorosityDG3DIPVariable* ipvprev; //= static_cast <const  nonLocalDamageDG3DIPVariable *> (ipvp);
    
      FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
      if(ipvtmp !=NULL)
      {
    
        ipvcur = static_cast<nonLocalPorosityDG3DIPVariable*>(ipvtmp->getIPvBulk());
        const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
        ipvprev = static_cast<const nonLocalPorosityDG3DIPVariable*>(ipvtmp2->getIPvBulk());
       }
      else
      {
        ipvcur = static_cast<nonLocalPorosityDG3DIPVariable*>(ipv);
        ipvprev = static_cast<const nonLocalPorosityDG3DIPVariable*>(ipvp);
      }
    
      /* strain xyz */
      const STensor3& F1 = ipvcur->getConstRefToDeformationGradient();
      const STensor3& F0 = ipvprev->getConstRefToDeformationGradient();
    
      /* data for J2 law */
      IPNonLocalPorosity* q1 = ipvcur->getIPNonLocalPorosity();
      const IPNonLocalPorosity* q0 = ipvprev->getIPNonLocalPorosity();
    
      // set additional kinematic variables
    	q1->setRefToCurrentOutwardNormal(ipvcur->getConstRefToCurrentOutwardNormal());
    	q1->setRefToReferenceOutwardNormal(ipvcur->getConstRefToReferenceOutwardNormal());
    
    
    
      /* compute stress */
      if (_nonLocalMethod == 1)
        _nldGursonlaw->constitutive(F0,F1,ipvcur->getRefToFirstPiolaKirchhoffStress(),q0,q1,ipvcur->getRefToTangentModuli(),
                            ipvcur->getRefToDLocalVariableDStrain()[0],ipvcur->getRefToDStressDNonLocalVariable()[0],
                            ipvcur->getRefToDLocalVariableDNonLocalVariable()(0,0),stiff);
      else if (_nonLocalMethod == 2)
        _nldGursonlaw->constitutive_NonLocalEqPlastic(F0,F1,ipvcur->getRefToFirstPiolaKirchhoffStress(),q0,q1,ipvcur->getRefToTangentModuli(),
                            ipvcur->getRefToDLocalVariableDStrain()[0],ipvcur->getRefToDStressDNonLocalVariable()[0],
                            ipvcur->getRefToDLocalVariableDNonLocalVariable()(0,0),stiff);
      else if (_nonLocalMethod == 3)
        _nldGursonlaw->constitutive_NonLocalVolumetricPlastic(F0,F1,ipvcur->getRefToFirstPiolaKirchhoffStress(),q0,q1,ipvcur->getRefToTangentModuli(),
                            ipvcur->getRefToDLocalVariableDStrain()[0],ipvcur->getRefToDStressDNonLocalVariable()[0],
                            ipvcur->getRefToDLocalVariableDNonLocalVariable()(0,0),stiff);
      else
        Msg::Fatal("nonlocal method must be correctly defined in NonLocalDamageGursonDG3DMaterialLaw");
    
      ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
    
    };
    
    
    
    
    
    
    //
    NonLocalPorousThomasonDG3DMaterialLaw::NonLocalPorousThomasonDG3DMaterialLaw(const int num,const double E,const double nu, const double rho, 
                          const double fVinitial, const double lambda0, const J2IsotropicHardening &j2IH, const CLengthLaw &cLLaw,
                          const double tol, const bool matrixByPert, const double pert):
                                                                  dG3DMaterialLaw(num,rho,true)
    {
      _nlpthomlaw = new mlawNonLocalPorousThomasonLaw(num,E,nu,rho,fVinitial,lambda0,j2IH,cLLaw,tol,matrixByPert,pert);
      fillElasticStiffness(E, nu, elasticStiffness);
    }
    
    NonLocalPorousThomasonDG3DMaterialLaw::NonLocalPorousThomasonDG3DMaterialLaw(const NonLocalPorousThomasonDG3DMaterialLaw &source) :
                                                                 dG3DMaterialLaw(source)
    {
    	_nlpthomlaw = NULL;
    	if (source._nlpthomlaw != NULL){
    		_nlpthomlaw = new mlawNonLocalPorousThomasonLaw(*(source._nlpthomlaw));
    	}
    }
    
    NonLocalPorousThomasonDG3DMaterialLaw::~NonLocalPorousThomasonDG3DMaterialLaw(){
    	if (_nlpthomlaw != NULL){
    		delete _nlpthomlaw;
    		_nlpthomlaw = NULL;
    	}
    }
    
    void NonLocalPorousThomasonDG3DMaterialLaw::setOrderForLogExp(const int newOrder){_nlpthomlaw->setOrderForLogExp(newOrder);};
    void NonLocalPorousThomasonDG3DMaterialLaw::setPredictorCorrectorParameters(const int InitialGuessMethod, const int maxite, const double theta)
      {_nlpthomlaw->setPredictorCorrectorParameters(InitialGuessMethod, maxite, theta);};
    void NonLocalPorousThomasonDG3DMaterialLaw::setViscosityLaw(const viscosityLaw& visco)
    {
      _nlpthomlaw->setViscosityLaw(visco);
    };
    void NonLocalPorousThomasonDG3DMaterialLaw::setNucleationLaw(const NucleationLaw& added_GdnLaw)
    {
      _nlpthomlaw->setNucleationLaw(added_GdnLaw);
    };
    void NonLocalPorousThomasonDG3DMaterialLaw::setScatterredInitialPorosity(const double fmin, const double fmax)
      {
        _nlpthomlaw->setScatterredInitialPorosity(fmin,fmax);
    };
    void NonLocalPorousThomasonDG3DMaterialLaw::setCoalescenceLaw(const CoalescenceLaw& added_coalesLaw)
    {
      _nlpthomlaw->setCoalescenceLaw(added_coalesLaw);
    };
      
    void NonLocalPorousThomasonDG3DMaterialLaw::setYieldSurfaceExponent(const double newN)
    {
      _nlpthomlaw->setYieldSurfaceExponent(newN);
    };
    
    void NonLocalPorousThomasonDG3DMaterialLaw::setSubStepping(const bool fl, const int maxNumStep){
      _nlpthomlaw->setSubStepping(fl,maxNumStep);
    }
    
    
    void NonLocalPorousThomasonDG3DMaterialLaw::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
    {
      // check interface element
      bool inter=true;
      const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
      if(iele==NULL) inter=false;
      double fvInit = _nlpthomlaw->getInitialPorosity();
      IPVariable* ipvi = new  nonLocalPorosityDG3DIPVariable(*_nlpthomlaw,fvInit,inter);
      IPVariable* ipv1 = new  nonLocalPorosityDG3DIPVariable(*_nlpthomlaw,fvInit,inter);
      IPVariable* ipv2 = new  nonLocalPorosityDG3DIPVariable(*_nlpthomlaw,fvInit,inter);
      if(ips != NULL) delete ips;
      ips = new IP3State(state_,ipvi,ipv1,ipv2);
      _nlpthomlaw->createIPState((static_cast <nonLocalPorosityDG3DIPVariable*> (ipvi))->getIPNonLocalPorosity(),
                             (static_cast <nonLocalPorosityDG3DIPVariable*> (ipv1))->getIPNonLocalPorosity(),
                 			 (static_cast <nonLocalPorosityDG3DIPVariable*> (ipv2))->getIPNonLocalPorosity());
    
    }
    
    
    void NonLocalPorousThomasonDG3DMaterialLaw::createIPVariable(IPVariable* &ipv,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
    {
      if(ipv !=NULL) delete ipv;
      bool inter=true;
      const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
      if(iele == NULL){
        inter=false;
      }
      double fvInit = _nlpthomlaw->getInitialPorosity();
      ipv = new  nonLocalPorosityDG3DIPVariable(*_nlpthomlaw,fvInit,inter);
      IPNonLocalPorosity * ipvnl = static_cast <nonLocalPorosityDG3DIPVariable*>(ipv)->getIPNonLocalPorosity();
      _nlpthomlaw->createIPVariable(ipvnl,ele,nbFF_,GP,gpt);
    }
    
    
    
    void NonLocalPorousThomasonDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac)
    {
      /* get ipvariable */
      nonLocalPorosityDG3DIPVariable* ipvcur; //= static_cast < nonLocalDamageDG3DIPVariable *> (ipv);
      const nonLocalPorosityDG3DIPVariable* ipvprev; //= static_cast <const  nonLocalDamageDG3DIPVariable *> (ipvp);
    
      FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
      if(ipvtmp !=NULL)
      {
    
        ipvcur = static_cast<nonLocalPorosityDG3DIPVariable*>(ipvtmp->getIPvBulk());
        const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
        ipvprev = static_cast<const nonLocalPorosityDG3DIPVariable*>(ipvtmp2->getIPvBulk());
       }
      else
      {
        ipvcur = static_cast<nonLocalPorosityDG3DIPVariable*>(ipv);
        ipvprev = static_cast<const nonLocalPorosityDG3DIPVariable*>(ipvp);
      }
    
      /* strain xyz */
      const STensor3& F1 = ipvcur->getConstRefToDeformationGradient();
      const STensor3& F0 = ipvprev->getConstRefToDeformationGradient();
    
      /* data for J2 law */
      IPNonLocalPorosity* q1 = ipvcur->getIPNonLocalPorosity();
      const IPNonLocalPorosity* q0 = ipvprev->getIPNonLocalPorosity();
    
    
      /* compute stress */
      _nlpthomlaw->constitutive(F0,F1,ipvcur->getRefToFirstPiolaKirchhoffStress(),q0,q1,ipvcur->getRefToTangentModuli(),
                            ipvcur->getRefToDLocalVariableDStrain()[0],ipvcur->getRefToDStressDNonLocalVariable()[0],
                            ipvcur->getRefToDLocalVariableDNonLocalVariable()(0,0),stiff);
    
      ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
    
    };
    
    
    
    
    
    
    
    
    
    
    
    //
    NonLocalPorousCoupledDG3DMaterialLaw::NonLocalPorousCoupledDG3DMaterialLaw(const int num,const double E,const double nu, const double rho, 
                          const double q1, const double q2, const double q3,
                          const double fVinitial, const double lambda0, const J2IsotropicHardening &j2IH, const CLengthLaw &cLLaw,
                          const double tol, const bool matrixbyPerturbation, const double pert):
                                    dG3DMaterialLaw(num,rho,true)
    {
      _nlpcoupledlaw = new mlawNonLocalPorousCoupledLaw(num, E, nu, rho, q1, q2, q3, fVinitial, lambda0,
                                                        j2IH, cLLaw, tol, matrixbyPerturbation, pert);
      fillElasticStiffness(E, nu, elasticStiffness);
    }
    
    NonLocalPorousCoupledDG3DMaterialLaw::NonLocalPorousCoupledDG3DMaterialLaw(const NonLocalPorousCoupledDG3DMaterialLaw &source) :
                                                                 dG3DMaterialLaw(source)
    {
    	_nlpcoupledlaw = NULL;
    	if (source._nlpcoupledlaw != NULL){
    		_nlpcoupledlaw = new mlawNonLocalPorousCoupledLaw(*(source._nlpcoupledlaw));
    	}
    }
    
    NonLocalPorousCoupledDG3DMaterialLaw::~NonLocalPorousCoupledDG3DMaterialLaw(){
    	if (_nlpcoupledlaw != NULL){
    		delete _nlpcoupledlaw;
    		_nlpcoupledlaw = NULL;
    	}
    }
    
    void NonLocalPorousCoupledDG3DMaterialLaw::setOrderForLogExp(const int newOrder){_nlpcoupledlaw->setOrderForLogExp(newOrder);};
    void NonLocalPorousCoupledDG3DMaterialLaw::setPredictorCorrectorParameters(const int InitialGuessMethod, const int maxite, const double theta)
      {_nlpcoupledlaw->setPredictorCorrectorParameters(InitialGuessMethod, maxite, theta);};
    void NonLocalPorousCoupledDG3DMaterialLaw::setViscosityLaw(const viscosityLaw& visco)
    {
      _nlpcoupledlaw->setViscosityLaw(visco);
    };
    void NonLocalPorousCoupledDG3DMaterialLaw::setNucleationLaw(const NucleationLaw& added_GdnLaw)
    {
      _nlpcoupledlaw->setNucleationLaw(added_GdnLaw);
    };
    void NonLocalPorousCoupledDG3DMaterialLaw::setScatterredInitialPorosity(const double fmin, const double fmax)
      {
        _nlpcoupledlaw->setScatterredInitialPorosity(fmin,fmax);
    };
    void NonLocalPorousCoupledDG3DMaterialLaw::setCoalescenceLaw(const CoalescenceLaw& added_coalesLaw)
      {_nlpcoupledlaw->setCoalescenceLaw(added_coalesLaw);};
      
    void NonLocalPorousCoupledDG3DMaterialLaw::setYieldSurfaceExponent(const double newN)
    {
      _nlpcoupledlaw->setYieldSurfaceExponent(newN);
    };
    
    void NonLocalPorousCoupledDG3DMaterialLaw::setFailureTolerance(const double tol, const double elTol, const double localF){
      _nlpcoupledlaw->setFailureTolerance(tol,elTol,localF);
    };
    
    void NonLocalPorousCoupledDG3DMaterialLaw::setSubStepping(const bool fl, const int maxNumStep){
      _nlpcoupledlaw->setSubStepping(fl,maxNumStep);
    }
    
    
    void NonLocalPorousCoupledDG3DMaterialLaw::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
    {
      // check interface element
      bool inter=true;
      const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
      if(iele==NULL) inter=false;
      double fvInit = _nlpcoupledlaw->getInitialPorosity();
      IPVariable* ipvi = new  nonLocalPorosityDG3DIPVariable(*_nlpcoupledlaw,fvInit,inter);
      IPVariable* ipv1 = new  nonLocalPorosityDG3DIPVariable(*_nlpcoupledlaw,fvInit,inter);
      IPVariable* ipv2 = new  nonLocalPorosityDG3DIPVariable(*_nlpcoupledlaw,fvInit,inter);
      if(ips != NULL) delete ips;
      ips = new IP3State(state_,ipvi,ipv1,ipv2);
      _nlpcoupledlaw->createIPState((static_cast <nonLocalPorosityDG3DIPVariable*> (ipvi))->getIPNonLocalPorosity(),
                             (static_cast <nonLocalPorosityDG3DIPVariable*> (ipv1))->getIPNonLocalPorosity(),
                 			 (static_cast <nonLocalPorosityDG3DIPVariable*> (ipv2))->getIPNonLocalPorosity());
    
    }
    
    
    void NonLocalPorousCoupledDG3DMaterialLaw::createIPVariable(IPVariable* &ipv,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
    {
      if(ipv !=NULL) delete ipv;
      bool inter=true;
      const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
      if(iele == NULL){
        inter=false;
      }
      double fvInit = _nlpcoupledlaw->getInitialPorosity();
      ipv = new  nonLocalPorosityDG3DIPVariable(*_nlpcoupledlaw,fvInit,inter);
      IPNonLocalPorosity * ipvnl = static_cast <nonLocalPorosityDG3DIPVariable*>(ipv)->getIPNonLocalPorosity();
      _nlpcoupledlaw->createIPVariable(ipvnl,ele,nbFF_,GP,gpt);
    }
    
    
    
    void NonLocalPorousCoupledDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac)
    {
      /* get ipvariable */
      nonLocalPorosityDG3DIPVariable* ipvcur; //= static_cast < nonLocalDamageDG3DIPVariable *> (ipv);
      const nonLocalPorosityDG3DIPVariable* ipvprev; //= static_cast <const  nonLocalDamageDG3DIPVariable *> (ipvp);
    
      FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
      if(ipvtmp !=NULL)
      {
    
        ipvcur = static_cast<nonLocalPorosityDG3DIPVariable*>(ipvtmp->getIPvBulk());
        const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
        ipvprev = static_cast<const nonLocalPorosityDG3DIPVariable*>(ipvtmp2->getIPvBulk());
       }
      else
      {
        ipvcur = static_cast<nonLocalPorosityDG3DIPVariable*>(ipv);
        ipvprev = static_cast<const nonLocalPorosityDG3DIPVariable*>(ipvp);
      }
    
      /* strain xyz */
      const STensor3& F1 = ipvcur->getConstRefToDeformationGradient();
      const STensor3& F0 = ipvprev->getConstRefToDeformationGradient();
    
      /* data for J2 law */
      IPNonLocalPorosity* q1 = ipvcur->getIPNonLocalPorosity();
      const IPNonLocalPorosity* q0 = ipvprev->getIPNonLocalPorosity();
    
    	// set additional kinematic variables
    	q1->setRefToCurrentOutwardNormal(ipvcur->getConstRefToCurrentOutwardNormal());
    	q1->setRefToReferenceOutwardNormal(ipvcur->getConstRefToReferenceOutwardNormal());
    
      /* compute stress */
      _nlpcoupledlaw->constitutive(F0,F1,ipvcur->getRefToFirstPiolaKirchhoffStress(),q0,q1,ipvcur->getRefToTangentModuli(),
                            ipvcur->getRefToDLocalVariableDStrain()[0],ipvcur->getRefToDStressDNonLocalVariable()[0],
                            ipvcur->getRefToDLocalVariableDNonLocalVariable()(0,0),stiff);
    
      ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
    
    };
    
    
    //
    LocalDamageGursonDG3DMaterialLaw::LocalDamageGursonDG3DMaterialLaw(const int num,const double rho,
                        const double E,const double nu,
                        const double q1, const double q2, const double q3, const double fVinitial,
                        const J2IsotropicHardening &_j2IH,
                         const double tol,
                        const bool matrixByPert, const double pert) :
                                                                  dG3DMaterialLaw(num,rho,true)
    {
      _localGursonlaw = new mlawNonLocalDamageGurson(num,E,nu,rho,q1, q2, q3,
                                                  fVinitial,_j2IH,tol,matrixByPert,pert);
      fillElasticStiffness(E, nu, elasticStiffness);
    }
    
    LocalDamageGursonDG3DMaterialLaw::LocalDamageGursonDG3DMaterialLaw(const LocalDamageGursonDG3DMaterialLaw &source) :
                                                                 dG3DMaterialLaw(source)
    {
    	_localGursonlaw= NULL;
    	if (source._localGursonlaw != NULL){
    		_localGursonlaw = dynamic_cast<mlawNonLocalDamageGurson*>(source._localGursonlaw->clone());
    	}
    }
    
    LocalDamageGursonDG3DMaterialLaw::~LocalDamageGursonDG3DMaterialLaw(){
    	if (_localGursonlaw != NULL){
    		delete _localGursonlaw;
    		_localGursonlaw = NULL;
    	}
    }
    
    void LocalDamageGursonDG3DMaterialLaw::setOrderForLogExp(const int newOrder){_localGursonlaw->setOrderForLogExp(newOrder);};
    void LocalDamageGursonDG3DMaterialLaw::setPredictorCorrectorParameters(const int InitialGuessMethod, const int maxite, const double theta,  const int sys)
      {_localGursonlaw->setPredictorCorrectorParameters(InitialGuessMethod, maxite, theta, sys);};
    void LocalDamageGursonDG3DMaterialLaw::setViscosityLaw(const viscosityLaw& visco)
    {
      _localGursonlaw->setViscosityLaw(visco);
    };
    void LocalDamageGursonDG3DMaterialLaw::setNucleationLaw(const NucleationLaw& added_GdnLaw)
    {
      _localGursonlaw->setNucleationLaw(added_GdnLaw);
    };
    void LocalDamageGursonDG3DMaterialLaw::setScatterredInitialPorosity(const double fmin, const double fmax)
      {
        _localGursonlaw->setScatterredInitialPorosity(fmin,fmax);
    };
    void LocalDamageGursonDG3DMaterialLaw::setCoalescenceLaw(const CoalescenceLaw& added_coalesLaw)
      {_localGursonlaw->setCoalescenceLaw(added_coalesLaw);};
    
    void LocalDamageGursonDG3DMaterialLaw::setFailureTolerance(const double tol, const double elTol, const double localF){
      _localGursonlaw->setFailureTolerance(tol,elTol,localF);
    };
    
    void LocalDamageGursonDG3DMaterialLaw::setSubStepping(const bool fl, const int maxNumStep){
    	_localGursonlaw->setSubStepping(fl,maxNumStep);
    }
    void LocalDamageGursonDG3DMaterialLaw::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
    {
      // check interface element
      bool inter=true;
      const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
      if(iele==NULL) inter=false;
      double fvInit = _localGursonlaw->getInitialPorosity();
      IPVariable* ipvi = new  nonLocalPorosityDG3DIPVariable(*_localGursonlaw,fvInit,inter,0);
      IPVariable* ipv1 = new  nonLocalPorosityDG3DIPVariable(*_localGursonlaw,fvInit,inter,0);
      IPVariable* ipv2 = new  nonLocalPorosityDG3DIPVariable(*_localGursonlaw,fvInit,inter,0);
      if(ips != NULL) delete ips;
      ips = new IP3State(state_,ipvi,ipv1,ipv2);
      _localGursonlaw->createIPState((static_cast <nonLocalPorosityDG3DIPVariable*> (ipvi))->getIPNonLocalPorosity(),
                             (static_cast <nonLocalPorosityDG3DIPVariable*> (ipv1))->getIPNonLocalPorosity(),
                 			 (static_cast <nonLocalPorosityDG3DIPVariable*> (ipv2))->getIPNonLocalPorosity());
    
    }
    
    void LocalDamageGursonDG3DMaterialLaw::createIPVariable(IPVariable* &ipv,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
    {
      if(ipv !=NULL) delete ipv;
      bool inter=true;
      const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
      if(iele == NULL){
        inter=false;
      }
      double fvInit = _localGursonlaw->getInitialPorosity();
      ipv = new  nonLocalPorosityDG3DIPVariable(*_localGursonlaw,fvInit,inter,0);
      IPNonLocalPorosity * ipvnl = static_cast <nonLocalPorosityDG3DIPVariable*>(ipv)->getIPNonLocalPorosity();
      _localGursonlaw->createIPVariable(ipvnl,ele,nbFF_,GP,gpt);
    
    }
    
    
    
    void LocalDamageGursonDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac)
    {
      /* get ipvariable */
      nonLocalPorosityDG3DIPVariable* ipvcur = NULL; 
      const nonLocalPorosityDG3DIPVariable* ipvprev = NULL;
    
      FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
      if(ipvtmp !=NULL)
      {
        ipvcur = static_cast<nonLocalPorosityDG3DIPVariable*>(ipvtmp->getIPvBulk());
        const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
        ipvprev = static_cast<const nonLocalPorosityDG3DIPVariable*>(ipvtmp2->getIPvBulk());
       }
      else
      {
        ipvcur = static_cast<nonLocalPorosityDG3DIPVariable*>(ipv);
        ipvprev = static_cast<const nonLocalPorosityDG3DIPVariable*>(ipvp);
      }
    
      /* strain xyz */
      const STensor3& F1 = ipvcur->getConstRefToDeformationGradient();
      const STensor3& F0 = ipvprev->getConstRefToDeformationGradient();
    
      /* data for J2 law */
      IPNonLocalPorosity* q1 = ipvcur->getIPNonLocalPorosity();
      const IPNonLocalPorosity* q0 = ipvprev->getIPNonLocalPorosity();
      
    	// set additional kinematic variables
    	q1->setRefToCurrentOutwardNormal(ipvcur->getConstRefToCurrentOutwardNormal());
    	q1->setRefToReferenceOutwardNormal(ipvcur->getConstRefToReferenceOutwardNormal());
    
      /* compute stress */
      _localGursonlaw->constitutive(F0,F1,ipvcur->getRefToFirstPiolaKirchhoffStress(),q0,q1,ipvcur->getRefToTangentModuli(),stiff);
      ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
    };