topology optimization of L-bracket using the pnorm of von-mises

Lbracket/topo.geo

+lc=0.004;
+mm=1e-3;
+
+Point(1) = {0, 0, 0, lc};
+Point(2) = {600*mm, 0, 0, lc};
+Point(3) = {600*mm, 240*mm, 0, lc};
+Point(4) = {240*mm, 240*mm, 0, lc};
+Point(5) = {240*mm, 600*mm, 0, lc};
+Point(6) = {0, 600*mm, 0, lc};
+Point(10) = {(600-10)*mm, 240*mm, 0, lc};
+Point(11) = {(600-10)*mm, (240-30)*mm, 0, lc};
+Point(12) = {600*mm, (240-30)*mm, 0, lc};
+
+Line(1) = {6, 1};
+Line(2) = {1, 2};
+Line(3) = {2, 12};
+Line(4) = {12, 3};
+Line(5) = {10, 3};
+Line(6) = {10, 4};
+Line(7) = {4, 5};
+Line(8) = {5, 6};
+Line(9) = {12, 11};
+Line(10) = {11, 10};
+
+Line Loop(1) = {1, 2, 3, 9, 10, 6, 7, 8};
+Plane Surface(1) = {1};
+
+Line Loop(2) = -{9, 10, 5, -4};
+Plane Surface(2) = {2};
+
+Physical Surface(1000) = {1};
+Physical Surface(1004) = {2};
+Physical Line(1001) = {8};
+Physical Line(1003) = {4};
+Physical Point(1002) = {3};
+
+Solver.AutoMesh = 0;
+Mesh.Algorithm = 1; 
+

Lbracket/topo.py

+# Open this file with Gmsh (interactively with File->Open->topo.py, or on the
+# command line with 'gmsh topo.py')
+#
+
+import numpy as np
+import optlab
+import onelab
+
+c = onelab.client(__file__)
+
+# get paths to the gmsh and getdp executable 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))
+
+# append correct path to model file names
+topo_geo = c.getPath('topo.geo')
+topo_msh = c.getPath('topo.msh')
+topo_pro = c.getPath('topo.pro')
+
+# build command strings with appropriate parameters and options
+mygetdp += ' -p 0 ' + topo_pro + ' -msh ' + topo_msh + ' '
+mygmsh += ' -p 0 -v 2 ' + topo_geo + ' '
+
+# load geometry in GUI to have something to look at
+c.openProject(topo_geo)
+
+# dry getdp run (without -solve or -pos option) to get model parameters in the GUI
+c.runSubClient('myGetDP', mygetdp)
+
+# define now optimization parameters
+# some of them as Onelab parameter, to be editable in the GUI
+maxIter = c.defineNumber('Optimization/00Max iterations', value=1600)
+maxChange = c.defineNumber('Optimization/01Max change', value=1e-2)
+filtering = c.defineNumber('Optimization/02Filtering/00Filter densities', value=1, choices=[0,1])
+filterContinuation = c.defineNumber('Optimization/02Filtering/01continuation?', value=0, choices=[0,1])
+nbLoopReduceFilterRadius = c.defineNumber('Optimization/02Filtering/02frequency of radius update', value=40)
+filterRadiusDecrease = c.defineNumber('Optimization/02Filtering/03reduce factor', value=0.7)
+filterRadiusMin = c.defineNumber('Optimization/02Filtering/04radius min', value=0.0015)
+filterRadius = c.getNumber('Optimization/3Density/2filter radius')
+
+# end of check phase. We're done if we do not start the optimization
+if c.action == 'check' :
+   exit(0)
+
+# mesh the geometry to create the density (background) mesh
+c.runSubClient('myGmsh', mygmsh + '-2')
+
+# initializes the density 'xe' field
+# This is done with a GetDP problem named 'DensityField'
+# that solves no FE equations
+# but only fills up the function space with an initial value. 
+c.runSubClient('myGetDP', mygetdp + '-solve DensityField')
+
+# read in the 'xe' field from GetDP into a python variable 'd'
+d = c.getNumbers('Optimization/Results/designVariable')
+
+# the variable d contains the following serialized data: [N, ele1, xe1, ..., eleN, xeN]
+# split it into homogeneous vectors
+numElements = d[0] # N
+elements = np.array(d[1::2]) # [ele1, ..., eleN ]
+#x = np.array(d[2::2]) # [xe1, ..., xeN]
+x = 0.5*np.ones(numElements)
+
+print 'number of elements in topology optimization:', numElements
+
+# This function serializes the python variables 'ele' and 'var' into
+# an ElementTable variable [N, ele1, xe1, ..., eleN, xeN]
+# and set it as a Onelab parameter named `var' 
+def setElementTable(var, ele, val):
+   data = np.ravel(np.column_stack((ele, val)))
+   data = np.insert(data, 0, numElements)
+   c.setNumber(var, values=data, visible=0)
+
+# The following function evaluates the optimization problem
+# at point `point' in the design space, i.e.
+# - the objective function,
+# - the constraints,
+# - the derivatives of the objective function w.r.t. the design variables,
+# - the derivatives of the constraints w.r.t. the design variables.
+def myOptimizationProblem(point=None):
+    #objective = c.getNumber('Optimization/Results/Energy')
+    objective = c.getNumber('Optimization/Results/pnorm of Von-Mises')
+    tmp = [c.getNumber('Optimization/Results/VolumeDomain')\
+      / (0.5*c.getNumber('Optimization/Results/VolumeDomainFull')) - 1.0]
+    constraints = np.array(tmp)
+
+    #grad_objective = np.array(c.getNumbers('Optimization/Results/GradOfCompliance')[2::2])
+    grad_objective = np.array(c.getNumbers('Optimization/Results/GradOfVonMises')[2::2])
+
+    tmp = list(np.asarray(c.getNumbers('Optimization/Results/GradOfVolume')[2::2])\
+         /(0.5*c.getNumber('Optimization/Results/VolumeDomainFull')))
+    grad_constraints = np.array(tmp)
+
+    return objective,constraints,grad_objective,grad_constraints
+
+# number of design variables and number of constraints
+numVariables = len(x)
+
+# lower bound for design variables (here all set to 0.001)
+lowerBound = 0.001*np.ones(numVariables)
+
+# upper bound for design variables (here all set to 1)
+upperBound = np.ones(numVariables)
+
+# Initialize the MMA optimizer
+optlab.mma.initialize(x, lowerBound, upperBound)
+
+# Set appropriate options for MMA
+optlab.mma.option.setNumber('General.Verbosity', 4)
+optlab.mma.option.setNumber('General.SaveHistory', 0)
+optlab.mma.option.setNumber('SubProblem.isRobustMoveLimit', 1)
+optlab.mma.option.setNumber('SubProblem.isRobustAsymptotes', 1)
+optlab.mma.option.setNumber('SubProblem.type', 0)
+optlab.mma.option.setNumber('SubProblem.addConvexity', 1)
+optlab.mma.option.setNumber('SubProblem.asymptotesRmax', 100.0)
+optlab.mma.option.setNumber('SubProblem.asymptotesRmin', 0.001)
+
+# Get iteration count (here it will be 1 - could be different in case of restart)
+it = optlab.mma.getOuterIteration()
+change = 1.0
+
+while change > maxChange and it <= maxIter and c.getString('topo/Action') != 'stop':
+
+    c.setNumber('Optimization/01Current iteration', value=it, readOnly=1,
+                attributes={'Highlight':'LightYellow'})
+    c.setNumber('Optimization/02Current change', value=change, readOnly=1,
+                attributes={'Highlight':'LightYellow'})
+
+    # get (copy of) current point
+    x = np.array(optlab.mma.getCurrentPoint())
+
+    if filtering:
+        setElementTable('Optimization/Results/filterInput', elements, x)
+        c.runSubClient('myGetDP', mygetdp + '-solve FilterDensity')
+        x = np.array(c.getNumbers('Optimization/Results/filterOutput')[2::2])
+
+    setElementTable('Optimization/Results/designVariable', elements, x)
+
+    # compute objective function and constraints as well as their sensitivity
+    # with respect to design variables at `x'
+    c.runSubClient('myGetDP', mygetdp + '-solve ComputePerformancesAndGradient')
+
+    # get the value of the objective function and of the constraints
+    # as well as their sensitivity with respect to design variables at `x'
+    objective,constraints,grad_objective,grad_constraints = myOptimizationProblem()
+    objective = c.getNumber('Optimization/Results/pnorm of Von-Mises')
+    constraints = np.array([c.getNumber('Optimization/Results/VolumeDomain')\
+      / (0.5*c.getNumber('Optimization/Results/VolumeDomainFull')) - 1.0])
+    grad_objective = np.array(c.getNumbers('Optimization/Results/GradOfVonMises')[2::2])
+    grad_constraints = np.asarray(c.getNumbers('Optimization/Results/GradOfVolume')[2::2])\
+         /(0.5*c.getNumber('Optimization/Results/VolumeDomainFull'))
+
+    if filtering:
+        # apply the chain rule for the sensitivity of objective
+        setElementTable('Optimization/Results/filterInput', elements, grad_objective)
+        c.runSubClient('myGetDP', mygetdp + '-solve FilterDensity')
+        grad_objective = np.array(c.getNumbers('Optimization/Results/filterOutput')[2::2])
+
+        # apply the chain rule for the sensitivity of constraints
+        d_con_dx = grad_constraints.reshape((len(constraints), numVariables))
+        for k,dfe in enumerate(d_con_dx):
+            setElementTable('Optimization/Results/filterInput', elements, dfe)
+            c.runSubClient('myGetDP', mygetdp + '-solve FilterDensity')
+            d_con_dx[k] = np.array(c.getNumbers('Optimization/Results/filterOutput')[2::2])
+        grad_constraints = d_con_dx.reshape((len(constraints) * numVariables,))
+
+    if it == 1: fscale = 1.0 / np.abs(objective)
+    objective *= fscale
+    grad_objective *= fscale
+
+    print '*'*50
+    print 'iteration: ', it
+    print 'change: ', change
+    print 'objective:', objective
+    print 'constraints:', constraints
+    c.sendInfo('Optimization: it. {}, obj. {}, constr. {}, change {}'.format(it,objective,constraints[0], change))
+
+    # call MMA and update the current point
+    optlab.mma.setMoveLimits(lowerBound, upperBound, 0.1)
+    optlab.mma.updateCurrentPoint(constraints, grad_objective, grad_constraints)
+    change = optlab.mma.getDesignChange()
+    it = optlab.mma.getOuterIteration()
+
+    if filtering and filterContinuation and \
+        (it % nbLoopReduceFilterRadius == 0.0 or change <= maxChange): 
+        filterRadius *= filterRadiusDecrease
+        print "decrease filter radius: ", filterRadius
+        c.setNumber('Optimization/3Density/2filter radius',value=filterRadius)
+        filtering = filterRadius <= filterRadiusMin
+        print "filtering: ",filtering
+        c.setNumber('Optimization/02Filtering/00Filter densities',value=filtering)
+        change = 1.0
+
+# This should be called at the end
+optlab.mma.finalize()