6
$\begingroup$

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]

Mathematica graphics

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}]

Mathematica graphics

We tetrahedralize it naively:

cavdel = DelaunayMesh[cavs];

Visualize using the same style from before:

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

Mathematica graphics

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:

Mathematica graphics

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

Mathematica graphics

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

Mathematica graphics

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.

$\endgroup$
  • $\begingroup$ I have some trouble understanding what you are looking for, perhaps you could rephrase a bit? $\endgroup$ – user21 Aug 1 '14 at 7:18
  • $\begingroup$ @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. $\endgroup$ – RunnyKine Aug 1 '14 at 7:21
  • 1
    $\begingroup$ Which dimensions are periodic ? all of them ? $\endgroup$ – lalmei Aug 1 '14 at 8:45
  • $\begingroup$ @lalmei. Yes, all of them. $\endgroup$ – RunnyKine Aug 1 '14 at 8:47
  • $\begingroup$ @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. $\endgroup$ – RunnyKine Aug 1 '14 at 8:54
4
$\begingroup$

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}] &

original data

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 = 
        MapThread[
                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}] &

shifted data

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:

shifted data for OP

$\endgroup$
  • $\begingroup$ @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?) $\endgroup$ – Silvia Aug 2 '14 at 0:32
  • $\begingroup$ @RunnyKine You're welcome. Thanks for acceptance. $\endgroup$ – Silvia Aug 2 '14 at 10:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.