import numpy as np from matplotlib.pyplot import * from scipy.spatial import Delaunay from scipy.linalg import det from numpy.random import uniform import random import time class Noeud: def __init__(self,x,y): self.x = x self.y = y self.triangles = [] def draw(self,style): plot([self.x],[self.y],style) class Triangle: def __init__(self,n1,n2,n3): self.noeuds = [n1,n2,n3] self.detT = n2.x*n3.y-n2.x*n1.y-n1.x*n3.y-n2.y*n3.x+n2.y*n1.x+n1.y*n3.x self.ST = 0.5*abs(self.detT) # aire du triangle if self.detT<0: (n2,n3) = (n3,n2) self.noeuds = [n1,n2,n3] self.detT = n2.x*n3.y-n2.x*n1.y-n1.x*n3.y-n2.y*n3.x+n2.y*n1.x+n1.y*n3.x self.voisins = [] self.virtuel = False # voir plus loin for i in range(3): self.noeuds[i].triangles.append(self) if self.detT==0: print("triangle avec 3 points alignés") def triangleVoisin(self,tr): # test si un triangle est voisin return (self.noeuds[0] in tr.noeuds) or (self.noeuds[1] in tr.noeuds) or (self.noeuds[2] in tr.noeuds) def ajouterVoisin(self,triangle): self.voisins.append(triangle) def enleverVoisin(self,triangle): self.voisins.remove(triangle) def dansCercle(self,n): n1 = self.noeuds[0] n2 = self.noeuds[1] n3 = self.noeuds[2] M = np.array([[n1.x,n1.y,n1.x**2+n1.y**2,1],[n2.x,n2.y,n2.x**2+n2.y**2,1],[n3.x,n3.y,n3.x**2+n3.y**2,1],[n.x,n.y,n.x**2+n.y**2,1]]) m = max(n1.x,n1.y,n2.x,n2.y,n3.x,n3.y,n.x,n.y) return det(M) > 1e-11*m*m def coteSegment(self,n,nA,nB): return (n.y-nA.y)*(nB.x-nA.x)-(n.x-nA.x)*(nB.y-nA.y) >= 0 def contientNoeud(self,n): return self.coteSegment(n,self.noeuds[0],self.noeuds[1]) and self.coteSegment(n,self.noeuds[1],self.noeuds[2]) and self.coteSegment(n,self.noeuds[2],self.noeuds[0]) def cercle(self): x1 = self.noeuds[0].x y1 = self.noeuds[0].y x2 = self.noeuds[1].x y2 = self.noeuds[1].y x3 = self.noeuds[2].x y3 = self.noeuds[2].y det = (x2-x1)*(y3-y1)-(y2-y1)*(x3-x1) xc = ((x2*x2-x1*x1+y2*y2-y1*y1)/2*(y3-y1)-(x3*x3-x1*x1+y3*y3-y1*y1)/2*(y2-y1))/det yc = ((x3*x3-x1*x1+y3*y3-y1*y1)/2*(x2-x1)-(x2*x2-x1*x1+y2*y2-y1*y1)/2*(x3-x1))/det R = np.sqrt((x1-xc)**2+(y1-yc)**2) return xc,yc,R 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==None: plot([self.noeuds[0].x,self.noeuds[1].x,self.noeuds[2].x,self.noeuds[0].x],[self.noeuds[0].y,self.noeuds[1].y,self.noeuds[2].y,self.noeuds[0].y],linewidth=linewidth) else: plot([self.noeuds[0].x,self.noeuds[1].x,self.noeuds[2].x,self.noeuds[0].x],[self.noeuds[0].y,self.noeuds[1].y,self.noeuds[2].y,self.noeuds[0].y],style,linewidth=linewidth) def drawCercle(self,style): xc,yc,R = self.cercle() theta = np.linspace(0,2*np.pi,360) x = xc+R*np.cos(theta) y = yc+R*np.sin(theta) if style!=None: plot(x,y,style) else: plot(x,y) def print(self): print("Triangle : ",self) for n in self.noeuds: print(n.x,n.y) class Maillage: def __init__(self): self.triangles = [] self.retire = [] def fromList(self,noeuds,simplices): triangles = [] for t in simplices: n1 = noeuds[t[0]] n2 = noeuds[t[1]] n3 = noeuds[t[2]] tr = Triangle(n1,n2,n3) triangles.append(tr) for i in range(len(triangles)): #recherche des voisins for j in range(len(triangles)): if i!=j and triangles[i].triangleVoisin(triangles[j]): triangles[i].ajouterVoisin(triangles[j]) self.triangles = triangles def draw(self,style,linewidth=1): for tr in self.triangles: tr.draw(style,linewidth) def drawCercles(self,style): for tr in self.triangles: tr.drawCercle(style) def verifierDelaunay(self): i = 0 test = True for tr in self.triangles: for n in tr.noeuds: for t in self.triangles: if t!=tr: if not(n in t.noeuds) and (n!=None) and t.dansCercle(n) : test=False i += 1 if t.virtuel: t.noeuds[0].draw("ro") t.noeuds[1].draw("ro") else: t.drawCercle("r-") return i 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])) points = np.array(points) tri = Delaunay(points) figure(figsize=(8,8)) maillage = Maillage() maillage.fromList(noeuds,tri.simplices) maillage.draw("k-") xlim(-10,10) ylim(-10,10) invalid = maillage.verifierDelaunay() figure(figsize=(8,8)) maillage.draw("k-") maillage.drawCercles(None) xlim(-10,10) ylim(-10,10) 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]) 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-") 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() 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) 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-") figure(figsize=(8,8)) maillage.ajouterNoeud(Noeud(0,0)) maillage.ajouterNoeud(Noeud(1,-5)) maillage.draw("k-") figure(figsize=(8,8)) maillage.ajouterNoeud(Noeud(5,5)) maillage.draw("k-") xlim(-10,10) ylim(-10,10) figure(figsize=(8,8)) maillage.ajouterNoeud(Noeud(9,6)) maillage.draw("k-") xlim(-10,10) ylim(-10,10) figure(figsize=(8,8)) maillage.ajouterNoeud(Noeud(5,6)) maillage.draw("k-") xlim(-10,10) ylim(-10,10) figure(figsize=(8,8)) maillage.ajouterNoeud(Noeud(-1,-10)) maillage.draw("k-") xlim(-10,10) ylim(-10,10) 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() 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) 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 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 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) 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) 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) 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) 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) 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) 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) 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) 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) 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) 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() 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() 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() 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) 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) 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) 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) 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) 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) 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) 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) 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)