Skip to content

Commit

Permalink
Switch from transpose to inverse & other changes
Browse files Browse the repository at this point in the history
  • Loading branch information
JWock82 committed Oct 12, 2021
1 parent d4d8b51 commit 08a84f0
Show file tree
Hide file tree
Showing 5 changed files with 74 additions and 56 deletions.
33 changes: 23 additions & 10 deletions PyNite/FEModel3D.py
Original file line number Diff line number Diff line change
Expand Up @@ -2367,14 +2367,14 @@ def analyze_PDelta(self, log=False, max_iter=30, tol=0.01, sparse=True):

#%%
def _calc_reactions(self, log=False):
'''
Calculates reactions once the model is solved.
"""
Calculates reactions internally once the model is solved.
Parameters
----------
log : bool, optional
Prints updates to the console if set to True. Default is False.
'''
"""

# Print a status update to the console
if log: print('- Calculating reactions')
Expand All @@ -2394,12 +2394,8 @@ def _calc_reactions(self, log=False):
node.RxnMZ[combo.name] = 0.0

# Determine if the node has any supports
if (node.support_DX == True) \
or (node.support_DY == True) \
or (node.support_DZ == True) \
or (node.support_RX == True) \
or (node.support_RY == True) \
or (node.support_RZ == True):
if (node.support_DX or node.support_DY or node.support_DZ
or node.support_RX or node.support_RY or node.support_RZ):

# Sum the spring end forces at the node
for spring in self.Springs.values():
Expand Down Expand Up @@ -2522,7 +2518,7 @@ def _calc_reactions(self, log=False):
# Get the quad's global force matrix
# Storing it as a local variable eliminates the need to rebuild it every time a term is needed
quad_F = quad.F(combo.name)

node.RxnFX[combo.name] += quad_F[0, 0]
node.RxnFY[combo.name] += quad_F[1, 0]
node.RxnFZ[combo.name] += quad_F[2, 0]
Expand Down Expand Up @@ -2714,3 +2710,20 @@ def _check_statics(self):
# Print the static check table
print(statics_table)
print('')

#%%
def _orphaned_nodes(self):

# Step through each node in the model
for node in self.Nodes.values():

orphaned = True

# Check to see if the node is attached to any quads
quads = [quad.Name for quad in self.Quads.values() if quad.i_node == node or quad.j_node == node or quad.m_node == node or quad.n_node == node]

if quads != []:
orphaned = False

if orphaned == True:
print('Node ' + node.Name + ' is orphaned.')
6 changes: 3 additions & 3 deletions PyNite/Member3D.py
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
# %%
from numpy import array, zeros, matrix, transpose, add, subtract, matmul, insert, cross, divide
from numpy import array, zeros, matrix, add, subtract, matmul, insert, cross, divide
from numpy.linalg import inv
from math import isclose
from PyNite.BeamSegZ import BeamSegZ
Expand Down Expand Up @@ -475,14 +475,14 @@ def T(self):
def K(self):

# Calculate and return the stiffness matrix in global coordinates
return matmul(matmul(transpose(self.T()), self.k()), self.T())
return matmul(matmul(inv(self.T()), self.k()), self.T())

#%%
# Member global geometric stiffness matrix
def Kg(self, P=0):

# Calculate and return the geometric stiffness matrix in global coordinates
return matmul(matmul(transpose(self.T()), self.kg(P)), self.T())
return matmul(matmul(inv(self.T()), self.kg(P)), self.T())

#%%
def F(self, combo_name='Combo 1'):
Expand Down
6 changes: 3 additions & 3 deletions PyNite/Plate3D.py
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
from numpy import zeros, matrix, array, matmul, transpose, cross, add
from numpy import zeros, matrix, array, matmul, cross, add
from numpy.linalg import inv, norm
from PyNite.LoadCombo import LoadCombo

Expand Down Expand Up @@ -106,7 +106,7 @@ def k_b(self):
# rotational stiffnesses in the matrix.
k_rz = min(abs(k[1, 1]), abs(k[2, 2]), abs(k[4, 4]), abs(k[5, 5]),
abs(k[7, 7]), abs(k[8, 8]), abs(k[10, 10]), abs(k[11, 11])
)/1000
)/1000000

# The matrix currently only holds terms related to bending action. We need to expand it to
# with placeholders for all the degrees of freedom so it can be directly added to the
Expand Down Expand Up @@ -422,7 +422,7 @@ def K(self):
"""

# Calculate and return the stiffness matrix in global coordinates
return matmul(matmul(transpose(self.T()), self.k()), self.T())
return matmul(matmul(inv(self.T()), self.k()), self.T())

def FER(self, combo_name='Combo 1'):
"""
Expand Down
81 changes: 43 additions & 38 deletions PyNite/Quad3D.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
# 1. "Finite Element Procedures, 2nd Edition", Klaus-Jurgen Bathe
# 2. "A First Course in the Finite Element Method, 4th Edition", Daryl L. Logan

from numpy import array, arccos, dot, cross, matmul, add, zeros, transpose
from numpy import array, arccos, dot, cross, matmul, add, zeros
from numpy.linalg import inv, det, norm
from math import sin, cos
from PyNite.LoadCombo import LoadCombo
Expand Down Expand Up @@ -319,7 +319,7 @@ def k_b(self):
# Calculate the stiffness of a weak spring for the drilling degree of freedom (rotation about local z)
k_rz = min(abs(k[1, 1]), abs(k[2, 2]), abs(k[4, 4]), abs(k[5, 5]),
abs(k[7, 7]), abs(k[8, 8]), abs(k[10, 10]), abs(k[11, 11])
)/1000
)/1000000

# Initialize the expanded stiffness matrix to all zeros
k_exp = zeros((24, 24))
Expand Down Expand Up @@ -434,13 +434,13 @@ def k(self):
self._local_coords()

# Sum the bending and membrane stiffness matrices
return self.k_b() + self.k_m()
return add(self.k_b(), self.k_m())

#%%
def f(self, combo_name='Combo 1'):
'''
"""
Returns the quad element's local end force vector
'''
"""

# Calculate and return the plate's local end force vector
return add(matmul(self.k(), self.d(combo_name)), self.fer(combo_name))
Expand Down Expand Up @@ -514,18 +514,23 @@ def fer(self, combo_name='Combo 1'):

#%%
def d(self, combo_name='Combo 1'):
'''
"""
Returns the quad element's local displacement vector
'''
"""

# Calculate and return the local displacement vector
return matmul(self.T(), self.D(combo_name))

#%%
def F(self, combo_name='Combo 1'):
'''
"""
Returns the quad element's global force vector
'''
Parameters
----------
combo_name : string
The load combination to get results for.
"""

# Calculate and return the global force vector
return matmul(inv(self.T()), self.f(combo_name))
Expand All @@ -547,33 +552,33 @@ def D(self, combo_name='Combo 1'):
D = zeros((24, 1))

# Read in the global displacements from the nodes
D.itemset((0, 0), self.m_node.DX[combo_name])
D.itemset((1, 0), self.m_node.DY[combo_name])
D.itemset((2, 0), self.m_node.DZ[combo_name])
D.itemset((3, 0), self.m_node.RX[combo_name])
D.itemset((4, 0), self.m_node.RY[combo_name])
D.itemset((5, 0), self.m_node.RZ[combo_name])

D.itemset((6, 0), self.n_node.DX[combo_name])
D.itemset((7, 0), self.n_node.DY[combo_name])
D.itemset((8, 0), self.n_node.DZ[combo_name])
D.itemset((9, 0), self.n_node.RX[combo_name])
D.itemset((10, 0), self.n_node.RY[combo_name])
D.itemset((11, 0), self.n_node.RZ[combo_name])

D.itemset((12, 0), self.i_node.DX[combo_name])
D.itemset((13, 0), self.i_node.DY[combo_name])
D.itemset((14, 0), self.i_node.DZ[combo_name])
D.itemset((15, 0), self.i_node.RX[combo_name])
D.itemset((16, 0), self.i_node.RY[combo_name])
D.itemset((17, 0), self.i_node.RZ[combo_name])

D.itemset((18, 0), self.j_node.DX[combo_name])
D.itemset((19, 0), self.j_node.DY[combo_name])
D.itemset((20, 0), self.j_node.DZ[combo_name])
D.itemset((21, 0), self.j_node.RX[combo_name])
D.itemset((22, 0), self.j_node.RY[combo_name])
D.itemset((23, 0), self.j_node.RZ[combo_name])
D[0, 0] = self.m_node.DX[combo_name]
D[1, 0] = self.m_node.DY[combo_name]
D[2, 0] = self.m_node.DZ[combo_name]
D[3, 0] = self.m_node.RX[combo_name]
D[4, 0] = self.m_node.RY[combo_name]
D[5, 0] = self.m_node.RZ[combo_name]

D[6, 0] = self.n_node.DX[combo_name]
D[7, 0] = self.n_node.DY[combo_name]
D[8, 0] = self.n_node.DZ[combo_name]
D[9, 0] = self.n_node.RX[combo_name]
D[10, 0] = self.n_node.RY[combo_name]
D[11, 0] = self.n_node.RZ[combo_name]

D[12, 0] = self.i_node.DX[combo_name]
D[13, 0] = self.i_node.DY[combo_name]
D[14, 0] = self.i_node.DZ[combo_name]
D[15, 0] = self.i_node.RX[combo_name]
D[16, 0] = self.i_node.RY[combo_name]
D[17, 0] = self.i_node.RZ[combo_name]

D[18, 0] = self.j_node.DX[combo_name]
D[19, 0] = self.j_node.DY[combo_name]
D[20, 0] = self.j_node.DZ[combo_name]
D[21, 0] = self.j_node.RX[combo_name]
D[22, 0] = self.j_node.RY[combo_name]
D[23, 0] = self.j_node.RZ[combo_name]

# Return the global displacement vector
return D
Expand All @@ -584,11 +589,11 @@ def K(self):
Returns the quad element's global stiffness matrix
'''

# Get the transpose matrix
# Get the transformation matrix
T = self.T()

# Calculate and return the stiffness matrix in global coordinates
return matmul(matmul(transpose(T), self.k()), T)
return matmul(matmul(inv(T), self.k()), T)

#%%
# Global fixed end reaction vector
Expand Down
4 changes: 2 additions & 2 deletions PyNite/Spring3D.py
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
# %%
from numpy import zeros, array, transpose, add, subtract, matmul, insert, cross, divide
from numpy import zeros, array, add, subtract, matmul, insert, cross, divide
from numpy.linalg import inv
from math import isclose
import PyNite.FixedEndReactions
Expand Down Expand Up @@ -187,7 +187,7 @@ def K(self):
'''

# Calculate and return the stiffness matrix in global coordinates
return matmul(matmul(transpose(self.T()), self.k()), self.T())
return matmul(matmul(inv(self.T()), self.k()), self.T())

#%%
def F(self, combo_name='Combo 1'):
Expand Down

0 comments on commit 08a84f0

Please sign in to comment.