Uncategorized

Integrals Over Arbitrary Triangular Regions for FEM

1 Comment

In this post I’ll present a recipe for taking an integral over an arbitrary triangular region using the SciPy integrate.dblquad() function. This is an important operation for implementing the Finite Elements method for solving partial differential equations. In school we are taught to perform a change of variables which involves splitting the triangle into two regions and performing the double integration on the simpler sub-domains after carefully calculating new limits of integration. This recipe maps the triangle to the unit square, and then calculates the double integral on the domain [0,1] \times [0,1]. I pieced this together after looking at this discussion on the MATLAB Central message board regarding the transformation of the triangle to the unit square, and this post on Paul’s Online Notes that touched on the calculation of the Jacobian, and this post by John D. Cook about choosing the correct error limits for quadrature integration.

The Jacobian

Suppose we have a transformation x=g(u,v) and y=h(u,v) for the independent variables of the function f(x,y). The Jacobian of this transformation is given by,

(1)   \begin{equation*} \dfrac{\partial (x,y)}{\partial (u,v)} = \mbox{det}\begin{vmatrix} \frac{\partial x}{\partial u} & \frac{\partial x}{\partial v} \\ \frac{\partial y}{\partial u} & \frac{\partial y}{\partial v} \end{vmatrix} \end{equation*}

Our transformation is given by,

(2)   \begin{align*} x &= g(u,v) = (1-u)x_{1} + u \bigl( (1-v)x_{2} + vx_{3} \bigr) \\ y &= h(u,v) = (1-u)y_{1} + u \bigl( (1-v)y_{2} + vy_{3} \bigr) \end{align*}

where p_{i}=(x_{i},y_{i}) are the three points that define the triangle. Then,

(3)   \begin{align*} \dfrac{\partial x}{\partial u} & = -x_{1} + (1-v)x_{2} + vx_{3} \\ \dfrac{\partial x}{\partial v} & = -ux_{2} + ux_{3} \\ \dfrac{\partial y}{\partial u} & = -y_{1} + (1-v)y_{2} + vy_{3} \\ \dfrac{\partial y}{\partial v} & = -uy_{2} + uy_{3} \\ \end{align*}

and so the determinant of the Jacobian of this transform becomes,

(4)   \begin{equation*} det \biggl( \dfrac{\partial(x,y)}{\partial(u,v)} \biggr) = \frac{\partial x}{\partial u}\frac{\partial y}{\partial v} - \frac{\partial x}{\partial v}\frac{\partial y}{\partial u} \\ \end{equation*}

If we resist the urge to write a lambda function, we can write this in a nice, almost descriptive Python function.

def dJ( u, v, p1, p2, p3 ):
    x1, y1 = p1
    x2, y2 = p2
    x3, y3 = p3
    dxdu = ( (1-v)*x2 + v*x3 - x1 )
    dxdv = ( u*x3 - u*x2 )
    dydu = ( (1-v)*y2 + v*y3 - y1 )
    dydv = ( u*y3 - u*y2 )
    return np.abs( dxdu*dydv - dxdv*dydu )

Integration

Once we transform the function over the triangular region to the unit square, we need to multiply by the absolute value of the determinant of the Jacobian of the transformation. We’ll perform the integration using the scipy.integrate function dblquad(), which is short for double quadrature. The first argument is the integrand, the next two arguments are the limits of the outer integral, and the next two arguments are functions for the limits of the inner integral. Since we’re integrating over the unit square, the limits of the inner integral will be given as lambda x: 0 and lambda x: 1, respectively. There are two more arguments specifying limits for the estimates of the relative and absolute error.

We can wrap all of this up into a function that takes the integrand and the points of the triangle, and then performs the transformation and the integration for us, saving us from working out a change of variables.

def tridblquad( integrand, p1, p2, p3 ):
    '''
    Perform double quadtrature integration on triangular domain.
    Input: function to integrate, points of triangle as tuples.
    Output: integral and estimated absolute error as a tuple.
    '''
    x1, y1 = p1 ; x2, y2 = p2 ; x3, y3 = p3
    # transformation to the unit square
    g = lambda u, v, c1, c2, c3: (1-u)*c1 + u*( (1-v)*c2 + v*c3 )
    # transformation for the integrand, 
    # including the Jacobian scaling factor
    def h( u, v ):
        x = g( u, v, x1, x2, x3 )
        y = g( u, v, y1, y2, y3 )
        I = integrand( x, y )
        I *= dJ( u, v, p1, p2, p3 )
        return I
    # perfrom the double integration using quadrature in the transformed space
    integral, error = scipy.integrate.dblquad( h, 0, 1, lambda x: 0, lambda x: 1, epsrel=1e-6, epsabs=0 )
    return integral, error

As a quick sanity-check, we can calculate the area of the (unit?) triangle with the vertices (0,0)-(0,1)-(1,0), which should be 0.5.

area, error = tridblquad( lambda x, y: 1., (0.,0.), (1.,0.), (0.,1.) )
print 'Area: {},   Error: {:.3e}'.format( area, error )

Area: 0.5,   Error: 5.551e-15

One thought on “Integrals Over Arbitrary Triangular Regions for FEM”

  1. Apart from mapping the triangle onto a quadrilateral, there are other ways of performing quadrature on the triangle too. Many publications are devoted to this topic; check out quadrature for a list of quadrature schemes for the triangle and their implementation.

Comments are closed.