Gauss integration of Finite Elements


In [1]:
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()

In [2]:
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.


In [3]:
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


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

Interpolation functions

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


In [4]:
Image(filename='../img_src/4node_element_enumeration.png', width=300)


Out[4]:

What leads to the following shape functions


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


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

Thus, the interpolation matrix renders


In [6]:
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


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

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


In [7]:
M_int = H.T*H

Derivatives interpolation matrix


In [8]:
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


Out[8]:
$$\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]$$

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


In [9]:
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*}$$


In [10]:
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


Out[10]:
$$\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]$$

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*}$$


In [11]:
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


Out[11]:
$$\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{2 \nu}{3} + \frac{1}{2} & - \frac{\nu}{2} + \frac{1}{8} & \frac{\nu}{6} & - \frac{1}{8} & \frac{\nu}{3} - \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}{6} & - \frac{1}{8} & - \frac{2 \nu}{3} + \frac{1}{2} & - \frac{\nu}{2} + \frac{1}{8} & \frac{\nu}{6} - \frac{1}{4} & \frac{1}{8} & \frac{\nu}{3} - \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}{3} - \frac{1}{4} & \frac{\nu}{2} - \frac{1}{8} & \frac{\nu}{6} - \frac{1}{4} & \frac{1}{8} & - \frac{2 \nu}{3} + \frac{1}{2} & - \frac{\nu}{2} + \frac{1}{8} & \frac{\nu}{6}\\\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}{3} - \frac{1}{4} & \frac{\nu}{2} - \frac{1}{8} & \frac{\nu}{6} & - \frac{1}{8} & - \frac{2 \nu}{3} + \frac{1}{2}\end{matrix}\right]$$

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


In [12]:
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


In [13]:
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;

}


In [14]:
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


In [15]:
(C_factor*K).subs([(E, S(8)/3), (nu, S(1)/3)])


Out[15]:
$$\left[\begin{matrix}\frac{5}{3} & \frac{3}{4} & - \frac{7}{6} & \frac{1}{4} & - \frac{5}{6} & - \frac{3}{4} & \frac{1}{3} & - \frac{1}{4}\\\frac{3}{4} & \frac{5}{3} & - \frac{1}{4} & \frac{1}{3} & - \frac{3}{4} & - \frac{5}{6} & \frac{1}{4} & - \frac{7}{6}\\- \frac{7}{6} & - \frac{1}{4} & \frac{5}{3} & - \frac{3}{4} & \frac{1}{3} & \frac{1}{4} & - \frac{5}{6} & \frac{3}{4}\\\frac{1}{4} & \frac{1}{3} & - \frac{3}{4} & \frac{5}{3} & - \frac{1}{4} & - \frac{7}{6} & \frac{3}{4} & - \frac{5}{6}\\- \frac{5}{6} & - \frac{3}{4} & \frac{1}{3} & - \frac{1}{4} & \frac{5}{3} & \frac{3}{4} & - \frac{7}{6} & \frac{1}{4}\\- \frac{3}{4} & - \frac{5}{6} & \frac{1}{4} & - \frac{7}{6} & \frac{3}{4} & \frac{5}{3} & - \frac{1}{4} & \frac{1}{3}\\\frac{1}{3} & \frac{1}{4} & - \frac{5}{6} & \frac{3}{4} & - \frac{7}{6} & - \frac{1}{4} & \frac{5}{3} & - \frac{3}{4}\\- \frac{1}{4} & - \frac{7}{6} & \frac{3}{4} & - \frac{5}{6} & \frac{1}{4} & \frac{1}{3} & - \frac{3}{4} & \frac{5}{3}\end{matrix}\right]$$

In [16]:
M.subs(rho, 1)


Out[16]:
$$\left[\begin{matrix}\frac{4}{9} & 0 & \frac{2}{9} & 0 & \frac{1}{9} & 0 & \frac{2}{9} & 0\\0 & \frac{4}{9} & 0 & \frac{2}{9} & 0 & \frac{1}{9} & 0 & \frac{2}{9}\\\frac{2}{9} & 0 & \frac{4}{9} & 0 & \frac{2}{9} & 0 & \frac{1}{9} & 0\\0 & \frac{2}{9} & 0 & \frac{4}{9} & 0 & \frac{2}{9} & 0 & \frac{1}{9}\\\frac{1}{9} & 0 & \frac{2}{9} & 0 & \frac{4}{9} & 0 & \frac{2}{9} & 0\\0 & \frac{1}{9} & 0 & \frac{2}{9} & 0 & \frac{4}{9} & 0 & \frac{2}{9}\\\frac{2}{9} & 0 & \frac{1}{9} & 0 & \frac{2}{9} & 0 & \frac{4}{9} & 0\\0 & \frac{2}{9} & 0 & \frac{1}{9} & 0 & \frac{2}{9} & 0 & \frac{4}{9}\end{matrix}\right]$$

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


In [17]:
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


In [18]:
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)

In [19]:
K_num = stiff_num(3)
K_num - K


Out[19]:
$$\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]$$

Best approach

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


In [20]:
from sympy.integrals.quadrature import gauss_legendre

In [21]:
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


In [22]:
from IPython.core.display import HTML
def css_styling():
    styles = open('../styles/custom_barba.css', 'r').read()
    return HTML(styles)
css_styling()


Out[22]: