NGsolveGeometry.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Example to show how to create CSG-geometry in Netgen and import as STL into exudyn
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2021-04-20
  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.utilities import *
 16from exudyn.FEM import *
 17
 18SC = exu.SystemContainer()
 19mbs = SC.AddSystem()
 20
 21import numpy as np
 22import time
 23
 24#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
 25#netgen/meshing part:
 26
 27#geometrical parameters:
 28L = 0.2         #Length of body
 29sy = 0.04       #width of body  (Y)
 30sz = 0.03       #height of body (Z)
 31dy = 0.005      #additional distance Y
 32r = 0.015       #radius of bolt
 33R = 0.025       #outer radius
 34x = 0.002
 35meshH = 0.005   #0.01 is default, 0.002 gives 100000 nodes and is fairly converged;
 36curvaturesafety = 5
 37
 38axis = 0.15
 39rotBasis = np.eye(3 )
 40
 41#steel: (not needed for drawing ...)
 42rho = 7850
 43Emodulus=2.1e11
 44nu=0.3
 45
 46#test high flexibility
 47Emodulus=2e8
 48# nModes = 32
 49
 50
 51#helper function for cylinder with netgen
 52def CSGcylinder(p0,p1,r):
 53    v = VSub(p1,p0)
 54    v = Normalize(v)
 55    cyl = Cylinder(Pnt(p0[0],p0[1],p0[2]), Pnt(p1[0],p1[1],p1[2]),
 56                   r) * Plane(Pnt(p0[0],p0[1],p0[2]), Vec(-v[0],-v[1],-v[2])) * Plane(Pnt(p1[0],p1[1],p1[2]), Vec(v[0],v[1],v[2]))
 57    return cyl
 58
 59meshCreated = False
 60
 61#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
 62import ngsolve as ngs
 63import netgen
 64from netgen.meshing import *
 65
 66from netgen.geom2d import unit_square
 67#import netgen.libngpy as libng
 68from netgen.csg import *
 69
 70
 71showCase = 'revolute'
 72showCase = 'spheric'
 73
 74if showCase == 'revolute':
 75    rotBasis = RotationMatrixX(-0.5*pi)
 76    #++++++++++++++++++++++++++++++++++++++++++++++++++++
 77    #first body
 78    geo0 = CSGeometry()
 79
 80    #plate
 81    block = OrthoBrick(Pnt(0, x, -0.5*sz),Pnt(L, sy-x, 0.5*sz))
 82    blockCyl = CSGcylinder(p0=[0,0,0], p1=[0,sy,0], r=R)
 83    bolt0 = CSGcylinder(p0=[0,sy,0], p1=[0,2*sy+dy,0], r=r)
 84    geo0.Add(block+blockCyl+bolt0)
 85
 86    mesh0 = ngs.Mesh( geo0.GenerateMesh(maxh=meshH, curvaturesafety=curvaturesafety))
 87    mesh0.Curve(1)
 88
 89    #++++++++++++++++++++++++++++++++++++++++++++++++++++
 90    #second body
 91    geo1 = CSGeometry()
 92
 93    #plate
 94    block = OrthoBrick(Pnt(0, x, -0.5*sz),Pnt(L, sy-x, 0.5*sz))
 95    blockCyl = CSGcylinder(p0=[0,0,0], p1=[0,sy,0], r=R)
 96    hole = CSGcylinder(p0=[0,-sy*0.1,0], p1=[0,1.1*sy,0], r=r)
 97    geo1.Add(block+blockCyl-hole)
 98
 99    mesh1 = ngs.Mesh( geo1.GenerateMesh(maxh=meshH, curvaturesafety=curvaturesafety))
100    mesh1.Curve(1)
101
102
103    if False: #set this to true, if you want to visualize the mesh inside netgen/ngsolve
104        import netgen.gui
105        ngs.Draw(mesh1)
106        for i in range(10000000):
107            netgen.Redraw() #this makes the netgen window interactive
108            time.sleep(0.05)
109
110    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
111    #import body with fem (this could be simplified in future ...)
112    #body0
113    fem0 = FEMinterface()
114    fem0.ImportMeshFromNGsolve(mesh0, density=rho, youngsModulus=Emodulus, poissonsRatio=nu)
115    graphics0 = GraphicsDataFromPointsAndTrigs(fem0.GetNodePositionsAsArray(), fem0.GetSurfaceTriangles(),
116                                               color=color4dodgerblue, )
117    graphics0 = AddEdgesAndSmoothenNormals(graphics0)
118
119    mbs.CreateRigidBody(referencePosition=[0,-sy*2-dy,0],
120                        referenceRotationMatrix=RotationMatrixX(0),
121                        inertia=InertiaCuboid(1000, [1,1,1]),
122                        graphicsDataList=[graphics0])
123
124
125    #+++++++++++++++++++++++++++++++++++++++++++++++++++++
126    #body1
127    fem1 = FEMinterface()
128    fem1.ImportMeshFromNGsolve(mesh1, density=rho, youngsModulus=Emodulus, poissonsRatio=nu)
129    graphics1 = GraphicsDataFromPointsAndTrigs(fem1.GetNodePositionsAsArray(), fem1.GetSurfaceTriangles(),
130                                               color=color4lightred)
131    graphics1 = AddEdgesAndSmoothenNormals(graphics1)
132
133    mbs.CreateRigidBody(referencePosition=[0,-sy,0],
134                        referenceRotationMatrix=RotationMatrixY(1.2*pi),
135                        inertia=InertiaCuboid(1000, [1,1,1]),
136                        graphicsDataList=[graphics1])
137
138if showCase == 'spheric':
139    #++++++++++++++++++++++++++++++++++++++++++++++++++++
140    #first body
141    geo0 = CSGeometry()
142    meshH *= 0.5
143
144    block = OrthoBrick(Pnt(0, -0.4*sy, -0.4*sz),Pnt(L, 0.4*sy, 0.4*sz))
145    sphere = Sphere( Pnt(0, 0, 0), R*1.2)
146    cutBlock = OrthoBrick(Pnt(-1, -sy, -sz),Pnt(-0.3*R, sy, sz))
147    cutSphere = Sphere( Pnt(0, 0, 0), R)
148    geo0.Add(block+sphere-cutBlock-cutSphere)
149
150    mesh0 = ngs.Mesh( geo0.GenerateMesh(maxh=meshH, curvaturesafety=curvaturesafety))
151    mesh0.Curve(1)
152
153    #++++++++++++++++++++++++++++++++++++++++++++++++++++
154    #second body
155    geo1 = CSGeometry()
156
157    #plate
158    block = OrthoBrick(Pnt(0, -0.4*sz, -0.4*sz),Pnt(L, 0.4*sz, 0.4*sz))
159    sphere = Sphere( Pnt(0, 0, 0), R)
160    geo1.Add(block+sphere)
161
162    mesh1 = ngs.Mesh( geo1.GenerateMesh(maxh=meshH, curvaturesafety=curvaturesafety))
163    mesh1.Curve(1)
164
165
166    if False: #set this to true, if you want to visualize the mesh inside netgen/ngsolve
167        import netgen.gui
168        ngs.Draw(mesh1)
169        for i in range(10000000):
170            netgen.Redraw() #this makes the netgen window interactive
171            time.sleep(0.05)
172
173    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
174    #import body with fem (this could be simplified in future ...)
175    #body0
176    fem0 = FEMinterface()
177    fem0.ImportMeshFromNGsolve(mesh0, density=rho, youngsModulus=Emodulus, poissonsRatio=nu)
178    graphics0 = GraphicsDataFromPointsAndTrigs(fem0.GetNodePositionsAsArray(), fem0.GetSurfaceTriangles(),
179                                               color=color4dodgerblue, )
180    graphics0 = AddEdgesAndSmoothenNormals(graphics0)
181
182    mbs.CreateRigidBody(referencePosition=[0,-sy,0],
183                        referenceRotationMatrix=RotationMatrixX(0),
184                        inertia=InertiaCuboid(1000, [1,1,1]),
185                        graphicsDataList=[graphics0])
186
187
188    #+++++++++++++++++++++++++++++++++++++++++++++++++++++
189    #body1
190    fem1 = FEMinterface()
191    fem1.ImportMeshFromNGsolve(mesh1, density=rho, youngsModulus=Emodulus, poissonsRatio=nu)
192    graphics1 = GraphicsDataFromPointsAndTrigs(fem1.GetNodePositionsAsArray(), fem1.GetSurfaceTriangles(),
193                                               color=color4lightred)
194    graphics1 = AddEdgesAndSmoothenNormals(graphics1)
195
196    mbs.CreateRigidBody(referencePosition=[0,-sy,0],
197                        referenceRotationMatrix=RotationMatrixZ(-0.12*pi)@RotationMatrixY(1.15*pi),
198                        inertia=InertiaCuboid(1000, [1,1,1]),
199                        graphicsDataList=[graphics1])
200
201
202#+++++++++++++++++++++++++++++++++++++++++++++++++
203#world basis
204gl=[]
205gl += [GraphicsDataText(point=rotBasis@[axis,0,0],text='X')]
206gl += [GraphicsDataText(point=rotBasis@[0,axis,0],text='Y')]
207gl += [GraphicsDataText(point=rotBasis@[0,0,axis],text='Z')]
208
209gl += [GraphicsDataBasis(origin=[0,0,0], rotationMatrix=rotBasis, length=axis)]
210
211mbs.CreateGround(graphicsDataList=gl)
212
213mbs.Assemble()
214#%%++++++++++++++++++++++++++++++
215SC.visualizationSettings.openGL.polygonOffset = 0.1 #to draw edges clearly
216SC.visualizationSettings.openGL.lineWidth = 2
217SC.visualizationSettings.openGL.multiSampling = 4
218# SC.visualizationSettings.general.drawWorldBasis = True
219# SC.visualizationSettings.general.worldBasisSize = axis
220SC.visualizationSettings.general.drawCoordinateSystem = False
221SC.visualizationSettings.general.textSize = 16
222SC.visualizationSettings.window.renderWindowSize = [1600,1200]
223
224mbs.SolveDynamic()
225
226mbs.SolutionViewer()