SwitchingConstraintsTest.py

You can view and download this file on Github: SwitchingConstraintsTest.py

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  A revolute joint and two coordinate constraints that are **switched on and off during the
  5#           integration** through `activeConnector`, over 2 s in 4000 steps: the sum of the final
  6#           coordinates.
  7#           The model compares against a reference value written into it, so its
  8#           result is that difference and its reference solution is 0 (#2632).
  9#
 10# Author:   Johannes Gerstmayr
 11# Date:     2019-11-01, reworked 2026-09-24
 12#
 13# 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.
 14#
 15#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 16
 17import exudyn as exu
 18from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 19import exudyn.graphics as graphics
 20
 21testIsActive = exu.sys.get('testIsActive', False)
 22exu.sys['testTolerance'] = 4e-13 #the tolerance RunAllModelUnitTests used for these ten
 23
 24SC = exu.SystemContainer()
 25mbs = SC.AddSystem()
 26
 27
 28rect = [-2.5,-1.5,0.5,1.5] #xmin,ymin,xmax,ymax
 29background = {'type':'Line', 'color':[0.1,0.1,0.8,1], 'data':[rect[0],rect[1],0, rect[2],rect[1],0, rect[2],rect[3],0, rect[0],rect[3],0, rect[0],rect[1],0]} #background
 30oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0], visualization=VObjectGround(graphicsData= [background])))
 31#nGround=mbs.AddNode(NodePointGround(referencePosition= [0,0,0]))
 32a = 0.5     #x-dim of pendulum
 33b = 0.05    #y-dim of pendulum
 34massRigid = 12
 35mass = 2 #of additional mass
 36inertiaRigid = massRigid/12*(2*a)**2
 37omega0 = 4
 38
 39graphics2 = {'type':'Line', 'color':[0.1,0.1,0.8,1], 'data':[-a,-b,0, a,-b,0, a,b,0, -a,b,0, -a,-b,0]} #background
 40nRigid = mbs.AddNode(Rigid2D(referenceCoordinates=[-0.5,0,0], initialVelocities=[0,omega0*a,omega0]));
 41oRigid = mbs.AddObject(RigidBody2D(mass=massRigid, inertia=inertiaRigid,nodeNumber=nRigid,visualization=VObjectRigidBody2D(graphicsData= [graphics2])))
 42
 43mR1 = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid, localPosition=[-0.5,0.,0.])) #support point
 44mG0 = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition=[-0.5-a,0.,0.]))
 45oRJoint = mbs.AddObject(RevoluteJoint2D(markerNumbers=[mG0,mR1],activeConnector=True))
 46
 47#mass point is attached with coordinate constraints:
 48mCoordR0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nRigid, coordinate=0)) 
 49mCoordR1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nRigid, coordinate=1)) 
 50
 51#additional mass point attached to COM of rigid body:
 52nMass = mbs.AddNode(Point2D(referenceCoordinates=[-0.5,0], initialVelocities=[0,omega0*a]));
 53oMass = mbs.AddObject(MassPoint2D(mass=mass, nodeNumber=nMass) )
 54
 55mCoordM0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMass, coordinate=0)) 
 56mCoordM1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMass, coordinate=1)) 
 57
 58oConstraint0 = mbs.AddObject(CoordinateConstraint(markerNumbers=[mCoordM0, mCoordR0],activeConnector=True))
 59oConstraint1 = mbs.AddObject(CoordinateConstraint(markerNumbers=[mCoordM1, mCoordR1],activeConnector=True))
 60
 61mbs.Assemble()
 62
 63simulationSettings = exu.SimulationSettings() #takes currently set values or default values
 64simulationSettings.timeIntegration.numberOfSteps = 4000
 65simulationSettings.timeIntegration.endTime = 2
 66simulationSettings.timeIntegration.newton.relativeTolerance = 1e-8
 67simulationSettings.timeIntegration.newton.absoluteTolerance = 1e-6
 68#simulationSettings.timeIntegration.verboseMode = 1
 69
 70#execute this python code in the same scope as this file at beginning of every time step:
 71#2019-12-18: subtract deltaT=1e-8 in order to avoid round off effects in 1e-16 regime
 72def UFswitchConnector(mbs, t):
 73    if t > (0.3 + 1e-8): 
 74        mbs.SetObjectParameter(oRJoint, 'activeConnector', False)
 75    if t > (0.1 + 1e-8): 
 76        mbs.SetObjectParameter(oConstraint0, 'activeConnector', False)
 77        mbs.SetObjectParameter(oConstraint1, 'activeConnector', False)
 78    return True #True, means that everything is alright, False=stop simulation
 79
 80mbs.SetPreStepUserFunction(UFswitchConnector)
 81
 82#simulationSettings.timeIntegration.verboseMode = 1
 83simulationSettings.timeIntegration.newton.useModifiedNewton = False
 84simulationSettings.timeIntegration.newton.numericalDifferentiation.minimumCoordinateSize = 1
 85#simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1
 86simulationSettings.timeIntegration.generalizedAlpha.useNewmark = True
 87simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = True
 88simulationSettings.solution.file.information = "Rigid pendulum with switching constraints"
 89simulationSettings.show.statistics = False
 90simulationSettings.solution.file.write=False
 91
 92#(not testIsActive) = True
 93if not testIsActive: 
 94    SC.renderer.Start()
 95
 96mbs.SolveDynamic(simulationSettings)
 97
 98if not testIsActive: 
 99    SC.renderer.DoIdleTasks()
100    SC.renderer.Stop() #safely close rendering window!
101
102u = mbs.systemData.GetODE2Coordinates()
103exu.Print('u =',u)
104exu.Print('sum(u) =',sum(u))
105
106testResult = abs(sum(u)-8.376384072333597) #2020-01-09: 376384072333597; before 13.12.2019: 8.342236349593959
107exu.Print('solution of SwitchingConstraintsTest=', testResult)
108exu.sys['testResult'] = testResult