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
area,b,c = gradient(x,y) # compute the numerical gradient
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);
    dSdx=( j22*dSdr-j12*dSds)/detJ
    return S,dSdx,dSdy,detJ

#Gradient Function
def gradient(x,y):
    #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

#Gauss point quadrature

#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)
area,b,c = gradient(x,y) # compute the numerical gradient
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
AK_quadrature = 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
    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)
    AK_quadrature = AK_quadrature + (dSdx.T.dot(dSdx) + dSdy.T.dot(dSdy))* wxarea

#Show the matrix
print AK_simplified
print " "
print AK_quadrature
  • $\begingroup$ 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)? $\endgroup$ Jul 16, 2017 at 20:57
  • $\begingroup$ 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 $\endgroup$
    – Electrode
    Jul 16, 2017 at 21:19
  • 2
    $\begingroup$ 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. $\endgroup$ Jul 16, 2017 at 21:33
  • $\begingroup$ Are your elements linear-elements?, i.e., three-nodes elements? $\endgroup$
    – nicoguaro
    Jul 16, 2017 at 23:36
  • $\begingroup$ Yes @nicoguaro, the elements are linear triangle (with 3 node). So they have no quadratic shape function $\endgroup$
    – Electrode
    Jul 17, 2017 at 9:43

1 Answer 1


It is not uncommon in a finite element model to have material properties or boundary conditions that vary as a function of the spatial coordinates. Often these are implemented as a "user-written" function that takes the coordinates as input values and returns the properties or boundary conditions at these coordinates. The input coordinates are typically defined in the global Cartesian system.

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.

  • $\begingroup$ Thank you! It is exactly what I have corrected in my code, and now it is working :D $\endgroup$
    – Electrode
    Jul 17, 2017 at 12:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.