Scientific Programming

Volume 2015 (2015), Article ID 797325, 22 pages

http://dx.doi.org/10.1155/2015/797325

## Finite Element Assembly Using an Embedded Domain Specific Language

^{1}Department of Mechanics, Royal Military Academy, Avenue de Renaissance 30, 1000 Brussels, Belgium^{2}von Karman Institute for Fluid Dynamics, Chaussée de Waterloo 72, 1640 Rhode-Saint-Genèse, Belgium^{3}LaSIE, La Rochelle University, Avenue Michel Crépeau, 17042 La Rochelle Cedex 1, France

Received 18 November 2013; Revised 22 October 2014; Accepted 19 January 2015

Academic Editor: Bormin Huang

Copyright © 2015 Bart Janssens et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

In finite element methods, numerical simulation of the problem requires the generation of a linear system based on an integral form of a problem. Using C++ meta-programming techniques, a method is developed that allows writing code that stays close to the mathematical formulation. We explain the specifics of our method, which relies on the Boost.Proto framework to simplify the evaluation of our language. Some practical examples are elaborated, together with an analysis of the performance. The abstraction overhead is quantified using benchmarks.

#### 1. Introduction

The application of the finite element method (FEM) requires the discretization of the integral form of the partial differential equations that govern the physics, thus transforming the problem into a set of algebraic equations that can be solved numerically. In essence, the resulting numerical model only depends on the governing equations and the set of basis and test functions. From the point of view of a model developer, it would therefore be ideal to only need to specify these parameters in a concise form that closely resembles the mathematical formulation of the problem, without sacrificing computational efficiency.

The current work is part of Coolfluid 3 [1], a C++ framework intended primarily for the solution of fluid dynamics problems and available as open source under the LGPL v3 license. The system is designed around a Component class that can provide a dynamic interface through Python scripting and a GUI. Model developers will typically write a set of components to build a solver for a given set of equations, using functionality provided by our framework and external libraries. Problem-dependent settings such as the mesh, boundary conditions, and model parameters can then be controlled through Python or a GUI, though these aspects are beyond the scope of the current paper. In this context, we provide an Embedded Domain Specific Language (EDSL) that can be used to implement finite element solver components. It uses a notation similar to the mathematical formulation of the weak form of a problem, allowing the programmer to focus on the physics rather than coding details. We assume a typical finite element workflow, where a global linear system is assembled from element matrix contributions. Our language can also describe boundary conditions and field arithmetic operations.

All language elements consist of standard C++ code, so they easily embed into the components of our framework. An extension mechanism is available, allowing any developer to implement his own language elements without touching the library code. The implementation uses the Boost.Proto [2] framework, which builds on template meta programming techniques for generating efficient code. The actual performance overhead will be quantified and compared to the global solution time for some applications.

The automatic generation of code based on an intuitive specification of the physical problem is of course a feature that is highly desirable for any numerical simulation framework. One example of previous work providing such functionality is OpenFOAM [3], a C++ library primarily intended for fluid mechanics. It allows for the easy expression of differential equations by providing an embedded tensor manipulation language, specific to finite volume discretizations. The FEniCS project [4] provides an easy and efficient method for developing finite element discretizations. The approach differs from ours in that the language is embedded into Python and compiled using a just-in-time compiler. It also offers the option of automated analytical evaluation of integrals [5]. Another example is the FEEL++ project [6], which also provides a domain specific language embedded into C++ using expression templates.

Notwithstanding the existence of all these excellent projects, we decided to implement our own framework for tight integration with the Coolfluid data structures and to support our workflow of building user-configurable components. The implementation relies on expression templates, just like FEEL++, but we use Boost.Proto to handle the expression template generation and to define a grammar for our language. This results in simpler code and has a limited impact on runtime performance, as we will show. Even though our language is not as feature-complete as the more established FEEL++ and FEniCS projects, we do believe that the easy integration of user-defined extensions into the language is a unique feature. The purpose of the current paper is to present how Proto can be used to construct a language for finite element modeling—based on a simple example building a mass matrix—and to show the capabilities and performance of the language we developed. Aspects that are new compared to previous work in the field are the use of Proto as a new way to construct the language and the possibility to add user-defined terminals to the language.

This paper is structured as follows. In Section 2 the mathematical formulation of the integral form will be laid out in order to clearly define the class of problems that is supported. Next, in Section 3 the mechanics for constructing and interpreting the EDSL are explained in detail, based on a simple application. Section 4 contains some application examples and Section 5 a performance analysis. Finally, in Section 6 we present our conclusions and suggest future work.

#### 2. Finite Element Discretization

In this section, we introduce the notation used for finite element discretizations, starting from the Poisson problem as an example. A more generic and much more thorough introduction to the finite element method is given in, for example, [7]. Considering a domain with boundary , the differential equations describing the Poisson problem are

To transform these continuous equations into a discrete problem, the unknown is interpolated using a linear combination of shape functions with the unknown coefficients :

The shape functions depend only on the spatial coordinates and to simplify the notation we will from here on just write instead of . The shape functions have a local support and the interpolation happens on a per-element basis, so is equal to the number of nodes in an element. The interpolation can also be written in vector form as . Here, is a row vector with the shape functions associated with unknown and is the column vector of the unknown coefficients for the element.

To obtain a linear system for the discrete problem, we need to multiply the equations with a weighting function and integrate over the domain. Using the Galerkin method (i.e., the weighting functions are the shape functions) as well as integrating per-element yields the following linear system for the Poisson problem, using the weak formulation here:

Due to the local support of the shape functions, the integral is written as a sum of integrals over all the elements. Each element integral results in an element matrix on the left hand side and an element vector on the right hand side. Summing up the contributions of all elements that share a given node assembles a global linear system that has the total number of nodes as dimension.

The indices are useful for problems with multiple unknowns. In this case, each unknown can be associated with its own shape function, and the element matrices and vectors are built up of blocks associated with each variable and equation.

The matrix assembly procedure is the same for all problems; only the values of the element matrices depend on the initial equation (1). This observation drives the current work: we will present a convenient language to write out the element equations and automate the steps that are common for all finite element discretizations. The code for the Poisson problem is presented in Listing 1, where lines and map directly to the element matrix and vector, as defined in (4), and lines and represent the assembly into the global system. Lines and form the core of our language, and the techniques for implementing this will be explained in detail in the next section. The remainder of Listing 1 will be further explained in Section 4.1.