Add some changes to tutorial

cubitpy/cubitpy.py

@@ -415,9 +415,10 @@ def create_dat(self, dat_path, pre_exodus=False):
         """
 
         # Check if output path exists.
-        dat_dir = os.path.dirname(dat_path)
-        if not os.path.exists(dat_dir):
-            raise ValueError("Path {} does not exist!".format(dat_dir))
+        if os.path.isabs(dat_path):
+            dat_dir = os.path.dirname(dat_path)
+            if not os.path.exists(dat_dir):
+                raise ValueError("Path {} does not exist!".format(dat_dir))
 
         if pre_exodus:
             # Create the dat file.

tutorial/tutorial.py

@@ -28,148 +28,195 @@
 # SOFTWARE.
 # -----------------------------------------------------------------------------
 """
-This script contains a tutorial for Cubitpy, following the Cubit tutorial. Most basic functionality is covered
-by this tutorial. For more information have a closer look at the test cases,
-as they cover all functionality.
+This file contains a tutorial for CubitPy, following the Cubit Step-By-Step tutorial
+from the cubit documentation. Most basic functionality is covered by this tutorial.
+For more information have a closer look at the test cases, as they cover all relevant
+functionality.
 """
 
-# Meshpy modules.
-import subprocess
-import cubitpy
-
-
-# Import cubitpy module.
+import numpy as np
 from cubitpy import CubitPy, cupy
-from cubitpy.mesh_creation_functions import create_brick
+from cubitpy.cubit_utility import get_surface_center
+
 
-from meshpy import mesh
+def cubit_step_by_step_tutorial_cli(input_file_path, *, display=True):
+    """This tutorial follows the Cubit step-by-step tutorial, which can be found
+    in the cubit documentation"""
 
+    # Geometry parameters
+    brick_size = 10.0
+    hole_radius = 3.0
 
-def cubitpy_tutorial(baci, output):
-    # In the first step, a cubitpy object is created to connect the
-    # code with the wrapper for Cubit. This also allows for direct insertion
-    # of commands, which will be shown at a later stage.
+    # The first step is to setup a CubitPy object. This object is derived from the
+    # cubit python interface and provides all functionality of the direct cubit
+    # python interface and also adds some additional functionality to make the
+    # creation of BACI input files easier.
     cubit = CubitPy()
 
-    # Once the cubitpy object is initialized, we can create our first brick
-    # object. For that we pass the cubitpy object, the dimensions of the brick
-    # in x, y, and z directions. the type of elements for the meshing phase,
-    # the meshing interval in [x, y, z] directions, and the element conditions
-    # (linear or nonlinear kinematics).
-    brick_1 = create_brick(cubit, 10, 10, 10, mesh=False)
+    # Once the CubitPy object is initialized, we can create our first brick
+    # object. To do so we use the cubit python interface function `brick`.
+    # This function returns a python object referencing the brick.
+    brick = cubit.brick(brick_size, brick_size, brick_size)
+
+    # Alternatively the brick could be created via the command line interface (CLI)
+    # (not all cubit commands are available through python functions).
+    # The `cubit.cmd` method allows to pass any cubit command as string.
+    #
+    # cubit.cmd(f"brick x {brick_size}")
 
     # The cube can be shown on your display using the following command
     # (You may have to click on the coordinate icon or refresh the display
     # for the item to show)
-    cubit.display_in_cubit()
+    if display:
+        cubit.display_in_cubit()
+
+    # You can also show the cube with ID labels on the geometry. The cupy object is
+    # a container holding all of cubits enums, such as the geometry type.
+    if display:
+        cubit.display_in_cubit(labels=[cupy.geometry.curve, cupy.geometry.surface])
+
+    # Now you must create the cylinder which will be used to cut the hole
+    # from the brick. This is accomplished with the CLI.
+    # (Alternatively we could also use the cubit python command `cubit.cylinder`)
+    cubit.cmd(f"create cylinder height {brick_size*1.2} radius {hole_radius}")
 
-    # Now you must form the cylinder which will be used to cut the hole
-    # from the brick. This is accomplished with the command:
-    cubit.cmd("create cylinder height 12 radius 3")
+    # Executing a command via the CLI does not return the geometry object. To get the
+    # cylinder object we can write
+    cylinder = cubit.volume(cubit.get_last_id(cupy.geometry.volume))
 
     # At this point you will see both a brick and a cylinder appear
     # in the CUBIT display window.
-    cubit.display_in_cubit()
+    if display:
+        cubit.display_in_cubit()
 
     # Now, the cylinder can be subtracted from the brick to form the hole in the block.
-    # Note that both original volumes are deleted in the Boolean operation and replaced
-    # with a new volume (with an id of 1) which is the result of the Boolean operation Subtract .
-    cubit.cmd("subtract 2 from 1")
-
-    # Now we start meshing by setting the intervals for the brick and defining the mesh
+    # Note that the after the subtraction, the cylinder volume is removed and the brick
+    # object now points to the result of the subtraction.
+    # For the subtraction, we need the ids of the respective volumes. There are two ways
+    # to obtain the IDs:
+    #   - Open Cubit and manually write down the volume IDs (in this case 1 and 2)
+    #         cubit.cmd("subtract 2 from 1")
+    #     This has the drawback that if something in the numbering changes, i.e., due to a
+    #     Cubit internal change or if more items are inserted before the current ones, that
+    #     the code has to be changed.
+    #   - The more sustainable way is to get the ID from the cubit objects directly, which can
+    #     be done with the `id()` method:
+    cubit.cmd(f"subtract volume {cylinder.id()} from volume {brick.id()}")
+
+    # Now we see, that the cylinder is removed from the brick
+    if display:
+        cubit.display_in_cubit()
+
+    # We now start with the meshing. For more details about this have a look at the
+    # Cubit tutorial.
+
+    # First setting the intervals for the brick and defining the mesh
     # size for the volume. we begin by setting an overall volume size interval.
     # Since the brick is 10 units in length on a side, this specifies that each
     # straight curve is to receive approximately 10 mesh elements.
-    cubit.cmd("volume 1 size 1.0")
-
-    # In order to better resolve the hole in the middle of the top surface,
-    # we set a smaller size for the curve bounding this hole. The IDs for the
-    # curves can be either obtained from the display or by selection methods
-    cubit.cmd(
-        """curve 15 interval size 0.5
-    curve 11 interval 5"""
-    )
+    cubit.cmd(f"volume {brick.id()} size 1.0")
 
-    # Now we can begin with meshing the surface with the hole. for that we
-    # define the mesh scheme by setting to pave. Again, the ID of the surface
-    # is obtained from the display or selection methods
-    cubit.cmd(
-        """surface 11 scheme pave
-    mesh surface 11"""
-    )
+    # We will use a sweep mesh, so we will first mesh one surface of the body. We want to
+    # first mesh the top surface. This surface can be selected by looking for the only
+    # surface with a normal vector in positive z-direction.
+    for surface in brick.surfaces():
+        # We use the following helper function to get the center of the surface
+        center = get_surface_center(surface)
 
-    # We can now see the meshed surface
-    # in the CUBIT display window.
-    cubit.display_in_cubit()
+        # Get the normal at the center (also works for the surface with the hole in it)
+        normal = surface.normal_at(center)
+
+        # Check if the normal is in positive z-direction
+        if np.dot([0, 0, 1], normal) > (1 - 1e-10):
+            mesh_surface = surface
+
+    # In order to better resolve the hole in the middle of the surface,
+    # we set a smaller size for the curve bounding this hole. We can either get the
+    # curve ID from the GUI, or select them in python.
+    # Here we select them in python. To do so, we loop over all curves of the surface and
+    # check if the curve lies on the surface of the cylinder and on the surface we want.
+    for curve in mesh_surface.curves():
+        # Position in the middle of the curve
+        curve_pos = curve.position_from_fraction(0.5)
+
+        # Radius w.r.t. the cylinder axis
+        radius = np.linalg.norm(curve_pos[:2])
+
+        if np.abs(radius - hole_radius) < 1e-10:
+            # The curve lies on the cylinder radius, now set the meshing interval
+            cubit.cmd(f"curve {curve.id()} interval size 0.5")
+
+    # Now we can mesh the surface
+    cubit.cmd(f"mesh surface {mesh_surface.id()}")
+
+    # We can now see the meshed surface in the CUBIT display window.
+    if display:
+        cubit.display_in_cubit()
 
     # The volume mesh can now be generated. Again, the first step is to specify
     # the type of meshing scheme to be used and the second step is to issue the
     # order to mesh. In certain cases, the scheme can be determined by CUBIT
     # automatically. For sweepable volumes, the automatic scheme detection
     # algorithm also identifies the source and target surfaces of the
     # sweep automatically.
+    cubit.cmd(f"volume {brick.id()} scheme auto")
+    cubit.cmd(f"mesh volume {brick.id()}")
 
-    cubit.cmd(
-        """volume 1 scheme auto
-    mesh volume 1"""
-    )
-
-    # We can now see the meshed volume
-    # in the CUBIT display window.
-    cubit.display_in_cubit()
+    # We can now see the meshed volume in the CUBIT display window.
+    if display:
+        cubit.display_in_cubit()
 
     # We can now define the boundary conditions, from supports to loading.
     # We will place a load on the top side of the cube, and fix the bottom.
     # It is worth noting that using coordinates to select the geometry
     # is preferred, as the element ID may change with different
-    # versions / runs of Cubit.
+    # versions / runs of Cubit (as mentioned earlier).
+    # The `cubit.add_node_set` method tracks the defined node sets, which allows
+    # an automated creation of the BCs in the input file (no need for a `.bc` file)
     cubit.add_node_set(
         cubit.group(add_value="add surface with y_coord < -4.99"),
         name="fix",
         bc_type=cupy.bc_type.dirichlet,
-        bc_description="NUMDOF 3 ONOFF 1 1 1 VAL 0 0 0 FUNCT 0 0 0 TAG monitor_reaction",
+        bc_description="NUMDOF 3 ONOFF 1 1 1 VAL 0 0 0 FUNCT 0 0 0",
     )
-
     cubit.add_node_set(
         cubit.group(add_value="add surface with y_coord > 4.99 "),
         name="load",
         bc_type=cupy.bc_type.neumann,
-        bc_description="NUMDOF 3 ONOFF 0 1 0 VAL 0 1E-1 0 FUNCT 0 1 0",
+        bc_description="NUMDOF 3 ONOFF 0 1 0 VAL 0 0.1 0 FUNCT 0 1 0",
         geometry_type=cupy.geometry.vertex,
     )
 
+    # Finally we have to set the element blocks.
+    cubit.add_element_type(brick.volumes()[0], el_type=cupy.element_type.hex8)
+
     # We can view the created mesh along with the node sets we defined
     # earlier with the boundary conditions.
-    cubit.display_in_cubit()
+    if display:
+        cubit.display_in_cubit()
 
     # Set the head string.
     cubit.head = """
     -------------------------------------------------------------------PROBLEM TYP
     PROBLEMTYP                      Structure
     ----------------------------------------------------------------------------IO
-    OUTPUT_BIN                      yes
+    OUTPUT_BIN                      no
     STRUCT_DISP                     yes
     FILESTEPS                       1000
     VERBOSITY                       Standard
     STRUCT_STRAIN                   gl
     STRUCT_STRESS                   cauchy
     OUTPUT_SPRING                   Yes
     WRITE_INITIAL_STATE             yes
-    --------------------------------------------------------IO/MONITOR STRUCTURE DBC
-    PRECISION_FILE                        16
-    PRECISION_SCREEN                      5
-    FILE_TYPE                             csv
-    INTERVAL_STEPS                        1
     ---------------------------------------------------------IO/RUNTIME VTK OUTPUT
-    OUTPUT_DATA_FORMAT              ascii
+    OUTPUT_DATA_FORMAT              binary
     INTERVAL_STEPS                  1
     EVERY_ITERATION                 no
     -----------------------------------------------IO/RUNTIME VTK OUTPUT/STRUCTURE
     OUTPUT_STRUCTURE                yes
     DISPLACEMENT                    yes
-    ELEMENT_OWNER                   no
+    ELEMENT_OWNER                   yes
     STRESS_STRAIN                   yes
-    GAUSS_POINT_DATA_OUTPUT_TYPE    none
     ----------------------------------------------------------------------SOLVER 1
     NAME                            Structure_Solver
     SOLVER                          Superlu
@@ -192,29 +239,16 @@ def cubitpy_tutorial(baci, output):
     Inner Iteration                 = No
     Outer Iteration StatusTest      = No
     ---------------------------------------------------------------------MATERIALS
-     MAT 1 MAT_Struct_StVenantKirchhoff  YOUNG 1.0E1  NUE 0.48 DENS 0 
+     MAT 1 MAT_Struct_StVenantKirchhoff  YOUNG 1.0E1  NUE 0.3 DENS 0 
     --------------------------------------------------------------------------FUNCT1
     SYMBOLIC_FUNCTION_OF_TIME t
     """
 
     # Write the input file.
-    cubit.create_dat(output + "/CubitTutorial.dat")
-
-    # Run Baci - get output
-    subprocess.run(
-        [
-            baci,
-            output + "/CubitTutorial.dat",
-            output + "/Test",
-        ]
-    )
+    cubit.create_dat(input_file_path)
 
 
 if __name__ == "__main__":
     """Execution part of script."""
 
-    # Set your baci-release directory and you designated output directory
-    baci = "/home/$USERNAME/Desktop/Baci/BuildBaci/baci-release"
-    output = "/home/$USERNAME/Desktop/CubitPy/Test"
-
-    cubitpy_tutorial(baci, output)
+    cubit_step_by_step_tutorial_cli("cubit_tutorial.dat", display=False)