Edi Liberato blog

This post is about some finite element method code that I wrote back in 2014. I am sharing it here, even though I never finished the data visualization part I wanted back then.

I won’t write about theory after all there’s a million books out there about this subject. I’m just sharing my take on hands on code, something that you don’t usually find in books. The code is written in python.

At the time I was heavily inspired by Sukhbinder Singh. Go check his blog if you didn’t yet, tons of interesting stuff there. I will also leave a big thanks! for him here, in case he ever happens to read this.

Trusses

2D Truss

Let’s start simple with a 2d application of the method. This example is taken from the following book:

Hutton, D. V. (2004). Fundamentals of finite element analysis. 1st ed. Boston: McGraw-Hill. Page 69.

The two-element truss in Figure is subjected to external loading as shown. Using the same node and element numbering as in Figure, determine the displacement components of node 3, the reaction force components at nodes 1 and 2, and the element displacements, stresses, and forces. The elements have modulus of elasticity E1=E2=10x10^6 lb/in^2 and cross-sectional areas A1=A2=1.5 in^2.

truss_2d

##############################################################
# Fem script to compute plane truss displacements and stresses 
# written by Eddie Liberato 16/11/16
##############################################################

import numpy as np

# constants
E=10e06
A=1.5

# structure geometry 
node_coord=np.array([[0,0],[0,40],[40,40]]) # nodes coordinates
elem_con=np.array([[0,2],[1,2]]) # elements connectivity

x=node_coord[:,0] # sliced array with only x node coordinates
y=node_coord[:,1] # sliced array with only y node coordinates
node_count=len(node_coord)
elem_count=len(elem_con)
struct_dof=2*node_count # (entire) structure degrees of freedom

# matrices initialization
displacement=np.zeros((struct_dof,1))
force=np.zeros((struct_dof,1))
sigma=np.zeros((elem_count,1))
stiffness=np.zeros((struct_dof,struct_dof))

# load assignments
force[4]=500
force[5]=300
 
# computations
for e in range(elem_count):
	index=elem_con[e]
	elem_dof=np.array([index[0]*2, index[0]*2+1, index[1]*2, index[1]*2+1])
	xl=x[index[1]]-x[index[0]]
	yl=y[index[1]]-y[index[0]]
	elem_length=np.sqrt(xl*xl+yl*yl)
	c=xl/elem_length
	s=yl/elem_length
	rot=np.array([[c*c, c*s, -c*c, -c*s],
                  [c*s, s*s, -c*s, -s*s],
                  [-c*c, -c*s, c*c, c*s],
                  [-c*s, -s*s, c*s, s*s]])
	k=(E*A/elem_length)*rot
	stiffness[np.ix_(elem_dof, elem_dof)] +=k

suppress_dof=np.array([0,1,2,3]) # constrained degrees of freedom
active_dof=np.setdiff1d(np.arange(struct_dof), suppress_dof)
displacement_aux=np.linalg.solve(stiffness[np.ix_(active_dof,active_dof)], force[np.ix_(active_dof)])
displacement[np.ix_(active_dof)]=displacement_aux
react=np.dot(stiffness, displacement)

for e in range(elem_count):
	index=elem_con[e]
	elem_dof=np.array([index[0]*2, index[0]*2+1, index[1]*2, index[1]*2+1])
	xl=x[index[1]]-x[index[0]]
	yl=y[index[1]]-y[index[0]]
	elem_length=np.sqrt(xl*xl+yl*yl)
	c=xl/elem_length
	s=yl/elem_length
	sigma[e]=(E/elem_length)*np.dot(np.array([-c,-s,c,s]), displacement[np.ix_(elem_dof)])

# emitting results to screen
print(f'displacements:\n {displacement}')
print(f'stress:\n {sigma}')

The output:

$ python truss2d.py
displacements:
 [[0.        ]
 [0.        ]
 [0.        ]
 [0.        ]
 [0.00053333]
 [0.00172941]]
stress:
 [[282.84271247]
 [133.33333333]]

All right, pretty cool. Short, dense script. A lot is being done in just a few lines of code. But even cooler is the fact that you can relatively easy expand the method (and the script) to solve 3d trusses, also called space trusses.

Space Truss

I’m taking the 3d example from the book of Darryl Logan, a book that I recommend as the first book to be read by someone studying the Finite Element Method.

Logan D. L. (2007). A First Course in the Finite Element Method. 4th ed. Thomson. Page 98.

Analyze the space truss shown in Figure. The truss is composed of four nodes, whose coordinates (in meters) are shown in the figure, and three elements, whose cross-sectional areas are all 10x10^-4 m^2. The modulus of elasticity E=210 GPa for all the elements. A load of 20 kN is applied at node 1 in the global x-direction. Nodes 2–4 are pin supported and thus constrained from movement in the x, y, and z directions.

truss_3d

###############################################################
# Fem script to compute space truss displacements and stresses 
# written by Eddie Liberato 16/11/16
###############################################################

import numpy as np

# constants
E=210e09
A=10e-04

# structure geometry
node_coord=np.array([[12,-3,-4],[0,0,0],[12,-3,-7],[14,6,0]]) # node cordenates
elem_con=np.array([[0,1],[0,2],[0,3]]) # elements connectivity

x=node_coord[:,0] # sliced array with only x node coord
y=node_coord[:,1] # sliced array with only y node coord
z=node_coord[:,2] # sliced array with only z node coord
node_count=len(node_coord)
elem_count=len(elem_con)
struct_dof=3*node_count # (entire) structure degrees of freedom

# matrices initialization
displacement=np.zeros((struct_dof,1))
force=np.zeros((struct_dof,1))
sigma=np.zeros((elem_count,1))
stiffness=np.zeros((struct_dof,struct_dof))
np.set_printoptions(precision=3)

# load assignments
force[0]=20e03

# computations
for e in range(elem_count):
	index=elem_con[e]
	elem_dof=np.array([index[0]*3, index[0]*3+1, index[0]*3+2,
                       index[1]*3, index[1]*3+1, index[1]*3+2])
	xl=x[index[1]]-x[index[0]]
	yl=y[index[1]]-y[index[0]]
	zl=z[index[1]]-z[index[0]]
	elem_length=np.sqrt(xl*xl+yl*yl+zl*zl)
	cx=xl/elem_length
	cy=yl/elem_length
	cz=zl/elem_length
	trans=np.array([[cx*cx,cx*cy,cx*cz,-cx*cx,-cx*cy,-cx*cz],
                    [cy*cx,cy*cy,cy*cz,-cy*cx,-cy*cy,-cy*cz],
                    [cz*cx,cz*cy,cz*cz,-cz*cx,-cz*cy,-cz*cz],
                    [-cx*cx,-cx*cy,-cx*cz,cx*cx,cx*cy,cx*cz],
                    [-cy*cx,-cy*cy,-cy*cz,cy*cx,cy*cy,cy*cz],
                    [-cz*cx,-cz*cy,-cz*cz,cz*cx,cz*cy,cz*cz]])
	k=(E*A/elem_length)*trans
	stiffness[np.ix_(elem_dof, elem_dof)] +=k

suppress_dof=np.array([3,4,5,6,7,8,9,10,11]) # constrained degrees of freedom
active_dof=np.setdiff1d(np.arange(struct_dof), suppress_dof)
displacement_aux=np.linalg.solve(stiffness[np.ix_(active_dof,active_dof)], force[np.ix_(active_dof)])
displacement[np.ix_(active_dof)]=displacement_aux
react=np.dot(stiffness, displacement)

for e in range(elem_count):
	index=elem_con[e]
	elem_dof=np.array([index[0]*3, index[0]*3+1, index[0]*3+2,
                       index[1]*3, index[1]*3+1, index[1]*3+2])
	xl=x[index[1]]-x[index[0]]
	yl=y[index[1]]-y[index[0]]
	zl=z[index[1]]-z[index[0]]
	elem_length=np.sqrt(xl*xl+yl*yl+zl*zl)
	cx=xl/elem_length
	cy=yl/elem_length
	cz=zl/elem_length
	sigma[e]=(E/elem_length)*np.dot(np.array([-cx,-cy,-cz,cx,cy,cz]), displacement[np.ix_(elem_dof)])

# emitting results to screen
print(f'reactions:\n {react.reshape(node_count,3)}')
print(f'displacements:\n {displacement}')
print(f'stress:\n {sigma}')

The output:

$ python truss3d.py
reactions:
 [[ 2.000e+04 -1.137e-13 -9.095e-13]
 [-1.895e+04  4.737e+03  6.316e+03]
 [ 0.000e+00  0.000e+00 -4.211e+03]
 [-1.053e+03 -4.737e+03 -2.105e+03]]
displacements:
 [[ 1.384e-03]
 [-5.157e-05]
 [ 6.015e-05]
 [ 0.000e+00]
 [ 0.000e+00]
 [ 0.000e+00]
 [ 0.000e+00]
 [ 0.000e+00]
 [ 0.000e+00]
 [ 0.000e+00]
 [ 0.000e+00]
 [ 0.000e+00]]
stress:
 [[20526315.789]
 [ 4210526.316]
 [-5289408.222]]

Frames

Portal Frame

A brief tale about how this started: My young self trying to design a bamboo bicycle frame. Unfortunately the bike itself never passed the drawing board stage, but in the process I learned a lot about mechanical structures and designing things using unusual materials. Here I share part of the effort, more specifically, the efforts to calculate the forces on members of the frame.

As you’ll see the examples are not about bike frames, but simpler structures. This is due to the fact that I don’t want to (re)write about theory here. I just want to share hands on code. The examples are taken from a book where you can find more detailed explanation about conceptual part of the method.

Logan D. L. (2007). A First Course in the Finite Element Method. 4th ed. Thomson. Page 218.

As the first example of rigid plane frame analysis, solve the simple ‘‘bent’’. The frame is fixed at nodes 1 and 4 and subjected to a positive horizontal force of 10,000 lb applied at node 2 and to a positive moment of 5000 lb-in. applied at node 3. The global-coordinate axes and the element lengths are shown in Figure. Let E=30x10^6 psi and A=10 in^2 for all elements, and let I=200 in^4 for elements 1 and 3, and I=100^4 for element 2.

frame_2d

###############################################################
# Fem script to compute plane frame forces and moments
# written by Eddie Liberato 16/11/16
###############################################################

import numpy as np

# constants
E=30e06         # psi
A=10            # in^2
I=[200,100,200] # in^4

# structure geometry
node_coord=np.array([[0,0],[0,10*12],[10*12,10*12],[10*12,0]]) # node coordinates
elem_con=np.array([[0,1],[1,2],[2,3]]) # elements connectivity

x=node_coord[:,0] # sliced array with only x node coord
y=node_coord[:,1] # sliced array with only y node coord
node_count=len(node_coord)
elem_count=len(elem_con)
struct_dof=3*node_count # (entire) structure degrees of freedom

# matrices initialization
displacement=np.zeros((struct_dof,1))
force=np.zeros((struct_dof,1))
stiffness=np.zeros((struct_dof,struct_dof))

# load assignments
force[3]=10000
force[8]=5000

# computations
for e in range(elem_count):
        index=elem_con[e]
        elem_dof=np.array([index[0]*3, index[0]*3+1, index[0]*3+2,
                           index[1]*3, index[1]*3+1, index[1]*3+2])
        xl=x[index[1]]-x[index[0]]
        yl=y[index[1]]-y[index[0]]
        el=np.sqrt(xl*xl+yl*yl)
        c=xl/el
        s=yl/el
        kglobal=(E/el)*np.array([[A*c**2+12*I[e]/el**2*s**2,(A-12*I[e]/el**2)*c*s,-6*I[e]/el*s,-(A*c**2+12*I[e]/el**2*s**2),-(A-12*I[e]/el**2)*c*s,-(6*I[e]/el)*s],
                                 [(A-12*I[e]/el**2)*c*s, A*s**2+12*I[e]/el**2*c**2, 6*I[e]/el*c, -(A-12*I[e]/el**2)*c*s, -(A*s**2+12*I[e]/el**2*c**2), (6*I[e]/el)*c],
                                 [-6*I[e]/el*s, 6*I[e]/el*c, 4*I[e], 6*I[e]/el*s, -(6*I[e]/el)*c, 2*I[e]],
                                 [-(A*c**2+12*I[e]/el**2*s**2),-(A-12*I[e]/el**2)*c*s, 6*I[e]/el*s, A*c**2+12*I[e]/el**2*s**2, (A-12*I[e]/el**2)*c*s, 6*I[e]/el*s],
                                 [-(A-12*I[e]/el**2)*c*s, -(A*s**2+12*I[e]/el**2*c**2), -(6*I[e]/el)*c, (A-12*I[e]/el**2)*c*s, A*s**2+12*I[e]/el**2*c**2, -(6*I[e]/el)*c],
                                 [-(6*I[e]/el)*s, (6*I[e]/el)*c, 2*I[e], 6*I[e]/el*s, -(6*I[e]/el)*c, 4*I[e]]])
        stiffness[np.ix_(elem_dof, elem_dof)] +=kglobal

suppress_dof=np.array([0,1,2,9,10,11])
active_dof=np.setdiff1d(np.arange(struct_dof), suppress_dof)
displacement_aux=np.linalg.solve(stiffness[np.ix_(active_dof,active_dof)], force[np.ix_(active_dof)])
displacement[np.ix_(active_dof)]=displacement_aux
print(f'displacements:\n {displacement}')

for e in range(elem_count):
	print(f'forces on element {e+1}')
	index=elem_con[e]
	elem_dof=np.array([index[0]*3, index[0]*3+1, index[0]*3+2,
                       index[1]*3, index[1]*3+1, index[1]*3+2])
	xl=x[index[1]]-x[index[0]]
	yl=y[index[1]]-y[index[0]]
	el=np.sqrt(xl*xl+yl*yl)
	c=xl/el
	s=yl/el
	trans=np.array([[c,s,0,0,0,0],[-s,c,0,0,0,0],[0,0,1,0,0,0],[0,0,0,c,s,0],[0,0,0,-s,c,0],[0,0,0,0,0,1]])
	Td=np.dot(trans, displacement[np.ix_(elem_dof)])
	c1=(A*E/el)
	c2=(E*I[e]/el**3)
	klocal=np.array([[c1,0,0,-c1,0,0],
                     [0,12*c2,6*c2*el,0,-12*c2,6*c2*el],[0,6*c2*el,4*c2*el**2,0,-6*c2*el,2*c2*el**2],
                     [-c1,0,0,c1,0,0],
                     [0,-12*c2,-6*c2*el,0,12*c2,-6*c2*el],
                     [0,6*c2*el,2*c2*el**2,0,-6*c2*el,4*c2*el**2]])
	elem_forces=np.dot(klocal,Td)
	print(f'{elem_forces}')

The output:

$ python frame2d.py
displacements:
 [[ 0.        ]
 [ 0.        ]
 [ 0.        ]
 [ 0.21136266]
 [ 0.00148133]
 [-0.00152603]
 [ 0.20935933]
 [-0.00148133]
 [-0.001486  ]
 [ 0.        ]
 [ 0.        ]
 [ 0.        ]]
forces on element 1
[[ -3703.31950207]
 [  4991.69435216]
 [375803.32156987]
 [  3703.31950207]
 [ -4991.69435216]
 [223200.00068927]]
forces on element 2
[[   5008.30564784]
 [  -3703.31950207]
 [-223200.00068927]
 [  -5008.30564784]
 [   3703.31950207]
 [-221198.3395597 ]]
forces on element 3
[[  3703.31950207]
 [  5008.30564784]
 [226198.3395597 ]
 [ -3703.31950207]
 [ -5008.30564784]
 [374798.33818117]]

3D Frame

Finally the script I used in practice when I was designing the bike, the 3d frame. Things get interesting and relatively complex here. You’ll need to care about directions and torsional parameters. Not for the faint of heart.

There’s a bit of bloat there because I failed the don’t repeat yourself mantra. I said to myself I’d come back later and rewrite that part but ended up never doing it.

Logan D. L. (2007). A First Course in the Finite Element Method. 4th ed. Thomson. Page 262.

Determine the displacements and rotations at the free node (node 1) and the element local forces and moments for the space frame shown in Figure. Also verify equilibrium at node 1. Let E=30 000 Ksi, G=10 000 Ksi, J=50 in^4, Iy=100 in^4, Iz=100 in^4, A=10 in^2, and L=100 inches for all three beam elements.

frame_3d

###############################################################
# Fem script to compute space frame forces and moments
# written by Eddie Liberato 16/11/16
###############################################################

import numpy as np

# constants
E=30e03     # psi
A=10        # in^2
Iy=100
Iz=100
Gs=10000
Jt=50

# structure geometry
node_coord=np.array([[0,0,0],[-100,0,0],[0,0,-100],[0,-100,0]]) # node coordinates
elem_con=np.array([[1,0],[2,0],[3,0]]) # elements connectivity
xx=node_coord[:,0] # sliced array with only x node coord
yy=node_coord[:,1] # sliced array with only y node coord
zz=node_coord[:,2] # sliced array with only z node coord
node_count=len(node_coord)
elem_count=len(elem_con)
struct_dof=6*node_count # (entire) structure degrees of freedom

# matrices initialization
displacement=np.zeros((struct_dof,1))
force=np.zeros((struct_dof,1))
stiffness=np.zeros((struct_dof,struct_dof))
np.set_printoptions(precision=3)

# load assignments
force[1]=-50
force[3]=-1000

# computations
for e in range(elem_count):
	index=elem_con[e]
	elem_dof=np.array([index[0]*6, index[0]*6+1, index[0]*6+2,
                       index[0]*6+3, index[0]*6+4, index[0]*6+5,
                       index[1]*6, index[1]*6+1, index[1]*6+2,
                       index[1]*6+3, index[1]*6+4, index[1]*6+5])
	xl=xx[index[1]]-xx[index[0]]
	yl=yy[index[1]]-yy[index[0]]
	zl=zz[index[1]]-zz[index[0]]
	el=np.sqrt(xl*xl+yl*yl+zl*zl)
	lx=xl/el
	mx=yl/el
	nx=zl/el
	Di=np.sqrt(lx*lx+mx*mx)
	if Di==0 and nx>0:
		transf=np.array([[0,0,1,0,0,0,0,0,0,0,0,0],
		                 [0,1,0,0,0,0,0,0,0,0,0,0],
                         [-1,0,0,0,0,0,0,0,0,0,0,0],
                         [0,0,0,0,0,1,0,0,0,0,0,0],
                         [0,0,0,0,1,0,0,0,0,0,0,0],
                         [0,0,0,-1,0,0,0,0,0,0,0,0],
                         [0,0,0,0,0,0,0,0,1,0,0,0],
                         [0,0,0,0,0,0,0,1,0,0,0,0],
                         [0,0,0,0,0,0,-1,0,0,0,0,0],
                         [0,0,0,0,0,0,0,0,0,0,0,1],
                         [0,0,0,0,0,0,0,0,0,0,1,0],
                         [0,0,0,0,0,0,0,0,0,-1,0,0]])
	elif Di==0 and nx<0:
		transf=np.array([[0,0,-1,0,0,0,0,0,0,0,0,0],
                         [0,1,0,0,0,0,0,0,0,0,0,0],
                         [1,0,0,0,0,0,0,0,0,0,0,0],
                         [0,0,0,0,0,-1,0,0,0,0,0,0],
                         [0,0,0,0,1,0,0,0,0,0,0,0],
                         [0,0,0,1,0,0,0,0,0,0,0,0],
                         [0,0,0,0,0,0,0,0,-1,0,0,0],
                         [0,0,0,0,0,0,0,1,0,0,0,0],
                         [0,0,0,0,0,0,1,0,0,0,0,0],
                         [0,0,0,0,0,0,0,0,0,0,0,-1],
                         [0,0,0,0,0,0,0,0,0,0,1,0],
                         [0,0,0,0,0,0,0,0,0,1,0,0]])
	else:
		ly=-mx/Di
		my=lx/Di
		ny=0
		lz=-lx*nx/Di
		mz=-mx*nx/Di
		nz=Di
		transf=np.array([[lx,mx,nx,0,0,0,0,0,0,0,0,0],
                         [ly,my,ny,0,0,0,0,0,0,0,0,0],
                         [lz,mz,nz,0,0,0,0,0,0,0,0,0],
                         [0,0,0,lx,mx,nx,0,0,0,0,0,0],
                         [0,0,0,ly,my,ny,0,0,0,0,0,0],
                         [0,0,0,lz,mz,nz,0,0,0,0,0,0],
                         [0,0,0,0,0,0,lx,mx,nx,0,0,0],
                         [0,0,0,0,0,0,ly,my,ny,0,0,0],
                         [0,0,0,0,0,0,lz,mz,nz,0,0,0],
                         [0,0,0,0,0,0,0,0,0,lx,mx,nx],
                         [0,0,0,0,0,0,0,0,0,ly,my,ny],
                         [0,0,0,0,0,0,0,0,0,lz,mz,nz]])

	klocal=np.array([[A*E/el,0,0,0,0,0,-A*E/el,0,0,0,0,0],
		[0,12*E*Iz/el**3,0,0,0,6*E*Iz/el**2,0,-12*E*Iz/el**3,0,0,0,6*E*Iz/el**2],
		[0,0,12*E*Iy/el**3,0,-6*E*Iy/el**2,0,0,0,-12*E*Iy/el**3,0,-6*E*Iy/el**2,0],
		[0,0,0,Gs*Jt/el,0,0,0,0,0,-Gs*Jt/el,0,0],
		[0,0,-6*E*Iy/el**2,0,4*E*Iy/el,0,0,0,6*E*Iy/el**2,0,2*E*Iy/el,0],
		[0,6*E*Iz/el**2,0,0,0,4*E*Iz/el,0,-6*E*Iz/el**2,0,0,0,2*E*Iz/el],
		[-A*E/el,0,0,0,0,0,A*E/el,0,0,0,0,0],
		[0,-12*E*Iz/el**3,0,0,0,-6*E*Iz/el**2,0,12*E*Iz/el**3,0,0,0,-6*E*Iz/el**2],
		[0,0,-12*E*Iy/el**3,0,6*E*Iy/el**2,0,0,0,12*E*Iy/el**3,0,6*E*Iy/el**2,0],
		[0,0,0,-Gs*Jt/el,0,0,0,0,0,Gs*Jt/el,0,0],
		[0,0,-6*E*Iy/el**2,0,2*E*Iy/el,0,0,0,6*E*Iy/el**2,0,4*E*Iy/el,0],
		[0,6*E*Iz/el**2,0,0,0,2*E*Iz/el,0,-6*E*Iz/el**2,0,0,0,4*E*Iz/el]])
	tptf=np.transpose(transf)
	tpklprod=np.dot(tptf,klocal)
	kglobal=np.dot(tpklprod,transf)
	stiffness[np.ix_(elem_dof, elem_dof)] +=kglobal

suppress_dof=np.array([6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23])
active_dof=np.setdiff1d(np.arange(struct_dof), suppress_dof)
displacement_aux=np.linalg.solve(stiffness[np.ix_(active_dof,active_dof)], force[np.ix_(active_dof)])
displacement[np.ix_(active_dof)]=displacement_aux
print(f'displacements:\n {displacement}')

for e in range(elem_count):
	print(f'forces on element {e+1}')
	index=elem_con[e]
	elem_dof=np.array([index[0]*6, index[0]*6+1, index[0]*6+2,
                       index[0]*6+3, index[0]*6+4, index[0]*6+5,
                       index[1]*6, index[1]*6+1, index[1]*6+2,
                       index[1]*6+3, index[1]*6+4, index[1]*6+5])
	xl=xx[index[1]]-xx[index[0]]
	yl=yy[index[1]]-yy[index[0]]
	zl=zz[index[1]]-zz[index[0]]
	el=np.sqrt(xl*xl+yl*yl+zl*zl)
	lx=xl/el
	mx=yl/el
	nx=zl/el
	Di=np.sqrt(lx*lx+mx*mx)
	if Di==0 and nx>0:
		transf=np.array([[0,0,1,0,0,0,0,0,0,0,0,0],
		                 [0,1,0,0,0,0,0,0,0,0,0,0],
                         [-1,0,0,0,0,0,0,0,0,0,0,0],
                         [0,0,0,0,0,1,0,0,0,0,0,0],
                         [0,0,0,0,1,0,0,0,0,0,0,0],
                         [0,0,0,-1,0,0,0,0,0,0,0,0],
                         [0,0,0,0,0,0,0,0,1,0,0,0],
                         [0,0,0,0,0,0,0,1,0,0,0,0],
                         [0,0,0,0,0,0,-1,0,0,0,0,0],
                         [0,0,0,0,0,0,0,0,0,0,0,1],
                         [0,0,0,0,0,0,0,0,0,0,1,0],
                         [0,0,0,0,0,0,0,0,0,-1,0,0]])
	elif Di==0 and nx<0:
		transf=np.array([[0,0,-1,0,0,0,0,0,0,0,0,0],
                         [0,1,0,0,0,0,0,0,0,0,0,0],
                         [1,0,0,0,0,0,0,0,0,0,0,0],
                         [0,0,0,0,0,-1,0,0,0,0,0,0],
                         [0,0,0,0,1,0,0,0,0,0,0,0],
                         [0,0,0,1,0,0,0,0,0,0,0,0],
                         [0,0,0,0,0,0,0,0,-1,0,0,0],
                         [0,0,0,0,0,0,0,1,0,0,0,0],
                         [0,0,0,0,0,0,1,0,0,0,0,0],
                         [0,0,0,0,0,0,0,0,0,0,0,-1],
                         [0,0,0,0,0,0,0,0,0,0,1,0],
                         [0,0,0,0,0,0,0,0,0,1,0,0]])
	else:
		ly=-mx/Di
		my=lx/Di
		ny=0
		lz=-lx*nx/Di
		mz=-mx*nx/Di
		nz=Di
		transf=np.array([[lx,mx,nx,0,0,0,0,0,0,0,0,0],
                         [ly,my,ny,0,0,0,0,0,0,0,0,0],
                         [lz,mz,nz,0,0,0,0,0,0,0,0,0],
                         [0,0,0,lx,mx,nx,0,0,0,0,0,0],
                         [0,0,0,ly,my,ny,0,0,0,0,0,0],
                         [0,0,0,lz,mz,nz,0,0,0,0,0,0],
                         [0,0,0,0,0,0,lx,mx,nx,0,0,0],
                         [0,0,0,0,0,0,ly,my,ny,0,0,0],
                         [0,0,0,0,0,0,lz,mz,nz,0,0,0],
                         [0,0,0,0,0,0,0,0,0,lx,mx,nx],
                         [0,0,0,0,0,0,0,0,0,ly,my,ny],
                         [0,0,0,0,0,0,0,0,0,lz,mz,nz]])

	Td=np.dot(transf, displacement[np.ix_(elem_dof)])
	klocal=np.array([[A*E/el,0,0,0,0,0,-A*E/el,0,0,0,0,0],
		[0,12*E*Iz/el**3,0,0,0,6*E*Iz/el**2,0,-12*E*Iz/el**3,0,0,0,6*E*Iz/el**2],
		[0,0,12*E*Iy/el**3,0,-6*E*Iy/el**2,0,0,0,-12*E*Iy/el**3,0,-6*E*Iy/el**2,0],
		[0,0,0,Gs*Jt/el,0,0,0,0,0,-Gs*Jt/el,0,0],
		[0,0,-6*E*Iy/el**2,0,4*E*Iy/el,0,0,0,6*E*Iy/el**2,0,2*E*Iy/el,0],
		[0,6*E*Iz/el**2,0,0,0,4*E*Iz/el,0,-6*E*Iz/el**2,0,0,0,2*E*Iz/el],
		[-A*E/el,0,0,0,0,0,A*E/el,0,0,0,0,0],
		[0,-12*E*Iz/el**3,0,0,0,-6*E*Iz/el**2,0,12*E*Iz/el**3,0,0,0,-6*E*Iz/el**2],
		[0,0,-12*E*Iy/el**3,0,6*E*Iy/el**2,0,0,0,12*E*Iy/el**3,0,6*E*Iy/el**2,0],
		[0,0,0,-Gs*Jt/el,0,0,0,0,0,Gs*Jt/el,0,0],
		[0,0,-6*E*Iy/el**2,0,2*E*Iy/el,0,0,0,6*E*Iy/el**2,0,4*E*Iy/el,0],
		[0,6*E*Iz/el**2,0,0,0,2*E*Iz/el,0,-6*E*Iz/el**2,0,0,0,4*E*Iz/el]])
	elem_forces=np.dot(klocal,Td)
	print(f'{elem_forces}')

The output:

$ python frame3d.py
displacements:
 [[ 7.098e-05]
 [-1.400e-02]
 [-2.352e-03]
 [-3.996e-03]
 [ 1.780e-05]
 [-1.033e-04]
 [ 0.000e+00]
 [ 0.000e+00]
 [ 0.000e+00]
 [ 0.000e+00]
 [ 0.000e+00]
 [ 0.000e+00]
 [ 0.000e+00]
 [ 0.000e+00]
 [ 0.000e+00]
 [ 0.000e+00]
 [ 0.000e+00]
 [ 0.000e+00]
 [ 0.000e+00]
 [ 0.000e+00]
 [ 0.000e+00]
 [ 0.000e+00]
 [ 0.000e+00]
 [ 0.000e+00]]
forces on element 1
[[ -0.213]
 [  0.318]
 [  0.053]
 [ 19.98 ]
 [ -3.165]
 [ 18.991]
 [  0.213]
 [ -0.318]
 [ -0.053]
 [-19.98 ]
 [ -2.097]
 [ 12.79 ]]
forces on element 2
[[ 7.056e+00]
 [ 7.697e+00]
 [-2.949e-02]
 [ 5.167e-01]
 [ 9.403e-01]
 [ 2.650e+02]
 [-7.056e+00]
 [-7.697e+00]
 [ 2.949e-02]
 [-5.167e-01]
 [ 2.008e+00]
 [ 5.047e+02]]
forces on element 3
[[ 4.199e+01]
 [-1.835e-01]
 [-7.108e+00]
 [-8.900e-02]
 [ 2.355e+02]
 [-6.073e+00]
 [-4.199e+01]
 [ 1.835e-01]
 [ 7.108e+00]
 [ 8.900e-02]
 [ 4.753e+02]
 [-1.227e+01]]

This last script craves for refactoring, but hey, it works. At the time I wanted to find some visualization library, rewrite and expand the code. I thought I’d post it only after it was finished. Fast forward to 2022, I realized I won’t do it anytime soon as I don’t design structural things anymore. There it is, anyway. Hopefully it can help someone.

Cheers,

Edi

Direct Stiffness Method in Python: Trusses and Frames

Table of Contents

Trusses

2D Truss

Space Truss

Frames

Portal Frame

3D Frame