Bug fixes and added tests for meshgroup and distmarch and related

examples/demo_distmarch.jl

@@ -6,7 +6,7 @@ using LinearAlgebra
 using ProgressMeter
 
 # Example geometry
-testCase = 7
+testCase = 3
 if testCase == 1
     s=1.0
     V=Vector{GeometryBasics.Point{3, Float64}}(undef,5)
@@ -18,6 +18,7 @@ if testCase == 1
     F[1 ] = TriangleFace{Int64}(1,2,3)
     # F,V=subtri(F,V,2)
 elseif testCase==2
+    # Bowtie mesh
     s=1.0
     V=Vector{GeometryBasics.Point{3, Float64}}(undef,5)
     V[1 ] = GeometryBasics.Point{3, Float64}( 0.0,    s, 0.0)
@@ -28,7 +29,7 @@ elseif testCase==2
 
     F = Vector{TriangleFace{Int64}}(undef,2)
     F[1 ] = TriangleFace{Int64}(1,2,3)
-    F[2 ] = TriangleFace{Int64}(3,4,5)
+    F[2 ] = TriangleFace{Int64}(5,4,3)
     # F,V=subtri(F,V,2)
 elseif testCase==3
     r = 1.0
@@ -64,7 +65,6 @@ elseif testCase==7 # Merged STL for single object
     F,V,_ = mergevertices(F,V)
 end
 
-
 z = [v[3] for v ∈ V]
 
 ind=[findmin(z)[2]]

examples/demo_quadplate.jl

@@ -2,7 +2,7 @@ using Comodo
 using GLMakie
 using GeometryBasics
 
-plateDim = [20,20]
+plateDim = [20.0,20.0]
 plateElem = [10,10]
 
 F,V = quadplate(plateDim,plateElem)

src/functions.jl

@@ -1403,40 +1403,48 @@ function con_edge_face(F,E_uni=nothing,indReverse=nothing)
 end
 
 function con_face_face(F,E_uni=nothing,indReverse=nothing,con_E2F=nothing,con_F2E=nothing)
-    if isnothing(E_uni)| isnothing(indReverse)
-        E = meshedges(F)
-        E_uni,_,indReverse = gunique(E; return_unique=true, return_index=true, return_inverse=true, sort_entries=true)    
-    end    
-    if isnothing(con_E2F) 
-        con_E2F = con_edge_face(F,E_uni)
-    end
-    if isnothing(con_F2E)
-        con_F2E = con_face_edge(F,E_uni,indReverse)                 
-    end
-    con_F2F = [Vector{Int64}() for _ ∈ 1:length(F)]
-    for i_f ∈ eachindex(F)
-        for i ∈ reduce(vcat,con_E2F[con_F2E[i_f]])    
-            if i!=i_f     
-                push!(con_F2F[i_f],i)
-            end 
+    if length(F)>1 # More than one face so compute connectivity
+        if isnothing(E_uni)| isnothing(indReverse)
+            E = meshedges(F)
+            E_uni,_,indReverse = gunique(E; return_unique=true, return_index=true, return_inverse=true, sort_entries=true)    
+        end    
+        if isnothing(con_E2F) 
+            con_E2F = con_edge_face(F,E_uni)
+        end
+        if isnothing(con_F2E)
+            con_F2E = con_face_edge(F,E_uni,indReverse)                 
+        end
+        con_F2F = [Vector{Int64}() for _ ∈ 1:length(F)]
+        for i_f ∈ eachindex(F)
+            for i ∈ reduce(vcat,con_E2F[con_F2E[i_f]])    
+                if i!=i_f     
+                    push!(con_F2F[i_f],i)
+                end 
+            end
         end
+        return con_F2F
+    else # Just one face, so return empty
+        return [Vector{Int64}() ]
     end
-    return con_F2F
 end
 
 function con_face_face_v(F,V=nothing,con_V2F=nothing)
-    if isnothing(con_V2F) 
-        con_V2F = con_vertex_face(F,V)  # VERTEX-FACE connectivity
-    end
-    con_F2F = [Vector{Int64}() for _ ∈ 1:length(F)]
-    for i_f ∈ eachindex(F)
-        for i ∈ reduce(vcat,con_V2F[F[i_f]])    
-            if i!=i_f     
-                push!(con_F2F[i_f],i)
-            end 
+    if length(F)>1 # More than one face so compute connectivity
+        if isnothing(con_V2F) 
+            con_V2F = con_vertex_face(F,V)  # VERTEX-FACE connectivity
+        end
+        con_F2F = [Vector{Int64}() for _ ∈ 1:length(F)]
+        for i_f ∈ eachindex(F)
+            for i ∈ reduce(vcat,con_V2F[F[i_f]])    
+                if i!=i_f     
+                    push!(con_F2F[i_f],i)
+                end 
+            end
         end
+        return con_F2F
+    else # Just one face, so return empty
+        return [Vector{Int64}() ]
     end
-    return con_F2F
 end
 
 function con_vertex_simplex(F,V=nothing)
@@ -1479,7 +1487,7 @@ end
 
 function con_vertex_vertex_f(F,V=nothing,con_V2F=nothing)
     if isnothing(V)
-        n = maximum(reduce(vcat,E))
+        n = maximum(reduce(vcat,F))
     else 
         n = length(V)
     end
@@ -1490,14 +1498,16 @@ function con_vertex_vertex_f(F,V=nothing,con_V2F=nothing)
 
     con_V2V = [Vector{Int64}() for _ ∈ 1:n]
     for i_v ∈ 1:n
-        for i ∈ unique(reduce(vcat,F[con_V2F[i_v]]))
-            if i_v!=i
-                push!(con_V2V[i_v],i)
+        if !isempty(con_V2F[i_v])
+            for i ∈ unique(reduce(vcat,F[con_V2F[i_v]]))
+                if i_v!=i
+                    push!(con_V2V[i_v],i)
+                end
             end
         end
     end
-
     return con_V2V
+
 end
 
 function con_vertex_vertex(E,V=nothing,con_V2E=nothing)
@@ -1512,9 +1522,11 @@ function con_vertex_vertex(E,V=nothing,con_V2E=nothing)
 
     con_V2V = [Vector{Int64}() for _ ∈ 1:n]
     for i_v ∈ 1:n
-        for i ∈ reduce(vcat,E[con_V2E[i_v]])
-            if i_v!=i
-                push!(con_V2V[i_v],i)
+        if !isempty(con_V2E[i_v])
+            for i ∈ reduce(vcat,E[con_V2E[i_v]])
+                if i_v!=i
+                    push!(con_V2V[i_v],i)
+                end
             end
         end
     end
@@ -2278,33 +2290,40 @@ function distmarch(F,V,indStart; d=nothing, dd=nothing, dist_tol=1e-3,con_V2V=no
     end
 
     # Set start distances to zero 
+    is_isolated =  isempty.(con_V2V)
+    d[is_isolated] .= NaN # Set isolated (non-connected) points to NaN
     d[indStart] .= 0.0
     l[indStart] .= 1:length(indStart)
 
-    c = false
-    ds = -1.0 # Set negative initially 
-
-    boolCheck = fill(true,length(V))
+    notGrowing = false
+    dist_sum_previous = -1.0 # Set negative initially 
+    count_inf_previous = length(d)-length(indStart) # number of Inf values currently
     while true                          
         for i ∈ eachindex(V) # For each point            
             for ii ∈ con_V2V[i] # Check umbrella neighbourhood
+                # Get closest point and distance from umbrella
                 minVal,minInd = findmin([d[ii],dd[sort([i,ii])]+d[i]])            
                 if minInd==2
-                    d[ii] = minVal                          
-                    l[ii] = l[i]
+                    d[ii] = minVal # Distance                          
+                    l[ii] = l[i] # Index
                 end
             end            
         end
-        if !any(isinf.(d)) # Start checking once all are no longer Inf
-            if c # If we were here before
-                if abs(sum(d)-ds)<dist_tol                                        
+        bool_inf = isinf.(d) # Booling to check number of points left at Inf
+        count_inf = count(bool_inf)
+        if count_inf == count_inf_previous #!any(isinf.(d)) # Start checking once all are no longer Inf
+            dist_sum = sum(d[.!is_isolated .&& .!bool_inf])
+            if notGrowing # If we were here before
+                if abs(dist_sum-dist_sum_previous)<dist_tol                                        
                     break                    
                 end
             end
-            c = true # Flip to denote we've been here           
-            ds = sum(d) # Now start computing the sum to check convergence
+            notGrowing = true # Flip to denote we've been here           
+            dist_sum_previous = dist_sum # Now start computing the sum to check convergence
         end
+        count_inf_previous = count_inf
     end
+    d[isinf.(d)] .= NaN # Change Inf to NaN
     return d,dd,l
 end
 

test/runtests.jl

@@ -1904,6 +1904,7 @@ end
 
 end
 
+
 @testset "meshgroup" verbose = true begin
     @testset "Single face" begin
         # Single triangle
@@ -1934,8 +1935,132 @@ end
         C = meshgroup(F)
         @test C == [1,2]
     end
+
+    @testset "Single group" begin
+        # Single tetrahedron
+        M = tetrahedron(1.0)
+        F = faces(M)
+        C = meshgroup(F)
+        @test C == ones(length(F))
+
+        # Single cube
+        M = cube(1.0)
+        F = faces(M)
+        C = meshgroup(F)
+        @test C == ones(length(F))
+    end
+
+    @testset "Two groups" begin
+        # Two tetrahedrons
+        M = tetrahedron(1.0)
+        F = faces(M)
+        V = coordinates(M)
+        n = length(F)
+        F2 = map(f-> f.+length(V),F)
+        V2 = map(v-> Point{3, Float64}(2.0+v[1],v[2],v[3]),V)
+        append!(F,F2)
+        append!(V,V2)
+        C = meshgroup(F)
+        @test C == repeat(1:2,inner=n)
+
+        # Two tetrahedrons
+        M = cube(1.0)
+        F = faces(M)
+        V = coordinates(M)
+        n = length(F)
+        F2 = map(f-> f.+length(V),F)
+        V2 = map(v-> Point{3, Float64}(2.0+v[1],v[2],v[3]),V)
+        append!(F,F2)
+        append!(V,V2)
+        C = meshgroup(F)
+        @test C == repeat(1:2,inner=n)
+    end
+end
+
+
+@testset "distmarch" verbose=true begin
+    eps_level = 0.001
+
+    @testset "Single face" begin
+        # Single triangle
+        F = TriangleFace{Int64}[[1,2,3]]
+        V = Point3{Float64}[[0.0,0.0,0.0],[1.0,0.0,0.0],[1.0,1.0,0.0]]
+        d,dd,l = distmarch(F,V,[1])
+        @test isapprox(d,[0.0,1.0,sqrt(2)],atol=eps_level)
+
+        # Single triangle, un-used nodes
+        F = TriangleFace{Int64}[[1,2,3]]
+        V = Point3{Float64}[[0.0,0.0,0.0],[1.0,0.0,0.0],[1.0,1.0,0.0],
+                            [1.0,0.0,0.0],[1.0,1.0,0.0]]
+        d,dd,l = distmarch(F,V,[1])
+        r = [0.0,1.0,sqrt(2),NaN,NaN]
+        b = .!isnan.(r)
+        @test isapprox(d[b],r[b],atol=eps_level) # Check reached entries
+        @test all(isnan.(d[.!b])) # Now check NaNs
+
+        # Single quad
+        F = QuadFace{Int64}[[1,2,3,4]]
+        V = Point3{Float64}[[0.0,0.0,0.0],[1.0,0.0,0.0],[1.0,1.0,0.0],[0.0,1.0,0.0]]
+        d,dd,l = distmarch(F,V,[1])
+        @test isapprox(d,[0.0,1.0,sqrt(2),1.0],atol=eps_level)
+    end
+
+    @testset "Multi-face meshes" begin
+        # Bowtie triangle set
+        V = Point3{Float64}[[0.0, 1.0, 0.0], [0.0, -1.0, 0.0], [1.0, 0.0, 0.0], 
+                            [2.0, 1.0, 0.0], [2.0, -1.0, 0.0]]
+        F = TriangleFace{Int64}[TriangleFace(1, 2, 3), TriangleFace(5, 4, 3)]
+        d,dd,l = distmarch(F,V,[1])
+        @test isapprox(d,[0.0,2.0,sqrt(2),2*sqrt(2),2*sqrt(2)],atol=eps_level)
+
+        # Two disconnected triangles
+        V = Point3{Float64}[[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [1.0, 1.0, 0.0], 
+                            [0.0, 0.0, 1.0], [1.0, 0.0, 1.0], [1.0, 1.0, 1.0]]
+        F = TriangleFace{Int64}[TriangleFace(1, 2, 3), TriangleFace(4, 5, 6)]
+        d,dd,l = distmarch(F,V,[1])
+        r = [0.0,1.0,sqrt(2),NaN,NaN,NaN]
+        b = .!isnan.(r)
+        @test isapprox(d[b],r[b],atol=eps_level) # Check reached entries
+        @test all(isnan.(d[.!b])) # Now check NaNs
+
+        # Two disconnected quads
+        V = Point3{Float64}[[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [1.0, 1.0, 0.0], [0.0,1.0,0.0],
+                            [0.0, 0.0, 1.0], [1.0, 0.0, 1.0], [1.0, 1.0, 1.0], [0.0,1.0,1.0]]
+        F = QuadFace{Int64}[QuadFace(1, 2, 3, 4), QuadFace(5, 6, 7, 8)]
+        d,dd,l = distmarch(F,V,[1])
+        r = [0.0,1.0,sqrt(2),1.0,NaN,NaN,NaN,NaN]
+        b = .!isnan.(r)
+        @test isapprox(d[b],r[b],atol=eps_level) # Check reached entries
+        @test all(isnan.(d[.!b])) # Now check NaNs
+
+        # Single cube
+        r = 2 * sqrt(3) / 2
+        M = cube(r)
+        F = faces(M)
+        V = coordinates(M)
+        d,dd,l = distmarch(F,V,[1])
+        @test isapprox(d,[0.0, 2.0, 2.8284271247461903, 2.0, 2.0,
+                         2.8284271247461903, 4.82842712474619, 2.8284271247461903],atol=eps_level) 
+
+        # Triangulated sphere, distance should approximate π 
+        r = 1.0
+        F,V = geosphere(4,r)
+        z = [v[3] for v ∈ V]
+        indStart =[findmin(z)[2]]
+        d,dd,l = distmarch(F,V,indStart)
+        @test isapprox(maximum(d),π,atol=0.01)
+
+        # Quadrangulated sphere, distance should approximate π 
+        r = 1.0
+        F,V = quadsphere(4,r)
+        z = [v[3] for v ∈ V]
+        indStart =[findmin(z)[2]]
+        d,dd,l = distmarch(F,V,indStart)
+        @test isapprox(maximum(d),π,atol=0.01)
+    end
 end
 
+
 @testset "separate vertices" begin
 
     eps_level = 0.001