diff --git a/NonLinearSolver/internalPoints/ipField.cpp b/NonLinearSolver/internalPoints/ipField.cpp
index bbbe9e24b9763e827506ec81c5ea158c46b246a5..2926a01bd5aec493934df7e26acaf16af9c6d49b 100644
--- a/NonLinearSolver/internalPoints/ipField.cpp
+++ b/NonLinearSolver/internalPoints/ipField.cpp
@@ -691,10 +691,12 @@ std::string IPField::ToString(const int i){
   else if (i == THERMALFLUX_X) return "THERMALFLUX_X";
   else if (i == THERMALFLUX_Y) return "THERMALFLUX_Y";
   else if (i == THERMALFLUX_Z) return "THERMALFLUX_Z";
+  else if (i == THERMALSOURCE) return "THERMALSOURCE";
   else if (i == VOLTAGE) return "VOLTAGE";
   else if (i == ELECTRICALFLUX_X) return "ELECTRICALFLUX_X";
   else if (i == ELECTRICALFLUX_Y) return "ELECTRICALFLUX_Y";
   else if (i == ELECTRICALFLUX_Z) return "ELECTRICALFLUX_Z";
+  else if (i == ELECTRICALSOURCE) return "ELECTRICALSOURCE";
   else if (i == MAGNETICVECTORPOTENTIAL_X) return "MAGNETICVECTORPOTENTIAL_X";
   else if (i == MAGNETICVECTORPOTENTIAL_Y) return "MAGNETICVECTORPOTENTIAL_Y";
   else if (i == MAGNETICVECTORPOTENTIAL_Z) return "MAGNETICVECTORPOTENTIAL_Z";
diff --git a/NonLinearSolver/internalPoints/ipField.h b/NonLinearSolver/internalPoints/ipField.h
index 1c6bdb3b7a97bca82976ae75363f54807d21185f..55ec5a4d416a01ea67f67fc502d98553aa11fa43 100644
--- a/NonLinearSolver/internalPoints/ipField.h
+++ b/NonLinearSolver/internalPoints/ipField.h
@@ -261,7 +261,7 @@ class IPField : public elementsField {
                     MTX_SVM,MTX_SIG_XX,MTX_SIG_YY,MTX_SIG_ZZ,MTX_SIG_XY,MTX_SIG_YZ,MTX_SIG_XZ,
                     INC_SVM,INC_SIG_XX,INC_SIG_YY,INC_SIG_ZZ,INC_SIG_XY,INC_SIG_YZ,INC_SIG_XZ,
                     TEMPERATURE,THERMALFLUX_X,THERMALFLUX_Y,THERMALFLUX_Z,THERMALSOURCE,
-                    VOLTAGE,ELECTRICALFLUX_X,ELECTRICALFLUX_Y,ELECTRICALFLUX_Z,
+                    VOLTAGE,ELECTRICALFLUX_X,ELECTRICALFLUX_Y,ELECTRICALFLUX_Z,ELECTRICALSOURCE,
                     LOCAL_POROSITY, NONLOCAL_POROSITY, CORRECTED_POROSITY, YIELD_POROSITY,
                     NUCLEATED_POROSITY_TOT, NUCLEATED_POROSITY_0, NUCLEATED_POROSITY_1, NUCLEATED_POROSITY_2,
                     NUCLEATION_ONSET_0, NUCLEATION_ONSET_1, NUCLEATION_ONSET_2,
diff --git a/NonLinearSolver/materialLaw/mlawElecMagGenericThermoMech.cpp b/NonLinearSolver/materialLaw/mlawElecMagGenericThermoMech.cpp
index 7e6878351c985d619d77b93b6fa41cce50bf1ba1..284fc7abda36abd53686be3bfb318b22dd29db50 100644
--- a/NonLinearSolver/materialLaw/mlawElecMagGenericThermoMech.cpp
+++ b/NonLinearSolver/materialLaw/mlawElecMagGenericThermoMech.cpp
@@ -262,6 +262,10 @@ void mlawElecMagGenericThermoMech::constitutive(
         double           &dwthdt,
         double           &dwthdV,
         STensor3         &dwthdF,
+        double           &wV,  //related to dV/dt
+        double           &dwVdt,
+        double           &dwVdV,
+        STensor3         &dwVdF,
         STensor3 &l10,
         STensor3 &l20,
         STensor3 &k10,
@@ -387,6 +391,13 @@ void mlawElecMagGenericThermoMech::constitutive(
     STensorOperation::zero(invP_TM);
     STensorOperation::inverseSTensor3(P_TM,invP_TM);
 
+    wV=0.;  //related to dV/dt
+    dwVdt=0.;
+    dwVdV=0.;
+    STensorOperation::zero(dwVdF);
+
+
+
     // if true then we solve EM problem and thus reevaluate Voltage field as well
     // else do not consider Voltage contributions -> TM problem
     if(_evaluateCurlField)
@@ -1397,4 +1408,4 @@ void mlawElecMagGenericThermoMech::constitutive(
             }
         }
     }
-}
\ No newline at end of file
+}
diff --git a/NonLinearSolver/materialLaw/mlawElecMagGenericThermoMech.h b/NonLinearSolver/materialLaw/mlawElecMagGenericThermoMech.h
index 7af4ba0d03695e532790e07af286e177266605ce..3c0bb566b3f6c7d049f252784232e5cfa49adead 100644
--- a/NonLinearSolver/materialLaw/mlawElecMagGenericThermoMech.h
+++ b/NonLinearSolver/materialLaw/mlawElecMagGenericThermoMech.h
@@ -31,6 +31,7 @@ protected:
     double _lz;
     // Initial electric scalar potential
     double _v0;
+    
     // Magnetic permeability directional components
     double _mu_x;
     double _mu_y;
@@ -263,6 +264,10 @@ public:
             double           &dwthdt,
             double           &dwthdV,
             STensor3         &dwthdF,
+            double           &wV,  //related to dV/dt
+            double           &dwVdt,
+            double           &dwVdV,
+            STensor3         &dwVdF,
             STensor3 &l10,
             STensor3 &l20,
             STensor3 &k10,
diff --git a/NonLinearSolver/materialLaw/mlawLinearElecMagTherMech.cpp b/NonLinearSolver/materialLaw/mlawLinearElecMagTherMech.cpp
index d24596f66cbd4e933dbafd5bcce2ee37f0758989..9299b5eba885fd0cc30dc3d01244f001f7fa88c8 100644
--- a/NonLinearSolver/materialLaw/mlawLinearElecMagTherMech.cpp
+++ b/NonLinearSolver/materialLaw/mlawLinearElecMagTherMech.cpp
@@ -157,6 +157,10 @@ void mlawLinearElecMagTherMech::constitutive(
                                             double &dw_Tdt,
                                             double &dw_TdV,
                                             STensor3 &dw_Tdf,
+                                            double &wV,  //related to dV/dt
+                                            double &dwVdt,
+                                            double &dwVdV,
+                                            STensor3 &dwVdF,
                                             STensor3 &l10,
                                             STensor3 &l20,
                                             STensor3 &k10,
@@ -313,6 +317,13 @@ void mlawLinearElecMagTherMech::constitutive(
     STensorOperation::zero(dVoltageEMFieldSourcedV);
     STensorOperation::zero(dVoltageEMFieldSourcedGradV);
 
+    wV=0.;  //related to dV/dt
+    dwVdt=0.;
+    dwVdV=0.;
+    STensorOperation::zero(dwVdF);
+
+
+
     static STensor3 dldT,dkdT,dnudT;
     double dseebeckdT = 0.0;
     STensor33 d_nu0dB(0.0);
diff --git a/NonLinearSolver/materialLaw/mlawLinearElecMagTherMech.h b/NonLinearSolver/materialLaw/mlawLinearElecMagTherMech.h
index 98697e2a821949a1e4d7544c2c8a368b4271ae7b..89aa2f3c74bf271300655094bc5d4b068f3086a6 100644
--- a/NonLinearSolver/materialLaw/mlawLinearElecMagTherMech.h
+++ b/NonLinearSolver/materialLaw/mlawLinearElecMagTherMech.h
@@ -109,6 +109,10 @@ class mlawLinearElecMagTherMech : public mlawLinearElecTherMech
                             double &dw_Tdt,
                             double &dw_TdV,
                             STensor3 &dw_Tdf,
+                            double &wV,  //related to dV/dt
+                            double &dwVdt,
+                            double &dwVdV,
+                            STensor3 &dwVdF,
                             STensor3 &l10,
                             STensor3 &l20,
                             STensor3 &k10,
diff --git a/dG3D/src/dG3DIPVariable.cpp b/dG3D/src/dG3DIPVariable.cpp
index b616ee77ed1fc0e0215627ab51c3881e67fbc603..a35210e82d7494d1574175dc3f54e699dc3829d6 100644
--- a/dG3D/src/dG3DIPVariable.cpp
+++ b/dG3D/src/dG3DIPVariable.cpp
@@ -5059,6 +5059,10 @@ double ElecMagTherMechDG3DIPVariableBase::get(const int i) const
   {
       return getConstRefToThermalSource();
   }
+  else if(i == IPField::ELECTRICALSOURCE)
+  {
+      return getConstRefToVoltageSource();
+  }
   else if(i == IPField::VOLTAGE)
   {
     return getConstRefToVoltage();
@@ -5179,6 +5183,10 @@ double & ElecMagTherMechDG3DIPVariableBase::getRef(const int i)
     {
         return getRefToThermalSource();
     }
+    else if(i == IPField::ELECTRICALSOURCE)
+    {
+        return getRefToVoltageSource();
+    }
     else if(i == IPField::VOLTAGE)
     {
         return getRefToVoltage();
diff --git a/dG3D/src/dG3DIPVariable.h b/dG3D/src/dG3DIPVariable.h
index ebb74105c74d4d6a8459d9653d74dcc5d8e9da1d..684b5cdb5304fb9b4f62c8c27455da3bb36db130 100644
--- a/dG3D/src/dG3DIPVariable.h
+++ b/dG3D/src/dG3DIPVariable.h
@@ -4159,6 +4159,9 @@ class ElecTherMechDG3DIPVariableBase : public dG3DIPVariable
   virtual const double getConstRefTodThermalSourcedT() const {return getConstRefTodFieldSourcedField()(0,0);}
   virtual double &getRefTodThermalSourcedT() {return  getRefTodFieldSourcedField()(0,0);}
 
+  virtual const double getConstRefTodThermalSourcedV() const {return getConstRefTodFieldSourcedField()(0,1);}
+  virtual double &getRefTodThermalSourcedV() {return  getRefTodFieldSourcedField()(0,1);}
+
   virtual const  STensor3 &getConstRefTodThermalSourcedF() const {return getConstRefTodFieldSourcedF()[0];}
   virtual  STensor3 &getRefTodThermalSourcedF() {return getRefTodFieldSourcedF()[0];}
 
@@ -4168,6 +4171,18 @@ class ElecTherMechDG3DIPVariableBase : public dG3DIPVariable
   virtual const  SVector3 &getConstRefTodThermalSourcedGradT() const {return  dThermalSourcedGradT;}//dwdgradT
   virtual  SVector3 &getRefTodThermalSourcedGradT() {return  dThermalSourcedGradT;}*/
 
+  virtual const double getConstRefToVoltageSource() const {return getConstRefToFieldSource()(1);}
+  virtual double &getRefToVoltageSource() {return  getRefToFieldSource()(1);}
+
+  virtual const double getConstRefTodVoltageSourcedT() const {return getConstRefTodFieldSourcedField()(1,0);}
+  virtual double &getRefTodVoltageSourcedT() {return  getRefTodFieldSourcedField()(1,0);}
+
+  virtual const double getConstRefTodVoltageSourcedV() const {return getConstRefTodFieldSourcedField()(1,1);}
+  virtual double &getRefTodVoltageSourcedV() {return  getRefTodFieldSourcedField()(1,1);}
+
+  virtual const  STensor3 &getConstRefTodVoltageSourcedF() const {return getConstRefTodFieldSourcedF()[1];}
+  virtual  STensor3 &getRefTodVoltageSourcedF() {return getRefTodFieldSourcedF()[1];}
+
   virtual const double getConstRefTofT() const {return  getConstRefToEnergyConjugatedField(0);}
   virtual double &getRefTofT() {return getRefToEnergyConjugatedField(0);}
 
@@ -4608,6 +4623,20 @@ class ElecMagTherMechDG3DIPVariableBase : public dG3DIPVariable
     virtual const  STensor3 &getConstRefTodThermalSourcedF() const {return getConstRefTodFieldSourcedF()[0];}
     virtual  STensor3 &getRefTodThermalSourcedF() {return getRefTodFieldSourcedF()[0];}
 
+    // VoltageSource
+    virtual const double getConstRefToVoltageSource() const {return getConstRefToFieldSource()(1);}
+    virtual double &getRefToVoltageSource() {return  getRefToFieldSource()(1);}
+
+    virtual const double getConstRefTodVoltageSourcedT() const {return getConstRefTodFieldSourcedField()(1,0);}
+    virtual double &getRefTodVoltageSourcedT() {return  getRefTodFieldSourcedField()(1,0);}
+
+    virtual const double getConstRefTodVoltageSourcedV() const {return getConstRefTodFieldSourcedField()(1,1);}
+    virtual double &getRefTodVoltageSourcedV() {return  getRefTodFieldSourcedField()(1,1);}
+
+    virtual const  STensor3 &getConstRefTodVoltageSourcedF() const {return getConstRefTodFieldSourcedF()[1];}
+    virtual  STensor3 &getRefTodVoltageSourcedF() {return getRefTodFieldSourcedF()[1];}
+
+
     // fT
     virtual const double getConstRefTofT() const {return  getConstRefToEnergyConjugatedField(0);}
     virtual double &getRefTofT() {return getRefToEnergyConjugatedField(0);}
diff --git a/dG3D/src/dG3DMaterialLaw.cpp b/dG3D/src/dG3DMaterialLaw.cpp
index 029a582049e60db9ac4edb919164f03f5368e327..caa48bd8cbac2fb4ee188aaddd17c963815c8bc0 100644
--- a/dG3D/src/dG3DMaterialLaw.cpp
+++ b/dG3D/src/dG3DMaterialLaw.cpp
@@ -8638,11 +8638,17 @@ const bool checkfrac, const bool dTangent)
     STensor33& dedF=ipvcur->getRefTodFluxjedF();
     STensor3& P=ipvcur->getRefToFirstPiolaKirchhoffStress();
     STensor43& Tangent=ipvcur->getRefToTangentModuli();
+    
     double &w_T = ipvcur->getRefToThermalSource();
     double &dw_TdT = ipvcur->getRefTodThermalSourcedT();
     double &dw_TdV = ipvcur->getRefTodThermalSourcedV();
     STensor3& dw_TdF = ipvcur->getRefTodThermalSourcedF();
 
+    double &w_V = ipvcur->getRefToThermalSource();
+    double &dw_VdT = ipvcur->getRefTodThermalSourcedT();
+    double &dw_VdV = ipvcur->getRefTodThermalSourcedV();
+    STensor3& dw_VdF = ipvcur->getRefTodThermalSourcedF();
+
     ipvcur->setRefToTherConductivity(this->lineark10);
     ipvcur->setRefToCrossTherConductivity(this->lineark20);
     ipvcur->setRefToElecConductivity(this->linearl10);
@@ -8783,7 +8789,7 @@ const bool checkfrac, const bool dTangent)
 
     _lawLinearEMTM->constitutive(F0,Fn,P, q0, q1, Tangent,T0, temperature, voltage,gradT, gradV,
     *fluxT,fluxje, dPdT,dPdV, *dqdgradT,*dqdT, *dqdgradV,*dqdV,dedV,dedgradV,dedT,dedgradT,*dqdF,dedF,
-    w_T,dw_TdT,dw_TdV,dw_TdF,linearl10,linearl20,lineark10,lineark20,dl10dT,dl20dT,dk10dT,dk20dT,dl10dv,
+    w_T,dw_TdT,dw_TdV,dw_TdF,w_V,dw_VdT,dw_VdV,dw_VdF,linearl10,linearl20,lineark10,lineark20,dl10dT,dl20dT,dk10dT,dk20dT,dl10dv,
     dl20dv,dk10dv,dk20dv,dl10dF,dl20dF,dk10dF,dk20dF,*fluxjy,*djydV,*djydgradV,*djydT,*djydgradT,
     *djydF,*djydA, *djydB, *djydgradVdT,*djydgradVdV,*djydgradVdF,*djydgradTdT,*djydgradTdV,*djydgradTdF,stiff,
     A0, An, B, H, js0, dBdT, dBdGradT, dBdV, dBdGradV, dBdF, dBdA, *dfluxTdA, dfluxjedA, *dfluxTdB, dfluxjedB, dw_TdA, dw_TdB, dmechSourcedB, sourceVectorField,
@@ -8951,6 +8957,11 @@ const bool checkfrac, const bool dTangent)
   double &dw_TdV = ipvcur->getRefTodThermalSourcedV();
   STensor3& dw_TdF = ipvcur->getRefTodThermalSourcedF();
 
+  double &w_V = ipvcur->getRefToThermalSource();
+  double &dw_VdT = ipvcur->getRefTodThermalSourcedT();
+  double &dw_VdV = ipvcur->getRefTodThermalSourcedV();
+  STensor3& dw_VdF = ipvcur->getRefTodThermalSourcedF();
+
   ipvcur->setRefToTherConductivity(this->lineark10);
   ipvcur->setRefToCrossTherConductivity(this->lineark20);
   ipvcur->setRefToElecConductivity(this->linearl10);
@@ -9091,7 +9102,7 @@ const bool checkfrac, const bool dTangent)
 
   _lawLinearEMInductor->constitutive(F0,Fn,P, q0, q1, Tangent,T0, temperature, voltage,gradT, gradV,
     *fluxT,fluxje, dPdT,dPdV, *dqdgradT,*dqdT, *dqdgradV,*dqdV,dedV,dedgradV,dedT,dedgradT,*dqdF,dedF,
-    w_T,dw_TdT,dw_TdV,dw_TdF,linearl10,linearl20,lineark10,lineark20,dl10dT,dl20dT,dk10dT,dk20dT,dl10dv,
+    w_T,dw_TdT,dw_TdV,dw_TdF,w_V,dw_VdT,dw_VdV,dw_VdF,linearl10,linearl20,lineark10,lineark20,dl10dT,dl20dT,dk10dT,dk20dT,dl10dv,
     dl20dv,dk10dv,dk20dv,dl10dF,dl20dF,dk10dF,dk20dF,*fluxjy,*djydV,*djydgradV,*djydT,*djydgradT,
     *djydF, *djydA, *djydB, *djydgradVdT,*djydgradVdV,*djydgradVdF,*djydgradTdT,*djydgradTdV,*djydgradTdF,stiff,
     A0, An, B, H, js0, dBdT, dBdGradT, dBdV, dBdGradV, dBdF, dBdA, *dfluxTdA, dfluxjedA, *dfluxTdB, dfluxjedB, dw_TdA, dw_TdB, dmechSourcedB,
@@ -9854,6 +9865,11 @@ void ElecMagGenericThermoMechanicsDG3DMaterialLaw::stress(IPVariable*ipv, const
     double &dw_TdV = ipvcur->getRefTodThermalSourcedV();
     STensor3& dw_TdF = ipvcur->getRefTodThermalSourcedF();
 
+    double &w_V = ipvcur->getRefToVoltageSource();
+    double &dw_VdT = ipvcur->getRefTodVoltageSourcedT();
+    double &dw_VdV = ipvcur->getRefTodVoltageSourcedV();
+    STensor3& dw_VdF = ipvcur->getRefTodVoltageSourcedF();
+
     double&  mechSource =  ipvcur->getRefToMechanicalSource()(0);
     double& dmechSourcedt = ipvcur->getRefTodMechanicalSourcedField()(0,0);
     double& dmechSourcedV = ipvcur->getRefTodMechanicalSourcedField()(0,1);
@@ -10003,7 +10019,7 @@ void ElecMagGenericThermoMechanicsDG3DMaterialLaw::stress(IPVariable*ipv, const
                                                dPdT,dPdV, dPdE, *dqdgradT,*dqdT, *dqdgradV,*dqdV,*dqdF,
                                                dedgradT,dedT,dedgradV,dedV,dedF,
                                                dDdgradT,dDdT,dDdgradV,dDdV,dDdF,
-                                               stiff, w_T,dw_TdT,dw_TdV,dw_TdF,linearl10,linearl20,
+                                               stiff, w_T,dw_TdT,dw_TdV,dw_TdF,w_V,dw_VdT,dw_VdV,dw_VdF,linearl10,linearl20,
                                                lineark10,lineark20,
                                                dl10dT,dl20dT,dk10dT,dk20dT,dl10dv,
                                                dl20dv,dk10dv,dk20dv,dl10dF,dl20dF,dk10dF,dk20dF,
@@ -10247,11 +10263,17 @@ void ElecMagInductorDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipv
     STensor3& dPdV=ipvcur->getRefTodPdV();
     STensor33& dPdE=ipvcur->getRefTodPdGradV();
     STensor43& Tangent=ipvcur->getRefToTangentModuli();
+
     double &w_T = ipvcur->getRefToThermalSource();
     double &dw_TdT = ipvcur->getRefTodThermalSourcedT();
     double &dw_TdV = ipvcur->getRefTodThermalSourcedV();
     STensor3& dw_TdF = ipvcur->getRefTodThermalSourcedF();
 
+    double &w_V = ipvcur->getRefToVoltageSource();
+    double &dw_VdT = ipvcur->getRefTodVoltageSourcedT();
+    double &dw_VdV = ipvcur->getRefTodVoltageSourcedV();
+    STensor3& dw_VdF = ipvcur->getRefTodVoltageSourcedF();
+
     double&  mechSource =  ipvcur->getRefToMechanicalSource()(0);
     double& dmechSourcedt = ipvcur->getRefTodMechanicalSourcedField()(0,0);
     double& dmechSourcedV = ipvcur->getRefTodMechanicalSourcedField()(0,1);
@@ -10403,7 +10425,7 @@ void ElecMagInductorDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipv
                                       dPdT,dPdV, dPdE, *dqdgradT,*dqdT, *dqdgradV,*dqdV,*dqdF,
                                       dedgradT,dedT,dedgradV,dedV,dedF,
                                       dDdgradT,dDdT,dDdgradV,dDdV,dDdF,
-                                      stiff, w_T,dw_TdT,dw_TdV,dw_TdF,linearl10,linearl20,
+                                      stiff, w_T,dw_TdT,dw_TdV,dw_TdF,w_V,dw_VdT,dw_VdV,dw_VdF,linearl10,linearl20,
                                       lineark10,lineark20,
                                       dl10dT,dl20dT,dk10dT,dk20dT,dl10dv,
                                       dl20dv,dk10dv,dk20dv,dl10dF,dl20dF,dk10dF,dk20dF,
diff --git a/dG3D/src/dG3DTerms.cpp b/dG3D/src/dG3DTerms.cpp
index aa79bc80a2235fe96a35bee7f7d9fea727a67862..e9c9b96025e4cfb58de6505d8652f1c4fba4c326 100644
--- a/dG3D/src/dG3DTerms.cpp
+++ b/dG3D/src/dG3DTerms.cpp
@@ -794,7 +794,7 @@ void dG3DForceBulk::get(MElement *ele,int npts,IntPt *GP,fullVector<double> &vFo
           {
             vFor(j+nbFFTotalLastRow) += getConstitutiveExtraDofDiffusionEqRatio()*ratio*(Grads[j+nbFFTotalLastRow][m]*flux(m));
           }
-          if(extraDOFField==0)
+          //if(extraDOFField==0)
           { // corresponds to temperature, but when considering electro-mechanics, ???
             vFor(j+nbFFTotalLastRow) += getConstitutiveExtraDofDiffusionEqRatio()*ratio*(w*Vals[j+nbFFTotalLastRow]);// related to cp
           }
@@ -1470,15 +1470,15 @@ void dG3DStiffnessBulk::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double>
               double dwdt = 0.;
               if (getConstitutiveExtraDofDiffusionAccountFieldSource())
               {
-                if(extraDOFField==0 && extraDOFField2==0){
+                //if(extraDOFField==0 && extraDOFField2==0){
                   dwdt += ipv->getConstRefTodFieldSourcedField()(extraDOFField,extraDOFField2);
-                }
+                //}
               }
               if (getConstitutiveExtraDofDiffusionAccountMecaSource())
               {
-                if(extraDOFField==0 && extraDOFField2==0){
+                //if(extraDOFField==0 && extraDOFField2==0){
                   dwdt += ipv->getConstRefTodMechanicalSourcedField()(extraDOFField,extraDOFField2);
-                }
+                //}
               }
               for(int j=0;j<nbFFRow; j++)
               {
@@ -1505,8 +1505,8 @@ void dG3DStiffnessBulk::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double>
                       }
                     }
                   }
-                  if(extraDOFField==0 && extraDOFField2==0) 
-                    mStiff(j+nbFFTotalLastRow,k+nbFFTotalLastCol)+= getConstitutiveExtraDofDiffusionEqRatio()*ratio*dwdt*dFielddEnergyConjugatedField1*
+                  //if(extraDOFField==0 && extraDOFField2==0) 
+                  mStiff(j+nbFFTotalLastRow,k+nbFFTotalLastCol)+= getConstitutiveExtraDofDiffusionEqRatio()*ratio*dwdt*dFielddEnergyConjugatedField1*
                                 Vals[j+nbFFTotalLastRow]*Vals[k+nbFFTotalLastCol];// related to cp
                 }
               }
@@ -1530,44 +1530,44 @@ void dG3DStiffnessBulk::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double>
               const STensor3  &dFluxdGradField = ipv->getConstRefTodFluxdGradField()[extraDOFField][extraDOFField2]; //in reference conf
               const SVector3  &dFluxdField     = ipv->getConstRefTodFluxdField()[extraDOFField][extraDOFField2]; //in reference conf
 
-              double dwdt = 0.;
-              double dwdV = 0.0;
+              double dwdField = 0.;
+              //double dwdV = 0.0;
               if (getConstitutiveExtraDofDiffusionAccountFieldSource())
               {
-                if(extraDOFField==0 && extraDOFField2==0)
+                //if(extraDOFField==0 && extraDOFField2==0)
                 {
-                  dwdt += ipv->getConstRefTodFieldSourcedField()(extraDOFField,extraDOFField2);
-                }
-                if (extraDOFField == 0 && extraDOFField2 == 1)
-                {
-                    dwdV += ipv->getConstRefTodFieldSourcedField()(extraDOFField,extraDOFField2);
+                  dwdField += ipv->getConstRefTodFieldSourcedField()(extraDOFField,extraDOFField2);
                 }
+                //if (extraDOFField == 0 && extraDOFField2 == 1)
+                //{
+                //    dwdV += ipv->getConstRefTodFieldSourcedField()(extraDOFField,extraDOFField2);
+                //}
               }
               if (getConstitutiveExtraDofDiffusionAccountMecaSource())
               {
-                if(extraDOFField==0 && extraDOFField2==0)
-                {
-                  dwdt += ipv->getConstRefTodMechanicalSourcedField()(extraDOFField,extraDOFField2);
-                }
-                if (extraDOFField == 0 && extraDOFField2 == 1)
-                {
-                  dwdV += ipv->getConstRefTodMechanicalSourcedField()(extraDOFField,extraDOFField2);
-                }
+                //if(extraDOFField==0 && extraDOFField2==0)
+                //{
+                  dwdField += ipv->getConstRefTodMechanicalSourcedField()(extraDOFField,extraDOFField2);
+                //}
+                //if (extraDOFField == 0 && extraDOFField2 == 1)
+                //{
+                //  dwdV += ipv->getConstRefTodMechanicalSourcedField()(extraDOFField,extraDOFField2);
+                //}
               }
-              SVector3 dw_AVdGradT(0.0);
-              SVector3 dw_AVdGradV(0.0);
-              double dw_AVdT = 0.0;
+              SVector3 dw_AVdGradField(0.0);
+              //SVector3 dw_AVdGradV(0.0);
+              double dw_AVdField = 0.0;
               if (getConstitutiveCurlAccountFieldSource())
               {
-                  if (extraDOFField == 0 && extraDOFField2 == 0)
-                  {
-                      dw_AVdT += ipv->getConstRefTodEMFieldSourcedField()(extraDOFField,extraDOFField2);
-                      dw_AVdGradT += ipv->getConstRefTodEMFieldSourcedGradField()[extraDOFField][extraDOFField2];
-                  }
-                  if (extraDOFField == 0 && extraDOFField2 == 1)
-                  {
-                      dw_AVdGradV += ipv->getConstRefTodEMFieldSourcedGradField()[extraDOFField][extraDOFField2];
-                  }
+                  //if (extraDOFField == 0 && extraDOFField2 == 0)
+                  //{
+                      dw_AVdField += ipv->getConstRefTodEMFieldSourcedField()(extraDOFField,extraDOFField2);
+                      dw_AVdGradField += ipv->getConstRefTodEMFieldSourcedGradField()[extraDOFField][extraDOFField2];
+                  //}
+                  //if (extraDOFField == 0 && extraDOFField2 == 1)
+                  //{
+                  //    dw_AVdGradV += ipv->getConstRefTodEMFieldSourcedGradField()[extraDOFField][extraDOFField2];
+                  //}
               }
               for(int j=0;j<nbFFRow; j++)
               {
@@ -1583,6 +1583,14 @@ void dG3DStiffnessBulk::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double>
                     }
                   }
                
+                  mStiff(j+nbFFTotalLastRow,k+nbFFTotalLastCol)+=getConstitutiveExtraDofDiffusionEqRatio()*ratio*dwdField*Vals[j+nbFFTotalLastRow]*Vals[k+nbFFTotalLastCol];
+                  mStiff(j+nbFFTotalLastRow,k+nbFFTotalLastCol) += getConstitutiveExtraDofDiffusionEqRatio()*ratio*dw_AVdField*Vals[j+nbFFTotalLastRow]*Vals[k+nbFFTotalLastCol];
+                  for (unsigned int m = 0; m < 3; ++m)
+                  {
+                     mStiff(j+nbFFTotalLastRow,k+nbFFTotalLastCol) += getConstitutiveExtraDofDiffusionEqRatio()*ratio*Vals[j+nbFFTotalLastRow]*
+                                                                           dw_AVdGradField.operator()(m)*Grads[k+nbFFTotalLastCol][m];
+                  }
+                  /*
                   if(extraDOFField==0 && extraDOFField2==0)
                   {
                       // related to cp
@@ -1605,7 +1613,7 @@ void dG3DStiffnessBulk::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double>
                           mStiff(j+nbFFTotalLastRow,k+nbFFTotalLastCol) += getConstitutiveExtraDofDiffusionEqRatio()*ratio*Vals[j+nbFFTotalLastRow]*
                                                                            dw_AVdGradV.operator()(m)*Grads[k+nbFFTotalLastCol][m];
                       }
-                  }
+                  }*/
                 }
               }
             }
@@ -1626,21 +1634,21 @@ void dG3DStiffnessBulk::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double>
 
           if (getConstitutiveExtraDofDiffusionAccountFieldSource())
           {
-            if(extraDOFField==0)
-            {
+            //if(extraDOFField==0)
+            //{
               dwdf += ipv->getConstRefTodFieldSourcedF()[extraDOFField];
-            }
+            //}
           }
           if (getConstitutiveExtraDofDiffusionAccountMecaSource())
           {
-            if(extraDOFField==0)
-            {
+            //if(extraDOFField==0)
+            //{
               dwdf += ipv->getConstRefTodMechanicalSourcedF()[extraDOFField];
-            }
+            //}
           }
           if (getConstitutiveCurlAccountFieldSource())
           {
-              if (extraDOFField==0)
+              //if (extraDOFField==0)
                   dwdf += ipv->getConstRefTodEMFieldSourcedF()[extraDOFField]; // related to EM heat source from eddy + hysteresis loss
           }
           
@@ -1765,7 +1773,7 @@ void dG3DStiffnessBulk::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double>
               SVector3 dw_AVdVectorPotential(0.0);
               SVector3 dw_AVdVectorCurl(0.0);
 
-              if (extraDOFField == 0)
+              //if (extraDOFField == 0)
               {
                   if (getConstitutiveExtraDofDiffusionAccountFieldSource())
                   {
@@ -1801,7 +1809,7 @@ void dG3DStiffnessBulk::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double>
                                                                                         CurlcurlVals[k+numEdgeShapeFuncTotalLastCol][n]);
                           }
 
-                          if (extraDOFField == 0)
+                          //if (extraDOFField == 0)
                           {
                               for (unsigned int l = 0; l < 3; ++l)
                               {