Merge pull request #228 from robbievanleeuwen/warping_lg

Use lagrangian multiplier for warping function

Author: robbievanleeuwen

sectionproperties/analysis/section.py

@@ -371,7 +371,7 @@ def warping_analysis(progress=None):
                     ),
                     total=len(warping_section.elements),
                 )
-                (k, k_lg, f_torsion) = warping_section.assemble_torsion(
+                (k_lg, f_torsion) = warping_section.assemble_torsion(
                     progress=progress, task=task
                 )
 
@@ -383,7 +383,7 @@ def warping_analysis(progress=None):
                 )
                 progress.update(0, advance=1)
             else:
-                (k, k_lg, f_torsion) = warping_section.assemble_torsion()
+                (k_lg, f_torsion) = warping_section.assemble_torsion()
 
             # ILU decomposition of stiffness matrices
             def ilu_decomp(progress=None, task=None):
@@ -424,9 +424,9 @@ def ilu_decomp(progress=None, task=None):
             # solve for warping function
             def solve_warping():
                 if solver_type == "cgs":
-                    omega = solver.solve_cgs(k, f_torsion, k_precond)
+                    omega = solver.solve_cgs_lagrange(k_lg, f_torsion, m=k_lg_precond)
                 elif solver_type == "direct":
-                    omega = solver.solve_direct(k, f_torsion)
+                    omega = solver.solve_direct_lagrange(k_lg, f_torsion)
 
                 return omega
 
@@ -457,7 +457,7 @@ def solve_warping():
             self.section_props.j = (
                 self.section_props.ixx_c
                 + self.section_props.iyy_c
-                - omega.dot(k.dot(np.transpose(omega)))
+                - omega.dot(k_lg[:-1, :-1].dot(np.transpose(omega)))
             )
 
             # assemble shear function load vectors
@@ -929,25 +929,25 @@ def calculate_frame_properties(self, solver_type="direct"):
             el.coords[0, :] -= self.section_props.cx
             el.coords[1, :] -= self.section_props.cy
 
-        (k, _, f) = warping_section.assemble_torsion(lg=False)
+        (k_lg, f) = warping_section.assemble_torsion()
 
         # if the cgs method is used, perform ILU decomposition
         if solver_type == "cgs":
-            k_precond = linalg.LinearOperator(
-                (self.num_nodes, self.num_nodes), linalg.spilu(k).solve
+            k_lg_precond = linalg.LinearOperator(
+                (self.num_nodes + 1, self.num_nodes + 1), linalg.spilu(k_lg).solve
             )
 
         # solve for warping function
         if solver_type == "cgs":
-            omega = solver.solve_cgs(k, f, k_precond)
+            omega = solver.solve_cgs_lagrange(k_lg, f, m=k_lg_precond)
         elif solver_type == "direct":
-            omega = solver.solve_direct(k, f)
+            omega = solver.solve_direct_lagrange(k_lg, f)
 
         # calculate the torsion constant
         self.section_props.j = (
             self.section_props.ixx_c
             + self.section_props.iyy_c
-            - omega.dot(k.dot(np.transpose(omega)))
+            - omega.dot(k_lg[:-1, :-1].dot(np.transpose(omega)))
         )
 
         return (
@@ -1223,18 +1223,15 @@ def calc_stress(progress=None):
         # return the stress_post object
         return stress_post
 
-    def assemble_torsion(self, lg=True, progress=None, task=None):
-        """Assembles stiffness matrices to be used for the computation of warping properties and
-        the torsion load vector (f_torsion). Both a regular (k) and Lagrangian multiplier (k_lg)
-        stiffness matrix are returned. The stiffness matrices are assembled using the sparse COO
-        format and returned in the sparse CSC format.
+    def assemble_torsion(self, progress=None, task=None):
+        """Assembles stiffness matrices to be used for the computation of warping
+        properties and the torsion load vector (f_torsion). A Lagrangian multiplier
+        (k_lg) stiffness matrix is returned. The stiffness matrix are assembled using
+        the sparse COO format and returned in the sparse CSC format.
 
-        :param bool lg: Whether or not to calculate the Lagrangian multiplier stiffness matrix
-
-        :return: Regular stiffness matrix, Lagrangian multiplier stiffness matrix and torsion load
-            vector *(k, k_lg, f_torsion)*
-        :rtype: tuple(:class:`scipy.sparse.csc_matrix`, :class:`scipy.sparse.csc_matrix`,
-            :class:`numpy.ndarray`)
+        :return: Lagrangian multiplier stiffness matrix and torsion load
+            vector *(k_lg, f_torsion)*
+        :rtype: tuple(:class:`scipy.sparse.csc_matrix`, :class:`numpy.ndarray`)
         """
 
         # initialise variables
@@ -1272,11 +1269,6 @@ def assemble_torsion(self, lg=True, progress=None, task=None):
             if progress is not None:
                 progress.update(task, advance=1)
 
-        k = coo_matrix((data, (row, col)), shape=(N, N))
-
-        if not lg:
-            return (csc_matrix(k), None, f_torsion)
-
         # construct Lagrangian multiplier matrix:
         # column vector of ones
         row = np.hstack((row, range(N)))
@@ -1295,7 +1287,7 @@ def assemble_torsion(self, lg=True, progress=None, task=None):
 
         k_lg = coo_matrix((data, (row, col)), shape=(N + 1, N + 1))
 
-        return (csc_matrix(k), csc_matrix(k_lg), f_torsion)
+        return (csc_matrix(k_lg), f_torsion)
 
     def plot_mesh(
         self,