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.