clean up

Ccore/README.txt

+Start the optimization environment for the C-core model 
+with the command
+
+$ gmsh ccore.py

Ccore/ccore.geo

@@ -102,73 +102,4 @@ Physical Surface("COILN", 4) = 4;
 
 Physical Line("DIR", 11) = {1 ... 4};
 
-
-/*
-  The remainder of the file is a generic block of commands
-  for the automatic computation of velocity fields.
-  This velocity is used for the computation of sensitivities [1].
-  It should be copied unmodified at the bottom 
-  of all shape optimisation .geo files.
-
-   [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[
-  PerturbMesh = 0
-  VelocityTag = -1
-  Perturbation = 1e-10
-  modelpath = CurrentDir
-];
-
-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;
-  Save View[0] StrCat(modelpath, Sprintf("view_%g.msh",0));
-
-  // 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;
-  Save View[1] StrCat(modelpath, Sprintf("view_%g.msh",1));
-
-  // 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";
-  Save View[0] StrCat(modelpath, Sprintf("view_%g.msh",2));
-  Delete View[1];
-  Save View[0] StrCat(modelpath, Sprintf("velocity_%g.msh",VelocityTag));
-  SendToServer View[0] Sprintf("Optimization/Results/velocity_%g",VelocityTag);
-EndIf
+Include "OnelabOptimize.geo";

Ccore/ccore.pro

@@ -175,14 +175,16 @@ Function {
   Flux = NbTurns*Mmf/(RelGap+RelCore); // analytic approximation
   SetNumber("Model/Results/Flux analytic", Flux);
   //Flux = 0;
-  If( ParamIndex == 0 ) // D
+  If( SweepOnParameter == 1 ) // D
     dFluxdTau = Flux/E*Flux/(NbTurns*Mmf)*(1/mu_analytic-1/mu0);
-  ElseIf( ParamIndex == 1 ) // E
+  ElseIf( SweepOnParameter == 2 ) // E
     dFluxdTau = Flux/E*(1.+4*Flux/(NbTurns*Mmf*mu_analytic));
     If( Flag_Jfixed )
       dFluxdTau -= 2*Flux*G/CoilSection;
     EndIf
-  ElseIf( ParamIndex == 2 ) // F
+  ElseIf( SweepOnParameter == 3 ) // F
+    dFluxdTau = 0;
+  Else
     dFluxdTau = 0;
   EndIf
 

Ccore/ccore.py

+
+import OnelabOptimize
+OnelabOptimize.minimize('ccore')

Ccore/ccore_common.pro

 
-// Onelab parameters related with the optimisation
 
-ParamIndex = 
-  DefineNumber[0, Name "Optimization/Sweep on parameter",
-	       Choices {0="Optimize on D and E", 
-			1="D: air gap length", 
-	       		2="E: C-core thickness"} ];
-Optimizer = 
-  DefineNumber[0, Name "Optimization/Optimizer",
-	       Choices {0="scipy", 
-			1="MMA"}];
-StepsInRange = 
-  DefineNumber[20, Name "Optimization/Steps in range", Visible ParamIndex];
-
-w = 
-  DefineNumber[0, Name "Optimization/Results/w", 
-	       Graph "030000000000000000000000000000000000"];
-dwdt =
-  DefineNumber[0, Name Sprintf["Optimization/Results/dwdtau_%g",0], 
-	       Graph "000300000000000000000000000000000000"];
-
-
-
-// Design parameters 'D' and 'E' are defined
-// within a 'DefineConstant[]' statement to be contrallable 
-// from the calling optimisation loop with a '-setnumber' line argument
-// for the computation of the velocity field.
+// Define first the candidate design variables of the problem
+// with appropriate bounds and a default first guess value. 
+// Explain in the "Help" option their meaning. 
+// These definitions must mandatorily be enclosed 
+// in a DefineConstant[] statement.
 
 DefineConstant[
  D = DefineNumber[0.01, Name "Optimization/Parameters/D", 
-		  Min 0.001, Max 0.050]
+		  Min 0.001, Max 0.050, Help "Air gap length"]
  E = DefineNumber[0.1 , Name "Optimization/Parameters/E", 
-		  Min 0.01 , Max 0.15]
+		  Min 0.01 , Max 0.15, Help "C-core thickness"]
+ F = DefineNumber[0.01, Name "Optimization/Parameters/F", 
+		  Min 0.001, Max 0.090, Help "Core-coil gap"]
 ];
 
-// A number of geometrical parameters of the physical model
-// are also used in the definition of the finite element model
-// and the optimisation problem.  
-// Definitions for 'D' and 'E' must be in a 'DefineConstant[]' statement
-// for the same reason as above 
-// (they also must be overriden by the '-setnumber' line argument)
-// They are however all gathered in the '_common.pro' file 
-// and in the 'DefineConstant[]' statement for the sake of clarity 
-// and a tidy organisation of the model files.
+// "OnelabOptimize.py" needs to know the Onelab name of the design variables.
+// They are communicated as the "Choices" of a Onelab string variable 
+// with the reserved name "Optimization/Design parameters"
+DesignParameters =
+  DefineString["", Name "Optimization/Design parameters", Visible 0,
+	       Choices {"Optimization/Parameters/D",
+			"Optimization/Parameters/E",
+			"Optimization/Parameters/F"}];
+
+// "OnelabOptimize.py" can work in two modes:
+// - parameter sweep and optimization over one individual design parameter
+// - optimization over a multidimensional design space
+// The mode is selected by the user with the parameter "SweepOnParameter".
+
+// If "SweepOnParameter" is positive, the objective function and the sensibility 
+// of the selected parameter are evaluated for "StepsInRange" points 
+// distributed evenly over the parameter range.
+// When completed, python prints in the terminal the results of the sweep
+// with in particular a finite difference approximation of the sensitivity,
+// in order to validate the sensitivity analysis for that parameter.
+// The parameter sweep is followed by a minimization over the selected parameter,
+// and skipped if "StepsInRange" is set to zero. 
+
+// If "SweepOnParameter" is zero, the optimization is performed
+// over the design space determined by the indices given in the "Choices"
+// of the variable "DesignSpace". 
+// The indices of the design parameter correspond to their order
+// in the Choices list of "DesignParameters". 
+
+SweepOnParameter = 
+  DefineNumber[0, Name "Optimization/Sweep on parameter", Min 1, Max 3,
+              Choices {0="Optimize on selected design space",
+		       1="D",
+		       2="E",
+		       3="F"}];
+StepsInRange = 
+  DefineNumber[20, Name "Optimization/Steps in range", Visible SweepOnParameter];
+
+DesignSpace =
+  DefineNumber[0, Name "Optimization/Design space", Visible 0,
+	       Choices {1,2}, Min 1, Max 3];
+
+// available optimization algorithms
+Optimizer = 
+  DefineNumber[0, Name "Optimization/Algorithm",
+               Choices {0="scipy", 
+                        1="mma"}];
+w = 
+  DefineNumber[0, Name "Optimization/Results/w", 
+	       Graph "030000000000000000000000000000000000"];
 
 
 DefineConstant[