Module: FEM

Support functions and helper classes for import of meshes, finite element models (ABAQUS, ANSYS, NETGEN) and for generation of FFRF (floating frame of reference) objects.

  • Author:
    Johannes Gerstmayr; Stefan Holzinger (Abaqus and Ansys import utilities); Joachim Schöberl (support for Netgen and NGsolve import and eigen computations)

  • Date: 2020-03-10 (created)

  • Notes: internal CSR matrix storage format contains 3 float numbers per row: [row, column, value], can be converted to scipy csr sparse matrices with function CSRtoScipySparseCSR(…)

Function: CompressedRowSparseToDenseMatrix

CompressedRowSparseToDenseMatrix(sparseData)

  • function description:
    convert zero-based sparse matrix data to dense numpy matrix
  • input:
    sparseData: format (per row): [row, column, value] ==> converted into dense format
  • output:
    a dense matrix as np.array

Function: MapSparseMatrixIndices

MapSparseMatrixIndices(matrix, sorting)

  • function description:
    resort a sparse matrix (internal CSR format) with given sorting for rows and columns; changes matrix directly! used for ANSYS matrix import

Function: VectorDiadicUnitMatrix3D

VectorDiadicUnitMatrix3D(v)

  • function description:
    compute diadic product of vector v and a 3D unit matrix = diadic(v,I\(_{3x3}\)); used for ObjectFFRF and CMS implementation

Function: CyclicCompareReversed

CyclicCompareReversed(list1, list2)

  • function description:
    compare cyclic two lists, reverse second list; return True, if any cyclic shifted lists are same, False otherwise

Function: AddEntryToCompressedRowSparseArray

AddEntryToCompressedRowSparseArray(sparseData, row, column, value)

  • function description:
    add entry to compressedRowSparse matrix, avoiding duplicates
    value is either added to existing entry (avoid duplicates) or a new entry is appended

Function: CSRtoRowsAndColumns

CSRtoRowsAndColumns(sparseMatrixCSR)

  • function description:
    compute rows and columns of a compressed sparse matrix and return as tuple: (rows,columns)

Function: CSRtoScipySparseCSR

CSRtoScipySparseCSR(sparseMatrixCSR)

  • function description:
    convert internal compressed CSR to scipy.sparse csr matrix

Relevant Examples (Ex) and TestModels (TM) with weblink to github:

Function: ScipySparseCSRtoCSR

ScipySparseCSRtoCSR(scipyCSR)

  • function description:
    convert scipy.sparse csr matrix to internal compressed CSR

Function: ResortIndicesOfCSRmatrix

ResortIndicesOfCSRmatrix(mXXYYZZ, numberOfRows)

  • function description:
    resort indices of given NGsolve CSR matrix in XXXYYYZZZ format to XYZXYZXYZ format; numberOfRows must be equal to columns
    needed for import from NGsolve

Function: ResortIndicesOfNGvector

ResortIndicesOfNGvector(vXXYYZZ)

  • function description:
    resort indices of given NGsolve vector in XXXYYYZZZ format to XYZXYZXYZ format

Function: ResortIndicesExudyn2NGvector

ResortIndicesExudyn2NGvector(vXYZXYZ)

  • function description:
    resort indices of given Exudyun vector XYZXYZXYZ to NGsolve vector in XXXYYYZZZ format

Function: ConvertHexToTrigs

ConvertHexToTrigs(nodeNumbers)

  • function description:
    convert list of Hex8/C3D8 element with 8 nodes in nodeNumbers into triangle-List
  • notes:
    works for Hex20 elements, but does only take the corner nodes for drawing!

Relevant Examples (Ex) and TestModels (TM) with weblink to github:

Function: ConvertTetToTrigs

ConvertTetToTrigs(nodeNumbers)

  • function description:
    convert list of Tet4/Tet10 element with 4 or 10 nodes in nodeNumbers into triangle-List
  • notes:
    works for Tet10 elements, but does only take the corner nodes for drawing!

Function: ConvertDenseToCompressedRowMatrix

ConvertDenseToCompressedRowMatrix(denseMatrix)

  • function description:
    convert numpy.array dense matrix to (internal) compressed row format

Function: ReadMatrixFromAnsysMMF

ReadMatrixFromAnsysMMF(fileName, verbose = False)

  • function description:
    This function reads either the mass or stiffness matrix from an Ansys
    Matrix Market Format (MMF). The corresponding matrix can either be exported
    as dense matrix or sparse matrix.
  • input:
    fileName of MMF file
  • output:
    internal compressed row sparse matrix (as (nrows x 3) numpy array)
  • author:
    Stefan Holzinger
  • notes:
    A MMF file can be created in Ansys by placing the following APDL code inside
    the solution tree in Ansys Workbench:
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! APDL code that exports sparse stiffnes and mass matrix in MMF format. If
    ! the dense matrix is needed, replace *SMAT with *DMAT in the following
    ! APDL code.
    ! Export the stiffness matrix in MMF format
    *SMAT,MatKD,D,IMPORT,FULL,file.full,STIFF
    *EXPORT,MatKD,MMF,fileNameStiffnessMatrix,,,
    ! Export the mass matrix in MMF format
    *SMAT,MatMD,D,IMPORT,FULL,file.full,MASS
    *EXPORT,MatMD,MMF,fileNameMassMatrix,,,
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    In case a lumped mass matrix is needed, place the following APDL Code inside
    the Modal Analysis Tree:
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! APDL code to force Ansys to use a lumped mass formulation (if available for
    ! used elements)
    LUMPM, ON, , 0
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

Function: ReadMatrixDOFmappingVectorFromAnsysTxt

ReadMatrixDOFmappingVectorFromAnsysTxt(fileName)

  • function description:
    read sorting vector for ANSYS mass and stiffness matrices and return sorting vector as np.array
    the file contains sorting for nodes and applies this sorting to the DOF (assuming 3 DOF per node!)
    the resulting sorted vector is already converted to 0-based indices

Function: ReadNodalCoordinatesFromAnsysTxt

ReadNodalCoordinatesFromAnsysTxt(fileName, verbose = False)

  • function description:
    This function reads the nodal coordinates exported from Ansys.
  • input:
    fileName (file name ending must be .txt!)
  • output:
    nodal coordinates as numpy array
  • author:
    Stefan Holzinger
  • notes:
    The nodal coordinates can be exported from Ansys by creating a named selection
    of the body whos mesh should to exported by choosing its geometry. Next,
    create a second named selcetion by using a worksheet. Add the named selection
    that was created first into the worksheet of the second named selection.
    Inside the working sheet, choose ‘convert’ and convert the first created
    named selection to ‘mesh node’ (Netzknoten in german) and click on generate
    to create the second named selection. Next, right click on the second
    named selection tha was created and choose ‘export’ and save the nodal
    coordinates as .txt file.

Function: ReadElementsFromAnsysTxt

ReadElementsFromAnsysTxt(fileName, verbose = False)

  • function description:
    This function reads the nodal coordinates exported from Ansys.
  • input:
    fileName (file name ending must be .txt!)
  • output:
    element connectivity as numpy array
  • author:
    Stefan Holzinger
  • notes:
    The elements can be exported from Ansys by creating a named selection
    of the body whos mesh should to exported by choosing its geometry. Next,
    create a second named selcetion by using a worksheet. Add the named selection
    that was created first into the worksheet of the second named selection.
    Inside the worksheet, choose ‘convert’ and convert the first created
    named selection to ‘mesh element’ (Netzelement in german) and click on generate
    to create the second named selection. Next, right click on the second
    named selection tha was created and choose ‘export’ and save the elements
    as .txt file.

Function: CMSObjectComputeNorm

CMSObjectComputeNorm(mbs, objectNumber, outputVariableType, norm = 'max', nodeNumberList = [])

  • function description:
    compute current (max, min, …) value for chosen ObjectFFRFreducedOrder object (CMSobject) with exu.OutputVariableType. The function operates on nodal values. This is a helper function, which can be used to conveniently compute output quantities of the CMSobject efficiently and to use it in sensors
  • input:
    mbs: MainSystem of objectNumber
    objectNumber: number of ObjectFFRFreducedOrder in mbs
    outputVariableType: a exu.OutputVariableType out of [StressLocal, DisplacementLocal, VelocityLocal]
    norm: string containing chosen norm to be computed, out of ‘Mises’, ‘maxNorm’, ‘min’, ‘max’; ‘max’ will return maximum of all components (component wise), ‘min’ does same but for minimum; ‘maxNorm’ computes np.linalg.norm for every node and then takes maximum of all norms; Mises computes von-Mises stress for every node and then takes maximum of all nodes
    nodeNumberList: list of mesh node numbers (from FEMinterface); if empty [], all nodes are used; otherwise, only given nodes are evaluated
  • output:
    return value or list of values according to chosen norm as np.array

Relevant Examples (Ex) and TestModels (TM) with weblink to github:

CLASS MaterialBaseClass (in module FEM)

class description:

material base class, e.g., for FiniteElement

CLASS KirchhoffMaterial(MaterialBaseClass) (in module FEM)

class description:

class for representation of Kirchhoff (linear elastic, 3D and 2D) material

  • notes:
    use planeStress=False for plane strain

Class function: Strain2Stress

Strain2Stress(self, strain)

  • classFunction:
    convert strain tensor into stress tensor using elasticity tensor

Class function: StrainVector2StressVector

StrainVector2StressVector(self, strainVector)

  • classFunction:
    convert strain vector into stress vector

Class function: StrainVector2StressVector2D

StrainVector2StressVector2D(self, strainVector2D)

  • classFunction:
    compute 2D stress vector from strain vector

Class function: LameParameters

LameParameters(self)

  • classFunction:
    compute Lame parameters from internal Young’s modulus and Poisson ratio
  • output:
    return vector [mu, lam] of Lame parameters

Relevant Examples (Ex) and TestModels (TM) with weblink to github:

CLASS FiniteElement (in module FEM)

class description:

finite element base class for lateron implementations of other finite elements

CLASS Tet4(FiniteElement) (in module FEM)

class description:

simplistic 4-noded tetrahedral interface to compute strain/stress at nodal points

CLASS ObjectFFRFinterface (in module FEM)

class description:

compute terms necessary for ObjectFFRF class used internally in FEMinterface to compute ObjectFFRF object this class holds all data for ObjectFFRF user functions

Class function: __init__

__init__(self, femInterface)

  • classFunction:
    initialize ObjectFFRFinterface with FEMinterface class
    initializes the ObjectFFRFinterface with nodes, modes, surface description and systemmatrices from FEMinterface
    data is then transfered to mbs object with classFunction AddObjectFFRF(…)

Class function: AddObjectFFRF

AddObjectFFRF(self, exu, mbs, positionRef = [0,0,0], eulerParametersRef = [1,0,0,0], initialVelocity = [0,0,0], initialAngularVelocity = [0,0,0], gravity = [0,0,0], constrainRigidBodyMotion = True, massProportionalDamping = 0, stiffnessProportionalDamping = 0, color = [0.1,0.9,0.1,1.])

  • classFunction:
    add according nodes, objects and constraints for FFRF object to MainSystem mbs; only implemented for Euler parameters
  • input:
    exu: the exudyn module
    mbs: a MainSystem object
    positionRef: reference position of created ObjectFFRF (set in rigid body node underlying to ObjectFFRF)
    eulerParametersRef: reference euler parameters of created ObjectFFRF (set in rigid body node underlying to ObjectFFRF)
    initialVelocity: initial velocity of created ObjectFFRF (set in rigid body node underlying to ObjectFFRF)
    initialAngularVelocity: initial angular velocity of created ObjectFFRF (set in rigid body node underlying to ObjectFFRF)
    gravity: set [0,0,0] if no gravity shall be applied, or to the gravity vector otherwise
    constrainRigidBodyMotion: set True in order to add constraint (Tisserand frame) in order to suppress rigid motion of mesh nodes
    color: provided as list of 4 RGBA values
    add object to mbs as well as according nodes

Class function: UFforce

UFforce(self, exu, mbs, t, q, q_t)

  • classFunction:
    optional forceUserFunction for ObjectFFRF (per default, this user function is ignored)

Class function: UFmassGenericODE2

UFmassGenericODE2(self, exu, mbs, t, q, q_t)

  • classFunction:
    optional massMatrixUserFunction for ObjectFFRF (per default, this user function is ignored)

Relevant Examples (Ex) and TestModels (TM) with weblink to github:

CLASS ObjectFFRFreducedOrderInterface (in module FEM)

class description:

compute terms necessary for ObjectFFRFreducedOrder class used internally in FEMinterface to compute ObjectFFRFreducedOrder dictionary this class holds all data for ObjectFFRFreducedOrder user functions

Class function: __init__

__init__(self, femInterface = None, rigidBodyNodeType = 'NodeType.RotationEulerParameters', roundMassMatrix = 1e-13, roundStiffnessMatrix = 1e-13)

  • classFunction:
    initialize ObjectFFRFreducedOrderInterface with FEMinterface class
    initializes the ObjectFFRFreducedOrderInterface with nodes, modes, surface description and reduced system matrices from FEMinterface
    data is then transfered to mbs object with classFunction AddObjectFFRFreducedOrderWithUserFunctions(…)
  • input:
    femInterface: must provide nodes, surfaceTriangles, modeBasis, massMatrix, stiffness; if femInterface=None, an empty ObjectFFRFreducedOrderInterface instance is created which may be used to load data with LoadFromFile()
    roundMassMatrix: use this value to set entries of reduced mass matrix to zero which are below the treshold
    roundStiffnessMatrix: use this value to set entries of reduced stiffness matrix to zero which are below the treshold

Class function: SaveToFile

SaveToFile(self, fileName, fileVersion = 1)

  • classFunction:
    save all data to a data filename; can be used to avoid loading femInterface and FE data
  • input:
    fileName: string for path and file name without ending ==> “.npy” will be added
    fileVersion: FOR EXPERTS: this allows to store in older format, will be recovered when loading; must be integer; version must by > 0; the default value will change in future!
  • output:
    stores file

Class function: LoadFromFile

LoadFromFile(self, fileName)

  • classFunction:
    load all data (nodes, elements, …) from a data filename previously stored with SaveToFile(…).
    this function is much faster than the text-based import functions
  • input:
    fileName: string for path and file name without ending ==> “.npy” will be added
  • output:
    loads data into fem (note that existing values are not overwritten!)

Class function: AddObjectFFRFreducedOrderWithUserFunctions

AddObjectFFRFreducedOrderWithUserFunctions(self, exu, mbs, positionRef = [0,0,0], initialVelocity = [0,0,0], rotationMatrixRef = [], initialAngularVelocity = [0,0,0], gravity = [0,0,0], UFforce = 0, UFmassMatrix = 0, massProportionalDamping = 0, stiffnessProportionalDamping = 0, color = [0.1,0.9,0.1,1.], eulerParametersRef = [])

  • classFunction:
    add according nodes, objects and constraints for ObjectFFRFreducedOrder object to MainSystem mbs; use this function with userfunctions=0 in order to use internal C++ functionality, which is approx. 10x faster; implementation of userfunctions also available for rotation vector (Lie group formulation), which needs further testing
  • input:
    exu: the exudyn module
    mbs: a MainSystem object
    positionRef: reference position of created ObjectFFRFreducedOrder (set in rigid body node underlying to ObjectFFRFreducedOrder)
    initialVelocity: initial velocity of created ObjectFFRFreducedOrder (set in rigid body node underlying to ObjectFFRFreducedOrder)
    rotationMatrixRef: reference rotation of created ObjectFFRFreducedOrder (set in rigid body node underlying to ObjectFFRFreducedOrder); if [], it becomes the unit matrix
    initialAngularVelocity: initial angular velocity of created ObjectFFRFreducedOrder (set in rigid body node underlying to ObjectFFRFreducedOrder)
    eulerParametersRef: DEPRECATED, use rotationParametersRef or rotationMatrixRef in future: reference euler parameters of created ObjectFFRFreducedOrder (set in rigid body node underlying to ObjectFFRFreducedOrder)
    gravity: set [0,0,0] if no gravity shall be applied, or to the gravity vector otherwise
    UFforce: (OPTIONAL, computation is slower) provide a user function, which computes the quadratic velocity vector and applied forces; see example
    UFmassMatrix: (OPTIONAL, computation is slower) provide a user function, which computes the quadratic velocity vector and applied forces; see example
    massProportionalDamping: Rayleigh damping factor for mass proportional damping (multiplied with reduced mass matrix), added to floating frame/modal coordinates only
    stiffnessProportionalDamping: Rayleigh damping factor for stiffness proportional damping, added to floating frame/modal coordinates only (multiplied with reduced stiffness matrix)
    color: provided as list of 4 RGBA values
  • example:
#example of a user function for forces:
def UFforceFFRFreducedOrder(mbs, t, itemIndex, qReduced, qReduced_t):
    return cms.UFforceFFRFreducedOrder(exu, mbs, t, qReduced, qReduced_t)
#example of a user function for mass matrix:
def UFmassFFRFreducedOrder(mbs, t, itemIndex, qReduced, qReduced\_t):
    return cms.UFmassFFRFreducedOrder(exu, mbs, t, qReduced, qReduced\_t)

Class function: UFmassFFRFreducedOrder

UFmassFFRFreducedOrder(self, exu, mbs, t, qReduced, qReduced_t)

  • classFunction:
    CMS mass matrix user function; qReduced and qReduced_t contain the coordiantes of the rigid body node and the modal coordinates in one vector!

Class function: UFforceFFRFreducedOrder

UFforceFFRFreducedOrder(self, exu, mbs, t, qReduced, qReduced_t)

  • classFunction:
    CMS force matrix user function; qReduced and qReduced_t contain the coordiantes of the rigid body node and the modal coordinates in one vector!

Class function: AddObjectFFRFreducedOrder

AddObjectFFRFreducedOrder(self, mbs, positionRef = [0,0,0], initialVelocity = [0,0,0], rotationMatrixRef = [], initialAngularVelocity = [0,0,0], massProportionalDamping = 0, stiffnessProportionalDamping = 0, gravity = [0,0,0], color = [0.1,0.9,0.1,1.])

  • classFunction:
    add according nodes, objects and constraints for ObjectFFRFreducedOrder object to MainSystem mbs; use this function in order to use internal C++ functionality, which is approx. 10x faster than AddObjectFFRFreducedOrderWithUserFunctions(…)
  • input:
    exu: the exudyn module
    mbs: a MainSystem object
    positionRef: reference position of created ObjectFFRFreducedOrder (set in rigid body node underlying to ObjectFFRFreducedOrder)
    initialVelocity: initial velocity of created ObjectFFRFreducedOrder (set in rigid body node underlying to ObjectFFRFreducedOrder)
    rotationMatrixRef: reference rotation of created ObjectFFRFreducedOrder (set in rigid body node underlying to ObjectFFRFreducedOrder); if [], it becomes the unit matrix
    initialAngularVelocity: initial angular velocity of created ObjectFFRFreducedOrder (set in rigid body node underlying to ObjectFFRFreducedOrder)
    massProportionalDamping: Rayleigh damping factor for mass proportional damping, added to floating frame/modal coordinates only
    stiffnessProportionalDamping: Rayleigh damping factor for stiffness proportional damping, added to floating frame/modal coordinates only
    gravity: set [0,0,0] if no gravity shall be applied, or to the gravity vector otherwise
    color: provided as list of 4 RGBA values

Relevant Examples (Ex) and TestModels (TM) with weblink to github:

CLASS HCBstaticModeSelection(Enum) (in module FEM)

class description:

helper calss for function ComputeHurtyCraigBamptonModes, declaring some computation options. It offers the following options:

  • allBoundaryNodes: compute a single static mode for every boundary coordinate

  • RBE2: static modes only for rigid body motion at boundary nodes; using rigid boundary surfaces (additional stiffening)

  • RBE3: static modes only for rigid body motion at boundary nodes; averaged rigid body motion at boundary surfaces (leads to deformation at boundaries)

  • noStaticModes: do not compute static modes, only eigen modes (not recommended; usually only for tests)

CLASS FEMinterface (in module FEM)

class description:

general interface to different FEM / mesh imports and export to EXUDYN functions use this class to import meshes from different meshing or FEM programs (NETGEN/NGsolve , ABAQUS, ANSYS, ..) and store it in a unique format do mesh operations, compute eigenmodes and reduced basis, etc. load/store the data efficiently with LoadFromFile(…), SaveToFile(…) if import functions are slow export to EXUDYN objects

Class function: __init__

__init__(self)

  • classFunction:
    initalize all data of the FEMinterface by, e.g., fem = FEMinterface()
  • example:
#**** this is not an example, just a description for internal variables ****
#default values for member variables stored internally in FEMinterface fem and typical structure:
fem.nodes = {}                 # {'Position':[[x0,y0,z0],...], 'RigidBodyRxyz':[[x0,y0,z0],...],  },...]                     #dictionary of different node lists
fem.elements = []              # [{'Name':'identifier', 'Tet4':[[n0,n1,n2,n3],...], 'Hex8':[[n0,...,n7],...],  },...]        #there may be several element sets
fem.massMatrix = np.zeros((0,0))    # np.array([[r0,c0,value0],[r1,c1,value1], ... ])                                #currently only in SparseCSR format allowed!
fem.stiffnessMatrix=np.zeros((0,0)) # np.array([[r0,c0,value0],[r1,c1,value1], ... ])                                #currently only in SparseCSR format allowed!
fem.surface = []               # [{'Name':'identifier', 'Trigs':[[n0,n1,n2],...], 'Quads':[[n0,...,n3],...],  },...]           #surface with faces
fem.nodeSets = []              # [{'Name':'identifier', 'NodeNumbers':[n_0,...,n_ns], 'NodeWeights':[w_0,...,w_ns]},...]     #for boundary conditions, etc.
fem.elementSets = []           # [{'Name':'identifier', 'ElementNumbers':[n_0,...,n_ns]},...]                                #for different volumes, etc.
fem.modeBasis = {}             # {'matrix':[[Psi_00,Psi_01, ..., Psi_0m],...,[Psi_n0,Psi_n1, ..., Psi_nm]],'type':'NormalModes'} #'NormalModes' are eigenmodes, 'HCBmodes' are Craig-Bampton modes including static modes
fem.eigenValues = []           # [ev0, ev1, ...]                                                                             #eigenvalues according to eigenvectors in mode basis
fem.postProcessingModes = {}   # {'matrix':<matrix containing stress components (xx,yy,zz,yz,xz,xy) in each column, rows are for every mesh node>,'outputVariableType':exudyn.OutputVariableType.StressLocal}

Class function: SaveToFile

SaveToFile(self, fileName, fileVersion = 13)

  • classFunction:
    save all data (nodes, elements, …) to a data filename; this function is much faster than the text-based import functions
  • input:
    fileName: string for path and file name without ending ==> “.npy” will be added
    fileVersion: FOR EXPERTS: this allows to store in older format, will be recovered when loading; must be integer; version must by > 0
  • output:
    stores file

Class function: LoadFromFile

LoadFromFile(self, fileName, forceVersion = None)

  • classFunction:
    load all data (nodes, elements, …) from a data filename previously stored with SaveToFile(…).
    this function is much faster than the text-based import functions
  • input:
    fileName: string for path and file name without ending ==> “.npy” will be added
    forceVersion: FOR EXPERTS: this allows to store in older format, will be recovered when loading; must be integer; for old files, use forceVersion=0
  • output:
    loads data into fem (note that existing values are not overwritten!); returns file version or None if version is not available

Class function: ImportFromAbaqusInputFile

ImportFromAbaqusInputFile(self, fileName, typeName = 'Part', name = 'Part-1', verbose = False, createSurfaceTrigs = True, surfaceTrigsAll = False)

  • classFunction:
    import nodes and elements from Abaqus input file and create surface elements;
    node numbers in elements are converted from 1-based indices to python’s 0-based indices;
    This function can only import one part or instance; this means that you have to merge all
    instances or parts in order to use this function for import of flexible bodies for order reduction methods
  • input:
    fileName: file name incl. path
    typeName: this is what is searched for regarding nodes and elements, see your .inp file
    name: if there are several parts, this name should address the according part name
    verbose: use True for some debug information
    createSurfaceTrigs: if True, triangles are created for visualization (triangles both for Tet and Hex elements)
    surfaceTrigsAll: if False, visualization triangles are created at the surface; if True, surface triangles are created also for interior elements
  • output:
    return node numbers as numpy array
  • notes:
    only works for Hex8, Hex20, Tet4 and Tet10 (C3D4, C3D8, C3D8R, C3D10, C3D20, C3D20R) elements; some functionality is untested and works in limited cases; only works for one single part or instance

Class function: ReadMassMatrixFromAbaqus

ReadMassMatrixFromAbaqus(self, fileName, type = 'SparseRowColumnValue')

  • classFunction:
    read mass matrix from compressed row text format (exported from Abaqus); in order to export system matrices, write the following lines in your Abaqus input file:
    *STEP
    *MATRIX GENERATE, STIFFNESS, MASS
    *MATRIX OUTPUT, STIFFNESS, MASS, FORMAT=COORDINATE
    *End Step

Class function: ReadStiffnessMatrixFromAbaqus

ReadStiffnessMatrixFromAbaqus(self, fileName, type = 'SparseRowColumnValue')

  • classFunction:
    read stiffness matrix from compressed row text format (exported from Abaqus)

Class function: ImportMeshFromNGsolve

ImportMeshFromNGsolve(self, mesh, density, youngsModulus, poissonsRatio, verbose = False, computeEigenmodes = False, meshOrder = 1, **kwargs)

  • classFunction:
    import mesh from NETGEN/NGsolve and setup mechanical problem
  • input:
    mesh: a previously created ngs.mesh (NGsolve mesh, see examples)
    youngsModulus: Young’s modulus used for mechanical model
    poissonsRatio: Poisson’s ratio used for mechanical model
    density: density used for mechanical model
    meshOrder: use 1 for linear elements and 2 for second order elements (recommended to use 2 for much higher accuracy!)
    verbose: set True to print out some status information
  • output:
    creates according nodes, elements, in FEM and returns [bfM, bfK, fes] which are the (mass matrix M, stiffness matrix K) bilinear forms and the finite element space fes
  • author:
    Johannes Gerstmayr, Joachim Schöberl
  • notes:
    setting ngsolve.SetNumThreads(nt) you can select the number of treads that are used for assemble or other functionality with NGsolve functionality

Class function: ComputeEigenmodesNGsolve

ComputeEigenmodesNGsolve(self, bfM, bfK, nModes, maxEigensolveIterations = 40, excludeRigidBodyModes = 0, verbose = False)

  • classFunction:
    compute nModes smallest eigenvalues and eigenmodes from mass and stiffnessMatrix; store mode vectors in modeBasis, but exclude a number of ‘excludeRigidBodyModes’ rigid body modes from modeBasis; uses scipy for solution of generalized eigenvalue problem
  • input:
    nModes: prescribe the number of modes to be computed; total computed modes are (nModes+excludeRigidBodyModes), but only nModes with smallest absolute eigenvalues are considered and stored
    excludeRigidBodyModes: if rigid body modes are expected (in case of free-free modes), then this number specifies the number of eigenmodes to be excluded in the stored basis (usually 6 modes in 3D)
    maxEigensolveIterations: maximum number of iterations for iterative eigensolver; default=40
    verbose: if True, output some relevant information during solving
  • output:
    eigenmodes are stored internally in FEMinterface as ‘modeBasis’ and eigenvalues as ‘eigenValues’
  • author:
    Johannes Gerstmayr, Joachim Schöberl

Class function: ComputeHurtyCraigBamptonModesNGsolve

ComputeHurtyCraigBamptonModesNGsolve(self, bfM, bfK, boundaryNodesList, nEigenModes, maxEigensolveIterations = 40, verbose = False)

  • classFunction:
    compute static and eigen modes based on Hurty-Craig-Bampton, for details see theory part Section Model order reduction and component mode synthesis. This function uses internal computational functionality of NGsolve and is often much faster than the scipy variant
  • input:
    bfM: bilinearform for mass matrix as retured in ImportMeshFromNGsolve(…)
    bfK: bilinearform for stiffness matrix as retured in ImportMeshFromNGsolve(…)
    boundaryNodesList: [nodeList0, nodeList1, …] a list of node lists, each of them representing a set of ‘Position’ nodes for which a rigid body interface (displacement/rotation and force/torque) is created; NOTE THAT boundary nodes may not overlap between the different node lists (no duplicated node indices!)
    nEigenModes: number of eigen modes in addition to static modes (may be zero for RBE2 computationMode); eigen modes are computed for the case where all rigid body motions at boundaries are fixed; only smallest nEigenModes absolute eigenvalues are considered
    maxEigensolveIterations: maximum number of iterations for iterative eigensolver; default=40
    verbose: if True, output some relevant information during solving
  • output:
    stores computed modes in self.modeBasis and abs(eigenvalues) in self.eigenValues
  • author:
    Johannes Gerstmayr, Joachim Schöberl

Class function: ComputePostProcessingModesNGsolve

ComputePostProcessingModesNGsolve(self, fes, material = 0, outputVariableType = 'OutputVariableType.StressLocal', verbose = False)

  • classFunction:
    compute special stress or strain modes in order to enable visualization of stresses and strains in ObjectFFRFreducedOrder; takes a NGsolve fes as input and uses internal NGsolve methods to efficiently compute stresses or strains
  • input:
    fes: finite element space as retured in ImportMeshFromNGsolve(…)
    material: specify material properties for computation of stresses, using a material class, e.g. material = KirchhoffMaterial(Emodulus, nu, rho); not needed for strains (material = 0)
    outputVariableType: specify either exudyn.OutputVariableType.StressLocal or exudyn.OutputVariableType.StrainLocal as the desired output variables
  • output:
    post processing modes are stored in FEMinterface in local variable postProcessingModes as a dictionary, where ‘matrix’ represents the modes and ‘outputVariableType’ stores the type of mode as a OutputVariableType
  • author:
    Johannes Gerstmayr, Joachim Schöberl
  • notes:
    This function is implemented in Python and rather slow for larger meshes; for NGsolve / Netgen meshes, see the according ComputePostProcessingModesNGsolve function, which is usually much faster

Class function: GetMassMatrix

GetMassMatrix(self, sparse = True)

  • classFunction:
    get sparse mass matrix in according format

Class function: GetStiffnessMatrix

GetStiffnessMatrix(self, sparse = True)

  • classFunction:
    get sparse stiffness matrix in according format

Class function: NumberOfNodes

NumberOfNodes(self)

  • classFunction:
    get total number of nodes

Class function: GetNodePositionsAsArray

GetNodePositionsAsArray(self)

  • classFunction:
    get node points as array; only possible, if there exists only one type of Position nodes
  • notes:
    in order to obtain a list of certain node positions, see example
  • example:
p=GetNodePositionsAsArray(self)[42] #get node 42 position
nodeList=[1,13,42]
pArray=GetNodePositionsAsArray(self)[nodeList] #get np.array with positions of node indices

Class function: GetNodePositionsMean

GetNodePositionsMean(self, nodeNumberList)

  • classFunction:
    get mean (average) position of nodes defined by list of node numbers

Class function: NumberOfCoordinates

NumberOfCoordinates(self)

  • classFunction:
    get number of total nodal coordinates

Class function: GetNodeAtPoint

GetNodeAtPoint(self, point, tolerance = 1e-5, raiseException = True)

  • classFunction:
    get node number for node at given point, e.g. p=[0.1,0.5,-0.2], using a tolerance (+/-) if coordinates are available only with reduced accuracy
    if not found, it returns an invalid index

Class function: GetNodesInPlane

GetNodesInPlane(self, point, normal, tolerance = 1e-5)

  • classFunction:
    get node numbers in plane defined by point p and (normalized) normal vector n using a tolerance for the distance to the plane
    if not found, it returns an empty list

Class function: GetNodesInCube

GetNodesInCube(self, pMin, pMax)

  • classFunction:
    get node numbers in cube, given by pMin and pMax, containing the minimum and maximum x, y, and z coordinates
  • output:
    returns list of nodes; if no nodes found, return an empty list
  • example:
nList = GetNodesInCube([-1,-0.2,0],[1,0.5,0.5])

Class function: GetNodesOnLine

GetNodesOnLine(self, p1, p2, tolerance = 1e-5)

  • classFunction:
    get node numbers lying on line defined by points p1 and p2 and tolerance, which is accepted for points slightly outside the surface

Class function: GetNodesOnCylinder

GetNodesOnCylinder(self, p1, p2, radius, tolerance = 1e-5)

  • classFunction:
    get node numbers lying on cylinder surface; cylinder defined by cylinder axes (points p1 and p2),
    cylinder radius and tolerance, which is accepted for points slightly outside the surface
    if not found, it returns an empty list

Class function: GetNodesOnCircle

GetNodesOnCircle(self, point, normal, r, tolerance = 1e-5)

  • classFunction:
    get node numbers lying on a circle, by point p, (normalized) normal vector n (which is the axis of the circle) and radius r
    using a tolerance for the distance to the plane
    if not found, it returns an empty list

Class function: GetNodeWeightsFromSurfaceAreas

GetNodeWeightsFromSurfaceAreas(self, nodeList, normalizeWeights = True)

  • classFunction:
    return list of node weights based on surface triangle areas; surface triangles are identified as such for which all nodes of a triangle are on the surface
    **nodes: requires that surface triangles have been already built during import of finite element mesh, or by calling VolumeToSurfaceElements!
  • input:
    nodeList: list of local (Position) node numbers
    normalizeWeights: if True, weights are normalized to sum(weights)==1; otherwise, returned list contains areas according to nodes per
  • output:
    numpy array with weights according to indices in node list

Class function: GetSurfaceTriangles

GetSurfaceTriangles(self)

  • classFunction:
    return surface trigs as node number list (for drawing in EXUDYN and for node weights)

Class function: VolumeToSurfaceElements

VolumeToSurfaceElements(self, verbose = False)

  • classFunction:
    generate surface elements from volume elements
    stores the surface in self.surface
    only works for one element list and only for element types ‘Hex8’, ‘Hex20’, ‘Tet4’ and ‘Tet10’

Class function: GetGyroscopicMatrix

GetGyroscopicMatrix(self, rotationAxis = 2, sparse = True)

  • classFunction:
    get gyroscopic matrix in according format; rotationAxis=[0,1,2] = [x,y,z]

Class function: ScaleMassMatrix

ScaleMassMatrix(self, factor)

  • classFunction:
    scale (=multiply) mass matrix with factor

Class function: ScaleStiffnessMatrix

ScaleStiffnessMatrix(self, factor)

  • classFunction:
    scale (=multiply) stiffness matrix with factor

Class function: AddElasticSupportAtNode

AddElasticSupportAtNode(self, nodeNumber, springStiffness = [1e8,1e8,1e8])

  • classFunction:
    modify stiffness matrix to add elastic support (joint, etc.) to a node; nodeNumber zero based (as everywhere in the code…)
    springStiffness must have length according to the node size

Class function: AddNodeMass

AddNodeMass(self, nodeNumber, addedMass)

  • classFunction:
    modify mass matrix by adding a mass to a certain node, modifying directly the mass matrix

Class function: CreateLinearFEMObjectGenericODE2

CreateLinearFEMObjectGenericODE2(self, mbs, color = [0.9,0.4,0.4,1.])

  • classFunction:
    create GenericODE2 object out of (linear) FEM model; uses always the sparse matrix mode, independent of the solver settings; this model can be directly used inside the multibody system as a static or dynamic FEM subsystem undergoing small deformations; computation is several magnitudes slower than ObjectFFRFreducedOrder
  • input:
    mbs: multibody system to which the GenericODE2 is added
  • output:
    return list [oGenericODE2, nodeList] containing object number of GenericODE2 as well as the list of mbs node numbers of all NodePoint nodes

Class function: CreateNonlinearFEMObjectGenericODE2NGsolve

CreateNonlinearFEMObjectGenericODE2NGsolve(self, mbs, mesh, density, youngsModulus, poissonsRatio, meshOrder = 1, color = [0.9,0.4,0.4,1.])

  • classFunction:
    create GenericODE2 object fully nonlinear FEM model using NGsolve; uses always the sparse matrix mode, independent of the solver settings; this model can be directly used inside the multibody system as a static or dynamic nonlinear FEM subsystem undergoing large deformations; computation is several magnitudes slower than ObjectFFRFreducedOrder
  • input:
    mbs: multibody system to which the GenericODE2 is added
    mesh: a previously created ngs.mesh (NGsolve mesh, see examples)
    youngsModulus: Young’s modulus used for mechanical model
    poissonsRatio: Poisson’s ratio used for mechanical model
    density: density used for mechanical model
    meshOrder: use 1 for linear elements and 2 for second order elements (recommended to use 2 for much higher accuracy!)
  • output:
    return list [oGenericODE2, nodeList] containing object number of GenericODE2 as well as the list of mbs node numbers of all NodePoint nodes
  • author:
    Johannes Gerstmayr, Joachim Schöberl
  • notes:
    The interface to NETGEN/NGsolve has been created together with Joachim Schöberl, main developer
    of NETGEN/NGsolve ; Thank’s a lot!
    download NGsolve at: https://ngsolve.org/
    NGsolve needs Python 3.7 (64bit) ==> use according EXUDYN version!
    note that node/element indices in the NGsolve mesh are 1-based and need to be converted to 0-base!

Class function: ComputeEigenmodes

ComputeEigenmodes(self, nModes, excludeRigidBodyModes = 0, useSparseSolver = True)

  • classFunction:
    compute nModes smallest eigenvalues and eigenmodes from mass and stiffnessMatrix; store mode vectors in modeBasis, but exclude a number of ‘excludeRigidBodyModes’ rigid body modes from modeBasis; uses scipy for solution of generalized eigenvalue problem
  • input:
    nModes: prescribe the number of modes to be computed; total computed modes are (nModes+excludeRigidBodyModes), but only nModes with smallest absolute eigenvalues are considered and stored
    excludeRigidBodyModes: if rigid body modes are expected (in case of free-free modes), then this number specifies the number of eigenmodes to be excluded in the stored basis (usually 6 modes in 3D)
    useSparseSolver: for larger systems, the sparse solver needs to be used, which iteratively solves the problem and uses a random number generator (internally in ARPACK): therefore, results are not fully repeatable!!!
  • output:
    eigenmodes are stored internally in FEMinterface as ‘modeBasis’ and eigenvalues as ‘eigenValues’
  • notes:
    for NGsolve / Netgen meshes, see the according ComputeEigenmodesNGsolve function, which is usually much faster

Class function: ComputeEigenModesWithBoundaryNodes

ComputeEigenModesWithBoundaryNodes(self, boundaryNodes, nEigenModes, useSparseSolver = True)

  • classFunction:
    compute eigenmodes, using a set of boundary nodes that are all fixed; very similar to ComputeEigenmodes, but with additional definition of (fixed) boundary nodes.
  • input:
    boundaryNodes: a list of boundary node indices, refering to ‘Position’ type nodes in FEMinterface; all coordinates of these nodes are fixed for the computation of the modes
    nEigenModes: prescribe the number of modes to be computed; only nEigenModes with smallest abs(eigenvalues) are considered and stored
    useSparseSolver: [yet NOT IMPLEMENTED] for larger systems, the sparse solver needs to be used, which iteratively solves the problem and uses a random number generator (internally in ARPACK): therefore, results are not fully repeatable!!!
  • output:
    eigenmodes are stored internally in FEMinterface as ‘modeBasis’ and eigenvalues as ‘eigenValues’

Class function: ComputeHurtyCraigBamptonModes

ComputeHurtyCraigBamptonModes(self, boundaryNodesList, nEigenModes, useSparseSolver = True, computationMode = HCBstaticModeSelection.RBE2, boundaryNodesWeights = [], excludeRigidBodyMotion = True, RBE3secondMomentOfAreaWeighting = True, verboseMode = False, timerTreshold = 20000)

  • classFunction:
    compute static and eigen modes based on Hurty-Craig-Bampton, for details see theory part Section Model order reduction and component mode synthesis. Note that this function may need significant time, depending on your hardware, but 50.000 nodes will require approx. 1-2 minutes and more nodes typically raise time more than linearly.
  • input:
    boundaryNodesList: [nodeList0, nodeList1, …] a list of node lists, each of them representing a set of ‘Position’ nodes for which a rigid body interface (displacement/rotation and force/torque) is created; NOTE THAT boundary nodes may not overlap between the different node lists (no duplicated node indices!)
    nEigenModes: number of eigen modes in addition to static modes (may be zero for RBE2/RBE3 computationMode); eigen modes are computed for the case where all rigid body motions at boundaries are fixed; only smallest nEigenModes absolute eigenvalues are considered
    useSparseSolver: for more than approx.~500 nodes, it is recommended to use the sparse solver; dense mode not available for RBE3
    computationMode: see class HCBstaticModeSelection for available modes; select RBE2 / RBE3 as standard, which is both efficient and accurate and which uses rigid-body-interfaces (6 independent modes) per boundary; RBE3 mode uses singular value decomposition, which requires full matrices for boundary nodes; this becomes slow in particular if the number of a single boundary node set gets larger than 500 nodes
    boundaryNodesWeights: according list of weights with same order as boundaryNodesList, as returned e.g. by FEMinterface.GetNodeWeightsFromSurfaceAreas(…)
    excludeRigidBodyMotion: if True (recommended), the first set of boundary modes is eliminated, which defines the reference conditions for the FFRF object
    RBE3secondMomentOfAreaWeighting: if True, the weighting of RBE3 boundaries is done according to second moment of area; if False, the more conventional (but less appropriate) quadratic distance to reference point weighting is used
    verboseMode: if True, some additional output is printed
    timerTreshold: for more DOF than this number, CPU times are printed even with verboseMode=False
  • output:
    stores computed modes in self.modeBasis and abs(eigenvalues) in self.eigenValues
  • notes:
    for NGsolve / Netgen meshes, see the according ComputeHurtyCraigBamptonModesNGsolve function, which is usually much faster - currently only implemented for RBE2 case

Class function: GetEigenFrequenciesHz

GetEigenFrequenciesHz(self)

  • classFunction:
    return list of eigenvalues in Hz of previously computed eigenmodes

Class function: ComputePostProcessingModes

ComputePostProcessingModes(self, material = 0, outputVariableType = 'OutputVariableType.StressLocal', numberOfThreads = 1)

  • classFunction:
    compute special stress or strain modes in order to enable visualization of stresses and strains in ObjectFFRFreducedOrder;
  • input:
    material: specify material properties for computation of stresses, using a material class, e.g. material = KirchhoffMaterial(Emodulus, nu, rho); not needed for strains
    outputVariableType: specify either exudyn.OutputVariableType.StressLocal or exudyn.OutputVariableType.StrainLocal as the desired output variables
    numberOfThreads: if numberOfThreads=1, it uses single threaded computation; if numberOfThreads>1, it uses the multiprocessing pools functionality, which requires that all code in your main file must be encapsulated within an if clause “if __name__ == ‘__main__’:”, see examples; if numberOfThreads==-1, it uses all threads/CPUs available
  • output:
    post processing modes are stored in FEMinterface in local variable postProcessingModes as a dictionary, where ‘matrix’ represents the modes and ‘outputVariableType’ stores the type of mode as a OutputVariableType
  • notes:
    This function is implemented in Python and rather slow for larger meshes; for NGsolve / Netgen meshes, see the according ComputePostProcessingModesNGsolve function, which is usually much faster

Class function: ComputeCampbellDiagram

ComputeCampbellDiagram(self, terminalFrequency, nEigenfrequencies = 10, frequencySteps = 25, rotationAxis = 2, plotDiagram = False, verbose = False, useCorotationalFrame = False, useSparseSolver = False)

  • classFunction:
    compute Campbell diagram for given mechanical system
    create a first order system Axd + Bx = 0 with x= [q,qd]’ and compute eigenvalues
    takes mass M, stiffness K and gyroscopic matrix G from FEMinterface
    currently only uses dense matrices, so it is limited to approx. 5000 unknowns!
  • input:
    terminalFrequency: frequency in Hz, up to which the campbell diagram is computed
    nEigenfrequencies: gives the number of computed eigenfrequencies(modes), in addition to the rigid body mode 0
    frequencySteps: gives the number of increments (gives frequencySteps+1 total points in campbell diagram)
    rotationAxis:[0,1,2] = [x,y,z] provides rotation axis
    plotDiagram: if True, plots diagram for nEigenfrequencies befor terminating
    verbose: if True, shows progress of computation; if verbose=2, prints also eigenfrequencies
    useCorotationalFrame: if False, the classic rotor dynamics formulation for rotationally-symmetric rotors is used, where the rotor can be understood in a Lagrangian-Eulerian manner: the rotation is represented by an additional (Eulerian) velocity in rotation direction; if True, the corotational frame is used, which gives a factor 2 in the gyroscopic matrix and can be used for non-symmetric rotors as well
    useSparseSolver: for larger systems, the sparse solver needs to be used for creation of system matrices and for the eigenvalue solver (uses a random number generator internally in ARPACK, therefore, results are not fully repeatable!!!)
  • output:
    [listFrequencies, campbellFrequencies]
    listFrequencies: list of computed frequencies
    campbellFrequencies: array of campbell frequencies per eigenfrequency of system

Class function: CheckConsistency

CheckConsistency(self)

  • classFunction:
    perform some consistency checks

Class function: ReadMassMatrixFromAnsys

ReadMassMatrixFromAnsys(self, fileName, dofMappingVectorFile, sparse = True, verbose = False)

  • classFunction:
    read mass matrix from CSV format (exported from Ansys)

Class function: ReadStiffnessMatrixFromAnsys

ReadStiffnessMatrixFromAnsys(self, fileName, dofMappingVectorFile, sparse = True, verbose = False)

  • classFunction:
    read stiffness matrix from CSV format (exported from Ansys)

Class function: ReadNodalCoordinatesFromAnsys

ReadNodalCoordinatesFromAnsys(self, fileName, verbose = False)

  • classFunction:
    read nodal coordinates (exported from Ansys as .txt-File)

Class function: ReadElementsFromAnsys

ReadElementsFromAnsys(self, fileName, verbose = False)

  • classFunction:
    read elements (exported from Ansys as .txt-File)

Relevant Examples (Ex) and TestModels (TM) with weblink to github: