diff --git a/src/DiscreteExteriorCalculus.jl b/src/DiscreteExteriorCalculus.jl
index 365d8f84..c05628f4 100644
--- a/src/DiscreteExteriorCalculus.jl
+++ b/src/DiscreteExteriorCalculus.jl
@@ -20,7 +20,7 @@ export DualSimplex, DualV, DualE, DualTri, DualChain, DualForm,
   ⋆, hodge_star, inv_hodge_star, δ, codifferential, ∇², laplace_beltrami, Δ, laplace_de_rham,
   ♭, flat, ♭_mat, ♯, ♯_mat, sharp, ∧, wedge_product, interior_product, interior_product_flat,
   ℒ, lie_derivative, lie_derivative_flat,
-  vertex_center, edge_center, triangle_center, dual_triangle_vertices,
+  vertex_center, edge_center, triangle_center, dual_triangle_vertices, dual_edge_vertices,
   dual_point, dual_volume, subdivide_duals!, DiagonalHodge, GeometricHodge,
   subdivide, PPSharp, AltPPSharp, DesbrunSharp, LLSDDSharp, de_sign,
   ♭♯, ♭♯_mat, flat_sharp, flat_sharp_mat
@@ -30,7 +30,7 @@ import Base: *
 import LinearAlgebra: mul!
 using LinearAlgebra: Diagonal, dot, norm, cross, pinv, qr, ColumnNorm
 using SparseArrays
-using StaticArrays: @SVector, SVector, SMatrix, MMatrix
+using StaticArrays: @SVector, SVector, SMatrix, MVector, MMatrix
 using GeometryBasics: Point2, Point3
 
 const Point2D = SVector{2,Float64}
@@ -174,6 +174,16 @@ elementary_duals(::Type{Val{0}}, s::AbstractDeltaDualComplex1D, v::Int) =
 elementary_duals(::Type{Val{1}}, s::AbstractDeltaDualComplex1D, e::Int) =
   SVector(edge_center(s,e))
 
+""" Boundary dual vertices of a dual edge.
+
+This accessor assumes that the simplicial identities for the dual hold.
+"""
+function dual_edge_vertices(s::HasDeltaSet1D, t...)
+    SVector(s[t..., :D_∂v0],
+            s[t..., :D_∂v1])
+end
+
+
 """ Boundary dual vertices of a dual triangle.
 
 This accessor assumes that the simplicial identities for the dual hold.
@@ -233,11 +243,37 @@ negatenz((I, V)) = (I, negate.(V))
 
 """ Construct 1D dual complex from 1D delta set.
 """
-function (::Type{S})(t::AbstractDeltaSet1D) where S <: AbstractDeltaDualComplex1D
-  s = S()
-  copy_parts!(s, t)
-  make_dual_simplices_1d!(s)
-  return s
+function (::Type{S})(s::AbstractDeltaSet1D) where S <: AbstractDeltaDualComplex1D
+  t = S()
+  copy_primal_1D!(t, s) # TODO: Revert to copy_parts! when performance is improved
+  make_dual_simplices_1d!(t)
+  return t
+end
+
+function copy_primal_1D!(t::HasDeltaSet1D, s::HasDeltaSet1D)
+
+  @assert nv(t) == 0
+  @assert ne(t) == 0
+
+  v_range = add_parts!(t, :V, nv(s))
+  e_range = add_parts!(t, :E, ne(s))
+
+  if has_subpart(s, :point)
+    @inbounds for v in v_range
+      t[v, :point] = s[v, :point]
+    end
+  end
+
+  @inbounds for e in e_range
+    t[e, :∂v0] = s[e, :∂v0]
+    t[e, :∂v1] = s[e, :∂v1]
+  end
+
+  if has_subpart(s, :edge_orientation)
+    @inbounds for e in e_range
+      t[e, :edge_orientation] = s[e, :edge_orientation]
+    end
+  end
 end
 
 make_dual_simplices_1d!(s::AbstractDeltaDualComplex1D) = make_dual_simplices_1d!(s, E)
@@ -364,27 +400,48 @@ Supports different methods of subdivision through the choice of geometric
 center, as defined by [`geometric_center`](@ref). In particular, barycentric
 subdivision and circumcentric subdivision are supported.
 """
-function subdivide_duals!(s::EmbeddedDeltaDualComplex1D, args...)
-  subdivide_duals_1d!(s, args...)
-  precompute_volumes_1d!(s)
+function subdivide_duals!(sd::EmbeddedDeltaDualComplex1D{_o, _l, point_type} where {_o, _l}, alg) where point_type
+  subdivide_duals_1d!(sd, point_type, alg)
+  precompute_volumes_1d!(sd, point_type)
 end
 
-function subdivide_duals_1d!(s::HasDeltaSet1D, alg)
-  for v in vertices(s)
-    s[vertex_center(s,v), :dual_point] = point(s, v)
+# TODO: Replace the individual accesses with vector accesses
+function subdivide_duals_1d!(sd::HasDeltaSet1D, ::Type{point_type}, alg) where point_type
+
+  point_arr = MVector{2, point_type}(undef)
+
+  @inbounds for v in vertices(sd)
+    sd[v, :dual_point] = sd[v, :point]
   end
-  for e in edges(s)
-    s[edge_center(s,e), :dual_point] = geometric_center(
-      point(s, edge_vertices(s, e)), alg)
+
+  @inbounds for e in edges(sd)
+    p1, p2 = edge_vertices(sd, e)
+    point_arr[1] = sd[p1, :point]
+    point_arr[2] = sd[p2, :point]
+
+    sd[sd[e, :edge_center], :dual_point] = geometric_center(point_arr, alg)
   end
 end
 
-function precompute_volumes_1d!(s::HasDeltaSet1D)
-  for e in edges(s)
-    s[e, :length] = volume(1,s,e,CayleyMengerDet())
+# TODO: Replace the individual accesses with vector accesses
+function precompute_volumes_1d!(sd::HasDeltaSet1D, ::Type{point_type}) where point_type
+
+  point_arr = MVector{2, point_type}(undef)
+
+  @inbounds for e in edges(sd)
+    p1, p2 = edge_vertices(sd, e)
+    point_arr[1] = sd[p1, :point]
+    point_arr[2] = sd[p2, :point]
+
+    sd[e, :length] = volume(point_arr)
   end
-  for e in parts(s, :DualE)
-    s[e, :dual_length] = dual_volume(1,s,e,CayleyMengerDet())
+
+  @inbounds for e in parts(sd, :DualE)
+    p1, p2 = dual_edge_vertices(sd, e)
+    point_arr[1] = sd[p1, :dual_point]
+    point_arr[2] = sd[p2, :dual_point]
+
+    sd[e, :dual_length] = volume(point_arr)
   end
 end
 
@@ -481,11 +538,31 @@ dual_derivative_nz(::Type{Val{1}}, s::AbstractDeltaDualComplex2D, x::Int) =
 
 """ Construct 2D dual complex from 2D delta set.
 """
-function (::Type{S})(t::AbstractDeltaSet2D) where S <: AbstractDeltaDualComplex2D
-  s = S()
-  copy_parts!(s, t)
-  make_dual_simplices_2d!(s)
-  return s
+function (::Type{S})(s::AbstractDeltaSet2D) where S <: AbstractDeltaDualComplex2D
+  t = S()
+  copy_primal_2D!(t, s) # TODO: Revert to copy_parts! when performance is improved
+  make_dual_simplices_2d!(t)
+  return t
+end
+
+function copy_primal_2D!(t::HasDeltaSet2D, s::HasDeltaSet2D)
+
+  @assert ntriangles(t) == 0
+
+  copy_primal_1D!(t, s)
+  tri_range = add_parts!(t, :Tri, ntriangles(s))
+
+  @inbounds for tri in tri_range
+    t[tri, :∂e0] = s[tri, :∂e0]
+    t[tri, :∂e1] = s[tri, :∂e1]
+    t[tri, :∂e2] = s[tri, :∂e2]
+  end
+
+  if has_subpart(s, :tri_orientation)
+    @inbounds for tri in tri_range
+      t[tri, :tri_orientation] = s[tri, :tri_orientation]
+    end
+  end
 end
 
 make_dual_simplices_1d!(s::AbstractDeltaDualComplex2D) = make_dual_simplices_1d!(s, Tri)
@@ -720,7 +797,7 @@ function get_orthogonal_vector(s::AbstractDeltaDualComplex2D, v::Int, e::Int)
   e2_vec - dot(e2_vec, e_vec)*e_vec
 end
 
-function ♯_assign!(♯_mat::AbstractSparseMatrix, s::AbstractDeltaDualComplex2D, 
+function ♯_assign!(♯_mat::AbstractSparseMatrix, s::AbstractDeltaDualComplex2D,
   v₀::Int, _::Int, t::Int, i::Int, tri_edges::SVector{3, Int}, tri_center::Int,
   out_vec, DS::DiscreteSharp)
   for e in deleteat(tri_edges, i)
@@ -733,7 +810,7 @@ function ♯_assign!(♯_mat::AbstractSparseMatrix, s::AbstractDeltaDualComplex2
   end
 end
 
-function ♯_assign!(♯_mat::AbstractSparseMatrix, s::AbstractDeltaDualComplex2D, 
+function ♯_assign!(♯_mat::AbstractSparseMatrix, s::AbstractDeltaDualComplex2D,
   _::Int, e₀::Int, t::Int, _::Int, _::SVector{3, Int}, tri_center::Int,
   out_vec, DS::DesbrunSharp)
   for v in edge_vertices(s, e₀)
@@ -832,26 +909,60 @@ function ∧(::Type{Tuple{1,1}}, s::HasDeltaSet2D, α, β, x::Int)
                       α[e1] * β[e0] - α[e0] * β[e1])) / 2
 end
 
-function subdivide_duals!(s::EmbeddedDeltaDualComplex2D, args...)
-  subdivide_duals_2d!(s, args...)
-  precompute_volumes_2d!(s)
+function subdivide_duals!(sd::EmbeddedDeltaDualComplex2D{_o, _l, point_type} where {_o, _l}, alg) where point_type
+  subdivide_duals_2d!(sd, point_type, alg)
+  precompute_volumes_2d!(sd, point_type)
 end
 
-function subdivide_duals_2d!(s::HasDeltaSet2D, alg)
-  subdivide_duals_1d!(s, alg)
-  for t in triangles(s)
-    s[triangle_center(s,t), :dual_point] = geometric_center(
-      point(s, triangle_vertices(s, t)), alg)
+# TODO: Replace the individual accesses with vector accesses
+function subdivide_duals_2d!(sd::HasDeltaSet2D, ::Type{point_type}, alg) where point_type
+  subdivide_duals_1d!(sd, point_type, alg)
+
+  point_arr = MVector{3, point_type}(undef)
+
+  @inbounds for t in triangles(sd)
+    p1, p2, p3 = triangle_vertices(sd, t)
+    point_arr[1] = sd[p1, :point]
+    point_arr[2] = sd[p2, :point]
+    point_arr[3] = sd[p3, :point]
+
+    sd[sd[t, :tri_center], :dual_point] = geometric_center(point_arr, alg)
   end
 end
 
-function precompute_volumes_2d!(s::HasDeltaSet2D)
-  precompute_volumes_1d!(s)
-  for t in triangles(s)
-    s[t, :area] = volume(2,s,t,CayleyMengerDet())
+function precompute_volumes_2d!(sd::HasDeltaSet2D, p::Type{point_type}) where point_type
+  precompute_volumes_1d!(sd, point_type)
+  set_volumes_2d!(Val{2}, sd, p)
+  set_dual_volumes_2d!(Val{2}, sd, p)
+end
+
+# TODO: Replace the individual accesses with vector accesses
+function set_volumes_2d!(::Type{Val{2}}, sd::HasDeltaSet2D, ::Type{point_type}) where point_type
+
+  point_arr = MVector{3, point_type}(undef)
+
+  @inbounds for t in triangles(sd)
+    p1, p2, p3 = triangle_vertices(sd, t)
+    point_arr[1] = sd[p1, :point]
+    point_arr[2] = sd[p2, :point]
+    point_arr[3] = sd[p3, :point]
+
+    sd[t, :area] = volume(point_arr)
   end
-  for t in parts(s, :DualTri)
-    s[t, :dual_area] = dual_volume(2,s,t,CayleyMengerDet())
+end
+
+# TODO: Replace the individual accesses with vector accesses
+function set_dual_volumes_2d!(::Type{Val{2}}, sd::HasDeltaSet2D, ::Type{point_type}) where point_type
+
+  point_arr = MVector{3, point_type}(undef)
+
+  @inbounds for t in parts(sd, :DualTri)
+    p1, p2, p3 = dual_triangle_vertices(sd, t)
+    point_arr[1] = sd[p1, :dual_point]
+    point_arr[2] = sd[p2, :dual_point]
+    point_arr[3] = sd[p3, :dual_point]
+
+    sd[t, :dual_area] = volume(point_arr)
   end
 end
 
diff --git a/test/Benchmarks.jl b/test/Benchmarks.jl
index 11b1b64e..97342456 100644
--- a/test/Benchmarks.jl
+++ b/test/Benchmarks.jl
@@ -10,8 +10,8 @@ using GeometryBasics: Point2, Point3
 Point2D = Point2{Float64}
 Point3D = Point3{Float64}
 
-begin 
-    # TODO: Support larger meshes like Ico7
+@info "Beginning DEC Operator Benchmarks"
+begin
     mesh_size = 5
     float_type::DataType = Float64
     primal_earth = loadmesh(Icosphere(mesh_size))
@@ -20,7 +20,7 @@ begin
     subdivide_duals!(earth, Barycenter());
 end
 
-begin 
+begin
     println("Mesh: " * "Blender Ico Sphere, $(mesh_size) Subdivisions")
     println("")
 
@@ -35,70 +35,58 @@ begin
     println("----------------------------------------------------------------")
 end
 
-suite = BenchmarkGroup()
+dec_op_suite = BenchmarkGroup()
 
 begin
-    suite["Exterior Derivative"] = BenchmarkGroup()
-
-    suite["Exterior Derivative"]["New Form-0"] = @benchmarkable dec_differential(0, $earth)
-    suite["Exterior Derivative"]["Old Form-0"] = @benchmarkable d(0, $earth)
+    dec_op_suite["Exterior Derivative"]["New Form-0"] = @benchmarkable dec_differential(0, $earth)
+    dec_op_suite["Exterior Derivative"]["Old Form-0"] = @benchmarkable d(0, $earth)
 
-    suite["Exterior Derivative"]["New Form-1"] = @benchmarkable dec_differential(1, $earth)
-    suite["Exterior Derivative"]["Old Form-1"] = @benchmarkable d(1, $earth)
+    dec_op_suite["Exterior Derivative"]["New Form-1"] = @benchmarkable dec_differential(1, $earth)
+    dec_op_suite["Exterior Derivative"]["Old Form-1"] = @benchmarkable d(1, $earth)
 end
 
 begin
-    suite["Boundary"] = BenchmarkGroup()
+    dec_op_suite["Boundary"]["New Form-1"] = @benchmarkable dec_boundary(1, $earth)
+    dec_op_suite["Boundary"]["Old Form-1"] = @benchmarkable ∂(1, $earth)
 
-    suite["Boundary"]["New Form-1"] = @benchmarkable dec_boundary(1, $earth)
-    suite["Boundary"]["Old Form-1"] = @benchmarkable ∂(1, $earth)
-
-    suite["Boundary"]["New Form-2"] = @benchmarkable dec_boundary(2, $earth)
-    suite["Boundary"]["Old Form-2"] = @benchmarkable ∂(2, $earth)
+    dec_op_suite["Boundary"]["New Form-2"] = @benchmarkable dec_boundary(2, $earth)
+    dec_op_suite["Boundary"]["Old Form-2"] = @benchmarkable ∂(2, $earth)
 end
 
 begin
-    suite["Dual Derivative"] = BenchmarkGroup()
-
-    suite["Dual Derivative"]["New Dual-Form-0"] = @benchmarkable dec_dual_derivative(0, $earth)
-    suite["Dual Derivative"]["Old Dual-Form-0"] = @benchmarkable dual_derivative(0, $earth)
+    dec_op_suite["Dual Derivative"]["New Dual-Form-0"] = @benchmarkable dec_dual_derivative(0, $earth)
+    dec_op_suite["Dual Derivative"]["Old Dual-Form-0"] = @benchmarkable dual_derivative(0, $earth)
 
-    suite["Dual Derivative"]["New Dual-Form-1"] = @benchmarkable dec_dual_derivative(1, $earth)
-    suite["Dual Derivative"]["Old Dual-Form-1"] = @benchmarkable dual_derivative(1, $earth)
+    dec_op_suite["Dual Derivative"]["New Dual-Form-1"] = @benchmarkable dec_dual_derivative(1, $earth)
+    dec_op_suite["Dual Derivative"]["Old Dual-Form-1"] = @benchmarkable dual_derivative(1, $earth)
 end
 
 begin
-    suite["Diagonal Hodge"] = BenchmarkGroup()
-
-    suite["Diagonal Hodge"]["New Form-0"] = @benchmarkable dec_hodge_star(0, $earth, DiagonalHodge())
-    suite["Diagonal Hodge"]["Old Form-0"] = @benchmarkable hodge_star(0, $earth, DiagonalHodge())
+    dec_op_suite["Diagonal Hodge"]["New Form-0"] = @benchmarkable dec_hodge_star(0, $earth, DiagonalHodge())
+    dec_op_suite["Diagonal Hodge"]["Old Form-0"] = @benchmarkable hodge_star(0, $earth, DiagonalHodge())
 
-    suite["Diagonal Hodge"]["New Form-1"] = @benchmarkable dec_hodge_star(1, $earth, DiagonalHodge())
-    suite["Diagonal Hodge"]["Old Form-1"] = @benchmarkable hodge_star(1, $earth, DiagonalHodge())
+    dec_op_suite["Diagonal Hodge"]["New Form-1"] = @benchmarkable dec_hodge_star(1, $earth, DiagonalHodge())
+    dec_op_suite["Diagonal Hodge"]["Old Form-1"] = @benchmarkable hodge_star(1, $earth, DiagonalHodge())
 
-    suite["Diagonal Hodge"]["New Form-2"] = @benchmarkable dec_hodge_star(2, $earth, DiagonalHodge())
-    suite["Diagonal Hodge"]["Old Form-2"] = @benchmarkable hodge_star(2, $earth, DiagonalHodge())
+    dec_op_suite["Diagonal Hodge"]["New Form-2"] = @benchmarkable dec_hodge_star(2, $earth, DiagonalHodge())
+    dec_op_suite["Diagonal Hodge"]["Old Form-2"] = @benchmarkable hodge_star(2, $earth, DiagonalHodge())
 
 end
 
 begin
-    suite["Geometric Hodge"] = BenchmarkGroup()
-
-    suite["Geometric Hodge"]["New Form-1"] = @benchmarkable dec_hodge_star(1, $earth, GeometricHodge())
-    suite["Geometric Hodge"]["Old Form-1"] = @benchmarkable hodge_star(1, $earth, GeometricHodge())
+    dec_op_suite["Geometric Hodge"]["New Form-1"] = @benchmarkable dec_hodge_star(1, $earth, GeometricHodge())
+    dec_op_suite["Geometric Hodge"]["Old Form-1"] = @benchmarkable hodge_star(1, $earth, GeometricHodge())
 end
 
 begin
-    suite["Inverse Diagonal Hodge"] = BenchmarkGroup()
-
-    suite["Inverse Diagonal Hodge"]["New Form-0"] = @benchmarkable dec_inv_hodge_star(0, $earth, DiagonalHodge())
-    suite["Inverse Diagonal Hodge"]["Old Form-0"] = @benchmarkable inv_hodge_star(0, $earth, DiagonalHodge())
+    dec_op_suite["Inverse Diagonal Hodge"]["New Form-0"] = @benchmarkable dec_inv_hodge_star(0, $earth, DiagonalHodge())
+    dec_op_suite["Inverse Diagonal Hodge"]["Old Form-0"] = @benchmarkable inv_hodge_star(0, $earth, DiagonalHodge())
 
-    suite["Inverse Diagonal Hodge"]["New Form-1"] = @benchmarkable dec_inv_hodge_star(1, $earth, DiagonalHodge())
-    suite["Inverse Diagonal Hodge"]["Old Form-1"] = @benchmarkable inv_hodge_star(1, $earth, DiagonalHodge())
+    dec_op_suite["Inverse Diagonal Hodge"]["New Form-1"] = @benchmarkable dec_inv_hodge_star(1, $earth, DiagonalHodge())
+    dec_op_suite["Inverse Diagonal Hodge"]["Old Form-1"] = @benchmarkable inv_hodge_star(1, $earth, DiagonalHodge())
 
-    suite["Inverse Diagonal Hodge"]["New Form-2"] = @benchmarkable dec_inv_hodge_star(2, $earth, DiagonalHodge())
-    suite["Inverse Diagonal Hodge"]["Old Form-2"] = @benchmarkable inv_hodge_star(2, $earth, DiagonalHodge())
+    dec_op_suite["Inverse Diagonal Hodge"]["New Form-2"] = @benchmarkable dec_inv_hodge_star(2, $earth, DiagonalHodge())
+    dec_op_suite["Inverse Diagonal Hodge"]["Old Form-2"] = @benchmarkable inv_hodge_star(2, $earth, DiagonalHodge())
 end
 
 begin
@@ -107,16 +95,16 @@ begin
     E_1, E_2 = rand(float_type, ne(earth)), rand(float_type, ne(earth))
     T_2 = rand(float_type, ntriangles(earth))
 
-    suite["Wedge Product"] = BenchmarkGroup()
+    dec_op_suite["Wedge Product"] = BenchmarkGroup()
 
-    suite["Wedge Product"]["New Form-0, Form-1"] = @benchmarkable dec_wedge_product(Tuple{0,1}, $earth)($V_1, $E_1)
-    suite["Wedge Product"]["Old Form-0, Form-1"] = @benchmarkable wedge_product(Tuple{0,1}, $earth, $V_1, $E_1)
+    dec_op_suite["Wedge Product"]["New Form-0, Form-1"] = @benchmarkable dec_wedge_product(Tuple{0,1}, $earth)($V_1, $E_1)
+    dec_op_suite["Wedge Product"]["Old Form-0, Form-1"] = @benchmarkable wedge_product(Tuple{0,1}, $earth, $V_1, $E_1)
 
-    suite["Wedge Product"]["New Form-1, Form-1"] = @benchmarkable dec_wedge_product(Tuple{1,1}, $earth)($E_1, $E_2)
-    suite["Wedge Product"]["Old Form-1, Form-1"] = @benchmarkable wedge_product(Tuple{1,1}, $earth, $E_1, $E_2)
+    dec_op_suite["Wedge Product"]["New Form-1, Form-1"] = @benchmarkable dec_wedge_product(Tuple{1,1}, $earth)($E_1, $E_2)
+    dec_op_suite["Wedge Product"]["Old Form-1, Form-1"] = @benchmarkable wedge_product(Tuple{1,1}, $earth, $E_1, $E_2)
 
-    suite["Wedge Product"]["New Form-0, Form-2"] = @benchmarkable dec_wedge_product(Tuple{0,2}, $earth)($V_1, $T_2)
-    suite["Wedge Product"]["Old Form-0, Form-2"] = @benchmarkable wedge_product(Tuple{0,2}, $earth, $V_1, $T_2)
+    dec_op_suite["Wedge Product"]["New Form-0, Form-2"] = @benchmarkable dec_wedge_product(Tuple{0,2}, $earth)($V_1, $T_2)
+    dec_op_suite["Wedge Product"]["Old Form-0, Form-2"] = @benchmarkable wedge_product(Tuple{0,2}, $earth, $V_1, $T_2)
 end
 
 begin
@@ -125,25 +113,25 @@ begin
     E_1, E_2 = rand(float_type, ne(earth)), rand(float_type, ne(earth))
     T_2 = rand(float_type, ntriangles(earth))
 
-    suite["Wedge Product Computation"] = BenchmarkGroup()
+    dec_op_suite["Wedge Product Computation"] = BenchmarkGroup()
 
     wdg01 = dec_wedge_product(Tuple{0,1}, earth)
     wdg11 = dec_wedge_product(Tuple{1,1}, earth)
     wdg02 = dec_wedge_product(Tuple{0,2}, earth)
 
-    suite["Wedge Product Computation"]["New Form-0, Form-1"] = @benchmarkable wdg01($V_1, $E_1)
-    suite["Wedge Product Computation"]["New Form-1, Form-1"] = @benchmarkable wdg11($E_1, $E_2)
-    suite["Wedge Product Computation"]["New Form-0, Form-2"] = @benchmarkable wdg02($V_1, $T_2)
+    dec_op_suite["Wedge Product Computation"]["New Form-0, Form-1"] = @benchmarkable wdg01($V_1, $E_1)
+    dec_op_suite["Wedge Product Computation"]["New Form-1, Form-1"] = @benchmarkable wdg11($E_1, $E_2)
+    dec_op_suite["Wedge Product Computation"]["New Form-0, Form-2"] = @benchmarkable wdg02($V_1, $T_2)
 end
 
+# tune!(dec_op_suite)
 
+@info "Running DEC Operator Benchmarks"
 
-# tune!(suite)
-
-results = run(suite, verbose = true, seconds = 1)
+dec_op_results = run(dec_op_suite, verbose = true, seconds = 1)
 
-for op in sort(collect(keys(results)))
-    test = median(results[op])
+for op in sort(collect(keys(dec_op_results)))
+    test = median(dec_op_results[op])
 
     println("Operator: $op")
     for k in sort(collect(keys(test)))
@@ -154,4 +142,44 @@ for op in sort(collect(keys(results)))
     println("----------------------------------------------------------------")
 end
 
-end
\ No newline at end of file
+@info "Beginning Dual Mesh Generation Benchmarks"
+begin
+    mesh_size = 100
+    grid_spacings = [1.0, 0.8]#, 0.5, 0.4, 0.25, 0.2]
+    point_type = Point2D
+    benchmark_dual_meshes = map(gs -> triangulated_grid(mesh_size, mesh_size, gs, gs, point_type), grid_spacings);
+end;
+@info "Generated Primal Meshes"
+
+function create_dual_mesh(s, point_type, center_type)
+    sd = EmbeddedDeltaDualComplex2D{Bool,Float64, point_type}(s)
+    subdivide_duals!(sd, center_type)
+    sd;
+end
+
+dual_mesh_suite = BenchmarkGroup()
+filler_name = "Grid Spacing"
+begin
+    for (i, grid_spacing) in enumerate(grid_spacings)
+        dual_mesh_suite["Barycenter"][filler_name][grid_spacing] = @benchmarkable create_dual_mesh($(benchmark_dual_meshes[i]), $point_type, $(Barycenter())) gcsample=true seconds=10
+        dual_mesh_suite["Circumcenter"][filler_name][grid_spacing] = @benchmarkable create_dual_mesh($(benchmark_dual_meshes[i]), $point_type, $(Circumcenter())) gcsample=true seconds=10
+    end
+end
+
+@info "Running Dual Mesh Benchmarks"
+
+dual_mesh_results = run(dual_mesh_suite, verbose = true)
+
+for center in sort(collect(keys(dual_mesh_results)))
+    trial = median(dual_mesh_results[center][filler_name])
+
+    println("Center: $center")
+    for k in sort(collect(keys(trial)), rev=true)
+        t = trial[k].time / 1e9
+        m = trial[k].memory / 1e9
+        println("$k, [$t s, $m GB]")
+    end
+    println("----------------------------------------------------------------")
+end
+
+end
diff --git a/test/DiscreteExteriorCalculus.jl b/test/DiscreteExteriorCalculus.jl
index 8d266b3f..d09e73ff 100644
--- a/test/DiscreteExteriorCalculus.jl
+++ b/test/DiscreteExteriorCalculus.jl
@@ -141,6 +141,7 @@ dual_vs = elementary_duals(2,s,2)
 @test s[elementary_duals(1,s,3), :D_∂v1] == repeat([edge_center(s,3)], 2)
 @test [length(elementary_duals(s, V(i))) for i in 1:4] == [4,2,4,2]
 @test dual_triangle_vertices(s, 1) == [1,7,10]
+@test dual_edge_vertices(s, 1) == [5,2]
 
 # 2D oriented dual complex
 #-------------------------
@@ -492,7 +493,7 @@ end
 onedx = eval_constant_form(s, @SVector [1.0,0.0,0.0])
 onedy = eval_constant_form(s, @SVector [0.0,1.0,0.0])
 @test ∧(s, onedx, onedy) == map(s[:tri_orientation], s[:area]) do o,a
-  # Note by the order of -1 and 1 here that 
+  # Note by the order of -1 and 1 here that
   a * (o ? -1 : 1)
 end
 
@@ -658,4 +659,3 @@ for (primal_s,s) in flat_meshes
 end
 
 end
-