X-Git-Url: https://melusine.eu.org/syracuse/G/git/?p=delaunay.git;a=blobdiff_plain;f=luamesh.lua;h=a1023aace9cb57f68b90d8fc0e32d05adfd5a7ce;hp=f165fe69389ca404a6720e4226d78913b402c285;hb=HEAD;hpb=99c7eb1c5fdd23d7617d3817806af416e4614d15 diff --git a/luamesh.lua b/luamesh.lua index f165fe6..a1023aa 100644 --- a/luamesh.lua +++ b/luamesh.lua @@ -1,18 +1,29 @@ +require "luamesh-polygon" +require "luamesh-tex" + +local function shallowCopy(original) + local copy = {} + for key, value in pairs(original) do + copy[key] = value + end + return copy +end + -- Bowyer and Watson algorithm -- Delaunay meshing function BowyerWatson (listPoints,bbox) local triangulation = {} - local lgth = #listPoints + lgth = #listPoints -- add four points to listPoints to have a bounding box listPoints = buildBoundingBox(listPoints) -- the first triangle - triangulation[1] = {lgth+1, lgth+2, lgth+3} + triangulation[1] = {lgth+1, lgth+2, lgth+3,type="bbox"} -- the second triangle - triangulation[2] = {lgth+1, lgth+3, lgth+4} + triangulation[2] = {lgth+1, lgth+3, lgth+4,type="bbox"} -- add points one by one for i=1,lgth do -- find the triangles which the circumcircle contained the point to add - badTriangles = buildBadTriangles(listPoints[i],triangulation) + badTriangles = buildBadTriangles(listPoints[i],triangulation,listPoints) -- build the polygon of the cavity containing the point to add polygon = buildCavity(badTriangles, triangulation) -- remove the bad triangles @@ -21,7 +32,11 @@ function BowyerWatson (listPoints,bbox) end -- build the new triangles and add them to triangulation for j=1,#polygon do - table.insert(triangulation,{polygon[j][1],polygon[j][2],i}) + if((polygon[j][1]>lgth) or (polygon[j][2]>lgth) or (i>lgth)) then + table.insert(triangulation,{polygon[j][1],polygon[j][2],i,type="bbox"}) + else + table.insert(triangulation,{polygon[j][1],polygon[j][2],i,type="in"}) + end end end -- end adding points of the listPoints -- remove bounding box @@ -64,13 +79,43 @@ function buildBoundingBox(listPoints) ymin = ymin - eps ymax = ymax + eps -- add points of the bounding box in last positions - table.insert(listPoints,{x=xmin,y=ymin}) - table.insert(listPoints,{x=xmin,y=ymax}) - table.insert(listPoints,{x=xmax,y=ymax}) - table.insert(listPoints,{x=xmax,y=ymin}) + table.insert(listPoints,{x=xmin,y=ymin,type="bbox"}) + table.insert(listPoints,{x=xmin,y=ymax,type="bbox"}) + table.insert(listPoints,{x=xmax,y=ymax,type="bbox"}) + table.insert(listPoints,{x=xmax,y=ymin,type="bbox"}) return listPoints end +function BoundingBox(listPoints) + -- listPoints : list of points + -- epsV : parameter for the distance of the bounding box + local xmin, xmax, ymin, ymax, eps + xmin = 1000 + ymin = 1000 + xmax = -1000 + ymax = -1000 + for i=1,#listPoints do + if (listPoints[i].x < xmin) then + xmin = listPoints[i].x + end + if (listPoints[i].x > xmax) then + xmax = listPoints[i].x + end + if (listPoints[i].y < ymin) then + ymin = listPoints[i].y + end + if (listPoints[i].y > ymax) then + ymax = listPoints[i].y + end + end + eps = math.max(math.abs(xmax-xmin),math.abs(ymax-ymin))*0.15 + xmin = xmin - eps + xmax = xmax + eps + ymin = ymin - eps + ymax = ymax + eps + return xmin, xmax, ymin, ymax +end + function removeBoundingBox(triangulation,lgth) -- build the four bounding box edge point1 = lgth+1 @@ -78,7 +123,7 @@ function removeBoundingBox(triangulation,lgth) point3 = lgth+3 point4 = lgth+4 -- for all triangle - newTriangulation = {} + local newTriangulation = {} for i=1,#triangulation do boolE1 = pointInTriangle(point1,triangulation[i]) boolE2 = pointInTriangle(point2,triangulation[i]) @@ -92,8 +137,8 @@ function removeBoundingBox(triangulation,lgth) end -function buildBadTriangles(point, triangulation) - badTriangles = {} +function buildBadTriangles(point, triangulation,listPoints) + local badTriangles = {} for j=1,#triangulation do -- for all triangles A = listPoints[triangulation[j][1]] B = listPoints[triangulation[j][2]] @@ -109,7 +154,7 @@ end -- construction of the cavity composed by the bad triangles around the point to add function buildCavity(badTriangles, triangulation) - polygon = {} + local polygon = {} for j=1,#badTriangles do -- for all bad triangles ind = badTriangles[j] for k=1,3 do -- for all edges @@ -203,13 +248,87 @@ function circoncircle(M, N, P) return Out, rad -- (center [Point], R [float]) end +-- compute the list of the circumcircle of a triangulation +function listCircumCenter(listPoints,triangulation) + local list = {} + for j=1,#triangulation do + A = listPoints[triangulation[j][1]] + B = listPoints[triangulation[j][2]] + C = listPoints[triangulation[j][3]] + center, radius = circoncircle(A,B,C) + table.insert(list,{x=center.x,y=center.y,r=radius}) + end + return list +end + +-- find the three neighbour triangles of T +function findNeighbour(T,i,triangulation) + -- T : triangle + -- i : index of T in triangualation + -- triangulation + + list = {} + -- define the three edge + e1 = {T[1],T[2]} + e2 = {T[2],T[3]} + e3 = {T[3],T[1]} + for j=1,#triangulation do + if j~= i then + if(edgeInTriangle(e1,triangulation[j])) then + table.insert(list,j) + end + if(edgeInTriangle(e2,triangulation[j])) then + table.insert(list,j) + end + if(edgeInTriangle(e3,triangulation[j])) then + table.insert(list,j) + end + end + end + return list +end + +-- test if edge are the same (reverse) +function equalEdge(e1,e2) + if(((e1[1] == e2[1]) and (e1[2] == e2[2])) or ((e1[1] == e2[2]) and (e1[2] == e2[1]))) then + return true + else + return false + end +end + +-- test if the edge belongs to the list +function edgeInList(e,listE) + output = false + for i=1,#listE do + if(equalEdge(e,listE[i])) then + output = true + end + end + return output +end + +-- build the edges of the Voronoi diagram with a given triangulation +function buildVoronoi(listPoints, triangulation) + local listCircumCircle = listCircumCenter(listPoints, triangulation) + local listVoronoi = {} + for i=1,#listCircumCircle do + listN = findNeighbour(triangulation[i],i,triangulation) + for j=1,#listN do + edge = {i,listN[j]} + if( not edgeInList(edge, listVoronoi)) then + table.insert(listVoronoi, edge) + end + end + end + return listVoronoi +end --------------------------- TeX -- build the list of points function buildList(chaine, mode) -- if mode = int : the list is given in the chaine string (x1,y1);(x2,y2);...;(xn,yn) -- if mode = ext : the list is given in a file line by line with space separation - listPoints = {} + local listPoints = {} if mode == "int" then local points = string.explode(chaine, ";") local lgth=#points @@ -233,330 +352,122 @@ function buildList(chaine, mode) end --- -function rectangleList(a,b,nbrA,nbrB) - stepA = a/nbrA - stepB = b/nbrB - listPoints = {} - k=1 - for i=1,(nbrA+1) do - for j=1,(nbrB+1) do - listPoints[k] = {x = (i-1)*stepA, y=(j-1)*stepB} - k=k+1 - end - end - return listPoints -end - --- trace a triangulation with TikZ -function traceMeshTikZ(listPoints, triangulation,points,color) - output = "" - for i=1,#triangulation do - PointI = listPoints[triangulation[i][1]] - PointJ = listPoints[triangulation[i][2]] - PointK = listPoints[triangulation[i][3]] - output = output .. "\\draw[color="..color.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;" - end - if(points=="points") then - for i=1,#listPoints do - output = output .. "\\draw[color=".."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};" +-- function to add points on a polygon to respect +-- the size of unit mesh +function addPointsPolygon(polygon,h) + local newPolygon = shallowCopy(polygon) + k=0 -- to follow in the newPolygon + for i=1,#polygon do + k = k+1 + ip = (i)%(#polygon)+1 + dist = math.sqrt(math.pow(polygon[i].x-polygon[ip].x,2) + math.pow(polygon[i].y-polygon[ip].y,2)) + -- if the distance between two ponits of the polygon is greater than 1.5*h + if(dist>=2*h) then + n = math.floor(dist/h) + step = dist/(n+1) + for j=1,n do + a = {x=polygon[i].x+j*step*(polygon[ip].x-polygon[i].x)/dist,y=polygon[i].y+j*step*(polygon[ip].y-polygon[i].y)/dist} + table.insert(newPolygon,k+j,a) + end + k=k+n end end - return output + return newPolygon end +-- function to build a gridpoints from the bounding box +-- with a prescribed +function buildGrid(listPoints,h,random) + -- listPoints : list of the points of the polygon, ordered + -- h : parameter for the grid + xmin, xmax, ymin, ymax = BoundingBox(listPoints) --- trace a triangulation with MP -function traceMeshMP(listPoints, triangulation,points,color) - output = ""; - output = output .. " pair MeshPoints[];" - for i=1,#listPoints do - output = output .. "MeshPoints[".. i .. "] = (" .. listPoints[i].x .. "," .. listPoints[i].y .. ");" - end - - for i=1,#triangulation do - PointI = listPoints[triangulation[i][1]] - PointJ = listPoints[triangulation[i][2]] - PointK = listPoints[triangulation[i][3]] - output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\mpcolor{"..color.. "};" - end - if(points=="points") then - for i=1,#listPoints do - output = output .. "dotlabel.llft( btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u) withcolor \\mpcolor{"..color.. "};" - end - end - return output + local grid = rectangleList(xmin,xmax,ymin,ymax,h,random) + return grid end - - --- buildMesh with TikZ -function buildMeshTikZ(chaine,mode,points,bbox,full,scale,color) - listPoints = buildList(chaine, mode) - triangulation = BowyerWatson(listPoints,bbox) - output = traceMeshTikZ(listPoints, triangulation,points,color) - if(full=="full") then - output = "\\noindent\\begin{tikzpicture}[x=" .. scale .. ",y=" .. scale .."]" .. output .."\\end{tikzpicture}" +-- function to build the list of points in the rectangle +function rectangleList(xmin,xmax,ymin,ymax,h,random) + -- for the random + math.randomseed( os.time() ) + nbrX = math.floor(math.abs(xmax-xmin)/h) + nbrY = math.floor(math.abs(ymax-ymin)/h) + local listPoints = {} + k=1 + for i=1,(nbrX+1) do + for j=1,(nbrY+1) do + rd = math.random() + if(random=="perturb") then + fact = 0.3*h + --print(fact) + else + fact = 0.0 + end + listPoints[k] = {x = xmin+(i-1)*h+rd*fact, y=ymin+(j-1)*h+rd*fact} + k=k+1 + end end - tex.sprint(output) + return listPoints end --- buildMesh with MP -function buildMeshMP(chaine,mode,points,bbox,full,scale,color) - listPoints = buildList(chaine, mode) - triangulation = BowyerWatson(listPoints,bbox) - output = traceMeshMP(listPoints, triangulation,points,color) - if(full=="full") then - output = "\\leavevmode\\begin{mplibcode}beginfig(0);u:="..scale.. ";" .. output .."endfig;\\end{mplibcode}" - else - output = "u:="..scale.. ";".. output +-- function to add points from a grid to the interior of a polygon +function addGridPoints(polygon, grid,h) + local listPoints = shallowCopy(polygon) + k = #polygon + for i=1, #grid do + --print(grid[i].x,grid[i].y) + --print(isInside(polygon,grid[i])) + if(isInside(polygon,grid[i],h)) then + k=k+1 + listPoints[k] = grid[i] + end end - tex.sprint(output) + return listPoints end --- buildMesh -function buildRect(largeur,a,b,nbrA, nbrB) - listPoints = rectangleList(a,b,nbrA,nbrB) - triangulation = BowyerWatson(listPoints,"none") - traceTikZ(listPoints, triangulation,largeur,"none") -end -- function give a real polygon without repeting points function cleanPoly(polygon) - polyNew = {} - e1 = polygon[1][1] - e2 = polygon[1][2] + local polyNew = {} + local polyCopy = shallowCopy(polygon) + e1 = polyCopy[1][1] + e2 = polyCopy[1][2] table.insert(polyNew, e1) table.insert(polyNew, e2) + table.remove(polyCopy,1) j = 2 - for i=2,#polygon do - bool1 = (polygon[i][1] == polyNew[j]) - bool2 = (polygon[i][2] == polyNew[j]) - if(bool1 or bool2) then -- the edge has a common point with polyNew[j] - if(not bool1) then - table.insert(polyNew, polygon[i][1]) - j = j+1 - elseif(not bool2) then - table.insert(polyNew, polygon[i][2]) - j = j+1 + while #polyCopy>1 do + i=1 + find = false + while (i<=#polyCopy and find==false) do + bool1 = (polyCopy[i][1] == polyNew[j]) + bool2 = (polyCopy[i][2] == polyNew[j]) + if(bool1 or bool2) then -- the edge has a common point with polyNew[j] + if(not bool1) then + table.insert(polyNew, polyCopy[i][1]) + find = true + table.remove(polyCopy,i) + j = j+1 + elseif(not bool2) then + table.insert(polyNew, polyCopy[i][2]) + find = true + table.remove(polyCopy,i) + j = j+1 + end end + i=i+1 end end return polyNew end --- -function TeXaddOnePointTikZ(chaine,point,step,color,colorBack, colorNew, colorCircle) - Sx,Sy=string.match(point,"%((.+),(.+)%)") - P = {x=Sx, y=Sy} - output = "" - listPoints = buildList(chaine, "int") - -- build the triangulation - triangulation = BowyerWatson(listPoints,"none") - badTriangles = buildBadTriangles(P,triangulation) - if(step == "badT") then - -- draw all triangle - for i=1,#triangulation do - PointI = listPoints[triangulation[i][1]] - PointJ = listPoints[triangulation[i][2]] - PointK = listPoints[triangulation[i][3]] - output = output .. "\\draw[color="..color.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;" - end - -- draw and fill the bad triangle - for i=1,#badTriangles do - PointI = listPoints[triangulation[badTriangles[i]][1]] - PointJ = listPoints[triangulation[badTriangles[i]][2]] - PointK = listPoints[triangulation[badTriangles[i]][3]] - output = output .. "\\draw[fill="..colorBack.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;" - end - -- draw the circoncircle - for i=1,#badTriangles do - PointI = listPoints[triangulation[badTriangles[i]][1]] - PointJ = listPoints[triangulation[badTriangles[i]][2]] - PointK = listPoints[triangulation[badTriangles[i]][3]] - center, radius = circoncircle(PointI, PointJ, PointK) - output = output .. "\\draw[dashed, color="..colorCircle.."] ("..center.x .. "," .. center.y .. ") circle ("..radius ..");" - end - -- mark the points - for i=1,#listPoints do - output = output .. "\\draw[color ="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};" - end - -- mark the point to add - output = output .. "\\draw[color="..colorNew.."] (" .. P.x ..",".. P.y .. ") node {$\\bullet$} node[anchor=north east] {$\\NewPoint$};" - elseif(step == "cavity") then - polygon = buildCavity(badTriangles, triangulation) - polyNew = cleanPoly(polygon) - -- remove the bad triangles - for j=1,#badTriangles do - table.remove(triangulation,badTriangles[j]-(j-1)) - end - -- draw the triangles - for i=1,#triangulation do - PointI = listPoints[triangulation[i][1]] - PointJ = listPoints[triangulation[i][2]] - PointK = listPoints[triangulation[i][3]] - output = output .. "\\draw[color="..color.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;" - end - -- fill and draw the cavity - path = "" - for i=1,#polyNew do - PointI = listPoints[polyNew[i]] - path = path .. "(".. PointI.x ..",".. PointI.y ..")--" - end - output = output .. "\\draw[color="..colorNew..",fill ="..colorBack..", thick] " .. path .. "cycle;" - -- mark the points of the mesh - for i=1,#listPoints do - output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};" - end - -- mark the adding point - output = output .. "\\draw[color="..colorNew.."] (" .. P.x ..",".. P.y .. ") node {$\\bullet$} node[anchor=north east] {$\\NewPoint$};" - elseif(step == "newT") then - polygon = buildCavity(badTriangles, triangulation) - polyNew = cleanPoly(polygon) - -- remove the bad triangles - for j=1,#badTriangles do - table.remove(triangulation,badTriangles[j]-(j-1)) - end - -- draw the triangle of the triangulation - for i=1,#triangulation do - PointI = listPoints[triangulation[i][1]] - PointJ = listPoints[triangulation[i][2]] - PointK = listPoints[triangulation[i][3]] - output = output .. "\\draw[color ="..color.."] (".. PointI.x ..",".. PointI.y ..")--("..PointJ.x..",".. PointJ.y ..")--("..PointK.x..",".. PointK.y ..")--cycle;" - end - -- fill and draw the cavity - path = "" - for i=1,#polyNew do - PointI = listPoints[polyNew[i]] - path = path .. "(".. PointI.x ..",".. PointI.y ..")--" - end - output = output .. "\\draw[color="..colorNew..",fill ="..colorBack..", thick] " .. path .. "cycle;" - -- draw the new triangles composed by the edges of the polygon and the added point - for i=1,#polygon do - 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 ..");" - output = output .. "\\draw[color="..colorNew..", thick]".."(".. listPoints[polygon[i][1]].x .. "," .. listPoints[polygon[i][1]].y .. ") -- (" .. P.x .. "," .. P.y ..");" - output = output .. "\\draw[color="..colorNew..", thick]".."(".. listPoints[polygon[i][2]].x .. "," .. listPoints[polygon[i][2]].y .. ") -- (" .. P.x .. "," .. P.y ..");" - end - -- mark points - for i=1,#listPoints do - output = output .. "\\draw[color="..color.."] (" .. listPoints[i].x ..",".. listPoints[i].y .. ") node {$\\bullet$} node[anchor=north east] {$\\MeshPoint_{" .. i .. "}$};" - end - -- mark the added point - output = output .. "\\draw[color="..colorNew.."] (" .. P.x ..",".. P.y .. ") node {$\\bullet$} node[anchor=north east] {$\\NewPoint$};" - end - return output -end - -function TeXaddOnePointMP(listPoints,P,step,color,colorBack, colorNew, colorCircle) - output = ""; - output = output .. "pair MeshPoints[];" - -- build the triangulation - triangulation = BowyerWatson(listPoints,"none") - badTriangles = buildBadTriangles(P,triangulation) - for i=1,#listPoints do - output = output .. "MeshPoints[".. i .. "] = (" .. listPoints[i].x .. "," .. listPoints[i].y .. ");" - end - if(step == "badT") then - -- draw all triangle - for i=1,#triangulation do - PointI = listPoints[triangulation[i][1]] - PointJ = listPoints[triangulation[i][2]] - PointK = listPoints[triangulation[i][3]] - output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\mpcolor{" .. color .."};" - end - -- draw and fill the bad triangle - for i=1,#badTriangles do - PointI = listPoints[triangulation[badTriangles[i]][1]] - PointJ = listPoints[triangulation[badTriangles[i]][2]] - PointK = listPoints[triangulation[badTriangles[i]][3]] - output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\mpcolor{" .. color .."};" - output = output .. "fill (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\mpcolor{" .. colorBack .."};" - end - -- draw the circoncircle - for i=1,#badTriangles do - PointI = listPoints[triangulation[badTriangles[i]][1]] - PointJ = listPoints[triangulation[badTriangles[i]][2]] - PointK = listPoints[triangulation[badTriangles[i]][3]] - center, radius = circoncircle(PointI, PointJ, PointK) - output = output .. "draw fullcircle scaled ("..radius .."*2u) shifted ("..center.x .. "*u," .. center.y .. "*u) dashed evenly withcolor \\mpcolor{" .. colorCircle .."};" - end - -- mark the points - for i=1,#listPoints do - output = output .. "dotlabel.llft (btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\mpcolor{" .. color .."};" - end - -- mark the point to add - output = output .. "dotlabel.llft (btex $\\NewPoint$ etex,(" .. P.x ..",".. P.y .. ")*u) withcolor \\mpcolor{" .. colorNew .."};" - elseif(step == "cavity") then - polygon = buildCavity(badTriangles, triangulation) - polyNew = cleanPoly(polygon) - -- remove the bad triangles - for j=1,#badTriangles do - table.remove(triangulation,badTriangles[j]-(j-1)) - end - -- draw the triangles - for i=1,#triangulation do - PointI = listPoints[triangulation[i][1]] - PointJ = listPoints[triangulation[i][2]] - PointK = listPoints[triangulation[i][3]] - output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\mpcolor{" .. color .."};" - end - -- fill and draw the cavity - path = "" - for i=1,#polyNew do - PointI = listPoints[polyNew[i]] - path = path .. "(".. PointI.x ..",".. PointI.y ..")*u--" - end - output = output .. "fill " .. path .. "cycle withcolor \\mpcolor{" .. colorBack .."};" - output = output .. "draw " .. path .. "cycle withcolor \\mpcolor{" .. colorNew .."} withpen pencircle scaled 1pt;" - -- mark the points of the mesh - for i=1,#listPoints do - output = output .. "dotlabel.llft (btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\mpcolor{" .. color .."};" - end - -- mark the adding point - output = output .. "dotlabel.llft (btex $\\NewPoint$ etex,(" .. P.x ..",".. P.y .. ")*u) withcolor \\mpcolor{" .. colorNew .."};" - elseif(step == "newT") then - polygon = buildCavity(badTriangles, triangulation) - polyNew = cleanPoly(polygon) - -- remove the bad triangles - for j=1,#badTriangles do - table.remove(triangulation,badTriangles[j]-(j-1)) - end - -- draw the triangle of the triangulation - for i=1,#triangulation do - PointI = listPoints[triangulation[i][1]] - PointJ = listPoints[triangulation[i][2]] - PointK = listPoints[triangulation[i][3]] - output = output .. "draw (".. PointI.x ..",".. PointI.y ..")*u--("..PointJ.x..",".. PointJ.y ..")*u--("..PointK.x..",".. PointK.y ..")*u--cycle withcolor \\mpcolor{" .. color .."};" - end - -- fill the cavity - path = "" - for i=1,#polyNew do - PointI = listPoints[polyNew[i]] - path = path .. "(".. PointI.x ..",".. PointI.y ..")*u--" - end - output = output .. "fill " .. path .. "cycle withcolor \\mpcolor{" .. colorBack .."};" - -- draw the new triangles composed by the edges of the polygon and the added point - for i=1,#polygon do - 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;" - 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;" - 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;" - end - -- mark points - for i=1,#listPoints do - output = output .. "dotlabel.llft (btex $\\MeshPoint_{" .. i .. "}$ etex, (" .. listPoints[i].x ..",".. listPoints[i].y .. ")*u ) withcolor \\mpcolor{" .. color .."};" - end - -- mark the added point - output = output .. "dotlabel.llft (btex $\\NewPoint$ etex,(" .. P.x ..",".. P.y .. ")*u) withcolor \\mpcolor{ " .. colorNew .."};" - end - return output -end -- build the list of points extern and stop at nbr function buildListExt(chaine, stop) - listPoints = {} + local listPoints = {} io.input(chaine) -- open the file text=io.read("*all") lines=string.explode(text,"\n+") -- all the lines @@ -569,23 +480,66 @@ function buildListExt(chaine, stop) return point, listPoints end - -function TeXFullOnePointTikZ(chaine,point,step,color,colorBack,colorNew,colorCircle,scale) - output = TeXaddOnePointTikZ(chaine,point,step,color,colorBack,colorNew,colorCircle) - output = "\\noindent\\begin{tikzpicture}[x="..scale..",y="..scale.."]".. output .. "\\end{tikzpicture}" - tex.sprint(output) +function split(pString, pPattern) + local Table = {} -- NOTE: use {n = 0} in Lua-5.0 + local fpat = "(.-)" .. pPattern + local last_end = 1 + local s, e, cap = pString:find(fpat, 1) + while s do + if s ~= 1 or cap ~= "" then + table.insert(Table,cap) + end + last_end = e+1 + s, e, cap = pString:find(fpat, last_end) + end + if last_end <= #pString then + cap = pString:sub(last_end) + table.insert(Table, cap) + end + return Table end -function TeXFullOnePointMP(chaine,point,step,color,colorBack,colorNew,colorCircle,scale,mode) - if(mode=="int") then - Sx,Sy=string.match(point,"%((.+),(.+)%)") - P = {x=Sx, y=Sy} - listPoints = buildList(chaine, mode) - else - -- point is a number - P, listPoints = buildListExt(chaine,tonumber(point)) +function readGmsh(file) + io.input(file) -- open the file + text=io.read("*all") + local lines = split(text,"\n+") -- all the lines + local listPoints={} + local triangulation ={} + boolNodes = false + Jnodes = 0 + boolElements = false + Jelements = 0 + J=0 + for i=1,#lines-J do + if(lines[i+J] == "$EndNodes") then + boolNodes = false + -- go to the next line + end + if(boolNodes) then -- we are in the Nodes environment + xy=split(lines[i+J]," +") + table.insert(listPoints,{x=tonumber(xy[2]),y=tonumber(xy[3])}) + end + if(lines[i+J] == "$Nodes") then + boolNodes = true + -- go to the next line + J=J+1 + end + if(lines[i+J] == "$EndElements") then + boolElements = false + -- go to the next line + end + if(boolElements) then -- we are in the Nodes environment + xy=split(lines[i+J]," +") + if(tonumber(xy[2]) == 2) then -- if the element is a triangle + nbrTags = xy[3]+1 + table.insert(triangulation,{tonumber(xy[2+nbrTags+1]),tonumber(xy[2+nbrTags+2]),tonumber(xy[2+nbrTags+3])}) + end + end + if(lines[i+J] == "$Elements") then + boolElements = true + -- go to the next line + J=J+1 + end end - output = TeXaddOnePointMP(listPoints,P,step,color,colorBack,colorNew,colorCircle) - output = "\\leavevmode\\begin{mplibcode}beginfig(0);u:="..scale..";".. output .. "endfig;\\end{mplibcode}" - tex.sprint(output) + return listPoints, triangulation end