NGsolvePistonEngine.py

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This is an EXUDYN example
#
# Details:  generate a piston engine with finite element mesh 
#           created with NGsolve and with variable number of pistons
#
# Author:   Johannes Gerstmayr
# Date:     2020-06-12
#
# 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 sys
import exudyn as exu

from exudyn.itemInterface import *
from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
import exudyn.graphics as graphics #only import if it does not conflict
from exudyn.rigidBodyUtilities import *
from exudyn.FEM import *

import time

import mkl
mkl.set_num_threads(20)

from ngsolve import *
from netgen.geom2d import unit_square

import netgen.libngpy as libng

netgenDrawing = False #set true, to show geometry and mesh in NETGEN
#if netgenDrawing, uncomment the following line and execute in external terminal, not in spyder (see preferences "Run"):
#import netgen.gui

from netgen.csg import *

import numpy as np
import timeit

verbose = True
meshSize = 0.005*2*2 #fast: 0.005*2; standard:0.005; fine: 0.0011: memory limit (96GB) for NGsolve; < 0.0015 makes problems with scipy eigensolver
meshOrder = 1 #2 for stresses!
showStresses = True #may take very long for large number of modes/nodes

#++++++++++++++++++++++++++++++++++++
#helper functions (copied from EXUDYN):
def RotationMatrixZ(angleRad):
    return np.array([ [np.cos(angleRad),-np.sin(angleRad), 0],
                      [np.sin(angleRad), np.cos(angleRad), 0],
                      [0,        0,        1] ]);
    
def VAdd(v0, v1):
    if len(v0) != len(v1): print("ERROR in VAdd: incompatible vectors!")
    n = len(v0)
    v = [0]*n
    for i in range(n):
        v[i] = v0[i]+v1[i]
    return v

def VSub(v0, v1):
    if len(v0) != len(v1): print("ERROR in VSub: incompatible vectors!")
    n = len(v0)
    v = [0]*n
    for i in range(n):
        v[i] = v0[i]-v1[i]
    return v

def NormL2(vector):
    value = 0
    for x in vector:
        value += x**2
    return value**0.5

def Normalize(v):
    v2=[0]*len(v)

    fact = NormL2(v)
    fact = 1./fact
    for i in range(len(v2)): 
        v2[i]=fact*v[i]
    return v2
#++++++++++++++++++++++++++++++++++++
startTotal = timeit.default_timer()
#parameters

#crank:
b1 = 0.012 #width of journal bearing
r1 = 0.012 #radius of journal bearing
dk = 0.015 #crank arm width (z)
bk = 0.032 #crank arm size (y)

l3 = 0.030
l4 = 0.040
#l4x= 0.005 #offset of counterweight
lk = 0.030 #l4*0.5+l3 #crank arm length (x)
bm = 0.065
dBevel = dk*0.5
#shaft:
r0 = 0.012 #0.012
d0 = 0.020 #shaft length at left/right support
d1 = 0.012 #shaft length at intermediate support

#distance rings:
db = 0.002          #width of distance ring
rdb0 = r0+db        #total radius of distance ring, shaft
rdb1 = r1+db        #total radius of distance ring, crank

#conrod:
bc = 0.024      #height of conrod
dc = 0.012      #width of conrod
lc = 0.080      #length of conrod (axis-axis)
r1o= r1+0.006   #outer radius of conrod at crank joint
r2 = 0.008      #radius of piston journal bearing
r2o= r2+0.006   #outer radius of conrod at piston joint

cylOffZ=0.010  #z-offset of cylinder cut out of conrod
cylR = 0.008    #radius of cylinder cut out of conrod

angC = 4*np.pi/180

#piston:
dpb = r2o-0.000   #axis inside piston
r2p = r2o+0.004   #0.018
lp = 0.034
bp = 0.050
lpAxis = dc+2*db
lOffCut = 0.011 #offset for cutout of big cylinder

#total length of one segment:
lTotal = db+dk+db+b1+db+dk+db+d1

#eps
eps = 5e-4 #added to faces, to avoid CSG-problems

#++++++++++++++++++++++++++++++++++++
#points
pLB = [0 ,0,-d0]
p0B = [0 ,0,0]
p1B = [0 ,0,db]
#p2B = [0, 0,db+dk]
p21B =[lk,0,db+dk]
p31B = [lk,0,db+dk+db]
p41B = [lk,0,db+dk+db+b1]
p51B =[lk,0,db+dk+db+b1+db]
p6B = [0 ,0,db+dk+db+b1+db+dk]
p7B = [0 ,0,db+dk+db+b1+db+dk+db]
p8B = [0 ,0,lTotal]

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

def CSGcube(pCenter,size):
    s2 = [0.5*size[0],0.5*size[1],0.5*size[2]]
    p0 = VSub(pCenter,s2)
    p1 = VAdd(pCenter,s2)
    brick = OrthoBrick(Pnt(p0[0],p0[1],p0[2]),Pnt(p1[0],p1[1],p1[2]))
    return brick


#transform points
def TransformCrank(p, zOff, zRot):
    p2 = RotationMatrixZ(zRot) @ p
    pOff=[0,0,zOff]
    return VAdd(p2,pOff)

#cube only in XY-plane, z infinite
def CSGcubeXY(pCenter,sizeX,sizeY,ex,ey):
    #print("pCenter=",pCenter)
    pl1 = Plane(Pnt(pCenter[0]-0.5*sizeX*ex[0],pCenter[1]-0.5*sizeX*ex[1],0),Vec(-ex[0],-ex[1],-ex[2]))
    pl2 = Plane(Pnt(pCenter[0]+0.5*sizeX*ex[0],pCenter[1]+0.5*sizeX*ex[1],0),Vec( ex[0], ex[1], ex[2]))

    pl3 = Plane(Pnt(pCenter[0]-0.5*sizeY*ey[0],pCenter[1]-0.5*sizeY*ey[1],0),Vec(-ey[0],-ey[1],-ey[2]))
    pl4 = Plane(Pnt(pCenter[0]+0.5*sizeY*ey[0],pCenter[1]+0.5*sizeY*ey[1],0),Vec( ey[0], ey[1], ey[2]))

    return pl1*pl2*pl3*pl4
    

#create one crank face at certain z-offset and rotation; side=1: left, side=-1: right
def GetCrankFace(zOff, zRot, side=1):
    ex = RotationMatrixZ(zRot) @ [1,0,0]
    ey = RotationMatrixZ(zRot) @ [0,1,0]
    #print("zOff=",zOff, "zRot=", zRot, "side=", side,"ex=", ex)
    pLeft = [0,0,zOff]
    pRight = [0,0,zOff+dk]
    pMid = [0,0,zOff+0.5*dk]

    pcLeft=VAdd(pLeft,lk*ex)
    pcRight=VAdd(pRight,lk*ex)
    f=0.5**0.5
    cyl1pl = Plane(Pnt(pcLeft[0],pcLeft[1],pcLeft[2]+0.5*dk-side*dk),Vec(f*ex[0],f*ex[1],f*ex[2]-side*f))        
    cyl1 = Cylinder(Pnt(pcLeft[0],pcLeft[1],pcLeft[2]-1), Pnt(pcRight[0],pcRight[1],pcRight[2]+1), 0.5*bk)*cyl1pl

    #cone2 = Cylinder(Pnt(pcLeft[0],pcLeft[1],pcLeft[2]-1), Pnt(pcRight[0],pcRight[1],pcRight[2]+1), lk+l4)
    cone2 = Cone(Pnt(pcLeft[0],pcLeft[1],pcLeft[2]-side*dBevel+0.5*dk), Pnt(pcLeft[0],pcLeft[1],pcLeft[2]+side*dBevel+0.5*dk), lk+l4-1.5*dBevel, lk+l4-0.5*dBevel)
    cube1 = CSGcubeXY(VAdd(pMid,0.49*l3*ex),1.02*l3,bk,ex,ey) #make l3 a little longer, to avoid bad edges
    cube2 = CSGcubeXY(VAdd(pMid,-0.5*l4*ex),1.0*l4,bm,ex,ey)*cone2

    pc3a = VAdd(pLeft,0.*l3*ex+(0.5*bk+0.4*l3)*ey)
    cyl3a = Cylinder(Pnt(pc3a[0],pc3a[1],pc3a[2]-1), Pnt(pc3a[0],pc3a[1],pc3a[2]+1), 0.42*l3)
    pc3b = VAdd(pLeft,0.*l3*ex+(-0.5*bk-0.4*l3)*ey)
    cyl3b = Cylinder(Pnt(pc3b[0],pc3b[1],pc3b[2]-1), Pnt(pc3b[0],pc3b[1],pc3b[2]+1), 0.42*l3)
    #cube3a = (CSGcubeXY(VAdd(pMid,0.26*l3*ex+(0.5*bk+0.26*l3)*ey),0.5*l3,0.5*l3,ex,ey)-cyl3a)
    
    return ((cube1+cube2+cyl1)-(cyl3a+cyl3b))*Plane(Pnt(0,0,pLeft[2]),Vec(0,0,-1))*Plane(Pnt(0,0,pRight[2]),Vec(0,0,1))
    #return (cube1+cube2+cyl1)*Plane(Pnt(0,0,pLeft[2]),Vec(0,0,-1))*Plane(Pnt(0,0,pRight[2]),Vec(0,0,1))

#generate one crank, rotated around z-axis in radiant
def GenerateCrank(zOff, zRot):
    pL = TransformCrank(pLB,zOff, zRot)
    p0 = TransformCrank(p0B,zOff, zRot)
    p1 = TransformCrank(p1B,zOff, zRot)

    p21 = TransformCrank(p21B,zOff, zRot)
    p31 = TransformCrank(p31B,zOff, zRot)
    p41 = TransformCrank(p41B,zOff, zRot)
    p51 = TransformCrank(p51B,zOff, zRot)

    p6 = TransformCrank(p6B,zOff, zRot)
    p7 = TransformCrank(p7B,zOff, zRot)
    p8 = TransformCrank(p8B,zOff, zRot)
    
    crank0 = CSGcylinder(pL,[p0[0],p0[1],p0[2]+eps],r0)
    crank1 = CSGcylinder(p0,[p1[0],p1[1],p1[2]+eps],rdb0)

    #conrod bearing:
    crank3 = CSGcylinder([p21[0],p21[1],p21[2]-eps],p31,rdb1)
    crank7 = CSGcylinder(p31,p41,r1)
    crank8 = CSGcylinder(p41,[p51[0],p51[1],p51[2]+eps],rdb1)
    
    crank9 = CSGcylinder([p6[0],p6[1],p6[2]-eps],p7,rdb0)
    crank10 = CSGcylinder([p7[0],p7[1],p7[2]-eps],p8,r0)

    #return crank0+crank1+crank3+crank4+crank5+crank6+crank7+crank8+crank4b+crank5b+crank6b+crank9+crank10
    if zOff==0:#add first shaft
        crank1 = crank1+crank0
    return crank1+GetCrankFace(db+zOff,zRot,1)+crank3+crank7+crank8+GetCrankFace(db+2*db+dk+b1+zOff,zRot,-1)+crank10+crank9


geoCrank = CSGeometry()

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#choose configuration for crankshaft:
#crankConfig = [0] #1-piston
#crankConfig = [np.pi/2] #1-piston
#crankConfig = [0,np.pi] #2-piston
#crankConfig = [0,np.pi*2./3.,2.*np.pi*2./3.] #3-piston
#crankConfig = [0,np.pi,np.pi,0] #4-piston
crankConfig = [0,np.pi*2./3.,2.*np.pi*2./3.,2.*np.pi*2./3.,np.pi*2./3.,0] #6-piston
#crankConfig = crankConfig*2 #12-piston

nPistons = len(crankConfig)

crank = GenerateCrank(0, crankConfig[0])
zPos = lTotal
for i in range(len(crankConfig)-1):
    angle = crankConfig[i+1]
    crank += GenerateCrank(zPos, angle)
    zPos += lTotal

# crank = (GenerateCrank(0, 0) + GenerateCrank(lTotal, np.pi*2./3.) + GenerateCrank(2*lTotal, np.pi*2.*2./3.)+
#           GenerateCrank(3*lTotal, np.pi*2.*2./3.) + GenerateCrank(4*lTotal, np.pi*2./3.))

geoCrank.Add(crank)

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#conrod model:
def GenerateConrod(zOff):
    ey0 = [0,1,0] #top/bottom face vector of conrod
    ey1 = [0,-1,0]

    ex0 = [1,0,0] #top/bottom face vector of conrod
    ex1 = [1,0,0]
    
    ey0 = RotationMatrixZ(-angC)@ey0
    ey1 = RotationMatrixZ(angC)@ey1
    ex0 = RotationMatrixZ(-angC)@ex0
    ex1 = RotationMatrixZ(angC)@ex1


    pl1 = Plane(Pnt(0, 0.5*bc,0),Vec(ey0[0],ey0[1],ey0[2]))
    pl2 = Plane(Pnt(0,-0.5*bc,0),Vec(ey1[0],ey1[1],ey1[2]))

    pl3 = Plane(Pnt(-0.5*lc,0,0),Vec(-1,0,0))
    pl4 = Plane(Pnt( 0.5*lc,0,0),Vec( 1,0,0))

    pl5 = Plane(Pnt( 0,0,-0.5*dc+zOff),Vec( 0,0,-1))
    pl6 = Plane(Pnt( 0,0, 0.5*dc+zOff),Vec( 0,0, 1))

    
    cylC1 = Cylinder(Pnt(-0.5*lc,0,-1), Pnt(-0.5*lc,0,1), r1)
    #cylC1o = Cylinder(Pnt(-0.5*lc,0,-1), Pnt(-0.5*lc,0,1), r1o)
    cylC1o = Sphere(Pnt(-0.5*lc,0,zOff), r1o) #in fact is a sphere

    cylC2 = Cylinder(Pnt( 0.5*lc,0,-1), Pnt( 0.5*lc,0,1), r2)
    #cylC2o = Cylinder(Pnt(0.5*lc,0,-1), Pnt( 0.5*lc,0,1), r2o)
    cylC2o = Sphere(Pnt(0.5*lc,0,zOff), r2o) #in fact is a sphere

    cylSideA = (Cylinder(Pnt(-0.5*lc+r1o,0,cylOffZ+zOff), Pnt(0.5*lc-r2o,0,cylOffZ+zOff), cylR)*
                Plane(Pnt(-0.5*lc+r1o-0.002,0,0),Vec(-1,0,0))*
                Plane(Pnt( 0.5*lc-r2o+0.002,0,0),Vec( 1,0,0)))

    cylSideB = (Cylinder(Pnt(-0.5*lc+r1o,0,-cylOffZ+zOff), Pnt(0.5*lc-r2o,0,-cylOffZ+zOff), cylR)*
                Plane(Pnt(-0.5*lc+r1o-0.002,0,0),Vec(-1,0,0))*
                Plane(Pnt( 0.5*lc-r2o+0.002,0,0),Vec( 1,0,0)))


    return ((pl1*pl2*pl3*pl4+cylC1o+cylC2o)-cylC1-cylC2)*pl5*pl6-cylSideA-cylSideB
    #return pl1*pl2*pl3*pl4*pl5*pl6

geoConrod = CSGeometry()
conrod = GenerateConrod(0)#db+dk+db+0.5*b1
geoConrod.Add(conrod)

# if netgenDrawing: 
#     Draw(geoCrank)

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#conrod model:
def GeneratePiston(zOff):
    p0 = [-dpb,0,zOff]
    p1 = [-dpb+lp,0,zOff]
    cylPo   = CSGcylinder(p0, p1, 0.5*bp) #piston outside
    cylPaxis= CSGcylinder([0,0,-0.5*lpAxis-eps+zOff],     [0,0, 0.5*lpAxis+eps+zOff], r2) #piston axis
    cylPaxis0= CSGcylinder([0,0,-0.5*lpAxis-eps+zOff],    [0,0,-0.5*lpAxis+db+zOff], r2+db) #piston axis
    cylPaxis1= CSGcylinder([0,0, 0.5*lpAxis-db+zOff], [0,0, 0.5*lpAxis+eps+zOff], r2+db) #piston axis
    cylPin  = CSGcylinder([0,0,-0.5*lpAxis+zOff], [0,0, 0.5*lpAxis+zOff], r2p) #piston inner cutout

    #box = CSGcube([0,0,zOff], [dpb+r2p,2*(r2p),lpAxis])
    box = CSGcube([-0.5*dpb,0,zOff], [dpb,2*(r2p)-0.002,lpAxis-0.000])

    cylCut  = CSGcylinder([-(l4+l3+lOffCut),0,-bp+zOff], [-(l4+l3+lOffCut),0, bp+zOff], l4+l3) #piston inner cutout

    return (cylPo-box-cylCut-cylPin)+cylPaxis+cylPaxis0+cylPaxis1

geoPiston = CSGeometry()
piston = GeneratePiston(0)#db+dk+db+0.5*b1
geoPiston.Add(piston)

if verbose: print("Generate meshes ...")
#do meshing, if geometry is successful
if True:
    meshCrank = Mesh( geoCrank.GenerateMesh(maxh=meshSize))
    meshCrank.Curve(1)
    if netgenDrawing: 
        Draw(meshCrank)
    
    if False:
        #save mesh to file:
        meshCrank.ngmesh.Export('testData/crankshaft.mesh','Neutral Format')

if True:
    meshConrod = Mesh( geoConrod.GenerateMesh(maxh=meshSize)) #in videos 0.003
    meshConrod.Curve(1)
    if netgenDrawing: 
        Draw(meshConrod)
    if False:
        meshConrod.ngmesh.Export('testData/conrod.mesh','Neutral Format')
    #+++++++++++++++++++++++++++++++++++++++++++++++++++++++

if True:
    meshPiston = Mesh( geoPiston.GenerateMesh(maxh=meshSize+0.001*0))
    meshPiston.Curve(1)
    if netgenDrawing: 
        Draw(meshPiston)
    if False:
        meshPiston.ngmesh.Export('testData/piston.mesh','Neutral Format')
    #+++++++++++++++++++++++++++++++++++++++++++++++++++++++

#here starts the EXUDYN part
if True:
    SC = exu.SystemContainer()
    mbs = SC.AddSystem()

    #crankshaft and piston mechanical parameters:
    density = 7850
    youngsModulus = 2.1e11 *1e-1
    poissonsRatio = 0.3
    fRotorStart = 20 #initial revolutions per second, only crankshaft

    totalFEcoordinates = 0 #accumulated FE-mesh coordinates
    #%%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    #import crankshaft mesh into EXUDYN FEMinterface
    femCrank = FEMinterface()
    eigenModesNGsolve=True
    nModes=8

    [bfM, bfK, fes] = femCrank.ImportMeshFromNGsolve(meshCrank, density, youngsModulus, poissonsRatio, 
                                                     verbose = True, meshOrder = meshOrder)
                          # computeEigenmodes=eigenModesNGsolve, excludeRigidBodyModes = 6,
                          # numberOfModes = nModes, maxEigensolveIterations=20)

    nModes = 20
    excludeRigidBodyModes = 6
    if verbose: print("number of coordinates crank =", femCrank.NumberOfCoordinates())
    if verbose: print("Compute eigenmodes crank ....")

    if not eigenModesNGsolve:
        startCrank = timeit.default_timer()
        femCrank.ComputeEigenmodes(nModes, excludeRigidBodyModes = excludeRigidBodyModes, useSparseSolver = True)
        stopCrank = timeit.default_timer()
        print("\ncrank eigen analysis time=", stopCrank-startCrank)
    else:
        start_time = time.time()
        femCrank.ComputeEigenmodesNGsolve(bfM, bfK, nModes=nModes, 
                                          excludeRigidBodyModes=excludeRigidBodyModes,  maxEigensolveIterations=20)
        print("NGsolve mode computation needed %.3f seconds" % (time.time() - start_time))
    
    totalFEcoordinates+=femCrank.NumberOfCoordinates()
    print("eigen freq. crank=", femCrank.GetEigenFrequenciesHz()[0:nModes])

    #+++++++++++++++++++++++++++++++++++++++++++++++++++++
    #compute stress modes:
    SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.Displacement
    mat = KirchhoffMaterial(youngsModulus, poissonsRatio, density)
    varType = exu.OutputVariableType.DisplacementLocal
    #varType = exu.OutputVariableType.StrainLocal
    if showStresses:
        print("ComputePostProcessingModes femCrank ... ")
        start_time = time.time()
        varType = exu.OutputVariableType.StressLocal
        femCrank.ComputePostProcessingModesNGsolve(fes, material=mat, 
                                       outputVariableType=varType)
        print("--- %s seconds ---" % (time.time() - start_time))
    
    SC.visualizationSettings.contour.outputVariable = varType
    
    #print("Create CMS object and matrices ....")
    cmsCrank = ObjectFFRFreducedOrderInterface(femCrank)
    
    #user functions should be defined outside of class:
    def UFmassFFRFreducedOrderCrank(mbs, t, itemIndex, qReduced, qReduced_t):
        return cmsCrank.UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    
    def UFforceFFRFreducedOrderCrank(mbs, t, itemIndex, qReduced, qReduced_t):
        return cmsCrank.UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    
    objFFRFcrank = cmsCrank.AddObjectFFRFreducedOrderWithUserFunctions(exu, mbs, 
                                                positionRef=[0,0,0], 
                                                eulerParametersRef=eulerParameters0, 
                                                initialVelocity=[0,0,0], initialAngularVelocity=[0,0,1*fRotorStart*2*pi],
                                                gravity = [0,-0*9.81,0],
                                                #UFforce=UFforceFFRFreducedOrderCrank, 
                                                #UFmassMatrix=UFmassFFRFreducedOrderCrank,
                                                color=[0.1,0.9,0.1,1.])
    mbs.SetObjectParameter(objFFRFcrank['oFFRFreducedOrder'],'VshowNodes',False)


    if False:#animate eigenmodes of crankshaft
        from exudyn.interactive import AnimateModes
        mbs.Assemble()

        SC.visualizationSettings.general.textSize = 16 #30 for cover figure
        SC.visualizationSettings.general.useGradientBackground = True
        SC.visualizationSettings.openGL.lineWidth = 2
        SC.visualizationSettings.openGL.showFaceEdges = True
        SC.visualizationSettings.openGL.showFaces = True
        SC.visualizationSettings.openGL.multiSampling = 4
        SC.visualizationSettings.nodes.show = False
        SC.visualizationSettings.window.renderWindowSize = [1600,1080]

        SC.visualizationSettings.contour.outputVariableComponent = 0

        SC.visualizationSettings.general.autoFitScene=False

        AnimateModes(SC, mbs, 1, period=0.2)
        exit()

    #%%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    #import conrod and piston mesh into EXUDYN FEMinterface and compute eigenmodes
    nModes = 8
    excludeRigidBodyModes = 6
    femConrod = FEMinterface()
    # femConrod.ImportMeshFromNGsolve(meshConrod, density, youngsModulus, poissonsRatio, verbose = False)
    [bfM, bfK, fes] = femConrod.ImportMeshFromNGsolve(meshConrod, density, youngsModulus, poissonsRatio, 
                                                      verbose = False, meshOrder = meshOrder)
                          # computeEigenmodes=eigenModesNGsolve, excludeRigidBodyModes = 6,
                          # numberOfModes = nModes, maxEigensolveIterations=20)
    if verbose: print("number of coordinates conrod =", femConrod.NumberOfCoordinates())
    if verbose: print("Compute eigenmodes conrod ....")

    if not eigenModesNGsolve:
        femConrod.ComputeEigenmodes(nModes, excludeRigidBodyModes = excludeRigidBodyModes, useSparseSolver = True)
    else:
        femConrod.ComputeEigenmodesNGsolve(bfM, bfK, nModes=nModes, excludeRigidBodyModes=excludeRigidBodyModes)

    totalFEcoordinates+=femConrod.NumberOfCoordinates()
    if verbose: print("eigen freq. conrod=", femConrod.GetEigenFrequenciesHz()[0:nModes])

    if showStresses:
        print("ComputePostProcessingModes femConrod ... ")
        start_time = time.time()
        femConrod.ComputePostProcessingModesNGsolve(fes, material=mat, 
                                       outputVariableType=varType)
        print("--- %s seconds ---" % (time.time() - start_time))

    #%%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    #import piston mesh into EXUDYN FEMinterface
    femPiston = FEMinterface()
    #femPiston.ImportMeshFromNGsolve(meshPiston, density, youngsModulus, poissonsRatio, verbose = False)
    [bfM, bfK, fes] = femPiston.ImportMeshFromNGsolve(meshPiston, density, youngsModulus, poissonsRatio, verbose = False, meshOrder = meshOrder)
    
    if verbose: print("number of coordinates piston =", femPiston.NumberOfCoordinates())
    if verbose: print("Compute eigenmodes piston ....")

    if not eigenModesNGsolve:
        femPiston.ComputeEigenmodes(nModes, excludeRigidBodyModes = excludeRigidBodyModes, useSparseSolver = True)
    else:
        femPiston.ComputeEigenmodesNGsolve(bfM, bfK, nModes=nModes, excludeRigidBodyModes=excludeRigidBodyModes)

    totalFEcoordinates+=femPiston.NumberOfCoordinates()
    if verbose: print("eigen freq. Piston=", femPiston.GetEigenFrequenciesHz()[0:nModes])

    if showStresses:
        print("ComputePostProcessingModes femPiston ... ")
        start_time = time.time()
        femPiston.ComputePostProcessingModesNGsolve(fes, material=mat, 
                                       outputVariableType=varType)
        print("--- %s seconds ---" % (time.time() - start_time))

    #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    #import multiple conrods and pistons

    #user functions should be defined outside of class:
    def UFmassFFRFreducedOrderConrod0(mbs, t, itemIndex, qReduced, qReduced_t):
        return cmsConrodList[0].UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFmassFFRFreducedOrderConrod1(mbs, t, itemIndex, qReduced, qReduced_t):
        return cmsConrodList[1].UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFmassFFRFreducedOrderConrod2(mbs, t, itemIndex, qReduced, qReduced_t):
        return cmsConrodList[2].UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFmassFFRFreducedOrderConrod3(mbs, t, itemIndex, qReduced, qReduced_t):
        return cmsConrodList[3].UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFmassFFRFreducedOrderConrod4(mbs, t, itemIndex, qReduced, qReduced_t):
        return cmsConrodList[4].UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFmassFFRFreducedOrderConrod5(mbs, t, itemIndex, qReduced, qReduced_t):
        return cmsConrodList[5].UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    
    def UFforceFFRFreducedOrderConrod0(mbs, t, itemIndex, qReduced, qReduced_t):
        return cmsConrodList[0].UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFforceFFRFreducedOrderConrod1(mbs, t, itemIndex, qReduced, qReduced_t):
        return cmsConrodList[1].UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFforceFFRFreducedOrderConrod2(mbs, t, itemIndex, qReduced, qReduced_t):
        return cmsConrodList[2].UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFforceFFRFreducedOrderConrod3(mbs, t, itemIndex, qReduced, qReduced_t):
        return cmsConrodList[3].UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFforceFFRFreducedOrderConrod4(mbs, t, itemIndex, qReduced, qReduced_t):
        return cmsConrodList[4].UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFforceFFRFreducedOrderConrod5(mbs, t, itemIndex, qReduced, qReduced_t):
        return cmsConrodList[5].UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)

    #user functions should be defined outside of class:
    def UFmassFFRFreducedOrderPiston0(mbs, t, itemIndex, qReduced, qReduced_t):
        return cmsPistonList[0].UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFmassFFRFreducedOrderPiston1(mbs, t, itemIndex, qReduced, qReduced_t):
        return cmsPistonList[1].UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFmassFFRFreducedOrderPiston2(mbs, t, itemIndex, qReduced, qReduced_t):
        return cmsPistonList[2].UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFmassFFRFreducedOrderPiston3(mbs, t, itemIndex, qReduced, qReduced_t):
        return cmsPistonList[3].UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFmassFFRFreducedOrderPiston4(mbs, t, itemIndex, qReduced, qReduced_t):
        return cmsPistonList[4].UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFmassFFRFreducedOrderPiston5(mbs, t, itemIndex, qReduced, qReduced_t):
        return cmsPistonList[5].UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    
    def UFforceFFRFreducedOrderPiston0(mbs, t, itemIndex, qReduced, qReduced_t):
        return cmsPistonList[0].UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFforceFFRFreducedOrderPiston1(mbs, t, itemIndex, qReduced, qReduced_t):
        return cmsPistonList[1].UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFforceFFRFreducedOrderPiston2(mbs, t, itemIndex, qReduced, qReduced_t):
        return cmsPistonList[2].UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFforceFFRFreducedOrderPiston3(mbs, t, itemIndex, qReduced, qReduced_t):
        return cmsPistonList[3].UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFforceFFRFreducedOrderPiston4(mbs, t, itemIndex, qReduced, qReduced_t):
        return cmsPistonList[4].UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFforceFFRFreducedOrderPiston5(mbs, t, itemIndex, qReduced, qReduced_t):
        return cmsPistonList[5].UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    
    #lists for multiple objects in conrods and pistons:
    UFmassFFRFreducedOrderConrodList=[UFmassFFRFreducedOrderConrod0,UFmassFFRFreducedOrderConrod1,
                                      UFmassFFRFreducedOrderConrod2,UFmassFFRFreducedOrderConrod3,
                                      UFmassFFRFreducedOrderConrod4,UFmassFFRFreducedOrderConrod5]
    UFforceFFRFreducedOrderConrodList=[UFforceFFRFreducedOrderConrod0,UFforceFFRFreducedOrderConrod1,
                                       UFforceFFRFreducedOrderConrod2,UFforceFFRFreducedOrderConrod3,
                                       UFforceFFRFreducedOrderConrod4,UFforceFFRFreducedOrderConrod5]
    objFFRFconrodList=[]
    cmsConrodList=[]
    UFmassFFRFreducedOrderPistonList=[UFmassFFRFreducedOrderPiston0,UFmassFFRFreducedOrderPiston1,
                                      UFmassFFRFreducedOrderPiston2,UFmassFFRFreducedOrderPiston3,
                                      UFmassFFRFreducedOrderPiston4,UFmassFFRFreducedOrderPiston5]
    UFforceFFRFreducedOrderPistonList=[UFforceFFRFreducedOrderPiston0,UFforceFFRFreducedOrderPiston1,
                                       UFforceFFRFreducedOrderPiston2,UFforceFFRFreducedOrderPiston3,
                                       UFforceFFRFreducedOrderPiston4,UFforceFFRFreducedOrderPiston5]
    objFFRFpistonList=[]
    cmsPistonList=[]
    pkList = []
    pcList = []
    ppList = []
    zOffsetList = []
    for iCrank in range(len(crankConfig)):
        zOffset = db+dk+db + lTotal*iCrank #left end of conrod, for multiple conrods in a loop
        zOffsetList.append(zOffset)
        #compute crank (pK), conrod (pC) and piston position (pP) for any crank angle:
        phi = crankConfig[iCrank]
        pK = np.array([lk*np.cos(phi),lk*np.sin(phi),0])
        alpha=np.arcsin(pK[1]/lc)
        pC = pK + np.array([0.5*lc*np.cos(alpha),-0.5*lc*np.sin(alpha),0])
        pP = pK + np.array([lc*np.cos(alpha),-lc*np.sin(alpha),0])
        pkList.append(pK)
        pcList.append(pC)
        ppList.append(pP)
        #print("pK=",pK)
        #print("pC=",pC)
        #print("pP=",pP)
        
        eulerParametersInit = RotationMatrix2EulerParameters(RotationMatrixZ(-alpha))
        #pRef = [lk+0.5*lc,0,zOffset+0.5*b1] #0-degree
        pRef = pC + [0,0,zOffset+0.5*b1]
                    
        #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        #import conrod CMS
        cmsConrod = ObjectFFRFreducedOrderInterface(femConrod)
        cmsConrodList.append(cmsConrod)
        objFFRFconrod = cmsConrod.AddObjectFFRFreducedOrderWithUserFunctions(exu, mbs, 
                                                    positionRef=pRef, 
                                                    eulerParametersRef=eulerParametersInit, 
                                                    initialVelocity=[0,0,0], 
                                                    initialAngularVelocity=[0,0,0*fRotorStart*2*pi],
                                                    gravity = [0,-0*9.81,0],
                                                    #UFforce=UFforceFFRFreducedOrderConrodList[iCrank], 
                                                    #UFmassMatrix=UFmassFFRFreducedOrderConrodList[iCrank],
                                                    color=[0.1,0.9,0.1,1.])
        mbs.SetObjectParameter(objFFRFconrod['oFFRFreducedOrder'],'VshowNodes',False)
        objFFRFconrodList.append(objFFRFconrod)
    
        #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        #import piston CMS
        cmsPiston = ObjectFFRFreducedOrderInterface(femPiston)
        cmsPistonList.append(cmsPiston)
 
        objFFRFpiston = cmsPiston.AddObjectFFRFreducedOrderWithUserFunctions(exu, mbs, 
                                                    positionRef=pP+[0,0,zOffset+0.5*b1], 
                                                    eulerParametersRef=eulerParameters0, 
                                                    initialVelocity=[0,0,0], initialAngularVelocity=[0,0,0*fRotorStart*2*pi],
                                                    gravity = [0,-0*9.81,0],
                                                    #UFforce=UFforceFFRFreducedOrderPistonList[iCrank], 
                                                    #UFmassMatrix=UFmassFFRFreducedOrderPistonList[iCrank],
                                                    color=[0.1,0.9,0.1,1.])
        mbs.SetObjectParameter(objFFRFpiston['oFFRFreducedOrder'],'VshowNodes',False)
        objFFRFpistonList.append(objFFRFpiston)

    #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    if True: #connect bodies:
        k = 1e6         #joint stiffness
        d = k*0.002     #joint damping
        nMarkerPerPiston = 10    #number of markers per crank/conrod/piston part

        genMarkerPos = [[0,0,-d0],[0,0,lTotal*nPistons]]
        genMarkerR   = [r0,r0]
        genMarkerFEM = [femCrank,femCrank]
        genMarkerObject = [objFFRFcrank,objFFRFcrank]

        for iCrank in range(len(crankConfig)):
            genMarkerPos += [pkList[iCrank]+[0,0,zOffsetList[iCrank]],pkList[iCrank]+[0,0,zOffsetList[iCrank]+b1],
                            [-0.5*lc,0,-0.5*dc],[-0.5*lc,0, 0.5*dc],[0.5*lc,0,-0.5*dc],[0.5*lc,0, 0.5*dc],
                            [0,0,-0.5*dc],[0,0,0.5*dc], [-dpb,0,0],[lp-dpb,0,0]]
            genMarkerR   += [r1,r1,
                            r1,r1,r2,r2,
                            r2,r2,0.5*bp,0.5*bp]
            genMarkerFEM += [femCrank,femCrank,
                            femConrod,femConrod,femConrod,femConrod,
                            femPiston,femPiston,femPiston,femPiston]
            genMarkerObject += [objFFRFcrank,objFFRFcrank,
                               objFFRFconrodList[iCrank],objFFRFconrodList[iCrank],objFFRFconrodList[iCrank],objFFRFconrodList[iCrank],
                               objFFRFpistonList[iCrank],objFFRFpistonList[iCrank],objFFRFpistonList[iCrank],objFFRFpistonList[iCrank]]

        markerList = []
        #generate markers for joints:
        for i in range(len(genMarkerPos)):
            p = genMarkerPos[i]
            nodeList=[]
            if p[2] != 0:
                nodeList= genMarkerFEM[i].GetNodesOnCircle(p, [0,0,1], genMarkerR[i])
            else:
                nodeList= genMarkerFEM[i].GetNodesOnCircle(p, [1,0,0], genMarkerR[i])
            #print("nodeList"+str(i)+":", nodeList)
            lenNodeList = len(nodeList)
            weights = np.array((1./lenNodeList)*np.ones(lenNodeList))
        
            markerList += [mbs.AddMarker(MarkerSuperElementPosition(bodyNumber=genMarkerObject[i]['oFFRFreducedOrder'], 
                                                            meshNodeNumbers=np.array(nodeList), #these are the meshNodeNumbers
                                                            weightingFactors=weights))]

        oGround = mbs.AddObject(ObjectGround(referencePosition= [0,0,0]))
        
        mGroundPosLeft = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition=genMarkerPos[0]))
        mGroundPosRight = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition=genMarkerPos[1]))
        

        #joints for crankshaft/ground
        oSJleft = mbs.AddObject(CartesianSpringDamper(markerNumbers=[mGroundPosLeft, markerList[0]],
                                            stiffness=[k,k,k], damping=[d,d,d]))
        oSJright = mbs.AddObject(CartesianSpringDamper(markerNumbers=[mGroundPosRight, markerList[1]],
                                            stiffness=[k,k,k], damping=[d,d,d]))

        for iCrank in range(len(crankConfig)):
            mOff = nMarkerPerPiston*iCrank
            #joints for crankshaft/conrod:
            oJointCCleft = mbs.AddObject(CartesianSpringDamper(markerNumbers=[markerList[mOff+2], markerList[mOff+4]],
                                                stiffness=[k,k,k], damping=[d,d,d]))
            oJointCCright= mbs.AddObject(CartesianSpringDamper(markerNumbers=[markerList[mOff+3], markerList[mOff+5]],
                                                stiffness=[k,k,k], damping=[d,d,d]))
    
            #joints for conrod/piston:
            oJointCPleft = mbs.AddObject(CartesianSpringDamper(markerNumbers=[markerList[mOff+6], markerList[mOff+8]],
                                                stiffness=[k,k,k], damping=[d,d,d]))
            oJointCPright= mbs.AddObject(CartesianSpringDamper(markerNumbers=[markerList[mOff+7], markerList[mOff+9]],
                                                stiffness=[k,k,k], damping=[d,d,d]))
    
            mGroundPosPiston = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, 
                                                                localPosition=[ppList[iCrank][0],0,zOffsetList[iCrank]+0.5*b1]))
            oJointPGleft = mbs.AddObject(CartesianSpringDamper(markerNumbers=[mGroundPosPiston, markerList[mOff+10]],
                                                stiffness=[0,k,k], damping=[0,d,d]))
            oJointPGright = mbs.AddObject(CartesianSpringDamper(markerNumbers=[mGroundPosPiston, markerList[mOff+11]],
                                                stiffness=[0,k,k], damping=[0,d,d]))


    stopTotal = timeit.default_timer()
    print("\ntotal elapsed time=", stopTotal-startTotal)
    mbs.Assemble()

    #now simulate model in exudyn:
    #%%+++++++++++++++++++++
    if True:
        print("totalFEcoordinates=",totalFEcoordinates)
        
        simulationSettings = exu.SimulationSettings()
        
        nodeDrawSize = 0.0005
        SC.visualizationSettings.general.textSize = 14 #30 for cover figure
        SC.visualizationSettings.general.useGradientBackground = True
        SC.visualizationSettings.openGL.lineWidth = 2

        SC.visualizationSettings.nodes.defaultSize = nodeDrawSize
        SC.visualizationSettings.nodes.drawNodesAsPoint = False
        SC.visualizationSettings.connectors.defaultSize = 2*nodeDrawSize
        SC.visualizationSettings.connectors.show = False
        
        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.bodies.deformationScaleFactor = 1 #use this factor to scale the deformation of modes
        
        SC.visualizationSettings.openGL.showFaceEdges = True
        SC.visualizationSettings.openGL.showFaces = True
        SC.visualizationSettings.openGL.multiSampling = 4
        
        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
        
        #SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.Displacement
        SC.visualizationSettings.contour.outputVariableComponent = -1
        SC.visualizationSettings.contour.reduceRange = True
        #SC.visualizationSettings.contour.automaticRange = False
        #SC.visualizationSettings.contour.maxValue = 3e7
        # SC.visualizationSettings.contour.minValue = -0.0003
        # SC.visualizationSettings.contour.maxValue =  0.0003
        
        simulationSettings.solutionSettings.solutionInformation = "NGsolve/NETGEN engine test"
        
        h=0.05e-3
        tEnd = 2
        
        simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
        simulationSettings.timeIntegration.endTime = tEnd
        simulationSettings.solutionSettings.solutionWritePeriod = h*10 #writing already costs much time
        simulationSettings.timeIntegration.verboseMode = 1
        #simulationSettings.timeIntegration.verboseModeFile = 3
        simulationSettings.timeIntegration.newton.useModifiedNewton = True
        
        simulationSettings.solutionSettings.sensorsWritePeriod = h
        #simulationSettings.solutionSettings.coordinatesSolutionFileName = "solution/coordinatesSolution.txt"
        simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse #faster, because system size already quite large
        
        simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5 #SHOULD work with 0.9 as well
        simulationSettings.displayStatistics = True
        #simulationSettings.displayComputationTime = True
        SC.visualizationSettings.general.autoFitScene = False #for reloading of renderState to work
        
        #create animation:
        if False:
            simulationSettings.solutionSettings.recordImagesInterval = 0.001
            SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
            SC.visualizationSettings.window.renderWindowSize=[1920,1080]

        exu.StartRenderer()
        if 'renderState' in exu.sys: SC.SetRenderState(exu.sys['renderState']) #load last model view

        mbs.WaitForUserToContinue() #press space to continue
        
        simulate = True #set false to show last stored solution
        if simulate:
            mbs.SolveDynamic(simulationSettings)
        else:
            SC.visualizationSettings.general.autoFitScene = False
            sol = LoadSolutionFile('coordinatesSolution.txt')
            if False: #directly show animation
                AnimateSolution(mbs, solution=sol, rowIncrement = 1, timeout=0.01, 
                                createImages = False, runLoop = True)
            else: #interact with animation
                
                mbs.SolutionViewer(sol, rowIncrement=1, timeout=0.02)


        if False: #draw with matplotlib, export as pdf
            SC.visualizationSettings.exportImages.saveImageFormat = "TXT"
            SC.visualizationSettings.exportImages.saveImageAsTextTriangles=True
            SC.RedrawAndSaveImage() #uses default filename
            
            from exudyn.plot import LoadImage, PlotImage

            # plot 2D
            # data = LoadImage('images/frame00000.txt', trianglesAsLines=True)
            # PlotImage(data, HT=HomogeneousTransformation(RotationMatrixZ(0.5*pi)@RotationMatrixX(0.5*pi), [0,0,0]), 
            #           lineWidths=0.5, lineStyles='-', title='', closeAll=True, plot3D=False,
            #           fileName='images/test.pdf')
            
            data = LoadImage('images/frame00000.txt', trianglesAsLines=False)
            PlotImage(data, HT=HomogeneousTransformation(2.5*RotationMatrixZ(0.5*pi)@RotationMatrixY(-0.5*pi), [0,1,0.25]), 
                      lineWidths=0.5, lineStyles='-', triangleEdgeColors='black', triangleEdgeWidths=0.25, title='', closeAll=True, plot3D=True,
                      fileName='images/test3D.pdf')
                    
        SC.WaitForRenderEngineStopFlag()
        exu.StopRenderer() #safely close rendering window!
        lastRenderState = SC.GetRenderState() #store model view for next simulation