diff --git a/NonLinearSolver/materialLaw/mlawLinearElecMagTherMech.cpp b/NonLinearSolver/materialLaw/mlawLinearElecMagTherMech.cpp
index 7fd097bf7ef8fb4305f36f466420d3153958be3d..7d6e5dc1c0f0716ce34ed3044d9f2ce774ab2f52 100644
--- a/NonLinearSolver/materialLaw/mlawLinearElecMagTherMech.cpp
+++ b/NonLinearSolver/materialLaw/mlawLinearElecMagTherMech.cpp
@@ -339,7 +339,7 @@ void mlawLinearElecMagTherMech::constitutive(
                 // j_e, q, j_y
                 fluxD(i) += (_l0(i, j) * (An(j) - A0(j)) * invdh);
                 fluxT(i) += ((_l0(i, j) * _seebeck * T) * (An(j) - A0(j)) * invdh);
-                fluxjy(i) += ((_l0(i, j) * V + _l0(i, j) * _seebeck * T) * (An(j) - A0(j)) * invdh);
+                fluxjy(i) = 0.0; //((_l0(i, j) * V + _l0(i, j) * _seebeck * T) * (An(j) - A0(j)) * invdh);
 
                 if (stiff)
                 {
@@ -356,9 +356,9 @@ 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) + _seebeck * T * dldT(i,j) + V * dldT(i,j)) * (An(j) - A0(j)) * invdh;
-                    djydV(i) += _l0(i, j) * (An(j) - A0(j)) * invdh;
-                    djydA(i, j) += (_l0(i, j) * V + _l0(i, j) * _seebeck * T) * invdh;
+                    djydT(i) = 0.0; // (_seebeck * _l0(i, j) + dseebeckdT * T * _l0(i, j) + _seebeck * T * dldT(i,j) + V * dldT(i,j)) * (An(j) - A0(j)) * invdh;
+                    djydV(i) = 0.0; // _l0(i, j) * (An(j) - A0(j)) * invdh;
+                    djydA(i, j) = 0.0; // (_l0(i, j) * V + _l0(i, j) * _seebeck * T) * invdh;
 
 
                     // derivatives of sourceVectorField
@@ -375,7 +375,7 @@ void mlawLinearElecMagTherMech::constitutive(
             // -j_e
             sourceVectorField(i) = -fluxD(i);
             // w_AV = j_e \cdot -dadt + h \cdot dbdt
-            emFieldSource += (sourceVectorField(i) * ((An(i) - A0(i)) * invdh) /*+ H(i) * (Bn(i) - B0(i)) * invdh*/);
+            emFieldSource += (sourceVectorField(i) * (((An(i) - A0(i)) * invdh) + gradV(i)) /*+ H(i) * (Bn(i) - B0(i)) * invdh*/);
         }
 
         // Need to check contribution to jy from source js0
@@ -383,7 +383,7 @@ void mlawLinearElecMagTherMech::constitutive(
         {
             fluxD(i) += 0.0;
             fluxT(i) += 0.0;
-            fluxjy(i) += (_seebeck * T + V) * js0(i);
+            fluxjy(i) = 0.0; // (_seebeck * T + V) * js0(i);
 
             if (stiff)
             {
@@ -393,7 +393,7 @@ void mlawLinearElecMagTherMech::constitutive(
                 dqdT(i) += 0.0;
                 dqdV(i) += 0.0;
 
-                djydT(i) += (_seebeck + dseebeckdT * T) * js0(i);
+                djydT(i) = 0.0; // (_seebeck + dseebeckdT * T) * js0(i);
 
                 // derivatives of sourceVectorField
                 dsourceVectorFielddV(i) += 0.0;
@@ -410,12 +410,12 @@ void mlawLinearElecMagTherMech::constitutive(
                     dsourceVectorFielddgradV(i, j) += 0.0;
                     dsourceVectorFielddgradT(i, j) += 0.0;
                 }
-                djydV(i) += js0(i);
+                djydV(i) = 0.0; // js0(i);
             }
             // -j_e
             sourceVectorField(i) = - js0(i);
             // w_AV = j_e \cdot -dadt + h \cdot dbdt
-            emFieldSource += (sourceVectorField(i) * ((An(i) - A0(i)) * invdh) /*+ H(i) * (Bn(i) - B0(i)) * invdh*/);
+            emFieldSource += (sourceVectorField(i) * (((An(i) - A0(i)) * invdh) + gradV(i)) /*+ H(i) * (Bn(i) - B0(i)) * invdh*/);
         }
 
         // For test with unit cube to check
@@ -470,7 +470,7 @@ void mlawLinearElecMagTherMech::constitutive(
             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)*/ 0.0;
             demFieldSourcedGradT(m) +=  (-dedgradT(m,n) * (An(n) - A0(n)) * invdh);
-            demFieldSourcedGradV(m) += (-dedgradV(m,n) * (An(n) - A0(n)) * invdh);
+            demFieldSourcedGradV(m) += (-dedgradV(m,n) * (An(n) - A0(n)) * invdh - dedgradV(m,n) * gradV(n) - fluxD(m));
         }
         demFieldSourcedT += (-dedT(m) * (An(m) - A0(m)) * invdh /*+ dnudT(m,n) * B(n) * (Bn(m) - B0(m)) * invdh*/);
         demFieldSourcedV += 0.0;
diff --git a/dG3D/src/dG3DIPVariable.h b/dG3D/src/dG3DIPVariable.h
index bcdee4f3a0bb26aa1984ba1f2d2855599a3faa63..774ce6070486ca8037d5efd5f8dae3ef3e654799 100644
--- a/dG3D/src/dG3DIPVariable.h
+++ b/dG3D/src/dG3DIPVariable.h
@@ -3126,9 +3126,11 @@ class ElecTherMechDG3DIPVariableBase : public dG3DIPVariable
   virtual SVector3 &getRefToGradT() {return getRefToGradField()[0];}
 
   virtual const SVector3 &getConstRefToFluxT() const
-      {Msg::Error("ElecTherMechDG3DIPVariableBase::getConstRefToFluxT should not be used; use getConstRefToFluxEnergy"); static SVector3 a; return a;}
+  { return getConstRefToFlux()[0];}
+      //{Msg::Error("ElecTherMechDG3DIPVariableBase::getConstRefToFluxT should not be used; use getConstRefToFluxEnergy"); static SVector3 a; return a;}
   virtual SVector3 &getRefToFluxT()
-      {Msg::Error("ElecTherMechDG3DIPVariableBase::getRefToFluxT should not be used; use getRefToFluxEnergy");static SVector3 a; return a;}
+      {return getRefToFlux()[0];}
+      //{Msg::Error("ElecTherMechDG3DIPVariableBase::getRefToFluxT should not be used; use getRefToFluxEnergy");static SVector3 a; return a;}
 
   virtual  const STensor3 &getConstRefTodPdT() const {return getConstRefTodPdField()[0];}
   virtual  STensor3 &getRefTodPdT() {return getRefTodPdField()[0];}
@@ -3138,13 +3140,15 @@ class ElecTherMechDG3DIPVariableBase : public dG3DIPVariable
 
   virtual const SVector3 & getConstRefTodFluxTdT() const
   {
-    Msg::Error("ElecTherMechDG3DIPVariableBase::getConstRefTodFluxTdT should not be used; use getConstRefTodFluxEnergydT");
-    static SVector3 a; return a;
+    //Msg::Error("ElecTherMechDG3DIPVariableBase::getConstRefTodFluxTdT should not be used; use getConstRefTodFluxEnergydT");
+    //static SVector3 a; return a;
+    return getConstRefTodFluxdField()[0][0];
   }
   virtual SVector3 & getRefTodFluxTdT()
   {
-    Msg::Error("ElecTherMechDG3DIPVariableBase::getRefTodFluxTdT should not be used; use getRefTodFluxEnergydT");
-    static SVector3 a; return a;
+    //Msg::Error("ElecTherMechDG3DIPVariableBase::getRefTodFluxTdT should not be used; use getRefTodFluxEnergydT");
+    //static SVector3 a; return a;
+    return getRefTodFluxdField()[0][0];
   }
 
   virtual const  SVector3 &getConstRefTodFluxTdV() const {return  getConstRefTodFluxdField()[0][1];}
@@ -3152,13 +3156,15 @@ class ElecTherMechDG3DIPVariableBase : public dG3DIPVariable
 
   virtual const  STensor3 & getConstRefTodFluxTdGradT() const
   {
-    Msg::Error("ElecTherMechDG3DIPVariableBase::getConstRefTodFluxTdGradT should not be used; use getConstRefTodFluxEnergydGradT");
-    static STensor3 a; return a;
+    //Msg::Error("ElecTherMechDG3DIPVariableBase::getConstRefTodFluxTdGradT should not be used; use getConstRefTodFluxEnergydGradT");
+    //static STensor3 a; return a;
+    return getConstRefTodFluxdGradField()[0][0];
   }//dqdgradT,
   virtual  STensor3 &getRefTodFluxTdGradT()
   {
-    Msg::Error("ElecTherMechDG3DIPVariableBase::getRefTodFluxTdGradT should not be used; use getRefTodFluxEnergydGradT");
-    static STensor3 a; return a;
+    //Msg::Error("ElecTherMechDG3DIPVariableBase::getRefTodFluxTdGradT should not be used; use getRefTodFluxEnergydGradT");
+    //static STensor3 a; return a;
+    return getRefTodFluxdGradField()[0][0];
   }
 
   virtual const  STensor3 &getConstRefTodFluxTdGradV() const {return  getConstRefTodFluxdGradField()[0][1];}//dqdgradV,
@@ -3166,13 +3172,15 @@ class ElecTherMechDG3DIPVariableBase : public dG3DIPVariable
 
   virtual const  STensor33 &getConstRefTodFluxTdF() const
   {
-    Msg::Error("ElecTherMechDG3DIPVariableBase::getConstRefTodFluxTdF should not be used; use getConstRefTodFluxEnergydF");
-    static STensor33 a; return a;
+    //Msg::Error("ElecTherMechDG3DIPVariableBase::getConstRefTodFluxTdF should not be used; use getConstRefTodFluxEnergydF");
+    //static STensor33 a; return a;
+    return getConstRefTodFluxdF()[0];
   }
   virtual  STensor33 &getRefTodFluxTdF()
   {
-    Msg::Error("ElecTherMechDG3DIPVariableBase::getRefTodFluxTdF should not be used; use getRefTodFluxEnergydF");
-    static STensor33 a; return a;
+    //Msg::Error("ElecTherMechDG3DIPVariableBase::getRefTodFluxTdF should not be used; use getRefTodFluxEnergydF");
+    //static STensor33 a; return a;
+    return getRefTodFluxdF()[0];
   }
 
 /* virtual const  SVector3 &getConstRefTodFluxTdT() const {return  getConstRefTodFluxdField()[0][0];}
@@ -3377,23 +3385,23 @@ class ElecTherMechDG3DIPVariableBase : public dG3DIPVariable
   virtual const SVector3  &getConstRefTodgradTdfV() const {return getConstRefTodGradFielddEnergyConjugatedField()[0][1];}
   virtual SVector3  &getRefTodgradTdfV() {return  getRefTodGradFielddEnergyConjugatedField()[0][1];}
 
-   virtual const SVector3  &getConstRefToFluxEnergy() const {return getConstRefToFlux()[0];}
-   virtual SVector3 &getRefToFluxEnergy(){return  getRefToFlux()[0];}
+   virtual const SVector3  &getConstRefToFluxEnergy() const {static SVector3 a; return a;}//return getConstRefToFlux()[0];}
+   virtual SVector3 &getRefToFluxEnergy(){static SVector3 a; return a;}//return  getRefToFlux()[0];}
 
-   virtual const SVector3  &getConstRefTodFluxEnergydV() const {return getConstRefTodFluxdField()[0][1];}
-   virtual SVector3 &getRefTodFluxEnergydV(){return  getRefTodFluxdField()[0][1];}
+   virtual const SVector3  &getConstRefTodFluxEnergydV() const {static SVector3 a; return a;}//return getConstRefTodFluxdField()[0][1];}
+   virtual SVector3 &getRefTodFluxEnergydV(){static SVector3 a; return a;}//return  getRefTodFluxdField()[0][1];}
 
-   virtual const STensor3 &getConstRefTodFluxEnergydGradV() const {return getConstRefTodFluxdGradField()[0][1];}
-   virtual STensor3 &getRefTodFluxEnergydGradV(){return  getRefTodFluxdGradField()[0][1];}
+   virtual const STensor3 &getConstRefTodFluxEnergydGradV() const {static STensor3 a; return a;}//return getConstRefTodFluxdGradField()[0][1];}
+   virtual STensor3 &getRefTodFluxEnergydGradV(){static STensor3 a; return a;}//return  getRefTodFluxdGradField()[0][1];}
 
-   virtual const SVector3  &getConstRefTodFluxEnergydT() const {return getConstRefTodFluxdField()[0][0];}
-   virtual SVector3 &getRefTodFluxEnergydT(){return  getRefTodFluxdField()[0][0];}
+   virtual const SVector3  &getConstRefTodFluxEnergydT() const {static SVector3 a; return a;}//return getConstRefTodFluxdField()[0][0];}
+   virtual SVector3 &getRefTodFluxEnergydT(){static SVector3 a; return a;}//return  getRefTodFluxdField()[0][0];}
 
-   virtual const STensor3 &getConstRefTodFluxEnergydGradT() const {return getConstRefTodFluxdGradField()[0][0];}
-   virtual STensor3 &getRefTodFluxEnergydGradT(){return  getRefTodFluxdGradField()[0][0];}
+   virtual const STensor3 &getConstRefTodFluxEnergydGradT() const {static STensor3 a; return a;}//return getConstRefTodFluxdGradField()[0][0];}
+   virtual STensor3 &getRefTodFluxEnergydGradT(){static STensor3 a; return a;}//return  getRefTodFluxdGradField()[0][0];}
 
-   virtual const STensor33 &getConstRefTodFluxEnergydF() const {return getConstRefTodFluxdF()[0];}
-   virtual STensor33 &getRefTodFluxEnergydF(){return  getRefTodFluxdF()[0];}
+   virtual const STensor33 &getConstRefTodFluxEnergydF() const {static STensor33 a; return a;}//return getConstRefTodFluxdF()[0];}
+   virtual STensor33 &getRefTodFluxEnergydF(){static STensor33 a; return a;}//return  getRefTodFluxdF()[0];}
 
  /*  virtual const STensor3 *getConstRefToEnergyConductivity() const {return getConstRefToEnergyK()[0][0];}
    virtual void setRefToEnergyConductivity(const STensor3 &eT){return  setConstRefToEnergyK(eT,0,0);}
@@ -3501,13 +3509,15 @@ class ElecMagTherMechDG3DIPVariableBase : public dG3DIPVariable
     // FluxT = q
     virtual const SVector3 &getConstRefToFluxT() const
     {
-        Msg::Error("ElecMagTherMechDG3DIPVariableBase::getConstRefToFluxT should not be used; use getConstRefToFluxEnergy");
-        static SVector3 a; return a;
+        //Msg::Error("ElecMagTherMechDG3DIPVariableBase::getConstRefToFluxT should not be used; use getConstRefToFluxEnergy");
+        //static SVector3 a; return a;
+        return getConstRefToFlux()[0];
     }
     virtual SVector3 &getRefToFluxT()
     {
-        Msg::Error("ElecMagTherMechDG3DIPVariableBase::getRefToFluxT should not be used; use getRefToFluxEnergy");
-        static SVector3 a; return a;
+        //Msg::Error("ElecMagTherMechDG3DIPVariableBase::getRefToFluxT should not be used; use getRefToFluxEnergy");
+        //static SVector3 a; return a;
+        return getRefToFlux()[0];
     }
 
     // dFluxTdV
@@ -3521,37 +3531,43 @@ class ElecMagTherMechDG3DIPVariableBase : public dG3DIPVariable
     // dFluxTdT
     virtual const SVector3 & getConstRefTodFluxTdT() const
     {
-        Msg::Error("ElecMagTherMechDG3DIPVariableBase::getConstRefTodFluxTdT should not be used; use getConstRefTodFluxEnergydT");
-        static SVector3 a; return a;
+        //Msg::Error("ElecMagTherMechDG3DIPVariableBase::getConstRefTodFluxTdT should not be used; use getConstRefTodFluxEnergydT");
+        //static SVector3 a; return a;
+        return getConstRefTodFluxdField()[0][0];
     }
     virtual SVector3 & getRefTodFluxTdT()
     {
-        Msg::Error("ElecMagTherMechDG3DIPVariableBase::getRefTodFluxTdT should not be used; use getRefTodFluxEnergydT");
-        static SVector3 a; return a;
+        //Msg::Error("ElecMagTherMechDG3DIPVariableBase::getRefTodFluxTdT should not be used; use getRefTodFluxEnergydT");
+        //static SVector3 a; return a;
+        return getRefTodFluxdField()[0][0];
     }
 
     // dFluxTdGradT
     virtual const STensor3 & getConstRefTodFluxTdGradT() const
     {
-        Msg::Error("ElecMagTherMechDG3DIPVariableBase::getConstRefTodFluxTdGradT should not be used; use getConstRefTodFluxEnergydGradT");
-        static STensor3 a; return a;
+        //Msg::Error("ElecMagTherMechDG3DIPVariableBase::getConstRefTodFluxTdGradT should not be used; use getConstRefTodFluxEnergydGradT");
+        //static STensor3 a; return a;
+        return getConstRefTodFluxdGradField()[0][0];
     }
     virtual STensor3 & getRefTodFluxTdGradT()
     {
-        Msg::Error("ElecMagTherMechDG3DIPVariableBase::getRefTodFluxTdGradT should not be used; use getRefTodFluxEnergydGradT");
-        static STensor3 a; return a;
+        //Msg::Error("ElecMagTherMechDG3DIPVariableBase::getRefTodFluxTdGradT should not be used; use getRefTodFluxEnergydGradT");
+        //static STensor3 a; return a;
+        return getRefTodFluxdGradField()[0][0];
     }
 
     // dFluxTdF
     virtual const STensor33 & getConstRefTodFluxTdF() const
     {
-        Msg::Error("ElecMagTherMechDG3DIPVariableBase::getConstRefTodFluxTdF should not be used; use getConstRefTodFluxEnergydF");
-        static STensor33 a; return a;
+        //Msg::Error("ElecMagTherMechDG3DIPVariableBase::getConstRefTodFluxTdF should not be used; use getConstRefTodFluxEnergydF");
+        //static STensor33 a; return a;
+        return getConstRefTodFluxdF()[0];
     }
     virtual STensor33 & getRefTodFluxTdF()
     {
-        Msg::Error("ElecMagTherMechDG3DIPVariableBase::getRefTodFluxTdF should not be used; use getRefTodFluxEnergydF");
-        static STensor33 a; return a;
+        //Msg::Error("ElecMagTherMechDG3DIPVariableBase::getRefTodFluxTdF should not be used; use getRefTodFluxEnergydF");
+        //static STensor33 a; return a;
+        return getRefTodFluxdF()[0];
     }
 
     // ElecCond
@@ -3787,36 +3803,36 @@ class ElecMagTherMechDG3DIPVariableBase : public dG3DIPVariable
     virtual SVector3  &getRefTodgradTdfV() {return  getRefTodGradFielddEnergyConjugatedField()[0][1];}
 
     // FluxEnergy
-    virtual const SVector3  &getConstRefToFluxEnergy() const {return getConstRefToFlux()[0];}
-    virtual SVector3 &getRefToFluxEnergy(){return  getRefToFlux()[0];}
+    virtual const SVector3  &getConstRefToFluxEnergy() const {static SVector3 a; return a;}//return getConstRefToFlux()[0];}
+    virtual SVector3 &getRefToFluxEnergy(){static SVector3 a; return a;}//return  getRefToFlux()[0];}
 
     // dFluxEnergydV
-    virtual const SVector3  &getConstRefTodFluxEnergydV() const {return getConstRefTodFluxdField()[0][1];}
-    virtual SVector3 &getRefTodFluxEnergydV(){return  getRefTodFluxdField()[0][1];}
+    virtual const SVector3  &getConstRefTodFluxEnergydV() const {static SVector3 a; return a;}//return getConstRefTodFluxdField()[0][1];}
+    virtual SVector3 &getRefTodFluxEnergydV(){static SVector3 a; return a;}//return  getRefTodFluxdField()[0][1];}
 
     // dFluxEnergydGradV
-    virtual const STensor3 &getConstRefTodFluxEnergydGradV() const {return getConstRefTodFluxdGradField()[0][1];}
-    virtual STensor3 &getRefTodFluxEnergydGradV(){return  getRefTodFluxdGradField()[0][1];}
+    virtual const STensor3 &getConstRefTodFluxEnergydGradV() const {static STensor3 a; return a;}//return getConstRefTodFluxdGradField()[0][1];}
+    virtual STensor3 &getRefTodFluxEnergydGradV(){static STensor3 a; return a;}//return  getRefTodFluxdGradField()[0][1];}
 
     // dFluxEnergydT
-    virtual const SVector3  &getConstRefTodFluxEnergydT() const {return getConstRefTodFluxdField()[0][0];}
-    virtual SVector3 &getRefTodFluxEnergydT(){return  getRefTodFluxdField()[0][0];}
+    virtual const SVector3  &getConstRefTodFluxEnergydT() const {static SVector3 a; return a;}//return getConstRefTodFluxdField()[0][0];}
+    virtual SVector3 &getRefTodFluxEnergydT(){static SVector3 a; return a;}//return  getRefTodFluxdField()[0][0];}
 
     // dFluxEnergydGradT
-    virtual const STensor3 &getConstRefTodFluxEnergydGradT() const {return getConstRefTodFluxdGradField()[0][0];}
-    virtual STensor3 &getRefTodFluxEnergydGradT(){return  getRefTodFluxdGradField()[0][0];}
+    virtual const STensor3 &getConstRefTodFluxEnergydGradT() const {static STensor3 a; return a;}//return getConstRefTodFluxdGradField()[0][0];}
+    virtual STensor3 &getRefTodFluxEnergydGradT(){static STensor3 a; return a;}//return  getRefTodFluxdGradField()[0][0];}
 
     // dFluxEnergydF
-    virtual const STensor33 &getConstRefTodFluxEnergydF() const {return getConstRefTodFluxdF()[0];}
-    virtual STensor33 &getRefTodFluxEnergydF(){return  getRefTodFluxdF()[0];}
+    virtual const STensor33 &getConstRefTodFluxEnergydF() const {static STensor33 a; return a;}//return getConstRefTodFluxdF()[0];}
+    virtual STensor33 &getRefTodFluxEnergydF(){static STensor33 a; return a;}//return  getRefTodFluxdF()[0];}
 
     // dFluxEnergydA
-    virtual const STensor3 & getConstRefTodFluxEnergydMagneticVectorPotential() const {return getConstRefTodFluxdVectorPotential()[0][0];}
-    virtual STensor3 & getRefTodFluxEnergydMagneticVectorPotential() {return getRefTodFluxdVectorPotential()[0][0];}
+    virtual const STensor3 & getConstRefTodFluxEnergydMagneticVectorPotential() const {static STensor3 a; return a;}//return getConstRefTodFluxdVectorPotential()[0][0];}
+    virtual STensor3 & getRefTodFluxEnergydMagneticVectorPotential() {static STensor3 a; return a;}//return getRefTodFluxdVectorPotential()[0][0];}
 
     // dFluxEnergydB
-    virtual const STensor3 & getConstRefTodFluxEnergydMagneticVectorCurl() const {return getConstRefTodFluxdVectorCurl()[0][0];}
-    virtual STensor3 & getRefTodFluxEnergydMagneticVectorCurl() {return getRefTodFluxdVectorCurl()[0][0];}
+    virtual const STensor3 & getConstRefTodFluxEnergydMagneticVectorCurl() const {static STensor3 a; return a;}//return getConstRefTodFluxdVectorCurl()[0][0];}
+    virtual STensor3 & getRefTodFluxEnergydMagneticVectorCurl() {static STensor3 a; return a;}//return getRefTodFluxdVectorCurl()[0][0];}
 
     // dFluxEnergydGradVdT
     virtual const STensor3 &getConstRefTodFluxEnergydGradVdT() const {return getConstRefTodFluxdGradFielddField()[0][1][0];}
diff --git a/dG3D/src/dG3DMaterialLaw.cpp b/dG3D/src/dG3DMaterialLaw.cpp
index 9568a3fd50ea8e8900143f1686b3049504304541..3d575c3e4631bc58563fcc834edfdbe80ef376ff 100644
--- a/dG3D/src/dG3DMaterialLaw.cpp
+++ b/dG3D/src/dG3DMaterialLaw.cpp
@@ -5477,17 +5477,17 @@ void LinearElecTherMechDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable
 
     STensor3& dPdT=ipvcur->getRefTodPdT();
 	STensor3& dPdV=ipvcur->getRefTodPdV();
-    SVector3 fluxT;
+    SVector3& fluxT=ipvcur->getRefToFluxT();
 	SVector3& fluxV=ipvcur->getRefToFluxV();
-    STensor3  dqdgradT;
-    STensor3  dqdgradV;
+    STensor3&  dqdgradT=ipvcur->getRefTodFluxTdGradT();
+    STensor3&  dqdgradV=ipvcur->getRefTodFluxTdGradV();
 	STensor3& dedgradV=ipvcur->getRefTodFluxVdGradV();
 	STensor3& dedgradT=ipvcur->getRefTodFluxVdGradT();
-    SVector3  dqdT;
-	SVector3  dqdV;
+    SVector3&  dqdT=ipvcur->getRefTodFluxTdT();
+	SVector3&  dqdV=ipvcur->getRefTodFluxTdV();
 	SVector3& dedV=ipvcur->getRefTodFluxVdV();
 	SVector3& dedT=ipvcur->getRefTodFluxVdT();
-    STensor33  dqdF;
+    STensor33&  dqdF=ipvcur->getRefTodFluxTdF();
 	STensor33& dedF=ipvcur->getRefTodFluxVdF();
         STensor3& P=ipvcur->getRefToFirstPiolaKirchhoffStress();
         STensor43& Tangent=ipvcur->getRefToTangentModuli();
@@ -5512,24 +5512,24 @@ void LinearElecTherMechDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable
 	STensor43& dk10dF=ipvcur->getRefTodTherConductivitydF();
 	STensor43& dk20dF=ipvcur->getRefTodCrossTherConductivitydF();
 
-	SVector3 &fluxjy=ipvcur->getRefToFluxEnergy();
-	SVector3 &djydV=ipvcur->getRefTodFluxEnergydV();
-	STensor3 &djydgradV=ipvcur->getRefTodFluxEnergydGradV();
-	SVector3 &djydT=ipvcur->getRefTodFluxEnergydT();
-	STensor3 &djydgradT=ipvcur->getRefTodFluxEnergydGradT();
-        STensor33 &djydF=ipvcur->getRefTodFluxEnergydF();
+	SVector3 fluxjy;
+	SVector3 djydV;
+	STensor3 djydgradV;
+	SVector3 djydT;
+	STensor3 djydgradT;
+        STensor33 djydF;
 
         /*ipvcur->setRefToEnergyConductivity(this->jy1);
 	STensor3 &djy1dT= ipvcur->getRefTodEnergyConductivitydT();
 	STensor3 &djy1dV= ipvcur->getRefTodEnergyConductivitydV();
 	STensor43 &djy1dF=ipvcur->getRefTodEnergyConductivitydF();*/
 
-	STensor3 &djydgradVdT= ipvcur->getRefTodFluxEnergydGradVdT();
-	STensor3 &djydgradVdV= ipvcur->getRefTodFluxEnergydGradVdV();
-	STensor43 &djydgradVdF=ipvcur->getRefTodFluxEnergydGradVdF();
-	STensor3 &djydgradTdT= ipvcur->getRefTodFluxEnergydGradTdT();
-	STensor3 &djydgradTdV= ipvcur->getRefTodFluxEnergydGradTdV();
-	STensor43 &djydgradTdF=ipvcur->getRefTodFluxEnergydGradTdF();
+	STensor3 djydgradVdT;
+	STensor3 djydgradVdV;
+	STensor43 djydgradVdF;
+	STensor3 djydgradTdT;
+	STensor3 djydgradTdV;
+	STensor43 djydgradTdF;
 
   /* data for J2 law */
   IPLinearElecTherMech* q1 = ipvcur->getIPLinearElecTherMech();
@@ -5667,17 +5667,17 @@ const bool checkfrac, const bool dTangent)
 
     STensor3& dPdT=ipvcur->getRefTodPdT();
     STensor3& dPdV=ipvcur->getRefTodPdV();
-    SVector3 fluxT;
+    SVector3& fluxT=ipvcur->getRefToFluxT();
     SVector3& fluxD=ipvcur->getRefToFluxD();
-    STensor3 dqdgradT;
-    STensor3 dqdgradV;
+    STensor3& dqdgradT=ipvcur->getRefTodFluxTdGradT();
+    STensor3& dqdgradV=ipvcur->getRefTodFluxTdGradV();
     STensor3& dedgradV=ipvcur->getRefTodFluxDdGradV();
     STensor3& dedgradT=ipvcur->getRefTodFluxDdGradT();
-    SVector3 dqdT;
-    SVector3 dqdV;
+    SVector3& dqdT=ipvcur->getRefTodFluxTdT();
+    SVector3& dqdV=ipvcur->getRefTodFluxTdV();
     SVector3& dedV=ipvcur->getRefTodFluxDdV();
     SVector3& dedT=ipvcur->getRefTodFluxDdT();
-    STensor33 dqdF;
+    STensor33& dqdF=ipvcur->getRefTodFluxTdF();
     STensor33& dedF=ipvcur->getRefTodFluxDdF();
     STensor3& P=ipvcur->getRefToFirstPiolaKirchhoffStress();
     STensor43& Tangent=ipvcur->getRefToTangentModuli();
@@ -5703,21 +5703,21 @@ const bool checkfrac, const bool dTangent)
     STensor43& dk10dF=ipvcur->getRefTodTherConductivitydF();
     STensor43& dk20dF=ipvcur->getRefTodCrossTherConductivitydF();
 
-    SVector3 &fluxjy=ipvcur->getRefToFluxEnergy();
-    SVector3 &djydV=ipvcur->getRefTodFluxEnergydV();
-    STensor3 &djydgradV=ipvcur->getRefTodFluxEnergydGradV();
-    SVector3 &djydT=ipvcur->getRefTodFluxEnergydT();
-    STensor3 &djydgradT=ipvcur->getRefTodFluxEnergydGradT();
-    STensor33 &djydF=ipvcur->getRefTodFluxEnergydF();
-    STensor3 & djydA = ipvcur->getRefTodFluxEnergydMagneticVectorPotential();
-    STensor3 & djydB = ipvcur->getRefTodFluxEnergydMagneticVectorCurl();
-
-    STensor3 &djydgradVdT= ipvcur->getRefTodFluxEnergydGradVdT();
-    STensor3 &djydgradVdV= ipvcur->getRefTodFluxEnergydGradVdV();
-    STensor43 &djydgradVdF=ipvcur->getRefTodFluxEnergydGradVdF();
-    STensor3 &djydgradTdT= ipvcur->getRefTodFluxEnergydGradTdT();
-    STensor3 &djydgradTdV= ipvcur->getRefTodFluxEnergydGradTdV();
-    STensor43 &djydgradTdF=ipvcur->getRefTodFluxEnergydGradTdF();
+    SVector3 fluxjy;
+    SVector3 djydV;
+    STensor3 djydgradV;
+    SVector3 djydT;
+    STensor3 djydgradT;
+    STensor33 djydF;
+    STensor3 djydA;
+    STensor3 djydB;
+
+    STensor3 djydgradVdT;
+    STensor3 djydgradVdV;
+    STensor43 djydgradVdF;
+    STensor3 djydgradTdT;
+    STensor3 djydgradTdV;
+    STensor43 djydgradTdF;
 
     /* data for J2 law */
     IPLinearElecMagTherMech* q1 = ipvcur->getIPLinearElecMagTherMech();
@@ -5734,9 +5734,9 @@ const bool checkfrac, const bool dTangent)
     STensor3 dBdGradV;
     STensor33 dBdF;
     STensor3 dBdA;
-    STensor3 dfluxTdA;
+    STensor3& dfluxTdA=ipvcur->getRefTodFluxTdMagneticVectorPotential();
     STensor3 &dfluxDdA = ipvcur->getRefTodFluxDdMagneticVectorPotential();
-    STensor3 dfluxTdB;
+    STensor3& dfluxTdB=ipvcur->getRefTodFluxTdMagneticVectorCurl();
     STensor3 & dfluxDdB = ipvcur->getRefTodFluxDdMagneticVectorCurl();
     SVector3 & dw_TdA = ipvcur->getRefTodFieldSourcedMagneticVectorPotential();
     SVector3 & dw_TdB = ipvcur->getRefTodFieldSourcedMagneticVectorCurl();
@@ -5922,17 +5922,17 @@ const bool checkfrac, const bool dTangent)
 
   STensor3& dPdT=ipvcur->getRefTodPdT();
   STensor3& dPdV=ipvcur->getRefTodPdV();
-  SVector3 fluxT;
+  SVector3& fluxT=ipvcur->getRefToFluxT();
   SVector3& fluxD=ipvcur->getRefToFluxD();
-  STensor3 dqdgradT;
-  STensor3 dqdgradV;
+  STensor3& dqdgradT=ipvcur->getRefTodFluxTdGradT();
+  STensor3& dqdgradV=ipvcur->getRefTodFluxTdGradV();
   STensor3& dedgradV=ipvcur->getRefTodFluxDdGradV();
   STensor3& dedgradT=ipvcur->getRefTodFluxDdGradT();
-  SVector3 dqdT;
-  SVector3 dqdV;
+  SVector3& dqdT=ipvcur->getRefTodFluxTdT();
+  SVector3& dqdV=ipvcur->getRefTodFluxTdV();
   SVector3& dedV=ipvcur->getRefTodFluxDdV();
   SVector3& dedT=ipvcur->getRefTodFluxDdT();
-  STensor33 dqdF;
+  STensor33& dqdF=ipvcur->getRefTodFluxTdF();
   STensor33& dedF=ipvcur->getRefTodFluxDdF();
   STensor3& P=ipvcur->getRefToFirstPiolaKirchhoffStress();
   STensor43& Tangent=ipvcur->getRefToTangentModuli();
@@ -5958,21 +5958,21 @@ const bool checkfrac, const bool dTangent)
   STensor43& dk10dF=ipvcur->getRefTodTherConductivitydF();
   STensor43& dk20dF=ipvcur->getRefTodCrossTherConductivitydF();
 
-  SVector3 &fluxjy=ipvcur->getRefToFluxEnergy();
-  SVector3 &djydV=ipvcur->getRefTodFluxEnergydV();
-  STensor3 &djydgradV=ipvcur->getRefTodFluxEnergydGradV();
-  SVector3 &djydT=ipvcur->getRefTodFluxEnergydT();
-  STensor3 &djydgradT=ipvcur->getRefTodFluxEnergydGradT();
-  STensor33 &djydF=ipvcur->getRefTodFluxEnergydF();
-  STensor3 & djydA = ipvcur->getRefTodFluxEnergydMagneticVectorPotential();
-  STensor3 & djydB = ipvcur->getRefTodFluxEnergydMagneticVectorCurl();
-
-  STensor3 &djydgradVdT= ipvcur->getRefTodFluxEnergydGradVdT();
-  STensor3 &djydgradVdV= ipvcur->getRefTodFluxEnergydGradVdV();
-  STensor43 &djydgradVdF=ipvcur->getRefTodFluxEnergydGradVdF();
-  STensor3 &djydgradTdT= ipvcur->getRefTodFluxEnergydGradTdT();
-  STensor3 &djydgradTdV= ipvcur->getRefTodFluxEnergydGradTdV();
-  STensor43 &djydgradTdF=ipvcur->getRefTodFluxEnergydGradTdF();
+  SVector3 fluxjy;
+  SVector3 djydV;
+  STensor3 djydgradV;
+  SVector3 djydT;
+  STensor3 djydgradT;
+  STensor33 djydF;
+  STensor3  djydA;
+  STensor3  djydB;
+
+  STensor3 djydgradVdT;
+  STensor3 djydgradVdV;
+  STensor43 djydgradVdF;
+  STensor3 djydgradTdT;
+  STensor3 djydgradTdV;
+  STensor43 djydgradTdF;
 
   /* data for J2 law */
   IPLinearElecMagInductor* q1 = ipvcur->getIPLinearElecMagInductor();
@@ -5989,9 +5989,9 @@ const bool checkfrac, const bool dTangent)
   STensor3 dBdGradV;
   STensor33 dBdF;
   STensor3 dBdA;
-  STensor3 dfluxTdA;
+  STensor3 &dfluxTdA=ipvcur->getRefTodFluxTdMagneticVectorPotential();
   STensor3 & dfluxDdA = ipvcur->getRefTodFluxDdMagneticVectorPotential();
-  STensor3 dfluxTdB;
+  STensor3 &dfluxTdB=ipvcur->getRefTodFluxTdMagneticVectorCurl();
   STensor3 & dfluxDdB = ipvcur->getRefTodFluxDdMagneticVectorCurl();
   SVector3 & dw_TdA = ipvcur->getRefTodFieldSourcedMagneticVectorPotential();
   SVector3 & dw_TdB = ipvcur->getRefTodFieldSourcedMagneticVectorCurl();