Merge pull request #27 from robbievanleeuwen/interaction_generalisation

Generalise moment interaction diagram

Author: robbievanleeuwen

concreteproperties/analysis_section.py

@@ -1,6 +1,7 @@
 from __future__ import annotations
 
 from dataclasses import dataclass
+from math import isinf
 from typing import TYPE_CHECKING, List, Optional, Tuple
 
 import numpy as np
@@ -161,7 +162,7 @@ def service_stress_analysis(
         theta: float,
         kappa: float,
         centroid: Tuple[float, float],
-    ) -> Tuple[float, float, float, float, float]:
+    ) -> Tuple[float, float, float, float]:
         r"""Performs a service stress analysis on the section.
 
         :param point_na: Point on the neutral axis
@@ -176,13 +177,12 @@ def service_stress_analysis(
 
         # initialise section actions
         n = 0
-        max_strain = 0
         m_x = 0
         m_y = 0
-        m_v = 0
+        max_strain = 0
 
         for el in self.elements:
-            el_n, el_m_x, el_m_y, el_m_v, el_max_strain = el.calculate_service_actions(
+            el_n, el_m_x, el_m_y, el_max_strain = el.calculate_service_actions(
                 point_na=point_na,
                 d_n=d_n,
                 theta=theta,
@@ -194,9 +194,8 @@ def service_stress_analysis(
             n += el_n
             m_x += el_m_x
             m_y += el_m_y
-            m_v += el_m_v
 
-        return n, m_x, m_y, m_v, max_strain
+        return n, m_x, m_y, max_strain
 
     def get_service_stress(
         self,
@@ -239,7 +238,7 @@ def get_service_stress(
             )
 
         # calculate total force
-        n, m_x, m_y, _, _ = self.service_stress_analysis(
+        n, m_x, m_y, _ = self.service_stress_analysis(
             point_na=point_na,
             d_n=d_n,
             theta=theta,
@@ -263,7 +262,7 @@ def ultimate_stress_analysis(
         d_n: float,
         theta: float,
         ultimate_strain: float,
-        pc: Tuple[float, float],
+        centroid: Tuple[float, float],
     ) -> Tuple[float, float, float]:
         r"""Performs an ultimate stress analysis on the section.
 
@@ -272,7 +271,7 @@ def ultimate_stress_analysis(
         :param theta: Angle (in radians) the neutral axis makes with the
             horizontal axis (:math:`-\pi \leq \theta \leq \pi`)
         :param ultimate_strain: Concrete strain at failure
-        :param pc: Location of the plastic centroid
+        :param centroid: Centroid about which to take moments
 
         :return: Axial force and resultant moments about the global axes
         """
@@ -288,7 +287,7 @@ def ultimate_stress_analysis(
                 d_n=d_n,
                 theta=theta,
                 ultimate_strain=ultimate_strain,
-                pc=pc,
+                centroid=centroid,
             )
 
             n += el_n
@@ -303,7 +302,7 @@ def get_ultimate_stress(
         point_na: Tuple[float, float],
         theta: float,
         ultimate_strain: float,
-        pc: Tuple[float, float],
+        centroid: Tuple[float, float],
     ) -> Tuple[np.ndarray, float, float, float]:
         r"""Given the neutral axis depth `d_n` and ultimate strain, determines the
         ultimate stresses with the section.
@@ -313,7 +312,7 @@ def get_ultimate_stress(
         :param theta: Angle (in radians) the neutral axis makes with the
             horizontal axis (:math:`-\pi \leq \theta \leq \pi`)
         :param ultimate_strain: Concrete strain at failure
-        :param pc: Location of the plastic centroid
+        :param centroid: Centroid about which to take moments
 
         :return: Ultimate stresses net force and distance from neutral axis to point of
             force action
@@ -325,13 +324,16 @@ def get_ultimate_stress(
         # loop through nodes
         for idx, node in enumerate(self.mesh_nodes):
             # get strain at node
-            strain = utils.get_ultimate_strain(
-                point=(node[0], node[1]),
-                point_na=point_na,
-                d_n=d_n,
-                theta=theta,
-                ultimate_strain=ultimate_strain,
-            )
+            if isinf(d_n):
+                strain = ultimate_strain
+            else:
+                strain = utils.get_ultimate_strain(
+                    point=(node[0], node[1]),
+                    point_na=point_na,
+                    d_n=d_n,
+                    theta=theta,
+                    ultimate_strain=ultimate_strain,
+                )
 
             # get stress at gauss point
             sig[idx] = self.geometry.material.ultimate_stress_strain_profile.get_stress(  # type: ignore
@@ -344,7 +346,7 @@ def get_ultimate_stress(
             d_n=d_n,
             theta=theta,
             ultimate_strain=ultimate_strain,
-            pc=pc,
+            centroid=centroid,
         )
 
         # calculate point of action
@@ -575,7 +577,7 @@ def calculate_service_actions(
         theta: float,
         kappa: float,
         centroid: Tuple[float, float],
-    ) -> Tuple[float, float, float, float, float]:
+    ) -> Tuple[float, float, float, float]:
         r"""Calculates service actions for the current finite element.
 
         :param point_na: Point on the neutral axis
@@ -590,10 +592,9 @@ def calculate_service_actions(
 
         # initialise element results
         force_e = 0
-        max_strain_e = 0
         m_x_e = 0
         m_y_e = 0
-        m_v_e = 0
+        max_strain_e = 0
 
         # get points for 1 point Gaussian integration
         gps = utils.gauss_points(n=1)
@@ -622,23 +623,20 @@ def calculate_service_actions(
             # calculate force (stress * area)
             force_gp = gp[0] * stress * j
 
-            # convert gauss point to local coordinates
-            u, _ = utils.global_to_local(theta=theta, x=x, y=y)
             # add force and moment
             force_e += force_gp
             m_x_e += force_e * (y - centroid[1])
             m_y_e += force_e * (x - centroid[0])
-            m_v_e += force_gp * u
 
-        return force_e, m_x_e, m_y_e, m_v_e, max_strain_e
+        return force_e, m_x_e, m_y_e, max_strain_e
 
     def calculate_ultimate_actions(
         self,
         point_na: Tuple[float, float],
         d_n: float,
         theta: float,
         ultimate_strain: float,
-        pc: Tuple[float, float],
+        centroid: Tuple[float, float],
     ) -> Tuple[float, float, float]:
         r"""Calculates ultimate actions for the current finite element.
 
@@ -647,7 +645,7 @@ def calculate_ultimate_actions(
         :param theta: Angle (in radians) the neutral axis makes with the
             horizontal axis (:math:`-\pi \leq \theta \leq \pi`)
         :param ultimate_strain: Concrete strain at failure
-        :param pc: Location of the plastic centroid
+        :param centroid: Centroid about which to take moments
 
         :return: Axial force and resultant moments about the global axes
         """
@@ -670,13 +668,16 @@ def calculate_ultimate_actions(
             y = np.dot(N, np.transpose(self.coords[1, :]))
 
             # get strain at gauss point
-            strain = utils.get_ultimate_strain(
-                point=(x, y),
-                point_na=point_na,
-                d_n=d_n,
-                theta=theta,
-                ultimate_strain=ultimate_strain,
-            )
+            if isinf(d_n):
+                strain = ultimate_strain
+            else:
+                strain = utils.get_ultimate_strain(
+                    point=(x, y),
+                    point_na=point_na,
+                    d_n=d_n,
+                    theta=theta,
+                    ultimate_strain=ultimate_strain,
+                )
 
             # get stress at gauss point
             stress = self.conc_material.ultimate_stress_strain_profile.get_stress(
@@ -688,7 +689,7 @@ def calculate_ultimate_actions(
 
             # add force and moment
             force_e += force_gp
-            m_x_e += force_gp * (y - pc[1])
-            m_y_e += force_gp * (x - pc[0])
+            m_x_e += force_gp * (y - centroid[1])
+            m_y_e += force_gp * (x - centroid[0])
 
         return force_e, m_x_e, m_y_e