Tutorial: Gmsh+CalculiX

Tutorial: Gmsh+CalculiX#

This tutorial mirrors the Abaqus-based core tutorials with a fully open-source, conda-forge workflow. The geometry is created and meshed with Gmsh. The finite element simulation is performed with CalculiX. Finally, the data is extracted with ccx2paraview and meshio for scripted quantity of interest (QoI) calculations and post-processing.

References#

Environment#

Warning

The Gmsh+CalculiX tutorial requires a different compute environment than the other tutorials. The following commands create a dedicated environment for the use of this tutorial. You can also use your existing tutorial environment with the conda install command while your environment is active.

All of the software required for this tutorial can be installed in a Conda environment with the Conda package manager. See the Conda installation and Conda environment management documentation for more details about using Conda.

  1. Create the Gmsh+ tutorial environment if it doesn’t exist

    $ conda create --name waves-gmsh-env --channel conda-forge waves 'scons>=4.6' python-gmsh calculix 'ccx2paraview>=3.2' meshio

  2. Activate the environment

    $ conda activate waves-gmsh-env

Directory Structure#

  1. Create and change to a new project root directory to house the tutorial files if you have not already done so. For example

$ mkdir -p ~/waves-tutorials
$ cd ~/waves-tutorials
$ pwd
/home/roppenheimer/waves-tutorials

  1. Create a new tutorial_gmsh directory with the waves fetch command below

$ pwd
/home/roppenheimer/waves-tutorials
$ waves fetch --destination tutorial_gmsh tutorials/tutorial_gmsh
$ ls tutorial_gmsh
SConscript  SConstruct  environment.yml  post_processing.py  rectangle.py  rectangle_compression.inp  strip_heading.py  vtu2xarray.py

  1. Make the new tutorial_gmsh directory the current working directory

$ pwd
/home/roppenheimer/waves-tutorials
$ cd tutorial_gmsh
$ pwd
/home/roppenheimer/waves-tutorials/tutorial_gmsh
$ ls
SConscript  SConstruct  environment.yml  post_processing.py  rectangle.py  rectangle_compression.inp  strip_heading.py  vtu2xarray.py

SConscript#

  1. Review the SConscript workflow configuration file.

The structure is sufficiently different from the core tutorials that a diff view is not as useful. Instead the contents of the new SConscript files are duplicated below.

waves-tutorials/tutorial_gmsh/SConscript

 1#! /usr/bin/env python
 2import waves
 3
 4Import("env")
 5
 6# Collect the target nodes to build concise target alias(es)
 7artifacts = []
 8workflow = []
 9
10# Geometry, Partition, Mesh
11artifacts.extend(
12    env.PythonScript(
13        target=["rectangle_gmsh.inp"],
14        source=["rectangle.py"],
15        subcommand_options="--output-file=${TARGET.abspath}",
16        action_suffix=env["action_suffix"],
17    )
18)
19
20artifacts.extend(
21    env.PythonScript(
22        target=["rectangle_mesh.inp"],
23        source=["strip_heading.py", "rectangle_gmsh.inp"],
24        subcommand_options="--input-file=${SOURCES[1].abspath} --output-file=${TARGET.abspath}",
25        action_suffix=env["action_suffix"],
26    )
27)
28
29# CalculiX Solve
30artifacts.extend(
31    env.CalculiX(
32        target=[f"rectangle_compression.{suffix}" for suffix in ("frd", "dat", "sta", "cvg", "12d")],
33        source=["rectangle_compression.inp"],
34        action_suffix=env["action_suffix"],
35    )
36)
37
38# Extract
39artifacts.extend(
40    env.Command(
41        target=["rectangle_compression.vtu", "rectangle_compression.vtu.stdout"],
42        source=["rectangle_compression.frd"],
43        action=["cd ${TARGET.dir.abspath} && ccx2paraview ${SOURCES[0].abspath} vtu ${action_suffix}"],
44        action_suffix=env["action_suffix"],
45    )
46)
47workflow.extend(
48    env.PythonScript(
49        target=["rectangle_compression.h5"],
50        source=["vtu2xarray.py", "rectangle_compression.vtu"],
51        subcommand_options="--input-file ${SOURCES[1].abspath} --output-file ${TARGET.abspath}",
52        action_suffix=env["action_suffix"],
53    )
54)
55
56# Post-processing
57script_options = "--input-file ${SOURCES[1:].abspath}"
58script_options += " --output-file ${TARGET.file} --x-units mm/mm --y-units MPa"
59workflow.extend(
60    env.PythonScript(
61        target=["stress_strain_comparison.pdf", "stress_strain_comparison.csv"],
62        source=["post_processing.py", "rectangle_compression.h5"],
63        subcommand_options=script_options,
64    )
65)
66artifacts.extend(workflow)
67
68# Collector alias named after the model simulation
69env.Alias("rectangle", workflow)
70
71if not env["unconditional_build"] and (not env["CCX_PROGRAM"] or not env["ccx2paraview"]):
72    print(
73        "Program 'CalculiX (ccx)' or 'ccx2paraview' was not found in construction environment. "
74        "Ignoring 'rectangle' target(s)"
75    )
76    Ignore([".", "rectangle"], workflow + artifacts)

Gmsh Python script#

  1. Review the Gmsh Python script rectangle.py.

The Gmsh script differs from the Abaqus journal files introduced in Tutorial 02: Partition and Mesh. Besides the differences in Abaqus and Gmsh commands, one major difference between the Abaqus and Gmsh scripts is the opportunity to use Python 3 with Gmsh, where Abaqus journal files must use the Abaqus controlled installation of Python 2. The Gmsh script creates the geometry, partition (sets), and mesh in a single script because the output formats used by Gmsh can not contain both the geometry definition and the mesh information.

waves-tutorials/tutorial_gmsh/rectangle.py

  1import typing
  2import pathlib
  3import argparse
  4
  5import numpy
  6import gmsh
  7
  8
  9script_name = pathlib.Path(__file__)
 10# Set default parameter values
 11default_output_file = pathlib.Path(script_name.with_suffix(".msh").name)
 12default_width = 1.0
 13default_height = 1.0
 14default_global_seed = 1.0
 15
 16
 17def main(
 18    output_file: pathlib.Path = default_output_file,
 19    width: float = default_width,
 20    height: float = default_height,
 21    global_seed: float = default_global_seed,
 22) -> None:
 23    """Create a simple rectangle geometry.
 24
 25    This script creates a simple Gmsh model with a single rectangle part.
 26
 27    :param output_file: The output file for the Gmsh model. Extension must match a supported Gmsh file type.
 28    :param width: The rectangle width
 29    :param height: The rectangle height
 30    :param global_seed: The global mesh seed size
 31
 32    :returns: writes ``output_file``
 33    """
 34    gmsh.initialize()
 35    gmsh.logger.start()
 36    gmsh.model.add("rectangle")
 37
 38    # Geometry
 39    rectangle_tag = gmsh.model.occ.addRectangle(0, 0, 0, width, height)
 40    gmsh.model.occ.synchronize()
 41
 42    # Partition
 43    # TODO: Move to separate function/script
 44    rectangle_group_tag = gmsh.model.addPhysicalGroup(2, [rectangle_tag], name="rectangle")
 45
 46    tolerance = 0.01 * min(width, height)
 47    bottom_left = get_entities_at_coordinates((0.0, 0.0, 0.0), 0)
 48    top_right = get_entities_at_coordinates((width, height, 0.0), 0)
 49    # fmt: off
 50    top = gmsh.model.getEntitiesInBoundingBox(
 51        0.0    - tolerance,  # noqa: E221 X-min
 52        height - tolerance,  # noqa: E221 Y-min
 53        0.0    - tolerance,  # noqa: E221 Z-min
 54        width  + tolerance,  # noqa: E221 X-max
 55        height + tolerance,  # noqa: E221 Y-max
 56        0.0    + tolerance,  # noqa: E221 Z-max
 57        0  # Entity dimension: points
 58    )
 59    bottom = gmsh.model.getEntitiesInBoundingBox(
 60        0.0    - tolerance,  # noqa: E221 X-min
 61        0.0    - tolerance,  # noqa: E221 Y-min
 62        0.0    - tolerance,  # noqa: E221 Z-min
 63        width  + tolerance,  # noqa: E221 X-max
 64        0.0    + tolerance,  # noqa: E221 Y-max
 65        0.0    + tolerance,  # noqa: E221 Z-max
 66        0  # Entity dimension: points
 67    )
 68    left = gmsh.model.getEntitiesInBoundingBox(
 69        0.0    - tolerance,  # noqa: E221 X-min
 70        0.0    - tolerance,  # noqa: E221 Y-min
 71        0.0    - tolerance,  # noqa: E221 Z-min
 72        0.0    + tolerance,  # noqa: E221 X-max
 73        height + tolerance,  # noqa: E221 Y-max
 74        0.0    + tolerance,  # noqa: E221 Z-max
 75        0  # Entity dimension: points
 76    )
 77    # fmt: on
 78    gmsh.model.addPhysicalGroup(0, tags_from_dimTags(bottom_left), name="bottom_left")
 79    gmsh.model.addPhysicalGroup(0, tags_from_dimTags(top_right), name="top_right")
 80    gmsh.model.addPhysicalGroup(0, tags_from_dimTags(top), name="top")
 81    gmsh.model.addPhysicalGroup(0, tags_from_dimTags(bottom), name="bottom")
 82    gmsh.model.addPhysicalGroup(0, tags_from_dimTags(left), name="left")
 83
 84    # Mesh
 85    # TODO: Move to separate function/script
 86    points = gmsh.model.getEntities(0)
 87    gmsh.model.mesh.setSize(points, global_seed)
 88    gmsh.model.mesh.generate(2)
 89    gmsh.model.mesh.recombine()
 90
 91    gmsh.option.setNumber("Mesh.SaveGroupsOfElements", 1)
 92    gmsh.option.setNumber("Mesh.SaveGroupsOfNodes", 1)
 93
 94    gmsh.write(str(output_file))
 95    gmsh.logger.stop()
 96    gmsh.finalize()
 97
 98
 99def tags_from_dimTags(dimTags: typing.List[typing.Tuple[int, int]]) -> typing.List[int]:
100    """Return tags from Gmsh entity ``dimTags`` list of tuples
101
102    :returns: list of tags
103    """
104    return [dimTag[1] for dimTag in dimTags]
105
106
107def get_entities_at_coordinates(
108    coordinates: typing.Tuple[float, float, float],
109    dimension: int,
110    tolerance: float = 1.0e-6,
111) -> typing.List[typing.Tuple[int, int]]:
112    """Return Gmsh ``dimTags`` of entities of dimension within bounding box determined by coordinates and tolerance
113
114    :param coordinates: 3D coordinates (X, Y, Z) for center of bounding box
115    :param dimension: Return entities matching dimension
116    :param tolerance: Distance from center of bounding box to edges of bounding box
117
118    :returns: Gmsh entity ``dimTags``
119    """
120    if len(coordinates) != 3:
121        raise RuntimeError("coordinates must have length 3")
122    center = numpy.array(coordinates)
123    tolerance_array = numpy.array([tolerance, tolerance, tolerance])
124    min_coordinate = center - tolerance_array
125    max_coordinate = center + tolerance_array
126    return gmsh.model.getEntitiesInBoundingBox(*min_coordinate, *max_coordinate, dimension)
127
128
129def get_parser():
130    prog = f"python {script_name.name} "
131    cli_description = "Create a simple rectangle geometry and write an ``output_file``.msh Gmsh model file."
132    parser = argparse.ArgumentParser(description=cli_description, prog=prog)
133    parser.add_argument(
134        "--output-file",
135        type=pathlib.Path,
136        default=default_output_file,
137        # fmt: off
138        help="The output file for the Gmsh model. Extension must match a supported Gmsh file type, e.g. "
139             "``output_file``.msh (default: %(default)s",
140        # fmt: on
141    )
142    parser.add_argument(
143        "--width",
144        type=float,
145        default=default_width,
146        help="The rectangle width",
147    )
148    parser.add_argument(
149        "--height",
150        type=float,
151        default=default_height,
152        help="The rectangle height",
153    )
154    parser.add_argument(
155        "--global-seed",
156        type=float,
157        default=default_global_seed,
158        help="The global mesh seed size (default: %(default)s)",
159    )
160    return parser
161
162
163if __name__ == "__main__":
164    parser = get_parser()
165    args = parser.parse_args()
166    main(
167        output_file=args.output_file,
168        width=args.width,
169        height=args.height,
170        global_seed=args.global_seed,
171    )

The other Python scripts will not be discussed in detail. They are small content handling and data conversion utilities or discussed in greater depth in other tutorials.

CalculiX Input File#

  1. Create or review the CalculiX input file from the contents below

waves-tutorials/tutorial_gmsh/rectangle_compression.inp

 1**
 2*HEADING
 3Compressing a rectangle
 4**
 5*INCLUDE, INPUT=rectangle_mesh.inp
 6*SOLID SECTION, ELSET=RECTANGLE, MATERIAL=mock
 7**
 8*MATERIAL, NAME=mock
 9*ELASTIC
10100, .3
11*DENSITY
12.27000E-08,
13*SPECIFIC HEAT
148.96E+08
15*CONDUCTIVITY
16132, 200
17*EXPANSION,ZERO=293.15
1823.58E-06
19**
20*************************************************************************
21**  STEP-1 LOAD STEP
22*************************************************************************
23*STEP, NLGEOM=NO, INC=100, AMPLITUDE=RAMP
24*STATIC
25.1, 1.0, 0.01, 0.5
26**
27*BOUNDARY
28bottom_left, 1,3
29bottom, 2,2
30top,2,2,-0.01
31**
32*EL FILE
33E, S
34*NODE FILE
35U
36**
37*END STEP

SConstruct#

Note that CalculiX differs from other solvers in the tutorials. CalculiX is deployed as a Conda package and is available in the launching Conda environment. It is still good practice to check if the executable is available and provide helpful feedback to developers about the excutable status and workflow configuration.

The structure has changed enough from the core tutorials that a diff view is not as useful. Instead the contents of the SConstruct file is duplicated below.

waves-tutorials/tutorial_gmsh/SConstruct

#! /usr/bin/env python
import os
import platform

import waves

# Add build directory CLI arg
AddOption(
    "--build-dir",
    dest="build_dir",
    default="build",
    nargs=1,
    type="string",
    action="store",
    metavar="DIR",
    help="SCons build (variant) root directory. Relative or absolute path. (default: '%default')",
)
AddOption(
    "--unconditional-build",
    dest="unconditional_build",
    default=False,
    action="store_true",
    help="Boolean to force building of conditionally ignored targets. (default: '%default')",
)
AddOption(
    "--verbose-tasks",
    dest="verbose_tasks",
    default=False,
    action="store_true",
    help="Boolean to 'tee' output for interactive and verbose task execution. (default: '%default')",
)

# Inherit user's full environment
env = waves.scons_extensions.WAVESEnvironment(
    ENV=os.environ.copy(),
    build_dir=GetOption("build_dir"),
    unconditional_build=GetOption("unconditional_build"),
    verbose_tasks=GetOption("verbose_tasks"),
)

# Action suffix
system = platform.system().lower()
if env["verbose_tasks"]:
    if system == "windows":  # Assume PowerShell
        env["action_suffix"] = "$(| Tee-Object -FilePath ${TARGETS[-1].abspath}$)"
    else:  # *Nix style tee
        env["action_suffix"] = "$(2>&1 | tee ${TARGETS[-1].abspath}$)"
else:
    # Default WAVES action suffix
    env["action_suffix"] = "$(> ${TARGETS[-1].abspath} 2>&1$)"

# Always copy (no sym-links) when duplicating
env.SetOption("duplicate", "copy")

# Always print failed task *.stdout files
env.PrintBuildFailures(print_stdout=True)

# Empty defaults list to avoid building all simulation targets by default
env.Default()

# Find required programs for conditional target ignoring and absolute path for use in target actions
env["CCX_PROGRAM"] = env.AddProgram(["ccx"])
env["ccx2paraview"] = env.AddProgram(["ccx2paraview"])

# Abaqus input file implicit dependency scanner
# Works for CalculiX because CalculiX uses the Abaqus include keyword semantics and input file extension
env.Append(SCANNERS=waves.scons_extensions.abaqus_input_scanner())

# Call SConscript file
SConscript("SConscript", variant_dir=env["build_dir"], exports={"env": env}, duplicate=True)

# List all aliases in help message.
# This must come *after* all expected Alias definitions and SConscript files.
env.ProjectHelp()

Build Targets#

  1. Build the new targets

$ pwd
/path/to/waves-tutorials/tutorial_gmsh
$ scons fierro
scons: Reading SConscript files ...
Checking whether 'ccx' program exists.../home/roppenheimer/anaconda3/envs/waves-gmsh-env/bin/ccx
Checking whether 'ccx2paraview' program exists.../home/roppenheimer/anaconda3/envs/waves-gmsh-env/bin/ccx2paraview
scons: done reading SConscript files.
scons: Building targets ...
cd /home/roppenheimer/waves-tutorials/tutorial_gmsh/build && python /home/roppenheimer/waves-tutorials/tutorial_gmsh/build/rectangle.py --output-file=/home/roppenheimer/waves-tutorials/tutorial_gmsh/build/rectangle_gmsh.inp > /home/roppenheimer/waves-tutorials/tutorial_gmsh/build/rectangle_gmsh.inp.stdout 2>&1
cd /home/roppenheimer/waves-tutorials/tutorial_gmsh/build && python /home/roppenheimer/waves-tutorials/tutorial_gmsh/build/strip_heading.py --input-file=/home/roppenheimer/waves-tutorials/tutorial_gmsh/build/rectangle_gmsh.inp --output-file=/home/roppenheimer/waves-tutorials/tutorial_gmsh/build/rectangle_mesh.inp > /home/roppenheimer/waves-tutorials/tutorial_gmsh/build/rectangle_mesh.inp.stdout 2>&1
cd /home/roppenheimer/waves-tutorials/tutorial_gmsh/build && /home/roppenheimer/anaconda3/envs/waves-gmsh-env/bin/ccx -i rectangle_compression > /home/roppenheimer/waves-tutorials/tutorial_gmsh/build/rectangle_compression.frd.stdout 2>&1
cd /home/roppenheimer/waves-tutorials/tutorial_gmsh/build && ccx2paraview /home/roppenheimer/waves-tutorials/tutorial_gmsh/build/rectangle_compression.frd > /home/roppenheimer/waves-tutorials/tutorial_gmsh/build/rectangle_compression.vtu.stdout 2>&1
cd /home/roppenheimer/waves-tutorials/tutorial_gmsh/build && python /home/roppenheimer/waves-tutorials/tutorial_gmsh/build/vtu2xarray.py --input-file /home/roppenheimer/waves-tutorials/tutorial_gmsh/build/rectangle_compression.vtu --output-file /home/roppenheimer/waves-tutorials/tutorial_gmsh/build/rectangle_compression.h5 > /home/roppenheimer/waves-tutorials/tutorial_gmsh/build/rectangle_compression.h5.stdout 2>&1
cd /home/roppenheimer/waves-tutorials/tutorial_gmsh/build && python /home/roppenheimer/waves-tutorials/tutorial_gmsh/build/post_processing.py --input-file /home/roppenheimer/waves-tutorials/tutorial_gmsh/build/rectangle_compression.h5 --output-file stress_strain_comparison.pdf --x-units mm/mm --y-units MPa > /home/roppenheimer/waves-tutorials/tutorial_gmsh/build/stress_strain_comparison.pdf.stdout 2>&1
scons: done building targets.

Output Files#

Explore the contents of the build directory using the find command against the build directory, as shown below.

$ pwd
/home/roppenheimer/waves-tutorials/tutorial_cubit
$ find build -type f
build/rectangle_compression.h5
build/rectangle_compression.inp
build/rectangle_compression.h5.stdout
build/rectangle_compression.sta
build/rectangle_compression.dat
build/rectangle_mesh.inp
build/rectangle_compression.vtu.stdout
build/rectangle_compression.12d
build/spooles.out
build/rectangle_gmsh.inp
build/rectangle_compression.cvg
build/stress_strain_comparison.pdf.stdout
build/rectangle_compression.frd.stdout
build/rectangle_compression.vtu
build/rectangle.py
build/stress_strain_comparison.pdf
build/post_processing.py
build/rectangle_gmsh.inp.stdout
build/SConscript
build/rectangle_compression.frd
build/strip_heading.py
build/vtu2xarray.py
build/stress_strain_comparison.csv
build/rectangle_mesh.inp.stdout
Contents