Skip to content
Snippets Groups Projects
Commit 7b2c9bfe authored by François Henrotte's avatar François Henrotte
Browse files

Merge branch 'Electrostatic' 

Second tutorial added and minor changes in tuto 1
parents b7642cb6 a78f46f2
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Electrostatic/microstrip.pro
+ 98
− 27
View file @ 7b2c9bfe
/* -------------------------------------------------------------------
File "microstrip.pro"
Tutorial 1 : electrostatic field of a microstrip
This file defines the problem dependent data structures for the
microstrip problem.
Features:
- Physical and abstract regions
- Scalar FunctionSpace with Dirichlet constraint
- Galerkin term for stiffness matrix
To compute the solution:
To compute the solution in a terminal:
getdp microstrip -solve EleSta_v
To compute post-processing results:
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)
The microstrip.pro file describes the finite element (FE) model.
------------------------------------------------------------------- */
Group {
/* One starts by giving explicit meaningful names to
......@@ -26,12 +30,12 @@ Group {
Electrode = Region[121];
SurfInf = Region[130];
/* We now define abstract regions to be used 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 = 0" is solved
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 )
Sur_Neu_Ele : Neumann bondary condition ( epsr[] n.Grad v = 0 )
Vol_xxx groups contain only volume elements of the mesh (triangles here).
Sur_xxx groups contain only surface elements of the mesh (lines here).
......@@ -43,14 +47,15 @@ Group {
}
Function {
/* Material laws (here the relative permittivity)
are defined piecewise in terms of the above defined physical regions */
/* The relative permittivity is defined piecewise using the above defined Groups */
epsr[Air] = 1.;
epsr[Diel1] = 9.8;
}
Constraint {
/* As for material laws, the Dirichlet boundary conditions can be defined piecewise
/* 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;
......@@ -61,22 +66,36 @@ Constraint {
}
}
Group{
/* The domain of definition of a FunctionSpace lists all regions
on which a field is defined.
/* The function space of the electric scalar potential v
is definied by
- a domain of definition, which is the Group "Dom_Hgrad_v_Ele",
- a type (Form0 means scalar field)
- a set of scalar basis functions (here nodal basis functions "BF_Node")
- constraints (here the Dirichlet boundary conditions)
Contrary to the above defined groups, which contain either volure or surface elements,
Dom_xxx groups contain both volume and surface elements */
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.*/
Group{
Dom_Hgrad_v_Ele = Region[ {Vol_Dielectric_Ele, Sur_Dir_Ele, Sur_Neu_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")
- a type ("Form0" means scalar field)
- a set of scalar basis functions ("BF_Node" means nodal basis functions)
- a constraint (here the Dirichlet boundary conditions)
The FE representation of the unknown field "v" reads
v = Sum_k sn_k vn_k
where the "vn_k" are the nodal values and "sn_k" the basis functions.
Not all nodal values 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.
*/
{ Name Hgrad_v_Ele; Type Form0;
BasisFunction {
{ Name sn; NameOfCoef vn; Function BF_Node;
......@@ -98,7 +117,7 @@ Jacobian {
}
Integration {
{ Name GradGrad ;
{ Name Int ;
Case { {Type Gauss ;
Case { { GeoElement Triangle ; NumberOfPoints 4 ; }
{ GeoElement Quadrangle ; NumberOfPoints 4 ; } }
......@@ -108,6 +127,41 @@ Integration {
}
Formulation {
/* The syntax of the Formulation{} section is a hard nut to crack.
So let's deal with it carefully.
The weak form of the electrostatic problem is particularly simple in this model,
as it has only one term:
(epsr[] Grad v, Grad vp)_Vol_Dielectric_Ele = 0 for all vp in S0(Vol_Dielectric_Ele).
The corresponding Euler-Lagrange are:
* Div ( epsr[] Grad v) = 0 on Vol_Dielectric_Ele
* epsr[] n.Grad v = 0 on Sur_Neu_Ele
The Galerkin{} statement is a symbolic representation of this weak formulation term.
The first argument represents the density to be integrated over finite elements,
the result of which being the bilinear form, i.e., the element-based matrix
that will be assembled in the linear system of FE equations.
What is a bit confusing is that the two comma-separated arguments of the bracket []
are not exactly interpreted the same way.
The symbol "d" represents the exterior derivative and is a synonym of "Grad".
The expression "d v" stands then for the gradient of the v field,
i.e., for the electric field up to a minus sign.
The Dof{} operator, finally, returns the array of basis functions of its argument.
As the Galerkin method uses as trial (test) functions the basis functions
of the unknown field, the density of the bilinear form should be
[ epsr[] * Dof{d v} , Dof{d v} ]
to highlight the symmetry.
But the second argument of the density is always reserved for the trial functions
and the Dof{} operator should therefore always be there.
It can thus be made implicit and, according to the GetDP syntax, it must be omitted.
Hence the final expression of the density below.
*/
{ Name Electrostatics_v; Type FemEquation;
Quantity {
{ Name v; Type Local; NameOfSpace Hgrad_v_Ele; }
......@@ -115,7 +169,7 @@ Formulation {
Equation {
Galerkin { [ epsr[] * Dof{d v} , {d v} ];
In Vol_Dielectric_Ele;
Jacobian Vol; Integration GradGrad; }
Jacobian Vol; Integration Int; }
}
}
}
......@@ -134,6 +188,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 scalar vector potential
- the electric fied
- 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;
......
Magnetostatic/electromagnet.geo 0 → 100644
+ 54
− 0
View file @ 7b2c9bfe
/* -------------------------------------------------------------------
File "electromagnet.geo"
This file is the geometrical description used by GMSH to produce
the file "electromagnet.msh".
------------------------------------------------------------------- */
dxCore = 50.e-3; dyCore = 100.e-3;
xInd = 75.e-3; dxInd = 25.e-3; dyInd = 50.e-3;
Include "electromagnet_common.pro";
s = DefineNumber[1, Name "Model parameters/Global mesh size",
Help "Reduce for finer mesh"];
p0 = 12.e-3 *s;
pCorex = 4.e-3 *s; pCorey0 = 8.e-3 *s; pCorey = 4.e-3 *s;
pIndx = 5.e-3 *s; pIndy = 5.e-3 *s;
pInt = 12.5e-3*s; pExt = 12.5e-3*s;
Point(1) = {0,0,0,p0};
Point(2) = {dxCore,0,0,pCorex};
Point(3) = {dxCore,dyCore,0,pCorey};
Point(4) = {0,dyCore,0,pCorey0};
Point(5) = {xInd,0,0,pIndx};
Point(6) = {xInd+dxInd,0,0,pIndx};
Point(7) = {xInd+dxInd,dyInd,0,pIndy};
Point(8) = {xInd,dyInd,0,pIndy};
Point(9) = {rInt,0,0,pInt};
Point(10) = {rExt,0,0,pExt};
Point(11) = {0,rInt,0,pInt};
Point(12) = {0,rExt,0,pExt};
Line(1) = {1,2}; Line(2) = {2,5}; Line(3) = {5,6};
Line(4) = {6,9}; Line(5) = {9,10}; Line(6) = {1,4};
Line(7) = {4,11}; Line(8) = {11,12}; Line(9) = {2,3};
Line(10) = {3,4}; Line(11) = {6,7}; Line(12) = {7,8};
Line(13) = {8,5};
Circle(14) = {9,1,11}; Circle(15) = {10,1,12};
Line Loop(16) = {-6,1,9,10}; Plane Surface(17) = {16};
Line Loop(18) = {11,12,13,3}; Plane Surface(19) = {18};
Line Loop(20) = {7,-14,-4,11,12,13,-2,9,10}; Plane Surface(21) = {-20}; // "-" for orientation
Line Loop(22) = {8,-15,-5,14}; Plane Surface(23) = {-22};
Physical Surface(101) = {21}; /* Air */
Physical Surface(102) = {17}; /* Core */
Physical Surface(103) = {19}; /* Ind */
Physical Surface(111) = {23}; /* AirInf */
Physical Line(1100) = {1,2,3,4,5}; /* Surface ht = 0 */
Physical Line(1101) = {6,7,8}; /* Surface bn = 0 */
Physical Line(1102) = {15}; /* Surface Inf */
Magnetostatic/electromagnet.pro 0 → 100644
+ 254
− 0
View file @ 7b2c9bfe
/* -------------------------------------------------------------------
Tutorial 2 : magnetostatic field of an electromagnet
Features:
- Infinite ring geometrical transformation
- Parameters shared by Gmsh and GetDp, and Onelab parameters
- FunctionSpaces for the 2D vector potential formulation
To compute the solution in a terminal:
getdp electromagnet -solve MagSta_a
getdp electromagnet -pos Map
To compute the solution interactively from the Gmsh GUI:
File > Open > electromagnet.pro
Run (button at the bottom of the left panel)
------------------------------------------------------------------- */
/* Electromagnetic fields expand to infinity.
The corresponding boundary condition can be imposed rigorously
by means of a gometrical transformation that maps a ring (or shell) of finite elements
to the complementary of its interior.
As this is a mere geometric transformation,
it is enough in the model description to attribute a special jacobian
to the ring region ("AirInf"). See Jacobian{} section below.
With this information, GetDP is able to deal with the correct transformation
of all quantities involved in the model.
The special jacobian "VolSphShell" has parameters.
There are 2 parameters in this case, "Val_Rint" and "Val_Rext",
which represent the inner and outer radii of the transformed ring region
and whose value must match those used
in the geometrical description of the model (.geo file).
This is a typical case where Gmsh and GetDP must consistently share parameter values.
To ensure consistency in all cases, common parameters are defined
is a specific file "electromagnet_common.pro",
which is included in both the .geo and .pro file of the model.
Besides sharing parameters between Gmsh and GetDP,
it is also useful to share some parameters (not all) with the user of the model,
i.e., to make them editable in the GUI before running the model.
Such variables are called Onelab variables (because the sharing mechanism
between the model and the GUI uses the Onelab interface).
Onelab parameters are defined with a "DefineNumber" statement,
which can be invoked in the .geo, .pro, or _common.pro files.
*/
Group {
// Physical regions:
Air = Region[ 101 ]; Core = Region[ 102 ];
Ind = Region[ 103 ]; AirInf = Region[ 111 ];
Surface_ht0 = Region[ 1100 ];
Surface_bn0 = Region[ 1101 ];
Surface_Inf = Region[ 1102 ];
/* Abstract regions :
The purpose of abstract regions is to allow a generic definition of
the FunctionSpace, Formulation and PostProcessing fields
with no reference to model-specific Physical regions.
We will show in a later tutorial how abstract formulations can then be isolated
in geometry independent template files, thanks to an appropriate declaration mechanism
(using DefineConstant[], DefineGroup[] and DefineFunction[]).
The abstract regions in this model have the following interpretation:
- Vol_Nu_Mag = region where the term [ nu[] * Dof{d a} , {d a} ] is assembled
- Vol_Js_Mag = region where the term [ - Dof{js} , {a} ] is assembled
- Vol_Inf = region where the infinite ring geometric transformation is applied
- Sur_Dir_Mag = Homogeneous Dirichlet part of the model's boundary;
- Sur_Neu_Mag = Homogeneous Neumann part of the model's boundary;
*/
Vol_Nu_Mag = Region[ {Air, AirInf, Core, Ind} ];
Vol_Js_Mag = Region[ Ind ];
Vol_Inf = Region[ {AirInf} ];
Sur_Dir_Mag = Region[ {Surface_bn0, Surface_Inf} ];
Sur_Neu_Mag = Region[ {Surface_ht0} ];
}
Function {
mu0 = 4.e-7 * Pi;
murCore = DefineNumber[100, Name "Model parameters/Mur core",
Help "Magnetic relative permeability of Core"];
nu [ Region[{Air, Ind, AirInf}] ] = 1. / mu0;
nu [ Core ] = 1. / (murCore * mu0);
NbTurns = 1000 ;
Current = DefineNumber[0.01, Name "Model parameters/Current",
Help "Current injected in coil [A]"];
Js_fct[ Ind ] = NbTurns*Current/SurfaceArea[];
}
/* In the 2D approximation, the magnetic vector potential A and the current density Js
are not scalars, but vectors with a z-component only:
A = Vector [ 0, 0, az(x,y,t) ]
Js = Vector [ 0, 0, jsz(x,y,t) ]
In order to compute derivatives and apply geometric transformations adequately,
GetDP needs this information.
Regarding discretization, now, A is node-based, whereas Js is region-wise constant.
Considering all this, one ends then up with the FunctionSpaces
"Hcurl_a_Mag_2D" and "Hregion_j_Mag_2D" as they are defined below.
The function space "Hregion_j_Mag_2D" provides one basis function,
and hence one degree of freedom, per physical region in the abstract region "Vol_Js_Mag".
The constraint "SourceCurrentDensityZ" fixes all these dofs,
so the FunctionSpace "Hregion_j_Mag_2D" is fully fixed and has no FE unknowns.
One could thus have replaced it by a simple function
and the Galerkin term would have been
Galerkin { [ Vector[ 0,0,-Js_fct[] ] , {a} ]; In Vol_Js_Mag;
Jacobian Vol; Integration Int; }
instead of
Galerkin { [ - Dof{js} , {a} ]; In Vol_Js_Mag;
Jacobian Vol; Integration Int; }
Thechosen implementation below is however more effeicient
as it avoids evaluating repeatedly the function Js_fct[] during assembly.
*/
Constraint {
{ Name Dirichlet_a_Mag;
Case {
{ Region Sur_Dir_Mag ; Value 0.; }
}
}
{ Name SourceCurrentDensityZ;
Case {
{ Region Vol_Js_Mag ; Value Js_fct[]; }
}
}
}
Group {
Dom_Hcurl_a_Mag_2D = Region[ {Vol_Nu_Mag, Sur_Neu_Mag} ];
}
FunctionSpace {
{ Name Hcurl_a_Mag_2D; Type Form1P; // Magnetic vector potential A
BasisFunction {
{ Name se; NameOfCoef ae; Function BF_PerpendicularEdge;
Support Dom_Hcurl_a_Mag_2D ; Entity NodesOf[ All ]; }
}
Constraint {
{ NameOfCoef ae; EntityType NodesOf;
NameOfConstraint Dirichlet_a_Mag; }
}
}
{ Name Hregion_j_Mag_2D; Type Vector; // Electric current density Js
BasisFunction {
{ Name sr; NameOfCoef jsr; Function BF_RegionZ;
Support Vol_Js_Mag; Entity Vol_Js_Mag; }
}
Constraint {
{ NameOfCoef jsr; EntityType Region;
NameOfConstraint SourceCurrentDensityZ; }
}
}
}
Include "electromagnet_common.pro";
Val_Rint = rInt; Val_Rext = rExt;
Jacobian {
{ Name Vol ;
Case { { Region Vol_Inf ;
Jacobian VolSphShell {Val_Rint, Val_Rext} ; }
{ Region All ; Jacobian Vol ; }
}
}
}
Integration {
{ Name Int ;
Case { {Type Gauss ;
Case { { GeoElement Triangle ; NumberOfPoints 4 ; }
{ GeoElement Quadrangle ; NumberOfPoints 4 ; }
}
}
}
}
}
Formulation {
{ Name Magnetostatics_a_2D; Type FemEquation;
Quantity {
{ Name a ; Type Local; NameOfSpace Hcurl_a_Mag_2D; }
{ Name js; Type Local; NameOfSpace Hregion_j_Mag_2D; }
}
Equation {
Galerkin { [ nu[] * Dof{d a} , {d a} ]; In Vol_Nu_Mag;
Jacobian Vol; Integration Int; }
Galerkin { [ - Dof{js} , {a} ]; In Vol_Js_Mag;
Jacobian Vol; Integration Int; }
}
}
}
Resolution {
{ Name MagSta_a;
System {
{ Name Sys_Mag; NameOfFormulation Magnetostatics_a_2D; }
}
Operation {
Generate[Sys_Mag]; Solve[Sys_Mag]; SaveSolution[Sys_Mag];
}
}
}
PostProcessing {
{ Name MagSta_a_2D; NameOfFormulation Magnetostatics_a_2D;
Quantity {
{ Name a;
Value {
Local { [ {a} ]; In Dom_Hcurl_a_Mag_2D; Jacobian Vol; }
}
}
{ Name az;
Value {
Local { [ CompZ[{a}] ]; In Dom_Hcurl_a_Mag_2D; Jacobian Vol; }
}
}
{ Name b;
Value {
Local { [ {d a} ]; In Dom_Hcurl_a_Mag_2D; Jacobian Vol; }
}
}
{ Name h;
Value {
Local { [ nu[] * {d a} ]; In Dom_Hcurl_a_Mag_2D; Jacobian Vol; }
}
}
}
}
}
e = 1.e-5;
p1 = {e,e,0};
p2 = {0.25-e,e,0}; // horizontal cut through model, just above x-axis.
PostOperation {
{ Name Map_a; NameOfPostProcessing MagSta_a_2D;
Operation {
Echo[ Str["l=PostProcessing.NbViews-1; View[l].IntervalsType = 3;"],
File "tmp.geo", LastTimeStepOnly] ;
Print[ az, OnElementsOf Dom_Hcurl_a_Mag_2D, File "az.pos" ];
Print[ b, OnLine{{List[p1]}{List[p2]}} {1000}, File "by.pos" ];
}
}
}
Magnetostatic/electromagnet_common.pro 0 → 100644
+ 4
− 0
View file @ 7b2c9bfe
// Parameters shared by Gmsh and GetDP
rInt = 200.e-3;
rExt = 250.e-3;
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment