Vaults and Steinmetz solids

A barrel vault is a half cylinder over a long hall. Where two halls of the same span cross, two barrel vaults meet over a square bay, and a builder can join them in two ways. A groin vault keeps the space of both barrels, their union, and rests on the bay’s four corners. A cloister vault keeps only the space they have in common, their intersection, and rests on the bay’s four walls.

Two whole cylinders of radius \(r\) crossed at a right angle intersect in a Steinmetz solid. It fills two thirds of the cube around it, \(16r^3/3\). Its surface is \(16r^2\), and three cylinders give \(8(2 - \sqrt{2})\,r^3\) and \(24(2 - \sqrt{2})\,r^2\). Cut at any height, the solid is a square, and the sphere of radius \(r\) inside it a circle inscribed in that square, so the sphere has \(\pi/4\) of the solid’s volume. The example builds the vaults and the solids from sketches and checks each of these numbers.

import math
import netgen.occ as occ
import planex as px

The barrel vault

The bay is 8 m square, so the barrels have a radius of \(r = 4\,\text{m}\). A barrel’s section is the arch, a half disc on the plane \(x = -r\), drawn from \((r, 0)\) up along a half circle and closed along the floor. Padded by \(2r\), it spans the bay along \(x\). Rotate of netgen.occ gives a copy turned by 90° about the \(z\) axis, which spans it along \(y\) (Figure 1). Each barrel holds \(\pi r^2/2 \cdot 2r = \pi r^3\):

r = 4.0
arch = px.QuickSketch(px.Plane(px.YZ, -r)).Move((r, 0), dir=(0, 1)).Arc(r, 180).Close()  # from (r, 0) upward
barrel = px.Pad(arch, 2 * r)
across = barrel.Rotate(occ.Axis((0, 0, 0), occ.Z), 90)
print("barrel", barrel.mass, "m³, exact", math.pi * r**3)

barrel 201.06192982974682 m³, exact 201.06192982974676

(a) A barrel vault along x

(b) Its copy along y

Figure 1: Two barrel vaults over one bay

Groin and cloister vaults

The groin vault is the union of the two barrels, the cloister vault their intersection (Figure 2). The cloister vault is the upper half of the Steinmetz solid, \(8r^3/3\), and the groin vault holds the two barrels less what they share, \(2\pi r^3 - 8r^3/3\):

groin = barrel + across
cloister = barrel * across
print("groin vault", groin.mass, "m³, exact", 2 * math.pi * r**3 - 8 * r**3 / 3)
print("cloister vault", cloister.mass, "m³, exact", 8 * r**3 / 3)

groin vault 231.45719299304196 m³, exact 231.45719299282686
cloister vault 170.6666666662947 m³, exact 170.66666666666666

(a) The groin vault, the union

(b) The cloister vault, the intersection

Figure 2: The two ways to join two barrel vaults

The two barrels’ surfaces meet along the diagonals of the bay. A point on the diagonal \(x = y\) lies in one barrel exactly when it lies in the other, so both vaults have the same section there, a half ellipse with the semi-axes \(r\sqrt{2}\) along the diagonal and \(r\) up. In the groin vault, these curves are the groins, the edges that give it its name. In the cloister vault, they are the inner corners where its four panels meet. Both sections have the area \(\pi r\sqrt{2} \cdot r / 2\):

diagonal = px.QuickSketch(px.Plane(origin=(0, 0, 0), normal=(1, -1, 0), axis=(1, 1, 0))).RectangleC(4 * r, 4 * r)
for name, vault in (("groin", groin), ("cloister", cloister)):
    cut = vault * diagonal.Face().ToOCCFace()
    print(name, "vault on the diagonal", sum(f.mass for f in cut.faces), "m², exact", math.pi * r * math.sqrt(2) * r / 2)

groin vault on the diagonal 35.543063505266936 m², exact 35.54306350526693
cloister vault on the diagonal 35.54306350526694 m², exact 35.54306350526693

Steinmetz solids

A whole disc on the same plane, padded by \(2r\), is a cylinder along \(x\). Its copies turned about \(z\) and about \(y\) lie along \(y\) and \(z\). Two of them intersect in the bicylinder, all three in the tricylinder (Figure 3). The intersection of two cylinders leaves each of the bicylinder’s curved faces in two pieces, joined along the square at \(z = 0\) where the surface runs on smoothly. UnifySameDomain of netgen.occ merges faces that lie on one surface into one. The volumes and surfaces are compared with the formulas:

cylinder = px.Pad(px.QuickSketch(px.Plane(px.YZ, -r)).Circle(r), 2 * r)
along_y = cylinder.Rotate(occ.Axis((0, 0, 0), occ.Z), 90)
along_z = cylinder.Rotate(occ.Axis((0, 0, 0), occ.Y), 90)
two = (cylinder * along_y).UnifySameDomain()  # four curved faces, not eight
three = two * along_z


def surface(solid):
    return sum(f.mass for f in solid.faces)


print("bicylinder", two.mass, "m³, exact", 16 * r**3 / 3, "| surface", surface(two), "m², exact", 16 * r**2)
k = 2 - math.sqrt(2)
print("tricylinder", three.mass, "m³, exact", 8 * k * r**3, "| surface", surface(three), "m², exact", 24 * k * r**2)

bicylinder 341.3333333334597 m³, exact 341.3333333333333 | surface 256.0000000000948 m², exact 256.0
tricylinder 299.9226563858374 m³, exact 299.9226560649753 | surface 224.94199228937805 m², exact 224.94199204873146

(a) Two cylinders

(b) Three cylinders

Figure 3: The Steinmetz solids of two and three cylinders

The sphere inside

At the height \(z\), both cylinders are a band of half width \(\sqrt{r^2 - z^2}\), so the bicylinder’s section is a square of side \(2\sqrt{r^2 - z^2}\). The sphere’s section is the circle of radius \(\sqrt{r^2 - z^2}\) inscribed in it. section cuts a solid with a large face on the plane at \(z\), as on the slicer page, and px.Project turns each cut face into a planex2d face (Figure 4). At every height, the circle has \(\pi/4\) of the square’s area:

ball = occ.Sphere(occ.Pnt(0, 0, 0), r)


def section(solid, z):
    plane = px.QuickSketch(px.Plane(px.XY, z)).RectangleC(4 * r, 4 * r).Face().ToOCCFace()
    return px.Unify([px.Face(px.Project(f)) for f in (solid * plane).faces])


for z in (0.0, 1.0, 2.0, 3.0, 3.9):
    square, circle = section(two, z), section(ball, z)
    print("z", z, "square", square.mass, "exact", 4 * (r**2 - z**2), "| circle / square", circle.mass / square.mass)
print("pi / 4", math.pi / 4)

z 0.0 square 64.0 exact 64.0 | circle / square 0.7853981633974483
z 1.0 square 60.00000000000001 exact 60.0 | circle / square 0.7853981633974483
z 2.0 square 47.99999999999999 exact 48.0 | circle / square 0.7853981633974483
z 3.0 square 28.000000000000007 exact 28.0 | circle / square 0.7853981633974476
z 3.9 square 3.1600000000000037 exact 3.1600000000000037 | circle / square 0.7853981633974514
pi / 4 0.7853981633974483

(a) At z = 0.0

(b) At z = 2.0

(c) At z = 3.0

Figure 4: Sections of the bicylinder, the sphere's circle shaded

Since every section of the sphere has \(\pi/4\) of the area of the bicylinder’s section at the same height, the volumes are in the same ratio, by Cavalieri’s principle. So \(\pi/4 \cdot 16r^3/3\) is the volume of the sphere, \(4\pi r^3/3\):

print("sphere", ball.mass, "m³ | pi / 4 of the bicylinder", math.pi / 4 * two.mass, "| exact", 4 / 3 * math.pi * r**3)

sphere 268.08257310632894 m³ | pi / 4 of the bicylinder 268.08257310642824 | exact 268.082573106329