Bassem Barhoumi

An improved axisymmetric Convected Boundary Element Method formulation with uniform flow

Keywords: Axisymmetric Convected Boundary Element Method, Convected Helmholtz Equation, singular integrals, Finite Element Method, axisymmetric infinite rigid duct

This paper presents an improved form of the convected Boundary Element Method (BEM) for axisymmetric problems in a subsonic uniform flow. The proposed formulation based on the axisymmetric Convected Helmholtz Equation (CHE) and its fundamental solution that describes the sound radiation from a monopole source. The variables in the new axisymmetric boundary integral formulation can be expressed explicitly in terms of the acoustic pressure and its particular normal derivative. Also, the constant coefficients are expressed only in terms of the axisymmetric convected Green's function and its convected normal derivative. The particular and convected derivatives reduce the flow effects of the normal and the flow direction derivatives incorporated in the conventional convected boundary integral formulas. The advanced form of the axisymmetric boundary integral representation with flow is a similar form of the axisymmetric boundary element method without flow. Precisely, the two new operators significantly reduce the computational burden of the classical BEM and then becomes the CPU time of BEM without flow. The formula is verified comparing to both analytical and Finite Element Methods (FEM) of an axisymmetric infinite rigid duct in a subsonic uniform flow.

Introduction

The numerical studies of some acoustic problems with fluid flow require only the special methods, such as the computational techniques of the Finite Element Method (FEM) and boundary element method (BEM) for solving any axisymmetric medium in the science and engineering. The finite element method is currently the most widely used numerical method for solving the interior medium with mean flow and that the boundary element method is used for solving the exterior medium with a constant flow [START_REF] Astley | Numerical methods for noise propagation in moving flows, with application to turbofan engines[END_REF].

The finite element method has been implemented for solving the Helmholtz equation for the acoustic wave diffraction problems without flow [START_REF] Kechroud | Preconditioning techniques for the solution of the Helmholtz equation by the finite element method[END_REF]. It has been implemented with a displacement formulation for the modelling and simulation of the fluid-structure interaction problems in which the coupling conditions are formulated [START_REF] Gabard | A numerical method for vibroacoustic problems with sheared mean flows[END_REF]. A numerical technique of the FEM for the modeling of the acoustic pressure inside an axisymmetric lined duct in the presence of a constant flow was developed and presented in Ref. [START_REF] Taktak | Numerical modelling of the acoustic pressure inside an axisymmetric lined Flow Duct[END_REF]. This method is based on the resolution of the convected Helmholtz equation in the modal representation which gives good results compared with the analytical ones. The FEM applied to determine the time-harmonic acoustic field in a 2D infinite waveguide with walls covered of an absorbing material, and in the presence of a mean flow assumed uniform far from the source [START_REF] Redon | Non-reecting boundary conditions for acoustic propagation in ducts with acoustic treatment and mean flow[END_REF]. The boundary conditions of the rigid and lined ducts are represented by transparent boundary conditions suitable for use in a finite element scheme which based on the Dirichlet-to-Neumann (DtN) boundary condition [START_REF] Redon | Non-reecting boundary conditions for acoustic propagation in ducts with acoustic treatment and mean flow[END_REF]. For particular acoustic domains, [START_REF] Balin | Boundary Element and Finite Element Coupling for Aeroacoustics Simulations[END_REF] has determined the acoustic propagation in a transformed domain by the Prandtl-Glauert transformation.

It is known that the finite element method has its limitations in modeling infinite domains. The unbounded medium contains the non-reflection condition which is represented by the Sommerfeld radiation condition in a medium at rest [START_REF] Sommerfeld | Partial differential equations in physics[END_REF] or in a subsonic uniform flow [START_REF] Balin | Boundary Element and Finite Element Coupling for Aeroacoustics Simulations[END_REF]. Using the boundary element method (BEM), which requires a discretization of only the generator of the acoustic domain and that this Sommerfeld condition is automatically fulfilled by the fundamental solution. This makes the BEM superior to the FEM in many cases, especially for the unbounded problems, such as the wave propagation in the infinite lined ducts.

The boundary element method can be formulated in the transformed acoustic medium [START_REF] Casenave | Coupled BEM-FEM for the convected Helmholtz equation with non-uniform flow in a bounded domain[END_REF] and in the original acoustic medium [START_REF] Wu | A direct boundary integral formulation for acoustic radiation in a subsonic uniform flow[END_REF]. In order to take advantage of both the FEM and the BEM, a coupled FEM/BEM approach has been proposed [START_REF] Zhang | A coupled FEM/BEM formulation for acoustic radiation in a subsonic non-uniform flow[END_REF]. A direct collocation boundary element formulation for the calculation of acoustic propagation in a subsonic uniform mean flow has been presented in Ref. [START_REF] Juhl | Non-axisymmetric acoustics propagation in and radiation from lined ducts in a subsonic uniform mean flow[END_REF] in which the model is directed towards the calculation of the propagation inside and radiation from axisymmetric lined ducts. An advanced boundary element/fast Fourier transform (BE/FFT) methodology for solving axisymmetric acoustic wave scattering and radiation problems with non-axisymmetric boundary conditions has been developed in Ref. [START_REF] Stephanos | An advanced boundary element/fast Fourier transform axisymmetric formulation for acoustic radiation and wave scattering problems[END_REF].

However, the main disadvantages of using boundary integral formulations in an original acoustic medium are the same as in the non-flow case [START_REF] Juhl | An axisymmetric boundary element formulation of sound wave propagation in fluids including viscous and thermal losses[END_REF]. In addition, the convection effects in the presence of flow substantially increase these difficulties and makes significantly more complicated the conventional integral representations [START_REF] Zhang | A coupled FEM/BEM formulation for acoustic radiation in a subsonic non-uniform flow[END_REF][START_REF] Juhl | Non-axisymmetric acoustics propagation in and radiation from lined ducts in a subsonic uniform mean flow[END_REF]. These axisymmetric formulations are derived from the three-dimensional boundary integral formulation developed in [START_REF] Wu | A direct boundary integral formulation for acoustic radiation in a subsonic uniform flow[END_REF] in which the convected terms due to the Green's function and its derivatives such as the normal derivative and the derivative in the ow direction.

Also, the authors in Ref. [START_REF] Beldi | A simple and effective axisymmetric convected Helmholtz integral equation[END_REF] have been used a reformulation of the threedimensional convected boundary element to solve the axisymmetric acoustic problems. The reformulation based on the non-standard normal derivative similar to the non-flow normal derivative of the Green's kernel. The non-standard derivative leads to a global form of the integral formula contains another supplementary convec-tive terms in which make difficulty the numerical implementation of the boundary element method. The acoustic variable of the classical and reformulation axisymmetric boundary integral representations in Refs. [START_REF] Zhang | A coupled FEM/BEM formulation for acoustic radiation in a subsonic non-uniform flow[END_REF][START_REF] Juhl | Non-axisymmetric acoustics propagation in and radiation from lined ducts in a subsonic uniform mean flow[END_REF][START_REF] Beldi | A simple and effective axisymmetric convected Helmholtz integral equation[END_REF] in terms of the acoustic pressure as well as its normal and tangential derivatives.

The aim of this paper is to contribute to reduce the complexity of convected boundary integral formulas for axisymmetric domains in a subsonic uniform flow. Thus, a new reformulation of the boundary element method for the axisymmetric mediums (BEMA) which based on the axisymmetric convected Helmholtz equation and its fundamental solution is presented. This new axisymmetric boundary element method requires only two operators which substantially reduce the flow effects of the normal and the flow direction derivatives of the axisymmetric convected Green's function incorporated in the classical convected boundary integral formulations. Through the use of these new convected terms, the new axisymmetric BEM formula has a similar form of the axisymmetric BEM without flow. The acoustic variables in the new axisymmetric boundary integral equation can be expressed explicitly in terms of the acoustic pressure and its particular normal derivative.

The numerical implementation of BEMA is developed for the acoustic propagation in an axisymmetric rigid duct and validated by comparison with FEMA method and the analytical solution in which the CPU times are presented.

An improved Axisymmetric BEM Formulation

Consider an acoustic medium Ω with a compressible fluid, inviscid, isentropic, which is characterized by the Mach flow vector M ∞ , of density ρ ∞ and speed of sound c ∞ in the axisymmetric physical space (o,r,z). In figure 1, ∂Ω = Γ ∞ ∪ Γ is the generator limiting the axisymmetric medium Ω, such as Γ ∞ and Γ are the exterior and interior generators of the outgoing normal vector n q (n rq , n zq ) and the tangential vector t q (t rq , t zq ) at the point q, respectively. The generator ∂Ω not perturbed the streamlines and unchanged the direction of the flow.

Axisymmetric Convected Helmholtz Equation

A combination of the Axisymmetric Linearized Euler Equations can be produced the Axisymmetric Convected Helmholtz Equation in a harmonic regime (exp(+iωt)) in which t is the time, ω is the angular frequency and i is the imaginary unit. The wave equation wearing on axisymmetric acoustic pressure p is decomposed into a non-flow part and a convection part [START_REF] Zhang | A coupled FEM/BEM formulation for acoustic radiation in a subsonic non-uniform flow[END_REF]

[∆ + k 2 ](p) -(M ∞ .∇) (2ik + (M ∞ .∇))(p) = 0, (1)
where ∆ and ∇ are the axisymmetric Laplace and Gradient operators. k = ω/c ∞ is the original wavenumber. The fundamental solution of Eq. ( 1) is the axisymmetric convected Green's function G k 0 = 2π 0 G k dβ, which is the integration of the threedimensional convected Green function G k over the azimuthal angle β. The Green's kernel G k (m, q) is given for the points m and q by

G k (m, q) = exp -ik * mq.M * ∞ + R * β (4παR * β ), (2)
where m(r m , z m ) and q(r q , z q ) are the axisymmetric source-observer points, α is the Prandtl-Glauert factor, k * = k/α is the convected wavenumber, and that M * ∞ = M ∞ /α is the convected Mach vector. The distance R * β is the convected radius which depends on the azimuthal order β, it's given for the points m and q by

R * β (m, q) = r * 2 + 4r q r m sin 2 (β/2), (3)
in which r * (m, q) = mq 2 + (mq.M * ∞ ) 2 is the axisymmetric convected radius and mq is the physical distance between source-observer. Also, the axisymmetric convected Green's function satisfies the following inhomogeneous convected Helmholtz equation

H G k 0 (m, q) = ∆ q + k 2 G k 0 (m, q) + (M ∞ .∇ q ) (2ik -(M ∞ .∇ q )) G k 0 (m, q) = - 1 2πr q δ (q -m) , (4)
in which (•) q is the operator at point q and δ is the Dirac delta function such as δ(q -m) = δ(r q -r m )δ(z q -z m ).

Axisymmetric convected Helmholtz integral equation

To obtain a boundary integral equation formulation for the acoustic problem, the Helmholtz equation Eq. ( 1) is multiplied by the axisymmetric convected Green's function G k 0 (m, q) and its fundamental equation Eq. ( 4) by p(q). Subtracting these two equations and integrating over the axisymmetric exterior medium Ω yields

p (m) = Ω G k 0 (m, q) H (p (q)) -p (q) H G k 0 (m, q) r q dΩ q , m ∈ Ω\Γ, (5)
Using the axisymmetric second Green's formula, the axisymmetric Gradient property and the axisymmetric divergence theorem (Appendix A), the axisymmetric acoustic pressure Eq. ( 5) takes the following form

p (m) = ∂Ω G k 0 (m, q) ∂p ∂n q (q) -p (q) ∂G k 0 ∂n q (m, q)] r q dΓ q -2ik ∂Ω M ∞n p (q) G k 0 (m, q)r q dΓ q - ∂Ω M ∞n G k 0 (m, q) (M ∞ .∇ q )p (q) + ∂Ω M ∞n p (q) (M ∞ .∇ q )G k 0 (m, q)r q dΓ q , m ∈ Ω\Γ, (6)
where ∇ q = σ n n q + σ t t q is the axisymmetric Gradient operator in which the couple (σ n = ∂/∂ nq , σ t = ∂/∂ tq ) designates the normal and tangential derivatives at point q, respectively. M ∞n = M ∞ .n q is the normal Mach number.

In order to simplify the boundary integral equation Eq. ( 6), the normal derivative ∂/∂n q and the flow direction derivative (M ∞ .∇ q ) can be converted into a particular normal derivative d/dn q , which is given by the following relation

d(•) dn q = ∂(•) ∂n q -M ∞n (M ∞ .∇ q ) (•), (7)
The particular normal derivative of the axisymmetric convected Green function dG k 0 /dn q = 2π 0 dG k /dn q dβ is the integration of the particular normal derivative of the three-dimensional Green kernel dG k /dn q over the azimuthal angle β. Using the normal and the flow direction derivatives of the three-dimensional Green's function in Refs. [START_REF] Zhang | A coupled FEM/BEM formulation for acoustic radiation in a subsonic non-uniform flow[END_REF][START_REF] Juhl | Non-axisymmetric acoustics propagation in and radiation from lined ducts in a subsonic uniform mean flow[END_REF][START_REF] Beldi | A simple and effective axisymmetric convected Helmholtz integral equation[END_REF], one obtains the following derivative function

dG k dn q (m, q) = -G k (m, q) 1 + ik * R * β (m, q) r nβ R * β 2 (m, q) + ikM ∞n ] , (8)
where r nβ = mq axi .n q is the normal distance and that the vector mq axi of coordinate (r q -r m cos(β), z q -z m ). The particular normal derivative of the convected Green's function Eq. ( 8) is decomposed into two terms. The first term generalizes the expression of the normal derivative in a non-flow case and that the second term explains the convective influence. Substituting the particular normal derivatives Eqs. ( 7) and ( 8) in the integral formula Eq. ( 6) yield

p (m) = ∂Ω G k 0 (m, q) dp dn q (q) -p (q) d k G k 0 dn q (m, q)] r q dΓ q , m ∈ Ω\Γ, (9)
The axisymmetric boundary integral formulation Eq. ( 9) is decomposed into two integrals in which the acoustic variables are expressed only in terms of the acoustic field as well as its particular normal derivative which interpreting as a condition applied to the interior generator Γ. Thus, the axisymmetric BEM Eq. ( 9) is similar to the boundary integral equation solution for solving the axisymmetric potential flows [START_REF] Lennon | Boundary Integral Equation solution to axisymmetric potential flows[END_REF] and the non-flow case (the right first term in the equation 6) [START_REF] Kirkup | The Boundary Element Method in Acoustics[END_REF]. From the particular normal derivative operator Eq. ( 7), we obtain the convected normal derivative operator

d k (•) dn q = d(•) dn q + 2ikM ∞n (•), (10)
Substituting the particular normal derivative of the three-dimensional Green's function Eq. ( 8) into the convected normal derivative operator Eq. ( 10), one obtains

d k G k dn q (m, q) = -G k (m, q) 1 + ik * R * β (m, q) r nβ R * β 2 (m, q) -ikM ∞n ] , (11)
The integration over the circular generator Γ ∞ given by Eq. ( 9), represents contributions waves which come in the interior region due to the source is very far from observer point; q → ∞. Thus, this integral is null. According to the threedimensional Green's function at infinity [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF], the axisymmetric Green's kernel G k 0 is proportional to 1/ R * β . Also, taking into account the term r nβ /R * β is equivalent to the tangential Prandtl-Glauert factor β * t = β/ 1 -M 2 ∞t , the convected normal derivative dG k 0 /dn q of Eq. ( 11) is proportional to ik(1/ 1 -M 2 ∞t -M ∞n )/ R * β . This leads to two convected radiation conditions which verified by the axisymmetric acoustic pressure and its Green's kernel

lim r→∞ p(r) = 0, (12)
lim r→∞ √ r dp dr (r) + i k 1 -M 2 ∞t -kM ∞n p(r) = 0, (13)
where M ∞t = M ∞ .t q is the tangential Mach number. The first condition Eq. ( 12) indicates the non-reflection of waves, and the second condition Eq. ( 13) generalizes the classical Sommerfeld radiation condition in a uniform flow subsonic [START_REF] Sommerfeld | Partial differential equations in physics[END_REF][START_REF] Casenave | Coupled BEM-FEM for the convected Helmholtz equation with non-uniform flow in a bounded domain[END_REF]. Substituting the non-reflection conditions into the boundary equation Eq. ( 9), we obtain

p (m) = ± Γ G k 0 (m, q) dp dn q (q) -p (q) d k G k 0 dn q (m, q)] r q dΓ q , m ∈ Ω\Γ, (14)
in which the signs (+, -) designates the exterior and interior acoustic domains.

When the point m is taken on the generator at point q, the Green's kernel and its convected normal derivative Eqs. ( 2) and ( 11) contain singular and regular parts. Thus, Eq. ( 14) requires to study of behavior of the kernels near the singularity.

Singular integrals

A classical procedure can be used to isolate the singularity problems of the axisymmetric Green's function and its convected derivative. In fact, the singularities can be isolated by decomposing the axisymmetric convected Green's function G k 0 into a singular static part G 0 0 (k = 0) and regular part g k 0 depending on the wavenumber. Then, the axisymmetric Green's function takes the following form

G k 0 (m, q) = G 0 0 (m, q) + g k 0 (m, q) , (15)
Thus, the convected normal derivative is given by

d k G k 0 dn q (m, q) = d k G 0 0 dn q (m, q) + d k g k 0 dn q (m, q) , (16)
Because the difference between the singular and regular behavior, the singular functions in Eqs. ( 15) and ( 16) are evaluated analytically and that the regular parts are evaluated numerically with standard Gauss-Legendre quadrature of high order. In fact, the static axisymmetric Green's function is given by [START_REF] Zhang | A coupled FEM/BEM formulation for acoustic radiation in a subsonic non-uniform flow[END_REF][START_REF] Juhl | Non-axisymmetric acoustics propagation in and radiation from lined ducts in a subsonic uniform mean flow[END_REF]]

G 0 0 (m, q) = K (τ ) πατ * , (17)
Using the normal and the flow direction derivatives of the static Green's kernel in Ref. [START_REF] Zhang | A coupled FEM/BEM formulation for acoustic radiation in a subsonic non-uniform flow[END_REF], one obtains the following static convected derivative

d k G 0 0 dn q (m, q) = - 1 πατ * r nπ/2 r * 2 - n rq 2r q E (τ ) + n rq 2r q -ikM ∞n K (τ ) , (18)
where K(τ ) and E(τ ) are the complete elliptic integrals of the first and second kinds, the parameters τ and τ * are given by τ 2 τ * 2 = 4r q r m ,τ * = r * 2 + 4r q r m ,

According to the complete elliptic integrals at r → 0 [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF], the static Green's function exhibits a logarithmic behavior near the singularity, as follows G 0 0 (m, q) ∼ -

1 2πα Log (r * ) 2πr q , r → 0, (20)
Thus, the static convected normal derivative Eq. ( 18) near the singularity is characterized by

d k G 0 0 dn q (m, q) ∼ - 1 2πα 1 2πr q r n0 r * 2 , r → 0, (21)
Since r n0 /r * = O(r * ), then the singularity of the static convected function d k G k 0 /dn q is thus only apparent. Following the exclusion procedure [START_REF] Bonnet | Boundary Integral Equation Methods for Solids and Fluids[END_REF][START_REF] Guiggiani | A General Algorithm for Multidimensional Cauchy Principal Value Integrals in the boundary Element Method[END_REF] to isolate the singularity problems in the boundary formula Eq. ( 14), which can be rewritten as

p (m) = ± lim ε→0 Γε G k 0 (m, q) dp dn q (q) -p (q) d k G k 0 dn q (m, q)] r q dΓ q ± lim ε→0 Γ ∓ ε G k 0 (m, q) dp dn q (q) -p (q) d k G k 0 dn q (m, q)] r q dΓ q , ( 22
)
where s ε is the disk of center m and of radius ε = r bounded by the circle Γ + ε ∪ Γ - dp dn q (q) r q dΓ q ∓ Γ p (q) d k G k 0 dn q (m, q) r q dΓ q , m ∈ Ω, [START_REF] Zhang | A hypersingular integral formulation for acoustic radiation in moving flows[END_REF] where c ± (m) is the free term at any point m, which is given by

c ± (m) = 1 - θ ∓ (m) 2π = 1 ± Γ - 1 2πβ mq.n q r * 2 (m, q) dΓ q , ( 24
)
in which θ ± (m) is the convected angle. When the point m in the exterior and interior domains, the convected angle θ ± (m) = 0 and that the coefficient c ± (m) = 1. If m is a regular point of the generator, θ ± (m) = π and that c ± (m) = 0.5.

Numerical implementation

The numerical implementation of the axisymmetric boundary integral equation Eq. ( 23) for a body of arbitrary shape is obtained by discretizing the generator Γ of the body with N one-dimensional quadratic isoparametric elements, according to the following convergence criterion

n e h e ≤ λα 2 , ( 25
)
where n is the number of finite elements for a wavelength λ and he is the diameter of the finite element. Then, by using the quadratic shape functions, the acoustic pressure p and its normal derivative σ n of any point q on the generator Γ are assumed to be given in terms of the nodal coordinates. Thus, the tangential derivative in the local coordinates system along the generator Γ at point q is given by [START_REF] Juhl | Non-axisymmetric acoustics propagation in and radiation from lined ducts in a subsonic uniform mean flow[END_REF][START_REF] Mancini | An integral formulation for wave propagation on weakly non-uniform potential flows[END_REF]

σ t (q) = σ ξ (q) /J i [ξ (q)] , (26)
with -1 ≤ ξ ≤ 1 is the local coordinate and J i is the Jacobien of reference-local transformation. Using the collocation method so that the point m coincides with the nodes q i of the discretized generator Γ, Eq. ( 14) can be written in complex matrix form with Robin condition

[A] p = [B n ] σ n + [B c ] p c , (27)
where p, p c and σ n are the nodal vectors containing the values of the two-dimensional acoustic pressure p, its Dirichlet acoustic pressure p c and its Neumann normal derivative σ n at each node, respectively, whereas [A] and [B n ] are the axisymmetric acoustic matrices with the elements of matrix [A] contain the convected angles and elementary integrals related to the axisymmetric convected Green's function G k 0 and its convected normal derivative d k G k 0 /dn q , while the elements of matrix [B n ] is composed only by the elementary integrals containing the integrands of axisymmetric Green's function and its convected derivative. [B c ] is the acoustic matrix containing the integrands G k 0 and d k G k 0 /dn q .

Axisymmetric Finite Element Method

The axisymmetric Finite Element Method is based on the weak variational formulation, which can be obtained by multiplying the axisymmetric wave equation Eq. ( 1) by a test function p * (the axisymmetric acoustic pressure p * satisfies the axisymmetric convected Helmholtz equation Eq. ( 1)). Then, integrating the result over the axisymmetric interior medium Ω i , one obtains

Ω i p * ∆pdΩ i + k 2 Ω i p * pdΩ i -2ikM ∞ Ω i p * ∂p ∂z dΩ i -M 2 ∞ Ω i p * ∂ 2 p ∂z 2 dΩ i , (28)
Using the axisymmetric second Green's formula, one obtains

Ω i [p * ∆p] rdΩ i = - Ω i ∇p * ∇prdΩ i + Γ i p * ∂p ∂n Γ i rdΓ i , (29)
Ω i p * ∂ 2 p ∂z 2 rdΩ i = - Ω i ∂p * ∂z ∂p ∂z rdΩ i + Γ i p * ∂p ∂n z rdΓ i , (30)
2

Ω i p * ∂p ∂z rdΩ i = Ω i p * ∂p ∂z rdΩ i - Ω i ∂p * ∂z prdΩ i + Γ i p * pn z rdΓ i , (31)
Substituting Eqs. ( 29)-(31) into Eq. ( 28), we obtain a developed form of the weak variational formulation in Refs. [START_REF] Taktak | Numerical modelling of the acoustic pressure inside an axisymmetric lined Flow Duct[END_REF][START_REF] Reddy | An Introduction to Finite Element Methods[END_REF] -

Ω i ∇p * ∇prdΩ i + k 2 Ω i p * prdΩ i + ikM ∞ Ω i p * ∂p ∂z rdΩ i -ikM ∞ Ω i ∂p * ∂z prdΩ i + M 2 ∞ Ω i ∂p * ∂z ∂p ∂z rdΩ i + Γ i p * ∂p ∂n Γ i rdΓ i + ikM ∞ Γ i p * pn z rdΓ i -M 2 ∞ Γ i p * ∂p ∂n z rdΓ i = 0, (32)
The same classical numerical discretization in Ref. [START_REF] Taktak | Numerical modelling of the acoustic pressure inside an axisymmetric lined Flow Duct[END_REF], the elementary integrals of the numerical variational formulation Eq. (32) over the reference triangular or interval elements T take the following form

T f * Pd T, (33)
in which f and P are the regular polynomials of order n (n < 3) containing the shape functions. Then, these integrals are evaluated analytically [START_REF] Dhatt | Presentation of the finite elements method[END_REF]. Thus, the solve equation is a similar form of the complex matrix form Eq. ( 27) in which [A], [B n ] and [B c ] are the axisymmetric acoustic matrices with the various integrals in Eq. (32).

Application: Axisymmetric rigid duct

We consider an infinite axisymmetric rigid duct of truncated domain Ω of radius R = 0.1m, of length L = 0.5m, of symmetry axis (oz), in a subsonic uniform flow of the Mach vector M ∞ = M ∞ e z in the (r-z) meridian plane. The truncated medium is limited by the rigid wall Γ R of the normal vector n Γ R and by two generators, the left boundary Γ -at z = 0 and the right boundary Γ + at z = L of normal vectors n Γ -and n Γ + , respectively (see Figure 3). The acoustic pressure in an axisymmetric rigid duct is [START_REF] Gabard | A numerical method for vibroacoustic problems with sheared mean flows[END_REF]

p ± n = A n J 0 (k r,n r)e ik ± z,n z , (34)
in which (0) is the azimuthal order and (n) is the radial order, A n is the axisymmetric amplitude of order n equal to

√ 2 R 1 J 2 0 (k r,n R) -J 2 1 (k r,n R)
, J 0 is the Bessel function of the first kind and that k ± z,n and k r,n are the axial and radial wavenumbers (± indicates the propagation mode in the increased (+) and decreased (-) z-direction) which satisfies the following relations ∂J 0 (k r,n r)

∂r r=R = 0, ( 35
) k ± z,n = -kM ∞ ± k 2 -α 2 k 2 r,n α 2 , ( 36
)
Since the term as the function of tangential derivative σ t over the axisymmetric generator duct is zero, then the new radiation condition is given by the following particular derivative σ at right and left generators Γ ±

σ ± = ±iα 2 k ± z,n p (37)
We note k ± n = α 2 k ± z,n is the new wavenumber in which its module reduces the module of the classical wavenumber k ± z,n for the axisymmetric excitations.

Numerical methodologies

• The two-dimensional rigid duct is discretized by triangular elements structured and regular with three nodes isoparametric of two straight perpendiculars.

• The axial and radial wavenumbers satisfies the convergence criterion of the discretization Eq. ( 25) in which the number of triangles is 1652 elements (see figure 4).

• The evaluation of the acoustic pressure at point m in an interior axisymmetric rigid duct Ω given by Eq. ( 23) is numerically where Γ = Γ -∪ Γ R ∪ Γ + is the generator of this duct, and taking into account that the free term at the corner points m = Γ -∩ Γ R and m = Γ R ∩ Γ + is c -(m) = 3/4.

Numerical results

We apply BEMA, FEMA techniques and the analytical solution (S.Anal) associated to the axisymmetric acoustic pressure which propagated in the increased z-direction of axisymmetric excitation (p + n /A n ). Figures 5 and6 give the cartographies of the real parts Re(P + n /A n ) for different values of the Mach vector by S.BEMA, S.FEMA and S.Anal methods. Figures 5 and6 are presented the real parts of the axisymmetric acoustic pressure Re(P + n /A n ) by FEMA and BEMA which are in very good concord with analytical solution of relative error given by the L 2 norm. For FEMA, the relative error is less than 0.5% and for the axisymmetric boundary integral method, the error is less than 0.8%.

We observe that for the initial mode (0), the axisymmetric acoustic pressure is uniform along the radial and axial waves. Also, when the subsonic uniform flow is increased, the wavelength increases, which is explained by the offset of the wave in the positive z-direction (see Figure 5).

In the other mode, the axisymmetric acoustic pressure along the radial and axial axis is non-uniform in the presence and absence of the flow. Indeed, in the upstream region, the acoustic field is radially amplified to be greater than that obtained in the downstream region in which we observe the Sommerfeld radiation condition Eq. ( 12), and we notice also that the uniform flow is proportional to the wavelength (see Figure 6).

In addition, the CPU times of the Numerical boundary element method, the finite element method and that the classical boundary element method in Ref. [START_REF] Beldi | A simple and effective axisymmetric convected Helmholtz integral equation[END_REF] with mean flow are obtained by a machine 2.93 GHz using MATLAB-ACOUSTIC codes, which are given by the following A comparison between the CPU times of FE and BE methods in table 1, show that the axisymmetric Convected boundary element formulation Eq. ( 23) reduces 80% of the computation time of Finite element method Eq. (32).

Compared to the classical axisymmetric boundary integral formulae, the proposed method reduces a CPU time equal to 60% of the classical axisymmetric BEM. Also, the relative error for the classical axisymmetric boundary integral equation is less than 2%.

We observe that when the number of the discretization increase in term of the numerical implementation of the axisymmetric boundary element method with mean flow, the gain time between BEM-FEM and BEM-Classical BEM increase.

For the BEM case, we observe that the numerical execution time for the non-flow case is the numerical execution time for the flow case. This is because the improved form of the axisymmetric boundary integral representation Eq. ( 23) with the mean flow is similar to the boundary element method without flow.

Conclusion and prospects

In this work, an improved form of the axisymmetric Helmholtz integral equation and conventional finite element method for an axisymmetric medium in a subsonic uniform flow have been developed and presented. The axisymmetric convected boundary element method is derived from two new operators such as the convected and particular derivatives, which substantially reduces the complexity due to the presence of flow in the conventional formulations. The formulation is derived to be easy to implement as a numerical tool for computational codes of the axisymmetric acoustic medium in a mean flow. Also, the improved boundary integral representation significantly reduces the computational burden of the convenctional axisymmetric finite element method. Precisely, for the rigid duct case, the axisymmetric BEM require three Matrices in terms of integrals over the generator and that axisymmetric FEM require six assembly matrices in terms of integrals over the medium and its generator. The proposed methodologies can be applied to simplified the modal boundary element method and its derivatives [START_REF] Juhl | Non-axisymmetric acoustics propagation in and radiation from lined ducts in a subsonic uniform mean flow[END_REF][START_REF] Zhang | A hypersingular integral formulation for acoustic radiation in moving flows[END_REF].

Finally, we think that the method could be extended for the analysis of physical phenomena in other noise pollution systems such as compressors, aircraft engines, and ventilation systems that contain different materials. It is hoped that this could have a significant impact in design optimization of flow duct systems in the industry sector.

Figure 1 :

 1 Figure 1: Axisymmetric acoustic medium.

  ∪ Γ - ε ) for the exterior domain and that Γ = lim ε→0 (Γ ε ∪ Γ + ε ) for the interior domain (see figure2).

Figure 2 :

 2 Figure 2: Analytical regular method.

Figure 3 :

 3 Figure 3: Axisymmetric rigid duct.

Figure 4 :

 4 Figure 4: Axisymmetric mesh.

Figure 5 :

 5 Figure 5: Real part of the axisymmetric acoustic pressure Re(P + 0 /A 0 ) by S.Anal (top), FEMA (middle) and BEMA (bottom) for M ∞ = 0 (left) and M ∞ = 0.4e z (right) with k = 30.

Figure 6 :

 6 Figure 6: Real part of the axisymmetric acoustic pressure Re(P + 2 /A 2 ) by S.Anal (top), FEMA (middle) and BEMA (bottom) for M ∞ = 0 (left) and M ∞ = 0.4e z (right) with k = 90.

table 1

 1

	Mach vector Presented BEM Classical BEM Presented FEM
	M ∞ = 0	2.71 min	3.08 min	8.45 min
	M ∞ = 0.4e	2.73 min	6.92 min	13.86 min

Table 1 :

 1 The CPU times of BEM, Classical BEM and FEM.

Tel.:+216 96 374 048, +33629522089, E-mail adresse: bassembarhoumi@hotmail.fr

APPENDIX A. AXISYMMETRIC ACOUSTIC PROPERTIES

The axisymmetric Green's identity formula as a function of the test function p * is given by Ω (∆ q p * (q) p (q) -∆ q p (q) p * (q)) dΩ q = ∂Ω [p * (q) ∂n q (q) -p (q) ∂p * ∂n q (q)]r q dΓ q , (A.1)

where ∂/∂ nq is the normal derivative and that ∆ q is the axisymmetric Laplace operator in an acoustics problem r-z

Using the axisymmetric Green's identity formula Eq. (A.1) for the axisymmetric Greens function G k 0 (m, q), yields

∂p ∂n q (q) -p (q) ∂G k 0 ∂n q (m, q) r q dΓ q , (A.

3)

The axisymmetric Gradient property is given by

And that the axisymmetric Divergence theorem is given by (p (q) M ∞ ) .∇ q (M ∞ .∇ q ) G k 0 (m, q) = ∇ q . (M ∞ .∇ q ) G k 0 (m, q) p (q) M ∞ -(M ∞ .∇ q ) G k 0 (m, q) ((M ∞ .∇ q ) p (q)) , (A.5)


C Nabti

H Seba

Subgraph Isomorphism Search in Massive Graph Databases

Keywords: Subgraph isomorphism, graph query, massive graph databases, graph summarizing, modular decomposition

Subgraph isomorphism search is a basic task in querying graph data. It consists to find all embeddings of a query graph in a data graph. It is encountered in many real world applications that require the management of structural data such as bioinformatics and chemistry. However, Subgraph isomorphism search is an NPcomplete problem which is prohibitively expensive in both memory and time in massive graph databases. To tackle this problem, we propose a new approach based on concepts widely different from existing works. Our approach relies on a summarized representation of the graph database that minimizes both the amount space required to store data graphs and the processing time of querying them. Experimental results show that our approach performs well compared to the most efficient approach of the literature.

INTRODUCTION

Graphs are naturally used to represent data in several domains and applications such as protein interactions, regions of images, 2D and 3D shape recognition, etc. Graphs are also a first class solution for massive data representation. In fact, graph databases enable us to build interesting models that map closely to the problem domain. In this context, querying graphs is a fundamental problem studied by a large research community. Subgraph isomorphism search is the basic type of graph queries. Given a query graph Q and a data graph G, the subgraph isomorphism problem is to find all embeddings of Q in G. However, this problem is NP-complete [START_REF] Garey | Computers and Intractability: A Guide to the Theory of NP-Completeness[END_REF] and the problem of finding practical solutions for massive graph databases is a challenge. Subgraph isomorphism knows a lot of works that consist in finding the best way to match a query graph with the data graphs. We review here the main existing algorithms. For an exhaustive and detailed list of solutions, several surveys are available on the topic such as [START_REF] Lee | An in-depth comparison of subgraph isomorphism algorithms in graph databases[END_REF] and [START_REF] Gallagher | Matching structure and semantics: A survey on graph-based pattern matching[END_REF]. Existing solutions are generally classified into two main categories:

1. Exact approaches: in this case, we are interested by returning the list of subgraphs of G that match exactly the query Q.

2. Inexact or error-tolerant approaches: in this case, we are interested by a ranked list of of subgraphs of G that are most similar to the query Q.

In both approaches, algorithms may be optimal or approximate. Optimal algorithms return a correct and complete solution but have, generally, an exponential time complexity. Approximate algorithms may not find the correct/complete solution but guarantee a polynomial time complexity. The inexact approach knows a flourishing research activity in several application domains such as databases and pattern recognition relying on different tools such as genetic algorithms [START_REF] Khoo | Multiple relational graphs mapping using genetic algorithms[END_REF], neural networks [START_REF] Micheli | Neural network for graphs : A contextual constructive approach[END_REF], etc. Within this approach, we generally compute a distance between the graphs. This distance reflects the degree of similarity or dissimilarity of the graphs. We focus here on exact subgraph isomorphism search where existing algorithms can be classified within 3 categories:

1. Backtracking-based algorithms such as Ullman's [START_REF] Ullmann | An Algorithm for Subgraph Isomorphism[END_REF], VF2 [START_REF] Cordella | A (Sub)Graph Isomorphism Algorithm for Matching Large Graphs[END_REF], QuickSI [START_REF] Shang | Taming verification hardness: An efficient algorithm for testing subgraph isomorphism[END_REF], GraphQL [START_REF] He | Graphs-at-a-time: Query language and access methods for graph databases[END_REF], GADDI [START_REF] Zhang | Gaddi: Distance index based subgraph matching in biological networks[END_REF], and SPath [START_REF] Zhao | On graph query optimization in large networks[END_REF]. This approach constructs a space search tree whose internal nodes correspond to partial solutions and leaves correspond to embeddings. The first exact subgraph isomorphism algorithm is due to Ullman (Ullmann, 1976). Ullmann's basic approach is to enumerate all possible mappings of vertices between the two graphs in a depth-first tree-search. Each node at level i of the search-tree maps a vertex of the query to a vertex of the data graph. Each path from the root to leaf in the search-tree represents a complete mapping between the query and a subgraph of the data graph. This search-space increases exponentially with the size of the input graphs, so all the solutions proposed in this approach are based on some pruning rules that prevent developing non necessary paths.

2. Indexing-based algorithms such as Grep [START_REF] Shasha | Algorithmics and applications of tree and graph searching[END_REF], gIndex [START_REF] Yan | Graph indexing: A frequent structure-based approach[END_REF], FG-Index [START_REF] Cheng | Fg-index: Towards verification-free query processing on graph databases[END_REF], Tree+△ [START_REF] Zhao | Graph indexing: Tree + delta >= graph[END_REF], gCode [START_REF] Zou | A novel spectral coding in a large graph database[END_REF], SwiftIndex [START_REF] Shang | Taming verification hardness: An efficient algorithm for testing subgraph isomorphism[END_REF], and C-Tree [START_REF] He | Closure-tree: An index structure for graph queries[END_REF]. In this approach, indexes are used to minimize the number of candidate graphs in the database. Then, a subgraph isomorphism search is launched on the candidates.

3. Candidate Region selecting algorithms such as Turbo ISO [START_REF] Han | Turboiso: Towards ultrafast and robust subgraph isomorphism search in large graph databases[END_REF]. In this approach, the idea is to target specified regions on the same graph for subgraph isomorphism search. These regions are selected according to the properties of the query.

A candidate region for a query graph Q is a subgraph of the data graph G which may contain embeddings of the query graph. So, performing subgraph isomorphism search on all candidate regions will ensure that all embeddings can be obtained. However, minimizing the number of candidate regions and the size of each region is obviously important for faster matching. The main solution within this approach is called Turbo ISO [START_REF] Han | Turboiso: Towards ultrafast and robust subgraph isomorphism search in large graph databases[END_REF]. In order to minimize the size of each candidate region, the authors of [START_REF] Han | Turboiso: Towards ultrafast and robust subgraph isomorphism search in large graph databases[END_REF] propose to :

(a) Rewrite the query Q into an equivalent NEC (Neighborhood Equivalence Class) tree Q ′ . In Q ′ each set of vertices that have the same label and the same set of adjacent query vertices are merged into one NEC vertex. So, a NEC vertex is a compressed form of a set of vertices. Consequently, using Q ′ instead of Q, will accelerate the candidate region exploration process, since the number of vertices is smaller. (b) Construct candidate regions for the query Q in the data graph G by constructing for each region a BFS search tree T G from the root node u ′ s of the NEC tree Q ′ so that each leaf is on the shortest path from u ′ s . Then, for the start vertex v s of each target candidate region, identify candidate data vertices for each query vertex by simply performing depth-first search using T G and starting from v s .

Minimizing the number of regions comes through a careful choice of the root of the NEC tree. For this, Turbo ISO ranks every query vertex u by Rank(u) = f req (G,L(u)) deg (u) , where f req (G, l) is the number of data vertices in G that have label l, and deg(u) means the degree of u. This ranking function favors lower frequencies and higher degrees which will minimize the number of regions.

When exploring candidate regions, Turbo ISO also minimizes the number of enumerated partial solutions by ordering the NEC vertices by increasing sizes. Thus, paths involving fewer vertices are explored first, the space is pruned rapidly if no isomorphism is possible. In [START_REF] Han | Turboiso: Towards ultrafast and robust subgraph isomorphism search in large graph databases[END_REF], Turbo ISO is compared to the other approaches and its superiority in processing queries is attested via extensive experimentations.

Another solution that can be classified in this category is STW proposed in [START_REF] Sun | Efficient subgraph matching on billion node graphs[END_REF]. In this solution, the authors propose a graph decomposition into STwigs. An STwig is a two level tree structure, q = (r, L), where r is the label of the root node and L is the set of labels of its child nodes. STwigs are non overlapping star structures, i.e., edge disjoint stars.

Given a query graph Q, [START_REF] Sun | Efficient subgraph matching on billion node graphs[END_REF] first decomposes Q into a set of STwigs, then it uses exploration to find matches to each STwig. Candidate region exploration concerns the graphs that contain these STWings.

In this paper,we present a new approach for subgraph isomorphism search on large graphs. The proposed approach is completely different from all the previous approaches, and outperforms the most efficient existing algorithm in our experimentation. The main idea of the proposed approach is to enhance both time and space requirements of subgraph isomorphism search on large graphs. This is achieved by working on summarized graphs that are simpler and smaller than the original graphs.

The rest of the paper is organized as follows: Section 2 defines our notations and presents the compression algorithm used to summarize the graphs. Section 3 describes the proposed algorithm for subgraph isomorphism search on summarized graphs. Section 4 presents an experimental evaluation to show the effective performance of the proposed technique over the existing solutions. Finally, Section 5 concludes the paper with a summary of our work and its perspectives.

Basics. We consider data graphs defined as simple 1 vertex labeled graphs. We rely on the terminology used in [START_REF] Basu | Data Complexity in Pattern Recognition[END_REF][START_REF] Fan | Graph pattern matching: From intractable to polynomial time[END_REF].

Definition 1. 2 A data graph G is a 3-tuple G = (V, E, ℓ)
, where V is a set of nodes (also called vertices), E ⊆ V × V is a set of edges connecting the nodes, ℓ : V → Σ is a function labeling the nodes where Σ is the sets of labels that can appear on the nodes.

In this paper, the notation G = (V, E), with ℓ omitted means that we actually do not need the labels of the vertices but just their identifiers (i.e., indexes).

An undirected edge between vertices u and v is denoted indifferently by (u, v) or (v, u). For each v ∈ V , d(v) denotes the degree of v, i.e., the number of neighbors of v, where a neighbor is a vertex adjacent to v. The label or set of labels of a vertex v is given by ℓ(v).

The number of vertices of a graph is called the order of the graph. The number of edges of a graph is called the size of the graph. A graph that is contained in another graph is called a subgraph and can be defined as follows:

Definition 2. A graph G 1 = (V 1 , E 1 , f V 1 ) is a sub- graph of a graph G 2 = (V 2 , E 2 , f V 2 ), denoted G 1 ⊆ G 2 , if V 1 ⊆ V 2 , E 1 ⊆ E 2 , f V 1 (x) = f V 2 (x) ∀x ∈ V 1 .
Graph isomorphism is defined as follows:

Definition 3. A graph G 1 = (V 1 , E 1 , f V 1 ) and a graph G 2 = (V 2 , E 2 , f V 2 , f E 2 ) are said to be isomorphic, de- noted G 1 ∼ = G2, if there exists a bijective function h : V 1 → V 2 such that the following conditions are met: 1. ∀x ∈ V 1 : f V 1 (x) = f V 2 (h(x)) 2. ∀(x, y) ∈ E 1 : (h(x), h(y)) ∈ E 2 3. ∀(h(x), h(y)) ∈ E 2 : (x, y) ∈ E 1
Given a query graph Q and a data graph G, the subgraph isomorphism search of Q in G consists to find all the subgraphs of G that are isomorphic to Q.

Graph Summarizing. Generally graph summarizing methods designate graph compression methods that aim to reduce the amount of storage space required for storing a graph so that the processing of the graph does not require its decompression. So, this summaries must retain an amount of the graph properties that are sufficient to the application. Actually, this 1 Simple graphs are graphs with no edges involving a single vertex.

2 For presentation simplicity, we do not consider edge labels.

Table 1: Notation

Symbol Description G = (V, E, ℓ) undirected vertex labeled graph, f is a labeling function V (G) vertex set of the graph G E(G) edge set of the graph G G the complement of the graph G d(v) degree of vertex v G[X]
the subgraph of G induced by the set of vertices X

C (G) compressed graph of G root(C(G))
root of the tree corresponding to C (G)

Father(x)
the module that contains vertex (or module) x Leaves(x) set of vertices contained in module x ℓ(x) set of label of vertex (or module) x line of research is very new and knows little works. We can cite for example [START_REF] Chen | Mining graph patterns efficiently via randomized summaries[END_REF] where the authors propose to summarize a graph by grouping the vertices that have the same label into supervertices. In [START_REF] Fan | Query preserving graph compression[END_REF], the authors summarize graphs in a manner that preserves a class of queries, i.e., a query of this class returns the same result when applied to a graph G and when applied to the compression of G. They considered mainly reachability queries. In this kind of queries, we are interested in verifying if there is a path between two vertices within the graph. In [START_REF] Lagraa | A distance measure for large graphs based on prime graphs[END_REF], the authors propose a similarity measure between two large graphs based on a similarity measure between compressed versions of these graphs. They use modular decomposition [START_REF] Gallai | Transitiv orientierbare graphen[END_REF]Möhring, 1985a) to compress the graphs. A triangle listing algorithm is also proposed on graphs compressed by modular decomposition in [START_REF] Lagraa | An efficient exact algorithm for triangle listing in large graphs[END_REF]. These two applications show that modular decomposition is a promising compression method for graphs. A modular decomposition of a graph consists to find within the graph all the sets of vertices that share the same neighboringhood. These sets of vertices are called modules. In our framework, we also rely on modular decomposition to compress graphs. Modular decomposition is a graph representation method introduced by Gallai [START_REF] Gallai | Transitiv orientierbare graphen[END_REF] to solve optimization problems. It was also used to recognize some graph classes (Möhring, 1985a;[START_REF] Möhring | Algorithmic aspects of comparability graphs and interval graphs[END_REF][START_REF] Spinrad | Efficient Graph Representation[END_REF]. For a survey of applications of modular decomposition see [START_REF] Gallai | Transitiv orientierbare graphen[END_REF]Möhring, 1985a;[START_REF] Dahlhaus | Efficient and practical modular decomposition[END_REF]. The basic concept of modular decomposition, used in the compression process, is the notion of module defined as follows:

Definition 4. A module of a graph G = (V, E) is a set M ⊆ V of
vertices where all vertices in M have the same neighbors in V M.

A module M of G can take one of the following types:

• Series: if G[M] is a clique. • Parallel: if G[M] is a clique. • Neighborhood: Both, G[M]
and G[M] are connected graphs.

Figure 1 presents a graph and its modules. To retain all the properties of the original graph with the obtained compact representation of the graph, adjacency information for neighborhood modules must be stored. Series and parallel modules need no information about adjacency. For example, The obtained compressed graph illustrated in Figure 2(e) is itself a neighborhood module that can be denoted : N(a, S(b, c), N(d, e, f , g), h, P(S(i, j), k), P(l, m)).

For this module, we retain the edges between the supervertices to keep adjacency information. This gives the final compressed graph. We also retain the edges that bind the vertices of the neighborhood module N(d, e, f , g).

We note also that each module is a tree whose leaves are the vertices of the original graph as illustrated in Figure 3 for the modules of our example. Given a vertex v, we denote by Father(v) the module that contains it and by root(C G) the module corresponding to the compressed graph. Given a module m, Leaves(m) gives the leaves, i.e., vertices of the module m. Also, we will use ℓ(x) to denote the set of labels of a module or vertex x. Table 1 summarizes our notations.

Modular decomposition has been the subject of extensive research for years [START_REF] Gallai | Transitiv orientierbare graphen[END_REF]Möhring, 1985a;[START_REF] Dahlhaus | Efficient and practical modular decomposition[END_REF][START_REF] Habib | A survey of the algorithmic aspects of modular decomposition[END_REF][START_REF] Quaddoura | Classical graphs decomposition and their totally 2010 decomposable graphs[END_REF]. Several algorithms that compute the modular decomposition of a graph are proposed in the literature [START_REF] Habib | A simple linear-time modular decomposition algorithm for graphs[END_REF]. The most efficient are linear time [START_REF] Capelle | Graph decompositions and factorizing permutations[END_REF][START_REF] Habib | A simple linear-time modular decomposition algorithm for graphs[END_REF][START_REF] Tedder | Simpler linear-time modular decomposition via recursive factorizing permutations[END_REF] and achieves in O(n + m) .

SUM ISO : SUBGRAPH ISOMORPHISM SEARCH ON SUMMARIZED GRAPHS

As mentioned before, we propose here an algorithm that finds all the embeddings of a query graph Q in a data graph G. Both Q and G are compressed as described in the previous section. We show that by doing this we enhance both memory requirement for storing the data graphs and the processing time of the search. The algorithm, called Sum ISO , takes in entry the compressed versions C (Q) and C (G) of Q and G respectively and reports all the embeddings of Q in G. The Algorithm operates in two phases: a candidate supervertex selection phase and a subgraph search phase. During the first phase, the compressed data graph is parsed to retain only regions of the graph that are likely to contain the query. This selection uses only the labels of the modules. During the second phase, a backtracking-like algorithm is used in each region to verify the embedding. In the following we detail both phases and show how we can find all the embeddings by parsing the compressed graph data.

Candidate Supervertex Selection

The aim of this phase is to determine the modules (supervertices) that are likely to match the query. With this step, we minimize the number of vertices of the data graph to be processed. For this, we explore the vertices of C (G) to get all those that contain at least one of the labels of the query. Let Cand denotes the obtained result with:

Cand = {m ∈ C (G) such that ℓ(m) ∩ ℓ(C (Q)) ̸ = ∅}.
After that the set of candidate modules is partitioned in several subsets where each of them is candidate for a single embedding. Each subset contains the minimum number of modules that satisfy all the labels of the query. Subgraph search is then invoked on each of these subsets. Algorithm 1 details candidate supervertex selection. Figure 4 [START_REF] Lagraa | A distance measure for large graphs based on prime graphs[END_REF].

the query. The partitioning of Cand yields to the set {{1, 2}}. This means that there is only one possible region of the graph to explore for subgraph isomorphism.

Algorithm 1: Supervertex Selection.

Data: A summarized data graph C (G) and a summarized query

C (Q). Result: A set of candidate supervertices of C (G) that match C (Q). begin Cand ← / 0; foreach m ∈ C (G) do if ℓ(Q) ∩ ℓ(m) ̸ = ∅ then Cand ← Cand ∪ {m}; end end C ← {s = {m 1 , m , • • • , m j }|ℓ(Q) ⊆ ℓ(s)}; foreach s ∈ C do P ← / 0; SubgraphSearch(C (Q), s, P); end end
Note that at this step, we have a set of candidates with no order. These candidates are selected solely on labels. No structural verification are done with the query. So, at the end of this step, we do not know if there is a subgraph in G that matches the query. The aim of the next step is to aggregate the candidate supervertices in order to verify if the structure of the query is preserved within them.

SUBGRAPH SEARCH

The subgraph search phase takes as inputs a query

C (Q) and a set s = {m 1 , m 2 , • • • , m j } of modules that
are likely to contain an embedding of the query. It returns all the embeddings of the query in these modules. An embedding is represented by a set P of pairs (u, v), where u is a query vertex and v is the data vertex that match u. Algorithm 2 shows all the details of this phase. For each vertex u in C (Q), Subgraph-Search first finds the set of candidate vertices C u from the vertices of the modules of the set s. A vertex v of the data graph matches u if it has the same label as u and all the neighbors of u are matched to neighbors of v. This is verified by a call to function IsJoinable (detailed in Algorithm 3). Given two vertices u ( from the query) and v (from the data graph) to be matched, function IsJoinable returns TRUE if the neighbors of vertex u are matched to neighbors of vertex v in the match P. To have the list of neighbors of a vertex in a compressed graph, we use function Neighbors that takes advantage from the tree structure of the compressed graph ( see Figure 3 to easily list the neighbors of a vertex as detailed in Algorithm 4. According to the type (series, parallel or neighborhood) of the module that contains the vertex we can easily determine its neighbors. Algorithm 4 parses the subtree of C (G) that contains u from the father of u upward to the root of C (G). If a visited vertex x is a series module, then all the leaves of its descendants that are not in the branch that contains u are neighbors of u. If the visited vertex is a neighborhood module, neighbors of u are determined according to the edges of the module.

When a match (u, v) is verified, in procedure SubgraphSearch, it is reported in P. As in any backtracking-based algorithm, SubgraphSearch uses recursion to complete the partial match until it meets the query. When a match fails, the procedure backtrack to the preceding state by removing the match.

Algorithm 2: SubgraphSearch.

Data: A set of modules from the data graph s = {m 1 , m 2 , • • • , m j }, the compressed query C (Q) and a partial embedding P. Algorithm 3: Verify that two vertices to be matched have the same adjacency (IsJoinable).

Result: All embeddings of Q in s. begin if |P| = |V (C (Q))| then Report P; else Choose a non matched vertex u from Leaves(m), m ∈ C (Q); C u ← { non matched v ∈ Leaves(m i ) such that m i ∈ s and ℓ(v) = ℓ(u) and IsJoinable(u, v, P)}; foreach v ∈ C u do P ← P ∪ {(u, v)}; SubgraphSearch(C (Q),
Data: Two vertices u and v to be matched. Result: True is the the vertices have the same adjacency. begin return

(∀u ′ ∈ Neighbors(u), if u ′ is matched to v ′ then v ′ ∈ Neignbors(v)); end

EXPERIMENTAL RESULTS

We evaluate the execution time performance of our algorithm, Sum ISO , over different type of graphs and size of queries. We also compared it with the most efficient state of the art algorithm, called Turbo ISO and presented in [START_REF] Han | Turboiso: Towards ultrafast and robust subgraph isomorphism search in large graph databases[END_REF]. We recall that Turbo ISO is itself compared to the other existing solutions in [START_REF] Han | Turboiso: Towards ultrafast and robust subgraph isomorphism search in large graph databases[END_REF] and showed to be superior to them.

We first describe the datasets used in the experiments, then we present our results.

Datasets

We use the same datasets considered in [START_REF] Han | Turboiso: Towards ultrafast and robust subgraph isomorphism search in large graph databases[END_REF] for proving the superiority of Turbo ISO against the other algorithm of the literature described in Section 1. These datasets are referred to as AIDS, NASA, and Human. Their description is as follows:

• AIDS database: This dataset consists of graphs representing molecular compounds. It contains 10, 000 small graphs of 27 edges. The number of unique labels in AIDS is 51. We present a summary of these graph databases in Table 2. Besides the average number of vertices and edges of the graphs in the dataset, we also give the average compression rate of each dataset. Given a graph G and its compressed graph C (G), the compression rate of G is given by: CR Graphs within the three datasets were preliminarily compressed using an extension of the algorithm proposed in [START_REF] Capelle | Graph decompositions and factorizing permutations[END_REF][START_REF] Habib | A simple linear-time modular decomposition algorithm for graphs[END_REF] that computes the modular decomposition of a graph in linear time. So, we compress an input graph in O(n + m) time, where n is the number of vertices and m the number of edges of the graph.

(G) = |E(C (G)|)| |E(G)| • 100%.
To show the storage saving obtained by compressing the datasets, Table 3 reports the size on disk of each dataset before and after compression. We also provide the time necessary to compress each dataset. The experiments are performed on a 2.40 GHz Intel(R) Core(T M) i5 -4210U 64 bits laptop with 8 GB of RAM running windows 7. The algorithm is implemented in C++.

We use the same query sets as in [START_REF] Han | Turboiso: Towards ultrafast and robust subgraph isomorphism search in large graph databases[END_REF]. These queries are constructed as follows [START_REF] Han | Turboiso: Towards ultrafast and robust subgraph isomorphism search in large graph databases[END_REF]:

• AIDS and NASA query sets: For each of these datasets, [START_REF] Han | Turboiso: Towards ultrafast and robust subgraph isomorphism search in large graph databases[END_REF] constructed 6 query sets (Q4, Q8, Q12, Q16, Q20, Q24), each of which contains 1,000 query graphs of the same size. Additionally, each query Q i is contained in a query Q i+1 . Each query is a subgraph of a graph in the dataset.

• Human Query sets: For this dataset, [START_REF] Han | Turboiso: Towards ultrafast and robust subgraph isomorphism search in large graph databases[END_REF] generated three kind of queries:

1. Subgraph queries as for the Aids and Nasa datasets. In this case, we have 10 query sets obtained by varying the number of query sizes from 1 to 10. 2. Clique queries where the query subgraph is a complete graph. For biological datasets, such as Human, a clique Query corresponds to a protein complex [START_REF] He | Graphs-at-a-time: Query language and access methods for graph databases[END_REF]. 3. Path queries where the query subgraph is a path. A path query corresponds to transcrip- tional or signaling pathways [START_REF] He | Graphs-at-a-time: Query language and access methods for graph databases[END_REF].

The time performance reported in the results is the average time computed over the sets of queries of the same size.

Results

Figure 5 shows the experimental results for AIDS. We can clearly see that the time performed by Turbo ISO decreases when the query size increases. This is explained in [START_REF] Han | Turboiso: Towards ultrafast and robust subgraph isomorphism search in large graph databases[END_REF] by the containment relationship among the query sets in AIDS. We can also observe the same behavior with Sum ISO which achieves better than Turbo ISO . In our case, this can be explained by important compression rate of AIDS that yields a small number of candidates to be considered.

Figure 6 shows the experimental results for NASA. For this dataset, Sum ISO achieves significantly better than Turbo ISO for all the queries.

Figure 7 shows the results of subgraph queries over the human dataset. The superiority of Sum ISO over Turbo ISO is clearly observable as soon as the query size is greater than 8.

Figure 8 shows the results of subgraph isomorphism search for path and clique queries over the Human dataset. For the clique queries, Sum ISO significantly outperforms Turbo ISO . This is mainly due to the fact that a clique is compressed to a single node in our approach. For path queries, we have also a better results than Turbo ISO even if not significantly. We explain this by the fact that paths are not summarized by modular decomposition.

CONCLUSIONS

In this paper, we presented a new approach to subgraph isomorphism search in massive graph databases. In our approach, data graphs are summarized to minimize storage requirement. Our subgraph isomorphism search algorithm, Sum ISO , finds all the embeddings of a query graph in a summarized data graph without decompressing the graph. Our experimentations show that the proposed approach achieves good performance on both time processing of queries and space storage of data graphs. However, several enhancement issues are possible and merit investigation. First, more experiments are needed to attest the efficiency of the approach. In fact, we used the same datasets as in [START_REF] Han | Turboiso: Towards ultrafast and robust subgraph isomorphism search in large graph databases[END_REF] for our evaluation to compare with the Turbo ISO algorithm. These datasets, namely AIDS, NASA and HUMAN are highly compressible with more than 40% of compression rate. However, not all the graphs are compressible with the same rate. So, it is interesting to study the behavior of the approach behaves with less compressible graphs and determine a compression rate threshold that conditions the use of compression. Second, it is interesting to compare the two approaches on larger graphs. Third, it is interesting to see if it is feasible to ran such an approach on a graph database such as Neo4j and also investigate how it can be implemented in a MapReduce framework. Finally, we have not used pruning methods in the SubgraphSearch phase and it may be possible to define some rules to prune the explored compressed paths by relying on the properties of the compression.

Acknowledgments

This work was supported by a grant from the ANR CAIR project.

Figure 1 :

 1 Figure1: A Graph and its Modules[START_REF] Lagraa | A distance measure for large graphs based on prime graphs[END_REF].The graph is compacted by replacing recursively each module by a supervertex as illustrated in Figure 2. To obtain a unique representation of the graph only the modules that do not overlap other modules are considered.To retain all the properties of the original graph with the obtained compact representation of the graph, adjacency information for neighborhood modules must be stored. Series and parallel modules need no information about adjacency. For example, The obtained compressed graph illustrated in Figure2(e) is itself a neighborhood module that can be denoted :

Figure 2 :Figure 3 :

 23 Figure 2: Compressing steps: S: series module. P: parallel module. N: neighborhood module (Lagraa et al., 2014).

Figure 4 :

 4 Figure 4: Example of Candidate Supervertex Selection.

  Figure 5: Aids dataset.

  Figure 7: Human dataset.

Figure 8 :

 8 Figure 8: Path and Clique Queries.

  Computing the set of neighbors of a vertex in a compressed graph (Neighbors).Data: A vertex u and a compressed graph C (G).

	Algorithm 4: Result: The set of neighbors of u in g.
	begin
	N ← / 0;
	z ← u;
	x ← Father(u);
	while x ̸ = root(C (G)) do
	switch type of x do
	case a series module
	foreach child y ̸ = z of x do
	N ← N ∪ Leaves(y);
	end
	case a Neighborhood module
	foreach edge (z, y) ∈ x do
	N ← N ∪ Leaves(y);
	end
	endsw
	z ← x;
	x ← Father(x);
	end
	return N
	end
	• HUMAN database: This dataset consists of one
	large graph representing a protein interaction net-
	work. This graph has 4, 675 vertices and 86, 282
	edges. The number of unique labels in the dataset
	is 90.
	• NASA database: This dataset contains 36, 790
	trees with an average size of 32, and a number
	of unique labels of 117, 302.

  It compares the number of edges in C (G) in respect to G.

	Dataset	Number of	avg|V |	avg|E|	CR
		graphs
	AIDS	10,000	26	27	56.8%
	NASA	36,790	94	32	44.2%
	HUMAN	1	4,675	86,282	61%

Table 2 :

 2 Graph Dataset Characteristics. avg|V |: average number of vertices. avg|E|: average number of edges. CR average compression rate.

Table 3 :

 3 Size on disk.

	Dataset	Size on disk (Mb)	Size on disk after	Compression
			compression (Mb)	time(ms)
	AIDS	4.59	2.18	230
	NASA	24	14.42	180
	HUMAN	1.15	0.24	195