Add get_stress_at_points

Author: normanrichardson

sectionproperties/analysis/section.py

@@ -27,6 +27,8 @@
 import sectionproperties.analysis.solver as solver
 import sectionproperties.post.post as post
 
+from shapely.geometry import asPoint, Polygon
+from shapely.strtree import STRtree
 
 class Section:
     """Class for structural cross-sections.
@@ -2129,7 +2131,7 @@ def get_sf_p(self):
     def get_stress_at_point(
         self, pt: List[float], N=0, Mxx=0, Myy=0, M11=0, M22=0, Mzz=0, Vx=0, Vy=0
     ) -> Tuple[float]:
-        """Calaculates the stress at a point wiithin an element for given design actions
+        """Calaculates the stress at a point within an element for given design actions
         and returns *(sigma_zz, tau_xz, tau_yz)*
 
         :param pt: The point. A list of the x and y coordinate
@@ -2196,6 +2198,90 @@ def get_stress_at_point(
         return sigs
 
 
+    def get_stress_at_points(
+        self, pts: List[List[float]], N=0, Mxx=0, Myy=0, M11=0, M22=0, Mzz=0, Vx=0, Vy=0, agg_func=np.average
+    ) -> List[Tuple[float]]:
+        """Calaculates the stress at a set of points within an element for given design actions
+        and returns *(sigma_zz, tau_xz, tau_yz)*
+
+        :param pts: The points. A list of several x and y coordinates
+        :type pts: list[list[float, float]]
+        :param float N: Axial force
+        :param float Vx: Shear force acting in the x-direction
+        :param float Vy: Shear force acting in the y-direction
+        :param float Mxx: Bending moment about the centroidal xx-axis
+        :param float Myy: Bending moment about the centroidal yy-axis
+        :param float M11: Bending moment about the centroidal 11-axis
+        :param float M22: Bending moment about the centroidal 22-axis
+        :param float Mzz: Torsion moment about the centroidal zz-axis
+        :return: Resultant normal and shear stresses list[(sigma_zz, tau_xz, tau_yz)]
+        :rtype: list[tuple(float, float, float)]
+        """
+
+        action = {
+            "N": N,
+            "Mxx": Mxx,
+            "Myy": Myy,
+            "M11": M11,
+            "M22": M22,
+            "Mzz": Mzz,
+            "Vx": Vx,
+            "Vy": Vy,
+        }
+
+        sect_prop = {
+            "ea": self.section_props.ea,
+            "cx": self.section_props.cx,
+            "cy": self.section_props.cy,
+            "ixx": self.section_props.ixx_c,
+            "iyy": self.section_props.iyy_c,
+            "ixy": self.section_props.ixy_c,
+            "i11": self.section_props.i11_c,
+            "i22": self.section_props.i22_c,
+            "phi": self.section_props.phi,
+            "j": self.section_props.j,
+            "Delta_s": self.section_props.Delta_s,
+            "nu": self.section_props.nu_eff,
+        }
+
+        #create the search tree
+        p_mesh = [Polygon(self.geometry.mesh["vertices"][tri][0:3]) for tri in self.geometry.mesh["triangles"]]
+        s_tree = STRtree(p_mesh)
+        index_by_id = dict((id(poly), i) for i, poly in enumerate(p_mesh))
+
+        for pt in pts:
+            query_geom = asPoint(pt)
+            tri_ids = [index_by_id[id(poly)] for poly in s_tree.query(query_geom) if poly.intersects(query_geom)]
+            if len(tri_ids)==1:
+                tri = self.elements[tri_ids[0]]
+                sig = tri.local_element_stress(
+                    p=pt,
+                    **action,
+                    **sect_prop,
+                    omega=self.section_props.omega[tri.node_ids],
+                    psi_shear=self.section_props.psi_shear[tri.node_ids],
+                    phi_shear=self.section_props.phi_shear[tri.node_ids],
+                )
+            else:
+                sigs = []
+                for idx in tri_ids:
+                    tri = self.elements[idx]
+                    sigs.append(tri.local_element_stress(
+                            p=pt,
+                            **action,
+                            **sect_prop,
+                            omega=self.section_props.omega[tri.node_ids],
+                            psi_shear=self.section_props.psi_shear[tri.node_ids],
+                            phi_shear=self.section_props.phi_shear[tri.node_ids],
+                        )
+                    )
+                sig = (
+                    agg_func([sig[0] for sig in sigs]),
+                    agg_func([sig[1] for sig in sigs]),
+                    agg_func([sig[2] for sig in sigs]),
+                )
+            yield sig
+
 class PlasticSection:
     """Class for the plastic analysis of cross-sections.