CMSexampleCourse.py

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This is an EXUDYN example
#
# Details:  Create flexible multibody system
#
# Author:   Johannes Gerstmayr 
# Date:     2021-09-10
# Python version: Python 3.7, 64bits, Anaconda3  + ngsolve + webgui_jupyter_widgets
# Jupyter:  requires upgrade of scipy and uninstall and install tk (tkinter)
# 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 #only import if it does not conflict
from exudyn.FEM import *
import numpy as np
from math import sqrt, sin, cos, pi

doMeshing = True #set false if mesh shall be loaded
useHCBmodes = True #Hurty-Craig-Bampton modes
computeStresses = False #takes some time
loadStresses = False #set True only if already computed previously; may lead to severe problems if wrong modes are loaded!!!!

fileName = 'testData/FMBStest1' #for load/save of FEM data


# # Parameter Definition
# ![example geometry](exampleBody.png "Geometry")

#flexible body dimensions:
femInterface = FEMinterface()
#geometrical parameters
ri = 0.010 #radius of hole/bolt
ra = 0.016 #outer radius
t  = 0.020 #part thickness, bolt length
f  = 0.004 #radius of part rounding
f2 = 0.002 #radius of bolt rounding
L = 0.250  #distance hole/bolt

hp = ra*sqrt(0.5)-f*(1-sqrt(0.5)) #height of bar
d1 = (ra+1*f)*sqrt(0.5) #length of rounded part
L1 = L-2*d1 #length of rectangular section

#number of eigenmodes
nModes = 8

rho = 1000
Emodulus=1e8
nu=0.3


# # Create FEM mesh in Netgen

if doMeshing: #needs netgen/ngsolve to be installed to compute mesh, see e.g.: https://github.com/NGSolve/ngsolve/releases

    from netgen.occ import *
    #from ngsolve.webgui import Draw #in Jupyter
    from ngsolve import Mesh, Draw

    #create part geometry
    wp = WorkPlane(gp_Ax3(p=(0,0,-1.5*t), n=Z, h=X))

    wp.Rotate(90).MoveTo(-ra,0).Arc(ra,-135).Arc(f,45).Line(0.5*L1).Line(0.5*L1).Arc(f,45).Arc(ra,-135)
    wp.Arc(ra,-135).Arc(f,45).Line(L1).Arc(f,45).Arc(ra,-135)
    wp.Close().Reverse()

    p1 = wp.Face().Extrude(t)
    
    #bolt:
    wp2 = WorkPlane(Axes(p=(0,0,-0.5*t), n=X, h=Y))
    f22 = np.sqrt(2)*f2
    #wp2.MoveTo(0,0).LineTo(ri,0).LineTo(ri,t).LineTo(0,t).LineTo(0,0) #without rounding
    wp2.MoveTo(0,0).Line(ri+f2).Rotate(180).Arc(f2,-90).Line(t-2*f2).Rotate(45).Line(f22).Rotate(45).Line(ri-f2).Rotate(90).Line(t).Close() #with rounding+chamfer
    axis2 = Axis((0,0,0),Z)
    p2 = wp2.Face().Revolve(axis2,360)

    #hole:
    wp3 = WorkPlane(Axes(p=(L,0,-1.5*t), n=X, h=Y))
    #wp2.MoveTo(0,0).LineTo(ri,0).LineTo(ri,t).LineTo(0,t).LineTo(0,0) #without rounding
    wp3.MoveTo(0,0).Line(ri+f2).Rotate(135).Line(f22).Rotate(-45).Line(t-2*f2).Rotate(-45).Line(f22).Rotate(135).Line(ri+f2).Rotate(90).Line(t).Close() #with rounding
    axis3 = Axis((L,0,0),Z)
    p3 = wp3.Face().Revolve(axis3,360)
    
    p1 = p1 + p2
    p1 = p1 - p3

    #for geometry check:
    #     box = Box((0,0,0), (ri,ri,ri)) #show (0,0,0)
    #     p1 = p1+ box
    geo = OCCGeometry( p1 )

    #Jupyter, webgui, draw geometry
    #NEEDS: pip install webgui_jupyter_widgets
    from netgen.webgui import Draw as DrawGeo
    #DrawGeo(geo.shape)

    #generate mesh:
    from ngsolve.webgui import Draw
    mesh = Mesh(geo.GenerateMesh(maxh=1.5*f2))
    #Jupyter, webgui, draw mesh
    Draw(mesh)



# # Import mesh into Exudyn
SC = exu.SystemContainer()
mbs = SC.AddSystem()

if doMeshing: #needs netgen/ngsolve to be installed to compute mesh, see e.g.: https://github.com/NGSolve/ngsolve/releases
    #save FEM mesh
    femInterface.ImportMeshFromNGsolve(mesh, density=rho, youngsModulus=Emodulus, poissonsRatio=nu)
    femInterface.SaveToFile(fileName)
else:
    femInterface.LoadFromFile(fileName)
print("number of nodes = ", femInterface.NumberOfNodes())


# In[12]:


femInterface.ComputeEigenmodes(nModes, excludeRigidBodyModes = 6, useSparseSolver = True)
if False: #activate to animate modes
    from exudyn.interactive import AnimateModes
    mbs.Reset()
    cms = ObjectFFRFreducedOrderInterface(femInterface)

    objFFRF = cms.AddObjectFFRFreducedOrder(mbs, positionRef=[0,0,0], 
                                                  initialVelocity=[0,0,0], 
                                                  initialAngularVelocity=[0,0,0],
                                                  color=[0.1,0.9,0.1,1.],
                                                  )
    mbs.Assemble()
    SC.visualizationSettings.nodes.show = False
    SC.visualizationSettings.openGL.showFaceEdges = True
    SC.visualizationSettings.openGL.multiSampling=4
    SC.visualizationSettings.openGL.lineWidth=2
    SC.visualizationSettings.window.renderWindowSize = [1600,1080]

    #%%+++++++++++++++++++++++++++++++++++++++
    SC.visualizationSettings.general.autoFitScene = False #otherwise, model may be difficult to be moved

    nodeNumber = objFFRF['nGenericODE2'] #this is the node with the generalized coordinates
    AnimateModes(SC, mbs, nodeNumber, period=0.1, 
                 scaleAmplitude = 0.02,
                 showTime=False, renderWindowText='Show modes\n',
                 runOnStart=True)


# # Define interfaces
addSensors = True
pLeft = [0,0,0] #midpoint of bolt
pRight = [L,0,-t] #midpoint of hole
pMid = [0.5*L,hp,-0.5*t] #midpoint of bar
pTip = [L+ra,0,-0.5*t] #midpoint of bar

#%%  
if addSensors:
    nMid = femInterface.GetNodeAtPoint(pMid, tolerance=1e-2) #tip node (do not use midpoint, as this may not be a mesh node ...)
    print("pMid=",pMid,", nMid=",nMid)
    nTip = femInterface.GetNodeAtPoint(pTip, tolerance=1e-2) #tip node (do not use midpoint, as this may not be a mesh node ...)
    print("pTip=",pTip,", nTip=",nTip)

tV = np.array([0,0,0.5*t])
nodesLeft = femInterface.GetNodesOnCylinder(pLeft-tV, pLeft+tV, ri)
# print('nodesLeft=',nodesLeft)
nodesRight = femInterface.GetNodesOnCylinder(pRight-tV, pRight+tV, ri)
# print('nodesRight=',nodesRight)

lenNodesLeft = len(nodesLeft)
weightsNodesLeft = np.array((1./lenNodesLeft)*np.ones(lenNodesLeft))

lenNodesRight = len(nodesRight)
weightsNodesRight = np.array((1./lenNodesRight)*np.ones(lenNodesRight))

boundaryList = [nodesLeft, nodesRight] #second boudary (right plane) not needed ...


# # Compute eigenmodes
#remark: ComputeEigenmodes requires upgrade of scipy (python -m pip install --upgrade scipy) as compared to Anaconda installation...
import time

print("compute modes... ")
start_time = time.time()

if useHCBmodes:
    femInterface.ComputeHurtyCraigBamptonModes(boundaryNodesList=boundaryList, 
                                  nEigenModes=nModes, 
                                  useSparseSolver=True,
                                  computationMode = HCBstaticModeSelection.RBE2)
else:
    femInterface.ComputeEigenmodes(nModes, excludeRigidBodyModes = 6, useSparseSolver = True)

print("computation of modes needed %.3f seconds" % (time.time() - start_time))
print("eigen freq.=", femInterface.GetEigenFrequenciesHz())




# # Compute stresses

femModesName = fileName+'modes'
if useHCBmodes:
    femModesName+='HCB'

varType = exu.OutputVariableType.StressLocal

if computeStresses:
    mat = KirchhoffMaterial(Emodulus, nu, rho)
    #varType = exu.OutputVariableType.StrainLocal
    print("ComputePostProcessingModes ... (may take a while)")
    start_time = time.time()
    femInterface.ComputePostProcessingModes(material=mat, 
                                   outputVariableType=varType)
    print("   ... needed %.3f seconds" % (time.time() - start_time))
    SC.visualizationSettings.contour.reduceRange=False
    SC.visualizationSettings.contour.outputVariable = varType
    SC.visualizationSettings.contour.outputVariableComponent = 0 #x-component
    #save modes + stresses
    femInterface.SaveToFile(femModesName)
else:
    if loadStresses:
        femInterface.LoadFromFile(femModesName)
        SC.visualizationSettings.contour.outputVariable = varType
        SC.visualizationSettings.contour.outputVariableComponent = 0 #x-component


# # Setup flexible body in exudyn
cms = ObjectFFRFreducedOrderInterface(femInterface)

objFFRF = cms.AddObjectFFRFreducedOrder(mbs, positionRef=[0,0,0], 
                                              initialVelocity=[0,0,0], 
                                              initialAngularVelocity=[0,0,0],
                                              color=[0.1,0.9,0.1,1.],
                                              )


# # Visualize modes
if False:
    from exudyn.interactive import AnimateModes
    mbs.Assemble()
    SC.visualizationSettings.nodes.show = False
    SC.visualizationSettings.openGL.showFaceEdges = True
    SC.visualizationSettings.openGL.multiSampling=4
    #SC.visualizationSettings.window.renderWindowSize = [1600,1080]
    SC.visualizationSettings.general.autoFitScene = False #otherwise, model may be difficult to be moved

    nodeNumber = objFFRF['nGenericODE2'] #this is the node with the generalized coordinates
    AnimateModes(SC, mbs, nodeNumber, scaleAmplitude = 0.1, runOnStart = True)


# # Add gravity

# In[11]:


#add gravity (not necessary if user functions used)
oFFRF = objFFRF['oFFRFreducedOrder']
mBody = mbs.AddMarker(MarkerBodyMass(bodyNumber=oFFRF))
mbs.AddLoad(LoadMassProportional(markerNumber=mBody, loadVector= [0,-9.81,0]))


# # Add joint constraint

# In[12]:


#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
#add markers and joints

#mRB = mbs.AddMarker(MarkerNodeRigid(nodeNumber=objFFRF['nRigidBody']))
oGround = mbs.AddObject(ObjectGround(referencePosition = [0,0,0]))
mGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber = oGround, 
                                        localPosition = pLeft))

#add marker
#NOTE: offset is added in order to compensate for small errors in average node position
#      because mesh is not fully symmetric and the average node position does not match 
#      the desired (body-fixed) joint position [0,0,0]; if not used, a small initial jump
#      happens during simulation when the body moves to the constrained position
mLeft = mbs.AddMarker(MarkerSuperElementRigid(bodyNumber=objFFRF['oFFRFreducedOrder'], 
                                              meshNodeNumbers=np.array(nodesLeft), #these are the meshNodeNumbers
                                              weightingFactors=weightsNodesLeft,
                                              offset=-femInterface.GetNodePositionsMean(nodesLeft),
                                             ))
oJoint = mbs.AddObject(GenericJoint(markerNumbers=[mGround, mLeft], 
                            constrainedAxes = [1,1,1,1,1,0],
                            visualization=VGenericJoint(axesRadius=0.05*ri, axesLength=1.5*t)))
# oJoint = mbs.AddObject(RevoluteJointZ(markerNumbers=[mGround, mLeft],
#                               visualization=VRevoluteJointZ(axisRadius=0.05*ri, axisLength=1.5*t)))

if False: #if this is used, remove offset in MarkerSuperElementRigid above
    #alternative to offset above: compensate joint offset by computation of current displacement in joint: (if not done in MarkerSuperElementRigid)
    mbs.Assemble() #initialize system to compute joint offset
    jointOffset = mbs.GetObjectOutput(oJoint,exu.OutputVariableType.DisplacementLocal)
    print('jointOffset=',jointOffset)

    mbs.SetMarkerParameter(mLeft.GetIndex(), 'offset', list(-jointOffset)) #compensate offset; mLeft.GetIndex() because of BUG751

    #now check new offset:
    mbs.Assemble() #initialize system to compute joint offset
    jointOffset = mbs.GetObjectOutput(oJoint,exu.OutputVariableType.DisplacementLocal)
    print('jointOffset=',jointOffset)


# # Add sensors
fileDir = 'solution/'
if addSensors:
    sMidDispl = mbs.AddSensor(SensorSuperElement(bodyNumber=objFFRF['oFFRFreducedOrder'], 
                              meshNodeNumber=nMid, #meshnode number!
                              fileName=fileDir+'uMid'+str(nModes)+'modes.txt', 
                              outputVariableType = exu.OutputVariableType.Displacement))
    sTipDispl = mbs.AddSensor(SensorSuperElement(bodyNumber=objFFRF['oFFRFreducedOrder'], 
                              meshNodeNumber=nTip, #meshnode number!
                              fileName=fileDir+'uTip'+str(nModes)+'modes.txt', 
                              outputVariableType = exu.OutputVariableType.Displacement))


# # Set up visualization
# (not needed)
nodeDrawSize = 0.0025 #for joint drawing
SC.visualizationSettings.nodes.defaultSize = nodeDrawSize
SC.visualizationSettings.nodes.drawNodesAsPoint = False
SC.visualizationSettings.connectors.defaultSize = nodeDrawSize

SC.visualizationSettings.nodes.show = False
SC.visualizationSettings.nodes.showBasis = True #of rigid body node of reference frame
SC.visualizationSettings.nodes.basisSize = t*4
SC.visualizationSettings.bodies.deformationScaleFactor = 1 #use this factor to scale the deformation of modes

SC.visualizationSettings.openGL.showFaceEdges = True
SC.visualizationSettings.openGL.showFaces = True

SC.visualizationSettings.sensors.show = True
SC.visualizationSettings.sensors.drawSimplified = False
SC.visualizationSettings.sensors.defaultSize = nodeDrawSize*2
SC.visualizationSettings.markers.drawSimplified = False
SC.visualizationSettings.markers.show = False
SC.visualizationSettings.markers.defaultSize = nodeDrawSize*2

SC.visualizationSettings.loads.drawSimplified = False
SC.visualizationSettings.loads.defaultSize = t*3
SC.visualizationSettings.loads.defaultRadius = 0.05*t

SC.visualizationSettings.window.renderWindowSize=[1280,720]
SC.visualizationSettings.openGL.multiSampling = 4

#create animation:
# simulationSettings.solutionSettings.recordImagesInterval = 0.005
# SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"


# # Set up simulation
mbs.Assemble() #initialize bodies, assemble system; necessary to simulate

simulationSettings = exu.SimulationSettings()
simulationSettings.solutionSettings.solutionInformation = "ObjectFFRFreducedOrder test"

h=1e-3
tEnd = 2

simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
simulationSettings.timeIntegration.endTime = tEnd
simulationSettings.solutionSettings.writeSolutionToFile = True
simulationSettings.solutionSettings.solutionWritePeriod = h
simulationSettings.timeIntegration.verboseMode = 1
#simulationSettings.timeIntegration.verboseModeFile = 3
simulationSettings.timeIntegration.newton.useModifiedNewton = True

simulationSettings.solutionSettings.sensorsWritePeriod = h

simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.8
#simulationSettings.displayStatistics = True
simulationSettings.displayComputationTime = True


# # Start renderer and Simulate

lifeVisualization = True

if lifeVisualization:
    SC.visualizationSettings.general.autoFitScene=False #if reloaded view settings
    exu.StartRenderer()
    if 'renderState' in exu.sys: SC.SetRenderState(exu.sys['renderState']) #load last model view
    mbs.WaitForUserToContinue() #press space to continue
        
mbs.SolveDynamic(#solverType=exu.DynamicSolverType.TrapezoidalIndex2, 
                  simulationSettings=simulationSettings)
            
if addSensors:
    uTip = mbs.GetSensorValues(sMidDispl)
    print("nModes=", nModes, ", mid displacement=", uTip)

if lifeVisualization:
    SC.WaitForRenderEngineStopFlag()
    exu.StopRenderer() #safely close rendering window!


# # 3D rendering of FMBS
if False: #use this to reload the solution and use SolutionViewer
    SC.visualizationSettings.general.autoFitScene=False #if reloaded view settings

    
    mbs.SolutionViewer() #can also be entered in IPython ...


# # Plot sensor

if addSensors:
    
    mbs.PlotSensor(sensorNumbers=[sMidDispl,sMidDispl,sMidDispl], components=[0,1,2])