diff --git a/dG3D/src/dG3DTerms.cpp b/dG3D/src/dG3DTerms.cpp
index d46e436d5449c58e670475726b7941318ab1dc60..a0f1d2f809104a463271f2cd69d05d5f82a36fe2 100644
--- a/dG3D/src/dG3DTerms.cpp
+++ b/dG3D/src/dG3DTerms.cpp
@@ -770,6 +770,11 @@ void dG3DForceBulk::get(MElement *ele,int npts,IntPt *GP,fullVector<double> &vFo
         if (getConstitutiveExtraDofDiffusionAccountMecaSource()){
            w += ipv->getConstRefToMechanicalSource()(extraDOFField);
         }
+        double w_AV = 0.0;
+        if (getConstitutiveCurlAccountFieldSource() && extraDOFField == 0)
+        {
+            w_AV += ipv->getConstRefToEMFieldSource()(extraDOFField);
+        }
         int indexField=3+getNumNonLocalVariable()+extraDOFField;
         int nbFFRow = nlspace->getNumShapeFunctions(ele,indexField);
         int nbFFTotalLastRow = nlspace->getShapeFunctionsIndex(ele,indexField);
@@ -782,6 +787,7 @@ void dG3DForceBulk::get(MElement *ele,int npts,IntPt *GP,fullVector<double> &vFo
           if(extraDOFField==0)
           { // corresponds to temperature, but when considering electro-mechanics, ???
             vFor(j+nbFFTotalLastRow) += getConstitutiveExtraDofDiffusionEqRatio()*ratio*(w*Vals[j+nbFFTotalLastRow]);// related to cp
+            vFor(j+nbFFTotalLastRow) += getConstitutiveExtraDofDiffusionEqRatio()*ratio*(w_AV*Vals[j+nbFFTotalLastRow]); // related to eddy current + hysteresis loss
           }
         }
       }
@@ -1510,6 +1516,21 @@ void dG3DStiffnessBulk::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double>
                   dwdV += ipv->getConstRefTodMechanicalSourcedField()(extraDOFField,extraDOFField2);
                 }
               }
+              SVector3 dw_AVdGradT(0.0);
+              SVector3 dw_AVdGradV(0.0);
+              double dw_AVdT = 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];
+                  }
+              }
               for(int j=0;j<nbFFRow; j++)
               {
                 for(int k=0;k<nbFFCol;k++)
@@ -1528,11 +1549,24 @@ void dG3DStiffnessBulk::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double>
                   {
                       // related to cp
                       mStiff(j+nbFFTotalLastRow,k+nbFFTotalLastCol)+=getConstitutiveExtraDofDiffusionEqRatio()*ratio*dwdt*Vals[j+nbFFTotalLastRow]*Vals[k+nbFFTotalLastCol];
+                      mStiff(j+nbFFTotalLastRow,k+nbFFTotalLastCol) -= getConstitutiveExtraDofDiffusionEqRatio()*ratio*dw_AVdT*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_AVdGradT.operator()(m)*Grads[k+nbFFTotalLastCol][m];
+                      }
                   }
                   if (extraDOFField == 0 && extraDOFField2 == 1)
                   {
                       mStiff(j+nbFFTotalLastRow,k+nbFFTotalLastCol)+=getConstitutiveExtraDofDiffusionEqRatio()*ratio*dwdV
                                                                         *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_AVdGradV.operator()(m)*Grads[k+nbFFTotalLastCol][m];
+                      }
                   }
                 }
               }
@@ -1683,6 +1717,8 @@ void dG3DStiffnessBulk::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double>
               const STensor3 & dFluxdVectorCurl = ipv->getConstRefTodFluxdVectorCurl()[extraDOFField][curlField];
 
               SVector3 dwTdVectorPotential(0.0);
+              SVector3 dw_AVdVectorPotential(0.0);
+              SVector3 dw_AVdVectorCurl(0.0);
 
               if (extraDOFField == 0)
               {
@@ -1695,6 +1731,12 @@ void dG3DStiffnessBulk::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double>
                   {
                       dwTdVectorPotential += ipv->getConstRefTodMechanicalSourcedVectorPotential()[extraDOFField][curlField];
                   }
+
+                  if (getConstitutiveCurlAccountFieldSource())
+                  {
+                      dw_AVdVectorPotential += ipv->getConstRefTodEMFieldSourcedVectorPotential()[extraDOFField][curlField];
+                      dw_AVdVectorCurl += ipv->getConstRefTodEMFieldSourcedVectorCurl()[extraDOFField][curlField];
+                  }
               }
 
               for (unsigned int j = 0; j < nbFFRow; ++j)
@@ -1723,6 +1765,16 @@ void dG3DStiffnessBulk::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double>
                                                                                            (Vals[j+nbFFTotalLastRow] *
                                                                                             dwTdVectorPotential.operator()(l) *
                                                                                             curlVals[k+numEdgeShapeFuncTotalLastCol][l]);
+                                  mStiff(j+nbFFTotalLastRow,
+                                         k+numEdgeShapeFuncTotalLastCol+numStandardDof) -= getConstitutiveExtraDofDiffusionEqRatio() * ratio *
+                                                                                           (Vals[j+nbFFTotalLastRow] *
+                                                                                            dw_AVdVectorPotential.operator()(l) *
+                                                                                            curlVals[k+numEdgeShapeFuncTotalLastCol][l]);
+                                  mStiff(j+nbFFTotalLastRow,
+                                         k+numEdgeShapeFuncTotalLastCol+numStandardDof) -= getConstitutiveExtraDofDiffusionEqRatio() * ratio *
+                                                                                           (Vals[j+nbFFTotalLastRow] *
+                                                                                            dw_AVdVectorCurl.operator()(l) *
+                                                                                            CurlcurlVals[k+numEdgeShapeFuncTotalLastCol][l]);
                               }
                           }
                       }