Integer matrix factorization for mesh defect detection

Francesca Rapetti<sup>a</sup>, François Dubois<sup>b</sup>, Alain Bossavit<sup>c</sup>

Received and accepted 4 February 2002

Note presented by Philippe G. Ciarlet.

**Abstract** The topological features of a given domain Ω in R³ are here analyzed by means of the homology groups of first and second order. Algebraic topology together with a particular QR type factorization in Z can be used to know whether Ω is connected and simply connected, as well as to check if a given discretization of Ω by means of simplices has been correctly realized. To cite this article: F. Rapetti et al., C. R. Acad. Sci. Paris, Ser. I 334 (2002) 717-720. © 2002 Académie des sciences/Éditions scientifiques et médicales Elsevier SAS

Factorisation matricielle à nombres entiers pour détecter des défauts
dans les maillages

Résumé
Les caractéristiques topologiques d'un domaine Ω de R³ sont analysées ici à l'aide
des groupes d'homologie du premier et second ordre. La topologie algébrique et une
factorization particulière de type QR dans Z peuvent être utilisées afin de savoir si Ω est
connexe et simplement connexe, de même que pour vérifier si une discrétisation de Ω par
éléments simpliciaux a été bien réalisée. Pour citer cet article: F. Rapetti et al., C. R. Acad.
Sci. Paris, Ser. I 334 (2002) 717-720. © 2002 Académie des sciences/Éditions scientifiques
et médicales Elsevier SAS

**1. Introduction**

A precise description of complicated industrial geometries relies on the use of computer assisted design tools. In this framework, the automatic generation of simplified meshes is often affected by topological mistakes (e.g., “holes”). We propose here an application of algebraic topology, graph theory and linear algebra to investigate the presence of mesh defects in two and three dimensions. Note that triangulating a domain reduces the calculation of homology groups to a finite procedure. The importance of simplified meshes relies on the fact that any triangulated domain is homeomorphic to one in which the elements have flat faces and straight edges. The remarkable thing is that these groups, in spite of being defined through a triangulation, do measure something intrinsic and geometric (they are topological invariants, i.e., they depend on the domain up to a homeomorphism) that does not depend on the mesh, if the latter is not too coarse. The homology groups of a domain have a direct link with the possibility of representing

<footnote><sup>a</sup> Laboratoire de mathématiques J.A. Dieudonné, CNRS & Université de Nice et Sophia-Antipolis, Parc Valrose, 06108 Nice cedex 02, France</footnote>

<footnote><sup>b</sup> Laboratoire des applications scientifiques du calcul intensif, CNRS, bat. 506, Université Paris Sud, 91403 Orsay, France</footnote>

<footnote><sup>c</sup> EdF, Division recherche et développement, 1, avenue du Général de Gaulle, 92141 Clamart cedex, France</footnote>

<footnote>*E-mail addresses*: frapetti@math.unice.fr (F. Rapetti); dubois@asci.fr (F. Dubois); Alain.Bossavit@der.edfgdf.fr
(A. Bossavit).</footnote>
curl-free and div-free vectors as gradients and curls, respectively. They contribute to build up suitable “belted” spanning trees of edges or faces [6]. This link is determinant in the construction of numerical algorithms for solving given problems in terms of scalar or vector potentials such as in electromagnetism or fluid dynamics. In a spirit similar to the GAP (Groups, Algorithms and Programming) free system for computational discrete mathematics, well known in the group theory research field, we propose an automatic procedure to compute a basis of generators for homology groups to investigate the presence of topological mesh defects. The very same procedure is used in [3] and for a recent application of similar tools to classify surfaces, see [5].

## 2. Cellular complexes and homology groups

Given a domain $\Omega \subset \mathbb{R}^3$ with boundary $\Gamma$, a simplicial mesh $m$ in $\Omega$ is a tessellation of $\bar{\Omega}$ by tetrahedra, under the condition that any two of them may intersect along a common face, edge or node, but in no other way. We denote by $N_m, E_m, F_m, T_m$ (nodes, edges, faces, and tetrahedra, respectively) the sets of simplices of dimension 0 to 3 thus obtained, and by $N_m, E_m, F_m, T_m$ their cardinalities. The notation $m^\Gamma$ stands for the trace of $m$ on $\Gamma$ and all quantities related to $m^\Gamma$ will carry the superscript $\Gamma$.

First we need to underline some combinatorial properties of the simplicial mesh. Let $M(r, s)$ denote the set of matrices $A$ whose elements are $A(i, j)$ with $1 \le i \le r$ and $1 \le j \le s$. Besides the list of nodes and their positions, the mesh data structure also contains incidence matrices, saying which node belongs to which edge, which edge bounds which face, etc. [1,4]. Moreover, there is a notion of orientation of the simplices that has to be taken into account. In short, an edge is not only a two-node subset of $N_m$, but an ordered set where the order implies an orientation. Let $e = \{\ell, n\}$ be an edge of the mesh oriented from the node $\ell$ to $n$. We can define the incidence numbers $G(e, n) = 1$, $G(e, \ell) = -1$ and $G(e, k) = 0$ for all nodes $k$ other than $\ell$ and $n$. These numbers form a rectangular matrix $G \in M(E_m, N_m)$, which describes how edges connect to nodes. A face $f = \{\ell, n, k\}$ has three vertices which are the nodes $\ell, n, k$. Note that $\{n, k, l\}$ and $\{k, l, n\}$ denote the same face $f$ whereas $\{n, l, k\}$ denotes an opposite oriented face, which is not supposed to belong to $F_m$ if $f$ does. An orientation of $f$ induces an orientation of its boundary. So, with respect to the orientation of the face $f$, the one of the edge $\{l, n\}$ is positive and negative the one of $\{k, n\}$. Then, we can define the incidence number $R(f, e) = 1$ or $-1$ if the orientation of $e$ matches or does not match, respectively the one on the boundary of $f$ and $R(f, e) = 0$ if $e$ is not an edge of $f$. These numbers form a matrix $R \in M(F_m, E_m)$. Finally, let us consider the tetrahedron $t = \{k, l, m, n\}$, positively oriented if the three vectors $\{k, l\}, \{k, m\}$ and $\{k, n\}$ define a positive frame ($t' = \{l, m, n, k\}$ has a negative orientation and does not belong to $T_m$ if $t$ does). A third matrix $D \in M(T_m, F_m)$ can be defined by setting $D(t, f) = \pm 1$ if face $f$ bounds the tetrahedron $t$, with the sign depending on whether the orientation of $f$ and of the boundary of $t$ match or not, and $D(t, f) = 0$ in case $f$ does not bound $t$. It can be easily proven that $RG = 0$ and $DR = 0$ [1].

Let $m$ be the simplified mesh on $\Omega \subset \mathbb{R}^3$. A $p$-chain $c$ is a function that associates a number $\alpha_i$ to each simplex of dimension $p$ in $m$ and we denote by $C_p(m)$ the set of all $p$-chains. The set $C_p(m)$ has a structure of Abelian group with respect to the addition of $p$-chains (two $p$-chains are added by adding the coefficients).

As an example, let us consider a path of edges of the mesh $m$ to go from a point $n_1$ to a point $n_2$: the path is an oriented line. Assigning an integer $\alpha_e$ equal to $+1$ or $-1$ when the edge $e$ belongs to the path and its orientation is in agreement or in disagreement with that of the path and 0 for all edges $e$ that do not belong to the path, we define a 1-chain. A circuit is a line plus a way to run along it; so, when the line is made of oriented edges, we need to tell the positive direction along each edge, which is precisely what the chain coefficient $\alpha_e$ does. In our case, we assume that all coefficients $\alpha_i$ are relative integers, i.e., $\alpha_i \in \mathbb{Z}$.

Next concept is the boundary operator $\partial_p : C_p(m) \to C_{p-1}(m), p > 0$. By definition, we have:

$$ \partial_1(e) = \sum_{n \in N_m} G(e, n)n, \quad \partial_2(f) = \sum_{e \in E_m} R(f, e)e, \quad \partial_3(t) = \sum_{f \in F_m} D(t, f)f. $$
Note that $\partial_p$ is represented by a matrix that is $G^t$, $R^t$ and $D^t$ depending on the dimension $p > 0$. We remark in particular that $\partial_{p+1} \circ \partial_p = 0$, i.e., the boundary of a boundary is zero.

We will say that a *p*-chain *c* is closed if $\partial_p(c) = 0$. Closed *p*-chains are called *p*-cycles and constitute the subset $Z_p(m) = \ker(\partial_p; C_p(m))$. A *p*-chain *c* is a boundary if there is a (*p* + 1)-chain $\gamma$ such that $c = \partial_{p+1}(\gamma)$. The *p*-boundaries constitute the subset $B_p(m) = \partial_{p+1}(C_{p+1}(m))$. Both $Z_p(m)$ and $B_p(m)$ are Abelian groups with respect to the addition of *p*-chains. Boundaries are cycles but not all cycles are boundaries: we have in fact that $B_p(m) \subset Z_p(m)$.

Our concern is to determine the 1-cycles and 2-cycles that are not boundaries, i.e., to define the generators of the homology group of order $p$ of the mesh, $H_p(m) = Z_p(m)/B_p(m)$, for $p=1$ and $2$. To understand the link between these groups and the detection of topological mesh defects, let us consider the situation where two separately meshed parts of a domain $\Omega$ are glued together along an interface $\mathcal{I}$ to give the global discretization of $\Omega$ (this occurs in general when $\Omega$ has complicated geometric features). We suppose that, after the gluing process, a couple of nodes, both on $\mathcal{I}$ and each of them on a different part of $\Omega$, come very close to each other but do not coincide, by mistake. In two dimensions, the 1-chain associated to the path of edges having one of these two nodes as extremity does not bound a cluster of triangles: it is an element of $H_1(m)$. In three dimensions, the 2-chain associated to the path of faces having one of these two nodes as extremity does not bound a cluster of tetrahedra: it is an element of $H_2(m)$.

A very first information on the connectivity genus of $\Omega$ and of its surface $\Gamma$ comes from the Euler-Poincaré characteristic [4]. For a given connected domain $\Omega$, the Euler-Poincaré characteristic of $\Omega$ is the integer $\chi(\Omega) = N_m - E_m + F_m - T_m$ and for a connected surface $\Gamma$, it is $\chi(\Gamma) = N_\Gamma - E_\Gamma + F_\Gamma$. The Euler-Poincaré characteristic is linked to the homology groups as follows: setting $b_p = \dim H_p(m)$ (resp. $b_p^\Gamma = \dim H_p(m^\Gamma)$), $p = 0, 1, 2, 3$ (resp. $p = 0, 1, 2$) for the Betti numbers of $\Omega$ (resp. of $\Gamma$), we have $\chi(\Gamma) = b_0^\Gamma - b_1^\Gamma + b_2^\Gamma$ and $\chi(\Omega) = b_0 - b_1 + b_2 - b_3$. The major point is that Betti numbers, and hence the Euler-Poincaré characteristic, are topological invariants: they do not depend on the adopted discretization. The constant $\chi(\Omega)$ is typically 0, 1 or 2. For the surface $\Gamma$ we have that $\chi(\Gamma) = 2(1 - r(\Gamma))$, where $r(\Gamma)$ is the genus of $\Gamma$ (i.e., the “number of loops” of $\Gamma$ or, more precisely, the maximum number of disjoint closed curves which can be drawn on $\Gamma$ without separating it). For more details on the subject, see [7].

### 3. An integer matrix factorization

In this section we present a matrix decomposition to compute the homology groups of order $p=1$ and $2$ of $\Omega \subset \mathbb{R}^3$. The basic idea is to make an integer $QR$ factorization of the matrices $G^t$, $R^t$ and $D^t$. Given a matrix $A \in \mathcal{M}(r, s)$, we compute a non-singular unimodular matrix $Q$ and a permutation matrix $P$ such that $R = QA P$ is upper triangular. As shown later, the two matrices $Q$ and $P$ are obtained as products of a certain number of local matrices $Q_{i,j}$ and $P_{i,j}$ and exhibit the row and column rank deficiency of $A$ [2].

To define $Q$ and $P$, we need two elementary operations. Firstly, the transformation $\Pi_1$ of a vector $v = (\varepsilon_i, \varepsilon_j)^t$ into the vector $\tilde{v} = (1, 0)^t$. To this purpose, let us introduce the matrices

$$Q_{i,j}^{\text{el}} = \begin{pmatrix} \varepsilon_i & 0 \\ -\varepsilon_i & \varepsilon_j \end{pmatrix}, \quad Q_{i,j}(\ell, q) = \begin{cases} \delta_{\ell,q}, & \ell \neq i, j, \quad q \neq i, j, \\ Q_{i,j}^{\text{el}}(1,1), & \ell = i, \quad q = i, \\ Q_{i,j}^{\text{el}}(1,2), & \ell = i, \quad q = j, \\ Q_{i,j}^{\text{el}}(2,1), & \ell = j, \quad q = i, \\ Q_{i,j}^{\text{el}}(2,2), & \ell = j, \quad q = j. \end{cases} \quad (Q_{i,j}^{\text{el}})^{-1} = \begin{pmatrix} \varepsilon_i & 0 \\ \varepsilon_j & \varepsilon_j \end{pmatrix}$$

In our case, $\varepsilon_i^2 = 1$, and the vector $\tilde{v} = \Pi_1(v) = Q_{i,j}^{\text{el}} v$. Secondly, we need the permutation $\Pi_2$ of a vector's components, i.e., the transformation of a vector $v = (\varepsilon_i, \varepsilon_j)^t$ into the vector $\tilde{v} = (\varepsilon_j, \varepsilon_i)^t$. To this purpose, we have $\tilde{v} = \Pi_2(v) = P_{i,j}^{\text{el}} v$, where $P_{i,j}^{\text{el}}$ is a permutation matrix; moreover, we introduce a matrix $P_{i,j}$ defined similarly to $Q_{i,j}$ (using $P_{i,j}^{\text{el}}$ instead of $Q_{i,j}^{\text{el}}$). We remark that $(P_{i,j}^{\text{el}})^{-1} = P_{i,j}^{\text{el}}$, owing to the fact
that $\mathcal{P}_{i,j}^{\text{el}}$ is a permutation matrix and that $(Q_{i,j})^{-1}$ is defined as $Q_{i,j}$ (using $(Q_{i,j}^{\text{el}})^{-1}$ instead of $Q_{i,j}^{\text{el}}$).

Now we describe the adopted procedure to build up $\mathcal{Q}$ and $\mathcal{P}$ for a given matrix $A \in \mathcal{M}(r, s)$.

*Procedure:* We set $\mathcal{Q} = \mathcal{Q}^0 \in \mathcal{M}(r, r)$, $\mathcal{P} = \mathcal{P}^0 \in \mathcal{M}(s, s)$. We loop on the column index $j$, $1 \le j \le s$:

1. We define $\mathcal{V}_j = \{i \mid \min\{j, r\} \le i \le \min\{s, r\}, A(i, j) \ne 0\}$ and we put $k$ equal to the cardinality of $\mathcal{V}_j$, $i_1$ equal to the smallest integer in $\mathcal{V}_j$ and $i_2$ equal to the smallest integer in $\mathcal{V}_j \setminus \{i_1\}$.

2. In case $k=0$: let $\mathcal{P}_{j,z}$ be the matrix of the transformation $\Pi_2$ that permutes the $j$-th column of $A$ with the $z$-th one. The $z$-th column is chosen to be the first column, starting from the last one in $A$, for which it exists a row index $s$ such that $A(s, z) \ne 0$. If the index $z$ exists, $\mathcal{P} \leftarrow \mathcal{P} \mathcal{P}_{j,z}$, $A \leftarrow A \mathcal{P}_{j,z}$ and we go back to step 1, otherwise we stop the *Procedure*.

3. In case $k \ne 0$ but $A(j, j) = 0$, we apply a partial pivot strategy. Let $Q_{j,i_1}$ be the matrix of the transformation $\Pi_2$ that permutes the $j$-th row with the $i_1$-th one; then, $Q \leftarrow Q_{j,i_1} Q$, $A \leftarrow Q_{j,i_1} A$, $i_1 \leftarrow j$ and we go to step 4.

4. In case $k \ge 2$ and $A(j, j) \ne 0$, let $Q_{i_1,i_2}^{\text{el}}$ be the matrix of the transformation $\Pi_1$ applied to the vector $(A(i_1, j), A(i_2, j))^t$ and $Q_{i_1,i_2}$ the associated matrix, then $Q \leftarrow Q_{i_1,i_2} Q$, $A \leftarrow Q_{i_1,i_2} A$ and we go back to step 1.

5. In case $k=1$ and $A(j, j) \ne 0$, then $j \leftarrow j+1$ and we go back to step 1.

Starting with $\mathcal{Q}^0 = I_r$ and $\mathcal{P}^0 = I_s$, at the end of the *Procedure*, the matrix $A$ has been replaced by $\mathcal{R}$, an upper triangular one. If this new matrix $\mathcal{R}$ does not contain zero rows, then $\text{dim}(\text{range}(\mathcal{R})) = r$. Otherwise $\text{dim}(\text{ker}(\mathcal{R}^t)) = r - \text{dim}(\text{range}(\mathcal{R}))$. We remark that the procedure converges and its computational cost is similar to the one of a $\mathcal{Q}\mathcal{R}$ decomposition by using Givens transformations. Now, the question is to apply the *Procedure* to compute the generators of the homology group of order $p=1$ and $2$.

(i) We apply the *Procedure* with $A = R^t$, $\mathcal{Q}^0 = I_{E_m}$, $\mathcal{P}^0 = I_{F_m}$ and we get two invertible matrices $\mathcal{Q}_R$ and $\mathcal{P}_R$ such that $\mathcal{R}_R = \mathcal{Q}_R R^t \mathcal{P}_R$ is upper triangular. The 1-cycles which bound a surface belong to the image of the matrix $R^t$ that is also the image of $\mathcal{R}_R$.

(ii) We define $\tilde{G}^t = G^t \mathcal{Q}_R^{-1}$: in this way we make a change of basis for the 1-chains. Looking at $\tilde{G}^t$ we see immediately from the presence of $n_c$ zero columns that the corresponding columns of $\mathcal{Q}_R^{-1}$ represent vectors that belong to the kernel of $G^t$. If $\text{dim}(\text{range}(R^t)) = \text{dim}(\text{ker}(G^t))$ then any 1-cycle bounds. In the other case, we apply the *Procedure* with $A = \tilde{G}^t$, $\mathcal{Q}^0 = I_{N_m}$, $\mathcal{P}^0 = I_{E_m - n_c}$. We then obtain two invertible matrices $\mathcal{Q}_G$ and $\tilde{\mathcal{P}}_G$ such that $\tilde{\mathcal{R}}_G = \mathcal{Q}_G \tilde{G}^t \tilde{\mathcal{P}}_G$ is upper triangular.

(iii) The rows in $\tilde{\mathcal{P}}_G$, corresponding to zero rows in $\tilde{\mathcal{R}}_G$, represent the vectors that complete the kernel of $G^t$. In fact, we are looking for $c$ such that $G^t c = 0$. This is equivalent to $\tilde{G}^t v = 0$ where $v$ has zero in the first $n_c$ components and, in the last $E_m - n_c$, any row in $\tilde{\mathcal{P}}_G$ corresponding to a new zero row in $\tilde{\mathcal{R}}_G$. Then, the components of $c = \mathcal{Q}_R^{-1} v$ are the coefficients of a 1-chain, generator of $H_1(m)$.
To determine the generators of $H_2(m)$, it is sufficient to perform parts (i), (ii), (iii) with $D^t$ at the place of $R^t$ and $R^t$ at the place of $G^t$.

## References

[1] A. Bossavit, Computational Electromagnetism: Variational Formulations, Complementarity, Edge Elements, Academic Press, New York, 1998.

[2] T.F. Chan, Rank revealing QR factorizations, Linear Algebra Appl. 88/89 (1987) 67-82.

[3] F. Dubois, F. Rapetti, Du tourbillon au champ de vitesse, in: Workshop at the Conservatoire des Arts et Métiers de Paris, on the subject "Modèles fluides et représentation en toubillons", Vol. 1, 2000, pp. 127-153, see also in MATAPLI 65 (2001) 87-88.

[4] P.J. Giblin, Graphs, Surfaces and Homology. An Introduction to Algebraic Topology, Chapman and Hall Math. Ser., 1977.

[5] R. Hiptmair, J. Ostrowski, Generators of $H_1(\Gamma_h, Z)$ for triangulated surfaces: Construction and classification, Sonderforschungsbereich 382, Universität Tübingen, Report 160, 2001.

[6] L. Kettunen, K. Forsman, A. Bossavit, Discrete spaces for divergence- and curl-free fields, IEEE Trans. Magn. 34 (5) (1998) 2851-2854.

[7] J. Stillwell, Classical Topology and Combinatorial Group Theory, Graduate Texts in Math., Vol. 72, Springer-Verlag, 1993.