diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index b047a8d..f316fd1 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -70,6 +70,41 @@ jobs:
       - run: |
           julia --color=yes --project=. --check-bounds=yes --depwarn=error -e '
             (1,) .== 1; include("test/GridapEmbeddedTests/runtests.jl")'
+  algoimutils:
+    name: AlgoimUtils ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        version:
+          - '1.8'
+        os:
+          - ubuntu-latest
+        arch:
+          - x64
+    steps:
+      - uses: actions/checkout@v2
+      - uses: julia-actions/setup-julia@v1
+        with:
+          version: ${{ matrix.version }}
+          arch: ${{ matrix.arch }}
+      - uses: actions/cache@v1
+        env:
+          cache-name: cache-artifacts
+        with:
+          path: ~/.julia/artifacts
+          key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
+          restore-keys: |
+            ${{ runner.os }}-test-${{ env.cache-name }}-
+            ${{ runner.os }}-test-
+            ${{ runner.os }}-
+      - uses: julia-actions/julia-buildpkg@v1
+      - run: |
+          julia --color=yes --project=. --check-bounds=yes --depwarn=error -e '
+            using Pkg; Pkg.instantiate()'
+      - run: |
+          julia --color=yes --project=. --check-bounds=yes --depwarn=error -e '
+            (1,) .== 1; include("test/AlgoimUtilsTests/runtests.jl")'
   docs:
     name: Documentation
     runs-on: ubuntu-latest
diff --git a/LICENSE.md b/LICENSE.md
index bc8fbce..652d476 100644
--- a/LICENSE.md
+++ b/LICENSE.md
@@ -1,7 +1,7 @@
 GridapEmbedded.jl Copyright and License
 ==
 
-Copyright (c) 2020 by [Francesc Verdugo](mailto:fverdugo@cimne.upc.edu) and [Santiago Badia](mailto:santiago.badia@monash.edu).
+Copyright (c) 2020 by [Francesc Verdugo](mailto:f.verdugo.rojano@vu.nl), [Eric Neiva](mailto:eric.neiva@college-de-france.fr) and [Santiago Badia](mailto:santiago.badia@monash.edu).
 
 GridapEmbedded.jl is licensed under the MIT Expat License
 
@@ -28,19 +28,33 @@ SOFTWARE.
 Citing `Gridap`
 ==
 
-In order to give credit to the `Gridap` contributors, we simply ask you to cite the refence below in any publication in which you have made use of `Gridap` packages.
-
-F. Verdugo and S. Badia. *A user-guide to Gridap -- grid-based approximation of partial differential equations in Julia*, 2019. [arXiv:1910.01412](https://arxiv.org/abs/1910.01412)
+In order to give credit to the `Gridap` contributors, we simply ask you to cite the references below in any publication in which you have made use of `Gridap` packages.
 
 ```
-@article{gridap_guide_2019,
-    author={Francesc Verdugo and Santiago Badia},
-    journal = {{arXiv}},
-    title = {{A user-guide to Gridap -- grid-based approximation of partial differential equations in Julia}},
-    year = {2019},
-    eprint={1910.01412},
-    archivePrefix={arXiv},
-    primaryClass={cs.MS},
+@article{Badia2020,
+  doi = {10.21105/joss.02520},
+  url = {https://doi.org/10.21105/joss.02520},
+  year = {2020},
+  publisher = {The Open Journal},
+  volume = {5},
+  number = {52},
+  pages = {2520},
+  author = {Santiago Badia and Francesc Verdugo},
+  title = {Gridap: An extensible Finite Element toolbox in Julia},
+  journal = {Journal of Open Source Software}
+}
+
+@article{Verdugo2022,
+  doi = {10.1016/j.cpc.2022.108341},
+  url = {https://doi.org/10.1016/j.cpc.2022.108341},
+  year = {2022},
+  month = jul,
+  publisher = {Elsevier {BV}},
+  volume = {276},
+  pages = {108341},
+  author = {Francesc Verdugo and Santiago Badia},
+  title = {The software design of Gridap: A Finite Element package based on the Julia {JIT} compiler},
+  journal = {Computer Physics Communications}
 }
 ```
 
@@ -92,4 +106,4 @@ By making a contribution to this project, I certify that:
 Contact
 ==
 
-Please, contact the project administrators, [Santiago Badia](mailto:santiago.badia@monash.edu) and [Francesc Verdugo](mailto:fverdugo@cimne.upc.edu), for further questions about licenses and terms of use.
+Please, contact the project administrators, [Francesc Verdugo](mailto:f.verdugo.rojano@vu.nl), [Eric Neiva](mailto:eric.neiva@college-de-france.fr) and [Santiago Badia](mailto:santiago.badia@monash.edu), for further questions about licenses and terms of use.
diff --git a/NEWS.md b/NEWS.md
index eeb9b68..82c962d 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -4,7 +4,12 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
-## [0.8.3] - 2024-08-02
+## [0.9.0] - 2024-01-22
+
+### Added
+- Interface to Algoim v0.2, which provides algoim's high-order quadrature algorithms for domains implicitly-defined by multivariate polynomials and high-order accurate algorithms for computing closest points on implicitly-defined surfaces. Since PR [#76](https://github.com/gridap/GridapEmbedded.jl/pull/76).
+
+## [0.8.3] - 2024-01-02
 
 ### Added
 - Support for FillArrays v1. Since PR [#75](https://github.com/gridap/GridapEmbedded.jl/pull/75).
diff --git a/Project.toml b/Project.toml
index 3e8715c..2a82563 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,11 +1,13 @@
 name = "GridapEmbedded"
 uuid = "8838a6a3-0006-4405-b874-385995508d5d"
-authors = ["Francesc Verdugo <fverdugo@cimne.upc.edu>"]
-version = "0.8.3"
+authors = ["Francesc Verdugo <f.verdugo.rojano@vu.nl>", "Eric Neiva <eric.neiva@college-de-france.fr>", "Santiago Badia <santiago.badia@monash.edu>"]
+version = "0.9.0"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
+Algoim = "0eb9048c-21de-4c7a-bfac-056de1940b74"
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
+CxxWrap = "1f15a43c-97ca-5a2a-ae31-89f07a497df4"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 Gridap = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
 LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d"
@@ -13,6 +15,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MiniQhull = "978d7f02-9e05-4691-894f-ae31a51d76ca"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+algoimWrapper_jll = "3c43aa7b-5398-51f3-8d75-8f051e6faa4d"
 
 [compat]
 AbstractTrees = "0.3.3, 0.4"
diff --git a/docs/make.jl b/docs/make.jl
index c09398c..3bad67a 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -8,7 +8,7 @@ makedocs(;
     ],
     repo="https://github.com/gridap/GridapEmbedded.jl/blob/{commit}{path}#L{line}",
     sitename="GridapEmbedded.jl",
-    authors="Francesc Verdugo <fverdugo@cimne.upc.edu>",
+    authors="Francesc Verdugo <f.verdugo.rojano@vu.nl>, Eric Neiva <eric.neiva@college-de-france.fr> and Santiago Badia <santiago.badia@monash.edu>",
 )
 
 deploydocs(;
diff --git a/src/AgFEM/AgFEM.jl b/src/AgFEM/AgFEM.jl
index 2bf2684..72fed37 100644
--- a/src/AgFEM/AgFEM.jl
+++ b/src/AgFEM/AgFEM.jl
@@ -18,6 +18,7 @@ using GridapEmbedded.Interfaces: compute_subcell_to_inout
 
 export aggregate
 export color_aggregates
+export aggregate_narrow_band
 export AggregateAllCutCells
 export compute_cell_bboxes
 export compute_cell_to_dface_bboxes
diff --git a/src/AgFEM/AggregateBoundingBoxes.jl b/src/AgFEM/AggregateBoundingBoxes.jl
index 4604144..80f26b0 100644
--- a/src/AgFEM/AggregateBoundingBoxes.jl
+++ b/src/AgFEM/AggregateBoundingBoxes.jl
@@ -19,6 +19,18 @@ function init_bboxes(cell_to_coords,cut::EmbeddedDiscretization;in_or_out=IN)
   bgcell_to_cbboxes
 end
 
+function init_bboxes(cell_to_coords,cut_measure::Measure)
+  bgcell_to_cbboxes = init_bboxes(cell_to_coords)
+  quad = get_cell_quadrature(cut_measure)
+  trian = get_triangulation(quad)
+  model = get_active_model(trian)
+  ccell_to_bgcell = get_cell_to_parent_cell(model)
+  for (cc,bc) in enumerate(ccell_to_bgcell)
+    bgcell_to_cbboxes[bc] = compute_subcell_bbox(quad.cell_point[cc])
+  end
+  bgcell_to_cbboxes
+end
+
 function init_cut_bboxes(cut,ccell_to_bgcell,in_or_out)
   subcell_to_inout   = compute_subcell_to_inout(cut,cut.geo)
   subcell_is_in      = lazy_map(i->i==in_or_out,subcell_to_inout)
diff --git a/src/AgFEM/CellAggregation.jl b/src/AgFEM/CellAggregation.jl
index c41df61..90df2fb 100644
--- a/src/AgFEM/CellAggregation.jl
+++ b/src/AgFEM/CellAggregation.jl
@@ -303,3 +303,60 @@ function _color_aggregates_barrier(cell_to_cellin,cell_to_faces,face_to_cells)
 
   cell_to_color
 end
+
+# Specialised methods for Algoim quadratures
+
+function aggregate(bgtrian,cell_to_is_active,cell_to_is_cut,in_or_out)
+  n_cells = length(cell_to_is_active)
+  @assert n_cells == length(cell_to_is_cut)
+
+  cell_to_unit_cut_meas = lazy_map(cell_to_is_active,cell_to_is_cut) do isa, isc
+    !isa ? 0.0 : (isc ? 0.0 : 1.0)
+  end
+
+  cell_to_inoutcut = lazy_map(cell_to_is_active,cell_to_is_cut) do isa, isc
+    !isa ? OUT : (isc ? CUT : IN)
+  end
+
+  cell_to_coords = get_cell_coordinates(bgtrian)
+  model = get_background_model(bgtrian)
+  topo = get_grid_topology(model)
+  D = num_cell_dims(model)
+  cell_to_faces = get_faces(topo,D,D-1)
+  face_to_cells = get_faces(topo,D-1,D)
+  # A hack follows to avoid constructing the actual facet_to_inoutcut array
+  facet_to_inoutcut = fill(in_or_out,num_faces(model,D-1)) 
+
+  threshold = 1.0
+  _aggregate_by_threshold_barrier(
+    threshold,cell_to_unit_cut_meas,facet_to_inoutcut,cell_to_inoutcut,
+    in_or_out,cell_to_coords,cell_to_faces,face_to_cells)
+end
+
+function aggregate_narrow_band(bgtrian,cell_to_is_in_narrow,cell_to_is_active,cell_to_is_cut,in_or_out)
+  n_cells = length(cell_to_is_in_narrow)
+  @assert n_cells == length(cell_to_is_active)
+  @assert n_cells == length(cell_to_is_cut)
+
+  cell_to_unit_cut_meas = lazy_map(cell_to_is_active,cell_to_is_cut) do isa, isc
+    ( isa & !isc ) ? 1.0 : 0.0
+  end
+
+  cell_to_inoutcut = lazy_map(cell_to_is_in_narrow,cell_to_is_active,cell_to_is_cut) do isn, isa, isc
+    !isn ? OUT : ( ( isa & !isc ) ? IN : CUT )
+  end
+
+  cell_to_coords = get_cell_coordinates(bgtrian)
+  model = get_background_model(bgtrian)
+  topo = get_grid_topology(model)
+  D = num_cell_dims(model)
+  cell_to_faces = get_faces(topo,D,D-1)
+  face_to_cells = get_faces(topo,D-1,D)
+  # A hack follows to avoid constructing the actual facet_to_inoutcut array
+  facet_to_inoutcut = fill(in_or_out,num_faces(model,D-1)) 
+
+  threshold = 1.0
+  _aggregate_by_threshold_barrier(
+    threshold,cell_to_unit_cut_meas,facet_to_inoutcut,cell_to_inoutcut,
+    in_or_out,cell_to_coords,cell_to_faces,face_to_cells)
+end
\ No newline at end of file
diff --git a/src/AlgoimUtils/AlgoimUtils.jl b/src/AlgoimUtils/AlgoimUtils.jl
new file mode 100644
index 0000000..c4d36d1
--- /dev/null
+++ b/src/AlgoimUtils/AlgoimUtils.jl
@@ -0,0 +1,410 @@
+module AlgoimUtils
+
+using CxxWrap
+using Algoim
+using Gridap
+
+import Base.prod
+@inline prod(x::Point{N,T}) where {N,T} = prod(x.data)
+
+import Algoim: to_array
+@inline to_array(x::Point{N,T}) where {N,T} = collect(x.data)
+
+using Gridap.Helpers
+using Gridap.ReferenceFEs
+using Gridap.Arrays: collect1d, CompressedArray, Table
+using Gridap.Adaptivity
+using Gridap.Fields: testitem
+using Gridap.CellData
+using Gridap.CellData: GenericCellField, get_data
+using Gridap.CellData: _point_to_cell_cache, _point_to_cell!
+using Gridap.Geometry
+
+using GridapEmbedded.Interfaces
+using GridapEmbedded.Interfaces: Simplex
+
+using MiniQhull
+using FillArrays
+
+import Algoim: CachedLevelSetValue
+import Algoim: CachedLevelSetGradient
+import Algoim: AlgoimCallLevelSetFunction
+import Algoim: normal
+import Gridap.ReferenceFEs: Quadrature
+
+export TriangulationAndMeasure
+export algoim
+export Quadrature
+export is_cell_active
+export restrict_measure
+export fill_cpp_data
+export fill_cpp_data_raw
+export compute_closest_point_projections
+export compute_normal_displacement
+export compute_distance_fe_function
+export delaunaytrian
+export convexhull
+
+struct Algoim <: QuadratureName end
+const algoim = Algoim()
+
+(φ::CachedLevelSetValue{<:CellField})(p,i::Float32) = begin
+  _p = Point(p)
+  (carr,cgetindex,ceval) = φ.c
+  evaluate!(ceval,getindex!(cgetindex,carr,Int(i)),_p)
+end
+
+(φ::CachedLevelSetGradient{<:CellField})(p,i::Float32) = begin
+  _p = Point(p)
+  (carr,cgetindex,ceval) = φ.c
+  _val = evaluate!(ceval,getindex!(cgetindex,carr,Int(i)),_p)
+  ConstCxxRef(to_uvector(to_array(_val)))
+end
+
+function AlgoimCallLevelSetFunction(φ::CellField,∇φ::CellField)
+  φtrian = get_triangulation(φ)
+  ∇φtrian = get_triangulation(∇φ)
+  xφ = testitem(testitem(get_cell_coordinates(φtrian)))
+  x∇φ = testitem(testitem(get_cell_coordinates(∇φtrian)))
+  cellφ = get_array(φ)
+  cell∇φ = get_array(∇φ)
+  cache_φ = (cellφ,array_cache(cellφ),return_cache(testitem(cellφ),xφ))
+  cache_∇φ = (cell∇φ,array_cache(cell∇φ),return_cache(testitem(cell∇φ),x∇φ))
+  AlgoimCallLevelSetFunction{typeof(φ),typeof(∇φ),typeof(cache_φ),typeof(cache_∇φ)}(φ,∇φ,cache_φ,cache_∇φ)
+end
+
+function normal(phi::AlgoimCallLevelSetFunction,x::AbstractVector{<:Point},cell_id::Int=1)
+  map(xi->normal(phi,xi,cell_id),x)
+end
+
+function normal(phi::AlgoimCallLevelSetFunction,trian::Triangulation)
+  f = [ x -> normal(phi,x,i) for i in 1:num_cells(trian) ]
+  GenericCellField(f,trian,PhysicalDomain())
+end
+
+function normal(ls::AlgoimCallLevelSetFunction{<:CellField,<:CellField},x::Point,cell_id::Int=1)
+  (carr,cgetindex,ceval) = ls.cache_∇φ
+  gx = evaluate!(ceval,getindex!(cgetindex,carr,cell_id),x)
+  gx/norm(gx)
+end
+
+function Quadrature(trian::Grid,::Algoim,phi::LevelSetFunction,degree::Int;kwargs...)
+  ctype_polytope = map(get_polytope,get_reffes(trian))
+  @notimplementedif !all(map(is_n_cube,ctype_polytope))
+  cell_to_coords = get_cell_coordinates(trian)
+  cell_to_bboxes = collect1d(lazy_map(a->(a[1],a[end]),cell_to_coords))
+  jls = JuliaFunctionLevelSet(phi,Val{num_dims(trian)}())
+  cell_to_quad = map(enumerate(cell_to_bboxes)) do (cell_id,bbox)
+    bbmin, bbmax = bbox
+    Quadrature(cell_id,bbmin,bbmax,jls,phi,degree;kwargs...)
+  end
+  CompressedArray(cell_to_quad,1:length(cell_to_quad))
+end
+
+function Quadrature(trian::Grid,::Algoim,
+                    phi1::LevelSetFunction,phi2::LevelSetFunction,
+                    degree::Int;kwargs...)
+  ctype_polytope = map(get_polytope,get_reffes(trian))
+  @notimplementedif !all(map(is_n_cube,ctype_polytope))
+  cell_to_coords = get_cell_coordinates(trian)
+  cell_to_bboxes = collect1d(lazy_map(a->(a[1],a[end]),cell_to_coords))
+  jls1 = JuliaFunctionLevelSet(phi1,Val{num_dims(trian)}())
+  jls2 = JuliaFunctionLevelSet(phi2,Val{num_dims(trian)}())
+  cell_to_quad = map(enumerate(cell_to_bboxes)) do (cell_id,bbox)
+    bbmin, bbmax = bbox
+    Quadrature(cell_id,bbmin,bbmax,jls1,jls2,phi1,phi2,degree;kwargs...)
+  end
+  CompressedArray(cell_to_quad,1:length(cell_to_quad))
+end
+
+function Quadrature(cell_id::Int,
+                    xmin::Point{N,T},
+                    xmax::Point{N,T},
+                    jls::LevelSetFunction,
+                    phi::LevelSetFunction,
+                    degree::Int;
+                    phase::Int=CUT) where {N,T}
+  coords, weights = fill_quad_data(jls,phi,xmin,xmax,phase,degree,cell_id)
+  GenericQuadrature(coords,weights,"Algoim quadrature of degree $degree")
+end
+
+function Quadrature(cell_id::Int,
+                    xmin::Point{N,T},
+                    xmax::Point{N,T},
+                    jls1::LevelSetFunction,
+                    jls2::LevelSetFunction,
+                    phi1::LevelSetFunction,
+                    phi2::LevelSetFunction,
+                    degree::Int;
+                    phase1::Int=IN,
+                    phase2::Int=IN) where {N,T}
+  coords, weights = fill_quad_data(jls1,jls2,phi1,phi2,xmin,xmax,
+                                   phase1,phase2,degree,cell_id)
+  GenericQuadrature(coords,weights,"Algoim quadrature of degree $degree")
+end
+
+function is_cell_active(meas::Measure)
+  has_non_empty_quad(x) = num_points(x) > 0
+  lazy_map(has_non_empty_quad,get_data(meas.quad))
+end
+
+function restrict_measure(meas::Measure,trian::Triangulation)
+  ocell_quad = meas.quad
+  cell_quad   = lazy_map(Reindex(ocell_quad.cell_quad),trian.tface_to_mface)
+  cell_point  = lazy_map(Reindex(ocell_quad.cell_point),trian.tface_to_mface)
+  cell_weight = lazy_map(Reindex(ocell_quad.cell_weight),trian.tface_to_mface)
+  dds = ocell_quad.data_domain_style
+  ids = ocell_quad.integration_domain_style
+  Measure(CellQuadrature(cell_quad,cell_point,cell_weight,trian,dds,ids))
+end
+
+function TriangulationAndMeasure(Ωbg::Triangulation,quad::Tuple)
+  msg = "TriangulationAndMeasure can only receive the background triangulation"
+  @notimplementedif num_cells(get_background_model(Ωbg)) != num_cells(Ωbg) msg
+  dΩbg = Measure(Ωbg,quad,data_domain_style=PhysicalDomain())
+  # RMK: This is a hack, but algoim interface does not let you 
+  # know if a (given) cell intersects the interior of the level 
+  # set. The hack consists in inferring this from the size of 
+  # each quadrature. I do not expect this hack implies a lot of 
+  # extra operations with regards to the proper way to do it.
+  cell_to_is_active = is_cell_active(dΩbg)
+  Ωᵃ = Triangulation(Ωbg,cell_to_is_active)
+  dΩᵃ = restrict_measure(dΩbg,Ωᵃ)
+  Ωᵃ,dΩᵃ,cell_to_is_active
+end
+
+using Gridap.Geometry: get_cell_to_parent_cell
+using Gridap.CellData: get_cell_quadrature
+
+function compute_closest_point_projections(Ω::Triangulation,φ;
+    cppdegree::Int=2,trim::Bool=false,limitstol::Float64=1.0e-8)
+  compute_closest_point_projections(get_active_model(Ω),φ,
+    cppdegree=cppdegree,trim=trim,limitstol=limitstol)
+end
+
+function compute_closest_point_projections(model::CartesianDiscreteModel,
+                                           φ::AlgoimCallLevelSetFunction;
+                                           cppdegree::Int=2,
+                                           trim::Bool=false,
+                                           limitstol::Float64=1.0e-8)
+  cdesc = get_cartesian_descriptor(model)
+  partition = Int32[cdesc.partition...]
+  xmin = cdesc.origin
+  xmax = xmin + Point(cdesc.sizes .* partition)
+  fill_cpp_data(φ,partition,xmin,xmax,cppdegree,trim,limitstol)
+end
+
+function compute_closest_point_projections(fespace::FESpace,
+                                           φ::AlgoimCallLevelSetFunction,
+                                           order::Int;
+                                           cppdegree::Int=2,
+                                           trim::Bool=false,
+                                           limitstol::Float64=1.0e-8)
+  trian = get_triangulation(fespace)
+  model = get_active_model(trian)
+  # TO-ANSWER: Do I need the rmodel in this scope or not?
+  rmodel = refine(model,order)
+  cps = compute_closest_point_projections(
+    rmodel,φ,order,cppdegree=cppdegree,trim=trim,limitstol=limitstol)
+  msg = "Is the FE space order the same as the input order?"
+  @assert length(cps) == num_free_dofs(fespace) msg
+  cps = node_to_dof_order(cps,fespace,rmodel,order)
+end
+
+function compute_closest_point_projections(model::AdaptedDiscreteModel,
+                                           φ::AlgoimCallLevelSetFunction,
+                                           order::Int;
+                                           cppdegree::Int=2,
+                                           trim::Bool=false,
+                                           limitstol::Float64=1.0e-8)
+  reffe = ReferenceFE(lagrangian,Float64,order)
+  rfespace = TestFESpace(model,reffe)
+  _rφ = interpolate_everywhere(φ.φ,rfespace)
+  rφ = AlgoimCallLevelSetFunction(_rφ,∇(_rφ))
+  cdesc = get_cartesian_descriptor(get_model(model))
+  partition = Int32[cdesc.partition...]
+  xmin = cdesc.origin
+  xmax = xmin + Point(cdesc.sizes .* partition)
+  fill_cpp_data(rφ,partition,xmin,xmax,cppdegree,trim,limitstol)
+end
+
+function node_to_dof_order(cps,
+                           fespace::FESpace,
+                           rmodel::AdaptedDiscreteModel,
+                           order::Int)
+  D = num_dims(rmodel)
+  cdesc = get_cartesian_descriptor(get_model(rmodel))
+  partition = cdesc.partition
+  orders = tfill(order,Val{D}())
+  ones = tfill(1,Val{D}())
+  range = CartesianIndices(orders.+1) .- CartesianIndex(ones) # 0-based
+  ldof_to_lnode = get_ldof_to_lnode(orders,D)
+  o2n_faces_map = rmodel.glue.o2n_faces_map
+  node_partition = partition .+ 1
+  ncells = num_cells(rmodel.parent)
+  cell_node_ids = lazy_map(1:ncells) do cellid
+    anchor_node = o2n_faces_map[cellid][1]
+    node_ijk = CartesianIndices(partition)[anchor_node]
+    range_cis = node_ijk .+ range
+    node_ids = LinearIndices(node_partition)[range_cis]
+    node_ids[ldof_to_lnode]
+  end
+  cell_dofs_ids = fespace.cell_dofs_ids
+  c1 = array_cache(fespace.cell_dofs_ids)
+  c2 = array_cache(cell_node_ids)
+  ncps = similar(cps)
+  for cellid in 1:ncells
+    dof_ids = getindex!(c1,cell_dofs_ids,cellid)
+    node_ids = getindex!(c2,cell_node_ids,cellid)
+    ncps[dof_ids] = cps[node_ids]
+  end
+  ncps
+end
+
+function get_ldof_to_lnode(orders,D)
+
+  # Generate indices of n-faces and order s.t.
+  # (1) dimension-increasing (2) lexicographic
+  bin_rang_nfaces = tfill(0:1,Val{D}())
+  bin_ids_nfaces = vec(collect(Iterators.product(bin_rang_nfaces...)))
+  sum_bin_ids_nfaces = sum.(bin_ids_nfaces)
+  bin_ids_nfaces = permute!(bin_ids_nfaces,sortperm(sum_bin_ids_nfaces))
+
+  # Generate LIs of basis funs s.t. order by n-faces
+  lids_b = LinearIndices(Tuple([orders[i]+1 for i=1:D]))
+
+  eet = eltype(eltype(bin_ids_nfaces))
+  f(x) = Tuple( x[i] == one(eet) ? (0:0) : (1:orders[i]:(orders[i]+1)) for i in 1:length(x) )
+  g(x) = Tuple( x[i] == one(eet) ? (2:orders[i]) : (0:0) for i in 1:length(x) )
+  rang_nfaces = map(f,bin_ids_nfaces)
+  rang_own_dofs = map(g,bin_ids_nfaces)
+
+  perm = Int64[]
+  for i = 1:length(bin_ids_nfaces)
+    cis_nfaces = CartesianIndices(rang_nfaces[i])
+    cis_own_dofs = CartesianIndices(rang_own_dofs[i])
+    for ci in cis_nfaces
+      ci = ci .+ cis_own_dofs
+      perm = vcat(perm,reshape(lids_b[ci],length(ci)))
+    end
+  end
+
+  perm
+end
+
+function compute_normal_displacement(
+    cps::AbstractVector{<:Point},
+    phi::AlgoimCallLevelSetFunction,
+    fun,
+    dt::Float64,
+    Ω::Triangulation)
+  # Note that cps must be (scalar) DoF-numbered, not lexicographic-numbered
+  searchmethod = KDTreeSearch()
+  cache1 = _point_to_cell_cache(searchmethod,Ω)
+  x_to_cell(x) = _point_to_cell!(cache1, x)
+  point_to_cell = lazy_map(x_to_cell, cps)
+  cell_to_points, _ = make_inverse_table(point_to_cell, num_cells(Ω))
+  cell_to_xs = lazy_map(Broadcasting(Reindex(cps)), cell_to_points)
+  cell_point_xs = CellPoint(cell_to_xs, Ω, PhysicalDomain())
+  fun_xs = evaluate(fun,cell_point_xs)
+  nΓ_xs = evaluate(normal(phi,Ω),cell_point_xs)
+  cell_point_disp = lazy_map(Broadcasting(⋅),fun_xs,nΓ_xs)
+  cache_vals = array_cache(cell_point_disp)
+  cache_ctop = array_cache(cell_to_points)
+  disps = zeros(Float64,length(cps))
+  for cell in 1:length(cell_to_points)
+    pts = getindex!(cache_ctop,cell_to_points,cell)
+    vals = getindex!(cache_vals,cell_point_disp,cell)
+    for (i,pt) in enumerate(pts)
+      val = vals[i]
+      disps[pt] = dt * val
+    end
+  end
+  disps
+end
+
+signed_distance(φ::Function,x,y) = sign(φ(y))*norm(x-y)
+
+signed_distance(φ::T,x,y) where {T<:Number} = sign(φ)*norm(x-y)
+
+function _compute_signed_distance(
+    φ::AlgoimCallLevelSetFunction{<:Function,<:Function},
+    cps::Vector{<:Point},cos::Matrix{<:Point})
+  _dist(x,y) = signed_distance(φ.φ,x,y)
+  lazy_map(_dist,cps,cos)
+end
+
+function _compute_signed_distance(
+    φ::AlgoimCallLevelSetFunction{<:CellField,<:CellField},
+    cps::Vector{<:Point},cos::Matrix{<:Point})
+  φs = get_free_dof_values(φ.φ)
+  lazy_map(signed_distance,φs,cps,cos)
+end
+
+function compute_distance_fe_function(
+    bgmodel::CartesianDiscreteModel,
+    fespace::FESpace,
+    φ::AlgoimCallLevelSetFunction,
+    order::Int;
+    cppdegree::Int=2)
+  cps = compute_closest_point_projections(
+    fespace,φ,order,cppdegree=cppdegree)
+  rmodel = refine(bgmodel,order)
+  cos = get_node_coordinates(rmodel)
+  cos = node_to_dof_order(cos,fespace,rmodel,order)
+  dists = _compute_signed_distance(φ,cps,cos) 
+  FEFunction(fespace,dists)
+end
+
+abstract type QhullType end
+
+struct DelaunayTrian <: QhullType end
+const delaunaytrian = DelaunayTrian()
+
+struct ConvexHull <: QhullType end
+const convexhull = ConvexHull()
+
+get_flags(::QhullType) = @abstractmethod
+get_flags(::DelaunayTrian) = "qhull d Qt Qbb Qc Qz"
+get_flags(::ConvexHull) = "qhull Qt Qc"
+
+get_dimension(::QhullType,dim) = @abstractmethod
+get_dimension(::DelaunayTrian,dim) = dim
+get_dimension(::ConvexHull,dim) = dim-1
+
+import Gridap.Visualization: visualization_data
+
+function visualization_data(meas::Measure,filename;cellfields=Dict(),qhulltype=DelaunayTrian())
+  node_coordinates = collect(Iterators.flatten(meas.quad.cell_point.values))
+  grid = _to_grid(node_coordinates,qhulltype)
+  ndata = Dict()
+  for (k,v) in cellfields
+    pts = get_cell_points(meas)
+    eval = evaluate(v,pts)
+    ndata[k] = collect(Iterators.flatten(eval))
+  end
+  visualization_data(grid,filename,nodaldata=ndata)
+end
+
+function visualization_data(meas::Vector{<:Measure},filename;cellfields=Dict(),qhulltype=DelaunayTrian())
+  node_coordinates = vcat(map(m->collect(Iterators.flatten(m.quad.cell_point.values)),meas)...)
+  grid = _to_grid(node_coordinates,qhulltype)
+  ndata = Dict()
+  for (k,v) in cellfields
+    pts = map(m->get_cell_points(m),meas)
+    eval = map(p->evaluate(v,p),pts)
+    ndata[k] = vcat(map(e->collect(Iterators.flatten(e)),eval)...)
+  end
+  visualization_data(grid,filename,nodaldata=ndata)
+end
+
+function _to_grid(node_coordinates::Vector{<:Point{Dp,Tp}},qhulltype) where {Dp,Tp}
+  d = get_dimension(qhulltype,Dp)
+  connectivity = delaunay(reinterpret(node_coordinates),get_flags(qhulltype))[1:(d+1),:]
+  cell_node_ids = Table(collect(eachcol(connectivity)))
+  reffes = [LagrangianRefFE(Float64,Simplex(Val{d}()),1)]
+  cell_types = collect(Fill(Int8(1),length(cell_node_ids)))
+  UnstructuredGrid(node_coordinates,cell_node_ids,reffes,cell_types)
+end
+
+end # module
\ No newline at end of file
diff --git a/src/Exports.jl b/src/Exports.jl
index 8a6f2a5..66974a1 100644
--- a/src/Exports.jl
+++ b/src/Exports.jl
@@ -44,9 +44,23 @@ end
 @publish AgFEM AgFEMSpace
 @publish AgFEM aggregate
 @publish AgFEM color_aggregates
+@publish AgFEM aggregate_narrow_band
 @publish AgFEM AggregateCutCellsByThreshold
 @publish AgFEM AggregateAllCutCells
 @publish AgFEM compute_cell_bboxes
 @publish AgFEM compute_cell_to_dface_bboxes
 
+@publish AlgoimUtils algoim
+@publish AlgoimUtils fill_quad_data
+@publish AlgoimUtils AlgoimCallLevelSetFunction
+@publish AlgoimUtils TriangulationAndMeasure
+@publish AlgoimUtils normal
+@publish AlgoimUtils fill_cpp_data
+@publish AlgoimUtils fill_cpp_data_raw
+@publish AlgoimUtils compute_closest_point_projections
+@publish AlgoimUtils compute_normal_displacement
+@publish AlgoimUtils compute_distance_fe_function
+@publish AlgoimUtils delaunaytrian
+@publish AlgoimUtils convexhull
+
 @publish MomentFittedQuadratures momentfitted
\ No newline at end of file
diff --git a/src/GridapEmbedded.jl b/src/GridapEmbedded.jl
index b84e57e..175cee0 100644
--- a/src/GridapEmbedded.jl
+++ b/src/GridapEmbedded.jl
@@ -10,6 +10,8 @@ include("AgFEM/AgFEM.jl")
 
 include("MomentFittedQuadratures/MomentFittedQuadratures.jl")
 
+include("AlgoimUtils/AlgoimUtils.jl")
+
 include("Exports.jl")
 
 end # module
diff --git a/test/AlgoimUtilsTests/AlgoimInterfaceTests.jl b/test/AlgoimUtilsTests/AlgoimInterfaceTests.jl
new file mode 100644
index 0000000..d1f64ec
--- /dev/null
+++ b/test/AlgoimUtilsTests/AlgoimInterfaceTests.jl
@@ -0,0 +1,108 @@
+module AlgoimInterfaceTests
+
+using Test
+using Gridap
+using Gridap.ReferenceFEs
+using GridapEmbedded
+using GridapEmbedded.Interfaces
+
+function run_case_raw_point(degree,phase)
+
+  u(x) = (x[1]*x[1]/(1.0*1.0)+x[2]*x[2]/(0.5*0.5)) - 1.0
+  gradu(x) = VectorValue(2.0*x[1]/(1.0*1.0),2.0*x[2]/(0.5*0.5))
+
+  phi = AlgoimCallLevelSetFunction(u,gradu)
+
+  xmin = Point(-1.1,-1.1)
+  xmax = Point(1.1,1.1)
+
+  _, weights = fill_quad_data(phi,xmin,xmax,phase,degree)
+
+  s = ∑(weights)
+  if phase == IN
+    @test s ≈ 1.6055602335693198
+  elseif phase == CUT
+    @test s ≈ 5.04463087731472
+  else
+    error()
+  end
+
+end
+
+function run_case_function(cells,degree,phase)
+
+  domain = (-1.1,1.1,-1.1,1.1)
+
+  model = CartesianDiscreteModel(domain,cells)
+  Ω = Triangulation(model)
+
+  u(x) = (x[1]*x[1]/(1.0*1.0)+x[2]*x[2]/(0.5*0.5)) - 1.0
+  gradu(x) = VectorValue(2.0*x[1]/(1.0*1.0),2.0*x[2]/(0.5*0.5))
+
+  phi = AlgoimCallLevelSetFunction(u,gradu)
+  quad = Quadrature(algoim,phi,degree,phase=phase)
+  dΩ = Measure(Ω,quad,data_domain_style=PhysicalDomain())
+
+  s = ∫(1)dΩ
+  if phase == IN
+    @test ∑(s) ≈ 1.570796326817732
+  elseif phase == CUT
+    @test ∑(s) ≈ 4.844224109684756
+  else
+    error()
+  end
+
+end
+
+function run_case_fe_function(cells,order,degree,phase)
+
+  domain = (-1.1,1.1,-1.1,1.1,-1.1,1.1)
+
+  model = CartesianDiscreteModel(domain,cells)
+  Ω = Triangulation(model)
+
+  reffe = ReferenceFE(lagrangian,Float64,order)
+  V = TestFESpace(Ω,reffe)
+  U = TrialFESpace(V)
+
+  u(x) = (x[1]*x[1]/(1.0*1.0)+x[2]*x[2]/(0.5*0.5)+x[3]*x[3]/(0.33*0.33)) - 1.0
+  gradu(x) = VectorValue(2.0*x[1]/(1.0*1.0),2.0*x[2]/(0.5*0.5),2.0*x[3]/(0.33*0.33))
+
+  uₕ = interpolate(u,U)
+  ∇uₕ = ∇(uₕ)
+
+  phi = AlgoimCallLevelSetFunction(uₕ,∇uₕ)
+  quad = Quadrature(algoim,phi,degree,phase=phase)
+  dΩ = Measure(Ω,quad,data_domain_style=PhysicalDomain())
+
+  s = ∫(1)dΩ
+  if phase == IN
+    @test ∑(s) ≈ 0.6911505144863921
+  elseif phase == CUT
+    @test ∑(s) ≈ 4.382658365204242
+  else
+    error()
+  end
+
+end
+
+f(x) = (x[1]*x[1]/(1.0*1.0)+x[2]*x[2]/(0.5*0.5)+x[3]*x[3]/(0.33*0.33)) - 1.0
+gradf(x) = VectorValue(2.0*x[1]/(1.0*1.0),2.0*x[2]/(0.5*0.5),2.0*x[3]/(0.33*0.33))
+
+phi = AlgoimCallLevelSetFunction(f,gradf)
+@test phi(VectorValue(1.0,0.5,0.33)) == 2.0
+@test all(phi.∇φ(VectorValue(1.0,0.5,0.33)) .≈ VectorValue(2.0, 4.0, 6.0606060606060606))
+
+run_case_raw_point(3,IN)
+run_case_raw_point(3,CUT)
+
+run_case_function((64,64),3,IN)
+run_case_function((64,66),3,IN)
+
+run_case_function((64,64),3,CUT)
+run_case_function((64,66),3,CUT)
+
+run_case_fe_function((32,32,32),2,3,IN)
+run_case_fe_function((32,32,32),2,3,CUT)
+
+end # module
\ No newline at end of file
diff --git a/test/AlgoimUtilsTests/ClosestPointTests.jl b/test/AlgoimUtilsTests/ClosestPointTests.jl
new file mode 100644
index 0000000..0779288
--- /dev/null
+++ b/test/AlgoimUtilsTests/ClosestPointTests.jl
@@ -0,0 +1,44 @@
+module ClosestPointTests
+
+using Test
+using CxxWrap
+using Gridap
+using GridapEmbedded
+
+const IN = -1
+const OUT = 1
+const CUT = 0
+
+function run_case(degree)
+
+  xmin = Point(-1.1,-1.1,-1.1)
+  xmax = Point(1.1,1.1,1.1)
+  partition = Int32[16,16,16]
+  
+  domain = (-1.1,1.1,-1.1,1.1,-1.1,1.1)
+  bgmodel = CartesianDiscreteModel(domain,partition)
+  Ω = Triangulation(bgmodel)
+  
+  reffeᵠ = ReferenceFE(lagrangian,Float64,1)
+  V = TestFESpace(Ω,reffeᵠ,conformity=:L2)
+  
+  f(x) = (x[1]*x[1]/(0.5*0.5)+x[2]*x[2]/(0.5*0.5)+x[3]*x[3]/(0.5*0.5)) - 1.0
+  function gradf(x::V) where {V}
+    V([2.0*x[1]/(0.5*0.5),2.0*x[2]/(0.5*0.5),2.0*x[3]/(0.5*0.5)])
+  end
+  
+  fₕ = interpolate_everywhere(f,V)
+  phi = AlgoimCallLevelSetFunction(fₕ,∇(fₕ))
+  
+  coords = fill_cpp_data(phi,partition,xmin,xmax,degree)
+  @test maximum(f.(coords)) < 1.0e-4
+
+end
+
+run_case(2)
+run_case(3)
+run_case(4)
+run_case(5)
+run_case(-1)
+
+end # module
\ No newline at end of file
diff --git a/test/AlgoimUtilsTests/DualQuadraturesTests.jl b/test/AlgoimUtilsTests/DualQuadraturesTests.jl
new file mode 100644
index 0000000..13f2fad
--- /dev/null
+++ b/test/AlgoimUtilsTests/DualQuadraturesTests.jl
@@ -0,0 +1,99 @@
+module DualQuadraturesTests
+
+using Test
+using Gridap
+using Gridap.ReferenceFEs
+using GridapEmbedded
+using GridapEmbedded.Interfaces
+
+# 3D funs
+
+u¹(x) = ( (x[1]-0.25) * (x[1]-0.25) +
+          (x[2]-0.25) * (x[2]-0.25) + 
+          (x[3]-0.25) * (x[3]-0.25) ) / ( 0.20 * 0.20 ) - 1.0
+
+u²(x) = ( (x[1]-0.75) * (x[1]-0.75) +
+          (x[2]-0.75) * (x[2]-0.75) + 
+          (x[3]-0.75) * (x[3]-0.75) ) / ( 0.22 * 0.22 ) - 1.0
+
+u³(x) = ( (x[1]-0.50) * (x[1]-0.50) +
+          (x[2]-0.50) * (x[2]-0.50) + 
+          (x[3]-0.50) * (x[3]-0.50) ) / ( 0.40 * 0.40 ) - 1.0
+
+# 2D funs
+
+v¹(x) = ( (x[1]-0.25) * (x[1]-0.25) +
+          (x[2]-0.25) * (x[2]-0.25) ) / ( 0.20 * 0.20 ) - 1.0
+
+v²(x) = ( (x[1]-0.75) * (x[1]-0.75) +
+          (x[2]-0.75) * (x[2]-0.75) ) / ( 0.22 * 0.22 ) - 1.0
+
+v³(x) = ( (x[1]-0.50) * (x[1]-0.50) +
+          (x[2]-0.50) * (x[2]-0.50) ) / ( 0.40 * 0.40 ) - 1.0
+
+function run_case_fe_function(domain,cells,order,degree,phase¹,phase²,u¹,u²)
+
+  model = CartesianDiscreteModel(domain,cells)
+  Ω = Triangulation(model)
+
+  reffe = ReferenceFE(lagrangian,Float64,order)
+  V = TestFESpace(Ω,reffe)
+  U = TrialFESpace(V)
+
+  uₕ¹  = interpolate(u¹,U)
+  ∇uₕ¹ = ∇(uₕ¹)
+
+  uₕ²  = interpolate(u²,U)
+  ∇uₕ² = ∇(uₕ²)
+
+  phi¹ = AlgoimCallLevelSetFunction(uₕ¹,∇uₕ¹)
+  phi² = AlgoimCallLevelSetFunction(uₕ²,∇uₕ²)
+
+  quad = Quadrature(algoim,phi¹,phi²,degree,phase1=phase¹,phase2=phase²)
+  dΩ = Measure(Ω,quad,data_domain_style=PhysicalDomain())
+
+  s = ∑(∫(1)dΩ)
+  if ( phase¹ == IN ) && ( phase² == OUT )
+    @test s ≈ 0.03351032164371535
+  elseif ( phase¹ == OUT ) && ( phase² == IN )
+    @test s ≈ 0.04460223810211001
+  elseif ( phase¹ == IN ) && ( phase² == CUT )
+    @test_skip s ≈ 0.11301756191551385 # Different precision in GitHub's CI workflow
+  elseif ( phase¹ == OUT ) && ( phase² == CUT )
+    @test_skip s ≈ 1.8976021283682756 # Different precision in GitHub's CI workflow
+  end
+
+end
+
+domain3D = (0.0,1.1,0.0,1.1,0.0,1.1)
+
+run_case_fe_function(domain3D,(24,24,24),2,5,IN,OUT,u¹,u²)
+# run_case_fe_function(domain3D,(24,24,24),2,5,OUT,IN,u¹,u²)
+
+run_case_fe_function(domain3D,(24,24,24),2,5,IN,CUT,u¹,u³)
+# run_case_fe_function(domain3D,(24,24,24),2,5,OUT,CUT,u¹,u³)
+
+# domain2D = (0.0,1.1,0.0,1.1)
+# domain3D = (0.0,1.1,0.0,1.1,0.0,1.1)
+# order  = 2
+# phaseA = IN
+# phaseB = OUT
+# u = u¹
+# v = u³
+
+# T = pi*16/25 # intersection 3D
+# # T = pi*4/5   # intersection 2D
+# # T = pi*(0.2^2+0.22^2) # volumes 2D
+# # T = (4/3)*pi*(0.2^3+0.22^3) # volumes 3D
+# # T = 4*pi*(0.2^2+0.22^2) # non-intersecting surfaces 3D
+
+# for deg in [3,4,5]
+#   @info "Case: degree = $deg"
+#   for n in [6,12,24,48]
+#     A = run_case_fe_function(domain3D,(n,n,n),2,deg,phaseA,CUT,u¹,u³)
+#     B = run_case_fe_function(domain3D,(n,n,n),2,deg,phaseB,CUT,u¹,u³)
+#     @show A+B-T
+#   end
+# end
+
+end # module
\ No newline at end of file
diff --git a/test/AlgoimUtilsTests/PoissonAlgoimTests.jl b/test/AlgoimUtilsTests/PoissonAlgoimTests.jl
new file mode 100644
index 0000000..e359b91
--- /dev/null
+++ b/test/AlgoimUtilsTests/PoissonAlgoimTests.jl
@@ -0,0 +1,64 @@
+module PoissonAlgoimTests
+
+using Test
+using Gridap
+using Gridap.ReferenceFEs
+using GridapEmbedded
+using GridapEmbedded.Interfaces
+
+function run_poisson(domain,cells,order,solution_degree)
+
+  φ(x) = x[1]-x[2]
+  ∇φ(x) = VectorValue(1.0,-1.0,0.0)
+  phi = AlgoimCallLevelSetFunction(φ,∇φ)
+
+  u(x) = sum(x)^solution_degree
+  f(x) = -Δ(u)(x)
+  ud(x) = u(x)
+
+  model = CartesianDiscreteModel(domain,cells)
+  Ω = Triangulation(model)
+
+  degree = order == 1 ? 3 : 2*order
+  vquad = Quadrature(algoim,phi,degree,phase=IN)
+  Ωᵃ,dΩᵃ = TriangulationAndMeasure(Ω,vquad)
+  squad = Quadrature(algoim,phi,degree)
+  _,dΓ = TriangulationAndMeasure(Ω,squad)
+  n_Γ = normal(phi,Ω)
+
+  reffe = ReferenceFE(lagrangian,Float64,order)
+  V = TestFESpace(Ωᵃ,reffe,conformity=:H1,dirichlet_tags="boundary")
+  U = TrialFESpace(V,ud)
+
+  # Nitsche method
+  γᵈ = 2.0*order^2
+  h = (domain[2]-domain[1])/cells[1]
+
+  a(u,v) =
+    ∫( ∇(v)⋅∇(u) )dΩᵃ +
+    ∫( (γᵈ/h)*v*u - v*(n_Γ⋅∇(u)) - (n_Γ⋅∇(v))*u )dΓ
+
+  l(v) =
+    ∫( v*f )dΩᵃ +
+    ∫( (γᵈ/h)*v*ud - (n_Γ⋅∇(v))*ud )dΓ
+
+  # FE problem
+  op = AffineFEOperator(a,l,U,V)
+  uₕ = solve(op)
+
+  eₕ = u - uₕ
+
+  l2(u) = √(∑( ∫( u*u )dΩᵃ ))
+  h1(u) = √(∑( ∫( u*u + ∇(u)⋅∇(u) )dΩᵃ ))
+
+  el2 = l2(eₕ)
+  eh1 = h1(eₕ)
+
+  @test el2 < 1.0e-014
+  @test eh1 < 1.0e-013
+  
+end
+
+run_poisson((-1.1,1.1,-1.1,1.1,-1.1,1.1),(8,8,8),1,1)
+
+end # module
\ No newline at end of file
diff --git a/test/AlgoimUtilsTests/QuadratureDegreeTests.jl b/test/AlgoimUtilsTests/QuadratureDegreeTests.jl
new file mode 100644
index 0000000..530a792
--- /dev/null
+++ b/test/AlgoimUtilsTests/QuadratureDegreeTests.jl
@@ -0,0 +1,167 @@
+module QuadratureDegreeTests
+
+using Test
+using Gridap
+using Gridap.ReferenceFEs
+using GridapEmbedded
+using GridapEmbedded.Interfaces
+
+φ(x) = x[1]-x[2]
+∇φ(x) = VectorValue(1.0,-1.0,0.0)
+phi = AlgoimCallLevelSetFunction(φ,∇φ)
+
+function run_diffusion_reaction_bulk(domain,cells,order,degree)
+
+  u(x) = sum(x)^order
+  f(x) = u(x)-Δ(u)(x)
+  ud(x) = u(x)
+
+  model = CartesianDiscreteModel(domain,cells)
+  Ω = Triangulation(model)
+
+  vquad = Quadrature(algoim,phi,degree,phase=IN)
+  Ωᵃ,dΩᵃ,is_a = TriangulationAndMeasure(Ω,vquad)
+  squad = Quadrature(algoim,phi,degree)
+  _,dΓ,is_c = TriangulationAndMeasure(Ω,squad)
+  n_Γ = normal(phi,Ω)
+
+  cell_to_cellin = aggregate(Ω,is_a,is_c,IN)
+
+  reffe = ReferenceFE(lagrangian,Float64,order)
+  Vstd = TestFESpace(Ωᵃ,reffe,conformity=:H1,dirichlet_tags="boundary")
+
+  V = AgFEMSpace(Vstd,cell_to_cellin)
+  U = TrialFESpace(V,ud)
+
+  # Nitsche method
+  γᵈ = 5.0*order^2
+  h = (domain[2]-domain[1])/cells[1]
+
+  a(u,v) =
+    ∫( v⋅u )dΩᵃ +
+    ∫( ∇(v)⋅∇(u) )dΩᵃ +
+    ∫( (γᵈ/h)*v*u - v*(n_Γ⋅∇(u)) - (n_Γ⋅∇(v))*u )dΓ
+
+  l(v) =
+    ∫( v*f )dΩᵃ +
+    ∫( (γᵈ/h)*v*ud - (n_Γ⋅∇(v))*ud )dΓ
+
+  # FE problem
+  op = AffineFEOperator(a,l,U,V)
+  uₕ = solve(op)
+
+  eₕ = u - uₕ
+
+  l2(u) = √(∑( ∫( u*u )dΩᵃ ))
+  h1(u) = √(∑( ∫( u*u + ∇(u)⋅∇(u) )dΩᵃ ))
+
+  el2 = l2(eₕ)
+  eh1 = h1(eₕ)
+
+  el2, eh1
+
+end
+
+left_project(u,n) = u - n⊗(n⋅u)
+∇ᵈ(u,n) = left_project(∇(u),n)
+
+function run_diffusion_reaction_surface(domain,cells,order,degree)
+
+  u(x) = sum(x)^order
+  f(x) = u(x)-Δ(u)(x)
+  ud(x) = u(x)
+
+  model = CartesianDiscreteModel(domain,cells)
+  Ω = Triangulation(model)
+
+  squad = Quadrature(algoim,phi,degree)
+  Ωc,dΓ = TriangulationAndMeasure(Ω,squad)
+
+  n_Γ = normal(phi,Ωc)
+  dΩc = Measure(Ωc,2*order)
+
+  reffe = ReferenceFE(lagrangian,Float64,order)
+
+  V = TestFESpace(Ωc,reffe,dirichlet_tags="boundary")
+  U = TrialFESpace(V,ud)
+
+  m(u,v) = ∫( v⋅u )dΓ
+  a(u,v) = ∫( ∇ᵈ(v,n_Γ)⋅∇ᵈ(u,n_Γ) )dΓ
+
+  # Γ-orthogonal derivative volume stabilisation
+  h = (domain[2]-domain[1])/cells[1]
+  γˢ = 10.0*h
+  s(u,v) = ∫( γˢ*((n_Γ⋅∇(v))⊙(n_Γ⋅∇(u))) )dΩc
+
+  aₛ(u,v) = m(u,v) + a(u,v) + s(u,v)
+  l(v) = ∫( v⊙f )dΓ
+
+  op = AffineFEOperator(aₛ,l,U,V)
+  uₕ = solve(op)
+  eₕ = uₕ - u
+
+  l2(v) = √(∑(∫(v⋅v)dΓ))
+  l2(eₕ),l2(∇ᵈ(eₕ,n_Γ))
+
+end
+
+domain = (-1.1,1.1,-1.1,1.1,-1.1,1.1)
+partition = (8,8,8)
+
+# BULK tests
+
+order = 1
+degree = 2
+el2, eh1 = run_diffusion_reaction_bulk(domain,partition,order,degree)
+@test el2 < 1.0e-015
+@test eh1 < 1.0e-014
+
+order = 2
+degree = 3
+el2, eh1 = run_diffusion_reaction_bulk(domain,partition,order,degree)
+@test el2 < 1.0e-013
+@test eh1 < 1.0e-012
+
+# For order > 2, expect accuracy loss with Lagrangian AgFEM
+
+# order = 3
+# degree = 5
+# el2, eh1 = run_diffusion_reaction_bulk(domain,partition,order,degree)
+# @test el2 < 1.0e-011
+# @test eh1 < 1.0e-009
+
+# order = 4
+# degree = 6
+# el2, eh1 = run_diffusion_reaction_bulk(domain,partition,order,degree)
+# @test el2 < 1.0e-007
+# @test eh1 < 1.0e-005
+
+# SURFACE tests
+
+order = 1
+degree = 2
+el2, eh1 = run_diffusion_reaction_surface(domain,partition,order,degree)
+@test el2 < 1.0e-015
+@test eh1 < 1.0e-014
+
+order = 2
+degree = 3
+el2, eh1 = run_diffusion_reaction_surface(domain,partition,order,degree)
+@test el2 < 1.0e-013
+@test eh1 < 1.0e-012
+
+# For order > 2, expect accuracy loss with Lagrangian AgFEM
+
+# order = 3
+# degree = 5
+# el2, eh1 = run_diffusion_reaction_surface(domain,partition,order,degree)
+# @test el2 < 1.0e-009
+# @test eh1 < 1.0e-008
+
+# order = 4
+# degree = 6
+# el2, eh1 = run_diffusion_reaction_surface(domain,partition,order,degree)
+# @test el2 < 1.0e-007
+# @test eh1 < 1.0e-005
+
+end # module
\ No newline at end of file
diff --git a/test/AlgoimUtilsTests/VisualizationTests.jl b/test/AlgoimUtilsTests/VisualizationTests.jl
new file mode 100644
index 0000000..9118178
--- /dev/null
+++ b/test/AlgoimUtilsTests/VisualizationTests.jl
@@ -0,0 +1,43 @@
+module VisualizationTests
+
+  using CxxWrap
+  using Gridap
+  using GridapEmbedded
+  using Gridap.ReferenceFEs
+
+  const IN = -1
+  const OUT = 1
+  const CUT = 0
+
+  order = 1
+  domain = (-1.1,1.1,-1.1,1.1)
+
+  f(x) = (x[1]*x[1]/(0.5*0.5)+x[2]*x[2]/(0.5*0.5)) - 1.0
+  function gradf(x::V) where {V}
+      V([2.0*x[1]/(0.5*0.5),2.0*x[2]/(0.5*0.5)])
+  end
+
+  n = 14
+  degree = 3
+  cppdegree = -1
+
+  partition = Int32[n,n]
+  bgmodel = CartesianDiscreteModel(domain,partition)
+  Ω = Triangulation(bgmodel)
+
+  reffeᵠ = ReferenceFE(lagrangian,Float64,order)
+  V = TestFESpace(Ω,reffeᵠ)
+
+  fₕ = interpolate_everywhere(f,V)
+  phi = AlgoimCallLevelSetFunction(fₕ,∇(fₕ))
+
+  squad = Quadrature(algoim,phi,degree,phase=CUT)
+  _,dΓ = TriangulationAndMeasure(Ω,squad)
+
+  vquad = Quadrature(algoim,phi,degree,phase=IN)
+  _,dΩ = TriangulationAndMeasure(Ω,vquad)
+
+  writevtk(dΓ,"res_sur",cellfields=["f"=>fₕ],qhulltype=convexhull)
+  writevtk([dΩ,dΓ],"res_vol",cellfields=["f"=>fₕ])
+
+end # module
\ No newline at end of file
diff --git a/test/AlgoimUtilsTests/VolumeConservationTests.jl b/test/AlgoimUtilsTests/VolumeConservationTests.jl
new file mode 100644
index 0000000..374d26d
--- /dev/null
+++ b/test/AlgoimUtilsTests/VolumeConservationTests.jl
@@ -0,0 +1,69 @@
+module VolumeConservationTests
+
+using Test
+using CxxWrap
+using Gridap
+using GridapEmbedded
+
+using Gridap.ReferenceFEs
+
+const IN = -1
+const OUT = 1
+const CUT = 0
+
+domain = (-1.1,1.1,-1.1,1.1)
+
+f(x) = (x[1]*x[1]/(0.5*0.5)+x[2]*x[2]/(0.5*0.5)) - 1.0
+function gradf(x::V) where {V}
+    V([2.0*x[1]/(0.5*0.5),2.0*x[2]/(0.5*0.5)])
+end
+
+function run_case(n::Int,order::Int,cppdegree::Int)
+
+  partition = Int32[n,n]
+  bgmodel = CartesianDiscreteModel(domain,partition)
+  Ω = Triangulation(bgmodel)
+
+  reffeᵠ = ReferenceFE(lagrangian,Float64,order)
+  V = TestFESpace(Ω,reffeᵠ)
+
+  fₕ = interpolate_everywhere(f,V)
+  phi = AlgoimCallLevelSetFunction(fₕ,∇(fₕ))
+
+  degree = order == 1 ? 3 : 2*order
+  vquad = Quadrature(algoim,phi,degree,phase=IN)
+  dΩ = Measure(Ω,vquad,data_domain_style=PhysicalDomain())
+  vol_1 = ∑(∫(1)dΩ)
+  # @show vol_1
+
+  cps = compute_closest_point_projections(V,phi,order,cppdegree=cppdegree)
+
+  reffeᶜ = ReferenceFE(lagrangian,VectorValue{2,Float64},order)
+  W = TestFESpace(Ω,reffeᶜ)
+
+  g(x) = VectorValue(1.0,0.0)
+  gₕ = interpolate_everywhere(g,W)
+
+  dt = 0.2
+  _phi₂ = compute_normal_displacement(cps,phi,gₕ,dt,Ω)
+  _phi₂ = get_free_dof_values(fₕ) - _phi₂
+  _phi₂ = FEFunction(V,_phi₂)
+  phi₂ = AlgoimCallLevelSetFunction(_phi₂,∇(_phi₂))
+
+  vquad = Quadrature(algoim,phi₂,degree,phase=IN)
+  dΩ₂ = Measure(Ω,vquad,data_domain_style=PhysicalDomain())
+  # writevtk(dΩ₂,"res_vol")
+  vol_2 = ∑(∫(1)dΩ₂)
+  # @show vol_2
+
+  abs(vol_1-vol_2)/abs(vol_1)
+
+end
+
+order = 2 # Affects CP accuracy _and_ volume conservation
+cppdegree = 2 # Affects CP accuracy, but not volume conservation
+
+@test run_case(12,order,cppdegree) < 1.0e-5
+# @test run_case(24,order,cppdegree) < 1.0e-6
+
+end # module
\ No newline at end of file
diff --git a/test/AlgoimUtilsTests/runtests.jl b/test/AlgoimUtilsTests/runtests.jl
new file mode 100644
index 0000000..f24d881
--- /dev/null
+++ b/test/AlgoimUtilsTests/runtests.jl
@@ -0,0 +1,19 @@
+module AlgoimTests
+
+using Test
+
+@time @testset "AlgoimInterface" begin include("AlgoimInterfaceTests.jl") end
+
+@time @testset "QuadratureDegree" begin include("QuadratureDegreeTests.jl") end
+
+@time @testset "DualQuadratures" begin include("DualQuadraturesTests.jl") end
+
+@time @testset "PoissonAlgoim" begin include("PoissonAlgoimTests.jl") end
+
+@time @testset "ClosestPoint" begin include("ClosestPointTests.jl") end
+
+@time @testset "VolumeConservation" begin include("VolumeConservationTests.jl") end
+
+@time @testset "Visualization" begin include("VisualizationTests.jl") end
+
+end # module
\ No newline at end of file