From 3bbdf96076f1d743dd527ff3da9b2d1303f36f55 Mon Sep 17 00:00:00 2001 From: steve Date: Mon, 14 Oct 2019 21:48:48 -0400 Subject: [PATCH 1/9] up fea for distmesh --- examples/fea.jl | 43 +++++++++++++++++++++++++++++++------------ 1 file changed, 31 insertions(+), 12 deletions(-) diff --git a/examples/fea.jl b/examples/fea.jl index f467597..0bb9734 100644 --- a/examples/fea.jl +++ b/examples/fea.jl @@ -1,20 +1,39 @@ using Descartes +using DistMesh +using StaticArrays +using Makie # params -beam_size = [50,10,10] -hole_ct = 5 -hole_d = 3 +function beam(beam_size = [50,10,10], + hole_ct = 5, + hole_d = 3) -# computations -hole_interval = beam_size[1]/(hole_ct + 1) + # computations + hole_interval = beam_size[1]/(hole_ct + 1) -c = Cuboid(beam_size) + c = Cuboid(beam_size) -for i = 1:hole_ct - h = translate([hole_interval*i, -1, beam_size[3]/2])* - rotate(-pi/2, [1,0,0])* - Cylinder(hole_d/2, beam_size[2]+2, center=false) - global c = diff(c, h) + for i = 1:hole_ct + h = translate([hole_interval*i, -1, beam_size[3]/2])* + rotate(-pi/2, [1,0,0])* + Cylinder(hole_d/2, beam_size[2]+2, center=false) + c = diff(c, h) + end + c end -save("fea.ply",HomogenousMesh(c)) +function deadmau5() + diff(union(Descartes.Sphere(5), + translate([4,2,0])Descartes.Sphere(4), + translate([-4,2,0])Descartes.Sphere(4)), + translate([5,2,3])Descartes.Sphere(3), + translate([-5,2,3])Descartes.Sphere(3)) +end + +b = deadmau5() +f(x) = FRep(b, x) +h = HyperRectangle(b) +@show f(SVector(1,1,1)) +p, t = distmeshnd(f, huniform, 1.5, origin=h.origin, widths=h.widths) + +#save("fea.ply",HomogenousMesh(c)) From fc643c0d489209217339195856cf2bddff9d5a1b Mon Sep 17 00:00:00 2001 From: steve Date: Sun, 20 Oct 2019 23:50:43 -0400 Subject: [PATCH 2/9] workspace --- examples/fea.jl | 41 +++++++++++++++++++++++++++++++++++++++-- 1 file changed, 39 insertions(+), 2 deletions(-) diff --git a/examples/fea.jl b/examples/fea.jl index 0bb9734..9790fe3 100644 --- a/examples/fea.jl +++ b/examples/fea.jl @@ -1,7 +1,9 @@ using Descartes using DistMesh +using GeometryTypes using StaticArrays using Makie +#using Colors # params function beam(beam_size = [50,10,10], @@ -16,6 +18,7 @@ function beam(beam_size = [50,10,10], for i = 1:hole_ct h = translate([hole_interval*i, -1, beam_size[3]/2])* rotate(-pi/2, [1,0,0])* + rotate(-pi/6, [0,1,0])* # adversarial case Cylinder(hole_d/2, beam_size[2]+2, center=false) c = diff(c, h) end @@ -30,10 +33,44 @@ function deadmau5() translate([-5,2,3])Descartes.Sphere(3)) end -b = deadmau5() +b = beam([10,10,10],1, 3) f(x) = FRep(b, x) +m = GLNormalMesh(b) +scene = mesh(m, color=:blue) +display(scene) +sleep(5) h = HyperRectangle(b) @show f(SVector(1,1,1)) -p, t = distmeshnd(f, huniform, 1.5, origin=h.origin, widths=h.widths) +p, t = distmeshnd(f, huniform, 0.7, origin=h.origin, widths=h.widths, vis=false) + +VertType = eltype(p) +pair = Tuple{Int32,Int32}[] # edge indices (Int32 since we use Tetgen) + +pair=resize!(pair, length(t)*6) +ls = Pair{VertType,VertType}[] + +for i in eachindex(t) + for ep in 1:6 + p1 = t[i][DistMesh.tetpairs[ep][1]] + p2 = t[i][DistMesh.tetpairs[ep][2]] + if p1 > p2 + pair[(i-1)*6+ep] = (p2,p1) + else + pair[(i-1)*6+ep] = (p1,p2) + end + end +end + +sort!(pair) +unique!(pair) + +# makie vis + +resize!(ls, length(pair)) +for i = 1:length(pair) + ls[i] = p[pair[i][1]] => p[pair[i][2]] +end +scene = Makie.linesegments(ls) +display(scene) #save("fea.ply",HomogenousMesh(c)) From 7ade50a94b7e01744159dfabe5db68f7838472e8 Mon Sep 17 00:00:00 2001 From: steve Date: Sun, 20 Oct 2019 23:51:34 -0400 Subject: [PATCH 3/9] clarify triangles --- docs/src/design.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/design.md b/docs/src/design.md index 606c193..7203f84 100644 --- a/docs/src/design.md +++ b/docs/src/design.md @@ -108,4 +108,4 @@ Any primitve or operation must implement two core operations; `FRep` and `HyperR In the next line we actually perform the meshing operation. We call the meshtype with the implicit function we created, `f`, in the bounds generated by `HyperRectangle`, uniformly sample the space by `samples`, and actually generate the triangular meshing using `algorithm`. In this case, the defaults are `samples=(128,128,128)` and `algorithm=MarchingCubes()`. -In the loop `for i = 2:length(primitives) ... end` we handle additional primitives passed to our meshing function so you can create a single mesh output from several different objects. In order to maintain performance and resolution these are not `union`ed. +In the loop `for i = 2:length(primitives) ... end` we handle additional primitives passed to our meshing function so you can create a single mesh output from several different objects. In order to maintain performance and resolution these are not `union`ed. This means if primitives are passed as seperate arguments ther will each get their own triangles. From 7447d8bd13c1565d330bf4e0639cd7c581665bea Mon Sep 17 00:00:00 2001 From: steve Date: Sun, 20 Oct 2019 23:52:27 -0400 Subject: [PATCH 4/9] try true-ish distance function for primitives --- src/constructors.jl | 4 ++++ src/frep.jl | 23 +++++++++++++++++------ 2 files changed, 21 insertions(+), 6 deletions(-) diff --git a/src/constructors.jl b/src/constructors.jl index e913894..1043094 100644 --- a/src/constructors.jl +++ b/src/constructors.jl @@ -54,6 +54,8 @@ end # # Union +@nospecialize + function CSGUnion(l::AbstractPrimitive{N1,T1}, r::AbstractPrimitive{N2,T2}) where {N1, N2, T1, T2} N1 == N2 || error("cannot create CSG between objects in R$N1 and R$N2") return CSGUnion{N1,T1, typeof(l), typeof(r)}(l,r) @@ -132,3 +134,5 @@ end function LinearExtrude(p::AbstractPrimitive{2,T}, d) where {T} LinearExtrude{3, T, typeof(p)}(p, convert(T, d), SMatrix{4,4}(1.0*I), SMatrix{4,4}(1.0*I)) end + +@specialize \ No newline at end of file diff --git a/src/frep.jl b/src/frep.jl index d8efa8b..8aecce9 100644 --- a/src/frep.jl +++ b/src/frep.jl @@ -1,5 +1,16 @@ # http://en.wikipedia.org/wiki/Function_representation +function fmax(x::T1,y::T2) where {T1, T2} + if x > zero(T1) && y > zero(T2) + return sqrt(x^2+y^2) + end + return max(x,y) +end + +fmax(x,y,z) = fmax(fmax(x,y),z) +fmax(a,b,c,d) = fmax(fmax(a,b),fmax(c,d)) +fmax(a,b,c,d,e,f) = fmax(fmax(a,b,c,d),fmax(e,f)) + function _radius(a,b,r) if abs(a-b) >= r return min(a,b) @@ -21,7 +32,7 @@ function FRep(p::Cylinder, v) x = v[1]*it[1,1]+v[2]*it[1,2]+v[3]*it[1,3]+it[1,4] y = v[1]*it[2,1]+v[2]*it[2,2]+v[3]*it[2,3]+it[2,4] z = v[1]*it[3,1]+v[2]*it[3,2]+v[3]*it[3,3]+it[3,4] - max(-z+p.bottom, z-p.height-p.bottom, sqrt(x*x + y*y) - p.radius) + fmax(-z+p.bottom, z-p.height-p.bottom, sqrt(x*x + y*y) - p.radius) end function FRep(p::Circle, v) @@ -39,7 +50,7 @@ function FRep(p::Cuboid, v) z = v[1]*it[3,1]+v[2]*it[3,2]+v[3]*it[3,3]+it[3,4] dx, dy, dz = p.dimensions lbx, lby,lbz = p.lowercorner - max(-x+lbx, x-dx-lbx, + fmax(-x+lbx, x-dx-lbx, -y+lby, y-dy-lby, -z+lbz, z-dz-lbz) end @@ -50,7 +61,7 @@ function FRep(p::Square, v) y = v[1]*it[2,1]+v[2]*it[2,2]+it[2,3] dx, dy = p.dimensions lbx, lby = p.lowercorner - max(-x+lbx, x-dx-lbx, + fmax(-x+lbx, x-dx-lbx, -y+lby, y-dy-lby) end @@ -59,11 +70,11 @@ function FRep(u::CSGUnion, v) end function FRep(u::CSGDiff, v) - max(FRep(u.left, v), -FRep(u.right, v)) + fmax(FRep(u.left, v), -FRep(u.right, v)) end function FRep(u::CSGIntersect, v) - max(FRep(u.left, v), FRep(u.right, v)) + fmax(FRep(u.left, v), FRep(u.right, v)) end function FRep(s::Shell, v) @@ -107,5 +118,5 @@ function FRep(p::LinearExtrude, v) y = v[1]*it[2,1]+v[2]*it[2,2]+v[3]*it[2,3]+it[2,4] z = v[1]*it[3,1]+v[2]*it[3,2]+v[3]*it[3,3]+it[3,4] r = FRep(p.primitive, v) - max(max(-z,z-p.distance), r) + fmax(fmax(-z,z-p.distance), r) end From da221ee4413966318d01f93de19a5a5a2ad1aa55 Mon Sep 17 00:00:00 2001 From: steve Date: Tue, 22 Oct 2019 23:53:13 -0400 Subject: [PATCH 5/9] test with packed spacing --- examples/fea.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/examples/fea.jl b/examples/fea.jl index 9790fe3..77dca8e 100644 --- a/examples/fea.jl +++ b/examples/fea.jl @@ -5,6 +5,8 @@ using StaticArrays using Makie #using Colors +Cylinder=Descartes.Cylinder + # params function beam(beam_size = [50,10,10], hole_ct = 5, @@ -41,7 +43,7 @@ display(scene) sleep(5) h = HyperRectangle(b) @show f(SVector(1,1,1)) -p, t = distmeshnd(f, huniform, 0.7, origin=h.origin, widths=h.widths, vis=false) +@time p, t = distmeshnd(f, huniform, 0.5, origin=h.origin, widths=h.widths, vis=false, distribution=:packed) VertType = eltype(p) pair = Tuple{Int32,Int32}[] # edge indices (Int32 since we use Tetgen) From cb635fe38c8330857046d5ffd1e8215f5aa87489 Mon Sep 17 00:00:00 2001 From: Steve Kelly Date: Fri, 1 Nov 2019 10:57:04 -0400 Subject: [PATCH 6/9] science --- examples/fea.jl | 23 +++++++++++++++++++++-- 1 file changed, 21 insertions(+), 2 deletions(-) diff --git a/examples/fea.jl b/examples/fea.jl index 77dca8e..3b7ae82 100644 --- a/examples/fea.jl +++ b/examples/fea.jl @@ -35,7 +35,8 @@ function deadmau5() translate([-5,2,3])Descartes.Sphere(3)) end -b = beam([10,10,10],1, 3) +#b = beam([10,10,10],1, 3) +b = Descartes.Sphere(5) f(x) = FRep(b, x) m = GLNormalMesh(b) scene = mesh(m, color=:blue) @@ -43,7 +44,8 @@ display(scene) sleep(5) h = HyperRectangle(b) @show f(SVector(1,1,1)) -@time p, t = distmeshnd(f, huniform, 0.5, origin=h.origin, widths=h.widths, vis=false, distribution=:packed) +statsdata = DistMesh.DistMeshStatistics() +@time p, t = distmesh(f, huniform, 1, origin=h.origin, widths=h.widths, stats=true, statsdata=statsdata, distribution=:regular) VertType = eltype(p) pair = Tuple{Int32,Int32}[] # edge indices (Int32 since we use Tetgen) @@ -75,4 +77,21 @@ end scene = Makie.linesegments(ls) display(scene) +qualities = DistMesh.triangle_qualities(p,t) + +statsdata +using Plots + +plt = Plots.histogram(qualities, title = "Quality", legend=false) +savefig(plt, "qual_hist.png") +plt = Plots.plot(statsdata.average_qual, title = "Average Quality", legend=false) +savefig(plt, "average_qual.png") +plt = Plots.plot(statsdata.median_qual, title = "Median Quality", legend=false) +savefig(plt, "median_qual.png") +plt = Plots.plot(statsdata.maxdp, title = "Max Displacement", legend=false) +savefig(plt, "max_displacement.png") +plt = Plots.plot(statsdata.maxmove, title = "Max Move", legend=false) +savefig(plt, "max_move.png") + + #save("fea.ply",HomogenousMesh(c)) From b2b91060b168b004dc8a3f17b80c7561a8e9eb84 Mon Sep 17 00:00:00 2001 From: Steve Kelly Date: Fri, 1 Nov 2019 18:56:24 -0400 Subject: [PATCH 7/9] better plots --- examples/fea.jl | 23 +++++++++++++---------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/examples/fea.jl b/examples/fea.jl index 3b7ae82..dd65ac2 100644 --- a/examples/fea.jl +++ b/examples/fea.jl @@ -82,16 +82,19 @@ qualities = DistMesh.triangle_qualities(p,t) statsdata using Plots -plt = Plots.histogram(qualities, title = "Quality", legend=false) -savefig(plt, "qual_hist.png") -plt = Plots.plot(statsdata.average_qual, title = "Average Quality", legend=false) -savefig(plt, "average_qual.png") -plt = Plots.plot(statsdata.median_qual, title = "Median Quality", legend=false) -savefig(plt, "median_qual.png") -plt = Plots.plot(statsdata.maxdp, title = "Max Displacement", legend=false) -savefig(plt, "max_displacement.png") -plt = Plots.plot(statsdata.maxmove, title = "Max Move", legend=false) -savefig(plt, "max_move.png") +qual_hist = Plots.histogram(qualities, title = "Quality", legend=false) +avg_plt = Plots.plot(statsdata.average_qual, title = "Average Quality", legend=false) +vline!(statsdata.retriangulations, line=(0.2, :dot, [:red])) +med_plt = Plots.plot(statsdata.median_qual, title = "Median Quality", legend=false) +vline!(statsdata.retriangulations, line=(0.2, :dot, [:red])) +maxdp_plt = Plots.plot(statsdata.maxdp, title = "Max Displacement", legend=false) +vline!(statsdata.retriangulations, line=(0.2, :dot, [:red])) +maxmove_plt = Plots.plot(statsdata.maxmove, title = "Max Move", legend=false) +vline!(statsdata.retriangulations, line=(0.2, :dot, [:red])) +plt = Plots.plot(avg_plt, med_plt,maxdp_plt,maxmove_plt,layout=(2,2)) +savefig(plt, "result_stat.svg") +savefig(qual_hist, "result_qual.svg") + #save("fea.ply",HomogenousMesh(c)) From 280e66a1b9f8d64d8a0d9ebb4bfd7eb54e3c20b7 Mon Sep 17 00:00:00 2001 From: Steve Kelly Date: Fri, 1 Nov 2019 20:00:37 -0400 Subject: [PATCH 8/9] plot tweaks --- examples/fea.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/examples/fea.jl b/examples/fea.jl index dd65ac2..885a40d 100644 --- a/examples/fea.jl +++ b/examples/fea.jl @@ -83,15 +83,15 @@ statsdata using Plots qual_hist = Plots.histogram(qualities, title = "Quality", legend=false) -avg_plt = Plots.plot(statsdata.average_qual, title = "Average Quality", 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) +med_plt = Plots.plot(statsdata.median_qual, title = "Median Quality", legend=false, ylabel="Quality") vline!(statsdata.retriangulations, line=(0.2, :dot, [:red])) -maxdp_plt = Plots.plot(statsdata.maxdp, title = "Max Displacement", legend=false) +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) +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)) +plt = Plots.plot(avg_plt, med_plt,maxdp_plt,maxmove_plt,layout=(2,2), xlabel="Iteration", ) savefig(plt, "result_stat.svg") savefig(qual_hist, "result_qual.svg") From e36ee88134c6a7cce23d0747d068bb1c5c082436 Mon Sep 17 00:00:00 2001 From: Steve Kelly Date: Sun, 22 Dec 2019 00:29:46 -0500 Subject: [PATCH 9/9] wip --- examples/fea.jl | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/examples/fea.jl b/examples/fea.jl index 885a40d..1b27a0e 100644 --- a/examples/fea.jl +++ b/examples/fea.jl @@ -6,6 +6,8 @@ using Makie #using Colors Cylinder=Descartes.Cylinder +translate=Descartes.translate +rotate=Descartes.rotate # params function beam(beam_size = [50,10,10], @@ -35,8 +37,8 @@ function deadmau5() translate([-5,2,3])Descartes.Sphere(3)) end -#b = beam([10,10,10],1, 3) -b = Descartes.Sphere(5) +b = beam([10,10,10],1, 3) +#b = Descartes.Sphere(5) f(x) = FRep(b, x) m = GLNormalMesh(b) scene = mesh(m, color=:blue) @@ -44,9 +46,8 @@ display(scene) sleep(5) h = HyperRectangle(b) @show f(SVector(1,1,1)) -statsdata = DistMesh.DistMeshStatistics() -@time p, t = distmesh(f, huniform, 1, origin=h.origin, widths=h.widths, stats=true, statsdata=statsdata, distribution=:regular) - +@time p, t, statsdata = distmesh(f, HUniform(), 0.9, DistMeshSetup(distribution=:packed,sort_interval=20,deltat=0.05,nonlinear=true,sort=true), origin=h.origin, widths=h.widths, stats=true) +@show length(p), length(t) VertType = eltype(p) pair = Tuple{Int32,Int32}[] # edge indices (Int32 since we use Tetgen)