Evolutionary Topology Optimisation with GSA

Inspired by a recent Arup Thoughts article about a radical minimum weight joint, and making use of some online references (Huang and Xie, 2010), I have developed a Python routine for running an Evolutionary Topology Optimisation (ETO) on a GSA mesh.

The basic idea of this method is quite simple: remove the bits of the structure that are doing the least work, then repeat until everything is working hard. Hard work in this case is defined by the ratio of the highest stress, or strain energy, to the lowest (I chose strain energy in this case), and switch off any element that falls beneath the Rejection Ratio (RR). The trick is to ramp up the RR slowly so that the remaining structure is stable. You do this by starting the RR low, say at 1% of peak stress or 0.1% strain energy, and then increasing it by an small amount, the Evolutionary Rate (ER), whenever all elements pass the RR test. You then stop when RR reaches a set amount, say 25%.

In this Python script (see Using Python Scripts with MassMotion for notes on getting started) I make use of Analysis stages, so nothing is actually deleted from the model and you can see how the structure optimises at each step. On the other hand, I delete the analysis results after each loop to keep the file size down, so you will need to analyse all if you want to see how the stresses develop over time.

As an example, here is a solid mesh with two fixed support points that evolves into a cross between an arch and a truss.

The routine will end when either the RR reaches the max value or when the analysis fails. If the analysis fails because of bits of the mesh coming loose, try halving the evolutionary rate, as this can allow the mesh time to remove the weak areas.

Note that when you run the script you will need to change the directory to something suitable on your machine, as well as possibly the file names.

Do please share any other interesting meshes and their required ETO settings.This routine should work on either 2d mesh or 1d elements, as both give a strain energy output. If you want to change the script to work on stresses you will need to be more careful as 1D and 2D elements will have different output formats.

Reference: Huang, X. and Xie, Y. (2010). Evolutionary topology optimization of continuum structures. Chichester, West Sussex: Wiley.

The Python 3.5 Script

## GSA Evolutionary Topology Optimisation ##
import win32com.client
import pythoncom

# variables - edit the directory to suit !
dir = "C:\\Users\\peter.debney\\Documents\\Visual Studio 2013\\Projects\\PythonTopologyOptimisation\\"
input = "two simple supports.gwb"
output = "two simple supports out.gwb"
initial_rejection_ratio = 0.0001  # initial Rejection Ratio
ER = 0.00005  # Evolutionary Rate
max_rejection_ratio = 0.025
RR = initial_rejection_ratio # Rejection Ratio

# open GSA model
gsa_obj = win32com.client.Dispatch("Gsa_8_7.ComAuto")  # early binding
gsa_obj.VersionString()
gsa_obj.Open(dir + input)

stage = 1
element_list = []
nElem = int(gsa_obj.GwaCommand("HIGHEST,EL"))
for iElem in range(nElem + 1):  # add existing elements to the list
    if int(gsa_obj.GwaCommand("EXIST,EL," + str(iElem))) == 1:
        element_list.append(iElem)

while RR <= max_rejection_ratio:
    print("RR = " + str('% 0.5f' % RR) + " Stage = " + str(stage) +
         " No. elements = " + str(len(element_list)))

     # clear existing results to keep file size down

    gsa_obj.DELETE("Results")

     # create analysis stage from element list
    gsa_obj.GwaCommand("ANAL_STAGE," + str(stage) + ",RR " + str('% 0.5f' % RR) +
             ",NO_RGB," + " ".join(map(str, element_list)) + " ,0")
     # create analysis task with stage
    gsa_obj.GwaCommand("TASK," + str(stage) + ",Static," + str(stage) +
             ",GSS,SOL_STATIC,1,0,128,SELF,none,none,DRCMEFNSU,MIN," +
             "AUTO,0.000000,0.000000,0.000000,NONE,ERROR,NONE,NONE," +
             "RAFT_LO,RESID_NO,0.000000,1.00000")
    gsa_obj.GwaCommand("GSS_PARAM," + str(stage) + ",SOLVER_PARALLEL_DIRECT," +
             "PRE_LINE_JACOBI,1.00000e-016,1.00000e-010,0.00100000," +
             "1.00000e-010,1.50000,RESIDUAL,AUTO_GEOM,1.00000e-012," +
             "32,1.00000e-012,10000,1.00000e-010,4096,0.000500000," +
             "DISP,1.00000,0,NODE,INT_NORM,STURM,MITC,ALLMAN_COOK," +
             "CONDITION,EIGEN_SUBSPACE_SHIFT,0,SHIFT_NONE,NO_TIED")
    gsa_obj.GwaCommand("ANAL," + str(stage) + ",Analysis Case " + str(stage) +
             "," + str(stage) + ",L1")
    gsa_obj.Analyse(stage)
    gsa_obj.SaveAs(dir + output)

     # read the strain energy results
    results = {}
    max_strain = 0
     for iElem in element_list:
         if int(gsa_obj.GwaCommand("EXIST,EL," + str(iElem))) == 1:
             # extract strain energy for each element
            strain = gsa_obj.GwaCommand("GET,STRAIN_ENERGY," + str(iElem) + ","
                + str(stage))

     # convert strain string to list
    strain = strain.split(',')
    strain = float(strain[3])
         if strain is not []:  ####
             if strain > max_strain:
                max_strain = strain
        results.update({str(iElem): strain})

     # loop through results and list any that have strain energy >= RR
    new_list = []
     for iElem in results:
         if results[iElem] >= RR * max_strain:
            new_list.append(iElem)

     if len(new_list) == len(element_list):  # increase RR if no elements removed
        RR = RR + ER

    element_list = new_list  # update element list
    stage += 1

    gsa_obj.SaveAs(dir + output)

gsa_obj.Close()
gsa_obj = None
print("Finished")