<!-- GENERATED by tools/generators/examplesDocsEmitter.py from python/ - do not edit -->
(testmodels-ngsolvecmstest)=
# NGsolveCMStest.py

You can view and download this file on Github: [NGsolveCMStest.py](https://github.com/jgerstmayr/EXUDYN/tree/master/python/TestModels/NGsolveCMStest.py)

```{code-block} python
:linenos:

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This is an EXUDYN example
#
# Details:  Test for Hurty-Craig-Bampton modes using a simple flexible pendulum meshed with Netgen
#
# Author:   Johannes Gerstmayr 
# Date:     2021-04-20
# Update:   2022-07-11: runs now with pip installed ngsolve on Python 3.10
#
# 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
import numpy as np
from exudyn.FEM import *
import os   #the mesh file is checked for existence, not loaded speculatively (#2469)
import time

testIsActive = exu.sys.get('testIsActive', False)


SC = exu.SystemContainer()
mbs = SC.AddSystem()

#for load/save of FEM data. The tracked file is an .npz of plain numpy arrays since 2026-09-17
#(#2469): the .pkl it replaces carried exudyn C++ types, so loading it under a second module
#(exudynCPPfast) raised 'type "Real" is already registered' and the model then regenerated the
#mesh - changing its own reference value. The .npz was CONVERTED from that .pkl, so the mesh and
#the reference value are unchanged. (The PKL/HDF5 round trip further below still tests those
#formats; it writes into the output directory, not here.)
fileName = 'testData/netgenTestMesh'

#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
#netgen/meshing part:
fem = FEMinterface()
#standard:
a = 0.1     #half height/width of beam
b = a
h = 2*a     #much too coarse, only for testing!!! 
L = 1       #Length of beam
nModes = 4

rho = 1000
Emodulus=1e6 #smaller to see some deformation
nu=0.3
meshCreated = False
meshOrder = 1 #use order 2 for higher accuracy, but more unknowns

hasNGsolve = False

try:
    #netgen/ngsolve: https://github.com/NGSolve/ngsolve/releases
    import ngsolve as ngs
    from netgen.meshing import *
    from netgen.csg import *
    hasNGsolve = True
except:
    exu.Print('NGsolve not installed; trying to load mesh')

#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
if hasNGsolve: 
    
    geo = CSGeometry()
    
    block = OrthoBrick(Pnt(0,-a,-b),Pnt(L,a,b))
    geo.Add(block)
    
    mesh = ngs.Mesh( geo.GenerateMesh(maxh=h))
    mesh.Curve(1)

    if False: #set this to true, if you want to visualize the mesh inside netgen/ngsolve
        # import netgen
        import netgen.gui
        ngs.Draw (mesh)
        for i in range(10000000):
            netgen.Redraw() #this makes the window interactive
            time.sleep(0.05)

    meshCreated = True
    #%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
    #Use fem to import FEM model and create FFRFreducedOrder object
    [bfM, bfK, fes] = fem.ImportMeshFromNGsolve(mesh, density=rho, 
                                                youngsModulus=Emodulus, poissonsRatio=nu,
                                                meshOrder=meshOrder)
    #if file does not exist, create it - otherwise don't change it!
    #if you want to replace it, delete the old file!
    #
    #DECIDE ON THE FILE, NOT ON AN EXCEPTION (#2469): this used to save whenever the LOAD failed,
    #and a load fails for reasons that have nothing to do with the file being absent - under a
    #second C++ module it raises 'type "Real" is already registered', and a new ngsolve version
    #would do the same. The mesh file is TRACKED and the reference value belongs to it, so an
    #overwrite silently moved the result by 2.4e-8 and left the working tree dirty.
    if not os.path.isfile(fileName + '.npz'):
        exu.Print('\nNGsolveCMStest: ' + fileName + '.npz does not exist; creating it from the '
                  'mesh just generated\n')
        fem.SaveToFile(fileName, mode='NPZ')

#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
#compute Hurty-Craig-Bampton modes
#the file is the reference: it is never rewritten here, so a failure to load has to be reported
#with its reason rather than repaired silently (#2469)
if not os.path.isfile(fileName + '.npz'):
    raise ValueError('NGsolveCMStest: ' + fileName + '.npz not found and ngsolve is not '
                     'available to create it')
try:
    fem.LoadFromFile(fileName, mode='NPZ')
except Exception as e:
    exu.Print(f'\nNGsolveCMStest: LoadFromFile(...) failed; {e}\n')
    raise ValueError('NGsolveCMStest: ' + fileName + '.npz exists but cannot be loaded (see the '
                     'message above); the file is tracked and holds the reference solution, so it '
                     'is NOT regenerated. Delete it deliberately to have it rebuilt.')

    
pLeft = [0,-a,-b]
pRight = [L,-a,-b]
nTip = fem.GetNodeAtPoint(pRight) #tip node (do not use midpoint, as this may not be a mesh node ...)

nodesLeftPlane = fem.GetNodesInPlane(pLeft, [-1,0,0])
weightsLeftPlane = fem.GetNodeWeightsFromSurfaceAreas(nodesLeftPlane)

boundaryList = [nodesLeftPlane] 

fem.ComputeHurtyCraigBamptonModes(boundaryNodesList=boundaryList, 
                              nEigenModes=nModes, 
                              useSparseSolver=True,
                              computationMode = HCBstaticModeSelection.RBE2)


#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
#compute stress modes for postprocessing (inaccurate for coarse meshes, just for visualization):
mat = KirchhoffMaterial(Emodulus, nu, rho)
varType = exu.OutputVariableType.StressLocal
start_time = time.time()
fem.ComputePostProcessingModes(material=mat, outputVariableType=varType)
SC.visualizationSettings.contour.reduceRange=False
SC.visualizationSettings.contour.outputVariable = varType
SC.visualizationSettings.contour.outputVariableComponent = 0 #x-component

#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
#now we test load/save for different formats!
#the files of this round trip are OUTPUT of the test: write them into
#exudyn.config.outputDirectory, which the test suite sets per model, so that two models cannot
#write the same file (#2418); the mesh loaded further above is INPUT and is read as given
from exudyn.basicUtilities import OutputFilePath
fileName2 = OutputFilePath(fileName+'2', 'NGsolveCMStest')
try:
    fem.SaveToFile(fileName2+'2.npz')
    fem.LoadFromFile(fileName2+'2', mode='NPZ')
except:
    #this should always work!
    raise ValueError('NGsolveCMStest: load/save with NPZ mode failed!')

try:
    fem.SaveToFile(fileName2+'.pkl', mode='PKL')
    fem.LoadFromFile(fileName2, mode='PKL')
except:
    #this should always work!
    raise ValueError('NGsolveCMStest: load/save with PKL mode failed!')

hasHDF5 = False
try:
    import h5py
    hasHDF5 = True
except:
    exu.Print('NGsolveCMStest: import of h5py failed; to test, requires to pip install h5py')
    
try:
    fem.SaveToFile(fileName2, mode='HDF5')
    fem.LoadFromFile(fileName2+'.hdf5')
except:
    if hasHDF5:
        raise ValueError('NGsolveCMStest: load/save with HDF5 mode failed even though h5py is installed!')
        #pass
    
#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
cms = ObjectFFRFreducedOrderInterface(fem)

objFFRF = cms.AddObjectFFRFreducedOrder(mbs, positionRef=[0,0,0], 
                                              initialVelocity=[0,0,0], 
                                              initialAngularVelocity=[0,0,0],
                                              gravity=[0,-9.81,0],
                                              color=[0.1,0.9,0.1,1.],
                                              )


#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
#add markers and joints
nodeDrawSize = 0.0025 #for joint drawing


mRB = mbs.AddMarker(MarkerNodeRigid(nodeNumber=objFFRF['nRigidBody']))
oGround = mbs.AddObject(ObjectGround(referencePosition= [0,0,0]))

if True:
    leftMidPoint = [0,0,0]
    
    mGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=leftMidPoint))

    mLeft = mbs.AddMarker(MarkerSuperElementRigid(bodyNumber=objFFRF['oFFRFreducedOrder'], 
                                                  meshNodeNumbers=np.array(nodesLeftPlane), #these are the meshNodeNumbers
                                                  weightingFactors=weightsLeftPlane))
    mbs.AddObject(GenericJoint(markerNumbers=[mGround, mLeft], 
                               constrainedAxes = [1,1,1,1,1,1*0],
                               visualization=VGenericJoint(axesRadius=0.1*a, axesLength=0.1*a)))

#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++
sensTipDispl = mbs.AddSensor(SensorSuperElement(bodyNumber=objFFRF['oFFRFreducedOrder'], 
                                                meshNodeNumber=nTip, #meshnode number!
                                                storeInternal=True,
                                                outputVariableType = exu.OutputVariableType.Displacement))
    
mbs.Assemble()

simulationSettings = exu.SimulationSettings()

SC.visualizationSettings.nodes.defaultSize = nodeDrawSize
SC.visualizationSettings.nodes.drawNodesAsPoint = False
SC.visualizationSettings.connectors.defaultSize = 2*nodeDrawSize

SC.visualizationSettings.nodes.show = False
SC.visualizationSettings.nodes.showBasis = True #of rigid body node of reference frame
SC.visualizationSettings.nodes.basisSize = 0.12
   
SC.visualizationSettings.sensors.show = True
SC.visualizationSettings.sensors.drawSimplified = False
SC.visualizationSettings.sensors.defaultSize = 0.01
SC.visualizationSettings.markers.drawSimplified = False
SC.visualizationSettings.markers.show = False
SC.visualizationSettings.markers.defaultSize = 0.01

SC.visualizationSettings.loads.drawSimplified = False


h=1e-3
tEnd = 0.1
if not testIsActive:
    tEnd = 4

simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
simulationSettings.timeIntegration.endTime = tEnd
simulationSettings.solution.file.write = False
simulationSettings.timeIntegration.verboseMode = 1
#simulationSettings.timeIntegration.verboseModeFile = 3
simulationSettings.timeIntegration.newton.useModifiedNewton = True

simulationSettings.solution.sensors.writePeriod = h

simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.8
simulationSettings.show.computationTime = True

#create animation:
SC.visualizationSettings.view0.window.renderWindowSize=[1920,1080]
SC.visualizationSettings.openGL.multiSampling = 4

if not testIsActive:
    SC.visualizationSettings.general.autoFitScene=False

    SC.renderer.Start()
    SC.renderer.RestoreSavedState() #load last model view

    SC.renderer.DoIdleTasks() #press space to continue

mbs.SolveDynamic(simulationSettings=simulationSettings)
    
uTip = mbs.GetSensorValues(sensTipDispl)
exu.Print("nModes=", nModes, ", tip displacement=", uTip)

    
if not testIsActive:
    SC.renderer.DoIdleTasks()
    SC.renderer.Stop() #safely close rendering window!

result = np.linalg.norm(uTip)
exu.Print('solution of NGsolveCMStest=',result)

#%%++++++++++++++++++++++++++++++++++++++++++++++++++++

exu.sys['testResult'] = result


#mbs.SolutionViewer()
```