more cleanups

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
 

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; }

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 ;
@@ -193,6 +202,25 @@ Integration {
   }
 }
 
+/* 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; }
     }
   }
 }