Implementation of the Immersed Boundary method in Code_Saturne
The following describes a proposal for implementation of an immersed boundary method into Code_Saturne. It is following the approach detailed in [1]. We consider here a simplified framework: 2D and fixed geometry (.i.e. fixed mesh and lagrangian points)
Variables and arrays used

nsupmax
maximum number of points in the support of a lagrangian point

nls
number of lagrangian points

xlag(nls)
, ylag(nls)
x and y coordinates of the lagrangian points

ds(nls)
arc step length

hx(nls)
, hy(nls)
dilatation parameters in the x and y direction

bmat(6,nls)
matrix of the RKPM method (6 because we assume we're in 2D)

nsup(nls)
number of points of the support of a given lagragian point

supelt(nsupmax,nls)
element index (i.e. Code_Saturne "iel") of the support

epsmat(nls)
coefficients
Steps of the method
 definition of the lagrangian points:
nls
, xlag
, ylag
, ds
 definition of the support for each lagragian point
 we find the nearest cell center I using findpt
 we define the neighborood as the set of cells which share a node with the cell I
 we define
hx
and hy
using the definition given in Eq.32, Eq.33 and Eq.34 of [1]
 we calculate the support using the Eq.35 of [1] (we choose and for now)
 we then define
nsup
and supelt
 calculation of the RKPM matrix
b
for each lagrangian points
 we calculate the moment matrix M (terms rescaled by the dilitation factors to avoid illconditionning), it needs the regular_delta function and the list of powers as defined in Eq.28 of [1]
 we calculate the determinant of M
 if det(M)=0 we set b=[1 0 0 0 0 0]
 if det(M)<>0 we perform a Cholesky decomposition of M = LL'
 we find the vector y1 so that L*y1=[1 0 0 0 0 0]
 we find the vector b so that L'*b=y1
 b is multiplied by H (rescaling matrix defined in Eq.37 of [1])
 calculation of the
epsilon
coefficients
 calculation of the terms of the matrix A (as defined in Eq.44 of [1]): we only consider the 5 righthand and lefthand neighbours of a given lagrangian point
 the system M*epsilon = [1 ..... 1] is solved. A Krylovtype iterative method is suggested for the resolution, no preconditionning required [1] (lagrangian points spacing ~ eulerian grid spacing).
Linear algebra
The following LAPACK subroutines are used:
 SDOT
 SPOFA (Cholesky decomposition)
 SPOSL (solves the real symmetric positive definite system using SPOFA)
 SAXPY (ax+y vector operation)
 SGEFA (LU decomposition)
 SGESL (solves a general system of linear equations, after factorization by SGEFA)
References
[1] Pinelli, A., I.Z. Naqavi, Piomelli, U. and Favier, J. Immersedboundary methods for general finitedifference and finitevolume NavierStokes solvers. Journal of Computational Physics, Vol. 229, pp. 90739091 (2010)