Doc: Avancées, apparence des commandes
[delaunay.git] / luamesh.lua
1 -- Bowyer and Watson algorithm
2 -- Delaunay meshing
3 function BowyerWatson (listPoints,bbox)
4    local triangulation = {}
5    local lgth = #listPoints
6    -- add four points to listPoints to have a bounding box
7    listPoints = buildBoundingBox(listPoints)
8    -- the first triangle
9    triangulation[1] = {lgth+1, lgth+2, lgth+3,type="bbox"}
10    -- the second triangle
11    triangulation[2] = {lgth+1, lgth+3, lgth+4,type="bbox"}
12    -- add points one by one
13    for i=1,lgth do
14       -- find the triangles which the circumcircle contained the point to add
15       badTriangles = buildBadTriangles(listPoints[i],triangulation)
16       -- build the polygon of the cavity containing the point to add
17       polygon = buildCavity(badTriangles, triangulation)
18       -- remove the bad triangles
19       for j=1,#badTriangles do
20          table.remove(triangulation,badTriangles[j]-(j-1))
21       end
22       -- build the new triangles and add them to triangulation
23       for j=1,#polygon do
24          if((polygon[j][1]>lgth) or (polygon[j][2]>lgth) or (i>lgth)) then
25             table.insert(triangulation,{polygon[j][1],polygon[j][2],i,type="bbox"})
26          else
27             table.insert(triangulation,{polygon[j][1],polygon[j][2],i,type="in"})
28          end
29       end
30    end  -- end adding points of the listPoints
31    -- remove bounding box
32    if(bbox ~= "bbox") then
33       triangulation = removeBoundingBox(triangulation,lgth)
34       table.remove(listPoints,lgth+1)
35       table.remove(listPoints,lgth+1)
36       table.remove(listPoints,lgth+1)
37       table.remove(listPoints,lgth+1)
38    end
39    return triangulation
40 end
41
42
43 function buildBoundingBox(listPoints)
44    -- listPoints : list of points
45    -- epsV        : parameter for the distance of the bounding box
46    local xmin, xmax, ymin, ymax, eps
47    xmin = 1000
48    ymin = 1000
49    xmax = -1000
50    ymax = -1000
51    for i=1,#listPoints do
52       if (listPoints[i].x < xmin) then
53          xmin = listPoints[i].x
54       end
55       if (listPoints[i].x > xmax) then
56          xmax = listPoints[i].x
57       end
58       if (listPoints[i].y < ymin) then
59          ymin = listPoints[i].y
60       end
61       if (listPoints[i].y > ymax) then
62          ymax = listPoints[i].y
63       end
64    end
65    eps = math.max(math.abs(xmax-xmin),math.abs(ymax-ymin))*0.15
66    xmin = xmin - eps
67    xmax = xmax + eps
68    ymin = ymin - eps
69    ymax = ymax + eps
70    -- add points of the bounding box in last positions
71    table.insert(listPoints,{x=xmin,y=ymin,type="bbox"})
72    table.insert(listPoints,{x=xmin,y=ymax,type="bbox"})
73    table.insert(listPoints,{x=xmax,y=ymax,type="bbox"})
74    table.insert(listPoints,{x=xmax,y=ymin,type="bbox"})
75    return listPoints
76 end
77
78 function removeBoundingBox(triangulation,lgth)
79    -- build the four bounding box edge
80    point1 = lgth+1
81    point2 = lgth+2
82    point3 = lgth+3
83    point4 = lgth+4
84    -- for all triangle
85    newTriangulation = {}
86    for i=1,#triangulation do
87       boolE1 = pointInTriangle(point1,triangulation[i])
88       boolE2 = pointInTriangle(point2,triangulation[i])
89       boolE3 = pointInTriangle(point3,triangulation[i])
90       boolE4 = pointInTriangle(point4,triangulation[i])
91       if((not boolE1) and (not boolE2) and (not boolE3) and (not boolE4)) then
92          table.insert(newTriangulation,triangulation[i])
93       end
94    end
95    return newTriangulation
96 end
97
98
99 function buildBadTriangles(point, triangulation)
100    badTriangles = {}
101    for j=1,#triangulation do -- for all triangles
102       A = listPoints[triangulation[j][1]]
103       B = listPoints[triangulation[j][2]]
104       C = listPoints[triangulation[j][3]]
105       center, radius = circoncircle(A,B,C)
106       CP = Vector(center,point)
107       if(VectorNorm(CP)<radius) then -- the point belongs to the circoncirle
108          table.insert(badTriangles,j)
109       end
110    end
111    return badTriangles
112 end
113
114 -- construction of the cavity composed by the bad triangles around the point to add
115 function buildCavity(badTriangles, triangulation)
116    polygon = {}
117    for j=1,#badTriangles do -- for all bad triangles
118       ind = badTriangles[j]
119       for k=1,3 do -- for all edges
120          edge = {triangulation[ind][k],triangulation[ind][k%3+1]}
121          edgeBord = false
122          for l = 1,#badTriangles do -- for all badtriangles
123             badInd = badTriangles[l]
124             if(badInd ~= ind) then -- if not the current one
125                edgeBord = edgeBord or edgeInTriangle(edge,triangulation[badInd])
126             end
127          end --
128          -- if the edge does not belong to another bad triangle
129          if(edgeBord == false) then
130             -- insert the edge to the cavity
131             table.insert(polygon,edge)
132          end
133       end --
134    end --
135    return polygon
136 end
137
138 function edgeInTriangle(e,t)
139    in1 = false
140    in2 = false
141    for i=1,3 do
142       if e[1] == t[i] then
143          in1 = true
144       end
145       if e[2] == t[i] then
146          in2 = true
147       end
148    end
149    out = (in1 and in2)
150    return out
151 end
152
153 function pointInTriangle(e,t)
154    in1 = false
155    for i=1,3 do
156       if e == t[i] then
157          in1 = true
158       end
159    end
160    return in1
161 end
162
163
164 function Vector(A,B)
165    local out = {x = B.x - A.x, y = B.y - A.y}
166    return out
167 end
168
169 function VectorNorm(NP)
170    return math.sqrt(NP.x*NP.x +NP.y*NP.y)
171 end
172
173 -- circoncircle
174 function circoncircle(M, N, P)
175    -- Compute center and radius of the circoncircle of the triangle M N P
176
177    -- return : (center [Point],radius [float])
178
179    local MN = Vector(M,N)
180    local NP = Vector(N,P)
181    local PM = Vector(P,M)
182    m = VectorNorm(NP)  -- |NP|
183    n = VectorNorm(PM)  -- |PM|
184    p = VectorNorm(MN)  -- |MN|
185
186    d = (m + n + p) * (-m + n + p) * (m - n + p) * (m + n - p)
187    if d > 0 then
188       rad = m * n * p / math.sqrt(d)
189    else
190       rad = 0
191    end
192    d = -2 * (M.x * NP.y + N.x * PM.y + P.x * MN.y)
193    O = {x=0, y=0}
194    OM = Vector(O, M)
195    ON = Vector(O, N)
196    OP = Vector(O, P)
197    om2 = math.pow(VectorNorm(OM),2)  -- |OM|**2
198    on2 = math.pow(VectorNorm(ON),2)  -- |ON|**2
199    op2 = math.pow(VectorNorm(OP),2)  -- |OP|**2
200    x0 = -(om2 * NP.y + on2 * PM.y + op2 * MN.y) / d
201    y0 = (om2 * NP.x + on2 * PM.x + op2 * MN.x) / d
202    if d == 0 then
203       Out = {nil, nil}
204    else
205       Out = {x=x0, y=y0}
206    end
207    return Out, rad  -- (center [Point], R [float])
208 end
209
210
211 -------------------------- TeX
212 -- build the list of points
213 function buildList(chaine, mode)
214    -- if mode = int : the list is given in the chaine string (x1,y1);(x2,y2);...;(xn,yn)
215    -- if mode = ext : the list is given in a file line by line with space separation
216    listPoints = {}
217    if mode == "int" then
218       local points = string.explode(chaine, ";")
219       local lgth=#points
220       for i=1,lgth do
221          Sx,Sy=string.match(points[i],"%((.+),(.+)%)")
222          listPoints[i]={x=tonumber(Sx),y=tonumber(Sy)}
223       end
224    elseif mode == "ext" then
225       io.input(chaine) -- open the file
226       text=io.read("*all")
227       lines=string.explode(text,"\n+") -- all the lines
228       tablePoints={}
229       for i=1,#lines do
230          xy=string.explode(lines[i]," +")
231          listPoints[i]={x=tonumber(xy[1]),y=tonumber(xy[2])}
232       end
233    else
234       print("Non existing mode")
235    end
236    return listPoints
237 end
238
239
240 --
241 function rectangleList(a,b,nbrA,nbrB)
242    stepA = a/nbrA
243    stepB = b/nbrB
244    listPoints = {}
245    k=1
246    for i=1,(nbrA+1) do
247       for j=1,(nbrB+1) do
248          listPoints[k] = {x = (i-1)*stepA, y=(j-1)*stepB}
249          k=k+1
250       end
251    end
252    return listPoints
253 end
254
255 -- trace a triangulation with TikZ
256 function traceMeshTikZ(listPoints, triangulation,points,color,colorBbox)
257    output = ""
258    for i=1,#triangulation do
259       PointI = listPoints[triangulation[i][1]]
260       PointJ = listPoints[triangulation[i][2]]
261       PointK = listPoints[triangulation[i][3]]
262       if(triangulation[i].type == "bbox") then
263          output = output .. "\\draw[color="..colorBbox.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
264       else
265          output = output .. "\\draw[color="..color.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
266       end
267    end
268    if(points=="points") then
269       for i=1,#listPoints do
270          if(listPoints[i].type == "bbox") then
271             output = output .. "\\draw[color="..colorBbox.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
272          else
273             output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
274          end
275       end
276    end
277    return output
278 end
279
280
281 -- trace a triangulation with MP
282 function traceMeshMP(listPoints, triangulation,points)
283    output = "";
284    output = output .. " pair MeshPoints[];"
285    for i=1,#listPoints do
286       output = output .. "MeshPoints[".. i .. "] = (" .. listPoints[i].x .. "," .. listPoints[i].y .. ");"
287    end
288
289    for i=1,#triangulation do
290       PointI = listPoints[triangulation[i][1]]
291       PointJ = listPoints[triangulation[i][2]]
292       PointK = listPoints[triangulation[i][3]]
293       if(triangulation[i].type == "bbox") then
294          output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolorBbox;"
295       else
296          output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolor;"
297       end
298    end
299    if(points=="points") then
300       j=1
301       for i=1,#listPoints do
302          if(listPoints[i].type == "bbox") then
303             output = output .. "dotlabel.llft (btex $\\MeshPoint^{*}_{"..j.."}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolorBbox ;"
304             j=j+1
305          else
306             output = output .. "dotlabel.llft (btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolor ;"
307          end
308       end
309    end
310    return output
311 end
312
313
314 -- buildMesh with MP
315 function buildMeshMPBW(chaine,mode,points,bbox,scale)
316    listPoints = buildList(chaine, mode)
317    triangulation = BowyerWatson(listPoints,bbox)
318    output = traceMeshMP(listPoints, triangulation,points)
319    output = "\\leavevmode\\begin{mplibcode}beginfig(0);u:="..scale.. ";" .. output .."endfig;\\end{mplibcode}"
320    tex.sprint(output)
321 end
322
323 -- buildMesh with MP include code
324 function buildMeshMPBWinc(chaine,beginning, ending,mode,points,bbox,scale)
325    listPoints = buildList(chaine, mode)
326    triangulation = BowyerWatson(listPoints,bbox)
327    output = traceMeshMP(listPoints, triangulation,points)
328    output = "\\leavevmode\\begin{mplibcode}u:="..scale..";"..beginning .. output .. ending .. "\\end{mplibcode}"
329    tex.sprint(output)
330 end
331
332 -- buildMesh with TikZ
333 function buildMeshTikZBW(chaine,mode,points,bbox,scale,color,colorBbox)
334    listPoints = buildList(chaine, mode)
335    triangulation = BowyerWatson(listPoints,bbox)
336    output = traceMeshTikZ(listPoints, triangulation,points,color,colorBbox)
337    output = "\\noindent\\begin{tikzpicture}[x=" .. scale .. ",y=" .. scale .."]" .. output .."\\end{tikzpicture}"
338    tex.sprint(output)
339 end
340
341 -- buildMesh with TikZ
342 function buildMeshTikZBWinc(chaine,beginning, ending,mode,points,bbox,scale,color,colorBbox)
343    listPoints = buildList(chaine, mode)
344    triangulation = BowyerWatson(listPoints,bbox)
345    output = traceMeshTikZ(listPoints, triangulation,points,color,colorBbox)
346    output = "\\noindent\\begin{tikzpicture}[x=" .. scale .. ",y=" .. scale .."]" ..beginning.. output..ending .."\\end{tikzpicture}"
347    tex.sprint(output)
348 end
349
350
351 -- print points of the mesh
352 function tracePointsMP(listPoints,points)
353    output = "";
354    output = output .. " pair MeshPoints[];"
355    for i=1,#listPoints do
356       output = output .. "MeshPoints[".. i .. "] = (" .. listPoints[i].x .. "," .. listPoints[i].y .. ");"
357    end
358    if(points=="points") then
359       j=1
360       for i=1,#listPoints do
361          if(listPoints[i].type == "bbox") then
362             output = output .. "dotlabel.llft (btex $\\MeshPoint^{*}_{" .. j .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolorBbox ;"
363             j=j+1
364          else
365             output = output .. "dotlabel.llft (btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolor ;"
366          end
367       end
368    else
369       for i=1,#listPoints do
370          if(listPoints[i].type == "bbox") then
371             output = output .. "drawdot  (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u  withcolor \\luameshmpcolorBbox withpen pencircle scaled 3;"
372          else
373             output = output .. "drawdot (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u  withcolor \\luameshmpcolor withpen pencircle scaled 3;"
374          end
375       end
376    end
377    return output
378 end
379
380 -- print points of the mesh
381 function tracePointsTikZ(listPoints,points,color,colorBbox)
382    output = "";
383    if(points=="points") then
384       for i=1,#listPoints do
385          if(listPoints[i].type == "bbox") then
386             output = output .. "\\draw[color="..colorBbox.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
387          else
388             output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
389          end
390       end
391    else
392       for i=1,#listPoints do
393          if(listPoints[i].type == "bbox") then
394             output = output .. "\\draw[color="..colorBbox.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} ;"
395          else
396             output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} ;"
397          end
398       end
399    end
400    return output
401 end
402
403 -- print points to mesh
404 function printPointsMP(chaine,mode,points,bbox,scale)
405    listPoints = buildList(chaine, mode)
406    if(bbox == "bbox" ) then
407       listPoints = buildBoundingBox(listPoints)
408    end
409    output = tracePointsMP(listPoints,points)
410    output = "\\leavevmode\\begin{mplibcode}beginfig(0);u:="..scale.. ";" .. output .."endfig;\\end{mplibcode}"
411    tex.sprint(output)
412 end
413
414
415 -- print points to mesh
416 function printPointsMPinc(chaine,beginning, ending, mode,points,bbox,scale)
417    listPoints = buildList(chaine, mode)
418    if(bbox == "bbox" ) then
419       listPoints = buildBoundingBox(listPoints)
420    end
421    output = tracePointsMP(listPoints,points)
422    output = "\\begin{mplibcode}u:="..scale..";"..beginning .. output .. ending .. "\\end{mplibcode}"
423    tex.sprint(output)
424 end
425
426 -- print points to mesh
427 function printPointsTikZ(chaine,mode,points,bbox,scale,color,colorBbox)
428    listPoints = buildList(chaine, mode)
429    if(bbox == "bbox" ) then
430       listPoints = buildBoundingBox(listPoints)
431    end
432    output = tracePointsTikZ(listPoints,points,color,colorBbox)
433    output = "\\noindent\\begin{tikzpicture}[x=" .. scale .. ",y=" .. scale .."]" .. output .."\\end{tikzpicture}"
434    tex.sprint(output)
435 end
436
437
438 -- print points to mesh
439 function printPointsTikZinc(chaine,beginning, ending, mode,points,bbox,scale,color,colorBbox)
440    listPoints = buildList(chaine, mode)
441    if(bbox == "bbox" ) then
442       listPoints = buildBoundingBox(listPoints)
443    end
444    output = tracePointsTikZ(listPoints,points,color,colorBbox)
445    output = "\\noindent\\begin{tikzpicture}[x=" .. scale .. ",y=" .. scale .."]" ..beginning.. output..ending .."\\end{tikzpicture}"
446    tex.sprint(output)
447 end
448
449
450 -- buildMesh
451 function buildRect(largeur,a,b,nbrA, nbrB)
452    listPoints = rectangleList(a,b,nbrA,nbrB)
453    triangulation = BowyerWatson(listPoints,"none")
454    traceTikZ(listPoints, triangulation,largeur,"none")
455 end
456
457 local function shallowCopy(original)
458    local copy = {}
459    for key, value in pairs(original) do
460       copy[key] = value
461    end
462    return copy
463 end
464
465 -- function give a real polygon without repeting points
466 function cleanPoly(polygon)
467    polyNew = {}
468    polyCopy = shallowCopy(polygon)
469    e1 = polyCopy[1][1]
470    e2 = polyCopy[1][2]
471    table.insert(polyNew, e1)
472    table.insert(polyNew, e2)
473    table.remove(polyCopy,1)
474    j = 2
475    while #polyCopy>1 do
476       i=1
477       find = false
478       while (i<=#polyCopy and find==false) do
479          bool1 = (polyCopy[i][1] == polyNew[j])
480          bool2 = (polyCopy[i][2] == polyNew[j])
481          if(bool1 or bool2) then -- the edge has a common point with polyNew[j]
482             if(not bool1) then
483                table.insert(polyNew, polyCopy[i][1])
484                find = true
485                table.remove(polyCopy,i)
486                j = j+1
487             elseif(not bool2) then
488                table.insert(polyNew, polyCopy[i][2])
489                find = true
490                table.remove(polyCopy,i)
491                j = j+1
492             end
493          end
494          i=i+1
495       end
496    end
497    return polyNew
498 end
499
500 --
501 function TeXaddOnePointTikZ(chaine,point,step,color,colorBack, colorNew, colorCircle,colorBbox)
502    output = ""
503    -- build the triangulation
504    triangulation = BowyerWatson(listPoints,"none")
505    badTriangles = buildBadTriangles(P,triangulation)
506    if(step == "badT") then
507       -- draw all triangle
508       for i=1,#triangulation do
509          PointI = listPoints[triangulation[i][1]]
510          PointJ = listPoints[triangulation[i][2]]
511          PointK = listPoints[triangulation[i][3]]
512          if(triangulation[i].type == "bbox") then
513             output = output .. "\\draw[color="..colorBbox.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
514          else
515             output = output .. "\\draw[color="..color.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
516          end
517       end
518       -- draw and fill the bad triangle
519       for i=1,#badTriangles do
520          PointI = listPoints[triangulation[badTriangles[i]][1]]
521          PointJ = listPoints[triangulation[badTriangles[i]][2]]
522          PointK = listPoints[triangulation[badTriangles[i]][3]]
523          output = output .. "\\draw[fill="..colorBack.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
524       end
525       -- draw the circoncircle
526       for i=1,#badTriangles do
527          PointI = listPoints[triangulation[badTriangles[i]][1]]
528          PointJ = listPoints[triangulation[badTriangles[i]][2]]
529          PointK = listPoints[triangulation[badTriangles[i]][3]]
530          center, radius = circoncircle(PointI, PointJ, PointK)
531          output = output .. "\\draw[dashed, color="..colorCircle.."] ("..center.x .. "," .. center.y .. ") circle ("..radius ..");"
532       end
533       -- mark the points
534       for i=1,#listPoints do
535          if(listPoints[i].type == "bbox") then
536             output = output .. "\\draw[color="..colorBbox.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
537          else
538             output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
539          end
540       end
541       -- mark the point to add
542       output = output .. "\\draw[color="..colorNew.."] (" .. P.x ..",".. P.y .. ") node {$\\bullet$} node[anchor=north east] {$\\NewPoint$};"
543    elseif(step == "cavity") then
544       polygon = buildCavity(badTriangles, triangulation)
545       polyNew = cleanPoly(polygon)
546       -- remove the bad triangles
547       for j=1,#badTriangles do
548          table.remove(triangulation,badTriangles[j]-(j-1))
549       end
550       -- draw the triangles
551       for i=1,#triangulation do
552          PointI = listPoints[triangulation[i][1]]
553          PointJ = listPoints[triangulation[i][2]]
554          PointK = listPoints[triangulation[i][3]]
555          if(triangulation[i].type == "bbox") then
556             output = output .. "\\draw[color="..colorBbox.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
557          else
558             output = output .. "\\draw[color="..color.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
559          end
560       end
561       -- fill and draw the cavity
562       path = ""
563       for i=1,#polyNew do
564          PointI = listPoints[polyNew[i]]
565          path = path .. "(".. PointI.x ..",".. PointI.y ..")--"
566       end
567       output = output .. "\\draw[color="..colorNew..",fill ="..colorBack..", thick] " .. path .. "cycle;"
568       -- mark the points of the mesh
569       for i=1,#listPoints do
570          if(listPoints[i].type == "bbox") then
571             output = output .. "\\draw[color="..colorBbox.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
572          else
573             output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
574          end
575       end
576       -- mark the adding point
577       output = output .. "\\draw[color="..colorNew.."] (" .. P.x ..",".. P.y .. ") node {$\\bullet$} node[anchor=north east] {$\\NewPoint$};"
578    elseif(step == "newT") then
579       polygon = buildCavity(badTriangles, triangulation)
580       polyNew = cleanPoly(polygon)
581       -- remove the bad triangles
582       for j=1,#badTriangles do
583          table.remove(triangulation,badTriangles[j]-(j-1))
584       end
585       -- draw the triangle of the triangulation
586       for i=1,#triangulation do
587          PointI = listPoints[triangulation[i][1]]
588          PointJ = listPoints[triangulation[i][2]]
589          PointK = listPoints[triangulation[i][3]]
590          if(triangulation[i].type == "bbox") then
591             output = output .. "\\draw[color="..colorBbox.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
592          else
593             output = output .. "\\draw[color="..color.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
594          end
595       end
596       -- fill and draw the cavity
597       path = ""
598       for i=1,#polyNew do
599          PointI = listPoints[polyNew[i]]
600          path = path .. "(".. PointI.x ..",".. PointI.y ..")--"
601       end
602       output = output .. "\\draw[color="..colorNew..",fill ="..colorBack..", thick] " .. path .. "cycle;"
603       -- draw the new triangles composed by the edges of the polygon and the added point
604       for i=1,#polygon do
605          output = output .. "\\draw[color=TeXCluaMeshNewTikZ, thick]".."(".. listPoints[polygon[i][1]].x .. "," .. listPoints[polygon[i][1]].y .. ") -- (" .. listPoints[polygon[i][2]].x .. "," .. listPoints[polygon[i][2]].y ..");"
606          output = output .. "\\draw[color="..colorNew..", thick]".."(".. listPoints[polygon[i][1]].x .. "," .. listPoints[polygon[i][1]].y .. ") -- (" .. P.x .. "," .. P.y ..");"
607          output = output .. "\\draw[color="..colorNew..", thick]".."(".. listPoints[polygon[i][2]].x .. "," .. listPoints[polygon[i][2]].y .. ") -- (" .. P.x .. "," .. P.y ..");"
608       end
609       -- mark points
610       for i=1,#listPoints do
611          if(listPoints[i].type == "bbox") then
612             output = output .. "\\draw[color="..colorBbox.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
613          else
614             output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
615          end
616       end
617       -- mark the added point
618       output = output .. "\\draw[color="..colorNew.."] (" .. P.x ..",".. P.y .. ") node {$\\bullet$} node[anchor=north east] {$\\NewPoint$};"
619    end
620    return output
621 end
622
623 function TeXaddOnePointMPBW(listPoints,P,step,bbox)
624    output = "";
625    output = output .. "pair MeshPoints[];"
626    -- build the triangulation
627    triangulation = BowyerWatson(listPoints,bbox)
628    badTriangles = buildBadTriangles(P,triangulation)
629    for i=1,#listPoints do
630       output = output .. "MeshPoints[".. i .. "] = (" .. listPoints[i].x .. "," .. listPoints[i].y .. ");"
631    end
632    if(step == "badT") then
633       -- draw all triangle
634       for i=1,#triangulation do
635          PointI = listPoints[triangulation[i][1]]
636          PointJ = listPoints[triangulation[i][2]]
637          PointK = listPoints[triangulation[i][3]]
638          if(triangulation[i].type == "bbox") then
639             output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolorBbox;"
640          else
641             output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolor;"
642          end
643       end
644       -- draw and fill the bad triangle
645       for i=1,#badTriangles do
646          PointI = listPoints[triangulation[badTriangles[i]][1]]
647          PointJ = listPoints[triangulation[badTriangles[i]][2]]
648          PointK = listPoints[triangulation[badTriangles[i]][3]]
649          output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolor;"
650          output = output .. "fill (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolorBack;"
651       end
652       -- draw the circoncircle
653       for i=1,#badTriangles do
654          PointI = listPoints[triangulation[badTriangles[i]][1]]
655          PointJ = listPoints[triangulation[badTriangles[i]][2]]
656          PointK = listPoints[triangulation[badTriangles[i]][3]]
657          center, radius = circoncircle(PointI, PointJ, PointK)
658          output = output .. "draw fullcircle scaled ("..radius .."*2u) shifted ("..center.x .. "*u," .. center.y .. "*u) dashed evenly withcolor \\luameshmpcolorCircle;"
659       end
660       -- mark the points
661       j=1
662       for i=1,#listPoints do
663          if(listPoints[i].type == "bbox") then
664             output = output .. "dotlabel.llft (btex $\\MeshPoint^{*}_{" .. j .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolorBbox ;"
665             j=j+1
666          else
667             output = output .. "dotlabel.llft (btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolor ;"
668          end
669       end
670       -- mark the point to add
671       output = output .. "dotlabel.llft (btex $\\NewPoint$ etex,(" .. P.x ..",".. P.y .. ")*u) withcolor \\luameshmpcolorNew;"
672    elseif(step == "cavity") then
673       polygon = buildCavity(badTriangles, triangulation)
674       polyNew = cleanPoly(polygon)
675       -- remove the bad triangles
676       for j=1,#badTriangles do
677          table.remove(triangulation,badTriangles[j]-(j-1))
678       end
679       -- draw the triangles
680       for i=1,#triangulation do
681          PointI = listPoints[triangulation[i][1]]
682          PointJ = listPoints[triangulation[i][2]]
683          PointK = listPoints[triangulation[i][3]]
684          if(triangulation[i].type == "bbox") then
685             output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolorBbox;"
686          else
687             output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolor;"
688          end
689       end
690       -- fill and draw the cavity
691       path = ""
692       for i=1,#polyNew do
693          PointI = listPoints[polyNew[i]]
694          path = path .. "(".. PointI.x ..",".. PointI.y ..")*u--"
695       end
696       output = output .. "fill " .. path .. "cycle withcolor \\luameshmpcolorBack;"
697       output = output .. "draw " .. path .. "cycle withcolor \\luameshmpcolorNew  withpen pencircle scaled 1pt;"
698       -- mark the points of the mesh
699       j=1
700       for i=1,#listPoints do
701          if(listPoints[i].type == "bbox") then
702             output = output .. "dotlabel.llft (btex $\\MeshPoint^{*}_{" .. j .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolorBbox ;"
703             j=j+1
704          else
705             output = output .. "dotlabel.llft (btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolor ;"
706          end
707       end
708       -- mark the adding point
709       output = output .. "dotlabel.llft (btex $\\NewPoint$ etex,(" .. P.x ..",".. P.y .. ")*u) withcolor \\luameshmpcolorNew ;"
710    elseif(step == "newT") then
711       polygon = buildCavity(badTriangles, triangulation)
712       polyNew = cleanPoly(polygon)
713       -- remove the bad triangles
714       for j=1,#badTriangles do
715          table.remove(triangulation,badTriangles[j]-(j-1))
716       end
717       -- draw the triangle of the triangulation
718       for i=1,#triangulation do
719          PointI = listPoints[triangulation[i][1]]
720          PointJ = listPoints[triangulation[i][2]]
721          PointK = listPoints[triangulation[i][3]]
722          if(triangulation[i].type == "bbox") then
723             output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolorBbox;"
724          else
725             output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolor;"
726          end
727       end
728       -- fill  the cavity
729       path = ""
730       for i=1,#polyNew do
731          PointI = listPoints[polyNew[i]]
732          path = path .. "(".. PointI.x ..",".. PointI.y ..")*u--"
733       end
734       output = output .. "fill " .. path .. "cycle withcolor \\luameshmpcolorBack;"
735       -- draw the new triangles composed by the edges of the polygon and the added point
736       for i=1,#polygon do
737          output = output .. "draw".."(".. listPoints[polygon[i][1]].x .. "," .. listPoints[polygon[i][1]].y .. ")*u -- (" .. listPoints[polygon[i][2]].x .. "," .. listPoints[polygon[i][2]].y ..")*u withcolor \\luameshmpcolorNew  withpen pencircle scaled 1pt;"
738          output = output .. "draw".."(".. listPoints[polygon[i][1]].x .. "," .. listPoints[polygon[i][1]].y .. ")*u -- (" .. P.x .. "," .. P.y ..")*u withcolor \\luameshmpcolorNew withpen pencircle scaled 1pt;"
739          output = output .. "draw".."(".. listPoints[polygon[i][2]].x .. "," .. listPoints[polygon[i][2]].y .. ")*u -- (" .. P.x .. "," .. P.y ..")*u withcolor \\luameshmpcolorNew withpen pencircle scaled 1pt;"
740       end
741       -- mark points
742       j=1
743       for i=1,#listPoints do
744          if(listPoints[i].type == "bbox") then
745             output = output .. "dotlabel.llft (btex $\\MeshPoint^{*}_{" .. j .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolorBbox ;"
746             j=j+1
747          else
748             output = output .. "dotlabel.llft (btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolor ;"
749          end
750       end
751       -- mark the added point
752       output = output .. "dotlabel.llft (btex $\\NewPoint$ etex,(" .. P.x ..",".. P.y .. ")*u) withcolor \\luameshmpcolorNew ;"
753    end
754    return output
755 end
756
757
758 -- build the list of points extern and stop at nbr
759 function buildListExt(chaine, stop)
760    listPoints = {}
761    io.input(chaine) -- open the file
762    text=io.read("*all")
763    lines=string.explode(text,"\n+") -- all the lines
764    for i=1,tonumber(stop) do
765       xy=string.explode(lines[i]," +")
766       table.insert(listPoints,{x=tonumber(xy[1]),y=tonumber(xy[2])})
767    end
768    xy=string.explode(lines[stop+1]," +")
769    point={x=tonumber(xy[1]),y=tonumber(xy[2])}
770    return point, listPoints
771 end
772
773
774 function TeXOnePointTikZBW(chaine,point,step,scale,mode,bbox,color,colorBack,colorNew,colorCircle,colorBbox)
775    if(mode=="int") then
776       Sx,Sy=string.match(point,"%((.+),(.+)%)")
777       P = {x=Sx, y=Sy}
778       listPoints = buildList(chaine, mode)
779    else
780       -- point is a number
781       P, listPoints = buildListExt(chaine,tonumber(point))
782    end
783    output = TeXaddOnePointTikZ(chaine,point,step,color,colorBack,colorNew,colorCircle,colorBbox)
784    output = "\\noindent\\begin{tikzpicture}[x="..scale..",y="..scale.."]".. output .. "\\end{tikzpicture}"
785    tex.sprint(output)
786 end
787
788 function TeXOnePointTikZBWinc(chaine,point,beginning, ending,step,scale,mode,bbox,color,colorBack,colorNew,colorCircle,colorBbox)
789    if(mode=="int") then
790       Sx,Sy=string.match(point,"%((.+),(.+)%)")
791       P = {x=Sx, y=Sy}
792       listPoints = buildList(chaine, mode)
793    else
794       -- point is a number
795       P, listPoints = buildListExt(chaine,tonumber(point))
796    end
797    output = TeXaddOnePointTikZ(chaine,point,step,color,colorBack,colorNew,colorCircle,colorBbox)
798    output = "\\noindent\\begin{tikzpicture}[x="..scale..",y="..scale.."]".. beginning..output ..ending.. "\\end{tikzpicture}"
799    tex.sprint(output)
800 end
801
802 function TeXOnePointMPBW(chaine,point,step,scale,mode,bbox)
803    if(mode=="int") then
804       Sx,Sy=string.match(point,"%((.+),(.+)%)")
805       P = {x=Sx, y=Sy}
806       listPoints = buildList(chaine, mode)
807    else
808       -- point is a number
809       P, listPoints = buildListExt(chaine,tonumber(point))
810    end
811    output = TeXaddOnePointMPBW(listPoints,P,step,bbox)
812    output = "\\leavevmode\\begin{mplibcode}beginfig(0);u:="..scale..";".. output .. "endfig;\\end{mplibcode}"
813    tex.sprint(output)
814 end
815
816 function TeXOnePointMPBWinc(chaine,point,beginning,ending,step,scale,mode,bbox)
817    if(mode=="int") then
818       Sx,Sy=string.match(point,"%((.+),(.+)%)")
819       P = {x=Sx, y=Sy}
820       listPoints = buildList(chaine, mode)
821    else
822       -- point is a number
823       P, listPoints = buildListExt(chaine,tonumber(point))
824    end
825    output = TeXaddOnePointMPBW(listPoints,P,step,bbox)
826    output = "\\begin{mplibcode}u:="..scale..";"..beginning .. output .. ending .. "\\end{mplibcode}"
827    tex.sprint(output)
828 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.