diff --git a/Project.toml b/Project.toml
index 002e7b4..67c8038 100644
--- a/Project.toml
+++ b/Project.toml
@@ -11,12 +11,13 @@ Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 Meshes = "eacbb407-ea5a-433e-ab97-5258b1ca43fa"
 PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc"
 
 [compat]
 DataFrames = "1"
 JuMP = "1.4"
-Meshes = "0.35,0.36"
+Meshes = "0.39,0.40"
 PrettyTables = "2"
 Rotations = "1.5"
 julia = "1.9"
diff --git a/docs/src/example.md b/docs/src/example.md
index ad766ab..b2c0c22 100644
--- a/docs/src/example.md
+++ b/docs/src/example.md
@@ -122,19 +122,19 @@ Note that a plane normal is given, that is only used later when visualizing the
 ## Rough geometry definitions
 
 fronthole_r = SimpleHole([82.5, 30, 40], 26)
-frontface_r = PlaneAndNormal([82.5, 30, 40], [1, 0, 0])
+frontface_r = SimplePlane([82.5, 30, 40], [1, 0, 0])
 
 righthole1_r = SimpleHole([66, 71.5, 55], 6)
 righthole2_r = SimpleHole([58, 74.5, 24], 4.905)
 righthole3_r = SimpleHole([21.5, 68.5, 40], 8)
-rightface1_r = PlaneAndNormal([66, 71.5, 55], [0, 1, 0])
-rightface2_r = PlaneAndNormal([58, 74.5, 24], [0, 1, 0])
-rightface3_r = PlaneAndNormal([21.5, 68.5, 40], [0, 1, 0])
+rightface1_r = SimplePlane([66, 71.5, 55], [0, 1, 0])
+rightface2_r = SimplePlane([58, 74.5, 24], [0, 1, 0])
+rightface3_r = SimplePlane([21.5, 68.5, 40], [0, 1, 0])
 
 backhole1_r = SimpleHole([-3, 44, 53.9], 6.2)
 backhole2_r = SimpleHole([-3, 16.1, 54], 6.25)
-backface1_r = PlaneAndNormal([-3, 44, 54], [-1, 0, 0])
-backface2_r = PlaneAndNormal([-3, 16, 54], [-1, 0, 0])
+backface1_r = SimplePlane([-3, 44, 54], [-1, 0, 0])
+backface2_r = SimplePlane([-3, 16, 54], [-1, 0, 0])
 ```
 
 ## Pairing the rough and machined features
diff --git a/examples/example.jl b/examples/example.jl
index 758e08c..fa7ef67 100644
--- a/examples/example.jl
+++ b/examples/example.jl
@@ -30,19 +30,19 @@ backface2_m = SimplePlane([14, 14, 0])
 ## Rough geometry definitions
 
 fronthole_r = SimpleHole([82.5, 30, 40], 26)
-frontface_r = PlaneAndNormal([82.5, 30, 40], [1, 0, 0])
+frontface_r = SimplePlane([82.5, 30, 40], [1, 0, 0])
 
 righthole1_r = SimpleHole([66, 71.5, 55], 6)
 righthole2_r = SimpleHole([58, 74.5, 24], 4.905)
 righthole3_r = SimpleHole([21.5, 68.5, 40], 8)
-rightface1_r = PlaneAndNormal([66, 71.5, 55], [0, 1, 0])
-rightface2_r = PlaneAndNormal([58, 74.5, 24], [0, 1, 0])
-rightface3_r = PlaneAndNormal([21.5, 68.5, 40], [0, 1, 0])
+rightface1_r = SimplePlane([66, 71.5, 55], [0, 1, 0])
+rightface2_r = SimplePlane([58, 74.5, 24], [0, 1, 0])
+rightface3_r = SimplePlane([21.5, 68.5, 40], [0, 1, 0])
 
 backhole1_r = SimpleHole([-3, 44, 53.9], 6.2)
 backhole2_r = SimpleHole([-3, 16.1, 54], 6.25)
-backface1_r = PlaneAndNormal([-3, 44, 54], [-1, 0, 0])
-backface2_r = PlaneAndNormal([-3, 16, 54], [-1, 0, 0])
+backface1_r = SimplePlane([-3, 44, 54], [-1, 0, 0])
+backface2_r = SimplePlane([-3, 16, 54], [-1, 0, 0])
 
 ## Geometry pairing and feature descriptors
 
diff --git a/src/BlankLocalizationCore.jl b/src/BlankLocalizationCore.jl
index 8d70613..50c8720 100644
--- a/src/BlankLocalizationCore.jl
+++ b/src/BlankLocalizationCore.jl
@@ -4,40 +4,40 @@ using JuMP
 using DataFrames: DataFrame, names, nrow
 using PrettyTables: pretty_table, ft_nonothing, tf_html_minimalist
 using Rotations: RotMatrix
+import Meshes
 using Meshes: SimpleMesh, vertices, boundingbox, connect, Point3, Vec3, Plane, Cylinder,
-    Rotate, Translate, Disk
+    Rotate, Translate, Disk, top, radius, normal, boundingbox
 using Logging: @warn
-using LinearAlgebra: norm, cross, normalize, normalize!
+using Statistics: mean
+using LinearAlgebra: norm, dot, cross, inv, normalize, normalize!
 using Printf: @sprintf
 
 export  PartZero,
         printpartzeropositions
 
-export  AbstractHoleGeometry,
-        AbstractPlaneGeometry,
-        GeometryStyle,
-        IsPrimitive,
-        IsFreeForm,
-        surfacepoints,
-        filteredsurfacepoints,
-        featurepoint,
-        featureradius,
-        visualizationgeometry,
+export  RepresentationStyle,
+        FeatureStyle,
+        AbstractLocalizationGeometry,
         SimpleHole,
-        SimplePlane,
-        PlaneAndNormal,
         MeshHole,
+        SimplePlane,
         MeshPlane,
-        FeatureDescriptor,
-        LocalizationFeature,
-        HoleLocalizationFeature,
-        PlaneLocalizationFeature,
-        localizationfeature,
-        OptimizationResult,
-        Tolerance,
+        featurepoint,
+        surfacepoints,
+        filteredsurfacepoints,
+        featureradius,
+        RoughFeature,
+        MachinedFeature,
+        LocalizationFeature
+        
+export  PositionTolerance,
+        ConcentrictyTolerance,
+        ProjectedDimensionTolerance,
+        addtolerance2model!
+        
+export  OptimizationResult,
         MultiOperationProblem,
-        setparameters!,
-        isoptimum
+        setparameters!
 
 export  createjumpmodel,
         setjumpresult!,
@@ -66,8 +66,10 @@ HV(v) = vcat(v, 1)
 
 include("partzeros.jl")
 include("geometries.jl")
+include("tolerances.jl")
+include("optimizationproblem.jl")
 include("optimization.jl")
-include("resultevaluation.jl")
+#include("resultevaluation.jl")
 include("visualization.jl")
 
 end
diff --git a/src/geometries.jl b/src/geometries.jl
index 1d58dfc..1500c8d 100644
--- a/src/geometries.jl
+++ b/src/geometries.jl
@@ -1,517 +1,235 @@
-"""Supertype for localization geometries."""
-abstract type AbstractLocalizationGeometry end
-
-"""Supertype of hole like localization geometries."""
-abstract type AbstractHoleGeometry <: AbstractLocalizationGeometry end
-
-"""Supertype of plane like geometries."""
-abstract type AbstractPlaneGeometry <: AbstractLocalizationGeometry end
-
-
-# Traits that all features need to have:
-# rough OR machined -> this can be eliminated by having two fields: rough and machined
-# primitive OR freeform
-# hole OR plane: trait/abstract type/type parameter???
-
-
 """
-    GeometryStyle
+    RepresentationStyle
 
 Trait that describes the "style" of an [`AbstractLocalizationGeometry`](@ref).
-Currently it can be either [`IsPrimitive`](@ref) or [`IsFreeForm`](@ref).
+If it can be described by a primitive, then it is [`Primitive`](@ref) and
+[`FreeForm`](@ref) otherwise.
 """
-abstract type GeometryStyle end
+abstract type RepresentationStyle end
 
 """Primitive geometries can be explicitly described, e.g. a box or sphere."""
-struct IsPrimitive <: GeometryStyle end
-"""Free form geometries are discrete representations, e.g. a mesh or a point cloud."""
-struct IsFreeForm <: GeometryStyle end
-
-# don't define a global default (currently it helps development and debugging)
-#GeometryStyle(::Type) = IsPrimitive()
-
-"""
-    visualizationgeometry(geom::AbstractLocalizationGeometry)
-
-Return a Meshes.jl object that can be visualized.
-"""
-function visualizationgeometry end
-
-"""
-    featurepoint()
-
-Return the feature point of an [`IsPrimitive`](@ref) geometry.
-Definition signature should look like: `featurepoint(::IsPrimitive, x)`.
-"""
-function featurepoint end
-
-"""
-    surfacepoints()
+struct Primitive <: RepresentationStyle end
 
-Return the points of the surface of an [`IsFreeForm`](@ref) geometry.
-Definition signature should look like: `surfacepoints(::IsFreeForm, x)`.
-"""
-function surfacepoints end
+"""Free form geometries are discrete representations, e.g. a mesh or a point cloud."""
+struct FreeForm <: RepresentationStyle end
 
-"""
-    filteredsurfacepoints()
+nicestr(x::Primitive) = "primitive"
+nicestr(x::FreeForm) = "free form"
 
-Return the filtered points of the surface of an [`IsFreeForm`](@ref) geometry, 
-that may define active constraints in the optimization task
-(for example convex hull of mesh).
-Definition signature should look like: `filteredsurfacepoints(::IsFreeForm, x)`.
-"""
-function filteredsurfacepoints end
+abstract type FeatureStyle end
+struct Planar <: FeatureStyle end
+struct Cylindrical <: FeatureStyle end
 
-"""
-    featureradius()
-
-Return the radius of a [`IsPrimitive`](@ref) geometry
-that is subtype of [`AbstractHoleGeometry`].
-There is a default implementation that can be used: `featureradius(::IsPrimitive, x) = x.r`.
-"""
-function featureradius end
+nicestr(x::Planar) = "planar"
+nicestr(x::Cylindrical) = "cylindrical"
 
-featurepoint(x::T) where {T} = featurepoint(GeometryStyle(T), x)
-featurepoint(::IsPrimitive, x) = x.p
-function featurepoint(::IsFreeForm, x)
-    error("Function `featurepoint` is not defined for `IsFreeForm`` features")
-end
+abstract type AbstractLocalizationGeometry end
 
+visualizationgeometry(g::AbstractLocalizationGeometry) = g.geom
 
-surfacepoints(x::T) where {T} = surfacepoints(GeometryStyle(T), x)
-function surfacepoints(::IsPrimitive, x)
-    error("Function `surfacepoints` is not defined for `IsFreeForm`` features")
+struct SimpleHole <: AbstractLocalizationGeometry
+    geom::Meshes.Disk
 end
 
+RepresentationStyle(g::SimpleHole) = Primitive()
+FeatureStyle(g::SimpleHole) = Cylindrical()
 
-filteredsurfacepoints(x::T) where {T} = filteredsurfacepoints(GeometryStyle(T), x)
-function filteredsurfacepoints(::IsPrimitive, x)
-    error("Function `surfacepoints` is not defined for `IsPrimitive`` features")
-end
+featurepoint(g::SimpleHole) = Meshes.plane(g.geom)(0,0)
+featureradius(g::SimpleHole) = radius(g.geom)
 
-featureradius(x::T) where {T<:AbstractHoleGeometry} = featureradius(GeometryStyle(T), x)
-featureradius(::IsPrimitive, x) = x.r
-function featureradius(::IsFreeForm, x)
-    error("Function `featureradius` is not defined for `IsFreeForm`` features")
+struct MeshHole <: AbstractLocalizationGeometry
+    geom #::Meshes.Mesh, but that is abstract
+    chull # vector of points / vector of vectors - TBD
 end
 
-"""
-SimpleHole <: AbstractHoleGeometry
+RepresentationStyle(g::MeshHole) = FreeForm()
+FeatureStyle(g::MeshHole) = Cylindrical()
 
-A simple "hole" structure with a center point and a radius.
-Axis of the hole is defined by its partzero taken from the feature descriptor.
-"""
-struct SimpleHole <: AbstractHoleGeometry
-    p::Vector{Float64}
-    r::Float64
-end
+surfacepoints(g::MeshHole) = vertices(g.geom)
+filteredsurfacepoints(g::MeshHole) = g.chull
 
-GeometryStyle(::Type{SimpleHole}) = IsPrimitive()
-
-function visualizationgeometry(hole::SimpleHole)
-    # p1: feature point
-    p1 = Point3(hole.p)
-    return Disk(Plane(p1, Vec3(0,0,1)), hole.r)
+struct SimplePlane <: AbstractLocalizationGeometry
+    geom::Meshes.Plane
 end
 
-"""
-SimplePlane <: AbstractPlaneGeometry
+RepresentationStyle(g::SimplePlane) = Primitive()
+FeatureStyle(g::SimplePlane) = Planar()
 
-A simple plane structure with one point.
-Normal vector of the plane is defined by its partzero taken from the feature descriptor.
-"""
-struct SimplePlane <: AbstractPlaneGeometry
-    p::Vector{Float64}
-end
+featurepoint(g::SimplePlane) = g.geom(0,0)
 
-GeometryStyle(::Type{SimplePlane}) = IsPrimitive()
+randnormal(v) = normalize(cross(v, rand(typeof(v))))
 
 function rectangleforplane(point, v1 ,v2, sidelength)
     c1 = point + sidelength/2*v1 + -1*sidelength/2*v2
     c2 = c1 + -1*sidelength*v1
     c3 = c2 + sidelength*v2
     c4 = c3 + sidelength*v1
-    g1 = Point3(c1)
-    g2 = Point3(c2)
-    g3 = Point3(c3)
-    g4 = Point3(c4)
-    return SimpleMesh([g1,g2,g3,g4], connect.([(1,2,3),(3,4,1)]))
+    return SimpleMesh([c1,c2,c3,c4], connect.([(1,2,3),(3,4,1)]))
 end
 
-function visualizationgeometry(plane::SimplePlane)
-    return rectangleforplane(plane.p, [1,0,0], [0,1,0], 20)
-end
-
-"""
-PlaneAndNormal <: AbstractPlaneGeometry
-
-A simple plane structure with one point and a normal vector.
-"""
-struct PlaneAndNormal <: AbstractPlaneGeometry
-    p::Vector{Float64}
-    n::Vector{Float64}
-end
-
-GeometryStyle(::Type{PlaneAndNormal}) = IsPrimitive()
-
-# should try first the 3 axes
-randnormal(v::Vector) = normalize(cross(v, rand(3)))
-
-function visualizationgeometry(plane::PlaneAndNormal)
-    o = plane.p
-    v1 = randnormal(plane.n)
-    v2 = cross(v1, plane.n)
+function rectangleforplane(plane::Meshes.Plane)
+    o = plane(0,0)
+    n = Meshes.normal(plane)
+    v1 = randnormal(n)
+    v2 = cross(v1, n)
     return rectangleforplane(o, v1, v2, 20)
 end
 
-"""
-MeshHole <: AbstractHoleGeometry
+visualizationgeometry(plane::SimplePlane) = rectangleforplane(plane.geom)
 
-A simple mesh hole geometry, that contains the mesh of the hole's surface and the convex
-hull of the points (see our paper for details).
-"""
-struct MeshHole <: AbstractHoleGeometry
-    surface::SimpleMesh
-    convexhull::Vector{Vector{Float64}}
+struct MeshPlane <: AbstractLocalizationGeometry
+    geom #::Meshes.Mesh, but that is abstract
 end
 
-GeometryStyle(::Type{MeshHole}) = IsFreeForm()
+RepresentationStyle(g::MeshPlane) = FreeForm()
+FeatureStyle(g::MeshPlane) = Planar()
 
-function visualizationgeometry(meshhole::MeshHole)
-    return meshhole.surface
-end
+surfacepoints(g::MeshPlane) = vertices(g.geom)
 
-function filteredsurfacepoints(::IsFreeForm, x::MeshHole)
-    return x.convexhull
+function filteredsurfacepoints(g::MeshPlane)
+    bb = boundingbox(vertices(g.geom))
+    return [[bb.min.coords.coords...], [bb.max.coords.coords...]]
 end
 
-"""
-MeshPlane <: AbstractPlaneGeometry
-
-A simple mesh plane geometry, that contains the mesh of a planar face.
-"""
-struct MeshPlane <: AbstractPlaneGeometry
-    surface::SimpleMesh
-end
-
-GeometryStyle(::Type{MeshPlane}) = IsFreeForm()
-
-function visualizationgeometry(meshplane::MeshPlane)
-    return meshplane.surface
-end
-
-function filteredsurfacepoints(::IsFreeForm, x::MeshPlane)
-    bbox = boundingbox(x.surface)
-    return [bbox.min.coords, bbox.max.coords]
-end
-
-
-"""
-Store description of a feature: its name, the corresponding part zero, if it has or has not
-a machined and a rough state.
-"""
-struct FeatureDescriptor
-    name::String
-    partzero::PartZero
-    hasmachined::Bool
-    hasrough::Bool
-end
-
-getfeaturename(f::FeatureDescriptor) = f.name
-getpartzero(f::FeatureDescriptor) = f.partzero
-getpartzeroname(f::FeatureDescriptor) = getpartzeroname(getpartzero(f))
-hasmachined(f::FeatureDescriptor) = f.hasmachined
-hasrough(f::FeatureDescriptor) = f.hasrough
-
-function Base.show(io::IO, fd::FeatureDescriptor)
-    print(io, fd.name, " in ", fd.partzero.name,
-    fd.hasmachined ? "; machined, " : "; ! machined, ",
-    fd.hasrough ? "rough" : "! rough")
-end
-
-"""
-    LocalizationFeature{R,M}
-
-Supertype of any localization features.
-A localization feature contains a feature
-descriptor ([`FeatureDescriptor`](@ref)) and a rough and machined geometry
-([`AbstractLocalizationGeometry`](@ref)).
-The two geometries must be of same type (hole, plane, etc.).
-If a feature doesn't have a rough of machined state, an empty object should be used
-(and the feature descriptor should also store this information).
-based on its rough geometry.
-"""
-abstract type LocalizationFeature{R,M} end
-
-function Base.show(io::IO, lf::LocalizationFeature)
-    print(io, typeof(lf), ": ", lf.descriptor.name)
-end
 
-"""
-    HoleLocalizationFeature(descriptor::FeatureDescriptor, rough::R, machined::M) where {R<:AbstractHoleGeometry,M<:AbstractHoleGeometry}
 
-A holelike localization feature. The rough and machined geometries don't necessarily
-have to be the same type.
-"""
-struct HoleLocalizationFeature{R<:AbstractHoleGeometry,M<:AbstractHoleGeometry} <: LocalizationFeature{R,M}
-    descriptor::FeatureDescriptor
-    rough::R
-    machined::M
-end
+"""Supertype for feature types."""
+abstract type AbstractFeature end
 
-"""
-    PlaneLocalizationFeature(descriptor::FeatureDescriptor, rough::R, machined::M) where {R<:AbstractPlaneGeometry,M<:AbstractPlaneGeometry}
+RepresentationStyle(f::AbstractFeature) = RepresentationStyle(f.geometry)
+FeatureStyle(f::AbstractFeature) = FeatureStyle(f.geometry)
+featurename(f::AbstractFeature) = f.name
+geometry(f::AbstractFeature) = f.geometry
 
-A planelike localization feature. The rough and machined geometries don't necessarily
-have to be the same type.
-"""
-struct PlaneLocalizationFeature{R<:AbstractPlaneGeometry,M<:AbstractPlaneGeometry} <: LocalizationFeature{R,M}
-    descriptor::FeatureDescriptor
-    rough::R
-    machined::M
-end
+isplanar(st::AbstractFeature) = FeatureStyle(st) === Planar()
+iscylindrical(st::AbstractFeature) = FeatureStyle(st) === Cylindrical()
 
 """
-    localizationfeature(descriptor::FeatureDescriptor, rough::AbstractHoleGeometry, machined::AbstractHoleGeometry)
+    featurepoint(f::AbstractFeature)
 
-Convenience constructor for [`HoleLocalizationFeature`](@ref).
+Return the feature point of an [`Primitive`](@ref) geometry.
 """
-function localizationfeature(descriptor, rough::AbstractHoleGeometry, machined::AbstractHoleGeometry)
-    return HoleLocalizationFeature(descriptor, rough, machined)
+featurepoint(f::AbstractFeature) = featurepoint(RepresentationStyle(f), f)
+featurepoint(::Primitive, f) = featurepoint(geometry(f))
+function featurepoint(::FreeForm, f)
+    error("Function `featurepoint` is not defined for `FreeForm` features.")
 end
 
 """
-    localizationfeature(descriptor::FeatureDescriptor, rough::AbstractPlaneGeometry, machined::AbstractPlaneGeometry)
+    surfacepoints(f::AbstractFeature)
 
-Convenience constructor for [`PlaneLocalizationFeature`](@ref).
+Return the points of the surface of a [`FreeForm`](@ref) geometry.
 """
-function localizationfeature(descriptor, rough::AbstractPlaneGeometry, machined::AbstractPlaneGeometry)
-    return PlaneLocalizationFeature(descriptor, rough, machined)
-end
+surfacepoints(f::AbstractFeature) = surfacepoints(RepresentationStyle(f), f)
 
-getfeaturename(f::LocalizationFeature) = getfeaturename(f.descriptor)
-getpartzero(f::LocalizationFeature) = getpartzero(f.descriptor)
-getpartzeroname(f::LocalizationFeature) = getpartzeroname(f.descriptor)
-hasmachined(f::LocalizationFeature) = hasmachined(f.descriptor)
-hasrough(f::LocalizationFeature) = hasrough(f.descriptor)
+surfacepoints(::FreeForm, f) = surfacepoints(geometry(f))
 
-getroughfeaturepoint(f::LocalizationFeature) = featurepoint(f.rough)
-getmachinedfeaturepoint(f::LocalizationFeature) = featurepoint(f.machined)
-getmachinedradius(f::LocalizationFeature) = featureradius(f.machined)
-getroughradius(f::LocalizationFeature) = featureradius(f.rough)
-
-getroughfilteredpoints(f::LocalizationFeature) = filteredsurfacepoints(f.rough)
-
-function getmachinedfeaturepointindatum(f::LocalizationFeature)
-    @assert hasmachined(f)
-    v = getmachinedfeaturepoint(f)
-    pz = getpartzero(f)
-    T = getpartzeroHM(pz)
-    v_indatum = T*HV(v)
-    return v_indatum[1:3]
+function surfacepoints(::Primitive, f)
+    error("Function `surfacepoints` is not defined for `Primitive` features.")
 end
 
 """
-    OptimizationResult
-
-Store the status (result) of an optimization run and the minimum allowance value.
-"""
-struct OptimizationResult
-    status::String
-    minallowance::Float64
-end
-
-function Base.show(io::IO, or::OptimizationResult)
-    print(io, or.status, ", minimum allowance: ", or.minallowance)
-end
-
-emptyor() = OptimizationResult("empty", 0.0)
+    filteredsurfacepoints(f::AbstractFeature)
 
+Return the filtered points of the surface of a [`FreeForm`](@ref) geometry.
 """
-    isoptimum(or::OptimizationResult)
+filteredsurfacepoints(f::AbstractFeature) = filteredsurfacepoints(RepresentationStyle(f), f)
 
-Tell if `or` is in an optimal solution state, either: `OPTIMAL` or `LOCALLY_SOLVED`.
-"""
-function isoptimum(or::OptimizationResult)
-    return (or.status == "OPTIMAL") | (or.status == "LOCALLY_SOLVED")
-end
+filteredsurfacepoints(::FreeForm, f) = filteredsurfacepoints(geometry(f))
 
-struct Tolerance
-    featurename1::String
-    ismachined1::Bool
-    projection::Function
-    featurename2::String
-    ismachined2::Bool
-    nominalvalue::Float64
-    lowervalue::Float64
-    uppervalue::Float64
-    note::String
+function filteredsurfacepoints(::Primitive, f)
+    error("Function `surfacepoints` is not defined for `Primitive` features.")
 end
 
 """
-    MultiOperationProblem
+    featureradius(f::AbstractFeature)
 
-Collect all data for a multi operation problem, including: part zeros, holes, planes,
-tolerances, parameters and optimization result.
+Return the radius of a [`Primitive`](@ref) geometry, that is `Cylindrical`.
 """
-mutable struct MultiOperationProblem
-    partzeros::Vector{PartZero}
-    holes::Vector{HoleLocalizationFeature}
-    planes::Vector{PlaneLocalizationFeature}
-    tolerances::Vector{Tolerance}
-    parameters::Dict{String,Any}
-    opresult::OptimizationResult
-end
-
-"""
-    MultiOperationProblem(partzeros, holes, planes, tolerances, parameters)
-
-Construct a multi operation problem.
-For usage, please see the example section in the documentation.
-The parameters for the optimization are also described there with greater details.
-
-# Arguments
+featureradius(f::AbstractFeature) = featureradius(RepresentationStyle(f), FeatureStyle(f), f)
 
-- `partzeros::Vector{PartZero}`: array of part zeros.
-- `holes::Vector{HoleLocalizationFeature}`: array of holes.
-- `planes::Vector{PlaneLocalizationFeature}`: array of planes.
-- `tolerances::Vector{Tolerance}`: array of tolerances.
-- `parameters::Dict{String,Any}`: parameters in the form of a dictionary. Keys include:
-    `minAllowance`, `OptimizeForToleranceCenter`, `UseTolerances`,
-    `SetPartZeroPosition`, `maxPlaneZAllowance`.
-"""
-function MultiOperationProblem(partzeros, holes, planes, tolerances, parameters)
-    return MultiOperationProblem(partzeros, holes, planes, tolerances, parameters, emptyor())
-end
+featureradius(::Primitive, ::Cylindrical, f) = featureradius(geometry(f))
 
-function problemtype(mop::MultiOperationProblem)
-    # problem type is depending on the rough geometries: IsPrimitive or IsFreeForm
-    # if there is at least one IsFreeForm rough geometry -> hybrid problem
-    holetypes = GeometryStyle.(typeof.(x.rough for x in mop.holes))
-    for ht in holetypes
-        ht === IsFreeForm() && return :HybridProblem
-    end
-    planetypes = GeometryStyle.(typeof.(x.rough for x in mop.planes))
-    for pt in planetypes
-        pt === IsFreeForm() && return :HybridProblem
-    end
-    return :PrimitiveProblem
+function featureradius(t1, t2, f)
+    error("Function `featureradius` is only defined for ::Primitive and ::Cylindrical features. Got $t1 and $t2")
 end
 
-function Base.show(io::IO, mop::MultiOperationProblem)
-    nh = size(mop.holes, 1)
-    np = size(mop.planes, 1)
-    npz = size(mop.partzeros, 1)
-    nts = size(mop.tolerances, 1)
-    sn = string(problemtype(mop))
-    print(io, sn,": ",
-    npz," part zero", npz > 1 ? "s, " : ", ",
-    nh," hole", nh > 1 ? "s, " : ", ",
-    np," plane", np > 1 ? "s, " : ", ",
-    nts," tolerance", nts > 1 ? "s" : "",
-    ", status: ", mop.opresult.status)
+struct RoughFeature <: AbstractFeature
+    name::String
+    geometry # <: AbstractLocalizationGeometry
 end
 
-"""
-    printpartzeropositions(mop::MultiOperationProblem)
-
-Print the positions of the part zeros of a `MultiOperationProblem`.
-"""
-printpartzeropositions(mop::MultiOperationProblem) = printpartzeropositions(mop.partzeros)
-
-"""
-    setparameters!(mop::MultiOperationProblem, pardict)
-
-Set parameter dictionary of a `MultiOperationProblem` to `pardict`.
-"""
-function setparameters!(mop::MultiOperationProblem, pardict)
-    mop.parameters = pardict
-    return mop
+struct MachinedFeature <: AbstractFeature
+    name::String
+    geometry # <: AbstractLocalizationGeometry
+    partzero::PartZero
 end
 
-"""
-    isoptimum(mop::MultiOperationProblem)
+partzero(f::MachinedFeature) = f.partzero
+partzeroname(f::MachinedFeature) = partzeroname(f.partzero)
 
-Tell if `mop`'s solution is in an optimal state, either: `OPTIMAL` or `LOCALLY_SOLVED`.
-"""
-isoptimum(mop::MultiOperationProblem) = isoptimum(mop.opresult)
 
 """
-    getfeaturebyname(mop::MultiOperationProblem, featurename)
+    LocalizationFeature
 
-Get a hole or plane feature by its name.
-It is assumed that all features have distinct names.
-Return `nothing`, if no feature is found with `featurename`.
+A feature that is machined and allowance can be computed for it.
+It has a name, a [`RoughFeature`](@ref) and a [`MachinedFeature`](@ref).
+The two geometries should be of same [`GeometryType`](@ref) (holelike or planelike),
+the constructor enforces this property.
 """
-function getfeaturebyname(mop::MultiOperationProblem, featurename)
-    function retbyname(array, name)
-        for f in array
-            if getfeaturename(f) == name
-                return f
-            end
+struct LocalizationFeature
+    name::String
+    roughfeature::RoughFeature
+    machinedfeature::MachinedFeature
+    function LocalizationFeature(n, r, m)
+        if FeatureStyle(r) === FeatureStyle(m)
+            new(n, r, m)
+        else
+            error("FeatureStyle of rough and machined features does not match!")
         end
-        return nothing
     end
-
-    hole_ = retbyname(mop.holes, featurename)
-    isnothing(hole_) || return hole_
-    # return plane even if it is nothing
-    return retbyname(mop.planes, featurename)
 end
 
-"""
-    collectholesbypartzero(mop::MultiOperationProblem, partzeroname)
+featurename(lf::LocalizationFeature) = lf.name
+partzero(lf::LocalizationFeature) = partzero(lf.machinedfeature)
+partzeroname(lf::LocalizationFeature) = partzeroname(partzero(lf))
+isplanar(lf::LocalizationFeature) = isplanar(lf.roughfeature)
+iscylindrical(lf::LocalizationFeature) = iscylindrical(lf.roughfeature)
+RepresentationStyle(lf::LocalizationFeature) = RepresentationStyle(lf.roughfeature)
+FeatureStyle(lf::LocalizationFeature) = FeatureStyle(lf.roughfeature)
 
-Collect holes that are grouped to part zero called `partzeroname`.
-"""
-function getholesbypartzero(mop::MultiOperationProblem, partzeroname)
-    return filter(x->getpartzeroname(x) == partzeroname, mop.holes)
+function Base.show(io::IO, lf::LocalizationFeature)
+    print(io, "LocalizationFeature: ", featurename(lf),
+    " ", nicestr(RepresentationStyle(lf)), ", ", nicestr(FeatureStyle(lf)))
 end
 
-"""
-    collectmachinedholes(mop::MultiOperationProblem)
+roughfeaturepoint(lf::LocalizationFeature) = featurepoint(lf.roughfeature)
+machinedfeaturepoint(lf::LocalizationFeature) = featurepoint(lf.machinedfeature)
+roughradius(lf::LocalizationFeature) = featureradius(lf.roughfeature)
+machinedradius(lf::LocalizationFeature) = featureradius(lf.machinedfeature)
+machinedfilteredsurfacepoints(lf::LocalizationFeature) = filteredsurfacepoints(lf.machinedfeature)
+roughfilteredsurfacepoints(lf::LocalizationFeature) = filteredsurfacepoints(lf.roughfeature)
 
-Collect holes that have a machined state.
-"""
-collectmachinedholes(mop::MultiOperationProblem) = filter(hasmachined, mop.holes)
-
-"""
-    collectmachinedplanes(mop::MultiOperationProblem)
-
-Collect planes that have a machined state.
-"""
-collectmachinedplanes(mop::MultiOperationProblem) = filter(hasmachined, mop.planes)
-
-"""
-    collectroughholes(mop::MultiOperationProblem)
-
-Collect those holes, that have rough stage.
-"""
-collectroughholes(mop::MultiOperationProblem) = filter(hasrough, mop.holes)
-
-"""
-    collectroughplanes(mop::MultiOperationProblem)
-
-Collect those planes, that have rough stage.
-"""
-collectroughplanes(mop::MultiOperationProblem) = filter(hasrough, mop.planes)
-
-
-"""
-    collectallowancedholes(mop::MultiOperationProblem)
-
-Collect holes that have machined and rough states, thus allowance should be calculated for.
-"""
-function collectallowancedholes(mop::MultiOperationProblem)
-    return [h for h in mop.holes if hasmachined(h) & hasrough(h)]
+#=
+function machinedfeaturepointindatum(f::LocalizationFeature)
+    v = machinedfeaturepoint(f)
+    pz = partzero(f)
+    T = getpartzeroHM(pz)
+    v_indatum = T*HV(v)
+    return v_indatum[1:3]
 end
 
 """
-    collectallowancedplanes(mop::MultiOperationProblem)
+    transformmachined2datum(feature, points)
 
-Collect planes that have machined and rough states, thus allowance should be calculated for.
+Transform a list of points with the part zero of `feature`.
 """
-function collectallowancedplanes(mop::MultiOperationProblem)
-    return [p for p in mop.planes if hasmachined(p) & hasrough(p)]
+function transformmachined2datum(feature, points)
+    pz = partzero(feature)
+    M = getpartzeroHM(pz)
+    newpoints = (M*HV(p) for p in points)
+    resultpoints = [p[1:3] for p in newpoints]
+    return resultpoints
 end
+=#
+
diff --git a/src/optimization.jl b/src/optimization.jl
index 5bc41e5..bd8d92f 100644
--- a/src/optimization.jl
+++ b/src/optimization.jl
@@ -1,20 +1,23 @@
-# dispatch on GeometryStyle trait
-function addhole2model!(model, hole::HoleLocalizationFeature{R,M}, ipzmatricedict) where {R,M}
-    return addhole2model!(GeometryStyle(R), model, hole, ipzmatricedict)
+function addallowancedfeature2model!(model, f, ipzmatricedict)
+    # dispatch on: primitive/free-form and planar/cylindrical
+    addallowancedfeature2model!(RepresentationStyle(f), FeatureStyle(f), model, f, ipzmatricedict)
 end
 
-function addhole2model!(::IsPrimitive, model, hole, ipzmatricedict)
+function addallowancedfeature2model!(::Primitive, ::Cylindrical, model, feature, ipzmatricedict)
     # access registered variables
     minAllowance = model[:minAllowance]
 
     # register distance variable:
-    dxy = @variable(model, base_name = string("d_xy_", getfeaturename(hole)), lower_bound = 0.0)
+    dxy = @variable(model, base_name = string("d_xy_", featurename(feature)), lower_bound = 0.0)
 
-    pzn = getpartzeroname(hole)
-    v_machined = getmachinedfeaturepoint(hole)
-    v_rough = getroughfeaturepoint(hole)
-    r_machined = getmachinedradius(hole)
-    r_rough = getroughradius(hole)
+    pzn = partzeroname(feature)
+    # featurepoints are Meshes.Point
+    v_machined_ = machinedfeaturepoint(feature)
+    v_machined = [v_machined_.coords.coords...]
+    v_rough_ = roughfeaturepoint(feature)
+    v_rough = [v_rough_.coords.coords...]
+    r_machined = machinedradius(feature)
+    r_rough = roughradius(feature)
     # equation (4)
     d_f = @expression(model, HV(v_machined)-ipzmatricedict[pzn]*HV(v_rough))
     # equation (5)
@@ -24,72 +27,75 @@ function addhole2model!(::IsPrimitive, model, hole, ipzmatricedict)
     return model
 end
 
-function addhole2model!(::IsFreeForm, model, hole, ipzmatricedict)
+function addallowancedfeature2model!(::Primitive, ::Planar, model, feature, ipzmatricedict)
     # access registered variables
     minAllowance = model[:minAllowance]
-
-    # filtered surface points of a free form surface
-    qs = getroughfilteredpoints(hole)
-    qiter = 1:length(qs)
+    maxPlaneZAllowance = model[:maxPlaneZAllowance]
 
     # register distance variable:
-    dxy = @variable(model, [qiter], base_name = string("d_xy_", getfeaturename(hole)), lower_bound = 0.0)
+    dz = @variable(model, base_name = string("d_z_", featurename(feature)))
+    
+    pzn = partzeroname(feature)
+    v_machined_ = machinedfeaturepoint(feature)
+    v_machined = [v_machined_.coords.coords...]
+    v_rough_ = roughfeaturepoint(feature)
+    v_rough = [v_rough_.coords.coords...]
 
-    pzn = getpartzeroname(hole)
-    v_machined = getmachinedfeaturepoint(hole)
-    r_machined = getmachinedradius(hole)
     # equation (4)
-    for (i, q) in enumerate(qs)
-        d_f = @expression(model, HV(v_machined)-ipzmatricedict[pzn]*HV(q))
-        # equation (5)
-        @constraint(model, dxy[i]*dxy[i] >= d_f[1]*d_f[1] + d_f[2]*d_f[2])
-        # equation (6)
-        @constraint(model, r_machined - dxy[i] >= minAllowance)
-    end
+    d_f = @expression(model, HV(v_machined)-ipzmatricedict[pzn]*HV(v_rough))
+    @constraint(model, dz == d_f[3])
+    # equation (7)
+    @constraint(model, -1*dz >= minAllowance)
+    @constraint(model, -1*dz <= maxPlaneZAllowance)
     return model
 end
 
-# dispatch on GeometryStyle trait
-function addplane2model!(model, plane::PlaneLocalizationFeature{R,M}, ipzmatricedict) where {R,M}
-    return addplane2model!(GeometryStyle(R), model, plane, ipzmatricedict)
-end
-
-function addplane2model!(::IsPrimitive, model, plane, ipzmatricedict)
+function addallowancedfeature2model!(::FreeForm, ::Cylindrical, model, feature, ipzmatricedict)
     # access registered variables
     minAllowance = model[:minAllowance]
-    maxPlaneZAllowance = model[:maxPlaneZAllowance]
-
+    
+    # filtered surface points of a free form surface
+    @show feature
+    qs = roughfilteredsurfacepoints(feature)
+    qiter = 1:length(qs)
+    
     # register distance variable:
-    dz = @variable(model, base_name = string("d_z_", getfeaturename(plane)))
+    dxy = @variable(model, [qiter], base_name = string("d_xy_", featurename(feature)), lower_bound = 0.0)
     
-    pzn = getpartzeroname(plane)
-    v_machined = getmachinedfeaturepoint(plane)
-    v_rough = getroughfeaturepoint(plane)
+    pzn = partzeroname(feature)
+    v_machined_ = machinedfeaturepoint(feature)
+    v_machined = [v_machined_.coords.coords...]
+    r_machined = machinedradius(feature)
     # equation (4)
-    d_f = @expression(model, HV(v_machined)-ipzmatricedict[pzn]*HV(v_rough))
-    @constraint(model, dz == d_f[3])
-    # equation (7)
-    @constraint(model, -1*dz >= minAllowance)
-    @constraint(model, -1*dz <= maxPlaneZAllowance)
+    for (i, q) in enumerate(qs)
+        #q = [qmeshes.coords.coords...]
+        d_f = @expression(model, HV(v_machined)-ipzmatricedict[pzn]*HV(q))
+        # equation (5)
+        @constraint(model, dxy[i]*dxy[i] >= d_f[1]*d_f[1] + d_f[2]*d_f[2])
+        # equation (6)
+        @constraint(model, r_machined - dxy[i] >= minAllowance)
+    end
     return model
 end
 
-function addplane2model!(::IsFreeForm, model, plane, ipzmatricedict)
+function addallowancedfeature2model!(::FreeForm, ::Planar, model, feature, ipzmatricedict)
     # access registered variables
     minAllowance = model[:minAllowance]
     maxPlaneZAllowance = model[:maxPlaneZAllowance]
 
     # filtered surface points of a free form surface
-    qs = getroughfilteredpoints(plane)
+    qs = roughfilteredsurfacepoints(feature)
     qiter = 1:length(qs)
 
     # register distance variable:
-    dz = @variable(model, [qiter], base_name = string("d_z_", getfeaturename(plane)))
+    dz = @variable(model, [qiter], base_name = string("d_z_", featurename(feature)))
 
-    pzn = getpartzeroname(plane)
-    v_machined = getmachinedfeaturepoint(plane)
+    pzn = partzeroname(feature)
+    v_machined_ = machinedfeaturepoint(feature)
+    v_machined = [v_machined_.coords.coords...]
     # equation (4)
     for (i, q) in enumerate(qs)
+        #q = [qmeshes.coords.coords...]
         d_f = @expression(model, HV(v_machined)-ipzmatricedict[pzn]*HV(q))
         # equation (5)
         @constraint(model, dz[i] == d_f[3])
@@ -102,31 +108,39 @@ end
 
 function addtolerances2model!(model, mop::MultiOperationProblem, pzmatricedict)
     for (i, t) in enumerate(mop.tolerances)
+        function fpointexpr(feature::RoughFeature)
+            v_ = featurepoint(feature)
+            v = [v_.coords.coords...]
+            return @expression(model, v)
+        end
+        function fpointexpr(feature::MachinedFeature)
+            pzn = partzeroname(feature)
+            v_ = featurepoint(feature)
+            v = [v_.coords.coords...]
+            return @expression(model, pzmatricedict[pzn]*HV(v))
+        end
+
         # access registered variables
         AbsValRelError = model[:AbsValRelError]
 
-        # get features and part zero names
-        f1 = getfeaturebyname(mop, t.featurename1)
-        f2 = getfeaturebyname(mop, t.featurename2)
-        @assert ! isnothing(f1) "Feature $(t.featurename1) does not exist!"
-        @assert ! isnothing(f2) "Feature $(t.featurename1) does not exist!"
-        pzn1 = getpartzeroname(f1)
-        pzn2 = getpartzeroname(f2)
-
         # equation (2)
-        v1 = @expression(model, t.ismachined1 ? pzmatricedict[pzn1]*HV(getmachinedfeaturepoint(f1)) : getroughfeaturepoint(f1))
-        v2 = @expression(model, t.ismachined2 ? pzmatricedict[pzn2]*HV(getmachinedfeaturepoint(f2)) : getroughfeaturepoint(f2))
+        v1 = fpointexpr(datumfeature(t))
+        v2 = fpointexpr(tolerancefeature(t))
         e_t = @expression(model, v1[1:3]-v2[1:3])
-        real_d = @expression(model, t.projection(e_t))
+        pa = projectionaxis(t)
+        real_d = @expression(model, dot(pa, e_t))
 
         # tolerance should be in interval if set to be used
+        # get lower and upper values
+        lv = minimumvalue(t)
+        uv = maximumvalue(t)
         if mop.parameters["UseTolerances"]
             # equation (3)
-            @constraint(model, t.lowervalue <= real_d <= t.uppervalue)
+            @constraint(model, lv <= real_d <= uv)
         end
         # equation (1)
-        abs_dist = @expression(model, real_d - (t.lowervalue+t.uppervalue)/2)
-        rel_dist = @expression(model, 2*abs_dist/(t.uppervalue-t.lowervalue))
+        abs_dist = @expression(model, real_d - (lv+uv)/2)
+        rel_dist = @expression(model, 2*abs_dist/(uv-lv))
         @constraint(model, rel_dist <= AbsValRelError[i])
         @constraint(model, rel_dist >= -1*AbsValRelError[i])
     end
@@ -157,10 +171,6 @@ function createjumpmodel(mop::MultiOperationProblem, optimizer; disable_string_n
     pzmatrices = [vcat(hcat(pz.rotation, pzpose[i,1:3]), [0 0 0 1]) for (i, pz) in enumerate(partzeros)]
     pzmatricedict = Dict([(partzeros[i].name, pzmatrices[i]) for i in pzr])
 
-    # holes that have machined state
-    machinedholes = collectallowancedholes(mop)
-    # planes that have machined state
-    machinedplanes = collectallowancedplanes(mop)
     # numer of tolerances
     ntolerances = length(mop.tolerances)
     
@@ -175,14 +185,11 @@ function createjumpmodel(mop::MultiOperationProblem, optimizer; disable_string_n
     ## tolerances
     addtolerances2model!(model, mop, pzmatricedict)
 
-    # allowance for holes
-    for h in machinedholes
-        addhole2model!(model, h, ipzmatricedict)
-    end
-    # allowance for planes
-    for p in machinedplanes
-        addplane2model!(model, p, ipzmatricedict)
+    # allowances
+    for f in mop.features
+        addallowancedfeature2model!(model, f, ipzmatricedict)
     end
+
     # if maximum allowance of planes is given, then set it
     if haskey(mop.parameters, "maxPlaneZAllowance")
         @assert mop.parameters["maxPlaneZAllowance"] > mop.parameters["minAllowance"] "Maximum plane z allowance must be larger, than minimum allowance!"
diff --git a/src/optimizationproblem.jl b/src/optimizationproblem.jl
new file mode 100644
index 0000000..fa86323
--- /dev/null
+++ b/src/optimizationproblem.jl
@@ -0,0 +1,134 @@
+"""
+    OptimizationResult
+
+Store the status (result) of an optimization run and the minimum allowance value.
+"""
+struct OptimizationResult
+    status::String
+    minallowance::Float64
+end
+
+function Base.show(io::IO, or::OptimizationResult)
+    print(io, or.status, ", minimum allowance: ", or.minallowance)
+end
+
+emptyor() = OptimizationResult("empty", 0.0)
+
+"""
+    isoptimum(or::OptimizationResult)
+
+Tell if `or` is in an optimal solution state, either: `OPTIMAL` or `LOCALLY_SOLVED`.
+"""
+function isoptimum(or::OptimizationResult)
+    return (or.status == "OPTIMAL") | (or.status == "LOCALLY_SOLVED")
+end
+
+"""
+    MultiOperationProblem
+
+A type storing all geometry data and parameters to generate a JuMP model.
+"""
+mutable struct MultiOperationProblem
+    partzeros::Vector{PartZero}
+    features::Vector{LocalizationFeature} # features that have rough and machined -> allowanced
+    tolerances # Vector{AbstractTolerance}
+    parameters::Dict{String,Real}
+    opresult::OptimizationResult
+end
+
+function MultiOperationProblem(partzeros, features, tolerances, parameters)
+    return MultiOperationProblem(partzeros, features, tolerances, parameters, emptyor())
+end
+
+function problemtype(mop::MultiOperationProblem)
+    # problem type is depending on the rough geometries: Primitive or FreeForm
+    # if there is at least one FreeForm rough geometry -> hybrid problem
+    featuretypes = RepresentationStyle.(mop.features)
+    for ft in featuretypes
+        ft === FreeForm() && return :HybridProblem
+    end
+    return :PrimitiveProblem
+end
+
+function cylindricalfeatures(mop::MultiOperationProblem)
+    filter(x->FeatureStyle(x) === Cylindrical(), mop.features)
+end
+
+function planarfeatures(mop::MultiOperationProblem)
+    filter(x->FeatureStyle(x) === Planar(), mop.features)
+end
+
+function Base.show(io::IO, mop::MultiOperationProblem)
+    nh = size(cylindricalfeatures(mop), 1)
+    np = size(planarfeatures(mop), 1)
+    npz = size(mop.partzeros, 1)
+    nts = size(mop.tolerances, 1)
+    sn = string(problemtype(mop))
+    print(io, sn,": ",
+    npz," part zero", npz > 1 ? "s, " : ", ",
+    nh," hole", nh > 1 ? "s, " : ", ",
+    np," plane", np > 1 ? "s, " : ", ",
+    nts," tolerance", nts > 1 ? "s" : "",
+    ", status: ", mop.opresult.status)
+end
+
+"""
+    printpartzeropositions(mop::MultiOperationProblem)
+
+Print the positions of the part zeros of a `MultiOperationProblem`.
+"""
+printpartzeropositions(mop::MultiOperationProblem) = printpartzeropositions(mop.partzeros)
+
+"""
+    setparameters!(mop::MultiOperationProblem, pardict)
+
+Set parameter dictionary of a `MultiOperationProblem` to `pardict`.
+"""
+function setparameters!(mop::MultiOperationProblem, pardict)
+    mop.parameters = pardict
+    return mop
+end
+
+"""
+    isoptimum(mop::MultiOperationProblem)
+
+Tell if `mop`'s solution is in an optimal state, either: `OPTIMAL` or `LOCALLY_SOLVED`.
+"""
+isoptimum(mop::MultiOperationProblem) = isoptimum(mop.opresult)
+
+
+"""
+    getfeaturebyname(features, featuresname)
+
+Get a hole or plane feature by its name.
+It is assumed that all features have distinct names.
+Return `nothing`, if no feature is found with `featuresname`.
+"""
+function getfeaturebyname(features, featuresname)
+    for f in features
+        if featurename(f) == featuresname
+            return f
+        end
+    end
+    return nothing
+end
+
+"""
+    getfeaturebyname(mop::MultiOperationProblem, featurename)
+
+Get a hole or plane feature by its name from a vector of features.
+It is assumed that all features have distinct names.
+Return `nothing`, if no feature is found with `featurename`.
+"""
+function getfeaturebyname(mop::MultiOperationProblem, featurename)
+    getfeaturebyname(mop.features, featurename)
+end
+
+"""
+    collectfeaturesbypartzero(mop::MultiOperationProblem, partzeroname)
+
+Collect features that are grouped to part zero called `partzeroname`.
+"""
+function collectfeaturesbypartzero(mop::MultiOperationProblem, partzeroname)
+    return filter(x->getpartzeroname(x) == partzeroname, mop.features)
+end
diff --git a/src/partzeros.jl b/src/partzeros.jl
index 1cfe860..60d6c59 100644
--- a/src/partzeros.jl
+++ b/src/partzeros.jl
@@ -21,13 +21,14 @@ end
 
 Base.show(io::IO, pz::PartZero) = print(io, "Part zero: \"", pz.name, "\"")
 
-Base.show(io::IO, ::MIME"text/plain", pz::PartZero) = print(io, "Part zero: \"", pz.name, "\"\n", getpartzeroHM(pz))
+Base.show(io::IO, ::MIME"text/plain", pz::PartZero) = print(io, "Part zero: \"",
+    pz.name, "\"\n", getpartzeroHM(pz))
 
 xaxis(partzero::PartZero) = partzero.rotation[:,1]
 yaxis(partzero::PartZero) = partzero.rotation[:,2]
 zaxis(partzero::PartZero) = partzero.rotation[:,3]
 
-getpartzeroname(pz::PartZero) = pz.name
+partzeroname(pz::PartZero) = pz.name
 
 """
     getpartzeroHM(partzero::PartZero)
@@ -68,7 +69,6 @@ function inverthomtr(M)
     return invtr
 end
 
-
 """
     getpartzerobyname(partzeros::Vector{PartZero}, partzeroname::AbstractString)
 
@@ -81,7 +81,6 @@ function getpartzerobyname(partzeros::Vector{PartZero}, partzeroname::AbstractSt
     return nothing
 end
 
-
 """
     printpartzeropositions(partzeros::Vector{PartZero})
 
diff --git a/src/resultevaluation.jl b/src/resultevaluation.jl
index 13812bd..3f8c2a8 100644
--- a/src/resultevaluation.jl
+++ b/src/resultevaluation.jl
@@ -11,8 +11,8 @@ function computeallowance(::IsPrimitive, hole::HoleLocalizationFeature)
 
     if hasmachined(hole) & hasrough(hole)
         v_mlocal = getmachinedfeaturepoint(hole)
-        pz = getpartzero(hole)
-        T_inv = getpartzeroinverseHM(pz)
+        pz = partzero(hole)
+        T_inv = partzeroinverseHM(pz)
         d_f = HV(v_mlocal) - T_inv*HV(v_r)
         xydist = norm(d_f[1:2])
         rallowance = r_m - r_r - xydist
@@ -35,8 +35,8 @@ function computeallowance(::IsFreeForm, hole::HoleLocalizationFeature)
         qs = getroughfilteredpoints(hole)
         rallowances = zeros(Float64, size(qs))
         
-        pz = getpartzero(hole)
-        T_inv = getpartzeroinverseHM(pz)
+        pz = partzero(hole)
+        T_inv = partzeroinverseHM(pz)
         v_mlocal = getmachinedfeaturepoint(hole)
         for (i, q) in enumerate(qs)
             d_f = HV(v_mlocal) - T_inv*HV(q)
@@ -64,8 +64,8 @@ function computeallowance(::IsPrimitive, plane::PlaneLocalizationFeature)
 
     if hasmachined(plane) & hasrough(plane)
         v_mlocal = getmachinedfeaturepoint(plane)
-        pz = getpartzero(plane)
-        T_inv = getpartzeroinverseHM(pz)
+        pz = partzero(plane)
+        T_inv = partzeroinverseHM(pz)
         d_f = HV(v_mlocal) - T_inv*HV(v_r)
         zdist = d_f[3]
         axallowance = -1*zdist
@@ -86,8 +86,8 @@ function computeallowance(::IsFreeForm, plane::PlaneLocalizationFeature)
         qs = getroughfilteredpoints(plane)
         axallowances = zeros(Float64, size(qs))
         
-        pz = getpartzero(plane)
-        T_inv = getpartzeroinverseHM(pz)
+        pz = partzero(plane)
+        T_inv = partzeroinverseHM(pz)
         v_mlocal = getmachinedfeaturepoint(plane)
         for (i, q) in enumerate(qs)
             d_f = HV(v_mlocal) - T_inv*HV(q)
@@ -125,7 +125,7 @@ function allowancetable(mop::MultiOperationProblem)
         xydist = htuple.xydistance
         rallowance = htuple.rallowance
 
-        push!(df, [getfeaturename(h), getpartzeroname(h), v_m[1], v_m[2], v_m[3], v_r[1],
+        push!(df, [getfeaturename(h), partzeroname(h), v_m[1], v_m[2], v_m[3], v_r[1],
             v_r[2], v_r[3], r_m, r_r, xydist, nothing, rallowance, nothing])
     end
     # planes
@@ -137,7 +137,7 @@ function allowancetable(mop::MultiOperationProblem)
         axallowance = ptuple.axallowance
         
 
-        push!(df, [getfeaturename(p), getpartzeroname(p), v_m[1], v_m[2], v_m[3], v_r[1],
+        push!(df, [getfeaturename(p), partzeroname(p), v_m[1], v_m[2], v_m[3], v_r[1],
             v_r[2], v_r[3], nothing, nothing, nothing, zdist, nothing, axallowance])
     end
     return df
@@ -213,9 +213,9 @@ function tolerancetable(mop::MultiOperationProblem)
         # get features and part zero names
         fname1, fname2 = fnameplusrORm(t)
         f1 = getfeaturebyname(mop, t.featurename1)
-        pzn1 = getpartzeroname(f1)
+        pzn1 = partzeroname(f1)
         f2 = getfeaturebyname(mop, t.featurename2)
-        pzn2 = getpartzeroname(f2)
+        pzn2 = partzeroname(f2)
 
         v1 = t.ismachined1 ? getmachinedfeaturepointindatum(f1) : getroughfeaturepoint(f1)
         v2 = t.ismachined2 ? getmachinedfeaturepointindatum(f2) : getroughfeaturepoint(f2)
diff --git a/src/tolerances.jl b/src/tolerances.jl
new file mode 100644
index 0000000..d7fe4a9
--- /dev/null
+++ b/src/tolerances.jl
@@ -0,0 +1,350 @@
+"""Supertype of any tolerance."""
+abstract type AbstractTolerance end
+
+"""
+Function to add any tolerance 
+"""
+function addtolerance2model! end
+
+# getters -> standardize package: do I write get_*** or just ***?
+# these are the defaults
+"""
+    tolerancefeature(t::AbstractTolerance)
+
+Return the feature that is toleranced relative to the datum feature.
+"""
+tolerancefeature(t::AbstractTolerance) = t.feature
+
+"""
+    datumfeature(t::AbstractTolerance)
+
+Return the datum feature.
+"""
+datumfeature(t::AbstractTolerance) = t.datumfeature
+
+"""
+    nominalvalue(t::AbstractTolerance)
+
+Return the nominal value of a tolerance.
+"""
+nominalvalue(t::AbstractTolerance) = t.nominalval
+
+"""
+    minimumvalue(t::AbstractTolerance)
+
+Return the minimum value of a tolerance.
+"""
+minimumvalue(t::AbstractTolerance) = t.minval
+
+"""
+    maximumvalue(t::AbstractTolerance)
+
+Return the maximum value of a tolerance.
+"""
+maximumvalue(t::AbstractTolerance) = t.maxval
+
+"""Supertype of location type tolerances."""
+abstract type LocationTolerance <: AbstractTolerance end
+
+struct PositionTolerance <: LocationTolerance
+    datumfeature
+    feature
+    nominalval
+    minval
+    maxval
+    function PositionTolerance(d, f, nv, miv, mav)
+        if (RepresentationStyle(d) !== Primitive()) || (RepresentationStyle(f) !== Primitive())
+            error("PositionTolerance is currently only defined between primitives!")
+        end
+        return new(d, f, nv, miv, mav)
+    end
+end
+
+struct ConcentrictyTolerance <: LocationTolerance
+    datumfeature
+    feature
+    nominalval
+    minval
+    maxval
+end
+
+"""
+    ProjectedDimensionTolerance <: AbstractTolerance
+
+Dimension tolerance between a feature and a datum feature. Projection means
+that a projection function is used for R^3->R.
+"""
+struct ProjectedDimensionTolerance <: AbstractTolerance
+    datumfeature
+    feature
+    nominalval
+    minval
+    maxval
+    projectionaxis
+end
+
+function Base.show(io::IO, t::ProjectedDimensionTolerance)
+    print(io, "ProjDimTol<", featurename(t.datumfeature), "-", featurename(t.feature),
+    " : ", t.minval, "-", t.maxval, ">",)
+end
+
+"""
+    projectionaxis(t::ProjectedDimensionTolerance)
+
+Return the projection axis of `t`.
+"""
+projectionaxis(t::ProjectedDimensionTolerance) = t.projectionaxis
+
+
+#=
+
+"""Supertype of tolerance types."""
+abstract type LocalizationToleranceType end
+
+struct PlanePlaneDistance <: LocalizationToleranceType end
+
+struct PlaneAxisDistance <: LocalizationToleranceType end
+
+struct AxisAxisDistance <: LocalizationToleranceType
+    projectionaxis::Vector{Float64}
+end
+
+struct AxisAxisConcentric <: LocalizationToleranceType end
+
+"""
+    toleranceddistance
+
+Should have signature like:
+`toleranceddistance(type::LocalizationToleranceType, feature1, machined1, feature2, machined2)`.
+"""
+function toleranceddistance end
+
+"""
+    getfeaturepoints(x::T)
+
+Get points of a geometry. If it's `IsPrimitive`,
+then return the feature point wrapped in an array.
+If it's `IsFreeForm`, then return the surface points in an array.
+"""
+getfeaturepoints(x::T) where {T} = getfeaturepoints(GeometryStyle(T), x)
+getfeaturepoints(::IsPrimitive, x) = [featurepoint(x)]
+getfeaturepoints(::IsFreeForm, x) = surfacepoints(x)
+
+"""Enum to store if a feature should be handled as rough or machined."""
+@enum IsMachined MACHINED ROUGH
+
+
+struct LocalizationTolerance
+    feature1::LocalizationFeature
+    machined1::IsMachined
+    feature2::LocalizationFeature
+    machined2::IsMachined
+    type::LocalizationToleranceType
+    nominalvalue::Float64
+    lowervalue::Float64
+    uppervalue::Float64
+    note::String
+end
+
+function toleranceddistance(type::PlanePlaneDistance, t::LocalizationTolerance)
+    f1 = t.feature1
+    m1 = t.machined1
+    f2 = t.feature2
+    m2 = t.machined2
+
+
+    #check if plane normals are parallel
+    zaxis1 = zaxis(getpartzero(f1))
+    zaxis2 = zaxis(getpartzero(f2))
+    if abs(dot(zaxis1, zaxis2)) < cosd(5)
+        error("Planes: $(getfeaturename(f1)) and $(getfeaturename(f2)) are not parallel!
+            Can't compute `PlanePlaneDistance`")
+    end
+
+    v_f1_ = m1 == MACHINED ? getfeaturepoints(f1.machined) : getfeaturepoints(f1.rough)
+    v_f2_ = m2 == MACHINED ? getfeaturepoints(f2.machined) : getfeaturepoints(f2.rough)
+    v_f1 = m1 == MACHINED ? transformmachined2datum(f1, v_f1_) : v_f1_
+    v_f2 = m2 == MACHINED ? transformmachined2datum(f2, v_f2_) : v_f2_
+
+    # pairwise distance
+    #TODO: not good for two IsFreeForm surfaces!!!!
+    difference_vectors = (v2-v1 for v1 in v_f1 for v2 in v_f2)
+    signed_distances = (dot(zaxis1, dv) for dv in difference_vectors)
+    d = mean(abs.(signed_distances))
+    return d
+end
+
+function toleranceddistance(type::PlaneAxisDistance, t::LocalizationTolerance)
+    f1 = t.feature1
+    m1 = t.machined1
+    f2 = t.feature2
+    m2 = t.machined2
+    # this is rather plane-point distance
+    # 1. get plane normal -> part zero z axis
+    # 2. get feature point of both features
+    # 3. compute the distance of the feature points
+    # 4. project the distance to the z axis
+    
+    # current implementation only handles IsPrimitive features
+    geom1 = m1 == MACHINED ? f1.machined : f1.rough
+    geom2 = m2 == MACHINED ? f2.machined : f2.rough
+    gs1 = GeometryStyle(typeof(geom1))
+    gs2 = GeometryStyle(typeof(geom2))
+
+    if (gs1) != IsPrimitive() || gs2 != IsPrimitive()
+        error("PlaneAxisDistance is only implemented when both features are IsPrimitive!
+        f1 is $gs1, f2 is $gs2")
+    end
+
+    if ! (geom1 isa AbstractPlaneGeometry) || ! (geom2 isa AbstractHoleGeometry)
+        error("PlaneAxisDistance is defined for a plane and a hole.
+        Got types: $(typeof(geom1)) and $(typeof(geom2))")
+    end
+
+    plane_n = zaxis(getpartzero(f1))
+    fp_plane = m1 == MACHINED ? getmachinedfeaturepointindatum(f1) : getroughfeaturepoint(f1)
+    fp_hole = m2 == MACHINED ? getmachinedfeaturepointindatum(f2) : getroughfeaturepoint(f2)
+
+    distancev = fp_plane - fp_hole
+    distance = abs(dot(plane_n, distancev))
+    
+    return distance
+end
+
+function toleranceddistance(type::AxisAxisDistance, t::LocalizationTolerance)
+    f1 = t.feature1
+    m1 = t.machined1
+    f2 = t.feature2
+    m2 = t.machined2
+    # this is rather point-point distance
+    # 1. get feature points
+    # 2. get axes of feature points
+    # 3. compute the distance of the feature points
+    # 4. project the distance along cross product of the two axes
+    
+    # current implementation only handles IsPrimitive features
+    geom1 = m1 == MACHINED ? f1.machined : f1.rough
+    geom2 = m2 == MACHINED ? f2.machined : f2.rough
+    gs1 = GeometryStyle(typeof(geom1))
+    gs2 = GeometryStyle(typeof(geom2))
+
+    if (gs1) != IsPrimitive() || gs2 != IsPrimitive()
+        error("AxisAxisDistance is only implemented when both features are IsPrimitive!
+        f1 is $gs1, f2 is $gs2")
+    end
+
+    if ! (geom1 isa AbstractHoleGeometry) || ! (geom2 isa AbstractHoleGeometry)
+        error("AxisAxisDistance is defined for a hole and a hole.
+        Got types: $(typeof(geom1)) and $(typeof(geom2))")
+    end
+
+    fp1 = m1 == MACHINED ? getmachinedfeaturepointindatum(f1) : getroughfeaturepoint(f1)
+    fp2 = m2 == MACHINED ? getmachinedfeaturepointindatum(f2) : getroughfeaturepoint(f2)
+
+    distancev = fp1 - fp2
+    distance = abs(dot(type.projectionaxis, distancev))
+
+    return distance
+end
+
+function aac_primitive_primitive(f1, m1, f2, m2)
+    #check if axes are parallel
+    zaxis1 = zaxis(getpartzero(f1))
+    zaxis2 = zaxis(getpartzero(f2))
+    if abs(dot(zaxis1, zaxis2)) < cosd(5)
+        error("Holes' axes are not parallel!
+        Can't compute `AxisAxisConcentric` tolerance.")
+    end
+
+    fp1 = m1 == MACHINED ? getmachinedfeaturepointindatum(f1) : getroughfeaturepoint(f1)
+    fp2 = m2 == MACHINED ? getmachinedfeaturepointindatum(f2) : getroughfeaturepoint(f2)
+
+    # distancev is in workpiece datum
+    distancev = fp1 - fp2
+    # needs to be in part zero orientation, so the "x-y distance" of the two feature points
+    # can be calculated
+    iR = inv(getpartzero(f1).rotation)
+    dv = iR*distancev
+    distance = norm(dv[1:2])
+
+    return distance
+end
+
+function aac_primitive_freeform(t::LocalizationTolerance, f1, m1, f2, m2)
+    #check if axes are parallel
+    # 1. get fp and radius of feature 1
+    # 2. get surface points of feature 2
+    # 3. calculate the distances from fp1 and surface points of 2
+    # 4. subtract t-nominalvalue
+    # 5. calculate the mean of those
+    zaxis1 = zaxis(getpartzero(f1))
+    zaxis2 = zaxis(getpartzero(f2))
+    if abs(dot(zaxis1, zaxis2)) < cosd(5)
+        error("Holes' axes are not parallel!
+        Can't compute `AxisAxisConcentric` tolerance.")
+    end
+    fp1 = m1 == MACHINED ? getmachinedfeaturepointindatum(f1) : getroughfeaturepoint(f1)
+    fp2s = m2 == MACHINED ? transformmachined2datum(f2, surfacepoints(f2.machined)) : surfacepoints(f2.rough)
+
+    # needs to be in part zero orientation, so the "x-y distance" of the two feature points
+    # can be calculated
+    iR = inv(getpartzero(f1).rotation)
+    # distancev is transformed into part zero orientation
+    distancevs = (iR*(fp1-v) for v in fp2s)
+    # x-y distance is calculated and subtract the t.nominalvalue
+    #distances = (norm(v[1:2])-t.nominalvalue for v in distancevs)
+    distances = (norm(v[1:2]) for v in distancevs)
+    distance = mean(distances)
+
+    return distance
+end
+
+function toleranceddistance(type::AxisAxisConcentric, t::LocalizationTolerance)
+    f1 = t.feature1
+    m1 = t.machined1
+    f2 = t.feature2
+    m2 = t.machined2
+    # this is rather point-point distance
+    # 1. get feature points
+    # 2. compute the distance of the feature points
+    # 4. project the distance to the normal of the axes' of the holes
+    
+    # current implementation only handles IsPrimitive features
+    geom1 = m1 == MACHINED ? f1.machined : f1.rough
+    geom2 = m2 == MACHINED ? f2.machined : f2.rough
+
+    if ! (geom1 isa AbstractHoleGeometry) || ! (geom2 isa AbstractHoleGeometry)
+        error("AxisAxisConcentric is defined for a hole and a hole.
+        Got types: $(typeof(geom1)) and $(typeof(geom2))")
+    end
+
+    gs1 = GeometryStyle(typeof(geom1))
+    gs2 = GeometryStyle(typeof(geom2))
+
+    if gs1 == IsPrimitive() && gs2 == IsPrimitive()
+        # this computes the distance between the axes of two holes
+        # ideally this should be zero, so optimizing to tolerance center doesnt really work
+        return aac_primitive_primitive(f1, m1, f2, m2)
+    elseif gs1 == IsPrimitive() && gs2 == IsFreeForm()
+        # this computes the distance between the axis of a hole (part zero z axis) and
+        # points of a freeform feature
+        return aac_primitive_freeform(t, f1, m1, f2, m2)
+    else
+        error("AxisAxisConcentric is not only implemented for this pair of geometries!
+        f1 is $gs1, f2 is $gs2")
+    end   
+end
+
+
+function toleranceddistance(t::LocalizationTolerance)
+    toleranceddistance(t.type, t)
+end
+
+function evaluatetolerance(t::LocalizationTolerance)
+    d = toleranceddistance(t)
+    abs_d = d - (t.lowervalue+t.uppervalue)/2
+    rel_d = 2*abs_d/(t.uppervalue-t.lowervalue)*100
+    return (d, rel_d)
+end
+
+=#
diff --git a/src/visualization.jl b/src/visualization.jl
index e821690..b1c56a5 100644
--- a/src/visualization.jl
+++ b/src/visualization.jl
@@ -4,14 +4,12 @@
 Get the rough state of `lf` and transform it according to its current part zero transformation.
 """
 function transformmachinedgeoms(lf::LocalizationFeature)
-    # for properly working, needs this to be merged:
-    # https://github.com/JuliaGeometry/Meshes.jl/pull/623
-    pz = getpartzero(lf)
+    pz = partzero(lf)
     R = pz.rotation
     rot = RotMatrix{3}(R)
     fr = Rotate(rot)
     ft = Translate(pz.position...)
-    geom = visualizationgeometry(lf.machined)
+    geom = visualizationgeometry(geometry(lf.machinedfeature))
     geom2 = fr(geom)
     return ft(geom2)
 end
@@ -22,7 +20,7 @@ end
 Generate Meshes object for each machined hole.
 """
 function genmachinedholes(mop::MultiOperationProblem)
-    holes = collectmachinedholes(mop)
+    holes = cylindricalfeatures(mop)
     return [transformmachinedgeoms(h) for h in holes]
 end
 
@@ -32,8 +30,8 @@ end
 Generate Meshes object for each rough hole.
 """
 function genroughholes(mop::MultiOperationProblem)
-    holes = collectroughholes(mop)
-    return [visualizationgeometry(h.rough) for h in holes]
+    holes = cylindricalfeatures(mop)
+    return [visualizationgeometry(geometry(h.roughfeature)) for h in holes]
 end
 
 """
@@ -42,7 +40,7 @@ end
 Generate Meshes object for each machined plane.
 """
 function genmachinedplanes(mop::MultiOperationProblem)
-    planes = collectmachinedplanes(mop)
+    planes = planarfeatures(mop)
     return [transformmachinedgeoms(p) for p in planes]
 end
 
@@ -52,6 +50,6 @@ end
 Generate Meshes object for each rough plane.
 """
 function genroughplanes(mop::MultiOperationProblem)
-    planes = collectroughplanes(mop)
-    return [visualizationgeometry(p.rough) for p in planes]
+    planes = planarfeatures(mop)
+    return [visualizationgeometry(geometry(p.roughfeature)) for p in planes]
 end
diff --git a/test/geometries.jl b/test/geometries.jl
index e199ee7..e7abd1c 100644
--- a/test/geometries.jl
+++ b/test/geometries.jl
@@ -1,47 +1,45 @@
 @testset "IsPrimitive geometries" begin
-    sh = SimpleHole([0, 0, 0], 29)
-    sp = SimplePlane([0, 0, 0])
-
-    @test featurepoint(sh) == [0, 0, 0]
-    @test featureradius(sh) == 29
-    @test featurepoint(sp) == [0, 0, 0]
-
-    @test_throws MethodError featureradius(sp)
-    @test_throws ErrorException surfacepoints(sp)
-    @test_throws ErrorException surfacepoints(sh)
-    @test_throws ErrorException filteredsurfacepoints(sp)
-    @test_throws ErrorException filteredsurfacepoints(sh)
-
-    pz1 = PartZero("pz1", [0, 0, 0], hcat([0, 1, 0], [0, 0, 1], [1, 0, 0]))
-    sh_r = SimpleHole([82.5, 30, 40], 26)
-    sp_r = PlaneAndNormal([82.5, 30, 40], [1, 0, 0])
-    fd_sh = FeatureDescriptor("simple-hole", pz1, true, true)
-    fd_sp = FeatureDescriptor("simple-plane", pz1, true, true)
-
-    h1 = HoleLocalizationFeature(fd_sh, sh_r, sh)
-    p1 = PlaneLocalizationFeature(fd_sp, sp_r, sp)
-
-    @test BLC.getfeaturename(h1) == "simple-hole"
-    @test BLC.getfeaturename(p1) == "simple-plane"
-    @test BLC.getpartzero(h1) === BLC.getpartzero(p1) === pz1
-    @test BLC.getpartzeroname(h1) == "pz1"
-    @test BLC.hasrough(h1)
-    @test BLC.hasrough(p1)
-    @test BLC.hasmachined(h1)
-    @test BLC.hasmachined(p1)
-
-    @test BLC.getroughfeaturepoint(h1) == [82.5, 30, 40]
-    @test BLC.getroughfeaturepoint(p1) == [82.5, 30, 40]
-    @test BLC.getmachinedfeaturepoint(h1) == [0, 0, 0]
-    @test BLC.getmachinedfeaturepoint(p1) == [0, 0, 0]
-    @test BLC.getmachinedradius(h1) == 29
-    @test_throws MethodError BLC.getmachinedradius(p1)
-    @test BLC.getroughradius(h1) == 26
-    @test_throws MethodError BLC.getroughradius(p1)
-
-    @test_throws ErrorException BLC.getroughfilteredpoints(h1)
-    @test_throws ErrorException BLC.getroughfilteredpoints(p1)
-
-    @test BLC.getmachinedfeaturepointindatum(h1) == [0, 0, 0]
-    @test BLC.getmachinedfeaturepointindatum(p1) == [0, 0, 0]
+    tpoints = [(0,0),(1,0),(1,1)]
+    tconnec = [connect((1,2,3))]
+    tmesh = SimpleMesh(tpoints, tconnec)
+    cyl1 = Cylinder(15.0)
+
+    pz1 = PartZero("pz1", [0,0,0], [1 0 0;0 1 0;0 0 1])
+
+    # IsFreeForm, HOLELIKE
+    r1 = RoughFeature("rh1", HOLELIKE, tmesh)
+    @test GeometryStyle(r1) == IsFreeForm()
+    @test isholelike(r1)
+    @test_throws ErrorException featureradius(r1)
+    @test_throws ErrorException featurepoint(r1)
+
+    # IsPrimitive, HOLELIKE
+    m1 = MachinedFeature("mh1", HOLELIKE, cyl1, pz1)
+    @test GeometryStyle(m1) == IsPrimitive()
+    @test isholelike(m1)
+    @test featurepoint(m1) == Point3(0,0,1)
+    @test featureradius(m1) == 15
+    @test_throws ErrorException surfacepoints(m1)
+
+    f1 = LocalizationFeature("h1", r1, m1)
+    @test GeometryStyle(f1) == IsFreeForm()
+
+    # IsPrimitive, PLANELIKE
+    r2 = RoughFeature("rp1", PLANELIKE, Plane(Point3(0,0,0), Vec3(0,0,1)))
+    @test GeometryStyle(r2) == IsPrimitive()
+    @test isplanelike(r2)
+    @test featurepoint(r2) == Point3(0,0,0)
+    @test_throws ErrorException featureradius(r2)
+    @test_throws ErrorException surfacepoints(r2)
+
+    # IsPrimitive, PLANELIKE
+    m2 = MachinedFeature("mp1", PLANELIKE, Plane(Point3(0,0,0), Vec3(0,0,1)), pz1)
+    @test GeometryStyle(m2) == IsPrimitive()
+    @test isplanelike(m2)
+    @test featurepoint(m2) == Point3(0,0,0)
+    @test_throws ErrorException featureradius(m2)
+    @test_throws ErrorException surfacepoints(m2)
+
+    f2 = LocalizationFeature("p1", r2, m2)
+    @test GeometryStyle(f2) == IsPrimitive()
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index b34fd96..806733e 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,11 +1,12 @@
 using BlankLocalizationCore
 using Test
-
 using Meshes
-using Ipopt
 
 const BLC = BlankLocalizationCore
 
+#using Ipopt
+
+#include("tolerances.jl")
 include("partzeros.jl")
 include("geometries.jl")
-include("testproblem.jl")
+#include("testproblem.jl")
diff --git a/test/testproblem.jl b/test/testproblem.jl
index f340601..a8774a5 100644
--- a/test/testproblem.jl
+++ b/test/testproblem.jl
@@ -12,36 +12,36 @@
     ## Machined geometry definitions
 
     fronthole_m = SimpleHole([0, 0, 0], 29)
-    frontface_m = SimplePlane([0, 0, 0])
+    frontface_m = SimplePlane([0, 0, 0], [1, 0, 0])
 
     righthole1_m = SimpleHole([16, 15, 0], 7.5)
     righthole2_m = SimpleHole([25, -16, 3], 9)
     righthole3_m = SimpleHole([60, 0, -3], 13.5)
-    rightface1_m = SimplePlane([16, 15, 0])
-    rightface2_m = SimplePlane([25, -16, 3])
-    rightface3_m = SimplePlane([60, 0, -3])
+    rightface1_m = SimplePlane([16, 15, 0], [0, 1, 0])
+    rightface2_m = SimplePlane([25, -16, 3], [0, 1, 0])
+    rightface3_m = SimplePlane([60, 0, -3], [0, 1, 0])
 
     backhole1_m = SimpleHole([-14, 14, 0], 9)
     backhole2_m = SimpleHole([14, 14, 0], 9)
-    backface1_m = SimplePlane([-14, 14, 0])
-    backface2_m = SimplePlane([14, 14, 0])
+    backface1_m = SimplePlane([-14, 14, 0], [-1, 0, 0])
+    backface2_m = SimplePlane([14, 14, 0], [-1, 0, 0])
 
     ## Rough geometry definitions
 
     fronthole_r = SimpleHole([82.5, 30, 40], 26)
-    frontface_r = PlaneAndNormal([82.5, 30, 40], [1, 0, 0])
+    frontface_r = SimplePlane([82.5, 30, 40], [1, 0, 0])
 
     righthole1_r = SimpleHole([66, 71.5, 55], 6)
     righthole2_r = SimpleHole([58, 74.5, 24], 4.905)
     righthole3_r = SimpleHole([21.5, 68.5, 40], 8)
-    rightface1_r = PlaneAndNormal([66, 71.5, 55], [0, 1, 0])
-    rightface2_r = PlaneAndNormal([58, 74.5, 24], [0, 1, 0])
-    rightface3_r = PlaneAndNormal([21.5, 68.5, 40], [0, 1, 0])
+    rightface1_r = SimplePlane([66, 71.5, 55], [0, 1, 0])
+    rightface2_r = SimplePlane([58, 74.5, 24], [0, 1, 0])
+    rightface3_r = SimplePlane([21.5, 68.5, 40], [0, 1, 0])
 
     backhole1_r = SimpleHole([-3, 44, 53.9], 6.2)
     backhole2_r = SimpleHole([-3, 16.1, 54], 6.25)
-    backface1_r = PlaneAndNormal([-3, 44, 54], [-1, 0, 0])
-    backface2_r = PlaneAndNormal([-3, 16, 54], [-1, 0, 0])
+    backface1_r = SimplePlane([-3, 44, 54], [-1, 0, 0])
+    backface2_r = SimplePlane([-3, 16, 54], [-1, 0, 0])
 
     ## Geometry pairing and feature descriptors
 
diff --git a/test/tolerances.jl b/test/tolerances.jl
new file mode 100644
index 0000000..0f63103
--- /dev/null
+++ b/test/tolerances.jl
@@ -0,0 +1,210 @@
+@testset "PlanePlaneDistance: tolerancedistance" begin
+    empty_plane = SimplePlane([0,0,0], [-1, 0, 0])
+
+    machined_plane1 = SimplePlane([0,0,0], [-1, 0, 0])
+    machined_plane2 = SimplePlane([5,6,7], [-1, 0, 0])
+
+    mpoints = Point3[(13,9,7.1), (-100,70,6.9), (5000,-2,7.)]
+    sm1 = SimpleMesh(mpoints, [connect((1,2,3))])
+    mesh_plane1 = MeshPlane(sm1)
+
+    pz1 = PartZero("testpz1", [0,0,0], [1 0 0;0 1 0;0 0 1])
+    pz2 = PartZero("testpz2", [0,0,0], [0 0 1;0 1 0;0 1 0])
+
+    # lf1 has machined_plane1 as machined feature
+    lf1 = BLC.LocalizationFeature("t1", pz1, empty_plane, machined_plane1)
+    # lf2 has machined_plane2 as machined feature
+    lf2 = BLC.LocalizationFeature("tf1", pz1, empty_plane, machined_plane2)
+    # lf3 has "wrong" z axis
+    lf3 = BLC.LocalizationFeature("t3", pz2, empty_plane, empty_plane)
+    # lf4 has both features a mesh, and "correct" part zero
+    lf4 = BLC.LocalizationFeature("t4", pz1, mesh_plane1, mesh_plane1)
+    
+    # testing machined-machined distance
+    # nominal, lower, upper, and note doesn't matter here
+    t1 = LocalizationTolerance(lf1, BLC.MACHINED, lf2, BLC.MACHINED, PlanePlaneDistance(), 0, 0, 0, "")
+    @test isapprox(toleranceddistance(t1), 7)
+
+    # testing if not parallel part zeros throw error
+    t2 = LocalizationTolerance(lf1, BLC.MACHINED, lf3, BLC.MACHINED, PlanePlaneDistance(), 0, 0, 0, "")
+    @test_throws ErrorException toleranceddistance(t2)
+
+    # testing IsFreeForm vs IsPrimitive
+    # both MACHINED and ROUGH
+    t3 = LocalizationTolerance(lf4, BLC.MACHINED, lf1, BLC.MACHINED, PlanePlaneDistance(), 0, 0, 0, "")
+    t4 = LocalizationTolerance(lf4, BLC.ROUGH, lf1, BLC.MACHINED, PlanePlaneDistance(), 0, 0, 0, "")
+    @test isapprox(toleranceddistance(t3), 7)
+    @test isapprox(toleranceddistance(t4), 7)
+
+    ## test if order of features doesn't matter
+    t1_r = LocalizationTolerance(lf2, BLC.MACHINED, lf1, BLC.MACHINED, PlanePlaneDistance(), 0, 0, 0, "")
+    @test isapprox(toleranceddistance(t1), toleranceddistance(t1_r))
+
+    t2_r = LocalizationTolerance(lf3, BLC.MACHINED, lf1, BLC.MACHINED, PlanePlaneDistance(), 0, 0, 0, "")
+    @test_throws ErrorException toleranceddistance(t2_r)
+
+    t3_r = LocalizationTolerance(lf1, BLC.MACHINED, lf4, BLC.MACHINED, PlanePlaneDistance(), 0, 0, 0, "")
+    t4_r = LocalizationTolerance(lf1, BLC.MACHINED, lf4, BLC.ROUGH, PlanePlaneDistance(), 0, 0, 0, "")
+    @test isapprox(toleranceddistance(t3), toleranceddistance(t3_r))
+    @test isapprox(toleranceddistance(t4), toleranceddistance(t4_r))
+end
+
+@testset "PlaneAxisDistance: tolerancedistance" begin
+    # test error throwing -> should be addressed by #7
+    machined_plane1 = SimplePlane([0,0,0], [-1, 0, 0])
+
+    mpoints = Point3[(13,9,7.1), (-100,70,6.9), (5000,-2,7.)]
+    sm1 = SimpleMesh(mpoints, [connect((1,2,3))])
+    mesh_plane1 = MeshPlane(sm1)
+
+    pz1_front = PartZero("testpz1", [0,0,0], [1 0 0;0 1 0;0 0 1])
+    pz2_right = PartZero("testpz2", [10,10,10], [-1 0 0;0 0 1;0 1 0])
+
+    # lf1 simple plane and simple plane
+    lf1 = BLC.LocalizationFeature("t1", pz1_front, machined_plane1, machined_plane1)
+    # lf2 simple plane and mash plane
+    lf2 = BLC.LocalizationFeature("tf1", pz1_front, machined_plane1, mesh_plane1)
+
+    t1_mr = LocalizationTolerance(lf1, BLC.ROUGH, lf2, BLC.MACHINED, PlaneAxisDistance(), 0, 0, 0, "")
+    t1_rm = LocalizationTolerance(lf2, BLC.MACHINED, lf1, BLC.ROUGH, PlaneAxisDistance(), 0, 0, 0, "")
+    @test_throws ErrorException toleranceddistance(t1_mr)
+    @test_throws ErrorException toleranceddistance(t1_rm)
+
+    # plane and hole test
+    t1_mm = LocalizationTolerance(lf1, BLC.MACHINED, lf2, BLC.ROUGH, PlaneAxisDistance(), 0, 0, 0, "")
+    @test_throws ErrorException toleranceddistance(t1_mm)
+
+    ## test distance calculation
+    sh1 = SimpleHole([10, 0, 0], 15)
+    sh2 = SimpleHole([0, -10, 0], 15)
+    empty_hole = SimpleHole([0,0,0], 0)
+
+    pl1 = SimplePlane([15, 0, 0], [-1, 0, 0])
+    pl2 = SimplePlane([0, 0, 15], [-1, 0, 0])
+    empty_plane = SimplePlane([0,0,0], [-1, 0, 0])
+
+    h1 = BLC.LocalizationFeature("fronth1", pz1_front, empty_hole, sh1)
+    h2 = BLC.LocalizationFeature("fronth1", pz1_front, empty_hole, sh2)
+
+    p1 = BLC.LocalizationFeature("rightp1", pz2_right, empty_plane, pl1)
+    p2 = BLC.LocalizationFeature("rightp2", pz2_right, empty_plane, pl2)
+
+    lt11 = LocalizationTolerance(p1, BLC.MACHINED, h1, BLC.MACHINED, PlaneAxisDistance(), 10, 10, 10, "")
+    lt12 = LocalizationTolerance(p1, BLC.MACHINED, h2, BLC.MACHINED, PlaneAxisDistance(), 20, 20, 20, "")
+
+    lt21 = LocalizationTolerance(p2, BLC.MACHINED, h1, BLC.MACHINED, PlaneAxisDistance(), 25, 25, 25, "")
+    lt22 = LocalizationTolerance(p2, BLC.MACHINED, h2, BLC.MACHINED, PlaneAxisDistance(), 35, 35, 35, "")
+
+    @test isapprox(toleranceddistance(lt11), 10)
+    @test isapprox(toleranceddistance(lt12), 20)
+    @test isapprox(toleranceddistance(lt21), 25)
+    @test isapprox(toleranceddistance(lt22), 35)
+end
+
+@testset "AxisAxisDistance: tolerancedistance" begin
+    ## test error throwing -> should be addressed by #7
+    machined_plane1 = SimplePlane([0,0,0], [-1, 0, 0])
+
+    mpoints = Point3[(13,9,7.1), (-100,70,6.9), (5000,-2,7.)]
+    sm1 = SimpleMesh(mpoints, [connect((1,2,3))])
+    mesh_plane1 = MeshPlane(sm1)
+
+    # part zeros are different of the previous ones
+    pz1_front = PartZero("frontpz1", [0,0,0], [0 0 1;1 0 0;0 1 0])
+    pz2_right = PartZero("rightpz2", [10,10,10], [-1 0 0;0 0 1;0 1 0])
+
+    # lf1 simple plane and simple plane
+    lf1 = BLC.LocalizationFeature("t1", pz1_front, machined_plane1, machined_plane1)
+    # lf2 simple plane and mesh plane
+    lf2 = BLC.LocalizationFeature("tf1", pz1_front, machined_plane1, mesh_plane1)
+
+    t1_mr = LocalizationTolerance(lf1, BLC.ROUGH, lf2, BLC.MACHINED, AxisAxisDistance([1,0,0]), 0, 0, 0, "")
+    t1_rm = LocalizationTolerance(lf2, BLC.MACHINED, lf1, BLC.ROUGH, AxisAxisDistance([1,0,0]), 0, 0, 0, "")
+    @test_throws ErrorException toleranceddistance(t1_mr)
+    @test_throws ErrorException toleranceddistance(t1_rm)
+
+    # plane and hole test
+    t1_mm = LocalizationTolerance(lf1, BLC.MACHINED, lf2, BLC.ROUGH, AxisAxisDistance([1,0,0]), 0, 0, 0, "")
+    @test_throws ErrorException toleranceddistance(t1_mm)
+
+    ## test distance calculation
+
+    sh1_machined = SimpleHole([10, 0, 0], 15)
+    sh2_machined = SimpleHole([5, 10, 15], 15)
+    sh1_rough = SimpleHole([0, 10, 0], 14.9)
+    sh2_rough = SimpleHole([5, 25, 20], 14.9)
+
+    fh1 = BLC.LocalizationFeature("fronth1", pz1_front, sh1_rough, sh1_machined)
+
+    rh1 = BLC.LocalizationFeature("righth1", pz2_right, sh2_rough, sh2_machined)
+
+    # all comibinations of rough-machined should be the same
+    aat1 = LocalizationTolerance(fh1, BLC.MACHINED, rh1, BLC.MACHINED, AxisAxisDistance([0, 0, 1]), 20, 20, 20, "")
+    aat2 = LocalizationTolerance(fh1, BLC.MACHINED, rh1, BLC.ROUGH, AxisAxisDistance([0, 0, 1]), 20, 20, 20, "")
+    aat3 = LocalizationTolerance(fh1, BLC.ROUGH, rh1, BLC.MACHINED, AxisAxisDistance([0, 0, 1]), 20, 20, 20, "")
+    aat4 = LocalizationTolerance(fh1, BLC.ROUGH, rh1, BLC.ROUGH, AxisAxisDistance([0, 0, 1]), 20, 20, 20, "")
+    @test isapprox(toleranceddistance(aat1), 20)
+    @test isapprox(toleranceddistance(aat2), 20)
+    @test isapprox(toleranceddistance(aat3), 20)
+    @test isapprox(toleranceddistance(aat4), 20)
+end
+
+@testset "AxisAxisConcentric: primitive-primitive tolerancedistance" begin
+    pz2_right = PartZero("rightpz2", [10,10,10], [-1 0 0;0 0 1;0 1 0])
+    pz3_right = PartZero("rightpz3", [0,0,0], [-1 0 0;0 0 1;0 1 0])
+    
+    aac1_m = SimpleHole([0,0,19], 15)
+    aac1_r = SimpleHole([10.1,28.9,10.1], 15)
+    
+    aac2_m = SimpleHole([-10, 10, 37], 15)
+    aac2_r = SimpleHole([9, 36, 9], 15)
+    
+    # distance aac1_m -> aac1_r: sqrt(2)/10
+    # distance aac1_m -> aac2_m: 0
+    
+    # distance aac2_m -> aac2_r: sqrt(2)
+    
+    aac_lf1 = BLC.LocalizationFeature("righth1", pz2_right, aac1_r, aac1_m)
+    aac_lf2 = BLC.LocalizationFeature("righth2", pz3_right, aac2_r, aac2_m)
+    
+    aac1 = LocalizationTolerance(aac_lf1, BLC.MACHINED, aac_lf1, BLC.ROUGH, AxisAxisConcentric(), sqrt(2)/10, 0,0, "")
+    aac2 = LocalizationTolerance(aac_lf1, BLC.MACHINED, aac_lf2, BLC.MACHINED, AxisAxisConcentric(), 0, 0,0, "")
+    aac3 = LocalizationTolerance(aac_lf2, BLC.ROUGH, aac_lf2, BLC.MACHINED, AxisAxisConcentric(), sqrt(2), 0,0, "")
+    
+    @test isapprox(toleranceddistance(aac1), sqrt(2)/10)
+    @test isapprox(toleranceddistance(aac2), 0)
+    @test isapprox(toleranceddistance(aac3), sqrt(2))
+end
+
+@testset "AxisAxisConcentric: primitive-freeform tolerancedistance" begin
+    pz_aac = PartZero("aac", [0,0,0], [1 0 0;0 1 0;0 0 1])
+
+    # unit circle, with center at (0,0,0)
+    aacpf_h = SimpleHole([0,0,0], 1)
+
+    aac_ps1 = Point3[(0,1,0),(1,0,1),(cosd(45), sind(45),-1)]
+    aac_c1 = [connect((1,2,3))]
+    aac_mesh1 = MeshHole(SimpleMesh(aac_ps1, aac_c1), [])
+
+    aacpf_lf1 = BLC.LocalizationFeature("h1", pz_aac, aac_mesh1, aacpf_h)
+    t_aacpf1 = LocalizationTolerance(aacpf_lf1, BLC.MACHINED, aacpf_lf1, BLC.ROUGH, AxisAxisConcentric(), 1,1,1,"")
+    @test isapprox(toleranceddistance(t_aacpf1), 1)
+
+    # every point is on radius 2 circle
+    aac_ps2 = Point3[(0,2,14), (-2,0,0), (2*cosd(30), 2*sind(30),-1), (2*cosd(130), 2*sind(130), 17)]
+    sm2 = SimpleMesh(aac_ps2, aac_c1)
+    aac_mesh2 = MeshHole(sm2, [])
+
+    aacpf_lf2 = BLC.LocalizationFeature("h1", pz_aac, aac_mesh2, aacpf_h)
+    t_aacpf2 = LocalizationTolerance(aacpf_lf2, BLC.MACHINED, aacpf_lf2, BLC.ROUGH, AxisAxisConcentric(), 1,0,2,"")
+    @test isapprox(toleranceddistance(t_aacpf2), 2)
+
+    # every point is on 0.5 radius circle
+    aac_ps3 = Point3[(0,0.5,0), (-0.5,0,4), (0.5*cosd(30), 0.5*sind(30),-6), (0.5*cosd(130), 0.5*sind(130), 17)]
+    sm3 = SimpleMesh(aac_ps3, aac_c1)
+    aac_mesh3 = MeshHole(sm3, [])
+
+    aacpf_lf3 = BLC.LocalizationFeature("h1", pz_aac, aac_mesh3, aacpf_h)
+    t_aacpf3 = LocalizationTolerance(aacpf_lf3, BLC.MACHINED, aacpf_lf3, BLC.ROUGH, AxisAxisConcentric(), 1,0.5,1.5,"")
+    @test isapprox(toleranceddistance(t_aacpf3), 0.5)
+end
\ No newline at end of file