Refactor RegularDiscretization

src/discretization/regular.jl

@@ -93,6 +93,43 @@ function appendtopo(::Sphere{3}, tg)
   SimpleTopology([middle; north; south])
 end
 
+function appendtopo(::CylinderSurface, tg)
+  sz = size(tg)
+  ip = isperiodic(tg)
+  np = @. sz + !ip
+  nx, ny = np
+
+  # connect quadrangles in the middle
+  middle = collect(elements(tg))
+
+  # connect south pole with triangles
+  south = map(1:(nx - 1)) do i
+    u = nx * ny + 1
+    v = i + 1
+    w = i
+    connect((u, v, w))
+  end
+  u = nx * ny + 1
+  v = 1
+  w = nx
+  push!(south, connect((u, v, w)))
+
+  # connect north pole with triangles
+  offset = nx * ny - nx
+  north = map(1:(nx - 1)) do i
+    u = nx * ny + 2
+    v = offset + i + 1
+    w = offset + i
+    connect((u, w, v))
+  end
+  u = nx * ny + 2
+  v = nx * ny - nx + 1
+  w = nx * ny
+  push!(north, connect((u, w, v)))
+
+  SimpleTopology([middle; north; south])
+end
+
 # --------------
 # SPECIAL CASES
 # --------------
@@ -134,55 +171,6 @@ function _rball(ball, method::RegularDiscretization)
   SimpleMesh(points, connec)
 end
 
-function discretize(cylsurf::CylinderSurface, method::RegularDiscretization)
-  nx, ny = fitdims(method.sizes, paramdim(cylsurf))
-
-  # sample points regularly
-  sampler = RegularSampling(nx, ny)
-  points = collect(sample(cylsurf, sampler))
-
-  # connect regular samples with quadrangles
-  topo = GridTopology((nx - 1, ny - 1))
-  middle = collect(elements(topo))
-  for j in 1:(ny - 1)
-    u = (j) * nx
-    v = (j - 1) * nx + 1
-    w = (j) * nx + 1
-    z = (j + 1) * nx
-    quad = connect((u, v, w, z))
-    push!(middle, quad)
-  end
-
-  # connect south pole with triangles
-  south = map(1:(nx - 1)) do i
-    u = nx * ny + 1
-    v = i + 1
-    w = i
-    connect((u, v, w))
-  end
-  u = nx * ny + 1
-  v = 1
-  w = nx
-  push!(south, connect((u, v, w)))
-
-  # connect north pole with triangles
-  offset = nx * ny - nx
-  north = map(1:(nx - 1)) do i
-    u = nx * ny + 2
-    v = offset + i + 1
-    w = offset + i
-    connect((u, w, v))
-  end
-  u = nx * ny + 2
-  v = nx * ny - nx + 1
-  w = nx * ny
-  push!(north, connect((u, w, v)))
-
-  connec = [middle; south; north]
-
-  SimpleMesh(points, connec)
-end
-
 function discretize(consurf::ConeSurface, method::RegularDiscretization)
   nx, ny = fitdims(method.sizes, paramdim(consurf))
 

test/discretization.jl

@@ -50,8 +50,8 @@
 
     cylsurf = CylinderSurface(Plane(P3(0, 0, 0), V3(0, 0, 1)), Plane(P3(1, 1, 1), V3(0, 0, 1)), T(1))
     mesh = discretize(cylsurf, RegularDiscretization(10))
-    @test nvertices(mesh) == 10 * 10 + 2
-    @test nelements(mesh) == 10 * (10 - 1) + 2 * 10
+    @test nvertices(mesh) == 10 * 11 + 2
+    @test nelements(mesh) == 10 * 10 + 2 * 10
     @test eltype(mesh) <: Ngon
     @test nvertices.(mesh) ⊆ [3, 4]
 
@@ -369,7 +369,7 @@
     mesh = discretize(cylsurf)
     @test !(eltype(mesh) <: Triangle)
     @test !(eltype(mesh) <: Quadrangle)
-    @test nelements(mesh) == 150
+    @test nelements(mesh) == 200
 
     grid = CartesianGrid(10)
     @test discretize(grid) == grid