# Foundation of FEM

“Spring is the foundation of FEM. Understand this well and you are through in this course.” said the IISC professor as we settled on the hard benches on the first day of our 6 month FEM course at IISC, bangalore.

The course proved very useful. Learnt a lot. Fortunately , being in NAL had access to matlab and so converted the learning into matlab scripts. Today those scripts must be lying in some system in NAL.

Two months back, working on a problem in office, revisited all those notes and had this impulse to revisit the spring FEM problem in python…

So here’s the result. Two problem defined and solved in python.

```
"""
Created on Tue April 06 18:46:01 2014

@author: Sukhbinder Singh

Finite Element Analysis of spring system

problem 1

#---^^^^^----#
#              #
#-----^^^^^----# ===> P
#              #
#---^^^^^----#

1              2            3,4

Spring Problem

P=100
ki=10
u1=u3=u4 =0

problem 2

#                             #
#------^^^^^^^^^^^^^^^^-------#
#                             #
#         |---^^^^--|         #
#---^^^^--|==>P     |--^^^^---#
#         |---^^^^--|         #
|                   #
|------^^^^^^-------#

1         2         3         4

k1=10,k2=15,k3=20,k4=25,k5=30,k6=35
P=100
u1=u4=0

Solution for problem 1
Displacement
0 0.0
1 3.33333333333
2 0.0
3 0.0

Reactions
0 -33.3333333333
2 -33.3333333333
3 -33.3333333333

Solution for problem 2
Displacement
0 0.0
1 1.54589371981
2 0.869565217391
3 0.0

Reactions
0 -23.1884057971
3 -76.8115942029

#usage of np.ix_
import numpy as np
A = np.array([ 1, 2, 3, 4, 4,
2, 3, 4, 5, 3,
4, 5, 6, 7, 2,
5, 6, 7, 8, 9,
6, 7, 8, 9, 0 ]).reshape(5, 5)

B = np.array([60, 70, 80, 90]).reshape(2, 2)
A[np.ix_([2, 4], [2, 4])] = B

# usage of np.setdiff1d(ar1,ar2)
# Return the sorted, unique values in `ar1` that are not in `ar2`.
a = np.array([1, 2, 3, 2, 4, 1])
b = np.array([3, 4, 5, 6])
np.setdiff1d(a, b)

"""
import numpy as np

p1 = True # True for problem 1

if p1: # problem 1
elemNodes= np.array([[0,1],[1,2],[1,3]])
numElem = elemNodes.shape[0]
numNodes=4
presDof=np.array([0,2,3])
k=np.array([10,10,10]) # Spring Stiffness
else: # problem 2
elemNodes= np.array([[0,3],[0,1],[1,2],[1,2],[1,3],[2,3]])
numElem = elemNodes.shape[0]
numNodes=4
presDof=np.array([0,3])
k=np.array([10,15,20,25,30,35]) # Spring Stiffness

disp=np.zeros((numNodes,1))
force=np.zeros((numNodes,1))
stiffness=np.zeros((numNodes,numNodes))

force[1]=100.0

for e in xrange(numElem):
stiffness[np.ix_(elemNodes[e,:],elemNodes[e,:])]+=np.array([[k[e],-k[e]],[-k[e],k[e]]])

actDof=np.setdiff1d(np.arange(numNodes),presDof)

disp1=np.linalg.solve(stiffness[np.ix_(actDof,actDof)],force[np.ix_(actDof)]);
disp[np.ix_(actDof)]=disp1

print "\nDisplacement"
for i,d in enumerate(disp):
print i,d[0]

# reactions
react = np.dot(stiffness,disp)
reaction = react[presDof]

print "\nReactions"
for (p,r) in zip(presDof,reaction):
print p,r[0]

```