Empowering Advanced Parametric Modes Clustering from Topological Data Analysis

Modal analysis is widely used for addressing NVH—Noise, Vibration, and Hardness—in automotive engineering. The so-called principal modes constitute an orthogonal basis, obtained from the eigenvectors related to the dynamical problem. When this basis is used for expressing the displacement field of a dynamical problem, the model equations become uncoupled. Moreover, a reduced basis can be defined according to the eigenvalues magnitude, leading to an uncoupled reduced model, especially appealing when solving large dynamical systems. However, engineering looks for optimal designs and therefore it focuses on parametric designs needing the efficient solution of parametric dynamical models. Solving parametrized eigenproblems remains a tricky issue, and, therefore, nonintrusive approaches are privileged. In that framework, a reduced basis consisting of the most significant eigenmodes is retained for each choice of the model parameters under consideration. Then, one is tempted to create a parametric reduced basis, by simply expressing the reduced basis parametrically by using an appropriate regression technique. However, an issue remains that limits the direct application of the just referred approach, the one related to the basis ordering. In order to order the modes before interpolating them, different techniques were proposed in the past, being the Modal Assurance Criterion—MAC—one of the most widely used. In the present paper, we proposed an alternative technique that, instead of operating at the eigenmodes level, classify the modes with respect to the deformed structure shapes that the eigenmodes induce, by invoking the so-called Topological Data Analysis—TDA—that ensures the invariance properties that topology ensure.


Introduction
Linear structural solid dynamics [1] expressed in the time domain results in the linear system of second order ordinary differential equations with the mass, damping, and stiffness matrices given by M, C, and K respectively, U the vector that contains the nodal displacements, and F the applied nodal forces. Its time integration can be performed by using any well experienced state of the art discretization technique, as [2] or [3].
In what follows, we will omit the damping term that results from the fact of assuming a proportional damping that expresses it as a combination of the mass and stiffness contributions.
To enhance the integration efficiency, mass lumping is usually considered, leading to a mass diagonal matrix. Model analysis looks also for enhancing the solution efficiency by decoupling the motion equation. For that purpose, the last extracts the basis {φ 1 , φ 2 , . . . , φ N } (N being the problem size, i.e., the number of degrees of freedom) by solving the eigenproblem −ω 2 M + K φ = 0, associated with the dynamical problem expressed in the Fourier space, where U and F refer to the Fourier transform of the nodal displacement U and forces F. The eigenmodes φ i , i = 1, · · · , N define an orthogonal basis, normalized with respect to the mass matrix, i.e., with δ the Kroenecker delta, and With P the matrix composed by the eigenmodes, i.e., P = (φ 1 , · · · , φ N ), the matrix form of the previous expressions reads P T MP = I and P T KP = K, with I the identity matrix and K the diagonal matrix with entries K ii = κ i .
In the modal basis U = Pϕ, the dynamical problem reads I d 2 ϕ(t) dt 2 + Kϕ(t) = P T F(t), (6) which constitutes a system of N uncoupled second order ordinary differential equations. The main limitation of modal analysis is the lack of validity of such basis in the case of parametric models. In the case of parametrized dynamical systems, with the model parameters grouped in vector µ, the model matrices will depend on those parameters, i.e., M(µ) and K(µ). The solution of parametric eigenproblems remains a tricky issue.
When one is not really interested in the transient regime, but much more in the forced regime, harmonic analysis represents a valuable route as considered in our former works where the so-called Proper Generalized Decomposition-PGD-enabled considering the frequency as a model extra-parameter as well as addressing general (non-proportional) damping and nonlinear dynamics, under the stringent real-time constraint, with even the inclusion of model parameters as extra-coordinates [4][5][6].
However, certain applications need accurate transient responses, and in that case the formulation and solution of the dynamical problem in the time domain is retained. Three routes are usually considered: 1.
The previously referred mass lumping that transforms the so-called consistent mass matrix into its diagonal counterpart, facilitating an explicit integration; 2.
In the context of model order reduction-MOR-, in [7], authors proposed a Proper Orthogonal Decomposition-POD-based reduced order modeling operating in the time domain. Ladeveze and coworkers proposed an extension of their radial approximation [8] for addressing mid-frequency dynamics, the so called TVCR (variational theory of complex rays) [9]. In our former works, we considered a PGD formulation for constructing a parametric transfer function [10] that allowed efficient solutions of transient dynamics operating in the time-domain. On the other hand, the separation of variables, at the heart of PGD [11], was extensively employed in the harmonic domain for solving multi-parametric dynamics, and was successfully extended to the nonlinear case, and then combined with modal analysis [4][5][6].

3.
Modal analysis is one the most widely used techniques for solving dynamical problems. Other than the benefits in the time integration, due to the dynamical system decoupling, the eigenmodes benefit from a physical interpretation, of great interest for the designer or structural analyst. However, when considering parametric models, as is always the case during the design stage, when the material and geometry are not totally defined, the dynamical modes depend on those parameters as previously discussed. Having a surrogate model expressing the parametric evolution of the eigenmodes is of great interest. Constructing those surrogate models is nowadays quite mature, by using usual and advanced nonlinear regressions [12], the last making use of sparsity and appropriate regularizations for operating in high-dimensional settings, while keeping as reduced as possible the number of data (eigenproblem solution), and leading to rich enough (nonlinear) regressions while avoiding overfitting.
Here, the trickiest issue is not the regression construction but the fact of ordering the different eigenmodes involved in the modal basis for each parameter choice, in order to create N clusters (or less in the reduced case), and putting in each one a mode of each modal basis, such that modes in each cluster remain close (in certain metrics). The main issue remains the metric to be used to successfully and efficiently accomplishing such clustering. In general, such clustering is performed by operating at the level of the eigenmodes, in the associated vector space, by using, for example, the Modal Assurance Criterion-MAC- [13] that proceed with comparing the modes resulting from each eigenproblem, by using the usual scalar product (modes similar to a given one should remain quite collinear).
The present paper focuses on the techniques based on modal analysis. As just described, usual techniques operate at the eignemodes level, defined in a vector space. When operating in high dimensional parametric spaces, sparsely sampled, the matrices involved in the resulting eigenproblem can vary a lot from one choice of the parameters to another, and consequently the scalar product criterion at the basis of the MAC could fail. On the other hand, the fact of proceeding in a vector space needs to carefully address the expression of the different modes by considering the same frame for all the analyzed mechanical systems.
For alleviating those difficulties, this paper proposes an alternative technique that instead of operating at the eigenmodes level, classifies the modes with respect to the deformed structures shape that the eigenmodes produce, by invoking the so-called Topological Data Analysis-TDA-that benefits from the invariance property of topology. The present paper does not aim at addressing the interpolation between the original and deformed structures that could be addressed by using optimal transport, or the parametrization of the manifold as in the parallel transport [14]. The present work aims only at addressing the deformed structures comparison using an appropriate metric, in particular the one based on the employment of TDA.
After the above short introduction, Section 2 addresses the problem and methodologies employed, Section 3 reports the obtained results, and Section 4 discusses the obtained results emphasizing the added value of the proposed approach.

Methods
As just indicated, the present paper aims at addressing the classification of a series of modal basis related to the eigenmodes of a thin structure equipped with a mesh consisting of shell elements, with displacement and rotation degrees of freedom at each node of the mesh.
The thickness of the structural part varies, with its consequent effect on the mass and stiffness matrices (damping is assumed proportional) and consequently on the eigenvalues and eigenvectors, the former defining the number of modes to be retained in the reduced basis. In the present study, the six rigid modes representing the whole structure translation (three modes) and rotation (three modes) will be discarded and among the remaining pairs eigenvalue-eigenvector, the most relevant six eigenvectors (corresponding to the six highest eigenvalues) retained in the reduced basis related to each choice of the model parameter (the thickness).
These six modes related to each structure (related to a thickness value) define a reduced basis that one would like render parametric. However, prior to constructing a regression able to define the reduced basis for each possible choice of the parameter (thickness), one should classify the six eigenmodes of each reduced basis associated with each structure into six clusters.
This task is compulsory to facilitate the interpolation in the parametric space and also to attach a physical sense to those modes. One could imagine that, for a given thickness, the most relevant deformation mode could be related to the extension, whereas, for another choice of the thickness, the most relevant deformation mode could be the bending. In such a case, one prefers to create a cluster grouping similar deformation modes, to evaluate how each of them depends on the parameter from one side, and on the other to facilitate the subsequent construction of the parametric modal reduced basis.
To preform such a clustering, we must employ an appropriate metric to compare those modes. In general, this comparison was traditionally performed by comparing the eigenmodes within the vector space to which they belong. In the present work, as announced previously, we prefer applying the deformation mode to the reference (undeformed) structure, that is, applying the eigenmode at the nodes location in the reference structure for obtaining the deformed structure related to each mode of each structure configuration (thickness) and then clusterizing the resulting deformed structures with respect to their shape.
Thus, we are employing a metric able to compare shapes, more than a metric for comparing the vectors (eigenmodes) that produced those shapes, the last being more intrinsic and inheriting invariance features. Moreover, in the present case study, eigenmodes are heterogeneous in the sense that they involve displacements and rotations, whereas the associated deformed surfaces are purely geometrical. For that purpose, the TDA will be applied.
In the remaining part of the present section, we will describe the available data and its organization, and then all the concepts enabling the use of a metric applying on the shapes, based on the employment of persistent homology, at the heart of the TDA, employed in our former works [15,16].

Data Description
As the different analyzed structures are shells of different thicknesses, from now on, we will describe these structures by their surfaces, each equipped of a node distribution and the associated mesh.
In this study, we consider a collection of M = 102 surfaces corresponding to the effect of a given deformation mode on the reference undeformed surface, as Figure 1 shows. Each surface M r , r = 1, . . . , M is equipped with a mesh associated with N = 3636 nodes, each node described by x n ≡ (x n , y n , z n ), n = 1, . . . , N and x n ∈ R 3 , all them making use of the same common coordinates frame.
The deformed structures consist of the nodes and elements resulting from the reference one by applying the associated deformation mode. There is neither nodes' redistribution nor refinement in the deformed surfaces. Figure 2 depicts the reference surface and the nodes distribution on it, from which all the deformed structures with their associated nodal distribution and deformed mesh will result. It is important to note that the undeformed and deformed meshes (elements connectivity) remain unchanged, facilitating the use of the proposed metrics discussed later. The M = 102 surfaces are associated with 17 different structures, each one of them having a different thickness, with the consequent effect on the mass and stiffness matrices, and therefore on the resulting eigenfrequencies and eigenmodes. For each of the 17 structures, the six most significant deformation modes (related to the six highest eigenvalues with the rigid modes excluded) are retained. As mentioned, by applying this 17 × 6 deformation modes to the original undeformed reference surface, the M = 102 deformed surfaces result.
The data are structured as a table, and, within each row, the six deformation modes related to a given structure (with its own thickness) are as follows: Our main aim in what follows is ordering the elements in the columns, such that each column will represent a similar deformation behavior.

On the Surface Topology
Consider a data-set M related to a given deformed surface defined from its N nodes, all of them in R 3 . We are interested in extracting the geometric features of M and how they are distributed across the different spatial scales.

Geometric Features
In order to describe the geometry of the data-set M, we first identify four types of geometrical features associated with M: where λ m , λ n and λ l are the barycentric coordinates, with λ m + λ n + λ l = 1; , and x m − x p are linearly independent, and then: where λ m , λ n , λ l and λ p are the barycentric coordinates, with λ m + λ n + λ l + λ p = 1. The vertices represent the dimension-0 features, segments the dimension-1 features, and triangles the dimension-2 features and tetrahedron dimension-3 features.

Features' Filtration
To describe the appearance and disappearance of the features of M across different scales, we consider the so-called Alpha Filtration. For that purpose, an interval [α min , α max ] is considered. It reflects the smallest features' scale, in our case α min = 0 (the vertex) and the largest one α max representing the largest distance between points of M.
The If two simplices in S d (M) have a common element σ, then there exists 0 ≤ l ≤ (d − 1) such that σ ∈ S l (M). Given the scale values (α j ) m j=0 , the Alpha Filtration is then a non-decreasing sequence describing the evolution of the features of the simplicial complex S(M) at each scale α j , and computed as follows: • First, the filtration value of each tetrahedron is computed as the circumradius of the tetrahedron if its circumsphere is empty, and as the minimum of the filtration values of the triangles that are within the circumsphere otherwise. • Similarly, the filtration value of each triangle is computed as the circumradius of the triangle if the circumcircle is empty, and as the minimum of the filtration values of the segments that are within the circumcircle otherwise. • Then, the filtration value of each segment is computed as its circumradius. • Finally, the filtration value of the vertices is set to 0.
The discrete values used for the radii are the α j , and all simplices that have a filtration value larger than α max are discarded.
The time complexity of the algorithm is O(n 2 ). The choice of the Alpha Filtration was motivated by its relatively much smaller size compared to other filtrations. A detailed definition and implementation are provided in [17,18].

Persistence Diagrams
In order to have a more exhaustive view on how the features are changing across different scales, the appearance and disappearance of each feature within the filtration are tracked and coded into the Homology Groups H k (M), where k = 0, 1, 2, 3 is the homology dimension.
The elements of a Homology Group H k (M) are classes of chains of simplices ("packets") σ ∈ S k (M), that is, simplices sharing faces, edges, or vertices. It can be seen as a connected component of the intersection of M with a k-dimensional linear subspace of R 3 . The use of Homology Groups allows for performing algebraic operations over their elements.
Given a Homology Group, we can now define how to track the appearance of the features across different scales, by defining the Homology Group at a scale α, H α k (M). It represents the classes of simplices as described previously, but taken from S α k (M). That is, the elements of S k (M) with a filtration value lesser that α.
This approach is known as the Persistent Homology. It allows for quantifying the appearance and disappearance of the features across the different scales and dimensions: • The birth scale b γ of the feature γ at homology k The birth scale represents the value at which the feature appeared in the filtration by combining lower dimensional simplices to form it. Conversely, the death scale represents the value at which the feature disappeared in the filtration by being combined into a higher dimensional feature. For example, if a vertex is part of a segment, then the death scale of the vertex is exactly the birth scale of the associated segment.
Note that, by definition, vertices always have a zero birth scale, while tetrahedrons always have an infinite death scale (in the numerical results, we removed the infinite values for computation purposes). Given that our data points M are embedded in R 3 , we will only track up to the dimension-2 features, thus the definition of S(M) with k = 0, 1, 2. More generally, if the data points are embedded in a d-dimensional manifold, the persistent homology can be computed up to dimension d − 1.
Finally, the persistence of the features throughout the scales can be represented by the so-called Persistence Diagram of M, defined at dimension-k from

Illustrating the Concepts on an Example
We illustrate the computational aspects of the Alpha Filtration on a simple example, consisting of six points in R 3 , as shown in The filtration values are computed and presented below in Table 1: We can then track the birth and death of the features and compute the persistence diagrams P D(M), as shown in Figure 4.

Matching Persistence Diagrams
Consider two data-sets M r and M s representing two deformed configurations of the same surface. A matching between two persistence diagrams with the same number of features, P D k (M r ) and P D k (M s ), for k = 0, 1, 2, is a bijective map ψ k that reads: The map ψ k associates each feature from P D k (M r ) to a feature from P D k (M s ). The Optimal Matching between P D k (M r ) and P D k (M s ) is a matchingψ k ψ k : P D k (M r ) −→ P D k (M s ), minimizing the transport cost C k to move the features from P D k (M r ) to P D k (M s ): When M r is the reference surface, and M s is any deformed surface resulting from M r , the optimal matchingψ k represents and quantifies the deformation from a topological viewpoint, at each dimension k = 0, 1, 2.
We note that, in our case considered here, the diagrams have all been reduced to their top 3000 persistence values, making the bijective matching possible. In the case of diagrams with a different number of points, a partial matching is rather considered.
The optimal matching is computed using a combinatorial optimization procedure, where the points in both diagrams are matched while minimizing the transport cost function defined above. A graphical representation of the matching is shown in Figure 5.

Multi-Scale Topological Measure of Surface Deformation
It is now possible to measure the degree of deformation from one data-set to another. For that purpose, consider two data-sets M r and M s representing two deformed states of the same surface, and a finite sequence of (α j ) m j=0 . Then, for k = 0, 1, 2, the k-distance between P D k (M r ) and P D k (M s ) reads whereψ k is the optimal matching between P D k (M r ) and P D k (M s ). An efficient computation of that distance W k , known as the Wasserstein Distance, is performed using the kernel linearization presented in Algorithm 2 of reference [19]. The Multi-Scale Topological Distance between M r and M s reads where ω k = W k P D k (M r ), P D k (M s ) , k = 0, 1, 2.

Comparing Tropological Descriptions of Deformed Surfaces
Consider now our collection {M 0 , M 1 , . . . , M M } of data-sets, consisting of the reference surface M 0 and M = 102 deformed surfaces M r , r = 1, . . . , M. By using the previously defined metric Ω, we can measure each surface deformation with respect to the reference one, such that ∀r ∈ {1, · · · , M} and we have where ∀r ∈ {1, . . . , M}, ∀k ∈ {0, 1, 2}, we denote The measure Ω r enables clustering the different surfaces (6 × 17, 17 being the number of structural configurations, each one related to a particular value of the thickness) into six clusters.

Modal Assurance Criterion
One of the most popular tools for the quantitative comparison of modal vectors is the Modal Assurance Criterion (MAC) [13]. The MAC criterion is a statistical indicator that is quite sensitive to large differences of the eigenmodes. In our case, each mode M r (r = 1, . . . , M) is decomposed in its linear and angular components (with respect to the three coordinate axes), resulting in the six vectors {Ψ r c } 1≤c≤6 . The MAC of two surfaces M r and M s is then computed according to The MAC takes values between 0 (representing no consistent correspondence) and 1 (representing a consistent correspondence). Values larger than 0.9 indicate consistent correspondence, whereas small values indicate poor resemblance of the two eigenmodes.
Thus, considering the six modes related to two different structures (with two different thicknesses), it is now possible to compute the so-called MAC matrix M to compare the modes and identify resemblances or discrepancies.
The MAC matrix M becomes diagonal dominant when modes are well ordered, whereas the loss of that diagonal dominancy informs on eventual permutations. In order to apply the MAC criterion in the case study addressed in the present paper, the first reduced modal basis consisting of the six modes related to the first thickness will be compared with the six modes of all the other configurations, the remaining 16 thickness choices.

Topological Modes Identification
By applying the methodology described in Section 2.2 to the surface data-sets, by first computing the Alpha Filtration, we obtain the persistence diagrams {P D k (M r )} M=102 r=0 , with the ones associated with the reference surface shown in Figure 6.
The multi-scale distance defined in Section 2.2.6 and here associated with each surface {Ω r } M=102 r=1 is then computed, and reported in Table 2. This multi-scale topological distance is then used to order the deformed surfaces to retain in each column of Table 3 those exhibiting a shape resemblance. From Tables 2 and 3,  the surfaces have been sorted using the values of Table 2, and their sorted order displayed in Table 3. The goal is to have the surfaces labeled from the least to the most deformed, according to our measure of deformation. In Table 3, numbers in red indicate surface (modes) permutations that have been made in order to classify all the shapes.    1  2  3  4  5  6  03  1  2  3  4  5  6  04  1  2  3  5  4  6  05  1  2  3  4  5  6  06  1  2  3  5  4  6  07  1  2  3  4  5  6  08  1  2  3  4  5  6  09  1  2  4  3  6  5  10  2  1  3  4  6  5  11 1 2 3 4 5 6

MAC Identification
Using the MAC criterion described in Section 2.3, we compute the MAC matrices comparing the model reduced basis (of the first thickness choice) with the remaining 16 reduced modal bases for the other thickness choices, and the results are reported in Figure 7.

Discussion and Conclusions
Labelling the surfaces as reported in Table 3 aims at classifying them according their shape induced by their deformation. The greater the value of the topological metric Ω, the more deformed the surfaces are. The surface discrepancies are quantified from the transport cost related to the matching of the topological features that the deformed meshes express, through the different scales and dimension.
The value of Ω can be interpreted as a level of topological deformation for a certain deformed mesh on the deformed surface compared with the undeformed mesh and surface. Thus, labels 1 to 6 in the case addressed here express the magnitude of the surface deformation. Figure 8 depicts the six ordered deformed surfaces for one particular case (structure with a given thickness).  Table 3, it can be noticed that the surface label usually matches the order of the eigenmodes provided by the eigenproblem solution. However, when modifying the structure thickness, some shapes that were important for a given thickness can be now more or less significant and the order of apparition in the eigenproblem differs. Thus, a permutation must be performed for ordering the modes with respect to their intrinsic shapes, here evaluated by using a metric based on topological concepts.The MAC matrices displayed in Figure 7 show similar tendencies, as the modes are globally consistent with their labels.

By inspecting
The presented topology-based framework for measuring surface deformations seems a very pertinent, powerful, and intrinsic way of quantifying, characterizing, and analyzing the deformation modes of structures. The strength of the framework relies on both the topology description of the surface at multiple scales, and the proposed measure based on the optimal matching of the features, to detect how each feature of the surface was deformed.