diff --git a/dG3D/src/dG3DTerms.cpp b/dG3D/src/dG3DTerms.cpp
index 925809c2efb0625f53b22119170c80096b25d2d3..6ef5b56e333f212d51b042be8bffe6281b718935 100644
--- a/dG3D/src/dG3DTerms.cpp
+++ b/dG3D/src/dG3DTerms.cpp
@@ -817,6 +817,7 @@ void dG3DForceBulk::get(MElement *ele,int npts,IntPt *GP,fullVector<double> &vFo
         const int nbFFRow =  nlspace->getNumShapeFunctions(ele,curveFieldIndex); // number of shape functinon for curveFieldIndex in curlVals and CurlcurlVals
         const int nbFFTotalLastRow = nlspace->getShapeFunctionsIndex(ele,curveFieldIndex); // the position of the first element of shape functinon for curveFieldIndex in curlVals and CurlcurlVals        
 
+        // SourceVectorField = -fluxD
         SVector3 w(0.0);
         if (getConstitutiveCurlAccountFieldSource())
         {
@@ -845,7 +846,7 @@ void dG3DForceBulk::get(MElement *ele,int npts,IntPt *GP,fullVector<double> &vFo
       }
     }
   }
-  vFor.print("force");
+  //vFor.print("force");
 }
 
 void dG3DStiffnessBulk::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &mStiff) const
@@ -1729,7 +1730,6 @@ void dG3DStiffnessBulk::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double>
             }
           }            
         }
-          std::cout << mStiff.operator()(32,25) << std::endl;
       }
     }
     
@@ -1776,21 +1776,15 @@ void dG3DStiffnessBulk::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double>
 
                 const STensor3 &dHdA = ipv->getConstRefTodVectorFielddVectorPotential()[curlField][curlField2];
                 const STensor3 &dHdB = ipv->getConstRefTodVectorFielddVectorCurl()[curlField][curlField2];
-                // only dFluxDdA, dFluxDdB contributions
-                const STensor3 &dFluxdVectorCurl = ipv->getConstRefTodFluxdVectorCurl()[1][curlField];
-                const STensor3 &dFluxdVectorPotential = ipv->getConstRefTodFluxdVectorPotential()[1][curlField];
 
-                /*STensor3 dSourceVectorFielddB(0.0);
-                STensor3 dMechFieldSourcedB(0.0);
+                STensor3 dSourceVectorFielddVectorCurl(0.0);
+                STensor3 dSourceVectorFielddVectorPotential(0.0);
 
                 if (getConstitutiveCurlAccountFieldSource())
                 {
-                    dSourceVectorFielddB += ipv->getConstRefTodSourceVectorFielddVectorCurl()[curlField][curlField2];
+                    dSourceVectorFielddVectorCurl += ipv->getConstRefTodSourceVectorFielddVectorCurl()[curlField][curlField2];
+                    dSourceVectorFielddVectorPotential += ipv->getConstRefTodSourceVectorFielddVectorPotential()[curlField][curlField2];
                 }
-                if (getConstitutiveCurlAccountMecaSource())
-                {
-                    dMechFieldSourcedB += ipv->getConstRefTodMechanicalFieldSourcedVectorCurl()[curlField][curlField2];
-                }*/
 
                 for (unsigned int j = 0; j < numEdgeShapeFuncRow; ++j)
                 {
@@ -1813,16 +1807,16 @@ void dG3DStiffnessBulk::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double>
                                          dHdB.operator()(m, n) *
                                          CurlcurlVals[k + numEdgeShapeFuncTotalLastCol][n]);
                                 mStiff(j + numEdgeShapeFuncTotalLastRow + numStandardDof,
-                                       k + numEdgeShapeFuncTotalLastCol + numStandardDof) +=
+                                       k + numEdgeShapeFuncTotalLastCol + numStandardDof) -=
                                         getConstitutiveCurlEqRatio() * ratio *
                                         (curlVals[j + numEdgeShapeFuncTotalLastRow][m] *
-                                         dFluxdVectorPotential.operator()(m, n) *
+                                        dSourceVectorFielddVectorPotential.operator()(m, n) *
                                          curlVals[k + numEdgeShapeFuncTotalLastCol][n]);
                                 mStiff(j + numEdgeShapeFuncTotalLastRow + numStandardDof,
-                                       k + numEdgeShapeFuncTotalLastCol + numStandardDof) +=
+                                       k + numEdgeShapeFuncTotalLastCol + numStandardDof) -=
                                         getConstitutiveCurlEqRatio() * ratio *
                                         (curlVals[j + numEdgeShapeFuncTotalLastRow][m] *
-                                         dFluxdVectorCurl.operator()(m, n) *
+                                        dSourceVectorFielddVectorCurl.operator()(m, n) *
                                          CurlcurlVals[k + numEdgeShapeFuncTotalLastCol][n]);
                             }
                         }
@@ -1890,17 +1884,12 @@ void dG3DStiffnessBulk::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double>
 
                       const SVector3 &dVectorFielddNonLocalVar = ipv->getConstRefTodVectorFielddNonLocalVariable()[curlField][nlk];
 
-                      /*SVector3 dSourceVectorFielddNonLocalVar(0.0);
-                      SVector3 dMechFieldSourcedNonLocalVar(0.0);
+                      SVector3 dSourceVectorFielddNonLocalVar(0.0);
 
                       if (getConstitutiveCurlAccountFieldSource())
                       {
                           dSourceVectorFielddNonLocalVar += ipv->getConstRefTodSourceVectorFielddNonLocalVariable()[curlField][nlk];
                       }
-                      if (getConstitutiveCurlAccountMecaSource())
-                      {
-                          dMechFieldSourcedNonLocalVar += ipv->getConstRefTodMechanicalFieldSourcedNonLocalVariable()[curlField][nlk];
-                      }*/
 
                       for (unsigned int j = 0; j < numEdgeShapeFuncRow; ++j)
                           for (unsigned int k = 0; k < nbFFCol; ++k)
@@ -1931,24 +1920,18 @@ void dG3DStiffnessBulk::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double>
                       const unsigned int nbFFCol = nlspace->getNumShapeFunctions(ele, IndexField);
                       const unsigned int nbFFTotalLastCol = nlspace->getShapeFunctionsIndex(ele, IndexField);
 
+                      // dHdV, dHdT, dHdGradV, dHdGradT
                       const SVector3 &dVectorFielddExtraDofField = ipv->getConstRefTodVectorFielddExtraDofField()[curlField][extraDOFField];
                       const STensor3 &dVectorFielddGradExtraDofField = ipv->getConstRefTodVectorFielddGradExtraDofField()[curlField][extraDOFField];
 
-                      // only dFluxDdField, dFluxDdGradField contributions
-                      const SVector3 &dFluxdExtraDofField = ipv->getConstRefTodFluxdField()[1][extraDOFField];
-                      const STensor3 &dFluxdGradExtraDofField = ipv->getConstRefTodFluxdGradField()[1][extraDOFField];
-
-                      /*SVector3 dSourceVectorFielddExtraDofField(0.0);
-                      SVector3 dMechFieldSourcedExtraDofField(0.0);
+                      SVector3 dSourceVectorFielddExtraDofField(0.0);
+                      STensor3 dSourceVectorFielddGradExtraDofField(0.0);
 
                       if (getConstitutiveCurlAccountFieldSource())
                       {
                           dSourceVectorFielddExtraDofField += ipv->getConstRefTodSourceVectorFielddExtraDofField()[curlField][extraDOFField];
+                          dSourceVectorFielddGradExtraDofField += ipv->getConstRefTodSourceVectorFielddGradExtraDofField()[curlField][extraDOFField];
                       }
-                      if (getConstitutiveCurlAccountMecaSource())
-                      {
-                          dMechFieldSourcedExtraDofField += ipv->getConstRefTodMechanicalFieldSourcedExtraDofField()[curlField][extraDOFField];
-                      }*/
 
                       for (unsigned int j = 0; j < numEdgeShapeFuncRow; ++j)
                           for (unsigned int k = 0; k < nbFFCol; ++k)
@@ -1960,7 +1943,7 @@ void dG3DStiffnessBulk::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double>
                                                                                                               CurlcurlVals[j+numEdgeShapeFuncTotalLastRow][m] *
                                                                                                               Vals[k+nbFFTotalLastCol]);
                                   mStiff(j+numEdgeShapeFuncTotalLastRow+numStandardDof,k+nbFFTotalLastCol) -= getConstitutiveCurlEqRatio() * ratio *
-                                                                                                              (dFluxdExtraDofField.operator()(m) *
+                                                                                                              (dSourceVectorFielddExtraDofField.operator()(m) *
                                                                                                                curlVals[j+numEdgeShapeFuncTotalLastRow][m] *
                                                                                                                Vals[k+nbFFTotalLastCol]);
                                   for (unsigned int n = 0; n < 3; ++n)
@@ -1971,7 +1954,7 @@ void dG3DStiffnessBulk::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double>
                                                                                                                   Grads[k+nbFFTotalLastCol][n]);
                                       mStiff(j+numEdgeShapeFuncTotalLastRow+numStandardDof,k+nbFFTotalLastCol) -= getConstitutiveCurlEqRatio() * ratio *
                                                                                                                   (curlVals[j+numEdgeShapeFuncTotalLastRow][m] *
-                                                                                                                   dFluxdGradExtraDofField.operator()(m,n) *
+                                                                                                                  dSourceVectorFielddGradExtraDofField.operator()(m,n) *
                                                                                                                    Grads[k+nbFFTotalLastCol][n]);
                                   }
                               }
@@ -1981,7 +1964,7 @@ void dG3DStiffnessBulk::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double>
           }
       }
     }
-  mStiff.print("bulk");
+  //mStiff.print("bulk");
   } // end if sym
   else
     printf("not implemented\n");