Rectification des fichiers et de la version pour le CTAN
[delaunay.git] / luamesh.lua
1 require "luamesh-polygon"
2 require "luamesh-tex"
3
4 local function shallowCopy(original)
5    local copy = {}
6    for key, value in pairs(original) do
7       copy[key] = value
8    end
9    return copy
10 end
11
12 -- Bowyer and Watson algorithm
13 -- Delaunay meshing
14 function BowyerWatson (listPoints,bbox)
15    local triangulation = {}
16    lgth = #listPoints
17    -- add four points to listPoints to have a bounding box
18    listPoints = buildBoundingBox(listPoints)
19    -- the first triangle
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
24    for i=1,lgth do
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))
32       end
33       -- build the new triangles and add them to triangulation
34       for j=1,#polygon do
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"})
37          else
38             table.insert(triangulation,{polygon[j][1],polygon[j][2],i,type="in"})
39          end
40       end
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)
49    end
50    return triangulation
51 end
52
53
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
58    xmin = 1000
59    ymin = 1000
60    xmax = -1000
61    ymax = -1000
62    for i=1,#listPoints do
63       if (listPoints[i].x < xmin) then
64          xmin = listPoints[i].x
65       end
66       if (listPoints[i].x > xmax) then
67          xmax = listPoints[i].x
68       end
69       if (listPoints[i].y < ymin) then
70          ymin = listPoints[i].y
71       end
72       if (listPoints[i].y > ymax) then
73          ymax = listPoints[i].y
74       end
75    end
76    eps = math.max(math.abs(xmax-xmin),math.abs(ymax-ymin))*0.15
77    xmin = xmin - eps
78    xmax = xmax + eps
79    ymin = ymin - eps
80    ymax = ymax + eps
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"})
86    return listPoints
87 end
88
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
93    xmin = 1000
94    ymin = 1000
95    xmax = -1000
96    ymax = -1000
97    for i=1,#listPoints do
98       if (listPoints[i].x < xmin) then
99          xmin = listPoints[i].x
100       end
101       if (listPoints[i].x > xmax) then
102          xmax = listPoints[i].x
103       end
104       if (listPoints[i].y < ymin) then
105          ymin = listPoints[i].y
106       end
107       if (listPoints[i].y > ymax) then
108          ymax = listPoints[i].y
109       end
110    end
111    eps = math.max(math.abs(xmax-xmin),math.abs(ymax-ymin))*0.15
112    xmin = xmin - eps
113    xmax = xmax + eps
114    ymin = ymin - eps
115    ymax = ymax + eps
116    return xmin, xmax, ymin, ymax
117 end
118
119 function removeBoundingBox(triangulation,lgth)
120    -- build the four bounding box edge
121    point1 = lgth+1
122    point2 = lgth+2
123    point3 = lgth+3
124    point4 = lgth+4
125    -- for all triangle
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])
134       end
135    end
136    return newTriangulation
137 end
138
139
140 function buildBadTriangles(point, triangulation,listPoints)
141    local badTriangles = {}
142    for j=1,#triangulation do -- for all triangles
143       A = listPoints[triangulation[j][1]]
144       B = listPoints[triangulation[j][2]]
145       C = listPoints[triangulation[j][3]]
146       center, radius = circoncircle(A,B,C)
147       CP = Vector(center,point)
148       if(VectorNorm(CP)<radius) then -- the point belongs to the circoncirle
149          table.insert(badTriangles,j)
150       end
151    end
152    return badTriangles
153 end
154
155 -- construction of the cavity composed by the bad triangles around the point to add
156 function buildCavity(badTriangles, triangulation)
157    local polygon = {}
158    for j=1,#badTriangles do -- for all bad triangles
159       ind = badTriangles[j]
160       for k=1,3 do -- for all edges
161          edge = {triangulation[ind][k],triangulation[ind][k%3+1]}
162          edgeBord = false
163          for l = 1,#badTriangles do -- for all badtriangles
164             badInd = badTriangles[l]
165             if(badInd ~= ind) then -- if not the current one
166                edgeBord = edgeBord or edgeInTriangle(edge,triangulation[badInd])
167             end
168          end --
169          -- if the edge does not belong to another bad triangle
170          if(edgeBord == false) then
171             -- insert the edge to the cavity
172             table.insert(polygon,edge)
173          end
174       end --
175    end --
176    return polygon
177 end
178
179 function edgeInTriangle(e,t)
180    in1 = false
181    in2 = false
182    for i=1,3 do
183       if e[1] == t[i] then
184          in1 = true
185       end
186       if e[2] == t[i] then
187          in2 = true
188       end
189    end
190    out = (in1 and in2)
191    return out
192 end
193
194 function pointInTriangle(e,t)
195    in1 = false
196    for i=1,3 do
197       if e == t[i] then
198          in1 = true
199       end
200    end
201    return in1
202 end
203
204
205 function Vector(A,B)
206    local out = {x = B.x - A.x, y = B.y - A.y}
207    return out
208 end
209
210 function VectorNorm(NP)
211    return math.sqrt(NP.x*NP.x +NP.y*NP.y)
212 end
213
214 -- circoncircle
215 function circoncircle(M, N, P)
216    -- Compute center and radius of the circoncircle of the triangle M N P
217
218    -- return : (center [Point],radius [float])
219
220    local MN = Vector(M,N)
221    local NP = Vector(N,P)
222    local PM = Vector(P,M)
223    m = VectorNorm(NP)  -- |NP|
224    n = VectorNorm(PM)  -- |PM|
225    p = VectorNorm(MN)  -- |MN|
226
227    d = (m + n + p) * (-m + n + p) * (m - n + p) * (m + n - p)
228    if d > 0 then
229       rad = m * n * p / math.sqrt(d)
230    else
231       rad = 0
232    end
233    d = -2 * (M.x * NP.y + N.x * PM.y + P.x * MN.y)
234    O = {x=0, y=0}
235    OM = Vector(O, M)
236    ON = Vector(O, N)
237    OP = Vector(O, P)
238    om2 = math.pow(VectorNorm(OM),2)  -- |OM|**2
239    on2 = math.pow(VectorNorm(ON),2)  -- |ON|**2
240    op2 = math.pow(VectorNorm(OP),2)  -- |OP|**2
241    x0 = -(om2 * NP.y + on2 * PM.y + op2 * MN.y) / d
242    y0 = (om2 * NP.x + on2 * PM.x + op2 * MN.x) / d
243    if d == 0 then
244       Out = {nil, nil}
245    else
246       Out = {x=x0, y=y0}
247    end
248    return Out, rad  -- (center [Point], R [float])
249 end
250
251 -- compute the list of the circumcircle of a triangulation
252 function listCircumCenter(listPoints,triangulation)
253    local list = {}
254    for j=1,#triangulation do
255       A = listPoints[triangulation[j][1]]
256       B = listPoints[triangulation[j][2]]
257       C = listPoints[triangulation[j][3]]
258       center, radius = circoncircle(A,B,C)
259       table.insert(list,{x=center.x,y=center.y,r=radius})
260    end
261    return list
262 end
263
264 -- find the three neighbour triangles of T
265 function findNeighbour(T,i,triangulation)
266    -- T : triangle
267    -- i : index of T in triangualation
268    -- triangulation
269
270    list = {}
271    -- define the three edge
272    e1 = {T[1],T[2]}
273    e2 = {T[2],T[3]}
274    e3 = {T[3],T[1]}
275    for j=1,#triangulation do
276       if j~= i then
277          if(edgeInTriangle(e1,triangulation[j])) then
278             table.insert(list,j)
279          end
280          if(edgeInTriangle(e2,triangulation[j])) then
281             table.insert(list,j)
282          end
283          if(edgeInTriangle(e3,triangulation[j])) then
284             table.insert(list,j)
285          end
286       end
287    end
288    return list
289 end
290
291 -- test if edge are the same (reverse)
292 function equalEdge(e1,e2)
293    if(((e1[1] == e2[1]) and (e1[2] == e2[2])) or ((e1[1] == e2[2]) and (e1[2] == e2[1]))) then
294       return true
295    else
296       return false
297    end
298 end
299
300 -- test if the edge belongs to the list
301 function edgeInList(e,listE)
302    output = false
303    for i=1,#listE do
304       if(equalEdge(e,listE[i])) then
305          output = true
306       end
307    end
308    return output
309 end
310
311 -- build the edges of the Voronoi diagram with a given triangulation
312 function buildVoronoi(listPoints, triangulation)
313    local listCircumCircle = listCircumCenter(listPoints, triangulation)
314    local listVoronoi = {}
315    for i=1,#listCircumCircle do
316       listN = findNeighbour(triangulation[i],i,triangulation)
317       for j=1,#listN do
318          edge = {i,listN[j]}
319          if( not edgeInList(edge, listVoronoi)) then
320             table.insert(listVoronoi, edge)
321          end
322       end
323    end
324    return listVoronoi
325 end
326
327 -- build the list of points
328 function buildList(chaine, mode)
329    -- if mode = int : the list is given in the chaine string (x1,y1);(x2,y2);...;(xn,yn)
330    -- if mode = ext : the list is given in a file line by line with space separation
331    local listPoints = {}
332    if mode == "int" then
333       local points = string.explode(chaine, ";")
334       local lgth=#points
335       for i=1,lgth do
336          Sx,Sy=string.match(points[i],"%((.+),(.+)%)")
337          listPoints[i]={x=tonumber(Sx),y=tonumber(Sy)}
338       end
339    elseif mode == "ext" then
340       io.input(chaine) -- open the file
341       text=io.read("*all")
342       lines=string.explode(text,"\n+") -- all the lines
343       tablePoints={}
344       for i=1,#lines do
345          xy=string.explode(lines[i]," +")
346          listPoints[i]={x=tonumber(xy[1]),y=tonumber(xy[2])}
347       end
348    else
349       print("Non existing mode")
350    end
351    return listPoints
352 end
353
354
355 -- function to add points on a polygon to respect
356 -- the size of unit mesh
357 function addPointsPolygon(polygon,h)
358    local newPolygon = shallowCopy(polygon)
359    k=0 -- to follow in the newPolygon
360    for i=1,#polygon do
361       k = k+1
362       ip = (i)%(#polygon)+1
363       dist = math.sqrt(math.pow(polygon[i].x-polygon[ip].x,2) + math.pow(polygon[i].y-polygon[ip].y,2))
364       -- if the distance between two ponits of the polygon is greater than 1.5*h
365       if(dist>=2*h) then
366          n = math.floor(dist/h)
367          step = dist/(n+1)
368          for j=1,n do
369             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}
370             table.insert(newPolygon,k+j,a)
371          end
372          k=k+n
373       end
374    end
375    return newPolygon
376 end
377
378 -- function to build a gridpoints from the bounding box
379 -- with a prescribed
380 function buildGrid(listPoints,h,random)
381    -- listPoints : list of the points of the polygon, ordered
382    -- h : parameter for the grid
383    xmin, xmax, ymin, ymax = BoundingBox(listPoints)
384
385    local grid = rectangleList(xmin,xmax,ymin,ymax,h,random)
386    return grid
387 end
388
389 -- function to build the list of points in the rectangle
390 function rectangleList(xmin,xmax,ymin,ymax,h,random)
391    -- for the random
392    math.randomseed( os.time() )
393    nbrX = math.floor(math.abs(xmax-xmin)/h)
394    nbrY = math.floor(math.abs(ymax-ymin)/h)
395    local listPoints = {}
396    k=1
397    for i=1,(nbrX+1) do
398       for j=1,(nbrY+1) do
399          rd = math.random()
400          if(random=="perturb") then
401             fact = 0.3*h
402             --print(fact)
403          else
404             fact = 0.0
405          end
406          listPoints[k] = {x = xmin+(i-1)*h+rd*fact, y=ymin+(j-1)*h+rd*fact}
407          k=k+1
408       end
409    end
410    return listPoints
411 end
412
413
414 -- function to add points from a grid to the interior of a polygon
415 function addGridPoints(polygon, grid,h)
416    local listPoints = shallowCopy(polygon)
417    k = #polygon
418    for i=1, #grid do
419       --print(grid[i].x,grid[i].y)
420       --print(isInside(polygon,grid[i]))
421       if(isInside(polygon,grid[i],h)) then
422          k=k+1
423          listPoints[k] = grid[i]
424       end
425    end
426    return listPoints
427 end
428
429
430
431 -- function give a real polygon without repeting points
432 function cleanPoly(polygon)
433    local polyNew = {}
434    local polyCopy = shallowCopy(polygon)
435    e1 = polyCopy[1][1]
436    e2 = polyCopy[1][2]
437    table.insert(polyNew, e1)
438    table.insert(polyNew, e2)
439    table.remove(polyCopy,1)
440    j = 2
441    while #polyCopy>1 do
442       i=1
443       find = false
444       while (i<=#polyCopy and find==false) do
445          bool1 = (polyCopy[i][1] == polyNew[j])
446          bool2 = (polyCopy[i][2] == polyNew[j])
447          if(bool1 or bool2) then -- the edge has a common point with polyNew[j]
448             if(not bool1) then
449                table.insert(polyNew, polyCopy[i][1])
450                find = true
451                table.remove(polyCopy,i)
452                j = j+1
453             elseif(not bool2) then
454                table.insert(polyNew, polyCopy[i][2])
455                find = true
456                table.remove(polyCopy,i)
457                j = j+1
458             end
459          end
460          i=i+1
461       end
462    end
463    return polyNew
464 end
465
466
467
468 -- build the list of points extern and stop at nbr
469 function buildListExt(chaine, stop)
470    local listPoints = {}
471    io.input(chaine) -- open the file
472    text=io.read("*all")
473    lines=string.explode(text,"\n+") -- all the lines
474    for i=1,tonumber(stop) do
475       xy=string.explode(lines[i]," +")
476       table.insert(listPoints,{x=tonumber(xy[1]),y=tonumber(xy[2])})
477    end
478    xy=string.explode(lines[stop+1]," +")
479    point={x=tonumber(xy[1]),y=tonumber(xy[2])}
480    return point, listPoints
481 end
482
483 function split(pString, pPattern)
484    local Table = {}  -- NOTE: use {n = 0} in Lua-5.0
485    local fpat = "(.-)" .. pPattern
486    local last_end = 1
487    local s, e, cap = pString:find(fpat, 1)
488    while s do
489       if s ~= 1 or cap ~= "" then
490          table.insert(Table,cap)
491       end
492       last_end = e+1
493       s, e, cap = pString:find(fpat, last_end)
494    end
495    if last_end <= #pString then
496       cap = pString:sub(last_end)
497       table.insert(Table, cap)
498    end
499    return Table
500 end
501
502 function readGmsh(file)
503    io.input(file) -- open the file
504    text=io.read("*all")
505    local lines = split(text,"\n+") -- all the lines
506    local listPoints={}
507    local triangulation ={}
508    boolNodes = false
509    Jnodes = 0
510    boolElements = false
511    Jelements = 0
512    J=0
513    for i=1,#lines-J do
514       if(lines[i+J] == "$EndNodes") then
515          boolNodes = false
516          -- go to the next line
517       end
518       if(boolNodes) then -- we are in the Nodes environment
519          xy=split(lines[i+J]," +")
520          table.insert(listPoints,{x=tonumber(xy[2]),y=tonumber(xy[3])})
521       end
522       if(lines[i+J] == "$Nodes") then
523          boolNodes = true
524          -- go to the next line
525          J=J+1
526       end
527       if(lines[i+J] == "$EndElements") then
528          boolElements = false
529          -- go to the next line
530       end
531       if(boolElements) then -- we are in the Nodes environment
532          xy=split(lines[i+J]," +")
533          if(tonumber(xy[2]) == 2) then -- if the element is a triangle
534             nbrTags = xy[3]+1
535             table.insert(triangulation,{tonumber(xy[2+nbrTags+1]),tonumber(xy[2+nbrTags+2]),tonumber(xy[2+nbrTags+3])})
536          end
537       end
538       if(lines[i+J] == "$Elements") then
539          boolElements = true
540          -- go to the next line
541          J=J+1
542       end
543    end
544    return listPoints, triangulation
545 end

Licence Creative Commons Les fichiers de Syracuse sont mis à disposition (sauf mention contraire) selon les termes de la
Licence Creative Commons Attribution - Pas d’Utilisation Commerciale - Partage dans les Mêmes Conditions 4.0 International.