---
jupyter:
  description: 'Vaults and Steinmetz solids with planex2d and OpenCASCADE: two barrel vaults over a square bay, united
    into a groin vault and intersected into a cloister vault, the solids of two and three crossed cylinders against
    16r³/3 and 8(2 - √2)r³, and the sphere inside, whose sections are pi/4 of the squares.'
---

# 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.

![A groin vault, two barrel vaults united over a square bay](../_static/gallery/vaults.png){.result}

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

```python nbsphinx="hidden"
from docs_figures import Row
```

## 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$:

```python
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)
```

```python nbsphinx="hide-input"
Row((barrel, "A barrel vault along x", {"views": [(barrel, "shade"), (barrel, "ink")], "angles": (-60, 0, 30)}),
    (across, "Its copy along y", {"views": [(across, "shade"), (across, "ink")], "angles": (-60, 0, 30)}),
    caption="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$:

```python
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)
```

```python nbsphinx="hide-input"
Row((groin, "The groin vault, the union", {"views": [(groin, "shade"), (groin, "ink")], "angles": (-60, 0, 30)}),
    (cloister, "The cloister vault, the intersection",
     {"views": [(cloister, "shade"), (cloister, "ink")], "angles": (-60, 0, 30)}),
    caption="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$:

```python
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)
```

## 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:

```python
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)
```

```python nbsphinx="hide-input"
Row((two, "Two cylinders", {"views": [(two, "shade"), (two, "ink")], "angles": (-60, 0, 30)}),
    (three, "Three cylinders", {"views": [(three, "shade"), (three, "ink")], "angles": (-60, 0, 30)}),
    caption="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](slicer.md) 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:

```python
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)
```

```python nbsphinx="hide-input"
Row(*[(section(two, z), f"At z = {z}", {"fills": [(section(ball, z), "#d98c2b")], "view": ((-r, r), (-r, r))})
      for z in (0.0, 2.0, 3.0)],
    caption="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$:

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

```python nbsphinx="hidden"
import sys
sys.path.insert(0, "..")
from docs_figures import ViewFigure
from gallery import save_figure
fig, aspect = ViewFigure([(groin, "shade"), (groin, "ink")], (-60, 0, 30))
save_figure(fig, "vaults", dpi=160)
```