INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING
Int. J. Numer. Meth. Engng 2005; 00:1–6
Prepared using nmeauth.cls [Version: 2002/09/18 v2.02]

Efficient Visualization of High Order Finite Elements
Jean-François Remacle1∗ , Nicolas Chevaugeon1 , Émilie Marchandise1,3 , Christophe
Geuzaine2,3
1 Département de Génie Civil et Environnemental, Université Catholique de Louvain, Bâtiment Vinci, Place

du Levant 1, B-1348 Louvain-la-Neuve, Belgium
2 Department of Applied Mathematics, California Institute of Technology, Pasadena CA
3 Fonds National de la Recherche Scientifique, rue d’Egmond, Bruxelles.

SUMMARY
A general method for the post-processing treatment of high order finite element fields is presented.
The method applies to general polynomial fields, including discontinuous finite element fields. The
technique uses error estimation and h-refinement to provide some optimal visualization grid. Somefiliering is added to the algorithm in order to focus the refinement on a visualization plane, one on thecomputation of one single iso-zero surface. Some 2D and 3D examples are provided that illustrate thepower of the technique. In addition, schemes and algorithms that are discussed in the paper are readily
available as part of an open source that is developped by the authors, namely Gmsh. Copyright c
2005 John Wiley & Sons, Ltd.
key words:

High order finite elements, Discontinuous Galerkin, Scientific visualization

1. Introduction
In the recent years, a large research effort has been devoted to the development of highorder finite element discretization techniques: spectral element methods [6], high order finite
elements [3, 2] or high order Discontinuous Galerkin Methods [12, 1, 5]. Our focus is high order
Discontinuous Galerkin Methods (DGM). Those are now used extensively for solving transient
aeroacoustics problems [10], electromagnetic problems [14], dynamics of compressible fluids[13, 12], Kortenweg-de Vries (KdV) equations [15] and many other relevant physical problems.
Using a quadrature free approach together with a very careful BLAS3 implementation, ourexperience show that, with 5th or 6th order polynomials on a tetrahedral mesh, we can run
a DGM code for which we use about 80 percent of the peek performance of an off-the-shelfdesktop computer. This makes high order DGM a very competitive approach.
In contrast, most of the finite element post-processing tools are unable to provide accurate
high order visualizations. Typically, quadratic elements are the highest available order for∗ Correspondence to: J.F. Remacle, Département de Génie Civil et Evironnemental, Université Catholique deLouvain, Bâtiment Vinci, Place du Levant 1, B-1348 Louvain-la-Neuve, Belgium.
† Email: remacle@gce.ucl.ac.be

Copyright c 2005 John Wiley & Sons, Ltd.

Received
Revised2

J.F. REMACLE, N. CHEVAUGEON E. MARCHANDISE C. GEUZAINE

visualization. In this paper, we propose a robust and efficient methodology for the visualizationof high order finite element solutions. This work includes the issues of contouring, iso-surfacing
and cutting. The methodology is inspired by h-adaptive finite element techniques [12, 13].
In all what follows, computer efficiency is always a central concern and we will be always
careful when it comes to cpu time issues. When plots and images of 2D and 3D results will beshown, we will always indicate computation time and number of rendered polygons.
Another concern that we have is compatibility. We think that it is the post-processor thathas to adapt to the solver and not the inverse. In our developments, it is the solver that dictates
how fields are interpolated and how elements are mapped. This means that we should be able
to visualize any kind of finite element solutions, continuous or not, using any sort of shapefunctions.
Finally, we think that providing some kind of source code is essential in this kind of business.
An implementation of our methodology is readily available as an open source code. Gmsh[8]is an automatic 3D finite element grid generator (primarily Delaunay) with a build-in CAD
engine and post-processor. Its design goal is to provide a simple meshing tool for academic
problems with parametric input and up to date visualization capabilities. Gmshis copyrightc 1997-2004 by C. Geuzaine and J.-F. Remacle and is distributed under the terms of the
GNU General Public License (GPL). The high order visualizations techniques described inthis paper have been added to Gmsh. All the illustrations shown below have been done withGmsh.

2. High order discontinuous finite elements
Finite element methods (FEMs) involve a double discretization. First, the physical domain Ω
is discretized into a collection of Ne elements
Te =

Ne
[

e

(1)

e=1

called a mesh. This first step is the one of geometrical discretization. Then, the continuousfunction spaces (infinite dimensional) are replaced by finite dimensional expansions. The
difference between the DGM and classical Finite Element Methods (FEMs) is that the solution
is approximated in each element separately: No a priori continuity requirements are needed.
The discrete solution may then be discontinuous at inter-element boundaries. Figure 1 showsa typical situation of three elements e1 , e2 and e3 . The approximated field u is smooth in eachelement but may be discontinuous at inter-element boundaries.
Consider a scalar field u(x, y, z). In each element, it is usually polynomial spaces that arechosen for approximating u. The approximation of u over element e, noted ue is written
ue =

d
X

φj (ξ, η, ζ)U (e, j)j=1

where the U (e, j) are the coefficients of the approximation or degrees of freedom, where d isthe size of the discrete function space and where φj (ξ, η, ζ) is the jth approximation function
defined in a common parametric space ξ, η, ζ. In each element, we have d coefficients, and,
Copyright c 2005 John Wiley & Sons, Ltd.
Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2005; 00:1–6

EFFICIENT VISUALIZATION OF HIGH ORDER FINITE ELEMENTS

3

u|e1
u|e3

u|e2

u

e1

e3
e2

y
x
Figure 1. Three elements e1 , e2 and e3 and the piecewise discontinuous solution u.

because we consider (this is the most general case) that all approximations are disconnected,
Ne ×d coefficients are required for describing u. In our implementation, we organize those data
in a Ne × d matrix that we call U . Uij is therefore the jth component of ui . We use the GNU
scientific library (GSL) for manipulating matrices and vectors. The GSL have the advantage
to provides direct and easy bindings to the Basic Linear Algebra Subroutines (BLAS).
The element e itself may be geometrically high order. We define the geometrical mappingas the mapping that transforms a reference element into the real element:xe = xe (ξ, η, ζ) , y e = y e (ξ, η, ζ) , z e = z e (ξ, η, ζ),
with
xe =

m
X

ψj (ξ, η, ζ)X(e, j),mXψj (ξ, η, ζ)Y (e, j),mXψj (ξ, η, ζ)Z(e, j).

j=1

e

y =

j=1

ze =

j=1

Matrices X, Y and Z are of size m × 3.
Our goal is to be able to treat (i.e. visualize) any piecewise polynomial approximation. For
that purpose, our visualization methodology has to• enable the represent of discontinuous fields,• distinguish geometric and functional discretization,
Copyright c 2005 John Wiley & Sons, Ltd.
Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2005; 00:1–6

4

J.F. REMACLE, N. CHEVAUGEON E. MARCHANDISE C. GEUZAINE• be able to use any set of polynomial interpolation function.
The fact that we have disconnected interpolations by defining as much coefficients as N e × dallow naturally to represent discontinuous fields. The fact that we consider the general case of
non-isoparametric elements, which means that d 6= m and φi 6= ψi , enables us to distinguishgeometric and functional discretizations.
Here, we only consider polynomial approximations. Usual FEM have limited choices for the
φi ’s due to the a priori continuity requirements of the approximation. Nodal, hierarchical,
Serendip or non-conforming basis are among the usual basis for classical FEMs. In the generalcase of a DGM, there are no limitations for the choice of the φi ’s. For giving the a generalform to our φi ’s, we define the two following matrices. The first matrix, named P , is of sized × 3. Matrix P allow the definition of a canonical polynomial basis:
P = {p1 (ξ, η, ζ), . . . , pd (ξ, η, ζ)}
wherepi = ξ Pi1 η Pi2 ζ Pi3 .
This basis is then used for building the φi ’s. For that, we define the coefficient matrix Φ ofsize d × d and we writedXφi =pk Φik .
k=1

For example, the a bilinear isoparametric quadrilateral element, we have




0 0 0
1 1
1
1
 1 0 0 1
 1 1 −1 −1 P =
 0 1 0  , Φ = 4  1 −1 −1 1  , P = {1, ξ, η, ξη} and
1 1 0
1 −1 1 −1
1{(1 + ξ)(1 + η), (1 + ξ)(1 − η), (1 − ξ)(1 − η), (1 − ξ)(1 + η)} .
4
This way of building the basis allow to define any polynomial basis, complete or not, of any
dimension. The advantage is that any FEM code can use the tool, eventually Gmsh, without
having to change its own representation to the one of the post-processor.
The same technique is used for describing geometrical mappings in the general form. A highorder finite element visualization structure is described in Algorithm 1. The DGM allow theuse of exotic polynomial approximations and are, therefore, a good support for testing ourvisualization methodology.
φ=

3. Contouring high order fields
One interresting approach to high order contouring can be found in [9]. This approach isbased on the definition of some kind of algebra of high order textures. Each basis function
φi correspond to one basis texture and the final visualization object is constructed by linear
combination of basis textures. Figure 2 show an example of this approach. Quadratic textures
s1 , s2 and s3 corresponding to hierarchical shape functions associated to the edges of a triangle
Copyright c 2005 John Wiley & Sons, Ltd.
Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2005; 00:1–6

EFFICIENT VISUALIZATION OF HIGH ORDER FINITE ELEMENTS

5

Algorithm 1 A general structure for the definition of high order finite element solution.
struct _high_order_fem {// coefficients for the element mappings
GSL_Matrix X,Y,Z;// coefficients for the function
GSL_Matrix U;// Canonical matrices for U and for the mapping
GSL_Matrix PU,PM;// Interpolation matrices
GSL_Matrix PHIU,PHIM;
};

s1

s2

s3

s1 + s 2 + s 3

Figure 2. Illustration of the visualization technique developed in [9] based on high order textures.

are shown. They are subsequently summed and the result s = s1 + s2 + s3 is computed as thegraphic sum of the basis textures. Disappointingly, this technique relies heavily on the graphic
hardware and very few off-the-shelf graphic cards allow to perform those kind of combinations
efficiently. Also, the most common graphical API, i.e. OpenGL, does not fully support linear
combination of textures.
More important, this approach has aspects that are way too restrictive for finite element
visualizations. The number of isocontours, for example, cannot be changed easily because basis
textures contain a fixed number of colors. Changing to a different scale, a logarithmic scale forexample, is difficult because the logarithm of the sum is not the sum of the logarithms. Also,
extracting iso-surfaces or plane cuts out of 3D data is not straightforward with this approach.
Here, we focus on what we believe to be a more convenient approach. Contouring high order
FEM fields is difficult because of the topological complexity of high order curves and surfaces.
If we consider 7th or 9th order approximations ue on a triangle, the shape of one iso-contour of
ue = C may be very complex (it could even be not connected). Only piecewise linear (p = 1)
interpolations allow straightforward representation: one iso-contour is a straight line, one isosurface is a polygon. We assume here that our visualization tool is able to deal with piecewise
linear simplices: lines, triangles and tetrahedra. It is also able to draw bi- and tri-linear fields
on quadrangles and hexahedra. In Gmsh, this is done by splitting quadrangles and hexahedra
into simplices and doing linear visualization.
One robust approach of high order visualization consist therefore in dividing the elements
into sub-elements and doing a linear visual approximation on every sub-element. The main
problem here is “how far do we have to divide in order to capture the complexity of the highCopyright c 2005 John Wiley & Sons, Ltd.
Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2005; 00:1–6

6

J.F. REMACLE, N. CHEVAUGEON E. MARCHANDISE C. GEUZAINE

order field”?. For addressing that issue, we use here h-refinement. The problem is made simpler
here by the fact that we know the exact visualization error i.e. the difference between u e andits “drawable” piecewise linear representation.

4. Visualization mesh
The technique we use here is based on classical AMR (Automatic Mesh Refinement)
methodologies, see for example [4, 12]. For each element type (triangle, quad, tetrahedron,hexahedron), we define a template for dividing the element itself recursively into sub-elementsof the same type. As an example, the refinement template is shown for the triangle on Figurer = 0, Nr = 3

r = 2, Nr = 15

r = 4, Nr = 153

Figure 3. Some AMR refinement templates for triangles at different resolution levels.

3 at different recursion levels r. The element sub-division pastern is performed in the reference
system of coordinates ξ, η, ζ. For a given maximal recursion level r, we obtain an array of N rvisualization points ξi , ηi and ζi with
Nr =

1 r
(2 + 1)(2r + 2)
2

for triangles,
Nr = (2r + 1)2
for quadrangles,
Nr =

1 r
(2 + 1)(2r + 2)(2r + 3)
6

for tetrahedra and
Nr = (2r + 1)3
for hexahedra. Using previous definitions, we construct the following interpolation matrix
Iij = φj (ξi , ηi , ζj )of size Nr × d. Matrix I is computed once and stored. The value of ue at visualization points
is computed efficiently using a simple matrix-vector productue (ξi , ηi , ζi ) = Iij U (e, j).
This product can be done efficiently using BLAS2 routines [7]. Typically, around one GigaFlop
can be obtained for p > 4 and r > 4 on a 2.4 GHz Pentium IV Xeon. Our experience has shown
Copyright c 2005 John Wiley & Sons, Ltd.
Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2005; 00:1–6

EFFICIENT VISUALIZATION OF HIGH ORDER FINITE ELEMENTS

7

us that it was faster to compute a large amount of points using fast linear algebra techniques
than only computing values when needed.
The number of triangles or quadrangles at recursion level r is 4r and the number oftetrahedron or hexahedron is 8r . All sub-elements should not always be visible: some kindof error analysis should be done in order to decide wether a given triangle is visible or not.
This decision has to be goal oriented i.e. it should depend on what has to be visualized:• A 2D colormap of the field in linear or logarithmic scale,• 3D iso-surfaces,• One given iso-contour or iso-surface,• The area delimited by the positive region of a levelset surface,• A zoom of a given area.
As an example, we represent one 9th order Lagrange shape function in the reference triangle.
In one case, we want to visualize the function using filled (colored) iso-values. In the othercase, we want to represent one given iso-value ue = C. Figure 4 show the two different adaptedvisualization grids relative to the two pre-defined visualization goals. We can do the same with

Figure 4. Visualization of one 9th order Lagrange shape function. On left, we see the visualization
mesh for the goal |u − uh | < 0.001. On center, we see the visualization mesh that has been optimizedin order to capture accurately one given iso-contour u = 0.5. On right, we see the corresponding
visualization results i.e a colormap and, in black, the iso-contour u = 0.5.

3D views on tetrahedral meshes and an example is shown at Figure 5.

5. Visualization error
We aim here to use classical h-refinement finite element techniques in order to optimize thevisualization mesh. For that, we should be able to estimate the visualization error. In a classicalposteriori error analysis, the “Quest of the Graal” is to find local error. Indeed, if the error
at a point is known, so is the exact solution. It is usually estimations of L2 or H 1 that areavailable. In this case, and contrary to the case of finite element a posteriori error estimation,
the exact error is known: the exact field, ue , being the one defined by the high order polynomialinterpolation and the approximate one, ueh being the piecewise linear field defined on thevisualization mesh. The exact local visualization error  is therefore simply defined as u e − ueh .
In this, our aim is to build an adapted visualization mesh while only using values computed
Copyright c 2005 John Wiley & Sons, Ltd.
Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2005; 00:1–6

8

J.F. REMACLE, N. CHEVAUGEON E. MARCHANDISE C. GEUZAINE

r=2

r=4

r=7

r = 7 plus vis. mesh

Figure 5. Visualization of one 4th order function on one single tetrahedron. Figures show one iso-surface
using different error thresholds. Last figure (bottom-right) show the visualization mesh together with
the iso-surface.

at the positions ξi , ηi , ζi defined by the interpolation matrix Iij . Let us consider the triangularelement depicted in Figure 6. We look if triangle e with nodal values u1 , u2 and u3 has tobe subdivided. At mid points, the exact solution u12 , u13 and u23 is available because itcorrespond to the ξi , ηi , ζi of child subdivisions. At those midpoints, linear approximations
would give average values 21 (u1 + u2 ), 12 (u1 + u3 ), and 21 (u2 + u3 ). We define the following error
indicator
e= max (|u1 + u2 − 2u12 | , |u1 + u3 − 2u13 | , |u2 + u3 − 2u23 |)= max (e12 , e13 , e23 )which is similar to what ca be found in [11]. Note that, for example,
e12 = |u1 + u2 − 2u12 | = h2η

∂2u+ O(h3η )∂η 2

so that the error indicator e is a measure of the maximal second derivative of the exact field.
We also see that, if we choose a maximal recursion level of r, then, element sizes are reduced
Copyright c 2005 John Wiley & Sons, Ltd.
Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2005; 00:1–6

9

EFFICIENT VISUALIZATION OF HIGH ORDER FINITE ELEMENTS

η
u2

u2
u223

u122
hη

u12

u1223
u23

u12

u23

u112
ξ
u1

u13

u3

u1233

u233

u1123
u1

u113 u13hξu133 u3hξFigure 6. One triangular element and its subdivision.

by a factor 2r and the visualization error is reduced by 22r . We proceed in the same way forother element types. Elements are subdivided recursively until e < ¯ where ¯ is a thresholdvalue.
The main advantage of this approach is its simplicity. Nodal values are computed once for
all and the error indicator only involves some few arithmetic operations. Hence, the approach
has a major drawback. If a feature has a characteristic size that is much smaller than the sizeof e, it can be missed by the procedure. This is especially true at low recursion levels. So, at
low recursion levels, more points are considered to compute the error indicator. The right partof Figure 6 shows the 2 recursion levels configuration. The error in element e is still computed
as the maximal numerical second derivatives at all interior points
e = max (e12 , e13 , e23 , e112 , e122 , e113 , e133 , e223 , e233 , e1123 , e1223 , e1233 )witheij = |2uij − ui − uj | ,13
eijj = 2uijj − ui − uj
2
2
and

1
eijjk = 2uijjk − uj − (ui + uj ) .
2
6. Visualization mesh adaptation

We now describe how we choose adaptively the sub-elements that are visible. The sub-element
data structure is described on Algorithm 2. We apply a recursive pattern for doing the meshrefinement and the subelement data structure is build like a multi dimensional tree. For
example, one triangular subelement has 4 children while a hexahedron subelement has 8.
Copyright c 2005 John Wiley & Sons, Ltd.
Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2005; 00:1–6

10

J.F. REMACLE, N. CHEVAUGEON E. MARCHANDISE C. GEUZAINE

There is only one instance of the subelement data structure : nodal values are updated foreach element using one BLAS2 matrix-vector multiplication. We call the root the initial
Algorithm 2 Visualization element class.
struct _subelement {staticint nbChilderen;
staticint nbNodes;
int
visible;int*node_numbers;_element *child;};unrefined element who is, in fact, the ancestor of all sub-elements. The mesh adaptation isdone in one recursion pass. On the way down of the recursion process, we compute the exact
error and decide weather it is useful to go to a downer level. On the way up of the recursion,
we may apply some filtering i.e. eliminate useless parts of the view. The filter may only allow
sub-elements whose nodal values are in a prescribed range (the case of the capture of one
given iso-surface). Another useful filter only allows to draw sub-elements that intersect a given
plane (the case of a planar cut in a 3D view). Algorithm 3 presents C++ code that is used todetermine which sub-elements of the root element are visible.

7. Two Dimensional Examples
7.1. Vortex-in-a-box
This example is usually chosen in order to test the ability of a numerical scheme to accurately
resolve thin filaments on the scale of the mesh which can occur in stretching and tearing flows.
We consider the following stream functionψ=1
sin2 (πx) sin2 (πy)π(2)

that defines the following velocity fieldvx = ∂ y ψ
vy = −∂x ψ=
cos(2πy) sin2 (πx)= − cos(2πx) sin2 (πy)(3)

We consider a disk of radius 0.15 placed at (0.5, 0.75) and the distance function u 20 =(x − 0.52 ) + (y − 0.75)2 − 0.152 . We aim to compute the following hyperbolic problem∂t u + v x ∂x u + v y ∂y u = 0
on a square domain of size [1, 1] and with u = u0 as initial conditions. The flow satisfies
vx = vy = 0 on the boundaries of the unit square. The resulting velocity field stretches out
the circle into a very long and thin filaments. For solving this problem, we use a 32 × 32
quadrilateral grid with a fourth-order (p = 4) DGM for space discretization. A fifth order
explicit Runge-Kutta is used for time stepping. Some visualization results are shown at time
step 1000.
Copyright c 2005 John Wiley & Sons, Ltd.
Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2005; 00:1–6

EFFICIENT VISUALIZATION OF HIGH ORDER FINITE ELEMENTS

11

Algorithm 3 Recursive algorithm that determines the set of visible sub-elements. The
algorithm is first called with the root element as first argument. It goes down in the subelement layers, stopping when the error is below a given threshold eps 0. If the algorithm
goes down to the deepest layer of the subelement tree, then the return value of the function
recur visible is false. More layers can be computed in order to ensure that the required
accuracy is reached.
bool recur_visible ( _subelement *e , double eps_0 ){if(!e->child[0]){// Max recursion has been attained without having// reached desired accuracy
e->visible = 1;
return false;}else// compute error using the 1- or the 2-level formula
double eps = compute_error ( e );if (eps > eps_0){
e->visible=0;
accuracy_reached = true;for (i=0;i<e->nbChilderen;++i)
accuracy_reached &= recur_visible (e->child[i],eps_0);
return accuracy_reached;}else{
e->visible=1;
return true;}}};

Figures 7 and 8 show an adaptive visualization. Figures 7–(a) and 7–(b) show the unrefinedmesh (1024 quads) and the unrefined visualization. Figures 7–(c) and 7–(d) show results fora visualization error of ¯ = 10−2 : the visualization mesh 7–(c) is made of 4021 quads andthe maximal recursion level required for obtaining the target error everywhere is r = 4. The
computation time for generating 7–(d) was 0.02 seconds. Figures 7–(e) and 7–(f) show results
for a visualization error of ¯ = 10−3 : the visualization mesh 7–(e) is made of 46747 quads andthe maximal recursion level required for obtaining the target error everywhere is r = 6. The
computation time for generating 7–(f) was 0.24 seconds.
Figure 8 shows an adaptive visualization that targets to capture accurately only the iso-zeroof function u. Figures 8–(a) and 8–(b) show the unrefined mesh (1024 quads) and the unrefinedvisualization. Figures 8–(c) and 8–(d) show results for a visualization error of ¯ = 10−2 : thevisualization mesh 8–(c) is made of 2032 quads. The computation time for generating 8–(d)was 0.03 seconds. Figures 8–(e) and 8–(f) show results for a visualization error of ¯ = 10−3 :Copyright c 2005 John Wiley & Sons, Ltd.
Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2005; 00:1–6

12

J.F. REMACLE, N. CHEVAUGEON E. MARCHANDISE C. GEUZAINE

the visualization mesh 8 (e) is made of 6454 quads and the computation time for generating
8–(f) was 0.08 seconds.
In the case ¯ = 10−3 , the maximal recursion level required for obtaining the desired accuracy
was r = 6. The equivalent uniformely refined mesh contains therefore 1024 × 4 6 quads, i.e. alittle more than four million quads. Our procedure allows to reduce this by a factor 650 in thecase of the iso-zero computation.
8. Three Dimensional Examples
8.1. Deformation of a sphereA sphere of radius .15 is placed within a unit computational domain at (.35, .35, .35) in avelocity field given by:vx= 2 sin2 (πx) sin(2πy) sin(2πz)g(t)vyvz= − sin(2πx) sin2 (πy) sin(2πz)g(t)= − sin(2πx) sin(2πy) sin2 (πz)g(t)(4)

The time dependence g(t) is given by:g(t) = cos(πt/T ),where the reversal time period is chosen to be T = 3. This three-dimensional incompressible
flow field combines a deformation in the x-y plane with a similar one in the x-z plane. The
sphere is subsequently deformed by the flow.
We have computed this problem using a non uniform hexahedral mesh composed of 32787
elements. Some results of visualization are presented in Figure 9. Left figures were computed
using a target error of ¯ = 10−3 while ¯ = 10−4 was used for the three right sub-figures of
Figure 9. The hehahedral refined mesh for the left and right figures is composed of 66172 and405784 elements respectively. Some 48426 triangles and 20585 quadrangles were used to draw
the iso-zero in the ¯ = 10−3 case. About 500000 triangles and 227000 quadrangles were usedin the ¯ = 10−4 case. Note that a maximal recursion level of r = 6 was necessary to obtain
the desired accuracy of ¯ = 10−4 . An element that has reached this level of recursion has been
cut into 86 = 262144 sub-elements. Only a fraction of these sub-elements are visible. A fullyrefined visualization mesh would have required to process 32787 × 262144 ' 8.5 10 9 elements.
Here, 405784 visualization cells correspond to about 21000 times less than the equivalent fullyrefined mesh. Figure 10 shows visualization results at different time steps for a target errorof ¯ = 10−4 . More than 4 million polygons are used at time step 2500. For that specific time
step, only 12 seconds of cpu time were necessary to build up the polygons, including theadaptive algorithm and the cutting of the 1084850 resulting hexahedra. Note that the visible
irregularities of the iso-surface at time step 2500 are caused by the numerical scheme (themesh is too coarse) and not to the accuracy of the visualization algorithm.
8.2. Propagation of acoustic modes in a duct.
We consider the problem of the propagation of acoustic modes in a quarter of a engine duct.
The geometry of the problem as well as the discretization mesh are represented in Figure
Copyright c 2005 John Wiley & Sons, Ltd.
Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2005; 00:1–6

EFFICIENT VISUALIZATION OF HIGH ORDER FINITE ELEMENTS

(a)

(b)

(c)

(d)

(e)

(f)

13

Figure 7. Visualization of function u at time step 1000 for the Vortex In a Box problem.

11. We have solved the Linearized Euler Equations (LEE) using 4th order polynomials anda Discontinuous Galerkin formulation. We have used our visualization algorithm for drawing
Copyright c 2005 John Wiley & Sons, Ltd.
Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2005; 00:1–6

14

J.F. REMACLE, N. CHEVAUGEON E. MARCHANDISE C. GEUZAINE

(a)

(b)

(c)

(d)

(e)

(f)

Figure 8. Visualization of the iso-zero of function u at time step 1000 for the Vortex In a Box problem.

contours of the acoustic pressure on both y = 0 and z = 0 planes with a target error of
Copyright c 2005 John Wiley & Sons, Ltd.
Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2005; 00:1–6

EFFICIENT VISUALIZATION OF HIGH ORDER FINITE ELEMENTS

non adaptive¯ = 10−3¯ = 10−4

non adaptive¯ = 10−3¯ = 10−4¯ = 10−315¯ = 10−4

Figure 9. Visualization of the iso-zero for the problem of the deformation of the sphere. Left figures
were computed using a target error of ¯ = 10−3 while ¯ = 10−4 was used for the three right sub-figures.
Top figures shows the iso-zero. Middle sub-figures show the visualization polygons. Bottom sub-figures
shows the volume visualization mesh.

¯ = 10−4 . Figure 13 shows visualization results in both refined and unrefined cases. Clearly,
the non-adaptive visualization strategy fails to obtain any relevant results. Figure 12 showsa zoom of the adaptive mesh in the z = 0 plane. Only 3 seconds of cpu time were needed
to generate the adaptive results. There are about one million polygons in each visualization
plane.
Copyright c 2005 John Wiley & Sons, Ltd.
Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2005; 00:1–6

16

J.F. REMACLE, N. CHEVAUGEON E. MARCHANDISE C. GEUZAINE

time step 0

time step 500

time step 1000

time step 1500

time step 2000

time step 2500

Figure 10. Adaptive visualization of the iso-zero for the problem of the deformation of the sphere atdifferent times for a target error of ¯ = 10−4 .

Copyright c 2005 John Wiley & Sons, Ltd.
Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2005; 00:1–6

EFFICIENT VISUALIZATION OF HIGH ORDER FINITE ELEMENTS

17

Figure 11. Geometry of the engine duct together with the computational mesh. The mesh is composedof 19395 tetrahedra.

Figure 12. View of the adaptive visualization plane in the z = 0 plane.

9. Conclusions
An adaptive technique for the visualization of high order finite element fields has beendeveloped. The technique is able to deal with general polynomial fields. Automatic mesh
Copyright c 2005 John Wiley & Sons, Ltd.
Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2005; 00:1–6

18

J.F. REMACLE, N. CHEVAUGEON E. MARCHANDISE C. GEUZAINE

Figure 13. Visualization of the acoustic pressure field on both y = 0 and z = 0 planes. Top sub-figure
shows the unrefined visualization results. Bottom sub-figure shows adaptive visualization results.

Copyright c 2005 John Wiley & Sons, Ltd.
Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2005; 00:1–6

EFFICIENT VISUALIZATION OF HIGH ORDER FINITE ELEMENTS

19

refinement (AMR) techniques have been used to generate optimal visualization grids. It has
been shown that the method was able to provide visualization results using only a small
fraction of the size that would have been required by an equivalent uniformely refined grid.
One of the main interests of the method developped here is its direct availability. An
implementation is provided in Gmshand this paper may be considered as a user’s guide.
As a future work, we will consider to extend the method the visualization of non polynomial
fields with a focus on the eXtended Finite Element Method (X-FEM). In X-FEM, nonpolynomial enrichments are used for the representation of the solution of linear elasticity
at the vivinity of the crack tip for example.

ACKNOWLEDGEMENTS

The various components of this work were supported by the European Community through Messiaen,
a Specific Target Research Programme of the 6th framework programme.

REFERENCES
1. M. Ainsworth. Dispersive and dissipative behavior of high order Discontinuous Galerkin finite element
methods. Journal of Computational Physics, 198(1):106–130, 2004.
2. Szabo B. and Babuska I. Finite Element Analysis. John Wiley and sons, 1991.
3. Ivo Babuska and Manil Suri. The p and h-p versions of the finite element method, basic principles andproperties. SIAM Rev., 36(4):578–632, 1994.
4. M.J. Berger and J. Oliger. Adaptive mesh refinement for hyperbolic partial differential equations. Journalof Computational Physics, 53:484–512, 1984.
5. B. Cockburn and C.-W. Shu. TVB Runge-Kutta local projection discontinuous Galerkin methods for
scalar conservation laws II: General framework. Mathematics of Computation, 52:411–435, 1989.
6. M. Deville, P. F. Fisher, and E. Mund. High-order methods for incompressible fluid flow. Cambridge
University Press, 2002.
7. J. J. Dongarra, J. Du Croz, I. S. Duff, and S. Hammarling. A set of Level 3 Basic Linear Algebra
Subprograms. ACM Trans. Math. Soft., 16:1–17, 1990.
8. C. Geuzaine and J.-F. Remacle. Gmsh web site. http://www.geuz.org/gmsh.
9. B. Haasdonk, M. Ohlberger, M. Rumpf, A. Schmidt, and K. G. Siebert. Multiresolution visualization ofhigher order adaptive finite element simulations. Computing, 70(3):181–204, 2003.
10. F. Hu and H. Atkins. Eigensolution analysis of the Discontinuous Galerkin Method with Nonuniform
Grids i One Space Dimension. Journal of Computational Physics, 182(2):516–545, 2002.
11. A. O. Leone, P. Marzano, E. Gobbetti, R. Scateni, and S. Pedinotti. Discontinuous finite element
visualization. In In the proceedings of th 8th International Symposium on Flow Visualisation, pages
1–6, 1998.
12. J.-F. Remacle, J.E. Flaherty, and M.S. Shephard. An adaptive discontinuous galerkin technique with an
orthogonal basis applied to compressible flow problems. SIAM Review, 45(1):53–72, 2003.
13. J.-F. Remacle, X. Li, J.E. Flaherty, and M.S. Shephard. Anisotropic adaptive simulation of transient
flows using discontinuous galerkin methods. International Journal for Numerical Methods in Engineering,62(7):899–923, 2005.
14. T.C. Warburton, S.J. Sherwin, and G.E. Karniadakis. Spectral basis functions for 2D hybrid hp elements.
SIAM Journal on Scientific Computing, 20:1671–1695, 1999.
15. Jue Yan and Chi-Wang Shu. A local discontinuous galerkin method for the kdv equations. SIAM Journalof Numerical Analysis, 40(2):769–791, 2002.

Copyright c 2005 John Wiley & Sons, Ltd.
Prepared using nmeauth.cls

Int. J. Numer. Meth. Engng 2005; 00:1–6