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

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

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

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This is an EXUDYN example
#
# Details:  Test model for MainSystem Create functions; the model itself is not very meaningful, but it uses many Create functions
#
# Author:   Johannes Gerstmayr
# Date:     2025-05-11
#
# 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)


#set up new multibody system to work with
SC = exu.SystemContainer()
mbs = SC.AddSystem()

graphicsGround = graphics.CheckerBoard(point=[0,2,-1], size=8)
oGround = mbs.CreateGround(referencePosition=[0,0,0], graphicsDataList=[graphicsGround])

oMass = mbs.CreateMassPoint(mass=5, referencePosition=[1,0,0], 
                            initialDisplacement=[0,0,0],
                            initialVelocity=[0,0.5,0],
                            gravity=[0,-9.81,0],
                            drawSize=0.3, color=graphics.color.orange,
                            )

oSpringDamper = mbs.CreateSpringDamper(itemNumbers=[oGround, oMass], #[body0,body1]
                                      localPosition0=[0,0,0], #locally on body0
                                      localPosition1=[0,0,0], #locally on body1
                                      referenceLength=None, #usually set to None (default) => takes the (relaxed) length in reference configuration
                                      stiffness=1e2, damping=1)

oDistance = mbs.CreateDistanceConstraint(itemNumbers=[oGround, oMass], #[body0,body1]
                                         localPosition0 = [0,0,0], #locally on body0
                                         localPosition1 = [0,0,0], #locally on body1
                                         distance=None)

#%%++++++++++++++++++++++++++++++++++++++++++++++++++

mass = 10 
outerRadius = 0.45
length = 2
volume = np.pi * outerRadius**2 * length
inertiaCylinder = InertiaCylinder(density=mass/volume, length=length, outerRadius=outerRadius, 
                                  axis=0, innerRadius=0.7*outerRadius)

oBody = mbs.CreateRigidBody(inertia=inertiaCylinder,
                            referencePosition=[3,0,0],
                            initialVelocity=[0,0,0],
                            initialAngularVelocity=[2*np.pi,0,0],
                            gravity=[0,-9.81,0],
                            #graphicsDataList=[graphicsCyl]
                            )


rbX = 2
rbY = 0.5
rbZ = 0.2

#create inertia instance for cuboid (brick) shape
mass=10 #kg
volume=rbX*rbY*rbZ
inertiaCube2 = InertiaCuboid(density=mass/volume, sideLengths=[rbX,rbY,rbZ])
inertiaCube2 = inertiaCube2.Translated([0.5,0,0])

oBody2 = mbs.CreateRigidBody(inertia = inertiaCube2,
                             referencePosition = [4,-0.5*rbX,0], 
                             referenceRotationMatrix = RotationVector2RotationMatrix([0,0,-0.5*np.pi]),
                             initialAngularVelocity = [0,0,0.2*2*np.pi],
                             initialVelocity = [0.2*np.pi*0.5*rbX,0,0],
                             gravity = [0,-9.81,0],
                             color=graphics.color.blue,
                             #graphicsDataList=[graphicsCube, graphics.Basis(origin=inertiaCube2.COM())],
                             )

loadMassPoint = mbs.CreateForce(itemNumber=oMass, loadVector=[10,0,0])
loadRigidBody = mbs.CreateForce(itemNumber=oBody, localPosition=[0,0,0], loadVector=[10,0,0])

def UFforce(mbs, t, loadVector):
    return (10+5*np.sin(t*10*2*np.pi))*np.array([0,5,0])

mbs.CreateForce(itemNumber=oBody,
                localPosition=[0,1.2,0.5], #position at body
                loadVectorUserFunction=UFforce)

def UserFunctionTorque(mbs, t, loadVector):
    return np.cos(t*2*np.pi)*np.array(loadVector)

#add torque with user function
mbs.CreateTorque(itemNumber=oBody, loadVector=[5,0,0], 
                 loadVectorUserFunction=UserFunctionTorque)

#+++++++++
inertiaSphere = InertiaSphere(mass=2, radius=0.25)
gSphere = inertiaSphere.GetGraphics(graphics.color.red)

oMass1 = mbs.CreateMassPoint(mass=2, referencePosition=[1,0,0],
                             gravity=[0,-9.81,0],
                             # drawSize=0.25, color=graphics.color.red,
                             graphicsDataList=[gSphere]
                             )

#create spherical joint between ground and mass point; could also be applied to two bodies; possible bodies: mass point or rigid body
mbs.CreateSphericalJoint(itemNumbers=[oGround, oMass1],
                         position=[1,0,0], #global position of joint (in reference configuration)
                         constrainedAxes=[1,0,0]) #x,y,z directions: 1 = fixed, 0 = free motion

#constrain Y and Z coordinate of mass point to move z=10*y
mbs.CreateCoordinateConstraint(itemNumbers=[oMass1, oMass1], 
                               coordinates=[2,1],
                               factor1=10)

mbs.CreatePrismaticJoint(itemNumbers=[oGround, oBody], 
                         position=[3,0,0], #global position of joint
                         axis=[1,0,0], #global axis of joint, can move in global x-direction
                         useGlobalFrame=True) #use local coordinates for joint definition

mbs.CreateRevoluteJoint(itemNumbers=[oBody, oBody2], 
                        position=[4,-0.5,0], #global position of joint
                        axis=[0,0,1], #rotation along global z-axis
                        useGlobalFrame=True)

mbs.CreateTorsionalSpringDamper(itemNumbers=[oBody, oBody2],
                                position=[4,-0.5,0], #global position of spring-damper
                                axis=[0,0,1],        #global rotation axis
                                stiffness=1000,
                                damping=20,
                                useGlobalFrame=True)

mass = 5 
mu = 0.05 #static friction
rDisc = 0.5
length = 0.1
volume = np.pi * rDisc**2 * length
inertiaCylinder = InertiaCylinder(density=mass/volume, length=length, outerRadius=rDisc, axis=0)

oDisc = mbs.CreateRigidBody(inertia=inertiaCylinder,
                            nodeType=exu.NodeType.RotationRxyz, #this allows to constrain Z-rotation
                            referencePosition=[0,2,rDisc-1], #reference position x/y/z of COM
                            initialVelocity=[0,2*np.pi*rDisc,0],   
                            initialAngularVelocity=[-2*np.pi,0.2,0],
                            gravity=[0,0,-9.81])

mbs.CreateRollingDiscPenalty(bodyNumbers=[oGround, oDisc], discRadius = rDisc, 
                             axisPosition=[0,0,0], axisVector=[1,0,0],
                             planePosition = [0,0,-1], planeNormal = [0,0,1],
                             contactStiffness = 1e4, contactDamping = 2e2, 
                             dryFriction = [mu,mu], color=graphics.color.red)

#constrain Z and Y-rotation of disc, so it rotates only around X
mbs.CreateCoordinateConstraint(itemNumbers=[oDisc, None], coordinates=[5,None])
mbs.CreateCoordinateConstraint(itemNumbers=[oDisc, oGround], coordinates=[4,None])#also works with ground

mass = 5 
rDisc = 0.5
length = 0.1
volume = np.pi * rDisc**2 * length
inertiaCylinder = InertiaCylinder(density=mass/volume, length=length, outerRadius=rDisc, axis=0)

#create a free rigid body using the defined inertia properties and applies gravity (y-direction).
oDisc1 = mbs.CreateRigidBody(inertia=inertiaCylinder,
                            referencePosition=[1,2,rDisc-1], #reference position x/y/z of COM
                            initialVelocity=[0,2*np.pi*rDisc,0],   
                            initialAngularVelocity=[-2*np.pi,0.2,0],
                            gravity=[0,0,-9.81])

mbs.CreateForce(itemNumber=oDisc1, loadVector=[10,0,0])

#create a 'rolling' joint between flat ground defined by plane, lying on oGround, and rigid body given as oDisc:
mbs.CreateRollingDisc(bodyNumbers=[oGround, oDisc1], discRadius = rDisc, 
                      axisPosition=[0,0,0], axisVector=[1,0,0], #relative to oDisc
                      planePosition = [0,0,-1], planeNormal = [0,0,1], #defines plane
                      constrainedAxes=[1,1,1] #constrain 3 axes: lateral motion, forward motion and normal contact
                      )

#constrain the two discs to have same X-motion
mbs.CreateCoordinateConstraint(itemNumbers=[oDisc,oDisc1],
                               velocityLevel=True,
                               coordinates=[0,0])

mbs.CreateCartesianSpringDamper(
    itemNumbers=[oGround, oBody2], #[body0,body1]
    localPosition0=[4,-0.5*rbX,0], #for body0
    localPosition1=[0,0,0],        #for body1
    stiffness = [100,10,10], #x,y,z stiffness
    damping = [5,2,2], #x,y,z damping
)

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

endTime = 0.5 #simulation time in seconds
if not testIsActive:
    endTime = 2.5

stepSize = 0.002 #should be small enough to achieve sufficient accuracy

#some simulation parameters:
simulationSettings = exu.SimulationSettings() #takes currently set values or default values
simulationSettings.timeIntegration.numberOfSteps = int(endTime/stepSize)
simulationSettings.timeIntegration.endTime = endTime
simulationSettings.timeIntegration.verboseMode = 1
simulationSettings.timeIntegration.newton.useModifiedNewton = True

#for redundant constraints, use these two settings
simulationSettings.linearSolver.ignoreSingularJacobian=True
simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenDense #use EigenSparse for larger systems alternatively
#simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
# simulationSettings.timeIntegration.newton.useModifiedNewton = True
SC.visualizationSettings.openGL.multiSampling = 4

mbs.SolveDynamic(simulationSettings)


#%%
#TAG: visualization
#DIFFICULTY: 10000
SC.visualizationSettings.nodes.drawNodesAsPoint = False
SC.visualizationSettings.nodes.showBasis = True

if not testIsActive:
    mbs.SolutionViewer()

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
ode2 = mbs.systemData.GetODE2Coordinates()

u = 0.01*np.linalg.norm(ode2)
exu.Print('solution of createFunctionsTest=',u) 

exu.sys['testResult'] = u
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
```