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