# How can one compute a 3D Periodic Delaunay Triangulation using DelaunayMesh?

Given a non-periodic point set, one can easily tetrahedralize it using the new (in V10) DelaunayMesh function. e.g.:

SeedRandom[0]
pts3d = RandomReal[4, {200, 3}];
del = DelaunayMesh[pts3d];
style = MapThread[Style[#1, Directive[#2]] &, {{1, 0, 2}, {{Thin, Purple},
{PointSize[0.02], Red}, {Opacity[0.5], Green}}}];
HighlightMesh[del, style]


Now consider the following periodic point set in 3D:

cavs = {{0., 1.199, 2.53}, {0., 1.265, 2.53}, {0.067, 1.199, 2.53}, {3.263,
1.199, 2.53}, {0.067, 1.265, 2.53}, {3.263, 1.265, 2.53}, {0.133,
1.199, 2.53}, {3.196, 1.199, 2.53}, {0.133, 1.265, 2.53}, {3.196,
1.265, 2.53}, {0.2, 1.199, 2.53}, {3.13, 1.199, 2.53}, {3.196,
1.199, 2.464}, {0.133, 1.265, 2.597}, {0.133, 1.332, 2.53}, {0.2,
1.265, 2.53}, {3.13, 1.265, 2.53}, {0.2, 1.199, 2.597}, {0.266,
1.199, 2.53}, {3.13, 1.199, 2.464}, {0.2, 1.265, 2.597}, {0.2,
1.332, 2.53}, {0.266, 1.265, 2.53}, {3.063, 1.265, 2.53}, {3.13,
1.265, 2.464}, {0.266, 1.132, 2.53}, {3.063, 1.199, 2.464}, {0.2,
1.265, 2.664}, {0.2, 1.332, 2.597}, {0.2, 1.398, 2.53}, {0.266,
1.332, 2.53}, {0.333, 1.265, 2.53}, {3.063, 1.265, 2.464}, {3.063,
1.332, 2.53}, {3.063, 1.199, 2.397}, {0.2, 1.398, 2.597}, {0.266,
1.332, 2.597}, {0.266, 1.398, 2.53}, {0.266, 1.332, 2.464}, {0.333,
1.332, 2.53}, {0.333, 1.265, 2.464}, {0.4, 1.265, 2.53}, {2.997,
1.265, 2.464}, {3.063, 1.265, 2.397}, {3.063, 1.332, 2.464}, {2.997,
1.332, 2.53}, {2.997, 1.199, 2.397}, {3.063, 1.132, 2.397}, {3.063,
1.199, 2.331}, {0.2, 1.465, 2.597}, {0.266, 1.398, 2.597}, {0.266,
1.332, 2.664}, {0.333, 1.332, 2.597}, {0.266, 1.465, 2.53}, {0.333,
1.398, 2.53}, {0.333, 1.332, 2.464}, {0.4, 1.332, 2.53}, {0.333,
1.199, 2.464}, {0.4, 1.265, 2.464}, {0.466, 1.265, 2.53}, {2.997,
1.265, 2.397}, {2.997, 1.332, 2.464}, {3.063, 1.265, 2.331}, {2.997,
1.332, 2.597}, {2.997, 1.398, 2.53}, {2.997, 1.199, 2.331}, {0.266,
1.465, 2.597}, {0.266, 1.398, 2.664}, {0.333, 1.398, 2.597}, {0.4,
1.398, 2.53}, {0.4, 1.332, 2.464}, {0.466, 1.332, 2.53}, {0.333,
1.132, 2.464}, {0.4, 1.199, 2.464}, {0.466, 1.265, 2.464}, {2.93,
1.265, 2.397}, {2.997, 1.332, 2.397}, {2.93, 1.332, 2.464}, {2.997,
1.398, 2.464}, {2.997, 1.398, 2.597}, {2.997, 1.132, 2.331}, {0.266,
1.465, 2.664}, {0.333, 1.465, 2.597}, {0.333, 1.398,
2.664}, {0.466, 1.332, 2.464}, {0.533, 1.332, 2.53}, {0.333, 1.065,
2.464}, {0.466, 1.265, 2.397}, {2.93, 1.332, 2.397}, {2.997, 1.332,
2.331}, {2.863, 1.332, 2.464}, {2.93, 1.132, 2.331}, {0.266, 1.465,
2.73}, {0.533, 1.332, 2.464}, {2.863, 1.332, 2.397}, {2.997, 1.398,
2.331}, {2.863, 1.132, 2.331}, {0.266, 1.465, 2.797}, {0.333, 1.465,
2.73}, {0.533, 1.332, 2.397}, {0.599, 1.332, 2.464}, {2.93, 1.398,
2.331}, {2.863, 1.199, 2.331}, {0.266, 1.465, 2.863}, {0.266, 1.532,
2.797}, {0.599, 1.332, 2.397}, {0.666, 1.332, 2.464}, {2.797,
1.199, 2.331}, {0.266, 1.532, 2.863}, {0.333, 1.532, 2.797}, {0.666,
1.332, 2.397}, {0.666, 1.332, 2.53}, {0.732, 1.332, 2.464}, {2.73,
1.199, 2.331}, {2.797, 1.265, 2.331}, {0.266, 1.598, 2.863}, {0.333,
1.532, 2.863}, {0.666, 1.398, 2.397}, {0.732, 1.332,
2.397}, {0.732, 1.332, 2.53}, {0.732, 1.398, 2.464}, {0.799, 1.332,
2.464}, {2.664, 1.199, 2.331}, {2.73, 1.132, 2.331}, {2.73, 1.265,
2.331}, {0.333, 1.598, 2.863}, {0.732, 1.398, 2.397}, {0.799, 1.332,
2.397}, {0.732, 1.332, 2.597}, {0.799, 1.398, 2.464}, {2.597,
1.199, 2.331}, {2.664, 1.199, 2.264}, {2.664, 1.199, 2.397}, {2.664,
1.265, 2.331}, {2.73, 1.132, 2.264}, {2.73, 1.265, 2.397}, {0.799,
1.398, 2.397}, {0.799, 1.398, 2.53}, {0.799, 1.465, 2.464}, {0.866,
1.398, 2.464}, {2.597, 1.199, 2.397}, {2.664, 1.199, 2.197}, {2.664,
1.265, 2.264}, {2.664, 1.265, 2.397}, {0.799, 1.398,
2.331}, {0.799, 1.465, 2.397}, {0.866, 1.398, 2.397}, {2.597, 1.132,
2.397}, {2.597, 1.265, 2.397}, {2.664, 1.332, 2.264}, {0.799,
1.465, 2.331}, {0.866, 1.398, 2.331}}


Here's what it looks like:

Graphics3D[{Red, PointSize[0.02], Point[cavs]}, BoxRatios -> {1, 1, 1}]


We tetrahedralize it naively:

cavdel = DelaunayMesh[cavs];


Visualize using the same style from before:

Show[HighlightMesh[cavdel, style], BoxRatios -> {1, 1, 1}]


Obviously, this is wrong as it's assuming the end of the box to be the end region and will lead to a larger volume than the true volume.

One can re-align the points in a box that gives them the minimum distance. When this is done the true arrangement of the points looks like this:

One can then easily use DelaunayMesh to tetrahedralize and obtain the following:

Which is the true Delaunay triangulation. The Delaunay tetrahedralization of the original periodic point set should look something like this:

Clearly, this is not quite right, since there are regions below and above the points that should fill up to the box edges and continue on the other side (wrap-around effect).

Since we can't give DelaunayMesh a Distance function AFAIK, my question is, given a set of periodic points in 3D how can one tetrahedralize it using DelaunayMesh?

Note: The length of the box is 3.2629 in each of x, y and z direction. Origin is (0, 0, 0) and Minimum image periodic boundary conditions were applied in all directions.

• I have some trouble understanding what you are looking for, perhaps you could rephrase a bit? – user21 Aug 1 '14 at 7:18
• @user21, Basically the distance between those points are not the normal EuclideanDistance so there is a wrap-around on points at the edge of the box. – RunnyKine Aug 1 '14 at 7:21
• Which dimensions are periodic ? all of them ? – lalmei Aug 1 '14 at 8:45
• @lalmei. Yes, all of them. – RunnyKine Aug 1 '14 at 8:47
• @Silvia. This is what I've done in the question, not just 6 directions but all 26 periodic boxes surrounding the center box. The problem with this approach is it's not feasible when you have 100's of such points in millions of configurations. – RunnyKine Aug 1 '14 at 8:54

The main idea is to find the gap position dimension by dimension. Take the first dimension for example, to accomplish that, we first project the hold points set to x axis, and do some very basic but fast statistics to locate the "gap", then we shift the coordinates according to the it.

Here is an example:

(* Size of the periodic cell: *)
ℒ = 10;

(* example data: *)
cavs = RandomVariate[
MultinormalDistribution[{3, 4, 5}, {{1, -(1/4), 1/3}, {-(1/4), 2/3, 1/5}, {1/3, 1/5, 1/2}}],
10^4] //
MapThread[#1@#2 &, {{Mod[#, ℒ, 4] &, Mod[#, ℒ, 2] &, Identity}, #}] &;

cavs // Graphics3D[{Blue, AbsolutePointSize[1], Point[#]},
PlotRange -> 2 {{0, ℒ}, {0, ℒ}, {0, ℒ}}, Axes -> True,
AxesLabel -> (Style[#, 15, Bold, Italic] & /@ {"x", "y", "z"}),
BoxRatios -> {1, 1, 1}] &


AbsoluteTiming[
(* coords is the coordinates of the points in any one dimension: *)
shiftParas = Module[{coords = #, crdSorted, threshold = 1, gapPos},
crdSorted = Sort[coords];
(* For sorted coords, gap position can be easily detected as a huge jump on a continuous distribution: *)
gapPos =
Differences[crdSorted] // Sign[threshold - #] & // FirstPosition[#, -1] &;
If[#, {#, crdSorted[[Join[gapPos, gapPos + 1]]] // Mean}, {#, crdSorted[[{1, -1}]] // Mean}] &@ListQ[gapPos]
] & /@ (cavs)
]

{0.007005, {{True, 8.14929}, {True, 8.73281}, {False, 5.22408}}}

(* The rest work is just shift the points according to the shiftParas: *)
cavsShifted =
If[#1[[1]],
Mod[#2, ℒ, #1[[2]]] - #1[[2]],
#2 + ℒ/2 - #1[[2]]
] &, {shiftParas, cavs}, 1];

cavsShifted //
Graphics3D[{Blue, AbsolutePointSize[1], Point[#]},
PlotRange -> {{0, ℒ}, {0, ℒ}, {0, ℒ}}, Axes -> True,
AxesLabel -> (Style[#, 15, Bold, Italic] & /@ {"x", "y", "z"}),
BoxRatios -> {1, 1, 1}] &


For $10^6$ points it will take about 1 sec on my laptop.

Note for a more sophisticated statistics, functions like HistogramList should be used rather than simply Differences. The timing will be a bit more, but I think is bearable.

## Update:

For OP's example, we have:

• @RunnyKine I believe it works, but please see my update and inform me my possible misunderstanding. (Have you set ℒ to your value i.e. 3.2629?) – Silvia Aug 2 '14 at 0:32
• @RunnyKine You're welcome. Thanks for acceptance. – Silvia Aug 2 '14 at 10:01