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