Ludovic Noels · Ludovic Noels · 4e6928a2 · 13bec0c9 · 13bec0c9 · 13bec0c9
--- a/NonLinearSolver/internalPoints/ipvariable.h
+++ b/NonLinearSolver/internalPoints/ipvariable.h
    return *this;
  };
  virtual double get(const int i) const{return 0.;} // Allow to get any component defined in the ipvariable
+  virtual double & getRef(const int i){ static double temp = 0.0; return temp; }; // Allow to get access to any component value defined in the ipvariable
--- a/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
+++ b/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
+                                avgValue += IPVariableOverTime.at(t).at(eleNum).at(j);
+                            }
+                        }
+                        avgValue /= count;
+                        tempMap.emplace(j, avgValue);
+                    }
+                    avgIPVariable.emplace(eleNum, tempMap);
+                }
+            }
+        }
+    }
+    Msg::Info("Done computing average IPVariable over time period");
+}
+
+void nonLinearMechSolver::rewriteIPDataOnPhysical(const std::map< std::size_t, std::map< int, double > > & avgIPVariable)
+{
--- a/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
+++ b/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
+
+            for (int idom = 0; idom < domainVector.size(); idom++)
+            {
+                partDomain *dom = domainVector[idom];
+                if (dom->g_find(ele))
+                {
+                    IntPt *GP;
+                    const int npts = dom->getBulkGaussIntegrationRule()->getIntPoints(ele, &GP);
+
+                    // loop over each Gauss point
+                    for (int j = 0; j < npts; j++)
+                    {
+                        // here can check GP against std::map coordinates for correctness
+
+                        double & ThermalEMFieldSource = ips[j]->getState(IPStateBase::current)->getRef(IPField::Output(ipval));
+                        ThermalEMFieldSource = avgIPVariable.at(eleNum).at(j);
--- a/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
+++ b/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
+            for (int idom = 0; idom < domainVector.size(); idom++)
+            {
+                partDomain *dom = domainVector[idom];
+                if (dom->g_find(ele))
+                {
+                    IntPt *GP;
+                    const int npts = dom->getBulkGaussIntegrationRule()->getIntPoints(ele, &GP);
+                    std::map< int, double > tempMap;
+
+                    // loop over each Gauss point
+                    for (int j = 0; j < npts; j++)
+                    {
+                        double avgValue = 0.0;
+                        int count = 0;
+
+                        // here can check GP against std::map coordinates for correctness
--- a/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
+++ b/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
                            if(t >= tstart && t <= tend)
                            {
                                count++;
+                                sum_dt += t;

                                // compute average of IPVariable over given time period
-                                avgValue += IPVariableOverTime.at(t).at(eleNum).at(j);
+                                //avgValue += IPVariableOverTime.at(t).at(eleNum).at(j);
+
+                                // sum(w_AV * dt)
+                                avgValue += (IPVariableOverTime.at(t).at(eleNum).at(j)) * t;
                            }
                        }
-                        avgValue /= count;
+                        //avgValue /= count;
+                        avgValue /= sum_dt;
--- a/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
+++ b/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
                        {
                            const double t = it->first;
+                            t_i = it->first;
                            if(t >= tstart && t <= tend)
                            {
                                count++;
-                                sum_dt += t;
+                                dt = (t_i - t_0);
+                                sum_dt += dt;

                                // compute average of IPVariable over given time period
                                //avgValue += IPVariableOverTime.at(t).at(eleNum).at(j);

                                // sum(w_AV * dt)
-                                avgValue += (IPVariableOverTime.at(t).at(eleNum).at(j)) * t;
+                                avgValue += (IPVariableOverTime.at(t).at(eleNum).at(j)) * dt;
--- a/dG3D/benchmarks/strongCouplingSMP/StrongCouplingEMTM.py
+++ b/dG3D/benchmarks/strongCouplingSMP/StrongCouplingEMTM.py
 # material law: vacuum/free space region
 lawnumvac = 3
 rhovac = 1.2 
-Gvac=156.e0 # Shear modulus
+Gvac=156.e1 # Shear modulus
 Muvac = 0.35 # Poisson ratio
 Evac= 2.0 * Gvac * (1.0 + Muvac) #youngs modulus
 alphavac = betavac = gammavac = 0. # parameter of anisotropy
-cpvac= 1012.*rhovac
-Kxvac=Kyvac=Kzvac= 26.e-3 # thermal conductivity tensor components
+cpvac= 1.e-12 #1012.*rhovac
--- a/dG3D/benchmarks/weakCouplingSMP/EM_weakCoupling.py
+++ b/dG3D/benchmarks/weakCouplingSMP/EM_weakCoupling.py
 Muvac = 0.35 # Poisson ratio (same as SMP)
 Evac= 2.0 * Gvac * (1.0 + Muvac) #youngs modulus
 alphavac = betavac = gammavac = 0. # parameter of anisotropy
-cpvac= 1012.*rhovac
-Kxvac=Kyvac=Kzvac= 26.e-3 # thermal conductivity tensor components
+cpvac= 0.0 #1012.*rhovac
--- a/dG3D/benchmarks/strongCouplingSMP/StrongCouplingEMTM.py
+++ b/dG3D/benchmarks/strongCouplingSMP/StrongCouplingEMTM.py
 sol = 2  # Gmm=0 (default) Taucs=1 PETsc=2
 soltype = 1 # StaticLinear=0 (default) StaticNonLinear=1
 #EMncycles = 5; # num of sinusoidal cycles of EM frequency application
-nstep = 80 # 8*EMncycles # number of step (used only if soltype=1)
+nstep = 200 # 8*EMncycles # number of step (used only if soltype=1)
 ftime = 5.e-3 # EMncycles/freq; # according to characteristic time
 #EMftime = EMncycles/freq # Final time (used only if soltype=1)
 tol=1.e-6   # relative tolerance for NR scheme (used only if soltype=1)
 nstepArch=1 # Number of step between 2 archiving (used only if soltype=1)
 fullDg = 0 #O = CG, 1 = DG
 beta1 = 1000.
-eqRatio =1.e8
+eqRatio =1.e0
--- a/dG3D/benchmarks/weakCouplingSMP/new_merge_mesh.geo 0 → 100644
+++ b/dG3D/benchmarks/weakCouplingSMP/new_merge_mesh.geo 0 → 100644
+  Rint = {200*cm, Min 0.15, Max 0.9, Step 0.1, Name StrCat[pp2, "7Inner [m]"],
+    Visible (Flag_Infinity==1), Highlight Str[colorpp] },
+  Rext = {280*cm, Min Rint, Max 1, Step 0.1, Name StrCat[pp2, StrCat["8", label_Rext]],
+    Visible 1, Highlight Str[colorpp] }
+];
+/*
+lc0 = Pi*Rint/6; // air
+Characteristic Length { PointsOf{ Volume{v4}; } } = lc0;
+
+lc1  = 2*wcoil/2; // coil
+Characteristic Length { PointsOf{ Volume{e(1)}; } } = lc1;
+*/
+//Mesh 2;
+//Mesh 3;
+//Mesh.SaveAll = 1;
+//Save "mergedMesh.msh";
--- a/NonLinearSolver/materialLaw/mlawLinearElecMagTherMech.cpp
+++ b/NonLinearSolver/materialLaw/mlawLinearElecMagTherMech.cpp
            ThermalEMFieldSource += (sourceVectorField(i) * ((An(i) - A0(i)) * invdh) /*+ H(i) * (Bn(i) - B0(i)) * invdh*/);

+            // to debug weak vs strong coupling results by
+            // using analytical EM thermal source in SMP only with
+            // A²/2 being mean value that is calculated
+            // in weak coupling for considered parameter values
+
+            /*
+              if(this->_num == 1) // if SMP
+                // A² cos²(2 pi f t)
+                // A²/2 = Mean value over one period = 3559017.23076498
+                ThermalEMFieldSource = 2.0 * 2227580. * std::cos(2.0 * M_PI * _freqEM * ctime) *
+                                        std::cos(2.0 * M_PI * _freqEM * ctime);
+            else
+                ThermalEMFieldSource = 0.0;
+            */
--- a/NonLinearSolver/materialLaw/mlawLinearElecMagTherMech.cpp
+++ b/NonLinearSolver/materialLaw/mlawLinearElecMagTherMech.cpp
    STensorOperation::zero(dmechFieldSourcedgradV);*/

    }
+/*
--- a/NonLinearSolver/materialLaw/mlawLinearElecTherMech.cpp
+++ b/NonLinearSolver/materialLaw/mlawLinearElecTherMech.cpp
        STensorOperation::zero(dqdF);

+        STensor3 Ft;
+        STensorOperation::unity(Ft);
+        //if (!evaluateCurlField)
+            //Ft = Fn;
+
        for(int i=0;i<3;i++)
        {
            fluxT(i)=0.;
            for(int j=0;j<3;j++)
+            //for(int k = 0; k < 3; ++k)
            {
-                fluxT(i)+=(_k(i,j)+_l0(i,j)*_seebeck*_seebeck*T)*gradT(j)+_l0(i,j)*_seebeck*T*gradV(j);
-                fluxjy(i)+=(_k(i,j)+_l0(i,j)*_seebeck*_seebeck*T+_l0(i,j)*_seebeck*V)*gradT(j)+(_l0(i,j)*V+_l0(i,j)*_seebeck*T)*gradV(j);
+                fluxT(i)+=(_k(i,j)+_l0(i,j)*_seebeck*_seebeck*T)*gradT(j)+_l0(i,j)*_seebeck*T/**Ft(k,j)*/*gradV(j);
--- a/dG3D/benchmarks/ThermoMechanicsInSMP/ThermoMechanicsInSMP.py 0 → 100644
+++ b/dG3D/benchmarks/ThermoMechanicsInSMP/ThermoMechanicsInSMP.py 0 → 100644
+#coding-Utf-8-*-
+from gmshpy import *
+
+
+#from dG3DpyDebug import*
+from dG3Dpy import*
+
+#script to launch unit test to debug strong vs weak coupling in SMP
+
+# Here strong coupled EM-TM in SMP only with analytically computed EM thermal source
+# No solution of EM fields but directly EM thermal source is imported to 
+# resolve TM in SMP
+
+# Need to uncomment in mlawLinearEMTM
+# wAV = A² cos² (2 pi f t)
+# with A²/2 = mean computed over one period for considered f, I from weak coupling test
--- a/NonLinearSolver/materialLaw/mlawElecMagGenericThermoMech.cpp
+++ b/NonLinearSolver/materialLaw/mlawElecMagGenericThermoMech.cpp
    const IPGenericTM & qTM0 = qEM0.getConstRefToIpThermoMech();
    IPGenericTM & qTM1 = qEM1.getRefToIpThermoMech();

+    if(_evaluateTemperature)
--- a/NonLinearSolver/materialLaw/mlawElecMagGenericThermoMech.cpp
+++ b/NonLinearSolver/materialLaw/mlawElecMagGenericThermoMech.cpp
+        // derivatives with deformation gradient
+        if(stiff)
+        {
+            for(int K = 0; K< 3; K++)
+            {
+                for(int m = 0; m< 3; m++)
+                {
+                    for(int N = 0; N< 3; N++)
+                    {
+                        for(int P = 0; P < 3; ++P)
+                        {
+                            for(int Q = 0; Q < 3; ++Q)
+                            {
+                                for(int L = 0; L< 3; L++)
+                                {
+                                    // missing mult with dFrefdFn
--- a/NonLinearSolver/materialLaw/mlawElecMagGenericThermoMech.cpp
+++ b/NonLinearSolver/materialLaw/mlawElecMagGenericThermoMech.cpp
+                                    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
+                                    dHdF(K,m,N) -= FRefInv(P,Q) * nut(K,L) * dFRefdFn(Q,P,m,N) * B(L);
--- a/NonLinearSolver/materialLaw/mlawElecMagGenericThermoMech.cpp
+++ b/NonLinearSolver/materialLaw/mlawElecMagGenericThermoMech.cpp
                {
                    for(int j = 0; j< 3; j++)
                    {
-                        lt(K,L) += Fn(K,i) * _l0(i,j) * Fn(L,j);
-                        nut(K,L) += Fninv(i,K) * _nu0(i,j) * Fninv(j,L); // <<-- Diff transformation
-                        mut(K,L) += Fn(K,i) * _mu0(i,j) * Fn(L,j);
+                        // considering pull-back to reference config
+                        lt(K,L) += FRefInv(K,i) * _l0(i,j) * FRefInv(L,j);
+                        nut(K,L) += FRef(i,K) * _nu0(i,j) * FRef(j,L); // <<-- Diff transformation
+                        mut(K,L) += FRefInv(K,i) * _mu0(i,j) * FRefInv(L,j);
+
+                        kt(K,L) += FRefInv(K,i) * k_(i,j) * FRefInv(L,j);
+
+                        // considering push-forward to intermediate config
--- a/dG3D/benchmarks/weakCouplingSMP/run.py
+++ b/dG3D/benchmarks/weakCouplingSMP/run.py
 os.system('mkdir '+EMresultsDir)

 initialDisp = 0.0;
-deltaDisp = 0.01;
+deltaDisp = 0.004;
--- a/dG3D/benchmarks/weakCouplingSMP/new_merge_mesh.geo
+++ b/dG3D/benchmarks/weakCouplingSMP/new_merge_mesh.geo
 Printf("surf = %g", e(9));
 */

-Transfinite Surface {e(2),e(3),e(4),e(5),e(6),e(7),e(8),e(9)};
+//#Transfinite Surface {e(2),e(3),e(4),e(5),e(6),e(7),e(8),e(9)};
--- a/dG3D/benchmarks/weakCouplingSMP/readAndPlotCSV.py 0 → 100644
+++ b/dG3D/benchmarks/weakCouplingSMP/readAndPlotCSV.py 0 → 100644
+import pandas as pd
--- a/dG3D/benchmarks/LargeDeformations_WeakCouplingEMGenericTM_SMP/EM_weakCoupling.py
+++ b/dG3D/benchmarks/LargeDeformations_WeakCouplingEMGenericTM_SMP/EM_weakCoupling.py
 # material law: vacuum/free space region
 lawnumvac = 7
 rhovac = 1.2 
-Gvac=156.e1 # Shear modulus (1.e-5 less than SMP)
-Muvac = 0.35 # Poisson ratio (same as SMP)
+Gvac=156.e3 # Shear modulus (1.e-5 less than SMP)
+Muvac = 0.45 # Poisson ratio (same as SMP)
--- a/dG3D/benchmarks/LargeDeformations_WeakCouplingEMGenericTM_SMP/EM_weakCoupling.py
+++ b/dG3D/benchmarks/LargeDeformations_WeakCouplingEMGenericTM_SMP/EM_weakCoupling.py

 # ElecMag Generic Thermo-Mech law for smp
 lawsmp = ElecMagGenericThermoMechanicsDG3DMaterialLaw(lawnum+2,rho,alpha,beta,gamma,lx,ly,lz,seebeck,v0,mag_mu_x,
-mag_mu_y,mag_mu_z,magpot0_x,magpot0_y, magpot0_z, Irms, freq, nTurnsCoil, coilLx, coilLy, coilLz, coilW,useFluxT,evaluateCurlField)
+mag_mu_y,mag_mu_z,magpot0_x,magpot0_y, magpot0_z, Irms, freq, nTurnsCoil, coilLx, coilLy, coilLz, coilW,useFluxT,evaluateCurlField,evaluateTemperature)

 lawsmp.setThermoMechanicalMaterialLaw(lawTM);
+lawsmp.setApplyReferenceF(False);
--- a/dG3D/benchmarks/LargeDeformations_WeakCouplingEMGenericTM_SMP/TM_weakCoupling.py
+++ b/dG3D/benchmarks/LargeDeformations_WeakCouplingEMGenericTM_SMP/TM_weakCoupling.py

 # material law for smp
 # Mechanical law
-sy0=1e10
+sy0=1e16
 h=1.e9
 harden = LinearExponentialJ2IsotropicHardening(1, sy0, h, 0., 10.)
 lawMech   = J2LinearDG3DMaterialLaw(lawnum,rho,E_matrix,Mu,harden)
+lawMech.setStrainOrder(-1)

 # Generic Thermo-Mech law
 lawTM   = GenericThermoMechanicsDG3DMaterialLaw(lawnum+1,rho,alpha,beta,gamma,tinitial,Kx,Ky,Kz,alphaTherm,alphaTherm,alphaTherm,cp)
 lawTM.setMechanicalMaterialLaw(lawMech);
-lawTM.setApplyReferenceF(True);
+lawTM.setApplyReferenceF(False);
--- a/dG3D/benchmarks/LargeDeformations_WeakCouplingEMGenericTM_SMP/new_merge_mesh.geo
+++ b/dG3D/benchmarks/LargeDeformations_WeakCouplingEMGenericTM_SMP/new_merge_mesh.geo
 Surface(s1) = {cl1, cl2};
 //+

-Transfinite Curve {c1,c2,c3,c4} = 11 Using Progression 1;
--- a/dG3D/benchmarks/LargeDeformations_WeakCouplingEMGenericTM_SMP/EM_weakCoupling.py
+++ b/dG3D/benchmarks/LargeDeformations_WeakCouplingEMGenericTM_SMP/EM_weakCoupling.py
 mysolver.addDomain(SMP_field)
 mysolver.addDomain(Inductor_field)
 mysolver.addDomain(Vacuum_field)
-mysolver.addMaterialLaw(lawsmp)
+mysolver.addMaterialLaw(pertLawsmp)
 mysolver.addMaterialLaw(lawind)
-mysolver.addMaterialLaw(lawvac)
+mysolver.addMaterialLaw(pertLawvac)
--- a/dG3D/benchmarks/LargeDeformations_WeakCouplingEMGenericTM_SMP/TM_weakCoupling.py
+++ b/dG3D/benchmarks/LargeDeformations_WeakCouplingEMGenericTM_SMP/TM_weakCoupling.py
 mysolver = nonLinearMechSolver(1000)
 mysolver.loadModel(meshfile)
 mysolver.addDomain(SMP_field)
-mysolver.addMaterialLaw(lawsmp)
+mysolver.addMaterialLaw(pertLawsmp)
--- a/dG3D/benchmarks/strongCouplingEMGenericTM_SMP/StrongCouplingEMGenericTM.py
+++ b/dG3D/benchmarks/strongCouplingEMGenericTM_SMP/StrongCouplingEMGenericTM.py
 mysolver.addDomain(SMP_field)
 mysolver.addDomain(Inductor_field)
 mysolver.addDomain(Vacuum_field)
-mysolver.addMaterialLaw(lawsmp)
+mysolver.addMaterialLaw(pertLawsmp)
 mysolver.addMaterialLaw(lawind)
-mysolver.addMaterialLaw(lawvac)
+mysolver.addMaterialLaw(pertLawvac)
--- a/NonLinearSolver/materialLaw/mlawElecMagGenericThermoMech.cpp
+++ b/NonLinearSolver/materialLaw/mlawElecMagGenericThermoMech.cpp
                {
                    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);
--- a/NonLinearSolver/materialLaw/mlawElecMagGenericThermoMech.cpp
+++ b/NonLinearSolver/materialLaw/mlawElecMagGenericThermoMech.cpp
        STensorOperation::multSTensor3(FRefTranspose,FRef,C);
        STensorOperation::multSTensor3(C,C,C_squared);
-        const double vacuum_permeability = 4.0 * M_PI * 1.e-7;

+        // 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/(2.0 * JacRef * vacuum_permeability) * C(i,j)
-                              - (rho)/(_mu_7) * I2(i,j)
-                              - (rho)/(_mu_8) * C(i,j)
-                              - (rho)/(_mu_9) * C_squared(i,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);
--- a/NonLinearSolver/materialLaw/mlawElecMagGenericThermoMech.cpp
+++ b/NonLinearSolver/materialLaw/mlawElecMagGenericThermoMech.cpp
                        {
                            // considering pull-back from intermediate to reference config
                            lt(K,L) += FnInv(K,i) * ltTemp(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);
-
                            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);
--- a/NonLinearSolver/materialLaw/mlawElecMagGenericThermoMech.cpp
+++ b/NonLinearSolver/materialLaw/mlawElecMagGenericThermoMech.cpp
                                    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
+                                    // 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);