From e65b845782693f1a5853cf1d5b35299f39643686 Mon Sep 17 00:00:00 2001
From: Daan Huybrechs <daan.huybrechs@kuleuven.be>
Date: Thu, 10 Jun 2021 16:27:17 +0200
Subject: [PATCH 1/5] improve boundary of cubes

---
 src/DomainSets.jl               |  1 +
 src/applications/coordinates.jl | 31 +++++++++++++++++
 src/domains/cube.jl             | 61 +++++++++++++++++++++++++++++++++
 src/domains/interval.jl         |  1 +
 src/domains/point.jl            | 12 +++----
 src/domains/simplex.jl          |  3 ++
 src/generic/domain.jl           | 24 ++++++++++++-
 src/generic/mapped.jl           | 15 +++++---
 src/generic/setoperations.jl    |  4 ++-
 src/maps/affine.jl              |  3 ++
 src/util/common.jl              |  1 +
 test/test_maps.jl               |  5 +++
 test/test_specific_domains.jl   | 10 +++++-
 13 files changed, 158 insertions(+), 13 deletions(-)
 create mode 100644 src/applications/coordinates.jl

diff --git a/src/DomainSets.jl b/src/DomainSets.jl
index 72b8a76..09dbb32 100644
--- a/src/DomainSets.jl
+++ b/src/DomainSets.jl
@@ -212,6 +212,7 @@ include("domains/cube.jl")
 include("domains/indicator.jl")
 include("domains/boundingbox.jl")
 
+include("applications/coordinates.jl")
 include("applications/rotation.jl")
 
 end # module
diff --git a/src/applications/coordinates.jl b/src/applications/coordinates.jl
new file mode 100644
index 0000000..eebd2f2
--- /dev/null
+++ b/src/applications/coordinates.jl
@@ -0,0 +1,31 @@
+
+"A `DomainPoint` is a point which is an element of a domain by construction."
+abstract type DomainPoint{T} end
+
+in(x::DomainPoint, d::Domain) = domain(x) == d || in(point(x), d)
+
+convert(::Type{T}, x::DomainPoint) where {T} = convert(T, point(x))
+
+
+## Points on a sphere
+
+"A point on the unit sphere."
+abstract type SpherePoint{T} <: DomainPoint{T} end
+
+domain(p::SpherePoint{T}) where {T<:StaticTypes} = UnitSphere{T}()
+domain(p::SpherePoint{T}) where {T<:AbstractVector} = UnitSphere{eltype(T)}(length(point(x)))
+
+"A point on the unit sphere represented by a standard Euclidean vector."
+struct EuclideanSpherePoint{T} <: SpherePoint{T}
+    x   ::  T
+end
+point(p::EuclideanSpherePoint) = p.x
+
+
+"A point on the unit sphere represented in spherical coordinates."
+struct SphericalCoordinate{T} <: SpherePoint{SVector{3,T}}
+    θ   ::  T   # inclination or polar angle
+    ϕ   ::  T   # azimuthal angle
+end
+
+point(p::SphericalCoordinate) = SVector(sin(p.θ)*cos(p.ϕ), sin(p.θ)*sin(p.ϕ), cos(p.θ))
diff --git a/src/domains/cube.jl b/src/domains/cube.jl
index 2a2c88c..c7149d5 100644
--- a/src/domains/cube.jl
+++ b/src/domains/cube.jl
@@ -2,6 +2,9 @@
 "A `HyperRectangle` is the cartesian product of intervals."
 abstract type HyperRectangle{T} <: ProductDomain{T} end
 
+convert(::Type{HyperRectangle}, d::ProductDomain) =
+	ProductDomain(map(t->convert(AbstractInterval, t), components(d)))
+
 boundingbox(d::HyperRectangle) = d
 
 "Compute all corners of the hyperrectangle."
@@ -19,6 +22,64 @@ function corners(d::HyperRectangle)
     corners
 end
 
+# map the interval [a,b] to a rectangle defined by c (bottom-left) and d (top-right)
+rectangle_map(a::Number, b::Number, c::SVector{2}, d::SVector{2}) = AffineMap((d-c)/(b-a), c - (d-c)/(b-a)*a)
+
+# map the rectangle defined by a and b, to the rectangle defined by c and d,
+# where the dimension of c/d is one larger, and one of the coordinates (dim) is fixed to dimval.
+function rectangle_map(a::SVector{M}, b::SVector{M}, c::SVector{N}, d::SVector{N}, dim, dimval) where {N,M}
+	@assert N == M+1
+	T = promote_type(eltype(a),eltype(c))
+	A = MMatrix{N,M,T}(undef)
+	B = MVector{N,T}(undef)
+	fill!(A, 0)
+	fill!(B, 0)
+	B[dim] = dimval
+	for m = 1:dim-1
+		# scalar map along the dimension "m"
+		mapdim = interval_map(a[m], b[m], c[m], d[m])
+		A[m,m] = mapdim.A
+		B[m] = mapdim.b
+	end
+	for m = dim+1:N
+		mapdim = interval_map(a[m-1], b[m-1], c[m], d[m])
+		A[m,m-1] = mapdim.A
+		B[m] = mapdim.b
+	end
+	AffineMap(SMatrix{N,M}(A), SVector{N}(B))
+end
+
+function boundary(d::HyperRectangle{SVector{2,T}}) where {T}
+	left = infimum(d)
+	right = supremum(d)
+	x1 = left[1]; y1 = left[1]; x2 = right[1]; y2 = right[2]
+	d_unit = UnitInterval{T}()
+	faces = [
+		ParametricDomain(rectangle_map(zero(T), one(T), SVector(x1,y1), SVector(x2,y1)), d_unit),
+		ParametricDomain(rectangle_map(zero(T), one(T), SVector(x2,y1), SVector(x2,y2)), d_unit),
+		ParametricDomain(rectangle_map(zero(T), one(T), SVector(x2,y2), SVector(x1,y2)), d_unit),
+		ParametricDomain(rectangle_map(zero(T), one(T), SVector(x1,y2), SVector(x1,y1)), d_unit)
+	]
+	UnionDomain(faces)
+end
+
+function boundary(d::HyperRectangle{SVector{N,T}}) where {N,T}
+	left2 = infimum(d)
+	right2 = supremum(d)
+	d_unit = UnitCube{SVector{N-1,T}}()
+	left1 = infimum(d_unit)
+	right1 = supremum(d_unit)
+
+	face1 = ParametricDomain(rectangle_map(left1, right1, left2, right2, 1, left2[1]), d_unit)
+	D = typeof(face1)
+	faces = D[]
+	for dim in 1:N
+		push!(faces, ParametricDomain(rectangle_map(left1, right1, left2, right2, dim, left2[dim]), d_unit))
+		push!(faces, ParametricDomain(rectangle_map(left1, right1, left2, right2, dim, right2[dim]), d_unit))
+	end
+	UnionDomain(faces)
+end
+
 
 "A `Cube` is a hyperrectangle with equal side lengths in each dimension."
 abstract type Cube{T} <: HyperRectangle{T} end
diff --git a/src/domains/interval.jl b/src/domains/interval.jl
index 8008b4f..0cfec84 100644
--- a/src/domains/interval.jl
+++ b/src/domains/interval.jl
@@ -31,6 +31,7 @@ isapprox(d1::AbstractInterval, d2::AbstractInterval) =
 
 
 boundary(d::AbstractInterval) = Point(leftendpoint(d)) ∪ Point(rightendpoint(d))
+corners(d::AbstractInterval) = [leftendpoint(d), rightendpoint(d)]
 
 boundingbox(d::AbstractInterval) = d
 
diff --git a/src/domains/point.jl b/src/domains/point.jl
index e71914b..9b5a0fd 100644
--- a/src/domains/point.jl
+++ b/src/domains/point.jl
@@ -37,17 +37,17 @@ isclosedset(d::Point) = true
 boundary(d::Point) = d
 boundingbox(d::Point) = d.x..d.x
 
+infimum(d::Point) = d.x
+supremum(d::Point) = d.x
+
 interior(d::Point{T}) where {T} = EmptySpace{T}()
 closure(d::Point) = d
 
 point_in_domain(d::Point) = d.x
 
-function map_domain(map, p::Point)
-    x = applymap(map, p.x)
-    Point(x)
-end
-
-mapped_domain(invmap, p::Point) = map_domain(inverse(invmap), p)
+mapped_domain(invmap, p::Point) = Point(inverse(invmap, p.x))
+map_domain(map, p::Point) = Point(applymap(map, p.x))
+parametric_domain(map, p::Point) = Point(applymap(map, p.x))
 
 for op in (:+,:-)
     @eval $op(a::Point, b::Point) = Point($op(a.x,b.x))
diff --git a/src/domains/simplex.jl b/src/domains/simplex.jl
index 318ef68..0261e09 100644
--- a/src/domains/simplex.jl
+++ b/src/domains/simplex.jl
@@ -3,6 +3,7 @@
 # An n-dimensional simplex
 ###########################
 
+"Supertype of an N-dimensional simplex."
 abstract type Simplex{T,C} <: Domain{T} end
 
 isclosedset(::Simplex{T,:closed}) where {T} = true
@@ -10,7 +11,9 @@ isclosedset(::Simplex{T,:open}) where {T} = false
 
 isopenset(d::Simplex) = !isclosedset(d)
 
+"The unit simplex is a polytope with the origin and all unit vectors as vertices."
 abstract type UnitSimplex{T,C} <: Simplex{T,C} end
+
 const ClosedUnitSimplex{T} = UnitSimplex{T,:closed}
 const OpenUnitSimplex{T} = UnitSimplex{T,:open}
 
diff --git a/src/generic/domain.jl b/src/generic/domain.jl
index 3c7d6f3..ce26cff 100644
--- a/src/generic/domain.jl
+++ b/src/generic/domain.jl
@@ -145,7 +145,29 @@ isreal(d::Domain) = isreal(eltype(d))
 infimum(d::Domain) = minimum(d)  # if the minimum exists, then it is also the infimum
 supremum(d::Domain) = maximum(d)  # if the maximum exists, then it is also the supremum
 
+
+"""
+Return a bounding box of the given domain.
+
+A bounding box is an interval, a hyperrectangle or the full space. It is such that
+each point in the domain also lies in the bounding box.
+"""
+boundingbox(d) = FullSpace{eltype(d)}()
+
 function boundary end
 const ∂ = boundary
 
-boundingbox(d) = FullSpace{eltype(d)}()
+
+# "Lazy representation of the boundary of a domain."
+# struct Boundary{T,D} <: Domain{T}
+#     domain  :: D
+# end
+#
+# domain(d::Boundary) = d.domain
+#
+# Boundary(domain) = Boundary{eltype(domain)}(domain)
+# Boundary{T}(domain) where {T} = Boundary{T,typeof(domain)}(domain)
+# Boundary{T}(domain::Domain{T}) where {T} = Boundary{T,typeof(domain)}(domain)
+# Boundary{T}(domain::Domain{S}) where {S,T} = Boundary{T}(convert(Domain{T}, domain))
+#
+# indomain(x, d::Boundary) = in(x, boundary(domain(d)))
diff --git a/src/generic/mapped.jl b/src/generic/mapped.jl
index b44fa5d..0597a92 100644
--- a/src/generic/mapped.jl
+++ b/src/generic/mapped.jl
@@ -31,6 +31,7 @@ isempty(d::AbstractMappedDomain) = isempty(superdomain(d))
 isopenset(d::AbstractMappedDomain) = isopenset(superdomain(d))
 isclosedset(d::AbstractMappedDomain) = isclosedset(superdomain(d))
 
+corners(d::AbstractMappedDomain) = [forward_map(d, x) for x in corners(superdomain(d))]
 
 ## I/O functionality
 
@@ -126,8 +127,11 @@ _mapped_domain2(invmap, domain, ::Type{S}, ::Type{T}) where {S,T} =
 # This assumes that the map types can be combined using \circ
 mapped_domain(invmap, d::MappedDomain) = mapped_domain(inverse_map(d) ∘ invmap, superdomain(d))
 
-boundary(d::MappedDomain) = _boundary(d, superdomain(d), inverse_map(d))
-_boundary(d::MappedDomain, superdomain, invmap) = MappedDomain(invmap, boundary(superdomain))
+boundary(d::MappedDomain) = _boundary(d, boundary(superdomain(d)), inverse_map(d))
+_boundary(d::MappedDomain, superbnd, invmap) = MappedDomain(invmap, superbnd)
+_boundary(d::MappedDomain, superbnd::UnionDomain, invmap) =
+    UnionDomain(map(t->mapped_domain(invmap, t), components(superbnd)))
+
 interior(d::MappedDomain) = _interior(d, superdomain(d), inverse_map(d))
 _interior(d::MappedDomain, superdomain, invmap) = MappedDomain(invmap, interior(superdomain))
 closure(d::MappedDomain) = _closure(d, superdomain(d), inverse_map(d))
@@ -161,8 +165,11 @@ inverse_map(d::ParametricDomain, y) = leftinverse(d.fmap, y)
 "Return the domain that results from mapping the given domain."
 parametric_domain(fmap, domain::Domain) = ParametricDomain(fmap, domain)
 
-boundary(d::ParametricDomain) = _boundary(d, superdomain(d), forward_map(d))
-_boundary(d::ParametricDomain, superdomain, fmap) = ParametricDomain(fmap, boundary(superdomain))
+boundary(d::ParametricDomain) = _boundary(d, boundary(superdomain(d)), forward_map(d))
+_boundary(d::ParametricDomain, superbnd, fmap) = ParametricDomain(fmap, superbnd)
+_boundary(d::ParametricDomain, superbnd::UnionDomain, fmap) =
+    UnionDomain(map(t -> parametric_domain(fmap, t), components(superbnd)))
+
 interior(d::ParametricDomain) = _interior(d, superdomain(d), forward_map(d))
 _interior(d::ParametricDomain, superdomain, fmap) = ParametricDomain(fmap, interior(superdomain))
 closure(d::ParametricDomain) = _closure(d, superdomain(d), forward_map(d))
diff --git a/src/generic/setoperations.jl b/src/generic/setoperations.jl
index a889cc7..2f65d61 100644
--- a/src/generic/setoperations.jl
+++ b/src/generic/setoperations.jl
@@ -76,7 +76,7 @@ _ud(d1, d2, d3, domains...) = uniondomain(uniondomain(d1, d2, d3), domains...)
 
 # avoid nested union domains
 uniondomain(d1::UnionDomain, d2::UnionDomain) =
-	d1 == d2 ? d1 : UnionDomain(components(d1)..., components(d2)...)
+	d1 == d2 ? d1 : UnionDomain(collect(Set(components(d1)) ∪ Set(components(d2))))
 uniondomain1(d1::UnionDomain, d2) = UnionDomain(components(d1)..., d2)
 uniondomain2(d1, d2::UnionDomain) = UnionDomain(d1, components(d2)...)
 
@@ -105,6 +105,8 @@ closure(d::UnionDomain) = UnionDomain(map(closure, components(d)))
 
 boundingbox(d::UnionDomain) = unionbox(map(boundingbox, components(d))...)
 
+boundary(d::UnionDomain) = uniondomain(map(boundary, components(d))...)
+
 Display.combinationsymbol(d::UnionDomain) = Display.Symbol('∪')
 Display.displaystencil(d::UnionDomain) = composite_displaystencil(d)
 show(io::IO, mime::MIME"text/plain", d::UnionDomain) = Display.composite_show(io, mime, d)
diff --git a/src/maps/affine.jl b/src/maps/affine.jl
index 1fe8058..f550633 100644
--- a/src/maps/affine.jl
+++ b/src/maps/affine.jl
@@ -370,6 +370,9 @@ AffineMap{T}(A, b) where {T} = GenericAffineMap{T}(A, b)
 
 similarmap(m::AffineMap, ::Type{T}) where {T} = AffineMap{T}(m.A, m.b)
 
+convert(::Type{AffineMap}, m) = (@assert isaffine(m); AffineMap(matrix(m), vector(m)))
+convert(::Type{AffineMap{T}}, m) where {T} = (@assert isaffine(m); AffineMap{T}(matrix(m), vector(m)))
+
 # If y = A*x+b, then x = inv(A)*(y-b) = inv(A)*y - inv(A)*b
 inverse(m::AffineMap) = (@assert issquaremap(m); AffineMap(inv(m.A), -inv(m.A)*m.b))
 inverse(m::AffineMap, x) = (@assert issquaremap(m); m.A \ (x-m.b))
diff --git a/src/util/common.jl b/src/util/common.jl
index 411d11d..200dd36 100644
--- a/src/util/common.jl
+++ b/src/util/common.jl
@@ -127,3 +127,4 @@ end
 # we use matrix_pinv rather than pinv to preserve static matrices
 matrix_pinv(A) = pinv(A)
 matrix_pinv(A::SMatrix{M,N}) where {M,N} = SMatrix{N,M}(pinv(A))
+matrix_pinv(A::SVector{N}) where {N} = convert(Transpose{Float64, SVector{N,Float64}}, pinv(A))
diff --git a/test/test_maps.jl b/test/test_maps.jl
index 2369c48..c0b4582 100644
--- a/test/test_maps.jl
+++ b/test/test_maps.jl
@@ -287,6 +287,11 @@ function test_affine_maps(T)
     @test DomainSets.to_vector(Vector{T}, A) == [0,0]
     @test DomainSets.to_vector(T, 2, 3) == 3
 
+    if T != BigFloat    # BigFloat's make pinv fail for StaticArrays
+        @test DomainSets.matrix_pinv(SMatrix{2,2}(rand(T),rand(T),rand(T),rand(T))) isa SMatrix{2,2}
+        @test DomainSets.matrix_pinv(SVector(rand(T),rand(T))) isa Transpose{T,SVector{2,T}}
+    end
+
     test_linearmap(T)
     test_translation(T)
     test_affinemap(T)
diff --git a/test/test_specific_domains.jl b/test/test_specific_domains.jl
index 4dd2abd..e8e0b88 100644
--- a/test/test_specific_domains.jl
+++ b/test/test_specific_domains.jl
@@ -129,6 +129,8 @@ end
         @test !isempty(d)
         @test boundary(d) == d
         @test boundingbox(d) == 1.0..1.0
+        @test infimum(d) == d.x
+        @test supremum(d) == d.x
         @test isclosedset(d)
         @test !isopenset(d)
         @test dimension(d) == 1
@@ -196,6 +198,8 @@ end
             @test leftendpoint(d) ∈ ∂(d)
             @test rightendpoint(d) ∈ ∂(d)
 
+            @test boundary(d) == uniondomain(Point(zero(T)), Point(one(T)))
+            @test corners(d) == [0,1]
             @test boundingbox(d) == d
 
             @test similar_interval(0..1, 0, big(1.0)) isa ClosedInterval{BigFloat}
@@ -1126,7 +1130,7 @@ end
         @test inverse_map(pd, m(0.4)) ≈ 0.4
         @test mapfrom_canonical(pd) == m
         @test canonicaldomain(pd) == 0..1
-        @test boundary(pd) isa ParametricDomain
+        @test boundary(pd) isa UnionDomain
         @test interior(pd) isa ParametricDomain
         @test closure(pd) isa ParametricDomain
     end
@@ -1538,6 +1542,9 @@ end
         @test UnitCube{SVector{2,BigFloat}}(Val(2)) isa EuclideanUnitCube{2,BigFloat}
 
         @test Set(corners((0..1)^2)) == Set([ [1,0], [0,1], [1,1], [0,0]])
+        @test Set(corners(UnitCube())) == Set([ [0,0,0], [1,0,0], [0,1,0], [0,0,1], [1,1,0], [1,0,1], [0,1,1], [1,1,1]])
+        @test boundary(boundary(UnitSquare())) == UnionDomain(map(t->Point(SVector{2}(t)), corners(UnitSquare())))
+        @test boundary(boundary(boundary(UnitCube()))) == UnionDomain(map(t->Point(SVector{3}(t)), corners(UnitCube())))
 
         @test StaticUnitCube() == UnitCube()
         @test StaticUnitCube(Val(3)) == UnitCube()
@@ -1579,6 +1586,7 @@ end
         @test !isempty(D)
         @test isclosedset(D)
         @test !isopenset(D)
+        @test convert(DomainSets.HyperRectangle, VcatDomain(UnitInterval(), UnitInterval())) isa StaticUnitCube
 
         @test approx_in(SA[-0.1,-0.1], D, 0.1)
         @test !approx_in(SA[-0.1,-0.1], D, 0.09)

From 7f49fbf371b2ddac4d1ddaf3447aa9e6e1759a90 Mon Sep 17 00:00:00 2001
From: Daan Huybrechs <daan.huybrechs@kuleuven.be>
Date: Fri, 11 Jun 2021 20:53:11 +0200
Subject: [PATCH 2/5] implement hash of domains, improve boundary of cubes

---
 src/DomainSets.jl             |  1 +
 src/domains/ball.jl           | 12 +++++
 src/domains/cube.jl           | 87 +++++++++++++++++++++++++++++------
 src/domains/indicator.jl      |  2 +
 src/domains/interval.jl       |  7 +++
 src/domains/levelset.jl       | 10 +++-
 src/domains/point.jl          |  6 ++-
 src/domains/simplex.jl        |  2 +
 src/domains/trivial.jl        |  6 +++
 src/generic/canonical.jl      | 17 ++++---
 src/generic/domain.jl         | 23 +++++++--
 src/generic/mapped.jl         | 16 ++++++-
 src/generic/productdomain.jl  | 10 ++--
 src/generic/setoperations.jl  |  4 +-
 src/maps/affine.jl            |  5 +-
 src/maps/basic.jl             |  4 +-
 src/maps/composite.jl         |  1 +
 src/maps/product.jl           |  1 +
 src/util/common.jl            |  4 ++
 test/test_specific_domains.jl | 10 +++-
 20 files changed, 191 insertions(+), 37 deletions(-)

diff --git a/src/DomainSets.jl b/src/DomainSets.jl
index 09dbb32..c26344e 100644
--- a/src/DomainSets.jl
+++ b/src/DomainSets.jl
@@ -101,6 +101,7 @@ export Domain, EuclideanDomain, VectorDomain,
     interior, closure,
     volume,
     point_in_domain,
+    normal, tangents,
     canonicaldomain, mapto_canonical, mapfrom_canonical, hascanonicaldomain,
     mapto,
     parameterdomain, parameterization, hasparameterization,
diff --git a/src/domains/ball.jl b/src/domains/ball.jl
index 0e1c205..84e1b41 100644
--- a/src/domains/ball.jl
+++ b/src/domains/ball.jl
@@ -54,6 +54,7 @@ isempty(d::OpenBall) = radius(d) == 0
 
 ==(d1::Ball, d2::Ball) = isclosedset(d1)==isclosedset(d2) &&
     radius(d1)==radius(d2) && center(d1)==center(d2)
+hash(d::Ball, h::UInt) = hashrec("Ball", isclosedset(d), radius(d), center(d), h)
 
 function issubset(d1::Ball, d2::Ball)
     if dimension(d1) == dimension(d2)
@@ -78,6 +79,10 @@ function issubset(d1::Ball, d2::Ball)
     end
 end
 
+normal(d::Ball, x) = normal(boundary(d), x)
+
+distance_to(d::Ball, x) = x ∈ d ? zero(prectype(d)) : norm(x-center(d))-radius(d)
+
 # We choose the center of the ball here. Concrete types should implement 'center'
 point_in_domain(d::Ball) = center(d)
 
@@ -110,6 +115,7 @@ approx_indomain(x, d::OpenUnitBall, tolerance) = norm(x) < 1+tolerance
 approx_indomain(x, d::ClosedUnitBall, tolerance) = norm(x) <= 1+tolerance
 
 
+
 ==(d1::UnitBall, d2::UnitBall) = isclosedset(d1)==isclosedset(d2) &&
     dimension(d1) == dimension(d2)
 
@@ -285,6 +291,10 @@ isempty(::Sphere) = false
 isclosedset(::Sphere) = true
 isopenset(::Sphere) = false
 
+normal(d::Sphere, x) = (x-center(x))/norm(x-center(x))
+
+distance_to(d::Sphere, x) = abs(norm(x-center(d))-radius(d))
+
 point_in_domain(d::Sphere) = center(d) + unitvector(d)
 
 ==(d1::Sphere, d2::Sphere) = radius(d1)==radius(d2) && center(d1)==center(d2)
@@ -302,6 +312,8 @@ abstract type UnitSphere{T} <: Sphere{T} end
 radius(::UnitSphere) = 1
 center(d::UnitSphere{T}) where {T<:StaticTypes} = zero(T)
 
+normal(d::UnitSphere, x) = x/norm(x)
+
 indomain(x, d::UnitSphere) = norm(x) == 1
 approx_indomain(x, d::UnitSphere, tolerance) =
     1-tolerance <= norm(x) <= 1+tolerance
diff --git a/src/domains/cube.jl b/src/domains/cube.jl
index c7149d5..0a322da 100644
--- a/src/domains/cube.jl
+++ b/src/domains/cube.jl
@@ -22,12 +22,16 @@ function corners(d::HyperRectangle)
     corners
 end
 
-# map the interval [a,b] to a rectangle defined by c (bottom-left) and d (top-right)
-rectangle_map(a::Number, b::Number, c::SVector{2}, d::SVector{2}) = AffineMap((d-c)/(b-a), c - (d-c)/(b-a)*a)
+# map the interval [a,b] to a cube defined by c (bottom-left) and d (top-right)
+cube_face_map(a::Number, b::Number, c::SVector{2}, d::SVector{2}) = AffineMap((d-c)/(b-a), c - (d-c)/(b-a)*a)
+function cube_face_map(a::Number, b::Number, c::Vector, d::Vector)
+	@assert length(c) == length(d) == 2
+	AffineMap((d-c)/(b-a), c - (d-c)/(b-a)*a)
+end
 
-# map the rectangle defined by a and b, to the rectangle defined by c and d,
-# where the dimension of c/d is one larger, and one of the coordinates (dim) is fixed to dimval.
-function rectangle_map(a::SVector{M}, b::SVector{M}, c::SVector{N}, d::SVector{N}, dim, dimval) where {N,M}
+# map the cube defined by a and b, to the cube defined by c and d,
+# where the dimension of c and d is one larger, and one of the coordinates (dim) is fixed to dimval.
+function cube_face_map(a::SVector{M}, b::SVector{M}, c::SVector{N}, d::SVector{N}, dim, dimval) where {N,M}
 	@assert N == M+1
 	T = promote_type(eltype(a),eltype(c))
 	A = MMatrix{N,M,T}(undef)
@@ -49,16 +53,42 @@ function rectangle_map(a::SVector{M}, b::SVector{M}, c::SVector{N}, d::SVector{N
 	AffineMap(SMatrix{N,M}(A), SVector{N}(B))
 end
 
+function cube_face_map(a::Vector, b::Vector, c::Vector, d::Vector, dim, dimval)
+	M = length(a)
+	N = length(c)
+	@assert length(a)==length(b)
+	@assert length(c)==length(d)
+	@assert N == M+1
+	T = promote_type(eltype(a),eltype(c))
+	A = Matrix{T}(undef, N, M)
+	B = Vector{T}(undef, N)
+	fill!(A, 0)
+	fill!(B, 0)
+	B[dim] = dimval
+	for m = 1:dim-1
+		# scalar map along the dimension "m"
+		mapdim = interval_map(a[m], b[m], c[m], d[m])
+		A[m,m] = mapdim.A
+		B[m] = mapdim.b
+	end
+	for m = dim+1:N
+		mapdim = interval_map(a[m-1], b[m-1], c[m], d[m])
+		A[m,m-1] = mapdim.A
+		B[m] = mapdim.b
+	end
+	AffineMap(A, B)
+end
+
 function boundary(d::HyperRectangle{SVector{2,T}}) where {T}
 	left = infimum(d)
 	right = supremum(d)
-	x1 = left[1]; y1 = left[1]; x2 = right[1]; y2 = right[2]
+	x1 = left[1]; y1 = left[2]; x2 = right[1]; y2 = right[2]
 	d_unit = UnitInterval{T}()
 	faces = [
-		ParametricDomain(rectangle_map(zero(T), one(T), SVector(x1,y1), SVector(x2,y1)), d_unit),
-		ParametricDomain(rectangle_map(zero(T), one(T), SVector(x2,y1), SVector(x2,y2)), d_unit),
-		ParametricDomain(rectangle_map(zero(T), one(T), SVector(x2,y2), SVector(x1,y2)), d_unit),
-		ParametricDomain(rectangle_map(zero(T), one(T), SVector(x1,y2), SVector(x1,y1)), d_unit)
+		ParametricDomain(cube_face_map(zero(T), one(T), SVector(x1,y1), SVector(x2,y1)), d_unit),
+		ParametricDomain(cube_face_map(zero(T), one(T), SVector(x2,y1), SVector(x2,y2)), d_unit),
+		ParametricDomain(cube_face_map(zero(T), one(T), SVector(x2,y2), SVector(x1,y2)), d_unit),
+		ParametricDomain(cube_face_map(zero(T), one(T), SVector(x1,y2), SVector(x1,y1)), d_unit)
 	]
 	UnionDomain(faces)
 end
@@ -70,16 +100,47 @@ function boundary(d::HyperRectangle{SVector{N,T}}) where {N,T}
 	left1 = infimum(d_unit)
 	right1 = supremum(d_unit)
 
-	face1 = ParametricDomain(rectangle_map(left1, right1, left2, right2, 1, left2[1]), d_unit)
+	face1 = ParametricDomain(cube_face_map(left1, right1, left2, right2, 1, left2[1]), d_unit)
 	D = typeof(face1)
 	faces = D[]
 	for dim in 1:N
-		push!(faces, ParametricDomain(rectangle_map(left1, right1, left2, right2, dim, left2[dim]), d_unit))
-		push!(faces, ParametricDomain(rectangle_map(left1, right1, left2, right2, dim, right2[dim]), d_unit))
+		push!(faces, ParametricDomain(cube_face_map(left1, right1, left2, right2, dim, left2[dim]), d_unit))
+		push!(faces, ParametricDomain(cube_face_map(left1, right1, left2, right2, dim, right2[dim]), d_unit))
 	end
 	UnionDomain(faces)
 end
 
+function boundary(d::HyperRectangle{Vector{T}}) where {T}
+	if dimension(d) == 2
+		left = infimum(d)
+		right = supremum(d)
+		x1 = left[1]; y1 = left[2]; x2 = right[1]; y2 = right[2]
+		d_unit = UnitInterval{T}()
+		faces = [
+			ParametricDomain(cube_face_map(zero(T), one(T), [x1,y1], [x2,y1]), d_unit),
+			ParametricDomain(cube_face_map(zero(T), one(T), [x2,y1], [x2,y2]), d_unit),
+			ParametricDomain(cube_face_map(zero(T), one(T), [x2,y2], [x1,y2]), d_unit),
+			ParametricDomain(cube_face_map(zero(T), one(T), [x1,y2], [x1,y1]), d_unit)
+		]
+		UnionDomain(faces)
+	else
+		left2 = infimum(d)
+		right2 = supremum(d)
+		d_unit = UnitCube(dimension(d)-1)
+		left1 = infimum(d_unit)
+		right1 = supremum(d_unit)
+
+		face1 = ParametricDomain(cube_face_map(left1, right1, left2, right2, 1, left2[1]), d_unit)
+		D = typeof(face1)
+		faces = D[]
+		for dim in 1:dimension(d)
+			push!(faces, ParametricDomain(cube_face_map(left1, right1, left2, right2, dim, left2[dim]), d_unit))
+			push!(faces, ParametricDomain(cube_face_map(left1, right1, left2, right2, dim, right2[dim]), d_unit))
+		end
+		UnionDomain(faces)
+	end
+end
+
 
 "A `Cube` is a hyperrectangle with equal side lengths in each dimension."
 abstract type Cube{T} <: HyperRectangle{T} end
diff --git a/src/domains/indicator.jl b/src/domains/indicator.jl
index 3afd13a..9e257c1 100644
--- a/src/domains/indicator.jl
+++ b/src/domains/indicator.jl
@@ -60,6 +60,8 @@ indomain(x, d::BoundedIndicatorFunction) = in(x, boundingdomain(d)) && d.f(x)
 
 ==(d1::BoundedIndicatorFunction, d2::BoundedIndicatorFunction) =
     indicatorfunction(d1)==indicatorfunction(d2) && boundingdomain(d1)==boundingdomain(d2)
+hash(d::BoundedIndicatorFunction, h::UInt) =
+    hashrec(indicatorfunction(d), boundingdomain(d), h)
 
 similardomain(d::BoundedIndicatorFunction, ::Type{T}) where {T} =
     BoundedIndicatorFunction(d.f, convert(Domain{T}, d.domain))
diff --git a/src/domains/interval.jl b/src/domains/interval.jl
index 0cfec84..8819138 100644
--- a/src/domains/interval.jl
+++ b/src/domains/interval.jl
@@ -5,6 +5,9 @@ iscompact(d::TypedEndpointsInterval) = false
 isinterval(d::Domain) = false
 isinterval(d::AbstractInterval) = true
 
+hash(d::AbstractInterval, h::UInt) =
+    hashrec(isleftopen(d), isrightopen(d), leftendpoint(d), rightendpoint(d), h)
+
 Display.object_parentheses(::AbstractInterval) = true
 
 approx_indomain(x, d::AbstractInterval, tolerance) =
@@ -33,6 +36,10 @@ isapprox(d1::AbstractInterval, d2::AbstractInterval) =
 boundary(d::AbstractInterval) = Point(leftendpoint(d)) ∪ Point(rightendpoint(d))
 corners(d::AbstractInterval) = [leftendpoint(d), rightendpoint(d)]
 
+normal(d::AbstractInterval, x) = (abs(minimum(d)-x) < abs(maximum(d)-x)) ? -one(eltype(d)) : one(eltype(d))
+
+distance_to(d::AbstractInterval, x) = x ∈ d ? zero(eltype(d)) : min(abs(x-supremum(d)), abs(x-infimum(d)))
+
 boundingbox(d::AbstractInterval) = d
 
 volume(d::AbstractInterval) = width(d)
diff --git a/src/domains/levelset.jl b/src/domains/levelset.jl
index 714c50f..50a5ac7 100644
--- a/src/domains/levelset.jl
+++ b/src/domains/levelset.jl
@@ -18,6 +18,7 @@ show(io::IO, d::AbstractLevelSet) =
 
 ==(d1::AbstractLevelSet, d2::AbstractLevelSet) = levelfun(d1)==levelfun(d2) &&
     level(d1)==level(d2)
+hash(d::AbstractLevelSet, h::UInt) = hashrec("AbstractLevelSet", levelfun(d), level(d), h)
 
 "The domain defined by `f(x)=0` for a given function `f`."
 struct ZeroSet{T,F} <: AbstractLevelSet{T}
@@ -60,6 +61,8 @@ show(io::IO, d::AbstractSublevelSet{T,:open}) where {T} =
 
 ==(d1::AbstractSublevelSet, d2::AbstractSublevelSet) = levelfun(d1)==levelfun(d2) &&
     level(d1)==level(d2)
+hash(d::AbstractSublevelSet, h::UInt) =
+    hashrec("AbstractSublevelSet", levelfun(d), level(d), h)
 
 "The domain where `f(x) <= 0` (or `f(x) < 0`)."
 struct SubzeroSet{T,C,F} <: AbstractSublevelSet{T,C}
@@ -110,8 +113,11 @@ show(io::IO, d::AbstractSuperlevelSet{T,:closed}) where {T} =
 show(io::IO, d::AbstractSuperlevelSet{T,:open}) where {T} =
     print(io, "superlevel set f(x) > $(level(d)) with f = $(levelfun(d))")
 
-==(d1::AbstractSuperlevelSet, d2::AbstractSuperlevelSet) = levelfun(d1)==levelfun(d2) &&
-    level(d1)==level(d2)
+==(d1::AbstractSuperlevelSet, d2::AbstractSuperlevelSet) =
+    levelfun(d1)==levelfun(d2) && level(d1)==level(d2)
+hash(d::AbstractSuperlevelSet, h::UInt) =
+    hashrec("AbstractSuperlevelSet", levelfun(d), level(d), h)
+
 
 "The domain where `f(x) >= 0` (or `f(x) > 0`)."
 struct SuperzeroSet{T,C,F} <: AbstractSuperlevelSet{T,C}
diff --git a/src/domains/point.jl b/src/domains/point.jl
index 9b5a0fd..6572cee 100644
--- a/src/domains/point.jl
+++ b/src/domains/point.jl
@@ -17,7 +17,9 @@ Number(d::Point) = convert(Number, d)
 convert(::Type{Domain}, c::Number) = Point(c)
 convert(::Type{Domain{T}}, c::Number) where T = Point{T}(c)
 
-==(a::Point,b::Point) = a.x == b.x
+==(d1::Point,d2::Point) = d1.x == d2.x
+hash(d::Point, h::UInt) = hashrec("Point", d.x, h)
+
 indomain(x, d::Point) = x == d.x
 isempty(::Point) = false
 
@@ -45,6 +47,8 @@ closure(d::Point) = d
 
 point_in_domain(d::Point) = d.x
 
+distance_to(d::Point, x) = norm(x-d.x)
+
 mapped_domain(invmap, p::Point) = Point(inverse(invmap, p.x))
 map_domain(map, p::Point) = Point(applymap(map, p.x))
 parametric_domain(map, p::Point) = Point(applymap(map, p.x))
diff --git a/src/domains/simplex.jl b/src/domains/simplex.jl
index 0261e09..02f1316 100644
--- a/src/domains/simplex.jl
+++ b/src/domains/simplex.jl
@@ -44,9 +44,11 @@ isempty(::UnitSimplex) = false
 
 ==(d1::UnitSimplex, d2::UnitSimplex) =
     isclosedset(d1)==isclosedset(d2) && dimension(d1)==dimension(d2)
+hash(d::UnitSimplex, h::UInt) = hashrec("UnitSimplex", isclosedset(d), dimension(d), h)
 
 boundingbox(d::UnitSimplex{T}) where {T} = UnitCube{T}(dimension(d))
 
+distance_to(d::UnitSimplex, x) = x ∈ d ? zero(prectype(d)) : minimum(distance_to(el, x) for el in components(boundary(d)))
 
 
 struct StaticUnitSimplex{T,C} <: UnitSimplex{T,C}
diff --git a/src/domains/trivial.jl b/src/domains/trivial.jl
index 7005444..a27b1be 100644
--- a/src/domains/trivial.jl
+++ b/src/domains/trivial.jl
@@ -28,6 +28,8 @@ interior(d::EmptySpace) = d
 closure(d::EmptySpace) = d
 boundingbox(d::EmptySpace) = d
 
+distance_to(d::EmptySpace, x) = convert(prectype(d), Inf)
+
 # Arithmetic operations
 
 issubset1(d1::EmptySpace, d2) = true
@@ -41,6 +43,7 @@ mapped_domain(map::Map, d::EmptySpace) = EmptySpace{codomaintype(map)}()
 ==(d1::EmptySpace, d2::EmptySpace) = true
 isequal1(d1::EmptySpace, d2) = isempty(d2)
 isequal2(d1, d2::EmptySpace) = isempty(d1)
+hash(d::EmptySpace, h::UInt) = hash("EmptySpace", h)
 
 
 "The full space of elements of type `T`."
@@ -79,6 +82,8 @@ boundary(d::FullSpace{T}) where {T} = EmptySpace{T}()
 interior(d::FullSpace) = d
 closure(d::FullSpace) = d
 
+distance_to(d::FullSpace, x) = zero(prectype(d))
+
 # Arithmetic operations
 
 issubset1(d1::FullSpace, d2) = isfullspace(d2)
@@ -89,6 +94,7 @@ map_domain(m::AbstractAffineMap{T}, d::FullSpace{T}) where {T} = d
 ==(d1::FullSpace, d2::FullSpace) = true
 isequal1(d1::FullSpace, d2) = isfullspace(d2)
 isequal2(d1, d2::FullSpace) = isfullspace(d1)
+hash(d::FullSpace, h::UInt) = hash("FullSpace", h)
 
 
 convert(::Type{Domain}, ::Type{T}) where T = FullSpace{T}()
diff --git a/src/generic/canonical.jl b/src/generic/canonical.jl
index 58d6f5a..36b68ee 100644
--- a/src/generic/canonical.jl
+++ b/src/generic/canonical.jl
@@ -12,12 +12,15 @@ For example, the canonical domain of an Interval `[a,b]` is the interval `[-1,1]
 """
 canonicaldomain(d::Domain) = d
 
+"""
+Check whether the given arguments are different. They are different when
+they have a different type, or if they have the same type but they are not equal.
+"""
+isdifferentfrom(d1::D, d2::D) where D = !(d1==d2)
+isdifferentfrom(d1, d2) = true
+
 "Does the domain have a canonical domain?"
-hascanonicaldomain(d) = _hascanonicaldomain(d, canonicaldomain(d))
-# the domain and the canonical domain have the same type: check equality
-_hascanonicaldomain(d::D, cd::D) where {D} = !(d==cd)
-# they have a different type: result is true even if the domains are equal
-_hascanonicaldomain(d, cd) = true
+hascanonicaldomain(d) = isdifferentfrom(d, canonicaldomain(d))
 
 identitymap(d) = IdentityMap{eltype(d)}(dimension(d))
 
@@ -34,7 +37,7 @@ mapto_canonical(d, x) = mapto_canonical(d)(x)
 abstract type CanonicalType end
 
 canonicaldomain(ctype::CanonicalType, d) = d
-hascanonicaldomain(ctype::CanonicalType, d) = _hascanonicaldomain(d, canonicaldomain(ctype, d))
+hascanonicaldomain(ctype::CanonicalType, d) = isdifferentfrom(d, canonicaldomain(ctype, d))
 
 mapfrom_canonical(ctype::CanonicalType, d, x) = mapfrom_canonical(ctype, d)(x)
 mapto_canonical(ctype::CanonicalType, d, x) = mapto_canonical(ctype, d)(x)
@@ -48,7 +51,7 @@ mapfrom_canonical(::Equal, d) = identitymap(d)
 mapto_canonical(::Equal, d) = leftinverse(mapto_canonical(Equal(), d))
 
 simplify(d) = canonicaldomain(Equal(), d)
-simplifies(d) = !(simplify(d)===d)
+simplifies(d) = hascanonicaldomain(Equal(), d)
 
 "A canonical domain that is isomorphic but may have different element type."
 struct Isomorphic <: CanonicalType end
diff --git a/src/generic/domain.jl b/src/generic/domain.jl
index ce26cff..a588579 100644
--- a/src/generic/domain.jl
+++ b/src/generic/domain.jl
@@ -16,6 +16,13 @@ Domain(d) = convert(Domain, d)
 convert(::Type{Domain{T}}, d::Domain{T}) where {T} = d
 convert(::Type{Domain{T}}, d::Domain{S}) where {S,T} = similardomain(d, T)
 
+"Can the domains be promoted without throwing an error?"
+promotable_domains(d1, d2) = promotable_eltypes(eltype(d1),eltype(d2))
+promotable_eltypes(::Type{S}, ::Type{T}) where {S,T} =
+    isconcretetype(promote_type(S, T))
+promotable_eltypes(::Type{S}, ::Type{T}) where {S<:AbstractVector,T<:AbstractVector} =
+    promotable_eltypes(eltype(S), eltype(T))
+
 "Promote the given domains to have a common element type."
 promote_domains() = ()
 promote_domains(domains...) = promote_domains(domains)
@@ -33,9 +40,6 @@ _convert_eltype(::Type{T}, d, ::Type{S}) where {S,T} =
 convert_eltype(::Type{T}, d::AbstractArray) where {T} = convert(AbstractArray{T}, d)
 convert_eltype(::Type{T}, d::Set) where {T} = convert(Set{T}, d)
 
-# TODO: rename and clarify this function
-compatible_eltype(d1, d2) = isconcretetype(promote_type(eltype(d1),eltype(d2)))
-
 promote(d1::Domain, d2::Domain) = promote_domains((d1,d2))
 promote(d1::Domain, d2) = promote_domains((d1,d2))
 promote(d1, d2::Domain) = promote_domains((d1,d2))
@@ -154,9 +158,22 @@ each point in the domain also lies in the bounding box.
 """
 boundingbox(d) = FullSpace{eltype(d)}()
 
+"Return the boundary of the given domain as a domain."
 function boundary end
 const ∂ = boundary
 
+"""
+Return the normal of the domain at the point `x`.
+
+It is assumed that `x` is a point on the boundary of the domain.
+"""
+function normal end
+
+"""
+Return the tangents of the domain at the point `x`. The tangents form a
+basis for the tangent plane, perpendicular to the normal direction at `x`.
+"""
+function tangents end
 
 # "Lazy representation of the boundary of a domain."
 # struct Boundary{T,D} <: Domain{T}
diff --git a/src/generic/mapped.jl b/src/generic/mapped.jl
index 0597a92..456d10f 100644
--- a/src/generic/mapped.jl
+++ b/src/generic/mapped.jl
@@ -159,11 +159,23 @@ similardomain(d::ParametricDomain, ::Type{T}) where {T} =
 forward_map(d::ParametricDomain) = d.fmap
 forward_map(d::ParametricDomain, x) = d.fmap(x)
 
-inverse_map(d::ParametricDomain) = leftinverse(d.fmap)
-inverse_map(d::ParametricDomain, y) = leftinverse(d.fmap, y)
+function indomain(x, d::ParametricDomain)
+    # To check for membership, we can't use the inverse map because it may not exist
+    # We assume a left inverse exists, but the left inverse may be many-to-one.
+    # So we also have to check whether the left-inverse-point maps back to x
+    y = leftinverse(d.fmap, x)
+    x2 = forward_map(d, y)
+    isapprox(x, x2)
+end
+
+==(d1::ParametricDomain, d2::ParametricDomain) =
+    forward_map(d1) == forward_map(d2) && superdomain(d1) == superdomain(d2)
+hash(d::ParametricDomain, h::UInt) = hashrec(forward_map(d), superdomain(d), h)
 
 "Return the domain that results from mapping the given domain."
 parametric_domain(fmap, domain::Domain) = ParametricDomain(fmap, domain)
+parametric_domain(fmap, domain::ParametricDomain) =
+    parametric_domain(fmap ∘ forward_map(domain), superdomain(domain))
 
 boundary(d::ParametricDomain) = _boundary(d, boundary(superdomain(d)), forward_map(d))
 _boundary(d::ParametricDomain, superbnd, fmap) = ParametricDomain(fmap, superbnd)
diff --git a/src/generic/productdomain.jl b/src/generic/productdomain.jl
index e2bbe2b..beb88b3 100644
--- a/src/generic/productdomain.jl
+++ b/src/generic/productdomain.jl
@@ -7,6 +7,7 @@ composition(d::ProductDomain) = Product()
 components(d::ProductDomain) = d.domains
 
 ==(d1::ProductDomain, d2::ProductDomain) = mapreduce(==, &, components(d1), components(d2))
+hash(d::ProductDomain, h::UInt) = hashrec("ProductDomain", collect(components(d)), h)
 
 isempty(d::ProductDomain) = any(isempty, components(d))
 isclosedset(d::ProductDomain) = all(isclosedset, components(d))
@@ -17,6 +18,8 @@ issubset(d1::ProductDomain, d2::ProductDomain) =
 
 volume(d::ProductDomain) = prod(map(volume, components(d)))
 
+distance_to(d::ProductDomain, x) = sqrt(sum(distance_to(component(d, i), x[i])^2 for i in 1:ncomponents(d)))
+
 compatibleproductdims(d1::ProductDomain, d2::ProductDomain) =
 	dimension(d1) == dimension(d2) &&
 		all(map(==, map(dimension, components(d1)), map(dimension, components(d2))))
@@ -29,10 +32,11 @@ show(io::IO, d::ProductDomain) = composite_show_compact(io, d)
 boundary_part(d::ProductDomain{T}, domains, i) where {T} =
 	ProductDomain{T}(domains[1:i-1]..., boundary(domains[i]), domains[i+1:end]...)
 
-boundary(d::ProductDomain) = _boundary(d, components(d))
-_boundary(d::ProductDomain, domains) =
+boundary(d::ProductDomain) = productboundary(d)
+productboundary(d) = productboundary(d, components(d))
+productboundary(d, domains) =
 	UnionDomain(boundary_part(d, domains, i) for i in 1:length(domains))
-_boundary(d::ProductDomain, domains::Tuple) =
+productboundary(d, domains::Tuple) =
 	UnionDomain(tuple((boundary_part(d, domains, i) for i in 1:length(domains))...))
 
 boundingbox(d::ProductDomain{T}) where {T} = ProductDomain{T}(map(boundingbox, components(d)))
diff --git a/src/generic/setoperations.jl b/src/generic/setoperations.jl
index 2f65d61..d6af007 100644
--- a/src/generic/setoperations.jl
+++ b/src/generic/setoperations.jl
@@ -1,6 +1,6 @@
 # The union, intersection and difference of domains are represented with lazy domains.
 
-issubset(d1::Domain, d2::Domain) = compatible_eltype(d1, d2) && issubset1(promote_domains(d1, d2)...)
+issubset(d1::Domain, d2::Domain) = promotable_domains(d1, d2) && issubset1(promote_domains(d1, d2)...)
 issubset(d1::Domain, d2) = issubset1(promote_domains(d1, d2)...)
 issubset(d1, d2::Domain) = issubset1(promote_domains(d1, d2)...)
 issubset1(d1, d2) = issubset2(d1, d2)
@@ -81,6 +81,7 @@ uniondomain1(d1::UnionDomain, d2) = UnionDomain(components(d1)..., d2)
 uniondomain2(d1, d2::UnionDomain) = UnionDomain(d1, components(d2)...)
 
 ==(a::UnionDomain, b::UnionDomain) = Set(components(a)) == Set(components(b))
+hash(d::UnionDomain, h::UInt) = hashrec("UnionDomain", Set(d.domains), h)
 
 
 convert(::Type{Domain}, v::AbstractVector{<:Domain}) = UnionDomain(v)
@@ -93,7 +94,6 @@ convert(::Type{Domain{T}}, s::AbstractSet) where {T} = UnionDomain{T}(map(Point,
 similardomain(d::UnionDomain, ::Type{T}) where {T} =
     UnionDomain(convert.(Domain{T}, components(d)))
 
-hash(d::UnionDomain, h::UInt) = hash(Set(d.domains), h)
 
 
 point_in_domain(d::UnionDomain) = convert(eltype(d), point_in_domain(component(d,1)))
diff --git a/src/maps/affine.jl b/src/maps/affine.jl
index f550633..388a0a0 100644
--- a/src/maps/affine.jl
+++ b/src/maps/affine.jl
@@ -56,6 +56,8 @@ isaffine(m::AbstractAffineMap) = true
     islinear(m1) && matrix(m1) == matrix(m2)
 ==(m1::IdentityMap, m2::AbstractAffineMap) = m2==m1
 
+hash(m::AbstractAffineMap, h::UInt) = hashrec(h, matrix(m), vector(m))
+
 mapsize(m::AbstractAffineMap) = _affine_mapsize(m, domaintype(m), unsafe_matrix(m), unsafe_vector(m))
 _affine_mapsize(m::AbstractAffineMap, T, A::AbstractArray, b) = size(A)
 _affine_mapsize(m::AbstractAffineMap, T, A::AbstractVector, b::AbstractVector) = (length(A),)
@@ -125,6 +127,7 @@ isreal(m::LinearMap) = _isreal(m, unsafe_matrix(m))
 _isreal(m::LinearMap, A) = isreal(A)
 
 ==(m1::LinearMap, m2::LinearMap) = matrix(m1) == matrix(m2)
+hash(m::LinearMap, h::UInt) = hash(matrix(m), h)
 
 # inverse should be called only on square maps, otherwise use
 # leftinverse or rightinverse in order to use pinv instead of inv
@@ -311,7 +314,7 @@ jacdet(m::Translation{T}, x) where {T} = one(eltype(T))
 isreal(m::Translation) = isreal(unsafe_vector(m))
 
 ==(m1::Translation, m2::Translation) = unsafe_vector(m1)==unsafe_vector(m2)
-
+hash(m::Translation, h::UInt) = hash(unsafe_vector(m), h)
 
 similarmap(m::Translation, ::Type{T}) where {T} = Translation{T}(m.b)
 
diff --git a/src/maps/basic.jl b/src/maps/basic.jl
index 6bfa1a4..d5859e1 100644
--- a/src/maps/basic.jl
+++ b/src/maps/basic.jl
@@ -62,7 +62,7 @@ similarmap(m::StaticIdentityMap, ::Type{T}) where {T} =
 convert(::Type{StaticIdentityMap{T}}, ::StaticIdentityMap) where {T} = StaticIdentityMap{T}()
 
 ==(m1::StaticIdentityMap, m2::StaticIdentityMap) = true
-
+hash(m::StaticIdentityMap, h::UInt) = hash("StaticIdentityMap", h)
 
 "Identity map with dynamic size determined by a dimension field."
 struct DynamicIdentityMap{T} <: IdentityMap{T}
@@ -83,6 +83,7 @@ similarmap(m::DynamicIdentityMap, ::Type{T}) where {T<:StaticTypes} =
     StaticIdentityMap{T}()
 
 ==(m1::DynamicIdentityMap, m2::DynamicIdentityMap) = m1.dimension == m2.dimension
+hash(m::DynamicIdentityMap, h::UInt) = hashrec("DynamicIdentityMap", dimension(m), h)
 
 
 "The supertype of constant maps from `T` to `U`."
@@ -114,6 +115,7 @@ determinantmap(m::ConstantMap{T}) where {T} = ConstantMap{T}(det(constant(m)))
 absmap(m::ConstantMap{T}) where {T} = ConstantMap{T}(abs(constant(m)))
 
 ==(m1::ConstantMap, m2::ConstantMap) = constant(m1)==constant(m2)
+hash(m::ConstantMap, h::UInt) = hashrec("ConstantMap", constant(m), h)
 
 similarmap(m::ConstantMap, ::Type{T}) where {T} = ConstantMap{T}(constant(m))
 similarmap(m::ConstantMap, ::Type{T}, ::Type{U}) where {T,U} = ConstantMap{T,U}(m.c)
diff --git a/src/maps/composite.jl b/src/maps/composite.jl
index 36f6ef2..ff30755 100644
--- a/src/maps/composite.jl
+++ b/src/maps/composite.jl
@@ -77,6 +77,7 @@ composedmap2(m1, m2::ComposedMap) = ComposedMap(m1, components(m2)...)
 
 ==(m1::ComposedMap, m2::ComposedMap) =
     ncomponents(m1) == ncomponents(m2) && all(map(isequal, components(m1), components(m2)))
+hash(m::ComposedMap, h::UInt) = hashrec("ComposedMap", ncomponents(m), components(m), h)
 
 Display.combinationsymbol(m::ComposedMap) = Display.Symbol('∘')
 Display.displaystencil(m::ComposedMap) =
diff --git a/src/maps/product.jl b/src/maps/product.jl
index 9960a18..c0bb227 100644
--- a/src/maps/product.jl
+++ b/src/maps/product.jl
@@ -75,6 +75,7 @@ end
 mapsize(m::ProductMap) = (sum(t->mapsize(t,1), components(m)), sum(t->mapsize(t,2), components(m)))
 
 ==(m1::ProductMap, m2::ProductMap) = all(map(isequal, components(m1), components(m2)))
+hash(m::ProductMap, h::UInt) = hashrec("ProductMap", components(m), h)
 
 Display.combinationsymbol(m::ProductMap) = Display.Symbol('⊗')
 Display.displaystencil(m::ProductMap) = composite_displaystencil(m)
diff --git a/src/util/common.jl b/src/util/common.jl
index 200dd36..1159444 100644
--- a/src/util/common.jl
+++ b/src/util/common.jl
@@ -20,6 +20,10 @@ function unitvector(d::Domain{T}) where {T<:AbstractVector}
 end
 unitvector(d::Domain{T}) where {T<:Number} = one(T)
 
+"Apply the `hash` function recursively to the given arguments."
+hashrec(h) = hash(h)
+hashrec(h, args...) = hash(h, hashrec(args...))
+
 #################
 # Precision type
 #################
diff --git a/test/test_specific_domains.jl b/test/test_specific_domains.jl
index e8e0b88..1803e36 100644
--- a/test/test_specific_domains.jl
+++ b/test/test_specific_domains.jl
@@ -717,6 +717,7 @@ end
         @test boundary(D) == UnitCircle()
         @test dimension(D) == 2
         @test boundingbox(D) == ProductDomain(ChebyshevInterval(), ChebyshevInterval())
+        @test normal(UnitDisk(), [sqrt(2)/2, sqrt(2)/2]) ≈ [sqrt(2)/2, sqrt(2)/2]
 
         @test boundingbox(UnitBall{Float64}()) == ChebyshevInterval()
 
@@ -1126,8 +1127,6 @@ end
         @test pd isa Domain{Vector{Float64}}
         @test forward_map(pd) == m
         @test forward_map(pd, 0.4) ≈ m(0.4)
-        @test inverse_map(pd)(m(0.4)) ≈ 0.4
-        @test inverse_map(pd, m(0.4)) ≈ 0.4
         @test mapfrom_canonical(pd) == m
         @test canonicaldomain(pd) == 0..1
         @test boundary(pd) isa UnionDomain
@@ -1545,6 +1544,13 @@ end
         @test Set(corners(UnitCube())) == Set([ [0,0,0], [1,0,0], [0,1,0], [0,0,1], [1,1,0], [1,0,1], [0,1,1], [1,1,1]])
         @test boundary(boundary(UnitSquare())) == UnionDomain(map(t->Point(SVector{2}(t)), corners(UnitSquare())))
         @test boundary(boundary(boundary(UnitCube()))) == UnionDomain(map(t->Point(SVector{3}(t)), corners(UnitCube())))
+        @test UnitCube() == UnitCube(3)
+        @test boundary(UnitCube()) == boundary(UnitCube(3))
+        @test boundary(boundary(UnitCube())) == boundary(boundary(UnitCube(3)))
+        @test boundary(boundary(boundary(UnitCube()))) == boundary(boundary(boundary(UnitCube(3))))
+        @test UnitSquare() == UnitCube(2)
+        @test boundary(UnitSquare()) == boundary(UnitCube(2))
+        @test boundary(boundary(UnitSquare())) == boundary(boundary(UnitCube(2)))
 
         @test StaticUnitCube() == UnitCube()
         @test StaticUnitCube(Val(3)) == UnitCube()

From 27c39780efc67ee96acbfaa9c6be10b038c39b6c Mon Sep 17 00:00:00 2001
From: Daan Huybrechs <daan.huybrechs@kuleuven.be>
Date: Tue, 22 Jun 2021 13:01:14 +0200
Subject: [PATCH 3/5] improve support for boundary, normals and hashes

---
 src/DomainSets.jl             |  2 +-
 src/domains/ball.jl           |  2 +-
 src/domains/cube.jl           | 53 ++++++++++++++------------
 src/domains/interval.jl       | 21 +++++++++++
 src/domains/simplex.jl        | 71 +++++++++++++++++++++++++++++++++--
 src/generic/setoperations.jl  |  1 +
 src/maps/affine.jl            |  5 +--
 src/maps/basic.jl             |  2 +-
 src/maps/composite.jl         |  2 +-
 src/maps/product.jl           |  4 +-
 src/util/common.jl            | 15 ++++++--
 test/test_maps.jl             | 11 +++++-
 test/test_specific_domains.jl |  2 +
 13 files changed, 147 insertions(+), 44 deletions(-)

diff --git a/src/DomainSets.jl b/src/DomainSets.jl
index c26344e..b995919 100644
--- a/src/DomainSets.jl
+++ b/src/DomainSets.jl
@@ -101,7 +101,7 @@ export Domain, EuclideanDomain, VectorDomain,
     interior, closure,
     volume,
     point_in_domain,
-    normal, tangents,
+    normal, tangents, distance_to,
     canonicaldomain, mapto_canonical, mapfrom_canonical, hascanonicaldomain,
     mapto,
     parameterdomain, parameterization, hasparameterization,
diff --git a/src/domains/ball.jl b/src/domains/ball.jl
index 84e1b41..83d8ee8 100644
--- a/src/domains/ball.jl
+++ b/src/domains/ball.jl
@@ -295,7 +295,7 @@ normal(d::Sphere, x) = (x-center(x))/norm(x-center(x))
 
 distance_to(d::Sphere, x) = abs(norm(x-center(d))-radius(d))
 
-point_in_domain(d::Sphere) = center(d) + unitvector(d)
+point_in_domain(d::Sphere) = center(d) + unitvector(d, 1)
 
 ==(d1::Sphere, d2::Sphere) = radius(d1)==radius(d2) && center(d1)==center(d2)
 
diff --git a/src/domains/cube.jl b/src/domains/cube.jl
index 0a322da..cb667db 100644
--- a/src/domains/cube.jl
+++ b/src/domains/cube.jl
@@ -12,14 +12,14 @@ function corners(d::HyperRectangle)
 	left = infimum(d)
 	right = supremum(d)
     N = length(left)
-    corners = [zeros(numtype(d), N) for i in 1:2^N]
+    corners = [similar(point_in_domain(d)) for i in 1:2^N]
     # All possible permutations of the corners
     for i=1:2^length(left)
         for j=1:N
             corners[i][j] = ((i>>(j-1))%2==0) ? left[j] : right[j]
         end
     end
-    corners
+    [convert(eltype(d), p) for p in corners]
 end
 
 # map the interval [a,b] to a cube defined by c (bottom-left) and d (top-right)
@@ -79,17 +79,21 @@ function cube_face_map(a::Vector, b::Vector, c::Vector, d::Vector, dim, dimval)
 	AffineMap(A, B)
 end
 
+# The boundary of a rectangle is a collection of mapped lower-dimensional rectangles.
+# Dimension 2 is a special case, because the lower dimension is 1 where we use
+# scalars instead of vectors
 function boundary(d::HyperRectangle{SVector{2,T}}) where {T}
 	left = infimum(d)
 	right = supremum(d)
 	x1 = left[1]; y1 = left[2]; x2 = right[1]; y2 = right[2]
 	d_unit = UnitInterval{T}()
-	faces = [
-		ParametricDomain(cube_face_map(zero(T), one(T), SVector(x1,y1), SVector(x2,y1)), d_unit),
-		ParametricDomain(cube_face_map(zero(T), one(T), SVector(x2,y1), SVector(x2,y2)), d_unit),
-		ParametricDomain(cube_face_map(zero(T), one(T), SVector(x2,y2), SVector(x1,y2)), d_unit),
-		ParametricDomain(cube_face_map(zero(T), one(T), SVector(x1,y2), SVector(x1,y1)), d_unit)
+	maps = [
+		cube_face_map(zero(T), one(T), SVector(x1,y1), SVector(x2,y1)),
+		cube_face_map(zero(T), one(T), SVector(x2,y1), SVector(x2,y2)),
+		cube_face_map(zero(T), one(T), SVector(x2,y2), SVector(x1,y2)),
+		cube_face_map(zero(T), one(T), SVector(x1,y2), SVector(x1,y1))
 	]
+	faces = map(m -> ParametricDomain(m, d_unit), maps)
 	UnionDomain(faces)
 end
 
@@ -100,13 +104,14 @@ function boundary(d::HyperRectangle{SVector{N,T}}) where {N,T}
 	left1 = infimum(d_unit)
 	right1 = supremum(d_unit)
 
-	face1 = ParametricDomain(cube_face_map(left1, right1, left2, right2, 1, left2[1]), d_unit)
-	D = typeof(face1)
-	faces = D[]
+	map1 = cube_face_map(left1, right1, left2, right2, 1, left2[1])
+	MAP = typeof(map1)
+	maps = MAP[]
 	for dim in 1:N
-		push!(faces, ParametricDomain(cube_face_map(left1, right1, left2, right2, dim, left2[dim]), d_unit))
-		push!(faces, ParametricDomain(cube_face_map(left1, right1, left2, right2, dim, right2[dim]), d_unit))
+		push!(maps, cube_face_map(left1, right1, left2, right2, dim, left2[dim]))
+		push!(maps, cube_face_map(left1, right1, left2, right2, dim, right2[dim]))
 	end
+	faces = map(m -> ParametricDomain(m, d_unit), maps)
 	UnionDomain(faces)
 end
 
@@ -116,13 +121,12 @@ function boundary(d::HyperRectangle{Vector{T}}) where {T}
 		right = supremum(d)
 		x1 = left[1]; y1 = left[2]; x2 = right[1]; y2 = right[2]
 		d_unit = UnitInterval{T}()
-		faces = [
-			ParametricDomain(cube_face_map(zero(T), one(T), [x1,y1], [x2,y1]), d_unit),
-			ParametricDomain(cube_face_map(zero(T), one(T), [x2,y1], [x2,y2]), d_unit),
-			ParametricDomain(cube_face_map(zero(T), one(T), [x2,y2], [x1,y2]), d_unit),
-			ParametricDomain(cube_face_map(zero(T), one(T), [x1,y2], [x1,y1]), d_unit)
+		maps = [
+			cube_face_map(zero(T), one(T), [x1,y1], [x2,y1]),
+			cube_face_map(zero(T), one(T), [x2,y1], [x2,y2]),
+			cube_face_map(zero(T), one(T), [x2,y2], [x1,y2]),
+			cube_face_map(zero(T), one(T), [x1,y2], [x1,y1])
 		]
-		UnionDomain(faces)
 	else
 		left2 = infimum(d)
 		right2 = supremum(d)
@@ -130,15 +134,16 @@ function boundary(d::HyperRectangle{Vector{T}}) where {T}
 		left1 = infimum(d_unit)
 		right1 = supremum(d_unit)
 
-		face1 = ParametricDomain(cube_face_map(left1, right1, left2, right2, 1, left2[1]), d_unit)
-		D = typeof(face1)
-		faces = D[]
+		map1 = cube_face_map(left1, right1, left2, right2, 1, left2[1])
+		MAP = typeof(map1)
+		maps = MAP[]
 		for dim in 1:dimension(d)
-			push!(faces, ParametricDomain(cube_face_map(left1, right1, left2, right2, dim, left2[dim]), d_unit))
-			push!(faces, ParametricDomain(cube_face_map(left1, right1, left2, right2, dim, right2[dim]), d_unit))
+			push!(maps, cube_face_map(left1, right1, left2, right2, dim, left2[dim]))
+			push!(maps, cube_face_map(left1, right1, left2, right2, dim, right2[dim]))
 		end
-		UnionDomain(faces)
 	end
+	faces = map(m -> ParametricDomain(m, d_unit), maps)
+	UnionDomain(faces)
 end
 
 
diff --git a/src/domains/interval.jl b/src/domains/interval.jl
index 8819138..188fd6b 100644
--- a/src/domains/interval.jl
+++ b/src/domains/interval.jl
@@ -28,6 +28,27 @@ function point_in_domain(d::AbstractInterval)
     mean(d)
 end
 
+# For an interval of integers, try to find an integer point
+function point_in_domain(d::AbstractInterval{T}) where {T<:Integer}
+    isempty(d) && throw(BoundsError())
+    x = round(T, mean(d))
+    if x ∈ d
+        x
+    else
+        # the mean is not inside the interval when rounded, this means
+        # we have to choose an endpoint
+        if isleftclosed(d)
+            leftendpoint(d)
+        elseif isrightclosed(d)
+            rightendpoint(d)
+        else
+            # the interval is open and contains no integers in its interior
+            throw(BoundsError())
+        end
+    end
+end
+
+
 isapprox(d1::AbstractInterval, d2::AbstractInterval) =
     isapprox(leftendpoint(d1), leftendpoint(d2); atol=default_tolerance(d1)) &&
     isapprox(rightendpoint(d1), rightendpoint(d2); atol=default_tolerance(d1))
diff --git a/src/domains/simplex.jl b/src/domains/simplex.jl
index 02f1316..9ae44c5 100644
--- a/src/domains/simplex.jl
+++ b/src/domains/simplex.jl
@@ -48,9 +48,25 @@ hash(d::UnitSimplex, h::UInt) = hashrec("UnitSimplex", isclosedset(d), dimension
 
 boundingbox(d::UnitSimplex{T}) where {T} = UnitCube{T}(dimension(d))
 
+boundary(d::UnitSimplex{T,:open}) where {T} = EmptySpace{T}()
+
+
 distance_to(d::UnitSimplex, x) = x ∈ d ? zero(prectype(d)) : minimum(distance_to(el, x) for el in components(boundary(d)))
 
+function normal(d::UnitSimplex, x)
+    z = similar(x)
+    fill!(z, 0)
+    if sum(x) ≈ 1
+        fill!(z, one(eltype(z))/sqrt(length(z)))
+    else
+        index = findmax(x)[2]
+        z[index] = -1
+    end
+    return convert(eltype(d), z/norm(z))
+end
+
 
+"A unit simplex whose dimension is determined by its element type."
 struct StaticUnitSimplex{T,C} <: UnitSimplex{T,C}
 end
 
@@ -72,7 +88,18 @@ const EuclideanUnitSimplex{N,T,C} = StaticUnitSimplex{SVector{N,T},C}
 
 EuclideanUnitSimplex{N}() where {N} = EuclideanUnitSimplex{N,Float64}()
 
+## A StaticUnitSimplex{<:Number} equals the interval [0,1]  (open or closed)
+convert(::Type{Interval}, d::StaticUnitSimplex{T,:closed}) where {T <: Number} =
+    UnitInterval{T}()
+convert(::Type{Interval}, d::StaticUnitSimplex{T,:open}) where {T <: Number} =
+    OpenInterval{T}(0, 1)
+
+canonicaldomain(::Equal, d::StaticUnitSimplex{T}) where {T<:Number} = convert(Interval, d)
+
+boundary(d::StaticUnitSimplex{T}) where {T<:Number} = boundary(convert(Interval, d))
+
 
+"A unit simplex with vector elements with variable dimension determined by a field."
 struct DynamicUnitSimplex{T,C} <: UnitSimplex{T,C}
     dimension   ::  Int
 
@@ -97,10 +124,7 @@ show(io::IO, d::VectorUnitSimplex{Float64,:closed}) = print(io, "UnitSimplex($(d
 center(d::EuclideanUnitSimplex{N,T}) where {N,T} = ones(SVector{N,T})/N
 center(d::VectorUnitSimplex{T}) where {T} = ones(T, dimension(d))/dimension(d)
 
-corners(::EuclideanUnitSimplex{N,T}) where {N,T} =
-    vcat([zero(SVector{N,T})], [ SVector{N,T}(ntuple( i -> i==j, N)) for j in 1:N])
-corners(d::VectorUnitSimplex{T}) where {T} =
-    vcat([zeros(T,dimension(d))], [ (z = zeros(T, dimension(d)); z[j]=1; z) for j in 1:dimension(d)])
+corners(d::UnitSimplex) = vcat([origin(d)], [ unitvector(d, i) for i in 1:dimension(d)])
 
 interior(d::EuclideanUnitSimplex{N,T}) where {N,T} = EuclideanUnitSimplex{N,T,:open}()
 closure(d::EuclideanUnitSimplex{N,T}) where {N,T} = EuclideanUnitSimplex{N,T,:closed}()
@@ -116,3 +140,42 @@ similardomain(d::StaticUnitSimplex{S,C}, ::Type{T}) where {S,T,C} =
     StaticUnitSimplex{T,C}()
 similardomain(d::DynamicUnitSimplex{S,C}, ::Type{T}) where {S,T,C} =
     DynamicUnitSimplex{T,C}(d.dimension)
+
+
+simplex_face_map(a::Number, b::Number, c::SVector{2}, d::SVector{2}) =
+	AffineMap((d-c)/(b-a), c - (d-c)/(b-a)*a)
+function simplex_face_map(a::Number, b::Number, c::Vector, d::Vector)
+	@assert length(c) == length(d) == 2
+	AffineMap((d-c)/(b-a), c - (d-c)/(b-a)*a)
+end
+
+function boundary(d::StaticUnitSimplex{SVector{2,T},:closed}) where {T}
+    d0 = UnitInterval{T}()
+	T0 = zero(T)
+	T1 = one(T)
+	maps = [
+		simplex_face_map(T0, T1, SVector(T0,T0), SVector(T1,T0)),
+		simplex_face_map(T0, T1, SVector(T1,T0), SVector(T0,T1)),
+		simplex_face_map(T0, T1, SVector(T0,T1), SVector(T0,T0))
+	]
+	faces = map(m -> ParametricDomain(m, d0), maps)
+	UnionDomain(faces)
+end
+
+# function boundary(d::StaticUnitSimplex{SVector{N,T},:closed}) where {N,T}
+# 	left2 = infimum(d)
+# 	right2 = supremum(d)
+# 	d0 = UnitSimplex{SVector{N-1,T},:closed}()
+# 	T0 = zero(T)
+# 	T1 = one(T)
+#
+# 	map1 = cube_face_map(left1, right1, left2, right2, 1, left2[1])
+# 	MAP = typeof(map1)
+# 	maps = MAP[]
+# 	for dim in 1:N
+# 		push!(maps, cube_face_map(left1, right1, left2, right2, dim, left2[dim]))
+# 		push!(maps, cube_face_map(left1, right1, left2, right2, dim, right2[dim]))
+# 	end
+# 	faces = map(m -> ParametricDomain(m, d_unit), maps)
+# 	UnionDomain(faces)
+# end
diff --git a/src/generic/setoperations.jl b/src/generic/setoperations.jl
index d6af007..dda05b5 100644
--- a/src/generic/setoperations.jl
+++ b/src/generic/setoperations.jl
@@ -248,6 +248,7 @@ similardomain(d::IntersectDomain, ::Type{T}) where {T} =
     IntersectDomain(convert.(Domain{T}, components(d)))
 
 ==(a::IntersectDomain, b::IntersectDomain) = Set(components(a)) == Set(components(b))
+hash(d::IntersectDomain, h::UInt) = hashrec("IntersectDomain", Set(components(d)), h)
 
 boundingbox(d::IntersectDomain) = intersectbox(map(boundingbox, components(d))...)
 
diff --git a/src/maps/affine.jl b/src/maps/affine.jl
index 388a0a0..fca3550 100644
--- a/src/maps/affine.jl
+++ b/src/maps/affine.jl
@@ -51,13 +51,12 @@ isaffine(m::AbstractAffineMap) = true
 
 ==(m1::AbstractAffineMap, m2::AbstractAffineMap) =
     matrix(m1) == matrix(m2) && vector(m1) == vector(m2)
+hash(m::AbstractAffineMap, h::UInt) = hashrec(h, matrix(m), vector(m))
 
 ==(m1::AbstractAffineMap, m2::IdentityMap) =
     islinear(m1) && matrix(m1) == matrix(m2)
 ==(m1::IdentityMap, m2::AbstractAffineMap) = m2==m1
 
-hash(m::AbstractAffineMap, h::UInt) = hashrec(h, matrix(m), vector(m))
-
 mapsize(m::AbstractAffineMap) = _affine_mapsize(m, domaintype(m), unsafe_matrix(m), unsafe_vector(m))
 _affine_mapsize(m::AbstractAffineMap, T, A::AbstractArray, b) = size(A)
 _affine_mapsize(m::AbstractAffineMap, T, A::AbstractVector, b::AbstractVector) = (length(A),)
@@ -127,7 +126,6 @@ isreal(m::LinearMap) = _isreal(m, unsafe_matrix(m))
 _isreal(m::LinearMap, A) = isreal(A)
 
 ==(m1::LinearMap, m2::LinearMap) = matrix(m1) == matrix(m2)
-hash(m::LinearMap, h::UInt) = hash(matrix(m), h)
 
 # inverse should be called only on square maps, otherwise use
 # leftinverse or rightinverse in order to use pinv instead of inv
@@ -314,7 +312,6 @@ jacdet(m::Translation{T}, x) where {T} = one(eltype(T))
 isreal(m::Translation) = isreal(unsafe_vector(m))
 
 ==(m1::Translation, m2::Translation) = unsafe_vector(m1)==unsafe_vector(m2)
-hash(m::Translation, h::UInt) = hash(unsafe_vector(m), h)
 
 similarmap(m::Translation, ::Type{T}) where {T} = Translation{T}(m.b)
 
diff --git a/src/maps/basic.jl b/src/maps/basic.jl
index d5859e1..c528961 100644
--- a/src/maps/basic.jl
+++ b/src/maps/basic.jl
@@ -83,7 +83,7 @@ similarmap(m::DynamicIdentityMap, ::Type{T}) where {T<:StaticTypes} =
     StaticIdentityMap{T}()
 
 ==(m1::DynamicIdentityMap, m2::DynamicIdentityMap) = m1.dimension == m2.dimension
-hash(m::DynamicIdentityMap, h::UInt) = hashrec("DynamicIdentityMap", dimension(m), h)
+hash(m::DynamicIdentityMap, h::UInt) = hashrec("DynamicIdentityMap", m.dimension, h)
 
 
 "The supertype of constant maps from `T` to `U`."
diff --git a/src/maps/composite.jl b/src/maps/composite.jl
index ff30755..3e02b78 100644
--- a/src/maps/composite.jl
+++ b/src/maps/composite.jl
@@ -77,7 +77,7 @@ composedmap2(m1, m2::ComposedMap) = ComposedMap(m1, components(m2)...)
 
 ==(m1::ComposedMap, m2::ComposedMap) =
     ncomponents(m1) == ncomponents(m2) && all(map(isequal, components(m1), components(m2)))
-hash(m::ComposedMap, h::UInt) = hashrec("ComposedMap", ncomponents(m), components(m), h)
+hash(m::ComposedMap, h::UInt) = hashrec("ComposedMap", collect(components(m)), h)
 
 Display.combinationsymbol(m::ComposedMap) = Display.Symbol('∘')
 Display.displaystencil(m::ComposedMap) =
diff --git a/src/maps/product.jl b/src/maps/product.jl
index c0bb227..85fbb77 100644
--- a/src/maps/product.jl
+++ b/src/maps/product.jl
@@ -75,7 +75,7 @@ end
 mapsize(m::ProductMap) = (sum(t->mapsize(t,1), components(m)), sum(t->mapsize(t,2), components(m)))
 
 ==(m1::ProductMap, m2::ProductMap) = all(map(isequal, components(m1), components(m2)))
-hash(m::ProductMap, h::UInt) = hashrec("ProductMap", components(m), h)
+hash(m::ProductMap, h::UInt) = hashrec("ProductMap", collect(components(m)), h)
 
 Display.combinationsymbol(m::ProductMap) = Display.Symbol('⊗')
 Display.displaystencil(m::ProductMap) = composite_displaystencil(m)
@@ -159,7 +159,7 @@ end
 # the Jacobian is a diagonal matrix
 toexternalmatrix(m::VectorProductMap, matrices) = Diagonal(matrices)
 
-dimension(m::VectorProductMap) = length(m.maps)
+mapsize(m::VectorProductMap) = (length(m.maps), length(m.maps))
 
 """
 A `TupleProductMap` is a product map with all components collected in a tuple.
diff --git a/src/util/common.jl b/src/util/common.jl
index 1159444..182575c 100644
--- a/src/util/common.jl
+++ b/src/util/common.jl
@@ -12,13 +12,20 @@ euclideandimension(::Type{T}) where {N,T <: StaticVector{N}} = N
 euclideandimension(::Type{T}) where {N,T <: NTuple{N,Any}} = N
 # Does not apply to Vector{T}: we don't know its dimension
 
-unitvector(d::Domain{T}) where {N,S,T<:SVector{N,S}} = SVector{N,S}(ntuple(i -> i==1, N))
-function unitvector(d::Domain{T}) where {T<:AbstractVector}
+unitvector(d::Domain{T}, dim) where {N,S,T<:SVector{N,S}} = SVector{N,S}(ntuple(i -> i==dim, N))
+function unitvector(d::Domain{T}, dim) where {T<:AbstractVector}
     p = zeros(eltype(T), dimension(d))
-    p[1] = 1
+    p[dim] = 1
     p
 end
-unitvector(d::Domain{T}) where {T<:Number} = one(T)
+unitvector(d::Domain{T}, dim) where {T<:Number} = (@assert dim==1; one(T))
+
+origin(d::Domain{T}) where {T <: StaticTypes} = zero(T)
+function origin(d::Domain{T}) where {T <: AbstractVector}
+	p = similar(point_in_domain(d))
+	fill!(p, 0)
+	convert(T, p)
+end
 
 "Apply the `hash` function recursively to the given arguments."
 hashrec(h) = hash(h)
diff --git a/test/test_maps.jl b/test/test_maps.jl
index c0b4582..7bc66b1 100644
--- a/test/test_maps.jl
+++ b/test/test_maps.jl
@@ -133,7 +133,7 @@ function test_generic_jacobian(m)
     end
 end
 
-# Test a map m with dimensions n
+# generic map test suite
 function test_generic_map(m)
     @test convert(Map{domaintype(m)}, m) == m
 
@@ -491,6 +491,7 @@ function test_identity_map(T)
     test_generic_map(i1)
     test_generic_map(i2)
     @test i1 == i2
+    @test hash(i1) == hash(i2)
     @test islinear(i1)
     @test isaffine(i1)
     @test convert(StaticIdentityMap{SVector{2,T}}, i1) === i2
@@ -509,12 +510,15 @@ function test_identity_map(T)
     test_generic_map(i3)
     r = rand(T, 10)
     @test i3(r) ≈ r
+    @test hash(i3) == hash(VectorIdentityMap{Int}(10))
 
     @test IdentityMap() ∘ LinearMap(2) == LinearMap(2.0)
 end
 
 function test_basic_maps(T)
     @test UnityMap{SVector{2,Float64}}() == UnityMap{SVector{2,Float64},Float64}()
+    @test hash(ConstantMap(2)) == hash(ConstantMap(2.0))
+    @test DomainSets.absmap(ConstantMap(-2)) == ConstantMap(2)
 end
 
 function test_composite_map(T)
@@ -602,6 +606,9 @@ function test_product_map(T)
     test_generic_map(m4)
     @test m4(SVector(r1,r2,r3,r4,r5)) ≈ SVector(m1(SVector(r1,r2))...,m2(SVector(r3,r4,r5))...)
 
+    @test ProductMap(ma, mb) == ProductMap([ma,mb])
+    @test hash(ProductMap(ma, mb)) == hash(ProductMap([ma,mb]))
+
     m5 = productmap(AffineMap(SMatrix{2,2,T}(1.0,2,3,4), SVector{2,T}(1,3)), LinearMap{T}(2.0))
     @test domaintype(component(m5,1)) == SVector{2,T}
     @test domaintype(component(m5,2)) == T
@@ -611,7 +618,7 @@ function test_product_map(T)
 
     m6 = ProductMap([ma,mb])
     @test m6 isa DomainSets.VectorProductMap
-    @test dimension(m6) == 2
+    @test mapsize(m6) == (2,2)
     @test convert(Map{SVector{2,T}}, m6) isa DomainSets.VcatMap
     test_generic_map(m6)
 
diff --git a/test/test_specific_domains.jl b/test/test_specific_domains.jl
index 1803e36..a2cd12c 100644
--- a/test/test_specific_domains.jl
+++ b/test/test_specific_domains.jl
@@ -139,6 +139,8 @@ end
         @test canonicaldomain(d) == Point(0.0)
         @test mapfrom_canonical(d) == Translation(1.0)
 
+        @test distance_to(d, 0.5) == abs(0.5-d.x)
+
         @test d .+ 1 == Domain(2.0)
         @test 1 .+ d == Domain(2.0)
         @test 1 .- d == Domain(0.0)

From 4e4004583bb501067e72436ad37c6b9ea52f84a0 Mon Sep 17 00:00:00 2001
From: Daan Huybrechs <daan.huybrechs@kuleuven.be>
Date: Tue, 22 Jun 2021 13:21:46 +0200
Subject: [PATCH 4/5] preserve ranges when changing eltype

---
 src/generic/domain.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/generic/domain.jl b/src/generic/domain.jl
index a588579..fbdefa0 100644
--- a/src/generic/domain.jl
+++ b/src/generic/domain.jl
@@ -38,6 +38,7 @@ _convert_eltype(::Type{T}, d, ::Type{S}) where {S,T} =
     error("Don't know how to convert the `eltype` of $(d).")
 # Some standard cases
 convert_eltype(::Type{T}, d::AbstractArray) where {T} = convert(AbstractArray{T}, d)
+convert_eltype(::Type{T}, d::AbstractRange) where {T} = map(T, d)
 convert_eltype(::Type{T}, d::Set) where {T} = convert(Set{T}, d)
 
 promote(d1::Domain, d2::Domain) = promote_domains((d1,d2))

From 6716483b9a620a527b2018a107ac1d8d06176459 Mon Sep 17 00:00:00 2001
From: Daan Huybrechs <daan.huybrechs@kuleuven.be>
Date: Sat, 14 Aug 2021 12:49:51 +0200
Subject: [PATCH 5/5] use left/rightendpoint instead of infimum and supremum

---
 src/domains/cube.jl | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/src/domains/cube.jl b/src/domains/cube.jl
index cb667db..b2ad15a 100644
--- a/src/domains/cube.jl
+++ b/src/domains/cube.jl
@@ -299,26 +299,26 @@ Rectangle(a::NTuple{N,T}, b::NTuple{N,T}) where {N,T} =
 Rectangle(domains::Tuple) = Rectangle(domains...)
 Rectangle(domains::ClosedInterval...) = Rectangle(promote_domains(domains)...)
 Rectangle(domains::ClosedInterval{T}...) where {T} =
-    Rectangle(map(infimum, domains), map(supremum, domains))
+    Rectangle(map(leftendpoint, domains), map(rightendpoint, domains))
 Rectangle(domains::AbstractVector{<:ClosedInterval}) =
-    Rectangle(map(infimum, domains), map(supremum, domains))
+    Rectangle(map(leftendpoint, domains), map(rightendpoint, domains))
 Rectangle(domains::Domain...) =
     error("The Rectangle constructor expects two points or a list of intervals (closed).")
 
 Rectangle{T}(domains::Tuple) where {T} = Rectangle{T}(domains...)
 Rectangle{T}(domains::ClosedInterval...) where {T} =
-	Rectangle{T}(map(infimum, domains), map(supremum, domains))
+	Rectangle{T}(map(leftendpoint, domains), map(rightendpoint, domains))
 Rectangle{T}(domains::AbstractVector{<:ClosedInterval}) where {T} =
-    Rectangle{T}(map(infimum, domains), map(supremum, domains))
+    Rectangle{T}(map(leftendpoint, domains), map(rightendpoint, domains))
 Rectangle{T}(domains::Domain...) where {T} =
     error("The Rectangle constructor expects two points or a list of intervals (closed).")
 
 ProductDomain(domains::ClosedInterval...) = Rectangle(domains...)
 ProductDomain(domains::AbstractVector{<:ClosedInterval}) =
-    Rectangle(map(infimum, domains), map(supremum, domains))
+    Rectangle(map(leftendpoint, domains), map(rightendpoint, domains))
 ProductDomain{T}(domains::ClosedInterval...) where {T} = Rectangle{T}(domains...)
 ProductDomain{T}(domains::AbstractVector{<:ClosedInterval}) where {T} =
-    Rectangle{T}(map(infimum, domains), map(supremum, domains))
+    Rectangle{T}(map(leftendpoint, domains), map(rightendpoint, domains))