Repository: sympy/sympy
File: sympy/physics/continuum_mechanics/column.py (Python)

Code before:
        self._integration_constants = dict(zip([C_N, C_u], integration_constants))

        self._is_solved = True

    def _solved_load(self):
        """
        Helper function that fills in the solved unknowns into
        the load equation.
        """
        solved_load = self._load
        solved_load = solved_load.subs(self._reaction_loads)
        solved_load = solved_load.subs(self._hinge_deflections)
        return solved_load

    def axial_force(self):
        """
        Returns a singularity function expression that
        represents the axial forces in the column.

        Examples
        ========
        A column with a length of 10 meters, elastic modulus E and
        area A is supported at both ends and in the middle. At x = 3
        the column is loaded by a point load of magnitude -1 kN. Between
        x = 5 and x = 10 the column is loaded by a uniform distributed
        load of -1 kN/m.

        >>> from sympy.physics.continuum_mechanics.column import Column
        >>> from sympy.core.symbol import symbols
        >>> E, A = symbols('E A')
        >>> c = Column(10, E, A)
        >>> c.apply_support(0)
        >>> c.apply_support(5)
        >>> c.apply_support(10)
        >>> c.apply_load(-5, 3, -1)
        >>> c.apply_load(-1, 5, 0, end=10)
        >>> c.axial_force()
        C_N - R_0*SingularityFunction(x, 0, 0) - R_10*SingularityFunction(x, 10, 0)
            - R_5*SingularityFunction(x, 5, 0) + 5*SingularityFunction(x, 3, 0)
            + SingularityFunction(x, 5, 1) - SingularityFunction(x, 10, 1)
        >>> c.solve_for_reaction_loads()
        >>> c.axial_force()
        -2*SingularityFunction(x, 0, 0) + 5*SingularityFunction(x, 3, 0)
            - 11*SingularityFunction(x, 5, 0)/2 + SingularityFunction(x, 5, 1)
            - 5*SingularityFunction(x, 10, 0)/2 - SingularityFunction(x, 10, 1)
        """
        load_equation = self._solved_load()
        x = self.variable
        C_N = self._integration_constants[Symbol('C_N')] if self._is_solved else Symbol('C_N')
        return -integrate(load_equation, x) + C_N


Review comment (refactor):
Duplicates between `reactions` and `self._applied_supports` are not handled now right?

Code after:
        {R_0: 31, R_10: 29}

        The same example when applied applied supports using manual method
        apply_load()

        >>> from sympy.physics.continuum_mechanics.column import Column
        >>> from sympy.core.symbol import symbols
        >>> E, A = symbols('E A')
        >>> R, G = symbols('R G')
        >>> c = Column(10, A, E)
        >>> c.apply_load(R,0,-1) # applied reaction load as a point load
        >>> c.apply_load(G,10,-1) # applied reaction load as a point load
        >>> c._bc_extension.append(0) # boundary conditions are added manually
        >>> c._bc_extension.append(10) # boundary conditions are added manually
        >>> c.apply_load(-5, 0, 0, end=10)
        >>> c.apply_load(-10, 4, -1)
        >>> c.solve_for_reaction_loads(R,G) # Reaction loads need to be passed.
        >>> c.reaction_loads
        {G: 29, R: 31}
        """
        x = self.variable
        qx = self._load
        L = self.length
        A = self.area
        E = self.elastic_modulus
        applied_supports = list(reactions) + self._applied_supports

        C_N, C_u = symbols('C_N, C_u')

        axial_force = -integrate(qx, x) + C_N
        extension = -(integrate(integrate(qx, x), x)) / (E * A) + C_N * x + C_u

        eq_axial_force = [
            limit(axial_force, x, 0, dir='-'),
            limit(axial_force, x, L, dir='+')
        ]

        eq_bc_displacement = [Eq(extension.subs(x, loc), 0) for loc in self._bc_extension]

        eq_bc_hinge = [Eq(limit(axial_force, x, loc, dir='+'), 0) for loc in self._bc_hinge] # Just right to avoid infinity

        total_eq = eq_axial_force + eq_bc_displacement + eq_bc_hinge
        unknowns = applied_supports + self._applied_hinges + [C_N, C_u]

        solution = list((linsolve(total_eq, unknowns).args)[0])
        solution = [nsimplify(s, rational=True) for s in solution] # To get rid of tiny residuals

        num_supports = len(applied_supports)
        reaction_solutions = solution[:num_supports]
        self._reaction_loads = dict(zip(applied_supports, reaction_solutions))