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. Note that since exudyn version 1.8.69 the mass and stiffness matrices in FEMinterface are either None or given in SciPy-sparse csr format (leading also to a new load/save fileVersion of FEMinterface)
- 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: OLD 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(…); the NEW format uses scipy’s internal sparse csr format! To switch to the old format, set exudyn.FEM.useOldCSRformat=True
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
- function description:compute diadic product of vector v and a 3D unit matrix = diadic(v,I
); 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 duplicatesvalue 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:DEPRECATED: convert internal compressed CSR to scipy.sparse csr matrix; should not be used and raises warning; use SparseTripletsToScipySparseCSR instead!
Relevant Examples (Ex) and TestModels (TM) with weblink to github:
shapeOptimization.py (Ex), ACFtest.py (TM)
Function: SparseTripletsToScipySparseCSR
SparseTripletsToScipySparseCSR(sparseTriplets
)
- function description:convert list of sparse triplets (or numpy array with one sparse triplet per row) into scipy.sparse csr_matrix
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 columnsneeded 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:
objectFFRFTest.py (TM)
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 OLD FEM internal compressed row sparse format; do not use!
Function: ReadMatrixFromAnsysMMF
ReadMatrixFromAnsysMMF(fileName
, verbose = False
)
- function description:This function reads either the mass or stiffness matrix from an AnsysMatrix Market Format (MMF). The corresponding matrix can either be exportedas 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 insidethe 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 insidethe 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.arraythe 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 selectionof the body whos mesh should to exported by choosing its geometry. Next,create a second named selcetion by using a worksheet. Add the named selectionthat was created first into the worksheet of the second named selection.Inside the working sheet, choose ‘convert’ and convert the first creatednamed selection to ‘mesh node’ (Netzknoten in german) and click on generateto create the second named selection. Next, right click on the secondnamed selection tha was created and choose ‘export’ and save the nodalcoordinates 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 selectionof the body whos mesh should to exported by choosing its geometry. Next,create a second named selcetion by using a worksheet. Add the named selectionthat was created first into the worksheet of the second named selection.Inside the worksheet, choose ‘convert’ and convert the first creatednamed selection to ‘mesh element’ (Netzelement in german) and click on generateto create the second named selection. Next, right click on the secondnamed selection tha was created and choose ‘export’ and save the elementsas .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 objectNumberobjectNumber
: number of ObjectFFRFreducedOrder in mbsoutputVariableType
: 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 nodesnodeNumberList
: 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:
netgenSTLtest.py (Ex), NGsolveCMStutorial.py (Ex)
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:
CMSexampleCourse.py (Ex), netgenSTLtest.py (Ex), NGsolveCMStutorial.py (Ex), NGsolveCraigBampton.py (Ex), NGsolvePistonEngine.py (Ex)
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 classinitializes the ObjectFFRFinterface with nodes, modes, surface description and systemmatrices from FEMinterfacedata 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 modulembs
: a MainSystem objectpositionRef
: 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 otherwiseconstrainRigidBodyMotion
: set True in order to add constraint (Tisserand frame) in order to suppress rigid motion of mesh nodescolor
: provided as list of 4 RGBA valuesadd 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:
objectFFRFTest2.py (TM)
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 classinitializes the ObjectFFRFreducedOrderInterface with nodes, modes, surface description and reduced system matrices from FEMinterfacedata 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 tresholdroundStiffnessMatrix
: 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 addedfileVersion
: 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 modulembs
: a MainSystem objectpositionRef
: 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 matrixinitialAngularVelocity
: 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 otherwiseUFforce
: (OPTIONAL, computation is slower) provide a user function, which computes the quadratic velocity vector and applied forces; see exampleUFmassMatrix
: (OPTIONAL, computation is slower) provide a user function, which computes the quadratic velocity vector and applied forces; see examplemassProportionalDamping
: Rayleigh damping factor for mass proportional damping (multiplied with reduced mass matrix), added to floating frame/modal coordinates onlystiffnessProportionalDamping
: 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 modulembs
: a MainSystem objectpositionRef
: 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 matrixinitialAngularVelocity
: 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 onlystiffnessProportionalDamping
: Rayleigh damping factor for stiffness proportional damping, added to floating frame/modal coordinates onlygravity
: set [0,0,0] if no gravity shall be applied, or to the gravity vector otherwisecolor
: provided as list of 4 RGBA values
Relevant Examples (Ex) and TestModels (TM) with weblink to github:
CMSexampleCourse.py (Ex), netgenSTLtest.py (Ex), NGsolveCMStutorial.py (Ex), NGsolveCraigBampton.py (Ex), NGsolvePistonEngine.py (Ex), abaqusImportTest.py (TM), NGsolveCrankShaftTest.py (TM), objectFFRFreducedOrderAccelerations.py (TM)
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 = None # np.array([[r0,c0,value0],[r1,c1,value1], ... ]) #currently only in SparseCSR format allowed!
fem.stiffnessMatrix= None # 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 = 2
)
- 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 addedfileVersion
: 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 addedforceVersion
: 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 allinstances or parts in order to use this function for import of flexible bodies for order reduction methods
- input:
fileName
: file name incl. pathtypeName
: this is what is searched for regarding nodes and elements, see your .inp filename
: if there are several parts, this name should address the according part nameverbose
: use True for some debug informationcreateSurfaceTrigs
: 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 createdngs.mesh
(NGsolve mesh, see examples)youngsModulus
: Young’s modulus used for mechanical modelpoissonsRatio
: Poisson’s ratio used for mechanical modeldensity
: density used for mechanical modelmeshOrder
: 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 storedexcludeRigidBodyModes
: 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=40verbose
: 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 consideredmaxEigensolveIterations
: maximum number of iterations for iterative eigensolver; default=40verbose
: 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 accuracyif 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 planeif 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 surfaceif 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 rusing a tolerance for the distance to the planeif 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 numbersnormalizeWeights
: 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 elementsstores the surface in self.surfaceonly 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 addedmesh
: a previously createdngs.mesh
(NGsolve mesh, see examples)youngsModulus
: Young’s modulus used for mechanical modelpoissonsRatio
: Poisson’s ratio used for mechanical modeldensity
: density used for mechanical modelmeshOrder
: 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 developerof 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 storedexcludeRigidBodyModes
: 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 modesnEigenModes
: prescribe the number of modes to be computed; only nEigenModes with smallest abs(eigenvalues) are considered and storeduseSparseSolver
: [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 considereduseSparseSolver
: for more than approx.~500 nodes, it is recommended to use the sparse solver; dense mode not available for RBE3computationMode
: 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 nodesboundaryNodesWeights
: 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 objectRBE3secondMomentOfAreaWeighting
: 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 usedverboseMode
: if True, some additional output is printedtimerTreshold
: 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 strainsoutputVariableType
: specify either exudyn.OutputVariableType.StressLocal or exudyn.OutputVariableType.StrainLocal as the desired output variablesnumberOfThreads
: 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 systemcreate a first order system Axd + Bx = 0 with x= [q,qd]’ and compute eigenvaluestakes mass M, stiffness K and gyroscopic matrix G from FEMinterfacecurrently only uses dense matrices, so it is limited to approx. 5000 unknowns!
- input:
terminalFrequency
: frequency in Hz, up to which the campbell diagram is computednEigenfrequencies
: gives the number of computed eigenfrequencies(modes), in addition to the rigid body mode 0frequencySteps
: gives the number of increments (gives frequencySteps+1 total points in campbell diagram)rotationAxis
:[0,1,2] = [x,y,z] provides rotation axisplotDiagram
: if True, plots diagram for nEigenfrequencies befor terminatingverbose
: if True, shows progress of computation; if verbose=2, prints also eigenfrequenciesuseCorotationalFrame
: 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 welluseSparseSolver
: 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 frequenciescampbellFrequencies
: 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:
CMSexampleCourse.py (Ex), netgenSTLtest.py (Ex), NGsolveCMStutorial.py (Ex), NGsolveCraigBampton.py (Ex), NGsolveGeometry.py (Ex), abaqusImportTest.py (TM), ACFtest.py (TM), compareAbaqusAnsysRotorEigenfrequencies.py (TM)