implement boundary and normal of domain

src/DomainSets.jl

@@ -101,6 +101,7 @@ export Domain, EuclideanDomain, VectorDomain,
     interior, closure,
     volume,
     point_in_domain,
+    normal, tangents, distance_to,
     canonicaldomain, mapto_canonical, mapfrom_canonical, hascanonicaldomain,
     mapto,
     parameterdomain, parameterization, hasparameterization,
@@ -212,6 +213,7 @@ include("domains/cube.jl")
 include("domains/indicator.jl")
 include("domains/boundingbox.jl")
 
+include("applications/coordinates.jl")
 include("applications/rotation.jl")
 
 end # module

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.θ))

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,7 +291,11 @@ isempty(::Sphere) = false
 isclosedset(::Sphere) = true
 isopenset(::Sphere) = false
 
-point_in_domain(d::Sphere) = center(d) + unitvector(d)
+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, 1)
 
 ==(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

src/domains/cube.jl

@@ -2,21 +2,148 @@
 "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."
 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)
+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 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)
+	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 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
+
+# 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}()
+	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
+
+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)
+
+	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
+
+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}()
+		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])
+		]
+	else
+		left2 = infimum(d)
+		right2 = supremum(d)
+		d_unit = UnitCube(dimension(d)-1)
+		left1 = infimum(d_unit)
+		right1 = supremum(d_unit)
+
+		map1 = cube_face_map(left1, right1, left2, right2, 1, left2[1])
+		MAP = typeof(map1)
+		maps = MAP[]
+		for dim in 1:dimension(d)
+			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
+	end
+	faces = map(m -> ParametricDomain(m, d_unit), maps)
+	UnionDomain(faces)
 end
 
 
@@ -172,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))
 
 
 

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))

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) =
@@ -25,12 +28,38 @@ 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))
 
 
 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
 

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}