Implement a nurbs cylindrical shell sector

meshpy/mesh_creation_functions/__init__.py

@@ -68,6 +68,7 @@
 # NURBS geometry functions
 from .nurbs_geometries import (
     create_nurbs_brick,
+    create_nurbs_cylindrical_shell_sector,
     create_nurbs_flat_plate_2d,
     create_nurbs_hemisphere_surface,
     create_nurbs_hollow_cylinder_segment_2d,
@@ -102,6 +103,7 @@
     "add_geomdl_nurbs_to_mesh",
     # NURBS geometry functions
     "create_nurbs_hollow_cylinder_segment_2d",
+    "create_nurbs_cylindrical_shell_sector",
     "create_nurbs_flat_plate_2d",
     "create_nurbs_brick",
     "create_nurbs_sphere_surface",

meshpy/mesh_creation_functions/nurbs_geometries.py

@@ -149,6 +149,107 @@ def create_nurbs_hollow_cylinder_segment_2d(
     return surf
 
 
+def create_nurbs_cylindrical_shell_sector(
+    radius, angle, length, *, n_ele_u=1, n_ele_v=1
+):
+    """Creates a patch of a surface of a 3-dimensional sector of a cylindrical
+    shell. The center of such cylindrical shell sector is located in [0, 0, 0]
+
+    Args
+    ----
+    radius: double
+        cylindrical shell radius
+    angle: double
+        angle of the cylindrical shell (radians)
+    length: double
+        length of the cylindrical shell
+    n_ele_u: int
+        number of elements in the parametric u-direction
+    n_ele_v: int
+        number of elements in the parametric v-direction
+
+
+    Return
+    ----
+    surf: geomdl object
+        geomdl object that contains the surface information
+    """
+
+    # Check the validity of the input values:
+    if (angle >= np.pi) or (angle < 0):
+        raise ValueError(
+            "The following algorithm for creating a cylindrical shell sector is only valid for 0 < angle <= pi."
+        )
+
+    # Create a NURBS surface instance
+    surf = NURBS.Surface()
+
+    # Set degrees
+    surf.degree_u = 2
+    surf.degree_v = 2
+
+    # Control points and set them to the surface
+    p_size_u = 3
+    p_size_v = 3
+
+    # Obtaining the control points
+    cp_1 = [-radius * np.sin(angle / 2), -length / 2, -radius * np.cos(angle / 2)]
+    cp_3 = [radius * np.sin(angle / 2), -length / 2, -radius * np.cos(angle / 2)]
+
+    # Calculating position of middle points. This is done by
+    # obtaining the tangents on the points cp_1 and cp_3 and
+    # calculating the intersection point between the tangents.
+    m_1 = -np.tan(-angle / 2)
+    m_3 = -np.tan(angle / 2)
+    b_1 = cp_1[2] + m_1 * cp_1[0]
+    b_3 = cp_3[2] + m_3 * cp_3[0]
+
+    inter_point_x = (b_3 - b_1) / (m_3 - m_1)
+    inter_point_z = m_1 * inter_point_x + b_1
+
+    # The intersection point is assigned to the middle points cp_2, cp_5 and cp_8
+    cp_2 = [inter_point_x, -length / 2, inter_point_z]
+
+    cp_4 = [-radius * np.sin(angle / 2), 0.0, -radius * np.cos(angle / 2)]
+    cp_5 = [inter_point_x, 0.0, inter_point_z]
+    cp_6 = [radius * np.sin(angle / 2), 0.0, -radius * np.cos(angle / 2)]
+
+    cp_7 = [
+        -radius * np.sin(angle / 2),
+        length / 2,
+        -radius * np.cos(angle / 2),
+    ]
+    cp_8 = [inter_point_x, length / 2, inter_point_z]
+    cp_9 = [radius * np.sin(angle / 2), length / 2, -radius * np.cos(angle / 2)]
+
+    ctrlpts = [cp_1, cp_4, cp_7, cp_2, cp_5, cp_8, cp_3, cp_6, cp_9]
+
+    weights = [
+        1.0,
+        1.0,
+        1.0,
+        np.cos(angle / 2),
+        np.cos(angle / 2),
+        np.cos(angle / 2),
+        1.0,
+        1.0,
+        1.0,
+    ]
+
+    t_ctrlptsw = compat.combine_ctrlpts_weights(ctrlpts, weights)
+
+    surf.ctrlpts_size_u = p_size_u
+    surf.ctrlpts_size_v = p_size_v
+    surf.ctrlptsw = t_ctrlptsw
+
+    surf.knotvector_u = [0.0, 0.0, 0.0, 1.0, 1.0, 1.0]
+    surf.knotvector_v = [0.0, 0.0, 0.0, 1.0, 1.0, 1.0]
+
+    do_uniform_knot_refinement_surface(surf, n_ele_u, n_ele_v)
+
+    return surf
+
+
 def create_nurbs_flat_plate_2d(width, length, *, n_ele_u=1, n_ele_v=1):
     """Creates a patch of a 2 dimensional flat plate.
 

tests/reference-files/test_nurbs_cylindrical_shell_sector.dat

@@ -0,0 +1,163 @@
+// -----------------------------------------------------------------------------
+// This input file was created with MeshPy.
+// Copyright (c) 2018-2024
+//     Ivo Steinbrecher
+//     Institute for Mathematics and Computer-Based Simulation
+//     Universitaet der Bundeswehr Muenchen
+//     https://www.unibw.de/imcs-en
+// -----------------------------------------------------------------------------
+-----------------------------------------------------------------------MATERIALS
+MAT 1 MAT_Struct_StVenantKirchhoff YOUNG 50 NUE 0.19 DENS 5.3e-07
+-----------------------------------------------------------STRUCTURE KNOTVECTORS
+NURBS_DIMENSION                       2
+BEGIN                                 NURBSPATCH
+ID                                    1
+NUMKNOTS                              9
+DEGREE                                2
+TYPE                                  Interpolated
+0.0
+0.0
+0.0
+0.25
+0.5
+0.75
+1.0
+1.0
+1.0
+NUMKNOTS                              8
+DEGREE                                2
+TYPE                                  Interpolated
+0.0
+0.0
+0.0
+0.3333333333333333
+0.6666666666666666
+1.0
+1.0
+1.0
+END                                   NURBSPATCH
+-------------------------------------------------------------DNODE-NODE TOPOLOGY
+NODE 30 DNODE 1
+NODE 6 DNODE 2
+NODE 25 DNODE 3
+NODE 1 DNODE 4
+-------------------------------------------------------------DLINE-NODE TOPOLOGY
+NODE 6 DLINE 1
+NODE 12 DLINE 1
+NODE 18 DLINE 1
+NODE 24 DLINE 1
+NODE 30 DLINE 1
+NODE 5 DLINE 2
+NODE 11 DLINE 2
+NODE 17 DLINE 2
+NODE 23 DLINE 2
+NODE 29 DLINE 2
+NODE 1 DLINE 3
+NODE 7 DLINE 3
+NODE 13 DLINE 3
+NODE 19 DLINE 3
+NODE 25 DLINE 3
+NODE 2 DLINE 4
+NODE 8 DLINE 4
+NODE 14 DLINE 4
+NODE 20 DLINE 4
+NODE 26 DLINE 4
+NODE 25 DLINE 5
+NODE 26 DLINE 5
+NODE 27 DLINE 5
+NODE 28 DLINE 5
+NODE 29 DLINE 5
+NODE 30 DLINE 5
+NODE 19 DLINE 6
+NODE 20 DLINE 6
+NODE 21 DLINE 6
+NODE 22 DLINE 6
+NODE 23 DLINE 6
+NODE 24 DLINE 6
+NODE 1 DLINE 7
+NODE 2 DLINE 7
+NODE 3 DLINE 7
+NODE 4 DLINE 7
+NODE 5 DLINE 7
+NODE 6 DLINE 7
+NODE 7 DLINE 8
+NODE 8 DLINE 8
+NODE 9 DLINE 8
+NODE 10 DLINE 8
+NODE 11 DLINE 8
+NODE 12 DLINE 8
+-------------------------------------------------------------DSURF-NODE TOPOLOGY
+NODE 1 DSURFACE 1
+NODE 2 DSURFACE 1
+NODE 3 DSURFACE 1
+NODE 4 DSURFACE 1
+NODE 5 DSURFACE 1
+NODE 6 DSURFACE 1
+NODE 7 DSURFACE 1
+NODE 8 DSURFACE 1
+NODE 9 DSURFACE 1
+NODE 10 DSURFACE 1
+NODE 11 DSURFACE 1
+NODE 12 DSURFACE 1
+NODE 13 DSURFACE 1
+NODE 14 DSURFACE 1
+NODE 15 DSURFACE 1
+NODE 16 DSURFACE 1
+NODE 17 DSURFACE 1
+NODE 18 DSURFACE 1
+NODE 19 DSURFACE 1
+NODE 20 DSURFACE 1
+NODE 21 DSURFACE 1
+NODE 22 DSURFACE 1
+NODE 23 DSURFACE 1
+NODE 24 DSURFACE 1
+NODE 25 DSURFACE 1
+NODE 26 DSURFACE 1
+NODE 27 DSURFACE 1
+NODE 28 DSURFACE 1
+NODE 29 DSURFACE 1
+NODE 30 DSURFACE 1
+---------------------------------------------------------------------NODE COORDS
+CP 1 COORD -0.5 -0.5 -0.866025403784 1.0
+CP 2 COORD -0.38799538113 -0.5 -0.930691300639 0.9665063509461097
+CP 3 COORD -0.133974596216 -0.5 -1 0.9330127018922194
+CP 4 COORD 0.133974596216 -0.5 -1 0.9330127018922194
+CP 5 COORD 0.38799538113 -0.5 -0.930691300639 0.9665063509461097
+CP 6 COORD 0.5 -0.5 -0.866025403784 1.0
+CP 7 COORD -0.5 -0.333333333333 -0.866025403784 1.0
+CP 8 COORD -0.38799538113 -0.333333333333 -0.930691300639 0.9665063509461098
+CP 9 COORD -0.133974596216 -0.333333333333 -1 0.9330127018922194
+CP 10 COORD 0.133974596216 -0.333333333333 -1 0.9330127018922194
+CP 11 COORD 0.38799538113 -0.333333333333 -0.930691300639 0.9665063509461098
+CP 12 COORD 0.5 -0.333333333333 -0.866025403784 1.0
+CP 13 COORD -0.5 0 -0.866025403784 1.0
+CP 14 COORD -0.38799538113 0 -0.930691300639 0.9665063509461098
+CP 15 COORD -0.133974596216 0 -1 0.9330127018922194
+CP 16 COORD 0.133974596216 0 -1 0.9330127018922194
+CP 17 COORD 0.38799538113 0 -0.930691300639 0.9665063509461098
+CP 18 COORD 0.5 0 -0.866025403784 1.0
+CP 19 COORD -0.5 0.333333333333 -0.866025403784 1.0
+CP 20 COORD -0.38799538113 0.333333333333 -0.930691300639 0.9665063509461097
+CP 21 COORD -0.133974596216 0.333333333333 -1 0.9330127018922194
+CP 22 COORD 0.133974596216 0.333333333333 -1 0.9330127018922194
+CP 23 COORD 0.38799538113 0.333333333333 -0.930691300639 0.9665063509461097
+CP 24 COORD 0.5 0.333333333333 -0.866025403784 1.0
+CP 25 COORD -0.5 0.5 -0.866025403784 1.0
+CP 26 COORD -0.38799538113 0.5 -0.930691300639 0.9665063509461097
+CP 27 COORD -0.133974596216 0.5 -1 0.9330127018922194
+CP 28 COORD 0.133974596216 0.5 -1 0.9330127018922194
+CP 29 COORD 0.38799538113 0.5 -0.930691300639 0.9665063509461097
+CP 30 COORD 0.5 0.5 -0.866025403784 1.0
+--------------------------------------------------------------STRUCTURE ELEMENTS
+1 WALLNURBS NURBS9 1 2 3 7 8 9 13 14 15  MAT 1 KINEM linear EAS none THICK 1.0 STRESS_STRAIN plane_strain GP 3 3
+2 WALLNURBS NURBS9 2 3 4 8 9 10 14 15 16  MAT 1 KINEM linear EAS none THICK 1.0 STRESS_STRAIN plane_strain GP 3 3
+3 WALLNURBS NURBS9 3 4 5 9 10 11 15 16 17  MAT 1 KINEM linear EAS none THICK 1.0 STRESS_STRAIN plane_strain GP 3 3
+4 WALLNURBS NURBS9 4 5 6 10 11 12 16 17 18  MAT 1 KINEM linear EAS none THICK 1.0 STRESS_STRAIN plane_strain GP 3 3
+5 WALLNURBS NURBS9 7 8 9 13 14 15 19 20 21  MAT 1 KINEM linear EAS none THICK 1.0 STRESS_STRAIN plane_strain GP 3 3
+6 WALLNURBS NURBS9 8 9 10 14 15 16 20 21 22  MAT 1 KINEM linear EAS none THICK 1.0 STRESS_STRAIN plane_strain GP 3 3
+7 WALLNURBS NURBS9 9 10 11 15 16 17 21 22 23  MAT 1 KINEM linear EAS none THICK 1.0 STRESS_STRAIN plane_strain GP 3 3
+8 WALLNURBS NURBS9 10 11 12 16 17 18 22 23 24  MAT 1 KINEM linear EAS none THICK 1.0 STRESS_STRAIN plane_strain GP 3 3
+9 WALLNURBS NURBS9 13 14 15 19 20 21 25 26 27  MAT 1 KINEM linear EAS none THICK 1.0 STRESS_STRAIN plane_strain GP 3 3
+10 WALLNURBS NURBS9 14 15 16 20 21 22 26 27 28  MAT 1 KINEM linear EAS none THICK 1.0 STRESS_STRAIN plane_strain GP 3 3
+11 WALLNURBS NURBS9 15 16 17 21 22 23 27 28 29  MAT 1 KINEM linear EAS none THICK 1.0 STRESS_STRAIN plane_strain GP 3 3
+12 WALLNURBS NURBS9 16 17 18 22 23 24 28 29 30  MAT 1 KINEM linear EAS none THICK 1.0 STRESS_STRAIN plane_strain GP 3 3

tests/testing_nurbs.py

@@ -37,13 +37,15 @@
 
 from meshpy import (
     InputFile,
+    InputSection,
     MaterialString,
     MaterialStVenantKirchhoff,
     Rotation,
 )
 from meshpy.mesh_creation_functions import (
     add_geomdl_nurbs_to_mesh,
     create_nurbs_brick,
+    create_nurbs_cylindrical_shell_sector,
     create_nurbs_flat_plate_2d,
     create_nurbs_hemisphere_surface,
     create_nurbs_hollow_cylinder_segment_2d,
@@ -86,6 +88,37 @@ def test_nurbs_hollow_cylinder_segment_2d(self):
         # Compare with the reference file
         compare_test_result(self, input_file.get_string(header=False))
 
+    def test_nurbs_cylindrical_shell_sector(self):
+        """Test the creation of a 3-dimensional cylindrical shell sector."""
+
+        # Create the surface of a quarter of a hollow cylinder
+        surf_obj = create_nurbs_cylindrical_shell_sector(
+            1, np.pi / 3, 1, n_ele_u=4, n_ele_v=3
+        )
+
+        # Create input file
+        input_file = InputFile()
+
+        # Add material
+        mat = MaterialStVenantKirchhoff(youngs_modulus=50, nu=0.19, density=5.3e-7)
+
+        # Create patch set
+        element_description = (
+            "KINEM linear EAS none THICK 1.0 STRESS_STRAIN plane_strain GP 3 3"
+        )
+
+        patch_set = add_geomdl_nurbs_to_mesh(
+            input_file,
+            surf_obj,
+            material=mat,
+            element_description=element_description,
+        )
+
+        input_file.add(patch_set)
+
+        # Compare with the reference file
+        compare_test_result(self, input_file.get_string(header=False))
+
     def test_nurbs_flat_plate_2d(self):
         """Test the creation of a two dimensional flat plate."""