From 3161cc74c888340e5c35364b9c41e1af72cef9cf Mon Sep 17 00:00:00 2001 From: Ivo Steinbrecher Date: Wed, 24 Jan 2024 09:39:23 +0100 Subject: [PATCH] Add some changes to tutorial --- cubitpy/cubitpy.py | 7 +- tutorial/tutorial.py | 210 +++++++++++++++++++++++++------------------ 2 files changed, 126 insertions(+), 91 deletions(-) diff --git a/cubitpy/cubitpy.py b/cubitpy/cubitpy.py index 36471b2..8be2e3b 100644 --- a/cubitpy/cubitpy.py +++ b/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. diff --git a/tutorial/tutorial.py b/tutorial/tutorial.py index a00edce..c325138 100644 --- a/tutorial/tutorial.py +++ b/tutorial/tutorial.py @@ -28,79 +28,130 @@ # 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 @@ -108,46 +159,48 @@ def cubitpy_tutorial(baci, output): # 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 @@ -155,21 +208,15 @@ def cubitpy_tutorial(baci, output): 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)