A fictitious domain finite element method for simulations of fluidstructure interactions: The NavierStokes equations coupled with a moving solid
Abstract
The paper extends a stabilized fictitious domain finite element method initially developed for the Stokes problem to the incompressible NavierStokes equations coupled with a moving solid. This method presents the advantage to predict an optimal approximation of the normal stress tensor at the interface. The dynamics of the solid is governed by the Newton’s laws and the interface between the fluid and the structure is materialized by a levelset which cuts the elements of the mesh. An algorithm is proposed in order to treat the time evolution of the geometry and numerical results are presented on a classical benchmark of the motion of a disk falling in a channel.
keywords:
Fluidstructure interactions, NavierStokes, Fictitious domain, eXtended Finite Element.S. COURT, M. FOURNIÉ. FICTITIOUS DOMAIN FOR FLUID STRUCTURE INTERACTION 0
1 Introduction
Fluidstructure interactions problems remain a challenge both for a comprehensive study of such problems as for the development of robust numerical methods (see a review in (1)). One class of numerical methods is based on meshes that are conformed to the interface where the physical boundary conditions are imposed (2); (3); (4). As the geometry of the fluid domain changes through the time, remeshing is needed, which is excessively timeconsuming, in particular for complex systems. An other class of numerical methods is based on nonconforming mesh with a fictitious domain approach where the mesh is cut by the boundary. Most of the nonconforming mesh methods are based on the immersed boundary methods where forceequivalent terms are added to the fluid equations in order to represent the fluid structure interaction (5); (6). Many related numerical methods have been developed, in particular the popular distributed Lagrange multiplier method, introduced for rigid bodies moving in an incompressible flow (7). In this method, the fluid domain is extended in order to cover the rigid domain where the fluid velocity is required to be equal to the rigid body velocity.
More recently, eXtended Finite Element Method introduced by Moës, Dolbow and Belytschko in (8) (see a review of such methods in (9)) has been adapted to fluid structure interactions problems in (10); (11); (12); (13). The idea is similar to the fictitious domain / Lagrange multiplier method aforementioned, but the fluid velocity is no longer extended inside the structure domain, and its value given by the structure velocity is enforced by a Lagrange multiplier only on the fluidstructure interface. One thus gets rid of unnecessary fluid unknowns. Besides, one easily recovers the normal trace of the Cauchy stress tensor on the interface. We note that this method has been originally developed for problems in structural mechanics mostly in the context of cracked domains, see for example (14); (15); (16); (17); (18). The specificity of the method is that it combines a levelset representation of the geometry of the crack with an enrichment of a finite element space by singular and discontinuous functions.
In the context of fluidstructure interactions, the difficulty related to the applications of such techniques lies in the choice of the Lagrange multiplier space used in order to take into account the interface, which is not trivial because of the fact that the interface cuts the mesh (see (19) for instance).
In particular, the natural mesh given by the points of intersection of the interface with the global mesh cannot be used directly. An algorithm to construct a multiplier space satisfying the infsup condition is developed in (19), but its implementation can be difficult in practice.
The method proposed in the present paper tackles this difficulty by using a stabilization technique proposed in (14). This method was adapted to contact problems in elastostatics in (20) and more recently to the Stokes problem in (21).
An important feature of this method (based on the eXtended Finite Element Method approach, similarly to (12); (13)) is that the Lagrange multiplier is identified with the normal trace of the Cauchy stress tensor at the interface.
Moreover, it is possible to obtain a good numerical approximation of (the proof is given in (21) for the Stokes problem). This property is crucial in fluidstructure interactions since this quantity gives the force exerted by the viscous fluid on the structure. In the present paper, we propose to extend this method to the NavierStokes equations coupled with a moving solid. Note that alternative methods based on the Nitsche’s work (22) (such as (23); (24) in the context of the Poisson problem and (25) in the context of the Stokes problems) do not introduce the Lagrange multiplier and thus do not necessarily provide a good numerical approximation of this force. Our method based on boundary forces is particular interesting for control flow around a structure. The control function can be localized on the boundary of the structure where we impose its local deformation. In order to perform direct numerical simulations of such a control, efficient tools based on accurate computations on the interface must be developed. The present approach is one brick in this research topic where recent development towards stabilized NavierStokes equations are proposed (like in (26)).
The outline of the paper is as follows. The continuous fluidstructure interactions problem is given in Section 2 and the weak formulation with the introduction of a Lagrange multiplier for imposing the boundary condition at the interface is given in Section 2.2. Next, in Section 3 the fictitious domain method is recalled with the introduction of the finite element method (Section 3.1) with a time discretization (Section 3.2). Section 4 is devoted to numerical tests and validation on a benchmark corresponding to the falling of a disk in a channel. The efficiency of the method is presented before conclusion.
2 The model
2.1 Fluidstructure interactions
We consider a moving solid which occupies a timedepending domain denoted by . The remaining domain corresponds to the fluid flow.
The displacement of a rigid solid can be given by the knowledge of , namely the position of its gravity center, and its rotation given by for , , where is the rotation angle of the solid (see Figure 1). Then at time the domain occupied by the structure is given by
Remark: This formulation can be extended in order to consider general deformations of the structure. Then we would have to define a mapping which corresponds to the deformation of the solid in its own frame of reference. Then, where , for .
The velocity of the incompressible viscous fluid of density is denoted by , the pressure by and is the dynamic viscosity. We denote by the outward unit normal vector to (the boundary of ), and the normal trace on the interface of the Cauchy stress tensor is given by
When gravity forces are considered (we denote by the gravity field), the fluid flow is modeled by the incompressible NavierStokes equations
(1) 
and the Newton’s laws are considered for the dynamics of the solid
(2) 
where is the mass of the solid, and is its moment of inertia.
At the interface , for the coupling between fluid and structure, we impose the continuity of the velocity
(3) 
The coupled system (1)–(3) has for unknowns , , and the angular velocity (a scalar function in 2D).
2.2 Weak formulation of the problem with stabilization terms
We consider the coupled system (1)(3) and we assume that the boundary condition imposed at the interface is sufficiently regular to make sense, and we introduce the following functional spaces (based on the classical Sobolev spaces , , and , see (27) for instance)
Due to the fact that we only consider boundary conditions of Dirichlet type, we impose to the pressure to have null average (this condition is taken into account in ). That variational formulation can be done in three steps:

Classical formulation of the NavierStokes and structure equations (in the formulation and are equal to );

Introduction of Lagrange multiplier in order to take into account the Dirichlet condition at the interface (in the formulation only is equal to );

Introduction of stabilization terms with a parameter .
This new unknown plays a critical role due to the fact that it is equal to the normal trace of the Cauchy stress tensor (this equality is described in (28)). The stabilization terms are associated with the constant parameter (chosen sufficiently small). The variational problem that we consider is the following:
Find such that  
where
Remark: The formulation can be justified by the introduction of an extended Lagrangian  à la BarbosaHughes, see (29)  whose a stationary point is a weak solution of the problem. The firstorder derivatives of this Lagrangian leads to forcing to reach the desired value corresponding to .
3 Fictitious domain approach
We refer to the article (21) for the details of the fictitious domain approach we consider here. In the following, we recall the method used for the present work.
3.1 Finite element discretization
The fictitious domain for the fluid is considered on the whole domain . Let us introduce three discrete finite element spaces, , and . Notice that the spaces introduced to define the weak formulation are included into those spaces defined all over the domain . In practice, is a simple domain, so that the construction of a unique mesh for all spaces is straightforward (the interface between the fluid and the structure is not considered). Let us consider for instance a rectangular domain where a structured uniform mesh can be constructed (see Figure 3). Classical finite element discretizations can be defined on the spaces , and . For , let us consider for instance a subspace of the continuous functions defined by
where is a finite dimensional space of regular functions, containing the polynom space of degree less or equal to an integer (). For more details, see (30) for instance. The mesh step stands for , where is the diameter of . In order to split the fluid domain and the structure domain, we define spaces on the fluid part and on the interface only, as
Notice that , , are respective natural discretizations of , and . It corresponds to cutting the basis functions of spaces , and , as shown in Figure 2. This approach is equivalent to the eXtended Finite Element Method, as proposed in (13) or (12), where the standard finite element method basis functions are multiplied by Heaviside functions ( for and for ), and the products are substituted in the variational formulation of the problem. Thus the degrees of freedom inside the fluid domain are used in the same way as in the standard finite element method, whereas the degrees of freedom in the solid domain at the vertexes of the elements cut by the interface (the so called virtual degrees of freedom) do not define the field variable at these nodes, but they are necessary to define the fields on and to compute the integrals over . The remaining degrees of freedom, corresponding to the basis functions with support completely outside of the fluid, are eliminated (see Figure 3). We refer to the papers aforementioned for more details.
The discrete problem consists in finding such that
This is a system of nonlinear differential algebraic equations which can be formulated into a compact form. We denote by , and the respective degrees of freedom of , and . After standard finite element discretization of the following bilinear forms
we define matrices like from , etc…, the vector from , from the gravity forces , the matrix computed by integration over of the basis functions and the matrix depending on the velocity and corresponding to the nonlinear convective term . Then the matrix formulation is given by
(4)  
(5)  
(6)  
(7)  
(8) 
At the interface represented by a levelset function which cuts the global mesh, the coupling between the fluid and the structure is imposed by a Dirichlet condition whose elements are determined through the computation of . The main advantage of our numerical method  mathematically justified in (21)  is to return an optimal approximation of at the interface. Getting a good approximation of this quantity is crucial for the dynamics of the system.
3.2 Time discretization and treatment of the nonlinearity
Classical methods like methods can be used for the time discretization. For a matter of unconditional stability of the scheme, we consider an implicit discretization based on the backward Euler method. We denote by the solution at the time level and is the time step. Particular attention must be done for a moving particle problem. Indeed, at the time level the solid occupies which is different from the previous time level . So, the field variable at the time level can become undefined near the interface since there was no fluid flow at the time level ( for the solid and for the fluid). In other words, some degrees of freedom for the fluid part which are not considered at the time level must be taken into account at the time level . In particular, the velocity field must be known in such nodes. In the present work, we impose the velocity to be equal to the motion of the solid. The validity of this approximation is justified as soon as time step is sufficiently small to ensure that the structure moves progressively across the mesh without jump of cells (when the levelset doesn’t cuts this cell). This constraint is not too strong and corresponds to the classical CFL condition for velocity of the structure.
In the following we present the algorithm we perform to compute at the time level the solution () on . To simplify, we assume that is constant. At the time level we have access to () on .

We update the geometry to determine . It corresponds to update the position of the levelset which is defined from and .

We complete the velocity defined on to the full domain by imposing the velocity on each node of to be equal to .
After this step, we know the Dirichlet condition for the velocity to impose at the interface . So we determine in (6) from . 
Finally, we compute such that
At this stage, the solution of the resulting nonlinear algebraic system is achieved by a Newton method. The initialization of the Newton algorithm is done with the solution at the previous time step (this solution is defined in item 4–).

We complete the velocity defined on to the full domain by imposing the velocity on each node of to be equal to .
Remark 1: In practice the midpoint method is used to update the geometry of the structure for the computation of .
Remark 2: Step 6– plays an important role to update the geometry. Indeed, after extension, we have access to the values of the solution at each node of the full domain. Thus new nodes that appear after update have some values and no interpolation is required.
4 Numerical tests: Free fall of a disk in a channel
For validation of our method, we consider the numerical simulation of the motion of a disk falling inside an incompressible Newtonian viscous fluid. The parameters used in the computation, for a disk of radius cm in a channel of dimension , are given in Table 1.
Parameter  

Unit  g/cm  g/cm  g/cm s  cm/s 
value  1  1.25  0.1  981 
This simulation is well documented in the literature and considered as a challenging benchmark.
We refer to the paper (7) where fictitious domain method is used and (31) for simulations with mesh adaptation.
For the finite element discretization, we consider classical Lagrange family with for respectively , , and , which is a choice that satisfies
the infsup condition required for such kind of problems (see (21) for more details). Uniform triangular meshes are used and defined by imposing a uniform repartition of points on the boundary of the domain. Two meshes are used, with points in direction and points in direction, and with points in direction and points in direction.
Numerical tests are performed with and without stabilization to underline the advantage of the method. When stabilization is considered, we choose where (see (21) for the justification of this choice). The parameter has to obey to a compromise between the coerciveness of the system and the weight of the stabilization term. The time discretization step is initialized to and adapted at each time iteration to satisfy a CFL condition. More precisely, we evaluate the norm of the velocity at each point of the structure and we deduce the maximum value . Then we impose . This condition is not restrictive and we observe that in all the tests.
In the literature, in order to study the fall of the disk, curves are given to show the evolution of the vertical velocity and the position of the center of the disk according to the time. We present the same analysis for different adaptation of our method. As expected, the disk reaches quickly a uniform fall velocity with slight moving on the right side of the vertical symmetry axis. This observation was already reported in the literature and is not specific to our method. One challenge is to propose robust methods that limit this breaking. In this section, we show that our method gives an answer to this question. Indeed, numerical simulations can be done with coarse meshes, even if it is not recommended with a fictitious domain approach (points on the interface can be far from the degrees of freedom introduced by the finite element method).
Contribution of the stabilization technique: Numerical tests are performed with the mesh , with and without performing the stabilization (). We compute the position of the disk according to the time and represent separately the vertical and horizontal positions.
In Figure 4 we represent the vertical velocity through the time. The results are similar, whether we perform stabilization or not. However, with the stabilization technique the method is more robust. If we zoom in (see Figure 4), without stabilization (red curve) some perturbations appear. When we compare the positions of the disk through the time, we do not observe difference on the vertical position in the Figure 5(a). However the difference is more important for the horizontal position and the rotation of the disk. The curves are plotted in red in Figure 6, with the mesh . With the stabilization technique the results are clearly improved. Without stabilization, the symmetry is broken even if it seems that the disk comes back around the symmetry axis of the cavity at the end of the simulation. This behavior related to the computation of the angular velocity which is represented in Figure 6(b).
Influence of the mesh size: With stabilization, we compare the simulations obtained with and . The results are given in Figure 5 (red curves for and blue curves for ). As expected, the smoothness of the solution is better when a sharper mesh is used. Besides, the result seems to be as good as results obtained in (7); (31) for instance. When a coarse mesh is used, the velocity is overestimated. This observation can be justified by the capability of the method for preserving the conservation of the mass. Indeed, with a coarse mesh, a numerical added mass appears in the system. This artificial mass is proportional to the stabilization terms (see and in the discrete problem) which are themselves proportional to the mesh size, and thus it can be neglected when the mesh size decreases.
The computation of the horizontal position of the disk is given in Figure 6(a). We observe that the disk tends to come back towards the symmetry axis of the cavity atthe end of the simulation with , unlike in the simulations with , where the symmetry breaking seems to growth. This phenomena can be observed in Figure 6(b) which represents the evolution of the rotation angle of the disk. With this angle always growths, unlike for the sharper mesh . This behavior can be justified by perturbation associated with our numerical method, in particular for the treatment of the nonlinear term. Moreover, a perturbation in horizontal velocity component is difficult to compensate during the simulation and contributes to the amplifying the phenomena. However, with stabilization, relevant values of the rotation are obtained and show that their influence is reduced compared with the translation.
5 Practical remarks on the numerical implementation.

All the numerical simulations were performed with the free generic library Getfem++ (32) (same source code for 2D and 3D) and implemented on High Performing Computers (parallel computations).

In order to compute properly the integrals over elements at the interface (during assembling procedure), external call to Qhull Library (33) is realized.

In the algorithm, steps 1– and 2– require computing of and , corresponding to integrations over the levelset. Such integrations require particular attention, in the sake of preserving a good accuracy. Indeed, the integrations must use nodes on levelset and accurate values on that nodes are required (no interpolation).

The method is very efficient in time computation, since it requires an update of the assembling matrices only locally near the interface.

As mentioned in (14), it is possible to define a reinforced stabilization technique in order to prevent difficulties that can occur when the intersection of the solid and the mesh over the whole domain introduce “very small” elements. The technique consists in selecting elements which are better to deduce the normal derivative on . A similar approach is given in (34). We think that this kind of reinforced stabilization technique can prevent the perturbations that appear during simulation (see zoom in Figure 4).
6 Conclusion
In this paper, we have considered a new fictitious domain method based on the extended finite element with stabilized term applied to the NavierStokes equations coupled with a moving solid. This method is quite simple to implement since all the variables (multipliers and primal variables) are defined on a single mesh independent of the computational domain. The algorithm leads to a robust method (good computation of the normal Cauchy stress tensor) whatever is the intersection of the domain with the  not necessarily sharp  mesh. The simulation of a falling disk with respect to the time confirms that our approach is able to predict well the interaction between the fluid and the structure. The stabilization must be considered to obtain more physical results preserving symmetry. Applications in 3D are in progress, in particular for control flow by acting on the boundary of the solid.
Acknowledgements
This work is partially supported by the foundation STAE in the context of the RTRA platform DYNAMORPH. It is based on the collaborative efforts with Alexei Lozinski and Yves Renard.
Footnotes
 michel.fournie@math.univtoulouse.fr
 journal: Journal of Fluids and Structures
References

G. Hou, J. Wang, A. Layton,
Numerical methods for
fluidstructure interaction—a review, Commun. Comput. Phys. 12 (2) (2012)
337–377.
doi:10.4208/cicp.291210.290411s.
URL http://dx.doi.org/10.4208/cicp.291210.290411s 
G. Legendre, T. Takahashi,
Convergence of a
LagrangeGalerkin method for a fluidrigid body system in ALE
formulation, M2AN Math. Model. Numer. Anal. 42 (4) (2008) 609–644.
doi:10.1051/m2an:2008020.
URL http://dx.doi.org/10.1051/m2an:2008020 
J. San Martín, J.F. Scheid, T. Takahashi, M. Tucsnak,
Convergence of the
LagrangeGalerkin method for the equations modelling the motion of a
fluidrigid system, SIAM J. Numer. Anal. 43 (4) (2005) 1536–1571.
doi:10.1137/S0036142903438161.
URL http://dx.doi.org/10.1137/S0036142903438161 
J. San Martín, L. Smaranda, T. Takahashi,
Convergence of a finite
element/ALE method for the Stokes equations in a domain depending on
time, J. Comput. Appl. Math. 230 (2) (2009) 521–545.
doi:10.1016/j.cam.2008.12.021.
URL http://dx.doi.org/10.1016/j.cam.2008.12.021 
C. S. Peskin, The immersed
boundary method, Acta Numer. 11 (2002) 479–517.
doi:10.1017/S0962492902000077.
URL http://dx.doi.org/10.1017/S0962492902000077 
R. Mittal, G. Iaccarino,
Immersed
boundary methods, in: Annual review of fluid mechanics. Vol. 37, Vol. 37
of Annu. Rev. Fluid Mech., Annual Reviews, Palo Alto, CA, 2005, pp. 239–261.
doi:10.1146/annurev.fluid.37.061903.175743.
URL http://dx.doi.org/10.1146/annurev.fluid.37.061903.175743  R. Glowinski, T.W. Pan, T. I. Hesla, D. D. Joseph, A distributed lagrange multiplier / fictitious domain method for particular flows, Int. J. of Multiphase Flow 25 (1999) 755–794.
 N. Moës, J. Dolbow, T. Belytschko, A finite element method for crack growth without remeshing, Internat. J. Numer. Methods Engrg. 46 (1) (1999) 131–150. doi:10.1002/(SICI)10970207(19990910)46:1<131::AIDNME726>3.0.CO;2J.

T.P. Fries, T. Belytschko, The
extended/generalized finite element method: an overview of the method and its
applications, Internat. J. Numer. Methods Engrg. 84 (3) (2010) 253–304.
doi:10.1002/nme.2914.
URL http://dx.doi.org/10.1002/nme.2914 
N. Moës, E. Béchet, M. Tourbier,
Imposing Dirichlet boundary
conditions in the extended finite element method, Internat. J. Numer.
Methods Engrg. 67 (12) (2006) 1641–1669.
doi:10.1002/nme.1675.
URL http://dx.doi.org/10.1002/nme.1675 
N. Sukumar, D. L. Chopp, N. Moës, T. Belytschko,
Modeling holes and
inclusions by level sets in the extended finiteelement method, Comput.
Methods Appl. Mech. Engrg. 190 (4647) (2001) 6183–6200.
doi:10.1016/S00457825(01)002158.
URL http://dx.doi.org/10.1016/S00457825(01)002158 
A. Gerstenberger, W. A. Wall,
An extended finite element
method/Lagrange multiplier based approach for fluidstructure interaction,
Comput. Methods Appl. Mech. Engrg. 197 (1920) (2008) 1699–1714.
doi:10.1016/j.cma.2007.07.002.
URL http://dx.doi.org/10.1016/j.cma.2007.07.002 
Y. J. Choi, M. A. Hulsen, H. E. H. Meijer,
An extended finite
element method for the simulation of particulate viscoelastic flows, J.
NonNewtonian Fluid Mech. 165 (1112) (2010) 607–624.
doi:doi:10.1016/j.jnnfm.2010.02.021.
URL http://dx.doi.org/doi:10.1016/j.jnnfm.2010.02.021 
J. Haslinger, Y. Renard, A new
fictitious domain approach inspired by the extended finite element method,
SIAM J. Numer. Anal. 47 (2) (2009) 1474–1499.
doi:10.1137/070704435.
URL http://dx.doi.org/10.1137/070704435 
N. Moës, A. Gravouil, T. Belytschko,
Nonplanar 3d crack growth by the
extended finite element and level sets, part i: Mechanical model, Internat.
J. Numer. Methods Engrg. 53 (11) (2002) 2549–2568.
doi:10.1002/nme.429.
URL http://dx.doi.org/10.1002/nme.429  F. L. Stazi, E. Budyn, J. Chessa, T. Belytschko, An extended finite element method with highorder elements for curved cracks, Comput. Mech. 31 (12) (2003) 38–48.
 N. Sukumar, N. Moës, B. Moran, T. Belytschko, Extended finite element method for threedimensional crack modelling, Int. J. Numer. Meth. Engng 48 (11) (2000) 1549–1570.

M. Stolarska, D. L. Chopp, N. Moës, T. Belytschko,
Modelling crack growth by level sets
in the extended finite element method, Int. J. Numer. Meth. Engng 51 (8)
(2001) 943–960.
doi:10.1002/nme.201.
URL http://dx.doi.org/10.1002/nme.201 
É. Béchet, N. Moës, B. Wohlmuth,
A stable Lagrange multiplier
space for stiff interface conditions within the extended finite element
method, Internat. J. Numer. Methods Engrg. 78 (8) (2009) 931–954.
doi:10.1002/nme.2515.
URL http://dx.doi.org/10.1002/nme.2515 
P. Hild, Y. Renard, A
stabilized Lagrange multiplier method for the finite element approximation
of contact problems in elastostatics, Numer. Math. 115 (1) (2010) 101–129.
doi:10.1007/s002110090273z.
URL http://dx.doi.org/10.1007/s002110090273z 
S. Court, M. Fournié, A. Lozinski,
A fictitious domain approach for
the Stokes problem based on the extended finite element method, Internat.
J. Numer. Methods Fluids 74 (2) (2014) 73–99.
doi:10.1002/fld.3839.
URL http://dx.doi.org/10.1002/fld.3839  J. Nitsche, Über ein Variationsprinzip zur Lösung von DirichletProblemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind, Abh. Math. Sem. Univ. Hamburg 36 (1971) 9–15, collection of articles dedicated to Lothar Collatz on his sixtieth birthday.

R. Becker, E. Burman, P. Hansbo, A
hierarchical NXFEM for fictitious domain simulations, Internat. J. Numer.
Methods Engrg. 86 (45) (2011) 549–559.
doi:10.1002/nme.3093.
URL http://dx.doi.org/10.1002/nme.3093 
E. Burman, P. Hansbo,
Fictitious domain finite
element methods using cut elements: II. A stabilized Nitsche method,
Appl. Numer. Math. 62 (4) (2012) 328–341.
doi:10.1016/j.apnum.2011.01.008.
URL http://dx.doi.org/10.1016/j.apnum.2011.01.008 
A. Massing, M. G. Larson, A. Logg, M. E. Rognes,
A stabilized Nitsche
fictitious domain method for the Stokes problem, J. Sci. Comput. 61 (3)
(2014) 604–628.
doi:10.1007/s1091501498389.
URL http://dx.doi.org/10.1007/s1091501498389  C. Airiau, J. M. Buchot, R. K. Dubey, M. Fournié, J. P. Raymond, Best actuator location for stabilization of the navierstokes equations, Submitted.
 L. C. Evans, Partial differential equations, 2nd Edition, Vol. 19 of Graduate Studies in Mathematics, American Mathematical Society, Providence, RI, 2010.

M. D. Gunzburger, S. L. Hou, Treating
inhomogeneous essential boundary conditions in finite element methods and the
calculation of boundary stresses, SIAM J. Numer. Anal. 29 (2) (1992)
390–424.
doi:10.1137/0729024.
URL http://dx.doi.org/10.1137/0729024 
H. J. C. Barbosa, T. J. R. Hughes,
The finite element
method with Lagrange multipliers on the boundary: circumventing the
BabuškaBrezzi condition, Comput. Methods Appl. Mech. Engrg. 85 (1)
(1991) 109–128.
doi:10.1016/00457825(91)90125P.
URL http://dx.doi.org/10.1016/00457825(91)90125P 
A. Ern, J.L. Guermond,
Theory and practice of
finite elements, Vol. 159 of Applied Mathematical Sciences, SpringerVerlag,
New York, 2004.
doi:10.1007/9781475743555.
URL http://dx.doi.org/10.1007/9781475743555  E. Hachem, S. Feghali, R. Codina, T. Coupez, Anisotropic adaptive meshing and monolithic variational multiscale method for fluid–structure interaction.

Y. Renard, J. Pommier, An open source
generic C++ library for finite element methods.
URL http://home.gna.org/getfem/ 
C. B. Barber, D. P. Dobkin, H. Huhdanpaa,
The quickhull algorithm for
convex hulls, ACM Trans. Math. Software 22 (4) (1996) 469–483.
doi:10.1145/235815.235821.
URL http://dx.doi.org/10.1145/235815.235821 
J. Pitkäranta, Local stability
conditions for the Babuška method of Lagrange multipliers, Math.
Comp. 35 (152) (1980) 1113–1129.
doi:10.2307/2006378.
URL http://dx.doi.org/10.2307/2006378