From 58476bd78530b848acb62cf681a003afde9646db Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Thu, 15 Mar 2018 21:42:14 +0100
Subject: [PATCH] more cleanups
---
Electrostatics/microstrip.pro | 6 +-
ElectrostaticsFloating/floating.pro | 10 +--
Magnetostatics/electromagnet.pro | 110 +++++++++++++++++-----------
3 files changed, 77 insertions(+), 49 deletions(-)
diff --git a/Electrostatics/microstrip.pro b/Electrostatics/microstrip.pro
index 7b152d6..3ac1462 100644
--- a/Electrostatics/microstrip.pro
+++ b/Electrostatics/microstrip.pro
@@ -54,7 +54,7 @@ Group {
Vol_Ele : volume where -Div(epsilon Grad v) = 0 is solved
Sur_Neu_Ele: surface where non homogeneous Neumann boundary conditions
- (on n.d == epsilon n.Grad v) are imposed
+ (on n.d = -n . (epsilon 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).
@@ -181,11 +181,11 @@ Formulation {
product over a domain D. 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')_Bnd_Vol_Ele = 0
+ (epsilon Grad v, Grad v')_Vol_Ele + (n . (epsilon Grad v), v')_Bnd_Vol_Ele = 0
holds for all v', where Bnd_Vol_Ele is the boundary of Vol_Ele. In our
microstrip example this surface term vanishes, as there is either no test
- function v' (on the Dirichlet boundary), or "epsilon n.Grad v" is zero
+ function v' (on the Dirichlet boundary), or "n. (epsilon Grad v)" is zero
(on the homogeneous Neumann boundary). We are thus eventually looking for
functions v in the function space Hgrad_v_Ele such that
diff --git a/ElectrostaticsFloating/floating.pro b/ElectrostaticsFloating/floating.pro
index 5b668d5..58f6b19 100644
--- a/ElectrostaticsFloating/floating.pro
+++ b/ElectrostaticsFloating/floating.pro
@@ -58,7 +58,7 @@ Group {
Vol_Ele : volume where -Div(epsilon Grad v) = 0 is solved
Sur_Neu_Ele : surface where non homogeneous Neumann boundary conditions
- (on n.d == epsilon n.Grad v) are imposed
+ (on n.d = -n . (epsilon Grad v)) are imposed
Sur_Electrodes_Ele : electrode regions */
Vol_Ele = Region[ {Air, Diel1} ];
@@ -140,15 +140,15 @@ FunctionSpace {
carried by that electrode. Indeed, let us consider the electrostatic weak
formulation derived in Tutorial 1: find v in Hgradv_Ele such that
- (epsilon Grad v, Grad v')_Vol_Ele + (epsilon n.Grad v, v')_Bnd_Vol_Ele = 0
+ (epsilon Grad v, Grad v')_Vol_Ele + (n . (epsilon Grad v), v')_Bnd_Vol_Ele = 0
holds for all test functions v'. When the test-function v' is BF_electrode,
the boundary term reduces to
- (epsilon n.Grad v, BF_electrode)_Sur_Electrodes_Ele.
+ (n . (epsilon Grad v), BF_electrode)_Sur_Electrodes_Ele.
Since BF_electrode == 1 on the electrode, the boundary term is actually
- simply equal to the integral of (epsilon n.Grad v) on the electrode,
+ simply equal to the integral of (n . epsilon Grad v) on the electrode,
i.e. the flux of the displacement field, which is by definition the
charge Q_electrode carried by the electrodes.
@@ -215,7 +215,7 @@ Formulation {
satisfies (just consider the equation corresponding to the test function
BF_electrode):
- Q_electrode = (-epsilon[] Grad v, Grad BF_electrode)_Vol_Ele */
+ Q_electrode = (-epsilon Grad v, Grad BF_electrode)_Vol_Ele */
{ Name Electrostatics_v; Type FemEquation;
Quantity {
{ Name v; Type Local; NameOfSpace Hgrad_v_Ele; }
diff --git a/Magnetostatics/electromagnet.pro b/Magnetostatics/electromagnet.pro
index 37293ee..d1956f3 100644
--- a/Magnetostatics/electromagnet.pro
+++ b/Magnetostatics/electromagnet.pro
@@ -81,9 +81,30 @@ Group {
Vol_S_Mag = Region[ Ind ];
Vol_Inf_Mag = Region[ AirInf ];
Sur_Dir_Mag = Region[ {Surface_bn0, Surface_Inf} ];
- Sur_Neu_Mag = Region[ {Surface_ht0} ];
+ Sur_Neu_Mag = Region[ Surface_ht0 ];
}
+/* The weak formulation for this problem is derived in a similar way to the
+ electrostatic weak formulation from Tutorial 1. The main difference is that
+ the fields are now vector-valued, and that we have one linear term in
+ addition to the bilinear term. The weak formulation reads: find a such that
+
+ (Curl(nu Curl a), a')_Vol_Mag = (js, a')_Vol_S_Mag
+
+ holds for all test functions a'. After integration by parts it reads: find a
+ such that
+
+ (nu Curl a, Curl a')_Vol_Mag + (n x (nu Curl a), a')_Bnd_Vol_Mag = (js, a')_Vol_S_Mag
+
+ In this electromagnet model the second (boundary term) vanishes, as there is
+ either no test function a' (on the Dirichlet boundary), or "n x (nu Curl a) =
+ n x h" is zero (on the homogeneous Neumann boundary). We are thus eventually
+ looking for functions a such that
+
+ (nu Curl a, Curl a')_Vol_Mag = (js, a')_Vol_S_Mag
+
+ holds for all a'. */
+
Function {
mu0 = 4.e-7 * Pi;
murCore = DefineNumber[100, Name "Model parameters/Mur core",
@@ -92,44 +113,14 @@ Function {
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]"];
+ NbTurns = 1000 ; // number of turns in the coil
js_fct[ Ind ] = -NbTurns*Current/SurfaceArea[];
/* The minus sign is to have the current in -e_z direction,
so that the magnetic induction field is in +e_y direction */
}
-/* 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_S_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 Integral term would have been
-
- Integral { [ Vector[ 0,0,-js_fct[] ] , {a} ]; In Vol_S_Mag;
- Jacobian Vol; Integration Int; }
-
- instead of
-
- Integral { [ - Dof{js} , {a} ]; In Vol_S_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 {
@@ -143,11 +134,28 @@ Constraint {
}
}
+/* In the 2D approximation, the magnetic vector potential a and the current
+ density js are vectors with a z-component only, i.e.:
+
+ a = Vector [ 0, 0, az(x,y) ]
+ js = Vector [ 0, 0, jsz(x,y) ]
+
+ These vector fields behave differently under derivation and geometrical
+ transformation though, and GetDP needs this information to perform these
+ operations correctly. This is reflected in the Type, BasisFunction and Entity
+ specified in the "Hcurl_a_Mag_2D" FunctionSpace for a ("Perpendicular 1-form"
+ with "BF_PerpendicularEdge" basis functions associated to nodes of the mesh)
+ and in the "Hregion_j_Mag_2D" FunctionSpace for js ("Vector" with
+ "BF_RegionZ" basis functions associated with the region Vol_S_Mag). Without
+ giving all the details, a is thus node-based, whereas js is region-wise
+ constant. */
+
Group {
Dom_Hcurl_a_Mag_2D = Region[ {Vol_Mag, Sur_Neu_Mag} ];
}
+
FunctionSpace {
- { Name Hcurl_a_Mag_2D; Type Form1P; // Magnetic vector potential A
+ { 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 ]; }
@@ -158,7 +166,7 @@ FunctionSpace {
}
}
- { Name Hregion_j_Mag_2D; Type Vector; // Electric current density Js
+ { Name Hregion_j_Mag_2D; Type Vector; // Electric current density js
BasisFunction {
{ Name sr; NameOfCoef jsr; Function BF_RegionZ;
Support Vol_S_Mag; Entity Vol_S_Mag; }
@@ -173,6 +181,7 @@ FunctionSpace {
Include "electromagnet_common.pro";
Val_Rint = rInt; Val_Rext = rExt;
+
Jacobian {
{ Name Vol ;
Case { { Region Vol_Inf_Mag ;
@@ -184,15 +193,34 @@ Jacobian {
Integration {
{ Name Int ;
- Case { {Type Gauss ;
- Case { { GeoElement Triangle ; NumberOfPoints 4 ; }
- { GeoElement Quadrangle ; NumberOfPoints 4 ; }
+ Case { { Type Gauss ;
+ Case { { GeoElement Triangle ; NumberOfPoints 4 ; }
+ { GeoElement Quadrangle ; NumberOfPoints 4 ; }
}
}
}
}
}
+/* 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_S_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 Integral term
+ below in the weak formulation could have been written as
+
+ Integral { [ Vector[ 0,0,-js_fct[] ] , {a} ];
+ In Vol_S_Mag; Jacobian Vol; Integration Int; }
+
+ instead of
+
+ Integral { [ - Dof{js} , {a} ];
+ In Vol_S_Mag; Jacobian Vol; Integration Int; }
+
+ The chosen implementation below is however more effeicient as it avoids
+ evaluating repeatedly the function js_fct[] during assembly.
+*/
+
Formulation {
{ Name Magnetostatics_a_2D; Type FemEquation;
Quantity {
@@ -200,10 +228,10 @@ Formulation {
{ Name js; Type Local; NameOfSpace Hregion_j_Mag_2D; }
}
Equation {
- Integral { [ nu[] * Dof{d a} , {d a} ]; In Vol_Mag;
- Jacobian Vol; Integration Int; }
- Integral { [ -Dof{js} , {a} ]; In Vol_S_Mag;
- Jacobian Vol; Integration Int; }
+ Integral { [ nu[] * Dof{d a} , {d a} ];
+ In Vol_Mag; Jacobian Vol; Integration Int; }
+ Integral { [ -Dof{js} , {a} ];
+ In Vol_S_Mag; Jacobian Vol; Integration Int; }
}
}
}
--
GitLab