# Computation of stiffness matrix with variable coefficient

I am implementing a finite element solver (in 2D) to solve the generic differential equation :

$$-\nabla(a(x) \nabla u) = f$$

# Brief explanation

By integrating and multipling by a test function, the weak form allows to integrate the stiffness matrix :

$$A_{ij} = \int a(x) \nabla \phi_i(x) \nabla \phi_j(x)\, dx\quad i=1,2,3\ .$$

More simply, we can write the local stiffness matrix (of a 3 node linear triangle) of the element $K$ (the assembly will be done in a next step):

$$A_{ij}^{K} = \int a(x) \nabla \phi_i(x) \nabla \phi_j(x)\, dx\quad i=1,2,3\ .$$

In first place, I use the central quadrature formula that allows to compute:

$$A_{ij}^{K} = a_\text{mean}(b_i b_j + c_i c_j) K_\text{area}$$

with $a_\text{mean}$ = central value of $a(x)$ in the center of the triangle. The coefficient are the one form the shape function: $\phi_i = a_i + b_i x1 + c_i x_2$

For one single local stiffness matrix, it represents this code in Python:

loc2glob = t[0:3,K] #local to global map
x = p[0,loc2glob] #get a list of the x-coordinates
y = p[1,loc2glob] #get a list of the y-coordinates
xc = np.mean(x), yc = np.mean(y)
a_centroid = a(xc,yc) #simplification
b = np.atleast_2d(b)
c = np.atleast_2d(c)
AK = (np.dot(b.T, b) + np.dot(c.T, c))*a_centroid*area

# And ...?

After the verification of this method, I end up with implementing a more general integration method (based on quadrature). The coordinates $x,y$ are transformed in $r,s$ (local integration, to have a general method). The Jacobian is computed numerically.

$$A_{ij}^{K} = \int_{el_K} a(r, s) \nabla \phi_i (r, s) \nabla \phi_i(r,s) |J_K(r,s)|\ dr\ ds$$

# A little more explanation on $(r,s)$ coordinate

In the domain of finite element method, computer implementation are used to switch the coordinate $(x,y)$ of a (triangle) element to an $(r,s)$ orthogonal base. The principle is very nice, it allows to do a variable change in the integral, and to have an easiest method to apply the quadrature :

$$\int_{\Omega_K} q(x,y) \phi_i(x,y)\phi_j(x,y)\ dx\ dy = \int_{\triangle} q(\xi, \eta) \psi_i(\xi, \eta)\psi_j(\xi, \eta) \left|\frac{\partial (x, y)}{\partial (\xi, \eta)}\right|\ d\xi\ d\eta$$

Graphically, it corresponds to :

So you can use the coordinate $(r,s)$ of a generic triangle for the integral, providing the Jacobian of the coordinate transformation.

In practice, I use a 3 gauss point method. The basics are here:

AK = np.zeros((3,3))
for q in range(len(qwgts)):
r = rspts[q,0] # r coordinate of the q_th quadrature point
s = rspts[q,1] # s coordinate of the q_th quadrature point
a_gauss_point = a(r,s)
S,dSdx,dSdy,detJ = isopmap(x,y,r,s,shapefcn)
dSdx =np.matrix(dSdx),dSdy =np.matrix(dSdy)
wxarea = a_gauss_point*qwgts[q]*detJ/2 #weight * area (corrected by 0.5)
AK = AK + (dSdx.T.dot(dSdx) + dSdy.T.dot(dSdy))* wxarea

I know that the method "isopmap" return a correct Jacobian, since I have tested with data found in the literature.

# My problem...

When I take $a(x,y) = 1$ , the two methods give the same result. But when I put $a(x,y) = x + y$ for example, I end with significant difference between the two matrix.

# And so?

Do you have any explanation? Of course the problem is about the value $a(r,s)$ in the second code, but I don't see where the thing fails...

# SOLUTION : EDIT (17.07.2017)

In fact it comes that I didn't understand properly the logic of the mapping (x,y) domain to (r,s) domain. In fact the integral can be written like this:

Then I just have to change in my python code the evaluation of the function $a(x,y)$:

x_physical = np.dot(x,S)
y_physical = np.dot(y,S)
a_gauss_point = a(x_physical,y_physical)

# If you want to reproduce this behaviour (EDIT : correct code uploaded)

You can download the Python coded that contains 100 lines of code in order to compute very easily the $A_K$ matrix.

The following is the code

from __future__ import division # avoid integer problem of division

#Import zone
import numpy as np

#Define shape function
def P1shapes(r,s):
S = np.array([1-r-s,r,s])
dSdr = np.array([-1,1,0])
dSds = np.array([-1,0,1])
return S,dSdr,dSds
#Jacobian function
def isopmap(x,y,r,s,shapefcn):
# x = vector of x coordinate of the element's point
#shapefcn = P1shapes or P2shapes
S,dSdr,dSds = shapefcn(r,s);
j11=np.dot(dSdr,x)
j12=np.dot(dSdr,y)
j21=np.dot(dSds,x)
j22=np.dot(dSds,y)
detJ=j11*j22-j12*j21
dSdx=( j22*dSdr-j12*dSds)/detJ
dSdy=(-j21*dSdr+j11*dSds)/detJ
return S,dSdx,dSdy,detJ

#for linear shape function (1st order)
area = polyArea(x,y)
b = np.array([y[1]-y[2],y[2]-y[0],y[0]-y[1]])/2/area
c = np.array([x[2]-x[1],x[0]-x[2],x[1]-x[0]])/2/area
return area,b,c

def polyArea(x,y):
#Shoelace formula => Compute Area of Polygons, x and y set of pointd given (vertices coordinates)
return 0.5*np.abs(np.dot(x,np.roll(y,1))-np.dot(y,np.roll(x,1)))

#Define a(x,y)
def a(x,y):
return x + y

qwgts=np.array([1/3,1/3,1/3])
rspts=np.array([[1/6,1/6],
[2/3,1/6],
[1/6,2/3]])

#Point of the 3 nodes triangles (test case)
x = [-21.68467035,-21.17695462 ,-22.18700401]
y = [-9.94204652 ,-8.91258056 ,-9.12362242]

#First method (simplifed one)
xc = np.mean(x)
yc = np.mean(y)
a_centroid = a(xc,yc) #simplification
b = np.atleast_2d(b)
c = np.atleast_2d(c)
AK_simplified = (np.dot(b.T, b) + np.dot(c.T, c))*a_centroid*area

#Second Method
for q in range(len(qwgts)):
r = rspts[q,0] # r coordinate of the q_th quadrature point
s = rspts[q,1] # s coordinate of the q_th quadrature point
S,dSdx,dSdy,detJ = isopmap(x,y,r,s,P1shapes)
#Map the gauss point to the physical domain (r,s)_gauss to (x,y)_gauss
x_physical = np.dot(x,S)
y_physical = np.dot(y,S)
a_gauss_point = a(x_physical,y_physical)

dSdx =np.matrix(dSdx)
dSdy =np.matrix(dSdy)
wxarea = a_gauss_point*qwgts[q]*detJ/2 #weight * area (corrected by 0.5)

#Show the matrix
print AK_simplified
print " "
• You say: "use a 2 gauss point method". Do you mean 3 gauss points? Why are you trying to pass (r,s) to a function that needs (x,y)? – Bill Greene Jul 16 '17 at 20:57
• Hi! Thank you for your reply! Yes I mean 3 gauss points, I have done the correction thanks :D I have added some explanation why the integral is done with the tupple (r,s), and not the coordinate (x,y). basically, it's a coordinate transformation – Electrode Jul 16 '17 at 21:19
• I do not see how you can implement the function $a(x,y)$ you show above by passing $r,s$. You need to do that coordinate transformation from $r,s$ to $x,y$ inside your integration point loop and pass those $x,y$ values to your $a$ function. – Bill Greene Jul 16 '17 at 21:33
• Are your elements linear-elements?, i.e., three-nodes elements? – nicoguaro Jul 16 '17 at 23:36
• Yes @nicoguaro, the elements are linear triangle (with 3 node). So they have no quadratic shape function – Electrode Jul 17 '17 at 9:43

If an isoparametric element formulation with numerical integration is used, the user-written property function will be called at each integration point. Given the parametric coordinates of the integration point, say $r$ and $s$ in two dimensions, the required global $x$ and $y$ coordinates are obtained from the basic isoparametric definition
$$x = \sum_{i=1}^{n_s} N_i(r,s) x_i$$ $$y = \sum_{i=1}^{n_s} N_i(r,s) y_i$$ where $N_i$ are the shape functions, $x_i$ and $y_i$ are the global coordinates of the element nodes, and $n_s$ is the number of nodes in the element.