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