⚾ Refactor meshbuilder and improve clarity in 2D and 3D mesh generation.

Author: ZenanH

src/utils.jl

@@ -43,7 +43,7 @@ function get_polygon(points::AbstractMatrix; ratio::Real=0.1, allow_holes::Bool=
 end
 
 """
-    meshbuilder(x::T, y::T; ϵ::String="FP64") where T <: AbstractRange
+    meshbuilder(x::AbstractRange, y::AbstractRange)
 
 Description:
 ---
@@ -54,18 +54,26 @@ Example:
 ```julia
 x = 0:0.1:1
 y = 0:0.1:1
-mesh = meshbuilder(x, y; ϵ="FP32") # returns a 2xN array of points
+mesh = meshbuilder(x, y) # returns a 2xN array of points
 ``` 
 """
-function meshbuilder(x::T, y::T; ϵ::String="FP64") where T <: AbstractRange
-    x_tmp = repeat(x', length(y), 1) |> vec
-    y_tmp = repeat(y , 1, length(x)) |> vec
-    T1 = ϵ == "FP32" ? Float32 : Float64
-    return Array{T1}(hcat(x_tmp, y_tmp)')
+function meshbuilder(x::AbstractRange, y::AbstractRange)::Array{Float64, 2}
+    nx, ny = length(x), length(y)
+    result = Array{Float64, 2}(undef, 2, nx * ny)
+    idx = 1
+    @inbounds for i = 1:nx
+        xi = x[i]
+        @simd for j = 1:ny
+            result[1, idx] = xi
+            result[2, idx] = y[j]
+            idx += 1
+        end
+    end
+    return result
 end
 
 """
-    meshbuilder(x::T, y::T, z::T; ϵ::String="FP64") where T <: AbstractRange
+    meshbuilder(x::AbstractRange, y::AbstractRange, z::AbstractRange)
 
 Description:
 ---
@@ -77,24 +85,25 @@ Example:
 x = 0:0.1:1
 y = 0:0.1:1
 z = 0:0.1:1
-mesh = meshbuilder(x, y, z; ϵ="FP32") # returns a 3xN array of points
+mesh = meshbuilder(x, y, z) # returns a 3xN array of points
 """
-function meshbuilder(x::T, y::T, z::T; ϵ::String="FP64") where T <: AbstractRange
-    vx      = x |> collect
-    vy      = y |> collect
-    vz      = z |> collect
-    m, n, o = length(vy), length(vx), length(vz)
-    vx      = reshape(vx, 1, n, 1)
-    vy      = reshape(vy, m, 1, 1)
-    vz      = reshape(vz, 1, 1, o)
-    om      = ones(Int, m)
-    on      = ones(Int, n)
-    oo      = ones(Int, o)
-    x_tmp   = vec(vx[om, :, oo])
-    y_tmp   = vec(vy[:, on, oo])
-    z_tmp   = vec(vz[om, on, :])
-    T1 = ϵ == "FP32" ? Float32 : Float64
-    return Array{T1}(hcat(x_tmp, y_tmp, z_tmp)')
+function meshbuilder(x::AbstractRange, y::AbstractRange, z::AbstractRange)::Array{Float64, 2}
+    nx, ny, nz = length(x), length(y), length(z)
+    result = Array{Float64, 2}(undef, 3, nx * ny * nz)
+    idx = 1
+    @inbounds for k = 1:nz
+        zk = z[k]
+        for j = 1:nx
+            xj = x[j]
+            @simd for i = 1:ny
+                result[1, idx] = xj
+                result[2, idx] = y[i]
+                result[3, idx] = zk
+                idx += 1
+            end
+        end
+    end
+    return result
 end
 
 """
@@ -149,8 +158,8 @@ If `xy` is false, sort the points by the z-, y-, and then x-coordinates (3D).
 function sortpts(pts::AbstractMatrix; xy::Bool=true)
     if size(pts, 1) == 2 || xy
         perm = sortperm(axes(pts, 2), by = i -> (pts[2, i], pts[1, i]))
-    elseif size(pts, 1) == 3
-        perm = sortperm(axes(pts, 2), by = i -> (pts[3, i], pts[2, i], pts[1, i]))
+    elseif size(pts, 1) == 3🍻
+        perm = sortperm(axes(pts, 2), by = i -> (pts[2, i], pts[1, i], pts[3, i]))
     else
         throw(ArgumentError("The input points should have 2 or 3 rows (2/3D)"))
     end