diff --git a/NonLinearSolver/materialLaw/failureCriterion.cpp b/NonLinearSolver/materialLaw/failureCriterion.cpp
index 6aa10876fa3bf37ad505f0ed5698984b6e7b8be7..d77d7e27c9a40db6ddafd84a1d5c3735528b3f54 100644
--- a/NonLinearSolver/materialLaw/failureCriterion.cpp
+++ b/NonLinearSolver/materialLaw/failureCriterion.cpp
@@ -10,13 +10,138 @@
 
 #include "failureCriterion.h"
 #include "ipField.h"
+#include "ipHyperelastic.h"
 
 TrivialFailureCriterion::TrivialFailureCriterion(const int num): FailureCriterionBase(num,true){};
 
 CriticalDamageCriterion::CriticalDamageCriterion(const int num, const double Dc): FailureCriterionBase(num,true),_Dc(Dc){};
 
-bool CriticalDamageCriterion::isFailed(const IPVariable* ipv) const{
-	double D = ipv->get(IPField::DAMAGE);
+double CriticalDamageCriterion::getFailureCriterion(const IPVariable* ipvprev, const IPVariable* ipv) const{
+  return ipv->get(IPField::DAMAGE);
+};
+
+bool CriticalDamageCriterion::isFailed(const IPVariable* ipvprev, const IPVariable* ipv) const{
+	double D = this->getFailureCriterion(ipvprev,ipv);
 	if (D > _Dc) return true;
 	else return false;
+};
+
+PowerFailureCriterion::PowerFailureCriterion(const int num, const RateFailureLaw& trac, const RateFailureLaw& comp, const double power): FailureCriterionBase(num),
+    _failurePower(power){
+  _compFailure = comp.clone();
+  _tracFailure = trac.clone();
+};
+
+PowerFailureCriterion::PowerFailureCriterion(const PowerFailureCriterion& src): FailureCriterionBase(src),_failurePower(src._failurePower){
+  _compFailure = NULL;
+  if (src._compFailure){
+    _compFailure= src._compFailure->clone();
+  }
+  _tracFailure = NULL;
+  if (src._tracFailure){
+    _tracFailure = src._tracFailure->clone();
+  }
+  
+};
+PowerFailureCriterion::~PowerFailureCriterion(){
+  if (_compFailure){
+    delete _compFailure;
+    _compFailure = NULL;
+  }
+  if (_tracFailure){
+    delete _tracFailure;
+    _tracFailure = NULL;
+  }
+}
+
+double PowerFailureCriterion::getFailureCriterion(const IPVariable* ipvprev, const IPVariable* ipv) const{
+  const IPHyperViscoElastoPlastic* viscoipv = dynamic_cast<const IPHyperViscoElastoPlastic*>(ipv);
+  if (viscoipv == NULL){
+    Msg::Error("PowerFailureCriterion::isFailed is implemented for IPHyperViscoElastoPlastic only");
+    return 1e10;
+  }
+  else{
+    // compute failure r and DrDF
+    const STensor3& K = viscoipv->_kirchhoff;
+    double p = K.trace()/3.;
+    STensor3 devK = K.dev();
+    double Keq = sqrt(1.5*devK.dotprod());
+    
+    double Xc, dXc, Xt, dXt;
+    const IPHyperViscoElastoPlastic* viscoipvprev = static_cast<const IPHyperViscoElastoPlastic*>(ipvprev);
+    _compFailure->get(viscoipvprev->_DgammaDt,Xc,dXc);
+    _tracFailure->get(viscoipvprev->_DgammaDt,Xt,dXt);
+
+    double m = Xt/Xc;
+    double mExal = pow(m,_failurePower);
+    double KeqXc = Keq/Xc;
+    double pXc = p/Xc;
+    double KeqExal = 0.;
+
+    if (KeqXc > 1e-20){
+      KeqExal = pow(KeqXc,_failurePower);
+    }
+
+    double b2 = (m+1.)/(mExal+m);
+    double b1 = (1.- mExal)/(mExal+m);
+    double b0 = 1.;
+
+    double Phif = b2*KeqExal+  b1* 3.*pXc - b0;
+    return Phif;
+  }
+};
+
+bool PowerFailureCriterion::isFailed(const IPVariable* ipvprev, const IPVariable* ipv) const{
+  const IPHyperViscoElastoPlastic* viscoipv = dynamic_cast<const IPHyperViscoElastoPlastic*>(ipv);
+  if (viscoipv == NULL){
+    Msg::Error("PowerFailureCriterion::isFailed is implemented for IPHyperViscoElastoPlastic only");
+    return false;
+  }
+  else{
+    double Phif = this->getFailureCriterion(ipvprev,ipv);
+    double r = viscoipv->getFailureOnset();  
+    double Ff =  Phif  - r;
+    if (Ff >=0) return true;
+    else return false;
+  }
+};
+
+CritialPlasticDeformationCriteration::CritialPlasticDeformationCriteration(const int num, const double a, const double b, const double c): FailureCriterionBase(num),
+  _a(a),_b(b),_c(c){};
+  
+CritialPlasticDeformationCriteration::CritialPlasticDeformationCriteration(const CritialPlasticDeformationCriteration& src): FailureCriterionBase(src),
+  _a(src._a),_b(src._b),_c(src._c){};
+
+double CritialPlasticDeformationCriteration::getFailureCriterion(const IPVariable* ipvprev, const IPVariable* ipv) const{
+  const IPHyperViscoElastoPlastic* viscoipv = dynamic_cast<const IPHyperViscoElastoPlastic*>(ipv);
+  if (viscoipv == NULL){
+    Msg::Error("PowerFailureCriterion::isFailed is implemented for IPHyperViscoElastoPlastic only");
+    return 1e10;
+  }
+  else{
+    // compute failure r and DrDF
+    const STensor3& K = viscoipv->_kirchhoff;
+    double pres = K.trace()/3.;
+    STensor3 devK = K.dev();
+    double Keq = sqrt(1.5*devK.dotprod());
+    double p = viscoipv->getConstRefToEqPlasticStrain();
+    
+    double T = pres/Keq;
+    if (T != T) T = 0.;
+    return (p - _a*(exp(-_b*T)) - _c);
+  }
+};
+bool CritialPlasticDeformationCriteration::isFailed(const IPVariable* ipvprev, const IPVariable* ipv) const{
+  const IPHyperViscoElastoPlastic* viscoipv = dynamic_cast<const IPHyperViscoElastoPlastic*>(ipv);
+  if (viscoipv == NULL){
+    Msg::Error("PowerFailureCriterion::isFailed is implemented for IPHyperViscoElastoPlastic only");
+    return false;
+  }
+  else{
+    double Phif = this->getFailureCriterion(ipvprev,ipv);
+    double r = viscoipv->getFailureOnset();  
+    double Ff =  Phif  - r;
+    if (Ff >=0) return true;
+    else return false;
+  }
 };
\ No newline at end of file
diff --git a/NonLinearSolver/materialLaw/failureCriterion.h b/NonLinearSolver/materialLaw/failureCriterion.h
index ea772e7d0f1c0431d18ea98324bf164c10f1f00b..4687df58678f9293df0e05f8f6a74282bea6256f 100644
--- a/NonLinearSolver/materialLaw/failureCriterion.h
+++ b/NonLinearSolver/materialLaw/failureCriterion.h
@@ -27,7 +27,8 @@ class FailureCriterionBase{
     FailureCriterionBase(const FailureCriterionBase& src): _num(src._num),_init(src._init){}
     virtual ~FailureCriterionBase(){};
     // failure check based on stress state tensor and internal variable
-    virtual bool isFailed(const IPVariable* ipv) const = 0;
+    virtual bool isFailed(const IPVariable* ipvprev, const IPVariable* ipv) const = 0;
+    virtual double getFailureCriterion(const IPVariable* ipvprev, const IPVariable* ipv) const = 0;
     virtual FailureCriterionBase* clone() const = 0;
   #endif // SWIG
 };
@@ -38,7 +39,8 @@ class TrivialFailureCriterion : public FailureCriterionBase{
     #ifndef SWIG
     TrivialFailureCriterion(const TrivialFailureCriterion& src): FailureCriterionBase(src){}
     virtual ~TrivialFailureCriterion(){}
-    virtual bool isFailed(const IPVariable* ipv) const {return false;};
+    virtual bool isFailed(const IPVariable* ipvprev, const IPVariable* ipv) const {return false;};
+    virtual double getFailureCriterion(const IPVariable* ipvprev, const IPVariable* ipv) const {return 1.;};
     virtual FailureCriterionBase* clone() const {return new TrivialFailureCriterion(*this);};
     #endif // SWIG
 };
@@ -52,10 +54,47 @@ class CriticalDamageCriterion : public FailureCriterionBase{
 		#ifndef SWIG
 		CriticalDamageCriterion(const CriticalDamageCriterion& src):FailureCriterionBase(src),_Dc(src._Dc){}
 		virtual ~CriticalDamageCriterion(){}
-		
-		virtual bool isFailed(const IPVariable* ipv) const;
+    virtual double getFailureCriterion(const IPVariable* ipvprev, const IPVariable* ipv) const;
+		virtual bool isFailed(const IPVariable* ipvprev, const IPVariable* ipv) const;
     virtual FailureCriterionBase* clone() const {return new CriticalDamageCriterion(*this);};
 		#endif // SWIG
 };
 
+
+class PowerFailureCriterion : public FailureCriterionBase{
+  protected:
+    #ifndef SWIG
+    RateFailureLaw* _compFailure;
+    RateFailureLaw* _tracFailure;
+    double _failurePower;
+    #endif //SWIG
+
+  public:
+    PowerFailureCriterion(const int num, const RateFailureLaw& trac, const RateFailureLaw& comp, const double power);
+    #ifndef SWIG
+    PowerFailureCriterion(const PowerFailureCriterion& src);
+    virtual ~PowerFailureCriterion();
+    virtual double getFailureCriterion(const IPVariable* ipvprev, const IPVariable* ipv) const;
+    virtual bool isFailed(const IPVariable* ipvprev, const IPVariable* ipv) const;
+    virtual FailureCriterionBase* clone() const {return new PowerFailureCriterion(*this);};
+    #endif //SWIG  
+};
+
+class CritialPlasticDeformationCriteration : public FailureCriterionBase{
+  protected:
+    #ifndef SWIG
+    double _a, _b, _c; // p_cr = a*(exp(-b*X))+c (X-stress-triaxiality = p/sigVM)
+    #endif //SWIG
+ 
+  public:
+    CritialPlasticDeformationCriteration(const int num, const double a, const double b, const double c);
+    #ifndef SWIG
+    CritialPlasticDeformationCriteration(const CritialPlasticDeformationCriteration& src);
+    virtual ~CritialPlasticDeformationCriteration(){};
+    virtual double getFailureCriterion(const IPVariable* ipvprev, const IPVariable* ipv) const;
+    virtual bool isFailed(const IPVariable* ipvprev, const IPVariable* ipv) const;
+    virtual FailureCriterionBase* clone() const {return new CritialPlasticDeformationCriteration(*this);};
+    #endif //SWIG
+  
+};
 #endif // FAILURECRITERION_H_
diff --git a/NonLinearSolver/materialLaw/mlawHyperelastic.cpp b/NonLinearSolver/materialLaw/mlawHyperelastic.cpp
index 88b3d59b48e117fe5bd8d582b23c0d40eb68d3ac..62ee704ab2378632ef263068fcc6d4e8574a8fe1 100644
--- a/NonLinearSolver/materialLaw/mlawHyperelastic.cpp
+++ b/NonLinearSolver/materialLaw/mlawHyperelastic.cpp
@@ -1771,21 +1771,24 @@ void mlawPowerYieldHyper::createIPState(IPStateBase* &ips,const bool* state_,con
 
 mlawPowerYieldHyperWithFailure::mlawPowerYieldHyperWithFailure(const int num,const double E,const double nu, const double rho,
                     const double tol, const bool matrixbyPerturbation, const double pert):
-                    mlawPowerYieldHyper(num,E,nu,rho,tol,matrixbyPerturbation,pert),_compFailure(NULL),_tracFailure(NULL){};
+                    mlawPowerYieldHyper(num,E,nu,rho,tol,matrixbyPerturbation,pert){
+  _failureCriterion = new TrivialFailureCriterion(num);
+};
 
 
 
 mlawPowerYieldHyperWithFailure::mlawPowerYieldHyperWithFailure(const mlawPowerYieldHyperWithFailure& src):
     mlawPowerYieldHyper(src){
-    _failurePower = src._failurePower;
-    _tracFailure = NULL;
-    _compFailure = NULL;
-    if (src._tracFailure != NULL) _tracFailure = src._tracFailure->clone();
-    if (src._compFailure != NULL) _compFailure = src._compFailure->clone();
+  _failureCriterion = NULL;
+  if (src._failureCriterion){
+    _failureCriterion = src._failureCriterion->clone();
+  }
 };
 mlawPowerYieldHyperWithFailure::~mlawPowerYieldHyperWithFailure(){
-  if (_compFailure != NULL) delete _compFailure;
-  if (_tracFailure != NULL) delete _tracFailure;
+  if (_failureCriterion != NULL) {
+    delete _failureCriterion;
+    _failureCriterion = NULL;
+  }
 };
 
 
@@ -1794,67 +1797,26 @@ void mlawPowerYieldHyperWithFailure::predictorCorrector(const STensor3& F, const
   mlawPowerYieldHyper::predictorCorrector(F,q0,q1,P,stiff,Tangent,dFedF,dFpdF);
 
   // compute failure r and DrDF
-  const STensor3& K = q1->_kirchhoff;
-  double p = K.trace()/3.;
-  STensor3 devK = K.dev();
-  double Keq = sqrt(1.5*devK.dotprod());
-
   double& r = q1->getRefToFailureOnset();
   double& gF = q1->getRefToFailurePlasticity();
-
   r = q0->getFailureOnset();
   gF = q0->getFailurePlasticity();
-
-  STensor3& dgFdF = q1->getRefToDFailurePlasticityDF();
-  dgFdF *= 0.;
-
-  double Xc, dXc, Xt, dXt;
-  _compFailure->get(q1->_DgammaDt,Xc,dXc);
-  _tracFailure->get(q1->_DgammaDt,Xt,dXt);
-
-  double m = Xt/Xc;
-  double mExal = pow(m,_failurePower);
-  double KeqXc = Keq/Xc;
-  double pXc = p/Xc;
-  double KeqExal = 0.;
-
-  if (KeqXc > 1e-20){
-    KeqExal = pow(KeqXc,_failurePower);
-  }
-
-  double b2 = (m+1.)/(mExal+m);
-  double b1 = (1.- mExal)/(mExal+m);
-  double b0 = 1.;
-
-  double Phif = b2*KeqExal+  b1* 3.*pXc - b0;
-
-
-  double Ff =  Phif  - r;
-
-  //Msg::Error("Xc =%e, Xt=%e, m = %e, Ff = %e Phif = %e, r= %e",Xc,Xt,m,Ff,Phif,r);
-
-  if ((Ff >0)){
-    //Msg::Error("Fracture onset Xc =%e, Xt=%e, m = %e, Ff = %e Phif = %e, r= %e",Xc,Xt,m,Ff,Phif,r);
-    r = Phif;
+  
+  if (_failureCriterion->isFailed(q0,q1)){
+    r = _failureCriterion->getFailureCriterion(q0,q1);
     gF += q1->_epspbarre - q0->_epspbarre;
     if (stiff){
-      dgFdF = q1->_DgammaDF;
+      q1->getRefToDFailurePlasticityDF() = q1->_DgammaDF;
     }
-  };
-
-};
-
-
-void mlawPowerYieldHyperWithFailure::setCompressiveFailureModel(const RateFailureLaw& compStrength){
-  if (_compFailure != NULL) delete _compFailure;
-  _compFailure = compStrength.clone();
-};
-
-void mlawPowerYieldHyperWithFailure::setTractionFailureModel(const RateFailureLaw& tracStrength){
-  if (_tracFailure != NULL) delete _tracFailure;
-  _tracFailure = tracStrength.clone();
+  }
+  else{
+    if (stiff){
+      q1->getRefToDFailurePlasticityDF() *= 0.;
+    }
+  }
 };
 
-void mlawPowerYieldHyperWithFailure::setFailurePower(const double power){
-  _failurePower = power;
-};
+void mlawPowerYieldHyperWithFailure::setFailureCriterion(const FailureCriterionBase& fCr){
+  if (_failureCriterion) delete _failureCriterion;
+  _failureCriterion = fCr.clone();
+};
\ No newline at end of file
diff --git a/NonLinearSolver/materialLaw/mlawHyperelastic.h b/NonLinearSolver/materialLaw/mlawHyperelastic.h
index 2ed3f00e78548906e09c6c5df1f124179d7c3040..b96dd9e2a16e66dec224fe5eaa84b81c863fc07d 100644
--- a/NonLinearSolver/materialLaw/mlawHyperelastic.h
+++ b/NonLinearSolver/materialLaw/mlawHyperelastic.h
@@ -16,6 +16,7 @@
 #include "ipHyperelastic.h"
 #include "viscosityLaw.h"
 #include "materialStrengthLaw.h"
+#include "failureCriterion.h"
 
 class mlawHyperViscoElastic : public materialLaw{
   #ifndef SWIG
@@ -211,25 +212,20 @@ class mlawPowerYieldHyper : public mlawHyperViscoElastic{
 
 class mlawPowerYieldHyperWithFailure : public mlawPowerYieldHyper{
   protected:
-    RateFailureLaw* _compFailure;
-    RateFailureLaw* _tracFailure;
-    double _failurePower;
-
-  #ifndef SWIG
+    #ifndef SWIG
+    FailureCriterionBase* _failureCriterion;
   protected:
     virtual void predictorCorrector(const STensor3& F, const IPHyperViscoElastoPlastic *q0_, IPHyperViscoElastoPlastic *q1,
                             STensor3&P, const bool stiff, STensor43& Tangent, STensor43& dFedF, STensor43& dFpdF) const;
-  #endif // SWIG
+    #endif // SWIG
   public:
     mlawPowerYieldHyperWithFailure(const int num,const double E,const double nu, const double rho,
                         const double tol = 1.e-6,
                         const bool matrixbyPerturbation = false, const double pert = 1e-8);
 
 
-    void setCompressiveFailureModel(const RateFailureLaw& compStrength);
-    void setTractionFailureModel(const RateFailureLaw& tracStrength);
-    void setFailurePower(const double power);
-
+    void setFailureCriterion(const FailureCriterionBase& fCr);
+    
   #ifndef SWIG
     mlawPowerYieldHyperWithFailure(const mlawPowerYieldHyperWithFailure& src);
     virtual ~mlawPowerYieldHyperWithFailure();
diff --git a/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp b/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
index b8535503a621f29fc886a7e9f07cb492682bcd36..811ad2ed6efc132357fdf22289568709aa529f01 100644
--- a/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
+++ b/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
@@ -1313,11 +1313,12 @@ void nonLinearMechSolver::checkElementErosion(IPField* ipf, const IPStateBase::w
         bool willBeEroded = false;
         
         for (int ip=0; ip < (*ips).size(); ip++){
+          const IPVariable* ipvprev = (*ips)[ip]->getState(IPStateBase::previous);
           const IPVariable* ipv = (*ips)[ip]->getState(ws);
           // check element
           if (_globalCheckFlag){
             // global check
-            if (_erosionGlobalCriterion->isFailed(ipv)){
+            if (_erosionGlobalCriterion->isFailed(ipvprev,ipv)){
               willBeEroded = true;
               break;
             }
@@ -1345,8 +1346,9 @@ void nonLinearMechSolver::checkElementErosion(IPField* ipf, const IPStateBase::w
         int numberNoFailedIP= 0;
         for (int ip=0; ip < (*ips).size(); ip++){
           const IPVariable* ipv = (*ips)[ip]->getState(ws);
+          const IPVariable* ipvprev = (*ips)[ip]->getState(IPStateBase::previous);
           if (_globalCheckFlag){
-            if (!_erosionGlobalCriterion->isFailed(ipv)){
+            if (!_erosionGlobalCriterion->isFailed(ipvprev,ipv)){
               numberNoFailedIP++;
               break;
             }
diff --git a/NonLinearSolver/nlsolver/scalarFunction.h b/NonLinearSolver/nlsolver/scalarFunction.h
index 9ad0c811e1395a46d14cf286ac6218fdeb7acc7e..5787c9ef27b7e03ea27276db776635df17b16f7e 100644
--- a/NonLinearSolver/nlsolver/scalarFunction.h
+++ b/NonLinearSolver/nlsolver/scalarFunction.h
@@ -19,6 +19,7 @@ class scalarFunction{
   scalarFunction(const scalarFunction& src){}
   virtual scalarFunction& operator=(const scalarFunction& src){return *this;}
   virtual ~scalarFunction(){};
+  virtual void setElement(MElement* el) const{}
   virtual double getVal(const double x) const = 0; // function value at x
   virtual double getDiff(const double x) const = 0; // function derivative at x
   virtual scalarFunction* clone() const = 0;
@@ -49,6 +50,42 @@ class constantScalarFunction : public scalarFunction{
     #endif // SWIG
 };
 
+class elementWiseScalarFunction : public constantScalarFunction{
+  protected:
+    std::map<int, double> _constantOnElements;
+    mutable MElement* _curELe;
+  public:
+    elementWiseScalarFunction(const double A) : constantScalarFunction(A),_curELe(NULL){}
+    void setValueOnElement(const int i, const double val) {_constantOnElements[i] = val;};
+    #ifndef SWIG
+    elementWiseScalarFunction(const elementWiseScalarFunction& src) : constantScalarFunction(src),_constantOnElements(src._constantOnElements),_curELe(NULL){}
+    virtual elementWiseScalarFunction& operator = (const scalarFunction& src){
+      constantScalarFunction::operator=(src);
+      const elementWiseScalarFunction* psrc = dynamic_cast<const elementWiseScalarFunction*>(&src);
+      if (psrc!= NULL){
+        _constantOnElements = psrc->_constantOnElements;
+        _curELe = NULL;
+      }
+      return *this;
+    }
+
+    virtual void setElement(MElement* el) const {_curELe = el;} // must be acessed during constobject
+    
+    virtual ~elementWiseScalarFunction(){};
+    virtual double getVal(const double x) const {
+      if (!_curELe) return _f;
+      std::map<int,double>::const_iterator it = _constantOnElements.find(_curELe->getNum());
+      if (it != _constantOnElements.end()) return it->second;
+      return _f;
+    };
+    virtual double getDiff(const double x) const {
+      Msg::Error("elementWiseScalarFunction::getDiff should not be used");
+      return 0.;
+    };
+    virtual scalarFunction* clone() const {return new elementWiseScalarFunction(*this);};
+    #endif //SWIG
+};
+
 class linearScalarFunction : public scalarFunction{
   protected:
     double _f0, _x0, _a; // f(x)= f0+a*(x-x0)
diff --git a/cm3apps/install.txt b/cm3apps/install.txt
index 0796462901587db6f2d32af33cbaf20597256497..2f6b95acdbb0d018fa821249f3f3cefb49c38ef1 100644
--- a/cm3apps/install.txt
+++ b/cm3apps/install.txt
@@ -483,6 +483,7 @@ C) Now you can install your project(s) (dgshell and/or dG3D and/or msch), NonLin
              CMAKE_BUILD_TYPE release (or debug)
              ENABLE_BUILD_DYNAMICS = ON (dynamic executable using the lib to not compile everything twice)
              ENABLE_BUILD_SHARED = ON (creates the Gmsh library you need)
+             DENABLE_INTERNAL_DEVELOPER_API=ON (to get the full gmsh code)
              ENABLE_MPI = ON (for // implementation)
              ENABLE_TAUCS=OFF (or install it but it's hard)
              ENABLE_SLEPC=ON
diff --git a/dG3D/benchmarks/interpolationPBC_FE1/idealHole.py b/dG3D/benchmarks/interpolationPBC_FE1/idealHole.py
index d0ea28e7e2d399840b9bbe915b3cf6e7491c2771..b2f819df37f0130f14e802bf663c2cda5739d1c0 100644
--- a/dG3D/benchmarks/interpolationPBC_FE1/idealHole.py
+++ b/dG3D/benchmarks/interpolationPBC_FE1/idealHole.py
@@ -67,7 +67,7 @@ microBC.setBCPhysical(1,3,5,2,4,6)
 method =3	# Periodic mesh = 0, Langrange interpolation = 1, Cubic spline interpolation =2,  FE linear= 3, FE Quad = 4
 degree = 7	# Order used for polynomial interpolation 
 addvertex = 1 # Polynomial interpolation by mesh vertex = 0, Polynomial interpolation by virtual vertex 
-
+microBC.forceInterpolationDegree(True)
 microBC.setPeriodicBCOptions(method, degree,bool(addvertex)) 
  # Deformation gradient
 microBC.setDeformationGradient(1.0,0.02,0.02,0.0,0.97,0,0,0,1.02)
@@ -108,10 +108,10 @@ mysolver.internalPointBuildView("Equivalent plastic strain",IPField.PLASTICSTRAI
 mysolver.solve()
 
 check = TestCheck()
-check.equal(-1.971206e+02,mysolver.getHomogenizedStress(0,0),1.e-4)
-check.equal(6.750333e+02,mysolver.getHomogenizedStress(0,1),1.e-4)
-check.equal(-3.06250E+03,mysolver.getHomogenizedStress(1,1),1.e-4)
-check.equal(1.383659e+05,mysolver.getHomogenizedTangent(0,0,0,0),1.e-4)
-check.equal(1.383432e+05,mysolver.getHomogenizedTangent(1,1,1,1),1.e-4)
-check.equal(3.381458e+04,mysolver.getHomogenizedTangent(0,1,0,1),1.e-4)
+check.equal(-1.968966e+02,mysolver.getHomogenizedStress(0,0),1.e-4)
+check.equal(6.710469e+02,mysolver.getHomogenizedStress(0,1),1.e-4)
+check.equal(-3.059271e+03,mysolver.getHomogenizedStress(1,1),1.e-4)
+check.equal(1.381866e+05,mysolver.getHomogenizedTangent(0,0,0,0),1.e-4)
+check.equal(1.381850e+05,mysolver.getHomogenizedTangent(1,1,1,1),1.e-4)
+check.equal(3.359562e+04,mysolver.getHomogenizedTangent(0,1,0,1),1.e-4)
 
diff --git a/dG3D/benchmarks/powerYieldViscoElastoPlasticFullFailure/model.py b/dG3D/benchmarks/powerYieldViscoElastoPlasticFullFailure/model.py
index 8a8e5a3ca42da684a2fb4e21a5d37a6b71f6896d..c4cff304cad4d94f4c0ed51caf5860cd1a05371e 100644
--- a/dG3D/benchmarks/powerYieldViscoElastoPlasticFullFailure/model.py
+++ b/dG3D/benchmarks/powerYieldViscoElastoPlasticFullFailure/model.py
@@ -50,9 +50,9 @@ lawMatrix.setViscoElasticData(4,48e6,0.073)
 Xc = ConstantFailureLaw(1,262.5e6)
 Xt = ConstantFailureLaw(2,100e6)
 
-lawMatrix.setCompressiveFailureModel(Xc)
-lawMatrix.setTractionFailureModel(Xt)
-lawMatrix.setFailurePower(2)
+failureCr = PowerFailureCriterion(1,Xt,Xc,2.)
+lawMatrix.setFailureCriterion(failureCr)
+
 eta = constantViscosityLaw(1,3e10)
 lawMatrix.setViscosityEffect(eta,0.21)
 
diff --git a/dG3D/src/NonLocalDamageHyperelasticDG3DMaterialLaw.cpp b/dG3D/src/NonLocalDamageHyperelasticDG3DMaterialLaw.cpp
index de57e1bb3e78b2c4cf3b97259c400f611d663a6b..a767eba5e4824623e897847589097cec6b4d4652 100644
--- a/dG3D/src/NonLocalDamageHyperelasticDG3DMaterialLaw.cpp
+++ b/dG3D/src/NonLocalDamageHyperelasticDG3DMaterialLaw.cpp
@@ -286,16 +286,9 @@ void LocalDamagePowerYieldHyperDG3DMaterialLawWithFailure::setViscoElasticData(c
   _localPowerYieldHyperlaw->setViscoElasticData(filename);
 }
 
-void LocalDamagePowerYieldHyperDG3DMaterialLawWithFailure::setCompressiveFailureModel(const RateFailureLaw& compStrength){
-   _localPowerYieldHyperlaw->setCompressiveFailureModel(compStrength);
+void LocalDamagePowerYieldHyperDG3DMaterialLawWithFailure::setFailureCriterion(const FailureCriterionBase& fCr){
+  _localPowerYieldHyperlaw->setFailureCriterion(fCr);
 };
-void LocalDamagePowerYieldHyperDG3DMaterialLawWithFailure::setTractionFailureModel(const RateFailureLaw& tracStrength){
-  _localPowerYieldHyperlaw->setTractionFailureModel(tracStrength);
-};
-void LocalDamagePowerYieldHyperDG3DMaterialLawWithFailure::setFailurePower(const double power){
-   _localPowerYieldHyperlaw->setFailurePower(power);
-};
-
 LocalDamagePowerYieldHyperDG3DMaterialLawWithFailure::~LocalDamagePowerYieldHyperDG3DMaterialLawWithFailure(){
   if (_localPowerYieldHyperlaw) delete _localPowerYieldHyperlaw;
 	_localPowerYieldHyperlaw = NULL;
@@ -485,14 +478,8 @@ void NonLocalDamagePowerYieldHyperDG3DMaterialLawWithFailure::setViscoElasticDat
   _nldPowerYieldHyperlaw->setViscoElasticData(filename);
 }
 
-void NonLocalDamagePowerYieldHyperDG3DMaterialLawWithFailure::setCompressiveFailureModel(const RateFailureLaw& compStrength){
-   _nldPowerYieldHyperlaw->setCompressiveFailureModel(compStrength);
-};
-void NonLocalDamagePowerYieldHyperDG3DMaterialLawWithFailure::setTractionFailureModel(const RateFailureLaw& tracStrength){
-  _nldPowerYieldHyperlaw->setTractionFailureModel(tracStrength);
-};
-void NonLocalDamagePowerYieldHyperDG3DMaterialLawWithFailure::setFailurePower(const double power){
-   _nldPowerYieldHyperlaw->setFailurePower(power);
+void NonLocalDamagePowerYieldHyperDG3DMaterialLawWithFailure::setFailureCriterion(const FailureCriterionBase& fCr){
+  _nldPowerYieldHyperlaw->setFailureCriterion(fCr);
 };
 
 NonLocalDamagePowerYieldHyperDG3DMaterialLawWithFailure::~NonLocalDamagePowerYieldHyperDG3DMaterialLawWithFailure(){
diff --git a/dG3D/src/NonLocalDamageHyperelasticDG3DMaterialLaw.h b/dG3D/src/NonLocalDamageHyperelasticDG3DMaterialLaw.h
index 3b9286ef6e6b40ff9ab71e37e6411712f98aebf7..c33b96139f44e6d3d9ee94ebc18c494a25b71269 100644
--- a/dG3D/src/NonLocalDamageHyperelasticDG3DMaterialLaw.h
+++ b/dG3D/src/NonLocalDamageHyperelasticDG3DMaterialLaw.h
@@ -99,11 +99,9 @@ class LocalDamagePowerYieldHyperDG3DMaterialLawWithFailure: public dG3DMaterialL
   void setViscoElasticData_Bulk(const int i, const double Ki, const double ki);
   void setViscoElasticData_Shear(const int i, const double Gi, const double gi);
   void setViscoElasticData(const std::string filename);
-
-  void setCompressiveFailureModel(const RateFailureLaw& compStrength);
-  void setTractionFailureModel(const RateFailureLaw& tracStrength);
-  void setFailurePower(const double power);
-
+  
+  void setFailureCriterion(const FailureCriterionBase& fCr);
+  
 #ifndef SWIG
   virtual ~LocalDamagePowerYieldHyperDG3DMaterialLawWithFailure();
   LocalDamagePowerYieldHyperDG3DMaterialLawWithFailure(const LocalDamagePowerYieldHyperDG3DMaterialLawWithFailure &source);
@@ -161,9 +159,7 @@ class NonLocalDamagePowerYieldHyperDG3DMaterialLawWithFailure: public dG3DMateri
   void setViscoElasticData_Shear(const int i, const double Gi, const double gi);
   void setViscoElasticData(const std::string filename);
 
-  void setCompressiveFailureModel(const RateFailureLaw& compStrength);
-  void setTractionFailureModel(const RateFailureLaw& tracStrength);
-  void setFailurePower(const double power);
+  void setFailureCriterion(const FailureCriterionBase& fCr);
 
 #ifndef SWIG
   virtual ~NonLocalDamagePowerYieldHyperDG3DMaterialLawWithFailure();
diff --git a/dG3D/src/dG3D1DDomain.cpp b/dG3D/src/dG3D1DDomain.cpp
index 8b28d0fc123333d5b9608edab7e851dde3c4c9e6..eabb5ebe7a9187360fac37f2403d8b1204de3676 100644
--- a/dG3D/src/dG3D1DDomain.cpp
+++ b/dG3D/src/dG3D1DDomain.cpp
@@ -16,33 +16,43 @@
 #include "nonLinearMechSolver.h"
 #include "pathFollowingTerms.h"
 
-
-dG3D1DDomain::dG3D1DDomain(const int tag, const int phys, const int sp_, const int lnum, const int fdg) : dG3DDomain(tag,phys,sp_,lnum,fdg,1),_crossSection(1.){};
-
-double dG3D1DDomain::getCrossSection(const int el) const{
-  std::map<int,double>::const_iterator it = _otherCrossSection.find(el);
-  if (it != _otherCrossSection.end()){
-    return it->second;
+dG3D1DDomain::dG3D1DDomain(const dG3D1DDomain& src) : dG3DDomain(src){
+  if (src._crossSectionDistribution){
+    _crossSectionDistribution = src._crossSectionDistribution->clone();
+  }
+  else{
+    _crossSectionDistribution = new constantScalarFunction(1.);
   }
-  else
-    return _crossSection;
 };
+dG3D1DDomain::~dG3D1DDomain(){
+  if (_crossSectionDistribution) {
+    delete _crossSectionDistribution;
+    _crossSectionDistribution = NULL;
+  }
+}
 
-void dG3D1DDomain::setCrossSection(const double A){
-  _crossSection = A;
+void dG3D1DDomain::setCrossSectionDistribution(const scalarFunction* fct){
+  if (_crossSectionDistribution) delete _crossSectionDistribution;
+  _crossSectionDistribution = fct->clone();
 };
-void dG3D1DDomain::setCrossSectionOnElement(const int el, const double A){
-  _otherCrossSection[el] = A;
+
+
+dG3D1DDomain::dG3D1DDomain(const int tag, const int phys, const int sp_, const int lnum, const int fdg) : dG3DDomain(tag,phys,sp_,lnum,fdg,1){
+  _crossSectionDistribution = new constantScalarFunction(1.);
+  #if _DEBUG
+  Msg::Warning("1D must be defined in x-coordinate only");
+  #endif //
 };
 
+
 void dG3D1DDomain::createTerms(unknownField *uf,IPField*ip)
 {
   FunctionSpace<double>* dgspace = static_cast<FunctionSpace<double>*>(_space);
-  ltermBulk = new dG3D1DForceBulk(this,ip);
+  ltermBulk = new dG3D1DForceBulk(this,ip,_crossSectionDistribution);
   if(_bmbp)
     btermBulk = new BilinearTermPerturbation<double>(ltermBulk,*dgspace,*dgspace,uf,ip,this,_eps);
   else
-   btermBulk = new dG3D1DStiffnessBulk(this,ip);
+   btermBulk = new dG3D1DStiffnessBulk(this,ip,_crossSectionDistribution);
 
   
   ltermBound = new nonLinearTermVoid();
@@ -91,33 +101,42 @@ void dG3D1DDomain::createTerms(unknownField *uf,IPField*ip)
 
 
 nonLocalDamageDG3D1DDDomain::nonLocalDamageDG3D1DDDomain(const int tag, const int phys, const int sp_, const int lnum, const int fdg, const double nonLocalEqRatio,const int numNonlocalVar):
-              nonLocalDamageDG3DDomain(tag,phys,sp_,lnum,fdg,nonLocalEqRatio,1,numNonlocalVar),_crossSection(1.){};
+              nonLocalDamageDG3DDomain(tag,phys,sp_,lnum,fdg,nonLocalEqRatio,1,numNonlocalVar){
+  _crossSectionDistribution = new constantScalarFunction(1.);        
+  #if _DEBUG
+  Msg::Warning("1D must be defined in x-coordinate only");
+  #endif //    
+};
               
-double nonLocalDamageDG3D1DDDomain::getCrossSection(const int el) const{
-  std::map<int,double>::const_iterator it = _otherCrossSection.find(el);
-  if (it != _otherCrossSection.end()){
-    return it->second;
+nonLocalDamageDG3D1DDDomain::nonLocalDamageDG3D1DDDomain(const nonLocalDamageDG3D1DDDomain& src): nonLocalDamageDG3DDomain(src){
+  if (src._crossSectionDistribution){
+    _crossSectionDistribution = src._crossSectionDistribution->clone();
+  }
+  else{
+    _crossSectionDistribution = new constantScalarFunction(1.);
   }
-  else
-    return _crossSection;
 };
 
-void nonLocalDamageDG3D1DDDomain::setCrossSection(const double A){
-  _crossSection = A;
+nonLocalDamageDG3D1DDDomain::~nonLocalDamageDG3D1DDDomain(){
+  if (_crossSectionDistribution) {
+    delete _crossSectionDistribution; _crossSectionDistribution = NULL;
+  }
 };
-void nonLocalDamageDG3D1DDDomain::setCrossSectionOnElement(const int el, const double A){
-  _otherCrossSection[el] = A;
+
+void nonLocalDamageDG3D1DDDomain::setCrossSectionDistribution(const scalarFunction* fct){
+  if (_crossSectionDistribution) delete _crossSectionDistribution;
+  _crossSectionDistribution = fct->clone();
 };
 
 
 void nonLocalDamageDG3D1DDDomain::createTerms(unknownField *uf,IPField*ip)
 {
   FunctionSpace<double>* dgspace = static_cast<FunctionSpace<double>*>(_space);
-  ltermBulk = new nonLocalDG3D1DForceBulk(this,ip,_nonLocalEqRatio);
+  ltermBulk = new nonLocalDG3D1DForceBulk(this,ip,_crossSectionDistribution,_nonLocalEqRatio);
   if(_bmbp)
     btermBulk = new BilinearTermPerturbation<double>(ltermBulk,*dgspace,*dgspace,uf,ip,this,_eps);
   else
-   btermBulk = new nonLocalDG3D1DStiffnessBulk(this,ip,_nonLocalEqRatio);
+   btermBulk = new nonLocalDG3D1DStiffnessBulk(this,ip,_crossSectionDistribution,_nonLocalEqRatio);
 
   
   ltermBound = new nonLinearTermVoid();
diff --git a/dG3D/src/dG3D1DDomain.h b/dG3D/src/dG3D1DDomain.h
index 9386e3ef04ced5f10413e6619fa5ce517d005b14..33808b4f6ffe8774e0524f0dfbaa137f86162e89 100644
--- a/dG3D/src/dG3D1DDomain.h
+++ b/dG3D/src/dG3D1DDomain.h
@@ -13,22 +13,19 @@
 #ifndef DG3D1DDOMAIN_H_
 #define DG3D1DDOMAIN_H_
 #include "dG3DDomain.h"
+#include "scalarFunction.h"
 
 class dG3D1DDomain : public dG3DDomain{
   protected:
     #ifndef SWIG
-    std::map<int,double> _otherCrossSection;
-    double _crossSection; // cross section
+    scalarFunction* _crossSectionDistribution;
     #endif //SWIG
   public:
     dG3D1DDomain(const int tag, const int phys, const int sp_, const int lnum, const int fdg);
-    
-    void setCrossSection(const double A);
-    void setCrossSectionOnElement(const int el, const double A);
+    void setCrossSectionDistribution(const scalarFunction* fct);
     #ifndef SWIG
-    dG3D1DDomain(const dG3D1DDomain& src) : dG3DDomain(src), _otherCrossSection(src._otherCrossSection), _crossSection(src._crossSection){};
-    virtual ~dG3D1DDomain(){}
-    double getCrossSection(const int el) const;
+    dG3D1DDomain(const dG3D1DDomain& src);
+    virtual ~dG3D1DDomain();
     virtual void createTerms(unknownField *uf,IPField*ip);
     virtual partDomain* clone() const {return new dG3D1DDomain(*this);};
     #endif //SWIG
@@ -37,20 +34,16 @@ class dG3D1DDomain : public dG3DDomain{
 class nonLocalDamageDG3D1DDDomain : public nonLocalDamageDG3DDomain{
   protected:
     #ifndef SWIG
-    std::map<int,double> _otherCrossSection;
-    double _crossSection; // cross section
+    scalarFunction* _crossSectionDistribution;
     #endif //SWIG
   
   public:
-    nonLocalDamageDG3D1DDDomain(const int tag, const int phys, const int sp_, const int lnum, const int fdg, const double nonLocalEqRatio,const int numNonlocalVar=1);
-              
-    void setCrossSection(const double A);
-    void setCrossSectionOnElement(const int el, const double A);
+    nonLocalDamageDG3D1DDDomain(const int tag, const int phys, const int sp_, const int lnum, const int fdg, const double nonLocalEqRatio,const int numNonlocalVar);
+    void setCrossSectionDistribution(const scalarFunction* fct);
     
     #ifndef SWIG
-    nonLocalDamageDG3D1DDDomain(const nonLocalDamageDG3D1DDDomain& src): nonLocalDamageDG3DDomain(src),_otherCrossSection(src._otherCrossSection), _crossSection(src._crossSection){};
-    virtual ~nonLocalDamageDG3D1DDDomain(){};
-    double getCrossSection(const int el) const;
+    nonLocalDamageDG3D1DDDomain(const nonLocalDamageDG3D1DDDomain& src);
+    virtual ~nonLocalDamageDG3D1DDDomain();
     virtual void createTerms(unknownField *uf,IPField*ip);
     virtual partDomain* clone() const {return new nonLocalDamageDG3D1DDDomain(*this);};
     #endif //SWIG
diff --git a/dG3D/src/dG3D1DTerms.cpp b/dG3D/src/dG3D1DTerms.cpp
index f996c126dfc9e3d2191a355f4d035277c55e236a..ff9611d533e63531b3686919e91aaa9766f15742 100644
--- a/dG3D/src/dG3D1DTerms.cpp
+++ b/dG3D/src/dG3D1DTerms.cpp
@@ -22,27 +22,25 @@ void dG3D1DForceBulk::get(MElement *ele,int npts,IntPt *GP,fullVector<double> &v
   
   const dG3D1DDomain* dom1D = dynamic_cast<const dG3D1DDomain*>(_dom);
   const nonLocalDamageDG3D1DDDomain* dom1DNonLocal = dynamic_cast<const nonLocalDamageDG3D1DDDomain*>(_dom);
-  double A = 1.;
-  if (dom1D != NULL){
-    A = dom1D->getCrossSection(ele->getNum());
-  }
-  else if (dom1DNonLocal != NULL){
-    A = dom1DNonLocal->getCrossSection(ele->getNum());
-  }
-  
   
   static SVector3 BT;
+  static SPoint3 point;
 
   // get gauss' point data
   _ipf->getIPv(ele,vipv);
-
-
+  _crossSection->setElement(ele);
+  
   for (int i = 0; i < npts; i++)
   {
     double weight = GP[i].weight;
     double & detJ = vipv[i]->getJacobianDeterminant(ele,GP[i]);
 
     const STensor3 *PK1 = &(vipv[i]->getConstRefToFirstPiolaKirchhoffStress());
+    
+    double u = GP[i].pt[0]; double v = GP[i].pt[1];  double w = GP[i].pt[2];
+    ele->pnt(u,v,w,point);
+    double A = _crossSection->getVal(point[0]); // crossection is in fuction of x only
+    
     double ratio = detJ*weight*A;
 
     const std::vector<TensorialTraits<double>::GradType> &Grads = vipv[i]->gradf(space,ele,GP[i]);
@@ -72,25 +70,22 @@ void dG3D1DStiffnessBulk::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double
 
   // get ipvariables
   _ipf->getIPv(ele,vipv);
+  _crossSection->setElement(ele);
 
   // sum on Gauss' points
   static SVector3 B;
   static SVector3 BT;
+  static SPoint3 point;
   
-  const dG3D1DDomain* dom1D = dynamic_cast<const dG3D1DDomain*>(_dom);
-  const nonLocalDamageDG3D1DDDomain* dom1DNonLocal = dynamic_cast<const nonLocalDamageDG3D1DDDomain*>(_dom);
-  double A = 1.;
-  if (dom1D != NULL){
-    A = dom1D->getCrossSection(ele->getNum());
-  }
-  else if (dom1DNonLocal != NULL){
-    A = dom1DNonLocal->getCrossSection(ele->getNum());
-  }
-
   for (int i = 0; i < npts; i++)
   {
     double weight = GP[i].weight;
     double& detJ = vipv[i]->getJacobianDeterminant(ele,GP[i]);
+    
+    double u = GP[i].pt[0]; double v = GP[i].pt[1];  double w = GP[i].pt[2];
+    ele->pnt(u,v,w,point);
+    double A = _crossSection->getVal(point[0]); // crossection is in fuction of x only
+    
     double ratio = detJ*weight*A;
     // x, y, z
     const STensor43 *H =&vipv[i]->getConstRefToTangentModuli();
@@ -133,17 +128,19 @@ void nonLocalDG3D1DForceBulk::get(MElement *ele,int npts,IntPt *GP,fullVector<do
   int nbFF = ele->getNumShapeFunctions();
   int threeTimesNbFF = 3*nbFF;
   
-  const nonLocalDamageDG3D1DDDomain* dom1D = dynamic_cast<const nonLocalDamageDG3D1DDDomain*>(_dom);
-  double A = 1.;
-  if (dom1D != NULL){
-    A = dom1D->getCrossSection(ele->getNum());
-  }
-  
+  static SPoint3 point;
+  _crossSection->setElement(ele);
+
   for (int i = 0; i < npts; i++)
   {
     double weight = GP[i].weight;
     double &detJ = vipv[i]->getJacobianDeterminant(ele,GP[i]);
-    double eqRatio = _nonLocalEqRatio;                    //by Wu Ling
+    double eqRatio = _nonLocalEqRatio;   
+    
+    double u = GP[i].pt[0]; double v = GP[i].pt[1];  double w = GP[i].pt[2];
+    ele->pnt(u,v,w,point);
+    double A = _crossSection->getVal(point[0]); // crossection is in fuction of x only
+    
     double ratio = detJ*weight*A;
 
     const std::vector<TensorialTraits<double>::GradType> &Grads = vipv[i]->gradf(space,ele,GP[i]);
@@ -177,19 +174,22 @@ void nonLocalDG3D1DStiffnessBulk::get(MElement *ele,int npts,IntPt *GP,fullMatri
   const FunctionSpaceBase* space  = _dom->getFunctionSpace();
   int nbdof = space->getNumKeys(ele);
   int nbFF = ele->getNumShapeFunctions();
-  const nonLocalDamageDG3D1DDDomain* dom1D = dynamic_cast<const nonLocalDamageDG3D1DDDomain*>(_dom);
-  double A = 1.;
-  if (dom1D != NULL){
-    A = dom1D->getCrossSection(ele->getNum());
-  }
+  
+  static SPoint3 point;
+  _crossSection->setElement(ele);
   
   for (int i = 0; i < npts; i++)
   {
 
     double weight = GP[i].weight;
     double &detJ = vipv[i]->getJacobianDeterminant(ele,GP[i]);
+    
+    double u = GP[i].pt[0]; double v = GP[i].pt[1];  double w = GP[i].pt[2];
+    ele->pnt(u,v,w,point);
+    double A = _crossSection->getVal(point[0]); // crossection is in fuction of x only
+    
     double ratio = detJ*weight*A;
-    double eqRatio = _nonLocalEqRatio;               //by Wu Ling
+    double eqRatio = _nonLocalEqRatio;        
 
     std::vector<TensorialTraits<double>::GradType> &Grads = vipv[i]->gradf(space,ele,GP[i]);
     std::vector<TensorialTraits<double>::ValType> &Vals = vipv[i]->f(space,ele,GP[i]);
diff --git a/dG3D/src/dG3D1DTerms.h b/dG3D/src/dG3D1DTerms.h
index b2ead47eeb8cdbccadfb2204a90d2c7039d23e2d..9826b9dff8f60a9bb14edcffc30d71fdfe37fbfd 100644
--- a/dG3D/src/dG3D1DTerms.h
+++ b/dG3D/src/dG3D1DTerms.h
@@ -17,16 +17,19 @@
 #include "ipField.h"
 #include "partDomain.h"
 #include "dG3DTerms.h"
+#include "scalarFunction.h"
 
 
 class dG3D1DForceBulk : public nonLinearTermBase<double>{
 	protected:
 		const partDomain* _dom;
     const IPField* _ipf;
+    const scalarFunction* _crossSection;
     mutable const dG3DIPVariableBase* vipv[256]; // max 256 Gauss' point ??
+    
 		
 	public:
-		dG3D1DForceBulk(const partDomain* d, const IPField *ip) : _dom(d),_ipf(ip){}
+		dG3D1DForceBulk(const partDomain* d, const IPField *ip, const scalarFunction* fct) : _dom(d),_ipf(ip),_crossSection(fct){}
 		virtual ~dG3D1DForceBulk(){}
 		virtual void get(MElement *ele,int npts,IntPt *GP,fullVector<double> &m) const;
 		virtual void get(MElement *ele, int npts, IntPt *GP, std::vector<fullVector<double> > &mv) const
@@ -38,7 +41,7 @@ class dG3D1DForceBulk : public nonLinearTermBase<double>{
 		
 		virtual LinearTermBase<double>* clone () const
 		{
-			return new dG3D1DForceBulk(_dom,_ipf);
+			return new dG3D1DForceBulk(_dom,_ipf,_crossSection);
 		}
 };
 
@@ -46,10 +49,12 @@ class dG3D1DStiffnessBulk : public BiNonLinearTermBase {
 	protected:
 		const partDomain* _dom;
     const IPField* _ipf;
+    const scalarFunction* _crossSection;
+    
     mutable const dG3DIPVariableBase* vipv[256]; // max 256 Gauss' point ??
 		
 	public:
-		dG3D1DStiffnessBulk(const partDomain* d, const IPField *ip) : _dom(d),_ipf(ip){}
+		dG3D1DStiffnessBulk(const partDomain* d, const IPField *ip, const scalarFunction* fct) : _dom(d),_ipf(ip),_crossSection(fct){}
 		virtual ~dG3D1DStiffnessBulk(){}									
 		virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m) const;
     
@@ -62,7 +67,7 @@ class dG3D1DStiffnessBulk : public BiNonLinearTermBase {
 		
 		virtual BilinearTermBase* clone () const
 		{
-			return new dG3D1DStiffnessBulk(_dom,_ipf);
+			return new dG3D1DStiffnessBulk(_dom,_ipf,_crossSection);
 		}
 };
 
@@ -71,7 +76,7 @@ class nonLocalDG3D1DForceBulk : public dG3D1DForceBulk{
     double _nonLocalEqRatio;
     
   public:
-		nonLocalDG3D1DForceBulk(const partDomain* d, const IPField *ip, const double eqR) : dG3D1DForceBulk(d,ip), _nonLocalEqRatio(eqR){}
+		nonLocalDG3D1DForceBulk(const partDomain* d, const IPField *ip, const scalarFunction* fct, const double eqR) : dG3D1DForceBulk(d,ip,fct), _nonLocalEqRatio(eqR){}
 		virtual ~nonLocalDG3D1DForceBulk(){}
 		virtual void get(MElement *ele,int npts,IntPt *GP,fullVector<double> &m) const;
 		virtual void get(MElement *ele, int npts, IntPt *GP, std::vector<fullVector<double> > &mv) const
@@ -83,7 +88,7 @@ class nonLocalDG3D1DForceBulk : public dG3D1DForceBulk{
 		
 		virtual LinearTermBase<double>* clone () const
 		{
-			return new nonLocalDG3D1DForceBulk(_dom,_ipf,_nonLocalEqRatio);
+			return new nonLocalDG3D1DForceBulk(_dom,_ipf,_crossSection,_nonLocalEqRatio);
 		}
 };
 
@@ -92,7 +97,7 @@ class nonLocalDG3D1DStiffnessBulk : public dG3D1DStiffnessBulk{
     double _nonLocalEqRatio;
     
   public:
-    nonLocalDG3D1DStiffnessBulk(const partDomain* d, const IPField *ip, const double eq) : dG3D1DStiffnessBulk(d,ip),_nonLocalEqRatio(eq){}
+    nonLocalDG3D1DStiffnessBulk(const partDomain* d, const IPField *ip, const scalarFunction* fct, const double eq) : dG3D1DStiffnessBulk(d,ip,fct),_nonLocalEqRatio(eq){}
 		virtual ~nonLocalDG3D1DStiffnessBulk(){}									
 		virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m) const;
     
@@ -105,7 +110,7 @@ class nonLocalDG3D1DStiffnessBulk : public dG3D1DStiffnessBulk{
 		
 		virtual BilinearTermBase* clone () const
 		{
-			return new nonLocalDG3D1DStiffnessBulk(_dom,_ipf,_nonLocalEqRatio);
+			return new nonLocalDG3D1DStiffnessBulk(_dom,_ipf,_crossSection,_nonLocalEqRatio);
 		}
 };