[Getdp] Magnetodynamics a-v formulation 3D ?

Thomas Jung Thomas.Jung at iisb.fraunhofer.de
Tue Apr 4 16:07:57 CEST 2006


On Friday 31 March 2006 15:45, Gilles Grégoire wrote:
> Le Vendredi 31 Mars 2006 15:20, Thomas Jung a écrit :
> > Hello everybody,
>
> Hello Thomas
>
> > I have spent some time trying to make a 3D example based on the example
> > in the documentation (a-v-formulation, 2D, magnetodynamics) - but am
> > apparently too stupid, cant get it running.
> >
> > So I thought maybe someone has done that already and could post an
> > example ?
>
> I haven't tryed this yet, but this could be a small mistake in the way you
> described your problem.
> Maybe you could post the files (*.geo and *.pro) describing your problem?
>
> > I am aware of the 3D magnetodynamics example in the wiki, but would like
> > to avoid the H-phi-formulation ....
> >
> > Thank you very much in advance !
> >
> > Thomas


Hello Gilles,

Here are the files (the .geo is a little lengthy, so I attach it).

I am trying to have a conductor ( "Coil"), with two ends where I want to 
prescribe electrical potential, get the current computed, and the induced 
current in the susceptor, which for the moment, is still not present 
(conductivity=0), because I still fail to compute the vector potential a.

I tried lots o f different combinations of solvers, preconditioners and other 
solver parameters, but get typically something like this:


Loading   : Problem definition 'InductionTestBox.pro'
Loading   : Problem definition 'Jacobian_Lib.pro'
Loading   : Problem definition 'Integration_Lib.pro'
Loading   : Problem definition 'MagDyn_av_3D.pro'
Info      : Selected Resolution 'MagDyn_vf'
Loading   : Geometric data 'InductionTestBox.msh'
Info      : System 'Sys_MagDyn' : Complex, Frequency = 10 Hz
P r e - P r o c e s s i n g . . .
Operation : Treatment Formulation 'MagDyn_avf'
Info      :   Generate ExtendedGroup '_CO_Entity_14' (EdgesOf)
Info      :   Generate ExtendedGroup '_CO_Entity_11' (NodesOf)
Info      :   Generate ExtendedGroup '_CO_Entity_12' (NodesOf)
Resources : cpu 5.812363 s / mem 0 kb
E n d   P r e - P r o c e s s i n g
P r o c e s s i n g . . .
Operation : Generate[Sys_MagDyn]
Solver    : Loading parameter file 'solver.par'
Info      : Setting System {A,b} to zero
Resources : cpu 16.5000 s / mem 0 kb
Operation : Solve[Sys_MagDyn]
Solver    : Scaling system of equations
Solver    : RCMK algebraic renumbering
Solver    : N: 264674, NZ: 4949840, BW max/avg: 110/18, SW max: 260437
Resources : cpu 21.325332 s / mem 0 kb
Solver    : ILUTP (Float, fill-in = 20)
Solver    : N: 264674, NZ: 10317084, BW max/avg: 40/38, SW max: 264252
Resources : cpu 1240.161505 s / mem 0 kb
Solver    : Direct version of Quasi Generalize Minimum RESidual (DQGMRES)
    1  5.8768961e+02  1.0000000e+00
Warning   : Iterative solver is facing a break-down
    2  5.8768961e+02  1.0000000e+00
Solver    : 2 Iterations / Residual: 587.69
Resources : cpu 1240.789544 s / mem 0 kb
Operation : SaveSolution[Sys_MagDyn]
Resources : cpu 1241.637597 s / mem 0 kb
E n d   P r o c e s s i n g
E n d

The problematic case is the resolution "MagDyn_vf"

Millions of thanks to anybody who should have a look at this !



Problem description (InductionTestBox.pro):

=========================================================================

/* InductionTestBox.pro */

Group {
   
  Air = Region[1];
  Coil = Region[2];
  Susceptor = Region[3];
  Ground1 = Region[12];
  Ground2 = Region[16];
  OuterBoundary = Region[{101,102,103,104,5,6}];

  /* conducting regions = DomainC_Mag */
  DomainC_Mag = Region[Coil];

  /* non-conducting regions = DomainCC_Mag */
  DomainCC_Mag = Region[{Air,Susceptor}];

  /* entire domain */
  Domain_Mag = Region[{DomainC_Mag, DomainCC_Mag}];

  /* domain with source current */
  DomainS_Mag = Region[Coil];
}

Function {

  epsr[Air] = 1.;
  epsr[Coil] = 1.;
  epsr[Susceptor] = 1.;

}

Function {

  sigma[Coil]=1.e6;
  sigma[Air]=0.;
  sigma[Susceptor]=0.;

}

Function {
    mu0 = 4.e-7 * Pi;

    nu [ Domain_Mag]  = 1. / mu0;

}

Function {

    js[Coil] = Vector[1000 ,0,0];

}

Constraint {

  { Name ElectricScalarPotential; Type Assign;
    Case {
      { Region Ground1; Value 0.; }
      { Region Ground2; Value 1.; }
    }
  }

  { Name MagneticVectorPotential_3D; Type Assign;
    Case {
      { Region OuterBoundary ; Value 0.; }
    }
  }

}

Include "Jacobian_Lib.pro"
Include "Integration_Lib.pro"
Include "MagDyn_av_3D.pro"

PostOperation {
  { Name MapE; NameOfPostProcessing EleSta_v;
     Operation {
       Print [ v, OnElementsOf DomainC_Mag, File "Coil_v.pos" ];
       Print [ e, OnElementsOf DomainC_Mag, File "Coil_e.pos" ];
       Print [ j, OnElementsOf DomainC_Mag, File "Coil_j.pos" ];
     }
  }

}

PostOperation {
  { Name MapA; NameOfPostProcessing MagSta;
     Operation {
       Print [ a, OnElementsOf Domain_Mag, File "MagSta_a.pos" ];
       Print [ h, OnElementsOf Domain_Mag, File "MagSta_h.pos" ];
     }
  }

}

PostOperation {
  { Name MapDyn; NameOfPostProcessing MagDyn_avf;
     Operation {
       Print [ a, OnElementsOf Domain_Mag, File "MagDyn_a.pos" ];
       Print [ h, OnElementsOf Domain_Mag, File "MagDyn_h.pos" ];
       Print [ j, OnElementsOf DomainC_Mag, File "MagDyn_j.pos" ];
     }
  }

}

=============================================================================

Function spaces, formulations, etc: 


/* -------------------------------------------------------------------
   File "MagDyn_av_3D.pro"
*/

eps0 = 8.854187818e-12;

Freq1=10;

Group {
  DefineGroup[ Domain_Mag, DomainS_Mag, DomainCC_Mag, DomainC_Mag ];
}

Function {
  DefineFunction[ epsr ];
}

Function {
  DefineFunction[ sigma ];
}

Function {
    DefineFunction[ nu ];
}

/* for electric potential */

FunctionSpace {
  { Name v_Ele; Type Form0;
    BasisFunction {
      // v = v  s   ,  for all nodes
      //      n  n
      { Name sn; NameOfCoef vn; Function BF_Node;
        Support DomainC_Mag; Entity NodesOf[ All ]; }
    }
    Constraint {
      { NameOfCoef vn; EntityType NodesOf; 
        NameOfConstraint ElectricScalarPotential; }
    }
  }

  // Magnetic vector potential a (b = curl a)
  { Name Hcurl_a_Mag_3D; Type Form1;
    BasisFunction {
	{ Name se; NameOfCoef ae; Function BF_Edge;
        Support Domain_Mag; Entity EdgesOf[ All ]; }
    }
    Constraint {
	{ NameOfCoef ae; EntityType EdgesOf;
        NameOfConstraint MagneticVectorPotential_3D; }
    }
  }
}

Formulation {
  { Name Electrostatics_v; Type FemEquation;
    Quantity {
      { Name v; Type Local; NameOfSpace v_Ele; }
    }
    Equation {
      Galerkin { [ epsr[] * Dof{d v} , {d v} ]; In DomainC_Mag; 
                 Jacobian Vol; Integration CurlCurl; }
    }
  }

  { Name Magnetostatics_a_3D; Type FemEquation;
    Quantity {
      { Name a ; Type Local; NameOfSpace Hcurl_a_Mag_3D; }
      { Name v ; Type Local; NameOfSpace v_Ele; }
      { Name j ; Type Local; NameOfSpace v_Ele; }
    }
    Equation {
      Galerkin { [ nu[] * Dof{d a} , {d a} ]; In Domain_Mag;
                 Jacobian Vol; Integration CurlCurl; }
      Galerkin { [ - js[] , {a} ]; In DomainS_Mag;
                 Jacobian Vol; Integration CurlCurl; }
      //Galerkin { [ -sigma[]*Dof{j} , {a} ]; In DomainS_Mag;
      //          Jacobian Vol; Integration CurlCurl; }
    }
  }

  { Name MagDyn_avf; Type FemEquation ;
    Quantity {
      { Name a ; Type Local  ; NameOfSpace Hcurl_a_Mag_3D ; }
      { Name v ; Type Local  ; NameOfSpace v_Ele; }
      //{ Name I ; Type Global ; NameOfSpace Hgrad_vf_Mag [I] ; }
      //{ Name V ; Type Global ; NameOfSpace Hgrad_vf_Mag [V] ; }
      // { Name vf; Type Local  ; NameOfSpace Hgrad_vf_Mag [vf] ; }
    }
    Equation {
      Galerkin { [ nu[] * Dof{d a} , {d a} ]; In Domain_Mag;
                 Jacobian Vol; Integration CurlCurl; }

      Galerkin { DtDof [ sigma[] * Dof{a} , {a} ]; In DomainC_Mag;
                 Jacobian Vol; Integration CurlCurl; }
      Galerkin { [ sigma[] * Dof{d v} , {a} ]; In DomainC_Mag;
                 Jacobian Vol; Integration CurlCurl; }
      
      Galerkin { DtDof [ sigma[] * Dof{a} , {d v} ]; In DomainC_Mag;
                 Jacobian Vol; Integration CurlCurl; }
      Galerkin { [ sigma[] * Dof{d v} , {d v} ]; In DomainC_Mag;
                 Jacobian Vol; Integration CurlCurl; }
      // GlobalTerm { [ Dof{I} , {V} ]; In Electrode_Mag; }
    }
  }
}

Resolution {
  { Name current;
    System {
      { Name Sys_Ele; NameOfFormulation Electrostatics_v; }
    }
    Operation { 
      Generate[Sys_Ele]; Solve[Sys_Ele]; SaveSolution[Sys_Ele];
    }
  }

  { Name MagSta_a_3D;
    System {
      { Name Sys_Mag; NameOfFormulation Magnetostatics_a_3D; }
    }
    Operation {
      Generate[Sys_Mag]; Solve[Sys_Mag]; SaveSolution[Sys_Mag];
    }
  }

  { Name MagDyn_vf;
    System {
      { Name Sys_MagDyn; NameOfFormulation MagDyn_avf;
        Type ComplexValue; Frequency {Freq1}; }
    }
    Operation { 
      Generate[Sys_MagDyn]; Solve[Sys_MagDyn]; SaveSolution[Sys_MagDyn];
    }
  }

  { Name both ;
    System {
      { Name A1 ; NameOfFormulation Electrostatics_v ; }
      { Name A2 ; NameOfFormulation Magnetostatics_a_3D; }
    }
    Operation { 
	Generate[A1] ; Solve[A1] ; SaveSolution[A1] ;
      Generate[A2] ; Solve[A2] ; SaveSolution[A2] ;
    }
  }
}


PostProcessing {
  { Name EleSta_v; NameOfFormulation Electrostatics_v;
    Quantity {
      { Name v; 
        Value { 
          Local { [ {v} ]; In DomainC_Mag; Jacobian Vol; } 
        }
      }
      { Name e; 
        Value { 
          Local { [ -{d v} ]; In DomainC_Mag; Jacobian Vol; }
        }
      }
      { Name d; 
        Value { 
          Local { [ -eps0*epsr[] * {d v} ]; In DomainC_Mag; 
                                             Jacobian Vol; } 
        } 
      }
      { Name j; 
        Value { 
          Local { [ sigma[] * {d v} ]; In DomainC_Mag; 
                                             Jacobian Vol; } 
        } 
      }


    }
  }
}

PostProcessing {
  { Name MagSta; NameOfFormulation Magnetostatics_a_3D;
    Quantity {
	{ Name h ; Value { Term { [ nu[] * {d a} ]     ; In Domain_Mag ; } } }
	{ Name a ; Value { Term { [ {a} ]              ; In Domain_Mag ; } } }
    }
  }
}

PostProcessing {
    { Name MagDyn_avf ; NameOfFormulation MagDyn_avf ;
    PostQuantity {
	{ Name v ; Value { Term { [ {v} ]              ; In DomainC_Mag ; } } }
	{ Name h ; Value { Term { [ nu[] * {d a} ]     ; In Domain_Mag ; } } }
	{ Name a ; Value { Term { [ {a} ]              ; In Domain_Mag ; } } }
	{ Name j ; Value { Term { [ - sigma[]*(Dt[{a}]+ {d v}) ]; In 
DomainC_Mag ; } } }
	//{ Name V ; Value { Term { [ {V} ]              ; In Electrode_Mag ; } } }
	//{ Name I ; Value { Term { [ {I} ]              ; In Electrode_Mag ; } } }
	//{ Name Z ; Value { Term { [ -{V}/{I} ]         ; In Electrode_Mag ; } } }

    }
    }
}

========================================================================

/* Jacobian_Lib.pro */

Jacobian {
  { Name Vol ;
    Case { { Region All; Jacobian Vol; }
    }
  }
  { Name Sur ;
    Case { { Region All ; Jacobian Sur ; }
    }
  }
}

=====================================================================

Integration {
    { Name CurlCurl ;
    Case { {Type Gauss ;
    Case { { GeoElement Triangle    ; NumberOfPoints  4 ; }
    { GeoElement Quadrangle  ; NumberOfPoints  4 ; }
    { GeoElement Tetrahedron ; NumberOfPoints  4 ; }
    { GeoElement Hexahedron  ; NumberOfPoints  6 ; }
    { GeoElement Prism       ; NumberOfPoints  9 ; } }
    }
    }
    }
}


Thomas Jung
Fraunhofer-Institut IISB
D-91058 Erlangen, Schottkystr. 10
+49 9131 761264
-------------- next part --------------
lc = 1;
Point(1) = {0, 0, 0, 0.075};
Point(2) = {0, 0, 1, 0.075};
Point(3) = {0, 1, 0, 0.075};
Point(4) = {0, 1, 1, 0.075};
Point(5) = {1, 0, 0, 0.075};
Point(6) = {1, 0, 1, 0.075};
Point(7) = {1, 1, 0, 0.075};
Point(8) = {1, 1, 1, 0.075};
Point(9) = {0.2, 0.2, 0.4, 0.03};
Point(10) = {0.2, 0.2, 0.6, 0.03};
Point(11) = {0.2, 0.8, 0.4, 0.03};
Point(12) = {0.2, 0.8, 0.6, 0.03};
Point(13) = {0.8, 0.2, 0.4, 0.03};
Point(14) = {0.8, 0.2, 0.6, 0.03};
Point(15) = {0.8, 0.3, 0.4, 0.03};
Point(16) = {0.7, 0.3, 0.4, 0.03};
Point(17) = {0.5, 0.3, 0.4, 0.03};
Point(18) = {0.3, 0.3, 0.4, 0.03};
Point(19) = {0.3, 0.7, 0.4, 0.03};
Point(20) = {0.5, 0.7, 0.4, 0.03};
Point(21) = {0.7, 0.7, 0.4, 0.03};
Point(22) = {0.8, 0.7, 0.4, 0.03};
Point(23) = {0.8, 0.8, 0.4, 0.03};
Point(24) = {0.8, 0.8, 0.6, 0.03};
Point(25) = {0.8, 0.3, 0.6, 0.03};
Point(26) = {0.7, 0.3, 0.6, 0.03};
Point(27) = {0.5, 0.3, 0.6, 0.03};
Point(28) = {0.3, 0.3, 0.6, 0.03};
Point(29) = {0.3, 0.7, 0.6, 0.03};
Point(30) = {0.5, 0.7, 0.6, 0.03};
Point(31) = {0.7, 0.7, 0.6, 0.03};
Point(32) = {0.8, 0.7, 0.6, 0.03};
Point(33) = {0.4, 0.4, 0.3, 0.03};
Point(34) = {0.4, 0.6, 0.3, 0.03};
Point(35) = {0.6, 0.6, 0.3, 0.03};
Point(36) = {0.6, 0.4, 0.3, 0.03};
Point(37) = {0.4, 0.4, 0.7, 0.03};
Point(38) = {0.4, 0.6, 0.7, 0.03};
Point(39) = {0.6, 0.6, 0.7, 0.03};
Point(40) = {0.6, 0.4, 0.7, 0.03};
Line(1) = {1,3};
Line(2) = {3,4};
Line(3) = {4,1};
Line(4) = {4,2};
Line(5) = {2,1};
Line(6) = {5,1};
Line(7) = {2,5};
Line(8) = {2,6};
Line(9) = {6,5};
Line(10) = {5,3};
Line(11) = {5,7};
Line(12) = {7,3};
Line(13) = {7,4};
Line(14) = {7,8};
Line(15) = {8,4};
Line(16) = {4,6};
Line(17) = {8,6};
Line(18) = {5,8};
Line(19) = {9,11};
Line(20) = {11,10};
Line(21) = {10,9};
Line(22) = {12,10};
Line(23) = {11,12};
Line(24) = {13,9};
Line(25) = {9,14};
Line(26) = {14,13};
Line(27) = {10,14};
Line(28) = {9,17};
Line(29) = {17,18};
Line(30) = {18,9};
Line(31) = {13,17};
Line(32) = {18,19};
Line(33) = {19,9};
Line(34) = {11,19};
Line(35) = {19,20};
Line(36) = {20,11};
Line(37) = {20,23};
Line(38) = {23,11};
Line(39) = {16,13};
Line(40) = {13,15};
Line(41) = {15,16};
Line(42) = {16,17};
Line(43) = {21,22};
Line(44) = {22,23};
Line(45) = {23,21};
Line(46) = {20,21};
Line(47) = {11,24};
Line(48) = {24,23};
Line(49) = {12,24};
Line(50) = {10,27};
Line(51) = {27,28};
Line(52) = {28,10};
Line(53) = {14,27};
Line(54) = {28,29};
Line(55) = {29,10};
Line(56) = {12,29};
Line(57) = {29,30};
Line(58) = {30,12};
Line(59) = {30,24};
Line(60) = {26,14};
Line(61) = {14,25};
Line(62) = {25,26};
Line(63) = {26,27};
Line(64) = {31,32};
Line(65) = {32,24};
Line(66) = {24,31};
Line(67) = {30,31};
Line(68) = {15,25};
Line(69) = {25,13};
Line(70) = {16,25};
Line(71) = {16,26};
Line(72) = {17,27};
Line(73) = {27,16};
Line(74) = {18,28};
Line(75) = {28,17};
Line(76) = {24,22};
Line(77) = {32,22};
Line(78) = {19,28};
Line(79) = {19,29};
Line(80) = {21,32};
Line(81) = {21,31};
Line(82) = {29,20};
Line(83) = {30,20};
Line(84) = {30,21};
Line(85) = {34,33};
Line(86) = {33,36};
Line(87) = {36,34};
Line(88) = {36,35};
Line(89) = {35,34};
Line(90) = {34,38};
Line(91) = {38,33};
Line(92) = {38,37};
Line(93) = {37,33};
Line(94) = {38,35};
Line(95) = {38,39};
Line(96) = {39,35};
Line(97) = {39,36};
Line(98) = {39,40};
Line(99) = {40,36};
Line(100) = {37,36};
Line(101) = {37,40};
Line(102) = {40,38};
Line Loop(1) = {1,2,3};
Line Loop(2) = {4,5,-3};
Line Loop(3) = {6,-5,7};
Line Loop(4) = {8,9,-7};
Line Loop(5) = {-1,-6,10};
Line Loop(6) = {11,12,-10};
Line Loop(7) = {-2,-12,13};
Line Loop(8) = {14,15,-13};
Line Loop(9) = {-8,-4,16};
Line Loop(10) = {-15,17,-16};
Line Loop(11) = {-14,-11,18};
Line Loop(12) = {-9,-17,-18};
Line Loop(13) = {19,20,21};
Line Loop(14) = {22,-20,23};
Line Loop(15) = {24,25,26};
Line Loop(16) = {27,-25,-21};
Line Loop(17) = {28,29,30};
Line Loop(18) = {-28,-24,31};
Line Loop(19) = {-30,32,33};
Line Loop(20) = {34,35,36};
Line Loop(21) = {-34,-19,-33};
Line Loop(22) = {-36,37,38};
Line Loop(23) = {39,40,41};
Line Loop(24) = {-39,42,-31};
Line Loop(25) = {43,44,45};
Line Loop(26) = {-45,-37,46};
Line Loop(27) = {38,47,48};
Line Loop(28) = {49,-47,23};
Line Loop(29) = {50,51,52};
Line Loop(30) = {-50,27,53};
Line Loop(31) = {-52,54,55};
Line Loop(32) = {56,57,58};
Line Loop(33) = {-56,22,-55};
Line Loop(34) = {-58,59,-49};
Line Loop(35) = {60,61,62};
Line Loop(36) = {-60,63,-53};
Line Loop(37) = {64,65,66};
Line Loop(38) = {-66,-59,67};
Line Loop(39) = {40,68,69};
Line Loop(40) = {-61,26,-69};
Line Loop(41) = {41,70,-68};
Line Loop(42) = {-62,-70,71};
Line Loop(43) = {42,72,73};
Line Loop(44) = {-63,-71,-73};
Line Loop(45) = {29,74,75};
Line Loop(46) = {-51,-72,-75};
Line Loop(47) = {44,-48,76};
Line Loop(48) = {-65,77,-76};
Line Loop(49) = {32,78,-74};
Line Loop(50) = {-54,-78,79};
Line Loop(51) = {-43,80,77};
Line Loop(52) = {64,-80,81};
Line Loop(53) = {-35,79,82};
Line Loop(54) = {57,83,-82};
Line Loop(55) = {-46,-83,84};
Line Loop(56) = {67,-81,-84};
Line Loop(57) = {85,86,87};
Line Loop(58) = {88,89,-87};
Line Loop(59) = {-85,90,91};
Line Loop(60) = {92,93,-91};
Line Loop(61) = {89,90,94};
Line Loop(62) = {95,96,-94};
Line Loop(63) = {88,-96,97};
Line Loop(64) = {98,99,-97};
Line Loop(65) = {-86,-93,100};
Line Loop(66) = {101,99,-100};
Line Loop(67) = {92,101,102};
Line Loop(68) = {-98,-95,-102};
Plane Surface(1) = {1};
Plane Surface(2) = {2};
Plane Surface(3) = {3};
Plane Surface(4) = {4};
Plane Surface(5) = {5};
Plane Surface(6) = {6};
Plane Surface(7) = {7};
Plane Surface(8) = {8};
Plane Surface(9) = {9};
Plane Surface(10) = {10};
Plane Surface(11) = {11};
Plane Surface(12) = {12};
Plane Surface(13) = {13};
Plane Surface(14) = {14};
Plane Surface(15) = {15};
Plane Surface(16) = {16};
Plane Surface(17) = {17};
Plane Surface(18) = {18};
Plane Surface(19) = {19};
Plane Surface(20) = {20};
Plane Surface(21) = {21};
Plane Surface(22) = {22};
Plane Surface(23) = {23};
Plane Surface(24) = {24};
Plane Surface(25) = {25};
Plane Surface(26) = {26};
Plane Surface(27) = {27};
Plane Surface(28) = {28};
Plane Surface(29) = {29};
Plane Surface(30) = {30};
Plane Surface(31) = {31};
Plane Surface(32) = {32};
Plane Surface(33) = {33};
Plane Surface(34) = {34};
Plane Surface(35) = {35};
Plane Surface(36) = {36};
Plane Surface(37) = {37};
Plane Surface(38) = {38};
Plane Surface(39) = {39};
Plane Surface(40) = {40};
Plane Surface(41) = {41};
Plane Surface(42) = {42};
Plane Surface(43) = {43};
Plane Surface(44) = {44};
Plane Surface(45) = {45};
Plane Surface(46) = {46};
Plane Surface(47) = {47};
Plane Surface(48) = {48};
Plane Surface(49) = {49};
Plane Surface(50) = {50};
Plane Surface(51) = {51};
Plane Surface(52) = {52};
Plane Surface(53) = {53};
Plane Surface(54) = {54};
Plane Surface(55) = {55};
Plane Surface(56) = {56};
Plane Surface(57) = {57};
Plane Surface(58) = {58};
Plane Surface(59) = {59};
Plane Surface(60) = {60};
Plane Surface(61) = {61};
Plane Surface(62) = {62};
Plane Surface(63) = {63};
Plane Surface(64) = {64};
Plane Surface(65) = {65};
Plane Surface(66) = {66};
Plane Surface(67) = {67};
Plane Surface(68) = {68};
Surface Loop(1) = {1,2,3,4,5,6,7,8,9,10,11,12};
Surface Loop(2) = {-13,-14,-15,-16,-17,-18,-19,-20,-21,-22,-23,-24,-25,-26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,-51,-52,-53,-54,-55,-56};
Surface Loop(3) = {-57,-58,-59,-60,61,62,63,64,-65,-66,67,68};
Surface Loop(4) = {57,58,59,60,-61,-62,-63,-64,65,66,-67,-68};
Surface Loop(5) = {13,14,15,16,17,18,19,20,21,22,23,24,25,26,-27,-28,-29,-30,-31,-32,-33,-34,-35,-36,-37,-38,-39,-40,-41,-42,-43,-44,-45,-46,-47,-48,-49,-50,51,52,53,54,55,56};
Volume(1)={1,2,3};
Volume(3)={5};
Volume(4)={4};
Physical Surface (1) = {1,2};
Physical Surface (2) = {3,4};
Physical Surface (3) = {5,6};
Physical Surface (4) = {7,8};
Physical Surface (5) = {9,10};
Physical Surface (6) = {11,12};
Physical Surface (7) = {13,14};
Physical Surface (8) = {15,16};
Physical Surface (9) = {17,18,19,20,21,22,23,24,25,26};
Physical Surface (10) = {27,28};
Physical Surface (11) = {29,30,31,32,33,34,35,36,37,38};
Physical Surface (12) = {39,40};
Physical Surface (13) = {41,42};
Physical Surface (14) = {43,44};
Physical Surface (15) = {45,46};
Physical Surface (16) = {47,48};
Physical Surface (17) = {49,50};
Physical Surface (18) = {51,52};
Physical Surface (19) = {53,54};
Physical Surface (20) = {55,56};
Physical Surface (21) = {57,58};
Physical Surface (22) = {59,60};
Physical Surface (23) = {61,62};
Physical Surface (24) = {63,64};
Physical Surface (25) = {65,66};
Physical Surface (26) = {67,68};
Physical Volume (1) = {1};
Physical Volume (2) = {3};
Physical Volume (3) = {4};