added comments and explanations

Ccore/ccore.geo

@@ -106,8 +106,15 @@ 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[

Ccore/ccore.pro

+/* -------------------------------------------------------------------
+   Tutorial: Optimization of ferromagnetic C-core
+
+   Features:
+   - Nonlinear magnetostatic model of a C-core
+   - Imposed current density or imposed current
+   - Air gap and core width as design parameters
+   - Prescribed By in the air gap as objective function 
+   - Direct approach of the shape sensitivity analysis
+   - Sensivity analysis with 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 tutorial deals with shape sensitivity analysis in Magnetostatics.
+  The model consists of a ferromagnetic C-core powered
+  by a coil with 'NbTurns' windings suplied with a DC current. 
+  The ferromagnetic behavior of the C-core material can be 
+  represented by either a linear of a nonlinear saturable law. 
+  The power supply is set to maintain either a constant total current
+  or a constant current density as the design parameters vary. 
+  The optimisation problem aims at having a prescribed value
+  for the By component of the magnetic flux density
+  in the middle of the air gap.
+  The two design variables selected for achieving this objective are
+  - the air gap width 'D',
+  - the core width 'E'.
+
+  The 'shp.py' is organised so as to allow either 
+  - to check the correctness of the sensitivity computation,
+  - effectively solve the optimisation problem.
+
+  The computation of sensitivities is a delicate step of the optimisation,
+  which can however be one hundred percent validated by finite difference.
+  This tutorial is thus dedicated to both the explanation and verification
+  of the direct computation of sensitivities in Magnetostatics.
+
+  Communication between the optimisation loop ('shp.py')
+  and the Onelab model is done via the Onelab database. 
+  Design variables are defined as Onelab variables whose name
+  has the reserved prefix 'Optimization/Parameters'.
+  Similarly, velocities and sensibilities are also exchanged as 
+  Onelab variables with reserved prefixes 'Optimization/Results/*'.
+
+  Finally, an extra mechanism is implemented 
+  to apply a temporary small perturbation of a design parameter 
+  in order to compute the associated velocity field
+  without interfering with the value stored in the Onelab database.
+  This affectation is done with a '-setnumber' line argument 
+  when calling Gmsh. 
+  In consequence, the definition of the design variables
+  as Onelab parameters must be enclosed in a 'DefineConstant[]' statement.
+  Remember that affectations within a 'DefineConstant[]'
+  are ignored if the variable already exists because,
+  and in particular if it has been set by a '-setnumber' line argument.
+ */
+
 Include "ccore_common.pro";
 
 DefineConstant[
@@ -8,6 +71,7 @@ DefineConstant[
   OptiIterNumber = 0
 ];
 
+// Onelab parameters of the C-core model 
 
 NL_tol_abs = 1e-8; 	// absolute tolerance on residual for noninear iterations
 NL_tol_relax = 1.0; 	// relaxation on residual for noninear iterations
@@ -37,9 +101,7 @@ Else
   CurrentDensity = Mmf / CoilSection;
 EndIf
 
-By_target = 1.;
-Px_target = CoreX + A - E/2 ;
-Py_target = CoreY + B/2 ;
+
 
 Group {
   // Physical regions (in capital letters):

Ccore/ccore_common.pro

 
-/*
-  This tutorial deals with shape sensitivity analysis in Magnetostatics.
-  The model consists of a ferromagnetic C-core powered
-  by a coil with 'NbTurns' windings suplied with a DC current. 
-  The ferromagnetic behavior of the C-core material can be 
-  represented by either a linear of a nonlinear saturable law. 
-  The power supply is set to maintain either a constant total current
-  or a constant current density as the design parameters vary. 
-
-  Three geometrical parametrs are 
-
- */
-
+// Onelab parameters related with the optimisation
 
 ParamIndex = 
   DefineNumber[0, Name "Optimization/Sweep on parameter",
-	       Choices {0="None", 
+	       Choices {0="Optimize on D and E", 
 			1="D: air gap length", 
-	       		2="E: C-core thickness",
-	       		3="F: Core-coil gap"} ];
+	       		2="E: C-core thickness"} ];
 StepsInRange = 
   DefineNumber[20, Name "Optimization/Steps in range", Visible ParamIndex];
 
-
 w = 
   DefineNumber[0, Name "Optimization/Results/w", 
 	       Graph "030000000000000000000000000000000000"];
@@ -32,15 +18,28 @@ dwdt =
 
 
 
+// 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.
+
 DefineConstant[
  D = DefineNumber[0.01, Name "Optimization/Parameters/D", 
 		  Min 0.001, Max 0.050]
  E = DefineNumber[0.1 , Name "Optimization/Parameters/E", 
 		  Min 0.01 , Max 0.15]
- F = DefineNumber[0.01, Name "Optimization/Parameters/F", 
-		  Min 0.001, Max 0.090]
 ];
 
+// 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.
+
 
 DefineConstant[
 	     // all dimensions in meter
@@ -60,3 +59,6 @@ DefineConstant[
   R=1        // refinement factor: R=5 far raw mesh, R=1 normal mesh
 ];
 
+By_target = 1.;
+Px_target = CoreX + A - E/2 ;
+Py_target = CoreY + B/2 ;

Ccore/shp.py

-# Open this file with Gmsh (interactively with File->Open->shape.py, or on the
-# command line with 'gmsh shape.py')
-#
+# Open this file with Gmsh
+# - interactively with File->Open->shp.py,
+# - in the command line with 'gmsh shape.py'
+
 
 import numpy as np
 import conveks
@@ -139,19 +140,18 @@ def test(n):
     return tab
 
 
-# x is the list of design parameters with format
+# List of design parameters with format
 # [ "Onelab name", value, Min, Max, "variable name" ]
 designParameters={
     0 : ['Optimization/Parameters/D', 0.01, 0.001, 0.050, 'D'],
     1 : ['Optimization/Parameters/E', 0.01 , 0.01, 0.15 , 'E']
-#    2 : ['Optimization/Parameters/F', 0.01, 0.001, 0.090, 'F']
 }
 
 
 x = {}
 index = int(c.getNumber('Optimization/Sweep on parameter'))-1
 
-if index >= 0:
+if index >= 0: # parameter sweep and optimisation over one design parameter
     x[0] = designParameters[index];
     test(int(c.getNumber('Optimization/Steps in range')))
     #callback=callbackF, 
@@ -160,7 +160,7 @@ if index >= 0:
                    options={'ftol': 1e-5, 'gtol': 1e-3, 'disp':True} )
     print res
 
-else:
+else: # optimisation over full parameter space
     for label, var in designParameters.iteritems():
         x[label] = var
     values = [designParameters[i][1] for i in range(len(designParameters))]