<!-- GENERATED by tools/generators/examplesDocsEmitter.py from python/ - do not edit -->
(testmodels-rigidbodyspringdamperintrinsic)=
# rigidBodySpringDamperIntrinsic.py

You can view and download this file on Github: [rigidBodySpringDamperIntrinsic.py](https://github.com/jgerstmayr/EXUDYN/tree/master/python/TestModels/rigidBodySpringDamperIntrinsic.py)

```{code-block} python
:linenos:

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This is an EXUDYN example
#
# Details:  Test rigid body formulation for different center of mass (COM)
#
# Author:   Johannes Gerstmayr
# Date:     2023-11-30
#
# Copyright:This file is part of Exudyn. Exudyn is free software. You can redistribute it and/or modify it under the terms of the Exudyn license. See 'LICENSE.txt' for more details.
#
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

import exudyn as exu
from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
import exudyn.graphics as graphics

import numpy as np

testIsActive = exu.sys.get('testIsActive', False)

SC = exu.SystemContainer()
mbs = SC.AddSystem()

#intrinsic is independent of order of markers; force/torque computed at mid-frame of joint
useIntrinsicFormulation=True

k=10000
d= k*0.05 * 0
kr = 500*0.5
dr = kr*0.05 * 0

L = 1
#++++++++++++++++++++++++++++++++++++++++++++++
#create two ground and two bodies, rigid body spring damper using different marker order
gGround = graphics.Brick(size = [L,0.1,0.25], color=graphics.color.blue)
oGround = mbs.CreateGround(referencePosition = [0,0,0], graphicsDataList=[gGround])

gBody1 = graphics.Brick(size = [L,0.1,0.2], color=graphics.color.red)
oBody1 = mbs.CreateRigidBody(referencePosition = [L,0.,0.],
                       referenceRotationMatrix = np.eye(3),
                       initialVelocity = [0.,1.,2.],
                       initialAngularVelocity = 5.*np.array([1.7,2.1,3.2]),
                       inertia=InertiaCuboid(density=1000,sideLengths=[L,0.1,0.2]),
                       gravity = [0.,0.,0.],
                       graphicsDataList = [gBody1],)
 
gBody2 = graphics.Brick(size = [L,0.1,0.2], color=graphics.color.green)
oBody2 = mbs.CreateRigidBody(referencePosition = [L,0.,0.],
                       referenceRotationMatrix = np.eye(3),
                       initialVelocity = [0.,1.,2.],
                       initialAngularVelocity = 5.*np.array([1.7,2.1,3.2]),
                       inertia=InertiaCuboid(density=1000,sideLengths=[L,0.1,0.2]),
                       gravity = [0.,0.,0.],
                       graphicsDataList = [gBody2],)
 
oRBSD1 = mbs.CreateRigidBodySpringDamper(itemNumbers = [oGround, oBody1],
                                        localPosition0 = [0.5*L,0.,0.], #global position
                                        localPosition1 = [-0.5*L,0.,0.], #global position
                                        stiffness = np.diag([k*0.4,k*0.5,k*0.7, kr,kr*2,kr*1.3]), 
                                        damping = np.diag([d,d,d, dr,dr,dr]), 
                                        offset = [0.,0.,0.,0.,0.,0.],
                                         useIntrinsicFormulation=useIntrinsicFormulation,
                                        )

oRBSD2 = mbs.CreateRigidBodySpringDamper(
                                         itemNumbers = [oBody2, oGround],
                                         localPosition0 = [-0.5*L,0.,0.], #global position
                                         localPosition1 = [0.5*L,0.,0.], #global position
                                         stiffness = np.diag([k*0.4,k*0.5,k*0.7, kr,kr*2,kr*1.3]), 
                                         damping = np.diag([d,d,d, dr,dr,dr]), 
                                         offset = [0.,0.,0.,0.,0.,0.],
                                         useIntrinsicFormulation=useIntrinsicFormulation,
                                        )

#++++++++++++++++++++++++++++++++++++++++++++++
#create two connected bodies, freely rotating:
oBody3a = mbs.CreateRigidBody(referencePosition = [0,1.,0.],
                       referenceRotationMatrix = np.eye(3),
                       initialVelocity = 0.*np.array([0.,1.,2.]),
                       initialAngularVelocity = 5*np.array([1.7,2.1,3.2]),
                       inertia=InertiaCuboid(density=1000,sideLengths=[L,0.1,0.2]),
                       gravity = [0.,0.,0.],
                       graphicsDataList = [gBody1],)
 
oBody3b = mbs.CreateRigidBody(referencePosition = [L,1.,0.],
                       referenceRotationMatrix = np.eye(3),
                       initialVelocity = 0.*np.array([0.,-1.,-2.]),
                       initialAngularVelocity = -5*np.array([1.7,2.1,3.2]),
                       inertia=InertiaCuboid(density=1000,sideLengths=[L,0.1,0.2]),
                       gravity = [0.,0.,0.],
                       graphicsDataList = [gBody2],)

oRBSD3 = mbs.CreateRigidBodySpringDamper(
                                         itemNumbers = [oBody3a, oBody3b],
                                         localPosition0 = [0.5*L,0.,0.], #global position
                                         localPosition1 = [-0.5*L,0.,0.], #global position
                                         stiffness = 1*np.diag([k*0.4,k*0.5,k*0.7, kr,kr*2,kr*1.3]), 
                                         damping = np.diag([d,d,d, dr,dr,dr]), 
                                         offset = [0.,0.,0.,0.,0.,0.],
                                         useIntrinsicFormulation=useIntrinsicFormulation,
                                        )

# mbs.SetObjectParameter(oRBSD1, 'useIntrinsicFormulation', True)
# mbs.SetObjectParameter(oRBSD2, 'useIntrinsicFormulation', True)
# mbs.SetObjectParameter(oRBSD3, 'useIntrinsicFormulation', True)

sBody1 = mbs.AddSensor(SensorBody(bodyNumber=oBody1, storeInternal=True, outputVariableType=exu.OutputVariableType.Position))
sBody2 = mbs.AddSensor(SensorBody(bodyNumber=oBody2, storeInternal=True, outputVariableType=exu.OutputVariableType.Position))
sBody3 = mbs.AddSensor(SensorBody(bodyNumber=oBody3a, storeInternal=True, outputVariableType=exu.OutputVariableType.AngularVelocity))


#+++++++++++++++++++++++++++++++++++++++++++++++++++++++
mbs.Assemble()

endTime = 1 #for stepSize=0.002: non-intrinsic formulation gets unstable around 7.5 seconds for body3 and for body1/2 around 73 seconds
if not testIsActive:
    endTime = 100

stepSize = 0.002
simulationSettings = exu.SimulationSettings()
#simulationSettings.show.computationTime = True
simulationSettings.timeIntegration.verboseMode = (not testIsActive)
simulationSettings.solution.file.write = (not testIsActive)

simulationSettings.timeIntegration.numberOfSteps = int(endTime/stepSize)
simulationSettings.timeIntegration.endTime = endTime
simulationSettings.timeIntegration.newton.useModifiedNewton = True


if not testIsActive:
    SC.renderer.Start()
    SC.renderer.DoIdleTasks()

mbs.SolveDynamic(simulationSettings)

if not testIsActive:
    SC.renderer.Stop() #safely close rendering window!


    mbs.PlotSensor([sBody1,sBody2], components=[1,1])
    mbs.PlotSensor([sBody1,sBody2], components=[2,2])
    mbs.PlotSensor([sBody3], components=[0,1,2])

    mbs.SolutionViewer()


p1=mbs.GetObjectOutputBody(oBody1, exu.OutputVariableType.Displacement)
exu.Print("p1=", p1)
p2=mbs.GetObjectOutputBody(oBody2, exu.OutputVariableType.Displacement)
exu.Print("p2=", p2)
omega3=mbs.GetObjectOutputBody(oBody3a, exu.OutputVariableType.AngularVelocity)
exu.Print("omega3=", omega3)


#+++++++++++++++++++++++++++++++++++++++++++++
u=np.linalg.norm(p1) + np.linalg.norm(p2) + np.linalg.norm(0.01*omega3)
exu.Print('solution of rigidBodySpringDamperIntrinsic test=',u)

exu.sys['testResult'] = u


if not testIsActive:
    SC.renderer.DoIdleTasks()
    SC.renderer.Stop() #safely close rendering window!
```