Local Grid Refinement: function unsuitable now can be used with only 2 arguments

CHANGELOG.md

@@ -1,4 +1,10 @@
 # Changes
+## v2.8.0 October 13, 2026
+
+### Added
+
+ - Local Grid Refinement: function `unsuitable` now can be used with only 2 arguments (barycentric coordinates, area/volume) in 2D and 3D. This unifies the signature for both 2D and 3D.
+ - New doc page: `"Local Grid Refinement"` and tests.
 
 ## v2.7.0 October 12, 2026
 

Project.toml

@@ -1,6 +1,6 @@
 name = "SimplexGridFactory"
 uuid = "57bfcd06-606e-45d6-baf4-4ba06da0efd5"
-version = "2.7.0"
+version = "2.8.0"
 authors = ["Juergen Fuhrmann <juergen.fuhrmann@wias-berlin.de>", "Patrick Jaap <patrick.jaap@wias-berlin.de>", "Jan Philipp Thiele  <jan-philipp.thiele@tu-braunschweig.de>"]
 
 [deps]
@@ -12,6 +12,7 @@ GridVisualize = "5eed8a63-0fb0-45eb-886d-8d5a387d12b8"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MeshIO = "7269a6da-0436-5bbc-96c2-40638cbb6118"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [weakdeps]
 TetGen = "c5d3f3f7-f850-59f6-8a2e-ffc6dc1317ea"
@@ -26,6 +27,7 @@ GridVisualize = "1.7.0"
 LinearAlgebra = "1.9"
 MeshIO = "0.4, 0.5"
 Printf = "1.6"
+StaticArrays = "1.9.23"
 TetGen = "1.5,2,3"
 Triangulate = "2.3.1, 3"
 julia = "1.9"

docs/make.jl

@@ -26,6 +26,7 @@ function mkdocs()
             "Home" => "index.md",
             "Changes" => "changes.md",
             "API" => "api.md",
+            "Local Grid Refinement" => "local_refinement.md",
             "Examples" => generated_examples,
             "Notebooks" => notebook_examples,
             "Internals" => "internals.md",

docs/src/local_refinement.md

@@ -0,0 +1,64 @@
+## Local refinement
+
+### 2D - Triangulate
+
+The `unsuitable` mechanism is used to refine the 2D triangulation on a domain, for instance surrounding a geometric object (e.g., a point or line). 
+Refinement is triggered when a condition for a given triangle is satisfied.
+For expamle when the distance between a point of interest and the center (barycentric coordinates) of the triangle falls below a threshold defined by the triangle’s area.
+#### Example usage:
+
+- One could define the `unsuitable` function with 6 points and the area of the triangle (x1, y1, x2, y2, x3, y3, area). This is the default signature in [Triangulate.jl](https://github.com/JuliaGeometry/Triangulate.jl):
+```julia
+    refinement_center = [0.5, 0.5]
+    function unsuitable(x1, y1, x2, y2, x3, y3, area)
+        bary = [(x1 + x2 + x3) / 3, (y2 + y2 + y3) / 3]
+        dist = norm(bary - refinement_center)
+        return area > 0.01 * dist ? 1 : 0
+    end
+```  
+
+- As a simpler alternative, one could directly use the barycentric coordinates and the area of the triangle (center, area): 
+```julia
+    refinement_center = [0.5, 0.5]
+    function unsuitable(center, area)
+        dist = norm(center - refinement_center)
+        return area > 0.01 * dist ? 1 : 0
+    end
+```
+To hand this function over to the builder:     
+```julia 
+    options!(builder; unsuitable = unsuitable)
+```
+
+
+### 3D - TetGen
+
+The  `unsuitable` mechanism is also used to refine the 3D triangulation on a domain, for instance surrounding a geometric object (e.g., a point or line). Refinement is triggered when a condition for a given tetrahedron is satisfied.
+For example when the distance between the point of interest and the barycentric coordinates of a tetrahedron falls below a threshold defined by the tetrahedron’s volume.   
+
+#### Example usage:
+
+- One could define the `unsuitable` function with 4 points (pa, pb, pc, pd). 
+This is the default signature in [TetGen.jl](https://github.com/JuliaGeometry/TetGen.jl):  
+```julia
+    refinement_center = [0.5, 0.5, 0.5]
+    function unsuitable(pa, pb, pc, pd)
+        vol = det(hcat(pb - pa, pc - pa, pd - pa)) / 6
+        center = 0.25 * (pa + pb + pc + pd) 
+        return vol > 0.05 * norm(center - refinement_center)^2.5
+    end
+```
+
+
+- As a simpler alternative, one could directly use the barycentric coordinates and the volume of the tetrahedron (center, volume):  
+```julia
+    refinement_center = [0.5, 0.5, 0.5]
+    function unsuitable(center, vol)
+        return vol > 0.05 * norm(center - refinement_center)^2.5
+    end
+
+```
+   To hand this function over to the builder:     
+```julia 
+    options!(builder; unsuitable = unsuitable)
+```

examples/examples2d.jl

@@ -168,9 +168,8 @@ function square_localref()
     facet!(builder, p4, p1)
 
     refinement_center = [0.5, 0.5]
-    function unsuitable(x1, y1, x2, y2, x3, y3, area)
-        bary = [(x1 + x2 + x3) / 3, (y2 + y2 + y3) / 3]
-        dist = norm(bary - refinement_center)
+    function unsuitable(center, area)
+        dist = norm(center - refinement_center)
         if area > 0.01 * dist
             return 1
         else

src/SimplexGridFactory.jl

@@ -5,16 +5,16 @@ $(read(joinpath(@__DIR__, "..", "README.md"), String))
 """
 module SimplexGridFactory
 
-using Printf: @sprintf
-using LinearAlgebra: norm
 
-using ElasticArrays: ElasticArray
-import ExtendableGrids
-using ExtendableGrids: dim_space, Coordinates, BFaceNodes, BFaceRegions, simplexgrid, BinnedPointList
 using DocStringExtensions: SIGNATURES, TYPEDEF, TYPEDSIGNATURES
-using GridVisualize: GridVisualize, GridVisualizer, plot_triangulateio!
+using ElasticArrays: ElasticArray
+using ExtendableGrids: ExtendableGrids, dim_space, Coordinates, BFaceNodes, BFaceRegions, simplexgrid, BinnedPointList
 using FileIO: load
-import MeshIO
+using GridVisualize: GridVisualize, GridVisualizer, plot_triangulateio!
+using LinearAlgebra: norm, det
+using MeshIO: MeshIO
+using Printf: @sprintf
+using StaticArrays: @MArray
 
 include("utils.jl")
 

src/plot.jl

@@ -59,7 +59,7 @@ function builderplot(
     end
 
     if !isnothing(opts[:unsuitable])
-        Triangulate.triunsuitable!(opts[:unsuitable])
+        set_unsuitable!(TriangulateType, Triangulate, opts[:unsuitable])
     end
 
     triout, vorout = Triangulate.triangulate(flags, triin)

src/tetgen.jl

@@ -8,15 +8,14 @@ See [`default_options`](@ref) for available `kwargs`.
 """
 function ExtendableGrids.simplexgrid(::Type{TetGenType}, TetGen, input; kwargs...)
     opts = blendoptions!(default_options(:tetgen); kwargs...)
-
     flags = makeflags(opts, :tetgen)
 
     if opts[:verbose]
         @show flags
     end
 
     if !isnothing(opts[:unsuitable])
-        TetGen.tetunsuitable!(opts[:unsuitable])
+        set_unsuitable!(TetGenType, TetGen, opts[:unsuitable])
     end
 
     tetout = TetGen.tetrahedralize(input, flags)

src/triangle.jl

@@ -14,7 +14,7 @@ function ExtendableGrids.simplexgrid(::Type{TriangulateType}, Triangulate, input
         @show flags
     end
     if !isnothing(opts[:unsuitable])
-        Triangulate.triunsuitable!(opts[:unsuitable])
+        set_unsuitable!(TriangulateType, Triangulate, opts[:unsuitable])
     end
 
     triout, vorout = Triangulate.triangulate(flags, input)

src/utils.jl

@@ -15,3 +15,50 @@ istriangulate(g) = typeof(g) == Module && nameof(g) == :Triangulate
 abstract type TetGenType end
 
 abstract type TriangulateType end
+
+function set_unsuitable!(::Type{TriangulateType}, Generator, unsuitable)
+
+    # We assume that the "Triangle" uses Cdouble internally
+    if hasmethod(unsuitable, Tuple{Vector{Float64}, Float64})
+
+        center = zeros(2)
+
+        new_unsuitable = (x1, y1, x2, y2, x3, y3, area) -> begin
+
+            # center of the triangle
+            center[1] = (x1 + x2 + x3) / 3
+            center[2] = (y1 + y2 + y3) / 3
+
+            return unsuitable(center, area)
+        end
+    else
+        new_unsuitable = unsuitable
+    end
+
+    return Generator.triunsuitable!(new_unsuitable)
+
+end
+
+function set_unsuitable!(::Type{TetGenType}, Generator, unsuitable)
+
+    # We assume that the "TetGen" uses Cdouble internally
+    if hasmethod(unsuitable, Tuple{Vector{Float64}, Float64})
+        center = zeros(3)
+        matrix = @MArray zeros(3, 3)
+
+        new_unsuitable = (pa, pb, pc, pd) -> begin
+            @. center = 0.25 * (pa + pb + pc + pd)
+            @. matrix[:, 1] = pb - pa
+            @. matrix[:, 2] = pc - pa
+            @. matrix[:, 3] = pd - pa
+
+            vol = abs(det(matrix)) / 6
+            return unsuitable(center, vol)
+        end
+    else
+        new_unsuitable = unsuitable
+    end
+
+    return Generator.tetunsuitable!(new_unsuitable)
+
+end

test/runtests.jl

@@ -107,7 +107,7 @@ end
 
 function test_triunsuitable(x1, y1, x2, y2, x3, y3, area)
     refinement_center = [0.5, 0.5]
-    bary = [(x1 + x2 + x3) / 3, (y2 + y2 + y3) / 3]
+    bary = [(x1 + x2 + x3) / 3, (y1 + y2 + y3) / 3]
     dist = norm(bary - refinement_center)
     if area > 0.01 * dist
         return 1
@@ -116,6 +116,12 @@ function test_triunsuitable(x1, y1, x2, y2, x3, y3, area)
     end
 end
 
+refinement_center = [0.5, 0.5]
+function test_triunsuitable_comfi(center, area)
+    dist = norm(center - refinement_center)
+    return area > 0.01 * dist ? 1 : 0
+end
+
 @testset "Simplexgrid (arrays 2d)" begin
     function test_simplesquare(; kwargs...)
         tio = SimplexGridFactory.triangulateio(
@@ -137,6 +143,7 @@ end
     @test testgrid(test_simplesquare(; minangle = 30), (91, 148, 32))
 
     @test testgrid(test_simplesquare(; unsuitable = test_triunsuitable), (299, 550, 46))
+    @test testgrid(test_simplesquare(; unsuitable = test_triunsuitable_comfi), (299, 550, 46))
 end
 
 @testset "SimplexGridBuilder 2d" begin
@@ -191,12 +198,18 @@ end
     @test testgrid(test_buildersquare1(), (299, 550, 46))
 end
 
+
 function test_tetunsuitable(pa, pb, pc, pd)
     vol = det(hcat(pb - pa, pc - pa, pd - pa)) / 6
     center = 0.25 * (pa + pb + pc + pd) - [0.5, 0.5, 0.5]
     return vol > 0.05 * norm(center)^2.5
 end
 
+refinement_center = [0.5, 0.5, 0.5]
+function test_tetunsuitable_comfi(center, vol)
+    return vol > 0.05 * norm(center - refinement_center)^2.5
+end
+
 @testset "Simplexgrid (arrays 3d)" begin
     function test_simplecube(; kwargs...)
         tio = SimplexGridFactory.tetgenio(
@@ -236,6 +249,7 @@ end
         @test testgrid(test_simplecube(; flags = "pAaqQD"), (63, 192, 96))
         @test testgrid(test_simplecube(; maxvolume = 0.05), (14, 24, 24))
         @test testgrid(test_simplecube(; unsuitable = test_tetunsuitable), (167, 788, 110))
+        @test testgrid(test_simplecube(; unsuitable = test_tetunsuitable_comfi), (167, 788, 110))
 
     end
 end