Merge pull request #19 from JuliaGeometry/sjk/edge_qual

Sjk/edge qual

Author: sjkelly

bench/benchmarks.jl

@@ -4,8 +4,8 @@ using Plots
 using Dates
 using Descartes
 using GeometryBasics
-using JLD
-using HDF5
+#using JLD
+#using HDF5
 
 include("util.jl")
 
@@ -60,31 +60,31 @@ println("Benchmarking DistMesh.jl...")
 
 timestamp = Dates.format(now(), "yyyy-mm-ddTHH_MM")
 # generate an orthogonal test suite
-for ttol in 0.01:0.01:0.05, deltat in 0.02:0.02:0.2, el = 0.05:0.05:0.2
+for ttol in 0.01:0.01:0.05, deltat in 0.05:0.05:0.1, el = 0.1:0.05:0.2
     rt = time()
-    p,t,s = distmesh(torus,
-                     huniform,
+    result = distmesh(torus,
+                     HUniform(),
                      el,
-                     DistMeshSetup(deltat=deltat, retriangulation_criteria=RetriangulateMaxMove(ttol)),
+                     DistMeshSetup(deltat=deltat,ttol=ttol,distribution=:packed),
                      origin = GeometryBasics.Point{3,Float64}(-2),
                      widths = GeometryBasics.Point{3,Float64}(4),
-                     stats=true,
-                     distribution=:packed)
+                     stats=true)
     running_time = time() - rt # approximate, since we mostly care about convergence factors
     item = "torus$timestamp"
-    folder = joinpath(@__DIR__, "output/$item")
+    folder = joinpath(@__DIR__, "output")
+    !isdir(folder) && mkdir(folder)
+    folder = joinpath(folder, "$item")
     !isdir(folder) && mkdir(folder)
     param_str = "_ttol=$(ttol)_deltat=$(deltat)_el=$(el)"
     # save plots
-    plotout(s, DistMesh.triangle_qualities(p,t), folder, param_str)
+    plotout(result.stats, DistMesh.triangle_qualities(result.points,result.tetrahedra), folder, param_str)
     # save dataset as JLD
-    jldopen("$folder/$item.jld", "w") do file
-        g = g_create(file, param_str)
-        g["points"] = p
-        g["tets"] = t
-        g["stats"] = s
-        g["running_time"] = running_time
-    end
+    # jldopen("$folder/$item.jld", "w") do file
+    #     g = g_create(file, param_str)
+    #     g["points"] = p
+    #     g["tets"] = t
+    #     g["stats"] = s
+    #     g["running_time"] = running_time
+    # end
     println(param_str)
 end
-

bench/util.jl

@@ -2,16 +2,30 @@
 function plotout(statsdata, qualities, folder, name)
 
     qual_hist = Plots.histogram(qualities, title = "Quality", bins=30, legend=false)
-    avg_plt = Plots.plot(statsdata.average_qual, title = "Average Quality", legend=false, ylabel="Quality")
-    vline!(statsdata.retriangulations, line=(0.2, :dot, [:red]))
-    med_plt = Plots.plot(statsdata.median_qual, title = "Median Quality", legend=false, ylabel="Quality")
+    # avg_plt = Plots.plot(statsdata.average_qual, title = "Average Tri Quality", legend=false, ylabel="Quality")
+    # vline!(statsdata.retriangulations, line=(0.2, :dot, [:red]))
+
+    # med_plt = Plots.plot(statsdata.median_qual, title = "Median Tri Quality", legend=false, ylabel="Quality")
+    # vline!(statsdata.retriangulations, line=(0.2, :dot, [:red]))
+
+    # min_plt = Plots.plot(statsdata.minimum_qual, title = "Minimum Tri Quality", legend=false, ylabel="Quality")
+    # vline!(statsdata.retriangulations, line=(0.2, :dot, [:red]))
+
+    # max_plt = Plots.plot(statsdata.maximum_qual, title = "Maximum Tri Quality", legend=false, ylabel="Quality")
+    # vline!(statsdata.retriangulations, line=(0.2, :dot, [:red]))
+    data = hcat(statsdata.average_volume_edge_ratio, statsdata.min_volume_edge_ratio, statsdata.max_volume_edge_ratio)
+
+    tq_plt = Plots.plot(data, title = "Tet Quality", legend=true, label=["Avg","Min","Max"], ylabel="Vol/Edge Ratio")
     vline!(statsdata.retriangulations, line=(0.2, :dot, [:red]))
+
     maxdp_plt = Plots.plot(statsdata.maxdp, title = "Max Displacement", legend=false, ylabel="Edge Displacement")
     vline!(statsdata.retriangulations, line=(0.2, :dot, [:red]))
+
     maxmove_plt = Plots.plot(statsdata.maxmove, title = "Max Move", legend=false, ylabel="Point Displacement")
     vline!(statsdata.retriangulations, line=(0.2, :dot, [:red]))
-    plt = Plots.plot(avg_plt, med_plt,maxdp_plt,maxmove_plt,layout=(2,2), xlabel="Iteration")
+
+    plt = Plots.plot(tq_plt,maxdp_plt,maxmove_plt,layout=(3,1), xlabel="Iteration")
 
     savefig(plt, "$folder/result_stat$name.svg")
     savefig(qual_hist, "$folder/result_qual$name.svg")
-end
+end

docs/src/index.md

@@ -1,5 +1,50 @@
 # DistMesh.jl
 
+
+DistMesh.jl implements simplex refinement on signed distance functions, or anything that
+has a sign, distance, and called like a function. The algorithm was first presented
+in 2004 by Per-Olof Persson, and was initially a port of the corresponding Matlab Code.
+
+## What is Simplex Refinement?
+
+In layman's terms, a simplex is either a triangle in the 2D case, or a tetrahedra in the 3D case.
+
+When simulating, you other want a few things from a mesh of simplices:
+    - Accurate approximation of boundaries and features
+    - Adaptive mesh sizes to improve accuracy
+    - Near-Regular Simplicies
+
+DistMesh is designed to address the above.
+
+## Algorithm Overview
+
+The basic processes is as follows:
+
+
+
+## Comparison to other refinements
+
+DistMesh generally has a very low memory footprint, and can refine without additional
+memory allocation. Similarly, since the global state of simplex qualities is accounted for
+in each refinement iteration, this leads to very high quality meshes.
+
+Aside from the above, since DistMesh works on signed distance functions it can handle
+complex and varied input data that are not in the form of surface meshes (Piecewise Linear Complicies).
+
+## Difference from the MatLab implementation
+
+Given the same parameters, the Julia implementation of DistMesh will generally perform
+4-60 times faster than the MatLab implementation. Delaunay Triangulation in MatLab uses
+QHull, whereas DistMesh.jl uses TetGen.
+
+## How do I get a Signed Distance Function?
+
+Here are some libraries that turn gridded and level set data into an approximate signed
+distance function:
+
+- Interpolations.jl
+- AdaptiveDistanceFields.jl
+
 ```@index
 ```
 

src/DistMesh.jl

@@ -13,7 +13,7 @@ abstract type AbstractDistMeshAlgorithm end
 """
     DistMeshSetup
 
-    Takes Keyword arguments as follows:
+Takes Keyword arguments as follows:
 
     iso (default: 0): Value of which to extract the isosurface, inside surface is negative
     deltat (default: 0.1): the fraction of edge displacement to apply each iteration
@@ -54,27 +54,39 @@ end
 """
     DistMeshQuality
 
+Use Tetrahedral quality analysis to control the meshing process
+
     iso (default: 0): Value of which to extract the iso surface, inside negative
     deltat (default: 0.1): the fraction of edge displacement to apply each iteration
 """
 struct DistMeshQuality{T} <: AbstractDistMeshAlgorithm
     iso::T
     deltat::T
-    minimum::T
+    filter_less_than::T # Remove tets less than the given quality
+    #allow_n_regressions::Int # Might want this
+    termination_quality::T # Once achieved, terminate
+    sort::Bool # use hilbert sort to cache-localize points
+    sort_interval::Int # retriangulations before resorting
+    nonlinear::Bool # uses nonlinear edge force
     distribution::Symbol # intial point distribution
 end
 
 function DistMeshQuality(;iso=0,
-                        ptol=.001,
                         deltat=0.05,
-                        ttol=0.02,
+                        filter_less_than=0.02,
+                        termination_quality=0.3,
+                        sort=false,
+                        sort_interval=20,
+                        nonlinear=true,
                         distribution=:regular)
-    T = promote_type(typeof(iso),typeof(ptol),typeof(deltat), typeof(ttol))
-    DistMeshQuality{T}(iso,
-                     deltat,
-                     ttol,
-                     ptol,
-                     distribution)
+    DistMeshQuality(iso,
+                    deltat,
+                    filter_less_than,
+                    termination_quality,
+                    sort,
+                    sort_interval,
+                    nonlinear,
+                    distribution)
 end
 
 """
@@ -85,21 +97,20 @@ end
 struct DistMeshStatistics{T}
     maxmove::Vector{T} # max point move in an iteration
     maxdp::Vector{T} # max displacmeent induced by an edge
-    average_qual::Vector{T}
-    median_qual::Vector{T}
-    minimum_qual::Vector{T}
-    maximum_qual::Vector{T}
+    min_volume_edge_ratio::Vector{T}
+    max_volume_edge_ratio::Vector{T}
+    average_volume_edge_ratio::Vector{T}
     retriangulations::Vector{Int} # Iteration num where retriangulation occured
 end
 
-DistMeshStatistics() = DistMeshStatistics{Float64}([],[],[],[],[],[],[])
+DistMeshStatistics() = DistMeshStatistics{Float64}([],[],[],[],[],[])
 
 """
 Uniform edge length function.
 """
 struct HUniform end
 
-include("common.jl")
+
 include("diff.jl")
 include("pointdistribution.jl")
 include("distmeshnd.jl")

src/common.jl

@@ -73,15 +73,10 @@ function distmesh(fdist::Function,
         p = VertType[]
     end
 
-    pt_dists = map(fdist, p) # cache to store point locations so we can minimize fdist calls
-
-    # add points to p based on the initial distribution
-    if setup.distribution === :regular
-        simplecubic!(fdist, p, pt_dists, h, setup.iso, origin, widths, VertType)
-    elseif setup.distribution === :packed
-        # face-centered cubic point distribution
-        facecenteredcubic!(fdist, p, pt_dists, h, setup.iso, origin, widths, VertType)
-    end
+    pt_dists = eltype(VertType)[]
+
+    # setup the initial point distribution specified in setup
+    point_distribution!(fdist,p,pt_dists,h, setup, origin, widths, VertType)
 
     # Result struct for holding points, simplices, and iteration statistics
     result = DistMeshResult(p,
@@ -97,11 +92,6 @@ function distmesh(fdist::Function,
     L0 = non_uniform ? eltype(VertType)[] : nothing # desired edge length computed by dh (edge length function)
     maxmove = typemax(eltype(VertType))         # stores an iteration max movement for retriangulation
 
-    # arrays for tracking quality metrics
-    tris = NTuple{3,Int32}[]        # triangles used for quality checks
-    triset = Set{NTuple{3,Int32}}() # set for triangles to ensure uniqueness
-    qualities = eltype(VertType)[]
-
     # information on each iteration
     lcount = 0 # iteration counter
     triangulationcount = 0 # triangulation counter
@@ -172,13 +162,10 @@ function distmesh(fdist::Function,
         if stats
             push!(result.stats.maxmove,maxmove)
             push!(result.stats.maxdp,maxdp)
-            triangle_qualities!(tris,triset,qualities,result.points,result.tetrahedra)
-            sort!(qualities) # sort for median calc and robust summation
-            mine, maxe = extrema(qualities)
-            push!(result.stats.average_qual, sum(qualities)/length(qualities))
-            push!(result.stats.median_qual, qualities[round(Int,length(qualities)/2)])
-            push!(result.stats.minimum_qual, mine)
-            push!(result.stats.maximum_qual, maxe)
+            min_v_edge, avg_v_edge, max_v_edge = volume_edge_stats(result.points,result.tetrahedra)
+            push!(result.stats.min_volume_edge_ratio, min_v_edge)
+            push!(result.stats.average_volume_edge_ratio, avg_v_edge)
+            push!(result.stats.max_volume_edge_ratio, max_v_edge)
         end
 
         # Termination criterion
@@ -226,3 +213,51 @@ function retriangulate!(fdist, result::DistMeshResult, geps, setup, triangulatio
     j <= lastindex(t) && resize!(t, j-1)
     nothing
 end
+
+
+function compute_displacements!(fh, dp, pair, L, L0, bars, p, setup,
+    ::Type{VertType}) where {VertType}
+
+    non_uniform = isa(typeof(L0), AbstractVector)
+
+    # compute edge lengths (L) and adaptive edge lengths (L0)
+    # Lp norm (p=3) is partially computed here
+    Lsum = zero(eltype(L))
+    L0sum = non_uniform ? zero(eltype(L0)) : length(pair)
+    for i in eachindex(pair)
+        pb = pair[i]
+        b1 = p[pb[1]]
+        b2 = p[pb[2]]
+        barvec = b1 - b2 # bar vector
+        bars[i] = barvec
+        L[i] = sqrt(sum(barvec.^2)) # length
+        non_uniform && (L0[i] = fh((b1+b2)./2))
+        Lsum = Lsum + L[i].^3
+        non_uniform && (L0sum = L0sum + L0[i].^3)
+    end
+
+    # zero out force at each node
+    for i in eachindex(dp)
+        dp[i] = zero(VertType)
+    end
+
+    # this is not hoisted correctly in the loop so we initialize here
+    # finish computing the Lp norm (p=3)
+    lscbrt = (1+(0.4/2^2))*cbrt(Lsum/L0sum)
+
+    # Move mesh points based on edge lengths L and forces F
+    for i in eachindex(pair)
+        if non_uniform && L[i] < L0[i]*lscbrt || L[i] < lscbrt
+            L0_f = non_uniform ? L0[i].*lscbrt : lscbrt
+            # compute force vectors
+            F = setup.nonlinear ? (L[i]+L0_f)*(L0_f-L[i])/(2*L0_f) : L0_f-L[i]
+            # edges are not allowed to pull, only repel
+            FBar = bars[i].*F./L[i]
+            # add the force vector to the node
+            b1 = pair[i][1]
+            b2 = pair[i][2]
+            dp[b1] = dp[b1] .+ FBar
+            dp[b2] = dp[b2] .- FBar
+        end
+    end
+end