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}
  10    -- the second triangle
  11    triangulation[2] = {lgth+1, lgth+3, lgth+4}
  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          table.insert(triangulation,{polygon[j][1],polygon[j][2],i})
  25       end
  26    end  -- end adding points of the listPoints
  27    -- remove bounding box
  28    if(bbox ~= "bbox") then
  29       triangulation = removeBoundingBox(triangulation,lgth)
  30       table.remove(listPoints,lgth+1)
  31       table.remove(listPoints,lgth+1)
  32       table.remove(listPoints,lgth+1)
  33       table.remove(listPoints,lgth+1)
  34    end
  35    return triangulation
  36 end
  37
  38
  39 function buildBoundingBox(listPoints)
  40    -- listPoints : list of points
  41    -- epsV        : parameter for the distance of the bounding box
  42    local xmin, xmax, ymin, ymax, eps
  43    xmin = 1000
  44    ymin = 1000
  45    xmax = -1000
  46    ymax = -1000
  47    for i=1,#listPoints do
  48       if (listPoints[i].x < xmin) then
  49          xmin = listPoints[i].x
  50       end
  51       if (listPoints[i].x > xmax) then
  52          xmax = listPoints[i].x
  53       end
  54       if (listPoints[i].y < ymin) then
  55          ymin = listPoints[i].y
  56       end
  57       if (listPoints[i].y > ymax) then
  58          ymax = listPoints[i].y
  59       end
  60    end
  61    eps = math.max(math.abs(xmax-xmin),math.abs(ymax-ymin))*0.15
  62    xmin = xmin - eps
  63    xmax = xmax + eps
  64    ymin = ymin - eps
  65    ymax = ymax + eps
  66    -- add points of the bounding box in last positions
  67    table.insert(listPoints,{x=xmin,y=ymin})
  68    table.insert(listPoints,{x=xmin,y=ymax})
  69    table.insert(listPoints,{x=xmax,y=ymax})
  70    table.insert(listPoints,{x=xmax,y=ymin})
  71    return listPoints
  72 end
  73
  74 function removeBoundingBox(triangulation,lgth)
  75    -- build the four bounding box edge
  76    point1 = lgth+1
  77    point2 = lgth+2
  78    point3 = lgth+3
  79    point4 = lgth+4
  80    -- for all triangle
  81    newTriangulation = {}
  82    for i=1,#triangulation do
  83       boolE1 = pointInTriangle(point1,triangulation[i])
  84       boolE2 = pointInTriangle(point2,triangulation[i])
  85       boolE3 = pointInTriangle(point3,triangulation[i])
  86       boolE4 = pointInTriangle(point4,triangulation[i])
  87       if((not boolE1) and (not boolE2) and (not boolE3) and (not boolE4)) then
  88          table.insert(newTriangulation,triangulation[i])
  89       end
  90    end
  91    return newTriangulation
  92 end
  93
  94
  95 function buildBadTriangles(point, triangulation)
  96    badTriangles = {}
  97    for j=1,#triangulation do -- for all triangles
  98       A = listPoints[triangulation[j][1]]
  99       B = listPoints[triangulation[j][2]]
 100       C = listPoints[triangulation[j][3]]
 101       center, radius = circoncircle(A,B,C)
 102       CP = Vector(center,point)
 103       if(VectorNorm(CP)<radius) then -- the point belongs to the circoncirle
 104          table.insert(badTriangles,j)
 105       end
 106    end
 107    return badTriangles
 108 end
 109
 110 -- construction of the cavity composed by the bad triangles around the point to add
 111 function buildCavity(badTriangles, triangulation)
 112    polygon = {}
 113    for j=1,#badTriangles do -- for all bad triangles
 114       ind = badTriangles[j]
 115       for k=1,3 do -- for all edges
 116          edge = {triangulation[ind][k],triangulation[ind][k%3+1]}
 117          edgeBord = false
 118          for l = 1,#badTriangles do -- for all badtriangles
 119             badInd = badTriangles[l]
 120             if(badInd ~= ind) then -- if not the current one
 121                edgeBord = edgeBord or edgeInTriangle(edge,triangulation[badInd])
 122             end
 123          end --
 124          -- if the edge does not belong to another bad triangle
 125          if(edgeBord == false) then
 126             -- insert the edge to the cavity
 127             table.insert(polygon,edge)
 128          end
 129       end --
 130    end --
 131    return polygon
 132 end
 133
 134 function edgeInTriangle(e,t)
 135    in1 = false
 136    in2 = false
 137    for i=1,3 do
 138       if e[1] == t[i] then
 139          in1 = true
 140       end
 141       if e[2] == t[i] then
 142          in2 = true
 143       end
 144    end
 145    out = (in1 and in2)
 146    return out
 147 end
 148
 149 function pointInTriangle(e,t)
 150    in1 = false
 151    for i=1,3 do
 152       if e == t[i] then
 153          in1 = true
 154       end
 155    end
 156    return in1
 157 end
 158
 159
 160 function Vector(A,B)
 161    local out = {x = B.x - A.x, y = B.y - A.y}
 162    return out
 163 end
 164
 165 function VectorNorm(NP)
 166    return math.sqrt(NP.x*NP.x +NP.y*NP.y)
 167 end
 168
 169 -- circoncircle
 170 function circoncircle(M, N, P)
 171    -- Compute center and radius of the circoncircle of the triangle M N P
 172
 173    -- return : (center [Point],radius [float])
 174
 175    local MN = Vector(M,N)
 176    local NP = Vector(N,P)
 177    local PM = Vector(P,M)
 178    m = VectorNorm(NP)  -- |NP|
 179    n = VectorNorm(PM)  -- |PM|
 180    p = VectorNorm(MN)  -- |MN|
 181
 182    d = (m + n + p) * (-m + n + p) * (m - n + p) * (m + n - p)
 183    if d > 0 then
 184       rad = m * n * p / math.sqrt(d)
 185    else
 186       rad = 0
 187    end
 188    d = -2 * (M.x * NP.y + N.x * PM.y + P.x * MN.y)
 189    O = {x=0, y=0}
 190    OM = Vector(O, M)
 191    ON = Vector(O, N)
 192    OP = Vector(O, P)
 193    om2 = math.pow(VectorNorm(OM),2)  -- |OM|**2
 194    on2 = math.pow(VectorNorm(ON),2)  -- |ON|**2
 195    op2 = math.pow(VectorNorm(OP),2)  -- |OP|**2
 196    x0 = -(om2 * NP.y + on2 * PM.y + op2 * MN.y) / d
 197    y0 = (om2 * NP.x + on2 * PM.x + op2 * MN.x) / d
 198    if d == 0 then
 199       Out = {nil, nil}
 200    else
 201       Out = {x=x0, y=y0}
 202    end
 203    return Out, rad  -- (center [Point], R [float])
 204 end
 205
 206
 207 -------------------------- TeX
 208 -- build the list of points
 209 function buildList(chaine, mode)
 210    -- if mode = int : the list is given in the chaine string (x1,y1);(x2,y2);...;(xn,yn)
 211    -- if mode = ext : the list is given in a file line by line with space separation
 212    listPoints = {}
 213    if mode == "int" then
 214       local points = string.explode(chaine, ";")
 215       local lgth=#points
 216       for i=1,lgth do
 217          Sx,Sy=string.match(points[i],"%((.+),(.+)%)")
 218          listPoints[i]={x=tonumber(Sx),y=tonumber(Sy)}
 219       end
 220    elseif mode == "ext" then
 221       io.input(chaine) -- open the file
 222       text=io.read("*all")
 223       lines=string.explode(text,"\n+") -- all the lines
 224       tablePoints={}
 225       for i=1,#lines do
 226          xy=string.explode(lines[i]," +")
 227          listPoints[i]={x=tonumber(xy[1]),y=tonumber(xy[2])}
 228       end
 229    else
 230       print("Non existing mode")
 231    end
 232    return listPoints
 233 end
 234
 235
 236 --
 237 function rectangleList(a,b,nbrA,nbrB)
 238    stepA = a/nbrA
 239    stepB = b/nbrB
 240    listPoints = {}
 241    k=1
 242    for i=1,(nbrA+1) do
 243       for j=1,(nbrB+1) do
 244          listPoints[k] = {x = (i-1)*stepA, y=(j-1)*stepB}
 245          k=k+1
 246       end
 247    end
 248    return listPoints
 249 end
 250
 251 -- trace a triangulation with TikZ
 252 function traceMeshTikZ(listPoints, triangulation,points,color)
 253    output = ""
 254    for i=1,#triangulation do
 255       PointI = listPoints[triangulation[i][1]]
 256       PointJ = listPoints[triangulation[i][2]]
 257       PointK = listPoints[triangulation[i][3]]
 258       output = output .. "\\draw[color="..color.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
 259    end
 260    if(points=="points") then
 261       for i=1,#listPoints do
 262          output = output .. "\\draw[color=".."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
 263       end
 264    end
 265    return output
 266 end
 267
 268
 269 -- trace a triangulation with MP
 270 function traceMeshMP(listPoints, triangulation,points,color)
 271    output = "";
 272    output = output .. " pair MeshPoints[];"
 273    for i=1,#listPoints do
 274       output = output .. "MeshPoints[".. i .. "] = (" .. listPoints[i].x .. "," .. listPoints[i].y .. ");"
 275    end
 276
 277    for i=1,#triangulation do
 278       PointI = listPoints[triangulation[i][1]]
 279       PointJ = listPoints[triangulation[i][2]]
 280       PointK = listPoints[triangulation[i][3]]
 281       output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\mpcolor{"..color.. "};"
 282    end
 283    if(points=="points") then
 284       for i=1,#listPoints do
 285          output = output .. "dotlabel.llft( btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u) withcolor \\mpcolor{"..color.. "};"
 286       end
 287    end
 288    return output
 289 end
 290
 291
 292
 293 -- buildMesh with TikZ
 294 function buildMeshTikZ(chaine,mode,points,bbox,full,scale,color)
 295    listPoints = buildList(chaine, mode)
 296    triangulation = BowyerWatson(listPoints,bbox)
 297    output = traceMeshTikZ(listPoints, triangulation,points,color)
 298    if(full=="full") then
 299       output = "\\noindent\\begin{tikzpicture}[x=" .. scale .. ",y=" .. scale .."]" .. output .."\\end{tikzpicture}"
 300    end
 301    tex.sprint(output)
 302 end
 303
 304
 305 -- buildMesh with MP
 306 function buildMeshMP(chaine,mode,points,bbox,full,scale,color)
 307    listPoints = buildList(chaine, mode)
 308    triangulation = BowyerWatson(listPoints,bbox)
 309    output = traceMeshMP(listPoints, triangulation,points,color)
 310    if(full=="full") then
 311       output = "\\leavevmode\\begin{mplibcode}beginfig(0);u:="..scale.. ";" .. output .."endfig;\\end{mplibcode}"
 312    else
 313       output = "u:="..scale.. ";".. output
 314    end
 315    tex.sprint(output)
 316 end
 317
 318 -- buildMesh
 319 function buildRect(largeur,a,b,nbrA, nbrB)
 320    listPoints = rectangleList(a,b,nbrA,nbrB)
 321    triangulation = BowyerWatson(listPoints,"none")
 322    traceTikZ(listPoints, triangulation,largeur,"none")
 323 end
 324
 325
 326 -- function give a real polygon without repeting points
 327 function cleanPoly(polygon)
 328    polyNew = {}
 329    e1 = polygon[1][1]
 330    e2 = polygon[1][2]
 331    table.insert(polyNew, e1)
 332    table.insert(polyNew, e2)
 333    j = 2
 334    for i=2,#polygon do
 335       bool1 = (polygon[i][1] == polyNew[j])
 336       bool2 = (polygon[i][2] == polyNew[j])
 337       if(bool1 or bool2) then -- the edge has a common point with polyNew[j]
 338          if(not bool1) then
 339             table.insert(polyNew, polygon[i][1])
 340             j = j+1
 341          elseif(not bool2) then
 342             table.insert(polyNew, polygon[i][2])
 343             j = j+1
 344          end
 345       end
 346    end
 347    return polyNew
 348 end
 349
 350 --
 351 function TeXaddOnePointTikZ(chaine,point,step,color,colorBack, colorNew, colorCircle)
 352    Sx,Sy=string.match(point,"%((.+),(.+)%)")
 353    P = {x=Sx, y=Sy}
 354    output = ""
 355    listPoints = buildList(chaine, "int")
 356    -- build the triangulation
 357    triangulation = BowyerWatson(listPoints,"none")
 358    badTriangles = buildBadTriangles(P,triangulation)
 359    if(step == "badT") then
 360       -- draw all triangle
 361       for i=1,#triangulation do
 362          PointI = listPoints[triangulation[i][1]]
 363          PointJ = listPoints[triangulation[i][2]]
 364          PointK = listPoints[triangulation[i][3]]
 365          output = output .. "\\draw[color="..color.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
 366       end
 367       -- draw and fill the bad triangle
 368       for i=1,#badTriangles do
 369          PointI = listPoints[triangulation[badTriangles[i]][1]]
 370          PointJ = listPoints[triangulation[badTriangles[i]][2]]
 371          PointK = listPoints[triangulation[badTriangles[i]][3]]
 372          output = output .. "\\draw[fill="..colorBack.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
 373       end
 374       -- draw the circoncircle
 375       for i=1,#badTriangles do
 376          PointI = listPoints[triangulation[badTriangles[i]][1]]
 377          PointJ = listPoints[triangulation[badTriangles[i]][2]]
 378          PointK = listPoints[triangulation[badTriangles[i]][3]]
 379          center, radius = circoncircle(PointI, PointJ, PointK)
 380          output = output .. "\\draw[dashed, color="..colorCircle.."] ("..center.x .. "," .. center.y .. ") circle ("..radius ..");"
 381       end
 382       -- mark the points
 383       for i=1,#listPoints do
 384          output = output .. "\\draw[color ="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
 385       end
 386       -- mark the point to add
 387       output = output .. "\\draw[color="..colorNew.."] (" .. P.x ..",".. P.y .. ") node {$\\bullet$} node[anchor=north east] {$\\NewPoint$};"
 388    elseif(step == "cavity") then
 389       polygon = buildCavity(badTriangles, triangulation)
 390       polyNew = cleanPoly(polygon)
 391       -- remove the bad triangles
 392       for j=1,#badTriangles do
 393          table.remove(triangulation,badTriangles[j]-(j-1))
 394       end
 395       -- draw the triangles
 396       for i=1,#triangulation do
 397          PointI = listPoints[triangulation[i][1]]
 398          PointJ = listPoints[triangulation[i][2]]
 399          PointK = listPoints[triangulation[i][3]]
 400          output = output .. "\\draw[color="..color.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
 401       end
 402       -- fill and draw the cavity
 403       path = ""
 404       for i=1,#polyNew do
 405          PointI = listPoints[polyNew[i]]
 406          path = path .. "(".. PointI.x ..",".. PointI.y ..")--"
 407       end
 408       output = output .. "\\draw[color="..colorNew..",fill ="..colorBack..", thick] " .. path .. "cycle;"
 409       -- mark the points of the mesh
 410       for i=1,#listPoints do
 411          output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
 412       end
 413       -- mark the adding point
 414       output = output .. "\\draw[color="..colorNew.."] (" .. P.x ..",".. P.y .. ") node {$\\bullet$} node[anchor=north east] {$\\NewPoint$};"
 415    elseif(step == "newT") then
 416       polygon = buildCavity(badTriangles, triangulation)
 417       polyNew = cleanPoly(polygon)
 418       -- remove the bad triangles
 419       for j=1,#badTriangles do
 420          table.remove(triangulation,badTriangles[j]-(j-1))
 421       end
 422       -- draw the triangle of the triangulation
 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          output = output .. "\\draw[color ="..color.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;"
 428       end
 429       -- fill and draw the cavity
 430       path = ""
 431       for i=1,#polyNew do
 432          PointI = listPoints[polyNew[i]]
 433          path = path .. "(".. PointI.x ..",".. PointI.y ..")--"
 434       end
 435       output = output .. "\\draw[color="..colorNew..",fill ="..colorBack..", thick] " .. path .. "cycle;"
 436       -- draw the new triangles composed by the edges of the polygon and the added point
 437       for i=1,#polygon do
 438          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 ..");"
 439          output = output .. "\\draw[color="..colorNew..", thick]".."(".. listPoints[polygon[i][1]].x .. "," .. listPoints[polygon[i][1]].y .. ") -- (" .. P.x .. "," .. P.y ..");"
 440          output = output .. "\\draw[color="..colorNew..", thick]".."(".. listPoints[polygon[i][2]].x .. "," .. listPoints[polygon[i][2]].y .. ") -- (" .. P.x .. "," .. P.y ..");"
 441       end
 442       -- mark points
 443       for i=1,#listPoints do
 444          output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};"
 445       end
 446       -- mark the added point
 447       output = output .. "\\draw[color="..colorNew.."] (" .. P.x ..",".. P.y .. ") node {$\\bullet$} node[anchor=north east] {$\\NewPoint$};"
 448    end
 449    return output
 450 end
 451
 452 function TeXaddOnePointMP(listPoints,P,step,color,colorBack, colorNew, colorCircle,bbox)
 453    output = "";
 454    output = output .. "pair MeshPoints[];"
 455    -- build the triangulation
 456    triangulation = BowyerWatson(listPoints,bbox)
 457    badTriangles = buildBadTriangles(P,triangulation)
 458    for i=1,#listPoints do
 459       output = output .. "MeshPoints[".. i .. "] = (" .. listPoints[i].x .. "," .. listPoints[i].y .. ");"
 460    end
 461    if(step == "badT") then
 462       -- draw all triangle
 463       for i=1,#triangulation do
 464          PointI = listPoints[triangulation[i][1]]
 465          PointJ = listPoints[triangulation[i][2]]
 466          PointK = listPoints[triangulation[i][3]]
 467          output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\mpcolor{" .. color .."};"
 468       end
 469       -- draw and fill the bad triangle
 470       for i=1,#badTriangles do
 471          PointI = listPoints[triangulation[badTriangles[i]][1]]
 472          PointJ = listPoints[triangulation[badTriangles[i]][2]]
 473          PointK = listPoints[triangulation[badTriangles[i]][3]]
 474          output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\mpcolor{" .. color .."};"
 475          output = output .. "fill (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\mpcolor{" .. colorBack .."};"
 476       end
 477       -- draw the circoncircle
 478       for i=1,#badTriangles do
 479          PointI = listPoints[triangulation[badTriangles[i]][1]]
 480          PointJ = listPoints[triangulation[badTriangles[i]][2]]
 481          PointK = listPoints[triangulation[badTriangles[i]][3]]
 482          center, radius = circoncircle(PointI, PointJ, PointK)
 483          output = output .. "draw fullcircle scaled ("..radius .."*2u) shifted ("..center.x .. "*u," .. center.y .. "*u) dashed evenly withcolor \\mpcolor{" .. colorCircle .."};"
 484       end
 485       -- mark the points
 486       for i=1,#listPoints do
 487          output = output .. "dotlabel.llft (btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\mpcolor{" .. color .."};"
 488       end
 489       -- mark the point to add
 490       output = output .. "dotlabel.llft (btex $\\NewPoint$ etex,(" .. P.x ..",".. P.y .. ")*u) withcolor \\mpcolor{" .. colorNew .."};"
 491    elseif(step == "cavity") then
 492       polygon = buildCavity(badTriangles, triangulation)
 493       polyNew = cleanPoly(polygon)
 494       -- remove the bad triangles
 495       for j=1,#badTriangles do
 496          table.remove(triangulation,badTriangles[j]-(j-1))
 497       end
 498       -- draw the triangles
 499       for i=1,#triangulation do
 500          PointI = listPoints[triangulation[i][1]]
 501          PointJ = listPoints[triangulation[i][2]]
 502          PointK = listPoints[triangulation[i][3]]
 503          output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\mpcolor{" .. color .."};"
 504       end
 505       -- fill and draw the cavity
 506       path = ""
 507       for i=1,#polyNew do
 508          PointI = listPoints[polyNew[i]]
 509          path = path .. "(".. PointI.x ..",".. PointI.y ..")*u--"
 510       end
 511       output = output .. "fill " .. path .. "cycle withcolor \\mpcolor{" .. colorBack .."};"
 512       output = output .. "draw " .. path .. "cycle withcolor \\mpcolor{" .. colorNew .."} withpen pencircle scaled 1pt;"
 513       -- mark the points of the mesh
 514       for i=1,#listPoints do
 515          output = output .. "dotlabel.llft (btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\mpcolor{" .. color .."};"
 516       end
 517       -- mark the adding point
 518       output = output .. "dotlabel.llft (btex $\\NewPoint$ etex,(" .. P.x ..",".. P.y .. ")*u) withcolor \\mpcolor{" .. colorNew .."};"
 519    elseif(step == "newT") then
 520       polygon = buildCavity(badTriangles, triangulation)
 521       polyNew = cleanPoly(polygon)
 522       -- remove the bad triangles
 523       for j=1,#badTriangles do
 524          table.remove(triangulation,badTriangles[j]-(j-1))
 525       end
 526       -- draw the triangle of the triangulation
 527       for i=1,#triangulation do
 528          PointI = listPoints[triangulation[i][1]]
 529          PointJ = listPoints[triangulation[i][2]]
 530          PointK = listPoints[triangulation[i][3]]
 531          output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\mpcolor{" .. color .."};"
 532       end
 533       -- fill  the cavity
 534       path = ""
 535       for i=1,#polyNew do
 536          PointI = listPoints[polyNew[i]]
 537          path = path .. "(".. PointI.x ..",".. PointI.y ..")*u--"
 538       end
 539       output = output .. "fill " .. path .. "cycle withcolor \\mpcolor{" .. colorBack .."};"
 540       -- draw the new triangles composed by the edges of the polygon and the added point
 541       for i=1,#polygon do
 542          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 \\mpcolor{" .. colorNew .."} withpen pencircle scaled 1pt;"
 543          output = output .. "draw".."(".. listPoints[polygon[i][1]].x .. "," .. listPoints[polygon[i][1]].y .. ")*u -- (" .. P.x .. "," .. P.y ..")*u withcolor \\mpcolor{" .. colorNew .."} withpen pencircle scaled 1pt;"
 544          output = output .. "draw".."(".. listPoints[polygon[i][2]].x .. "," .. listPoints[polygon[i][2]].y .. ")*u -- (" .. P.x .. "," .. P.y ..")*u withcolor \\mpcolor{" .. colorNew .."} withpen pencircle scaled 1pt;"
 545       end
 546       -- mark points
 547       for i=1,#listPoints do
 548          output = output .. "dotlabel.llft (btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\mpcolor{" .. color .."};"
 549       end
 550       -- mark the added point
 551       output = output .. "dotlabel.llft (btex $\\NewPoint$ etex,(" .. P.x ..",".. P.y .. ")*u) withcolor \\mpcolor{ " .. colorNew .."};"
 552    end
 553    return output
 554 end
 555
 556
 557 -- build the list of points extern and stop at nbr
 558 function buildListExt(chaine, stop)
 559    listPoints = {}
 560    io.input(chaine) -- open the file
 561    text=io.read("*all")
 562    lines=string.explode(text,"\n+") -- all the lines
 563    for i=1,tonumber(stop) do
 564       xy=string.explode(lines[i]," +")
 565       table.insert(listPoints,{x=tonumber(xy[1]),y=tonumber(xy[2])})
 566    end
 567    xy=string.explode(lines[stop+1]," +")
 568    point={x=tonumber(xy[1]),y=tonumber(xy[2])}
 569    return point, listPoints
 570 end
 571
 572
 573 function TeXOnePointTikZ(chaine,point,step,color,colorBack,colorNew,colorCircle,scale)
 574    output = TeXaddOnePointTikZ(chaine,point,step,color,colorBack,colorNew,colorCircle)
 575    output = "\\noindent\\begin{tikzpicture}[x="..scale..",y="..scale.."]".. output .. "\\end{tikzpicture}"
 576    tex.sprint(output)
 577 end
 578
 579 function TeXOnePointMP(chaine,point,step,color,colorBack,colorNew,colorCircle,scale,mode,picture,bbox)
 580    if(mode=="int") then
 581       Sx,Sy=string.match(point,"%((.+),(.+)%)")
 582       P = {x=Sx, y=Sy}
 583       listPoints = buildList(chaine, mode)
 584    else
 585       -- point is a number
 586       P, listPoints = buildListExt(chaine,tonumber(point))
 587    end
 588    output = TeXaddOnePointMP(listPoints,P,step,color,colorBack,colorNew,colorCircle,bbox)
 589    if(picture=="full") then
 590       output = "\\leavevmode\\begin{mplibcode}beginfig(0);u:="..scale..";".. output .. "endfig;\\end{mplibcode}"
 591    else
 592       output = "u:="..scale..";".. output
 593    end
 594    tex.sprint(output)
 595 end