Merge pull request #36 from chethana-rao/main

Hemisphere

examples/demo_hemisphere_mesh.jl

@@ -0,0 +1,32 @@
+using Comodo
+using GLMakie
+using GeometryBasics
+using Rotations
+using Statistics
+using LinearAlgebra
+
+
+# Demo for building a  hemissphere from a 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
+F0,V0 = hemisphere(0,r)
+F1,V1 = hemisphere(1,r)
+F2,V2 = hemisphere(2,r)
+F3,V3 = hemisphere(3,r)
+
+# Visualization
+fig = Figure(size=(1200,800))
+ax1 = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "n=0")
+hp1 = poly!(ax1,GeometryBasics.Mesh(V0,F0), strokewidth=1,color=:white, strokecolor=:black, shading = FastShading, transparency=false)
+
+ax2 = Axis3(fig[1, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "n=1")
+hp2 = poly!(ax2,GeometryBasics.Mesh(V1,F1), strokewidth=1,color=:white, strokecolor=:black, shading = FastShading, transparency=false)
+ax3 = Axis3(fig[2, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "n=2")
+hp3 = poly!(ax3,GeometryBasics.Mesh(V2,F2), strokewidth=1,color=:white, strokecolor=:black, shading = FastShading, transparency=false)
+ax4 = Axis3(fig[2, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "n=3")
+hp4 = poly!(ax4,GeometryBasics.Mesh(V3,F3), strokewidth=1,color=:white, strokecolor=:black, shading = FastShading, transparency=false)
+fig

src/Comodo.jl

@@ -31,7 +31,7 @@ export gridpoints, gridpoints_equilateral
 export interp_biharmonic_spline, interp_biharmonic, nbezier
 export ind2sub, sub2ind, meshedges, edgelengths 
 export gunique, unique_simplices, unique_dict, occursonce
-export subtri, subquad, geosphere, quad2tri 
+export subtri, subquad, geosphere,hemisphere,quad2tri 
 export icosahedron, tetrahedron, cube, dodecahedron, octahedron, platonicsolid
 export tofaces, topoints, togeometrybasics_mesh
 export quadplate, quadsphere, smoothmesh_laplacian, smoothmesh_hc

src/functions.jl

@@ -1847,6 +1847,11 @@ function subquad(F::Vector{NgonFace{4,TF}},V::Vector{Point{ND,TV}},n::Int; metho
 end
 
 
+function pushtoradius_(v,r=1.0)    
+    return v.*(r/sqrt(sum(v.^2))) # Push to sphere       
+end
+
+
 """
     geosphere(n::Int,r::T) where T <: Real
 
@@ -1865,17 +1870,74 @@ function geosphere(n::Int,r::T) where T <: Real
     M = platonicsolid(4,r)
     V = coordinates(M)
     F = faces(M)
-    for _ = 1:n
+    for _ = 1:n        
+        numPointsNow = length(V)
         F,V = subtri(F,V,1)
-        for q in eachindex(V)
-            v = V[q]
-            rn = sqrt(sum(v.^2))
-            V[q] = v .* (r/rn)
+        @inbounds for i in numPointsNow+1:length(V)
+            V[i] = pushtoradius_(V[i],r) # Overwrite points
         end
     end
     return F,V
 end
 
+"""
+    geosphere(n::Int,r::T) where T <: Real
+
+Returns a geodesic sphere triangulation
+
+# Description
+
+This function returns a geodesic hemispherephere triangulation based on the number of
+refinement iterations `n` and the radius `r`. Geodesic spheres (aka Buckminster-Fuller
+ spheres) are triangulations of a sphere that have near uniform edge lenghts. 
+The algorithm starts with a regular icosahedron that is adjusted to generate a half icosahedron. 
+Next this icosahedron is refined  `n` times, while nodes are pushed to a sphere surface with radius `r` at each
+iteration. 
+"""
+function hemisphere(n::Int,r::T) where T <: Real
+    if n == 0 #creates a half a isosaheddron
+        F,V = geosphere(0,r) # Creating the faces and vertices of a full sphere
+
+        VC = simplexcenter(F,V) # Finding triangle centre coordiantes
+        F = [F[i] for i in findall(map(v-> v[3]>=0,VC))] # Remove some faces using z of central coordinates
+        F,V = remove_unused_vertices(F,V) # Cleanup/remove unused vertices after faces were removed
+
+        for (i,v) in enumerate(V)
+            if v[3]<0
+                v = Point3([v[1],v[2],0.0]) # Make z zero
+                V[i] = pushtoradius_(v,r) # Overwrite point        
+            end
+        end
+    elseif n>0
+        F,V = geosphere(1,r) # Creating the faces and vertices of a full sphere
+
+        ϕ = Base.MathConstants.golden # (1.0+sqrt(5.0))/2.0, Golden ratio
+        s = r/sqrt(ϕ + 2.0)
+        t = ϕ*s
+        α = atan(t/s) # angle needed for roating the sphere to 90 deg
+        Q = RotXYZ(0.0,α,0.0)
+        V = [GeometryBasics.Point{3, Float64}(Q*v) for v in V] 
+
+        VC = simplexcenter(F,V) # Finding triangle centre coordiantes
+        F = [F[i] for i in findall(map(v-> v[3]>=0,VC))] # Remove some faces at central coordinates z=0
+        F,V = remove_unused_vertices(F,V) # Cleanup/remove unused vertices after faces were removed
+
+        if n>1
+            for _ = 1:n-1
+                numPointsNow = length(V)
+                F,V = subtri(F,V,1)
+                @inbounds for i in numPointsNow+1:length(V)
+                    V[i] = pushtoradius_(V[i],r) # Overwrite points
+                end
+            end
+        end
+
+    else
+        # error message here
+    end
+    return F,V
+end
+
 """
     hexbox(boxDim,boxEl)