Up Team32 with updated hysteresis Functions (..._Vinch_K becoming ..._EB) +...

Up Team32 with updated hysteresis Functions (..._Vinch_K becoming ..._EB) + new options: Flag_UpdateSeparated + Flag_AdaptRelax_Version choices

Up Team32 with updated hysteresis Functions (..._Vinch_K becoming ..._EB) +...
3a8f2306 · Kevin Jacques · b48bd603 · 3a8f2306 · 3a8f2306 · 3a8f2306
Commit 3a8f2306 authored 6 years ago by Kevin Jacques
--- a/Team32/BH_team32materialdata.pro
+++ b/Team32/BH_team32materialdata.pro
+Function{
+  /*
+  // The version NON-any below as a form like this (harder to deal with):
+  // (form like this:		         ___    )
+  // (                            __|       )
+  // (                         __|          ) 
+  */
+/*List_H={0,42.9013,51.4518,73.9502,108.958,202.746,251.765,348.401,492.594,727.547,906.266,1232.82,5441.59,1e+09};
+List_H2={0,1840.52,2647.29,5468.63,11871.8,41105.8,63385.5,121383,242649,529325,821319,1.51984e+06,2.9611e+07,1e+18};
+List_Janhy={0,0.706605,0.750429,0.834449,0.927388,1.09949,1.16511,1.26115,1.34964,1.4243,1.45561,1.48941,1.56197,1.58323};
+List_Banhy={0,0.706659,0.750493,0.834542,0.927525,1.09974,1.16542,1.26158,1.35025,1.42521,1.45675,1.49096,1.56881,1258.22};
+List_Banhy2={0,0.499366,0.56324,0.696461,0.860303,1.20943,1.35821,1.59159,1.82319,2.03123,2.12212,2.22296,2.46117,1.58312e+06};
+List_Muanhy={0.0479842,0.0164717,0.0145863,0.0112852,0.00851268,0.00542425,0.00462902,0.00362107,0.00274111,0.00195893,0.00160742,0.00120939,0.0002883,1.25822e-06};
+List_Nuanhy={20.8402,60.7101,68.5573,88.6116,117.472,184.357,216.029,276.161,364.816,510.484,622.115,826.861,3468.61,794773};
+List_J={0,0.150626,0.249049,0.499732,0.750528,1.0153,1.095,1.2004,1.29992,1.39605,1.44988,1.4999,1.56197,1.58323};
+List_B={0,0.15068,0.249113,0.499825,0.750665,1.01556,1.09531,1.20084,1.30054,1.39697,1.45102,1.50145,1.56881,1258.22};
+List_B2={0,0.0227043,0.0620575,0.249825,0.563498,1.03136,1.19971,1.44202,1.6914,1.95152,2.10545,2.25434,2.46117,1.58312e+06};
+List_Mu={0.000572824,0.00351224,0.00484168,0.00675895,0.00688949,0.00500902,0.00435054,0.00344672,0.00264018,0.00192011,0.00160109,0.0012179,0.0002883,1.25822e-06};
+List_Nu={1745.74,284.719,206.54,147.952,145.149,199.64,229.857,290.131,378.762,520.804,624.574,821.086,3468.61,794773};
+*/
+// More points
+List_H={0,1,5.06122,9.12245,13.1837,17.2449,21.3061,25.3673,29.4286,33.4898,37.551,41.6122,42.9013,45.6735,49.7347,51.4518,53.7959,57.8571,61.9184,65.9796,70.0408,73.9502,74.102,78.1633,82.2245,86.2857,90.3469,94.4082,98.4694,102.531,106.592,108.958,110.653,114.714,118.776,122.837,126.898,130.959,135.02,139.082,143.143,147.204,151.265,155.327,159.388,163.449,167.51,171.571,175.633,179.694,183.755,187.816,191.878,195.939,200,202.746,251.765,348.401,492.594,727.547,906.266,1232.82,5441.59,1e+09};
+List_H2={0,1,25.616,83.2191,173.809,297.387,453.951,643.502,866.041,1121.57,1410.08,1731.58,1840.52,2086.07,2473.54,2647.29,2894,3347.45,3833.88,4353.31,4905.72,5468.63,5491.11,6109.5,6760.87,7445.22,8162.57,8912.9,9696.22,10512.5,11361.8,11871.8,12244.1,13159.4,14107.6,15088.9,16103.1,17150.3,18230.5,19343.7,20489.9,21669,22881.2,24126.3,25404.5,26715.6,28059.7,29436.8,30846.8,32289.9,33765.9,35275,36817,38392,40000,41105.8,63385.5,121383,242649,529325,821319,1.51984e+06,2.9611e+07,1e+18};
+List_Janhy={0,0.0309761,0.153936,0.266504,0.363752,0.444473,0.510097,0.563225,0.606585,0.642523,0.672877,0.699018,0.706605,0.721948,0.742395,0.750429,0.760892,0.777831,0.793504,0.808137,0.821899,0.834449,0.834925,0.84732,0.859169,0.87054,0.88149,0.892064,0.902301,0.912232,0.921884,0.927388,0.93128,0.94044,0.94938,0.958114,0.966655,0.975014,0.9832,0.991222,0.999086,1.0068,1.01437,1.0218,1.02909,1.03626,1.0433,1.05021,1.05701,1.06368,1.07025,1.0767,1.08304,1.08927,1.0954,1.09949,1.16511,1.26115,1.34964,1.4243,1.45561,1.48941,1.56197,1.58323};
+List_Banhy={0,0.0309774,0.153942,0.266515,0.363768,0.444494,0.510123,0.563257,0.606622,0.642565,0.672924,0.69907,0.706659,0.722005,0.742458,0.750493,0.76096,0.777903,0.793582,0.80822,0.821987,0.834542,0.835018,0.847418,0.859272,0.870649,0.881604,0.892183,0.902425,0.912361,0.922018,0.927525,0.931419,0.940584,0.949529,0.958268,0.966814,0.975178,0.98337,0.991397,0.999266,1.00699,1.01456,1.02199,1.02929,1.03646,1.04351,1.05043,1.05723,1.06391,1.07048,1.07693,1.08328,1.08952,1.09565,1.09974,1.16542,1.26158,1.35025,1.42521,1.45675,1.49096,1.56881,1258.22};
+List_Banhy2={0,0.000959599,0.0236982,0.0710303,0.132327,0.197575,0.260226,0.317259,0.367991,0.41289,0.452827,0.488699,0.499366,0.521292,0.551244,0.56324,0.57906,0.605133,0.629773,0.653219,0.675662,0.696461,0.697254,0.718117,0.738349,0.758029,0.777225,0.795991,0.81437,0.832402,0.850117,0.860303,0.867542,0.884698,0.901605,0.918278,0.93473,0.950973,0.967016,0.982867,0.998533,1.01402,1.02933,1.04447,1.05945,1.07426,1.08891,1.1034,1.11773,1.1319,1.14592,1.15978,1.17349,1.18705,1.20046,1.20943,1.35821,1.59159,1.82319,2.03123,2.12212,2.22296,2.46117,1.58312e+06};
+List_Muanhy={0.0309975,0.0309774,0.030416,0.0292153,0.0275923,0.0257754,0.0239426,0.022204,0.0206134,0.0191869,0.0179203,0.0167996,0.0164717,0.015808,0.0149284,0.0145863,0.0141453,0.0134452,0.0128166,0.0122495,0.0117358,0.0112852,0.0112685,0.0108416,0.0104503,0.0100903,0.00975798,0.00945027,0.00916452,0.00889842,0.00864999,0.00851268,0.00841747,0.00819936,0.00799431,0.00780115,0.00761883,0.00744643,0.00728312,0.00712816,0.0069809,0.00684074,0.00670716,0.00657965,0.0064578,0.00634121,0.00622951,0.00612239,0.00601953,0.00592067,0.00582556,0.00573396,0.00564568,0.00556051,0.00547827,0.00542425,0.00462902,0.00362107,0.00274111,0.00195893,0.00160742,0.00120939,0.0002883,1.25822e-06};
+List_Nuanhy={32.2607,32.2816,32.8775,34.2286,36.242,38.7967,41.7666,45.0369,48.5122,52.1189,55.8027,59.5251,60.7101,63.2592,66.9865,68.5573,70.6948,74.3758,78.0239,81.6357,85.2092,88.6116,88.7431,92.237,95.6909,99.1051,102.48,105.817,109.116,112.379,115.607,117.472,118.801,121.961,125.089,128.186,131.254,134.293,137.304,140.289,143.248,146.183,149.094,151.984,154.851,157.699,160.526,163.335,166.126,168.9,171.657,174.399,177.127,179.84,182.539,184.357,216.029,276.161,364.816,510.484,622.115,826.861,3468.61,794773};
+List_J={0,0.000102733,0.00257973,0.00821228,0.0168,0.0281424,0.0420391,0.0582897,0.0766939,0.0970511,0.119161,0.142823,0.150626,0.175788,0.227968,0.249049,0.27606,0.32408,0.372237,0.418886,0.46238,0.499732,0.501078,0.536706,0.571417,0.60484,0.636603,0.666336,0.693668,0.718228,0.739646,0.750528,0.757812,0.774727,0.790909,0.80639,0.821202,0.835378,0.848948,0.861945,0.874401,0.886347,0.897815,0.908838,0.919446,0.929673,0.93955,0.949108,0.95838,0.967397,0.976192,0.984796,0.993241,1.00156,1.00978,1.0153,1.095,1.2004,1.29992,1.39605,1.44988,1.4999,1.56197,1.58323};
+List_B={0,0.00010399,0.00258609,0.00822374,0.0168165,0.028164,0.0420659,0.0583216,0.0767308,0.0970931,0.119208,0.142875,0.15068,0.175845,0.228031,0.249113,0.276127,0.324153,0.372315,0.418969,0.462468,0.499825,0.501171,0.536805,0.571521,0.604948,0.636716,0.666454,0.693791,0.718357,0.73978,0.750665,0.757951,0.774871,0.791058,0.806545,0.821362,0.835542,0.849118,0.86212,0.87458,0.886532,0.898005,0.909033,0.919647,0.929879,0.93976,0.949324,0.958601,0.967623,0.976423,0.985032,0.993482,1.00181,1.01003,1.01556,1.09531,1.20084,1.30054,1.39697,1.45102,1.50145,1.56881,1258.22};
+List_B2={0,1.08138e-08,6.68786e-06,6.76299e-05,0.000282796,0.000793213,0.00176954,0.00340141,0.00588762,0.00942708,0.0142106,0.0204134,0.0227043,0.0309215,0.0519981,0.0620575,0.0762462,0.105075,0.138619,0.175535,0.213877,0.249825,0.251172,0.288159,0.326636,0.365962,0.405408,0.444161,0.481346,0.516036,0.547274,0.563498,0.57449,0.600425,0.625773,0.650514,0.674635,0.698131,0.721001,0.743251,0.764891,0.785938,0.806413,0.826341,0.84575,0.864674,0.883149,0.901215,0.918915,0.936295,0.953402,0.970288,0.987007,1.00361,1.02017,1.03136,1.19971,1.44202,1.6914,1.95152,2.10545,2.25434,2.46117,1.58312e+06};
+List_Mu={7.43501e-05,0.00010399,0.000510961,0.000901484,0.00127556,0.00163318,0.00197436,0.00229908,0.00260736,0.00289919,0.00317456,0.00343349,0.00351224,0.00385005,0.00458495,0.00484168,0.00513286,0.00560264,0.006013,0.00634998,0.00660284,0.00675895,0.00676325,0.00686774,0.00695074,0.00701099,0.00704746,0.00705929,0.00704576,0.00700627,0.0069403,0.00688949,0.0068498,0.00675479,0.00666011,0.00656599,0.00647262,0.00638017,0.00628881,0.00619866,0.00610984,0.00602247,0.00593662,0.0058524,0.00576987,0.00568911,0.00561017,0.00553311,0.00545799,0.00538484,0.00531372,0.00524466,0.00517769,0.00511285,0.00505017,0.00500902,0.00435054,0.00344672,0.00264018,0.00192011,0.00160109,0.0012179,0.0002883,1.25822e-06};
+List_Nu={13449.9,9616.35,1957.1,1109.28,783.971,612.302,506.494,434.956,383.53,344.924,315.004,291.249,284.719,259.737,218.105,206.54,194.823,178.487,166.306,157.481,151.45,147.952,147.858,145.608,143.87,142.633,141.895,141.657,141.929,142.729,144.086,145.149,145.99,148.043,150.148,152.3,154.497,156.736,159.013,161.325,163.67,166.045,168.446,170.87,173.314,175.775,178.248,180.73,183.218,185.706,188.192,190.67,193.136,195.586,198.013,199.64,229.857,290.131,378.762,520.804,624.574,821.086,3468.61,794773};
+any_nu_b2  = ListAlt[List_B2, List_Nu] ;
+//any_nu_b2  = ListAlt[List_Banhy2, List_Nuanhy] ; // take the any version (see header)
+// Akima Interpolation problematic with Nu ...
+nu_any[] = InterpolationAkima[$1]{List[any_nu_b2]} ;
+dnudb2_any[] = dInterpolationAkima[$1]{List[any_nu_b2]} ;
+// Try Linear Interpolation instead with Nu ...
+//nu_any[] = InterpolationLinear[$1]{List[any_nu_b2]} ;
+//dnudb2_any[] = dInterpolationLinear[$1]{List[any_nu_b2]} ;
+h_any[] = nu_any[(SquNorm[$1])] * $1 ;
+dhdb_any[] = TensorDiag[1,1,1] * nu_any[SquNorm[$1]] + 2*dnudb2_any[SquNorm[$1]] * SquDyadicProduct[$1]  ;  
+dhdb_any_NL[] = 2*dnudb2_any[SquNorm[$1]] * SquDyadicProduct[$1]  ;  
+//nu[DomainNL]=nu_any[(SquNorm[$1])]; // already in magstadyn_a.pro
+//dnudbxb[] = 2*dnudb2_any[(SquNorm[$1])] * SquDyadicProduct[$1]; // defined here abobe as dhdb_any_NL
+}
--- a/Team32/MyCircuit.pro
+++ b/Team32/MyCircuit.pro
@@ -69,16 +69,16 @@ Constraint {
      If (Flag_TestCase==3)
      Case Circuit1 {
-        { Region Ein1 ; Branch {1,2} ; }
+        { Region Ein1 ; Branch {2,1} ; } // 2 1 since 25/8/18 !! (E1Inv)
        { Region R1 ;   Branch {2,3} ; }
        { Region IndB~{1} ; Branch {3,4} ; }
-        { Region IndA~{1} ; Branch {4,1} ; }
+        { Region IndA~{1} ; Branch {1,4} ; } // 1 4 since 25/8/18 !!
      }
      Case Circuit2 {
-        { Region Ein2 ; Branch {1,2} ; }
+        { Region Ein2 ; Branch {1,2} ; } //(E2dir)
        { Region R2   ; Branch {2,3} ; }
        { Region IndB~{2} ; Branch {3,4} ; }
-        { Region IndA~{2} ; Branch {4,1} ; }
+        { Region IndA~{2} ; Branch {1,4} ; } // 4 1 since 25/8/18 !!
      }
      EndIf
      If (Flag_TestCase==4)
@@ -108,7 +108,7 @@ Constraint {
      If (Flag_TestCase==3)
      Case Circuit1 {
-        { Region Ein1 ; Branch {1,2} ; }
+        { Region Ein1 ; Branch {2,1} ; } // 2 1 since 12/12/18 !! (E1Inv)
        { Region R1 ;   Branch {2,3} ; }
        { Region Ind~{1} ; Branch {3,1} ; }
      }

--- a/Team32/MyItLoopPro.pro
+++ b/Team32/MyItLoopPro.pro
@@ -42,8 +42,9 @@ Macro MyItLoopPro
  Test[Flag_AdaptRelax==1]{
    CopySolution[A,'x_Prev']; //SET : save solution from the previous converged time step at iter 0
-    Evaluate[ $res_Prev=$res,
+    Evaluate[ $res_Opt=$res,
-              $relax_Opt=$relax];   
+              $relax_Opt=$relax,
+              $inc_Opt=0];   
  }
  {
    Generate[A]; 
@@ -60,7 +61,12 @@ Macro MyItLoopPro
    // if adapted relaxation is activated 
    Test[Flag_AdaptRelax==1]
    {
+      If(Flag_AdaptRelax_Version ==1 )
        Call MySolveJac_AdaptRelax;
+      EndIf
+      If(Flag_AdaptRelax_Version ==2 )
+        Call MySolveJac_AdaptRelax2;
+      EndIf
    }
    // ... otherwise, if adapted relaxation is not activated 
@@ -69,12 +75,13 @@ Macro MyItLoopPro
      Evaluate[ $syscount = $syscount + 1, 
                $relaxcount=1, 
                $relaxcounttot=$relaxcounttot+1 ];
-      Generate[A]; GetResidual[A, $res]; 
+      Generate[A]; 
+      GetNormIncrement[A, $inc]; //Check the new increment
      // TODO: GetIncrement h_only   
      //       " PostOperation { { h_only , Reltol_Mag, Abstol_Mag,  MeanL2Norm }} "
+      GetResidual[A, $res]; 
+      Test[$iter==1]{Evaluate[$res0=$res];} // compute first residual at step $TimeStep
    }
-    GetNormIncrement[A, $inc]; //Check the new increment
    GetNormSolution[A, $sol]; //Check the new solution
    Call DisplayRunTimeVar;
@@ -90,51 +97,55 @@ Macro MySolveJac_AdaptRelax2
  //NB:solution from the previous converged time step at iter 0 is saved in x_Prev at this point
  Evaluate[ $redo=1, 
            $relax=relax_max, 
-            $res_Prev=$res, 
+            $res_Opt=$res, 
            $relaxcount=0];
  /*  get increment dx by solving the system */
  Generate[A];
  Solve[A];
  Evaluate[ $syscount = $syscount + 1];
-  CopySolution[A,'x_New']; //SET : save increment dx (found without relaxation $relax=1)
+  CopySolution[A,'x_New']; 
-  //TODO: "DX = x_New - x_Prev" 
+  //SET : save increment dx (found without relaxation $relax=relax_max)
-  // or TODO: CopyIncrement[A,'DX'];
+  AddVector[A, 1/relax_max, 'x_New', -1/relax_max, 'x_Prev', 'Delta_x']; 
  While[ $relax>relax_min-0.5*relax_step &&
         $redo>0 ]
  // while all relaxation factor have not been tested yet 
  {
    /* new trial solution = x + Frelax * dx */
-    // TODO:"x_New = x_Prev + $relax * DX"
+    AddVector[A, 1, 'x_Prev', $relax, 'Delta_x', 'x_New'];
    CopySolution['x_New', A]; // RESET: reset new trial solution in the system
    /* calculate residual with trial solution res= b-A(x_New)*x_New*/
    Generate[A];  // regenerate A with new trial solution (x_New) 
-    GetResidual[A, $res]; //Check the new residual  
    // TODO: GetIncrement h_only   
    //       " PostOperation { { h_only , Reltol_Mag, Abstol_Mag,  MeanL2Norm }} "
+    AddVector[A, 1, 'x_New', -1, 'x_Prev', 'weighted_Delta_x'];     
+    CopyIncrement['weighted_Delta_x',A];
+    GetNormIncrement[A, $inc]; //Check the new increment
+    GetResidual[A, $res]; //Check the new residual
    Evaluate[ $relaxcount=$relaxcount+1, 
              $relaxcounttot=$relaxcounttot+1 ];
    // compute first residual at step $TimeStep
    Test[$iter==1 && $relax==relax_max]{
      Evaluate[ $res0=$res, 
-                $res_Prev=$res, 
+                $res_Opt=$res, 
-                $relax_Opt=$relax ];
+                $relax_Opt=$relax,
+                $inc_Opt=$inc ];
      CopySolution[A,'x_Opt']; //SET : save actual solution as the optimal one  
    }        
    /* // If one wants to be more verbose
-    Print[{$TimeStep, $iter, $res, $res/$res0, $inc, $relax, ($res < $res_Prev ), $relaxcount}, Format
+    Print[{$TimeStep, $iter, $res, $res/$res0, $inc, $relax, ($res < $res_Opt ), $relaxcount}, Format
-      "ts=%05g, it=%03g, abs=%1.5e, rel=%1.5e, dx=%1.5e, relax = %1.1g, (better? %g), relaxcount=%02g"];
+      "ts=%05g, it=%03g, abs=%1.5e, rel=%1.5e, dx=%1.5e, relax = %1.5g, (better? %g), relaxcount=%02g (V2)"];
    //*/
    // if the residual has decreased
-    Test[ ($res < $res_Prev ) || $relax>relax_max-0.5*relax_step ]{    
+    Test[ ($res < $res_Opt ) || $relax==relax_max ]{    
      CopySolution[A,'x_Opt']; //SET : save actual solution as the optimal one  'x_New is best'
-      Evaluate[$res_Prev=$res,
+      Evaluate[$res_Opt=$res,
-               $relax_Opt=$relax];                     
+               $relax_Opt=$relax,
+               $inc_Opt=$inc];                     
    }
    { // otherwise, if the residual has not decreased... 
      // and if you don't want to check All Factors: Break the try loop prematurely
@@ -145,8 +156,9 @@ Macro MySolveJac_AdaptRelax2
  }
  CopySolution['x_Opt',A]; //RESET : reload optimal solution
  CopySolution[A,'x_Prev']; //SET : save solution from the previous converged iteration              
-  Evaluate[$relax=$relax_Opt,
+  Evaluate[$res=$res_Opt,
-           $res=$res_Prev];
+           $relax=$relax_Opt,
+           $inc=$inc_Opt];
 Return
 //...............................................................
@@ -158,7 +170,7 @@ Macro MySolveJac_AdaptRelax
  Evaluate[ $redo=1, 
            $relax=relax_max, 
-            $res_Prev=$res, 
+            $res_Opt=$res, 
            $relaxcount=0];
  While[ $relax>relax_min-0.5*relax_step &&
@@ -172,26 +184,29 @@ Macro MySolveJac_AdaptRelax
              $relaxcount=$relaxcount+1, 
              $relaxcounttot=$relaxcounttot+1 ];
    Generate[A]; 
+    GetNormIncrement[A, $inc]; //Check the new increment
    GetResidual[A, $res]; //Check the new residual 
    // compute first residual at step $TimeStep
    Test[$iter==1 && $relax==relax_max]{
      Evaluate[ $res0=$res, 
-                $res_Prev=$res, 
+                $res_Opt=$res, 
-                $relax_Opt=$relax ];
+                $relax_Opt=$relax,
+                $inc_Opt=$inc ];
      CopySolution[A,'x_Opt']; //SET : save actual solution as the optimal one}  
    }             
    /* // If one wants to be more verbose
-    Print[{$TimeStep, $iter, $res, $res/$res0, $inc, $relax, ($res < $res_Prev ), $relaxcount}, Format
+    Print[{$TimeStep, $iter, $res, $res/$res0, $inc, $relax, ($res < $res_Opt ), $relaxcount}, Format
-      "ts=%05g, it=%03g, abs=%1.5e, rel=%1.5e, dx=%1.5e, relax = %1.1g, (better? %g), relaxcount=%02g"];
+      "ts=%05g, it=%03g, abs=%1.5e, rel=%1.5e, dx=%1.5e, relax = %1.5g, (better? %g), relaxcount=%02g (V1)"];
    //*/
    // if the residual has decreased
-    Test[ ($res < $res_Prev ) || $relax>relax_max-0.5*relax_step ]{    
+    Test[ ($res < $res_Opt ) || $relax==relax_max ]{    
      CopySolution[A,'x_Opt']; //SET : save actual solution as the optimal one
-      Evaluate[$res_Prev=$res,
+      Evaluate[$res_Opt=$res,
-               $relax_Opt=$relax];                     
+               $relax_Opt=$relax,
+               $inc_Opt=$inc];                     
    }
    { // otherwise, if the residual has not decreased... 
      // and if you don't want to check All Factors: Break the try loop prematurely
@@ -202,8 +217,9 @@ Macro MySolveJac_AdaptRelax
  }
  CopySolution['x_Opt',A]; //RESET : reload optimal solution
  CopySolution[A,'x_Prev']; //SET : save solution from the previous converged iteration              
-  Evaluate[$relax=$relax_Opt,
+  Evaluate[$res=$res_Opt,
-           $res=$res_Prev];
+           $relax=$relax_Opt,
+           $inc=$inc_Opt];
 Return
 //...............................................................
@@ -215,12 +231,12 @@ Macro DisplayRunTimeVar
    If(Flag_NLRes!=NLRES_ITERATIVELOOPPRO) 
      Print[{$TimeStep, $iter, $res, $resL, $resN, $relax, $relaxcount}, Format
-      "*** ts=%05g, it=%03g, res_used=%1.5e, (res = %1.5e, resN=%1.5e), relaxopt = %1.3g, relaxcount=%02g ***"];
+      "*** ts=%05g, it=%03g, res_used=%1.5e, (res = %1.5e, resN=%1.5e), relaxopt = %1.5g, relaxcount=%02g ***"];
    EndIf
    If(Flag_NLRes==NLRES_ITERATIVELOOPPRO) 
      Print[{$TimeStep, $iter, $res, $res/$res0, $inc, $relax, $relaxcount}, Format
-        "*** ts=%05g, it=%03g, abs=%1.5e, rel=%1.5e, dx=%1.5e, relaxopt = %1.3g, relaxcount=%02g ***"];
+        "*** ts=%05g, it=%03g, abs=%1.5e, rel=%1.5e, dx=%1.5e, relaxopt = %1.5g, relaxcount=%02g ***"];
    EndIf
    PostOperation[PostOpItLoop_a];

--- a/Team32/magstadyn_a.pro
+++ b/Team32/magstadyn_a.pro
--- a/Team32/magstadyn_avs_3d.pro
+++ b/Team32/magstadyn_avs_3d.pro
--- a/Team32/param_EnergHyst.dat
+++ b/Team32/param_EnergHyst.dat
 ///*//default for t32 #CHECK here
 dim00 = (Flag_3Dmodel==1)?3:2;
 N00   = 3; 
-FLAG_TANORLANG00 = 1; 
+FLAG_ANHYLAW00 = 2; 
-FLAG_VARORDIFF00 = 2; 
+FLAG_APPROACH00 = 2; 
 FLAG_INVMETHOD00 = 1; 
-FLAG_ROOTFINDING1D00 = 0; 
+FLAG_JACEVAL00 = 1; 
 FLAG_MINMETHOD00 = 1; 
 FLAG_WARNING00 = 0 ;     
 TOLERANCE_JS00 = 1.e-3;  
@@ -23,12 +23,12 @@ FLAG_HOMO00    = 1;
 /*
 dim00 = (Flag_3Dmodel==1)?3:2;
 N00   = 3; 
-FLAG_TANORLANG00 = 1; 
+FLAG_ANHYLAW00 = 1; 
-FLAG_VARORDIFF00 = 2; 
+FLAG_APPROACH00 = 2; 
 FLAG_INVMETHOD00 = 1; 
-FLAG_ROOTFINDING1D00 = 0; 
+FLAG_JACEVAL00 = 1; 
 FLAG_MINMETHOD00 = 1; 
-FLAG_WARNING00 = 3 ;     
+FLAG_WARNING00 = 0 ;     
 TOLERANCE_JS00 = 1.e-3;  
 TOLERANCE_000  = 1.e-8; 
 TOLERANCE_NR00 = 1.e-8;  

--- a/Team32/param_EnergHyst.pro
+++ b/Team32/param_EnergHyst.pro
@@ -5,39 +5,27 @@ Function {
  dim = dim00;
  N   = N00; //nbr de cercle
-  FLAG_TANORLANG = FLAG_TANORLANG00; 
+  FLAG_ANHYLAW = FLAG_ANHYLAW00; 
 /* 
-FLAG_TANORLANG = 1 --> hyperbolic tangent
+FLAG_ANHYLAW = 1 --> hyperbolic tangent
-FLAG_TANORLANG = 2 --> double langevin function
+FLAG_ANHYLAW = 2 --> double langevin function
 */
-  FLAG_VARORDIFF = FLAG_VARORDIFF00; 
+  FLAG_APPROACH = FLAG_APPROACH00; 
 /*
-FLAG_VARORDIFF = 1 --> variational approach (Jk)
+FLAG_APPROACH = 1 --> variational approach (Jk)
-FLAG_VARORDIFF = 2 --> simple differential approach (hrk)
+FLAG_APPROACH = 2 --> vector play model approach (hrk)
-FLAG_VARORDIFF = 3 --> full differential approach (hrk)
+FLAG_APPROACH = 3 --> full differential approach (hrk)
 */
  FLAG_INVMETHOD = FLAG_INVMETHOD00; 
 /*
 FLAG_INVMETHOD = 1 --> NR_ana (homemade)
 FLAG_INVMETHOD = 2 --> NR_num (homemade)
 FLAG_INVMETHOD = 3 --> bfgs (homemade)
-FLAG_INVMETHOD = 4 --> hybridsj (gsl)
-FLAG_INVMETHOD = 5 --> hybridj (gsl)
-FLAG_INVMETHOD = 6 --> newton (gsl)
-FLAG_INVMETHOD = 7 --> gnewton (gsl)
 */
-  FLAG_ROOTFINDING1D = FLAG_ROOTFINDING1D00; 
+  FLAG_JACEVAL = FLAG_JACEVAL00; 
 /*
-FLAG_ROOTFINDING1D  = 0 --> Not done
+FLAG_JACEVAL = 1 --> analytical Jacobian
-FLAG_ROOTFINDING1D  = 1 --> (1D minimization) golden section (gsl)
+FLAG_JACEVAL = 2 --> numerical Jacobian
-FLAG_ROOTFINDING1D  = 2 --> (1D minimization) brent (gsl)
-FLAG_ROOTFINDING1D  = 3 --> (1D minimization) quad golden (gsl)
-FLAG_ROOTFINDING1D  = 4 --> (1D root bracketing) bisection (gsl)
-FLAG_ROOTFINDING1D  = 5 --> (1D root bracketing) brent (gsl)
-FLAG_ROOTFINDING1D  = 6 --> (1D root bracketing) falsepos (gsl)
-FLAG_ROOTFINDING1D  = 7 --> (1D root finding with derivatives) newton (gsl)
-FLAG_ROOTFINDING1D  = 8 --> (1D root finding with derivatives) secant (gsl)
-FLAG_ROOTFINDING1D  = 9 --> (1D root finding with derivatives) steffenson (gsl)
 */
  FLAG_MINMETHOD = FLAG_MINMETHOD00; 
 /*
@@ -47,18 +35,21 @@ FLAG_MINMETHOD = 3 --> conjugate pr (gsl)
 FLAG_MINMETHOD = 4 --> bfgs2 (gsl)
 FLAG_MINMETHOD = 5 --> bfgs (gsl)
 FLAG_MINMETHOD = 6 --> steepest descent (gsl)
+FLAG_MINMETHOD = 11   --> steepest descent+ (homemade)\n"
+FLAG_MINMETHOD = 22   --> conjugate Fletcher-Reeves (homemade)\n"
+FLAG_MINMETHOD = 33   --> conjugate Polak-Ribiere (homemade)\n"
+FLAG_MINMETHOD = 333  --> conjugate Polak-Ribiere+ (homemade)\n"
+FLAG_MINMETHOD = 1999 --> conjugate Dai Yuan 1999 (p.85) (homemade)\n"
+FLAG_MINMETHOD = 2005 --> conjugate Hager Zhang 2005 (p.161) (homemade)\n"
+FLAG_MINMETHOD = 77   --> newton (homemade)\n", ::FLAG_MINMETHOD);
 */
  FLAG_WARNING = FLAG_WARNING00 ;     // SHOW WARNING 
 /*
 #define FLAG_WARNING_INFO_INV         1
-#define FLAG_WARNING_INFO_ROOTFINDING 2
+#define FLAG_WARNING_INFO_APPROACH    2
-#define FLAG_WARNING_INFO_MIN         3
+#define FLAG_WARNING_STOP_INV         10
-#define FLAG_WARNING_DISP_INV         11
+#define FLAG_WARNING_DISPABOVEITER    1
-#define FLAG_WARNING_DISP_ROOTFINDING 12
-#define FLAG_WARNING_STOP_INV         101
-#define FLAG_WARNING_STOP_ROOTFINDING 102
-#define FLAG_WARNING_ITER             100
 */
  TOLERANCE_JS = TOLERANCE_JS00;  // SENSITIVE_PARAM (1.e-3) // 1.e-4 
  TOLERANCE_0  = TOLERANCE_000;  // SENSITIVE_PARAM (1.e-7)
@@ -73,14 +64,14 @@ FLAG_MINMETHOD = 6 --> steepest descent (gsl)
                        //                  1.e0 for VinchT & transfo)
  FLAG_HOMO    = FLAG_HOMO00;
-  If(FLAG_TANORLANG==1)
+  If(FLAG_ANHYLAW==1)
  // Material M25050A
    Ja   = 1.22;
    ha   = 65;
    Jb   = 0.;
    hb   = 0.;   
  EndIf
-  If(FLAG_TANORLANG==2)
+  If(FLAG_ANHYLAW==2)
 /* 
    // Material M23535A (default for square) #CHECK here
    Ja   = 0.595;
@@ -113,86 +104,73 @@ Jb   = 0.572282990517
 hb   = 399.474553851
 */
  EndIf
-  If(N==3)
-///*
+  If(N==1)
+    kappa_1 = 0.;
+    w_1   = 1.;
+  ElseIf (N==3)
+/*
  // Material M25050A (default for square) #CHECK here
    Js_1  = 0.11;
    Js_2  = 0.8;
    Js_3  = 0.31;
-    chi_1 = 0;
+    kappa_1 = 0;
-    chi_2 = 16;
+    kappa_2 = 16;
-    chi_3 = 47;
+    kappa_3 = 47;
    w_1   = 0.090163934426230;
    w_2   = 0.655737704918033;
    w_3   = 0.254098360655738;
 /*/
 /*    //team32_RD_z_anhyshift (abstract ISEM 2017)
-    chi_1 = 0*11.6413;
+    kappa_1 = 0*11.6413;
-    chi_2 = 50.1413;
+    kappa_2 = 50.1413;
-    chi_3 = 107.518;
+    kappa_3 = 107.518;
    w_1   = 0.173524;
    w_2   = 0.591535;
    w_3   = 0.234941;
-    Js_1  = w_1*(Ja+Jb);
-    Js_2  = w_2*(Ja+Jb);
-    Js_3  = w_3*(Ja+Jb);
 //*/ 
 /*    //FH(xini=400) // THIS IS PROBLEMATIC (CONVERGENCE PROBLEM ok now) (new try)
-    chi_1 = 0.;
+    kappa_1 = 0.;
-    chi_2 = 53.779614004;
+    kappa_2 = 53.779614004;
-    chi_3 = 233.384447699;
+    kappa_3 = 233.384447699;
    w_1   = 0.1048; //+0.1 to converge
    w_2   = 0.817943263499; // -0.1 to converge
    w_3   = 0.0772371911006;
-    Js_1  = w_1*(Ja+Jb);
-    Js_2  = w_2*(Ja+Jb);
-    Js_3  = w_3*(Ja+Jb);
 //*/
 /*    //FH(xini=400) // NEW 1/9/2017 (chamonix) (default for t32) #CHECK here
-    chi_1 = 0.;
+    kappa_1 = 0.;
-    chi_2 = 53.789872;
+    kappa_2 = 53.789872;
-    chi_3 = 233.910892;
+    kappa_3 = 233.910892;
    w_1   = 0.0280702473932;
    w_2   = 0.896565361832;
    w_3   = 0.0753643907746;
-    Js_1  = w_1*(Ja+Jb);
-    Js_2  = w_2*(Ja+Jb);
-    Js_3  = w_3*(Ja+Jb);
 //*/
 /*//FH(xini=400) TD // NEW 1/9/2017 (chamonix)
 0 0.0565084014449 0.0
 1 0.839372763961 65.418004
 2 0.104118834594 246.564466
 */
-  EndIf
+///*   // NEW 8/8/2018  for team32
-  If(N==1)
+    kappa_1 = 0.;
-    Js_1  = 1;
+    kappa_2 = 54.705531;
-    Js_2  = 0.;
+    kappa_3 = 199.248421;
-    Js_3  = 0.;    
+    w_1   = 0.028745;
-    chi_1 = 0.;
+    w_2   = 0.921797;
-    chi_2 = 0.;
+    w_3   = 0.049458;
-    chi_3 = 0.;
+//*/
-    w_1   = 1.;
+ElseIf (N==5)   //FH(xini=400) // NEW 1/9/2017
-    w_2   = 0.;
+    kappa_1 = 0.;
-    w_3   = 0.;
+    kappa_2 = 14.81413;
-  EndIf
+    kappa_3 = 52.317547;
-  If(N==5)   //FH(xini=400) // NEW 1/9/2017
+    kappa_4 = 102.345013;
-    chi_1 = 0.;
+    kappa_5 = 247.012801;
-    chi_2 = 14.81413;
-    chi_3 = 52.317547;
-    chi_4 = 102.345013;
-    chi_5 = 247.012801;
    w_1   = 0.028258057159;
    w_2   = 0.087100091228;
    w_3   = 0.702372713403;
    w_4   = 0.134380357677;
    w_5   = 0.0478887805338;
-    Js_1  = w_1*(Ja+Jb);
-    Js_2  = w_2*(Ja+Jb);
-    Js_3  = w_3*(Ja+Jb);   
-    Js_4  = w_4*(Ja+Jb);
-    Js_5  = w_5*(Ja+Jb);
  /*//FH(xini=400) TD // NEW 1/9/2017
  0 0.0574009326857 0.0
@@ -202,17 +180,15 @@ hb   = 399.474553851
  4 0.0646106266745 259.419756*/
 EndIf
-  If (N==3)
+  If (N==1)
  param_EnergHyst={ dim, N,  
                    Ja, ha, Jb, hb, 
-                    w_1, chi_1,
+                    w_1, kappa_1,
-                    w_2, chi_2,
-                    w_3, chi_3,
                    FLAG_INVMETHOD,
-                    FLAG_ROOTFINDING1D, 
+                    FLAG_JACEVAL, 
-                    FLAG_VARORDIFF, 
+                    FLAG_APPROACH, 
                    FLAG_MINMETHOD,
-                    FLAG_TANORLANG,
+                    FLAG_ANHYLAW,
                    FLAG_WARNING,
                    TOLERANCE_JS,
                    TOLERANCE_0,
@@ -225,21 +201,42 @@ hb   = 399.474553851
                    DELTA_0,
                    FLAG_HOMO
                   };
-  EndIf
+  ElseIf (N==3)
+  param_EnergHyst={ dim, N,  
-  If (N==5)
+                    Ja, ha, Jb, hb, 
+                    w_1, kappa_1,
+                    w_2, kappa_2,
+                    w_3, kappa_3,
+                    FLAG_INVMETHOD,
+                    FLAG_JACEVAL, 
+                    FLAG_APPROACH, 
+                    FLAG_MINMETHOD,
+                    FLAG_ANHYLAW,
+                    FLAG_WARNING,
+                    TOLERANCE_JS,
+                    TOLERANCE_0,
+                    TOLERANCE_NR,
+                    MAX_ITER_NR,
+                    TOLERANCE_OM,
+                    MAX_ITER_OM,
+                    FLAG_ANA,
+                    TOLERANCE_NJ,
+                    DELTA_0,
+                    FLAG_HOMO
+                   };
+  ElseIf (N==5)
  param_EnergHyst={ dim, N,  
                    Ja, ha, Jb, hb, 
-                    w_1, chi_1,
+                    w_1, kappa_1,
-                    w_2, chi_2,
+                    w_2, kappa_2,
-                    w_3, chi_3,
+                    w_3, kappa_3,
-                    w_4, chi_4,
+                    w_4, kappa_4,
-                    w_5, chi_5,
+                    w_5, kappa_5,
                    FLAG_INVMETHOD,
-                    FLAG_ROOTFINDING1D, 
+                    FLAG_JACEVAL, 
-                    FLAG_VARORDIFF, 
+                    FLAG_APPROACH, 
                    FLAG_MINMETHOD,
-                    FLAG_TANORLANG,
+                    FLAG_ANHYLAW,
                    FLAG_WARNING,
                    TOLERANCE_JS,
                    TOLERANCE_0,
@@ -252,16 +249,18 @@ hb   = 399.474553851
                    DELTA_0,
                    FLAG_HOMO
                   };
+  Else
+    Printf["N=%g is not valid",N];
+    Error["The parameters (w_k, kappa_k) for the asked number of cells (N) is not written yet in 'param_EnergHyst.pro'"];
  EndIf
 //*******************************************************************
 // DISPLAY PARAMETERS SETTING
 //*******************************************************************
 Label_FLAG_INVMETHOD="...";
-Label_FLAG_ROOTFINDING1D="...";
+Label_FLAG_JACEVAL="...";
-Label_FLAG_TANORLANG="...";
+Label_FLAG_ANHYLAW="...";
-Label_FLAG_VARORDIFF="...";
+Label_FLAG_APPROACH="...";
 Label_FLAG_MINMETHOD="...";
 Label_FLAG_ANA="...";
  If (FLAG_INVMETHOD == 1)
@@ -273,102 +272,66 @@ Label_FLAG_ANA="...";
  If (FLAG_INVMETHOD == 3)
  Label_FLAG_INVMETHOD="FLAG_INVMETHOD=3    --> bfgs (homemade)";
  EndIf
-  If (FLAG_INVMETHOD == 4)
+  If(FLAG_JACEVAL == 1)
-  Label_FLAG_INVMETHOD="FLAG_INVMETHOD=4     --> hybridsj (gsl)";
+  Label_FLAG_JACEVAL="FLAG_JACEVAL=1      --> analytical Jacobian";
-  EndIf
-  If (FLAG_INVMETHOD == 5)
-  Label_FLAG_INVMETHOD="FLAG_INVMETHOD=5     --> hybridj (gsl)";
-  EndIf
-  If (FLAG_INVMETHOD == 6)
-  Label_FLAG_INVMETHOD="FLAG_INVMETHOD=6     --> newton (gsl)";
-  EndIf
-  If (FLAG_INVMETHOD == 7)
-  Label_FLAG_INVMETHOD="FLAG_INVMETHOD=7     --> gnewton (gsl)";
-  EndIf
-  If(FLAG_ROOTFINDING1D == 0 && FLAG_INVMETHOD != 0)
-  Label_FLAG_ROOTFINDING1D="FLAG_ROOTFINDING1D=0 --> not done";
-  EndIf
-  If(FLAG_ROOTFINDING1D == 1 && FLAG_INVMETHOD != 0)
-  Label_FLAG_ROOTFINDING1D="FLAG_ROOTFINDING1D=1 --> (1D minimization) golden section (gsl)";
-  EndIf
-  If(FLAG_ROOTFINDING1D == 2 && FLAG_INVMETHOD != 0)
-  Label_FLAG_ROOTFINDING1D="FLAG_ROOTFINDING1D=2 --> (1D minimization) brent (gsl)";
-  EndIf
-  If(FLAG_ROOTFINDING1D == 3 && FLAG_INVMETHOD != 0)
-  Label_FLAG_ROOTFINDING1D="FLAG_ROOTFINDING1D=3 --> (1D minimization) quad golden (gsl)";
-  EndIf
-  If(FLAG_ROOTFINDING1D == 4 && FLAG_INVMETHOD != 0)
-  Label_FLAG_ROOTFINDING1D="FLAG_ROOTFINDING1D=4 --> (1D root bracketing) bisection (gsl)";
-  EndIf
-  If(FLAG_ROOTFINDING1D == 5 && FLAG_INVMETHOD != 0)
-  Label_FLAG_ROOTFINDING1D="FLAG_ROOTFINDING1D=5 --> (1D root bracketing) brent (gsl)";
-  EndIf
-  If(FLAG_ROOTFINDING1D == 6 && FLAG_INVMETHOD != 0)
-  Label_FLAG_ROOTFINDING1D="FLAG_ROOTFINDING1D=6 --> (1D root bracketing) falsepos (gsl)";
-  EndIf
-  If(FLAG_ROOTFINDING1D == 7 && FLAG_INVMETHOD != 0)
-  Label_FLAG_ROOTFINDING1D="FLAG_ROOTFINDING1D=7 --> (1D root finding with derivatives) newton (gsl)";
-  EndIf
-  If(FLAG_ROOTFINDING1D == 8 && FLAG_INVMETHOD != 0)
-  Label_FLAG_ROOTFINDING1D="FLAG_ROOTFINDING1D=8 --> (1D root finding with derivatives) secant (gsl)";
  EndIf
-  If(FLAG_ROOTFINDING1D == 9 && FLAG_INVMETHOD != 0)
+  If(FLAG_JACEVAL == 2)
-  Label_FLAG_ROOTFINDING1D="FLAG_ROOTFINDING1D=9 --> (1D root finding with derivatives) steffenson (gsl)";
+  Label_FLAG_JACEVAL="FLAG_JACEVAL=2      --> numerical Jacobian";
  EndIf
-  If (FLAG_TANORLANG == 1)
+  If (FLAG_ANHYLAW == 1)
-  Label_FLAG_TANORLANG="FLAG_TANORLANG=1     --> tangent hyperbolic";
+  Label_FLAG_ANHYLAW="FLAG_ANHYLAW=1     --> tangent hyperbolic";
  EndIf
-  If (FLAG_TANORLANG == 2)
+  If (FLAG_ANHYLAW == 2)
-  Label_FLAG_TANORLANG="FLAG_TANORLANG=2     --> double langevin function";
+  Label_FLAG_ANHYLAW="FLAG_ANHYLAW=2     --> double langevin function";
  EndIf
-  If (FLAG_VARORDIFF == 1)
+  If (FLAG_APPROACH == 1)
-  Label_FLAG_VARORDIFF="FLAG_VARORDIFF=1     --> variational approach (Jk)";
+  Label_FLAG_APPROACH="FLAG_APPROACH=1     --> variational approach (Jk)";
  EndIf
-  If (FLAG_VARORDIFF == 2)
+  If (FLAG_APPROACH == 2)
-  Label_FLAG_VARORDIFF="FLAG_VARORDIFF=2     --> simple differential approach (hrk)";
+  Label_FLAG_APPROACH="FLAG_APPROACH=2     --> vector play model approach (hrk)";
  Label_FLAG_MINMETHOD ="...";
  EndIf
-  If (FLAG_VARORDIFF == 3)
+  If (FLAG_APPROACH == 3)
-  Label_FLAG_VARORDIFF="FLAG_VARORDIFF=2     --> full differential approach (hrk)";
+  Label_FLAG_APPROACH="FLAG_APPROACH=3     --> full differential approach (hrk)";
  Label_FLAG_MINMETHOD ="...";
  EndIf
-  If (FLAG_MINMETHOD == 1 && FLAG_VARORDIFF==1)
+  If (FLAG_MINMETHOD == 1 && FLAG_APPROACH==1)
  Label_FLAG_MINMETHOD ="FLAG_MINMETHOD=1     --> steepest descent (homemade)";
  EndIf
-  If (FLAG_MINMETHOD == 11 && FLAG_VARORDIFF==1)
+  If (FLAG_MINMETHOD == 11 && FLAG_APPROACH==1)
  Label_FLAG_MINMETHOD ="FLAG_MINMETHOD=11    --> NEW steepest descent naive (homemade)";
  EndIf
-  If (FLAG_MINMETHOD == 22 && FLAG_VARORDIFF==1)
+  If (FLAG_MINMETHOD == 22 && FLAG_APPROACH==1)
  Label_FLAG_MINMETHOD ="FLAG_MINMETHOD=11    --> NEW conjugate Fletcher-Reeves (homemade)";
  EndIf
-  If (FLAG_MINMETHOD == 33 && FLAG_VARORDIFF==1)
+  If (FLAG_MINMETHOD == 33 && FLAG_APPROACH==1)
  Label_FLAG_MINMETHOD ="FLAG_MINMETHOD=11    --> NEW conjugate Polak-Ribiere (homemade)";
  EndIf
-  If (FLAG_MINMETHOD == 333 && FLAG_VARORDIFF==1)
+  If (FLAG_MINMETHOD == 333 && FLAG_APPROACH==1)
  Label_FLAG_MINMETHOD ="FLAG_MINMETHOD=11    --> NEW conjugate Polak-Ribiere+ (homemade)";
  EndIf
-  If (FLAG_MINMETHOD == 1999 && FLAG_VARORDIFF==1)
+  If (FLAG_MINMETHOD == 1999 && FLAG_APPROACH==1)
  Label_FLAG_MINMETHOD ="FLAG_MINMETHOD=11    --> NEW conjugate Dai Yuan 1999 (p.85) (homemade)";
  EndIf
-  If (FLAG_MINMETHOD == 2005 && FLAG_VARORDIFF==1)
+  If (FLAG_MINMETHOD == 2005 && FLAG_APPROACH==1)
  Label_FLAG_MINMETHOD ="FLAG_MINMETHOD=11    --> NEW conjugate Hager Zhang 2005 (p.161) (homemade)";
  EndIf
-  If (FLAG_MINMETHOD == 77 && FLAG_VARORDIFF==1)
+  If (FLAG_MINMETHOD == 77 && FLAG_APPROACH==1)
  Label_FLAG_MINMETHOD ="FLAG_MINMETHOD=11    --> NEW newton (homemade)";
  EndIf
-  If (FLAG_MINMETHOD == 2 && FLAG_VARORDIFF==1)
+  If (FLAG_MINMETHOD == 2 && FLAG_APPROACH==1)
  Label_FLAG_MINMETHOD ="FLAG_MINMETHOD=2     --> conjugate fr (gsl)";
  EndIf
-  If (FLAG_MINMETHOD == 3 && FLAG_VARORDIFF==1)
+  If (FLAG_MINMETHOD == 3 && FLAG_APPROACH==1)
  Label_FLAG_MINMETHOD ="FLAG_MINMETHOD=3     --> conjugate pr (gsl)";
  EndIf
-  If (FLAG_MINMETHOD == 4 && FLAG_VARORDIFF==1)
+  If (FLAG_MINMETHOD == 4 && FLAG_APPROACH==1)
  Label_FLAG_MINMETHOD ="FLAG_MINMETHOD=4     --> bfgs2 (gsl)";
  EndIf
-  If (FLAG_MINMETHOD == 5 && FLAG_VARORDIFF==1)
+  If (FLAG_MINMETHOD == 5 && FLAG_APPROACH==1)
  Label_FLAG_MINMETHOD ="FLAG_MINMETHOD=5     --> bfgs (gsl)";
  EndIf
-  If (FLAG_MINMETHOD == 6 && FLAG_VARORDIFF==1)
+  If (FLAG_MINMETHOD == 6 && FLAG_APPROACH==1)
  Label_FLAG_MINMETHOD ="FLAG_MINMETHOD=6     --> steepest descent (gsl)";
  EndIf
  If (FLAG_ANA ==0)
@@ -383,16 +346,16 @@ Label_FLAG_ANA="...";
  Printf["Number of cells      --> %g",N];
  Printf["k   omega    kappa    "];
  For iSub In {1:N}
-    Printf["%g   %g   %g   ",iSub,w~{iSub},chi~{iSub}];
+    Printf["%g   %g   %g   ",iSub,w~{iSub},kappa~{iSub}];
  EndFor
  Printf[""];
-  Printf[StrCat[Label_FLAG_TANORLANG]];
+  Printf[StrCat[Label_FLAG_ANHYLAW]];
  Printf["Ja   ha   Jb   hb"];
  Printf["%g   %g   %g   %g",Ja,ha,Jb,hb];
  Printf[""];
-  Printf[StrCat[Label_FLAG_VARORDIFF]];
+  Printf[StrCat[Label_FLAG_APPROACH]];
  Printf[StrCat[Label_FLAG_INVMETHOD]];
-  Printf[StrCat[Label_FLAG_ROOTFINDING1D]];
+  Printf[StrCat[Label_FLAG_JACEVAL]];
  Printf[StrCat[Label_FLAG_MINMETHOD]];
  Printf[StrCat[Label_FLAG_ANA]];
  Printf["FLAG_HOMO            --> %g",FLAG_HOMO];
@@ -402,7 +365,7 @@ Label_FLAG_ANA="...";
  Printf[""];
  Printf["TOLERANCE_NR         --> %g",TOLERANCE_NR];
  Printf["MAX_ITER_NR          --> %g",MAX_ITER_NR];
-  If(FLAG_VARORDIFF==1)
+  If(FLAG_APPROACH==1)
  Printf["TOLERANCE_NJ         --> %g",TOLERANCE_NJ];
  Printf["DELTA_0              --> %g",DELTA_0];
  Printf[""];

--- a/Team32/square.geo
+++ b/Team32/square.geo
@@ -36,6 +36,7 @@ Physical Line(15) = {1,2,3,4};
 Physical Surface(16)={6};
+Recursive Color LightGrey { Surface{ 6}; }
 // Post-processing point
@@ -51,3 +52,4 @@ For k In {1:num_postop_points}
 EndFor
 EndIf
--- a/Team32/square.pro
+++ b/Team32/square.pro
@@ -24,7 +24,7 @@ getdp square.pro -pos Get_AllTS
 // To launch computation and Then
 // do PostOp after the end of the simulation (create res/*.pos res/*dat)
-getdp square.pro -solve Analysis -pos Get_AllTS 
+getdp square.pro -v 3 -solve Analysis -pos Get_AllTS 
 //---------
@@ -67,7 +67,11 @@ Include "square_data.geo";
 // Output Directory
 //-------------------------------------------------------------------------
+If (Exists(simulabel))
+    Dir = Sprintf("res%g/",simulabel);
+Else
    Dir = 'res/';
+EndIf
 //Include "param_EnergHyst.dat";
 //Dir   = Sprintf("res_N%g_C%g/",N00,Flag_TestCase_00);
@@ -131,8 +135,8 @@ Function {
  hx[] = Frelax[] * (  hax     * Cos[2*Pi*Freq    *$Time]
                     + haharmx * Cos[2*Pi*freqharm*$Time]) ;
  hy[] = Frelax[] * (  hay0 
-                     + hay     * (Cos[phase]*Cos[2*Pi*Freq    *$Time] + Sin[phase]*Sin[2*Pi*Freq    *$Time])
+                     + hay     * (Cos[phase]*Cos[2*Pi*Freq    *$Time] - Sin[phase]*Sin[2*Pi*Freq    *$Time])
-                     + haharmy * (Cos[phase]*Cos[2*Pi*freqharm*$Time] + Sin[phase]*Sin[2*Pi*freqharm*$Time]) );
+                     + haharmy * (Cos[phase]*Cos[2*Pi*freqharm*$Time] - Sin[phase]*Sin[2*Pi*freqharm*$Time]) );
 /*
  hx[] = Frelax[] * (  hax     * Complex_MH[1,0]{Freq} 
                     + haharmx * Complex_MH[1,0]{freqharm}) ;

--- a/Team32/square_data.dat
+++ b/Team32/square_data.dat
+simulabel = 0;
+THICKNESS_00 = 1;
+NumEle_00 = 1;
+mm = 0.001;
+Flag_3Dmodel = 0;
+Flag_AnalysisType_00 = 1;
+freq_00 = 1;
+nstep_00 = 200;
+NbT_00 = 2;
+tmax_00 = 1;
+ha_00 = 150;
+Flag_TestCase_00 = 0;
+hax_00 = 150;
+haharmx_00 = 0;
+hay0_00 = 0;
+hay_00 = 150;
+haharmy_00 = 0;
+freqharm_00 = 0;
+phase_00 = (20/180)*3.141592653589793;
+Flag_NL_law00 = 4;
+Flag_NLRes00 = 1;
+Flag_UpdateSeparated= 1;
+Flag_ItLoopNRes00 = 31;
+Flag_ItLoopNDoFirstIter00 = 0;
+Nb_max_iter00 = 50;
+stop_criterion00 = 1e-05;
+relaxation_factor00 = 1;
+Flag_AdaptRelax00 = 1;
+relax_max00 = 1;
+relax_min00 = 0.1;
+relax_numtest00 = 10;
+TestAllFactors00 = 1;
+Flag_AdaptRelax_Version = 2;
+Flag_ExtrapolationOrder = 0;
--- a/Team32/square_data.geo
+++ b/Team32/square_data.geo
 //-------------------------------------------------------------------------
 // Default Parameters
 //-------------------------------------------------------------------------
+LOAD_DATA_FILE=0;
+If (LOAD_DATA_FILE && (FileExists("square_data.dat")))
+        Printf('square_data.dat is used');
+        Include "square_data.dat";
+Else
 // GEO parameters .............................
 THICKNESS_00 = 1; // thickness of the sample
@@ -14,13 +19,13 @@ Flag_AnalysisType_00 = 1; // 0: "static"
                          // 1: "time domain"
                          // 2: "frequency domain"
 freq_00 = 1;    // frequency (no impact as there is no dynamic effect)
-nstep_00 = 200; // number of time step per period (without Flexible DT)
+nstep_00 = 100; // number of time step per period (without Flexible DT)
-NbT_00 = 2; //number of periods
+NbT_00 = 1; //number of periods
 tmax_00 = 1; // final time of simulation
 // SOURCE parameters ............................
 ha_00 = 150; // maximal amplitude of the imposed magnetic field
-Flag_TestCase_00 = 1; // 0: CASE 0 (Default)
+Flag_TestCase_00 = 6; // 0: CASE 0 (Default)
                      // 1: CASE 1 (1Dx)
                      // 2: CASE 2 (1Dx+harm)
                      // 3: CASE 3 (1Dxy)
@@ -45,25 +50,44 @@ Flag_NL_law00 = 4;  // 0: "linear",
                    // 3: "Jiles-Atherton hysteresis model",
                    // 4: "EnergHyst model"
-Flag_NLRes00 = 0; // 0: use classical IterativeLoop to solve the NL (non linear) problem
+Flag_ExtrapolationOrder = 0;  // degree 0 ==> initialization at the previous Time Solution
+                              // degree 1 ==> linear extrapolation from the 2 previous Time Solutions
+                              // degree 2 ==> quadratic extrapolation from the 3 previous Time Solutions
+Flag_NLRes00 = 2; // 0: use classical IterativeLoop to solve the NL (non linear) problem
                  // 1: use IterativeLoopN to solve the NL (non linear) problem
                  // 2: solve the NL (non linear) problem "by hand"
 // non linear loop default parameters
+Flag_UpdateSeparated= 1;
 Nb_max_iter00       = 50; // maximum number of NL (non linear) iterations
 stop_criterion00    = 1e-5;  // strop criterion for the NL (non linear) iteration
 relaxation_factor00 = 1; // default relaxation factor (between ]0,1] )
 Flag_AdaptRelax00   = 1; // set to 1 to test different relaxation factors
+Flag_AdaptRelax_Version = 2; // allow to chose between MySolveJac_AdaptRelax or MySolveJac_AdaptRelax2 in 'MyItLoopPro.pro' 
+                             // (used when Flag_NLRes00==2 and Flag_AdaptRelax00==1)
 relax_max00         = 1 ;
 relax_min00         = 0.1;
 relax_numtest00     = 10;
-TestAllFactors00    = 2;  // 0 : try first relaxation factors (from the list) and stops when the residual goes up 
+TestAllFactors00    = 1;  // 0 : try first relaxation factors (from the list) and stops when the residual goes up 
                          // 1 : try every relaxation factors (from the list) and keep the optimal one
                          // 2 : find the maximum relaxation factor that decreases the residual:
                          //     - start with the relaxation factor from the previous time step or from the previous iteration
                          //     - the relaxation factor is multiplied by a ratio as long as the residual decreases
                          //     - the relaxation factor is decreased by a ratio until a decreasing residual is found
                          // [3 : Build a parabola based on three relaxation factors and deduce a minimizing relaxation factor (NOT WORKING!!)]
+Flag_ItLoopNRes00 = 1;// 0  : Solution (dx=x-xp; x=x)
+                       // 1  : Residual (dx=res=b-Ax; x=b) (default for square ???)
+                       // 2  : RecalcResidual (dx=res=b-Ax; x=b) 
+                       // 31 : PostOp unknown field (az or f)
+                       // 32 : PostOp b field
+                       // 33 : PostOp h field (default for t32)
+Flag_ItLoopNDoFirstIter00 = 0;// set to 1 if one wants to force the first NL iteration before evaluating the error with IterativeLoopN
+                              // helpfull for square.pro with EnergHyst or JA because PostOp h field = 0 at iter 1, which stops IterativeLoopN prematurely
+EndIf
 Reltol_Mag00        = stop_criterion00; // 0 before with IterativeLoopN
 Abstol_Mag00        = stop_criterion00;
 Reltoldx_Mag00      = 1e-5*stop_criterion00;
@@ -86,6 +110,13 @@ NLRES_ITERATIVELOOP    = 0;
 NLRES_ITERATIVELOOPN   = 1;
 NLRES_ITERATIVELOOPPRO = 2;
+NLITLOOPN_SOLUTION       = 0;
+NLITLOOPN_RESIDUAL       = 1;
+NLITLOOPN_RECALCRESIDUAL = 2;
+NLITLOOPN_POSTOPX        = 31;
+NLITLOOPN_POSTOPB        = 32;
+NLITLOOPN_POSTOPH        = 33;
 AN_STATIC    = 0;
 AN_TIME      = 1;
 AN_FREQUENCY = 2;
@@ -183,6 +214,23 @@ DefineConstant[
      Help Str["- Use 'IterativeLoop' to solve the NL (non linear) problem with classical IterativeLoop Operation",
               "- Use 'IterativeLoopN' to solve the NL (non linear) problem with IterativeLoopN Operation",
               "- Use 'IterativeLoopPro' to solve the NL (non linear) problem 'by hand'"]},
+  Flag_ItLoopNRes = { (Flag_NLRes!=NLRES_ITERATIVELOOPN)? NLITLOOPN_SOLUTION:Flag_ItLoopNRes00  , Choices {
+      NLITLOOPN_SOLUTION ="Solution",
+      NLITLOOPN_RESIDUAL ="RecalcResidual",
+      NLITLOOPN_RECALCRESIDUAL ="Residual",
+      NLITLOOPN_POSTOPX ="PostOp Potential",
+      NLITLOOPN_POSTOPB ="PostOp b field",
+      NLITLOOPN_POSTOPH ="PostOp h field"}, 
+    Name StrCat[ppNL, "1Error Calculation based on ..."], Highlight Str[colNL1], Visible Flag_NL ,ReadOnly (Flag_NLRes!=NLRES_ITERATIVELOOPN),
+      Help Str["- the solution (dx=x-xp; x=x)",
+               "- the residual (dx=res=b-Ax; x=b)",
+               "- the recalculated residual (dx=res=b-Ax; x=b)",
+               "- the post-operated potential field (az or phi)",
+               "- the post-operated flux density field (b)",
+               "- the post-operated magnetic field (h)"]},
+  Flag_ItLoopNDoFirstIter = {Flag_ItLoopNDoFirstIter00, Choices{0,1}, Name StrCat[ppNL, "1First iteration done before error calculation"], 
+    Highlight Str[colNL1], Visible (Flag_NL && Flag_NLRes==NLRES_ITERATIVELOOPN),
+    Help Str["Set to 1 to force the first NL iteration to be done before evaluating the error with IterativeLoopN"]},
  Nb_max_iter = {Nb_max_iter00, Name StrCat[ppNL, "1Nb max iter"], 
    Highlight Str[colNL4], Visible Flag_NL}

--- a/Team32/t32.geo
+++ b/Team32/t32.geo
@@ -19,10 +19,10 @@ lc2  = wcore/5;
 /*
-lc0=lc0/2;
+lc0=lc0/3;
-lc1=lc1/2;
+lc1=lc1/3;
-lc2=lc2/2;
+lc2=lc2/3;
-*/
+//*/
 lcri = lc0*4; // Pi*Rint/30;
 lcro = lc0*4; // Pi*Rext/30;
@@ -271,32 +271,32 @@ If(Flag_3Dmodel==0)
  //=================================================
  // Physical regions for FE analysis (2D)
  //=================================================
-    Physical Surface("Core", CORE)   = surf_ECore[];
+  Physical Surface(CORE)   = surf_ECore[];
  nb_surf_coil = #surf_Coil[];
-  Physical Surface("Coil right (left side)", COIL)   = surf_Coil[{0:nb_surf_coil-1:4}];
+  Physical Surface(COIL)   = surf_Coil[{0:nb_surf_coil-1:4}];
-  Physical Surface("Coil right (right side)", COIL+1) = surf_Coil[{1:nb_surf_coil-1:4}];
+  Physical Surface(COIL+1) = surf_Coil[{1:nb_surf_coil-1:4}];
-  Physical Surface("Coil left (right side)", COIL+2) = surf_Coil[{2:nb_surf_coil-1:4}];
+  Physical Surface(COIL+2) = surf_Coil[{2:nb_surf_coil-1:4}];
-  Physical Surface("Coil left (left side)", COIL+3) = surf_Coil[{3:nb_surf_coil-1:4}];
+  Physical Surface(COIL+3) = surf_Coil[{3:nb_surf_coil-1:4}];
-  Physical Surface("Air", AIR)    = surf_Air[];
+  Physical Surface(AIR)    = surf_Air[];
-  Physical Surface("Air infinity shell", AIRINF) = surf_AirInf[];
+  Physical Surface(AIRINF) = surf_AirInf[];
-  Physical Line("Outer boundary", SURF_AIROUT) = ln_bnd[];
+  Physical Line(SURF_AIROUT) = ln_bnd[];
  //=================================================
  // Some colors... for aesthetics :-)
  //=================================================
  Reverse Surface {surf_ECore[{2,3}], surf_Coil[{2,3}],surf_Air[{2,3}],surf_AirInf[1]};
-  Recursive Color SkyBlue { Surface{surf_Air[]};}
+  Recursive Color White { Surface{surf_Air[]};}
-  Recursive Color Blue { Surface{surf_AirInf[]};}
+  Recursive Color White { Surface{surf_AirInf[]};}
-  Recursive Color SteelBlue { Surface{ surf_ECore[]}; }
+  Recursive Color LightGrey { Surface{ surf_ECore[]}; }
-  Recursive Color Red    { Surface{surf_Coil[{0,1}]}; }
+  Recursive Color LightBlue    { Surface{surf_Coil[{0,1}]}; }
-  Recursive Color Yellow { Surface{surf_Coil[{2,3}]}; }
+  Recursive Color Pink { Surface{surf_Coil[{2,3}]}; }
  If(!Flag_Symmetry)
    Reverse Surface {surf_ECore[{4,5}], surf_Coil[{4,5}],surf_Air[{4,5}],surf_AirInf[2]};
-    Recursive Color Red    { Surface{surf_Coil[{4,5}]}; }
+    Recursive Color LightBlue    { Surface{surf_Coil[{4,5}]}; }
-    Recursive Color Yellow { Surface{surf_Coil[{6,7}]}; }
+    Recursive Color Pink { Surface{surf_Coil[{6,7}]}; }
  EndIf
 EndIf
@@ -338,7 +338,7 @@ If(Flag_3Dmodel==1)
  bnd_CoilR_plus_hole[] = Abs(CombinedBoundary{Volume{vol_CoilR[], vol_inCoilR[]};});
  bnd_inCoilR[]  = Abs(CombinedBoundary{Volume{vol_inCoilR[]};});
  bnd_inCoilR[] -= bnd_CoilR_plus_hole[];
-Physical Surface("Skin Coil right- only around cut", SKINCOIL_CUT+0) = bnd_inCoilR[];
+  Physical Surface("Skin Coil right - only cut",SKINCOIL_CUT+0) = bnd_inCoilR[];
  bnd_CoilR_plus_hole[] -= Abs(CombinedBoundary{Volume{vol_CoilR[]};});
  bnd_CoilR_plus_hole[] -= surf_cut_xy[];
@@ -362,7 +362,7 @@ Physical Surface("Skin Coil right- only around cut", SKINCOIL_CUT+0) = bnd_inCoi
  bnd_CoilL_plus_hole[] = Abs(CombinedBoundary{Volume{vol_CoilL[], vol_inCoilL[]};});
  bnd_inCoilL[]  = Abs(CombinedBoundary{Volume{vol_inCoilL[]};});
  bnd_inCoilL[] -= bnd_CoilL_plus_hole[];
-  Physical Surface("Skin Coil left- only around cut", SKINCOIL_CUT+1) = bnd_inCoilL[];
+  Physical Surface("Skin Coil left - only cut",SKINCOIL_CUT+1) = bnd_inCoilL[];
  bnd_CoilL_plus_hole[] -= Abs(CombinedBoundary{Volume{vol_CoilL[]};});
  bnd_CoilL_plus_hole[] -= surf_cut_xy[];
@@ -373,8 +373,8 @@ Physical Surface("Skin Coil right- only around cut", SKINCOIL_CUT+0) = bnd_inCoi
  EndIf
  Physical Volume("Air",AIR)         = vol_Air[];
-  Physical Surface("Symmetry cut XY", SURF_CUTXY) = surf_cut_xy[];
+  Physical Surface("cut XY", SURF_CUTXY) = surf_cut_xy[];
-  Physical Surface("Symmetry cut XZ", SURF_CUTXZ) = surf_cut_xz[];
+  Physical Surface("cut XZ", SURF_CUTXZ) = surf_cut_xz[];
  all_vol[] = Volume '*';
  bnd_all[] = Abs(CombinedBoundary{Volume{all_vol[]};});
@@ -385,14 +385,14 @@ Physical Surface("Skin Coil right- only around cut", SKINCOIL_CUT+0) = bnd_inCoi
  //=================================================
  // Some colors... for aesthetics :-)
  //=================================================
-  Recursive Color SkyBlue { Volume{vol_Air[]};}
+  Recursive Color White { Volume{vol_Air[]};}
-  Recursive Color SteelBlue { Volume{ vol_Core[]}; }
+  Recursive Color LightGrey { Volume{ vol_Core[]}; }
  Recursive Color Red    { Volume{vol_CoilR[]}; }
  Recursive Color Yellow { Volume{vol_CoilL[]}; }
 EndIf
+aaa=newp;
 // Post-processing point
 For k In {1:num_postop_points}
  Point(newp) = {xpos~{k},ypos~{k}, 0 }; // for visu

--- a/Team32/t32.pro
+++ b/Team32/t32.pro
@@ -25,6 +25,7 @@ getdp t32.pro -pos Get_AllTS
 // To launch computation and Then
 // do PostOp after the end of the simulation (create res/*.pos res/*dat)
 getdp t32.pro -solve Analysis -pos Get_AllTS 
+getdp t32.pro -solve Analysis -pos Get_AllTS_3D
 //---------
@@ -157,7 +158,11 @@ Function {
    xc2_2 = -xc1_2;
    SurfCoil[] = SurfaceArea[]{ELECCOIL+0} ;
+/*
+    vDir[Ind_2] = Vector[0,0,1];
+    vDir[Ind_1] = Vector[0,0,-1];
+*/
+///*
    vDir[Ind_2] = (Fabs[Z[]]>=Lz/2 && Fabs[X[]-x_1]<=wcore/2 ) ? Vector[-1,0,0]:(
      (Fabs[Z[]]<Lz/2)? Vector[0, 0, (X[]>x_1)?-1:1 ]:
      -Vector[ -Sin[Atan2[Z[]-Lz/2, X[]-((X[]>x_1)?xc1_1:xc1_2) ]#1], 0, Cos[#1]] ) ;
@@ -165,7 +170,7 @@ Function {
    vDir[Ind_1] = (Fabs[Z[]]>=Lz/2 && Fabs[X[]-x_2]<=wcore/2 ) ? Vector[-1,0,0]:(
      (Fabs[Z[]]<Lz/2)? Vector[0, 0, (X[]>x_2)?-1:1 ]:
      -Vector[ -Sin[Atan2[Z[]-Lz/2, X[]-((X[]>x_2)?xc2_2:xc2_1) ]#1], 0, Cos[#1]] ) ;
+//*/
  EndIf
  js0[] = NbWires[]/SurfCoil[]*vDir[] ;
@@ -215,8 +220,8 @@ Function {
  //pA = (Flag_AnalysisType==AN_STATIC) ? Pi/2: 0.; // Not Used
  //IA[] = F_Sin_wt_p[]{2*Pi*Freq, pA}; // DomainB
-  I1[] = ((Flag_NL_law==NL_JA ||  Flag_NL_law==NL_ENERGHYST)? Frelax[]:1) * (F_Sin_wt_p[]{2*Pi*Freq, Phase_1}+V5[]);
+  I1[] = ((Flag_AnalysisType==AN_TIME)? Frelax[]:1) * (F_Sin_wt_p[]{2*Pi*Freq, Phase_1}+V5[]);
-  I2[] = ((Flag_NL_law==NL_JA ||  Flag_NL_law==NL_ENERGHYST)? Frelax[]:1) * (F_Sin_wt_p[]{2*Pi*Freq, Phase_2}+V5[]);
+  I2[] = ((Flag_AnalysisType==AN_TIME)? Frelax[]:1) * (F_Sin_wt_p[]{2*Pi*Freq, Phase_2}+V5[]);
 }

--- a/Team32/t32_data.geo
+++ b/Team32/t32_data.geo
@@ -27,20 +27,30 @@ Flag_TestCase00 = 3; //1: - case 1: windings series connected with sinusoidal su
                     //4: - case 4: secondary winding on load  
 // RESOLUTION parameters .........................
-Flag_NL_law00 = 4;  // 0: "linear",
+Flag_NL_law00 = 1;  // 0: "linear",
                    // 1: "analytical",
                    // 2: "anhysteretic part of JA-model",
                    // 3: "Jiles-Atherton hysteresis model",
                    // 4: "EnergHyst model"
+Flag_ExtrapolationOrder = 2;  // degree 0 ==> initialization at the previous Time Solution
+                              // degree 1 ==> linear extrapolation from the 2 previous Time Solutions
+                              // degree 2 ==> quadratic extrapolation from the 3 previous Time Solutions
 // non linear loop default parameters
 Flag_NLRes00 = 1; // 0: use classical IterativeLoop to solve the NL (non linear) problem
                  // 1: use IterativeLoopN to solve the NL (non linear) problem
                  // 2: solve the NL (non linear) problem "by hand"
+Flag_DomainLam00= 0;
+Flag_UpdateSeparated= 1;
 Nb_max_iter00       = 50; // maximum number of NL (non linear) iterations
 stop_criterion00    = 1e-5;  // strop criterion for the NL (non linear) iteration
 relaxation_factor00 = 1; // default relaxation factor (between ]0,1] )
 Flag_AdaptRelax00   = 1; // set to 1 to test different relaxation factors
+Flag_AdaptRelax_Version = 2; // allow to chose between MySolveJac_AdaptRelax or MySolveJac_AdaptRelax2 in 'MyItLoopPro.pro' 
+                             // (used when Flag_NLRes00==2 and Flag_AdaptRelax00==1)
 relax_max00         = 1 ;
 relax_min00         = 0.1;
 relax_numtest00     = 10;
@@ -51,6 +61,15 @@ TestAllFactors00    = 1;  // 0 : try first relaxation factors (from the list) an
                          //     - the relaxation factor is multiplied by a ratio as long as the residual decreases
                          //     - the relaxation factor is decreased by a ratio until a decreasing residual is found
                          // [3 : Build a parabola based on three relaxation factors and deduce a minimizing relaxation factor (NOT WORKING!!)]
+Flag_ItLoopNRes00 = 33;// 0  : Solution (dx=x-xp; x=x)
+                       // 1  : Residual (dx=res=b-Ax; x=b) (default for square ???)
+                       // 2  : RecalcResidual (dx=res=b-Ax; x=b) 
+                       // 31 : PostOp unknown field (az or f)
+                       // 32 : PostOp b field
+                       // 33 : PostOp h field (default for t32)
+Flag_ItLoopNDoFirstIter00 = 0;// set to 1 if one wants to force the first NL iteration before evaluating the error with IterativeLoopN
+                              // helpfull for square.pro with EnergHyst or JA because PostOp h field = 0 at iter 1, which stops IterativeLoopN prematurely
 Reltol_Mag00        = stop_criterion00; // 0 before with IterativeLoopN
 Abstol_Mag00        = stop_criterion00;
 Reltoldx_Mag00      = 1e-5*stop_criterion00;
@@ -73,6 +92,13 @@ NLRES_ITERATIVELOOP    = 0;
 NLRES_ITERATIVELOOPN   = 1;
 NLRES_ITERATIVELOOPPRO = 2;
+NLITLOOPN_SOLUTION       = 0;
+NLITLOOPN_RESIDUAL       = 1;
+NLITLOOPN_RECALCRESIDUAL = 2;
+NLITLOOPN_POSTOPX        = 31;
+NLITLOOPN_POSTOPB        = 32;
+NLITLOOPN_POSTOPH        = 33;
 AN_STATIC    = 0;
 AN_TIME      = 1;
 AN_FREQUENCY = 2;
@@ -129,6 +155,8 @@ DefineConstant[
  //Printf("====> SymmetryFactor=%g", SymmetryFactor);
 ];
+  Printf("====> SymmetryFactor=%g", SymmetryFactor);
 // Geometric dimensions........................ (subMenu ppDIM)
 close_menu = 1;
 DefineConstant[
@@ -144,7 +172,9 @@ DefineConstant[
    Highlight Str[colorro]},
  hcore = {hcoil+2*wcore, Name StrCat[ppDIM, "2Core height [m]"], ReadOnly 1,
    Highlight Str[colorro]},
-  AxialLength = {cteaux*2.4*mm, Name StrCat[ppDIM, "0Length along z-axis [m]"], Highlight Str[colorpp]}
+  AxialLength = {cteaux*2.4*mm, Name StrCat[ppDIM, "0Length along z-axis [m]"], Highlight Str[colorpp]},
+  Thickness = {0.48*mm,  Name StrCat[ppDIM, "1Thickness of lamination [m]"],
+    Highlight Str[colorpp],Closed close_menu}
 ];
 // Excitation Source .......................... (SubMenu ppTC)
@@ -183,6 +213,8 @@ DefineConstant[
    Highlight Str[colAC2],Visible (Flag_AnalysisType==AN_TIME)}
  NbSteps = {nstep_00, Name StrCat[ppAC, "Number of steps"], 
    Highlight Str[colAC2], Visible (Flag_AnalysisType==AN_TIME)}
+  Flag_DomainLam = {Flag_DomainLam00, Choices{0,1}, 
+    Name  StrCat[ppAC, "1Activate DomainLam"], Visible 1}
 ];
 // Material Law ............................... (SubMenu ppML)
@@ -209,6 +241,25 @@ DefineConstant[
      Help Str["- Use 'IterativeLoop' to solve the NL (non linear) problem with classical IterativeLoop Operation",
               "- Use 'IterativeLoopN' to solve the NL (non linear) problem with IterativeLoopN Operation",
               "- Use 'IterativeLoopPro' to solve the NL (non linear) problem 'by hand'"]},
+  Flag_ItLoopNRes = { (Flag_NLRes!=NLRES_ITERATIVELOOPN)? NLITLOOPN_SOLUTION:Flag_ItLoopNRes00  , Choices {
+      NLITLOOPN_SOLUTION ="Solution",
+      NLITLOOPN_RESIDUAL ="RecalcResidual",
+      NLITLOOPN_RECALCRESIDUAL ="Residual",
+      NLITLOOPN_POSTOPX ="PostOp Potential",
+      NLITLOOPN_POSTOPB ="PostOp b field",
+      NLITLOOPN_POSTOPH ="PostOp h field"}, 
+    Name StrCat[ppNL, "1Error Calculation based on ..."], Highlight Str[colNL1], Visible Flag_NL ,ReadOnly (Flag_NLRes!=NLRES_ITERATIVELOOPN),
+      Help Str["- the solution (dx=x-xp; x=x)",
+               "- the residual (dx=res=b-Ax; x=b)",
+               "- the recalculated residual (dx=res=b-Ax; x=b)",
+               "- the post-operated potential field (az or phi)",
+               "- the post-operated flux density field (b)",
+               "- the post-operated magnetic field (h)"]},
+  Flag_ItLoopNDoFirstIter = {Flag_ItLoopNDoFirstIter00, Choices{0,1}, Name StrCat[ppNL, "1First iteration done before error calculation"], 
+    Highlight Str[colNL1], Visible (Flag_NL && Flag_NLRes==NLRES_ITERATIVELOOPN),
+    Help Str["Set to 1 to force the first NL iteration to be done before evaluating the error with IterativeLoopN"]},
  Nb_max_iter = {Nb_max_iter00, Name StrCat[ppNL, "1Nb max iter"], 
    Highlight Str[colNL4], Visible Flag_NL}
  stop_criterion = {stop_criterion00, Name StrCat[ppNL, "2stop criterion"], 
@@ -276,7 +327,8 @@ mu0 = 4.e-7 * Pi ;
 sigma_al = 3.72e7 ; // conductivity of aluminum [S/m]
 sigma_cu = 5.77e7  ; // conductivity of copper [S/m]
-sigma_core = 2.5e7/1000;
+//sigma_core = 2.5e7/1000;
+sigma_core = 1.78e6; // for team 32
 sigma_coil = 5.9e7;
 mur_fe     = 1000; // linear case