You can view and download this file on Github: NGsolvePostProcessingStresses.py
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This is an EXUDYN example
#
# Details: Test for meshing with NETGEN and import of FEM model;
# Model is a simple flexible pendulum meshed with tet elements;
# Note that the model is overly flexible (linearized strain assumption not valid),
# but it should serve as a demonstration of the FFRFreducedOrder modeling
#
# Author: Johannes Gerstmayr
# Date: 2021-02-05
#
# 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.itemInterface import *
from exudyn.utilities import *
from exudyn.FEM import *
from exudyn.graphicsDataUtilities import *
SC = exu.SystemContainer()
mbs = SC.AddSystem()
import numpy as np
import timeit
import exudyn.basicUtilities as eb
import exudyn.rigidBodyUtilities as rb
import exudyn.utilities as eu
import numpy as np
useGraphics = True
fileName = 'testData/netgenBrick' #for load/save of FEM data
if __name__ == '__main__': #needed to use multiprocessing for mode computation
#+++++++++++++++++++++++++++++++++++++++++++++++++++++
#netgen/meshing part:
fem = FEMinterface()
#standard:
a = 0.025 #height/width of beam
b = a
h = 0.3*a
L = 1 #Length of beam
nModes = 10
#plate:
# a = 0.025 #height/width of beam
# b = 0.4
# L = 1 #Length of beam
# h = 0.6*a
# nModes = 40
rho = 1000
Emodulus=1e7
nu=0.3
meshCreated = False
#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
if True: #needs netgen/ngsolve to be installed to compute mesh, see e.g.: https://github.com/NGSolve/ngsolve/releases
import ngsolve as ngs
from netgen.geom2d import unit_square
import netgen.libngpy as libng
from netgen.csg import *
geo = CSGeometry()
block = OrthoBrick(Pnt(0,-a,-b),Pnt(L,a,b))
geo.Add(block)
#Draw (geo)
mesh = ngs.Mesh( geo.GenerateMesh(maxh=h))
mesh.Curve(1)
if False: #set this to true, if you want to visualize the mesh inside netgen/ngsolve
import netgen.gui
Draw (mesh)
netgen.Redraw()
#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
#Use fem to import FEM model and create FFRFreducedOrder object
fem.ImportMeshFromNGsolve(mesh, density=rho, youngsModulus=Emodulus, poissonsRatio=nu)
meshCreated = True
if (h==a): #save only if it has smaller size
fem.SaveToFile(fileName)
#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
if not meshCreated: fem.LoadFromFile(fileName)
print("nNodes=",fem.NumberOfNodes())
fem.ComputeEigenmodes(nModes, excludeRigidBodyModes = 6, useSparseSolver = True)
#print("eigen freq.=", fem.GetEigenFrequenciesHz())
#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
#compute stress modes:
mat = KirchhoffMaterial(Emodulus, nu, rho)
varType = exu.OutputVariableType.StressLocal
#varType = exu.OutputVariableType.StrainLocal
print("ComputePostProcessingModes ... (may take a while)")
start_time = time.time()
fem.ComputePostProcessingModes(material=mat,
outputVariableType=varType,
numberOfThreads=5)
print("--- %s seconds ---" % (time.time() - start_time))
#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
print("create CMS element ...")
cms = ObjectFFRFreducedOrderInterface(fem)
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.SetObjectParameter(objectNumber=objFFRF['oFFRFreducedOrder'],
# parameterName='outputVariableModeBasis',
# value=stressModes)
# mbs.SetObjectParameter(objectNumber=objFFRF['oFFRFreducedOrder'],
# parameterName='outputVariableTypeModeBasis',
# value=exu.OutputVariableType.StressLocal) #type=stress modes ...
#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 markers and joints
nodeDrawSize = 0.0025 #for joint drawing
pLeft = [0,-a,-b]
pRight = [0,-a,b]
nTip = fem.GetNodeAtPoint([L,-a,-b]) #tip node
#print("nMid=",nMid)
mRB = mbs.AddMarker(MarkerNodeRigid(nodeNumber=objFFRF['nRigidBody']))
oGround = mbs.AddObject(ObjectGround(referencePosition= [0,0,0]))
mGroundPosLeft = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition=pLeft))
mGroundPosRight = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition=pRight))
#++++++++++++++++++++++++++++++++++++++++++
#find nodes at left and right surface:
#nodeListLeft = fem.GetNodesInPlane(pLeft, [0,0,1])
#nodeListRight = fem.GetNodesInPlane(pRight, [0,0,1])
nodeListLeft = [fem.GetNodeAtPoint(pLeft)]
nodeListRight = [fem.GetNodeAtPoint(pRight)]
lenLeft = len(nodeListLeft)
lenRight = len(nodeListRight)
weightsLeft = np.array((1./lenLeft)*np.ones(lenLeft))
weightsRight = np.array((1./lenRight)*np.ones(lenRight))
addSupports = True
if addSupports:
k = 10e8 #joint stiffness
d = k*0.01 #joint damping
useSpringDamper = True
mLeft = mbs.AddMarker(MarkerSuperElementPosition(bodyNumber=objFFRF['oFFRFreducedOrder'],
meshNodeNumbers=np.array(nodeListLeft), #these are the meshNodeNumbers
weightingFactors=weightsLeft))
mRight = mbs.AddMarker(MarkerSuperElementPosition(bodyNumber=objFFRF['oFFRFreducedOrder'],
meshNodeNumbers=np.array(nodeListRight), #these are the meshNodeNumbers
weightingFactors=weightsRight))
if useSpringDamper:
oSJleft = mbs.AddObject(CartesianSpringDamper(markerNumbers=[mLeft, mGroundPosLeft],
stiffness=[k,k,k], damping=[d,d,d]))
oSJright = mbs.AddObject(CartesianSpringDamper(markerNumbers=[mRight,mGroundPosRight],
stiffness=[k,k,0], damping=[d,d,d]))
else:
oSJleft = mbs.AddObject(SphericalJoint(markerNumbers=[mGroundPosLeft,mLeft], visualization=VObjectJointSpherical(jointRadius=nodeDrawSize)))
oSJright= mbs.AddObject(SphericalJoint(markerNumbers=[mGroundPosRight,mRight], visualization=VObjectJointSpherical(jointRadius=nodeDrawSize)))
#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
fileDir = 'solution/'
mbs.AddSensor(SensorSuperElement(bodyNumber=objFFRF['oFFRFreducedOrder'], meshNodeNumber=nTip, #meshnode number!
fileName=fileDir+'nMidDisplacementCMS'+str(nModes)+'Test.txt',
outputVariableType = exu.OutputVariableType.Displacement))
mbs.Assemble()
simulationSettings = exu.SimulationSettings()
SC.visualizationSettings.nodes.defaultSize = nodeDrawSize
SC.visualizationSettings.nodes.drawNodesAsPoint = False
SC.visualizationSettings.connectors.defaultSize = 2*nodeDrawSize
SC.visualizationSettings.nodes.show = False
SC.visualizationSettings.nodes.showBasis = True #of rigid body node of reference frame
SC.visualizationSettings.nodes.basisSize = 0.12
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 = 0.01
SC.visualizationSettings.markers.drawSimplified = False
SC.visualizationSettings.markers.show = False
SC.visualizationSettings.markers.defaultSize = 0.01
SC.visualizationSettings.loads.drawSimplified = False
# SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.DisplacementLocal
# SC.visualizationSettings.contour.outputVariableComponent = 0 #x-component
SC.visualizationSettings.contour.reduceRange=False
SC.visualizationSettings.contour.outputVariable = varType
SC.visualizationSettings.contour.outputVariableComponent = 0 #x-component
simulationSettings.solutionSettings.solutionInformation = "ObjectFFRFreducedOrder test"
h=0.25e-3
tEnd = 0.05
simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
simulationSettings.timeIntegration.endTime = tEnd
simulationSettings.solutionSettings.writeSolutionToFile = False
simulationSettings.timeIntegration.verboseMode = 1
#simulationSettings.timeIntegration.verboseModeFile = 3
simulationSettings.timeIntegration.newton.useModifiedNewton = True
simulationSettings.solutionSettings.sensorsWritePeriod = h
simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.8 #SHOULD work with 0.9 as well
#simulationSettings.displayStatistics = True
#simulationSettings.displayComputationTime = True
#create animation:
# simulationSettings.solutionSettings.recordImagesInterval = 0.005
# SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
SC.visualizationSettings.window.renderWindowSize=[1920,1080]
SC.visualizationSettings.openGL.multiSampling = 4
if True:
if useGraphics:
SC.visualizationSettings.general.autoFitScene=False
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 useGraphics:
SC.WaitForRenderEngineStopFlag()
exu.StopRenderer() #safely close rendering window!