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):
        if not(tr in self.triangles): return
        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):
        if not(tr in self.triangles): return
        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.373493975903615

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.9672544080604535

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.979848866498741

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.

Nous ajoutons aussi une fonction qui permet d'extraire les nœuds du bord à partir des triangles virtuels, ordonnés dans le sens direct.

class Maillage5(Maillage):
    def __init__(self):
        Maillage.__init__(self)
    def enleverTriangle(self,tr):
        if not(tr in self.triangles): return
        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 rechercheTrianglesConflit(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
        trCavite.append(tr)            
        rechercheCavite(n,tr,trCavite) # recherche récursive des autres triangles
        return trCavite
        
    def ajouterNoeud(self,n):
        trCavite = self.rechercheTrianglesConflit(n)
        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)
        
    def polyedreBord(self):
        virtuels = []
        for tr in self.triangles:
            if tr.virtuel:
                triangle = tr
                break
        polyedre = []
        n1 = triangle.noeuds[0]
        premier = n1
        n2 = triangle.noeuds[1]
        polyedre.append(n1)
        recherche = True
        while recherche:
            for tr in triangle.voisins:
                if tr.virtuel and tr.noeuds[0]==n2:
                    triangle = tr
                    n1 = tr.noeuds[0]
                    n2 = tr.noeuds[1]
                    polyedre.append(n1)
                    break
            if n2==premier:
                recherche = False
        return polyedre 
                            

 
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.028177833437696

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

print(t1)
--> 0.5186221599578857

print(t2)
--> 0.9765517711639404

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)
--> 2.2342886924743652

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. Lors du retrait des triangles, il faut aussi retirer les virtuels associés aux côtés de bord de ces triangles puisque ceux-ci sont regénérés..

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 intersectionCotes(self,A,B,C,D,eps=0):
        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) > eps
        test2 = (x-B.x)*(A.x-B.x)+(y-B.y)*(A.y-B.y) > eps
        test3 = (x-C.x)*(D.x-C.x)+(y-C.y)*(D.y-C.y) > eps
        test4 = (x-D.x)*(C.x-D.x)+(y-D.y)*(C.y-D.y) > eps
        test = test1 and test2 and test3 and test4
        if test: return Noeud(x,y) 
        return None
    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.intersectionCotes(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)
        
        
        virtuels =  []
        for tr in retrait:
            for tv in tr.voisins:
                if tv.virtuel and (tv.noeuds[0] in tr.noeuds) and (tv.noeuds[1] in tr.noeuds): # triangle virtuel ajdacent
                    virtuels.append(tv)
            
        for tr in virtuels:
            tr.draw("b-",linewidth=3)
            self.enleverTriangle(tr)
        
        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. Définition

Les maillages utilisées pour la méthode des éléments finis comportent des frontières sur lesquelles des conditions limites sont imposées, en particulier la condition limite de Neumann. Des nœuds doivent évidemment se trouver sur ces frontières et il faut qu'elles soient représentés par des côtés de triangles, c'est-à-dire qu'aucun triangle ne doit traverser une frontière. On appelera segment un côté de la frontière constitué de deux nœuds consécutifs. Chaque segment doit appartenir à un triangle, comme l'illuste la figure suivante, où la frontière est tracée en trait plus épais :

Figure pleine page

Une triangulation de Delaunay contrainte doit respecter ces contraintes sur les segments. L'ensemble d'éléments à triangulariser est donc constitué d'un ensemble de nœuds et d'un ensemble de segments (constitués de nœuds appartenant au premier).

Il est plus simple d'obtenir cette triangulation si l'on accepte que les triangles adjacents aux segments ne sont pas nécessairement Delaunay. En effet, dans le cas contraire, on doit ajouter des nœuds à l'ensemble introduit initialement. On considère donc une triangulation de Delaunay contrainte dans laquelle les triangles adjacents aux segments ne sont pas nécessairement Delaunay (la triangulation n'est donc pas Delaunay). Il s'en suit que les segments ne sont pas nécessairement Delaunay ni localement Delaunay.

Les triangles adjacents aux segments doivent être Delaunay contraints (Delaunay-C), une propriété moins forte que Delaunay. On dit qu'un nœud A est invivisible d'un nœud B si le segment AB coupe un segment de la frontière mais n'en fait pas partie, autrement dit si les deux nœuds sont séparés par la frontière. Le disque ouvert circonscrit d'un triangle Delaunay-C peut contenir un nœud à condition qu'il soit invisible de tous les points de l'intérieur du triangle.

Figure pleine page

Sur la figure ci-dessus, N₁N₂ est un segment de la frontière. Le triangle adjacent T, constitué des nœuds N₁,N₂,N₃ est Delaunay-C : son disque ouvert circonscrit contient le nœud N mais celui-ci est invisible depuis tout point situé dans le triangle T. Le côté N₁N₂ n'est pas localement Delaunay : il n'a pas de disque ouvert circonscrit ne contenant aucun nœud des deux triangles adjacents.

4.b. Ajout d'un nœud

Considérons l'ajout d'un nœud N dans un triangulation Delaunay-C, qui peut comporter déjà des segments de frontière. Dans l'algorithme de Bowyer-Watson appliqué à une triangulation Delaunay, on élimine tous les triangles qui ne sont plus Delaunay, c'est-à-dire dont le disque ouvert circonscrit contient N. Dans le cas présent, on doit éliminer seulement les triangles qui ne sont plus Delaunay-C. Cela est obtenu avec une variante de l'algorithme de Bowyer-Watson :

• (1) Recherche d'un triangle contenant N.
• (2) Recherche récursive de tous les triangles dont le disque ouvert ouvert circonscrit contient N mais sans franchir un segment de frontière. 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é.

Les étapes (1) et (2) sont modifiées. Pour l'étape (1), on recherche un triangle contenant N, pas seulement un triangle dont le disque ouvert circonscrit contient N. Dans l'étape (2), la recherche récursive s'arrête lorsqu'on rencontre un triangle dont le côté commun avec le précédent est un segment de la frontière.

Figure pleine page

Si les nœuds ajoutés sont préalablement enregistrés dans un graphe de conflits, le premier triangle trouvé a simplement un disque ouvert circonscrit contenant N mais ne contient pas nécessairement N. Supposons que ce premier triangle soit T₁ sur la figure ci-dessus. Une recherche récursive par voisinage permet de trouver le triangle qui contient N (triangle T₂). Si N se trouve exactement sur un côté commun à deux triangles, on obtient l'un des deux triangles auxquels N appartient. À partir de T₂, on recherche récursivement les triangles dont le disque ouvert circonscrit contient N mais sans jamais franchir la frontière. On obtient T₁ et T₄. Lorsqu'on parvient à T₅ à partir de T₄, ces deux triangles ont un côté commun qui appartient à la frontière donc T₅ n'est pas retenu et on ne recherche pas parmi ses voisins. On enlève donc les triangles T₁,T₂,T₄ et on construit 4 nouveaux triangles à partir de N et des côtés de la cavité. Remarquons que le triangle T₅ n'a pas été enlevé bien que son disque ouvert circonscrit contient N. Si on l'avait enlevé, le segment de frontière qu'il contient aurait évidemment été enlevé, ce que précisément on ne veut pas faire. Les nouveaux triangles sont Delaunay-C.

4.c. Ajout d'un segment de frontière

Un segment de frontière peut être ajouté dès que ses deux nœuds ont été ajoutés. Dans un premier temps, on peut par exemple obtenir une triangulation de Delaunay avec tous les nœuds (y compris ceux des frontières) puis ajouter les segments. Après cela, il sera toujours possible d'ajouter d'autres nœuds avec l'algorithme décrit ci-dessus.

La première étape de l'ajout d'un segment consiste à éliminer les triangles coupés par ce segment. La figure suivante montre un exemple avec un segment AB à ajouter :

Figure pleine page

Les triangles T₁,T₂,T₃ sont retirés car il sont coupés par le segment. Le retrait de ces triangles puis l'ajout du segment donne deux cavités C₁ et C₂ que l'on triangularise comme expliqué dans la partie Retrait d'un nœud. Sur cet exemple, il n'y a rien à faire pour la cavité C₂ puisqu'elle ne contient que 3 nœuds. À la différence du remplissage d'une cavité dans une triangulation de Delaunay, les nouveaux triangles peuvent être simplement Delaunay-C (il est évident sur la figure que les triangles créés dans C₁ et C₂ ne sont pas Delaunay mais Delaunay-C).

Il faut détecter efficacement les triangles coupés par le segment AB. Un triangle est coupé par le segment si et seulement si au moins un de ses côtés est coupé par le segment (voir en annexe le test d'intersection de deux triangles). On commence à tester les triangles auquels le nœud A appartient. Si une intersection est trouvée, il faut tester le triangle adjacent à ce dernier et rechercher une intersection sur ses deux autres côtés et ainsi de suite jusqu'à rencontrer un triangle pour lequel il n'y a qu'une intersection (le triangle contenant le nœud B. Il peut arriver qu'une partie du segment soit déjà un côté d'un triangle, comme le montre la figure suivante :

Figure pleine page

Il y a deux triangles, T₁ et T'₁, qui ont deux nœuds situés sur le segment. Ces deux triangles sont retirés. Il faut continuer la recherche sur les triangles auxquels le second nœud appartient.

La recherche de l'intersection d'un triangle avec le segment AB a trois réponses possibles :

• (1) Aucune intersection.
• (2) Une seule intersection sur un sommet.
• (3) Deux intersections sur deux sommets.
• (4) Une intersection sur un côté et une intersection sur un sommet.
• (5) Une intersection sur un côté et une intersection sur un côté.

Le triangle est enlevé dans les cas (3),(4) et (5). Lorsqu'il y a deux intersections, on les range par distance croissante au nœud A. La deuxième intersection permet de savoir comment rechercher le triangle suivant :

• (a) S'il s'agit d'un sommet, on arrête la recherche si ce sommet est le nœud B, sinon on recherche le triangle contenant ce nœud qui a deux intersections avec le segment.
• (b) S'il s'agit d'un côté, on recherche le triangle adjacent, c'est-à-dire le triangle qui a ce côté mais qui n'est pas le premier.

La recherche se fait ainsi de manière itérative en parcourant les triangles dans l'ordre d'intersection par le segment orienté AB. Pour démarrer la recherche, il faut trouver le triangle contenant A qui a deux intersections avec le segment.

Lorsque les triangles copués par le segment sont trouvés, on récupère leurs sommets et on génère les nœuds des deux cavités, constituées respectivement des nœuds situés à gauche et à droite du segment orienté AB. La triangulation de ces cavités se fait comme expliqué dans la partie Retrait d'un nœud. Cette triangularisation générée comporte évidemment un triangle ayant le segment AB pour côté. Contrairement au cas du retrait d'un nœud, le rangement des nœuds d'une cavité dans le sens direct, nécessaire pour traiter les bords, n'est plus trivial.

La classe Maillage6 est complétée pour permettre l'ajout de segments. Ceux-ci sont stockés dans une liste.

La fonction intersectionCoteSegment renvoie l'intersection éventuelle entre un côté de triangle et un segment.

La fonction rechercheTrianglesConflit, qui constitue les étapes 1 et 2 de l'algorithme de Bowyer-Watson, est réécrite afin de pouvoir fonctionner lorsque des nœuds sont ajoutés alors qu'une frontière a déjà été définie.

class Maillage7(Maillage6):
    def __init__(self):
        Maillage6.__init__(self)
        self.segments = []
        
    def rechercheTrianglesConflit(self,n):
        def rechercheTriangleContenant(n,tr,premier,trCavite):
            for t in tr.voisins:
                
                if not(t.virtuel) and t.contientNoeud(n):
                    premier.append(t)
                    return
                elif not(t in trCavite) and t.dansCercle(n):
                    
                    trCavite.append(t)
                    rechercheTriangleContenant(n,t,premier)
                    if len(premier)>0: return
        tr = None
        if not(n.triangleConflit):
            for t in self.triangles:
                if t.dansCercle(n):
                    tr = t
                    break
        else:
            tr = n.triangleConflit
        if tr.dansCercle(n):
            premier = [tr]
        else:
            premier = []
            trCavite = []
            rechercheTriangleContenant(n,tr,premier,trCavite)

        

        def rechercheCavite(n,tr,trCavite):
            for t in tr.voisins:
                frontiere = False
                for seg in self.segments:
                    if not(t.virtuel) and not(tr.virtuel) and t.contientCote(seg[0],seg[1]) and tr.contientCote(seg[0],seg[1]):
                        frontiere=True
                        break
                if not(frontiere) and not(t in trCavite) and t.dansCercle(n):
                    trCavite.append(t)
                    rechercheCavite(n,t,trCavite)
                  
        trCavite = []
        trCavite.append(premier[0])
        rechercheCavite(n,premier[0],trCavite) # recherche récursive des autres triangles
        return trCavite
        
    def alignement(self,nA,nB,n):
        return (n.x-nA.x)*(nB.y-nA.y)-(n.y-nA.y)*(nB.x-nA.x) == 0
    
    def trierIntersect(self,A,intersect):
        n1 = intersect[0][0]
        n2 = intersect[1][0]
        d1 = (n1.x-A.x)**2+(n1.y-A.y)**2
        d2 = (n2.x-A.x)**2+(n2.y-A.y)**2
        if d1 > d2:
            return (intersect[1],intersect[0])
        else:
            return intersect 
     
    def intersectionTriangleSegment(self,tr,nA,nB):
        intersect = []
        if tr.virtuel: return intersect
        #liste d'intersections
        # intersection : (P,[n]), P : point d'intersection, [n] : liste des sommets du triangles (côté ou sommet coupé)
        for n in tr.noeuds:
            if n!=None and self.alignement(nA,nB,n):
                #print("alignement :",n.x,n.y)
                intersect.append((n,[n]))
        
        for p in [(0,1),(1,2),(2,0)]:
            nC = tr.noeuds[p[0]] 
            nD = tr.noeuds[p[1]]
            P = self.intersectionCotes(nA,nB,nC,nD,eps=1e-11)
            if P: intersect.append((P,[nC,nD])) 
        if len(intersect) != 2: return intersect
        return self.trierIntersect(nA,intersect) 
        
    
    def ajouterSegment(self,nA,nB):
        self.segments.append([nA,nB])
        retrait = []
        for tr in nA.triangles:
            intersect = self.intersectionTriangleSegment(tr,nA,nB)
            if len(intersect)==2:   
                retrait.append(tr)
                triangle = tr
                inter = intersect
        recherche = True
        intersect = inter
        
        while recherche:
            noeuds = intersect[1][1] #intersection la plus éloignée de A
            if len(noeuds)==2: # intersection avec un côté du Triangle
                n1 = noeuds[0]
                n2 = noeuds[1] 
                #recherche du triangle adjacent 
                for tr in n1.triangles:
                    if tr!=triangle and tr in n2.triangles:
                        retrait.append(tr) 
                        suivant = tr     
                         
            elif len(noeuds)==1: # intersection avec un sommet 
                n1 = noeuds[0]
                if n1==nB:
                    recherche = False
                    break
                for tr in n1.triangles: # recherche des triangles voisins coupés par le segment
                    if tr!=triangle:
                        intersect = self.intersectionTriangleSegment(tr,nA,nB)
                        if len(intersect)==2:
                            retrait.append(tr)
                            suivant = tr  
            triangle = suivant    
            intersect =  self.intersectionTriangleSegment(triangle,nA,nB)
            if nB in triangle.noeuds: recherche=True
                    
            
        polyedre1 = []
        polyedre2 = []
        
        for tr in retrait:
            for n in tr.noeuds:
                if not(n in polyedre1) and (n.x-nA.x)*(nB.y-nA.y)-(n.y-nA.y)*(nB.x-nA.x) >=0:
                    polyedre1.append(n)
                    #n.draw("bo")
                if not(n in polyedre2) and (n.x-nA.x)*(nB.y-nA.y)-(n.y-nA.y)*(nB.x-nA.x) <=0:
                    polyedre2.append(n)
                    #n.draw("go")
        voisins1 = []
        for n in polyedre1:
            for tr in n.triangles:
                if not(tr in voisins1) and not(tr in retrait):
                    voisins1.append(tr)
        voisins2 = []
        for n in polyedre2:
            for tr in n.triangles:
                if not(tr in voisins2) and not(tr in retrait):
                    voisins2.append(tr)
        for tr in retrait: 
            tr.draw("r-") 
            self.enleverTriangle(tr)
        insert1 = []
        retraitCavite1 = []
        bordsCavite1 = []
        insert2 = []
        retraitCavite2 = []
        bordsCavite2 = [] 
        
        if len(polyedre1) >2:
            cavite1 = Maillage7()
            cavite1.triangleInitial(polyedre1[0],polyedre1[1],polyedre1[2])
            try:
                if (len(polyedre1)>3):
                    cavite1.ajouterListeNoeuds(polyedre1[3::],len(polyedre1)-3)
            except:
                print("échec de triangulation de la cavité 1")  
                return  # recherche des triangles générés à conserver
            insert1 = []
            retraitCavite1 = []
            bordsCavite1 = [] 
            for tr in cavite1.triangles:
                insertion = True
                if tr.virtuel: 
                    insertion=False
                    bordsCavite1.append([tr.noeuds[0],tr.noeuds[1]])
                else:
                    for tv in voisins1: 
                        if self.intersectionTriangles(tv,tr): 
                            insertion = False  
                            break  
                if insertion:  
                    insert1.append(tr) 
                else: 
                    retraitCavite1.append(tr)
        
        if len(polyedre2)>2:
            cavite2 = Maillage7()
            cavite2.triangleInitial(polyedre2[0],polyedre2[1],polyedre2[2])
            try:
                if (len(polyedre2)>3):
                    cavite2.ajouterListeNoeuds(polyedre2[3::],len(polyedre2)-3)
            except:
                print("échec de triangulation de la cavité 2")  
                return
            for tr in cavite2.triangles:
                insertion = True
                if tr.virtuel: 
                    insertion=False
                    bordsCavite2.append([tr.noeuds[0],tr.noeuds[1]])
                else:
                    for tv in voisins2: 
                        if self.intersectionTriangles(tv,tr): 
                            insertion = False  
                            break  
                if insertion:  
                    insert2.append(tr) 
                else: 
                    retraitCavite2.append(tr)
        
        
        # on enlève les triangles inutiles des cavités (pour la mise à jour des noeuds)        
        for tr in retraitCavite1:
            cavite1.enleverTriangle(tr)
        for tr in retraitCavite2:
            cavite2.enleverTriangle(tr)


        # insertion des nouveaux triangles dans la triangulation
        for tr in insert1:     
            self.triangles.append(tr)   
            for tv in voisins1: 
                if tv.triangleVoisin(tr): 
                    tr.ajouterVoisin(tv)
                    tv.ajouterVoisin(tr)
        for tr in insert2:     
            self.triangles.append(tr)   
            for tv in voisins2: 
                if tv.triangleVoisin(tr): 
                    tr.ajouterVoisin(tv)
                    tv.ajouterVoisin(tr)
    
    def drawSegments(self,style="r-",linewidth=1):
        for seg in self.segments:
            plot([seg[0].x,seg[1].x],[seg[0].y,seg[1].y],style,linewidth=linewidth)
              

Voici un exemple, comportant deux nœuds qui vont servir à définir un segment de frontière :

    
maillage = Maillage7() 
n1 = Noeud(-10,-10)
n2 = Noeud(10,-5)  
n3 = Noeud(0,10)  
maillage.triangleInitial(n1,n2,n3)
N = 50
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)
nA = Noeud(-2,-2) 
noeuds.append(nA)    
nB = Noeud(6,5)
noeuds.append(nB)
maillage.ajouterListeNoeuds(noeuds,len(noeuds))
figure(figsize=(8,8))
maillage.draw("k-")
nA.draw("ro")
nB.draw("ro") 
xlim(-12,12)
ylim(-12,12)  

fig22.pdf

figure(figsize=(8,8))
maillage.draw("k-")
maillage.ajouterSegment(nA,nB)
plot([nA.x,nB.x],[nA.y,nB.y],"k--")
#nA.draw("ro") 
#nB.draw("ro")
xlim(-12,12) 
ylim(-12,12)  
        

fig23.pdf

figure(figsize=(8,8))
maillage.draw("k-")
plot([nA.x,nB.x],[nA.y,nB.y],"k-",linewidth=3)
#nA.draw("ro") 
#nB.draw("ro")
xlim(-12,12) 
ylim(-12,12)  
        

fig24.pdf

L'algorithme fonctionne aussi si les deux nœuds du segment sont sur les bords de la triangulation (les triangles virtuels ne sont pas affectés) :

figure(figsize=(8,8))
maillage.draw("k-")
bord = maillage.polyedreBord()
ia = random.randint(0,len(bord)-1)
ib = (ia+2)%len(bord)
nA = bord[ia] 
nB = bord[ib]
maillage.ajouterSegment(nA,nB)
plot([nA.x,nB.x],[nA.y,nB.y],"k--")
xlim(-12,12) 
ylim(-12,12)
        

fig25.pdf

 
figure(figsize=(8,8))
maillage.draw("k-")
plot([nA.x,nB.x],[nA.y,nB.y],"k-",linewidth=3)
#nA.draw("ro") 
#nB.draw("ro")
xlim(-12,12) 
ylim(-12,12)  
        

fig26.pdf

Voici un exemple où une frontière circulaire est générée à l'intérieur d'un maillage carré :

N = 30
x = np.linspace(-10,10,N)
y = np.linspace(-10,10,N)
maillage = Maillage7() 
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
R = 5
Nf = 30
dtheta = np.pi*2/Nf
frontiere = [] 
for i in range(Nf):
    theta = i*dtheta
    n = Noeud(R*np.cos(theta),R*np.sin(theta))
    noeuds.append(n)
    frontiere.append(n)
P += Nf
figure(figsize=(8,8))
maillage.ajouterListeNoeuds(noeuds,P)
maillage.draw("k-") 
xlim(-11,11) 
ylim(-11,11)        
        

fig27.pdf

Dans cette triangulation de Delaunay, les segments de la frontière ne coïncident pas nécessairement avec un côté de triangle. On ajoute les segments :

for i in range(Nf-1):
    maillage.ajouterSegment(frontiere[i],frontiere[i+1])
maillage.ajouterSegment(frontiere[Nf-1],frontiere[0])
figure(figsize=(8,8))
maillage.draw("k-")
maillage.drawSegments("r-") 
xlim(-11,11) 
ylim(-11,11)       
        

fig28.pdf

Voici un test de l'ajout d'un nœud après que la frontière soit créée avec une triangulation Delaunay-C. Lorsque le nœud ajouté n'est pas dans le graphe des conflits, une recherche itérative d'un triangle en conflit est faite dans la fonction rechercheTrianglesConflit.

figure(figsize=(8,8))
n1 = Noeud(R+0.3,0)
maillage.ajouterNoeud(n1)
n2 = Noeud(R-0.3,0)
maillage.ajouterNoeud(n2)
maillage.draw("k-") 
maillage.drawSegments("r-")  
n1.draw("ro")
n2.draw("ro")
xlim(2,8) 
ylim(-3,3)       
        

fig29.pdf

5.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 côtés se coupent si et seulement si le point d'intersection M appartient aux deux côtés (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 côtés ouverts). Il faut donc exclure du test ce cas.