diff --git a/src/DistMesh.jl b/src/DistMesh.jl index cefeea6..20c8749 100644 --- a/src/DistMesh.jl +++ b/src/DistMesh.jl @@ -19,45 +19,36 @@ const tettriangles = ((1,2,3),(1,2,4),(2,3,4),(1,3,4)) 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,RT} +struct DistMeshSetup{T} iso::T deltat::T - retriangulation_criteria::RT ptol::T + ttol::T + maxmove_delta::Int + maxmove_delta_delay::Int droptets::Bool # drop tetrahedra with centroids outside of the boundary + initial_points::Symbol end -function DistMeshSetup(;iso=0, +function DistMeshSetup(;iso=0.0, ptol=.001, - deltat=0.05, - retriangulation_criteria=RetriangulateMaxMove(0.02), - droptets=true) + deltat=0.2, + ttol=0.1, + maxmove_delta=6, # number of prior samples to consider + maxmove_delta_delay=2, + droptets=true, + initial_points=:regular) T = promote_type(typeof(iso),typeof(ptol),typeof(deltat)) - DistMeshSetup{T,typeof(retriangulation_criteria)}(iso, - deltat, - retriangulation_criteria, - ptol, - droptets) -end - -abstract type AbstractRetriangulationCriteria end - -""" - retriangulate if a move over a given tolerance occurs -""" -struct RetriangulateMaxMove{T} <: AbstractRetriangulationCriteria - ttol::T -end - -""" - retriangulate on a positive delta over N moves, after M force iterations -""" -struct RetriangulateMaxMoveDelta <: AbstractRetriangulationCriteria - move_count::Int - iterations::Int + DistMeshSetup{T}(iso, + deltat, + ptol, + ttol, + maxmove_delta, + maxmove_delta_delay, + droptets, + initial_points) end - """ DistMeshStatistics @@ -88,7 +79,7 @@ include("decompositions.jl") #export distmeshsurface export distmesh, DistMeshSetup, DistMeshStatistics -export RetriangulateMaxMove +export RetriangulateMaxMove, RetriangulateMaxMoveDelta export huniform end # module diff --git a/src/distmeshnd.jl b/src/distmeshnd.jl index 1f2a7d6..e4fda73 100644 --- a/src/distmeshnd.jl +++ b/src/distmeshnd.jl @@ -15,11 +15,10 @@ d(p) = sqrt(sum(p.^2))-1 p,t = distmeshnd(d,huniform,0.2) """ -function distmesh(fdist::Function,fh::Function,h::Number, setup::DistMeshSetup{T,ReTri}=DistMeshSetup(), ::Type{VertType}=GeometryBasics.Point{3,Float64}; origin=VertType(-1,-1,-1), +function distmesh(fdist::Function,fh::Function,h::Number, setup::DistMeshSetup{T}=DistMeshSetup(), ::Type{VertType}=GeometryBasics.Point{3,Float64}; origin=VertType(-1,-1,-1), widths=VertType(2,2,2), fix::Vector{VertType}=VertType[], - stats=false, - distribution=:regular) where {VertType, T, ReTri} + stats=false) where {VertType, T} ptol=setup.ptol L0mult=1+.4/2^2 @@ -40,9 +39,9 @@ function distmesh(fdist::Function,fh::Function,h::Number, setup::DistMeshSetup{T # 'fix' points that do not move p = copy(fix) - if distribution == :regular + if setup.initial_points === :regular simplecubic!(fdist, p, h, setup.iso, origin, widths, VertType) - elseif distribution == :packed + elseif setup.intial_points === :packed # face-centered cubic point distribution facecenteredcubic!(fdist, p, h, setup.iso, origin, widths, VertType) end @@ -62,21 +61,24 @@ function distmesh(fdist::Function,fh::Function,h::Number, setup::DistMeshSetup{T tris = NTuple{3,Int32}[] # array to store triangles used for quality checks triset = Set{NTuple{3,Int32}}() # set for triangles to ensure uniqueness qualities = eltype(VertType)[] - #maxmoves = eltype(VertType)[] + maxmoves = eltype(VertType)[] # information on each iteration statsdata = DistMeshStatistics() + last_retri = 0 # iterations since last retriangulation @inbounds while true # Retriangulation by Delaunay - + if length(maxmoves) == setup.maxmove_delta + popfirst!(maxmoves) + end + lcount > 0 && push!(maxmoves, maxmove) # dont include the very first iteration # if large move, retriangulation - if ReTri <: RetriangulateMaxMove && maxmove>setup.retriangulation_criteria.ttol*h + if lcount < 0 && maxmove > setup.ttol*h || last_retri > setup.maxmove_delta_delay && maxmove > sum(maxmoves)/length(maxmoves) triangulation = delaunayn(p) t_d = triangulation.tetrahedra resize!(t, length(t_d)) copyto!(t, t_d) # we need to copy since we have a shared reference with tetgen - sort!(t) # sort tetrahedra so points are closer in mem # average points to get mid point of each tetrahedra # if the mid point of the tetrahedra is outside of @@ -115,9 +117,13 @@ function distmesh(fdist::Function,fh::Function,h::Number, setup::DistMeshSetup{T resize!(L, length(pair)) resize!(L0, length(pair)) + empty!(maxmoves) + last_retri = 0 stats && push!(statsdata.retriangulations, lcount) end + last_retri = last_retri + 1 + # 6. Move mesh points based on edge lengths L and forces F Lsum = zero(eltype(L)) L0sum = zero(eltype(L0)) diff --git a/test/runtests.jl b/test/runtests.jl index 977065b..bb83bcd 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -11,7 +11,7 @@ using StaticArrays @test length(p) == 485 @test length(t) == 2207 - p,t,_ = distmesh(d,huniform,0.2, distribution=:packed) + p,t,_ = distmesh(d,huniform,0.2, DistMeshSetup(initial_points=:packed)) @test length(p) == 742 @test length(t) == 3472 end