Update plastic centroid algorithm to search over entire section depth and remove pc_region

Author: robbievanleeuwen

docs/build/doctrees/environment.pickle

@@ -158,7 +158,7 @@ <h1>Source code for sectionproperties.pre.pre</h1><div class="highlight"><pre>
 <span class="kn">import</span> <span class="nn">meshpy.triangle</span> <span class="k">as</span> <span class="nn">triangle</span>
 
 
-<div class="viewcode-block" id="Material"><a class="viewcode-back" href="../../../rst/api.html#sectionproperties.pre.pre.Material">[docs]</a><span class="k">class</span> <span class="nc">Material</span><span class="p">:</span>
+<div class="viewcode-block" id="Material"><a class="viewcode-back" href="../../../rst/geom_mesh.html#sectionproperties.pre.pre.Material">[docs]</a><span class="k">class</span> <span class="nc">Material</span><span class="p">:</span>
     <span class="sd">&quot;&quot;&quot;Class for structural materials.</span>
 
 <span class="sd">    Provides a way of storing material properties related to a specific</span>

docs/build/doctrees/rst/analysis.doctree

@@ -403,7 +403,7 @@ <h2>Plastic Analysis<a class="headerlink" href="#plastic-analysis" title="Permal
 </div>
 <dl class="method">
 <dt>
-<code class="descclassname">CrossSection.</code><code class="descname">calculate_plastic_properties</code><span class="sig-paren">(</span><em>pc_region=[0.95, 0.95, 0.95, 0.95], time_info=False, verbose=False, debug=False</em><span class="sig-paren">)</span><a class="reference internal" href="../_modules/sectionproperties/analysis/cross_section.html#CrossSection.calculate_plastic_properties"><span class="viewcode-link">[source]</span></a></dt>
+<code class="descclassname">CrossSection.</code><code class="descname">calculate_plastic_properties</code><span class="sig-paren">(</span><em>time_info=False</em>, <em>verbose=False</em>, <em>debug=False</em><span class="sig-paren">)</span><a class="reference internal" href="../_modules/sectionproperties/analysis/cross_section.html#CrossSection.calculate_plastic_properties"><span class="viewcode-link">[source]</span></a></dt>
 <dd><p>Calculates the plastic properties of the cross-section and stores
 the, in the
 <a class="reference internal" href="api.html#sectionproperties.analysis.cross_section.SectionProperties" title="sectionproperties.analysis.cross_section.SectionProperties"><code class="xref py py-class docutils literal notranslate"><span class="pre">SectionProperties</span></code></a>
@@ -413,12 +413,6 @@ <h2>Plastic Analysis<a class="headerlink" href="#plastic-analysis" title="Permal
 <col class="field-body" />
 <tbody valign="top">
 <tr class="field-odd field"><th class="field-name">Parameters:</th><td class="field-body"><ul class="first last simple">
-<li><strong>pc_region</strong> (<em>list</em><em>[</em><em>float</em><em>, </em><em>float</em><em>, </em><em>float</em><em>, </em><em>float</em><em>]</em>) – Region of cross-section in which to search for the
-plastic centroids. The value represents the fraction of the depth
-to search within perpendicular to the axis being considered
-[x-axis, y-axis, 11-axis, 22-axis]. Ensure 0 &lt; pc_region[i] &lt; 1.
-The default scheme searches 95% of the section depth (ie. centroid
-<span class="math notranslate nohighlight">\(\pm\)</span> 47.5%)</li>
 <li><strong>time_info</strong> (<em>bool</em>) – If set to True, a detailed description of the
 computation and the time cost is printed to the terminal.</li>
 <li><strong>verbose</strong> (<em>bool</em>) – If set to True, the number of iterations required

docs/build/doctrees/rst/api.doctree

@@ -1728,7 +1728,7 @@ <h4>CrossSection Class<a class="headerlink" href="#crosssection-class" title="Pe
 
 <dl class="method">
 <dt id="sectionproperties.analysis.cross_section.CrossSection.calculate_plastic_properties">
-<code class="descname">calculate_plastic_properties</code><span class="sig-paren">(</span><em>pc_region=[0.95, 0.95, 0.95, 0.95], time_info=False, verbose=False, debug=False</em><span class="sig-paren">)</span><a class="reference internal" href="../_modules/sectionproperties/analysis/cross_section.html#CrossSection.calculate_plastic_properties"><span class="viewcode-link">[source]</span></a><a class="headerlink" href="#sectionproperties.analysis.cross_section.CrossSection.calculate_plastic_properties" title="Permalink to this definition">¶</a></dt>
+<code class="descname">calculate_plastic_properties</code><span class="sig-paren">(</span><em>time_info=False</em>, <em>verbose=False</em>, <em>debug=False</em><span class="sig-paren">)</span><a class="reference internal" href="../_modules/sectionproperties/analysis/cross_section.html#CrossSection.calculate_plastic_properties"><span class="viewcode-link">[source]</span></a><a class="headerlink" href="#sectionproperties.analysis.cross_section.CrossSection.calculate_plastic_properties" title="Permalink to this definition">¶</a></dt>
 <dd><p>Calculates the plastic properties of the cross-section and stores
 the, in the
 <a class="reference internal" href="#sectionproperties.analysis.cross_section.SectionProperties" title="sectionproperties.analysis.cross_section.SectionProperties"><code class="xref py py-class docutils literal notranslate"><span class="pre">SectionProperties</span></code></a>
@@ -1738,12 +1738,6 @@ <h4>CrossSection Class<a class="headerlink" href="#crosssection-class" title="Pe
 <col class="field-body" />
 <tbody valign="top">
 <tr class="field-odd field"><th class="field-name">Parameters:</th><td class="field-body"><ul class="first last simple">
-<li><strong>pc_region</strong> (<em>list</em><em>[</em><em>float</em><em>, </em><em>float</em><em>, </em><em>float</em><em>, </em><em>float</em><em>]</em>) – Region of cross-section in which to search for the
-plastic centroids. The value represents the fraction of the depth
-to search within perpendicular to the axis being considered
-[x-axis, y-axis, 11-axis, 22-axis]. Ensure 0 &lt; pc_region[i] &lt; 1.
-The default scheme searches 95% of the section depth (ie. centroid
-<span class="math notranslate nohighlight">\(\pm\)</span> 47.5%)</li>
 <li><strong>time_info</strong> (<em>bool</em>) – If set to True, a detailed description of the
 computation and the time cost is printed to the terminal.</li>
 <li><strong>verbose</strong> (<em>bool</em>) – If set to True, the number of iterations required
@@ -2724,7 +2718,7 @@ <h4>PlasticSection Class<a class="headerlink" href="#plasticsection-class" title
 
 <dl class="method">
 <dt id="sectionproperties.analysis.cross_section.PlasticSection.calculate_plastic_properties">
-<code class="descname">calculate_plastic_properties</code><span class="sig-paren">(</span><em>cross_section</em>, <em>pc_region</em>, <em>verbose</em><span class="sig-paren">)</span><a class="reference internal" href="../_modules/sectionproperties/analysis/cross_section.html#PlasticSection.calculate_plastic_properties"><span class="viewcode-link">[source]</span></a><a class="headerlink" href="#sectionproperties.analysis.cross_section.PlasticSection.calculate_plastic_properties" title="Permalink to this definition">¶</a></dt>
+<code class="descname">calculate_plastic_properties</code><span class="sig-paren">(</span><em>cross_section</em>, <em>verbose</em><span class="sig-paren">)</span><a class="reference internal" href="../_modules/sectionproperties/analysis/cross_section.html#PlasticSection.calculate_plastic_properties"><span class="viewcode-link">[source]</span></a><a class="headerlink" href="#sectionproperties.analysis.cross_section.PlasticSection.calculate_plastic_properties" title="Permalink to this definition">¶</a></dt>
 <dd><p>Calculates the location of the plastic centroid with respect to the
 centroidal and principal bending axes, the plastic section moduli and
 shape factors and stores the results to the supplied
@@ -2736,8 +2730,6 @@ <h4>PlasticSection Class<a class="headerlink" href="#plasticsection-class" title
 <tr class="field-odd field"><th class="field-name">Parameters:</th><td class="field-body"><ul class="first last simple">
 <li><strong>cross_section</strong> (<a class="reference internal" href="#sectionproperties.analysis.cross_section.CrossSection" title="sectionproperties.analysis.cross_section.CrossSection"><code class="xref py py-class docutils literal notranslate"><span class="pre">CrossSection</span></code></a>) – Cross section object that uses the same geometry
 and materials specified in the class constructor</li>
-<li><strong>pc_region</strong> (<em>list</em><em>[</em><em>float</em><em>, </em><em>float</em><em>, </em><em>float</em><em>, </em><em>float</em><em>]</em>) – Region of cross-section in which to search for the
-plastic centroids.</li>
 <li><strong>verbose</strong> (<em>bool</em>) – If set to True, the number of iterations required
 for each plastic axis is printed to the terminal.</li>
 </ul>
@@ -2843,10 +2835,10 @@ <h4>PlasticSection Class<a class="headerlink" href="#plasticsection-class" title
 
 <dl class="method">
 <dt id="sectionproperties.analysis.cross_section.PlasticSection.pc_algorithm">
-<code class="descname">pc_algorithm</code><span class="sig-paren">(</span><em>u</em>, <em>dlim</em>, <em>axis</em>, <em>pc_region</em>, <em>verbose</em><span class="sig-paren">)</span><a class="reference internal" href="../_modules/sectionproperties/analysis/cross_section.html#PlasticSection.pc_algorithm"><span class="viewcode-link">[source]</span></a><a class="headerlink" href="#sectionproperties.analysis.cross_section.PlasticSection.pc_algorithm" title="Permalink to this definition">¶</a></dt>
+<code class="descname">pc_algorithm</code><span class="sig-paren">(</span><em>u</em>, <em>dlim</em>, <em>axis</em>, <em>verbose</em><span class="sig-paren">)</span><a class="reference internal" href="../_modules/sectionproperties/analysis/cross_section.html#PlasticSection.pc_algorithm"><span class="viewcode-link">[source]</span></a><a class="headerlink" href="#sectionproperties.analysis.cross_section.PlasticSection.pc_algorithm" title="Permalink to this definition">¶</a></dt>
 <dd><p>An algorithm used for solving for the location of the plastic
 centroid. The algorithm searches for the location of the axis, defined
-by unit vector <em>u</em> and within limits <em>pc_region</em>, that satisfies force
+by unit vector <em>u</em> and within the section depth, that satisfies force
 equilibrium.</p>
 <table class="docutils field-list" frame="void" rules="none">
 <col class="field-name" />
@@ -2858,8 +2850,6 @@ <h4>PlasticSection Class<a class="headerlink" href="#plasticsection-class" title
 centroid to the extreme fibres perpendicular to the axis</li>
 <li><strong>axis</strong> (<em>int</em>) – The current axis direction: 1 (e.g. x or 11) or
 2 (e.g. y or 22)</li>
-<li><strong>pc_region</strong> (<em>list</em><em>[</em><em>float</em><em>, </em><em>float</em><em>, </em><em>float</em><em>, </em><em>float</em><em>]</em>) – Region of cross-section in which to search for the
-plastic centroids.</li>
 <li><strong>verbose</strong> (<em>bool</em>) – If set to True, the number of iterations required
 for each plastic axis is printed to the terminal.</li>
 </ul>

docs/build/html/_modules/sectionproperties/analysis/cross_section.html

@@ -694,21 +694,13 @@ def calculate_frame():
             self.section_props.iyy_c, self.section_props.ixy_c,
             self.section_props.j, self.section_props.phi)
 
-    def calculate_plastic_properties(self, pc_region=[0.95, 0.95, 0.95, 0.95],
-                                     time_info=False,
-                                     verbose=False, debug=False):
+    def calculate_plastic_properties(self, time_info=False, verbose=False,
+                                     debug=False):
         """Calculates the plastic properties of the cross-section and stores
         the, in the
         :class:`~sectionproperties.analysis.cross_section.SectionProperties`
         object contained in the ``section_props`` class variable.
 
-        :param pc_region: Region of cross-section in which to search for the
-            plastic centroids. The value represents the fraction of the depth
-            to search within perpendicular to the axis being considered
-            [x-axis, y-axis, 11-axis, 22-axis]. Ensure 0 < pc_region[i] < 1.
-            The default scheme searches 95% of the section depth (ie. centroid
-            :math:`\pm` 47.5%)
-        :type pc_region: list[float, float, float, float]
         :param bool time_info: If set to True, a detailed description of the
             computation and the time cost is printed to the terminal.
         :param bool verbose: If set to True, the number of iterations required
@@ -750,11 +742,11 @@ def calc_plastic():
 
             # calculate plastic properties
             try:
-                plastic_section.calculate_plastic_properties(
-                    self, pc_region, verbose)
+                plastic_section.calculate_plastic_properties(self, verbose)
             except ValueError:
                 str = "Plastic section properties calculation failed. "
-                str += "Consider increasing pc_region bounds."
+                str += "Contact robbie.vanleeuwen@gmail.com "
+                str += "with your analysis parameters."
                 raise RuntimeError(str)
 
         if time_info:
@@ -1815,7 +1807,7 @@ def calculate_centroid(self, elements):
 
         return (qy / ea, qx / ea)
 
-    def calculate_plastic_properties(self, cross_section, pc_region, verbose):
+    def calculate_plastic_properties(self, cross_section, verbose):
         """Calculates the location of the plastic centroid with respect to the
         centroidal and principal bending axes, the plastic section moduli and
         shape factors and stores the results to the supplied
@@ -1825,9 +1817,6 @@ def calculate_plastic_properties(self, cross_section, pc_region, verbose):
             and materials specified in the class constructor
         :type cross_section:
             :class:`~sectionproperties.analysis.cross_section.CrossSection`
-        :param pc_region: Region of cross-section in which to search for the
-            plastic centroids.
-        :type pc_region: list[float, float, float, float]
         :param bool verbose: If set to True, the number of iterations required
             for each plastic axis is printed to the terminal.
         """
@@ -1838,7 +1827,7 @@ def calculate_plastic_properties(self, cross_section, pc_region, verbose):
 
         # 1a) Calculate x-axis plastic centroid
         (y_pc, r, f, c_top, c_bot) = self.pc_algorithm(
-            np.array([1, 0]), fibres[2:], 1, pc_region[1], verbose)
+            np.array([1, 0]), fibres[2:], 1, verbose)
 
         self.check_convergence(r, 'x-axis')
         cross_section.section_props.y_pc = y_pc
@@ -1849,7 +1838,7 @@ def calculate_plastic_properties(self, cross_section, pc_region, verbose):
 
         # 1b) Calculate y-axis plastic centroid
         (x_pc, r, f, c_top, c_bot) = self.pc_algorithm(
-            np.array([0, 1]), fibres[0:2], 2, pc_region[0], verbose)
+            np.array([0, 1]), fibres[0:2], 2, verbose)
 
         self.check_convergence(r, 'y-axis')
         cross_section.section_props.x_pc = x_pc
@@ -1872,7 +1861,7 @@ def calculate_plastic_properties(self, cross_section, pc_region, verbose):
 
         # 2a) Calculate 11-axis plastic centroid
         (y22_pc, r, f, c_top, c_bot) = self.pc_algorithm(
-            ux, fibres[2:], 1, pc_region[3], verbose)
+            ux, fibres[2:], 1, verbose)
 
         # calculate the centroids in the principal coordinate system
         c_top_p = fea.principal_coordinate(
@@ -1889,7 +1878,7 @@ def calculate_plastic_properties(self, cross_section, pc_region, verbose):
 
         # 2b) Calculate 22-axis plastic centroid
         (x11_pc, r, f, c_top, c_bot) = self.pc_algorithm(
-            uy, fibres[0:2], 2, pc_region[2], verbose)
+            uy, fibres[0:2], 2, verbose)
 
         # calculate the centroids in the principal coordinate system
         c_top_p = fea.principal_coordinate(
@@ -2036,10 +2025,10 @@ def evaluate_force_eq(self, d, u, u_p, verbose):
         # return the force norm
         return f_norm
 
-    def pc_algorithm(self, u, dlim, axis, pc_region, verbose):
+    def pc_algorithm(self, u, dlim, axis, verbose):
         """An algorithm used for solving for the location of the plastic
         centroid. The algorithm searches for the location of the axis, defined
-        by unit vector *u* and within limits *pc_region*, that satisfies force
+        by unit vector *u* and within the section depth, that satisfies force
         equilibrium.
 
         :param u: Unit vector defining the direction of the axis
@@ -2049,9 +2038,6 @@ def pc_algorithm(self, u, dlim, axis, pc_region, verbose):
         :type dlim: list[float, float]
         :param int axis: The current axis direction: 1 (e.g. x or 11) or
             2 (e.g. y or 22)
-        :param pc_region: Region of cross-section in which to search for the
-            plastic centroids.
-        :type pc_region: list[float, float, float, float]
         :param bool verbose: If set to True, the number of iterations required
             for each plastic axis is printed to the terminal.
         :return: The distance to the plastic centroid axis *d*, the result
@@ -2062,17 +2048,14 @@ def pc_algorithm(self, u, dlim, axis, pc_region, verbose):
             list[float, float], list[float, float])
         """
 
-        D = dlim[1] - dlim[0]  # depth of the section perpendicular to the axis
-
         # calculate vector perpendicular to u
         if axis == 1:
             u_p = np.array([-u[1], u[0]])
         else:
             u_p = np.array([u[1], -u[0]])
 
-        x = (1 - pc_region) / 2
-        a = dlim[0] + x * D
-        b = dlim[1] - x * D
+        a = dlim[0]
+        b = dlim[1]
 
         (d, r) = brentq(self.evaluate_force_eq, a, b, args=(u, u_p, verbose),
                         full_output=True, disp=False, xtol=1e-6, rtol=1e-6)
@@ -2117,6 +2100,14 @@ def calculate_plastic_force(self, elements, u, p):
                 qx_bot += qx_el
                 qy_bot += qy_el
 
+        # if there are no elements in the top/bottom prevent division by zero
+        # N.B. the algorithm will never converge at this point, this is purely
+        # done to ensure a 100% search range
+        if ea_top == 0:
+            ea_top = 1
+        if ea_bot == 0:
+            ea_bot = 1
+
         # calculate the centroid of the top and bottom segments and save
         self.c_top = [qy_top / ea_top, qx_top / ea_top]
         self.c_bot = [qy_bot / ea_bot, qx_bot / ea_bot]
@@ -2220,6 +2211,11 @@ def add_line(self, geometry, line):
                     int_pts.append(new_pt)
                     fct_idx.append(idx)
 
+        # if less than 2 intersection points are found, we are at the edge of
+        # the section, therefore no line to add
+        if len(int_pts) < 2:
+            return
+
         # sort intersection points and facet list first by x, then by y
         int_pts = np.array(int_pts)
         idx_sort = np.lexsort((int_pts[:, 0], int_pts[:, 1]))