Foundation of FEM

Spring_FEM_in_Python
“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]


Advertisements

3 thoughts on “Foundation of FEM

  1. Pingback: Class version of Spring FEM in Python | SukhbinderSingh.com

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s