Skewed Quads Not Producing Results

[bjhowie]

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)

[JWock82]

Thanks for reporting this. I’ll make fixing this my next priority.

[JWock82]

I'm not seeing the same issue. I skewed the grid on the rectangular wall example 10% and the results look great (see below). The deflections, moments, and shears are all close to the values for a non-skewed wall, which is to be expected for a small skew.

Has your fork pulled in the latest changes to the Quad element? Version 0.0.24 is the latest version. Have you perhaps changed something in your fork? Is there something else going on with your code? Perhaps with the way you are assembling your mesh or applying your supports?

.

[bjhowie]

Interesting. Yes my fork has pulled all of the latest changes.

I've modified the code to run on your version 0.0.24 and am still getting the same issue. It works perfectly without any skew (which i've validated with a Strand7 model), but gives me nothing with the skew. If I change the order of the quad node definition (to be anti-clockwise when viewed from the default view position in vtk (+Z?) ) it will give results (instead of just 'nan') but they look to be quite incorrect.

Are you able to try running the below code to see if you can replicate the issue?

# 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
distort=0

# Create a test model with a grid of quad elements
node_array=[]
plate_array=[]
support_nodes=[]
nc,qc = 0,0
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*distort) + skew*yi

        nn = str(nc)
        plate_test.AddNode(nn, x2, y2, 0)
        current_row.append(nn)
        nc+=1

        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]

            print("Plate Definition:",n4,n3,n2,n1)
            nq=str(qc)
            plate_test.AddQuad(nq,n1,n2,n3,n4,0.01,2e11,0.25)
            #plate_test.AddQuad(nq,n4,n3,n2,n1,0.01,2e11,0.25)
            qc+=1

            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)

[JWock82]

I’m fairly certain PyNite is handling this correctly. I’m having a hard time following your code’s logic without any comments in it. There are some red flags in your code that jump out at me. Why are you using the mod operator % with a value of 8? There will be very few instances where that is zero, and I don’t think it’s at the boundaries of your mesh. Also, comparing calculated values to exact zero ignores computer numerical precision limits and is generally a bad practice. I suggest using isclose() instead, or at least rounding the values to some reasonable decimal place. Why is the one dimensional array node_array being sliced twice using the [ ][ ] notation? Is your first slicer supposed to include a range instead of an index? Once you slice it once there’s only one value in it.

I'll run your code this weekend and see what I get. My suspicion is that it's not placing nodes and supports where you think it is.

[JWock82]

The problem was on my side, not yours. I think I solved it. It looks like it was only working when the local x-axis was parallel to the global X-axis. It should be resolved now. Let me know if it's working as expected or not. Thanks for finding this issue. I have not pushed it to PyPI just yet, so you'll need to use the github version to verify for now. I'll get around to that later today.

[JWock82]

Also it shouldn’t matter whether you define the nodes clockwise or counter clockwise now. It should just flip the plate’s local axes when you do that.

[bjhowie]

Apologies, I'd have commented the code better but i was only using it for testing and was pretty confident it was correct. The intent is that the test is representative of a slab in bending, so i'm using the mod 8 operator to add a point support at roughly every 8 plates, the actual value doesn't matter that much. Unsure to where you're referring with the float comparison to zero, the only place i'm checking for equality with zero is after the mod operator which returns an integer anyway. node_arrayand plate_arrayare both 2-dimensional.

I've pulled your updates and am now getting deflection results, instead of just 'nan', but they dont look correct.

This is the deflection plot with a skew of 0. It matches what i get in Strand7 pretty well.

This is the deflection plot with a skew of 0.01.

With such a small skew, the results should be very similar, but as you can see they are not.

[JWock82]

Strange behavior for sure. I’ll keep looking into it.

[JWock82]

It’s got to have something to do with the local x-axis not being parallel to the global one. When they are parallel it seems to work fine for skewed plates.

[JWock82]

I've narrowed the issue down to the MITC4 B matrix. I'll have to go back through the textbook I used to derive it. It may be that the closed form equations in the book were developed for the nodes to be defined in a specific order. If I reorder how your nodes are defined to be plate_test.AddQuad(nq,n3,n2,n1,n4,0.01,2e11,0.25) I get the correct displacement plot. I'll keep digging into this.

[bjhowie]

Glad you're making progress. It's far beyond my understanding unfortunately. I never managed to get past frame elements in 2D space.

Good luck!

[JWock82]

Update: I haven't forgot about this. Just working through it slowly. I ran my quad example problem yesterday with a skew applied to it, and the nodes defined in a counter-clockwise fashion (similar to your model), and I got good results, which is really confusing me. It suggests to me that the order the nodes are defined in and the local axes orientation is not the issue. Yet the order of the nodes causes issues in your file. I'm going to have a closer look at your file to understand it better because I haven't been able to identify any obvious issues in PyNite.

[JWock82]

I went over your file very closely today and I don't see any problems with it. I think this problem is going to take a considerable effort to root out. Probably a lot of math. I'll continue working on it.

[bjhowie]

That's good to hear i guess.

What happens if you rotate your quad example problem so that the quad local axes don't align with the global axes? That seems to be a contributing factor to the issue.

[JWock82]

I've reviewed my derivation of this element closely multiple times, and I'm certain I've followed the procedures and equations in the reference text exactly. This problem keeps evading me and is slowing down production of other features. I'm going to sideline it for now. The effort required to find the bug seems to be more than the feature is worth when there are so many other important features I'd like to add and refine. I'll tag this issue as "help-wanted" for now. Perhaps someone with a little more knowledge of these elements will be able to spot the error in the meantime.

[bjhowie]

An error in the reference material perhaps? Must have been very frustrating, hopefully somebody can figure it out eventually.

I'm sure you're sick of looking at matrix derivations, but i'd love to see triangular elements added in the future. These will work much easier with most meshing algorithms for use on wall and slab type structures.

[JWock82]

Perhaps. I found an older version of the equations online, and it had negative signs in different places in the equations and used a different sign convention for plate section curvatures. I applied those equations and got different answers, but still erroneous, so I went back to the newer equations.

I've placed my derivation in a Jupyter Notebook in the "Derivations" folder for anyone who wants to look over my solution more closely: MITC 4 Quad Element.ipynb. You'll need to install Jupyter or Jupyter Lab to open and view it properly. It doesn't display the equations correctly in GitHub's viewer.

[bjhowie]

So some progress; i decided i'd take a look at this problem again and the it appears to have at least partially resolved itself. When I run the code above with the latest PyNite version and a skew applied, the results look good..... I even applied a random jiggle to the plate x and y coords to stress test the quads with in-plane distortion.

However, If i apply even a very small random variation in the Z coords of the nodes, the results collapse again.

Were you seeing the same thing when looking into it previously?

[JWock82]

I did fix an issue with the local coordinate system early on when I was looking into this. That's probably why your model is working now.

As far as a random variation in the Z goes, I would expect the results to collapse if the plates are warped in any way. They are not meant to function that way. It takes 3 points to define a plane. If the 4th point is not in the same plane as the other 3 the plate is warped.

Try skewing your plates the opposite direction. I still seem to have trouble depending which direction the plates are skewed.

[JWock82]

The problem is fixed. There were two issues going on here. First, there was an error in defining the local coordinate system. I fixed that one early on, but I was still getting bad results in some instances. Yesterday I discovered that some membrane stiffness terms were being placed in the wrong positions in the quad element's global stiffness matrix. I've fixed the issue and tested the results. I've released the correction in version 0.0.28.

[bjhowie]

Excellent. They also seem to be more resilient to out-of-plane distortion now. Good work!

[JWock82]

That took way more effort than it should have to find. By far the toughest bug I've had to locate so far. Thank you for reporting this and for your help in testing the fix.