NGsolveOCCboundaries.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Test for NGsolve interface with fem, using OCC Extrude() and interface names
  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 = 1.5
 32length = 4
 33thickness = 0.1
 34arc = 0.5
 35maxh = 0.5
 36
 37#simple solid of revolution geometry:
 38wp = occ.WorkPlane(occ.Axes(p=(0,0,0), n=occ.X, h=occ.Y))
 39#wp.MoveTo(0,0).Line(0.5*width*2).Rotate(90).Line(height).Rotate(45).Line(0.5).Close()
 40wp.MoveTo(0,0).Line(0.5*width-arc).Arc(arc,90).Line(height-arc).Rotate(90).Line(thickness).\
 41        Rotate(90).Line(height-arc).Arc(arc-thickness,-90).Line(width-2*arc).Arc(arc-thickness,-90).\
 42        Line(height-arc).Rotate(90).Line(thickness).Rotate(90).Line(height-arc).\
 43        Arc(arc,90).Close()
 44
 45frame = wp.Face().Extrude(length)
 46
 47bDim = 0.1
 48vBox = [0.5*2*bDim,0.5*bDim,0.5*bDim]
 49
 50boxList = []
 51interfaceNameList = []
 52for fy in [-1,1]:
 53    for fx in [0,0.5,1]:
 54        name = 'box_y'+str(fy)+'_x'+str(fx)
 55        interfaceNameList.append(name)
 56        pBox = np.array((0.5*2*bDim+(length-2*bDim)*fx,fy*(-0.5*width+thickness+0.5*bDim),0.5*height))
 57        box = occ.Box(tuple(pBox-vBox),
 58                      tuple(pBox+vBox))
 59        if fy > 0:
 60            box.faces.Max((0, 1, 0)).name = name
 61        else:
 62            box.faces.Min((0, 1, 0)).name = name
 63        boxList.append(box)
 64
 65geo = occ.OCCGeometry(frame+boxList[0]+boxList[1]+boxList[2]+boxList[3]+boxList[4]+boxList[5])
 66exu.Print('meshing ...')
 67
 68if False:
 69    import netgen.gui #this starts netgen gui; Press button "Visual" and activate "Auto-redraw after (sec)"; Then select "Mesh"
 70
 71
 72geoMesh = geo.GenerateMesh(maxh=maxh,
 73                           curvaturesafety=1.5,
 74                           )
 75
 76mesh = Mesh(geoMesh)
 77
 78
 79gFloor = graphics.CheckerBoard(point=[0,0,0.],size=8)
 80oGround = mbs.CreateGround(graphicsDataList=[gFloor])
 81
 82#+++++++++++++++++++++++++++++++++++++++++++++
 83#create FEM
 84rho = 2800
 85Emodulus = 8e11
 86nu = 0.3
 87femInterface = FEMinterface()
 88[bfM, bfK, fes] = femInterface.ImportMeshFromNGsolve(mesh,
 89                                                     density=rho, youngsModulus=Emodulus, poissonsRatio=nu,
 90                                                     boundaryNamesList=interfaceNameList,
 91                                                     meshOrder=2
 92                                                     )
 93
 94[boundaryNodesList, boundaryWeightsList] = femInterface.GetBoundaryNodeSetsAsLists()
 95femInterface.ComputeHurtyCraigBamptonModes(boundaryNodesList=boundaryNodesList,
 96                                           nEigenModes=8,
 97                                           excludeRigidBodyMotion=True,
 98                                           boundaryNodesWeights=boundaryWeightsList,
 99                                           computationMode=HCBstaticModeSelection.RBE2)
100
101exu.Print('eigenfrequencies (Hz):\n',femInterface.GetEigenFrequenciesHz(),sep='')
102
103cms = ObjectFFRFreducedOrderInterface(femInterface)
104
105objFFRF = cms.AddObjectFFRFreducedOrder(mbs, positionRef=[0,0,0],
106                                              initialVelocity=[0,0,0],
107                                              initialAngularVelocity=[0,0,0],
108                                              color=[0.9,0.9,0.9,1.],
109                                              )
110
111mbs.Assemble()
112
113
114
115#+++++++++++++++++++++++++++++++++++++++++++++
116if True: #activate to animate modes
117    from exudyn.interactive import AnimateModes
118    SC.visualizationSettings.nodes.show = False
119    SC.visualizationSettings.view0.scene.showFaceEdges = True
120    SC.visualizationSettings.openGL.multiSampling=2
121    SC.visualizationSettings.openGL.lineWidth=2
122    SC.visualizationSettings.view0.window.renderWindowSize = [1600,1080]
123
124
125    nodeNumber = objFFRF['nGenericODE2'] #this is the node with the generalized coordinates
126
127    SC.renderer.Start()              #start graphics visualization
128    SC.renderer.SetModelView(zoom=2.011357,rotationVector=[-0.9999199,0.6224784,0.9588339],centerPoint=[0.7874073,1.334914,0])
129    AnimateModes(SC, mbs, nodeNumber, period=0.1, showTime=False, renderWindowText='Hurty-Craig-Bampton: 2 x 6 static modes and 8 eigenmodes\n',
130                 runOnStart=True)
131    import sys
132    sys.exit()
133
134#+++++++++++++++++++++++++++++++++++++++++++++
135
136SC.visualizationSettings.view0.window.renderWindowSize=[1200,800]
137SC.visualizationSettings.general.autoFitScene=False
138
139#SC.visualizationSettings.view0.scene.drawCoordinateSystem = False
140SC.visualizationSettings.openGL.lineWidth = 2
141#SC.visualizationSettings.openGL.light0.position = [-2.0, 4.0, 1.0, 1.0]
142SC.visualizationSettings.loads.show = False
143SC.visualizationSettings.openGL.light0.position=[2,2,10,1]
144
145#raytracing options
146SC.visualizationSettings.openGL.multiSampling = 2
147SC.visualizationSettings.openGL.light0.shadow = 0.2
148SC.visualizationSettings.openGL.light1.shadow = 0.2
149SC.visualizationSettings.raytracer.numberOfThreads = 64
150
151SC.visualizationSettings.view0.camera.useRaytracer = False #set True for raytracing
152SC.visualizationSettings.raytracer.keepWindowActive= True
153SC.visualizationSettings.raytracer.advanced.searchTreeFactor = 8
154
155
156#visualize in Exudyn:
157SC.renderer.Start()              #start graphics visualization
158
159SC.renderer.SetModelView(zoom=2.011357,rotationVector=[-0.9999199,0.6224784,0.9588339],centerPoint=[0.7874073,1.334914,0])
160SC.renderer.DoIdleTasks() #press space to continue
161
162SC.renderer.Stop() #safely close rendering window!