tutorial 2

Magnetostatic/electromagnet.geo

+/* -------------------------------------------------------------------
+   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

+/* -------------------------------------------------------------------
+   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

+// Parameters shared by Gmsh and GetDP
+
+rInt   = 200.e-3; 
+rExt   = 250.e-3;