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])