I finally got to write my own FEA code. Always wanted to do so but dreadful because I’m easily distracted and need a good time to focus on this.
Without a background in engineering yet working in an engineering company, I hope writing FEA can help me understand the engineering math.
I mainly follow this article right here and translate it to python. I will eventually upload my code to git but for now I will copy paste the pieces of my code here.
https://podgorskiy.com/spblog/304/writing-a-fem-solver-in-less-the-180-lines-of-code
I use the platform I’m most familiar with, Rhino Grasshopper which have existing methods to define and use mesh. However, I have to use Rhino 5 32bit this time because I want to use numpy for matrix calculation. I have covered the installation in a post below
Inputs
To start a FEA we need boundaries conditions:
A good triangulated mesh
Nodes which are mesh vertices
Support Conditions (nodes that are restrained)
Loading Conditions (loading vectors, loadPts)
Material Info (Poisson Ratio and Young Modulus)
Below is my set up in grasshopper
I position the mesh at origin
Set Up
In order to calculate efficiently, we need to do some “prepping”. We will set up list of nodes, edges, list of loads, list of constraints and list of elements (an element is a triangulated face)
Nodes: The input nodes currently are point3d. However, for convenience, we often refer to the nodes using its index. Whenever the situation arise, I use index to call out the nodes.
Edges: I store the indexes of the line’s start point and end point in a value pair
def GetEdges(mesh):
edges = [] #list of vertex indices that make up the edges
for i in range(mesh.TopologyEdges.Count):
vPair = mesh.TopologyEdges.GetTopologyVertices(i)
vPairPython = (vPair.I, vPair.J)
edges.append(vPairPython)
return edges
Loads: For loads, I need to store them in an array twice the size of the nodes array. Alternate items of the array are X and Y coordinates of the load vector
def GetLoads(nodes, loadPts, loadVs):
loads = [0] * (len(nodes) * 2)
for i in range(len(nodes)):
for j in range(len(loadPts)):
loadV = loadVs[j]
dist = rs.Distance(loadPts[j],nodes[i])
if dist < 0.001:
loads[2 * i + 0] = loadV.X
loads[2 * i + 1] = loadV.Y
return loads
Constraints: I define a class for constraint where type is indicated by integer (1 is for constraint, x, 2 is for y, 3 is for both), node is point3d.
class Constraint:
def __init__(self, type, node):
self.Type = type
self.Node = node
def GetConstraints(supportBoundaries, supportConds, nodes):
constraints = []
for i in range(len(nodes)):
for j in range(len(supportBoundaries)):
if rs.PointInPlanarClosedCurve(nodes[i], supportBoundaries[j]) == 1:
con = Constraint(supportConds[j], i)
constraints.append(con)
return constraints
Elements: I created Element object woth properties are Nodes and NodeIds
class Element:
def __init__(self, elemNodes, nodeIds):
self.Nodes = elemNodes
self.NodeIds = nodeIds
def GetElements(mesh):
faceCt = rs.MeshFaceCount(mesh)
elements = []
for i in range(faceCt):
faceVerticeIDs = rs.MeshFaceVertices(mesh)[i]
nodeIds = faceVerticeIDs[:-1]
elemNodes = [nodes[nodeIds[0]], nodes[nodeIds[1]], nodes[nodeIds[2]]]
element = Element(elemNodes, nodeIds)
elements.append(element)
return elements
To list out all of them we have:
edges = GetEdges(mesh) loads = GetLoads(nodes, loadPts, loadVs) constraints = GetConstraints(supportBoundaries,supportConds, nodes ) #(restraintType, id) elements = GetElements(mesh)
And…you’re done setting up the items to start on doing some math. I will cover the calculation/matrix manipulations in the next post
BitByDesign 