luamesh.lua

   1 require "polygon-luamesh"
   2
   3 -- Bowyer and Watson algorithm
   4 -- Delaunay meshing
   5 function BowyerWatson (listPoints,bbox)
   6    local triangulation = {}
   7    local lgth = #listPoints
   8    -- add four points to listPoints to have a bounding box
   9    listPoints = buildBoundingBox(listPoints)
  10    -- the first triangle
  11    triangulation[1] = {lgth+1, lgth+2, lgth+3,type="bbox"}
  12    -- the second triangle
  13    triangulation[2] = {lgth+1, lgth+3, lgth+4,type="bbox"}
  14    -- add points one by one
  15    for i=1,lgth do
  16       -- find the triangles which the circumcircle contained the point to add
  17       badTriangles = buildBadTriangles(listPoints[i],triangulation)
  18       -- build the polygon of the cavity containing the point to add
  19       polygon = buildCavity(badTriangles, triangulation)
  20       -- remove the bad triangles
  21       for j=1,#badTriangles do
  22          table.remove(triangulation,badTriangles[j]-(j-1))
  23       end
  24       -- build the new triangles and add them to triangulation
  25       for j=1,#polygon do
  26          if((polygon[j][1]>lgth) or (polygon[j][2]>lgth) or (i>lgth)) then
  27             table.insert(triangulation,{polygon[j][1],polygon[j][2],i,type="bbox"})
  28          else
  29             table.insert(triangulation,{polygon[j][1],polygon[j][2],i,type="in"})
  30          end
  31       end
  32    end  -- end adding points of the listPoints
  33    -- remove bounding box
  34    if(bbox ~= "bbox") then
  35       triangulation = removeBoundingBox(triangulation,lgth)
  36       table.remove(listPoints,lgth+1)
  37       table.remove(listPoints,lgth+1)
  38       table.remove(listPoints,lgth+1)
  39       table.remove(listPoints,lgth+1)
  40    end
  41    return triangulation
  42 end
  43
  44
  45 function buildBoundingBox(listPoints)
  46    -- listPoints : list of points
  47    -- epsV        : parameter for the distance of the bounding box
  48    local xmin, xmax, ymin, ymax, eps
  49    xmin = 1000
  50    ymin = 1000
  51    xmax = -1000
  52    ymax = -1000
  53    for i=1,#listPoints do
  54       if (listPoints[i].x < xmin) then
  55          xmin = listPoints[i].x
  56       end
  57       if (listPoints[i].x > xmax) then
  58          xmax = listPoints[i].x
  59       end
  60       if (listPoints[i].y < ymin) then
  61          ymin = listPoints[i].y
  62       end
  63       if (listPoints[i].y > ymax) then
  64          ymax = listPoints[i].y
  65       end
  66    end
  67    eps = math.max(math.abs(xmax-xmin),math.abs(ymax-ymin))*0.15
  68    xmin = xmin - eps
  69    xmax = xmax + eps
  70    ymin = ymin - eps
  71    ymax = ymax + eps
  72    -- add points of the bounding box in last positions
  73    table.insert(listPoints,{x=xmin,y=ymin,type="bbox"})
  74    table.insert(listPoints,{x=xmin,y=ymax,type="bbox"})
  75    table.insert(listPoints,{x=xmax,y=ymax,type="bbox"})
  76    table.insert(listPoints,{x=xmax,y=ymin,type="bbox"})
  77    return listPoints
  78 end
  79
  80 function removeBoundingBox(triangulation,lgth)
  81    -- build the four bounding box edge
  82    point1 = lgth+1
  83    point2 = lgth+2
  84    point3 = lgth+3
  85    point4 = lgth+4
  86    -- for all triangle
  87    newTriangulation = {}
  88    for i=1,#triangulation do
  89       boolE1 = pointInTriangle(point1,triangulation[i])
  90       boolE2 = pointInTriangle(point2,triangulation[i])
  91       boolE3 = pointInTriangle(point3,triangulation[i])
  92       boolE4 = pointInTriangle(point4,triangulation[i])
  93       if((not boolE1) and (not boolE2) and (not boolE3) and (not boolE4)) then
  94          table.insert(newTriangulation,triangulation[i])
  95       end
  96    end
  97    return newTriangulation
  98 end
  99
 100
 101 function buildBadTriangles(point, triangulation)
 102    badTriangles = {}
 103    for j=1,#triangulation do -- for all triangles
 104       A = listPoints[triangulation[j][1]]
 105       B = listPoints[triangulation[j][2]]
 106       C = listPoints[triangulation[j][3]]
 107       center, radius = circoncircle(A,B,C)
 108       CP = Vector(center,point)
 109       if(VectorNorm(CP)<radius) then -- the point belongs to the circoncirle
 110          table.insert(badTriangles,j)
 111       end
 112    end
 113    return badTriangles
 114 end
 115
 116 -- construction of the cavity composed by the bad triangles around the point to add
 117 function buildCavity(badTriangles, triangulation)
 118    polygon = {}
 119    for j=1,#badTriangles do -- for all bad triangles
 120       ind = badTriangles[j]
 121       for k=1,3 do -- for all edges
 122          edge = {triangulation[ind][k],triangulation[ind][k%3+1]}
 123          edgeBord = false
 124          for l = 1,#badTriangles do -- for all badtriangles
 125             badInd = badTriangles[l]
 126             if(badInd ~= ind) then -- if not the current one
 127                edgeBord = edgeBord or edgeInTriangle(edge,triangulation[badInd])
 128             end
 129          end --
 130          -- if the edge does not belong to another bad triangle
 131          if(edgeBord == false) then
 132             -- insert the edge to the cavity
 133             table.insert(polygon,edge)
 134          end
 135       end --
 136    end --
 137    return polygon
 138 end
 139
 140 function edgeInTriangle(e,t)
 141    in1 = false
 142    in2 = false
 143    for i=1,3 do
 144       if e[1] == t[i] then
 145          in1 = true
 146       end
 147       if e[2] == t[i] then
 148          in2 = true
 149       end
 150    end
 151    out = (in1 and in2)
 152    return out
 153 end
 154
 155 function pointInTriangle(e,t)
 156    in1 = false
 157    for i=1,3 do
 158       if e == t[i] then
 159          in1 = true
 160       end
 161    end
 162    return in1
 163 end
 164
 165
 166 function Vector(A,B)
 167    local out = {x = B.x - A.x, y = B.y - A.y}
 168    return out
 169 end
 170
 171 function VectorNorm(NP)
 172    return math.sqrt(NP.x*NP.x +NP.y*NP.y)
 173 end
 174
 175 -- circoncircle
 176 function circoncircle(M, N, P)
 177    -- Compute center and radius of the circoncircle of the triangle M N P
 178
 179    -- return : (center [Point],radius [float])
 180
 181    local MN = Vector(M,N)
 182    local NP = Vector(N,P)
 183    local PM = Vector(P,M)
 184    m = VectorNorm(NP)  -- |NP|
 185    n = VectorNorm(PM)  -- |PM|
 186    p = VectorNorm(MN)  -- |MN|
 187
 188    d = (m + n + p) * (-m + n + p) * (m - n + p) * (m + n - p)
 189    if d > 0 then
 190       rad = m * n * p / math.sqrt(d)
 191    else
 192       rad = 0
 193    end
 194    d = -2 * (M.x * NP.y + N.x * PM.y + P.x * MN.y)
 195    O = {x=0, y=0}
 196    OM = Vector(O, M)
 197    ON = Vector(O, N)
 198    OP = Vector(O, P)
 199    om2 = math.pow(VectorNorm(OM),2)  -- |OM|**2
 200    on2 = math.pow(VectorNorm(ON),2)  -- |ON|**2
 201    op2 = math.pow(VectorNorm(OP),2)  -- |OP|**2
 202    x0 = -(om2 * NP.y + on2 * PM.y + op2 * MN.y) / d
 203    y0 = (om2 * NP.x + on2 * PM.x + op2 * MN.x) / d
 204    if d == 0 then
 205       Out = {nil, nil}
 206    else
 207       Out = {x=x0, y=y0}
 208    end
 209    return Out, rad  -- (center [Point], R [float])
 210 end
 211
 212 -- compute the list of the circumcircle of a triangulation
 213 function listCircumCenter(listPoints,triangulation)
 214    list = {}
 215    for j=1,#triangulation do
 216       A = listPoints[triangulation[j][1]]
 217       B = listPoints[triangulation[j][2]]
 218       C = listPoints[triangulation[j][3]]
 219       center, radius = circoncircle(A,B,C)
 220       table.insert(list,{x=center.x,y=center.y,r=radius})
 221    end
 222    return list
 223 end
 224
 225 -- find the three neighbour triangles of T
 226 function findNeighbour(T,i,triangulation)
 227    -- T : triangle
 228    -- i : index of T in triangualation
 229    -- triangulation
 230
 231    list = {}
 232    -- define the three edge
 233    e1 = {T[1],T[2]}
 234    e2 = {T[2],T[3]}
 235    e3 = {T[3],T[1]}
 236    for j=1,#triangulation do
 237       if j~= i then
 238          if(edgeInTriangle(e1,triangulation[j])) then
 239             table.insert(list,j)
 240          end
 241          if(edgeInTriangle(e2,triangulation[j])) then
 242             table.insert(list,j)
 243          end
 244          if(edgeInTriangle(e3,triangulation[j])) then
 245             table.insert(list,j)
 246          end
 247       end
 248    end
 249    return list
 250 end
 251
 252 -- test if edge are the same (reverse)
 253 function equalEdge(e1,e2)
 254    if(((e1[1] == e2[1]) and (e1[2] == e2[2])) or ((e1[1] == e2[2]) and (e1[2] == e2[1]))) then
 255       return true
 256    else
 257       return false
 258    end
 259 end
 260
 261 -- test if the edge belongs to the list
 262 function edgeInList(e,listE)
 263    output = false
 264    for i=1,#listE do
 265       if(equalEdge(e,listE[i])) then
 266          output = true
 267       end
 268    end
 269    return output
 270 end
 271
 272 -- build the edges of the Voronoi diagram with a given triangulation
 273 function buildVoronoi(listPoints, triangulation)
 274    listCircumCircle = listCircumCenter(listPoints, triangulation)
 275    listVoronoi = {}
 276    for i=1,#listCircumCircle do
 277       listN = findNeighbour(triangulation[i],i,triangulation)
 278       for j=1,#listN do
 279          edge = {i,listN[j]}
 280          if( not edgeInList(edge, listVoronoi)) then
 281             table.insert(listVoronoi, edge)
 282          end
 283       end
 284    end
 285    return listVoronoi
 286 end
 287
 288 -------------------------- TeX
 289 -- build the list of points
 290 function buildList(chaine, mode)
 291    -- if mode = int : the list is given in the chaine string (x1,y1);(x2,y2);...;(xn,yn)
 292    -- if mode = ext : the list is given in a file line by line with space separation
 293    listPoints = {}
 294    if mode == "int" then
 295       local points = string.explode(chaine, ";")
 296       local lgth=#points
 297       for i=1,lgth do
 298          Sx,Sy=string.match(points[i],"%((.+),(.+)%)")
 299          listPoints[i]={x=tonumber(Sx),y=tonumber(Sy)}
 300       end
 301    elseif mode == "ext" then
 302       io.input(chaine) -- open the file
 303       text=io.read("*all")
 304       lines=string.explode(text,"\n+") -- all the lines
 305       tablePoints={}
 306       for i=1,#lines do
 307          xy=string.explode(lines[i]," +")
 308          listPoints[i]={x=tonumber(xy[1]),y=tonumber(xy[2])}
 309       end
 310    else
 311       print("Non existing mode")
 312    end
 313    return listPoints
 314 end
 315
 316
 317 --
 318 function rectangleList(a,b,nbrA,nbrB)
 319    stepA = a/nbrA
 320    stepB = b/nbrB
 321    listPoints = {}
 322    k=1
 323    for i=1,(nbrA+1) do
 324       for j=1,(nbrB+1) do
 325          listPoints[k] = {x = (i-1)*stepA, y=(j-1)*stepB}
 326          k=k+1
 327       end
 328    end
 329    return listPoints
 330 end
 331
 332
 333 -- trace Voronoi with MP
 334 function traceVoronoiMP(listPoints, triangulation,listVoronoi, points, tri,styleD,styleV)
 335    if(styleD == "dashed") then
 336       sDelaunay = "dashed evenly"
 337    else
 338       sDelaunay = ""
 339    end
 340    if(styleV == "dashed") then
 341       sVoronoi = "dashed evenly"
 342    else
 343       sVoronoi = ""
 344    end
 345    listCircumC = listCircumCenter(listPoints,triangulation)
 346    output = "";
 347    output = output .. " pair MeshPoints[];"
 348    for i=1,#listPoints do
 349       output = output .. "MeshPoints[".. i .. "] = (" .. listPoints[i].x .. "," .. listPoints[i].y .. ")*u;"
 350    end
 351    output = output .. " pair CircumCenters[];"
 352    for i=1,#listCircumC do
 353       output = output .. "CircumCenters[".. i .. "] = (" .. listCircumC[i].x .. "," .. listCircumC[i].y .. ")*u;"
 354    end
 355    if(tri=="show") then
 356       for i=1,#triangulation do
 357          PointI = listPoints[triangulation[i][1]]
 358          PointJ = listPoints[triangulation[i][2]]
 359          PointK = listPoints[triangulation[i][3]]
 360          if(triangulation[i].type == "bbox") then
 361             output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle "..sDelaunay.." withcolor \\luameshmpcolorBbox;"
 362          else
 363             output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle "..sDelaunay.." withcolor \\luameshmpcolor;"
 364          end
 365       end
 366    end
 367    for i=1,#listVoronoi do
 368       PointI = listCircumC[listVoronoi[i][1]]
 369       PointJ = listCircumC[listVoronoi[i][2]]
 370       output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u "..sVoronoi.." withcolor \\luameshmpcolorVoronoi;"
 371    end
 372    if(points=="points") then
 373       j=1
 374       for i=1,#listPoints do
 375          if(listPoints[i].type == "bbox") then
 376             output = output .. "dotlabel.llft (btex $\\MeshPoint^{*}_{"..j.."}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolorBbox ;"
 377             j=j+1
 378          else
 379             output = output .. "dotlabel.llft (btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolor ;"
 380          end
 381       end
 382       for i=1,#listCircumC do
 383          output = output .. "dotlabel.llft (btex $\\CircumPoint_{" .. i .. "}$ etex, (" .. listCircumC[i].x ..",".. listCircumC[i].y .. ")*u ) withcolor \\luameshmpcolorVoronoi ;"
 384       end
 385    else
 386       j=1
 387       for i=1,#listPoints do
 388          if(listPoints[i].type == "bbox") then
 389             output = output .. "drawdot  (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u  withcolor \\luameshmpcolorBbox withpen pencircle scaled 3;"
 390             j=j+1
 391          else
 392             output = output .. "drawdot  (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u  withcolor \\luameshmpcolor withpen pencircle scaled 3;"
 393          end
 394       end
 395       for i=1,#listCircumC do
 396          output = output .. "drawdot  (" .. listCircumC[i].x ..",".. listCircumC[i].y .. ")*u  withcolor \\luameshmpcolorVoronoi withpen pencircle scaled 3;"
 397       end
 398    end
 399
 400    return output
 401 end
 402
 403
 404 -- trace Voronoi with TikZ
 405 function traceVoronoiTikZ(listPoints, triangulation,listVoronoi, points, tri,color,colorBbox,colorVoronoi,styleD,styleV)
 406    if(styleD == "dashed") then
 407       sDelaunay = ",dashed"
 408    else
 409       sDelaunay = ""
 410    end
 411    if(styleV == "dashed") then
 412       sVoronoi = ",dashed"
 413    else
 414       sVoronoi = ""
 415    end
 416    listCircumC = listCircumCenter(listPoints,triangulation)
 417     output = ""
 418    for i=1,#listPoints do
 419       output = output .. "\\coordinate (MeshPoints".. i .. ") at  (" .. listPoints[i].x .. "," .. listPoints[i].y .. ");"
 420    end
 421    for i=1,#listCircumC do
 422       output = output .. "\\coordinate (CircumPoints".. i .. ") at  (" .. listCircumC[i].x .. "," .. listCircumC[i].y .. ");"
 423    end
 424    if(tri=="show") then
 425       for i=1,#triangulation do
 426          PointI = listPoints[triangulation[i][1]]
 427          PointJ = listPoints[triangulation[i][2]]
 428          PointK = listPoints[triangulation[i][3]]
 429          if(triangulation[i].type == "bbox") then
 430             output = output .. "\\draw[color="..colorBbox..sDelaunay.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
 431          else
 432             output = output .. "\\draw[color="..color..sDelaunay.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
 433          end
 434       end
 435    end
 436    for i=1,#listVoronoi do
 437       PointI = listCircumC[listVoronoi[i][1]]
 438       PointJ = listCircumC[listVoronoi[i][2]]
 439       output = output .. "\\draw[color="..colorVoronoi..sVoronoi.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..");"
 440    end
 441    if(points=="points") then
 442       j=1
 443       for i=1,#listPoints do
 444          if(listPoints[i].type == "bbox") then
 445             output = output .. "\\draw[color="..colorBbox.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint^*_{" .. j .. "}$};"
 446             j=j+1
 447          else
 448             output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
 449          end
 450       end
 451       for i=1,#listCircumC do
 452          output = output .. "\\draw[color="..colorVoronoi.."] (" .. listCircumC[i].x ..",".. listCircumC[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\CircumPoint_{" .. i .. "}$};"
 453       end
 454    else
 455       j=1
 456       for i=1,#listPoints do
 457          if(listPoints[i].type == "bbox") then
 458             output = output .. "\\draw[color="..colorBbox.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} ;"
 459             j=j+1
 460          else
 461             output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} ;"
 462          end
 463       end
 464       for i=1,#listCircumC do
 465          output = output .. "\\draw[color="..colorVoronoi.."] (" .. listCircumC[i].x ..",".. listCircumC[i].y .. ") node {$\\bullet$};"
 466       end
 467    end
 468    return output
 469 end
 470
 471
 472
 473 -- buildVoronoi with MP
 474 function buildVoronoiMPBW(chaine,mode,points,bbox,scale,tri,styleD,styleV)
 475    listPoints = buildList(chaine, mode)
 476    triangulation = BowyerWatson(listPoints,bbox)
 477    listVoronoi = buildVoronoi(listPoints, triangulation)
 478    output = traceVoronoiMP(listPoints,triangulation,listVoronoi,points,tri,styleD,styleV)
 479    output = "\\leavevmode\\begin{mplibcode}beginfig(0);u:="..scale.. ";" .. output .."endfig;\\end{mplibcode}"
 480    tex.sprint(output)
 481 end
 482
 483
 484 -- buildVoronoi with TikZ
 485 function buildVoronoiTikZBW(chaine,mode,points,bbox,scale,tri,color,colorBbox,colorVoronoi,styleD,styleV)
 486    listPoints = buildList(chaine, mode)
 487    triangulation = BowyerWatson(listPoints,bbox)
 488    listVoronoi = buildVoronoi(listPoints, triangulation)
 489    output = traceVoronoiTikZ(listPoints,triangulation,listVoronoi,points,tri,color,colorBbox,colorVoronoi,styleD,styleV)
 490    output = "\\noindent\\begin{tikzpicture}[x=" .. scale .. ",y=" .. scale .."]" .. output .."\\end{tikzpicture}"   tex.sprint(output)
 491 end
 492
 493
 494 -- buildVoronoi with MP
 495 function buildVoronoiMPBWinc(chaine,beginning, ending,mode,points,bbox,scale,tri,styleD,styleV)
 496    listPoints = buildList(chaine, mode)
 497    triangulation = BowyerWatson(listPoints,bbox)
 498    listVoronoi = buildVoronoi(listPoints, triangulation)
 499    output = traceVoronoiMP(listPoints,triangulation,listVoronoi,points,tri,styleD,styleV)
 500    output = "\\leavevmode\\begin{mplibcode}u:="..scale..";"..beginning .. output .. ending .. "\\end{mplibcode}"
 501    tex.sprint(output)
 502 end
 503
 504
 505 -- buildVoronoi with TikZ
 506 function buildVoronoiTikZBWinc(chaine,beginning, ending,mode,points,bbox,scale,tri,color,colorBbox,colorVoronoi)
 507    listPoints = buildList(chaine, mode,styleD,styleV)
 508    triangulation = BowyerWatson(listPoints,bbox)
 509    listVoronoi = buildVoronoi(listPoints, triangulation)
 510    output = traceVoronoiTikZ(listPoints,triangulation,listVoronoi,points,tri,color,colorBbox,colorVoronoi,styleD,styleV)
 511    output = "\\noindent\\begin{tikzpicture}[x=" .. scale .. ",y=" .. scale .."]" ..beginning.. output..ending .."\\end{tikzpicture}"
 512    tex.sprint(output)
 513 end
 514
 515
 516
 517 -- trace a triangulation with TikZ
 518 function traceMeshTikZ(listPoints, triangulation,points,color,colorBbox)
 519    output = ""
 520    for i=1,#listPoints do
 521       output = output .. "\\coordinate (MeshPoints".. i .. ") at  (" .. listPoints[i].x .. "," .. listPoints[i].y .. ");"
 522    end
 523    for i=1,#triangulation do
 524       PointI = listPoints[triangulation[i][1]]
 525       PointJ = listPoints[triangulation[i][2]]
 526       PointK = listPoints[triangulation[i][3]]
 527       if(triangulation[i].type == "bbox") then
 528          output = output .. "\\draw[color="..colorBbox.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
 529       else
 530          output = output .. "\\draw[color="..color.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
 531       end
 532    end
 533    if(points=="points") then
 534       j=1
 535       for i=1,#listPoints do
 536          if(listPoints[i].type == "bbox") then
 537             output = output .. "\\draw[color="..colorBbox.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint^*_{" .. j .. "}$};"
 538             j=j+1
 539          else
 540             output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
 541          end
 542       end
 543    end
 544    return output
 545 end
 546
 547
 548 -- trace a triangulation with MP
 549 function traceMeshMP(listPoints, triangulation,points)
 550    output = "";
 551    output = output .. " pair MeshPoints[];"
 552    for i=1,#listPoints do
 553       output = output .. "MeshPoints[".. i .. "] = (" .. listPoints[i].x .. "," .. listPoints[i].y .. ")*u;"
 554    end
 555    for i=1,#triangulation do
 556       PointI = listPoints[triangulation[i][1]]
 557       PointJ = listPoints[triangulation[i][2]]
 558       PointK = listPoints[triangulation[i][3]]
 559       if(triangulation[i].type == "bbox") then
 560          output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolorBbox;"
 561       else
 562          output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolor;"
 563       end
 564    end
 565    if(points=="points") then
 566       j=1
 567       for i=1,#listPoints do
 568          if(listPoints[i].type == "bbox") then
 569             output = output .. "dotlabel.llft (btex $\\MeshPoint^{*}_{"..j.."}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolorBbox ;"
 570             j=j+1
 571          else
 572             output = output .. "dotlabel.llft (btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolor ;"
 573          end
 574       end
 575    end
 576    return output
 577 end
 578
 579
 580 -- buildMesh with MP
 581 function buildMeshMPBW(chaine,mode,points,bbox,scale)
 582    listPoints = buildList(chaine, mode)
 583    triangulation = BowyerWatson(listPoints,bbox)
 584    output = traceMeshMP(listPoints, triangulation,points)
 585    output = "\\leavevmode\\begin{mplibcode}beginfig(0);u:="..scale.. ";" .. output .."endfig;\\end{mplibcode}"
 586    tex.sprint(output)
 587 end
 588
 589 -- buildMesh with MP include code
 590 function buildMeshMPBWinc(chaine,beginning, ending,mode,points,bbox,scale)
 591    listPoints = buildList(chaine, mode)
 592    triangulation = BowyerWatson(listPoints,bbox)
 593    output = traceMeshMP(listPoints, triangulation,points)
 594    output = "\\leavevmode\\begin{mplibcode}u:="..scale..";"..beginning .. output .. ending .. "\\end{mplibcode}"
 595    tex.sprint(output)
 596 end
 597
 598 -- buildMesh with TikZ
 599 function buildMeshTikZBW(chaine,mode,points,bbox,scale,color,colorBbox)
 600    listPoints = buildList(chaine, mode)
 601    triangulation = BowyerWatson(listPoints,bbox)
 602    output = traceMeshTikZ(listPoints, triangulation,points,color,colorBbox)
 603    output = "\\noindent\\begin{tikzpicture}[x=" .. scale .. ",y=" .. scale .."]" .. output .."\\end{tikzpicture}"
 604    tex.sprint(output)
 605 end
 606
 607 -- buildMesh with TikZ
 608 function buildMeshTikZBWinc(chaine,beginning, ending,mode,points,bbox,scale,color,colorBbox)
 609    listPoints = buildList(chaine, mode)
 610    triangulation = BowyerWatson(listPoints,bbox)
 611    output = traceMeshTikZ(listPoints, triangulation,points,color,colorBbox)
 612    output = "\\noindent\\begin{tikzpicture}[x=" .. scale .. ",y=" .. scale .."]" ..beginning.. output..ending .."\\end{tikzpicture}"
 613    tex.sprint(output)
 614 end
 615
 616
 617 -- print points of the mesh
 618 function tracePointsMP(listPoints,points)
 619    output = "";
 620    output = output .. " pair MeshPoints[];"
 621    for i=1,#listPoints do
 622       output = output .. "MeshPoints[".. i .. "] = (" .. listPoints[i].x .. "," .. listPoints[i].y .. ")*u;"
 623    end
 624    if(points=="points") then
 625       j=1
 626       for i=1,#listPoints do
 627          if(listPoints[i].type == "bbox") then
 628             output = output .. "dotlabel.llft (btex $\\MeshPoint^{*}_{" .. j .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolorBbox ;"
 629             j=j+1
 630          else
 631             output = output .. "dotlabel.llft (btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolor ;"
 632          end
 633       end
 634    else
 635       for i=1,#listPoints do
 636          if(listPoints[i].type == "bbox") then
 637             output = output .. "drawdot  (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u  withcolor \\luameshmpcolorBbox withpen pencircle scaled 3;"
 638          else
 639             output = output .. "drawdot (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u  withcolor \\luameshmpcolor withpen pencircle scaled 3;"
 640          end
 641       end
 642    end
 643    return output
 644 end
 645
 646 -- print points of the mesh
 647 function tracePointsTikZ(listPoints,points,color,colorBbox)
 648    output = "";
 649    for i=1,#listPoints do
 650       output = output .. "\\coordinate (MeshPoints".. i .. ") at  (" .. listPoints[i].x .. "," .. listPoints[i].y .. ");"
 651    end
 652    if(points=="points") then
 653       j=1
 654       for i=1,#listPoints do
 655          if(listPoints[i].type == "bbox") then
 656             output = output .. "\\draw[color="..colorBbox.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint^*_{" .. j .. "}$};"
 657             j = j+1
 658          else
 659             output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
 660          end
 661       end
 662    else
 663       for i=1,#listPoints do
 664          if(listPoints[i].type == "bbox") then
 665             output = output .. "\\draw[color="..colorBbox.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} ;"
 666          else
 667             output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} ;"
 668          end
 669       end
 670    end
 671    return output
 672 end
 673
 674 -- print points to mesh
 675 function printPointsMP(chaine,mode,points,bbox,scale)
 676    listPoints = buildList(chaine, mode)
 677    if(bbox == "bbox" ) then
 678       listPoints = buildBoundingBox(listPoints)
 679    end
 680    output = tracePointsMP(listPoints,points)
 681    output = "\\leavevmode\\begin{mplibcode}beginfig(0);u:="..scale.. ";" .. output .."endfig;\\end{mplibcode}"
 682    tex.sprint(output)
 683 end
 684
 685
 686 -- print points to mesh
 687 function printPointsMPinc(chaine,beginning, ending, mode,points,bbox,scale)
 688    listPoints = buildList(chaine, mode)
 689    if(bbox == "bbox" ) then
 690       listPoints = buildBoundingBox(listPoints)
 691    end
 692    output = tracePointsMP(listPoints,points)
 693    output = "\\leavevmode\\begin{mplibcode}u:="..scale..";"..beginning .. output .. ending .. "\\end{mplibcode}"
 694    tex.sprint(output)
 695 end
 696
 697 -- print points to mesh
 698 function printPointsTikZ(chaine,mode,points,bbox,scale,color,colorBbox)
 699    listPoints = buildList(chaine, mode)
 700    if(bbox == "bbox" ) then
 701       listPoints = buildBoundingBox(listPoints)
 702    end
 703    output = tracePointsTikZ(listPoints,points,color,colorBbox)
 704    output = "\\noindent\\begin{tikzpicture}[x=" .. scale .. ",y=" .. scale .."]" .. output .."\\end{tikzpicture}"
 705    tex.sprint(output)
 706 end
 707
 708
 709 -- print points to mesh
 710 function printPointsTikZinc(chaine,beginning, ending, mode,points,bbox,scale,color,colorBbox)
 711    listPoints = buildList(chaine, mode)
 712    if(bbox == "bbox" ) then
 713       listPoints = buildBoundingBox(listPoints)
 714    end
 715    output = tracePointsTikZ(listPoints,points,color,colorBbox)
 716    output = "\\noindent\\begin{tikzpicture}[x=" .. scale .. ",y=" .. scale .."]" ..beginning.. output..ending .."\\end{tikzpicture}"
 717    tex.sprint(output)
 718 end
 719
 720
 721 -- buildMesh
 722 function buildRect(largeur,a,b,nbrA, nbrB)
 723    listPoints = rectangleList(a,b,nbrA,nbrB)
 724    triangulation = BowyerWatson(listPoints,"none")
 725    traceTikZ(listPoints, triangulation,largeur,"none")
 726 end
 727
 728 local function shallowCopy(original)
 729    local copy = {}
 730    for key, value in pairs(original) do
 731       copy[key] = value
 732    end
 733    return copy
 734 end
 735
 736 -- function give a real polygon without repeting points
 737 function cleanPoly(polygon)
 738    polyNew = {}
 739    polyCopy = shallowCopy(polygon)
 740    e1 = polyCopy[1][1]
 741    e2 = polyCopy[1][2]
 742    table.insert(polyNew, e1)
 743    table.insert(polyNew, e2)
 744    table.remove(polyCopy,1)
 745    j = 2
 746    while #polyCopy>1 do
 747       i=1
 748       find = false
 749       while (i<=#polyCopy and find==false) do
 750          bool1 = (polyCopy[i][1] == polyNew[j])
 751          bool2 = (polyCopy[i][2] == polyNew[j])
 752          if(bool1 or bool2) then -- the edge has a common point with polyNew[j]
 753             if(not bool1) then
 754                table.insert(polyNew, polyCopy[i][1])
 755                find = true
 756                table.remove(polyCopy,i)
 757                j = j+1
 758             elseif(not bool2) then
 759                table.insert(polyNew, polyCopy[i][2])
 760                find = true
 761                table.remove(polyCopy,i)
 762                j = j+1
 763             end
 764          end
 765          i=i+1
 766       end
 767    end
 768    return polyNew
 769 end
 770
 771 --
 772 function TeXaddOnePointTikZ(listPoints,P,step,bbox,color,colorBack, colorNew, colorCircle,colorBbox)
 773    output = ""
 774    -- build the triangulation
 775    triangulation = BowyerWatson(listPoints,bbox)
 776    badTriangles = buildBadTriangles(P,triangulation)
 777    for i=1,#listPoints do
 778       output = output .. "\\coordinate (MeshPoints".. i .. ") at  (" .. listPoints[i].x .. "," .. listPoints[i].y .. ");"
 779    end
 780    if(step == "badT") then
 781       -- draw all triangle
 782       for i=1,#triangulation do
 783          PointI = listPoints[triangulation[i][1]]
 784          PointJ = listPoints[triangulation[i][2]]
 785          PointK = listPoints[triangulation[i][3]]
 786          if(triangulation[i].type == "bbox") then
 787             output = output .. "\\draw[color="..colorBbox.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
 788          else
 789             output = output .. "\\draw[color="..color.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
 790          end
 791       end
 792       -- draw and fill the bad triangle
 793       for i=1,#badTriangles do
 794          PointI = listPoints[triangulation[badTriangles[i]][1]]
 795          PointJ = listPoints[triangulation[badTriangles[i]][2]]
 796          PointK = listPoints[triangulation[badTriangles[i]][3]]
 797          output = output .. "\\draw[fill="..colorBack.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
 798       end
 799       -- draw the circoncircle
 800       for i=1,#badTriangles do
 801          PointI = listPoints[triangulation[badTriangles[i]][1]]
 802          PointJ = listPoints[triangulation[badTriangles[i]][2]]
 803          PointK = listPoints[triangulation[badTriangles[i]][3]]
 804          center, radius = circoncircle(PointI, PointJ, PointK)
 805          output = output .. "\\draw[dashed, color="..colorCircle.."] ("..center.x .. "," .. center.y .. ") circle ("..radius ..");"
 806       end
 807       -- mark the points
 808       j=1
 809       for i=1,#listPoints do
 810          if(listPoints[i].type == "bbox") then
 811             output = output .. "\\draw[color="..colorBbox.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint^*_{" .. j .. "}$};"
 812             j = j+1
 813          else
 814             output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
 815          end
 816       end
 817       -- mark the point to add
 818       output = output .. "\\draw[color="..colorNew.."] (" .. P.x ..",".. P.y .. ") node {$\\bullet$} node[anchor=north east] {$\\NewPoint$};"
 819    elseif(step == "cavity") then
 820       polygon = buildCavity(badTriangles, triangulation)
 821       polyNew = cleanPoly(polygon)
 822       -- remove the bad triangles
 823       for j=1,#badTriangles do
 824          table.remove(triangulation,badTriangles[j]-(j-1))
 825       end
 826       -- draw the triangles
 827       for i=1,#triangulation do
 828          PointI = listPoints[triangulation[i][1]]
 829          PointJ = listPoints[triangulation[i][2]]
 830          PointK = listPoints[triangulation[i][3]]
 831          if(triangulation[i].type == "bbox") then
 832             output = output .. "\\draw[color="..colorBbox.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
 833          else
 834             output = output .. "\\draw[color="..color.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
 835          end
 836       end
 837       -- fill and draw the cavity
 838       path = ""
 839       for i=1,#polyNew do
 840          PointI = listPoints[polyNew[i]]
 841          path = path .. "(".. PointI.x ..",".. PointI.y ..")--"
 842       end
 843       output = output .. "\\draw[color="..colorNew..",fill ="..colorBack..", thick] " .. path .. "cycle;"
 844       -- mark the points of the mesh
 845       j=1
 846       for i=1,#listPoints do
 847          if(listPoints[i].type == "bbox") then
 848             output = output .. "\\draw[color="..colorBbox.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint^*_{" .. j .. "}$};"
 849             j=j+1
 850          else
 851             output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
 852          end
 853       end
 854       -- mark the adding point
 855       output = output .. "\\draw[color="..colorNew.."] (" .. P.x ..",".. P.y .. ") node {$\\bullet$} node[anchor=north east] {$\\NewPoint$};"
 856    elseif(step == "newT") then
 857       polygon = buildCavity(badTriangles, triangulation)
 858       polyNew = cleanPoly(polygon)
 859       -- remove the bad triangles
 860       for j=1,#badTriangles do
 861          table.remove(triangulation,badTriangles[j]-(j-1))
 862       end
 863       -- draw the triangle of the triangulation
 864       for i=1,#triangulation do
 865          PointI = listPoints[triangulation[i][1]]
 866          PointJ = listPoints[triangulation[i][2]]
 867          PointK = listPoints[triangulation[i][3]]
 868          if(triangulation[i].type == "bbox") then
 869             output = output .. "\\draw[color="..colorBbox.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
 870          else
 871             output = output .. "\\draw[color="..color.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
 872          end
 873       end
 874       -- fill and draw the cavity
 875       path = ""
 876       for i=1,#polyNew do
 877          PointI = listPoints[polyNew[i]]
 878          path = path .. "(".. PointI.x ..",".. PointI.y ..")--"
 879       end
 880       output = output .. "\\draw[color="..colorNew..",fill ="..colorBack..", thick] " .. path .. "cycle;"
 881       -- draw the new triangles composed by the edges of the polygon and the added point
 882       for i=1,#polygon do
 883          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 ..");"
 884          output = output .. "\\draw[color="..colorNew..", thick]".."(".. listPoints[polygon[i][1]].x .. "," .. listPoints[polygon[i][1]].y .. ") -- (" .. P.x .. "," .. P.y ..");"
 885          output = output .. "\\draw[color="..colorNew..", thick]".."(".. listPoints[polygon[i][2]].x .. "," .. listPoints[polygon[i][2]].y .. ") -- (" .. P.x .. "," .. P.y ..");"
 886       end
 887       -- mark points
 888       j=1
 889       for i=1,#listPoints do
 890          if(listPoints[i].type == "bbox") then
 891             output = output .. "\\draw[color="..colorBbox.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint^*_{" .. j .. "}$};"
 892             j=j+1
 893          else
 894             output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
 895          end
 896       end
 897       -- mark the added point
 898       output = output .. "\\draw[color="..colorNew.."] (" .. P.x ..",".. P.y .. ") node {$\\bullet$} node[anchor=north east] {$\\NewPoint$};"
 899    end
 900    return output
 901 end
 902
 903 function TeXaddOnePointMPBW(listPoints,P,step,bbox)
 904    output = "";
 905    output = output .. "pair MeshPoints[];"
 906    -- build the triangulation
 907    triangulation = BowyerWatson(listPoints,bbox)
 908    badTriangles = buildBadTriangles(P,triangulation)
 909    for i=1,#listPoints do
 910       output = output .. "MeshPoints[".. i .. "] = (" .. listPoints[i].x .. "," .. listPoints[i].y .. ")*u;"
 911    end
 912    if(step == "badT") then
 913       -- draw all triangle
 914       for i=1,#triangulation do
 915          PointI = listPoints[triangulation[i][1]]
 916          PointJ = listPoints[triangulation[i][2]]
 917          PointK = listPoints[triangulation[i][3]]
 918          if(triangulation[i].type == "bbox") then
 919             output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolorBbox;"
 920          else
 921             output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolor;"
 922          end
 923       end
 924       -- draw and fill the bad triangle
 925       for i=1,#badTriangles do
 926          PointI = listPoints[triangulation[badTriangles[i]][1]]
 927          PointJ = listPoints[triangulation[badTriangles[i]][2]]
 928          PointK = listPoints[triangulation[badTriangles[i]][3]]
 929          output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolor;"
 930          output = output .. "fill (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolorBack;"
 931       end
 932       -- draw the circoncircle
 933       for i=1,#badTriangles do
 934          PointI = listPoints[triangulation[badTriangles[i]][1]]
 935          PointJ = listPoints[triangulation[badTriangles[i]][2]]
 936          PointK = listPoints[triangulation[badTriangles[i]][3]]
 937          center, radius = circoncircle(PointI, PointJ, PointK)
 938          output = output .. "draw fullcircle scaled ("..radius .."*2u) shifted ("..center.x .. "*u," .. center.y .. "*u) dashed evenly withcolor \\luameshmpcolorCircle;"
 939       end
 940       -- mark the points
 941       j=1
 942       for i=1,#listPoints do
 943          if(listPoints[i].type == "bbox") then
 944             output = output .. "dotlabel.llft (btex $\\MeshPoint^{*}_{" .. j .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolorBbox ;"
 945             j=j+1
 946          else
 947             output = output .. "dotlabel.llft (btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolor ;"
 948          end
 949       end
 950       -- mark the point to add
 951       output = output .. "dotlabel.llft (btex $\\NewPoint$ etex,(" .. P.x ..",".. P.y .. ")*u) withcolor \\luameshmpcolorNew;"
 952    elseif(step == "cavity") then
 953       polygon = buildCavity(badTriangles, triangulation)
 954       polyNew = cleanPoly(polygon)
 955       -- remove the bad triangles
 956       for j=1,#badTriangles do
 957          table.remove(triangulation,badTriangles[j]-(j-1))
 958       end
 959       -- draw the triangles
 960       for i=1,#triangulation do
 961          PointI = listPoints[triangulation[i][1]]
 962          PointJ = listPoints[triangulation[i][2]]
 963          PointK = listPoints[triangulation[i][3]]
 964          if(triangulation[i].type == "bbox") then
 965             output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolorBbox;"
 966          else
 967             output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolor;"
 968          end
 969       end
 970       -- fill and draw the cavity
 971       path = ""
 972       for i=1,#polyNew do
 973          PointI = listPoints[polyNew[i]]
 974          path = path .. "(".. PointI.x ..",".. PointI.y ..")*u--"
 975       end
 976       output = output .. "fill " .. path .. "cycle withcolor \\luameshmpcolorBack;"
 977       output = output .. "draw " .. path .. "cycle withcolor \\luameshmpcolorNew  withpen pencircle scaled 1pt;"
 978       -- mark the points of the mesh
 979       j=1
 980       for i=1,#listPoints do
 981          if(listPoints[i].type == "bbox") then
 982             output = output .. "dotlabel.llft (btex $\\MeshPoint^{*}_{" .. j .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolorBbox ;"
 983             j=j+1
 984          else
 985             output = output .. "dotlabel.llft (btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolor ;"
 986          end
 987       end
 988       -- mark the adding point
 989       output = output .. "dotlabel.llft (btex $\\NewPoint$ etex,(" .. P.x ..",".. P.y .. ")*u) withcolor \\luameshmpcolorNew ;"
 990    elseif(step == "newT") then
 991       polygon = buildCavity(badTriangles, triangulation)
 992       polyNew = cleanPoly(polygon)
 993       -- remove the bad triangles
 994       for j=1,#badTriangles do
 995          table.remove(triangulation,badTriangles[j]-(j-1))
 996       end
 997       -- draw the triangle of the triangulation
 998       for i=1,#triangulation do
 999          PointI = listPoints[triangulation[i][1]]
1000          PointJ = listPoints[triangulation[i][2]]
1001          PointK = listPoints[triangulation[i][3]]
1002          if(triangulation[i].type == "bbox") then
1003             output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolorBbox;"
1004          else
1005             output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\luameshmpcolor;"
1006          end
1007       end
1008       -- fill  the cavity
1009       path = ""
1010       for i=1,#polyNew do
1011          PointI = listPoints[polyNew[i]]
1012          path = path .. "(".. PointI.x ..",".. PointI.y ..")*u--"
1013       end
1014       output = output .. "fill " .. path .. "cycle withcolor \\luameshmpcolorBack;"
1015       -- draw the new triangles composed by the edges of the polygon and the added point
1016       for i=1,#polygon do
1017          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;"
1018          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;"
1019          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;"
1020       end
1021       -- mark points
1022       j=1
1023       for i=1,#listPoints do
1024          if(listPoints[i].type == "bbox") then
1025             output = output .. "dotlabel.llft (btex $\\MeshPoint^{*}_{" .. j .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolorBbox ;"
1026             j=j+1
1027          else
1028             output = output .. "dotlabel.llft (btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\luameshmpcolor ;"
1029          end
1030       end
1031       -- mark the added point
1032       output = output .. "dotlabel.llft (btex $\\NewPoint$ etex,(" .. P.x ..",".. P.y .. ")*u) withcolor \\luameshmpcolorNew ;"
1033    end
1034    return output
1035 end
1036
1037
1038 -- build the list of points extern and stop at nbr
1039 function buildListExt(chaine, stop)
1040    listPoints = {}
1041    io.input(chaine) -- open the file
1042    text=io.read("*all")
1043    lines=string.explode(text,"\n+") -- all the lines
1044    for i=1,tonumber(stop) do
1045       xy=string.explode(lines[i]," +")
1046       table.insert(listPoints,{x=tonumber(xy[1]),y=tonumber(xy[2])})
1047    end
1048    xy=string.explode(lines[stop+1]," +")
1049    point={x=tonumber(xy[1]),y=tonumber(xy[2])}
1050    return point, listPoints
1051 end
1052
1053
1054 function TeXOnePointTikZBW(chaine,point,step,scale,mode,bbox,color,colorBack,colorNew,colorCircle,colorBbox)
1055    if(mode=="int") then
1056       Sx,Sy=string.match(point,"%((.+),(.+)%)")
1057       P = {x=Sx, y=Sy}
1058       listPoints = buildList(chaine, mode)
1059    else
1060       -- point is a number
1061       P, listPoints = buildListExt(chaine,tonumber(point))
1062    end
1063    output = TeXaddOnePointTikZ(listPoints,P,step,bbox,color,colorBack,colorNew,colorCircle,colorBbox)
1064    output = "\\noindent\\begin{tikzpicture}[x="..scale..",y="..scale.."]".. output .. "\\end{tikzpicture}"
1065    tex.sprint(output)
1066 end
1067
1068 function TeXOnePointTikZBWinc(chaine,point,beginning, ending,step,scale,mode,bbox,color,colorBack,colorNew,colorCircle,colorBbox)
1069    if(mode=="int") then
1070       Sx,Sy=string.match(point,"%((.+),(.+)%)")
1071       P = {x=Sx, y=Sy}
1072       listPoints = buildList(chaine, mode)
1073    else
1074       -- point is a number
1075       P, listPoints = buildListExt(chaine,tonumber(point))
1076    end
1077    output = TeXaddOnePointTikZ(listPoints,P,step,bbox,color,colorBack,colorNew,colorCircle,colorBbox)
1078    output = "\\noindent\\begin{tikzpicture}[x="..scale..",y="..scale.."]".. beginning..output ..ending.. "\\end{tikzpicture}"
1079    tex.sprint(output)
1080 end
1081
1082 function TeXOnePointMPBW(chaine,point,step,scale,mode,bbox)
1083    if(mode=="int") then
1084       Sx,Sy=string.match(point,"%((.+),(.+)%)")
1085       P = {x=Sx, y=Sy}
1086       listPoints = buildList(chaine, mode)
1087    else
1088       -- point is a number
1089       P, listPoints = buildListExt(chaine,tonumber(point))
1090    end
1091    output = TeXaddOnePointMPBW(listPoints,P,step,bbox)
1092    output = "\\leavevmode\\begin{mplibcode}beginfig(0);u:="..scale..";".. output .. "endfig;\\end{mplibcode}"
1093    tex.sprint(output)
1094 end
1095
1096 function TeXOnePointMPBWinc(chaine,point,beginning,ending,step,scale,mode,bbox)
1097    if(mode=="int") then
1098       Sx,Sy=string.match(point,"%((.+),(.+)%)")
1099       P = {x=Sx, y=Sy}
1100       listPoints = buildList(chaine, mode)
1101    else
1102       -- point is a number
1103       P, listPoints = buildListExt(chaine,tonumber(point))
1104    end
1105    output = TeXaddOnePointMPBW(listPoints,P,step,bbox)
1106    output = "\\begin{mplibcode}u:="..scale..";"..beginning .. output .. ending .. "\\end{mplibcode}"
1107    tex.sprint(output)
1108 end
1109
1110
1111 function split(pString, pPattern)
1112    local Table = {}  -- NOTE: use {n = 0} in Lua-5.0
1113    local fpat = "(.-)" .. pPattern
1114    local last_end = 1
1115    local s, e, cap = pString:find(fpat, 1)
1116    while s do
1117       if s ~= 1 or cap ~= "" then
1118          table.insert(Table,cap)
1119       end
1120       last_end = e+1
1121       s, e, cap = pString:find(fpat, last_end)
1122    end
1123    if last_end <= #pString then
1124       cap = pString:sub(last_end)
1125       table.insert(Table, cap)
1126    end
1127    return Table
1128 end
1129
1130 function readGmsh(file)
1131    io.input(file) -- open the file
1132    text=io.read("*all")
1133    local lines = split(text,"\n+") -- all the lines
1134    listPoints={}
1135    triangulation ={}
1136    boolNodes = false
1137    Jnodes = 0
1138    boolElements = false
1139    Jelements = 0
1140    J=0
1141    for i=1,#lines-J do
1142       if(lines[i+J] == "$EndNodes") then
1143          boolNodes = false
1144          -- go to the next line
1145       end
1146       if(boolNodes) then -- we are in the Nodes environment
1147          xy=split(lines[i+J]," +")
1148          table.insert(listPoints,{x=tonumber(xy[2]),y=tonumber(xy[3])})
1149       end
1150       if(lines[i+J] == "$Nodes") then
1151          boolNodes = true
1152          -- go to the next line
1153          J=J+1
1154       end
1155       if(lines[i+J] == "$EndElements") then
1156          boolElements = false
1157          -- go to the next line
1158       end
1159       if(boolElements) then -- we are in the Nodes environment
1160          xy=split(lines[i+J]," +")
1161          if(tonumber(xy[2]) == 2) then -- if the element is a triangle
1162             nbrTags = xy[3]+1
1163             table.insert(triangulation,{tonumber(xy[2+nbrTags+1]),tonumber(xy[2+nbrTags+2]),tonumber(xy[2+nbrTags+3])})
1164          end
1165       end
1166       if(lines[i+J] == "$Elements") then
1167          boolElements = true
1168          -- go to the next line
1169          J=J+1
1170       end
1171    end
1172    return listPoints, triangulation
1173 end
1174
1175 function drawGmshMP(file,points,scale)
1176    listPoints,triangulation = readGmsh(file)
1177    output = traceMeshMP(listPoints,triangulation,points)
1178    output = "\\leavevmode\\begin{mplibcode}beginfig(0);u:="..scale.. ";" .. output .."endfig;\\end{mplibcode}"
1179    tex.sprint(output)
1180 end
1181
1182 function drawGmshMPinc(file,beginning,ending,points,scale)
1183    listPoints,triangulation = readGmsh(file)
1184    output = traceMeshMP(listPoints,triangulation,points)
1185    output = "\\leavevmode\\begin{mplibcode}u:="..scale..";"..beginning .. output .. ending .. "\\end{mplibcode}"
1186    tex.sprint(output)
1187 end
1188
1189
1190
1191 --
1192 function drawGmshTikZ(file,points,scale,color)
1193    listPoints,triangulation = readGmsh(file)
1194    output = traceMeshTikZ(listPoints, triangulation,points,color,colorBbox)
1195    output = "\\noindent\\begin{tikzpicture}[x=" .. scale .. ",y=" .. scale .."]" .. output .."\\end{tikzpicture}"
1196    tex.sprint(output)
1197 end
1198
1199 --
1200 function drawGmshTikZinc(file,beginning, ending,points,scale,color)
1201    listPoints,triangulation = readGmsh(file)
1202    output = traceMeshTikZ(listPoints, triangulation,points,color,colorBbox)
1203    output = "\\noindent\\begin{tikzpicture}[x=" .. scale .. ",y=" .. scale .."]" ..beginning.. output..ending .."\\end{tikzpicture}"
1204    tex.sprint(output)
1205 end
1206
1207
1208 -- buildVoronoi with MP
1209 function gmshVoronoiMP(file,points,scale,tri,styleD,styleV)
1210    listPoints,triangulation = readGmsh(file)
1211    listVoronoi = buildVoronoi(listPoints, triangulation)
1212    output = traceVoronoiMP(listPoints,triangulation,listVoronoi,points,tri,styleD,styleV)
1213    output = "\\leavevmode\\begin{mplibcode}beginfig(0);u:="..scale.. ";" .. output .."endfig;\\end{mplibcode}"
1214    tex.sprint(output)
1215 end
1216
1217
1218 -- buildVoronoi with TikZ
1219 function gmshVoronoiTikZ(file,points,scale,tri,color,colorVoronoi,styleD,styleV)
1220    listPoints,triangulation = readGmsh(file)
1221    listVoronoi = buildVoronoi(listPoints, triangulation)
1222    output = traceVoronoiTikZ(listPoints,triangulation,listVoronoi,points,tri,color,colorBbox,colorVoronoi,styleD,styleV)
1223    output = "\\noindent\\begin{tikzpicture}[x=" .. scale .. ",y=" .. scale .."]" .. output .."\\end{tikzpicture}"   tex.sprint(output)
1224 end
1225
1226
1227 -- buildVoronoi with MP
1228 function gmshVoronoiMPinc(file,beginning, ending,points,scale,tri,styleD,styleV)
1229    listPoints,triangulation = readGmsh(file)
1230    listVoronoi = buildVoronoi(listPoints, triangulation)
1231    output = traceVoronoiMP(listPoints,triangulation,listVoronoi,points,tri,styleD,styleV)
1232    output = "\\leavevmode\\begin{mplibcode}u:="..scale..";"..beginning .. output .. ending .. "\\end{mplibcode}"
1233    tex.sprint(output)
1234 end
1235
1236
1237 -- buildVoronoi with TikZ
1238 function gmshVoronoiTikZinc(file,beginning, ending,points,scale,tri,color,colorVoronoi,styleD,styleV)
1239    listPoints,triangulation = readGmsh(file)
1240    listVoronoi = buildVoronoi(listPoints, triangulation)
1241    output = traceVoronoiTikZ(listPoints,triangulation,listVoronoi,points,tri,color,colorBbox,colorVoronoi,styleD,styleV)
1242    output = "\\noindent\\begin{tikzpicture}[x=" .. scale .. ",y=" .. scale .."]" ..beginning.. output..ending .."\\end{tikzpicture}"
1243    tex.sprint(output)
1244 end