1 require "luamesh-polygon"
4 local function shallowCopy(original)
6 for key, value in pairs(original) do
12 -- Bowyer and Watson algorithm
14 function BowyerWatson (listPoints,bbox)
15 local triangulation = {}
16 local lgth = #listPoints
17 -- add four points to listPoints to have a bounding box
18 listPoints = buildBoundingBox(listPoints)
20 triangulation[1] = {lgth+1, lgth+2, lgth+3,type="bbox"}
21 -- the second triangle
22 triangulation[2] = {lgth+1, lgth+3, lgth+4,type="bbox"}
23 -- add points one by one
25 -- find the triangles which the circumcircle contained the point to add
26 badTriangles = buildBadTriangles(listPoints[i],triangulation,listPoints)
27 -- build the polygon of the cavity containing the point to add
28 polygon = buildCavity(badTriangles, triangulation)
29 -- remove the bad triangles
30 for j=1,#badTriangles do
31 table.remove(triangulation,badTriangles[j]-(j-1))
33 -- build the new triangles and add them to triangulation
35 if((polygon[j][1]>lgth) or (polygon[j][2]>lgth) or (i>lgth)) then
36 table.insert(triangulation,{polygon[j][1],polygon[j][2],i,type="bbox"})
38 table.insert(triangulation,{polygon[j][1],polygon[j][2],i,type="in"})
41 end -- end adding points of the listPoints
42 -- remove bounding box
43 if(bbox ~= "bbox") then
44 triangulation = removeBoundingBox(triangulation,lgth)
45 table.remove(listPoints,lgth+1)
46 table.remove(listPoints,lgth+1)
47 table.remove(listPoints,lgth+1)
48 table.remove(listPoints,lgth+1)
54 function buildBoundingBox(listPoints)
55 -- listPoints : list of points
56 -- epsV : parameter for the distance of the bounding box
57 local xmin, xmax, ymin, ymax, eps
62 for i=1,#listPoints do
63 if (listPoints[i].x < xmin) then
64 xmin = listPoints[i].x
66 if (listPoints[i].x > xmax) then
67 xmax = listPoints[i].x
69 if (listPoints[i].y < ymin) then
70 ymin = listPoints[i].y
72 if (listPoints[i].y > ymax) then
73 ymax = listPoints[i].y
76 eps = math.max(math.abs(xmax-xmin),math.abs(ymax-ymin))*0.15
81 -- add points of the bounding box in last positions
82 table.insert(listPoints,{x=xmin,y=ymin,type="bbox"})
83 table.insert(listPoints,{x=xmin,y=ymax,type="bbox"})
84 table.insert(listPoints,{x=xmax,y=ymax,type="bbox"})
85 table.insert(listPoints,{x=xmax,y=ymin,type="bbox"})
89 function BoundingBox(listPoints)
90 -- listPoints : list of points
91 -- epsV : parameter for the distance of the bounding box
92 local xmin, xmax, ymin, ymax, eps
97 for i=1,#listPoints do
98 if (listPoints[i].x < xmin) then
99 xmin = listPoints[i].x
101 if (listPoints[i].x > xmax) then
102 xmax = listPoints[i].x
104 if (listPoints[i].y < ymin) then
105 ymin = listPoints[i].y
107 if (listPoints[i].y > ymax) then
108 ymax = listPoints[i].y
111 eps = math.max(math.abs(xmax-xmin),math.abs(ymax-ymin))*0.15
116 return xmin, xmax, ymin, ymax
119 function removeBoundingBox(triangulation,lgth)
120 -- build the four bounding box edge
126 local newTriangulation = {}
127 for i=1,#triangulation do
128 boolE1 = pointInTriangle(point1,triangulation[i])
129 boolE2 = pointInTriangle(point2,triangulation[i])
130 boolE3 = pointInTriangle(point3,triangulation[i])
131 boolE4 = pointInTriangle(point4,triangulation[i])
132 if((not boolE1) and (not boolE2) and (not boolE3) and (not boolE4)) then
133 table.insert(newTriangulation,triangulation[i])
136 return newTriangulation
140 function buildBadTriangles(point, triangulation,listPoints)
141 local badTriangles = {}
142 print(#triangulation)
143 for j=1,#triangulation do -- for all triangles
144 A = listPoints[triangulation[j][1]]
145 B = listPoints[triangulation[j][2]]
146 C = listPoints[triangulation[j][3]]
147 center, radius = circoncircle(A,B,C)
148 CP = Vector(center,point)
149 if(VectorNorm(CP)<radius) then -- the point belongs to the circoncirle
150 table.insert(badTriangles,j)
156 -- construction of the cavity composed by the bad triangles around the point to add
157 function buildCavity(badTriangles, triangulation)
159 for j=1,#badTriangles do -- for all bad triangles
160 ind = badTriangles[j]
161 for k=1,3 do -- for all edges
162 edge = {triangulation[ind][k],triangulation[ind][k%3+1]}
164 for l = 1,#badTriangles do -- for all badtriangles
165 badInd = badTriangles[l]
166 if(badInd ~= ind) then -- if not the current one
167 edgeBord = edgeBord or edgeInTriangle(edge,triangulation[badInd])
170 -- if the edge does not belong to another bad triangle
171 if(edgeBord == false) then
172 -- insert the edge to the cavity
173 table.insert(polygon,edge)
180 function edgeInTriangle(e,t)
195 function pointInTriangle(e,t)
207 local out = {x = B.x - A.x, y = B.y - A.y}
211 function VectorNorm(NP)
212 return math.sqrt(NP.x*NP.x +NP.y*NP.y)
216 function circoncircle(M, N, P)
217 -- Compute center and radius of the circoncircle of the triangle M N P
219 -- return : (center [Point],radius [float])
221 local MN = Vector(M,N)
222 local NP = Vector(N,P)
223 local PM = Vector(P,M)
224 m = VectorNorm(NP) -- |NP|
225 n = VectorNorm(PM) -- |PM|
226 p = VectorNorm(MN) -- |MN|
228 d = (m + n + p) * (-m + n + p) * (m - n + p) * (m + n - p)
230 rad = m * n * p / math.sqrt(d)
234 d = -2 * (M.x * NP.y + N.x * PM.y + P.x * MN.y)
239 om2 = math.pow(VectorNorm(OM),2) -- |OM|**2
240 on2 = math.pow(VectorNorm(ON),2) -- |ON|**2
241 op2 = math.pow(VectorNorm(OP),2) -- |OP|**2
242 x0 = -(om2 * NP.y + on2 * PM.y + op2 * MN.y) / d
243 y0 = (om2 * NP.x + on2 * PM.x + op2 * MN.x) / d
249 return Out, rad -- (center [Point], R [float])
252 -- compute the list of the circumcircle of a triangulation
253 function listCircumCenter(listPoints,triangulation)
255 for j=1,#triangulation do
256 A = listPoints[triangulation[j][1]]
257 B = listPoints[triangulation[j][2]]
258 C = listPoints[triangulation[j][3]]
259 center, radius = circoncircle(A,B,C)
260 table.insert(list,{x=center.x,y=center.y,r=radius})
265 -- find the three neighbour triangles of T
266 function findNeighbour(T,i,triangulation)
268 -- i : index of T in triangualation
272 -- define the three edge
276 for j=1,#triangulation do
278 if(edgeInTriangle(e1,triangulation[j])) then
281 if(edgeInTriangle(e2,triangulation[j])) then
284 if(edgeInTriangle(e3,triangulation[j])) then
292 -- test if edge are the same (reverse)
293 function equalEdge(e1,e2)
294 if(((e1[1] == e2[1]) and (e1[2] == e2[2])) or ((e1[1] == e2[2]) and (e1[2] == e2[1]))) then
301 -- test if the edge belongs to the list
302 function edgeInList(e,listE)
305 if(equalEdge(e,listE[i])) then
312 -- build the edges of the Voronoi diagram with a given triangulation
313 function buildVoronoi(listPoints, triangulation)
314 local listCircumCircle = listCircumCenter(listPoints, triangulation)
315 local listVoronoi = {}
316 for i=1,#listCircumCircle do
317 listN = findNeighbour(triangulation[i],i,triangulation)
320 if( not edgeInList(edge, listVoronoi)) then
321 table.insert(listVoronoi, edge)
328 -- build the list of points
329 function buildList(chaine, mode)
330 -- if mode = int : the list is given in the chaine string (x1,y1);(x2,y2);...;(xn,yn)
331 -- if mode = ext : the list is given in a file line by line with space separation
332 local listPoints = {}
333 if mode == "int" then
334 local points = string.explode(chaine, ";")
337 Sx,Sy=string.match(points[i],"%((.+),(.+)%)")
338 listPoints[i]={x=tonumber(Sx),y=tonumber(Sy)}
340 elseif mode == "ext" then
341 io.input(chaine) -- open the file
343 lines=string.explode(text,"\n+") -- all the lines
346 xy=string.explode(lines[i]," +")
347 listPoints[i]={x=tonumber(xy[1]),y=tonumber(xy[2])}
350 print("Non existing mode")
356 -- function to add points on a polygon to respect
357 -- the size of unit mesh
358 function addPointsPolygon(polygon,h)
359 local newPolygon = shallowCopy(polygon)
360 k=0 -- to follow in the newPolygon
363 ip = (i)%(#polygon)+1
364 dist = math.sqrt(math.pow(polygon[i].x-polygon[ip].x,2) + math.pow(polygon[i].y-polygon[ip].y,2))
365 -- if the distance between two ponits of the polygon is greater than 1.5*h
367 n = math.floor(dist/h)
368 print(polygon[i].x,polygon[i].y,polygon[ip].x,polygon[ip].y)
374 a = {x=polygon[i].x+j*step*(polygon[ip].x-polygon[i].x)/dist,y=polygon[i].y+j*step*(polygon[ip].y-polygon[i].y)/dist}
375 print("new = "..a.x.." "..a.y)
376 table.insert(newPolygon,k+j,a)
384 -- function to build a gridpoints from the bounding box
386 function buildGrid(listPoints,h)
387 -- listPoints : list of the points of the polygon, ordered
388 -- h : parameter for the grid
389 xmin, xmax, ymin, ymax = BoundingBox(listPoints)
391 local grid = rectangleList(xmin,xmax,ymin,ymax,h)
395 -- function to build the list of points in the rectangle
396 function rectangleList(xmin,xmax,ymin,ymax,h)
397 nbrX = math.floor(math.abs(xmax-xmin)/h)
398 nbrY = math.floor(math.abs(ymax-ymin)/h)
399 local listPoints = {}
403 listPoints[k] = {x = xmin+(i-1)*h, y=ymin+(j-1)*h}
411 -- function to add points from a grid to the interior of a polygon
412 function addGridPoints(polygon, grid,h)
413 local listPoints = shallowCopy(polygon)
416 --print(grid[i].x,grid[i].y)
417 --print(isInside(polygon,grid[i]))
418 if(isInside(polygon,grid[i],h)) then
420 listPoints[k] = grid[i]
428 -- function give a real polygon without repeting points
429 function cleanPoly(polygon)
431 polyCopy = shallowCopy(polygon)
434 table.insert(polyNew, e1)
435 table.insert(polyNew, e2)
436 table.remove(polyCopy,1)
441 while (i<=#polyCopy and find==false) do
442 bool1 = (polyCopy[i][1] == polyNew[j])
443 bool2 = (polyCopy[i][2] == polyNew[j])
444 if(bool1 or bool2) then -- the edge has a common point with polyNew[j]
446 table.insert(polyNew, polyCopy[i][1])
448 table.remove(polyCopy,i)
450 elseif(not bool2) then
451 table.insert(polyNew, polyCopy[i][2])
453 table.remove(polyCopy,i)
465 -- build the list of points extern and stop at nbr
466 function buildListExt(chaine, stop)
468 io.input(chaine) -- open the file
470 lines=string.explode(text,"\n+") -- all the lines
471 for i=1,tonumber(stop) do
472 xy=string.explode(lines[i]," +")
473 table.insert(listPoints,{x=tonumber(xy[1]),y=tonumber(xy[2])})
475 xy=string.explode(lines[stop+1]," +")
476 point={x=tonumber(xy[1]),y=tonumber(xy[2])}
477 return point, listPoints
480 function split(pString, pPattern)
481 local Table = {} -- NOTE: use {n = 0} in Lua-5.0
482 local fpat = "(.-)" .. pPattern
484 local s, e, cap = pString:find(fpat, 1)
486 if s ~= 1 or cap ~= "" then
487 table.insert(Table,cap)
490 s, e, cap = pString:find(fpat, last_end)
492 if last_end <= #pString then
493 cap = pString:sub(last_end)
494 table.insert(Table, cap)
499 function readGmsh(file)
500 io.input(file) -- open the file
502 local lines = split(text,"\n+") -- all the lines
511 if(lines[i+J] == "$EndNodes") then
513 -- go to the next line
515 if(boolNodes) then -- we are in the Nodes environment
516 xy=split(lines[i+J]," +")
517 table.insert(listPoints,{x=tonumber(xy[2]),y=tonumber(xy[3])})
519 if(lines[i+J] == "$Nodes") then
521 -- go to the next line
524 if(lines[i+J] == "$EndElements") then
526 -- go to the next line
528 if(boolElements) then -- we are in the Nodes environment
529 xy=split(lines[i+J]," +")
530 if(tonumber(xy[2]) == 2) then -- if the element is a triangle
532 table.insert(triangulation,{tonumber(xy[2+nbrTags+1]),tonumber(xy[2+nbrTags+2]),tonumber(xy[2+nbrTags+3])})
535 if(lines[i+J] == "$Elements") then
537 -- go to the next line
541 return listPoints, triangulation