Notebook

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¶

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

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]: