You can view and download this file on Github: NGsolveGeometry.py
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This is an EXUDYN example
#
# Details: Example to show how to create CSG-geometry in Netgen and import as STL into exudyn
#
# Author: Johannes Gerstmayr
# Date: 2021-04-20
#
# 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.
#
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
import exudyn as exu
from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
import exudyn.graphics as graphics #only import if it does not conflict
from exudyn.FEM import *
SC = exu.SystemContainer()
mbs = SC.AddSystem()
import numpy as np
import time
#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
#netgen/meshing part:
#geometrical parameters:
L = 0.2 #Length of body
sy = 0.04 #width of body (Y)
sz = 0.03 #height of body (Z)
dy = 0.005 #additional distance Y
r = 0.015 #radius of bolt
R = 0.025 #outer radius
x = 0.002
meshH = 0.005 #0.01 is default, 0.002 gives 100000 nodes and is fairly converged;
curvaturesafety = 5
axis = 0.15
rotBasis = np.eye(3 )
#steel: (not needed for drawing ...)
rho = 7850
Emodulus=2.1e11
nu=0.3
#test high flexibility
Emodulus=2e8
# nModes = 32
#helper function for cylinder with netgen
def CSGcylinder(p0,p1,r):
v = VSub(p1,p0)
v = Normalize(v)
cyl = Cylinder(Pnt(p0[0],p0[1],p0[2]), Pnt(p1[0],p1[1],p1[2]),
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]))
return cyl
meshCreated = False
#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
import ngsolve as ngs
import netgen
from netgen.meshing import *
from netgen.geom2d import unit_square
#import netgen.libngpy as libng
from netgen.csg import *
showCase = 'revolute'
showCase = 'spheric'
if showCase == 'revolute':
rotBasis = RotationMatrixX(-0.5*pi)
#++++++++++++++++++++++++++++++++++++++++++++++++++++
#first body
geo0 = CSGeometry()
#plate
block = OrthoBrick(Pnt(0, x, -0.5*sz),Pnt(L, sy-x, 0.5*sz))
blockCyl = CSGcylinder(p0=[0,0,0], p1=[0,sy,0], r=R)
bolt0 = CSGcylinder(p0=[0,sy,0], p1=[0,2*sy+dy,0], r=r)
geo0.Add(block+blockCyl+bolt0)
mesh0 = ngs.Mesh( geo0.GenerateMesh(maxh=meshH, curvaturesafety=curvaturesafety))
mesh0.Curve(1)
#++++++++++++++++++++++++++++++++++++++++++++++++++++
#second body
geo1 = CSGeometry()
#plate
block = OrthoBrick(Pnt(0, x, -0.5*sz),Pnt(L, sy-x, 0.5*sz))
blockCyl = CSGcylinder(p0=[0,0,0], p1=[0,sy,0], r=R)
hole = CSGcylinder(p0=[0,-sy*0.1,0], p1=[0,1.1*sy,0], r=r)
geo1.Add(block+blockCyl-hole)
mesh1 = ngs.Mesh( geo1.GenerateMesh(maxh=meshH, curvaturesafety=curvaturesafety))
mesh1.Curve(1)
if False: #set this to true, if you want to visualize the mesh inside netgen/ngsolve
import netgen.gui
ngs.Draw(mesh1)
for i in range(10000000):
netgen.Redraw() #this makes the netgen window interactive
time.sleep(0.05)
#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
#import body with fem (this could be simplified in future ...)
#body0
fem0 = FEMinterface()
fem0.ImportMeshFromNGsolve(mesh0, density=rho, youngsModulus=Emodulus, poissonsRatio=nu)
graphics0 = graphics.FromPointsAndTrigs(fem0.GetNodePositionsAsArray(), fem0.GetSurfaceTriangles(),
color=graphics.color.dodgerblue, )
graphics0 = AddEdgesAndSmoothenNormals(graphics0)
mbs.CreateRigidBody(referencePosition=[0,-sy*2-dy,0],
referenceRotationMatrix=RotationMatrixX(0),
inertia=InertiaCuboid(1000, [1,1,1]),
graphicsDataList=[graphics0])
#+++++++++++++++++++++++++++++++++++++++++++++++++++++
#body1
fem1 = FEMinterface()
fem1.ImportMeshFromNGsolve(mesh1, density=rho, youngsModulus=Emodulus, poissonsRatio=nu)
graphics1 = graphics.FromPointsAndTrigs(fem1.GetNodePositionsAsArray(), fem1.GetSurfaceTriangles(),
color=graphics.color.lightred)
graphics1 = AddEdgesAndSmoothenNormals(graphics1)
mbs.CreateRigidBody(referencePosition=[0,-sy,0],
referenceRotationMatrix=RotationMatrixY(1.2*pi),
inertia=InertiaCuboid(1000, [1,1,1]),
graphicsDataList=[graphics1])
if showCase == 'spheric':
#++++++++++++++++++++++++++++++++++++++++++++++++++++
#first body
geo0 = CSGeometry()
meshH *= 0.5
block = OrthoBrick(Pnt(0, -0.4*sy, -0.4*sz),Pnt(L, 0.4*sy, 0.4*sz))
sphere = Sphere( Pnt(0, 0, 0), R*1.2)
cutBlock = OrthoBrick(Pnt(-1, -sy, -sz),Pnt(-0.3*R, sy, sz))
cutSphere = Sphere( Pnt(0, 0, 0), R)
geo0.Add(block+sphere-cutBlock-cutSphere)
mesh0 = ngs.Mesh( geo0.GenerateMesh(maxh=meshH, curvaturesafety=curvaturesafety))
mesh0.Curve(1)
#++++++++++++++++++++++++++++++++++++++++++++++++++++
#second body
geo1 = CSGeometry()
#plate
block = OrthoBrick(Pnt(0, -0.4*sz, -0.4*sz),Pnt(L, 0.4*sz, 0.4*sz))
sphere = Sphere( Pnt(0, 0, 0), R)
geo1.Add(block+sphere)
mesh1 = ngs.Mesh( geo1.GenerateMesh(maxh=meshH, curvaturesafety=curvaturesafety))
mesh1.Curve(1)
if False: #set this to true, if you want to visualize the mesh inside netgen/ngsolve
import netgen.gui
ngs.Draw(mesh1)
for i in range(10000000):
netgen.Redraw() #this makes the netgen window interactive
time.sleep(0.05)
#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
#import body with fem (this could be simplified in future ...)
#body0
fem0 = FEMinterface()
fem0.ImportMeshFromNGsolve(mesh0, density=rho, youngsModulus=Emodulus, poissonsRatio=nu)
graphics0 = graphics.FromPointsAndTrigs(fem0.GetNodePositionsAsArray(), fem0.GetSurfaceTriangles(),
color=graphics.color.dodgerblue, )
graphics0 = AddEdgesAndSmoothenNormals(graphics0)
mbs.CreateRigidBody(referencePosition=[0,-sy,0],
referenceRotationMatrix=RotationMatrixX(0),
inertia=InertiaCuboid(1000, [1,1,1]),
graphicsDataList=[graphics0])
#+++++++++++++++++++++++++++++++++++++++++++++++++++++
#body1
fem1 = FEMinterface()
fem1.ImportMeshFromNGsolve(mesh1, density=rho, youngsModulus=Emodulus, poissonsRatio=nu)
graphics1 = graphics.FromPointsAndTrigs(fem1.GetNodePositionsAsArray(), fem1.GetSurfaceTriangles(),
color=graphics.color.lightred)
graphics1 = AddEdgesAndSmoothenNormals(graphics1)
mbs.CreateRigidBody(referencePosition=[0,-sy,0],
referenceRotationMatrix=RotationMatrixZ(-0.12*pi)@RotationMatrixY(1.15*pi),
inertia=InertiaCuboid(1000, [1,1,1]),
graphicsDataList=[graphics1])
#+++++++++++++++++++++++++++++++++++++++++++++++++
#world basis
gl=[]
gl += [graphics.Text(point=rotBasis@[axis,0,0],text='X')]
gl += [graphics.Text(point=rotBasis@[0,axis,0],text='Y')]
gl += [graphics.Text(point=rotBasis@[0,0,axis],text='Z')]
gl += [graphics.Basis(origin=[0,0,0], rotationMatrix=rotBasis, length=axis)]
mbs.CreateGround(graphicsDataList=gl)
mbs.Assemble()
#%%++++++++++++++++++++++++++++++
SC.visualizationSettings.openGL.polygonOffset = 0.1 #to draw edges clearly
SC.visualizationSettings.openGL.lineWidth = 2
SC.visualizationSettings.openGL.multiSampling = 4
# SC.visualizationSettings.general.drawWorldBasis = True
# SC.visualizationSettings.general.worldBasisSize = axis
SC.visualizationSettings.general.drawCoordinateSystem = False
SC.visualizationSettings.general.textSize = 16
SC.visualizationSettings.window.renderWindowSize = [1600,1200]
mbs.SolveDynamic()
mbs.SolutionViewer()