improve support for boundary, gradient and normal

[daanhb]

Some domains already define the boundary function, e.g., the boundary of the unit disk is the unit circle. It would be good to extend this to all domains that can be used in geometry. At the same time, we can define gradients and normal directions when they make sense.

Right now we have

julia> using DomainSets; boundary(UnitCube())
D₄ ∪ D₂ ∪ D₃

D₁ = Point{Float64}(0.0) ∪ Point{Float64}(1.0)
D₂ = (0.0..1.0 (Unit)) × D₁ × (0.0..1.0 (Unit))
D₃ = (0.0..1.0 (Unit)) × (0.0..1.0 (Unit)) × D₁
D₄ = D₁ × (0.0..1.0 (Unit)) × (0.0..1.0 (Unit))

This is not wrong, but it may be better to support hyperfaces, so that the boundary of all polytopes can be a collection of hyperfaces (of the same type).

[daanhb]

Closing in favour of #100