Faster Triangulated Grid Generation

src/Meshes.jl

@@ -69,7 +69,7 @@ loadmesh_icosphere_helper(obj_file_name) = EmbeddedDeltaSet2D(
 
 
 # This function was once the gridmeshes.jl file from Decapodes.jl.
-"""    function triangulated_grid(max_x, max_y, dx, dy, point_type, compress=false)
+"""    function triangulated_grid(max_x, max_y, dx, dy, point_type, compress=true)
 
 Triangulate the rectangle [0,max_x] x [0,max_y] by approximately equilateral
 triangles of width `dx` and height `dy`.
@@ -80,53 +80,60 @@ otherwise, keep `dx` as is.
 function triangulated_grid(max_x, max_y, dx, dy, point_type, compress=true)
   s = EmbeddedDeltaSet2D{Bool, point_type}()
 
-  # Place equally-spaced points in a max_x by max_y rectangle.
-  coords = point_type == Point3D ? map(x -> point_type(x..., 0), Iterators.product(0:dx:max_x, 0:dy:max_y)) : map(x -> point_type(x...), Iterators.product(0:dx:max_x, 0:dy:max_y))
-  # Perturb every other row right by half a dx.
-  coords[:, 2:2:end] = map(coords[:, 2:2:end]) do row
-    if point_type == Point3D
-      row .+ [dx/2, 0, 0]
-    else
-      row .+ [dx/2, 0]
+  scale = max_x/(max_x+dx/2)
+  shift = dx/2
+
+  nx = length(0:dx:max_x)
+  ny = length(0:dy:max_y)
+  add_vertices!(s, nx * ny)
+
+  for (y_idx, raw_y) in enumerate(0:dy:max_y)
+    for (x_idx, raw_x) in enumerate(0:dx:max_x)
+      x_coord = raw_x + shift * iseven(y_idx)
+      if compress
+        x_coord *= scale
+      end
+
+      i = x_idx + nx * (y_idx - 1)
+      s[i, :point] = if point_type <: Point3
+        point_type(x_coord, raw_y, 0)
+      else
+        point_type(x_coord, raw_y)
+      end
+
     end
   end
 
-  if compress
-    # The perturbation moved the right-most points past max_x, so compress along x.
-    map!(coords, coords) do coord
-      if point_type == Point3D
-        diagm([max_x/(max_x+dx/2), 1, 1]) * coord
+  for y in 1:ny-1, x in 1:nx-1
+      i = x + nx * (y - 1)
+      if iseven(y)
+        glue_triangle!(s, i, i+nx, i+nx+1)
+        glue_triangle!(s, i, i+1, i+nx+1)
       else
-        diagm([max_x/(max_x+dx/2), 1]) * coord
-      end
+        glue_triangle!(s, i, i+1, i+nx)
+        glue_triangle!(s, i+1, i+nx, i+nx+1)
     end
   end
-
-  add_vertices!(s, length(coords), point = vec(coords))
 
-  nx = length(0:dx:max_x)
+  # Orient and return.
+  s[:edge_orientation] = false
 
-  # Matrix that stores indices of points.
-  idcs = reshape(eachindex(coords), size(coords))
-  # Only grab vertices that will be the bottom-left corner of a subdivided square.
-  idcs = idcs[begin:end-1, begin:end-1]
-
-  # Subdivide every other row along the opposite diagonal.
-  for i in idcs[:, begin+1:2:end]
-    glue_sorted_triangle!(s, i, i+nx, i+nx+1)
-    glue_sorted_triangle!(s, i, i+1, i+nx+1)
-  end
-  for i in idcs[:, begin:2:end]
-    glue_sorted_triangle!(s, i, i+1, i+nx)
-    glue_sorted_triangle!(s, i+1, i+nx, i+nx+1)
+  nxtri = 2 * (nx - 1)
+  nytri = ny - 1
+
+  for y in 1:nytri
+    tri_orient = !iseven(y)
+    for x in 1:2:nxtri
+      i = x + nxtri * (y - 1)
+      s[i, :tri_orientation] = tri_orient
+      s[i + 1, :tri_orientation] = !tri_orient
+    end
   end
 
-  # Orient and return.
-  s[:edge_orientation]=true
-  orient!(s)
   s
 end
 
+
 # This function was once the sphericalmeshes.jl file from Decapodes.jl.
 """
     makeSphere(minLat, maxLat, dLat, minLong, maxLong, dLong, radius)

test/Meshes.jl

@@ -3,6 +3,7 @@ module TestMeshes
 using Catlab, Catlab.CategoricalAlgebra
 using CombinatorialSpaces
 using Test
+using GeometryBasics: Point2, Point3
 
 magnitude = (sqrt ∘ (x -> foldl(+, x*x)))
 unit_radius = 1
@@ -94,6 +95,49 @@ thermo_icosphere5 = loadmesh(Icosphere(5, thermosphere_radius))
 @test ρ == thermosphere_radius
 @test euler_characteristic(thermo_icosphere5) == 2
 
+# Testing the Triangulated Grid
+function test_trigrid_size(s, max_x, max_y, dx, dy)
+  nx = length(0:dx:max_x)
+  ny = length(0:dy:max_y)
+
+  @test nparts(s, :V) == nx * ny
+  @test nparts(s, :Tri) == 2 * (nx - 1) * (ny - 1)
+end
+
+s = triangulated_grid(1, 1, 1, 1, Point2{Float64}, false)
+test_trigrid_size(s, 1, 1, 1, 1)
+@test s[:point] == [[0.0, 0.0], [1.0, 0.0], [0.5, 1.0], [1.5, 1.0]]
+@test s[:tri_orientation] == [true, false]
+@test orient!(s)
+
+s = triangulated_grid(1, 1, 1, 1, Point2{Float64}, true)
+test_trigrid_size(s, 1, 1, 1, 1)
+@test s[:point] == [[0.0, 0.0], [2/3, 0.0], [1/3, 1.0], [1.0, 1.0]]
+@test s[:tri_orientation] == [true, false]
+@test orient!(s)
+
+s = triangulated_grid(2, 2, 1, 1, Point2{Float64}, false)
+test_trigrid_size(s, 2, 2, 1, 1)
+@test s[:point] == [[0.0, 0.0], [1.0, 0.0], [2.0, 0.0],
+                    [0.5, 1.0], [1.5, 1.0], [2.5, 1.0],
+                    [0.0, 2.0], [1.0, 2.0], [2.0, 2.0]]
+@test s[:tri_orientation] == [true, false, true, false, false, true, false, true]
+@test orient!(s)
+
+s = triangulated_grid(1, 1, 0.49, 0.78, Point2{Float64}, false)
+test_trigrid_size(s, 1, 1, 0.49, 0.78)
+@test orient!(s)
+
+s = triangulated_grid(1.6, 3.76, 0.49, 0.78, Point2{Float64}, true)
+test_trigrid_size(s, 1.6, 3.76, 0.49, 0.78)
+@test orient!(s)
+
+lx = 25
+s = triangulated_grid(lx, 2, 1, 0.99, Point2{Float64}, true)
+test_trigrid_size(s, lx, 2, 1, 0.99)
+@test maximum(getindex.(s[:point], 1)) == lx
+
+
 # Tests for the SphericalMeshes
 ρ = 6371+90
 s, npi, spi = makeSphere(0, 180, 5, 0, 360, 5, ρ)

test/Operators.jl

@@ -67,20 +67,20 @@ subdivide_duals!(rect, Barycenter());
 flat_meshes = [tri_345()[2], tri_345_false()[2], right_scalene_unit_hypot()[2], grid_345()[2], tg, rect];
 
 @testset "Exterior Derivative" begin
-    for i in 0:0 
+    for i in 0:0
         for sd in dual_meshes_1D
             @test all(dec_differential(i, sd) .== d(i, sd))
         end
     end
-    for i in 0:1 
+    for i in 0:1
         for sd in dual_meshes_2D
             @test all(dec_differential(i, sd) .== d(i, sd))
         end
     end
 end
 
 @testset "Boundary" begin
-    for i in 1:1 
+    for i in 1:1
         for sd in dual_meshes_1D
             @test all(dec_boundary(i, sd) .== ∂(i, sd))
         end
@@ -94,7 +94,7 @@ end
 end
 
 @testset "Dual Derivative" begin
-    for i in 0:0 
+    for i in 0:0
         for sd in dual_meshes_1D
             @test all(dec_dual_derivative(i, sd) .== dual_derivative(i, sd))
         end
@@ -406,9 +406,9 @@ end
 
   X♯ = SVector{3,Float64}(3,3,0)
   mag_selfadv, mag_dp, mag_∂ₜu  = euler_equation_test(X♯, tg)
-  @test .60 < (count(mag_selfadv .< 1e-1) / length(mag_selfadv))
-  @test .60 < (count(mag_dp .< 1e-1) / length(mag_dp))
-  @test .60 < (count(mag_∂ₜu .< 1e-1) / length(mag_∂ₜu))
+  @test 97/162 <= (count(mag_selfadv .< 1e-1) / length(mag_selfadv))
+  @test 97/162 <= (count(mag_dp .< 1e-1) / length(mag_dp))
+  @test 97/162 <= (count(mag_∂ₜu .< 1e-1) / length(mag_∂ₜu))
 
   # u := ⋆xdx
   # ιᵤu = x²