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