diff --git a/NEWS.md b/NEWS.md
index 0d5f02c88..e8c1b5787 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -14,6 +14,9 @@
 - `PointLocationHistory` was not marked as public. This has been fixed. See [#198](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/198).
 - Fixed an issue with missing docstrings and duplicate docstrings in the documentation. See [#198](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/198).
 - `copy` and `deepcopy` are now correctly implemented for `PolygonTree`s and `PolygonHierarchy`s. See [#199](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/199)
+- Implemented `copy` and `deepcopy` for `Triangulation` and `VoronoiTessellation`.  See [#201](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/201). 
+- Fixed a bug with `Triangulation`s `polygon_hierarchy` not being correctly aliased with the `polygon_hierarchy` from the `BoundaryEnricher`, and similarly for the `boundary_edge_map`.  See [#201](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/201).
+- Implemented `==` for `VoronoiTessellation`. See [#201](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/201).
 
 ## 1.5.0
 
diff --git a/src/algorithms/triangulation/constrained_triangulation.jl b/src/algorithms/triangulation/constrained_triangulation.jl
index 4cc653888..3f02b0817 100644
--- a/src/algorithms/triangulation/constrained_triangulation.jl
+++ b/src/algorithms/triangulation/constrained_triangulation.jl
@@ -316,7 +316,13 @@ function remake_triangulation_with_constraints(tri::Triangulation, segments, bou
     boundary_curves = get_boundary_curves(tri)
     I = integer_type(tri)
     E = edge_type(tri)
-    boundary_edge_map = construct_boundary_edge_map(boundary_nodes, I, E)
+    if is_curve_bounded(tri)
+        boundary_edge_map = get_boundary_edge_map(get_boundary_enricher(tri))
+        empty!(boundary_edge_map)
+        _construct_boundary_edge_map!(boundary_edge_map, boundary_nodes)
+    else
+        boundary_edge_map = construct_boundary_edge_map(boundary_nodes, I, E)
+    end
     ghost_vertex_map = get_ghost_vertex_map(tri)
     new_ghost_vertex_map = construct_ghost_vertex_map(boundary_nodes, I) # Delay putting these in until we are done with the triangulation
     ghost_vertex_ranges = get_ghost_vertex_ranges(tri)
@@ -349,7 +355,7 @@ function remake_triangulation_with_constraints(tri::Triangulation, segments, bou
 end
 
 """
-    replace_ghost_vertex_information(tri::Triangulation, ghost_vertex_map, ghost_vertex_ranges) -> Triangulation
+    replace_ghost_vertex_information(tri::Triangulation, ghost_vertex_map, ghost_vertex_ranges, polygon_hierarchy) -> Triangulation
 
 Replaces the ghost vertex information in `tri` with `ghost_vertex_map` and `ghost_vertex_ranges`, using the results from 
 [`remake_triangulation_with_constraints`](@ref).
@@ -358,11 +364,12 @@ Replaces the ghost vertex information in `tri` with `ghost_vertex_map` and `ghos
 - `tri::Triangulation`: The triangulation to remake.
 - `ghost_vertex_map`: The ghost vertex map to add to the triangulation.
 - `ghost_vertex_ranges`: The ghost vertex ranges to add to the triangulation.
+- `polygon_hierarchy`: The polygon hierarchy to add to the triangulation.
 
 # Outputs
 - `new_tri::Triangulation`: The new triangulation, now containing `ghost_vertex_map` in the `ghost_vertex_map` field and `ghost_vertex_ranges` in the `ghost_vertex_ranges` field.
 """
-function replace_ghost_vertex_information(tri::Triangulation, ghost_vertex_map, ghost_vertex_ranges)
+function replace_ghost_vertex_information(tri::Triangulation, ghost_vertex_map, ghost_vertex_ranges, polygon_hierarchy = get_polygon_hierarchy(tri))
     points = get_points(tri)
     triangles = get_triangles(tri)
     boundary_nodes = get_boundary_nodes(tri)
@@ -376,7 +383,6 @@ function replace_ghost_vertex_information(tri::Triangulation, ghost_vertex_map,
     boundary_edge_map = get_boundary_edge_map(tri)
     ch = get_convex_hull(tri)
     representative_point_list = get_representative_point_list(tri)
-    polygon_hierarchy = get_polygon_hierarchy(tri)
     boundary_enricher = get_boundary_enricher(tri)
     cache = get_cache(tri)
     return Triangulation(
@@ -828,13 +834,11 @@ function constrained_triangulation!(tri::Triangulation, segments, boundary_nodes
     for e in each_edge(all_segments)
         add_segment!(new_tri, e; predicates, rng)
     end
-    new_tri_2 = replace_ghost_vertex_information(new_tri, ghost_vertex_map, ghost_vertex_ranges)
+    new_tri_2 = replace_ghost_vertex_information(new_tri, ghost_vertex_map, ghost_vertex_ranges, full_polygon_hierarchy)
     if !(isnothing(boundary_nodes) || !has_boundary_nodes(boundary_nodes)) && delete_holes
         delete_holes!(new_tri_2)
         add_boundary_information!(new_tri_2)
         add_ghost_triangles!(new_tri_2) # fix the ghost triangles 
-        polygon_hierarchy = get_polygon_hierarchy(new_tri_2)
-        copyto!(polygon_hierarchy, full_polygon_hierarchy)
     end
     return new_tri_2
 end
diff --git a/src/data_structures/convex_hull.jl b/src/data_structures/convex_hull.jl
index 2f7692529..74e5c82fd 100644
--- a/src/data_structures/convex_hull.jl
+++ b/src/data_structures/convex_hull.jl
@@ -9,7 +9,7 @@ Struct for representing a convex hull. See also [`convex_hull`](@ref).
 
 # Constructors 
     ConvexHull(points, vertices)
-    convex_hull(points; IntegerType=Int)
+    convex_hull(points; predicates=AdaptiveKernel(), IntegerType=Int)
 """
 struct ConvexHull{P, I}
     points::P
@@ -32,6 +32,12 @@ function Base.show(io::IO, m::MIME"text/plain", ch::ConvexHull)
 end
 Base.sizehint!(ch::ConvexHull, n) = Base.sizehint!(get_vertices(ch), n)
 
+function Base.copy(ch::ConvexHull)
+    p = get_points(ch)
+    v = get_vertices(ch)
+    return ConvexHull(copy(p), copy(v))
+end
+
 @doc """
     get_points(convex_hull::ConvexHull) -> Points
 
diff --git a/src/data_structures/mesh_refinement/boundary_enricher.jl b/src/data_structures/mesh_refinement/boundary_enricher.jl
index ca13df746..1b45ab05d 100644
--- a/src/data_structures/mesh_refinement/boundary_enricher.jl
+++ b/src/data_structures/mesh_refinement/boundary_enricher.jl
@@ -15,6 +15,8 @@ Base.:(==)(member₁::SmallAngleComplexMember, member₂::SmallAngleComplexMembe
 Base.show(io::IO, member::SmallAngleComplexMember) = print(io, "SmallAngleComplexMember with parent curve ", get_parent_curve(member), " and next edge ", get_next_edge(member))
 Base.show(io::IO, ::MIME"text/plain", member::SmallAngleComplexMember) = print(io, "SmallAngleComplexMember with parent curve ", get_parent_curve(member), " and next edge ", get_next_edge(member), ".")
 
+Base.copy(member::SmallAngleComplexMember) = member
+
 """
     replace_next_edge(member::SmallAngleComplexMember{I}, next_edge) where {I} -> SmallAngleComplexMember{I}
 
@@ -45,6 +47,8 @@ Base.:(==)(complex₁::SmallAngleComplex, complex₂::SmallAngleComplex) = get_a
 Base.show(io::IO, complex::SmallAngleComplex) = print(io, "SmallAngleComplex with apex ", get_apex(complex), " and ", length(get_members(complex)), " members")
 Base.show(io::IO, ::MIME"text/plain", complex::SmallAngleComplex) = print(io, "SmallAngleComplex with apex ", get_apex(complex), " and ", length(get_members(complex)), " members.")
 
+Base.copy(complex::SmallAngleComplex) = SmallAngleComplex(get_apex(complex), copy(get_members(complex)))
+
 """
     replace_next_edge!(complex::SmallAngleComplex, member_id, next_edge)
 
@@ -105,8 +109,8 @@ get_members(complex::SmallAngleComplex) = complex.members
 
 Returns a map from an apex vertex to a list of all curves that define a small angle complex associated with that apex vertex.
 """
-function get_small_angle_complexes(points, boundary_nodes, boundary_curves, segments = nothing; IntegerType = Int)
-    d = Dict{IntegerType, Vector{SmallAngleComplex{IntegerType}}}()
+function get_small_angle_complexes(points, boundary_nodes, boundary_curves, segments=nothing; IntegerType=Int)
+    d = Dict{IntegerType,Vector{SmallAngleComplex{IntegerType}}}()
     if has_multiple_curves(boundary_nodes)
         _get_small_angle_complexes_multiple_curves!(d, boundary_nodes, boundary_curves, IntegerType)
     elseif has_multiple_sections(boundary_nodes)
@@ -208,7 +212,7 @@ Segments are stored twice. The vertices associated with a vertex are sorted coun
 to define the coordinates.
 """
 function construct_segment_map(segments, points, ::Type{I}) where {I}
-    segment_map = Dict{I, Vector{I}}()
+    segment_map = Dict{I,Vector{I}}()
     for e in each_edge(segments)
         i, j = edge_vertices(e)
         iset = get!(Vector{I}, segment_map, i)
@@ -225,12 +229,12 @@ function construct_segment_map(segments, points, ::Type{I}) where {I}
         qx, qy = getxy(q)
         base = (qx - px, qy - py)
         sort!(
-            vertices, by = _vertex -> begin
+            vertices, by=_vertex -> begin
                 _q = get_point(points, _vertex)
                 _qx, _qy = getxy(_q)
                 next_base = (_qx - px, _qy - py)
                 return angle_between(base, next_base)
-            end, rev = false,
+            end, rev=false,
         )
     end
     return segment_map
@@ -251,12 +255,12 @@ function sort_members!(complex::SmallAngleComplex, points)
     qx, qy = getxy(q)
     base = (qx - px, qy - py)
     sort!(
-        members, by = member -> begin
+        members, by=member -> begin
             _q = get_point(points, get_next_edge(member))
             _qx, _qy = getxy(_q)
             next_base = (_qx - px, _qy - py)
             return angle_between(base, next_base)
-        end, rev = false,
+        end, rev=false,
     )
     return complex
 end
@@ -367,38 +371,59 @@ field will no longer aliased with the input `boundary_nodes`, although `points`
 The argument `n` is used in [`polygonise`](@ref) for filling out the boundary temporarily in order to construct the [`PolygonHierarchy`](@ref). The argument `coarse_n` defines the initial coarse discretisation 
 through [`coarse_discretisation!`](@ref); the default `n=0` means that the coarse discretisation will be performed until the maximum total variation of a subcurve is less than π/2.
 """
-struct BoundaryEnricher{P, B, C, I, BM, S, E}
+struct BoundaryEnricher{P,B,C,I,BM,S,E}
     points::P
     boundary_nodes::B
     segments::S
     boundary_curves::C
     polygon_hierarchy::PolygonHierarchy{I}
-    parent_map::Dict{E, I}
-    curve_index_map::Dict{I, I}
+    parent_map::Dict{E,I}
+    curve_index_map::Dict{I,I}
     boundary_edge_map::BM
     spatial_tree::BoundaryRTree{P}
     queue::Queue{I}
-    small_angle_complexes::Dict{I, Vector{SmallAngleComplex{I}}}
+    small_angle_complexes::Dict{I,Vector{SmallAngleComplex{I}}}
+end
+
+Base.copy(enricher::BoundaryEnricher) = enrcopy(enricher)
+enrcopy(::Nothing; kwargs...) = nothing
+
+function enrcopy(enricher::BoundaryEnricher; 
+        points = copy(get_points(enricher)), 
+        boundary_nodes  = copy(get_boundary_nodes(enricher)),
+        segments = copy(get_segments(enricher)),
+        boundary_curves = _plcopy.(get_boundary_curves(enricher); points),
+        polygon_hierarchy = copy(get_polygon_hierarchy(enricher)),
+        boundary_edge_map = copy(get_boundary_edge_map(enricher))) 
+    parent_map = copy(get_parent_map(enricher))
+    curve_index_map = copy(get_curve_index_map(enricher))
+    spatial_tree = copy(get_spatial_tree(enricher))
+    queue = copy(get_queue(enricher))
+    small_angle_complexes = copy(get_small_angle_complexes(enricher))
+    return BoundaryEnricher(points, boundary_nodes, segments, boundary_curves,
+        polygon_hierarchy, parent_map, curve_index_map,
+        boundary_edge_map, spatial_tree, queue, small_angle_complexes)
 end
+
 function BoundaryEnricher(
-        points::P, boundary_nodes::B, segments = nothing; n = 4096, coarse_n = 0,
-        IntegerType::Type{I} = Int,
-        EdgeType::Type{E} = isnothing(segments) ? NTuple{2, IntegerType} : (edge_type ∘ typeof)(segments),
-        EdgesType::Type{Es} = isnothing(segments) ? Set{EdgeType} : typeof(segments),
-    ) where {P, B, I, E, Es}
+    points::P, boundary_nodes::B, segments=nothing; n=4096, coarse_n=0,
+    IntegerType::Type{I}=Int,
+    EdgeType::Type{E}=isnothing(segments) ? NTuple{2,IntegerType} : (edge_type ∘ typeof)(segments),
+    EdgesType::Type{Es}=isnothing(segments) ? Set{EdgeType} : typeof(segments),
+) where {P,B,I,E,Es}
     boundary_curves, new_boundary_nodes = convert_boundary_curves!(points, boundary_nodes, I)
     polygon_hierarchy = construct_polygon_hierarchy(points, new_boundary_nodes, boundary_curves; IntegerType, n)
     return _construct_boundary_enricher(points, new_boundary_nodes, boundary_curves, polygon_hierarchy, segments, n, coarse_n, I, E, Es)
 end
-function _construct_boundary_enricher(points, boundary_nodes, boundary_curves, polygon_hierarchy, segments, n, coarse_n, ::Type{I}, ::Type{E}, ::Type{Es}) where {I, E, Es}
+function _construct_boundary_enricher(points, boundary_nodes, boundary_curves, polygon_hierarchy, segments, n, coarse_n, ::Type{I}, ::Type{E}, ::Type{Es}) where {I,E,Es}
     expand_bounds!(polygon_hierarchy, ε(number_type(points)))
-    coarse_discretisation!(points, boundary_nodes, boundary_curves; n = coarse_n)
+    coarse_discretisation!(points, boundary_nodes, boundary_curves; n=coarse_n)
     boundary_edge_map = construct_boundary_edge_map(boundary_nodes, I)
-    parent_map = Dict{E, I}()
-    curve_index_map = Dict{I, I}()
+    parent_map = Dict{E,I}()
+    curve_index_map = Dict{I,I}()
     spatial_tree = BoundaryRTree(points)
     queue = Queue{I}()
-    small_angle_complexes = get_small_angle_complexes(points, boundary_nodes, boundary_curves, segments; IntegerType = I)
+    small_angle_complexes = get_small_angle_complexes(points, boundary_nodes, boundary_curves, segments; IntegerType=I)
     _segments = isnothing(segments) ? Es() : segments
     enricher = BoundaryEnricher(points, boundary_nodes, _segments, boundary_curves, polygon_hierarchy, parent_map, curve_index_map, boundary_edge_map, spatial_tree, queue, small_angle_complexes)
     construct_parent_map!(enricher)
@@ -579,7 +604,7 @@ end
 
 Returns the parent of the edge `(i, j)` in `boundary_enricher`. If the edge is not in the parent map, then `$∅` is returned.
 """
-get_parent(boundary_enricher::BoundaryEnricher{P, B, C, I}, i, j) where {P, B, C, I} = get(get_parent_map(boundary_enricher), (i, j), I(∅))
+get_parent(boundary_enricher::BoundaryEnricher{P,B,C,I}, i, j) where {P,B,C,I} = get(get_parent_map(boundary_enricher), (i, j), I(∅))
 
 """
     set_parent!(boundary_enricher::BoundaryEnricher, i, j, k)
@@ -654,7 +679,7 @@ function _construct_parent_map_multiple_curves!(enricher::BoundaryEnricher)
     end
     return enricher
 end
-function _construct_parent_map_multiple_sections!(enricher::BoundaryEnricher, boundary_nodes = get_boundary_nodes(enricher), ctr = 1)
+function _construct_parent_map_multiple_sections!(enricher::BoundaryEnricher, boundary_nodes=get_boundary_nodes(enricher), ctr=1)
     ns = num_sections(boundary_nodes)
     for i in 1:ns
         section_nodes = get_boundary_nodes(boundary_nodes, i)
@@ -663,7 +688,7 @@ function _construct_parent_map_multiple_sections!(enricher::BoundaryEnricher, bo
     end
     return enricher
 end
-function _construct_parent_map_contiguous!(enricher::BoundaryEnricher, boundary_nodes = get_boundary_nodes(enricher), ctr = 1)
+function _construct_parent_map_contiguous!(enricher::BoundaryEnricher, boundary_nodes=get_boundary_nodes(enricher), ctr=1)
     n = num_boundary_edges(boundary_nodes)
     for i in 1:n
         u = get_boundary_nodes(boundary_nodes, i)
@@ -725,7 +750,7 @@ Returns `true` if the `curve_index`th curve in `enricher` is piecewise linear, a
     boundary_curves = get_boundary_curves(enricher)
     return is_piecewise_linear(boundary_curves, curve_index)
 end
-@inline function is_piecewise_linear(boundary_curves::C, curve_index) where {C <: Tuple}
+@inline function is_piecewise_linear(boundary_curves::C, curve_index) where {C<:Tuple}
     isempty(boundary_curves) && return true
     return eval_fnc_at_het_tuple_element(is_piecewise_linear, boundary_curves, curve_index)
 end
@@ -739,7 +764,7 @@ Returns `true` if the `curve_index`th curve in `enricher` is a [`LineSegment`](@
     boundary_curves = get_boundary_curves(enricher)
     return is_linear(boundary_curves, curve_index)
 end
-@inline function is_linear(boundary_curves::C, curve_index) where {C <: Tuple}
+@inline function is_linear(boundary_curves::C, curve_index) where {C<:Tuple}
     isempty(boundary_curves) && return true
     return eval_fnc_at_het_tuple_element(is_linear, boundary_curves, curve_index)
 end
@@ -753,7 +778,7 @@ Returns the inverse of the `curve_index`th curve at `q`.
     boundary_curves = get_boundary_curves(enricher)
     return get_inverse(boundary_curves, curve_index, q)
 end
-@inline function get_inverse(boundary_curves::C, curve_index, q) where {C <: Tuple}
+@inline function get_inverse(boundary_curves::C, curve_index, q) where {C<:Tuple}
     return eval_fnc_at_het_tuple_element_with_arg(get_inverse, boundary_curves, (q,), curve_index)
 end
 
@@ -766,7 +791,7 @@ Returns the equivariation split of the `curve_index`th curve between `t₁` and
     boundary_curves = get_boundary_curves(enricher)
     return get_equivariation_split(boundary_curves, curve_index, t₁, t₂)
 end
-@inline function get_equivariation_split(boundary_curves::C, curve_index, t₁, t₂) where {C <: Tuple}
+@inline function get_equivariation_split(boundary_curves::C, curve_index, t₁, t₂) where {C<:Tuple}
     return eval_fnc_at_het_tuple_element_with_arg(get_equivariation_split, boundary_curves, (t₁, t₂), curve_index)
 end
 
@@ -779,7 +804,7 @@ Returns the equidistant split of the `curve_index`th curve between `t₁` and `t
     boundary_curves = get_boundary_curves(enricher)
     return get_equidistant_split(boundary_curves, curve_index, t₁, t₂)
 end
-@inline function get_equidistant_split(boundary_curves::C, curve_index, t₁, t₂) where {C <: Tuple}
+@inline function get_equidistant_split(boundary_curves::C, curve_index, t₁, t₂) where {C<:Tuple}
     return eval_fnc_at_het_tuple_element_with_arg(get_equidistant_split, boundary_curves, (t₁, t₂), curve_index)
 end
 
@@ -793,7 +818,7 @@ Returns the `curve_index`th boundary curve at `t`.
     boundary_curves = get_boundary_curves(enricher)
     return eval_boundary_curve(boundary_curves, curve_index, t)
 end
-@inline function eval_boundary_curve(boundary_curves::C, curve_index, t) where {C <: Tuple}
+@inline function eval_boundary_curve(boundary_curves::C, curve_index, t) where {C<:Tuple}
     return eval_fnc_in_het_tuple(boundary_curves, t, curve_index)
 end
 
@@ -815,10 +840,10 @@ See the documentation for more information about these choices.
     return point_position_relative_to_curve(kernel, boundary_curves, curve_index, p)
 end
 @inline point_position_relative_to_curve(enricher::BoundaryEnricher, curve_index, p) = point_position_relative_to_curve(AdaptiveKernel(), enricher, curve_index, p)
-@inline function point_position_relative_to_curve(kernel::AbstractPredicateKernel, boundary_curves::C, curve_index, p) where {C <: Tuple}
+@inline function point_position_relative_to_curve(kernel::AbstractPredicateKernel, boundary_curves::C, curve_index, p) where {C<:Tuple}
     return eval_fnc_at_het_tuple_element_with_arg_and_prearg(point_position_relative_to_curve, boundary_curves, kernel, (p,), curve_index)
 end
-@inline point_position_relative_to_curve(boundary_curves::C, curve_index, p) where {C <: Tuple} = point_position_relative_to_curve(AdaptiveKernel(), boundary_curves, curve_index, p)
+@inline point_position_relative_to_curve(boundary_curves::C, curve_index, p) where {C<:Tuple} = point_position_relative_to_curve(AdaptiveKernel(), boundary_curves, curve_index, p)
 
 """
     angle_between(enricher::BoundaryEnricher, curve_index1, curve_index2) -> Float64 
@@ -874,7 +899,7 @@ Fills out a set of points for a curve-bounded domain for use with [`PolygonHiera
     If the boundary is not curve bounded, then `new_points` and `new_boundary_nodes` remain aliased 
     with the input `points` and `boundary_nodes`.
 """
-function polygonise(points, boundary_nodes, boundary_curves; n = 4096)
+function polygonise(points, boundary_nodes, boundary_curves; n=4096)
     new_points = deepcopy(points)
     new_boundary_nodes = deepcopy(boundary_nodes)
     coarse_discretisation!(new_points, new_boundary_nodes, boundary_curves; n)
@@ -923,16 +948,19 @@ end
 
 Updates the fields of `enricher` after splitting a boundary edge `(i, j)` at the `r`th vertex. The `update_boundary_nodes` argument
 can be used to avoid inserting an additional boundary node when `boundary_nodes` was already updated somewhere else (e.g., we need this for 
-mesh refinement which already updates the `boundary_nodes` which is aliased with the same field in the enricher).
+mesh refinement which already updates the `boundary_nodes` which is aliased with the same field in the enricher). Note that not updating the 
+boundary nodes also implies that the `boundary_edge_map` will not be updated.
 """
-function split_boundary_edge!(enricher::BoundaryEnricher, i, j, r, update_boundary_nodes = Val(true))
+function split_boundary_edge!(enricher::BoundaryEnricher, i, j, r, update_boundary_nodes=Val(true))
     boundary_nodes = get_boundary_nodes(enricher)
     boundary_edge_map = get_boundary_edge_map(enricher)
     spatial_tree = get_spatial_tree(enricher)
-    pos = get_boundary_edge_map(enricher, i, j)
-    new_pos = (pos[1], pos[2] + 1)
-    is_true(update_boundary_nodes) && insert_boundary_node!(boundary_nodes, new_pos, r)
-    split_boundary_edge_map!(boundary_edge_map, boundary_nodes, pos, i, j)
+    if is_true(update_boundary_nodes)
+        pos = get_boundary_edge_map(enricher, i, j)
+        new_pos = (pos[1], pos[2] + 1)
+        insert_boundary_node!(boundary_nodes, new_pos, r)
+        split_boundary_edge_map!(boundary_edge_map, boundary_nodes, pos, i, j)
+    end
     split_edge!(spatial_tree, i, j, r)
     update_parent_map!(enricher, i, j, r)
     return enricher
@@ -945,7 +973,7 @@ Updates the fields of `enricher` after splitting an interior segment `(i, j)` at
 The `update_segments` argument can be used to avoid inserting an additional segment when `segments` was already updated somewhere else (e.g., we need this for
 mesh refinement which already updates the `interior_segments` which is aliased with the `segments` field in the enricher).
 """
-function split_interior_segment!(enricher::BoundaryEnricher, i, j, r, update_segments = Val(true))
+function split_interior_segment!(enricher::BoundaryEnricher, i, j, r, update_segments=Val(true))
     segments = get_segments(enricher)
     spatial_tree = get_spatial_tree(enricher)
     E = edge_type(segments)
@@ -971,7 +999,7 @@ The `is_interior` argument can be used to specify whether the edge is an interio
 
 See also [`split_boundary_edge!`](@ref) and [`split_interior_segment!`](@ref).
 """
-function split_edge!(enricher::BoundaryEnricher, i, j, r, update_boundary_nodes = Val(true), update_segments = Val(true), is_interior = is_segment(enricher, i, j))
+function split_edge!(enricher::BoundaryEnricher, i, j, r, update_boundary_nodes=Val(true), update_segments=Val(true), is_interior=is_segment(enricher, i, j))
     if is_interior
         split_interior_segment!(enricher, i, j, r, update_segments)
     else
diff --git a/src/data_structures/mesh_refinement/curves/abstract.jl b/src/data_structures/mesh_refinement/curves/abstract.jl
index 0d9ba0a9c..29f5edb76 100644
--- a/src/data_structures/mesh_refinement/curves/abstract.jl
+++ b/src/data_structures/mesh_refinement/curves/abstract.jl
@@ -1028,3 +1028,5 @@ Returns an [`AbstractParametricCurve`](@ref) that reverses the orientation of `c
 function Base.reverse(c::AbstractParametricCurve)
     return _reverse(c)
 end
+
+_plcopy(c::AbstractParametricCurve; kwargs...) = copy(c) # for handling PiecewiseLinear
\ No newline at end of file
diff --git a/src/data_structures/mesh_refinement/curves/beziercurve.jl b/src/data_structures/mesh_refinement/curves/beziercurve.jl
index d8dffa6be..81858b936 100644
--- a/src/data_structures/mesh_refinement/curves/beziercurve.jl
+++ b/src/data_structures/mesh_refinement/curves/beziercurve.jl
@@ -45,6 +45,15 @@ function Base.:(==)(b₁::BezierCurve, b₂::BezierCurve)
     return true
 end
 
+function Base.copy(b::BezierCurve)
+    return BezierCurve(
+        copy(b.control_points),
+        copy(b.cache),
+        copy(b.lookup_table),
+        copy(b.orientation_markers)
+    )
+end
+
 function BezierCurve(control_points::Vector{NTuple{2,Float64}}; lookup_steps=5000, kwargs...)
     cache = similar(control_points) # will be copyto! later
     lookup_table = similar(control_points, lookup_steps)
diff --git a/src/data_structures/mesh_refinement/curves/bspline.jl b/src/data_structures/mesh_refinement/curves/bspline.jl
index a1e7fa9fe..dafaf6dbb 100644
--- a/src/data_structures/mesh_refinement/curves/bspline.jl
+++ b/src/data_structures/mesh_refinement/curves/bspline.jl
@@ -45,6 +45,16 @@ function Base.:(==)(b₁::BSpline, b₂::BSpline)
     return true
 end
 
+function Base.copy(b::BSpline)
+    return BSpline(
+        copy(b.control_points),
+        copy(b.knots),
+        copy(b.cache),
+        copy(b.lookup_table),
+        copy(b.orientation_markers)
+    )
+end
+
 function BSpline(control_points::Vector{NTuple{2,Float64}}; degree=3, lookup_steps=5000, kwargs...)
     nc = length(control_points)
     @assert degree ≥ 1 "Degree must be at least 1, got $degree."
diff --git a/src/data_structures/mesh_refinement/curves/catmullromspline.jl b/src/data_structures/mesh_refinement/curves/catmullromspline.jl
index d37927b19..4523d24a2 100644
--- a/src/data_structures/mesh_refinement/curves/catmullromspline.jl
+++ b/src/data_structures/mesh_refinement/curves/catmullromspline.jl
@@ -46,6 +46,8 @@ struct CatmullRomSplineSegment <: AbstractParametricCurve
     p₂::NTuple{2,Float64}
 end
 
+Base.copy(c::CatmullRomSplineSegment) = c
+
 function (c::CatmullRomSplineSegment)(t)
     if iszero(t)
         return c.p₁
@@ -203,6 +205,21 @@ struct CatmullRomSpline <: AbstractParametricCurve
     orientation_markers::Vector{Float64}
 end
 
+function Base.copy(c::CatmullRomSpline)
+    return CatmullRomSpline(
+        copy(c.control_points),
+        copy(c.knots),
+        copy(c.lookup_table),
+        c.alpha,
+        c.tension,
+        c.left,
+        c.right,
+        copy(c.lengths),
+        copy(c.segments),
+        copy(c.orientation_markers)
+    )
+end
+
 function _reverse(c::CatmullRomSpline)
     return CatmullRomSpline(
         reverse(c.control_points),
diff --git a/src/data_structures/mesh_refinement/curves/circulararc.jl b/src/data_structures/mesh_refinement/curves/circulararc.jl
index cf2c3b447..f6d29704e 100644
--- a/src/data_structures/mesh_refinement/curves/circulararc.jl
+++ b/src/data_structures/mesh_refinement/curves/circulararc.jl
@@ -44,6 +44,8 @@ function Base.:(==)(c₁::CircularArc, c₂::CircularArc)
     return true
 end
 
+Base.copy(c::CircularArc) = c
+
 function CircularArc(p, q, c; positive=true)
     px, py = getxy(p)
     qx, qy = getxy(q)
diff --git a/src/data_structures/mesh_refinement/curves/ellipticalarc.jl b/src/data_structures/mesh_refinement/curves/ellipticalarc.jl
index 45fc7dc89..3f211ab32 100644
--- a/src/data_structures/mesh_refinement/curves/ellipticalarc.jl
+++ b/src/data_structures/mesh_refinement/curves/ellipticalarc.jl
@@ -42,6 +42,8 @@ function Base.:(==)(e₁::EllipticalArc, e₂::EllipticalArc)
     return true
 end
 
+Base.copy(e::EllipticalArc) = e
+
 function EllipticalArc(p, q, c, α, β, θ°; positive=true)
     px, py = getxy(p)
     qx, qy = getxy(q)
diff --git a/src/data_structures/mesh_refinement/curves/linesegment.jl b/src/data_structures/mesh_refinement/curves/linesegment.jl
index cd68b996a..cc884e654 100644
--- a/src/data_structures/mesh_refinement/curves/linesegment.jl
+++ b/src/data_structures/mesh_refinement/curves/linesegment.jl
@@ -21,6 +21,8 @@ struct LineSegment <: AbstractParametricCurve # line segment
 end
 Base.:(==)(L₁::LineSegment, L₂::LineSegment) = L₁.first == L₂.first && L₁.last == L₂.last
 
+Base.copy(L::LineSegment) = L
+
 function LineSegment(p₀, p₁)
     return LineSegment(p₀, p₁, dist(p₀, p₁))
 end
diff --git a/src/data_structures/mesh_refinement/curves/piecewiselinear.jl b/src/data_structures/mesh_refinement/curves/piecewiselinear.jl
index 2684c35a7..f536ebf87 100644
--- a/src/data_structures/mesh_refinement/curves/piecewiselinear.jl
+++ b/src/data_structures/mesh_refinement/curves/piecewiselinear.jl
@@ -22,6 +22,11 @@ end
 Base.show(io::IO, ::PiecewiseLinear) = print(io, "PiecewiseLinear()")
 Base.show(io::IO, ::MIME"text/plain", L::PiecewiseLinear) = Base.show(io, L)
 Base.:(==)(L1::PiecewiseLinear, L2::PiecewiseLinear) = get_points(L1) == get_points(L2) && get_boundary_nodes(L1) == get_boundary_nodes(L2)
+
+Base.copy(L::PiecewiseLinear) = _plcopy(L) # so we can check the aliasing when we copy from a triangulation
+_plcopy(L::PiecewiseLinear; points=copy(get_points(L)), boundary_nodes=copy(get_boundary_nodes(L))) = PiecewiseLinear(points, boundary_nodes)
+
+
 is_piecewise_linear(::PiecewiseLinear) = true
 get_points(pl::PiecewiseLinear) = pl.points
 get_boundary_nodes(pl::PiecewiseLinear) = pl.boundary_nodes
diff --git a/src/data_structures/queue/queue.jl b/src/data_structures/queue/queue.jl
index 1961f5c6f..c49ac25cd 100644
--- a/src/data_structures/queue/queue.jl
+++ b/src/data_structures/queue/queue.jl
@@ -14,6 +14,8 @@ end
 Queue{T}() where {T} = Queue{T}(T[])
 Base.:(==)(q1::Queue, q2::Queue) = q1.data == q2.data
 
+Base.copy(queue::Queue) = Queue(copy(queue.data))
+
 """
     isempty(queue::Queue) -> Bool 
 
diff --git a/src/data_structures/representative_coordinates/representative_coordinates.jl b/src/data_structures/representative_coordinates/representative_coordinates.jl
index c3f34e6ac..82bd865b2 100644
--- a/src/data_structures/representative_coordinates/representative_coordinates.jl
+++ b/src/data_structures/representative_coordinates/representative_coordinates.jl
@@ -32,6 +32,10 @@ function Base.convert(::Type{RepresentativeCoordinates{I, T}}, c::Representative
     return c
 end
 
+function Base.copy(c::RepresentativeCoordinates)
+    return RepresentativeCoordinates(getx(c), gety(c), getn(c))
+end
+
 """
     getx(c::RepresentativeCoordinates) -> Number
 
diff --git a/src/data_structures/trees/rtree.jl b/src/data_structures/trees/rtree.jl
index fc0e09019..8f3dddc61 100644
--- a/src/data_structures/trees/rtree.jl
+++ b/src/data_structures/trees/rtree.jl
@@ -45,6 +45,8 @@ struct BoundingInterval <: AbstractBoundingShape
 end
 Base.show(io::IO, I::BoundingInterval) = print(io, "[$(I.a) .. $(I.b)]")
 
+Base.copy(I::BoundingInterval) = I
+
 """
     InvalidBoundingInterval
 
@@ -138,6 +140,8 @@ struct BoundingBox <: AbstractBoundingShape
 end
 Base.show(io::IO, r::BoundingBox) = print(io, "[$(r.x.a), $(r.x.b)] × [$(r.y.a), $(r.y.b)]")
 
+Base.copy(r::BoundingBox) = r
+
 """
     InvalidBoundingBox
 
@@ -279,7 +283,7 @@ function bounding_box(points)
     ymin = Inf
     ymax = -Inf
     for p in each_point(points)
-        px, py = getxy(p)
+        px, py = _getxy(p)
         xmin = min(xmin, px)
         xmax = max(xmax, px)
         ymin = min(ymin, py)
@@ -331,6 +335,8 @@ struct DiametralBoundingBox
     edge::NTuple{2, Int}
 end
 
+Base.copy(bbox::DiametralBoundingBox) = bbox
+
 """
     bounding_box(points, i, j) -> DiametralBoundingBox
 
@@ -508,6 +514,24 @@ function Base.:(==)(leaf1::Leaf, leaf2::Leaf)
     return true
 end
 
+function Base.deepcopy_internal(leaf::Leaf{Branch}, dict::IdDict) where {Branch}
+    haskey(dict, leaf) && return dict[leaf]
+    new_parent = has_parent(leaf) ? Base.deepcopy_internal(get_parent(leaf), dict) : nothing
+    bbox = get_bounding_box(leaf)
+    cchildren = copy(get_children(leaf))
+    new_leaf = Leaf{Branch}(new_parent, bbox, cchildren)
+    dict[leaf] = new_leaf
+    return new_leaf
+end
+
+function Base.copy(leaf::Leaf{B}) where {B}
+    parent = get_parent(leaf)
+    cparent = isnothing(parent) ? nothing : copy(parent)
+    bbox = get_bounding_box(leaf)
+    cchildren = copy(get_children(leaf))
+    return Leaf{B}(cparent, bbox, cchildren)
+end
+
 """
     mutable struct Branch <: AbstractNode
 
@@ -540,6 +564,24 @@ function Base.:(==)(branch1::Branch, branch2::Branch)
     return true
 end
 
+function Base.deepcopy_internal(branch::Branch, dict::IdDict)
+    haskey(dict, branch) && return dict[branch]
+    new_parent = has_parent(branch) ? Base.deepcopy_internal(get_parent(branch), dict) : nothing
+    bbox = get_bounding_box(branch)
+    cchildren = copy(get_children(branch))
+    new_branch = Branch(new_parent, bbox, cchildren, get_level(branch))
+    dict[branch] = new_branch
+    return new_branch
+end
+
+function Base.copy(branch::Branch)
+    parent = get_parent(branch)
+    cparent = isnothing(parent) ? nothing : copy(parent)
+    bbox = get_bounding_box(branch)
+    cchildren = copy(get_children(branch))
+    return Branch(cparent, bbox, cchildren, get_level(branch))
+end
+
 function Base.show(io::IO, ::MIME"text/plain", leaf::Leaf)
     parent = get_parent(leaf)
     nchildren = num_children(leaf)
@@ -615,6 +657,13 @@ struct NodeCache{Node, Child} # similar to why we use TriangulationCache. Think
     end
 end
 
+Base.:(==)(cache1::NodeCache, cache2::NodeCache) = get_size_limit(cache1) == get_size_limit(cache2) && cache1.cache == cache2.cache
+
+function Base.copy(cache::NodeCache{Node, Child}) where {Node, Child}
+    # Just generates a new cache
+    return NodeCache{Node, Child}(get_size_limit(cache))
+end
+
 """
     BranchCache
 
@@ -713,6 +762,10 @@ end
 RTreeIntersectionCache() = RTreeIntersectionCache(Vector{Int}(), BitVector())
 Base.sizehint!(cache::RTreeIntersectionCache, n) = (sizehint!(get_node_indices(cache), n); sizehint!(get_need_tests(cache), n))
 
+Base.:(==)(cache1::RTreeIntersectionCache, cache2::RTreeIntersectionCache) = get_node_indices(cache1) == get_node_indices(cache2) && get_need_tests(cache1) == get_need_tests(cache2)
+
+Base.copy(cache::RTreeIntersectionCache) = RTreeIntersectionCache() # Just generates a new cache
+
 """
     get_node_indices(cache::RTreeIntersectionCache) -> Vector{Int}
 
@@ -761,32 +814,45 @@ mutable struct RTree # linear
     @const free_cache::BitVector
     @const detached_cache::Vector{Union{Branch, Leaf{Branch}}}
     @const intersection_cache::NTuple{2, RTreeIntersectionCache}
-    function RTree(; size_limit = 100, fill_factor = 0.7) # https://en.wikipedia.org/wiki/R-tree: "however best performance has been experienced with a minimum fill of 30%–40%)
-        branch_cache = BranchCache(size_limit)
-        twig_cache = TwigCache(size_limit)
-        leaf_cache = LeafCache(size_limit)
-        root = spawn_node!(leaf_cache)
-        num_elements = 0
-        free_cache = BitVector()
-        sizehint!(free_cache, size_limit)
-        detached_cache = Vector{Union{Branch, Leaf{Branch}}}()
-        sizehint!(detached_cache, size_limit)
-        cache1, cache2 = RTreeIntersectionCache(), RTreeIntersectionCache()
-        sizehint!(cache1, ceil(Int, log2(size_limit)))
-        sizehint!(cache2, ceil(Int, log2(size_limit)))
-        return new(
-            root,
-            num_elements,
-            branch_cache,
-            twig_cache,
-            leaf_cache,
-            fill_factor,
-            free_cache,
-            detached_cache,
-            (cache1, cache2),
-        )
-    end
 end
+function RTree(; size_limit = 100, fill_factor = 0.7) # https://en.wikipedia.org/wiki/R-tree: "however best performance has been experienced with a minimum fill of 30%–40%)
+    branch_cache = BranchCache(size_limit)
+    twig_cache = TwigCache(size_limit)
+    leaf_cache = LeafCache(size_limit)
+    root = spawn_node!(leaf_cache)
+    num_elements = 0
+    free_cache = BitVector()
+    sizehint!(free_cache, size_limit)
+    detached_cache = Vector{Union{Branch, Leaf{Branch}}}()
+    sizehint!(detached_cache, size_limit)
+    cache1, cache2 = RTreeIntersectionCache(), RTreeIntersectionCache()
+    sizehint!(cache1, ceil(Int, log2(size_limit)))
+    sizehint!(cache2, ceil(Int, log2(size_limit)))
+    return RTree(
+        root,
+        num_elements,
+        branch_cache,
+        twig_cache,
+        leaf_cache,
+        fill_factor,
+        free_cache,
+        detached_cache,
+        (cache1, cache2),
+    )
+end
+
+function Base.copy(tree::RTree)
+    root = copy(get_root(tree))
+    branch_cache = copy(get_branch_cache(tree))
+    twig_cache = copy(get_twig_cache(tree))
+    leaf_cache = copy(get_leaf_cache(tree))
+    fill_factor = get_fill_factor(tree)
+    free_cache = copy(get_free_cache(tree))
+    detached_cache = copy(get_detached_cache(tree))
+    intersection_cache = copy.(get_intersection_cache(tree))
+    return RTree(root, num_elements(tree), branch_cache, twig_cache, leaf_cache, fill_factor, free_cache, detached_cache, intersection_cache)
+end
+
 function Base.:(==)(tree1::RTree, tree2::RTree)
     num_elements(tree1) ≠ num_elements(tree2) && return false
     get_root(tree1) ≠ get_root(tree2) && return false
@@ -1071,6 +1137,12 @@ struct BoundaryRTree{P}
 end
 BoundaryRTree(points) = BoundaryRTree(RTree(), points)
 
+function Base.copy(tree::BoundaryRTree)
+    ctree = copy(tree.tree)
+    cpoints = copy(tree.points)
+    return BoundaryRTree(ctree, cpoints)
+end
+
 function Base.:(==)(tree1::BoundaryRTree, tree2::BoundaryRTree)
     tree1.points ≠ tree2.points && return false
     tree1.tree ≠ tree2.tree && return false
diff --git a/src/data_structures/triangulation/adjacent.jl b/src/data_structures/triangulation/adjacent.jl
index f66c4031d..e9a63e4e3 100644
--- a/src/data_structures/triangulation/adjacent.jl
+++ b/src/data_structures/triangulation/adjacent.jl
@@ -23,6 +23,8 @@ function Base.show(io::IO, m::MIME"text/plain", adj::Adjacent{I, E}) where {I, E
 end
 Base.sizehint!(adj::Adjacent, n) = sizehint!(get_adjacent(adj), n)
 
+Base.copy(adj::Adjacent) = Adjacent(copy(get_adjacent(adj)))
+
 """
     get_adjacent(adj::Adjacent) -> Dict 
 
diff --git a/src/data_structures/triangulation/adjacent2vertex.jl b/src/data_structures/triangulation/adjacent2vertex.jl
index c29610c25..b9621311a 100644
--- a/src/data_structures/triangulation/adjacent2vertex.jl
+++ b/src/data_structures/triangulation/adjacent2vertex.jl
@@ -24,6 +24,8 @@ function Base.show(io::IO, m::MIME"text/plain", adj2v::Adjacent2Vertex{I, Es}) w
 end
 Base.sizehint!(adj2v::Adjacent2Vertex, n) = Base.sizehint!(get_adjacent2vertex(adj2v), n)
 
+Base.copy(adj2v::Adjacent2Vertex) = Adjacent2Vertex(copy(get_adjacent2vertex(adj2v)))
+
 """
     get_adjacent2vertex(adj2v::Adjacent2Vertex) -> Dict
 
diff --git a/src/data_structures/triangulation/graph.jl b/src/data_structures/triangulation/graph.jl
index 4bd8e82de..666be051b 100644
--- a/src/data_structures/triangulation/graph.jl
+++ b/src/data_structures/triangulation/graph.jl
@@ -54,6 +54,13 @@ function Base.sizehint!(graph::Graph, n1, n2, n3)
     return graph
 end
 
+function Base.copy(graph::Graph)
+    V = get_vertices(graph)
+    E = get_edges(graph)
+    N = get_neighbours(graph)
+    return Graph(copy(V), copy(E), copy(N))
+end
+
 """
     get_vertices(graph::Graph) -> Set{Vertex}
 
diff --git a/src/data_structures/triangulation/methods/weights.jl b/src/data_structures/triangulation/methods/weights.jl
index 38909e14e..893cdc9a5 100644
--- a/src/data_structures/triangulation/methods/weights.jl
+++ b/src/data_structures/triangulation/methods/weights.jl
@@ -27,6 +27,8 @@ ZeroWeight() = ZeroWeight{Float64}()
 get_weight(weights::ZeroWeight{T}, i) where {T} = is_point3(i) ? i[3] : zero(T)
 get_weight(weights::ZeroWeight{T}, i::Integer) where {T} = zero(T)
 
+Base.copy(w::ZeroWeight) = w
+
 """
     add_weight!(tri::Triangulation, w)
 
diff --git a/src/data_structures/triangulation/triangulation.jl b/src/data_structures/triangulation/triangulation.jl
index 7da4a6718..06e0d9916 100644
--- a/src/data_structures/triangulation/triangulation.jl
+++ b/src/data_structures/triangulation/triangulation.jl
@@ -108,6 +108,36 @@ struct Triangulation{P,T,BN,W,I,E,Es,BC,BEM,GVM,GVR,BPL,C,BE}
     boundary_enricher::BE
     cache::C
 end
+
+function Base.copy(tri::Triangulation)
+    points = copy(get_points(tri))
+    triangles = copy(get_triangles(tri))
+    boundary_nodes = copy(get_boundary_nodes(tri))
+    interior_segments = copy(get_interior_segments(tri))
+    all_segments = copy(get_all_segments(tri))
+    weights = copy(get_weights(tri))
+    adjacent = copy(get_adjacent(tri))
+    adjacent2vertex = copy(get_adjacent2vertex(tri))
+    graph = copy(get_graph(tri))
+    boundary_curves = _plcopy.(get_boundary_curves(tri); points)
+    boundary_edge_map = _bemcopy(get_boundary_edge_map(tri); boundary_nodes)
+    ghost_vertex_map = copy(get_ghost_vertex_map(tri))
+    ghost_vertex_ranges = copy(get_ghost_vertex_ranges(tri))
+    convex_hull = copy(get_convex_hull(tri))
+    representative_point_list = copy(get_representative_point_list(tri))
+    polygon_hierarchy = copy(get_polygon_hierarchy(tri))
+    boundary_enricher = enrcopy(get_boundary_enricher(tri); points, 
+        boundary_nodes, segments=interior_segments,boundary_curves,
+        polygon_hierarchy,boundary_edge_map) 
+    cache = _copy_cache(get_cache(tri); weights) 
+    return Triangulation(
+        points, triangles, boundary_nodes, interior_segments, all_segments, weights,
+        adjacent, adjacent2vertex, graph, boundary_curves, boundary_edge_map,
+        ghost_vertex_map, ghost_vertex_ranges, convex_hull, representative_point_list,
+        polygon_hierarchy, boundary_enricher, cache
+    )
+end
+
 function Base.show(io::IO, ::MIME"text/plain", tri::Triangulation)
     println(io, "Delaunay Triangulation.")
     println(io, "   Number of vertices: ", num_solid_vertices(tri))
diff --git a/src/data_structures/triangulation/triangulation_cache.jl b/src/data_structures/triangulation/triangulation_cache.jl
index bcb9f2737..2c1786ad6 100644
--- a/src/data_structures/triangulation/triangulation_cache.jl
+++ b/src/data_structures/triangulation/triangulation_cache.jl
@@ -52,6 +52,21 @@ function TriangulationCache()
     return TriangulationCache(nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing)
 end
 
+Base.copy(cache::EmptyTriangulationCache) = cache
+Base.copy(cache::TriangulationCache) = _copy_cache(cache)
+function _copy_cache(cache::TriangulationCache; weights = copy(get_weights(get_triangulation(cache))))
+    # Doesn't actually copy the values in the cache. Just generates a new cache.
+    tri = get_triangulation(cache)
+    points = copy(get_points(tri))
+    I = integer_type(tri)
+    E = edge_type(tri)
+    V = triangle_type(tri)
+    Es = edges_type(tri)
+    T = triangles_type(tri)
+    new_cache = _build_cache(points, I, E, V, Es, T, weights, Val(true))
+    return new_cache
+end
+
 """
     get_triangulation(cache::TriangulationCache) -> Triangulation
 
diff --git a/src/data_structures/voronoi.jl b/src/data_structures/voronoi.jl
index 9c85d79b1..0580918f5 100644
--- a/src/data_structures/voronoi.jl
+++ b/src/data_structures/voronoi.jl
@@ -49,6 +49,35 @@ struct VoronoiTessellation{Tr <: Triangulation, P, I, T, S, E}
     boundary_polygons::Set{I}
 end
 
+function Base.copy(vorn::VoronoiTessellation)
+    return VoronoiTessellation(
+        copy(get_triangulation(vorn)),
+        copy(get_generators(vorn)),
+        copy(get_polygon_points(vorn)),
+        copy(get_polygons(vorn)),
+        copy(get_circumcenter_to_triangle(vorn)),
+        copy(get_triangle_to_circumcenter(vorn)),
+        copy(get_unbounded_polygons(vorn)),
+        copy(get_cocircular_circumcenters(vorn)),
+        copy(get_adjacent(vorn)),
+        copy(get_boundary_polygons(vorn))
+    )
+end
+
+function Base.:(==)(vorn1::VoronoiTessellation, vorn2::VoronoiTessellation)
+    get_triangulation(vorn1) ≠ get_triangulation(vorn2) && return false 
+    get_generators(vorn1) ≠ get_generators(vorn2) && return false
+    get_polygon_points(vorn1) ≠ get_polygon_points(vorn2) && return false
+    get_polygons(vorn1) ≠ get_polygons(vorn2) && return false
+    get_circumcenter_to_triangle(vorn1) ≠ get_circumcenter_to_triangle(vorn2) && return false
+    get_triangle_to_circumcenter(vorn1) ≠ get_triangle_to_circumcenter(vorn2) && return false
+    get_unbounded_polygons(vorn1) ≠ get_unbounded_polygons(vorn2) && return false
+    get_cocircular_circumcenters(vorn1) ≠ get_cocircular_circumcenters(vorn2) && return false
+    get_adjacent(vorn1) ≠ get_adjacent(vorn2) && return false
+    get_boundary_polygons(vorn1) ≠ get_boundary_polygons(vorn2) && return false
+    return true
+end
+
 """
     get_triangulation(vorn::VoronoiTessellation) -> Triangulation 
 
diff --git a/src/geometric_primitives/boundary_nodes.jl b/src/geometric_primitives/boundary_nodes.jl
index 85109fc16..76b2a5558 100644
--- a/src/geometric_primitives/boundary_nodes.jl
+++ b/src/geometric_primitives/boundary_nodes.jl
@@ -390,6 +390,9 @@ Dict{Tuple{Int64, Int64}, Tuple{Tuple{Int64, Int64}, Int64}} with 12 entries:
 end
 function _construct_boundary_edge_map_multiple_curves(boundary_nodes, IntegerType::Type{I}, EdgeType::Type{E}) where {I, E}
     dict = Dict{E, Tuple{NTuple{2, I}, I}}()
+    return _construct_boundary_edge_map!(dict, boundary_nodes)
+end
+function _construct_boundary_edge_map!(dict::Dict{E, Tuple{NTuple{2, I}, I}}, boundary_nodes) where {E, I}
     nc = num_curves(boundary_nodes)
     for m in 1:nc
         bn_m = get_boundary_nodes(boundary_nodes, m)
@@ -409,6 +412,9 @@ function _construct_boundary_edge_map_multiple_curves(boundary_nodes, IntegerTyp
 end
 function _construct_boundary_edge_map_multiple_sections(boundary_nodes, IntegerType::Type{I}, EdgeType::Type{E}) where {I, E}
     dict = Dict{E, NTuple{2, I}}()
+    return _construct_boundary_edge_map!(dict, boundary_nodes)
+end
+function _construct_boundary_edge_map!(dict::Dict{E, NTuple{2, I}}, boundary_nodes) where {E, I} 
     ns = num_sections(boundary_nodes)
     for n in 1:ns
         bn_n = get_boundary_nodes(boundary_nodes, n)
@@ -424,6 +430,9 @@ function _construct_boundary_edge_map_multiple_sections(boundary_nodes, IntegerT
 end
 function _construct_boundary_edge_map_contiguous(boundary_nodes::A, IntegerType::Type{I}, EdgeType::Type{E}) where {A, I, E}
     dict = Dict{E, Tuple{A, I}}()
+    return _construct_boundary_edge_map!(dict, boundary_nodes)
+end
+function _construct_boundary_edge_map!(dict::Dict{E, Tuple{A, I}}, boundary_nodes) where {A, I, E}
     ne = num_boundary_edges(boundary_nodes)
     for ℓ in 1:ne
         u = get_boundary_nodes(boundary_nodes, ℓ)
@@ -434,6 +443,16 @@ function _construct_boundary_edge_map_contiguous(boundary_nodes::A, IntegerType:
     return dict
 end
 
+function _bemcopy(dict::Dict{E, Tuple{A, I}}; boundary_nodes) where {A, I, E} # so that boundary_nodes remains aliased with this dict when copying in Triangulation
+    new_dict = Dict{E, Tuple{A, I}}()
+    for (e, (pos, ℓ)) in dict
+        new_dict[e] = (boundary_nodes, ℓ)
+    end
+    return new_dict
+end
+_bemcopy(dict::Dict{E, NTuple{2, I}}; boundary_nodes) where {I, E} = copy(dict)
+_bemcopy(dict::Dict{E, Tuple{NTuple{2, I}, I}}; boundary_nodes) where {E, I} = copy(dict)
+
 """
     insert_boundary_node!(boundary_nodes, pos, node)
 
diff --git a/test/data_structures/adjacent.jl b/test/data_structures/adjacent.jl
index b4d77a3bc..1c63ffe88 100644
--- a/test/data_structures/adjacent.jl
+++ b/test/data_structures/adjacent.jl
@@ -4,12 +4,12 @@ using DataStructures
 using StaticArrays
 
 global def_adj = DT.∅
-global default_1 = Dict{NTuple{2, Int}, Int}()
-global default_2 = Dict{NTuple{2, Int32}, Int32}()
-global default_3 = Dict{Vector{Int}, Int}()
-global adj_1 = DT.Adjacent{Int, NTuple{2, Int}}()
-global adj_2 = DT.Adjacent{Int32, NTuple{2, Int32}}()
-global adj_3 = DT.Adjacent{Int, Vector{Int}}()
+global default_1 = Dict{NTuple{2,Int},Int}()
+global default_2 = Dict{NTuple{2,Int32},Int32}()
+global default_3 = Dict{Vector{Int},Int}()
+global adj_1 = DT.Adjacent{Int,NTuple{2,Int}}()
+global adj_2 = DT.Adjacent{Int32,NTuple{2,Int32}}()
+global adj_3 = DT.Adjacent{Int,Vector{Int}}()
 
 @testset "Constructors and getters" begin
     @test adj_1.adjacent == default_1
@@ -28,7 +28,7 @@ global dict_2 = Dict(
     @SVector[1, 2] => 4, @SVector[2, 3] => 10, @SVector[5, 6] => 15,
     @SVector[20, 5] => 72,
 )
-global dict_3 = Dict{NTuple{2, Int32}, Int32}(
+global dict_3 = Dict{NTuple{2,Int32},Int32}(
     (1, 2) => 4, (2, 3) => 10, (5, 6) => 15,
     (20, 5) => 72,
 )
@@ -112,3 +112,12 @@ end
     empty!(adj)
     @test isempty(get_adjacent(adj))
 end
+
+@testset "copy/deepcopy" begin
+    tri = triangulate(rand(2, 50))
+    adj = get_adjacent(tri)
+    adj2 = copy(adj)
+    @test adj == adj2 && !(adj === adj2)
+    adj2 = deepcopy(adj)
+    @test adj == adj2 && !(adj === adj2)
+end
\ No newline at end of file
diff --git a/test/data_structures/adjacent2vertex.jl b/test/data_structures/adjacent2vertex.jl
index 99d31f216..5f04284c7 100644
--- a/test/data_structures/adjacent2vertex.jl
+++ b/test/data_structures/adjacent2vertex.jl
@@ -3,12 +3,12 @@ const DT = DelaunayTriangulation
 using DataStructures
 using StaticArrays
 
-global dict_1 = Dict{Int, Set{NTuple{2, Int}}}()
-global dict_2 = Dict{Int, Vector{NTuple{2, Int}}}()
-global dict_3 = Dict{Int32, Set{SVector{2, Int32}}}()
-global adj2v_1 = DT.Adjacent2Vertex{Int, Set{NTuple{2, Int}}}()
-global adj2v_2 = DT.Adjacent2Vertex{Int, Vector{NTuple{2, Int}}}()
-global adj2v_3 = DT.Adjacent2Vertex{Int32, Set{SVector{2, Int32}}}()
+global dict_1 = Dict{Int,Set{NTuple{2,Int}}}()
+global dict_2 = Dict{Int,Vector{NTuple{2,Int}}}()
+global dict_3 = Dict{Int32,Set{SVector{2,Int32}}}()
+global adj2v_1 = DT.Adjacent2Vertex{Int,Set{NTuple{2,Int}}}()
+global adj2v_2 = DT.Adjacent2Vertex{Int,Vector{NTuple{2,Int}}}()
+global adj2v_3 = DT.Adjacent2Vertex{Int32,Set{SVector{2,Int32}}}()
 
 @testset "Constructors and getters" begin
     @test adj2v_1.adjacent2vertex == dict_1
@@ -32,28 +32,28 @@ global dict_2 = Dict(
     2 => [(5, 7), (10, 14), (2, 3), (5, 9)],
     3 => [(10, 25), (23, 29)],
 )
-global dict_3 = Dict{Int32, Set{SVector{2, Int32}}}(
-    1 => Set{SVector{2, Int32}}(
+global dict_3 = Dict{Int32,Set{SVector{2,Int32}}}(
+    1 => Set{SVector{2,Int32}}(
         (
-            @SVector[1, 2],
-            @SVector[3, 4],
-            @SVector[10, 15],
-            @SVector[2, 5],
-        ),
+        @SVector[1, 2],
+        @SVector[3, 4],
+        @SVector[10, 15],
+        @SVector[2, 5],
     ),
-    2 => Set{SVector{2, Int32}}(
+    ),
+    2 => Set{SVector{2,Int32}}(
         (
-            @SVector[5, 7],
-            @SVector[10, 14],
-            @SVector[2, 3],
-            @SVector[5, 9],
-        ),
+        @SVector[5, 7],
+        @SVector[10, 14],
+        @SVector[2, 3],
+        @SVector[5, 9],
+    ),
     ),
-    3 => Set{SVector{2, Int32}}(
+    3 => Set{SVector{2,Int32}}(
         (
-            @SVector[10, 25],
-            @SVector[23, 29],
-        ),
+        @SVector[10, 25],
+        @SVector[23, 29],
+    ),
     ),
 )
 global adj2v_1 = DT.Adjacent2Vertex(dict_1)
@@ -75,20 +75,20 @@ global adj2v_3 = DT.Adjacent2Vertex(dict_3)
             @test w3 == [(10, 25), (23, 29)]
         elseif adj2v === adj2v_3
             @test w1 ==
-                Set{SVector{2, Int32}}(
+                  Set{SVector{2,Int32}}(
                 (
-                    @SVector[1, 2], @SVector[3, 4], @SVector[10, 15],
-                    @SVector[2, 5],
-                ),
+                @SVector[1, 2], @SVector[3, 4], @SVector[10, 15],
+                @SVector[2, 5],
+            ),
             )
             @test w2 ==
-                Set{SVector{2, Int32}}(
+                  Set{SVector{2,Int32}}(
                 (
-                    @SVector[5, 7], @SVector[10, 14], @SVector[2, 3],
-                    @SVector[5, 9],
-                ),
+                @SVector[5, 7], @SVector[10, 14], @SVector[2, 3],
+                @SVector[5, 9],
+            ),
             )
-            @test w3 == Set{SVector{2, Int32}}((@SVector[10, 25], @SVector[23, 29]))
+            @test w3 == Set{SVector{2,Int32}}((@SVector[10, 25], @SVector[23, 29]))
         end
 
         DT.add_adjacent2vertex!(adj2v, 1, 23, 50)
@@ -103,24 +103,24 @@ global adj2v_3 = DT.Adjacent2Vertex(dict_3)
             @test w3 == [(10, 25), (23, 29), (38, 173)]
         elseif adj2v === adj2v_3
             @test w1 ==
-                Set{SVector{2, Int32}}(
+                  Set{SVector{2,Int32}}(
                 (
-                    @SVector[1, 2], @SVector[3, 4], @SVector[10, 15],
-                    @SVector[2, 5], @SVector[23, 50],
-                ),
+                @SVector[1, 2], @SVector[3, 4], @SVector[10, 15],
+                @SVector[2, 5], @SVector[23, 50],
+            ),
             )
             @test w2 ==
-                Set{SVector{2, Int32}}(
+                  Set{SVector{2,Int32}}(
                 (
-                    @SVector[5, 7], @SVector[10, 14], @SVector[2, 3],
-                    @SVector[5, 9],
-                ),
+                @SVector[5, 7], @SVector[10, 14], @SVector[2, 3],
+                @SVector[5, 9],
+            ),
             )
-            @test w3 == Set{SVector{2, Int32}}(
+            @test w3 == Set{SVector{2,Int32}}(
                 (
-                    @SVector[10, 25], @SVector[23, 29],
-                    @SVector[38, 173],
-                ),
+                @SVector[10, 25], @SVector[23, 29],
+                @SVector[38, 173],
+            ),
             )
         end
 
@@ -137,19 +137,19 @@ global adj2v_3 = DT.Adjacent2Vertex(dict_3)
                 @test w3 == [(10, 25), (23, 29), (38, 173)]
             elseif adj2v === adj2v_3
                 @test w1 ==
-                    Set{SVector{2, Int32}}(
+                      Set{SVector{2,Int32}}(
                     (
-                        @SVector[1, 2], @SVector[10, 15], @SVector[2, 5],
-                        @SVector[23, 50],
-                    ),
+                    @SVector[1, 2], @SVector[10, 15], @SVector[2, 5],
+                    @SVector[23, 50],
+                ),
                 )
                 @test w2 ==
-                    Set{SVector{2, Int32}}((@SVector[5, 7], @SVector[10, 14], @SVector[2, 3]))
-                @test w3 == Set{SVector{2, Int32}}(
+                      Set{SVector{2,Int32}}((@SVector[5, 7], @SVector[10, 14], @SVector[2, 3]))
+                @test w3 == Set{SVector{2,Int32}}(
                     (
-                        @SVector[10, 25], @SVector[23, 29],
-                        @SVector[38, 173],
-                    ),
+                    @SVector[10, 25], @SVector[23, 29],
+                    @SVector[38, 173],
+                ),
                 )
             end
         end
@@ -164,14 +164,14 @@ global adj2v_3 = DT.Adjacent2Vertex(dict_3)
             @test w2 == [(5, 7), (10, 14), (2, 3), (5, 9)]
         elseif adj2v === adj2v_3
             @test w1 ==
-                Set{SVector{2, Int32}}(
+                  Set{SVector{2,Int32}}(
                 (
-                    @SVector[1, 2], @SVector[10, 15], @SVector[2, 5],
-                    @SVector[23, 50],
-                ),
+                @SVector[1, 2], @SVector[10, 15], @SVector[2, 5],
+                @SVector[23, 50],
+            ),
             )
             @test w2 ==
-                Set{SVector{2, Int32}}((@SVector[5, 7], @SVector[10, 14], @SVector[2, 3]))
+                  Set{SVector{2,Int32}}((@SVector[5, 7], @SVector[10, 14], @SVector[2, 3]))
         end
 
         DT.add_triangle!(adj2v, 51, 52, 53)
@@ -218,7 +218,7 @@ global adj2v_3 = DT.Adjacent2Vertex(dict_3)
 end
 
 @testset "Seeing if Adjacent2Vertex is empty and clearing empty sets" begin
-    adj2v = DT.Adjacent2Vertex{Int, Set{NTuple{2, Int}}}()
+    adj2v = DT.Adjacent2Vertex{Int,Set{NTuple{2,Int}}}()
     DT.add_adjacent2vertex!(adj2v, 2, 5, 7)
     DT.add_adjacent2vertex!(adj2v, 2, 7, 13)
     DT.add_adjacent2vertex!(adj2v, 13, 5, 23)
@@ -239,3 +239,18 @@ end
     empty!(adj2v)
     @test isempty(get_adjacent2vertex(adj2v))
 end
+
+@testset "copy/deepcopy" begin
+    tri = triangulate(rand(2, 50))
+    adj2v = get_adjacent2vertex(tri)
+    adj2v2 = copy(adj2v)
+    @test adj2v == adj2v2 && !(adj2v === adj2v2)
+    S = get_adjacent2vertex(adj2v, 1)
+    S2 = get_adjacent2vertex(adj2v2, 1)
+    @test S === S2 # shallow
+    adj2v3 = deepcopy(adj2v)
+    @test adj2v == adj2v3 && !(adj2v === adj2v3)
+    S = get_adjacent2vertex(adj2v, 1)
+    S2 = get_adjacent2vertex(adj2v3, 1)
+    @test !(S === S2) # deep
+end
\ No newline at end of file
diff --git a/test/data_structures/convex_hull.jl b/test/data_structures/convex_hull.jl
index 07afff0a3..4c0d45d38 100644
--- a/test/data_structures/convex_hull.jl
+++ b/test/data_structures/convex_hull.jl
@@ -77,3 +77,11 @@ end
     points = rand(2, 50)
     @test convex_hull(points).vertices == convex_hull(vcat(points, points)).vertices
 end
+
+@testset "copy/deepcopy" begin
+    ch = convex_hull(rand(2, 50))
+    ch2 = copy(ch)
+    @test ch == ch2 && !(ch === ch2)
+    ch2 = deepcopy(ch)
+    @test ch == ch2 && !(ch === ch2)
+end
diff --git a/test/data_structures/curves.jl b/test/data_structures/curves.jl
index 967b06ee6..6de9878c9 100644
--- a/test/data_structures/curves.jl
+++ b/test/data_structures/curves.jl
@@ -1148,7 +1148,7 @@ end
     @test revspl.control_points == _revspl.control_points
     for t in LinRange(0, 1, 1000)
         @test spl(t) ⪧ _revspl(1 - t)
-        @test revspl(t) ⪧ _revspl(t) 
+        @test revspl(t) ⪧ _revspl(t)
     end
     @test revspl.lookup_table ⪧ _revspl.lookup_table
     @test revspl.orientation_markers ≈ _revspl.orientation_markers
@@ -1756,7 +1756,7 @@ end
     for t in LinRange(0, 1, 1500)
         @test spl(t) ⪧ _revspl(1 - t)
         @test revspl(t) ⪧ _revspl(t)
-    end 
+    end
     @test _revspl.control_points == revspl.control_points
     @test _revspl.knots == revspl.knots
     @test _revspl.cache == revspl.cache
@@ -2330,7 +2330,7 @@ end
     @test CatmullRomSpline(ctrl1, _tension=0.2) ≠ CatmullRomSpline(ctrl1)
 
     ## reverse 
-     control_points = [(-9.0, 5.0), (-7.0, 6.0), (-6.0, 4.0), (-3.0, 5.0), (0.0, 0.0), (2.0, -5.0), (2.0, 10.0), (5.0, 12.0), (-9.0, 10.0)]
+    control_points = [(-9.0, 5.0), (-7.0, 6.0), (-6.0, 4.0), (-3.0, 5.0), (0.0, 0.0), (2.0, -5.0), (2.0, 10.0), (5.0, 12.0), (-9.0, 10.0)]
     spl = CatmullRomSpline(control_points)
     revspl = CatmullRomSpline(reverse(control_points))
     _revspl = reverse(spl)
@@ -2349,7 +2349,7 @@ end
         @test DT.differentiate(revspl, t) ⪧ DT.differentiate(_revspl, t)
         @test DT.twice_differentiate(revspl, t) ⪧ DT.twice_differentiate(_revspl, t)
         @test DT.thrice_differentiate(revspl, t) ⪧ DT.thrice_differentiate(_revspl, t)
-    end 
+    end
 end
 
 @testset "angle_between" begin
@@ -2558,3 +2558,97 @@ end
     @test θ56 ≈ 360 - 101.3099324740202
     @test θ61 ≈ 360 - 348.6900675259798
 end
+
+@testset "copy/deepcopy" begin
+    @testset "LineSegment" begin
+        L = LineSegment((0.0, 0.0), (1.0, 1.0))
+        L2 = copy(L)
+        @test L === L
+        L2 = deepcopy(L)
+        @test L === L2
+    end
+
+    @testset "PiecewiseLinear" begin
+        points = [(0.0, 0.0), (10.0, 0.0), (10.0, 10.0)]
+        boundary = [[1, 2, 3]]
+        PL = DT.PiecewiseLinear(points, boundary)
+        PL2 = copy(PL)
+        @test PL == PL2 && !(PL === PL2)
+        @test get_boundary_nodes(PL)[1] === get_boundary_nodes(PL2)[1]
+        PL2 = deepcopy(PL)
+        @test PL == PL2 && !(PL === PL2)
+        @test !(get_boundary_nodes(PL)[1] === get_boundary_nodes(PL2)[1])
+    end
+
+    @testset "BezierCurve" begin
+        control_points = [(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)]
+        spl = BezierCurve(control_points)
+        spl2 = copy(spl)
+        @test spl == spl2 && !(spl === spl2)
+        @test spl.cache == spl2.cache && !(spl.cache === spl2.cache)
+        @test spl.lookup_table == spl2.lookup_table && !(spl.lookup_table === spl2.lookup_table)
+        @test spl.orientation_markers == spl2.orientation_markers && !(spl.orientation_markers === spl2.orientation_markers)
+        spl2 = deepcopy(spl)
+        @test spl == spl2 && !(spl === spl2)
+    end
+
+    @testset "CircularArc" begin
+        p, q, c = (0.0, 0.0), (1.0, 0.0), (0.5, 0.5)
+        cir = CircularArc(p, q, c)
+        cir2 = copy(cir)
+        @test cir === cir2
+        cir2 = deepcopy(cir)
+        @test cir === cir2
+    end
+
+    @testset "EllipticalArc" begin
+        p, q, c, A, B, rot = (0.0, 0.0), (1.0, 0.0), (0.5, 0.5), 1.0, 1.0, 0.0
+        ell = DT.EllipticalArc(p, q, c, A, B, rot)
+        ell2 = copy(ell)
+        @test ell === ell2
+        ell2 = deepcopy(ell)
+        @test ell === ell2
+    end
+
+    @testset "CatmullRomSplineSegment" begin
+        p0, p1, p2, p3, a, t = (0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0), 0.5, 0.5
+        spl = DT.catmull_rom_spline_segment(p0, p1, p2, p3, a, t)
+        spl2 = copy(spl)
+        @test spl === spl2
+        spl2 = deepcopy(spl)
+        @test spl === spl2
+    end
+
+    @testset "CatmullRomSpline" begin
+        control_points = [(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)]
+        spl = CatmullRomSpline(control_points)
+        spl2 = copy(spl)
+        @test spl == spl2 && !(spl === spl2)
+        @test spl.control_points == spl2.control_points && !(spl.control_points === spl2.control_points)
+        @test spl.knots == spl2.knots && !(spl.knots === spl2.knots)
+        @test spl.lookup_table == spl2.lookup_table && !(spl.lookup_table === spl2.lookup_table)
+        @test spl.alpha === spl2.alpha
+        @test spl.tension === spl2.tension
+        @test spl.left === spl2.left
+        @test spl.right === spl2.right
+        @test spl.lengths == spl2.lengths && !(spl.lengths === spl2.lengths)
+        @test spl.segments == spl2.segments && !(spl.segments === spl2.segments)
+        @test spl.orientation_markers == spl2.orientation_markers && !(spl.orientation_markers === spl2.orientation_markers)
+        spl2 = deepcopy(spl)
+        @test spl == spl2 && !(spl === spl2)
+    end
+
+    @testset "BSpline" begin
+        control_points = [(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)]
+        spl = BSpline(control_points)
+        spl2 = copy(spl)
+        @test spl == spl2 && !(spl === spl2)
+        @test spl.control_points == spl2.control_points && !(spl.control_points === spl2.control_points)
+        @test spl.knots == spl2.knots && !(spl.knots === spl2.knots)
+        @test spl.lookup_table == spl2.lookup_table && !(spl.lookup_table === spl2.lookup_table)
+        @test spl.orientation_markers == spl2.orientation_markers && !(spl.orientation_markers === spl2.orientation_markers)
+        @test spl.cache == spl2.cache && !(spl.cache === spl2.cache)
+        spl2 = deepcopy(spl)
+        @test spl == spl2 && !(spl === spl2)
+    end
+end
\ No newline at end of file
diff --git a/test/data_structures/graph.jl b/test/data_structures/graph.jl
index 7991ce168..f243e3e62 100644
--- a/test/data_structures/graph.jl
+++ b/test/data_structures/graph.jl
@@ -2,7 +2,6 @@ using ..DelaunayTriangulation
 const DT = DelaunayTriangulation
 using DataStructures
 
-
 @testset "Constructing empty graphs" begin
     g1 = DT.Graph{Int64}()
     g2 = DT.Graph{Int32}()
@@ -164,4 +163,19 @@ end
     @test isempty(get_neighbours(graph))
     @test isempty(DT.get_edges(graph))
     @test isempty(DT.get_vertices(graph))
+end
+
+@testset "copy/deepcopy" begin
+    tri = triangulate(rand(2, 50))
+    graph = get_graph(tri)
+    graph2 = copy(graph)
+    @test graph == graph2 && !(graph === graph2)
+    S = get_neighbours(tri, 1)
+    S2 = get_neighbours(graph2, 1)
+    @test S === S2
+    graph2 = deepcopy(graph)
+    @test graph == graph2 
+    S = get_neighbours(tri, 1)
+    S2 = get_neighbours(graph2, 1)
+    @test S == S2 && !(S === S2)
 end
\ No newline at end of file
diff --git a/test/data_structures/queue.jl b/test/data_structures/queue.jl
index 3fa25fbc1..6ed807739 100644
--- a/test/data_structures/queue.jl
+++ b/test/data_structures/queue.jl
@@ -32,3 +32,11 @@ push!(q1, 1)
 @test q1 ≠ q2
 push!(q2, 1)
 @test q1 == q2
+
+@testset "copy/deepcopy" begin
+    queue = DT.Queue{Int}([1, 2, 3, 4])
+    queue2 = copy(queue)
+    @test queue == queue2 && !(queue === queue2)
+    queue2 = deepcopy(queue)
+    @test queue == queue2 && !(queue === queue2)
+end
\ No newline at end of file
diff --git a/test/data_structures/representative.jl b/test/data_structures/representative.jl
index aa46992b8..ccc0eb8d3 100644
--- a/test/data_structures/representative.jl
+++ b/test/data_structures/representative.jl
@@ -2,9 +2,8 @@ using ..DelaunayTriangulation
 const DT = DelaunayTriangulation
 using StatsBase
 
-
 @testset "Initialise" begin
-    c = DT.RepresentativeCoordinates{Int, Float64}()
+    c = DT.RepresentativeCoordinates{Int,Float64}()
     @test c == DT.RepresentativeCoordinates(0.0, 0.0, 0)
 end
 
@@ -17,7 +16,7 @@ end
 
 @testset "Conversion" begin
     c = DT.RepresentativeCoordinates(0.39182, 0.5919198, 5)
-    newc = convert(DT.RepresentativeCoordinates{Int32, Float32}, c)
+    newc = convert(DT.RepresentativeCoordinates{Int32,Float32}, c)
     @test newc.x isa Float32
     @test newc.y isa Float32
     @test newc.n isa Int32
@@ -33,7 +32,7 @@ end
 end
 
 @testset "Adding and deleting points to representative coordinates" begin
-    c = DT.RepresentativeCoordinates{Int, Float64}()
+    c = DT.RepresentativeCoordinates{Int,Float64}()
     pts = [15randn(2) for _ in 1:200]
     foreach(p -> DT.add_point!(c, p), pts)
     mx, my = mean(pts)
@@ -79,7 +78,7 @@ end
 @testset "Polylabel, representative points, and distance" begin
     x, y, x1, x2, x3, x4, x5, y1, y2, y3, y4, y5 = complicated_geometry()
     boundary_nodes, points = convert_boundary_points_to_indices(x, y)
-    tri = triangulate(points; boundary_nodes, delete_ghosts = false)
+    tri = triangulate(points; boundary_nodes, delete_ghosts=false)
     for i in (1, 2, 3)
         local x1, y1, x2, y2, x3, y3, x4, y4, x5, y5
         if i == 3
@@ -126,3 +125,11 @@ end
         @test collect(get_point(tri, -11)) ≈ [x5, y5]
     end
 end
+
+@testset "copy/deepcopy" begin
+    c = DT.RepresentativeCoordinates(0.39182, 0.5919198, 5)
+    c2 = copy(c)
+    @test c == c2 && !(c === c2)
+    c2 = deepcopy(c)
+    @test c == c2 && !(c === c2)
+end
\ No newline at end of file
diff --git a/test/data_structures/rtree.jl b/test/data_structures/rtree.jl
index b5c687db9..2747a7412 100644
--- a/test/data_structures/rtree.jl
+++ b/test/data_structures/rtree.jl
@@ -6,6 +6,10 @@ using Random
 const DT = DelaunayTriangulation
 const SI = SpatialIndexing
 
+@testset "copy for BoundingInterval" begin
+    I = DT.BoundingInterval(1.0, 2.0)
+    @test I === copy(I)
+end
 
 @testset "Expanding a BoundingBox" begin
     bbox = DT.BoundingBox(2, 5, 13, 20)
@@ -14,6 +18,11 @@ const SI = SpatialIndexing
     @test bbox2 == DT.BoundingBox(2 - 0.3, 5 + 0.3, 13 - 0.7, 20 + 0.7)
 end
 
+@testset "copy for BoundingBox" begin
+    bbox = DT.BoundingBox(2, 5, 13, 20)
+    @test copy(bbox) === bbox
+end
+
 @testset "BoundingBox for a set of points" begin
     for _ in 1:100
         points = rand(2, 50000)
@@ -54,6 +63,7 @@ end
         for i in axes(points, 2)
             for j in axes(points, 2)
                 bbox = DT.bounding_box(points, i, j)
+                @test copy(bbox) === bbox
                 _bbox = DT.get_bounding_box(bbox)
                 _edge = DT.get_edge(bbox)
                 @test _edge == (i, j)
@@ -78,6 +88,11 @@ end
     leaf1 = DT.Leaf{DT.Branch}(parent, bbox, [dbx1, dbx2])
     leaf2 = DT.Leaf{DT.Branch}(parent, bbox, [dbx1, dbx2])
     @test leaf1 == leaf2
+    @test leaf1 == leaf2
+    @test leaf1 == copy(leaf1) && !(leaf1 === copy(leaf1))
+    @test typeof(copy(leaf1)) === typeof(leaf1)
+    @test leaf1 == deepcopy(leaf1) && !(leaf1 === deepcopy(leaf1))
+    @test typeof(deepcopy(leaf1)) === typeof(leaf1)
     leaf2 = DT.Leaf{DT.Branch}(parent, bbox1, [dbx1, dbx2])
     @test leaf1 ≠ leaf2
     leaf2 = DT.Leaf{DT.Branch}(parent, bbox, [dbx1])
@@ -85,12 +100,21 @@ end
     leaf1 = DT.Leaf{DT.Branch}(DT.Branch(), bbox, [dbx1, dbx2])
     leaf2 = DT.Leaf{DT.Branch}(DT.Branch(), bbox, [dbx1, dbx2])
     @test leaf1 == leaf2
+    @test leaf1 == copy(leaf1) && !(leaf1 === copy(leaf1))
+    @test typeof(copy(leaf1)) === typeof(leaf1)
+    @test leaf1 == deepcopy(leaf1) && !(leaf1 === deepcopy(leaf1))
+    @test typeof(deepcopy(leaf1)) === typeof(leaf1)
 
     leaf1 = DT.Leaf{DT.Branch}(parent, bbox, [dbx1, dbx2])
     leaf2 = DT.Leaf{DT.Branch}(DT.Branch(), bbox, [dbx1, dbx2])
     branch1 = DT.Branch(parent, bbox, [leaf1, leaf2], 1)
     branch2 = DT.Branch(parent, bbox, [leaf1, leaf2], 1)
     @test branch1 == branch2
+    @test leaf1 == leaf2
+    @test branch1 == copy(branch1) && !(branch1 === copy(branch1))
+    @test typeof(copy(branch1)) === typeof(branch1)
+    @test branch1 == deepcopy(branch1) && !(branch1 === deepcopy(branch1))
+    @test typeof(deepcopy(branch1)) === typeof(branch1)
     branch2 = DT.Branch(parent, bbox, [leaf1, leaf2], 2)
     @test branch1 ≠ branch2
     branch2 = DT.Branch(parent, bbox1, [leaf1, leaf2], 1)
@@ -99,6 +123,11 @@ end
     @test branch1 ≠ branch2
     branch1 = DT.Branch(DT.Branch(), bbox, [leaf1, leaf2], 1)
     branch2 = DT.Branch(DT.Branch(), bbox, [leaf1, leaf2], 1)
+    @test leaf1 == leaf2
+    @test branch1 == copy(branch1) && !(branch1 === copy(branch1))
+    @test typeof(copy(branch1)) === typeof(branch1)
+    @test branch1 == deepcopy(branch1) && !(branch1 === deepcopy(branch1))
+    @test typeof(deepcopy(branch1)) === typeof(branch1)
     @test branch1 == branch2
 end
 
@@ -278,3 +307,55 @@ end
     insert!(tree2, 1, 2)
     @test tree1 == tree2
 end
+
+@testset "RTree copy/deepcopy" begin
+    for f in (copy, deepcopy)
+        tree1 = DT.RTree()
+        tree2 = f(tree1)
+        @test tree1 == tree2 && !(tree1 === tree2)
+        @test typeof(tree1) === typeof(tree2)
+        @test DT.get_branch_cache(tree1) == DT.get_branch_cache(tree2) && !(DT.get_branch_cache(tree1) === DT.get_branch_cache(tree2))
+        @test DT.get_twig_cache(tree1) == DT.get_twig_cache(tree2) && !(DT.get_twig_cache(tree1) === DT.get_twig_cache(tree2))
+        @test DT.get_leaf_cache(tree1) == DT.get_leaf_cache(tree2) && !(DT.get_leaf_cache(tree1) === DT.get_leaf_cache(tree2))
+        @test DT.get_fill_factor(tree1) == DT.get_fill_factor(tree2) 
+        @test DT.get_free_cache(tree1) == DT.get_free_cache(tree2) && !(DT.get_free_cache(tree1) === DT.get_free_cache(tree2))
+        @test DT.get_detached_cache(tree1) == DT.get_detached_cache(tree2) && !(DT.get_detached_cache(tree1) === DT.get_detached_cache(tree2))
+        @test DT.get_intersection_cache(tree1) == DT.get_intersection_cache(tree2) && !(DT.get_intersection_cache(tree1) === DT.get_intersection_cache(tree2))
+
+        insert!(tree1, DT.DiametralBoundingBox(DT.BoundingBox((0.0, 0.0)), (1, 2)))
+        insert!(tree1, DT.DiametralBoundingBox(DT.BoundingBox((0.3, 0.0)), (1, 3)))
+        @test tree1 ≠ tree2
+        tree2 = f(tree1)
+        @test tree1 == tree2 && !(tree1 === tree2)
+        @test typeof(tree1) === typeof(tree2)
+        @test DT.get_branch_cache(tree1) == DT.get_branch_cache(tree2) && !(DT.get_branch_cache(tree1) === DT.get_branch_cache(tree2))
+        @test DT.get_twig_cache(tree1) == DT.get_twig_cache(tree2) && !(DT.get_twig_cache(tree1) === DT.get_twig_cache(tree2))
+        @test DT.get_leaf_cache(tree1) == DT.get_leaf_cache(tree2) && !(DT.get_leaf_cache(tree1) === DT.get_leaf_cache(tree2))
+        @test DT.get_fill_factor(tree1) == DT.get_fill_factor(tree2)
+        @test DT.get_free_cache(tree1) == DT.get_free_cache(tree2) && !(DT.get_free_cache(tree1) === DT.get_free_cache(tree2))
+        @test DT.get_detached_cache(tree1) == DT.get_detached_cache(tree2) && !(DT.get_detached_cache(tree1) === DT.get_detached_cache(tree2))
+        @test DT.get_intersection_cache(tree1) == DT.get_intersection_cache(tree2) && !(DT.get_intersection_cache(tree1) === DT.get_intersection_cache(tree2))
+    end
+end
+
+@testset "BoundaryRTree copy/deepcopy" begin
+    for f in (copy, deepcopy)
+        n = 25
+        points = 5 .* randn(2, n)
+        tree1 = DT.BoundaryRTree(points)
+        tree2 = f(tree1)
+        @test tree1 == tree2 && !(tree1 === tree2)
+        @test typeof(tree1) === typeof(tree2)
+        @test tree1.points == tree2.points && !(tree1.points === tree2.points)
+        @test tree1.tree == tree2.tree && !(tree1.tree === tree2.tree)
+
+        insert!(tree1, 1, 2)
+        insert!(tree1, 1, 3)
+        @test tree1 ≠ tree2
+        tree2 = f(tree1)
+        @test tree1 == tree2 && !(tree1 === tree2)
+        @test typeof(tree1) === typeof(tree2)
+        @test tree1.points == tree2.points && !(tree1.points === tree2.points)
+        @test tree1.tree == tree2.tree && !(tree1.tree === tree2.tree)
+    end
+end
\ No newline at end of file
diff --git a/test/data_structures/triangulation.jl b/test/data_structures/triangulation.jl
index a9825f42e..4e99d23f3 100644
--- a/test/data_structures/triangulation.jl
+++ b/test/data_structures/triangulation.jl
@@ -131,7 +131,7 @@ refine!(tri_4; max_area=1.0e-1A, use_circumcenter=true, use_lens=false)
     @test DT.get_ghost_vertex_ranges(tri) == tri.ghost_vertex_ranges == DT.construct_ghost_vertex_ranges(get_boundary_nodes(tri))
     @test DT.get_boundary_edge_map(tri) == tri.boundary_edge_map == DT.construct_boundary_edge_map(get_boundary_nodes(tri))
     @test DT.get_cache(tri) == tri.cache
-    @test DT.get_incircle_cache(tri) === tri.cache.incircle_cache 
+    @test DT.get_incircle_cache(tri) === tri.cache.incircle_cache
     @test DT.get_orient3_cache(tri) === tri.cache.orient3_cache
     @test DT.get_insphere_cache(tri) === tri.cache.insphere_cache
     @inferred DT.get_points(tri)
@@ -1487,3 +1487,120 @@ end
     @test tri1 == tri2
     @test !DT.has_boundary_nodes(tri2)
 end
+
+@testset "copy/deepcopy" begin 
+    tri = triangulate(rand(2, 50))
+    _tri1 = copy(tri)
+    _tri2 = deepcopy(tri)
+    @test typeof(tri) == typeof(_tri1) == typeof(_tri2)
+    for tri2 in (_tri1, _tri2)
+        @test tri == tri2 && !(tri === tri2)
+        for f in fieldnames(typeof(tri))
+            f in (:weights, :boundary_enricher) && continue
+            @test getfield(tri, f) == getfield(tri2, f)
+            f == :boundary_curves && continue
+            @test !(getfield(tri, f) === getfield(tri2, f))
+        end
+    end
+
+    a, b, c, d = [(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)]
+    points = [a, b, c, d]
+    append!(points, tuple.(rand(100), rand(100)))
+    tri = triangulate(points; boundary_nodes=[1, 2, 3, 4, 1], segments=Set([(1, 3)]))
+    _tri1 = copy(tri)
+    _tri2 = deepcopy(tri)
+    @test typeof(tri) == typeof(_tri1) == typeof(_tri2)
+    for tri2 in (_tri1, _tri2)
+        @test tri == tri2 && !(tri === tri2)
+        for f in fieldnames(typeof(tri))
+            f in (:weights, :boundary_enricher, :cache) && continue
+            @test getfield(tri, f) == getfield(tri2, f)
+            f == :boundary_curves && continue
+            @test !(getfield(tri, f) === getfield(tri2, f))
+        end
+        bem = DT.get_boundary_edge_map(tri)
+        _bem = DT.get_boundary_edge_map(tri2)
+        @test bem == _bem && !(bem === _bem)
+        bn = DT.get_boundary_nodes(tri)
+        _bn = DT.get_boundary_nodes(tri2)
+        cbn = first.(values(bem))
+        _cbn = first.(values(_bem))
+        @test all(x -> x === bn, cbn)
+        @test all(x -> x === _bn, _cbn)
+        @test all(x -> !(x === _bn), cbn)
+    end
+
+    curve_XI = [
+        [
+            [1, 2, 3], [DT.EllipticalArc((2.0, 0.0), (-2.0, 0.0), (0.0, 0.0), 2, 1 / 2, 0.0)],
+        ],
+        [
+            [BSpline([(0.0, 0.4), (1.0, 0.2), (0.0, 0.1), (-1.0, 0.2), (0.0, 0.4)])],
+        ],
+        [
+            [4, 5, 6, 7, 4],
+        ],
+        [
+            [BezierCurve(reverse([(-1.0, -3.0), (-1.0, -2.5), (0.0, -2.5), (0.0, -2.0)]))], [CatmullRomSpline(reverse([(0.0, -2.0), (1.0, -3.0), (0.0, -4.0), (-1.0, -3.0)]))],
+        ],
+        [
+            [12, 11, 10, 12],
+        ],
+        [
+            [CircularArc((1.1, -3.0), (1.1, -3.0), (0.0, -3.0), positive=false)],
+        ],
+    ]
+    points_XI = [(-2.0, 0.0), (0.0, 0.0), (2.0, 0.0), (-2.0, -5.0), (2.0, -5.0), (2.0, -1 / 10), (-2.0, -1 / 10), (-1.0, -3.0), (0.0, -4.0), (0.0, -2.3), (-0.5, -3.5), (0.9, -3.0)]
+    fpoints_XI = flatten_boundary_nodes(points_XI, curve_XI)
+    points_XI_extra = copy(points_XI)
+    push!(points_XI_extra, (-1.0, -4.0), (-1.0, -2.0), (0.0, -4.1), (1.0, -4.0), (1.0, -2.0), (1.0, 0.4), (0.0, 0.49), (1.99, -2.0), (-1.99, -2.0), (-1.99, -3.71))
+    tri = triangulate(points_XI_extra; boundary_nodes=curve_XI)
+    _tri1 = copy(tri)
+    _tri2 = deepcopy(tri)
+    @test typeof(tri) == typeof(_tri1) == typeof(_tri2)
+    for tri2 in (_tri1, _tri2)
+        @test tri == tri2 && !(tri === tri2)
+        for f in fieldnames(typeof(tri))
+            f in (:weights, :cache) && continue
+            @test getfield(tri, f) == getfield(tri2, f)
+            @test !(getfield(tri, f) === getfield(tri2, f))
+        end
+        enricher = DT.get_boundary_enricher(tri)
+        _enricher = DT.get_boundary_enricher(tri2)
+        tree = DT.get_spatial_tree(enricher)
+        _tree = DT.get_spatial_tree(_enricher)
+        @test tree == _tree && !(tree === _tree)
+        @test enricher == _enricher && !(enricher === _enricher)
+        @test DT.get_points(enricher) === DT.get_points(tri)
+        @test DT.get_boundary_nodes(enricher) === DT.get_boundary_nodes(tri)
+        @test DT.get_boundary_curves(enricher) === DT.get_boundary_curves(tri)
+        @test DT.get_segments(enricher) === DT.get_interior_segments(tri)
+        @test DT.get_polygon_hierarchy(enricher) === DT.get_polygon_hierarchy(tri)
+        @test DT.get_boundary_edge_map(enricher) === DT.get_boundary_edge_map(tri)
+        @test DT.get_points(_enricher) === DT.get_points(tri2)
+        @test DT.get_boundary_nodes(_enricher) === DT.get_boundary_nodes(tri2)
+        @test DT.get_boundary_curves(_enricher) === DT.get_boundary_curves(tri2)
+        @test DT.get_segments(_enricher) === DT.get_interior_segments(tri2)
+        @test DT.get_polygon_hierarchy(_enricher) === DT.get_polygon_hierarchy(tri2)
+        @test DT.get_boundary_edge_map(_enricher) === DT.get_boundary_edge_map(tri2)
+    end
+
+    tri = triangulate(rand(2, 50); weights = rand(50))
+    _tri1 = copy(tri)
+    _tri2 = deepcopy(tri)
+    @test typeof(tri) == typeof(_tri1) == typeof(_tri2)
+    for tri2 in (_tri1, _tri2)
+        @test tri == tri2 && !(tri === tri2)
+        for f in fieldnames(typeof(tri))
+            f in (:boundary_enricher, :cache) && continue
+            @test getfield(tri, f) == getfield(tri2, f)
+            f in (:boundary_curves,) && continue
+            @test !(getfield(tri, f) === getfield(tri2, f))
+        end
+        @test get_weights(tri) == get_weights(tri2) && !(get_weights(tri) === get_weights(tri2))
+        cache = DT.get_cache(tri)
+        _cache = DT.get_cache(tri2)
+        @test DT.get_weights(DT.get_triangulation(cache)) === DT.get_weights(DT.get_triangulation_2(cache)) === DT.get_weights(tri)
+        @test DT.get_weights(DT.get_triangulation(_cache)) === DT.get_weights(DT.get_triangulation_2(_cache)) === DT.get_weights(tri2)
+    end
+end
\ No newline at end of file
diff --git a/test/data_structures/triangulation_cache.jl b/test/data_structures/triangulation_cache.jl
index 0fb721d04..7a569ddf8 100644
--- a/test/data_structures/triangulation_cache.jl
+++ b/test/data_structures/triangulation_cache.jl
@@ -20,7 +20,7 @@ cache = DT.TriangulationCache()
 @test isnothing(DT.get_surrounding_polygon(cache))
 @test typeof(cache) == DT.EmptyTriangulationCache
 
-tri = triangulate(rand(2, 50); weights = rand(50))
+tri = triangulate(rand(2, 50); weights=rand(50))
 cache = DT.get_cache(tri)
 @test get_weights(DT.get_triangulation(DT.get_cache(tri))) === get_weights(tri)
 DT.unconstrained_triangulation!(DT.get_triangulation(cache))
@@ -88,4 +88,21 @@ empty!(cache)
 @test isempty(surrounding_polygon)
 @test length.(DT.get_incircle_cache(cache)) == length.(AdaptivePredicates.incircleadapt_cache(Float64))
 @test length.(DT.get_orient3_cache(cache)) == length.(AdaptivePredicates.orient3adapt_cache(Float64))
-@test length.(DT.get_insphere_cache(cache)) == length.(AdaptivePredicates.insphereexact_cache(Float64))
\ No newline at end of file
+@test length.(DT.get_insphere_cache(cache)) == length.(AdaptivePredicates.insphereexact_cache(Float64))
+
+@testset "copy/deepcopy" begin
+    cache = DT.TriangulationCache()
+    @test copy(cache) === cache
+
+    tri = triangulate(rand(2, 50))
+    cache = DT.get_cache(tri)
+    cache2 = copy(cache)
+    @inferred copy(cache2)
+    @test typeof(cache2) == typeof(cache) && !(cache2 === cache)
+
+    tri = triangulate(rand(2, 50); weights=rand(50))
+    cache = DT.get_cache(tri)
+    cache2 = copy(cache)
+    @inferred copy(cache2)
+    @test typeof(cache2) == typeof(cache) && !(cache2 === cache)
+end
\ No newline at end of file
diff --git a/test/interfaces/boundary_nodes.jl b/test/interfaces/boundary_nodes.jl
index 980ff6617..cd1f61d43 100644
--- a/test/interfaces/boundary_nodes.jl
+++ b/test/interfaces/boundary_nodes.jl
@@ -155,6 +155,13 @@ end
         @test get_boundary_nodes(S, k) == ij[1]
         @test get_boundary_nodes(S, k + 1) == ij[2]
     end
+    cbn = copy(bn)
+    _bn_map = DT._bemcopy(bn_map; boundary_nodes=cbn)
+    @test bn_map == _bn_map 
+    _bn = first.(values(_bn_map))
+    @test all(x -> x === cbn, _bn)
+    @test all(x -> !(x === bn), _bn)
+    @test all(==(bn), _bn)
     bn = [[1, 2, 3, 4], [4, 5, 6, 7, 8], [8, 9, 10, 1]]
     bn_map = DT.construct_boundary_edge_map(bn)
     for (ij, (index, k)) in bn_map
@@ -162,6 +169,8 @@ end
         @test get_boundary_nodes(S, k) == ij[1]
         @test get_boundary_nodes(S, k + 1) == ij[2]
     end
+    _bn_map = DT._bemcopy(bn_map; boundary_nodes=bn)
+    @test _bn_map == bn_map && !(bn_map === _bn_map)
     bn = [
         [[1, 2, 3, 4, 5], [5, 6, 7], [7, 8], [8, 9, 10, 1]],
         [[13, 14, 15, 16, 17], [17, 18, 19, 20], [20, 13]],
@@ -172,6 +181,8 @@ end
         @test get_boundary_nodes(S, k) == ij[1]
         @test get_boundary_nodes(S, k + 1) == ij[2]
     end
+    _bn_map = DT._bemcopy(bn_map; boundary_nodes=bn)
+    @test _bn_map == bn_map && !(bn_map === _bn_map)
     bn = Int[]
     bn_map = DT.construct_boundary_edge_map(bn)
     @test bn_map == Dict{Tuple{Int32, Int32}, Tuple{Vector{Int}, Int}}()
diff --git a/test/refinement/curve_bounded.jl b/test/refinement/curve_bounded.jl
index 30609ee25..d1879b929 100644
--- a/test/refinement/curve_bounded.jl
+++ b/test/refinement/curve_bounded.jl
@@ -496,6 +496,25 @@ end
     @test DT.reorient_edge(enricher, 23, 14) == (14, 23)
 end
 
+@testset "copy/deepcopy" begin 
+    ppoints_XI = deepcopy(points_XI)
+    ccurve_XI = deepcopy(curve_XI)
+    tri = triangulate(ppoints_XI; boundary_nodes = ccurve_XI)
+    enricher = DT.get_boundary_enricher(tri)
+    enricher2 = copy(enricher)
+    @test enricher == enricher2 && !(enricher === enricher2)
+    @test get_points(enricher) === get_points(get_boundary_curves(enricher)[1])
+    @test get_points(enricher) === get_points(get_boundary_curves(enricher)[4])
+    @test get_points(enricher) === get_points(get_boundary_curves(enricher)[7])
+    @test get_points(enricher) === DT.get_spatial_tree(enricher).points
+    enricher2 = deepcopy(enricher)
+    @test enricher == enricher2 && !(enricher === enricher2)
+    @test get_points(enricher) === get_points(get_boundary_curves(enricher)[1])
+    @test get_points(enricher) === get_points(get_boundary_curves(enricher)[4])
+    @test get_points(enricher) === get_points(get_boundary_curves(enricher)[7])
+    @test get_points(enricher) === DT.get_spatial_tree(enricher).points
+end
+
 @testset "coarse_discretisation" begin
     # I 
     points = deepcopy(points_I)
@@ -1018,6 +1037,9 @@ end
     member = SACM(1, 7)
     nmember = DT.replace_next_edge(member, 12)
     @test DT.get_parent_curve(nmember) == 1 && DT.get_next_edge(nmember) == 12
+    membercopy = copy(member)
+    @test member === membercopy
+    @test member === deepcopy(member)
 
     points, boundary_nodes = deepcopy(points_I), deepcopy(curve_I)
     enricher = DT.BoundaryEnricher(points, boundary_nodes)
@@ -1058,6 +1080,11 @@ end
         1 => [SAC(1, [SACM(3, 11), SACM(1, 2)])],
     )
     @test DT.get_small_angle_complexes(enricher) == _complexes
+    cp1 = DT.get_small_angle_complexes(enricher)
+    cp2 = copy(cp1)
+    @test cp1 == cp2 && !(cp1 === cp2)
+    cp2 = deepcopy(cp1)
+    @test cp1 == cp2 && !(cp1 === cp2)
 
     A, B, C, D, E, F, G, H, I, J, K = (0.0, 0.0), (0.2, 1.4), (0.6, 1.2),
     (1.2, 0.2), (1.2, -0.2), (-1.4, -0.2),
diff --git a/test/triangulation/weighted.jl b/test/triangulation/weighted.jl
index e743f9765..ccf6baa2f 100644
--- a/test/triangulation/weighted.jl
+++ b/test/triangulation/weighted.jl
@@ -12,6 +12,7 @@ using ..DelaunayTriangulation: add_weight!, get_weight, get_weights
     zw = DT.ZeroWeight{Float32}()
     @inferred DT.ZeroWeight{Float32}()
     @test get_weight(zw, 1) ⊢ 0.0f0
+    @test copy(zw) === zw
 end
 
 @testset "get_weight" begin
diff --git a/test/voronoi/voronoi.jl b/test/voronoi/voronoi.jl
index c9ff9367d..b2cd1f7a3 100644
--- a/test/voronoi/voronoi.jl
+++ b/test/voronoi/voronoi.jl
@@ -1713,4 +1713,64 @@ end
         c = DT.get_centroid(vorn, index)
         @test c ⪧ get_generator(vorn, index) rtol=1e-2 atol=1e-4
     end
+end
+
+@testset "==" begin
+    tri = triangulate(rand(2, 100))
+    vorn = voronoi(tri, clip=true)
+    vorn2 = voronoi(tri, clip=true)
+    @test vorn == vorn
+    @test vorn == vorn2
+    g1 = vorn.generators[1]
+    vorn.generators[1] = vorn.generators[2] 
+    @test vorn != vorn2
+    vorn.generators[1] = g1
+    @test vorn == vorn2
+    p1 = vorn.polygon_points[1]
+    vorn.polygon_points[1] = vorn.polygon_points[2]
+    @test vorn != vorn2
+    vorn.polygon_points[1] = p1
+    @test vorn == vorn2
+    poly1 = vorn.polygons[1]
+    vorn.polygons[1] = vorn.polygons[2]
+    @test vorn != vorn2
+    vorn.polygons[1] = poly1
+    @test vorn == vorn2
+    T1 = vorn.circumcenter_to_triangle[1]
+    vorn.circumcenter_to_triangle[1] = vorn.circumcenter_to_triangle[2]
+    @test vorn != vorn2
+    vorn.circumcenter_to_triangle[1] = T1
+    @test vorn == vorn2
+    T2 = vorn.circumcenter_to_triangle[2]
+    p1 = vorn.triangle_to_circumcenter[T1]
+    vorn.triangle_to_circumcenter[T1] = vorn.triangle_to_circumcenter[T2]
+    @test vorn != vorn2
+    vorn.triangle_to_circumcenter[T1] = p1
+    @test vorn == vorn2
+    push!(vorn.unbounded_polygons, 5)
+    @test vorn != vorn2
+    delete!(vorn.unbounded_polygons, 5)
+    (u, v), w = first(vorn.adjacent.adjacent)
+    vorn.adjacent.adjacent[(u, v)] = w + 1
+    @test vorn != vorn2
+    vorn.adjacent.adjacent[(u, v)] = w
+    @test vorn == vorn2
+    push!(vorn.boundary_polygons, 5000)
+    @test vorn != vorn2
+    delete!(vorn.boundary_polygons, 5000)
+    @test vorn == vorn2
+end
+
+@testset "copy/deepcopy" begin 
+    vorn = voronoi(triangulate(rand(2, 100)), clip = true, smooth = true)
+    vorn1 = deepcopy(vorn)
+    vorn2 = copy(vorn)
+    @test typeof(vorn1) == typeof(vorn) == typeof(vorn2)
+    for _vorn in (vorn1, vorn2)
+        @test vorn == _vorn && !(vorn === _vorn)
+        for f in fieldnames(typeof(vorn))
+            @test getfield(vorn, f) == getfield(_vorn, f)
+            @test !(getfield(vorn, f) === getfield(_vorn, f))
+        end
+    end
 end
\ No newline at end of file