<jupyter_start><jupyter_text>Gauss integration of Finite Elements<jupyter_code>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()
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')<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>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<jupyter_output><empty_output><jupyter_text>Interpolation functions The enumeration that we are using for the elements is shown below<jupyter_code>Image(filename='../img_src/4node_element_enumeration.png', width=300)<jupyter_output><empty_output><jupyter_text>What leads to the following shape functions<jupyter_code>N = S(1)/4*Matrix([(1 + r)*(1 + s),
         (1 - r)*(1 + s),
         (1 - r)*(1 - s),
         (1 + r)*(1 - s)])
N<jupyter_output><empty_output><jupyter_text>Thus, the interpolation matrix renders<jupyter_code>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<jupyter_output><empty_output><jupyter_text>And the mass matrix integrand is$$M_\text{int}=H^TH$$<jupyter_code>M_int = H.T*H<jupyter_output><empty_output><jupyter_text>Derivatives interpolation matrix<jupyter_code>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<jupyter_output><empty_output><jupyter_text>Being the stiffness matrix integrand$$K_\text{int} = B^T C B$$<jupyter_code>K_int = B.T*C*B<jupyter_output><empty_output><jupyter_text>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*}$$<jupyter_code>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<jupyter_output><empty_output><jupyter_text>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*}$$<jupyter_code>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<jupyter_output><empty_output><jupyter_text>We can generate automatically code for `Fortran`, `C` or `Octave/Matlab`, although it will be useful just for non-distorted elements.<jupyter_code>K_local = MatrixSymbol('K_local', 8, 8)
code = codegen(("local_stiff", Eq(K_local, simplify(K))), "f95")
print code[0][1]
code = codegen(("local_stiff", Eq(K_local, simplify(K))), "C")
print code[0][1]
code = codegen(("local_stiff", Eq(K_local, simplify(K))), "Octave")
print code[0][1]<jupyter_output>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 +[...]<jupyter_text>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<jupyter_code>(C_factor*K).subs([(E, S(8)/3), (nu, S(1)/3)])
M.subs(rho, 1)<jupyter_output><empty_output><jupyter_text>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<jupyter_code>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))]]<jupyter_output><empty_output><jupyter_text>And the numerical integral is computed as<jupyter_code>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)
    
K_num = stiff_num(3)
K_num - K<jupyter_output><empty_output><jupyter_text>Best approach A best approach is to use Python built-in functions for computing the Gauss-Legendre nodes and weights<jupyter_code>from sympy.integrals.quadrature import gauss_legendre
x, w = gauss_legendre(5, 15)
print x
print w<jupyter_output>[-0.906179845938664, -0.538469310105683, 0, 0.538469310105683, 0.906179845938664]
[0.236926885056189, 0.478628670499366, 0.568888888888889, 0.478628670499366, 0.236926885056189]<jupyter_text>References [1] http://en.wikipedia.org/wiki/Gaussian_quadrature<jupyter_code>from IPython.core.display import HTML
def css_styling():
    styles = open('../styles/custom_barba.css', 'r').read()
    return HTML(styles)
css_styling()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17842702491196036, 'fraction_numerical': 0.06782313812442937, 'mean_word_length': 3.286193404136389, 'pattern_counts': {'":': 0, '<': 63, '<?xml version=': 0, '>': 64, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '308574', 'license': 'MIT', 'n_tokens_mistral': 3364, 'n_tokens_neox': 3121, 'n_words': 740, 'path': 'notebooks/gauss_integration.ipynb', 'repo_name': 'jomorlier/FEM-Notes', 'size': '7667'}