Ludovic Noels · 4e6928a2 · 13bec0c9 · 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