Made Beam class accept cross-section geometries

sympy/physics/continuum_mechanics/beam.py

@@ -13,6 +13,7 @@
 from sympy.integrals import integrate
 from sympy.series import limit
 from sympy.plotting import plot, PlotGrid
+from sympy.geometry.entity import GeometryEntity
 from sympy.external import import_module
 from sympy.utilities.decorator import doctest_depends_on
 from sympy import lambdify
@@ -80,29 +81,42 @@ def __init__(self, length, elastic_modulus, second_moment, variable=Symbol('x'),
 
         Parameters
         ==========
+
         length : Sympifyable
             A Symbol or value representing the Beam's length.
+
         elastic_modulus : Sympifyable
             A SymPy expression representing the Beam's Modulus of Elasticity.
             It is a measure of the stiffness of the Beam material. It can
             also be a continuous function of position along the beam.
-        second_moment : Sympifyable
-            A SymPy expression representing the Beam's Second moment of area.
-            It is a geometrical property of an area which reflects how its
-            points are distributed with respect to its neutral axis. It can
-            also be a continuous function of position along the beam.
+
+        second_moment : Sympifyable or Geometry object
+            Describes the cross-section of the beam via a SymPy expression
+            representing the Beam's second moment of area. It is a geometrical
+            property of an area which reflects how its points are distributed
+            with respect to its neutral axis. It can also be a continuous
+            function of position along the beam. Alternatively ``second_moment``
+            can be a shape object such as a ``Polygon`` from the geometry module
+            representing the shape of the cross-section of the beam. The second moment
+            of area will be computed from the shape object internally.
+
         variable : Symbol, optional
             A Symbol object that will be used as the variable along the beam
             while representing the load, shear, moment, slope and deflection
             curve. By default, it is set to ``Symbol('x')``.
+
         base_char : String, optional
             A String that will be used as base character to generate sequential
             symbols for integration constants in cases where boundary conditions
             are not sufficient to solve them.
         """
         self.length = length
         self.elastic_modulus = elastic_modulus
-        self.second_moment = second_moment
+        if isinstance(second_moment, GeometryEntity):
+            self.cross_section = second_moment
+        else:
+            self.cross_section = None
+            self.second_moment = second_moment
         self.variable = variable
         self._base_char = base_char
         self._boundary_conditions = {'deflection': [], 'slope': []}
@@ -113,7 +127,8 @@ def __init__(self, length, elastic_modulus, second_moment, variable=Symbol('x'),
         self._hinge_position = None
 
     def __str__(self):
-        str_sol = 'Beam({}, {}, {})'.format(sstr(self._length), sstr(self._elastic_modulus), sstr(self._second_moment))
+        shape_description = self._cross_section if self._cross_section else self._second_moment
+        str_sol = 'Beam({}, {}, {})'.format(sstr(self._length), sstr(self._elastic_modulus), sstr(shape_description))
         return str_sol
 
     @property
@@ -180,7 +195,22 @@ def second_moment(self):
 
     @second_moment.setter
     def second_moment(self, i):
-        self._second_moment = sympify(i)
+        self._cross_section = None
+        if isinstance(i, GeometryEntity):
+            raise ValueError("To update cross-section geometry use `cross_section` attribute")
+        else:
+            self._second_moment = sympify(i)
+
+    @property
+    def cross_section(self):
+        """Cross-section of the beam"""
+        return self._cross_section
+
+    @cross_section.setter
+    def cross_section(self, s):
+        if s:
+            self._second_moment = s.second_moment_of_area()[0]
+        self._cross_section = s
 
     @property
     def boundary_conditions(self):

sympy/physics/continuum_mechanics/tests/test_beam.py

@@ -1,9 +1,10 @@
-from sympy import Symbol, symbols, S, simplify, Interval
+from sympy import Symbol, symbols, S, simplify, Interval, pi
 from sympy.physics.continuum_mechanics.beam import Beam
 from sympy.functions import SingularityFunction, Piecewise, meijerg, Abs, log
 from sympy.utilities.pytest import raises, slow
 from sympy.physics.units import meter, newton, kilo, giga, milli
 from sympy.physics.continuum_mechanics.beam import Beam3D
+from sympy.geometry import Circle, Polygon, Point2D, Triangle
 
 x = Symbol('x')
 y = Symbol('y')
@@ -607,3 +608,69 @@ def test_parabolic_loads():
     loading = beam.load.xreplace({L: 10, E: 20, I: 30, P: 40})
     assert loading.xreplace({x: 5}) == 40 * 5**8
     assert loading.xreplace({x: 15}) == 0
+
+
+def test_cross_section():
+    I = Symbol('I')
+    l = Symbol('l')
+    E = Symbol('E')
+    C3, C4 = symbols('C3, C4')
+    a, c, g, h, r, n = symbols('a, c, g, h, r, n')
+
+    # test for second_moment and cross_section setter
+    b0 = Beam(l, E, I)
+    assert b0.second_moment == I
+    assert b0.cross_section == None
+    b0.cross_section = Circle((0, 0), 5)
+    assert b0.second_moment == 625*pi/4
+    assert b0.cross_section == Circle((0, 0), 5)
+    b0.second_moment = 2*n - 6
+    assert b0.second_moment == 2*n-6
+    assert b0.cross_section == None
+    with raises(ValueError):
+        b0.second_moment = Circle((0, 0), 5)
+
+    # beam with a circular cross-section
+    b1 = Beam(50, E, Circle((0, 0), r))
+    assert b1.cross_section == Circle((0, 0), r)
+    assert b1.second_moment == pi*r*Abs(r)**3/4
+
+    b1.apply_load(-10, 0, -1)
+    b1.apply_load(R1, 5, -1)
+    b1.apply_load(R2, 50, -1)
+    b1.apply_load(90, 45, -2)
+    b1.solve_for_reaction_loads(R1, R2)
+    assert b1.load == (-10*SingularityFunction(x, 0, -1) + 82*SingularityFunction(x, 5, -1)/S(9)
+                         + 90*SingularityFunction(x, 45, -2) + 8*SingularityFunction(x, 50, -1)/9)
+    assert b1.bending_moment() == (-10*SingularityFunction(x, 0, 1) + 82*SingularityFunction(x, 5, 1)/9
+                                     + 90*SingularityFunction(x, 45, 0) + 8*SingularityFunction(x, 50, 1)/9)
+    q = (-5*SingularityFunction(x, 0, 2) + 41*SingularityFunction(x, 5, 2)/S(9)
+           + 90*SingularityFunction(x, 45, 1) + 4*SingularityFunction(x, 50, 2)/S(9))/(pi*E*r*Abs(r)**3)
+    assert b1.slope() == C3 + 4*q
+    q = (-5*SingularityFunction(x, 0, 3)/3 + 41*SingularityFunction(x, 5, 3)/27 + 45*SingularityFunction(x, 45, 2)
+           + 4*SingularityFunction(x, 50, 3)/27)/(pi*E*r*Abs(r)**3)
+    assert b1.deflection() == C3*x + C4 + 4*q
+
+    # beam with a recatangular cross-section
+    b2 = Beam(20, E, Polygon((0, 0), (a, 0), (a, c), (0, c)))
+    assert b2.cross_section == Polygon((0, 0), (a, 0), (a, c), (0, c))
+    assert b2.second_moment == a*c**3/12
+    # beam with a triangular cross-section
+    b3 = Beam(15, E, Triangle((0, 0), (g, 0), (g/2, h)))
+    assert b3.cross_section == Triangle(Point2D(0, 0), Point2D(g, 0), Point2D(g/2, h))
+    assert b3.second_moment == g*h**3/36
+
+    # composite beam
+    b = b2.join(b3, "fixed")
+    b.apply_load(-30, 0, -1)
+    b.apply_load(65, 0, -2)
+    b.apply_load(40, 0, -1)
+    b.bc_slope = [(0, 0)]
+    b.bc_deflection = [(0, 0)]
+
+    assert b.second_moment == Piecewise((a*c**3/12, x <= 20), (g*h**3/36, x <= 35))
+    assert b.cross_section == None
+    assert b.length == 35
+    assert b.slope().subs(x, 7) == 8400/(E*a*c**3)
+    assert b.slope().subs(x, 25) == 52200/(E*g*h**3) + 39600/(E*a*c**3)
+    assert b.deflection().subs(x, 30) == 537000/(E*g*h**3) + 712000/(E*a*c**3)