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