# Analytic integration of 4-nodes Finite Element



In [1]:

from __future__ import division
from sympy import *




In [2]:

init_session()




IPython console for SymPy 1.2 (Python 3.6.8-64-bit) (ground types: gmpy)

These commands were executed:
>>> from __future__ import division
>>> from sympy import *
>>> x, y, z, t = symbols('x y z t')
>>> k, m, n = symbols('k m n', integer=True)
>>> f, g, h = symbols('f g h', cls=Function)
>>> init_printing()

Documentation can be found at http://docs.sympy.org/1.2/




In [3]:

E, rho = symbols('E rho', positive=True)
nu, r, s = symbols('nu r s')



Constitutive equations in Voigt's notation



In [4]:

factor = E/((1+nu)*(1-2*nu))
factor = 1
C = factor*Matrix([
[1 - nu,nu,0],
[nu,1,0],
[0,0,(1-2*nu)/2]])
display(C)




$$\left[\begin{matrix}- \nu + 1 & \nu & 0\\\nu & 1 & 0\\0 & 0 & - \nu + \frac{1}{2}\end{matrix}\right]$$



Interpolation functions

2-----1
|     |
|     |
3-----4


In [5]:

N = S(1)/4*Matrix([(1 + r) * (1 + s),
(1 - r) * (1 + s),
(1 - r) * (1 - s),
(1 + r) * (1 - s)])



Interpolation matrix



In [6]:

H = zeros(2, 8)
for i in range(4):
H[0, 2*i] = N[i]
H[1, 2*i + 1] = N[i]

display(H)




$$\left[\begin{matrix}\frac{\left(r + 1\right) \left(s + 1\right)}{4} & 0 & \frac{\left(- r + 1\right) \left(s + 1\right)}{4} & 0 & \frac{\left(- r + 1\right) \left(- s + 1\right)}{4} & 0 & \frac{\left(r + 1\right) \left(- s + 1\right)}{4} & 0\\0 & \frac{\left(r + 1\right) \left(s + 1\right)}{4} & 0 & \frac{\left(- r + 1\right) \left(s + 1\right)}{4} & 0 & \frac{\left(- r + 1\right) \left(- s + 1\right)}{4} & 0 & \frac{\left(r + 1\right) \left(- s + 1\right)}{4}\end{matrix}\right]$$



Mass matrix



In [7]:

M_inte = H.T * H
M = integrate(integrate(rho*M_inte, [r, -1, 1]), [s, -1, 1])
display(M)




$$\left[\begin{matrix}\frac{4 \rho}{9} & 0 & \frac{2 \rho}{9} & 0 & \frac{\rho}{9} & 0 & \frac{2 \rho}{9} & 0\\0 & \frac{4 \rho}{9} & 0 & \frac{2 \rho}{9} & 0 & \frac{\rho}{9} & 0 & \frac{2 \rho}{9}\\\frac{2 \rho}{9} & 0 & \frac{4 \rho}{9} & 0 & \frac{2 \rho}{9} & 0 & \frac{\rho}{9} & 0\\0 & \frac{2 \rho}{9} & 0 & \frac{4 \rho}{9} & 0 & \frac{2 \rho}{9} & 0 & \frac{\rho}{9}\\\frac{\rho}{9} & 0 & \frac{2 \rho}{9} & 0 & \frac{4 \rho}{9} & 0 & \frac{2 \rho}{9} & 0\\0 & \frac{\rho}{9} & 0 & \frac{2 \rho}{9} & 0 & \frac{4 \rho}{9} & 0 & \frac{2 \rho}{9}\\\frac{2 \rho}{9} & 0 & \frac{\rho}{9} & 0 & \frac{2 \rho}{9} & 0 & \frac{4 \rho}{9} & 0\\0 & \frac{2 \rho}{9} & 0 & \frac{\rho}{9} & 0 & \frac{2 \rho}{9} & 0 & \frac{4 \rho}{9}\end{matrix}\right]$$



Derivatives interpolation matrix



In [8]:

B = zeros(3, 8)
for i in range(4):
B[0, 2*i] = diff(N[i], r)
B[1, 2*i + 1] = diff(N[i], s)
B[2, 2*i] = diff(N[i], s)
B[2, 2*i + 1] = diff(N[i], r)

display(B)




$$\left[\begin{matrix}\frac{s}{4} + \frac{1}{4} & 0 & - \frac{s}{4} - \frac{1}{4} & 0 & \frac{s}{4} - \frac{1}{4} & 0 & - \frac{s}{4} + \frac{1}{4} & 0\\0 & \frac{r}{4} + \frac{1}{4} & 0 & - \frac{r}{4} + \frac{1}{4} & 0 & \frac{r}{4} - \frac{1}{4} & 0 & - \frac{r}{4} - \frac{1}{4}\\\frac{r}{4} + \frac{1}{4} & \frac{s}{4} + \frac{1}{4} & - \frac{r}{4} + \frac{1}{4} & - \frac{s}{4} - \frac{1}{4} & \frac{r}{4} - \frac{1}{4} & \frac{s}{4} - \frac{1}{4} & - \frac{r}{4} - \frac{1}{4} & - \frac{s}{4} + \frac{1}{4}\end{matrix}\right]$$



Stiffness matrix



In [9]:

K_inte = B.T * C * B
K = simplify(integrate(integrate(K_inte, [r, -1, 1]), [s, -1, 1]))




In [10]:

display(K)




$$\left[\begin{matrix}- \frac{2 \nu}{3} + \frac{1}{2} & \frac{1}{8} & \frac{\nu}{6} - \frac{1}{4} & \frac{\nu}{2} - \frac{1}{8} & \frac{\nu}{3} - \frac{1}{4} & - \frac{1}{8} & \frac{\nu}{6} & - \frac{\nu}{2} + \frac{1}{8}\\\frac{1}{8} & - \frac{\nu}{3} + \frac{1}{2} & - \frac{\nu}{2} + \frac{1}{8} & \frac{\nu}{3} & - \frac{1}{8} & \frac{\nu}{6} - \frac{1}{4} & \frac{\nu}{2} - \frac{1}{8} & - \frac{\nu}{6} - \frac{1}{4}\\\frac{\nu}{6} - \frac{1}{4} & - \frac{\nu}{2} + \frac{1}{8} & - \frac{2 \nu}{3} + \frac{1}{2} & - \frac{1}{8} & \frac{\nu}{6} & \frac{\nu}{2} - \frac{1}{8} & \frac{\nu}{3} - \frac{1}{4} & \frac{1}{8}\\\frac{\nu}{2} - \frac{1}{8} & \frac{\nu}{3} & - \frac{1}{8} & - \frac{\nu}{3} + \frac{1}{2} & - \frac{\nu}{2} + \frac{1}{8} & - \frac{\nu}{6} - \frac{1}{4} & \frac{1}{8} & \frac{\nu}{6} - \frac{1}{4}\\\frac{\nu}{3} - \frac{1}{4} & - \frac{1}{8} & \frac{\nu}{6} & - \frac{\nu}{2} + \frac{1}{8} & - \frac{2 \nu}{3} + \frac{1}{2} & \frac{1}{8} & \frac{\nu}{6} - \frac{1}{4} & \frac{\nu}{2} - \frac{1}{8}\\- \frac{1}{8} & \frac{\nu}{6} - \frac{1}{4} & \frac{\nu}{2} - \frac{1}{8} & - \frac{\nu}{6} - \frac{1}{4} & \frac{1}{8} & - \frac{\nu}{3} + \frac{1}{2} & - \frac{\nu}{2} + \frac{1}{8} & \frac{\nu}{3}\\\frac{\nu}{6} & \frac{\nu}{2} - \frac{1}{8} & \frac{\nu}{3} - \frac{1}{4} & \frac{1}{8} & \frac{\nu}{6} - \frac{1}{4} & - \frac{\nu}{2} + \frac{1}{8} & - \frac{2 \nu}{3} + \frac{1}{2} & - \frac{1}{8}\\- \frac{\nu}{2} + \frac{1}{8} & - \frac{\nu}{6} - \frac{1}{4} & \frac{1}{8} & \frac{\nu}{6} - \frac{1}{4} & \frac{\nu}{2} - \frac{1}{8} & \frac{\nu}{3} & - \frac{1}{8} & - \frac{\nu}{3} + \frac{1}{2}\end{matrix}\right]$$



## Gauss integration



In [11]:

wts = [2, [1,1], [S(5)/9, S(8)/9, S(5)/9],
[(18+sqrt(30))/36,(18+sqrt(30))/36, (18-sqrt(30))/36, (18-sqrt(30))/36]];
pts = [0, [-sqrt(S(1)/3), sqrt(S(1)/3)], [-sqrt(S(3)/5), 0, sqrt(S(3)/5)],
[-sqrt(S(3)/7 - S(2)/7*sqrt((6)/5)), sqrt(S(3)/7 - S(2)/7*sqrt(S(6)/5)),
-sqrt(S(3)/7 + S(2)/7*sqrt(S(6)/5)), sqrt(S(3)/7 + S(2)/7*sqrt(S(6)/5))]]




In [12]:

K_1pts = simplify(K_inte.subs({r:0, s:0})*wts[0])
display(K_1pts)




$$\left[\begin{matrix}- \frac{\nu}{4} + \frac{3}{16} & \frac{1}{16} & - \frac{1}{16} & \frac{\nu}{4} - \frac{1}{16} & \frac{\nu}{4} - \frac{3}{16} & - \frac{1}{16} & \frac{1}{16} & - \frac{\nu}{4} + \frac{1}{16}\\\frac{1}{16} & - \frac{\nu}{8} + \frac{3}{16} & - \frac{\nu}{4} + \frac{1}{16} & \frac{\nu}{8} + \frac{1}{16} & - \frac{1}{16} & \frac{\nu}{8} - \frac{3}{16} & \frac{\nu}{4} - \frac{1}{16} & - \frac{\nu}{8} - \frac{1}{16}\\- \frac{1}{16} & - \frac{\nu}{4} + \frac{1}{16} & - \frac{\nu}{4} + \frac{3}{16} & - \frac{1}{16} & \frac{1}{16} & \frac{\nu}{4} - \frac{1}{16} & \frac{\nu}{4} - \frac{3}{16} & \frac{1}{16}\\\frac{\nu}{4} - \frac{1}{16} & \frac{\nu}{8} + \frac{1}{16} & - \frac{1}{16} & - \frac{\nu}{8} + \frac{3}{16} & - \frac{\nu}{4} + \frac{1}{16} & - \frac{\nu}{8} - \frac{1}{16} & \frac{1}{16} & \frac{\nu}{8} - \frac{3}{16}\\\frac{\nu}{4} - \frac{3}{16} & - \frac{1}{16} & \frac{1}{16} & - \frac{\nu}{4} + \frac{1}{16} & - \frac{\nu}{4} + \frac{3}{16} & \frac{1}{16} & - \frac{1}{16} & \frac{\nu}{4} - \frac{1}{16}\\- \frac{1}{16} & \frac{\nu}{8} - \frac{3}{16} & \frac{\nu}{4} - \frac{1}{16} & - \frac{\nu}{8} - \frac{1}{16} & \frac{1}{16} & - \frac{\nu}{8} + \frac{3}{16} & - \frac{\nu}{4} + \frac{1}{16} & \frac{\nu}{8} + \frac{1}{16}\\\frac{1}{16} & \frac{\nu}{4} - \frac{1}{16} & \frac{\nu}{4} - \frac{3}{16} & \frac{1}{16} & - \frac{1}{16} & - \frac{\nu}{4} + \frac{1}{16} & - \frac{\nu}{4} + \frac{3}{16} & - \frac{1}{16}\\- \frac{\nu}{4} + \frac{1}{16} & - \frac{\nu}{8} - \frac{1}{16} & \frac{1}{16} & \frac{\nu}{8} - \frac{3}{16} & \frac{\nu}{4} - \frac{1}{16} & \frac{\nu}{8} + \frac{1}{16} & - \frac{1}{16} & - \frac{\nu}{8} + \frac{3}{16}\end{matrix}\right]$$




In [13]:

K_2pts = zeros(8, 8)
for row in range(8):
for col in range(8):
K_2pts[row, col] = sum(wts[1][kx]*wts[1][ky]*
K_inte[row, col].subs({r: pts[1][kx], s:pts[1][ky]})
for kx in range(2) for ky in range(2))




In [14]:

simplify(K_2pts - K)




Out[14]:

$$\left[\begin{matrix}0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\end{matrix}\right]$$




In [15]:

K_3pts = zeros(8, 8)
for row in range(8):
for col in range(8):
K_3pts[row, col] = sum(wts[2][kx]*wts[2][ky]*
K_inte[row, col].subs({r: pts[2][kx], s:pts[2][ky]})
for kx in range(3) for ky in range(3))




In [16]:

simplify(K_3pts - K)




Out[16]:

$$\left[\begin{matrix}0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\end{matrix}\right]$$




In [ ]:




In [ ]:



# References

1. Bathe, Klaus-Jürgen. Finite element procedures. Klaus-Jurgen Bathe, 2006.


In [17]:

from IPython.core.display import HTML
def css_styling():
return HTML(styles)
css_styling()




Out[17]:

/* Based on Lorena Barba template available at: https://github.com/barbagroup/AeroPython/blob/master/styles/custom.css*/
@font-face {
font-family: "Computer Modern";
src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');
}
div.cell{
width:800px;
margin-left:16% !important;
margin-right:auto;
}
h1 {
font-family: 'Alegreya Sans', sans-serif;
}
h2 {
font-family: 'Fenix', serif;
}
h3{
font-family: 'Fenix', serif;
margin-top:12px;
margin-bottom: 3px;
}
h4{
font-family: 'Fenix', serif;
}
h5 {
font-family: 'Alegreya Sans', sans-serif;
}
div.text_cell_render{
font-family: 'Alegreya Sans',Computer Modern, "Helvetica Neue", Arial, Helvetica, Geneva, sans-serif;
line-height: 135%;
font-size: 120%;
width:600px;
margin-left:auto;
margin-right:auto;
}
.CodeMirror{
font-family: "Source Code Pro";
font-size: 90%;
}
/* .prompt{
display: None;
}*/
.text_cell_render h1 {
font-weight: 200;
font-size: 50pt;
line-height: 100%;
color:#CD2305;
margin-bottom: 0.5em;
margin-top: 0.5em;
display: block;
}
.text_cell_render h5 {
font-weight: 300;
font-size: 16pt;
color: #CD2305;
font-style: italic;
margin-bottom: .5em;
margin-top: 0.5em;
display: block;
}
.warning{
color: rgb( 240, 20, 20 )
}

MathJax.Hub.Config({
TeX: {
extensions: ["AMSmath.js"]
},
tex2jax: {
inlineMath: [ ['$','$'], ["\$","\$"] ],
displayMath: [ ['$$','$$'], ["\$","\$"] ]
},
displayAlign: 'center', // Change this to 'center' to center equations.
"HTML-CSS": {
styles: {'.MathJax_Display': {"margin": 4}}
}
});




In [ ]: