Skip to content
Snippets Groups Projects
Commit 4857d9a1 authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

rework explanations in electrostatic tutorial

parent 4e8a8e10
Branches
No related tags found
No related merge requests found
Pipeline #
Show whitespace changes
Inline Side-by-side
Electrostatics/microstrip.geo
+ 15
− 15
View file @ 4857d9a1
/* -------------------------------------------------------------------
File "microstrip.geo"
This file is the geometrical description used by Gmsh to produce
the mesh file "microstrip.msh" that is read by GetDP.
This file is the geometrical description used by Gmsh to produce the mesh
file "microstrip.msh", which is read by GetDP when analysing the file
"microstrip.pro".
------------------------------------------------------------------- */
/* Definition of some parameters for geometrical dimensions, i.e. h (height of
......@@ -40,19 +41,18 @@ Line(11) = {5,7};
/* Definition of geometrical surfaces */
Line Loop(12) = {1, 2, -8, -9, 7}; Plane Surface(13) = {12};
Line Loop(14) = {10,-11,8,3,4,5}; Plane Surface(15) = {14};
Curve Loop(12) = {1, 2, -8, -9, 7}; Plane Surface(13) = {12};
Curve Loop(14) = {10,-11,8,3,4,5}; Plane Surface(15) = {14};
/* Definition of Physical entities (surfaces, lines). The definition of
Physical entities (Surfaces and Lines) tells Gmsh the elements and associated
region numbers that have to be saved in the mesh file 'microstrip.msh'. For
example, Region 111 is made of the triangle elements of the geometric surface
13, whereas Region 121 is made of line elements of the geometric lines 9, 10
and 11. */
/* Definition of Physical entities. The definition of Physical entities
(Surfaces and Curves) tells Gmsh the elements and associated region numbers
that have to be saved in the mesh file 'microstrip.msh'. For example, Region
111 is made of the triangle elements of the geometric surface 13, whereas
Region 121 is made of line elements of the geometric curves 9, 10 and 11. */
Physical Surface("Air", 101) = {15};
Physical Surface("Dielectric", 111) = {13};
Physical Line ("Ground", 120) = {1} ;
Physical Line ("Electrode", 121) = {9,10,11} ;
Physical Line ("Surface infinity", 130) = {2,3,4} ;
Physical Curve("Ground", 120) = {1} ;
Physical Curve("Electrode", 121) = {9,10,11} ;
Physical Curve("Surface infinity", 130) = {2,3,4} ;
Electrostatics/microstrip.pro
+ 151
− 124
View file @ 4857d9a1
......@@ -16,48 +16,67 @@
Run (button at the bottom of the left panel)
------------------------------------------------------------------- */
/* In this first tutorial we consider the solution of electric fields given a
static distribution of electric potentials. This corresponds to a so-called
"electrostatic" physical model, obtained by combining the time-invariant
Maxwell-Ampere equation (Curl e = 0, with e the electric field) with Gauss'
law (Div d = rho, with d the displacement field and rho the charge density)
and the dielectric constitutive law (d = epsilon e, with epsilon the
dielectric permittivity).
Since Curl e = 0, e can be derived from a scalar electric potential v, such
that e = -Grad v. Plugging this potential in Gauss' law and using the
constitutive law leads to a scalar Poisson equation in terms of the electric
potential: -Div(epsilon Grad v) = rho.
We consider here the special case where rho = 0, to model a conducting
microstrip line on top of a dielectric substrate. A Dirichlet boundary
condition sets the potential to 1 mV on the boundary of the line (called
"Electrode" below) and 0 V on the ground. A homogeneous Neumann boundary
condition (zero flux of the displacement field, i.e. n.d = 0) is imposed on a
surface truncating the modelling domain. */
Group {
/* One starts by giving explicit meaningful names to
the Physical regions defined in the "microstrip.msh" mesh file.
There are 2 volume regions and 3 surface regions in this model. */
/* One starts by giving explicit meaningful names to the Physical regions
defined in the "microstrip.msh" mesh file. There are 2 volume regions and 3
surface regions in this model. */
Air = Region[101];
Diel1 = Region[111];
Ground = Region[120];
Electrode = Region[121];
SurfInf = Region[130];
/* We now define abstract regions to be used below
in the definition of the scalar electric potential formulation:
/* We now define abstract regions to be used below in the definition of the
scalar electric potential formulation:
Vol_Dielectric_Ele : dielectric volume regions where
"Div ( epsr[] Grad v)" is solved
Sur_Dir_Ele : Dirichlet boundary condition (v imposed)
Sur_Neu_Ele : Neumann bondary condition ( epsr[] n.Grad v = 0 )
Vol_Ele : volume where -Div(epsilon Grad v) = 0 is solved
SurNeumann_Ele : surface where non homogeneous Neumann boundary conditions
(on n.d == epsilon n.Grad v) are imposed
Vol_xxx groups contain only volume elements of the mesh (triangles here).
Sur_xxx groups contain only surface elements of the mesh (lines here).
*/
Vol* groups contain only volume elements of the mesh (triangles here).
Sur* groups contain only surface elements of the mesh (lines here). */
Vol_Dielectric_Ele = Region[ {Air, Diel1} ];
Sur_Dir_Ele = Region[ {Ground, Electrode} ];
Sur_Neu_Ele = Region[ {SurfInf} ];
Vol_Ele = Region[ {Air, Diel1} ];
SurNeumann_Ele = Region[ {} ]; // empty for this model
}
Function {
/* Material laws (here the relative permittivity)
are defined piecewise in terms of the above defined physical regions */
/* Constants are evaluated (only once) when the .pro file is parsed by GetDP,
before any finite element calculation takes place. */
eps0 = 8.854187818e-12; // permittivity of empty space
epsr[Air] = 1.;
epsr[Diel1] = 9.8;
/* Material laws (here the dielectric permittivity) are defined as piecewise
functions (note the square brackets), in terms of the above defined
physical regions. */
epsilon[Air] = 1. * eps0;
epsilon[Diel1] = 9.8 * eps0;
}
Constraint {
/* As for material laws, the Dirichlet boundary condition
is defined piecewise.
The constraint "Dirichlet_Ele" is invoked in the FunctionSpace below */
/* As for material laws, the Dirichlet boundary condition is defined
piecewise. The constraint "Dirichlet_Ele" is invoked in the FunctionSpace
below. */
{ Name Dirichlet_Ele; Type Assign;
Case {
{ Region Ground; Value 0.; }
......@@ -67,53 +86,60 @@ Constraint {
}
Group{
/* The domain of definition of a FunctionSpace lists all regions
on which a field is defined.
/* The domain of definition of a FunctionSpace lists all regions on which a
field is defined.
Contrary to Vol_xxx and Sur_xxx regions,
which contain only volume or surface regions, resp.,
domains of definition Dom_xxx regions may contain
both volume and surface regions.
Hence the use of the prefixes Vol_, Sur_ and Dom_ to avoid confusions.*/
Contrary to Vol_xxx and Sur_xxx regions, which contain only volume or
surface regions, resp., domains of definition Dom_xxx regions may contain
both volume and surface regions. Hence the use of the prefixes Vol_, Sur_
and Dom_ to avoid confusions.*/
Dom_Hgrad_v_Ele = Region[ {Vol_Dielectric_Ele, Sur_Dir_Ele, Sur_Neu_Ele} ];
Dom_Hgrad_v_Ele = Region[ {Vol_Ele, SurNeumann_Ele} ];
}
FunctionSpace {
/* The function space in which we shall pick
the electric scalar potential "v" solution
is definied by
- a domain of definition ("Dom_Hgrad_v_Ele")
/* The function space in which we shall pick the electric scalar potential "v"
solution is definied by
- a domain of definition (the "Support": "Dom_Hgrad_v_Ele")
- a type ("Form0" means scalar field)
- a set of scalar basis functions ("BF_Node" means nodal basis functions)
- a set of entities to which the basis functions are associated ("Entity":
here all the nodes of the domain of definition "NodesOf[All]")
- a constraint (here the Dirichlet boundary conditions)
The FE expansion of the unknown field "v" reads
v = Sum_k vn_k sn_k
v(x,y) = Sum_k vn_k sn_k(x,y)
where the "vn_k" are the nodal values (connectors)
and "sn_k" the basis functions.
Not all connectors are unknowns of the FE problem,
due to the "Constraint", which assigns particular values
to the nodes of the region "Sur_Dir_Ele".
GetDP deals with that automatically
on basis of the definition of the FunctionSpace.
where the "vn_k" are the nodal values (connectors) and "sn_k(x,y)" the
nodal basis functions. Not all connectors are unknowns of the FE problem,
due to the "Constraint", which assigns particular values to the nodes of
the Ground and Electrode regions. GetDP deals with that automatically on
basis of the definition of the FunctionSpace.
*/
{ Name Hgrad_v_Ele; Type Form0;
BasisFunction {
{ Name sn; NameOfCoef vn; Function BF_Node;
Support Dom_Hgrad_v_Ele; Entity NodesOf[ All ]; }
// using "NodesOf[All]" instead of "NodesOf[Dom_Hgrad_v_Ele]" is an
// optimization, which allows GetDP to not explicitly build the list of
// all nodes
}
Constraint {
{ NameOfCoef vn; EntityType NodesOf;
NameOfConstraint Dirichlet_Ele; }
{ NameOfCoef vn; EntityType NodesOf; NameOfConstraint Dirichlet_Ele; }
}
}
}
Jacobian {
/* Jacobians are used to specify the mapping between elements in the mesh and
the reference elements (in a standard unit-type cell) over which
integration is performed. "Vol" represents the classical 1-to-1 mapping
between elements of identical geometrical dimension, i.e. in this case
triangles/quadrangles in the plane onto the reference triangle/quadrangle
(2D to 2D). "Sur" would be used to map triangles/quadrangles in a 3D volume
onto the reference triangle/quadrangle (3D to 2D). "Lin" would be used to
map segments in a volume to the unit segment (3D to 1D). */
{ Name Vol ;
Case {
{ Region All ; Jacobian Vol ; }
......@@ -122,6 +148,7 @@ Jacobian {
}
Integration {
/* A Gauss quadrature rule with 4 points is used for all integrations below. */
{ Name Int ;
Case { { Type Gauss ;
Case { { GeoElement Triangle ; NumberOfPoints 4 ; }
......@@ -132,85 +159,93 @@ Integration {
}
Formulation {
/* The syntax of the Formulation{} section is a hard nut to crack.
So let's deal with it carefully.
/* The syntax of the Formulation section is a harder nut to crack. So let's
deal with it carefully.
A GetDP Formulation encodes a so-called weak formulation of the original
partial differential equation, i.e. of -Div(epsilon Grad v) = 0. This weak
formulation involves finding v such that
(-Div(epsilon Grad v) , v')_Vol_Ele = 0
holds for all so-called "test-functions" v', where (.,.) denotes an inner
product. If the test-functions v' are differentiable, integration by parts
using Green's identity leads to finding v such that
(epsilon Grad v, Grad v')_Vol_Ele + (epsilon n.Grad v, v')_SurNeumann_Ele = 0
The weak form of the electrostatic problem is particularly simple
in this model, as it has only one term:
holds for all v'.
(epsr[] Grad v, Grad vp)_Vol_Dielectric_Ele = 0
for all vp in S0(Vol_Dielectric_Ele).
In this particular example the Neumann boundary term vanishes and we are
looking for functions v in the function space Hgrad_v_Ele such that
The corresponding Euler-Lagrange equations are:
* Div ( epsr[] Grad v) = 0 on Vol_Dielectric_Ele
* epsr[] n.Grad v = 0 on Sur_Neu_Ele
(epsilon Grad v, Grad v')_Vol_Ele = 0
holds for all v', where the test functions v' are the same basis functions
("sn_k") as the ones used to interpolate the unknown potential v.
The Galerkin statement in the Formulation is a symbolic representation of
this weak formulation.
The Galerkin{} statement is a symbolic representation
of this weak formulation term.
It has got 4 semicolon separated arguments:
* the density [,] to be integrated,
* the density [.,.] to be integrated, with the test-functions (always) on
the right side of the comma
* the domain of integration,
* the jacobian of the transformation reference element -> real element,
* the Jacobian of the transformation reference element <-> real element,
* the integration method to be used.
The symbol "d" represents the exterior derivative,
and it is a synonym of "Grad" when applied to a scalar function.
The expression "d v" stands thus for the gradient of the v field,
i.e., for the electric field up to a minus sign.
What is a bit confusing is that the two comma-separated terms
of the bracket [,], in the first argument,
are not interpreted the same way. Let us unravel this in detail.
As the Galerkin method uses as trial (test) functions
the basis functions "sn_k" of the unknown field "v",
the density should be something like this
In the density, braces around a quantity (such as "{v}") indicate that this
quantity belongs to a FunctionSpace. Differential operators can be applied
within braces (such as "{Grad v}"); in particular the symbol "d" represents
the exterior derivative, and it is a synonym of "Grad" when applied to a
scalar function (such as "{d v}").
[ epsr[] * {d v} , basis_functions_of {d v} ].
What can be a bit confusing is that the two comma-separated terms of the
density are not interpreted exactly in the same way. Let us unravel this in
detail. As the Galerkin method uses as test functions the basis functions
"sn_k" of the unknown field "v", the density should be something like this
Since the second argument is devoted to the trial functions,
the operator "basis_functions_of" would always be there.
It can therefore be made implicit, and,
according to the GetDP syntax, it is omitted.
So, one writes simply "{d v}".
[ epsilon[] * {d v} , basis_functions_of {d v} ].
The first argument, on the other hand, can contain
a much wider variety of expressions than the second one.
In this problem, it contains
Since the second argument is always devoted to test functions, the operator
"basis_functions_of" would always be there. It can therefore be made
implicit, and, according to the GetDP syntax, it is omitted. So, one
writes simply "{d v}".
epsr[] * {d v} = Sum_k vn_k epsr[]*{d sn_k}
The first argument, on the other hand, can contain a much wider variety of
expressions than the second one. In this problem, it contains
which is the eletric displacement vector, up to a sign.
Here, we have two valid syntaxes, with very different meanings.
epsilon[] * {d v} = Sum_k vn_k epsilon[]*{d sn_k}
The first argument can be expressed in terms of
the FE expansion of "v" at the present system solution.
This is indicated by invoking the Dof{} operator:
which is the eletric displacement vector, up to a sign. Here, we have two
valid syntaxes, with very different meanings.
[ epsr[] * Dof{d v} , {d v} ].
The first argument can be expressed in terms of the FE expansion of "v" at
the present system solution, i.e. when the coefficients vn_k are unknown.
This is indicated by prefixing the braces with "Dof":
Another option, which would not work here,
is to evaluate the first argument
with the last available already computed solution.
For this the Dof{} operator is omitted:
[ epsilon[] * Dof{d v} , {d v} ].
[ epsr[] * {d v} , {d v} ].
Another option, which would not work here, is to evaluate the first
argument with the last available already computed solution, i.e. simply
perform the interpolation with known coefficients vn_k. For this the Dof
prefix operator would be omitted:
Both choices are commonly used in different contexts,
and we shall come back on this often in subsequent tutorials.
[ epsilon[] * {d v} , {d v} ].
To express the stiffness matrix of the electrostatic problem at hand,
we have to take the first option.
Hence the final expression of the density below.
Both choices are commonly used in different contexts, and we shall come
back on this often in subsequent tutorials.
*/
To express the stiffness matrix of the electrostatic problem at hand, we
have to take the first option. Hence the final expression of the density
below. */
{ Name Electrostatics_v; Type FemEquation;
Quantity {
{ Name v; Type Local; NameOfSpace Hgrad_v_Ele; }
}
Equation {
Galerkin { [ epsr[] * Dof{d v} , {d v} ];
In Vol_Dielectric_Ele;
Jacobian Vol; Integration Int; }
Galerkin { [ epsilon[] * Dof{d v} , {d v} ];
In Vol_Ele; Jacobian Vol; Integration Int; }
}
}
}
......@@ -226,27 +261,23 @@ Resolution {
}
}
eps0 = 8.854187818e-12; // permittivity of empty space
/* Post-processing is done in two parts.
The first part defines, in terms of the Formulation,
which itself refers to the FunctionSpace,
a number of quantities that can be evaluated at the postprocessing stage.
The three quantities defined here are :
The first part defines, in terms of the Formulation, which itself refers to
the FunctionSpace, a number of quantities that can be evaluated at the
postprocessing stage. The three quantities defined here are:
- the scalar vector potential,
- the electric field,
- the electric displacement.
The second part consists in defining groups of post-processing operations,
which can be invoked separately.
The first group is invoked by default when Gmsh is run interactively.
Each Operation specifies
which can be invoked separately. The first group is invoked by default when
Gmsh is run interactively. Each Operation specifies
- a quantity to be eveluated,
- the region on which the evaluation is done,
- the name of the output file.
The generated post-processing files are automatically displayed by Gmsh
if the "Merge result automatically" option is enabled
(which is the default). */
The generated post-processing files are automatically displayed by Gmsh if
the "Merge result automatically" option is enabled (which is the default). */
PostProcessing {
{ Name EleSta_v; NameOfFormulation Electrostatics_v;
......@@ -265,7 +296,7 @@ PostProcessing {
}
{ Name d;
Value {
Local { [ -eps0*epsr[] * {d v} ];
Local { [ -epsilon[] * {d v} ];
In Dom_Hgrad_v_Ele; Jacobian Vol; }
}
}
......@@ -281,10 +312,6 @@ PostOperation {
Operation {
Print [ v, OnElementsOf Dom_Hgrad_v_Ele, File "mStrip_v.pos" ];
Print [ e, OnElementsOf Dom_Hgrad_v_Ele, File "mStrip_e.pos" ];
Echo[ StrCat["l=PostProcessing.NbViews-1;",
"View[l].IntervalsType = 3;",
"View[l].NbIso = 40;"],
File "tmp.geo", LastTimeStepOnly] ;
Print [ e, OnLine {{e,h,0}{14.e-3,h,0}}{60}, File "Cut_e.pos" ];
}
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment