Repository: esa/polyhedral-gravity-model (language: python)

Issue:
Create pypi publish workflow
# Changelog

_53 Commit-Tortures later, finally, the CI runs for pypy, CPython on Windows, OSX, and Manylinux (due to psychological reasons, all 53 commits were rebased to one since the commit messages degraded to single letters)_

- Add workflow which builds the `sdist` and `wheels` for pypy, CPython on Windows, OSX, and Manylinux for all Python versions 3.6 - 3.11
- Add automatic upload to PyPi
- *Trusted Publishing* is enabled on the PyPi-side, so no secrets or similar need to be added to the repository



## Details

- Windows
  - We use OpenMP parallelization due to compile issues with GNU
  - We use GNU since MSVC is somehow unavailable in the Docker Image (?)
  - Because Windows GNU `cmath` header has a weird include bug, we do some magic with `CXXFLAGS=-D_hypot=hypot` to fix the namespace issue
- OSX
  - We use TBB
  - We don't use Apple Clang since the Docker Image uses an old MacOS, which has compile issues with C++17
  - We use GCC12
- Linux
  - We use TBB


Discussion / hints:
@gomezzz 

Yeah 😌, Windows made me a bit mad. OSX was, for once, pretty straightforward once I got the `manylinux` thing running (but that took a while, as well).

Definitely, a good idea to upload to testpypi anyway and import checking - a massive benefit for little to no work.
I updated the PR accordingly.

So you mean making a release from the branch of the Pull Request (this does not work, I think, since the PR comes from a fork)? Or similarly?
Oh no, already pays of. The Windows import fails
@gomezzz 

With the exception of the upload to real `PyPi`, everything is tested, working (compare https://github.com/schuhmaj/polyhedral-gravity-model-cibuild/actions/runs/5941981799), and ready for review again :)

Fix patch:
diff --git a/.github/workflows/wheels.yml b/.github/workflows/wheels.yml
index 2db1eed..39f5b87 100644
--- a/.github/workflows/wheels.yml
+++ b/.github/workflows/wheels.yml
@@ -7,9 +7,11 @@ on:
     # See https://docs.github.com/en/actions/reference/events-that-trigger-workflows#release
     # for more information on the release event
     prerelease: false
+  # Allows you to run this workflow manually from the Actions tab
+  workflow_dispatch:
 
 jobs:
-  # 1.a Buidl the wheels on a matrix of Windows, MacOS, and Linux plattforms usinbg cibuildwheel
+  # 1.a Buidl the wheels on a matrix of Windows, MacOS, and Linux platforms using cibuildwheel
   build_wheels:
     name: Build wheels on ${{ matrix.os }}
     runs-on: ${{ matrix.os }}
@@ -44,7 +46,7 @@ jobs:
           package-dir: .
           output-dir: dist
         if: matrix.os == 'macos-latest'
-      # Set up the Visual Studio environment on Windows (required, so that CNaje finds the compiler)
+      # Set up the Visual Studio environment on Windows (required, so that CMake finds the compiler)
       - uses: ilammy/msvc-dev-cmd@v1
         if: matrix.os == 'windows-latest'
       - name: Build wheels (Windows)
@@ -93,7 +95,7 @@ jobs:
           repository-url: https://test.pypi.org/legacy/
 
   # 3. Check if the package can be installed from testpypi
-  #    Notice, that this is more of an install test since the
+  #    Notice, that this is more of an installation test since the
   #    import check has already been done in the build_wheels job
   check_testpypi:
     needs: [upload_testpypi]
@@ -126,7 +128,7 @@ jobs:
       - name: Check import
         run: python -c "import polyhedral_gravity"
 
-  # 4. Upload the wheels to the actualy Python Package Index
+  # 4. Upload the wheels to the actually Python Package Index
   #    using trusted publishing
   upload_pypi:
     needs: [build_wheels, make_sdist, check_testpypi]
diff --git a/environment.yml b/environment.yml
new file mode 100644
index 0000000..6fcf810
--- /dev/null
+++ b/environment.yml
@@ -0,0 +1,7 @@
+name: polyhedral-gravity-env
+channels:
+  - conda-forge
+dependencies:
+  - numpy==1.23.4
+  - pytest
+  - pytest-xdist
\ No newline at end of file
diff --git a/script/time.py b/script/time.py
index 23a045c..db56e47 100755
--- a/script/time.py
+++ b/script/time.py
@@ -4,6 +4,7 @@
 import matplotlib.pyplot as plt
 import mesh_utility
 import timeit
+import pickle
 
 vertices, faces = mesh_utility.read_pk_file("/Users/schuhmaj/Programming/polyhedral-gravity-model/script/mesh/Eros.pk")
 vertices, faces = np.array(vertices), np.array(faces)
@@ -67,3 +68,24 @@
 print(f"--> Time taken: {delta:.3f} microseconds per point")
 
 
+########################################################################################################################
+# PICKLE SUPPORT
+########################################################################################################################
+
+
+with open("evaluable.pk", "wb") as f:
+    pickle.dump(evaluable, f, pickle.HIGHEST_PROTOCOL)
+    
+with open("evaluable.pk", "rb") as f:
+    evaluable2 = pickle.load(f)
+    
+
+start_time = timeit.default_timer()
+evaluable2(computation_points)
+end_time = timeit.default_timer()
+
+delta = (end_time - start_time) / N * 1e6
+# Print the time in milliseconds
+print("######## Pickle ########")
+print(f"--> 1 time {N} points with GravityEvaluable")
+print(f"--> Time taken: {delta:.3f} microseconds per point")
diff --git a/setup.py b/setup.py
index cd9453a..d7e73aa 100644
--- a/setup.py
+++ b/setup.py
@@ -147,7 +147,7 @@ def build_extension(self, ext):
 # --------------------------------------------------------------------------------
 setup(
     name="polyhedral_gravity",
-    version="2.0.3",
+    version="2.1",
     author="Jonas Schuhmacher",
     author_email="jonas.schuhmacher@tum.de",
     description="Package to compute full gravity tensor of a given constant density polyhedron for arbitrary points",
diff --git a/src/polyhedralGravity/calculation/GravityEvaluable.cpp b/src/polyhedralGravity/calculation/GravityEvaluable.cpp
index be399ed..36343af 100644
--- a/src/polyhedralGravity/calculation/GravityEvaluable.cpp
+++ b/src/polyhedralGravity/calculation/GravityEvaluable.cpp
@@ -52,10 +52,10 @@ namespace polyhedralGravity {
 
         if constexpr (Parallelization) {
             result = thrust::transform_reduce(thrust::device, zip1, zip2, &GravityEvaluable::evaluateFace, result,
-                                              util::operator+<double, Array3, Array6>);
+                                              util::operator+ < double, Array3, Array6 > );
         } else {
             result = thrust::transform_reduce(thrust::host, zip1, zip2, &GravityEvaluable::evaluateFace, result,
-                                              util::operator+<double, Array3, Array6>);
+                                              util::operator+ < double, Array3, Array6 > );
         }
 
         SPDLOG_LOGGER_DEBUG(PolyhedralGravityLogger::DEFAULT_LOGGER.getLogger(),
@@ -203,7 +203,11 @@ namespace polyhedralGravity {
             // indicators for numerical magnitudes appearing during the calculation:
             // planeDistance gets very big when far away, sum2 remains independently very small
             SPDLOG_LOGGER_WARN(PolyhedralGravityLogger::DEFAULT_LOGGER.getLogger(),
-                               "The results of point [{}, {}, {}] may be wrong due to floating point arithmetic");
+                               "While evaluating the plane with coordinates v1 = [{}, {}, {}], v2 = [{}, {}, {}], "
+                               "v3 = [{}, {}, {}] (with computation point re-located at the origin) a "
+                               "significant difference of magnitudes occurred during the evaluation. "
+                               "This may lead to numerically unstable results!", face[0][0], face[0][1], face[0][2],
+                               face[1][0], face[1][1], face[1][2], face[2][0], face[2][1], face[2][2]);
         }
 
         //6. Step: Sum for tensor
@@ -229,10 +233,14 @@ namespace polyhedralGravity {
 
     std::string GravityEvaluable::toString() const {
         std::stringstream ss;
-        ss  << "<polyhedral_gravity.GravityEvaluable, density=" << _density
-            << ", vertices= " << _polyhedron.countVertices()
-            << ", faces= " << _polyhedron.countFaces() << ">";
+        ss << "<polyhedral_gravity.GravityEvaluable, density=" << _density << ", vertices= "
+           << _polyhedron.countVertices() << ", faces= " << _polyhedron.countFaces() << ">";
         return ss.str();
     }
 
+    std::tuple<Polyhedron, double, std::vector<Array3Triplet>, std::vector<Array3>, std::vector<Array3Triplet>>
+    GravityEvaluable::getState() const {
+        return std::make_tuple(_polyhedron, _density, _segmentVectors, _planeUnitNormals, _segmentUnitNormals);
+    }
+
 }
diff --git a/src/polyhedralGravity/calculation/GravityEvaluable.h b/src/polyhedralGravity/calculation/GravityEvaluable.h
index 1b053ca..2c847f8 100644
--- a/src/polyhedralGravity/calculation/GravityEvaluable.h
+++ b/src/polyhedralGravity/calculation/GravityEvaluable.h
@@ -59,16 +59,32 @@ namespace polyhedralGravity {
          * @param polyhedralSource - the vertices & faces of the polyhedron as tuple or the filenames of the files
          * @param density - the constant density in [kg/m^3]
          */
-        GravityEvaluable(
-                const PolyhedralSource &polyhedralSource,
-                double density) : _polyhedron{
+        GravityEvaluable(const PolyhedralSource &polyhedralSource, double density) : _polyhedron{
                 std::holds_alternative<std::tuple<std::vector<Array3>, std::vector<IndexArray3>>>(polyhedralSource)
                 ? Polyhedron{std::get<std::tuple<std::vector<Array3>, std::vector<IndexArray3>>>(polyhedralSource)}
-                : TetgenAdapter(std::get<std::vector<std::string>>(polyhedralSource)).getPolyhedron()
-        }, _density{density} {
+                : TetgenAdapter(std::get<std::vector<std::string>>(polyhedralSource)).getPolyhedron()},
+                                                                                     _density{density} {
             this->prepare();
         }
 
+        /**
+         * Instantiates a GravityEvaluable with a given constant density polyhedron and caches.
+         * This is for restoring a GravityEvaluable from a previous state.
+         * @param polyhedron - the polyhedron
+         * @param density - the constant density in [kg/m^3]
+         * @param segmentVectors - the segment vectors
+         * @param planeUnitNormals - the plane unit normals
+         * @param segmentUnitNormals - the segment unit normals
+         */
+        GravityEvaluable(const Polyhedron &polyhedron, double density, const std::vector<Array3Triplet> &segmentVectors,
+                         const std::vector<Array3> &planeUnitNormals,
+                         const std::vector<Array3Triplet> &segmentUnitNormals) : _polyhedron{polyhedron},
+                                                                                 _density{density},
+                                                                                 _segmentVectors{segmentVectors},
+                                                                                 _planeUnitNormals{planeUnitNormals},
+                                                                                 _segmentUnitNormals{
+                                                                                         segmentUnitNormals} {}
+
         /**
          * Evaluates the polyhedrale gravity model for a given constant density polyhedron at computation
          * point P. Wrapper for evaluate<parallelization>.
@@ -100,6 +116,14 @@ namespace polyhedralGravity {
          */
         std::string toString() const;
 
+        /**
+         * Returns the polyhedron, the density and the internal caches.
+         *
+         * @return tuple of polyhedron, density, segmentVectors, planeUnitNormals and segmentUnitNormals
+         */
+        std::tuple<Polyhedron, double, std::vector<Array3Triplet>, std::vector<Array3>, std::vector<Array3Triplet>>
+        getState() const;
+
     private:
 
         /**
@@ -120,6 +144,7 @@ namespace polyhedralGravity {
         template<bool Parallelization = true>
         GravityModelResult evaluate(const Array3 &computationPoint) const;
 
+
         /**
          * Evaluates the polyhedrale gravity model for a given constant density polyhedron at computation
          * at multiple computation points.
@@ -130,7 +155,6 @@ namespace polyhedralGravity {
         template<bool Parallelization = true>
         std::vector<GravityModelResult> evaluate(const std::vector<Array3> &computationPoints) const;
 
-
         /**
          * Evaluates the polyhedrale gravity model for a given constant density polyhedron at computation a certain face.
          * @param tuple - consisting of face, segmentVectors, planeUnitNormal and segmentUnitNormals
diff --git a/src/polyhedralGravityPython/PolyhedralGravityPython.cpp b/src/polyhedralGravityPython/PolyhedralGravityPython.cpp
index 4cd4e5a..c7c6600 100644
--- a/src/polyhedralGravityPython/PolyhedralGravityPython.cpp
+++ b/src/polyhedralGravityPython/PolyhedralGravityPython.cpp
@@ -83,7 +83,29 @@ PYBIND11_MODULE(polyhedral_gravity, m) {
                     Returns:
                         Either a tuple of potential, acceleration and second derivatives at the computation points or
                         if multiple computation points are given, the result is a list of tuples
-                )mydelimiter", py::arg("computation_points"), py::arg("parallel") = true);
+                )mydelimiter", py::arg("computation_points"), py::arg("parallel") = true)
+            .def(py::pickle(
+                    [](const GravityEvaluable &evaluable) {
+                        const auto&[polyhedron, density, segmentVectors, planeUnitNormals, segmentUnitNormals] =
+                                evaluable.getState();
+                        return py::make_tuple(polyhedron.getVertices(), polyhedron.getFaces(), density, segmentVectors,
+                                              planeUnitNormals, segmentUnitNormals);
+                    },
+                    [](const py::tuple &tuple) {
+                        constexpr size_t tupleSize = 6;
+                        if (tuple.size() != tupleSize) {
+                            throw std::runtime_error("Invalid state!");
+                        }
+                        Polyhedron polyhedron {
+                            tuple[0].cast<std::vector<Array3>>(), tuple[1].cast<std::vector<IndexArray3>>()
+                        };
+                        GravityEvaluable evaluable{
+                            polyhedron, tuple[2].cast<double>(), tuple[3].cast<std::vector<Array3Triplet>>(),
+                            tuple[4].cast<std::vector<Array3>>(), tuple[5].cast<std::vector<Array3Triplet>>()
+                        };
+                        return evaluable;
+                    }
+            ));
 
     py::module_ utility = m.def_submodule("utility",
                                           "This submodule contains useful utility functions like parsing meshes "


Test patch:
diff --git a/.github/workflows/build-and-test.yml b/.github/workflows/build-and-test.yml
new file mode 100644
index 0000000..548883e
--- /dev/null
+++ b/.github/workflows/build-and-test.yml
@@ -0,0 +1,65 @@
+name: Build & Test
+on:
+  # Triggers the workflow on push or pull request events but only for the master branch
+  push:
+    branches: [ main ]
+  pull_request:
+    branches: [ main ]
+  # Allows you to run this workflow manually from the Actions tab
+  workflow_dispatch:
+
+# A workflow run is made up of one or more jobs that can run sequentially or in parallel
+jobs:
+  # Build & Test on Linux
+  build-and-test-linux:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v3
+      - name: Build
+        run: |
+          mkdir build
+          cd build
+          cmake -DCMAKE_BUILD_TYPE=Release ..
+          cmake --build .
+      - name: Test
+        run: cd build/test && ./polyhedralGravity_test
+
+  # Install the python interface by building from source and run pytest
+  pytest-linux:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v3
+      - name: Install Ninja Build
+        run: |
+          sudo apt-get install ninja-build -y
+      - name: Install Conda environment from environment.yml
+        uses: mamba-org/provision-with-micromamba@main
+        with:
+          environment-file: environment.yml
+          cache-downloads: true
+          cache-env: true
+      - name: Install & Test polyhedral-gravity
+        shell: bash -l {0}
+        run: |
+          pip install . -vv --no-build-isolation
+          pytest -n 3
+          
+
+  # Builds the polyhedral gravity python interface the interface
+  # No testing on Windows since it takes too long
+  # Enables the long paths feature on Windows (newer thrust versions require this)
+  build-windows:
+    runs-on: windows-latest
+    steps:
+      - uses: actions/checkout@v3
+      - uses: ilammy/msvc-dev-cmd@v1
+      - name: Build
+        run: |
+          New-ItemProperty -Path "HKLM:\SYSTEM\CurrentControlSet\Control\FileSystem" -Name "LongPathsEnabled" -Value 1 -PropertyType DWORD -Force
+          git config --system core.longpaths true
+          mkdir build
+          cd build
+          cmake -DCMAKE_BUILD_TYPE=Release -DBUILD_POLYHEDRAL_PYTHON_INTERFACE=ON ..
+          cmake --build . --config Release
+
+      
diff --git a/.github/workflows/ctest.yml b/.github/workflows/ctest.yml
deleted file mode 100644
index 0913fd3..0000000
--- a/.github/workflows/ctest.yml
+++ /dev/null
@@ -1,60 +0,0 @@
-# This is a basic workflow to help you get started with Actions
-
-name: CI
-
-# Controls when the action will run. 
-on:
-  # Triggers the workflow on push or pull request events but only for the master branch
-  push:
-    branches: [ main ]
-  pull_request:
-    branches: [ main ]
-
-  # Allows you to run this workflow manually from the Actions tab
-  workflow_dispatch:
-
-# A workflow run is made up of one or more jobs that can run sequentially or in parallel
-jobs:
-  # This workflow contains a single job called "build-and-test"
-  build-and-test-linux:
-    # The type of runner that the job will run on
-    runs-on: ubuntu-latest
-
-    # Steps represent a sequence of tasks that will be executed as part of the job
-    steps:
-      # Checks-out your repository under $GITHUB_WORKSPACE, so your job can access it
-      - uses: actions/checkout@v2
-
-      - name: Build
-        run: |
-          mkdir build
-          cd build
-          cmake -DCMAKE_BUILD_TYPE=Release ..
-          cmake --build .
-          
-      - name: Test
-        run: cd build/test && ./polyhedralGravity_test
-
-  build-windows:
-    # The type of runner that the job will run on
-    runs-on: windows-latest
-
-    # Steps represent a sequence of tasks that will be executed as part of the job
-    steps:
-      # Checks-out your repository under $GITHUB_WORKSPACE, so your job can access it
-      - uses: actions/checkout@v2
-
-      # Enables Long Paths on Windows (newer thrust versions require this)
-      - name: Build
-        run: |
-          New-ItemProperty -Path "HKLM:\SYSTEM\CurrentControlSet\Control\FileSystem" -Name "LongPathsEnabled" -Value 1 -PropertyType DWORD -Force
-          git config --system core.longpaths true
-          mkdir build
-          cd build
-          cmake -DCMAKE_BUILD_TYPE=Release ..
-          cmake --build . --config Release
-
-# No testing on Windows GitHub since it takes too long
-#      - name: Test
-#        run: cd build && ctest
-      
diff --git a/pytest.ini b/pytest.ini
new file mode 100644
index 0000000..e6a7c51
--- /dev/null
+++ b/pytest.ini
@@ -0,0 +1,4 @@
+[pytest]
+minversion = 6.0
+testpaths =
+    test
\ No newline at end of file
diff --git a/test/python/test_polyhedral_gravity.py b/test/python/test_polyhedral_gravity.py
new file mode 100644
index 0000000..6f9fa8d
--- /dev/null
+++ b/test/python/test_polyhedral_gravity.py
@@ -0,0 +1,136 @@
+from typing import Tuple, List, Union
+import polyhedral_gravity
+import numpy as np
+import pickle
+import pytest
+from pathlib import Path
+
+CUBE_VERTICES = np.array([
+    [-1, -1, -1],
+    [1, -1, -1],
+    [1, 1, -1],
+    [-1, 1, -1],
+    [-1, -1, 1],
+    [1, -1, 1],
+    [1, 1, 1],
+    [-1, 1, 1]
+])
+
+CUBE_FACES = np.array([
+    [1, 3, 2],
+    [0, 3, 1],
+    [0, 1, 5],
+    [0, 5, 4],
+    [0, 7, 3],
+    [0, 4, 7],
+    [1, 2, 6],
+    [1, 6, 5],
+    [2, 3, 6],
+    [3, 7, 6],
+    [4, 5, 6],
+    [4, 6, 7]
+])
+
+DENSITY = 1.0
+
+CUBE_VERTICES_FILE = Path("test/resources/cube.node")
+CUBE_FACE_FILE = Path("test/resources/cube.face")
+
+
+@pytest.fixture(scope="session")
+def reference_solution() -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
+    """Reads the cube reference solution from a file and returns it as a tuple"""
+    reference_file_path = Path("test/resources/analytic_cube_solution.txt")
+    reference = np.loadtxt(reference_file_path)
+    points = reference[:, :3]
+    expected_potential = reference[:, 3].flatten()
+    expected_acceleration = reference[:, 4:]
+    return points, expected_potential, expected_acceleration
+
+
+@pytest.mark.parametrize(
+    "polyhedral_source", [(CUBE_VERTICES, CUBE_FACES), [str(CUBE_VERTICES_FILE), str(CUBE_FACE_FILE)]],
+    ids=["with_arrays", "with_file_input"]
+)
+def test_polyhedral_gravity(
+        polyhedral_source: Union[Tuple[np.ndarray, np.ndarray], List[str]],
+        reference_solution: Tuple[np.ndarray, np.ndarray, np.ndarray]) -> None:
+    """Checks that the evaluate function the correct results and
+    is callable with file/ array inputs.
+    """
+    points, expected_potential, expected_acceleration = reference_solution
+    sol = polyhedral_gravity.evaluate(
+        polyhedral_source=polyhedral_source,
+        density=DENSITY,
+        computation_points=points,
+        parallel=True
+    )
+    sol = np.array([np.array(x) for x in sol])
+    actual_potential = sol[:, 0].flatten()
+    assert np.testing.assert_array_almost_equal(actual_potential, expected_potential) is None
+
+    actual_acceleration = np.array([np.array(x) for x in sol[:, 1]])
+    assert np.testing.assert_array_almost_equal(actual_acceleration, expected_acceleration) is None
+
+
+@pytest.mark.parametrize(
+    "polyhedral_source", [(CUBE_VERTICES, CUBE_FACES), [str(CUBE_VERTICES_FILE), str(CUBE_FACE_FILE)]],
+    ids=["with_arrays", "with_file_input"]
+)
+def test_polyhedral_gravity_evaluable(
+        polyhedral_source: Union[Tuple[np.ndarray, np.ndarray], List[str]],
+        reference_solution: Tuple[np.ndarray, np.ndarray, np.ndarray]) -> None:
+    """Checks that the evaluable produces the correct results and
+    is instantiable with file/ array inputs.
+    """
+    points, expected_potential, expected_acceleration = reference_solution
+    evaluable = polyhedral_gravity.GravityEvaluable(
+        polyhedral_source=polyhedral_source,
+        density=DENSITY,
+    )
+    sol = evaluable(
+        computation_points=points,
+        parallel=True
+    )
+    sol = np.array([np.array(x) for x in sol])
+    actual_potential = sol[:, 0].flatten()
+    assert np.testing.assert_array_almost_equal(actual_potential, expected_potential) is None
+
+    actual_acceleration = np.array([np.array(x) for x in sol[:, 1]])
+    assert np.testing.assert_array_almost_equal(actual_acceleration, expected_acceleration) is None
+
+
+@pytest.mark.parametrize(
+    "polyhedral_source", [(CUBE_VERTICES, CUBE_FACES), [str(CUBE_VERTICES_FILE), str(CUBE_FACE_FILE)]],
+    ids=["with_arrays", "with_file_input"]
+)
+def test_polyhedral_evaluable_pickle(
+        polyhedral_source: Union[Tuple[np.ndarray, np.ndarray], List[str]],
+        reference_solution: Tuple[np.ndarray, np.ndarray, np.ndarray],
+        tmp_path: Path,) -> None:
+    """Tests that the evaluable can be pickled and unpicked and that the results
+    are still correct afterward (i.e. that the internal cache is correctly
+    pickled and unpicked as well).
+    """
+    points, expected_potential, expected_acceleration = reference_solution
+    initial_evaluable = polyhedral_gravity.GravityEvaluable(
+        polyhedral_source=polyhedral_source,
+        density=DENSITY
+    )
+    pickle_output = tmp_path.joinpath("evaluable.pk")
+    with open(pickle_output, "wb") as f:
+        pickle.dump(initial_evaluable, f, pickle.HIGHEST_PROTOCOL)
+
+    with open(pickle_output, "rb") as f:
+        read_evaluable = pickle.load(f)
+
+    sol = read_evaluable(
+        computation_points=points,
+        parallel=True
+    )
+    sol = np.array([np.array(x) for x in sol])
+    actual_potential = sol[:, 0].flatten()
+    assert np.testing.assert_array_almost_equal(actual_potential, expected_potential) is None
+
+    actual_acceleration = np.array([np.array(x) for x in sol[:, 1]])
+    assert np.testing.assert_array_almost_equal(actual_acceleration, expected_acceleration) is None
diff --git a/test/resources/cube.face b/test/resources/cube.face
new file mode 100644
index 0000000..3174a32
--- /dev/null
+++ b/test/resources/cube.face
@@ -0,0 +1,15 @@
+# Number of faces, boundary marker off
+12 0
+# Face index, nodes of face
+0    1  3  2
+1    0  3  1
+2    0  1  5
+3    0  5  4
+4    0  7  3
+5    0  4  7
+6    1  2  6
+7    1  6  5
+8    2  3  6
+9    3  7  6
+10   4  5  6
+11   4  6  7
diff --git a/test/resources/cube.node b/test/resources/cube.node
new file mode 100644
index 0000000..9a77cfb
--- /dev/null
+++ b/test/resources/cube.node
@@ -0,0 +1,11 @@
+# Node count, 3 dimensions, no attribute, no boundary marker
+8 3 0 0
+# Node index, node coordinates
+0 -1 -1 -1
+1  1 -1 -1
+2  1  1 -1
+3 -1  1 -1
+4 -1 -1  1
+5  1 -1  1
+6  1  1  1
+7 -1  1  1