Add support for predicates in TriRefinement (#1070)

* Add support for predicates in 'TriRefinement'

* Apply suggestions

* Apply suggestions

* Apply suggestions

* Apply suggestions

* Rename variables

* Fix typo

* Rename variables

* Update src/refinement/tri.jl

Co-authored-by: Júlio Hoffimann <[email protected]>

* Add tests

* Add more tests

* Update code order

* Update code order

---------

Co-authored-by: Júlio Hoffimann <[email protected]>

Author: eliascarv

src/refinement/tri.jl

@@ -3,14 +3,20 @@
 # ------------------------------------------------------------------
 
 """
-    TriRefinement()
+    TriRefinement([pred])
 
 Refinement of polygonal meshes into triangles.
-A n-gon is subdivided into n triangles.
+A n-gon for which the predicate `pred` holds true
+is subdivided into n triangles. The method refines all
+n-gons if the `pred` is ommited.
 """
-struct TriRefinement <: RefinementMethod end
+struct TriRefinement{F} <: RefinementMethod
+  pred::F
+end
+
+TriRefinement() = TriRefinement(nothing)
 
-function refine(mesh, ::TriRefinement)
+function refine(mesh, method::TriRefinement)
   assertion(paramdim(mesh) == 2, "TriRefinement only defined for surface meshes")
   (eltype(mesh) <: Triangle) || return simplexify(mesh)
 
@@ -21,9 +27,19 @@ function refine(mesh, ::TriRefinement)
   # convert to half-edge structure
   t = convert(HalfEdgeTopology, connec)
 
+  # indices to refine
+  rinds = if isnothing(method.pred)
+    1:nelements(t)
+  else
+    filter(i -> method.pred(mesh[i]), 1:nelements(t))
+  end
+
+  # indices to preserve
+  pinds = setdiff(1:nelements(t), rinds)
+
   # add centroids of elements
   ∂₂₀ = Boundary{2,0}(t)
-  epts = map(1:nelements(t)) do elem
+  rpts = map(rinds) do elem
     is = ∂₂₀(elem)
     cₒ = sum(i -> to(points[i]), is) / length(is)
     withcrs(mesh, cₒ)
@@ -33,23 +49,31 @@ function refine(mesh, ::TriRefinement)
   vpts = points
 
   # new points in refined mesh
-  newpoints = [vpts; epts]
+  newpoints = [vpts; rpts]
 
+  # new connectivities in refined mesh
+  newconnec = Connectivity{Triangle,3}[]
+
+  # offset to new vertex indices
   offset = length(vpts)
 
-  # connect vertices into new triangles
-  newconnec = Connectivity{Triangle,3}[]
-  for elem in 1:nelements(t)
+  # connectivities of new triangles
+  for (i, elem) in enumerate(rinds)
     verts = ∂₂₀(elem)
     nv = length(verts)
-    for i in 1:nv
-      u = elem + offset
-      v = verts[mod1(i, nv)]
-      w = verts[mod1(i + 1, nv)]
+    for j in 1:nv
+      u = i + offset
+      v = verts[mod1(j, nv)]
+      w = verts[mod1(j + 1, nv)]
       tri = connect((u, v, w))
       push!(newconnec, tri)
     end
   end
 
+  # connectivities of preserved elements
+  for elem in pinds
+    push!(newconnec, element(t, elem))
+  end
+
   SimpleMesh(newpoints, newconnec)
 end

test/refinement.jl

@@ -15,6 +15,17 @@
   grid = CartesianGrid((3, 3), merc(0, 0), (T(1), T(1)))
   ref = refine(grid, TriRefinement())
   @test crs(ref) === crs(grid)
+
+  # predicate
+  points = cart.([(0, 0), (4, 0), (8, 0), (3, 1), (5, 1), (2, 2), (4, 2), (6, 2), (4, 4)])
+  connec = connect.([(1, 2, 6), (2, 3, 8), (6, 8, 9), (2, 5, 4), (4, 5, 7), (4, 7, 6), (5, 8, 7)])
+  mesh = SimpleMesh(points, connec)
+  ref = refine(mesh, TriRefinement(e -> measure(e) > T(1) * u"m^2"))
+  @test nelements(ref) == 13
+  @test nvertices(ref) == 12
+  ref = refine(mesh, TriRefinement(e -> measure(e) ≤ T(1) * u"m^2"))
+  @test nelements(ref) == 15
+  @test nvertices(ref) == 13
 end
 
 @testitem "QuadRefinement" setup = [Setup] begin