Merge pull request #51 from ModiaSim/an_mprWithDouble64

An mpr with double64

Author: GerhardHippmann

Project.toml

@@ -6,6 +6,7 @@ version = "0.5.1-dev"
 [deps]
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -21,6 +22,7 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 [compat]
 Colors = "0.12, 0.11, 0.10"
 DataFrames = "0.22, 0.21, 0.20, 0.19"
+DoubleFloats = "1.1"
 JSON = "0.21"
 ModiaLang = "0.8.1"
 OrderedCollections = "1.4, 1.3, 1.2, 1.1"

src/Basics/_module.jl

@@ -17,7 +17,7 @@ import Modia3D
 
 export trailingPartOfName
 
-export neps, sign_eps, radToDeg
+export nepsType, neps, sign_eps, radToDeg
 export getAndCheckFullLibraryPath, getEnvironmentVariable
 
 export normalizeVector, BoundingBox

src/Basics/constantsAndFunctions.jl

@@ -7,19 +7,21 @@
 
 
 # Epsilon and sign
-const neps = sqrt( eps(Modia3D.MPRFloatType) )
+const neps = sqrt( eps() )
 
-function sign_eps(value::T; ) where {T}
-    seps::T = 100.0*neps
+nepsType(::Type{T}) where {T} = sqrt( eps(T) ) # mpr and bounding box calculation use this
+
+function sign_eps(value::T) where {T}
+    seps::T = 100.0*nepsType(T)
     return value > seps ? T(1.0) : (value < -seps ? T(-1.0) : T(0.0))
 end
 
 function normalizeVector(n::SVector{3,T}) where {T}
     nabs = norm(n)
-    if nabs <= neps
-        println("neps ", neps)
+    if nabs <= nepsType(T)
+        println("neps ", nepsType(T))
         println("nabs ", nabs)
-        @assert(nabs > neps) # && norm(vec) > eps()
+        @assert(nabs > nepsType(T)) # && norm(vec) > eps()
         # return nothing
     end
     return n/nabs

src/Frames/interpolation.jl

@@ -4,7 +4,6 @@
 # This file is part of module
 #   Modia3D.Frames (Modia3D/Frames/_module.jl)
 #
-const seps = sqrt(eps())
 
 
 """
@@ -76,7 +75,7 @@ struct Path
 
     function Path(r::AbstractVector, q::AbstractVector=Quaternion[];
                   v::AbstractVector=ones(length(r)),
-                  seps=sqrt(eps()) )
+                  seps=Modia3D.neps )
         nframes = size(r, 1)
         @assert(seps > 0.0)
         @assert(nframes > 1)

src/Modia3D.jl

@@ -96,6 +96,7 @@ convertAndStripUnit(TargetType, requiredUnit, value) =
             convert(TargetType, ustrip.( uconvert.(requiredUnit, value))) : convert(TargetType, value)
 
 # MPRFloatType is used to change betweeen Double64 and Float64 for mpr calculations
+using DoubleFloats
 const MPRFloatType = Float64
 
 # Include sub-modules

src/Shapes/boundingBoxes.jl

@@ -78,21 +78,21 @@ end
         halfDiameter = T(0.5*shape.diameter)
         if shape.axis == 1
             enorm = norm(SVector(e_abs[2], e_abs[3]))
-            if enorm <= Modia3D.neps
+            if enorm <= Modia3D.nepsType(T)
                 return Basics.sign_eps(e_abs[1])*SVector(halfLength, 0.0, 0.0)
             else
                 return Basics.sign_eps(e_abs[1])*SVector(halfLength, 0.0, 0.0) + SVector(0.0, (halfDiameter)*e_abs[2], (halfDiameter)*e_abs[3])/enorm
             end
         elseif shape.axis == 2
             enorm = norm(SVector(e_abs[3], e_abs[1]))
-            if enorm <= Modia3D.neps
+            if enorm <= Modia3D.nepsType(T)
                 return Basics.sign_eps(e_abs[2])*SVector(0.0, halfLength, 0.0)
             else
                 return Basics.sign_eps(e_abs[2])*SVector(0.0, halfLength, 0.0) + SVector((halfDiameter)*e_abs[1], 0.0, (halfDiameter)*e_abs[3])/enorm
             end
         else
             enorm = norm(SVector(e_abs[1], e_abs[2]))
-            if enorm <= Modia3D.neps
+            if enorm <= Modia3D.nepsType(T)
                 return Basics.sign_eps(e_abs[3])*SVector(0.0, 0.0, halfLength)
             else
                 return Basics.sign_eps(e_abs[3])*SVector(0.0, 0.0, halfLength) + SVector((halfDiameter)*e_abs[1], (halfDiameter)*e_abs[2], 0.0)/enorm
@@ -137,15 +137,15 @@ end
                     return SVector(shapeLength, 0.0, 0.0)  # apex is support point
                 else  # frustum of a cone
                     enorm = norm(SVector(e_abs[2], e_abs[3]))
-                    if enorm > Modia3D.neps
+                    if enorm > Modia3D.nepsType(T)
                         return SVector(shapeLength, 0.0, 0.0) + SVector(0.0, topRadius*e_abs[2], topRadius*e_abs[3]) / enorm  # point on top circle is support point
                     else
                         return SVector(shapeLength, 0.0, 0.0)  # top circle center is support point
                     end
                 end
             else
                 enorm = norm(SVector(e_abs[2], e_abs[3]))
-                if enorm > Modia3D.neps
+                if enorm > Modia3D.nepsType(T)
                     return SVector(0.0, baseRadius*e_abs[2], baseRadius*e_abs[3]) / enorm  # point on base circle is support point
                 else
                     return SVector{3,T}(0.0, 0.0, 0.0)  # base circle center is support point
@@ -158,15 +158,15 @@ end
                     return SVector(0.0, shapeLength, 0.0)  # apex is support point
                 else  # frustum of a cone
                     enorm = norm(SVector(e_abs[3], e_abs[1]))
-                    if enorm > Modia3D.neps
+                    if enorm > Modia3D.nepsType(T)
                         return SVector(0.0, shapeLength, 0.0) + SVector(topRadius*e_abs[1], 0.0, topRadius*e_abs[3]) / enorm  # point on top circle is support point
                     else
                         return SVector(0.0, shapeLength, 0.0)  # top circle center is support point
                     end
                 end
             else
                 enorm = norm(SVector(e_abs[3], e_abs[1]))
-                if enorm > Modia3D.neps
+                if enorm > Modia3D.nepsType(T)
                     return SVector(baseRadius*e_abs[1], 0.0, baseRadius*e_abs[3]) / enorm  # point on base circle is support point
                 else
                     return SVector{3,T}(0.0, 0.0, 0.0)  # base circle center is support point
@@ -179,15 +179,15 @@ end
                     return SVector(0.0, 0.0, shapeLength)  # apex is support point
                 else  # frustum of a cone
                     enorm = norm(SVector(e_abs[1], e_abs[2]))
-                    if enorm > Modia3D.neps
+                    if enorm > Modia3D.nepsType(T)
                         return SVector(0.0, 0.0, shapeLength) + SVector(topRadius*e_abs[1], topRadius*e_abs[2], 0.0) / enorm  # point on top circle is support point
                     else
                         return SVector(0.0, 0.0, shapeLength)  # top circle center is support point
                     end
                 end
             else
                 enorm = norm(SVector(e_abs[1], e_abs[2]))
-                if enorm > Modia3D.neps
+                if enorm > Modia3D.nepsType(T)
                     return SVector(baseRadius*e_abs[1], baseRadius*e_abs[2], 0.0) / enorm  # point on base circle is support point
                 else
                     return SVector{3,T}(0.0, 0.0, 0.0)  # base circle center is support point
@@ -205,21 +205,21 @@ end
         halfThickness = T(0.5*shape.thickness)
         if shape.axis == 1
             enorm = norm(SVector(e_abs[1], e_abs[2]))
-            if enorm <= Modia3D.neps
+            if enorm <= Modia3D.nepsType(T)
                 return SVector(Basics.sign_eps(e_abs[1])*halfLength, 0.0, Basics.sign_eps(e_abs[3])*halfThickness)
             else
                 return SVector(Basics.sign_eps(e_abs[1])*halfLength, 0.0, Basics.sign_eps(e_abs[3])*halfThickness) + SVector(halfWidth*e_abs[1], halfWidth*e_abs[2], 0.0)/enorm
             end
         elseif shape.axis == 2
             enorm = norm(SVector(e_abs[2], e_abs[3]))
-            if enorm <= Modia3D.neps
+            if enorm <= Modia3D.nepsType(T)
                 return SVector(Basics.sign_eps(e_abs[1])*halfThickness, Basics.sign_eps(e_abs[2])*halfLength, 0.0)
             else
                 return SVector(Basics.sign_eps(e_abs[1])*halfThickness, Basics.sign_eps(e_abs[2])*halfLength, 0.0) + SVector(0.0, halfWidth*e_abs[2], halfWidth*e_abs[3])/enorm
             end
         else
             enorm = norm(SVector(e_abs[3], e_abs[1]))
-            if enorm <= Modia3D.neps
+            if enorm <= Modia3D.nepsType(T)
                 return SVector(0.0, Basics.sign_eps(e_abs[2])*halfThickness, Basics.sign_eps(e_abs[3])*halfLength)
             else
                 return SVector(0.0, Basics.sign_eps(e_abs[2])*halfThickness, Basics.sign_eps(e_abs[3])*halfLength) + SVector(halfWidth*e_abs[1], 0.0, halfWidth*e_abs[3])/enorm

src/Shapes/computePropertiesFileMes.jl

@@ -14,7 +14,6 @@
 #                         All faces should be specified in right-handed/counter-clockwise order!
 
 
-neps = sqrt(eps())
 
 function helpingFunc(w0,w1,w2)
     tmp0 = w0 + w1
@@ -70,7 +69,7 @@ function computeMassProperties(vertices::Vector{SVector{3,Float64}}, triangle_in
     end
 
     volume = integral[1]
-    if volume > neps
+    if volume > Modia3D.neps
         # center of volume
         centroid = [integral[2], integral[3], integral[4]]/volume
 

src/contactDetection/ContactDetectionMPR/ContactDetectionMPR_handler.jl

@@ -94,7 +94,7 @@ Base.:isequal(key1::DistanceKey, key2::DistanceKey) =
 
 
 """
-    handler = ContactDetectionMPR_handler(;tol_rel = 1e-4, niter_max=100, neps=sqrt(eps()))
+    handler = ContactDetectionMPR_handler(;tol_rel = 1e-7, niter_max=100)
 
 Generate a new contact handler for usage of the MPR algorithm
 The handler instance contains all information
@@ -105,8 +105,6 @@ about the contact situation.
 - `tol_rel`: Relative tolerance to compute the contact point (> 0.0)
 - `niter_max`: Maximum number of iterations of the MPR algorithm. If this number is reached,
                an error occurs (> 0).
-- `neps`: Small number used to check whether a floating number is close to zero (> 0.0).
-
 """
 mutable struct ContactDetectionMPR_handler{T} <: Modia3D.AbstractContactDetection
     distanceComputed::Bool
@@ -117,7 +115,6 @@ mutable struct ContactDetectionMPR_handler{T} <: Modia3D.AbstractContactDetectio
 
     tol_rel::T
     niter_max::Int
-    neps::T
 
     contactPairs::Composition.ContactPairs
 
@@ -127,12 +124,10 @@ mutable struct ContactDetectionMPR_handler{T} <: Modia3D.AbstractContactDetectio
     defaultContactSphereDiameter::Float64
 
     function ContactDetectionMPR_handler{T}(;tol_rel  = 1.0e-7,
-                                            niter_max = 100 ,
-                                            neps      = sqrt(eps(T))) where {T}
+                                            niter_max = 100) where {T}
         @assert(tol_rel > 0.0)
         @assert(niter_max > 0)
-        @assert(neps > 0.0)
-        new(false, Dict{PairID,ContactPair}(), Dict{PairID,ContactPair}(), 42.0, tol_rel, niter_max, neps)
+        new(false, Dict{PairID,ContactPair}(), Dict{PairID,ContactPair}(), 42.0, tol_rel, niter_max)
     end
 end
 ContactDetectionMPR_handler() = ContactDetectionMPR_handler{Modia3D.MPRFloatType}()

src/contactDetection/ContactDetectionMPR/analyzeMPR.jl

@@ -20,12 +20,13 @@ mutable struct Plane
 end
 
 
-function intersect3DSegmentPlane(seg::Segment, plane::Plane, neps)
+function intersect3DSegmentPlane(seg::Segment, plane::Plane)
    u = seg.P1 - seg.P0
    w = seg.P0 - plane.V0
 
    dividend = -dot(plane.n, w)
    divisor = dot(plane.n, u)
+   neps = Modia3D.nepsType(T)
 
    # checks if segment and plane are parallel
    if abs(divisor) < neps  # segment is parallel to plane
@@ -53,11 +54,11 @@ end
 pointInTriangle(p,a,b,c) = (sameSideTriangle(p,a,b,c) && sameSideTriangle(p,b,a,c) && sameSideTriangle(p,c,a,b))
 
 
-function doesRayIntersectPortal(r1,r2,r3,point,neps)
+function doesRayIntersectPortal(r1,r2,r3,point)
    plane = Plane(r1,r2,r3)
    segment = Segment(point, SVector(0.0, 0.0, 0.0))
 
-   (value, intersectionPoint) = intersect3DSegmentPlane(segment, plane, neps)
+   (value, intersectionPoint) = intersect3DSegmentPlane(segment, plane)
    if value == 1
       if !pointInTriangle(intersectionPoint, r1, r2, r3)
          error("Ray = ", point ," does not intersect portal (r1,r2,r3) = ", r1, r2, r3)
@@ -68,10 +69,10 @@ function doesRayIntersectPortal(r1,r2,r3,point,neps)
    return true
 end
 
-function analyzeFinalPortal(r1, r2, r3, r4, neps)
+function analyzeFinalPortal(r1, r2, r3, r4)
    plane = Plane(r1,r2,r3)
    segment = Segment(r4, SVector(0.0, 0.0, 0.0))
-   (value, intersectionPoint) = intersect3DSegmentPlane(segment, plane, neps)
+   (value, intersectionPoint) = intersect3DSegmentPlane(segment, plane)
    if value == 1
       if !pointInTriangle(intersectionPoint, r1, r2, r3)
          println("Ray r4 = ", r4 ," is not intersecting last portal (r1,r2,r3) = ", r1, r2, r3)