# Gauss integration of Finite Elements


```python
from __future__ import division
from sympy.utilities.codegen import codegen
from sympy import *
from sympy import init_printing
from IPython.display import Image
init_printing()
```


```python
r, s, t, x, y, z = symbols('r s t x y z') 
k, m, n = symbols('k m n', integer=True)
rho, nu, E = symbols('rho, nu, E')
```

## Predefinition

The constitutive model tensor in Voigt notation (plane strain) is
$$C = \frac{(1 - \nu) E}{(1 - 2\nu) (1 + \nu) }
\begin{pmatrix}
1 & \frac{\nu}{1-\nu} & 0\\
\frac{\nu}{1-\nu} & 1 & 0\\
0 & 0 & \frac{1 - 2\nu}{2(1 - \nu)}
\end{pmatrix}$$

But for simplicity we are going to use
$$\hat{C} = \frac{C (1 - 2\nu) (1 + \nu)}{E} = 
\begin{pmatrix}
1-\nu & \nu & 0\\
\nu & 1-\nu & 0\\
0 & 0 & \frac{1 - 2\nu}{2}
\end{pmatrix} \enspace ,$$
since we can always multiply by that factor afterwards to obtain the correct stiffness matrices.


```python
C = Matrix([[1 - nu, nu, 0],
           [nu, 1 - nu, 0],
           [0, 0, (1 - 2*nu)/2]])
C_factor = E/(1-2*nu)/(1 + nu)
C
```

## Interpolation functions

The enumeration that we are using for the elements is shown below


```python
Image(filename='../img_src/4node_element_enumeration.png', width=300)
```

What leads to the following shape functions


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

Thus, the interpolation matrix renders


```python
H = zeros(2,8)
for i in range(4):
    H[0, 2*i] = N[i]
    H[1, 2*i + 1] = N[i]
    
H.T  # Transpose of the interpolation matrix
```

And the mass matrix integrand is
$$M_\text{int}=H^TH$$


```python
M_int = H.T*H
```

## Derivatives interpolation matrix


```python
dHdr = zeros(2,4)
for i in range(4):
    dHdr[0,i] = diff(N[i],r)
    dHdr[1,i] = diff(N[i],s)

jaco = eye(2)  # Jacobian matrix, identity for now
dHdx = jaco*dHdr

B = zeros(3,8)
for i in range(4):
    B[0, 2*i] = dHdx[0, i]
    B[1, 2*i+1] = dHdx[1, i]
    B[2, 2*i] = dHdx[1, i]
    B[2, 2*i+1] = dHdx[0, i]
    
B
```

Being the stiffness matrix integrand
$$K_\text{int} = B^T C B$$


```python
K_int = B.T*C*B
```

## Analytic integration

The mass matrix is obtained integrating the product of the interpolator matrix with itself, i.e.
$$\begin{align*}
M &= \int\limits_{-1}^{1}\int\limits_{-1}^{1} M_\text{int} dr\, ds\\
 &= \int\limits_{-1}^{1}\int\limits_{-1}^{1} H^T H\, dr\, ds \enspace .
\end{align*}$$


```python
M = zeros(8,8)
for i in range(8):
    for j in range(8):
        M[i,j] = rho*integrate(M_int[i,j],(r,-1,1), (s,-1,1))
        
M
```

The stiffness matrix is obtained integrating the product of the interpolator-derivatives (displacement-to-strains) matrix with the constitutive tensor and itself, i.e.
$$\begin{align*}
K &= \int\limits_{-1}^{1}\int\limits_{-1}^{1} K_\text{int} dr\, ds\\
 &= \int\limits_{-1}^{1}\int\limits_{-1}^{1} B^T C\, B\, dr\, ds \enspace .
\end{align*}$$


```python
K = zeros(8,8)
for i in range(8):
    for j in range(8):
        K[i,j] = integrate(K_int[i,j], (r,-1,1), (s,-1,1))

K
```

We can generate automatically code for `Fortran`, `C` or `Octave/Matlab`, although it will be useful just for non-distorted elements.


```python
K_local = MatrixSymbol('K_local', 8, 8)
code = codegen(("local_stiff", Eq(K_local, simplify(K))), "f95")
print code[0][1]
```

    !******************************************************************************
    !*                      Code generated with sympy 0.7.6                       *
    !*                                                                            *
    !*              See http://www.sympy.org/ for more information.               *
    !*                                                                            *
    !*                       This file is part of 'project'                       *
    !******************************************************************************
    
    subroutine local_stiff(nu, K_local)
    implicit none
    REAL*8, intent(in) :: nu
    REAL*8, intent(out), dimension(1:8, 1:8) :: K_local
    
    K_local(1, 1) = -2.0d0/3.0d0*nu + 1.0d0/2.0d0
    K_local(2, 1) = 1.0d0/8.0d0
    K_local(3, 1) = (1.0d0/6.0d0)*nu - 1.0d0/4.0d0
    K_local(4, 1) = (1.0d0/2.0d0)*nu - 1.0d0/8.0d0
    K_local(5, 1) = (1.0d0/3.0d0)*nu - 1.0d0/4.0d0
    K_local(6, 1) = -1.0d0/8.0d0
    K_local(7, 1) = (1.0d0/6.0d0)*nu
    K_local(8, 1) = -1.0d0/2.0d0*nu + 1.0d0/8.0d0
    K_local(1, 2) = 1.0d0/8.0d0
    K_local(2, 2) = -2.0d0/3.0d0*nu + 1.0d0/2.0d0
    K_local(3, 2) = -1.0d0/2.0d0*nu + 1.0d0/8.0d0
    K_local(4, 2) = (1.0d0/6.0d0)*nu
    K_local(5, 2) = -1.0d0/8.0d0
    K_local(6, 2) = (1.0d0/3.0d0)*nu - 1.0d0/4.0d0
    K_local(7, 2) = (1.0d0/2.0d0)*nu - 1.0d0/8.0d0
    K_local(8, 2) = (1.0d0/6.0d0)*nu - 1.0d0/4.0d0
    K_local(1, 3) = (1.0d0/6.0d0)*nu - 1.0d0/4.0d0
    K_local(2, 3) = -1.0d0/2.0d0*nu + 1.0d0/8.0d0
    K_local(3, 3) = -2.0d0/3.0d0*nu + 1.0d0/2.0d0
    K_local(4, 3) = -1.0d0/8.0d0
    K_local(5, 3) = (1.0d0/6.0d0)*nu
    K_local(6, 3) = (1.0d0/2.0d0)*nu - 1.0d0/8.0d0
    K_local(7, 3) = (1.0d0/3.0d0)*nu - 1.0d0/4.0d0
    K_local(8, 3) = 1.0d0/8.0d0
    K_local(1, 4) = (1.0d0/2.0d0)*nu - 1.0d0/8.0d0
    K_local(2, 4) = (1.0d0/6.0d0)*nu
    K_local(3, 4) = -1.0d0/8.0d0
    K_local(4, 4) = -2.0d0/3.0d0*nu + 1.0d0/2.0d0
    K_local(5, 4) = -1.0d0/2.0d0*nu + 1.0d0/8.0d0
    K_local(6, 4) = (1.0d0/6.0d0)*nu - 1.0d0/4.0d0
    K_local(7, 4) = 1.0d0/8.0d0
    K_local(8, 4) = (1.0d0/3.0d0)*nu - 1.0d0/4.0d0
    K_local(1, 5) = (1.0d0/3.0d0)*nu - 1.0d0/4.0d0
    K_local(2, 5) = -1.0d0/8.0d0
    K_local(3, 5) = (1.0d0/6.0d0)*nu
    K_local(4, 5) = -1.0d0/2.0d0*nu + 1.0d0/8.0d0
    K_local(5, 5) = -2.0d0/3.0d0*nu + 1.0d0/2.0d0
    K_local(6, 5) = 1.0d0/8.0d0
    K_local(7, 5) = (1.0d0/6.0d0)*nu - 1.0d0/4.0d0
    K_local(8, 5) = (1.0d0/2.0d0)*nu - 1.0d0/8.0d0
    K_local(1, 6) = -1.0d0/8.0d0
    K_local(2, 6) = (1.0d0/3.0d0)*nu - 1.0d0/4.0d0
    K_local(3, 6) = (1.0d0/2.0d0)*nu - 1.0d0/8.0d0
    K_local(4, 6) = (1.0d0/6.0d0)*nu - 1.0d0/4.0d0
    K_local(5, 6) = 1.0d0/8.0d0
    K_local(6, 6) = -2.0d0/3.0d0*nu + 1.0d0/2.0d0
    K_local(7, 6) = -1.0d0/2.0d0*nu + 1.0d0/8.0d0
    K_local(8, 6) = (1.0d0/6.0d0)*nu
    K_local(1, 7) = (1.0d0/6.0d0)*nu
    K_local(2, 7) = (1.0d0/2.0d0)*nu - 1.0d0/8.0d0
    K_local(3, 7) = (1.0d0/3.0d0)*nu - 1.0d0/4.0d0
    K_local(4, 7) = 1.0d0/8.0d0
    K_local(5, 7) = (1.0d0/6.0d0)*nu - 1.0d0/4.0d0
    K_local(6, 7) = -1.0d0/2.0d0*nu + 1.0d0/8.0d0
    K_local(7, 7) = -2.0d0/3.0d0*nu + 1.0d0/2.0d0
    K_local(8, 7) = -1.0d0/8.0d0
    K_local(1, 8) = -1.0d0/2.0d0*nu + 1.0d0/8.0d0
    K_local(2, 8) = (1.0d0/6.0d0)*nu - 1.0d0/4.0d0
    K_local(3, 8) = 1.0d0/8.0d0
    K_local(4, 8) = (1.0d0/3.0d0)*nu - 1.0d0/4.0d0
    K_local(5, 8) = (1.0d0/2.0d0)*nu - 1.0d0/8.0d0
    K_local(6, 8) = (1.0d0/6.0d0)*nu
    K_local(7, 8) = -1.0d0/8.0d0
    K_local(8, 8) = -2.0d0/3.0d0*nu + 1.0d0/2.0d0
    
    end subroutine
    



```python
code = codegen(("local_stiff", Eq(K_local, simplify(K))), "C")
print code[0][1]
```

    /******************************************************************************
     *                      Code generated with sympy 0.7.6                       *
     *                                                                            *
     *              See http://www.sympy.org/ for more information.               *
     *                                                                            *
     *                       This file is part of 'project'                       *
     ******************************************************************************/
    #include "local_stiff.h"
    #include <math.h>
    
    void local_stiff(double nu, double *K_local) {
    
       K_local[0] = -2.0L/3.0L*nu + 1.0L/2.0L;
       K_local[1] = 1.0L/8.0L;
       K_local[2] = (1.0L/6.0L)*nu - 1.0L/4.0L;
       K_local[3] = (1.0L/2.0L)*nu - 1.0L/8.0L;
       K_local[4] = (1.0L/3.0L)*nu - 1.0L/4.0L;
       K_local[5] = -1.0L/8.0L;
       K_local[6] = (1.0L/6.0L)*nu;
       K_local[7] = -1.0L/2.0L*nu + 1.0L/8.0L;
       K_local[8] = 1.0L/8.0L;
       K_local[9] = -2.0L/3.0L*nu + 1.0L/2.0L;
       K_local[10] = -1.0L/2.0L*nu + 1.0L/8.0L;
       K_local[11] = (1.0L/6.0L)*nu;
       K_local[12] = -1.0L/8.0L;
       K_local[13] = (1.0L/3.0L)*nu - 1.0L/4.0L;
       K_local[14] = (1.0L/2.0L)*nu - 1.0L/8.0L;
       K_local[15] = (1.0L/6.0L)*nu - 1.0L/4.0L;
       K_local[16] = (1.0L/6.0L)*nu - 1.0L/4.0L;
       K_local[17] = -1.0L/2.0L*nu + 1.0L/8.0L;
       K_local[18] = -2.0L/3.0L*nu + 1.0L/2.0L;
       K_local[19] = -1.0L/8.0L;
       K_local[20] = (1.0L/6.0L)*nu;
       K_local[21] = (1.0L/2.0L)*nu - 1.0L/8.0L;
       K_local[22] = (1.0L/3.0L)*nu - 1.0L/4.0L;
       K_local[23] = 1.0L/8.0L;
       K_local[24] = (1.0L/2.0L)*nu - 1.0L/8.0L;
       K_local[25] = (1.0L/6.0L)*nu;
       K_local[26] = -1.0L/8.0L;
       K_local[27] = -2.0L/3.0L*nu + 1.0L/2.0L;
       K_local[28] = -1.0L/2.0L*nu + 1.0L/8.0L;
       K_local[29] = (1.0L/6.0L)*nu - 1.0L/4.0L;
       K_local[30] = 1.0L/8.0L;
       K_local[31] = (1.0L/3.0L)*nu - 1.0L/4.0L;
       K_local[32] = (1.0L/3.0L)*nu - 1.0L/4.0L;
       K_local[33] = -1.0L/8.0L;
       K_local[34] = (1.0L/6.0L)*nu;
       K_local[35] = -1.0L/2.0L*nu + 1.0L/8.0L;
       K_local[36] = -2.0L/3.0L*nu + 1.0L/2.0L;
       K_local[37] = 1.0L/8.0L;
       K_local[38] = (1.0L/6.0L)*nu - 1.0L/4.0L;
       K_local[39] = (1.0L/2.0L)*nu - 1.0L/8.0L;
       K_local[40] = -1.0L/8.0L;
       K_local[41] = (1.0L/3.0L)*nu - 1.0L/4.0L;
       K_local[42] = (1.0L/2.0L)*nu - 1.0L/8.0L;
       K_local[43] = (1.0L/6.0L)*nu - 1.0L/4.0L;
       K_local[44] = 1.0L/8.0L;
       K_local[45] = -2.0L/3.0L*nu + 1.0L/2.0L;
       K_local[46] = -1.0L/2.0L*nu + 1.0L/8.0L;
       K_local[47] = (1.0L/6.0L)*nu;
       K_local[48] = (1.0L/6.0L)*nu;
       K_local[49] = (1.0L/2.0L)*nu - 1.0L/8.0L;
       K_local[50] = (1.0L/3.0L)*nu - 1.0L/4.0L;
       K_local[51] = 1.0L/8.0L;
       K_local[52] = (1.0L/6.0L)*nu - 1.0L/4.0L;
       K_local[53] = -1.0L/2.0L*nu + 1.0L/8.0L;
       K_local[54] = -2.0L/3.0L*nu + 1.0L/2.0L;
       K_local[55] = -1.0L/8.0L;
       K_local[56] = -1.0L/2.0L*nu + 1.0L/8.0L;
       K_local[57] = (1.0L/6.0L)*nu - 1.0L/4.0L;
       K_local[58] = 1.0L/8.0L;
       K_local[59] = (1.0L/3.0L)*nu - 1.0L/4.0L;
       K_local[60] = (1.0L/2.0L)*nu - 1.0L/8.0L;
       K_local[61] = (1.0L/6.0L)*nu;
       K_local[62] = -1.0L/8.0L;
       K_local[63] = -2.0L/3.0L*nu + 1.0L/2.0L;
    
    }
    



```python
code = codegen(("local_stiff", Eq(K_local, simplify(K))), "Octave")
print code[0][1]
```

    function K_local = local_stiff(nu)
      %LOCAL_STIFF  Autogenerated by sympy
      %   Code generated with sympy 0.7.6
      %
      %   See http://www.sympy.org/ for more information.
      %
      %   This file is part of 'project'
    
      K_local = [-2*nu/3 + 1/2           1/8    nu/6 - 1/4    nu/2 - 1/8    nu/3 - 1/4          -1/8          nu/6   -nu/2 + 1/8;
      1/8 -2*nu/3 + 1/2   -nu/2 + 1/8          nu/6          -1/8    nu/3 - 1/4    nu/2 - 1/8    nu/6 - 1/4;
      nu/6 - 1/4   -nu/2 + 1/8 -2*nu/3 + 1/2          -1/8          nu/6    nu/2 - 1/8    nu/3 - 1/4           1/8;
      nu/2 - 1/8          nu/6          -1/8 -2*nu/3 + 1/2   -nu/2 + 1/8    nu/6 - 1/4           1/8    nu/3 - 1/4;
      nu/3 - 1/4          -1/8          nu/6   -nu/2 + 1/8 -2*nu/3 + 1/2           1/8    nu/6 - 1/4    nu/2 - 1/8;
      -1/8    nu/3 - 1/4    nu/2 - 1/8    nu/6 - 1/4           1/8 -2*nu/3 + 1/2   -nu/2 + 1/8          nu/6;
      nu/6    nu/2 - 1/8    nu/3 - 1/4           1/8    nu/6 - 1/4   -nu/2 + 1/8 -2*nu/3 + 1/2          -1/8;
      -nu/2 + 1/8    nu/6 - 1/4           1/8    nu/3 - 1/4    nu/2 - 1/8          nu/6          -1/8 -2*nu/3 + 1/2];
    
    end
    


We can check some numerical vales for $E=8/3$ Pa, $\nu=1/3$ and $\rho=1$ kg\m$^3$, where we can multiply by the factor that we took away from the stiffness tensor


```python
(C_factor*K).subs([(E, S(8)/3), (nu, S(1)/3)])
```


```python
M.subs(rho, 1)
```

## Gauss integration

As stated before, the analytic expressions for the mass and stiffness matrices is useful for non-distorted elements. In the general case, a mapping between distorted elements and these _canonical_ elements is used to simplify the integration domain. When this transformation is done, the functions to be integrated are more convoluted and we should use numerical integration like _Gauss-Legendre quadrature_.

The Gauss-Legendre quadrature approximates the integral:
$$     \int_{-1}^1 f(x)\,dx \approx \sum_{i=1}^n w_i f(x_i)$$

The nodes $x_i$ of an order $n$ quadrature rule are the roots of $P_n$
and the weights $w_i$ are given by:
$$w_i = \frac{2}{\left(1-x_i^2\right) \left(P'_n(x_i)\right)^2}$$

For the first four orders, the weights and nodes are


```python
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(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)), sqrt(S(3)/7 + S(2)/7*sqrt(S(6)/5))]]
```

And the numerical integral is computed as


```python
def stiff_num(n):
    """Compute the stiffness matrix using Gauss quadrature
    
    Parameters
    ----------
    n : int
        Order of the polynomial.
        
    """
    if n>4:
        raise Exception("Number of points not valid")
    K_num = zeros(8,8)
    for x_i, w_i in zip(pts[n-1], wts[n-1]):
        for y_j, w_j in zip(pts[n-1], wts[n-1]):
            K_num = K_num + w_i*w_j*K_int.subs([(r,x_i), (s,y_j)])

    return simplify(K_num)
    
```


```python
K_num = stiff_num(3)
K_num - K
```

### Best approach

A best approach is to use Python built-in functions for computing the Gauss-Legendre nodes and weights


```python
from sympy.integrals.quadrature import gauss_legendre
```


```python
x, w = gauss_legendre(5, 15)
print x
print w
```

    [-0.906179845938664, -0.538469310105683, 0, 0.538469310105683, 0.906179845938664]
    [0.236926885056189, 0.478628670499366, 0.568888888888889, 0.478628670499366, 0.236926885056189]


## References

[1] http://en.wikipedia.org/wiki/Gaussian_quadrature


```python
from IPython.core.display import HTML
def css_styling():
    styles = open('../styles/custom_barba.css', 'r').read()
    return HTML(styles)
css_styling()
```





<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>
<style>
/* 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 )
}
</style>






d50044ddb978e876bfa666d6ba145ebb32f85c8d
ipynb
Jupyter Notebook
notebooks/gauss_integration.ipynb
jomorlier/FEM-Notes
3b81053aee79dc59965c3622bc0d0eb6cfc7e8ae
MIT
2020-04-15T01:53:14.000Z
2020-04-15T01:53:14.000Z
notebooks/gauss_integration.ipynb
jomorlier/FEM-Notes
3b81053aee79dc59965c3622bc0d0eb6cfc7e8ae
MIT
notebooks/gauss_integration.ipynb
jomorlier/FEM-Notes
3b81053aee79dc59965c3622bc0d0eb6cfc7e8ae
MIT
2020-05-25T17:19:53.000Z
2020-05-25T17:19:53.000Z
Qwen/Qwen-72B
1. YES
2. YES