Bug fix for openings in meshes

PyNite/Mesh.py

@@ -545,20 +545,20 @@ def generate(self):
 
                 if ((round(opng.y_bott + opng.height, 10) >= round(top, 10))
                 and (round(opng.y_bott, 10) <= round(bott))
-                and (round(opng.x_left, 10) >= round(left, 10))
-                and (round(opng.x_left + opng.width, 10) <= round(right, 10))):
+                and (round(opng.x_left, 10) <= round(left, 10))
+                and (round(opng.x_left + opng.width, 10) >= round(right, 10))):
 
                     # Mark the element for deletion if it's not already marked
                     if element.Name not in element_del_list:
                         element_del_list.append(element.Name)
 
-        # Delete the nodes marked for deletion
-        for node_name in node_del_list:
-            del self.nodes[node_name]
-
         # Delete the elements marked for deletion
         for element_name in element_del_list:
             del self.elements[element_name]
+
+        # Delete the nodes marked for deletion
+        for node_name in node_del_list:
+            del self.nodes[node_name]
 
         # Find any remaining orphaned nodes around the perimeter of the mesh
         node_del_list = []

README.md

@@ -62,8 +62,9 @@ PyNite depends on the following packages:
 * sympy: Only needed if you want to view the derivations used to build PyNite.
 
 # What's New?
-v0.0.57 & 0.0.59
-* Bug fix for small openings that weren't being added to meshes.
+v0.0.57 & 0.0.60
+* Fixed a stubborn bug that wouldn't create openings if they didn't have a node inside them. This prevented openings from showing up in small meshes.
+* Added unit testing to help identify similar problems with mesh openings in the future.
 
 v0.0.56
 * Bug fix for orthotropic rectangular plates. The stiffness was slightly off on rectangular plates in version 0.0.55. Prior versions were not affected.

Testing/test_shear_wall.py

@@ -10,6 +10,7 @@
 from PyNite.Mesh import CylinderMesh, RectangleMesh
 import sys
 from io import StringIO
+from math import isclose
 
 class TestShearWalls(unittest.TestCase):
 
@@ -130,4 +131,104 @@ def test_cracked_rect_shear_wall(self):
         delta2 = V*H**3/(3*E*I) + 1.2*V*H/(G*A)
 
         # Check that the solution matches the theoretical solution within 2%
-        self.assertLess(abs(1 - delta1/delta2), 0.02, 'Failed cracked rect plate shear wall test.')
+        self.assertLess(abs(1 - delta1/delta2), 0.02, 'Failed cracked rect plate shear wall test.')
+
+    def test_shear_wall_openings(self):
+        # This example demonstrates how to analyze a shear wall with openings. It
+        # follows Section 10.5.3 of "Masonry Structures - Behavior and Design, 2nd
+        # Edition" by Robert G. Drysdale, Ahmad A. Hamid, and Lawrie R. Baker. The
+        # solution given in that text is obtained using an approximation method that
+        # isn't nearly as accurate as the finite element method, so some differences
+        # in the final results are expected.
+
+        # Set material properties for the wall (2 ksi masonry)
+        f_m = 2000        # Masonry compressive strength (psi)
+        E = 900*f_m/1000  # Masonry modulus of elasticity (ksi)
+        nu = 0.17         # Poisson's ratio for masonry
+
+        # Choose a desired mesh size. The program will try to stick to this as best as it can.
+        mesh_size = 6  # in
+
+        # Set the wall's dimensions
+        width = 26*12   # Wall overall width (in)
+        height = 16*12  # Wall overall height (in)
+        t = 8           # Masonry thickness (in)
+
+        # Generate the rectangular mesh
+        # The effects of cracked masonry can be modeled by adjusting the `ky_mod` factor. For this example
+        # uncracked masonry will be used to match the textbook problem.
+        mesh = RectangleMesh(mesh_size, width, height, t, E, nu, kx_mod=1, ky_mod=1, origin=[0, 0, 0],
+                            plane='XY', start_node='N1', start_element='R1', element_type='Rect')
+
+        # Add a 4' wide x 12' tall door opening to the mesh
+        mesh.add_rect_opening(name='Door 1', x_left=2*12, y_bott=0*12, width=4*12, height=12*12)
+
+        # Add a 4' wide x 4' tall window opening to the mesh
+        mesh.add_rect_opening(name='Window 1', x_left=8*12, y_bott=8*12, width=4*12, height=4*12)
+
+        # Add another 4' wide x 4' tall window opening to the mesh
+        mesh.add_rect_opening(name='Window 2', x_left=14*12, y_bott=8*12, width=4*12, height=4*12)
+
+        # Add another 4' wide x 12' tall door opening to the mesh
+        mesh.add_rect_opening(name='Door 2', x_left=20*12, y_bott=0*12, width=4*12, height=12*12)
+
+        # Generate the mesh now that we've defined all the openings
+        mesh.generate()
+
+        # Create a finite element model
+        model = FEModel3D()
+
+        # Add the mesh to the model
+        model.add_mesh(mesh)
+
+        # Shear at the top of the wall
+        V = 100  # kip
+
+        # The shear at the top of the wall will be distributed evently to all the
+        # nodes at the top of the wall. Determine how many nodes are at the top of the
+        # wall.
+        n = len([node for node in model.Nodes.values() if isclose(node.Y, height)])
+        v = V/n
+
+        # Add supports and loads to the nodes
+        for node in model.Nodes.values():
+
+            # Determine if the node is at the base of the wall
+            if isclose(node.Y, 0):
+                # Fix the base of the wall
+                model.def_support(node.Name, True, True, True, True, True, True)
+            # Determine if the node is at the top of the wall
+            elif isclose(node.Y, height):
+                # Add out-of-plane support (provided by the diaphragm)
+                model.def_support(node.Name, False, False, True, False, False, False)
+                # Add a seismic shear load to the top of the wall (applied by the diaphragm)
+                model.add_node_load(node.Name, 'FX', v, case='E')
+
+        # Add a load combination named 'Seismic' with a factor of 1.0 applied to any loads designated as
+        # 'E'.
+        model.add_load_combo('Seismic', {'E': 1.0})
+
+        # Analyze the model
+        model.analyze(log=True, check_statics=True)
+
+        # Render the model and plot the `Txy` shears.
+        # window = render_model(model, text_height=1, render_loads=True, deformed_shape=True,
+        #                       deformed_scale=200, color_map='Txy', scalar_bar=False,
+        #                       combo_name='Seismic', labels=False, screenshot='console')
+        from PyNite.Visualization import Renderer
+        renderer = Renderer(model)
+        renderer.combo_name = 'Seismic'
+        renderer.color_map = 'Txy'
+        renderer.annotation_size = 1
+        renderer.deformed_shape = True
+        renderer.deformed_scale = 200
+        renderer.scalar_bar = True
+        # renderer.render_model()
+        renderer.screenshot()
+
+        # Print the maximum displacement
+        d_max = max([node.DX['Seismic'] for node in model.Nodes.values()])
+        print('Max displacement: ', d_max, 'in')
+        print('Expected displacement from reference text: ', 7.623/E*t, 'in')
+        print('Wall rigidity: ', V/d_max, 'kips/in')
+        print('Expected rigidity from reference text: ', 1/7.623*E*t, 'kips/in')

setup.py

@@ -5,7 +5,7 @@
 
 setuptools.setup(
     name="PyNiteFEA",
-    version="0.0.59",
+    version="0.0.60",
     author="D. Craig Brinck, PE, SE",
     author_email="[email protected]",
     description="A simple elastic 3D structural finite element library for Python.",