diff --git a/Conductors3D/w2d.pro b/Conductors3D/w2d.pro
index f43b73edcc2904a974724a4765d6df47a223925d..b09e70de2d01273f532f581113d9e2207062460d 100644
--- a/Conductors3D/w2d.pro
+++ b/Conductors3D/w2d.pro
@@ -36,13 +36,13 @@ Group{
   LWIRES = Region[ {} ];
   For i In {1:NumWires}
     LWIRE~{i} = Region[ (50+i) ];
-	LWIRES += Region[ LWIRE~{i} ];
-	If( Flag_Thin == 0 )
+    LWIRES += Region[ LWIRE~{i} ];
+    If( Flag_Thin == 0 )
       VWIRE~{i} = Region[ (10+i) ];
       VWIRES += Region[ VWIRE~{i} ];
-	  SWIRE~{i} = Region[ {} ];
-	Else
-	  VWIRE~{i} = Region[ {} ];
+      SWIRE~{i} = Region[ {} ];
+    Else
+      VWIRE~{i} = Region[ {} ];
       SWIRE~{i} = ElementsOf[ AIR, OnOneSideOf LWIRE~{i} ];
     EndIf
   EndFor
@@ -100,9 +100,9 @@ Function{
   For i In {1:NumWires}
     // Current density
     J[ Region[ {VWIRE~{i}, LWIRE~{i}} ] ] = 
-	  Vector[ 0, 0, WI~{i}/A_c * Exp[ i[]*WP~{i}*deg ] ];
+    Vector[ 0, 0, WI~{i}/A_c * Exp[ i[]*WP~{i}*deg ] ];
 
-	// local radial coordinate with origin on wire i
+    // local radial coordinate with origin on wire i
     R~{i}[] = Hypot[ X[]-WX~{i} , Y[]-WY~{i} ];
   EndFor
 
@@ -111,10 +111,10 @@ Function{
 
   //radial analytical solutions with a zero flux condition imposed at $1(=r)=rs
   Analytic_A[] = // per Amp
-	(($1>rw) ?  mu0/(2*Pi) * Log[rs/$1] :
-				mu0/(2*Pi) * Log[rs/rw] + i[]*( wI[$1]-wI[rw] )/(sigma*omega) ); 
+                (($1>rw) ?  mu0/(2*Pi) * Log[rs/$1] :
+                  mu0/(2*Pi) * Log[rs/rw] + i[]*( wI[$1]-wI[rw] )/(sigma*omega) ); 
   AnalyticStatic_A[] = mu0/(2*Pi) * // per Amp
-	(($1>rw) ? Log[rs/$1] : Log[rs/rw] + mur*(1-($1/rw)^2)/2 );
+                (($1>rw) ? Log[rs/$1] : Log[rs/rw] + mur*(1-($1/rw)^2)/2 );
  
   // Impedance of thin wire p.u. length
   R_DC = 1./(sigma*A_c);
@@ -136,7 +136,7 @@ Function{
     + WI_3*Analytic_A[R_3[]] + (WI_1+WI_2+WI_3)*mu0/(2*Pi) * NumWires * Log[rout/rs] ;
     Correction_A[] = ((R_1[]<rs) ? WI_1*Analytic_A[R_1[]] : 
                      ((R_2[]<rs) ? WI_2*Analytic_A[R_2[]] : 
-					 ((R_3[]<rs) ? WI_3*Analytic_A[R_3[]] : 0 )));
+                     ((R_3[]<rs) ? WI_3*Analytic_A[R_3[]] : 0 )));
   EndIf
 }
 
@@ -155,14 +155,14 @@ Integration {
     Case {
       { Type Gauss ;
         Case {
-	  { GeoElement Point       ; NumberOfPoints  1 ; }
-	  { GeoElement Line        ; NumberOfPoints  3 ; }
-	  { GeoElement Triangle    ; NumberOfPoints  3 ; }
-	  { GeoElement Quadrangle  ; NumberOfPoints  4 ; }
-	  { GeoElement Prism       ; NumberOfPoints 21 ; }
-	  { GeoElement Tetrahedron ; NumberOfPoints  4 ; }
-	  { GeoElement Hexahedron  ; NumberOfPoints  6 ; }
-	  { GeoElement Pyramid     ; NumberOfPoints  8 ; }
+          { GeoElement Point       ; NumberOfPoints  1 ; }
+          { GeoElement Line        ; NumberOfPoints  3 ; }
+          { GeoElement Triangle    ; NumberOfPoints  3 ; }
+          { GeoElement Quadrangle  ; NumberOfPoints  4 ; }
+          { GeoElement Prism       ; NumberOfPoints 21 ; }
+          { GeoElement Tetrahedron ; NumberOfPoints  4 ; }
+          { GeoElement Hexahedron  ; NumberOfPoints  6 ; }
+          { GeoElement Pyramid     ; NumberOfPoints  8 ; }
         }
       }
     }
@@ -182,7 +182,7 @@ Constraint{
     Case {
       For i In {1:NumWires}
         { Region VWIRE~{i} ; Value  WI~{i} ; }
-	    { Region LWIRE~{i} ; Value  WI~{i} ; }
+        { Region LWIRE~{i} ; Value  WI~{i} ; }
       EndFor
     }
   }
@@ -206,9 +206,9 @@ FunctionSpace {
   { Name Hcurl_a_2D; Type Form1P;
     BasisFunction {
       { Name sc; NameOfCoef ac; Function BF_PerpendicularEdge;
-		Support Dom_Hcurl_a; Entity NodesOf[All]; }
+        Support Dom_Hcurl_a; Entity NodesOf[All]; }
       { Name sw; NameOfCoef aw; Function BF_PerpendicularEdge;
-		Support Dom_Hthin_a; Entity NodesOf[lVol_C]; }
+        Support Dom_Hthin_a; Entity NodesOf[lVol_C]; }
     }
     SubSpace {
       { Name ac ; NameOfBasisFunction { sc } ; }
@@ -216,7 +216,7 @@ FunctionSpace {
     }
     Constraint {
       { NameOfCoef ac;
-	EntityType Auto; NameOfConstraint Hcurl_a_2D_ac; }
+        EntityType Auto; NameOfConstraint Hcurl_a_2D_ac; }
     }
   }
 
@@ -253,9 +253,9 @@ FunctionSpace {
   { Name Hcurl_a_3D; Type Form1;
     BasisFunction {
       { Name sc; NameOfCoef ac; Function BF_Edge;
-	Support Dom_Hcurl_a; Entity EdgesOf[All]; }
+        Support Dom_Hcurl_a; Entity EdgesOf[All]; }
       { Name sw; NameOfCoef aw; Function BF_Edge;
-	Support Dom_Hthin_a; Entity EdgesOf[lVol_C]; }
+        Support Dom_Hthin_a; Entity EdgesOf[lVol_C]; }
     }
     SubSpace {
       { Name ac ; NameOfBasisFunction { sc } ; }
@@ -263,16 +263,16 @@ FunctionSpace {
     }
     Constraint {
       { NameOfCoef ac;
-	EntityType Auto; NameOfConstraint Hcurl_a_3D_ac; }
+        EntityType Auto; NameOfConstraint Hcurl_a_3D_ac; }
     }
   }
 
   { Name Hgrad_u_3D; Type Form0; 
     BasisFunction {
       { Name sn; NameOfCoef un; Function BF_Node;
-	Support Dom_Hgrad_u; Entity NodesOf[ Vol_C, Not Electrodes]; }
+        Support Dom_Hgrad_u; Entity NodesOf[ Vol_C, Not Electrodes]; }
       { Name sn2; NameOfCoef un2; Function BF_GroupOfNodes;
-	Support Dom_Hgrad_u; Entity GroupsOfNodesOf[ Electrodes ]; }  
+        Support Dom_Hgrad_u; Entity GroupsOfNodesOf[ Electrodes ]; }  
    }
    GlobalQuantity {
       { Name U; Type AliasOf       ; NameOfCoef un2; }
@@ -287,7 +287,6 @@ FunctionSpace {
 
 Formulation {
   { Name MagnetoDynamics; Type FemEquation;
-  If( !Flag_3D )
     Quantity {
       { Name a;  Type Local;  NameOfSpace Hcurl_a_2D; }
       { Name ac; Type Local;  NameOfSpace Hcurl_a_2D[ac]; }
@@ -331,56 +330,19 @@ Formulation {
       GlobalTerm { [ Dof{I} , {U} ]; In Vol_C; }
 
       GlobalTerm { [ Dof{Us} , {Is} ]; 
-		In lVol_C; }
-	  GlobalTerm { [ Analytic_R[] * Dof{Is} , {Is} ]; 
-	  	In lVol_C;}
+        In lVol_C; }
+      GlobalTerm { [ Analytic_R[] * Dof{Is} , {Is} ]; 
+        In lVol_C;}
       Integral { DtDof [ Dof{ac}/A_c , {i} ];
         In lVol_C; Jacobian Vol; Integration I1; }
+
       If( Flag_Thin != 3)
         Integral { DtDof [ -Dof{aw}/A_c , {i} ];
-		  In lVol_C; Jacobian Vol; Integration I1;}
-	    GlobalTerm { DtDof [ Analytic_L[] * Dof{Is} , {Is} ]; 
-	  	  In lVol_C;}
+          In lVol_C; Jacobian Vol; Integration I1;}
+        GlobalTerm { DtDof [ Analytic_L[] * Dof{Is} , {Is} ]; 
+          In lVol_C;}
       EndIf
     }
-	Else // ###############    3D 
-    Quantity {
-      { Name a;  Type Local;  NameOfSpace Hcurl_a_3D; }
-      { Name ac; Type Local;  NameOfSpace Hcurl_a_3D[ac]; }
-      { Name aw; Type Local;  NameOfSpace Hcurl_a_3D[aw]; }
-      // { Name v;  Type Local;  NameOfSpace Hgrad_u_3D; }
-      // { Name U;  Type Global; NameOfSpace Hgrad_u_3D [U]; }
-      // { Name I;  Type Global; NameOfSpace Hgrad_u_3D [I]; }
-    }
-
-    Equation {
-    /*
-      Integral { [ nu[] * Dof{d ac} , {d ac} ];
-        In Vol_nu; Jacobian Vol; Integration I1; }
-      Integral { [ nu[] * Dof{d aw} , {d aw} ];
-        In Vol_nu; Jacobian Vol; Integration I1; }
-
-      Integral { [ -Vector[ 0, 0, I0/A_c] , {ac} ];
-        In lVol_C; Jacobian Vol; Integration I1; }
-      Integral { [ -Vector[ 0, 0, I0/A_c] , {aw} ];
-        In lVol_C; Jacobian Vol; Integration I1; }
-
-	  Integral { [ -Vector[ 0, 0, I0/A_c] , {a} ];
-		In Vol_C; Jacobian Vol; Integration I1; }
-
-      Integral { DtDof[ sigma[] * Dof{a} , {a} ];
-        In Vol_C; Jacobian Vol; Integration I1; }
-      Integral { [ sigma[] * Dof{d v} , {a} ];
-        In Vol_C; Jacobian Vol; Integration I1; }
-      Integral { DtDof[ sigma[] * Dof{a} , {d v} ];
-        In Vol_C; Jacobian Vol; Integration I1; }      
-      Integral { [ sigma[] * Dof{d v} , {d v} ];
-        In Vol_C; Jacobian Vol; Integration I1; }   
-      GlobalTerm { [Dof{U} , {I} ]; In Electrodes; }
-      */
-     
-    }
-  EndIf
   }
 }
 
@@ -388,126 +350,125 @@ Resolution {
  { Name Dynamic;
     System {
       { Name S; NameOfFormulation MagnetoDynamics;
-	Type ComplexValue; Frequency Freq; }
+        Type ComplexValue; Frequency Freq; }
     }
     Operation {
       CreateDir[DataDir];
 
       // SOLVER OPTION
       If(type_solver == S::SOLVER.gmres)
-		SetGlobalSolverOptions[
-							   "-pc_type ilu -ksp_type gmres -ksp_max_it 1000 -ksp_rtol 1e-6 -ksp_atol 1e-6"];
-      ElseIf(type_solver == S::SOLVER.fgmres)
-		SetGlobalSolverOptions[
-							   "-pc_type ilu -ksp_type fgmres -ksp_max_it 1000 -ksp_rtol 1e-6 -ksp_atol 1e-6"];
-      ElseIf(type_solver == S::SOLVER.bcgs)
-		SetGlobalSolverOptions[
-							   "-pc_type ilu -ksp_type bcgs -ksp_max_it 1000  -ksp_rtol 1e-6 -ksp_atol 1e-6"];
+        SetGlobalSolverOptions[
+          "-pc_type ilu -ksp_type gmres -ksp_max_it 1000 -ksp_rtol 1e-6 -ksp_atol 1e-6"];
+        ElseIf(type_solver == S::SOLVER.fgmres)
+        SetGlobalSolverOptions[
+          "-pc_type ilu -ksp_type fgmres -ksp_max_it 1000 -ksp_rtol 1e-6 -ksp_atol 1e-6"];
+        ElseIf(type_solver == S::SOLVER.bcgs)
+        SetGlobalSolverOptions[
+          "-pc_type ilu -ksp_type bcgs -ksp_max_it 1000  -ksp_rtol 1e-6 -ksp_atol 1e-6"];
       EndIf
-		SetGlobalSolverOptions["-petsc_prealloc 900"];
-
+      SetGlobalSolverOptions["-petsc_prealloc 900"];
+      
       Generate[S]; Solve[S]; SaveSolution[S]; 
 
-	  PostOperation[integaz];
-
-	  If( Flag_Maps == 1 )
-	    PostOperation[mapaz];
-	    PostOperation[cutaz];
-
-		DeleteFile["U_full.dat"];
-		DeleteFile["U_raw.dat"];
-		DeleteFile["U_reg.dat"];
-		DeleteFile["U_naive.dat"];
-		DeleteFile["joule_full.dat"];
-		DeleteFile["joule_raw.dat"];
-		DeleteFile["joule_reg.dat"];
-		DeleteFile["joule_naive.dat"];
-	  EndIf
-
-	  If( Flag_Thin != 0 )
-  	    For i In {1:NumWires}
-	       Print[ {i, rw, Sqrt[ $sarea~{i}/Pi], rs, $sarea~{i}, Pi*rs^2, skin_depth, $intby~{i}/$sarea~{i}}, 
-				 Format "WIRE %2g: rw = %7.3e rs = %7.3e (%7.3e) as = %7.3e (%7.3e) delta = %7.3e bave = %7.3e" ] ;
-		EndFor
-	  EndIf
+      PostOperation[integaz];
+
+      If( Flag_Maps == 1 )
+        PostOperation[mapaz];
+        PostOperation[cutaz];
+
+        DeleteFile["U_full.dat"];
+        DeleteFile["U_raw.dat"];
+        DeleteFile["U_reg.dat"];
+        DeleteFile["U_naive.dat"];
+        DeleteFile["joule_full.dat"];
+        DeleteFile["joule_raw.dat"];
+        DeleteFile["joule_reg.dat"];
+        DeleteFile["joule_naive.dat"];
+      EndIf
+
+      If( Flag_Thin != 0 )
+        For i In {1:NumWires}
+          Print[ {i, rw, Sqrt[ $sarea~{i}/Pi], rs, $sarea~{i}, Pi*rs^2, skin_depth, $intby~{i}/$sarea~{i}}, 
+            Format "WIRE %2g: rw = %7.3e rs = %7.3e (%7.3e) as = %7.3e (%7.3e) delta = %7.3e bave = %7.3e" ] ;
+          //Print[ {$sarea~{i}}, Format "DEBUG %7.3e" ] ;
+        EndFor
+      EndIf
     }
- }
+  }
 }
 
 PostProcessing {
   { Name PostProcessings; NameOfFormulation MagnetoDynamics;
     PostQuantity {
       { Name az;
-		Value { Local { [ CompZ[{ac}] - (Flag_Thin!=0)*(CompZ[{aw}]-Correction_A[])];
-			In  Dom_Hcurl_a; Jacobian Vol; }}}
+        Value { Local { [ CompZ[{ac}] - (Flag_Thin!=0)*(CompZ[{aw}]-Correction_A[])];
+            In  Dom_Hcurl_a; Jacobian Vol; }}}
       { Name acz;
-		Value { Local { [ CompZ[ {ac}] ];
-			In  Dom_Hcurl_a; Jacobian Vol; }}}
+        Value { Local { [ CompZ[ {ac}] ];
+            In  Dom_Hcurl_a; Jacobian Vol; }}}
       { Name awz;
-		Value { Local { [ CompZ[ {aw}] ];
-			In Dom_Hthin_a; Jacobian Vol; }}}
+        Value { Local { [ CompZ[ {aw}] ];
+            In Dom_Hthin_a; Jacobian Vol; }}}
       { Name acwz;
-		Value { Local { [ CompZ[ {ac} - {aw} ] ];
-			In Dom_Hcurl_a; Jacobian Vol; }}}
+        Value { Local { [ CompZ[ {ac} - {aw} ] ];
+            In Dom_Hcurl_a; Jacobian Vol; }}}
       { Name exact;
-		Value { Local { [ Exact_A[] ] ;
-			In Dom_Hcurl_a; Jacobian Vol; }}}
+        Value { Local { [ Exact_A[] ] ;
+            In Dom_Hcurl_a; Jacobian Vol; }}}
       { Name corr;
-		Value { Local { [ Correction_A[] ] ;
-			In Dom_Hcurl_a; Jacobian Vol; }}}
+        Value { Local { [ Correction_A[] ] ;
+            In Dom_Hcurl_a; Jacobian Vol; }}}
       { Name bcw;
-	   	Value { Local { [ {d ac} - {d aw} ];
-	   		In Dom_Hcurl_a; Jacobian Vol; }}}
+        Value { Local { [ {d ac} - {d aw} ];
+            In Dom_Hcurl_a; Jacobian Vol; }}}
 
       { Name j;
-		Value { Local { [ -sigma[] * ( Dt[{a}] + {dv} ) ]; 
-			In Vol_C; Jacobian Vol; }}}
-
-	  { Name un;
-		Value{ Local { [ 1 ] ;
-			In Vol_nu; Integration I1 ; Jacobian Vol; }}}  
-	  { Name surf;
-		Value{ Integral{ [ 1 ] ;
-			In Vol_nu; Integration I1 ; Jacobian Vol; }}}
-
-
-	  If( Flag_Thin == 0 )
-		{ Name flux;
-		  Value { Integral { [ CompZ[ {a} ] / A_c  ]; 
-			  In Vol_C; Integration I1 ; Jacobian Vol; }}}
+        Value { Local { [ -sigma[] * ( Dt[{a}] + {dv} ) ]; 
+            In Vol_C; Jacobian Vol; }}}
+
+      { Name un;
+        Value{ Local { [ 1 ] ;
+            In Vol_nu; Integration I1 ; Jacobian Vol; }}}  
+      { Name surf;
+        Value{ Integral{ [ 1 ] ;
+            In Vol_nu; Integration I1 ; Jacobian Vol; }}}
+
+      If( Flag_Thin == 0 )
+        { Name flux;
+          Value { Integral { [ CompZ[ {a} ] / A_c  ]; 
+              In Vol_C; Integration I1 ; Jacobian Vol; }}}
         { Name JouleLosses;
-		  Value { Integral { [ sigma[]*SquNorm[ Dt[{a}] + {dv} ] ] ;
-			  In Vol_C ; Integration I1 ; Jacobian Vol ; }}}
-		{ Name U;
+          Value { Integral { [ sigma[]*SquNorm[ Dt[{a}] + {dv} ] ] ;
+              In Vol_C ; Integration I1 ; Jacobian Vol ; }}}
+        { Name U;
           Value { Term { [  Cart2Pol [ {U} ] ]; In Vol_C; }}}
-		{ Name L;
-		  Value { Term { [  Im [ {U}/{I}/omega ] ]; In Vol_C; }}}
+        { Name L;
+          Value { Term { [  Im [ {U}/{I}/omega ] ]; In Vol_C; }}}
 
-	  Else
-		{ Name flux; // naive flux
-		  Value { Integral { [ CompZ[ {ac} ] ]; 
-			  In lVol_C; Integration I1 ; Jacobian Vol; }}}
+        Else
+        { Name flux; // naive flux
+          Value { Integral { [ CompZ[ {ac} ] ]; 
+              In lVol_C; Integration I1 ; Jacobian Vol; }}}
 
-	    { Name JouleLosses; 
-		  Value { Term { [ Analytic_R[]*{Is}*{Is} ] ; In lVol_C; }}}
-		{ Name U;
+        { Name JouleLosses; 
+          Value { Term { [ Analytic_R[]*{Is}*{Is} ] ; In lVol_C; }}}
+        { Name U;
           Value { Term { [ Cart2Pol [ {Us} ] ] ; In lVol_C; }}}
         { Name L; // only valid with one current 
           Value { Term { [ Im [ {Us}/{Is}/omega ] ] ; In lVol_C; }}}
-
-	    For i In {1:NumWires}
-		  { Name area~{i};
-		    Value{ Integral{ [ 1 ] ;
-				In SWIRE~{i}; Integration I1 ; Jacobian Vol; }}}
-		  { Name intbx~{i}; // integral of b, 1st step for computing b average 
-			 Value { Integral { [ Norm [ CompX[ {d ac} - {d aw} ] ] ] ; 
-				 In SWIRE~{i}; Integration I1 ; Jacobian Vol; }}}
-		  { Name intby~{i}; // integral of b, 1st step for computing b average 
-			 Value { Integral { [ Norm [ CompY[ {d ac} - {d aw} ] ] ] ; 
-				 In SWIRE~{i}; Integration I1 ; Jacobian Vol; }}}
-	    EndFor
-		
-	  EndIf
+        
+        For i In {1:NumWires}
+          { Name area~{i};
+            Value{ Integral{ [ 1 ] ;
+                In SWIRE~{i}; Integration I1 ; Jacobian Vol; }}}
+          { Name intbx~{i}; // integral of b, 1st step for computing b average 
+            Value { Integral { [ Norm [ CompX[ {d ac} - {d aw} ] ] ] ; 
+                In SWIRE~{i}; Integration I1 ; Jacobian Vol; }}}
+          { Name intby~{i}; // integral of b, 1st step for computing b average 
+            Value { Integral { [ Norm [ CompY[ {d ac} - {d aw} ] ] ] ; 
+                In SWIRE~{i}; Integration I1 ; Jacobian Vol; }}}
+        EndFor
+      EndIf
     }
   }
 }
@@ -515,22 +476,22 @@ PostProcessing {
 PostOperation mapaz UsingPost PostProcessings {
   Print [az, OnElementsOf Vol_nu, File "az.pos"];
   Echo[ StrCat["l=PostProcessing.NbViews-1;",
-			   "View[l].IntervalsType = 3;",
-			   "View[l].NbIso = 30;",
-			   "View[l].RaiseZ = 15000;",
-			   "View[l].NormalRaise = 0;"],
-		File "tmp.geo", LastTimeStepOnly];
+      "View[l].IntervalsType = 3;",
+      "View[l].NbIso = 30;",
+      "View[l].RaiseZ = 15000;",
+      "View[l].NormalRaise = 0;"],
+    File "tmp.geo", LastTimeStepOnly];
  
   If(Flag_Thin != 0)
     Print [bcw, OnElementsOf Vol_Sleeve, File "b.pos"];
 
     Print [acwz, OnElementsOf Vol_nu, File "acz.pos"];
-	Echo[ StrCat["l=PostProcessing.NbViews-1;",
-				 "View[l].IntervalsType = 3;",
-				 "View[l].NbIso = 30;",
-				 "View[l].RaiseZ = 15000;",
-				 "View[l].NormalRaise = 0;"],
-		   File "tmp.geo", LastTimeStepOnly] ;
+    Echo[ StrCat["l=PostProcessing.NbViews-1;",
+        "View[l].IntervalsType = 3;",
+        "View[l].NbIso = 30;",
+        "View[l].RaiseZ = 15000;",
+        "View[l].NormalRaise = 0;"],
+      File "tmp.geo", LastTimeStepOnly] ;
   EndIf
 }
 
@@ -542,53 +503,53 @@ PostOperation cutaz UsingPost PostProcessings {
   // cut of az displayed in GUI
   If( Flag_Thin == 0 )
     Print [ az, OnLine { {-Xmax,0,0} {Xmax,0,0} }{NbPoints},
-			Format Gmsh, File "Cut_az.pos" ];
+      Format Gmsh, File "Cut_az.pos" ];
     Echo[ StrCat["l=PostProcessing.NbViews-1;",
-				 "View[l].Name = 'cut az full';",
-				 "View[l].Axes = 3;",
-				 "View[l].LineWidth = 3;",
-				 "View[l].Type = 2;"],
-	File "tmp.geo", LastTimeStepOnly];
+        "View[l].Name = 'cut az full';",
+        "View[l].Axes = 3;",
+        "View[l].LineWidth = 3;",
+        "View[l].Type = 2;"],
+      File "tmp.geo", LastTimeStepOnly];
 
-  Else
+    Else
 
-	Print [ acz, OnLine { {-Xmax,0,0} {Xmax,0,0} }{NbPoints},
-			Format Gmsh, File "Cut_az.pos" ];
+    Print [ acz, OnLine { {-Xmax,0,0} {Xmax,0,0} }{NbPoints},
+      Format Gmsh, File "Cut_az.pos" ];
     Echo[ StrCat["l=PostProcessing.NbViews-1;",
-				 "View[l].Name = 'cut az not corrected';",
-				 "View[l].Axes = 3;",
-				 "View[l].LineWidth = 3;",
-				 "View[l].Type = 2;"],
-		  File "tmp.geo", LastTimeStepOnly];
+        "View[l].Name = 'cut az not corrected';",
+        "View[l].Axes = 3;",
+        "View[l].LineWidth = 3;",
+        "View[l].Type = 2;"],
+      File "tmp.geo", LastTimeStepOnly];
 
     Print [ az, OnLine { {-Xmax,0,0} {Xmax,0,0} }{NbPoints},
-			Format Gmsh, File "Cut_az.pos" ];
+      Format Gmsh, File "Cut_az.pos" ];
     Echo[ StrCat["l=PostProcessing.NbViews-1;",
-				 "View[l].Name = 'cut az corrected';",
-				 "View[l].Axes = 3;",
-				 "View[l].LineWidth = 3;",
-				 "View[l].Type = 2;"],
-	File "tmp.geo", LastTimeStepOnly];
+        "View[l].Name = 'cut az corrected';",
+        "View[l].Axes = 3;",
+        "View[l].LineWidth = 3;",
+        "View[l].Type = 2;"],
+      File "tmp.geo", LastTimeStepOnly];
   EndIf
 
   // cuts in txt format for Gnuplot
   Print [ exact, OnLine { {0,0,0} {Xmax,0,0} } {NbPoints},
-	  Format SimpleTable, File "Cut_analytic.txt" ];
+    Format SimpleTable, File "Cut_analytic.txt" ];
 
   If( Flag_Thin == 0 )
     Print [ az, OnLine { {0,0,0} {Xmax,0,0} } {NbPoints},
-	    Format SimpleTable, File "Cut_fullfem.txt" ];
+      Format SimpleTable, File "Cut_fullfem.txt" ];
   Else
     Print [ az, OnLine { {0,0,0} {Xmax,0,0} } {NbPoints},
-	    Format SimpleTable, File "Cut_thinfem.txt" ];
+      Format SimpleTable, File "Cut_thinfem.txt" ];
     Print [ acz, OnLine { {0,0,0} {Xmax,0,0} } {NbPoints},
-	    Format SimpleTable, File "Cut_ac.txt" ];
+      Format SimpleTable, File "Cut_ac.txt" ];
     Print [ awz, OnLine { {0,0,0} {Xmax,0,0} } {NbPoints},
-	    Format SimpleTable, File "Cut_aw.txt" ];
+      Format SimpleTable, File "Cut_aw.txt" ];
     Print [ acwz, OnLine { {0,0,0} {Xmax,0,0} } {NbPoints},
-	    Format SimpleTable, File "Cut_acw.txt" ];
+      Format SimpleTable, File "Cut_acw.txt" ];
     Print [ corr, OnLine { {0,0,0} {Xmax,0,0} } {NbPoints},
-	    Format SimpleTable, File "Cut_correction.txt" ];
+      Format SimpleTable, File "Cut_correction.txt" ];
   EndIf
 }
 
@@ -597,111 +558,53 @@ PostOperation integaz UsingPost PostProcessings {
   For i In {1:NumWires}
 
     If( Flag_Thin == 0 )
-	  Print[ U, OnRegion Vol_C, File > "U_full.dat", Format Table,
-			 SendToServer Sprintf["}Voltage/Wire %g/1Full",i]];
+      Print[ U, OnRegion Vol_C, File > "U_full.dat", Format Table,
+        SendToServer Sprintf["}Voltage/Wire %g/1Full",i]];
       Print[ L, OnRegion Vol_C, File > "L_full.dat", Format Table,
         SendToServer Sprintf["}Inductance/Wire %g/1Full",i]];
-	  Print[ JouleLosses[ VWIRE~{i}], OnGlobal,
-			 Format Table, File > "joule_full.dat", 
-	 		 SendToServer Sprintf["}Joule Losses/Wire %g/1Full",i]];
+      Print[ JouleLosses[ VWIRE~{i}], OnGlobal,
+        Format Table, File > "joule_full.dat", 
+        SendToServer Sprintf["}Joule Losses/Wire %g/1Full",i]];
     EndIf
-	If( Flag_Thin == 1 )
-	  Print[ U, OnRegion lVol_C, File > "U_raw.dat", Format Table,
+    
+    If( Flag_Thin == 1 )
+      Print[ U, OnRegion lVol_C, File > "U_raw.dat", Format Table,
         SendToServer Sprintf["}Voltage/Wire %g/2Raw",i]];
       Print[ L, OnRegion lVol_C, File > "L_raw.dat", Format Table,
         SendToServer Sprintf["}Inductance/Wire %g/2Raw",i]];
-	  Print[ JouleLosses, OnRegion LWIRE~{i},
+      Print[ JouleLosses, OnRegion LWIRE~{i},
         Format Table, File > "joule_raw.dat", 
         SendToServer Sprintf["}Joule Losses/Wire %g/2Raw",i]];
     EndIf
-	If( Flag_Thin == 2 )
-	  Print[ U, OnRegion lVol_C, File > "U_reg.dat", Format Table,
-			 SendToServer Sprintf["}Voltage/Wire %g/3Reg",i]];
+    
+    If( Flag_Thin == 2 )
+      Print[ U, OnRegion lVol_C, File > "U_reg.dat", Format Table,
+        SendToServer Sprintf["}Voltage/Wire %g/3Reg",i]];
       Print[ L, OnRegion lVol_C, File > "L_reg.dat", Format Table,
         SendToServer Sprintf["}Inductance/Wire %g/3Reg",i]];
-  	  Print[ JouleLosses, OnRegion LWIRE~{i},
-			 Format Table, File > "joule_reg.dat", 
-			 SendToServer Sprintf["}Joule Losses/Wire %g/3Reg",i]];
+      Print[ JouleLosses, OnRegion LWIRE~{i},
+        Format Table, File > "joule_reg.dat", 
+        SendToServer Sprintf["}Joule Losses/Wire %g/3Reg",i]];
     EndIf
-	If( Flag_Thin == 3 )
-	  Print[ U, OnRegion lVol_C, File > "U_naive.dat", Format Table,
-			 SendToServer Sprintf["}Voltage/Wire %g/4Naive",i]];
+
+    If( Flag_Thin == 3 )
+      Print[ U, OnRegion lVol_C, File > "U_naive.dat", Format Table,
+        SendToServer Sprintf["}Voltage/Wire %g/4Naive",i]];
       Print[ L, OnRegion lVol_C, File > "L_naive.dat", Format Table,
         SendToServer Sprintf["}Inductance/Wire %g/4naive",i]];
-  	  Print[ JouleLosses, OnRegion LWIRE~{i},
-			 Format Table, File > "joule_naive.dat", 
-			 SendToServer Sprintf["}Joule Losses/Wire %g/4Naiv",i]];
+      Print[ JouleLosses, OnRegion LWIRE~{i},
+        Format Table, File > "joule_naive.dat", 
+        SendToServer Sprintf["}Joule Losses/Wire %g/4Naiv",i]];
     EndIf
 
-	If( Flag_Thin != 0 )
-      Print[ area~{i}[SWIRE~{i}], OnGlobal,
-			 Format Table, File >"zzz",
-			 StoreInVariable $sarea~{i}];
+    If( Flag_Thin != 0 )
+      Print[ surf[SWIRE~{i}], OnGlobal,
+        Format Table, File >"zzz", StoreInVariable $sarea~{i}];
       Print[ intby~{i}[SWIRE~{i}], OnGlobal,
-	  		 Format Table, File > "zzz",
-	  		 StoreInVariable $intby~{i}];	  
-	EndIf
+        Format Table, File > "zzz", StoreInVariable $intby~{i}];
+    EndIf
 
   EndFor
 }
 
 
-
-
-	  /*{ Name flux_corr; // corrected flux 
-		  Value { 
-			Integral { [ CompZ[ {ac} - {aw} ] ]; 
-			  In lVol_C; Integration I1 ; Jacobian Vol; }
-			Integral { [ Analytic_L[]*{Is} ]; 
-			  In lVol_C; Integration I1 ; Jacobian Vol; }
-		  }
-		  }*/
-
-    /* Print [awz, OnElementsOf Vol_nu, File "awz.pos"];
-    Echo[ StrCat["l=PostProcessing.NbViews-1;",
-				 "View[l].IntervalsType = 3;",
-				 "View[l].NbIso = 30;",
-				 "View[l].RaiseZ = 15000;",
-				 "View[l].NormalRaise = 0;"],
-		  File "tmp.geo", LastTimeStepOnly] ;
-
-    Print [acz, OnElementsOf Vol_nu, File "acz.pos"];
-    Echo[ StrCat["l=PostProcessing.NbViews-1;",
-				 "View[l].IntervalsType = 3;",
-				 "View[l].NbIso = 30;",
-				 "View[l].RaiseZ = 15000;",
-				 "View[l].NormalRaise = 0;"],
-			File "tmp.geo", LastTimeStepOnly] ; 
-
-
-
-      Print[ flux[ VWIRE~{i} ], OnGlobal,
-			 Format Table, File > "flux_full.dat", 
-			 SendToServer Sprintf["}Flux/Wire %g/1Full",i]];
-      Print[ flux[ LWIRE~{i} ], OnGlobal,
-			 Format Table, File > "flux_raw.dat", 
-			 SendToServer Sprintf["}Flux/Wire %g/2Raw",i]];
-      Print[ flux[ LWIRE~{i} ], OnGlobal,
-			 Format Table, File > "flux_reg.dat", 
-			 SendToServer Sprintf["}Flux/Wire %g/4Reg",i]];
-*/
-  /*
-  X = rw/skin_depth;
-  pI[] = (X/2)*Re[(1+i[])*(J_0[(1+i[])*X]/J_1[(1+i[])*X])];
-  Resist[] = Re[pI[]*R_dc];
- // radial analytical solution with a zero flux condition imposed at r=rs
-    Analytic_a~{i}[] = mu0/(2*Pi) * WI~{i} * 
-	    ( (R~{i}[]>WR~{i}) ? Log[rs/R~{i}[]] :
-		Log[rs/WR~{i}] + (J0[tau[]*R~{i}[]/rw] - J0[tau[]]) / J1[tau[]] / tau[] ); 
-  
-    Analytic_L~{i}[] = mu0 / (2*Pi) * 
-		   ( 2/tau[]/tau[]- J0[tau[]] / J1[tau[]] / tau[] - Log[rs/WR~{i}] ) ; 
-
-    // Idem in static case (thus particular case of the above)
-    AnalyticStatic_a~{i}[] = mu0 * WI~{i} / (2*Pi) * 
-       ( (R~{i}[]>WR~{i}) ? Log[rs/R~{i}[]] : Log[rs/WR~{i}] + (1-(R~{i}[]/WR~{i})^2)/2. );
-    Analytic_e~{i}[] = mu0 * WI~{i} / (2*Pi) * i[] * omega
-						    * J0[tau[]*R~{i}[]/WR~{i}] / J1[tau[]] / tau[] ; 
-    Analytic_flux~{i}[] = mu0 * WI~{i} / (2*Pi) * 
-      ( i[] * (skin_depth/WR~{i})^2 + Log[rs/WR~{i}] - J0[tau[]] / J1[tau[]] / tau[] ) ; 
-  */
diff --git a/Conductors3D/w3d.geo b/Conductors3D/w3d.geo
index cdc2a4c36b18c3234fda10407ffbe39caaab55dd..521752c6a744034f83eecd986a19f23e142539bd 100644
--- a/Conductors3D/w3d.geo
+++ b/Conductors3D/w3d.geo
@@ -116,7 +116,10 @@ ElseIf( Flag_RegSleeve ) // regular sleeve
   s1 = news; Plane Surface(s1) = {ll1,-llWires[]};
   e[] = Extrude {0,0,Lz} { Surface{ s1 } ; /*Layers{1}; Recombine;*/ } ;
   Physical Volume ("AIR", 1) = { e[1], Wires[] };
-  Physical Volume ("BLA", 5) = { e[1] };
+  // auxiliary Physical Volume to test the tree.
+  // Note that two Physical Volumes defined on a same region
+  // yield duplicate finite element in GetDP. 
+  //Physical Volume ("BLA", 5) = { e[1] }; this 
 
 Else
 
@@ -148,4 +151,5 @@ If( !Flag_Thin ) // Wires with finite radius
     Electrodes[] += {20+i, 30+i};
   EndFor
 EndIf  
-Physical Line("LINTREE", 4) = { Boundary { Physical Surface { Electrodes[] }; } };
+
+Physical Line("LINTREE", 4) = { CombinedBoundary { Physical Surface { Electrodes[] }; } };
diff --git a/Conductors3D/w3d.pro b/Conductors3D/w3d.pro
index 69d45783ca8360459587665e02adc9f181eb4b4c..95e4db75a56e269a6c270bfe041573acfe4e025b 100644
--- a/Conductors3D/w3d.pro
+++ b/Conductors3D/w3d.pro
@@ -33,8 +33,7 @@ Group{
   AIR = Region[ 1 ];
   INF = Region[ 2 ];
   SKIN = Region[ 3 ];
-  LINTREE = Region[ 4 ]; // not used
-  BLA = Region[ 5 ]; 
+  LINTREE = Region[ 4 ];
 
   VWIRES = Region[ {} ];
   LWIRES = Region[ {} ];
@@ -62,7 +61,6 @@ Group{
     CATHODES += Region[ CATHODE~{i} ];
   EndFor
 
-
   // Abstract regions
 
   If( !Flag_Thin )
@@ -76,7 +74,7 @@ Group{
   Vol_nu = Region[ {Vol_C, Vol_CC} ];
   Sur_Dirichlet_a = Region[ INF ];
   Electrodes = Region[ {ANODES, CATHODES} ];   
-  
+      
   Dom_Hcurl_a = Region[ Vol_nu ];
   Dom_Hgrad_u = Region[ Vol_C ];
   Dom_Hregion_i = Region[ Vol_C ];
@@ -84,7 +82,6 @@ Group{
   // additional Groups for the semi_analytic approach
   Dom_Hthin_a = ElementsOf[ Vol_nu, OnOneSideOf LWIRES ];
   Vol_Tree = ElementsOf[ Vol_nu, DisjointOf LWIRES ];
-  //Vol_Tree = Region[ { BLA } ];
   Sur_Tree = Region[ { Sur_Dirichlet_a /*, SKIN*/ } ];
   If( !Flag_SemiAnalytic )
     Lin_Tree = Region[ {} ];
@@ -92,7 +89,6 @@ Group{
     Lin_Tree = Region[ { LWIRES } ];
   EndIf
 }
-  
 
 Function{
   omega = 2*Pi*Freq;
@@ -267,7 +263,7 @@ FunctionSpace {
       { NameOfCoef ac;
 		EntityType Auto; NameOfConstraint Hcurl_a_3D; }
 	  { NameOfCoef ac; EntityType EdgesOfTreeIn ; EntitySubType StartingOn ;
-		NameOfConstraint GaugeCondition_a ; }
+	    NameOfConstraint GaugeCondition_a ; }
     }
   }
   
@@ -576,16 +572,13 @@ PostOperation map UsingPost MagnetoDynamics {
     Print[ J, OnElementsOf Vol_C, File "j.pos"];
   Else
     Print[ bs, OnElementsOf Vol_nu, File "bs.pos"];
-    PrintGroup[ EdgesOf[ Vol_nu, ConnectedTo LWIRES ],
-      In Vol_nu, File "group.pos"];
-    PrintGroup[ _CO_Entity_71, In Vol_nu, File "Tree.pos"];
+    // PrintGroup[ EdgesOf[ Vol_nu, ConnectedTo LWIRES ],
+    //   In Vol_nu, File "group.pos"];
+    //PrintGroup[ _CO_Entity_71, In Vol_nu, File "Tree.pos"];
   EndIf
-  
-  // PrintGroup[ EdgesOfTreeIn[ { Vol_Tree }, StartingOn { Sur_Tree, Lin_Tree } ],
-  //   In Vol_nu, File "Tree.pos"];
-
-  
-
+  //PrintGroup[ _CO_Entity_30, In Vol_nu, File "Tree.pos"];
+    PrintGroup[ EdgesOfTreeIn[ Vol_Tree , StartingOn { Sur_Tree, Lin_Tree } ],
+      In Vol_nu, File "Tree.pos"];
 }
 
 PostOperation integaz UsingPost MagnetoDynamics {
diff --git a/Conductors3D/w3d_common.pro b/Conductors3D/w3d_common.pro
index 8bb6be1e77fc806bf42ed068eab59d570a9a0c7b..3f3af2a4c871556ab3031503cd2bfabd6ee00c1c 100644
--- a/Conductors3D/w3d_common.pro
+++ b/Conductors3D/w3d_common.pro
@@ -22,8 +22,8 @@ DefineConstant[
   
   Flag_Stranded = {0, Name "Parameters/04Stranded conductors", 
                    Choices {0,1}, Visible !Flag_SemiAnalytic || !Flag_Thin}
-  Flag_Dual = {0, Name "Parameters/05Dual approach", 
-    Choices {0,1}, Visible Flag_SemiAnalytic && Flag_Thin}
+  // Flag_Dual = {0, Name "Parameters/05Dual approach", 
+  //   Choices {0,1}, Visible Flag_SemiAnalytic && Flag_Thin}
   Flag_U = {0, Name "Parameters/06Impose U", Choices {0,1}, Visible 1}
   
   WireRadius   = {1, Name "Parameters/10Wire radius [mm]"}