A logistical way to recursively subdivide triangles (into 6 sub-triangles), such that all detail is captured and no triangles are wasted

Question

I have a surface composed of triangles. I want to recursively subdivide these triangles into 6 sub-triangles, which are then themselves further subdivided into 6 more and so on.

I want for this subdivision to carry on only until enough detail is captured, and no more after that.

I do not know what the heck I should measure for this. I tried using the laplacian, but it didn't seem to work. Maybe there is a better measurement, or perhaps either my laplaican function could be wrong.

There are two functions called append_surface in the MWE. swap their order to see each's effect. In particular, I've included the current function, along with my experimental function which I made in an attempt to solve this question.

The current function produces:

The experimental function produces:

mwe.tex

\documentclass[ tikz ,border = 1cm ]{standalone} \makeatletter \tikzset{ /lua-tikz3dtools/.is family % abbreviated "td" ,/lua-tikz3dtools/.cd ,/lua-tikz3dtools/parametric/.cd % appreviated "p" ,/lua-tikz3dtools/parametric/matrix/.cd % abbreviated "m" ,object/.code = {\protected@edef\luatikztdtools@p@m@transformation{#1}} ,object = {identity_matrix()} ,name/.code = {\protected@edef\luatikztdtools@p@m@name{#1}} ,name = {Change me!} ,/lua-tikz3dtools/parametric/surface/.cd % abbreviated "s" ,ustart/.estore in = \luatikztdtools@p@s@ustart ,ustop/.estore in = \luatikztdtools@p@s@ustop ,usamples/.estore in = \luatikztdtools@p@s@usamples ,vstart/.estore in = \luatikztdtools@p@s@vstart ,vstop/.estore in = \luatikztdtools@p@s@vstop ,vsamples/.estore in = \luatikztdtools@p@s@vsamples ,x/.code = {\protected@edef\luatikztdtools@p@s@x{#1}} ,y/.code = {\protected@edef\luatikztdtools@p@s@y{#1}} ,z/.code = {\protected@edef\luatikztdtools@p@s@z{#1}} ,transformation/.code = {\protected@edef\luatikztdtools@p@s@transformation{#1}} ,fill options/.code = {\protected@edef\luatikztdtools@p@s@filloptions{#1}} ,ustart = {0} ,ustop = {1} ,usamples = {10} ,vstart = {0} ,vstop = {1} ,vsamples = {10} ,x = {sphere(tau*u, pi*v)[1][1]} ,y = {sphere(tau*u, pi*v)[1][2]} ,z = {sphere(tau*u, pi*v)[1][3]} ,transformation = {identity_matrix()} ,fill options = {fill, fill opacity = 0.3} } \ExplSyntaxOn \lua_load_module:n { mwe-lua } \NewDocumentCommand { \appendsurface } {o} { \group_begin: \tikzset{ /lua-tikz3dtools/parametric/surface/.search~also = {/tikz} ,/lua-tikz3dtools/parametric/surface/.cd ,#1 } \__lua_tikztdtools_appendsurface: \group_end: } \NewDocumentCommand { \displaysegments } {} { \group_begin: \__lua_tikztdtools_displaysegments: \group_end: } \NewDocumentCommand { \setobject } {o} { \group_begin: \tikzset{ /lua-tikz3dtools/parametric/matrix/.search~also = {/tikz} ,/lua-tikz3dtools/parametric/matrix/.cd ,#1 } \__lua_tikztdtools_setobject: \group_end: } \ExplSyntaxOff \makeatother \begin{document} \begin{tikzpicture} \setobject[ name = {T} ,object = { matrix_multiply( euler(pi/2,pi/3,5*pi/6) ,translate(0,0,-15) ) } ] \appendsurface[ ustart = {-7} ,ustop = {7} ,usamples = {7} ,vstart = {-7} ,vstop = {7} ,vsamples = {7} ,transformation = {T} ,x = {u} ,y = {v} ,z = { 3 / ((u)^2+(v)^2+1/2) } ,fill options = { preaction = { fill = red!70!black ,fill opacity = 1 } ,postaction = { draw = black ,ultra thin ,line join = round ,line cap = round } } ] \displaysegments \end{tikzpicture} \end{document}

mwe-lua.lua

 local eps = 0.0000001 --- scalar multiplication --- @param a number the scalar coefficient --- @param v table<table<number>> the vector --- @return table<table<number>> the scaled vector local function scalar_multiplication(a, v) local tmp = v[1] return { {a*tmp[1], a*tmp[2], a*tmp[3], 1} } end --- vector subtraction --- @param u table<table<number>> the first vector --- @param v table<table<number>> the second vector --- @return table<table<number>> their difference local function vector_subtraction(u, v) local U = u[1] local V = v[1] return { {U[1]-V[1], U[2]-V[2], U[3]-V[3], 1} } end --- vector addition --- @param u table<table<number>> the first vector --- @param v table<table<number>> the second vector --- @return table<table<number>> their sum local function vector_addition(u, v) local U = u[1] local V = v[1] return { {U[1]+V[1], U[2]+V[2], U[3]+V[3], 1} } end --- length of a vector --- @param u table<table<number>> the vector to be measured --- @return number the vectors length local function length(u) local result = math.sqrt((u[1][1])^2 + (u[1][2])^2 + (u[1][3])^2) return result end --- distance between two points --- @param P1 table<table<number>> a point --- @param P2 table<table<number>> another point --- @return number their distance local function distance(P1, P2) local A = P1[1] local B = P2[1] return math.sqrt( (A[1]-B[1])^2 + (A[2]-B[2])^2 + (A[3]-B[3])^2 ) end -- https://tex.stackexchange.com/a/747040 --- Creates a TeX command that evaluates a Lua function --- --- @param name string The name of the `\csname` to define --- @param func function --- @param args table<string> The TeX types of the function arguments --- @param protected boolean|nil Define the command as `\protected` --- @return nil local function register_tex_cmd(name, func, args, protected) -- The extended version of this function uses `N` and `w` where appropriate, -- but only using `n` is good enough for exposition purposes. name = "__lua_tikztdtools_" .. name .. ":" .. ("n"):rep(#args) -- Push the appropriate scanner functions onto the scanning stack. local scanners = {} for _, arg in ipairs(args) do scanners[#scanners+1] = token['scan_' .. arg] end -- An intermediate function that properly "scans" for its arguments -- in the TeX side. local scanning_func = function() local values = {} for _, scanner in ipairs(scanners) do values[#values+1] = scanner() end func(table.unpack(values)) end local index = luatexbase.new_luafunction(name) lua.get_functions_table()[index] = scanning_func if protected then token.set_lua(name, index, "protected") else token.set_lua(name, index) end end local segments = {} local mm = {} mm.tau = 2*math.pi local pi = math.pi local cos, sin = math.cos, math.sin --- matrix multiplication --- --- @param A table<table<number>> left matrix --- @param B table<table<number>> right matrix --- @return table<table<number>> the product function mm.matrix_multiply(A, B) local rows_A = #A local columns_A = #A[1] local rows_B = #B local columns_B = #B[1] assert( columns_A == rows_B, string.format( [[ Wrong size matrices for multiplication. Size A: %d,%d Size B: %d,%d ]], rows_A, columns_A, rows_B, columns_B ) ) local product = {} for row = 1, rows_A do product[row] = {} for column = 1, columns_B do product[row][column] = 0 for dot_product_step = 1, columns_A do local a = A[row][dot_product_step] local b = B[dot_product_step][column] assert(type(a) == "number", string.format("Expected number but got %s in A[%d][%d]", type(a), row, dot_product_step)) assert(type(b) == "number", string.format("Expected number but got %s in B[%d][%d]", type(b), dot_product_step, column)) product[row][column] = product[row][column] + a * b end end end return product end function mm.yrotation(angle) local c = cos(angle) local s = sin(angle) return { {c,0,-s,0} ,{0,1,0,0} ,{s,0,c,0} ,{0,0,0,1} } end function mm.translate(x,y,z) return { {1,0,0,0} ,{0,1,0,0} ,{0,0,1,0} ,{x,y,z,1} } end function mm.zrotation(angle) local c = cos(angle) local s = sin(angle) return { {c,s,0,0} ,{-s,c,0,0} ,{0,0,1,0} ,{0,0,0,1} } end function mm.euler(alpha,beta,gamma) return mm.matrix_multiply( mm.zrotation(gamma) ,mm.matrix_multiply( mm.yrotation(beta) ,mm.zrotation(alpha) ) ) end function mm.midpoint(triangle) local P,Q,R = table.unpack(triangle) local x = (P[1]+Q[1]+R[1])/3 local y = (P[2]+Q[2]+R[2])/3 local z = (P[3]+Q[3]+R[3])/3 return {{x,y,z,1}} end local lua_tikz3dtools_parametric = {} lua_tikz3dtools_parametric.math = {} for k, v in pairs(_G) do lua_tikz3dtools_parametric.math[k] = v end for k, v in pairs(mm) do lua_tikz3dtools_parametric.math[k] = v end for k, v in pairs(math) do lua_tikz3dtools_parametric.math[k] = v end local function single_string_expression(str) return load(("return %s"):format(str), "expression", "t", lua_tikz3dtools_parametric.math)() end local function double_string_function(str) return load(("return function(u,v) return %s end"):format(str), "expression", "t", lua_tikz3dtools_parametric.math)() end local function append_surface(hash) local ustart = hash.ustart local ustop = hash.ustop local usamples = hash.usamples local vstart = hash.vstart local vstop = hash.vstop local vsamples = hash.vsamples local x = hash.x local y = hash.y local z = hash.z local drawoptions = hash.drawoptions local filloptions = hash.filloptions local transformation = hash.transformation x = double_string_function(x) y = double_string_function(y) z = double_string_function(z) ustart = single_string_expression(ustart) ustop = single_string_expression(ustop) usamples = single_string_expression(usamples) vstart = single_string_expression(vstart) vstop = single_string_expression(vstop) vsamples = single_string_expression(vsamples) transformation = single_string_expression(transformation) local ustep = (ustop - ustart) / (usamples - 1) local vstep = (vstop - vstart) / (vsamples - 1) local function parametric_surface(u, v) return { x(u,v), y(u,v), z(u,v), 1 } end local function distance(P1, P2) return math.sqrt( (P1[1]-P2[1])^2 + (P1[2]-P2[2])^2 + (P1[3]-P2[3])^2 ) end local eps = 0.0000001 for i = 0, usamples - 2 do local u = ustart + i * ustep for j = 0, vsamples - 2 do local v = vstart + j * vstep local A = parametric_surface(u, v) local B = parametric_surface(u + ustep, v) local C = parametric_surface(u, v + vstep) local D = parametric_surface(u + ustep, v + vstep) local the_segment1 = mm.matrix_multiply({ A, B, D },transformation) local the_segment2 = mm.matrix_multiply({ A, C, D },transformation) if not ( distance(A, B) < eps or distance(B, D) < eps or distance(A, D) < eps ) then table.insert( segments, { segment = the_segment1, filloptions = filloptions, type = "triangle" } ) end if not ( distance(A, C) < eps or distance(C, D) < eps or distance(A, D) < eps ) then table.insert( segments, { segment = the_segment2, filloptions = filloptions, type = "triangle" } ) end end end end local function append_surface(hash) local ustart = hash.ustart local ustop = hash.ustop local vstart = hash.vstart local vstop = hash.vstop local x = hash.x local y = hash.y local z = hash.z local filloptions = hash.filloptions local transformation = hash.transformation x = double_string_function(x) y = double_string_function(y) z = double_string_function(z) ustart = single_string_expression(ustart) ustop = single_string_expression(ustop) vstart = single_string_expression(vstart) vstop = single_string_expression(vstop) transformation = single_string_expression(transformation) --- parametric definition of our surface This function maps --- a point from the 2D domain onto the 3D codomain. --- @param u number the first coordinate component --- @param v number the second coordinate component --- @return table<number> the map of the given coordinates onto R3 local function parametric_surface(u, v) return { { x(u,v), y(u,v), z(u,v), 1 } } end --- distance between two DOMAIN points --- @param A table<table<number>> the first point in the DOMAIN --- @param B table<table<number>> the second point in the DOMAIN --- @return number their Euclidean distance in the DOMAIN local function dist2(A,B) local Ax, Ay = A[1], A[2] local Bx, By = B[1], B[2] return math.sqrt((Ax-Bx)^2+(Ay-By)^2) end --- gives the incentre of a triangle on the 2D domain --- @param T table<table<number>> --- @return table<number> the incentre (u,v) local function incentre2(T) local A,B,C = T[1],T[2],T[3] local DAB = dist2(A,B) local DBC = dist2(B,C) local DCA = dist2(C,A) local S = DAB + DBC + DCA if S == 0 then -- degenerate triangle, return centroid as fallback return { (A[1]+B[1]+C[1]) / 3, (A[2]+B[2]+C[2]) / 3 } end -- weighted average: (a*A + b*B + c*C) / (a+b+c) local dotx = (A[1]*DAB + B[1]*DBC + C[1]*DCA) / S local doty = (A[2]*DAB + B[2]*DBC + C[2]*DCA) / S return { dotx, doty } end --- partial derivatives of f(x,y) at the point (x=a,y=b) local function first_partial_derivatives(a,b) local eps = 1e-6 local dfdu = scalar_multiplication( 1/eps ,vector_subtraction( parametric_surface(a+eps, b) ,parametric_surface(a,b) ) ) local dfdv = scalar_multiplication( 1/eps ,vector_subtraction( parametric_surface(a, b+eps) ,parametric_surface(a,b) ) ) return {dfdu, dfdv} end --- second partial derivatives --- @param a number the first coordinate component --- @param b number the second coordinate component --- @return table<nnumber> the second partial derivates local function second_partial_derivatives(a, b) local eps = 1e-6 local FPD1, FPD2 = first_partial_derivatives(a, b)[1], first_partial_derivatives(a, b)[2] local FPD1h = first_partial_derivatives(a+eps, b)[1] local FPD2h = first_partial_derivatives(a, b+eps)[2] local d2fdu2 = scalar_multiplication( 1/eps ,vector_subtraction( FPD1h ,FPD1 ) ) local d2fdv2 = scalar_multiplication( 1/eps ,vector_subtraction( FPD2h ,FPD2 ) ) return {d2fdu2, d2fdv2} end local function laplacian(a,b) local P, Q = second_partial_derivatives(a, b)[1], second_partial_derivatives(a, b)[2] return scalar_multiplication(1/2,vector_addition(P, Q)) end -- local LL = parametric_surface(ustart, vstart) -- local LR = parametric_surface(ustop, vstart) -- local UL = parametric_surface(ustart, vstop) -- local UR = parametric_surface(ustop, vstop) local TRI1 = {{ustart, vstart}, {ustop, vstart}, {ustart, vstop}} local TRI2 = {{ustop, vstop}, {ustop, vstart}, {ustart, vstop}} local tris = {TRI1, TRI2} local max_depth = 4 -- you can tune this local iteration = 0 local min_DT = 0.05 local max_DT = 0.5 local tiny = 1e-12 while iteration < max_depth do local newtris = {} local refined = false for _, tri in ipairs(tris) do local v1, v2, v3 = tri[1], tri[2], tri[3] local P = incentre2(tri) -- Adaptive threshold based on Laplacian local PL = laplacian(P[1], P[2]) local px, py, pz = PL[1][1], PL[1][2], PL[1][3] local Lmag = math.sqrt(px*px + py*py + pz*pz) local invL = 1 / (Lmag + tiny) local DT = math.max(min_DT, math.min(invL, max_DT)) -- Distances to control refinement local D1 = dist2(P, v1) local D2 = dist2(P, v2) local D3 = dist2(P, v3) if (D1 < DT) or (D2 < DT) or (D3 < DT) then table.insert(newtris, tri) else -- Midpoints of edges (not P→vertex midpoints) local e12 = { (v1[1] + v2[1]) * 0.5, (v1[2] + v2[2]) * 0.5 } local e23 = { (v2[1] + v3[1]) * 0.5, (v2[2] + v3[2]) * 0.5 } local e31 = { (v3[1] + v1[1]) * 0.5, (v3[2] + v1[2]) * 0.5 } -- Six subtriangles around P that fill the parent: local t1 = { P, v1, e12 } local t2 = { P, e12, v2 } local t3 = { P, v2, e23 } local t4 = { P, e23, v3 } local t5 = { P, v3, e31 } local t6 = { P, e31, v1 } table.insert(newtris, t1) table.insert(newtris, t2) table.insert(newtris, t3) table.insert(newtris, t4) table.insert(newtris, t5) table.insert(newtris, t6) refined = true end end tris = newtris iteration = iteration + 1 if not refined then break end end for _, tri in ipairs(tris) do -- map UV to 3D local A = parametric_surface(tri[1][1], tri[1][2])[1] local B = parametric_surface(tri[2][1], tri[2][2])[1] local C = parametric_surface(tri[3][1], tri[3][2])[1] -- matrix-multiply by transformation local the_segment = mm.matrix_multiply({A, B, C}, transformation) if not ( distance({A}, {B}) < eps or distance({B}, {C}) < eps or distance({A}, {C}) < eps ) then table.insert( segments, { segment = the_segment, filloptions = filloptions, type = "triangle" } ) end end end register_tex_cmd( "appendsurface", function() append_surface{ ustart = token.get_macro("luatikztdtools@p@s@ustart"), ustop = token.get_macro("luatikztdtools@p@s@ustop"), usamples = token.get_macro("luatikztdtools@p@s@usamples"), vstart = token.get_macro("luatikztdtools@p@s@vstart"), vstop = token.get_macro("luatikztdtools@p@s@vstop"), vsamples = token.get_macro("luatikztdtools@p@s@vsamples"), x = token.get_macro("luatikztdtools@p@s@x"), y = token.get_macro("luatikztdtools@p@s@y"), z = token.get_macro("luatikztdtools@p@s@z"), transformation = token.get_macro("luatikztdtools@p@s@transformation"), filloptions = token.get_macro("luatikztdtools@p@s@filloptions"), } end, { } ) local function make_object(name, tbl) local obj = {} obj[name] = tbl return obj end local function set_matrix(hash) local transformation = single_string_expression(hash.transformation) local name = hash.name local tbl = make_object(name, transformation) for k, v in pairs(tbl) do lua_tikz3dtools_parametric.math[k] = v end end register_tex_cmd( "setobject", function() set_matrix{ name = token.get_macro("luatikztdtools@p@m@name"), transformation = token.get_macro("luatikztdtools@p@m@transformation"), } end, { } ) local function display_segments() local function comparator(a,b) return mm.midpoint(a.segment)[1][3] < mm.midpoint(b.segment)[1][3] end -- Perform sorting after filtering out skipped segments table.sort(segments, comparator) -- Iterate through the sorted segments for rendering for _, segment in ipairs(segments) do -- Handle triangle segments if segment.type == "triangle" then local P, Q, R = segment.segment[1], segment.segment[2], segment.segment[3] tex.sprint(string.format( "\\path[%s] (%f,%f) -- (%f,%f) -- (%f,%f) -- cycle;", segment.filloptions, P[1], P[2], Q[1], Q[2], R[1], R[2] )) end end -- Clear the segments array after processing segments = {} end register_tex_cmd("displaysegments", function() display_segments() end, { })

The key is the funciton to which this is a 3D surface plot. If you have that function, use it. Otherwise, you need an approximation, preferably one which will produce a sharp peak in the center (alas, not a polynomial). — John Kormylo
– John Kormylo, Commented Oct 20 at 15:18

kabenyuk · Accepted Answer · 2025-10-25 11:52:59Z

I apologize for posting a new answer instead of editing my previous one.

Nevertheless, I believe this presents a substantially different approach to adaptive mesh generation based on Delaunay triangulation. The main goal is to construct a high‑quality triangular mesh that is refined near user‑specified feature points and gradually coarsens with distance from them.

The algorithm works as follows:

Boundary corners, a central point, and a set of “refinement” points are placed inside a rectangular domain. Around each refinement point, an initial ring of auxiliary points is generated.
On each iteration, a Delaunay triangulation of the current point set is built (using the Bowyer–Watson algorithm with a super‑triangle). For every triangle, the circumcenter and circumradius are computed. If the radius exceeds a locally allowed cell size—determined by a smooth function of the distance to the nearest refinement point—the circumcenter is inserted as a new mesh point, provided it lies inside the domain and is not too close to existing points.
The process stops either when the maximum number of iterations is reached or when no new points are added (i.e., the mesh has converged).
The resulting triangles together with the function values z(x,y) are written to surface.dat, which is then visualized as a 3D surface using pgfplots via \addplot3[patch, patch type=triangle].
The figure shows the output for the test function
z(u,v) = 3/((u+2)^2 + (v+2)^2 + 0.5) + 4/((u-2)^2 + (v-2)^2 + 0.5).

\documentclass[tikz, border=24pt]{standalone} \usepackage{pgfplots} \usepgfplotslibrary{patchplots} \pgfplotsset{compat=1.18} \usepackage{luacode} % ------------------------------------------------------------ % ADAPTIVE MESH PARAMETERS % ------------------------------------------------------------ \def\AdaptiveDomain{{-7, 7, -7, 7}} % Singular/refinement points and initial kernel radius r around each point, format: {{x, y, r}, ...} \def\AdaptiveRefinementPoints{{{-2, -2, 1.5}, {2, 2, 1.0}}} \def\AdaptiveMaxIter{6} % number of refinement iterations \def\AdaptiveHmin{0.01} % fine mesh size near refinement points \def\AdaptiveHmax{2.0} % coarse mesh size far from features \def\AdaptiveDecay{4.0} % smoothness of transition (smaller = sharper transition) \def\AdaptiveBoundaryPoints{6} % number of points per edge of the rectangular domain boundary \def\AdaptiveDebug{true} % enable debug output \begin{luacode*} -- Helper functions -- Parse LaTeX array string into Lua table local function from_latex_array(str) str = str:gsub("^%s*", ""):gsub("%s*$", "") if str:sub(1,1) == "{" and str:sub(-1) == "}" then str = "return " .. str local f, err = load(str) if f then return f() else error("Parsing error: " .. tostring(err)) end else error("Expected a table in {...}") end end local ADAPTIVE_PARAMS = { domain = from_latex_array(token.get_macro("AdaptiveDomain")), refinement_points = from_latex_array(token.get_macro("AdaptiveRefinementPoints")), max_iter = tonumber(token.get_macro("AdaptiveMaxIter")), h_min = tonumber(token.get_macro("AdaptiveHmin")), h_max = tonumber(token.get_macro("AdaptiveHmax")), decay = tonumber(token.get_macro("AdaptiveDecay")), boundary_points_per_edge = tonumber(token.get_macro("AdaptiveBoundaryPoints")), debug = (token.get_macro("AdaptiveDebug") == "true"), show_tri_ids = (token.get_macro("AdaptiveShowTriIDs") == "true") } -- === Geometric and meshing utility functions === -- Parse LaTeX array string into Lua table local function from_latex_array(str) str = str:gsub("^%s*", ""):gsub("%s*$", "") if str:sub(1,1) == "{" and str:sub(-1) == "}" then str = "return " .. str local f, err = load(str) if f then return f() else error("Parsing error: " .. tostring(err)) end else error("Expected a table in {...}") end end -- Check if point p lies inside the circumcircle of triangle (a,b,c) function inCircumcircle(points, a, b, c, p) local ax, ay = points[a][1], points[a][2] local bx, by = points[b][1], points[b][2] local cx, cy = points[c][1], points[c][2] local px, py = points[p][1], points[p][2] local dx = ax - px local dy = ay - py local ex = bx - px local ey = by - py local fx = cx - px local fy = cy - py local ap = dx * dx + dy * dy local bp = ex * ex + ey * ey local cp = fx * fx + fy * fy local det = dx * (ey * cp - bp * fy) - dy * (ex * cp - bp * fx) + ap * (ex * fy - ey * fx) local orient = (bx - ax) * (cy - ay) - (by - ay) * (cx - ax) if math.abs(orient) < 1e-14 then return false -- degenerate triangle end return det * orient > 0 end -- Compute circumcenter and circumradius of triangle (a,b,c) function circumcenter_and_radius(points, a, b, c) local ax, ay = points[a][1], points[a][2] local bx, by = points[b][1], points[b][2] local cx, cy = points[c][1], points[c][2] local d = 2 * (ax*(by - cy) + bx*(cy - ay) + cx*(ay - by)) if math.abs(d) < 1e-14 then return nil, nil end local a2 = ax*ax + ay*ay local b2 = bx*bx + by*by local c2 = cx*cx + cy*cy local ux = (a2*(by - cy) + b2*(cy - ay) + c2*(ay - by)) / d local uy = (a2*(cx - bx) + b2*(ax - cx) + c2*(bx - ax)) / d local R = math.sqrt((ux - ax)^2 + (uy - ay)^2) return {ux, uy}, R end -- Compute target edge length at (x,y) based on distance to refinement points function target_h(x, y, ref_points) if #ref_points == 0 then return ADAPTIVE_PARAMS.h_max end local min_d = math.huge for _, p in ipairs(ref_points) do local dx = x - p[1] local dy = y - p[2] local d = math.sqrt(dx*dx + dy*dy) if d < min_d then min_d = d end end local h_min = ADAPTIVE_PARAMS.h_min local h_max = ADAPTIVE_PARAMS.h_max local decay = ADAPTIVE_PARAMS.decay return h_min + (h_max - h_min) * (1 - math.exp(-min_d / decay)) end -- Create a large enclosing super-triangle and return indices of its vertices function makeSuperTriangleIndices(points) local minx, miny = points[1][1], points[1][2] local maxx, maxy = minx, miny for _, p in ipairs(points) do minx = math.min(minx, p[1]) maxx = math.max(maxx, p[1]) miny = math.min(miny, p[2]) maxy = math.max(maxy, p[2]) end local d = math.max(maxx - minx, maxy - miny) * 10 local cx = (minx + maxx) / 2 local cy = (miny + maxy) / 2 -- Add super-triangle vertices to point list and return their indices table.insert(points, {cx - d, cy - d}) -- s1 table.insert(points, {cx + d, cy - d}) -- s2 table.insert(points, {cx, cy + d}) -- s3 return #points - 2, #points - 1, #points end -- Perform Delaunay triangulation using the Bowyer–Watson incremental algorithm function delaunay(points) -- Add super-triangle and get its vertex indices local s1, s2, s3 = makeSuperTriangleIndices(points) local triangles = {{s1, s2, s3}} -- Process only original points (before super-triangle was added) local original_n = #points - 3 for p_idx = 1, original_n do -- 1) find "bad" triangles whose circumcircle contains p_idx local badTriangles = {} for i, tri in ipairs(triangles) do if inCircumcircle(points, tri[1], tri[2], tri[3], p_idx) then table.insert(badTriangles, i) end end -- 2) count edge occurrences (canonical ordering) local edgeCount = {} for _, i in ipairs(badTriangles) do local tri = triangles[i] local edges = { {tri[1], tri[2]}, {tri[2], tri[3]}, {tri[3], tri[1]} } for _, e in ipairs(edges) do local u, v = e[1], e[2] if u > v then u, v = v, u end local key = string.format("%d-%d", u, v) edgeCount[key] = (edgeCount[key] or 0) + 1 end end -- 3) extract boundary edges (appear only once) local boundary = {} for key, cnt in pairs(edgeCount) do if cnt == 1 then local u, v = string.match(key, "(%d+)-(%d+)") if u then table.insert(boundary, {tonumber(u), tonumber(v)}) end end end table.sort(boundary, function(a,b) if a[1] ~= b[1] then return a[1] < b[1] end return a[2] < b[2] end) -- 4) remove bad triangles (in reverse order to preserve indices) table.sort(badTriangles, function(a,b) return a > b end) for _, i in ipairs(badTriangles) do table.remove(triangles, i) end -- 5) re-triangulate cavity with new point (ensure CCW orientation) for _, e in ipairs(boundary) do local a, b = e[1], e[2] local orient = (points[b][1]-points[a][1])*(points[p_idx][2]-points[a][2]) - (points[b][2]-points[a][2])*(points[p_idx][1]-points[a][1]) if orient > 0 then table.insert(triangles, {a, b, p_idx}) else table.insert(triangles, {b, a, p_idx}) end end end -- Remove triangles that include super-triangle vertices local filtered = {} for _, tri in ipairs(triangles) do local keep = true for _, v in ipairs(tri) do if v == s1 or v == s2 or v == s3 then keep = false break end end if keep then table.insert(filtered, tri) end end return filtered end -- === Adaptive Delaunay mesh generation === -- Generate an adaptive triangular mesh over a rectangular domain function adaptive_delaunay(ref_points, domain, max_iter) local xmin, xmax = domain[1], domain[2] local ymin, ymax = domain[3], domain[4] local n_edge = ADAPTIVE_PARAMS.boundary_points_per_edge local points = { {xmin, ymin}, {xmax, ymin}, {xmax, ymax}, {xmin, ymax}, {(xmin + xmax) / 2, (ymin + ymax) / 2} } local function inside(x, y) return (x >= xmin and x <= xmax and y >= ymin and y <= ymax) end local function exists_close(pt, set, tol) for _, q in ipairs(set) do local dx, dy = pt[1]-q[1], pt[2]-q[2] if dx*dx + dy*dy < tol*tol then return true end end return false end -- (1) seed refinement points and their hexagonal rings for idx, p in ipairs(ref_points) do local x0, y0, r if #p >= 3 then x0, y0, r = p[1], p[2], p[3] else x0, y0 = p[1], p[2] r = ADAPTIVE_PARAMS.h_max * 0.3 end if not inside(x0, y0) then if ADAPTIVE_PARAMS.debug then print("DEBUG: Refinement point", idx, "outside domain, skipped") end goto continue_point end if not exists_close({x0, y0}, points, r * 0.3) then table.insert(points, {x0, y0}) if ADAPTIVE_PARAMS.debug then print("DEBUG: Added center at (", x0, ",", y0, ") with r =", r) end end for k = 0, 5 do local angle = k * math.pi / 3 local x = x0 + r * math.cos(angle) local y = y0 + r * math.sin(angle) if inside(x, y) and not exists_close({x, y}, points, r * 0.3) then table.insert(points, {x, y}) end end ::continue_point:: end -- Add boundary helper points along each edge for i = 1, n_edge do local t = i / (n_edge + 1) table.insert(points, {xmin + t*(xmax - xmin), ymin}) table.insert(points, {xmax, ymin + t*(ymax - ymin)}) table.insert(points, {xmin + t*(xmax - xmin), ymax}) table.insert(points, {xmin, ymin + t*(ymax - ymin)}) end -- Iterative refinement loop for iter = 1, max_iter do -- Triangulate current point set (copy to avoid modifying original) local work_points = {} for i, p in ipairs(points) do work_points[i] = {p[1], p[2]} end local triangles = delaunay(work_points) local new_points = {} for _, tri in ipairs(triangles) do -- Skip triangles involving super-triangle vertices local a, b, c = tri[1], tri[2], tri[3] if a > #points or b > #points or c > #points then goto continue end local center, R = circumcenter_and_radius(points, a, b, c) if not center or R > 1e6 then goto continue end if center[1] < xmin or center[1] > xmax or center[2] < ymin or center[2] > ymax then goto continue end local h_loc = target_h(center[1], center[2], ref_points) local R_target = h_loc / math.sqrt(3) local merge_tol = math.max(ADAPTIVE_PARAMS.h_min * 0.5, h_loc * 0.6) -- Insert circumcenter if triangle is too large if R > R_target then if (not exists_close(center, points, merge_tol)) and (not exists_close(center, new_points, merge_tol)) then table.insert(new_points, center) end end ::continue:: end if #new_points == 0 then if ADAPTIVE_PARAMS.debug then print("Converged at iter", iter) end break end for _, p in ipairs(new_points) do table.insert(points, p) end if ADAPTIVE_PARAMS.debug then print("Iter", iter, "added:", #new_points) end end -- Final clean triangulation and filtering local n_final = #points local work_points = {} for i = 1, n_final do work_points[i] = {points[i][1], points[i][2]} end local raw_triangles = delaunay(work_points) -- Keep only valid triangles inside domain and with reasonable coordinates local final_triangles = {} for _, tri in ipairs(raw_triangles) do if tri[1] <= n_final and tri[2] <= n_final and tri[3] <= n_final then local a, b, c = points[tri[1]], points[tri[2]], points[tri[3]] if math.abs(a[1]) < 1e4 and math.abs(a[2]) < 1e4 and math.abs(b[1]) < 1e4 and math.abs(b[2]) < 1e4 and math.abs(c[1]) < 1e4 and math.abs(c[2]) < 1e4 then table.insert(final_triangles, tri) end end end return final_triangles, points end -- === Export to pgfplots-compatible file === -- Export adaptive mesh as a 3D surface (x y z) to 'surface.dat' function exportSurfaceToPgfplots() local domain = ADAPTIVE_PARAMS.domain local ref_points = ADAPTIVE_PARAMS.refinement_points local max_iter = ADAPTIVE_PARAMS.max_iter print("Generating adaptive mesh for 3D surface...") local triangles, points = adaptive_delaunay(ref_points, domain, max_iter) print("Points:", #points, "Triangles:", #triangles) local function z(u, v) local z1 = 3.0 / ((u + 2.0)^2 + (v + 2.0)^2 + 0.5) local z2 = 4.0 / ((u - 2.0)^2 + (v - 2.0)^2 + 0.5) return z1 + z2 end local file = io.open("surface.dat", "w") if not file then error("Failed to create surface.dat") end for _, tri in ipairs(triangles) do for _, idx in ipairs(tri) do local pt = points[idx] local x, y = pt[1], pt[2] local zz = z(x, y) file:write(string.format("%.6f %.6f %.6f\n", x, y, zz)) end end file:close() print("File surface.dat written successfully.") end \end{luacode*} % Генерация данных \directlua{exportSurfaceToPgfplots()} \begin{document} \begin{tikzpicture} \begin{axis}[ view={60}{30}, hide axis, % colormap/viridis, width=12cm, height=9cm, zmin=0, zmax=6.5, % grid=major, % ticklabel style={font=\footnotesize}, ] \addplot3[ patch, patch type=triangle, shader=flat, draw=black!40, line width=0.05pt, fill=red!60!black, fill opacity=0.85 ] table {surface.dat}; \end{axis} \end{tikzpicture} \end{document}

P.S. Excellent external meshing tools such as Triangle, Gmsh, etc., certainly exist and can be used to generate adaptive meshes for LaTeX plots. However, that workflow typically involves two separate stages: (1) mesh generation with an external program, and (2) plotting in LaTeX. In contrast, the solution presented here performs both mesh generation and visualization entirely within the LaTeX ecosystem, using only Lua and pgfplots.

Edit
Another example. This is a plot of the surface defined by the equation
z = u * e^{-u^2 - v^2}. This is first-order derivative of the 2D Gaussian in the u-direction.

To obtain the surface shown in the figure, the following modifications must be made to the code above:

Replace the adaptive mesh parameters block with the following:

\def\AdaptiveDomain{{-2.5, 2.5, -2.5, 2.5}} % Singular/refinement points and initial kernel radius r around each point, format: {{x, y, r}, ...} \def\AdaptiveRefinementPoints{{{0, 0, 1.5},{0.7, 0.0, 1.5},{-0.7, 0.0, 1.5}}} \def\AdaptiveMaxIter{6} % number of refinement iterations \def\AdaptiveHmin{0.08} % fine mesh size near refinement points \def\AdaptiveHmax{1.0} % coarse mesh size far from features \def\AdaptiveDecay{4.0} % smoothness of transition (smaller = sharper transition) \def\AdaptiveBoundaryPoints{6} % number of points per edge of the rectangular domain boundary \def\AdaptiveDebug{true} % enable debug output

Replace the function definition local function z(u, v)...end with:

local function z(u, v) return u * math.exp(-u^2 - v^2) end

Replace the \begin{tikzpicture}...\end{tikzpicture} block with:

\begin{tikzpicture} \begin{axis}[ view={150}{20}, hide axis, width=10cm, height=8cm, zmin=-0.5, zmax=0.5, ] \addplot3[ patch, patch type=triangle, shader=flat, draw=black!40, line width=0.05pt, fill=blue!60!black, fill opacity=0.85 ] table {surface.dat}; \end{axis} \end{tikzpicture}

I love it, but is there a way to not have to manually set the refinement points? What if it were a surface that was bumpy in many places? — Jasper
– Jasper, Commented Oct 25 at 14:31
@Jasper Thank you so much. Of course, I was initially moving in this direction, as evident from my first reply, but the problem is that this is a rather complex task in computational mathematics and programming, and I’m not sufficiently skilled to implement it. I did dabble a bit in FEM some time ago, which is why I realized the best approach is interactive: the user selects special points, and the computer handles everything else. And if there are many special points, a sufficiently fine mesh everywhere becomes inevitable. — kabenyuk
– kabenyuk, Commented Oct 25 at 15:11

nidarfp · Accepted Answer · 2025-10-21 09:50:10Z

I don't have a general method to suggest. But in this example we can easily have facets of variable sizes because the surface has an axis of symmetry, so we can calculate a parametric curve with points that are not regularly distributed, and rotate this curve around the axis. I drew the curve on the surface as well as its points to show their distribution:

\documentclass[border=5pt]{standalone}% compile with lualatex only \usepackage[svgnames]{xcolor} \usepackage[3d]{luadraw}%https://github.com/pfradin/luadraw %https://tex.stackexchange.com/questions/752634/a-logistical-way-to-recursively-subdivide-triangles-into-6-sub-triangles-such \begin{document} \begin{luadraw}{name=surface} local g = graph3d:new{window3d={-5,5,-5,5,-1,6},adjust2d=true,bbox=false,size={10,10}, viewdir={-30,50}} g:Labelsize("footnotesize") local curve = parametric3d(function(t) return M(t,0,3/(t^2+0.5)) end, 0,7, 15)[1] -- The points on this curve are not regularly spaced. local S = rotline(curve,{Origin,vecK},-180,180) -- We rotate the curve around the axis (O,k) to obtain the surface. g:Dboxaxes3d({grid=true,gridcolor="gray",fillcolor="LightGray",xyzstep=2}) g:Dfacet(S, {usepalette={palAutumn,"z"},clip=true}) g:Dpolyline3d(curve,false,"DarkGreen,line width=0.8pt",true) g:Ddots3d(curve,"scale=0.75,DarkGreen",true) g:Dlabel3d("curve $M(t,0,\\frac3{t^2+0.5})$", M(5.5,0,0),{dir={vecJ,vecK},node_options="DarkGreen"}) g:Show() \end{luadraw} \end{document}

Why here need parametric3d(...)[1], what is [1] for? I only found that parametric3d return a table of 3D segments. BTW, glad to see some v2.2's new features here😀. — Explorer
– Explorer, Commented Oct 21 at 14:01
@Explorer parametric3d() returns a llst (table) of lists of 3d points, because a curve can be made up of several pieces (several connected components), so [1] means that we take the first component of the curve (and as the function is continuous, there is no other). — nidarfp
– nidarfp, Commented Oct 21 at 14:12
A little doubt about the mmanual "The function parametric3d(p,t1,t2,nbdots,discont,nbdiv) calculates the points of the curve and returns a 3D polygonal line (no drawing)." How can we know that "a 3D polygonal line" is actually "a list (table) of lists of 3d points"? I would just regard "a 3D polygonal line" as "a table of 3D points, which consist a whole 3D polygonal line"... — Explorer
– Explorer, Commented Oct 21 at 14:19
Doc. p. 11: A polygonal line is a list (table) of connected components, and a connected component is a list (table) of complex numbers representing the affixes of the points. For example, the statement: local L = { {Z(-4,0), Z(0,2), Z(1,3)}, {Z(0,0), Z(4,-2), Z(1,-1)}} creates a polygonal line with two components in a variable L. — nidarfp
– nidarfp, Commented Oct 21 at 15:18
More clearly now! I neglected that just now.... Shame on me! — Explorer
– Explorer, Commented Oct 21 at 16:10

kabenyuk · Accepted Answer · 2025-10-20 16:31:01Z

It seems to me that refining the mesh by splitting a triangle into six triangles is not a good solution. With this approach, as you get closer to surface features (e.g., a function peak), the triangle angles tend to zero, which leads to strong element degeneracy and forces increasingly high refinement density even in regions where the function varies only moderately.

Instead, I propose the following strategy. The initial parametric domain is first partitioned into a uniform triangulation (as in your approach), with its density controlled by the initdens parameter. Then each triangle is checked for the need of further refinement. The approach is inspired by ideas from the finite element method (FEM), where the quality of the numerical solution directly depends on the shape and local density of elements: in regions of high curvature or large gradients a finer mesh is required, whereas on smooth regions larger elements are sufficient.

When refining, a triangle is split into four triangles by its midlines. As a result, four triangles are obtained, each similar to the original and, crucially, with angles preserved—this prevents element degeneration under recursive subdivision.

The adaptation criterion is based on estimating the local approximation error—which, in my view, is fairly reliable. For each triangle, the true position of the point corresponding to the barycenter of the triangle in parameter space is computed and compared to the linear interpolation of its vertex coordinates. Importantly, the comparison is performed not only in the vertical z coordinate, but using the full three-dimensional Euclidean distance in R^3.

The process is controlled by three key parameters:

initdens — sets the initial density of the uniform mesh;
tolerance — the error threshold below which no refinement is required;
steps — the maximum depth of recursive subdivision.

Despite its theoretical appeal, this approach has a significant drawback: local refinement breaks mesh conformity. At the boundaries between refined and unrefined elements, so-called T-junctions arise, which leads to microscopic gaps—“lacunae”—between triangles and, in some cases, to partial overlaps. The magnitude of these artifacts depends on the relationship among the parameters above; for the illustration presented, they were chosen so that the gaps are barely noticeable.

I tried not to alter the original architecture of your code. A new implementation of local function append_surface() was added to the Lua file, and only minimal edits were made to the rest—both in LaTeX and Lua—to support the new variant. However, I had to remove both of your local function append_surface() definitions due to response length limitations.

mwe.tex

\documentclass[ tikz ,border = 1cm ]{standalone} \makeatletter \tikzset{ /lua-tikz3dtools/.is family % abbreviated "td" ,/lua-tikz3dtools/.cd ,/lua-tikz3dtools/parametric/.cd % appreviated "p" ,/lua-tikz3dtools/parametric/matrix/.cd % abbreviated "m" ,object/.code = {\protected@edef\luatikztdtools@p@m@transformation{#1}} ,object = {identity_matrix()} ,name/.code = {\protected@edef\luatikztdtools@p@m@name{#1}} ,name = {Change me!} ,/lua-tikz3dtools/parametric/surface/.cd % abbreviated "s" ,ustart/.estore in = \luatikztdtools@p@s@ustart ,ustop/.estore in = \luatikztdtools@p@s@ustop ,usamples/.estore in = \luatikztdtools@p@s@usamples ,vstart/.estore in = \luatikztdtools@p@s@vstart ,vstop/.estore in = \luatikztdtools@p@s@vstop ,vsamples/.estore in = \luatikztdtools@p@s@vsamples ,x/.code = {\protected@edef\luatikztdtools@p@s@x{#1}} ,y/.code = {\protected@edef\luatikztdtools@p@s@y{#1}} ,z/.code = {\protected@edef\luatikztdtools@p@s@z{#1}} ,transformation/.code = {\protected@edef\luatikztdtools@p@s@transformation{#1}} ,fill options/.code = {\protected@edef\luatikztdtools@p@s@filloptions{#1}} ,ustart = {0} ,ustop = {1} ,usamples = {10} ,vstart = {0} ,vstop = {1} ,vsamples = {10} ,x = {sphere(tau*u, pi*v)[1][1]} ,y = {sphere(tau*u, pi*v)[1][2]} ,z = {sphere(tau*u, pi*v)[1][3]} ,transformation = {identity_matrix()} ,fill options = {fill, fill opacity = 0.3} ,initdens/.estore in = \luatikztdtools@p@s@initdens ,tolerance/.estore in = \luatikztdtools@p@s@tolerance ,steps/.estore in = \luatikztdtools@p@s@steps } \ExplSyntaxOn \lua_load_module:n {mwe-lua} \NewDocumentCommand { \appendsurface } {o} { \group_begin: \tikzset{ /lua-tikz3dtools/parametric/surface/.search~also = {/tikz} ,/lua-tikz3dtools/parametric/surface/.cd ,#1 } \__lua_tikztdtools_appendsurface: \group_end: } \NewDocumentCommand { \displaysegments } {} { \group_begin: \__lua_tikztdtools_displaysegments: \group_end: } \NewDocumentCommand { \setobject } {o} { \group_begin: \tikzset{ /lua-tikz3dtools/parametric/matrix/.search~also = {/tikz} ,/lua-tikz3dtools/parametric/matrix/.cd ,#1 } \__lua_tikztdtools_setobject: \group_end: } \ExplSyntaxOff \makeatother \begin{document} \begin{tikzpicture} \setobject[ name = {T} ,object = { matrix_multiply( euler(pi/2,pi/3,5*pi/6) ,translate(0,0,-15) ) } ] \appendsurface[ ustart = -7 ,ustop = 7 % ,usamples = 7 ,vstart = -7 ,vstop = 7 % ,vsamples = 7 ,transformation = T ,initdens = 10 ,steps = 4 ,tolerance = 0.004 ,x = {u} ,y = {v} ,z = { 3 / ((u)^2+(v)^2+1/2) } ,fill options = { preaction = { fill = red!70!black ,fill opacity = 1 } ,postaction = { draw = black ,ultra thin ,line join = round ,line cap = round } } ] \displaysegments \end{tikzpicture} \end{document}

mve-lua.lua

 local eps = 0.0000001 --- scalar multiplication --- @param a number the scalar coefficient --- @param v table<table<number>> the vector --- @return table<table<number>> the scaled vector local function scalar_multiplication(a, v) local tmp = v[1] return { {a*tmp[1], a*tmp[2], a*tmp[3], 1} } end --- vector subtraction --- @param u table<table<number>> the first vector --- @param v table<table<number>> the second vector --- @return table<table<number>> their difference local function vector_subtraction(u, v) local U = u[1] local V = v[1] return { {U[1]-V[1], U[2]-V[2], U[3]-V[3], 1} } end --- vector addition --- @param u table<table<number>> the first vector --- @param v table<table<number>> the second vector --- @return table<table<number>> their sum local function vector_addition(u, v) local U = u[1] local V = v[1] return { {U[1]+V[1], U[2]+V[2], U[3]+V[3], 1} } end --- length of a vector --- @param u table<table<number>> the vector to be measured --- @return number the vectors length local function length(u) local result = math.sqrt((u[1][1])^2 + (u[1][2])^2 + (u[1][3])^2) return result end --- distance between two points --- @param P1 table<table<number>> a point --- @param P2 table<table<number>> another point --- @return number their distance local function distance(P1, P2) local A = P1[1] local B = P2[1] return math.sqrt( (A[1]-B[1])^2 + (A[2]-B[2])^2 + (A[3]-B[3])^2 ) end -- https://tex.stackexchange.com/a/747040 --- Creates a TeX command that evaluates a Lua function --- --- @param name string The name of the `\csname` to define --- @param func function --- @param args table<string> The TeX types of the function arguments --- @param protected boolean|nil Define the command as `\protected` --- @return nil local function register_tex_cmd(name, func, args, protected) -- The extended version of this function uses `N` and `w` where appropriate, -- but only using `n` is good enough for exposition purposes. name = "__lua_tikztdtools_" .. name .. ":" .. ("n"):rep(#args) -- Push the appropriate scanner functions onto the scanning stack. local scanners = {} for _, arg in ipairs(args) do scanners[#scanners+1] = token['scan_' .. arg] end -- An intermediate function that properly "scans" for its arguments -- in the TeX side. local scanning_func = function() local values = {} for _, scanner in ipairs(scanners) do values[#values+1] = scanner() end func(table.unpack(values)) end local index = luatexbase.new_luafunction(name) lua.get_functions_table()[index] = scanning_func if protected then token.set_lua(name, index, "protected") else token.set_lua(name, index) end end local segments = {} local mm = {} mm.tau = 2*math.pi local pi = math.pi local cos, sin = math.cos, math.sin --- matrix multiplication --- --- @param A table<table<number>> left matrix --- @param B table<table<number>> right matrix --- @return table<table<number>> the product function mm.matrix_multiply(A, B) local rows_A = #A local columns_A = #A[1] local rows_B = #B local columns_B = #B[1] assert( columns_A == rows_B, string.format( [[ Wrong size matrices for multiplication. Size A: %d,%d Size B: %d,%d ]], rows_A, columns_A, rows_B, columns_B ) ) local product = {} for row = 1, rows_A do product[row] = {} for column = 1, columns_B do product[row][column] = 0 for dot_product_step = 1, columns_A do local a = A[row][dot_product_step] local b = B[dot_product_step][column] assert(type(a) == "number", string.format("Expected number but got %s in A[%d][%d]", type(a), row, dot_product_step)) assert(type(b) == "number", string.format("Expected number but got %s in B[%d][%d]", type(b), dot_product_step, column)) product[row][column] = product[row][column] + a * b end end end return product end function mm.yrotation(angle) local c = cos(angle) local s = sin(angle) return { {c,0,-s,0} ,{0,1,0,0} ,{s,0,c,0} ,{0,0,0,1} } end function mm.translate(x,y,z) return { {1,0,0,0} ,{0,1,0,0} ,{0,0,1,0} ,{x,y,z,1} } end function mm.zrotation(angle) local c = cos(angle) local s = sin(angle) return { {c,s,0,0} ,{-s,c,0,0} ,{0,0,1,0} ,{0,0,0,1} } end function mm.euler(alpha,beta,gamma) return mm.matrix_multiply( mm.zrotation(gamma) ,mm.matrix_multiply( mm.yrotation(beta) ,mm.zrotation(alpha) ) ) end function mm.midpoint(triangle) local P,Q,R = table.unpack(triangle) local x = (P[1]+Q[1]+R[1])/3 local y = (P[2]+Q[2]+R[2])/3 local z = (P[3]+Q[3]+R[3])/3 return {{x,y,z,1}} end local lua_tikz3dtools_parametric = {} lua_tikz3dtools_parametric.math = {} for k, v in pairs(_G) do lua_tikz3dtools_parametric.math[k] = v end for k, v in pairs(mm) do lua_tikz3dtools_parametric.math[k] = v end for k, v in pairs(math) do lua_tikz3dtools_parametric.math[k] = v end local function single_string_expression(str) return load(("return %s"):format(str), "expression", "t", lua_tikz3dtools_parametric.math)() end local function double_string_function(str) return load(("return function(u,v) return %s end"):format(str), "expression", "t", lua_tikz3dtools_parametric.math)() end -- my append_surface local function append_surface(hash) local ustart = hash.ustart local ustop = hash.ustop local vstart = hash.vstart local vstop = hash.vstop local x = hash.x local y = hash.y local z = hash.z local filloptions = hash.filloptions local transformation = hash.transformation local max_depth = tonumber(hash.steps) or 1 local tolerance = tonumber(hash.tolerance) or 0.01 local initdens = tonumber(hash.initdens) or 2 x = double_string_function(x) y = double_string_function(y) z = double_string_function(z) ustart = single_string_expression(ustart) ustop = single_string_expression(ustop) vstart = single_string_expression(vstart) vstop = single_string_expression(vstop) transformation = single_string_expression(transformation) -- Parametric surface in homogeneous coords local function parametric_surface(u, v) return { x(u,v), y(u,v), z(u,v), 1 } end -- Estimate error: difference between true center value and linear interp local function needs_refinement(tri) local A, B, C = tri[1], tri[2], tri[3] local uA, vA = A[1], A[2] local uB, vB = B[1], B[2] local uC, vC = C[1], C[2] local P = parametric_surface(uA, vA) local Q = parametric_surface(uB, vB) local R = parametric_surface(uC, vC) local uG = (uA + uB + uC) / 3 local vG = (vA + vB + vC) / 3 local G_true = parametric_surface(uG, vG) local G_interp = { (P[1] + Q[1] + R[1]) / 3, (P[2] + Q[2] + R[2]) / 3, (P[3] + Q[3] + R[3]) / 3 } local err = distance({G_true}, {G_interp}) return err > tolerance end -- Initial triangulation of the domain rectangle local tris = {} local n = math.max(2, initdens) local du = (ustop - ustart) / (n - 1) local dv = (vstop - vstart) / (n - 1) for i = 0, n - 2 do for j = 0, n - 2 do local u0 = ustart + i * du local v0 = vstart + j * dv local u1 = u0 + du local v1 = v0 + dv -- Два треугольника на квадратную ячейку table.insert(tris, {{u0, v0}, {u1, v0}, {u0, v1}}) table.insert(tris, {{u1, v1}, {u1, v0}, {u0, v1}}) end end for depth = 1, max_depth do local newtris = {} local refined = false for _, tri in ipairs(tris) do if depth == max_depth or not needs_refinement(tri) then table.insert(newtris, tri) else local v1, v2, v3 = tri[1], tri[2], tri[3] -- Midpoints local e12 = { (v1[1] + v2[1]) * 0.5, (v1[2] + v2[2]) * 0.5 } local e23 = { (v2[1] + v3[1]) * 0.5, (v2[2] + v3[2]) * 0.5 } local e31 = { (v3[1] + v1[1]) * 0.5, (v3[2] + v1[2]) * 0.5 } -- 4-subtriangle refinement table.insert(newtris, {v1, e12, e31}) table.insert(newtris, {v2, e23, e12}) table.insert(newtris, {v3, e31, e23}) table.insert(newtris, {e12, e23, e31}) refined = true end end tris = newtris if not refined then break end end -- Project to 3D, transform, and store for _, tri in ipairs(tris) do local A = parametric_surface(tri[1][1], tri[1][2]) local B = parametric_surface(tri[2][1], tri[2][2]) local C = parametric_surface(tri[3][1], tri[3][2]) local the_segment = mm.matrix_multiply({A, B, C}, transformation) if distance({A}, {B}) >= eps and distance({B}, {C}) >= eps and distance({A}, {C}) >= eps then table.insert(segments, { segment = the_segment, filloptions = filloptions, type = "triangle" }) end end end register_tex_cmd( "appendsurface", function() append_surface{ ustart = token.get_macro("luatikztdtools@p@s@ustart"), ustop = token.get_macro("luatikztdtools@p@s@ustop"), usamples = token.get_macro("luatikztdtools@p@s@usamples"), vstart = token.get_macro("luatikztdtools@p@s@vstart"), vstop = token.get_macro("luatikztdtools@p@s@vstop"), vsamples = token.get_macro("luatikztdtools@p@s@vsamples"), x = token.get_macro("luatikztdtools@p@s@x"), y = token.get_macro("luatikztdtools@p@s@y"), z = token.get_macro("luatikztdtools@p@s@z"), transformation = token.get_macro("luatikztdtools@p@s@transformation"), filloptions = token.get_macro("luatikztdtools@p@s@filloptions"), steps = token.get_macro("luatikztdtools@p@s@steps"), tolerance = token.get_macro("luatikztdtools@p@s@tolerance"), initdens = token.get_macro("luatikztdtools@p@s@initdens"), } end, { } ) local function make_object(name, tbl) local obj = {} obj[name] = tbl return obj end local function set_matrix(hash) local transformation = single_string_expression(hash.transformation) local name = hash.name local tbl = make_object(name, transformation) for k, v in pairs(tbl) do lua_tikz3dtools_parametric.math[k] = v end end register_tex_cmd( "setobject", function() set_matrix{ name = token.get_macro("luatikztdtools@p@m@name"), transformation = token.get_macro("luatikztdtools@p@m@transformation"), } end, { } ) local function display_segments() local function comparator(a,b) return mm.midpoint(a.segment)[1][3] < mm.midpoint(b.segment)[1][3] end -- Perform sorting after filtering out skipped segments table.sort(segments, comparator) -- Iterate through the sorted segments for rendering for _, segment in ipairs(segments) do -- Handle triangle segments if segment.type == "triangle" then local P, Q, R = segment.segment[1], segment.segment[2], segment.segment[3] tex.sprint(string.format( "\\path[%s] (%f,%f) -- (%f,%f) -- (%f,%f) -- cycle;", segment.filloptions, P[1], P[2], Q[1], Q[2], R[1], R[2] )) end end -- Clear the segments array after processing segments = {} end register_tex_cmd("displaysegments", function() display_segments() end, { })

Wow! I am very impressed with your work, as usual. Expect a bounty when the timer allows me to. If I may ask, why does the more dense tessellation seem to occur along one of the long diagonals, and not the other? Wouldn't you expect it to be more radially centred around the origin? — Jasper
– Jasper, Commented Oct 20 at 17:22
One problem I've noticed: Triangles which are adjacent to more subdivided ones can have gaps when the subdivided ones interpolate an interior point along a line segment. — Jasper
– Jasper, Commented Oct 20 at 18:20
Thanks for the kind words! The diagonal bias comes from the rectangular initial grid—it doesn’t perfectly capture the radial symmetry of the function. The gaps appear because refinement is non‑conforming: when one triangle is split but its neighbour isn’t, T‑junctions form. A fully conforming mesh (where every shared vertex belongs to all adjacent triangles) would fix this, but it requires maintaining mesh topology, which greatly complicates the code. I’m still thinking about these issues. It would be worth looking into how this is handled (or isn’t) in FEM. — kabenyuk
– kabenyuk, Commented Oct 21 at 4:36
+1. Sometimes I wonder if people will implement a computer algebra system entirely in Lua, LaTeX and TikZ!! — Apoorv Potnis
– Apoorv Potnis, Commented Oct 25 at 12:01
@ApoorvPotnis Thanks! Well, according to the Church–Turing thesis, if something is computable at all, it can be done in Lua—after all, Lua is Turing-complete. Whether it should be done in TikZ… that’s a question for the typesetting gods. 😅 — kabenyuk
– kabenyuk, Commented Oct 25 at 12:14

Stack Exchange Network

A logistical way to recursively subdivide triangles (into 6 sub-triangles), such that all detail is captured and no triangles are wasted

mwe.tex

mwe-lua.lua

3 Answers 3

You must log in to answer this question.

Linked

Hot Network Questions