From 7141b93a7b595b1b15f5725660f153641ce26f38 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tam=C3=A1s=20Cserteg?= <csertegt@gmail.com>
Date: Thu, 6 Jul 2023 16:05:49 +0200
Subject: [PATCH 01/27] clean up PartZero functions a little bit

---
 src/partzeros.jl | 20 ++------------------
 1 file changed, 2 insertions(+), 18 deletions(-)

diff --git a/src/partzeros.jl b/src/partzeros.jl
index 2fe28db..d0c2d36 100644
--- a/src/partzeros.jl
+++ b/src/partzeros.jl
@@ -8,7 +8,8 @@ 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]
@@ -55,7 +56,6 @@ function inverthomtr(M)
     return invtr
 end
 
-
 """
     getpartzerobyname(partzeros::Vector{PartZero}, partzeroname::AbstractString)
 
@@ -68,22 +68,6 @@ function getpartzerobyname(partzeros::Vector{PartZero}, partzeroname::AbstractSt
     return nothing
 end
 
-"""
-    getcoordinateindatumorientation(partzero::PartZero, v, axis::Int)
-
-Get the `axis`th coordinate of vector `v` in the orientation of the workpiece datum.
-(Inverse transformed with the part zero, the part zero's position being zero.)
-"""
-function getcoordinateindatumorientation(partzero::PartZero, v, axis::Int)
-    M = getpartzeroHM(partzero)
-    M[1:3,4] = [0,0,0]
-    #iM = inv(M)
-    fv = M*vcat(v, 1)
-    # not sure why, but part zero is needed, not its inverted part
-    #fv = iM*vcat(v, 1)
-    return fv[axis]
-end
-
 """
     printpartzeropositions(partzeros::Vector{PartZero})
 

From 5fa9b71f34bb3771234ece15c71c4f11aeeba991 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tam=C3=A1s=20Cserteg?= <csertegt@gmail.com>
Date: Thu, 6 Jul 2023 16:10:07 +0200
Subject: [PATCH 02/27] split geometries and optimization problem definition

---
 src/BlankLocalizationCore.jl |   6 +-
 src/geometries.jl            | 146 -----------------------------------
 src/optimizationproblem.jl   | 145 ++++++++++++++++++++++++++++++++++
 3 files changed, 149 insertions(+), 148 deletions(-)
 create mode 100644 src/optimizationproblem.jl

diff --git a/src/BlankLocalizationCore.jl b/src/BlankLocalizationCore.jl
index 17cea00..4eac9b0 100644
--- a/src/BlankLocalizationCore.jl
+++ b/src/BlankLocalizationCore.jl
@@ -17,8 +17,9 @@ export  SimpleHole,
         MeshPlane,
         FeatureDescriptor,
         HoleLocalizationFeature,
-        PlaneLocalizationFeature,
-        OptimizationResult,
+        PlaneLocalizationFeature
+
+export  OptimizationResult,
         Tolerance,
         MultiOperationProblem
 
@@ -42,6 +43,7 @@ HV(v) = vcat(v, 1)
 
 include("partzeros.jl")
 include("geometries.jl")
+include("optimizationproblem.jl")
 include("optimization.jl")
 include("resultevaluation.jl")
 
diff --git a/src/geometries.jl b/src/geometries.jl
index a1dfdfb..e265b3b 100644
--- a/src/geometries.jl
+++ b/src/geometries.jl
@@ -196,149 +196,3 @@ function getmachinedfeaturepointindatum(f::LocalizationFeature)
     v_indatum = T*HV(v)
     return v_indatum[1:3]
 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)
-
-struct Tolerance
-    featurename1::String
-    ismachined1::Bool
-    projection::Function
-    featurename2::String
-    ismachined2::Bool
-    nominalvalue::Float64
-    lowervalue::Float64
-    uppervalue::Float64
-    note::String
-end
-
-mutable struct MultiOperationProblem
-    partzeros::Vector{PartZero}
-    holes::Vector{HoleLocalizationFeature}
-    planes::Vector{PlaneLocalizationFeature}
-    tolerances::Vector{Tolerance}
-    parameters::Dict{String,Real}
-    opresult::OptimizationResult
-end
-
-function MultiOperationProblem(partzeros, holes, planes, tolerances, parameters)
-    return MultiOperationProblem(partzeros, holes, planes, tolerances, parameters, emptyor())
-end
-
-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
-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)
-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
-
-"""
-    getfeaturebyname(mop::MultiOperationProblem, featurename)
-
-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`.
-"""
-function getfeaturebyname(mop::MultiOperationProblem, featurename)
-    function retbyname(array, name)
-        for f in array
-            if getfeaturename(f) == name
-                return f
-            end
-        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)
-
-Collect holes that are grouped to part zero called `partzeroname`.
-"""
-function getholesbypartzero(mop::MultiOperationProblem, partzeroname)
-    return filter(x->getpartzeroname(x) == partzeroname, mop.holes)
-end
-
-"""
-    collectmachinedholes(mop::MultiOperationProblem)
-
-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)
diff --git a/src/optimizationproblem.jl b/src/optimizationproblem.jl
new file mode 100644
index 0000000..a4d8fe6
--- /dev/null
+++ b/src/optimizationproblem.jl
@@ -0,0 +1,145 @@
+"""
+    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)
+
+struct Tolerance
+    featurename1::String
+    ismachined1::Bool
+    projection::Function
+    featurename2::String
+    ismachined2::Bool
+    nominalvalue::Float64
+    lowervalue::Float64
+    uppervalue::Float64
+    note::String
+end
+
+mutable struct MultiOperationProblem
+    partzeros::Vector{PartZero}
+    holes::Vector{HoleLocalizationFeature}
+    planes::Vector{PlaneLocalizationFeature}
+    tolerances::Vector{Tolerance}
+    parameters::Dict{String,Real}
+    opresult::OptimizationResult
+end
+
+function MultiOperationProblem(partzeros, holes, planes, tolerances, parameters)
+    return MultiOperationProblem(partzeros, holes, planes, tolerances, parameters, emptyor())
+end
+
+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
+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)
+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
+
+"""
+    getfeaturebyname(mop::MultiOperationProblem, featurename)
+
+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`.
+"""
+function getfeaturebyname(mop::MultiOperationProblem, featurename)
+    function retbyname(array, name)
+        for f in array
+            if getfeaturename(f) == name
+                return f
+            end
+        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)
+
+Collect holes that are grouped to part zero called `partzeroname`.
+"""
+function getholesbypartzero(mop::MultiOperationProblem, partzeroname)
+    return filter(x->getpartzeroname(x) == partzeroname, mop.holes)
+end
+
+"""
+    collectmachinedholes(mop::MultiOperationProblem)
+
+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)

From 6ef1b62cec704d6f865da2c8c9883a325e9eea7d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tam=C3=A1s=20Cserteg?= <csertegt@gmail.com>
Date: Fri, 7 Jul 2023 11:44:05 +0200
Subject: [PATCH 03/27] rework geometry storage

---
 src/BlankLocalizationCore.jl | 12 ++---
 src/geometries.jl            | 90 +++++++++++-------------------------
 2 files changed, 33 insertions(+), 69 deletions(-)

diff --git a/src/BlankLocalizationCore.jl b/src/BlankLocalizationCore.jl
index 4eac9b0..4f0719a 100644
--- a/src/BlankLocalizationCore.jl
+++ b/src/BlankLocalizationCore.jl
@@ -15,10 +15,10 @@ export  SimpleHole,
         SimplePlane,
         MeshHole,
         MeshPlane,
-        FeatureDescriptor,
+        LocalizationFeature,
         HoleLocalizationFeature,
         PlaneLocalizationFeature
-
+#=
 export  OptimizationResult,
         Tolerance,
         MultiOperationProblem
@@ -31,7 +31,7 @@ export  allowancetable,
         printallowancetable,
         tolerancetable,
         printtolerancetable
-
+=#
 """Union type for `Float64` and `Nothing`."""
 const FON = Union{Nothing,Float64}
 
@@ -43,8 +43,8 @@ HV(v) = vcat(v, 1)
 
 include("partzeros.jl")
 include("geometries.jl")
-include("optimizationproblem.jl")
-include("optimization.jl")
-include("resultevaluation.jl")
+#include("optimizationproblem.jl")
+#include("optimization.jl")
+#include("resultevaluation.jl")
 
 end
diff --git a/src/geometries.jl b/src/geometries.jl
index e265b3b..4c6c5c2 100644
--- a/src/geometries.jl
+++ b/src/geometries.jl
@@ -1,4 +1,4 @@
-"""Supertype for localization geometries."""
+"""Supertype of localization geometries."""
 abstract type AbstractLocalizationGeometry end
 
 """Supertype of hole like localization geometries."""
@@ -7,18 +7,16 @@ 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???
-
-
-"""Trait that describes the "style" of an [`AbstractLocalizationGeometry`](@ref)."""
+"""
+Trait that describes the "style" of an [`AbstractLocalizationGeometry`](@ref).
+If it can be described by a primitive, then it is [`IsPrimitive`](@ref) and
+[`IsFreeForm`](@ref) otherwise.
+"""
 abstract type GeometryStyle 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
 
@@ -37,7 +35,6 @@ end
 # free form geometries have surfacepoints
 # return all points of a free form surface
 surfacepoints(x::T) where {T} = surfacepoints(GeometryStyle(T), x)
-#surfacepoints(::IsFreeForm, x) = x.p
 function surfacepoints(::IsPrimitive, x)
     error("Function `surfacepoints` is not defined for `IsFreeForm`` features")
 end
@@ -115,71 +112,39 @@ end
 GeometryStyle(::Type{MeshPlane}) = IsFreeForm()
 
 """
-Store description of a feature: its name, the corresponding part zero, if it has or has not
-a machined and a rough state.
+    LocalizationFeature{R,M}
+
+A feature that is machined and allowance can be computed for it.
+It has a name, a [`PartZero`](@ref), and a rough and machined geometry
+([`AbstractLocalizationGeometry`](@ref)).
+The two geometries should be of same type (hole, plane, etc.),
+but [`HoleLocalizationFeature`](@ref) and [`PlaneLocalizationFeature`](@ref) enforce this
+property
 """
-struct FeatureDescriptor
+struct LocalizationFeature{R,M}
     name::String
     partzero::PartZero
-    hasmachined::Bool
-    hasrough::Bool
+    rough::R
+    machined::M
 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
+const HoleLocalizationFeature{R,M} = LocalizationFeature{R,M} where {R<:AbstractHoleGeometry, M<:AbstractHoleGeometry}
 
-function Base.show(io::IO, fd::FeatureDescriptor)
-    print(io, fd.name, " in ", fd.partzero.name,
-    fd.hasmachined ? "; machined, " : "; ! machined, ",
-    fd.hasrough ? "rough" : "! rough")
-end
+HoleLocalizationFeature(n, p, r, m) = LocalizationFeature(n, p, r, m)
 
-"""
-    LocalizationFeature{R,M}
+const PlaneLocalizationFeature{R,M} = LocalizationFeature{R,M} where {R<:AbstractPlaneGeometry, M<:AbstractPlaneGeometry}
 
-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).
-A `LocalizationFeature` must define if it is [`IsPrimitive`](@ref) or [`IsFreeForm`](@ref)
-based on its rough geometry.
-"""
-abstract type LocalizationFeature{R,M} end
+PlaneLocalizationFeature(n, p, r, m) = LocalizationFeature(n, p, r, m)
 
-# I thought this should cover all subtypes, but it doesn't. But I don't know why
-#GeometryStyle(::Type{LocalizationFeature{R,M}}) where {R,M} = GeometryStyle(R)
+GeometryStyle(::Type{LocalizationFeature{R,M}}) where {R,M} = GeometryStyle(R)
 
 function Base.show(io::IO, lf::LocalizationFeature)
-    print(io, typeof(lf), ": ", lf.descriptor.name)
-end
-
-struct HoleLocalizationFeature{R<:AbstractHoleGeometry,M<:AbstractHoleGeometry} <: LocalizationFeature{R,M}
-    descriptor::FeatureDescriptor
-    rough::R
-    machined::M
+    print(io, typeof(lf), ": ", getfeaturename(lf))
 end
 
-#GeometryStyle(::Type{HoleLocalizationFeature{R,M}}) where {R,M} = GeometryStyle(R)
-
-struct PlaneLocalizationFeature{R<:AbstractPlaneGeometry,M<:AbstractPlaneGeometry} <: LocalizationFeature{R,M}
-    descriptor::FeatureDescriptor
-    rough::R
-    machined::M
-end
-
-#GeometryStyle(::Type{PlaneLocalizationFeature{R,M}}) where {R,M} = GeometryStyle(R)
-
-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)
+getfeaturename(f::LocalizationFeature) = f.name
+getpartzero(f::LocalizationFeature) = getpartzero(f.partzero)
+getpartzeroname(f::LocalizationFeature) = getpartzeroname(f.partzero)
 
 getroughfeaturepoint(f::LocalizationFeature) = featurepoint(f.rough)
 getmachinedfeaturepoint(f::LocalizationFeature) = featurepoint(f.machined)
@@ -189,7 +154,6 @@ 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)

From 09d823ff6b7ce828fed1474ec87331c5e8b1803e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tam=C3=A1s=20Cserteg?= <csertegt@gmail.com>
Date: Fri, 7 Jul 2023 15:06:44 +0200
Subject: [PATCH 04/27] fix error message

---
 src/geometries.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/geometries.jl b/src/geometries.jl
index 4c6c5c2..83e11b5 100644
--- a/src/geometries.jl
+++ b/src/geometries.jl
@@ -36,7 +36,7 @@ end
 # return all points of a free form surface
 surfacepoints(x::T) where {T} = surfacepoints(GeometryStyle(T), x)
 function surfacepoints(::IsPrimitive, x)
-    error("Function `surfacepoints` is not defined for `IsFreeForm`` features")
+    error("Function `surfacepoints` is not defined for `IsPrimitive`` features")
 end
 
 # free form geometries have surfacepoints

From 5d5a6fd88f299766d47d330c0509a44d47c153d0 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tam=C3=A1s=20Cserteg?= <csertegt@gmail.com>
Date: Fri, 7 Jul 2023 15:07:03 +0200
Subject: [PATCH 05/27] get started with tolerance implementation

---
 Project.toml                 |  3 +-
 src/BlankLocalizationCore.jl |  4 +-
 src/tolerances.jl            | 90 ++++++++++++++++++++++++++++++++++++
 3 files changed, 95 insertions(+), 2 deletions(-)
 create mode 100644 src/tolerances.jl

diff --git a/Project.toml b/Project.toml
index 74cecd8..ce8e89f 100644
--- a/Project.toml
+++ b/Project.toml
@@ -11,10 +11,11 @@ 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"
 
 [compat]
-JuMP = "1.4"
 DataFrames = "1"
+JuMP = "1.4"
 PrettyTables = "2"
 julia = "1.6"
 
diff --git a/src/BlankLocalizationCore.jl b/src/BlankLocalizationCore.jl
index 4f0719a..1d66304 100644
--- a/src/BlankLocalizationCore.jl
+++ b/src/BlankLocalizationCore.jl
@@ -5,7 +5,8 @@ using DataFrames: DataFrame, names, nrow
 using PrettyTables: pretty_table, ft_nonothing, tf_html_minimalist
 using Meshes: SimpleMesh, vertices, boundingbox
 using Logging: @warn
-using LinearAlgebra: norm
+using LinearAlgebra: norm, dot
+using Statistics: mean
 using Printf: @sprintf
 
 export  PartZero,
@@ -43,6 +44,7 @@ HV(v) = vcat(v, 1)
 
 include("partzeros.jl")
 include("geometries.jl")
+include("tolerances.jl")
 #include("optimizationproblem.jl")
 #include("optimization.jl")
 #include("resultevaluation.jl")
diff --git a/src/tolerances.jl b/src/tolerances.jl
new file mode 100644
index 0000000..aea3268
--- /dev/null
+++ b/src/tolerances.jl
@@ -0,0 +1,90 @@
+"""Supertype of tolerance types."""
+abstract type LocalizationToleranceType end
+
+struct PlanePlaneDistance <: LocalizationToleranceType end
+
+struct PlaneAxisDistance <: LocalizationToleranceType end
+
+struct AxisAxisDistance <: LocalizationToleranceType 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
+
+"""
+    transformmachined2datum(feature, points)
+
+Transform a list of points with the part zero of `feature`.
+"""
+function transformmachined2datum(feature, points)
+    pz = getpartzero(feature)
+    M = getpartzeroHM(pz)
+    newpoints = (M*HV(p) for p in points)
+    resultpoints = [p[1:3] for p in newpoints]
+    return resultpoints
+end
+
+function toleranceddistance(type::PlanePlaneDistance, f1, m1, f2, m2)
+    #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
+
+
+    # sokszor fog kelleni az, hogy v_feature = ismachined ? v_m : v_r
+    # 1. kelleni fog pont vagy pontok -> mindig pontlista vs pontlista legyen
+    # 2. kelleni fog, hogy kell-e
+    v_f1_ = getfeaturepoints(f1)
+    v_f2_ = getfeaturepoints(f2)
+    v_f1 = m1 == MACHINED ? transformmachined2datum(f1, v_f1_) : v_f1_
+    v_f2 = m2 == MACHINED ? transformmachined2datum(f2, v_f2_) : v_f2_
+
+    # pairwise distance
+    distance_vectors = (v2-v1 for v1 in v_f1 for v2 in v_f2)
+    mean_vector = mean(distance_vectors)
+    d = abs(dot(mean_vector, zaxis1))
+    return d
+end
+
+
+
+
+struct LocalizationTolerance
+    feature1::LocalizationFeature
+    machined1::IsMachined
+    feature2::LocalizationFeature
+    machined2::IsMachined
+    type::LocalizationToleranceType
+    nominalvalue::Float64
+    lowervalue::Float64
+    uppervalue::Float64
+    note::String
+end
+
+
+function toleranceddistance(t::LocalizationTolerance)
+    toleranceddistance(t.tolerancetype, t.feature1, t.machined1, t.feature2, t.machined2)
+end
\ No newline at end of file

From c0c726a966cc2bb8b1c45290043fd5858b3e7781 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tam=C3=A1s=20Cserteg?= <csertegt@gmail.com>
Date: Fri, 7 Jul 2023 16:11:34 +0200
Subject: [PATCH 06/27] fix bugs in tolerance calculation

---
 src/geometries.jl | 17 +++++++++++++++--
 src/tolerances.jl | 17 +++++++----------
 2 files changed, 22 insertions(+), 12 deletions(-)

diff --git a/src/geometries.jl b/src/geometries.jl
index 83e11b5..76d0da2 100644
--- a/src/geometries.jl
+++ b/src/geometries.jl
@@ -91,6 +91,12 @@ end
 
 GeometryStyle(::Type{MeshHole}) = IsFreeForm()
 
+function surfacepoints(::IsFreeForm, x::MeshHole)
+    points = vertices(x.surface)
+    verts = [x.coords for x in points]
+    return verts
+end
+
 function filteredsurfacepoints(::IsFreeForm, x::MeshHole)
     return x.convexhull
 end
@@ -104,12 +110,19 @@ struct MeshPlane <: AbstractPlaneGeometry
     surface::SimpleMesh
 end
 
+GeometryStyle(::Type{MeshPlane}) = IsFreeForm()
+
+function surfacepoints(::IsFreeForm, x::MeshPlane)
+    points = vertices(x.surface)
+    verts = [x.coords for x in points]
+    return verts
+end
+
 function filteredsurfacepoints(::IsFreeForm, x::MeshPlane)
     bbox = boundingbox(x.surface)
     return [bbox.min.coords, bbox.max.coords]
 end
 
-GeometryStyle(::Type{MeshPlane}) = IsFreeForm()
 
 """
     LocalizationFeature{R,M}
@@ -143,7 +156,7 @@ function Base.show(io::IO, lf::LocalizationFeature)
 end
 
 getfeaturename(f::LocalizationFeature) = f.name
-getpartzero(f::LocalizationFeature) = getpartzero(f.partzero)
+getpartzero(f::LocalizationFeature) = f.partzero
 getpartzeroname(f::LocalizationFeature) = getpartzeroname(f.partzero)
 
 getroughfeaturepoint(f::LocalizationFeature) = featurepoint(f.rough)
diff --git a/src/tolerances.jl b/src/tolerances.jl
index aea3268..3f2abfa 100644
--- a/src/tolerances.jl
+++ b/src/tolerances.jl
@@ -53,19 +53,16 @@ function toleranceddistance(type::PlanePlaneDistance, f1, m1, f2, m2)
             Can't compute `PlanePlaneDistance`")
     end
 
-
-    # sokszor fog kelleni az, hogy v_feature = ismachined ? v_m : v_r
-    # 1. kelleni fog pont vagy pontok -> mindig pontlista vs pontlista legyen
-    # 2. kelleni fog, hogy kell-e
-    v_f1_ = getfeaturepoints(f1)
-    v_f2_ = getfeaturepoints(f2)
+    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
-    distance_vectors = (v2-v1 for v1 in v_f1 for v2 in v_f2)
-    mean_vector = mean(distance_vectors)
-    d = abs(dot(mean_vector, zaxis1))
+    #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
 
@@ -86,5 +83,5 @@ end
 
 
 function toleranceddistance(t::LocalizationTolerance)
-    toleranceddistance(t.tolerancetype, t.feature1, t.machined1, t.feature2, t.machined2)
+    toleranceddistance(t.type, t.feature1, t.machined1, t.feature2, t.machined2)
 end
\ No newline at end of file

From 325be70b1022f33d28ab9ae7d54cc895acf4dca3 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tam=C3=A1s=20Cserteg?= <csertegt@gmail.com>
Date: Fri, 7 Jul 2023 16:11:52 +0200
Subject: [PATCH 07/27] get started with testing tolerances

---
 Project.toml                 |  3 ++-
 src/BlankLocalizationCore.jl |  8 ++++++
 test/runtests.jl             |  7 +++--
 test/tolerances.jl           | 50 ++++++++++++++++++++++++++++++++++++
 4 files changed, 65 insertions(+), 3 deletions(-)
 create mode 100644 test/tolerances.jl

diff --git a/Project.toml b/Project.toml
index ce8e89f..f8c34f7 100644
--- a/Project.toml
+++ b/Project.toml
@@ -21,6 +21,7 @@ julia = "1.6"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Meshes = "eacbb407-ea5a-433e-ab97-5258b1ca43fa"
 
 [targets]
-test = ["Test"]
+test = ["Test", "Meshes"]
diff --git a/src/BlankLocalizationCore.jl b/src/BlankLocalizationCore.jl
index 1d66304..f47f046 100644
--- a/src/BlankLocalizationCore.jl
+++ b/src/BlankLocalizationCore.jl
@@ -19,6 +19,14 @@ export  SimpleHole,
         LocalizationFeature,
         HoleLocalizationFeature,
         PlaneLocalizationFeature
+
+export  PlanePlaneDistance,
+        PlaneAxisDistance,
+        AxisAxisDistance,
+        AxisAxisConcentric,
+        LocalizationTolerance,
+        toleranceddistance
+
 #=
 export  OptimizationResult,
         Tolerance,
diff --git a/test/runtests.jl b/test/runtests.jl
index 33787e3..8fbc2ca 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,6 +1,9 @@
 using BlankLocalizationCore
 using Test
+using Meshes
 
-@testset "BlankLocalizationCore.jl" begin
-    # Write your tests here.
+const BLC = BlankLocalizationCore
+
+@testset "BlankLocalizationCore.jl: Tolerances" begin
+    include("tolerances.jl")
 end
diff --git a/test/tolerances.jl b/test/tolerances.jl
new file mode 100644
index 0000000..5a06927
--- /dev/null
+++ b/test/tolerances.jl
@@ -0,0 +1,50 @@
+@testset "PlanePlaneDistance: tolerancedistance" begin
+    empty_plane = SimplePlane([0,0,0])
+
+    machined_plane1 = SimplePlane([0,0,0])
+    machined_plane2 = SimplePlane([5,6,7])
+
+    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 = LocalizationFeature("t1", pz1, empty_plane, machined_plane1)
+    # lf2 has machined_plane2 as machined feature
+    lf2 = LocalizationFeature("tf1", pz1, empty_plane, machined_plane2)
+    # lf3 has "wrong" z axis
+    lf3 = LocalizationFeature("t3", pz2, empty_plane, empty_plane)
+    # lf4 has both features a mesh, and "correct" part zero
+    lf4 = 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
\ No newline at end of file

From 50aa346e7a28a050b9989d20e69db1118b109336 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tam=C3=A1s=20Cserteg?= <csertegt@gmail.com>
Date: Mon, 10 Jul 2023 14:23:02 +0200
Subject: [PATCH 08/27] implement PlaneAxisDistance type tolerance

---
 src/tolerances.jl  | 33 +++++++++++++++++++++++++++++
 test/tolerances.jl | 52 ++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 85 insertions(+)

diff --git a/src/tolerances.jl b/src/tolerances.jl
index 3f2abfa..9bb83b1 100644
--- a/src/tolerances.jl
+++ b/src/tolerances.jl
@@ -66,6 +66,39 @@ function toleranceddistance(type::PlanePlaneDistance, f1, m1, f2, m2)
     return d
 end
 
+function toleranceddistance(type::PlaneAxisDistance, f1, m1, f2, m2)
+    # 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
+
 
 
 
diff --git a/test/tolerances.jl b/test/tolerances.jl
index 5a06927..20e3a0f 100644
--- a/test/tolerances.jl
+++ b/test/tolerances.jl
@@ -47,4 +47,56 @@
     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])
+
+    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 = LocalizationFeature("t1", pz1_front, machined_plane1, machined_plane1)
+    # lf2 simple plane and mash plane
+    lf2 = 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])
+    pl2 = SimplePlane([0, 0, 15])
+    empty_plane = SimplePlane([0,0,0])
+
+    h1 = LocalizationFeature("fronth1", pz1_front, empty_hole, sh1)
+    h2 = LocalizationFeature("fronth1", pz1_front, empty_hole, sh2)
+
+    p1 = LocalizationFeature("rightp1", pz2_right, empty_plane, pl1)
+    p2 = 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
\ No newline at end of file

From 6bb530691365c658fd0861be9d9262f9c422147f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tam=C3=A1s=20Cserteg?= <csertegt@gmail.com>
Date: Mon, 10 Jul 2023 14:57:49 +0200
Subject: [PATCH 09/27] implement AxisAxisDistance tolerance

---
 src/BlankLocalizationCore.jl |  2 +-
 src/tolerances.jl            | 39 +++++++++++++++++++++++++++++
 test/tolerances.jl           | 48 ++++++++++++++++++++++++++++++++++++
 3 files changed, 88 insertions(+), 1 deletion(-)

diff --git a/src/BlankLocalizationCore.jl b/src/BlankLocalizationCore.jl
index f47f046..1b0cd4f 100644
--- a/src/BlankLocalizationCore.jl
+++ b/src/BlankLocalizationCore.jl
@@ -5,7 +5,7 @@ using DataFrames: DataFrame, names, nrow
 using PrettyTables: pretty_table, ft_nonothing, tf_html_minimalist
 using Meshes: SimpleMesh, vertices, boundingbox
 using Logging: @warn
-using LinearAlgebra: norm, dot
+using LinearAlgebra: norm, dot, cross
 using Statistics: mean
 using Printf: @sprintf
 
diff --git a/src/tolerances.jl b/src/tolerances.jl
index 9bb83b1..4159970 100644
--- a/src/tolerances.jl
+++ b/src/tolerances.jl
@@ -99,6 +99,45 @@ function toleranceddistance(type::PlaneAxisDistance, f1, m1, f2, m2)
     return distance
 end
 
+function toleranceddistance(type::AxisAxisDistance, f1, m1, f2, m2)
+    # 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)
+
+    # TODO
+    # what if the two axes are parallel?
+    zaxis1 = zaxis(getpartzero(f1))
+    zaxis2 = zaxis(getpartzero(f2))
+
+    perpend_v = cross(zaxis1, zaxis2)
+
+    distancev = fp1 - fp2
+    distance = abs(dot(perpend_v, distancev))
+
+    return distance
+end
+
 
 
 
diff --git a/test/tolerances.jl b/test/tolerances.jl
index 20e3a0f..d59f480 100644
--- a/test/tolerances.jl
+++ b/test/tolerances.jl
@@ -99,4 +99,52 @@ end
     @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])
+
+    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 = LocalizationFeature("t1", pz1_front, machined_plane1, machined_plane1)
+    # lf2 simple plane and mesh plane
+    lf2 = LocalizationFeature("tf1", pz1_front, machined_plane1, mesh_plane1)
+
+    t1_mr = LocalizationTolerance(lf1, BLC.ROUGH, lf2, BLC.MACHINED, AxisAxisDistance(), 0, 0, 0, "")
+    t1_rm = LocalizationTolerance(lf2, BLC.MACHINED, lf1, BLC.ROUGH, AxisAxisDistance(), 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(), 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 = LocalizationFeature("fronth1", pz1_front, sh1_rough, sh1_machined)
+
+    rh1 = 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(), 20, 20, 20, "")
+    aat2 = LocalizationTolerance(fh1, BLC.MACHINED, rh1, BLC.ROUGH, AxisAxisDistance(), 20, 20, 20, "")
+    aat3 = LocalizationTolerance(fh1, BLC.ROUGH, rh1, BLC.MACHINED, AxisAxisDistance(), 20, 20, 20, "")
+    aat4 = LocalizationTolerance(fh1, BLC.ROUGH, rh1, BLC.ROUGH, AxisAxisDistance(), 20, 20, 20, "")
+    @test isapprox(toleranceddistance(aat1), 20)
+    @test isapprox(toleranceddistance(aat2), 20)
+    @test isapprox(toleranceddistance(aat3), 20)
+    @test isapprox(toleranceddistance(aat4), 20)
 end
\ No newline at end of file

From dc722c0b46308915620939c88a5ca2d64da002b7 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tam=C3=A1s=20Cserteg?= <csertegt@gmail.com>
Date: Mon, 10 Jul 2023 15:47:11 +0200
Subject: [PATCH 10/27] add AxisAxisConcentric tolerance

---
 src/BlankLocalizationCore.jl |  2 +-
 src/tolerances.jl            | 42 ++++++++++++++++++++++++++++++++++++
 test/tolerances.jl           | 29 ++++++++++++++++++++++++-
 3 files changed, 71 insertions(+), 2 deletions(-)

diff --git a/src/BlankLocalizationCore.jl b/src/BlankLocalizationCore.jl
index 1b0cd4f..3e2f120 100644
--- a/src/BlankLocalizationCore.jl
+++ b/src/BlankLocalizationCore.jl
@@ -5,7 +5,7 @@ using DataFrames: DataFrame, names, nrow
 using PrettyTables: pretty_table, ft_nonothing, tf_html_minimalist
 using Meshes: SimpleMesh, vertices, boundingbox
 using Logging: @warn
-using LinearAlgebra: norm, dot, cross
+using LinearAlgebra: norm, dot, cross, inv
 using Statistics: mean
 using Printf: @sprintf
 
diff --git a/src/tolerances.jl b/src/tolerances.jl
index 4159970..e6191e0 100644
--- a/src/tolerances.jl
+++ b/src/tolerances.jl
@@ -138,7 +138,49 @@ function toleranceddistance(type::AxisAxisDistance, f1, m1, f2, m2)
     return distance
 end
 
+function toleranceddistance(type::AxisAxisConcentric, f1, m1, f2, m2)
+    # 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
+    gs1 = GeometryStyle(typeof(geom1))
+    gs2 = GeometryStyle(typeof(geom2))
+
+    if (gs1) != IsPrimitive() || gs2 != IsPrimitive()
+        error("AxisAxisConcentric is only implemented when both features are IsPrimitive!
+        f1 is $gs1, f2 is $gs2")
+    end
+
+    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
+
+    #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
 
 
 struct LocalizationTolerance
diff --git a/test/tolerances.jl b/test/tolerances.jl
index d59f480..ae4c63c 100644
--- a/test/tolerances.jl
+++ b/test/tolerances.jl
@@ -147,4 +147,31 @@ end
     @test isapprox(toleranceddistance(aat2), 20)
     @test isapprox(toleranceddistance(aat3), 20)
     @test isapprox(toleranceddistance(aat4), 20)
-end
\ No newline at end of file
+end
+
+@testset "AxisAxisConcentric: 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 = LocalizationFeature("righth1", pz2_right, aac1_r, aac1_m)
+    aac_lf2 = 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

From 19aa14dd9eb4144f44148bc8eedbc7c024b77319 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tam=C3=A1s=20Cserteg?= <csertegt@gmail.com>
Date: Wed, 12 Jul 2023 11:28:59 +0200
Subject: [PATCH 11/27] add primitive-freeform axis-axis tolerance

---
 src/geometries.jl  |  13 ++++
 src/tolerances.jl  | 146 ++++++++++++++++++++++++++++++---------------
 test/tolerances.jl |  35 ++++++++++-
 3 files changed, 146 insertions(+), 48 deletions(-)

diff --git a/src/geometries.jl b/src/geometries.jl
index 76d0da2..2d51a31 100644
--- a/src/geometries.jl
+++ b/src/geometries.jl
@@ -173,3 +173,16 @@ function getmachinedfeaturepointindatum(f::LocalizationFeature)
     v_indatum = T*HV(v)
     return v_indatum[1:3]
 end
+
+"""
+    transformmachined2datum(feature, points)
+
+Transform a list of points with the part zero of `feature`.
+"""
+function transformmachined2datum(feature, points)
+    pz = getpartzero(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/tolerances.jl b/src/tolerances.jl
index e6191e0..3c169b7 100644
--- a/src/tolerances.jl
+++ b/src/tolerances.jl
@@ -31,20 +31,26 @@ getfeaturepoints(::IsFreeForm, x) = surfacepoints(x)
 """Enum to store if a feature should be handled as rough or machined."""
 @enum IsMachined MACHINED ROUGH
 
-"""
-    transformmachined2datum(feature, points)
 
-Transform a list of points with the part zero of `feature`.
-"""
-function transformmachined2datum(feature, points)
-    pz = getpartzero(feature)
-    M = getpartzeroHM(pz)
-    newpoints = (M*HV(p) for p in points)
-    resultpoints = [p[1:3] for p in newpoints]
-    return resultpoints
+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, f1, m1, f2, m2)
+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))
@@ -66,7 +72,11 @@ function toleranceddistance(type::PlanePlaneDistance, f1, m1, f2, m2)
     return d
 end
 
-function toleranceddistance(type::PlaneAxisDistance, f1, m1, f2, m2)
+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
@@ -99,7 +109,11 @@ function toleranceddistance(type::PlaneAxisDistance, f1, m1, f2, m2)
     return distance
 end
 
-function toleranceddistance(type::AxisAxisDistance, f1, m1, f2, m2)
+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
@@ -138,28 +152,7 @@ function toleranceddistance(type::AxisAxisDistance, f1, m1, f2, m2)
     return distance
 end
 
-function toleranceddistance(type::AxisAxisConcentric, f1, m1, f2, m2)
-    # 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
-    gs1 = GeometryStyle(typeof(geom1))
-    gs2 = GeometryStyle(typeof(geom2))
-
-    if (gs1) != IsPrimitive() || gs2 != IsPrimitive()
-        error("AxisAxisConcentric is only implemented when both features are IsPrimitive!
-        f1 is $gs1, f2 is $gs2")
-    end
-
-    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
-
+function aac_primitive_primitive(f1, m1, f2, m2)
     #check if axes are parallel
     zaxis1 = zaxis(getpartzero(f1))
     zaxis2 = zaxis(getpartzero(f2))
@@ -182,20 +175,79 @@ function toleranceddistance(type::AxisAxisConcentric, f1, m1, f2, m2)
     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
 
-struct LocalizationTolerance
-    feature1::LocalizationFeature
-    machined1::IsMachined
-    feature2::LocalizationFeature
-    machined2::IsMachined
-    type::LocalizationToleranceType
-    nominalvalue::Float64
-    lowervalue::Float64
-    uppervalue::Float64
-    note::String
+    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.feature1, t.machined1, t.feature2, t.machined2)
-end
\ No newline at end of file
+    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/test/tolerances.jl b/test/tolerances.jl
index ae4c63c..586451f 100644
--- a/test/tolerances.jl
+++ b/test/tolerances.jl
@@ -149,7 +149,7 @@ end
     @test isapprox(toleranceddistance(aat4), 20)
 end
 
-@testset "AxisAxisConcentric: tolerancedistance" begin
+@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])
     
@@ -175,3 +175,36 @@ end
     @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 = 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 = 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 = 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

From e4b447bd45d45948fc8bdf94e65871e4ca127574 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tam=C3=A1s=20Cserteg?= <csertegt@gmail.com>
Date: Thu, 13 Jul 2023 09:45:35 +0200
Subject: [PATCH 12/27] allow to override projection axis in AxisAxisDistance

---
 src/tolerances.jl | 14 ++++++++++----
 1 file changed, 10 insertions(+), 4 deletions(-)

diff --git a/src/tolerances.jl b/src/tolerances.jl
index 3c169b7..af0421b 100644
--- a/src/tolerances.jl
+++ b/src/tolerances.jl
@@ -5,7 +5,11 @@ struct PlanePlaneDistance <: LocalizationToleranceType end
 
 struct PlaneAxisDistance <: LocalizationToleranceType end
 
-struct AxisAxisDistance <: LocalizationToleranceType end
+struct AxisAxisDistance <: LocalizationToleranceType
+    projectionaxis::Union{Nothing,Vector{Float64}}
+end
+
+AxisAxisDistance() = AxisAxisDistance(nothing)
 
 struct AxisAxisConcentric <: LocalizationToleranceType end
 
@@ -109,6 +113,10 @@ function toleranceddistance(type::PlaneAxisDistance, t::LocalizationTolerance)
     return distance
 end
 
+#TODO
+# what happens if both axes are parallel???
+# then which axis is used to check the distance?
+# -> an axis can be defined
 function toleranceddistance(type::AxisAxisDistance, t::LocalizationTolerance)
     f1 = t.feature1
     m1 = t.machined1
@@ -139,12 +147,10 @@ function toleranceddistance(type::AxisAxisDistance, t::LocalizationTolerance)
     fp1 = m1 == MACHINED ? getmachinedfeaturepointindatum(f1) : getroughfeaturepoint(f1)
     fp2 = m2 == MACHINED ? getmachinedfeaturepointindatum(f2) : getroughfeaturepoint(f2)
 
-    # TODO
-    # what if the two axes are parallel?
     zaxis1 = zaxis(getpartzero(f1))
     zaxis2 = zaxis(getpartzero(f2))
 
-    perpend_v = cross(zaxis1, zaxis2)
+    perpend_v = isnothing(type.projectionaxis) ? cross(zaxis1, zaxis2) : type.projectionaxis
 
     distancev = fp1 - fp2
     distance = abs(dot(perpend_v, distancev))

From cce4065ebade0f89bb62dbf3ae2c06ab9f13702d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tam=C3=A1s=20Cserteg?= <csertegt@gmail.com>
Date: Thu, 13 Jul 2023 10:44:36 +0200
Subject: [PATCH 13/27] adjust optimization problem to recent changes

---
 src/BlankLocalizationCore.jl |  8 ++-
 src/optimizationproblem.jl   | 94 ++++++++++--------------------------
 2 files changed, 28 insertions(+), 74 deletions(-)

diff --git a/src/BlankLocalizationCore.jl b/src/BlankLocalizationCore.jl
index 3e2f120..e2cc70c 100644
--- a/src/BlankLocalizationCore.jl
+++ b/src/BlankLocalizationCore.jl
@@ -27,11 +27,9 @@ export  PlanePlaneDistance,
         LocalizationTolerance,
         toleranceddistance
 
-#=
 export  OptimizationResult,
-        Tolerance,
         MultiOperationProblem
-
+#=
 export  createjumpmodel,
         setjumpresult!,
         optimizeproblem!
@@ -53,8 +51,8 @@ HV(v) = vcat(v, 1)
 include("partzeros.jl")
 include("geometries.jl")
 include("tolerances.jl")
-#include("optimizationproblem.jl")
-#include("optimization.jl")
+include("optimizationproblem.jl")
+include("optimization.jl")
 #include("resultevaluation.jl")
 
 end
diff --git a/src/optimizationproblem.jl b/src/optimizationproblem.jl
index a4d8fe6..1fc12d2 100644
--- a/src/optimizationproblem.jl
+++ b/src/optimizationproblem.jl
@@ -14,48 +14,39 @@ end
 
 emptyor() = OptimizationResult("empty", 0.0)
 
-struct Tolerance
-    featurename1::String
-    ismachined1::Bool
-    projection::Function
-    featurename2::String
-    ismachined2::Bool
-    nominalvalue::Float64
-    lowervalue::Float64
-    uppervalue::Float64
-    note::String
-end
-
 mutable struct MultiOperationProblem
     partzeros::Vector{PartZero}
-    holes::Vector{HoleLocalizationFeature}
-    planes::Vector{PlaneLocalizationFeature}
-    tolerances::Vector{Tolerance}
+    features::Vector{LocalizationFeature} # features that have rough and machined -> allowanced
+    tolerances::Vector{LocalizationTolerance}
     parameters::Dict{String,Real}
     opresult::OptimizationResult
 end
 
-function MultiOperationProblem(partzeros, holes, planes, tolerances, parameters)
-    return MultiOperationProblem(partzeros, holes, planes, tolerances, parameters, emptyor())
+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: 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
+    featuretypes = GeometryStyle.(typeof.(x.rough for x in mop.features))
+    for ft in featuretypes
+        ft === IsFreeForm() && return :HybridProblem
     end
     return :PrimitiveProblem
 end
 
+function collectholefeatures(mop::MultiOperationProblem)
+    filter(x->x isa HoleLocalizationFeature, mop.features)
+end
+
+function collectplanefeatures(mop::MultiOperationProblem)
+    filter(x->x isa PlaneLocalizationFeature, mop.features)
+end
+
 function Base.show(io::IO, mop::MultiOperationProblem)
-    nh = size(mop.holes, 1)
-    np = size(mop.planes, 1)
+    nh = size(collectholefeatures(mop), 1)
+    np = size(collectplanefeatures(mop), 1)
     npz = size(mop.partzeros, 1)
     nts = size(mop.tolerances, 1)
     sn = string(problemtype(mop))
@@ -92,54 +83,19 @@ It is assumed that all features have distinct names.
 Return `nothing`, if no feature is found with `featurename`.
 """
 function getfeaturebyname(mop::MultiOperationProblem, featurename)
-    function retbyname(array, name)
-        for f in array
-            if getfeaturename(f) == name
-                return f
-            end
+    for f in mop.features
+        if getfeaturename(f) == featurename
+            return f
         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)
+    return nothing
 end
 
 """
-    collectholesbypartzero(mop::MultiOperationProblem, partzeroname)
+    collectfeaturesbypartzero(mop::MultiOperationProblem, partzeroname)
 
-Collect holes that are grouped to part zero called `partzeroname`.
+Collect features that are grouped to part zero called `partzeroname`.
 """
-function getholesbypartzero(mop::MultiOperationProblem, partzeroname)
-    return filter(x->getpartzeroname(x) == partzeroname, mop.holes)
+function collectfeaturesbypartzero(mop::MultiOperationProblem, partzeroname)
+    return filter(x->getpartzeroname(x) == partzeroname, mop.features)
 end
-
-"""
-    collectmachinedholes(mop::MultiOperationProblem)
-
-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)

From 52e819f0eba381d40c497e74cc17e05e396ed219 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tam=C3=A1s=20Cserteg?= <csertegt@gmail.com>
Date: Mon, 13 Nov 2023 18:40:48 +0100
Subject: [PATCH 14/27] remplace PlaneAndNormal with SimplePlane

---
 docs/src/example.md          | 12 ++++++------
 examples/example.jl          | 12 ++++++------
 src/BlankLocalizationCore.jl |  1 -
 src/geometries.jl            | 23 ++---------------------
 test/testproblem.jl          | 12 ++++++------
 5 files changed, 20 insertions(+), 40 deletions(-)

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 a7bc590..a9b382d 100644
--- a/src/BlankLocalizationCore.jl
+++ b/src/BlankLocalizationCore.jl
@@ -26,7 +26,6 @@ export  AbstractHoleGeometry,
         visualizationgeometry,
         SimpleHole,
         SimplePlane,
-        PlaneAndNormal,
         MeshHole,
         MeshPlane,
         LocalizationFeature,
diff --git a/src/geometries.jl b/src/geometries.jl
index 46f7a2b..6098797 100644
--- a/src/geometries.jl
+++ b/src/geometries.jl
@@ -113,11 +113,11 @@ end
 """
 SimplePlane <: AbstractPlaneGeometry
 
-A simple plane structure with one point.
-Normal vector of the plane is defined by its partzero taken from the feature descriptor.
+A simple plane structure with one point and a normal vector.
 """
 struct SimplePlane <: AbstractPlaneGeometry
     p::Vector{Float64}
+    n::Vector{Float64}
 end
 
 GeometryStyle(::Type{SimplePlane}) = IsPrimitive()
@@ -135,25 +135,6 @@ function rectangleforplane(point, v1 ,v2, sidelength)
 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)
diff --git a/test/testproblem.jl b/test/testproblem.jl
index f340601..f593bf8 100644
--- a/test/testproblem.jl
+++ b/test/testproblem.jl
@@ -29,19 +29,19 @@
     ## 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
 

From 9fdc58016dff2af46311b7c09e5544848900a735 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tam=C3=A1s=20Cserteg?= <csertegt@gmail.com>
Date: Mon, 13 Nov 2023 18:46:46 +0100
Subject: [PATCH 15/27] don't export LocalizationFeature

---
 src/BlankLocalizationCore.jl |  1 -
 test/tolerances.jl           | 38 ++++++++++++++++++------------------
 2 files changed, 19 insertions(+), 20 deletions(-)

diff --git a/src/BlankLocalizationCore.jl b/src/BlankLocalizationCore.jl
index a9b382d..4e26241 100644
--- a/src/BlankLocalizationCore.jl
+++ b/src/BlankLocalizationCore.jl
@@ -28,7 +28,6 @@ export  AbstractHoleGeometry,
         SimplePlane,
         MeshHole,
         MeshPlane,
-        LocalizationFeature,
         HoleLocalizationFeature,
         PlaneLocalizationFeature
         
diff --git a/test/tolerances.jl b/test/tolerances.jl
index 586451f..c810054 100644
--- a/test/tolerances.jl
+++ b/test/tolerances.jl
@@ -12,13 +12,13 @@
     pz2 = PartZero("testpz2", [0,0,0], [0 0 1;0 1 0;0 1 0])
 
     # lf1 has machined_plane1 as machined feature
-    lf1 = LocalizationFeature("t1", pz1, empty_plane, machined_plane1)
+    lf1 = BLC.LocalizationFeature("t1", pz1, empty_plane, machined_plane1)
     # lf2 has machined_plane2 as machined feature
-    lf2 = LocalizationFeature("tf1", pz1, empty_plane, machined_plane2)
+    lf2 = BLC.LocalizationFeature("tf1", pz1, empty_plane, machined_plane2)
     # lf3 has "wrong" z axis
-    lf3 = LocalizationFeature("t3", pz2, empty_plane, empty_plane)
+    lf3 = BLC.LocalizationFeature("t3", pz2, empty_plane, empty_plane)
     # lf4 has both features a mesh, and "correct" part zero
-    lf4 = LocalizationFeature("t4", pz1, mesh_plane1, mesh_plane1)
+    lf4 = BLC.LocalizationFeature("t4", pz1, mesh_plane1, mesh_plane1)
     
     # testing machined-machined distance
     # nominal, lower, upper, and note doesn't matter here
@@ -61,9 +61,9 @@ end
     pz2_right = PartZero("testpz2", [10,10,10], [-1 0 0;0 0 1;0 1 0])
 
     # lf1 simple plane and simple plane
-    lf1 = LocalizationFeature("t1", pz1_front, machined_plane1, machined_plane1)
+    lf1 = BLC.LocalizationFeature("t1", pz1_front, machined_plane1, machined_plane1)
     # lf2 simple plane and mash plane
-    lf2 = LocalizationFeature("tf1", pz1_front, machined_plane1, mesh_plane1)
+    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, "")
@@ -83,11 +83,11 @@ end
     pl2 = SimplePlane([0, 0, 15])
     empty_plane = SimplePlane([0,0,0])
 
-    h1 = LocalizationFeature("fronth1", pz1_front, empty_hole, sh1)
-    h2 = LocalizationFeature("fronth1", pz1_front, empty_hole, sh2)
+    h1 = BLC.LocalizationFeature("fronth1", pz1_front, empty_hole, sh1)
+    h2 = BLC.LocalizationFeature("fronth1", pz1_front, empty_hole, sh2)
 
-    p1 = LocalizationFeature("rightp1", pz2_right, empty_plane, pl1)
-    p2 = LocalizationFeature("rightp2", pz2_right, empty_plane, pl2)
+    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, "")
@@ -114,9 +114,9 @@ end
     pz2_right = PartZero("rightpz2", [10,10,10], [-1 0 0;0 0 1;0 1 0])
 
     # lf1 simple plane and simple plane
-    lf1 = LocalizationFeature("t1", pz1_front, machined_plane1, machined_plane1)
+    lf1 = BLC.LocalizationFeature("t1", pz1_front, machined_plane1, machined_plane1)
     # lf2 simple plane and mesh plane
-    lf2 = LocalizationFeature("tf1", pz1_front, machined_plane1, mesh_plane1)
+    lf2 = BLC.LocalizationFeature("tf1", pz1_front, machined_plane1, mesh_plane1)
 
     t1_mr = LocalizationTolerance(lf1, BLC.ROUGH, lf2, BLC.MACHINED, AxisAxisDistance(), 0, 0, 0, "")
     t1_rm = LocalizationTolerance(lf2, BLC.MACHINED, lf1, BLC.ROUGH, AxisAxisDistance(), 0, 0, 0, "")
@@ -134,9 +134,9 @@ end
     sh1_rough = SimpleHole([0, 10, 0], 14.9)
     sh2_rough = SimpleHole([5, 25, 20], 14.9)
 
-    fh1 = LocalizationFeature("fronth1", pz1_front, sh1_rough, sh1_machined)
+    fh1 = BLC.LocalizationFeature("fronth1", pz1_front, sh1_rough, sh1_machined)
 
-    rh1 = LocalizationFeature("righth1", pz2_right, sh2_rough, sh2_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(), 20, 20, 20, "")
@@ -164,8 +164,8 @@ end
     
     # distance aac2_m -> aac2_r: sqrt(2)
     
-    aac_lf1 = LocalizationFeature("righth1", pz2_right, aac1_r, aac1_m)
-    aac_lf2 = LocalizationFeature("righth2", pz3_right, aac2_r, aac2_m)
+    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, "")
@@ -186,7 +186,7 @@ end
     aac_c1 = [connect((1,2,3))]
     aac_mesh1 = MeshHole(SimpleMesh(aac_ps1, aac_c1), [])
 
-    aacpf_lf1 = LocalizationFeature("h1", pz_aac, aac_mesh1, aacpf_h)
+    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)
 
@@ -195,7 +195,7 @@ end
     sm2 = SimpleMesh(aac_ps2, aac_c1)
     aac_mesh2 = MeshHole(sm2, [])
 
-    aacpf_lf2 = LocalizationFeature("h1", pz_aac, aac_mesh2, aacpf_h)
+    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)
 
@@ -204,7 +204,7 @@ end
     sm3 = SimpleMesh(aac_ps3, aac_c1)
     aac_mesh3 = MeshHole(sm3, [])
 
-    aacpf_lf3 = LocalizationFeature("h1", pz_aac, aac_mesh3, aacpf_h)
+    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

From 999e5717c84a08faa70c4c0648803b1caa878f35 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tam=C3=A1s=20Cserteg?= <csertegt@gmail.com>
Date: Mon, 13 Nov 2023 18:48:42 +0100
Subject: [PATCH 16/27] polishing look and fixing typo

---
 src/geometries.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/geometries.jl b/src/geometries.jl
index 6098797..88ecd16 100644
--- a/src/geometries.jl
+++ b/src/geometries.jl
@@ -82,7 +82,7 @@ end
 
 filteredsurfacepoints(x::T) where {T} = filteredsurfacepoints(GeometryStyle(T), x)
 function filteredsurfacepoints(::IsPrimitive, x)
-    error("Function `surfacepoints` is not defined for `IsPrimitive`` features")
+    error("Function `filteredsurfacepoints` is not defined for `IsPrimitive`` features")
 end
 
 featureradius(x::T) where {T<:AbstractHoleGeometry} = featureradius(GeometryStyle(T), x)
@@ -232,8 +232,8 @@ getpartzeroname(f::LocalizationFeature) = getpartzeroname(f.partzero)
 
 getroughfeaturepoint(f::LocalizationFeature) = featurepoint(f.rough)
 getmachinedfeaturepoint(f::LocalizationFeature) = featurepoint(f.machined)
-getmachinedradius(f::LocalizationFeature) = featureradius(f.machined)
 getroughradius(f::LocalizationFeature) = featureradius(f.rough)
+getmachinedradius(f::LocalizationFeature) = featureradius(f.machined)
 
 getroughfilteredpoints(f::LocalizationFeature) = filteredsurfacepoints(f.rough)
 

From 3e01c288a32078111409af91085501b47084ac7e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tam=C3=A1s=20Cserteg?= <csertegt@gmail.com>
Date: Mon, 13 Nov 2023 18:52:10 +0100
Subject: [PATCH 17/27] AxisAxisDistance must have a projectionaxis

---
 src/tolerances.jl | 15 ++-------------
 1 file changed, 2 insertions(+), 13 deletions(-)

diff --git a/src/tolerances.jl b/src/tolerances.jl
index af0421b..543ab5f 100644
--- a/src/tolerances.jl
+++ b/src/tolerances.jl
@@ -6,11 +6,9 @@ struct PlanePlaneDistance <: LocalizationToleranceType end
 struct PlaneAxisDistance <: LocalizationToleranceType end
 
 struct AxisAxisDistance <: LocalizationToleranceType
-    projectionaxis::Union{Nothing,Vector{Float64}}
+    projectionaxis::Vector{Float64}
 end
 
-AxisAxisDistance() = AxisAxisDistance(nothing)
-
 struct AxisAxisConcentric <: LocalizationToleranceType end
 
 """
@@ -113,10 +111,6 @@ function toleranceddistance(type::PlaneAxisDistance, t::LocalizationTolerance)
     return distance
 end
 
-#TODO
-# what happens if both axes are parallel???
-# then which axis is used to check the distance?
-# -> an axis can be defined
 function toleranceddistance(type::AxisAxisDistance, t::LocalizationTolerance)
     f1 = t.feature1
     m1 = t.machined1
@@ -147,13 +141,8 @@ function toleranceddistance(type::AxisAxisDistance, t::LocalizationTolerance)
     fp1 = m1 == MACHINED ? getmachinedfeaturepointindatum(f1) : getroughfeaturepoint(f1)
     fp2 = m2 == MACHINED ? getmachinedfeaturepointindatum(f2) : getroughfeaturepoint(f2)
 
-    zaxis1 = zaxis(getpartzero(f1))
-    zaxis2 = zaxis(getpartzero(f2))
-
-    perpend_v = isnothing(type.projectionaxis) ? cross(zaxis1, zaxis2) : type.projectionaxis
-
     distancev = fp1 - fp2
-    distance = abs(dot(perpend_v, distancev))
+    distance = abs(dot(type.projectionaxis, distancev))
 
     return distance
 end

From cdc443506c1610d260f11e203c8863e60478a284 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tam=C3=A1s=20Cserteg?= <csertegt@gmail.com>
Date: Mon, 13 Nov 2023 19:06:59 +0100
Subject: [PATCH 18/27] fix tests after changing SimplePlane and
 AxisAxisDistance signature

---
 test/geometries.jl  |  2 +-
 test/testproblem.jl | 12 ++++++------
 test/tolerances.jl  | 30 +++++++++++++++---------------
 3 files changed, 22 insertions(+), 22 deletions(-)

diff --git a/test/geometries.jl b/test/geometries.jl
index 4c93c01..c48357b 100644
--- a/test/geometries.jl
+++ b/test/geometries.jl
@@ -1,6 +1,6 @@
 @testset "IsPrimitive geometries" begin
     sh = SimpleHole([0, 0, 0], 29)
-    sp = SimplePlane([0, 0, 0])
+    sp = SimplePlane([0, 0, 0], [0, 0, 1])
 
     @test featurepoint(sh) == [0, 0, 0]
     @test featureradius(sh) == 29
diff --git a/test/testproblem.jl b/test/testproblem.jl
index f593bf8..a8774a5 100644
--- a/test/testproblem.jl
+++ b/test/testproblem.jl
@@ -12,19 +12,19 @@
     ## 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
 
diff --git a/test/tolerances.jl b/test/tolerances.jl
index c810054..0f63103 100644
--- a/test/tolerances.jl
+++ b/test/tolerances.jl
@@ -1,8 +1,8 @@
 @testset "PlanePlaneDistance: tolerancedistance" begin
-    empty_plane = SimplePlane([0,0,0])
+    empty_plane = SimplePlane([0,0,0], [-1, 0, 0])
 
-    machined_plane1 = SimplePlane([0,0,0])
-    machined_plane2 = SimplePlane([5,6,7])
+    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))])
@@ -51,7 +51,7 @@ end
 
 @testset "PlaneAxisDistance: tolerancedistance" begin
     # test error throwing -> should be addressed by #7
-    machined_plane1 = SimplePlane([0,0,0])
+    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))])
@@ -79,9 +79,9 @@ end
     sh2 = SimpleHole([0, -10, 0], 15)
     empty_hole = SimpleHole([0,0,0], 0)
 
-    pl1 = SimplePlane([15, 0, 0])
-    pl2 = SimplePlane([0, 0, 15])
-    empty_plane = SimplePlane([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)
@@ -103,7 +103,7 @@ end
 
 @testset "AxisAxisDistance: tolerancedistance" begin
     ## test error throwing -> should be addressed by #7
-    machined_plane1 = SimplePlane([0,0,0])
+    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))])
@@ -118,13 +118,13 @@ end
     # 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(), 0, 0, 0, "")
-    t1_rm = LocalizationTolerance(lf2, BLC.MACHINED, lf1, BLC.ROUGH, AxisAxisDistance(), 0, 0, 0, "")
+    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(), 0, 0, 0, "")
+    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
@@ -139,10 +139,10 @@ end
     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(), 20, 20, 20, "")
-    aat2 = LocalizationTolerance(fh1, BLC.MACHINED, rh1, BLC.ROUGH, AxisAxisDistance(), 20, 20, 20, "")
-    aat3 = LocalizationTolerance(fh1, BLC.ROUGH, rh1, BLC.MACHINED, AxisAxisDistance(), 20, 20, 20, "")
-    aat4 = LocalizationTolerance(fh1, BLC.ROUGH, rh1, BLC.ROUGH, AxisAxisDistance(), 20, 20, 20, "")
+    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)

From 4aa351e59defd3dc0a467482d5c61828b4e57234 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tam=C3=A1s=20Cserteg?= <csertegt@gmail.com>
Date: Mon, 13 Nov 2023 21:33:28 +0100
Subject: [PATCH 19/27] doc updates

---
 src/geometries.jl          |  5 ++---
 src/optimizationproblem.jl | 10 ++++++++++
 2 files changed, 12 insertions(+), 3 deletions(-)

diff --git a/src/geometries.jl b/src/geometries.jl
index 88ecd16..6fa7b97 100644
--- a/src/geometries.jl
+++ b/src/geometries.jl
@@ -201,9 +201,8 @@ end
 A feature that is machined and allowance can be computed for it.
 It has a name, a [`PartZero`](@ref), and a rough and machined geometry
 ([`AbstractLocalizationGeometry`](@ref)).
-The two geometries should be of same type (hole, plane, etc.),
-but [`HoleLocalizationFeature`](@ref) and [`PlaneLocalizationFeature`](@ref) enforce this
-property
+The two geometries should be of same type (hole, plane, etc.), but only
+[`HoleLocalizationFeature`](@ref) and [`PlaneLocalizationFeature`](@ref) enforce this property.
 """
 struct LocalizationFeature{R,M}
     name::String
diff --git a/src/optimizationproblem.jl b/src/optimizationproblem.jl
index 52a02f3..fd51549 100644
--- a/src/optimizationproblem.jl
+++ b/src/optimizationproblem.jl
@@ -23,6 +23,16 @@ 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.
+
+The elements of `features` can only be `HoleLocalizationFeature`s and `PlaneLocalizationFeature`s.
+Otherwise the JuMP model building fails.
+
+The features in tolerances ([`LocalizationTolerance`](@ref)) can be `nothing`, if needed.
+"""
 mutable struct MultiOperationProblem
     partzeros::Vector{PartZero}
     features::Vector{LocalizationFeature} # features that have rough and machined -> allowanced

From a5326a4bced333a843b20b803d9008ea0fdce013 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tam=C3=A1s=20Cserteg?= <csertegt@gmail.com>
Date: Fri, 19 Jan 2024 11:57:42 +0100
Subject: [PATCH 20/27] rename getpartzeroname function

---
 src/optimization.jl | 12 ++++++------
 src/partzeros.jl    |  2 +-
 2 files changed, 7 insertions(+), 7 deletions(-)

diff --git a/src/optimization.jl b/src/optimization.jl
index 5bc41e5..1ca9811 100644
--- a/src/optimization.jl
+++ b/src/optimization.jl
@@ -10,7 +10,7 @@ function addhole2model!(::IsPrimitive, model, hole, ipzmatricedict)
     # register distance variable:
     dxy = @variable(model, base_name = string("d_xy_", getfeaturename(hole)), lower_bound = 0.0)
 
-    pzn = getpartzeroname(hole)
+    pzn = partzeroname(hole)
     v_machined = getmachinedfeaturepoint(hole)
     v_rough = getroughfeaturepoint(hole)
     r_machined = getmachinedradius(hole)
@@ -35,7 +35,7 @@ function addhole2model!(::IsFreeForm, model, hole, ipzmatricedict)
     # register distance variable:
     dxy = @variable(model, [qiter], base_name = string("d_xy_", getfeaturename(hole)), lower_bound = 0.0)
 
-    pzn = getpartzeroname(hole)
+    pzn = partzeroname(hole)
     v_machined = getmachinedfeaturepoint(hole)
     r_machined = getmachinedradius(hole)
     # equation (4)
@@ -62,7 +62,7 @@ function addplane2model!(::IsPrimitive, model, plane, ipzmatricedict)
     # register distance variable:
     dz = @variable(model, base_name = string("d_z_", getfeaturename(plane)))
     
-    pzn = getpartzeroname(plane)
+    pzn = partzeroname(plane)
     v_machined = getmachinedfeaturepoint(plane)
     v_rough = getroughfeaturepoint(plane)
     # equation (4)
@@ -86,7 +86,7 @@ function addplane2model!(::IsFreeForm, model, plane, ipzmatricedict)
     # register distance variable:
     dz = @variable(model, [qiter], base_name = string("d_z_", getfeaturename(plane)))
 
-    pzn = getpartzeroname(plane)
+    pzn = partzeroname(plane)
     v_machined = getmachinedfeaturepoint(plane)
     # equation (4)
     for (i, q) in enumerate(qs)
@@ -110,8 +110,8 @@ function addtolerances2model!(model, mop::MultiOperationProblem, pzmatricedict)
         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)
+        pzn1 = partzeroname(f1)
+        pzn2 = partzeroname(f2)
 
         # equation (2)
         v1 = @expression(model, t.ismachined1 ? pzmatricedict[pzn1]*HV(getmachinedfeaturepoint(f1)) : getroughfeaturepoint(f1))
diff --git a/src/partzeros.jl b/src/partzeros.jl
index a5c86d2..60d6c59 100644
--- a/src/partzeros.jl
+++ b/src/partzeros.jl
@@ -28,7 +28,7 @@ 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)

From 54fe9a602ce8948060470edc53c12228228e063d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tam=C3=A1s=20Cserteg?= <csertegt@gmail.com>
Date: Fri, 19 Jan 2024 11:58:32 +0100
Subject: [PATCH 21/27] rework geometry storage

---
 src/BlankLocalizationCore.jl |  34 +++--
 src/geometries.jl            | 270 ++++++++++++-----------------------
 test/geometries.jl           |  55 +++++--
 test/runtests.jl             |   2 +-
 4 files changed, 150 insertions(+), 211 deletions(-)

diff --git a/src/BlankLocalizationCore.jl b/src/BlankLocalizationCore.jl
index 4e26241..de6fe1b 100644
--- a/src/BlankLocalizationCore.jl
+++ b/src/BlankLocalizationCore.jl
@@ -4,8 +4,9 @@ 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
 using Logging: @warn
 using Statistics: mean
 using LinearAlgebra: norm, dot, cross, inv, normalize, normalize!
@@ -14,22 +15,19 @@ using Printf: @sprintf
 export  PartZero,
         printpartzeropositions
 
-export  AbstractHoleGeometry,
-        AbstractPlaneGeometry,
-        GeometryStyle,
+export  GeometryStyle,
+        PLANELIKE,
+        HOLELIKE,
+        isplanelike,
+        isholelike,
         IsPrimitive,
         IsFreeForm,
-        surfacepoints,
-        filteredsurfacepoints,
         featurepoint,
+        surfacepoints,
         featureradius,
-        visualizationgeometry,
-        SimpleHole,
-        SimplePlane,
-        MeshHole,
-        MeshPlane,
-        HoleLocalizationFeature,
-        PlaneLocalizationFeature
+        RoughFeature,
+        MachinedFeature,
+        LocalizationFeature
         
 export  PlanePlaneDistance,
         PlaneAxisDistance,
@@ -69,10 +67,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("visualization.jl")
+#include("tolerances.jl")
+#include("optimizationproblem.jl")
+#include("optimization.jl")
+#include("resultevaluation.jl")
+#include("visualization.jl")
 
 end
diff --git a/src/geometries.jl b/src/geometries.jl
index 6fa7b97..c215032 100644
--- a/src/geometries.jl
+++ b/src/geometries.jl
@@ -1,12 +1,3 @@
-"""Supertype of 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
-
 """
     GeometryStyle
 
@@ -22,223 +13,142 @@ 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()
+"""Enum to store if a feature is hole or plane like."""
+@enum GeometryType PLANELIKE HOLELIKE
 
-Return the points of the surface of an [`IsFreeForm`](@ref) geometry.
-Definition signature should look like: `surfacepoints(::IsFreeForm, x)`.
-"""
-function surfacepoints end
+"""Supertype for feature types."""
+abstract type AbstractFeature{G} end
 
-"""
-    filteredsurfacepoints()
+GeometryStyle(f::AbstractFeature) = GeometryStyle(f.geometry)
+featurename(f::AbstractFeature) = f.name
+geometrytype(f::AbstractFeature) = f.gtype
+geometry(f::AbstractFeature) = f.geometry
 
-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
+isplanelike(st::AbstractFeature) = geometrytype(st) == PLANELIKE
+isholelike(st::AbstractFeature) = geometrytype(st) == HOLELIKE
 
 """
-    featureradius()
+    featurepoint(f::AbstractFeature)
 
-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`.
+Return the feature point of an [`IsPrimitive`](@ref) geometry.
 """
-function featureradius end
-
-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
-
-
-surfacepoints(x::T) where {T} = surfacepoints(GeometryStyle(T), x)
-function surfacepoints(::IsPrimitive, x)
-    error("Function `surfacepoints` is not defined for `IsPrimitive`` features")
-end
-
-
-filteredsurfacepoints(x::T) where {T} = filteredsurfacepoints(GeometryStyle(T), x)
-function filteredsurfacepoints(::IsPrimitive, x)
-    error("Function `filteredsurfacepoints` is not defined for `IsPrimitive`` features")
-end
-
-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")
+featurepoint(f::AbstractFeature) = featurepoint(GeometryStyle(f), f)
+featurepoint(::IsPrimitive, f) = featurepoint(geometry(f))
+function featurepoint(::IsFreeForm, f)
+    error("Function `featurepoint` is not defined for `IsFreeForm` features.")
 end
 
 """
-SimpleHole <: AbstractHoleGeometry
+    surfacepoints(f::AbstractFeature)
 
-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.
+Return the points of the surface of an [`IsFreeForm`](@ref) geometry.
 """
-struct SimpleHole <: AbstractHoleGeometry
-    p::Vector{Float64}
-    r::Float64
-end
+surfacepoints(f::AbstractFeature) = surfacepoints(GeometryStyle(f), f)
 
-GeometryStyle(::Type{SimpleHole}) = IsPrimitive()
+surfacepoints(::IsFreeForm, f) = surfacepoints(geometry(f))
 
-function visualizationgeometry(hole::SimpleHole)
-    # p1: feature point
-    p1 = Point3(hole.p)
-    return Disk(Plane(p1, Vec3(0,0,1)), hole.r)
+function surfacepoints(::IsPrimitive, f)
+    error("Function `surfacepoints` is not defined for `IsPrimitive` features.")
 end
 
 """
-SimplePlane <: AbstractPlaneGeometry
+    featureradius(f::AbstractFeature)
 
-A simple plane structure with one point and a normal vector.
+Return the radius of a [`IsPrimitive`](@ref) geometry, that is `HOLELIKE`.
 """
-struct SimplePlane <: AbstractPlaneGeometry
-    p::Vector{Float64}
-    n::Vector{Float64}
-end
+featureradius(f::AbstractFeature) = featureradius(GeometryStyle(f), f)
 
-GeometryStyle(::Type{SimplePlane}) = IsPrimitive()
-
-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)]))
+function featureradius(::IsPrimitive, f)
+    if isholelike(f)
+        return featureradius(geometry(f))
+    else
+        error("Function `featureradius` is only defined for `HOLELIKE` features.")
+    end
 end
 
-function visualizationgeometry(plane::SimplePlane)
-    o = plane.p
-    v1 = randnormal(plane.n)
-    v2 = cross(v1, plane.n)
-    return rectangleforplane(o, v1, v2, 20)
-end
-
-"""
-MeshHole <: AbstractHoleGeometry
-
-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}}
+function featureradius(::IsFreeForm, f)
+    error("Function `featureradius` is not defined for `IsFreeForm` features.")
 end
 
-GeometryStyle(::Type{MeshHole}) = IsFreeForm()
+## default functions for Cylinders, Planes and Meshes
+featurepoint(g::Meshes.Plane) = g(0,0)
 
-function surfacepoints(::IsFreeForm, x::MeshHole)
-    points = vertices(x.surface)
-    verts = [x.coords for x in points]
-    return verts
-end
+featurepoint(g::Meshes.Cylinder) = featurepoint(top(g))
 
-function visualizationgeometry(meshhole::MeshHole)
-    return meshhole.surface
-end
+surfacepoints(g::Meshes.Mesh) = vertices(g)
 
-function filteredsurfacepoints(::IsFreeForm, x::MeshHole)
-    return x.convexhull
-end
+featureradius(g::Meshes.Cylinder) = radius(g)
 
-"""
-MeshPlane <: AbstractPlaneGeometry
-
-A simple mesh plane geometry, that contains the mesh of a planar face.
-"""
-struct MeshPlane <: AbstractPlaneGeometry
-    surface::SimpleMesh
+struct RoughFeature{G} <: AbstractFeature{G}
+    name::String
+    gtype::GeometryType
+    geometry::G
 end
 
-GeometryStyle(::Type{MeshPlane}) = IsFreeForm()
+GeometryStyle(f::Meshes.Primitive) = IsPrimitive()
+GeometryStyle(f::Meshes.Domain) = IsFreeForm()
+#GeometryStyle(f) = IsFreeForm()
 
-function surfacepoints(::IsFreeForm, x::MeshPlane)
-    points = vertices(x.surface)
-    verts = [x.coords for x in points]
-    return verts
+#=
+for ST = subtypes(Meshes.Primitive)
+    eval(quote
+        GeometryStyle(::Type{$ST}) = IsPrimitive()    
+    end)
 end
+GeometryStyle(::Type) = IsFreeForm()
+=#
 
-function visualizationgeometry(meshplane::MeshPlane)
-    return meshplane.surface
+struct MachinedFeature{G<:Meshes.Primitive} <: AbstractFeature{G}
+    name::String
+    gtype::GeometryType
+    geometry::G
+    partzero::PartZero
 end
 
-function filteredsurfacepoints(::IsFreeForm, x::MeshPlane)
-    bbox = boundingbox(x.surface)
-    return [bbox.min.coords, bbox.max.coords]
-end
+partzero(f::MachinedFeature) = f.partzero
+partzeroname(f::MachinedFeature) = partzeroname(f.partzero)
 
 
 """
-    LocalizationFeature{R,M}
+    LocalizationFeature
 
 A feature that is machined and allowance can be computed for it.
-It has a name, a [`PartZero`](@ref), and a rough and machined geometry
-([`AbstractLocalizationGeometry`](@ref)).
-The two geometries should be of same type (hole, plane, etc.), but only
-[`HoleLocalizationFeature`](@ref) and [`PlaneLocalizationFeature`](@ref) enforce this property.
+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.
 """
-struct LocalizationFeature{R,M}
+struct LocalizationFeature{G}
     name::String
-    partzero::PartZero
-    rough::R
-    machined::M
-end
-
-const HoleLocalizationFeature{R,M} = LocalizationFeature{R,M} where {R<:AbstractHoleGeometry, M<:AbstractHoleGeometry}
-
-HoleLocalizationFeature(n, p, r, m) = LocalizationFeature(n, p, r, m)
-
-const PlaneLocalizationFeature{R,M} = LocalizationFeature{R,M} where {R<:AbstractPlaneGeometry, M<:AbstractPlaneGeometry}
-
-PlaneLocalizationFeature(n, p, r, m) = LocalizationFeature(n, p, r, m)
-
-GeometryStyle(::Type{LocalizationFeature{R,M}}) where {R,M} = GeometryStyle(R)
+    roughfeature::RoughFeature{G}
+    machinedfeature::MachinedFeature
+    function LocalizationFeature(n, r::RoughFeature{T}, m) where {T}
+        if isplanelike(r) == isplanelike(m)
+            new{T}(n, r, m)
+        else
+            error("GeometryType of rough and machined does not match!")
+        end
+    end
+end
+
+featurename(lf::LocalizationFeature) = lf.name
+partzero(lf::LocalizationFeature) = partzero(lf.machinedfeature)
+partzeroname(lf::LocalizationFeature) = partzeroname(lf.partzero)
+isplanelike(lf::LocalizationFeature) = isplanelike(lf.roughfeature)
+isholelike(lf::LocalizationFeature) = isholelike(lf.roughfeature)
+GeometryStyle(lf::LocalizationFeature) = GeometryStyle(lf.roughfeature)
 
 function Base.show(io::IO, lf::LocalizationFeature)
-    print(io, typeof(lf), ": ", getfeaturename(lf))
+    print(io, typeof(lf), ": ", featurename(lf))
 end
 
-getfeaturename(f::LocalizationFeature) = f.name
-getpartzero(f::LocalizationFeature) = f.partzero
-getpartzeroname(f::LocalizationFeature) = getpartzeroname(f.partzero)
-
-getroughfeaturepoint(f::LocalizationFeature) = featurepoint(f.rough)
-getmachinedfeaturepoint(f::LocalizationFeature) = featurepoint(f.machined)
-getroughradius(f::LocalizationFeature) = featureradius(f.rough)
-getmachinedradius(f::LocalizationFeature) = featureradius(f.machined)
-
-getroughfilteredpoints(f::LocalizationFeature) = filteredsurfacepoints(f.rough)
+roughfeaturepoint(lf::LocalizationFeature) = featurepoint(lf.rough)
+machinedfeaturepoint(lf::LocalizationFeature) = featurepoint(lf.machined)
+roughradius(lf::LocalizationFeature) = featureradius(lf.rough)
+machinedradius(lf::LocalizationFeature) = featureradius(lf.machined)
 
-function getmachinedfeaturepointindatum(f::LocalizationFeature)
-    v = getmachinedfeaturepoint(f)
-    pz = getpartzero(f)
+function machinedfeaturepointindatum(f::LocalizationFeature)
+    v = machinedfeaturepoint(f)
+    pz = partzero(f)
     T = getpartzeroHM(pz)
     v_indatum = T*HV(v)
     return v_indatum[1:3]
@@ -250,7 +160,7 @@ end
 Transform a list of points with the part zero of `feature`.
 """
 function transformmachined2datum(feature, points)
-    pz = getpartzero(feature)
+    pz = partzero(feature)
     M = getpartzeroHM(pz)
     newpoints = (M*HV(p) for p in points)
     resultpoints = [p[1:3] for p in newpoints]
diff --git a/test/geometries.jl b/test/geometries.jl
index c48357b..e7abd1c 100644
--- a/test/geometries.jl
+++ b/test/geometries.jl
@@ -1,14 +1,45 @@
 @testset "IsPrimitive geometries" begin
-    sh = SimpleHole([0, 0, 0], 29)
-    sp = SimplePlane([0, 0, 0], [0, 0, 1])
-
-    @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)
+    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 e8f08a3..806733e 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -6,7 +6,7 @@ const BLC = BlankLocalizationCore
 
 #using Ipopt
 
-include("tolerances.jl")
+#include("tolerances.jl")
 include("partzeros.jl")
 include("geometries.jl")
 #include("testproblem.jl")

From b1ee3dbd4354380e4d1033ed510133acbb344f95 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tam=C3=A1s=20Cserteg?= <csertegt@gmail.com>
Date: Wed, 24 Jan 2024 15:37:17 +0100
Subject: [PATCH 22/27] rework feature and tolerance storage

---
 src/BlankLocalizationCore.jl |  29 +++---
 src/geometries.jl            | 171 ++++++++++++++++++++---------------
 src/optimizationproblem.jl   |  44 +++++----
 src/resultevaluation.jl      |  24 ++---
 src/tolerances.jl            | 102 +++++++++++++++++++++
 5 files changed, 255 insertions(+), 115 deletions(-)

diff --git a/src/BlankLocalizationCore.jl b/src/BlankLocalizationCore.jl
index de6fe1b..87406d3 100644
--- a/src/BlankLocalizationCore.jl
+++ b/src/BlankLocalizationCore.jl
@@ -15,26 +15,25 @@ using Printf: @sprintf
 export  PartZero,
         printpartzeropositions
 
-export  GeometryStyle,
-        PLANELIKE,
-        HOLELIKE,
-        isplanelike,
-        isholelike,
-        IsPrimitive,
-        IsFreeForm,
+export  RepresentationStyle,
+        FeatureStyle,
+        AbstractLocalizationGeometry,
+        SimpleHole,
+        MeshHole,
+        SimplePlane,
+        MeshPlane,
         featurepoint,
         surfacepoints,
+        filteredsurfacepoints,
         featureradius,
         RoughFeature,
         MachinedFeature,
         LocalizationFeature
         
-export  PlanePlaneDistance,
-        PlaneAxisDistance,
-        AxisAxisDistance,
-        AxisAxisConcentric,
-        LocalizationTolerance,
-        toleranceddistance
+export  PositionTolerance,
+        ConcentrictyTolerance,
+        ProjectedDimensionTolerance,
+        addtolerance2model!
         
 export  OptimizationResult,
         MultiOperationProblem,
@@ -67,8 +66,8 @@ HV(v) = vcat(v, 1)
 
 include("partzeros.jl")
 include("geometries.jl")
-#include("tolerances.jl")
-#include("optimizationproblem.jl")
+include("tolerances.jl")
+include("optimizationproblem.jl")
 #include("optimization.jl")
 #include("resultevaluation.jl")
 #include("visualization.jl")
diff --git a/src/geometries.jl b/src/geometries.jl
index c215032..804e4da 100644
--- a/src/geometries.jl
+++ b/src/geometries.jl
@@ -1,107 +1,129 @@
 """
-    GeometryStyle
+    RepresentationStyle
 
 Trait that describes the "style" of an [`AbstractLocalizationGeometry`](@ref).
-If it can be described by a primitive, then it is [`IsPrimitive`](@ref) and
-[`IsFreeForm`](@ref) otherwise.
+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
+struct Primitive <: RepresentationStyle end
 
 """Free form geometries are discrete representations, e.g. a mesh or a point cloud."""
-struct IsFreeForm <: GeometryStyle end
+struct FreeForm <: RepresentationStyle end
+
+nicestr(x::Primitive) = "primitive"
+nicestr(x::FreeForm) = "free form"
+
+abstract type FeatureStyle end
+struct Planar <: FeatureStyle end
+struct Cylindrical <: FeatureStyle end
+
+nicestr(x::Planar) = "planar"
+nicestr(x::Cylindrical) = "cylindrical"
+
+abstract type AbstractLocalizationGeometry end
+
+struct SimpleHole <: AbstractLocalizationGeometry
+    geom::Meshes.Disk
+end
+
+RepresentationStyle(g::SimpleHole) = Primitive()
+FeatureStyle(g::SimpleHole) = Cylindrical()
+
+featurepoint(g::SimpleHole) = plane(g.geom)(0,0)
+featureradius(g::SimpleHole) = radius(g.geom)
+
+struct MeshHole <: AbstractLocalizationGeometry
+    geom #::Meshes.Mesh, but that is abstract
+    chull # vector of points / vector of vectors - TBD
+end
+
+RepresentationStyle(g::MeshHole) = FreeForm()
+FeatureStyle(g::MeshHole) = Cylindrical()
+
+surfacepoints(g::MeshHole) = vertices(g.geom)
+filteredsurfacepoints(g::MeshHole) = g.chull
+
+struct SimplePlane <: AbstractLocalizationGeometry
+    geom::Meshes.Plane
+end
+
+RepresentationStyle(g::SimplePlane) = Primitive()
+FeatureStyle(g::SimplePlane) = Planar()
+
+featurepoint(g::SimplePlane) = g.geom(0,0)
+
+struct MeshPlane <: AbstractLocalizationGeometry
+    geom #::Meshes.Mesh, but that is abstract
+end
+
+RepresentationStyle(g::MeshPlane) = FreeForm()
+FeatureStyle(g::MeshPlane) = Planar()
+
+surfacepoints(g::MeshPlane) = vertices(g.geom)
+filteredsurfacepoints(g::MeshPlane) = vertices(g.geom)
+
+
 
-"""Enum to store if a feature is hole or plane like."""
-@enum GeometryType PLANELIKE HOLELIKE
 
 """Supertype for feature types."""
-abstract type AbstractFeature{G} end
+abstract type AbstractFeature end
 
-GeometryStyle(f::AbstractFeature) = GeometryStyle(f.geometry)
+RepresentationStyle(f::AbstractFeature) = RepresentationStyle(f.geometry)
+FeatureStyle(f::AbstractFeature) = FeatureStyle(f.geometry)
 featurename(f::AbstractFeature) = f.name
-geometrytype(f::AbstractFeature) = f.gtype
 geometry(f::AbstractFeature) = f.geometry
 
-isplanelike(st::AbstractFeature) = geometrytype(st) == PLANELIKE
-isholelike(st::AbstractFeature) = geometrytype(st) == HOLELIKE
+isplanar(st::AbstractFeature) = FeatureStyle(st) === Planar()
+iscylindrical(st::AbstractFeature) = FeatureStyle(st) === Cylindrical()
 
 """
     featurepoint(f::AbstractFeature)
 
-Return the feature point of an [`IsPrimitive`](@ref) geometry.
+Return the feature point of an [`Primitive`](@ref) geometry.
 """
-featurepoint(f::AbstractFeature) = featurepoint(GeometryStyle(f), f)
-featurepoint(::IsPrimitive, f) = featurepoint(geometry(f))
-function featurepoint(::IsFreeForm, f)
-    error("Function `featurepoint` is not defined for `IsFreeForm` features.")
+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
 
 """
     surfacepoints(f::AbstractFeature)
 
-Return the points of the surface of an [`IsFreeForm`](@ref) geometry.
+Return the points of the surface of an [`FreeForm`](@ref) geometry.
 """
-surfacepoints(f::AbstractFeature) = surfacepoints(GeometryStyle(f), f)
+surfacepoints(f::AbstractFeature) = surfacepoints(RepresentationStyle(f), f)
 
-surfacepoints(::IsFreeForm, f) = surfacepoints(geometry(f))
+surfacepoints(::FreeForm, f) = surfacepoints(geometry(f))
 
-function surfacepoints(::IsPrimitive, f)
-    error("Function `surfacepoints` is not defined for `IsPrimitive` features.")
+function surfacepoints(::Primitive, f)
+    error("Function `surfacepoints` is not defined for `Primitive` features.")
 end
 
 """
     featureradius(f::AbstractFeature)
 
-Return the radius of a [`IsPrimitive`](@ref) geometry, that is `HOLELIKE`.
+Return the radius of a [`Primitive`](@ref) geometry, that is `HOLELIKE`.
 """
-featureradius(f::AbstractFeature) = featureradius(GeometryStyle(f), f)
+featureradius(f::AbstractFeature) = featureradius(RepresentationStyle(f), FeatureStyle(f), f)
 
-function featureradius(::IsPrimitive, f)
-    if isholelike(f)
-        return featureradius(geometry(f))
-    else
-        error("Function `featureradius` is only defined for `HOLELIKE` features.")
-    end
-end
+featureradius(::Primitive, ::Cylindrical, f) = featureradius(geometry(f))
 
-function featureradius(::IsFreeForm, f)
-    error("Function `featureradius` is not defined for `IsFreeForm` features.")
+function featureradius(t1, t2, f)
+    error("Function `featureradius` is only defined for ::Primitive and ::Cylindrical features. Got $t1 and $t2")
 end
 
-## default functions for Cylinders, Planes and Meshes
-featurepoint(g::Meshes.Plane) = g(0,0)
-
-featurepoint(g::Meshes.Cylinder) = featurepoint(top(g))
-
-surfacepoints(g::Meshes.Mesh) = vertices(g)
-
-featureradius(g::Meshes.Cylinder) = radius(g)
-
-struct RoughFeature{G} <: AbstractFeature{G}
+struct RoughFeature <: AbstractFeature
     name::String
-    gtype::GeometryType
-    geometry::G
+    geometry # <: AbstractLocalizationGeometry
 end
 
-GeometryStyle(f::Meshes.Primitive) = IsPrimitive()
-GeometryStyle(f::Meshes.Domain) = IsFreeForm()
-#GeometryStyle(f) = IsFreeForm()
-
-#=
-for ST = subtypes(Meshes.Primitive)
-    eval(quote
-        GeometryStyle(::Type{$ST}) = IsPrimitive()    
-    end)
-end
-GeometryStyle(::Type) = IsFreeForm()
-=#
-
-struct MachinedFeature{G<:Meshes.Primitive} <: AbstractFeature{G}
+struct MachinedFeature <: AbstractFeature
     name::String
-    gtype::GeometryType
-    geometry::G
+    geometry # <: AbstractLocalizationGeometry
     partzero::PartZero
 end
 
@@ -117,15 +139,15 @@ 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.
 """
-struct LocalizationFeature{G}
+struct LocalizationFeature
     name::String
-    roughfeature::RoughFeature{G}
+    roughfeature::RoughFeature
     machinedfeature::MachinedFeature
-    function LocalizationFeature(n, r::RoughFeature{T}, m) where {T}
-        if isplanelike(r) == isplanelike(m)
-            new{T}(n, r, m)
+    function LocalizationFeature(n, r, m)
+        if FeatureStyle(r) === FeatureStyle(m)
+            new(n, r, m)
         else
-            error("GeometryType of rough and machined does not match!")
+            error("FeatureStyle of rough and machined features does not match!")
         end
     end
 end
@@ -133,12 +155,14 @@ end
 featurename(lf::LocalizationFeature) = lf.name
 partzero(lf::LocalizationFeature) = partzero(lf.machinedfeature)
 partzeroname(lf::LocalizationFeature) = partzeroname(lf.partzero)
-isplanelike(lf::LocalizationFeature) = isplanelike(lf.roughfeature)
-isholelike(lf::LocalizationFeature) = isholelike(lf.roughfeature)
-GeometryStyle(lf::LocalizationFeature) = GeometryStyle(lf.roughfeature)
+isplanar(lf::LocalizationFeature) = isplanar(lf.roughfeature)
+iscylindrical(lf::LocalizationFeature) = iscylindrical(lf.roughfeature)
+RepresentationStyle(lf::LocalizationFeature) = RepresentationStyle(lf.roughfeature)
+FeatureStyle(lf::LocalizationFeature) = FeatureStyle(lf.roughfeature)
 
 function Base.show(io::IO, lf::LocalizationFeature)
-    print(io, typeof(lf), ": ", featurename(lf))
+    print(io, "LocalizationFeature: ", featurename(lf),
+    " ", nicestr(RepresentationStyle(lf)), ", ", nicestr(FeatureStyle(lf)))
 end
 
 roughfeaturepoint(lf::LocalizationFeature) = featurepoint(lf.rough)
@@ -146,6 +170,7 @@ machinedfeaturepoint(lf::LocalizationFeature) = featurepoint(lf.machined)
 roughradius(lf::LocalizationFeature) = featureradius(lf.rough)
 machinedradius(lf::LocalizationFeature) = featureradius(lf.machined)
 
+#=
 function machinedfeaturepointindatum(f::LocalizationFeature)
     v = machinedfeaturepoint(f)
     pz = partzero(f)
@@ -166,3 +191,5 @@ function transformmachined2datum(feature, points)
     resultpoints = [p[1:3] for p in newpoints]
     return resultpoints
 end
+=#
+
diff --git a/src/optimizationproblem.jl b/src/optimizationproblem.jl
index fd51549..72f20e0 100644
--- a/src/optimizationproblem.jl
+++ b/src/optimizationproblem.jl
@@ -36,7 +36,7 @@ The features in tolerances ([`LocalizationTolerance`](@ref)) can be `nothing`, i
 mutable struct MultiOperationProblem
     partzeros::Vector{PartZero}
     features::Vector{LocalizationFeature} # features that have rough and machined -> allowanced
-    tolerances::Vector{LocalizationTolerance}
+    tolerances # Vector{AbstractTolerance}
     parameters::Dict{String,Real}
     opresult::OptimizationResult
 end
@@ -46,26 +46,26 @@ function MultiOperationProblem(partzeros, features, tolerances, parameters)
 end
 
 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
-    featuretypes = GeometryStyle.(typeof.(x.rough for x in mop.features))
+    # 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 === IsFreeForm() && return :HybridProblem
+        ft === FreeForm() && return :HybridProblem
     end
     return :PrimitiveProblem
 end
 
-function collectholefeatures(mop::MultiOperationProblem)
-    filter(x->x isa HoleLocalizationFeature, mop.features)
+function cylindricalfeatures(mop::MultiOperationProblem)
+    filter(x->FeatureStyle(x) === Cylindrical(), mop.features)
 end
 
-function collectplanefeatures(mop::MultiOperationProblem)
-    filter(x->x isa PlaneLocalizationFeature, mop.features)
+function planarfeatures(mop::MultiOperationProblem)
+    filter(x->FeatureStyle(x) === Planar(), mop.features)
 end
 
 function Base.show(io::IO, mop::MultiOperationProblem)
-    nh = size(collectholefeatures(mop), 1)
-    np = size(collectplanefeatures(mop), 1)
+    nh = size(cylindricalfeatures(mop), 1)
+    np = size(planarfeatures(mop), 1)
     npz = size(mop.partzeros, 1)
     nts = size(mop.tolerances, 1)
     sn = string(problemtype(mop))
@@ -101,22 +101,34 @@ Tell if `mop`'s solution is in an optimal state, either: `OPTIMAL` or `LOCALLY_S
 """
 isoptimum(mop::MultiOperationProblem) = isoptimum(mop.opresult)
 
+
 """
-    getfeaturebyname(mop::MultiOperationProblem, featurename)
+    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 `featurename`.
+Return `nothing`, if no feature is found with `featuresname`.
 """
-function getfeaturebyname(mop::MultiOperationProblem, featurename)
-    for f in mop.features
-        if getfeaturename(f) == featurename
+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)
 
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
index 543ab5f..1b31393 100644
--- a/src/tolerances.jl
+++ b/src/tolerances.jl
@@ -1,3 +1,103 @@
+"""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
+    projection
+end
+
+function Base.show(io::IO, t::ProjectedDimensionTolerance)
+    print(io, "ProjDimTol<", featurename(t.datumfeature), "-", featurename(t.feature),
+    " : ", t.minval, "-", t.maxval, ">",)
+end
+
+"""
+    projection(t::ProjectedDimensionTolerance)
+
+Return the projection function of `t`.
+"""
+projection(t::ProjectedDimensionTolerance) = t.projection
+
+
+#=
+
 """Supertype of tolerance types."""
 abstract type LocalizationToleranceType end
 
@@ -246,3 +346,5 @@ function evaluatetolerance(t::LocalizationTolerance)
     rel_d = 2*abs_d/(t.uppervalue-t.lowervalue)*100
     return (d, rel_d)
 end
+
+=#

From e445ef6e76fefd9cacebe26ffd9271899f985d0f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tam=C3=A1s=20Cserteg?= <csertegt@gmail.com>
Date: Thu, 25 Jan 2024 10:23:18 +0100
Subject: [PATCH 23/27] update Meshes to 0.39

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 9abb3f0..1d13794 100644
--- a/Project.toml
+++ b/Project.toml
@@ -17,7 +17,7 @@ Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc"
 [compat]
 DataFrames = "1"
 JuMP = "1.4"
-Meshes = "0.35,0.36"
+Meshes = "0.39"
 PrettyTables = "2"
 Rotations = "1.5"
 julia = "1.9"

From f938e0ede21c26e5eac3b67cb95676a25b79fdcd Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tam=C3=A1s=20Cserteg?= <csertegt@gmail.com>
Date: Thu, 25 Jan 2024 11:10:00 +0100
Subject: [PATCH 24/27] update jump model construction

---
 src/BlankLocalizationCore.jl |   2 +-
 src/geometries.jl            |  16 ++--
 src/optimization.jl          | 156 ++++++++++++++++++-----------------
 src/optimizationproblem.jl   |   5 --
 src/tolerances.jl            |   8 +-
 5 files changed, 95 insertions(+), 92 deletions(-)

diff --git a/src/BlankLocalizationCore.jl b/src/BlankLocalizationCore.jl
index 87406d3..90e5257 100644
--- a/src/BlankLocalizationCore.jl
+++ b/src/BlankLocalizationCore.jl
@@ -68,7 +68,7 @@ include("partzeros.jl")
 include("geometries.jl")
 include("tolerances.jl")
 include("optimizationproblem.jl")
-#include("optimization.jl")
+include("optimization.jl")
 #include("resultevaluation.jl")
 #include("visualization.jl")
 
diff --git a/src/geometries.jl b/src/geometries.jl
index 804e4da..0a6a22f 100644
--- a/src/geometries.jl
+++ b/src/geometries.jl
@@ -32,7 +32,7 @@ end
 RepresentationStyle(g::SimpleHole) = Primitive()
 FeatureStyle(g::SimpleHole) = Cylindrical()
 
-featurepoint(g::SimpleHole) = plane(g.geom)(0,0)
+featurepoint(g::SimpleHole) = Meshes.plane(g.geom)(0,0)
 featureradius(g::SimpleHole) = radius(g.geom)
 
 struct MeshHole <: AbstractLocalizationGeometry
@@ -106,7 +106,7 @@ end
 """
     featureradius(f::AbstractFeature)
 
-Return the radius of a [`Primitive`](@ref) geometry, that is `HOLELIKE`.
+Return the radius of a [`Primitive`](@ref) geometry, that is `Cylindrical`.
 """
 featureradius(f::AbstractFeature) = featureradius(RepresentationStyle(f), FeatureStyle(f), f)
 
@@ -154,7 +154,7 @@ end
 
 featurename(lf::LocalizationFeature) = lf.name
 partzero(lf::LocalizationFeature) = partzero(lf.machinedfeature)
-partzeroname(lf::LocalizationFeature) = partzeroname(lf.partzero)
+partzeroname(lf::LocalizationFeature) = partzeroname(partzero(lf))
 isplanar(lf::LocalizationFeature) = isplanar(lf.roughfeature)
 iscylindrical(lf::LocalizationFeature) = iscylindrical(lf.roughfeature)
 RepresentationStyle(lf::LocalizationFeature) = RepresentationStyle(lf.roughfeature)
@@ -165,10 +165,12 @@ function Base.show(io::IO, lf::LocalizationFeature)
     " ", nicestr(RepresentationStyle(lf)), ", ", nicestr(FeatureStyle(lf)))
 end
 
-roughfeaturepoint(lf::LocalizationFeature) = featurepoint(lf.rough)
-machinedfeaturepoint(lf::LocalizationFeature) = featurepoint(lf.machined)
-roughradius(lf::LocalizationFeature) = featureradius(lf.rough)
-machinedradius(lf::LocalizationFeature) = featureradius(lf.machined)
+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)
 
 #=
 function machinedfeaturepointindatum(f::LocalizationFeature)
diff --git a/src/optimization.jl b/src/optimization.jl
index 1ca9811..e3533f2 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)
-
-    pzn = partzeroname(hole)
-    v_machined = getmachinedfeaturepoint(hole)
-    v_rough = getroughfeaturepoint(hole)
-    r_machined = getmachinedradius(hole)
-    r_rough = getroughradius(hole)
+    dxy = @variable(model, base_name = string("d_xy_", featurename(feature)), lower_bound = 0.0)
+
+    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,74 @@ 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 = partzeroname(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
+    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_", getfeaturename(feature)), lower_bound = 0.0)
     
-    pzn = partzeroname(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, qmeshes) 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 = getroughfilteredpoints(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 = partzeroname(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)
+    for (i, qmeshes) 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 +107,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 = partzeroname(f1)
-        pzn2 = partzeroname(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 +170,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 +184,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
index 72f20e0..fa86323 100644
--- a/src/optimizationproblem.jl
+++ b/src/optimizationproblem.jl
@@ -27,11 +27,6 @@ end
     MultiOperationProblem
 
 A type storing all geometry data and parameters to generate a JuMP model.
-
-The elements of `features` can only be `HoleLocalizationFeature`s and `PlaneLocalizationFeature`s.
-Otherwise the JuMP model building fails.
-
-The features in tolerances ([`LocalizationTolerance`](@ref)) can be `nothing`, if needed.
 """
 mutable struct MultiOperationProblem
     partzeros::Vector{PartZero}
diff --git a/src/tolerances.jl b/src/tolerances.jl
index 1b31393..d7fe4a9 100644
--- a/src/tolerances.jl
+++ b/src/tolerances.jl
@@ -80,7 +80,7 @@ struct ProjectedDimensionTolerance <: AbstractTolerance
     nominalval
     minval
     maxval
-    projection
+    projectionaxis
 end
 
 function Base.show(io::IO, t::ProjectedDimensionTolerance)
@@ -89,11 +89,11 @@ function Base.show(io::IO, t::ProjectedDimensionTolerance)
 end
 
 """
-    projection(t::ProjectedDimensionTolerance)
+    projectionaxis(t::ProjectedDimensionTolerance)
 
-Return the projection function of `t`.
+Return the projection axis of `t`.
 """
-projection(t::ProjectedDimensionTolerance) = t.projection
+projectionaxis(t::ProjectedDimensionTolerance) = t.projectionaxis
 
 
 #=

From 5f60c6b9db24e84b2a37c3628a692b8f1c8d8c0d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tam=C3=A1s=20Cserteg?= <csertegt@gmail.com>
Date: Sun, 28 Jan 2024 20:48:29 +0100
Subject: [PATCH 25/27] uptade Meshes compat to 0.40

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 1d13794..67c8038 100644
--- a/Project.toml
+++ b/Project.toml
@@ -17,7 +17,7 @@ Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc"
 [compat]
 DataFrames = "1"
 JuMP = "1.4"
-Meshes = "0.39"
+Meshes = "0.39,0.40"
 PrettyTables = "2"
 Rotations = "1.5"
 julia = "1.9"

From 0a7c7724edf60024d5c7aeb9f954589049bbf745 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tam=C3=A1s=20Cserteg?= <csertegt@gmail.com>
Date: Sun, 28 Jan 2024 20:49:33 +0100
Subject: [PATCH 26/27] update visualization code for current state

---
 src/BlankLocalizationCore.jl |  4 ++--
 src/geometries.jl            | 22 ++++++++++++++++++++++
 src/visualization.jl         | 18 ++++++++----------
 3 files changed, 32 insertions(+), 12 deletions(-)

diff --git a/src/BlankLocalizationCore.jl b/src/BlankLocalizationCore.jl
index 90e5257..8d724a9 100644
--- a/src/BlankLocalizationCore.jl
+++ b/src/BlankLocalizationCore.jl
@@ -6,7 +6,7 @@ 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, top, radius
+    Rotate, Translate, Disk, top, radius, normal
 using Logging: @warn
 using Statistics: mean
 using LinearAlgebra: norm, dot, cross, inv, normalize, normalize!
@@ -70,6 +70,6 @@ include("tolerances.jl")
 include("optimizationproblem.jl")
 include("optimization.jl")
 #include("resultevaluation.jl")
-#include("visualization.jl")
+include("visualization.jl")
 
 end
diff --git a/src/geometries.jl b/src/geometries.jl
index 0a6a22f..2d87a3b 100644
--- a/src/geometries.jl
+++ b/src/geometries.jl
@@ -25,6 +25,8 @@ nicestr(x::Cylindrical) = "cylindrical"
 
 abstract type AbstractLocalizationGeometry end
 
+visualizationgeometry(g::AbstractLocalizationGeometry) = g.geom
+
 struct SimpleHole <: AbstractLocalizationGeometry
     geom::Meshes.Disk
 end
@@ -55,6 +57,26 @@ FeatureStyle(g::SimplePlane) = Planar()
 
 featurepoint(g::SimplePlane) = g.geom(0,0)
 
+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
+    return SimpleMesh([c1,c2,c3,c4], connect.([(1,2,3),(3,4,1)]))
+end
+
+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
+
+visualizationgeometry(plane::SimplePlane) = rectangleforplane(plane.geom)
+
 struct MeshPlane <: AbstractLocalizationGeometry
     geom #::Meshes.Mesh, but that is abstract
 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

From cba9713f4feb0a30252fb7305c33859a0cde493b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tam=C3=A1s=20Cserteg?= <csertegt@gmail.com>
Date: Sun, 28 Jan 2024 21:42:45 +0100
Subject: [PATCH 27/27] update mesh optimization code

---
 src/BlankLocalizationCore.jl |  2 +-
 src/geometries.jl            | 20 ++++++++++++++++++--
 src/optimization.jl          | 13 +++++++------
 3 files changed, 26 insertions(+), 9 deletions(-)

diff --git a/src/BlankLocalizationCore.jl b/src/BlankLocalizationCore.jl
index 8d724a9..50c8720 100644
--- a/src/BlankLocalizationCore.jl
+++ b/src/BlankLocalizationCore.jl
@@ -6,7 +6,7 @@ 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, top, radius, normal
+    Rotate, Translate, Disk, top, radius, normal, boundingbox
 using Logging: @warn
 using Statistics: mean
 using LinearAlgebra: norm, dot, cross, inv, normalize, normalize!
diff --git a/src/geometries.jl b/src/geometries.jl
index 2d87a3b..1500c8d 100644
--- a/src/geometries.jl
+++ b/src/geometries.jl
@@ -85,8 +85,11 @@ RepresentationStyle(g::MeshPlane) = FreeForm()
 FeatureStyle(g::MeshPlane) = Planar()
 
 surfacepoints(g::MeshPlane) = vertices(g.geom)
-filteredsurfacepoints(g::MeshPlane) = vertices(g.geom)
 
+function filteredsurfacepoints(g::MeshPlane)
+    bb = boundingbox(vertices(g.geom))
+    return [[bb.min.coords.coords...], [bb.max.coords.coords...]]
+end
 
 
 
@@ -115,7 +118,7 @@ end
 """
     surfacepoints(f::AbstractFeature)
 
-Return the points of the surface of an [`FreeForm`](@ref) geometry.
+Return the points of the surface of a [`FreeForm`](@ref) geometry.
 """
 surfacepoints(f::AbstractFeature) = surfacepoints(RepresentationStyle(f), f)
 
@@ -125,6 +128,19 @@ function surfacepoints(::Primitive, f)
     error("Function `surfacepoints` is not defined for `Primitive` features.")
 end
 
+"""
+    filteredsurfacepoints(f::AbstractFeature)
+
+Return the filtered points of the surface of a [`FreeForm`](@ref) geometry.
+"""
+filteredsurfacepoints(f::AbstractFeature) = filteredsurfacepoints(RepresentationStyle(f), f)
+
+filteredsurfacepoints(::FreeForm, f) = filteredsurfacepoints(geometry(f))
+
+function filteredsurfacepoints(::Primitive, f)
+    error("Function `surfacepoints` is not defined for `Primitive` features.")
+end
+
 """
     featureradius(f::AbstractFeature)
 
diff --git a/src/optimization.jl b/src/optimization.jl
index e3533f2..bd8d92f 100644
--- a/src/optimization.jl
+++ b/src/optimization.jl
@@ -55,19 +55,20 @@ function addallowancedfeature2model!(::FreeForm, ::Cylindrical, model, feature,
     minAllowance = model[:minAllowance]
     
     # filtered surface points of a free form surface
+    @show feature
     qs = roughfilteredsurfacepoints(feature)
     qiter = 1:length(qs)
     
     # register distance variable:
-    dxy = @variable(model, [qiter], base_name = string("d_xy_", getfeaturename(feature)), lower_bound = 0.0)
+    dxy = @variable(model, [qiter], base_name = string("d_xy_", featurename(feature)), lower_bound = 0.0)
     
     pzn = partzeroname(feature)
     v_machined_ = machinedfeaturepoint(feature)
     v_machined = [v_machined_.coords.coords...]
     r_machined = machinedradius(feature)
     # equation (4)
-    for (i, qmeshes) in enumerate(qs)
-        q = [qmeshes.coords.coords...]
+    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])
@@ -83,7 +84,7 @@ function addallowancedfeature2model!(::FreeForm, ::Planar, model, feature, ipzma
     maxPlaneZAllowance = model[:maxPlaneZAllowance]
 
     # filtered surface points of a free form surface
-    qs = getroughfilteredpoints(feature)
+    qs = roughfilteredsurfacepoints(feature)
     qiter = 1:length(qs)
 
     # register distance variable:
@@ -93,8 +94,8 @@ function addallowancedfeature2model!(::FreeForm, ::Planar, model, feature, ipzma
     v_machined_ = machinedfeaturepoint(feature)
     v_machined = [v_machined_.coords.coords...]
     # equation (4)
-    for (i, qmeshes) in enumerate(qs)
-        q = [qmeshes.coords.coords...]
+    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])