Select Git revision
OCCEdge.cpp
microstrip.pro 12.28 KiB
/* -------------------------------------------------------------------
Tutorial 1 : electrostatic field of a microstrip
Features:
- Physical and abstract regions
- Scalar FunctionSpace with Dirichlet constraint
- Galerkin term for stiffness matrix
To compute the solution in a terminal:
getdp microstrip -solve EleSta_v
getdp microstrip -pos Map
getdp microstrip -pos Cut
To compute the solution interactively from the Gmsh GUI:
File > Open > microstrip.pro
Run (button at the bottom of the left panel)
------------------------------------------------------------------- */
/* In this first tutorial we consider the calculation of the electric field
given a static distribution of electric potential. This corresponds to an
"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 microstrip 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 the left boundary of the domain to account for the symmetry of the
problem, as well as on the top and right boundaries that truncate 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. */
Air = Region[101];
Diel1 = Region[111];
Ground = Region[120];
Electrode = Region[121];
/* We now define abstract regions to be used below in the definition of the
scalar electric potential formulation:
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* groups contain only volume elements of the mesh (triangles here).
Sur* groups contain only surface elements of the mesh (lines here). */
Vol_Ele = Region[ {Air, Diel1} ];
SurNeumann_Ele = Region[ {} ]; // empty for this model
}
Function {
/* 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
/* 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. */
{ Name Dirichlet_Ele; Type Assign;
Case {
{ Region Ground; Value 0.; }
{ Region Electrode; Value 1.e-3; }
}
}
}
Group{
/* 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.*/
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 (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(x,y) = Sum_k vn_k sn_k(x,y)
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; }
}
}
}
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 ; }
}
}
}
Integration {
/* A Gauss quadrature rule with 4 points is used for all integrations below. */
{ Name Int ;
Case { { Type Gauss ;
Case { { GeoElement Triangle ; NumberOfPoints 4 ; }
{ GeoElement Quadrangle ; NumberOfPoints 4 ; } }
}
}
}
}
Formulation {
/* 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
holds for all v'.
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
(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.
It has got 4 semicolon separated arguments:
* 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 integration method to be used.
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}").
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
[ epsilon[] * {d v} , basis_functions_of {d v} ].
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}".
The first argument, on the other hand, can contain a much wider variety of
expressions than the second one. In this problem, it contains
epsilon[] * {d v} = Sum_k vn_k epsilon[]*{d sn_k}
which is the eletric displacement vector, up to a sign. Here, we have two
valid syntaxes, with very different meanings.
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":
[ epsilon[] * Dof{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:
[ epsilon[] * {d v} , {d v} ].
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 { [ epsilon[] * Dof{d v} , {d v} ];
In Vol_Ele; Jacobian Vol; Integration Int; }
}
}
}
Resolution {
{ Name EleSta_v;
System {
{ Name Sys_Ele; NameOfFormulation Electrostatics_v; }
}
Operation {
Generate[Sys_Ele]; Solve[Sys_Ele]; SaveSolution[Sys_Ele];
}
}
}
/* 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 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
- 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). */
PostProcessing {
{ Name EleSta_v; NameOfFormulation Electrostatics_v;
Quantity {
{ Name v;
Value {
Local { [ {v} ];
In Dom_Hgrad_v_Ele; Jacobian Vol; }
}
}
{ Name e;
Value {
Local { [ -{d v} ];
In Dom_Hgrad_v_Ele; Jacobian Vol; }
}
}
{ Name d;
Value {
Local { [ -epsilon[] * {d v} ];
In Dom_Hgrad_v_Ele; Jacobian Vol; }
}
}
}
}
}
e = 1.e-7; // tolerance to ensure that the cut is inside the simulation domain
h = 2.e-3; // vertical position of the cut
PostOperation {
{ Name Map; NameOfPostProcessing EleSta_v;
Operation {
Print [ v, OnElementsOf Dom_Hgrad_v_Ele, File "mStrip_v.pos" ];
Print [ e, OnElementsOf Dom_Hgrad_v_Ele, File "mStrip_e.pos" ];
Print [ e, OnLine {{e,h,0}{14.e-3,h,0}}{60}, File "Cut_e.pos" ];
}
}
{ Name Cut; NameOfPostProcessing EleSta_v;
// same cut as above, with more points and exported in raw text format
Operation {
Print [ e, OnLine {{e,e,0}{14.e-3,e,0}} {500}, Format TimeTable, File "Cut_e.txt" ];
}
}
}