Christophe Droz

Régis Boukadia

Wim Desmet

A multi-scale model order reduction scheme for transient modelling of periodic structures

Keywords: Periodic structure, Transient Modelling, Wave Finite Element, Model Order Reduction, Multiscale, Floquet-Bloch theory

In this paper, a Floquet Model Order Reduction (MOR) method is proposed for the modelling of finite periodic structures subjected to harmonic or transient loads. This Floquet MOR generates frequency-independent projection matrices of the full-scale structure using a Floquet expansion of the local state-vectors computed on a single unit-cell with the Forced Wave Finite Element framework. This two-step MOR method can handle large multi-scale waveguides regardless of the scale and the number of unit-cells. A comparison with a standard sub-structuring approach based on unit-cell wave-mode reduction demonstrates that Floquet MOR is able to create reduced models of finite periodic structures with an unequaled computational efficiency. The method's ability to produce small and accurate models suitable for time-domain simulations presents a high potential for virtual sensing, real-time Structural Health Monitoring (SHM) and digital twins for large periodic structures and metamaterials.

Introduction

The requirement for high-fidelity models describing the structural dynamics and acoustic behavior of engineering structures increases with our society's growing concern for human safety and comfort. Indeed, the development of scalable materials and meta-structures capable of combining high structural integrity with vibroacoustic performances has been exerting a considerable pressure on the computational requirements in fields such as model-based Noise and Vibration Harshness (NVH) and predictive Structural Health Monitoring (SHM). To address these computational performances issues, Model Order Reduction (MOR) and sub-structuring techniques have been used for dynamic systems as early as in the 60s with the works of Hurty [START_REF] Hurty | Dynamic analysis of structural systems using component modes[END_REF] on modal decomposition in order to produce displacement transformation matrices for Component Mode Synthesis (CMS). Reviews by Besselink et al. [START_REF] Besselink | A comparison of model reduction techniques from structural dynamics, numerical mathematics and systems and control[END_REF] and Klerk et al. [START_REF] De Klerk | General framework for dynamic substructuring: history, review and classification of techniques[END_REF] provide a history and relevant classification techniques for general sub-structuring and MOR techniques.

Email address: christophe.droz@gmail.com (Christophe Droz) Nowadays, an emphasis is being put on time domain applications of MOR and dynamic sub-structuring techniques. A detailed survey of projection-based MOR was proposed by Lulf et al. [START_REF] Lülf | Reduced bases for nonlinear structural dynamic systems: A comparative study[END_REF] for discrete non-linear problems.

A typical example is the interaction of a dynamic system with a 'hard' non-linearity (e.g. contact, friction), which may induce chaotic responses that cannot generally be studied with non-linear frequency methods (e.g harmonic balance method (HBM)). In damage detection and broader SHM applications, Fourier transform can sometimes be used to unlock efficient frequency-domain MOR, but remains mostly limited to steady-state linear (or linearizable) systems with time invariant parameters. The same limitations exist for the development and exploitation of digital twins. Kapteyn et al. [START_REF] Kapteyn | Data-driven physics-based digital twins via a library of component-based reducedorder models[END_REF] used component-based MOR to build a physics-based near-real-time digital twin, aiming for a drone's ability to perform autonomous SHM. Time-domain MOR is also critical in the field of virtual sensing, were real-time state/input estimations methods are needed for certain vibroacoustic applications (e.g. see [START_REF] Van Ophem | Model based virtual intensity measurements for exterior vibro-acoustic radiation[END_REF][START_REF] Van De Walle | Virtual microphone sensing through vibro-acoustic modelling and kalman filtering[END_REF][START_REF] Huang | Improved kalman filter damage detection approach based on lp regularization[END_REF]). Such applications are particularly challenging when the models' outputs require a high level of spatial resolution across a large-scale structure. It is then crucial for a MOR strategy to produce small and accurate models which are actionable via direct time integration.

Modal projection methods such as CMS are popular for reducing large dynamic systems, as the resulting Reduced-Order Model (ROM) is generally exploitable in both frequency-and time-domains. When a large structure is made of an assembly of many sub-components, multi-scale MOR methods are needed to circumvent the computation of global modes and solve smaller eigenvalue problems. In model-based SHM and vibro-acoustic applications, multi-scale problems are often due to the co-existence of a macroscopic structure and a periodic architecture at the micro-or mesoscale (e.g. meta-materials, lattices, thin-walled or stiffened plates). To exploit their repetitive nature, periodic structure theory has been used to compute the dynamic response of periodic media based on a refined finite element (FE) description of their unit-cell [START_REF] Mead | A general theory of harmonic wave propagation in linear periodic systems with multiple coupling[END_REF]. The Wave Finite Element Method (WFEM) is a popular implementation allowing the computation of the finite response of loaded waveguides under different boundary conditions. Its interest lies in the reduction of the overall system into a considerably smaller set of propagating and evanescent wave amplitudes, which are resulting from a Floquet expansion of the displacements and forces state-vectors. The final system's dimension is independent of the number of unit-cells (UCs) in the structure. A reduced formulation of this Forced WFEM was investigated in Boukadia et al. [START_REF] Boukadia | A bloch wave reduction scheme for ultrafast band diagram and dynamic response computation in periodic structures[END_REF] to explore wave localization effects in a stiffened plate. However the main drawback of this Forced WFEM framework is the underlying periodicity assumption of the Floquet expansion, which makes the WFEM incompatible with most valuable time-domain applications. Recently some authors [START_REF] Mallouli | Damage detection on composite beam under transverse impact using the wave finite element method[END_REF][START_REF] Malik | Transient ultrasonic guided wave simulation in layered composite structures using a hybrid wave and finite element scheme[END_REF] used Inverse Fourier Transform to achieve time domain analyses of periodic structures, thus using the method's outstanding performances for harmonic analyses to predict their transient response. Unfortunately, this approach is not compatible with real-time model estimation and updating, which are crucial for digital twin or virtual sensing applications.

Modal projection techniques have been used at the unit-cell scale in combination with the Periodic structure theory, to reduce the size of spectral eigenvalue problems and accelerate the computation of wave dispersion characteristics.

Different techniques were proposed to build the local projection matrices: a combination of a component modes with a subset of least evanescent wave vectors were used by Droz et al. [START_REF] Droz | A reduced formulation for the free-wave propagation analysis in composite structures[END_REF][START_REF] Droz | A hybrid wave-mode formulation for the vibro-acoustic analysis of 2d periodic structures[END_REF] to reduce the interior and boundaries of the periodic unit-cells. As an alternative to wave vectors, Krattinger and Hussein [START_REF] Krattiger | Generalized bloch mode synthesis for accelerated calculation of elastic band structures[END_REF] adapted the CMS to Floquet-Bloch conditions and extended interior CMS reduction with boundary modes. A variety of similar techniques were proposed to generate a ROM for periodic UCs in order to compute their dispersion characteristics with 1D or 2D periodic finite element frameworks [START_REF] Palermo | A reduced bloch operator finite element method for fast calculation of elastic complex band structures[END_REF][START_REF] Beli | A projection-based model reduction strategy for the wave and vibration analysis of rotating periodic structures[END_REF][START_REF] Fan | Model reduction schemes for the wave and finite element method using the free modes of a unit cell[END_REF]. These methods are often qualified as 'UC MOR', as the local modes and waves derived from the UCs are used to build a UC projection matrix, but do not take into account the global loads or boundary conditions of the overall structure. Still, a local MOR approach can be used to assemble a global projection matrix (see for example [START_REF] Krattiger | Interface reduction for hurty/craig-bampton substructured models: Review and improvements[END_REF][START_REF] Pinault | A general superelement generation strategy for piecewise periodic media[END_REF]), later used to reduce the finite structure and compute its dynamic response under various boundary conditions. For example, Hu and Oskay [START_REF] Hu | Spectral variational multiscale model for transient dynamics of phononic crystals and acoustic metamaterials[END_REF] used a truncated local modal basis to create a frequency independent projection matrix. The global system built by assembling the local projection matrices was then used to conduct transient simulations in a large-scale phononic meta-material. Although these approaches, referred as 'Local MOR' in this paper, only require a single UC to build the projection matrix, the periodic nature of the structure is not fully exploited at the macro-scale, and the dimension of the global ROM always grows with the number of UCs.

To summarize, two approaches based on UC analysis have been investigated to build the ROM of periodic waveguides, suitable for time-domain simulations: (i) the Periodic Structure Theory, which is a highly efficient way to reduce a periodic system regardless the number of UCs. However it relies on a frequency-dependent projection inherent to the state vector's expansion onto Floquet vectors. (ii) the 'Local MOR' approach, which takes advantage of local modes and wave solutions to build a reduced UC ROM, later assembled to create a global projection matrix.

However the dimension of the global system depends on the number of UCs hence it does not fully take advantage of the structure's periodicity.

In this paper, a two-step strategy is proposed for generating a macro-scale ROM, actionable in time-domain while being based on Floquet theory and reduced UC analysis. Sections 2.1 and 2.2 remind the formalism and highlight the methodological and conceptual differences between Floquet expansion and local MOR methods for computing the dynamic response of finite periodic structures. In section 2.3 a Floquet MOR approach is introduced to generate global projection matrices from the expansion of forced state vectors computed using reduced UCs. The pros and cons of local MOR and Floquet MOR are discussed in section 3, with a thorough comparison of their numerical performances with regards to model size and UC scaling. In section 4 the Floquet MOR and an implicit integration scheme are used to study the impulse response of a periodic structure with 2D elliptic inclusions. The possible limitations of the method when it comes to high-frequency contents are stressed and conclusions are drawn in section 5.

Formulation of the model-order reduction scheme for periodic structures

Wave Finite Element Framework for dynamic analyses

Let us first remind the numerical formalism used to model a finite periodic structure based on Floquet theory and UC finite elements. Consider a linear elastic structure made of repetitive cells along one direction. The governing equation of the overall system writes:

M ∂ 2 U(t) ∂t 2 + C ∂U(t) ∂t + K U(t) = F (t) (1)
where cursive M, C and K stand for the assembled FE mass, damping and stiffness matrices of the system, while U(t) and F (t) are the global nodal displacements and external forces. Assuming the structure is made of 'N' identical cells with compatible meshes, the FE model of a single sub-structure of index 'n' can be used to describe the force equilibrium in the absence of inner loads:

                 M LL M LI M LR M IL M II M IR M RL M RI M RR                  ∂ 2 ∂t 2                  U (n) L U (n) I U (n) R                  (t)+                  C LL C LI C LR C IL C II C IR C RL C RI C RR                  ∂ ∂t                  U (n) L U (n) I U (n) R                  (t)+                  K LL K LI K LR K IL K II K IR K RL K RI K RR                                   U (n) L U (n) I U (n) R                  (t) =                  F (n) L 0 F (n) R                  (t) (2)
where M, C and K are the assembled FE mass, damping and stiffness matrices of the UC, while the sub-scripts 'L', 'I', 'R' stand for left-, inner-and right-sided degrees of freedom (DOF), respectively. Uppercase U(t) and F(t) denote the cell's nodal displacements and external forces.

In harmonic regime, the governing equation of the global system writes D(ω)U(ω) = F (ω), where

D = K + jωC -ω 2 M.
Using Floquet theory provides the solution U(x) of a linear system of equations H(x)U(x) = 0 where H is a periodic operator such as H(x + d) = H(x). The solution can decomposed with Floquet vectors as:

U(x) = V(x)e κx
where κ is a constant matrix and V is an invertible d-periodic matrix. In the following, Ψ denotes the wave matrices and λ = e κd are the Floquet coefficients (i.e. propagation constants), used to relate the consecutive state-vectors of the UC's displacements:

         u (n+1) L u (n+1) I          = λ          u (n) L u (n) I          (3)
Introducing the periodicity conditions between two consecutive UCs, using the above mentioned state-vectors yields the following spectral eigenproblem:

D(λ)          u L u I          =          0 0          ( 4
)
where D is the cell's dynamic stiffness matrix and:

D(λ) : D -→           D RL D RI 0 0           1 λ +           D RR + D LL D LI D IL D II           +           D LR 0 D IR 0           λ
Note that with the WFEM, a dynamic condensation is usually achieved to simplify the quadratic eigenvalue problem (EVP) into the common expression: A variety of numerical schemes [START_REF] Patel | On computing the eigenvalues of a symplectic pencil[END_REF] can be used to solve Eq.( 5). One of the major features of WFEM lies in the straightforwardness of boundary conditions implementation for determining the finite waveguide's dynamic response.

λD LR + (D LL + D RR ) + 1 λ D RL u L = 0 (
Indeed, the multi-modal wave expansion can be written to decompose any state vector's displacements u (n) L using the solutions of Eqs.( 5) or (4), sorted into positive-and negative-going (ψ + , ψ -) eigenvectors and their respective propagation constants as:

u (n) = S s=1 λ n s ψ + s q + s + λ N-n s ψ - s q - s (6)
where q + and q -are the complex wave amplitudes and 'S' the number of solutions retained in the decomposition.

Λ ω Λ -1 (ω)

Ψ -(ω) Ψ + (ω) F 0 F N Floquet expansion: u (n) (ω) q + (ω), q -(ω)
Dispersion analysis: D(ω) Ψ(ω), Λ(ω) Considering the free-free waveguide configuration, the force equilibrium on both sides subjected respectively to F 0 and F N yields the following system of equations, where Ψ are the concatenated eigenvectors and Λ the diagonal eigenvalues matrix:

u (0) u (n) u (N) u (N-1) u (n+1) u (1)
          D LL Ψ + + D LR Ψ + Λ D LL Ψ -Λ N + D LR Ψ -Λ N-1 D RL Ψ + Λ N-1 + D RR Ψ + Λ N D RL Ψ -Λ + D RR Ψ -                    q + q -          =          F 0 F N          (7)
Solving Eq.( 7), nodal displacements can be derived at any cell (n) using Eq.( 6). Note that the size of Eq.( 7) is twice the number of DOF in a state vector and does not depend on the number of UCs in the overall structure. It is also worth mentioning that the computation effort mainly lies in the determination of the frequency-dependent wavemodes Ψ and propagation constants Λ. Then, the dynamic response of the periodic structure with various boundary conditions or number of UCs can be determined.

Local MOR based on reduced unit-cell analysis

The above formulation allows a drastic reduction of the dynamic system's size. However, the Floquet expansion is based on the spatial periodicity of the dynamic operator, which is only valid in harmonic regime. The projection is therefore based on frequency dependent entities hence is not suitable for the reduction of the system in time domain.

In the following, the local modes of the UC are used to build a global projection matrix of the periodic structure.

This corresponds to the so-called Local MOR and can be seen as a multi-scale sub-structuring technique taking advantage of the periodicity to reuse a same reduction basis for each UC. The governing equation of each UC (i.e.

Eq.( 2)) is reduced by creating a local projection matrix Γ such that Eq.( 2) becomes:

                                 Γ * MΓ ∂ 2 p (n) (t) ∂t 2 + Γ * CΓ ∂p (n) (t) ∂t + Γ * KΓ p (n) (t) = Γ * F (n) (t)                      u (n) L u (n) I u (n+1) L                      = Γ p (n) (t) (8)
Several reduction techniques were proposed to find a suitable local projection matrix to reduce the UCs, which can be seen a specific application of more general modal sub-structuring approaches. The periodic nature of the media provides a straightforward dividing where either normal and constrained modes or even wave-modes can be used to describe the inner cells' dynamics.

In the following a reduction of the UCs is achieved using a two-step approach. In harmonic regime, Eq.( 2) writes:

                 D LL D LI D LR D IL D II D IR D RL D RI D RR                  (ω)                  u (n) L u (n) I u (n+1) L                  (ω) =                  f (n) L 0 -f (n+1) L                  (ω) (9)
The inner displacement DOF u (n) I of the UC can be approximated using a Component Mode Synthesis as follows:

u (n) I (ω) ≈ Φ s p (n) s (ω) + Φ c p (n) c (ω) (10)
where the static modes Φ s = Φ s| L Φ s| R (i.e. constrained interface modes) are calculated by:

Φ s = -K -1 II K IL K IR (11)
while the component modes are obtained selecting a truncated set of modal eigenvalues as shown in Eq.( 12). Note that the specification of a frequency limit criterion ω max is necessary:

                 K II -ν j M II φ j = 0 Φ c p c = j∈Ω φ j p c, j ; Ω = j, ν j ≤ ω 2 max ( 12
)
Krylov decomposition [START_REF] Stewart | A krylov-schur algorithm for large eigenproblems[END_REF] can be advantageously used to implement the Arnoldi algorithm required to solve Eq.( 12), given the problem's sparsity. Solutions are retained up to ν i ≥ ω 2 p , the maximal frequency above which modes are truncated. This frequency-dependent criterion is an important parameter impacting the overall size of the reduced model. The following (frequency-independent) projection matrix can be used to reduce the inner DOF of each UC's state vectors:

                 u (n) L u (n) I u (n+1) L                  (ω) =                  I 0 0 Φ s| L Φ c Φ s| R 0 0 I                                   u (n) L p c (n) u (n+1) L                  (ω) (13)
This approach was used effectively in several sub-structuring and homogenization techniques. However, the interfaces between the UCs are not affected by the CMS, although a variety of modal descriptions have been investigated

to reduce the cells' edges. In the context of a wave-based modelling of the structure, Floquet solutions (i.e. wave shapes computed by WFEM) can be used as an efficient way to describe the state vectors of a UC. The concept of Floquet state-vector reduction for dispersion analyses was introduced in [START_REF] Droz | A reduced formulation for the free-wave propagation analysis in composite structures[END_REF] then improved in [START_REF] Boukadia | A bloch wave reduction scheme for ultrafast band diagram and dynamic response computation in periodic structures[END_REF] by using the inverse WFEM formulation to speed up solution sub-set computation and include evanescent waves' contribution in the reduced model. The goal is to propose a frequency-independent projection matrix Ψ w to replace the edge state-vector

u (n)
L by a reduced one p (n) L . Interested reader can also refer to [START_REF] Beli | A projection-based model reduction strategy for the wave and vibration analysis of rotating periodic structures[END_REF] for other examples of MOR for the WFEM. For a sake of completeness, a state-vector reduction procedure (detailed in [START_REF] Boukadia | A bloch wave reduction scheme for ultrafast band diagram and dynamic response computation in periodic structures[END_REF]) is summarized below:

• The first step is to identify the partial contours of the Brillouin zone. It can be done using one direct (for λ(ω max )) and one inverse (for ω({|λ| = 1})) resolution of the spectral eigenproblem, using the UC element matrices,

• Then, a wavenumber sampling can be defined and used to solve the inverse (wavenumber-dependent) eigenproblems,

• The real and imaginary parts of each state-vector is collected to build a reduced solution subspace, using a SVD procedure,

• Eventually, one can build the wave projection basis Ψ w for each cell's edge u (n) L , so that the displacement statevectors can be written:

u (n) L (ω) = Ψ w p (n) L (ω) (14)
The local projection matrix Γ defined in Eq.( 8) for the overall UC is then written:

Γ =                  Ψ w 0 0 Φ s| L Φ c Φ s| R 0 0 Ψ w                  (15)
The global projection matrix can eventually be assembled using the local scale projection matrices as detailed in Eq.( 16). This local MOR approach is synthesized in Fig. 2.

U(t) =                                                            U (0) L U (0) I . . . U (n) L U (n) I U (n+1) L . . . U (N) L                                                            (t) =                                                           Ψ w 0 0 Φ s| L Φ c Φ s| R . . . Ψ w 0 Φ s| L Φ c Φ s| R 0 Ψ w . . . 0 Ψ w                                                                                                                      P (0) w P (0) c . . . P (n) w P (n) c P (n+1)
w . . .

P (N) w                                                            (t) = B P(t) (16)
The main advantages of Local MOR are its ease of computation (size of a single UC) and the fact that the resulting basis is intrinsically frequency-independent. However the drawback of Local MOR reduction arises when the number of local modes, waves or UCs is high: the dimension of the reduced model is the product of the number of local modes and waves by the number of UCs. Consequently, a large number of UCs may produce an excessively large reduction basis, while an insufficient number of UCs may be responsible for an inaccurate model. Moreover, Local MOR does not take into account the specific loads or boundary conditions hence generates a sub-optimal reduction basis.

Formulation of a global MOR based on the Floquet expansion of forced WFEM state-vectors

Note however that Local MOR can be used prior to the Floquet expansion described in Eq. [START_REF] Van De Walle | Virtual microphone sensing through vibro-acoustic modelling and kalman filtering[END_REF]. In this case, the modal basis can be over-determined without consequences on the final number of unknown wave amplitudes q + (ω) and q -(ω). The main advantage of WFEM remains, which is to provide fast dynamic response to a given load or BCs, disregarding the number of UCs composing the periodic structure. Its limitation lies in the frequency-dependence of the wave-modes derived from the spectral eigenproblem. This drawback will be tackled using a snapshot-type Proper Orthogonal Decomposition technique, to generate a low dimension model approximating the waveguide's dynamics from a subset of full state vectors' responses computed using WFEM. An illustration of the strategy is displayed in Fig. 3. The local projection matrix is first re-arranged from Eq.( 15), so that the state vectors write:

                 u (n) L u (n) I u (n+1) L                  (ω) =                  Ψ 0 0 Φ s| L Φ s| R Φ c 0 Ψ 0                                   p (n) L p (n+1) L p (n) I                  (ω) (17)
The dynamic stiffness matrices are then reduced using the above projection: D = Γ * DΓ. Note that neither the number of modes used to describe the inner cell's dynamic, nor the number of waves used to describe the cell's edges contribute to the final ROM dimension. Both selection constraints can be loosen with minor consequences to the overall performances of the MOR algorithm. A condensation of the inner DOF (i.e. component modes of the UC) can be achieved to define the dynamic stiffness matrix D:

D =           DLL DLR DRL DRR           -           DLI DRI           DII -1 DIL DIR (18)
The condensed form of the direct spectral problem defined in Eq.( 5) is solved to assemble the (frequencydependent) reduced wave solutions Ψ + (ω), Ψ -(ω) and the diagonal propagation constants' matrix Λ(ω). The forced response system is then built using the WFEM framework (i.e. Eq.( 7) for a free end waveguide). The solution is derived as a couple of complex wave amplitudes q + (ω) and q -(ω). The selection of appropriate snapshots to be later used for a computationally efficient proper orthogonal decomposition can be achieved using different heuristics (e.g. modal energy, error minimization). In this paper, the poles of the frequency response function averaged on each node subjected to external load are used to generate the snapshots. The key advantage of WFEM lies in its extremely low computational effort for punctual FRF determination, which allows a direct identification of resonances and antiresonances of the finite structure. Let us denote Ω = {ω i } 0≤i≤N f the frequency snapshots' subset collecting the poles of the loaded waveguide. Note that the global state vectors do not have to be computed to identify the poles (nor does the global dynamic stiffness matrix). The small number of snapshots N f allows a fast collection of a equal number of reduced edge state vectors p L using the Floquet expansion, defined for UC index 'n' at each frequency ω i by:

p (n) L (ω i ) = Ψ + (ω i ) Λ n (ω i ) q + (ω i ) + Ψ -(ω i ) Λ N-n (ω i ) q -(ω i ) (19)
while the inner part p I of the reduced state vectors can be obtained using Eq.( 18):

p (n) I (ω i ) = DII -1 DIL p (n) L (ω i ) DIR p (n+1) L (ω i ) (20)
The local displacement state vectors u (n) can then be retrieved using Eq.( 17). The organised global state vectors are collected into the expanded snapshot matrix U Ω as follows:

U Ω =                                                                     u (0) L u (0) I . . . u (N-1) I u (N) L                                  ω 1 ,                                  u (0) L u (0) I . . . u (N-1) I u (N) L                                  ω 2 , . . . ,                                  u (0) L u (0) I . . . u (N-1) I u (N) L                                  ω N f                                    (21)
The dimensions of matrix U Ω generally make it advantageous to use Lanczos bi-diagonalization with a partial re-orthogonalization [START_REF] Larsen | Lanczos bidiagonalization with partial reorthogonalization[END_REF] to achieve an computationally efficient singular value decomposition where the snapshot matrix writes:

U Ω = V Σ W H ( 22
)
where Σ is a diagonal matrix containing the sorted singular values of the decomposition. A truncation is achieved to a limited number 'r' of singular values using a residue minimum defined by: R i=1

σ i > (1 -) Tr (Σ) (23)
Note that in the general case, one has r ≤ N f , although the combination of an accurate pole identification (which produces reliable snapshots) and a low specified residue -→ 0 can yield r = N f , which becomes the final dimension of the reduced macro-scale model. The projection matrix is therefore assembled using the r columns-elements of V as:

B = V σ 1 , . . . , V σ r (24)
The projection matrix B should be recomputed if the boundary conditions are modified. The force boundary condition (described in Eq.( 7)) is the most general case, since other boundary conditions can be written similarly after achieving a dynamic condensation of the constrained edge element to the subsequent (n = 1) or second to last (n = N -1) cell. The use of explicit displacement-or force-based edge conditions can however prove advantageous when large UCs are considered. The global system defined in time domain by Eq.( 1) becomes after Floquet MOR:

               B * MB ∂ 2 p(t) ∂t 2 + B * CB ∂p(t) ∂t + B * KB p(t) = B * F (t) U(t) = B p(t) (25)
where p(t) is a reduced variable of dimension 'r'. The reduction parameters and numerical performances of the two methods are discussed in next section.

Discussion on numerical aspects and performances of Local and Floquet MOR

Comparison of modeling strategies and reduction parameters

Before analysing the method's accuracy for time-domain applications, let us discuss some sensitive numerical aspects of the proposed reduction scheme. Indeed, UC reduction is classically meant to accelerate the computation of dispersion curves, hence its performances are often solely discussed in terms of accuracy of the wave dispersion characteristics (e.g. wavenumbers) and compared for different UCs' meshes or geometries. These do not constitute a relevant benchmark for addressing finite structures. In that sense, Local MOR described in Sec.2.2 can be seen as a direct extension of existing UC reduction techniques for the construction of a global frequency-independent projection matrix, to be later used for time domain analysis. The local MOR and the FEM will both be compared to the Floquet MOR introduced in Sec.2.3, which also uses local UC reduction prior to the creation of a global projection basis, but now resulting from a Floquet expansion onto the entire structure. The finite nature of the structure means that Floquet MOR is affected differently by the UC reduction parameters than the Local MOR, as the number of UCs is not linked to the minimum number of required Floquet expanders. A reasonable hypothesis is that the local MOR will prove more efficient than Floquet MOR for structures constituted of a small number of large UCs, while Floquet MOR will become advantageous as the structure involves an increasingly large number of small UCs. To conduct a fair comparison between these MOR strategies and evaluate our method, let us consider the following parametric benchmarks where the two above-mentioned design parameters are distinguished in Sec.3.2 and 3.3: (i) the number of UCs and (ii) their scaling (or segmentation).

Influence of the model's dimensions on the computational performances

To evaluate the effects of the number of UCs on the Local and Floquet MOR techniques, one can apply these methods to the five structures displayed in Fig. solvers. This is due to the increasing time required to operate the projection matrices -which grow consistently with the number of UCs -due notably to the high number of state vectors required to reduce the UCs' interfaces.

To summarize, when the structure's size increases (i.e. the number of UCs, DOF, modes), the proposed Floquet MOR remains efficient without loss of benefits over FEM, while Local MOR becomes almost counter-productive as the length increases.

Influence of unit-cell's scaling on the computational performances

To evaluate the influence of the modelling scale on MOR efficiency, let us now consider the five models displayed in Fig. 6. Given the waveguide's homogeneity one can define an arbitrary UC dimension to divide the 1000 mm-long aluminium structure into identical segments. It is therefore divided into 5, 20, 100, 250 and 500 UCs of respective lengths 200 mm, 50 mm, 10 mm, 4 mm and 2 mm. Note that in all five cases the assembled FE model is strictly identical, namely a 34034 DOF model. The smallest constitutive UC is 1-element long (2 mm) and displayed in Fig. 6.

The dynamic analysis is conducted between 0 and 50 kHz with 1000 frequency samples and the load configuration, MOR residue parameters and CPU are the same as in Sec. while the resolution constitutes most of the computational effort for higher levels of the waveguide's segmentation.

Modelling To summarize, this parametric analysis shows that the Floquet MOR is always significantly faster than the alternatives for both reduction and resolution, reaching over 99% reduction compared with accelerated FEM (i.e. with SLU-solver). This conclusion stands, notwithstanding the segmentation of the structure. In fact, the Floquet MOR efficiency further increases with segmentation, in opposition to Local MOR which is more suited to lower levels of multiscality.

4. Application: transient response of a periodic structure

Description of the model

In this section the performances of both methods are studied for the computation of the transient dynamic response of a simple example of elastic metamaterial: the structure is a slender aluminium phononic-type waveguide with 30 elliptic inclusions producing both Bragg and locally-resonant bandgaps. A detailed description of the example's geometry, FE model and load configuration is provided in Fig. 8. The structure is subjected to three different force impulses applied along the vertical direction on the left boundary while the right end is free. The input signal is a Gaussian pulse F(t) of amplitude F 0 , center frequency ω 0 , cut-off frequency t c at truncating amplitude -100dB defined by:

F(t) = F 0 cos (ω 0 (t -t c )) exp -(t -t c ) 2 2σ 2 (26)
where σ = 1 2π f v and f v are respectively the Gaussian pulse's temporal and spectrum bandwidth variances. Displacements are then monitored at both ends over a duration of 2 ms.

Time integration and resolution

A direct time integration is used to solve the second order differential equation, based on Newmark method [START_REF] Newmark | A method of computation for structural dynamics[END_REF].

The numerical scheme used below is implicit and a constant acceleration scheme (i.e. unconditionally stable) is used by specifying α = 1 4 and δ = 1 2 in the integration constants defined as:

               a 0 = 1 α∆ 2 , a 1 = δ α∆ , a 2 = 1 α∆ , a 3 = 1 2α -1 a 4 = δ α -1, a 5 = ∆ δ 2α -1 , a 6 = ∆(1 -δ), a 7 = ∆δ ( 27
)
where ∆ is the selected time step.

Defining K eff the global effective stiffness matrix defined by: K eff = a 0 M + a 1 C + K and denoting U T , V T , A T the approximations of respective displacement state vectors' derivatives ∂ 2 U(t) ∂t 2 , ∂U(t) ∂t and U(t), the following iterative procedure can be used to solve the linear system subjected to the time-dependent load:

• The effective load vector is computed using:

R T +∆ = F (T + ∆) + M [a 0 U T + a 2 V T + a 3 A T ] + C [a 1 U T + a 4 V T + a 5 A T ] , (28)
• The displacement state vector at time T + ∆ is computed by:

U T +∆ = K-1 eff R T +∆ , ( 29
)
where Keff is a factorisation of K eff into lower -upper triangular matrices.

• The velocities and accelerations are derived using:

               A T +∆ = a 0 (U T +∆ -U T ) -a 2 V T -a 3 A T V T +∆ = V T + a 6 A T + a 7 A T +∆ , (30)
• The previous steps are resumed using the time step increment:

T = T + ∆
The time resolution is adapted to the different load cases described in the following, based on the highest frequency content of their normalized input spectrum above a 10 -4 threshold. 3 is computed using FEM and both MOR methods. A detailed description of each input signal's center frequency, fractional bandwidth ratio and time sampling is provided in the caption of Tab.3. Since the accuracy of the reduced model tends to decrease with frequency, it means that the model is expected to produce highly accurate results for both P1 and P2 and a slightly increased error for P3, whose spectrum extents closer to the upper validity limit of the ROM.

Results and accuracy

The computational efforts and DOF reduction produced by the two methods are detailed in Tab. (1) 0 =10 kHz and a fractional bandwidth ratio of 20%, in order to have a spectrum covering the most part of the first pass-band. Pulse P2 is centered on ω (2) 0 =31 kHz with a bandwidth ratio of 5% to be fully contained in the stop-band. Pulse P3 is centered on ω (3) 0 =47.5 kHz with a 2.5% bandwidth ratio, corresponding to the 2 nd pass-band. In each case, the input amplitude is F 0 =100 N and the simulation duration is 2 ms. The time resolution is set to guarantee a sampling frequency of 3× the center frequency of the signal, i.e. 30 kHz for P1 (∆ (1) =1.67µs), 93 kHz for P2 (∆ (2) =0.54µs) and 142.5 kHz for P3 (∆ (3) =0.35µs).

where • is the scalar product and the L 2 -norm. The errors produced by both methods on P1, P2 and P3 analyses are displayed in Fig. 10. In the P3 input case, an additional analysis was conducted, extending the duration up to 2 t max in order to observe the time convergence of the error.

Firstly, one can note the outstanding accuracy of Floquet MOR for P1 and P2 analyses. The error on P3 is slightly higher over time as expected, although it remains accurate for the first 1 ms and still converges around 3% at 4 ms.

Second, one clearly identifies regular peaks on all MAC plots. These times are corresponding to the reflection of the pulses on the edges of the structure. A possible explanation for this phenomenon is that each reflection on a boundary generates a variety of evanescent or highly decaying waves which the ROM fails to describe accurately. The absence of such error peaks for P2 (i.e. inside the bandgap) tends to support this interpretation. These discrepancies however vanish with the far field, probably due to the attenuation. One can thirdly observe that both Floquet MOR and Local MOR produce very similar errors. Interestingly, this can be interpreted as the consequence of their shared wave-mode decomposition prerequisite, inherent to these method (i.e. based on a frequency truncation) and possibly responsible for a poor projection of the waveguide's interfaces. The proposed Floquet MOR thus does not produce any additional error.

The time-dependent upper surface nodal displacements on the two locations (defined in Fig. 8) U Y (0, t) and U Y (L, t) are displayed for each input signal in Fig. 11. The reduced model produces a highly accurate response without any observable discrepancy over FEM, either on the source point and 30 UC further, at the end of the structure. The low frequency components, which are sometimes inadequately described by a reduced wave decomposition and here appear a translation motion (e.g. top-right in Fig. 11), are also perfectly approximated by the Floquet MOR method.

Noteworthy, the response to signal P2 (i.e. inside the stop-band) at x = L (middle-right in Fig. 11) is also accurately Convergence is shown to be attained with higher simulation duration for Pulse 3, and is established around -3%. described with the reduced model, although this is a residual displacement resulting from second-order and converted waves propagating through the 30 UCs. The fact that the Floquet MOR maintains high accuracy above the first Brillouin zone, both inside and outside bandgaps is remarkable and should be emphasized. Moreover, for the smallest signal's wavelength (P3, bottom in Fig. 11) which is the upper limit of the Floquet MOR validity bandwidth, the discrepancy with FEM remains invisible, despite the increased MAC error mentioned earlier. One can reasonably deduce that the method is still efficient and remains highly accurate up to the upper bandwidth used to generate the reduced model.

Conclusions

A Floquet MOR method was investigated in this paper, with the objective to create small and accurate ROM of periodic structures which can be used for time domain applications. A global frequency-independent projection matrix was built from a POD-like basis construction strategy and a Floquet expansion of local state-vectors. By leveraging the outstanding performances of the forced WFEM framework, the proposed Floquet MOR quickly generates an accurate ROM. Comparing an extension of a local reduction strategy (i.e. local MOR) to the Floquet MOR, the authors demonstrated a clear added value of the latter, which can be summarized as follows:

• Floquet MOR is considerably faster than local MOR, notwithstanding sparsity of the system or the solver used.

• The resulting reduced model is extremely small and nearly minimal given the number of modes in the bandwidth.

• Floquet MOR is less sensitive, if at all, to a change of the modeling scale or the number of unit-cells in the structure.

Still, some aspects of the method could be further investigated, as it was observed for example that the Floquet MOR is not strictly insensitive to the scaling. This means that a better trade-off could be sought between the snapshots sampling procedure and the SVD residue threshold. In time-domain applications, small spikes of the MAC error were also identified during interface reflections. Excluding the boundary cells from the reduction could therefore contribute to reduce near-field ROM discrepancies.

It is expected that Floquet MOR will be highly valuable in fields such as virtual sensing, real-time SHM and digital twins, by providing a unit-cell-based MOR technique, exploitable for the time-domain analysis of large-scaled finite periodic structures and metamaterials.

5 )

 5 where blackboard bold D = D b,b -D b,I D -1 II D I,b stands in the following for the condensed cell's matrix with sub-script notation b = {L, R}.

Figure 1 :

 1 Figure 1: Concept diagram of the forced WFEM.

FE

Figure 2 :

 2 Figure 2: Concept diagram for the Local MOR method. Local ROM is replicated for each UC in the structure and does not depend on the boundary conditions or load case.

FFigure 3 :

 3 Figure 3: Concept diagram for the Floquet MOR, exploiting Forced WFEM with a POD-like technique to generate global projection matrix from a Floquet expansion of the UC state-vectors.

Figure 4 :

 4 Figure 4: Description of the 5 FE models used for comparison between local and Floquet MOR techniques and detailed model of a single UC. The cell is made of aluminium, has a length d x =20 mm, a thickness h=10 mm and is meshed using 49 nodes and 72 linear 2D triangular plane-strain elements with 2 translations per node, resulting in a total of 98 DOF.

4 .

 4 The structures are built by concatenating 20, 30, 40, 50 and 100 identical UCs, resulting in respectively: 400 mm, 600 mm, 800 mm, 1000 mm and 2000 mm-long structures. The constitutive UC is displayed in Fig.4. The size of each model is displayed in the 3 rd row of Tab.1. For each of these structures, five different strategies are used to conduct the same dynamic analysis: (a) the FEM with a full solver; (b) the FEM with an accelerated Sparse LU/Cholesky solver (SLU) (to exploit matrices' sparsity); the Local MOR using full (c) and SLU (d) solvers; and finally (e) the proposed Floquet MOR 1 . The dynamic analysis is a frequency response computation of the beam subjected to a left-side vertical force between 0 and 200 kHz. Identical reduction parameters are used for the local reduction, waves and component modes are kept up to ω max = 400 kHz, the residues for the local UC and Floquet reductions are respectively set to ε UC = 10 -6 and Floq = 10 -9

Figure 5 :

 5 Figure 5: Reduction + 2000 FRF computations time vs. FEM full resolution. Number of DOF is proportionate to initial number of UCs.

3 . 2 .

 32 Evanescent waves (filtered to |λ| > 0.2) are collected up to 100 kHz and component modes up to ω max =150 kHz. Results are presented in Tab.2 and Fig.7.

Figure 6 :

 6 Figure 6: Description of a same structure modelled using 5 different segmentation. The cell has a length (1 element) d x =2 mm, a thickness (8 elements) h=10 mm and is meshed using 51 nodes and 16 quadratic 2D triangular plane-strain (X,Y) elements, resulting in a total of 102 DOF. The number of section elements in the respective UCs (from top to bottom) is: 1, 2, 5, 25 and 100, respectively.

Figure 7 :

 7 Figure 7: (a) Reduction + resolution time of local and Floquet MOR methods vs. FEM (right axis). (b) Bar diagram indicating the proportion of time dedicated to build reduced model vs. resolution (left axis).

Figure 8 :

 8 Figure 8: Application example of the 30-UC structure and detailed FE model of a single UC. The structure is made of aluminium (E=70 GPa, ν=0.33, ρ=2700 kg.m -3 ), has a length d x =20 mm, a thickness h=10 mm and is meshed using 1197 nodes and 537 quadratic 2D triangular planestrain (X,Y) elements, resulting in 2394 DOF per UC. Each periodic cell has a centered 16 mm per 4 mm elliptic inclusion. The overall model is made of 30 UCs and 70254 DOF.

4. 3 . 4 Figure 9 :

 349 Figure 9: (a) Frequency response function at input force location (co-localized), and at the end of the structure (x=L). Frequency spectra of pulses 1, 2 and 3 are overlaid. (b) Wave dispersion characteristics (real wavenumber parts), showing 4 propagating waves in the considered bandwidth, and corresponding bandwidths of the three pulses in the various pass-and stop-bands. Wave 1 is mainly actuated, although all waves (including evanescent) are produced by load.

Figure 10 :

 10 Figure 10: Time dependence estimation of the global error produced by Floquet and Local MOR methods using full field MAC over FEM results.

Figure 11 :

 11 Figure 11: Time signals comparison at two different locations for the three pulse's inputs.

  Original concept, Method, Software, Investigation, Writing-Original draft, Funding acquisition. Regis Boukadia: Method, Investigation, Writing-review & editing. Wim Desmet: Critical review, Project supervision

Table 1 :

 1 . Influence of the number of UCs (i.e. model size) on the performances of Local and

		Modelling/Solver		Number of UCs
	Method Step duration (ms)	20	30	40	50	100
		Model size (DOF)	1694 2534 3374 4214	8414		10 -1	ROM + FRF Computation
	FEM	Solving (full-ms/step)	197	562 1157 2295 15328
	Local MOR	Solving (SLU-ms/step) 7.92 Model size (DOF) 494 ROM building (s) 0.61 Solving (full-ms/step) 7.99 20.62 46.34 86.0 477.85 12.2 16.7 21.1 46.4 734 974 1214 2414 1.50 3.16 6.2 43.15 Solving (SLU-ms/step) 1.25 2.80 3.60 4.70 11.19	Total duration (vs full FEM)	10 -3 10 -2
	Floquet MOR	Model size (DOF) ROM building (s) Solving (full-ms/step)	92 2.05 0.50	134 2.20 0.79	179 2.60 2.70 220 1.32 1.86	353 4.24 4.15		10 -4	FEM (LU) Local ROM Local ROM (LU) Number of DOFs 2000 3000 4000 5000 6000 7000 8000 Floquet ROM

Floquet MOR using full and reduced solvers. For each method, final number of DOF is indicated as well as the duration needed to reduce and solve the system. Note that the SLU-solver is not used for the Floquet MOR, since it produces full matrices. All computations are made on the same 2.90 GHz CPU limited to 16 Go RAM. The resolution time is averaged on 2000 frequency samples and expressed in sec./step.

Table 2 :

 2 Detailed performances of each of the MOR algorithms and FEM simulation. The three

		/Solver		Number of divisions
	Method Step	5	20	100	250	500
		Model size (DOF)			34034
	FEM	Solving-SLU (s)			460.9
		Model size (DOF)	394	874 2423 6024 13026
	Local	UC reduction (s)	2.20	0.28	2.10	0.20	0.17
	MOR	Global projection	5.64	3.28	3.07	3.14	18.29
		Solving-SLU (s)	9.02 11.67 16.72 39.94 102.31
		Model size (DOF)	77	70	62	60	58
		UC reduction (s)	1.85	0.29	1.93	0.20	0.18
	Floquet
		Poles identification (s) 5.28	4.55	2.77	3.21	2.98
	MOR
		Floquet expansion (s)	1.56	1.26	1.16	1.24	1.42
		Projection (s)	0.31	0.28	0.25	0.25	0.24
		Solving (s)	0.40	0.35	0.26	0.25	0.26

methods (FEM with SLU resolution algorithm, Local MOR with SLU resolution algorithm, Floquet MOR) are implemented, resulting in the given model size (in DOF) while the ROM step duration is indicated for each method and UC scaling. Given the size of this model, only SLUsolvers could be used for FEM and local MOR comparisons, so as not to exceed computational and memory resources.

  [START_REF] De Klerk | General framework for dynamic substructuring: history, review and classification of techniques[END_REF]. The Floquet MOR produces a system of 160 DOF against the 70254 DOF of the initial FEM. Noteworthy, both reduction techniques are preferable to FEM in the three analyses, both in terms of computation time and model dimension.Interestingly, both Local and Floquet MOR take about the same time to create a reduced model. In fact, although the UC reduction is identical in both methods, the additional cost of Floquet MOR (i.e. for pole identification and Floquet expansion) is counterbalanced by its fast projection time on a much smaller basis. Overall, is it clear that the performances of Floquet MOR are far superior to FEM and to Local MOR. For the more expensive (i.e. Pulse 3) analysis, a resolution with Floquet MOR runs over 35× faster than Local MOR and 550× faster than FEM.

				Pulse 1 Pulse 2 Pulse 3 Model size
	Method	Steps and duration	1201	3721	5701	(DOF)
	FEM	Resolution (s) 137.41 425.89 656.71	70254
		Building ROM (s)		17.25
	Local MOR	Resolution (s)	9.05	27.35	42.31	2064
		Building ROM (s)		16.96
	Floquet MOR	Resolution (s)	0.42	1.02	1.15	160
	Accuracy is a crucial criterion to assess the performances of the reduced-order modelling techniques. Denoting
	U FEM (t) and U ROM (t) the global nodal displacement fields computed using FEM and Local or Floquet MOR respec-
	tively, the overall time-dependent reduction error is evaluated using the Modal Assurance Criterion residue defined
	by:
		ε(t) = 1 -	U FEM (t) • U ROM (t) U FEM (t) U ROM (t)		(31)

Table 3 :

 3 Computational effort of the different modelling techniques. For Local and Floquet reduction, the time dedicated to build the ROM is distinguished from the resolution time. Pulse P1 has a center-frequency ω

Although full solver is used as a benchmark, it is usually preferable to use SLU for FEM. However, this is not always true for Local MOR, which tends to produce full matrices.

Acknowledgements

This work was funded by the European Union's Horizon 2020 research and innovation programme, through the individual Marie Skłodowska-Curie fellowship of Christophe Droz (WIDEA) [grant number 797034]. The research of R. F. Boukadia is funded by an Early Stage Researcher grant within the European Project VIPER Marie Curie Initial Training Network [grant number 675441] and the KU Leuven Industrial Research Fund (IOF).


Johann Mortara
email: johann.mortara@univ-cotedazur.fr

Philippe Collet
email: philippe.collet@univ-cotedazur.fr

Anne-Marie Pinna-Dery
email: anne-marie.pinna@univ-cotedazur.fr

Extending the Identification of Object-Oriented Variability Implementations using Usage Relationships

Keywords: Software and its engineering → Software product lines, Software reverse engineering, Object oriented architectures variability identification, variability visualization, variability-rich object-oriented software systems, symfinder

Many variability-rich object-oriented systems rely on multiple traditional techniques (inheritance, patterns) to implement their variability in a single codebase. These variability implementation places are neither explicit nor documented, hampering their detection and variability comprehension. Based on the identification of symmetry property in seven implementation techniques, a first approach was proposed with symfinder to automatically identify and display the variability of a system in a graph-based visualization structured by inheritance. However, composition, or more generally the usage relationship, is extensively used to implement the variability in object-oriented systems, and without this information, comprehending the large amount of variability identified by symfinder is not trivial. In this paper, we present symfinder-2, an extension of the former approach that incorporates the usage relationships to better identify potential variability implementations. We provide two ways to mark classes as entry points, user-defined and automatic, so that the visualization is filtered and enables users to have a better focus when they identify variability. We also report on the evaluation of this extension to ten open-source Java-based systems.

INTRODUCTION

Whatever is their size, most modern object-oriented software systems are variability-rich [START_REF] Galster | Variability in Software Systems -A Systematic Literature Review[END_REF]. The variability among the software variants of these systems is often implemented in a single codebase using traditional techniques, such as inheritance, parameters, overloading, or software design patterns [START_REF] Gacek | Implementing Product Line Variabilities[END_REF][START_REF] Tërnava | On the Diversity of Capturing Variability at the Implementation Level[END_REF]. As they do not follow a fully-fledged software product line approach [START_REF] Capilla | Systems and software variability management[END_REF], their domain variability (i.e., features) is hardly documented or made explicit to code assets. This directly hampers the comprehension and maintainability of the implemented variability [START_REF] Pohl | Software Product Line Engineering: Foundations, Principles and Techniques[END_REF].

Several manual and automatic approaches can be found in the literature for locating, identifying, and visualizing domain features in code [START_REF] Roberto | A Systematic Mapping Study of Information Visualization for Software Product Line Engineering[END_REF][START_REF] Mistrík | Software Engineering for Variability Intensive Systems: Foundations and Applications[END_REF][START_REF] Rubin | A survey of feature location techniques[END_REF]. But, in these approaches, features are either required to be known in advance or to be identified through reengineering a set of clone-and-own or legacy software systems into a product line. Actually, considering the single codebase systems, there was a lack of approaches to identify variation points (vp-s) with variants 1 implemented by object-oriented techniques [START_REF] Lozano | An Overview of Techniques for Detecting Software Variability Concepts in Source Code[END_REF], but also a lack of solutions to visualize them [START_REF] Roberto | A Systematic Mapping Study of Information Visualization for Software Product Line Engineering[END_REF]. Complimentary to existing reverse engineering approaches, a new variability identification support was proposed with symfinder [START_REF] Tërnava | Identifying and visualizing variability in object-oriented variability-rich systems[END_REF], which handles variability-rich object-oriented systems developed in a single codebase. More recently, Michelon et al. [START_REF] Gabriela K Michelon | A Hybrid Feature Location Technique for Re-engineering Single Systems into Software Product Lines[END_REF] also proposed a hybrid technique using runtime trace generation on variants from the same codebase to locate features, showing the interest of the research community in the context of single codebase systems.

On its side, symfinder provides automatic identification and visualization of potential vp-s with variants in object-oriented code assets, with support for the Java and C ++ languages [START_REF] Mortara | Identifying and Mapping Implemented Variabilities in Java and C++ Systems using symfinder[END_REF][START_REF] Mortara | symfinder: a toolchain for the identification and visualization of object-oriented variability implementations[END_REF]. Its identification approach of potential vp-s with variants is using the property of symmetry in their implementation techniques, for example, inheritance defines a substitution symmetry between the immutable part of the superclass (the vp) and the possible changes in its subclasses (the variants). symfinder also provides a visual graph-based representation of the identified variabilities, in which vp-s and variants are nodes related by their inheritance links while visual parameters, such as size and colour intensity, vary according to the number of symmetries in the method level. The aim is to make zones with many symmetries easily distinguishable. symfinder has been successfully applied to several real Java and C ++ variability-rich systems [START_REF] Mortara | Identifying and Mapping Implemented Variabilities in Java and C++ Systems using symfinder[END_REF][START_REF] Mortara | Mapping Features to Automatically Identified Object-Oriented Variability Implementations: The Case of ArgoUML-SPL[END_REF][START_REF] Tërnava | Identifying and visualizing variability in object-oriented variability-rich systems[END_REF], showing that a large part of the identified potential vp-s and variants actually implemented domain features [START_REF] Mortara | Identifying and Mapping Implemented Variabilities in Java and C++ Systems using symfinder[END_REF], that is, they are actual vp-s with variants.

However, as more experiments and applications were conducted with symfinder, some drawbacks also appeared, especially on large systems. First, the visualization provided by symfinder exhibits vp-s and variants with heavy usage of inheritance, while some of the variability relies also on other mechanisms and are thus less visible. When some vp-s or variants are using others (e.g., through attributes or parameters), this usage relationship is not exploited while it should create denser zones of variability and help in better structuring the variability visualization of large systems. Second, the user feedback shows that one tends to start an analysis through some well-known entry point classes, being either exposed facade classes or a well-defined application programming interface (API) [START_REF] Martin P Robillard | Automated API property inference techniques[END_REF] that publish the customizable functionalities. However, symfinder did not provide any means to browse the identified variability from these entry points.

In the following, we first provide some background on the concept of symmetry used and the symfinder principles (Section 2). We then analyse the issues to be tackled (Section 3) and describe the different extensions made on the identification and visualization parts to build a new version of symfinder named symfinder-2 (Section 4). During the identification of potential variation points with variants, we take into account their usage relationships so to display them in the visualization. We also refine the visualization so that entry point classes of a targeted system, such as API classes, can be used to find more easily the important zones of variability. A parameterized density metric can also help in filtering the visualization when no entry points are given. We applied symfinder-2 to ten Java-based variability-rich systems and observed the impact of our changes and done improvements (Section 5). We notably evaluate the visualized graph, the remaining vp-s and variants, and the scalability of our extension. We finally discuss related work (Section 6), and conclude while evoking future work (Section 7).

BACKGROUND

In this section, we introduce the concept of symmetry in implementation techniques, vp-s with variants, their density, and symfinder.

Symmetry in object-oriented techniques

From the natural sciences, the symmetry of an object is defined as the immunity to a possible change [START_REF] Rosen | Symmetry Rules: How Science and Nature are Founded on Symmetry[END_REF]. Several studies have shown that most of the object-oriented techniques, such as inheritance, overloading, including software design patterns, can also be interpreted in terms of symmetry [START_REF] Coplien | Toward a General Formal Foundation of Design. Symmetry and Broken Symmetry[END_REF][START_REF] Zhao | Understanding symmetry in objectoriented languages[END_REF]. Considering a whole codebase, these techniques can be seen as local symmetries [START_REF] Tërnava | Identifying and visualizing variability in object-oriented variability-rich systems[END_REF], which allow a part of code to change while another part remains unchanged.

To illustrate these symmetries, let us consider an illustrative example from JFreeChart2 , a variability-rich object-oriented library that provides a family of charts, such as pie and meter charts, given in Figure 1. In this library, the abstract class Plot is the common part between the related variations of PiePlot and MeterPlot. These are some of the variations that make it possible to differentiate the variants of JFreeChart objects at the implementation level. The common and variation parts are commonly abstracted in terms of variation points (vp-s) and variants, respectively [START_REF] Czarnecki | Cool features and tough decisions: a comparison of variability modeling approaches[END_REF][START_REF] Jacobson | Software reuse: architecture process and organization for business success[END_REF]. In our example, the class Plot is a vp with two variants, PiePlot and MeterPlot. They can be abstracted as vp_Plot, v_PiePlot, and v_MeterPlot, respectively. Following the symmetry definition, inheritance defines a substitution symmetry for its subtypes. Here, the possibility of a change in the superclass Plot materializes in its potential different subtypes, such as PiePlot and MeterPlot, which vary regarding the way how they draw a chart. Still, they also preserve and conform to the common behaviour of their superclass.

Considering this way of reasoning, which is also illustrated in Figure 1, most of the object-oriented variability implementation techniques, including design patterns, can be described in [START_REF] Coplien | Toward a General Formal Foundation of Design. Symmetry and Broken Symmetry[END_REF][START_REF] Zhao | Understanding symmetry in objectoriented languages[END_REF]. At an abstract level, a vp represents the unchanged part while its variants are the changed parts in code assets of a system [START_REF] Tërnava | Identifying and visualizing variability in object-oriented variability-rich systems[END_REF].

The symfinder approach

In object-oriented systems, the domain variability is hardly documented and traced to code assets (e.g., in terms of features in a feature model [START_REF] Kyo | FORM: A feature-; oriented reuse method with domain-; specific reference architectures[END_REF]), as they are mainly not organized as fullyfledged product lines [START_REF] Apel | Feature-Oriented Software Product Lines[END_REF]. Consequently, their existing vp-s with variants in code assets, as in Figure 1, are neither explicit nor documented. In the literature, there are several manual and automatic approaches for locating, identifying, and visualizing domain features in code [START_REF] Dit | Feature Location in Source Code: A Taxonomy and Survey[END_REF][START_REF] Roberto | A Systematic Mapping Study of Information Visualization for Software Product Line Engineering[END_REF][START_REF] Gabriela K Michelon | A Hybrid Feature Location Technique for Re-engineering Single Systems into Software Product Lines[END_REF][START_REF] Mistrík | Software Engineering for Variability Intensive Systems: Foundations and Applications[END_REF]. In these approaches, features are either required to be known in advance or to be identified through reengineering a set of clone-and-own or legacy software systems into a product line. Unlike them, our recent tooled approach named symfinder [START_REF] Mortara | Identifying and Mapping Implemented Variabilities in Java and C++ Systems using symfinder[END_REF][START_REF] Mortara | symfinder: a toolchain for the identification and visualization of object-oriented variability implementations[END_REF][START_REF] Tërnava | Identifying and visualizing variability in object-oriented variability-rich systems[END_REF] focuses on identifying the variability implementations in single codebase object-oriented systems. It provides automatic identification and visualization of potential vp-s with variants in code assets of a Java or C ++ -based variability-rich system 3 . The identification approach of potential vp-s with variants is based on the identification of local symmetries in seven objectoriented variability implementation techniques, namely class as type, class subtyping, method and constructor overloading, strategy, template, decorator, and factory patterns.

symfinder provides a visual representation of the potential vp-s and variants. Figure 2 shows the visualization of identified potential vp-s with variants for the illustrative example of JFreeChart given in Figure 1. Each class level vp or variant is visualized as a circle node that points out the used implementation technique, edges show their inheritance relationship, then the size and shades of the red colour of nodes indicate the number of method level vp-s with variants 4 . We follow here the visual principles of preattentive perception [START_REF] Christopher G Healey | High-speed visual estimation using preattentive processing[END_REF], using some of the seven parameters that can vary in visualization in order to represent data. The visualized variability then forms a disconnected graph based on inheritance links, while the visual representation of the metrics associated with each node (e.g., size, colour intensity) creates zones that can be easily distinguished by their different density of symmetries 5 . For example, the left part of Figure 2 is denser than the right part 6 .

The reader can note that we use vp-s and variants at the conceptual level while using sometimes classes when we discuss the technical identification parts, as well as nodes as the representation of classes in the visualization. In addition, symfinder provides several other options for facilitating the comprehension of the variability of a system: filtering out the solitary nodes, showing or hiding class variants, metrics on identified vp-s and variants at different levels, coloring vp-s with variants that belong to the same package, and information on vp or variant code assets when hovering 7 .

PROBLEM STATEMENT

symfinder has been applied to several Java and C ++ variability-rich systems [START_REF] Mortara | Identifying and Mapping Implemented Variabilities in Java and C++ Systems using symfinder[END_REF][START_REF] Tërnava | Identifying and visualizing variability in object-oriented variability-rich systems[END_REF], showing that it can successfully identify the potential variability in their code assets. It also helped to comprehend this variability by using the generated visualization. Moreover, an application on the ArgoUML-SPL [START_REF] Martinez | Feature location benchmark with argoUML SPL[END_REF] has demonstrated that a large part of the identified vp-s and variants by symfinder actually implement the ArgoUML's reverse-engineered domain features [START_REF] Mortara | Identifying and Mapping Implemented Variabilities in Java and C++ Systems using symfinder[END_REF]. However, as more experiments and applications were conducted with symfinder, two issues also appeared.

Issue 1. Identifying inheritance relationships is not enough. By construction, the identified vp-s with variants with heavy usage of inheritance are the most visible in the symfinder visualization, that is, those implemented by classes that are represented as nodes. They are simply grouped with this relationship and the other visualized indicators, such as the size and colour intensity of the nodes. While the experiments show that real vp-s and variants were successfully identified with symfinder, it also appeared that their visualized density based only on the inheritance relationship is not enough to comprehend the variability of a system. In all studied systems, we identified many cases in which subclasses being variants were using other vp-s and variants. For example, in JFreeChart, CompassPlot, a variant of the vp Plot, uses every variant of the vp MeterNeedle, but they are displayed as disconnected in the visualization.

Besides, as studied by Zhao and Coplien [START_REF] Zhao | Understanding symmetry in objectoriented languages[END_REF], reusability in object-oriented systems is about the instantiation of templates, composition of the instances, and substitution of the instances. symfinder clearly identifies substitution of the instances, as for instance inheritance has a substitution symmetry [START_REF] Tërnava | Identifying and visualizing variability in object-oriented variability-rich systems[END_REF], but the composition of other instances is not taken into account while it is used for reuse of code assets and implementation of design patterns. In the following, we will consider this relationship in the broad sense of the term, as a usage relationship, encompassing cases where a class uses another class in attributes and parameters of its methods.

Issue 2. Entry points are missing. With its zooming, filtering, and hovering capabilities, symfinder's visualization naturally relies on the Shneiderman visual information seeking mantra [START_REF] Shneiderman | The eyes have it: A task by data type taxonomy for information visualizations[END_REF]: overview first, zoom and filter, then details on demand. Nevertheless, in many systems and especially in large ones, users gave us feedback that they were missing some clear entry points to start browsing the visualization. The vast majority of studied systems indeed exposed facade classes, which were natural entry points expected by the users, or even a well-defined application programming interface (API) [START_REF] Martin P Robillard | Automated API property inference techniques[END_REF] where their reusable and customizable functionalities are textually exposed. The second issue we thus identify in this paper is the need for entry points to be specified and exploited in the visualization to facilitate variability comprehension.

symfinder-2

To address the identified issues, we extended the variability identification and visualization of symfinder, creating symfinder-2. Its sources and conducted experiments are publicly available 8 9 .

Handling the usage relationships

To address Issue 1, we extended the identification step by also identifying the usage relationships among the identified vp-s with variants in a codebase. At parsing time, each usage relationship is identified and added to the symfinder-2 graph database. Specifically, we consider that class A uses another class B if class B is used as a field type or method parameter type in class A [START_REF] Freeman | Head First Design Patterns[END_REF]. For instance, in the example of JFreeChart in Figure 1, class Plot is used by class JFreeChart as it is referenced by a field inside the Plot class. Therefore, each of these classes will be identified by symfinder-2 as either vp or variant and will be visualized with their inheritance and usage relationships.

Visually, we show in symfinder-2 both inheritance and usage relationships, where inheritance relationships are grey arrows and usage relationships are represented as dashed (light) blue arrows. For example, the new visualization for the JFreeChart example, given in Figure 1, is shown in Figure 3. In comparison with its prior visualization in symfinder, given in Figure 2, now both relationships are represented. In particular, the usage relationship between Plot and JFreeChart is now explicit in Figure 3.

Handling the entry points

To address Issue 2, we extended the visualization to provide some entry points from which users can begin the exploration of the identified variability in a targeted system. In the following, we present two kinds of entry points that are available in symfinder-2.

User-defined entry points. We first added the possibility to filter and visualize only those potential vp-s with variants of a system that are under analysis by a user. As shown in Figure 3, one can choose a specific classpath and use it to refine the identified variability in its visualization. symfinder-2 will show only those potential vp-s with variants and their interconnected vp-s or variants that are under the inheritance and usage relationships. For example, from the shown data in Figure 3, 926 potential vp-s at class and method levels are identified in JFreeChart.

Secondly, we added the possibility to filter in the potential vp-s with variants based on the direction or level of their usage relationships (cf. 'usage-type' and 'usage-level' boxes as in Figure 3). For example, the usage type OUT filters in only those class level vp-s and variants that use other vp-s or variants. Then, a usage level of 4 will display all vp-s and variants that have a usage relationship to one of the vp-s or variants through at most 4 identified usages transitively.

Lastly, we added the possibility to filter out the less dense zones with variability in the visualization. To achieve this filtering, a user needs to set a density threshold of potential vp-s with variants, which density is based on their usage and inheritance relationships 10 . For instance, the visualization from JFreeChart given in Figure 2 shows two places with different densities of potential vp-s and variants. The left one seems denser than the right one with the inheritance-only symfinder, whereas when usage relationships are considered using symfinder-2, both of these places are interrelated and create a new denser zone of potential vp-s and variants given in Figure 3. The density threshold is configured through setting two parameters in symfinder-2's configuration file: (1) the minimum number of variants that a potential vp has, and (2) the maximum number of usage relationships between potential vp-s. These two values are used in queries over the graph database to identify the local symmetries that are under or above the given threshold. Finally, symfinder-2 generates a new visualization where the potential vp-s and variants that are under the threshold are coloured in light grey and can be further excluded ('Show only dense zones' in Figure 3). As an example, in JFreeChart we set a threshold with the number of variants in a potential vp ≥ 5 and the usage relationships between the potential vp-s or variants ≤ 3. Hence, the v_DefaultDrawingSupplier variant is highlighted in grey as its vp_DrawingSupplier vp has less than 5 variants.

Automatically defined entry points. Another kind of entry points that we implemented in symfinder-2 is the automatic filtering of potential vp-s and variants that are advertised as being part of an API. Actually, symfinder-2 identifies classes annotated by the API Guardian library 11 as specific API entry points. The visualization then automatically adds all identified vp-s and variants as entry points so that they can be used to refine it as described in the previous paragraph. These added capabilities are made available to help users to chose interesting entry points in the visualization and to comprehend progressively the whole identified variability of a given system.

EVALUATION

In this section, we first introduce the subject systems and then define the research questions to evaluate the extended approach of symfinder-2. The symfinder-2 and the used data to obtain the presented results are available online in a reproduction package 12 .

Subject systems

To evaluate symfinder-2, we chose ten popular variability-rich subject systems, being Java applications, frameworks, or libraries (cf. Table 1). For the time frame of up to twelve last years, they have received between 150 and 8, 000 stars in GitHub, but we particularly considered them because of the following criteria. The first six ones were already used to evaluate the first version of symfinder [START_REF] Tërnava | Identifying and visualizing variability in object-oriented variability-rich systems[END_REF], namely Java AWT, Apache CXF, JUnit, Apache Maven, JFreeChart, and ArgoUML. Then, we chose the other three as they use in their codebase a form of API annotations, namely the API Guardian library, to annotate each code unit that constitutes their API. These new chosen systems are Cucumber -a framework for BDD testing, Logbook -a library to enable logging for different client-and server-side technologies, and Riptide -a library based on Spring to implement client-side response routing. Finally, we selected the NetBeans IDE because of its size with about 5 𝑀 lines of code (LoC), which helps in evaluating the scalability issues of both the approach and the prototyped toolchain.

𝑅𝑄 1 : Improved visualization

The visualization generated by symfinder shows the class level vp-s with variants as nodes, which are linked together through inheritance relationships. Whereas with symfinder-2, they are linked together through inheritance and usage relationships. In both cases, we observed that they form disconnected graph structures, which looks slightly different. Therefore, we defined 𝑅𝑄 1 :does the identification of usage relationships have changed the variability visualization of a given system by symfinder-2? To this end, we applied symfinder and symfinder-2 to all ten subject systems and compared their respective visualizations. Results are given in Table 2. On all studied systems, we notice a smaller number of disconnected graphs and isolated nodes by symfinder-2 for the same number of nodes displayed, meaning that zones in the visualization 12 https://zenodo.org/record/4946730 that seemed previously uncorrelated are now linked through usage relationships and appear as such on the new visualization. We observe that the difference between the number of disconnected graphs is not proportional to the size of the studied system. For instance, the number of NetBeans' graphs are reduced by 61% whereas JFreeChart's graphs are reduced by 90%. However, their number could be related to the architecture of the project. A project of an important size may have an architecture in layers, limiting the number of interactions between classes, and therefore exhibit fewer usage relationships. Besides, we notice that some isolated nodes still appear on symfinder-2's visualization. This may suggest that other usage mechanisms are present in the studied systems [START_REF] Lau | A taxonomy of software composition mechanisms[END_REF].

Taking into account these specific types of usage relationships is part of our future work. Further, it is important to emphasize that although the visualizations by symfinder and symfinder-2 have different numbers of disconnected graphs, their overall number and kinds of identified vp-s with variants remain unchanged. This indicates that symfinder-2 is an extension of symfinder with an intact variability identification. To conclude, the reduced number of disconnected graphs by symfinder-2 shows an improved and denser visualization of the identified vp-s and variants for a given system.

𝑅𝑄 2 : Starting density threshold

In addition to the two basic user-defined variability filtering capabilities that are added in the visualization and can be activated interactively, the third filtering capability is based on the density of potential vp-s and variants. As explained in Section 4.2, our filtering by density has for objective to reduce the number of potential vp-s and variants (a.k.a., nodes) visible on the visualization through identifying those vp-s that have a minimum number of variants or those that are linked to another vp through a maximum number of usage relationships. This threshold is set before their identification, and setting it is not trivial as the user needs to know which is the right threshold to start with. Therefore, we defined 𝑅𝑄 2 : what is the starting density threshold to begin with the comprehension of the visualized variability by symfinder-2?

To this end, we run symfinder-2 on each subject system with three different density thresholds, which is set by two parameters: the minimum number of variants of a potential vp and the maximum number of usage relationships between potential vp-s and variants (i.e., usage hops). Namely, the thresholds with (A) ≥ 5 variants and ≤ 3 usage hops, (B) ≥ 10 variants and ≤ 3 usage hops, and (C) ≥ 30 variants and ≤ 2 usage hops. We have carefully chosen these parameters values, based on a previous empirical evaluation with the ArgoUML, JFreeChart, and Java AWT, and for which we have the best knowledge to manually evaluate the impact of the density threshold. By increasing the threshold on the number of variants, we aim to consider only highly-dense vp-s , whereas by decreasing the threshold on the usage hops between such vp-s, we aim to consider only highly-dense vp-s which are close in terms of usage relationships. The obtained results are given in Table 3. It can be observed that, in all subject systems, fewer nodes are displayed when using any of the three threshold values. For instance, JFreeChart has in total 578 identified vp-s with variants at class level. After applying the three threshold values, 34, 15, and 3 vp-s and variants (i.e., entry points) remain, respectively. Moreover, the number of remaining vp-s decreases as we increase the minimum number of variants for a vp, and decreases the number of usage relationships between both of them. However, these values are not adapted to such projects like Logbook or Riptide, for which no vp-s remain with these three default thresholds.

These results suggest that determining a set of appropriate values for the parameters when setting the density threshold is highly dependent on the studied project's characteristics. That is, even some large projects in terms of lines of code, such as Maven, can have a considerable number of potential vp-s but with few variants. For such reason, setting a high density threshold may filter out most or all potential vp-s with variants (i.e., there will be no entry points). Based on our experiments with ten subjects, we conclude that the first threshold value (i.e., ≥ 5 variants and ≤ 3 usage hops) can be used as a good starting point to begin the exploration and comprehension of the identified variability by symfinder-2.

𝑅𝑄 3 : Usefulness of API-based filtering

As is shown in Table 3, the number of identified potential vp-s with variants is extensive (cf. nodes). For example, the Cucumber framework has 238 potential vp-s with 282 variants at class and method levels. Analysing and comprehending the Cucumber's variability from the provided visualization 13 , which exhibits 331 nodes, can 13 Cucumber's visualization: https://deathstar3.github.io/symfinder2-demo/ visualizations/cucumber-v6.8.0-usage.html be really difficult. But, as a testing framework, Cucumber has an API that exposes classes for defining the dependency steps and an object factory for customizing the dependency injections. Therefore, we defined 𝑅𝑄 3 : is the API information of a given system useful to simplify its identified variability by symfinder-2?

To answer it, we ran symfinder-2 on three subject systems, namely Cucumber, Logbook, and Riptide, while taking into account the code units annotated by developers using the API Guardian library.

We then compared the number of nodes that are displayed in their visualizations before and after using their respective API to filter in the related classes. Results are given in Figure 4. It can be observed that in all three subjects the visualization with the applied API has notably fewer nodes than the original visualization by symfinder-2. We manually checked that, while reducing the number of potential vp-s with variants, it always shows those that are considered as the most API relevant. They can help to comprehend each system's variability, that is, to give us an insight on the implemented domain variability. In the three cases, these vp-s with variants can be used as entry points for users in order to begin with the variability comprehension of the systems. We interpret this filtering by an API as facilitation in the overview and zooming parts of the Shneiderman mantra [START_REF] Shneiderman | The eyes have it: A task by data type taxonomy for information visualizations[END_REF]: overview first, zoom and filter, then details on demand. In the longer term, we believe that this should be extended by the integration of the symfinder-2 toolchain with other sources that contain variability information for a given system, or simply to show how different variability information sources could be blended together.

𝑅𝑄 4 : Scalability

Lastly, we defined 𝑅𝑄 4 : does the identification of usage relationships impact the scalability of symfinder-2? This question aims at determining whether the additional functionalities in symfinder-2 harm its scalability. The impact could be located in the identification phase, as the usage relationship is identified in the source code representation, as well as in the visualization phase, where new elements are computed to filter all displayed elements. To this end, we measured the time taken to identify and visualize the variability in all ten subject systems. We conducted our experiments on a Linux environment running Arch Linux 5.11.12-arch1-1 x64 with Intel i7-9850H (12 cores) @ 4.6GHz and 32Go memory. The visualizations are tested using Mozilla Firefox 87.0 and Google Chrome 89.0.4389.114. We noticed that the computation needed to render and display the visualization is not very time consuming, with 700 ms on Chrome and 850 ms on Firefox for the NetBeans system. The identification step is clearly the most time-consuming activity. We hence measured and compared this time on the ten subject systems with symfinder and symfinder-2. Figure 5 summarizes the obtained execution times for both versions (with density calculation based on thresholds of 5 variants and 3 hops). Although the execution time with symfinder-2 is higher than with symfinder for every system, we observe that the difference increases with the size of a system and number of identified vp-s and variants, for instance, 20% of difference for Riptide (23 sec → 28 sec) and 85% of difference for NetBeans (01:02:10 → 01:55:04). This can be explained by the fact that a higher number of relationships between classes needs to be parsed and treated by symfinder-2 in its database. While there seems to have an exponential evolution w.r.t. LoC with systems of the size of NetBeans (5 M LoC), we believe the analysis step is still adapted to large systems, as symfinder was also successfully applied to Firefox and its 25 M LoC [START_REF] Mortara | symfinder: a toolchain for the identification and visualization of object-oriented variability implementations[END_REF]. Then, waiting for around 1 or 2 hours to run symfinder-2 only on the new releases of a project, for example, every 6 months, is affordable.

Discussion and threats to validity

Summary of 𝑅𝑄 1 -𝑅𝑄 4 . symfinder-2 provides a more focused identification and visualization of relationships among the potential vp-s with variants than symfinder. Depending on the system's size, symfinder-2 can take between 30 seconds to 2 hours to identify between 250 and 11 K potential vp-s with variants. Besides, it supports users with up to four ways to begin the variability comprehension of a given system from its visualization. In particular, our experiments suggest using the density threshold with ≥ 5 variants and ≤ 3 usage hops or, if available, the API-based filtering.

Internal threats to validity. A first internal threat concerns the distinction of real vp-s and variants from the potential ones proposed by the different versions of symfinder. For the considered subject systems (except NetBeans that was too large), we manually determined whether the remained vp-s and variants after applying the thresholds represent some variability implementations. We thus did a sample verification by examining identified classes, checking for their documentation on the project website, and devising the kind of variability that was implemented. We could be partially wrong in our interpretation, but we believe it has a limited impact. Then, determining whether an identified vp or variant actually implements some domain variability was hard to conclude as non of the subject systems, except ArgoUML, had a ground truth.

External threats to validity. To address the research questions, we used up to 10 subject systems, which vary across domains, size, type, and developers. While the dataset is still small, we have good confidence that the obtained results also apply to other Java-based variability-rich systems of mid-size. Besides, our experiments over NetBeans show that while the toolchain is likely to scale on very large systems, the proposed improvements in symfinder-2 are not sufficient to comprehend all the implemented variabilities. Entry points and the usage links enable to provide a better overview and better filtering over the system, but it is still difficult to browse effectively towards a comprehension of the system variability. We expect future work on exploiting additional information (e.g., the architecture of the system, component definitions) and on providing a further improved visualization to facilitate this activity.

RELATED WORK

Multiple approaches for feature location [START_REF] Roberto | A Systematic Mapping Study of Information Visualization for Software Product Line Engineering[END_REF][START_REF] Mistrík | Software Engineering for Variability Intensive Systems: Foundations and Applications[END_REF][START_REF] Rubin | A survey of feature location techniques[END_REF] and feature identification have been proposed [START_REF] Wesley Kg Assunção | Reengineering legacy applications into software product lines: a systematic mapping[END_REF], but they mainly rely on multiple systems [START_REF] Martinez | Bottom-up technologies for reuse: automated extractive adoption of software product lines[END_REF][START_REF] Ziadi | Feature identification from the source code of product variants[END_REF] or use a form of feature annotations in their codebase [START_REF] Hunsen | Preprocessor-based variability in opensource and industrial software systems: An empirical study[END_REF][START_REF] Duc | Validating consistency between a feature model and its implementation[END_REF][START_REF] Liebig | An analysis of the variability in forty preprocessor-based software product lines[END_REF]. Unlike them, we use the property of symmetry in object-oriented constructs and design patterns, to identify the implemented variability in a single codebase. Outside the software variability domain, some other works rely on the identification of inheritance and composition to define hotspots, being zones of an object-oriented design that are exposed to client software and hence have to be comprehended to ensure reuse [START_REF] Flores | HotSpotter: a JavaML-based approach to discover Framework's HotSpots[END_REF][START_REF] Schauer | Hot spot recovery in object-oriented software with inheritance and composition template methods[END_REF]. While this identification relies on design patterns detection used to implement reusable interfaces, it was not related to variability implementations.

Previous studies on software API comprehension focus on the extraction of usage patterns relying on unit tests analysis [START_REF] Zhu | Mining api usage examples from test code[END_REF], client code analysis [START_REF] Zhong | MAPO: Mining and recommending API usage patterns[END_REF], or approaches combining client and library code analysis [START_REF] Saied | A cooperative approach for combining clientbased and library-based API usage pattern mining[END_REF], as well as their evolution in time [START_REF] Huppe | Mining complex temporal API usage patterns: an evolutionary approach[END_REF][START_REF] Aymen Saied | Towards assisting developers in API usage by automated recovery of complex temporal patterns[END_REF] to help developers in library reuse. In the variability domain, APIs have been studied for the variability of their evolution [START_REF] Alrubaye | Variability in Library Evolution: An Exploratory Study on Open-Source Java Libraries[END_REF], but not to facilitate the comprehension of variability implementations through a visualization approach as in our work.

CONCLUSION

As object-oriented variability-rich software systems are often not organized as software product lines, identifying their variability implementations is difficult. While proposing a first solution, the identification method provided by symfinder is incomplete as it does not take into account usage relationships between classes and cannot provide appropriate means to help users start the comprehension activity. In symfinder-2 we extended this method with usage relationships and provided filtering capabilities based on both automatically detected API-related entry points, and a userdefined density threshold. Application to ten systems has shown that symfinder-2 provides a more focused identification and visualization of object-oriented variability implementations.

In the future, we plan to include in our tooled approach additional object-oriented metrics and usage mechanisms [START_REF] Lau | A taxonomy of software composition mechanisms[END_REF] and to exploit architectural information of the systems to improve the comprehension of variability implementations and their evolution over time. We also intend to study other types of visualization that would provide a more intuitive comprehension of the system's implemented variability.

Figure 1 :

 1 Figure 1: Inheritance as a local symmetry in JFreeChart

-

  Class as type (vp or variant), -Class variant with inner vp-s -Abstract class (vp), -Interface (vp), -Constructor overloads (vp) -Method overloads (vp), -→ Inheritance relationship Node with 'F', 'S', 'T', 'D' symbol -Factory, Strategy, Template, or Decorator pattern

Figure 2 :

 2 Figure 2: Identified vp-s and variants using symfinder for the excerpt of JFreeChart given in Figure 1.

Figure 3 :

 3 Figure3: The identified vp-s and variants using symfinder-2 for the excerpt of JFreeChart given in Figure1

Figure 4 :

 4 Figure 4: Number of vp-s and variants (nodes) displayed on the visualization before and after refinement by API

Figure 5 :

 5 Figure 5: Execution time in symfinder and symfinder-2

Table 1 :

 1 The ten variability-rich subject systems. API: Ddocumented or A -annotated. Type: L -library, F -framework, and A -application

	Subject system	Commit	LoC	# vp-s	# variants	API / Type
	Java AWT	3319fcb	69,974	795	1,706 D / L
	Apache CXF	4da7b71	48,655	3,403	7,625 D / F
	JUnit	60aaf96	7,717	109	245 D / F
	Maven	97c98ec	105,342	612	1,147 D / A
	JFreeChart	1f6a91f	94,384	926	1,923 D / L
	ArgoUML	d135342	134,367	776	1,959 D / A
	Cucumber	323f724	42,662	238	282 A / F
	Logbook	f0f36e7	16,210	96	162 A / L
	Riptide	48b03a7	12,626	102	218 A / L
	NetBeans	cade258	5,058,448	3,621	6,736 D / A

Table 2 :

 2 Comparison of the number of disconnected graphs and isolated nodes with symfinder and symfinder-2

	Subject	Nodes	symfinder Graphs Isolated Graphs Isolated symfinder-2
	Java AWT	431	55	142	2	20
	Apache CXF	3085	473	1149	105	500
	JUnit	118	23	36	6	18
	Maven	616	177	172	21	79
	JFreeChart	578	54	167	5	51
	ArgoUML	1270	123	460	38	183
	Cucumber	331	45	122	14	50
	Logbook	117	19	40	4	16
	Riptide	89	20	37	8	19
	NetBeans	3498	504	1666	195	836

Table 3 :

 3 #nodes identified as being part of dense zones compared to the total number of nodes in all subject systems

				symfinder-2
	Project	symfinder	≥ 5 v-s	≥ 10 v-s	≥ 30 v-s
		nodes ≤ 3 hops ≤ 3 hops ≤ 2 hops
	Java AWT	431	28	22	3
	Apache CXF	3086	98	32	4
	JUnit	118	5	0	0
	Maven	616	8	1	0
	JFreeChart	578	34	15	3
	ArgoUML	1258	40	15	3
	Cucumber	331	4	0	0
	Logbook	117	0	0	0
	Riptide	89	0	0	0
	NetBeans	3494	58	22	2

http://www.jfree.org/jfreechart/

In this paper, we extend only the support for Java-based systems.

PiePlot is a variant of Plot, however it is also a vp as it has two variants, thus is described as vp_PiePlot on the figure.

It should be noted that the identification is based on symmetry in implementation techniques while the visualization on their density. Hence, one should not confuse the two and expect any kind of symmetry in the visualization.

The whole visualization for the JFreeChart 1.5.0 system is available online: https: //deathstar3.github.io/symfinder2-demo/visualizations/jfreechart-v1.5.0.html

Blue names and arrows in Figure2have been manually added for illustration.

symfinder-2's website: https://deathstar3.github.io/symfinder2-demo/

https://doi.org/10.5281/zenodo.4946729

Due to page limitations, a detailed definition and implementation of the density threshold is given here: https://deathstar3.github.io/symfinder2-demo/density.html

https://github.com/apiguardian-team/apiguardian

ACKNOWLEDGMENTS

We thank Florian Ainadou, Paul-Marie Djekinnou, and Djotiham Nabagou for their contribution in the development of symfinder-2.