# Stiffness matrix computation for 4 node quadrilateral element

I am writing a finite element code for heat transfer (scalar field problem) and starting from simple 4 node quadrilateral element. I tried computing conductance (stiffness) matrix in the physical coordinate systems and comparing the answer with isoparametric system.

Element looks like:

 |Y
|
4 _____________3
|             |
|             |
|             |
|             |
|             |
|_____________|  ---->X
1             2


(Side 1-2 is not straight.) Coordinates of nodes 1,2,3,4 are :

1: (0,0)
2: (2,0.5)
3: (2,1)
4: (0,1)


Then I make isoparametric transformation to bi-unit (e,n):

       |n
|
4 _____________3
|     |       |
|     |       |
|     |------ | ----->e
|             |
|             |
|_____________|
1             2


So, Coordinates of nodes 1,2,3,4 become :

1: (-1,-1)
2: (1,-1)
3: (1,1)
4: (-1,1)


I created a python simpy script to compute stiffness matrix.

import sympy as sp
import numpy as np
from sympy import init_printing

n,e = sp.symbols('n e')
# shape functions in isoparametric
N1_iso = (1/4)*(e-1)*(n-1)
N2_iso = (1/4)*(e+1)*(-n+1)
N3_iso = (1/4)*(e+1)*(n+1)
N4_iso = (1/4)*(-e+1)*(n+1)

#gradient of shape function matrix in isoparametric
GN_mat = sp.Matrix([ [sp.diff(N1_iso,e),sp.diff(N2_iso,e), sp.diff(N3_iso,e), sp.diff(N4_iso,e)],
[sp.diff(N1_iso,n),sp.diff(N2_iso,n), sp.diff(N3_iso,n), sp.diff(N4_iso,n)]])

#point matrix
P_mat = sp.Matrix([[0.0,0.0],[2.0,0.5],[2.0,1.0],[0.0,1.0]])
#Jacobian matrix
J_mat = (1.0/4.0)*(GN_mat*P_mat)

#Jacobian matrix determinant and inverse
J_det=J_mat.det()
J_inv=J_mat.inv()

#"B" matrix, assuming unit conductivity matrix
B_mat=J_inv*GN_mat

#conductance (stiffness) matrix
expr1 = B_mat.T*B_mat*J_det

#Double integrating conductance (stiffness) matrix
expr2 = sp.integrate(expr1,(n,-1,1))
#Stiffness (conductance) matrix using isoparametric transformation
K_mat=sp.integrate( (expr2),(e,-1,1))
print('K_mat = \n {}'.format(K_mat ))

#Now in physical coordinate system
x,y = sp.symbols('x y')
#Area of the element = 2
A=1.5

#shape function
N1 = (1/A)* (x-2)*(y-1)
N2 = (1/A)* (x-0)*(1-y)
N3 = (1/A)* (x-0)*(y-0.5)
N4 = (1/A)* (2-x)*(y-0)

B_mat2 = sp.Matrix([[sp.diff(N1,x), sp.diff(N2,x), sp.diff(N3,x), sp.diff(N4,x)],
[sp.diff(N1,y), sp.diff(N2,y), sp.diff(N3,y), sp.diff(N4,y)]])

expr10 = sp.integrate(B_mat2.T*B_mat2,(y,0,1))
#Stiffness (conductance) matrix using physical coordinates
K_mat2 = sp.integrate(expr10,(x,0,2))
print('     K_mat2 = \n {}'.format(K_mat2))


If I make point 2: (2,0.0), domain becomes rectangle and I get same answers. But, if it is not the case, I get different answers for K_mat and K_mat2. Any suggestions what is going wrong?

• Could you clarify your question? What is K_mat and what K_mat2? – nicoguaro Nov 14 '17 at 2:20
• I have comments in the code snippet, but anyway edited the question. I got K_mat by isoparametric transformation and K_mat2 without using it - just using physical coordinate system and its shape functions. – vcx34178 Nov 14 '17 at 9:45
• @v34178 At the moment, your big problem is that your shape functions don't seem to be consistent with your transformation. As a side note, you also need to think about the limits of your integrations in the physical case. – origimbo Nov 14 '17 at 10:58
• @origimbo, I got your suggestion about integration limits in physical domain. Can you explain what you mean by shape functions don't seem to be consistent with your transformation? I found that shape functions, in physical and isoparametric domains, match to nodes (node position and numbering in the shown figure). Appreciate your feedback. – vcx34178 Nov 15 '17 at 1:51

As I'm sure you're already aware, the Galerkin finite element method expands both test and trial function spaces in terms of a finite basis, $$T=\sum_i T_i N_i(\mathbf{x}).$$ Hopefully you're also aware that different choices of the $N$s will give different discretizations, and thus different matrix problems.
The isoparametric bilinear element for quadrilaterals, Q4, is defined by the shape functions you give for the reference quadrilateral (your $N^\mbox{iso}$), along with a spatial transformation from the reference to physical space, $$x(e,n) = x_1 N_1(e,n) + x_2 N_2(e,n) +x_3 N_3(e,n) + x_4 N_4(e,n),$$ $$y(e,n) = y_1 N_1(e,n) + y_2 N_2(e,n) +y_3 N_3(e,n) + y_4 N_4(e,n).$$ This means that for a given element we can trivially write $$N_i(x(e,n),y(e,n)) = N_i(e,n),$$ but to give an explicit definition in terms of $x$ and $y$ we have to invert the relation above.
At this point it should be obvious that you've chosen a different physical basis, since (1) your physical $N_1$ shape function isn't equal to unity at vertex 1 and (2) your shape functions don't all vanish on the edges not defined by "their" vertex. The last point is important, since it's what guarantees that the Q4 representation is continuous across elements. The reason that your code works when you change the vertex position is that you then have a physical quadrilateral which is a linear translation of the reference quadrilateral. In that case, the (trivial) invertion does give back what you've written.