1 -- Bowyer and Watson algorithm
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)
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
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))
22 -- build the new triangles and add them to triangulation
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"})
27 table.insert(triangulation,{polygon[j][1],polygon[j][2],i,type="in"})
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)
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
51 for i=1,#listPoints do
52 if (listPoints[i].x < xmin) then
53 xmin = listPoints[i].x
55 if (listPoints[i].x > xmax) then
56 xmax = listPoints[i].x
58 if (listPoints[i].y < ymin) then
59 ymin = listPoints[i].y
61 if (listPoints[i].y > ymax) then
62 ymax = listPoints[i].y
65 eps = math.max(math.abs(xmax-xmin),math.abs(ymax-ymin))*0.15
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"})
78 function removeBoundingBox(triangulation,lgth)
79 -- build the four bounding box edge
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])
95 return newTriangulation
99 function buildBadTriangles(point, triangulation)
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)
114 -- construction of the cavity composed by the bad triangles around the point to add
115 function buildCavity(badTriangles, triangulation)
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]}
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])
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)
138 function edgeInTriangle(e,t)
153 function pointInTriangle(e,t)
165 local out = {x = B.x - A.x, y = B.y - A.y}
169 function VectorNorm(NP)
170 return math.sqrt(NP.x*NP.x +NP.y*NP.y)
174 function circoncircle(M, N, P)
175 -- Compute center and radius of the circoncircle of the triangle M N P
177 -- return : (center [Point],radius [float])
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|
186 d = (m + n + p) * (-m + n + p) * (m - n + p) * (m + n - p)
188 rad = m * n * p / math.sqrt(d)
192 d = -2 * (M.x * NP.y + N.x * PM.y + P.x * MN.y)
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
207 return Out, rad -- (center [Point], R [float])
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
217 if mode == "int" then
218 local points = string.explode(chaine, ";")
221 Sx,Sy=string.match(points[i],"%((.+),(.+)%)")
222 listPoints[i]={x=tonumber(Sx),y=tonumber(Sy)}
224 elseif mode == "ext" then
225 io.input(chaine) -- open the file
227 lines=string.explode(text,"\n+") -- all the lines
230 xy=string.explode(lines[i]," +")
231 listPoints[i]={x=tonumber(xy[1]),y=tonumber(xy[2])}
234 print("Non existing mode")
241 function rectangleList(a,b,nbrA,nbrB)
248 listPoints[k] = {x = (i-1)*stepA, y=(j-1)*stepB}
255 -- trace a triangulation with TikZ
256 function traceMeshTikZ(listPoints, triangulation,points,color,colorBbox)
258 for i=1,#listPoints do
259 output = output .. "\\coordinate (MeshPoints".. i .. ") at (" .. listPoints[i].x .. "," .. listPoints[i].y .. ");"
261 for i=1,#triangulation do
262 PointI = listPoints[triangulation[i][1]]
263 PointJ = listPoints[triangulation[i][2]]
264 PointK = listPoints[triangulation[i][3]]
265 if(triangulation[i].type == "bbox") then
266 output = output .. "\\draw[color="..colorBbox.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
268 output = output .. "\\draw[color="..color.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
271 if(points=="points") then
273 for i=1,#listPoints do
274 if(listPoints[i].type == "bbox") then
275 output = output .. "\\draw[color="..colorBbox.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint^*_{" .. j .. "}$};"
278 output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
286 -- trace a triangulation with MP
287 function traceMeshMP(listPoints, triangulation,points)
289 output = output .. " pair MeshPoints[];"
290 for i=1,#listPoints do
291 output = output .. "MeshPoints[".. i .. "] = (" .. listPoints[i].x .. "," .. listPoints[i].y .. ")*u;"
294 for i=1,#triangulation do
295 PointI = listPoints[triangulation[i][1]]
296 PointJ = listPoints[triangulation[i][2]]
297 PointK = listPoints[triangulation[i][3]]
298 if(triangulation[i].type == "bbox") then
299 output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolorBbox;"
301 output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolor;"
304 if(points=="points") then
306 for i=1,#listPoints do
307 if(listPoints[i].type == "bbox") then
308 output = output .. "dotlabel.llft (btex $\\MeshPoint^{*}_{"..j.."}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolorBbox ;"
311 output = output .. "dotlabel.llft (btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolor ;"
320 function buildMeshMPBW(chaine,mode,points,bbox,scale)
321 listPoints = buildList(chaine, mode)
322 triangulation = BowyerWatson(listPoints,bbox)
323 output = traceMeshMP(listPoints, triangulation,points)
324 output = "\\leavevmode\\begin{mplibcode}beginfig(0);u:="..scale.. ";" .. output .."endfig;\\end{mplibcode}"
328 -- buildMesh with MP include code
329 function buildMeshMPBWinc(chaine,beginning, ending,mode,points,bbox,scale)
330 listPoints = buildList(chaine, mode)
331 triangulation = BowyerWatson(listPoints,bbox)
332 output = traceMeshMP(listPoints, triangulation,points)
333 output = "\\leavevmode\\begin{mplibcode}u:="..scale..";"..beginning .. output .. ending .. "\\end{mplibcode}"
337 -- buildMesh with TikZ
338 function buildMeshTikZBW(chaine,mode,points,bbox,scale,color,colorBbox)
339 listPoints = buildList(chaine, mode)
340 triangulation = BowyerWatson(listPoints,bbox)
341 output = traceMeshTikZ(listPoints, triangulation,points,color,colorBbox)
342 output = "\\noindent\\begin{tikzpicture}[x=" .. scale .. ",y=" .. scale .."]" .. output .."\\end{tikzpicture}"
346 -- buildMesh with TikZ
347 function buildMeshTikZBWinc(chaine,beginning, ending,mode,points,bbox,scale,color,colorBbox)
348 listPoints = buildList(chaine, mode)
349 triangulation = BowyerWatson(listPoints,bbox)
350 output = traceMeshTikZ(listPoints, triangulation,points,color,colorBbox)
351 output = "\\noindent\\begin{tikzpicture}[x=" .. scale .. ",y=" .. scale .."]" ..beginning.. output..ending .."\\end{tikzpicture}"
356 -- print points of the mesh
357 function tracePointsMP(listPoints,points)
359 output = output .. " pair MeshPoints[];"
360 for i=1,#listPoints do
361 output = output .. "MeshPoints[".. i .. "] = (" .. listPoints[i].x .. "," .. listPoints[i].y .. ")*u;"
363 if(points=="points") then
365 for i=1,#listPoints do
366 if(listPoints[i].type == "bbox") then
367 output = output .. "dotlabel.llft (btex $\\MeshPoint^{*}_{" .. j .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolorBbox ;"
370 output = output .. "dotlabel.llft (btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolor ;"
374 for i=1,#listPoints do
375 if(listPoints[i].type == "bbox") then
376 output = output .. "drawdot (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u withcolor \\luameshmpcolorBbox withpen pencircle scaled 3;"
378 output = output .. "drawdot (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u withcolor \\luameshmpcolor withpen pencircle scaled 3;"
385 -- print points of the mesh
386 function tracePointsTikZ(listPoints,points,color,colorBbox)
388 for i=1,#listPoints do
389 output = output .. "\\coordinate (MeshPoints".. i .. ") at (" .. listPoints[i].x .. "," .. listPoints[i].y .. ");"
391 if(points=="points") then
393 for i=1,#listPoints do
394 if(listPoints[i].type == "bbox") then
395 output = output .. "\\draw[color="..colorBbox.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint^*_{" .. j .. "}$};"
398 output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
402 for i=1,#listPoints do
403 if(listPoints[i].type == "bbox") then
404 output = output .. "\\draw[color="..colorBbox.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} ;"
406 output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} ;"
413 -- print points to mesh
414 function printPointsMP(chaine,mode,points,bbox,scale)
415 listPoints = buildList(chaine, mode)
416 if(bbox == "bbox" ) then
417 listPoints = buildBoundingBox(listPoints)
419 output = tracePointsMP(listPoints,points)
420 output = "\\leavevmode\\begin{mplibcode}beginfig(0);u:="..scale.. ";" .. output .."endfig;\\end{mplibcode}"
425 -- print points to mesh
426 function printPointsMPinc(chaine,beginning, ending, mode,points,bbox,scale)
427 listPoints = buildList(chaine, mode)
428 if(bbox == "bbox" ) then
429 listPoints = buildBoundingBox(listPoints)
431 output = tracePointsMP(listPoints,points)
432 output = "\\leavevmode\\begin{mplibcode}u:="..scale..";"..beginning .. output .. ending .. "\\end{mplibcode}"
436 -- print points to mesh
437 function printPointsTikZ(chaine,mode,points,bbox,scale,color,colorBbox)
438 listPoints = buildList(chaine, mode)
439 if(bbox == "bbox" ) then
440 listPoints = buildBoundingBox(listPoints)
442 output = tracePointsTikZ(listPoints,points,color,colorBbox)
443 output = "\\noindent\\begin{tikzpicture}[x=" .. scale .. ",y=" .. scale .."]" .. output .."\\end{tikzpicture}"
448 -- print points to mesh
449 function printPointsTikZinc(chaine,beginning, ending, mode,points,bbox,scale,color,colorBbox)
450 listPoints = buildList(chaine, mode)
451 if(bbox == "bbox" ) then
452 listPoints = buildBoundingBox(listPoints)
454 output = tracePointsTikZ(listPoints,points,color,colorBbox)
455 output = "\\noindent\\begin{tikzpicture}[x=" .. scale .. ",y=" .. scale .."]" ..beginning.. output..ending .."\\end{tikzpicture}"
461 function buildRect(largeur,a,b,nbrA, nbrB)
462 listPoints = rectangleList(a,b,nbrA,nbrB)
463 triangulation = BowyerWatson(listPoints,"none")
464 traceTikZ(listPoints, triangulation,largeur,"none")
467 local function shallowCopy(original)
469 for key, value in pairs(original) do
475 -- function give a real polygon without repeting points
476 function cleanPoly(polygon)
478 polyCopy = shallowCopy(polygon)
481 table.insert(polyNew, e1)
482 table.insert(polyNew, e2)
483 table.remove(polyCopy,1)
488 while (i<=#polyCopy and find==false) do
489 bool1 = (polyCopy[i][1] == polyNew[j])
490 bool2 = (polyCopy[i][2] == polyNew[j])
491 if(bool1 or bool2) then -- the edge has a common point with polyNew[j]
493 table.insert(polyNew, polyCopy[i][1])
495 table.remove(polyCopy,i)
497 elseif(not bool2) then
498 table.insert(polyNew, polyCopy[i][2])
500 table.remove(polyCopy,i)
511 function TeXaddOnePointTikZ(listPoints,P,step,bbox,color,colorBack, colorNew, colorCircle,colorBbox)
513 -- build the triangulation
514 triangulation = BowyerWatson(listPoints,bbox)
515 badTriangles = buildBadTriangles(P,triangulation)
516 for i=1,#listPoints do
517 output = output .. "\\coordinate (MeshPoints".. i .. ") at (" .. listPoints[i].x .. "," .. listPoints[i].y .. ");"
519 if(step == "badT") then
521 for i=1,#triangulation do
522 PointI = listPoints[triangulation[i][1]]
523 PointJ = listPoints[triangulation[i][2]]
524 PointK = listPoints[triangulation[i][3]]
525 if(triangulation[i].type == "bbox") then
526 output = output .. "\\draw[color="..colorBbox.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
528 output = output .. "\\draw[color="..color.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
531 -- draw and fill the bad triangle
532 for i=1,#badTriangles do
533 PointI = listPoints[triangulation[badTriangles[i]][1]]
534 PointJ = listPoints[triangulation[badTriangles[i]][2]]
535 PointK = listPoints[triangulation[badTriangles[i]][3]]
536 output = output .. "\\draw[fill="..colorBack.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
538 -- draw the circoncircle
539 for i=1,#badTriangles do
540 PointI = listPoints[triangulation[badTriangles[i]][1]]
541 PointJ = listPoints[triangulation[badTriangles[i]][2]]
542 PointK = listPoints[triangulation[badTriangles[i]][3]]
543 center, radius = circoncircle(PointI, PointJ, PointK)
544 output = output .. "\\draw[dashed, color="..colorCircle.."] ("..center.x .. "," .. center.y .. ") circle ("..radius ..");"
548 for i=1,#listPoints do
549 if(listPoints[i].type == "bbox") then
550 output = output .. "\\draw[color="..colorBbox.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint^*_{" .. j .. "}$};"
553 output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
556 -- mark the point to add
557 output = output .. "\\draw[color="..colorNew.."] (" .. P.x ..",".. P.y .. ") node {$\\bullet$} node[anchor=north east] {$\\NewPoint$};"
558 elseif(step == "cavity") then
559 polygon = buildCavity(badTriangles, triangulation)
560 polyNew = cleanPoly(polygon)
561 -- remove the bad triangles
562 for j=1,#badTriangles do
563 table.remove(triangulation,badTriangles[j]-(j-1))
565 -- draw the triangles
566 for i=1,#triangulation do
567 PointI = listPoints[triangulation[i][1]]
568 PointJ = listPoints[triangulation[i][2]]
569 PointK = listPoints[triangulation[i][3]]
570 if(triangulation[i].type == "bbox") then
571 output = output .. "\\draw[color="..colorBbox.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
573 output = output .. "\\draw[color="..color.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
576 -- fill and draw the cavity
579 PointI = listPoints[polyNew[i]]
580 path = path .. "(".. PointI.x ..",".. PointI.y ..")--"
582 output = output .. "\\draw[color="..colorNew..",fill ="..colorBack..", thick] " .. path .. "cycle;"
583 -- mark the points of the mesh
585 for i=1,#listPoints do
586 if(listPoints[i].type == "bbox") then
587 output = output .. "\\draw[color="..colorBbox.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint^*_{" .. j .. "}$};"
590 output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
593 -- mark the adding point
594 output = output .. "\\draw[color="..colorNew.."] (" .. P.x ..",".. P.y .. ") node {$\\bullet$} node[anchor=north east] {$\\NewPoint$};"
595 elseif(step == "newT") then
596 polygon = buildCavity(badTriangles, triangulation)
597 polyNew = cleanPoly(polygon)
598 -- remove the bad triangles
599 for j=1,#badTriangles do
600 table.remove(triangulation,badTriangles[j]-(j-1))
602 -- draw the triangle of the triangulation
603 for i=1,#triangulation do
604 PointI = listPoints[triangulation[i][1]]
605 PointJ = listPoints[triangulation[i][2]]
606 PointK = listPoints[triangulation[i][3]]
607 if(triangulation[i].type == "bbox") then
608 output = output .. "\\draw[color="..colorBbox.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
610 output = output .. "\\draw[color="..color.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
613 -- fill and draw the cavity
616 PointI = listPoints[polyNew[i]]
617 path = path .. "(".. PointI.x ..",".. PointI.y ..")--"
619 output = output .. "\\draw[color="..colorNew..",fill ="..colorBack..", thick] " .. path .. "cycle;"
620 -- draw the new triangles composed by the edges of the polygon and the added point
622 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 ..");"
623 output = output .. "\\draw[color="..colorNew..", thick]".."(".. listPoints[polygon[i][1]].x .. "," .. listPoints[polygon[i][1]].y .. ") -- (" .. P.x .. "," .. P.y ..");"
624 output = output .. "\\draw[color="..colorNew..", thick]".."(".. listPoints[polygon[i][2]].x .. "," .. listPoints[polygon[i][2]].y .. ") -- (" .. P.x .. "," .. P.y ..");"
628 for i=1,#listPoints do
629 if(listPoints[i].type == "bbox") then
630 output = output .. "\\draw[color="..colorBbox.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint^*_{" .. j .. "}$};"
633 output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
636 -- mark the added point
637 output = output .. "\\draw[color="..colorNew.."] (" .. P.x ..",".. P.y .. ") node {$\\bullet$} node[anchor=north east] {$\\NewPoint$};"
642 function TeXaddOnePointMPBW(listPoints,P,step,bbox)
644 output = output .. "pair MeshPoints[];"
645 -- build the triangulation
646 triangulation = BowyerWatson(listPoints,bbox)
647 badTriangles = buildBadTriangles(P,triangulation)
648 for i=1,#listPoints do
649 output = output .. "MeshPoints[".. i .. "] = (" .. listPoints[i].x .. "," .. listPoints[i].y .. ")*u;"
651 if(step == "badT") then
653 for i=1,#triangulation do
654 PointI = listPoints[triangulation[i][1]]
655 PointJ = listPoints[triangulation[i][2]]
656 PointK = listPoints[triangulation[i][3]]
657 if(triangulation[i].type == "bbox") then
658 output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolorBbox;"
660 output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolor;"
663 -- draw and fill the bad triangle
664 for i=1,#badTriangles do
665 PointI = listPoints[triangulation[badTriangles[i]][1]]
666 PointJ = listPoints[triangulation[badTriangles[i]][2]]
667 PointK = listPoints[triangulation[badTriangles[i]][3]]
668 output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolor;"
669 output = output .. "fill (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolorBack;"
671 -- draw the circoncircle
672 for i=1,#badTriangles do
673 PointI = listPoints[triangulation[badTriangles[i]][1]]
674 PointJ = listPoints[triangulation[badTriangles[i]][2]]
675 PointK = listPoints[triangulation[badTriangles[i]][3]]
676 center, radius = circoncircle(PointI, PointJ, PointK)
677 output = output .. "draw fullcircle scaled ("..radius .."*2u) shifted ("..center.x .. "*u," .. center.y .. "*u) dashed evenly withcolor \\luameshmpcolorCircle;"
681 for i=1,#listPoints do
682 if(listPoints[i].type == "bbox") then
683 output = output .. "dotlabel.llft (btex $\\MeshPoint^{*}_{" .. j .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolorBbox ;"
686 output = output .. "dotlabel.llft (btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolor ;"
689 -- mark the point to add
690 output = output .. "dotlabel.llft (btex $\\NewPoint$ etex,(" .. P.x ..",".. P.y .. ")*u) withcolor \\luameshmpcolorNew;"
691 elseif(step == "cavity") then
692 polygon = buildCavity(badTriangles, triangulation)
693 polyNew = cleanPoly(polygon)
694 -- remove the bad triangles
695 for j=1,#badTriangles do
696 table.remove(triangulation,badTriangles[j]-(j-1))
698 -- draw the triangles
699 for i=1,#triangulation do
700 PointI = listPoints[triangulation[i][1]]
701 PointJ = listPoints[triangulation[i][2]]
702 PointK = listPoints[triangulation[i][3]]
703 if(triangulation[i].type == "bbox") then
704 output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolorBbox;"
706 output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolor;"
709 -- fill and draw the cavity
712 PointI = listPoints[polyNew[i]]
713 path = path .. "(".. PointI.x ..",".. PointI.y ..")*u--"
715 output = output .. "fill " .. path .. "cycle withcolor \\luameshmpcolorBack;"
716 output = output .. "draw " .. path .. "cycle withcolor \\luameshmpcolorNew withpen pencircle scaled 1pt;"
717 -- mark the points of the mesh
719 for i=1,#listPoints do
720 if(listPoints[i].type == "bbox") then
721 output = output .. "dotlabel.llft (btex $\\MeshPoint^{*}_{" .. j .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolorBbox ;"
724 output = output .. "dotlabel.llft (btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolor ;"
727 -- mark the adding point
728 output = output .. "dotlabel.llft (btex $\\NewPoint$ etex,(" .. P.x ..",".. P.y .. ")*u) withcolor \\luameshmpcolorNew ;"
729 elseif(step == "newT") then
730 polygon = buildCavity(badTriangles, triangulation)
731 polyNew = cleanPoly(polygon)
732 -- remove the bad triangles
733 for j=1,#badTriangles do
734 table.remove(triangulation,badTriangles[j]-(j-1))
736 -- draw the triangle of the triangulation
737 for i=1,#triangulation do
738 PointI = listPoints[triangulation[i][1]]
739 PointJ = listPoints[triangulation[i][2]]
740 PointK = listPoints[triangulation[i][3]]
741 if(triangulation[i].type == "bbox") then
742 output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolorBbox;"
744 output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolor;"
750 PointI = listPoints[polyNew[i]]
751 path = path .. "(".. PointI.x ..",".. PointI.y ..")*u--"
753 output = output .. "fill " .. path .. "cycle withcolor \\luameshmpcolorBack;"
754 -- draw the new triangles composed by the edges of the polygon and the added point
756 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;"
757 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;"
758 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;"
762 for i=1,#listPoints do
763 if(listPoints[i].type == "bbox") then
764 output = output .. "dotlabel.llft (btex $\\MeshPoint^{*}_{" .. j .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolorBbox ;"
767 output = output .. "dotlabel.llft (btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolor ;"
770 -- mark the added point
771 output = output .. "dotlabel.llft (btex $\\NewPoint$ etex,(" .. P.x ..",".. P.y .. ")*u) withcolor \\luameshmpcolorNew ;"
777 -- build the list of points extern and stop at nbr
778 function buildListExt(chaine, stop)
780 io.input(chaine) -- open the file
782 lines=string.explode(text,"\n+") -- all the lines
783 for i=1,tonumber(stop) do
784 xy=string.explode(lines[i]," +")
785 table.insert(listPoints,{x=tonumber(xy[1]),y=tonumber(xy[2])})
787 xy=string.explode(lines[stop+1]," +")
788 point={x=tonumber(xy[1]),y=tonumber(xy[2])}
789 return point, listPoints
793 function TeXOnePointTikZBW(chaine,point,step,scale,mode,bbox,color,colorBack,colorNew,colorCircle,colorBbox)
795 Sx,Sy=string.match(point,"%((.+),(.+)%)")
797 listPoints = buildList(chaine, mode)
800 P, listPoints = buildListExt(chaine,tonumber(point))
802 output = TeXaddOnePointTikZ(listPoints,P,step,bbox,color,colorBack,colorNew,colorCircle,colorBbox)
803 output = "\\noindent\\begin{tikzpicture}[x="..scale..",y="..scale.."]".. output .. "\\end{tikzpicture}"
807 function TeXOnePointTikZBWinc(chaine,point,beginning, ending,step,scale,mode,bbox,color,colorBack,colorNew,colorCircle,colorBbox)
809 Sx,Sy=string.match(point,"%((.+),(.+)%)")
811 listPoints = buildList(chaine, mode)
814 P, listPoints = buildListExt(chaine,tonumber(point))
816 output = TeXaddOnePointTikZ(listPoints,P,step,bbox,color,colorBack,colorNew,colorCircle,colorBbox)
817 output = "\\noindent\\begin{tikzpicture}[x="..scale..",y="..scale.."]".. beginning..output ..ending.. "\\end{tikzpicture}"
821 function TeXOnePointMPBW(chaine,point,step,scale,mode,bbox)
823 Sx,Sy=string.match(point,"%((.+),(.+)%)")
825 listPoints = buildList(chaine, mode)
828 P, listPoints = buildListExt(chaine,tonumber(point))
830 output = TeXaddOnePointMPBW(listPoints,P,step,bbox)
831 output = "\\leavevmode\\begin{mplibcode}beginfig(0);u:="..scale..";".. output .. "endfig;\\end{mplibcode}"
835 function TeXOnePointMPBWinc(chaine,point,beginning,ending,step,scale,mode,bbox)
837 Sx,Sy=string.match(point,"%((.+),(.+)%)")
839 listPoints = buildList(chaine, mode)
842 P, listPoints = buildListExt(chaine,tonumber(point))
844 output = TeXaddOnePointMPBW(listPoints,P,step,bbox)
845 output = "\\begin{mplibcode}u:="..scale..";"..beginning .. output .. ending .. "\\end{mplibcode}"