linearFEMgenericODE2.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Import Abaqus FEM mesh and perform linear finite element analysis, for demonstration and example purpose;
  5#           example can be used to model linear FEM bodies directly in exudyn (may be slow; consider ObjectFFRFreducedOrder instead!!!)
  6#           2 cases considered: 1) use K and M directly; 2) use M and implement jacobian manually (to show how this would look like!)
  7#
  8# Author:   Johannes Gerstmayr, Andreas Zwölfer
  9# Date:     2024-10-06
 10#
 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.
 12#
 13#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 14
 15import exudyn as exu
 16from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 17import exudyn.graphics as graphics
 18from exudyn.FEM import *
 19
 20testIsActive = exu.sys.get('testIsActive', False)
 21#testIsActive = True #without test
 22
 23
 24import numpy as np
 25import time
 26import sys
 27import copy
 28
 29SC = exu.SystemContainer()
 30mbs = SC.AddSystem() 
 31
 32#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
 33#Use FEMinterface to import FEM model and create FFRFreducedOrder object
 34fileDir = 'testData/abaqus/block'
 35
 36fem = FEMinterface()
 37inputFileName = fileDir+'C3D8'
 38nodes=fem.ImportFromAbaqusInputFile(inputFileName+'.inp', typeName='Instance', name='rotor-1')
 39    
 40fem.ReadMassMatrixFromAbaqus(inputFileName+'_MASS1.mtx')
 41fem.ReadStiffnessMatrixFromAbaqus(inputFileName+'_STIF1.mtx')
 42
 43#scalse stiffness matrix; just done here as example to get larger deformations!!!
 44fem.ScaleStiffnessMatrix(1e-3) 
 45
 46nodePositions = fem.GetNodePositionsAsArray()
 47if False:
 48    exu.Print('number of nodes:', len(nodePositions))
 49    exu.Print('number of elements:', len(fem.elements[0]['Hex8']) )
 50    exu.Print('non-zero elements of massMatrix:', len(fem.massMatrix) )
 51    exu.Print('non-zero elements of stiffnessMatrix:', len(fem.stiffnessMatrix) )
 52    exu.Print('===================')
 53    
 54pLeft = [0,0,0]
 55pRight = [4,0,0]
 56nTip = fem.GetNodeAtPoint(pRight) #tip node (do not use midpoint, as this may not be a mesh node ...)
 57
 58nodesLeftPlane = fem.GetNodesInPlane(pLeft, [-1,0,0])
 59weightsLeftPlane = fem.GetNodeWeightsFromSurfaceAreas(nodesLeftPlane)
 60
 61nodesRightPlane = fem.GetNodesInPlane(pRight, [-1,0,0])
 62weightsRightPlane = fem.GetNodeWeightsFromSurfaceAreas(nodesRightPlane)
 63
 64boundaryList = [nodesLeftPlane] 
 65
 66exu.Print('nodesLeftPlane',nodesLeftPlane)
 67exu.Print('nodesRightPlane',nodesRightPlane)
 68
 69 
 70#compute the according sparse matrices
 71nRows = fem.NumberOfCoordinates()
 72nCols = nRows
 73nNodes = fem.NumberOfNodes()
 74
 75Mfem = fem.GetMassMatrix(sparse=True)
 76Kfem = fem.GetStiffnessMatrix(sparse=True)
 77#for test purposes, use stiffness-proportional damping
 78Dfem = 5e-3*Kfem #this is the damping factor alpha, using D=alpha*K
 79
 80#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 81
 82#create copy of matrix, only used for jacobian, as we need to multiply with fODE2 factor:
 83jac = exu.MatrixContainer()
 84
 85#add nodes to mbs model (if you have several bodies, you have to do this for each body):
 86allNodeList = [] #create node list
 87for node in nodePositions:
 88    allNodeList += [mbs.AddNode(NodePoint(referenceCoordinates=node))]
 89 
 90
 91#create user function; use Scipy's csr matrices for fast multiplication!
 92#as these terms are on the right-hand-side, the sign is negative
 93def UFforce(mbs, t, itemNumber, c, c_t):
 94    return -Kfem @ c - Dfem @ c_t
 95
 96#compute the jacobian; this is for the left-hand-side and therefore has positive K and D
 97def UFjacobian(mbs, t, itemNumber, c, c_t, fODE2, fODE2_t):
 98    #KDfem needs to include factors fODE2 and fODE2_t for jacobian
 99    #note that adding both components of Kfem and Dfem only works here because they have identical sparse format!
100    #OLD (exudyn < 1.8.69): KDfem[:,2] = fODE2 * Kfem[:,2] + fODE2_t * Dfem[:,2]
101    KDfem= fODE2 * Kfem + fODE2_t * Dfem
102    jac.SetWithSparseMatrix(KDfem) 
103    return jac
104    
105testUserFunction = False #set False to switcht to (faster) case without user function
106
107#now add generic body  
108if testUserFunction:
109    oGenericODE2 = mbs.AddObject(ObjectGenericODE2(nodeNumbers = allNodeList, 
110                                                    massMatrix=exu.MatrixContainer(Mfem), 
111                                                    forceUserFunction=UFforce,
112                                                    jacobianUserFunction=UFjacobian,
113                                                    visualization=VObjectGenericODE2(triangleMesh = fem.GetSurfaceTriangles(), color=graphics.color.lightred)
114                                                    ))
115else:
116    oGenericODE2 = mbs.AddObject(ObjectGenericODE2(nodeNumbers = allNodeList, 
117                                                   massMatrix=exu.MatrixContainer(Mfem), 
118                                                   stiffnessMatrix=exu.MatrixContainer(Kfem),
119                                                   dampingMatrix=exu.MatrixContainer(Dfem),
120                                                   visualization=VObjectGenericODE2(triangleMesh = fem.GetSurfaceTriangles(), color=graphics.color.lightred)
121                                                   ))
122 
123
124#%%++++++++++++++++++++++++++++++++++++++++
125#add constraints:
126oGround = mbs.AddObject(ObjectGround()) 
127mGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=[0,0.5,0.5]))
128 
129mBeam = mbs.AddMarker(MarkerSuperElementRigid(bodyNumber=oGenericODE2,
130                                              meshNodeNumbers=np.array(nodesLeftPlane), 
131                                              weightingFactors=weightsLeftPlane))
132 
133
134cGroundBeam = mbs.AddObject(GenericJoint(markerNumbers=[mGround, mBeam], constrainedAxes=[1,1,1,1,1,1]))
135     
136
137  
138#apply forces to right end, using all nodes equally:
139for node in nodesRightPlane:
140    mNode = mbs.AddMarker(MarkerNodePosition(nodeNumber=node))
141    mbs.AddLoad(Force(markerNumber=mNode, loadVector=[0,-2e4,0]))
142
143 
144 
145#add sensor to view results
146sTip = mbs.AddSensor(SensorNode(nodeNumber=nTip, storeInternal=True,
147                                outputVariableType=exu.OutputVariableType.Displacement))
148
149
150
151
152#%%++++++++++++++++++++++++++++++++++++++++
153mbs.Assemble()
154
155simulationSettings = exu.SimulationSettings()
156
157SC.visualizationSettings.nodes.defaultSize = 0.02
158SC.visualizationSettings.nodes.drawNodesAsPoint = False
159SC.visualizationSettings.loads.drawSimplified = False
160SC.visualizationSettings.openGL.multiSampling=4
161SC.visualizationSettings.openGL.lineWidth=2
162    
163h=1e-2 #default: 5e-4
164tEnd = 0.1
165if not testIsActive:
166    tEnd = 10
167
168simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
169simulationSettings.timeIntegration.endTime = tEnd
170simulationSettings.solution.file.write = (not testIsActive)
171simulationSettings.solution.file.writePeriod = 0.04
172simulationSettings.timeIntegration.verboseMode = 1
173# simulationSettings.timeIntegration.stepInformation = 255#8192-1
174#simulationSettings.timeIntegration.verboseModeFile = 3
175simulationSettings.timeIntegration.newton.useModifiedNewton = True
176simulationSettings.solution.sensors.writePeriod = h
177simulationSettings.timeIntegration.generalizedAlpha.computeInitialAccelerations = False
178#simulationSettings.show.statistics = True
179simulationSettings.show.computationTime = True
180simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
181
182 
183if not testIsActive:
184    SC.renderer.Start()
185    #SC.renderer.DoIdleTasks() #press space to continueq
186
187mbs.SolveDynamic(simulationSettings=simulationSettings)
188sensorValues = mbs.GetSensorValues(sTip)
189exu.Print('uTip=',list(sensorValues))
190#regular:       uTip= [-0.07561723475265847, -0.42046930823607603, 0.0001934540937925318]
191#user function: uTip= [-0.07561723475266069, -0.42046930823609246, 0.00019345409379260147]
192
193result = np.linalg.norm(sensorValues)
194exu.Print('solution of linearFEMgenericODE2=',result)
195
196exu.sys['testResult'] = result
197
198
199if not testIsActive:
200    #SC.renderer.DoIdleTasks()
201    SC.renderer.Stop() #safely close rendering window!
202
203    mbs.SolutionViewer()
204
205    # plot results:
206    if True: 
207        from exudyn.plot import PlotSensor
208        PlotSensor(mbs, sTip, components=[0,1,2])