diff --git a/Project.toml b/Project.toml
index fa0a9c2..fa90aa1 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "DomainSets"
 uuid = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
-version = "0.6.7"
+version = "0.7.0"
 
 [deps]
 CompositeTypes = "b152e2b5-7a66-4b01-a709-34e65c35f657"
@@ -11,15 +11,15 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
+julia = "1.6"
 CompositeTypes = "0.1.2"
 IntervalSets = "0.7.4"
-StaticArrays = "0.12.2, 1"
 StableRNGs = "1"
-julia = "1.6"
+StaticArrays = "0.12.2, 1"
 
 [extras]
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
 test = ["Test", "StableRNGs"]
diff --git a/docs/src/interface.md b/docs/src/interface.md
new file mode 100644
index 0000000..e69de29
diff --git a/src/DomainSets.jl b/src/DomainSets.jl
index 8fcb99f..a895ebd 100644
--- a/src/DomainSets.jl
+++ b/src/DomainSets.jl
@@ -43,7 +43,6 @@ export ..
 
 import CompositeTypes: component, components
 
-
 ################################
 ## Exhaustive list of exports
 ################################
@@ -81,6 +80,8 @@ export AffineMap, Translation, LinearMap,
 
 ## Generic domains
 
+@deprecate point_in_domain(d) choice(d)
+
 # from generic/domain.jl
 export Domain, EuclideanDomain, VectorDomain,
     dimension,
@@ -88,7 +89,7 @@ export Domain, EuclideanDomain, VectorDomain,
     isopenset, isclosedset, iscompact,
     boundary, ∂,
     interior, closure,
-    point_in_domain
+    choice
 
 # from generic/geometry.jl
 export boundingbox,
@@ -192,6 +193,7 @@ include("maps/basic.jl")
 include("maps/affine.jl")
 include("maps/arithmetics.jl")
 
+include("generic/core.jl")
 include("generic/domain.jl")
 include("generic/geometry.jl")
 include("generic/canonical.jl")
diff --git a/src/domains/ball.jl b/src/domains/ball.jl
index f68f51e..272dd53 100644
--- a/src/domains/ball.jl
+++ b/src/domains/ball.jl
@@ -56,7 +56,7 @@ isequaldomain(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)
+function issubset_domain(d1::Ball, d2::Ball)
     if dimension(d1) == dimension(d2)
         if center(d1) == center(d2)
             if radius(d1) < radius(d2)
@@ -84,7 +84,7 @@ 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)
+choice(d::Ball) = center(d)
 
 
 "The unit ball."
@@ -119,7 +119,7 @@ approx_indomain(x, d::ClosedUnitBall, tolerance) = norm(x) <= 1+tolerance
 isequaldomain(d1::UnitBall, d2::UnitBall) = isclosedset(d1)==isclosedset(d2) &&
     dimension(d1) == dimension(d2)
 
-issubset(d1::UnitBall, d2::UnitBall) =
+issubset_domain(d1::UnitBall, d2::UnitBall) =
     dimension(d1) == dimension(d2) && (isclosedset(d2) || isopenset(d1))
 
 convert(::Type{SublevelSet}, d::UnitBall{T,C}) where {T,C} =
@@ -304,7 +304,7 @@ normal(d::Sphere, x) = (x-center(d))/norm(x-center(d))
 
 distance_to(d::Sphere, x) = abs(norm(x-center(d))-radius(d))
 
-point_in_domain(d::Sphere) = center(d) + unitvector(d, 1)
+choice(d::Sphere) = center(d) + unitvector(d, 1)
 
 isequaldomain(d1::Sphere, d2::Sphere) =
     radius(d1)==radius(d2) && center(d1)==center(d2)
diff --git a/src/domains/boundingbox.jl b/src/domains/boundingbox.jl
index bc30d11..e5ea390 100644
--- a/src/domains/boundingbox.jl
+++ b/src/domains/boundingbox.jl
@@ -3,12 +3,10 @@ boundingbox(d::Vector{T}) where {T <: Number} = minimum(d)..maximum(d)
 boundingbox(d::Set{T}) where {T<:Number} = minimum(d)..maximum(d)
 
 "Return the bounding box of the union of two or more bounding boxes."
-unionbox(d::Domain) = d
-unionbox(d1::Domain, d2::Domain) = unionbox(promote_domains(d1, d2)...)
-unionbox(d1::Domain, d2::Domain, domains...) =
+unionbox(d) = d
+unionbox(d1, d2) = unionbox1(promote_domains(d1, d2)...)
+unionbox(d1, d2, domains...) =
     unionbox(unionbox(d1,d2), domains...)
-
-unionbox(d1::Domain{T}, d2::Domain{T}) where {T} = unionbox1(d1, d2)
 unionbox1(d1, d2) = unionbox2(d1, d2)
 unionbox2(d1, d2) = fullspace(d1)
 unionbox1(d1::EmptySpace, d2) = d2
@@ -16,26 +14,28 @@ unionbox1(d1::FullSpace, d2) = d1
 unionbox2(d1, d2::EmptySpace) = d1
 unionbox2(d1, d2::FullSpace) = d2
 
-unionbox(d1::D, d2::D) where {D<:FixedInterval} = d1
+unionbox1(d1::D, d2::D) where {D<:FixedInterval} = d1
 
-function unionbox(d1::AbstractInterval{T}, d2::AbstractInterval{T}) where {T}
+function unionbox1(d1::AbstractInterval, d2::AbstractInterval)
     a, b = endpoints(d1)
     c, d = endpoints(d2)
     A = min(a, c)
     B = max(b, d)
-    isinf(A) && isinf(B) ? fullspace(T) : A..B
+    isinf(A) && isinf(B) ? fullspace(d1) : A..B
 end
 
-unionbox(d1::HyperRectangle{T}, d2::HyperRectangle{T}) where {T} =
+function unionbox(d1::HyperRectangle, d2::HyperRectangle)
+    @assert dimension(d1) == dimension(d2)
+    T = promote_type(eltype(d1),eltype(d2))
     Rectangle{T}(map(unionbox, components(d1), components(d2)))
+end
 
 "Return the bounding box of the intersection of two or more bounding boxes."
-intersectbox(d::Domain) = d
-intersectbox(d1::Domain, d2::Domain) = intersectbox(promote_domains(d1, d2)...)
-intersectbox(d1::Domain, d2::Domain, domains...) =
+intersectbox(d) = d
+intersectbox(d1, d2) = intersectbox1(promote_domains(d1, d2)...)
+intersectbox(d1, d2, domains...) =
     intersectbox(intersectbox(d1,d2), domains...)
 
-intersectbox(d1::Domain{T}, d2::Domain{T}) where {T} = intersectbox1(d1, d2)
 intersectbox1(d1, d2) = intersectbox2(d1, d2)
 intersectbox2(d1, d2) = fullspace(d1)
 intersectbox1(d1::EmptySpace, d2) = d1
@@ -43,12 +43,14 @@ intersectbox1(d1::FullSpace, d2) = d2
 intersectbox2(d1, d2::EmptySpace) = d2
 intersectbox2(d1, d2::FullSpace) = d1
 
-intersectbox(d1::D, d2::D) where {D<:FixedInterval} = d1
+intersectbox1(d1::D, d2::D) where {D<:FixedInterval} = d1
 
-intersectbox(d1::AbstractInterval{T}, d2::AbstractInterval{T}) where {T} =
+intersectbox1(d1::AbstractInterval, d2::AbstractInterval) =
     intersectdomain(d1, d2)
 
-function intersectbox(d1::HyperRectangle{T}, d2::HyperRectangle{T}) where {T}
+function intersectbox1(d1::HyperRectangle, d2::HyperRectangle)
+    @assert dimension(d1) == dimension(d2)
+    T = eltype(d1)
     d = Rectangle{T}(map(intersectbox, components(d1), components(d2)))
     isempty(d) ? emptyspace(d) : d
 end
diff --git a/src/domains/cube.jl b/src/domains/cube.jl
index a34cd63..bace51c 100644
--- a/src/domains/cube.jl
+++ b/src/domains/cube.jl
@@ -351,8 +351,8 @@ ProductDomain(domains::NTuple{N,D}) where {N,D <: FixedInterval} =
 ProductDomain(domains::SVector{N,<:FixedInterval}) where {N} =
 	FixedIntervalProduct(domains)
 ProductDomain{T}(domains::D...) where {N,S,T<:SVector{N,S},D <: FixedInterval} =
-	FixedIntervalProduct(convert.(Ref(Domain{S}), domains))
+	FixedIntervalProduct(convert_eltype.(Ref(S), domains))
 ProductDomain{T}(domains::NTuple{N,D}) where {N,S,T<:SVector{N,S},D <: FixedInterval} =
-	FixedIntervalProduct(convert.(Ref(Domain{S}), domains))
+	FixedIntervalProduct(convert_eltype.(Ref(S), domains))
 ProductDomain{T}(domains::SVector{N,<:FixedInterval}) where {N,T<:SVector{N}} =
 	FixedIntervalProduct(domains)
diff --git a/src/domains/indicator.jl b/src/domains/indicator.jl
index 905c74b..4e631cc 100644
--- a/src/domains/indicator.jl
+++ b/src/domains/indicator.jl
@@ -13,7 +13,7 @@ implementing `in` directly.
 abstract type AbstractIndicatorFunction{T} <: Domain{T} end
 
 "The indicator function of a domain is the function `f(x) = x ∈ D`."
-indicatorfunction(d::Domain) = x -> x ∈ d
+indicatorfunction(d) = x -> x ∈ d
 
 indomain(x, d::AbstractIndicatorFunction) = _indomain(x, d, indicatorfunction(d))
 _indomain(x, d::AbstractIndicatorFunction, f) = f(x)
@@ -35,8 +35,8 @@ indicatorfunction(d::IndicatorFunction) = d.f
 similardomain(d::IndicatorFunction, ::Type{T}) where {T} = IndicatorFunction{T}(d.f)
 
 convert(::Type{IndicatorFunction}, d::AbstractIndicatorFunction) = d
-convert(::Type{IndicatorFunction}, d::Domain{T}) where {T} =
-    IndicatorFunction{T}(indicatorfunction(d))
+convert(::Type{IndicatorFunction}, d) =
+    IndicatorFunction{domaineltype(d)}(indicatorfunction(checkdomain(d)))
 
 isequaldomain(d1::IndicatorFunction, d2::IndicatorFunction) =
     indicatorfunction(d1)==indicatorfunction(d2)
@@ -45,18 +45,19 @@ intersectdomain1(d1::IndicatorFunction, d2) = BoundedIndicatorFunction(d1.f, d2)
 intersectdomain2(d1, d2::IndicatorFunction) = BoundedIndicatorFunction(d2.f, d1)
 
 "An indicator function with a known bounding domain."
-struct BoundedIndicatorFunction{F,D,T} <: AbstractIndicatorFunction{T}
+struct BoundedIndicatorFunction{T,F,D} <: AbstractIndicatorFunction{T}
     f       ::  F
     domain  ::  D
 end
 
-BoundedIndicatorFunction(f, domain::Domain{T}) where {T} =
-    BoundedIndicatorFunction{typeof(f),typeof(domain),T}(f, domain)
-
-function BoundedIndicatorFunction(f, domain)
-    T = eltype(domain)
-    BoundedIndicatorFunction{typeof(f),typeof(domain),T}(f, domain)
-end
+BoundedIndicatorFunction(f, domain) =
+    BoundedIndicatorFunction{domaineltype(domain)}(f, domain)
+BoundedIndicatorFunction{T}(f, domain::Domain{T}) where {T} =
+    BoundedIndicatorFunction{T,typeof(f),typeof(domain)}(f, domain)
+BoundedIndicatorFunction{T}(f, domain) where {T} =
+    _BoundedIndicatorFunction(T, f, checkdomain(domain))
+_BoundedIndicatorFunction(::Type{T}, f, domain) where {T} =
+    BoundedIndicatorFunction{T,typeof(f),typeof(domain)}(f, domain)
 
 indicatorfunction(d::BoundedIndicatorFunction) = d.f
 
@@ -70,7 +71,7 @@ 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))
+    BoundedIndicatorFunction(d.f, convert_eltype(T, d.domain))
 
 boundingbox(d::BoundedIndicatorFunction) = boundingbox(boundingdomain(d))
 
diff --git a/src/domains/interval.jl b/src/domains/interval.jl
index 67bcede..096368d 100644
--- a/src/domains/interval.jl
+++ b/src/domains/interval.jl
@@ -2,7 +2,7 @@
 iscompact(d::TypedEndpointsInterval{:closed,:closed}) = true
 iscompact(d::TypedEndpointsInterval) = false
 
-isinterval(d::Domain) = false
+isinterval(d) = false
 isinterval(d::AbstractInterval) = true
 
 hash(d::AbstractInterval, h::UInt) =
@@ -23,7 +23,7 @@ approx_indomain(x, d::TypedEndpointsInterval{:open,:open}, tolerance) =
     (x > leftendpoint(d)-tolerance) && (x < rightendpoint(d)+tolerance)
 
 
-function point_in_domain(d::AbstractInterval)
+function choice(d::AbstractInterval)
     isempty(d) && throw(BoundsError())
     mean(d)
 end
@@ -31,7 +31,7 @@ end
 center(d::AbstractInterval) = mean(d)
 
 # For an interval of integers, try to find an integer point
-function point_in_domain(d::AbstractInterval{T}) where {T<:Integer}
+function choice(d::AbstractInterval{T}) where {T<:Integer}
     isempty(d) && throw(BoundsError())
     x = round(T, mean(d))
     if x ∈ d
@@ -54,9 +54,9 @@ end
 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))
+normal(d::AbstractInterval, x) = (abs(minimum(d)-x) < abs(maximum(d)-x)) ? -one(domaineltype(d)) : one(domaineltype(d))
 
-distance_to(d::AbstractInterval, x) = x ∈ d ? zero(eltype(d)) : min(abs(x-supremum(d)), abs(x-infimum(d)))
+distance_to(d::AbstractInterval, x) = x ∈ d ? zero(domaineltype(d)) : min(abs(x-supremum(d)), abs(x-infimum(d)))
 
 boundingbox(d::AbstractInterval) = d
 
@@ -79,6 +79,10 @@ and `ChebyshevInterval`.
 abstract type FixedInterval{L,R,T} <: TypedEndpointsInterval{L,R,T} end
 const ClosedFixedInterval{T} = FixedInterval{:closed,:closed,T}
 
+convert(::Type{AbstractInterval{T}}, d::FixedInterval{L,R,T}) where {L,R,T} = d
+convert(::Type{AbstractInterval{T}}, d::FixedInterval{L,R,S}) where {L,R,S,T} =
+    similardomain(d, T)
+
 closure(d::AbstractInterval) = ClosedInterval(endpoints(d)...)
 closure(d::ClosedInterval) = d
 
@@ -244,8 +248,8 @@ function similar_interval(d::HalfLine{T,C}, a::S, b::S) where {T,S,C}
     HalfLine{promote_type(float(T),S),C}()
 end
 
-point_in_domain(d::NonnegativeRealLine) = zero(eltype(d))
-point_in_domain(d::PositiveRealLine) = one(eltype(d))
+choice(d::NonnegativeRealLine) = zero(domaineltype(d))
+choice(d::PositiveRealLine) = one(domaineltype(d))
 
 
 "The negative halfline `(-∞,0]` or `(-∞,0)`, right-closed or right-open."
@@ -276,15 +280,15 @@ function similar_interval(d::NegativeHalfLine{T,C}, a::S, b::S) where {S,T,C}
     NegativeHalfLine{promote_type(S,float(T)),C}()
 end
 
-point_in_domain(d::NegativeRealLine) = -one(eltype(d))
-point_in_domain(d::NonpositiveRealLine) = zero(eltype(d))
+choice(d::NegativeRealLine) = -one(domaineltype(d))
+choice(d::NonpositiveRealLine) = zero(domaineltype(d))
 
 
 "The real line `(-∞,∞)`."
 struct RealLine{T} <: FixedInterval{:open,:open,T} end
 RealLine() = RealLine{Float64}()
 
-point_in_domain(d::RealLine) = zero(eltype(d))
+choice(d::RealLine) = zero(eltype(d))
 
 similardomain(::RealLine, ::Type{T}) where {T} = RealLine{T}()
 
@@ -318,11 +322,11 @@ intersect(d1::FixedInterval, d2::FixedInterval) = intersectdomain(d1, d2)
 
 # Promotion to joint type T
 uniondomain(d1::TypedEndpointsInterval, d2::TypedEndpointsInterval) =
-    uniondomain(promote(d1,d2)...)
+    uniondomain(promote_domains(d1,d2)...)
 intersectdomain(d1::TypedEndpointsInterval, d2::TypedEndpointsInterval) =
-    intersectdomain(promote(d1,d2)...)
+    intersectdomain(promote_domains(d1,d2)...)
 setdiffdomain(d1::AbstractInterval, d2::AbstractInterval) =
-    setdiffdomain(promote(d1,d2)...)
+    setdiffdomain(promote_domains(d1,d2)...)
 
 
 
diff --git a/src/domains/point.jl b/src/domains/point.jl
index 4d2a54b..fed2e6e 100644
--- a/src/domains/point.jl
+++ b/src/domains/point.jl
@@ -48,7 +48,7 @@ supremum(d::Point) = pointval(d)
 interior(d::Point) = emptyspace(d)
 closure(d::Point) = d
 
-point_in_domain(d::Point) = pointval(d)
+choice(d::Point) = pointval(d)
 
 distance_to(d::Point, x) = norm(x-pointval(d))
 
diff --git a/src/domains/simplex.jl b/src/domains/simplex.jl
index ed23f52..059bcbc 100644
--- a/src/domains/simplex.jl
+++ b/src/domains/simplex.jl
@@ -143,7 +143,7 @@ closure(d::VectorUnitSimplex{T}) where {T} = VectorUnitSimplex{T,:closed}(dimens
 
 # We pick the center point, because it belongs to the domain regardless of
 # whether it is open or closed.
-point_in_domain(d::UnitSimplex) = center(d)
+choice(d::UnitSimplex) = center(d)
 
 similardomain(d::StaticUnitSimplex{S,C}, ::Type{T}) where {S,T,C} =
     StaticUnitSimplex{T,C}()
diff --git a/src/domains/trivial.jl b/src/domains/trivial.jl
index 763513b..75283b4 100644
--- a/src/domains/trivial.jl
+++ b/src/domains/trivial.jl
@@ -13,7 +13,7 @@ EmptySpace(::Type{T}) where {T} = EmptySpace{T}()
 similardomain(::EmptySpace, ::Type{T}) where {T} = EmptySpace{T}()
 
 "Return the empty space with the same element type as the given domain."
-emptyspace(d) = emptyspace(eltype(d))
+emptyspace(d) = emptyspace(domaineltype(d))
 emptyspace(::Type{T}) where {T} = EmptySpace{T}()
 
 indomain(x::T, d::EmptySpace{T}) where {T} = false
@@ -64,14 +64,17 @@ struct FullSpace{T} <: Domain{T} end
 const AnyFullSpace = FullSpace{Any}
 
 FullSpace() = FullSpace{Float64}()
-FullSpace(d) = FullSpace{eltype(d)}()
+FullSpace(d) = FullSpace{domaineltype(d)}()
 
 "Return the full space with the same element type as the given domain."
-fullspace(d) = fullspace(eltype(d))
+fullspace(d) = fullspace(domaineltype(d))
 fullspace(::Type{T}) where {T} = FullSpace{T}()
 
 isfullspace(d::FullSpace) = true
 isfullspace(d::Domain) = false
+isfullspace(d) = _isfullspace(d, DomainStyle(d))
+_isfullspace(d, ::IsDomain) = false
+_isfullspace(d, ::NotDomain) = error("isfullspace invoked on non-domain type")
 
 similardomain(::FullSpace, ::Type{T}) where {T} = FullSpace{T}()
 
@@ -84,7 +87,7 @@ approx_indomain(x, d::FullSpace, tolerance) = in(x, d)
 show(io::IO, d::FullSpace) = print(io, "{x} (full space)")
 
 # We choose the origin as a point in the full space
-point_in_domain(d::FullSpace) = zero(eltype(d))
+choice(d::FullSpace) = zero(eltype(d))
 
 isempty(::FullSpace) = false
 
@@ -125,7 +128,7 @@ type `T`.
 struct TypeDomain{T} <: Domain{T} end
 
 "Return the domain for the element type of the given domain."
-typedomain(d) = typedomain(eltype(d))
+typedomain(d) = typedomain(domaineltype(d))
 typedomain(::Type{T}) where {T} = TypeDomain{T}()
 
 iscompatiblepair(x::T, d::TypeDomain{T}) where {T} = true
diff --git a/src/generic/broadcast.jl b/src/generic/broadcast.jl
index d0acc05..96178e5 100644
--- a/src/generic/broadcast.jl
+++ b/src/generic/broadcast.jl
@@ -6,8 +6,10 @@
 struct DomainSetStyle <: Base.Broadcast.BroadcastStyle end
 
 Base.BroadcastStyle(::Type{<:Domain}) = DomainSetStyle()
+Base.BroadcastStyle(::Type{<:AsDomain}) = DomainSetStyle()
 
 Base.broadcastable(d::Domain) = d
+Base.broadcastable(d::AsDomain) = d
 
 # DomainSetStyle doesn't mix with ArrayStyle
 Base.BroadcastStyle(::DomainSetStyle, ::Base.Broadcast.AbstractArrayStyle) = DomainSetStyle()
@@ -15,34 +17,35 @@ Base.BroadcastStyle(::Base.Broadcast.AbstractArrayStyle, ::DomainSetStyle) = Dom
 
 import Base.Broadcast: broadcasted
 
-broadcasted(::DomainSetStyle, ::typeof(+), a::Union{Number,AbstractArray}, d::Domain) =
-    map_domain(Translation(a), d)
-broadcasted(::DomainSetStyle, ::typeof(+), d::Domain, a::Union{Number,AbstractArray}) =
-    map_domain(Translation(a), d)
+broadcasted(::DomainSetStyle, ::typeof(+), a::Union{Number,AbstractArray}, d::AnyDomain) =
+    map_domain(Translation(a), domain(d))
+broadcasted(::DomainSetStyle, ::typeof(+), d::AnyDomain, a::Union{Number,AbstractArray}) =
+    map_domain(Translation(a), domain(d))
 
-broadcasted(::DomainSetStyle, ::typeof(-), a::Union{Number,AbstractArray}, d::Domain) =
-    map_domain(AffineMap(-1, a), d)
-broadcasted(::DomainSetStyle, ::typeof(-), d::Domain, a::Union{Number,AbstractArray}) =
-    map_domain(Translation(-a), d)
-broadcasted(::DomainSetStyle, ::typeof(-), d::Domain{T}) where {T} =
-    map_domain(LinearMap{T}(-1), d)
+broadcasted(::DomainSetStyle, ::typeof(-), a::Union{Number,AbstractArray}, d::AnyDomain) =
+    map_domain(AffineMap(-1, a), domain(d))
+broadcasted(::DomainSetStyle, ::typeof(-), d::AnyDomain, a::Union{Number,AbstractArray}) =
+    map_domain(Translation(-a), domain(d))
+broadcasted(::DomainSetStyle, ::typeof(-), d::AnyDomain) =
+    map_domain(LinearMap{domaineltype(d)}(-1), domain(d))
 
-broadcasted(::DomainSetStyle, ::typeof(*), a::Number, d::Domain{T}) where {T} =
-    map_domain(LinearMap{T}(a), d)
-broadcasted(::DomainSetStyle, ::typeof(*), d::Domain{T}, a::Number) where {T} =
-    map_domain(LinearMap{T}(a), d)
+broadcasted(::DomainSetStyle, ::typeof(*), a::Number, d::AnyDomain) =
+    map_domain(LinearMap{domaineltype(d)}(a), domain(d))
+broadcasted(::DomainSetStyle, ::typeof(*), d::AnyDomain, a::Number) =
+    map_domain(LinearMap{domaineltype(d)}(a), domain(d))
 
 
-broadcasted(::DomainSetStyle, ::typeof(/), d::Domain, a::Number) =
-    mapped_domain(LinearMap(a), d)
+broadcasted(::DomainSetStyle, ::typeof(/), d::AnyDomain, a::Number) =
+    mapped_domain(LinearMap(a), domain(d))
 
-broadcasted(::DomainSetStyle, ::typeof(\), a::Number, d::Domain) =
-    mapped_domain(LinearMap(a), d)
+broadcasted(::DomainSetStyle, ::typeof(\), a::Number, d::AnyDomain) =
+    mapped_domain(LinearMap(a), domain(d))
 
-broadcasted(::DomainSetStyle, m::AbstractMap, d::Domain) = map_domain(m, d)
+broadcasted(::DomainSetStyle, m::AbstractMap, d::AnyDomain) =
+    map_domain(m, domain(d))
 
-broadcasted(::DomainSetStyle, fun::Function, d::Domain{T}) where {T} =
-    convert(Map{T}, fun).(d)
+broadcasted(::DomainSetStyle, fun::Function, d::AnyDomain) =
+    convert(Map{domaineltype(d)}, fun).(d)
 
 # Intercept broadcast applied to `in`, e.g. in.(A, d).
 # This gives domains an opportunity to provide a more efficient implementation
@@ -50,31 +53,34 @@ broadcasted(::DomainSetStyle, fun::Function, d::Domain{T}) where {T} =
 # has particular structure. A common case would be a raster of points,
 # for plotting purposes.
 # This call can be avoided by typing in.(A, Ref(d)) instead.
-broadcasted(::DomainSetStyle, ::typeof(in), A, d::Domain) = vectorized_in(A, d)
-broadcasted(::DomainSetStyle, ::typeof(approx_in), A, d::Domain) = vectorized_approx_in(A, d)
-broadcasted(::DomainSetStyle, ::typeof(approx_in), A, d::Domain, tol) = vectorized_approx_in(A, d, tol)
+broadcasted(::DomainSetStyle, ::typeof(in), A, d::AnyDomain) =
+    vectorized_in(A, domain(d))
+broadcasted(::DomainSetStyle, ::typeof(approx_in), A, d::AnyDomain) =
+    vectorized_approx_in(A, domain(d))
+broadcasted(::DomainSetStyle, ::typeof(approx_in), A, d::AnyDomain, tol) =
+    vectorized_approx_in(A, domain(d), tol)
 
-broadcasted(::DomainSetStyle, ::typeof(∉), A, d::Domain) = A .∉ Ref(d)
+broadcasted(::DomainSetStyle, ::typeof(∉), A, d::AnyDomain) = A .∉ Ref(domain(d))
 
 @deprecate broadcast_in(A, d::Domain) vectorized_in(A, d)
 @deprecate broadcast_approx_in(A, d::Domain) vectorized_approx_in(A, d)
 @deprecate broadcast_approx_in(A, d::Domain, tol) vectorized_approx_in(A, d, tol)
 
 "Vectorized version of `in`: apply `x ∈ d` to all elements of `A`."
-vectorized_in(A, d::Domain) = in.(A, Ref(d))
+vectorized_in(A, d) = in.(A, Ref(d))
 
 "Vectorized version of `approx_in`: apply `x ∈ d` to all elements of `A`."
-vectorized_approx_in(A, d::Domain) = approx_in.(A, Ref(d))
-vectorized_approx_in(A, d::Domain, tol) = approx_in.(A, Ref(d), tol)
+vectorized_approx_in(A, d) = approx_in.(A, Ref(d))
+vectorized_approx_in(A, d, tol) = approx_in.(A, Ref(d), tol)
 
 
 ## Some arithmetics
 
 # Allow unary minus, but use broadcast for the implementation
--(d::Domain) = (-).(d)
+-(d::AnyDomain) = (-).(d)
 
 # Allow multiplication and division by numbers, like for vectors
-*(a::Number, domain::Domain) = a .* domain
-*(domain::Domain, a::Number) = domain .* a
-/(domain::Domain, a::Number) = domain ./ a
-\(a::Number, domain::Domain) = a .\ domain
+*(a::Number, domain::AnyDomain) = a .* domain
+*(domain::AnyDomain, a::Number) = domain .* a
+/(domain::AnyDomain, a::Number) = domain ./ a
+\(a::Number, domain::AnyDomain) = a .\ domain
diff --git a/src/generic/canonical.jl b/src/generic/canonical.jl
index 7787b13..6b3bb82 100644
--- a/src/generic/canonical.jl
+++ b/src/generic/canonical.jl
@@ -1,6 +1,6 @@
 
 """
-    canonicaldomain([ctype::CanonicalType, ]d::Domain)
+    canonicaldomain([ctype::CanonicalType, ]domain)
 
 Return an associated canonical domain, if any, of the given domain.
 
@@ -18,7 +18,7 @@ canonicaldomain(d) = d
 "Does the domain have a canonical domain?"
 hascanonicaldomain(d) = !(d === canonicaldomain(d))
 
-identitymap(d) = IdentityMap{eltype(d)}(dimension(d))
+identitymap(d) = IdentityMap{domaineltype(d)}(dimension(d))
 
 """
     mapfrom_canonical(d[, x])
@@ -57,6 +57,11 @@ mapto_canonical(::Equal, d) = leftinverse(mapfrom_canonical(Equal(), d))
 simplify(d) = canonicaldomain(Equal(), d)
 simplifies(d) = hascanonicaldomain(Equal(), d)
 
+"Convert the given domain to a domain defined in DomainSets.jl."
+todomainset(d::Domain) = d
+todomainset(d::AsDomain) = todomainset(domain(d))
+
+
 "A canonical domain that is isomorphic but may have different element type."
 struct Isomorphic <: CanonicalType end
 
@@ -75,6 +80,7 @@ mapfrom_canonical(::Isomorphic, d::Domain{NTuple{N,T}}) where {N,T} =
     VectorToTuple{N,T}()
 
 
+
 "A parameter domain that can be mapped to the domain."
 struct Parameterization <: CanonicalType end
 
@@ -119,7 +125,6 @@ no_known_mapto(d1, d2) = d1 == d2 ? identitymap(d1) : error("No map known betwee
 # We generically define the equality of Domains, using the framework
 # of canonical domains.
 
-
 @deprecate isequal1(d1,d2) isequaldomain1(d1,d2)
 @deprecate isequal2(d1,d2) isequaldomain2(d1,d2)
 
@@ -130,4 +135,6 @@ isequaldomain1(d1, d2) = simplifies(d1) ? simplify(d1)==d2 : isequaldomain2(d1,
 # simplify the second argument
 isequaldomain2(d1, d2) = simplifies(d2) ? d1==simplify(d2) : d1===d2
 
-==(d1::Domain, d2::Domain) = isequaldomain(d1, d2)
+==(d1::AnyDomain, d2::AnyDomain) = isequaldomain(domain(d1), domain(d2))
+==(d1::AnyDomain, d2) = isequaldomain(domain(d1), d2)
+==(d1, d2::AnyDomain) = isequaldomain(d1, domain(d2))
diff --git a/src/generic/core.jl b/src/generic/core.jl
new file mode 100644
index 0000000..682c8c4
--- /dev/null
+++ b/src/generic/core.jl
@@ -0,0 +1,113 @@
+export Domain,
+    domain,
+    domaineltype,
+    DomainStyle,
+    IsDomain,
+    NotDomain,
+    AsDomain,
+    AnyDomain,
+    checkdomain
+
+# """
+# A domain is a set of elements that is possibly continuous.
+#
+# Examples may be intervals and triangles. These are geometrical shapes, but more
+# generally a domain can be any type that supports `in`. Conceptually, a domain
+# is the set of all elements `x` for which `in(x, domain)` returns true.
+#
+# A `Domain{T}` is a domain with elements of type `T`, in analogy with
+# `AbstractSet{T}` and `AbstractVector{T}`. Although domains may be defined by a
+# mathematical condition such as `a <= x <= b`, irrespective of the type of `x`,
+# points generated to belong to the domain have type `T`.
+# """
+# abstract type Domain{T} end
+
+import IntervalSets: Domain
+Base.eltype(::Type{<:Domain{T}}) where {T} = T
+
+"""
+    domaineltype(d)
+
+The `domaineltype` of a continuous domain is the element type of any
+discretization of that domain.
+"""
+domaineltype(d) = domaineltype(typeof(d))
+domaineltype(::Type{D}) where D = eltype(D)
+domaineltype(d::Domain{T}) where T = T      # shortcut definition
+
+
+abstract type DomainStyle end
+
+"""
+    IsDomain()
+
+indicates an object implements the domain interface.
+"""
+struct IsDomain <: DomainStyle end
+
+"""
+    NotDomain()
+
+indicates an object does not implement the domain interface.
+"""
+struct NotDomain <: DomainStyle end
+
+
+"""
+    DomainStyle(d)
+
+The domain style of `d` is a trait to indicate whether or not `d` implements
+the domain interface. If so, the object supports `in(x,d)` and optionally also
+implements `domaineltype(d)`.
+"""
+DomainStyle(d) = DomainStyle(typeof(d))
+DomainStyle(::Type) = NotDomain()
+DomainStyle(::Type{<:Domain}) = IsDomain()
+DomainStyle(::Type{<:Number}) = IsDomain()
+DomainStyle(::Type{<:AbstractSet}) = IsDomain()
+DomainStyle(::Type{<:AbstractArray}) = IsDomain()
+
+"""
+A reference to a domain.
+
+In a function call, `AsDomain(x)` can be used to indicate that `x` should be
+treated as a domain in the function, e.g., `foo(x, AsDomain(d))`.
+"""
+abstract type AsDomain end
+
+"A reference to a specific given domain."
+struct DomainRef{D} <: AsDomain
+    domain  ::  D
+end
+domain(d::DomainRef) = d.domain
+domain(d::Domain) = d
+
+Base.eltype(::Type{<:DomainRef{D}}) where D = domaineltype(D)
+
+domaineltype(d::AsDomain) = domaineltype(domain(d))
+domaineltype(::Type{<:DomainRef{D}}) where D = domaineltype(D)
+
+AsDomain(d) = DomainRef(d)
+AsDomain(d::Domain) = d
+
+"""
+`AnyDomain` can be a concrete domain or a reference to a domain.
+
+In both cases `domain(d::AnyDomain)` returns the domain itself.
+"""
+const AnyDomain = Union{Domain,AsDomain}
+
+"""
+   checkdomain(d)
+
+Checks that `d` is a domain or refers to a domain and if so returns that domain,
+throws an error otherwise.
+"""
+checkdomain(d::Domain) = d
+# we trust the explicit intention of a user providing a domain reference
+checkdomain(d::AsDomain) = domain(d)
+# for other objects we check DomainStyle
+checkdomain(d) = _checkdomain(d, DomainStyle(d))
+_checkdomain(d, ::IsDomain) = d
+_checkdomain(d, ::NotDomain) =
+    error("Domain does not implement domain interface as indicated by DomainStyle.")
diff --git a/src/generic/domain.jl b/src/generic/domain.jl
index a9c1bcd..001162f 100644
--- a/src/generic/domain.jl
+++ b/src/generic/domain.jl
@@ -1,15 +1,15 @@
 # Definition of the abstract Domain type and its interface
 
-# The type Domain{T} is defined in IntervalSets.jl
-
-eltype(::Type{<:Domain{T}}) where {T} = T
-prectype(::Type{<:Domain{T}}) where {T} = prectype(T)
-numtype(::Type{<:Domain{T}}) where {T} = numtype(T)
+Domain(d) = convert(Domain, d)
 
-convert_numtype(d::Domain{T}, ::Type{U}) where {T,U} = convert(Domain{to_numtype(T,U)}, d)
-convert_prectype(d::Domain{T}, ::Type{U}) where {T,U} = convert(Domain{to_prectype(T,U)}, d)
+prectype(::Type{<:Domain{T}}) where T = prectype(T)
+numtype(::Type{<:Domain{T}}) where T = numtype(T)
 
-Domain(d) = convert(Domain, d)
+# Domain-specific prectype and numtype default to using domaineltype
+domain_prectype(d) = prectype(domaineltype(d))
+domain_numtype(d) = numtype(domaineltype(d))
+prectype(d::AsDomain) = prectype(domaineltype(d))
+numtype(d::AsDomain) = numtype(domaineltype(d))
 
 # Concrete types can implement similardomain(d, ::Type{T}) where {T}
 # to support convert(Domain{T}, d) functionality.
@@ -17,7 +17,7 @@ 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(domains...) = promotable_eltypes(map(eltype, domains)...)
+promotable_domains(domains...) = promotable_eltypes(map(domaineltype, domains)...)
 promotable_eltypes(types...) = isconcretetype(promote_type(types...))
 promotable_eltypes(::Type{S}, ::Type{T}) where {S<:AbstractVector,T<:AbstractVector} =
     promotable_eltypes(eltype(S), eltype(T))
@@ -25,20 +25,17 @@ promotable_eltypes(::Type{S}, ::Type{T}) where {S<:AbstractVector,T<:AbstractVec
 "Promote the given domains to have a common element type."
 promote_domains() = ()
 promote_domains(domains...) = promote_domains(domains)
-promote_domains(domains) = convert_eltype.(mapreduce(eltype, promote_type, domains), domains)
-
-promote_domains(domains::AbstractSet{<:Domain{T}}) where {T} = domains
-promote_domains(domains::AbstractSet{<:Domain}) = Set(promote_domains(collect(domains)))
+promote_domains(domains) = convert_eltype.(mapreduce(domaineltype, promote_type, domains), domains)
 
 convert_eltype(::Type{T}, d::Domain) where {T} = convert(Domain{T}, d)
-convert_eltype(::Type{T}, d) where {T} = _convert_eltype(T, d, eltype(d))
-_convert_eltype(::Type{T}, d, ::Type{T}) where {T} = d
-_convert_eltype(::Type{T}, d, ::Type{S}) where {S,T} =
+convert_eltype(::Type{T}, d) where {T} = _convert_eltype(T, d, domaineltype(d), DomainStyle(d))
+_convert_eltype(::Type{T}, d, ::Type{T}, ::IsDomain) where {T} = d
+_convert_eltype(::Type{T}, d, ::Type{S}, ::IsDomain) where {S,T} =
     error("Don't know how to convert the `eltype` of $(d).")
+_convert_eltype(::Type{T}, d, ::Type{T}, ::NotDomain) where {T} = d
+_convert_eltype(::Type{T}, d, ::Type{S}, ::NotDomain) where {S,T} =
+    error("Convert eltype: argument given is not a domain.")
 
-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))
 
 "A `EuclideanDomain` is any domain whose eltype is `<:StaticVector{N,T}`."
 const EuclideanDomain{N,T} = Domain{<:StaticVector{N,T}}
@@ -52,10 +49,10 @@ const AbstractVectorDomain{T} = Domain{<:AbstractVector{T}}
 CompositeTypes.Display.displaysymbol(d::Domain) = 'D'
 
 "What is the Euclidean dimension of the domain?"
-dimension(::Domain{T}) where {T} = euclideandimension(T)
+dimension(d) = euclideandimension(domaineltype(d))
 
 "Is the given combination of point and domain compatible?"
-iscompatiblepair(x, d) = _iscompatiblepair(x, d, typeof(x), eltype(d))
+iscompatiblepair(x, d) = _iscompatiblepair(x, d, typeof(x), domaineltype(d))
 _iscompatiblepair(x, d, ::Type{S}, ::Type{T}) where {S,T} =
     _iscompatiblepair(x, d, S, T, promote_type(S,T))
 _iscompatiblepair(x, d, ::Type{S}, ::Type{T}, ::Type{U}) where {S,T,U} = true
@@ -70,15 +67,15 @@ iscompatiblepair(x::AbstractVector, ::EuclideanDomain{N}) where {N} = length(x)=
 # Note: there are cases where this warning reveals a bug, and cases where it is
 # annoying. In cases where it is annoying, the domain may want to specialize `in`.
 compatible_or_false(x, domain) =
-    iscompatiblepair(x, domain) ? true : (@warn "`in`: incompatible combination of point: $(typeof(x)) and domain eltype: $(eltype(domain)). Returning false."; false)
+    iscompatiblepair(x, domain) ? true : (@warn "`in`: incompatible combination of point: $(typeof(x)) and domain eltype: $(domaineltype(domain)). Returning false."; false)
 
 compatible_or_false(x::AbstractVector, domain::AbstractVectorDomain) =
     iscompatiblepair(x, domain) ? true : (@warn "`in`: incompatible combination of vector with length $(length(x)) and domain '$(domain)' with dimension $(dimension(domain)). Returning false."; false)
 
 
 "Promote point and domain to compatible types."
-promote_pair(x, d) = _promote_pair(x, d, promote_type(typeof(x),eltype(d)))
-_promote_pair(x, d, ::Type{T}) where {T} = convert(T, x), convert(Domain{T}, d)
+promote_pair(x, d) = _promote_pair(x, d, promote_type(typeof(x), domaineltype(d)))
+_promote_pair(x, d, ::Type{T}) where {T} = convert(T, x), convert_eltype(T, d)
 _promote_pair(x, d, ::Type{Any}) = x, d
 # Some exceptions:
 # - matching types: avoid promotion just in case it is expensive
@@ -108,9 +105,14 @@ Return a suitable tolerance to use for verifying whether a point is close to
 a domain. Typically, the tolerance is close to the precision limit of the numeric
 type associated with the domain.
 """
-default_tolerance(d::Domain) = default_tolerance(prectype(d))
-default_tolerance(::Type{T}) where {T <: AbstractFloat} = 100eps(T)
+domain_tolerance(d) = domain_tolerance(domain_prectype(d))
+domain_tolerance(::Type{T}) where {T <: AbstractFloat} = 100eps(T)
 
+# a version with a tolerance, for use in approx_in
+function compatible_or_false(x, d, tol)
+    tol >= 0 || error("Tolerance has to be positive in `approx_in`.")
+    compatible_or_false(x, d)
+end
 
 """
 `approx_in(x, domain::Domain [, tolerance])`
@@ -127,19 +129,21 @@ Up to inexact computations due to floating point numbers, it should also be
 the case that `approx_in(x, d, 0) == in(x,d)`. This implies that `approx_in`
 reflects whether a domain is open or closed.
 """
-approx_in(x, d::Domain) = approx_in(x, d, default_tolerance(d))
-
-function compatible_or_false(x, d, tol)
-    tol >= 0 || error("Tolerance has to be positive in `approx_in`.")
-    compatible_or_false(x, d)
-end
-
+approx_in(x, d) = approx_in(x, d, domain_tolerance(d))
 approx_in(x, d::Domain, tol) =
     compatible_or_false(x, d, tol) && approx_indomain(promote_pair(x, d)..., tol)
+approx_in(x, d::AsDomain, tol) =
+    compatible_or_false(x, domain(d), tol) && approx_indomain(promote_pair(x, domain(d))..., tol)
 
 # Fallback to `in`
-approx_indomain(x, d::Domain, tol) = in(x, d)
+approx_indomain(x, d, tol) = in(x, d)
 
 isapprox(d1::Domain, d2::Domain; kwds...) = d1 == d2
 
 isreal(d::Domain) = isreal(eltype(d))
+
+"Return a point from the given domain."
+function choice(d) end
+
+choice(d::AbstractSet) = first(d)
+choice(d::AbstractArray) = first(d)
diff --git a/src/generic/generator.jl b/src/generic/generator.jl
index 1c103d9..d95096b 100644
--- a/src/generic/generator.jl
+++ b/src/generic/generator.jl
@@ -5,9 +5,9 @@ Domain(gen::Base.Generator) = generator_domain(gen)
 generator_domain(gen::Base.Generator) = generator_domain(gen.f, gen.iter)
 
 # Example: Domain(x for x in 0..1)
-generator_domain(f::typeof(identity), iter::Domain) = iter
+generator_domain(f::typeof(identity), iter) = checkdomain(iter)
 # Example: Domain(x>1 for x in 0..2)
-generator_domain(f, iter::Domain) = BoundedIndicatorFunction(f, iter)
+generator_domain(f, iter) = BoundedIndicatorFunction(f, checkdomain(iter))
 
 # Example: Domain(x for x in 0..2 if x > 1)
 generator_domain(f::typeof(identity), iter::Base.Iterators.Filter) =
diff --git a/src/generic/lazy.jl b/src/generic/lazy.jl
index 0afd0b9..0c05c1e 100644
--- a/src/generic/lazy.jl
+++ b/src/generic/lazy.jl
@@ -81,8 +81,8 @@ _approx_indomain(x, d, tolerance, ::Combination, domains) =
 _approx_indomain(x, d, tolerance, ::Product, domains) =
     mapreduce((u,v)->approx_in(u, v, tolerance), &, x, domains)
 
-point_in_domain(d::SimpleLazyDomain) = toexternalpoint(d, point_in_domain(superdomain(d)))
-point_in_domain(d::CompositeDomain) = toexternalpoint(d, map(point_in_domain, components(d)))
+choice(d::SimpleLazyDomain) = toexternalpoint(d, choice(superdomain(d)))
+choice(d::CompositeDomain) = toexternalpoint(d, map(choice, components(d)))
 
 isequaldomain(a::D, b::D) where {D<:CompositeDomain} = components(a) == components(b)
 
@@ -128,11 +128,12 @@ struct WrappedDomain{T,D} <: DerivedDomain{T}
 end
 
 WrappedDomain(domain::Domain{T}) where {T} = WrappedDomain{T}(domain)
-WrappedDomain(domain) = WrappedDomain{eltype(domain)}(domain)
+WrappedDomain(domain) = WrappedDomain{domaineltype(domain)}(domain)
 
 WrappedDomain{T}(domain::D) where {T,D<:Domain{T}} = WrappedDomain{T,D}(domain)
 WrappedDomain{T}(domain::Domain) where {T} = WrappedDomain{T}(convert(Domain{T}, domain))
-WrappedDomain{T}(domain) where {T} = WrappedDomain{T,typeof(domain)}(domain)
+WrappedDomain{T}(domain) where {T} = _WrappedDomain(convert_eltype(T, domain))
+_WrappedDomain(domain) = WrappedDomain{domaineltype(domain),typeof(domain)}(domain)
 
 similardomain(d::WrappedDomain, ::Type{T}) where {T} = WrappedDomain{T}(d.domain)
 
diff --git a/src/generic/mapped.jl b/src/generic/mapped.jl
index 81d5e7d..7a28177 100644
--- a/src/generic/mapped.jl
+++ b/src/generic/mapped.jl
@@ -52,18 +52,23 @@ struct MappedDomain{T,F,D} <: AbstractMappedDomain{T}
 end
 
 # In the constructor, we have to decide which T to use for the MappedDomain.
-# - we don't know anything about invmap: deduce T from the given domain
-MappedDomain(invmap, domain::Domain{T}) where {T} = MappedDomain{T}(invmap, domain)
+# - if we don't know anything about invmap: deduce T from the given domain
+MappedDomain(invmap, domain) =
+    MappedDomain{domaineltype(domain)}(invmap, domain)
 # - if the map is a Map{T}, use that T for the MappedDomain
-MappedDomain(invmap::Map{T}, domain::Domain{S}) where {S,T} = MappedDomain{T}(invmap, domain)
-# - if T is given in the constructor, by all means we use that
-MappedDomain{T}(invmap, domain::Domain) where {T} =
-    MappedDomain{T,typeof(invmap),typeof(domain)}(invmap, domain)
+MappedDomain(invmap::Map{T}, domain) where {T} =
+    MappedDomain{T}(invmap, domain)
+# If T is given in the constructor, by all means we use that:
+MappedDomain{T}(invmap, domain) where {T} =
+    _MappedDomain(T, invmap, checkdomain(domain))
 # - in that case, if the map is a Map{S}, make sure that S matches T
-MappedDomain{T}(invmap::Map{T}, domain::Domain) where {T} =
-    MappedDomain{T,typeof(invmap),typeof(domain)}(invmap, domain)
-MappedDomain{T}(invmap::Map{S}, domain::Domain) where {S,T} =
+MappedDomain{T}(invmap::Map{T}, domain) where {T} =
+    _MappedDomain(T, invmap, checkdomain(domain))
+MappedDomain{T}(invmap::Map{S}, domain) where {S,T} =
     MappedDomain{T}(convert(Map{T}, invmap), domain)
+# invoke the constructor
+_MappedDomain(::Type{T}, invmap, domain) where {T} =
+    MappedDomain{T,typeof(invmap),typeof(domain)}(invmap, domain)
 
 similardomain(d::MappedDomain, ::Type{T}) where {T} =
     MappedDomain{T}(d.invmap, d.domain)
@@ -75,7 +80,7 @@ inverse_map(d::MappedDomain) = d.invmap
 inverse_map(d::MappedDomain, y) = d.invmap(y)
 
 "Map a domain with the inverse of the given map"
-map_domain(map, domain::Domain) = _map_domain(map, domain)
+map_domain(map, domain) = _map_domain(map, checkdomain(domain))
 
 # Fallback: we don't know anything about map, just try to invert
 _map_domain(map, domain) = mapped_domain(inverse(map), domain)
@@ -83,11 +88,11 @@ _map_domain(map::Map{T}, domain::Domain{T}) where {T} =
     mapped_domain(inverse(map), domain)
 # If map is a Map{T}, then verify and if necessary update T
 function _map_domain(map::Map, domain)
-    U = codomaintype(map, eltype(domain))
+    U = codomaintype(map, domaineltype(domain))
     if U == Union{}
         error("incompatible types of $(map) and $(domain)")
     end
-    mapped_domain(inverse(convert(Map{U}, map)), convert(Domain{U}, domain))
+    mapped_domain(inverse(convert(Map{U}, map)), convert_eltype(U, domain))
 end
 
 isequaldomain(a::MappedDomain, b::MappedDomain) =
@@ -95,7 +100,7 @@ isequaldomain(a::MappedDomain, b::MappedDomain) =
 
 
 "Make a mapped domain with the given inverse map"
-mapped_domain(invmap, domain::Domain) = _mapped_domain(invmap, domain)
+mapped_domain(invmap, domain) = _mapped_domain(invmap, checkdomain(domain))
 
 # We face the same task as in the constructor: attempt to identify T
 # Here, we are more flexible, and attempt to do more conversions. We assume
@@ -107,11 +112,11 @@ _mapped_domain(invmap, domain) = MappedDomain(invmap, domain)
 # - invmap is a Map{T}: its codomaintype should match the eltype of the domain
 # -- first, update the numtype
 _mapped_domain(invmap::Map, domain) =
-    _mapped_domain(invmap, domain, promote_type(numtype(invmap),numtype(domain)))
-_mapped_domain(invmap::Map{T}, domain::Domain{S}, ::Type{U}) where {S,T,U} =
+    _mapped_domain(invmap, domain, promote_type(numtype(invmap),domain_numtype(domain)))
+_mapped_domain(invmap, domain, ::Type{U}) where {U} =
     _mapped_domain2(convert_numtype(invmap,U), convert_numtype(domain,U))
 # -- then, ensure the codomaintype of the map equals the element type of the domain
-_mapped_domain2(invmap, domain) = _mapped_domain2(invmap, domain, codomaintype(invmap), eltype(domain))
+_mapped_domain2(invmap, domain) = _mapped_domain2(invmap, domain, codomaintype(invmap), domaineltype(domain))
 # --- it's okay
 _mapped_domain2(invmap, domain, ::Type{T}, ::Type{T}) where {T} =
     MappedDomain(invmap, domain)
@@ -150,9 +155,14 @@ struct ParametricDomain{T,F,D} <: AbstractMappedDomain{T}
     domain  ::  D
 end
 
-ParametricDomain(fmap, domain::Domain) = ParametricDomain{codomaintype(fmap)}(fmap, domain)
+ParametricDomain(fmap, domain) =
+    ParametricDomain{codomaintype(fmap)}(fmap, domain)
 ParametricDomain{T}(fmap, domain::Domain) where {T} =
     ParametricDomain{T,typeof(fmap),typeof(domain)}(fmap, domain)
+ParametricDomain{T}(fmap, domain) where {T} =
+    _ParametricDomain(T, fmap, checkdomain(domain))
+_ParametricDomain(::Type{T}, fmap, domain) where {T} =
+    ParametricDomain{T,typeof(fmap),typeof(domain)}(fmap, domain)
 
 similardomain(d::ParametricDomain, ::Type{T}) where {T} =
     ParametricDomain{T}(d.fmap, d.domain)
@@ -174,7 +184,7 @@ isequaldomain(d1::ParametricDomain, d2::ParametricDomain) =
 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(fmap, domain)
 parametric_domain(fmap, domain::ParametricDomain) =
     parametric_domain(fmap ∘ forward_map(domain), superdomain(domain))
 
diff --git a/src/generic/productdomain.jl b/src/generic/productdomain.jl
index 0c09cce..821cfad 100644
--- a/src/generic/productdomain.jl
+++ b/src/generic/productdomain.jl
@@ -8,15 +8,15 @@ components(d::ProductDomain) = d.domains
 factors(d::ProductDomain) = components(d)
 
 isequaldomain(d1::ProductDomain, d2::ProductDomain) =
-	mapreduce(==, &, components(d1), components(d2))
+	compatibleproductdims(d1,d2) && 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))
 isopenset(d::ProductDomain) = all(isopenset, components(d))
 
-issubset(d1::ProductDomain, d2::ProductDomain) =
-	compatibleproductdims(d1, d2) && all(map(issubset, components(d1), components(d2)))
+issubset_domain(d1::ProductDomain, d2::ProductDomain) =
+	compatibleproductdims(d1, d2) && all(map(issubset, factors(d1), factors(d2)))
 
 volume(d::ProductDomain) = prod(map(volume, components(d)))
 
@@ -24,7 +24,7 @@ distance_to(d::ProductDomain, x) = sqrt(sum(distance_to(component(d, i), x[i])^2
 
 compatibleproductdims(d1::ProductDomain, d2::ProductDomain) =
 	dimension(d1) == dimension(d2) &&
-		all(map(==, map(dimension, components(d1)), map(dimension, components(d2))))
+		all(map(==, map(dimension, factors(d1)), map(dimension, factors(d2))))
 
 Display.combinationsymbol(d::ProductDomain) = Display.Times()
 Display.displaystencil(d::ProductDomain) = composite_displaystencil(d)
@@ -54,10 +54,11 @@ closure(d::ProductDomain) = ProductDomain(map(closure, components(d)))
 center(d::ProductDomain) = toexternalpoint(d, map(center, components(d)))
 
 VcatDomainElement = Union{Domain{<:Number},EuclideanDomain}
+VcatEltype = Union{Type{<:Number},Type{<:SVector}}
 
-ProductDomain(domains...) = _ProductDomain(map(Domain, domains)...)
-_ProductDomain(domains...) = TupleProductDomain(domains...)
-_ProductDomain(domains::VcatDomainElement...) = VcatDomain(domains...)
+ProductDomain(domains...) = _ProductDomain(domains, map(eltype, domains)...)
+_ProductDomain(domains, types...) = TupleProductDomain(domains...)
+_ProductDomain(domains, types::VcatEltype...) = VcatDomain(domains...)
 ProductDomain(domains::AbstractVector) = VectorProductDomain(domains)
 # To create a tuple product domain, invoke ProductDomain{T}. Here, we splat
 # and this may end up creating a VcatDomain instead.
@@ -85,11 +86,12 @@ productdomain2(d1, d2::ProductDomain) = ProductDomain(d1, factors(d2)...)
 
 # Only override cross for variables of type Domain, it may have a different
 # meaning for other variables (like the vector cross product)
-cross(x::Domain...) = productdomain(x...)
+cross(d::Domain...) = productdomain(d...)
 
 ^(d::Domain, n::Int) = productdomain(ntuple(i->d, n)...)
 
-similardomain(d::ProductDomain, ::Type{T}) where {T} = ProductDomain{T}(components(d))
+similardomain(d::ProductDomain, ::Type{T}) where {T} =
+	ProductDomain{T}(factors(d))
 
 canonicaldomain(d::ProductDomain) = any(map(hascanonicaldomain, factors(d))) ?
 	ProductDomain(map(canonicaldomain, components(d))) : d
@@ -150,26 +152,26 @@ struct VectorProductDomain{V<:AbstractVector,DD<:AbstractVector} <: ProductDomai
 	domains	::	DD
 
 	function VectorProductDomain{V,DD}(domains::DD) where {V,DD}
-		@assert eltype(eltype(domains)) == eltype(V)
+		@assert eltype(domaineltype(domains[1])) == eltype(V)
 		new(domains)
 	end
 end
 
 VectorProductDomain(domains::AbstractVector) =
-	VectorProductDomain{Vector{eltype(eltype(domains))}}(domains)
+	VectorProductDomain{Vector{eltype(domaineltype(domains[1]))}}(domains)
 
 VectorProductDomain{V}(domains::AbstractVector{<:Domain{T}}) where {T,V<:AbstractVector{T}} =
 	VectorProductDomain{V,typeof(domains)}(domains)
 function VectorProductDomain{V}(domains::AbstractVector) where {T,V<:AbstractVector{T}}
-	Tdomains = convert.(Domain{T}, domains)
-	VectorProductDomain{V}(Tdomains)
+	Tdomains = convert_eltype.(T, domains)
+	VectorProductDomain{V,typeof(Tdomains)}(Tdomains)
 end
 
 # Convenience: allow constructor to be called with multiple arguments, or with
 # a container that is not a vector
-VectorProductDomain(domains::Domain...) = VectorProductDomain(domains)
+VectorProductDomain(domains...) = VectorProductDomain(promote_domains(domains...))
 VectorProductDomain(domains) = VectorProductDomain(collect(domains))
-VectorProductDomain{V}(domains::Domain...) where {V} = VectorProductDomain{V}(domains)
+VectorProductDomain{V}(domains...) where {V} = VectorProductDomain{V}(promote_domains(domains...))
 VectorProductDomain{V}(domains) where {V} = VectorProductDomain{V}(collect(domains))
 
 # the dimension equals the number of composite elements
@@ -193,8 +195,7 @@ struct TupleProductDomain{T,DD} <: ProductDomain{T}
 end
 
 TupleProductDomain(domains::Vector) = TupleProductDomain(domains...)
-TupleProductDomain(domains::Domain...) = TupleProductDomain(domains)
-TupleProductDomain(domains...) = TupleProductDomain(map(Domain, domains)...)
+TupleProductDomain(domains...) = TupleProductDomain(map(checkdomain, domains))
 function TupleProductDomain(domains::Tuple)
 	T = Tuple{map(eltype, domains)...}
 	TupleProductDomain{T}(domains)
@@ -203,7 +204,7 @@ end
 TupleProductDomain{T}(domains::Vector) where {T} = TupleProductDomain{T}(domains...)
 TupleProductDomain{T}(domains...) where {T} = TupleProductDomain{T}(domains)
 function TupleProductDomain{T}(domains::Tuple) where {T <: Tuple}
-	Tdomains = map((t,d) -> convert(Domain{t},d), tuple(T.parameters...), domains)
+	Tdomains = map((t,d) -> convert_eltype(t, d), tuple(T.parameters...), domains)
 	TupleProductDomain{T,typeof(Tdomains)}(Tdomains)
 end
 TupleProductDomain{T}(domains::Tuple) where {T} =
diff --git a/src/generic/setoperations.jl b/src/generic/setoperations.jl
index f64abb7..d86e763 100644
--- a/src/generic/setoperations.jl
+++ b/src/generic/setoperations.jl
@@ -1,21 +1,24 @@
 # The union, intersection and difference of domains are represented with lazy domains.
 
-issubset(d1::Domain, d2::Domain) = promotable_domains(d1, d2) && issubset1(promote_domains(d1, d2)...)
-issubset(d1, d2::Domain) = issubset1(promote_domains(d1, d2)...)
-issubset1(d1, d2) = issubset2(d1, d2)
-issubset2(d1, d2) = d1 == d2
+issubset(d1::AnyDomain, d2::AnyDomain) = issubset_domain(domain(d1), domain(d2))
+issubset(d1::AnyDomain, d2) = issubset_domain(domain(d1), d2)
+
+issubset_domain(d1, d2) =
+	promotable_domains(d1, d2) && issubset1(promote_domains(d1, d2)...)
+
+issubset1(d1, d2) = simplifies(d1) ? issubset(simplify(d1),d2) : issubset2(d1, d2)
+issubset2(d1, d2) = simplifies(d2) ? issubset(d1, simplify(d2)) : d1 == d2
 # this last fallback is only an approximation of the truth: if d1 equals d2, then
 # d1 is a subset of d2, but the reverse is not true. So we might be returning false
 # even when the correct mathematical answer is true.
 # What `issubset` means for the optimizations below is:
 # - if true: we are sure that d1 is a subset of d2
 # - if false: either it really is false, or we don't really know
-# On top of that, we only invoke issubset(d1,d2) when the arguments are both
-# Domains, or when d2 is a Domain, anticipating that issubset(d1::Domain,d2)
-# is often harder to implement.
 
-issubset(d1::AbstractArray, d2::Domain) = all(in(d2), d1)
-issubset(d1::AbstractSet, d2::Domain) = all(in(d2), d1)
+issubset1(d1::Number, d2) = in(d1, d2)
+issubset1(d1::AbstractArray, d2) = all(in(d2), d1)
+issubset1(d1::AbstractSet, d2) = all(in(d2), d1)
+
 
 ############################
 # The union of two domains
@@ -31,16 +34,15 @@ end
 """
 The `UnionDomain` and `UnionDomain{T}` constructors can be invoked in three ways:
 - with a list of arguments: `UnionDomain(d1, d2, ...)`
-- with a single domain: `UnionDomain(d::Domain)`
-- or with any iterable list of domains: `UnionDomain(domains)`
+- or with an iterable list of domains: `UnionDomain(domains)`
 """
 UnionDomain(domains...) = UnionDomain(domains)
-UnionDomain(d::Domain) = UnionDomain((d,))
+@deprecate UnionDomain(domain::Domain) UnionDomain((domain,))
 UnionDomain(domains) = _UnionDomain(promote_domains(domains))
-_UnionDomain(domains) = _UnionDomain(eltype(first(domains)), domains)
+_UnionDomain(domains) = _UnionDomain(domaineltype(first(domains)), domains)
 
 UnionDomain{T}(domains...) where {T} = UnionDomain{T}(domains)
-UnionDomain{T}(d::Domain) where {T} = UnionDomain{T}((d,))
+@deprecate UnionDomain{T}(domain::Domain) where {T} UnionDomain{T}((domain,))
 UnionDomain{T}(domains) where {T} = _UnionDomain(T, convert_eltype.(T, domains))
 _UnionDomain(::Type{T}, domains) where {T} = UnionDomain{T,typeof(domains)}(domains)
 
@@ -48,21 +50,19 @@ _UnionDomain(::Type{T}, domains) where {T} = UnionDomain{T,typeof(domains)}(doma
 composition(d::UnionDomain) = Combination()
 combine(d::UnionDomain, results) = reduce(|, results)
 
-# Make d1 ∪ d2 invoke `uniondomain` if one of the first two arguments is a Domain
-union(d1::Domain, d2::Domain, domains...) = uniondomain(d1, d2, domains...)
-union(d1::Domain, d2, domains...) = uniondomain(d1, d2, domains...)
-union(d1, d2::Domain, domains...) = uniondomain(d1, d2, domains...)
+# Make d1 ∪ d2 invoke `uniondomain` if one of the first two arguments is a domain
+union(d1::AnyDomain, d2::AnyDomain, domains...) =
+	uniondomain(domain(d1), domain(d2), domains...)
+union(d1::AnyDomain, d2, domains...) = uniondomain(domain(d1), d2, domains...)
+union(d1, d2::AnyDomain, domains...) = uniondomain(d1, domain(d2), domains...)
 
 uniondomain() = emptyspace(Any)
 uniondomain(d1) = d1
 uniondomain(d1, d2) = uniondomain1(promote_domains(d1, d2)...)
 
 # simplification: in case the domains are equal, we don't create a union
-uniondomain1(d1, d2) = d1 == d2 ? d1 : uniondomain2(d1, d2)
-# if d1 is a Domain, we go one step further and invoke `issubset`
-uniondomain1(d1::Domain, d2) = issubset(d2, d1) ? d1 : uniondomain2(d1, d2)
-uniondomain2(d1, d2) = UnionDomain(d1, d2)
-uniondomain2(d1, d2::Domain) = issubset(d1, d2) ? d2 : UnionDomain(d1, d2)
+uniondomain1(d1, d2) = issubset_domain(d2, d1) ? d1 : uniondomain2(d1, d2)
+uniondomain2(d1, d2) = issubset_domain(d1, d2) ? d2 : UnionDomain(d1, d2)
 
 uniondomain(d1, d2, d3) = _ud3(promote_domains(d1, d2, d3)...)
 _ud3(d1, d2, d3) =
@@ -105,11 +105,11 @@ convert(::Type{Domain{T}}, v::AbstractSet{<:Domain}) where {T} = UnionDomain{T}(
 convert(::Type{Domain{T}}, s::AbstractSet) where {T} = UnionDomain{T}(map(Point,collect(s)))
 
 similardomain(d::UnionDomain, ::Type{T}) where {T} =
-    UnionDomain(convert.(Domain{T}, components(d)))
+    UnionDomain(convert_eltype.(T, components(d)))
 
 
 
-point_in_domain(d::UnionDomain) = convert(eltype(d), point_in_domain(component(d,1)))
+choice(d::UnionDomain) = convert(eltype(d), choice(component(d,1)))
 
 isempty(d::UnionDomain) = all(isempty, d.domains)
 
@@ -173,16 +173,15 @@ end
 """
 The `IntersectDomain` constructor can be invoked in one of three ways:
 - with a list of arguments: IntersectDomain(d1, d2, ...)
-- with a single domain: IntersectDomain(d::Domain)
 - or with any iterable list of domains: IntersectDomain(domains)
 """
 IntersectDomain(domains...) = IntersectDomain(domains)
-IntersectDomain(d::Domain) = IntersectDomain((d,))
+@deprecate IntersectDomain(domain::Domain) IntersectDomain((domain,))
 IntersectDomain(domains) = _IntersectDomain(promote_domains(domains))
-_IntersectDomain(domains) = IntersectDomain{eltype(first(domains))}(domains)
+_IntersectDomain(domains) = IntersectDomain{domaineltype(first(domains))}(domains)
 
 IntersectDomain{T}(domains...) where {T} = IntersectDomain{T}(domains)
-IntersectDomain{T}(d::Domain) where {T} = IntersectDomain{T}((d,))
+@deprecate IntersectDomain{T}(domain::Domain) where T IntersectDomain((domain,))
 IntersectDomain{T}(domains) where {T} = _IntersectDomain(T, convert_eltype.(T, domains))
 _IntersectDomain(::Type{T}, domains) where {T} = IntersectDomain{T,typeof(domains)}(domains)
 
@@ -192,18 +191,21 @@ combine(d::IntersectDomain, results) = reduce(&, results)
 
 
 # Make d1 ∩ d2 invoke `intersectdomain` if one of the first two arguments is a Domain
-intersect(d1::Domain, d2::Domain, domains...) = intersectdomain(d1, d2, domains...)
-intersect(d1::Domain, d2, domains...) = intersectdomain(d1, d2, domains...)
-intersect(d1, d2::Domain, domains...) = intersectdomain(d1, d2, domains...)
+intersect(d1::AnyDomain, d2::AnyDomain, domains...) =
+	intersectdomain(domain(d1), domain(d2), domains...)
+# if a domain is combined with another argument that has a domain interpretation,
+# proceed with intersectdomain
+intersect(d1::AnyDomain, d2, domains...) = _intersect(domain(d1), d2, DomainStyle(d2), domains...)
+intersect(d1, d2::AnyDomain, domains...) = _intersect(d1, domain(d2), DomainStyle(d1), domains...)
+_intersect(d1, d2, ::IsDomain, domains...) = intersectdomain(d1, d2, domains...)
+_intersect(d1, d2, ::NotDomain, domains...) = intersectdomain(d1, d2, domains...)
 
 intersectdomain() = emptyspace(Any)
 intersectdomain(d1) = d1
 intersectdomain(d1, d2) = intersectdomain1(promote_domains(d1, d2)...)
 
-intersectdomain1(d1, d2) = d1 == d2 ? d1 : intersectdomain2(d1, d2)
-intersectdomain1(d1::Domain, d2) = issubset(d2, d1) ? d2 : intersectdomain2(d1, d2)
-intersectdomain2(d1, d2) = IntersectDomain(d1, d2)
-intersectdomain2(d1, d2::Domain) = issubset(d1, d2) ? d1 : IntersectDomain(d1, d2)
+intersectdomain1(d1, d2) = issubset_domain(d2, d1) ? d2 : intersectdomain2(d1, d2)
+intersectdomain2(d1, d2) = issubset_domain(d1, d2) ? d1 : IntersectDomain(d1, d2)
 
 intersectdomain(d1, d2, d3) = _id3(promote_domains(d1, d2, d3)...)
 _id3(d1, d2, d3) =
@@ -239,10 +241,12 @@ function intersectdomain(d1::UnionDomain, d2::UnionDomain)
     d1 == d2 && return d1
     uniondomain(intersectdomain.(Ref(d1), components(d2))...)
 end
-intersectdomain(d1::UnionDomain, d2::Domain) = uniondomain(intersectdomain.(d1.domains, Ref(d2))...)
-intersectdomain(d1::Domain, d2::UnionDomain) = uniondomain(intersectdomain.(Ref(d1), d2.domains)...)
+intersectdomain1(d1::UnionDomain, d2) = uniondomain(intersectdomain.(d1.domains, Ref(d2))...)
+intersectdomain2(d1, d2::UnionDomain) = uniondomain(intersectdomain.(Ref(d1), d2.domains)...)
 
-(&)(d1::Domain, d2::Domain) = intersectdomain(d1,d2)
+(&)(d1::AnyDomain, d2::AnyDomain) = intersectdomain(domain(d1),domain(d2))
+(&)(d1::AnyDomain, d2) = intersectdomain(domain(d1), d2)
+(&)(d1, d2::AnyDomain) = intersectdomain(d1, domain(d2))
 
 function intersectdomain(d1::ProductDomain, d2::ProductDomain)
 	if compatibleproductdims(d1, d2)
@@ -253,7 +257,7 @@ function intersectdomain(d1::ProductDomain, d2::ProductDomain)
 end
 
 similardomain(d::IntersectDomain, ::Type{T}) where {T} =
-    IntersectDomain(convert.(Domain{T}, components(d)))
+    IntersectDomain(convert_eltype.(T, components(d)))
 
 isequaldomain(a::IntersectDomain, b::IntersectDomain) = Set(components(a)) == Set(components(b))
 hash(d::IntersectDomain, h::UInt) = hashrec("IntersectDomain", Set(components(d)), h)
@@ -281,7 +285,7 @@ struct SetdiffDomain{T,DD} <: CompositeDomain{T}
 end
 
 SetdiffDomain(d1, d2) = _SetdiffDomain(promote_domains((d1, d2))...)
-_SetdiffDomain(d1, d2) = SetdiffDomain{eltype(d1)}((d1,d2))
+_SetdiffDomain(d1, d2) = SetdiffDomain{domaineltype(d1)}((d1,d2))
 SetdiffDomain{T}(domains) where {T} = SetdiffDomain{T,typeof(domains)}(domains)
 
 # The difference between two domains corresponds to a logical AND NOT of their characteristic functions
@@ -294,18 +298,18 @@ _approx_indomain(x, d::SetdiffDomain, comp::Combination, domains, tolerance) =
     approx_in(x, domains[1], tolerance) & !in(x, domains[2])
 
 similardomain(d::SetdiffDomain, ::Type{T}) where {T} =
-    SetdiffDomain(convert(Domain{T}, d.domains[1]), convert(Domain{T}, d.domains[2]))
+    SetdiffDomain(convert_eltype(T, d.domains[1]), convert_eltype(T, d.domains[2]))
 
 # use \ as a synomym for setdiff, in the context of domains (though, generically,
 # \ means left division in Julia)
-\(d1::Domain, d2::Domain) = setdiffdomain(d1, d2)
-\(d1::Domain, d2) = setdiffdomain(d1, d2)
-\(d1, d2::Domain) = setdiffdomain(d1, d2)
+\(d1::AnyDomain, d2::AnyDomain) = setdiffdomain(domain(d1), domain(d2))
+\(d1::AnyDomain, d2) = setdiffdomain(domain(d1), d2)
+\(d1, d2::AnyDomain) = setdiffdomain(d1, domain(d2))
 
 # Make setdiff invoke `setdiffdomain` if one of the arguments is a Domain
-setdiff(d1::Domain, d2::Domain) = setdiffdomain(d1, d2)
-setdiff(d1::Domain, d2) = setdiffdomain(d1, d2)
-setdiff(d1, d2::Domain) = setdiffdomain(d1, d2)
+setdiff(d1::AnyDomain, d2::AnyDomain) = setdiffdomain(domain(d1), domain(d2))
+setdiff(d1::AnyDomain, d2) = setdiffdomain(domain(d1), d2)
+setdiff(d1, d2::AnyDomain) = setdiffdomain(d1, domain(d2))
 
 setdiffdomain(d1, d2) = setdiffdomain1(promote_domains(d1, d2)...)
 setdiffdomain1(d1, d2) = setdiffdomain2(d1, d2)
@@ -313,17 +317,7 @@ setdiffdomain1(d1, d2) = setdiffdomain2(d1, d2)
 function setdiffdomain2(d1, d2)
 	if isempty(d2)
 		d1
-	elseif d1 == d2
-		emptyspace(d1)
-	else
-		SetdiffDomain(d1, d2)
-	end
-end
-# similar to uniondomain, if d2 is a Domain we use `issubset` as an optimization
-function setdiffdomain2(d1, d2::Domain)
-	if isempty(d2)
-		d1
-	elseif issubset(d1,d2)
+	elseif issubset_domain(d1,d2)
 		emptyspace(d1)
 	else
 		SetdiffDomain(d1, d2)
diff --git a/src/maps/lazy.jl b/src/maps/lazy.jl
index cbd14cb..f0cd612 100644
--- a/src/maps/lazy.jl
+++ b/src/maps/lazy.jl
@@ -41,6 +41,7 @@ similarmap(m::WrappedMap, ::Type{T}) where {T} = WrappedMap{T}(m)
 convert(::Type{Map}, m::Map) = m
 convert(::Type{Map}, m) = WrappedMap(m)
 convert(::Type{Map{T}}, m) where {T} = WrappedMap{T}(m)
+convert(::Type{Map{T}}, m::WrappedMap) where {T} = WrappedMap{T}(m.map)
 
 ==(m1::WrappedMap, m2::Function) = m1.map == m2
 ==(m1::Function, m2::WrappedMap) = m1 == m2.map
diff --git a/src/maps/map.jl b/src/maps/map.jl
index 2c149d0..9377d8e 100644
--- a/src/maps/map.jl
+++ b/src/maps/map.jl
@@ -33,13 +33,13 @@ codomaintype(::Type{M}, ::Type{T}) where {M,T} = Any
 codomaintype(M::Type{<:AbstractMap}, ::Type{T}) where {T} = Base.promote_op(applymap, M, T)
 codomaintype(M::Type{<:TypedMap{T,U}}, ::Type{T}) where {T,U} = U
 
+prectype(::Type{<:Map{T}}) where T = prectype(T)
+numtype(::Type{<:Map{T}}) where T = numtype(T)
+
 isreal(m::AbstractMap) = isreal(domaintype(m)) && isreal(codomaintype(m))
 isreal(::UniformScaling{T}) where {T} = isreal(T)
 isreal(::Type{UniformScaling{T}}) where {T} = isreal(T)
 
-numtype(::Type{<:Map{T}}) where {T} = numtype(T)
-prectype(::Type{<:Map{T}}) where {T} = prectype(T)
-
 convert(::Type{AbstractMap}, m::AbstractMap) = m
 convert(::Type{Map{T}}, m::Map{T}) where {T} = m
 convert(::Type{Map{T}}, m::Map{S}) where {S,T} = similarmap(m, T)
diff --git a/src/util/common.jl b/src/util/common.jl
index ea6377a..dd219d7 100644
--- a/src/util/common.jl
+++ b/src/util/common.jl
@@ -10,7 +10,8 @@ const StaticTypes = Union{Number,<:StaticVector{N} where N,<:NTuple{N,Any} where
 euclideandimension(::Type{T}) where {T <: Number} = 1
 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
+euclideandimension(::Type{T}) where {T} =
+	error("Don't know the euclidean dimension of $(T).")
 
 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}
@@ -22,7 +23,7 @@ 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))
+	p = similar(choice(d))
 	fill!(p, 0)
 	convert(T, p)
 end
@@ -62,8 +63,11 @@ convert_eltype(::Type{T}, d::Number) where {T} = convert(T, d)
 
 "The floating point precision type associated with the argument."
 prectype(x) = prectype(typeof(x))
+prectype(::Type{T}) where T = Any		# fallback
+prectype(::Type{T}) where {T<:AbstractFloat} = T
+prectype(::Type{T}) where {T<:Number} = prectype(float(T))
+prectype(::Type{<:AbstractArray{T}}) where T = prectype(T)
 prectype(::Type{<:Complex{T}}) where {T} = prectype(T)
-prectype(::Type{<:AbstractArray{T}}) where {T} = prectype(T)
 prectype(::Type{NTuple{N,T}}) where {N,T} = prectype(T)
 prectype(::Type{Tuple{A}}) where {A} = prectype(A)
 prectype(::Type{Tuple{A,B}}) where {A,B} = prectype(A,B)
@@ -71,18 +75,12 @@ prectype(::Type{Tuple{A,B,C}}) where {A,B,C} = prectype(A,B,C)
 @generated function prectype(T::Type{<:Tuple})
     quote $(promote_type(map(prectype, T.parameters[1].parameters)...)) end
 end
-prectype(::Type{T}) where {T<:AbstractFloat} = T
-prectype(::Type{T}) where {T<:Number} = prectype(float(T))
 
 prectype(a...) = promote_type(map(prectype, a)...)
 
 "Convert `x` such that its `prectype` equals `U`."
-convert_prectype(x, ::Type{U}) where {U} = convert(to_prectype(typeof(x),U), x)
-
-# function convert_prectype(x, ::Type{U}) where {U}
-#     @warn "The order of arguments of convert_prectype has changed."
-#     prectype(U, x)
-# end
+convert_prectype(x, ::Type{U}) where {U} =
+	convert_eltype(to_prectype(eltype(x), U), x)
 
 "Return the type to which `U` can be converted, such that the `prectype` becomes `T`."
 to_prectype(::Type{T}, ::Type{U}) where {T,U} = error("Don't know how to convert the numtype of $(T) to $(U).")
@@ -104,10 +102,11 @@ _promote_prectype(U, a, b, c...) =
 # Numeric type
 #################
 
-"The numeric element type of x in a Euclidean space."
+"The numeric element type of `x` in a Euclidean space."
 numtype(x) = numtype(typeof(x))
+numtype(::Type{T}) where T = Any
 numtype(::Type{T}) where {T<:Number} = T
-numtype(::Type{T}) where {T} = eltype(T)
+numtype(::Type{<:AbstractArray{T}}) where T = T
 numtype(::Type{NTuple{N,T}}) where {N,T} = T
 numtype(::Type{Tuple{A,B}}) where {A,B} = promote_type(numtype(A), numtype(B))
 numtype(::Type{Tuple{A,B,C}}) where {A,B,C} = promote_type(numtype(A), numtype(B), numtype(C))
@@ -119,7 +118,8 @@ end
 numtype(a...) = promote_type(map(numtype, a)...)
 
 "Convert `x` such that its `numtype` equals `U`."
-convert_numtype(x, ::Type{U}) where {U} = convert(to_numtype(typeof(x),U), x)
+convert_numtype(x, ::Type{U}) where {U} =
+	convert_eltype(to_numtype(eltype(x), U), x)
 
 to_numtype(::Type{T}, ::Type{U}) where {T,U} = error("Don't know how to convert the numtype of $(T) to $(U).")
 to_numtype(::Type{T}, ::Type{U}) where {T <: Number,U <: Number} = U
diff --git a/test/runtests.jl b/test/runtests.jl
index 4567da4..57f5144 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -3,6 +3,9 @@ using Test, LinearAlgebra, StaticArrays, Random, StableRNGs
 using DomainSets
 using CompositeTypes.Indexing
 
+const io = IOBuffer()
+const textmime = MIME"text/plain"()
+
 include("test_common.jl")
 include("test_maps.jl")
 include("test_generic_domain.jl")
diff --git a/test/test_applications.jl b/test/test_applications.jl
index ecc564b..1337e9a 100644
--- a/test/test_applications.jl
+++ b/test/test_applications.jl
@@ -103,7 +103,7 @@ function test_rand(T)
         rs = rand(rng, b, n)
         n_1 = sum(r in region_1 for r in rs)
         n_2 = sum(r in region_2 for r in rs)
-        @test isapprox(n_1, n_2, rtol=0.1) 
+        @test isapprox(n_1, n_2, rtol=0.1)
     end
 end
 
diff --git a/test/test_common.jl b/test/test_common.jl
index 839b834..d8175e6 100644
--- a/test/test_common.jl
+++ b/test/test_common.jl
@@ -1,16 +1,23 @@
-
-using DomainSets: convert_numtype, convert_prectype,
-    promote_numtype, promote_prectype
+using DomainSets:
+    convert_eltype,
+    convert_numtype,
+    convert_prectype,
+    promote_numtype,
+    promote_prectype,
+    euclideandimension,
+    domaineltype,
+    domain_numtype,
+    domain_prectype
 
 function test_dimension()
-    @test DomainSets.euclideandimension(Int) == 1
-    @test DomainSets.euclideandimension(Float64) == 1
-    @test DomainSets.euclideandimension(ComplexF64) == 1
-    @test DomainSets.euclideandimension(SVector{2,Float64}) == 2
-    @test DomainSets.euclideandimension(MVector{2,Float64}) == 2
-    @test DomainSets.euclideandimension(Tuple{Int,Int}) == 2
-    @test DomainSets.euclideandimension(Tuple{Int,Float64}) == 2
-    @test_throws MethodError DomainSets.euclideandimension(Vector{Float64})
+    @test euclideandimension(Int) == 1
+    @test euclideandimension(Float64) == 1
+    @test euclideandimension(ComplexF64) == 1
+    @test euclideandimension(SVector{2,Float64}) == 2
+    @test euclideandimension(MVector{2,Float64}) == 2
+    @test euclideandimension(Tuple{Int,Int}) == 2
+    @test euclideandimension(Tuple{Int,Float64}) == 2
+    @test_throws ErrorException euclideandimension(Vector{Float64})
 end
 
 function test_components()
@@ -29,13 +36,13 @@ function test_prectype()
     @test prectype(1.0+2.0im) == Float64
     @test prectype([1.0+2.0im, 3.0]) == Float64
     @test prectype(NTuple{2,Int}) == Float64
-    @test prectype((1.0,)) == Float64
-    @test prectype((1.0, 2.0)) == Float64
-    @test prectype((1.0, 2.0, 3.0)) == Float64
-    @test prectype((1.0, big(2.0), 3.0+im)) == BigFloat
+    @test prectype(typeof((1.0,))) == Float64
+    @test prectype(typeof((1.0, 2.0))) == Float64
+    @test prectype(typeof((1.0, 2.0, 3.0))) == Float64
+    @test prectype(typeof((1.0, big(2.0), 3.0+im))) == BigFloat
     @test prectype(NTuple{4,Int}) == Float64
     @test @inferred(prectype(1, 2.0)) == Float64
-    @test @inferred(prectype((1, 2.0, 3, 40+im))) == Float64
+    @test @inferred(prectype(typeof((1, 2.0, 3, 40+im)))) == Float64
 
     @test convert_prectype(2, Float64) == 2
     @test convert_prectype(2, Float64) isa Float64
@@ -63,14 +70,14 @@ function test_numtype()
     @test numtype(1.0+2.0im) == Complex{Float64}
     @test numtype([1.0+2.0im, 3.0]) == Complex{Float64}
     @test numtype(NTuple{2,Complex{Int}}) == Complex{Int}
-    @test numtype((1.0,)) == Float64
-    @test numtype((1.0, 2.0)) == Float64
+    @test numtype(typeof((1.0,))) == Float64
+    @test numtype(typeof((1.0, 2.0))) == Float64
     @test numtype((1.0, 2.0, 3.0)) == Float64
     @test numtype((1.0, 2.0, 3.0, 4.0)) == Float64
     @test numtype(1.0, big(2.0), 3.0+im) == Complex{BigFloat}
-    @test numtype((1.0, big(2.0), 3.0+im)) == Complex{BigFloat}
+    @test numtype(typeof((1.0, big(2.0), 3.0+im))) == Complex{BigFloat}
     @test @inferred(numtype(1, 2.0)) == Float64
-    @test @inferred(numtype((1, 2.0, 3, 40+im))) == Complex{Float64}
+    @test @inferred(numtype(typeof((1, 2.0, 3, 40+im)))) == Complex{Float64}
 
     @test convert_numtype(2, Float64) == 2
     @test convert_numtype(2, Float64) isa Float64
@@ -83,7 +90,6 @@ function test_numtype()
     @test promote_numtype(2, 3.0+im, big(4)) isa Tuple{Complex{BigFloat},Complex{BigFloat},Complex{BigFloat}}
 end
 
-using DomainSets: convert_eltype
 function test_eltype()
     @test convert_eltype(Float64, Point(0)) isa Point{Float64}
     @test convert_eltype(Float64, Point(0)) == Point(0)
diff --git a/test/test_domain_interval.jl b/test/test_domain_interval.jl
index 10a962b..921a019 100644
--- a/test/test_domain_interval.jl
+++ b/test/test_domain_interval.jl
@@ -232,8 +232,8 @@ function test_intervals()
         @test isopenset(d)
         @test DomainSets.isfullspace(d)
         @test 1 ∈ d
-        @test point_in_domain(d) ∈ d
-        @test point_in_domain(d) isa T
+        @test choice(d) ∈ d
+        @test choice(d) isa T
         @test isempty(boundary(d))
         @test_throws ArgumentError minimum(d)
         @test_throws ArgumentError maximum(d)
diff --git a/test/test_domain_product.jl b/test/test_domain_product.jl
index 28e5a47..f7369db 100644
--- a/test/test_domain_product.jl
+++ b/test/test_domain_product.jl
@@ -42,9 +42,9 @@ function test_product_domains()
         @test_logs (:warn, "`in`: incompatible combination of vector with length 3 and domain '($(-1.0..1.0)) × ($(-1.0..1.0))' with dimension 2. Returning false.") [0.0,0.0,0.0] ∉ d1
         @test_logs (:warn, "`in`: incompatible combination of vector with length 1 and domain '($(-1.0..1.0)) × ($(-1.0..1.0))' with dimension 2. Returning false.") [0.0] ∉ d1
 
-        d3 = VcatDomain(-1.0 .. 1.0, -1.5 .. 2.5)
-        @test SA[0.5,0.5] ∈ d3
-        @test SA[-1.1,0.3] ∉ d3
+        d2 = VcatDomain(-1.0 .. 1.0, -1.5 .. 2.5)
+        @test SA[0.5,0.5] ∈ d2
+        @test SA[-1.1,0.3] ∉ d2
 
         d3 = VcatDomain(1.05 * UnitDisk(), -1.0 .. 1.0)
         @inferred(cross(1.05 * UnitDisk(), -1.0 .. 1.0)) === d3
@@ -52,7 +52,7 @@ function test_product_domains()
         @test eltype(d3) == SVector{3,Float64}
         @test SA[0.5,0.5,0.8] ∈ d3
         @test SA[-1.1,0.3,0.1] ∉ d3
-        @test point_in_domain(d3) ∈ d3
+        @test choice(d3) ∈ d3
     end
     @testset "mixed intervals" begin
         d = (0..1) × (0.0..1)
@@ -63,7 +63,7 @@ function test_product_domains()
         # Make sure promotion of domains happened
         @test eltype(component(d,1)) == Float64
         @test eltype(component(d,2)) == Float64
-        @test point_in_domain(d) ∈ d
+        @test choice(d) ∈ d
     end
     @testset "vector domains" begin
         d1 = VectorProductDomain([0..1.0, 0..2.0])
@@ -72,7 +72,7 @@ function test_product_domains()
         @test dimension(d1) == 2
         @test [0.1,0.2] ∈ d1
         @test SA[0.1,0.2] ∈ d1
-        @test point_in_domain(d1) ∈ d1
+        @test choice(d1) ∈ d1
         @test convert(Domain{Vector{BigFloat}}, d1) == d1
         d1big = convert(Domain{Vector{BigFloat}}, d1)
         @test eltype(d1big) == Vector{BigFloat}
@@ -81,7 +81,7 @@ function test_product_domains()
         d2 = VectorProductDomain([0..1, 0..3])
         @test dimension(d2) == 2
         @test [0.1,0.2] ∈ d2
-        @test point_in_domain(d2) ∈ d2
+        @test choice(d2) ∈ d2
 
         # other constructor calls
         @test VectorProductDomain(0..1, 1..2) isa VectorProductDomain{Vector{Int}}
@@ -110,7 +110,7 @@ function test_product_domains()
         @test rand(10) ∈ d4
         @test 2 .+ rand(10) ∉ d4
 
-        @test VectorProductDomain{SVector{2,Float64}}(SVector(0..1, 0..2)).domains[1] isa Domain{Float64}
+        @test eltype(VectorProductDomain{SVector{2,Float64}}(SVector(0..1, 0..2)).domains[1]) == Float64
     end
     @testset "Tuple product domain" begin
         # Use the constructor ProductDomain{T} directly
@@ -240,7 +240,7 @@ function test_product_domains()
         @test d4 isa Rectangle
         @test SA[0.5,0.5,0.8] ∈ d4
         @test SA[-1.1,0.3,0.1] ∉ d4
-        @test point_in_domain(d4) ∈ d4
+        @test choice(d4) ∈ d4
         @test center(d4) == [0.0,0.0,0.0]
 
         @test d1[Component(1)] == -1..1
@@ -251,13 +251,13 @@ function test_product_domains()
         @test d5 isa Rectangle
         @test SA[0.,0.5,0.5] ∈ d5
         @test SA[0.,-1.1,0.3] ∉ d5
-        @test point_in_domain(d5) ∈ d5
+        @test choice(d5) ∈ d5
 
         d6 = d1 × d1
         @test d6 isa Rectangle
         @test SA[0.,0.,0.5,0.5] ∈ d6
         @test SA[0.,0.,-1.1,0.3] ∉ d6
-        @test point_in_domain(d6) ∈ d6
+        @test choice(d6) ∈ d6
 
         @test Rectangle( SA[1,2], SA[2.0,3.0]) isa Rectangle{SVector{2,Float64}}
         @test Rectangle([0..1, 2..3]) isa Rectangle{Vector{Int}}
diff --git a/test/test_domain_simplex.jl b/test/test_domain_simplex.jl
index f8b520f..4af60cc 100644
--- a/test/test_domain_simplex.jl
+++ b/test/test_domain_simplex.jl
@@ -70,7 +70,7 @@ function test_simplex()
     @test !isopenset(d)
     @test isopenset(interior(d))
     @test closure(d) == d
-    @test point_in_domain(d) ∈ d
+    @test choice(d) ∈ d
 
     # open/closed
     d2 = EuclideanUnitSimplex{2,Float64,:open}()
@@ -83,7 +83,7 @@ function test_simplex()
     @test !approx_in(SA[-0.01,0.0], d2, 0.001)
 
     d3 = EuclideanUnitSimplex{3,BigFloat}()
-    @test point_in_domain(d3) ∈ d3
+    @test choice(d3) ∈ d3
     x0 = big(0.0)
     x1 = big(1.0)
     x2 = big(0.3)
@@ -119,7 +119,7 @@ function test_simplex()
     @test SA[0.2,-0.2] ∉ D
     @test convert(Domain{Vector{BigFloat}}, D) == VectorUnitSimplex{BigFloat}(2)
     @test corners(D) == [ [0.0,0.0], [1.0,0.0], [0.0,1.0]]
-    @test boundingbox(D) == UnitCube(4)
+    @test boundingbox(D) == UnitCube(2)
 
     d4 = EuclideanUnitSimplex{4,Float64}()
     @test corners(d4) isa SVector{5,SVector{4,Float64}}
diff --git a/test/test_generic_domain.jl b/test/test_generic_domain.jl
index b7af73c..380c43f 100644
--- a/test/test_generic_domain.jl
+++ b/test/test_generic_domain.jl
@@ -1,4 +1,4 @@
-import DomainSets: factors, nfactors, factor
+using DomainSets: factors, nfactors, factor
 
 
 widen_eltype(::Type{T}) where {T<:Number} = widen(T)
@@ -9,16 +9,21 @@ widen_eltype(::Type{Vector{T}}) where {T<:Number} = Vector{widen(T)}
 # We test the generic functionality of a domain.
 # These tests check whether the given domain correctly implements the
 # interface of a domain.
-function test_generic_domain(d::Domain)
-    @test isreal(d) == isreal(eltype(d))
-    @test isreal(d) == isreal(numtype(d))
+function test_generic_domain(d)
+    @test isreal(d) == isreal(domaineltype(d))
+    @test isreal(d) == isreal(domain_numtype(d))
 
-    @test convert(Domain{eltype(d)}, d) == d
-    @test convert(Domain{widen_eltype(eltype(d))}, d) == d
+    if d isa Domain
+        @test convert(Domain{eltype(d)}, d) == d
+        @test convert(Domain{widen_eltype(eltype(d))}, d) == d
+    elseif DomainStyle(d) isa IsDomain
+        @test convert_eltype(eltype(d), d) == d
+        @test convert_eltype(widen_eltype(eltype(d)), d) == d
+    end
     @test prectype(convert_prectype(d, BigFloat)) == BigFloat
 
     if !isempty(d)
-        x = point_in_domain(d)
+        x = choice(d)
         @test x ∈ d
         @test approx_in(x, d, 0.01)
         @test_throws ErrorException approx_in(x, d, -1)
@@ -26,17 +31,17 @@ function test_generic_domain(d::Domain)
         @test approx_in.([x,x], d, 0.01) == approx_in.([x,x], Ref(d), 0.01)
     else
         try
-            x = point_in_domain(d)
+            x = choice(d)
             @test false
         catch
         end
     end
-    @test canonicaldomain(DomainSets.Equal(), d) == d
+    @test isequaldomain(canonicaldomain(DomainSets.Equal(), d), d)
     if hascanonicaldomain(d)
         cd = canonicaldomain(d)
         @test mapfrom_canonical(d) == mapto(cd, d)
         @test mapto_canonical(d) == mapto(d, cd)
-        x1 = point_in_domain(cd)
+        x1 = choice(cd)
         @test mapfrom_canonical(d, x1) ∈ d
         @test mapto_canonical(d, x) ∈ cd
     else
@@ -46,7 +51,7 @@ function test_generic_domain(d::Domain)
     if hasparameterization(d)
         par = parameterdomain(d)
         @test mapfrom_parameterdomain(d) == mapto(par, d)
-        xp = point_in_domain(par)
+        xp = choice(par)
         @test approx_in(mapfrom_parameterdomain(d, xp), d)
     end
     if iscomposite(d)
@@ -109,24 +114,19 @@ end
         # functionality using broadcast
         @test 2 * (1..2) == 2 .* (1..2)
         @test (1..2) * 2 == (1..2) .* 2
-        @test (1..2) / 2 ≈ (0.5..1)
-        @test 2 \ (1..2) ≈ (0.5..1)
+        @test (1..2) / 2 ≈ 0.5..1
+        @test 2 \ (1..2) ≈ 0.5..1
         @test all(rand(4) .∈ (-1..1))
-        @test all(approx_in.(rand(4), -1..1))
-        @test all(approx_in.([1.005,1.0005], -1..1, [1e-2,1e-3]))
+        @test all(approx_in.(rand(4), (-1..1)))
+        @test all(approx_in.([1.005,1.0005], (-1..1), [1e-2,1e-3]))
         @test all(rand(4) .∉ (2..3))
         @test_throws MethodError (0..1) + 0.4
 
         # promotion
         @test DomainSets.promote_domains() == ()
-        s1 = Set([0..1,2..3,3..4.0])
-        @test s1 isa Set{<:Domain{Float64}}
-        @test DomainSets.promote_domains(s1) == s1
-        s2 = Set([0..1,2..3,3..4.0, Point(2)])
-        @test s2 isa Set{Domain}
-        @test DomainSets.promote_domains(s2) isa Set{<:Domain{Float64}}
-        @test DomainSets.promote(0..1.0, [1,2,3]) isa Tuple{Interval,Vector{Float64}}
-        @test DomainSets.promote([1,2,3], 0..1.0) isa Tuple{Vector{Float64},Interval}
+        @test DomainSets.promote_domains(0..1, 2..4.0) isa Tuple{ClosedInterval{Float64},ClosedInterval{Float64}}
+        @test DomainSets.promote_domains(0..1.0, [1,2,3]) isa Tuple{Interval,Vector{Float64}}
+        @test DomainSets.promote_domains([1,2,3], 0..1.0) isa Tuple{Vector{Float64},Interval}
 
         # compatible point-domain pairs
         @test DomainSets.iscompatiblepair(0.5, 0..1)
@@ -143,4 +143,7 @@ end
         @test_throws ErrorException DomainSets.convert_eltype(Float64, (1,2)) isa NTuple{2,Float64}
         @test DomainSets.convert_eltype(Int, (1,2)) == (1,2)
     end
+
+    @test choice(Set([1,2,3])) ∈ Set([1,2,3])
+    @test choice([1,2,3]) == 1
 end
diff --git a/test/test_maps.jl b/test/test_maps.jl
index 3fea8bc..b4b21a6 100644
--- a/test/test_maps.jl
+++ b/test/test_maps.jl
@@ -29,9 +29,9 @@ maps_to_test(T) = [
     AffineMap(rand(T, 3, 2), rand(T, 3)), # vector map, rectangular
     AffineMap(LinearAlgebra.I, one(T)/2),  # use UniformScaling object as A
     AffineMap(LinearAlgebra.I, randvec(T,2)),  # use UniformScaling object as A
-    DomainSets.GenericAffineMap(randvec(T, 2, 2), randvec(T, 2)),
-    DomainSets.GenericAffineMap(T(1.2), randvec(T, 2)),
-    DomainSets.GenericAffineMap(randvec(T, 3, 2), randvec(T, 3)),
+    GenericAffineMap(randvec(T, 2, 2), randvec(T, 2)),
+    GenericAffineMap(T(1.2), randvec(T, 2)),
+    GenericAffineMap(randvec(T, 3, 2), randvec(T, 3)),
     Translation(randvec(T, 3)),
     LinearMap(randvec(T, 2, 2)),
     LinearMap(randvec(T, 2)),
diff --git a/test/test_readme.jl b/test/test_readme.jl
index 54edadb..6cc21fc 100644
--- a/test/test_readme.jl
+++ b/test/test_readme.jl
@@ -4,7 +4,7 @@
     using DomainSets, StaticArrays
     @test repr(UnitInterval()) == "$(0.0..1.0) (Unit)"
     @test repr(ChebyshevInterval()) == "$(-1.0..1.0) (Chebyshev)"
-    @test repr(HalfLine()) == "$(0.0..Inf) (closed-open) (HalfLine)"
+    @test repr(HalfLine())[1:10] == "$(0.0..Inf)"
 
     using DomainSets: ×
     @test repr((-1..1) × (0..3) × (4.0..5.0)) == "($(-1.0..1.0)) × ($(0.0..3.0)) × ($(4.0..5.0))"
diff --git a/test/test_setoperations.jl b/test/test_setoperations.jl
index 3efc36b..b8469f3 100644
--- a/test/test_setoperations.jl
+++ b/test/test_setoperations.jl
@@ -29,7 +29,7 @@
         ũ2 = UnionDomain([d1,d2])
         @test ũ2 == ũ2
         @test u1 == ũ2
-        @test UnionDomain{SVector{2,Float64}}(d1) isa UnionDomain
+        @test UnionDomain{SVector{2,Float64}}((d1,)) isa UnionDomain
 
         # Don't create a union with two identical elements
         @test UnitDisk() ∪ UnitDisk() isa UnitDisk
@@ -74,14 +74,12 @@
         @test UnionDomain(d1,d2) == UnionDomain(d2,d1)
 
         @test UnionDomain((d1,d2)) == UnionDomain(d1,d2)
-        @test UnionDomain(d1) isa UnionDomain
         @test UnionDomain(UnionDomain(d1,d2),d3) == UnionDomain(d3,UnionDomain(d1,d2))
 
-        @test convert(Domain, [0..1,2..3,3..4]) isa UnionDomain{Int}
-        @test convert(Domain{Float64}, [0..1,2..3,3..4]) isa UnionDomain{Float64}
-        @test convert(Domain, Set([0..1,2..3,3..4])) isa UnionDomain{Int}
-        @test convert(Domain{Float64}, Set([0..1,2..3,3..4])) isa UnionDomain{Float64}
-        @test convert(Domain, Set([1,2,3])) isa UnionDomain{Int}
+        @test convert(Domain, [Point(0),Point(2),Point(3)]) isa UnionDomain{Int}
+        @test convert(Domain{Float64}, [Point(1),Point(2),Point(3)]) isa UnionDomain{Float64}
+        @test convert(Domain, Set([Point(0),Point(2),Point(3)])) isa UnionDomain{Int}
+        @test convert(Domain{Float64}, Set([Point(0),Point(2),Point(3)])) isa UnionDomain{Float64}
         @test 2 ∈ convert(Domain, Set([1,2,3]))
         @test 2 ∈ convert(Domain{Float64}, Set([1,2,3]))
         @test 4 ∉ convert(Domain, Set([1,2,3]))
@@ -97,14 +95,12 @@
         @test String(take!(io)) == "UnitDisk() ∪ (($(-0.9..0.9)) × ($(-0.9..0.9)))"
 
         # repeated union
-        @test ncomponents(uniondomain(UnitBall{Float64}(), UnitInterval(), UnitInterval())) == 2
-        @test ncomponents(uniondomain(UnitInterval(), UnitBall{Float64}(), UnitInterval())) == 2
-        @test ncomponents(uniondomain(UnitInterval(), UnitInterval(), UnitBall{Float64}())) == 2
+        @test ncomponents(uniondomain(-2 .. -1, UnitInterval(), UnitInterval())) == 2
+        @test ncomponents(uniondomain(UnitInterval(), -2 .. -1, UnitInterval())) == 2
+        @test ncomponents(uniondomain(UnitInterval(), UnitInterval(), -2 .. -1)) == 2
     end
 
     @testset "intersect" begin
-        @test IntersectDomain(0..1) == IntersectDomain((0..1))
-        @test IntersectDomain{Float64}(0..1) == IntersectDomain(0.0..1.0)
         @test IntersectDomain{Float64}(0..1, 1..2) == IntersectDomain((0..1, 1..2))
         @test intersectdomain(0..1, 1..2) == Point(1)
 
@@ -154,9 +150,9 @@
         @test intersectdomain(2..4, uniondomain(0..1, 2..3)) == 2..3
 
         # repeated intersection
-        @test ncomponents(intersectdomain(UnitBall{Float64}(), UnitInterval(), UnitInterval())) == 2
-        @test ncomponents(intersectdomain(UnitInterval(), UnitBall{Float64}(), UnitInterval())) == 2
-        @test ncomponents(intersectdomain(UnitInterval(), UnitInterval(), UnitBall{Float64}())) == 2
+        @test ncomponents(intersectdomain(UnitBall{Float64}(), UnitInterval(), UnitInterval())) == 0
+        @test ncomponents(intersectdomain(UnitInterval(), UnitBall{Float64}(), UnitInterval())) == 0
+        @test ncomponents(intersectdomain(UnitInterval(), UnitInterval(), UnitBall{Float64}())) == 0
         # larger intersection expressions
         @test intersectdomain(0..1, 1..3, Point(0.4), 2..5, FullSpace(), Point(-0.2)) isa EmptySpace
         @test intersectdomain(0..1, 1..3, Point(1.0)) == Point(1.0)
@@ -198,7 +194,7 @@
         @test setdiff(0.5, 0..1) == setdiffdomain(0.5, 0..1)
 
         @test setdiff(0..1, EmptySpace()) == 0..1
-        @test setdiff(0..1, 0.0..1.0) == EmptySpace()
+        @test setdiffdomain(0..1, 0.0..1.0) == EmptySpace()
 
         @test (0..1)^2 \ UnitCircle() == UnitInterval()^2 \ UnitCircle()
     end
@@ -206,25 +202,25 @@
     @testset "arithmetic" begin
         d1 = (0..1)
         d2 = (2..3)
-        d = UnionDomain(d1) ∪ UnionDomain(d2)
+        d = UnionDomain((d1,)) ∪ UnionDomain((d2,))
 
-        @test d .+ 1 == UnionDomain(d1 .+ 1) ∪ (d2 .+ 1)
-        @test d .- 1 == UnionDomain(d1 .- 1) ∪ (d2 .- 1)
-        @test 2 * d  == UnionDomain(2 * d1)  ∪ (2 * d2)
-        @test d * 2 == UnionDomain(d1 * 2) ∪ (d2 * 2)
-        @test d / 2 == UnionDomain(d1 / 2) ∪ (d2 / 2)
-        @test 2 \ d == UnionDomain(2 \ d1) ∪ (2 \ d2)
+        @test d .+ 1 == UnionDomain((d1 .+ 1,)) ∪ (d2 .+ 1)
+        @test d .- 1 == UnionDomain((d1 .- 1,)) ∪ (d2 .- 1)
+        @test 2 * d  == UnionDomain((2*d1,))  ∪ 2 * d2
+        @test d * 2 == UnionDomain((d1 * 2,)) ∪ (d2 * 2)
+        @test d / 2 == UnionDomain((d1 / 2,)) ∪ (d2 / 2)
+        @test 2 \ d == UnionDomain((2 \ d1,)) ∪ (2 \ d2)
 
         @test infimum(d) == minimum(d) == 0
         @test supremum(d) == maximum(d) == 3
     end
 
     @testset "different types" begin
-        d̃1 = (0..1)
-        d1 = (0f0.. 1f0)
-        d2 = (2..3)
+        d1 = 0..1
+        d2 = 0f0..1f0
+        d3 = 2..3
 
-        @test UnionDomain(d1) ∪ d2 == UnionDomain(d̃1) ∪ d2
+        @test UnionDomain(d1,d2) ∪ d3 == d1 ∪ UnionDomain(d2, d3)
     end
 
     @testset "disk × interval" begin
diff --git a/test/test_specific_domains.jl b/test/test_specific_domains.jl
index 463aaa6..539ae0b 100644
--- a/test/test_specific_domains.jl
+++ b/test/test_specific_domains.jl
@@ -9,9 +9,6 @@ struct Basis3Vector <: StaticVector{3,Float64} end
 
 Base.getindex(::Basis3Vector, k::Int) = k == 1 ? 1.0 : 0.0
 
-const io = IOBuffer()
-const textmime = MIME"text/plain"()
-
 struct NamedBall <: DomainSets.DerivedDomain{SVector{2,Float64}}
     domain  ::  Domain{SVector{2,Float64}}
 
@@ -84,7 +81,7 @@ include("test_domain_simplex.jl")
         @test convert(Domain{BigFloat}, d1) === FullSpace{BigFloat}()
         @test DomainSets.euclideanspace(Val{2}()) == FullSpace{SVector{2,Float64}}()
         @test 0.5 ∈ d1
-        @test point_in_domain(d1) == 0
+        @test choice(d1) == 0
         @test d1 ∪ d1 == d1
         @test d1 ∩ d1 == d1
         @test isempty(d1) == false
@@ -259,7 +256,7 @@ include("test_domain_simplex.jl")
         @test Number(Point(1)) ≡ convert(Number, Point(1)) ≡ convert(Int, Point(1)) ≡ 1
         @test convert(Domain{Float64}, 1) isa Point{Float64}
 
-        @test point_in_domain(Point(1)) == 1
+        @test choice(Point(1)) == 1
 
         @test Point(1) + Point(2) == Point(3)
         @test Point(1) - Point(2) == Point(-1)
@@ -287,9 +284,9 @@ include("test_domain_simplex.jl")
         @test repr(Point(3.0)) == "Point(3.0)"
     end
 
-    @testset "intervals" begin
-        test_intervals()
-    end
+    # @testset "intervals" begin
+    #     test_intervals()
+    # end
 
     @testset "balls" begin
         test_balls()
@@ -422,7 +419,7 @@ include("test_domain_simplex.jl")
         D = rotate(UnitInterval()^2, π/2)
         @test SA[-0.9, 0.9] ∈ D
         @test SA[-1.1, -1.1] ∉ D
-        x = point_in_domain(D)
+        x = choice(D)
         @test forward_map(D)(x) ≈ forward_map(D, x)
         @test DomainSets.toexternalpoint(D, x) ≈ forward_map(D, x)
 
@@ -606,7 +603,7 @@ include("test_domain_simplex.jl")
         @test [0.4,-0.2] ∉ d3
 
         d4 = Domain( x+y+z > 0 for (x,y) in UnitDisk(), z in 0..1.0)
-        @test d4 isa DomainSets.BoundedIndicatorFunction{F,<:TupleProductDomain} where F
+        @test d4 isa DomainSets.BoundedIndicatorFunction{T,F,<:TupleProductDomain} where {T,F}
         @test eltype(d4) == Tuple{SVector{2, Float64}, Float64}
         @test DomainSets.indicatorfunction(d4) isa Function
         @test DomainSets.boundingdomain(d4) == TupleProductDomain(UnitDisk(), 0..1.0)