2.7.2: A Simplified ESO-Procedure on Shells

The example file “SimpleShellEso.gh” contains also a Python script which performs a simple evolutionary structural optimization (ESO) procedure on shell elements. As it applies no filters for calculating the fitness of individual shell triangles checkerboard patterns result (see fig. 2.7.2.1). Alas the script can be easily extended to include more elaborate fitness calculation schemes.

Fig. 2.7.2.1: SimpleShellEso.gh

The script shows how to work directly with the C++ model in order to avoid costly mappings to and from the C#-model:

import clr

clr.AddReferenceToFileAndPath("C:\Program Files\Rhino 6\Plug-ins\Karamba\Karamba.gha")
clr.AddReferenceToFileAndPath("C:\Program Files\Rhino 6\Plug-ins\Karamba\KarambaCommon.dll")

import Karamba.Models.Model as Model
import Karamba.Elements.ModelShell as Shell
import Karamba.Materials.FemMaterial_Isotrop as FemMaterial
import feb.ShellMesh as ShellMesh
import feb.TriShell3D as TriShell3D
import feb.VectSurface3DSigEps as TriStates
import feb.Deform as Deform
import feb.Response as Response
import feb.EnergyVisitor as EnergyVisitor
import Rhino.Geometry as Rh

from operator import attrgetter

# encapsulate ESO properties of shell elements
class EsoItem:
    def __init__(self, shell_elem, elem_ind):
       self.active = True
       self.fitness = 0
       self.shell_elem = shell_elem
       self.area = shell_elem.area()
       self.ind = elem_ind

    def update(self, energy_visitor):
       self.fitness = energy_visitor.elasticEnergy(self.ind) / self.area
        

# clone model to avoid side effects
model = Model_in.Clone()
model.deepCloneFEModel()

# generate ESO properties of each triangular shell element
eso_items = []
for elem in model.elems:
    if type(elem) != Shell:
       continue
    tri_mesh = model.febmodel.triMesh(elem.fe_id)
    for i in xrange(tri_mesh.numberOfElems()):
        eso_items.append(EsoItem(tri_mesh.elem(i), i))

nremove_per_iter = int(NRemove/NIter+1)
n_removed = 0

# do the ESO iterations
for iter in xrange(NIter):
    analysis = Deform(model.febmodel)
    response = Response(analysis)

    try:
        response.updateNodalDisplacements()
        response.updateMemberForces()
    except:
        raise Exception("singular stiffness")

    energy_visitor = EnergyVisitor(model.febmodel, model.febmodel.state(0), 0);
    energy_visitor.visit(model.febmodel);
   
    for eso_item in eso_items:
        eso_item.update(energy_visitor)

    eso_items = sorted(eso_items, key = attrgetter("fitness"))

    n_removed_per_iter = 0
    has_changed = False
    for eso_item in eso_items:
        if (n_removed >= NRemove): break
        if (n_removed_per_iter >= nremove_per_iter): break
        if (eso_item.active == False):
            continue
        eso_item.shell_elem.softKilled(True)
        eso_item.active = False
        n_removed +=1
        n_removed_per_iter +=1
     
    has_changed = True
    if (has_changed == False):
        break
    model.febmodel.touch()

    # create active and inactive mesh for output
    active_mesh = Rh.Mesh()
    inactive_mesh = Rh.Mesh()
    for i in xrange(model.febmodel.numberOfNodes()):
        feb_pos = model.febmodel.node(i).pos()
        active_mesh.Vertices.Add(Rh.Point3d(feb_pos.x(), feb_pos.y(), feb_pos.z()))
        inactive_mesh.Vertices.Add(Rh.Point3d(feb_pos.x(), feb_pos.y(), feb_pos.z()))
        
    for eso_item in eso_items:
        ind0 = eso_item.shell_elem.node(0).ind()
        ind1 = eso_item.shell_elem.node(1).ind()
        ind2 = eso_item.shell_elem.node(2).ind()
        if (eso_item.active):
            active_mesh.Faces.AddFace(Rh.MeshFace(ind0, ind1, ind2))
        else:
            inactive_mesh.Faces.AddFace(Rh.MeshFace(ind0, ind1, ind2))
                
    activeMesh = active_mesh
    inactiveMesh = inactive_mesh

In the above code the class “ESOItem” handles the book-keeping necessary in the optimization steps. it contains the activation-state, the fitness, the element’s area, a reference to the C++-element and the element’s index in the C#-model. The “update”-method calculates the specific elastic energy of the underlying shell-element.

Activation and deactivation of model elements works via setting the soft-kill status of C++-elements to “True” or “False” (see line 74). On model-assembly the stiffness of the corresponding element will be multiplied with the soft-kill-factor which is 1.0×10^−10. This factor can be set on the C++-model via “softKillFactor(new_factor)” if necessary.

The last part of the script categorizes the shell-faces into active or in-active adding their geometry to the corresponding output-meshes.