All test models changed to Model3D(..) and partially adapted. Some minor fixes in the code. Release notes added. Version number changed to 0.9.0.

Author: MartinOtter

docs/src/index.md

@@ -48,6 +48,28 @@ Download and install the free DLR SimVis Community Edition, e.g. with [https://v
 
 ## Release Notes
 
+### Version 0.9.0
+
+Non-backwards compatible changes
+
+- Operator `buildModia3D(..)` is removed. Instead, the new constructor `Model3D(..)` must be used at the top level of a
+  Modia3D definition. It is now possible to define several, independent multibody systems 
+  (currently, only one of them can have animation and animation export).
+
+- If init/start vectors are defined (e.g. initial state of a FreeMotion joint), they must be defined as SVector{3,Float64}(..).
+  Otherwise, errors occur during compilation.
+  
+- Requires ModiaLang Version 0.11
+
+
+Other changes
+
+- Code generation changed significantly, in order that the interface to Modia3D functions is type stable.
+  As a result, simulation is more efficient and much less memory is allocated during simulation.
+  
+- All test models changed, due to the non-backwards compatible change
+
+
 ### Version 0.8.2
 
 - Fix Cone contact detection

docs/src/internal/Profiling.md

@@ -14,7 +14,7 @@ The simplest technique is to put `@time` in front of the statement that should b
 
 module BouncingSphereSimulation_time
 
-@time begin BouncingSphere = Model(
+@time begin BouncingSphere = Model3D(
     ...
     )
 end
@@ -30,7 +30,7 @@ This gives the following output:
 ```julia
   0.000116 seconds (199 allocations: 12.344 KiB)
 
-Instantiating model Main.BouncingSphereSimulation.buildModia3D(BouncingSphere)
+Instantiating model Main.BouncingSphereSimulation.BouncingSphere
   1.106054 seconds (2.40 M allocations: 150.386 MiB, 3.92% gc time, 62.47% compilation time)
   0.322050 seconds (579.29 k allocations: 32.774 MiB, 5.81% compilation time)
   0.077613 seconds (1.74 k allocations: 97.844 KiB)
@@ -91,7 +91,7 @@ function `getDerivatives!(..)` is called twice to force compilation of this func
 The `logTiming=true` flag generates the following output in this case:
 
 ```
-Instantiating model Main.MobileWithLogTiming.buildModia3D(Mobile)
+Instantiating model Main.MobileWithLogTiming.Mobile
 
 Execute getDerivatives
 First executions of getDerivatives
@@ -101,9 +101,9 @@ First executions of getDerivatives
 ... first simulation:
   8.341341 seconds (37.34 M allocations: 2.236 GiB, 3.36% gc time)
 ... second simulation:
-... Simulate model buildModia3D(Mobile)
+... Simulate model Mobile
       Initialization at time = 0.0 s
-      Termination of buildModia3D(Mobile) at time = 5.0 s
+      Termination of Mobile at time = 5.0 s
         cpuTime         = 11.7 s
         allocated       = 2290.0 MiB
         algorithm       = CVODE_BDF
@@ -123,7 +123,7 @@ First executions of getDerivatives
         nStateEvents    = 0
         nRestartEvents  = 0
 
-... Timings for simulation of buildModia3D(Mobile):
+... Timings for simulation of Mobile:
  ──────────────────────────────────────────────────────────────────────────────────────────────
                                                        Time                   Allocations
                                                ──────────────────────   ───────────────────────

docs/src/tutorial/CollisionHandling.md

@@ -11,7 +11,7 @@ module BouncingSphere3D
 
 using Modia3D
 
-BouncingSphere = Model(
+BouncingSphere = Model3D(
     boxHeigth = 0.1,
     world      = Object3D(feature=Scene(enableContactDetection = true,
                     animationFile = "BouncingSphere.json")),
@@ -28,7 +28,7 @@ BouncingSphere = Model(
     free       = FreeMotion(obj1=:world, obj2=:sphere, r=Var(init=[0.0, 1.0, 0.0]))
 )
 
-bouncingSphere = @instantiateModel(buildModia3D(BouncingSphere), unitless=true)
+bouncingSphere = @instantiateModel(BouncingSphere, unitless=true)
 simulate!(bouncingSphere, stopTime=2.2, dtmax=0.1)
 
 @usingModiaPlot

docs/src/tutorial/GettingStarted.md

@@ -15,14 +15,14 @@ module Pendulum1
 
 using Modia3D
 
-Pendulum = Model(
+Pendulum = Model3D(
     world     = Object3D(feature=Scene()),
     body      = Object3D(feature=Solid(massProperties=MassProperties(mass=1.0))),
     bodyFrame = Object3D(parent=:body, translation=[-0.5, 0.0, 0.0]),
     rev       = Revolute(obj1=:world, obj2=:bodyFrame)
 )
 
-pendulum = @instantiateModel(buildModia3D(Pendulum), unitless=true)
+pendulum = @instantiateModel(Pendulum, unitless=true)
 simulate!(pendulum, stopTime=3.0)
 
 @usingModiaPlot     # use the plot package defined by ENV["MODIA_PLOT"]
@@ -38,7 +38,7 @@ julia> include("$(Modia3D.path)/test/Tutorial/Pendulum1.jl")
 
 The `world` [Object3D](@ref) has feature [Scene](@ref) and is therefore the inertial system. The `body` Object3D is a [Solid](@ref) and defines the pendulum as a mass point with `mass = 1.0`. The `bodyFrame` Object3D defines a coordinate system on the `body` that is translated along the x-axis. A revolute joint connects `world` with `bodyFrame`.
 
-With command `buildModia3D(Pendulum)`, the model definition is inspected and a few lines of code included, in order that joint variables are communicated between the Modia equations and the Modia3D multibody program. Keyword `unitless=true` defines that code generation is performed without units (because Modia3D does not yet fully support units in all components).
+With command `Model3D(..)`, a few lines of code are included, in order that joint variables are communicated between the Modia equations and the Modia3D multibody program. Keyword `unitless=true` defines that code generation is performed without units (because Modia3D does not yet fully support units in all components).
 
 The commands above generate an instance of the model, simulate it and generate the following plot:
 
@@ -54,7 +54,7 @@ module Pendulum2
 
 using Modia3D
 
-Pendulum = Model(
+Pendulum = Model3D(
     world = Object3D(feature=Scene(animationFile="Pendulum2.json")),
     obj1  = Object3D(feature=Solid(shape=Beam(axis=1, length=1.0, width=0.2, thickness=0.2),
                 solidMaterial="Steel", visualMaterial=VisualMaterial(color="Blue"))),
@@ -63,7 +63,7 @@ Pendulum = Model(
     rev   = Revolute(obj1=:world, obj2=:obj2)
 )
 
-pendulum = @instantiateModel(buildModia3D(Pendulum), unitless=true)
+pendulum = @instantiateModel(Pendulum, unitless=true)
 simulate!(pendulum, stopTime=3.0)
 
 @usingModiaPlot
@@ -97,7 +97,7 @@ using Modia3D
 
 include("$(Modia3D.modelsPath)/AllModels.jl")
 
-Pendulum = Model(
+Pendulum = Model3D(
     world = Object3D(feature=Scene(animationFile="Pendulum3.json")),
     obj1  = Object3D(feature=Solid(shape=Beam(axis=1, length=1.0, width=0.2, thickness=0.2),
                 solidMaterial="Steel", visualMaterial=VisualMaterial(color="Blue"))),
@@ -111,7 +111,7 @@ Pendulum = Model(
                 (damper.flange_a, fixed.flange)]
 )
 
-pendulum = @instantiateModel(buildModia3D(Pendulum), unitless=true)
+pendulum = @instantiateModel(Pendulum, unitless=true)
 simulate!(pendulum, stopTime=3.0)
 
 @usingModiaPlot

src/Composition/ForceElements/Bushing.jl

@@ -119,7 +119,7 @@ Bushing(; kwargs...) = Bushing{Float64}(; kwargs...)
 
 
 # Compute deformation angles from rotation matrix
-function anglesFromRotation(largeAngles::Bool, R12::SMatrix{3,3,F}, w12::SVector{3,F}) where F <: Modia3D.VarFloatType
+function anglesFromRotation(largeAngles::Bool, R12::SMatrix{3,3,F,9}, w12::SVector{3,F})::Tuple{SVector{3,F},SVector{3,F},SMatrix{2,2,F,4}} where F <: Modia3D.VarFloatType
     if largeAngles
         sbe = clamp(R12[3,1], F(-1.0), F(1.0))
         cbe = sqrt(F(1.0) - sbe*sbe)
@@ -132,10 +132,10 @@ function anglesFromRotation(largeAngles::Bool, R12::SMatrix{3,3,F}, w12::SVector
             ald = w12[1] + (sal*w12[2] - cal*w12[3])*sbe/cbe
             bed = cal*w12[2] + sal*w12[3]
             gad = (-sal*w12[2] + cal*w12[3])/cbe
-            return (SVector{3,F}(al, be, ga), SVector{3,F}(ald, bed, gad), SMatrix{2,2,F}(sal, cal, sbe, cbe))
+            return (SVector{3,F}(al, be, ga), SVector{3,F}(ald, bed, gad), SMatrix{2,2,F,4}(sal, cal, sbe, cbe))
         else
             @error("Gimbal lock of Bushing transformation.")
-            return (SVector{3,F}(0.0, 0.0, 0.0), SVector{3,F}(0.0, 0.0, 0.0), SMatrix{2,2,F}(0.0, 0.0, 0.0, 0.0))
+            return (SVector{3,F}(0.0, 0.0, 0.0), SVector{3,F}(0.0, 0.0, 0.0), SMatrix{2,2,F,4}(0.0, 0.0, 0.0, 0.0))
         end
     else
         al = R12[2,3]
@@ -147,12 +147,12 @@ function anglesFromRotation(largeAngles::Bool, R12::SMatrix{3,3,F}, w12::SVector
         if (max(abs(al), abs(be), abs(ga)) > 0.174)
             @warn("Bushing angle exceeds 10 deg.")
         end
-        return (SVector{3,F}(al, be, ga), SVector{3,F}(ald, bed, gad), SMatrix{2,2,F}(0.0, 0.0, 0.0, 0.0))
+        return (SVector{3,F}(al, be, ga), SVector{3,F}(ald, bed, gad), SMatrix{2,2,F,4}(0.0, 0.0, 0.0, 0.0))
     end
 end
 
 # Compute torque vector from force law moments
-function torqueFromMoments(largeAngles::Bool, moments::SVector{3,F}, sico::SMatrix{2,2,F}) where F <: Modia3D.VarFloatType
+function torqueFromMoments(largeAngles::Bool, moments::SVector{3,F}, sico::SMatrix{2,2,F,4})::SVector{3,F} where F <: Modia3D.VarFloatType
     if largeAngles
         tx = moments[1] + sico[1,2]*moments[3]
         ty = sico[2,1]*moments[2] - sico[1,1]*sico[2,2]*moments[3]
@@ -164,13 +164,13 @@ function torqueFromMoments(largeAngles::Bool, moments::SVector{3,F}, sico::SMatr
 end
 
 
-function initializeForceElement(force::Bushing{F}) where F <: Modia3D.VarFloatType
+function initializeForceElement(force::Bushing{F})::Nothing where F <: Modia3D.VarFloatType
     force.obj1.hasForceElement = true
     force.obj2.hasForceElement = true
     return nothing
 end
 
-function evaluateForceElement(force::Bushing{F}) where F <: Modia3D.VarFloatType
+function evaluateForceElement(force::Bushing{F})::Nothing where F <: Modia3D.VarFloatType
     R12 = measFrameRotation(force.obj2; frameOrig=force.obj1)
     r12 = measFramePosition(force.obj2; frameOrig=force.obj1, frameCoord=force.obj1)
     w12 = measFrameRotVelocity(force.obj2; frameOrig=force.obj1, frameCoord=force.obj1)
@@ -190,6 +190,6 @@ function evaluateForceElement(force::Bushing{F}) where F <: Modia3D.VarFloatType
     return nothing
 end
 
-function terminateForceElement(force::Bushing{F}) where F <: Modia3D.VarFloatType
+function terminateForceElement(force::Bushing{F})::Nothing where F <: Modia3D.VarFloatType
     return nothing
 end

src/Composition/dynamics.jl

@@ -370,7 +370,6 @@ function computeGeneralizedForces!(mbs::MultibodyData{F}, _leq)::MultibodyData{F
         exportAnimation = scene.exportAnimation
 
         if isTerminal  #if Modia.isTerminalOfAllSegments(sim)
-            println("... isTerminal = true")
             TimerOutputs.@timeit instantiatedModel.timer "Modia3D_4 isTerminal" begin
                 for force in forceElements
                     terminateForceElement(force)

src/Composition/joints/joints.jl

@@ -559,7 +559,7 @@ Copy accelerations of prismatic joints into mbs.
     @assert(NJOINTS == length(mbs.prismaticObjects))
     scene = mbs.scene
     @inbounds for (i,obj) in enumerate(mbs.prismaticObjects)
-        scene.prismatic[obj.jointIndex] = args[i]
+        scene.prismatic[obj.jointIndex].a = args[i]
     end
     return mbs
 end

src/ModiaInterface/_module.jl

@@ -25,6 +25,7 @@ import OrderedCollections
 import Modia3D
 import Modia3D.Composition
 using  ModiaLang
+using  StaticArrays
 
 
 # ModiaLang models

src/ModiaInterface/buildModia3D.jl

@@ -1,7 +1,8 @@
 using OrderedCollections
 
-appendSymbol(path::Nothing, name::Symbol) = name
-appendSymbol(path         , name::Symbol) = :( $path.$name )
+appendSymbol(path::Nothing, name::Union{Symbol,Expr}) = name
+appendSymbol(path         , name::Union{Symbol,Expr}) = :( $path.$name )
+appendSymbol(path::String , name::Union{Symbol,Expr}) = path == "" ? string(name) : path*"."*string(name)
 
 derSymbol(path::Nothing, name::Symbol) = :(der($name))
 derSymbol(path         , name::Symbol) = :(der($path.$name))
@@ -106,7 +107,7 @@ function buildModia3D!(model::AbstractDict, FloatType::Type, TimeType::Type,
         if jointType == :Revolute || jointType == :RevoluteWithFlange
             ndofTotal += 1
             NRevolute += 1
-            revoluteIndices[modelPathAsString*"."*string(path)] = NRevolute
+            revoluteIndices[appendSymbol(modelPathAsString,path)] = NRevolute
             push!(jointStatesRevolute, appendSymbol(path, :phi))
             push!(jointStatesRevolute, appendSymbol(path, :w))
             if jointType == :RevoluteWithFlange
@@ -125,7 +126,7 @@ function buildModia3D!(model::AbstractDict, FloatType::Type, TimeType::Type,
         elseif jointType == :Prismatic || jointType == :PrismaticWithFlange
             ndofTotal  += 1
             NPrismatic += 1
-            prismaticIndices[modelPathAsString*"."*string(path)] = NPrismatic            
+            prismaticIndices[appendSymbol(modelPathAsString,path)] = NPrismatic            
             push!(jointStatesPrismatic, appendSymbol(path, :s))
             push!(jointStatesPrismatic, appendSymbol(path, :v))
             if jointType == :PrismaticWithFlange
@@ -144,7 +145,7 @@ function buildModia3D!(model::AbstractDict, FloatType::Type, TimeType::Type,
         elseif jointType == :FreeMotion
             ndofTotal += 6
             NFreeMotion += 1
-            freeMotionIndices[modelPathAsString*"."*string(path)] = NFreeMotion           
+            freeMotionIndices[appendSymbol(modelPathAsString,path)] = NFreeMotion           
             push!(jointStatesFreeMotion, appendSymbol(path, :r))
             push!(jointStatesFreeMotion, appendSymbol(path, :v))
             push!(jointStatesFreeMotion, appendSymbol(path, :rot))            

src/ModiaInterface/model3D.jl

@@ -94,29 +94,21 @@ PrismaticWithFlange(; obj1, obj2, axis=1, s=Var(init=0.0), v=Var(init=0.0), canC
 
 Return joint rot. kinematics matrix `J` for Cardan angles `rot123` (rotation sequence x-y-z).
 """
-function J123(rot123::AbstractVector)
-
+function J123(rot123::SVector{3,F})::SMatrix{3,3,F,9} where {F}
     (sbe, cbe) = sincos(rot123[2])
     (sga, cga) = sincos(rot123[3])
-    return [     cga     -sga  0.0 ;
-             cbe*sga  cbe*cga  0.0 ;
-            -sbe*cga  sbe*sga  cbe ] / cbe
-
+    return SMatrix{3,3,F,9}(cga/cbe, sga, -sbe*cga/cbe, -sga/cbe, cga, sbe*sga/cbe, F(0.0), F(0.0), F(1.0))
 end
 
 """
     J = J132(rot132::AbstractVector)
 
 Return joint rot. kinematics matrix `J` for Cardan angles `rot132` (rotation sequence x-z-y).
 """
-function J132(rot132::AbstractVector)
-
+function J132(rot132::SVector{3,F})::SMatrix{3,3,F,9} where {F}
     (sga, cga) = sincos(rot132[2])
     (sbe, cbe) = sincos(rot132[3])
-    return [ cbe      0.0  sbe     ;
-            -sbe*cga  0.0  cbe*cga ;
-             cbe*sga  cga  sbe*sga ] / cga
-
+    return SMatrix{3,3,F,9}(cbe/cga, -sbe, cbe*sga/cga, F(0.0), F(0.0), F(1.0), sbe/cga, cbe, sbe*sga/cga)
 end