init quality-based meshing algo

Author: sjkelly

src/DistMesh.jl

@@ -8,14 +8,15 @@ using GeometryBasics: Triangle, Tetrahedron, Mesh, Polytope, Point
 const tetpairs = ((1,2),(1,3),(1,4),(2,3),(2,4),(3,4))
 const tettriangles = ((1,2,3),(1,2,4),(2,3,4),(1,3,4))
 
+abstract type AbstractDistMeshAlgorithm end
 
 """
     DistMeshSetup
 
     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 DistMeshSetup{T}
+struct DistMeshSetup{T} <: AbstractDistMeshAlgorithm
     iso::T
     deltat::T
     ttol::T
@@ -36,6 +37,32 @@ function DistMeshSetup(;iso=0,
                      distribution)
 end
 
+"""
+    DistMeshQuality
+
+    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
+    distribution::Symbol # intial point distribution
+end
+
+function DistMeshQuality(;iso=0,
+                        ptol=.001,
+                        deltat=0.05,
+                        ttol=0.02,
+                        distribution=:regular)
+    T = promote_type(typeof(iso),typeof(ptol),typeof(deltat), typeof(ttol))
+    DistMeshQuality{T}(iso,
+                     deltat,
+                     ttol,
+                     ptol,
+                     distribution)
+end
+
 """
     DistMeshStatistics
 

src/distmeshnd.jl

@@ -18,7 +18,7 @@
 function distmesh(fdist::Function,
                   fh::Union{Function,HUniform},
                   h::Number,
-                  setup::DistMeshSetup=DistMeshSetup();
+                  setup::AbstractDistMeshAlgorithm=DistMeshSetup();
                   origin=GeometryBasics.Point{3,Float64}(-1,-1,-1),
                   widths=GeometryBasics.Point{3,Float64}(2,2,2),
                   fix=nothing,
@@ -209,3 +209,177 @@ function distmesh(fdist::Function,
         end
     end
 end
+
+function distmesh(fdist::Function,
+                  fh,
+                  h::Number,
+                  setup::DistMeshQuality,
+                  origin,
+                  widths,
+                  fix,
+                  ::Val{stats},
+                  ::Type{VertType}) where {VertType, stats}
+
+    geps=1e-1*h+setup.iso # parameter for filtering tets outside bounds and considering for max displacment of a node
+
+    # static parameter info
+    non_uniform = isa(fh, Function) # so we can elide fh calls
+    have_fixed = !isa(fix, Nothing)
+
+    #ptol=.001; ttol=.1; L0mult=1+.4/2^(dim-1); deltat=.2; geps=1e-1*h;
+
+    # initialize Vertex Arrays
+    # the first N points in the array correspond to'fix' points that do not move
+    if have_fixed
+        p = copy(fix)
+    else
+        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
+
+    # initialize arrays
+    pair_set = Set{Tuple{Int32,Int32}}()        # set used for ensure we have a unique set of edges
+    pair = Tuple{Int32,Int32}[]                 # edge indices (Int32 since we use Tetgen)
+    dp = zeros(VertType, length(p))             # displacement at each node
+    bars = VertType[]                           # the vector of each edge
+    L = eltype(VertType)[]                      # vector length of each edge
+    non_uniform && (L0 = eltype(VertType)[])    # desired edge length computed by dh (edge length function)
+    t = GeometryBasics.SimplexFace{4,Int32}[]   # tetrahedra indices from delaunay triangulation
+    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
+    statsdata = DistMeshStatistics()
+    lcount = 0 # iteration counter
+
+    @inbounds while true
+        # if large move, retriangulation
+        if maxmove>setup.ttol*h
+
+            delaunayn!(fdist, p, t, geps) # compute a new delaunay triangulation
+
+            tet_to_edges!(pair, pair_set, t) # Describe each edge by a unique pair of nodes
+
+            # resize arrays for new pair counts
+            resize!(bars, length(pair))
+            resize!(L, length(pair))
+            non_uniform && resize!(L0, length(pair))
+
+            stats && push!(statsdata.retriangulations, lcount)
+        end
+
+        # 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)
+            L0_f = non_uniform ? L0[i].*lscbrt : lscbrt
+            # compute force vectors
+            F = L0_f-L[i]
+            # edges are not allowed to pull, only repel
+            if F > zero(eltype(L))
+                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
+
+        # Zero out forces on fix points
+        if have_fixed
+            for i in eachindex(fix)
+                dp[i] = zero(VertType)
+            end
+        end
+
+        # apply point forces and
+        # bring outside points back to the boundary
+        maxdp = typemin(eltype(VertType))
+        maxmove = typemin(eltype(VertType))
+        for i in eachindex(p)
+
+            p0 = p[i] # store original point location
+            p[i] = p[i].+setup.deltat.*dp[i] # apply displacements to points
+
+            # Check if we are verifiably within the bounds and use this value
+            # to avoid recomputing fdist. This increases performance greatly on
+            # complex distance functions and large node counts
+            move = sqrt(sum((p[i]-p0).^2))          # compute movement from the displacement
+            d_est = pt_dists[i] + move              # apply the movement to our cache point
+            d = d_est < -geps ? d_est : fdist(p[i]) # determine if we need correct or approximate distance
+            pt_dists[i] = d                         # store distance
+
+            if d < -geps
+                maxdp = max(maxdp, setup.deltat*sqrt(sum(dp[i].^2)))
+            end
+
+            if d <= setup.iso
+                maxmove = max(move,maxmove)
+                continue
+            end
+
+            # bring points back to boundary if outside using central difference
+            p[i] = p[i] .- centraldiff(fdist,p[i]).*(d+setup.iso)
+            maxmove = max(sqrt(sum((p[i]-p0).^2)), maxmove)
+            pt_dists[i] = setup.iso # ideally
+        end
+
+        # increment iteration counter
+        lcount = lcount + 1
+
+        # save iteration stats
+        if stats
+            push!(statsdata.maxmove,maxmove)
+            push!(statsdata.maxdp,maxdp)
+            triangle_qualities!(tris,triset,qualities,p,t)
+            sort!(qualities) # sort for median calc and robust summation
+            mine, maxe = extrema(qualities)
+            push!(statsdata.average_qual, sum(qualities)/length(qualities))
+            push!(statsdata.median_qual, qualities[round(Int,length(qualities)/2)])
+            push!(statsdata.minimum_qual, mine)
+            push!(statsdata.maximum_qual, maxe)
+        end
+
+        # Termination criterion
+        if maxdp<setup.ptol*h
+            return p, t, statsdata
+        end
+    end
+end