[Getdp] @D axisymmetric electrostatic problem

Nacho Andres nacho.andres at deos.tudelft.nl
Wed Apr 6 15:37:06 CEST 2005


Dear All,

I had a similar problem with an axisimmetric problem and what I imposed
is that the partial derivative of the magnitude to solve (in my case
this was the electric potential) should be zero at the axis of symmetry.
This can be easily done by adding the symmetry axis to the Basis
Functions in your Function Space, which would cause a null normal
derivative at that line to be imposed. I got a solution which looked
perfectly correct to me .. Am I right?

Good Luck,
Nacho

On Wed, 2005-04-06 at 15:16, Patrick Dular wrote:
> Gilles,
> 
> You are right, the radial field (x) component along the (y) vertical
> axis is not zero. Nevertheless, this component is negligible in
> comparison with the y-component. It is a good average of what happens
> in the vicinity of the axis.
> 
> The reason is the following. In your computation, you use first order
> finite elements for approximating the electric scalar potential. The
> electric field is computed as (minus) the gradient of the scalar
> potential and is therefore constant in each finite element. Each
> constant field, being an average of the electric field  in each
> element, can then have a radial component in some regions in case a
> significant radial exists next to the axis.
> 
> You can try to use the files in the attached  archive, in which a
> second order approximation is used for the scalar potential
> (homogeneous boundary conditions have been added for the additional
> degrees of freedom located on Dirichlet boundaries). A better accuracy
> will be obtained. With second order elements, you can reduce the
> number of elements in your mesh, otherwise the number of degrees of
> freedom can increase a lot (it will be then given by the number of
> nodes plus the number of edges). The bandwidth of the system matrix
> will be higher as well, so you can have to modify some parameters in
> your 'solver.par' file (e.g., to increase 'Nb_Fill') in order to allow
> convergence in the solving.
> 
> Patrick
> 
> 
> Gilles Quemener wrote:
> > Patrick,
> > 
> > Thanks a lot for looking into this problem. Due to the cylindrical
> > symmetry
> > around the vertical axis (Y), I am expecting to get the field
> > component along 
> > the horizonthal X (the radial component) such that E_r(@r=0) = 0 ?
> > But when 
> > I look in the file "Er" which contains the radial field component
> > within the
> >  chamber, I obtain  the following values for X = 0 :
> > 
> > # GetDP 1.0.0, binary
> > # PostData 'Er'
> > # Type Num  X Y Z  N1 N2 N3  Values  <Values>...
> > 15 66639  0 -0.04375 0  0 0 1  -23296.40465766675
> > 15 70830  0 -0.04275 0  0 0 1  -23443.70018452443
> > 15 64469  0 -0.04175 0  0 0 1  -23516.95331158148
> > 15 64469  0 -0.04074999999999999 0  0 0 1  -23516.95331158148
> > 15 65867  0 -0.03975 0  0 0 1  -21410.95426836388
> > 15 65867  0 -0.03875 0  0 0 1  -21410.95426836388
> > 15 62465  0 -0.03775 0  0 0 1  -19915.45461649752
> > 15 62465  0 -0.03675 0  0 0 1  -19915.45461649752
> > 15 58233  0 -0.03575 0  0 0 1  -16345.95241338546
> > 15 58233  0 -0.03475 0  0 0 1  -16345.95241338546
> > 15 67260  0 -0.03375 0  0 0 1  -21486.27224563957
> > 15 67260  0 -0.03274999999999999 0  0 0 1  -21486.27224563957
> > 15 59075  0 -0.03175 0  0 0 1  -20806.32373356267
> > 15 59075  0 -0.03075 0  0 0 1  -20806.32373356267
> > 15 69192  0 -0.02975 0  0 0 1  -17121.91130877584
> > 15 69192  0 -0.02875 0  0 0 1  -17121.91130877584
> > 15 69321  0 -0.02775 0  0 0 1  -13255.12227938629
> > 15 69321  0 -0.02675 0  0 0 1  -13255.12227938629
> > 15 62927  0 -0.02575 0  0 0 1  -13854.12868189462
> > 15 62927  0 -0.02475 0  0 0 1  -13854.12868189462
> > 15 55140  0 -0.02375 0  0 0 1  -14258.16683221585
> > 15 65487  0 -0.02275 0  0 0 1  -11406.89200786646
> > 15 65487  0 -0.02175 0  0 0 1  -11406.89200786646
> > 15 69991  0 -0.02075 0  0 0 1  -9523.366909591761
> > 15 69991  0 -0.01975 0  0 0 1  -9523.366909591761
> > 15 71120  0 -0.01875 0  0 0 1  -9543.905133518228
> > 15 71120  0 -0.01775 0  0 0 1  -9543.905133518228
> > 15 59775  0 -0.01675 0  0 0 1  -8941.018697394928
> > 15 59775  0 -0.01575 0  0 0 1  -8941.018697394928
> > .......................
> > etc....
> > This is a little puzzling for me !
> > 
> >        Gilles
> > 
> > Patrick Dular wrote:
> > > Gilles,
> > > 
> > > The solution I get for your problem looks correct regarding the
> > > distribution of the electric field.
> > > 
> > > In which region do you have an undesirable nonzero radial
> > > component of the electric field?
> > > 
> > > Patrick
> > > 
> > > Gilles Quemener wrote:
> > > > Hello Gmsh/GetDP users,
> > > > 
> > > > Currently I am trying to solve a simple 2D-axysymmetric
> > > > electrostatic problem.
> > > > In this problem,  due to the symmetry, I should obtain a zero
> > > > radial component
> > > > of the elctric field, but this is not what I get. could someone
> > > > help me finding my
> > > > mistake(s) ?
> > > > 
> > > > The problem is the following :  I have a cylindrical chamber
> > > > with a top electrode
> > > > set at a voltage of -1.8e05 V and a bottom one which is
> > > > grounded. The cylindrical
> > > > wall is made of a dielectric (quartz). See the attached picture
> > > > for a detailled geometry
> > > > where only half of the chamber is drawn and should be rotated
> > > > around the left vertical 
> > > > axis (corresponding to r=0) I also attached the .geo and
> > > > different .pro files that I took 
> > > > from some 2D examples provided on the website.
> > > > Thanks a lot for any help,
> > > > 
> > > >                    Gilles
> > > > 
> > > > 
> > > > ________________________________________________________________
> > > > /* -------------------------------------------------------------------
> > > >    File "mStrip.geo"
> > > > 
> > > >    This file is the geometrical description used by GMSH to produce
> > > >    the file "mStrip.msh".
> > > >    ------------------------------------------------------------------- */
> > > > 
> > > > /* 
> > > >    Definition of some parameters for geometrical dimensions, i.e.
> > > >    hd (height of 'Diel1'), td (thickness of dielectric 'Diel1')
> > > >    wge (width of grounded electrode 'Elec0'), tge (thickness of 'Elec0')
> > > >    wve (width of charged electrode 'ElecV'),  tve (thickness of 'ElecV')
> > > >    xBox (width of the air box) and yBox (height of the air box) 
> > > > */
> > > > cm = 1.e-02;
> > > > 
> > > > /* ********************************************************** */
> > > > /* ********************************************************** */
> > > > /* ********************************************************** */
> > > > /* chamber quartz wall */
> > > > hd = 15.*cm ; td = 1.5*cm ; red = 25.*cm; 
> > > > rid = red - td;
> > > > 
> > > > /* +V electrode (1) */
> > > > wve = 35.5*cm; tve = 9.85*cm ; 
> > > > epve = 3.*cm; hbve = 7.625*cm;
> > > > rcv = tve / 2. ;
> > > > 
> > > > /* Grounded electrode */
> > > > wge = 37.0*cm; tge = 16.1*cm ; rtge = 3.35*cm; 
> > > > epge = 3.*cm; hhge = 4.375*cm;
> > > > rcg = tge / 6. ;
> > > > retraitge = wge - wve;
> > > > rhge = wge - rcg - retraitge;
> > > > rbge = (rhge + red) / 2.;
> > > > 
> > > > /* +V electrode (2) */
> > > > rbve = rhge;
> > > > rhve = rbge;
> > > > 
> > > > /* curvature radius of all electrodes */
> > > > rcurv = td / 4.;
> > > > 
> > > > /* profondeur des gorges */
> > > > pg = 1.5*cm;
> > > > 
> > > > /* epaisseur alu */
> > > > epalu = 0.3*cm;
> > > > 
> > > > /* External volume */
> > > > xBox = (wge + wve)/2. * 4.5 ;  yBox = hd / 2. * 20. ;
> > > > 
> > > > /* Definition of parameters for local mesh dimensions */
> > > > 
> > > > s = 1. ;
> > > > p0 = hd / 10. * s ;
> > > > pGE = wge/2. / 10. * s ;  pCE = wve/2. / 10. * s ;
> > > > pmiddle = (pCE + pGE) / 2. * 10.;
> > > > pGEcorner = wge/2. / 100. * s ;  pCEcorner = wve/2. / 50. * s ;
> > > > pxBox = xBox / 10. * s ;  pyBox = yBox / 8. * s ;
> > > > 
> > > > /* ********************************************************** */
> > > > /* ********************************************************** */
> > > > /* ********************************************************** */
> > > > /* Definition of gemetrical points */
> > > > 
> > > > /* external box */
> > > > Point(1) = { 0              , -yBox        , 0, pyBox} ;
> > > > Point(2) = { xBox           , -yBox        , 0, pyBox} ;
> > > > Point(3) = { xBox           ,  0.          , 0, pyBox} ;
> > > > Point(4) = { xBox           ,  yBox        , 0, pyBox} ;
> > > > Point(5) = { 0              ,  yBox        , 0, pyBox} ;
> > > > Point(6) = { 0              ,  hbve + epve , 0, pCE} ;
> > > > Point(7) = { 0              ,  hbve        , 0, pGEcorner} ;
> > > > Point(8) = { 0              ,  0.          , 0, pGEcorner} ;
> > > > Point(9) = { 0              , -hhge        , 0, pGEcorner} ;
> > > > Point(10) = { 0             , -hhge - epge , 0, pGEcorner} ;
> > > > 
> > > > /* grounded electrode (1) */
> > > > Point(11) = { rbge          , -hhge - tge  , 0, pGEcorner} ;
> > > > Point(12) = { wge - rcg     , -hhge - tge  , 0, pGEcorner} ;
> > > > Point(13) = { wge - rcg     , -hhge - tge + rcg , 0, pGEcorner} ;
> > > > Point(14) = { wge           , -hhge - tge + rcg , 0, pGEcorner} ;
> > > > Point(15) = { wge           , -hhge - rcg  , 0, pGEcorner} ;
> > > > Point(16) = { wge - rcg     , -hhge - rcg  , 0, pGEcorner} ;
> > > > Point(17) = { wge - rcg     , -hhge        , 0, pGEcorner} ;
> > > > Point(18) = { rhge          , -hhge        , 0, pGEcorner} ;
> > > > Point(19) = { red + rcurv   , -hhge        , 0, pGEcorner} ;
> > > > Point(20) = { red + rcurv   , -hhge - rcurv, 0, pGEcorner} ;
> > > > Point(21) = { red           , -hhge - rcurv, 0, pGEcorner} ;
> > > > Point(22) = { red           , -hhge - pg   , 0, pGEcorner} ;
> > > > Point(23) = { red - td      , -hhge - pg   , 0, pGEcorner} ;
> > > > Point(24) = { red - td      , -hhge - rcurv, 0, pGEcorner} ;
> > > > Point(25) = { red - td - rcurv , -hhge - rcurv , 0, pGEcorner} ;
> > > > Point(26) = { red - td - rcurv , -hhge     , 0, pGEcorner} ;
> > > > Point(27) = { (red +rtge)/2., -hhge        , 0, pGEcorner} ;
> > > > Point(28) = { rtge + rcurv  , -hhge        , 0, pGEcorner} ;
> > > > Point(29) = { rtge + rcurv  , -hhge - rcurv, 0, pGEcorner} ;
> > > > Point(30) = { rtge          , -hhge - rcurv, 0, pGEcorner} ;
> > > > Point(31) = { rtge          , -hhge - epge , 0, pGEcorner} ;
> > > > 
> > > > /* chamber wall external middle point */
> > > > Point(32) = { red           ,  0   , 0, pGEcorner} ;
> > > > 
> > > > /* +V electode (1/3) */
> > > > Point(33) = { red           , hbve + rcurv, 0, pGEcorner} ;
> > > > Point(34) = { red           , hbve + pg   , 0, pGEcorner} ;
> > > > Point(35) = { red - td      , hbve + pg   , 0, pGEcorner} ;
> > > > Point(36) = { red - td      , hbve + rcurv, 0, pGEcorner} ;
> > > > 
> > > > /* chamber wall internal middle point */
> > > > Point(37) = { red - td      ,  0   , 0, pGEcorner} ;
> > > > 
> > > > /* +V electode (2/3) */
> > > > Point(38) = { red - td - rcurv , hbve + rcurv, 0, pGEcorner} ;
> > > > Point(39) = { red - td - rcurv , hbve        , 0, pGEcorner} ;
> > > > Point(40) = { (red - td)/2.    , hbve        , 0, pGEcorner} ;
> > > > Point(41) = { red + rcurv      , hbve        , 0, pGEcorner} ;
> > > > Point(42) = { red + rcurv      , hbve + rcurv, 0, pGEcorner} ;
> > > > Point(43) = { rbve             , hbve        , 0, pGEcorner} ;
> > > > Point(44) = { rbve             , hbve + rcv  , 0, pGEcorner} ;
> > > > Point(45) = { rbve + rcv       , hbve + rcv  , 0, pGEcorner} ;
> > > > Point(46) = { rbve             , hbve+ 2.*rcv, 0, pGEcorner} ;
> > > > Point(47) = { rhve             , hbve+ epve  , 0, pGEcorner} ;
> > > > 
> > > > /* grounded electrode (2) */
> > > > Point(48) = { rbge          , -hhge - epge  , 0, pGEcorner} ;
> > > > Point(49) = { rhge          , -hhge - epalu , 0, pGEcorner} ;
> > > > Point(50) = { wge - rcg     , -hhge - epalu , 0, pGEcorner} ;
> > > > Point(51) = { wge - epalu   , -hhge - rcg  , 0, pGEcorner} ;
> > > > Point(52) = { wge - epalu   , -hhge - tge + rcg , 0, pGEcorner} ;
> > > > Point(53) = { wge - rcg     , -hhge - tge + epalu , 0, pGEcorner} ;
> > > > Point(54) = { rbge          , -hhge - tge + epalu , 0, pGEcorner} ;
> > > > 
> > > > /* +V electode (3/3) */
> > > > Point(55) = { rbve                , hbve + epalu     , 0, pGEcorner} ;
> > > > Point(56) = { rbve + rcv - epalu  , hbve + rcv       , 0, pGEcorner} ;
> > > > Point(57) = { rbve             , hbve+ 2.*rcv - epalu, 0, pGEcorner} ;
> > > > 
> > > > /* grounded valve */
> > > > Point(58) = { 1.5*rtge         , -hhge - epge     , 0, pGEcorner};
> > > > Point(59) = { 1.5*rtge         , -hhge - 1.5*epge , 0, pGEcorner};
> > > > Point(60) = { 0.               , -hhge - 1.5*epge , 0, pGEcorner};
> > > > 
> > > > /* Additional points for meshing characteristic length */
> > > > Point(61) = { 0             , -hhge - 2.*epge , 0, pGEcorner} ;
> > > > Point(62) = { 0             , -hhge - 4.*epge , 0, pGEcorner} ;
> > > > Point(63) = { 0             , -hhge - 6.*epge , 0, pGEcorner} ;
> > > > 
> > > > 
> > > > /* ********************************************************** */
> > > > /* ********************************************************** */
> > > > /* ********************************************************** */
> > > > /* Definition of gemetrical lines */
> > > > 
> > > > Line(1) = {1,2};
> > > > Line(2) = {2,3};
> > > > Line(3) = {3,4};
> > > > Line(4) = {4,5};
> > > > Line(5) = {5,6};
> > > > Line(7) = {7,8};
> > > > Line(8) = {8,9};
> > > > Line(9) = {9,10};
> > > > Line(10) = {60,61,63,1};
> > > > 
> > > > Line(11) = {60,59};
> > > > Line(12) = {59,58};
> > > > Line(13) = {58,48};
> > > > Line(14) = {48,49};
> > > > Line(15) = {49,50};
> > > > Circle(16) = {50,16,51};
> > > > Line(17) = {51,52};
> > > > Circle(18) = {52,13,53};
> > > > Line(19) = {53,54};
> > > > Line(20) = {54,11};
> > > > Line(21) = {11,12};
> > > > Circle(22) = {12,13,14};
> > > > Line(23) = {14,15};
> > > > Circle(24) = {15,16,17};
> > > > Line(25) = {17,18};
> > > > Line(26) = {18,19};
> > > > Circle(27) = {19,20,21};
> > > > Line(28) = {21,22};
> > > > Line(29) = {22,23};
> > > > Line(30) = {23,24};
> > > > Circle(31) = {24,25,26};
> > > > Line(32) = {26,27};
> > > > Line(33) = {27,28};
> > > > Circle(34) = {28,29,30};
> > > > Line(35) = {30,31};
> > > > Line(36) = {31,10};
> > > > 
> > > > Line(37) = {21,32};
> > > > Line(38) = {32,33};
> > > > Line(39) = {36,37};
> > > > Line(40) = {37,24};
> > > > 
> > > > Line(41) = {6,47};
> > > > Line(42) = {47,55};
> > > > Circle(43) = {55,44,56};
> > > > Circle(44) = {56,44,57};
> > > > Line(45) = {57,46};
> > > > Circle(46) = {46,44,45};
> > > > Circle(47) = {45,44,43};
> > > > Line(48) = {43,41};
> > > > Circle(49) = {41,42,33};
> > > > Line(50) = {33,34};
> > > > Line(51) = {34,35};
> > > > Line(52) = {35,36};
> > > > Circle(53) = {36,38,39};
> > > > Line(54) = {39,40};
> > > > Line(55) = {40,7};
> > > > 
> > > > 
> > > > /* ********************************************************** */
> > > > /* ********************************************************** */
> > > > /* ********************************************************** */
> > > > /* Definition of geometrical surfaces */
> > > > 
> > > > /* external vacuum */
> > > > Line Loop(30) = {1,2,3,4,5,41,42,43,44,45,46,47,48,49,-38,-37,-27,-26,-25,-24,-23,-22,-21,-20,-19,-18,-17,-16,-15,-14,-13,-12,-11,10};
> > > > Plane Surface(31) = {30};
> > > > 
> > > > /* quartz wall */
> > > > Line Loop(32) = {37,38,50,51,52,39,40,-30,-29,-28};
> > > > Plane Surface(33) = {32};
> > > > 
> > > > /* internal vacuum */
> > > > Line Loop(34) = {7,8,9,-36,-35,-34,-33,-32,-31,-40,-39,53,54,55};
> > > > Plane Surface(35) = {34};
> > > > 
> > > > /* ********************************************************** */
> > > > /* ********************************************************** */
> > > > /* ********************************************************** */
> > > > /* Definition of Physical entities (surfaces, lines). The Physical
> > > >    entities tell GMSH the elements and their associated region numbers
> > > >    to save in the file 'mStrip.msh'. For example, the Region 
> > > >    111 is made of elements of surface 13, while the Region 121 is 
> > > >    made of elements of lines 9, 10 and 11 */
> > > > 
> > > > Physical Surface (101) = {31} ;   /* External Vacuum */
> > > > Physical Surface (102) = {33} ;   /* Quartz */
> > > > Physical Surface (103) = {35} ;   /* Internal Vacuum */
> > > > 
> > > > 
> > > > Physical Line (120) = {11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36} ;  /* Ground */
> > > > Physical Line (121) = {41,42,43,44,45,46,47,48,49,50,51,52,53,54,55} ; /* Elctrode V */
> > > > Physical Line (130) = {1,2,3,4} ;   /* SurfInf */
> > > > 
> > > > ________________________________________________________________
> > > > /* -------------------------------------------------------------------
> > > >    File "mStrip.pro"
> > > > 
> > > >    This file defines the problem dependent data structures for the
> > > >    microstrip problem.
> > > >    
> > > >    To compute the solution: getdp mStrip -solve EleSta_v
> > > >    To compute post-results: getdp mStrip -pos Map
> > > >                          or getdp mStrip -pos Cut
> > > >    ------------------------------------------------------------------- */
> > > > 
> > > > Group {
> > > >    
> > > >   /* Let's start by defining the interface (i.e. elementary groups)
> > > >      between the mesh file and GetDP (no mesh object is defined, so 
> > > >      the default mesh will be assumed to be in GMSH format and located 
> > > >      in "mChambre.msh") */
> > > > 
> > > >   EVacuum = Region[101] ; Quartz = Region[102] ; IVacuum = Region[103] ;
> > > >   Ground = Region[120] ; ElectrodeV = Region[121]; 
> > > >   SurfInf = Region[130];
> > > > 
> > > >   /* We can then define a global group (used in "EleSta_vl.pro",
> > > >      the file containing the function spaces and formulations) */
> > > > 
> > > >   DomainCC_Ele = Region[{EVacuum, Quartz, IVacuum}] ;
> > > > 
> > > > }
> > > > 
> > > > Function {
> > > > 
> > > >   /* The relative permittivity (needed in the formulation) is piecewise
> > > >      defined in elementary groups */
> > > > 
> > > >   epsr[EVacuum] = 1. ;
> > > >   epsr[IVacuum] = 1. ;
> > > > /*  epsr[Quartz] = 4.0 ;*/
> > > >   epsr[Quartz] = 3.8 ;
> > > > /*  epsr[Quartz] = 4.9 ;*/
> > > > 
> > > > }
> > > > 
> > > > Constraint {
> > > > 
> > > >   /* Now, some Dirichlet conditions are defined. The name 
> > > >      'ElectricScalarPotential' refers to the constraint name given in
> > > >      the function space */
> > > > 
> > > >   { Name ElectricScalarPotential ; Type Assign ;
> > > >     Case {
> > > >       { Region Region[{Ground, SurfInf}] ; Value 0. ; }
> > > >       { Region ElectrodeV ; Value -1.8e05 ; }
> > > >     }
> > > >   }
> > > > 
> > > > }
> > > > 
> > > > /* The formulation used and its tools, considered as being 
> > > >    in a black box, can now be included */
> > > > 
> > > > Include "Jacobian_Lib.pro"
> > > > Include "Integration_Lib.pro"
> > > > Include "EleSta_vl.pro"
> > > > 
> > > > /* Finally, we can define some operations to output results */
> > > > 
> > > > /* e = 1.e-7 ;*/
> > > > PostOperation {
> > > >   { Name Map ; NameOfPostProcessing EleSta_v ;
> > > >      Operation {
> > > >        Print [ v, OnElementsOf DomainCC_Ele, File "mChambre_v.pos" ] ;
> > > >        Print [ E, OnElementsOf DomainCC_Ele, File "mChambre_e.pos" ] ;
> > > >        Print [ modE, OnElementsOf DomainCC_Ele, File "mChambre_mode.pos" ] ;
> > > >        Print [ Er, OnElementsOf DomainCC_Ele, File "mChambre_ex.pos" ] ;
> > > >        Print [ Ez, OnElementsOf DomainCC_Ele, File "mChambre_ey.pos" ] ;
> > > >        Print [ th_EB, OnElementsOf DomainCC_Ele, File "mChambre_theb.pos" ] ;
> > > >      }
> > > >   }
> > > >   { Name Cut ; NameOfPostProcessing EleSta_v ;
> > > >      Operation {
> > > >        Print [ E, OnLine {{0,0,0}{0.5,0,0}} {500}, File "Cut_e", Format Gnuplot ] ;
> > > >        Print [ v, OnLine {{0.007,-0.5,0}{0.007,0.5,0}} {500}, File "Cut_v", Format Gnuplot ] ;
> > > >        Print [ E, OnLine {{0.007,-0.5,0}{0.007,0.5,0}} {500}, File "Cut_e2", Format Gnuplot ] ;
> > > >        Print [ Er, OnLine {{0.007,-0.5,0}{0.007,0.5,0}} {500}, File "Cut_Er", Format Gnuplot ] ;
> > > >        Print [ Ez, OnLine {{0.007,-0.5,0}{0.007,0.5,0}} {500}, File "Cut_Ez", Format Gnuplot ] ;
> > > >        Print [ Er, OnPlane {{0.,-0.04375,0}{0.235,-0.04375,0}{0.,0.07625,0}} {235,120}, File "Er", Format Gnuplot ] ;
> > > >        Print [ Ez, OnPlane {{0.,-0.04375,0}{0.235,-0.04375,0}{0.,0.07625,0}} {235,120}, File "Ez", Format Gnuplot ] ;
> > > >        Print [ th_EB, OnPlane {{0.,-0.04375,0}{0.235,-0.04375,0}{0.,0.07625,0}} {235,120}, File "th_EB", Format Gnuplot ] ;
> > > >      }
> > > >   }
> > > > 
> > > > }
> > > >   
> > > > 
> > > > ________________________________________________________________
> > > > /* -------------------------------------------------------------------
> > > >    File "EleSta_v.pro"
> > > > 
> > > >    Electrostatics - Electric scalar potential v formulation
> > > >    ------------------------------------------------------------------- 
> > > > 
> > > >    I N P U T
> > > >    ---------
> > > > 
> > > >    Global Groups :  (Extension '_Ele' is for Electric problem)
> > > >    -------------
> > > >    Domain_Ele               Whole electric domain (not used)
> > > >    DomainCC_Ele             Nonconducting regions
> > > >    DomainC_Ele              Conducting regions (not used)
> > > > 
> > > >    Function :
> > > >    --------
> > > >    epsr[]                   Relative permittivity
> > > > 
> > > >    Constraint :
> > > >    ----------
> > > >    ElectricScalarPotential  Fixed electric scalar potential
> > > >                             (classical boundary condition)
> > > > 
> > > >    Physical constants :
> > > >    ------------------                                               */
> > > > 
> > > >    eps0 = 8.854187818e-12 ;
> > > > 
> > > > /* O U T P U T
> > > >    -----------
> > > > 
> > > >   PostQuantities :
> > > >   --------------
> > > >    v : Electric scalar potential
> > > >    e : Electric field
> > > >    d : Electric flux density
> > > > 
> > > > */
> > > > 
> > > > Group {
> > > >   DefineGroup[ Domain_Ele, DomainCC_Ele, DomainC_Ele ] ;
> > > > }
> > > > 
> > > > Function {
> > > >   DefineFunction[ epsr ] ;
> > > > }
> > > > 
> > > > FunctionSpace {
> > > >   { Name Hgrad_v_Ele ; Type Form0 ;
> > > >     BasisFunction {
> > > >       // v = v  s   ,  for all nodes
> > > >       //      n  n
> > > >       { Name sn ; NameOfCoef vn ; Function BF_Node ;
> > > >         Support DomainCC_Ele ; Entity NodesOf[ All ] ; }
> > > >     }
> > > >     Constraint {
> > > >       { NameOfCoef vn ; EntityType NodesOf ; 
> > > >         NameOfConstraint ElectricScalarPotential ; }
> > > >     }
> > > >   }
> > > > }
> > > > 
> > > > 
> > > > Formulation {
> > > >   { Name Electrostatics_v ; Type FemEquation ;
> > > >     Quantity {
> > > >       { Name v ; Type Local ; NameOfSpace Hgrad_v_Ele ; }
> > > >     }
> > > >     Equation {
> > > >       Galerkin { [ epsr[] * Dof{d v} , {d v} ] ; In DomainCC_Ele ; 
> > > >                  Jacobian VolAxi ; Integration GradGrad ; }
> > > >     }
> > > >   }
> > > > }
> > > > 
> > > > 
> > > > Resolution {
> > > >   { Name EleSta_v ;
> > > >     System {
> > > >       { Name Sys_Ele ; NameOfFormulation Electrostatics_v ; }
> > > >     }
> > > >     Operation { 
> > > >       Generate Sys_Ele ; Solve Sys_Ele ; SaveSolution Sys_Ele ; 
> > > >     }
> > > >   }
> > > > }
> > > > 
> > > > 
> > > > PostProcessing {
> > > >   { Name EleSta_v ; NameOfFormulation Electrostatics_v ;
> > > >     Quantity {
> > > >       { Name v ; 
> > > >         Value { 
> > > >           Local { [ {v} ] ; In DomainCC_Ele ; Jacobian VolAxi ; } 
> > > >         }
> > > >       }
> > > >       { Name E ; 
> > > >         Value { 
> > > >           Local { [ -{d v} ] ; In DomainCC_Ele ; Jacobian VolAxi ; }
> > > >         }
> > > >       }
> > > >       { Name Er ; 
> > > >         Value { 
> > > >           Local { [CompX[ -{d v} ]] ; In DomainCC_Ele ; Jacobian VolAxi ; }
> > > >         }
> > > >       }
> > > >       { Name Ez ; 
> > > >         Value { 
> > > >           Local { [CompY[ -{d v} ]] ; In DomainCC_Ele ; Jacobian VolAxi ; }
> > > >         }
> > > >       }
> > > >       { Name d ; 
> > > >         Value { 
> > > >           Local { [ -eps0*epsr[] * {d v} ] ; In DomainCC_Ele ; Jacobian VolAxi ; } 
> > > >         } 
> > > >       }
> > > >       { Name th_EB ; 
> > > >         Value { 
> > > >           Local { [Atan2[CompX[ -{d v} ],CompY[ -{d v} ]]] ; In DomainCC_Ele ; Jacobian VolAxi ; }
> > > >         }
> > > >       }
> > > >       { Name modE ; 
> > > >         Value { 
> > > >           Local { [Norm[ -{d v} ]] ; In DomainCC_Ele ; Jacobian VolAxi ; }
> > > >         }
> > > >       }
> > > > /*
> > > >       { Name th_EB ; 
> > > >         Value { 
> > > >           Local { [Atan2[Fabs[CompX[ -{d v} ]],Fabs[CompY[ -{d v} ]]]] ; In DomainCC_Ele ; Jacobian VolAxi ; }
> > > >         }
> > > >       }
> > > > */
> > > >     }
> > > >   }
> > > > }
> > > >   
> > > > 
> > > > ________________________________________________________________
> > > > /*
> > > >   Jacobian methods
> > > >     VolAxi
> > > > */
> > > > 
> > > > /* I N P U T
> > > >    ---------
> > > > 
> > > >   GlobalGroup :
> > > >   -----------
> > > >     DomainInf                Regions with Spherical Shell Transformation
> > > >     DomainInfRectX           Regions with Rectangular Transformation in X direction
> > > >     DomainInfRectY           Regions with Rectangular Transformation in Y direction
> > > >     DomainInfRectZ           Regions with Rectangular Transformation in Z direction
> > > > 
> > > >   Parameters :
> > > >   ----------
> > > >     Val_Rint, Val_Rext       Inner and outer radius of the Spherical Shell
> > > >                              of DomainInf
> > > >     Val_Xint, Val_Xext       Inner and outer coordinates for 
> > > >                              rectangular transformation with axis X
> > > >     Val_Yint, Val_Yext       idem axis Y
> > > >     Val_Zint, Val_Zext       idem axis Z
> > > > 
> > > > */
> > > > 
> > > > /* --------------------------------------------------------------------------*/
> > > > 
> > > > Group {
> > > >   DefineGroup[ DomainInf ] ;
> > > >   DefineVariable[ Val_Rint, Val_Rext ] ;
> > > >   DefineGroup[ DomainInfRectX, DomainInfRectY, DomainInfRectZ ] ;
> > > >   DefineVariable[ Val_Xint, Val_Xext, Val_Yint, Val_Yext, Val_Zint, Val_Zext ] ;
> > > > }
> > > > 
> > > > /* --------------------------------------------------------------------------*/
> > > > 
> > > > Jacobian {
> > > >   {
> > > >     Name VolAxi ;
> > > >     Case {
> > > >       { Region DomainInf ;
> > > >         Jacobian VolAxiSphShell {Val_Rint, Val_Rext} ; }
> > > >       { Region DomainInfRectX ;
> > > >         Jacobian VolAxiRectShell {Val_Xint, Val_Xext, 1} ; }
> > > >       { Region DomainInfRectY ;
> > > >         Jacobian VolAxiRectShell {Val_Yint, Val_Yext, 2} ; }
> > > >       { Region DomainInfRectZ ;
> > > >         Jacobian VolAxiRectShell {Val_Zint, Val_Zext, 3} ; }
> > > >       { Region All ; Jacobian VolAxi ; }
> > > >     }
> > > >   }
> > > >   {
> > > >     Name SurAxi ;
> > > >     Case {
> > > >       { Region All ; Jacobian SurAxi ; }
> > > >     }
> > > >   }
> > > > }
> > > > 
> > > > /* --------------------------------------------------------------------------*/
> > > >   
> > > > 
> > > > ________________________________________________________________
> > > > /*
> > > >   Integration method
> > > >     GradGrad
> > > >     CurlCurl
> > > > */
> > > > 
> > > > 
> > > > Integration {
> > > >   { Name GradGrad ;
> > > >     Case { {Type Gauss ;
> > > >             Case { { GeoElement Triangle    ; NumberOfPoints  4 ; }
> > > >                    { GeoElement Quadrangle  ; NumberOfPoints  4 ; }
> > > >                    { GeoElement Tetrahedron ; NumberOfPoints  4 ; }
> > > >                    { GeoElement Hexahedron  ; NumberOfPoints  6 ; }
> > > >                    { GeoElement Prism       ; NumberOfPoints  9 ; } }
> > > >            }
> > > >          }
> > > >   }
> > > >   { 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 ; } }
> > > >            }
> > > >          }
> > > >   }
> > > > 
> > > >   {
> > > >     Name Sur ;
> > > >     Case {
> > > >       {
> > > > 	Type Gauss ;
> > > > 	Case {
> > > > 	  { GeoElement Line ; NumberOfPoints  3 ; }
> > > > 	}
> > > >       }
> > > >     }
> > > >   }
> > > > 
> > > > }
> > > > 
> > > > 
> > > > /* --------------------------------------------------------------------------*/
> > > > /* --------------------------------------------------------------------------*/
> > > >   
> > > > 
> > > > ________________________________________________________________
> > > > _______________________________________________
> > > > getdp mailing list
> > > > getdp at geuz.org
> > > > http://www.geuz.org/mailman/listinfo/getdp
> > > >   
> > > 
> > > -- 
> > > Patrick Dular, Dr. Ir., Research associate, F.N.R.S.
> > > Department of Electrical Engineering and Computer Science
> > > Unit of Applied Electricity
> > > University of Liege - Montefiore Institute - B28 - Parking P32
> > > B-4000 Liege - Belgium - Tel. +32-4 3663710 - Fax +32-4 3662910
> > > E-mail: Patrick.Dular at ulg.ac.be
> 
> -- 
> Patrick Dular, Dr. Ir., Research associate, F.N.R.S.
> Department of Electrical Engineering and Computer Science
> Unit of Applied Electricity
> University of Liege - Montefiore Institute - B28 - Parking P32
> B-4000 Liege - Belgium - Tel. +32-4 3663710 - Fax +32-4 3662910
> E-mail: Patrick.Dular at ulg.ac.be
> 
> ______________________________________________________________________
> _______________________________________________
> getdp mailing list
> getdp at geuz.org
> http://www.geuz.org/mailman/listinfo/getdp