To follow on from the recent series on using hyperbolic tessellation to generate patterns that might be used for 3D printing, I decided to research a slightly different approach.
While I found hyperbolic tessellation reasonably straightforward for generating 2D patterns, it was much harder to adapt to 3D, mainly because we’d need to create irregular polyhedra rather than just irregular 2D polygons. I had enough trouble getting my head around the 2D side of things that I decided not to go down that path for now, at least. The good news from Alex Fielder is that he’s been making good progress getting the hyperbolic tessellation code working inside Inventor (very cool), focusing – as far as I know – on using the 2D pattern to generate ribs rather than polyhedra (a very pragmatic approach).
I do think the type of pattern created through hyperbolic tessellation is probably more interesting from a structural perspective (having smaller versions at the edges with larger ones at the middle just makes a lot of sense), especially if working with polyhedra rather than ribs.
I expect to come back to something similar, in due course, but for now have been looking into a different approach: thinking about more regular (and therefore simpler) geometric forms by looking at circle packing before moving on to the very much related sphere packing.
There are a number of interesting strategies for circle packing: the one I like best is called the Apollonian gasket (which has been covered on a number of different Wikipedia pages, as well as on Wolfram MathWorld), which is also known as Leibniz packing. A related, specific case is that of Ford circles. The idea is that after creating three initial circles (each tangent to each other and to the outer circle), you draw the largest possible circle in the remaining space, continuing to sub-divide fractally as far as you’d like.
[As an aside… while researching this topic, I came across this interesting site by Takaya Iwamoto, a former colleague at Autodesk, who clearly shares my interest in fractals.]
And this of course also applies to 3D, although it may require a slightly different approach to get there.
I really like this pattern: while it may not create an internal structure optimised for strength, it should still be pretty interesting from a material savings perspective (short of having a completely hollow shell, of course). And I fully expect to looking into more structural filling techniques, further down the line.
The main problem to solve when generating the above pattern is to find a fourth mutually tangent circle (in 2D) or a fifth mutually tangent sphere (in 3D). It’s possible to use Descartes’ theorem in 2D or the Soddy-Gosset theorem in 3D (and beyond, for the truly fearless).
My mathematics being as they are (i.e rusty at best), I did my best to find someone else’s code to adapt for the various posts in this series. I started with a C# implementation (it also has a VB.NET sibling), and reworked it to generate AutoCAD geometry.
using Autodesk.AutoCAD.ApplicationServices;
using Autodesk.AutoCAD.DatabaseServices;
using Autodesk.AutoCAD.EditorInput;
using Autodesk.AutoCAD.Geometry;
using Autodesk.AutoCAD.Runtime;
using System.Collections.Generic;
using System;
namespace CirclePacking
{
public class Commands
{
static double _minRad = 0.01;
static short _color = 5;
public class MyCircle
{
public Point2d Center;
public double Radius;
public short ColorIndex;
}
[CommandMethod("AGC")]
public static void ApollonianGasket()
{
Document doc =
Application.DocumentManager.MdiActiveDocument;
Database db = doc.Database;
Editor ed = doc.Editor;
PromptEntityOptions peo =
new PromptEntityOptions("\nSelect circle: ");
peo.AllowNone = true;
peo.SetRejectMessage("\nMust be a circle.");
peo.AddAllowedClass(typeof(Circle), false);
PromptEntityResult per = ed.GetEntity(peo);
if (per.Status != PromptStatus.OK &&
per.Status != PromptStatus.None
)
return;
Transaction tr =
doc.TransactionManager.StartTransaction();
using (tr)
{
BlockTableRecord btr =
(BlockTableRecord)tr.GetObject(
db.CurrentSpaceId, OpenMode.ForWrite
);
// If the user selected a circle, use it
Circle cir;
if (per.Status != PromptStatus.None && per.ObjectId != null)
{
cir =
tr.GetObject(per.ObjectId, OpenMode.ForRead) as Circle;
}
else
{
// Otherwise create a new one at the default location
cir = new Circle(Point3d.Origin, Vector3d.ZAxis, 500);
btr.AppendEntity(cir);
tr.AddNewlyCreatedDBObject(cir, true);
}
List<MyCircle> circles =
FindApollonianPacking(cir.Radius * 2);
// Offset each of the returned circles to fit our outer one
Vector3d offset =
cir.Center - Point3d.Origin -
new Vector3d(cir.Radius, cir.Radius, 0);
// Add each of the returned circles to the current space
foreach (MyCircle myc in circles)
{
Point3d cen =
new Point3d(myc.Center.X, myc.Center.Y, 0) + offset;
Circle c = new Circle(cen, Vector3d.ZAxis, myc.Radius);
c.ColorIndex = myc.ColorIndex;
btr.AppendEntity(c);
tr.AddNewlyCreatedDBObject(c, true);
}
tr.Commit();
}
}
// Find the Apollonian packing
private static List<MyCircle> FindApollonianPacking(
double width
)
{
List<MyCircle> circles = new List<MyCircle>();
// Create the three central tangent circles
double radius = width * 0.232f,
x = width / 2,
gaskHeight = 2 * (radius + 2 * radius / Math.Sqrt(3)),
y = (width - gaskHeight) / 2 + radius;
MyCircle cir0 = CreateCircle(x, y, radius);
circles.Add(cir0);
x -= radius;
y += radius * Math.Sqrt(3);
MyCircle cir1 = CreateCircle(x, y, radius);
circles.Add(cir1);
x += 2 * radius;
MyCircle cir2 = CreateCircle(x, y, radius);
circles.Add(cir2);
// Find the circle that contains them all
MyCircle outer =
FindApollonianCircle(cir0, cir1, cir2, -1, -1, -1);
circles.Add(outer);
// Set level to smaller values such as 3 to see partially
// drawn gaskets
int level = 10000;
// Find the central circle
FindCircleOutsideAll(circles, level, cir0, cir1, cir2);
// Find circles tangent to the big circle
FindCircleOutsideTwo(circles, level, cir0, cir1, outer);
FindCircleOutsideTwo(circles, level, cir1, cir2, outer);
FindCircleOutsideTwo(circles, level, cir2, cir0, outer);
return circles;
}
// Draw a circle tangent to these three circles and that is
// outside all three
private static void FindCircleOutsideAll(
List<MyCircle> circles,
int level, MyCircle cir0, MyCircle cir1, MyCircle cir2
)
{
MyCircle newCir =
FindApollonianCircle(cir0, cir1, cir2, 1, 1, 1);
if (newCir == null || newCir.Radius < _minRad)
return;
circles.Add(newCir);
if (--level > 0)
{
FindCircleOutsideAll(circles, level, cir0, cir1, newCir);
FindCircleOutsideAll(circles, level, cir0, cir2, newCir);
FindCircleOutsideAll(circles, level, cir1, cir2, newCir);
}
}
// Draw a circle tangent to these three circles and that is
// outside two of them
private static void FindCircleOutsideTwo(
List<MyCircle> circles, int level,
MyCircle cir0, MyCircle cir1, MyCircle contains
)
{
MyCircle newCir =
FindApollonianCircle(cir0, cir1, contains, 1, 1, -1);
if (newCir == null || newCir.Radius < _minRad)
return;
circles.Add(newCir);
if (--level > 0)
{
FindCircleOutsideTwo(circles, level, newCir, cir0, contains);
FindCircleOutsideTwo(circles, level, newCir, cir1, contains);
FindCircleOutsideAll(circles, level, cir0, cir1, newCir);
}
}
private static MyCircle FindApollonianCircle(
MyCircle c1, MyCircle c2, MyCircle c3, int s1, int s2, int s3
)
{
// Make sure c2 doesn't have the same X or Y coordinate as
// the others
double tiny = Tolerance.Global.EqualPoint;
if (
Math.Abs(c2.Center.X - c1.Center.X) < tiny ||
Math.Abs(c2.Center.Y - c1.Center.Y) < tiny
)
{
MyCircle tempCir = c2;
c2 = c3;
c3 = tempCir;
int temp_s = s2;
s2 = s3;
s3 = temp_s;
}
if (
Math.Abs(c2.Center.X - c3.Center.X) < tiny ||
Math.Abs(c2.Center.Y - c3.Center.Y) < tiny
)
{
MyCircle tempCir = c2;
c2 = c1;
c1 = tempCir;
int temp_s = s2;
s2 = s1;
s1 = temp_s;
}
double x1 = c1.Center.X,
y1 = c1.Center.Y,
r1 = c1.Radius,
x2 = c2.Center.X,
y2 = c2.Center.Y,
r2 = c2.Radius,
x3 = c3.Center.X,
y3 = c3.Center.Y,
r3 = c3.Radius;
double v11 = 2 * x2 - 2 * x1,
v12 = 2 * y2 - 2 * y1,
v13 =
x1 * x1 - x2 * x2 + y1 * y1 -
y2 * y2 - r1 * r1 + r2 * r2,
v14 = 2 * s2 * r2 - 2 * s1 * r1;
double v21 = 2 * x3 - 2 * x2,
v22 = 2 * y3 - 2 * y2,
v23 =
x2 * x2 - x3 * x3 + y2 * y2 - y3 * y3 -
r2 * r2 + r3 * r3,
v24 = 2 * s3 * r3 - 2 * s2 * r2;
double w12 = v12 / v11,
w13 = v13 / v11,
w14 = v14 / v11;
double w22 = v22 / v21 - w12,
w23 = v23 / v21 - w13,
w24 = v24 / v21 - w14;
double P = -w23 / w22,
Q = w24 / w22,
M = -w12 * P - w13,
N = w14 - w12 * Q;
double a = N * N + Q * Q - 1,
b =
2 * M * N - 2 * N * x1 + 2 * P * Q -
2 * Q * y1 + 2 * s1 * r1,
c =
x1 * x1 + M * M - 2 * M * x1 + P * P +
y1 * y1 - 2 * P * y1 - r1 * r1;
// Find roots of a quadratic equation
double[] solutions = QuadraticSolutions(a, b, c);
if (solutions.Length < 1)
return null;
double rs = solutions[0],
xs = M + N * rs,
ys = P + Q * rs;
if (
Math.Abs(xs) < tiny ||
Math.Abs(ys) < tiny ||
Math.Abs(rs) < tiny
)
return null;
return CreateCircle(xs, ys, Math.Abs(rs));
}
// Return solutions to a quadratic equation
private static double[] QuadraticSolutions(
double a, double b, double c
)
{
double discriminant = b * b - 4 * a * c;
// See if there are no real solutions
if (discriminant < 0)
return new double[] { };
// See if there is one solution
if (discriminant < Tolerance.Global.EqualPoint)
return new double[] { -b / (2 * a) };
// There are two solutions
return new double[]
{
(-b + Math.Sqrt(discriminant)) / (2 * a),
(-b - Math.Sqrt(discriminant)) / (2 * a),
};
}
private static MyCircle CreateCircle(
double x, double y, double radius
)
{
MyCircle cir = new MyCircle();
cir.Center = new Point2d(x, y);
cir.Radius = radius;
cir.ColorIndex = _color;
return cir;
}
}
}
The AGC command allows you to select an existing circle to fill or it creates one at the origin and fills that. One thing to bear in mind: with this implementation the size of the circle does matter: as the fractal stops at circles of radius 0.01, an outer circle of 500 radius is going to lead to more geometry that one with a radius of 10 (for instance).
The implementation is reasonably straightforward – the main function, FindApollonianCircle(), is probably the hardest to understand – but I tried not to change anything very much with the algorithm (which means, just like last time, there may well be more elegant ways to solve the problem).
Throughout this series of investigations, I’ve been itching to switch across to F# to see what’s possible, there. In the next post or two in the series, we’ll take a look at implementing an Apollonian gasket in F# before moving on to the ultimately more relevant world of 3D.
