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