Skip to content

Commit d750083

Browse files
Add result element ContactResult
- Add result element for MPR contact results. - Function `getElasticContactPair` based on implementation of [email protected]. - Add comments and documentation.
1 parent 504a0bd commit d750083

File tree

9 files changed

+230
-27
lines changed

9 files changed

+230
-27
lines changed

docs/make.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -22,6 +22,7 @@ makedocs(
2222
"Components/Materials.md"
2323
"Components/GravityField.md"
2424
"Components/ForceElements.md"
25+
"Components/ResultElements.md"
2526
],
2627
"Functions" => "Functions.md",
2728
"Internal" => [
Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,11 @@
1+
# Result Elements
2+
3+
```@meta
4+
CurrentModule = Modia3D.Composition
5+
```
6+
7+
## ContactResult
8+
9+
```@docs
10+
ContactResult
11+
```
Lines changed: 132 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,132 @@
1+
2+
"""
3+
result = ContactResult(; object1, object2, objectCoordinateRef)
4+
5+
Return a `result` providing results of elastic MPR contacts between
6+
`object1::`[`Object3D`](@ref) and `object2::`[`Object3D`](@ref).
7+
If `objectCoordinateRef::`[`Object3D`](@ref) is defined, all vector
8+
results are resolved in `objectCoordinateRef`, otherwise in world.
9+
10+
# Results
11+
12+
- `penetration` is the normal contact penetration (positive in case of contact).
13+
- `penetrationVelocity` is the normal contact penetration velocity (positive for compression).
14+
- `tangentialVelocity` is the absolute value of the tangential relative velocity.
15+
- `angularVelocity` it the absolute value of the relative angular velocity.
16+
- `normalForce` is the normal contact force (positive for pressure).
17+
- `tangentialForce` is the absolute value of the tangential contact force.
18+
- `torque` is the absolute value of the contact torque.
19+
- `positionVector` is the absolute position vector of the contact point, resolved in `objectCoordinateRef`.
20+
- `normalVector` is the unit vector in contact normal direction, pointing into `object1`, resolved in `objectCoordinateRef`.
21+
- `forceVector` is the total contact force vector acting at the contact point on object1, resolved in `objectCoordinateRef`.
22+
On `object2` the same force vector is applied in inverse direction.
23+
- `torqueVector` is the total contact torque vector acting at the contact point on object1, resolved in `objectCoordinateRef`.
24+
On `object2` the same torque vector is applied in inverse direction.
25+
"""
26+
mutable struct ContactResult{F <: Modia3D.VarFloatType} <: Modia3D.AbstractResultElement
27+
28+
path::String
29+
30+
object1::Object3D{F}
31+
object2::Object3D{F}
32+
objectCoordinateRef::Union{Object3D{F}, Nothing}
33+
34+
penetrationResultIndex::Int
35+
penetrationVelocityResultIndex::Int
36+
tangentialVelocityResultIndex::Int
37+
angularVelocityResultIndex::Int
38+
normalForceResultIndex::Int
39+
tangentialForceResultIndex::Int
40+
torqueResultIndex::Int
41+
positionVectorResultIndex::Int
42+
normalVectorResultIndex::Int
43+
forceVectorResultIndex::Int
44+
torqueVectorResultIndex::Int
45+
46+
function ContactResult{F}(; path::String = "",
47+
object1::Object3D{F},
48+
object2::Object3D{F},
49+
objectCoordinateRef::Union{Object3D{F}, Nothing}=nothing ) where F <: Modia3D.VarFloatType
50+
return new(path, object1, object2, objectCoordinateRef)
51+
end
52+
end
53+
ContactResult(; kwargs...) = ContactResult{Float64}(; kwargs...)
54+
55+
56+
function initializeResultElement(model::Modia.InstantiatedModel{F,TimeType}, result::ContactResult{F}) where {F <: Modia3D.VarFloatType, TimeType <: AbstractFloat}
57+
58+
result.penetrationResultIndex = Modia.new_w_segmented_variable!(model, result.path*".penetration" , F(0), "m")
59+
result.penetrationVelocityResultIndex = Modia.new_w_segmented_variable!(model, result.path*".penetrationVelocity", F(0), "m/s")
60+
result.tangentialVelocityResultIndex = Modia.new_w_segmented_variable!(model, result.path*".tangentialVelocity" , F(0), "m/s")
61+
result.angularVelocityResultIndex = Modia.new_w_segmented_variable!(model, result.path*".angularVelocity" , F(0), "rad/s")
62+
result.normalForceResultIndex = Modia.new_w_segmented_variable!(model, result.path*".normalForce" , F(0), "N")
63+
result.tangentialForceResultIndex = Modia.new_w_segmented_variable!(model, result.path*".tangentialForce" , F(0), "N")
64+
result.torqueResultIndex = Modia.new_w_segmented_variable!(model, result.path*".torque" , F(0), "N*m")
65+
result.positionVectorResultIndex = Modia.new_w_segmented_variable!(model, result.path*".positionVector" , SVector{3,F}(0, 0, 0), "m")
66+
result.normalVectorResultIndex = Modia.new_w_segmented_variable!(model, result.path*".normalVector" , SVector{3,F}(0, 0, 0))
67+
result.forceVectorResultIndex = Modia.new_w_segmented_variable!(model, result.path*".forceVector" , SVector{3,F}(0, 0, 0), "N")
68+
result.torqueVectorResultIndex = Modia.new_w_segmented_variable!(model, result.path*".torqueVector" , SVector{3,F}(0, 0, 0), "N*m")
69+
70+
return nothing
71+
end
72+
73+
function evaluateResultElement(model::Modia.InstantiatedModel{F,TimeType}, scene::Modia3D.Composition.Scene{F}, result::ContactResult{F}, time::TimeType) where {F <: Modia3D.VarFloatType, TimeType <: AbstractFloat}
74+
75+
(contactPair, converse) = getElasticContactPair(scene, result.object1, result.object2)
76+
if !isnothing(contactPair)
77+
penetration = contactPair.results.penetration
78+
penetrationVelocity = contactPair.results.penetrationVelocity
79+
tangentialVelocity = contactPair.results.tangentialVelocity
80+
angularVelocity = contactPair.results.angularVelocity
81+
normalForce = contactPair.results.normalForce
82+
tangentialForce = contactPair.results.tangentialForce
83+
torque = contactPair.results.torque
84+
positionVector = contactPair.results.positionVector
85+
normalVector = contactPair.results.normalVector
86+
forceVector = contactPair.results.forceVector
87+
torqueVector = contactPair.results.torqueVector
88+
if converse
89+
normalVector = -normalVector
90+
forceVector = -forceVector
91+
torqueVector = -torqueVector
92+
end
93+
if !isnothing(result.objectCoordinateRef)
94+
positionVector = result.objectCoordinateRef.R_abs * positionVector
95+
normalVector = result.objectCoordinateRef.R_abs * normalVector
96+
forceVector = result.objectCoordinateRef.R_abs * forceVector
97+
torqueVector = result.objectCoordinateRef.R_abs * torqueVector
98+
end
99+
else
100+
penetration = F(0)
101+
penetrationVelocity = F(0)
102+
tangentialVelocity = F(0)
103+
angularVelocity = F(0)
104+
normalForce = F(0)
105+
tangentialForce = F(0)
106+
torque = F(0)
107+
positionVector = SVector{3,F}(0, 0, 0)
108+
normalVector = SVector{3,F}(0, 0, 0)
109+
forceVector = SVector{3,F}(0, 0, 0)
110+
torqueVector = SVector{3,F}(0, 0, 0)
111+
end
112+
113+
if Modia.storeResults(model)
114+
Modia.copy_w_segmented_value_to_result(model, result.penetrationResultIndex, penetration)
115+
Modia.copy_w_segmented_value_to_result(model, result.penetrationVelocityResultIndex, penetrationVelocity)
116+
Modia.copy_w_segmented_value_to_result(model, result.tangentialVelocityResultIndex, tangentialVelocity)
117+
Modia.copy_w_segmented_value_to_result(model, result.angularVelocityResultIndex, angularVelocity)
118+
Modia.copy_w_segmented_value_to_result(model, result.normalForceResultIndex, normalForce)
119+
Modia.copy_w_segmented_value_to_result(model, result.tangentialForceResultIndex, tangentialForce)
120+
Modia.copy_w_segmented_value_to_result(model, result.torqueResultIndex, torque)
121+
Modia.copy_w_segmented_value_to_result(model, result.positionVectorResultIndex, positionVector)
122+
Modia.copy_w_segmented_value_to_result(model, result.normalVectorResultIndex, normalVector)
123+
Modia.copy_w_segmented_value_to_result(model, result.forceVectorResultIndex, forceVector)
124+
Modia.copy_w_segmented_value_to_result(model, result.torqueVectorResultIndex, torqueVector)
125+
end
126+
127+
return nothing
128+
end
129+
130+
function terminateResultElement(result::ContactResult{F}) where F <: Modia3D.VarFloatType
131+
return nothing
132+
end

src/Composition/_module.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -117,6 +117,7 @@ include("scene.jl") # must be included after superObjects.jl
117117
include(joinpath("joints", "joints.jl"))
118118
include(joinpath("joints", "changeJointState.jl"))
119119

120+
include(joinpath("ResultElements", "ContactResult.jl"))
120121
include("sensors.jl")
121122
include("frameMeasurements.jl")
122123
include("frameForceTorque.jl")

src/Composition/dynamicCollision.jl

Lines changed: 21 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -49,8 +49,8 @@ function dealWithContacts!(sim::Modia.InstantiatedModel{F, T}, scene::Scene{F},
4949

5050
for (pairID::Int64, pair::ContactPair{F}) in ch.contactDict
5151
obj1 = pair.obj1
52-
obj2= pair.obj2
53-
rContact= (pair.contactPoint1 + pair.contactPoint2)/F(2.0)
52+
obj2 = pair.obj2
53+
rContact = (pair.contactPoint1 + pair.contactPoint2)/F(2.0)
5454
contactNormal = pair.contactNormal
5555
if Modia.isEvent(sim)
5656
# println("$(sim.time): ", obj1.path, " ", obj2.path)
@@ -65,7 +65,7 @@ function dealWithContacts!(sim::Modia.InstantiatedModel{F, T}, scene::Scene{F},
6565
pair.distanceWithHysteresis, time, file, sim)
6666
elseif pair.pairKind == Modia3D.ElasticContactPairKind
6767
elasticContactPairMaterial::Composition.ElasticContactPairResponseMaterial = pair.contactPairMaterial
68-
(f1,f2,t1,t2) = responseCalculation(elasticContactPairMaterial, obj1, obj2,
68+
(f1,f2,t1,t2,pair.results) = responseCalculation(elasticContactPairMaterial, obj1, obj2,
6969
rContact, contactNormal,
7070
pair.distanceWithHysteresis, time, file, sim)
7171
elseif pair.pairKind == Modia3D.ObserverContactPairKind
@@ -245,3 +245,21 @@ function setVisualizationContactProperties!(obj::Composition.Object3D{F}, transp
245245
obj.feature.visualMaterial.transparency = transparency
246246
return nothing
247247
end
248+
249+
250+
function getElasticContactPair(scene::Scene{F}, obj1::Composition.Object3D{F}, obj2::Composition.Object3D{F}) where F <: Modia3D.VarFloatType
251+
252+
if scene.options.enableContactDetection
253+
for (pairID::Int64, pair::ContactPair{F}) in scene.options.contactDetection.contactDict
254+
if pair.pairKind == Modia3D.ElasticContactPairKind
255+
if (obj1 === pair.obj1 && obj2 === pair.obj2)
256+
return (pair, false) # object 1/2 assignment consistent
257+
elseif (obj2 === pair.obj1 && obj1 === pair.obj2)
258+
return (pair, true) # object 1/2 assignment converse
259+
end
260+
end
261+
end
262+
end
263+
264+
return (nothing, false)
265+
end

src/Composition/responseCalculation/elasticCollisionResponse.jl

Lines changed: 34 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -126,7 +126,7 @@ at time `time`.
126126
"""
127127
function responseCalculation(material::ElasticContactPairResponseMaterial, obj1::Object3D{F}, obj2::Object3D{F},
128128
rContact::SVector{3,F}, e_n::SVector{3,F},
129-
s::F, time, file, sim)::Tuple{SVector{3,F},SVector{3,F},SVector{3,F},SVector{3,F}} where F <: Modia3D.VarFloatType
129+
s::F, time, file, sim)::Tuple{SVector{3,F},SVector{3,F},SVector{3,F},SVector{3,F},ContactResults{F}} where F <: Modia3D.VarFloatType
130130
# Material
131131
c_res = F(material.c_res)
132132
c_geo = F(material.c_geo)
@@ -140,23 +140,32 @@ function responseCalculation(material::ElasticContactPairResponseMaterial, obj1:
140140

141141
### signed velocity and relative velocity ####
142142
# Contact points and distances to local part frame (in world frame)
143+
# r_rel1 ... position vector from obj1 to average contact point, resolved in world
144+
# r_rel2 ... position vector from obj2 to average contact point, resolved in world
143145
r_rel1 = rContact - obj1.r_abs
144146
r_rel2 = rContact - obj2.r_abs
145147

146148
# Velocities and angular velocities of contact frames in world frame
149+
147150
w1 = obj1.R_abs'*obj1.w
148151
w2 = obj2.R_abs'*obj2.w
149152
v1 = obj1.v0 + cross(w1,r_rel1)
150153
v2 = obj2.v0 + cross(w2,r_rel2)
151154

152155
# Velocities and angular velocities in normal and tangential direction
156+
# w_rel ... angular velocity vector of obj2 relative to obj1, resolved in world
157+
# e_w_reg ... regularised direction vector of w_rel, resolved in world
158+
# v_rel ... velocity vector at average contact point on obj2 relative to obj1, resolved in world
153159
w_rel = w2 - w1
154160
e_w_reg = w_rel/Modia3D.regularize(norm(w_rel),wsmall)
155161
v_rel = v2 - v1
156162

157-
# delta_dot ... signed relative velocity in normal direction
158-
# v_t ... relative velocity vector in tangential direction
159-
# delta ... signed distance
163+
# e_n ... points from contact point on obj2 to contact point on obj1, i.e. normal direction into obj2
164+
# s ... signed distance of the contact points, negative during contact
165+
# delta_dot ... signed relative velocity in normal direction, positive for expansion (sign inconsistent with delta!)
166+
# v_t ... tangential velocity vector at average contact point on obj2 relative to obj1, resolved in world
167+
# e_t_reg ... regularised direction vector of v_t, resolved in world
168+
# delta ... penetration, positive during contact
160169
delta_dot = dot(v_rel,e_n)
161170
v_t = v_rel - delta_dot*e_n
162171
e_t_reg = v_t/Modia3D.regularize(norm(v_t),vsmall)
@@ -165,16 +174,27 @@ function responseCalculation(material::ElasticContactPairResponseMaterial, obj1:
165174
delta_comp = delta * sqrt(abs(delta))
166175

167176
#fn = -max(F(0.0), c_res * c_geo * delta_comp * (1 - d_res*delta_dot) )
168-
fn = -c_res * c_geo * delta_comp * (1 - d_res*delta_dot)
169-
ft = -mu_k*fn*e_t_reg
170-
f1 = fn * e_n + ft
171-
f2 = -f1
172-
173-
tau = -mu_r * mu_r_geo * fn * e_w_reg
174-
t1 = cross(r_rel1,f1) + tau
175-
t2 = cross(r_rel2,f2) - tau
176-
177-
return (f1,f2,t1,t2)
177+
fn = -c_res * c_geo * delta_comp * (1 - d_res*delta_dot) # normal force, negative for pressure!
178+
ft = -mu_k*fn*e_t_reg # tangential force vector acting on obj1, resolved in world
179+
f1 = fn * e_n + ft # total force vector acting on obj1, resolved in world
180+
f2 = -f1 # total force vector acting on obj2, resolved in world
181+
182+
tau = -mu_r * mu_r_geo * fn * e_w_reg # torque vector acting at average contact point on obj1, resolved in world
183+
t1 = cross(r_rel1,f1) + tau # torque vector acting at obj1, resolved in world
184+
t2 = cross(r_rel2,f2) - tau # torque vector acting at obj2, resolved in world
185+
186+
return (f1, f2, t1, t2,
187+
ContactResults(delta, # normal contact penetration (positive during contact)
188+
-delta_dot, # normal contact penetration velocity (positive for compression)
189+
norm(v_t), # absolute value of the tangential relative velocity
190+
norm(w_rel), # absolute value of the relative angular velocity
191+
-fn, # normal contact force (positive for pressure)
192+
norm(ft), # absolute value of the tangential contact force
193+
norm(tau), # absolute value of the contact torque
194+
rContact, # absolute position vector of the contact point, resolved in world
195+
-e_n, # unit vector in contact normal direction, pointing into obj1, resolved in world
196+
f1, # total contact force vector acting at the contact point on obj1, resolved in world
197+
tau)) # total contact torque vector acting at the contact point on obj1, resolved in world
178198
end
179199

180200
responseCalculation(material::Nothing, obj1::Object3D{F}, obj2::Object3D{F},

src/ModiaInterface/_module.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,7 @@ export Revolute, RevoluteWithFlange
1919
export Prismatic, PrismaticWithFlange
2020
export singularRem, FreeMotion, change_rotSequenceInNextIteration!
2121
export WorldForce, WorldTorque, Bushing, SpringDamperPtP, PolygonalContactModel
22+
export ContactResult
2223
export FreeMotion2
2324

2425
export ModelActions, ActionAttach, ActionRelease, ActionReleaseAndAttach, ActionDelete, EventAfterPeriod, ActionWait

src/ModiaInterface/model3D.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -30,6 +30,7 @@ WorldTorque( ; kwargs...) = Par(; _constructor = :(Modia3D.Composition.
3030
Bushing( ; kwargs...) = Par(; _constructor = :(Modia3D.Composition.Bushing{FloatType}) , _path = true, kwargs...)
3131
SpringDamperPtP( ; kwargs...) = Par(; _constructor = :(Modia3D.Composition.SpringDamperPtP{FloatType}) , _path = true, kwargs...)
3232
PolygonalContactModel(; kwargs...) = Par(; _constructor = :(Modia3D.Composition.PolygonalContactModel{FloatType}), _path = true, kwargs...)
33+
ContactResult( ; kwargs...) = Par(; _constructor = :(Modia3D.Composition.ContactResult{FloatType}) , _path = true, kwargs...)
3334

3435
MassPropertiesFromShape() = Par(; _constructor = :(Modia3D.Shapes.MassPropertiesFromShape{FloatType}))
3536
MassPropertiesFromShapeAndMass(; mass) = Par(; _constructor = :(Modia3D.Shapes.MassPropertiesFromShapeAndMass{FloatType}), mass = mass)

src/contactDetection/ContactDetectionMPR/ContactDetectionMPR_handler.jl

Lines changed: 28 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,22 @@ using OrderedCollections
1111
const PairID = Int64
1212

1313

14+
# Struct for contact results, used by ContactResult element
15+
struct ContactResults{F <: Modia3D.VarFloatType}
16+
penetration::F # normal contact penetration (positive during contact)
17+
penetrationVelocity::F # normal contact penetration velocity (positive for compression)
18+
tangentialVelocity::F # absolute value of the tangential relative velocity
19+
angularVelocity::F # absolute value of the relative angular velocity
20+
normalForce::F # normal contact force (positive for pressure)
21+
tangentialForce::F # absolute value of the tangential contact force
22+
torque::F # absolute value of the contact torque
23+
positionVector::SVector{3,F} # absolute position vector of the contact point, resolved in world
24+
normalVector::SVector{3,F} # unit vector in contact normal direction, pointing into obj1, resolved in world
25+
forceVector::SVector{3,F} # total contact force vector acting at the contact point on obj1, resolved in world
26+
torqueVector::SVector{3,F} # total contact torque vector acting at the contact point on obj1, resolved in world
27+
end
28+
29+
1430
"""
1531
ContactPair(contactPointA, contactPointB,
1632
contactNormal, objA, objB,
@@ -21,24 +37,26 @@ const PairID = Int64
2137
Generate a new `ContactPair` object of two objects that have `contact = true`.
2238
"""
2339
mutable struct ContactPair{F <: Modia3D.VarFloatType}
24-
contactPoint1::SVector{3,F}
25-
contactPoint2::SVector{3,F}
26-
contactNormal::SVector{3,F}
40+
contactPoint1::SVector{3,F} # absolute position of contact point on object1 resolved in world
41+
contactPoint2::SVector{3,F} # absolute position of contact point on object2 resolved in world
42+
contactNormal::SVector{3,F} # vector pointing from contactPoint2 to contactPoint1 resolved in world
2743
obj1::Object3D{F}
2844
obj2::Object3D{F}
29-
distanceWithHysteresis::F
45+
distanceWithHysteresis::F # signed distance: positive: Euclidean distance (no contact), negative: penetration depth (contact)
3046

3147
supportPointsDefined::Bool
32-
support1A::SVector{3,F}
33-
support2A::SVector{3,F}
34-
support3A::SVector{3,F}
35-
support1B::SVector{3,F}
36-
support2B::SVector{3,F}
37-
support3B::SVector{3,F}
48+
support1A::SVector{3,F} # absolute position of support point 1 on obj1 resolved in world
49+
support2A::SVector{3,F} # absolute position of support point 2 on obj1 resolved in world
50+
support3A::SVector{3,F} # absolute position of support point 3 on obj1 resolved in world
51+
support1B::SVector{3,F} # absolute position of support point 1 on obj2 resolved in world
52+
support2B::SVector{3,F} # absolute position of support point 2 on obj2 resolved in world
53+
support3B::SVector{3,F} # absolute position of support point 3 on obj2 resolved in world
3854

3955
pairKind::Shapes.PairMaterialKind
4056
contactPairMaterial::Union{Modia3D.AbstractContactPairMaterial,Nothing} # only if contact = true, otherwise not defined
4157

58+
results::ContactResults{F}
59+
4260
ContactPair{F}(contactPoint1::SVector{3,F}, contactPoint2::SVector{3,F},
4361
contactNormal::SVector{3,F},
4462
obj1::Object3D{F}, obj2::Object3D{F}, distanceWithHysteresis::F, supportPointsDefined::Bool,

0 commit comments

Comments
 (0)