Hi. First off, PyNite is a great project and i'm very excited to see the pace at which it is improving.

Describe the bug Quad elements don't seem to produce a valid result if they're skewed from rectangular, but no interpreter errors are raised

To Reproduce The below code produces a rectangular test grid in order to test out the quad functionality. It seems to work fine when i set the skew to zero (ie. all quads are rectangular), but as soon as the skew is non-zero all of the reactions and displacements become 'nan' and the displacement plot no longer correctly plots the quad displacements.

Apologies as the code runs on my own fork, in which i've added element auto-naming (when 'None' is passed to the constructor function)...

# Import 'FEModel3D' and 'Visualization' from 'PyNite'
from PyNite import FEModel3D
from PyNite import Visualization

# Create a new model
plate_test = FEModel3D()

hplates=10
vplates=20
plate_dim=0.5
skew=0.01

# Create a test model with a grid of quad elements
node_array=[]
plate_array=[]
support_nodes=[]
for yi in range(vplates+1):
    y2=(yi+1)*plate_dim 
    
    current_row=[]
    current_row_plates=[]
    for xi in range(hplates+1):
        x2=(xi+1)*plate_dim * (1+yi*skew)
        
        current_row.append(plate_test.AddNode(None, x2, y2, 0))
        
        if (xi==2 or (xi>0 and xi % 8==0)) and (yi==2 or (yi>0 and yi % 8 ==0)):
            plate_test.DefineSupport(current_row[-1], True, True, True, False, False, False)
            support_nodes.append(current_row[-1])
        
        if yi>0 and xi>0:
            n1 = node_array[yi-1][xi-1]
            n2 = current_row[xi-1]
            n3 = current_row[xi]
            n4 = node_array[yi-1][xi]
            
            nq=plate_test.AddQuad(None,n1,n2,n3,n4,0.01,2e11,0.25)

            plate_test.AddQuadSurfacePressure(nq,1000)
            current_row_plates.append(nq)
            
                                      
    node_array.append(current_row)
    plate_array.append(current_row_plates)
        

# # Analyze the model
plate_test.Analyze()

#Print critical node reactions
for i,n in enumerate(support_nodes):
    print("Node",i,"reaction is",plate_test.GetNode(n).RxnFZ["Combo 1"])
  
#Print node displacements
i=0
for row in node_array:
    for nd in row:
        print("Node",i,"displacement is",plate_test.GetNode(nd).DZ["Combo 1"])
        i+=1

# # Render the model for viewing.
Visualization.RenderModel(plate_test, text_height=0.05, deformed_scale=0.001, deformed_shape=False, color_map='dz', render_loads=True)

In order to obtaine 3D deformed shape, I've introduced axial deformation and axial translation in PyNite. I've defined two new functions in BeamSegZ.py:

"Translation1" is defined in Member3D.py forsì SegmentsZ[0] and SegmentsZ[i].

With these new functions, deformed shape of 2D frame shown in pull request Update Visualization.py #30 now it's correct.

There are some problems with axial distributed load (p1 and p2); I think because I have to add other terms (terms that includes p1 and p2) in Axial translation's formula. In this moment, using this formulation, with a simple supported beam with axial distributed load, I otained the follow axial translations:

I've tried to understand formula used in Deflection function in BeamSegZ but I don't remeber how get terms who account for distributed load p1 and p2.

Thaks for supporting!

I'd love to see functionality where the solver automatically recognizes the presence of an intermediate node within the length of a member, and internally divides the member into smaller sub-members. With that functionality, an intermediate framing member will behave correctly, without having to subdivide the member by hand. This functionality is really nice to have in commercial solvers.

Hi,

I just installed PyNite with pip.

Installation ended correctly.

But when I want to use it, if I type "import pynite", I get "No module named 'PyNite'".

Other imports such as "import numpy" or whatever work just fine.

I work with WinPython-64-bit-3.6.0.1.Qt5.

I'm still trying to get my head around the maths of it, but:

1. A mechanism for orientating beams by specifying a reference point in the beam's ZX-plane.
2. Don't create the transformation matrix, T(), more than once.
3. I think T() needs to be transposed when transforming back to the global coordinate frame.

Hi,

I've been putting PyNite through its paces recently and it seems to be holding up really well, so congrats on creating such a robust FEA solver. I have a few ideas for improvements i'd like to raise for discussion:

1. The Scipy sparse solver doesnt appear to throw an exception when the input stiffness matrix is singular, it simply prints a warning. This isn't ideal for the operation of PyNite, where we want to be able to catch this condition to resolve the instability. I'm not too familiar with Scipy so not sure if there's a way to force it to throw an exception with a singular input matrix, but as a work around maybe PyNite could check the first element in the returned displacement Matrix and if it is 'NaN' then raise the exception manually
2. The solver is very slow for large models with multiple load combinations. The reason for this seems to be that the stiffness matrix is regenerated (and more critically, re-inversed) for each load combination. For linear models (no tension/compression only elements) this seems entirely unnecessry, the stiffness matrix should be able to be generated once, then re-used for every load combination. As an even further optimisation, the stiffness matrix could be pre-inversed then just multiplied with the force matrix for each load combo
3. The above point is not possible for non-linear models however as a thought, potentially the full stiffness matrix could be generated once then just have the stiffness components for the inactive compression/tension only elements subtracted in each iteration/load combo. Unsure what sort of performance increase this would provide as you'd still have to do the matrix inversion in each iteration.

Would be interested to hear your thoughts on the above.

Happy to contribute some code, but I am not particularly familiar with the implementation specifics of the matrix operations I've discussed.

Description: The calculated values for the moment and deflection do not match the expected values printed at the end.

To Reproduce: Steps to reproduce the behavior:

1. Run Testing & Debuging/AISC 2nd Order Benchmark Check.py
2. See value mismatch in the console output.

Output Displayed:

PyNite Calculated Moment:  -99.99499999998429
PyNite Calculated Displacement:  7.944827586205967
AISC Benchmark Problem Moment:  435.6165910296658
AISC Benchmark Problem Displacement:  40.27399092355989

I introduced the possibility to apply support's settlements and nodal displacements/rotations at a free node. Because of this new feature, solving method is changed because partioning of matrix K and vectors P and FER is needed. I hope you like my work and you take it into account for PyNiteFEA develope.

Thanks!

I would like to introduce settlement of a support or, in general, the possibility to choose a nodal displacement. Pratically it rappresents a new load for the model, but it rappresent also a non homogeneus boundary condition.

Could you give me some suggestions (guidelines)?

Dear Craig and hello everyone,

as described in FEModel3D, scipy is used to make pynite work faster with a huge amount of nodes. It would be great, if there is a possibility to use pynite without scipy also. Pynite allows to pass sparse=False as argument to the function analyze in FEModel3D:

truss.analyze(check_statics=False, sparse=False)

Anyway: Scipy is still imported in FEModel3D.py Furthermore coo_matrix and lil_matrix is is called in the code afterwards. Would it be possible to make conditionals in order to avoid to load and to use scipy? (I would like to make a simple add-on for pynite in blender 3D. When running blender 3D on linux, it is possible to use pip in order to install additional modules to blender python. When switching to other operating systems, this would be a limitation. But maybe someone has also another idea how to make this work.)

Thanks for your answers in advance!

Greetings, Chris

Quads still seem to have some issues when they are not aligned with the global coordinate system. Currently working on a fix to the element formulation.

Hi Craig,

I can see the code for assembling and solving FE are all written in FEModel3D.py. It might be beneficial to have the solution part of the code written in a separate script to maintain the code tidy if more functionalities are added (different nonlinear solvers, dynamic solver, etc) in the future.

BR Roger

Sorry to come again with problems and not solutions, but i'm having some trouble with a proof-of-concept lateral stability model that i have been trying to analyse with PyNite. It is a model of an RC building core and shear wall, with some rigid members connecting the two to represent the rigid floor diaphragm.

The model runs fine, and the stress plots on the quad elements look good, but the model fails statics checks. The global reactions are larger than the sum of the applied forces. The 'check_statics' function built in to the 'analyze' method confirms this.

Any thoughts on what could be causing this?

The code to generate the model is complex so i've pickled the model and attached it in a zip file below:

stability_test.zip

I have subclassed PyNite.FEModel3D so you'll need to use the below code to read the model.

import pickle

infile = ...

class customUnpickler(pickle.Unpickler):
    def find_class(self, module, name):
        if module == "Analysis.AnalysisBase":
            module = "PyNite.FEModel3D"
            name = "FEModel3D"

        return super().find_class(module, name)
    
with open(infile,mode="rb") as f: model=customUnpickler(f).load()

Also, i've noticed that Visualization._PrepContour fails when there are nodes in the model not connected to Quad or Plate elements. node.contour = sum(node.contour)/len(node.contour) results in a divide by zero error.

I know it's a big one, but code checks would be useful for various material types: steel, concrete, timber, and masonry. Building a whole Excel spreadsheet to do it, then manually import the results into it, is complex and time-consuming.

First I just want to say I am very impressed with this project! I too have been writing/developing an FEA package currently as well. My primary use case is for doing tilt up/slender wall panel analysis design. After finding your project I am considering now just forking/contributing here. The main things that I see that would still need to be implemented are:

1. Support for self-weight of the finite elements. I could have missed this but I don't see any density functions etc.

2. Extend the meshing routines to support quad/rect FE generations with voids/openings

3. Support for Cracked sections/moment of inertias. Most of the design codes require assuming a fully cracked section. Therefore I think I would need to add a few new variables to allow some sort of coefficients to reduce the gross moment of inertia based purely off the section geometry.

4. Composite Sections? The cracked moment of inertia is a function of reinforcement and axial load at the section. So I suppose I may also want to add support for maybe a "Reinforced Concrete" material and add some logic to modify the stiffness matrix based on this.

5. Easy way to add line loads + eccentricity to the meshed structure and convert these to lumped nodal loads/moments at the individual FE plates.

6. Second order analysis for geometric stiffness to handle the P-delta contribution to ultimate load moments.

I am pretty sure I know how I would code most of these as for implementation but I am curious of your inputs if you have already implemented any of the above and I missed it. Or if you have any thoughts or see any fundamental issues to implementing these.

Again great work and I wish I had found your library sooner!

Hi Graig,

Hope you are doing well and hope you remember me. I have released the standardsolver beam analysis and it is getting some good reviews. It is not very popular yet but I am spreading the message. Got some good reviews in Reddit.

https://www.reddit.com/r/StructuralEngineering/comments/kgyzcr/free_beam_analysis_software/

I am now working on building a similar calculation for 2D frame analysis.

I have a feature request if it is possible.

For nodal supports, they can be only aligned in global X, Y. Is there a way to specify an angle to them with respect to the global axis, For example, if you want to apply only lateral restrain to an inclined member and allow the member to move axially free, At present, it seems it is not possible from reading the wiki.

Regards Lino

It will be usefull if Plate3D element can support non-uniform (trapezodial) loads. Is there any way to load Plate3D element non-uniformly that i can't see?