setminus is now part of the API

src/CompScienceMeshes.jl

@@ -21,7 +21,7 @@ export index
 export euclidianbasis, point
 
 # default mesh creation
-export mesh, readmesh, writemesh, meshgeo
+export mesh, readmesh, writemesh, meshgeo, setminus
 export meshsegment, meshrectangle, meshcircle, meshdisk, meshsphere, meshcuboid
 export tetmeshsphere, tetmeshcuboid
 export subdMesh

src/charts.jl

@@ -145,7 +145,7 @@ of vertices supplied minus one.
     end
 end
 
-simplex(vertices...) = simplex(SVector((vertices...,)))
+simplex(vertices::SVector{N,T}...) where {N,T<:Number} = simplex(SVector((vertices...,)))
 
 @generated function simplex(vertices, ::Type{Val{D}}) where D
     P = eltype(vertices)
@@ -320,7 +320,7 @@ function edges(s::Simplex{3,2})
         simplex(s.vertices[1], s.vertices[2])]
 end
 
-function faces(c)
+function faces(c::Simplex{3,2})
     @SVector[
         simplex(c[2], c[3]),
         simplex(c[3], c[1]),
@@ -337,6 +337,10 @@ function faces(c::CompScienceMeshes.Simplex{3,3,0,4,T}) where {T}
     ]
 end
 
+function faces(ch::Simplex{3,2}, ::Type{Val{0}})
+    simplex.(ch.vertices)
+end
+
 
 """
     ReferenceSimplex{Dimension, CoordType, NumVertices}
@@ -377,13 +381,26 @@ function faces(c::ReferenceSimplex{2})
         simplex(p[1], p[2]))
 end
 
+function faces(ch::ReferenceSimplex{2}, ::Type{Val{0}})
+    simplex.(vertices(ch))
+end
 
 barytocart(ch::ReferenceSimplex, u) = u # barytocart(ch.simplex, u)
 carttobary(ch::ReferenceSimplex, p) = SVector{2}(p[1:2]) # carttobary(ch.simplex, p)
 
-domain(ch::Simplex{U,D,C,N,T}) where {U,D,C,T,N} = ReferenceSimplex{D,T,N}()
+function domain(ch::Simplex{U,D,C,N,T}) where {U,D,C,T,N} ReferenceSimplex{D,T,N}() end
+function domain(ch::ReferenceSimplex) ch end
+
 neighborhood(ch::ReferenceSimplex, u) = SVector(u)
 neighborhood(ch::ReferenceSimplex, u::AbstractVector) = SVector{length(u)}(u)
+
+function vertices(ch::ReferenceSimplex{1,T}) where {T}
+    SVector(
+        point(T,1),
+        point(T,0),
+    )
+end
+
 function vertices(ch::ReferenceSimplex{2,T}) where {T}
     SVector(
         point(T,1,0),

src/tpredicates.jl

@@ -8,6 +8,22 @@ function in(mesh::CompScienceMeshes.AbstractMesh)
     return f
 end
 
+function setminus(m1::AbstractMesh, m2::AbstractMesh)
+    m = submesh(!in(m2), m1)
+end
+
+@testitem "setminus" begin
+    m1 = meshrectangle(1.0, 1.0, 0.5, 3)
+    m2 = submesh(m1) do m,c
+        ch = chart(m,c)
+        x = cartesian(center(ch))
+        x[1] < 0.5
+    end
+    m3 = setminus(m1,m2)
+    @test length(m2) == 4
+    @test length(m3) == 4
+end
+
 """
 Creates a predicate that can be used to check wheter an edge is interior to
 a surface (true) or on its boundary (false). This predicate is based on combinatorics.