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

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

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

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This is an EXUDYN example
#
# Details:  Momentum conservation of ObjectContactSphereTorus (#2127): a free ball in the groove of a free
#           ring, both rigid bodies, no gravity; the ball hits the groove and the ring spins. The contact
#           forces are internal, so the total linear and angular momentum of the two bodies are conserved,
#           with and without friction. Measured: the angular momentum about the origin stays at its initial
#           value to 3e-9 (without friction) and 2e-9 (with friction), the linear momentum to 1e-15 - the
#           torque of the normal force on the ring, which acts at the contact point and not at the center of
#           the ring, was missing without friction (angular momentum off by 6e-3 of 0.09).
#
# Author:   Johannes Gerstmayr
# Date:     2026-09-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 numpy as np

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

def Run(mu, stepSize=1e-5, tEnd=0.05):
    """the ball and the ring; returns the times, the total angular and linear momentum"""
    SC=exu.SystemContainer(); mbs=SC.AddSystem()
    iB=InertiaSphere(mass=0.1, radius=0.01)
    iR=InertiaCylinder(density=7800, length=0.02, outerRadius=0.12, innerRadius=0.1, axis=2)
    nB=mbs.AddNode(NodeRigidBodyEP(referenceCoordinates=[0.1005,0,0]+eulerParameters0, initialVelocities=[0.5,1,0.2]+list(AngularVelocity2EulerParameters_t([0,0,50],eulerParameters0))))
    nR=mbs.AddNode(NodeRigidBodyEP(referenceCoordinates=[0,0,0]+eulerParameters0, initialVelocities=[0,0,0]+list(AngularVelocity2EulerParameters_t([0,2,3],eulerParameters0))))
    b=mbs.AddObject(ObjectRigidBody(nodeNumber=nB, mass=iB.Mass(), inertia=iB.GetInertia6D()))
    r=mbs.AddObject(ObjectRigidBody(nodeNumber=nR, mass=iR.Mass(), inertia=iR.GetInertia6D(), centerOfMass=iR.COM()))
    mB=mbs.AddMarker(MarkerBodyRigid(bodyNumber=b)); mR=mbs.AddMarker(MarkerBodyRigid(bodyNumber=r))
    nD=mbs.AddNode(NodeGenericData(numberOfDataCoordinates=4, initialCoordinates=[0,0,0,0]))
    mbs.AddObject(ObjectContactSphereTorus(markerNumbers=[mB,mR], nodeNumber=nD, sphereRadius=0.01, torusMajorRadius=0.1, torusMinorRadius=0.011,
        torusAxis=[0,0,1], contactStiffness=1e6, contactDamping=1e2, dynamicFriction=mu, frictionProportionalZone=1e-3))
    sens={}
    for name,n in [('B',nB),('R',nR)]:
        for ov in ['Position','Velocity','AngularVelocity','RotationMatrix']:
            sens[name+ov]=mbs.AddSensor(SensorNode(nodeNumber=n,outputVariableType=getattr(exu.OutputVariableType,ov),storeInternal=True))
    mbs.Assemble()
    s=exu.SimulationSettings(); s.timeIntegration.numberOfSteps=int(tEnd/stepSize); s.timeIntegration.endTime=tEnd
    s.timeIntegration.newton.useModifiedNewton = False #Just for the test; modified Newton is usually faster
    s.solution.file.write=False; s.solution.sensors.writePeriod=stepSize*10
    s.timeIntegration.generalizedAlpha.spectralRadius=1
    mbs.SolveDynamic(s)
    D={k:mbs.GetSensorStoredData(v)[:,1:] for k,v in sens.items()}
    t=mbs.GetSensorStoredData(sens['BPosition'])[:,0]
    def L(name, inertia, mass, com):
        J=inertia.Inertia()
        out=[]; P=[]
        for k in range(len(t)):
            A=D[name+'RotationMatrix'][k].reshape(3,3); w=D[name+'AngularVelocity'][k]
            p=D[name+'Position'][k]+A@com; v=D[name+'Velocity'][k]+np.cross(w,A@com)
            #inertia about COM
            Jc=J - mass*(np.dot(com,com)*np.eye(3)-np.outer(com,com))
            out.append(np.cross(p,mass*v)+A@Jc@A.T@w); P.append(mass*v)
        return np.array(out), np.array(P)
    LB,PB=L('B',iB,iB.Mass(),np.zeros(3)); LR,PR=L('R',iR,iR.Mass(),np.array(iR.COM()))
    Lt=LB+LR; Pt=PB+PR
    return t, Lt, Pt, np.abs(LB).max()


testResult = 0
for mu in [0, 0.3]:
    t, L, P, _ = Run(mu)
    errorL = np.abs(L-L[0]).max()
    errorP = np.abs(P-P[0]).max()
    exu.Print('friction', mu, ': angular momentum', L[-1].round(8), ', largest change', errorL, ', of the linear momentum', errorP)
    testResult += sum(L[-1]) + (errorL < 1e-7) + (errorP < 1e-12)

exu.Print('solution of contactSphereTorusMomentumTest=', testResult)
exu.sys['testResult'] = testResult
```