diff --git a/NonLinearSolver/internalPoints/ipField.cpp b/NonLinearSolver/internalPoints/ipField.cpp
index 9ed837c93526b5331ba1374201704a37862da2f6..b96547a01867042d1f1f30f36f97a737909edbae 100644
--- a/NonLinearSolver/internalPoints/ipField.cpp
+++ b/NonLinearSolver/internalPoints/ipField.cpp
@@ -175,6 +175,7 @@ std::string IPField::ToString(const int i){
   else if (i == INDUCTORSOURCEVECTORFIELD_X) return "INDUCTORSOURCEVECTORFIELD_X";
   else if (i == INDUCTORSOURCEVECTORFIELD_Y) return "INDUCTORSOURCEVECTORFIELD_Y";
   else if (i == INDUCTORSOURCEVECTORFIELD_Z) return "INDUCTORSOURCEVECTORFIELD_Z";
+  else if (i == EMFIELDSOURCE) return "EMFIELDSOURCE";
   else if (i == LOCAL_POROSITY) return "LOCAL_POROSITY";
   else if (i == NONLOCAL_POROSITY) return "NONLOCAL_POROSITY";
   else if (i == CORRECTED_POROSITY) return "CORRECTED_POROSITY";
diff --git a/NonLinearSolver/internalPoints/ipField.h b/NonLinearSolver/internalPoints/ipField.h
index a2b493d0434effa08a6a98df2c4543258605bef7..df60f2db640eae12faca8e97f8ca3d9db9295f3f 100644
--- a/NonLinearSolver/internalPoints/ipField.h
+++ b/NonLinearSolver/internalPoints/ipField.h
@@ -138,7 +138,7 @@ class IPField : public elementsField {
                     MAGNETICINDUCTION_X, MAGNETICINDUCTION_Y, MAGNETICINDUCTION_Z,
                     MAGNETICFIELD_X, MAGNETICFIELD_Y, MAGNETICFIELD_Z,
                     EMSOURCEVECTORFIELD_X,EMSOURCEVECTORFIELD_Y,EMSOURCEVECTORFIELD_Z,
-                    INDUCTORSOURCEVECTORFIELD_X,INDUCTORSOURCEVECTORFIELD_Y,INDUCTORSOURCEVECTORFIELD_Z,
+                    INDUCTORSOURCEVECTORFIELD_X,INDUCTORSOURCEVECTORFIELD_Y,INDUCTORSOURCEVECTORFIELD_Z, EMFIELDSOURCE,
                     BROKEN};
     enum  Operator { MEAN_VALUE=1, MIN_VALUE, MAX_VALUE, CRUDE_VALUE};      
     static std::string ToString(const int i);
diff --git a/dG3D/src/FractureCohesiveDG3DIPVariable.h b/dG3D/src/FractureCohesiveDG3DIPVariable.h
index bc81eac4b98a5da0afe0c312fd7062f198a82f80..318e0c8b7dc13545613970da9dffef0b28ac6a3d 100644
--- a/dG3D/src/FractureCohesiveDG3DIPVariable.h
+++ b/dG3D/src/FractureCohesiveDG3DIPVariable.h
@@ -262,6 +262,8 @@ class FractureCohesive3DIPVariable : public IPVariable2ForFracture<dG3DIPVariabl
   virtual std::vector< std::vector< SVector3 > > & getRefTodMechanicalFieldSourcedNonLocalVariable() { return _ipvbulk->getRefTodMechanicalFieldSourcedNonLocalVariable();}
   virtual const SVector3 & getConstRefTodSourceVectorFielddt(const int idex) const {return _ipvbulk->getConstRefTodSourceVectorFielddt(idex);}
   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();}
 
   //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 af87419d8ad42f3df9be191683c8162aa1a5f4ba..09e267746fe023099abb5ee8e4c7bee61329dcb6 100644
--- a/dG3D/src/dG3DIPVariable.cpp
+++ b/dG3D/src/dG3DIPVariable.cpp
@@ -370,6 +370,7 @@ 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);
 
         std::vector< STensor3 > StdSTen3Tmp;
         StdSTen3Tmp.resize(numCurlVar, STensor3(0.0));
@@ -413,6 +414,7 @@ constitutiveCurlData::constitutiveCurlData(const constitutiveCurlData& src)
   inductorSourceVectorField(src.inductorSourceVectorField),
   sourceVectorField(src.sourceVectorField),
   dSourceVectorFielddt(src.dSourceVectorFielddt),
+  emFieldSource(src.emFieldSource),
   dSourceVectorFielddVectorCurl(src.dSourceVectorFielddVectorCurl),
   dSourceVectorFielddVectorPotential(src.dSourceVectorFielddVectorPotential),
   dVectorFielddVectorPotential(src.dVectorFielddVectorPotential),
@@ -447,6 +449,7 @@ constitutiveCurlData & constitutiveCurlData::operator=(const constitutiveCurlDat
     inductorSourceVectorField = src.inductorSourceVectorField;
     sourceVectorField = src.sourceVectorField;
     dSourceVectorFielddt = src.dSourceVectorFielddt;
+    emFieldSource = src.emFieldSource;
     dSourceVectorFielddVectorCurl = src.dSourceVectorFielddVectorCurl;
     dSourceVectorFielddVectorPotential = src.dSourceVectorFielddVectorPotential;
     dVectorFielddVectorPotential = src.dVectorFielddVectorPotential;
@@ -481,6 +484,7 @@ void constitutiveCurlData::restart()
     restartManager::restart(inductorSourceVectorField);
     restartManager::restart(sourceVectorField);
     restartManager::restart(dSourceVectorFielddt);
+    restartManager::restart(emFieldSource);
     restartManager::restart(dSourceVectorFielddVectorCurl);
     restartManager::restart(dSourceVectorFielddVectorPotential);
     restartManager::restart(dVectorFielddVectorPotential);
@@ -3602,6 +3606,10 @@ double ElecMagTherMechDG3DIPVariableBase::get(const int i) const
   {
       return getConstRefToinductorSourceVectorField()[2];
   }
+  else if(i == IPField::EMFIELDSOURCE)
+  {
+      return getConstRefToemFieldSource();
+  }
   else
   {
     return dG3DIPVariable::get(i);
diff --git a/dG3D/src/dG3DIPVariable.h b/dG3D/src/dG3DIPVariable.h
index 0c4bcdebe6b9aed0e47e6b0369e41a8009c7278a..ecc0933cb8ab96b77d10b7d5246e8ebedc70118f 100644
--- a/dG3D/src/dG3DIPVariable.h
+++ b/dG3D/src/dG3DIPVariable.h
@@ -310,6 +310,8 @@ class dG3DIPVariableBase : public ipFiniteStrain
   virtual std::vector< std::vector< SVector3 > > & getRefTodMechanicalFieldSourcedNonLocalVariable()=0;
   virtual const SVector3 & getConstRefTodSourceVectorFielddt(const int idex) const=0;
   virtual SVector3 & getRefTodSourceVectorFielddt(const int idex)=0;
+  virtual const fullVector<double> & getConstRefToEMFieldSource() const=0;
+  virtual fullVector<double> & getRefToEMFieldSource()=0;
   virtual const std::vector< std::vector< STensor3 > > & getConstRefTodSourceVectorFielddVectorCurl() const=0;
   virtual std::vector< std::vector< STensor3 > > & getRefTodSourceVectorFielddVectorCurl()=0;
   virtual const std::vector< std::vector< STensor3 > > & getConstRefTodSourceVectorFielddVectorPotential() const=0;
@@ -472,6 +474,7 @@ 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
   std::vector< std::vector< STensor3 > > dSourceVectorFielddVectorCurl;
   std::vector< std::vector< STensor3 > > dSourceVectorFielddVectorPotential;
 
@@ -1424,6 +1427,20 @@ protected:
     return _constitutiveCurlData->dSourceVectorFielddt[idex];
   }
 
+  virtual const fullVector<double> & getConstRefToEMFieldSource() const
+  {
+      if (_constitutiveCurlData == NULL)
+          Msg::Error("getConstRefToEMFieldSource(): _constitutiveCurlData not created");
+      return _constitutiveCurlData->emFieldSource;
+  }
+
+  virtual fullVector<double> & getRefToEMFieldSource()
+  {
+      if (_constitutiveCurlData == NULL)
+          Msg::Error("getRefToEMFieldSource(): _constitutiveCurlData not created");
+      return _constitutiveCurlData->emFieldSource;
+  }
+
     virtual const std::vector< std::vector< STensor3 > > & getConstRefTodSourceVectorFielddVectorPotential() const
     {
         if (_constitutiveCurlData == NULL)
@@ -3920,6 +3937,16 @@ class ElecMagTherMechDG3DIPVariableBase : public dG3DIPVariable
       return getRefTodSourceVectorFielddt(0);
     }
 
+    // EMFieldSource : w_AV
+    virtual const double getConstRefToemFieldSource() const
+    {
+        return getConstRefToEMFieldSource()(0);
+    }
+    virtual double & getRefToemFieldSource()
+    {
+        return getRefToEMFieldSource()(0);
+    }
+
     // dsourceVectorFielddA
     virtual const STensor3 & getConstRefTodSourceVectorFielddMagneticVectorPotential() const
     {
diff --git a/dG3D/src/nonLocalDamageDG3DIPVariable.h b/dG3D/src/nonLocalDamageDG3DIPVariable.h
index 373c32a40bec442669785bcd665c36281500094d..882acd4672e57e37014b8b73bbce67daaacda7e5 100644
--- a/dG3D/src/nonLocalDamageDG3DIPVariable.h
+++ b/dG3D/src/nonLocalDamageDG3DIPVariable.h
@@ -954,7 +954,19 @@ virtual SVector3 & getRefToInductorSourceVectorField(const int idex)
     if (_constitutiveCurlData == NULL)
       Msg::Error("getRefTodSourceVectorFielddt : _constitutiveCurlData not created");
     return _constitutiveCurlData->dSourceVectorFielddt[idex];
-  }  
+  }
+  virtual const fullVector<double> & getConstRefToEMFieldSource() const
+  {
+      if (_constitutiveCurlData == NULL)
+          Msg::Error("getConstRefToEMFieldSource : _constitutiveCurlData not created");
+      return _constitutiveCurlData->emFieldSource;
+  }
+  virtual fullVector<double> & getRefToEMFieldSource()
+  {
+    if (_constitutiveCurlData == NULL)
+        Msg::Error("getRefToEMFieldSource : _constitutiveCurlData not created");
+    return _constitutiveCurlData->emFieldSource;
+  }
 
   //Energy Conjugated Field (fT & fv)