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.
diff --git a/examples/fea.jl b/examples/fea.jl
index f467597..1b27a0e 100644
--- a/examples/fea.jl
+++ b/examples/fea.jl
@@ -1,20 +1,101 @@
 using Descartes
+using DistMesh
+using GeometryTypes
+using StaticArrays
+using Makie
+#using Colors
+
+Cylinder=Descartes.Cylinder
+translate=Descartes.translate
+rotate=Descartes.rotate
 
 # 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)
+
+    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])*
+                rotate(-pi/6, [0,1,0])* # adversarial case
+                    Cylinder(hole_d/2, beam_size[2]+2, center=false)
+        c = diff(c, h)
+    end
+    c
+end
+
+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 = beam([10,10,10],1, 3)
+#b = Descartes.Sphere(5)
+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))
+@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)
 
-# computations
-hole_interval = beam_size[1]/(hole_ct + 1)
+pair=resize!(pair, length(t)*6)
+ls = Pair{VertType,VertType}[]
 
-c = Cuboid(beam_size)
+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)
 
-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)
+# 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)
+
+qualities = DistMesh.triangle_qualities(p,t)
+
+statsdata
+using Plots
+
+qual_hist = Plots.histogram(qualities, title = "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, ylabel="Quality")
+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", )
+savefig(plt, "result_stat.svg")
+savefig(qual_hist, "result_qual.svg")
+
+
 
-save("fea.ply",HomogenousMesh(c))
+#save("fea.ply",HomogenousMesh(c))
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