shapeOptimization.py
You can view and download this file on Github: shapeOptimization.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Example for simple shape optimization with NGsolve/Netgen
5#
6# Author: Johannes Gerstmayr
7# Date: 2023-05-23
8#
9# Copyright:This file is part of Exudyn. Exudyn is free software. You can redistribute it and/or modify it under the terms of the Exudyn license. See 'LICENSE.txt' for more details.
10#
11#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
12
13import exudyn as exu
14from exudyn.utilities import *
15from exudyn.FEM import *
16from exudyn.graphicsDataUtilities import *
17from exudyn.processing import GeneticOptimization, ParameterVariation, PlotOptimizationResults2D, Minimize
18
19import numpy as np
20import time
21#import timeit
22
23import exudyn.basicUtilities as eb
24import exudyn.rigidBodyUtilities as rb
25import exudyn.utilities as eu
26
27import numpy as np
28
29useGraphics = True
30
31
32
33#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
34#netgen/meshing part:
35fem = FEMinterface()
36#standard:
37L = 1 #Length of beam
38wy = 0.25*L
39#wzNominal = 0.25*L
40wzNominal = 0.25*L
41
42meshSize = wzNominal*0.25 #*0.25
43nModes = 8
44
45rho = 1000 #polymer or similar
46Emodulus=5e9
47nu=0.3
48meshCreated = False
49verbose = True
50
51import ngsolve as ngs
52import netgen
53from netgen.meshing import *
54
55from netgen.geom2d import unit_square
56#import netgen.libngpy as libng
57from netgen.csg import *
58
59
60#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
61#if True: #needs netgen/ngsolve to be installed with pip install; to compute mesh, see e.g.: https://github.com/NGSolve/ngsolve/releases
62# rangeMax = 4
63# for param in range(rangeMax):
64def ParameterFunction(parameterSet):
65 SC = exu.SystemContainer()
66 mbs = SC.AddSystem()
67
68 try:
69 #++++++++++++++++++++++++++++++++++++++++++++++
70 #++++++++++++++++++++++++++++++++++++++++++++++
71 #store default parameters in structure (all these parameters can be varied!)
72 class P: pass #create emtpy structure for parameters; simplifies way to update parameters
73 P.x0 = 0.4*L
74 P.x1 = 0.6*L
75 P.x2 = 0.85*L
76 P.r0 = wy*0.2
77 P.r1 = wy*0.2
78 P.r2 = wy*0.2
79
80 P.computationIndex = 'Ref'
81
82 # #now update parameters with parameterSet (will work with any parameters in structure P)
83 for key,value in parameterSet.items():
84 setattr(P,key,value)
85
86 verbose = P.computationIndex == 'Ref'
87
88 geo = CSGeometry()
89
90 Vblock = L*wy*wzNominal
91 Vcyls = P.r0**2*pi*wzNominal + P.r1**2*pi*wzNominal + P.r2**2*pi*wzNominal
92 wz = wzNominal*Vblock/(Vblock-Vcyls) #increase thickness
93
94 block = OrthoBrick(Pnt(0,-0.5*wy,-0.5*wz),Pnt(L,0.5*wy,0.5*wz))
95 cyl0 = Cylinder(Pnt(P.x0,0,-0.5*wz),Pnt(P.x0,0,0.5*wz),P.r0)
96 cyl1 = Cylinder(Pnt(P.x1,0,-0.5*wz),Pnt(P.x1,0,0.5*wz),P.r1)
97 cyl2 = Cylinder(Pnt(P.x2,0,-0.5*wz),Pnt(P.x2,0,0.5*wz),P.r2)
98
99 part = block - (cyl0+cyl1+cyl2)
100 geo.Add(part)
101
102 #Draw (geo)
103
104 mesh = ngs.Mesh( geo.GenerateMesh(maxh=meshSize, curvaturesafety=1.25)) #smaller values give coarser mesh
105 mesh.Curve(2)
106
107 if False: #set this to true, if you want to visualize the mesh inside netgen/ngsolve
108 # import netgen
109 #+++++++++++++++++++
110 import netgen.gui
111 ngs.Draw (mesh)
112 for i in range(40):
113 netgen.Redraw() #this makes the window interactive
114 time.sleep(0.05)
115
116 meshCreated = True
117 #+++++++++++++++++++++++++++++++++++++++++++++++++++++
118 #Use fem to import FEM model and create FFRFreducedOrder object
119 fem.ImportMeshFromNGsolve(mesh, density=rho, youngsModulus=Emodulus, poissonsRatio=nu, meshOrder=2)
120
121 # print("nNodes=",fem.NumberOfNodes())
122 #check total mass:
123 if verbose:
124 print("nNodes=",fem.NumberOfNodes())
125 nNodes = fem.NumberOfNodes()
126 M = CSRtoScipySparseCSR(fem.GetMassMatrix(sparse=True))
127 Phit = np.kron(np.ones(nNodes),np.eye(3)).T
128 totalMass = (Phit.T @ M @ Phit)[0,0]
129
130 print("total mass=",totalMass)
131
132 #+++++++++++++++++++++++++++++++++++++++++++++++++++++
133 #compute Hurty-Craig-Bampton modes
134
135 if True:
136 pLeft = [0,-0.5*wy,-0.5*wz]
137 pRight = [L,-0.5*wy,-0.5*wz]
138 nTip = fem.GetNodeAtPoint(pRight) #tip node (do not use midpoint, as this may not be a mesh node ...)
139 #print("nMid=",nMid)
140 nodesLeftPlane = fem.GetNodesInPlane(pLeft, [-1,0,0])
141 lenNodesLeftPlane = len(nodesLeftPlane)
142 weightsLeftPlane = np.array((1./lenNodesLeftPlane)*np.ones(lenNodesLeftPlane))
143
144 nodesRightPlane = fem.GetNodesInPlane(pRight, [-1,0,0])
145 lenNodesRightPlane = len(nodesRightPlane)
146 weightsRightPlane = np.array((1./lenNodesRightPlane)*np.ones(lenNodesRightPlane))
147
148 #boundaryList = [nodesLeftPlane, nodesRightPlane] #second boudary (right plane) not needed ...
149 boundaryList = [nodesLeftPlane]
150
151 # if verbose:
152 # print("compute HCB modes... ")
153
154 start_time = time.time()
155 fem.ComputeHurtyCraigBamptonModes(boundaryNodesList=boundaryList,
156 nEigenModes=nModes,
157 useSparseSolver=True,
158 computationMode = HCBstaticModeSelection.RBE2,
159 verboseMode=False,
160 timerTreshold=100000)
161
162 fMin = min(fem.GetEigenFrequenciesHz())
163 if verbose:
164 # print("HCB modes needed %.3f seconds" % (time.time() - start_time))
165 print('smallest Eigenfrequency=',fMin)
166
167
168 #alternatives:
169 #fem.ComputeEigenModesWithBoundaryNodes(boundaryNodes=nodesLeftPlane, nEigenModes=nModes, useSparseSolver=False)
170 #fem.ComputeEigenmodes(nModes, excludeRigidBodyModes = 6, useSparseSolver = True)
171 #print("eigen freq.=", fem.GetEigenFrequenciesHz())
172
173
174 #+++++++++++++++++++++++++++++++++++++++++++++++++++++
175 #animate modes
176 if verbose:
177 print("create CMS element ...")
178 cms = ObjectFFRFreducedOrderInterface(fem)
179
180 objFFRF = cms.AddObjectFFRFreducedOrder(mbs, positionRef=[0,0,0],
181 initialVelocity=[0,0,0],
182 initialAngularVelocity=[0,0,0],
183 gravity=[0,-9.81,0],
184 color=[0.1,0.9,0.1,1.],
185 )
186
187 from exudyn.interactive import AnimateModes
188 mbs.Assemble()
189 SC.visualizationSettings.nodes.show = False
190 SC.visualizationSettings.openGL.showFaceEdges = True
191 SC.visualizationSettings.openGL.multiSampling=4
192
193 #+++++++++++++++++++++++++++++++++++++++
194 #animate modes of ObjectFFRFreducedOrder (only needs generic node containing modal coordinates)
195 SC.visualizationSettings.general.autoFitScene = False #otherwise, model may be difficult to be moved
196
197 nodeNumber = objFFRF['nGenericODE2'] #this is the node with the generalized coordinates
198 AnimateModes(SC, mbs, nodeNumber, runOnStart= True)
199 # import sys
200 # sys.exit()
201
202 #do not show values larger than -120 to have nicer plots ...
203 return 140 + min(-100,-fMin) #maximize eigenfrequency => minimize ...
204 except:
205 return 140-99
206
207
208
209#%%
210#now perform parameter variation
211if __name__ == '__main__': #include this to enable parallel processing
212 import time
213
214 # refval = ParameterFunction({'computationIndex':'Ref'}) # compute reference solution
215 # import sys
216 # sys.exit()
217 #%%++++++++++++++++++++++++++++++++++++++++++++++++++++
218 #GeneticOptimization
219 start_time = time.time()
220 rMax = 0.5*wy
221 rRange = (0.3*rMax, 0.95*rMax)
222
223 x0Range = (1*rMax, 3*rMax)
224 x1Range = (4.5*rMax, 5*rMax)
225 x2Range = (6*rMax, 7*rMax)
226
227
228 if False:
229 #this option does not work here, as it would require the ParameterFunction to restrict the ranges to feasible values
230 [pOpt, vOpt, pList, values] = Minimize(
231 objectiveFunction = ParameterFunction,
232 parameters = {'r0':rRange, 'r1':rRange, 'r2':rRange,
233 'x0':x0Range, 'x1':x1Range, 'x2':x2Range, }, #parameters provide search range
234 showProgress=True,
235 enforceBounds=True,
236 addComputationIndex = True,
237 options={'maxiter':100},
238 resultsFile='solution/shapeOptimization.txt'
239 )
240 else:
241 [pOpt, vOpt, pList, values] = GeneticOptimization(
242 objectiveFunction = ParameterFunction,
243 parameters = {'r0':rRange, 'r1':rRange, 'r2':rRange,
244 'x0':x0Range, 'x1':x1Range, 'x2':x2Range, }, #parameters provide search range
245 numberOfGenerations = 10,
246 populationSize = 100,
247 elitistRatio = 0.1,
248 # crossoverProbability = 0, #makes no sense!
249 rangeReductionFactor = 0.7,
250 randomizerInitialization=0, #for reproducible results
251 #distanceFactor = 0.1,
252 distanceFactor = 0.,
253 debugMode=False,
254 useMultiProcessing=True, #may be problematic for test
255 numberOfThreads = 20,
256 addComputationIndex=True,
257 showProgress=True,
258 resultsFile = 'solution/shapeOptimization.txt',
259 )
260
261 exu.Print("--- %s seconds ---" % (time.time() - start_time))
262
263 exu.Print("[pOpt, vOpt]=", [pOpt, vOpt])
264 u = vOpt
265 exu.Print("optimized eigenfrequency=",-u)
266
267 if True:
268 # from mpl_toolkits.mplot3d import Axes3D # noqa: F401 unused import
269 import matplotlib.pyplot as plt
270
271 plt.close('all')
272 [figList, axList] = PlotOptimizationResults2D(pList, values, yLogScale=False)
273
274 refval = ParameterFunction({**pOpt,'computationIndex':'Ref'} ) # compute reference solution
275
276#i9, 14 cores, numberOfThreads = 20: 575.85 seconds (ngen=10, pposize=100)
277#with distanceFactor: 117.16 (ngen=12, popsize=200)
278#without distanceFactor: 117.5287 (ngen=10, popsize=100)
279
280#pOpt = {'r0': 0.04012197560525166, 'r1': 0.08694009640170029, 'r2': 0.11845643292562649, 'x0': 0.2828892936420842, 'x1': 0.6170637302367472, 'x2': 0.8749125935436654}
281#optimized eigenfrequency= 117.52874179749946
282#p0 = {'r0': 0.0625, 'r1': 0.0625, 'r2': 0.0625, 'x0': 0.25, 'x1': 0.59375, 'x2': 0.8125}