diff --git a/Team25/shape.geo b/Team25/shape.geo
new file mode 100644
index 0000000000000000000000000000000000000000..71d2bab2f61c3410987a2266fe0885da52c9d7df
--- /dev/null
+++ b/Team25/shape.geo
@@ -0,0 +1,152 @@
+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
diff --git a/Team25/shape.pro b/Team25/shape.pro
new file mode 100644
index 0000000000000000000000000000000000000000..6cdaf8f3f1a13fcfb29669cf248001e7e3902b75
--- /dev/null
+++ b/Team25/shape.pro
@@ -0,0 +1,421 @@
+/* -------------------------------------------------------------------
+ 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];
+ }
+ }
+}
diff --git a/Team25/shape.py b/Team25/shape.py
new file mode 100644
index 0000000000000000000000000000000000000000..6a01cc95144b2efda7e138ef8e9d9a37aabaddfe
--- /dev/null
+++ b/Team25/shape.py
@@ -0,0 +1,181 @@
+# 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()