Introduction
automesh is an automatic mesh generation tool that converts a segmentation
into a hexahedral finite element mesh.
Segmentation
Segmentation is the process of categorizing pixels that compose a digital image into a class that represents some subject of interest. For example, in the image below, the image pixels are classified into classes of sky, trees, cat, grass, and cow.
Figure: Example of semantic segmentation, from Li et al.1
- Semantic segmentation does not differentiate between objects of the same class.
- Instance segmentation does differentiate between objects of the same class.
These two concepts are shown below:
Figure: Distinction between semantic segmentation and instance segmentation, from Lin et al.2
Both segmentation types, semantic and instance, can be used with automesh. However, automesh operates on a 3D segmentation, not a 2D segmentation, as present in a digital image. To obtain a 3D segmentation, two or more images are stacked to compose a volume.
The structured volume of a stacked of pixel composes a volumetric unit called a voxel. A voxel, in the context of this work, will have the same dimensionality in the x and y dimension as the pixel in the image space, and will have the z dimensionality that is the stack interval distance between each image slice. All pixels are rectangular, and all voxels are cuboid.
The figure below illustrates the concept of stacked images:
Figure: Example of stacking several images to create a 3D representation, from Bit et al.3
The digital image sources are frequently medical images, obtained by CT or MR, though automesh can be used for any subject that can be represented as a stacked segmentation. Anatomical regions are classified into categories. For example, in the image below, ten unique integers have been used to represent
bone, disc, vasculature, airway/sinus, membrane, cerebral spinal fluid, white matter, gray matter, muscle, and skin.
Figure: Example of a 3D voxel model, segmented into 10 categories, from Terpsma et al.4
Given a 3D segmentation, for any image slice that composes it, the pixels have been classified
into categories that are designated with unique, non-negative integers.
The range of integer values is limited to 256 = 2^8, since the uint8 data type is specified.
A practical example of a range could be [0, 1, 2, 3, 4]. The integers do not need to be sequential,
so a range of [4, 501, 2, 0, 42] is also valid, but not conventional.
Segmentations are frequently serialized (saved to disc) as either a Numpy (.npy) file
or a SPN (.spn) file.
A SPN file is a text (human-readable) file that contains a single a column of non-negative integer values. Each integer value defines a unique category of a segmentation.
Axis order (for example,
x, y, then z; or, z, y, x, etc.) is not implied by the SPN structure;
so additional data, typically provided through a configuration file, is
needed to uniquely interpret the pixel tile and voxel stack order
of the data in the SPN file.
For subjects that the human anatomy, we use the Patient Coordinate System (PCS), which directs the
x, y, and z axes to the left, posterior, and superior, as shown below:
| Patient Coordinate System: | Left, Posterior, Superior (x, y, z) |
|---|---|
 |  |
Figure: Illustration of the patient coordinate system, left figure from Terpsma et al.4 and right figure from Sharma.5
Finite Element Mesh
ABAQUS
To come.
EXODUS II
EXODUS II is a model developed to store and retrieve data for finite element analyses. It is used for preprocesing (problem definition), postprocessing (results visualization), as well as code to code data transfer. An EXODUS II data file is a random access, machine independent binary file.6
EXODUS II depends on the Network Common Data Form (NetCDF) library.
NetCDF is a public domain database library that provides low-level data storage. The NetCDF library stores data in eXternal Data Representation (XDR) format, which provides machine independency.
EXODUS II library functions provide a map between finite element data objects and NetCDF dimensions, attributes, and variables.
EXODUS II data objects:
- Initialization Data
- Number of nodes
- Number of elements
- optional informational text
- et cetera
- Model - static objects (i.e., objects that do not change over time)
- Nodal coordinates
- Element connectivity
- Node sets
- Side sets
- optional Results
- Nodal results
- Element results
- Global results
Note:
automeshwill use Initialization Data and Model sections; it will not use the Results section.
We use the Exodus II convention for a hexahedral element local node numbering:
Figure: Exodus II hexahedral local finite element numbering scheme, from Schoof et al.6
References
1
Li FF, Johnson J, Yeung S. Lecture 11: Dection and Segmentation, CS 231n, Stanford Unveristy, 2017. link
2
Lin TY, Maire M, Belongie S, Hays J, Perona P, Ramanan D, Dollár P, Zitnick CL. Microsoft coco: Common objects in context. In Computer Vision–ECCV 2014: 13th European Conference, Zurich, Switzerland, September 6-12, 2014, Proceedings, Part V 13 2014 (pp. 740-755). Springer International Publishing. link
3
Bit A, Ghagare D, Rizvanov AA, Chattopadhyay H. Assessment of influences of stenoses in right carotid artery on left carotid artery using wall stress marker. BioMed research international. 2017;2017(1):2935195. link
4
Terpsma RJ, Hovey CB. Blunt impact brain injury using cellular injury criterion. Sandia National Lab. (SNL-NM), Albuquerque, NM (United States); 2020 Oct 1. link
5
Sharma S. DICOM Coordinate Systems — 3D DICOM for computer vision engineers, Medium, 2021-12-22. link
6
Schoof LA, Yarberry VR. EXODUS II: a finite element data model. Sandia National Lab. (SNL-NM), Albuquerque, NM (United States); 1994 Sep 1. link
Installation
The command line interface (CLI) and Python interface are installed separately.
Command Line Interface
cargo install automesh
Python Interface
pip install automesh
Examples
Following are examples created with automesh.
Tests
The following illustates a subset of the unit tests used to validate the code implementation. For a complete listing of the unit sets, see voxels.rs and voxel.py.
The Python script used to generate the figures is included below.
Remark: We use the convention np when importing numpy as follows:
import numpy as np
Single
The minimum working example (MWE) is a single voxel, used to create a single mesh consisting of one block consisting of a single element. The NumPy input single.npy contains the following segmentation:
segmentation = np.array(
[
[
[ 11, ],
],
],
dtype=np.uint8,
)
where the segmentation 11 denotes block 11 in the finite element mesh.
Remark: Serialization (write and read)
| Write | Read |
|---|---|
Use the np.save command to serialize the segmentation a .npy file | Use the np.load command to deserialize the segmentation from a .npy file |
Example: Write the data in segmentation to a file called seg.npy | |
np.save("seg.npy", segmentation) | Example: Read the data from the file seg.npy to a variable called loaded_array |
loaded_array = np.load("seg.npy") |
Equivalently, the single.spn contains a single integer:
11 # x:1 y:1 z:1
The resulting finite element mesh is visualized is shown in the following figure:
Figure: The single.png visualization, (left) lattice node numbers, (right)
mesh node numbers. Lattice node numbers appear in gray, with (x, y, z)
indices in parenthesis. The right-hand rule is used. Lattice coordinates
start at (0, 0, 0), and proceed along the x-axis, then
the y-axis, and then the z-axis.
The finite element mesh local node numbering map to the following global node numbers identically, and :
[1, 2, 4, 3, 5, 6, 8, 7]
->
[1, 2, 4, 3, 5, 6, 8, 7]
which is a special case not typically observed, as shown in more complex examples below.
Double
The next level of complexity example is a two-voxel domain, used to create
a single block composed of two finite elements. We test propagation in
both the x and y directions. The figures below show these two
meshes.
Double X
11 # x:1 y:1 z:1
11 # 2 1 1
where the segmentation 11 denotes block 11 in the finite element mesh.
Figure: Mesh composed of a single block with two elements, propagating along
the x-axis, (left) lattice node numbers, (right) mesh node numbers.
Double Y
11 # x:1 y:1 z:1
11 # 1 2 1
where the segmentation 11 denotes block 11 in the finite element mesh.
Figure: Mesh composed of a single block with two elements, propagating along
the y-axis, (left) lattice node numbers, (right) mesh node numbers.
Triple
11 # x:1 y:1 z:1
11 # 2 1 1
11 # 3 1 1
where the segmentation 11 denotes block 11 in the finite element mesh.
Figure: Mesh composed of a single block with three elements, propagating along
the x-axis, (left) lattice node numbers, (right) mesh node numbers.
Quadruple
11 # x:1 y:1 z:1
11 # 2 1 1
11 # 3 1 1
11 # 4 1 1
where the segmentation 11 denotes block 11 in the finite element mesh.
Figure: Mesh composed of a single block with four elements, propagating along
the x-axis, (left) lattice node numbers, (right) mesh node numbers.
Quadruple with Voids
99 # x:1 y:1 z:1
0 # 2 1 1
0 # 3 1 1
99 # 4 1 1
where the segmentation 99 denotes block 99 in the finite element mesh, and segmentation 0 is excluded from the mesh.
Figure: Mesh composed of a single block with two elements, propagating along
the x-axis and two voids, (left) lattice node numbers, (right) mesh node
numbers.
Quadruple with Two Blocks
100 # x:1 y:1 z:1
101 # 2 1 1
101 # 3 1 1
100 # 4 1 1
where the segmentation 100 and 101 denotes block 100 and 101,
respectively in the finite element mesh.
Figure: Mesh composed of two blocks with two elements elements each,
propagating along the x-axis, (left) lattice node numbers, (right) mesh
node numbers.
Quadruple with Two Blocks and Void
102 # x:1 y:1 z:1
103 # 2 1 1
0 # 3 1 1
102 # 4 1 1
where the segmentation 102 and 103 denotes block 102 and 103,
respectively, in the finite element mesh, and segmentation 0 will be included from the finite element mesh.
Figure: Mesh composed of one block with two elements, a second block with one
element, and a void, propagating along the x-axis, (left) lattice node
numbers, (right) mesh node numbers.
Cube
11 # x:1 y:1 z:1
11 # _ 2 _ 1 1
11 # 1 2 1
11 # _ 2 _ 2 _ 1
11 # 1 1 2
11 # _ 2 _ 1 2
11 # 1 2 2
11 # _ 2 _ 2 _ 2
where the segmentation 11 denotes block 11 in the finite element mesh.
Figure: Mesh composed of one block with eight elements, (left) lattice node numbers, (right) mesh node numbers.
Cube with Multi Blocks and Void
82 # x:1 y:1 z:1
2 # _ 2 _ 1 1
2 # 1 2 1
2 # _ 2 _ 2 _ 1
0 # 1 1 2
31 # _ 2 _ 1 2
0 # 1 2 2
44 # _ 2 _ 2 _ 2
where the segmentation 82, 2, 31 and 44 denotes block 82, 2, 31
and 44, respectively, in the finite element mesh, and segmentation 0 will
be included from the finite element mesh.
Figure: Mesh composed of four blocks (block 82 has one element, block 2
has three elements, block 31 has one element, and block 44 has one
element), (left) lattice node numbers, (right) mesh node numbers.
Cube with Inclusion
11 # x:1 y:1 z:1
11 # 2 1 1
11 # _ 3 _ 1 1
11 # 1 2 1
11 # 2 2 1
11 # _ 3 _ 2 1
11 # 1 3 1
11 # 2 3 1
11 # _ 3 _ 3 _ 1
11 # 1 1 2
11 # 2 1 2
11 # _ 3 _ 1 2
11 # 1 2 2
88 # 2 2 2
11 # _ 3 _ 2 2
11 # 1 3 2
11 # 2 3 2
11 # _ 3 _ 3 _ 2
11 # 1 1 3
11 # 2 1 3
11 # _ 3 _ 1 3
11 # 1 2 3
11 # 2 2 3
11 # _ 3 _ 2 3
11 # 1 3 3
11 # 2 3 3
11 # _ 3 _ 3 _ 3
Figure: Mesh composed of 26 voxels of (block 11) and one voxel inslusion
(block 88), (left) lattice node numbers, (right) mesh node numbers.
Letter F
11 # x:1 y:1 z:1
0 # 2 1 1
0 # _ 3 _ 1 1
11 # 1 2 1
0 # 2 2 1
0 # _ 3 _ 2 1
11 # 1 3 1
11 # 2 3 1
0 # _ 3 _ 3 1
11 # 1 4 1
0 # 2 4 1
0 # _ 3 _ 4 1
11 # 1 5 1
11 # 2 5 1
11 # _ 3 _ 5 _ 1
where the segmentation 11 denotes block 11 in the finite element mesh.
Figure: Mesh composed of a single block with eight elements, (left) lattice node numbers, (right) mesh node numbers.
Letter F in 3D
1 # x:1 y:1 z:1
1 # 2 1 1
1 # 3 1 1
1 # _ 4 _ 1 1
1 # 1 2 1
1 # 2 2 1
1 # 3 2 1
1 # _ 4 _ 2 1
1 # 1 3 1
1 # 2 3 1
1 # 3 3 1
1 # _ 4 _ 3 1
1 # 1 4 1
1 # 2 4 1
1 # 3 4 1
1 # _ 4 _ 4 1
1 # 1 5 1
1 # 2 5 1
1 # 3 5 1
1 # _ 4 _ 5 _ 1
1 # x:1 y:1 z:2
0 # 2 1 2
0 # 3 1 2
0 # _ 4 _ 1 2
1 # 1 2 2
0 # 2 2 2
0 # 3 2 2
0 # _ 4 _ 2 2
1 # 1 3 2
1 # 2 3 2
1 # 3 3 2
1 # _ 4 _ 3 2
1 # 1 4 2
0 # 2 4 2
0 # 3 4 2
0 # _ 4 _ 4 2
1 # 1 5 2
1 # 2 5 2
1 # 3 5 2
1 # _ 4 _ 5 _ 2
1 # x:1 y:1 z:3
0 # 2 1 j
0 # 3 1 2
0 # _ 4 _ 1 2
1 # 1 2 3
0 # 2 2 3
0 # 3 2 3
0 # _ 4 _ 2 3
1 # 1 3 3
0 # 2 3 3
0 # 3 3 3
0 # _ 4 _ 3 3
1 # 1 4 3
0 # 2 4 3
0 # 3 4 3
0 # _ 4 _ 4 3
1 # 1 5 3
1 # 2 5 3
1 # 3 5 3
1 # _ 4 _ 5 _ 3
which corresponds to --nelx 4, --nely 5, and --nelz 3 in the
command line interface.
Figure: Mesh composed of a single block with thirty-nine elements, (left) lattice node numbers, (right) mesh node numbers.
The shape of the solid segmentation is more easily seen without the lattice and element nodes, and with decreased opacity, as shown below:
Figure: Mesh composed of a single block with thirty-nine elements, shown with decreased opacity and without lattice and element node numbers.
Sparse
0 # x:1 y:1 z:1
0 # 2 1 1
0 # 3 1 1
0 # 4 1 1
2 # _ 5 _ 1 1
0 # 1 2 1
1 # 2 2 1
0 # 3 2 1
0 # 4 2 1
2 # _ 5 _ 2 1
1 # 1 3 1
2 # 2 3 1
0 # 3 3 1
2 # 4 3 1
0 # _ 5 _ 3 1
0 # 1 4 1
1 # 2 4 1
0 # 3 4 1
2 # 4 4 1
0 # _ 5 _ 4 1
1 # 1 5 1
0 # 2 5 1
0 # 3 5 1
0 # 4 5 1
1 # _ 5 _ 5 _ 1
2 # x:1 y:1 z:2
0 # 2 1 2
2 # 3 1 2
0 # 4 1 2
0 # _ 5 _ 1 2
1 # 1 2 2
1 # 2 2 2
0 # 3 2 2
2 # 4 2 2
2 # _ 5 _ 2 2
2 # 1 3 2
0 # 2 3 2
0 # 3 3 2
0 # 4 3 2
0 # _ 5 _ 3 2
1 # 1 4 2
0 # 2 4 2
0 # 3 4 2
2 # 4 4 2
0 # _ 5 _ 4 2
2 # 1 5 2
0 # 2 5 2
2 # 3 5 2
0 # 4 5 2
2 # _ 5 _ 5 _ 2
0 # x:1 y:1 z:3
0 # 2 1 3
1 # 3 1 3
0 # 4 1 3
2 # _ 5 _ 1 3
0 # 1 2 3
0 # 2 2 3
0 # 3 2 3
1 # 4 2 3
2 # _ 5 _ 2 3
0 # 1 3 3
0 # 2 3 3
2 # 3 3 3
2 # 4 3 3
2 # _ 5 _ 3 3
0 # 1 4 3
0 # 2 4 3
1 # 3 4 3
0 # 4 4 3
1 # _ 5 _ 4 3
0 # 1 5 3
1 # 2 5 3
0 # 3 5 3
1 # 4 5 3
0 # _ 5 _ 5 _ 3
0 # x:1 y:1 z:4
1 # 2 1 4
2 # 3 1 4
1 # 4 1 4
2 # _ 5 _ 1 4
2 # 1 2 4
0 # 2 2 4
2 # 3 2 4
0 # 4 2 4
1 # _ 5 _ 2 4
1 # 1 3 4
2 # 2 3 4
2 # 3 3 4
0 # 4 3 4
0 # _ 5 _ 3 4
2 # 1 4 4
1 # 2 4 4
1 # 3 4 4
1 # 4 4 4
1 # _ 5 _ 4 4
0 # 1 5 4
0 # 2 5 4
1 # 3 5 4
0 # 4 5 4
0 # _ 5 _ 5 _ 4
0 # x:1 y:1 z:5
1 # 2 1 5
0 # 3 1 5
2 # 4 1 5
0 # _ 5 _ 1 5
1 # 1 2 5
0 # 2 2 5
0 # 3 2 5
0 # 4 2 5
2 # _ 5 _ 2 5
0 # 1 3 5
1 # 2 3 5
0 # 3 3 5
0 # 4 3 5
0 # _ 5 _ 3 5
1 # 1 4 5
0 # 2 4 5
0 # 3 4 5
0 # 4 4 5
0 # _ 5 _ 4 5
0 # 1 5 5
0 # 2 5 5
1 # 3 5 5
2 # 4 5 5
1 # _ 5 _ 5 _ 5
where the segmentation 1 denotes block 1 and segmentation 2 denotes block 2 in the finite eelement mesh (with segmentation 0 excluded).
Figure: Sparse mesh composed of two materials at random voxel locations.
Figure: Sparse mesh composed of two materials at random voxel locations, shown with decreased opactity and without lattice and element node numbers.
Source
The figures were created with figures.py and tested with test_figures.py:
"""This module demonstrates creating a pixel slice in the (x, y)
plane, and then appending layers in the z axis, to create a 3D
voxel lattice, as a precursor for a hexahedral finite element mesh.
This module assumes the virtual environment has been loaded.
Example:
cd ~/autotwin/automesh
source .venv/bin/activate
cd book/examples/unit_tests
pip install matplotlib
python figures.py
Ouput:
The `output_npy` file data structure
The `output_png` file visualization
"""
# standard library
import datetime
from pathlib import Path
from typing import Final, NamedTuple
# third-party libary
import matplotlib.pyplot as plt
from matplotlib.colors import LightSource
import numpy as np
from numpy.typing import NDArray
# module library
# import sandbox.figures_data as fd # why doesn't this work?
def hello_world() -> str:
"""Simple example of a function hooked to a command line entry point.
Returns:
The canonical "Hello world!" string.
"""
return "Hello world!"
class Example(NamedTuple):
"""A base class that has all of the fields required to specialize into a
specific example."""
figure_title: str = "Figure Title"
file_stem: str = "filename"
segmentation = np.array(
[
[
[
1,
],
],
],
dtype=np.uint8,
)
included_ids = (1,)
gold_lattice = None
gold_mesh_lattice_connectivity = None
gold_mesh_element_connectivity = None
COMMON_TITLE: Final[str] = "Lattice Index and Coordinates: "
class Single(Example):
"""A specific example of a single voxel."""
figure_title: str = COMMON_TITLE + "Single"
file_stem: str = "single"
segmentation = np.array(
[
[
[
11,
],
],
],
dtype=np.uint8,
)
included_ids = (11,)
gold_lattice = ((1, 2, 4, 3, 5, 6, 8, 7),)
gold_mesh_lattice_connectivity = (
(
11,
(1, 2, 4, 3, 5, 6, 8, 7),
),
)
gold_mesh_element_connectivity = (
(
11,
(1, 2, 4, 3, 5, 6, 8, 7),
),
)
class DoubleX(Example):
"""A specific example of a double voxel, coursed along the x-axis."""
figure_title: str = COMMON_TITLE + "DoubleX"
file_stem: str = "double_x"
segmentation = np.array(
[
[
[
11,
11,
],
],
],
dtype=np.uint8,
)
included_ids = (11,)
gold_lattice = (
(1, 2, 5, 4, 7, 8, 11, 10),
(2, 3, 6, 5, 8, 9, 12, 11),
)
gold_mesh_lattice_connectivity = (
(
11,
(1, 2, 5, 4, 7, 8, 11, 10),
(2, 3, 6, 5, 8, 9, 12, 11),
),
)
gold_mesh_element_connectivity = (
(
11,
(1, 2, 5, 4, 7, 8, 11, 10),
(2, 3, 6, 5, 8, 9, 12, 11),
),
)
class DoubleY(Example):
"""A specific example of a double voxel, coursed along the y-axis."""
figure_title: str = COMMON_TITLE + "DoubleY"
file_stem: str = "double_y"
segmentation = np.array(
[
[
[
11,
],
[
11,
],
],
],
dtype=np.uint8,
)
included_ids = (11,)
gold_lattice = (
(1, 2, 4, 3, 7, 8, 10, 9),
(3, 4, 6, 5, 9, 10, 12, 11),
)
gold_mesh_lattice_connectivity = (
(
11,
(1, 2, 4, 3, 7, 8, 10, 9),
(3, 4, 6, 5, 9, 10, 12, 11),
),
)
gold_mesh_element_connectivity = (
(
11,
(1, 2, 4, 3, 7, 8, 10, 9),
(3, 4, 6, 5, 9, 10, 12, 11),
),
)
class TripleX(Example):
"""A triple voxel lattice, coursed along the x-axis."""
figure_title: str = COMMON_TITLE + "Triple"
file_stem: str = "triple_x"
segmentation = np.array(
[
[
[
11,
11,
11,
],
],
],
dtype=np.uint8,
)
included_ids = (11,)
gold_lattice = (
(1, 2, 6, 5, 9, 10, 14, 13),
(2, 3, 7, 6, 10, 11, 15, 14),
(3, 4, 8, 7, 11, 12, 16, 15),
)
gold_mesh_lattice_connectivity = (
(
11,
(1, 2, 6, 5, 9, 10, 14, 13),
(2, 3, 7, 6, 10, 11, 15, 14),
(3, 4, 8, 7, 11, 12, 16, 15),
),
)
gold_mesh_element_connectivity = (
(
11,
(1, 2, 6, 5, 9, 10, 14, 13),
(2, 3, 7, 6, 10, 11, 15, 14),
(3, 4, 8, 7, 11, 12, 16, 15),
),
)
class QuadrupleX(Example):
"""A quadruple voxel lattice, coursed along the x-axis."""
figure_title: str = COMMON_TITLE + "Quadruple"
file_stem: str = "quadruple_x"
segmentation = np.array(
[
[
[
11,
11,
11,
11,
],
],
],
dtype=np.uint8,
)
included_ids = (11,)
gold_lattice = (
(1, 2, 7, 6, 11, 12, 17, 16),
(2, 3, 8, 7, 12, 13, 18, 17),
(3, 4, 9, 8, 13, 14, 19, 18),
(4, 5, 10, 9, 14, 15, 20, 19),
)
gold_mesh_lattice_connectivity = (
(
11,
(1, 2, 7, 6, 11, 12, 17, 16),
(2, 3, 8, 7, 12, 13, 18, 17),
(3, 4, 9, 8, 13, 14, 19, 18),
(4, 5, 10, 9, 14, 15, 20, 19),
),
)
gold_mesh_element_connectivity = (
(
11,
(1, 2, 7, 6, 11, 12, 17, 16),
(2, 3, 8, 7, 12, 13, 18, 17),
(3, 4, 9, 8, 13, 14, 19, 18),
(4, 5, 10, 9, 14, 15, 20, 19),
),
)
class Quadruple2VoidsX(Example):
"""A quadruple voxel lattice, coursed along the x-axis, with two
intermediate voxels in the segmentation being void.
"""
figure_title: str = COMMON_TITLE + "Quadruple2VoidsX"
file_stem: str = "quadruple_2_voids_x"
segmentation = np.array(
[
[
[
99,
0,
0,
99,
],
],
],
dtype=np.uint8,
)
included_ids = (99,)
gold_lattice = (
(1, 2, 7, 6, 11, 12, 17, 16),
(2, 3, 8, 7, 12, 13, 18, 17),
(3, 4, 9, 8, 13, 14, 19, 18),
(4, 5, 10, 9, 14, 15, 20, 19),
)
gold_mesh_lattice_connectivity = (
(
99,
(1, 2, 7, 6, 11, 12, 17, 16),
(4, 5, 10, 9, 14, 15, 20, 19),
),
)
gold_mesh_element_connectivity = (
(
99,
(1, 2, 6, 5, 9, 10, 14, 13),
(3, 4, 8, 7, 11, 12, 16, 15),
),
)
class Quadruple2Blocks(Example):
"""A quadruple voxel lattice, with the first intermediate voxel being
the second block and the second intermediate voxel being void.
"""
figure_title: str = COMMON_TITLE + "Quadruple2Blocks"
file_stem: str = "quadruple_2_blocks"
segmentation = np.array(
[
[
[
100,
101,
101,
100,
],
],
],
dtype=np.uint8,
)
included_ids = (
100,
101,
)
gold_lattice = (
(1, 2, 7, 6, 11, 12, 17, 16),
(2, 3, 8, 7, 12, 13, 18, 17),
(3, 4, 9, 8, 13, 14, 19, 18),
(4, 5, 10, 9, 14, 15, 20, 19),
)
gold_mesh_lattice_connectivity = (
(
100,
(1, 2, 7, 6, 11, 12, 17, 16),
(4, 5, 10, 9, 14, 15, 20, 19),
),
(
101,
(2, 3, 8, 7, 12, 13, 18, 17),
(3, 4, 9, 8, 13, 14, 19, 18),
),
)
gold_mesh_element_connectivity = (
(
100,
(1, 2, 7, 6, 11, 12, 17, 16),
(4, 5, 10, 9, 14, 15, 20, 19),
),
(
101,
(2, 3, 8, 7, 12, 13, 18, 17),
(3, 4, 9, 8, 13, 14, 19, 18),
),
)
class Quadruple2BlocksVoid(Example):
"""A quadruple voxel lattice, with the first intermediate voxel being
the second block and the second intermediate voxel being void.
"""
figure_title: str = COMMON_TITLE + "Quadruple2BlocksVoid"
file_stem: str = "quadruple_2_blocks_void"
segmentation = np.array(
[
[
[
102,
103,
0,
102,
],
],
],
dtype=np.uint8,
)
included_ids = (
102,
103,
)
gold_lattice = (
(1, 2, 7, 6, 11, 12, 17, 16),
(2, 3, 8, 7, 12, 13, 18, 17),
(3, 4, 9, 8, 13, 14, 19, 18),
(4, 5, 10, 9, 14, 15, 20, 19),
)
gold_mesh_lattice_connectivity = (
(
102,
(1, 2, 7, 6, 11, 12, 17, 16),
(4, 5, 10, 9, 14, 15, 20, 19),
),
(
103,
(2, 3, 8, 7, 12, 13, 18, 17),
),
)
gold_mesh_element_connectivity = (
(
102,
(1, 2, 7, 6, 11, 12, 17, 16),
(4, 5, 10, 9, 14, 15, 20, 19),
),
(
103,
(2, 3, 8, 7, 12, 13, 18, 17),
),
)
class Cube(Example):
"""A (2 x 2 x 2) voxel cube."""
figure_title: str = COMMON_TITLE + "Cube"
file_stem: str = "cube"
segmentation = np.array(
[
[
[
11,
11,
],
[
11,
11,
],
],
[
[
11,
11,
],
[
11,
11,
],
],
],
dtype=np.uint8,
)
included_ids = (11,)
gold_lattice = (
(1, 2, 5, 4, 10, 11, 14, 13),
(2, 3, 6, 5, 11, 12, 15, 14),
(4, 5, 8, 7, 13, 14, 17, 16),
(5, 6, 9, 8, 14, 15, 18, 17),
(10, 11, 14, 13, 19, 20, 23, 22),
(11, 12, 15, 14, 20, 21, 24, 23),
(13, 14, 17, 16, 22, 23, 26, 25),
(14, 15, 18, 17, 23, 24, 27, 26),
)
gold_mesh_lattice_connectivity = (
(
11,
(1, 2, 5, 4, 10, 11, 14, 13),
(2, 3, 6, 5, 11, 12, 15, 14),
(4, 5, 8, 7, 13, 14, 17, 16),
(5, 6, 9, 8, 14, 15, 18, 17),
(10, 11, 14, 13, 19, 20, 23, 22),
(11, 12, 15, 14, 20, 21, 24, 23),
(13, 14, 17, 16, 22, 23, 26, 25),
(14, 15, 18, 17, 23, 24, 27, 26),
),
)
gold_mesh_element_connectivity = (
(
11,
(1, 2, 5, 4, 10, 11, 14, 13),
(2, 3, 6, 5, 11, 12, 15, 14),
(4, 5, 8, 7, 13, 14, 17, 16),
(5, 6, 9, 8, 14, 15, 18, 17),
(10, 11, 14, 13, 19, 20, 23, 22),
(11, 12, 15, 14, 20, 21, 24, 23),
(13, 14, 17, 16, 22, 23, 26, 25),
(14, 15, 18, 17, 23, 24, 27, 26),
),
)
class CubeMulti(Example):
"""A (2 x 2 x 2) voxel cube with two voids and six elements."""
figure_title: str = COMMON_TITLE + "CubeMulti"
file_stem: str = "cube_multi"
segmentation = np.array(
[
[
[
82,
2,
],
[
2,
2,
],
],
[
[
0,
31,
],
[
0,
44,
],
],
],
dtype=np.uint8,
)
included_ids = (
82,
2,
31,
44,
)
gold_lattice = (
(1, 2, 5, 4, 10, 11, 14, 13),
(2, 3, 6, 5, 11, 12, 15, 14),
(4, 5, 8, 7, 13, 14, 17, 16),
(5, 6, 9, 8, 14, 15, 18, 17),
(10, 11, 14, 13, 19, 20, 23, 22),
(11, 12, 15, 14, 20, 21, 24, 23),
(13, 14, 17, 16, 22, 23, 26, 25),
(14, 15, 18, 17, 23, 24, 27, 26),
)
gold_mesh_lattice_connectivity = (
# (
# 0,
# (10, 11, 14, 13, 19, 20, 23, 22),
# ),
# (
# 0,
# (13, 14, 17, 16, 22, 23, 26, 25),
(
2,
(2, 3, 6, 5, 11, 12, 15, 14),
(4, 5, 8, 7, 13, 14, 17, 16),
(5, 6, 9, 8, 14, 15, 18, 17),
),
(
31,
(11, 12, 15, 14, 20, 21, 24, 23),
),
(
44,
(14, 15, 18, 17, 23, 24, 27, 26),
),
(
82,
(1, 2, 5, 4, 10, 11, 14, 13),
),
)
gold_mesh_element_connectivity = (
(
2,
(2, 3, 6, 5, 11, 12, 15, 14),
(4, 5, 8, 7, 13, 14, 17, 16),
(5, 6, 9, 8, 14, 15, 18, 17),
),
(
31,
(11, 12, 15, 14, 19, 20, 22, 21),
),
(
44,
(14, 15, 18, 17, 21, 22, 24, 23),
),
(
82,
(1, 2, 5, 4, 10, 11, 14, 13),
),
)
class CubeWithInclusion(Example):
"""A (3 x 3 x 3) voxel cube with a single voxel inclusion
at the center.
"""
figure_title: str = COMMON_TITLE + "CubeWithInclusion"
file_stem: str = "cube_with_inclusion"
segmentation = np.array(
[
[
[
11,
11,
11,
],
[
11,
11,
11,
],
[
11,
11,
11,
],
],
[
[
11,
11,
11,
],
[
11,
88,
11,
],
[
11,
11,
11,
],
],
[
[
11,
11,
11,
],
[
11,
11,
11,
],
[
11,
11,
11,
],
],
],
dtype=np.uint8,
)
included_ids = (
11,
88,
)
gold_lattice = (
(1, 2, 6, 5, 17, 18, 22, 21),
(2, 3, 7, 6, 18, 19, 23, 22),
(3, 4, 8, 7, 19, 20, 24, 23),
(5, 6, 10, 9, 21, 22, 26, 25),
(6, 7, 11, 10, 22, 23, 27, 26),
(7, 8, 12, 11, 23, 24, 28, 27),
(9, 10, 14, 13, 25, 26, 30, 29),
(10, 11, 15, 14, 26, 27, 31, 30),
(11, 12, 16, 15, 27, 28, 32, 31),
(17, 18, 22, 21, 33, 34, 38, 37),
(18, 19, 23, 22, 34, 35, 39, 38),
(19, 20, 24, 23, 35, 36, 40, 39),
(21, 22, 26, 25, 37, 38, 42, 41),
(22, 23, 27, 26, 38, 39, 43, 42),
(23, 24, 28, 27, 39, 40, 44, 43),
(25, 26, 30, 29, 41, 42, 46, 45),
(26, 27, 31, 30, 42, 43, 47, 46),
(27, 28, 32, 31, 43, 44, 48, 47),
(33, 34, 38, 37, 49, 50, 54, 53),
(34, 35, 39, 38, 50, 51, 55, 54),
(35, 36, 40, 39, 51, 52, 56, 55),
(37, 38, 42, 41, 53, 54, 58, 57),
(38, 39, 43, 42, 54, 55, 59, 58),
(39, 40, 44, 43, 55, 56, 60, 59),
(41, 42, 46, 45, 57, 58, 62, 61),
(42, 43, 47, 46, 58, 59, 63, 62),
(43, 44, 48, 47, 59, 60, 64, 63),
)
gold_mesh_lattice_connectivity = (
(
11,
(1, 2, 6, 5, 17, 18, 22, 21),
(2, 3, 7, 6, 18, 19, 23, 22),
(3, 4, 8, 7, 19, 20, 24, 23),
(5, 6, 10, 9, 21, 22, 26, 25),
(6, 7, 11, 10, 22, 23, 27, 26),
(7, 8, 12, 11, 23, 24, 28, 27),
(9, 10, 14, 13, 25, 26, 30, 29),
(10, 11, 15, 14, 26, 27, 31, 30),
(11, 12, 16, 15, 27, 28, 32, 31),
(17, 18, 22, 21, 33, 34, 38, 37),
(18, 19, 23, 22, 34, 35, 39, 38),
(19, 20, 24, 23, 35, 36, 40, 39),
(21, 22, 26, 25, 37, 38, 42, 41),
(23, 24, 28, 27, 39, 40, 44, 43),
(25, 26, 30, 29, 41, 42, 46, 45),
(26, 27, 31, 30, 42, 43, 47, 46),
(27, 28, 32, 31, 43, 44, 48, 47),
(33, 34, 38, 37, 49, 50, 54, 53),
(34, 35, 39, 38, 50, 51, 55, 54),
(35, 36, 40, 39, 51, 52, 56, 55),
(37, 38, 42, 41, 53, 54, 58, 57),
(38, 39, 43, 42, 54, 55, 59, 58),
(39, 40, 44, 43, 55, 56, 60, 59),
(41, 42, 46, 45, 57, 58, 62, 61),
(42, 43, 47, 46, 58, 59, 63, 62),
(43, 44, 48, 47, 59, 60, 64, 63),
),
(
88,
(22, 23, 27, 26, 38, 39, 43, 42),
),
)
gold_mesh_element_connectivity = (
(
11,
(1, 2, 6, 5, 17, 18, 22, 21),
(2, 3, 7, 6, 18, 19, 23, 22),
(3, 4, 8, 7, 19, 20, 24, 23),
(5, 6, 10, 9, 21, 22, 26, 25),
(6, 7, 11, 10, 22, 23, 27, 26),
(7, 8, 12, 11, 23, 24, 28, 27),
(9, 10, 14, 13, 25, 26, 30, 29),
(10, 11, 15, 14, 26, 27, 31, 30),
(11, 12, 16, 15, 27, 28, 32, 31),
(17, 18, 22, 21, 33, 34, 38, 37),
(18, 19, 23, 22, 34, 35, 39, 38),
(19, 20, 24, 23, 35, 36, 40, 39),
(21, 22, 26, 25, 37, 38, 42, 41),
(23, 24, 28, 27, 39, 40, 44, 43),
(25, 26, 30, 29, 41, 42, 46, 45),
(26, 27, 31, 30, 42, 43, 47, 46),
(27, 28, 32, 31, 43, 44, 48, 47),
(33, 34, 38, 37, 49, 50, 54, 53),
(34, 35, 39, 38, 50, 51, 55, 54),
(35, 36, 40, 39, 51, 52, 56, 55),
(37, 38, 42, 41, 53, 54, 58, 57),
(38, 39, 43, 42, 54, 55, 59, 58),
(39, 40, 44, 43, 55, 56, 60, 59),
(41, 42, 46, 45, 57, 58, 62, 61),
(42, 43, 47, 46, 58, 59, 63, 62),
(43, 44, 48, 47, 59, 60, 64, 63),
),
(
88,
(22, 23, 27, 26, 38, 39, 43, 42),
),
)
class LetterF(Example):
"""A minimal letter F example."""
figure_title: str = COMMON_TITLE + "LetterF"
file_stem: str = "letter_f"
segmentation = np.array(
[
[
[
11,
0,
0,
],
[
11,
0,
0,
],
[
11,
11,
0,
],
[
11,
0,
0,
],
[
11,
11,
11,
],
],
],
dtype=np.uint8,
)
included_ids = (11,)
gold_lattice = (
(1, 2, 6, 5, 25, 26, 30, 29),
(2, 3, 7, 6, 26, 27, 31, 30),
(3, 4, 8, 7, 27, 28, 32, 31),
(5, 6, 10, 9, 29, 30, 34, 33),
(6, 7, 11, 10, 30, 31, 35, 34),
(7, 8, 12, 11, 31, 32, 36, 35),
(9, 10, 14, 13, 33, 34, 38, 37),
(10, 11, 15, 14, 34, 35, 39, 38),
(11, 12, 16, 15, 35, 36, 40, 39),
(13, 14, 18, 17, 37, 38, 42, 41),
(14, 15, 19, 18, 38, 39, 43, 42),
(15, 16, 20, 19, 39, 40, 44, 43),
(17, 18, 22, 21, 41, 42, 46, 45),
(18, 19, 23, 22, 42, 43, 47, 46),
(19, 20, 24, 23, 43, 44, 48, 47),
)
gold_mesh_lattice_connectivity = (
(
11,
(1, 2, 6, 5, 25, 26, 30, 29),
# (2, 3, 7, 6, 26, 27, 31, 30),
# (3, 4, 8, 7, 27, 28, 32, 31),
(5, 6, 10, 9, 29, 30, 34, 33),
# (6, 7, 11, 10, 30, 31, 35, 34),
# (7, 8, 12, 11, 31, 32, 36, 35),
(9, 10, 14, 13, 33, 34, 38, 37),
(10, 11, 15, 14, 34, 35, 39, 38),
# (11, 12, 16, 15, 35, 36, 40, 39),
(13, 14, 18, 17, 37, 38, 42, 41),
# (14, 15, 19, 18, 38, 39, 43, 42),
# (15, 16, 20, 19, 39, 40, 44, 43),
(17, 18, 22, 21, 41, 42, 46, 45),
(18, 19, 23, 22, 42, 43, 47, 46),
(19, 20, 24, 23, 43, 44, 48, 47),
),
)
gold_mesh_element_connectivity = (
(
11,
(1, 2, 4, 3, 19, 20, 22, 21),
#
#
(3, 4, 6, 5, 21, 22, 24, 23),
#
#
(5, 6, 9, 8, 23, 24, 27, 26),
(6, 7, 10, 9, 24, 25, 28, 27),
#
(8, 9, 12, 11, 26, 27, 30, 29),
#
#
(11, 12, 16, 15, 29, 30, 34, 33),
(12, 13, 17, 16, 30, 31, 35, 34),
(13, 14, 18, 17, 31, 32, 36, 35),
),
)
class LetterF3D(Example):
"""A three dimensional variation of the letter F, in a non-standard
orientation.
"""
figure_title: str = COMMON_TITLE + "LetterF3D"
file_stem: str = "letter_f_3d"
segmentation = np.array(
[
[
[1, 1, 1, 1],
[1, 1, 1, 1],
[1, 1, 1, 1],
[1, 1, 1, 1],
[1, 1, 1, 1],
],
[
[1, 0, 0, 0],
[1, 0, 0, 0],
[1, 1, 1, 1],
[1, 0, 0, 0],
[1, 1, 1, 1],
],
[
[1, 0, 0, 0],
[1, 0, 0, 0],
[1, 0, 0, 0],
[1, 0, 0, 0],
[1, 1, 1, 1],
],
],
dtype=np.uint8,
)
included_ids = (1,)
gold_lattice = (
(1, 2, 7, 6, 31, 32, 37, 36),
(2, 3, 8, 7, 32, 33, 38, 37),
(3, 4, 9, 8, 33, 34, 39, 38),
(4, 5, 10, 9, 34, 35, 40, 39),
(6, 7, 12, 11, 36, 37, 42, 41),
(7, 8, 13, 12, 37, 38, 43, 42),
(8, 9, 14, 13, 38, 39, 44, 43),
(9, 10, 15, 14, 39, 40, 45, 44),
(11, 12, 17, 16, 41, 42, 47, 46),
(12, 13, 18, 17, 42, 43, 48, 47),
(13, 14, 19, 18, 43, 44, 49, 48),
(14, 15, 20, 19, 44, 45, 50, 49),
(16, 17, 22, 21, 46, 47, 52, 51),
(17, 18, 23, 22, 47, 48, 53, 52),
(18, 19, 24, 23, 48, 49, 54, 53),
(19, 20, 25, 24, 49, 50, 55, 54),
(21, 22, 27, 26, 51, 52, 57, 56),
(22, 23, 28, 27, 52, 53, 58, 57),
(23, 24, 29, 28, 53, 54, 59, 58),
(24, 25, 30, 29, 54, 55, 60, 59),
(31, 32, 37, 36, 61, 62, 67, 66),
(32, 33, 38, 37, 62, 63, 68, 67),
(33, 34, 39, 38, 63, 64, 69, 68),
(34, 35, 40, 39, 64, 65, 70, 69),
(36, 37, 42, 41, 66, 67, 72, 71),
(37, 38, 43, 42, 67, 68, 73, 72),
(38, 39, 44, 43, 68, 69, 74, 73),
(39, 40, 45, 44, 69, 70, 75, 74),
(41, 42, 47, 46, 71, 72, 77, 76),
(42, 43, 48, 47, 72, 73, 78, 77),
(43, 44, 49, 48, 73, 74, 79, 78),
(44, 45, 50, 49, 74, 75, 80, 79),
(46, 47, 52, 51, 76, 77, 82, 81),
(47, 48, 53, 52, 77, 78, 83, 82),
(48, 49, 54, 53, 78, 79, 84, 83),
(49, 50, 55, 54, 79, 80, 85, 84),
(51, 52, 57, 56, 81, 82, 87, 86),
(52, 53, 58, 57, 82, 83, 88, 87),
(53, 54, 59, 58, 83, 84, 89, 88),
(54, 55, 60, 59, 84, 85, 90, 89),
(61, 62, 67, 66, 91, 92, 97, 96),
(62, 63, 68, 67, 92, 93, 98, 97),
(63, 64, 69, 68, 93, 94, 99, 98),
(64, 65, 70, 69, 94, 95, 100, 99),
(66, 67, 72, 71, 96, 97, 102, 101),
(67, 68, 73, 72, 97, 98, 103, 102),
(68, 69, 74, 73, 98, 99, 104, 103),
(69, 70, 75, 74, 99, 100, 105, 104),
(71, 72, 77, 76, 101, 102, 107, 106),
(72, 73, 78, 77, 102, 103, 108, 107),
(73, 74, 79, 78, 103, 104, 109, 108),
(74, 75, 80, 79, 104, 105, 110, 109),
(76, 77, 82, 81, 106, 107, 112, 111),
(77, 78, 83, 82, 107, 108, 113, 112),
(78, 79, 84, 83, 108, 109, 114, 113),
(79, 80, 85, 84, 109, 110, 115, 114),
(81, 82, 87, 86, 111, 112, 117, 116),
(82, 83, 88, 87, 112, 113, 118, 117),
(83, 84, 89, 88, 113, 114, 119, 118),
(84, 85, 90, 89, 114, 115, 120, 119),
)
gold_mesh_lattice_connectivity = (
(
1,
(1, 2, 7, 6, 31, 32, 37, 36),
(2, 3, 8, 7, 32, 33, 38, 37),
(3, 4, 9, 8, 33, 34, 39, 38),
(4, 5, 10, 9, 34, 35, 40, 39),
(6, 7, 12, 11, 36, 37, 42, 41),
(7, 8, 13, 12, 37, 38, 43, 42),
(8, 9, 14, 13, 38, 39, 44, 43),
(9, 10, 15, 14, 39, 40, 45, 44),
(11, 12, 17, 16, 41, 42, 47, 46),
(12, 13, 18, 17, 42, 43, 48, 47),
(13, 14, 19, 18, 43, 44, 49, 48),
(14, 15, 20, 19, 44, 45, 50, 49),
(16, 17, 22, 21, 46, 47, 52, 51),
(17, 18, 23, 22, 47, 48, 53, 52),
(18, 19, 24, 23, 48, 49, 54, 53),
(19, 20, 25, 24, 49, 50, 55, 54),
(21, 22, 27, 26, 51, 52, 57, 56),
(22, 23, 28, 27, 52, 53, 58, 57),
(23, 24, 29, 28, 53, 54, 59, 58),
(24, 25, 30, 29, 54, 55, 60, 59),
(31, 32, 37, 36, 61, 62, 67, 66),
(36, 37, 42, 41, 66, 67, 72, 71),
(41, 42, 47, 46, 71, 72, 77, 76),
(42, 43, 48, 47, 72, 73, 78, 77),
(43, 44, 49, 48, 73, 74, 79, 78),
(44, 45, 50, 49, 74, 75, 80, 79),
(46, 47, 52, 51, 76, 77, 82, 81),
(51, 52, 57, 56, 81, 82, 87, 86),
(52, 53, 58, 57, 82, 83, 88, 87),
(53, 54, 59, 58, 83, 84, 89, 88),
(54, 55, 60, 59, 84, 85, 90, 89),
(61, 62, 67, 66, 91, 92, 97, 96),
(66, 67, 72, 71, 96, 97, 102, 101),
(71, 72, 77, 76, 101, 102, 107, 106),
(76, 77, 82, 81, 106, 107, 112, 111),
(81, 82, 87, 86, 111, 112, 117, 116),
(82, 83, 88, 87, 112, 113, 118, 117),
(83, 84, 89, 88, 113, 114, 119, 118),
(84, 85, 90, 89, 114, 115, 120, 119),
),
)
gold_mesh_element_connectivity = (
(
1,
(1, 2, 7, 6, 31, 32, 37, 36),
(2, 3, 8, 7, 32, 33, 38, 37),
(3, 4, 9, 8, 33, 34, 39, 38),
(4, 5, 10, 9, 34, 35, 40, 39),
(6, 7, 12, 11, 36, 37, 42, 41),
(7, 8, 13, 12, 37, 38, 43, 42),
(8, 9, 14, 13, 38, 39, 44, 43),
(9, 10, 15, 14, 39, 40, 45, 44),
(11, 12, 17, 16, 41, 42, 47, 46),
(12, 13, 18, 17, 42, 43, 48, 47),
(13, 14, 19, 18, 43, 44, 49, 48),
(14, 15, 20, 19, 44, 45, 50, 49),
(16, 17, 22, 21, 46, 47, 52, 51),
(17, 18, 23, 22, 47, 48, 53, 52),
(18, 19, 24, 23, 48, 49, 54, 53),
(19, 20, 25, 24, 49, 50, 55, 54),
(21, 22, 27, 26, 51, 52, 57, 56),
(22, 23, 28, 27, 52, 53, 58, 57),
(23, 24, 29, 28, 53, 54, 59, 58),
(24, 25, 30, 29, 54, 55, 60, 59),
(31, 32, 37, 36, 61, 62, 64, 63),
(36, 37, 42, 41, 63, 64, 66, 65),
(41, 42, 47, 46, 65, 66, 71, 70),
(42, 43, 48, 47, 66, 67, 72, 71),
(43, 44, 49, 48, 67, 68, 73, 72),
(44, 45, 50, 49, 68, 69, 74, 73),
(46, 47, 52, 51, 70, 71, 76, 75),
(51, 52, 57, 56, 75, 76, 81, 80),
(52, 53, 58, 57, 76, 77, 82, 81),
(53, 54, 59, 58, 77, 78, 83, 82),
(54, 55, 60, 59, 78, 79, 84, 83),
(61, 62, 64, 63, 85, 86, 88, 87),
(63, 64, 66, 65, 87, 88, 90, 89),
(65, 66, 71, 70, 89, 90, 92, 91),
(70, 71, 76, 75, 91, 92, 94, 93),
(75, 76, 81, 80, 93, 94, 99, 98),
(76, 77, 82, 81, 94, 95, 100, 99),
(77, 78, 83, 82, 95, 96, 101, 100),
(78, 79, 84, 83, 96, 97, 102, 101),
),
)
class Sparse(Example):
"""A radomized 5x5x5 segmentation."""
figure_title: str = COMMON_TITLE + "Sparse"
file_stem: str = "sparse"
segmentation = np.array(
[
[
[0, 0, 0, 0, 2],
[0, 1, 0, 0, 2],
[1, 2, 0, 2, 0],
[0, 1, 0, 2, 0],
[1, 0, 0, 0, 1],
],
[
[2, 0, 2, 0, 0],
[1, 1, 0, 2, 2],
[2, 0, 0, 0, 0],
[1, 0, 0, 2, 0],
[2, 0, 2, 0, 2],
],
[
[0, 0, 1, 0, 2],
[0, 0, 0, 1, 2],
[0, 0, 2, 2, 2],
[0, 0, 1, 0, 1],
[0, 1, 0, 1, 0],
],
[
[0, 1, 2, 1, 2],
[2, 0, 2, 0, 1],
[1, 2, 2, 0, 0],
[2, 1, 1, 1, 1],
[0, 0, 1, 0, 0],
],
[
[0, 1, 0, 2, 0],
[1, 0, 0, 0, 2],
[0, 1, 0, 0, 0],
[1, 0, 0, 0, 0],
[0, 0, 1, 2, 1],
],
],
dtype=np.uint8,
)
included_ids = (
1,
2,
)
gold_lattice = (
(1, 2, 8, 7, 37, 38, 44, 43),
(2, 3, 9, 8, 38, 39, 45, 44),
(3, 4, 10, 9, 39, 40, 46, 45),
(4, 5, 11, 10, 40, 41, 47, 46),
(5, 6, 12, 11, 41, 42, 48, 47),
(7, 8, 14, 13, 43, 44, 50, 49),
(8, 9, 15, 14, 44, 45, 51, 50),
(9, 10, 16, 15, 45, 46, 52, 51),
(10, 11, 17, 16, 46, 47, 53, 52),
(11, 12, 18, 17, 47, 48, 54, 53),
(13, 14, 20, 19, 49, 50, 56, 55),
(14, 15, 21, 20, 50, 51, 57, 56),
(15, 16, 22, 21, 51, 52, 58, 57),
(16, 17, 23, 22, 52, 53, 59, 58),
(17, 18, 24, 23, 53, 54, 60, 59),
(19, 20, 26, 25, 55, 56, 62, 61),
(20, 21, 27, 26, 56, 57, 63, 62),
(21, 22, 28, 27, 57, 58, 64, 63),
(22, 23, 29, 28, 58, 59, 65, 64),
(23, 24, 30, 29, 59, 60, 66, 65),
(25, 26, 32, 31, 61, 62, 68, 67),
(26, 27, 33, 32, 62, 63, 69, 68),
(27, 28, 34, 33, 63, 64, 70, 69),
(28, 29, 35, 34, 64, 65, 71, 70),
(29, 30, 36, 35, 65, 66, 72, 71),
(37, 38, 44, 43, 73, 74, 80, 79),
(38, 39, 45, 44, 74, 75, 81, 80),
(39, 40, 46, 45, 75, 76, 82, 81),
(40, 41, 47, 46, 76, 77, 83, 82),
(41, 42, 48, 47, 77, 78, 84, 83),
(43, 44, 50, 49, 79, 80, 86, 85),
(44, 45, 51, 50, 80, 81, 87, 86),
(45, 46, 52, 51, 81, 82, 88, 87),
(46, 47, 53, 52, 82, 83, 89, 88),
(47, 48, 54, 53, 83, 84, 90, 89),
(49, 50, 56, 55, 85, 86, 92, 91),
(50, 51, 57, 56, 86, 87, 93, 92),
(51, 52, 58, 57, 87, 88, 94, 93),
(52, 53, 59, 58, 88, 89, 95, 94),
(53, 54, 60, 59, 89, 90, 96, 95),
(55, 56, 62, 61, 91, 92, 98, 97),
(56, 57, 63, 62, 92, 93, 99, 98),
(57, 58, 64, 63, 93, 94, 100, 99),
(58, 59, 65, 64, 94, 95, 101, 100),
(59, 60, 66, 65, 95, 96, 102, 101),
(61, 62, 68, 67, 97, 98, 104, 103),
(62, 63, 69, 68, 98, 99, 105, 104),
(63, 64, 70, 69, 99, 100, 106, 105),
(64, 65, 71, 70, 100, 101, 107, 106),
(65, 66, 72, 71, 101, 102, 108, 107),
(73, 74, 80, 79, 109, 110, 116, 115),
(74, 75, 81, 80, 110, 111, 117, 116),
(75, 76, 82, 81, 111, 112, 118, 117),
(76, 77, 83, 82, 112, 113, 119, 118),
(77, 78, 84, 83, 113, 114, 120, 119),
(79, 80, 86, 85, 115, 116, 122, 121),
(80, 81, 87, 86, 116, 117, 123, 122),
(81, 82, 88, 87, 117, 118, 124, 123),
(82, 83, 89, 88, 118, 119, 125, 124),
(83, 84, 90, 89, 119, 120, 126, 125),
(85, 86, 92, 91, 121, 122, 128, 127),
(86, 87, 93, 92, 122, 123, 129, 128),
(87, 88, 94, 93, 123, 124, 130, 129),
(88, 89, 95, 94, 124, 125, 131, 130),
(89, 90, 96, 95, 125, 126, 132, 131),
(91, 92, 98, 97, 127, 128, 134, 133),
(92, 93, 99, 98, 128, 129, 135, 134),
(93, 94, 100, 99, 129, 130, 136, 135),
(94, 95, 101, 100, 130, 131, 137, 136),
(95, 96, 102, 101, 131, 132, 138, 137),
(97, 98, 104, 103, 133, 134, 140, 139),
(98, 99, 105, 104, 134, 135, 141, 140),
(99, 100, 106, 105, 135, 136, 142, 141),
(100, 101, 107, 106, 136, 137, 143, 142),
(101, 102, 108, 107, 137, 138, 144, 143),
(109, 110, 116, 115, 145, 146, 152, 151),
(110, 111, 117, 116, 146, 147, 153, 152),
(111, 112, 118, 117, 147, 148, 154, 153),
(112, 113, 119, 118, 148, 149, 155, 154),
(113, 114, 120, 119, 149, 150, 156, 155),
(115, 116, 122, 121, 151, 152, 158, 157),
(116, 117, 123, 122, 152, 153, 159, 158),
(117, 118, 124, 123, 153, 154, 160, 159),
(118, 119, 125, 124, 154, 155, 161, 160),
(119, 120, 126, 125, 155, 156, 162, 161),
(121, 122, 128, 127, 157, 158, 164, 163),
(122, 123, 129, 128, 158, 159, 165, 164),
(123, 124, 130, 129, 159, 160, 166, 165),
(124, 125, 131, 130, 160, 161, 167, 166),
(125, 126, 132, 131, 161, 162, 168, 167),
(127, 128, 134, 133, 163, 164, 170, 169),
(128, 129, 135, 134, 164, 165, 171, 170),
(129, 130, 136, 135, 165, 166, 172, 171),
(130, 131, 137, 136, 166, 167, 173, 172),
(131, 132, 138, 137, 167, 168, 174, 173),
(133, 134, 140, 139, 169, 170, 176, 175),
(134, 135, 141, 140, 170, 171, 177, 176),
(135, 136, 142, 141, 171, 172, 178, 177),
(136, 137, 143, 142, 172, 173, 179, 178),
(137, 138, 144, 143, 173, 174, 180, 179),
(145, 146, 152, 151, 181, 182, 188, 187),
(146, 147, 153, 152, 182, 183, 189, 188),
(147, 148, 154, 153, 183, 184, 190, 189),
(148, 149, 155, 154, 184, 185, 191, 190),
(149, 150, 156, 155, 185, 186, 192, 191),
(151, 152, 158, 157, 187, 188, 194, 193),
(152, 153, 159, 158, 188, 189, 195, 194),
(153, 154, 160, 159, 189, 190, 196, 195),
(154, 155, 161, 160, 190, 191, 197, 196),
(155, 156, 162, 161, 191, 192, 198, 197),
(157, 158, 164, 163, 193, 194, 200, 199),
(158, 159, 165, 164, 194, 195, 201, 200),
(159, 160, 166, 165, 195, 196, 202, 201),
(160, 161, 167, 166, 196, 197, 203, 202),
(161, 162, 168, 167, 197, 198, 204, 203),
(163, 164, 170, 169, 199, 200, 206, 205),
(164, 165, 171, 170, 200, 201, 207, 206),
(165, 166, 172, 171, 201, 202, 208, 207),
(166, 167, 173, 172, 202, 203, 209, 208),
(167, 168, 174, 173, 203, 204, 210, 209),
(169, 170, 176, 175, 205, 206, 212, 211),
(170, 171, 177, 176, 206, 207, 213, 212),
(171, 172, 178, 177, 207, 208, 214, 213),
(172, 173, 179, 178, 208, 209, 215, 214),
(173, 174, 180, 179, 209, 210, 216, 215),
)
gold_mesh_lattice_connectivity = (
(
1,
(8, 9, 15, 14, 44, 45, 51, 50),
(13, 14, 20, 19, 49, 50, 56, 55),
(20, 21, 27, 26, 56, 57, 63, 62),
(25, 26, 32, 31, 61, 62, 68, 67),
(29, 30, 36, 35, 65, 66, 72, 71),
(43, 44, 50, 49, 79, 80, 86, 85),
(44, 45, 51, 50, 80, 81, 87, 86),
(55, 56, 62, 61, 91, 92, 98, 97),
(75, 76, 82, 81, 111, 112, 118, 117),
(82, 83, 89, 88, 118, 119, 125, 124),
(93, 94, 100, 99, 129, 130, 136, 135),
(95, 96, 102, 101, 131, 132, 138, 137),
(98, 99, 105, 104, 134, 135, 141, 140),
(100, 101, 107, 106, 136, 137, 143, 142),
(110, 111, 117, 116, 146, 147, 153, 152),
(112, 113, 119, 118, 148, 149, 155, 154),
(119, 120, 126, 125, 155, 156, 162, 161),
(121, 122, 128, 127, 157, 158, 164, 163),
(128, 129, 135, 134, 164, 165, 171, 170),
(129, 130, 136, 135, 165, 166, 172, 171),
(130, 131, 137, 136, 166, 167, 173, 172),
(131, 132, 138, 137, 167, 168, 174, 173),
(135, 136, 142, 141, 171, 172, 178, 177),
(146, 147, 153, 152, 182, 183, 189, 188),
(151, 152, 158, 157, 187, 188, 194, 193),
(158, 159, 165, 164, 194, 195, 201, 200),
(163, 164, 170, 169, 199, 200, 206, 205),
(171, 172, 178, 177, 207, 208, 214, 213),
(173, 174, 180, 179, 209, 210, 216, 215),
),
(
2,
(5, 6, 12, 11, 41, 42, 48, 47),
(11, 12, 18, 17, 47, 48, 54, 53),
(14, 15, 21, 20, 50, 51, 57, 56),
(16, 17, 23, 22, 52, 53, 59, 58),
(22, 23, 29, 28, 58, 59, 65, 64),
(37, 38, 44, 43, 73, 74, 80, 79),
(39, 40, 46, 45, 75, 76, 82, 81),
(46, 47, 53, 52, 82, 83, 89, 88),
(47, 48, 54, 53, 83, 84, 90, 89),
(49, 50, 56, 55, 85, 86, 92, 91),
(58, 59, 65, 64, 94, 95, 101, 100),
(61, 62, 68, 67, 97, 98, 104, 103),
(63, 64, 70, 69, 99, 100, 106, 105),
(65, 66, 72, 71, 101, 102, 108, 107),
(77, 78, 84, 83, 113, 114, 120, 119),
(83, 84, 90, 89, 119, 120, 126, 125),
(87, 88, 94, 93, 123, 124, 130, 129),
(88, 89, 95, 94, 124, 125, 131, 130),
(89, 90, 96, 95, 125, 126, 132, 131),
(111, 112, 118, 117, 147, 148, 154, 153),
(113, 114, 120, 119, 149, 150, 156, 155),
(115, 116, 122, 121, 151, 152, 158, 157),
(117, 118, 124, 123, 153, 154, 160, 159),
(122, 123, 129, 128, 158, 159, 165, 164),
(123, 124, 130, 129, 159, 160, 166, 165),
(127, 128, 134, 133, 163, 164, 170, 169),
(148, 149, 155, 154, 184, 185, 191, 190),
(155, 156, 162, 161, 191, 192, 198, 197),
(172, 173, 179, 178, 208, 209, 215, 214),
),
)
gold_mesh_element_connectivity = (
(
1,
(3, 4, 9, 8, 35, 36, 42, 41),
(7, 8, 14, 13, 40, 41, 47, 46),
(14, 15, 20, 19, 47, 48, 53, 52),
(18, 19, 25, 24, 51, 52, 58, 57),
(22, 23, 27, 26, 55, 56, 62, 61),
(34, 35, 41, 40, 69, 70, 76, 75),
(35, 36, 42, 41, 70, 71, 77, 76),
(46, 47, 52, 51, 81, 82, 88, 87),
(65, 66, 72, 71, 100, 101, 107, 106),
(72, 73, 79, 78, 107, 108, 114, 113),
(83, 84, 90, 89, 118, 119, 125, 124),
(85, 86, 92, 91, 120, 121, 127, 126),
(88, 89, 95, 94, 123, 124, 129, 128),
(90, 91, 97, 96, 125, 126, 131, 130),
(99, 100, 106, 105, 132, 133, 139, 138),
(101, 102, 108, 107, 134, 135, 141, 140),
(108, 109, 115, 114, 141, 142, 148, 147),
(110, 111, 117, 116, 143, 144, 150, 149),
(117, 118, 124, 123, 150, 151, 157, 156),
(118, 119, 125, 124, 151, 152, 158, 157),
(119, 120, 126, 125, 152, 153, 159, 158),
(120, 121, 127, 126, 153, 154, 160, 159),
(124, 125, 130, 129, 157, 158, 162, 161),
(132, 133, 139, 138, 165, 166, 171, 170),
(137, 138, 144, 143, 169, 170, 176, 175),
(144, 145, 151, 150, 176, 177, 182, 181),
(149, 150, 156, 155, 180, 181, 184, 183),
(157, 158, 162, 161, 185, 186, 190, 189),
(159, 160, 164, 163, 187, 188, 192, 191),
),
(
2,
(1, 2, 6, 5, 32, 33, 39, 38),
(5, 6, 12, 11, 38, 39, 45, 44),
(8, 9, 15, 14, 41, 42, 48, 47),
(10, 11, 17, 16, 43, 44, 50, 49),
(16, 17, 22, 21, 49, 50, 55, 54),
(28, 29, 35, 34, 63, 64, 70, 69),
(30, 31, 37, 36, 65, 66, 72, 71),
(37, 38, 44, 43, 72, 73, 79, 78),
(38, 39, 45, 44, 73, 74, 80, 79),
(40, 41, 47, 46, 75, 76, 82, 81),
(49, 50, 55, 54, 84, 85, 91, 90),
(51, 52, 58, 57, 87, 88, 94, 93),
(53, 54, 60, 59, 89, 90, 96, 95),
(55, 56, 62, 61, 91, 92, 98, 97),
(67, 68, 74, 73, 102, 103, 109, 108),
(73, 74, 80, 79, 108, 109, 115, 114),
(77, 78, 84, 83, 112, 113, 119, 118),
(78, 79, 85, 84, 113, 114, 120, 119),
(79, 80, 86, 85, 114, 115, 121, 120),
(100, 101, 107, 106, 133, 134, 140, 139),
(102, 103, 109, 108, 135, 136, 142, 141),
(104, 105, 111, 110, 137, 138, 144, 143),
(106, 107, 113, 112, 139, 140, 146, 145),
(111, 112, 118, 117, 144, 145, 151, 150),
(112, 113, 119, 118, 145, 146, 152, 151),
(116, 117, 123, 122, 149, 150, 156, 155),
(134, 135, 141, 140, 167, 168, 173, 172),
(141, 142, 148, 147, 173, 174, 179, 178),
(158, 159, 163, 162, 186, 187, 191, 190),
),
)
def lattice_connectivity(ex: Example) -> NDArray[np.uint8]:
"""Given an Example, prints the lattice connectivity."""
offset = 0
nz, ny, nx = ex.segmentation.shape
nzp, nyp, nxp = nz + 1, ny + 1, nx + 1
# Generate the lattice nodes
lattice_nodes = []
lattice_node = 0
for k in range(nzp):
for j in range(nyp):
for i in range(nxp):
lattice_node += 1
lattice_nodes.append([lattice_node, i, j, k])
# connectivity for each voxel
cvs = []
offset = 0
# print("processing indices...")
for iz in range(nz):
for iy in range(ny):
for ix in range(nx):
# print(f"(ix, iy, iz) = ({ix}, {iy}, {iz})")
cv = offset + np.array(
[
(iz + 0) * (nxp * nyp) + (iy + 0) * nxp + ix + 1,
(iz + 0) * (nxp * nyp) + (iy + 0) * nxp + ix + 2,
(iz + 0) * (nxp * nyp) + (iy + 1) * nxp + ix + 2,
(iz + 0) * (nxp * nyp) + (iy + 1) * nxp + ix + 1,
(iz + 1) * (nxp * nyp) + (iy + 0) * nxp + ix + 1,
(iz + 1) * (nxp * nyp) + (iy + 0) * nxp + ix + 2,
(iz + 1) * (nxp * nyp) + (iy + 1) * nxp + ix + 2,
(iz + 1) * (nxp * nyp) + (iy + 1) * nxp + ix + 1,
]
)
cvs.append(cv)
cs = np.vstack(cvs)
# voxel by voxel comparison
# breakpoint()
vv = ex.gold_lattice == cs
assert np.all(vv)
return cs
def mesh_lattice_connectivity(
ex: Example,
lattice: np.ndarray,
) -> tuple:
"""Given an Example (in particular, the Example's voxel data structure,
a segmentation) and the `lattice_connectivity`, create the connectivity
for the mesh with lattice node numbers. A voxel with a segmentation id not
in the Example's included ids tuple is excluded from the mesh.
"""
# segmentation = ex.segmentation.flatten().squeeze()
segmentation = ex.segmentation.flatten()
# breakpoint()
# assert that the list of included ids is equal
included_set_unordered = set(ex.included_ids)
included_list_ordered = sorted(included_set_unordered)
# breakpoint()
seg_set = set(segmentation)
for item in included_list_ordered:
assert (
item in seg_set
), f"Error: `included_ids` item {item} is not in the segmentation"
# Create a list of finite elements from the lattice elements. If the
# lattice element has a segmentation id that is not in the included_ids,
# exlude the voxel element from the collected list to create the finite
# element list
blocks = () # empty tuple
# breakpoint()
for bb in included_list_ordered:
# included_elements = []
elements = () # empty tuple
elements = elements + (bb,) # insert the block number
for i, element in enumerate(lattice):
if bb == segmentation[i]:
# breakpoint()
elements = elements + (tuple(element.tolist()),) # overwrite
blocks = blocks + (elements,) # overwrite
# breakpoint()
# return np.array(blocks)
return blocks
def renumber(source: tuple, old: tuple, new: tuple) -> tuple:
"""Given a source tuple, composed of a list of positive integers,
a tuple of `old` numbers that maps into `new` numbers, return the
source tuple with the `new` numbers."""
# the old and the new tuples musts have the same length
err = "Tuples `old` and `new` must have equal length."
assert len(old) == len(new), err
result = ()
for item in source:
idx = old.index(item)
new_value = new[idx]
result = result + (new_value,)
return result
def mesh_element_connectivity(mesh_with_lattice_connectivity: tuple):
"""Given a mesh with lattice connectivity, return a mesh with finite
element connectivity.
"""
# create a list of unordered lattice node numbers
ln = []
for item in mesh_with_lattice_connectivity:
# print(f"item is {item}")
# The first item is the block number
# block = item[0]
# The second and onward items are the elements
elements = item[1:]
for element in elements:
ln += list(element)
ln_set = set(ln) # sets are not necessarily ordered
ln_ordered = tuple(sorted(ln_set)) # now these unique integers are ordered
# and they will map into the new compressed unique interger list `mapsto`
mapsto = tuple(range(1, len(ln_ordered) + 1))
# now build a mesh_with_element_connectivity
mesh = () # empty tuple
# breakpoint()
for item in mesh_with_lattice_connectivity:
# The first item is the block number
block_number = item[0]
block_and_elements = () # empty tuple
# insert the block number
block_and_elements = block_and_elements + (block_number,)
for element in item[1:]:
new_element = renumber(source=element, old=ln_ordered, new=mapsto)
# overwrite
block_and_elements = block_and_elements + (new_element,)
mesh = mesh + (block_and_elements,) # overwrite
return mesh
def flatten_tuple(t):
"""Uses recursion to convert nested tuples into a single-sevel tuple.
Example:
nested_tuple = (1, (2, 3), (4, (5, 6)), 7)
flattened_tuple = flatten_tuple(nested_tuple)
print(flattened_tuple) # Output: (1, 2, 3, 4, 5, 6, 7)
"""
flat_list = []
for item in t:
if isinstance(item, tuple):
flat_list.extend(flatten_tuple(item))
else:
flat_list.append(item)
# breakpoint()
return tuple(flat_list)
def elements_without_block_ids(mesh: tuple) -> tuple:
"""Given a mesh, removes the block ids and returns only just the
element connectivities.
"""
aa = ()
for item in mesh:
bb = item[1:]
aa = aa + bb
return aa
def main():
"""The main program."""
# Create an instance of a specific example
# user input begin
examples = [
Single(),
DoubleX(),
# DoubleY(),
# TripleX(),
# QuadrupleX(),
# Quadruple2VoidsX(),
# Quadruple2Blocks(),
# Quadruple2BlocksVoid(),
# Cube(),
# CubeMulti(),
# CubeWithInclusion(),
# LetterF(),
# LetterF3D(),
# Sparse(),
]
# output_dir: Final[str] = "~/scratch"
output_dir: Final[Path] = Path(__file__).parent
dpi: Final[int] = 300 # resolution, dots per inch
for ex in examples:
# computation
output_npy: Path = (
Path(output_dir).expanduser().joinpath(ex.file_stem + ".npy")
)
# visualization
visualize: bool = True # True performs post-processing visualization
output_png_short = ex.file_stem + ".png"
output_png: Path = (
Path(output_dir).expanduser().joinpath(output_png_short)
)
# el, az, roll = 25, -115, 0
# el, az, roll = 28, -115, 0
# el, az, roll = 63, -110, 0 # used for most visuals
el, az, roll = 11, -111, 0 # used for CubeWithInclusion
# el, az, roll = 60, -121, 0
# el, az, roll = 42, -120, 0
#
# colors
# cmap = cm.get_cmap("viridis") # viridis colormap
# cmap = plt.get_cmap(name="viridis")
cmap = plt.get_cmap(name="tab10")
# number of discrete colors
num_colors = len(ex.included_ids)
colors = cmap(np.linspace(0, 1, num_colors))
# breakpoint()
# azimuth (deg):
# 0 is east (from +y-axis looking back toward origin)
# 90 is north (from +x-axis looking back toward origin)
# 180 is west (from -y-axis looking back toward origin)
# 270 is south (from -x-axis looking back toward origin)
# elevation (deg): 0 is horizontal, 90 is vertical (+z-axis up)
lightsource = LightSource(azdeg=325, altdeg=45) # azimuth, elevation
nodes_shown: bool = True
# nodes_shown: bool = False
voxel_alpha: float = 0.1
# voxel_alpha: float = 0.7
# io: if the output directory does not already exist, create it
output_path = Path(output_dir).expanduser()
if not output_path.exists():
print(f"Could not find existing output directory: {output_path}")
Path.mkdir(output_path)
print(f"Created: {output_path}")
assert output_path.exists()
nelz, nely, nelx = ex.segmentation.shape
lc = lattice_connectivity(ex=ex)
mesh_w_lattice_conn = mesh_lattice_connectivity(ex=ex, lattice=lc)
# breakpoint()
err = "Calculated lattice connectivity error."
assert mesh_w_lattice_conn == ex.gold_mesh_lattice_connectivity, err
mesh_w_element_conn = mesh_element_connectivity(mesh_w_lattice_conn)
err = "Calcualted element connectivity error." # overwrite
assert mesh_w_element_conn == ex.gold_mesh_element_connectivity, err
# save the numpy data as a .npy file
np.save(output_npy, ex.segmentation)
print(f"Saved: {output_npy}")
# to load the array back from the .npy file,
# use the numpy.load function:
loaded_array = np.load(output_npy)
# verify the loaded array
print(f"segmentation loaded from saved file: {loaded_array}")
assert np.all(loaded_array == ex.segmentation)
# now that the .npy file has been created and verified,
# move it to the repo at ~/autotwin/automesh/tests/input
if not visualize:
return
# visualization
# Define the dimensions of the lattice
nxp, nyp, nzp = (nelx + 1, nely + 1, nelz + 1)
# Create a figure and a 3D axis
# fig = plt.figure()
fig = plt.figure(figsize=(10, 5)) # Adjust the figure size
# fig = plt.figure(figsize=(8, 4)) # Adjust the figure size
# ax = fig.add_subplot(111, projection="3d")
# figure with 1 row, 2 columns
ax = fig.add_subplot(1, 2, 1, projection="3d") # r1, c2, 1st subplot
ax2 = fig.add_subplot(1, 2, 2, projection="3d") # r1, c2, 2nd subplot
# For 3D plotting of voxels in matplotlib, we must swap the 'x' and the
# 'z' axes. The original axes in the segmentation are (z, y, x) and
# are numbered (0, 1, 2). We want new exists as (x, y, z) and thus
# with numbering (2, 1, 0).
vox = np.transpose(ex.segmentation, (2, 1, 0))
# add voxels for each of the included materials
for i, block_id in enumerate(ex.included_ids):
# breakpoint()
solid = vox == block_id
# ax.voxels(solid, facecolors=voxel_color, alpha=voxel_alpha)
# ax.voxels(solid, facecolors=colors[i], alpha=voxel_alpha)
ax.voxels(
solid,
facecolors=colors[i],
edgecolor=colors[i],
alpha=voxel_alpha,
lightsource=lightsource,
)
# plot the same voxels on the 2nd axis
ax2.voxels(
solid,
facecolors=colors[i],
edgecolor=colors[i],
alpha=voxel_alpha,
lightsource=lightsource,
)
# breakpoint()
# Generate the lattice points
x = []
y = []
z = []
labels = []
# Generate the element points
xel = []
yel = []
zel = []
# generate a set from the element connectivity
# breakpoint()
# ec_set = set(flatten_tuple(mesh_w_lattice_conn)) # bug!
# bug fix:
ec_set = set(
flatten_tuple(elements_without_block_ids(mesh_w_lattice_conn))
)
# breakpoint()
lattice_ijk = 0
# gnn = global node number
gnn = 0
gnn_labels = []
for k in range(nzp):
for j in range(nyp):
for i in range(nxp):
x.append(i)
y.append(j)
z.append(k)
if lattice_ijk + 1 in ec_set:
gnn += 1
xel.append(i)
yel.append(j)
zel.append(k)
gnn_labels.append(f" {gnn}")
lattice_ijk += 1
labels.append(f" {lattice_ijk}: ({i},{j},{k})")
if nodes_shown:
# Plot the lattice coordinates
ax.scatter(
x,
y,
z,
s=20,
facecolors="red",
edgecolors="none",
)
# Label the lattice coordinates
for n, label in enumerate(labels):
ax.text(x[n], y[n], z[n], label, color="darkgray", fontsize=8)
# Plot the nodes included in the finite element connectivity
ax2.scatter(
xel,
yel,
zel,
s=30,
facecolors="blue",
edgecolors="blue",
)
# Label the global node numbers
for n, label in enumerate(gnn_labels):
ax2.text(
xel[n], yel[n], zel[n], label, color="darkblue", fontsize=8
)
# Set labels for the axes
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")
# repeat for the 2nd axis
ax2.set_xlabel("x")
ax2.set_ylabel("y")
ax2.set_zlabel("z")
x_ticks = list(range(nxp))
y_ticks = list(range(nyp))
z_ticks = list(range(nzp))
ax.set_xticks(x_ticks)
ax.set_yticks(y_ticks)
ax.set_zticks(z_ticks)
# repeat for the 2nd axis
ax2.set_xticks(x_ticks)
ax2.set_yticks(y_ticks)
ax2.set_zticks(z_ticks)
ax.set_xlim(float(x_ticks[0]), float(x_ticks[-1]))
ax.set_ylim(float(y_ticks[0]), float(y_ticks[-1]))
ax.set_zlim(float(z_ticks[0]), float(z_ticks[-1]))
# repeat for the 2nd axis
ax2.set_xlim(float(x_ticks[0]), float(x_ticks[-1]))
ax2.set_ylim(float(y_ticks[0]), float(y_ticks[-1]))
ax2.set_zlim(float(z_ticks[0]), float(z_ticks[-1]))
# Set the camera view
ax.set_aspect("equal")
ax.view_init(elev=el, azim=az, roll=roll)
# repeat for the 2nd axis
ax2.set_aspect("equal")
ax2.view_init(elev=el, azim=az, roll=roll)
# Adjust the distance of the camera. The default value is 10.
# Increasing/decreasing this value will zoom in/out, respectively.
# ax.dist = 5 # Change the distance of the camera
# Doesn't seem to work, and the title is clipping the uppermost node
# and lattice numbers, so suppress the titles for now.
# Set the title
# ax.set_title(ex.figure_title)
# Add a footnote
# Get the current date and time in UTC
now_utc = datetime.datetime.now(datetime.UTC)
# Format the date and time as a string
timestamp_utc = now_utc.strftime("%Y-%m-%d %H:%M:%S UTC")
fn = f"Figure: {output_png_short} "
fn += f"created with {__file__}\non {timestamp_utc}."
fig.text(0.5, 0.01, fn, ha="center", fontsize=8)
# Show the plot
plt.show()
# plt.show()
fig.savefig(output_png, dpi=dpi)
print(f"Saved: {output_png}")
if __name__ == "__main__":
main()
r"""This module tests functionality of the included module.
Example:
cd ~/autotwin/automesh
Activate the venv with one of the following:
source .venv/bin/activate # for bash shell
source .venv/bin/activate.csh # for c shell
source .venv/bin/activate.fish # for fish shell
.\.venv\Scripts\activate # for powershell
cd book/examples/unit_tests
python -m pytest test_figures.py -v # -v is for verbose
to run a single test in this module, for example `test_hello` function:
python -m pytest test_figures.py::test_hello -v
"""
import pytest
import figures as ff
def test_hello():
"""This is a minimum working example (MWE) test."""
assert ff.hello_world() == "Hello world!"
def test_renumber():
"""Tests that the renumber function works as expected."""
source = (300, 22, 1)
old = (1, 22, 300, 40)
new = (42, 2, 9, 1000)
result = ff.renumber(source=source, old=old, new=new)
assert result == (9, 2, 42)
# Assure that tuples old and new of unequal length raise an AssertionError
new = (42, 2) # overwrite
err = "Tuples `old` and `new` must have equal length."
with pytest.raises(AssertionError, match=err):
_ = ff.renumber(source=source, old=old, new=new)
def test_mesh_with_element_connectivity():
"""Test CubeMulti by hand."""
gold_mesh_lattice_connectivity = (
(
2,
(2, 3, 6, 5, 11, 12, 15, 14),
(4, 5, 8, 7, 13, 14, 17, 16),
(5, 6, 9, 8, 14, 15, 18, 17),
),
(
31,
(11, 12, 15, 14, 20, 21, 24, 23),
),
(
44,
(14, 15, 18, 17, 23, 24, 27, 26),
),
(
82,
(1, 2, 5, 4, 10, 11, 14, 13),
),
)
gold_mesh_element_connectivity = (
(
2,
(2, 3, 6, 5, 11, 12, 15, 14),
(4, 5, 8, 7, 13, 14, 17, 16),
(5, 6, 9, 8, 14, 15, 18, 17),
),
(31, (11, 12, 15, 14, 19, 20, 22, 21)),
(44, (14, 15, 18, 17, 21, 22, 24, 23)),
(82, (1, 2, 5, 4, 10, 11, 14, 13)),
)
result = ff.mesh_element_connectivity(
mesh_with_lattice_connectivity=gold_mesh_lattice_connectivity
)
assert result == gold_mesh_element_connectivity
def test_elements_no_block_ids():
"""Given a mesh, strips the block ids from the"""
known_input = (
(
2,
(2, 3, 6, 5, 11, 12, 15, 14),
(4, 5, 8, 7, 13, 14, 17, 16),
(5, 6, 9, 8, 14, 15, 18, 17),
),
(31, (11, 12, 15, 14, 20, 21, 24, 23)),
(44, (14, 15, 18, 17, 23, 24, 27, 26)),
(82, (1, 2, 5, 4, 10, 11, 14, 13)),
)
gold_output = (
(2, 3, 6, 5, 11, 12, 15, 14),
(4, 5, 8, 7, 13, 14, 17, 16),
(5, 6, 9, 8, 14, 15, 18, 17),
(11, 12, 15, 14, 20, 21, 24, 23),
(14, 15, 18, 17, 23, 24, 27, 26),
(1, 2, 5, 4, 10, 11, 14, 13),
)
result = ff.elements_without_block_ids(mesh=known_input)
assert result == gold_output
Spheres
We segment a sphere into coarse voxel meshes.
Segmentation
Using spheres.py,
"""This module creates a voxelized sphere and exports it as a .npy file.
Example
-------
source ~/autotwin/automesh/.venv/bin/activate
python spheres.py
"""
from pathlib import Path
from typing import Final
from matplotlib.colors import LightSource
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
def sphere(radius: int, dtype=np.uint8) -> np.ndarray:
"""Generate a 3D voxelized representation of a sphere.
Parameters
----------
radius: int
The radius of the sphere. Minimum value is 1.
dtype: data-type, optional
The data type of the output array. Default is np.uint8.
Returns
-------
np.ndarray
A 3D numpy array of shape (2*radius+1, 2*radius+1, 2*radius+1)
representing the voxelized sphere. Voxels within the sphere are
set to 1, and those outside are set to 0.
Raises
------
ValueError
If the radius is less than 1.
Example
-------
>>> sphere(radius=1) returns
array(
[
[[0, 0, 0], [0, 1, 0], [0, 0, 0]],
[[0, 1, 0], [1, 1, 1], [0, 1, 0]],
[[0, 0, 0], [0, 1, 0], [0, 0, 0]]
],
dtype=uint8
)
Reference
---------
Adapted from:
https://github.com/scikit-image/scikit-image/blob/v0.24.0/skimage/morphology/footprints.py#L763-L833
"""
if radius < 1:
raise ValueError("Radius must be >= 1")
n_voxels_per_side = 2 * radius + 1
vox_z, vox_y, vox_x = np.mgrid[
-radius:radius:n_voxels_per_side * 1j,
-radius:radius:n_voxels_per_side * 1j,
-radius:radius:n_voxels_per_side * 1j,
]
voxel_radius_squared = vox_x**2 + vox_y**2 + vox_z**2
result = np.array(voxel_radius_squared <= radius * radius, dtype=dtype)
return result
# User input begin
spheres = {
"radius_1": sphere(radius=1),
"radius_3": sphere(radius=3),
"radius_5": sphere(radius=5),
}
aa = Path(__file__)
bb = aa.with_suffix(".png")
# Visualize the elements.
# width, height = 8, 4
width, height = 6, 3
fig = plt.figure(figsize=(width, height))
# fig = plt.figure(figsize=(8, 8))
el, az, roll = 63, -110, 0
cmap = plt.get_cmap(name="tab10")
num_colors = len(spheres)
voxel_alpha: Final[float] = 0.9
colors = cmap(np.linspace(0, 1, num_colors))
lightsource = LightSource(azdeg=325, altdeg=45) # azimuth, elevation
# lightsource = LightSource(azdeg=325, altdeg=90) # azimuth, elevation
dpi: Final[int] = 300 # resolution, dots per inch
visualize: Final[bool] = False # turn to True to show the figure on screen
serialize: Final[bool] = False # turn to True to save .png and .npy files
# User input end
N_SUBPLOTS = len(spheres)
IDX = 1
for index, (key, value) in enumerate(spheres.items()):
ax = fig.add_subplot(1, N_SUBPLOTS, index+1, projection=Axes3D.name)
ax.voxels(
value,
facecolors=colors[index],
edgecolor=colors[index],
alpha=voxel_alpha,
lightsource=lightsource)
ax.set_title(key)
IDX += 1
# Set labels for the axes
ax.set_xlabel("x (voxels)")
ax.set_ylabel("y (voxels)")
ax.set_zlabel("z (voxels)")
# Set the camera view
ax.set_aspect("equal")
ax.view_init(elev=el, azim=az, roll=roll)
if serialize:
cc = aa.with_stem("spheres_" + key)
dd = cc.with_suffix(".npy")
# Save the data in .npy format
np.save(dd, value)
print(f"Saved: {dd}")
fig.tight_layout()
if visualize:
plt.show()
if serialize:
fig.savefig(bb, dpi=dpi)
print(f"Saved: {bb}")
create very coarse spheres of varying
resolution (radius=1, radius=3, and radius=5), as shown below:
Figure: Sphere segmentations at selected resolutions, shown in the voxel domain.
For the radius=1 case, the underyling data structure appears as:
spheres["radius_1"]
array([[[0, 0, 0],
[0, 1, 0],
[0, 0, 0]],
[[0, 1, 0],
[1, 1, 1],
[0, 1, 0]],
[[0, 0, 0],
[0, 1, 0],
[0, 0, 0]]], dtype=uint8)
For the radius=3 case, the underyling data structure appears as:
spheres["sradius_3"]
array([[[0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 1, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0]],
[[0, 0, 0, 0, 0, 0, 0],
[0, 0, 1, 1, 1, 0, 0],
[0, 1, 1, 1, 1, 1, 0],
[0, 1, 1, 1, 1, 1, 0],
[0, 1, 1, 1, 1, 1, 0],
[0, 0, 1, 1, 1, 0, 0],
[0, 0, 0, 0, 0, 0, 0]],
[[0, 0, 0, 0, 0, 0, 0],
[0, 1, 1, 1, 1, 1, 0],
[0, 1, 1, 1, 1, 1, 0],
[0, 1, 1, 1, 1, 1, 0],
[0, 1, 1, 1, 1, 1, 0],
[0, 1, 1, 1, 1, 1, 0],
[0, 0, 0, 0, 0, 0, 0]],
[[0, 0, 0, 1, 0, 0, 0],
[0, 1, 1, 1, 1, 1, 0],
[0, 1, 1, 1, 1, 1, 0],
[1, 1, 1, 1, 1, 1, 1],
[0, 1, 1, 1, 1, 1, 0],
[0, 1, 1, 1, 1, 1, 0],
[0, 0, 0, 1, 0, 0, 0]],
[[0, 0, 0, 0, 0, 0, 0],
[0, 1, 1, 1, 1, 1, 0],
[0, 1, 1, 1, 1, 1, 0],
[0, 1, 1, 1, 1, 1, 0],
[0, 1, 1, 1, 1, 1, 0],
[0, 1, 1, 1, 1, 1, 0],
[0, 0, 0, 0, 0, 0, 0]],
[[0, 0, 0, 0, 0, 0, 0],
[0, 0, 1, 1, 1, 0, 0],
[0, 1, 1, 1, 1, 1, 0],
[0, 1, 1, 1, 1, 1, 0],
[0, 1, 1, 1, 1, 1, 0],
[0, 0, 1, 1, 1, 0, 0],
[0, 0, 0, 0, 0, 0, 0]],
[[0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 1, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0]]], dtype=uint8)
Because of its large size, the data structure for sphere_5 is not shown here.
These segmentations are saved to
Autotwin
Autotwin is used to convert the .npy segmentations into .inp meshes.
automesh -i spheres_radius_1.npy -o spheres_radius_1.inp -x 3 -y 3 -z 3
automesh -i spheres_radius_3.npy -o spheres_radius_3.inp -x 7 -y 7 -z 7
automesh -i spheres_radius_5.npy -o spheres_radius_5_.inp -x 11 -y 11 -z 11
Mesh
The spheres_radius_1.inp file:
*HEADING
autotwin.automesh
version 0.1.6
autogenerated on 2024-09-23 20:42:20.599041 UTC
**
*NODE, NSET=ALLNODES
1, 1.000000e0, 1.000000e0, 0.000000e0
2, 2.000000e0, 1.000000e0, 0.000000e0
3, 1.000000e0, 2.000000e0, 0.000000e0
4, 2.000000e0, 2.000000e0, 0.000000e0
5, 1.000000e0, 0.000000e0, 1.000000e0
6, 2.000000e0, 0.000000e0, 1.000000e0
7, 0.000000e0, 1.000000e0, 1.000000e0
8, 1.000000e0, 1.000000e0, 1.000000e0
9, 2.000000e0, 1.000000e0, 1.000000e0
10, 3.000000e0, 1.000000e0, 1.000000e0
11, 0.000000e0, 2.000000e0, 1.000000e0
12, 1.000000e0, 2.000000e0, 1.000000e0
13, 2.000000e0, 2.000000e0, 1.000000e0
14, 3.000000e0, 2.000000e0, 1.000000e0
15, 1.000000e0, 3.000000e0, 1.000000e0
16, 2.000000e0, 3.000000e0, 1.000000e0
17, 1.000000e0, 0.000000e0, 2.000000e0
18, 2.000000e0, 0.000000e0, 2.000000e0
19, 0.000000e0, 1.000000e0, 2.000000e0
20, 1.000000e0, 1.000000e0, 2.000000e0
21, 2.000000e0, 1.000000e0, 2.000000e0
22, 3.000000e0, 1.000000e0, 2.000000e0
23, 0.000000e0, 2.000000e0, 2.000000e0
24, 1.000000e0, 2.000000e0, 2.000000e0
25, 2.000000e0, 2.000000e0, 2.000000e0
26, 3.000000e0, 2.000000e0, 2.000000e0
27, 1.000000e0, 3.000000e0, 2.000000e0
28, 2.000000e0, 3.000000e0, 2.000000e0
29, 1.000000e0, 1.000000e0, 3.000000e0
30, 2.000000e0, 1.000000e0, 3.000000e0
31, 1.000000e0, 2.000000e0, 3.000000e0
32, 2.000000e0, 2.000000e0, 3.000000e0
**
*ELEMENT, TYPE=C3D8R, ELSET=EB1
1, 1, 2, 4, 3, 8, 9, 13, 12
2, 5, 6, 9, 8, 17, 18, 21, 20
3, 7, 8, 12, 11, 19, 20, 24, 23
4, 8, 9, 13, 12, 20, 21, 25, 24
5, 9, 10, 14, 13, 21, 22, 26, 25
6, 12, 13, 16, 15, 24, 25, 28, 27
7, 20, 21, 25, 24, 29, 30, 32, 31
**
*SOLID SECTION, ELSET=EB1, MATERIAL=Default-Steel
The spheres_radius_3.inp file:
*HEADING
autotwin.automesh
version 0.1.6
autogenerated on 2024-09-23 20:44:01.593367 UTC
**
*NODE, NSET=ALLNODES
1, 3.000000e0, 3.000000e0, 0.000000e0
2, 4.000000e0, 3.000000e0, 0.000000e0
3, 3.000000e0, 4.000000e0, 0.000000e0
4, 4.000000e0, 4.000000e0, 0.000000e0
5, 2.000000e0, 1.000000e0, 1.000000e0
6, 3.000000e0, 1.000000e0, 1.000000e0
7, 4.000000e0, 1.000000e0, 1.000000e0
8, 5.000000e0, 1.000000e0, 1.000000e0
9, 1.000000e0, 2.000000e0, 1.000000e0
10, 2.000000e0, 2.000000e0, 1.000000e0
11, 3.000000e0, 2.000000e0, 1.000000e0
12, 4.000000e0, 2.000000e0, 1.000000e0
13, 5.000000e0, 2.000000e0, 1.000000e0
14, 6.000000e0, 2.000000e0, 1.000000e0
15, 1.000000e0, 3.000000e0, 1.000000e0
16, 2.000000e0, 3.000000e0, 1.000000e0
17, 3.000000e0, 3.000000e0, 1.000000e0
18, 4.000000e0, 3.000000e0, 1.000000e0
19, 5.000000e0, 3.000000e0, 1.000000e0
20, 6.000000e0, 3.000000e0, 1.000000e0
21, 1.000000e0, 4.000000e0, 1.000000e0
22, 2.000000e0, 4.000000e0, 1.000000e0
23, 3.000000e0, 4.000000e0, 1.000000e0
24, 4.000000e0, 4.000000e0, 1.000000e0
25, 5.000000e0, 4.000000e0, 1.000000e0
26, 6.000000e0, 4.000000e0, 1.000000e0
27, 1.000000e0, 5.000000e0, 1.000000e0
28, 2.000000e0, 5.000000e0, 1.000000e0
29, 3.000000e0, 5.000000e0, 1.000000e0
30, 4.000000e0, 5.000000e0, 1.000000e0
31, 5.000000e0, 5.000000e0, 1.000000e0
32, 6.000000e0, 5.000000e0, 1.000000e0
33, 2.000000e0, 6.000000e0, 1.000000e0
34, 3.000000e0, 6.000000e0, 1.000000e0
35, 4.000000e0, 6.000000e0, 1.000000e0
36, 5.000000e0, 6.000000e0, 1.000000e0
37, 1.000000e0, 1.000000e0, 2.000000e0
38, 2.000000e0, 1.000000e0, 2.000000e0
39, 3.000000e0, 1.000000e0, 2.000000e0
40, 4.000000e0, 1.000000e0, 2.000000e0
41, 5.000000e0, 1.000000e0, 2.000000e0
42, 6.000000e0, 1.000000e0, 2.000000e0
43, 1.000000e0, 2.000000e0, 2.000000e0
44, 2.000000e0, 2.000000e0, 2.000000e0
45, 3.000000e0, 2.000000e0, 2.000000e0
46, 4.000000e0, 2.000000e0, 2.000000e0
47, 5.000000e0, 2.000000e0, 2.000000e0
48, 6.000000e0, 2.000000e0, 2.000000e0
49, 1.000000e0, 3.000000e0, 2.000000e0
50, 2.000000e0, 3.000000e0, 2.000000e0
51, 3.000000e0, 3.000000e0, 2.000000e0
52, 4.000000e0, 3.000000e0, 2.000000e0
53, 5.000000e0, 3.000000e0, 2.000000e0
54, 6.000000e0, 3.000000e0, 2.000000e0
55, 1.000000e0, 4.000000e0, 2.000000e0
56, 2.000000e0, 4.000000e0, 2.000000e0
57, 3.000000e0, 4.000000e0, 2.000000e0
58, 4.000000e0, 4.000000e0, 2.000000e0
59, 5.000000e0, 4.000000e0, 2.000000e0
60, 6.000000e0, 4.000000e0, 2.000000e0
61, 1.000000e0, 5.000000e0, 2.000000e0
62, 2.000000e0, 5.000000e0, 2.000000e0
63, 3.000000e0, 5.000000e0, 2.000000e0
64, 4.000000e0, 5.000000e0, 2.000000e0
65, 5.000000e0, 5.000000e0, 2.000000e0
66, 6.000000e0, 5.000000e0, 2.000000e0
67, 1.000000e0, 6.000000e0, 2.000000e0
68, 2.000000e0, 6.000000e0, 2.000000e0
69, 3.000000e0, 6.000000e0, 2.000000e0
70, 4.000000e0, 6.000000e0, 2.000000e0
71, 5.000000e0, 6.000000e0, 2.000000e0
72, 6.000000e0, 6.000000e0, 2.000000e0
73, 3.000000e0, 0.000000e0, 3.000000e0
74, 4.000000e0, 0.000000e0, 3.000000e0
75, 1.000000e0, 1.000000e0, 3.000000e0
76, 2.000000e0, 1.000000e0, 3.000000e0
77, 3.000000e0, 1.000000e0, 3.000000e0
78, 4.000000e0, 1.000000e0, 3.000000e0
79, 5.000000e0, 1.000000e0, 3.000000e0
80, 6.000000e0, 1.000000e0, 3.000000e0
81, 1.000000e0, 2.000000e0, 3.000000e0
82, 2.000000e0, 2.000000e0, 3.000000e0
83, 3.000000e0, 2.000000e0, 3.000000e0
84, 4.000000e0, 2.000000e0, 3.000000e0
85, 5.000000e0, 2.000000e0, 3.000000e0
86, 6.000000e0, 2.000000e0, 3.000000e0
87, 0.000000e0, 3.000000e0, 3.000000e0
88, 1.000000e0, 3.000000e0, 3.000000e0
89, 2.000000e0, 3.000000e0, 3.000000e0
90, 3.000000e0, 3.000000e0, 3.000000e0
91, 4.000000e0, 3.000000e0, 3.000000e0
92, 5.000000e0, 3.000000e0, 3.000000e0
93, 6.000000e0, 3.000000e0, 3.000000e0
94, 7.000000e0, 3.000000e0, 3.000000e0
95, 0.000000e0, 4.000000e0, 3.000000e0
96, 1.000000e0, 4.000000e0, 3.000000e0
97, 2.000000e0, 4.000000e0, 3.000000e0
98, 3.000000e0, 4.000000e0, 3.000000e0
99, 4.000000e0, 4.000000e0, 3.000000e0
100, 5.000000e0, 4.000000e0, 3.000000e0
101, 6.000000e0, 4.000000e0, 3.000000e0
102, 7.000000e0, 4.000000e0, 3.000000e0
103, 1.000000e0, 5.000000e0, 3.000000e0
104, 2.000000e0, 5.000000e0, 3.000000e0
105, 3.000000e0, 5.000000e0, 3.000000e0
106, 4.000000e0, 5.000000e0, 3.000000e0
107, 5.000000e0, 5.000000e0, 3.000000e0
108, 6.000000e0, 5.000000e0, 3.000000e0
109, 1.000000e0, 6.000000e0, 3.000000e0
110, 2.000000e0, 6.000000e0, 3.000000e0
111, 3.000000e0, 6.000000e0, 3.000000e0
112, 4.000000e0, 6.000000e0, 3.000000e0
113, 5.000000e0, 6.000000e0, 3.000000e0
114, 6.000000e0, 6.000000e0, 3.000000e0
115, 3.000000e0, 7.000000e0, 3.000000e0
116, 4.000000e0, 7.000000e0, 3.000000e0
117, 3.000000e0, 0.000000e0, 4.000000e0
118, 4.000000e0, 0.000000e0, 4.000000e0
119, 1.000000e0, 1.000000e0, 4.000000e0
120, 2.000000e0, 1.000000e0, 4.000000e0
121, 3.000000e0, 1.000000e0, 4.000000e0
122, 4.000000e0, 1.000000e0, 4.000000e0
123, 5.000000e0, 1.000000e0, 4.000000e0
124, 6.000000e0, 1.000000e0, 4.000000e0
125, 1.000000e0, 2.000000e0, 4.000000e0
126, 2.000000e0, 2.000000e0, 4.000000e0
127, 3.000000e0, 2.000000e0, 4.000000e0
128, 4.000000e0, 2.000000e0, 4.000000e0
129, 5.000000e0, 2.000000e0, 4.000000e0
130, 6.000000e0, 2.000000e0, 4.000000e0
131, 0.000000e0, 3.000000e0, 4.000000e0
132, 1.000000e0, 3.000000e0, 4.000000e0
133, 2.000000e0, 3.000000e0, 4.000000e0
134, 3.000000e0, 3.000000e0, 4.000000e0
135, 4.000000e0, 3.000000e0, 4.000000e0
136, 5.000000e0, 3.000000e0, 4.000000e0
137, 6.000000e0, 3.000000e0, 4.000000e0
138, 7.000000e0, 3.000000e0, 4.000000e0
139, 0.000000e0, 4.000000e0, 4.000000e0
140, 1.000000e0, 4.000000e0, 4.000000e0
141, 2.000000e0, 4.000000e0, 4.000000e0
142, 3.000000e0, 4.000000e0, 4.000000e0
143, 4.000000e0, 4.000000e0, 4.000000e0
144, 5.000000e0, 4.000000e0, 4.000000e0
145, 6.000000e0, 4.000000e0, 4.000000e0
146, 7.000000e0, 4.000000e0, 4.000000e0
147, 1.000000e0, 5.000000e0, 4.000000e0
148, 2.000000e0, 5.000000e0, 4.000000e0
149, 3.000000e0, 5.000000e0, 4.000000e0
150, 4.000000e0, 5.000000e0, 4.000000e0
151, 5.000000e0, 5.000000e0, 4.000000e0
152, 6.000000e0, 5.000000e0, 4.000000e0
153, 1.000000e0, 6.000000e0, 4.000000e0
154, 2.000000e0, 6.000000e0, 4.000000e0
155, 3.000000e0, 6.000000e0, 4.000000e0
156, 4.000000e0, 6.000000e0, 4.000000e0
157, 5.000000e0, 6.000000e0, 4.000000e0
158, 6.000000e0, 6.000000e0, 4.000000e0
159, 3.000000e0, 7.000000e0, 4.000000e0
160, 4.000000e0, 7.000000e0, 4.000000e0
161, 1.000000e0, 1.000000e0, 5.000000e0
162, 2.000000e0, 1.000000e0, 5.000000e0
163, 3.000000e0, 1.000000e0, 5.000000e0
164, 4.000000e0, 1.000000e0, 5.000000e0
165, 5.000000e0, 1.000000e0, 5.000000e0
166, 6.000000e0, 1.000000e0, 5.000000e0
167, 1.000000e0, 2.000000e0, 5.000000e0
168, 2.000000e0, 2.000000e0, 5.000000e0
169, 3.000000e0, 2.000000e0, 5.000000e0
170, 4.000000e0, 2.000000e0, 5.000000e0
171, 5.000000e0, 2.000000e0, 5.000000e0
172, 6.000000e0, 2.000000e0, 5.000000e0
173, 1.000000e0, 3.000000e0, 5.000000e0
174, 2.000000e0, 3.000000e0, 5.000000e0
175, 3.000000e0, 3.000000e0, 5.000000e0
176, 4.000000e0, 3.000000e0, 5.000000e0
177, 5.000000e0, 3.000000e0, 5.000000e0
178, 6.000000e0, 3.000000e0, 5.000000e0
179, 1.000000e0, 4.000000e0, 5.000000e0
180, 2.000000e0, 4.000000e0, 5.000000e0
181, 3.000000e0, 4.000000e0, 5.000000e0
182, 4.000000e0, 4.000000e0, 5.000000e0
183, 5.000000e0, 4.000000e0, 5.000000e0
184, 6.000000e0, 4.000000e0, 5.000000e0
185, 1.000000e0, 5.000000e0, 5.000000e0
186, 2.000000e0, 5.000000e0, 5.000000e0
187, 3.000000e0, 5.000000e0, 5.000000e0
188, 4.000000e0, 5.000000e0, 5.000000e0
189, 5.000000e0, 5.000000e0, 5.000000e0
190, 6.000000e0, 5.000000e0, 5.000000e0
191, 1.000000e0, 6.000000e0, 5.000000e0
192, 2.000000e0, 6.000000e0, 5.000000e0
193, 3.000000e0, 6.000000e0, 5.000000e0
194, 4.000000e0, 6.000000e0, 5.000000e0
195, 5.000000e0, 6.000000e0, 5.000000e0
196, 6.000000e0, 6.000000e0, 5.000000e0
197, 2.000000e0, 1.000000e0, 6.000000e0
198, 3.000000e0, 1.000000e0, 6.000000e0
199, 4.000000e0, 1.000000e0, 6.000000e0
200, 5.000000e0, 1.000000e0, 6.000000e0
201, 1.000000e0, 2.000000e0, 6.000000e0
202, 2.000000e0, 2.000000e0, 6.000000e0
203, 3.000000e0, 2.000000e0, 6.000000e0
204, 4.000000e0, 2.000000e0, 6.000000e0
205, 5.000000e0, 2.000000e0, 6.000000e0
206, 6.000000e0, 2.000000e0, 6.000000e0
207, 1.000000e0, 3.000000e0, 6.000000e0
208, 2.000000e0, 3.000000e0, 6.000000e0
209, 3.000000e0, 3.000000e0, 6.000000e0
210, 4.000000e0, 3.000000e0, 6.000000e0
211, 5.000000e0, 3.000000e0, 6.000000e0
212, 6.000000e0, 3.000000e0, 6.000000e0
213, 1.000000e0, 4.000000e0, 6.000000e0
214, 2.000000e0, 4.000000e0, 6.000000e0
215, 3.000000e0, 4.000000e0, 6.000000e0
216, 4.000000e0, 4.000000e0, 6.000000e0
217, 5.000000e0, 4.000000e0, 6.000000e0
218, 6.000000e0, 4.000000e0, 6.000000e0
219, 1.000000e0, 5.000000e0, 6.000000e0
220, 2.000000e0, 5.000000e0, 6.000000e0
221, 3.000000e0, 5.000000e0, 6.000000e0
222, 4.000000e0, 5.000000e0, 6.000000e0
223, 5.000000e0, 5.000000e0, 6.000000e0
224, 6.000000e0, 5.000000e0, 6.000000e0
225, 2.000000e0, 6.000000e0, 6.000000e0
226, 3.000000e0, 6.000000e0, 6.000000e0
227, 4.000000e0, 6.000000e0, 6.000000e0
228, 5.000000e0, 6.000000e0, 6.000000e0
229, 3.000000e0, 3.000000e0, 7.000000e0
230, 4.000000e0, 3.000000e0, 7.000000e0
231, 3.000000e0, 4.000000e0, 7.000000e0
232, 4.000000e0, 4.000000e0, 7.000000e0
**
*ELEMENT, TYPE=C3D8R, ELSET=EB1
1, 1, 2, 4, 3, 17, 18, 24, 23
2, 5, 6, 11, 10, 38, 39, 45, 44
3, 6, 7, 12, 11, 39, 40, 46, 45
4, 7, 8, 13, 12, 40, 41, 47, 46
5, 9, 10, 16, 15, 43, 44, 50, 49
6, 10, 11, 17, 16, 44, 45, 51, 50
7, 11, 12, 18, 17, 45, 46, 52, 51
8, 12, 13, 19, 18, 46, 47, 53, 52
9, 13, 14, 20, 19, 47, 48, 54, 53
10, 15, 16, 22, 21, 49, 50, 56, 55
11, 16, 17, 23, 22, 50, 51, 57, 56
12, 17, 18, 24, 23, 51, 52, 58, 57
13, 18, 19, 25, 24, 52, 53, 59, 58
14, 19, 20, 26, 25, 53, 54, 60, 59
15, 21, 22, 28, 27, 55, 56, 62, 61
16, 22, 23, 29, 28, 56, 57, 63, 62
17, 23, 24, 30, 29, 57, 58, 64, 63
18, 24, 25, 31, 30, 58, 59, 65, 64
19, 25, 26, 32, 31, 59, 60, 66, 65
20, 28, 29, 34, 33, 62, 63, 69, 68
21, 29, 30, 35, 34, 63, 64, 70, 69
22, 30, 31, 36, 35, 64, 65, 71, 70
23, 37, 38, 44, 43, 75, 76, 82, 81
24, 38, 39, 45, 44, 76, 77, 83, 82
25, 39, 40, 46, 45, 77, 78, 84, 83
26, 40, 41, 47, 46, 78, 79, 85, 84
27, 41, 42, 48, 47, 79, 80, 86, 85
28, 43, 44, 50, 49, 81, 82, 89, 88
29, 44, 45, 51, 50, 82, 83, 90, 89
30, 45, 46, 52, 51, 83, 84, 91, 90
31, 46, 47, 53, 52, 84, 85, 92, 91
32, 47, 48, 54, 53, 85, 86, 93, 92
33, 49, 50, 56, 55, 88, 89, 97, 96
34, 50, 51, 57, 56, 89, 90, 98, 97
35, 51, 52, 58, 57, 90, 91, 99, 98
36, 52, 53, 59, 58, 91, 92, 100, 99
37, 53, 54, 60, 59, 92, 93, 101, 100
38, 55, 56, 62, 61, 96, 97, 104, 103
39, 56, 57, 63, 62, 97, 98, 105, 104
40, 57, 58, 64, 63, 98, 99, 106, 105
41, 58, 59, 65, 64, 99, 100, 107, 106
42, 59, 60, 66, 65, 100, 101, 108, 107
43, 61, 62, 68, 67, 103, 104, 110, 109
44, 62, 63, 69, 68, 104, 105, 111, 110
45, 63, 64, 70, 69, 105, 106, 112, 111
46, 64, 65, 71, 70, 106, 107, 113, 112
47, 65, 66, 72, 71, 107, 108, 114, 113
48, 73, 74, 78, 77, 117, 118, 122, 121
49, 75, 76, 82, 81, 119, 120, 126, 125
50, 76, 77, 83, 82, 120, 121, 127, 126
51, 77, 78, 84, 83, 121, 122, 128, 127
52, 78, 79, 85, 84, 122, 123, 129, 128
53, 79, 80, 86, 85, 123, 124, 130, 129
54, 81, 82, 89, 88, 125, 126, 133, 132
55, 82, 83, 90, 89, 126, 127, 134, 133
56, 83, 84, 91, 90, 127, 128, 135, 134
57, 84, 85, 92, 91, 128, 129, 136, 135
58, 85, 86, 93, 92, 129, 130, 137, 136
59, 87, 88, 96, 95, 131, 132, 140, 139
60, 88, 89, 97, 96, 132, 133, 141, 140
61, 89, 90, 98, 97, 133, 134, 142, 141
62, 90, 91, 99, 98, 134, 135, 143, 142
63, 91, 92, 100, 99, 135, 136, 144, 143
64, 92, 93, 101, 100, 136, 137, 145, 144
65, 93, 94, 102, 101, 137, 138, 146, 145
66, 96, 97, 104, 103, 140, 141, 148, 147
67, 97, 98, 105, 104, 141, 142, 149, 148
68, 98, 99, 106, 105, 142, 143, 150, 149
69, 99, 100, 107, 106, 143, 144, 151, 150
70, 100, 101, 108, 107, 144, 145, 152, 151
71, 103, 104, 110, 109, 147, 148, 154, 153
72, 104, 105, 111, 110, 148, 149, 155, 154
73, 105, 106, 112, 111, 149, 150, 156, 155
74, 106, 107, 113, 112, 150, 151, 157, 156
75, 107, 108, 114, 113, 151, 152, 158, 157
76, 111, 112, 116, 115, 155, 156, 160, 159
77, 119, 120, 126, 125, 161, 162, 168, 167
78, 120, 121, 127, 126, 162, 163, 169, 168
79, 121, 122, 128, 127, 163, 164, 170, 169
80, 122, 123, 129, 128, 164, 165, 171, 170
81, 123, 124, 130, 129, 165, 166, 172, 171
82, 125, 126, 133, 132, 167, 168, 174, 173
83, 126, 127, 134, 133, 168, 169, 175, 174
84, 127, 128, 135, 134, 169, 170, 176, 175
85, 128, 129, 136, 135, 170, 171, 177, 176
86, 129, 130, 137, 136, 171, 172, 178, 177
87, 132, 133, 141, 140, 173, 174, 180, 179
88, 133, 134, 142, 141, 174, 175, 181, 180
89, 134, 135, 143, 142, 175, 176, 182, 181
90, 135, 136, 144, 143, 176, 177, 183, 182
91, 136, 137, 145, 144, 177, 178, 184, 183
92, 140, 141, 148, 147, 179, 180, 186, 185
93, 141, 142, 149, 148, 180, 181, 187, 186
94, 142, 143, 150, 149, 181, 182, 188, 187
95, 143, 144, 151, 150, 182, 183, 189, 188
96, 144, 145, 152, 151, 183, 184, 190, 189
97, 147, 148, 154, 153, 185, 186, 192, 191
98, 148, 149, 155, 154, 186, 187, 193, 192
99, 149, 150, 156, 155, 187, 188, 194, 193
100, 150, 151, 157, 156, 188, 189, 195, 194
101, 151, 152, 158, 157, 189, 190, 196, 195
102, 162, 163, 169, 168, 197, 198, 203, 202
103, 163, 164, 170, 169, 198, 199, 204, 203
104, 164, 165, 171, 170, 199, 200, 205, 204
105, 167, 168, 174, 173, 201, 202, 208, 207
106, 168, 169, 175, 174, 202, 203, 209, 208
107, 169, 170, 176, 175, 203, 204, 210, 209
108, 170, 171, 177, 176, 204, 205, 211, 210
109, 171, 172, 178, 177, 205, 206, 212, 211
110, 173, 174, 180, 179, 207, 208, 214, 213
111, 174, 175, 181, 180, 208, 209, 215, 214
112, 175, 176, 182, 181, 209, 210, 216, 215
113, 176, 177, 183, 182, 210, 211, 217, 216
114, 177, 178, 184, 183, 211, 212, 218, 217
115, 179, 180, 186, 185, 213, 214, 220, 219
116, 180, 181, 187, 186, 214, 215, 221, 220
117, 181, 182, 188, 187, 215, 216, 222, 221
118, 182, 183, 189, 188, 216, 217, 223, 222
119, 183, 184, 190, 189, 217, 218, 224, 223
120, 186, 187, 193, 192, 220, 221, 226, 225
121, 187, 188, 194, 193, 221, 222, 227, 226
122, 188, 189, 195, 194, 222, 223, 228, 227
123, 209, 210, 216, 215, 229, 230, 232, 231
**
*SOLID SECTION, ELSET=EB1, MATERIAL=Default-Steel
Because of its large size, the mesh structure for sphere_5 is not shown here.
Spheres - Continued
Use the fundamentals learned in the previous example to create a more sophisticated example: Concentric, high-resolution spheres consisting of three materials.
Problem Statement
Given
Given three concentric spheres of radius 10, 11, and 12 cm, as shown in the figure below,
Figure: Schematic cross-section of three concentric spheres of radius 10, 11, and 12 cm. Grid spacing is 1 cm.
Find
Use the following segmentation resolutions,
| resolution (vox/cm) | element side length (cm) | nelx | # voxels |
|---|---|---|---|
| 1 | 1.0 | 24 | 13,824 |
| 2 | 0.5 | 48 | 110,592 |
| 4 | 0.25 | 96 | 884,736 |
| 10 | 0.1 | 240 | 13,824,000 |
with a cubic domain (nelx = nely = nelz),
to create finite element meshes.
Solution
Python Segmentation
Use spheres_cont.py to create segmentations,
"""This module builds on the `spheres.py` module to create high resolution,
three-material, concentric spheres and export the voxelization as a .npy
file.
Example
-------
source ~/autotwin/automesh/.venv/bin/activate
cd ~/autotwin/automesh/book/examples/spheres_cont
python spheres_cont.py
Output
------
Saved: ~/autotwin/automesh/book/examples/spheres_cont/spheres_resolution_1.npy
Saved: ~/autotwin/automesh/book/examples/spheres_cont/spheres_resolution_2.npy
Saved: ~/autotwin/automesh/book/examples/spheres_cont/spheres_resolution_3.npy
Saved: ~/autotwin/automesh/book/examples/spheres_cont/spheres_resolution_4.npy
"""
from pathlib import Path
from typing import Final
from matplotlib.colors import LightSource
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
def sphere(resolution: int, dtype=np.uint8) -> np.ndarray:
"""Generate a 3D voxelized representation of three concentric spheres
of 10, 11, and 12 cm, at a given resolution.
Parameters
----------
resolution : int
The resolution as voxels per centimeter. Minimum value is 1.
dtype: data-type, optional
The data type of the output array. Default is np.uint8.
Returns
-------
np.ndarray
A 3D numpy array representing the voxelized spheres. Voxels within
the inner sphere are set to 1, the intermediate shell are set to 2,
and the outer shell are set to 3. Voxels outside the spheres are
set to 0.
Raises
------
ValueError
If the resolution is less than 1.
"""
print(f"Creating sphere with resolution: {resolution}")
if resolution < 1:
raise ValueError("Resolution must be >= 1")
r10 = 10 # cm
r11 = 11 # cm
r12 = 12 # cm
# We change the algorithm a bit here so we can exactly match the radius:
# number of voxels per side length (nvps)
# nvps = 2 * r12 * resolution + 1
nvps = 2 * r12 * resolution
vox_z, vox_y, vox_x = np.mgrid[
-r12: r12: nvps * 1j,
-r12: r12: nvps * 1j,
-r12: r12: nvps * 1j,
]
domain = vox_x**2 + vox_y**2 + vox_z**2
mask_10_in = np.array(domain <= r10 * r10, dtype=dtype)
mask_11_in = np.array(domain <= r11 * r11, dtype=dtype)
mask_12_in = np.array(domain <= r12 * r12, dtype=dtype)
mask_10_11 = mask_11_in - mask_10_in
mask_11_12 = mask_12_in - mask_11_in
shell_10_11 = 2 * mask_10_11
shell_11_12 = 3 * mask_11_12
result = mask_10_in + shell_10_11 + shell_11_12
# breakpoint()
print(f"Completed: Sphere with resolution: {resolution}")
return result
rr = (1, 2, 4, 10) # resolutions (voxels per cm)
lims = tuple(map(lambda x: [0, 24*x], rr)) # limits
tt = tuple(map(lambda x: [0, 12*x, 24*x], rr)) # ticks
# User input begin
spheres = {
"resolution_1": sphere(resolution=rr[0]),
"resolution_2": sphere(resolution=rr[1]),
"resolution_3": sphere(resolution=rr[2]),
"resolution_4": sphere(resolution=rr[3]),
}
aa = Path(__file__)
bb = aa.with_suffix(".png")
# Visualize the elements.
# width, height = 8, 4
width, height = 6, 3
fig = plt.figure(figsize=(width, height))
# fig = plt.figure(figsize=(8, 8))
el, az, roll = 63, -110, 0
cmap = plt.get_cmap(name="tab10")
num_colors = len(spheres)
voxel_alpha: Final[float] = 0.9
colors = cmap(np.linspace(0, 1, num_colors))
lightsource = LightSource(azdeg=325, altdeg=45) # azimuth, elevation
# lightsource = LightSource(azdeg=325, altdeg=90) # azimuth, elevation
dpi: Final[int] = 300 # resolution, dots per inch
visualize: Final[bool] = False # turn to True to show the figure on screen
serialize: Final[bool] = True # turn to True to save .png and .npy files
# User input end
N_SUBPLOTS = len(spheres)
for index, (key, value) in enumerate(spheres.items()):
if visualize:
print(f"index: {index}")
print(f"key: {key}")
print(f"value: {value}")
ax = fig.add_subplot(1, N_SUBPLOTS, index+1, projection=Axes3D.name)
ax.voxels(
value,
facecolors=colors[index],
edgecolor=colors[index],
alpha=voxel_alpha,
lightsource=lightsource)
ax.set_title(key)
# Set labels for the axes
ax.set_xlabel("x (voxels)")
ax.set_ylabel("y (voxels)")
ax.set_zlabel("z (voxels)")
ax.set_xticks(ticks=tt[index])
ax.set_yticks(ticks=tt[index])
ax.set_zticks(ticks=tt[index])
ax.set_xlim(lims[index])
ax.set_ylim(lims[index])
ax.set_zlim(lims[index])
# Set the camera view
ax.set_aspect("equal")
ax.view_init(elev=el, azim=az, roll=roll)
if serialize:
cc = aa.with_stem("spheres_" + key)
dd = cc.with_suffix(".npy")
# Save the data in .npy format
np.save(dd, value)
print(f"Saved: {dd}")
# fig.tight_layout() # don't use as it clips the x-axis label
if visualize:
plt.show()
if serialize:
fig.savefig(bb, dpi=dpi)
print(f"Saved: {bb}")
Figure: Sphere segmentations at selected resolutions, shown in the voxel domain. Because plotting large domains with Matplotlib is slow, only the first two resolutions are shown.
automesh
Use automesh to convert the segmentations into finite element meshes.
alias automesh='/Users/chovey/autotwin/automesh/target/release/automesh'
cd ~/autotwin/automesh/book/examples/spheres_cont/
automesh mesh -i spheres_resolution_1.npy \
-o spheres_resolution_1.inp \
-x 24 -y 24 -z 24 \
--xtranslate -12 --ytranslate -12 --ztranslate -12
automesh mesh -i spheres_resolution_2.npy \
-o spheres_resolution_2.inp \
-x 48 -y 48 -z 48 \
--xscale 0.5 --yscale 0.5 --yzscale 0.5 \
--xtranslate -12 --ytranslate -12 --ztranslate -12
automesh mesh -i spheres_resolution_3.npy \
-o spheres_resolution_3.inp \
-x 96 -y 96 -z 96 \
--xscale 0.25 --yscale 0.25 --zscale 0.25 \
--xtranslate -12 --ytranslate -12 --ztranslate -12
automesh mesh -i spheres_resolution_4.npy \
-o spheres_resolution_4.inp \
-x 240 -y 240 -z 240 \
--xscale 0.1 --yscale 0.1 --yzscale 0.1 \
--xtranslate -12 --ytranslate -12 --ztranslate -12
| resolution | 1 vox/cm | 2 vox/cm | 4 vox/cm | 10 vox/cm |
|---|---|---|---|---|
| midline |  |  |  |  |
| isometric |  |  |  |  |
Figure: Finite element meshes at various resolutions, shown with half-symmetric cut plane, in front view and isometric view.
Table: Summary of vital results with automesh version 0.1.10 on macOS laptop.
| resolution (vox/cm) | processing time | .npy file size | .inp file size | .g file size |
|---|---|---|---|---|
| 1 | 11.839625ms | 14 kB | 0.962 MB | 0.557 MB |
| 2 | 89.120459ms | 111 kB | 8.5 MB | 4.5 MB |
| 4 | 662.89025ms | 885 kB | 73.6 MB | 36.8 MB |
| 10 | 10.01070525s | 13.8 MB | 1.23 GB | 577 MB |
Cubit is used for the visualizations with the following recipe:
reset
import abaqus "/Users/chovey/autotwin/automesh/book/examples/spheres_cont/spheres_resolution_1.inp"
set exodus netcdf4 off
set large exodus file on
export mesh "/Users/chovey/autotwin/automesh/book/examples/spheres_cont/spheres_resolution_1.g" overwrite
reset
import mesh "/Users/chovey/autotwin/automesh/book/examples/spheres_cont/spheres_resolution_1.g" lite
graphics scale off
graphics scale on
graphics clip off
view iso
graphics clip on plane location 0 -1.0 0 direction 0 1 0
view up 0 0 1
view from 100 -100 100
graphics clip manipulation off
view bottom
Comparison
Set up reference to the Sculpt binary,
alias sculpt='/Applications/Cubit-16.14/Cubit.app/Contents/MacOS/sculpt'
automesh convert -i spheres_resolution_1.npy -o spheres_resolution_1.spn
automesh convert -i spheres_resolution_2.npy -o spheres_resolution_2.spn
automesh convert -i spheres_resolution_3.npy -o spheres_resolution_3.spn
automesh convert -i spheres_resolution_4.npy -o spheres_resolution_4.spn
Run Sculpt
Test case with letter_f_3d.spn (renamed to test_f.spn for the test below):
sculpt --num_procs 1 --input_spn "test_f.spn" \
--nelx 4 --nely 5 --nelz 3 \
--spn_xyz_order 5 \
--stair 1
Elapsed Time 0.025113 sec. (0.000419 min.)
cd ~/autotwin/automesh/book/examples/spheres_cont/
sculpt --num_procs 1 --input_spn "spheres_resolution_1.spn" \
-x 24 -y 24 -z 24 \
--xtranslate -24 --ytranslate -24 --ztranslate -24 \
--spn_xyz_order 5 \
--stair 1
sculpt --num_procs 1 --input_spn "spheres_resolution_2.spn" \
-x 48 -y 48 -z 48 \
--xscale 0.5 --yscale 0.5 --zscale 0.5 \
--xtranslate -12 --ytranslate -12 --ztranslate -12
--spn_xyz_order 5
--stair 1
sculpt --num_procs 1 --input_spn "spheres_resolution_3.spn" \
-x 96 -y 96 -z 96
--xscale 0.25 --yscale 0.25 --zscale 0.25 \
--xtranslate -12 --ytranslate -12 --ztranslate -12 \
--spn_xyz_order 5 \
--stair 1
sculpt --num_procs 1 --input_spn "spheres_resolution_4.spn" \
-x 240 -y 240 -z 240 \
--xscale 0.1 --yscale 0.1 --zscale 0.1 \
--xtranslate -12 --ytranslate -12 --ztranslate -12 \
--spn_xyz_order 5 \
--stair 1
| test | nelx | lines | time automesh | time Sculpt | automesh speed up multiple |
|---|---|---|---|---|---|
| 1 | 24 | 13,824 | 11.839625ms | 1.101862s | 93x |
| 2 | 48 | 110,592 | 89.120459ms | 3.246166s | 36x |
| 3 | 96 | 884,736 | 662.89025ms | 24.414653s | 36x |
| 4 | 240 | 13,824,000 | 10.01070525s | 449.339395s | 44x |
Smoothing
All degrees of freedom in the mesh must be in one, and only one, of the following smoothing categories:
- Prescribed
- Homogeneous
- Inhomogeneous
- Free
- Exterior
- Interface
- Interior
Figure: Two element test problem.
Table: Nodal coordinates 1-12, with x, y, z, degrees of freedom.
|node|x|y|z|->|||dof|
|:--:|:-:|:-:|:-:|:--:|::|::|:-:|
|1|0.0|0.0|0.0||1|2|3|
|2|1.0|0.0|0.0||4|5|6|
|3|2.0|0.0|0.0||7|8|9|
|4|0.0|1.0|0.0||10|11|12|
|5|1.0|1.0|0.0||13|14|15|
|6|2.0|1.0|0.0||16|17|18|
|7|0.0|0.0|1.0||19|20|21|
|8|1.0|0.0|1.0||22|23|24|
|9|2.0|0.0|1.0||25|26|27|
|10|0.0|1.0|1.0||28|29|30|
|11|1.0|1.0|1.0||31|32|33|
|12|2.0|1.0|1.0||34|35|36|
Table. The node neighbors.
| node | neighbor node(s) |
|---|---|
| 1 | 2, 4, 7 |
| 2 | 1, 3, 5, 8 |
| 3 | 2, 6, 9 |
| 4 | 1, 5, 10 |
| 5 | 2, 4, 6, 11 |
| 6 | 3, 5, 12 |
| 7 | 1, 8, 10 |
| 8 | 2, 7, 9, 11 |
| 9 | 3, 8, 12 |
| 10 | 4, 7, 11 |
| 11 | 5, 8, 10, 12 |
| 12 | 6, 9, 11 |
All Free
Following is a test where all degrees of freedom are and
hierarchical smoothing is OFF.
class DofType(Enum):
"""All degrees of freedom must belong to one, and only one, of the
following smoothing categories.
"""
PRESCRIBED_HOMOGENEOUS = 0
PRESCRIBED_INHOMOGENEOUS = 1
FREE_EXTERIOR = 2
FREE_INTERFACE = 3
FREE_INTERIOR = 4
dofset: DofSet = (
(4, 4, 4),
(4, 4, 4),
(4, 4, 4),
(4, 4, 4),
(4, 4, 4),
(4, 4, 4),
(4, 4, 4),
(4, 4, 4),
(4, 4, 4),
(4, 4, 4),
(4, 4, 4),
(4, 4, 4),
)
Table: Smoothed configuration (x, y, z).
| node | x | y | z |
|---|---|---|---|
| 1 | 0.1 | 0.1 | 0.1 |
| 2 | 1.0 | 0.075 | 0.075 |
| 3 | 1.9 | 0.1 | 0.1 |
| 4 | 0.1 | 0.9 | 0.1 |
| 5 | 1.0 | 0.925 | 0.075 |
| 6 | 1.9 | 0.9 | 0.1 |
| 7 | 0.1 | 0.1 | 0.9 |
| 8 | 1.0 | 0.075 | 0.925 |
| 9 | 1.9 | 0.1 | 0.9 |
| 10 | 0.1 | 0.9 | 0.9 |
| 11 | 1.0 | 0.925 | 0.925 |
| 12 | 1.9 | 0.9 | 0.9 |
Figure: Two element test problem (left) original configuration, (right) subject to one iteration of Laplace smoothing.
Temporary
Should move extra details from Rust API documentation to book somewhere about scaling/translation/numbering since most people will not be reading the Rust API documentation.
"The voxel data can be scaled and translated (in that order)."
Interfaces
There are three main ways to use autotwin:
- The command line interface
- The Python interface
- The Rust interface
Command Line Interface
automesh --help
@@@@@@@@@@@@@@@@
@@@@ @@@@@@@@@@
@@@@ @@@@@@@@@@@
@@@@ @@@@@@@@@@@@ automesh: Automatic mesh generation
@@ @@ @@ Chad B. Hovey <chovey@sandia.gov>
@@ @@ @@ Michael R. Buche <mrbuche@sandia.gov>
@@@@@@@@@@@@ @@@
@@@@@@@@@@@ @@@@ Notes:
@@@@@@@@@@ @@@@@ @ - Input/output file types are inferred.
@@@@@@@@@@@@@@@@ - Scaling is applied before translation.
Usage: automesh [COMMAND]
Commands:
convert Converts between segmentation input file types
mesh Creates a finite element mesh from a segmentation
smooth Applies smoothing to an existing mesh file
help Print this message or the help of the given subcommand(s)
Options:
-h, --help Print help
-V, --version Print version
Example
Convert a Numpy segmentation file to an Abaqus input file:
automesh mesh --input single.npy --output single.inp
The terminal output:
automesh 0.1.10
Reading single.npy
Done 141.917µs
Meshing single.inp
Done 155.628µs
Writing single.inp
Done 150.88µs
The resulting Abaqus input file:
*HEADING
autotwin.automesh
version 0.1.10
autogenerated on 2024-10-15 18:00:00.989487426 UTC
**
*NODE, NSET=ALLNODES
1, 0.000000e0, 0.000000e0, 0.000000e0
2, 1.000000e0, 0.000000e0, 0.000000e0
3, 0.000000e0, 1.000000e0, 0.000000e0
4, 1.000000e0, 1.000000e0, 0.000000e0
5, 0.000000e0, 0.000000e0, 1.000000e0
6, 1.000000e0, 0.000000e0, 1.000000e0
7, 0.000000e0, 1.000000e0, 1.000000e0
8, 1.000000e0, 1.000000e0, 1.000000e0
**
*ELEMENT, TYPE=C3D8R, ELSET=EB1
1, 1, 2, 4, 3, 5, 6, 8, 7
**
*SOLID SECTION, ELSET=EB1, MATERIAL=Default-Steel
Python Interface
from automesh import FiniteElements, Voxels
Example
Convert a NumPy segmentation file to an Abaqus input file:
from automesh import Voxels
voxels = Voxels.from_npy("single.npy")
fem = voxels.as_finite_elements()
fem.write_inp("single.inp")
The resulting Abaqus input file:
*HEADING
autotwin.automesh
version 0.1.10
autogenerated on 2024-10-15 18:00:01.357929060 UTC
**
*NODE, NSET=ALLNODES
1, 0.000000e0, 0.000000e0, 0.000000e0
2, 1.000000e0, 0.000000e0, 0.000000e0
3, 0.000000e0, 1.000000e0, 0.000000e0
4, 1.000000e0, 1.000000e0, 0.000000e0
5, 0.000000e0, 0.000000e0, 1.000000e0
6, 1.000000e0, 0.000000e0, 1.000000e0
7, 0.000000e0, 1.000000e0, 1.000000e0
8, 1.000000e0, 1.000000e0, 1.000000e0
**
*ELEMENT, TYPE=C3D8R, ELSET=EB1
1, 1, 2, 4, 3, 5, 6, 8, 7
**
*SOLID SECTION, ELSET=EB1, MATERIAL=Default-Steel
Rust Interface
#![allow(unused)] fn main() { extern crate automesh; use automesh::{Abaqus, FiniteElements, Voxels}; }
Example
Convert a Numpy segmentation file to an Abaqus input file:
use automesh::{Abaqus, Voxels};
fn main() {
let voxels = Voxels::from_npy("single.npy");
let scale = [1.0, 1.0, 1.0];
let translation = [0.0, 0.0, 0.0];
let fem = voxels.into_finite_elements(&scale, &translation);
fem.write_inp("single.inp");
}The resulting Abaqus input file:
*HEADING
autotwin.automesh
version 0.1.10
autogenerated on 2024-10-15 18:00:00.989487426 UTC
**
*NODE, NSET=ALLNODES
1, 0.000000e0, 0.000000e0, 0.000000e0
2, 1.000000e0, 0.000000e0, 0.000000e0
3, 0.000000e0, 1.000000e0, 0.000000e0
4, 1.000000e0, 1.000000e0, 0.000000e0
5, 0.000000e0, 0.000000e0, 1.000000e0
6, 1.000000e0, 0.000000e0, 1.000000e0
7, 0.000000e0, 1.000000e0, 1.000000e0
8, 1.000000e0, 1.000000e0, 1.000000e0
**
*ELEMENT, TYPE=C3D8R, ELSET=EB1
1, 1, 2, 4, 3, 5, 6, 8, 7
**
*SOLID SECTION, ELSET=EB1, MATERIAL=Default-Steel
Theory
This section contains a description of some of the theory used for implementation.
Smoothing
Both Laplacian smoothing and Taubin smoothing1 2 are smoothing operations that adjust the positions of the nodes in a finite element mesh.
Laplacian smoothing, based on the Laplacian operator, computes the average position of a point's neighbors and moves the point toward the average. This reduces high-frequency noise, but can result in a loss of shape and detail, with overall shrinkage.
Taubin smoothing is an extension of Laplacian smoothing that seeks to overcome the shrinkage drawback associated with the Laplacian approach. Taubin is a two-pass approach. The first pass smooths the mesh. The second pass re-expands the mesh.
Laplacian Smoothing
Consider a subject node with position . The subject node connects to neighbor points for $i \in [1, n]$ through edges.
For concereteness, consider a node with four neighbors, shown in the figure below.
Figure: The subject node with edge connections (dotted lines) to neighbor nodes with $i \in [1, n]$ (withouth loss of generality, the specific example of is shown). The average position of all neighbors of is denoted , and the gap (dashed line) originates at and terminates at .
Define as the average position of all neighbors of ,
Define the gap vector as originating at and terminating at (viz., ),
Let be the positive scaling factor for the gap .
Since
subdivision of this relationship into several substeps gives rise to an iterative approach. We typically select $\lambda < 1$ to avoid overshoot of the update, .
At iteration , we update the position of by an amount to as
with
Thus
and finally
The formulation above, based on the average position of the neighbors, is a special case of the more generalized presentation of Laplace smoothing, wherein a normalized weighting factor, , is used:
When all weights are equal and normalized by the number of neighbors, , the special case presented in the box above is recovered.
Example
For a 1D configuration, consider a node with initial position with two neighbors (that never move) with positions and (). With , the table below shows updates for for position .
Table: Iteration updates of a 1D example.
| 0 | 0.5 | 1.5 | -1 | -0.3 |
| 1 | 0.5 | 1.2 | -0.7 | -0.21 |
| 2 | 0.5 | 0.99 | -0.49 | -0.147 |
| 3 | 0.5 | 0.843 | -0.343 | -0.1029 |
| 4 | 0.5 | 0.7401 | -0.2401 | -0.07203 |
| 5 | 0.5 | 0.66807 | -0.16807 | -0.050421 |
| 6 | 0.5 | 0.617649 | -0.117649 | -0.0352947 |
| 7 | 0.5 | 0.5823543 | -0.0823543 | -0.02470629 |
| 8 | 0.5 | 0.55764801 | -0.05764801 | -0.017294403 |
| 9 | 0.5 | 0.540353607 | -0.040353607 | -0.012106082 |
| 10 | 0.5 | 0.528247525 | -0.028247525 | -0.008474257 |
Figure: Convergence of position toward as a function of iteration .
Taubin Smoothing
Taubin smoothing is a two-parameter, two-pass iterative variation of Laplace smoothing. Specifically with the definitions used in Laplacian smoothing, a second negative parameter is used, where
$$ \lambda \in \mathbb{R}^+ \subset (0, 1) \hspace{0.5cm} \rm{and} \hspace{0.5cm} \lambda < -\mu. $$
The first parameter, , tends to smooth (and shrink) the domain. The second parameter, , tends to expand the domain.
Taubin smoothing is written as, for , $k < k_{\rm{max}}$, , with an even number,
- First pass (if is even):
- Second pass (if is odd):
In any second pass (any pass with odd), the algorithm uses the updated positions from the previous (even) iteration to compute the new positions. So, the average is taken from the updated neighbor positions rather than the original neighbor positions. Some presentation of Taubin smoothing do not carefully state the second pass update, and so we emphasize it here.
Hierarchical Control
Hierarchical control classifies all nodes in a mesh as belonging to a surface , interface , or interior . These categories are mutually exclusive. Any and all nodes must belong to one, and only one, of these three categories. For a given node , let
- the set of surface neighbors be denoted ,
- the set of interface neighbors be denoted , and
- the set of interior neighbors be denoted .
Hierarchical control states that
- for any surface node , only surface nodes connected via edges to are considered neighbors,
- for any interface node , only surface nodes and interface nodes connected via edges to are neighbors, and
- for any interior node , surface nodes , interface nodes , and interior nodes connecting via edges to are neighbors.
The following figure shows this concept:
Figure: Classification of nodes into categories of surface nodes , interface nodes , and interior nodes . Hierarchical relationship: a surface node's neighbors are other other edge-connected surface nodes, an interface node's neighbors are other edge-connected interface nodes or surface nodes, and an interior node's neighbors are edge-connected nodes of any category.
Chen Example
Chen3 used medical image voxel data to create a structured hexahedral mesh. They noded that the approach generated a mesh with "jagged edges on mesh surface and material interfaces," which can cause numerical artifacts.
Chen used hierarchical Taubin mesh smoothing for eight (8) iterations, with and to smooth the outer and inner surfaces of the mesh.
References
1
Taubin G. Curve and surface smoothing without shrinkage. In Proceedings of IEEE international conference on computer vision 1995 Jun 20 (pp. 852-857). IEEE. paper
2
Taubin G. A signal processing approach to fair surface design. In Proceedings of the 22nd annual conference on Computer graphics and interactive techniques 1995 Sep 15 (pp. 351-358). paper
3
Chen Y, Ostoja-Starzewski M. MRI-based finite element modeling of head trauma: spherically focusing shear waves. Acta mechanica. 2010 Aug;213(1):155-67. paper
Development
Prerequisites
Optional
- VS Code with the following extensions
- GitHub CLI
Development Cycle Overview
- Branch
- Develop
cargo build- develop:
- tests
- implementation
- document:
mdbook build- output: automesh/book/build
mdbook serve- interactive mode
- on local machine, with Firefox, open the
index.htmlfile., e.g., file:///Users/chovey/autotwin/automesh/book/build/index.html
- test:
cargo testcargo run// test without required input and output flagscargo run --release -- -i tests/input/f.npy -o foo.exocargo run -- --help
- precommit:
pre-commit run --all-files
cargo doc --open
- Test
maturin develop --release --features python
- Merge request
Rust Development
Install Rust
Install Rust using rustup with the default standard installation:
curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh
Rust updates occur every six week. To update Rust:
rustup update
Clone Repository
git clone git@github.com:autotwin/automesh.git
Build and Test the Source Code
cd automesh
cargo test
Lint the Source Code
cargo clippy
Build and Open the API Documentation
cargo rustdoc --open -- --html-in-header docs/katex.html
Python Development
Install Rust prior to completing items below.
Install Python
Install a Python version supported by automesh.
Create a Virtual Environment
Note: If a virtual environment folder automesh/.venv already exists from previous installs, then remove it as follows:
cd automesh # from the automesh directory
(.venv) deactivate # if the virtual environment is currently active
rm -rf .venv # remove the virtual environment folder
# with `rm -rf .venv/`.
python -m venv .venv # create the virtual environment
# activate the venv with one of the following:
source .venv/bin/activate # for bash shell
source .venv/bin/activate.csh # for c shell
source .venv/bin/activate.fish # for fish shell
.\.venv\Scripts/activate # for powershell
pip install --upgrade pip
pip install maturin
Build and Test the Source Code
maturin develop --features python --extras dev
pytest
Lint the Source Code
cargo clippy --features python
pycodestyle --exclude=.venv,book,docs .
Build and Open the API Documentation
pdoc automesh --math --no-show-source --template-directory docs/
Encrypted
We create an encrypted .npy file for the sake of testing. We wish to test the case where the .npy file exists, it can be found, but it cannot be opened.
from cryptography.fernet import Fernet
import numpy as np
# Generate a key for encryption
key = Fernet.generate_key()
cipher_suite = Fernet(key)
# Create a valid .npy file
data = np.array([1, 2, 3, 4, 5])
np.save('valid.npy', data)
# Read the file content
with open('valid.npy', 'rb') as f:
file_data = f.read()
# Encrypt the file content
encrypted_data = cipher_suite.encrypt(file_data)
# Write the encrypted data to a new file
with open('encrypted.npy', 'wb') as f:
f.write(encrypted_data)
print(f"Encryption key: {key.decode()}")