The Constructive Solid Geometry (CSG) package

The csg package is intended for the universal creation of geometry models based on regularizing set operations. Solid geometry modeling is quite common in Monte Carlo neutronics and shielding codes, and the basic principles should be familiar to anybody with experience in codes such as MCNP and SERPENT. The main difference in our approach is that we emphasize solids over surfaces. More specifically, the geometry is encoded directly in the volumes, and there is no need for tedious book keeping of interface surfaces.

Note

From a practical stand point, this means that you do not need to have a defined pool of surfaces (or primitives), but can construct them as you need them, even if they are used by many cells. The system will automatically deduce shared surfaces. Thus, regions (or volumes) can be constructed completely independently from one another, allowing one to focus on a particular part of an assembly, without having to take the entire construction into consideration.

However, if you are more comfortable with the traditional approach, you can still make use of a single pool of surfaces (or primitives).

The package has two basic layers: The first simply provides a convenient, code independent interface for creating and storing CSG geometries. This layer does not perform any geometric operations, it merely stores the construction tree.

The second layer uses a polyhedra package to perform the actual geometric operations. This layer is optional, but is currently required for the geometry visualization pipeline, as well as numerous geometric algorithms (e.g. cell simplifications, meshing, volume calculations etc).

This package creates basic volumes (bounded or unbounded regions in 3D space), which can later be used to construct physical cells (the primary building blocks of components), by adding data describing its purpose (e.g fill it with a material, or burnable bundle). Components of the csg package are typically used in the following order:

Choose a set of primitive regions that best capture the region one wishes to construct.

Use transformations (rotations and translations) to correctly orient the primitives.

Combine the primitives using set operations.

Primitive regions

Primitive regions constitute the basic building blocks out of which more complicated regions are constructed. Since these regions are bounded by simple (and often parametric) surfaces, they play the same role as surface definitions in Monte Carlo codes.

All basic primitives are defined in the csg.primitives module, while some extensions are defined in the csg.region_macros module.

Note

All primitives have a show() method, which opens a basic three dimensional view of the region.

Planes

These are infinite, half space regions bounded by a plane. The following plane constructions are available:

primitives.Px(position=p)

Region below a plane perpendicular the \(x\)-axis.

Parameters:: position (length) – Intersection point with the \(x\)-axis.

Parametric set definition:

\[\{(x,y,z) : x \leq p\}\]

Example:

>>> p = primitives.Px(2 * units.cm)

primitives.Py(position=p)

Region below a plane perpendicular the \(y\)-axis.

Parameters:: position (length) – Intersection point with the \(y\)-axis.

Parametric set definition:

\[\{(x,y,z) : y \leq p\}\]

Example:

>>> p = primitives.Py(1 * units.cm)

primitives.Pz(position=p)

Region below a plane perpendicular the \(z\)-axis.

Parameters:: position (length) – Intersection point with the \(z\)-axis.

Parametric set definition:

\[\{(x,y,z) : z \leq p\}\]

Example:

>>> p = primitives.Pz(3 * units.cm)

primitives.Plane(a=a, b=b, c=c, d=p)

Region below an arbitrary plane.

Parameters:

a (float) – \(x\) component of plane normal.
b (float) – \(y\) component of plane normal.
c (float) – \(z\) component of plane normal.
d (length) – Position of plane.

Parametric set definition:

\[\{(x,y,z) : ax + by + cz \leq p\}\]

Example:

>>> p = primitives.Plane(a=1.0, b=1.0, c=0.0, d=-1.0 * units.cm)

primitives.create_plane(p0, p1, p2)

Create a plane through 3 (non co-linear) points.

Parameters:

p0 (point) – First point.
p1 (point) – Second point.
p2 (point) – Third point.

Example:

>>> p = primitives.create_plane((1.0 * units.cm, 0, 0), (0, 1 * units.cm, 0), (0, 0, 1 * units.cm))
>>> print(p.a, p.b, p.c)
1.0 1.0 1.0

Cylindrical regions

These are regions that are bounded in one coordinate plane by circular shapes, and is infinite in the perpendicular dimension.

primitives.Cx(radius=r (OR diameter=d), center=(y0, z0))

Region inside an infinite cylinder parallel to the \(x\)-axis.

Parameters:

radius (length) – Radius of the cylinder.
diameter (length) – Diameter of the cylinder (instead of radius, but not both).
center (point) – Center of the cylinder in the \(yz\)-plane. Default value is \((0,0)\).

Parametric set definition:

\[\left\{(x,y,z) : (y-y_0)^2 + (z-z_0)^2 \leq r^2\right\}\]

Examples:

>>> p0 = primitives.Cx(radius=1.0 * units.cm) # centered at (0,0)
>>> p1 = primitives.Cx(diameter=2.0 * units.cm, center=(-1.0 * units.cm, 0.0))

primitives.Cy(radius=r (OR diameter=d), center=(x0, z0))

Region inside an infinite cylinder parallel to the \(y\)-axis.

Parameters:

radius (length) – Radius of the cylinder.
diameter (length) – Diameter of the cylinder (instead of radius, but not both).
center (point) – Center of the cylinder in the \(xz\)-plane. Default value is \((0,0)\).

Parametric set definition:

\[\left\{(x,y,z) : (x-x_0)^2 + (z-z_0)^2 \leq r^2\right\}\]

Examples:

>>> p0 = primitives.Cy(radius=1.0 * units.cm) # centered at (0,0)
>>> p1 = primitives.Cy(diameter=2.0 * units.cm, center=(-1.0 * units.cm, 0.0))

primitives.Cz(radius=r (OR diameter=d), center=(x0, y0))

Region inside an infinite cylinder parallel to the \(z\)-axis.

Parameters:

radius (length) – Radius of the cylinder.
diameter (length) – Diameter of the cylinder (instead of radius, but not both).
center (point) – Center of the cylinder in the \(xy\)-plane. Default value is \((0,0)\).

Parametric set definition:

\[\left\{(x,y,z) : (x-x_0)^2 + (y-y_0)^2 \leq r^2\right\}\]

Examples:

>>> p0 = primitives.Cz(radius=1.0 * units.cm) # centered at (0,0)
>>> p1 = primitives.Cz(diameter=2.0 * units.cm, center=(-1.0 * units.cm, 0.0))

Note

There is no arbitrary orientated cylinder primitive. Choose one of the above and add the appropriate transformation.

region_macros.pipe(diameter, width, center=(0, 0), alignment=primitives.ALIGNMENT.Z_AXIS)

Cylindrical pipe, that is, region between two cylinders.

Parameters:

diameter (length) – Diameter of the outer (bounding) cylinder.
width (length) – Width of the pipe sleeve. The inner diameter is calculated as diameter less two times this value.
center (point) – Center of the pipe in the \(uv\)-plane.
alignment – Axial axis of the pipe.

Parametric set definition:

\[\left\{(u,v,a) : \left(\frac{d}{2} - w\right)^2 \leq (u-u_0)^2 + (v-v_0)^2 \leq \frac{d}{2}^2\right\},\]

where \(d\) is diameter, and \(w\) is width.

Note

The pipe is unbounded in the axial direction.

primitives.wedge(radius, arc_length, angle=0.0 * units.degrees, alignment=primitives.ALIGNMENT.Z_AXIS, center=(0, 0), no_cylindrical_bound=False)

Region bounded by a circular segment in the \(uv\)-plane.

Parameters:

radius (length) – Radius of the circle.
arc_length (length) – Total length of the segment arc.
angle – Center angle of the arc.
alignment – Axial axis of the wedge.
center (point) – Center of the circle in the \(uv\)-plane.
no_cylindrical_bound (bool) – Flag indicating if the circular bound should be included.

In polar coordinates,

\[u-u_0 = r \cos(\theta), v-v_0 = r \sin(\theta)\]

the region is described by

\[\left\{(r, \theta) : 0 \leq r \leq r_b, |\theta - \theta_0| \leq \frac{l}{2 r_0}\right\},\]

where \(r_0\) is radius, \(\theta_0\) is angle, abd \(l\) is arc_length. The radial upper bound \(r_b\) is radius if no_cylindrical_bound is False, \(\infty\) otherwise.

Rectangular regions

These are axially aligned regions that are bounded in one coordinate plane by a rectangle, and either infinite or bounded in the perpendicular directions. The alignment parameter determines the axes the primitive is parallel to (i.e. the infinite direction), and is set by choosing the appropriate value from the ALIGNMENT container. It always defaults to ALIGNMENT.Z_AXIS if not specified.

In all cases \((u,v)\) is used to denote coordinates in the perpendicular plane:

Alignment

Perpendicular Plane

X_AXIS (\(a=x\))

\(yz\)-plane (\(u=y, v=z\))

Y_AXIS (\(a=y\))

\(xz\)-plane (\(u=x, v=z\))

Z_AXIS (\(a=z\))

\(xy\)-plane (\(u=x, v=y\))

primitives.SquareCylinder(radius=r, center=(u0, v0), alignment=ALIGNMENT.Z_AXIS)

Region bounded by a square in the \(uv\)-plane.

Parameters:

radius (length) – Half the length (or width) of the square.
center (point) – Center of the square in the \(uv\)-plane. Default is \((0,0\)).
alignment – Axis parallel to infinite extension.

Parametric set definition:

\[\{(u,v,a) : |u-u_0| \leq r, |v-v_0| \leq r\}\]

Examples:

>>> p0 = primitives.SquareCylinder(radius=1.0 * units.cm, alignment=ALIGNMENT.X_AXIS)
>>> p1 = primitives.SquareCylinder(radius=1.0 * units.cm, alignment=ALIGNMENT.Z_AXIS, center=(0.0, -1.0 * units.mm))

primitives.RectangularCylinder(width=du, length=dv, center=(u0, v0), alignment=ALIGNMENT.Z_AXIS)

Region bounded by a rectangle in the \(uv\)-plane.

Parameters:

width (length) – Length of the square in the \(u\) coordinate direction.
length (length) – Length of the square in the \(v\) coordinate direction.
center (point) – Center of the square in the \(uv\)-plane. Default is \((0,0\)).
alignment – Axis parallel to infinite extension.

Parametric set definition:

\[\left\{(u,v,a) : |u-u_0| \leq \frac{du}{2}, |v-v_0| \leq \frac{dv}{2}\right\}\]

Examples:

>>> p0 = primitives.RectangularCylinder(width=1.0 * units.cm, length=2.0 * units.cm, alignment=ALIGNMENT.X_AXIS)

primitives.Cube(radius=r, center=(x0, y0, z0))

Region bounded by an axis aligned cube.

Parameters:

radius (length) – Half the length of any side.
center – Center of the cube. Default value in \((0,0,0)\).

\[\{(x,y,z) : |x-x_0| \leq r, |y-y_0| \leq r, |z-z_0| \leq r\}\]

Examples:

>>> p = primitives.Cube(radius=1.0 * units.cm, center=(0.0, 1.0 * units.mm, 0.0))

primitives.Cuboid(x_min, x_max, y_min, y_max, z_min, z_max)

Region bounded by a cuboid.

Parameters:

x_min (length) – Minimum \(x\)-coordinate.
x_max (length) – Maximum \(x\)-coordinate.
y_min (length) – Minimum \(y\)-coordinate.
y_max (length) – Maximum \(y\)-coordinate.
z_min (length) – Minimum \(z\)-coordinate.
z_max (length) – Maximum \(z\)-coordinate.

Parametric set definition:

\[\{(x,y,z) : x_{min} \leq x \leq x_{max}, y_{min} \leq y \leq y_{max}, z_{min} \leq z \leq z_{max}\}\]

Examples:

>>> p = primitives.Cuboid(-1.0 * units.cm, 0.0 * units.cm, 0.0 * units.cm, 1.0 * units.cm, -1.0 * units.cm, 1.0 * units.cm)

Conic regions

primitives.Cone(scaling=t, focus=a0, center=(u0, v0), alignment=ALIGNMENT.Z_AXIS, side=0)

Basic single or two sided cone.

Parameters:

scaling (float) – Rate at which the cone decreases.
focus (length) – Cone focus point in the axial coordinate direction.
center (point) – Center of cone in the \(uv\)-plane. Default value in \((0,0)\).
alignment – Axis parallel to cone.
side (int) – Flag indicating which side(s) of the cone should be included. Use 0 for both sides, -1 for the side below the focus, and +1 for the side above the focus. The default value is 0.

Parametric set definition:

\[\left\{(u, v, a): \sqrt{(u-u0)^2 + (v-v0)^2} \leq \sqrt{t} (a-a0)\right\}\]

primitives.conic_section(large_circle_radius, small_circle_radius, height, center=(0, 0, 0), alignment=ALIGNMENT.Z_AXIS)

Conic section between two circles.

Parameters:

large_circle_radius (length) – Radius of the larger (lower) circle.
small_circle_radius (length) – Radius of the small (upper) circle.
height (length) – Total height of the section.
center (point) – Center of mass of the section.
alignment – Axis parallel to the height of the section.

Example:

>>> p = primitives.conic_section(2.0 * units.cm, 1.5 * units.cm, 2.0 * units.cm)
>>> p.show()

primitives.rectangular_cone(base_rectangle_width, base_rectangle_length, upper_rectangle_width, upper_rectangle_length, height, center=(0, 0, 0), alignment=ALIGNMENT.Z_AXIS)

Pyramid like region connecting two rectangles.

Parameters:

base_rectangle_width (length) – Length along the \(u\) coordinate of the base (lower) rectangle.
base_rectangle_length (length) – Length along the \(v\) coordinate of the base (lower) rectangle.
upper_rectangle_width (length) – Length along the \(u\) coordinate of the upper rectangle.
upper_rectangle_length (length) – Length along the \(v\) coordinate of the upper rectangle.
height (length) – Distance between lower and upper rectangle.
center (point) – Center of mass of the entire region.
alighnment – Axis parallel to the height of the pyramid.

Examples:

>>> p = primitives.rectangular_cone(5 * units.cm, 4 * units.cm, 3 * units.cm, 2 * units.cm, height=3 * units.cm)
>>> p.show()

Spherical regions

primitives.Sphere(radius=r, center=(x0, y0, z0))

Three dimensional spherical region.

Parameters:

radius (length) – Radius of the sphere.
center (point) – Center of the sphere. Default is \((0,0,0)\).

Parametric set definition:

\[\left\{(x,y,z) : (x-x_0)^2 + (y-y_0)^2 + (z-z_0)^2 \leq r^2\right\}\]

primitives.Ellipsoid(small_radius, large_radius, center=(0, 0, 0), orientation=ALIGNMENT.Z_AXIS)

Three dimensional region bounded by an ellipsoid.

Parameters:

small_radius (length) – Small radius of the ellipsoid (short axis).
large_radius (length) – Large radius of the ellipsoid (long axis).
center (point) – Center of the ellipsoid.
orientation – Direction of the long axis, that is, where the scaling should take place. See the description below.

The orientation parameter can have any of the following values:

A coordinate axis (i.e ALIGNMENT.X_AXIS, ALIGNMENT.Y_AXIS or ALIGNMENT.Z_AXIS) in which case the long axis will be aligned with this axis.

A string specifying a coordinate plane (‘xy’, ‘xz’ or ‘yz’). In this case, the short axis will align with the coordinate axis perpendicular to the specified plane.

An arbitrary point, in which case the long axis will be aligned with this vector.

Note

Currently, an ellipsoid is not an independent primitive, but merely a transformed (scaled) version of primitives.Sphere().

Polyhedral regions

primitives.Hexagon(diameter=None, radius=None, inner_radius=None, orientation=ALIGNMENT.X_AXIS, center=(0, 0), alignment=ALIGNMENT.Z_AXIS)

Create a region bounded by a hexagonal prism.

Parameters:

diameter (length) – Outer diameter of the hexagon (that is, all vertices are on this circle).
outer_radius (length) – Outer radius of the hexagon (that is, all vertices are on this circle).
inner_radius (length) – Inner radius of the hexagon (that is, all edges are tangent to this circle)
orientation – On what axis a vertex will be placed. If ALIGNMENT.X_AXIS, a vertex will be placed on the \(u\)-axis creating a flat top hexagon. Otherwise, if ALIGNMENT.Y_AXIS, a vertex will be placed on the \(v\)-axis creating a pointy top hexagon.
alignment – Axial alignment of the hexagonal region.
center – Center of the region in the \(uv\)-plane.

Note

Only one of the parameters diameter, radius or inner_radius should be specified.

region_macros.extruded(points, height=1.0 * units.cm, direction=primitives.ALIGNMENT.Z_AXIS, include_top=True, include_bottom=True)

Create a three dimensional region by extruding a set of co-planar point enclosing a convex polygon.

Parameters:

points (list (point)) – List of co-planar points.
height (length) – Height of the extrusion.
direction – Direction along which the points should be extruded. Either ALIGNMENT.X_AXIS, ALIGNMENT.Y_AXIS, or ALIGNMENT.Z_AXIS for extrusion in a coordinate direction, or an arbitrary point.
include_top (bool) – Flag indicating if the extrusion should be bounded from above. If False, the region will be infinite in the positive alignment direction.
include_bottom (bool) – Flag indicating if the extrusion should be bounded from above. If False, the region will be infinite in the negative alignment direction.

Attention

The points must be ordered clockwise with respect to direction (left hand rule).
The last point in the list will be connected to the first point to create a closed cycle.
For the extrusion to be perpendicular, all points must be in the plane perpendicular to direction.
If direction is a point, the height parameter is ignored, and the height of the extrusion is equal to the magnitude of the direction.

Example:

>>> points = [(-1 * units.cm, 0), (1 * units.cm, 0), (0.5 * units.cm, 1 * units.cm), (-0.5 * units.cm, 1 * units.cm)]
>>> points.reverse()     # correct orientation
>>> p = region_macros.extruded(points, height=2* units.cm)
>>> p.show()

Strips

Basic semi-infinite regions bounded by axially aligned planes.

region_macros.axial_strip(lower=None, upper=None)

Infinite region bounded by two axial (\(z\)-axes) planes.

Parameters:

lower (length or None) – Lower axial bound. If None, region will not be bounded from below.
upper (length or None) – Upper axial bound. If None, region will not be bounded from above.

Parametric set definition:

\[\left\{(x,y,z) : l < z < u \right\}\]

where \(l\) is lower or \(-\infty\) if lower is None, and \(u\) is upper or \(+\infty\) if upper is None.

region_macros.x_strip(xmin=None, xmax=None)

Infinite region bounded by two \(x\)-axes planes.

Parameters:

xmin (length or None) – Lower bound. If None, region will not be bounded from below.
xmax (length or None) – Upper bound. If None, region will not be bounded from above.

Parametric set definition:

\[\left\{(x,y,z) : l < x < u \right\}\]

where \(l\) is xmin or \(-\infty\) if xmin is None, and \(u\) is xmax or \(+\infty\) if xmax is None.

region_macros.y_strip(ymin=None, ymax=None)

Infinite region bounded by two \(y\)-axes planes.

Parameters:

ymin (length or None) – Lower bound. If None, region will not be bounded from below.
ymax (length or None) – Upper bound. If None, region will not be bounded from above.

Parametric set definition:

\[\left\{(x,y,z) : l < y < u \right\}\]

where \(l\) is ymin or \(-\infty\) if ymin is None, and \(u\) is ymax or \(+\infty\) if ymax is None.

region_macros.rectangle(xmin=None, xmax=None, ymin=None, ymax=None)

Rectangular region in the radial \(xy\)-plane. Unbounded in the axial direction.

Parameters:

xmin (length or None) – Minimum \(x\)-coordinate. If None, the region will be unbounded on the left.
xmax (length or None) – Maximum \(x\)-coordinate. If None, the region will be unbounded on the right.
ymin (length or None) – Minimum \(y\)-coordinate. If None, the region will be unbounded on the bottom.
ymax (length or None) – Maximum \(y\)-coordinate. If None, the region will be unbounded on the top.

Parametric set definition:

\[\left\{(x,y,z) : x_{min} < x < x_{max}, y_{min} < y < y_{max} \right\}\]

where \(x_{min}\), \(y_{min}\) is \(-\infty\) if xmin, ymin is None respectively, and \(x_{max}\), \(y_{max}\) is \(+\infty\) if xmax, ymax is None respectively.

Other radial regions

These regions are bounded by specified shapes in the \(xy\) plane, but are infinite in the axial direction.

region_macros.clipped_box(diameter, cut)

Square in the \(xy\)-plane with clipped corners.

Parameters:

diameter (length) – Total width of the square.
cut (length) – Piece that should be clipped from each side, so that the remaining edge has length diameter less two times cut.

Note

The square is centered at \(x=0, y=0\).

Example:

>>> p = region_macros.clipped_box(10 * units.cm, 2.0 * units.cm)
>>> p.show()

region_macros.square_cross(radius, width, center=(0, 0))

Square cruciform region.

Parameters:

radius (length) – Distance from center to cross end.
width (length) – Width of the cross pieces.
center (point) – Center of the cross in the \(xy\)-plane.

Example:

>>> p = region_macros.square_cross(5 * units.cm, 2.0 * units.cm)
>>> p.show()

Trivial regions

primitives.Space(): Region representing the entire space.

primitive.Empty(): Region representing an empty space (null set).

Transformations

Since most of the available primitives are aligned to the coordinate axis, transformations (rotations and translations) must be added in order to capture the correct region. This is achieved with the add_transformation method, common to all regions (primitive or not). The parameter passed to this method must be a valid affine transformation.

In it’s most general from, and affine transformation is a \(4\times4\) matrix, with the first \(3\times 3\) entries \(R_{ij}\) specifying the rotation matrix, and the translation \(b_j\) specified in column 4:

\[\begin{split}\begin{bmatrix} R_{11} & R_{12} & R_{13} & b_1 \\ R_{21} & R_{22} & R_{23} & b_2 \\ R_{31} & R_{32} & R_{33} & b_3 \\ 0 & 0 & 0 & 1 \\ \end{bmatrix}\end{split}\]

Affine transformations transform points \(r=(x,y,z)\) to \(r'=(x',y',z')\) as follows:

\[r' = R r + b.\]

Note that translations are applied after the rotation. The package allows the user full freedom in specifying the affine transformation matrix, but the core.geometry module contains a number of convenience functions for creating the most commonly used ones.

Basic translations

geometry.translation(x=0.0, y=0.0, z=0.0)

Creates a rotation free translation matrix.

Parameters:

x (length) – The \(x\)-coordinate of the translation vector.
y (length) – The \(y\)-coordinate of the translation vector.
z (length) – The \(z\)-coordinate of the translation vector.

Example:

>>> p = primitives.Cx(radius=1.0 * units.cm)
>>> p.add_transformation(geometry.translation(1.0 * units.cm, 2.0 * units.cm, 1.0 * units.cm))

Rotation around an arbitrary vector

geometry.rotation_around_vector(angle, vector=(0, 0, 1))

Create a pure rotation transformation, which rotates all points by a angle around vector.

Parameters:

angle – Angle of rotation (in degrees or radians).
vector – Either a three coordinate vector, or an axis (ALIGNMENT.X_AXIS, ALIGNMENT.Y_AXIS or ALIGNMENT.Z_AXIS).

Example:

>>> p = primitives.Cx(radius=1.0 * units.cm)
>>> p.add_transformation(geometry.rotation_around_vector(33.0 * units.degrees, ALIGNMENT.Z_AXIS))

Note

The rotation is anchored at the origin \(x=0, y=0, z=0\).

Attention

It is important to use units.degrees when specifying angles in degrees. Otherwise the system will assume it is in radians!

Euler rotations

geometry.euler_rotations(roll=0.0, pitch=0.0, yaw=0.0)

Creates a rotation based on euler operations.

Parameters:

roll – Roll angle (rotation around \(x\)-axis).
pith – Pitch angle (rotation around \(y\)-axis).
yaw – Yaw angle (rotation around \(z\)-axis).

Example:

>>> p = primitives.Cx(radius=1.0 * units.cm)
>>> p.add_transformation(geometry.euler_rotations(pitch=10.0 * units.degrees, yaw=45.0 * units.degrees))

Combining affine transformations

Affine transformations are automatically combined when multiple transformations are added to a region. For instance, the following will rotate and then translate a region p:

p.add_transformation(geometry.rotation_around_vector(33.0 * units.degrees, [1,0,0]))
p.add_transformation(geometry.translation(1.0 * units.cm, 2.0 * units.cm, 1.0 * units.cm))

Attention

Affine transformations are, in general, non-commutative. Thus, the order in which transformations are added is important!

Regularizing set operations

Regularizing set operations (complement, intersection and union) are used to construct more complicated regions out of primitives. The CSG package allows two notation styles when specifying set operations:

Set operation style: ~ for complement, & for intersections and | for unions.

MCNP style: + for complement, space for intersections and : for unions.

Note

The MCNP style notation is provided mostly for people with years of experience with MCNP input decks. However, we recommend adapting the set operation notation, as it is more expressive, and the entire validation example set was constructed using this notation.

Complement

The complement of a region denotes everything not in the interior of the region.

inside = primitives.Cylinder(radius=1.0 * units.cm)
outside = ~inside

Or, equivalently, in the MCNP notation:

inside = primitives.Cylinder(radius=1.0 * units.cm)
outside = +inside

Intersection

The intersection between two regions r1 and r2 denotes everything that is in r1 and r2. For example, the following describe the region inside a cube and outside a cylinder:

r1 = primitives.Cube(radius=2.0 * units.cm)
r2 = primitives.Cylinder(radius=1.0 * units.cm)
region = r1 & (~r2)      #reads: inside cube AND NOT inside cylinder
region.show()

Alternatively, in the MCNP style notation, intersections is not specified by a symbol:

r1 = primitives.Cube(radius=2.0 * units.cm)
r2 = primitives.Cylinder(radius=1.0 * units.cm)
region = -r1 +r2        # reads: inside cube AND outside cylinder
region.show()

Attention

Although the + symbol is optional in MCNP cell cards, it cannot be omitted here!

Union

The union of two regions r1 and r2 denotes everything that is in r1 or r2. The following example illustrates how region outside a rectangle can be constructed as the union of half plane regions:

outside_rect = primitives.Px(-1.0 * units.cm) | ~primitives.Px(1.0 * units.cm)  \
               | primitives.Py(-1.0 * units.cm) | ~primitives.Py(1.0 * units.cm)

Equivalently, in MCNP style notation:

outside_rect = -primitives.Px(-1.0 * units.cm) : +primitives.Px(1.0 * units.cm) \
                : -primitives.Py(-1.0 * units.cm) : +primitives.Py(1.0 * units.cm)

Set operations on composite regions

The above set operations are not restricted to primitives, and can be applied to any region. For example:

>>> from csg import *
>>> r1 = primitives.Cube(radius=2.0 * units.cm) & (~primitives.Cylinder(radius=1.0 * units.cm))
>>> co = (-1.5 * units.cm, 0.0 * units.cm, 0.0 * units.cm)
>>> r2 = primitives.Cube(radius=2.0 * units.cm, center=co) & ~primitives.Cylinder(radius=1.0 * units.cm, center=co)
>>> r3 = r1 & ~r2
>>> r3.show()

Additional region functions

The following are some useful functions when working with regions:

region.intersects(r1, r2)

Check if two regions intersects.

Parameters:

r1 – First region (either a primitive or composite).
r2 – Second region (either a primitive or composite).

Returns:

True if the intersection of r1 and r2 has non-empty interior, False otherwise.

region.volume(r)

Computes the volume of th region.

Parameters:: r – Region to calculate volume of. Either a primitive or region constructed through set operations.
Returns:: The volume of the region (a volume quantity).

Note

Infinite regions will be clipped to the universe box, so the returned volume will be finite. Thus, a meaningful value can only be expected for bounded regions.

region.decompose(r)

Decompose a region into disjoint pieces. This function can be used to create simpler regions from a region constructed through a complicated sequence of set operations.

Parameters:: r – Region to decompose.
Returns:: A list of mutually disjoint regions.

Repeated structures

The csg package has no predefined set of repeated structures. However, the entire python machinery is at your disposal. Using loops that update transformations, one can recreate any lattice type structure (and much more!)

The following example builds a complicated region by repeatedly adding rotated cylinders:

from csg import *

axial = ~primitives.Pz(-2.0 * units.cm) & primitives.Pz(2.0 * units.cm)

prim = primitives.Cylinder(radius=0.75 * units.cm, center=(0.0 * units.cm, 1.0 * units.cm)) & axial

# Initiate with a copy of the first cylinder
region = prim.clone()

for i in range(0, 7):
    prim.add_transformation(geometry.rotation_around_vector(45 * units.degrees, [0, 0, 1]))
    region |= prim

region.show()

Working with cells

Cells form the building blocks of all assemblies. They are essentially regions (bounded or unbounded) in three dimensional space, with additional data connecting them to the physical structure of an assembly.

cell.Cell(structure, material=None, facility=None, bundle=None, description=None, part=None, outside=False)

Create a physical cell.

Parameters:

structure – Any region (either a primitive, or composite region constructed using set operations).
material (material) – Material that fills the region.
facility (str) – Mark this cell as containing the specified facility. This means that the contents of this cell is controlled by the assembly.
bundle (bundle_tag) – Indicate that the cell contains a burn bundle. This means that the contents of the cell is controlled by the specified burnable material manager.
description (str) – Short description of what the cell is modeling.
part (str) – Assembly part this cell belongs to. Parts are used to group cells into real physical parts, and is required when using the CAD link when documenting models.
outside (bool) – Flag indicating if the cell falls outside the region of interest. That is, cell with outside True is considered to lie outside the boundary of the problem domain.

Note

All these keyword parameters (except structure) are also attributes of the cell container itself, and can be assigned after construction, e.g

>>> c = cell.Cell(my_region)
>>> c.material = my_material

The following method is used to retrieve the underlying structure of the cell:

cell.region()

Returns:: The underlying structure of the cell. Either a primitive or composite (made from CSG operations) region.

This method can be used to emulate the MCNP # operator. For example, assuming c1 is an existing cell,

>>> r = ~c1.region()

gives the region not in c1.

Cell operations

cell.reflect(axis, *cells)

Reflect a number of cells through a coordinate plane.

Parameters:

axis (int) –
Reflection axis. Use
- ALIGNMENT.X_AXES for reflection through the \(yz\)-plane,
- ALIGNMENT.Y_AXES for reflection through the \(xz\)-plane, and
- ALIGNMENT.Z_AXES for reflection through the \(xy\)-plane.
cells – List of cells to reflect.

Returns:

A new list of reflected cells (equal to the number of input cells).

This function is very useful when modeling symmetric components, and effectively reduces the number of cells that need to be defined by half.

Example:

>>> from csg import *
>>> r = region_macros.rectangle(0.5 * units.mm, 1.5 * units.mm, -1.0 * units.mm, 1.0 * units.mm)
>>> c1 = cell.Cell(r)
>>> rc = cell.reflect(ALIGNMENT.X_AXIS, c1)
>>> (c1 | rc[0]).show()

Composites

The csg.composites and csg.adapters contain a number of structures which can create sets of cells. These, can be viewed as a pre-defined collection of parts which can be used in your models. Since the composites are usually defined by a large number parameters, they are used by first instantiating the data container, setting its parameters, then calling the build() method to construct the cells.

class composites.RoundedSquare

Square region with rounded (circular) corners.

alignment : integer = ALIGNMENT.Z_AXIS: Coordinate axis along which the height of the cells extend. Either ALIGNMENT.X_AXIS, ALIGNMENT.Y_AXIS or ALIGNMENT.Z_AXIS.

center : point = :math:`(0,0,0)`: Center of the part.

height : length !: Height of the part along the axes defined by alignment.

square_radius : length !: Total radius of the square (that is, half the width).

corner_radius : length !: Radius of curvature of the corners.

Example:

>>> from csg import *
>>> rs = composites.RoundedSquare()
>>> rs.square_radius = 2.0 * units.cm
>>> rs.corner_radius = 0.5 * units.cm
>>> rs.height = 5.0 * units.cm
>>> cells = rs.build()
>>> cells.show()

class composites.RoundedSquareTube

Strip between two square regions with rounded (circular) corners. Accepts all the parameters defined in composites.RoundedSquare as well as:

thickness : length !: Width of the tubular region.

Note

The parameters square_radius, corner_radius refers to the dimensions of the outer square.

Example:

>>> rs = composites.RoundedSquareTube()
>>> rs.square_radius = 2.0 * units.cm
>>> rs.corner_radius = 0.5 * units.cm
>>> rs.height = 5.0 * units.cm
>>> rs.thickness = 1.0 * units.mm
>>> cells = rs.build()
>>> cells.show()

class composites.RoundedSquareCutOut

Region between a square and a square with rounded corners. Accepts all the parameters defined in composites.RoundedSquare as well as:

outside_square_radius : length !: Half the width of the outer (boundary) square region.

Note

In this case, the parameters square_radius, corner_radius refers to the dimensions of the inner square.

Example:

>>> rs = composites.RoundedSquareCutOut()
>>> rs.square_radius = 2.0 * units.cm
>>> rs.corner_radius = 0.5 * units.cm
>>> rs.height = 5.0 * units.cm
>>> rs.outside_square_radius = 2.0 * units.cm + 1.0 * units.mm
>>> cells = rs.build()
>>> cells.show()

class composites.PinCell

Two dimensional pin cell definition. This composites adds consecutive rings, creating a nested set of tubular cells.

The composite accepts the following parameter:

fill_material : material !: Material that fills space outside the last ring.

Consecutive rings are added using the following three methods:

add_cylindrical_shell(radius, material=None, bundle=None, facility=None, part=None):

Adds a cylindrical layer.

Parameters:

radius (length) – Radius of the layer.
material (material) – Material that fills the region between the cylinder with radius and the previously added layer.
facility (str) – Facility that the region between the cylinder with radius and the previously added layer contain.
bundle (bundle_tag) – Indicate that the current cell contains a burn bundle.
part (str) – Assembly part the current cell belongs to.

Note

If this is the first layer added, the region will be a complete cylinder.

add_square_shell(radius, material=None, bundle=None, facility=None, part=None):

Adds a square layer.

Parameters:

radius (length) – Radius of the layer.
material (material) – Material that fills the region between the square with radius and the previously added layer.
facility (str) – Facility that the region between the square with radius and the previously added layer contain.
bundle (bundle_tag) – Indicate that the current cell contains a burn bundle.
part (str) – Assembly part the current cell belongs to.

Note

If this is the first layer added, the region will be a complete square.

add_generic_shell(region, material=None, bundle=None, facility=None, part=None):

Adds an arbitrarily shaped shell.

Parameters:

region (region) – Outer boundary of the shell.
material (material) – Material that fills the region between the boundary region and the previously added layer.
facility (str) – Facility that the region between the boundary region and the previously added layer contain.
bundle (bundle_tag) – Indicate that the current cell contains a burn bundle.
part (str) – Assembly part the current cell belongs to.

Attention

The order in which layers are added is important, as the first layer will be used to construct the inner cell, the second layer adds the tubular cell between it and the previous cell etc. In particular, the radius parameter of subsequent calls must be increasing.

Note

The object also allows a short hand method using the | operator. See the examples below.

build(center=(0, 0), axial=primitives.Space(), radial=primitives.Space())

Create all the cells defined by this composite.

Parameters:

center (point) – Center (in the \(xy\)-plane) at which cells should be constructed.
axial (region) – Region that should be used to clip cells axially.
radial – Region that should be used to clip cells radially. This clip is only applied to the region outside the last layer.

Returns:

A list of cells.

Examples:

>>> from csg import *
>>> from core import *
>>> pc = composites.PinCell()
>>> pc.add_cylindrical_shell(0.8 * units.cm, bundle=bundel_tags.fuel, part='Fuel')
>>> from material_library.gasses import Helium
>>> pc.add_cylindrical_shell(0.9 * units.cm, material=Helium())
>>> from material_library.structural import Al6061
>>> pc.add_cylindrical_shell(1.2 * units.cm, material=Al6061(), part='Cladding')
>>> from material_library.moderators import BoratedWater
>>> pc.fill_material = BoratedWater()

Instead of using explicit function calls to add shells, the following equivalent short hand notation can also be used:

>>> pc = PinCell() | (0.8 * units.cm, bundle_tags.fuel) | (0.9 * units.cm, Helium()) | (1.2 * units.cm, Al6061())

To add a crate like structure, add the following to the above:

>>> pc.add_square_shell(3 * units.cm, material=LightWater())
>>> pc.fill_material = Al6061()

Adapters

class adapters.CylindricalAdapter

Models a sequence of axially aligned, stacked cylinders, with conic sections connecting them.

alignment
: integer
= ALIGNMENT.Z_AXIS

Coordinate axis along which the cylinders extend. Either ALIGNMENT.X_AXIS, ALIGNMENT.Y_AXIS or ALIGNMENT.Z_AXIS.

center
: point
= :math:`(0,0,0)`

Center of the part.

structural_material
: material
!

Material that should be used to fill structural cells.

fill_material
: material


Material that fills all other cells.

holes
: list (region)


Optional holes that should be punched into the structure of the adapter. This is either a single list of regions, whose interiors describe the hole (or cavity) structure that should be removed from all cylindrical segment, or a list equal to the number of cylindrical sections, giving the hole structures for each segment.

Attention

These regions are used exactly as specified, and are not transformed to any particular segment position.

Note

Holes are filled with the fill_material.

cylinder_diameters
: list (length)
!

List giving the outer diameter of each cylinder in the stack. The length of this list is also used to determine the total number of cylinders in the stack.

wall_thickness
: list (length)
!

List giving the thickness of each cylinder in the stack, that is, the inner diameter is the outer diameter, minus two times this value.

Must have the same number of entries as cylinder_diameters.

ranges
: list (tuple)
!

List giving the axial extent and positioning of each cylindrical segment in the stack. Entries are a tuple of two values giving the upper and lower bound respectively.

Must have the same number of entries as cylinder_diameters.

Attention

The intervals must be mutually disjoint, and either increasing or decreasing.

build(include_fill=True)

Creates all the cells making up the structure of the adapter.

Parameters:

include_fill (bool) – Add all non-structural cells. These cells fill be filled with material fill_material.

Returns:

A list of cells that can be added to your component.

Each segment will be connected with a conic structure, whose outer boundary connects the outer diameters of the two cylinders, and whose inner region connects the inner diameters. If there is no distance between subsequent sections, this connection piece will be missing.

Note

If fill_material is True, the method will also add all outside cells, effectively filling the entire strip between the minimum and maximum value in ranges.

Example:
>>> from csg import *
>>> from material_library.structural import Al6061
>>> adp = adapters.CylindricalAdapter()
>>> adp.cylinder_diameters = [2.0 * units.cm, 1.5 * units.cm]
>>> adp.wall_thickness = [2.0 * units.mm, 1.0 * units.mm]
>>> adp.ranges = [(0.0 * units.cm, -1.0 * units.cm), (-1.5 * units.cm, -3.0 * units.cm)]
>>> p = adp.build()
>>> p.show()
To punch a 4 mm hole through the top cylinder add:
>>> h = primitives.Cx(radius=2.0 * units.mm, center=(0, -0.5 * units.cm))
>>> adp.holes = [h, None]
>>> p = adp.build()
>>> p.show()

class adapters.RectangleToCircleAdapter

Part that connects a rectangular shaped base to a cylinder, with an (optional) cylindrical nozzle piece at the end.

alignment
: integer
= ALIGNMENT.Z_AXIS

Coordinate axis along which the adapter extend. Either ALIGNMENT.X_AXIS, ALIGNMENT.Y_AXIS or ALIGNMENT.Z_AXIS.

center
: point
= :math:`(0,0,0)`

Center of the part.

Note

The axial coordinate is taken as the center of the main part, that is, the rectangular box without the nozzle.

height
: length
!

Total height of the box adapter. This excludes the nozzle piece.

structural_material
: material
!

Material that should be used to fill structural cells.

fill_material
: material


Material that fills all other cells.

holes
: list (region)


Optional holes that should be punched into the structure of the adapter. This is a single list of regions, whose interiors describe the hole (or cavity) structure that should be removed from the main adapter part.

Note

Holes are filled with the fill_material.

outer_rectangle_width
: length
!

Total width (\(u\) dimension) of the adapter’s rectangular base.

outer_rectangle_length
: length
!

Total length (\(v\) dimension) of the adapter’s rectangular base.

inner_rectangle_width
: length
!

Width (\(u\) dimension) of the adapter’s rectangular base’s inner region. Must be less than outer_rectangle_width, so that the side wall thickness is the half the difference between these two values.

inner_rectangle_length
: length
!

Length (\(v\) dimension) of the adapter’s rectangular base’s inner region. Must be less than outer_rectangle_length, so that the side wall thickness is the half the difference between these two values.

target_outer_radius
: length
!

Outer radius of the target cylindrical structure.

target_inner_radius
: length
!

Inner radius of the target cylindrical structure.

nozzle_height
: length
= 0.0

Optional height of the cylindrical nozzle that will be added to the target cylindrical part of the adapter. If larger than 0.0, a cylindrical tube of this height with outer radius target_outer_radius and inner radius target_inner_radius will be added to this part.

nozzle_holes
: list (region)


List of holes that should be cut from the nozzle piece.

build(include_fill=True)

Creates all the cells making up the structure of the adapter.

Parameters:

include_fill (bool) – Add all non-structural cells. These cells fill be filled with material fill_material.

Returns:

A list of cells that can be added to your component.

Note

If fill_material is True, the method will also add all outside cells, effectively filling the entire strip between center.a - 0.5 * height and center.a + 0.5 * height + nozzle_height.

Examples:
>>> from csg import *
>>> adp = adapters.RectangleToCircleAdapter()
>>> adp.outer_rectangle_width = 0.5 * (73.2 + 73.05) * units.mm
>>> adp.outer_rectangle_length = 0.5 * (74.50 + 74.35) * units.mm
>>> adp.inner_rectangle_width = 63.5 * units.mm
>>> adp.inner_rectangle_length = 66 * units.mm
>>> adp.target_outer_radius = 0.5 * 70.55 * units.mm
>>> adp.target_inner_radius = 0.5 * 63.5 * units.mm
>>> adp.height = 63.5 * units.mm + 19.25 * units.mm
>>> from material_library.structural import Al6061
>>> from material_library.moderators import H2O
>>> adp.structural_material = Al6061()
>>> adp.fill_material = H2O()
>>> p = adp.build()
>>> p.show()
To add 5 cm of nozzle tubing:
>>> adp.nozzle_height = 5 * units.cm
>>> p = adp.build()
>>> p.show()