NGsolveOCCboundaries2.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Test for NGsolve interface with fem, using OCC Extrude() and Revolve()
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2026-02-02
  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 *
 16import exudyn.graphics as graphics
 17from exudyn.FEM import HCBstaticModeSelection, FEMinterface, ObjectFFRFreducedOrderInterface
 18
 19SC = exu.SystemContainer()
 20mbs = SC.AddSystem()
 21
 22import numpy as np
 23
 24import time
 25
 26import netgen.occ as occ
 27from ngsolve import Mesh, Draw
 28
 29
 30width = 3
 31height = 2
 32length = 4
 33thickness = 0.08
 34arc = 0.5
 35bDim = 0.1
 36hBox = 0.8*height
 37arc2 = 0.4
 38
 39#create frame
 40wp = occ.WorkPlane(occ.Axes(p=(0,0,0), n=occ.X, h=occ.Y)) #axis=x, plane: (h=Y, v=Z)
 41#wp.MoveTo(0,0).Line(0.5*width*2).Rotate(90).Line(height).Rotate(45).Line(0.5).Close()
 42wp.MoveTo(0,0).Line(0.5*width-arc).Arc(arc,90).Line(height-arc).Rotate(90).Line(thickness).\
 43        Rotate(90).Line(height-arc).Arc(arc-thickness,-90).Line(width-2*arc).Arc(arc-thickness,-90).\
 44        Line(height-arc).Rotate(90).Line(thickness).Rotate(90).Line(height-arc).\
 45        Arc(arc,90).Close()
 46
 47frame = wp.Face().Extrude(length)
 48interfaceNameList = ['ground']
 49frame.faces.Min((0,0,1)).name = 'ground'
 50
 51#create tube
 52rTube = (0.5*width-thickness-bDim*0.5)
 53rTubeInner = 0.5*rTube
 54lTube = length
 55wp1 = occ.WorkPlane(occ.Axes(p=(0,0,0), n=occ.Y, h=occ.Z)) #axis=X, plane: (h=X, v=Z)
 56
 57#MoveTo(h,v): move cursor to position in workplane
 58#LineTo(h,v): create wire to to position in workplane
 59#Line(dx,dy): create wire in direction given
 60#Arc( radius, angleDegree): create arc in current direction with radius and angle
 61#ArcTo:  #destination h, v and destination tangent
 62#Rotate(angleDegree): change current direction relative to previous one (affects Line() and Arc() )
 63wp1.MoveTo( rTube, 0 )\
 64   .LineTo( rTube, lTube)\
 65   .LineTo( rTube-thickness, lTube)\
 66   .LineTo( rTube-thickness, thickness)\
 67   .Line( -thickness*1.5,0).Rotate(180)\
 68   .Arc( arc2, -90)\
 69   .LineTo( rTubeInner+thickness, lTube-thickness)\
 70   .LineTo( rTube-thickness*1.5, lTube-thickness)\
 71   .LineTo( rTube-thickness*1.5, lTube)\
 72   .LineTo( rTubeInner+thickness, lTube)\
 73   .LineTo( rTubeInner, lTube)\
 74   .LineTo( rTubeInner, arc2 )\
 75   .ArcTo(rTubeInner+arc2, 0, (0,-1))\
 76   .Close()
 77
 78
 79axis1 = occ.Axis((0,0,0),occ.X)
 80tube = wp1.Face().Revolve(axis1,360).Move((0,0,hBox))
 81
 82#boxes between frame and tube:
 83vBox = [2*0.5*bDim,2*0.5*bDim,0.5*bDim]
 84
 85boxList = []
 86for fy in [-1,1]:
 87    for fx in [0,0.5,1]:
 88        pBox = np.array((0.5*2*bDim+(length-2*bDim)*fx,fy*(-0.5*width+thickness+vBox[1]),hBox))
 89        box = occ.Box(tuple(pBox-vBox),
 90                      tuple(pBox+vBox))
 91        boxList.append(box)
 92
 93geo = occ.OCCGeometry(frame+tube+boxList[0]+boxList[1]+boxList[2]+boxList[3]+boxList[4]+boxList[5])
 94print('meshing ...')
 95maxh=15
 96
 97if False:
 98    import netgen.gui #this starts netgen gui; Press button "Visual" and activate "Auto-redraw after (sec)"; Then select "Mesh"
 99
100
101geoMesh = geo.GenerateMesh(maxh=maxh,
102                           curvaturesafety=2,
103                           )
104
105mesh = Mesh(geoMesh)
106
107
108gFloor = graphics.CheckerBoard(point=[1,0,0.],size=8)
109oGround = mbs.CreateGround(graphicsDataList=[gFloor])
110
111if False: #just show mesh:
112    [points, triangles, normals] = graphics.NGsolveMesh2PointsAndTrigs(mesh=mesh,
113                                                                       scale=1,
114                                                                       meshOrder=2,
115                                                                       addNormals=True,
116                                                                       )
117    color = graphics.color.steelblue
118    meshColor=graphics.color.lawngreen[0:3]+[graphics.material.indexChrome]
119
120    gMesh = graphics.FromPointsAndTrigs(points, triangles, normals=normals,
121                                        color=meshColor)
122    mbs.CreateGround(graphicsDataList=[gMesh])
123    SC.visualizationSettings.view0.scene.showFaceEdges = True
124
125
126rho = 2800
127Emodulus = 8e11
128nu = 0.3
129femInterface = FEMinterface()
130[bfM, bfK, fes] = femInterface.ImportMeshFromNGsolve(mesh,
131                                                     density=rho, youngsModulus=Emodulus, poissonsRatio=nu,
132                                                     boundaryNamesList=interfaceNameList,
133                                                     meshOrder=2
134                                                     )
135
136[boundaryNodesList, boundaryWeightsList] = femInterface.GetBoundaryNodeSetsAsLists()
137femInterface.ComputeHurtyCraigBamptonModes(boundaryNodesList=boundaryNodesList,
138                                           nEigenModes=24,
139                                           excludeRigidBodyMotion=True,
140                                           boundaryNodesWeights=boundaryWeightsList,
141                                           computationMode=HCBstaticModeSelection.RBE2)
142
143print('eigenfrequencies (Hz):\n',femInterface.GetEigenFrequenciesHz(),sep='')
144
145#invert boundary faces (BUG in surface creation in FEM?)
146trigs = femInterface.surface[0]['Trigs']
147femInterface.surface[0]['Trigs'][:,[0,2,1]] = trigs
148
149cms = ObjectFFRFreducedOrderInterface(femInterface)
150
151objFFRF = cms.AddObjectFFRFreducedOrder(mbs, positionRef=[0,0,0],
152                                              initialVelocity=[0,0,0],
153                                              initialAngularVelocity=[0,0,0],
154                                              color=[0.9,0.9,0.9,1.],
155                                              )
156
157mbs.Assemble()
158
159
160#%%
161if True: #activate to animate modes
162    from exudyn.interactive import AnimateModes
163
164    SC.visualizationSettings.nodes.show = False
165    SC.visualizationSettings.view0.scene.showFaceEdges = False
166    SC.visualizationSettings.openGL.multiSampling = 4
167    SC.visualizationSettings.openGL.light0.shadow = 0.2
168    SC.visualizationSettings.openGL.light0.position=[2,2,10,1]
169    SC.visualizationSettings.openGL.light0.diffuse = 0.4
170    SC.visualizationSettings.openGL.lineWidth=2
171    SC.visualizationSettings.openGL.lightModelAmbient = [0.6,0.6,0.6,1]
172    SC.visualizationSettings.view0.window.renderWindowSize = [1400,1024]
173    #SC.visualizationSettings.view0.camera.clippingPlaneNormal = [0,-1,0]
174
175
176    nodeNumber = objFFRF['nGenericODE2'] #this is the node with the generalized coordinates
177
178    SC.renderer.Start()              #start graphics visualization
179    SC.renderer.SetModelView(zoom=2.011357,rotationVector=[-0.9999199,0.6224784,0.9588339],centerPoint=[0.7874073,1.334914,0])
180    AnimateModes(SC, mbs, nodeNumber, period=0.1,
181                 showTime=False, scaleAmplitude=20,
182                 runOnStart=True)
183    import sys
184    sys.exit()
185
186
187SC.visualizationSettings.view0.window.renderWindowSize=[1200,800]
188SC.visualizationSettings.general.autoFitScene=False
189
190#SC.visualizationSettings.view0.scene.drawCoordinateSystem = False
191SC.visualizationSettings.openGL.lineWidth = 2
192#SC.visualizationSettings.openGL.light0.position = [-2.0, 4.0, 1.0, 1.0]
193SC.visualizationSettings.loads.show = False
194SC.visualizationSettings.openGL.light0.position=[2,2,10,1]
195
196#raytracing options
197SC.visualizationSettings.openGL.multiSampling = 2
198SC.visualizationSettings.openGL.light0.shadow = 0.2
199SC.visualizationSettings.openGL.light1.shadow = 0.2
200SC.visualizationSettings.raytracer.numberOfThreads = 64
201
202SC.visualizationSettings.view0.camera.useRaytracer = False #set True for raytracing
203SC.visualizationSettings.raytracer.keepWindowActive= True
204SC.visualizationSettings.raytracer.advanced.searchTreeFactor = 8
205
206
207#visualize in Exudyn:
208SC.renderer.Start()              #start graphics visualization
209
210SC.renderer.SetModelView(zoom=2.011357,rotationVector=[-0.9999199,0.6224784,0.9588339],centerPoint=[0.7874073,1.334914,0])
211SC.renderer.DoIdleTasks() #press space to continue
212
213SC.renderer.Stop() #safely close rendering window!