NGsolveLinearFEM.py
You can view and download this file on Github: NGsolveLinearFEM.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Linear FEM model using NGsolve and ObjectGenericODE2
5#
6# Author: Johannes Gerstmayr, Joachim Schöberl
7# Date: 2021-10-05
8#
9# 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.
10#
11#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
12
13
14import exudyn as exu
15from exudyn.itemInterface import *
16from exudyn.utilities import *
17from exudyn.FEM import *
18from exudyn.graphicsDataUtilities import *
19
20SC = exu.SystemContainer()
21mbs = SC.AddSystem()
22
23import numpy as np
24import sys
25import time
26
27
28#import netgen.geom2d as geom2d
29from netgen.occ import *
30import ngsolve as ngs
31
32# from ngsolve.webgui import Draw
33# from netgen.webgui import Draw as DrawGeo
34
35#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
36# define geometry and mesh
37L = 1
38wy = 0.1
39wz = 0.12
40body = Box((0,0,0), (L,wy, wz))
41#body.bc("all")
42
43faces = body.SubShapes(FACE)
44faces[0].bc("left")
45faces[0].col=(1,0,0)
46
47geo = OCCGeometry(body)
48mesh = ngs.Mesh(geo.GenerateMesh(maxh=0.05*1)) #0.05*0.25 gives quite fine mesh (13GB)
49#DrawGeo(geo.shape)
50#Draw(mesh)
51#print(mesh.dim)
52
53#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
54# define material parameters and energy
55meshOrder = 1
56youngsModulus = 210
57nu = 0.2
58mu = youngsModulus / 2 / (1+nu)
59lam = youngsModulus * nu / ((1+nu)*(1-2*nu))
60density = 1
61
62
63fem=FEMinterface()
64fem.ImportMeshFromNGsolve(mesh, density, youngsModulus, nu, meshOrder=meshOrder)
65
66#%%++++++++++++++++++++++++++++++++++++++++
67[oGenericODE2, allNodeList] = fem.CreateLinearFEMObjectGenericODE2(mbs, color=color4dodgerblue)
68
69#%%++++++++++++++++++++++++++++++++++++++++
70#add forces on right side and fix on left side:
71nLists = 2
72nodeLists = [[]]*nLists
73nNodes = [0]*nLists
74
75nodeLists[0] = fem.GetNodesInPlane(point=[0,0,0], normal=[1,0,0])
76nodeLists[1] = fem.GetNodesInPlane(point=[L,0,0], normal=[1,0,0])
77
78for i in range(nLists):
79 nNodes[i] = len(nodeLists[i])
80
81#apply force to right end:
82fLoad = 1/nNodes[1] * np.array([0,-1e-3,0])
83for i in nodeLists[1]:
84 mNode = mbs.AddMarker(MarkerNodePosition(nodeNumber=i))
85 mbs.AddLoad(Force(markerNumber=mNode, loadVector=fLoad))
86
87oGround = mbs.AddObject(ObjectGround())
88
89if False:
90 #apply single sphereical constraints to left end:
91 for i in nodeLists[0]:
92 mNode = mbs.AddMarker(MarkerNodePosition(nodeNumber=i))
93 mGroundI = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround,
94 localPosition=fem.GetNodePositionsAsArray()[i]))
95 mbs.AddObject(ObjectJointSpherical(markerNumbers = [mNode, mGroundI],
96 # constrainedAxes=[0,0,0],
97 visualization=VSphericalJoint(jointRadius=0.015)))
98else: #use superelement marker
99 #pMid = [0,wy*0.5,wz*0.5]
100 pMid = fem.GetNodePositionsMean(nodeLists[0])
101 mGroundI = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround,
102 localPosition=pMid))
103 mLeft = mbs.AddMarker(MarkerSuperElementRigid(bodyNumber=oGenericODE2,
104 meshNodeNumbers=nodeLists[0],
105 useAlternativeApproach=False,
106 weightingFactors=[1/nNodes[0]]*nNodes[0],
107 offset = [0,0,0]))
108 #mbs.AddObject(ObjectJointSpherical(markerNumbers = [mLeft, mGroundI],
109 # visualization=VSphericalJoint(jointRadius=0.015)))
110 mbs.AddObject(GenericJoint(markerNumbers = [mLeft, mGroundI],
111 visualization=VGenericJoint(axesRadius=0.015, axesLength=0.02)))
112
113
114
115
116#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
117
118mbs.Assemble()
119
120simulationSettings = exu.SimulationSettings()
121
122nodeDrawSize = 0.01
123
124SC.visualizationSettings.nodes.defaultSize = nodeDrawSize
125SC.visualizationSettings.nodes.drawNodesAsPoint = False
126SC.visualizationSettings.connectors.defaultSize = 1.25*nodeDrawSize
127
128SC.visualizationSettings.nodes.show = False
129SC.visualizationSettings.nodes.showBasis = False #of rigid body node of reference frame
130SC.visualizationSettings.nodes.basisSize = 0.12
131SC.visualizationSettings.bodies.deformationScaleFactor = 1 #use this factor to scale the deformation of modes
132
133SC.visualizationSettings.openGL.showFaceEdges = True
134SC.visualizationSettings.openGL.showFaces = True
135
136SC.visualizationSettings.sensors.show = True
137SC.visualizationSettings.sensors.drawSimplified = False
138SC.visualizationSettings.sensors.defaultSize = 0.01
139
140SC.visualizationSettings.markers.show = True
141SC.visualizationSettings.markers.defaultSize=1.2*nodeDrawSize
142SC.visualizationSettings.markers.drawSimplified = False
143
144SC.visualizationSettings.loads.show = False
145SC.visualizationSettings.loads.drawSimplified = False
146SC.visualizationSettings.loads.defaultSize=0.1
147SC.visualizationSettings.loads.defaultRadius = 0.002
148
149SC.visualizationSettings.openGL.multiSampling=4
150SC.visualizationSettings.openGL.lineWidth=2
151
152h=1e-3*0.5
153tEnd = 2
154
155simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
156simulationSettings.timeIntegration.endTime = tEnd
157simulationSettings.solutionSettings.writeSolutionToFile = False
158simulationSettings.timeIntegration.verboseMode = 1
159#simulationSettings.timeIntegration.verboseModeFile = 3
160simulationSettings.timeIntegration.newton.useModifiedNewton = True
161
162simulationSettings.solutionSettings.sensorsWritePeriod = h
163
164simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.7
165#simulationSettings.displayStatistics = True
166simulationSettings.displayComputationTime = True
167simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
168#create animation:
169# simulationSettings.solutionSettings.recordImagesInterval = 0.005
170# SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
171SC.visualizationSettings.window.renderWindowSize=[1920,1080]
172SC.visualizationSettings.openGL.multiSampling = 4
173# SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.Displacement
174# SC.visualizationSettings.contour.outputVariableComponent = 1 #y-component
175
176useGraphics=True
177if True:
178 if useGraphics:
179 SC.visualizationSettings.general.autoFitScene=False
180
181 exu.StartRenderer()
182 if 'renderState' in exu.sys: SC.SetRenderState(exu.sys['renderState']) #load last model view
183
184 mbs.WaitForUserToContinue() #press space to continue
185
186 #SC.RedrawAndSaveImage()
187 if True:
188 # mbs.SolveDynamic(solverType=exu.DynamicSolverType.TrapezoidalIndex2,
189 # simulationSettings=simulationSettings)
190 mbs.SolveDynamic(simulationSettings=simulationSettings)
191 else:
192 mbs.SolveStatic(simulationSettings=simulationSettings)
193
194 # uTip = mbs.GetSensorValues(sensTipDispl)[1]
195 # print("nModes=", nModes, ", tip displacement=", uTip)
196
197 if useGraphics:
198 SC.WaitForRenderEngineStopFlag()
199 exu.StopRenderer() #safely close rendering window!
200
201 if False:
202
203 mbs.PlotSensor(sensorNumbers=[sensBushingVel], components=[1])