#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This is an EXUDYN example
#
# Details:  generate a piston engine with variable number of pistons
#
# Author:   <NAME>
# 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
sys.path.append('C:/ProgramData/ngsolve/lib/site-packages')

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:
import netgen.gui

from netgen.csg import *

import numpy as np
import timeit

verbose = True

#++++++++++++++++++++++++++++++++++++
#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
meshSize = 0.0025 #standard: 0.005; 0.0015 does not always work

#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

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)
    #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)
    meshConrod.ngmesh.Export('testData/conrod.mesh','Neutral Format')
    #+++++++++++++++++++++++++++++++++++++++++++++++++++++++

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

#here starts the EXUDYN part
if True:
    import sys
    #C:\DATA\cpp\EXUDYN_git\main\bin\WorkingRelease64P37
    #sys.path.append('C:/DATA/cpp/EXUDYN_git/main/bin/WorkingRelease64P37') #for exudyn, itemInterface and from exudyn.utilities import *
    import exudyn as exu
    
    from exudyn.itemInterface import *
    from exudyn.utilities import *
    from exudyn.rigidBodyUtilities import *
    from exudyn.FEM import *
    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()
    femCrank.ImportMeshFromNGsolve(meshCrank, density, youngsModulus, poissonsRatio, verbose = True)

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

    startCrank = timeit.default_timer()
    femCrank.ComputeEigenmodes(nModes, excludeRigidBodyModes = excludeRigidBodyModes, useSparseSolver = True)
    totalFEcoordinates+=femCrank.NumberOfCoordinates()
    stopCrank = timeit.default_timer()
    print("\ncrank eigen analysis time=", stopCrank-startCrank)

    print("eigen freq. crank=", femCrank.GetEigenFrequenciesHz()[0:nModes])
    
    #print("Create CMS object and matrices ....")
    cmsCrank = ObjectFFRFreducedOrderInterface(femCrank)
    
    #user functions should be defined outside of class:
    def UFmassFFRFreducedOrderCrank(t, qReduced, qReduced_t):
        return cmsCrank.UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    
    def UFforceFFRFreducedOrderCrank(t, 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)


    #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    #import conrod and piston mesh into EXUDYN FEMinterface and compute eigenmodes
    femConrod = FEMinterface()
    femConrod.ImportMeshFromNGsolve(meshConrod, density, youngsModulus, poissonsRatio, verbose = False)
    if verbose: print("number of coordinates conrod =", femConrod.NumberOfCoordinates())
    if verbose: print("Compute eigenmodes conrod ....")

    nModes = 8
    excludeRigidBodyModes = 6
    #print("Compute eigenmodes Scipy ....")

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

    #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    #import piston mesh into EXUDYN FEMinterface
    femPiston = FEMinterface()
    femPiston.ImportMeshFromNGsolve(meshPiston, density, youngsModulus, poissonsRatio, verbose = False)

    nModes = 8
    excludeRigidBodyModes = 6

    if verbose: print("number of coordinates piston =", femPiston.NumberOfCoordinates())
    if verbose: print("Compute eigenmodes piston ....")

    femPiston.ComputeEigenmodes(nModes, excludeRigidBodyModes = excludeRigidBodyModes, useSparseSolver = True)
    totalFEcoordinates+=femPiston.NumberOfCoordinates()
    if verbose: print("eigen freq. Piston=", femPiston.GetEigenFrequenciesHz()[0:nModes])
    #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    #import multiple conrods and pistons

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

    #user functions should be defined outside of class:
    def UFmassFFRFreducedOrderPiston0(t, qReduced, qReduced_t):
        return cmsPistonList[0].UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFmassFFRFreducedOrderPiston1(t, qReduced, qReduced_t):
        return cmsPistonList[1].UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFmassFFRFreducedOrderPiston2(t, qReduced, qReduced_t):
        return cmsPistonList[2].UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFmassFFRFreducedOrderPiston3(t, qReduced, qReduced_t):
        return cmsPistonList[3].UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFmassFFRFreducedOrderPiston4(t, qReduced, qReduced_t):
        return cmsPistonList[4].UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFmassFFRFreducedOrderPiston5(t, qReduced, qReduced_t):
        return cmsPistonList[5].UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    
    def UFforceFFRFreducedOrderPiston0(t, qReduced, qReduced_t):
        return cmsPistonList[0].UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFforceFFRFreducedOrderPiston1(t, qReduced, qReduced_t):
        return cmsPistonList[1].UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFforceFFRFreducedOrderPiston2(t, qReduced, qReduced_t):
        return cmsPistonList[2].UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFforceFFRFreducedOrderPiston3(t, qReduced, qReduced_t):
        return cmsPistonList[3].UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFforceFFRFreducedOrderPiston4(t, qReduced, qReduced_t):
        return cmsPistonList[4].UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFforceFFRFreducedOrderPiston5(t, 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)
    #now simulate model in exudyn:
    if True:
        mbs.Assemble()
        print("totalFEcoordinates=",totalFEcoordinates)
        
        simulationSettings = exu.SimulationSettings()
        
        nodeDrawSize = 0.0005
        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 = 0
        SC.visualizationSettings.contour.reduceRange = False
        # SC.visualizationSettings.contour.minValue = -0.0003
        # SC.visualizationSettings.contour.maxValue =  0.0003
        
        simulationSettings.solutionSettings.solutionInformation = "NGsolve/NETGEN engine test"
        
        h=0.25e-3
        tEnd = 0.5
        
        simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
        simulationSettings.timeIntegration.endTime = tEnd
        simulationSettings.solutionSettings.solutionWritePeriod = h
        simulationSettings.timeIntegration.verboseMode = 1
        #simulationSettings.timeIntegration.verboseModeFile = 3
        simulationSettings.timeIntegration.newton.useModifiedNewton = True
        
        simulationSettings.solutionSettings.sensorsWritePeriod = h
        #simulationSettings.solutionSettings.coordinatesSolutionFileName = "solution/coordinatesSolution.txt"
        
        simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5 #SHOULD work with 0.9 as well
        simulationSettings.displayStatistics = True
        simulationSettings.displayComputationTime = True
        
        #create animation:
        #simulationSettings.solutionSettings.recordImagesInterval = 0.002
        #SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
        #SC.visualizationSettings.window.renderWindowSize=[1920,1080]
        
        exu.StartRenderer()
        lastRenderState = {'centerPoint': [-0.1251770257949829, -0.07375780493021011, 0.0],
                             'maxSceneSize': 0.10000000149011612,
                             'zoom': 0.15295113623142242,
                             'currentWindowSize': [1024, 768],
                             'modelRotation': np.array([[ 0.484963  ,  0.04052661,  0.87359524],
                                    [-0.3886355 ,  0.9048576 ,  0.17376813],
                                    [-0.783437  , -0.42378116,  0.4545725 ]])}
        SC.SetRenderState(lastRenderState) #load last model view
    
        mbs.WaitForUserToContinue() #press space to continue
        mbs.WaitForUserToContinue() #press space to continue
        
        SC.TimeIntegrationSolve(mbs, 'GeneralizedAlpha', simulationSettings)
        
        SC.WaitForRenderEngineStopFlag()
        exu.StopRenderer() #safely close rendering window!
        lastRenderState = SC.GetRenderState() #store model view for next simulation
    













#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This is an EXUDYN example
#
# Details:  generate a piston engine with variable number of pistons
#
# Author:   <NAME>
# 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
sys.path.append('C:/ProgramData/ngsolve/lib/site-packages')

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:
import netgen.gui

from netgen.csg import *

import numpy as np
import timeit

verbose = True

#++++++++++++++++++++++++++++++++++++
#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
meshSize = 0.0025 #standard: 0.005; 0.0015 does not always work

#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

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)
    #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)
    meshConrod.ngmesh.Export('testData/conrod.mesh','Neutral Format')
    #+++++++++++++++++++++++++++++++++++++++++++++++++++++++

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

#here starts the EXUDYN part
if True:
    import sys
    #C:\DATA\cpp\EXUDYN_git\main\bin\WorkingRelease64P37
    #sys.path.append('C:/DATA/cpp/EXUDYN_git/main/bin/WorkingRelease64P37') #for exudyn, itemInterface and from exudyn.utilities import *
    import exudyn as exu
    
    from exudyn.itemInterface import *
    from exudyn.utilities import *
    from exudyn.rigidBodyUtilities import *
    from exudyn.FEM import *
    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()
    femCrank.ImportMeshFromNGsolve(meshCrank, density, youngsModulus, poissonsRatio, verbose = True)

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

    startCrank = timeit.default_timer()
    femCrank.ComputeEigenmodes(nModes, excludeRigidBodyModes = excludeRigidBodyModes, useSparseSolver = True)
    totalFEcoordinates+=femCrank.NumberOfCoordinates()
    stopCrank = timeit.default_timer()
    print("\ncrank eigen analysis time=", stopCrank-startCrank)

    print("eigen freq. crank=", femCrank.GetEigenFrequenciesHz()[0:nModes])
    
    #print("Create CMS object and matrices ....")
    cmsCrank = ObjectFFRFreducedOrderInterface(femCrank)
    
    #user functions should be defined outside of class:
    def UFmassFFRFreducedOrderCrank(t, qReduced, qReduced_t):
        return cmsCrank.UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    
    def UFforceFFRFreducedOrderCrank(t, 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)


    #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    #import conrod and piston mesh into EXUDYN FEMinterface and compute eigenmodes
    femConrod = FEMinterface()
    femConrod.ImportMeshFromNGsolve(meshConrod, density, youngsModulus, poissonsRatio, verbose = False)
    if verbose: print("number of coordinates conrod =", femConrod.NumberOfCoordinates())
    if verbose: print("Compute eigenmodes conrod ....")

    nModes = 8
    excludeRigidBodyModes = 6
    #print("Compute eigenmodes Scipy ....")

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

    #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    #import piston mesh into EXUDYN FEMinterface
    femPiston = FEMinterface()
    femPiston.ImportMeshFromNGsolve(meshPiston, density, youngsModulus, poissonsRatio, verbose = False)

    nModes = 8
    excludeRigidBodyModes = 6

    if verbose: print("number of coordinates piston =", femPiston.NumberOfCoordinates())
    if verbose: print("Compute eigenmodes piston ....")

    femPiston.ComputeEigenmodes(nModes, excludeRigidBodyModes = excludeRigidBodyModes, useSparseSolver = True)
    totalFEcoordinates+=femPiston.NumberOfCoordinates()
    if verbose: print("eigen freq. Piston=", femPiston.GetEigenFrequenciesHz()[0:nModes])
    #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    #import multiple conrods and pistons

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

    #user functions should be defined outside of class:
    def UFmassFFRFreducedOrderPiston0(t, qReduced, qReduced_t):
        return cmsPistonList[0].UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFmassFFRFreducedOrderPiston1(t, qReduced, qReduced_t):
        return cmsPistonList[1].UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFmassFFRFreducedOrderPiston2(t, qReduced, qReduced_t):
        return cmsPistonList[2].UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFmassFFRFreducedOrderPiston3(t, qReduced, qReduced_t):
        return cmsPistonList[3].UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFmassFFRFreducedOrderPiston4(t, qReduced, qReduced_t):
        return cmsPistonList[4].UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFmassFFRFreducedOrderPiston5(t, qReduced, qReduced_t):
        return cmsPistonList[5].UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    
    def UFforceFFRFreducedOrderPiston0(t, qReduced, qReduced_t):
        return cmsPistonList[0].UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFforceFFRFreducedOrderPiston1(t, qReduced, qReduced_t):
        return cmsPistonList[1].UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFforceFFRFreducedOrderPiston2(t, qReduced, qReduced_t):
        return cmsPistonList[2].UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFforceFFRFreducedOrderPiston3(t, qReduced, qReduced_t):
        return cmsPistonList[3].UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFforceFFRFreducedOrderPiston4(t, qReduced, qReduced_t):
        return cmsPistonList[4].UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
    def UFforceFFRFreducedOrderPiston5(t, 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)
    #now simulate model in exudyn:
    if True:
        mbs.Assemble()
        print("totalFEcoordinates=",totalFEcoordinates)
        
        simulationSettings = exu.SimulationSettings()
        
        nodeDrawSize = 0.0005
        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 = 0
        SC.visualizationSettings.contour.reduceRange = False
        # SC.visualizationSettings.contour.minValue = -0.0003
        # SC.visualizationSettings.contour.maxValue =  0.0003
        
        simulationSettings.solutionSettings.solutionInformation = "NGsolve/NETGEN engine test"
        
        h=0.25e-3
        tEnd = 0.5
        
        simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
        simulationSettings.timeIntegration.endTime = tEnd
        simulationSettings.solutionSettings.solutionWritePeriod = h
        simulationSettings.timeIntegration.verboseMode = 1
        #simulationSettings.timeIntegration.verboseModeFile = 3
        simulationSettings.timeIntegration.newton.useModifiedNewton = True
        
        simulationSettings.solutionSettings.sensorsWritePeriod = h
        #simulationSettings.solutionSettings.coordinatesSolutionFileName = "solution/coordinatesSolution.txt"
        
        simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5 #SHOULD work with 0.9 as well
        simulationSettings.displayStatistics = True
        simulationSettings.displayComputationTime = True
        
        #create animation:
        #simulationSettings.solutionSettings.recordImagesInterval = 0.002
        #SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
        #SC.visualizationSettings.window.renderWindowSize=[1920,1080]
        
        exu.StartRenderer()
        lastRenderState = {'centerPoint': [-0.1251770257949829, -0.07375780493021011, 0.0],
                             'maxSceneSize': 0.10000000149011612,
                             'zoom': 0.15295113623142242,
                             'currentWindowSize': [1024, 768],
                             'modelRotation': np.array([[ 0.484963  ,  0.04052661,  0.87359524],
                                    [-0.3886355 ,  0.9048576 ,  0.17376813],
                                    [-0.783437  , -0.42378116,  0.4545725 ]])}
        SC.SetRenderState(lastRenderState) #load last model view
    
        mbs.WaitForUserToContinue() #press space to continue
        mbs.WaitForUserToContinue() #press space to continue
        
        SC.TimeIntegrationSolve(mbs, 'GeneralizedAlpha', simulationSettings)
        
        SC.WaitForRenderEngineStopFlag()
        exu.StopRenderer() #safely close rendering window!
        lastRenderState = SC.GetRenderState() #store model view for next simulation
    













#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # This is an EXUDYN example # # Details:  generate a piston engine with variable number of pistons # # Author:   <NAME> # 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. # #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #set true, to show geometry and mesh in NETGEN #if netgenDrawing: #++++++++++++++++++++++++++++++++++++ #helper functions (copied from EXUDYN): #++++++++++++++++++++++++++++++++++++ #parameters #standard: 0.005; 0.0015 does not always work #crank: #width of journal bearing #radius of journal bearing #crank arm width (z) #crank arm size (y) #l4x= 0.005 #offset of counterweight #l4*0.5+l3 #crank arm length (x) #shaft: #0.012 #shaft length at left/right support #shaft length at intermediate support #distance rings: #width of distance ring #total radius of distance ring, shaft #total radius of distance ring, crank #conrod: #height of conrod #width of conrod #length of conrod (axis-axis) #outer radius of conrod at crank joint #radius of piston journal bearing #outer radius of conrod at piston joint #z-offset of cylinder cut out of conrod #radius of cylinder cut out of conrod #piston: #axis inside piston #0.018 #offset for cutout of big cylinder #total length of one segment: #eps #added to faces, to avoid CSG-problems #++++++++++++++++++++++++++++++++++++ #points #p2B = [0, 0,db+dk] #transform points #cube only in XY-plane, z infinite #print("pCenter=",pCenter) #create one crank face at certain z-offset and rotation; side=1: left, side=-1: right #print("zOff=",zOff, "zRot=", zRot, "side=", side,"ex=", ex) #cone2 = Cylinder(Pnt(pcLeft[0],pcLeft[1],pcLeft[2]-1), Pnt(pcRight[0],pcRight[1],pcRight[2]+1), lk+l4) #make l3 a little longer, to avoid bad edges #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)*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 #conrod bearing: #return crank0+crank1+crank3+crank4+crank5+crank6+crank7+crank8+crank4b+crank5b+crank6b+crank9+crank10 #add first shaft #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #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 #6-piston # 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.)) #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #conrod model: #top/bottom face vector of conrod #top/bottom face vector of conrod #cylC1o = Cylinder(Pnt(-0.5*lc,0,-1), Pnt(-0.5*lc,0,1), r1o) #in fact is a sphere #cylC2o = Cylinder(Pnt(0.5*lc,0,-1), Pnt( 0.5*lc,0,1), r2o) #in fact is a sphere #return pl1*pl2*pl3*pl4*pl5*pl6 #db+dk+db+0.5*b1 # if netgenDrawing: #     Draw(geoCrank) #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #conrod model: #piston outside #piston axis #piston axis #piston axis #piston inner cutout #box = CSGcube([0,0,zOff], [dpb+r2p,2*(r2p),lpAxis]) #piston inner cutout #db+dk+db+0.5*b1 #do meshing, if geometry is successful #save mesh to file: #in videos 0.003 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++ #+++++++++++++++++++++++++++++++++++++++++++++++++++++++ #here starts the EXUDYN part #C:\DATA\cpp\EXUDYN_git\main\bin\WorkingRelease64P37 #sys.path.append('C:/DATA/cpp/EXUDYN_git/main/bin/WorkingRelease64P37') #for exudyn, itemInterface and from exudyn.utilities import * #crankshaft and piston mechanical parameters: #initial revolutions per second, only crankshaft #accumulated FE-mesh coordinates #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #import crankshaft mesh into EXUDYN FEMinterface #print("Create CMS object and matrices ....") #user functions should be defined outside of class: #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #import conrod and piston mesh into EXUDYN FEMinterface and compute eigenmodes #print("Compute eigenmodes Scipy ....") #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #import piston mesh into EXUDYN FEMinterface #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #import multiple conrods and pistons #user functions should be defined outside of class: #user functions should be defined outside of class: #lists for multiple objects in conrods and pistons: #left end of conrod, for multiple conrods in a loop #compute crank (pK), conrod (pC) and piston position (pP) for any crank angle: #print("pK=",pK) #print("pC=",pC) #print("pP=",pP) #pRef = [lk+0.5*lc,0,zOffset+0.5*b1] #0-degree #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #import conrod CMS #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #import piston CMS #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #connect bodies: #joint stiffness #joint damping #number of markers per crank/conrod/piston part #generate markers for joints: #print("nodeList"+str(i)+":", nodeList) #these are the meshNodeNumbers #joints for crankshaft/ground #joints for crankshaft/conrod: #joints for conrod/piston: #now simulate model in exudyn: #of rigid body node of reference frame #use this factor to scale the deformation of modes # SC.visualizationSettings.contour.minValue = -0.0003 # SC.visualizationSettings.contour.maxValue =  0.0003 #simulationSettings.timeIntegration.verboseModeFile = 3 #simulationSettings.solutionSettings.coordinatesSolutionFileName = "solution/coordinatesSolution.txt" #SHOULD work with 0.9 as well #create animation: #simulationSettings.solutionSettings.recordImagesInterval = 0.002 #SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame" #SC.visualizationSettings.window.renderWindowSize=[1920,1080] #load last model view #press space to continue #press space to continue #safely close rendering window! #store model view for next simulation