diff --git a/Project.toml b/Project.toml
index eb5766a7..66362225 100644
--- a/Project.toml
+++ b/Project.toml
@@ -23,7 +23,7 @@ Downloads = "1.4"
 FileIO = "1.6"
 FillArrays = "0.11, 0.13"
 GraphRecipes = "0.5"
-Gridap = "0.17"
+Gridap = "0.17.16"
 GridapEmbedded = "0.8"
 IterativeSolvers = "0.9"
 MeshIO = "0.4"
diff --git a/src/Embedded.jl b/src/Embedded.jl
index f0980723..1e9a654b 100644
--- a/src/Embedded.jl
+++ b/src/Embedded.jl
@@ -73,37 +73,38 @@ function _cut_stl(model::DiscreteModel,geom::STLGeometry;kwargs...)
   (bgcell_to_ioc,subcells,cell_to_io,subfacets,face_to_io,oid_to_ls),bgface_to_ioc
 end
 
+function send_to_ref_space(
+  model::DiscreteModel,
+  cell_to_bgcell::Vector,
+  subgrid::Grid)
+
+  send_to_ref_space(get_grid(model),cell_to_bgcell,subgrid)
+end
+
 function send_to_ref_space(grid::Grid,cell_to_bgcell::Vector,subgrid::Grid)
-  send_to_ref_space(
-    grid,
-    cell_to_bgcell,
-    get_cell_node_ids(subgrid),
-    get_node_coordinates(subgrid))
+  bgcell_map = get_cell_map(grid)
+  bgcell_invmap = lazy_map(inverse_map,bgcell_map)
+  cell_invmap = lazy_map(Reindex(bgcell_invmap),cell_to_bgcell)
+  cell_nodes = get_cell_node_ids(subgrid)
+  node_coords = get_node_coordinates(subgrid)
+  node_cell = _first_inverse_index_map(cell_nodes,num_nodes(subgrid))
+  node_invmap = lazy_map(Reindex(cell_invmap),node_cell)
+  node_rcoords = lazy_map(evaluate,node_invmap,node_coords)
+  _collect(node_rcoords)
 end
 
-function send_to_ref_space(
-  grid::Grid,
-  cell_to_bgcell::Vector,
-  cell_to_points::AbstractVector,
-  points::Vector)
-
-  rpoints = zero(points)
-  cache = array_cache(cell_to_points)
-  bgcache = array_cache(get_cell_node_ids(grid))
-  for i in 1:length(cell_to_points)
-    for p in getindex!(cache,cell_to_points,i)
-      rpoints[p] = send_to_ref_space!(bgcache,grid,cell_to_bgcell[i],points[p])
+function _first_inverse_index_map(a_to_b,nb)
+  na = length(a_to_b)
+  T = eltype(eltype(a_to_b))
+  b_to_a = ones(T,nb)
+  c = array_cache(a_to_b)
+  for a in 1:na
+    bs = getindex!(c,a_to_b,a)
+    for b in bs
+      b_to_a[b] = a
     end
   end
-  rpoints
-end
-
-function send_to_ref_space!(cache,grid::Grid,cell::Integer,point::Point)
-  cell_nodes = getindex!(cache,get_cell_node_ids(grid),cell)
-  node_coordinates = get_node_coordinates(grid)
-  pmin = node_coordinates[cell_nodes[1]]
-  pmax = node_coordinates[cell_nodes[end]]
-  Point(Tuple(point-pmin)./Tuple(pmax-pmin))
+  b_to_a
 end
 
 function _normals(geom::STLGeometry,face_to_stlface)
@@ -130,3 +131,13 @@ function check_requisites(geo::STLGeometry,bgmodel::DiscreteModel;kwargs...)
   stl = get_stl(geo)
   check_requisites(stl,bgmodel)
 end
+
+function _collect(a::LazyArray)
+  c = array_cache(a)
+  T = eltype(a)
+  b = zeros(T,size(a))
+  for i in eachindex(a)
+    b[i] = getindex!(c,a,i)
+  end
+  b
+end
diff --git a/src/STLs.jl b/src/STLs.jl
index 74c3efed..162f0589 100644
--- a/src/STLs.jl
+++ b/src/STLs.jl
@@ -223,6 +223,7 @@ function get_cell!(c,grid::Grid,i)
   p = get_polytope( only(get_reffes(grid)) )
   get_cell!(c,T,X,p,i)
 end
+
 function get_facet_cache(model::DiscreteModel)
   Dc = num_dims(model)
   get_dface_cache(model,Dc-1)
@@ -240,9 +241,11 @@ end
 
 function get_dface!(c,model::DiscreteModel,i,::Val{d}) where d
   topo = get_grid_topology(model)
+  Dc = num_dims(model)
+  PT = ExtrusionPolytope{Dc}
   T = get_face_vertices(topo,d)
   X = get_vertex_coordinates(topo)
-  p = get_polytope( only(get_reffes(model)) )
+  p = get_polytope( only(get_reffes(model)) )::PT
   get_dface!(c,T,X,p,i,Val{d}())
 end
 
@@ -867,6 +870,18 @@ end
 
 min_height(model::DiscreteModel) = min_height(get_grid(model))
 
+function max_length(grid::Grid)
+  max_len = 0.0
+  c = get_cell_cache(grid)
+  for i in 1:num_cells(grid)
+    cell = get_cell!(c,grid,i)
+    max_len = max(max_len,max_length(cell))
+  end
+  max_len
+end
+
+max_length(model::DiscreteModel) = max_length(get_grid(model))
+
 # Cartesian Grid
 
 function get_bounding_box(grid::CartesianGrid)
@@ -886,8 +901,13 @@ function measure(a::CartesianGrid{D,T,typeof(identity)},mask) where {D,T}
   cell_vol * count(mask)
 end
 
+function compute_cell_to_facets(model_a::DiscreteModel,grid_b)
+  grid_a = get_grid(model_a)
+  compute_cell_to_facets(grid_a,grid_b)
+end
 
-function compute_cell_to_facets(grid::CartesianGrid,stl::STL)
+
+function compute_cell_to_facets(grid::CartesianGrid,stl)
   CELL_EXPANSION_FACTOR = 1e-3
   desc = get_cartesian_descriptor(grid)
   @assert length(get_reffes(grid)) == 1
@@ -920,13 +940,112 @@ function compute_cell_to_facets(grid::CartesianGrid,stl::STL)
       end
     end
   end
-  cell_to_stl_facets = [ Int32[] for _ in 1:num_cells(grid) ]
-  for (cells,stl_facets) in zip(thread_to_cells,thread_to_stl_facets)
-    for (cell,stl_facet) in zip(cells,stl_facets)
-      push!(cell_to_stl_facets[cell],stl_facet)
+  ncells = num_cells(grid)
+  assemble_threaded_sparse_map(thread_to_cells,thread_to_stl_facets,ncells)
+end
+
+function compute_cell_to_facets(a::UnstructuredGrid,b)
+  tmp = cartesian_bounding_model(a)
+  tmp_to_a = compute_cell_to_facets(tmp,a)
+  tmp_to_b = compute_cell_to_facets(tmp,b)
+  a_to_tmp = inverse_index_map(tmp_to_a,num_cells(a))
+  a_to_b = compose_index_map(a_to_tmp,tmp_to_b)
+  a_to_b = filter_intersections(a,b,a_to_b)
+  a_to_b
+end
+
+function filter_intersections(a::Grid,b,a_to_b)
+  CELL_EXPANSION_FACTOR = 1e-3
+  n = Threads.nthreads()
+  t_as = map(_->Int32[],1:n)
+  t_bs = map(_->Int32[],1:n)
+  t_c = map(_->(get_cell_cache(a),get_cell_cache(b),array_cache(a_to_b)),1:n)
+  Threads.@threads for ai in 1:num_cells(a)
+    i = Threads.threadid()
+    ca,cb,c = t_c[i]
+    as = t_as[i]
+    bs = t_bs[i]
+    cell_a = get_cell!(ca,a,ai)
+    pmin,pmax = get_bounding_box(cell_a)
+    Δ = norm(pmin-pmax) * CELL_EXPANSION_FACTOR
+    cell_a = expand_face(cell_a,Δ)
+    for bi in getindex!(c,a_to_b,ai)
+      cell_b = get_cell!(cb,b,bi)
+      if has_intersection(cell_a,cell_b)
+        push!(as,ai)
+        push!(bs,bi)
+      end
+    end
+  end
+  assemble_threaded_sparse_map(t_as,t_bs,num_cells(a))
+end
+
+function cartesian_bounding_model(grid::Grid)
+  h = max_length(grid)
+  pmin,pmax = get_bounding_box(grid)
+  s = Tuple(pmax-pmin)
+  partition = Int.( ceil.( s ./ h  ))
+  CartesianDiscreteModel(pmin,pmax,partition)
+end
+
+function inverse_index_map(
+  a_to_b::AbstractVector{<:AbstractVector},
+  nb=maximum(maximum(a_to_b)))
+
+  na = length(a_to_b)
+  as = Int32[]
+  bs = Int32[]
+  cache = array_cache(a_to_b)
+  for a in 1:na
+    for b in getindex!(cache,a_to_b,a)
+      push!(as,a)
+      push!(bs,b)
+    end
+  end
+  assemble_sparse_map(bs,as,nb)
+end
+
+function compose_index_map(
+  a_to_b::AbstractVector{<:AbstractVector},
+  b_to_c::AbstractVector{<:AbstractVector})
+
+  na = length(a_to_b)
+  as = Int32[]
+  cs = Int32[]
+  cache_ab = array_cache(a_to_b)
+  cache_bc = array_cache(b_to_c)
+  for a in 1:na
+    for b in getindex!(cache_ab,a_to_b,a)
+      for c in getindex!(cache_bc,b_to_c,b)
+        push!(as,a)
+        push!(cs,c)
+      end
+    end
+  end
+  assemble_sparse_map(as,cs,na)
+end
+
+function assemble_threaded_sparse_map(tas,tbs,na=maximum(maximum,tas))
+  a_to_b = map(_->Int32[],1:na)
+  for (as,bs) in zip(tas,tbs)
+    assemble_sparse_map!(a_to_b,as,bs)
+  end
+  a_to_b
+end
+
+function assemble_sparse_map(as,bs,na=maximum(as))
+  a_to_b = map(_->Int32[],1:na)
+  assemble_sparse_map!(a_to_b,as,bs)
+  a_to_b
+end
+
+function assemble_sparse_map!(a_to_b,as,bs)
+  for (a,b) in zip(as,bs)
+    if b ∉ a_to_b[a]
+      push!(a_to_b[a],b)
     end
   end
-  cell_to_stl_facets
+  a_to_b
 end
 
 function get_cells_around(desc::CartesianDescriptor{D},pmin::Point,pmax::Point) where D
@@ -934,7 +1053,7 @@ function get_cells_around(desc::CartesianDescriptor{D},pmin::Point,pmax::Point)
   _,cmax = get_cell_bounds(desc,pmax)
   cmin = CartesianIndices(desc.partition)[cmin]
   cmax = CartesianIndices(desc.partition)[cmax]
-  ranges = ntuple( i -> cmin.I[i] : cmax.I[i], Val{D}() )
+  ranges = ntuple( i -> cmin.I[i]::Int : cmax.I[i]::Int, Val{D}() )
   CartesianIndices( ranges )
 end
 
diff --git a/src/SimplexFaces.jl b/src/SimplexFaces.jl
index 00973b50..f0280372 100644
--- a/src/SimplexFaces.jl
+++ b/src/SimplexFaces.jl
@@ -149,6 +149,11 @@ is_n_cube(a::Face) = is_n_cube(get_polytope(a))
 
 simplex_face(v::Point...) = simplex_face(v)
 
+function simplex_face(coords::AbstractVector{<:Point{D}}) where D
+  tcoords = ntuple( i -> coords[i], Val{D+1}() )
+  simplex_face(tcoords)
+end
+
 @inline function simplex_face(v::NTuple{N,<:Point}) where N
   if N == 1
     v[1]
@@ -303,6 +308,15 @@ function min_height(f::Face{Df}) where Df
   min_dist
 end
 
+function max_length(f::Face{Df}) where Df
+  max_len = 0.0
+  for e in 1:num_edges(f)
+    edge = get_edge(f,e)
+    max_len = max(max_len,length(edge))
+  end
+  max_len
+end
+
 get_bounding_box(p::Point) = p,p
 
 function get_bounding_box(f::Face)
@@ -435,6 +449,118 @@ function contains_projection(f::Face,p::Point)
   true
 end
 
+function intersection_point(f::Face{1,D},p::AbstractPlane{D}) where D
+  v1 = f[1]
+  v2 = f[2]
+  d1 = signed_distance(v1,p)
+  d2 = signed_distance(v2,p)
+  if d1-d2 == 0
+    d1,d2 = 1,10
+  end
+  (v1*d2-v2*d1)/(d2-d1)
+end
+
+intersection_point(a::Face{1},b::Face) = _intersection_point(a,b)
+
+intersection_point(a::Face{1},b::Face{1}) = _intersection_point(a,b)
+
+intersection_point(a::Face,b::Face{1}) = _intersection_point(b,a)
+
+function _intersection_point(e::Face{1},f::Face)
+  @notimplementedif num_dims(f) + 1 != num_point_dims(f)
+  n = normal(f)
+  c = center(f)
+  plane = Plane(c,n)
+  intersection_point(e,plane)
+end
+
+function expand_face(f::Face,dist::Real)
+  Dc = num_dims(f)
+  if is_simplex(f)
+    verts = ntuple(i->_pull_vertex(f,i,dist),Val{Dc+1}())
+    simplex_face( verts )
+  else
+    @notimplemented
+  end
+end
+
+function _pull_vertex(f::Face,v::Integer,d::Real)
+  dir = _vertex_direction(f,v)
+  f[v] + dir*d
+end
+
+function _vertex_direction(f::Face,v::Integer)
+  p = get_polytope(f)
+  v_to_e = get_faces(p,0,1)
+  e_to_v = get_faces(p,1,0)
+  dir = zero(f[1])
+  n = length(v_to_e[v])
+  for e in v_to_e[v]
+    for vneig in e_to_v[e]
+      if vneig != v
+        dv = f[v] - f[vneig]
+        dir += dv / n
+      end
+    end
+  end
+  @assert norm(dir) > 0
+  dir / norm(dir)
+end
+
+
+# General predicates
+
+function has_intersection(a::Face,b::Face)
+  Dp = num_point_dims(a)
+  pa = get_polytope(a)
+  pb = get_polytope(b)
+  for da in 0:num_dims(a)
+    db = Dp-da
+    if num_dims(b) ≥ db
+      for fa in 1:num_faces(pa,da), fb in 1:num_faces(pb,db)
+        if has_intersection_point(a,da,fa,b,db,fb)
+          return true
+        end
+      end
+    end
+  end
+  false
+end
+
+function has_intersection_point(p::Point,f::Face{D,D}) where D
+  contains_projection(f,p)
+end
+
+has_intersection_point(f::Face,p::Point) = has_intersection_point(p,f)
+
+function has_intersection_point(a::Face,b::Face)
+  p = intersection_point(a,b)
+  contains_projection(a,p) && contains_projection(b,p)
+end
+
+function has_intersection_point(
+  f::Face,
+  d::Integer,
+  df::Integer,
+  args::Vararg{<:Any,N}) where N
+
+  if d == 0
+    f0 = get_dface(f,df,Val{0}())
+    has_intersection_point(args...,f0)
+  elseif d == 1
+    f1 = get_dface(f,df,Val{1}())
+    has_intersection_point(args...,f1)
+  elseif d == 2
+    f2 = get_dface(f,df,Val{2}())
+    has_intersection_point(args...,f2)
+  elseif d == 3
+    f3 = get_dface(f,df,Val{3}())
+    has_intersection_point(args...,f3)
+  else
+    @notimplemented
+  end
+end
+
 # Voxel predicates
 
 function voxel_intersection(f::Face{1,2},pmin::Point,pmax::Point,p::Polytope)
@@ -510,6 +636,26 @@ function voxel_intersection(e::Face{1,D},pmin::Point{D},pmax::Point{D}) where D
   t_min < t_max
 end
 
+function voxel_intersection(
+  c::Face{D,D},
+  pmin::Point{D},
+  pmax::Point{D},
+  p::Polytope) where D
+
+  @check all( pmin .≤ pmax ) "pmin  > pmax"
+  for i in 1:num_facets(c)
+    f = get_facet(c,i)
+    voxel_intersection(f,pmin,pmax,p) && return true
+  end
+  for i in 1:num_facets(c)
+    f = get_facet(c,i)
+    plane = Plane( center(f), normal(c,i) )
+    v_to_dists = _compute_distances(plane,pmin,pmax,p)
+    all( d -> d < 0, v_to_dists) || return false
+  end
+  true
+end
+
 function _compute_distances(
   Π::Plane,
   pmin::Point,
@@ -579,6 +725,69 @@ function get_cell_planes(p::Polytope,pmin::Point,pmax::Point)
   Π_cell, Π_ids, Π_inout
 end
 
+function get_cell_planes(p::Polytope,coords)
+  if is_simplex(p)
+    f = simplex_face(coords)
+    get_cell_planes(f)
+  elseif is_n_cube(p)
+    pmin,pmax = coords[1],coords[end]
+    get_cell_planes(p,pmin,pmax)
+  end
+end
+
+function get_cell_planes(face::Face)
+  n_facets = num_facets(face)
+  facets = lazy_map(get_facet,Fill(face,n_facets),1:n_facets)
+  centers = lazy_map(center,facets)
+  normals = lazy_map(normal,Fill(face,n_facets),1:n_facets)
+  mask = lazy_map(n-> n > -n, normals)
+  facs = lazy_map(i->(-1)^(i+1),mask)
+  _normals = lazy_map(*,normals,facs)
+  planes = lazy_map(Plane,centers,_normals)
+  ids = -(1:n_facets)
+  planes,ids,mask
+end
+
+function displace(plane::Plane,dist,oriented=true)
+  c0 = center(plane)
+  n0 = normal(plane)
+  n = oriented ? n0 : -n0
+  c = c0 + n*dist
+  Plane(c,n0)
+end
+
+function displace(plane::CartesianPlane,dist,oriented=true)
+  c0 = origin(plane)
+  d = plane.d
+  or = plane.positive
+  diff = or == oriented ? dist : -dist
+  c = c0 + diff
+  dir = or ? 1 : -1
+  CartesianPlane(c,d,dir)
+end
+
+function expand_planes(
+  planes::AbstractVector{<:AbstractPlane},
+  inout::AbstractVector,
+  dist)
+
+  n_planes = length(planes)
+  lazy_map(displace,planes,Fill(dist,n_planes),inout)
+end
+
+function _visualization_grid(face::Face)
+  X = collect(face.vertices)
+  T = Table( [1:length(X)] )
+  ctype = [1]
+  reffes = [LagrangianRefFE(get_polytope(face))]
+  UnstructuredGrid(X,T,reffes,ctype)
+end
+
+function writevtk(face::Face,filename;kwargs...)
+  grid = _visualization_grid(face)
+  writevtk(grid,filename;kwargs...)
+end
+
 ## Orthogonal
 
 function orthogonal(a::VectorValue{2})
diff --git a/src/SubTriangulations.jl b/src/SubTriangulations.jl
index 1b4919c3..4c356f62 100644
--- a/src/SubTriangulations.jl
+++ b/src/SubTriangulations.jl
@@ -16,7 +16,7 @@ struct SubtriangulationLabels
 end
 
 function subtriangulate(
-  bgmodel::CartesianDiscreteModel,
+  bgmodel::DiscreteModel,
   stlmodel::DiscreteModel;
   threading=:spawn,
   kdtree=false,
@@ -89,7 +89,7 @@ function subtriangulate(
   cell_grid, face_grid, labels
 end
 
-function subtriangulate(bgmodel::CartesianDiscreteModel,args...;kwargs...)
+function subtriangulate(bgmodel::DiscreteModel,args...;kwargs...)
   stlmodel = _to_stl_model(args...)
   subtriangulate(bgmodel,stlmodel;kwargs...)
 end
@@ -119,14 +119,12 @@ function compute_polyhedra!(caches,Γ0,stl,p,f_to_isempty,Πf,Πr,
   c_to_stlf,node_to_coords,cell_to_nodes,cell;atol,kdtree)
 
   cell_coords = _get_cell_coordinates!(caches,node_to_coords,cell_to_nodes,cell)
-  pmin = cell_coords[1]
-  pmax = cell_coords[end]
 
    _faces = _get_cell_faces!(caches,stl,c_to_stlf,cell,f_to_isempty)
   facets,empty_facets,reflex_faces = _faces
 
-  Πk,Πk_ids,Πk_io = get_cell_planes(p,pmin,pmax)
-  Πk⁺,Πk⁺_ids,Πk⁺_io = get_cell_planes(p,pmin-atol,pmax+atol)
+  Πk,Πk_ids,Πk_io = get_cell_planes(p,cell_coords)
+  Πk⁺,Πk⁺_ids,Πk⁺_io = expand_planes(Πk,Πk_io,atol), Πk_ids, Πk_io
   Πk⁺_ids = Πk⁺_ids .+ Πk_ids[end]
 
   Πkf,Πkf_ids = filter_face_planes(stl,Πr,reflex_faces,Πf,facets)
@@ -274,8 +272,12 @@ function delete_small_subfacets!(bgmodel,T,X,arrays...)
   delete_small_subfaces!(bgmodel,T,X,TRI,arrays...)
 end
 
+function max_length(model::CartesianDiscreteModel)
+  float(get_cartesian_descriptor(model).sizes[1])
+end
+
 function delete_small_subfaces!(bgmodel,T,X,p::Polytope{D},arrays...) where D
-  h = float(get_cartesian_descriptor(bgmodel).sizes[1])
+  h = max_length(bgmodel)
   c = array_cache(T)
   ids = findall( i -> measure(get_cell!(c,T,X,p,i)) < eps(h^D), 1:length(T) )
   deleteat!(T,ids)
@@ -905,18 +907,9 @@ function complete_nodes_to_inout!(node_to_inout,polys,inout,p::Polytope)
     for v in 1:num_vertices(poly)
       isactive(poly,v) || continue
     #  v = inout == FACE_CUT ? v : v_to_v[v]
-      if count(Π -> -2*D ≤ Π < 0 ,v_to_Π[v]) == D
-        i = 0
-        node = 0
-        for _ in 1:D
-          i = findnext(Π -> Π < 0,v_to_Π[v],i+1)
-          f = -v_to_Π[v][i]
-          d = D - ((f-1)>>1)
-          ud = iseven(f)
-          node |= ud<<(d-1)
-        end
-        node += 1
-        #@assert  node_to_inout[node] ∈ (inout,UNSET)
+      node = find_polytope_vertex(v_to_Π[v],p)
+      if node != UNSET
+      #@assert  node_to_inout[node] ∈ (inout,UNSET)
         _inout = inout
         if node_to_inout[node] ∉ (inout,UNSET)
           _inout = FACE_CUT
@@ -928,6 +921,52 @@ function complete_nodes_to_inout!(node_to_inout,polys,inout,p::Polytope)
   node_to_inout
 end
 
+function find_polytope_vertex(planes,p::Polytope)
+  if is_n_cube(p)
+    _find_n_cube_vertex(planes,p)
+  else
+    _find_polytope_vertex(planes,p)
+  end
+end
+
+function _find_polytope_vertex(planes,p::Polytope)
+  D = num_dims(p)
+  n_facets = num_faces(p,D-1)
+  if count(Π -> -n_facets ≤ Π < 0 ,planes) == D
+    f_to_v = get_faces(p,D-1,0)
+    v_to_f = get_faces(p,0,D-1)
+    i = findfirst(Π -> Π < 0,planes)
+    f = -planes[i]
+    for v in f_to_v[f]
+      if all(f->(-f) ∈ planes,v_to_f[v])
+        return v
+      end
+    end
+  end
+  UNSET
+end
+
+
+function _find_n_cube_vertex(planes,p::Polytope)
+  @notimplementedif !is_n_cube(p)
+  D = num_dims(p)
+  if count(Π -> -2*D ≤ Π < 0 , planes) == D
+    i = 0
+    node = 0
+    for _ in 1:D
+      i = findnext(Π -> Π < 0, planes,i+1)
+      f = -planes[i]
+      d = D - ((f-1)>>1)
+      ud = iseven(f)
+      node |= ud<<(d-1)
+    end
+    node += 1
+  else
+    node = UNSET
+  end
+  node
+end
+
 function complete_facets_to_inoutcut!(facet_to_inoutcut,polys,inout,p::Polytope)
   facet_list = Int32[]
   for poly in polys
diff --git a/test/Integration.jl b/test/Integration.jl
index 7da373b1..3506f0e6 100644
--- a/test/Integration.jl
+++ b/test/Integration.jl
@@ -7,7 +7,7 @@ using GridapEmbedded
 using Test
 
 using STLCutters: compute_stl_model
-using STLCutters: read_stl, merge_nodes, get_bounding_box 
+using STLCutters: read_stl, merge_nodes, get_bounding_box
 
 # Manufactured solution
 u0(x) = x[1] + x[2] - x[3]
@@ -69,4 +69,55 @@ a = sum( ∫( ∇(v)⋅∇(u) ) * dΩ )
 b = sum( ∫( v⋅n_Γd⋅∇(u) ) * dΓd )
 @test abs( a-b ) < 1e-9
 
+
+# Simplex background
+#
+
+bgmodel = CartesianDiscreteModel(pmin,pmax,partition)
+bgmodel = simplexify(bgmodel,positive=true)
+
+# Cut the background model
+
+cutgeo, = cut(bgmodel,geo)
+
+# Setup integration meshes
+Ω = Triangulation(cutgeo,PHYSICAL)
+Γd = EmbeddedBoundary(cutgeo)
+Γg = GhostSkeleton(cutgeo)
+
+# Setup normal vectors
+n_Γd = get_normal_vector(Γd)
+n_Γg = get_normal_vector(Γg)
+
+#writevtk(Ω,"trian_O")
+#writevtk(Γd,"trian_Gd",cellfields=["normal"=>n_Γd])
+#writevtk(Γg,"trian_Gg",cellfields=["normal"=>n_Γg])
+#writevtk(Triangulation(bgmodel),"bgtrian")
+
+# Setup Lebesgue measures
+order = 1
+degree = 2*order
+dΩ = Measure(Ω,degree)
+dΓd = Measure(Γd,degree)
+dΓg = Measure(Γg,degree)
+
+vol = sum( ∫(1)*dΩ  )
+surf = sum( ∫(1)*dΓd )
+
+@test vol ≈ 1
+@test surf ≈ 6
+
+# Setup FESpace
+Ω_act = Triangulation(cutgeo,ACTIVE,geo)
+V = TestFESpace(Ω_act,ReferenceFE(lagrangian,Float64,order),conformity=:H1)
+v = FEFunction(V,rand(num_free_dofs(V)))
+u = interpolate(u0,V)
+
+# Check divergence theorem
+a = sum( ∫( ∇(v)⋅∇(u) ) * dΩ )
+b = sum( ∫( v⋅n_Γd⋅∇(u) ) * dΓd )
+@test abs( a-b ) < 1e-9
+
+
+
 end
diff --git a/test/SampleGeometriesTests.jl b/test/SampleGeometriesTests.jl
index 293bea93..1f54d72d 100644
--- a/test/SampleGeometriesTests.jl
+++ b/test/SampleGeometriesTests.jl
@@ -13,7 +13,7 @@ using STLCutters: compute_cartesian_descriptor
 download = download_thingi10k
 
 function main(filename;
-  nmin=10,nmax=100,δ=0.2,tolfactor=1000,kdtree=false,output=nothing)
+  nmin=10,nmax=100,δ=0.2,tolfactor=1000,kdtree=false,simplex=false,output=nothing)
 
   println("Running: $(basename(filename))")
 
@@ -23,6 +23,9 @@ function main(filename;
   Δ = (pmax-pmin)*δ
   desc = compute_cartesian_descriptor(pmin-Δ,pmax+Δ;nmin,nmax)
   model = CartesianDiscreteModel(desc)
+  if simplex
+    model = simplexify(model,positive=true)
+  end
 
   @test check_requisites(stl,model)
 
@@ -72,6 +75,7 @@ end
 filename = joinpath(@__DIR__,"../test/data/47076.stl")
 main(filename,nmax=50)
 main(filename,nmax=10,nmin=10,kdtree=true)
+main(filename,nmax=50,nmin=10,simplex=true)
 
 filename = joinpath(@__DIR__,"../test/data/47076_sf.obj")
 main(filename,nmax=50)
diff --git a/test/SimplexFacesTests.jl b/test/SimplexFacesTests.jl
index 708540c4..76f23a58 100644
--- a/test/SimplexFacesTests.jl
+++ b/test/SimplexFacesTests.jl
@@ -13,9 +13,12 @@ using STLCutters: normal
 using STLCutters: origin
 using STLCutters: signed_distance
 using STLCutters: contains_projection
+using STLCutters: expand_face
+using STLCutters: has_intersection
 using STLCutters: voxel_intersection
 using STLCutters: simplex_face
 using STLCutters: min_height
+using STLCutters: displace
 using STLCutters: Plane
 using STLCutters: CartesianPlane
 
@@ -277,6 +280,81 @@ b = Point(1.0,1.0,1.0)
 @test contains_projection(f,a)
 @test !contains_projection(f,b)
 
+# Displace planes
+
+x = Point(1.0,2.0,3.0)
+n = VectorValue(0,0,1.0)
+
+plane0 = Plane(x,n)
+dist = 0.5
+plane = displace(plane0,dist)
+@test normal(plane) == normal(plane0) == n
+@test origin(plane) == Point(1,2,3.5)
+
+plane = displace(plane0,dist,false)
+@test normal(plane) == normal(plane0) == n
+@test origin(plane) == Point(1,2,2.5)
+
+x = Point(1,2,3)
+d = 3
+plane0 = CartesianPlane(x,d,1)
+plane = displace(plane0,dist)
+dist = 0.5
+@test origin(plane)⋅normal(plane) == origin(plane0)⋅normal(plane0)+dist == 3.5
+
+x = Point(1,2,3)
+d = 3
+plane0 = CartesianPlane(x,d,-1)
+plane = displace(plane0,dist)
+dist = 0.5
+@test origin(plane)⋅normal(plane) == origin(plane0)⋅normal(plane0)+dist == -2.5
+
+x = Point(1,2,3)
+d = 3
+plane0 = CartesianPlane(x,d,1)
+plane = displace(plane0,dist,false)
+dist = 0.5
+@test origin(plane)⋅normal(plane) == origin(plane0)⋅normal(plane0)-dist == 2.5
+
+x = Point(1,2,3)
+d = 3
+plane0 = CartesianPlane(x,d,-1)
+plane = displace(plane0,dist,false)
+dist = 0.5
+@test origin(plane)⋅normal(plane) == origin(plane0)⋅normal(plane0)-dist == -3.5
+
+# General intersection
+
+p1 = Point(0.0,0.0,0.0)
+p2 = Point(1.0,0.0,0.0)
+p3 = Point(0.0,1.0,0.0)
+p4 = Point(0.0,0.0,1.0)
+f1 = simplex_face(p1,p2,p3,p4)
+
+offset = 0.5
+p1 = Point(0.0,0.0,0.0) + offset
+p2 = Point(1.0,0.0,0.0) + offset
+p3 = Point(0.0,1.0,0.0) + offset
+p4 = Point(0.0,0.0,1.0) + offset
+f2 = simplex_face(p1,p2,p3,p4)
+
+f3 = expand_face(f1,0.5)
+f4 = expand_face(f2,0.5)
+
+@test !has_intersection(f1,f2)
+@test has_intersection(f1,f3)
+@test has_intersection(f3,f1)
+@test !has_intersection(f2,f3)
+@test has_intersection(f1,f4)
+@test has_intersection(f2,f4)
+@test has_intersection(f3,f4)
+@test has_intersection(f4,f3)
+
+#writevtk(f1,"f1")
+#writevtk(f2,"f2")
+#writevtk(f3,"f3")
+#writevtk(f4,"f4")
+
 # Voxel intersection Face{2}
 
 p = QUAD
@@ -351,5 +429,38 @@ e = simplex_face(p1,p2,p3)
 @test !voxel_intersection(d,pmin,pmax,p)
 @test !voxel_intersection(e,pmin,pmax,p)
 
+# Voxel intersection Cell
+
+t = simplex_face(get_vertex_coordinates(TET)...)
+pmin,pmax = get_bounding_box(t)
+p0 = center(t)
+
+@test voxel_intersection(t,pmin,pmax,HEX)
+@test voxel_intersection(t,pmin,pmin+0.1,HEX)
+@test !voxel_intersection(t,pmax,pmax+0.1,HEX)
+@test voxel_intersection(t,p0,p0+0.1,HEX)
+
+f, = writevtk(t,"tet")
+rm.(f)
+
+
+pmin = Point(0.0,0.0,0.0)
+pmax = Point(12.0,12.0,12.0)
+
+tverts =
+( Point(19.0, 0.0, 1.0),
+  Point(17.0, 5.0, -1.0),
+  Point(10.0, 5.0, -1.0),
+  Point(10.0, 5.0, 7.0))
+t = STLCutters.simplex_face(tverts)
+
+@test voxel_intersection(t,pmin,pmax,HEX)
+
+vverts = map(i->STLCutters._get_vertex(HEX,pmin,pmax,i),1:2^3)
+v = STLCutters.Polyhedron(HEX,vverts)
+writevtk(t,"t")
+writevtk(v,"v")
+
+
 end # module
 
diff --git a/test/SubTriangulationsTests.jl b/test/SubTriangulationsTests.jl
index 835ac2aa..2eda3619 100644
--- a/test/SubTriangulationsTests.jl
+++ b/test/SubTriangulationsTests.jl
@@ -258,6 +258,101 @@ println("Num subcells: $(num_cells(subcells))")
 @test in_volume + out_volume ≈ volume(grid)
 @test in_volume ≈ 74.12595970214333
 
+# Simplex background grid
+
+X,T,N = read_stl(joinpath(@__DIR__,"data/Bunny-LowPoly.stl"))
+stl = compute_stl_model(T,X)
+stl = merge_nodes(stl)
+
+#writevtk(get_grid(stl),"stl")
+
+p = HEX
+δ = 0.2
+n = 10
+D = 3
+
+pmin,pmax = get_bounding_box(stl)
+Δ = (pmax-pmin)*δ
+pmin = pmin - Δ
+pmax = pmax + Δ
+partition = (n,n,n)
+
+model = CartesianDiscreteModel(pmin,pmax,partition)
+model = simplexify(model,positive=true)
+grid = get_grid(model)
+
+@time subcells, subfaces, labels = subtriangulate(model,stl,kdtree=false)
+
+#writevtk(submesh,"submesh",cellfields=["inout"=>k_to_io,"bgcell"=>k_to_bgcell])
+#writevtk(grid,"bgmesh",cellfields=["inoutcut"=>bgcell_to_ioc])
+
+bgmesh_in_vol = volume(grid,labels.bgcell_to_ioc,:in)
+bgmesh_out_vol  = volume(grid,labels.bgcell_to_ioc,:out)
+bgmesh_cut_vol  = volume(grid,labels.bgcell_to_ioc,:cut)
+submesh_in_vol  = volume(subcells,labels.cell_to_io,:in)
+submesh_out_vol = volume(subcells,labels.cell_to_io,:out)
+
+in_volume = bgmesh_in_vol + submesh_in_vol
+out_volume = bgmesh_out_vol + submesh_out_vol
+cut_volume = bgmesh_cut_vol
+
+
+#celldata = [ "inoutcut" => labels.bgcell_to_ioc ]
+#writevtk( grid, "bgcells"; celldata )
+#
+#celldata = [ "inout" => labels.cell_to_io, "bgcell" => labels.cell_to_bgcell ]
+#writevtk( subcells, "subcells"; celldata )
+#
+#celldata = [ "bgcell" => labels.face_to_bgcell ]
+#writevtk( subfaces, "subfaces" )
+
+println("Num subcells: $(num_cells(subcells))")
+@test surface(get_grid(stl)) ≈ surface(subfaces)
+@test submesh_in_vol + submesh_out_vol ≈ cut_volume
+@test in_volume + out_volume ≈ volume(grid)
+@test in_volume ≈ 273280.03374196636
+
+X,T,N = read_stl(joinpath(@__DIR__,"data/wine_glass.stl"))
+stl = compute_stl_model(T,X)
+stl = merge_nodes(stl)
+
+#writevtk(get_grid(stl),"stl")
+
+p = HEX
+δ = 0.2
+n = 20
+D = 3
+
+pmin,pmax = get_bounding_box(stl)
+Δ = (pmax-pmin)*δ
+pmin = pmin - Δ
+pmax = pmax + Δ
+partition = (n,n,n)
+
+model = CartesianDiscreteModel(pmin,pmax,partition)
+model = simplexify(model,positive=true)
+grid = get_grid(model)
+
+@time subcells, subfaces, labels = subtriangulate(model,stl,kdtree=false)
+
+#writevtk(submesh,"submesh",cellfields=["inout"=>k_to_io,"bgcell"=>k_to_bgcell])
+#writevtk(grid,"bgmesh",cellfields=["inoutcut"=>bgcell_to_ioc])
+
+bgmesh_in_vol = volume(grid,labels.bgcell_to_ioc,:in)
+bgmesh_out_vol  = volume(grid,labels.bgcell_to_ioc,:out)
+bgmesh_cut_vol  = volume(grid,labels.bgcell_to_ioc,:cut)
+submesh_in_vol  = volume(subcells,labels.cell_to_io,:in)
+submesh_out_vol = volume(subcells,labels.cell_to_io,:out)
+
+in_volume = bgmesh_in_vol + submesh_in_vol
+out_volume = bgmesh_out_vol + submesh_out_vol
+cut_volume = bgmesh_cut_vol
+
+println("Num subcells: $(num_cells(subcells))")
+@test surface(get_grid(stl)) ≈ surface(subfaces)
+@test submesh_in_vol + submesh_out_vol ≈ cut_volume
+@test in_volume + out_volume ≈ volume(grid)
+@test in_volume ≈ 74.12595970214333
 ## Kd-Tree
 
 X,T,N = read_stl(joinpath(@__DIR__,"data/Bunny-LowPoly.stl"))