Repository: sympy/sympy

Issue:
Truss - compute forces using sympy functions rather math module
Feature request: the `Truss` class is a very nice addition. I think it would be important for it to compute the internal forces using symbolic functions instead of the math module.

For example, the current output is:

```python
from sympy import *
from sympy.physics.continuum_mechanics import Truss

t = Truss()
t.add_node("A", 0, 0)
t.add_node("B", 4, 0)
t.add_node("C", 8, 0)
t.add_node("D", 12, 0)
t.add_node("E", 4, 4)
t.add_node("F", 8, 4)
t.add_member("AB", "A", "B")
t.add_member("BC", "B", "C")
t.add_member("CD", "C", "D")
t.add_member("AE", "A", "E")
t.add_member("BF", "B", "F")
t.add_member("DF", "D", "F")
t.add_member("BE", "B", "E")
t.add_member("CF", "C", "F")
t.add_member("EF", "E", "F")
t.apply_load("B", magnitude=1, direction=270)
t.apply_support("A", type="pinned")
t.apply_support("D", type="roller")
t.solve()
t.internal_forces
# out:
# {'AB': 0.666666666666666,
#  'BC': 0.333333333333333,
#  'CD': 0.333333333333333,
#  'AE': -0.666666666666667*sqrt(2),
#  'BF': 0.333333333333333*sqrt(2),
#  'DF': -0.333333333333333*sqrt(2),
#  'BE': 0.666666666666667,
#  'CF': 0,
#  'EF': -0.666666666666667}
```

I'd rather have it like this:

```text
{'AB': 2/3,
 'BC': 1/3,
 'CD': 1/3,
 'AE': -2*sqrt(2)/3,
 'BF': sqrt(2)/3,
 'DF': -sqrt(2)/3,
 'BE': 2/3,
 'CF': 0,
 'EF': -2/3}
```

Switching from symbolic to numeric is actually simple.

Discussion:

Patch:
diff --git a/sympy/physics/continuum_mechanics/truss.py b/sympy/physics/continuum_mechanics/truss.py
index 5c2d0d0517..449ae52a7f 100644
--- a/sympy/physics/continuum_mechanics/truss.py
+++ b/sympy/physics/continuum_mechanics/truss.py
@@ -11,7 +11,7 @@
 from sympy import Matrix, pi
 from sympy.functions.elementary.miscellaneous import sqrt
 from sympy.matrices.dense import zeros
-import math
+from sympy import sin, cos
 
 
 
@@ -657,9 +657,9 @@ def solve(self):
         >>> t.apply_support("node_2", type="roller")
         >>> t.solve()
         >>> t.reaction_loads
-        {'R_node_1_x': 0, 'R_node_1_y': 6.66666666666667, 'R_node_2_y': 3.33333333333333}
+        {'R_node_1_x': 0, 'R_node_1_y': 20/3, 'R_node_2_y': 10/3}
         >>> t.internal_forces
-        {'member_1': 6.66666666666666, 'member_2': 6.66666666666667, 'member_3': -6.66666666666667*sqrt(2), 'member_4': -3.33333333333333*sqrt(5), 'member_5': 10.0}
+        {'member_1': 20/3, 'member_2': 20/3, 'member_3': -20*sqrt(2)/3, 'member_4': -10*sqrt(5)/3, 'member_5': 10}
         """
         count_reaction_loads = 0
         for node in self._nodes:
@@ -677,8 +677,8 @@ def solve(self):
             if node[0] in list(self._loads):
                 for load in self._loads[node[0]]:
                     if load[0]!=Symbol('R_'+str(node[0])+'_x') and load[0]!=Symbol('R_'+str(node[0])+'_y'):
-                        load_matrix[load_matrix_row] -= load[0]*math.cos(pi*load[1]/180)
-                        load_matrix[load_matrix_row + 1] -= load[0]*math.sin(pi*load[1]/180)
+                        load_matrix[load_matrix_row] -= load[0]*cos(pi*load[1]/180)
+                        load_matrix[load_matrix_row + 1] -= load[0]*sin(pi*load[1]/180)
             load_matrix_row += 2
         cols = 0
         row = 0


Test patch:
diff --git a/sympy/physics/continuum_mechanics/tests/test_truss.py b/sympy/physics/continuum_mechanics/tests/test_truss.py
index addeaa46fb..addd6ce022 100644
--- a/sympy/physics/continuum_mechanics/tests/test_truss.py
+++ b/sympy/physics/continuum_mechanics/tests/test_truss.py
@@ -1,5 +1,6 @@
 from sympy.core.symbol import Symbol, symbols
 from sympy.physics.continuum_mechanics.truss import Truss
+from sympy import sqrt
 
 
 def test_truss():
@@ -99,10 +100,9 @@ def test_truss():
 
     # testing the solve method
     t.solve()
-    assert t.reaction_loads['R_A_x']/P - (-1.4142135623731) < 1e-10
-    assert t.reaction_loads['R_A_y']/P - (-2.41421356237309) < 1e-10
-    assert t.reaction_loads['R_D_y']/P - (-0.5) < 1e-10
-    assert t.internal_forces[AB]/P < 1e-10
-    assert t.internal_forces[CD] < 1e-10
-    assert t.internal_forces[AC] < 1e-10
-    # assert t.internal_forces == {AB: 3.06161699786838e-17*P, CD: 0, AC: 0}
+    assert t.reaction_loads['R_A_x'] == -sqrt(2)*P
+    assert t.reaction_loads['R_A_y'] == -sqrt(2)*P - P
+    assert t.reaction_loads['R_D_y'] == -P/2
+    assert t.internal_forces[AB]/P == 0
+    assert t.internal_forces[CD] == 0
+    assert t.internal_forces[AC] == 0