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,#listPoints do
 259       output = output .. "\\coordinate (MeshPoints".. i .. ") at  (" .. listPoints[i].x .. "," .. listPoints[i].y .. ");"
 260    end
 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;"
 267       else
 268          output = output .. "\\draw[color="..color.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
 269       end
 270    end
 271    if(points=="points") then
 272       j=1
 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 .. "}$};"
 276             j=j+1
 277          else
 278             output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
 279          end
 280       end
 281    end
 282    return output
 283 end
 284
 285
 286 -- trace a triangulation with MP
 287 function traceMeshMP(listPoints, triangulation,points)
 288    output = "";
 289    output = output .. " pair MeshPoints[];"
 290    for i=1,#listPoints do
 291       output = output .. "MeshPoints[".. i .. "] = (" .. listPoints[i].x .. "," .. listPoints[i].y .. ")*u;"
 292    end
 293
 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;"
 300       else
 301          output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolor;"
 302       end
 303    end
 304    if(points=="points") then
 305       j=1
 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 ;"
 309             j=j+1
 310          else
 311             output = output .. "dotlabel.llft (btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolor ;"
 312          end
 313       end
 314    end
 315    return output
 316 end
 317
 318
 319 -- buildMesh with MP
 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}"
 325    tex.sprint(output)
 326 end
 327
 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}"
 334    tex.sprint(output)
 335 end
 336
 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}"
 343    tex.sprint(output)
 344 end
 345
 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}"
 352    tex.sprint(output)
 353 end
 354
 355
 356 -- print points of the mesh
 357 function tracePointsMP(listPoints,points)
 358    output = "";
 359    output = output .. " pair MeshPoints[];"
 360    for i=1,#listPoints do
 361       output = output .. "MeshPoints[".. i .. "] = (" .. listPoints[i].x .. "," .. listPoints[i].y .. ")*u;"
 362    end
 363    if(points=="points") then
 364       j=1
 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 ;"
 368             j=j+1
 369          else
 370             output = output .. "dotlabel.llft (btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolor ;"
 371          end
 372       end
 373    else
 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;"
 377          else
 378             output = output .. "drawdot (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u  withcolor \\luameshmpcolor withpen pencircle scaled 3;"
 379          end
 380       end
 381    end
 382    return output
 383 end
 384
 385 -- print points of the mesh
 386 function tracePointsTikZ(listPoints,points,color,colorBbox)
 387    output = "";
 388    for i=1,#listPoints do
 389       output = output .. "\\coordinate (MeshPoints".. i .. ") at  (" .. listPoints[i].x .. "," .. listPoints[i].y .. ");"
 390    end
 391    if(points=="points") then
 392       j=1
 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 .. "}$};"
 396             j = j+1
 397          else
 398             output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
 399          end
 400       end
 401    else
 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$} ;"
 405          else
 406             output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} ;"
 407          end
 408       end
 409    end
 410    return output
 411 end
 412
 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)
 418    end
 419    output = tracePointsMP(listPoints,points)
 420    output = "\\leavevmode\\begin{mplibcode}beginfig(0);u:="..scale.. ";" .. output .."endfig;\\end{mplibcode}"
 421    tex.sprint(output)
 422 end
 423
 424
 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)
 430    end
 431    output = tracePointsMP(listPoints,points)
 432    output = "\\leavevmode\\begin{mplibcode}u:="..scale..";"..beginning .. output .. ending .. "\\end{mplibcode}"
 433    tex.sprint(output)
 434 end
 435
 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)
 441    end
 442    output = tracePointsTikZ(listPoints,points,color,colorBbox)
 443    output = "\\noindent\\begin{tikzpicture}[x=" .. scale .. ",y=" .. scale .."]" .. output .."\\end{tikzpicture}"
 444    tex.sprint(output)
 445 end
 446
 447
 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)
 453    end
 454    output = tracePointsTikZ(listPoints,points,color,colorBbox)
 455    output = "\\noindent\\begin{tikzpicture}[x=" .. scale .. ",y=" .. scale .."]" ..beginning.. output..ending .."\\end{tikzpicture}"
 456    tex.sprint(output)
 457 end
 458
 459
 460 -- buildMesh
 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")
 465 end
 466
 467 local function shallowCopy(original)
 468    local copy = {}
 469    for key, value in pairs(original) do
 470       copy[key] = value
 471    end
 472    return copy
 473 end
 474
 475 -- function give a real polygon without repeting points
 476 function cleanPoly(polygon)
 477    polyNew = {}
 478    polyCopy = shallowCopy(polygon)
 479    e1 = polyCopy[1][1]
 480    e2 = polyCopy[1][2]
 481    table.insert(polyNew, e1)
 482    table.insert(polyNew, e2)
 483    table.remove(polyCopy,1)
 484    j = 2
 485    while #polyCopy>1 do
 486       i=1
 487       find = false
 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]
 492             if(not bool1) then
 493                table.insert(polyNew, polyCopy[i][1])
 494                find = true
 495                table.remove(polyCopy,i)
 496                j = j+1
 497             elseif(not bool2) then
 498                table.insert(polyNew, polyCopy[i][2])
 499                find = true
 500                table.remove(polyCopy,i)
 501                j = j+1
 502             end
 503          end
 504          i=i+1
 505       end
 506    end
 507    return polyNew
 508 end
 509
 510 --
 511 function TeXaddOnePointTikZ(listPoints,P,step,bbox,color,colorBack, colorNew, colorCircle,colorBbox)
 512    output = ""
 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 .. ");"
 518    end
 519    if(step == "badT") then
 520       -- draw all triangle
 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;"
 527          else
 528             output = output .. "\\draw[color="..color.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
 529          end
 530       end
 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;"
 537       end
 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 ..");"
 545       end
 546       -- mark the points
 547       j=1
 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 .. "}$};"
 551             j = j+1
 552          else
 553             output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
 554          end
 555       end
 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))
 564       end
 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;"
 572          else
 573             output = output .. "\\draw[color="..color.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
 574          end
 575       end
 576       -- fill and draw the cavity
 577       path = ""
 578       for i=1,#polyNew do
 579          PointI = listPoints[polyNew[i]]
 580          path = path .. "(".. PointI.x ..",".. PointI.y ..")--"
 581       end
 582       output = output .. "\\draw[color="..colorNew..",fill ="..colorBack..", thick] " .. path .. "cycle;"
 583       -- mark the points of the mesh
 584       j=1
 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 .. "}$};"
 588             j=j+1
 589          else
 590             output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
 591          end
 592       end
 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))
 601       end
 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;"
 609          else
 610             output = output .. "\\draw[color="..color.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
 611          end
 612       end
 613       -- fill and draw the cavity
 614       path = ""
 615       for i=1,#polyNew do
 616          PointI = listPoints[polyNew[i]]
 617          path = path .. "(".. PointI.x ..",".. PointI.y ..")--"
 618       end
 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
 621       for i=1,#polygon do
 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 ..");"
 625       end
 626       -- mark points
 627       j=1
 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 .. "}$};"
 631             j=j+1
 632          else
 633             output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
 634          end
 635       end
 636       -- mark the added point
 637       output = output .. "\\draw[color="..colorNew.."] (" .. P.x ..",".. P.y .. ") node {$\\bullet$} node[anchor=north east] {$\\NewPoint$};"
 638    end
 639    return output
 640 end
 641
 642 function TeXaddOnePointMPBW(listPoints,P,step,bbox)
 643    output = "";
 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;"
 650    end
 651    if(step == "badT") then
 652       -- draw all triangle
 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;"
 659          else
 660             output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolor;"
 661          end
 662       end
 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;"
 670       end
 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;"
 678       end
 679       -- mark the points
 680       j=1
 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 ;"
 684             j=j+1
 685          else
 686             output = output .. "dotlabel.llft (btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolor ;"
 687          end
 688       end
 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))
 697       end
 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;"
 705          else
 706             output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolor;"
 707          end
 708       end
 709       -- fill and draw the cavity
 710       path = ""
 711       for i=1,#polyNew do
 712          PointI = listPoints[polyNew[i]]
 713          path = path .. "(".. PointI.x ..",".. PointI.y ..")*u--"
 714       end
 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
 718       j=1
 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 ;"
 722             j=j+1
 723          else
 724             output = output .. "dotlabel.llft (btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolor ;"
 725          end
 726       end
 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))
 735       end
 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;"
 743          else
 744             output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolor;"
 745          end
 746       end
 747       -- fill  the cavity
 748       path = ""
 749       for i=1,#polyNew do
 750          PointI = listPoints[polyNew[i]]
 751          path = path .. "(".. PointI.x ..",".. PointI.y ..")*u--"
 752       end
 753       output = output .. "fill " .. path .. "cycle withcolor \\luameshmpcolorBack;"
 754       -- draw the new triangles composed by the edges of the polygon and the added point
 755       for i=1,#polygon do
 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;"
 759       end
 760       -- mark points
 761       j=1
 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 ;"
 765             j=j+1
 766          else
 767             output = output .. "dotlabel.llft (btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolor ;"
 768          end
 769       end
 770       -- mark the added point
 771       output = output .. "dotlabel.llft (btex $\\NewPoint$ etex,(" .. P.x ..",".. P.y .. ")*u) withcolor \\luameshmpcolorNew ;"
 772    end
 773    return output
 774 end
 775
 776
 777 -- build the list of points extern and stop at nbr
 778 function buildListExt(chaine, stop)
 779    listPoints = {}
 780    io.input(chaine) -- open the file
 781    text=io.read("*all")
 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])})
 786    end
 787    xy=string.explode(lines[stop+1]," +")
 788    point={x=tonumber(xy[1]),y=tonumber(xy[2])}
 789    return point, listPoints
 790 end
 791
 792
 793 function TeXOnePointTikZBW(chaine,point,step,scale,mode,bbox,color,colorBack,colorNew,colorCircle,colorBbox)
 794    if(mode=="int") then
 795       Sx,Sy=string.match(point,"%((.+),(.+)%)")
 796       P = {x=Sx, y=Sy}
 797       listPoints = buildList(chaine, mode)
 798    else
 799       -- point is a number
 800       P, listPoints = buildListExt(chaine,tonumber(point))
 801    end
 802    output = TeXaddOnePointTikZ(listPoints,P,step,bbox,color,colorBack,colorNew,colorCircle,colorBbox)
 803    output = "\\noindent\\begin{tikzpicture}[x="..scale..",y="..scale.."]".. output .. "\\end{tikzpicture}"
 804    tex.sprint(output)
 805 end
 806
 807 function TeXOnePointTikZBWinc(chaine,point,beginning, ending,step,scale,mode,bbox,color,colorBack,colorNew,colorCircle,colorBbox)
 808    if(mode=="int") then
 809       Sx,Sy=string.match(point,"%((.+),(.+)%)")
 810       P = {x=Sx, y=Sy}
 811       listPoints = buildList(chaine, mode)
 812    else
 813       -- point is a number
 814       P, listPoints = buildListExt(chaine,tonumber(point))
 815    end
 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}"
 818    tex.sprint(output)
 819 end
 820
 821 function TeXOnePointMPBW(chaine,point,step,scale,mode,bbox)
 822    if(mode=="int") then
 823       Sx,Sy=string.match(point,"%((.+),(.+)%)")
 824       P = {x=Sx, y=Sy}
 825       listPoints = buildList(chaine, mode)
 826    else
 827       -- point is a number
 828       P, listPoints = buildListExt(chaine,tonumber(point))
 829    end
 830    output = TeXaddOnePointMPBW(listPoints,P,step,bbox)
 831    output = "\\leavevmode\\begin{mplibcode}beginfig(0);u:="..scale..";".. output .. "endfig;\\end{mplibcode}"
 832    tex.sprint(output)
 833 end
 834
 835 function TeXOnePointMPBWinc(chaine,point,beginning,ending,step,scale,mode,bbox)
 836    if(mode=="int") then
 837       Sx,Sy=string.match(point,"%((.+),(.+)%)")
 838       P = {x=Sx, y=Sy}
 839       listPoints = buildList(chaine, mode)
 840    else
 841       -- point is a number
 842       P, listPoints = buildListExt(chaine,tonumber(point))
 843    end
 844    output = TeXaddOnePointMPBW(listPoints,P,step,bbox)
 845    output = "\\begin{mplibcode}u:="..scale..";"..beginning .. output .. ending .. "\\end{mplibcode}"
 846    tex.sprint(output)
 847 end