Improvements related to hemisphere and disc meshing functions

examples/demo_geosphere.jl

@@ -16,18 +16,29 @@ n2 = 1 # Number of refinement iterations
 F2,V2 = geosphere(n2,r)
 
 n3 = 3 # Number of refinement iterations
-F3,V3 = geosphere(n3,r)
+method3 = :linear
+F3,V3 = geosphere(n3,r; method=method3)
+
+n4 = 3 # Number of refinement iterations
+method4 = :Loop
+F4,V4 = geosphere(n4,r; method=method4)
+
 
 #Visualize mesh
+lineWidth = 1
+
 fig = Figure(size=(1600,800))
 
-ax1=Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "A geodesic sphere n=0")
-hp1=poly!(ax1,V1,F1, strokewidth=2,color=:white, shading = FastShading)
+ax1=Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "A geodesic sphere n=$n1")
+hp1=poly!(ax1,V1,F1, strokewidth=lineWidth,color=:white, shading = FastShading)
+
+ax2=Axis3(fig[1, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "A geodesic sphere n=$n2")
+hp2=poly!(ax2,V2,F2, strokewidth=lineWidth,color=:white, shading = FastShading)
 
-ax2=Axis3(fig[1, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "A geodesic sphere n=1")
-hp2=poly!(ax2,V2,F2, strokewidth=2,color=:white, shading = FastShading)
+ax3=Axis3(fig[2, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "A geodesic sphere n=$n3, method=$method3")
+hp3=poly!(ax3,V3,F3, strokewidth=lineWidth,color=:white, shading = FastShading)
 
-ax3=Axis3(fig[1, 3], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "A geodesic sphere n=3")
-hp3=poly!(ax3,V3,F3, strokewidth=2,color=:white, shading = FastShading)
+ax4=Axis3(fig[2, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "A geodesic sphere n=$n4, method=$method4")
+hp4=poly!(ax4,V4,F4, strokewidth=lineWidth,color=:white, shading = FastShading)
 
 fig

examples/demo_hemisphere.jl

@@ -6,36 +6,44 @@ using Statistics
 using LinearAlgebra
 
 
-# Demo for building a  hemissphere from a isosaheddron .
+# Demo for building a  hemissphere from an isosaheddron.
 
 # Example geometry for a sphere that is cut so some edges are boundary edges
 n = 1# Number of refinement steps of the geodesic sphere
-r = 1.0 # Sphere radius
+r = 1.0 # hemisphere radius
 
-r = 1.0
-F1,V1 = hemisphere(1,r; face_type=:tri)
-F2,V2 = hemisphere(2,r; face_type=:tri)
-F3,V3 = hemisphere(3,r; face_type=:tri)
+F1,V1 = hemisphere(0,r; face_type=:tri,closed=true)
+F2,V2 = hemisphere(1,r; face_type=:tri)
+F3,V3 = hemisphere(2,r; face_type=:tri,closed=true)
+F4,V4 = hemisphere(3,r; face_type=:tri)
 
-F4,V4 = hemisphere(1,r; face_type=:quad)
-F5,V5 = hemisphere(2,r; face_type=:quad)
-F6,V6 = hemisphere(3,r; face_type=:quad)
+F5,V5 = hemisphere(0,r; face_type=:quad,closed=true)
+F6,V6 = hemisphere(1,r; face_type=:quad)
+F7,V7 = hemisphere(2,r; face_type=:quad,closed=true)
+F8,V8 = hemisphere(2,r; face_type=:quad)
 
 # Visualization
 fig = Figure(size=(1200,800))
 
 ax1 = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "face_type=:tri, n=0")
 hp1 = poly!(ax1,GeometryBasics.Mesh(V1,F1), strokewidth=1,color=:white, strokecolor=:black, shading = FastShading, transparency=false)
-ax2 = Axis3(fig[1, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "face_type=:tri, n=2")
-hp2 = poly!(ax2,GeometryBasics.Mesh(V2,F2), strokewidth=1,color=:white, strokecolor=:black, shading = FastShading, transparency=false)
-ax3 = Axis3(fig[1, 3], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "face_type=:tri, n=3")
-hp4 = poly!(ax3,GeometryBasics.Mesh(V3,F3), strokewidth=1,color=:white, strokecolor=:black, shading = FastShading, transparency=false)
+# normalplot(ax1,F1,V1)
 
-ax4 = Axis3(fig[2, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "face_type=:quad, n=0")
+ax2 = Axis3(fig[1, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "face_type=:tri, n=1")
+hp2 = poly!(ax2,GeometryBasics.Mesh(V2,F2), strokewidth=1,color=:white, strokecolor=:black, shading = FastShading, transparency=false)
+ax3 = Axis3(fig[1, 3], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "face_type=:tri, n=2")
+hp3 = poly!(ax3,GeometryBasics.Mesh(V3,F3), strokewidth=1,color=:white, strokecolor=:black, shading = FastShading, transparency=false)
+ax4 = Axis3(fig[1, 4], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "face_type=:tri, n=3")
 hp4 = poly!(ax4,GeometryBasics.Mesh(V4,F4), strokewidth=1,color=:white, strokecolor=:black, shading = FastShading, transparency=false)
-ax5 = Axis3(fig[2, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "face_type=:quad, n=2")
+
+ax5 = Axis3(fig[2, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "face_type=:quad, n=0")
 hp5 = poly!(ax5,GeometryBasics.Mesh(V5,F5), strokewidth=1,color=:white, strokecolor=:black, shading = FastShading, transparency=false)
-ax6 = Axis3(fig[2, 3], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "face_type=:quad, n=3")
+# normalplot(ax5,F5,V5)
+ax6 = Axis3(fig[2, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "face_type=:quad, n=1")
 hp6 = poly!(ax6,GeometryBasics.Mesh(V6,F6), strokewidth=1,color=:white, strokecolor=:black, shading = FastShading, transparency=false)
+ax7 = Axis3(fig[2, 3], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "face_type=:quad, n=2")
+hp7 = poly!(ax7,GeometryBasics.Mesh(V7,F7), strokewidth=1,color=:white, strokecolor=:black, shading = FastShading, transparency=false)
+ax8 = Axis3(fig[2, 4], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "face_type=:quad, n=3")
+hp8 = poly!(ax8,GeometryBasics.Mesh(V8,F8), strokewidth=1,color=:white, strokecolor=:black, shading = FastShading, transparency=false)
 
 fig

examples/demo_perlin.jl

@@ -0,0 +1,130 @@
+using GLMakie
+using Random
+
+function randangle(siz)
+    A = Matrix{Float64}(undef,siz)    
+    for i in eachindex(A)
+        A[i] = rand()*pi*rand((-1,1))
+    end
+    return A
+end
+
+Random.seed!(1)
+
+# Define grid
+size_grid = (25,25) # grid size
+sampleFactor = 35 # Pixel sample factor wrt grid
+pixelSize = 1/sampleFactor # Pixel size assuming grid has unit steps
+
+# Create grid vectors 
+A = randangle(size_grid) # Random angles
+
+# for w = 0:1:19
+# A .+= (2*pi/20)
+
+# A = pi/2 .* ones(size_grid)
+# for i in eachindex(A)
+#     A[i] *= rand((-1,1))
+#     # println(i)
+# end
+# A[7]*=-1
+
+Ux = cos.(A) # Unit vector x-component
+Uy = sin.(A) # Unit vector y-component
+
+# Initialise image
+size_image = (size_grid .- 1) .* sampleFactor # image size
+M = Matrix{Float64}(undef,size_image) # Start as undef Float64
+
+# Pre-compute grid cell quantities
+xy = range(0+pixelSize/2,1-pixelSize/2,sampleFactor) # x or y coordinates within a grid cell
+
+X = [x for x in 0:sampleFactor:(size_grid[2]-1)*sampleFactor, j in 0:sampleFactor:(size_grid[1]-1)*sampleFactor]'
+Y = [y for i in 0:sampleFactor:(size_grid[2]-1)*sampleFactor, y in 0:sampleFactor:(size_grid[1]-1)*sampleFactor]'
+
+# ig = ceil(Int,i./sampleFactor) # Grid row index
+                # ip = mod1(i,sampleFactor) # Cell row index
+
+                # jg = ceil(Int,j./sampleFactor) # Grid column index
+                # jp = mod1(j,sampleFactor) # Cell column index
+
+# ceil.(collect(1:1:25)./sampleFactor)
+
+function fade_perlin(a,b,t)
+    # q = 3.0*t^3 - 2.0*t^3 # Smoothstep
+    q = 0.5 .- 0.5.*cos.(t*pi) # Cosine
+    # q = t # Linear 
+
+    # Fade function q = 6t⁵-15t⁴+10t³    
+    # q = t^3 * (t * (6.0 * t - 15.0) + 10.0) # Perlin fade function
+
+    return (1.0-q)*a + q*b
+end
+
+xc = [0,1,1,0]
+yc = [0,0,1,1]
+for ip in 1:sampleFactor # For each cell row
+    for jp in 1:sampleFactor # For each cell column
+        for ig in 1:size_grid[1]-1 # For each grid row    
+            for jg in 1:size_grid[2]-1 # For each grid column
+                i = (ig-1)*sampleFactor + ip # Pixel row index
+                j = (jg-1)*sampleFactor + jp # Pixel column index
+
+                # Current pixel cell coordinates
+                px = xy[jp]
+                py = xy[ip]
+
+                # Offset vector components
+                xc1 = px # -xc[1] Offset vector 1 x
+                xc2 = px-xc[2] # Offset vector 2 x
+                xc3 = px-xc[3] # Offset vector 3 x
+                xc4 = px-xc[4] # Offset vector 4 x
+
+                yc1 = py # -yc[2] Offset vector 1 y                
+                yc2 = py-yc[2] # Offset vector 2 y               
+                yc3 = py-yc[3] # Offset vector 3 y                
+                yc4 = py-yc[4] # Offset vector 4 y
+
+                u1x = Ux[ig  ,jg]
+                u2x = Ux[ig  ,jg+1]
+                u3x = Ux[ig+1,jg+1]
+                u4x = Ux[ig+1,jg]
+
+                u1y = Uy[ig  ,jg]
+                u2y = Uy[ig  ,jg+1]
+                u3y = Uy[ig+1,jg+1]
+                u4y = Uy[ig+1,jg]
+
+                d1 = xc1.*u1x + yc1.*u1y
+                d2 = xc2.*u2x + yc2.*u2y
+                d3 = xc3.*u3x + yc3.*u3y
+                d4 = xc4.*u4x + yc4.*u4y
+
+                # Interpolation 
+                # Fade function 6t⁵-15t⁴+10t³
+                d12 = fade_perlin(d1,d2,px)
+                d34 = fade_perlin(d4,d3,px)
+                d = fade_perlin(d12,d34,py)
+
+                M[i,j] = d
+
+            end
+        end
+    end
+end
+
+fig = Figure(size=(1500,1500))
+# ax1 = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Edge angles")
+ax1 = Axis(fig[1, 1], aspect = DataAspect(), title = "Perlin noise",limits=(-sampleFactor,size_image[2]+sampleFactor,-sampleFactor,size_image[1]+sampleFactor) )
+hm = image!(ax1, M',interpolate=false,colormap = Makie.Reverse(:Spectral),colorrange=(-0.5,0.5)) #
+
+arrows!(ax1,reduce(vcat,X), reduce(vcat,Y), reduce(vcat,Ux), reduce(vcat,Uy), arrowsize = 10, lengthscale = sampleFactor/3, color = :black, linewidth=1)
+
+Colorbar(fig[1, 2], hm)
+
+fig
+
+# save("/home/kevin/DATA/Julia/Comodo.jl/assets/img/perlin_noise_c"*string(w)*".jpg",fig,px_per_unit = 1)
+# save("/home/kevin/DATA/Julia/Comodo.jl/assets/img/perlin_noise_large.jpg",fig,px_per_unit = 2)
+
+# end

examples/demo_quaddisc.jl

@@ -0,0 +1,48 @@
+using Comodo
+using GLMakie
+using GeometryBasics
+using LinearAlgebra
+
+#=
+This is a demonstration of the capabilities of the `quaddisc` function which 
+generates the faces `F` and vertices `V` for a quadrangulated disc (circle).
+=#
+
+# Define input parameters
+
+
+r = 1.0 # Radius
+
+n1 = 0
+method1 = :linear
+F1,V1 = quaddisc(r,n1; method=method1)
+
+n2 = 1
+method2 = :Catmull_Clark
+F2,V2 = quaddisc(r,n2; method=method2)
+
+n3 = 2
+method3 = :Catmull_Clark
+F3,V3 = quaddisc(r,n3; method=method3)
+
+n4 = 3
+method4 = :Catmull_Clark
+F4,V4 = quaddisc(r,n4; method=method4, orientation=:down)
+
+# Visualization
+fig = Figure(size=(800,800))
+
+ax1 = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Quandrangulated disc, n=$n1, method=$method1")
+hp1 = poly!(ax1,GeometryBasics.Mesh(V1,F1), strokewidth=1,color=:white,shading=FastShading,transparency=false)
+
+ax2 = Axis3(fig[1, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Quandrangulated disc, n=$n2, method=$method2")
+hp2 = poly!(ax2,GeometryBasics.Mesh(V2,F2), strokewidth=1,color=:white,shading=FastShading,transparency=false)
+
+ax3 = Axis3(fig[2, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Quandrangulated disc, n=$n3, method=$method3")
+hp3 = poly!(ax3,GeometryBasics.Mesh(V3,F3), strokewidth=1,color=:white,shading=FastShading,transparency=false)
+
+ax4 = Axis3(fig[2, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Quandrangulated disc, n=$n4, method=$method4")
+hp4 = poly!(ax4,GeometryBasics.Mesh(V4,F4), strokewidth=1,color=:white,shading=FastShading,transparency=false)
+
+
+fig

examples/demo_quadplate.jl

@@ -2,18 +2,15 @@ using Comodo
 using GLMakie
 using GeometryBasics
 
-plateDim = [20.0,24.0]
-plateElem = [11,16]
-
-F,V = quadplate(plateDim,plateElem)
-N = facenormal(F,V)
-
+plateDim1 = [20.0,24.0]
+plateElem1 = [11,16]
+orientation1 = :up
+F1,V1 = quadplate(plateDim1,plateElem1; orientation=orientation1)
 
 ## Visualize mesh
-fig = Figure(size=(1000,1000))
+fig = Figure(size=(1200,800))
+
+ax1 = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Quadrilateral mesh plate")
+hp2 = poly!(ax1,GeometryBasics.Mesh(V1,F1), strokewidth=3,color=:white, shading = FastShading)
 
-ax1=Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "A quadrilateral mesh of a plate",azimuth=-pi/2,elevation=pi/2)
-hp2=poly!(ax1,GeometryBasics.Mesh(V,F), strokewidth=3,color=:white, shading = FastShading)
-# hp2 = scatter!(ax1, V,markersize=15,color=:orange)
-# normalplot(ax1,F,V)
 fig

examples/demo_subquad.jl

@@ -51,11 +51,11 @@ ax1 = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z",
 wireframe!(ax1,M, linewidth=lineWidth,color=:red, overdraw=false)
 poly!(ax1,GeometryBasics.Mesh(Vn1,Fn1), strokewidth=strokewidth1,color=:white,shading=FastShading,transparency=false)
 
-ax2 = Axis3(fig[1, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Linearly, n=2")
+ax2 = Axis3(fig[1, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Linear, n=2")
 wireframe!(ax2,M, linewidth=lineWidth,color=:red, overdraw=false)
 poly!(ax2,GeometryBasics.Mesh(Vn2,Fn2), strokewidth=strokewidth1,color=:white,shading=FastShading,transparency=false)
 
-ax3 = Axis3(fig[1, 3], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Linearly, n=3")
+ax3 = Axis3(fig[1, 3], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Linear, n=3")
 hp1 = wireframe!(ax3,M, linewidth=lineWidth,color=:red, overdraw=false)
 hp2 = poly!(ax3,GeometryBasics.Mesh(Vn3,Fn3), strokewidth=strokewidth1,color=:white,shading=FastShading,transparency=false)
 Legend(fig[1, 4],[hp1,hp2],["Initial","Refined"])

examples/demo_tridisc.jl

@@ -9,17 +9,40 @@ generates the faces `F` and vertices `V` for a triangulated disc (circle).
 =#
 
 # Define input parameters
-r = 2.0 # Radius
-n = 3 # Number of refinement iterations
+r = 1.0 # Radius
 
 # Create disc mesh 
-F,V = tridisc(r,n)
+n1 = 0
+F1,V1 = tridisc(r,n1)
 
-## Visualization
+n2 = 2
+F2,V2 = tridisc(r,n2)
+
+n3 = 3
+ngon3 = 5
+method3 = :linear
+F3,V3 = tridisc(r,n3; ngon=ngon3, method=method3)
+
+n4 = 3
+ngon4 = 6
+method4 = :Loop
+F4,V4 = tridisc(r,n4; ngon=ngon4, method=method4)
+
+# Visualization
 fig = Figure(size=(800,800))
 
-ax1 = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "A triangulated mesh of a circular domain (disc)")
-# scatter!(ax1,V,markersize=10,color = :black)
-hp2 = poly!(ax1,GeometryBasics.Mesh(V,F), strokewidth=1,color=:white,shading=FastShading,transparency=false)
-# normalplot(ax1,F,V)
+ax1 = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Triangulated disc, n=$n1")
+hp1 = poly!(ax1,GeometryBasics.Mesh(V1,F1), strokewidth=1,color=:white,shading=FastShading,transparency=false)
+# normalplot(ax1,F1,V1)
+
+ax2 = Axis3(fig[1, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Triangulated disc, n=$n2")
+hp2 = poly!(ax2,GeometryBasics.Mesh(V2,F2), strokewidth=1,color=:white,shading=FastShading,transparency=false)
+
+ax3 = Axis3(fig[2, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Triangulated disc, n=$n3, ngon=$ngon3, method=$method3")
+hp3 = poly!(ax3,GeometryBasics.Mesh(V3,F3), strokewidth=1,color=:white,shading=FastShading,transparency=false)
+
+ax4 = Axis3(fig[2, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Triangulated disc, n=$n4, ngon=$ngon4, method=$method4")
+hp4 = poly!(ax4,GeometryBasics.Mesh(V4,F4), strokewidth=1,color=:white,shading=FastShading,transparency=false)
+
+
 fig

examples/demo_triplate.jl

@@ -0,0 +1,18 @@
+using Comodo
+using GLMakie
+using GeometryBasics
+
+plateDim1 = [20.0,24.0]
+pointSpacing1 = 2.0
+
+orientation1 = :up
+F1,V1 = triplate(plateDim1,pointSpacing1; orientation=orientation1)
+
+## Visualization
+linewidth = 1
+fig = Figure(size=(1200,800))
+
+ax1 = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Triangulated mesh plate")
+hp2 = poly!(ax1,GeometryBasics.Mesh(V1,F1), strokewidth=linewidth,color=:white, shading = FastShading)
+# normalplot(ax1,F1,V1)
+fig