Started implementing 3D moving mesh algorithm. Not really coupled to SWIFT yet.

.gitignore

@@ -79,7 +79,7 @@ tests/testRiemannTRRS
 tests/testRiemannHLLC
 tests/testMatrixInversion
 tests/testVoronoi1D
-tests/threadpool_test
+tests/testVoronoi3D
 
 theory/latex/swift.pdf
 theory/kernel/kernels.pdf

src/hydro/Shadowswift/voronoi3d_cell.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Bert Vandenbroucke (bert.vandenbroucke@gmail.com).
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+#ifndef SWIFT_VORONOI1D_CELL_H
+#define SWIFT_VORONOI1D_CELL_H
+
+/* Maximal number of neighbours that can be stored in a voronoi_cell struct */
+#define VORONOI3D_MAXNUMNGB 100
+/* Maximal number of vertices that can be stored in a voronoi_cell struct */
+#define VORONOI3D_MAXNUMVERT 500
+/* Maximal number of edges that can be stored in a voronoi_cell struct */
+#define VORONOI3D_MAXNUMEDGE 1500
+/* Maximal number of faces that can be stored in a voronoi_cell struct */
+#define VORONOI3D_MAXFACE 100
+
+/* 3D Voronoi cell */
+struct voronoi_cell {
+
+  /* The position of the generator of the cell. */
+  double x[3];
+
+  /* The volume of the 3D cell. */
+  float volume;
+
+  /* The centroid of the cell. */
+  float centroid[3];
+
+  /* Number of cell vertices. */
+  int nvert;
+
+  /* Vertex coordinates. */
+  float vertices[3 * VORONOI3D_MAXNUMVERT];
+
+  /* Number of edges for every vertex. */
+  int orders[VORONOI3D_MAXNUMVERT];
+
+  /* Offsets of the edges, edgeindices and neighbours corresponding to a
+     particular vertex in the internal arrays */
+  int offsets[VORONOI3D_MAXNUMVERT];
+
+  /* Edge information. */
+  int edges[VORONOI3D_MAXNUMEDGE];
+
+  /* Additional edge information. */
+  int edgeindices[VORONOI3D_MAXNUMEDGE];
+
+  /* Neighbour information. */
+  unsigned long long ngbs[VORONOI3D_MAXNUMEDGE];
+
+  /* Number of faces of the cell. */
+  int nface;
+
+  /* Surface areas of the cell faces. */
+  float face_areas[VORONOI3D_MAXFACE];
+
+  /* Midpoints of the cell faces. */
+  float face_midpoints[VORONOI3D_MAXFACE][3];
+};
+
+/**
+ * @brief Copy the contents of the 3D Voronoi cell pointed to by source into the
+ * 3D Voronoi cell pointed to by destination
+ *
+ * @param source Pointer to a 3D Voronoi cell to read from.
+ * @param destination Pointer to a 3D Voronoi cell to write to.
+ */
+__attribute__((always_inline)) INLINE void voronoi3d_cell_copy(
+    struct voronoi_cell *source, struct voronoi_cell *destination) {
+
+  /* Copy the position of the generator of the cell. */
+  destination->x[0] = source->x[0];
+  destination->x[1] = source->x[1];
+  destination->x[2] = source->x[2];
+
+  /* Copy the volume of the 3D cell. */
+  destination->volume = source->volume;
+
+  /* Copy the centroid of the cell. */
+  destination->centroid[0] = source->centroid[0];
+  destination->centroid[1] = source->centroid[1];
+  destination->centroid[2] = source->centroid[2];
+
+  /* Copy the number of cell vertices. */
+  destination->nvert = source->nvert;
+
+  /* Copy the vertex coordinates. We only copy the 3*nvert first coordinates. */
+  for (int i = 0; i < 3 * source->nvert; ++i) {
+    destination->vertices[i] = source->vertices[i];
+  }
+
+  /* Copy the number of edges for every vertex. Again, we only copy the nvert
+     first values. */
+  for (int i = 0; i < source->nvert; ++i) {
+    destination->orders[i] = source->orders[i];
+  }
+
+  /* Copy the nvert first values of the offsets. */
+  for (int i = 0; i < source->nvert; ++i) {
+    destination->offsets[i] = source->offsets[i];
+  }
+
+  /* Copy the edge information. No idea how many edges we have, so we copy
+     everything. */
+  for (int i = 0; i < VORONOI3D_MAXNUMEDGE; ++i) {
+    destination->edges[i] = source->edges[i];
+  }
+
+  /* Copy all additional edge information. */
+  for (int i = 0; i < VORONOI3D_MAXNUMEDGE; ++i) {
+    destination->edgeindices[i] = source->edgeindices[i];
+  }
+
+  /* Copy neighbour information. Since neighbours are stored per edge, the total
+     number of neighbours in this list is larger than numngb and we copy
+     everything. */
+  for (int i = 0; i < VORONOI3D_MAXNUMEDGE; ++i) {
+    destination->ngbs[i] = source->ngbs[i];
+  }
+
+  /* Copy the number of faces of the cell. */
+  destination->nface = source->nface;
+
+  /* Copy the surface areas of the cell faces. We only copy the nface first
+     values. */
+  for (int i = 0; i < source->nface; ++i) {
+    destination->face_areas[i] = source->face_areas[i];
+  }
+
+  /* Copy the midpoints of the cell faces. */
+  for (int i = 0; i < source->nface; ++i) {
+    destination->face_midpoints[i][0] = source->face_midpoints[i][0];
+    destination->face_midpoints[i][1] = source->face_midpoints[i][1];
+    destination->face_midpoints[i][2] = source->face_midpoints[i][2];
+  }
+}
+
+#endif  // SWIFT_VORONOI1D_CELL_H

src/hydro/Shadowswift/voronoi_algorithm.h

@@ -26,7 +26,6 @@
 #warning "2D moving mesh not implemented yet!"
 #include "voronoi2d_algorithm.h"
 #elif defined(HYDRO_DIMENSION_3D)
-#warning "3D moving mesh not implemented yet!"
 #include "voronoi3d_algorithm.h"
 #else
 #error "You have to select a dimension for the hydro!"

src/hydro/Shadowswift/voronoi_cell.h

@@ -26,7 +26,6 @@
 #warning "2D moving mesh not implemented yet!"
 #include "voronoi2d_cell.h"
 #elif defined(HYDRO_DIMENSION_3D)
-#warning "3D moving mesh not implemented yet!"
 #include "voronoi3d_algorithm.h"
 #else
 #error "You have to select a dimension for the hydro!"

tests/Makefile.am

@@ -25,7 +25,7 @@ TESTS = testGreetings testMaths testReading.sh testSingle testKernel testSymmetr
         testPair.sh testPairPerturbed.sh test27cells.sh test27cellsPerturbed.sh  \
         testParser.sh testSPHStep test125cells.sh testKernelGrav testFFT \
         testAdiabaticIndex testRiemannExact testRiemannTRRS testRiemannHLLC \
-        testMatrixInversion testVoronoi1D testThreadpool
+        testMatrixInversion testVoronoi1D testVoronoi3D testThreadpool
 
 # List of test programs to compile
 check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
@@ -33,7 +33,8 @@ check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
                  testKernel testKernelGrav testFFT testInteractions testMaths \
                  testSymmetry testThreadpool \
                  testAdiabaticIndex testRiemannExact testRiemannTRRS \
-                 testRiemannHLLC testMatrixInversion testVoronoi1D
+                 testRiemannHLLC testMatrixInversion testVoronoi1D \
+                 testVoronoi3D
 
 # Sources for the individual programs
 testGreetings_SOURCES = testGreetings.c
@@ -78,6 +79,8 @@ testMatrixInversion_SOURCES = testMatrixInversion.c
 
 testVoronoi1D_SOURCES = testVoronoi1D.c
 
+testVoronoi3D_SOURCES = testVoronoi3D.c
+
 testThreadpool_SOURCES = testThreadpool.c
 
 # Files necessary for distribution

tests/testVoronoi3D.c

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (C) 2016 Bert Vandenbroucke (bert.vandenbroucke@gmail.com).
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+#include <stdlib.h>
+#include "error.h"
+#include "hydro/Shadowswift/voronoi3d_algorithm.h"
+
+/**
+ * @brief Check if voronoi_volume_tetrahedron() works
+ */
+void test_voronoi_volume_tetrahedron() {
+  float v1[3] = {0., 0., 0.};
+  float v2[3] = {0., 0., 1.};
+  float v3[3] = {0., 1., 0.};
+  float v4[3] = {1., 0., 0.};
+
+  float V = voronoi_volume_tetrahedron(v1, v2, v3, v4);
+  assert(V == 1.0f / 6.0f);
+}
+
+/**
+ * @brief Check if voronoi_centroid_tetrahedron() works
+ */
+void test_voronoi_centroid_tetrahedron() {
+  float v1[3] = {0., 0., 0.};
+  float v2[3] = {0., 0., 1.};
+  float v3[3] = {0., 1., 0.};
+  float v4[3] = {1., 0., 0.};
+
+  float centroid[3];
+  voronoi_centroid_tetrahedron(centroid, v1, v2, v3, v4);
+  assert(centroid[0] == 0.25f);
+  assert(centroid[1] == 0.25f);
+  assert(centroid[2] == 0.25f);
+}
+
+/**
+ * @brief Check if voronoi_calculate_cell() works
+ */
+void test_calculate_cell() {
+  struct voronoi_cell cell;
+
+  cell.x[0] = 0.5f;
+  cell.x[1] = 0.5f;
+  cell.x[2] = 0.5f;
+
+  /* Initialize the cell to a large cube. */
+  voronoi_initialize(&cell);
+  /* Calculate the volume and centroid of the large cube. */
+  voronoi_calculate_cell(&cell);
+
+  /* Update these values if you ever change to another large cube! */
+  assert(cell.volume == 27.0f);
+  assert(cell.centroid[0] = 0.5f);
+  assert(cell.centroid[1] = 0.5f);
+  assert(cell.centroid[2] = 0.5f);
+}
+
+int main() {
+
+  /* Check basic Voronoi cell functions */
+  test_voronoi_volume_tetrahedron();
+  test_voronoi_centroid_tetrahedron();
+  test_calculate_cell();
+
+  /* Create a Voronoi cell */
+  double x[3] = {0.5f, 0.5f, 0.5f};
+  struct voronoi_cell cell;
+  voronoi_cell_init(&cell, x);
+
+  /* Interact with neighbours */
+  float x0[3] = {0.5f, 0.0f, 0.0f};
+  float x1[3] = {-0.5f, 0.0f, 0.0f};
+  float x2[3] = {0.0f, 0.5f, 0.0f};
+  float x3[3] = {0.0f, -0.5f, 0.0f};
+  float x4[3] = {0.0f, 0.0f, 0.5f};
+  float x5[3] = {0.0f, 0.0f, -0.5f};
+  voronoi_cell_interact(&cell, x0, 1);
+  voronoi_cell_interact(&cell, x1, 2);
+  voronoi_cell_interact(&cell, x2, 3);
+  voronoi_cell_interact(&cell, x3, 4);
+  voronoi_cell_interact(&cell, x4, 5);
+  voronoi_cell_interact(&cell, x5, 6);
+
+  /* Interact with some more neighbours to check if they are properly ignored */
+  float xE0[3] = {0.6f, 0.0f, 0.1f};
+  float xE1[3] = {-0.7f, 0.2f, 0.04f};
+  voronoi_cell_interact(&cell, xE0, 7);
+  voronoi_cell_interact(&cell, xE1, 8);
+
+  /* Finalize cell and check results */
+  voronoi_cell_finalize(&cell);
+
+  if (cell.volume != 0.125f) {
+    error("Wrong volume: %g!", cell.volume);
+  }
+  if (fabs(cell.centroid[0] - 0.5f) > 1.e-5f ||
+      fabs(cell.centroid[1] - 0.5f) > 1.e-5f ||
+      fabs(cell.centroid[2] - 0.5f) > 1.e-5f) {
+    error("Wrong centroid: %g %g %g!", cell.centroid[0], cell.centroid[1],
+          cell.centroid[2]);
+  }
+  /* Check neighbour order */
+  // TODO
+
+  /* Check face method */
+  float A[3];
+  float midpoint[3];
+  voronoi_get_face(&cell, 1, A, midpoint);
+  // TODO
+
+  return 0;
+}