---
jupyter:
  description: A heater track along a Hilbert curve, built with planex2d and checked against its exact area; then
    the current, the Joule heat and the temperature of the plate.
---

# Hilbert curve heater track

A heater track along a Hilbert curve fills its square evenly and therefore heats it evenly. The example has two
parts. The first builds the track of order 7, with 16384 points on a `pitch` = 2 mm grid and `width` = 1 mm wide. The
second puts a heater of order 4 on an aluminum plate and computes the current through the track, its Joule heat and
the plate's temperature.

![The temperature of an aluminum plate heated by a Hilbert curve heater track](../_static/gallery/hilbert-heater.png){.result}

```python
import math
import numpy as np
import ngsolve as ngs
from ngsolve_webgpu.jupyter import Draw
import planex as px
```

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

`hilbert(order)` gives the points of the curve, on the grid of the pitch:

```python
pitch = 2.0
width = 1.0


def hilbert(order):
    n, points = 2**order, []
    for d in range(n * n):
        x = y = 0
        s, t = 1, d
        while s < n:
            rx = 1 & (t // 2)
            ry = 1 & (t ^ rx)
            if ry == 0:
                if rx == 1:
                    x, y = s - 1 - x, s - 1 - y
                x, y = y, x
            x, y = x + s * rx, y + s * ry
            t //= 4
            s *= 2
        points.append((pitch * x, pitch * y))
    return points


points = hilbert(7)
print(len(points), "points")
```

## The track

The path is a sketch of lines through the points, and the track is its wire offset by half the width:

```python
path = px.QuickSketch(points[0])
for p in points[1:]:
    path.Line(p)
track = path.Wire().Offset(width / 2)
```

**Check:** the track's area is the width times the path length, plus the two round ends. At every 90° turn, the
square (w/2)² inside the corner is counted twice, and the quarter disc outside it is added.

```python
turns = sum(1 for p, q, r in zip(points, points[1:], points[2:])
            if (q[0] - p[0]) * (r[1] - q[1]) != (q[1] - p[1]) * (r[0] - q[0]))
length = pitch * (len(points) - 1)
exact = width * length + math.pi * width**2 / 4 + turns * (width / 2) ** 2 * (math.pi / 4 - 1)
print(turns, "turns")
print("area", track.mass, "exact", exact)
```

## A heater on a plate

The track of order 7 is 33 m long, far too long to mesh here. The heater is the same curve of order 4 on the same
grid, with 256 points in a 30 mm square, made of constantan foil 25 µm thick, on a 40 mm aluminum plate 1.5 mm
thick. Leads 3 long take both ends to the plate's lower edge, 2 below the heater, whose two halves are named as the
contacts. The plate reaches 5 beyond the heater on the other sides. Cutting the track to the plate gives its ends
those names, and gluing it into the plate gives one mesh for both (Figure 1).

```python
heater_points = hilbert(4)
side = heater_points[-1][0]
path = px.QuickSketch((0, -3)).Line(heater_points[0])
for p in heater_points[1:]:
    path.Line(p)
path.Line((side, -3))
plate = px.QuickSketch((-5, -2)).Line((side / 2, -2), name="minus").Line((side + 5, -2), name="plus")
plate = plate.Line((side + 5, side + 5)).Line((-5, side + 5)).Close().Face()
plate.name = "plate"
heater = path.Wire().Offset(width / 2) * plate
heater.name = "track"
heater.vertices.hpref = 1
board = px.Glue([plate, heater])
```

```python nbsphinx="hide-input"
Row((board, "The plate and the track of order 4"), caption="A heater of constantan foil on an aluminum plate")
```

## Current

In the foil the potential solves Laplace's equation with the supply voltage on one contact and zero on the other. The
current comes from the leftover equations at the fixed contact, so the electric power V·I must equal the Joule heat
integrated over the track. The resistance in squares (the resistance over the foil's sheet resistance ρ/t) is less than
the center line's length over the width, because the current cuts across at a corner. A sharp 90° corner counts as about
0.56 of a square instead of one, and these corners have a rounded outside. The inside corners of the turns are
re-entrant, so the current density is singular there, and the mesh is refined towards them, with the potential of order
3. The color scale ends at 1.2 times the mean field along the track, and at the inside corners the field goes far above
it.

```python
rho = 4.9e-4
foil = 0.025
volts = 5.0
mesh = ngs.Mesh(px.OCCGeometry(board).GenerateMesh(maxh=0.5))
mesh.RefineHP(levels=2)
mesh.Curve(3)
E = ngs.H1(mesh, order=3, definedon="track", dirichlet="plus|minus")
phi, psi = E.TnT()
conduct = ngs.BilinearForm(foil / rho * ngs.grad(phi) * ngs.grad(psi) * ngs.dx("track", bonus_intorder=4)).Assemble()
potential = ngs.GridFunction(E)
potential.Set(volts, definedon=mesh.Boundaries("plus"))
potential.vec.data += conduct.mat.Inverse(E.FreeDofs()) * (-conduct.mat * potential.vec)

reaction = potential.vec.CreateVector()
reaction.data = conduct.mat * potential.vec
contact = ngs.GridFunction(E)
contact.Set(1, definedon=mesh.Boundaries("plus"))
current = ngs.InnerProduct(contact.vec, reaction)
joule = foil / rho * ngs.grad(potential) * ngs.grad(potential)
center_line = [(0, -2)] + heater_points + [(side, -2)]
corners = sum(1 for p, q, r in zip(center_line, center_line[1:], center_line[2:])
              if (q[0] - p[0]) * (r[1] - q[1]) != (q[1] - p[1]) * (r[0] - q[0]))
strip = sum(math.dist(p, q) for p, q in zip(center_line, center_line[1:])) / width
squares = volts / current * foil / rho
heat = ngs.Integrate(joule, mesh, definedon=mesh.Materials("track"), order=10)
print(mesh.ne, "elements")
print("current", round(current, 4), "A, resistance", round(volts / current, 3), "Ohm")
print("power V I", volts * current, "W, Joule heat", heat, "W")
print(round(squares, 1), "squares instead of", round(strip), "along the center line")
print("each of the", corners, "corners counts", round(1 - (strip - squares) / corners, 2), "of a square")
Draw(-ngs.grad(potential), mesh, min=0, max=1.2 * volts / (squares * width), autoscale=False)
```

## Temperature

The plate conducts the heat sideways and gives it to still air on both faces. With the air at 20 °C, the heat the plate
gives off must equal the electric power; the refined elements at the corners are not affine, so the heat source is
integrated with the same rule as the Joule heat above. The temperature over the heated square shows how even the heater
is.

```python
k = 0.237
plate_t = 1.5
h = 1e-5
T_air = 20.0
T = ngs.H1(mesh, order=2)
theta, v = T.TnT()
plate_form = k * plate_t * ngs.grad(theta) * ngs.grad(v) + 2 * h * theta * v
heat = ngs.BilinearForm(plate_form * ngs.dx(bonus_intorder=4)).Assemble()
rules = {ngs.ET.TRIG: ngs.IntegrationRule(ngs.TRIG, 10), ngs.ET.QUAD: ngs.IntegrationRule(ngs.QUAD, 10)}
source = ngs.LinearForm(joule * v * ngs.dx("track", intrules=rules)).Assemble()
rise = ngs.GridFunction(T)
rise.vec.data = heat.mat.Inverse(T.FreeDofs()) * source.vec

print("heat given off", ngs.Integrate(2 * h * rise, mesh, order=10), "W, electric power", volts * current, "W")
points3 = {ngs.ET.TRIG: ngs.IntegrationRule(ngs.TRIG, 3), ngs.ET.QUAD: ngs.IntegrationRule(ngs.QUAD, 3)}
samples = mesh.MapToAllElements(points3, ngs.VOL)
xy = ngs.CF((ngs.x, ngs.y))(samples)
inside = (xy[:, 0] > 0) & (xy[:, 0] < side) & (xy[:, 1] > 0) & (xy[:, 1] < side)
temperature = T_air + np.asarray(rise(samples)).ravel()[inside]
print("heated square from", round(temperature.min(), 2), "to", round(temperature.max(), 2), "°C")
Draw(T_air + rise, mesh)
```

```python nbsphinx="hidden"
from ngsolve.webgui import Draw as WebguiDraw
import sys
sys.path.insert(0, "..")
from gallery import save_scene, save_figure, HIDDEN, CLEAN
save_scene(WebguiDraw(T_air + rise, mesh, settings={**CLEAN, "Misc": {"subdivision": 2}, "Colormap": {"ncolors": 32}}, show=False), "hilbert-heater")
save_scene(WebguiDraw(ngs.Norm(ngs.grad(potential)), mesh, min=0, max=1.2 * volts / (squares * width),
                      settings={**CLEAN, "Misc": {"subdivision": 2}, "Colormap": {"ncolors": 32}}, show=False),
           "hilbert-current")
```