From f0cb67780234e637edf735f6597627290831cb2a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Fran=C3=A7ois=20Henrotte?= <francois.henrotte@uclouvain.be>
Date: Fri, 12 Oct 2018 17:07:39 +0200
Subject: [PATCH] tuto shape optimization in Magnetostatics
---
Ccore/ccore.geo | 167 +++++++++++++++
Ccore/ccore.pro | 457 +++++++++++++++++++++++++++++++++++++++++
Ccore/ccore_common.pro | 62 ++++++
Ccore/shp.plt | 53 +++++
Ccore/shp.py | 245 ++++++++++++++++++++++
5 files changed, 984 insertions(+)
create mode 100644 Ccore/ccore.geo
create mode 100644 Ccore/ccore.pro
create mode 100644 Ccore/ccore_common.pro
create mode 100644 Ccore/shp.plt
create mode 100644 Ccore/shp.py
diff --git a/Ccore/ccore.geo b/Ccore/ccore.geo
new file mode 100644
index 0000000..348dff3
--- /dev/null
+++ b/Ccore/ccore.geo
@@ -0,0 +1,167 @@
+
+/*
+ All geometrical parameters of the model are used in both
+ Gmsh (CAO) and GetDP (shape optimisation).
+ They are therefore defined in the file "ccore_common.pro".
+*/
+
+Include "ccore_common.pro";
+
+// The mesh refinement factor is rather a global user parameter,
+// therefore defined as a Onelab variable.
+R = DefineNumber[1, Name "Model/Parameters/Mesh refinement factor",
+ Help "R=5 far raw mesh, R=1 normal mesh"];
+
+// Center C-core in air box
+CoreX=(L-A)/2; // X-position of C-core's bottom left corner
+CoreY=(L-B)/2; // Y-position of C-core's bottom left corner
+
+// external boundary
+lc1=L/10.*R;
+
+Point(1) = {0,0,0,lc1};
+Point(2) = {L,0,0,lc1};
+Point(3) = {L,L,0,lc1};
+Point(4) = {0,L,0,lc1};
+Line(1) = {1,2};
+Line(2) = {2,3};
+Line(3) = {3,4};
+Line(4) = {4,1};
+Curve Loop(1) = {1 ... 4};
+
+// magnetic C-core
+
+lc2=A/10.*R;
+lc3=D/2.*R;
+
+Point(5) = {CoreX,CoreY,0,lc2};
+Point(6) = {CoreX+A,CoreY,0,lc2};
+Point(7) = {CoreX+A,CoreY+(B-D)/2.,0,lc3};
+Point(8) = {CoreX+A-E,CoreY+(B-D)/2.,0,lc3};
+Point(9) = {CoreX+A-E,CoreY+E,0,lc2};
+Point(10) = {CoreX+E,CoreY+E,0,lc2};
+Point(11) = {CoreX+E,CoreY+B-E,0,lc2};
+Point(12) = {CoreX+A-E,CoreY+B-E,0,lc2};
+Point(13) = {CoreX+A-E,CoreY+(B+D)/2.,0,lc3};
+Point(14) = {CoreX+A,CoreY+(B+D)/2.,0,lc3};
+Point(15) = {CoreX+A,CoreY+B,0,lc2};
+Point(16) = {CoreX,CoreY+B,0,lc2};
+
+Line(5) = {5,6};
+Line(6) = {6,7};
+Line(7) = {7,8};
+Line(8) = {8,9};
+Line(9) = {9,10};
+Line(10) = {10,11};
+Line(11) = {11,12};
+Line(12) = {12,13};
+Line(13) = {13,14};
+Line(14) = {14,15};
+Line(15) = {15,16};
+Line(16) = {16,5};
+
+Curve Loop(2) = {5 ... 16};
+
+// inductors
+
+lc4=lc2; //F/2.*R;
+
+Point(17) = {CoreX+E+F,CoreY+E+F,0,lc4};
+Point(18) = {CoreX+E+F+G,CoreY+E+F,0,lc4};
+Point(19) = {CoreX+E+F+G,CoreY+B-E-F,0,lc4};
+Point(20) = {CoreX+E+F,CoreY+B-E-F,0,lc4};
+Line(17) = {17,18};
+Line(18) = {18,19};
+Line(19) = {19,20};
+Line(20) = {20,17};
+
+Curve Loop(3) = {17 ... 20};
+
+Point(21) = {CoreX-F-G,CoreY+E+F,0,lc4};
+Point(22) = {CoreX-F,CoreY+E+F,0,lc4};
+Point(23) = {CoreX-F,CoreY+B-E-F,0,lc4};
+Point(24) = {CoreX-F-G,CoreY+B-E-F,0,lc4};
+Line(21) = {21,22};
+Line(22) = {22,23};
+Line(23) = {23,24};
+Line(24) = {24,21};
+
+Curve Loop(4) = {21 ... 24};
+
+
+Plane Surface(1) = {1,2,3,4};
+Plane Surface(2) = {2};
+Plane Surface(3) = {3};
+Plane Surface(4) = {4};
+
+
+Physical Surface("AIR", 1) = 1;
+Physical Surface("CORE", 2) = 2;
+Physical Surface("COILP", 3) = 3;
+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.
+ It should be copied unmodified at the bottom
+ of all shape optimisation .geo files.
+ */
+
+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
diff --git a/Ccore/ccore.pro b/Ccore/ccore.pro
new file mode 100644
index 0000000..e32cb7f
--- /dev/null
+++ b/Ccore/ccore.pro
@@ -0,0 +1,457 @@
+Include "ccore_common.pro";
+
+DefineConstant[
+ R_ = {"GetPerformances", Name "GetDP/1ResolutionChoices", Visible 0}
+ C_ = {"-solve -v 5 -v2", Name "GetDP/9ComputeCommand", Visible 0}
+ P_ = {"", Name "GetDP/2PostOperationChoices", Visible 0}
+ VelocityTag = -1
+ OptiIterNumber = 0
+];
+
+
+NL_tol_abs = 1e-8; // absolute tolerance on residual for noninear iterations
+NL_tol_relax = 1.0; // relaxation on residual for noninear iterations
+NL_iter_max = 50; // maximum number of noninear iterations
+Flag_NL =
+ DefineNumber[ 1, Name "Model/Parameters/Non linear core material", Choices {0,1}];
+murCore =
+ DefineNumber[300, Name "Model/Parameters/Mur core", Visible 1,//!Flag_NL,
+ Help "Magnetic relative permeability of Core"];
+Flag_Jfixed =
+ DefineNumber[ 1, Name "Model/Parameters/Fixed current density", Choices {0,1}];
+CurrentDensity = 1e6 *
+ DefineNumber[1, Name "Model/Parameters/Current Density", Visible Flag_Jfixed,
+ Help "Current density in coil [A/mm2]"];
+NbTurns = 300;
+Current =
+ DefineNumber[CurrentDensity*CoilSection_ref/NbTurns,
+ Name "Model/Parameters/Current [A]", Visible !Flag_Jfixed,
+ Help "Current in coil [A]"];
+
+CoilSection = G*(B-2*E-2*F);
+If(Flag_Jfixed)
+ Mmf = CurrentDensity * CoilSection;
+ Current = Mmf / NbTurns;
+Else
+ Mmf = NbTurns * Current;
+ CurrentDensity = Mmf / CoilSection;
+EndIf
+
+By_target = 1.;
+Px_target = CoreX + A - E/2 ;
+Py_target = CoreY + B/2 ;
+
+Group {
+ // Physical regions (in capital letters):
+ AIR = Region[ 1 ];
+ CORE = Region[ 2 ];
+ COILP = Region[ 3 ];
+ COILN = Region[ 4 ];
+ NOFLUX = Region[ 11 ];
+
+ // Abstract regions
+ Vol_Mag = Region[ {AIR, CORE, COILP, COILN} ];
+ Vol_S_Mag = Region[ {COILP, COILN} ];
+ Sur_Dir_Mag = Region[ {NOFLUX} ];
+ Sur_Neu_Mag = Region[ {} ];
+
+ Vol_NL_Mag = Region[ {} ];
+ If(Flag_NL)
+ Vol_NL_Mag = Region[ {CORE} ];
+ EndIf
+ Vol_L_Mag = Region[ {Vol_Mag,-Vol_NL_Mag} ];
+}
+
+Function {
+ mu0 = 4.e-7 * Pi ;
+ nu0 = 1. / mu0 ;
+
+ nu [ Region[{AIR, COILP, COILN}] ] = nu0;
+
+ Mat_h_r = { 0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 125, 150, 175, 200, 250,
+ 300, 400, 500, 600, 700, 800, 900, 1000, 1250, 1500, 2000, 2500, 5000,
+ 7500, 10000, 15000, 20000, 59000, 174000, 514000, 1520000, 4470000,
+ 13200000, 38900000, 115000000, 339000000, 1000000000 } ;
+ Mat_b_r = { 0.0, 0.194880829963, 0.377143018857, 0.537767739762, 0.672888260835,
+ 0.783043000477, 0.871342430831,0.941778611986, 0.998183303557, 1.04378111223,
+ 1.08110469369, 1.14963767549, 1.19607212343, 1.22964695907, 1.25515221835,
+ 1.29162498935, 1.31678879432, 1.35015120537, 1.37220092877, 1.38859114656,
+ 1.4017440574, 1.41287024565, 1.42264180514, 1.43146158921, 1.45082466146,
+ 1.46784549989, 1.49819370601, 1.52578650709, 1.64314027719, 1.73458485332,
+ 1.8039068939,1.89568786291, 1.95213815187, 2.1390774927, 2.45827909293,
+ 3.32303272825, 5.85485500678, 13.2701832298, 35.2114648741, 99.8027446541,
+ 291.062951228, 854.036370229, 2515.3105707 } ;
+
+ Mat_b2_r = Mat_b_r()^2;
+ Mat_nu_r = Mat_h_r()/Mat_b_r();
+ Mat_nu_r(0) = Mat_nu_r(1);
+ Mat_nu_b2_r = ListAlt[Mat_b2_r(), Mat_nu_r()];
+ nu_interp_r[] = InterpolationLinear[ SquNorm[$1] ]{Mat_nu_b2_r()};
+ dnudb2_interp_r[] = dInterpolationLinear[SquNorm[$1]]{Mat_nu_b2_r()};
+ h_interp_r[] = nu_interp_r[$1] * $1 ;
+
+ If(!Flag_NL)
+ nu [ CORE ] = nu0/murCore;
+ dhdb_NL[ CORE ] = 0;
+ mu_analytic = mu0*murCore;
+ Else
+ nu [ Vol_NL_Mag ] = nu_interp_r[$1] ;
+ dhdb_NL[ Vol_NL_Mag ] = 2*dnudb2_interp_r[$1#2]*SquDyadicProduct[#2];
+ mu_analytic = 1/Mat_nu_r(0);
+ EndIf
+
+ // shape function of J : js[] = nI * Jshape[]
+ Jshape[ COILP ] = Vector[0,0,+1]/CoilSection;
+ Jshape[ COILN ] = Vector[0,0,-1]/CoilSection;
+ js[] = Mmf * Jshape[];
+
+ // Analytical solution
+ RelGap = D/(mu0*E);
+ RelCore = (2*(A+B-2*E)-D)/(mu_analytic*E);
+ Flux = NbTurns*Mmf/(RelGap+RelCore); // analytic approximation
+ SetNumber("Model/Results/Flux analytic", Flux);
+ //Flux = 0;
+ If( ParamIndex == 0 ) // D
+ dFluxdTau = Flux/E*Flux/(NbTurns*Mmf)*(1/mu_analytic-1/mu0);
+ ElseIf( ParamIndex == 1 ) // E
+ dFluxdTau = Flux/E*(1.+4*Flux/(NbTurns*Mmf*mu_analytic));
+ If( Flag_Jfixed )
+ dFluxdTau -= 2*Flux*G/CoilSection;
+ EndIf
+ ElseIf( ParamIndex == 2 ) // F
+ dFluxdTau = 0;
+ EndIf
+
+
+ VV[] = Vector[$1,$2,0*$3] ;
+ dV[] = TensorV[$1,$2,0*$3] ;
+ Lie2form[] = TTrace[$2]*$1 - Transpose[$2]*$1; // $1=2-form (vector field), $2=dV (tensor)
+
+ If(Flag_Jfixed)
+ LieOf_js[] = Lie2form[js[], $1]; // $1=dV
+ Else
+ LieOf_js[] = 0;
+ EndIf
+
+ // 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;
+}
+
+
+
+Group {
+ Dom_Hcurl_a_Mag_2D = Region[ {Vol_Mag, Sur_Neu_Mag} ];
+}
+
+Function{
+ l_a = ListFromServer["Optimization/Results/a"];
+ aFromServer[] = ValueFromIndex[]{l_a()};
+ 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 {
+ { Name Dirichlet_a_Mag;
+ Case {
+ { Region NOFLUX; Type Assign; Value 0; }
+ }
+ }
+ { Name aFromServer;
+ Case {
+ { Region Vol_Mag; Type Assign; Value aFromServer[]; }
+ { Region NOFLUX; Type Assign; Value 0; }
+ }
+ }
+ For i In {1:3}
+ { Name velocity~{i} ;
+ Case {
+ { Region Vol_Mag ; Value velocity~{i}[]; }
+ }
+ }
+ EndFor
+}
+
+
+
+FunctionSpace {
+ { Name Hcurl_a_2D; Type Form1P;
+ BasisFunction {
+ { Name se; NameOfCoef ae; Function BF_PerpendicularEdge;
+ Support Dom_Hcurl_a_Mag_2D ; Entity NodesOf[ All ]; }
+ }
+ Constraint {
+ { NameOfCoef ae; EntityType NodesOf;
+ NameOfConstraint Dirichlet_a_Mag; }
+ }
+ }
+
+ { Name Hcurl_a_2D_fullyfixed; Type Form1P;
+ BasisFunction{
+ { Name se1; NameOfCoef ae1; Function BF_PerpendicularEdge;
+ Support Vol_Mag; Entity NodesOf[All]; }
+ }
+ Constraint{
+ { NameOfCoef ae1; EntityType NodesOf; NameOfConstraint aFromServer; }
+ }
+ }
+ For i In {1:3}
+ { Name H_v~{i} ; Type Form0;
+ BasisFunction {
+ { Name sn ; NameOfCoef un ; Function BF_Node ;
+ Support Vol_Mag; Entity NodesOf[ All ] ; }
+ }
+ Constraint {
+ { NameOfCoef un ; EntityType NodesOf ; NameOfConstraint velocity~{i}; }
+ }
+ }
+ EndFor
+}
+
+Jacobian {
+ { Name Vol ;
+ Case {
+ { Region All ; Jacobian Vol ; }
+ }
+ }
+ { Name Sur;
+ Case {
+ { Region All; Jacobian Sur; }
+ }
+ }
+}
+
+Integration {
+ { Name Int ;
+ Case { { Type Gauss ;
+ Case {
+ { GeoElement Point; NumberOfPoints 1; }
+ { GeoElement Line; NumberOfPoints 5; }
+ { GeoElement Triangle; NumberOfPoints 4; }
+ { GeoElement Quadrangle; NumberOfPoints 4; }
+ }
+ }
+ }
+ }
+}
+
+
+// -------------------------------------------------------------------------
+// This resolution solves the direct problem to compute the sensitivity
+// of the induction flux with respect to a given design variable.
+
+
+Formulation {
+ { Name MagSta_a; Type FemEquation;
+ Quantity {
+ { Name a; Type Local; NameOfSpace Hcurl_a_2D; }
+ }
+ Equation {
+ Integral { [ nu[{d a}] * Dof{d a} , {d a} ];
+ In Vol_Mag ; Jacobian Vol; Integration Int; }
+ Integral { JacNL [ dhdb_NL[{d a}] * Dof{d a} , {d a} ];
+ In Vol_NL_Mag; Jacobian Vol; Integration Int; }
+ Integral { [ -js[] , {a} ];
+ In Vol_S_Mag; Jacobian Vol; Integration Int; }
+ }
+ }
+}
+
+Resolution {
+ { Name GetPerformances;
+ System {
+ { Name SYS; NameOfFormulation MagSta_a;}
+ }
+ Operation {
+ If( OptiIterNumber == 1 )
+ CreateDir[ResDir];
+ DeleteFile[StrCat[ResDir,"w.txt"]];
+ DeleteFile[StrCat[ResDir,"Lie_w.txt"]];
+ EndIf
+ InitSolution[SYS];
+ IterativeLoop[NL_iter_max, NL_tol_abs, NL_tol_relax]{
+ GenerateJac[SYS];
+ SolveJac[SYS];
+ }
+ PostOperation[Get_ObjectiveConstraints];
+ }
+ }
+}
+
+PostProcessing {
+ { Name MagSta_a_2D; NameOfFormulation MagSta_a;
+ Quantity {
+ { Name az;
+ Value {
+ Term { [ CompZ[{a}] ]; In Dom_Hcurl_a_Mag_2D; Jacobian Vol; }
+ }
+ }
+ { Name b;
+ Value {
+ Term { [ {d a} ]; In Dom_Hcurl_a_Mag_2D; Jacobian Vol; }
+ }
+ }
+ { Name by;
+ Value {
+ Term { [ CompY[{d a}] ]; In Dom_Hcurl_a_Mag_2D; Jacobian Vol; }
+ }
+ }
+ { Name js;
+ Value {
+ Term { [ js[] ]; In Vol_S_Mag; Jacobian Vol; }
+ }
+ }
+ { Name Flux ; Value { Integral {
+ [ {a} * NbTurns * Jshape[] ] ;
+ In Vol_S_Mag ; Jacobian Vol ; Integration Int ; } }
+ }
+ { Name fobj_b_at_point;
+ Value {
+ Term { [ SquNorm[ CompY[{d a}] - By_target ] ];
+ In Dom_Hcurl_a_Mag_2D; Jacobian Vol; }
+ }
+ }
+ }
+ }
+}
+
+PostOperation {
+ { Name Get_ObjectiveConstraints; NameOfPostProcessing MagSta_a_2D;
+ Operation{
+
+ CreateDir["res"];
+ //Print[js, OnElementsOf Vol_S_Mag, File "res/js.pos"];
+ Print[b, OnElementsOf Vol_Mag, File "res/b.pos"];
+
+ Print[az, OnElementsOf Vol_Mag, File "res/az.pos"];
+ Echo[ Str["k= PostProcessing.NbViews-1;",
+ "View[k].IntervalsType = 3;",
+ "View[0].NbIso = 30;"],
+ File "tmp.geo", LastTimeStepOnly] ;
+
+ // In memory transfer of the a field for assembling sensitivity equations
+ Print[az, OnElementsOf Vol_Mag, Format NodeTable, File "",
+ SendToServer "Optimization/Results/a", Hidden 1];
+
+ Print[ Flux[ Vol_S_Mag ], OnGlobal,
+ Format Table, File > "",
+ SendToServer "Model/Results/Flux computed"];
+
+ Print[ by, OnPoint { Px_target, Py_target, 0},
+ Format Table, File "",
+ SendToServer "Model/Results/b probe"];
+
+ Print[ fobj_b_at_point, OnPoint { Px_target, Py_target, 0},
+ Format Table, File > "res/w.txt",
+ SendToServer "Optimization/Results/w"];
+ }
+ }
+}
+
+// -------------------------------------------------------------------------
+// This resolution solves the direct problem to compute the sensitivity
+// of the induction flux with respect to a given design variable.
+
+Formulation {
+ { Name DirectFormulationOf_MagSta_a_2D; Type FemEquation ;
+ Quantity {
+ { Name Lva; Type Local; NameOfSpace Hcurl_a_2D; }
+ { Name a; Type Local; NameOfSpace Hcurl_a_2D_fullyfixed; }
+ 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 Vol_Mag; Jacobian Vol; Integration Int; }
+ Galerkin { [ dhdb_NL[{d a}] * Dof{d Lva}, {d Lva} ];
+ In Vol_NL_Mag; Jacobian Vol; Integration Int; }
+ Galerkin { [LieOf_HB[{d a},dV[{d v_1},{d v_2},{d v_3}]], {d Lva}];
+ In Vol_Mag; Jacobian Vol; Integration Int; }
+ Galerkin { [LieOf_HB_NL[{d a},dV[{d v_1},{d v_2},{d v_3}]], {d Lva}];
+ In Vol_NL_Mag; Jacobian Vol; Integration Int; }
+
+ Integral { [ -LieOf_js[dV[{d v_1},{d v_2},{d v_3}]], {Lva} ] ;
+ In Vol_S_Mag ; Jacobian Vol; Integration Int; }
+
+ Integral { [ 0*Dof{a}, {a} ] ;
+ In Vol_Mag; Jacobian Vol ; Integration Int ; }
+ For i In {1:3}
+ Integral { [ 0*Dof{v~{i}}, {v~{i}} ] ;
+ In Vol_Mag; Jacobian Vol ; Integration Int ; }
+ EndFor
+ }
+ }
+}
+
+Resolution {
+ { Name GetGradient_wrt_dv;
+ System {
+ { Name DIR; NameOfFormulation DirectFormulationOf_MagSta_a_2D; }
+ }
+ Operation {
+ InitSolution[DIR];
+ Generate[DIR];
+ Solve[DIR];
+ PostOperation[Get_GradOf_w];
+ }
+ }
+}
+
+PostProcessing{
+ { Name Direct_MagSta; NameOfFormulation DirectFormulationOf_MagSta_a_2D;
+ PostQuantity {
+ { Name az;
+ Value{
+ Term{ [ CompZ[{a}] ]; In Vol_Mag; Jacobian Vol; }
+ }
+ }
+ { Name VV;
+ Value{
+ Term{ [ Vector[{v_1},{v_2},{v_3}] ]; In Vol_Mag ; Jacobian Vol;}
+ }
+ }
+ { Name Lie_az;
+ Value {
+ Term { [ CompZ[ {Lva}] ] ; In Vol_Mag ; Jacobian Vol; }
+ }
+ }
+ { Name Lie_b_at_point ; Value {
+ Term { [ 2 * (CompY[{d a}] - By_target) * CompY[{d Lva}] ] ;
+ In Vol_Mag; Jacobian Vol; }
+ }
+ }
+ }
+ }
+}
+
+PostOperation{
+ { Name Get_GradOf_w; NameOfPostProcessing Direct_MagSta;
+ Operation {
+ CreateDir["res"];
+
+ Print[ VV, OnElementsOf Vol_Mag,File "res/velocity.pos" ];
+ Echo[ Str["View[PostProcessing.NbViews-1].GlyphLocation = 2;"],
+ File "tmp.geo", LastTimeStepOnly] ;
+
+ //Print[az, OnElementsOf Vol_Mag, File "res/azstar.pos"];
+ Echo[ Str["k= PostProcessing.NbViews-1;",
+ "View[k].IntervalsType = 3;",
+ "View[0].NbIso = 30;"],
+ File "tmp.geo", LastTimeStepOnly] ;
+
+ Print[ Lie_az, OnElementsOf Vol_Mag, File "res/Lie_az.pos" ];
+ Echo[ Str["k= PostProcessing.NbViews-1;",
+ "View[k].IntervalsType = 3;",
+ "View[0].NbIso = 30;"],
+ File "tmp.geo", LastTimeStepOnly] ;
+
+ Print[ Lie_b_at_point, OnPoint{ Px_target, Py_target, 0},
+ Format Table, File > "res/Lie_w.txt",
+ SendToServer Sprintf["Optimization/Results/dwdtau_%g",VelocityTag]];
+
+ }
+ }
+}
diff --git a/Ccore/ccore_common.pro b/Ccore/ccore_common.pro
new file mode 100644
index 0000000..3898872
--- /dev/null
+++ b/Ccore/ccore_common.pro
@@ -0,0 +1,62 @@
+
+/*
+ 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
+
+ */
+
+
+ParamIndex =
+ DefineNumber[0, Name "Optimization/Sweep on parameter",
+ Choices {0="None",
+ 1="D: air gap length",
+ 2="E: C-core thickness",
+ 3="F: Core-coil gap"} ];
+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"];
+
+
+
+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]
+];
+
+
+DefineConstant[
+ // all dimensions in meter
+ L=1 // sidelength of the air box square
+
+ CoreX=0.3 // X-position of C-core's bottom left corner
+ CoreY=0.3 // Y-position of C-core's bottom left corner
+ A=0.4 // C-core length along X-axis
+ B=0.4 // C-core length along Y-axis
+ D=0.01 // air gap length
+ E=0.1 // C-core thickness
+ F=0.01 // gap between C-core and coil
+ G=0.05 // coil width along X-axis
+
+ CoilSection_ref = G*(B-2*E-2*F)
+
+ R=1 // refinement factor: R=5 far raw mesh, R=1 normal mesh
+];
+
diff --git a/Ccore/shp.plt b/Ccore/shp.plt
new file mode 100644
index 0000000..7fee2e8
--- /dev/null
+++ b/Ccore/shp.plt
@@ -0,0 +1,53 @@
+#set terminal pdf
+
+set grid
+set ytics nomirror
+set y2tics
+set style data lines
+
+set xlabel "tau"
+set y2label "w"
+set ylabel "dw/dtau"
+
+set output "aaa.pdf"
+plot "aaa" u 2:3 axis x1y2 t"w",\
+ "aaa" u 2:4 t"dwdt",\
+ "aaa" u 2:5 w linesp t"dwdt FD"
+
+pause -1
+
+set output "optishape1.pdf"
+plot "Torque_linTH.txt" u 2:3 axis x1y2 t"w",\
+ "Torque_linTH.txt" u 2:4 t"dwdt",\
+ "Torque_linTH.txt" u 2:5 w linesp t"dwdt FD"
+
+pause -1
+
+set output "optishape2.pdf"
+plot "LieDevSmoother.txt" u 2:3 axis x1y2 t"w",\
+ "LieDevSmoother.txt" u 2:4 t"dwdt",\
+ "LieDevSmoother.txt" u 2:5 w linesp t"dwdt FD"
+
+pause -1
+
+set output "optishape3.pdf"
+plot "Torque_linTH2.txt" u 2:3 axis x1y2 t"w",\
+ "Torque_linTH2.txt" u 2:4 t"dwdt",\
+ "Torque_linTH2.txt" u 2:5 w linesp t"dwdt FD"
+
+pause -1
+
+set output "optishape4.pdf"
+plot "LieDevSmoother2.txt" u 2:3 axis x1y2 t"w",\
+ "LieDevSmoother2.txt" u 2:4 t"dwdt",\
+ "LieDevSmoother2.txt" u 2:5 w linesp t"dwdt FD"
+
+pause -1
+
+set output "optishape5.pdf"
+plot "Torque_TH.txt" u 2:3 axis x1y2 t"w",\
+ "Torque_TH.txt" u 2:4 t"dwdt",\
+ "Torque_TH.txt" u 2:5 w linesp t"dwdt FD"
+
+pause -1
+
diff --git a/Ccore/shp.py b/Ccore/shp.py
new file mode 100644
index 0000000..32aa8d2
--- /dev/null
+++ b/Ccore/shp.py
@@ -0,0 +1,245 @@
+# Open this file with Gmsh (interactively with File->Open->shape.py, or on the
+# command line with 'gmsh shape.py')
+#
+
+import numpy as np
+import conveks
+import onelab
+from scipy.optimize import minimize
+
+
+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))
+
+modelName = 'ccore'
+
+# append correct path to model file names
+file_geo = c.getPath(modelName+'.geo')
+file_msh = c.getPath(modelName+'.msh')
+file_pro = c.getPath(modelName+'.pro')
+
+# build command strings with appropriate parameters and options
+mygmsh += ' -p 0 -v 2 -parametric ' + file_geo + ' '
+mygetdp += ' -p 0 ' + file_pro + ' -msh ' + file_msh + ' '
+
+# load geometry in GUI to have something to look at
+c.openProject(file_geo)
+
+c.sendCommand('Solver.AutoMesh = -1;')
+
+# dry getdp run (without -solve or -pos option) to get model parameters in the GUI
+c.runSubClient('myGetDP', mygetdp)
+
+# define optimization parameters as Onelab parameter (editable in the GUI)
+maxIter = c.defineNumber('Optimization/00Max iterations', value=100)
+maxChange = c.defineNumber('Optimization/01Max change', value=1e-0)
+
+# end of check phase. We're done if we do not start the optimization
+if c.action == 'check' :
+ exit(0)
+
+## 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-7):
+ for id, var in xvar.iteritems():
+ c.runNonBlockingSubClient('myGmsh', mygmsh \
+ + ' -setnumber PerturbMesh 1'\
+ + ' -setnumber Perturbation '+str(perturb)\
+ + ' -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])
+
+def fobj(xFromOpti):
+ # send the current point to GetDP model
+ for label, var in x.iteritems():
+ x[label][1] = xFromOpti[label]
+ c.setNumber(var[0],value=x[label][1])
+
+ # reload the geometry in GUI to see the geometrical changes
+ c.openProject(file_geo)
+
+ c.runSubClient('myGmsh', mygmsh + ' -2')
+ getVelocityField(x)
+ c.runSubClient('myGetDP', mygetdp + '-solve GetPerformances' \
+ + ' - setnumber OuterIteration ' + str(0)) # it !!
+ return c.getNumber('Optimization/Results/w')
+
+
+def grad(xFromOpti):
+ # send the current point to GetDP model
+ for label, var in x.iteritems():
+ x[label][1] = xFromOpti[label]
+ c.setNumber(var[0],value=x[label][1])
+
+ # 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()
+ grads = [c.getNumber('Optimization/Results/dwdtau_0'), c.getNumber('Optimization/Results/dwdtau_1')]
+ return np.asarray(grads)
+
+Nfeval = 1
+def callbackF(Xi):
+ global Nfeval
+ print '{0:4d} {1: 3.6f} {2: 3.6f} {3: 3.6f}'.format(Nfeval, Xi[0], fobj(Xi), grad(Xi)[0])
+ Nfeval += 1
+
+def test(n):
+ tab = np.zeros((n+1,4))
+ for step in range(n+1):
+ xstep = x[0][2] + (x[0][3]-x[0][2])/float(n)*step
+ f,g,d = fobj([xstep]), grad([xstep]), 0
+ tab[step] = [xstep, f, g, d]
+ if step >= 2 :
+ tab[step-1][3] = (tab[step][1]-tab[step-2][1])/(tab[step][0]-tab[step-2][0])
+ if c.getString(c.name+'/Action') == 'stop':
+ exit(1)
+ print
+ for step in range(n+1):
+ ind1, ind2 = max(step-1,0), min(step+1,n)
+ tab[step][3] = (tab[ind2][1]-tab[ind1][1])/(tab[ind2][0]-tab[ind1][0])
+ print "%4d" % step,
+ for item in tab[step] : print "%3.6f" % item,
+ print
+
+ return tab
+
+
+# x is the 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:
+ x[0] = designParameters[index];
+ test(int(c.getNumber('Optimization/Steps in range')))
+ #callback=callbackF,
+ res = minimize(fobj, [x[0][1]], jac=grad, bounds=[(x[0][2],x[0][3])],\
+ callback=callbackF, method = 'L-BFGS-B', tol=None,\
+ options={'ftol': 1e-5, 'gtol': 1e-3, 'disp':True} )
+ print res
+
+else:
+ for label, var in designParameters.iteritems():
+ x[label] = var
+ values = [designParameters[i][1] for i in range(len(designParameters))]
+ bounds = [(designParameters[i][2], designParameters[i][3]) for i in range(len(designParameters))]
+ ## print x
+ ## print values, bounds
+ ## print fobj(values)
+ ## print grad(values)
+ ## exit(1)
+
+ res = minimize(fobj, values, jac=grad, bounds=bounds, callback=callbackF,\
+ method = 'L-BFGS-B', tol=None, options={'ftol': 1e-5, 'gtol': 1e-3, 'disp':True} )
+ print res
+
+
+
+exit(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
+conveks.mma.initialize(initialPoint, lowerBound, upperBound)
+
+# Set some options for MMA
+conveks.mma.option.setNumber('General.Verbosity', 0) # def 4
+conveks.mma.option.setNumber('General.SaveHistory', 0)
+conveks.mma.option.setNumber('SubProblem.isRobustMoveLimit', 1) # def 1
+conveks.mma.option.setNumber('SubProblem.isRobustAsymptotes', 1) # def 1
+conveks.mma.option.setNumber('SubProblem.type', 0)
+conveks.mma.option.setNumber('SubProblem.addConvexity', 0)
+conveks.mma.option.setNumber('SubProblem.asymptotesRmax', 100.0)
+conveks.mma.option.setNumber('SubProblem.asymptotesRmin', 0.001) # def 0.001
+
+# Get iteration count (here it will be 1 - could be different in case of restart)
+it = conveks.mma.getOuterIteration()
+change = 1.0
+
+while it <= maxIter and c.getString(c.name+'/Action') != 'stop': #and change > 1e-6 :
+ c.setNumber('Optimization/01Current iteration', value=it, readOnly=1,
+ attributes={'Highlight':'LightYellow'})
+
+ # get (copy of) current point
+ xFromMMA = np.array(conveks.mma.getCurrentPoint())
+ print 'xFromMMA = ', xFromMMA
+
+ # get the value of the objective function and of the constraints
+ # as well as their sensitivity with respect to design variables at `x'
+ objective = fobj(xFromMMA)
+ grad_objective = grad(xFromMMA)
+ constraints = np.array([np.sum(xFromMMA)/100.0-1.0])
+ grad_constraints = np.ones(numVariables)/100.0
+
+ if it == 1: fscale = 1.0 / np.abs(objective)
+ objective *= fscale
+ grad_objective *= fscale
+
+ print it, xFromMMA, objective, grad_objective
+ c.sendInfo('Optimization: it. {}, obj. {}, constr. {}'.format(it,objective,constraints[0]))
+
+ # call MMA and update the current point
+ conveks.mma.setMoveLimits(lowerBound, upperBound, 1e-0) # parametre interessant
+ conveks.mma.updateCurrentPoint(constraints,grad_objective,grad_constraints)
+ change = conveks.mma.getDesignChange()
+ it = conveks.mma.getOuterIteration()
+
+# This should be called at the end
+conveks.mma.finalize()
+
--
GitLab