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