new tuto 9 on magnetic forces

MagneticForces/InfiniteBox.geo

+Function Cuboid
+// This function receives a vector pnt[] with 8 points
+// and generates 12 lines, 6 line loops, 1 volume,
+// 1 physical volume and one physicale surface
+
+l12 = newl; Line(l12) = {pnt[0], pnt[1]};
+l23 = newl; Line(l23) = {pnt[1], pnt[2]};
+l34 = newl; Line(l34) = {pnt[2], pnt[3]};
+l41 = newl; Line(l41) = {pnt[3], pnt[0]};
+l56 = newl; Line(l56) = {pnt[4], pnt[5]};
+l67 = newl; Line(l67) = {pnt[5], pnt[6]};
+l78 = newl; Line(l78) = {pnt[6], pnt[7]};
+l85 = newl; Line(l85) = {pnt[7], pnt[4]};
+l15 = newl; Line(l15) = {pnt[0], pnt[4]};
+l26 = newl; Line(l26) = {pnt[1], pnt[5]};
+l37 = newl; Line(l37) = {pnt[2], pnt[6]};
+l48 = newl; Line(l48) = {pnt[3], pnt[7]};
+
+ll1 = newll; Line Loop(ll1) = { l12, l23, l34, l41 };
+ll2 = newll; Line Loop(ll2) = { l56, l67, l78, l85 };
+ll3 = newll; Line Loop(ll3) = { l12, l26, -l56, -l15 };
+ll4 = newll; Line Loop(ll4) = { l23, l37, -l67, -l26 };
+ll5 = newll; Line Loop(ll5) = { l34, l48, -l78, -l37 };
+ll6 = newll; Line Loop(ll6) = { l41, l15, -l85, -l48 };
+
+s1 = news; Plane Surface(s1) = { ll1 } ;
+s2 = news; Plane Surface(s2) = { ll2 } ;
+s3 = news; Plane Surface(s3) = { ll3 } ;
+s4 = news; Plane Surface(s4) = { ll4 } ;
+s5 = news; Plane Surface(s5) = { ll5 } ;
+s6 = news; Plane Surface(s6) = { ll6 } ;
+
+sl = newsl; Surface Loop(sl) = { -s1, s2, s3, s4, s5, s6 };
+v = newv; Volume(v) = { sl };
+
+If( Flag_TransfInf )
+//Mesh.Algorithm3D = 4; // (4=Frontal, 5=Frontal Delaunay, 6=Frontal Hex, 7=MMG3D, 9=R-tree)
+  For num In { l12:l85 } 
+    Transfinite Line{ num } = 10;
+  EndFor
+  For num In { l15:l48 } 
+    Transfinite Line{ num } = 5;
+  EndFor
+  For num In { s1:s6 } 
+    Transfinite Surface{ num } ;
+  EndFor
+  Transfinite Volume{ v } ;
+EndIf
+Return
+
+
+DefineConstant[
+  Flag_InfiniteBox = {1, Choices{0,1}, Name "Infinite box/Add infinite box"}
+  Flag_TransfInf = {0, Choices{0,1}, Name "Infinite box/Transfinite mesh", Visible 0}
+  ratioInf = {1.5, Name "Infinite box/Ratio ext-int", Visible Flag_InfiniteBox}
+  //ratioLc = {25, Name "Infinite box/ration int-Lc", Visible Flag_InfiniteBox}
+  xInt = {1, Name "Infinite box/xInt", Visible 0}
+  yInt = {1, Name "Infinite box/yInt", Visible 0}
+  zInt = {1, Name "Infinite box/zInt", Visible 0}
+  xExt = {xInt*ratioInf, Name "Infinite box/xExt", Visible 0}
+  yExt = {yInt*ratioInf, Name "Infinite box/yExt", Visible 0}
+  zExt = {zInt*ratioInf, Name "Infinite box/zExt", Visible 0}
+  xCnt = {0, Name "Infinite box/xCenter", Visible 0}
+  yCnt = {0, Name "Infinite box/yCenter", Visible 0}
+  zCnt = {0, Name "Infinite box/zCenter", Visible 0}
+];
+
+If(!Flag_InfiniteBox)
+  Abort; // exit
+EndIf
+
+// Compute parameters related to the Infinite box
+
+BoundingBox;
+xCnt = (General.MinX + General.MaxX) / 2.;
+yCnt = (General.MinY + General.MaxY) / 2.;
+zCnt = (General.MinZ + General.MaxZ) / 2.;
+xInt = (General.MaxX - General.MinX)/2. * ratioBox;
+yInt = (General.MaxY - General.MinZ)/2. * ratioBox;
+zInt = (General.MaxZ - General.MinZ)/2. * ratioBox;
+xExt = xInt * ratioInf;
+yExt = yInt * ratioInf;
+zExt = zInt * ratioInf;
+SetNumber("Infinite box/xInt",xInt);
+SetNumber("Infinite box/yInt",yInt);
+SetNumber("Infinite box/zInt",zInt);
+SetNumber("Infinite box/xExt",xExt);
+SetNumber("Infinite box/yExt",yExt);
+SetNumber("Infinite box/zExt",zExt);
+SetNumber("Infinite box/xCenter",xCnt);
+SetNumber("Infinite box/yCenter",yCnt);
+SetNumber("Infinite box/zCenter",zCnt);
+
+lc1inf = lc2;
+lc2inf = lc1inf;
+
+p1 = newp; Point (p1) = {xCnt-xInt, yCnt-yInt, zCnt-zInt, lc1inf};
+p2 = newp; Point (p2) = {xCnt+xInt, yCnt-yInt, zCnt-zInt, lc1inf};
+p3 = newp; Point (p3) = {xCnt+xInt, yCnt+yInt, zCnt-zInt, lc1inf};
+p4 = newp; Point (p4) = {xCnt-xInt, yCnt+yInt, zCnt-zInt, lc1inf};
+p5 = newp; Point (p5) = {xCnt-xInt, yCnt-yInt, zCnt+zInt, lc1inf};
+p6 = newp; Point (p6) = {xCnt+xInt, yCnt-yInt, zCnt+zInt, lc1inf};
+p7 = newp; Point (p7) = {xCnt+xInt, yCnt+yInt, zCnt+zInt, lc1inf};
+p8 = newp; Point (p8) = {xCnt-xInt, yCnt+yInt, zCnt+zInt, lc1inf};
+
+pp1 = newp; Point (pp1) = {xCnt-xExt, yCnt-yExt, zCnt-zExt, lc2inf};
+pp2 = newp; Point (pp2) = {xCnt+xExt, yCnt-yExt, zCnt-zExt, lc2inf};
+pp3 = newp; Point (pp3) = {xCnt+xExt, yCnt+yExt, zCnt-zExt, lc2inf};
+pp4 = newp; Point (pp4) = {xCnt-xExt, yCnt+yExt, zCnt-zExt, lc2inf};
+pp5 = newp; Point (pp5) = {xCnt-xExt, yCnt-yExt, zCnt+zExt, lc2inf};
+pp6 = newp; Point (pp6) = {xCnt+xExt, yCnt-yExt, zCnt+zExt, lc2inf};
+pp7 = newp; Point (pp7) = {xCnt+xExt, yCnt+yExt, zCnt+zExt, lc2inf};
+pp8 = newp; Point (pp8) = {xCnt-xExt, yCnt+yExt, zCnt+zExt, lc2inf};
+
+SurfInf[] = {} ;
+SurfExt[] = {} ;
+VolInf[] = {};
+
+pnt[]={p1,p2,p3,p4,pp1,pp2,pp3,pp4};
+Call Cuboid; 
+SurfInt[] += {s1} ; SurfExt[] += {s2} ; VolInf[] += {v};
+
+pnt[]={p5,p6,p7,p8,pp5,pp6,pp7,pp8};
+Call Cuboid;
+SurfInt[] += {s1} ; SurfExt[] += {s2} ; VolInf[] += {v};
+
+pnt[]={p1,p2,p6,p5,pp1,pp2,pp6,pp5};
+Call Cuboid;
+SurfInt[] += {s1} ; SurfExt[] += {s2} ; VolInf[] += {v};
+
+pnt[]={p3,p4,p8,p7,pp3,pp4,pp8,pp7};
+Call Cuboid;
+SurfInt[] += {s1} ; SurfExt[] += {s2} ; VolInf[] += {v};
+
+pnt[]={p2,p3,p7,p6,pp2,pp3,pp7,pp6};
+Call Cuboid;
+SurfInt[] += {s1} ; SurfExt[] += {s2} ; VolInf[] += {v};
+
+pnt[]={p4,p1,p5,p8,pp4,pp1,pp5,pp8};
+Call Cuboid;
+SurfInt[] += {s1} ; SurfExt[] += {s2} ; VolInf[] += {v};
+
+Coherence;
+
+Physical Volume("InfiniteX", 1) = { VolInf[4], VolInf[5] };
+Physical Volume("InfiniteY", 2) = { VolInf[2], VolInf[3] };
+Physical Volume("InfiniteZ", 3) = { VolInf[0], VolInf[1] };
+
+

MagneticForces/magnets.geo

+// include parameters common to geometry and solver
+
+Include "magnets_common.pro";
+
+// define geometry-specific parameters
+DefineConstant[
+  lc1 = {TotalMemory <= 2048 ? 10*mm : 5*mm, Name "Parameters/2Mesh size on magnets [m]"}
+  lc2 = {TotalMemory <= 2048 ? 20*mm : 10*mm, Name "Parameters/2Mesh size at infinity [m]"}
+  ratioBox = {2.0, Name "Parameters/2Ratio Box to model size"}
+];
+
+// change global Gmsh options
+Mesh.Optimize = 1; // optimize quality of tetrahedra
+Mesh.VolumeEdges = 0; // hide volume edges
+Geometry.ExactExtrusion = 0; // to allow rotation of extruded shapes
+Solver.AutoMesh = 2; // always remesh if necessary (don't reuse mesh on disk)
+
+// create magnets
+For i In {1:NumMagnets}
+  If(M~{i} == 0) // cylinder
+    p1 = newp; Point(p1) = {X~{i}, Y~{i}-L~{i}/2, Z~{i}, lc1};
+    p2 = newp; Point(p2) = {X~{i}+R~{i}, Y~{i}-L~{i}/2, Z~{i}, lc1};
+    p3 = newp; Point(p3) = {X~{i}, Y~{i}+L~{i}/2, Z~{i}, lc1};
+    p4 = newp; Point(p4) = {X~{i}+R~{i}, Y~{i}+L~{i}/2, Z~{i}, lc1};
+    l1 = newl; Line(l1) = {p1,p2}; l2 = newl; Line(l2) = {p2,p4};
+    l3 = newl; Line(l3) = {p4,p3}; l4 = newl; Line(l4) = {p3,p1};
+    ll1 = newll; Line Loop(ll1) = {l1,l2,l3,l4};
+    s1 = news; Plane Surface(s1) = {ll1};
+    e1[] = Extrude {{0, 1, 0}, {X~{i}, Y~{i}, Z~{i}}, Pi/2} { Surface{s1}; };
+    e2[] = Extrude {{0, 1, 0}, {X~{i}, Y~{i}, Z~{i}}, Pi/2} { Surface{e1[0]}; };
+    e3[] = Extrude {{0, 1, 0}, {X~{i}, Y~{i}, Z~{i}}, Pi/2} { Surface{e2[0]}; };
+    e4[] = Extrude {{0, 1, 0}, {X~{i}, Y~{i}, Z~{i}}, Pi/2} { Surface{e3[0]}; };
+    Magnet~{i}[] = {e1[1], e2[1], e3[1], e4[1]};
+  EndIf
+  If(M~{i} == 1) // parallelepiped
+    p1 = newp; Point(p1) = {X~{i}-Lx~{i}/2, Y~{i}-Ly~{i}/2, Z~{i}-Lz~{i}/2, lc1};
+    p2 = newp; Point(p2) = {X~{i}+Lx~{i}/2, Y~{i}-Ly~{i}/2, Z~{i}-Lz~{i}/2, lc1};
+    p3 = newp; Point(p3) = {X~{i}+Lx~{i}/2, Y~{i}+Ly~{i}/2, Z~{i}-Lz~{i}/2, lc1};
+    p4 = newp; Point(p4) = {X~{i}-Lx~{i}/2, Y~{i}+Ly~{i}/2, Z~{i}-Lz~{i}/2, lc1};
+    l1 = newl; Line(l1) = {p1,p2}; l2 = newl; Line(l2) = {p2,p3};
+    l3 = newl; Line(l3) = {p3,p4}; l4 = newl; Line(l4) = {p4,p1};
+    ll1 = newll; Line Loop(ll1) = {l1,l2,l3,l4};
+    s1 = news; Plane Surface(s1) = {ll1};
+    e1[] = Extrude {0, 0, Lz~{i}} { Surface{s1}; };
+    Magnet~{i}[] = {e1[1]};
+  EndIf
+  Rotate { {0,0,1}, {X~{i},Y~{i},Z~{i}}, deg*Rz~{i} } { Volume{Magnet~{i}[]}; }
+  Rotate { {0,1,0}, {X~{i},Y~{i},Z~{i}}, deg*Ry~{i} } { Volume{Magnet~{i}[]}; }
+  Rotate { {1,0,0}, {X~{i},Y~{i},Z~{i}}, deg*Rx~{i} } { Volume{Magnet~{i}[]}; }
+
+  Physical Volume(Sprintf("Magnet_%g",i),10*i) = Magnet~{i}[]; // magnet volume
+  skin~{i}[] = CombinedBoundary{ Volume{Magnet~{i}[]}; };
+  Physical Surface(Sprintf("SkinMagnet_%g",i),10*i+1) = -skin~{i}[]; // magnet skin
+EndFor
+
+
+If (Flag_InfiniteBox )
+   Include "InfiniteBox.geo";
+Else
+// create air box around magnets
+  BoundingBox; // recompute model bounding box
+  cx = (General.MinX + General.MaxX) / 2;
+  cy = (General.MinY + General.MaxY) / 2;
+  cz = (General.MinZ + General.MaxZ) / 2;
+  lx = (General.MaxX - General.MinX) * 2 * ratioBox;
+  ly = (General.MaxY - General.MinZ) * ratioBox; 
+  lz = (General.MaxZ - General.MinZ) * 2 * ratioBox;
+  p1 = newp; Point (p1) = {cx-lx/2, cy-ly/2, cz-lz/2, lc2};
+  p2 = newp; Point (p2) = {cx+lx/2, cy-ly/2, cz-lz/2, lc2};
+  l1 = newl; Line(l1) = {p1, p2};
+  e1[] = Extrude {0, ly, 0} { Line{l1}; };
+  e2[] = Extrude {0, 0, lz} { Surface{e1[1]}; };
+  Delete { Volume{e2[1]}; }
+  SurfInt[] = { e1[1],e2[0],e2[2],e2[3],e2[4],e2[5] };
+  SurfExt[] = SurfInt[];
+EndIf
+
+
+
+For num In {0:#SurfInt[]-1}
+  Printf("%g %g", num, SurfInt[num]);
+EndFor
+
+temp = newsl; Surface Loop(temp) = { SurfInt[] };
+vv[] = { temp};
+For i In {1:NumMagnets}
+  temp = newsl; Surface Loop(temp) = { skin~{i}[] };
+  vv[] += temp;
+EndFor
+v1 = newv; Volume(v1) = { vv[] };
+Physical Volume("AirBox", 4) = v1;
+
+Physical Surface("OuterSurface", 5) = { SurfExt[] };
+Physical Surface("InnerSurface", 6) = { SurfInt[] };

MagneticForces/magnets.pro

+/* -------------------------------------------------------------------
+   Tutorial 9 : 3D magnetostatic dual formulations and magnetic forces
+
+   Features:
+   - Dual 3D magnetostatic formulations
+   - Boundary condition at infinity with infinite elements
+   - Maxwell stress tensor and rigid-body magnetic forces
+
+   To compute the solution in a terminal:
+       getdp magnets.pro
+
+   To compute the solution interactively from the Gmsh GUI:
+       File > Open > magnets.pro
+       Run (button at the bottom of the left panel)
+   ------------------------------------------------------------------- */
+
+/*
+ This rather didactic tutorial solves the electromagnetic field 
+ and the rigid-body forces acting on a set of magnetic pieces
+ of parallelepipedic or cylindrical shape.
+ Besides position and dimension, each piece is attributed 
+ a (constant) magnetic permability and/or a remanence field.
+ The pieces are all (a bit abusively) generically called "Magnet" 
+ in the problem decription below, irresective of whether they are
+ truly permanent magnets or ferromagnetic barrels. 
+
+ The tutorial model proposes the both dual 3D magnetostatic formulations:
+ the magnetic vector potential formulation with spanning-tree gauging,
+ and the scalar magnetic potential formulation.
+ The later is rather simple in this case since, as there are no conductors,
+ the known coercive field hc[] is the only source field "hs" one needs 
+ to represens the magnetic field:
+  
+   h = hs - grad phi   ,  hs = -hc.
+
+ As in tutorial 2 (magnetostatic field of an electromagnet), a shell
+ of so-called infinite elements is used here to impose the exact
+ zero-field boundary condition at infinity.
+ The shell is generated automatically by including "InfiniteBox.geo"
+ at the end of the geometrical description of the model. 
+ It can be placed rather close of the magnets without loss of accuracy.
+
+ The preferred way to compute electromagnetic forces in GetDP
+ is as an explicit by-product of the Maxwell stress tensor "TM[{b}]",
+ which is a material dependent function of the magnetic induction "b" field. 
+ Exactly like we computed the heat flux "q(S)" through a surface "S"
+ using a special auxiliary function "g(S)" associated with that surface 
+ in the tutorial "Tutorial 5 : thermal problem with contact resistances",
+ the magnetic force acting on a rigid body in empty space can be evaluated
+ as the flux of the Maxwell stress tensor through that surface.
+ There is one auxiliary function "g(SkinMagnet~{i}) = un~{i}"
+ for each magnet and the resultant magnetic force acting on "Magnet~{i}"
+ is given by the integral:
+
+ f~{i} = Integral [ TM[{b}] * {-grad un~{i}} ] ;
+
+ It should be insisted on the fact that the Maxwell stress is discontinuous
+ on material discontinuities, and that magnetic forces on rigid bodies
+ only depend on the Maxwell stress tensor of empty space and 
+ on the "b" and "h" field distribution on the outer side of "SkinMagnet~{i}".
+
+ "{-grad un~{i}}" in the above formula can be regarded 
+ as the normal vector to "SkinMagnet~{i}"
+ in the one element thick layer "layer~{i}" of finite elements 
+ around "Magnet~{i}", and "f~{i}", is thus indeed the flux of "TM[]"
+ through the surface of "Magnet~{i}".
+ 
+ The support of "{-grad un~{i}}" is limited to "layer~{i}",
+ which is much smaller than "AirBox".
+ To speed up the computation of forces, a special domain "Vol_Force"
+ for force integrations is defined, which contains only
+ the layers  "layer~{i}" of alla magnets.  
+*/
+
+Include "magnets_common.pro"
+
+DefineConstant[
+  // preset all getdp options and make them invisible
+  R_ = {"MagSta_a", Name "GetDP/1ResolutionChoices", Visible 1,
+	Choices {"MagSta_a", "MagSta_phi"}},
+  C_ = {"-solve -v 5 -v2 -bin", Name "GetDP/9ComputeCommand", Visible 0}
+  P_ = {"", Name "GetDP/2PostOperationChoices", Visible 0}
+];
+
+Group{
+  // Geometrical regions (give litteral labels to geometrical region numbers)
+  domInfX = Region[1];
+  domInfY = Region[2];
+  domInfZ = Region[3];
+  AirBox  = Region[4];
+  Outer   = Region[5];
+
+  For i In {1:NumMagnets}
+    Magnet~{i} = Region[ {(10*i)}]; 
+    SkinMagnet~{i} = Region[ {(10*i+1)} ];
+    Layer~{i} =  Region[AirBox, OnOneSideOf SkinMagnet~{i}] ;
+  EndFor
+
+  // Abstract Groups (group geometrical regions into formulation relevant groups)
+  Vol_Inf = Region[ {domInfX, domInfY, domInfZ} ]; 
+  Vol_Air = Region[ {AirBox, Vol_Inf} ];
+
+  Vol_Magnet = Region[{}];
+  Sur_Magnet = Region[{}];
+  Vol_Force = Region[{}];
+  For i In {1:NumMagnets}
+    Sur_Magnet += Region[SkinMagnet~{i}]; 
+    Vol_Magnet += Region[Magnet~{i}];
+    Vol_Layer += Region[Layer~{i}];
+  EndFor
+
+  Vol_mu = Region[ {Vol_Air, Vol_Magnet}];
+
+  Sur_Dirichlet_phi = Region[ Outer ];
+  Sur_Dirichlet_a   = Region[ Outer ];
+
+  Dom_Hgrad_phi = Region[ {Vol_Air, Vol_Magnet, Sur_Dirichlet_phi} ];
+  Dom_Hcurl_a = Region[ {Vol_Air, Vol_Magnet, Sur_Dirichlet_a} ];
+  Vol_Force = Region [ Vol_Layer ];
+  //Vol_Force = Region [ Vol_Air ];
+}
+
+Function{
+  mu0 = 4*Pi*1e-7;
+  mu[ Vol_Air ] = mu0;
+
+  For i In {1:NumMagnets}
+    // coercive field of magnets
+    DefineConstant[
+      HC~{i} = {-BR~{i}/mu0,
+        Name Sprintf("Parameters/Magnet %g/0Coercive magnetic field [Am^-1]", i), Visible 0}
+    ];
+    hc[Magnet~{i}] = Rotate[Vector[0, HC~{i}, 0], Rx~{i}, Ry~{i}, Rz~{i}];
+    br[Magnet~{i}] = Rotate[Vector[0, BR~{i}, 0], Rx~{i}, Ry~{i}, Rz~{i}];
+    mu[Magnet~{i}] = mu0*MUR~{i};
+  EndFor
+
+  nu[] = 1.0/mu[];
+
+  // Maxwell stress tensor (valid for both formulations and linear materials
+  TM[] = ( SquDyadicProduct[$1] - SquNorm[$1] * TensorDiag[0.5, 0.5, 0.5] ) / mu[] ;
+}
+
+Jacobian {
+  { Name Vol ;
+    Case {
+      { Region All ; Jacobian Vol ; }
+      {Region domInfX; Jacobian VolRectShell {xInt,xExt,1,xCnt,yCnt,zCnt};}      
+      {Region domInfY; Jacobian VolRectShell {yInt,yExt,2,xCnt,yCnt,zCnt};}
+      {Region domInfZ; Jacobian VolRectShell {zInt,zExt,3,xCnt,yCnt,zCnt};}
+    }
+  }
+}
+
+Integration {
+  { Name Int ; 
+    Case {
+      { Type Gauss ;
+        Case {
+	  { GeoElement Triangle    ; NumberOfPoints 4 ; }
+	  { GeoElement Quadrangle  ; NumberOfPoints 4 ; }
+          { GeoElement Tetrahedron ; NumberOfPoints 4 ; }
+	  { GeoElement Hexahedron  ; NumberOfPoints  6 ; }
+	  { GeoElement Prism       ; NumberOfPoints  6 ; } 
+	}
+      }
+    }
+  }
+}
+
+Constraint {
+  { Name phi ;
+    Case {
+      { Region Sur_Dirichlet_phi ; Value 0. ; }
+    }
+  }
+  { Name a ;
+    Case {
+      { Region Sur_Dirichlet_a ; Value 0. ; }
+    }
+  }
+  { Name GaugeCondition_a ; Type Assign ;
+    Case {
+      { Region Dom_Hcurl_a ; SubRegion Sur_Dirichlet_a ; Value 0. ; }
+    }
+  }
+  For i In {1:NumMagnets}
+    { Name Magnet~{i} ;
+      Case {
+        { Region SkinMagnet~{i} ; Value 1. ; }
+      }
+    }
+  EndFor
+}
+
+FunctionSpace {
+  { Name Hgrad_phi ; Type Form0 ; // scalar magnetic potential
+    BasisFunction {
+      { Name sn ; NameOfCoef phin ; Function BF_Node ;
+        Support Dom_Hgrad_phi ; Entity NodesOf[ All ] ; }
+    }
+    Constraint {
+      { NameOfCoef phin ; EntityType NodesOf ; NameOfConstraint phi ; }
+    }
+  }
+  { Name Hcurl_a; Type Form1; // vector magnetic potential
+    BasisFunction {
+      { Name se;  NameOfCoef ae;  Function BF_Edge; 
+	Support Dom_Hcurl_a ;Entity EdgesOf[ All ]; }
+    }
+    Constraint {
+      { NameOfCoef ae;  EntityType EdgesOf ; NameOfConstraint a; }
+      { NameOfCoef ae;  EntityType EdgesOfTreeIn ; EntitySubType StartingOn ;
+        NameOfConstraint GaugeCondition_a ; }
+    }
+  }
+  // auxiliary field on layer of elements touching each magnet, for the
+  // accurate integration of the Maxwell stress tensor (using the gradient of
+  // this field)
+  For i In {1:NumMagnets}
+    { Name Magnet~{i} ; Type Form0 ;
+      BasisFunction {
+        { Name sn ; NameOfCoef un ; Function BF_GroupOfNodes ;
+          Support Vol_Air ; Entity GroupsOfNodesOf[ SkinMagnet~{i} ] ; }
+      }
+      Constraint {
+        { NameOfCoef un ; EntityType GroupsOfNodesOf ; NameOfConstraint Magnet~{i} ; }
+      }
+    }
+  EndFor
+}
+
+Formulation {
+  { Name MagSta_phi ; Type FemEquation ;
+    Quantity {
+      { Name phi ; Type Local ; NameOfSpace Hgrad_phi ; }
+      For i In {1:NumMagnets}
+        { Name un~{i} ; Type Local ; NameOfSpace Magnet~{i} ; }
+      EndFor
+    }
+    Equation {
+      Galerkin { [ - mu[] * Dof{d phi} , {d phi} ] ;
+        In Vol_mu ; Jacobian Vol ; Integration Int ; }
+      Galerkin { [ - mu[] * hc[] , {d phi} ] ;
+        In Vol_Magnet ; Jacobian Vol ; Integration Int ; }
+      For i In {1:NumMagnets} // dummy term to define dofs for fully fixed space
+        Galerkin { [ 0 * Dof{un~{i}} , {un~{i}} ] ;
+          In Vol_Air ; Jacobian Vol ; Integration Int ; }
+      EndFor
+    }
+  }
+  { Name MagSta_a; Type FemEquation ;
+    Quantity {
+      { Name a  ; Type Local  ; NameOfSpace Hcurl_a ; }
+      For i In {1:NumMagnets}
+        { Name un~{i} ; Type Local ; NameOfSpace Magnet~{i} ; }
+      EndFor
+    }
+    Equation {
+      Galerkin { [ nu[] * Dof{d a} , {d a} ] ;
+        In Vol_mu ; Jacobian Vol ; Integration Int ; }
+      Galerkin { [ nu[] * br[] , {d a} ] ;
+        In Vol_Magnet ; Jacobian Vol ; Integration Int ; }
+      For i In {1:NumMagnets} 
+      // dummy term to define dofs for fully fixed space
+        Galerkin { [ 0 * Dof{un~{i}} , {un~{i}} ] ;
+          In Vol_Air ; Jacobian Vol ; Integration Int ; }
+      EndFor
+    }
+  }
+}
+
+Resolution {
+  { Name MagSta_phi ;
+    System {
+      { Name A ; NameOfFormulation MagSta_phi ; }
+    }
+    Operation {
+      Generate[A] ; Solve[A] ; SaveSolution[A] ;
+      PostOperation[MagSta_phi] ;
+    }
+  }
+  { Name MagSta_a ;
+    System {
+      { Name A ; NameOfFormulation MagSta_a ; }
+    }
+    Operation {
+      Generate[A] ; Solve[A] ; SaveSolution[A] ;
+      PostOperation[MagSta_a] ;
+    }
+  }
+}
+
+PostProcessing {
+  { Name MagSta_phi ; NameOfFormulation MagSta_phi ;
+    Quantity {
+      { Name b   ; 
+	Value { Local { [ - mu[] * {d phi} ] ; In Dom_Hgrad_phi ; Jacobian Vol ; }
+	        Local { [ - mu[] * hc[] ]    ; In Vol_Magnet ; Jacobian Vol ; } } }
+      { Name h   ; 
+	Value { Local { [ - {d phi} ]        ; In Dom_Hgrad_phi ; Jacobian Vol ; } } }
+      { Name hc  ; Value { Local { [ hc[] ]  ; In Vol_Magnet ; Jacobian Vol ; } } }
+      { Name phi ; Value { Local { [ {phi} ] ; In Dom_Hgrad_phi ; Jacobian Vol ; } } }
+      For i In {1:NumMagnets}
+        { Name un~{i} ; Value { Local { [ {un~{i}} ] ; In Vol_Force ; Jacobian Vol ; } } }
+        { Name f~{i} ; Value { Integral { [ - TM[-mu[] * {d phi}] * {d un~{i}} ] ;
+              In Vol_Force ; Jacobian Vol ; Integration Int ; } } }
+        { Name fx~{i} ; Value { Integral { [ CompX[- TM[-mu[] * {d phi}] * {d un~{i}} ] ] ;
+              In Vol_Force ; Jacobian Vol ; Integration Int ; } } }
+        { Name fy~{i} ; Value { Integral { [ CompY[- TM[-mu[] * {d phi}] * {d un~{i}} ] ] ;
+              In Vol_Force ; Jacobian Vol ; Integration Int ; } } }
+        { Name fz~{i} ; Value { Integral { [ CompZ[- TM[-mu[] * {d phi}] * {d un~{i}} ] ] ;
+              In Vol_Force ; Jacobian Vol ; Integration Int ; } } }
+      EndFor
+    }
+  }
+  { Name MagSta_a ; NameOfFormulation MagSta_a ;
+    PostQuantity {
+      { Name b ; Value { Local { [ {d a} ]; In Dom_Hcurl_a ; Jacobian Vol; } } }
+      { Name a ; Value { Local { [ {a} ]; In Dom_Hcurl_a ; Jacobian Vol; } } }
+      { Name br ; Value { Local { [ br[] ]; In Vol_Magnet ; Jacobian Vol; } } }
+      For i In {1:NumMagnets}
+        { Name un~{i} ; Value { Local { [ {un~{i}} ] ; In Dom_Hcurl_a ; Jacobian Vol ; } } }
+        { Name f~{i} ; Value { Integral { [ - TM[{d a}] * {d un~{i}} ] ;
+              In Vol_Force ; Jacobian Vol ; Integration Int ; } } }
+        { Name fx~{i} ; Value { Integral { [ CompX[- TM[{d a}] * {d un~{i}} ] ] ;
+              In Vol_Force ; Jacobian Vol ; Integration Int ; } } }
+        { Name fy~{i} ; Value { Integral { [ CompY[- TM[{d a}] * {d un~{i}} ] ] ;
+              In Vol_Force ; Jacobian Vol ; Integration Int ; } } }
+        { Name fz~{i} ; Value { Integral { [ CompZ[- TM[{d a}] * {d un~{i}} ] ] ;
+              In Vol_Force ; Jacobian Vol ; Integration Int ; } } }
+      EndFor
+    }
+  }
+}
+
+PostOperation {
+  { Name MagSta_phi ; NameOfPostProcessing MagSta_phi;
+    Operation {
+      Print[ b, OnElementsOf Vol_mu, File "b.pos" ] ;
+      Echo[ Str["l=PostProcessing.NbViews-1;",
+		"View[l].ArrowSizeMax = 100;",
+		"View[l].CenterGlyphs = 1;",
+		"View[l].VectorType = 1;" ] ,
+        File "tmp.geo", LastTimeStepOnly] ;
+      For i In {1:NumMagnets}
+        Print[ f~{i}[Vol_Air], OnGlobal, Format Table, File > "F.dat"  ];
+        Print[ fx~{i}[Vol_Air], OnGlobal, Format Table, File > "Fx.dat",
+          SendToServer Sprintf("Output/Magnet %g/X force [N]", i), Color "Ivory"  ];
+        Print[ fy~{i}[Vol_Air], OnGlobal, Format Table, File > "Fy.dat",
+          SendToServer Sprintf("Output/Magnet %g/Y force [N]", i), Color "Ivory"  ];
+        Print[ fz~{i}[Vol_Air], OnGlobal, Format Table, File > "Fz.dat",
+          SendToServer Sprintf("Output/Magnet %g/Z force [N]", i), Color "Ivory"  ];
+      EndFor
+    }
+  }
+  { Name MagSta_a ; NameOfPostProcessing MagSta_a ;
+    Operation {
+      Print[ b,  OnElementsOf Vol_mu,  File "b.pos" ];
+      Echo[ Str["l=PostProcessing.NbViews-1;",
+		"View[l].ArrowSizeMax = 100;",
+		"View[l].CenterGlyphs = 1;",
+		"View[l].VectorType = 1;" ] ,
+	    File "tmp.geo", LastTimeStepOnly] ;
+      //Print[ br,  OnElementsOf Vol_Magnet,  File "br.pos" ];
+      //Print[ a,  OnElementsOf Vol_mu,  File "a.pos" ];
+      For i In {1:NumMagnets}
+      //Print[ un~{i}, OnElementsOf Domain, File "un.pos"  ];
+        Print[ f~{i}[Vol_Air], OnGlobal, Format Table, File > "F.dat"  ];
+        Print[ fx~{i}[Vol_Air], OnGlobal, Format Table, File > "Fx.dat",
+          SendToServer Sprintf("Output/Magnet %g/X force [N]", i), Color "Ivory"  ];
+        Print[ fy~{i}[Vol_Air], OnGlobal, Format Table, File > "Fy.dat",
+          SendToServer Sprintf("Output/Magnet %g/Y force [N]", i), Color "Ivory"  ];
+        Print[ fz~{i}[Vol_Air], OnGlobal, Format Table, File > "Fz.dat",
+          SendToServer Sprintf("Output/Magnet %g/Z force [N]", i), Color "Ivory"  ];
+      EndFor 
+    }
+  }
+}
+

MagneticForces/magnets_common.pro

+mm = 1.e-3;
+deg = Pi/180.;
+DefineConstant[
+  NumMagnets = {2, Min 1, Max 20, Step 1, Name "Parameters/0Number of magnets"}
+  Flag_InfiniteBox = {1, Choices{0,1}, Name "Infinite box/Add infinite box"}
+  Flag_FullMenu = {0, Choices{0,1}, Name "Parameters/Show all parameters"}
+  //DistToBox = {50*mm, Name "Infinite box/Distance to box [m]", Visible Flag_InfiniteBox}
+  //lcInf = {DistToBox, Name "Infinite box/Mesh size", Visible Flag_InfiniteBox}
+];
+
+For i In {1:NumMagnets}
+  DefineConstant[
+    X~{i} = {0, Min -100*mm, Max 100*mm, Step mm, Visible Flag_FullMenu,
+      Name Sprintf("Parameters/Magnet %g/0X position [m]", i) },
+    Y~{i} = { (i-1)*60*mm, Min -100*mm, Max 100*mm, Step mm,
+      Name Sprintf("Parameters/Magnet %g/0Y position [m]", i) },
+    Z~{i} = {0, Min -100*mm, Max 100*mm, Step mm, Visible Flag_FullMenu,
+      Name Sprintf("Parameters/Magnet %g/0Z position [m]", i) },
+
+    M~{i} = {(i==1)?0:1, Choices{0="Cylinder",1="Cube"},
+      Name Sprintf("Parameters/Magnet %g/00Shape", i)},
+
+    R~{i} = {20*mm, Min mm, Max 100*mm, Step mm,
+      Name Sprintf("Parameters/Magnet %g/1Radius [m]", i),
+      Visible (M~{i} == 0) },
+    L~{i} = {50*mm, Min mm, Max 100*mm, Step mm,
+      Name Sprintf("Parameters/Magnet %g/1Length [m]", i),
+      Visible (M~{i} == 0) },
+
+    Lx~{i} = {50*mm, Min mm, Max 100*mm, Step mm,
+      Name Sprintf("Parameters/Magnet %g/1X length [m]", i),
+      Visible (M~{i} == 1) },
+    Ly~{i} = {50*mm, Min mm, Max 100*mm, Step mm, Visible Flag_FullMenu,
+      Name Sprintf("Parameters/Magnet %g/1XY aspect ratio", i),
+      Visible (M~{i} == 1) },
+    Lz~{i} = {50*mm, Min mm, Max 100*mm, Step mm, Visible Flag_FullMenu,
+      Name Sprintf("Parameters/Magnet %g/1XZ aspect ration", i),
+      Visible (M~{i} == 1) },
+
+    Rx~{i} = {0, Min -Pi, Max Pi, Step Pi/180, Visible Flag_FullMenu,
+      Name Sprintf("Parameters/Magnet %g/2X rotation [deg]", i) },
+    Ry~{i} = {0, Min -Pi, Max Pi, Step Pi/180, Visible Flag_FullMenu,
+      Name Sprintf("Parameters/Magnet %g/2Y rotation [deg]", i) },
+    Rz~{i} = {0, Min -Pi, Max Pi, Step Pi/180,
+      Name Sprintf("Parameters/Magnet %g/2Z rotation [deg]", i) },
+
+    MUR~{i} = {(i==1)?1.:1000.,
+      Name Sprintf("Parameters/Magnet %g/3Mu relative []", i)},
+    BR~{i} = {(i==1)?1.0:0.0,
+      Name Sprintf("Parameters/Magnet %g/3Br [T]", i)}
+  ];
+EndFor
+
+//The geometrical parameters of the Infinite box. 
+DefineConstant[
+  xInt = {1, Name "Infinite box/xInt", Visible 0}
+  yInt = {1, Name "Infinite box/yInt", Visible 0}
+  zInt = {1, Name "Infinite box/zInt", Visible 0}
+  xExt = {xInt*2, Name "Infinite box/xExt", Visible 0}
+  yExt = {yInt*2, Name "Infinite box/yExt", Visible 0}
+  zExt = {zInt*2, Name "Infinite box/zExt", Visible 0}
+  xCnt = {0, Name "Infinite box/xCenter", Visible 0}
+  yCnt = {0, Name "Infinite box/yCenter", Visible 0}
+  zCnt = {0, Name "Infinite box/zCenter", Visible 0}
+];