ss = SystemElements()
ss.add_element([2, 2], g = 1)
ss.add_element([4, 0], g=1)
ss.add_element([6, 2], g=1)
ss.add_element([8, 0], g=1)
ss.add_support_hinged(1)
ss.add_support_hinged(5)
ss.solve()
ss.show_axial_force()

Does not lead to symmetrical loading

I have made some test for nonlinear analysis with anastruct.
First AISC 360 benchmark was checked.
https://notebooks.azure.com/lukasz-laba1/projects/ansastruct/html/AISC%20360%20benchmark-20190301.ipynb
All considered cases agree with AISC benchmarks. It's Great !!
Then I made some test with imperfection modeling (big and small delta).
https://notebooks.azure.com/lukasz-laba1/projects/ansastruct/html/DAM-test_20190228.ipynb
The idea was to create and calculate one by one deformated structures and finally get envelope of results.
I think about some new features in anstruct to make analysis like this more easy:

• method to make curved member with small delta parameter is needed,
• parameter for big delta effect is needed,
• some idea for making envelop is needed.

I thinking about anaStruct contribution. Does this idea for direct nonlinear analysis could be right way to anaStruct upgrade?

Congratulations on your V1.0 Beta release.

Hi, thanks for this awesome project!

Is there a way to change the number of decimal digits after the comma in the plot? Right now it displays only the integer part.

Hi Ritchie,
Congratulations on anaStruct. It's a really great tool.
I have a script that uses the .get_element_results() method and it worked last time I used it a few months ago (sorry, but I don't remember which version of anaStruct I had at the time).
However, now the .get_element_results() method returns

AttributeError: 'Element' object has no attribute 'ai'

This behaviour happens with any model, including the examples in the documentation.
All the best,
Pedro

I have made first test with importing geometry from dxf. EZDXF package was used to grab data from dxf file. I have used this package a lot before.
Here is DxfLoader class draft
https://github.com/lukaszlaba/anaStruct/blob/dxfloader/anastruct/fem/dxfloader/dxfloader.py
It is also possible to make that section be loaded from dxf to.
Is it possible to change section or other section proprieties of element after the element was crated? Some proprieties must be still defined with script code - like hinges for example, so some existing elements must be modified. I think only structure geometry and sections is what we want to import from dxf to be happy.

The following script for a simple king post truss fails with error anastruct.basic.FEMException: ('StabilityError', 'The eigenvalues of the stiffness matrix are non zero, which indicates a instable structure. Check your support conditions'):

Example Script

from anastruct.fem.system import SystemElements

sys = SystemElements(mesh=50)
sys.add_element(
    location=[[0.0, 0.0], [2.5, 0.0]], EA=68300, EI=128, g=0, spring={1: 0, 2: 0}
)
sys.add_element(
    location=[[0.0, 0.0], [2.5, 2.0]], EA=68300, EI=128, g=0, spring={1: 0, 2: 0}
)
sys.add_element(
    location=[[2.5, 0.0], [5.0, 0.0]], EA=68300, EI=128, g=0, spring={1: 0, 2: 0}
)
sys.add_element(
    location=[[2.5, 2.0], [2.5, 0.0]], EA=68300, EI=128, g=0, spring={1: 0, 2: 0}
)
sys.add_element(
    location=[[2.5, 2.0], [5.0, 0.0]], EA=68300, EI=128, g=0, spring={1: 0, 2: 0}
)
sys.add_support_hinged(node_id=1)
sys.add_support_hinged(node_id=4)
sys.point_load(Fx=0, Fy=-20.0, node_id=3)
sys.solve()
sys.show_results()

Discussion

I see that the function ensure_single_hinge() tries to ensure that no node ends up with an unrestrained rotational degree of freedom, by ensuring that all nodes have at least one element with a non-hinged connection to the node. However, that function is broken due to a minor syntax error:
if spring is not None and 0 in spring:
should instead read
if spring is not None and 0 in spring.values():

I've tested a fix for this (along with a couple consequential downstream fixes that come up when this function is fixed). See below diffs to util.py, system.py, node.py, and postprocess.py

Diff of Proposed Fix

diff --git a/anastruct/fem/node.py b/anastruct/fem/node.py
index f905ca0..9dea756 100644
--- a/anastruct/fem/node.py
+++ b/anastruct/fem/node.py
@@ -1,5 +1,5 @@
 class Node:
-    def __init__(self, id, Fx=0, Fz=0, Ty=0, ux=0, uz=0, phi_y=0, vertex=None):
+    def __init__(self, id, Fx=0, Fz=0, Ty=0, ux=0, uz=0, phi_y=0, vertex=None, hinge=False):
         """
         :param id: ID of the node, integer
         :param Fx: Value of Fx
@@ -20,7 +20,7 @@ class Node:
         self.uz = uz
         self.phi_y = phi_y
         self.vertex = vertex
-        self.hinge = False
+        self.hinge = hinge
         self.elements = {}
 
     @property
diff --git a/anastruct/fem/postprocess.py b/anastruct/fem/postprocess.py
index d6b2d93..8750ee3 100644
--- a/anastruct/fem/postprocess.py
+++ b/anastruct/fem/postprocess.py
@@ -20,7 +20,7 @@ class SystemLevel:
 
         for el in self.system.element_map.values():
             # post processor element level
-            self.post_el.node_results(el)
+            self.post_el.node_results(self, el)
 
     def node_results_system(self):
         for k, v in self.system.node_element_map.items():
@@ -88,11 +88,15 @@ class ElementLevel:
         self.system = system
 
     @staticmethod
-    def node_results(element):
+    def node_results(self, element):
         """
         Determine node results on the element level.
 
         """
+        # Check for hinges
+        hinge1 = self.system.node_map[element.node_id1].hinge
+        hinge2 = self.system.node_map[element.node_id2].hinge
+        
         # Global coordinates system
         element.node_map[element.node_id1] = Node(
             id=element.node_id1,
@@ -101,10 +105,15 @@ class ElementLevel:
             Fz=element.element_force_vector[1]
             + element.element_primary_force_vector[1],
             Ty=element.element_force_vector[2]
-            + element.element_primary_force_vector[2],
+            + element.element_primary_force_vector[2]
+            if not hinge1
+            else 0,
             ux=element.element_displacement_vector[0],
             uz=element.element_displacement_vector[1],
-            phi_y=element.element_displacement_vector[2],
+            phi_y=element.element_displacement_vector[2]
+            if not hinge1
+            else 0,
+            hinge=hinge1,
         )
 
         element.node_map[element.node_id2] = Node(
@@ -114,10 +123,15 @@ class ElementLevel:
             Fz=element.element_force_vector[4]
             + element.element_primary_force_vector[4],
             Ty=element.element_force_vector[5]
-            + element.element_primary_force_vector[5],
+            + element.element_primary_force_vector[5]
+            if not hinge2
+            else 0,
             ux=element.element_displacement_vector[3],
             uz=element.element_displacement_vector[4],
-            phi_y=element.element_displacement_vector[5],
+            phi_y=element.element_displacement_vector[5]
+            if not hinge2
+            else 0,
+            hinge=hinge2,
         )
 
         # Local coordinate system. With inclined supports
diff --git a/anastruct/fem/system.py b/anastruct/fem/system.py
index 62d240c..41451c5 100644
--- a/anastruct/fem/system.py
+++ b/anastruct/fem/system.py
@@ -248,7 +248,7 @@ class SystemElements:
         system_components.util.append_node_id(
             self, point_1, point_2, node_id1, node_id2
         )
-        system_components.util.ensure_single_hinge(self, spring, node_id1, node_id2)
+        spring = system_components.util.ensure_single_hinge(self, spring, node_id1, node_id2)
 
         # add element
         element = Element(
@@ -634,7 +634,7 @@ class SystemElements:
 
         for id_ in node_id:
             id_ = _negative_index_to_id(id_, self.node_map.keys())
-            system_components.util.support_check(self, id_)
+            #system_components.util.support_check(self, id_)
 
             # add the support to the support list for the plotter
             self.supports_hinged.append(self.node_map[id_])
@@ -653,7 +653,7 @@ class SystemElements:
 
         for id_ in node_id:
             id_ = _negative_index_to_id(id_, self.node_map.keys())
-            system_components.util.support_check(self, id_)
+            #system_components.util.support_check(self, id_)
 
             if direction == "x":
                 direction = 2
diff --git a/anastruct/fem/system_components/util.py b/anastruct/fem/system_components/util.py
index e757a7e..8e2e212 100644
--- a/anastruct/fem/system_components/util.py
+++ b/anastruct/fem/system_components/util.py
@@ -4,31 +4,68 @@ from anastruct.basic import FEMException, angle_x_axis
 
 
 def ensure_single_hinge(system, spring, node_id1, node_id2):
-    if spring is not None and 0 in spring:
+    if spring is not None and 0 in spring.values():
         """
         Must be one rotational fixed element per node. Thus keep track of the hinges (k == 0).
         """
 
         for node in range(1, 3):
-            if spring[node] == 0:  # node is a hinged node
-                if node == 1:
-                    node_id = node_id1
-                else:
-                    node_id = node_id2
-
-                system.node_map[node_id].hinge = True
-
+            if node == 1:
+                node_id = node_id1
+            else:
+                node_id = node_id2
+            
+            if node in spring.keys() and spring[node] == 0:  # node is a hinged node
                 if len(system.node_map[node_id].elements) > 0:
-                    pass_hinge = not all(
-                        [
-                            el.hinge == node
-                            for el in system.node_map[node_id].elements.values()
-                        ]
-                    )
+                    hinges = []
+                    for el in system.node_map[node_id].elements.values():
+                        if node_id == el.node_id1:
+                            hinges.append(1 in el.springs.keys() and el.springs[1] == 0)
+                        elif node_id == el.node_id2:
+                            hinges.append(2 in el.springs.keys() and el.springs[2] == 0)
+                    pass_hinge = not all(hinges)
                 else:
                     pass_hinge = True
                 if not pass_hinge:
+                    system.node_map[node_id].hinge = True
                     del spring[node]  # too many hinges at that element.
+            
+            elif node not in spring.keys() or (node in spring.keys() and spring[node] != 0):
+                """
+                If a fixed element is added after a hinged element,
+                then add the removed spring release back in and set node as not hinged
+                """
+                system.node_map[node_id].hinge = False
+                if len(system.node_map[node_id].elements) > 0:
+                    for el in system.node_map[node_id].elements.values():
+                        if node_id == el.node_id1:
+                            system.element_map[el.id].springs.update({
+                                1: 0
+                            })
+                        elif node_id == el.node_id2:
+                            system.element_map[el.id].springs.update({
+                                2: 0
+                            })
+    else:
+        """
+        If a fixed element is added after a hinged element,
+        then add the removed spring release back in and set node as not hinged
+        """
+        if system.node_map[node_id1].hinge:
+            system.node_map[node_id1].hinge = False
+            for el in system.node_map[node_id1].elements.values():
+                system.element_map[el.id].springs.update({
+                    1: 0
+                })
+
+        if system.node_map[node_id2].hinge:
+            system.node_map[node_id2].hinge = False
+            for el in system.node_map[node_id2].elements.values():
+                system.element_map[el.id].springs.update({
+                    2: 0
+                })
+    
+    return spring
 
 
 def append_node_id(self, point_1, point_2, node_id1, node_id2):
@@ -75,7 +112,7 @@ def det_node_ids(system, point_1, point_2):
 def support_check(system, node_id):
     if system.node_map[node_id].hinge:
         raise FEMException(
-            "Flawed inputs", "You cannot add a support to a hinged node."
+            "Flawed inputs", "You cannot add a rotation-restraining support to a hinged node."
         )

A

Hello there.

I came across anaStruct and I couldn't help myself to try it out. One thing that's struck me is that when I tried running a simple cantilever simulation, the deflection given in the charts (ss.show_displacement()) is not the maximum displacement? It seems to give me the displacement at the midpoint between the nodes?

After some digging and troubleshooting I found that the nodal displacement information gives me results equal to what I was expecting. Am I doing something wrong?

``E = 210 * 1e9 # GPa
F= -50 # kN
Ix = 360000 # mm^4
A = 1200 # mm^2

Expected results:

print((F1e3 * 1**3)/(3Esection.Ix1e-12))

anaStruct:

ss = SystemElements(figsize=(15,15))
ss.add_element(location=[[0,0],[1,0]],EI=(EIx1e-12 * 1e-3),EA=(EA1e-9))
ss.add_support_fixed(node_id=1)
ss.point_load(node_id=2,Fy=F)

ss.solve()
ss.show_displacement()
print(ss.get_node_displacements(node_id=2))
print(ss.get_node_result_range('uy'))
print(ss.get_element_results(element_id=1))``

Anastruct output:
{'id': 2, 'ux': 0.0, 'uy': -0.22045855379188709, 'phi_y': 0.3306878306878306}
[-0.0, -0.22045855379188709]
{'id': 1, 'length': 1.0, 'alpha': -0.0, 'u': 0.0, 'N': 0.0, 'wmax': -0.03995224242642931, 'wmin': 0.002646510934973756, 'w': None, 'Mmin': 0.0, 'Mmax': 49.99999999999999, 'M': None, 'q': 0}

I've got a snippet of some demo code + outputs below. You can see that "get_element_results" and "get_element_result_range" result in different outputs for Max Moment. It seems that the get_element_result_range function is actually pulling in node results.

`# Solve and display results
ss.solve()
ss.show_results()

moment_results = ss.get_element_result_range("moment")
shear_results = ss.get_element_result_range("shear")
axial_results = ss.get_element_result_range("axial")
displacement_results = ss.get_node_result_range("uy")


element_result_string = str(ss.get_element_results(0))
for results in ss.get_element_results(0):
    for data in results.items():
        print(data)
    print("\n")

print(moment_results)
print(shear_results)

`

Can I just have the values appear on the show_axial_force()?

Also can I combine show_structure with show_reaction_force and axial_force?
So that I can see internal forces and external forces without the need of shaded diagrams which obstruct everything.

D

Is it possible to have a load case combination (e.g. 1.25 Dead, 1.5 Live, 0.5 Snow), where the 'dead load' case is the only load case where self weight is considered?

From my understanding, self weight appears in all load cases.

Hi,
Is there a way to define a new member using the node numbers already defined?
eg
ss.add_element(location=[ss.node[1], ss.node[2])

Only work if the direction is similar (y). Make independent.

How can I make a support that only resists vertical loads?
I get that fixed can do vertical and horizontal, but why isn't there a vertical only option?

90 degr > angle < 270 degr result into errors and flawed results

It would be nice to have more section defined in database and have section name parameter for method that create element. Also section mark on show_structure figure would be great. Please find strupy project where already big database of section exist
https://bitbucket.org/struthonteam/strupy/src/master/
There are US, UK and EU sections
https://bitbucket.org/struthonteam/strupy/src/master/strupy/steel/database_sections/
Strupy package can be used as new dependence, or similar solution can be implement in anastruct what is better idea I think. Dose it is good way for anaStruct upgrade?

NumPy's latest release has implemented a new function for assembling arrays. It reminded me of how stiffness matrixes are built and (though I have not reviewed your code) thought you might find it useful in a refactoring pass.

Problem Description

This displacement plot of a uniformly-loaded propped beam has a sharp corner at the second support that's obviously not correct:

Code for example above:

from anastruct.fem.system import SystemElements
sys = SystemElements(mesh=250)
sys.add_element(location=[(0.0, 0), (1.3333333333333333, 0)],
                EA=356000.0,
                EI=1332.0000000000002)
sys.add_element(location=[(1.3333333333333333, 0), (2.0, 0)],
                EA=356000.0,
                EI=1332.0000000000002)
sys.add_support_hinged(node_id=sys.find_node_id((0, 0)))
sys.add_support_hinged(node_id=sys.find_node_id((1.3333333333333333, 0)))
sys.q_load(q=-1.0, element_id=1, direction="y")
sys.q_load(q=-1.0, element_id=2, direction="y")
sys.solve()
sys.show_displacement()

Discussion

I saw that a few months ago, the elemental deflections were changed to calculate by means of averaging integrations from both the left and the right. However, there seems to be a flaw in the logic of this, which results in some odd deflection graph shapes and slightly incorrect deflections. The offending code appears to be the following in postprocess.py/determine_displacements():

# Take the average of cumulative integration from both sides. This is due to numerical differences, that
            # would be observable in symmetrical structures.
            phi_neg = (
                -0.5
                * (
                    integrate_array(element.bending_moment, dx)
                    + integrate_array(element.bending_moment[::-1], dx)
                )
                / element.EI
            )
            w = integrate_array(phi_neg, dx)

The first integrate_array() returns a cumulative sum of the moment from, say, left to right. The second integrate_array(), however, flips the moment array such that the moments themselves are right to left, and it never flips them back. In effect, therefore, the deflections from the left side of the beam are partially averaged with the deflection from the right side of the beam.

Since maximum deflections in a beam are typically near to the beginning, middle, or end of the beam, the maximum deflections are essentially unaffected by this problem - but deflections near the quarter points of asymmetric beams in particular are notably off.

Proposed solution

I like the idea of integrating in both directions, and that definitely does clean up some asymmetries that shouldn't be there. However, I would suggest running through the full integration in each direction entirely independently, and then flip the reversed array back before averaging both sets of integrations back together:

            dx = element.l / (len(element.bending_moment) - 1)
            lx = np.linspace(0, element.l, con)

            # Take the average of cumulative integration from both sides. This is due to numerical differences, that
            # would be observable in symmetrical structures.
            phi_neg1 = -integrate_array(element.bending_moment,
                                        dx) / element.EI
            w1 = integrate_array(phi_neg1, dx)

            # Angle between last w and elements axis. The w array will be corrected so that this angle == 0.
            alpha1 = np.arctan(w1[-1] / element.l)
            w1 = w1 - lx * np.tan(alpha1)

            phi_neg2 = -integrate_array(element.bending_moment[::-1],
                                        dx) / element.EI
            w2 = integrate_array(phi_neg2, dx)

            # Angle between last w and elements axis. The w array will be corrected so that this angle == 0.
            alpha2 = np.arctan(w2[-1] / element.l)
            w2 = w2[::-1] - lx[::-1] * np.tan(alpha2)

            element.deflection = -(w1 + w2) / 2
            element.max_deflection = np.max(np.abs(element.deflection))

Output of example script with proposed solution

Observe that the slope at that second support is now nearly flat, and there's now a smooth transition at that support.

Verification

I put together a quick Excel spreadsheet for checking the deflections of a propped beam, based upon the analytical equations for its moments, shears, and deflections. The results are significantly off in the current anaStruct version, but from spot-checking, seem to be accurate to about the third decimal place with the above proposed solution.
Propped Beam Check Sheet.xlsx

Hi,
How to clear information from memory?
I want to model one structure after another, I'm using anaStruct in a loop to check if a certain bar meets an established parameter. However, an error occurs due to an overlap of data from the second structure with the first.

I was wondering if it was possible to change the member properties (EA, EI) after the model has been defined and solved? I would like to update the members based on displacement.
Currently I am redefining the structure each time.

The following two routines should give identical results (after transforming axes), but they do not. With inclined roller supports and q_loads, incorrect reactions are being calculated, and sum of loads does not equal sum of reactions. Point loads appear to work correctly; the problem appears to only be with q_loads.

I tried tracing this back to find the problem. It looks like it's something to do with the element_primary_force_vector not being transformed correctly with inclined roller supports, but I wasn't able to pin down exactly where the problem is.

Inclined Beam with Y-Axis Roller

from anastruct import SystemElements
ss = SystemElements(EA=356000, EI=1330)
ss.add_element(location=[[0,0],[7.07,7.07]])
ss.add_support_hinged(node_id=1)
ss.add_support_roll(node_id=2, direction="y")
ss.q_load(q=-1, element_id=1, direction="element")
ss.solve()

Note that sum of reactions perpendicular to the beam are correctly equal to the total applied load of 10, and an axial load is present.

Inclined Roller Support

from anastruct import SystemElements
ss = SystemElements(EA=356000, EI=1330)
ss.add_element(location=[[0,0],[10,0]])
ss.add_support_hinged(node_id=1)
ss.add_support_roll(node_id=2, angle=45)
ss.q_load(q=-1, element_id=1, direction="element")
ss.solve()

Note that sum of reactions perpendicular to the beam are incorrectly equal to 8.54, and no axial load is present.

If an axial load is applied at the same node at which an inclined roller support is located, then the results are blatantly wrong, with sum of loads != sum of reactions.

Example script:

from anastruct.fem.system import SystemElements
sys = SystemElements(mesh=250)
sys.add_element(location=[(0.0, 0), (1.0, 0)], EA=356000.0, EI=1332.0000000000002)
sys.add_support_hinged(node_id=sys.find_node_id((0, 0)))
sys.add_support_roll(node_id=sys.find_node_id((1.0, 0)), angle=45)
sys.point_load(Fx=10.0, Fy=0, node_id=2)
sys.solve()
sys.validate()
sys.show_reaction_force()

Which results in the following:

Note that sum of reactions = 30, though total applied load = 10.

Hi @ritchie46 & @lukaszlaba & @rodrigo1392 ,

I worked on two days about this issue. Can you please help me. This is very important. How to fix this Value error. When we using reaction force this value error came. axial, bending, deflection forces working. Can you please help me. show_reaction_force(show=True) is working when we did not used load combinations.

from anastruct.fem.system import SystemElements
from anastruct.fem.util.load import LoadCase, LoadCombination
import matplotlib.pyplot as plt
ss = SystemElements(EA=14000, EI=5.36)
ss.add_element(location=[[0, 0], [6, 0]])
ss.discretize(n=6)
ss.add_support_fixed(node_id=1)
ss.add_support_fixed(node_id=7)
lc_dead = LoadCase('dead')
lc_dead.q_load(q=-10, element_id=2)
lc_dead.point_load(Fx=30, Fy=-50, node_id=2)
combination1 = LoadCombination('SLS')
combination1.add_load_case(lc_dead, factor=1.0)
results1 = combination1.solve(ss)
for k, ss in results1.items():
    if k == 'combination':
         results1[k].show_reaction_force(show=True)
         results1[k].show_bending_moment(show=True)
         results1[k].show_displacement(show=True)
         results1[k].show_axial_force(show=True)

File "E:/main/fem_analysis/l_profile.py", line 72, in l_profile_loads
    results1[k].show_reaction_force(show=True)
  File "E:\models\fem_model\fem\system.py", line 910, in show_reaction_force
    return self.plotter.reaction_force(figsize, verbosity, scale, offset, show)
  File "E:\models\fem_model\fem\plotter\mpl.py", line 669, in reaction_force
    max_force = max(
ValueError: max() arg is an empty sequence

While the following code is an extreme example, it gives an obviously nonsense displacement plot. Note that the maximum displacement (at the midpoint) is correct, but there's an incorrect discontinuity in the plot.

from anastruct import SystemElements
ss = SystemElements(EA=356000, EI=1330)
ss.add_element(location=[[0,0],[10,0]])
ss.add_element(location=[[10,0],[20,0]])
ss.add_support_hinged(node_id=1)
ss.add_support_hinged(node_id=3)
ss.point_load(Fy=-500, node_id=2)
ss.solve()
ss.show_displacement()

I believe the problem is on line 221 of fem/postprocess.py, where
w = w - lx * np.sin(alpha)
should actually be
w = w - lx * np.tan(alpha)

Making this change results in the following displacement graph instead (and the same maximum displacement):

Hi,
See the following model:

from anastruct.fem.system import SystemElements
ss = SystemElements()
ss.add_element(location=[[0, 0], [5, 0]], EA=5e9, EI=8000, spring={2:0})
ss.add_element(location=[[5, 0], [10, 0]], EA=5e9, EI=8000)

ss.q_load(element_id=1, q=-9)
ss.q_load(element_id=2, q=-9)

ss.add_support_fixed(node_id=1)
ss.add_support_fixed(node_id=3)

ss.solve()
ss.show_structure()
ss.show_reaction_force()
ss.show_shear_force()
ss.show_bending_moment()
ss.show_displacement()

It should give at least symmetrical results with support loads of 45 and bending moments of 112.5. Adding a spring with a numerical small value will also lead to different results.

Small request; the use of a kinematic matrix and a constitutive matrix is interesting, is there anyone who has some literature about that?

Marcel

Hi,
I'm using anaStruct in my thesis, I would like to cite your project. How do I should cite your project?
Do you have the BibTex format?

Thanks.

Hi Ritchie,

I just installed anaStruct, and so far I like it a lot, so thanks for sharing! However, when playing around with a simple cantilever beam case, I'm running into a number of errors that I don't quite understand, so I was hoping you could give some explanation.

So here's my toy problem:

from anastruct import SystemElements

Ne = 10     # number of elements
F_y = 10    # point load

ss = SystemElements(EI=5e3, EA=1e5)
ss.add_multiple_elements([[0, 0], [10, 0]], Ne)
# clamp first node
ss.add_support_fixed(1)
# point load on last node
ss.point_load(ss.id_last_node, Fy = -F_y)

# solve, linear
ss.solve()

# solve with geometrical non-linearity
ss.solve(geometrical_non_linear = True)

With the values of Ne and F_y used above, everything goes fine with the linear solver (I'll get to the geometrical_non_linear later on). However, if I increase Ne to, say, 100, ss.solve() gives me this error:

Traceback (most recent call last):

  File "<ipython-input-16-a7e1225dab88>", line 1, in <module>
    ss.solve()

  File "/Applications/miniconda3/lib/python3.6/site-packages/anastruct/fem/system.py", line 532, in solve
    "The eigenvalues of the stiffness matrix are non zero, "

FEMException: ('StabilityError', 'The eigenvalues of the stiffness matrix are non zero, which indicates a instable structure. Check your support conditions')

Solving the same problem (i.e. Ne = 10, Fy = 10) with geometrical_non_linear set to True results in the following error (independent of whether I add extra discretisation, e.g. discretize_kwargs=dict(n=20)):

Traceback (most recent call last):

  File "<ipython-input-47-3c8b776ea164>", line 1, in <module>
    ss.solve(geometrical_non_linear = True)

  File "/Applications/miniconda3/lib/python3.6/site-packages/anastruct/fem/system.py", line 554, in solve
    self, verbosity, discretize_kwargs=discretize_kwargs

  File "/Applications/miniconda3/lib/python3.6/site-packages/anastruct/fem/system_components/solver.py", line 132, in geometrically_non_linear
    buckling_factor = det_linear_buckling(buckling_system)

  File "/Applications/miniconda3/lib/python3.6/site-packages/anastruct/fem/system_components/solver.py", line 104, in det_linear_buckling
    kg = system.reduced_system_matrix - k0

ValueError: operands could not be broadcast together with shapes (20,20) (30,30) 

Can you explain what's causing these errors?

Many thanks,
Sita

P.S. I just noticed my code might be a bit unclear/misleading, so to clarify: I'm using either ss.solve() or ss.solve(geometrical_non_linear = True), not both in one go.

@ritchie46 what breaks with Python 3.5?

Is it possible to add a roll support not in x or y direction, but, for example in 45° of x-axis?

results = combination.solve(ss)
for k, ss in results.items():
results1[k].show_reaction_force(show=True)

When i used load comination it show error.

File "E:/main/fem_analysis/t_profile.py", line 66, in t_profile_loads
#results1[k].show_displacement(show=True)

File "E:\fem_model\fem\system.py", line 910, in show_reaction_force
return self.plotter.reaction_force(figsize, verbosity, scale, offset, show)
File "E:\fem_model\fem\plotter\mpl.py", line 669, in reaction_force
max_force = max(
ValueError: max() arg is an empty sequence

Hi @ritchie46!

Thanks for creating this library! So far it has been really interesting to play around with. I'm looking to use it to test my current bridge design before building it out. However, I'm pretty new to structural analysis + load simulation in Python, and would love any insights you have! The goal would be to visualize points of tension and the load distribution in our bridge, alongside the maximum load (in kilograms the bridge can support.

Here are a couple of diagrams of our design so far! Any help is gladly appreciated!

Hi! I have a university project where I must calculate the forces and strains on a fixed beam using the Python language and would like to represent them in 2D. I created the code below, but it is not running on Spyder. Can you help me?

a=float(input("Enter a value for a (mm):"))
P=float(input("Enter a value for P (KN)"))
from anastruct import SystemElements
ss = SystemElements()
ss.add_element(location=[700,0])
ss.add_support_fixed(node_id=1)
ss.point_load(ss.id=[a], Fy=-1*P)
ss.show_structure()

Hi @ritchie46 ,

If possible, can u please upload sample axial and reaction forces with load combination examples

HI @ritchie46 ,@lukaszlaba, @rodrigo1392
I want display results(values or json or object info) with out plot(graphs) when we use Show_results().

Thanks for advance.

The system matrix is changed when a spring is added. In the non linear solver this matrix is removed.

42 for el in ss.element_map.values():
43 # apply wind load on elements that are vertical
---> 44 if np.isclose(np.sin(el.ai), 1):
45 ss.q_load(
46 q=1,

AttributeError: 'Element' object has no attribute 'ai'

I am getting attribute error issue intermediate program. How to solve this