Computing cross fields - A PDE approach based on the Ginzburg-Landau theory

06/02/2017 ∙ by Pierre-Alexandre Beaufort, et al. ∙ Université catholique de Louvain 0

Cross fields are auxiliary in the generation of quadrangular meshes. A method to generate cross fields on surface manifolds is presented in this paper. Algebraic topology constraints on quadrangular meshes are first discussed. The duality between quadrangular meshes and cross fields is then outlined, and a generalization to cross fields of the Poincaré-Hopf theorem is proposed, which highlights some fundamental and important topological constraints on cross fields. A finite element formulation for the computation of cross fields is then presented, which is based on Ginzburg-Landau equations and makes use of edge-based Crouzeix-Raviart interpolation functions. It is first presented in the planar case, and then extended to a general surface manifold. Finally, application examples are solved and discussed.



There are no comments yet.


page 10

page 11

page 12

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

In the finite element community quadrangular elements are usually considered to be better than triangular elements, because there are twice as much triangles than quadrangles in a mesh for the same number of vertices. Moreover, quadrangular elements support tensorial operations, and ease the tracking of preferential directions in mesh refinement. Producing high quality quadrangular meshes is however not a trivial task. An appropriate approach (Lee and Lo, 1994) is the indirect one. It consists of (i) spawning points, (ii) triangulating them, and finally (iii) combining triangles into quadrangles.

Even if the two last steps are performed by efficient algorithms Remacle et al. (2012), the final mesh quality strongly depends on the adequacy of the position of the spawned points with a potential quadrangular grid. One strategy to place the points relies on a precomputed cross field. The ideal cross field should be smooth, topologically consistent, and aligned with boundaries and/or with geometrical features of the surface like curvature. Many methods have been developed to produce cross fields Bommes et al. (2009); Kälberer et al. (2007); Kowalski et al. (2013); Lai et al. (2010); Palacios and Zhang (2007); Ray et al. (2008); Cohen and Desbrun (2006); Vaxman et al. (2016). They either use an optimization process, or solve a PDE. All those papers provide deep insights into the modeling of direction fields.

In this paper, we propose to build cross fields based on Ginzburg-Landau theory. We provide then step by step the required mathematical background (§24). The partial differential cross field problem is formulated in variational form by the Ginzburg-Landau functional (§4.3). This functional (15) consists of a smoothness term that minimizes the gradient of the cross field and a penalty term that ensures its norm remains closed to unity. The asymptotic behavior (17) of Ginzburg-Landau functional naturally gives fields with well-distributed minimum critical points of minimum index. Crouzeix-Raviart finite elements are used for the interpolation (§5), easing the computation thanks to global/local representations of cross fields (Fig. 7). Finally, the nonlinear problem is solved using a Newton-Raphson scheme (19). Results on a unit sphere and oceans are provided in §6.

2 Topology of triangular and quadrilateral meshes

Assume an orientable surface embedded in . Let be the number of handles of the surface. The topological characteristic , which is also called the genus of the surface, is the maximum number of cuttings along non-intersecting closed curves that won’t make the surface disconnected. Let also be the number of connected components of the boundary of the surface. The Euler characteristic of is then the integer

One has for a sphere, whereas for a disk (), and for a torus () or a cylinder ().

Consider now a mesh on with nodes (also called vertices), edges and facets. The Euler formula


provides a general relationship betweeen the numbers of nodes, edges and facets in the mesh Eppstein (2009). If nodes (and hence edges) are on the boundary , and if the number of edges (or nodes) per facet is noted ( for triangulations and for quadrangulations, meshes mixing triangles with quadrangles being excluded), the following identity holds : all facets have edges, edges have two adjacent facets and edges have one adjacent facet. Hence the relationship


Elimination of between (2) and (1) yields


which is true for any triangulation or quadrangulation.

A regular mesh has only regular vertices. An internal vertex is regular if it has exactly adjacent triangles or adjacent quadrangles, whereas a boundary vertex is regular if it has exactly adjacent triangles or adjacent quadrangles. One has then




respectively for a regular triangulation and a regular quadrangulation. Substitution of (4) and (5) into (3) shows that only surfaces with a zero Euler characteristic can be paved with a regular mesh. If , irregular vertices will necessarily be present in the mesh.

The number and the index of the irregular vertices is tightly linked to the Euler characteristic , which is a topological invariant of the surface. We call valence of a vertex the number of facets adjacent to the vertex in the mesh. In a regular mesh, all vertices have the same valence . In a non regular mesh, on the other hand, a number of irregular vertices have a valence , and one notes the integer the valence mismatch of a vertex.

Assume a quadrangulation with irregular internal vertices of valence , and irregular boundary vertices of valence , given. All other vertices are regular. There are then regular internal vertices of valence , and regular boundary vertices of valence , so that one can write


and the substraction of (3) with yields

showing that, in a quadrangulation, each irregular vertex counts for in the Euler characteristic, a quantity called the indice of the irregular vertex .

Summing up now on different possible values for , one can establish that a quadrangulation of a surface with Euler characteristic verifies


Consider, for instance, the quadrangulation of a disk, which is a surface with . A minimum of irregular vertices of index must be present. They can be located either on the boundary (vertices of valence 1) or inside the disk (vertices of valence 3), Fig. 1.

Figure 1: A quadrilateral mesh of a circle. Four irregular vertices of index (in red) are required to obtain such a mesh. The irregular vertices may be inside the disk (left) or on its boundary (right)
Figure 2: Different quadrangulations of a L-shaped domain. Irregular vertices of index are displayed in red, whereas ones of index are displayed in blue. The sum of the indices of the irregular vertices is equal to in all cases.

Fig. 2 shows three different quadangulations of a L-shaped domain (). Regular boundary nodes should all have a valence of 2. The mesh on the left has irregular vertices located at the corners of the domain : five with index , and one with index . The central mesh, on the other hand, has the minimum amount of irregular vertices, i.e., four ones of index . The right mesh generated by recombination of a standard Delaunay triangular mesh Remacle et al. (2012) has eight vertices of index , and twelve vertices of index , both on the boundary and inside the domain. Quality meshes should have as few irregular vertices as possible. In what follows, a general method allowing to position irregular vertices before meshing the surface is presented.

3 Why cross fields?

Cross fields are auxiliary in the generation of quadrangular meshes. We shall show that nonregular vertices defined in the previous section are precisely the critical points of a cross field, and that these critical points of the cross field can also be related to the Euler characteristic of the meshed surface. This result represents an important theoretical limit on the regularity of quadrangular meshes.

3.1 Continuity

A cross field is a field defined on a surface with values in the quotient space , where is the circle group and is the group of quadrilateral symmetry. Pictorially, it associates to each point of the surface

to be meshed a cross made of four unit vectors orthogonal with each others in the tangent plane

of the surface.

Figure 3: Diffential function .

A surface can be identified with its tangent space in any neighborhood that is sufficiently small to have curvature effects negligible. This local identification of the surface with a vector space endows it with a natural parallel transport rule, so that the angular differential can be defined as the minimal angle, with its sign, between the branches of and any of the branches of for any pair of points where is defined, Fig. 3. Taking now as reference the cross , an angular coordinate


can be defined for crosses in . The cross field is deemed continuous (regular) in if the limit


exists (i.e., is unique). It is then equal to . Isolated points , of where the limit (9) does not exist are called critical points or zeros of the cross field.

3.2 Index and degree

Although defined locally, the notion of continuity gives unexpectedly valuable information about the topology of , which is a nonlocal concept. To see this, consider a cross field defined on a quadrangular element delimited by four (possibly curvilinear) edges. Assume the cross field is parallel to the four edges (i.e., one of the four branches of the cross is parallel to the tangent vector of the edge at each point of the edge, except the extremities) and prolongates smoothly inside the quadrangle. This field is discontinuous at corners where edges do not meet at right angle, but it is continous everywhere else. Making the same construction for all elements of a quadrangular mesh, one obtains a cross field topologically identified with the quadrangular mesh, and that is continuous everywhere except at the vertices of the mesh. This field has thus got isolated critical points at mesh vertices, but not all critical points have the same significance. Some critical points have a specific topological value, associated with the notion of index.

To introduce the notion of index, an angular coordinate needs to be defined for points in a neighborhood of a critical point . Picking up an arbitrary regular point , , the local unit vector basis

with the normal to , is constructed, and hence a local polar coordinate system


can be defined for points in .

A circular curve of infinitesimal radius centered around the vertex is now considered. As the angles (10) and (8) are precisely the elements of the groups and , respectively, the cross field on can be regarded as a mapping


The mapping is continuous, since circles around the critical point , but does not cross it. The index of at is the degree of the mapping (11), i.e., the number of times the codomain wraps around the domain under the mapping. Its algebraic expression is easily expressed in terms of the angles and as

where is . In case of a vertex of valence , i.e., a vertex adjacent to quadrangular elements, the integral evaluates as


where the ’s are the angles of the quadrangular elements adjacent to the considered vertex , and where the obvious relationship has been used. The cross field has index 0 at vertices adjacent to four quadrangular elements, whereas it has index , resp. at vertices adjacent to 3, resp. 5, quadrangular elements meet, Fig. 4. As one sees, the index is a topological characteristic of the cross field at the critical point . It does not depend on the choice of the curve , nor on the choice of an angular reference for the angles and .

Figure 4: Illustration of vertices where the indices of the cross field represented in red are respectively , from left to right. The index only depends on the number of quadrangles adjacent to the vertex, independently of the values of the angles , which need not be identical as they are in the figure.

3.3 Poincaré-Hopf theorem

Equation (12), relates the index of the cross field at a critical point with one forth of valence of the corresponding mesh vertex. This result can be combined with the algebraic topology result of previous section (7) that each internal irregular vertex of valence counts for in the Euler characteristic of the underlying surface. This yields the relationship


for the critical points of a cross field defined on a surface .

This is a generalization Poincaré-Hopf theorem, which states that the sum of the indices of the critical points of a vector field defined on a surface without boundary is equal to the Euler characteristic of the surface. This famous theorem draws an unexpected and profound link between two apparently distinct areas of mathematics, topology and analysis. Whereas vector fields have integer indices at critical points, cross fields have indices that are multiples of 1/4. Still the topological relationship (13) of Poincaré-Hopf holds in both cases. Actually, our developments reach same inferences than Ray et al. (2008).

4 Cross field computation : the planar case

4.1 Vector representation of cross fields

Only scalar quantities can be compared at different points of a manifold. For the comparison or, more generally, for differential calculus with nonscalar quantities like cross fields, a parallel transport rule needs to be defined on the manifold. On a surface manifold, this rule can take the form of a regular vector field which gives at each point the direction of the reference angle 0. Poincaré-Hopf theorem says that such a field does not exist in general, and in particular on manifolds whose Euler characteristic is not zero. The situation is however easier in the planar case. A global Cartesian coordinate frame can always be defined over the plane, and be used to evaluate the orientation of the cross field. We shall therefore expose the cross field computation method in the planar case first, and generalize to nonplanar surfaces, where we will have to deal with local reference frames, in a subsequent section.

A cross is an element of the group , which can be represented by the angle it forms with the local reference frame. Yet, due to the quadrilateral symmetry, four different angles in represent the same cross field . Let for instance the angles and represent the same cross. The average represents another cross, whereas the difference is not zero. So, we have and , which clearly indicates that the values of the cross field do not live in a linear (affine) space. This makes the representation by improper for finite element interpolation. The solution is two-fold. First, the angle is multiplied by four, so that the group be mapped on the unit circle , and the cross be therefore represented by a unit norm vector . Then, the vector is represented in components in the reference frame as

4.2 Laplacian smoothing

Computing the cross field consists thus now of computing the vector field representation , which obviously lives in a linear space (a 2D plane). The components of are fixed on the boundaries of so that the crosses are parallel with the exterior normal vector i.e.,

Propagating inside is here be done by solving a Laplacian problem. Even though the vector representation is unitary on , it tends to drift away from inside the domain. The computed finite element solution lies therefore outside the unit circle and must be projected back on to recover the angle

Due to the multiplication by 4, the indices of the critical points of the vector field verify


4.3 The Ginzburg-Landau model

Numerical experiments show that the norm of the vector field computed by Laplacian smoothing (See previous section) decreases quite rapidly as one moves away from the boundary , leaving in practice large zones in the bulk of the computation domain where the solution is small, and the computed cross field inaccurate. A more satisfactory formulation consists of ensuring that the norm of remains unitary over the whole computation domain. This problem can be formulated in variational form in terms of the Ginzburg-Landau functional


The first term minimizes the gradient of the cross field and is therefore responsible for the laplacian smoothing introduced in the previous section. The second term is a penality term that vanishes when . The penality parameter , called coherence length, has the dimension of a length. The Euler-Lagrange equations of the functional (15) are the quasi-linear PDE’s


called Ginzburg-Landau equations. If is small (enough) with respect to the dimension of , then is of norm everywhere but in the vicinity of the isolated critical points .

The asymptotic behavior of Ginzburg-Landau energy can be written as




as (see Bethuel et al. (2012), Introduction, Formulae 11 and 12).

In asymptotic regime, the energy is thus composed of three terms. The first term of (17) blows up as , i.e., energy becomes unbounded if critical points are present. When is small, this first term dominates, and one is essentially minimizing with the constraint (14). This indicates that a critical point of index has a cost of in terms of energy, whereas 2 critical points of index have a cost of . All critical points should therefore be of index , and their number should be . This is indeed good news for our purpose : good cross fields should have few critical points of lower indices.

The second term of (17) is the renormalized energy (18). It remains bounded when tends to . The double sum in reveals the existence of a logarithmic force between critical points. The force is attractive between critical points with indices of opposite signs, and a repulsive between critical points with indices of the same signs. The second term in (18) is more complicated and is detailed in Bethuel et al. (2012). Basically, represents a repulsing force that forbids critical points to approach the boundaries.

Finally, the third term in (17) vanishes as . At the limit, all energy is thus carried by the critical points of the field. All this together allows to believe that Ginzburg-Landau model is a good choice for computing cross fields. It produces few critical points, which are moreover well-distributed over the domain.

5 Computation of cross fields: nonplanar generalization

The finite element computation method for cross fields is now generalized to the case of nonplanar surfaces. Consider the conformal triangulation of a nonplanar surface manifold , each triangle being defined by the vertices , and . Since no global reference frame exists on a nonplanar surface, a local reference frame is associated to each edge of the triangulation. Let be the edge of the mesh, joining nodes and , and be the average of the normals vectors of the two triangles adjacent to . The vectors

form a local frame in which to represent the connector values of the discretized cross field ,

which are attached to the center of the edges of the triangulation. Actually, is assumed to be the same along within both planes of triangles sharing . This assumption eases computation and gives a planar-like representation, Fig. (a)a.

(a) Global representation.

(b) Local representation.
Figure 7: Cross field over the edge of a mesh.

As the connector values are attached to the edges of the mesh, and not to the nodes, Crouzeix-Raviart interpolation functions are used instead of conventional Lagrange shape functions Crouzeix and Raviart (1973). The Crouzeix-Raviart shape functions equal on corresponding edge , and on the opposite vertices (Fig. 8) in the two adjacent triangular elements. They are polynomial and their analytic expression in the reference triangle reads

where indices and enclosed in parenthesis denote the local edge numbering in the considered triangular element.

Figure 8: Third Crouzeix-Raviart function shape (shaded in grey) over reference triangle (in blue).

Each of the three edges of a triangle has its own local reference frame. If one is to interpolate expressions involving the vector field over this element, the three edge-based reference frames have to be appropriately related with each other Ray et al. (2016). We arbitrarily take the reference frame of the first edge of the element as reference, and express the angular coordinate of the two other edges in function of this one with the relationships (Fig. (b)b)

Thus, the 6 local unknowns of triangle can be expressed as a function of the 6 edge unknowns by

and we have the interpolation

for the vector field in the triangle .

A Newton scheme is proposed to converge to the solution. The Newton iteration at stage for solving (16) consists of solving:


The elementary matrix and the elementary vector of element are then given by




It is then necessary to transform those elementary matrix and vector in the reference frames of the edges as

Then, standard finite element assembly can be performed. Boundary conditions are simply

on every edge of . This nice simplification is due to the fact that unknowns are defined on the reference frame of the edges.

6 Results

6.1 The sphere

As a first example, let us compute the cross field on a unit sphere. The sphere has no boundary so we choose randomly one edge of the mesh and fix the cross field for this specific edge. The mesh of the sphere is made of 2960 triangles (see Fig. 9). A value of was chosen for the computation, which corresponds roughly to twice of the mesh size. A total of 29 Newton iterations were necessary to converge, by reducing the residual norm to .

Figure 9: Mesh of the sphere. Colors correspond to the 2-norm of the cross field. The 8 critical points are located on two squares of side which correspond to the size of the inscribed cube. The two squares are tilted by 45 degrees.

The location of the critical points is indeed not what we expected: our initial intuition was that critical points would be located at the corners of an inscribed cube of side . In all our computations i.e. while changing the mesh and , critical points are located on two squares of side , those two squares being tilded by degrees around their common axe (see Fig. 9). Equilateral triangle patterns are formed between critical points that belong to both squares. In reality, our solution is the right solution. In the asymptotic regime, the location of the critical points tends to minimize (see Equations (17) and (18)). We have thus computed for tilting angles ranging from to . Fig. 10 shows clearly that the minimum of the energy corresponds to an angle of , which is exactly what is found by the finite element formulation.

Figure 10: Energy vs. tilting angle for the sphere. The minimum corresponds to a titing angle of .

Fig. 11 shows the cross field as well as the separatrices. The separatrices were computed “by hand”.

Figure 11: Separatrices from cross fields.

The solution that has been found is related to what is called the Whyte’s problem Saff and Kuijlaars (1997); Dragnev et al. (2002) that consists of finding points on the sphere which positions maximizes the product of their distances. The critical points are called logarithmic extreme points or elliptic Fekete points Fekete (1923).

The specific configuration that corresponds to is called an anticube (or square antiprism) and is exactly the one that was found numerically.

At this point, one could wonder what happens if the symmetry group of the hexagon was used instead of the symmetry group of the square in the definition of the cross field. Pictorially, the cross field is then an “asterisk” field. In this case, we expect to find critical points distributed in such a way that Ginzburg-Landau energy is minimized. Changing the group of symmetry simply amounts in our code in changing the rotation matrix . The results for the spherical surface are depicted in Fig. 12. Critical points are the summits of an icosahedron, which is the solution of Whyte’s problem for . This superb result shows that it is indeed possible to use cross fields not only for building quadrangles but also to build equilateral triangles.

Figure 12: Asterisk field (6 turns) with their critical points on the corners of an icosahedron.

6.2 World oceans and seas

As it has been said previously, direction fields are used to generate points which are adequate for quadrangles (from cross fields), either equilateral triangles (from asterisk fields). Here, we show meshes that are built from such direction fields, using algorithms of Georgiadis et al. (2017). Some oceans and seas (up to geometry resolution and mesh size) are meshed over the Earth with only quadrangles, Fig. 13.

Figure 13: Quadrangular mesh of oceans and seas, built from cross fields.

On Fig. 16, we give results of cross (Fig. (a)a) and asterisk (Fig. (b)b) fields. Observe that it is a zoom onto center of Fig. 13. Branches which color tends to blue (and thus zero norm) fairly correspond to critical points. We notice that critical points lay near boundaries (i.e. coasts). It means that repulsion of critical points is stronger than repulsion of boundary in (18). In other words, critical points tend to be close to boundaries (if any). There are more critical points within asterisk fields, than within cross fields. It is due to the fact that a nonregular vertex in a triangular mesh has an index of , while in a quadrangular mesh it is . Finally, we see that direction fields in both cases are smooth and enable production of high quality meshes.

(a) From cross fields.
(b) From asterisk fields.
Figure 16: Output meshes from direction fields. Colorbar represents norm of direction branches.

7 Conclusion

This article has mainly underlined the consistency of using Ginzburg-Landau equations to compute cross fields for indirect quadrangles meshing, and even so-called asterisk fields in order to produce equilateral triangles. The interest of using Ginzburg-Landau functional is related to its asymptotic behavior as the coherence length tends to zero: it gives optimal critical points of direction fields, i.e. those of lower index as it is expected by continuous and discrete topology. Besides, Ginzburg-Landau functional is intuitive with its two complementary terms: smoothing and penalty; it draws direction fields with aimed characteristics, as it has be shown by our results.


The present study was carried out in the framework of the project ”Large Scale Simulation of Waves in Complex Media”, which is funded by the Communauté Française de Belgique under contract ARC WAVES 15/19-03 with the aim of developing and using Gmsh Geuzaine and Remacle (2009).


  • Lee and Lo (1994) C. K. Lee, S. Lo, A new scheme for the generation of a graded quadrilateral mesh, Computers & structures 52 (1994) 847–857.
  • Remacle et al. (2012) J.-F. Remacle, J. Lambrechts, B. Seny, E. Marchandise, A. Johnen, C. Geuzaine, Blossom-Quad: A non-uniform quadrilateral mesh generator using a minimum-cost perfect-matching algorithm, International Journal for Numerical Methods in Engineering 89 (2012) 1102–1119.
  • Bommes et al. (2009) D. Bommes, H. Zimmer, L. Kobbelt, H. Zimmer, L. Kobbelt, Mixed-integer quadrangulation, ACM Transactions on Graphics 28 (2009) 1.
  • Kälberer et al. (2007) F. Kälberer, M. Nieser, K. Polthier, Quadcover-surface parameterization using branched coverings, in: Computer graphics forum, volume 26, Wiley Online Library, 2007, pp. 375–384.
  • Kowalski et al. (2013) N. Kowalski, F. Ledoux, P. Frey, A pde based approach to multidomain partitioning and quadrilateral meshing, Proceedings of the 21st international meshing roundtable (2013) 137–154.
  • Lai et al. (2010) Y.-K. Lai, M. Jin, X. Xie, Y. He, J. Palacios, E. Zhang, S.-M. Hu, X. Gu, Metric-driven rosy field design and remeshing, IEEE Transactions on Visualization and Computer Graphics 16 (2010) 95–108.
  • Palacios and Zhang (2007) J. Palacios, E. Zhang, Rotational symmetry field design on surfaces, ACM Transactions on Graphics (TOG) 26 (2007) 55.
  • Ray et al. (2008) N. Ray, B. Vallet, W. C. Li, B. Lévy, N-symmetry direction field design, ACM Transactions on Graphics 27 (2008) 1–13.
  • Cohen and Desbrun (2006) Y. T. P. A. D. Cohen, S. M. Desbrun, Designing quadrangulations with discrete harmonic forms, in: Eurographics symposium on geometry processing, 2006, pp. 1–10.
  • Vaxman et al. (2016) A. Vaxman, M. Campen, O. Diamanti, D. Panozzo, D. Bommes, K. Hildebrandt, M. Ben-Chen, Directional field synthesis, design, and processing, in: Computer Graphics Forum, volume 35, Wiley Online Library, 2016, pp. 545–572.
  • Eppstein (2009) D. Eppstein, Nineteen proofs of Euler’s formula: , Information and Computer Sciences, University of California, Irvine (2009).
  • Bethuel et al. (2012) F. Bethuel, H. Brezis, F. Hélein, Ginzburg-Landau Vortices, volume 13, Springer Science & Business Media, 2012.
  • Crouzeix and Raviart (1973) M. Crouzeix, P.-A. Raviart, Conforming and nonconforming finite element methods for solving the stationary stokes equations i, Revue française d’automatique, informatique, recherche opérationnelle. Mathématique 7 (1973) 33–75.
  • Ray et al. (2016) N. Ray, D. Sokolov, B. Lévy, Practical 3d frame field generation, ACM Transactions on Graphics (TOG) 35 (2016) 233.
  • Saff and Kuijlaars (1997) E. B. Saff, A. B. Kuijlaars, Distributing many points on a sphere, The mathematical intelligencer 19 (1997) 5–11.
  • Dragnev et al. (2002) P. D. Dragnev, D. Legg, D. Townsend, On the separation of logarithmic points on the sphere, in: Approximation Theory X: Abstract and Classical Analysis, Vanderbilt University Press, Nashville, TN, 2002, pp. 137–144.
  • Fekete (1923) M. Fekete, Über die verteilung der wurzeln bei gewissen algebraischen gleichungen mit ganzzahligen koeffizienten, Mathematische Zeitschrift 17 (1923) 228–249.
  • Georgiadis et al. (2017) C. Georgiadis, P.-A. Beaufort, J. Lambrechts, J.-F. Remacle, Mesh generation on the sphere using cross and asterisk fields (2017). Submitted as a research note for the IMR26.
  • Geuzaine and Remacle (2009) C. Geuzaine, J.-F. Remacle, Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities, International Journal for Numerical Methods in Engineering 79 (2009) 1309–1331.