diff --git a/MagneticForces/InfiniteBox.geo b/MagneticForces/InfiniteBox.geo
new file mode 100644
index 0000000000000000000000000000000000000000..0d9cb8523269bf9301be47852689a96739e39592
--- /dev/null
+++ b/MagneticForces/InfiniteBox.geo
@@ -0,0 +1,149 @@
+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] };
+
+
diff --git a/MagneticForces/magnets.geo b/MagneticForces/magnets.geo
new file mode 100644
index 0000000000000000000000000000000000000000..89143ab0700d57d224c62d711b7d2cae74503a83
--- /dev/null
+++ b/MagneticForces/magnets.geo
@@ -0,0 +1,94 @@
+// 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[] };
diff --git a/MagneticForces/magnets.pro b/MagneticForces/magnets.pro
new file mode 100644
index 0000000000000000000000000000000000000000..3421498848cd9eb7e1c424d4b8056ba1c6ff251a
--- /dev/null
+++ b/MagneticForces/magnets.pro
@@ -0,0 +1,380 @@
+/* -------------------------------------------------------------------
+ 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
+ }
+ }
+}
+
diff --git a/MagneticForces/magnets_common.pro b/MagneticForces/magnets_common.pro
new file mode 100644
index 0000000000000000000000000000000000000000..3c75dce75aa16ece8ea75c1329382fc82947e957
--- /dev/null
+++ b/MagneticForces/magnets_common.pro
@@ -0,0 +1,65 @@
+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}
+];