diff --git a/Project.toml b/Project.toml
index 2b4ade63..b64c49c8 100644
--- a/Project.toml
+++ b/Project.toml
@@ -9,6 +9,7 @@ ACSets = "227ef7b5-1206-438b-ac65-934d6da304b8"
 Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 Catlab = "134e5e36-593f-5add-ad60-77f754baafbe"
+DataMigrations = "0c4ad18d-0c49-4bc2-90d5-5bca8f00d6ae"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
diff --git a/src/DiscreteExteriorCalculus.jl b/src/DiscreteExteriorCalculus.jl
index e827bfc5..92dc7076 100644
--- a/src/DiscreteExteriorCalculus.jl
+++ b/src/DiscreteExteriorCalculus.jl
@@ -18,21 +18,30 @@ export DualSimplex, DualV, DualE, DualTri, DualChain, DualForm,
   SimplexCenter, Barycenter, Circumcenter, Incenter, geometric_center,
   subsimplices, primal_vertex, elementary_duals, dual_boundary, dual_derivative,
   ⋆, hodge_star, inv_hodge_star, δ, codifferential, ∇², laplace_beltrami, Δ, laplace_de_rham,
-  ♭, flat, ♯, sharp, ∧, wedge_product, interior_product, interior_product_flat,
+  ♭, 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,
-  dual_point, dual_volume, subdivide_duals!, DiagonalHodge, GeometricHodge
+  dual_point, dual_volume, subdivide_duals!, DiagonalHodge, GeometricHodge,
+  subdivide, PPSharp, AltPPSharp, DesbrunSharp, LLSDDSharp, de_sign,
+  ♭♯, ♭♯_mat, flat_sharp, flat_sharp_mat
 
 import Base: ndims
-using LinearAlgebra: Diagonal, dot, norm, cross
+import Base: *
+import LinearAlgebra: mul!
+using LinearAlgebra: Diagonal, dot, norm, cross, pinv
 using SparseArrays
-using StaticArrays: @SVector, SVector
+using StaticArrays: @SVector, SVector, SMatrix
+using GeometryBasics: Point2, Point3
+
+const Point2D = SVector{2,Float64}
+const Point3D = SVector{3,Float64}
 
 using ACSets.DenseACSets: attrtype_type
 using Catlab, Catlab.CategoricalAlgebra.CSets
 using Catlab.CategoricalAlgebra.FinSets: deleteat
 import Catlab.CategoricalAlgebra.CSets: ∧
 import Catlab.Theories: Δ
+using DataMigrations: @migrate
 
 using ..ArrayUtils, ..SimplicialSets
 using ..SimplicialSets: CayleyMengerDet, operator_nz, ∂_nz, d_nz,
@@ -44,11 +53,65 @@ struct DPPFlat <: DiscreteFlat end
 
 abstract type DiscreteSharp end
 struct PPSharp <: DiscreteSharp end
+struct AltPPSharp <: DiscreteSharp end
+struct DesbrunSharp <: DiscreteSharp end
+struct LLSDDSharp <: DiscreteSharp end
 
 abstract type DiscreteHodge end
 struct GeometricHodge <: DiscreteHodge end
 struct DiagonalHodge  <: DiscreteHodge end
 
+# Euclidean geometry
+####################
+
+""" A notion of "geometric center" of a simplex.
+
+See also: [`geometric_center`](@ref).
+"""
+abstract type SimplexCenter end
+
+""" Calculate the center of simplex spanned by given points.
+
+The first argument is a list of points and the second specifies the notion of
+"center", via an instance of [`SimplexCenter`](@ref).
+"""
+function geometric_center end
+
+""" Barycenter, aka centroid, of a simplex.
+"""
+struct Barycenter <: SimplexCenter end
+
+function geometric_center(points, ::Barycenter)
+  sum(points) / length(points)
+end
+
+""" Circumcenter, or center of circumscribed circle, of a simplex.
+
+The circumcenter is calculated by inverting the Cayley-Menger matrix, as
+explained by
+[Westdendorp](https://westy31.home.xs4all.nl/Circumsphere/ncircumsphere.htm).
+This method of calculation is also used in the package
+[AlphaShapes.jl](https://github.com/harveydevereux/AlphaShapes.jl).
+"""
+struct Circumcenter <: SimplexCenter end
+
+function geometric_center(points, ::Circumcenter)
+  CM = cayley_menger(points...)
+  barycentric_coords = inv(CM)[1,2:end]
+  mapreduce(*, +, barycentric_coords, points)
+end
+
+""" Incenter, or center of inscribed circle, of a simplex.
+"""
+struct Incenter <: SimplexCenter end
+
+function geometric_center(points, ::Incenter)
+  length(points) > 2 || return geometric_center(points, Barycenter())
+  face_volumes = map(i -> volume(deleteat(points, i)), eachindex(points))
+  barycentric_coords = face_volumes / sum(face_volumes)
+  mapreduce(*, +, barycentric_coords, points)
+end
+
 # 1D dual complex
 #################
 
@@ -121,6 +184,29 @@ function dual_triangle_vertices(s::HasDeltaSet1D, t...)
           s[s[t..., :D_∂e0], :D_∂v0])
 end
 
+""" Subdivide a 1D delta set.
+"""
+function subdivide(s::HasDeltaSet1D)
+  @migrate typeof(s) s begin
+    V => @cases begin
+      v::V
+      e::E
+    end
+    E => @cases begin
+      e₁::E
+      e₂::E
+    end
+    ∂v1 => begin
+      e₁ => e
+      e₂ => e
+    end
+    ∂v0 => begin
+      e₁ => (v∘∂v1)
+      e₂ => (v∘∂v0)
+    end
+  end
+end
+
 # 1D oriented dual complex
 #-------------------------
 
@@ -195,6 +281,41 @@ function make_dual_simplices_1d!(s::HasDeltaSet1D, ::Type{Simplex{n}}) where n
   D_edges
 end
 
+# TODO: Instead of copying-and-pasting the DeltaSet1D version:
+# - Use metaprogramming, or
+# - Don't use the migration DSL, but rather the lower-level functor interface.
+# TODO: When Catlab PR #823 "Data migrations with Julia functions on attributes"
+# is merged, ensure that oriented-ness is preserved. (Flip one of the
+# orientations.)
+""" Subdivide an oriented 1D delta set.
+
+Note that this function does NOT currently guarantee that if the input is
+oriented, then the output will be.
+"""
+function subdivide(s::OrientedDeltaSet1D{T}) where T
+  @migrate typeof(s) s begin
+    V => @cases begin
+      v::V
+      e::E
+    end
+    E => @cases begin
+      e₁::E
+      e₂::E
+    end
+    ∂v1 => begin
+      e₁ => e
+      e₂ => e
+    end
+    ∂v0 => begin
+      e₁ => (v∘∂v1)
+      e₂ => (v∘∂v0)
+    end
+    Orientation => Orientation
+    # TODO: One of these edge orientations must be flipped. (e₂?)
+    edge_orientation => (e₁ => edge_orientation; e₂ => edge_orientation)
+  end
+end
+
 # 1D embedded dual complex
 #-------------------------
 
@@ -267,6 +388,18 @@ function precompute_volumes_1d!(s::HasDeltaSet1D)
   end
 end
 
+# TODO: When Catlab PR #823 "Data migrations with Julia functions on attributes"
+# is merged, encode subdivision like so:
+#function subdivide(s::EmbeddedDeltaSet1D{T,U}, alg::V) where {T,U,V <: SimplexCenter}
+#  @migrate typeof(s) s begin
+#    ...
+#    edge_orientation => (e₁ => edge_orientation; e₂ => !(edge_orientation))
+#    Point => Point
+#    point => (v => point; e => geometric_center([e₁ ⋅ point, e₂ ⋅ point], alg))
+#    ...
+#  end
+#end
+
 # 2D dual complex
 #################
 
@@ -435,7 +568,7 @@ relative_sign(x::Bool, y::Bool) = (x && y) || (!x && !y)
   dual_area::Attr(DualTri, Real)
 end
 
-""" Embedded dual complex of an embedded 12 delta set.
+""" Embedded dual complex of an embedded 2D delta set.
 
 Although they are redundant information, the lengths and areas of the
 primal/dual edges and triangles are precomputed and stored.
@@ -473,21 +606,62 @@ function ♭(s::AbstractDeltaDualComplex2D, X::AbstractVector, ::DPPFlat)
   # happens to be inconvenient.
   tri_map = Dict{Int,Int}(triangle_center(s,t) => t for t in triangles(s))
 
+  # TODO: Remove these comments before merging.
+  # For each primal edge:
   map(edges(s)) do e
+    # Get the vector from src to tgt (oriented correctly).
     e_vec = (point(s, tgt(s,e)) - point(s, src(s,e))) * sign(1,s,e)
+    # Grab all the dual edges of this primal edge.
     dual_edges = elementary_duals(1,s,e)
+    # And the corresponding lengths.
     dual_lengths = dual_volume(1, s, dual_edges)
+    # For each of these dual edges:
     mapreduce(+, dual_edges, dual_lengths) do dual_e, dual_length
+      # Get the vector at the center of the triangle this edge is pointing at.
       X_vec = X[tri_map[s[dual_e, :D_∂v0]]]
+      # Take their dot product and multiply by the length of this dual edge.
       dual_length * dot(X_vec, e_vec)
+      # When done, sum these weights up and divide by the total length.
     end / sum(dual_lengths)
   end
 end
 
-function ♯(s::AbstractDeltaDualComplex2D, α::AbstractVector, ::PPSharp)
+function ♭_mat(s::AbstractDeltaDualComplex2D)
+  ♭_mat(s, ∂(2,s))
+end
+
+function ♭_mat(s::AbstractDeltaDualComplex2D, p2s)
+  mat_type = SMatrix{1, length(eltype(s[:point])), eltype(eltype(s[:point])), length(eltype(s[:point]))}
+  ♭_mat = spzeros(mat_type, ne(s), ntriangles(s))
+  for e in edges(s)
+    # The vector associated with this primal edge.
+    e_vec = (point(s, tgt(s,e)) - point(s, src(s,e))) * sign(1,s,e)
+    # The triangles associated with this primal edge.
+    tris = p2s[e,:].nzind
+    # The dual vertex at the center of this primal edge.
+    center = edge_center(s, e)
+    # The centers of the triangles associated with this primal edge.
+    dvs = triangle_center(s, tris)
+    # The dual edge pointing to each triangle.
+    des = map(dvs) do dv
+      # (This is the edges(s,src,tgt) function.)
+      only(de for de in incident(s, dv, :D_∂v0) if s[de, :D_∂v1] == center)
+    end
+    # The lengths of those dual edges.
+    dels = volume(s, DualE(des))
+    # The sum of the lengths of the dual edges at each primal edge.
+    dels_sum = sum(dels)
+
+    for (tri, del) in zip(tris, dels)
+      ♭_mat[e, tri] = del * mat_type(e_vec) / dels_sum
+    end
+  end
+  ♭_mat
+end
+
+function ♯(s::AbstractDeltaDualComplex2D, α::AbstractVector, DS::DiscreteSharp)
   α♯ = zeros(attrtype_type(s, :Point), nv(s))
   for t in triangles(s)
-    area = volume(2,s,t)
     tri_center, tri_edges = triangle_center(s,t), triangle_edges(s,t)
     tri_point = dual_point(s, tri_center)
     for (i, (v₀, e₀)) in enumerate(zip(triangle_vertices(s,t), tri_edges))
@@ -501,6 +675,7 @@ function ♯(s::AbstractDeltaDualComplex2D, α::AbstractVector, ::PPSharp)
         v, sgn = src(s,e) == v₀ ? (tgt(s,e), -1) : (src(s,e), +1)
         dual_area = sum(dual_volume(2,s,d) for d in elementary_duals(0,s,v)
                         if s[s[d, :D_∂e0], :D_∂v0] == tri_center)
+        area = ♯_denominator(s, v, t, DS)
         α♯[v] += sgn * sign(1,s,e) * α[e] * (dual_area / area) * out_vec
       end
     end
@@ -508,6 +683,132 @@ function ♯(s::AbstractDeltaDualComplex2D, α::AbstractVector, ::PPSharp)
   α♯
 end
 
+function ♯(s::AbstractDeltaDualComplex2D, α::AbstractVector, ::LLSDDSharp)
+  ♯_m = ♯_mat(s, LLSDDSharp())
+  ♯_m * α
+end
+
+""" Divided weighted normals by | σⁿ | .
+
+This weighting is that used in equation 5.8.1 from Hirani.
+
+See Hirani §5.8.
+"""
+♯_denominator(s::AbstractDeltaDualComplex2D, _::Int, t::Int, ::DiscreteSharp) =
+  volume(2,s,t)
+
+""" Divided weighted normals by | ⋆v | .
+
+This weighting is NOT that of equation 5.8.1, but a different weighting scheme.
+We essentially replace the denominator in equation 5.8.1 with | ⋆v | . This
+may be what Hirani intended, and perhaps the denominator | σⁿ | in that equation
+is either a mistake or clerical error.
+
+See Hirani §5.8.
+"""
+♯_denominator(s::AbstractDeltaDualComplex2D, v::Int, _::Int, ::AltPPSharp) =
+  sum(dual_volume(2,s, elementary_duals(0,s,v)))
+
+"""    function get_orthogonal_vector(s::AbstractDeltaDualComplex2D, v::Int, e::Int)
+
+Find a vector orthogonal to e pointing into the triangle shared with v.
+"""
+function get_orthogonal_vector(s::AbstractDeltaDualComplex2D, v::Int, e::Int)
+  e_vec = point(s, tgt(s, e)) - point(s, src(s, e))
+  e_vec /= norm(e_vec)
+  e2_vec = point(s, v) - point(s, src(s, e))
+  e2_vec - dot(e2_vec, e_vec)*e_vec
+end
+
+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)
+    v, sgn = src(s,e) == v₀ ? (tgt(s,e), -1) : (src(s,e), +1)
+    # | ⋆vₓ ∩ σⁿ |
+    dual_area = sum(dual_volume(2,s,d) for d in elementary_duals(0,s,v)
+                    if s[s[d, :D_∂e0], :D_∂v0] == tri_center)
+    area = ♯_denominator(s, v, t, DS)
+    ♯_mat[v,e] += sgn * sign(1,s,e) * (dual_area / area) * out_vec
+  end
+end
+
+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₀)
+    sgn = v == tgt(s,e₀) ? -1 : +1
+    # | ⋆vₓ ∩ σⁿ |
+    dual_area = sum(dual_volume(2,s,d) for d in elementary_duals(0,s,v)
+                    if s[s[d, :D_∂e0], :D_∂v0] == tri_center)
+    area = ♯_denominator(s, v, t, DS)
+    ♯_mat[v,e₀] += sgn * sign(1,s,e₀) * (dual_area / area) * out_vec
+  end
+end
+
+"""    function ♯_mat(s::AbstractDeltaDualComplex2D, DS::DiscreteSharp)
+
+Sharpen a 1-form into a vector field.
+
+3 primal-primal methods are supported. See [`♯_denominator`](@ref) for the distinction between Hirani's method and and an "Alternative" method. Desbrun's definition is selected with `DesbrunSharp`, and is like Hirani's, save for dividing by the norm twice.
+
+A dual-dual method which uses linear least squares to estimate a vector field is selected with `LLSDDSharp`.
+"""
+function ♯_mat(s::AbstractDeltaDualComplex2D, DS::DiscreteSharp)
+  ♯_mat = spzeros(attrtype_type(s, :Point), (nv(s), ne(s)))
+  for t in triangles(s)
+    tri_center, tri_edges = triangle_center(s,t), triangle_edges(s,t)
+    for (i, (v₀, e₀)) in enumerate(zip(triangle_vertices(s,t), tri_edges))
+      out_vec = get_orthogonal_vector(s, v₀, e₀)
+      h = norm(out_vec)
+      out_vec /= DS == DesbrunSharp() ? h : h^2
+      ♯_assign!(♯_mat, s, v₀, e₀, t, i, tri_edges, tri_center, out_vec, DS)
+    end
+  end
+  ♯_mat
+end
+
+de_sign(s,de) = s[de, :D_edge_orientation] ? +1 : -1
+
+"""    function ♯_mat(s::AbstractDeltaDualComplex2D, ::LLSDDSharp)
+
+Sharpen a dual 1-form into a DualVectorField, using linear least squares.
+
+Up to floating point error, this method perfectly produces fields which are constant over any triangle in the domain. Assume that the contribution of each half-edge to the value stored on the entire dual edge is proportional to their lengths. Since this least squares method does not perform pre-normalization, the contribution of each half-edge value is proportional to its length on the given triangle. Satisfying the continuous exterior calculus, sharpened vectors are constrained to lie on their triangle (i.e. they are indeed tangent).
+
+It is not known whether this method has been exploited previously in the DEC literature, or defined in code elsewhere.
+"""
+function ♯_mat(s::AbstractDeltaDualComplex2D, ::LLSDDSharp)
+  # TODO: Grab point information out of s at the type level.
+  pt = attrtype_type(s, :Point)
+  ♯_m = spzeros(SVector{length(pt), eltype(pt)},
+                findnz(d(1,s))[[1,2]]...)
+  for t in triangles(s)
+    tri_center, tri_edges = triangle_center(s,t), sort(triangle_edges(s,t))
+    # | ⋆eₓ ∩ σⁿ |
+    star_e_cap_t = map(tri_edges) do e
+      only(filter(elementary_duals(1,s,e)) do de
+        s[de, :D_∂v0] == tri_center
+      end)
+    end
+    de_vecs = map(star_e_cap_t) do de
+      de_sign(s,de) *
+        (dual_point(s,s[de, :D_∂v0]) - dual_point(s,s[de, :D_∂v1]))
+    end
+    weights = s[star_e_cap_t, :dual_length] ./
+      map(tri_edges) do e
+        sum(s[elementary_duals(1,s,e), :dual_length])
+      end
+    # TODO: Move around ' as appropriate to minimize transposing.
+    X = stack(de_vecs)'
+    LLS = pinv(X'*(X))*(X')
+    for (i,e) in enumerate(tri_edges)
+      ♯_m[t, e] = LLS[:,i]'*weights[i]
+    end
+  end
+  ♯_m
+end
+
 function ∧(::Type{Tuple{1,1}}, s::HasDeltaSet2D, α, β, x::Int)
   # XXX: This calculation of the volume coefficients is awkward due to the
   # design decision described in `SchDeltaDualComplex1D`.
@@ -599,6 +900,36 @@ this correspondence, a basis for primal ``n``-chains defines the basis for dual
 """
 @vector_struct DualVectorField
 
+# When the user wraps a vector of SVectors in DualVectorField(), automatically
+# reinterpret the result of multiplication from a vector of length 1 SVectors
+# to a vector of floats.
+function *(A::AbstractMatrix{T}, x::DualVectorField) where T
+  reinterpret(Float64, invoke(*, Tuple{AbstractMatrix{T} where T, AbstractVector{S} where S}, A, x))
+end
+
+function mul!(C::Vector{H}, A::AbstractVecOrMat,
+  B::DualVectorField, α::Number, β::Number) where {H <: Number}
+  size(A, 2) == size(B, 1) || throw(DimensionMismatch())
+  size(A, 1) == size(C, 1) || throw(DimensionMismatch())
+  size(B, 2) == size(C, 2) || throw(DimensionMismatch())
+  nzv = nonzeros(A)
+  rv = rowvals(A)
+  if β != 1
+      #β != 0 ? rmul!(C, β) : fill!(C, zero(eltype(C)))
+      β != 0 ? rmul!(C, β) : fill!(C, zero(H))
+  end
+  for k in 1:size(C, 2)
+      @inbounds for col in 1:size(A, 2)
+          αxj = B[col,k] * α
+          for j in nzrange(A, col)
+              #C[rv[j], k] += nzv[j]*αxj
+              C[rv[j], k] += only(nzv[j]*αxj)
+          end
+      end
+  end
+  C
+end
+
 ndims(s::AbstractDeltaDualComplex1D) = 1
 ndims(s::AbstractDeltaDualComplex2D) = 2
 
@@ -609,7 +940,7 @@ volume(s::HasDeltaSet, x::DualSimplex{n}, args...) where n =
 
 """ List of dual simplices comprising the subdivision of a primal simplex.
 
-A primal ``n``-simplex is always subdivided into ``n!`` dual ``n``-simplices,
+A primal ``n``-simplex is always subdivided into ``(n+1)!`` dual ``n``-simplices,
 not be confused with the [`elementary_duals`](@ref) which have complementary
 dimension.
 
@@ -641,7 +972,7 @@ In 2D dual complexes, the elementary duals of...
 
 - primal vertices are dual triangles
 - primal edges are dual edges
-- primal triangles are (single) dual triangles
+- primal triangles are (single) dual vertices
 """
 elementary_duals(s::HasDeltaSet, x::Simplex{n}) where n =
   DualSimplex{ndims(s)-n}(elementary_duals(Val{n}, s, x.data))
@@ -969,11 +1300,9 @@ See also: the sharp operator [`♯`](@ref).
 """
 const flat = ♭
 
-""" Sharp operator for converting 1-forms to vector fields.
+""" Sharp operator for converting primal 1-forms to primal vector fields.
 
-A generic function for discrete sharp operators. Currently only the
-primal-primal flat from (Hirani 2003, Definition 5.8.1 and Remark 2.7.2) is
-implemented.
+This the primal-primal sharp from Hirani 2003, Definition 5.8.1 and Remark 2.7.2.
 
 !!! note
 
@@ -987,14 +1316,49 @@ implemented.
     knowledge, our implementation is the correct one and agrees with Hirani's
     description, if not his figure.
 
-See also: the flat operator [`♭`](@ref).
+See also: [`♭`](@ref) and [`♯_mat`](@ref), which returns a matrix that encodes this operator.
 """
 ♯(s::HasDeltaSet, α::EForm) = PrimalVectorField(♯(s, α.data, PPSharp()))
 
+""" Sharp operator for converting dual 1-forms to dual vector fields.
+
+This dual-dual sharp uses a method of local linear least squares to provide a
+tangent vector field.
+
+See also: [`♯_mat`](@ref), which returns a matrix that encodes this operator.
+"""
+♯(s::HasDeltaSet, α::DualForm{1}) = DualVectorField(♯(s, α.data, LLSDDSharp()))
+
 """ Alias for the sharp operator [`♯`](@ref).
 """
 const sharp = ♯
 
+"""    ♭♯_mat(s::HasDeltaSet)
+
+Make a dual 1-form primal by chaining ♭ᵈᵖ♯ᵈᵈ.
+
+This returns a matrix which can be multiplied by a dual 1-form.
+See also [`♭♯`](@ref).
+"""
+♭♯_mat(s::HasDeltaSet) = ♭_mat(s) * ♯_mat(s, LLSDDSharp())
+
+"""    ♭♯(s::HasDeltaSet, α::SimplexForm{1})
+
+Make a dual 1-form primal by chaining ♭ᵈᵖ♯ᵈᵈ.
+
+This returns the given dual 1-form as a primal 1-form.
+See also [`♭♯_mat`](@ref).
+"""
+♭♯(s::HasDeltaSet, α::SimplexForm{1}) = only.(♭♯_mat(s) * α)
+
+""" Alias for the flat-sharp dual-to-primal interpolation operator [`♭♯`](@ref).
+"""
+const flat_sharp = ♭♯
+
+""" Alias for the flat-sharp dual-to-primal interpolation matrix [`♭♯_mat`](@ref).
+"""
+const flat_sharp_mat = ♭♯_mat
+
 """ Wedge product of discrete forms.
 
 The wedge product of a ``k``-form and an ``l``-form is a ``(k+l)``-form.
@@ -1002,7 +1366,8 @@ The wedge product of a ``k``-form and an ``l``-form is a ``(k+l)``-form.
 The DEC and related systems have several flavors of wedge product. This one is
 the discrete primal-primal wedge product introduced in (Hirani, 2003, Chapter 7)
 and (Desbrun et al 2005, Section 8). It depends on the geometric embedding and
-requires the dual complex.
+requires the dual complex. Note that we diverge from Hirani in that his
+formulation explicitly divides by (k+1)!. We do not do so in this computation.
 """
 ∧(s::HasDeltaSet, α::SimplexForm{k}, β::SimplexForm{l}) where {k,l} =
   SimplexForm{k+l}(∧(Tuple{k,l}, s, α.data, β.data))
@@ -1027,7 +1392,7 @@ function wedge_product_zero(::Type{Val{k}}, s::HasDeltaSet,
   subs = subsimplices(k, s, x)
   vs = primal_vertex(k, s, subs)
   coeffs = map(x′ -> dual_volume(k,s,x′), subs) / volume(k,s,x)
-  dot(coeffs, f[vs]) * α[x] / factorial(k)
+  dot(coeffs, f[vs]) * α[x]
 end
 
 """ Alias for the wedge product operator [`∧`](@ref).
@@ -1097,55 +1462,26 @@ function lie_derivative_flat(::Type{Val{2}}, s::HasDeltaSet,
   dual_derivative(1, s, interior_product_flat(2, s, X♭, α; kw...))
 end
 
-# Euclidean geometry
-####################
-
-""" A notion of "geometric center" of a simplex.
-
-See also: [`geometric_center`](@ref).
-"""
-abstract type SimplexCenter end
-
-""" Calculate the center of simplex spanned by given points.
-
-The first argument is a list of points and the second specifies the notion of
-"center", via an instance of [`SimplexCenter`](@ref).
-"""
-function geometric_center end
-
-""" Barycenter, aka centroid, of a simplex.
-"""
-struct Barycenter <: SimplexCenter end
-
-function geometric_center(points, ::Barycenter)
-  sum(points) / length(points)
+function eval_constant_primal_form(s::EmbeddedDeltaDualComplex2D{Bool, Float64, T} where T<:Union{Point3D, Point3{Float64}}, α::SVector{3,Float64})
+  EForm(map(edges(s)) do e
+          dot(α, point(s, tgt(s,e)) - point(s, src(s,e))) * sign(1,s,e)
+        end)
 end
-
-""" Circumcenter, or center of circumscribed circle, of a simplex.
-
-The circumcenter is calculated by inverting the Cayley-Menger matrix, as
-explained by
-[Westdendorp](https://westy31.home.xs4all.nl/Circumsphere/ncircumsphere.htm).
-This method of calculation is also used in the package
-[AlphaShapes.jl](https://github.com/harveydevereux/AlphaShapes.jl).
-"""
-struct Circumcenter <: SimplexCenter end
-
-function geometric_center(points, ::Circumcenter)
-  CM = cayley_menger(points...)
-  barycentric_coords = inv(CM)[1,2:end]
-  mapreduce(*, +, barycentric_coords, points)
+function eval_constant_primal_form(s::EmbeddedDeltaDualComplex2D{Bool, Float64, T} where T<:Union{Point2D, Point2{Float64}}, α::SVector{3,Float64})
+  α = SVector{2,Float64}(α[1],α[2])
+  EForm(map(edges(s)) do e
+          dot(α, point(s, tgt(s,e)) - point(s, src(s,e))) * sign(1,s,e)
+        end)
 end
 
-""" Incenter, or center of inscribed circle, of a simplex.
-"""
-struct Incenter <: SimplexCenter end
-
-function geometric_center(points, ::Incenter)
-  length(points) > 2 || return geometric_center(points, Barycenter())
-  face_volumes = map(i -> volume(deleteat(points, i)), eachindex(points))
-  barycentric_coords = face_volumes / sum(face_volumes)
-  mapreduce(*, +, barycentric_coords, points)
+# Evaluate a constant dual form
+# XXX: This "left/right-hand-rule" trick only works when z=0.
+# XXX: So, do not use this function to test e.g. curved surfaces.
+function eval_constant_dual_form(s::EmbeddedDeltaDualComplex2D, α::SVector{3,Float64})
+  EForm(
+    hodge_star(1,s) *
+      eval_constant_primal_form(s, SVector{3,Float64}(α[2], -α[1], α[3])))
 end
 
+
 end
diff --git a/src/FastDEC.jl b/src/FastDEC.jl
index 96a90f98..854e35ce 100644
--- a/src/FastDEC.jl
+++ b/src/FastDEC.jl
@@ -1,14 +1,19 @@
 module FastDEC
 using LinearAlgebra: Diagonal, dot, norm, cross
 using StaticArrays: SVector
-using SparseArrays: sparse
+using SparseArrays: sparse, spzeros
 using LinearAlgebra
 using Base.Iterators
 using ACSets.DenseACSets: attrtype_type
 using Catlab, Catlab.CategoricalAlgebra.CSets
 using ..SimplicialSets, ..DiscreteExteriorCalculus
+import ..DiscreteExteriorCalculus: ∧
 
-export dec_boundary, dec_differential, dec_dual_derivative, dec_hodge_star, dec_inv_hodge_star, dec_wedge_product, dec_c_wedge_product, dec_p_wedge_product, dec_c_wedge_product!
+export dec_boundary, dec_differential, dec_dual_derivative, dec_hodge_star, dec_inv_hodge_star, dec_wedge_product, dec_c_wedge_product, dec_p_wedge_product, dec_c_wedge_product!,
+       dec_wedge_product_pd, dec_wedge_product_dp, ∧,
+       interior_product_dd, ℒ_dd,
+       dec_wedge_product_dd,
+       Δᵈ
 
 """
     dec_p_wedge_product(::Type{Tuple{0,1}}, sd::EmbeddedDeltaDualComplex1D)
@@ -230,6 +235,150 @@ function dec_wedge_product(::Type{Tuple{1,1}}, sd::HasDeltaSet2D)
     (α, β) -> dec_c_wedge_product(Tuple{1,1}, α, β, val_pack)
 end
 
+"""
+    function wedge_dd_01_mat(sd::HasDeltaSet)
+
+Returns a matrix that can be multiplied to a dual 0-form, before being
+elementwise-multiplied by a dual 1-form, encoding the wedge product.
+"""
+function wedge_dd_01_mat(sd::HasDeltaSet)
+  m = spzeros(ne(sd), ntriangles(sd))
+  for e in edges(sd)
+    des = elementary_duals(1,sd,e)
+    dvs = sd[des, :D_∂v0]
+    tris = only.(incident(sd, dvs, :tri_center))
+    ws = sd[des, :dual_length] ./ sum(sd[des, :dual_length])
+    for (w,t) in zip(ws,tris)
+      m[e,t] = w
+    end
+  end
+  m
+end
+
+"""
+    dec_wedge_product_dd(::Type{Tuple{0,1}}, sd::HasDeltaSet)
+
+Returns a cached function that computes the wedge product between a dual
+0-form and a dual 1-form.
+"""
+function dec_wedge_product_dd(::Type{Tuple{0,1}}, sd::HasDeltaSet)
+  m = wedge_dd_01_mat(sd)
+  (f,g) -> (m * f) .* g
+end
+
+"""
+    dec_wedge_product_dd(::Type{Tuple{1,0}}, sd::HasDeltaSet)
+
+Returns a cached function that computes the wedge product between a dual
+1-form and a dual 0-form.
+"""
+function dec_wedge_product_dd(::Type{Tuple{1,0}}, sd::HasDeltaSet)
+  m = wedge_dd_01_mat(sd)
+  (f,g) -> f .* (m * g)
+end
+
+"""
+    function wedge_pd_01_mat(sd::HasDeltaSet)
+
+Returns a matrix that can be multiplied to a primal 0-form, before being
+elementwise-multiplied by a dual 1-form, encoding the wedge product.
+
+This function assumes barycentric means and performs bilinear interpolation. It
+is not known if this definition has appeared in the literature or any code.
+"""
+function wedge_pd_01_mat(sd::HasDeltaSet)
+  m = spzeros(ne(sd), nv(sd))
+  for e in edges(sd)
+    α, β = edge_vertices(sd,e)
+    des = elementary_duals(1,sd,e)
+    dvs = sd[des, :D_∂v0]
+    tris = only.(incident(sd, dvs, :tri_center))
+    γδ = map(tris) do t
+      only(filter(x -> x ∉ [α,β], triangle_vertices(sd,t)))
+    end
+    ws = sd[des, :dual_length] ./ sum(sd[des, :dual_length])
+    for (w,l) in zip(ws,γδ)
+      m[e,α] += w*5/12
+      m[e,β] += w*5/12
+      m[e,l] += w*2/12
+    end
+  end
+  m
+end
+
+"""
+    dec_wedge_product_pd(::Type{Tuple{1,0}}, sd::HasDeltaSet)
+
+Returns a cached function that computes the wedge product between a dual
+1-form and a primal 0-form.
+
+This function assumes barycentric means and performs bilinear interpolation. It
+is not known if this definition has appeared in the literature or any code.
+"""
+function dec_wedge_product_dp(::Type{Tuple{1,0}}, sd::HasDeltaSet)
+  m = wedge_pd_01_mat(sd)
+  (f,g) -> f .* (m * g)
+end
+
+"""
+    function dec_wedge_product_pd(::Type{Tuple{0,1}}, sd::HasDeltaSet)
+
+Returns a cached function that computes the wedge product between a primal
+0-form and a dual 1-form.
+
+This function assumes barycentric means and performs bilinear interpolation. It
+is not known if this definition has appeared in the literature or any code.
+"""
+function dec_wedge_product_pd(::Type{Tuple{0,1}}, sd::HasDeltaSet)
+  m = wedge_pd_01_mat(sd)
+  (g,f) -> (m * g) .* f
+end
+
+"""
+    dec_wedge_product_pd(::Type{Tuple{1,1}}, sd::HasDeltaSet)
+
+Returns a cached function that computes the wedge product between a primal
+1-form and a dual 1-form.
+"""
+function dec_wedge_product_pd(::Type{Tuple{1,1}}, sd::HasDeltaSet)
+  ♭♯_m = ♭♯_mat(sd)
+  ♭♯_cached(x) = only.(♭♯_m * x)
+  Λ_cached = dec_wedge_product(Tuple{1, 1}, sd)
+  (f, g) -> Λ_cached(f, ♭♯_cached(g))
+end
+
+"""
+    dec_wedge_product_dp(::Type{Tuple{1,1}}, sd::HasDeltaSet)
+
+Returns a cached function that computes the wedge product between a dual 1-form
+and a primal 1-form.
+"""
+function dec_wedge_product_dp(::Type{Tuple{1,1}}, sd::HasDeltaSet)
+  ♭♯_m = ♭♯_mat(sd)
+  ♭♯_cached(x) = only.(♭♯_m * x)
+  Λ_cached = dec_wedge_product(Tuple{1, 1}, sd)
+  (f, g) -> Λ_cached(♭♯_cached(f), g)
+end
+
+""" Wedge product of a primal 1-form and a dual 1-form.
+
+Chain the musical isomorphisms to interpolate the dual 1-form to a primal
+1-form, using the linear least squares ♯. Then use the CombinatorialSpaces
+version of the Hirani primal-primal weddge.
+"""
+∧(s::HasDeltaSet, α::SimplexForm{1}, β::DualForm{1}) =
+  dec_wedge_product_pd(Tuple{1,1}, s)(α, β)
+
+""" Wedge product of a dual 1-form and a primal 1-form.
+
+Chain the musical isomorphisms to interpolate the dual 1-form to a primal
+1-form. Then use the CombinatorialSpaces version of the Hirani primal-primal
+weddge (without explicitly dividing by 2.)
+"""
+∧(s::HasDeltaSet, α::DualForm{1}, β::SimplexForm{1}) =
+  dec_wedge_product_dp(Tuple{1,1}, s)(α, β)
+
+
 # Boundary Operators
 dec_boundary(n::Int, sd::HasDeltaSet) = sparse(dec_p_boundary(Val{n}, sd)...)
 
@@ -532,4 +681,82 @@ end
 
 dec_inv_hodge_star(::Type{Val{2}}, sd::EmbeddedDeltaDualComplex2D, ::GeometricHodge) =
     dec_inv_hodge_star(Val{2}, sd, DiagonalHodge())
+
+"""    function interior_product_dd(::Type{Tuple{1,1}}, s::SimplicialSets.HasDeltaSet)
+
+Given a dual 1-form and a dual 1-form, return their interior product as a dual 0-form.
+"""
+function interior_product_dd(::Type{Tuple{1,1}}, s::SimplicialSets.HasDeltaSet)
+  ihs1 = dec_inv_hodge_star(Val{1}, s, GeometricHodge())
+  Λ11 = dec_wedge_product_pd(Tuple{1,1}, s)
+  hs2 = dec_hodge_star(Val{2}, s, GeometricHodge())
+
+  (f,g) -> hs2 * Λ11(ihs1(g), f)
+end
+
+"""    function interior_product_dd(::Type{Tuple{1,1}}, s::SimplicialSets.HasDeltaSet)
+
+Given a dual 1-form and a dual 2-form, return their interior product as a dual 1-form.
+"""
+function interior_product_dd(::Type{Tuple{1,2}}, s::SimplicialSets.HasDeltaSet)
+  ihs0 = dec_inv_hodge_star(Val{0}, s, GeometricHodge())
+  hs1 = dec_hodge_star(Val{1}, s, GeometricHodge())
+  ♭♯_m = ♭♯_mat(s)
+  Λ01_m = wedge_pd_01_mat(s)
+  (f,g) -> hs1 * only.(♭♯_m * ((Λ01_m * ihs0 * g) .* f))
+end
+
+"""    function ℒ_dd(::Type{Tuple{1,1}}, s::SimplicialSets.HasDeltaSet)
+
+Given a dual 1-form and a dual 1-form, return their lie derivative as a dual 1-form.
+"""
+function ℒ_dd(::Type{Tuple{1,1}}, s::SimplicialSets.HasDeltaSet)
+  # TODO: Check signs.
+  # ℒ := -diuv - iduv
+  d0 = dec_dual_derivative(0, s)
+  d1 = dec_dual_derivative(1, s)
+  i1 = interior_product_dd(Tuple{1,1}, s)
+  i2 = interior_product_dd(Tuple{1,2}, s)
+
+  (f,g) ->
+    -(d0 * i1(f,g)) -
+      i2(f,d1 * g)
+end
+
+const lie_derivative_dd = ℒ_dd
+
+"""    function Δᵈ_mat(::Type{Val{0}}, s::SimplicialSets.HasDeltaSet)
+
+Return a function matrix encoding the dual 0-form Laplacian.
+"""
+function Δᵈ(::Type{Val{0}}, s::SimplicialSets.HasDeltaSet)
+  dd0 = dec_dual_derivative(0, s);
+  ihs1 = dec_inv_hodge_star(1, s, GeometricHodge());
+  d1 = dec_differential(1,s);
+  hs2 = dec_hodge_star(2, s, GeometricHodge());
+  m = hs2 * d1
+  x -> hs2 * d1 * ihs1(dd0 * x)
+end
+
+"""    function Δᵈ_mat(::Type{Val{2}}, s::SimplicialSets.HasDeltaSet)
+
+Return a function matrix encoding the dual 1-form Laplacian.
+"""
+function Δᵈ(::Type{Val{1}}, s::SimplicialSets.HasDeltaSet)
+  dd0 = dec_dual_derivative(0, s);
+  ihs1 = dec_inv_hodge_star(1, s, GeometricHodge());
+  d1 = dec_differential(1,s);
+  hs2 = dec_hodge_star(2, s, GeometricHodge());
+  dd1 = dec_dual_derivative(1, s);
+  ihs0 = dec_inv_hodge_star(0, s, GeometricHodge());
+  d0 = dec_differential(0,s);
+  hs1 = dec_hodge_star(1, s, GeometricHodge());
+  m = hs1 * d0 * ihs0 * dd1
+  n = dd0 * hs2 * d1
+  x -> begin
+    m * x +
+    n * ihs1(x)
+  end
+end
+
 end
diff --git a/src/Meshes.jl b/src/Meshes.jl
index a1ec90ce..a76333ed 100644
--- a/src/Meshes.jl
+++ b/src/Meshes.jl
@@ -73,10 +73,15 @@ Load in a mesh from the Artifacts.jl system and automatically subdivide it.
 """
 loadmesh(s, subdivision=Circumcenter()) = subdivide_duals!(loadmesh(s), subdivision)
 
-# Once part of GridMeshes
 
-# Note that dx will be slightly less than what is given, since points are
-# compressed to lie within max_x.
+# This function was once the gridmeshes.jl file from Decapodes.jl.
+"""    function triangulated_grid(max_x, max_y, dx, dy, point_type)
+
+Return a grid of (near) equilateral triangles.
+
+Note that dx will be slightly less than what is given, since points are
+compressed to lie within max_x.
+"""
 function triangulated_grid(max_x, max_y, dx, dy, point_type)
 
   s = EmbeddedDeltaSet2D{Bool, point_type}()
@@ -125,8 +130,7 @@ function triangulated_grid(max_x, max_y, dx, dy, point_type)
   s
 end
 
-# Once part of SphericalMeshes
-
+# This function was once the sphericalmeshes.jl file from Decapodes.jl.
 """
     makeSphere(minLat, maxLat, dLat, minLong, maxLong, dLong, radius)
 
@@ -134,6 +138,10 @@ Construct a spherical mesh (inclusively) bounded by the given latitudes and
 longitudes, discretized at dLat and dLong intervals, at the given radius from
 Earth's center.
 
+Note that this construction returns a UV-sphere. DEC simulations are more
+accurate on meshes with (near) equilateral triangles, such as the icospheres
+available through [`loadmesh`](@ref).
+
 We say that:
 - 90°N is 0
 - 90°S is 180
@@ -284,4 +292,75 @@ function makeSphere(minLat, maxLat, dLat, minLong, maxLong, dLong, radius)
   return s, north_pole_idx, south_pole_idx
 end
 
-end
\ No newline at end of file
+"""    function tri_345()
+
+Return the primal and dual mesh of a triangle with edge lengths 3,4,5 and with true orientation.
+See also: [`tri_345_false`](@ref)
+"""
+function tri_345()
+  primal_s = EmbeddedDeltaSet2D{Bool,Point3D}()
+  add_vertices!(primal_s, 3, point=[Point3D(0,0,0), Point3D(3,0,0), Point3D(3,4,0)])
+  glue_triangle!(primal_s, 1, 2, 3, tri_orientation=true)
+  primal_s[:edge_orientation] = true
+  s = EmbeddedDeltaDualComplex2D{Bool,Float64,Point3D}(primal_s)
+  subdivide_duals!(s, Barycenter())
+  (primal_s, s)
+end
+
+"""    function tri_345_false()
+
+Return the primal and dual mesh of a triangle with edge lengths 3,4,5 and with false orientation.
+See also: [`tri_345`](@ref)
+"""
+function tri_345_false()
+  primal_s = EmbeddedDeltaSet2D{Bool,Point3D}()
+  add_vertices!(primal_s, 3, point=[Point3D(0,0,0), Point3D(3,0,0), Point3D(3,4,0)])
+  glue_triangle!(primal_s, 1, 2, 3, tri_orientation=false)
+  primal_s[:edge_orientation] = true
+  s = EmbeddedDeltaDualComplex2D{Bool,Float64,Point3D}(primal_s)
+  subdivide_duals!(s, Barycenter())
+  (primal_s, s)
+end
+
+"""    function grid_345()
+
+Return the primal and dual mesh of a grid of 3-4-5 triangles.
+See also: [`tri_345`](@ref)
+"""
+function grid_345()
+  primal_s = EmbeddedDeltaSet2D{Bool,Point3D}()
+  add_vertices!(primal_s, 9,
+    point=[Point3D(0,+4,0), Point3D(3,+4,0), Point3D(6,+4,0),
+          Point3D(0, 0,0), Point3D(3, 0,0), Point3D(6, 0,0),
+          Point3D(0,-4,0), Point3D(3,-4,0), Point3D(6,-4,0)])
+  glue_sorted_triangle!(primal_s, 1, 2, 4)
+  glue_sorted_triangle!(primal_s, 5, 2, 4)
+  glue_sorted_triangle!(primal_s, 5, 2, 3)
+  glue_sorted_triangle!(primal_s, 5, 6, 3)
+  glue_sorted_triangle!(primal_s, 5, 7, 4)
+  glue_sorted_triangle!(primal_s, 5, 7, 8)
+  glue_sorted_triangle!(primal_s, 5, 6, 8)
+  glue_sorted_triangle!(primal_s, 9, 6, 8)
+  primal_s[:edge_orientation] = true
+  s = EmbeddedDeltaDualComplex2D{Bool,Float64,Point3D}(primal_s)
+  subdivide_duals!(s, Barycenter())
+  (primal_s, s)
+end
+
+"""    function right_scalene_unit_hypot()
+
+Return the primal and dual mesh of a right scalene triangle with unit hypotenuse.
+"""
+function right_scalene_unit_hypot()
+  primal_s = EmbeddedDeltaSet2D{Bool,Point2D}()
+  add_vertices!(primal_s, 3,
+    point=[Point2D(0,0), Point2D(1/√2,0), Point2D(1/√2,1/√2)])
+  glue_sorted_triangle!(primal_s, 1, 2, 3)
+  primal_s[:edge_orientation] = true
+  s = EmbeddedDeltaDualComplex2D{Bool,Float64,Point2D}(primal_s)
+  subdivide_duals!(s, Barycenter())
+  (primal_s, s)
+end
+
+
+end
diff --git a/src/SimplicialSets.jl b/src/SimplicialSets.jl
index 5211acd3..d406723e 100644
--- a/src/SimplicialSets.jl
+++ b/src/SimplicialSets.jl
@@ -34,7 +34,8 @@ export Simplex, V, E, Tri, Tet, SimplexChain, VChain, EChain, TriChain, TetChain
   add_triangle!, glue_triangle!, glue_sorted_triangle!,
   tetrahedron_triangles, tetrahedron_edges, tetrahedron_vertices, ntetrahedra,
   tetrahedra, add_tetrahedron!, glue_tetrahedron!, glue_sorted_tetrahedron!,
-  is_manifold_like, nonboundaries, glue_sorted_tet_cube!
+  glue_sorted_tet_cube!, is_manifold_like, nonboundaries,
+  star, St, closed_star, St̄, link, Lk
 
 using LinearAlgebra: det
 using SparseArrays
@@ -905,4 +906,93 @@ function nonboundaries(s::HasDeltaSet, ::Type{Simplex{n}}) where n
   end
 end
 
+# Topological helper functions
+##############################
+
+star(s::AbstractDeltaSet1D, v::Int) = star(s, v, E)
+star(s::AbstractDeltaSet2D, v::Int) = star(s, v, Tri)
+star(s::AbstractDeltaSet3D, v::Int) = star(s, v, Tet)
+
+""" Star of a vertex in a delta set.
+
+Munkres §2 ≈ "The union of the interiors of those simplices of s that have v as
+a vertex."
+
+Return a vector of simplices of dimensions 0 to n.
+
+Recall that interior(σ) = σ - boundary(σ), Munkres §1.
+
+Note that we are returning interiors alone. This means, e.g. a triangle may be
+returned without one or more of its edges. Consequentially, the output of this
+function may not be storable in an ACSet.
+
+This is not the Hodge star [`⋆`](@ref).
+
+See also [`closed_star`](@ref), [`link`](@ref).
+"""
+function star(s::HasDeltaSet, v::Int, ::Type{Simplex{n}}) where n
+  # Recursively compute cofaces, incrementing dimension.
+  cofaces_1n = accumulate(1:n; init=[v]) do c, p
+    Simplex{p}(union([Iterators.flatten(coface(p,i,s,c)) for i in 0:p]...))
+  end
+  pushfirst!(cofaces_1n, V([v]))
+end
+
+""" Alias for the star operator [`star`](@ref), not the Hodge star.
+"""
+St = star
+
+closed_star(s::AbstractDeltaSet1D, v::Int) = closed_star(s, v, star(s, v), E)
+closed_star(s::AbstractDeltaSet2D, v::Int) = closed_star(s, v, star(s, v), Tri)
+closed_star(s::AbstractDeltaSet3D, v::Int) = closed_star(s, v, star(s, v), Tet)
+
+""" Closed star of a vertex in a delta set.
+
+Munkres §2 ≈ "The union of all simplices of s having v as a vertex."
+
+Return a vector of simplices of dimensions 0 to n.
+
+Note that we do not return polytopes, but rather the simplices which together
+form such polytopes, in no particular order.
+
+This is not the Hodge star [`⋆`](@ref).
+
+See also [`star`](@ref), [`link`](@ref).
+"""
+function closed_star(s::HasDeltaSet, v::Int, Sts::Vector{Simplex}, ::Type{Simplex{n}}) where n
+  faces_0nminus1 = map(1:n, Sts, Sts[begin+1:end]) do p, cₚ, cₚ₊₁
+    Simplex{p-1}(union(cₚ, [∂(p,i,s,cₚ₊₁) for i in 0:p]...))
+  end
+  push!(faces_0nminus1, last(Sts))
+end
+
+""" Alias for the closed star operator [`closed_star`](@ref), not the Hodge star.
+"""
+St̄ = closed_star
+
+link(s::AbstractDeltaSet1D, v::Int) = link(s, v, E)
+link(s::AbstractDeltaSet2D, v::Int) = link(s, v, Tri)
+link(s::AbstractDeltaSet3D, v::Int) = link(s, v, Tet)
+
+""" Link of a vertex in a delta set.
+
+Munkres §2 ≈ "The set St̄(v) - St(v)."
+
+Return a vector of simplices of dimensions 0 to n.
+
+These are the simplices which are in the closed star of v, but not in the star
+of v.
+
+See also [`star`](@ref), [`closed_star`](@ref).
+"""
+function link(s::HasDeltaSet, v::Int, ::Type{Simplex{n}}) where n
+  map(closed_star(s,v), star(s,v)) do closed, interior
+    setdiff(closed, interior)
+  end
+end
+
+""" Alias for the link operator [`link`](@ref).
+"""
+Lk = link
+
 end
diff --git a/test/DiscreteExteriorCalculus.jl b/test/DiscreteExteriorCalculus.jl
index fc494e41..09857a7b 100644
--- a/test/DiscreteExteriorCalculus.jl
+++ b/test/DiscreteExteriorCalculus.jl
@@ -1,11 +1,16 @@
 module TestDiscreteExteriorCalculus
 using Test
 
-using LinearAlgebra: Diagonal
+using LinearAlgebra: Diagonal, mul!, norm, dot
 using SparseArrays, StaticArrays
 
 using Catlab.CategoricalAlgebra.CSets
+using ACSets
+using ACSets.DenseACSets: attrtype_type
 using CombinatorialSpaces
+using CombinatorialSpaces.Meshes: tri_345, tri_345_false, grid_345, right_scalene_unit_hypot
+using CombinatorialSpaces.DiscreteExteriorCalculus: eval_constant_primal_form, eval_constant_dual_form
+using GeometryBasics: Point2, Point3
 
 const Point2D = SVector{2,Float64}
 const Point3D = SVector{3,Float64}
@@ -29,6 +34,10 @@ dual_es = elementary_duals(0,s,5)
 @test s[dual_es, :D_∂v0] == edge_center(s, 1:4)
 @test elementary_duals(s, V(5)) == DualE(dual_es)
 
+primal_s′ = subdivide(primal_s)
+@test nv(primal_s′) == nv(primal_s) + ne(primal_s)
+@test ne(primal_s′) == 2*ne(primal_s)
+
 # 1D oriented dual complex
 #-------------------------
 
@@ -44,6 +53,11 @@ s = OrientedDeltaDualComplex1D{Bool}(primal_s)
 @test dual_boundary(1,s) == ∂(1,s)'
 @test dual_derivative(0,s) == -d(0,s)'
 
+primal_s′ = subdivide(primal_s)
+@test nv(primal_s′) == nv(primal_s) + ne(primal_s)
+@test ne(primal_s′) == 2*ne(primal_s)
+@test orient!(primal_s′)
+
 # 1D embedded dual complex
 #-------------------------
 
@@ -163,6 +177,56 @@ end
 # 2D embedded dual complex
 #-------------------------
 
+unit_vector(θ) = cos(θ), sin(θ), 0
+
+function get_regular_polygon(n::Int)
+  n < 3 && error("Cannot construct a polygon with fewer than 3 points.")
+  primal_s = EmbeddedDeltaSet2D{Bool,Point3D}()
+  exterior_points = map(Point3D∘unit_vector, range(0, 2pi-(pi/(n/2)), length=n))
+  add_vertices!(primal_s, n+1, point=[exterior_points..., Point3D(0, 0, 0)])
+  foreach(1:n-1) do x
+    glue_sorted_triangle!(primal_s, n+1, x, x+1)
+  end
+  glue_sorted_triangle!(primal_s, n+1, n, 1)
+  primal_s
+end
+
+# Triangulated hexagon in ℝ³,
+# from Hirani §5.5 Figure 5.5, showing ♭ᵈᵖᵖ(X) to be 0 for a non-trivial X.
+primal_s = get_regular_polygon(6)
+# Rotate counter-clockwise by pi/6 to match the Hirani figure.
+θ = -pi/6
+primal_s[:point] = [[[cos(θ), -sin(θ), 0];; [sin(θ), cos(θ), 0];; [0,0,1]] * p for p in primal_s[:point]]
+s = EmbeddedDeltaDualComplex2D{Bool,Float64,Point3D}(primal_s)
+subdivide_duals!(s, Circumcenter())
+X = map([8,4,0,8,4,0]) do i
+  SVector(unit_vector(i*(2pi/12)))
+end
+@test all(♭(s, DualVectorField(X))) do i
+  isapprox(i, 0.0; atol=1e-15)
+end
+♭_m = ♭_mat(s)
+@test all(reinterpret(Float64, ♭_m * X)) do i
+  isapprox(i, 0.0; atol=1e-15)
+end
+# Half of Proposition 5.5.3: ♭ᵈᵖᵖ is not injective:
+@test all(map(♭_m * X, ♭_m * 2X) do x♭, x2♭
+  isapprox(x♭, x2♭; atol=1e-15)
+end)
+
+#using Random
+#Random.seed!(1)
+#Y = map(1:6) do i
+#  10 .* SVector(randn(), randn(), randn())
+#end
+Y = SVector{3, Float64}[[-0.7058313895389791, 5.314767537831963, -8.06852326006714], [24.56991333983293, 11.648740735275195, 2.6756441862888507], [17.499336925282453, -8.260207919192975, -10.427524178910968], [-3.291338458564041, -4.822519154091156, 11.822427034001585], [4.806914145449797, -0.025581805543075972, 14.346125517139171], [5.320444346188966, 2.4384867352385236, 0.20245933893191853]]
+Y♭ = zeros(SVector{1, Float64}, ne(s))
+mul!(Y♭, ♭_m, Y)
+Y♭_floats = reinterpret(Float64, Y♭)
+@test all(map(♭(s, DualVectorField(Y)), Y♭_floats) do orig, new
+  isapprox(orig, new; atol=20*eps(Float64))
+end)
+
 # Single triangle: numerical example from Gillette's notes on DEC, §2.13.
 #
 # Compared with Gillette, edges #2 and #3 are swapped in the ordering, which
@@ -259,6 +323,73 @@ subdivide_duals!(s, Incenter())
                         -1.657  1.172 -1.172;
                         -1.172 -1.172 -1.657], atol=1e-3)
 
+# This plots a primal vector field over a simplicial complex.
+# This is helpful for debugging operators like sharp and flat.
+#using CairoMakie
+#""" Plot the primal vector field X♯ over the simplicial complex primal_s.
+#"""
+#function plot_pvf(primal_s, X♯; ls=1f0, title="Primal Vector Field")
+#  # Makie will throw an error if the colorrange end points are not unique:
+#  extX = extrema(norm.(X♯))
+#  range = extX[1] ≈ extX[2] ? (0,extX[2]) : extX
+#  f = Figure()
+#  ax = Axis(f[1, 1], title=title)
+#  wireframe!(ax, primal_s, color=:gray95)
+#  scatter!(ax, getindex.(primal_s[:point],1), getindex.(primal_s[:point],2), color = norm.(X♯), colorrange=range)
+#  arrows!(ax, getindex.(primal_s[:point],1), getindex.(primal_s[:point],2), getindex.(X♯,1), getindex.(X♯,2), lengthscale=ls)
+#  Colorbar(f[1,2], limits=range)
+#  hidedecorations!(ax)
+#  f
+#end
+#plot_pvf(primal_s, ♯(s, E))
+#plot_pvf(primal_s, ♯(s, E, AltPPSharp()))
+
+# Evaluate a constant 1-form α assuming Euclidean space. (Inner product is dot)
+eval_constant_form(s, α::SVector) = map(edges(s)) do e
+  dot(α, point(s, tgt(s,e)) - point(s, src(s,e))) * sign(1,s,e)
+end |> EForm
+
+function test_♯(s, covector::SVector; atol=1e-8)
+  X = eval_constant_form(s, covector)
+  # Test that the Hirani field is approximately parallel to the given field.
+  X♯ = ♯(s, X)
+  @test all(map(X♯) do x
+    isapprox(dot(x, covector), norm(x) * norm(covector), atol=atol)
+  end)
+  # Test that the Alternate field approximately equals the given field.
+  X♯ = ♯(s, X, AltPPSharp())
+  @test all(isapprox.(X♯, [covector]))
+  # Test that the matrix and non-matrix versions yield the same result.
+  @test all(isapprox.(♯_mat(s, PPSharp()) * X, ♯(s, X, PPSharp())))
+  @test all(isapprox.(♯_mat(s, AltPPSharp()) * X, ♯(s, X, AltPPSharp())))
+end
+#using Random
+#Random.seed!(1)
+#Point3D(randn(), randn(), 0)
+vfs = [Point3D(1,0,0), Point3D(1,1,0), Point3D(-3,2,0), Point3D(0,0,0),
+  Point3D(-0.07058313895389791, 0.5314767537831963, 0.0)]
+
+# 3,4,5 triangle.
+primal_s, s = tri_345()
+foreach(vf -> test_♯(s, vf), vfs)
+♯_m = ♯_mat(s, DesbrunSharp())
+X = eval_constant_form(s, Point3D(1,0,0))
+X♯ = ♯_m * X
+@test all(X♯ .≈ [Point3D(.8,.4,0), Point3D(0,-1,0), Point3D(-.8,.6,0)])
+
+# Grid of 3,4,5 triangles.
+primal_s, s = grid_345()
+foreach(vf -> test_♯(s, vf), vfs)
+# TODO: Compute results for Desbrun's ♯ by hand.
+
+# Triangulated regular dodecagon.
+primal_s = get_regular_polygon(12)
+primal_s[:point] = [Point3D(1/4,1/5,0) + p for p in primal_s[:point]]
+s = EmbeddedDeltaDualComplex2D{Bool,Float64,Point3D}(primal_s)
+subdivide_duals!(s, Circumcenter())
+foreach(vf -> test_♯(s, vf), vfs)
+# TODO: Compute results for Desbrun's ♯ by hand.
+
 # Triangulated square with consistent orientation.
 primal_s = EmbeddedDeltaSet2D{Bool,Point2D}()
 add_vertices!(primal_s, 4, point=[Point2D(-1,+1), Point2D(+1,+1),
@@ -268,18 +399,23 @@ glue_triangle!(primal_s, 1, 3, 4, tri_orientation=true)
 primal_s[:edge_orientation] = true
 s = EmbeddedDeltaDualComplex2D{Bool,Float64,Point2D}(primal_s)
 subdivide_duals!(s, Barycenter())
+♭_m = ♭_mat(s)
 
-x̂, ŷ, zero = @SVector([1,0]), @SVector([0,1]), @SVector([0,0])
+x̂, ŷ, ẑ = @SVector([1,0]), @SVector([0,1]), @SVector([0,0])
 @test ♭(s, DualVectorField([x̂, -x̂])) ≈ EForm([2,0,0,2,0])
 @test ♭(s, DualVectorField([ŷ, -ŷ])) ≈ EForm([0,-2,0,0,2])
 @test ♭(s, DualVectorField([(x̂-ŷ)/√2, (x̂-ŷ)/√2]))[3] ≈ 2*√2
-@test ♭(s, DualVectorField([(x̂-ŷ)/√2, zero]))[3] ≈ √2
+@test ♭(s, DualVectorField([(x̂-ŷ)/√2, ẑ]))[3] ≈ √2
+@test reinterpret(Float64, ♭_m * [x̂, -x̂]) ≈ EForm([2,0,0,2,0])
+@test reinterpret(Float64, ♭_m * [ŷ, -ŷ]) ≈ EForm([0,-2,0,0,2])
+@test reinterpret(Float64, ♭_m * [(x̂-ŷ)/√2, (x̂-ŷ)/√2])[3] ≈ 2*√2
+@test reinterpret(Float64, ♭_m * [(x̂-ŷ)/√2, ẑ])[3] ≈ √2
 X = ♯(s, EForm([2,0,0,2,0]))::PrimalVectorField
 @test X[2][1] > 0 && X[4][1] < 0
 X = ♯(s, EForm([0,-2,0,0,2]))
 @test X[2][2] > 0 && X[4][2] < 0
 
-@test ∧(s, VForm([2,2,2,2]), TriForm([2.5, 5]))::TriForm ≈ TriForm([2.5, 5])
+@test all(∧(s, VForm([2,2,2,2]), TriForm([2.5, 5]))::TriForm .≈ TriForm([5.0, 10.0]))
 vform, triform = VForm([1.5, 2, 2.5, 3]), TriForm([5, 7.5])
 @test ∧(s, vform, triform) ≈ ∧(s, triform, vform)
 eform1, eform2 = EForm([1.5, 2, 2.5, 3, 3.5]), EForm([3, 7, 10, 11, 15])
@@ -325,4 +461,200 @@ subdivide_duals!(s, Barycenter())
 @test isapprox(Δ(2, s), reshape([-24.0], (1,1)), atol=1e-3)
 @test isapprox(Δ(2, s; hodge=DiagonalHodge()), reshape([-24.0], (1,1)), atol=1e-3)
 
+# 3,4,5 triangle of unit area.
+# Unit area makes computing integrals by hand simple.
+primal_s = EmbeddedDeltaSet2D{Bool,Point2D}()
+add_vertices!(primal_s, 3, point=[Point2D(0,0), Point2D(6/8,0), Point2D(6/8,8/3)])
+glue_triangle!(primal_s, 1, 2, 3, tri_orientation=true)
+primal_s[:edge_orientation] = true
+s = EmbeddedDeltaDualComplex2D{Bool,Float64,Point2D}(primal_s)
+subdivide_duals!(s, Barycenter())
+#@assert only(s[:area]) == 1.0
+
+# The wedge product of 0-form and a 2-form "is multiplication."
+# α = 3
+# β = 5 dx ∧ dy
+# α ∧ β = 15 dx ∧ dy
+# We also write just αβ, when α is a 0-form.
+@test only(∧(s, VForm([3,3,3]), TriForm([only(s[:area])*5]))) ≈ 15
+
+# Grid of 3,4,5 triangles.
+primal_s, s = grid_345()
+
+# ∀ βᵏ ∧(α,βᵏ) = id(βᵏ), where α = 1.
+α = VForm(ones(nv(s)))
+for k = 0:2
+  βᵏ = SimplexForm{k}(collect(1:nsimplices(k,s)))
+  @test all(∧(s, α, βᵏ) .≈ βᵏ)
+end
+
+# 1dx ∧ 1dy = 1 dx∧dy
+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 
+  a * (o ? -1 : 1)
+end
+
+# 1dx∧dy = -1dy∧dx
+@test ∧(s, onedx, onedy) == -∧(s, onedy, onedx)
+
+# 3dx ∧ 2dy = 2 dx ∧ 3dy
+@test ∧(s, EForm(3*onedx), EForm(2*onedy)) == ∧(s, EForm(2*onedx), EForm(3*onedy))
+
+# A triangulated quadrilateral where edges are all of distinct length.
+primal_s = EmbeddedDeltaSet2D{Bool,Point2D}()
+add_vertices!(primal_s, 4, point=[Point2D(0,0), Point2D(1,0), Point2D(0,2), Point2D(-2,5)])
+glue_triangle!(primal_s, 1, 2, 3)
+glue_triangle!(primal_s, 1, 3, 4)
+s = EmbeddedDeltaDualComplex2D{Bool,Float64,Point2D}(primal_s)
+subdivide_duals!(s, Barycenter())
+X = [SVector(2,3), SVector(5,7)]
+♭_m = ♭_mat(s)
+X♭ = zeros(ne(s))
+mul!(X♭, ♭_m,  DualVectorField(X))
+@test all(map(♭(s, DualVectorField(X)), X♭) do orig, new
+  isapprox(orig, new; atol=20*eps(Float64))
+end)
+@test all(map(♭(s, DualVectorField(X)), ♭_m * DualVectorField(X)) do orig, new
+  isapprox(orig, new; atol=20*eps(Float64))
+end)
+
+# Test whether the implementation of ♭ assumes that the DualVectorField vector
+# is indexed in the same order as Tri i.e. 1:ntriangles(s).
+primal_s = EmbeddedDeltaSet2D{Bool,Point2D}()
+add_vertices!(primal_s, 4, point=[Point2D(0,0), Point2D(1,0), Point2D(0,2), Point2D(-2,5)])
+glue_triangle!(primal_s, 1, 2, 3)
+glue_triangle!(primal_s, 1, 3, 4)
+s = EmbeddedDeltaDualComplex2D{Bool,Float64,Point2D}(primal_s)
+subdivide_duals!(s, Barycenter())
+
+primal_s′ = EmbeddedDeltaSet2D{Bool,Point2D}()
+add_vertices!(primal_s′, 4, point=[Point2D(0,0), Point2D(1,0), Point2D(0,2), Point2D(-2,5)])
+glue_triangle!(primal_s′, 1, 2, 3)
+glue_triangle!(primal_s′, 1, 3, 4)
+s′ = EmbeddedDeltaDualComplex2D{Bool,Float64,Point2D}(primal_s′)
+s′[1, :tri_center] = 11
+s′[2, :tri_center] = 10
+s′[[11,13,15,17,19,21], :D_∂v0] = 11
+s′[[12,14,16,18,20,22], :D_∂v0] = 10
+subdivide_duals!(s′, Barycenter())
+#@assert is_isomorphic(s,s′)
+
+X = [SVector(2,3), SVector(5,7)]
+
+@test ♭(s, DualVectorField(X)) == ♭(s′, DualVectorField(X))
+@test ♭_mat(s) * DualVectorField(X) == ♭_mat(s′) * DualVectorField(X)
+
+tg′ = triangulated_grid(100,100,10,10,Point2D);
+tg = EmbeddedDeltaDualComplex2D{Bool,Float64,Point2D}(tg′);
+subdivide_duals!(tg, Barycenter());
+
+rect′ = loadmesh(Rectangle_30x10());
+rect = EmbeddedDeltaDualComplex2D{Bool,Float64,Point3D}(rect′);
+subdivide_duals!(rect, Barycenter());
+
+flat_meshes = [tri_345(), tri_345_false(), right_scalene_unit_hypot(), grid_345(), (tg′, tg), (rect′, rect)];
+
+# Test that the technique for evaluating 1-forms is consistent across primal
+# and dual forms.
+-# ♭ᵈ_discrete(f_continuous) = ⋆₁_discrete∘♭ᵖ_discrete(⋆_continuous f_continuous)
+for (primal_s,s) in flat_meshes
+  α = SVector(1/√2,1/√2,0)
+  left_hand_α = SVector(1/√2,-1/√2,0)
+  f′ = eval_constant_primal_form(s, left_hand_α)
+  f̃ = eval_constant_dual_form(s, α)
+  @test all(map((x,y) -> isapprox( x,y,atol=1e-15), hodge_star(1,s) * f′, f̃))
+end
+
+function all_are_approx(f::Vector{SVector{3,Float64}}, g::SVector{3,Float64}; atol)
+  all(map(f) do f̂
+    all(isapprox.(f̂, g, atol=atol))
+  end)
+end
+function all_are_approx(f::Vector{SVector{2,Float64}}, g::SVector{3,Float64}; atol)
+  all(map(f) do f̂
+        all(isapprox.(f̂, SVector{2,Float64}(g[1],g[2]), atol=atol))
+  end)
+end
+
+for (primal_s,s) in flat_meshes
+  # Test that the least-squares ♯ from dual vector fields to dual 1-forms
+  # preserves constant fields.
+  ♯_m = ♯_mat(s, LLSDDSharp())
+
+  # This tests that numerically raising indices is equivalent to analytically
+  # raising indices.
+  # X_continuous = ♯ᵖ_discrete∘♭ᵖ_discrete(X_continuous)
+  X♯ = SVector{3,Float64}(1/√2,1/√2,0)
+  X♭ = eval_constant_dual_form(s, X♯)
+  @test all_are_approx(♯_m * X♭,  X♯; atol=1e-13) ||
+        all_are_approx(♯_m * X♭, -X♯; atol=1e-13)
+
+  # This tests that numerically raising-then-lowering indices is equivalent to
+  # evaluating a 1-form directly.
+  # Demonstrate how to cache the chained musical isomorphisms.
+  ♭_m = ♭_mat(s)
+  ♭_♯_m = ♭_m * ♯_m
+  # Note that we check explicitly both cases of signedness, because orient!
+  # picks in/outside without regard to the embedding.
+  # This is the "right/left-hand-rule trick."
+  @test all(isapprox.(only.(♭_♯_m * X♭),  eval_constant_primal_form(s, X♯), atol=1e-12)) ||
+        all(isapprox.(only.(♭_♯_m * X♭), -eval_constant_primal_form(s, X♯), atol=1e-12))
+
+  # This test shows how the musical isomorphism chaining lets you further
+  # define a primal-dual wedge that preserves properties from the continuous
+  # exterior calculus.
+  Λpd = dec_wedge_product_pd(Tuple{1,1}, s)
+  Λdp = dec_wedge_product_dp(Tuple{1,1}, s)
+
+  f_def = SVector{3,Float64}(2,7,0)
+  g_def = SVector{3,Float64}(8,1,0)
+  f′ = eval_constant_primal_form(s, f_def)
+  g̃  = eval_constant_dual_form(s, g_def)
+  f̃  = eval_constant_dual_form(s, f_def)
+  g′ = eval_constant_primal_form(s, g_def)
+
+  # Antisymmetry:
+  @test all(Λpd(f′, g̃) .≈ -1 * Λpd(g′, f̃))
+  @test all(Λdp(f̃, g′) .≈ -1 * Λdp(g̃, f′))
+
+  # Antisymmetry (across Λpd and Λdp!):
+  @test all(Λpd(f′, g̃) .== -1 * Λdp(g̃, f′))
+  @test all(Λdp(f̃, g′) .== -1 * Λpd(g′, f̃))
+
+  # f∧f = 0 (implied by antisymmetry):
+  @test all(isapprox.(Λpd(f′, f̃), 0, atol=1e-10))
+  @test all(isapprox.(Λpd(g′, g̃), 0, atol=1e-10))
+
+  # Test and demonstrate the convenience functions:
+  @test all(∧(s, SimplexForm{1}(f′), DualForm{1}(g̃)) .≈ -1 * ∧(s, SimplexForm{1}(g′), DualForm{1}(f̃)))
+  @test all(∧(s, DualForm{1}(f̃), SimplexForm{1}(g′)) .≈ -1 * ∧(s, DualForm{1}(g̃), SimplexForm{1}(f′)))
+
+  # TODO: Test the Leibniz rule.
+end
+
+# Discretely test the continuous property:
+# u :=  fdx + gdy
+# ⋆u = -gdx + fdy
+# u∧⋆u = (fdx + gdy)∧(-gdx + fdy)
+#      = ffdxdy - ggdydx
+#      = (ff+gg)dxdy
+# ⋆(u∧⋆u) = ff+gg
+for (primal_s,s) in flat_meshes
+  f = 2
+  g = 7
+  ff_gg = f*f + g*g
+
+  u_def = SVector{3,Float64}(f,g,0)
+
+  u = eval_constant_primal_form(s, u_def)
+  u_star = hodge_star(1,s) * u
+
+  @test all(isapprox.(
+    sign(2,s) .* hodge_star(2,s) * ∧(s, SimplexForm{1}(u), DualForm{1}(u_star)),
+    ff_gg,
+    atol=1e-10))
+end
+
 end
diff --git a/test/Operators.jl b/test/Operators.jl
index a6f6add5..52016fca 100644
--- a/test/Operators.jl
+++ b/test/Operators.jl
@@ -3,9 +3,14 @@ module TestOperators
 using Test
 using SparseArrays
 using LinearAlgebra
+using Catlab
 using CombinatorialSpaces
+using CombinatorialSpaces.Meshes: tri_345, tri_345_false, grid_345, right_scalene_unit_hypot
+using CombinatorialSpaces.DiscreteExteriorCalculus: eval_constant_primal_form
 using Random
 using GeometryBasics: Point2, Point3
+using StaticArrays: SVector
+using Statistics: mean, var
 
 Point2D = Point2{Float64}
 Point3D = Point3{Float64}
@@ -50,6 +55,16 @@ dual_meshes_2D = [(generate_dual_mesh ∘ loadmesh ∘ Icosphere).(1:2)...,
                (generate_dual_mesh).([triangulated_grid(10,10,8,8,Point3D), makeSphere(5, 175, 5, 0, 360, 5, 6371+90)[1]])...,
                (loadmesh)(Torus_30x10())];
 
+tg′ = triangulated_grid(100,100,10,10,Point2D);
+tg = EmbeddedDeltaDualComplex2D{Bool,Float64,Point2D}(tg′);
+subdivide_duals!(tg, Barycenter());
+
+rect′ = loadmesh(Rectangle_30x10());
+rect = EmbeddedDeltaDualComplex2D{Bool,Float64,Point3D}(rect′);
+subdivide_duals!(rect, Barycenter());
+
+flat_meshes = [tri_345()[2], tri_345_false()[2], right_scalene_unit_hypot()[2], grid_345()[2], tg, rect];
+
 @testset "Exterior Derivative" begin
     for i in 0:0 
         for sd in dual_meshes_1D
@@ -186,10 +201,170 @@ end
         @test all(isapprox.(wdg01(V_1, E_2), ∧(Tuple{0, 1}, sd, V_1, E_2); rtol = 1e-14))
         @test all(wdg01(V_ones, E_ones) .== E_ones)
 
-        @test all(dec_wedge_product(Tuple{0, 2}, sd)(V_1, T_2) / 2 .== ∧(Tuple{0, 2}, sd, V_1, T_2))
+        @test all(dec_wedge_product(Tuple{0, 2}, sd)(V_1, T_2) .== ∧(Tuple{0, 2}, sd, V_1, T_2))
 
         @test all(dec_wedge_product(Tuple{1, 1}, sd)(E_1, E_2) .== ∧(Tuple{1, 1}, sd, E_1, E_2))
     end
 end
 
-end
\ No newline at end of file
+# TODO: Move all boundary helper functions into CombinatorialSpaces.
+function boundary_inds(::Type{Val{0}}, s)
+  ∂1_inds = boundary_inds(Val{1}, s)
+  unique(vcat(s[∂1_inds,:∂v0],s[∂1_inds,:∂v1]))
+end
+function boundary_inds(::Type{Val{1}}, s)
+  collect(findall(x -> x != 0, boundary(Val{2},s) * fill(1,ntriangles(s))))
+end
+function boundary_inds(::Type{Val{2}}, s)
+  ∂1_inds = boundary_inds(Val{1}, s)
+  inds = map([:∂e0, :∂e1, :∂e2]) do esym
+    vcat(incident(s, ∂1_inds, esym)...)
+  end
+  unique(vcat(inds...))
+end
+
+@testset "Primal-Dual Wedge Product 0-1" begin
+    for sd in flat_meshes
+      # Allocate the cached wedge operator.
+      Λ10 = dec_wedge_product_dp(Tuple{1,0}, sd)
+      Λ01 = dec_wedge_product_pd(Tuple{0,1}, sd)
+      ♯_m = ♯_mat(sd, LLSDDSharp())
+
+      # Define test data
+      X♯ = SVector{3,Float64}(1/√2,1/√2,0)
+      f = hodge_star(1,sd) * eval_constant_primal_form(sd, X♯)
+      g = fill(5.0, nv(sd))
+
+      # f := 1/√2dx + 1/√2dy:
+      # g := 5
+      # ⋆f = -1/√2dx + 1/√2dy
+      # ⋆f∧g = 5(-1/√2dx + 1/√2dy) = -5/√2dx + 5/√2dy
+      @test all(Λ10(f,g) .≈ hodge_star(1,sd) * eval_constant_primal_form(sd, SVector{3,Float64}(5/√2,5/√2,0)))
+
+      # Test symmetry across Λ10 and Λ01.
+      @test all(Λ10(f,g) .== Λ01(g,f))
+    end
+    for sd in flat_meshes
+      # Here, We test only for matching values on interior edges, because the
+      # numerical solution assumes values past the boundary are 0, whereas the
+      # analytic solution has no such artifacting.  i.e. The numerical solution
+      # does not hold on boundary edges.
+      interior_edges = setdiff(edges(sd), boundary_inds(Val{1}, sd))
+      length(interior_edges) == 0 && continue
+
+      # Allocate the cached wedge operator.
+      Λ10 = dec_wedge_product_dp(Tuple{1,0}, sd)
+      ♯_m = ♯_mat(sd, LLSDDSharp())
+      # Define test data and the analytic solution.
+      a = map(point(sd)) do p
+        p[1] + 4*p[2]
+      end
+      f = hodge_star(1,sd) * d(0,sd) * a
+      g = map(point(sd)) do p
+        4*p[1] + 16*p[2]
+      end
+      h = map(point(sd)) do p
+        -p[1] - 4*p[2]
+      end
+      dx = eval_constant_primal_form(sd, SVector{3,Float64}(1,0,0))
+      dy = eval_constant_primal_form(sd, SVector{3,Float64}(0,1,0))
+      fΛa_analytic = hodge_star(1,sd) * (-dec_wedge_product(Tuple{0,1}, sd)(h, dx) .+ dec_wedge_product(Tuple{0,1}, sd)(g, dy))
+
+      # a := x + 4y
+      # f := ⋆da
+      # f = ⋆(∂a/∂x dx + ∂a/∂y dy)
+      #   = ⋆(dx + 4dy)
+      # f = 4dx - dy
+      # f∧a = (4dx - dy) ∧ (x + 4y)
+      #     = 4(x + 4y)dx -(x + 4y)dy
+      #     = (4x + 16y)dx + (-x - 4y)dy
+      @test all(isapprox.(Λ10(f,a)[interior_edges], fΛa_analytic[interior_edges], atol=1e-8))
+    end
+end
+
+@testset "Dual-Dual Wedge Product 0-1" begin
+    for sd in flat_meshes
+      # Allocate the cached wedge operator.
+      Λ10 = dec_wedge_product_dd(Tuple{1,0}, sd)
+      Λ01 = dec_wedge_product_dd(Tuple{0,1}, sd)
+
+      # Define test data
+      X♯ = SVector{3,Float64}(1/√2,1/√2,0)
+      f = hodge_star(1,sd) * eval_constant_primal_form(sd, X♯)
+      g = fill(5.0, ntriangles(sd))
+
+      # f := 1/√2dx + 1/√2dy:
+      # g := 5
+      # ⋆f = -1/√2dx + 1/√2dy
+      # ⋆f∧g = 5(-1/√2dx + 1/√2dy) = -5/√2dx + 5/√2dy
+      @test all(Λ10(f,g) .≈ hodge_star(1,sd) * eval_constant_primal_form(sd, SVector{3,Float64}(5/√2,5/√2,0)))
+
+      # Test symmetry across Λ10 and Λ01.
+      @test all(Λ10(f,g) .== Λ01(g,f))
+    end
+end
+
+function euler_equation_test(X♯, sd)
+  # We are not interested in boundary conditions.
+  interior_edges = setdiff(edges(sd), boundary_inds(Val{1}, sd))
+  interior_points = setdiff(vertices(sd), boundary_inds(Val{0}, sd))
+  interior_tris = setdiff(triangles(sd), boundary_inds(Val{2}, sd))
+  # Allocate the cached operators.
+  d0 = dec_dual_derivative(0, sd)
+  ι1 = interior_product_dd(Tuple{1,1}, sd)
+  ι2 = interior_product_dd(Tuple{1,2}, sd)
+  ℒ1 = ℒ_dd(Tuple{1,1}, sd)
+
+  # This is a uniform, constant flow.
+  u = hodge_star(1,sd) * eval_constant_primal_form(sd, X♯)
+
+  # Recall Euler's Equation:
+  # ∂ₜu = ℒᵤu - 0.5dιᵤu = -1/ρdp + b.
+  # We expect for a uniform flow then that ∂ₜu = 0.
+  # (We ignore the pressure term, and only investigate the interior of the
+  # domain.)
+
+  dtu = ℒ1(u,u) - 0.5*d0*ι1(u,u)
+  error = (sign(2,sd) .* ι1(dtu,dtu))[interior_tris]
+  # Note that "error" accumulates in the region around the boundaries.
+  # That is not really error, but rather the effect of boundary conditions and
+  # the influence of pressure.
+  #scatter(sd[sd[:tri_center], :dual_point][interior_tris], color=abs.(error))
+  error
+end
+
+@testset "Dual-Dual Interior Product and Lie Derivative" begin
+  X♯ = SVector{3,Float64}(1/√2,1/√2,0)
+  error = euler_equation_test(X♯, rect)
+  @test abs(mean(error)) < 8e-2
+  @test var(error) < 9e-2
+  # Note that boundary conditions "bleed into" the interior of the domain
+  # somewhat. So this is not all "error" in the usual sense.
+  # TODO: Grab points that are distance X away from the interior.
+
+  X♯ = SVector{3,Float64}(1/√2,1/√2,0)
+  error = euler_equation_test(X♯, tg)
+  @test abs(mean(error)) < 7e-4
+  @test var(error) < 7e-6
+
+  X♯ = SVector{3,Float64}(3,3,0)
+  error = euler_equation_test(X♯, tg)
+  @test abs(mean(error)) < 3e-1
+  @test var(error) < 8e-1
+
+  # u := ⋆xdx
+  # ιᵤu = -x²
+  sd = rect
+  f = map(point(sd)) do p
+    p[1]
+  end
+  dx = eval_constant_primal_form(sd, SVector{3,Float64}(1,0,0))
+  u = hodge_star(1,sd) * dec_wedge_product(Tuple{0,1}, sd)(f, dx)
+  ι1 = interior_product_dd(Tuple{1,1}, sd)
+  interior_tris = setdiff(triangles(sd), boundary_inds(Val{2}, sd))
+  @test all(<(8e-3), (sign(2,sd) .* ι1(u,u) .- map(sd[sd[:tri_center], :dual_point]) do (x,_,_)
+    -x*x
+  end)[interior_tris])
+end
+
+end
diff --git a/test/Project.toml b/test/Project.toml
index 47800dee..7fbc7e0a 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,4 +1,5 @@
 [deps]
+ACSets = "227ef7b5-1206-438b-ac65-934d6da304b8"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Catlab = "134e5e36-593f-5add-ad60-77f754baafbe"
 CombinatorialSpaces = "b1c52339-7909-45ad-8b6a-6e388f7c67f2"
@@ -7,4 +8,5 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/test/SimplicialSets.jl b/test/SimplicialSets.jl
index 64e0897c..3d604bca 100644
--- a/test/SimplicialSets.jl
+++ b/test/SimplicialSets.jl
@@ -417,6 +417,71 @@ for t in tetrahedra(s)
 end
 @test sum([volume(s, Tet(t)) for t in tetrahedra(s)]) == 8*2
 
+# Six tetrahedra of equal volume forming a cube with edge length 1.
+# The connectivity is that of Blessent's 2009 thesis, Figure 3.2.
+s = EmbeddedDeltaSet3D{Bool,Point3D}()
+add_vertices!(s, 8, point=[
+  Point3D(0,1,0), Point3D(0,0,0), Point3D(1,1,0), Point3D(1,0,0),
+  Point3D(0,1,1), Point3D(0,0,1), Point3D(1,1,1), Point3D(1,0,1)])
+# See Table 3.1 "Mesh connectivity".
+glue_sorted_tetrahedron!(s, 3, 5, 4, 2)
+glue_sorted_tetrahedron!(s, 7, 6, 8, 4)
+glue_sorted_tetrahedron!(s, 5, 6, 7, 4)
+glue_sorted_tetrahedron!(s, 3, 5, 7, 4)
+glue_sorted_tetrahedron!(s, 5, 6, 4, 2)
+glue_sorted_tetrahedron!(s, 1, 5, 3, 2)
+@test ntetrahedra(s) == 6
+@test tetrahedra(s) == 1:6
+@test ntriangles(s) == 18
+@test triangles(s) == 1:18 # (2 * num cube faces) + (2 * num "cuts" into cube)
+@test ne(s) == 19 # (num cube edges) + (num cube faces) + (internal diagonal)
+@test is_semi_simplicial(s, 2)
+@test is_semi_simplicial(s, 3)
+s[:edge_orientation] = true
+s[:tri_orientation] = true
+s[:tet_orientation] = true
+orient!(s)
+for t in tetrahedra(s)
+  @test volume(s, Tet(t)) ≈ 1/6
+end
+@test is_manifold_like(s)
+for i in 1:3
+  @test isempty(nonboundaries(s)[i])
+end
+@test d(1, s) * d(0, s) * vertices(s) == zeros(ntriangles(s))
+@test d(2, s) * d(1, s) * edges(s) == zeros(ntetrahedra(s))
+
+# Five tetrahedra example from Letniowski 1992, as given by Blessent Table 3.3b.
+s = EmbeddedDeltaSet3D{Bool,Point3D}()
+add_vertices!(s, 6, point=[
+  # See Table 3.3a "Nodal coordinates"
+  Point3D(-2, -2,   0.5),
+  Point3D( 0, -2,   0.1),
+  Point3D(-2,  0,   0.1),
+  Point3D( 0,  0.1, 0),
+  Point3D(-2, -2,  -0.25),
+  Point3D(-2, -2,   1.5)])
+# See Table 3.3b "Connectivity list for Letniowski's example"
+glue_sorted_tetrahedron!(s, 1, 2, 4, 6)
+glue_sorted_tetrahedron!(s, 1, 3, 4, 6)
+glue_sorted_tetrahedron!(s, 1, 2, 3, 5)
+glue_sorted_tetrahedron!(s, 2, 3, 4, 5)
+glue_sorted_tetrahedron!(s, 1, 2, 3, 4)
+@test ntetrahedra(s) == 5
+@test tetrahedra(s) == 1:5
+@test is_semi_simplicial(s, 2)
+@test is_semi_simplicial(s, 3)
+s[:edge_orientation] = true
+s[:tri_orientation] = true
+s[:tet_orientation] = true
+orient!(s)
+@test is_manifold_like(s)
+for i in 1:3
+  @test isempty(nonboundaries(s)[i])
+end
+@test d(1, s) * d(0, s) * vertices(s) == zeros(ntriangles(s))
+@test d(2, s) * d(1, s) * edges(s) == zeros(ntetrahedra(s))
+
 # Euclidean geometry
 ####################
 
@@ -431,4 +496,61 @@ p1, p2, p3 = Point3D(1,0,0), Point3D(0,1,0), Point3D(0,0,1)
 p1, p2, p3, p4 = SVector(1,0,0,0), SVector(0,1,0,0), SVector(0,0,1,0), SVector(0,0,0,1)
 @test volume([p1, p2, p3, p4]) ≈ std_simplex_volume(3)
 
+# Topological helper functions
+##############################
+
+# §62 Example 1 Figure 62.1 from Munkres 1984:
+s = DeltaSet2D()
+# 6 adjacent triangles forming a hexagon, with a stand-alone edge at the center.
+g, h = 7, 8
+add_vertices!(s, 8)
+foreach(1:5, 2:6) do x,y
+  glue_sorted_triangle!(s, g,x,y)
+end
+glue_sorted_triangle!(s, g,6,1)
+add_sorted_edge!(s, g,h)
+
+Stg = star(s, g)
+@test issetequal(Stg[1], [g])
+@test issetequal(Stg[2], union(coface(1,0,s,g), union(coface(1,1,s,g))))
+@test issetequal(Stg[3], triangles(s))
+
+St̄g = closed_star(s, g)
+@test issetequal(St̄g[1], vertices(s))
+@test issetequal(St̄g[2], edges(s))
+@test issetequal(St̄g[3], triangles(s))
+
+# "The link of the vertex g consists of the hexagon ... and the vertex h."
+Lkg = link(s, g)
+@test issetequal(Lkg[1], [1,2,3,4,5,6, 8])
+@test issetequal(Lkg[2], setdiff(St̄g[2], Stg[2]))
+@test isempty(Lkg[3])
+
+# "The link of the vertex h is the vertex g."
+Lkh = link(s, h)
+@test issetequal(Lkh[1], [g])
+@test isempty(Lkh[2])
+@test isempty(Lkh[3])
+
+# §62 Example 1 Figure 62.2 from Munkres 1984:
+s = DeltaSet3D()
+# 5 adjacent tetrahedra forming a pentagonal bipyramid.
+add_vertices!(s, 7)
+foreach(1:4, 2:5) do x,y
+  glue_sorted_tetrahedron!(s, 6,7, x,y)
+end
+glue_sorted_tetrahedron!(s, 6,7, 5,1)
+
+# "The link of the vertex a is the union of the 2-simplices bfg and efg."
+# bfg=[267], efg=[567]
+@test link(s,1)[3] == [1, 11]
+
+# "The link of the vertex f is the cone: abcdea * g
+Lkf = link(s,6)
+es = union(∂(0, s, Tri(Lkf[3])), ∂(1, s, Tri(Lkf[3])), ∂(2, s, Tri(Lkf[3])))
+vs = union(∂(0, s, E(es)), ∂(1, s, E(es)))
+@test Set(es) == Set(Lkf[2])
+@test Set(vs) == Set(Lkf[1])
+@test Set(Lkf[1]) == Set([1,2,3,4,5,7])
+
 end