diff --git a/NonLinearSolver/materialLaw/mlawElecMagGenericThermoMech.cpp b/NonLinearSolver/materialLaw/mlawElecMagGenericThermoMech.cpp
index 4c6357676423b8cecbc9603228c361e582ad7a8a..dcf9e7ce29f166d8fda56c2ed39d0aaadd0085c3 100644
--- a/NonLinearSolver/materialLaw/mlawElecMagGenericThermoMech.cpp
+++ b/NonLinearSolver/materialLaw/mlawElecMagGenericThermoMech.cpp
@@ -359,7 +359,8 @@ void mlawElecMagGenericThermoMech::constitutive(
         double dseebeckdT;
         static STensor3 dldT,dkdT,dnudT,dmudT;//,k1,k2,l10,l20;
 
-        STensor33 d_nu0dB(0.0);
+        static STensor33 d_nu0dB;
+        STensorOperation::zero(d_nu0dB);
 
         STensorOperation::zero(dldT);
         STensorOperation::zero(dkdT);
@@ -475,180 +476,186 @@ void mlawElecMagGenericThermoMech::constitutive(
         STensorOperation::zero(dmutdF);
         STensorOperation::zero(dktdF);
 
-        static STensor3 FRef,FRefInv;
+        static STensor3 FRef,FRefInv,FRefTranspose,FRefInvTranspose;
         static STensor43 dFRefdFn;
         STensorOperation::zero(FRef);
         STensorOperation::zero(FRefInv);
+        STensorOperation::zero(FRefTranspose);
+        STensorOperation::zero(FRefInvTranspose);
         STensorOperation::zero(dFRefdFn);
-        double JacRef(0.0);
-
-        static STensor3 ltTemp, nutTemp, mutTemp, ktTemp;
-        STensorOperation::zero(ltTemp);
-        STensorOperation::zero(nutTemp);
-        STensorOperation::zero(mutTemp);
-        STensorOperation::zero(ktTemp);
+        static double JacRef;
+        STensorOperation::zero(JacRef);
 
         const STensor3 I2(1.0);
         const STensor43 I4(2.0,0.0);
-        const double vacuum_permeability = 4.0 * M_PI * 1.e-7;
-        const double rho = this->_genericThermoMechLaw->getConstRefMechanicalMaterialLaw().density();
+        static const double vacuum_permeability = 4.0 * M_PI * 1.e-7;
+        static const double rho = this->_genericThermoMechLaw->getConstRefMechanicalMaterialLaw().density();
 
-        // strong coupling: FRef = Fn
-        // weak coupling: FRef != Fn (when Fn != I)
-        // FRef = F_TM, Fn = F_EM (= I possibly)
+        // strong coupling: Ftot = Fn Fref (with Fref = I)
+        // weak coupling: Ftot = Fn Fref (with FRef = F_TM, Fn = I)
         if(applyReferenceF)
         {
+            //FRef =Fn; // here we have Fn identity in any case but I would leave it to be correct
             FRef = qTM1.getRefToReferenceF();
-            JacRef = FRef.determinant();
-            STensorOperation::inverseSTensor3(FRef, FRefInv);
-            //FRef *=Fn;
-            //dFRefdFn=..;
-            for (int i=0; i<3; i++){
-                for (int j=0; j<3; j++){
-                    for (int k=0; k<3; k++){
-                        for (int m=0; m<3; m++){
-                            for (int n=0; n<3; n++){
-                                dFRefdFn(i,k,m,n) += qTM1.getRefToReferenceF()(i,j) * I2(j,m) * I2(k,n);
-                            }
-                        }
-                    }
-                }
-            }
         }
         else
         {
-            FRef = Fn;
-            JacRef = FRef.determinant();
-            STensorOperation::inverseSTensor3(FRef, FRefInv);
-            dFRefdFn = I4;
+            FRef = I2;
         }
-
-        if(applyReferenceF)
-        {
-            for(int K = 0; K< 3; K++)
-            {
-                for(int L = 0; L< 3; L++)
-                {
-                    for(int i = 0; i< 3; i++)
-                    {
-                        nutTemp(K,L) += 1.0/(JacRef * vacuum_permeability) * FRefInv(K,i) * FRefInv(i,L) +
-                                            (2.0 * rho)/(_mu_8) * FRefInv(K,i) * FRefInv(i,L);
-                        for(int j = 0; j< 3; j++)
-                        {
-                            // considering push-forward to intermediate config
-                            ltTemp(K,L) += FRef(K,i) * _l0(i,j) * FRef(L,j);
-                            ktTemp(K,L) += FRef(K,i) * k_(i,j) * FRef(L,j);
-                            //--nutTemp(K,L) += FRefInv(i,K) * _nu0(i,j) * FRefInv(j,L); // <<-- Diff transformation
-                            //--mutTemp(K,L) += FRef(K,i) * _mu0(i,j) * FRef(L,j);
+        JacRef = FRef.determinant();
+        STensorOperation::inverseSTensor3(FRef, FRefInv);
+        STensorOperation::transposeSTensor3(FRef,FRefTranspose);
+        STensorOperation::transposeSTensor3(FRefInv,FRefInvTranspose);
+        dFRefdFn = I4;
+
+        STensor3 Ftot(Fn);
+        Ftot*=FRef;
+
+        const double Jactot = Ftot.determinant();
+        static STensor43 dFTotdFn;
+        STensorOperation::zero(dFTotdFn);
+        static STensor3 FtotInv, FtotTranspose, FtotInvTranspose;
+        STensorOperation::zero(FtotInv);
+        STensorOperation::zero(FtotTranspose);
+        STensorOperation::zero(FtotInvTranspose);
+        STensorOperation::inverseSTensor3(Ftot,FtotInv);
+        STensorOperation::transposeSTensor3(Ftot,FtotTranspose);
+        STensorOperation::transposeSTensor3(FtotInv,FtotInvTranspose);
+
+        static STensor3 Ctot, Ctot_squared;
+        STensorOperation::zero(Ctot);
+        STensorOperation::zero(Ctot_squared);
+        STensorOperation::multSTensor3(FtotTranspose,Ftot,Ctot);
+        STensorOperation::multSTensor3(Ctot,Ctot,Ctot_squared);
+
+        static STensor3 FnInv, FnTranspose, FnInvTranspose;
+        STensorOperation::zero(FnInv);
+        STensorOperation::zero(FnTranspose);
+        STensorOperation::zero(FnInvTranspose);
+        STensorOperation::inverseSTensor3(Fn,FnInv);
+        STensorOperation::transposeSTensor3(Fn,FnTranspose);
+        STensorOperation::transposeSTensor3(FnInv,FnInvTranspose);
+
+        for (int i=0; i<3; i++){
+            for (int j=0; j<3; j++){
+                for (int k=0; k<3; k++){
+                    for (int m=0; m<3; m++){
+                        for (int n=0; n<3; n++){
+                            dFTotdFn(i,k,m,n) += I2(i,m) * I2(j,n) * FRef(j,k);
                         }
                     }
                 }
             }
+        }
 
-            // considering push-forward to intermediate config
-            ltTemp *= 1.0/JacRef;
-            ktTemp *= 1.0/JacRef;
-            //--nutTemp *= JacRef; // <<-- Diff transformation
-            //--mutTemp *= 1.0/JacRef;
+        if(applyReferenceF)
+        {
+            // in weak EM problem with Fn = I
+            // d*dFn = 0 always because Fn is constant
+            STensorOperation::zero(dFTotdFn);
         }
 
+        //from HERE NO MORE applyReferenceFq
+
         // Magnetic permeability tensor from constitutive model
         // I_7 = B \otimes B : I, I_8 = B \otimes B : C, I_9 = B \otimes B : C^2
         // W_mag + W_mag-mech = mu_7^{-1} I_7 + mu_8^{-1} I_8 + mu_9^{-1} I_9
-        STensor3 mu_ref(0.0);
-        STensor3 nu_ref(0.0);
-        static STensor3 dnu_refdT;
-        static STensor33 dnu_refdB;
-        STensorOperation::zero(dnu_refdT);
-        STensorOperation::zero(dnu_refdB);
-        STensor3 C(0.0),FRefTranspose(0.0),C_squared(0.0);
-
-        STensorOperation::transposeSTensor3(FRef, FRefTranspose);
-        STensorOperation::multSTensor3(FRefTranspose,FRef,C);
-        STensorOperation::multSTensor3(C,C,C_squared);
-
-        // in case of strong coupling
-        for(unsigned int i = 0; i < 3; ++i)
-            for(unsigned int j = 0; j < 3; ++j)
-            {
-                nu_ref(i,j) = 1.0/(JacRef * vacuum_permeability) * C(i,j)
-                              /*+ (2.0 * rho)/(_mu_7) * I2(i,j)*/
-                              + (2.0 * rho)/(_mu_8) * C(i,j)
-                              /*+ (2.0 * rho)/(_mu_9) * C_squared(i,j)*/;
-            }
-
-        STensorOperation::inverseSTensor3(nu_ref,mu_ref);
-
-        // need to update kt etc...
-        STensor3 FRefInvTranspose;
-        STensorOperation::transposeSTensor3(FRefInv, FRefInvTranspose);
-        const STensor3 dJacRefdF= JacRef * FRefInvTranspose;
-
-        static STensor3 FnInv,FnInvTranspose;
-        STensorOperation::zero(FnInv);
-        STensorOperation::zero(FnInvTranspose);
-        STensorOperation::inverseSTensor3(Fn, FnInv);
-        STensorOperation::transposeSTensor3(FnInv, FnInvTranspose);
-        const double Jacn = Fn.determinant();
-        const STensor3 dJacndF= Jacn * FnInvTranspose;
-
-        // transformation to updated reference configuration
-        for(int K = 0; K< 3; K++)
-        {
-            for(int L = 0; L< 3; L++)
-            {
-                for(int i = 0; i< 3; i++)
-                {
-                    for(int j = 0; j< 3; j++)
-                    {
-                        if(applyReferenceF)
-                        {
-                            // considering pull-back from intermediate to reference config
-                            lt(K,L) += FnInv(K,i) * ltTemp(i,j) * FnInv(L,j);
-                            kt(K,L) += FnInv(K,i) * ktTemp(i,j) * FnInv(L,j);
-                            //--nut(K,L) += Fn(i,K) * nutTemp(i,j) * Fn(j,L); // <<-- Diff transformation
-                            //--mut(K,L) += FnInv(K,i) * mutTemp(i,j) * FnInv(L,j);
-                            nu_ref(K,L) += Fn(i,K) * nutTemp(i,j) * Fn(j,L);
+        for(int K = 0; K< 3; K++) {
+            for (int L = 0; L < 3; L++) {
+                for (int i = 0; i < 3; i++) {
+                    for (int j = 0; j < 3; j++) {
+                        for(int X = 0; X< 3; X++) {
+                            for (int Y = 0; Y < 3; Y++) {
+                                lt(K,L) += FRef(K,X) * FtotInv(X,i) * _l0(i,j) * FtotInv(Y,j) * FRef(L,Y);
+                                kt(K,L) += FRef(K,X) * FtotInv(X,i) * k_(i,j) * FtotInv(Y,j) * FRef(L,Y);
+                            }
                         }
+                        nut(K,L) += FRefInv(i,K) * Ctot(i,j) * (1.0/vacuum_permeability) * FRefInv(j,L);
+                        if(_mu_7 != 0.0)
+                            nut(K,L) += FRefInv(i,K) * Jactot * rho * I2(i,j) * (1.0/_mu_7) * FRefInv(j,L);
+                        if(_mu_8 != 0.0)
+                            nut(K,L) += FRefInv(i,K) * Jactot * rho * Ctot(i,j) * (1.0/_mu_8) * FRefInv(j,L);
+                        if(_mu_9 != 0.0)
+                            nut(K,L) += FRefInv(i,K) * Jactot * rho * Ctot_squared(i,j) * (1.0/_mu_9) * FRefInv(j,L);
+                    }
+                }
+            }
+        }
 
-                        else{
-                        // considering pull-back from current to reference config
-                        lt(K,L) += FRefInv(K,i) * _l0(i,j) * FRefInv(L,j);
-                        kt(K,L) += FRefInv(K,i) * k_(i,j) * FRefInv(L,j);
+        lt *= (Jactot/JacRef);
+        kt *= (Jactot/JacRef);
+        nut *= (JacRef/Jactot);
+        STensorOperation::inverseSTensor3(nut,mut);
+
+        static STensor43 dFTotdFTot,dFTotInvdFTot,dFTotInvTransposedFTot;
+        static STensor43 dCTotdFTot, dCTot_squareddFTot;
+        STensorOperation::zero(dFTotdFTot);
+        STensorOperation::zero(dFTotInvdFTot);
+        STensorOperation::zero(dFTotInvTransposedFTot);
+        STensorOperation::zero(dCTotdFTot);
+        STensorOperation::zero(dCTot_squareddFTot);
+
+        for (int i=0; i<3; i++){
+            for (int j=0; j<3; j++){
+                for (int k=0; k<3; k++){
+                    for(int l=0; l<3; ++l){
+                        dFTotInvdFTot(i,j,k,l) -= FtotInv(i,k) * FtotInv(l,j);
+                        dFTotInvTransposedFTot(i,j,k,l) -= FtotInv(i,l) * FtotInv(k,j);
+                        dCTotdFTot(i,j,k,l) += I2(i,l) * Ftot(k,j);
+                        dCTotdFTot(i,j,k,l) += I2(j,l) * Ftot(k,i);
+                        for(int r=0; r<3; ++r){
+                            dCTot_squareddFTot(i,j,k,l) += I2(i,l) * Ftot(k,r) * Ctot(r,j);
+                            dCTot_squareddFTot(i,j,k,l) += I2(r,l) * Ftot(k,i) * Ctot(r,j);
+                            dCTot_squareddFTot(i,j,k,l) += I2(r,l) * Ftot(k,j) * Ctot(i,r);
+                            dCTot_squareddFTot(i,j,k,l) += I2(j,l) * Ftot(k,r) * Ctot(i,r);
                         }
-
-                        // considering push-forward to intermediate config
-                        /*lt(K,L) += FRef(K,i) * _l0(i,j) * FRef(L,j);
-                        nut(K,L) += FRefInv(i,K) * _nu0(i,j) * FRefInv(j,L); // <<-- Diff transformation
-                        mut(K,L) += FRef(K,i) * _mu0(i,j) * FRef(L,j);
-
-                        kt(K,L) += FRef(K,i) * k_(i,j) * FRef(L,j);*/
                     }
                 }
             }
         }
 
-        if(applyReferenceF)
-        {
-            // considering pull-back from intermediate to reference config
-            lt *= Jacn;
-            kt *= Jacn;
-            //--nut *= 1.0/Jacn; // <<-- Diff transformation
-            //--mut *= Jacn;
-            STensorOperation::inverseSTensor3(nu_ref,mu_ref);
-        }
+        static STensor43 dltdFTot, dktdFTot, dnutdFTot;
+        STensorOperation::zero(dltdFTot);
+        STensorOperation::zero(dktdFTot);
+        STensorOperation::zero(dnutdFTot);
+
+        for(int K = 0; K< 3; K++) {
+            for (int L = 0; L < 3; L++) {
+                for (int i = 0; i < 3; i++) {
+                    for (int j = 0; j < 3; j++) {
+                        dltdFTot(K,L,i,j) += lt(K,L) * FtotInv(j,i);
+                        dktdFTot(K,L,i,j) += kt(K,L) * FtotInv(j,i);
+                        dnutdFTot(K,L,i,j) -= nut(K,L) * FtotInv(j,i);
+                        for(int X = 0; X< 3; X++) {
+                            for (int Y = 0; Y < 3; Y++) {
+                                for(int p = 0; p < 3; ++p){
+                                    for(int Q = 0; Q < 3; ++Q){
+                                    dltdFTot(K,L,p,Q) -= (Jactot/JacRef) * FRef(K,X) * FtotInv(X,p) * FtotInv(Q,i) * _l0(i,j) * FtotInv(Y,j) * FRef(L,Y);
+                                    dltdFTot(K,L,p,Q) -= (Jactot/JacRef) * FRef(K,X) * FtotInv(X,i) * _l0(i,j) * FtotInv(Y,p) * FtotInv(Q,j) * FRef(L,Y);
+
+                                    dktdFTot(K,L,p,Q) -= (Jactot/JacRef) * FRef(K,X) * FtotInv(X,p) * FtotInv(Q,i) * k_(i,j) * FtotInv(Y,j) * FRef(L,Y);
+                                    dktdFTot(K,L,p,Q) -= (Jactot/JacRef) * FRef(K,X) * FtotInv(X,i) * k_(i,j) * FtotInv(Y,p) * FtotInv(Q,j) * FRef(L,Y);
+                                    }
+                                }
 
-        else{
-        // considering pull-back from current to reference config
-        lt *= JacRef;
-        kt *= JacRef;
-        }
+                                dnutdFTot(K,L,X,Y) += (JacRef/Jactot) * FRefInv(i,K) * (1.0/vacuum_permeability) * dCTotdFTot(i,j,X,Y) * FRefInv(j,L);
+                                if(_mu_7 != 0.0)
+                                    dnutdFTot(K,L,X,Y) += (JacRef/Jactot) * FRefInv(i,K) * Jactot * rho * I2(i,j) * FtotInv(Y,X) * (1.0/_mu_7) * FRefInv(j,L);
 
-        // considering push-forward to intermediate config
-        /*lt *= 1.0/JacRef;
-        nut *= JacRef; // <<-- Diff transformation
-        mut *= 1.0/JacRef;
-        kt *= 1.0/JacRef;*/
+                                if(_mu_8 != 0.0) {
+                                    dnutdFTot(K, L, X, Y) += (JacRef / Jactot) * FRefInv(i, K) * Jactot * rho * Ctot(i, j) * FtotInv(Y, X) * (1.0 / _mu_8) * FRefInv(j, L);
+                                    dnutdFTot(K, L, X, Y) += (JacRef / Jactot) * FRefInv(i, K) * Jactot * rho * dCTotdFTot(i, j, X, Y) * (1.0 / _mu_8) * FRefInv(j, L);
+                                }
+                                if(_mu_9 != 0.0) {
+                                    dnutdFTot(K,L,X,Y) += (JacRef/Jactot) * FRefInv(i,K) * Jactot * rho * Ctot_squared(i,j) * FtotInv(Y,X) * (1.0/_mu_9) * FRefInv(j,L);
+                                    dnutdFTot(K,L,X,Y) += (JacRef/Jactot) * FRefInv(i,K) * Jactot * rho * dCTot_squareddFTot(i,j,X,Y) * (1.0/_mu_9) * FRefInv(j,L);
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
 
         // for DG only
         for(int i = 0; i< 3; i++)
@@ -736,15 +743,12 @@ void mlawElecMagGenericThermoMech::constitutive(
                     djydgradTdT(i,j)+=lt(i,j)*_seebeck*_seebeck;
                     djydgradTdV(i,j)+=lt(i,j)*_seebeck;
 
-                    //--dHdT(i) += dnudT(i,j) * B(j);
-                    dHdT(i) += dnu_refdT(i,j) * B(j);
+                    dHdT(i) += dnudT(i,j) * B(j);
 
                     for (unsigned int k = 0; k < 3; ++k)
-                        //--dHdB(i,j) += (d_nu0dB(i,j,k) * B(k));
-                        dHdB(i,j) += (dnu_refdB(i,j,k) * B(k));
+                        dHdB(i,j) += (d_nu0dB(i,j,k) * B(k));
 
-                    //--dHdB(i,j) += nut(i,j);
-                    dHdB(i,j) += nu_ref(i,j);
+                    dHdB(i,j) += nut(i,j);
                 }
             }
         }
@@ -756,8 +760,7 @@ void mlawElecMagGenericThermoMech::constitutive(
             H(i) = 0.0;
             for (unsigned int j = 0; j < 3; ++j)
             {
-                //--H(i) += nut(i,j) * B(j); // old nu empirical model
-                H(i) += nu_ref(i,j) * B(j); // New nu using constitutive model
+                H(i) += nut(i,j) * B(j);
             }
         }
 
@@ -902,30 +905,16 @@ void mlawElecMagGenericThermoMech::constitutive(
             }
         }
 
-        // Double Pull-back transformation in weak EM problem
-        // first: from current to intermediate with Fn (=I possibly)
-        // second: from intermediate to reference config with FRef
-        // to correctly employ it in weak TM problem
-        if(applyReferenceF)
-        {
-            ThermalEMFieldSource *= JacRef * Jacn;
-            dThermalEMFieldSourcedA *= JacRef * Jacn;
-            dThermalEMFieldSourcedB *= JacRef * Jacn;
-            dThermalEMFieldSourcedGradT *= JacRef * Jacn;
-            dThermalEMFieldSourcedGradV *= JacRef * Jacn;
-            dThermalEMFieldSourcedT *= JacRef * Jacn;
-            dThermalEMFieldSourcedV *= JacRef * Jacn;
-        }
-        else // Pull-back transformation in strong coupling
-        {
-            ThermalEMFieldSource *= JacRef;
-            dThermalEMFieldSourcedA *= JacRef;
-            dThermalEMFieldSourcedB *= JacRef;
-            dThermalEMFieldSourcedGradT *= JacRef;
-            dThermalEMFieldSourcedGradV *= JacRef;
-            dThermalEMFieldSourcedT *= JacRef;
-            dThermalEMFieldSourcedV *= JacRef;
-        }
+        // to use in weak TM with transformation
+        // because em source computed in intermediate config
+        // but needs to be used in reference config in weak TM
+        ThermalEMFieldSource *= JacRef;
+        dThermalEMFieldSourcedA *= JacRef;
+        dThermalEMFieldSourcedB *= JacRef;
+        dThermalEMFieldSourcedGradT *= JacRef;
+        dThermalEMFieldSourcedGradV *= JacRef;
+        dThermalEMFieldSourcedT *= JacRef;
+        dThermalEMFieldSourcedV *= JacRef;
 
         // derivatives with deformation gradient
         if(stiff)
@@ -936,100 +925,59 @@ void mlawElecMagGenericThermoMech::constitutive(
                 {
                     for(int N = 0; N< 3; N++)
                     {
-                        for(int P = 0; P < 3; ++P)
+                        for(int L = 0; L< 3; L++)
                         {
-                            for(int Q = 0; Q < 3; ++Q)
+                            //dHdF(K,m,N) += nut(K,L) * dBdF(L,m,N);
+                            for(int p = 0; p < 3; ++p)
                             {
-                                for(int L = 0; L< 3; L++)
+                                for(int Q = 0; Q < 3; ++Q)
                                 {
-                                    // missing mult with dFrefdFn
-                                    dedF(K,m,N) -= FRefInv(K,Q) * lt(P,L) * dFRefdFn(Q,P,m,N) * gradV(L);
-                                    dedF(K,m,N) -= FRefInv(L,Q) * lt(K,P) * dFRefdFn(Q,P,m,N) * gradV(L);
-                                    dedF(K,m,N) += FRefInv(P,Q) * lt(K,L) * dFRefdFn(Q,P,m,N) * gradV(L);
-                                    dedF(K,m,N) -= FRefInv(K,Q) * lt(P,L) * dFRefdFn(Q,P,m,N) * _seebeck * gradT(L);
-                                    dedF(K,m,N) -= FRefInv(L,Q) * lt(K,P) * dFRefdFn(Q,P,m,N) * _seebeck * gradT(L);
-                                    dedF(K,m,N) += FRefInv(P,Q) * lt(K,L) * dFRefdFn(Q,P,m,N) * _seebeck * gradT(L);
-                                    dedF(K,m,N) -= FRefInv(K,Q) * lt(P,L) * dFRefdFn(Q,P,m,N) * ((An(L) - A0(L)) * invdh);
-                                    dedF(K,m,N) -= FRefInv(L,Q) * lt(K,P) * dFRefdFn(Q,P,m,N) * ((An(L) - A0(L)) * invdh);
-                                    dedF(K,m,N) += FRefInv(P,Q) * lt(K,L) * dFRefdFn(Q,P,m,N) * ((An(L) - A0(L)) * invdh);
-
-                                    dqdF(K,m,N) -= FRefInv(K,Q) * lt(P,L) * dFRefdFn(Q,P,m,N) * _seebeck * T * gradV(L);
-                                    dqdF(K,m,N) -= FRefInv(L,Q) * lt(K,P) * dFRefdFn(Q,P,m,N) * _seebeck * T * gradV(L);
-                                    dqdF(K,m,N) += FRefInv(P,Q) * lt(K,L) * dFRefdFn(Q,P,m,N) * _seebeck * T * gradV(L);
-                                    dqdF(K,m,N) -= FRefInv(K,Q) * lt(P,L) * dFRefdFn(Q,P,m,N) * _seebeck * _seebeck * T * gradT(L);
-                                    dqdF(K,m,N) -= FRefInv(L,Q) * lt(K,P) * dFRefdFn(Q,P,m,N) * _seebeck * _seebeck * T * gradT(L);
-                                    dqdF(K,m,N) += FRefInv(P,Q) * lt(K,L) * dFRefdFn(Q,P,m,N) * _seebeck * _seebeck * T * gradT(L);
-                                    dqdF(K,m,N) -= FRefInv(K,Q) * lt(P,L) * dFRefdFn(Q,P,m,N) * _seebeck * T * ((An(L) - A0(L)) * invdh);
-                                    dqdF(K,m,N) -= FRefInv(L,Q) * lt(K,P) * dFRefdFn(Q,P,m,N) * _seebeck * T * ((An(L) - A0(L)) * invdh);
-                                    dqdF(K,m,N) += FRefInv(P,Q) * lt(K,L) * dFRefdFn(Q,P,m,N) * _seebeck * T * ((An(L) - A0(L)) * invdh);
-
-                                    // need to check this after using nu_ref
-                                    dHdF(K,m,N) -= FRefInv(P,Q) * nut(K,L) * dFRefdFn(Q,P,m,N) * B(L);
-
-                                    for(unsigned int a = 0; a < 3; ++a)
-                                        for(unsigned int b = 0; b < 3; ++b)
-                                        {
-                                            dHdF(K,m,N) += (1.0/Jacn) * I2(a,P)* I2(K,Q) * _nu0(a,b) * FRef(b,L) * dFRefdFn(P,Q,m,N) * B(L);
-                                            dHdF(K,m,N) += (1.0/Jacn) * FRef(a,K) * _nu0(a,b) * I2(b,P)* I2(L,Q) * dFRefdFn(P,Q,m,N) * B(L);
-                                        }
-
-                                    dsourceVectorFielddF(K,m,N) -= -1.0 * FRefInv(K,Q) * lt(P,L) * dFRefdFn(Q,P,m,N) * gradV(L);
-                                    dsourceVectorFielddF(K,m,N) -= -1.0 * FRefInv(L,Q) * lt(K,P) * dFRefdFn(Q,P,m,N) * gradV(L);
-                                    dsourceVectorFielddF(K,m,N) += -1.0 * FRefInv(P,Q) * lt(K,L) * dFRefdFn(Q,P,m,N) * gradV(L);
-                                    dsourceVectorFielddF(K,m,N) -= -1.0 * FRefInv(K,Q) * lt(P,L) * dFRefdFn(Q,P,m,N) * _seebeck * gradT(L);
-                                    dsourceVectorFielddF(K,m,N) -= -1.0 * FRefInv(L,Q) * lt(K,P) * dFRefdFn(Q,P,m,N) * _seebeck * gradT(L);
-                                    dsourceVectorFielddF(K,m,N) += -1.0 * FRefInv(P,Q) * lt(K,L) * dFRefdFn(Q,P,m,N) * _seebeck * gradT(L);
-                                    dsourceVectorFielddF(K,m,N) -= -1.0 * FRefInv(K,Q) * lt(P,L) * dFRefdFn(Q,P,m,N) * ((An(L) - A0(L)) * invdh);
-                                    dsourceVectorFielddF(K,m,N) -= -1.0 * FRefInv(L,Q) * lt(K,P) * dFRefdFn(Q,P,m,N) * ((An(L) - A0(L)) * invdh);
-                                    dsourceVectorFielddF(K,m,N) += -1.0 * FRefInv(P,Q) * lt(K,L) * dFRefdFn(Q,P,m,N) * ((An(L) - A0(L)) * invdh);
+                                    dedF(K,m,N) += dltdFTot(K,L,p,Q) * dFTotdFn(p,Q,m,N) * gradV(L);
+                                    dedF(K,m,N) += dltdFTot(K,L,p,Q) * dFTotdFn(p,Q,m,N) * _alpha * gradT(L);
+                                    dedF(K,m,N) += dltdFTot(K,L,p,Q)* dFTotdFn(p,Q,m,N) * ((An(L) - A0(L)) * invdh);
+
+                                    dqdF(K,m,N) += dltdFTot(K,L,p,Q) * dFTotdFn(p,Q,m,N) * _seebeck * _seebeck * T * gradT(L);
+                                    dqdF(K,m,N) += dltdFTot(K,L,p,Q) * dFTotdFn(p,Q,m,N) * _seebeck * T * gradV(L);
+                                    dqdF(K,m,N) += dltdFTot(K,L,p,Q)* dFTotdFn(p,Q,m,N) * _seebeck * T * ((An(L) - A0(L)) * invdh);
+                                    dqdF(K,m,N) += dktdFTot(K,L,p,Q) * dFTotdFn(p,Q,m,N) * gradT(L);
+
+                                    djydF(K,m,N) += dktdFTot(K,L,p,Q) * dFTotdFn(p,Q,m,N) * gradT(L);
+                                    djydF(K,m,N) += dltdFTot(K,L,p,Q) * dFTotdFn(p,Q,m,N) * _seebeck * _seebeck * T * gradT(L);
+                                    djydF(K,m,N) += dltdFTot(K,L,p,Q) * dFTotdFn(p,Q,m,N) * _seebeck * V * gradT(L);
+                                    djydF(K,m,N) += dltdFTot(K,L,p,Q) * dFTotdFn(p,Q,m,N) * V * gradV(L);
+                                    djydF(K,m,N) += dltdFTot(K,L,p,Q) * dFTotdFn(p,Q,m,N) * _seebeck * T * gradV(L);
+                                    djydF(K,m,N) += dltdFTot(K,L,p,Q) * dFTotdFn(p,Q,m,N) * V * ((An(L) - A0(L)) * invdh);
+                                    djydF(K,m,N) += dltdFTot(K,L,p,Q) * dFTotdFn(p,Q,m,N) * _seebeck * T * ((An(L) - A0(L)) * invdh);
+
+                                    dHdF(K,m,N) += dnutdFTot(K,L,p,Q) * dFTotdFn(p,Q,m,N) * B(L);
+
+                                    dsourceVectorFielddF(K,m,N) += -1.0 * dltdFTot(K,L,p,Q) * dFTotdFn(p,Q,m,N) * gradV(L);
+                                    dsourceVectorFielddF(K,m,N) += -1.0 * dltdFTot(K,L,p,Q) * dFTotdFn(p,Q,m,N) * _alpha * gradT(L);
+                                    dsourceVectorFielddF(K,m,N) += -1.0 * dltdFTot(K,L,p,Q)* dFTotdFn(p,Q,m,N) * ((An(L) - A0(L)) * invdh);
                                 }
                             }
                         }
+                    }
+                }
+            }
 
-                        // need to check index
-                        // double Pull-back transformation
-                        // first: weak EM problem with Fn (=I possibly) from current to intermediate config
-                        // second: weak EM in intermediate to weak TM in reference config with FRef
-                        if(applyReferenceF)
-                        {
-                            dThermalEMFieldSourcedF(m,N) += JacRef * Jacn * dsourceVectorFielddF(K, m,N) * ((An(K) - A0(K)) * invdh);
-
-                            if (_useFluxT)
-                                dThermalEMFieldSourcedF(m,N) += JacRef * Jacn * dsourceVectorFielddF(K, m,N) * gradV(K);
-
-                            // Magnetic hysteresis loss contribution
-                            /*dThermalEMFieldSourcedF(m,N) += dHdF(K,m,N) * (Bn(K) - B0(K)) * invdh;
-                            dThermalEMFieldSourcedF(m,N) += H(K) * dBdF(K,m,N) * invdh;*/
-
-                            // Missing terms of dJacdF for ThermalEMFieldSource to
-                            // correctly employ it in TM problem
-                            dThermalEMFieldSourcedF(m,N) += JacRef * dJacndF(m,N) * sourceVectorField(K) * ((An(K) - A0(K)) * invdh);
-                            dThermalEMFieldSourcedF(m,N) += JacRef * dJacndF(m,N) * sourceVectorField(K) * gradV(K);
-                            dThermalEMFieldSourcedF(m,N) += Jacn * dJacRefdF(m,N) * sourceVectorField(K) * ((An(K) - A0(K)) * invdh);
-                            dThermalEMFieldSourcedF(m,N) += Jacn * dJacRefdF(m,N) * sourceVectorField(K) * gradV(K);
-                        }
-                        else //strong coupling case
-                        {
-                            dThermalEMFieldSourcedF(m,N) += JacRef * dsourceVectorFielddF(K, m,N) * ((An(K) - A0(K)) * invdh);
-
-                            if (_useFluxT)
-                                dThermalEMFieldSourcedF(m,N) += JacRef * dsourceVectorFielddF(K, m,N) * gradV(K);
-
-                            // Magnetic hysteresis loss contribution
-                            /*dThermalEMFieldSourcedF(m,N) += dHdF(K,m,N) * (Bn(K) - B0(K)) * invdh;
-                            dThermalEMFieldSourcedF(m,N) += H(K) * dBdF(K,m,N) * invdh;*/
+            for(int K = 0; K< 3; K++)
+            {
+                for(int m = 0; m< 3; m++)
+                {
+                    for(int N = 0; N< 3; N++)
+                    {
+                        dThermalEMFieldSourcedF(m,N) += JacRef * dsourceVectorFielddF(K,m,N) * ((An(K) - A0(K)) * invdh);
 
-                            // Missing terms of dJacdF for ThermalEMFieldSource to
-                            // correctly employ it in TM problem
-                            dThermalEMFieldSourcedF(m,N) += dJacRefdF(m,N) * sourceVectorField(K) * ((An(K) - A0(K)) * invdh);
-                            dThermalEMFieldSourcedF(m,N) += dJacRefdF(m,N) * sourceVectorField(K) * gradV(K);
-                        }
+                        if (_useFluxT)
+                            dThermalEMFieldSourcedF(m,N) += JacRef * dsourceVectorFielddF(K,m,N) * gradV(K);
 
+                        // Magnetic hysteresis loss contribution
+                        /*dThermalEMFieldSourcedF(m,N) += JacRef * dHdF(K,m,N) * (Bn(K) - B0(K)) * invdh;
+                        dThermalEMFieldSourcedF(m,N) += JacRef * H(K) * dBdF(K,m,N) * invdh;*/
                     }
                 }
             }
-
         }
-
     }
 }
\ No newline at end of file
diff --git a/NonLinearSolver/materialLaw/mlawGenericTM.cpp b/NonLinearSolver/materialLaw/mlawGenericTM.cpp
index c55c5b85798f1f78bc83f872ef55aead5506b955..b360013dc7dc6a07f0d172eaf5bfaa2360dc619e 100644
--- a/NonLinearSolver/materialLaw/mlawGenericTM.cpp
+++ b/NonLinearSolver/materialLaw/mlawGenericTM.cpp
@@ -258,15 +258,14 @@ void mlawGenericTM::constitutive(const STensor3   &F0,
 
     if(applyReferenceF)
     {
+       //FRef = Fn; // here we have Fn identity in any case but I would leave it to be correct
        FRef = qTM1.getRefToReferenceF();
-       FRef *=Fn;
-       //dFRefdFn=..;
         for (int i=0; i<3; i++){
             for (int j=0; j<3; j++){
                 for (int k=0; k<3; k++){
                     for (int m=0; m<3; m++){
                         for (int n=0; n<3; n++){
-                            dFRefdFn(i,k,m,n) += qTM1.getRefToReferenceF()(i,j) * I2(j,m) * I2(k,n);
+                            dFRefdFn(i,k,m,n) += qTM1.getRefToReferenceF()(i,j) * I2(j,m) * I2(k,n); /// to chnage to reflect order
                         }
                     }
                 }
@@ -276,7 +275,6 @@ void mlawGenericTM::constitutive(const STensor3   &F0,
     else
     {
         FRef = Fn;
-        //STensorOperation::unity(dFRefdFn);
         dFRefdFn = I4;
     }
     // need to updqte kt etc...