Skip to content

Commit 450e0cc

Browse files
jgerstmayrjgerstmayr
committed
final adaptations for version 1.9.0 - sparse matrix conversion in FEM, warnings, and print leftovers
1 parent 435099f commit 450e0cc

File tree

5 files changed

+71
-55
lines changed

5 files changed

+71
-55
lines changed

‎main/pythonDev/TestModels/runTestExamples.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -329,7 +329,7 @@ def GetFileNames(dirPath, fileEnding=''):
329329
else:
330330
exu.Print(str(len(examplesFailed)) + ' Examples OUT OF '+ str(totalExamples) + ' FAILED: ')
331331
for ef in examplesFailed:
332-
exu.Print(' TestModel ' + ef + ' FAILED')
332+
exu.Print(' Example ' + ef + ' FAILED')
333333

334334
exu.Print('Skipped '+str(nSkipped)+' examples')
335335
exu.Print('******************************************')

‎main/pythonDev/exudyn/FEM.py‎

Lines changed: 59 additions & 43 deletions
Original file line numberDiff line numberDiff line change
@@ -132,7 +132,7 @@ def CSRtoRowsAndColumns(sparseMatrixCSR):
132132

133133

134134
warnedCSRtoScipySparseCSR = False #add warning if this function is used with old format!
135-
#**function: convert internal compressed CSR to scipy.sparse csr matrix
135+
#**function: DEPRECATED: convert internal compressed CSR to scipy.sparse csr matrix; should not be used and raises warning; use SparseTripletsToScipySparseCSR instead!
136136
def CSRtoScipySparseCSR(sparseMatrixCSR):
137137
if IsListOrArray(sparseMatrixCSR):
138138
global warnedCSRtoScipySparseCSR
@@ -149,6 +149,17 @@ def CSRtoScipySparseCSR(sparseMatrixCSR):
149149
CheckForSciPyMatrix(sparseMatrixCSR, 'CSRtoScipySparseCSR')
150150
return sparseMatrixCSR
151151

152+
#**function: convert list of sparse triplets (or numpy array with one sparse triplet per row) into scipy.sparse csr_matrix
153+
def SparseTripletsToScipySparseCSR(sparseTriplets):
154+
if IsListOrArray(sparseTriplets):
155+
if type(sparseTriplets) == list:
156+
sparseTriplets = np.array(sparseTriplets)
157+
158+
from scipy.sparse import csr_matrix
159+
X = csr_matrix((sparseTriplets[:,2],(sparseTriplets[:,0].astype(int),sparseTriplets[:,1].astype(int))))
160+
return X
161+
else:
162+
raise ValueError('SparseTripletsToScipySparseCSR: requires sparse triplets as input')
152163

153164
#**function: convert scipy.sparse csr matrix to internal compressed CSR
154165
def ScipySparseCSRtoCSR(scipyCSR):
@@ -806,10 +817,10 @@ def __init__(self, femInterface):
806817
self.nodeArray = femInterface.GetNodePositionsAsArray()
807818
self.trigList = femInterface.GetSurfaceTriangles()
808819

820+
#self.massMatrixCSR = CSRtoScipySparseCSR(femInterface.GetMassMatrix(sparse=True)) #for multiplications
809821
#stiffnessMatrixCSR = CSRtoScipySparseCSR(femInterface.GetStiffnessMatrix(sparse=True))
810-
self.massMatrixCSR = CSRtoScipySparseCSR(femInterface.GetMassMatrix(sparse=True)) #for multiplications
811-
self.stiffnessMatrixSparse = femInterface.GetStiffnessMatrix(sparse=True)
812822
self.massMatrixSparse = femInterface.GetMassMatrix(sparse=True)
823+
self.stiffnessMatrixSparse = femInterface.GetStiffnessMatrix(sparse=True)
813824

814825
#new coordinates:
815826
self.nNodes = len(self.nodeArray) #stored in nNodes x 3 np-array
@@ -826,12 +837,12 @@ def __init__(self, femInterface):
826837
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++
827838
#FFRFreduced constant matrices:
828839
self.Phit = np.kron(np.ones(self.nNodes),np.eye(3)).T
829-
self.PhitTM = self.Phit.T @ self.massMatrixCSR #LARGE MATRIX COMPUTATION
840+
self.PhitTM = self.Phit.T @ self.massMatrixSparse #LARGE MATRIX COMPUTATION
830841
self.xRef = self.nodeArray.flatten() #node reference values in single vector (can be added then to q[7:])
831842
self.xRefTilde = ComputeSkewMatrix(self.xRef) #rfTilde without q
832843

833844
#not needed, but may be interesting for checks:
834-
self.inertiaLocal = self.xRefTilde.T @ self.massMatrixCSR @ self.xRefTilde #LARGE MATRIX COMPUTATION
845+
self.inertiaLocal = self.xRefTilde.T @ self.massMatrixSparse @ self.xRefTilde #LARGE MATRIX COMPUTATION
835846

836847

837848
#**classFunction: add according nodes, objects and constraints for FFRF object to MainSystem mbs; only implemented for Euler parameters
@@ -871,16 +882,19 @@ def AddObjectFFRF(self, exu, mbs,
871882

872883
stiffnessMatrixMC = exu.MatrixContainer()
873884
#stiffnessMatrixMC.SetWithSparseMatrixCSR(self.nODE2FF, self.nODE2FF, self.stiffnessMatrixSparse,useDenseMatrix=False)
874-
stiffnessMatrixMC.SetWithDenseMatrix(CompressedRowSparseToDenseMatrix(self.stiffnessMatrixSparse),useDenseMatrix=False)
885+
#OLD: stiffnessMatrixMC.SetWithDenseMatrix(CompressedRowSparseToDenseMatrix(self.stiffnessMatrixSparse),useDenseMatrix=False)
886+
stiffnessMatrixMC.SetWithSparseMatrix(self.stiffnessMatrixSparse,useDenseMatrix=False)
875887

876888
massMatrixMC = exu.MatrixContainer()
877889
#massMatrixMC.SetWithSparseMatrixCSR(self.nODE2FF, self.nODE2FF, self.massMatrixSparse,useDenseMatrix=False)
878-
massMatrixMC.SetWithDenseMatrix(CompressedRowSparseToDenseMatrix(self.massMatrixSparse),useDenseMatrix=False)
890+
#OLD: massMatrixMC.SetWithDenseMatrix(CompressedRowSparseToDenseMatrix(self.massMatrixSparse),useDenseMatrix=False)
891+
massMatrixMC.SetWithSparseMatrix(self.massMatrixSparse,useDenseMatrix=False)
879892

880893
if (massProportionalDamping != 0 or massProportionalDamping != 0):
881894
dampingMatrixMC = exu.MatrixContainer()
882-
dampingMatrixMC.SetWithDenseMatrix(massProportionalDamping*CompressedRowSparseToDenseMatrix(self.massMatrixSparse)+
883-
stiffnessProportionalDamping*CompressedRowSparseToDenseMatrix(self.stiffnessMatrixSparse),useDenseMatrix=False)
895+
dampingMatrixMC.SetWithSparseMatrix(massProportionalDamping*self.massMatrixSparse+
896+
stiffnessProportionalDamping*self.stiffnessMatrixSparse,
897+
useDenseMatrix=False)
884898
else:
885899
dampingMatrixMC=[]
886900

@@ -907,7 +921,7 @@ def AddObjectFFRF(self, exu, mbs,
907921
mGroundCoordinate = mbs.AddMarker(eii.MarkerNodeCoordinates(nodeNumber = nGround)) #Ground node ==> no action
908922

909923
X1 = np.zeros((6, self.nODE2FFRF))
910-
X1rot = ComputeSkewMatrix(self.xRef).T @ self.massMatrixCSR
924+
X1rot = ComputeSkewMatrix(self.xRef).T @ self.massMatrixSparse
911925
X1[0:3,self.nODE2rigid:] = self.PhitTM
912926
X1[3:6,self.nODE2rigid:] = X1rot
913927
offset = np.zeros(6)
@@ -956,11 +970,11 @@ def UFforce(self, exu, mbs, t, q, q_t):
956970
fTrans = A @ (omega3Dtilde @ self.PhitTM @ rfTilde @ omega3D + 2*self.PhitTM @ cF_tTilde @ omega3D)
957971
force[0:self.dim3D] = fTrans
958972

959-
fRot = -G.T@(omega3Dtilde @ rfTilde.T @ self.massMatrixCSR @ rfTilde @ omega3D +
960-
2*rfTilde.T @ self.massMatrixCSR @ cF_tTilde @ omega3D)
973+
fRot = -G.T@(omega3Dtilde @ rfTilde.T @ self.massMatrixSparse @ rfTilde @ omega3D +
974+
2*rfTilde.T @ self.massMatrixSparse @ cF_tTilde @ omega3D)
961975
force[self.dim3D:self.nODE2rigid] = fRot
962976

963-
fFlex = -self.massMatrixCSR @ (omegaTilde2 @ rF + 2*(omegaTilde @ cF_t))
977+
fFlex = -self.massMatrixSparse @ (omegaTilde2 @ rF + 2*(omegaTilde @ cF_t))
964978
force[self.nODE2rigid:] = fFlex
965979

966980
#add gravity:
@@ -992,7 +1006,7 @@ def UFmassGenericODE2(self, exu, mbs, t, q, q_t):
9921006
Mnew[0:self.dim3D, self.nODE2rigid:] = Mtf
9931007
Mnew[self.nODE2rigid:, 0:self.dim3D] = Mtf.T
9941008
#Mrf:
995-
Mrf = -G.T @ rfTilde.T @ self.massMatrixCSR
1009+
Mrf = -G.T @ rfTilde.T @ self.massMatrixSparse
9961010
Mnew[self.dim3D:self.dim3D+self.nODE2rot, self.nODE2rigid:] = Mrf
9971011
Mnew[self.nODE2rigid:, self.dim3D:self.dim3D+self.nODE2rot] = Mrf.T
9981012
#Mrr:
@@ -1091,12 +1105,14 @@ def __init__(self, femInterface=None, rigidBodyNodeType = 'NodeType.RotationEule
10911105
self.trigList = femInterface.GetSurfaceTriangles()
10921106
self.postProcessingModes = femInterface.postProcessingModes
10931107

1094-
stiffnessMatrixCSR = CSRtoScipySparseCSR(femInterface.GetStiffnessMatrix(sparse=True))
1095-
massMatrixCSR = CSRtoScipySparseCSR(femInterface.GetMassMatrix(sparse=True))
1108+
# stiffnessMatrixCSR = CSRtoScipySparseCSR(femInterface.GetStiffnessMatrix(sparse=True))
1109+
# massMatrixCSR = CSRtoScipySparseCSR(femInterface.GetMassMatrix(sparse=True))
1110+
stiffnessMatrixSparse = femInterface.GetStiffnessMatrix(sparse=True)
1111+
massMatrixSparse = femInterface.GetMassMatrix(sparse=True)
10961112

10971113
#compute reduced mass and stiffness matrices, only flexible coordinates:
1098-
self.massMatrixReduced = self.modeBasis.T @ massMatrixCSR @ self.modeBasis #LARGE MATRIX COMPUTATION
1099-
self.stiffnessMatrixReduced = self.modeBasis.T @ stiffnessMatrixCSR @ self.modeBasis #LARGE MATRIX COMPUTATION
1114+
self.massMatrixReduced = self.modeBasis.T @ massMatrixSparse @ self.modeBasis #LARGE MATRIX COMPUTATION
1115+
self.stiffnessMatrixReduced = self.modeBasis.T @ stiffnessMatrixSparse @ self.modeBasis #LARGE MATRIX COMPUTATION
11001116
RoundMatrix(self.massMatrixReduced, roundMassMatrix*abs(self.massMatrixReduced).max()) #erase off-diagonal terms for higher efficiency ...
11011117
RoundMatrix(self.stiffnessMatrixReduced, roundStiffnessMatrix*abs(self.stiffnessMatrixReduced).max())
11021118

@@ -1121,10 +1137,10 @@ def __init__(self, femInterface=None, rigidBodyNodeType = 'NodeType.RotationEule
11211137
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++
11221138
#FFRFreduced constant matrices:
11231139
Phit = np.kron(np.ones(self.nNodes),np.eye(3)).T
1124-
PhitTM = Phit.T @ massMatrixCSR #LARGE MATRIX COMPUTATION
1140+
PhitTM = Phit.T @ massMatrixSparse #LARGE MATRIX COMPUTATION
11251141
self.xRef = nodeArray.flatten() #node reference values in single vector (can be added then to q[7:])
11261142
xRefTilde = ComputeSkewMatrix(self.xRef) #rfTilde without q
1127-
self.inertiaLocal = xRefTilde.T @ massMatrixCSR @ xRefTilde #LARGE MATRIX COMPUTATION
1143+
self.inertiaLocal = xRefTilde.T @ massMatrixSparse @ xRefTilde #LARGE MATRIX COMPUTATION
11281144

11291145
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++
11301146
#prepare CMS matrices (some already given):
@@ -1139,12 +1155,12 @@ def __init__(self, femInterface=None, rigidBodyNodeType = 'NodeType.RotationEule
11391155
self.PsiTilde = ComputeSkewMatrix(self.modeBasis)
11401156
self.chiU = 1./self.totalMass*(PhitTM @ self.xRef) #center of mass
11411157
self.chiUtilde = ComputeSkewMatrix(self.chiU)
1142-
self.mPsiTildePsi = self.PsiTilde.T @ massMatrixCSR @ self.modeBasis #LARGE MATRIX COMPUTATION
1143-
self.mPsiTildePsiTilde = self.PsiTilde.T @ massMatrixCSR @ self.PsiTilde#LARGE MATRIX COMPUTATION
1144-
self.mPhitTPsi = Phit.T @ massMatrixCSR @ self.modeBasis #LARGE MATRIX COMPUTATION
1145-
self.mPhitTPsiTilde = Phit.T @ massMatrixCSR @ self.PsiTilde #LARGE MATRIX COMPUTATION
1146-
self.mXRefTildePsi = xRefTilde.T @ massMatrixCSR @ self.modeBasis #LARGE MATRIX COMPUTATION
1147-
self.mXRefTildePsiTilde = xRefTilde.T @ massMatrixCSR @ self.PsiTilde #LARGE MATRIX COMPUTATION
1158+
self.mPsiTildePsi = self.PsiTilde.T @ massMatrixSparse @ self.modeBasis #LARGE MATRIX COMPUTATION
1159+
self.mPsiTildePsiTilde = self.PsiTilde.T @ massMatrixSparse @ self.PsiTilde#LARGE MATRIX COMPUTATION
1160+
self.mPhitTPsi = Phit.T @ massMatrixSparse @ self.modeBasis #LARGE MATRIX COMPUTATION
1161+
self.mPhitTPsiTilde = Phit.T @ massMatrixSparse @ self.PsiTilde #LARGE MATRIX COMPUTATION
1162+
self.mXRefTildePsi = xRefTilde.T @ massMatrixSparse @ self.modeBasis #LARGE MATRIX COMPUTATION
1163+
self.mXRefTildePsiTilde = xRefTilde.T @ massMatrixSparse @ self.PsiTilde #LARGE MATRIX COMPUTATION
11481164

11491165
#fill already parts of the mass matrix in order to avoid copying this constant data during user function calls:
11501166
self.massMatrixFFRFreduced[0:3,0:3] = self.Mtt
@@ -1698,8 +1714,8 @@ def LoadFromFile(self, fileName, forceVersion=None):
16981714
self.massMatrix = MK['massMatrix']
16991715
self.stiffnessMatrix = MK['stiffnessMatrix']
17001716
else:
1701-
self.massMatrix = CSRtoScipySparseCSR( np.load(f) )
1702-
self.stiffnessMatrix = CSRtoScipySparseCSR( np.load(f) )
1717+
self.massMatrix = SparseTripletsToScipySparseCSR( np.load(f) )
1718+
self.stiffnessMatrix = SparseTripletsToScipySparseCSR( np.load(f) )
17031719
self.surface = list(np.load(f, allow_pickle=True))
17041720
self.nodeSets = list(np.load(f, allow_pickle=True))
17051721
self.elementSets = list(np.load(f, allow_pickle=True))
@@ -1907,7 +1923,7 @@ def ReadMassMatrixFromAbaqus(self, fileName, type='SparseRowColumnValue'):
19071923
self.massMatrix[:,1] -= 1
19081924

19091925
if not useOldCSRformat:
1910-
self.massMatrix = CSRtoScipySparseCSR(self.massMatrix)
1926+
self.massMatrix = SparseTripletsToScipySparseCSR(self.massMatrix)
19111927

19121928
#**classFunction: read stiffness matrix from compressed row text format (exported from Abaqus)
19131929
def ReadStiffnessMatrixFromAbaqus(self, fileName, type='SparseRowColumnValue'):
@@ -1916,7 +1932,7 @@ def ReadStiffnessMatrixFromAbaqus(self, fileName, type='SparseRowColumnValue'):
19161932
self.stiffnessMatrix[:,1] -= 1
19171933

19181934
if not useOldCSRformat:
1919-
self.stiffnessMatrix = CSRtoScipySparseCSR(self.stiffnessMatrix)
1935+
self.stiffnessMatrix = SparseTripletsToScipySparseCSR(self.stiffnessMatrix)
19201936

19211937

19221938
#%%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -2067,8 +2083,8 @@ def Flip3D(v):
20672083

20682084
self.nodes = {'Position':nodes}
20692085
self.elements = [{'Name':'NGsolve','Tet4':elements}]
2070-
self.massMatrix = CSRtoScipySparseCSR(M1)
2071-
self.stiffnessMatrix = CSRtoScipySparseCSR(K1)
2086+
self.massMatrix = SparseTripletsToScipySparseCSR(M1)
2087+
self.stiffnessMatrix = SparseTripletsToScipySparseCSR(K1)
20722088
self.surface = [{'Name':'meshSurface','Trigs':trigList}]
20732089

20742090
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -2740,7 +2756,7 @@ def GetGyroscopicMatrix(self, rotationAxis=2, sparse=True):
27402756
if sparse:
27412757
from scipy import sparse
27422758
xBlock = sparse.kron(sparse.eye(nNodes), X) #create big block-diagonal matrix
2743-
G=np.dot(xBlock,CSRtoScipySparseCSR(self.massMatrix))
2759+
G=np.dot(xBlock,self.massMatrix)
27442760
else:
27452761
xBlock = np.kron(np.eye(nNodes), X) #create big block-diagonal matrix
27462762
G=np.dot(xBlock,CompressedRowSparseToDenseMatrix(self.massMatrix))
@@ -2815,9 +2831,9 @@ def CreateLinearFEMObjectGenericODE2(self, mbs, color=[0.9,0.4,0.4,1.]):
28152831

28162832
nRows = self.NumberOfCoordinates()
28172833
Mcsr = exu.MatrixContainer()
2818-
Mcsr.SetWithSparseMatrixCSR(nRows,nRows,self.GetMassMatrix(sparse=True), useDenseMatrix=False)
2834+
Mcsr.SetWithSparseMatrix(self.GetMassMatrix(sparse=True), nRows,nRows,useDenseMatrix=False)
28192835
Kcsr = exu.MatrixContainer()
2820-
Kcsr.SetWithSparseMatrixCSR(nRows,nRows,self.GetStiffnessMatrix(sparse=True), useDenseMatrix=False)
2836+
Kcsr.SetWithSparseMatrix(self.GetStiffnessMatrix(sparse=True), nRows,nRows,useDenseMatrix=False)
28212837

28222838
#now add generic body built from FEM model with mass and stiffness matrix (optional damping could be added):
28232839
oGenericODE2 = mbs.AddObject(eii.ObjectGenericODE2(nodeNumbers = allNodeList,
@@ -2967,8 +2983,8 @@ def ComputeEigenmodes(self, nModes, excludeRigidBodyModes = 0, useSparseSolver =
29672983
#sorted, sparse eigen vectors
29682984
from scipy.sparse.linalg import eigsh #eigh for symmetric matrices, positive definite
29692985

2970-
K = CSRtoScipySparseCSR(self.GetStiffnessMatrix(sparse=True))
2971-
M = CSRtoScipySparseCSR(self.GetMassMatrix(sparse=True))
2986+
K = self.GetStiffnessMatrix(sparse=True)
2987+
M = self.GetMassMatrix(sparse=True)
29722988
#optional, using shift-invert mode; DOES NOT WORK:
29732989
#guess for smallest eigenvalue:
29742990
#n=self.NumberOfCoordinates()
@@ -3104,8 +3120,8 @@ def ComputeHurtyCraigBamptonModes(self,
31043120
if useSparseSolver:
31053121
from scipy.sparse.linalg import eigsh, factorized #eigh for symmetric matrices, positive definite
31063122

3107-
K = CSRtoScipySparseCSR(self.GetStiffnessMatrix(sparse=True))
3108-
M = CSRtoScipySparseCSR(self.GetMassMatrix(sparse=True))
3123+
K = self.GetStiffnessMatrix(sparse=True)
3124+
M = self.GetMassMatrix(sparse=True)
31093125

31103126
else: #not recommended for more than 2000 nodes (6000 DOF)!
31113127
if computationMode == HCBstaticModeSelection.RBE3:
@@ -3755,8 +3771,8 @@ def ComputeCampbellDiagram(self, terminalFrequency, nEigenfrequencies=10, freque
37553771
#d(x) = [ ] * x
37563772
# [-Minv*K, -omega*Minv*G]
37573773

3758-
M = CSRtoScipySparseCSR(self.GetMassMatrix(sparse=True))
3759-
K = CSRtoScipySparseCSR(self.GetStiffnessMatrix(sparse=True))
3774+
M = self.GetMassMatrix(sparse=True)
3775+
K = self.GetStiffnessMatrix(sparse=True)
37603776
G = self.GetGyroscopicMatrix(rotationAxis=rotationAxis, sparse=True) #already in scypi sparse format
37613777

37623778
nODE = self.NumberOfCoordinates()
@@ -3885,7 +3901,7 @@ def ReadMassMatrixFromAnsys(self, fileName, dofMappingVectorFile, sparse=True, v
38853901

38863902
MapSparseMatrixIndices(self.massMatrix, sorting)
38873903
if not useOldCSRformat:
3888-
self.massMatrix = CSRtoScipySparseCSR(self.massMatrix)
3904+
self.massMatrix = SparseTripletsToScipySparseCSR(self.massMatrix)
38893905

38903906
#**classFunction: read stiffness matrix from CSV format (exported from Ansys)
38913907
def ReadStiffnessMatrixFromAnsys(self, fileName, dofMappingVectorFile, sparse=True, verbose=False):
@@ -3901,7 +3917,7 @@ def ReadStiffnessMatrixFromAnsys(self, fileName, dofMappingVectorFile, sparse=Tr
39013917

39023918
MapSparseMatrixIndices(self.stiffnessMatrix, sorting)
39033919
if not useOldCSRformat:
3904-
self.stiffnessMatrix = CSRtoScipySparseCSR(self.stiffnessMatrix)
3920+
self.stiffnessMatrix = SparseTripletsToScipySparseCSR(self.stiffnessMatrix)
39053921

39063922
#**classFunction: read nodal coordinates (exported from Ansys as .txt-File)
39073923
def ReadNodalCoordinatesFromAnsys(self, fileName, verbose=False):

‎main/src/Linalg/SymbolicUtilities.h‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -508,8 +508,8 @@ class PySymbolicUserFunction : public SymbolicFunction
508508
}
509509
else if (py::isinstance<py::int_>(pyItemIndex)) //MainSystem
510510
{
511-
Index itemIndex = py::cast<Index>(pyItemIndex);
512-
CHECKandTHROW(itemIndex == -1, "Symbolic::SymbolicUserFunction: invalid item type (must be ObjectIndex, LoadIndex or -1 for MainSystem)");
511+
//Index itemIndex = py::cast<Index>(pyItemIndex);
512+
CHECKandTHROW(py::cast<Index>(pyItemIndex) == -1, "Symbolic::SymbolicUserFunction: invalid item type (must be ObjectIndex, LoadIndex or -1 for MainSystem)");
513513
indexType = "None";
514514
itemTypeName = "MainSystem";
515515
itemNumber = -1;

0 commit comments

Comments
 (0)