Penrose Kite and Dart Tilings with Haskell Diagrams
Infinite non-periodic tessellations of Roger Penrose’s kite and dart tiles.
As part of a collaboration with Stephen Huggett, working on some mathematical properties of Penrose tilings, I recognised the need for quick renderings of tilings. I thought Haskell diagrams would be helpful here, and that turned out to be an excellent choice. Two dimensional vectors were well-suited to describing tiling operations and these are included as part of the diagrams package.
This literate Haskell uses the Haskell diagrams package to draw tilings with kites and darts. It also implements the main operations of inflate and decompose which are essential for constructing tilings (explained below).
Firstly, these 5 lines are needed in Haskell to use the diagrams package:
{-# LANGUAGE NoMonomorphismRestriction #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE TypeFamilies #-}import Diagrams.Prelude
import Diagrams.Backend.SVG.CmdLineLegal tilings
These are the kite and dart tiles.
The red line marking here on the right hand copies, is purely to illustrate rules about how tiles can be put together for legal (non-periodic) tilings. Obviously edges can only be put together when they have the same length. If all the tiles are marked with red lines as illustrated on the right, the vertices where tiles meet must all have a red line or none must have a red line at that vertex. This prevents us from forming a simple rhombus by placing a kite top at the base of a dart and thus enabling periodic tilings.
All edges are powers of the golden section 
phi.
phi::Double
phi = (1.0 + sqrt 5.0) / 2.0So if the shorter edges are unit length, then the longer edges have length phi. We also have the interesting property of the golden section that 
All angles in the figures are multiples of tt which is 36 deg or 1/10 turn. We use ttangle to express such angles (e.g 180 degrees is ttangle 5).
ttangle:: Int -> Angle Double
ttangle n = (fromIntegral (n `mod` 10))*^tt
where tt = 1/10 @@ turnComponents
In order to implement inflate and decompose, we need to work with half tiles. These can be represented with a simple 2D vector to provide orientation and scale. We will call these components : Left Dart, Right Dart, Left Kite, Right Kite and use the type
data Component = LD (V2 Double)
| RD (V2 Double)
| LK (V2 Double)
| RK (V2 Double)The vector represents the join edge of each half tile where halves come together. The origin for a dart is the tip, and the origin for a kite is the acute angle tip (marked in the figure with a red dot).
These are the only 4 components we use (oriented along the x axis)
ldart,rdart,lkite,rkite:: Component
ldart = LD unitX
rdart = RD unitX
lkite = LK (phi*^unitX)
rkite = RK (phi*^unitX)
Perhaps confusingly we regard left and right of a dart differently from left and right of a kite when viewed from the origin. The diagram shows the left dart before the right dart and the left kite before the right kite. Thus in a complete tile, going clockwise round the origin the right dart comes before left dart, but the left kite comes before right kite.
When it comes to drawing components, for the simplest case, we just want to show the two tile edges of each component (and not the join edge). These edges are calculated as a list of 2 new vectors, using the join edge vector v. They are ordered clockwise from the origin of each component
compEdges:: Component -> [V2 Double]
compEdges (LD v) = [v',v ^-^ v'] where v' = phi*^rotate (ttangle 9) v
compEdges (RD v) = [v',v ^-^ v'] where v' = phi*^rotate (ttangle 1) v
compEdges (RK v) = [v',v ^-^ v'] where v' = rotate (ttangle 9) v
compEdges (LK v) = [v',v ^-^ v'] where v' = rotate (ttangle 1) vNow drawing lines for the 2 outer edges of a component is simply
drawComp:: Component -> Diagram B
drawComp = strokeLine . fromOffsets . compEdgesIt is also useful to calculate a list of the 4 tile edges of a completed half-tile component clockwise from the origin of the tile. (This is useful for colour filling a tile)
tileEdges:: Component -> [V2 Double]
tileEdges (LD v) =
compEdges (RD v) ++ map (zero ^-^) (reverse (compEdges (LD v)))
tileEdges (RD v) = tileEdges (LD v)
tileEdges (LK v) =
compEdges (LK v) ++ map (zero ^-^) (reverse (compEdges (RK v)))
tileEdges (RK v) = tileEdges (LK v)To fill whole tiles with colours, darts with dcol and kites with kcol we use fillDK. This uses only the left components to identify the whole tile and ignores right components so that a tile is not filled twice.
fillDK:: Colour Double -> Colour Double -> Component -> Diagram B
fillDK dcol kcol c = case c of
(LD _) -> (strokeLoop $ glueLine $ fromOffsets $ tileEdges c) # fc dcol
(LK _) -> (strokeLoop $ glueLine $ fromOffsets $ tileEdges c) # fc kcol
_ -> memptyTo show half tiles separately, we will use drawJComp which adds the join edge of each half tile to make loops with 3 edges.
drawJComp:: Component -> Diagram B
drawJComp = strokeLoop . closeLine . fromOffsets . compEdgesApart from drawing, we can also scale a component
scaleComp:: Double -> Component -> Component
scaleComp = compMap . scale where compMap applies a function to the vector of a component
compMap f (LD v) = LD (f v)
compMap f (RD v) = RD (f v)
compMap f (LK v) = LK (f v)
compMap f (RK v) = RK (f v)and rotate a component by an angle
rotateComp :: Angle Double -> Component -> Component
rotateComp = compMap . rotate (Positive rotations are in the anticlockwise direction.)
Patches
A patch is a list of components each with an offset represented by a vector
type Patch = [(V2 Double, Component)]To turn a whole patch into a diagram using some function cd for drawing the components, we use
patchWith cd patch = position $ fmap offset patch
where offset (v,c) = (origin .+^ v, cd c)Here the offset vector is used to calculate a point from the origin for each component to be positioned and drawn (with cd).
The common special case drawPatch uses drawComp on each component
drawPatch = patchWith drawCompOther operations on patches include scaling a patch
scalePatch :: Double -> Patch -> Patch
scalePatch r = fmap (\(v,c) -> (scale r v, scaleComp r c))and rotating a patch by an angle
rotatePatch :: Angle Double -> Patch -> Patch
rotatePatch a = fmap (\(v,c) -> (rotate a v, rotateComp a c))To move a patch by a vector, we simply add the vector to each offset vector
translatePatch:: V2 Double -> Patch -> Patch
translatePatch v0 = fmap (\(v,c) -> (v0 ^+^ v, c))As an aid to creating patches with 5-fold rotational symmetry, we combine 5 copies of a basic patch (rotated by multiples of ttangle 2 successively).
penta:: Patch -> Patch
penta p = concatMap copy [0..4]
where copy n = rotatePatch (ttangle (2*n)) pThis must be used with care to avoid nonsense patches. But two special cases are
sun = penta [(zero, rkite), (zero, lkite)]
star = penta [(zero, rdart), (zero, ldart)]This figure shows some example patches, drawn with drawPatch The first is a star and the second is a sun.
The tools so far for creating patches may seem limited (and do not help with ensuring legal tilings), but there is an even bigger problem.
Correct Tilings
Unfortunately, correct tilings – that is, tilings which can be extended to infinity – are not as simple as just legal tilings. It is not enough to have a legal tiling, because an apparent (legal) choice of placing one tile can have non-local consequences, causing a conflict with a choice made far away in a patch of tiles, resulting in a patch which cannot be extended. This suggests that constructing correct patches is far from trivial.
The infinite number of possible infinite tilings do have some remarkable properties. Any finite patch from one of them, will occur in all the others (infinitely many times) and within a relatively small radius of any point in an infinite tiling. (For details of this see links at the end)
This is why we need a different approach to constructing larger patches. There are two significant processes used for creating patches, namely inflate (also called compose) and decompose.
To understand these processes, take a look at the following figure.
Here the small components have been drawn in an unusual way. The edges have been drawn with dashed lines, but long edges of kites have been emphasised with a solid line and the join edges of darts marked with a red line. From this you may be able to make out a patch of larger scale kites and darts. This is an inflated patch arising from the smaller scale patch. Conversely, the larger kites and darts decompose to the smaller scale ones.
Decomposition
Since the rule for decomposition is uniquely determined, we can express it as a simple function on patches.
decompose :: Patch -> Patch
decompose = concatMap decompCwhere the function decompC acts on components and produces a list of the smaller components contained in the component. For example, a larger right dart will produce both a smaller right dart and a smaller left kite. Decomposing a component with offset v also takes care of the positioning, scale and rotation of the new components.
decompC (v, RD vd) = [(v, LK vd )
,(v ^+^ v', RD vd' )
] where v' = phi*^rotate (ttangle 1) vd
vd' = (2-phi) *^ (zero ^-^ v')
-- (2-phi) = 1/phi^2
decompC (v, LD vd) = [(v, RK vd )
,(v ^+^ v', LD vd' )
] where v' = phi*^rotate (ttangle 9) vd
vd' = (2-phi) *^ (zero ^-^ v')
-- (2-phi) = 1/phi^2
decompC (v, RK vk) = [(v, RD vd' )
,(v ^+^ v', LK vk' )
,(v ^+^ v', RK vk' )
] where v' = rotate (ttangle 9) vk
vd' = (2-phi) *^ v'
-- v'/phi^2
vk' = ((phi-1) *^ vk) ^-^ v'
-- (phi-1) = 1/phi
decompC (v, LK vk) = [(v, LD vd' )
,(v ^+^ v', RK vk' )
,(v ^+^ v', LK vk' )
] where v' = rotate (ttangle 1) vk
vd' = (2-phi) *^ v'
-- v'/phi^2
vk' = ((phi-1) *^ vk) ^-^ v'
-- (phi-1) = 1/phiThis is illustrated in the following figure for the cases of a right dart and a right kite.
The symmetric diagrams for left components are easy to work out from these, so they are not illustrated.
With the decompose operation we can start with a simple correct patch, and decompose repeatedly to get more and more detailed patches. (Each decomposition scales the tiles down by a factor of 
This figure illustrates how each component decomposes with 4 decomposition steps below each one.
components = [ldart, rdart, lkite, rkite]
fourDecomps = hsep 1 $ fmap decomps components # lw thin where
decomps cpt = vsep 1 $
fmap drawPatch $ take 5 $ decompositions [(zero,cpt)] We have made use of the fact that we can create an infinite list of finer and finer decompositions of any patch, using:
decompositions:: Patch -> [Patch]
decompositions p = inf where inf = p:fmap decompose infWe could get the n-fold decomposition of a patch as just the nth item in such a list
multiDecomp :: Int -> Patch -> Patch
multiDecomp n p = decompositions p !! nFor example, here is an infinite list of decomposed versions of sun.
suns = decompositions sunThe coloured tiling shown at the beginning is simply 6 decompositions of sun displayed using fillDK
sun6 = suns!!6
filledSun6 = patchWith (fillDK red blue) sun6 # lw ultraThinThe earlier figure illustrating larger kites and darts emphasised from the smaller ones is also sun6 but this time drawn with
experimentFig = patchWith experiment sun6 # lw thinwhere components are drawn with
experiment:: Component -> Diagram B
experiment c = emph c <> (drawComp c # dashingN [0.002,0.002] 0 # lw ultraThin)
where emph c = case c of
(LD v) -> (strokeLine . fromOffsets) [v] # lc red
-- emphasise join edge of darts in red
(RD v) -> (strokeLine . fromOffsets) [v] # lc red
(LK v) -> (strokeLine . fromOffsets) [rotate (ttangle 1) v]
-- emphasise long edge for kites
(RK v) -> (strokeLine . fromOffsets) [rotate (ttangle 9) v]Inflation
You might expect inflation (also called composition) to be a kind of inverse to decomposition, but it is a bit more complicated than that. With our current representation of components, we can only inflate single components. This amounts to embedding the component into a larger component that matches how the larger component decomposes. There is thus a choice at each inflation step as to which of several possibilities we select as the larger half-tile. We represent this choice as a list of alternatives. This list should not be confused with a patch. It only makes sense to select one of the alternatives giving a new single component.
The earlier diagram illustrating how decompositions are calculated also shows the two choices for inflating a right dart into either a right kite or a larger right dart. There will be two symmetric choices for a left dart, and three choices for left and right kites.
Once again we work with offsets so we can calculate the offset for any resulting component (to ensure the larger component contains the original in its original position in a decomposition).
inflate :: (V2 Double, Component) -> [(V2 Double, Component)]
inflate (v, RD vd) = [(v ^+^ v', RD vd')
,(v, RK vk )
] where v' = (phi+1) *^ vd
-- vd*phi^2
vd' = rotate (ttangle 9) (vd ^-^ v')
vk = rotate (ttangle 1) v'
inflate (v, LD vd) = [(v ^+^ v', LD vd')
,(v, LK vk )
] where v' = (phi+1) *^ vd
-- vd*phi^2
vd' = rotate (ttangle 1) (vd ^-^ v')
vk = rotate (ttangle 9) v'
inflate (v, RK vk) = [(v, LD vk )
,(v ^+^ lv', LK lvk')
,(v ^+^ rv', RK rvk')
] where lv' = phi*^rotate (ttangle 9) vk
rv' = phi*^rotate (ttangle 1) vk
rvk' = phi*^rotate (ttangle 7) vk
lvk' = phi*^rotate (ttangle 3) vk
inflate (v, LK vk) = [(v, RD vk )
,(v ^+^ rv', RK rvk')
,(v ^+^ lv', LK lvk')
] where v0 = rotate (ttangle 1) vk
lv' = phi*^rotate (ttangle 9) vk
rv' = phi*^rotate (ttangle 1) vk
rvk' = phi*^rotate (ttangle 7) vk
lvk' = phi*^rotate (ttangle 3) vkAs the result is a list of alternatives, we need to select one to do further inflations. We can express all the alternatives after n steps as inflations n where
inflations :: Int -> (V2 Double, Component) -> [(V2 Double, Component)]
inflations 0 vc = [vc]
inflations n vc = [vc'' | vc' <- inflate vc, vc'' <- inflations (n-1) vc'](That last line could be written more succinctly with do notation)
This figure illustrates 5 consecutive choices for inflating a left dart to produce a left kite. On the left, the finishing component is shown with the starting component embedded, and on the right the 5-fold decomposition of the result is shown.
fiveInflate = hsep 1 $ fmap drawPatch [[ld,lk'], multiDecomp 5 [lk']]
-- two seperate patches
where
ld = (zero,ldart)
lk = inflate ld !!1
rk = inflate lk !!1
rk' = inflate rk !!2
ld' = inflate rk' !!0
lk' = inflate ld' !!1Finally, at the end of this literate haskell program we choose which figure to draw as output.
fig::Diagram B
fig = filledSun6main = mainWith figThat’s it. But, What about inflating whole patches?, I hear you ask. Unfortunately we need to answer questions like what components are adjacent to a component in a patch and whether there is a corresponding other half for a component. These cannot be done with our simple vector representations. We would need some form of planar graph representation, which is much more involved. That is another story.
Many thanks to Stephen Huggett for his inspirations concerning the tilings. A library version of the above code is available on GitHub
Further reading on Penrose Tilings
As well as the Wikipedia entry Penrose Tilings I recommend two articles in Scientific American from 2005 by David Austin Penrose Tiles Talk Across Miles and Penrose Tilings Tied up in Ribbons.
There is also a very interesting article by Roger Penrose himself: Penrose R Tilings and quasi-crystals; a non-local growth problem? in Aperiodicity and Order 2, edited by Jarich M, Academic Press, 1989.
More information about the diagrams package can be found from the home page Haskell diagrams
is a linear operator (such as the laplace operator 
or higher order or lower order partial derivatives or any linear combination of these), 
is given and 
is the function we want to solve for, and both are defined on a region 
. The multigrid scheme starts with a fine grid 
(where 
is the grid spacing) on which we want to obtain an approximate solution by finite difference methods. The scheme involves the use of a coarser grid (e.g. 
) and, recursively, a stack of coarser and coarser grids to apply error correction on approximations. The full multigrid scheme also uses coarser grids to improve initial approximations used by multigrid.
on the fine grid (
.
are called the residuals (
) and these we can calculate.
,
and get a corrected guess
to get a final approximation.
we want a (linear) mapping that is full rank. E.g identity plus averaging neighbours.
or take some weighting of neighbours as well. A common requirement is that this is also linear and a constant times the transpose of the interpolation when these are expressed as matrix operations.
for some function 
on the boundary 
. We will take the standard approach of representing the boundary conditions using a mask array (with 0’s indicating boundary positions and 1’s at non-boundary positions) along with a boundary values array which has 0’s at non-boundary positions. This will enable us to adapt these two arrays for different coarsenesses of grid.
. We can thus assume two dimensional arrays represent 
, Poisson’s equation has the form
subject to the Dirichlet boundary condition 
for 
, then we divide by 4 and reset the boundary with the mask and values.
, we can use the same formula as above to approximate 
. Rearranging we have 
by 
between 0 and 1. This slightly increases runtime, but does not improve errors for 
, where 
is a scalar potential (such as temperature) over a smooth surface.
in the 
direction and 
in the 
direction). Approximating solutions numerically amounts to approximating a fixedpoint for the grid values satisfying the boundary values and also, for non-boundary points (i,j) satisfying
then this simplifies to
from previous values 
until the iterations are close enough to converging.
to 
for calculated values in the previous column and in the previouus row).
as the weight we calculate first 
(substituting 
