diff options
Diffstat (limited to 'graphics/pgf/contrib/pgfplots/tex/prepcontour.lua')
-rwxr-xr-x | graphics/pgf/contrib/pgfplots/tex/prepcontour.lua | 708 |
1 files changed, 708 insertions, 0 deletions
diff --git a/graphics/pgf/contrib/pgfplots/tex/prepcontour.lua b/graphics/pgf/contrib/pgfplots/tex/prepcontour.lua new file mode 100755 index 0000000000..cfed282c24 --- /dev/null +++ b/graphics/pgf/contrib/pgfplots/tex/prepcontour.lua @@ -0,0 +1,708 @@ +#!/usr/bin/luatex + +--[[ + +prepcontour [Lua variant] - prepare contour lines (for pgfplots) + +Version: 1.4 (2021-02-22) + +Copyright (C) 2020-2021 Francesco Poli <invernomuto@paranoici.org> + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +--]] + +-- class creation function +function FpolyCreateClass() + local newclass = {} -- this table is the created class (which will act + -- as metatable for the instantiated objects) + + newclass.__index = newclass -- set the __index metamethod (so that + -- members missing from objects are searched + -- for in the class) + + newclass.new = function (...) -- member function that will be used to + -- instantiate a new object of this class + local object = setmetatable({}, newclass) -- a table with this class + -- set as metatable + + if object.init then + object:init(...) -- run the init method, if present + end + + return object + end + + return newclass +end + +-- class for the coordinates (and meta value) of a point in 3D space +Coord = FpolyCreateClass() + +function Coord:init() + self.x = { nil, nil, nil } -- point x,y,z coordinates + self.meta = nil -- point meta value +end + +-- class for a 2D mesh of points in 3D space +PrepcMesh = FpolyCreateClass() + +function PrepcMesh:init(yvaries, nblocks, nlines, copylines, + input, output) + if yvaries then + self.ni = nblocks - 1 -- number of cells in the i direction + self.nj = nlines - 1 -- number of cells in the j direction + else + self.ni = nlines - 1 -- number of cells in the i direction + self.nj = nblocks - 1 -- number of cells in the j direction + end + self.st = self.nj + 1 -- stride (number of nodes in the j direction) + self.si = self.nj + 1 -- stride for done_i + self.sj = self.nj -- stride for done_j + local cl = tonumber(copylines) or 0 + self.is = input or io.stdin + self.os = output or io.stdout + + self.debug = false + + local st = self.st -- shorter name + + self.coords = {} -- mesh of nodes + -- use single-index table for efficiency + -- (i,j) will be mapped to [i*st+j] + --[[ + The mesh of nodes represented by this table is: + + nj o o o o o o o o o o o o + + ^ o o o o o o o o o o o o + | + 2| o o o o o o o o o o o o + | + 1| o o o o o o o o o o o o + | + j=0+ o o o o o o o o o o o o + + +-----------------------------> + i=0 1 2 ni + --]] + + -- copy the first "cl" text lines from self.is to self.os + for l = 1, cl do + self.os:write(self.is:read('L')) + end + + -- read all the nodes + for bl = 0, nblocks - 1 do + if bl > 0 then + self.is:read('L') -- this text line should be empty + end + for ln = 0, nlines - 1 do + if yvaries then + i, j = bl, ln + else + i, j = ln, bl + end + n = Coord.new() + n.x[1], n.x[2], n.x[3], n.meta = self.is:read('n', 'n', 'n', 'n') + self.coords[i*st+j] = n + end + end +end + + +--[[ begin core of the program logic ]]-- + +-- build N contour lines between meta_min and meta_max +function PrepcMesh:autocontour(N, meta_min, meta_max, corners, tolerance) + -- subdivide the meta_min÷meta_max interval into N equal sub-intervals + -- and pick the midpoints of those sub-intervals + meta_min = tonumber(meta_min) + meta_max = tonumber(meta_max) + local step = (meta_max - meta_min)/N + local meta_mid = (meta_max + meta_min)/2 + local n_mid = (N + 1)/2 + for n = 1, N do + local isoval = meta_mid + (n - n_mid)*step + self:contour(isoval, corners, tolerance) + end +end + +-- build contour lines for meta==isoval +function PrepcMesh:contour(isoval, corners, tolerance) + -- set relative tolerance + local tol = tolerance or 1e-3 + + -- short names + local st = self.st + local si = self.si + local sj = self.sj + + -- local variables + local i, j + + if self.done_i == nil then + self.done_i = {} -- table of markers for i-sides + -- again single-index table + -- (i,j) will be mapped to [i*si+j] + --[[ + The corresponding i-sides are: + + nj o--o--o--o--o--o--o--o--o--o--o-- + + ^ o--o--o--o--o--o--o--o--o--o--o-- + | + 2| o--o--o--o--o--o--o--o--o--o--o-- + | + 1| o--o--o--o--o--o--o--o--o--o--o-- + | + j=0+ o--o--o--o--o--o--o--o--o--o--o-- + + +--------------------------> + i=0 1 2 ni-1 + --]] + end + + if self.done_j == nil then + self.done_j = {} -- table of markers for j-sides + -- again single-index table + -- (i,j) will be mapped to [i*sj+j] + --[[ + The corresponding j-sides are: + + | | | | | | | | | | | | + nj-1 o o o o o o o o o o o o + | | | | | | | | | | | | + ^ o o o o o o o o o o o o + | | | | | | | | | | | | | + 1| o o o o o o o o o o o o + | | | | | | | | | | | | | + j=0+ o o o o o o o o o o o o + + +-----------------------------> + i=0 1 2 ni + --]] + end + + if self.newcoords == nil then + self.newcoords = {} -- nodes for contour lines will be placed here + end + + -- scan all the i-sides, searching for + -- intersections with the level surface + for i = 0, self.ni - 1 do + for j = 0, self.nj do + self.done_i[i*si+j] = ((self.coords[ i*st+j].meta > isoval) == + (self.coords[(i+1)*st+j].meta > isoval)) + --[[ + This is equivalent to testing the two nodes of the i-side: + + j o-----o + i i+1 + + If meta > isoval in both nodes, or in none of them, then + set done_i to true (never look again at this i-side). + Otherwise, set done_i to false (we have not finished with + this i-side, since it contains one point of a contour line). + --]] + end + end + + -- scan all the j-sides, searching for + -- intersections with the level surface + for i = 0, self.ni do + for j = 0, self.nj - 1 do + self.done_j[i*sj+j] = ((self.coords[i*st+j ].meta > isoval) == + (self.coords[i*st+j+1].meta > isoval)) + --[[ + This is equivalent to testing the two nodes of the j-side: + + j+1 o + | + | + j o + i + + Similarly, set done_j accordingly... + --]] + end + end + + if self.debug then + self:show_sides(isoval) + end + + --[[ + We now have all the sides containing points of the meta==isoval + contour lines. We need to connect the points into distinct lines. + + First of all, scan all the boundary sides, searching for sides + with done==false. For each boundary side with done==false, set + its done to true and build one contour line, starting from + the side itself and the real cell adjacent to it. + + How can we specify a cell and one of its sides? + We can specify cell (i,j) and a side (South, West, North, East) + + N + j o-------o + | | + W | (i,j) | E + | | + j-1 o-------o + i-1 S i + --]] + + j = 0 -- i-sides with j=0 + for i = 0, self.ni - 1 do + if not self.done_i[i*si+j] then + self.done_i[i*si+j] = true + self:build_line(isoval, tol, i+1, j+1, 'S', corners) + if self.debug then + self:show_sides(isoval) + end + end + end + + i = 0 -- j-sides with i=0 + for j = 0, self.nj - 1 do + if not self.done_j[i*sj+j] then + self.done_j[i*sj+j] = true + self:build_line(isoval, tol, i+1, j+1, 'W', corners) + if self.debug then + self:show_sides(isoval) + end + end + end + + j = self.nj -- i-sides with j=nj + for i = 0, self.ni - 1 do + if not self.done_i[i*si+j] then + self.done_i[i*si+j] = true + self:build_line(isoval, tol, i+1, j , 'N', corners) + if self.debug then + self:show_sides(isoval) + end + end + end + + i = self.ni -- j-sides with i=ni + for j = 0, self.nj - 1 do + if not self.done_j[i*sj+j] then + self.done_j[i*sj+j] = true + self:build_line(isoval, tol, i , j+1, 'E', corners) + if self.debug then + self:show_sides(isoval) + end + end + end + + --[[ + Finally, scan all the internal sides, again searching for sides + with done==false. For each internal side with done==false, do _not_ + alter its done value and build one contour line, starting from + the side itself and one of the two real cells adjacent to it. + --]] + + -- i-sides with j=1,...,nj-1 (i-sides with j=0 are already done) + -- j-sides with i=1,...,ni-1 (j-sides with i=0 are already done) + for j = 0, self.nj - 1 do + for i = 0, self.ni - 1 do + if not self.done_i[i*si+j] then + self:build_line(isoval, tol, i+1, j+1, 'S', corners) + if self.debug then + self:show_sides(isoval) + end + end + if not self.done_j[i*sj+j] then + self:build_line(isoval, tol, i , j+1, 'E', corners) + if self.debug then + self:show_sides(isoval) + end + end + end + end +end + +-- build a single contour line +function PrepcMesh:build_line(isoval, tol, ic, jc, side, corners) + -- short names + local st = self.st + local si = self.si + local sj = self.sj + + -- local variables + local ia, ja, na, wa + local ib, jb, nb, wb + local ie, je, ne + local id, jd, nd + local pt, k, count, next_side + local cc, xi_k, eta_k, ck, above, sign_ck + local xi_v, eta_v, xi_0, xi_1, eta_0, eta_1 + local next_done, next_done_idx + + while true do + --[[ + Start from cell (ic,jc) and its South, West, North, or East side + + e d b d a N b d b + o-------o o-------o o-*-----o o-------o + | | | | | | | * + | c | W | c | | c | | c | E + | | * | | | | | + o----*--o o-------o o-------o o-------o + a S b a e e d e a + --]] + + -- find nodes a and b + if side == 'E' then ia = ic else ia = ic - 1 end + if side == 'N' then ja = jc else ja = jc - 1 end + if side == 'W' then ib = ic - 1 else ib = ic end + if side == 'S' then jb = jc - 1 else jb = jc end + if ia < 0 or ib > self.ni or ja < 0 or jb > self.nj then + self.newcoords[#self.newcoords + 1] = Coord.new() + break -- out of range: abort + end + na = self.coords[ia*st+ja] -- node a + nb = self.coords[ib*st+jb] -- node b + + -- compute the intersection point between the side + -- and the meta==isoval level surface + wa = (isoval - nb.meta)/(na.meta - nb.meta) -- weights + wb = 1 - wa + + pt = Coord.new() + for k = 1, 3 do + pt.x[k] = wa*na.x[k] + wb*nb.x[k] -- intersection point + end + pt.meta = isoval + -- add the intersection point to the contour line + self.newcoords[#self.newcoords + 1] = pt + + -- check the cell: if it is a phantom cell, then stop + if ic < 1 or ic > self.ni or jc < 1 or jc > self.nj then + self.newcoords[#self.newcoords + 1] = Coord.new() + break -- the contour line is finished + end + + -- find nodes e and d + if side == 'W' then ie = ic else ie = ic - 1 end + if side == 'S' then je = jc else je = jc - 1 end + if side == 'E' then id = ic - 1 else id = ic end + if side == 'N' then jd = jc - 1 else jd = jc end + ne = self.coords[ie*st+je] -- node e + nd = self.coords[id*st+jd] -- node d + + -- look at the other three sides of the cell: + -- how many of them have done==false ? + count = 0 + if side ~= 'S' and not self.done_i[(ic-1)*si+(jc-1)] then + count = count + 1 + next_side = 'N' -- next side could be the North side of the + -- adjacent south cell + end + if side ~= 'W' and not self.done_j[(ic-1)*sj+(jc-1)] then + count = count + 1 + next_side = 'E' -- next side could be the East side of the + -- adjacent west cell + end + if side ~= 'N' and not self.done_i[(ic-1)*si+ jc ] then + count = count + 1 + next_side = 'S' -- next side could be the South side of the + -- adjacent north cell + end + if side ~= 'E' and not self.done_j[ ic *sj+(jc-1)] then + count = count + 1 + next_side = 'W' -- next side could be the West side of the + -- adjacent east cell + end + + --[[ + if zero other sides have done==false, then stop + + o o + + + o---*-o + --]] + if count == 0 then + self.newcoords[#self.newcoords + 1] = Coord.new() + break -- the contour line is finished + end + + --[[ + if exactly one other side has done==false, then + next_side already specifies the next side to start from + + o-*---o + | + \ + o---*-o + + nothing to be done for the time being... + --]] + + --[[ + the surface inside the cell is thought as the bilinear + interpolation of the four nodes and hence mapped onto + a unit square in the xi,eta plane: + + eta ^ + | pt.x[k] = (1-xi)*(1-eta)*na.x[k] + + 1+------+ xi *(1-eta)*nb.x[k] + + |e d| (1-xi)* eta *ne.x[k] + + | | xi * eta *nd.x[k] + |a b| + 0+------+------> + 0 1 xi + + the equation of the contour line is therefore: + cc*xi*eta - eta_k*xi - xi_k*eta == ck + --]] + cc = na.meta - nb.meta + nd.meta - ne.meta + eta_k = na.meta - nb.meta + xi_k = na.meta - ne.meta + ck = isoval - na.meta + + if math.abs(cc) <= + math.max(math.abs(tol*eta_k), math.abs(tol*xi_k)) then + --[[ + cc is negligible and the contour line is + basically a straight line: + eta_k*xi + xi_k*eta + ck == 0 + --]] + cc = 0 + else + --[[ + cc is non-zero (and non-negligible) and we can divide + both sides of the equation by cc: + xi*eta - (eta_k/cc)*xi - (xi_k/cc)*eta == ck/cc + --]] + xi_k = xi_k /cc -- let's rename the coefficients + eta_k = eta_k/cc + ck = xi_k*eta_k + ck/cc + --[[ + the contour line is an equilateral hyperbola (in the + xi,eta plane) with asymptotes xi == xi_k and eta == eta_k + and constant product ck: + (xi - xi_k)*(eta - eta_k) == ck + --]] + if math.abs(ck) <= math.abs(tol/2) then + --[[ + ck is negligible and the hyperbola basically + degenerates into its asymptotes + --]] + ck = 0 + end + end + + --[[ + if two or three other sides have done==false, then + choose the next side + + ? + o-*---o + | * ? + ? * | + o---*-o + + --]] + if count >= 2 then + -- how can we choose the next side? + --[[ + if there's more than one side to choose from, cc is + necessarily non-zero (a straight line could not cross + more than two sides of the unit square) and the + center k of the hyperbola is necessarily inside the + cell (otherwise one hyperbola branch would be completely + outside the cell and the other branch could not cross + more than two sides) + + hence, choose the next side by checking the quadrant + where the hyperbola branch crossing side a-b lies + + eta ^ + | : + 1+--:---+ + |e : d| + ······+k······ + |a : b| + 0+--:---+------> + 0 : 1 xi + --]] + if ck == 0 then -- math.abs(ck) <= math.abs(tol/2) + -- choose side e-d + next_side = side + elseif ck < 0 then -- ck < -tol/2 + -- choose side b-d + if side == 'N' or side == 'S' then + next_side = 'W' + else + next_side = 'S' + end + else -- ck > +tol/2 + -- choose side a-e + if side == 'N' or side == 'S' then + next_side = 'E' + else + next_side = 'N' + end + end + end + + if corners then -- enhanced corner algorithm + + --[[ + if, within the cell, the contour line is not straight + we can improve its representation by computing an + additional point: the vertex v of the hyperbola branch, + as long as it lies inside the cell + --]] + if cc ~= 0 then + -- depending on where the center k is, consider the + -- vertex above or below k + if eta_k < 0 then above = 1 else above = -1 end + -- also check the sign of ck + if ck < 0 then sign_ck = -1 else sign_ck = 1 end + -- compute the vertex v + xi_v = xi_k + sign_ck*above*math.sqrt(math.abs(ck)) + eta_v = eta_k + above*math.sqrt(math.abs(ck)) + + --[[ + the vertex will be considered inside the cell, as long as + xi_0 <= xi_v <= xi_1 and eta_0 <= eta_v <= eta_1 + --]] + xi_0 = - math.abs(tol/10) + xi_1 = 1 + math.abs(tol/10) + eta_0 = math.abs(tol) + eta_1 = 1 + math.abs(tol/10) + -- eta_0 is meant to reject a vertex too close to side a-b + -- (which already has a contour line point!) + -- we should also reject a vertex too close to the next side + if next_side == side then + -- next side is e-d + eta_1 = 1 - math.abs(tol) + elseif next_side == 'W' or next_side == 'S' then + -- next side is b-d + xi_1 = 1 - math.abs(tol) + else -- next_side == 'E' or next_side == 'N' + -- next side is a-e + xi_0 = math.abs(tol) + end + + if xi_0 <= xi_v and xi_v <= xi_1 and + eta_0 <= eta_v and eta_v <= eta_1 then + -- v is inside the cell: compute its actual coordinates + pt = Coord.new() + for k = 1, 3 do + pt.x[k] = (1-xi_v)*(1-eta_v)*na.x[k] + + xi_v *(1-eta_v)*nb.x[k] + + (1-xi_v)* eta_v *ne.x[k] + + xi_v * eta_v *nd.x[k] + end + pt.meta = isoval + -- add the vertex to the contour line + self.newcoords[#self.newcoords + 1] = pt + end + end + + end + + -- point to done value for the next side and + -- determine indexes for the next cell + -- (the other cell adjacent to the next side) + if next_side == 'W' then ic = ic + 1 end + if next_side == 'S' then jc = jc + 1 end + if next_side == 'N' or next_side == 'S' then + next_done = self.done_i + next_done_idx = (ic-1)*si+(jc-1) + else + next_done = self.done_j + next_done_idx = (ic-1)*sj+(jc-1) + end + if next_side == 'E' then ic = ic - 1 end + if next_side == 'N' then jc = jc - 1 end + + -- if next side is already done, then stop + if next_done[next_done_idx] then + self.newcoords[#self.newcoords + 1] = Coord.new() + break -- the contour line is finished + end + + -- set next side done to true and + -- iterate starting from this new side and the next cell + next_done[next_done_idx] = true + side = next_side + end +end + +--[[ end core of the program logic ]]-- + + +-- print out all the contour lines +function PrepcMesh:printcontours() + for n = 1, #self.newcoords do + if self.newcoords[n].x[1] then + self.os:write(string.format("%14.6g%14.6g%14.6g%14.6g\n", + self.newcoords[n].x[1], + self.newcoords[n].x[2], + self.newcoords[n].x[3], + self.newcoords[n].meta)) + else + self.os:write('\n') + end + end +end + +-- print a debug representation of mesh and not-done sides +function PrepcMesh:show_sides(isoval) + -- short names + local st = self.st + local si = self.si + local sj = self.sj + + self.os:write(string.format("# contour level = %.6g \n#\n", isoval)) + for j = self.nj, 0, -1 do + if j <= self.nj - 1 then + self.os:write("# ") + for i = 0, self.ni do + if self.done_j[i*sj+j] then + self.os:write(" ") + else + self.os:write(" *") + end + end + self.os:write("\n") + end + self.os:write("# ") + for i = 0, self.ni do + if self.coords[i*st+j].meta > isoval then + self.os:write("•") + else + self.os:write("o") + end + if i <= self.ni - 1 then + if self.done_i[i*si+j] then + self.os:write(" ") + else + self.os:write("*") + end + end + end + self.os:write("\n") + end + self.os:write("#\n") +end + |