2.7.2: A Simplified ESO-Procedure on Shells

import clr
clr.AddReferenceToFileAndPath("C:\Program Files\Rhino 8\Plug-ins\karamba\Karamba.gha")
clr.AddReferenceToFileAndPath("C:\Program Files\Rhino 8\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 Karamba.Utilities.Utils as Utils
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)
loadCaseCombinationIndex = 0
Utils.handleError(response.update(loadCaseCombinationIndex), analysis);
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_meshLast updated