diff --git a/NonLinearSolver/materialLaw/mlawLinearElecMagTherMech.cpp b/NonLinearSolver/materialLaw/mlawLinearElecMagTherMech.cpp
index 2a1febb72c4d51d531b85bf828f59f7a4673c544..e170d8394766a9b6fe16396983255ef19db408ef 100644
--- a/NonLinearSolver/materialLaw/mlawLinearElecMagTherMech.cpp
+++ b/NonLinearSolver/materialLaw/mlawLinearElecMagTherMech.cpp
@@ -225,7 +225,14 @@ void mlawLinearElecMagTherMech::constitutive(
                                             SVector3 & dmechFieldSourcedT,
                                             STensor3 & dmechFieldSourcedgradT,
                                             SVector3 & dmechFieldSourcedV,
-                                            STensor3 & dmechFieldSourcedgradV
+                                            STensor3 & dmechFieldSourcedgradV,
+                                            SVector3 & demFieldSourcedA,
+                                            SVector3 & demFieldSourcedB,
+                                            STensor3 & demFieldSourcedF,
+                                            double & demFieldSourcedT,
+                                            SVector3 & demFieldSourcedGradT,
+                                            double & demFieldSourcedV,
+                                            SVector3 & demFieldSourcedGradV
 
                            )const
 {
@@ -272,6 +279,13 @@ void mlawLinearElecMagTherMech::constitutive(
     STensorOperation::zero(dsourceVectorFielddF);
     STensorOperation::zero(djydA);
     STensorOperation::zero(djydB);
+    STensorOperation::zero(demFieldSourcedA);
+    STensorOperation::zero(demFieldSourcedB);
+    STensorOperation::zero(demFieldSourcedF);
+    STensorOperation::zero(demFieldSourcedT);
+    STensorOperation::zero(demFieldSourcedGradT);
+    STensorOperation::zero(demFieldSourcedV);
+    STensorOperation::zero(demFieldSourcedGradV);
 
     static STensor3 dldT,dkdT;
     double dseebeckdT = 0.0;
@@ -343,7 +357,7 @@ void mlawLinearElecMagTherMech::constitutive(
                     dfluxTdA(i, j) += _l0(i, j) * _seebeck * T * invdh;
 
                     // derivatives of fluxjy
-                    djydT(i) += (_seebeck * _l0(i, j) + dseebeckdT * T * _l0(i, j)) * (An(j) - A0(j)) * invdh +
+                    djydT(i) += (_seebeck * _l0(i, j) + dseebeckdT * T * _l0(i, j) + _seebeck * T * dldT(i,j) + V * dldT(i,j)) * (An(j) - A0(j)) * invdh +
                                 (_seebeck + dseebeckdT * T) * js0(i);
                     djydV(i) += _l0(i, j) * (An(j) - A0(j)) * invdh + js0(i);
                     djydA(i, j) += (_l0(i, j) * V + _l0(i, j) * _seebeck * T) * invdh;
@@ -359,9 +373,10 @@ void mlawLinearElecMagTherMech::constitutive(
                     dsourceVectorFielddA(i, j) += -_l0(i, j) * invdh;
                 }
             }
+            // -j_e
             sourceVectorField(i) = -fluxD(i) - js0(i);
-            // w_AV = j_e \cdot e
-            emFieldSource += sourceVectorField(i) * (-gradV(i) - (An(i) - A0(i)) * invdh);
+            // w_AV = j_e \cdot -dadt + h \cdot dbdt
+            emFieldSource += (sourceVectorField(i) * ((An(i) - A0(i)) * invdh) /*+ H(i) * (Bn(i) - B0(i)) * invdh*/);
         }
 
         // Need to check contribution to jy from source js0
@@ -398,9 +413,10 @@ void mlawLinearElecMagTherMech::constitutive(
                     dsourceVectorFielddgradT(i, j) += 0.0;
                 }
             }
+            // -j_e
             sourceVectorField(i) = -fluxD(i) - js0(i);
-            // w_AV = j_e \cdot e
-            emFieldSource += sourceVectorField(i) * (-gradV(i) - (An(i) - A0(i)) * invdh);
+            // w_AV = j_e \cdot -dadt + h \cdot dbdt
+            emFieldSource += (sourceVectorField(i) * ((An(i) - A0(i)) * invdh) /*+ H(i) * (Bn(i) - B0(i)) * invdh*/);
         }
 
         // For test with unit cube to check
@@ -448,7 +464,19 @@ void mlawLinearElecMagTherMech::constitutive(
 
     dw_TdV += 0.0;
     for (unsigned  int m = 0; m < 3; ++m)
+    {
         dw_TdA(m) += 0.0;
+        for (unsigned int n = 0; n < 3; ++n)
+        {
+            demFieldSourcedA(m) += (-dfluxDdA(m,n) * ((An(n) - A0(n)) * invdh)) - (fluxD(m) * invdh);
+            demFieldSourcedB(m) += (/*dHdB(m,n) * (Bn(n) - B0(n)) * invdh +*/ H(m) * invdh);
+            demFieldSourcedGradT(m) +=  -dedgradT(m,n) * (An(n) - A0(n)) * invdh;
+            demFieldSourcedGradV(m) += -dedgradV(m,n) * (An(n) - A0(n)) * invdh;
+        }
+        demFieldSourcedT += -dedT(m) * (An(m) - A0(m)) * invdh;
+        demFieldSourcedV += 0.0;
+    }
+
     }
 
     // Later think of flux terms and higher order terms d*d*
diff --git a/NonLinearSolver/materialLaw/mlawLinearElecMagTherMech.h b/NonLinearSolver/materialLaw/mlawLinearElecMagTherMech.h
index 1700068243258b0115b167fedb034312af67e24d..de0ee39def5a7c10e188f52e37b703b16ee5d98b 100644
--- a/NonLinearSolver/materialLaw/mlawLinearElecMagTherMech.h
+++ b/NonLinearSolver/materialLaw/mlawLinearElecMagTherMech.h
@@ -172,7 +172,14 @@ class mlawLinearElecMagTherMech : public mlawLinearElecTherMech
                             SVector3 & dmechFieldSourcedT,
                             STensor3 & dmechFieldSourcedgradT,
                             SVector3 & dmechFieldSourcedV,
-                            STensor3 & dmechFieldSourcedgradV
+                            STensor3 & dmechFieldSourcedgradV,
+                            SVector3 & demFieldSourcedA,
+                            SVector3 & demFieldSourcedB,
+                            STensor3 & demFieldSourcedF,
+                            double & demFieldSourcedT,
+                            SVector3 & demFieldSourcedGradT,
+                            double & demFieldSourcedV,
+                            SVector3 & demFieldSourcedGradV
                             ) const;
     double getInitialTemperature() const {return _t0;}
     double getInitialVoltage()     const {return _v0;}
diff --git a/dG3D/src/FractureCohesiveDG3DIPVariable.h b/dG3D/src/FractureCohesiveDG3DIPVariable.h
index 318e0c8b7dc13545613970da9dffef0b28ac6a3d..abe85b4b3b74947533a2e4b2fd4b2f4d92257661 100644
--- a/dG3D/src/FractureCohesiveDG3DIPVariable.h
+++ b/dG3D/src/FractureCohesiveDG3DIPVariable.h
@@ -264,6 +264,16 @@ class FractureCohesive3DIPVariable : public IPVariable2ForFracture<dG3DIPVariabl
   virtual SVector3 & getRefTodSourceVectorFielddt(const int idex) {return _ipvbulk->getRefTodSourceVectorFielddt(idex);}
   virtual const fullVector<double> & getConstRefToEMFieldSource() const {return _ipvbulk->getConstRefToEMFieldSource();}
   virtual fullVector<double> & getRefToEMFieldSource() {return _ipvbulk->getRefToEMFieldSource();}
+  virtual const std::vector< std::vector< SVector3 > > & getConstRefTodEMFieldSourcedVectorPotential() const {return _ipvbulk->getConstRefTodEMFieldSourcedVectorPotential();}
+  virtual std::vector< std::vector< SVector3 > > & getRefTodEMFieldSourcedVectorPotential() {return _ipvbulk->getRefTodEMFieldSourcedVectorPotential();}
+  virtual const std::vector< std::vector< SVector3 > > & getConstRefTodEMFieldSourcedVectorCurl() const {return _ipvbulk->getConstRefTodEMFieldSourcedVectorCurl();}
+  virtual std::vector< std::vector< SVector3 > > & getRefTodEMFieldSourcedVectorCurl() {return _ipvbulk->getRefTodEMFieldSourcedVectorCurl();}
+  virtual const std::vector< STensor3 > & getConstRefTodEMFieldSourcedF() const { return _ipvbulk->getConstRefTodEMFieldSourcedF();}
+  virtual std::vector< STensor3 > & getRefTodEMFieldSourcedF() { return _ipvbulk->getRefTodEMFieldSourcedF();}
+  virtual const fullMatrix<double> & getConstRefTodEMFieldSourcedField() const {return _ipvbulk->getConstRefTodEMFieldSourcedField();}
+  virtual fullMatrix<double> & getRefTodEMFieldSourcedField() {return _ipvbulk->getRefTodEMFieldSourcedField();}
+  virtual const std::vector< std::vector< SVector3 > > & getConstRefTodEMFieldSourcedGradField() const {return _ipvbulk->getConstRefTodEMFieldSourcedGradField();}
+  virtual std::vector< std::vector< SVector3 > > & getRefTodEMFieldSourcedGradField() {return _ipvbulk->getRefTodEMFieldSourcedGradField();}
 
   //Energy Conjugated Field (fT & fv)
   virtual double getConstRefToEnergyConjugatedField(const int idex) const {return _ipvbulk->getConstRefToEnergyConjugatedField(idex);};
diff --git a/dG3D/src/dG3DIPVariable.cpp b/dG3D/src/dG3DIPVariable.cpp
index c32a1655c90e281e6e25f491350fc457c6641328..38ae0d2979f77f391e436407d95e5408cd39274a 100644
--- a/dG3D/src/dG3DIPVariable.cpp
+++ b/dG3D/src/dG3DIPVariable.cpp
@@ -370,7 +370,10 @@ constitutiveCurlData::constitutiveCurlData(const int numCurlVar, const int numNo
         inductorSourceVectorField.resize(numCurlVar, SVector3(0.0));
         sourceVectorField.resize(numCurlVar, SVector3(0.0));
         dSourceVectorFielddt.resize(numCurlVar, SVector3(0.0));
-        emFieldSource.resize(numCurlVar); emFieldSource.setAll(0.0);
+        emFieldSource.resize(numExtraDofDiffusion); emFieldSource.setAll(0.0);
+
+        std::vector< SVector3 > StdSVec3Tmp;
+        StdSVec3Tmp.resize(numCurlVar, SVector3(0.0));
 
         std::vector< STensor3 > StdSTen3Tmp;
         StdSTen3Tmp.resize(numCurlVar, STensor3(0.0));
@@ -386,6 +389,11 @@ constitutiveCurlData::constitutiveCurlData(const int numCurlVar, const int numNo
         std::vector< STensor3 > StdSTen3TmpCouplingNonLocal;
         StdSTen3TmpCouplingNonLocal.resize(numNonLocal, STensor3(0.0));
 
+        demFieldSourcedVectorPotential.resize(numExtraDofDiffusion, StdSVec3Tmp);
+        demFieldSourcedVectorCurl.resize(numExtraDofDiffusion, StdSVec3Tmp);
+        demFieldSourcedF.resize(numExtraDofDiffusion, STensor3(0.0));
+        demFieldSourcedField.resize(numExtraDofDiffusion, numExtraDofDiffusion); demFieldSourcedField.setAll(0.0);
+        demFieldSourcedGradField.resize(numExtraDofDiffusion, StdSVec3Tmp);
         dVectorFielddVectorPotential.resize(numCurlVar, StdSTen3Tmp);
         dVectorFielddVectorCurl.resize(numCurlVar, StdSTen3Tmp);
         dPdVectorCurl.resize(numCurlVar, STensor33(0.0));
@@ -417,6 +425,8 @@ constitutiveCurlData::constitutiveCurlData(const constitutiveCurlData& src)
   emFieldSource(src.emFieldSource),
   dSourceVectorFielddVectorCurl(src.dSourceVectorFielddVectorCurl),
   dSourceVectorFielddVectorPotential(src.dSourceVectorFielddVectorPotential),
+  demFieldSourcedVectorPotential(src.demFieldSourcedVectorPotential),
+  demFieldSourcedVectorCurl(src.demFieldSourcedVectorCurl),
   dVectorFielddVectorPotential(src.dVectorFielddVectorPotential),
   dVectorFielddVectorCurl(src.dVectorFielddVectorCurl),
   dPdVectorCurl(src.dPdVectorCurl),
@@ -425,12 +435,15 @@ constitutiveCurlData::constitutiveCurlData(const constitutiveCurlData& src)
   mechanicalFieldSource(src.mechanicalFieldSource),
   dMechanicalFieldSourcedVectorCurl(src.dMechanicalFieldSourcedVectorCurl),
   dMechanicalFieldSourcedF(src.dMechanicalFieldSourcedF),
+  demFieldSourcedF(src.demFieldSourcedF),
   dVectorFielddExtraDofField(src.dVectorFielddExtraDofField),
   dVectorFielddGradExtraDofField(src.dVectorFielddGradExtraDofField),
   dSourceVectorFielddExtraDofField(src.dSourceVectorFielddExtraDofField),
   dSourceVectorFielddGradExtraDofField(src.dSourceVectorFielddGradExtraDofField),
   dMechanicalFieldSourcedExtraDofField(src.dMechanicalFieldSourcedExtraDofField),
   dMechanicalFieldSourcedGradExtraDofField(src.dMechanicalFieldSourcedGradExtraDofField),
+  demFieldSourcedField(src.demFieldSourcedField),
+  demFieldSourcedGradField(src.demFieldSourcedGradField),
   dVectorFielddNonLocalVariable(src.dVectorFielddNonLocalVariable),
   dSourceVectorFielddNonLocalVariable(src.dSourceVectorFielddNonLocalVariable),
   dMechanicalFieldSourcedNonLocalVariable(src.dMechanicalFieldSourcedNonLocalVariable)
@@ -452,6 +465,8 @@ constitutiveCurlData & constitutiveCurlData::operator=(const constitutiveCurlDat
     emFieldSource = src.emFieldSource;
     dSourceVectorFielddVectorCurl = src.dSourceVectorFielddVectorCurl;
     dSourceVectorFielddVectorPotential = src.dSourceVectorFielddVectorPotential;
+    demFieldSourcedVectorPotential = src.demFieldSourcedVectorPotential;
+    demFieldSourcedVectorCurl = src.demFieldSourcedVectorCurl;
     dVectorFielddVectorPotential = src.dVectorFielddVectorPotential;
     dVectorFielddVectorCurl = src.dVectorFielddVectorCurl;
     dPdVectorCurl = src.dPdVectorCurl;
@@ -460,12 +475,15 @@ constitutiveCurlData & constitutiveCurlData::operator=(const constitutiveCurlDat
     mechanicalFieldSource = src.mechanicalFieldSource;
     dMechanicalFieldSourcedVectorCurl = src.dMechanicalFieldSourcedVectorCurl;
     dMechanicalFieldSourcedF = src.dMechanicalFieldSourcedF;
+    demFieldSourcedF = src.demFieldSourcedF;
     dVectorFielddExtraDofField = src.dVectorFielddExtraDofField;
     dVectorFielddGradExtraDofField = src.dVectorFielddGradExtraDofField;
     dSourceVectorFielddExtraDofField = src.dSourceVectorFielddExtraDofField;
     dSourceVectorFielddGradExtraDofField = src.dSourceVectorFielddGradExtraDofField;
     dMechanicalFieldSourcedExtraDofField = src.dMechanicalFieldSourcedExtraDofField;
     dMechanicalFieldSourcedGradExtraDofField = src.dMechanicalFieldSourcedGradExtraDofField;
+    demFieldSourcedField = src.demFieldSourcedField;
+    demFieldSourcedGradField = src.demFieldSourcedGradField;
     dVectorFielddNonLocalVariable = src.dVectorFielddNonLocalVariable;
     dSourceVectorFielddNonLocalVariable = src.dSourceVectorFielddNonLocalVariable;
     dMechanicalFieldSourcedNonLocalVariable = src.dMechanicalFieldSourcedNonLocalVariable;
@@ -484,9 +502,11 @@ void constitutiveCurlData::restart()
     restartManager::restart(inductorSourceVectorField);
     restartManager::restart(sourceVectorField);
     restartManager::restart(dSourceVectorFielddt);
-    restartManager::restart(emFieldSource.getDataPtr(),numConstitutiveCurlVariable);
+    restartManager::restart(emFieldSource.getDataPtr(),numConstitutiveExtraDofDiffusionVariable);
     restartManager::restart(dSourceVectorFielddVectorCurl);
     restartManager::restart(dSourceVectorFielddVectorPotential);
+    restartManager::restart(demFieldSourcedVectorPotential);
+    restartManager::restart(demFieldSourcedVectorCurl);
     restartManager::restart(dVectorFielddVectorPotential);
     restartManager::restart(dVectorFielddVectorCurl);
     restartManager::restart(dVectorFielddF);
@@ -495,12 +515,15 @@ void constitutiveCurlData::restart()
     restartManager::restart(mechanicalFieldSource);
     restartManager::restart(dMechanicalFieldSourcedVectorCurl);
     restartManager::restart(dMechanicalFieldSourcedF);
+    restartManager::restart(demFieldSourcedF);
     restartManager::restart(dVectorFielddGradExtraDofField);
     restartManager::restart(dVectorFielddNonLocalVariable);
     restartManager::restart(dSourceVectorFielddExtraDofField);
     restartManager::restart(dSourceVectorFielddGradExtraDofField);
     restartManager::restart(dMechanicalFieldSourcedExtraDofField);
     restartManager::restart(dMechanicalFieldSourcedGradExtraDofField);
+    restartManager::restart(demFieldSourcedField.getDataPtr(),numConstitutiveExtraDofDiffusionVariable*numConstitutiveExtraDofDiffusionVariable);
+    restartManager::restart(demFieldSourcedGradField);
     restartManager::restart(dVectorFielddExtraDofField);
     restartManager::restart(dSourceVectorFielddNonLocalVariable);
     restartManager::restart(dMechanicalFieldSourcedNonLocalVariable);
diff --git a/dG3D/src/dG3DIPVariable.h b/dG3D/src/dG3DIPVariable.h
index ecc0933cb8ab96b77d10b7d5246e8ebedc70118f..bcdee4f3a0bb26aa1984ba1f2d2855599a3faa63 100644
--- a/dG3D/src/dG3DIPVariable.h
+++ b/dG3D/src/dG3DIPVariable.h
@@ -316,6 +316,16 @@ class dG3DIPVariableBase : public ipFiniteStrain
   virtual std::vector< std::vector< STensor3 > > & getRefTodSourceVectorFielddVectorCurl()=0;
   virtual const std::vector< std::vector< STensor3 > > & getConstRefTodSourceVectorFielddVectorPotential() const=0;
   virtual std::vector< std::vector< STensor3 > > & getRefTodSourceVectorFielddVectorPotential()=0;
+  virtual const std::vector< std::vector< SVector3 > > & getConstRefTodEMFieldSourcedVectorPotential() const=0;
+  virtual std::vector< std::vector< SVector3 > > & getRefTodEMFieldSourcedVectorPotential()=0;
+  virtual const std::vector< std::vector< SVector3 > > & getConstRefTodEMFieldSourcedVectorCurl() const=0;
+  virtual std::vector< std::vector< SVector3 > > & getRefTodEMFieldSourcedVectorCurl()=0;
+  virtual const std::vector< STensor3 > & getConstRefTodEMFieldSourcedF() const=0;
+  virtual std::vector< STensor3 > & getRefTodEMFieldSourcedF()=0;
+  virtual const fullMatrix<double> & getConstRefTodEMFieldSourcedField() const=0;
+  virtual fullMatrix<double> & getRefTodEMFieldSourcedField()=0;
+  virtual const std::vector< std::vector< SVector3 > > & getConstRefTodEMFieldSourcedGradField() const=0;
+  virtual std::vector< std::vector< SVector3 > > & getRefTodEMFieldSourcedGradField()=0;
 
   //Energy Conjugated Field (fT & fv)
 
@@ -474,9 +484,11 @@ class constitutiveCurlData{
   std::vector< SVector3 >           inductorSourceVectorField; // current density applied in inductor js0
   std::vector< SVector3 >           sourceVectorField; // Vector magnetic field source (T, grad V, F)
   std::vector< SVector3 >           dSourceVectorFielddt; // time derivative
-  fullVector<double>                emFieldSource; // w_AV = j_e \cdot e
+  fullVector<double>                emFieldSource; // w_AV Eddy current + hysteresis loss
   std::vector< std::vector< STensor3 > > dSourceVectorFielddVectorCurl;
   std::vector< std::vector< STensor3 > > dSourceVectorFielddVectorPotential;
+  std::vector< std::vector< SVector3 > > demFieldSourcedVectorPotential;
+  std::vector< std::vector< SVector3 > > demFieldSourcedVectorCurl;
 
   std::vector<std::vector< STensor3 > > dVectorFielddVectorPotential;
   std::vector<std::vector< STensor3 > > dVectorFielddVectorCurl;
@@ -486,6 +498,7 @@ class constitutiveCurlData{
   std::vector< STensor33 > dVectorFielddF;
   std::vector< STensor33 >  dSourceVectorFielddF;
   std::vector< STensor33 >  dMechanicalFieldSourcedF;
+  std::vector< STensor3 >  demFieldSourcedF;
 
   std::vector< SVector3 >  mechanicalFieldSource;   // EM force proportional to freq (V, B, F)
   std::vector< std::vector< STensor3 > > dMechanicalFieldSourcedVectorCurl;
@@ -497,6 +510,8 @@ class constitutiveCurlData{
   std::vector< std::vector< STensor3 > > dSourceVectorFielddGradExtraDofField;
   std::vector< std::vector< SVector3 > > dMechanicalFieldSourcedExtraDofField;
   std::vector< std::vector< STensor3 > > dMechanicalFieldSourcedGradExtraDofField;
+  fullMatrix< double >                demFieldSourcedField;
+  std::vector< std::vector< SVector3 > > demFieldSourcedGradField;
 
   //coupling with nonLocal
   std::vector<std::vector< SVector3 > >   dVectorFielddNonLocalVariable;
@@ -1441,6 +1456,76 @@ protected:
       return _constitutiveCurlData->emFieldSource;
   }
 
+  virtual const std::vector< std::vector< SVector3 > > & getConstRefTodEMFieldSourcedVectorPotential() const
+  {
+      if (_constitutiveCurlData == NULL)
+          Msg::Error("getConstRefTodEMFieldSourcedVectorPotential(): _constitutiveCurlData not created");
+      return _constitutiveCurlData->demFieldSourcedVectorPotential;
+  }
+
+    virtual std::vector< std::vector< SVector3 > > & getRefTodEMFieldSourcedVectorPotential()
+    {
+        if (_constitutiveCurlData == NULL)
+            Msg::Error("getRefTodEMFieldSourcedVectorPotential(): _constitutiveCurlData not created");
+        return _constitutiveCurlData->demFieldSourcedVectorPotential;
+    }
+
+    virtual const std::vector< std::vector< SVector3 > > & getConstRefTodEMFieldSourcedVectorCurl() const
+    {
+        if (_constitutiveCurlData == NULL)
+            Msg::Error("getConstRefTodEMFieldSourcedVectorCurl(): _constitutiveCurlData not created");
+        return _constitutiveCurlData->demFieldSourcedVectorCurl;
+    }
+
+    virtual std::vector< std::vector< SVector3 > > & getRefTodEMFieldSourcedVectorCurl()
+    {
+        if (_constitutiveCurlData == NULL)
+            Msg::Error("getRefTodEMFieldSourcedVectorCurl(): _constitutiveCurlData not created");
+        return _constitutiveCurlData->demFieldSourcedVectorCurl;
+    }
+
+    virtual const std::vector< STensor3 > & getConstRefTodEMFieldSourcedF() const
+    {
+        if (_constitutiveCurlData == NULL)
+            Msg::Error("getConstRefTodEMFieldSourcedF(): _constitutiveCurlData not created");
+        return _constitutiveCurlData->demFieldSourcedF;
+    }
+
+    virtual std::vector< STensor3 > & getRefTodEMFieldSourcedF()
+    {
+        if (_constitutiveCurlData == NULL)
+            Msg::Error("getRefTodEMFieldSourcedF(): _constitutiveCurlData not created");
+        return _constitutiveCurlData->demFieldSourcedF;
+    }
+
+    virtual const fullMatrix<double> & getConstRefTodEMFieldSourcedField() const
+    {
+        if (_constitutiveCurlData == NULL)
+            Msg::Error("getConstRefTodEMFieldSourcedField(): _constitutiveCurlData not created");
+        return _constitutiveCurlData->demFieldSourcedField;
+    }
+
+    virtual fullMatrix<double> & getRefTodEMFieldSourcedField()
+    {
+        if (_constitutiveCurlData == NULL)
+            Msg::Error("getRefTodEMFieldSourcedField(): _constitutiveCurlData not created");
+        return _constitutiveCurlData->demFieldSourcedField;
+    }
+
+    virtual const std::vector< std::vector< SVector3 > > & getConstRefTodEMFieldSourcedGradField() const
+    {
+        if (_constitutiveCurlData == NULL)
+            Msg::Error("getConstRefTodEMFieldSourcedGradField(): _constitutiveCurlData not created");
+        return _constitutiveCurlData->demFieldSourcedGradField;
+    }
+
+    virtual std::vector< std::vector< SVector3 > > & getRefTodEMFieldSourcedGradField()
+    {
+        if (_constitutiveCurlData == NULL)
+            Msg::Error("getRefTodEMFieldSourcedGradField(): _constitutiveCurlData not created");
+        return _constitutiveCurlData->demFieldSourcedGradField;
+    }
+
     virtual const std::vector< std::vector< STensor3 > > & getConstRefTodSourceVectorFielddVectorPotential() const
     {
         if (_constitutiveCurlData == NULL)
@@ -3947,6 +4032,76 @@ class ElecMagTherMechDG3DIPVariableBase : public dG3DIPVariable
         return getRefToEMFieldSource()(0);
     }
 
+    virtual const SVector3 & getConstRefTodemFieldSourcedMagneticVectorPotential() const
+    {
+        return getConstRefTodEMFieldSourcedVectorPotential()[0][0];
+    }
+
+    virtual SVector3 & getRefTodemFieldSourcedMagneticVectorPotential()
+    {
+        return getRefTodEMFieldSourcedVectorPotential()[0][0];
+    }
+
+    virtual const SVector3 & getConstRefTodemFieldSourcedMagneticVectorCurl() const
+    {
+        return getConstRefTodEMFieldSourcedVectorCurl()[0][0];
+    }
+
+    virtual SVector3 & getRefTodemFieldSourcedMagneticVectorCurl()
+    {
+        return getRefTodEMFieldSourcedVectorCurl()[0][0];
+    }
+
+    virtual const STensor3 & getConstRefTodemFieldSourcedF() const
+    {
+        return getConstRefTodEMFieldSourcedF()[0];
+    }
+
+    virtual STensor3 & getRefTodemFieldSourcedF()
+    {
+        return getRefTodEMFieldSourcedF()[0];
+    }
+
+    virtual const double getConstRefTodemFieldSourcedT() const
+    {
+        return getConstRefTodEMFieldSourcedField()(0,0);
+    }
+
+    virtual double & getRefTodemFieldSourcedT()
+    {
+        return getRefTodEMFieldSourcedField()(0,0);
+    }
+
+    virtual const double getConstRefTodemFieldSourcedV() const
+    {
+        return getConstRefTodEMFieldSourcedField()(0,1);
+    }
+
+    virtual double & getRefTodemFieldSourcedV()
+    {
+        return getRefTodEMFieldSourcedField()(0,1);
+    }
+
+    virtual const SVector3 & getConstRefTodemFieldSourcedGradT() const
+    {
+        return getConstRefTodEMFieldSourcedGradField()[0][0];
+    }
+
+    virtual SVector3 & getRefTodemFieldSourcedGradT()
+    {
+        return getRefTodEMFieldSourcedGradField()[0][0];
+    }
+
+    virtual const SVector3 & getConstRefTodemFieldSourcedGradV() const
+    {
+        return getConstRefTodEMFieldSourcedGradField()[0][1];
+    }
+
+    virtual SVector3 & getRefTodemFieldSourcedGradV()
+    {
+        return getRefTodEMFieldSourcedGradField()[0][1];
+    }
+
     // dsourceVectorFielddA
     virtual const STensor3 & getConstRefTodSourceVectorFielddMagneticVectorPotential() const
     {
diff --git a/dG3D/src/dG3DMaterialLaw.cpp b/dG3D/src/dG3DMaterialLaw.cpp
index be0294421159a75785113cd48bb318370e3270c3..9568a3fd50ea8e8900143f1686b3049504304541 100644
--- a/dG3D/src/dG3DMaterialLaw.cpp
+++ b/dG3D/src/dG3DMaterialLaw.cpp
@@ -5766,6 +5766,13 @@ const bool checkfrac, const bool dTangent)
     STensor3 & dmechFieldSourcedgradT = ipvcur->getRefTodMechanicalFieldSourcedGradT();
     SVector3 & dmechFieldSourcedV = ipvcur->getRefTodMechanicalFieldSourcedV();
     STensor3 & dmechFieldSourcedgradV = ipvcur->getRefTodMechanicalFieldSourcedGradV();
+    SVector3 & demFieldSourcedA = ipvcur->getRefTodemFieldSourcedMagneticVectorPotential();
+    SVector3 & demFieldSourcedB = ipvcur->getRefTodemFieldSourcedMagneticVectorCurl();
+    STensor3 & demFieldSourcedF = ipvcur->getRefTodemFieldSourcedF();
+    double & demFieldSourcedT = ipvcur->getRefTodemFieldSourcedT();
+    SVector3 & demFieldSourcedGradT = ipvcur->getRefTodemFieldSourcedGradT();
+    double & demFieldSourcedV = ipvcur->getRefTodemFieldSourcedV();
+    SVector3 & demFieldSourcedGradV = ipvcur->getRefTodemFieldSourcedGradV();
 
     _lawLinearEMTM->constitutive(F0,Fn,P, q0, q1, Tangent,T0, temperature, voltage,gradT, gradV,
     fluxT,fluxD, dPdT,dPdV, dqdgradT,dqdT, dqdgradV,dqdV,dedV,dedgradV,dedT,dedgradT,dqdF,dedF,
@@ -5775,7 +5782,8 @@ const bool checkfrac, const bool dTangent)
     A0, An, B, H, js0, dBdT, dBdGradT, dBdV, dBdGradV, dBdF, dBdA, dfluxTdA, dfluxDdA, dfluxTdB, dfluxDdB, dw_TdA, dw_TdB, dmechSourcedB, sourceVectorField,
     dsourceVectorFielddt, emFieldSource, dHdA, dHdB, dPdB, dHdF, dsourceVectorFielddF, mechFieldSource, dmechFieldSourcedB, dmechFieldSourcedF,
     dHdT, dHdgradT, dHdV, dHdgradV, dsourceVectorFielddT, dsourceVectorFielddgradT, dsourceVectorFielddV, dsourceVectorFielddgradV, dsourceVectorFielddA,
-    dsourceVectorFielddB, dmechFieldSourcedT, dmechFieldSourcedgradT, dmechFieldSourcedV, dmechFieldSourcedgradV);
+    dsourceVectorFielddB, dmechFieldSourcedT, dmechFieldSourcedgradT, dmechFieldSourcedV, dmechFieldSourcedgradV,
+    demFieldSourcedA, demFieldSourcedB, demFieldSourcedF, demFieldSourcedT, demFieldSourcedGradT, demFieldSourcedV, demFieldSourcedGradV);
 
     ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
 }
@@ -6013,6 +6021,13 @@ const bool checkfrac, const bool dTangent)
   STensor3 & dmechFieldSourcedgradT = ipvcur->getRefTodMechanicalFieldSourcedGradT();
   SVector3 & dmechFieldSourcedV = ipvcur->getRefTodMechanicalFieldSourcedV();
   STensor3 & dmechFieldSourcedgradV = ipvcur->getRefTodMechanicalFieldSourcedGradV();
+  SVector3 & demFieldSourcedA = ipvcur->getRefTodemFieldSourcedMagneticVectorPotential();
+  SVector3 & demFieldSourcedB = ipvcur->getRefTodemFieldSourcedMagneticVectorCurl();
+  STensor3 & demFieldSourcedF = ipvcur->getRefTodemFieldSourcedF();
+  double & demFieldSourcedT = ipvcur->getRefTodemFieldSourcedT();
+  SVector3 & demFieldSourcedGradT = ipvcur->getRefTodemFieldSourcedGradT();
+  double & demFieldSourcedV = ipvcur->getRefTodemFieldSourcedV();
+  SVector3 & demFieldSourcedGradV = ipvcur->getRefTodemFieldSourcedGradV();
   
   _lawLinearEMInductor->constitutive(F0,Fn,P, q0, q1, Tangent,T0, temperature, voltage,gradT, gradV,
     fluxT,fluxD, dPdT,dPdV, dqdgradT,dqdT, dqdgradV,dqdV,dedV,dedgradV,dedT,dedgradT,dqdF,dedF,
@@ -6022,7 +6037,8 @@ const bool checkfrac, const bool dTangent)
     A0, An, B, H, js0, dBdT, dBdGradT, dBdV, dBdGradV, dBdF, dBdA, dfluxTdA, dfluxDdA, dfluxTdB, dfluxDdB, dw_TdA, dw_TdB, dmechSourcedB,
     sourceVectorField, dsourceVectorFielddt, emFieldSource, dHdA, dHdB, dPdB, dHdF, dsourceVectorFielddF, mechFieldSource, dmechFieldSourcedB, dmechFieldSourcedF,
     dHdT, dHdgradT, dHdV, dHdgradV, dsourceVectorFielddT, dsourceVectorFielddgradT, dsourceVectorFielddV, dsourceVectorFielddgradV, dsourceVectorFielddA,
-    dsourceVectorFielddB, dmechFieldSourcedT, dmechFieldSourcedgradT, dmechFieldSourcedV, dmechFieldSourcedgradV);
+    dsourceVectorFielddB, dmechFieldSourcedT, dmechFieldSourcedgradT, dmechFieldSourcedV, dmechFieldSourcedgradV,
+    demFieldSourcedA, demFieldSourcedB, demFieldSourcedF, demFieldSourcedT, demFieldSourcedGradT, demFieldSourcedV, demFieldSourcedGradV);
 
   ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
 }
diff --git a/dG3D/src/nonLocalDamageDG3DIPVariable.h b/dG3D/src/nonLocalDamageDG3DIPVariable.h
index 882acd4672e57e37014b8b73bbce67daaacda7e5..0474d00d6e26c7e4b0d7fd10b58e0dd4af5f3ddb 100644
--- a/dG3D/src/nonLocalDamageDG3DIPVariable.h
+++ b/dG3D/src/nonLocalDamageDG3DIPVariable.h
@@ -967,6 +967,66 @@ virtual SVector3 & getRefToInductorSourceVectorField(const int idex)
         Msg::Error("getRefToEMFieldSource : _constitutiveCurlData not created");
     return _constitutiveCurlData->emFieldSource;
   }
+  virtual const std::vector< std::vector< SVector3 > > & getConstRefTodEMFieldSourcedVectorPotential() const
+  {
+      if (_constitutiveCurlData == NULL)
+          Msg::Error("getConstRefTodEMFieldSourcedVectorPotential : _constitutiveCurlData not created");
+      return _constitutiveCurlData->demFieldSourcedVectorPotential;
+  }
+  virtual std::vector< std::vector< SVector3 > > & getRefTodEMFieldSourcedVectorPotential()
+  {
+      if (_constitutiveCurlData == NULL)
+          Msg::Error("getRefTodEMFieldSourcedVectorPotential : _constitutiveCurlData not created");
+      return _constitutiveCurlData->demFieldSourcedVectorPotential;
+  }
+  virtual const std::vector< std::vector< SVector3 > > & getConstRefTodEMFieldSourcedVectorCurl() const
+  {
+      if (_constitutiveCurlData == NULL)
+          Msg::Error("getConstRefTodEMFieldSourcedVectorCurl : _constitutiveCurlData not created");
+      return _constitutiveCurlData->demFieldSourcedVectorCurl;
+  }
+  virtual std::vector< std::vector< SVector3 > > & getRefTodEMFieldSourcedVectorCurl()
+  {
+      if (_constitutiveCurlData == NULL)
+          Msg::Error("getRefTodEMFieldSourcedVectorCurl : _constitutiveCurlData not created");
+      return _constitutiveCurlData->demFieldSourcedVectorCurl;
+  }
+  virtual const std::vector< STensor3 > & getConstRefTodEMFieldSourcedF() const
+  {
+      if (_constitutiveCurlData == NULL)
+          Msg::Error("getConstRefTodEMFieldSourcedF : _constitutiveCurlData not created");
+      return _constitutiveCurlData->demFieldSourcedF;
+  }
+  virtual std::vector< STensor3 > & getRefTodEMFieldSourcedF()
+  {
+      if (_constitutiveCurlData == NULL)
+          Msg::Error("getRefTodEMFieldSourcedF : _constitutiveCurlData not created");
+      return _constitutiveCurlData->demFieldSourcedF;
+  }
+  virtual const fullMatrix<double> & getConstRefTodEMFieldSourcedField() const
+  {
+      if (_constitutiveCurlData == NULL)
+          Msg::Error("getConstRefTodEMFieldSourcedField : _constitutiveCurlData not created");
+      return _constitutiveCurlData->demFieldSourcedField;
+  }
+  virtual fullMatrix<double> & getRefTodEMFieldSourcedField()
+  {
+      if (_constitutiveCurlData == NULL)
+          Msg::Error("getRefTodEMFieldSourcedField : _constitutiveCurlData not created");
+      return _constitutiveCurlData->demFieldSourcedField;
+  }
+  virtual const std::vector< std::vector< SVector3 > > & getConstRefTodEMFieldSourcedGradField() const
+  {
+      if (_constitutiveCurlData == NULL)
+          Msg::Error("getConstRefTodEMFieldSourcedGradField : _constitutiveCurlData not created");
+      return _constitutiveCurlData->demFieldSourcedGradField;
+  }
+  virtual std::vector< std::vector< SVector3 > > & getRefTodEMFieldSourcedGradField()
+  {
+      if (_constitutiveCurlData == NULL)
+          Msg::Error("getConstRefTodEMFieldSourcedGradField : _constitutiveCurlData not created");
+      return _constitutiveCurlData->demFieldSourcedGradField;
+  }
 
   //Energy Conjugated Field (fT & fv)