Team25

Team25/shape.geo

+DefineConstant[
+  mm = 1e-3
+  deg = Pi/180
+  R1 = 6*mm
+  L2 = 18*mm
+  L3 = 15*mm
+  L4 = L2/L3*Sqrt[L3^2-(12.5*mm)^2]
+  angleMeasure = 50*deg
+  lc = 12*mm
+  lcd = lc/15
+
+  // mesh perturb
+  PerturbMesh = 0
+  VelocityTag = -1
+  Perturbation = 1e-6
+  modelpath = CurrentDir
+];
+
+Point(1) = {0, 0, 0, lcd};
+
+For k In {0:10}
+  theta = k * 90/10 * deg;
+  DefineConstant[Rs~{k} = {R1, Name Sprintf["Constructive parameters/spline radius %g",k+1]}];
+  Point(2020+k) = {Rs~{k}*Cos[theta], Rs~{k}*Sin[theta], 0, lcd};  
+EndFor
+
+BSpline(4020) = {2020,2021, 2022, 2023};
+BSpline(4021) = {2023, 2024, 2025, 2026};
+BSpline(4022) = {2026, 2027, 2028, 2029, 2030};
+
+Point(26) = {L2/L3*Sqrt[L3^2-((12.5-2)*mm)^2], (12.5-2)*mm, 0, lcd};
+
+For k In {0:10}
+  theta = k * 35/10 * deg;
+  DefineConstant[Rso~{k} = {L2*L3/Sqrt[(L3*Cos[theta])^2+(L2*Sin[theta])^2], 
+    Name Sprintf["Constructive parameters/outer spline radius %g",k+1]}];
+  Point(6020+k) = {Rso~{k}*Cos[theta], Rso~{k}*Sin[theta], 0, lcd};  
+EndFor
+
+BSpline(4023) = {6020, 6021, 6022, 6023};
+BSpline(4024) = {6023, 6024, 6025, 6026};
+BSpline(4025) = {6026, 6027, 6028};
+BSpline(4026) = {6028, 6029, 6030, 26};
+
+Point(6) = {20*mm, 0, 0, lcd};
+Point(8) = {0, L3, 0, lc};
+Point(9) = {25*mm, 0, 0, lc};
+Point(10) = {163*mm, 0, 0, lc};
+Point(11) = {25*mm, 50*mm, 0, lc};
+Point(12) = {113*mm, 50*mm, 0, lc};
+Point(13) = {113*mm, (50+80)*mm, 0, lc};
+Point(14) = {0, (50+80)*mm, 0, lc};
+Point(15) = {0, 180*mm, 0, lc};
+Point(16) = {163*mm, 180*mm, 0, lc};
+Point(17) = {25*mm, (50+7.5)*mm, 0, lc};
+Point(18) = {113*mm, (50+7.5)*mm, 0, lc};
+Point(19) = {25*mm, (50+7.5+39)*mm, 0, lc};
+Point(20) = {113*mm, (50+7.5+39)*mm, 0, lc};
+
+Point(1021) = {(9.5+2.25)*Cos[angleMeasure]*mm, (9.5+2.25)*Sin[angleMeasure]*mm, 0, lcd};
+Point(1022) = {(9.5+2.25)*mm, 0, 0, lcd};
+
+Point(23) = {20*mm, 12.5*mm, 0, lcd};
+Point(24) = {20*mm-L4, 12.5*mm, 0, lcd};
+Point(25) = {20*mm-L4, (12.5-2)*mm, 0, lcd};
+
+Circle(1016) = {1022, 1, 1021};
+Line(2) = {9, 10};
+Line(3) = {10, 16};
+Line(4) = {16, 15};
+Line(5) = {15, 14};
+Line(6) = {14, 13};
+Line(8) = {12, 11};
+Line(9) = {11, 9};
+Line(10) = {17, 18};
+Line(11) = {12, 18};
+Line(12) = {18, 20};
+Line(13) = {20, 19};
+Line(14) = {19, 17};
+Line(15) = {20, 13};
+Line(18) = {1, 2020};
+Line(24) = {2020, 6020};
+Line(22) = {1, 2030};
+Line(26) = {6, 23};
+Line(27) = {24, 23};
+Line(28) = {24, 25};
+Line(29) = {25, 26};
+Line(30) = {6020, 6};
+Line(31) = {6, 9};
+Line(32) = {2030, 14};
+
+Curve Loop(1) = {-22, 18, 4020, 4021, 4022};
+Plane Surface(1) = {1};
+Curve Loop(2) = {4023,4024,4025,4026, -29, -28, 27, -26, -30};
+Plane Surface(2) = {2};
+Curve Loop(3) = {9, 2, 3, 4, 5, 6, -15, -12, -11, 8};
+Plane Surface(3) = {3};
+Curve Loop(4) = {13, 14, 10, 12};
+Plane Surface(4) = {4};
+Curve Loop(6) = {24, 4023,4024,4025,4026, -29, -28, 27, -26, 31, -9, -8, 11, -10, -14, -13, 15, -6, -32, -4020, -4021, -4022};
+Plane Surface(6) = {6};
+
+Physical Surface("pole", 1) = {3};
+Physical Surface("inductor", 2) = {4};
+Physical Surface("die molds", 3) = {1,2};
+Physical Surface("air", 5) = {6};
+Physical Line("dirichlet", 6) = {3,4,18,20,24,30,31,2};
+Physical Line("neuman", 7) = {22,32};
+
+If(PerturbMesh == 1)
+  Printf("Computing velocity field ...");
+  modelName = StrPrefix(StrRelative(General.FileName));  
+
+  SyncModel;
+  
+  // Merge the original mesh
+  Merge StrCat(modelpath, modelName, ".msh");
+
+  // create a view with the original node coordinates as a vector dataset
+  Plugin(NewView).NumComp = 3;
+  Plugin(NewView).Run;
+  Plugin(ModifyComponents).View = 0;
+  Plugin(ModifyComponents).Expression0 = "x";
+  Plugin(ModifyComponents).Expression1 = "y";
+  Plugin(ModifyComponents).Expression2 = "z";
+  Plugin(ModifyComponents).Run;
+
+  // relocate the vertices of the original mesh on the perturbed geometry
+  RelocateMesh Point "*";
+  RelocateMesh Line "*";
+  RelocateMesh Surface "*";
+
+  // Create a view with the perturbed node coordinates as vector dataset
+  Plugin(NewView).NumComp = 3;
+  Plugin(NewView).Run;
+  Plugin(ModifyComponents).View = 1;
+  Plugin(ModifyComponents).Expression0 = "x";
+  Plugin(ModifyComponents).Expression1 = "y";
+  Plugin(ModifyComponents).Expression2 = "z";
+  Plugin(ModifyComponents).Run;
+
+  // compute the velocity field by subtracting the two vector datasets
+  Plugin(ModifyComponents).View = 0;
+  Plugin(ModifyComponents).OtherView = 1;
+  Plugin(ModifyComponents).Expression0 = Sprintf("(w0 - v0)/(%g)", Perturbation);
+  Plugin(ModifyComponents).Expression1 = Sprintf("(w1 - v1)/(%g)", Perturbation);
+  Plugin(ModifyComponents).Expression2 = Sprintf("(w2 - v2)/(%g)", Perturbation);
+  Plugin(ModifyComponents).Run;
+  View.Name = "velocity";
+  Delete View[1];
+  SendToServer View[0] Sprintf("Optimization/Results/velocity_%g",VelocityTag); // Hidden
+EndIf

Team25/shape.pro

+/* -------------------------------------------------------------------
+   Tutorial: Optimization of Die Press Model (TEAM25)
+
+   Features:
+   - Nonlinear magnetostatic model of the die press
+   - CAD parameters as ONELAB parameters
+   - Direct approach of the shape sensitivity analysis 
+     based on the Lie derivative
+
+   To compute the solution in a terminal:
+       getdp shape.pro -solve GetPerformancesAndGradients
+
+   To compute the solution interactively from the Gmsh GUI:
+       File > Open > shape.pro
+       Run (button at the bottom of the left panel)
+   ------------------------------------------------------------------- */
+
+/* This model computes the static magnetic field produced by a DC current. This
+   corresponds to a "magnetostatic" physical model, obtained by combining the
+   time-invariant Maxwell-Ampere equation (curl h = js, with h the magnetic
+   field and js the source current density) with Gauss' law (Div b = 0, with b
+   the magnetic flux density) and the magnetic constitutive law (b = mu h, with
+   mu the magnetic permeability).
+
+   Since Div b = 0, b can be derived from a vector magnetic potential a, such
+   that b = curl a. Plugging this potential in Maxwell-Ampere's law and using
+   the constitutive law leads to a vector Poisson equation in terms of the
+   magnetic vector potential: curl(nu curl a) = js, where nu = 1/mu is
+   the reluctivity.
+
+   It is useful to share parameters 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 "DefineConstant" statement, 
+   which can be invoked in the .geo and .pro files.
+
+   The optimization part is organized as follows:
+   
+   (1) We replaced both curves g-h and i-j by nurbs which have a certain 
+   of control points. We set the shape design variables as the distance 
+   of these control points from the origin point. Their actual value is
+   provided by shape.py.
+   
+   (2) We make use of the solution of the nonlinear magnetostatic problem,
+   for a given CAD configuration, to compute the objective function, 
+     w = Sum_{i=1}^{10} ((bpx_i-box_o)^2 + (bpy_i-boy_o)^2),
+   as defined in TEAM25.
+   
+   (3) We compute the derivative of the flux density along the velocity 
+   field generated by the perturbation of each design variable, by means
+   of the Lie derivative. The derivative of he objective, w, is therefore
+   obtained by applying the chain rule. This is known as the direct approach
+   of the sensitivity analysis, see [1].
+
+   The value of the objective function as well as is derivative with respect 
+   to each design variable are provided to the ONELAB database. They can be 
+   accessed by shape.py.
+
+   [1]:Erin Kuci, François Henrotte, Pierre Duysinx, and Christophe Geuzaine. 
+   ``Design sensitivity analysis for shape optimization based on the Lie derivative'',
+   Computer Methods in Applied Mechanics and Engineering 317 (2017), pp. 702 –722. 
+   issn: 0045-7825.  
+*/
+
+DefineConstant [
+  Opt_ResDir = "res_opt/"
+  Opt_ResDir_Onelab = "Optimization/Results/"
+  Model_smallAT = {1, Name "Model/Small Ampere-Turn", Choices{0,1}} 
+  Flag_PrintLevel = 1
+  VelocityTag = -1
+  // Nonlinear iterations
+  Nb_max_iter = 30
+  relaxation_factor = 1
+  stop_criterion = 1e-5
+];
+
+Group {
+  // One starts by giving explicit meaningful names to the Physical regions
+  // defined in the "shape.geo" mesh file. 
+  Coil = Region[2];
+  Core = Region[{1,3}];
+
+  // Then one difines abstract regions so as to allow a generic definition of the
+  // FunctionSpace, Formulation and PostProcessing fields with no reference to
+  // model-specific Physical regions.
+  Domain = Region[{Core,Coil,5}]; 
+  NoFlux = Region[6];
+}
+
+Function{
+  // Material law (here the magnetic reluctivity) is defined as piecewise
+  // function (note the square brackets), in terms of the above defined
+  // physical regions. 
+
+  mu0 = 4 * Pi * 1e-7;
+  Mat1_b = {0.0,0.11,0.18,0.28,0.35,0.74,0.82,0.91,0.98,1.02,1.08,1.15,
+    1.27,1.32,1.36,1.39,1.42,1.47,1.51,1.54,1.56,1.60,1.64,1.72};
+  Mat1_h = {0.0,140.0,178.0,215.0,253.0,391.0,452.0,529.0,596.0,677.0,774.0,902.0,
+    1164.0,1299.0,1462.0,1640.0,1851.0,2262.0,2685.0,3038.0,3395.0,4094.0,4756.0,7079.0};
+  Mat1_b2 = List[Mat1_b]^2;
+  Mat1_nu = List[Mat1_h]/List[Mat1_b];
+  Mat1_nu(0) = Mat1_nu(1);
+  Mat1_nu_b2  = ListAlt[Mat1_b2, Mat1_nu];
+  nu[Core] = InterpolationLinear[ SquNorm[$1] ]{List[Mat1_nu_b2]} ;
+  nu[Region[{Domain,-Core}]] = 1/mu0; // linear
+
+  dnudb2[] = dInterpolationLinear[SquNorm[$1]]{List[Mat1_nu_b2]} ;
+  dnudb_1[] = 2.0*dInterpolationLinear[SquNorm[$1]]{List[Mat1_nu_b2]}*SquDyadicProduct[$1];
+  dhdb_NL[Core] = 2*dnudb2[$1#1] * SquDyadicProduct[#1];
+
+  // This is the current density which feeds the inductor.
+  js[Coil] = Vector[0, 0, (Model_smallAT==1?4253:17500)/SurfaceArea[]{2}];
+
+  // This is the desired x and y components of the induction field 
+  // in the cavity as defined in TEAM25. 
+  bx[] = (Model_smallAT==1?0.35:1.5) * X[]/Sqrt[X[]^2+Y[]^2];
+  by[] = (Model_smallAT==1?0.35:1.5) * Y[]/Sqrt[X[]^2+Y[]^2];
+}
+
+Jacobian {
+  { Name Vol; Case { { Region All ; Jacobian Vol; } } }
+  { Name Sur; Case { { Region All ; Jacobian Sur; } } }
+}
+
+Integration{
+  { Name Int2;
+    Case{
+      { Type Gauss;
+        Case{
+          { GeoElement Line;       NumberOfPoints 3; }
+          { GeoElement Triangle;   NumberOfPoints 4; }
+          { GeoElement Quadrangle; NumberOfPoints 4; }
+        }
+      }
+    }
+  }
+}
+
+// -------------------------------------------------------------------------
+// The weak formulation for this problem is derived. The fields are 
+// vector-valued, and we have one linear (source) term in addition 
+// to the bilinear term. 
+
+Constraint{
+  { Name Dirichlet;
+    Case{
+      { Region NoFlux; Type Assign; Value 0; }
+    }
+  }
+}
+
+FunctionSpace {
+  { Name HCurl2D; Type Form1P;
+    BasisFunction{
+      { Name se1; NameOfCoef ae1; Function BF_PerpendicularEdge;
+        Support Region[Domain]; Entity NodesOf[All]; }
+    }
+    Constraint{
+      { NameOfCoef ae1; EntityType NodesOf; NameOfConstraint Dirichlet; }
+    }
+  }
+}
+
+Formulation{
+  { Name MVP2D; Type FemEquation;
+    Quantity{
+      { Name a; Type Local; NameOfSpace HCurl2D; }
+    }
+    Equation{
+      Galerkin{ [nu[{d a}] * Dof{d a}, {d a}];
+        In Domain;  Jacobian Vol; Integration Int2; }
+      Galerkin { JacNL [ dhdb_NL[{d a}] * Dof{d a} , {d a} ];
+        In Core; Jacobian Vol ; Integration Int2 ; }
+      Galerkin { [-js[], {a}];
+        In Coil; Jacobian Vol; Integration Int2; }
+    }
+  }
+}
+
+// -------------------------------------------------------------------------
+// Handling of the velocity field (through a fully-fixed function space 
+// defined on Domain). 
+
+Function{
+  For i In {1:3}
+    l_v~{i} = ListFromServer[Sprintf["Optimization/Results/velocity_%g_%g",VelocityTag,i]];
+    velocity~{i}[] = ValueFromIndex[]{l_v~{i}()};
+  EndFor
+}
+
+Constraint {
+  For i In {1:3}
+    { Name velocity~{i} ;
+      Case {
+        { Region Domain ; Value velocity~{i}[]; }
+      }
+    }
+  EndFor
+}
+
+FunctionSpace {
+  For i In {1:3}
+    { Name H_v~{i} ; Type Form0;  
+      BasisFunction {
+        { Name sn ; NameOfCoef un ; Function BF_Node ; 
+          Support Domain; Entity NodesOf[ All ] ; }
+      }
+      Constraint {
+        { NameOfCoef un ; EntityType NodesOf ; NameOfConstraint velocity~{i}; }
+      }
+    }
+  EndFor
+}
+
+// -------------------------------------------------------------------------
+// Magnetic vector potential handling (through a fully-fixed function space 
+// defined on Domain). 
+
+Function{
+  l_a = ListFromServer[StrCat[Opt_ResDir_Onelab,"a"]];
+  aFromServer[] = ValueFromIndex[]{l_a()};
+}
+
+Constraint {
+  { Name aFromServer;
+    Case {
+      { Region Domain ; Value aFromServer[]; }
+      { Region NoFlux; Type Assign; Value 0; }
+    }
+  }
+}
+
+FunctionSpace{
+  { Name HCurl2D_a_fullyfixed; Type Form1P;
+    BasisFunction{
+      { Name se1; NameOfCoef ae1; Function BF_PerpendicularEdge;
+        Support Region[Domain]; Entity NodesOf[All]; }
+    }
+    Constraint{
+      { NameOfCoef ae1; EntityType NodesOf; NameOfConstraint aFromServer; }
+    }
+  }
+}
+
+// -------------------------------------------------------------------------
+// Handling of the direct sensitivity problem (depending on a given
+// design variable represeted by the velocity field).
+
+Function{
+  dV[] = Transpose[Tensor[CompX[$1],CompX[$2],0,
+            CompY[$1],CompY[$2],0,
+            CompZ[$1],CompZ[$2],0]];
+
+  // Lie derivative of H(B) where B is a 2-form and H a 1-form ($1:{d a}, $2:dV)
+  LieOf_HB[] = nu[$1] * (Transpose[$2] * $1 - TTrace[$2] * $1 + $2 * $1) ;
+  LieOf_HB_NL[] = dhdb_NL[$1] * (Transpose[$2] - TTrace[$2] * TensorDiag[1,1,1]) * $1;
+
+  // Lie derivative of the 2-form Js ($1:Js[], $2:dV)
+  LieOf_Js[] = TTrace[$2]*js[] - Transpose[$2]*js[];
+}
+
+FunctionSpace{
+  { Name HCurl2D_Lva; Type Form1P ;
+    BasisFunction {
+      { Name se1 ; NameOfCoef ae1 ; Function BF_PerpendicularEdge ;
+        Support Region[{ Domain}] ; Entity NodesOf [ All ] ; }
+    }
+    Constraint {
+      { NameOfCoef ae1 ; EntityType NodesOf ; NameOfConstraint Dirichlet; }
+    }
+  }
+}
+
+Formulation {
+  { Name DirectFormulationOf_MagSta_a_2D; Type FemEquation ;
+    Quantity {
+      { Name a; Type Local; NameOfSpace HCurl2D_a_fullyfixed; }
+      { Name Lva; Type Local; NameOfSpace HCurl2D_Lva; }
+      For i In {1:3}
+        { Name v~{i} ; Type Local ; NameOfSpace H_v~{i}; }
+      EndFor
+    }
+    Equation {
+      Galerkin { [ nu[{d a}] * Dof{d Lva}, {d Lva} ];
+        In Domain; Jacobian Vol; Integration Int2; }
+      Galerkin { [LieOf_HB[{d a},dV[{d v_1},{d v_2},{d v_3}]],{d Lva}];
+        In Domain; Jacobian Vol; Integration Int2; }
+      Galerkin { [LieOf_HB_NL[{d a},dV[{d v_1},{d v_2},{d v_3}]],{d Lva}];
+        In Core; Jacobian Vol; Integration Int2; }
+      Galerkin { [ LieOf_Js[js[], dV[{d v_1},{d v_2},{d v_3}]], {Lva} ] ;
+        In Coil ; Jacobian Vol; Integration Int2; }
+
+      Galerkin { [ 0*Dof{a}, {a} ] ;
+        In Domain; Jacobian Vol ; Integration Int2 ; }
+
+      For i In {1:3}
+        Galerkin { [ 0*Dof{v~{i}}, {v~{i}} ] ;
+          In Domain; Jacobian Vol ; Integration Int2 ; }
+      EndFor
+    }
+  }
+}
+
+// -------------------------------------------------------------------------
+// In the following, we create the post-operations necessary
+// for the computation of objective and constraints, 
+// as well as their sensitivities. 
+
+PostProcessing {
+  { Name ObjectiveConstraints; NameOfFormulation MVP2D;
+    Quantity {
+      { Name w; Value{ Term{[(CompX[{d a}]-bx[])^2+(CompY[{d a}]-by[])^2 ]; In Domain; Jacobian Vol; } } }
+      { Name az; Value{ Term{ [CompZ[{a}]]; In Domain; Jacobian Vol; } } }
+      { Name b; Value{ Term{ [{d a}]; In Domain; Jacobian Vol; } } }
+      { Name bMag; Value{ Term{ [Norm[{d a}]]; In Domain; Jacobian Vol; } } }
+      { Name js; Value{ Term{ [js[]]; In Coil; Jacobian Vol; } } }
+      { Name mur; Value{ Term{ [1/(nu[{d a}]*mu0)]; In Domain; Jacobian Vol; } } }
+    }
+  }
+}
+
+PostOperation Get_ObjectiveConstraints UsingPost ObjectiveConstraints {
+  CreateDir[Opt_ResDir];
+  Print[w, 
+    OnGrid {(9.5e-3+2.25e-3)*Cos[$A],(9.5e-3+2.25e-3)*Sin[$A],0}{0:50*Pi/180:5*Pi/180,0,0},
+    Format SimpleTable, File StrCat[Opt_ResDir,"w.txt"]];
+  Print[bMag,  
+    OnGrid {(9.5e-3+2.25e-3)*Cos[$A],(9.5e-3+2.25e-3)*Sin[$A],0}{0:50*Pi/180:5*Pi/180,0,0},
+    Format SimpleTable, File StrCat[Opt_ResDir,"bMag.txt"]];
+  Print[az, OnElementsOf Domain, Format NodeTable, File "", 
+    SendToServer StrCat[Opt_ResDir_Onelab,"a"], Hidden 1];
+  Print[bMag, OnElementsOf Domain, File StrCat[Opt_ResDir,"az.pos"]];
+  //Print[az, OnElementsOf Domain, File StrCat[Opt_ResDir,"az.pos"]];
+  If(Flag_PrintLevel>5)
+    Print[mur, OnElementsOf Domain, File StrCat[Opt_ResDir,"mur.pos"]];
+    Print[az, OnElementsOf Domain, File StrCat[Opt_ResDir,"az.pos"]];
+    Print[a, OnElementsOf Domain, File StrCat[Opt_ResDir,"a.pos"]];
+    Print[b, OnElementsOf Domain, File StrCat[Opt_ResDir,"b.pos"]];
+    Print[bMag, OnElementsOf Domain, File StrCat[Opt_ResDir,"bMag.pos"]];
+    Print[js, OnElementsOf Coil, File StrCat[Opt_ResDir,"js.pos"]];
+  EndIf
+}
+
+// -------------------------------------------------------------------------
+// Sensitivity of w based on a direct method
+
+PostProcessing{
+  { Name Direct_MagSta; NameOfFormulation DirectFormulationOf_MagSta_a_2D;
+    PostQuantity {
+      { Name Lie_w; Value { Term { [ 2*(CompX[{d a}]-bx[])*CompX[{d Lva}]
+        +2*(CompY[{d a}]-by[])*CompY[{d Lva}] ] ; In Domain; Jacobian Vol; }}}
+    }
+  }
+}
+
+PostOperation Get_GradOf_w UsingPost Direct_MagSta {
+  Print[Lie_w, OnGrid {(9.5e-3+2.25e-3)*Cos[$A],(9.5e-3+2.25e-3)*Sin[$A], 0} { 0:50*Pi/180:5*Pi/180, 0, 0 },
+    Format SimpleTable, File StrCat[Opt_ResDir,Sprintf["Grad_w_wrt_dv_%g.txt",VelocityTag]]];
+}
+
+// Show useful data
+PostProcessing {
+  { Name Show_shape; NameOfFormulation DirectFormulationOf_MagSta_a_2D;
+    PostQuantity {
+      { Name VV; Value{Term{[ Vector[{v_1},{v_2},{v_3}] ];In Domain ; Jacobian Vol;}}}
+      { Name Lie_az; Value { Term { [ CompZ[{Lva}] ] ; In Domain ; Jacobian Vol; }}}
+    }
+  }
+}
+
+PostOperation Show_shape UsingPost Show_shape{
+  CreateDir[Opt_ResDir];
+  Print[ VV, OnElementsOf Domain,File StrCat[Opt_ResDir, "velocity.pos"] ]; 
+  Print[ Lie_az, OnElementsOf Domain,File StrCat[Opt_ResDir, "Lie_az.pos"] ];
+}
+
+// -------------------------------------------------------------------------
+// The following resolution solves the physical problem
+// to obtain the objective and the constraints.
+
+Resolution {
+  { Name GetPerformances;
+    System {
+      { Name SYS; NameOfFormulation MVP2D;}
+    }
+    Operation {
+      // Initialization of the systems
+      InitSolution[SYS];
+
+      // Solve the nonlinear magnetostatic problem
+      IterativeLoop[Nb_max_iter, stop_criterion, relaxation_factor]{
+        GenerateJac[SYS]; 
+        SolveJac[SYS];
+        Evaluate[ SetNumberRunTime[$Iteration]{"iterNL"} ]; 
+        Evaluate[ SetNumberRunTime[$Residual]{"residualNL"} ]; 
+      }
+
+      // Compute the objective function and the constraints
+      PostOperation[Get_ObjectiveConstraints];
+    }
+  }
+}
+
+// -------------------------------------------------------------------------
+// This resolution solves the direct problem to compute the sensitivity
+// of the induction flux with respect to a given design variable.
+
+Resolution {
+  { Name GetGradient_wrt_dv;
+    System {
+      { Name DIR; NameOfFormulation DirectFormulationOf_MagSta_a_2D; }
+    }
+    Operation {
+      InitSolution[DIR];
+      Generate[DIR];
+      Solve[DIR]; 
+      PostOperation[Get_GradOf_w];
+    }
+  }
+}

Team25/shape.py

+# Open this file with Gmsh (interactively with File->Open->penning.py, or on the
+# command line with 'gmsh penning.py')
+#
+from shutil import copyfile
+import numpy as np
+import optlab
+import onelab
+
+c = onelab.client(__file__)
+
+# get gmsh and getdp locations from Gmsh options
+mygmsh = c.getString('General.ExecutableFileName')
+mygetdp = ''
+for s in range(9):
+   n = c.getString('Solver.Name' + str(s))
+   if(n == 'GetDP'):
+      mygetdp = c.getString('Solver.Executable' + str(s))
+      break
+if(not len(mygetdp)):
+   c.sendError('This appears to be the first time you are trying to run GetDP')
+   c.sendError('Please run a GetDP model interactively once with Gmsh to ' +
+               'initialize the solver location')
+   exit(0)
+
+c.sendInfo('Will use gmsh={0} and getdp={1}'.format(mygmsh, mygetdp))
+
+# append correct path to model file names
+file_geo = c.getPath('shape.geo')
+file_msh = c.getPath('shape.msh')
+file_pro = c.getPath('shape.pro')
+
+# build command strings with appropriate parameters and options
+mygetdp += ' -p 0 ' + file_pro + ' -msh ' + file_msh + ' '
+mygmsh += ' -p 0 -v 2 -parametric ' + file_geo + ' '
+
+# load geometry in GUI to have something to look at
+c.openProject(file_geo)
+
+# dry getdp run (without -solve or -pos option) to get model parameters in the GUI
+c.runSubClient('myGetDP', mygetdp)
+
+# define now optimization parameters
+# some of them as Onelab parameter, to be editable in the GUI
+maxIter = c.defineNumber('Optimization/00Max iterations', value=100)
+maxChange = c.defineNumber('Optimization/01Max change', value=1e-5)
+
+# end of check phase. We're done if we do not start the optimization
+if c.action == 'check' :
+   exit(0)
+
+x = {}
+for k in xrange(11):
+    xe = 6.0e-3
+    x[k] = ['Constructive parameters/spline radius '+str(k+1),xe, xe-4.0e-3, xe+4.0e-3, 'Rs_'+str(k)]
+sizeX = len(x.keys())
+
+L2 = 18.0e-3;L3 = 15.0e-3
+for k in xrange(11):
+    theta = float(k) * 35.0/10.0 * np.pi/180.0
+    xe = L2*L3/((L3*np.cos(theta))**2.0+(L2*np.sin(theta))**2.0)**0.5
+    x[sizeX+k] = ['Constructive parameters/outer spline radius '+str(k+1), xe, xe-4.0e-3, xe+1.5e-3, 'Rso_'+str(k)]
+
+def readSimpleTable(path):
+    with open(path) as FileObj:
+        return np.array([float(lines.split()[-1]) for lines in FileObj])
+
+def setNodeTable(var, nodes, val):
+    data = np.ravel(np.column_stack((nodes, val)))
+    data = np.insert(data, 0, float(len(nodes)))
+    c.setNumber(var, values=data, visible=0)
+
+def getVelocityField(xvar, perturb=1e-6):
+    for id, var in xvar.iteritems():
+        c.runNonBlockingSubClient('myGmsh', mygmsh \
+          + ' -setnumber PerturbMesh 1'\
+          + ' -setnumber VelocityTag '+str(id)\
+          + ' -setnumber '+var[-1]+' '+str(var[1]+perturb)\
+          + ' -run')
+    c.waitOnSubClients()
+    # send the x,y,z components of each velocity field
+    for label, var in x.iteritems():
+        d = c.getNumbers('Optimization/Results/velocity_'+str(label))
+        if label==0:nodes = np.array(d[1::4])
+        setNodeTable('Optimization/Results/velocity_'+str(label)+'_1',nodes,np.array(d[2::4]))
+        setNodeTable('Optimization/Results/velocity_'+str(label)+'_2',nodes,np.array(d[3::4]))
+        setNodeTable('Optimization/Results/velocity_'+str(label)+'_3',nodes,np.array(d[4::4]))
+        c.setNumber(var[0],value=x[label][1])
+
+# number of design variables and number of constraints
+numVariables = len(x.keys())
+
+# initial value of the variables in the design space,
+# the lower bound for design variables,
+# as well as the upper bound for design variables.
+initialPoint = np.zeros(numVariables)
+lowerBound = np.zeros(numVariables)
+upperBound = np.zeros(numVariables)
+for label, var in x.iteritems():
+  initialPoint[label] = var[1]
+  lowerBound[label] = var[2]
+  upperBound[label] = var[3]
+
+# Initialize the MMA optimizer
+optlab.mma.initialize(initialPoint, lowerBound, upperBound)
+
+# Set some options for MMA
+optlab.mma.option.setNumber('General.Verbosity', 4)
+optlab.mma.option.setNumber('General.SaveHistory', 0)
+optlab.mma.option.setNumber('SubProblem.isRobustMoveLimit', 1)
+optlab.mma.option.setNumber('SubProblem.isRobustAsymptotes', 1)
+optlab.mma.option.setNumber('SubProblem.type', 0)
+optlab.mma.option.setNumber('SubProblem.addConvexity', 0)
+optlab.mma.option.setNumber('SubProblem.asymptotesRmax', 100.0)
+optlab.mma.option.setNumber('SubProblem.asymptotesRmin', 0.001)
+
+# Get iteration count (here it will be 1 - could be different in case of restart)
+it = optlab.mma.getOuterIteration()
+change = 1.0
+
+while it <= maxIter and c.getString('shape/Action') != 'stop':
+
+    c.setNumber('Optimization/01Current iteration', value=it, readOnly=1,
+                attributes={'Highlight':'LightYellow'})
+
+    # get (copy of) current point
+    xFromMMA = np.array(optlab.mma.getCurrentPoint())
+
+    # send the current point to GetDP model
+    for label, var in x.iteritems():
+        x[label][1] = xFromMMA[label]
+        c.setNumber(var[0],value=x[label][1])
+
+    # reload the geometry in GUI to see the geometrical changes
+    c.openProject(file_geo)
+
+    # mesh the geometry
+    c.runSubClient('myGmsh', mygmsh + ' -2')
+
+    # generate the velocity field of each design variable
+    getVelocityField(x)
+
+    # compute objective function and constraints 
+    c.runSubClient('myGetDP', mygetdp + '-solve GetPerformances')
+
+    # as well as their sensitivity with respect to design variables at `x'
+    for dv, var in x.iteritems():
+        c.runNonBlockingSubClient('myGetDP', mygetdp \
+          + ' -setnumber VelocityTag '+str(dv)\
+          + ' -solve GetGradient_wrt_dv')
+    c.waitOnSubClients()
+
+    # get the value of the objective function and of the constraints
+    # as well as their sensitivity with respect to design variables at `x'
+    objective = np.sum(readSimpleTable('res_opt/w.txt'))
+    constraints = np.array([np.sum(xFromMMA)/100.0-1.0])
+    grad_objective = np.asarray([np.sum(readSimpleTable('res_opt/Grad_w_wrt_dv_'+str(dv)+'.txt'))\
+        for dv in xrange(numVariables)])
+    grad_constraints = np.ones(numVariables)/100.0	     
+
+    if it == 1: fscale = 1.0 / objective
+    objective *= fscale
+    grad_objective *= fscale
+
+    print '*'*50
+    print 'iteration: ', it
+    print 'change: ', change
+    print 'current point:', xFromMMA
+    print 'objective:', objective
+    print 'constraints:', constraints
+    c.sendInfo('Optimization: it. {}, obj. {}, constr. {}'.format(it,objective,constraints[0]))
+    #print 'gradient of objective:', grad_objective
+    #print 'gradient of constraints:', grad_constraints
+
+    # call MMA and update the current point
+    optlab.mma.setMoveLimits(lowerBound, upperBound, 1.0e-4)
+    optlab.mma.updateCurrentPoint(constraints,grad_objective,grad_constraints)
+    change = optlab.mma.getDesignChange()
+    it = optlab.mma.getOuterIteration()
+
+# This should be called at the end
+optlab.mma.finalize()