Triangulation de Delaunay

3.a. Insertion d'un nœud interne

L'algorithme de Bowyer-Watson permet d'ajouter à une triangulation un nouveau nœud et de mettre à jour la triangulation de Delaunay de manière efficace. Dans un premier temps, on traite le cas où le nœud ajouté se trouve dans la triangulation, c'est-à-dire dans l'enveloppe convexe de l'ensemble des nœuds.

Soit N le nouveau nœud. La première étape consiste à repérer un triangle dont le disque ouvert circonscrit contient N. Un tel triangle existe évidemment puisque N se trouve dans l'enveloppe convexe des nœuds.

La deuxième étape consiste à trouver tous les triangles dont le disque ouvert circonscrit contient N. Pour obtenir ces triangles on les recherche tout d'abord parmi les voisins du premier triangle, puis parmi les voisins des voisins et ainsi de suite. La recherche se fait donc de manière récursive.

Les triangles dont le disque ouvert circonscrit contient N forment une cavité polyhédrique, comme illusté sur la figure suivante. Plus précisément, enlever ces triangles donne une cavité polyhédrique.

La troisième étape consiste à générer des nouveaux triangles. Il s'agit des triangles formés du nœud N et des nœuds de la cavité polyhédrique. Ces triangles sont Delaunay. La nouvelle triangulation est donc Delaunay. Pour générer ces triangles, il suffit de repérer le côté des triangles enlevés (les triangles de la cavité) qui ne sont partagés par aucun autre côté des autres triangles enlevés. Ces côtés sont en effet ceux du polyhèdre.

La figure suivante montre une triangulation de Delaunay avec un nœud ajouté N. Les triangles dont le disque ouvert circonscrit contient N sont grisés et leur cercle tracé. Les bords de la partie grisée forment la cavité.

Figure pleine page

La figure suivante montre les nouveaux triangles (grisés) et leur cercle circonscrit, qui effectivement ne contiennent aucuns nœuds.

Figure pleine page

Pour résumer, voici les quatre parties de l'algorithme :

• (1) Recherche d'un triangle dont le disque ouvert circonscrit contient N.
• (2) Recherche récursive de tous les triangles dont le disque ouvert circonscrit contient N. Ces triangles forment une cavité polyhédrique.
• (3) Retrait des triangles de la cavité.
• (4) Génération des triangles constitués de N et des côtés non partagés des triangles de la cavité.

La partie la moins efficace de l'algorithme est (1) car elle se fait a priori en itérant sur tous les nœuds de la triangulation. Il sera donc souhaitable d'accélérer cette étape en recherchant les triangles dans une zone particulière (mais cette optimisation dépendera de la méthode de maillage). Si le nouveau nœud est inséré dans un triangle particulier connu, cette recherche n'est évidemment plus nécessaire. On verra plus loin une méthode d'accélération de cette recherche dans le cas où les nœuds ajoutés sont connus à l'avance.

Nous ajoutons à la classe Maillage la fonction enleverTriangle, qui permet d'enlever un triangle de la liste des triangles, de la listes de voisins des triangles voisins et de la liste des triangles des nœuds sommets de ce triangle. La fonction ajouterNoeud permet d'ajouter un nœud et de mettre à jour la triangulation au moyen de l'algorithme de Bowyer-Watson.

class Maillage2(Maillage):
    def __init__(self):
        Maillage.__init__(self)
    def enleverTriangle(self,tr):
        self.retire.append(tr)
        self.triangles.remove(tr)
        for trv in tr.voisins:
            trv.enleverVoisin(tr)
        for i in range(3):
            if tr.noeuds[i]!=None:
                tr.noeuds[i].triangles.remove(tr)
    def ajouterNoeud(self,n):
        def rechercheCavite(n,tr,trCavite):
            for t in tr.voisins:
                if not(t in trCavite) and t.dansCercle(n):
                    trCavite.append(t)
                    rechercheCavite(n,t,trCavite)
                
        trCavite = []
        for tr in self.triangles: # recherche d'un triangle dont le cercle circonscrit contient le nouveau noeud
            if tr.dansCercle(n):
                trCavite.append(tr)            
                rechercheCavite(n,tr,trCavite) # recherche récursive des autres triangles
                break
        nouveaux = []
        
        for i in range(len(trCavite)):
            tr = trCavite[i]
            self.enleverTriangle(tr)
            
            for c in [[0,1],[0,2],[1,2]]: # recherche du côté non partagé
                na = tr.noeuds[c[0]]
                nb = tr.noeuds[c[1]]
                ajout = True
                for j in range(len(trCavite)):
                    if i!=j:
                        if trCavite[j].contientCote(na,nb):
                            ajout = False
                            break
                if ajout:
                    nouvTr = Triangle(n,na,nb)
                    nouveaux.append(nouvTr)
                    self.triangles.append(nouvTr)
                    
                    # recherche des voisins du nouveau triangle parmi les voisins des triangles enlevés
                    voisins = []
                    for t in trCavite:
                        for trv in t.voisins:
                            if not(trv in trCavite) and not(trv in voisins) and trv.triangleVoisin(nouvTr):
                                nouvTr.ajouterVoisin(trv)
                                trv.ajouterVoisin(nouvTr)
                                voisins.append(trv)
        # les nouveaux triangles sont voisins
        for i in range(len(nouveaux)):
            for j in range(len(nouveaux)):
                if i!=j:
                    nouveaux[i].ajouterVoisin(nouveaux[j])

    
                   

La recherche des voisins des nouveaux triangles se fait en deux étapes : on recherche tout d'abord les voisins parmi les voisins des triangles enlevés puis parmi les nouveaux triangles. Il est important qu'aucun voisin n'apparaisse plusieurs fois dans une liste de voisins.

Voici un exemple : on génére tout d'abord une triangulation de Delaunay aléatoire avec une enveloppe convexe carrée puis on ajute quelques nœuds de positions aléatoires :

  
points = []
noeuds = []
N = 10
x = uniform(-10,10,N)
y = uniform(-10,10,N)
for i in range(N):
    points.append([x[i],y[i]])
    noeuds.append(Noeud(x[i],y[i]))
for x,y in [(-10,-10),(-10,10),(10,10),(10,-10)]:
    points.append([x,y])
    noeuds.append(Noeud(x,y))
    
points = np.array(points)
tri = Delaunay(points)

figure(figsize=(8,8))
maillage = Maillage2()
maillage.fromList(noeuds,tri.simplices)
maillage.draw("k-")

            

fig3.pdf

 
N = 5
x = uniform(-10,10,N)
y = uniform(-10,10,N)
noeuds = []
for i in range(N):
    n = Noeud(x[i],y[i])
    maillage.ajouterNoeud(n)
    noeuds.append(n)
figure(figsize=(8,8))
maillage.draw("k-")
for n in noeuds:
    n.draw("ro")
invalid = maillage.verifierDelaunay()
            

fig4.pdf

print(invalid)
--> 0

3.b. Insertion d'un nœud externe

Un nœud externe n'est pas dans l'enveloppe convexe de l'ensemble initial. Dans ce cas, il s'agit d'ajouter des triangles qui comportent le nœud N ajouté et un segment de la frontière de la triangulation initiale (si le disque ouvert circonscrit du triangle associé à ce segment ne contient pas N). De plus, il faut supprimer les triangles dont le disque ouvert circonscrit contient N.

La figure suivante représente la situation dans le cas où aucun disque ouvert circonscrit ne contient N :

Figure pleine page

Les triangles à ajouter sont tracés en trait pointillé. Pour savoir avec quels segments de la frontière les nouveaux triangles doivent être créés, il suffit de déterminer de quel côté de la droite formé par le segment le nœud N se situe. Par exemple pour le segment N₁N₂, cette droite est tracée en trait pointillé rouge sur la figure. La droite en trait pointillé bleu est celle d'un segment avec lequel le nœud N ne formera pas un nouveau triangle.

Un nouveau triangle doit être créé si N se trouve à droite de la droite orientée (N₁N₂). La condition s'écrit :

$$\left(\overrightarrow{N_{1}N}\wedge \overrightarrow{N_1{N}_2}\right)\cdot \overrightarrow{u_z}>0\text{(3)}$$

Si N se trouve sur le segment N₁N₂ alors il se trouve aussi dans disque ouvert circonscrit du triangle auquel le segment appartient et ce triangle doit être enlevé (c'est ausi le cas si N est très proche du segment).

Le cas du nœud ajouté externe peut être intégré à l'algoritme de Bowyer-Watson en introduisant des triangles virtuels. À chaque segment de la frontière on associe un triangle virtuel dont le troisième sommet est un nœud virtuel. Ce nœud n'a pas de position particulière, son rôle est de permettre d'intégrer le segment dans un triangle qui sera traité comme tous les triangles par l'algorithme. La condition d'appartenance au disque ouvert circonscrit, qui permet de savoir si un triangle doit être enlevé, est remplacée par la condition $\mathrm{(3)}$. Si cette condition est vérifiée, le triangle virtuel est retiré. Les deux côtés virtuels de ce triangles deviennent des segments de bord s'il ne sont communs avec aucun autre côté de triangles enlevés. Si le nœud N se trouve sur le segment N₁N₂, il faut enlever le triangle réel auquel ce segment appartient mais aussi le triangle virtuel (car il aura deux segments à la place d'un seul). La condition pour que N appartienne au segment est :

$$\begin{array}{cc}\hfill & \left(\overrightarrow{N_{1}N}\wedge \overrightarrow{N_1{N}_2}\right)\cdot \overrightarrow{u_z}=0\hfill \\ \hfill & \overrightarrow{N_1{N}_2}\cdot \overrightarrow{N_{1}N}>0\hfill \\ \hfill & \overrightarrow{N_2{N}_1}\cdot \overrightarrow{N_{2}N}>0\hfill \end{array}$$

C'est cependant une condition théorique qui a peu de change de se produire à cause des erreurs d'arrondis. Il est bien plus probable que N se trouve à proximité du segment. La figure suivante montre cette situation après retrait des triangles de la cavité :

Figure pleine page

Dans ce cas, deux triangles sont enlevés, un triangle virtuel et un triangle réel, et quatre sont créés : deux triangles virtuels (associés au deux nouveaux segments de la frontière) et deux triangles réels (qui remplissent l'espace laissé par le triangle réel enlevé).

Nous devons implémenter une classe pour les triangles virtuels, qui comporte les mêmes fonctions que celle du triangle réel, en particulier la fonction dansCercle mais celle-ci fera le test décrit ci-dessus. On définit donc une nouvelle classe avec la même interface (en C++, il faudrait faire une classe virtuelle de base et deux classes dérivées) :

 
class TriangleVirtuel(Triangle):
    def __init__(self,n1,n2):
        self.noeuds = [n1,n2,None] # noeud virtuel = None
        self.voisins = []
        self.virtuel = True
        n1.triangles.append(self)
        n2.triangles.append(self)
    def triangleVoisin(self,tr):
        # test si un triangle est voisin
        return (self.noeuds[0] in tr.noeuds) or (self.noeuds[1] in tr.noeuds) 
    def ajouterVoisin(self,triangle):
        self.voisins.append(triangle)
    def enleverVoisin(self,triangle):
        self.voisins.remove(triangle)
    def dansCercle(self,n):
        x1 = self.noeuds[0].x
        y1 = self.noeuds[0].y
        x2 = self.noeuds[1].x 
        y2 = self.noeuds[1].y 
        test = (n.x-x1)*(y2-y1)-(n.y-y1)*(x2-x1)
        m = max(x1,x2,y1,y2)
        eps = 1e-11*m*m
        if test > eps :
            return True
        surSeg = -eps < test < eps and (x2-x1)*(n.x-x1)+(y2-y1)*(n.y-y1) > 0 and (x1-x2)*(n.x-x2)+(y1-y2)*(n.y-y2) > 0
        #print(surSeg,x1,y1,x2,y2,n.x,n.y)
        return surSeg
            
    def cercle(self):
        return None
    def contientCote(self,na,nb):
        return (na in self.noeuds) and (nb in self.noeuds)
    def coteCommun(self,tr):
        # test si un autre triangle comporte un côté commun avec le triangle
        return self.contientCote(tr.noeuds[0],tr.noeuds[1]) or self.contientCote(tr.noeuds[0],tr.noeuds[2]) or self.contientCote(tr.noeuds[1],tr.noeuds[2])
    def draw(self,style=None,linewidth=1):
        if style:
            plot([self.noeuds[0].x,self.noeuds[1].x],[self.noeuds[0].y,self.noeuds[1].y],style,linewidth=linewidth)
        else:
            plot([self.noeuds[0].x,self.noeuds[1].x],[self.noeuds[0].y,self.noeuds[1].y],linewidth=linewidth)
    def drawCercle(self,style):
        pass
    def print(self):
        print("Triangle : ",self)
        for n in self.noeuds[0:2]:
            print(n.x,n.y) 
                     

Pour mettre au point la nouvelle version de la classe Maillage, nous lui ajoutons une fonction triangleInitial, qui permet de créer un maillage comportant seulement un triangle réel et trois triangles virtuels, un pour chaque côté. Bien évidemment, on pourra plus tard définir une forme initiale quelconque.

Bien que l'algorithme soit le même, la fonction ajouterNoeud doit être légèrement modifiée car certains nouveaux triangles sont virtuels. L'algorithme permet de générer automatiquement les triangles virtuels associés aux nouveaux segments de la frontière. La figure suivante représente la situation où 3 segments de la frontière, et seulement 3, ne sont plus sur la frontière après la transformation. Le triangle virtuel associé à un segment de la frontière est représenté sous la forme d'une flèche, le sens de celle-ci indiquant l'ordre des nœuds dans le triangle virtuel (le nœud virtuel est toujours le troisième).

Figure pleine page

Par exemple, le triangle virtuel A est constitué des nœuds N₁,N₂,V où V désigne le nœud virtuel (représenté par l'objet None). Les triangles virtuels A,B et C sont dans la cavité et sont donc enlevés, indépendamment du fait que le segment correspondant reste ou pas dans la triangulation (il est retiré si le triangle réel auquel il appartient est dans la cavité). Un nouveau triangle est obtenu à partir d'un segment non partagé par d'autres triangles de la cavité auquel on ajoute le nœud N. Dans le cas du triangle virtuel A, les côtés N₁V et N₁N₂ ne sont pas partagés. Il y a donc un nouveau triangle réel N₁,N,N₂ et un nouveau triangle virtuel N₁,N,V (triangle D). Dans le cas du triangle B, les deux côtés N₂V et N₃V sont partagés respectivement avec les triangles A et C; seul son côté N₂N₂ n'est pas partagé donc il y a un nouveau triangle réel N₂,N,N₃ mais aucun triangle virtuel. Dans le cas du triangle C, les côtés non partagés sont N₃N₄ et N₃V donc il y a un nouveau triangle réel N,N₄,N₃ et un nouveau triangle virtuel N,N₄,V. Nous voyons donc que l'algorithme génère bien un triangle virtuel pour chaque nouveau segment de la frontière. Il faut remarquer que le nœud virtuel doit être considéré comm un seul et même nœud pour tous les triangles virtuels, ce qui est facile à obtenir si le nœud virtuel est toujours représenté par l'objet None. Lors de la création d'un triangle virtuel, il faut faire en sorte que ses deux premiers nœuds définissent un segment de la frontière parcouru dans le sens direct. Pour cela, il faut que le nœud du triangle enlevé (autre que V) qui appartient au nouveau triangle garde sa position dans le nouveau triangle, la troisième position étant toujours occupée par V et la position restante par N. Par exemple pour le nouveau triangle D, le premier nœud est N₁ puisque ce nœud avait la première position dans le triangle A. Pour le nouveau triangle E, le nœud N₄ est en deuxième position puisque c'était sa position dans le triangle C.

 
class Maillage3(Maillage):
    def __init__(self):
        Maillage.__init__(self)
    def enleverTriangle(self,tr):
        self.retire.append(tr)
        self.triangles.remove(tr)
        for trv in tr.voisins:
            trv.enleverVoisin(tr)
        for i in range(3):
            if tr.noeuds[i]!=None:
                tr.noeuds[i].triangles.remove(tr)
    def triangleInitial(self,n1,n2,n3):
        tr = Triangle(n1,n2,n3)
        self.triangles.append(tr)
        bord1 = TriangleVirtuel(tr.noeuds[0],tr.noeuds[1])
        bord2 = TriangleVirtuel(tr.noeuds[1],tr.noeuds[2])
        bord3 = TriangleVirtuel(tr.noeuds[2],tr.noeuds[0])
        bords = [bord1,bord2,bord3]
        for b in bords:
            self.triangles.append(b)
            b.ajouterVoisin(tr)
            tr.ajouterVoisin(b)
        bord1.ajouterVoisin(bord2)
        bord2.ajouterVoisin(bord1)
        bord1.ajouterVoisin(bord3)
        bord3.ajouterVoisin(bord1)
        bord2.ajouterVoisin(bord3)
        bord3.ajouterVoisin(bord2)
    def ajouterNoeud(self,n):
        def rechercheCavite(n,tr,trCavite):
            for t in tr.voisins:
                if not(t in trCavite) and t.dansCercle(n):
                    trCavite.append(t)
                    rechercheCavite(n,t,trCavite) 
                
        trCavite = []
        
        for tr in self.triangles: # recherche d'un triangle dont le cercle circonscrit contient le nouveau noeud
            if tr.dansCercle(n):
                trCavite.append(tr)            
                rechercheCavite(n,tr,trCavite) # recherche récursive des autres triangles
                break
        nouveaux = []
        for i in range(len(trCavite)):
            tr = trCavite[i]
            
            self.enleverTriangle(tr)
            
            for c in [[0,1],[0,2],[1,2]]: # recherche du côté non partagé
                na = tr.noeuds[c[0]]
                nb = tr.noeuds[c[1]]
                ajout = True
                for j in range(len(trCavite)):
                    if i!=j:
                        if trCavite[j].contientCote(na,nb):
                            ajout = False
                            break
                if ajout:
                    if nb!=None:
                        nouvTr = Triangle(n,na,nb) 
                    else :
                        if tr.noeuds[0]==na:
                            nouvTr = TriangleVirtuel(na,n)
                        else:
                            nouvTr = TriangleVirtuel(n,na)
                    nouveaux.append(nouvTr)
                    self.triangles.append(nouvTr)
                    
                    # recherche des voisins du nouveau triangle parmi les voisins des triangles enlevés
                    voisins = []
                    for t in trCavite:
                        for trv in t.voisins:
                            if not(trv in trCavite) and not(trv in voisins) and trv.triangleVoisin(nouvTr):
                                nouvTr.ajouterVoisin(trv)
                                trv.ajouterVoisin(nouvTr)
                                voisins.append(trv)
        # les nouveaux triangles sont voisins
        for i in range(len(nouveaux)):
            for j in range(len(nouveaux)):
                if i!=j:
                    nouveaux[i].ajouterVoisin(nouveaux[j])
        return len(nouveaux)

                      

 
maillage = Maillage3()
n1 = Noeud(-10,-10)
n2 = Noeud(10,-10)
n3 = Noeud(0,10)
maillage.triangleInitial(n1,n2,n3)

maillage.draw("k-")
                      

fig5.pdf

On commence par ajouter des nœuds à l'intérieur :

 
figure(figsize=(8,8))
maillage.ajouterNoeud(Noeud(0,0))
maillage.ajouterNoeud(Noeud(1,-5))
maillage.draw("k-")
                      

fig6.pdf

puis un nœud à l'extérieur :

figure(figsize=(8,8))
maillage.ajouterNoeud(Noeud(5,5))
maillage.draw("k-")
xlim(-10,10)
ylim(-10,10)

                      

fig7.pdf

figure(figsize=(8,8))
maillage.ajouterNoeud(Noeud(9,6))
maillage.draw("k-")
xlim(-10,10)
ylim(-10,10)

                      

fig8.pdf

Pour vérifier que tout fonctionne bien, on ajoute un nœd à l'intérieur d'un des nouveaux triangles :

figure(figsize=(8,8))
maillage.ajouterNoeud(Noeud(5,6))
maillage.draw("k-")
xlim(-10,10)
ylim(-10,10)

                      

fig9.pdf

Un cas particulier à vérifier est celui d'un nouveau nœud situé exactement sur un segment de la frontière. Ajoutons un nœud sur la frontière horizontale (sans erreur d'arrondi possible dans ce cas) :

figure(figsize=(8,8))
maillage.ajouterNoeud(Noeud(-1,-10))
maillage.draw("k-")
xlim(-10,10)
ylim(-10,10)

                      

fig10.pdf

Finalement, l'algorithme permet de générer une triangulation de Delaunay d'un ensemble de points aléatoires quelconque, à condition de partir d'un triangle initial, qui sont évidemment les trois premiers nœuds :

 
maillage = Maillage3()
n1 = Noeud(-10,-10)
n2 = Noeud(10,-10)
n3 = Noeud(0,10)
maillage.triangleInitial(n1,n2,n3)
N = 10
x = uniform(-10,10,N)
y = uniform(-10,10,N)
noeuds = [n1,n2,n3]
for i in range(N):
    n = Noeud(x[i],y[i])
    maillage.ajouterNoeud(n)
    noeuds.append(n)
figure(figsize=(8,8))
maillage.draw("k-")
for n in noeuds:
    n.draw("ro")   
invalid = maillage.verifierDelaunay()
            

fig11.pdf

print(invalid)
--> 0

C'est vraisemblablement l'algorithme implémenté par la fonction scipy.spatial.Delaunay mais nous disposons ici d'un algorithme qui permet de construire un maillage de Delaunay nœud par nœud.

3.c. Insertion aléatoire

L'efficacité de la génération des maillages est une question importante car les maillages pour la méthode des éléments finis peuvent comporter plusieurs centaines de milliers de nœuds.

Outre l'étape de localisation du premier triangle dont le disque ouvert circonscrit contient le nouveau nœud, la destruction et l'ajout de triangles peut être très inefficace lorsque leur nombre est trop grand. Cela se produit pour l'ajout de nœuds externes lorsqu'il sont répartis sur une grille (ou à peu près). Voici un exemple :

 
N = 20
x = np.linspace(-10,10,N)
y = np.linspace(-10,10,N)
maillage = Maillage3()
noeuds = []
n1 = Noeud(x[0],y[0])
n2 = Noeud(x[1],y[0])
n3 = Noeud(x[0],y[1])
maillage.triangleInitial(n1,n2,n3)
for j in range(N):
    for i in range(N):
        if not((i,j) in [(0,0),(1,0),(0,1)]): 
            n = Noeud(x[i],y[j])
            noeuds.append(n)
            

Voici un état intermédiaire du maillage :

P = 4*N+3
nNouv = 0
for k in range(P):
    nNouv += maillage.ajouterNoeud(noeuds[k])
figure(figsize=(8,8))
maillage.draw("k-")
xlim(-10,10)
ylim(-10,10)
nNouv /= P
             

fig12.pdf

print(nNouv)
--> 10.795180722891565

Le nombre de nouveaux triangles créés à chaque insertion est très grand (le nombre de triangles enlevés est quasiment le même). Dans le cas présent, il est en moyenne égale au nombre de nœuds par rangée divisé par 2.

Une solution simple pour augmenter l'efficacité est d'ajouter les nœuds dans un ordre aléatoire. Bien évidemment, cela n'est possible que pour les nœuds dont les positions sont connues dès le début. Nous ajoutons à la classe Maillage une fonction prenant en argument une liste de nœuds à ajouter et qui les ajoute après avoir effectué des permuations aléatoires (autant que le nombre de nœuds de la liste).

class Maillage4(Maillage3):
    def __init__(self):
        Maillage3.__init__(self)
    def ajouterListeNoeuds(self,noeuds,P):
        # P : nombre de noeuds ajoutés 
        N = len(noeuds)
        for p in range(N):
            i = random.randint(0,N-1)
            j = random.randint(0,N-1)
            noeuds[i],noeuds[j] = noeuds[j],noeuds[i]
        nNouv = 0
        for k in range(P):
            nNouv += self.ajouterNoeud(noeuds[k])
        return nNouv/P
             

Voici un état intermédiaire du maillage :

 
N = 20
x = np.linspace(-10,10,N)
y = np.linspace(-10,10,N)
maillage = Maillage4()
noeuds = []
n1 = Noeud(x[0],y[0])
n2 = Noeud(x[1],y[0])
n3 = Noeud(x[0],y[1])
maillage.triangleInitial(n1,n2,n3)
for j in range(N):
    for i in range(N):
        if not((i,j) in [(0,0),(1,0),(0,1)]): 
            n = Noeud(x[i],y[j])
            noeuds.append(n)
P = 4*N+3
nNouv = maillage.ajouterListeNoeuds(noeuds,P)
figure(figsize=(8,8))
maillage.draw("k-")
xlim(-10,10)
ylim(-10,10)
             

fig13.pdf

print(nNouv)
--> 5.397590361445783

On comprend aisément l'intérêt de l'ajout aléatoire : la plupart des nœuds ajoutés sont internes.

Voici le maillage complet :

 
N = 20
x = np.linspace(-10,10,N)
y = np.linspace(-10,10,N)
maillage = Maillage4()
noeuds = []
n1 = Noeud(x[0],y[0])
n2 = Noeud(x[1],y[0])
n3 = Noeud(x[0],y[1])
maillage.triangleInitial(n1,n2,n3)
for j in range(N):
    for i in range(N):
        if not((i,j) in [(0,0),(1,0),(0,1)]): 
            n = Noeud(x[i],y[j])
            noeuds.append(n)
P = N*N-3
nNouv = maillage.ajouterListeNoeuds(noeuds,P)
figure(figsize=(8,8))
maillage.draw("k-")
xlim(-10,10)
ylim(-10,10)
             

fig14.pdf

print(nNouv)
--> 4.944584382871536

Le nombre moyen de triangles créés à chaque insertion est réduit. On est ici dans un cas où il existe plusieurs triangulations de Delaunay : les cercles circonscrits des triangles non situés près des bords contiennent en effet 4 nœuds (sur le cercle).

Il peut être montré que si les nœuds sont insérés dans un ordre aléatoire, le nombre moyen de triangles créés à chaque insertion est inférieur à 6. Ainsi, si le nombre de nœuds par rangée est beaucoup plus grand que 12, le gain en efficacité est très important. Voici un exemple avec 40 nœuds par rangée :

 
N = 20#N = 40
x = np.linspace(-10,10,N)
y = np.linspace(-10,10,N)
maillage = Maillage4()
noeuds = []
n1 = Noeud(x[0],y[0])
n2 = Noeud(x[1],y[0])
n3 = Noeud(x[0],y[1])
maillage.triangleInitial(n1,n2,n3)
for j in range(N):
    for i in range(N):
        if not((i,j) in [(0,0),(1,0),(0,1)]): 
            n = Noeud(x[i],y[j])
            noeuds.append(n)
P = N*N-3
t1 = time.time()
nNouv = maillage.ajouterListeNoeuds(noeuds,P)   
t1 = time.time()-t1
figure(figsize=(8,8))
maillage.draw("k-")
xlim(-10,10)
ylim(-10,10)
             

fig15.pdf

print(nNouv)
--> 4.9773299748110835

Cette méthode a cependant un inconvénient : il est plus difficile d'accélérer l'étape de recherche du premier triangle car la position de chaque nœud ajouté est aléatoire. Pour accélérer la première étape, il faut localiser le nœud inséré dans la triangulation, c'est-à-dire déterminer efficacement un triangle dont le cercle circonscrit contient ce nœud. Si on exclut cette étape de localisation (si on la suppose O(1)), la complexité temporelle de l'insertion aléatoire est O(N) où N est le nombres de nœuds insérés. Si l'insertion n'est pas aléatoire, la complexité peut être O(N²) dans le pire des cas (nœuds sur une grille).

3.d. Accélération de la recherche du premier triangle

Lorsqu'un nouveau nœud N est ajouté, il faut tout d'abord trouver un triangle dont le disque ouvert circonscrit contienne N. La recherche par parcours systématique de tous les triangles est très inefficace. L'utilisation d'un graphe de conflits permet d'accélérer cette recherche. Il s'agit d'un graphe biparti qui relie les nœuds non encore ajoutés aux triangles de la triangulation (réels ou virtuels). Chaque nœud non encore ajouté (mais connu à l'avance) est relié à un et un seul triangle, avec lequel il est en conflit, c'est-à-dire que son disque ouvert circonscrit contient le nœud.

Figure pleine page

Un nœud est connecté à un seul triangle mais un triangle peut avoir plusieurs nœuds avec lequel il est en conflit. Par exemple le triangle T₂ est en conflit avec les nœuds N₁ et N₃. Lorsqu'un nœud est ajouté à la triangulation, le graphe des conflits donne immédiatement un triangle dont le disque circonscrit contient ce nœud. Par exemple si le nœud N₇ est ajouté, ce triangle est T₆.

Le graphe doit être mis à jour après chaque insertion d'un nouveau nœud N_i. La première mise à jour évidente consiste à enlever ce nœud du graphe. Il faut aussi enlever les triangles qui sont supprimés, ceux dont le disque circonscrit contient N_i et qui sont déterminé à partir du premier par recherche récursive par voisinage, comme expliqué plus haut. Enfin il faut ajouter les nouveaux triangles. Enlever un nœud est une opération très simple, consistant à enlever une seule liaison. Enlever un triangle consiste à enlever les liaisons de tous les nœuds avec lesquels ce triangle est en conflit mais cette opération laisse des nœuds sans conflit. Il faut donc attribuer un nouveau triangle à ces nœuds orphelins. Par ailleurs, de nouveaux triangles sont ajoutés au graphe et ces triangles sont des canditats pour les conflits attribués aux nœuds orphelins.

Pour reprendre la figure ci-dessus, supposons qu'on enlève le nœud du graphe N₆ (c'est-à-dire que ce nœud est ajouté à la triangulation) et que les triangles enlevés soient T₃ (déterminé d'après le graphe), T₄ et T₆. On suppose par ailleurs que 3 triangles sont ajoutés, les triangles T₁₀,T₁₁,T₁₂. L'état du graphe devient :

Figure pleine page

Il faut ensuite attribuer un triangle à chacun des nœuds N₄,N₅,N₇,N₈. Soit N_k l'un de ces nœud. Il était en conflit avec un des triangles de la cavité (qui a été enlevé), de sommets A,B,C. Les points N_i et N_k appartiennent au disque ouvert circonscrit de ce triangle. Considérons la droite N_iN_k: elle coupe un côté de la cavité polyhédrique, qui est un côté d'un nouveau triangle et aussi un côté non partagé des triangles enlevés. Notons T le nouveau triangle construit à partir de ce côté, colorié en gris sur la figure suivante :

Figure pleine page

Si N_k est à l'intérieur de la cavité, alors il est dans le triangle T, donc dans son disque circonscrit. C'est le cas du point N_k sur la figure ci-dessus (à gauche). Le nœud N_k est donc en conflit avec le triangle T. Si le nœud N_k n'est pas à l'intérieur de la cavité (point N_k' sur la figure à droite), alors ce nœuds est aussi en conflit à avec le triangle T (démonstration ?). Nous voyons donc que, pour chaque nœud qui n'a plus de conflit, il est possible d'en trouver un parmi les triangles ajoutés dans la cavité, ce qui permet de mettre à jours de graphe des conflits.

Un nœud peut avoir un triangle virtuel enregistré comme conflit s'il se trouve à droite du segment de frontière associé à ce triangle (le demi-plan est l'équivalent du cercle circonscrit). Lorsqu'une cavité contient des triangles virtuels (c.-à-d. des bords de la triangulation), il peut arriver qu'un nouveau triangle virtuel soit attribué à un nœud orphelin (qui étaient en conflit avec un triangle de la cavité, réel ou virtuel). La figure suivante montre un exemple où la cavité comporte trois triangles virtuels et un triangle réel :

Figure pleine page

Le cercle circonscrit du triangle réel est tracé en pointillé. Les trois triangles virtuels sont représentés par une flèche bleue. Les nouveaux triangles sont représentés à droite, avec en vert les côtés qui partagent le nœud ajouté N_i. Comme précédemment, N_k est un nœud qui était en conflit avec un triangle (figure de gauche) et, puisque ce triangle est supprimé, il faut lui attribuer un nouveau triangle avec lequel il est en conflit. Prenons l'exemple du nœud N_k représenté sur la figure. Avant suppression des triangles, ce nœud avait peut-être comme triangle associé dans le graphe des conflits le triangle virtuel situé en haut (mais cela peut être un autre triangle). Considérons, pour chaque nœud réel du triangle virtuel, la droite passant par un nœud du segment et par un point quelconque P situé dans la triangulation (le même point pour tous les triangles). Ces droites sont représentées en trait pointillé sur la figure ci-dessus. Elles permettent de partager l'extérieur de la triangulation en domaines, chacun étant associé à un trianle virtuel. Le côté virtuel du triangle est ainsi une demi-droite. Si le nœud N_k est dans le domaine d'un triangle virtuel, celui-ci est en conflit avec lui. Conclusion : si aucun triangle réel créé n'est en conflit avec le nœud N_k alors il existe un nouveau nœud virtuel qui l'est. L'attribution dans le graphe des conflits d'un triangle au nœud N_k consiste donc à en chercher un parmi les nouveaux nœuds. Bien qu'il soit possible d'accélérer cette recherche, nous la ferons par un parcours linéaire des nouveaux triangles car le nombre de ceux-ci est en moyenne inférieur à 6.

Pour implémenter l'algorithme, il n'est pas nécessaire de disposer explicitement d'un graphe biparti. Il suffit que chaque nœud dispose d'une référence d'un triangle avec lequel il est en conflit et que chaque triangle dispose d'une liste de références de nœuds qui sont en conflit avec lui. Voici la nouvelle classe Noeud :

class Noeud:
    def __init__(self,x,y):
        self.x = x 
        self.y = y 
        self.triangleConflit = None
        self.triangles = []
    def definirConflit(self,tr):
        self.triangleConflit = tr
    def enleverConflit(self):
        self.triangleConflit = None 
    def draw(self,style):
        plot([self.x],[self.y],style)
                             

Voici la nouvelle classe Triangle, dérivée de la première :

class Triangle2(Triangle):
    def __init__(self,n1,n2,n3):
        Triangle.__init__(self,n1,n2,n3)
        self.noeudsConflit = []
    def ajouterConflit(self,n):
        self.noeudsConflit.append(n)
    def enleverConflit(self,n):
        self.noeudsConflit.remove(n)
    
                             

et la nouvelle classe pour les triangles virtuels :

class TriangleVirtuel2(TriangleVirtuel):
    def __init__(self,n1,n2):
        TriangleVirtuel.__init__(self,n1,n2)
        self.noeudsConflit = []
    def ajouterConflit(self,n):
        self.noeudsConflit.append(n)
    def enleverConflit(self,n):
        self.noeudsConflit.remove(n)
                             

La classe Maillage3 doit être revue entièrement. Dans la fonction ajouterListNoeuds, il faut initialiser le graphe des conflits (cette fonction est appelée après que le premier triangle réel soit créé, avec trois triangles virtuels). Dans la fonction ajouterNoeud, la recherche des triangles en conflit avec le nœud ajouté est beaucoup plus rapide que précédemment puisque le premier triangle en conflit est mémorisé dans le graphe des conflits, c'est-à-dire dans l'object Noeud. les autres triangles sont trouvés par recherche récursive comme précédemment. Une liste de nœuds en conflit avec les triangles de la cavité, qui sont enlevé, est établie (noeudsSansConflit). La création des nouveaux triangles se fait comme précédemment mais, dès qu'un triangle est créé, on parcours la liste noeudsSansConflit pour attribuer éventuellement ce triangle à un des nœuds de cette liste, à condition que l'attribution n'ait pas déjà été faite. Le nombre de triangles créés est en moyenne inférieur à 6 donc cette mise à jour du graphe des conflits se fait en temps constant.

class Maillage5(Maillage):
    def __init__(self):
        Maillage.__init__(self)
    def enleverTriangle(self,tr):
        self.retire.append(tr)
        self.triangles.remove(tr)
        for trv in tr.voisins:
            trv.enleverVoisin(tr)
        for i in range(3):
            if tr.noeuds[i]!=None:
                tr.noeuds[i].triangles.remove(tr)
    def triangleInitial(self,n1,n2,n3):
        tr = Triangle2(n1,n2,n3)
        self.triangles.append(tr)
        bord1 = TriangleVirtuel2(tr.noeuds[0],tr.noeuds[1])
        bord2 = TriangleVirtuel2(tr.noeuds[1],tr.noeuds[2])
        bord3 = TriangleVirtuel2(tr.noeuds[2],tr.noeuds[0])
        bords = [bord1,bord2,bord3]
        for b in bords:
            self.triangles.append(b)
            b.ajouterVoisin(tr)
            tr.ajouterVoisin(b)
        bord1.ajouterVoisin(bord2)
        bord2.ajouterVoisin(bord1)
        bord1.ajouterVoisin(bord3)
        bord3.ajouterVoisin(bord1)
        bord2.ajouterVoisin(bord3)
        bord3.ajouterVoisin(bord2)
        self.P = Noeud((n1.x+n2.x+n3.x)/3,(n1.y+n2.y+n3.y)/3)
        
        
    def ajouterListeNoeuds(self,noeuds,P):
        # P : nombre de noeuds ajoutés
        self.noeuds = noeuds  
        for n in noeuds: #initialisation du graphe des conflits avec les 4 premiers triangles
            for i in range(4):
                if self.triangles[i].dansCercle(n):
                    n.definirConflit(self.triangles[i])
                    self.triangles[i].ajouterConflit(n)
                    break
        N = len(noeuds)
        for p in range(N):
            i = random.randint(0,N-1)
            j = random.randint(0,N-1)
            noeuds[i],noeuds[j] = noeuds[j],noeuds[i]
        nNouv = 0
        for k in range(P):
            nNouv += self.ajouterNoeud(noeuds[k])
        
        return nNouv/P 
    def printConflits(self):
        for i in range(len(self.noeuds)):
            print("noeud ",i," conflit :")
            self.noeuds[i].triangleConflit.print()
        for i in range(len(self.triangles)):
            print("triangle ",i)
            for n in self.triangles[i].noeudsConflit:
                print("noeud : ",n.x,n.y)
            
        
    def ajouterNoeud(self,n):
        def rechercheCavite(n,tr,trCavite):
            for t in tr.voisins:
                if not(t in trCavite) and t.dansCercle(n):
                    trCavite.append(t)
                    rechercheCavite(n,t,trCavite)
                
        trCavite = []
        tr = n.triangleConflit
        # bug : noeud sans conflit ?
        trCavite.append(tr)            
        rechercheCavite(n,tr,trCavite) # recherche récursive des autres triangles
        nouveaux = []
        # mise à jour du graphe des conflits
        noeudsSansConflit = [] # noeuds n'ayant plus de triangle associé dans le graphe des conflits
        for tr in trCavite:
            for no in tr.noeudsConflit:
                if no!=n:
                    noeudsSansConflit.append(no)
                    no.enleverConflit()
        
        for i in range(len(trCavite)):
            tr = trCavite[i]
            self.enleverTriangle(tr)
            
            for c in [[0,1],[0,2],[1,2]]: # recherche du côté non partagé
                na = tr.noeuds[c[0]]
                nb = tr.noeuds[c[1]]
                ajout = True
                for j in range(len(trCavite)):
                    if i!=j:
                        if trCavite[j].contientCote(na,nb):
                            ajout = False
                            break
                if ajout:
                    if nb!=None:
                        nouvTr = Triangle2(n,na,nb) 
                    else :
                        if tr.noeuds[0]==na:
                            nouvTr = TriangleVirtuel2(na,n)
                        else:
                            nouvTr = TriangleVirtuel2(n,na)
                    nouveaux.append(nouvTr)
                    self.triangles.append(nouvTr)
                    
                    # recherche d'un triangle en conflit pour chaque noeud qui n'en a plus
                    # nouvTr est-il en conflit avec l'un de ces noeuds ? 
                    for nsc in noeudsSansConflit:
                        if nouvTr.virtuel==False:
                            if nsc.triangleConflit==None and nouvTr.dansCercle(nsc):
                                nsc.definirConflit(nouvTr)
                                nouvTr.ajouterConflit(nsc)
                                
                                
                        else:
                            if nsc.triangleConflit==None and nouvTr.dansCercle(nsc) :
                                nsc.definirConflit(nouvTr)
                                nouvTr.ajouterConflit(nsc)
                                
                                 
                    # recherche des voisins du nouveau triangle parmi les voisins des triangles enlevés
                    voisins = []
                    for t in trCavite:
                        for trv in t.voisins:
                            if not(trv in trCavite) and not(trv in voisins) and trv.triangleVoisin(nouvTr):
                                nouvTr.ajouterVoisin(trv)
                                trv.ajouterVoisin(nouvTr)
                                voisins.append(trv)
        # les nouveaux triangles sont voisins
        for i in range(len(nouveaux)):
            for j in range(len(nouveaux)):
                if i!=j:
                    nouveaux[i].ajouterVoisin(nouveaux[j])
                    
        
        return len(nouveaux)
                            

 
N = 40
x = np.linspace(0,10,N)
y = np.linspace(0,10,N)
maillage = Maillage5() 
noeuds = []
n1 = Noeud(x[0],y[0])
n2 = Noeud(x[1],y[0])
n3 = Noeud(x[0],y[1])
maillage.triangleInitial(n1,n2,n3)
for j in range(N):
    for i in range(N):
        if not((i,j) in [(0,0),(1,0),(0,1)]): 
            n = Noeud(x[i],y[j])
            noeuds.append(n)
P = N*N-3
figure(figsize=(8,8))
t2 = time.time()
nNouv = maillage.ajouterListeNoeuds(noeuds,P)   
t2 = time.time()-t2
maillage.draw("k-") 
xlim(-1,11) 
ylim(-1,11)
             

fig16.pdf

print(nNouv)
--> 5.028804007514089

Voici les temps de calcul sans l'accélération de la recherche du premier triangle et avec :

print(t1)
--> 0.4074852466583252

print(t2)
--> 0.807161808013916

Voici la génération d'une triangulation avec deux fois plus de nœuds que ci-dessus :

 
N = int(N*np.sqrt(2))
x = np.linspace(0,10,N)
y = np.linspace(0,10,N)
maillage = Maillage5() 
noeuds = []
n1 = Noeud(x[0],y[0])
n2 = Noeud(x[1],y[0])
n3 = Noeud(x[0],y[1])
maillage.triangleInitial(n1,n2,n3)
for j in range(N):
    for i in range(N):
        if not((i,j) in [(0,0),(1,0),(0,1)]): 
            n = Noeud(x[i],y[j])
            noeuds.append(n)
P = N*N-3
figure(figsize=(8,8))
t3 = time.time()
nNouv = maillage.ajouterListeNoeuds(noeuds,P)   
t3 = time.time()-t3
maillage.draw("k-") 
xlim(-1,11) 
ylim(-1,11)
             

fig17.pdf

print(t3)
--> 1.6960170269012451

3.e. Retrait d'un nœud

Lorsque le maillage n'est pas satisfaisant, on peut être amené à supprimer un nœud. Cette suppression a pour conséquence immédiate la suppression de tous les triangles qui ont ce nœud comme sommet. Ceux-ci sont immédiatement disponibles car il sont stockés dans une liste associée à chaque nœud. La suppression de ces triangles conduit à une cavité polyhédrique (pas nécessairement convexe) définie par les nœuds des triangles enlevés autres que le nœud enlevé. Ces nœuds sont rangés afin que le polyhèdre soit parcouru dans le sens direct. La figure suivante représente le nœud N à enlever et les triangles auxquels il appartient, qu'il faut aussi enlever.

Figure pleine page

On part du premier triangle rencontré dans la liste des triangles du nœud N, le triangle T₁ sur la figure. Dans ce triangle, on considère le nœud N₁ qui suit N et le nœud suivant N₂. Les nœuds N₁ et N₂ constituent les deux premiers nœuds du polyèdre. Parmis les triangles auxquels N appartient, on recherche celui qui contient N₂, le triangle T₂, dans lequel on trouve le nœud N₃. La recherche est ainsi répété itérativement jusqu'à obtenir à nouveau le nœud N₁. On obtient ainsi la liste des nœuds du polyèdre rangés dans le sens direct. Si le nœud N se trouve sur la frontière de la triangulation, il y a exactement deux triangles virtuels auxquels ce nœud appartient (figure ci-dessus à droite). Lorsqu'on parvient à un triangle virtuel (T₄), qui correspond à un segment de bord, on identifie l'indice dans le triangle virtuel du nœuds de ce segment autre que N (nœud N₃). Si cet indice vaut 0, comme c'est le cas pour T₄ alors le triangle suivant est l'autre triangle virtuel (triangle T₅) de la liste des triangles auxquels N appartient. Pour le triangle T₅, cet indice vaut 1 donc le triangle suivant est un triangle réel contenant le nœud N₄.

Il faut ensuite triangulariser la cavité avec l'algorithme de Bowyer-Watson en traitant les nœuds du polyhèdre comme des nouveaux nœuds, bien qu'ils ne soient pas globalement nouveaux puisqu'ils appartiennent à d'autres triangles (éventuellement à des triangles virtuels). Les trois premiers nœuds du polyèdre définissent le premier triangle (à condition qu'ils ne soient pas alignés).

La triangulation par l'algorithme de Bowyer-Watson, à partir de ce premier triangle et par ajout des nœuds restant du polyèdre, permet d'obtenir une triangulation de Delaunay de la cavité. Cependant, la triangulation obtenue est convexe alors que la cavité ne l'est pas toujours. Cette situation est illustée sur la figure suivante, dans le cas où aucun nœud de la cavité n'est sur un bord de la triangulation :

Figure pleine page

Les deux nouveaux triangles coloriés en rouge doivent être éliminés. Il s'agit de triangles qui ont une intersection avec au moins un des triangles de la triangulation après retrait de ceux de la cavité, ou bien de triangles identiques à un de ces triangles. Il faut détecter efficacement ces triangles. Deux triangles ont une intersection si et seulement si un côté de l'un coupe un côté de l'autre (voir annexe). Pour chaque nouveau triangle, il faut faire ce test avec les triangles auxquels les nœud du triangle appartiennent (qui sont mémorisés dans les nœuds). Par ailleurs, les triangles virtuels générés lors de la triangulation de la cavité, représentés par une flèche rouge dur la figure, doivent aussi être élminés (on verra plus loin qu'on les élimine même s'ils sont au bord).

Il faut aussi considérer le cas où le nœud enlevé se trouve sur la frontière :

Figure pleine page

L'intérieur de la triangulation se trouve à droite. Les triangles enlevés sont représentés en trait pointillé. Les triangles virtuels qui définissaient les bords sont aussi enlevé puisque le nœud enlevé en fait partie. Le polyèdre définissant la cavité contient donc les nœuds 1,2,3 et 4. Les nouveaux triangles apparaissant lors de la triangulation de la cavité sont, par exemple, les triangles 1,2,3 et 2,4,3. Le premier, tracé en bleu, doit être conservé. Le second, tracé en rouge, doit être éliminé. On a aussi représenté les nouveaux triangles virtuels qui doivent être éliminés en rouge et ceux qui doivent être conservés en bleu.

Nous ne voyons pas de méthode simple pour décider si un triangle virtuel doit être gardé ou pas. De plus, ces triangles virtuels ne sont pas toujours orientés convenablement. La stratégie la plus simple est donc de supprimer tous les triangles virtuels générés par la triangulation de la cavité puis de générer des triangles virtuels pour les bords.

Le premier triangle est similaire au premier triangle introduit précédemment pour l'algorithme de Bowyer-Watson, avec un triangle virtuel associé à chacun de ses trois côtés. Il reste à introduire les autres nœuds (nombre de nœuds du polyèdre moins 3) en suivant l'algorithme exposé ci-dessus et implémenté dans la classe Maillage5. La cavité est donc triangularisé indépendamment du reste de la triangulation. Lorsqu'un nœud du polyèdre est ajouté, il n'appartient évidemment à aucun disque circonscrit des triangles situés hors de la cavité. Grace à l'algorithme de Bowyer-Watson, il n'appartient à aucun disque circonscrit des nouveaux triangles. De plus, aucun nœud à l'extérieur de la cavité n'est dans le disque ouvert circonscrit d'un des nouveaux triangles (démonstration ?). En conséquence, la triangulation de la cavité conduit bien à une triangulation globalement Delaunay.

L'implémentation de la triangulation de la cavité est immédiate : il suffit de créer une instance de la classe Maillage5 et de l'utiliser pour générer les triangles, à la suite de quoi ceux-ci sont récupérés et insérés dans la triangulation globale s'ils sont valides. Il faut faire attention au fait que les nœuds du polyèdre sont modifiés par la triangulation de la cavité : il faut donc enlever les triangles invalides de la triangulation de la cavité. Les potentiels voisins des nouveaux triangles sont parmi ces nouveaux triangles et parmi ceux qui étaient voisins des triangles enlevés.

Les triangles virtuels générés par la triangulation de la cavité permettent d'obtenir les nœuds des bords de cette cavité. Lorsqu'un segment de bord appartient à un et un seul triangle, un nouveau triangle virtuel doit être généré pour ce bord, car il s'agit d'un nouveau bord de la triangulation.

Lorsque le nœud est enlevé, les triangles contenant ce nœud sont enlevés et une liste de leurs voisins est générée. Après création des nouveaux triangles, les voisins de ceux-ci sont recherchés dans cette liste.

La classe Maillage6, dérivée de la classe Maillage5, comporte en plus la fonction pour enlever un nœud. Elle comporte aussi une fonction pour détecter l'intersection de deux triangles, bien que celle-ci devrait plutôt être placée dans la classe Triangle.

class Maillage6(Maillage5):
    def __init__(self):
        Maillage5.__init__(self)
    def intersectionSegments(self,A,B,C,D):
        det = (B.y-A.y)*(C.x-D.x)-(D.y-C.y)*(A.x-B.x)
        if det==0: return False
        x = (((B.y-A.y)*A.x+(A.x-B.x)*A.y)*(C.x-D.x)-((D.y-C.y)*C.x+(C.x-D.x)*C.y)*(A.x-B.x))/det
        y = (((D.y-C.y)*C.x+(C.x-D.x)*C.y)*(B.y-A.y)-((B.y-A.y)*A.x+(A.x-B.x)*A.y)*(D.y-C.y))/det
        test1 = (x-A.x)*(B.x-A.x)+(y-A.y)*(B.y-A.y) > 0
        test2 = (x-B.x)*(A.x-B.x)+(y-B.y)*(A.y-B.y) > 0
        test3 = (x-C.x)*(D.x-C.x)+(y-C.y)*(D.y-C.y) > 0
        test4 = (x-D.x)*(C.x-D.x)+(y-D.y)*(C.y-D.y) > 0
        test = test1 and test2 and test3 and test4
        return test
    def intersectionTriangles(self,tr1,tr2):
        if tr1.virtuel or tr2.virtuel:
            return False
        if tr1==tr2: return False
        paires = [(0,1),(0,2),(1,2)]      
        if (tr1.noeuds[0] in tr2.noeuds) and (tr1.noeuds[1] in tr2.noeuds) and (tr1.noeuds[2] in tr2.noeuds): # triangles identiques
            return True
         
        for p in paires:
            for q in paires:
                A = tr1.noeuds[p[0]]
                B = tr1.noeuds[p[1]]
                C = tr2.noeuds[q[0]]
                D = tr2.noeuds[q[1]]
                if A!=C and A!=D and B!=C and B!=D and self.intersectionSegments(A,B,C,D): return True
        return False
    def enleverNoeud(self,n):
        voisins = [] # voisins des triangles enlevés
        retrait = [] # triangles à enlever
        for tr in n.triangles:
            for trv in tr.voisins:
                if not(n in trv.noeuds) and not(trv in voisins):
                    voisins.append(trv)  
            retrait.append(tr)   
                   
        
        # recherche des noeuds du polyèdre dans le sens direct
        polyedre = []
        recherche = True
        premier = None
        k = 0
        while recherche:
            i = 0
            while tr.noeuds[i]!=n: i = (i+1)%3
            # i : indice du noeud dans le triangle 
            k += 1
            if not(tr.virtuel):
                i1 = (i+1)%3
                i2 = (i1+1)%3
                n1 = tr.noeuds[i1]
                if n1==premier:
                    recherche = False 
                    break
                if premier==None : premier = n1
                n2 = tr.noeuds[i2]
                polyedre.append(n1)
                suivant = False
                for t in n.triangles: # recherche du triangle suivant
                    if t!=tr and n2 in t.noeuds: 
                        tr = t
                        suivant = True 
                        break
            else:
                i1 = abs(i-1) 
                n1 = tr.noeuds[i1] # l'autre noeud réel
                if i1==1:
                    for t in n.triangles:  
                        if t!=tr and not(t.virtuel) and n1 in t.noeuds: 
                            tr = t
                            break
                else:
                    polyedre.append(n1)   
                    for t in n.triangles: 
                        if t!=tr and t.virtuel: 
                            tr = t  
                            break
        N = len(polyedre)
        
        
        for tr in retrait:
            self.enleverTriangle(tr) 
              
        cavite = Maillage6()
        cavite.triangleInitial(polyedre[0],polyedre[1],polyedre[2])
        try:
            if len(polyedre) > 3:
                cavite.ajouterListeNoeuds(polyedre[3::],len(polyedre)-3)
        except:
            print("échec de triangulation de la cavité")  
            return
        #cavite.draw("g-",linewidth=3)
         
        # recherche des triangles générés à conserver
        insert = []
        retraitCavite = []
        bordsCavite = [] 
        for tr in cavite.triangles:
            
            insertion = True
            if tr.virtuel: 
                insertion=False
                bordsCavite.append([tr.noeuds[0],tr.noeuds[1]])
            else:
                for tv in voisins: 
                    if self.intersectionTriangles(tv,tr): 
                        insertion = False  
                        break  
            if insertion:  
                insert.append(tr) 
            else: 
                retraitCavite.append(tr) 
        
        # on enlève les triangles inutiles de la cavité (pour la mise à jour des noeuds)        
        for tr in retraitCavite:
            cavite.enleverTriangle(tr)  
        
        #recherche des bords pour créer les nouveaux triangles virtuels.
        trianglesVirtuels = []
        for seg in bordsCavite:
            n1 = seg[0]
            n2 = seg[1]
            num = 0
            for tr in n1.triangles:
                if  n2 in tr.noeuds:
                    num += 1
            if num==1:  # segment appartenant à un seul triangle
                tv = TriangleVirtuel2(n1,n2)
                #if not(tv.dansCercle(n)):
                #    tv.noeuds[0],tv.noeuds[1] = tv.noeuds[1],tv.noeuds[0]
                trianglesVirtuels.append(tv)
                
                  
        for tr in trianglesVirtuels:  
            insert.append(tr)    
        
        self.nouvTriangles = insert
        self.voisins = voisins     
        self.polyedre = polyedre
        self.trianglesVirtuels = trianglesVirtuels 
        
        # insertion des nouveaux triangles dans la triangulation
        for tr in insert:     
            self.triangles.append(tr)   
            for tv in voisins: 
                if tv.triangleVoisin(tr): 
                    tr.ajouterVoisin(tv)
                    tv.ajouterVoisin(tr)
                  
              
                

Voici un exemple, comportant deux nœuds qui vont être enlevés :

    
maillage = Maillage6()
n1 = Noeud(-10,-10)
n2 = Noeud(10,-10)  
n3 = Noeud(0,10)  
maillage.triangleInitial(n1,n2,n3)
N = 30
x = uniform(-10,10,N)
y = uniform(-10,10,N) 
noeuds = []
for i in range(N):
    n = Noeud(x[i],y[i])
    noeuds.append(n)
n1 = Noeud(0,0) # noeud central à enlever
noeuds.append(n1)    
n2 = Noeud(-11,11) # noeud de bord à enlever
noeuds.append(n2)
maillage.ajouterListeNoeuds(noeuds,len(noeuds))
figure(figsize=(8,8))
maillage.draw("k-")
n1.draw("ro")
n2.draw("ro")
xlim(-12,12)
ylim(-12,12) 

fig18.pdf

Les nœuds à enlever sont marqués en rouge. On enlève le nœud situé au centre :

figure(figsize=(8,8))
noeuds.remove(n1)
n1.draw("ro")
maillage.enleverNoeud(n1)
maillage.draw("k-") 
for tr in maillage.nouvTriangles: 
    tr.draw("r-")

xlim(-12,12)
ylim(-12,12)
invalid = maillage.verifierDelaunay()
            

fig19.pdf

print(invalid)
--> 0

Les nouveaux triangles sont tracés en rouge. On retire le nœud du bord :

figure(figsize=(8,8))
noeuds.remove(n2) 
n2.draw("ro")
maillage.enleverNoeud(n2)
maillage.draw("k-") 
for tr in maillage.nouvTriangles: 
    tr.draw("r-") 
 
xlim(-12,12)  
ylim(-12,12)
invalid = maillage.verifierDelaunay()          
            

fig20.pdf

print(invalid)
--> 0

Parmi les nouveaux triangles, il y a les nouveaux triangles virtuels, correspondant aux nouveaux bords. On enlève quelques nœuds sélectionnés au hasard

 

for k in range(10):
    i = random.randint(0,len(noeuds)-1)
    n = noeuds[i] 
    noeuds.remove(n)
    maillage.enleverNoeud(n)
figure(figsize=(8,8)) 
maillage.draw("k-") 
xlim(-12,12) 
ylim(-12,12)
invalid = maillage.verifierDelaunay()
            

fig21.pdf

print(invalid)
--> 0

4.a. Test d'intersection de deux triangles

Deux triangles ont une intersection si et seulement si un côté de l'un coupe un côté de l'autre. Il faut donc détecter l'intersection entre un segment AB et un segment CD. On commence par rechercher le point d'intersection des droites (AB) et (BC), dont les coordonnées sont solutions du système suivant :

$$\begin{array}{cc}\hfill & (y_B-y_A)x+(x_A-x_B)y=(y_B-y_A)x_A+(x_A-x_B)y_A\hfill \\ \hfill & (y_D-y_C)x+(x_C-x_D)y=(y_D-y_C)x_C+(x_C-x_D)y_C\hfill \end{array}$$

Le déterminant de ce système est :

$$det=(y_B-y_A)(x_C-x_D)-(y_D-y_C)(x_A-x_B)$$

S'il est nul, les deux droites sont parallèles. C'est en particulier le cas lorsqu'on teste deux côtés communs de deux triangles. En conséquence, si le déterminant est nul les deux côtés ne se coupent pas. S'il est non nul, on obtient les coordonnées du point d'intersection avec la règle de Cramer :

$$\begin{array}{cc}\hfill & x=\frac{((y_B-y_A)x_A+(x_A-x_B)y_A)(x_C-x_D)-((y_D-y_C)x_C+(x_C-x_D)y_C)(x_A-x_B)}{det}\hfill \\ \hfill & y=\frac{((y_D-y_C)x_C+(x_C-x_D)y_C)(y_B-y_A)-((y_B-y_A)x_A+(x_A-x_B)y_A)(y_D-y_C)}{det}\hfill \end{array}$$

Les deux segments se coupent si et seulement si le point d'intersection M appartient aux deux segments (ouverts) :

$$\begin{array}{cc}\hfill & \overrightarrow{AM}\cdot \overrightarrow{AB}>0\hfill \\ \hfill & \overrightarrow{BM}\cdot \overrightarrow{BA}>0\hfill \\ \hfill & \overrightarrow{CM}\cdot \overrightarrow{CD}>0\hfill \\ \hfill & \overrightarrow{DM}\cdot \overrightarrow{DC}>0\hfill \end{array}$$

c'est-à-dire :

$$\begin{array}{cc}\hfill & (x-x_A)(x_B-x_A)+(y-y_A)(y_B-y_A)>0\hfill \\ \hfill & (x-x_B)(x_A-x_B)+(y-y_B)(y_A-y_B)>0\hfill \\ \hfill & (x-x_C)(x_D-x_C)+(y-y_C)(y_D-y_C)>0\hfill \\ \hfill & (x-x_D)(x_C-x_D)+(y-y_D)(y_C-y_D)>0\hfill \end{array}$$

Si un des nœuds du segment AB est aussi un nœud du segment DC, par exemple si A=C, le test risque de conclure à une intersection à cause des erreurs d'arrondi, alors qu'il n'y a pas d'intersection dans ce cas (puisqu'on considère les segments ouverts). Il faut donc exclure du test ce cas.