Skip to content

Commit

Permalink
Optimize
Browse files Browse the repository at this point in the history
  • Loading branch information
@TLCFEM
TLCFEM committed Sep 29, 2023
1 parent b72b12e commit 84f5ad8
Show file tree
Hide file tree
Showing 2 changed files with 49 additions and 46 deletions.
40 changes: 16 additions & 24 deletions src/sectionproperties/analysis/fea.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -100,17 +100,14 @@ def geometric_properties(
# determine shape function, shape function derivative and jacobian
n, _, j = shape_function(coords=self.coords, gauss_point=gp)

nx, ny = self.coords @ n

area += gp[0] * j
qx += gp[0] * np.dot(n, np.transpose(self.coords[1, :])) * j
qy += gp[0] * np.dot(n, np.transpose(self.coords[0, :])) * j
ixx += gp[0] * np.dot(n, np.transpose(self.coords[1, :])) ** 2 * j
iyy += gp[0] * np.dot(n, np.transpose(self.coords[0, :])) ** 2 * j
ixy += (
gp[0]
* np.dot(n, np.transpose(self.coords[1, :]))
* np.dot(n, np.transpose(self.coords[0, :]))
* j
)
qx += gp[0] * ny * j
qy += gp[0] * nx * j
ixx += gp[0] * ny**2 * j
iyy += gp[0] * nx**2 * j
ixy += gp[0] * ny * nx * j

return (
area,
Expand Down Expand Up @@ -142,12 +139,11 @@ def torsion_properties(self) -> tuple[np.ndarray, np.ndarray]:
n, b, j = shape_function(coords=self.coords, gauss_point=gp)

# determine x and y position at Gauss point
nx = np.dot(n, np.transpose(self.coords[0, :]))
ny = np.dot(n, np.transpose(self.coords[1, :]))
nx, ny = self.coords @ n

# calculated modulus weighted stiffness matrix and load vector
k_el += (
gp[0] * np.dot(np.transpose(b), b) * j * (self.material.elastic_modulus)
gp[0] * np.dot(np.transpose(b), b) * j * self.material.elastic_modulus
)
f_el += (
gp[0]
Expand Down Expand Up @@ -188,8 +184,7 @@ def shear_load_vectors(
n, b, j = shape_function(coords=self.coords, gauss_point=gp)

# determine x and y position at Gauss point
nx = np.dot(n, np.transpose(self.coords[0, :]))
ny = np.dot(n, np.transpose(self.coords[1, :]))
nx, ny = self.coords @ n

# determine shear parameters
r = nx**2 - ny**2
Expand Down Expand Up @@ -260,9 +255,9 @@ def shear_warping_integrals(
n, _, j = shape_function(coords=self.coords, gauss_point=gp)

# determine x and y position at Gauss point
nx = np.dot(n, np.transpose(self.coords[0, :]))
ny = np.dot(n, np.transpose(self.coords[1, :]))
n_omega = np.dot(n, np.transpose(omega))
nx, ny = self.coords @ n

n_omega = np.dot(n, omega)

sc_xint += (
gp[0]
Expand Down Expand Up @@ -320,8 +315,7 @@ def shear_coefficients(
n, b, j = shape_function(coords=self.coords, gauss_point=gp)

# determine x and y position at Gauss point
nx = np.dot(n, np.transpose(self.coords[0, :]))
ny = np.dot(n, np.transpose(self.coords[1, :]))
nx, ny = self.coords @ n

# determine shear parameters
r = nx**2 - ny**2
Expand Down Expand Up @@ -385,8 +379,7 @@ def monosymmetry_integrals(
n, _, j = shape_function(coords=self.coords, gauss_point=gp)

# determine x and y position at Gauss point
nx = np.dot(n, np.transpose(self.coords[0, :]))
ny = np.dot(n, np.transpose(self.coords[1, :]))
nx, ny = self.coords @ n

# determine 11 and 22 position at Gauss point
nx_11, ny_22 = principal_coordinate(phi=phi, x=nx, y=ny)
Expand Down Expand Up @@ -519,8 +512,7 @@ def element_stress(
n_shape, b, _ = shape_function(coords=coords_c, gauss_point=gp)

# determine x and y position at Gauss point
nx = np.dot(n_shape, np.transpose(coords_c[0, :]))
ny = np.dot(n_shape, np.transpose(coords_c[1, :]))
nx, ny = coords_c @ n_shape

# determine 11 and 22 position at Gauss point
nx_11, ny_22 = principal_coordinate(phi=phi, x=nx, y=ny)
Expand Down
55 changes: 33 additions & 22 deletions src/sectionproperties/analysis/section.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1442,38 +1442,49 @@ def assemble_torsion(
# assemble the torsion load vector
f_torsion[el.node_ids] += f_el

# create row index vector
r = np.repeat(el.node_ids, n)

# create column index vector
c = np.tile(el.node_ids, n)

# flatten element stiffness matrix
k = k_el.flatten()

# add to global arrays
row = np.hstack((row, r))
col = np.hstack((col, c))
data = np.hstack((data, k))
# # create row index vector
# r = np.repeat(el.node_ids, n)
#
# # create column index vector
# c = np.tile(el.node_ids, n)
#
# # flatten element stiffness matrix
# k = k_el.flatten()
#
# # add to global arrays
# row = np.hstack((row, r))
# col = np.hstack((col, c))
# data = np.hstack((data, k))

for node_id in el.node_ids:
row.extend([node_id] * n)
col.extend(list(el.node_ids) * n)
data.extend(k_el.flatten())

if progress and task:
progress.update(task_id=task, advance=1)

# construct Lagrangian multiplier matrix:
# column vector of ones
row = np.hstack((row, range(n_size)))
col = np.hstack((col, np.repeat(n_size, n_size)))
data = np.hstack((data, np.repeat(1, n_size)))
row.extend(range(n_size))
col.extend([n_size] * n_size)
data.extend([1] * n_size)
# row = np.hstack((row, range(n_size)))
# col = np.hstack((col, np.repeat(n_size, n_size)))
# data = np.hstack((data, np.repeat(1, n_size)))

# row vector of ones
row = np.hstack((row, np.repeat(n_size, n_size)))
col = np.hstack((col, range(n_size)))
data = np.hstack((data, np.repeat(1, n_size)))
row.extend([n_size] * n_size)
col.extend(range(n_size))
data.extend([1] * n_size)
# row = np.hstack((row, np.repeat(n_size, n_size)))
# col = np.hstack((col, range(n_size)))
# data = np.hstack((data, np.repeat(1, n_size)))

# zero in bottom right corner
row = np.hstack((row, n_size))
col = np.hstack((col, n_size))
data = np.hstack((data, 0))
# row = np.hstack((row, n_size))
# col = np.hstack((col, n_size))
# data = np.hstack((data, 0))

k_lg = coo_matrix(
(data, (row, col)), shape=(n_size + 1, n_size + 1), dtype=float
Expand Down

0 comments on commit 84f5ad8

Please sign in to comment.