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

Issue:
Allow to chose parallel or serial execution at runtime
The current implementation is compiled with either `CPP` (serial), `OMP` or `TBB` backend.

For E.g. the package on `conda` is always parallelized with `OMP`. This needs to change.

Thus, enable the user to choose `serial` or `parallel` execution at runtime via an argument given to the evaluated parameter. I.e. make the `CPP` backend always usable.


Discussion / hints:

Fix patch:
diff --git a/.gitignore b/.gitignore
index 75f77c2..3fe54c9 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,4 +1,5 @@
 .idea
 cmake-build-*
 build
-polyhedral_gravity.egg-info
\ No newline at end of file
+polyhedral_gravity.egg-info
+dist
\ No newline at end of file
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 2d5434e..4148467 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -74,62 +74,10 @@ endif ()
 thrust_create_target(Thrust HOST CPP DEVICE ${POLYHEDRAL_GRAVITY_PARALLELIZATION})
 message(STATUS "Set Parallelization: HOST CPP DEVICE ${POLYHEDRAL_GRAVITY_PARALLELIZATION}")
 
-################################
-# Building the Polyhedral Libray
-################################
-if (BUILD_POLYHEDRAL_GRAVITY_LIBRARY)
-    # YAML is only required for the library target, but not for e.g. the Python interface
-    include(yaml)
-
-    # Including the actual Polyhedral source files
-    file(GLOB_RECURSE SRC
-        "${PROJECT_SOURCE_DIR}/src/polyhedralGravity/*.h"
-        "${PROJECT_SOURCE_DIR}/src/polyhedralGravity/*.cpp")
-
-    add_library(${PROJECT_NAME}_lib OBJECT ${SRC})
-
-    # Adds the include Path PROJECT_SOURCE_DIR/src to the target polyhedralGravity
-    target_include_directories(${PROJECT_NAME}_lib PUBLIC
-            "${PROJECT_SOURCE_DIR}/src"
-            )
-
-    # Link libaries
-    target_link_libraries(${PROJECT_NAME}_lib
-            spdlog::spdlog
-            tetgen
-            yaml-cpp
-            xsimd
-            )
-
-    # Adds the include Path PROJECT_SOURCE_DIR/src to the target polyhedralGravity
-    target_include_directories(${PROJECT_NAME}_lib PUBLIC
-            "${PROJECT_SOURCE_DIR}/src"
-            )
-
-    # Link libraries needed in all targets
-    target_link_libraries(${PROJECT_NAME}_lib
-            spdlog::spdlog
-            yaml-cpp
-            tetgen
-            xsimd
-            Thrust
-            )
-
-endif()
-
-#####################################
-# Building the Polyhedral Executable
-#####################################
-if (BUILD_POLYHEDRAL_GRAVITY_EXECUTABLE)
-    # Main file, containing the gravity model's driver
-    set(mainFile ${PROJECT_SOURCE_DIR}/src/main.cpp)
-
-    # Building the standalone Executable
-    add_executable(${PROJECT_NAME} ${mainFile})
-
-    # Link executable with library
-    target_link_libraries(${PROJECT_NAME} PUBLIC ${PROJECT_NAME}_lib)
-endif()
+############################################################
+# Polyhedral Gravity Library & Executable & Python Interface
+############################################################
+add_subdirectory(${PROJECT_SOURCE_DIR}/src)
 
 ##############################
 # Building the Polyhedral Docs
@@ -149,11 +97,3 @@ if (BUILD_POLYHEDRAL_GRAVITY_TESTS)
     # Subdirectory where the tests are located
     add_subdirectory(${PROJECT_SOURCE_DIR}/test)
 endif ()
-
-##########################################
-# Building the Polyhedral Python Interface
-##########################################
-if(BUILD_POLYHEDRAL_GRAVITY_PYTHON_INTERFACE)
-    message(STATUS "Building the Polyhedral Gravity Python Interface")
-    add_subdirectory(${PROJECT_SOURCE_DIR}/python-binding)
-endif()
\ No newline at end of file
diff --git a/LICENSE b/LICENSE
index f288702..a97f251 100644
--- a/LICENSE
+++ b/LICENSE
@@ -672,3 +672,25 @@ may consider it more useful to permit linking proprietary applications with
 the library.  If this is what you want to do, use the GNU Lesser General
 Public License instead of this License.  But first, please read
 <https://www.gnu.org/licenses/why-not-lgpl.html>.
+
+
+##############################################################################
+-- NOTE: Third party dependency used by this software --
+##############################################################################
+This software depends on the following third party software
+and users must comply to its license terms.
+
+- spdlog (MIT License)
+  https://github.com/gabime/spdlog/blob/76dfc7e7c0d3c69d3cdaa3399b63545235ccbb02/LICENSE
+- yaml-cpp (MIT Lucense)
+  https://github.com/jbeder/yaml-cpp/blob/b8882652fcbeba4c00dec019a39da91e702e474e/LICENSE
+- thrust (Apache License Version 2.0)
+  https://github.com/NVIDIA/thrust/blob/24486a169a62a58ef8f824d3dc9613c006b6f5a7/LICENSE
+- oneTBB (Apache License Version 2.0)
+  https://github.com/oneapi-src/oneTBB/blob/37a89d47f2bca3700f9777732c41f08ac44d2d0a/LICENSE.txt
+- xsimd (BSD 3 License)
+  https://github.com/xtensor-stack/xsimd/blob/e6fa5aca6320d6ccaf24c123ab2af9b0f2f09cc1/LICENSE
+- pybind11 (BSD 3 License)
+  https://github.com/pybind/pybind11/blob/8d08dc64ca300342852ceaa7b1d65fe9f69dab06/LICENSE
+- Tetgen (GNU AGPL v3)
+  https://github.com/libigl/tetgen/blob/4f3bfba3997f20aa1f96cfaff604313a8c2c85b6/LICENSE
diff --git a/MANIFEST.in b/MANIFEST.in
new file mode 100644
index 0000000..5b8c9c0
--- /dev/null
+++ b/MANIFEST.in
@@ -0,0 +1,3 @@
+include CMakeLists.txt
+recursive-include src *
+recursive-include cmake *
\ No newline at end of file
diff --git a/README.md b/README.md
index 40db0d7..15d8633 100644
--- a/README.md
+++ b/README.md
@@ -38,13 +38,13 @@ on [readthedocs](https://polyhedral-gravity-model-cpp.readthedocs.io/en/latest/)
 The project uses the following dependencies,
 all of them are **automatically** set-up via CMake:
 
-- GoogleTest 1.11.0 (only required for testing)
-- spdlog 1.11.0 (required for logging)
-- tetgen 1.6 (required for I/O)
-- yaml-cpp 0.7.0 (required for I/O)
-- thrust 1.16.0 (required for parallelization and utility)
-- xsimd 8.1.0 (required for vectorization of the `atan(..)`)
-- pybind11 2.10.3 (required for the Python interface, but not the C++ standalone)
+- GoogleTest (1.13.0 or compatible), only required for testing
+- spdlog (1.11.0 or compatible), required for logging
+- tetgen (1.6 or compatible), required for I/O
+- yaml-cpp (0.7.0 or compatible), required for I/O
+- thrust (2.1.0 or compatible), required for parallelization and utility
+- xsimd (11.1.0 or compatible), required for vectorization of the `atan(..)`
+- pybind11 (2.10.4 or compatible), required for the Python interface, but not the C++ standalone
 
 ## Python interface
 
diff --git a/cmake/pybind11.cmake b/cmake/pybind11.cmake
index 200f2ed..774e98c 100644
--- a/cmake/pybind11.cmake
+++ b/cmake/pybind11.cmake
@@ -2,7 +2,7 @@ include(FetchContent)
 
 message(STATUS "Setting up pybind11")
 
-find_package(pybind11 2.10.3 QUIET)
+find_package(pybind11 2.10.4 QUIET)
 
 if (${pybind11_FOUND})
 
@@ -12,10 +12,10 @@ else()
 
     message(STATUS "Using pybind11 from git repository")
 
-    #Fetches the version 2.10.3 from the official github of pybind11
+    #Fetches the version 2.10.4 from the official github of pybind11
     FetchContent_Declare(pybind11
             GIT_REPOSITORY https://github.com/pybind/pybind11
-            GIT_TAG v2.10.3
+            GIT_TAG v2.10.4
             )
 
     FetchContent_MakeAvailable(pybind11)
diff --git a/cmake/thrust.cmake b/cmake/thrust.cmake
index a7ef8bf..4c6ebb8 100644
--- a/cmake/thrust.cmake
+++ b/cmake/thrust.cmake
@@ -11,7 +11,7 @@ set(THRUST_ENABLE_EXAMPLES "OFF")
 # Set standard CPP Dialect to 17 (default of thrust would be 14)
 set(THRUST_CPP_DIALECT 17)
 
-find_package(Thrust 1.16.0 QUIET)
+find_package(Thrust 2.1.0 QUIET)
 
 if (${Thrust_FOUND})
 
@@ -19,10 +19,10 @@ if (${Thrust_FOUND})
 
 else()
     message(STATUS "Using thrust from git repository")
-    # Fetches the version 1.16.0 of the official NVIDIA Thrust repository
+    # Fetches the version 2.1.0 of the official NVIDIA Thrust repository
     FetchContent_Declare(thrust
             GIT_REPOSITORY https://github.com/NVIDIA/thrust.git
-            GIT_TAG 1.16.0
+            GIT_TAG 2.1.0
             )
     FetchContent_MakeAvailable(thrust)
 endif()
diff --git a/cmake/xsimd.cmake b/cmake/xsimd.cmake
index 3318d70..dd3e820 100644
--- a/cmake/xsimd.cmake
+++ b/cmake/xsimd.cmake
@@ -3,7 +3,7 @@ include(FetchContent)
 message(STATUS "Setting up xsimd via CMake")
 
 
-find_package(xsimd 8.1.0 QUIET)
+find_package(xsimd 11.1.0 QUIET)
 
 if (${xsimd_FOUND})
 
@@ -13,10 +13,10 @@ else()
 
     message(STATUS "Using xsimd from git repository")
 
-    #Fetches the version 8.1.0 from the official github of tbb
+    #Fetches the version 11.1.0 from the official github of tbb
     FetchContent_Declare(xsimd
             GIT_REPOSITORY https://github.com/xtensor-stack/xsimd.git
-            GIT_TAG 8.1.0
+            GIT_TAG 11.1.0
             )
 
     FetchContent_MakeAvailable(xsimd)
diff --git a/docs/api/calculation.rst b/docs/api/calculation.rst
index 1c3e9a9..1ff2680 100644
--- a/docs/api/calculation.rst
+++ b/docs/api/calculation.rst
@@ -9,9 +9,19 @@ since it contains the important function :code:`evaluate(..)`
 used evaluation of the polyhedral gravity model at computation
 point(s) P.
 
+Additionally, the class :code:`GravityEvaluable` is provided as a way to
+perform the evaluation of the polyhedral gravity model repeatedly
+without the need to re-initialize the polyhedron and the gravity model as
+some caching is performed. It is recommended to use this feature when evaluating the model repeatedly on a small number of points as the overhead of setting up the models is larger than the computational cost.
+
 Further, it implements the `Möller–Trumbore intersection algorithm <https://en.wikipedia.org/wiki/Möller–Trumbore_intersection_algorithm>`__
 capable of checking if the plane unit normals of the polyhedron are pointing outwards.
 
+GravityEvaluable
+----------------
+
+.. doxygenclass:: polyhedralGravity::GravityEvaluable
+
 
 GravityModel
 ------------
@@ -20,6 +30,11 @@ GravityModel
 
 
 MeshChecking
------------
+------------
+
+.. doxygennamespace:: polyhedralGravity::MeshChecking
+
+Other
+-----
 
-.. doxygennamespace:: polyhedralGravity::MeshChecking
\ No newline at end of file
+.. doxygenfunction:: polyhedralGravity::transformPolyhedron
diff --git a/docs/api/python.rst b/docs/api/python.rst
index b773fa7..289f44a 100644
--- a/docs/api/python.rst
+++ b/docs/api/python.rst
@@ -7,67 +7,75 @@ polyhedral_gravity
 Computes the full gravity tensor for a given constant density polyhedron which consists of some
 vertices and triangular faces at given computation points
 
-.. py::module:: polyhedral_gravity
+.. py:module:: polyhedral_gravity
 
-Computes the full gravity tensor for a given constant density polyhedron which consists of some vertices and
-triangular faces at a given computation point :math:`P`
+    .. py:function:: evaluate(polyhedral_source, density, computation_points, parallel=True)
 
-.. py:function:: evaluate(vertices, faces, density, computation_point)
-   :noindex:
+        Evaluates the polyhedral gravity model for a given constant density polyhedron at a given computation point.
 
-   Evaluate the full gravity tensor for a given constant density polyhedron which consists of some vertices and
-   triangular faces at a given computation point P
+        :param Union[Tuple[List[List[float[3]]], List[List[int[3]]]], List[str]] polyhedral_source:
+            The vertices & faces of the polyhedron as tuple or
+            the filenames of the files containing the vertices & faces as list of strings
+            Supports numpy arrays!
+            Supported are .node/.face, mesh, .ply, .off, .stl files.
+        :param float density:
+            The constant density of the polyhedron in :math:`\frac{kg}{m^3}`
+        :param List[float[3]] computation_points:
+            The cartesian computation points as tuple or list of points
+        :param bool parallel:
+            If true, the computation is done in parallel (default: true)
 
-   :param List[List[float[3]]] vertices: vertices of the polyhedron
-   :param List[List[int[3]]] faces: faces of the polyhedron
-   :param float density: constant density in :math:`\frac{kg}{m^3}`
-   :param List[float[3]] computation_point: cartesian computation point :math:`P`
-   :return: tuple of potential, acceleration, and second derivative tensor
-   :rtype: Tuple[float, List[float[3]], List[float[6]]]
+        :rtype: Union[Tuple[float, List[float[3]], List[float[6]]], List[Tuple[float, List[float[3]], List[float[6]]]]]
+        :return:
+            Either a tuple of potential :math:`V`, acceleration (:math:`V_x, V_y, V_z`)
+            and second derivatives (:math:`V_{xx}, V_{yy}, V_{zz}, V_{xy}, V_{xz}, V_{yz}`)
+            at the computation points or if multiple computation points are given, the result is a list of tuples
 
 
-.. py:function:: evaluate(vertices, faces, density, computation_points)
-   :noindex:
+    .. py:class:: GravityEvaluable
 
-   Evaluate the full gravity tensor for a given constant density polyhedron which consists of some vertices and
-   triangular faces at multiple given computation points
+        A class to evaluate the polyhedral gravity model for a given constant density polyhedron at a given computation point.
+        It provides a :code:`__call__` method to evaluate the polyhedral gravity model for computation points while
+        also caching the polyhedron data over the lifetime of the object.
 
-   :param List[List[float[3]]] vertices: vertices of the polyhedron
-   :param List[List[int[3]]] faces: faces of the polyhedron
-   :param float density: constant density in :math:`\frac{kg}{m^3}`
-   :param List[List[float[3]]] computation_points: multiple cartesian computation points :math:`P`
-   :return: list of tuple of potential, acceleration, and second derivative tensor
-   :rtype: List[Tuple[float, List[float[3]], List[float[6]]]]
+        .. py:method:: __init__(polyhedral_source, density)
 
-.. py:function:: evaluate(input_files, density, computation_point)
-   :noindex:
+            :param Union[Tuple[List[List[float[3]]], List[List[int[3]]]], List[str]] polyhedral_source:
+                The vertices & faces as (N, 3) array-like of the polyhedron as tuple or
+                the filenames of the files containing the vertices & faces as list of strings
+                Supports numpy arrays!
+                Supported are .node/.face, mesh, .ply, .off, .stl files.
+            :param float density:
+                The constant density of the polyhedron in :math:`\frac{kg}{m^3}`
 
-   Evaluate the full gravity tensor for a given constant density polyhedron which consists of some vertices and
-   triangular faces at a given computation point :math:`P`
+        .. py:method:: __call__(computation_points, parallel=True)
 
-   :param List[str] input_files: polyhedral source files
-   :param float density: constant density in :math:`\frac{kg}{m^3}`
-   :param List[float[3]] computation_point: cartesian computation point :math:`P`
-   :return: tuple of potential, acceleration, and second derivative tensor
-   :rtype: Tuple[float, List[float[3]], List[float[6]]]
+            Evaluates the polyhedral gravity model for a given constant density polyhedron at a given computation point.
 
-.. py:function:: evaluate(input_files, density, computation_points)
-   :noindex:
+            :param List[float[3]] computation_points:
+                The cartesian computation point
+            :param bool parallel:
+                If true, the computation is done in parallel (default: true)
 
-    Evaluate the full gravity tensor for a given constant density polyhedron which consists of some vertices and
-    triangular faces at multiple given computation points
+            :rtype: Tuple[float, List[float[3]], List[float[6]]]
+            :return:
+                A tuple of potential :math:`V`, acceleration (:math:`V_x, V_y, V_z`)
+                and second derivatives (:math:`V_{xx}, V_{yy}, V_{zz}, V_{xy}, V_{xz}, V_{yz}`)
+                at the computation point.
 
-   :param List[str] input_files: polyhedral source files
-   :param List[List[int[3]]] faces: faces of the polyhedron
-   :param float density: constant density in :math:`\frac{kg}{m^3}`
-   :param List[List[float[3]]] computation_points: multiple cartesian computation points :math:`P`
-   :return: list of tuple of potential, acceleration, and second derivative tensor
-   :rtype: List[Tuple[float, List[float[3]], List[float[6]]]]
+        .. py:method:: __repr__()
+
+            Returns a string representation of the polyhedral gravity evaluable.
+
+            :rtype: str
+            :return: A string representation of the polyhedron
 
 
 utility
 ~~~~~~~
 
+.. py:module:: polyhedral_gravity.utility
+
 This submodule contains useful utility functions like parsing meshes
 or checking if the polyhedron's mesh plane unit normals point outwards
 like it is required by the polyhedral-gravity model.
diff --git a/docs/build.rst b/docs/build.rst
index c002ecf..ead5b0c 100644
--- a/docs/build.rst
+++ b/docs/build.rst
@@ -6,13 +6,13 @@ Requirements
 
 Requirements (all of them are automatically set-up via :code:`CMake`):
 
-- GoogleTest 1.11.0 (only required for testing)
-- spdlog 1.11.0 (required for logging)
-- tetgen 1.6 (required for I/O)
-- yaml-cpp 0.7.0 (required for I/O)
-- thrust 1.16.0 (required for parallelization and utility)
-- xsimd 8.1.0 (required for vectorization of the `atan(..)`)
-- pybind11 2.10.3 (required for the Python interface, but not the C++ standalone)
+- GoogleTest (1.13.0 or compatible), only required for testing
+- spdlog (1.11.0 or compatible), required for logging
+- tetgen (1.6 or compatible), required for I/O
+- yaml-cpp (0.7.0 or compatible), required for I/O
+- thrust (2.1.0 or compatible), required for parallelization and utility
+- xsimd (11.1.0 or compatible), required for vectorization of the :code:`atan(..)`
+- pybind11 (2.10.4 or compatible), required for the Python interface, but not the C++ standalone
 
 
 Build with CMake
@@ -54,7 +54,12 @@ The module will be build using CMake. Just execute in repository root:
     pip install .
 
 To modify the build options (like parallelization) have a look
-at the :code:`setupy.py`.
+at the :code:`setupy.py`. As simple example, to modify the parallelization,
+just set the environment variable like below:
+
+.. code-block::
+
+    export POLYHEDRAL_GRAVITY_PARALLELIZATION="TBB"
 
 Installation with conda
 -----------------------
diff --git a/docs/quick_start_cpp.rst b/docs/quick_start_cpp.rst
index 2860319..b4b1feb 100644
--- a/docs/quick_start_cpp.rst
+++ b/docs/quick_start_cpp.rst
@@ -53,38 +53,44 @@ The following paragraph gives some examples on how to
 use the polyhedral model as library from within source code
 (All examples with :code:`using namespace polyhedralGravity`).
 
+Individual Function (without caching)
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
 
 **Example 1:** Evaluating the gravity model for a given polyhedron
 defined from within source code for a specific point and density.
+Further, we disable the parallelization using the optional fourth parameter (which defaults to true).
 
 .. code-block:: cpp
 
         // Defining every input parameter in the source code
         std::vector<std::array<double, 3>> vertices = ...
         std::vector<std::array<size_t, 3>> faces = ...
+        Polyhedron polyhedron{vertices, faces};
         double density = ...
         std::array<double, 3> point = ...
 
-        // Main method, notice that the last argument can also be a
-        // vector of points
-        GravityResult result = GravityModel::evaluate({vertices, faces}, density, point);
+        // Returns either a single of vector of results
+        // Here, we only have one point. Thus we get a single result
+        const auto[pot, acc, tensor] = GravityModel::evaluate(polyhedron, density, point, false);
 
 
 **Example 2:** Evaluating the gravity model for a given polyhedron
 in some source files for a specific point and density.
+Further, we explicitly enable the parallelization using the optional fourth parameter
+(which defaults to true).
 
 .. code-block:: cpp
 
         // Reading just the polyhedral source from file,
         // whereas the rest is defined within source-code
-        TetgenAdapter input{{"tsoulis.node", "tsoulis.face"}};
-        Polyhedron poly = input.getPolyhedron();
+        auto files = std::vector<std::string>{"tsoulis.node", "tsoulis.face"};
         double density = ...
         std::array<double, 3> point = ...
 
-        // Main method, notice that the last argument can also be a
-        // vector of points
-        GravityResult result = GravityModel::evaluate(poly, density, point);
+        // Returns either a single of vector of results
+        // Here, we only have one point. Thus we get a single result
+        const auto[pot, acc, tensor] = GravityModel::evaluate(files, density, point, true);
 
 
 **Example 3:** Evaluating the gravity model for a given configuration
@@ -96,11 +102,12 @@ from a .yaml file.
         std::shared_ptr<ConfigSource> config = std::make_shared<YAMLConfigReader>("config.yaml");
         Polyhedron poly = config->getDataSource()->getPolyhedron();
         double density = config->getDensity();
-        std::array<double, 3> point = config->getPointsOfInterest()[0];
+        // This time, we use multiple points
+        std::vector<std::array<double, 3>> points = config->getPointsOfInterest();
 
-        // Main method, notice that the last argument can also be a
-        // vector of points
-        GravityResult result = GravityModel::evaluate(poly, density, point);
+        // Returns either a single of vector of results
+        // Here, we have multiple point. Thus we get a vector of results!
+        const results = GravityModel::evaluate(poly, density, points);
 
 **Example 4:** A guard statement checks that the plane unit
 normals are pointing outwards and no triangle is degenerated.
@@ -120,3 +127,30 @@ about the vertices' ordering due to its quadratic complexity!
             GravityResult result = GravityModel::evaluate(poly, density, point);
         }
 
+
+
+Evaluation with caching
+~~~~~~~~~~~~~~~~~~~~~~~
+
+Evaluating the gravity model for a given polyhedron
+defined from within source code for a specific point and density.
+
+.. code-block:: cpp
+
+        // Defining every input parameter in the source code
+        std::vector<std::array<double, 3>> vertices = ...
+        std::vector<std::array<size_t, 3>> faces = ...
+        Polyhedron polyhedron{vertices, faces};
+        double density = ...
+        std::array<double, 3> point = ...
+        std::vector<std::array<double, 3>> points = ...
+
+        // Instantiation of the GravityEvaluable object
+        GravityEvaluable evaluable{polyhedron, density};
+
+        // From now, we can evaluate the gravity model for any point with
+        const auto[pot, acc, tensor] = evaluable(point);
+        // or for multiple points with
+        const auto results = evaluable(points);
+        // and we can also disable e.g. the parallelization like for the free function
+        const auto singleResultTuple = evaluable(point, false);
diff --git a/docs/quick_start_python.rst b/docs/quick_start_python.rst
index fa3edc9..9eb96e7 100644
--- a/docs/quick_start_python.rst
+++ b/docs/quick_start_python.rst
@@ -2,7 +2,47 @@ Quick Start Python
 ==================
 
 The use of the Python interface is pretty straight-forward since
-there is only one method: :code:`evaluate(..)`.
+there is only one method: :code:`evaluate(..)` or the alternative
+class :code:`GravityEvaluable` caching the polyhedron data.
+
+The method calls follow the same pattern as the C++ interface. Thus it is always:
+
+.. code-block:: python
+
+    import polyhedral_gravity as model
+
+    #############################################
+    # Free function call
+    #############################################
+    # Define the polyhedral_source which is either (vertices, faces) or a list of files
+    # Define the density (a float)
+    # Define the computation_points which is either a single point or a list of points
+    # Define if the computation is parallel or not (the default is parallel, which corresponds to True)
+    results = model.evaluate(
+        polyhedral_source=polyhedral_source,
+        density=density
+        computation_points=computation_points
+        parallel=True
+    )
+
+    #############################################
+    # With the evaluable class
+    # This allows multiple evaluations on the same polyhedron
+    # without the overhead of setting up the normals etc.
+    #############################################
+    # Parameters are the same as for the free function call
+    evaluable = model.GravityEvaluable(
+        polyhedral_source=polyhedral_source,
+        density=density
+    )
+    results = evaluable(
+        computation_points=computation_points,
+        parallel=True
+    )
+
+
+Free function
+-------------
 
 **Example 1:** Evaluating the gravity model for a given polyhedron
 defined from within source code for a specific point and density.
@@ -18,13 +58,14 @@ defined from within source code for a specific point and density.
         computation_point = ... # [] or np.array of length 3
 
         # Evaluate the gravity model
-        # Notice that the last argument could also be a list of points
-        # Returns (list of) tuple of potential, acceleration and tensor
-        potential, acceleration, tensor = model.evaluate(vertices, faces, density, computation_point)
+        # Notice that the third argument could also be a list of points
+        # Returns a tuple of potential, acceleration and tensor
+        potential, acceleration, tensor = model.evaluate((vertices, faces), density, computation_point, parallel=True)
 
 
 **Example 2a:** Evaluating the gravity model for a given polyhedron
-in some source files for a specific point and density.
+in some source files for several points and density.
+Of course, you can also use keyword arguments for the parameters.
 
 .. code-block:: python
 
@@ -34,12 +75,17 @@ in some source files for a specific point and density.
         file_vertices = '___.node'      # str, path to file
         file_nodes = '___.face'         # str, path to file
         density = ...                   # float
-        computation_point = ...         # [] or np.array of length 3
+        computation_points = ...        # [] or np.array of length 3
 
         # Evaluate the gravity model
         # Notice that the last argument could also be a list of points
-        # Returns (list of) tuple of potential, acceleration and tensor
-        potential, acceleration, tensor = model.evaluate([file_vertices, file_nodes], density, computation_point)
+        # Returns a list of tuple of potential, acceleration and tensor
+        results = model.evaluate(
+            polyhedral_source=[file_vertices, file_nodes],
+            density=density,
+            computation_points=computation_points
+            parallel=True
+        )
 
 
 **Example 2b:** Evaluating the gravity model for a given polyhedron
@@ -56,7 +102,7 @@ in some source files for a specific point and density.
 
         # Evaluate the gravity model
         # Notice that the last argument could also be a list of points
-        # Returns (list of) tuple of potential, acceleration and tensor
+        # Returns a tuple of potential, acceleration and tensor
         potential, acceleration, tensor = model.evaluate([mesh], density, computation_point)
 
 
@@ -85,5 +131,49 @@ about the vertices' ordering due to its quadratic complexity!
 
     # Additional guard statement
     if mesh_sanity.check_mesh(vertices, faces):
-        potential, acceleration, tensor = model.evaluate(vertices, faces, density, computation_point)
+        potential, acceleration, tensor = model.evaluate((vertices, faces), density, computation_point)
+
 
+Evaluable class
+---------------
+
+Use the :code:`GravityEvaluable` class to cache the polyhedron data over multiple calls.
+This drastically improves the performance, as the polyhedral data is "stored" on the C++ side,
+rather than being converted from Python to C++ for every call.
+
+This approach is especially useful one wants to calculate multiple points for the same polyhedron, but
+the points are not known in advance (e.g. when propagating a spacecraft).
+Have a look at the example below to see how to use the :code:`GravityEvaluable` class.
+
+.. code-block:: python
+
+        import polyhedral_gravity as model
+
+        # Defining every input parameter in the source code
+        vertices = ...           # [] of [] or np.array of length 3 and type float
+        faces = ...              # [] of [] or np.array of length 3 and type int
+        density = ...            # float
+        computation_points = ... # [] or np.array of length 3
+
+        # Create the evaluable object
+        evaluable = model.GravityEvaluable(polyhedral_source, density)
+
+        for point in computation_points:
+            # Evaluate the gravity model
+            potential, acceleration, tensor = evaluable(point, parallel=True)
+
+Below is a comparison of the performance of the free function and the evaluable class.
+Basically, as soon as you have more than one point to evaluate, the evaluable class is faster and
+thus recommended.
+
++----------------------------------------+-------------------------------+
+| Test                                   | Time Per Point (microseconds) |
++========================================+===============================+
+| Free Function (1000 times 1 point)     | 7765.073                      |
++----------------------------------------+-------------------------------+
+| Free Function (1 time N points)        | 275.917                       |
++----------------------------------------+-------------------------------+
+| GravityEvaluable (1000 times 1 Point)  | 313.408                       |
++----------------------------------------+-------------------------------+
+| GravityEvaluable (1 time 1000 Points)  | 253.031                       |
++----------------------------------------+-------------------------------+
diff --git a/python-binding/python-interface.cpp b/python-binding/python-interface.cpp
deleted file mode 100644
index 8f6dd9b..0000000
--- a/python-binding/python-interface.cpp
+++ /dev/null
@@ -1,187 +0,0 @@
-#include <tuple>
-#include "pybind11/pybind11.h"
-#include "pybind11/stl.h"
-
-#include "polyhedralGravity/model/Polyhedron.h"
-#include "polyhedralGravity/model/GravityModelData.h"
-#include "polyhedralGravity/calculation/GravityModel.h"
-#include "polyhedralGravity/calculation/MeshChecking.h"
-#include "polyhedralGravity/input/TetgenAdapter.h"
-
-
-namespace py = pybind11;
-
-/**
- * This method converts a named tuples GravityModelResult to a tuple
- * @param result a named tuple of type GravityModelResult
- * @return tuple of potential, acceleration and gradiometric tensor
- */
-std::tuple<double, std::array<double, 3>, std::array<double, 6>> convertToTuple(
-        const polyhedralGravity::GravityModelResult &result) {
-    return std::make_tuple(result.gravitationalPotential,
-                           result.acceleration,
-                           result.gradiometricTensor);
-}
-
-/**
- * This method converts a vector of named tuples GravityModelResult to a vector of tuples
- * @param resultVector a vector of named tuples of type GravityModelResult
- * @return vector of tuples of potential, acceleration and gradiometric tensor
- */
-std::vector<std::tuple<double, std::array<double, 3>, std::array<double, 6>>> convertToTupleVector(
-        const std::vector<polyhedralGravity::GravityModelResult> &resultVector) {
-    std::vector<std::tuple<double, std::array<double, 3>, std::array<double, 6>>> resVectorTuple{resultVector.size()};
-    std::transform(resultVector.cbegin(), resultVector.cend(), resVectorTuple.begin(), &convertToTuple);
-    return resVectorTuple;
-}
-
-PYBIND11_MODULE(polyhedral_gravity, m) {
-    using namespace polyhedralGravity;
-    m.doc() = "Computes the full gravity tensor for a given constant density polyhedron which consists of some "
-              "vertices and triangular faces at a given computation points";
-
-    /*
-     * Methods for vertices and faces from the script context
-     */
-
-    m.def("evaluate",
-          [](const std::vector<std::array<double, 3>> &vertices, const std::vector<std::array<size_t, 3>> &faces,
-             double density, const std::array<double, 3> &computationPoint) {
-              return convertToTuple(GravityModel::evaluate({vertices, faces}, density, computationPoint));
-          }, R"mydelimiter(
-            Evaluate the full gravity tensor for a given constant density polyhedron which consists of some vertices
-            and triangular faces at a given computation point P.
-
-            Args:
-                vertices (2-D array-like): (N, 3) array of vertex coordinates (floats).
-                faces (2-D array-like): (N, 3) array of faces, vertex-indices (ints).
-                density (float): the constant density of the polyhedron.
-                computation_point (2-D array-like): (N, 3) cartesian points.
-            Returns:
-                tuple of potential, acceleration and second derivative gravity tensor.
-
-          )mydelimiter",
-          py::arg("vertices"), py::arg("faces"), py::arg("density"), py::arg("computation_point"));
-
-    m.def("evaluate",
-          [](const std::vector<std::array<double, 3>> &vertices, const std::vector<std::array<size_t, 3>> &faces,
-             double density, const std::vector<std::array<double, 3>> &computationPoints) {
-              return convertToTupleVector(GravityModel::evaluate({vertices, faces}, density, computationPoints));
-          }, R"mydelimiter(
-            Evaluate the full gravity tensor for a given constant density polyhedron which consists of some vertices
-            and triangular faces at multiple given computation points.
-
-            Args:
-                vertices (2-D array-like): (N, 3) array of vertex coordinates (floats).
-                faces (2-D array-like): (N, 3) array of faces, vertex-indices (ints).
-                density (float): the constant density of the polyhedron.
-                computation_point (array-like): cartesian point.
-            Returns:
-                array of tuples of potential, acceleration and second derivative gravity tensor.
-
-          )mydelimiter",
-          py::arg("vertices"), py::arg("faces"), py::arg("density"), py::arg("computation_points"));
-
-    /*
-     * Methods for vertices and faces from files via TetgenAdapter
-     */
-
-    m.def("evaluate",
-          [](const std::vector<std::string> &filenames,
-             double density, const std::array<double, 3> &computationPoint) {
-              TetgenAdapter tetgen{filenames};
-              return convertToTuple(GravityModel::evaluate(tetgen.getPolyhedron(), density, computationPoint));
-          }, R"mydelimiter(
-            Evaluate the full gravity tensor for a given constant density polyhedron which consists of some vertices
-            and triangular faces at a given computation point P. The vertices and faces are read from input
-            files (either .node/.face, mesh, .ply, .off, .stl). File-Order matters in case of the first option!
-
-            Args:
-                input_files (List[str]): list of input files.
-                density (float): the constant density of the polyhedron.
-                computation_point (array-like): cartesian point.
-            Returns:
-                tuple of potential, acceleration and second derivative gravity tensor.
-
-          )mydelimiter",
-          py::arg("input_files"), py::arg("density"), py::arg("computation_point"));
-
-    m.def("evaluate",
-          [](const std::vector<std::string> &filenames,
-             double density, const std::vector<std::array<double, 3>> &computationPoints) {
-              TetgenAdapter tetgen{filenames};
-              return convertToTupleVector(GravityModel::evaluate(tetgen.getPolyhedron(), density, computationPoints));
-          }, R"mydelimiter(
-            Evaluate the full gravity tensor for a given constant density polyhedron which consists of some vertices
-            and triangular faces multiple given computation points. The vertices and faces are read from input
-            files (either .node/.face, mesh, .ply, .off, .stl). File-Order matters in case of the first option!
-
-            Args:
-                input_files (List[str]): list of input files.
-                density (float): the constant density of the polyhedron.
-                computation_point (2-D array-like): (N, 3) cartesian points.
-            Returns:
-                array of tuples of potential, acceleration and second derivative gravity tensor.
-
-          )mydelimiter",
-          py::arg("input_files"), py::arg("density"), py::arg("computation_points"));
-
-    py::module_ utility = m.def_submodule("utility",
-                                          "This submodule contains useful utility functions like parsing meshes "
-                                          "or checking if the polyhedron's mesh plane unit normals point outwards "
-                                          "like it is required by the polyhedral-gravity model.");
-
-    utility.def("read",
-                [](const std::vector<std::string> &filenames) {
-                    TetgenAdapter tetgen{filenames};
-                    auto polyhedron = tetgen.getPolyhedron();
-                    return std::make_tuple(polyhedron.getVertices(), polyhedron.getFaces());
-                }, R"mydelimiter(
-            Reads a polyhedron from a mesh file. The vertices and faces are read from input
-            files (either .node/.face, mesh, .ply, .off, .stl). File-Order matters in case of the first option!
-
-            Args:
-                input_files (List[str]): list of input files.
-            Returns:
-                tuple of vertices (N, 3) (floats) and faces (N, 3) (ints).
-          )mydelimiter",
-                py::arg("input_files"));
-
-    utility.def("check_mesh",
-                [](const std::vector<std::array<double, 3>> &vertices,
-                   const std::vector<std::array<size_t, 3>> &faces) {
-                    Polyhedron poly {vertices, faces};
-                    return MeshChecking::checkTrianglesNotDegenerated(poly) && MeshChecking::checkNormalsOutwardPointing(poly);
-                }, R"mydelimiter(
-                Checks if no triangles of the polyhedral mesh are degenerated by checking that their surface area
-                is greater zero.
-                Further, Checks if all the polyhedron's plane unit normals are pointing outwards.
-
-                Args:
-                    vertices (2-D array-like): (N, 3) array of vertex coordinates (floats).
-                    faces (2-D array-like): (N, 3) array of faces, vertex-indices (ints).
-                Returns:
-                    True if no triangle is degenerate and the polyhedron's plane unit normals are all pointing outwards.
-          )mydelimiter",
-                py::arg("vertices"), py::arg("faces"));
-
-
-    utility.def("check_mesh",
-                [](const std::vector<std::string> &filenames) {
-                    TetgenAdapter tetgen{filenames};
-                    Polyhedron poly = tetgen.getPolyhedron();
-                    return MeshChecking::checkTrianglesNotDegenerated(poly) && MeshChecking::checkNormalsOutwardPointing(poly);
-                }, R"mydelimiter(
-                Checks if no triangles of the polyhedral mesh are degenerated by checking that their surface area
-                is greater zero.
-                Further, Checks if all the polyhedron's plane unit normals are pointing outwards.
-                Reads a polyhedron from a mesh file. The vertices and faces are read from input
-                files (either .node/.face, mesh, .ply, .off, .stl). File-Order matters in case of the first option!
-
-                Args:
-                    input_files (List[str]): list of input files.
-                Returns:
-                    True if no triangle is degenerate and the polyhedron's plane unit normals are all pointing outwards.
-          )mydelimiter",
-                py::arg("input_files"));
-}
diff --git a/script/time.py b/script/time.py
new file mode 100755
index 0000000..23a045c
--- /dev/null
+++ b/script/time.py
@@ -0,0 +1,69 @@
+#!python3
+import polyhedral_gravity
+import numpy as np
+import matplotlib.pyplot as plt
+import mesh_utility
+import timeit
+
+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)
+
+# Generate 1000 random cartesian points
+N = 1000
+computation_points = np.random.uniform(-2, 2, (N, 3))
+
+DENSITY = 1.0
+
+########################################################################################################################
+# OLD
+########################################################################################################################
+start_time = timeit.default_timer()
+for i in range(N):
+    polyhedral_gravity.evaluate((vertices, faces), DENSITY, computation_points[i])
+end_time = timeit.default_timer()
+
+delta = (end_time - start_time) / N * 1e6
+# Print the time in milliseconds
+print("######## Old Single-Point ########")
+print(f"--> {N} times 1 point with old interface")
+print(f"--> Time taken: {delta:.3f} microseconds per point")
+
+start_time = timeit.default_timer()
+polyhedral_gravity.evaluate((vertices, faces), DENSITY, computation_points)
+end_time = timeit.default_timer()
+
+delta = (end_time - start_time) / N * 1e6
+# Print the time in milliseconds
+print("######## Old Multi-Point ########")
+print("--> 1 time N points with old interface")
+print(f"--> Time taken: {delta:.3f} microseconds per point")
+
+########################################################################################################################
+# NEW
+########################################################################################################################
+evaluable = polyhedral_gravity.GravityEvaluable((vertices, faces), DENSITY)
+
+start_time = timeit.default_timer()
+for i in range(N):
+    evaluable(computation_points[i])
+end_time = timeit.default_timer()
+
+delta = (end_time - start_time) / N * 1e6
+# Print the time in milliseconds
+print("######## GravityEvaluable (Single-Point) #########")
+print(f"--> {N} times 1 point with GravityEvaluable")
+print(f"--> Time taken: {delta:.3f} microseconds per point")
+
+evaluable = polyhedral_gravity.GravityEvaluable((vertices, faces), DENSITY)
+
+start_time = timeit.default_timer()
+evaluable(computation_points)
+end_time = timeit.default_timer()
+
+delta = (end_time - start_time) / N * 1e6
+# Print the time in milliseconds
+print("######## GravityEvaluable (Multi-Point) ########")
+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 bd8c3c7..980b820 100644
--- a/setup.py
+++ b/setup.py
@@ -16,8 +16,8 @@
 CMAKE_OPTIONS = {
     # The Build Type (Should be release!)
     "CMAKE_BUILD_TYPE": "Release",
-    # Modify to change the parallelization (Default value: CPP)
-    "POLYHEDRAL_GRAVITY_PARALLELIZATION": "CPP",
+    # Modify to change the parallelization (Default value: TBB)
+    "POLYHEDRAL_GRAVITY_PARALLELIZATION": "TBB",
     # Default value (INFO=2)
     "LOGGING_LEVEL": 2,
     # Default value (OFF)
@@ -32,8 +32,8 @@
     "BUILD_POLYHEDRAL_GRAVITY_EXECUTABLE": "OFF",
     # Should be of course ON!
     "BUILD_POLYHEDRAL_GRAVITY_PYTHON_INTERFACE": "ON",
-    # Build shared libs by default
-    "BUILD_SHARED_LIBS": "ON"
+    # Build static libs by default (On conda-forge we build shared libs by setting this to ON)
+    "BUILD_SHARED_LIBS": "OFF"
 }
 # ---------------------------------------------------------------------------------
 
@@ -141,7 +141,7 @@ def build_extension(self, ext):
 # --------------------------------------------------------------------------------
 setup(
     name="polyhedral_gravity",
-    version="1.2.2",
+    version="2.0",
     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",
@@ -152,7 +152,9 @@ def build_extension(self, ext):
     """,
     ext_modules=[CMakeExtension("polyhedral_gravity")],
     cmdclass={"build_ext": CMakeBuild},
+    license_file="LICENSE",
     zip_safe=False,
-    python_requires=">=3.6"
+    python_requires=">=3.6",
+    include_package_data=True,
 )
 # --------------------------------------------------------------------------------
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
new file mode 100644
index 0000000..5875530
--- /dev/null
+++ b/src/CMakeLists.txt
@@ -0,0 +1,67 @@
+###############################
+# Building the Polyhedral Libray
+################################
+if (BUILD_POLYHEDRAL_GRAVITY_LIBRARY)
+    message(STATUS "Building the Polyhedral Gravity Library")
+
+    # YAML is only required for the library target, but not for e.g. the Python interface
+    include(yaml)
+
+    add_subdirectory(polyhedralGravity)
+
+    # Adds the include Path PROJECT_SOURCE_DIR/src to the target polyhedralGravity
+    target_include_directories(${PROJECT_NAME}_lib PUBLIC
+            "${CMAKE_CURRENT_SOURCE_DIR}"
+            )
+
+    # Link libraries needed in all targets
+    target_link_libraries(${PROJECT_NAME}_lib
+            spdlog::spdlog
+            yaml-cpp
+            tetgen
+            xsimd
+            Thrust
+            )
+endif ()
+#####################################
+# Building the Polyhedral Executable
+#####################################
+if (BUILD_POLYHEDRAL_GRAVITY_EXECUTABLE)
+    message(STATUS "Building the Polyhedral Gravity Executable")
+
+    # Building the standalone Executable
+    add_executable(${PROJECT_NAME} main.cpp)
+
+    # Link executable with library
+    target_link_libraries(${PROJECT_NAME} PUBLIC ${PROJECT_NAME}_lib)
+
+    # Place the executable in the top level directory
+    set_target_properties(${PROJECT_NAME} PROPERTIES
+            ARCHIVE_OUTPUT_DIRECTORY ${PROJECT_BINARY_DIR}
+            LIBRARY_OUTPUT_DIRECTORY ${PROJECT_BINARY_DIR}
+            RUNTIME_OUTPUT_DIRECTORY ${PROJECT_BINARY_DIR}
+            )
+endif()
+
+##########################################
+# Building the Polyhedral Python Interface
+##########################################
+if(BUILD_POLYHEDRAL_GRAVITY_PYTHON_INTERFACE)
+    message(STATUS "Building the Polyhedral Gravity Python Interface")
+
+    # Include pybind11
+    include(pybind11)
+
+    add_subdirectory(polyhedralGravityPython)
+
+    target_include_directories(polyhedral_gravity PUBLIC
+            "${CMAKE_CURRENT_SOURCE_DIR}"
+            )
+
+    target_link_libraries(polyhedral_gravity PUBLIC
+            spdlog::spdlog
+            tetgen
+            xsimd
+            Thrust
+            )
+endif()
\ No newline at end of file
diff --git a/src/main.cpp b/src/main.cpp
index 99d2bde..42ef01a 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -10,8 +10,8 @@ int main(int argc, char *argv[]) {
     using namespace polyhedralGravity;
     if (argc != 2) {
         SPDLOG_LOGGER_INFO(PolyhedralGravityLogger::DEFAULT_LOGGER.getLogger(), "Wrong program call! "
-                                                                                 "Please use the program like that:\n"
-                                                                                 "./polyhedralGravity [YAML-Configuration-File]\n");
+                                                                                "Please use the program like this:\n"
+                                                                                "./polyhedralGravity [YAML-Configuration-File]\n");
         return 0;
     }
 
diff --git a/src/polyhedralGravity/CMakeLists.txt b/src/polyhedralGravity/CMakeLists.txt
new file mode 100644
index 0000000..efd3fb1
--- /dev/null
+++ b/src/polyhedralGravity/CMakeLists.txt
@@ -0,0 +1,6 @@
+# Including the actual Polyhedral source files
+file(GLOB_RECURSE SRC
+        "${CMAKE_CURRENT_SOURCE_DIR}/*.h"
+        "${CMAKE_CURRENT_SOURCE_DIR}/*.cpp")
+
+add_library(${PROJECT_NAME}_lib OBJECT ${SRC})
diff --git a/src/polyhedralGravity/calculation/GravityEvaluable.cpp b/src/polyhedralGravity/calculation/GravityEvaluable.cpp
new file mode 100644
index 0000000..be399ed
--- /dev/null
+++ b/src/polyhedralGravity/calculation/GravityEvaluable.cpp
@@ -0,0 +1,238 @@
+#include "GravityEvaluable.h"
+
+
+namespace polyhedralGravity {
+
+    void GravityEvaluable::prepare() const {
+        using namespace GravityModel::detail;
+        // Initialize the vectors and allocate the required memory
+        const size_t n = _polyhedron.countFaces();
+        const auto &vertices = _polyhedron.getVertices();
+        const auto &faces = _polyhedron.getFaces();
+        _segmentVectors.resize(n);
+        _planeUnitNormals.resize(n);
+        _segmentUnitNormals.resize(n);
+
+        // Create the iterators for the for_each loop over the polyhedral faces
+        thrust::counting_iterator<size_t> begin{0};
+        thrust::counting_iterator<size_t> end{n};
+
+        // Compute the segment vectors, the plane unit normals and the segment unit normals
+        thrust::for_each(thrust::device, begin, end, [this, &faces, &vertices](size_t index) {
+            Array3Triplet face{vertices[faces[index][0]], vertices[faces[index][1]], vertices[faces[index][2]]};
+            //1-01 Step: Compute Segment Vectors G_pq which describe each one the edge between two vertices
+            _segmentVectors[index] = buildVectorsOfSegments(face[0], face[1], face[2]);
+            //1-02 Step: Compute the Plane Unit Normals N_p (pointing outside the polyhedron)
+            _planeUnitNormals[index] = buildUnitNormalOfPlane(_segmentVectors[index][0], _segmentVectors[index][1]);
+            //1-03 Step: Compute Segment Unit Normals n_pq (normal pointing away from each segment)
+            _segmentUnitNormals[index] = buildUnitNormalOfSegments(_segmentVectors[index], _planeUnitNormals[index]);
+        });
+    }
+
+    template<bool Parallelization>
+    GravityModelResult GravityEvaluable::evaluate(const Array3 &computationPoint) const {
+        using namespace GravityModel::detail;
+        using namespace util;
+        SPDLOG_LOGGER_DEBUG(PolyhedralGravityLogger::DEFAULT_LOGGER.getLogger(),
+                            "Evaluation for computation point P = [{}, {}, {}] started, given density = {} kg/m^3",
+                            computationPoint[0], computationPoint[1], computationPoint[2], _density);
+        /*
+         * Calculate V and Vx, Vy, Vz and Vxx, Vyy, Vzz, Vxy, Vxz, Vyz
+         */
+        auto polyhedronIterator = transformPolyhedron(_polyhedron, computationPoint);
+        auto zip1 = util::zip(polyhedronIterator.first, _segmentVectors.begin(), _planeUnitNormals.begin(),
+                              _segmentUnitNormals.begin());
+        auto zip2 = util::zip(polyhedronIterator.second, _segmentVectors.end(), _planeUnitNormals.end(),
+                              _segmentUnitNormals.end());
+
+        SPDLOG_LOGGER_DEBUG(PolyhedralGravityLogger::DEFAULT_LOGGER.getLogger(),
+                            "Starting to iterate over the planes...");
+        GravityModelResult result{};
+        auto &[potential, acceleration, gradiometricTensor] = result;
+
+        if constexpr (Parallelization) {
+            result = thrust::transform_reduce(thrust::device, zip1, zip2, &GravityEvaluable::evaluateFace, result,
+                                              util::operator+<double, Array3, Array6>);
+        } else {
+            result = thrust::transform_reduce(thrust::host, zip1, zip2, &GravityEvaluable::evaluateFace, result,
+                                              util::operator+<double, Array3, Array6>);
+        }
+
+        SPDLOG_LOGGER_DEBUG(PolyhedralGravityLogger::DEFAULT_LOGGER.getLogger(),
+                            "Finished the sums. Applying final prefix and eliminating rounding errors.");
+
+        //9. Step: Compute prefix consisting of GRAVITATIONAL_CONSTANT * density
+        const double prefix = util::GRAVITATIONAL_CONSTANT * _density;
+
+        //10. Step: Final expressions after application of the prefix (and a division by 2 for the potential)
+        potential = (potential * prefix) / 2.0;
+        acceleration = acceleration * prefix;
+        gradiometricTensor = gradiometricTensor * prefix;
+        return result;
+    }
+
+    // Explicit template instantiation of the single point evaluate method
+    template GravityModelResult GravityEvaluable::evaluate<true>(const Array3 &computationPoints) const;
+
+    template GravityModelResult GravityEvaluable::evaluate<false>(const Array3 &computationPoints) const;
+
+    template<bool Parallelization>
+    std::vector<GravityModelResult> GravityEvaluable::evaluate(const std::vector<Array3> &computationPoints) const {
+        std::vector<GravityModelResult> result{computationPoints.size()};
+        if constexpr (Parallelization) {
+            thrust::transform(thrust::device, computationPoints.begin(), computationPoints.end(), result.begin(),
+                              [this](const Array3 &computationPoint) {
+                                  return this->evaluate<false>(computationPoint);
+                              });
+        } else {
+            thrust::transform(thrust::host, computationPoints.begin(), computationPoints.end(), result.begin(),
+                              [this](const Array3 &computationPoint) {
+                                  return this->evaluate<false>(computationPoint);
+                              });
+        }
+        return result;
+    }
+
+    // Explicit template instantiation of the multipoint evaluate method
+    template std::vector<GravityModelResult>
+    GravityEvaluable::evaluate<true>(const std::vector<Array3> &computationPoints) const;
+
+    template std::vector<GravityModelResult>
+    GravityEvaluable::evaluate<false>(const std::vector<Array3> &computationPoints) const;
+
+    GravityModelResult
+    GravityEvaluable::evaluateFace(const thrust::tuple<Array3Triplet, Array3Triplet, Array3, Array3Triplet> &tuple) {
+        using namespace util;
+        using namespace GravityModel::detail;
+        const auto &face = thrust::get<0>(tuple);
+        const auto &segmentVectors = thrust::get<1>(tuple);
+        const auto &planeUnitNormal = thrust::get<2>(tuple);
+        const auto &segmentUnitNormals = thrust::get<3>(tuple);
+        SPDLOG_LOGGER_TRACE(PolyhedralGravityLogger::DEFAULT_LOGGER.getLogger(),
+                            "Evaluating the plane with vertices: v1 = [{}, {}, {}], v2 = [{}, {}, {}], "
+                            "v3 = [{}, {}, {}]", 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]);
+        //1. Step: Compute ingredients for current plane which were not computed before
+        //1-04 Step: Compute Plane Normal Orientation sigma_p (direction of N_p in relation to P)
+        double planeNormalOrientation = computeUnitNormalOfPlaneDirection(planeUnitNormal, face[0]);
+        //1-05 Step: Compute Hessian Normal Plane Representation
+        HessianPlane hessianPlane = computeHessianPlane(face[0], face[1], face[2]);
+        //1-06 Step: Compute distance h_p between P and P'
+        double planeDistance = distanceBetweenOriginAndPlane(hessianPlane);
+        //1-07 Step: Compute the actual position of P' (projection of P on the plane)
+        Array3 orthogonalProjectionPointOnPlane = projectPointOrthogonallyOntoPlane(planeUnitNormal, planeDistance,
+                                                                                    hessianPlane);
+        //1-08 Step: Compute the segment normal orientation sigma_pq (direction of n_pq in relation to P')
+        Array3 segmentNormalOrientations = computeUnitNormalOfSegmentsDirections(face, orthogonalProjectionPointOnPlane,
+                                                                                 segmentUnitNormals);
+        //1-09 Step: Compute the orthogonal projection point P'' of P' on each segment
+        Array3Triplet orthogonalProjectionPointsOnSegmentsForPlane = projectPointOrthogonallyOntoSegments(
+                orthogonalProjectionPointOnPlane, segmentNormalOrientations, face);
+        //1-10 Step: Compute the segment distances h_pq between P'' and P'
+        Array3 segmentDistances = distancesBetweenProjectionPoints(orthogonalProjectionPointOnPlane,
+                                                                   orthogonalProjectionPointsOnSegmentsForPlane);
+        //1-11 Step: Compute the 3D distances l1, l2 (between P and vertices)
+        // and 1D distances s1, s2 (between P'' and vertices)
+        std::array<Distance, 3> distances = distancesToSegmentEndpoints(segmentVectors,
+                                                                        orthogonalProjectionPointsOnSegmentsForPlane,
+                                                                        face);
+        //1-12 Step: Compute the euclidian Norms of the vectors consisting of P and the vertices
+        // they are later used for determining the position of P in relation to the plane
+        Array3 projectionPointVertexNorms = computeNormsOfProjectionPointAndVertices(orthogonalProjectionPointOnPlane,
+                                                                                     face);
+        //1-13 Step: Compute the transcendental Expressions LN_pq and AN_pq
+        std::array<TranscendentalExpression, 3> transcendentalExpressions = computeTranscendentalExpressions(distances,
+                                                                                                             planeDistance,
+                                                                                                             segmentDistances,
+                                                                                                             segmentNormalOrientations,
+                                                                                                             projectionPointVertexNorms);
+        //1-14 Step: Compute the singularities sing A and sing B if P' is located in the plane,
+        // on any vertex, or on one segment (G_pq)
+        std::pair<double, Array3> singularities = computeSingularityTerms(segmentVectors, segmentNormalOrientations,
+                                                                          projectionPointVertexNorms, planeUnitNormal,
+                                                                          planeDistance, planeNormalOrientation);
+        //2. Step: Compute Sum 1 used for potential and acceleration (first derivative)
+        // sum over: sigma_pq * h_pq * LN_pq
+        // --> Equation 11/12 the first summation in the brackets
+        auto zipIteratorSum1PotentialAcceleration = util::zipPair(segmentNormalOrientations, segmentDistances,
+                                                                  transcendentalExpressions);
+        const double sum1PotentialAcceleration = std::accumulate(zipIteratorSum1PotentialAcceleration.first,
+                                                                 zipIteratorSum1PotentialAcceleration.second, 0.0,
+                                                                 [](double acc, const auto &tuple) {
+                                                                     const double &segmentOrientation = thrust::get<0>(
+                                                                             tuple);
+                                                                     const double &segmentDistance = thrust::get<1>(
+                                                                             tuple);
+                                                                     const TranscendentalExpression &transcendentalExpressions = thrust::get<2>(
+                                                                             tuple);
+                                                                     return acc + segmentOrientation * segmentDistance *
+                                                                                  transcendentalExpressions.ln;
+                                                                 });
+
+        //3. Step: Compute Sum 1 used for the gradiometric tensor (second derivative)
+        // sum over: n_pq * LN_pq
+        // --> Equation 13 the first summation in the brackets
+        auto zipIteratorSum1Tensor = util::zipPair(segmentUnitNormals, transcendentalExpressions);
+        const Array3 sum1Tensor = std::accumulate(zipIteratorSum1Tensor.first, zipIteratorSum1Tensor.second,
+                                                  Array3{0.0, 0.0, 0.0}, [](const Array3 &acc, const auto &tuple) {
+                    const Array3 &segmentNormal = thrust::get<0>(tuple);
+                    const TranscendentalExpression &transcendentalExpressions = thrust::get<1>(tuple);
+                    return acc + (segmentNormal * transcendentalExpressions.ln);
+                });
+
+        //4. Step: Compute Sum 2 which is the same for every result parameter
+        // sum over: sigma_pq * AN_pq
+        // --> Equation 11/12/13 the second summation in the brackets
+        auto zipIteratorSum2 = util::zipPair(segmentNormalOrientations, transcendentalExpressions);
+        const double sum2 = std::accumulate(zipIteratorSum2.first, zipIteratorSum2.second, 0.0,
+                                            [](double acc, const auto &tuple) {
+                                                const double &segmentOrientation = thrust::get<0>(tuple);
+                                                const TranscendentalExpression &transcendentalExpressions = thrust::get<1>(
+                                                        tuple);
+                                                return acc + segmentOrientation * transcendentalExpressions.an;
+                                            });
+
+        //5. Step: Sum for potential and acceleration
+        // consisting of: sum1 + h_p * sum2 + sing A
+        // --> Equation 11/12 the total sum of the brackets
+        const double planeSumPotentialAcceleration =
+                sum1PotentialAcceleration + planeDistance * sum2 + singularities.first;
+
+        if (isCriticalDifference(planeDistance, sum2)) {
+            // The multiplication planeDistance * sum2 is not the root cause, but both numbers are good
+            // 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");
+        }
+
+        //6. Step: Sum for tensor
+        // consisting of: sum1 + sigma_p * N_p * sum2 + sing B
+        // --> Equation 13 the total sum of the brackets
+        const Array3 subSum = (sum1Tensor + (planeUnitNormal * (planeNormalOrientation * sum2))) + singularities.second;
+        // first component: trivial case Vxx, Vyy, Vzz --> just N_p * subSum
+        // 00, 11, 22 --> xx, yy, zz with x as 0, y as 1, z as 2
+        const Array3 first = planeUnitNormal * subSum;
+        // second component: reordering required to build Vxy, Vxz, Vyz
+        // 01, 02, 12 --> xy, xz, yz with x as 0, y as 1, z as 2
+        const Array3 reorderedNp = {planeUnitNormal[0], planeUnitNormal[0], planeUnitNormal[1]};
+        const Array3 reorderedSubSum = {subSum[1], subSum[2], subSum[2]};
+        const Array3 second = reorderedNp * reorderedSubSum;
+
+        //7. Step: Multiply with prefix
+        // Equation (11): sigma_p * h_p * sum
+        // Equation (12): N_p * sum
+        // Equation (13): already done above, just concat the two components for later summation
+        return std::make_tuple(planeNormalOrientation * planeDistance * planeSumPotentialAcceleration,
+                               planeUnitNormal * planeSumPotentialAcceleration, concat(first, second));
+    }
+
+    std::string GravityEvaluable::toString() const {
+        std::stringstream ss;
+        ss  << "<polyhedral_gravity.GravityEvaluable, density=" << _density
+            << ", vertices= " << _polyhedron.countVertices()
+            << ", faces= " << _polyhedron.countFaces() << ">";
+        return ss.str();
+    }
+
+}
diff --git a/src/polyhedralGravity/calculation/GravityEvaluable.h b/src/polyhedralGravity/calculation/GravityEvaluable.h
new file mode 100644
index 0000000..9ebd704
--- /dev/null
+++ b/src/polyhedralGravity/calculation/GravityEvaluable.h
@@ -0,0 +1,145 @@
+#pragma once
+
+#include <tuple>
+#include <variant>
+#include <string>
+#include <sstream>
+
+#include "thrust/transform.h"
+#include "polyhedralGravity/calculation/GravityModel.h"
+
+#include "polyhedralGravity/input/TetgenAdapter.h"
+#include "polyhedralGravity/model/GravityModelData.h"
+#include "polyhedralGravity/model/Polyhedron.h"
+#include "thrust/execution_policy.h"
+
+
+namespace polyhedralGravity {
+
+    /**
+     * Class for evaluating the polyhedrale gravity model for a given constant density polyhedron.
+     * Caches the polyhedron and data which is independent of the computation point P.
+     * Provides an operator() for evaluating the polyhedrale gravity model for a given constant density polyhedron
+     * at computation point P and choosing between parallel and serial evaluation.
+     */
+    class GravityEvaluable {
+
+        /** The constant density polyhedron consisting of vertices and triangular faces */
+        const Polyhedron _polyhedron;
+
+        /** The constant density of the polyhedron in [kg/m^3] */
+        const double _density;
+
+        /** Cache for the segment vectors (segments between vertices of a polyhedral face) */
+        mutable std::vector<Array3Triplet> _segmentVectors;
+
+        /** Cache for the plane unit normals (unit normals of the polyhedral faces) */
+        mutable std::vector<Array3> _planeUnitNormals;
+
+        /** Cache for the segment unit normals (unit normals of each the polyhedral faces' segments) */
+        mutable std::vector<Array3Triplet> _segmentUnitNormals;
+
+
+    public:
+
+        /**
+         * Instantiates a GravityEvaluable with a given constant density polyhedron.
+         * @param polyhedron - the polyhedron
+         * @param density - the constant density in [kg/m^3]
+         *
+         * @note This is a separate constructor since the polyhedron as a class it not exposed to the user via
+         * the Python Interface. Thus, this constructor is only available via the C++ interface.
+         */
+        GravityEvaluable(const Polyhedron &polyhedron, double density) : _polyhedron{polyhedron}, _density{density} {
+            this->prepare();
+        }
+
+        /**
+         * Instantiates a GravityEvaluable with a given constant density polyhedron.
+         * @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{
+                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} {
+            this->prepare();
+        }
+
+        /**
+         * Evaluates the polyhedrale gravity model for a given constant density polyhedron at computation
+         * point P. Wrapper for evaluate<parallelization>.
+         * @param computationPoints - the computation point P or multiple computation points in a vector
+         * @param parallelization - if true, the calculation is parallelized
+         * @return the GravityModelResult containing the potential, acceleration and second derivative
+         */
+        inline std::variant<GravityModelResult, std::vector<GravityModelResult>>
+        operator()(const std::variant<Array3, std::vector<Array3>> &computationPoints,
+                   bool parallelization = true) const {
+            if (parallelization) {
+                if (std::holds_alternative<Array3>(computationPoints)) {
+                    return this->evaluate<true>(std::get<Array3>(computationPoints));
+                } else {
+                    return this->evaluate<true>(std::get<std::vector<Array3>>(computationPoints));
+                }
+            } else {
+                if (std::holds_alternative<Array3>(computationPoints)) {
+                    return this->evaluate<false>(std::get<Array3>(computationPoints));
+                } else {
+                    return this->evaluate<false>(std::get<std::vector<Array3>>(computationPoints));
+                }
+            }
+        }
+
+        /**
+         * Returns a string representation of the GravityEvaluable.
+         * @return string representation of the GravityEvaluable
+         */
+        std::string toString() const;
+
+    private:
+
+        /**
+         * Prepares the polyhedron for the evaluation by calculating the segment vectors, the plane unit normals
+         * and the segment unit normals.
+         * Called by the constructor once.
+         */
+        void prepare() const;
+
+        /**
+        * Evaluates the polyhedrale gravity model for a given constant density polyhedron at computation
+        * point P.
+        * @tparam Parallelization - if true, the calculation is parallelized
+        * @param computationPoint - the computation Point P
+        * @return the GravityModelResult containing the potential, the acceleration and the change of acceleration
+        * at computation Point P
+        */
+        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.
+         * @tparam Parallelization - if true, the calculation is parallelized
+         * @param computationPoints - the computation Points
+         * @return vector of GravityModelResults containing the potential, the acceleration and the change of acceleration
+         */
+        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
+         * @return the GravityModelResult containing the potential, the acceleration and the change of acceleration which
+         * this face contributes to the computation point
+         */
+        static GravityModelResult
+        evaluateFace(const thrust::tuple<Array3Triplet, Array3Triplet, Array3, Array3Triplet> &tuple);
+
+    };
+
+} // namespace polyhedralGravity
\ No newline at end of file
diff --git a/src/polyhedralGravity/calculation/GravityModel.cpp b/src/polyhedralGravity/calculation/GravityModel.cpp
index 7477dc1..a88fcdf 100644
--- a/src/polyhedralGravity/calculation/GravityModel.cpp
+++ b/src/polyhedralGravity/calculation/GravityModel.cpp
@@ -1,226 +1,40 @@
 #include "GravityModel.h"
 
-namespace polyhedralGravity {
-
-    GravityModelResult GravityModel::evaluate(
-            const Polyhedron &polyhedron, double density, const Array3 &computationPoint) {
-        using namespace util;
-        using namespace detail;
-        SPDLOG_LOGGER_DEBUG(PolyhedralGravityLogger::DEFAULT_LOGGER.getLogger(),
-                            "Evaluation for computation point P = [{}, {}, {}] started, given density = {} kg/m^3",
-                            computationPoint[0], computationPoint[1], computationPoint[2], density);
-        /*
-         * Calculate V and Vx, Vy, Vz and Vxx, Vyy, Vzz, Vxy, Vxz, Vyz
-         */
-        auto polyhedronIterator = transformPolyhedron(polyhedron, computationPoint);
-
-        SPDLOG_LOGGER_DEBUG(PolyhedralGravityLogger::DEFAULT_LOGGER.getLogger(),
-                            "Starting to iterate over the planes...");
-        GravityModelResult result{};
-        result = thrust::transform_reduce(
-                thrust::device,
-                polyhedronIterator.first,
-                polyhedronIterator.second,
-                [&computationPoint](const Array3Triplet &face) {
-                    using namespace util;
-                    using namespace detail;
-                    SPDLOG_LOGGER_TRACE(PolyhedralGravityLogger::DEFAULT_LOGGER.getLogger(),
-                                        "Evaluating the plane with vertices: v1 = [{}, {}, {}], v2 = [{}, {}, {}], "
-                                        "v3 = [{}, {}, {}]",
-                                        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]);
-                    //1. Step: Compute ingredients for current plane
-                    //1-01 Step: Compute Segment Vectors G_pq which describe each one the edge between two vertices
-                    Array3Triplet segmentVectors = buildVectorsOfSegments(face[0], face[1], face[2]);
-                    //1-02 Step: Compute the Plane Unit Normals N_p (pointing outside the polyhedron)
-                    Array3 planeUnitNormal = buildUnitNormalOfPlane(segmentVectors[0], segmentVectors[1]);
-                    //1-03 Step: Compute Segment Unit Normals n_pq (normal pointing away from each segment)
-                    Array3Triplet segmentUnitNormals = buildUnitNormalOfSegments(segmentVectors, planeUnitNormal);
-                    //1-04 Step: Compute Plane Normal Orientation sigma_p (direction of N_p in relation to P)
-                    double planeNormalOrientation = computeUnitNormalOfPlaneDirection(planeUnitNormal, face[0]);
-                    //1-05 Step: Compute Hessian Normal Plane Representation
-                    HessianPlane hessianPlane = computeHessianPlane(face[0], face[1], face[2]);
-                    //1-06 Step: Compute distance h_p between P and P'
-                    double planeDistance = distanceBetweenOriginAndPlane(hessianPlane);
-                    //1-07 Step: Compute the actual position of P' (projection of P on the plane)
-                    Array3 orthogonalProjectionPointOnPlane = projectPointOrthogonallyOntoPlane(
-                            planeUnitNormal,
-                            planeDistance,
-                            hessianPlane);
-                    //1-08 Step: Compute the segment normal orientation sigma_pq (direction of n_pq in relation to P')
-                    Array3 segmentNormalOrientations = computeUnitNormalOfSegmentsDirections(
-                            face,
-                            orthogonalProjectionPointOnPlane,
-                            segmentUnitNormals);
-                    //1-09 Step: Compute the orthogonal projection point P'' of P' on each segment
-                    Array3Triplet orthogonalProjectionPointsOnSegmentsForPlane = projectPointOrthogonallyOntoSegments(
-                            orthogonalProjectionPointOnPlane,
-                            segmentNormalOrientations,
-                            face);
-                    //1-10 Step: Compute the segment distances h_pq between P'' and P'
-                    Array3 segmentDistances = distancesBetweenProjectionPoints(
-                            orthogonalProjectionPointOnPlane,
-                            orthogonalProjectionPointsOnSegmentsForPlane);
-                    //1-11 Step: Compute the 3D distances l1, l2 (between P and vertices)
-                    // and 1D distances s1, s2 (between P'' and vertices)
-                    std::array<Distance, 3> distances = distancesToSegmentEndpoints(
-                            segmentVectors,
-                            orthogonalProjectionPointsOnSegmentsForPlane,
-                            face);
-                    //1-12 Step: Compute the euclidian Norms of the vectors consisting of P and the vertices
-                    // they are later used for determining the position of P in relation to the plane
-                    Array3 projectionPointVertexNorms = computeNormsOfProjectionPointAndVertices(
-                            orthogonalProjectionPointOnPlane,
-                            face);
-                    //1-13 Step: Compute the transcendental Expressions LN_pq and AN_pq
-                    std::array<TranscendentalExpression, 3> transcendentalExpressions = computeTranscendentalExpressions(
-                            distances,
-                            planeDistance,
-                            segmentDistances,
-                            segmentNormalOrientations,
-                            projectionPointVertexNorms);
-                    //1-14 Step: Compute the singularities sing A and sing B if P' is located in the plane,
-                    // on any vertex, or on one segment (G_pq)
-                    std::pair<double, Array3> singularities = computeSingularityTerms(
-                            segmentVectors,
-                            segmentNormalOrientations,
-                            projectionPointVertexNorms,
-                            planeUnitNormal, planeDistance,
-                            planeNormalOrientation);
-                    //2. Step: Compute Sum 1 used for potential and acceleration (first derivative)
-                    // sum over: sigma_pq * h_pq * LN_pq
-                    // --> Equation 11/12 the first summation in the brackets
-                    auto zipIteratorSum1PotentialAcceleration =
-                            util::zipPair(segmentNormalOrientations, segmentDistances, transcendentalExpressions);
-                    const double sum1PotentialAcceleration = std::accumulate(
-                            zipIteratorSum1PotentialAcceleration.first,
-                            zipIteratorSum1PotentialAcceleration.second,
-                            0.0,
-                            [](double acc, const auto &tuple) {
-                                const double &segmentOrientation = thrust::get<0>(tuple);
-                                const double &segmentDistance = thrust::get<1>(tuple);
-                                const TranscendentalExpression &transcendentalExpressions = thrust::get<2>(tuple);
-                                return acc + segmentOrientation * segmentDistance * transcendentalExpressions.ln;
-                            });
-
-                    //3. Step: Compute Sum 1 used for the gradiometric tensor (second derivative)
-                    // sum over: n_pq * LN_pq
-                    // --> Equation 13 the first summation in the brackets
-                    auto zipIteratorSum1Tensor = util::zipPair(segmentUnitNormals, transcendentalExpressions);
-                    const Array3 sum1Tensor = std::accumulate(
-                            zipIteratorSum1Tensor.first,
-                            zipIteratorSum1Tensor.second,
-                            Array3{0.0, 0.0, 0.0},
-                            [](const Array3 &acc, const auto &tuple) {
-                                const Array3 &segmentNormal = thrust::get<0>(tuple);
-                                const TranscendentalExpression &transcendentalExpressions = thrust::get<1>(tuple);
-                                return acc + (segmentNormal * transcendentalExpressions.ln);
-                            });
-
-                    //4. Step: Compute Sum 2 which is the same for every result parameter
-                    // sum over: sigma_pq * AN_pq
-                    // --> Equation 11/12/13 the second summation in the brackets
-                    auto zipIteratorSum2 = util::zipPair(segmentNormalOrientations, transcendentalExpressions);
-                    const double sum2 = std::accumulate(
-                            zipIteratorSum2.first,
-                            zipIteratorSum2.second,
-                            0.0,
-                            [](double acc, const auto &tuple) {
-                                const double &segmentOrientation = thrust::get<0>(tuple);
-                                const TranscendentalExpression &transcendentalExpressions = thrust::get<1>(tuple);
-                                return acc + segmentOrientation * transcendentalExpressions.an;
-                            });
-
-                    //5. Step: Sum for potential and acceleration
-                    // consisting of: sum1 + h_p * sum2 + sing A
-                    // --> Equation 11/12 the total sum of the brackets
-                    const double planeSumPotentialAcceleration =
-                            sum1PotentialAcceleration + planeDistance * sum2 + singularities.first;
-
-                    if (isCriticalDifference(planeDistance, sum2)) {
-                        // The multiplication planeDistance * sum2 is not the root cause, but both numbers are good
-                        // 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",
-                                computationPoint[0], computationPoint[1], computationPoint[2]);
-                    }
-
-                    //6. Step: Sum for tensor
-                    // consisting of: sum1 + sigma_p * N_p * sum2 + sing B
-                    // --> Equation 13 the total sum of the brackets
-                    const Array3 subSum = (sum1Tensor + (planeUnitNormal * (planeNormalOrientation * sum2))) + singularities.second;
-                    // first component: trivial case Vxx, Vyy, Vzz --> just N_p * subSum
-                    // 00, 11, 22 --> xx, yy, zz with x as 0, y as 1, z as 2
-                    const Array3 first = planeUnitNormal * subSum;
-                    // second component: reordering required to build Vxy, Vxz, Vyz
-                    // 01, 02, 12 --> xy, xz, yz with x as 0, y as 1, z as 2
-                    const Array3 reorderedNp = {planeUnitNormal[0], planeUnitNormal[0], planeUnitNormal[1]};
-                    const Array3 reorderedSubSum = {subSum[1], subSum[2], subSum[2]};
-                    const Array3 second = reorderedNp * reorderedSubSum;
-
-                    //7. Step: Multiply with prefix
-                    // Equation (11): sigma_p * h_p * sum
-                    // Equation (12): N_p * sum
-                    // Equation (13): already done above, just concat the two components for later summation
-                    return GravityModelResult{
-                            planeNormalOrientation * planeDistance * planeSumPotentialAcceleration,
-                            planeUnitNormal * planeSumPotentialAcceleration,
-                            concat(first, second)
-                    };
-                },
-                result,
-                [](const GravityModelResult &a, const GravityModelResult &b) {
-                    //8. Step: Accumulate the partial sums
-                    return GravityModelResult{
-                            a.gravitationalPotential + b.gravitationalPotential,
-                            a.acceleration + b.acceleration,
-                            a.gradiometricTensor + b.gradiometricTensor
-                    };
-                });
-
-        SPDLOG_LOGGER_DEBUG(PolyhedralGravityLogger::DEFAULT_LOGGER.getLogger(),
-                            "Finished the sums. Applying final prefix and eliminating rounding errors.");
-
-        //9. Step: Eliminate rounding errors in the result and set those tiny values to "really" zero
-        result.eliminateRoundingErrors();
-
-        //10. Step: Compute prefix consisting of GRAVITATIONAL_CONSTANT * density
-        const double prefix = util::GRAVITATIONAL_CONSTANT * density;
-
-        //11. Step: Final expressions after application of the prefix (and a division by 2 for the potential)
-        result.gravitationalPotential = (result.gravitationalPotential * prefix) / 2.0;
-        result.acceleration = result.acceleration * prefix;
-        result.gradiometricTensor = result.gradiometricTensor * prefix;
-        return result;
+namespace polyhedralGravity::GravityModel {
+
+    GravityModelResult evaluate(
+            const PolyhedronOrSource &polyhedron,
+            double density, const Array3 &computationPoint, bool parallel) {
+        return std::get<GravityModelResult>(
+                std::visit([&parallel, &density, &computationPoint](const auto &polyhedron) {
+                    GravityEvaluable evaluable{polyhedron, density};
+                    return evaluable(computationPoint, parallel);
+                }, polyhedron));
     }
 
-    std::vector<GravityModelResult> GravityModel::evaluate(
-            const Polyhedron &polyhedron, double density, const std::vector<Array3> &computationPoints) {
-        std::vector<GravityModelResult> result{computationPoints.size()};
-        thrust::transform(computationPoints.begin(), computationPoints.end(), result.begin(),
-                          [&polyhedron, density](const Array3 &computationPoint) {
-                              return evaluate(polyhedron, density, computationPoint);
-                          });
-        return result;
+    std::vector<GravityModelResult>
+    evaluate(const PolyhedronOrSource &polyhedron, double density, const std::vector<Array3> &computationPoints,
+             bool parallel) {
+        return std::get<std::vector<GravityModelResult>>(std::visit([&parallel, &density, &computationPoints](const auto &polyhedron) {
+            GravityEvaluable evaluable{polyhedron, density};
+            return evaluable(computationPoints, parallel);
+        }, polyhedron));
     }
 
-    Array3Triplet GravityModel::detail::buildVectorsOfSegments(
-            const Array3 &vertex0, const Array3 &vertex1, const Array3 &vertex2) {
+
+    Array3Triplet detail::buildVectorsOfSegments(const Array3 &vertex0, const Array3 &vertex1, const Array3 &vertex2) {
         using util::operator-;
         //Calculate G_ij
         return {vertex1 - vertex0, vertex2 - vertex1, vertex0 - vertex2};
     }
 
-    Array3 GravityModel::detail::buildUnitNormalOfPlane(const Array3 &segmentVector1, const Array3 &segmentVector2) {
+    Array3 detail::buildUnitNormalOfPlane(const Array3 &segmentVector1, const Array3 &segmentVector2) {
         //Calculate N_i as (G_i1 * G_i2) / |G_i1 * G_i2| with * being the cross product
         return util::normal(segmentVector1, segmentVector2);
     }
 
-    Array3Triplet GravityModel::detail::buildUnitNormalOfSegments(
-            const Array3Triplet &segmentVectors, const Array3 &planeUnitNormal) {
+    Array3Triplet
+    detail::buildUnitNormalOfSegments(const Array3Triplet &segmentVectors, const Array3 &planeUnitNormal) {
         Array3Triplet segmentUnitNormal{};
         //Calculate n_ij as (G_ij * N_i) / |G_ig * N_i| with * being the cross product
         std::transform(segmentVectors.cbegin(), segmentVectors.end(), segmentUnitNormal.begin(),
@@ -230,8 +44,7 @@ namespace polyhedralGravity {
         return segmentUnitNormal;
     }
 
-    double
-    GravityModel::detail::computeUnitNormalOfPlaneDirection(const Array3 &planeUnitNormal, const Array3 &vertex0) {
+    double detail::computeUnitNormalOfPlaneDirection(const Array3 &planeUnitNormal, const Array3 &vertex0) {
         using namespace util;
         //Calculate N_i * -G_i1 where * is the dot product and then use the inverted sgn
         //We abstain on the double multiplication with -1 in the line above and beyond since two
@@ -239,7 +52,7 @@ namespace polyhedralGravity {
         return sgn(dot(planeUnitNormal, vertex0), util::EPSILON);
     }
 
-    HessianPlane GravityModel::detail::computeHessianPlane(const Array3 &p, const Array3 &q, const Array3 &r) {
+    HessianPlane detail::computeHessianPlane(const Array3 &p, const Array3 &q, const Array3 &r) {
         using namespace util;
         constexpr Array3 origin{0.0, 0.0, 0.0};
         const auto crossProduct = cross(p - q, p - r);
@@ -249,18 +62,14 @@ namespace polyhedralGravity {
         return {crossProduct[0], crossProduct[1], crossProduct[2], d};
     }
 
-    double GravityModel::detail::distanceBetweenOriginAndPlane(const HessianPlane &hessianPlane) {
+    double detail::distanceBetweenOriginAndPlane(const HessianPlane &hessianPlane) {
         //Compute h_p as D/sqrt(A^2 + B^2 + C^2)
         return std::abs(hessianPlane.d / std::sqrt(
-                hessianPlane.a * hessianPlane.a +
-                hessianPlane.b * hessianPlane.b +
-                hessianPlane.c * hessianPlane.c));
+                hessianPlane.a * hessianPlane.a + hessianPlane.b * hessianPlane.b + hessianPlane.c * hessianPlane.c));
     }
 
-    Array3 GravityModel::detail::projectPointOrthogonallyOntoPlane(
-            const Array3 &planeUnitNormal,
-            double planeDistance,
-            const HessianPlane &hessianPlane) {
+    Array3 detail::projectPointOrthogonallyOntoPlane(const Array3 &planeUnitNormal, double planeDistance,
+                                                     const HessianPlane &hessianPlane) {
         using namespace util;
         //Calculate the projection point by (22) P'_ = N_i / norm(N_i) * h_i
         // norm(N_i) is always 1 since N_i is a "normed" vector --> we do not need this division
@@ -294,10 +103,9 @@ namespace polyhedralGravity {
         return orthogonalProjectionPoint;
     }
 
-    Array3 GravityModel::detail::computeUnitNormalOfSegmentsDirections(
-            const Array3Triplet &vertices,
-            const Array3 &projectionPointOnPlane,
-            const Array3Triplet &segmentUnitNormalsForPlane) {
+    Array3
+    detail::computeUnitNormalOfSegmentsDirections(const Array3Triplet &vertices, const Array3 &projectionPointOnPlane,
+                                                  const Array3Triplet &segmentUnitNormalsForPlane) {
         using namespace util;
         std::array<double, 3> segmentNormalOrientations{};
         //Equation (23)
@@ -305,26 +113,24 @@ namespace polyhedralGravity {
         //i.e. the coordinates of the training-planes' nodes --> projectionPointOnPlane - vertex
         //Calculate n_ij * x_ij with * being the dot product and use the inverted sgn to determine the value of sigma_pq
         //--> sgn((dot(segmentNormalOrientation, projectionPointOnPlane - vertex)), util::EPSILON) * -1.0
-        std::transform(segmentUnitNormalsForPlane.cbegin(), segmentUnitNormalsForPlane.cend(),
-                       vertices.cbegin(), segmentNormalOrientations.begin(),
+        std::transform(segmentUnitNormalsForPlane.cbegin(), segmentUnitNormalsForPlane.cend(), vertices.cbegin(),
+                       segmentNormalOrientations.begin(),
                        [&projectionPointOnPlane](const Array3 &segmentUnitNormal, const Array3 &vertex) {
                            using namespace util;
-                           return sgn((dot(segmentUnitNormal, projectionPointOnPlane - vertex)), util::EPSILON) *
-                                  -1.0;
+                           return sgn((dot(segmentUnitNormal, projectionPointOnPlane - vertex)), util::EPSILON) * -1.0;
                        });
         return segmentNormalOrientations;
     }
 
-    Array3Triplet GravityModel::detail::projectPointOrthogonallyOntoSegments(
-            const Array3 &projectionPointOnPlane,
-            const Array3 &segmentNormalOrientations,
-            const Array3Triplet &face) {
+    Array3Triplet detail::projectPointOrthogonallyOntoSegments(const Array3 &projectionPointOnPlane,
+                                                               const Array3 &segmentNormalOrientations,
+                                                               const Array3Triplet &face) {
         auto counterJ = thrust::counting_iterator<unsigned int>(0);
         Array3Triplet orthogonalProjectionPointOnSegmentPerPlane{};
 
         //Running over the segments of this plane
-        thrust::transform(segmentNormalOrientations.begin(), segmentNormalOrientations.end(),
-                          counterJ, orthogonalProjectionPointOnSegmentPerPlane.begin(),
+        thrust::transform(segmentNormalOrientations.begin(), segmentNormalOrientations.end(), counterJ,
+                          orthogonalProjectionPointOnSegmentPerPlane.begin(),
                           [&](const double &segmentNormalOrientation, const unsigned int j) {
                               //We actually only accept +0.0 or -0.0, so the equal comparison is ok
                               if (segmentNormalOrientation == 0.0) {
@@ -335,35 +141,31 @@ namespace polyhedralGravity {
                                   //with the endpoints v1 and v2
                                   const auto &vertex1 = face[j];
                                   const auto &vertex2 = face[(j + 1) % 3];
-                                  return projectPointOrthogonallyOntoSegment(vertex1, vertex2,
-                                                                             projectionPointOnPlane);
+                                  return projectPointOrthogonallyOntoSegment(vertex1, vertex2, projectionPointOnPlane);
                               }
                           });
         return orthogonalProjectionPointOnSegmentPerPlane;
     }
 
-    Array3 GravityModel::detail::projectPointOrthogonallyOntoSegment(const Array3 &vertex1, const Array3 &vertex2,
-                                                                     const Array3 &orthogonalProjectionPointOnPlane) {
+    Array3 detail::projectPointOrthogonallyOntoSegment(const Array3 &vertex1, const Array3 &vertex2,
+                                                       const Array3 &orthogonalProjectionPointOnPlane) {
         using namespace util;
         //Preparing our the planes/ equations in matrix form
         const Array3 matrixRow1 = vertex2 - vertex1;
         const Array3 matrixRow2 = cross(vertex1 - orthogonalProjectionPointOnPlane, matrixRow1);
         const Array3 matrixRow3 = cross(matrixRow2, matrixRow1);
         const Array3 d = {dot(matrixRow1, orthogonalProjectionPointOnPlane),
-                          dot(matrixRow2, orthogonalProjectionPointOnPlane),
-                          dot(matrixRow3, vertex1)};
+                          dot(matrixRow2, orthogonalProjectionPointOnPlane), dot(matrixRow3, vertex1)};
         Matrix<double, 3, 3> columnMatrix = transpose(Matrix<double, 3, 3>{matrixRow1, matrixRow2, matrixRow3});
         //Calculation and solving the equations of above
         const double determinant = det(columnMatrix);
-        return Array3{
-                det(Matrix<double, 3, 3>{d, columnMatrix[1], columnMatrix[2]}),
-                det(Matrix<double, 3, 3>{columnMatrix[0], d, columnMatrix[2]}),
-                det(Matrix<double, 3, 3>{columnMatrix[0], columnMatrix[1], d})
-        } / determinant;
+        return Array3{det(Matrix<double, 3, 3>{d, columnMatrix[1], columnMatrix[2]}),
+                      det(Matrix<double, 3, 3>{columnMatrix[0], d, columnMatrix[2]}),
+                      det(Matrix<double, 3, 3>{columnMatrix[0], columnMatrix[1], d})} / determinant;
     }
 
-    Array3 GravityModel::detail::distancesBetweenProjectionPoints(const Array3 &orthogonalProjectionPointOnPlane,
-                                                                  const Array3Triplet &orthogonalProjectionPointOnSegments) {
+    Array3 detail::distancesBetweenProjectionPoints(const Array3 &orthogonalProjectionPointOnPlane,
+                                                    const Array3Triplet &orthogonalProjectionPointOnSegments) {
         std::array<double, 3> segmentDistances{};
         //The inner loop with the running j --> iterating over the segments
         //Using the values P'_i and P''_ij for the calculation of the distance
@@ -375,164 +177,160 @@ namespace polyhedralGravity {
         return segmentDistances;
     }
 
-    std::array<Distance, 3> GravityModel::detail::distancesToSegmentEndpoints(
-            const Array3Triplet &segmentVectorsForPlane,
-            const Array3Triplet &orthogonalProjectionPointsOnSegmentForPlane,
-            const Array3Triplet &face) {
+    std::array<Distance, 3> detail::distancesToSegmentEndpoints(const Array3Triplet &segmentVectorsForPlane,
+                                                                const Array3Triplet &orthogonalProjectionPointsOnSegmentForPlane,
+                                                                const Array3Triplet &face) {
         std::array<Distance, 3> distancesForPlane{};
         auto counter = thrust::counting_iterator<unsigned int>(0);
         auto zip = util::zipPair(segmentVectorsForPlane, orthogonalProjectionPointsOnSegmentForPlane);
 
-        thrust::transform(zip.first, zip.second, counter,
-                          distancesForPlane.begin(), [&face](const auto &tuple, unsigned int j) {
-                    using namespace util;
-                    Distance distance{};
-                    //segment vector G_pq
-                    const Array3 &segmentVector = thrust::get<0>(tuple);
-                    //orthogonal projection point on segment P'' for plane p and segment q
-                    const Array3 &orthogonalProjectionPointsOnSegment = thrust::get<1>(tuple);
-
-                    //Calculate the 3D distances between P (0, 0, 0) and
-                    // the segment endpoints face[j] and face[(j + 1) % 3])
-                    distance.l1 = euclideanNorm(face[j]);
-                    distance.l2 = euclideanNorm(face[(j + 1) % 3]);
-                    //Calculate the 1D distances between P'' (every segment has its own) and
-                    // the segment endpoints face[j] and face[(j + 1) % 3])
-                    distance.s1 = euclideanNorm(orthogonalProjectionPointsOnSegment - face[j]);
-                    distance.s2 = euclideanNorm(orthogonalProjectionPointsOnSegment - face[(j + 1) % 3]);
-
-                    /*
-                     * Additional remark:
-                     * Details on these conditions are in the second paper referenced in the README.md (Tsoulis, 2021)
-                     * The numbering of these conditions is equal to the numbering scheme of the paper
-                     * Assign a sign to those magnitudes depending on the relative position of P'' to the two
-                     * segment endpoints
-                     */
-
-                    //4. Option: |s1 - l1| == 0 && |s2 - l2| == 0 Computation point P is located from the beginning on
-                    // the direction of a specific segment (P coincides with P' and P'')
-                    if (std::abs(distance.s1 - distance.l1) < EPSILON &&
-                        std::abs(distance.s2 - distance.l2) < EPSILON) {
-                        //4. Option - Case 2: P is located on the segment from its right side
-                        // s1 = -|s1|, s2 = -|s2|, l1 = -|l1|, l2 = -|l2|
-                        if (distance.s2 < distance.s1) {
-                            distance.s1 *= -1.0;
-                            distance.s2 *= -1.0;
-                            distance.l1 *= -1.0;
-                            distance.l2 *= -1.0;
-                            return distance;
-                        } else if (std::abs(distance.s2 - distance.s1) < EPSILON) {
-                            //4. Option - Case 1: P is located inside the segment (s2 == s1)
-                            // s1 = -|s1|, s2 = |s2|, l1 = -|l1|, l2 = |l2|
-                            distance.s1 *= -1.0;
-                            distance.l1 *= -1.0;
-                            return distance;
-                        }
-                        //4. Option - Case 3: P is located on the segment from its left side
-                        // s1 = |s1|, s2 = |s2|, l1 = |l1|, l2 = |l2| --> Nothing to do!
-                    } else {
-                        const double norm = euclideanNorm(segmentVector);
-                        if (distance.s1 < norm && distance.s2 < norm) {
-                            //1. Option: |s1| < |G_ij| && |s2| < |G_ij| Point P'' is situated inside the segment
-                            // s1 = -|s1|, s2 = |s2|, l1 = |l1|, l2 = |l2|
-                            distance.s1 *= -1.0;
-                            return distance;
-                        } else if (distance.s2 < distance.s1) {
-                            //2. Option: |s2| < |s1| Point P'' is on the right side of the segment
-                            // s1 = -|s1|, s2 = -|s2|, l1 = |l1|, l2 = |l2|
-                            distance.s1 *= -1.0;
-                            distance.s2 *= -1.0;
-                            return distance;
-                        }
-                        //3. Option: |s1| < |s2| Point P'' is on the left side of the segment
-                        // s1 = |s1|, s2 = |s2|, l1 = |l1|, l2 = |l2| --> Nothing to do!
-                    }
-                    return distance;
-                });
+        thrust::transform(zip.first, zip.second, counter, distancesForPlane.begin(),
+                          [&face](const auto &tuple, unsigned int j) {
+                              using namespace util;
+                              Distance distance{};
+                              //segment vector G_pq
+                              const Array3 &segmentVector = thrust::get<0>(tuple);
+                              //orthogonal projection point on segment P'' for plane p and segment q
+                              const Array3 &orthogonalProjectionPointsOnSegment = thrust::get<1>(tuple);
+
+                              //Calculate the 3D distances between P (0, 0, 0) and
+                              // the segment endpoints face[j] and face[(j + 1) % 3])
+                              distance.l1 = euclideanNorm(face[j]);
+                              distance.l2 = euclideanNorm(face[(j + 1) % 3]);
+                              //Calculate the 1D distances between P'' (every segment has its own) and
+                              // the segment endpoints face[j] and face[(j + 1) % 3])
+                              distance.s1 = euclideanNorm(orthogonalProjectionPointsOnSegment - face[j]);
+                              distance.s2 = euclideanNorm(orthogonalProjectionPointsOnSegment - face[(j + 1) % 3]);
+
+                              /*
+                               * Additional remark:
+                               * Details on these conditions are in the second paper referenced in the README.md (Tsoulis, 2021)
+                               * The numbering of these conditions is equal to the numbering scheme of the paper
+                               * Assign a sign to those magnitudes depending on the relative position of P'' to the two
+                               * segment endpoints
+                               */
+
+                              //4. Option: |s1 - l1| == 0 && |s2 - l2| == 0 Computation point P is located from the beginning on
+                              // the direction of a specific segment (P coincides with P' and P'')
+                              if (std::abs(distance.s1 - distance.l1) < EPSILON &&
+                                  std::abs(distance.s2 - distance.l2) < EPSILON) {
+                                  //4. Option - Case 2: P is located on the segment from its right side
+                                  // s1 = -|s1|, s2 = -|s2|, l1 = -|l1|, l2 = -|l2|
+                                  if (distance.s2 < distance.s1) {
+                                      distance.s1 *= -1.0;
+                                      distance.s2 *= -1.0;
+                                      distance.l1 *= -1.0;
+                                      distance.l2 *= -1.0;
+                                      return distance;
+                                  } else if (std::abs(distance.s2 - distance.s1) < EPSILON) {
+                                      //4. Option - Case 1: P is located inside the segment (s2 == s1)
+                                      // s1 = -|s1|, s2 = |s2|, l1 = -|l1|, l2 = |l2|
+                                      distance.s1 *= -1.0;
+                                      distance.l1 *= -1.0;
+                                      return distance;
+                                  }
+                                  //4. Option - Case 3: P is located on the segment from its left side
+                                  // s1 = |s1|, s2 = |s2|, l1 = |l1|, l2 = |l2| --> Nothing to do!
+                              } else {
+                                  const double norm = euclideanNorm(segmentVector);
+                                  if (distance.s1 < norm && distance.s2 < norm) {
+                                      //1. Option: |s1| < |G_ij| && |s2| < |G_ij| Point P'' is situated inside the segment
+                                      // s1 = -|s1|, s2 = |s2|, l1 = |l1|, l2 = |l2|
+                                      distance.s1 *= -1.0;
+                                      return distance;
+                                  } else if (distance.s2 < distance.s1) {
+                                      //2. Option: |s2| < |s1| Point P'' is on the right side of the segment
+                                      // s1 = -|s1|, s2 = -|s2|, l1 = |l1|, l2 = |l2|
+                                      distance.s1 *= -1.0;
+                                      distance.s2 *= -1.0;
+                                      return distance;
+                                  }
+                                  //3. Option: |s1| < |s2| Point P'' is on the left side of the segment
+                                  // s1 = |s1|, s2 = |s2|, l1 = |l1|, l2 = |l2| --> Nothing to do!
+                              }
+                              return distance;
+                          });
         return distancesForPlane;
     }
 
     std::array<TranscendentalExpression, 3>
-    GravityModel::detail::computeTranscendentalExpressions(
-            const std::array<Distance, 3> &distancesForPlane,
-            double planeDistance, const Array3 &segmentDistancesForPlane,
-            const Array3 &segmentNormalOrientationsForPlane,
-            const Array3 &projectionPointVertexNorms) {
+    detail::computeTranscendentalExpressions(const std::array<Distance, 3> &distancesForPlane, double planeDistance,
+                                             const Array3 &segmentDistancesForPlane,
+                                             const Array3 &segmentNormalOrientationsForPlane,
+                                             const Array3 &projectionPointVertexNorms) {
         std::array<TranscendentalExpression, 3> transcendentalExpressionsForPlane{};
 
         //Zip iterator consisting of 3D and 1D distances l1/l2 and s1/2 for this plane | h_pq | sigma_pq for this plane
         auto zip = util::zipPair(distancesForPlane, segmentDistancesForPlane, segmentNormalOrientationsForPlane);
         auto counter = thrust::counting_iterator<unsigned int>(0);
 
-        thrust::transform(zip.first, zip.second, counter,
-                          transcendentalExpressionsForPlane.begin(), [&](const auto &tuple, const unsigned int j) {
-                    using namespace util;
-                    //distances l1, l2, s1, s1 for this segment q of plane p
-                    const Distance &distance = thrust::get<0>(tuple);
-                    //segment distance h_pq for this segment q of plane p
-                    const double segmentDistance = thrust::get<1>(tuple);
-                    //segment normal orientation sigma_pq for this segment q of plane p
-                    const double segmentNormalOrientation = thrust::get<2>(tuple);
-
-                    //Result for this segment
-                    TranscendentalExpression transcendentalExpressionPerSegment{};
-
-                    //Computation of the norm of P' and segment endpoints
-                    // If the one of the norms == 0 then P' lies on the corresponding vertex and coincides with P''
-                    const double r1Norm = projectionPointVertexNorms[(j + 1) % 3];
-                    const double r2Norm = projectionPointVertexNorms[j];
-
-                    //Compute LN_pq according to (14)
-                    // If sigma_pq == 0 && either of the distances of P' to the two segment endpoints == 0 OR
-                    // the 1D and 3D distances are smaller than some EPSILON
-                    // then LN_pq can be set to zero
-                    if ((segmentNormalOrientation == 0.0 && (r1Norm < EPSILON || r2Norm < EPSILON)) ||
-                        (std::abs(distance.s1 + distance.s2) < EPSILON &&
-                         std::abs(distance.l1 + distance.l2) < EPSILON)) {
-                        transcendentalExpressionPerSegment.ln = 0.0;
-                    } else {
-                        //Implementation of
-                        // log((s2_pq + l2_pq) / (s1_pq + l1_pq))
-                        transcendentalExpressionPerSegment.ln =
-                                std::log((distance.s2 + distance.l2) / (distance.s1 + distance.l1));
-                    }
-
-                    //Compute AN_pq according to (15)
-                    // If h_p == 0 or h_pq == 0 then AN_pq is zero, too (distances are always positive!)
-                    if (planeDistance < EPSILON || segmentDistance < EPSILON) {
-                        transcendentalExpressionPerSegment.an = 0.0;
-                    } else {
-                        //Implementation of:
-                        // atan(h_p * s2_pq / h_pq * l2_pq) - atan(h_p * s1_pq / h_pq * l1_pq)
-                        // in a vectorized manner to reduce the number of atan(..) and divisions(..)
-
-                        //The last subtraction is implemented via -distance.l1 (the minus/ inversion
-                        // is sustained through all operations, even the atan(..)
-                        xsimd::batch<double> reg1({distance.s2, distance.s1});
-                        xsimd::batch<double> reg2({distance.l2, -distance.l1});
-
-                        reg1 = xsimd::mul(reg1, planeDistance);
-                        reg2 = xsimd::mul(reg2, segmentDistance);
-
-                        reg1 = xsimd::div(reg1, reg2);
-                        reg1 = xsimd::atan(reg1);
-
-                        // "Subtraction" masked as addition (l1 was previously inverted with a minus)
-                        transcendentalExpressionPerSegment.an = xsimd::hadd(reg1);
-                    }
-
-                    return transcendentalExpressionPerSegment;
-                });
+        thrust::transform(zip.first, zip.second, counter, transcendentalExpressionsForPlane.begin(),
+                          [&](const auto &tuple, const unsigned int j) {
+                              using namespace util;
+                              //distances l1, l2, s1, s1 for this segment q of plane p
+                              const Distance &distance = thrust::get<0>(tuple);
+                              //segment distance h_pq for this segment q of plane p
+                              const double segmentDistance = thrust::get<1>(tuple);
+                              //segment normal orientation sigma_pq for this segment q of plane p
+                              const double segmentNormalOrientation = thrust::get<2>(tuple);
+
+                              //Result for this segment
+                              TranscendentalExpression transcendentalExpressionPerSegment{};
+
+                              //Computation of the norm of P' and segment endpoints
+                              // If the one of the norms == 0 then P' lies on the corresponding vertex and coincides with P''
+                              const double r1Norm = projectionPointVertexNorms[(j + 1) % 3];
+                              const double r2Norm = projectionPointVertexNorms[j];
+
+                              //Compute LN_pq according to (14)
+                              // If sigma_pq == 0 && either of the distances of P' to the two segment endpoints == 0 OR
+                              // the 1D and 3D distances are smaller than some EPSILON
+                              // then LN_pq can be set to zero
+                              if ((segmentNormalOrientation == 0.0 && (r1Norm < EPSILON || r2Norm < EPSILON)) ||
+                                  (std::abs(distance.s1 + distance.s2) < EPSILON &&
+                                   std::abs(distance.l1 + distance.l2) < EPSILON)) {
+                                  transcendentalExpressionPerSegment.ln = 0.0;
+                              } else {
+                                  //Implementation of
+                                  // log((s2_pq + l2_pq) / (s1_pq + l1_pq))
+                                  transcendentalExpressionPerSegment.ln = std::log(
+                                          (distance.s2 + distance.l2) / (distance.s1 + distance.l1));
+                              }
+
+                              //Compute AN_pq according to (15)
+                              // If h_p == 0 or h_pq == 0 then AN_pq is zero, too (distances are always positive!)
+                              if (planeDistance < EPSILON || segmentDistance < EPSILON) {
+                                  transcendentalExpressionPerSegment.an = 0.0;
+                              } else {
+                                  //Implementation of:
+                                  // atan(h_p * s2_pq / h_pq * l2_pq) - atan(h_p * s1_pq / h_pq * l1_pq)
+                                  // in a vectorized manner to reduce the number of atan(..) and divisions(..)
+
+                                  //The last subtraction is implemented via -distance.l1 (the minus/ inversion
+                                  // is sustained through all operations, even the atan(..)
+                                  xsimd::batch<double> reg1(distance.s2, distance.s1);
+                                  xsimd::batch<double> reg2(distance.l2, -distance.l1);
+
+                                  reg1 = xsimd::mul(reg1, planeDistance);
+                                  reg2 = xsimd::mul(reg2, segmentDistance);
+
+                                  reg1 = xsimd::div(reg1, reg2);
+                                  reg1 = xsimd::atan(reg1);
+
+                                  // "Subtraction" masked as addition (l1 was previously inverted with a minus)
+                                  transcendentalExpressionPerSegment.an = xsimd::reduce_add(reg1);
+                              }
+
+                              return transcendentalExpressionPerSegment;
+                          });
         return transcendentalExpressionsForPlane;
     }
 
-    std::pair<double, Array3> GravityModel::detail::computeSingularityTerms(
-            const Array3Triplet &segmentVectorsForPlane,
-            const Array3 &segmentNormalOrientationForPlane,
-            const Array3 &projectionPointVertexNorms,
-            const Array3 &planeUnitNormal,
-            double planeDistance,
-            double planeNormalOrientation) {
+    std::pair<double, Array3> detail::computeSingularityTerms(const Array3Triplet &segmentVectorsForPlane,
+                                                              const Array3 &segmentNormalOrientationForPlane,
+                                                              const Array3 &projectionPointVertexNorms,
+                                                              const Array3 &planeUnitNormal, double planeDistance,
+                                                              double planeNormalOrientation) {
         //1. Case: If all sigma_pq for a given plane p are 1.0 then P' lies inside the plane S_p
         if (std::all_of(segmentNormalOrientationForPlane.cbegin(), segmentNormalOrientationForPlane.cend(),
                         [](const double sigma) { return sigma == 1.0; })) {
@@ -565,9 +363,9 @@ namespace polyhedralGravity {
                    projectionPointVertexNorms[j] < segmentVectorNorm;
         })) {
             using namespace util;
-            return std::make_pair(
-                    -1.0 * util::PI * planeDistance,                                //sing alpha = -pi*h_p
-                    planeUnitNormal * (-1.0 * util::PI * planeNormalOrientation));  //sing beta  = -pi*sigma_p*N_p
+            return std::make_pair(-1.0 * util::PI * planeDistance,                                //sing alpha = -pi*h_p
+                                  planeUnitNormal *
+                                  (-1.0 * util::PI * planeNormalOrientation));  //sing beta  = -pi*sigma_p*N_p
         }
         //3. Case If sigma_pq == 0 AND norm(P' - v1) < 0 || norm(P' - v2) < 0
         // then P' is located at one of G_p's vertices
@@ -598,21 +396,18 @@ namespace polyhedralGravity {
             const Array3 &g2 = r1Norm == 0.0 ? segmentVectorsForPlane[(j + 1) % 3] : segmentVectorsForPlane[j];
             // theta = arcos((G_2 * -G_1) / (|G_2| * |G_1|))
             const double gdot = dot(g1 * -1.0, g2);
-            const double theta =
-                    gdot == 0.0 ? util::PI_2 : std::acos(gdot / (euclideanNorm(g1) * euclideanNorm(g2)));
-            return std::make_pair(
-                    -1.0 * theta * planeDistance,                               //sing alpha = -theta*h_p
-                    planeUnitNormal * (-1.0 * theta * planeNormalOrientation)); //sing beta  = -theta*sigma_p*N_p
+            const double theta = gdot == 0.0 ? util::PI_2 : std::acos(gdot / (euclideanNorm(g1) * euclideanNorm(g2)));
+            return std::make_pair(-1.0 * theta * planeDistance,                               //sing alpha = -theta*h_p
+                                  planeUnitNormal *
+                                  (-1.0 * theta * planeNormalOrientation)); //sing beta  = -theta*sigma_p*N_p
         }
         //4. Case Otherwise P' is located outside the plane S_p and then the singularity equals zero
-        return std::make_pair(
-                0.0,                                                            //sing alpha = 0
-                Array3{0.0, 0.0, 0.0});                                         //sing beta  = 0
+        return std::make_pair(0.0,                                                            //sing alpha = 0
+                              Array3{0.0, 0.0, 0.0});                                         //sing beta  = 0
     }
 
-    Array3 GravityModel::detail::computeNormsOfProjectionPointAndVertices(
-            const Array3 &orthogonalProjectionPointOnPlane,
-            const Array3Triplet &face) {
+    Array3 detail::computeNormsOfProjectionPointAndVertices(const Array3 &orthogonalProjectionPointOnPlane,
+                                                            const Array3Triplet &face) {
         using namespace util;
         return {euclideanNorm(orthogonalProjectionPointOnPlane - face[0]),
                 euclideanNorm(orthogonalProjectionPointOnPlane - face[1]),
diff --git a/src/polyhedralGravity/calculation/GravityModel.h b/src/polyhedralGravity/calculation/GravityModel.h
index f464d46..021d71d 100644
--- a/src/polyhedralGravity/calculation/GravityModel.h
+++ b/src/polyhedralGravity/calculation/GravityModel.h
@@ -4,8 +4,13 @@
 #include <array>
 #include <vector>
 #include <algorithm>
+#include <variant>
+
+#include "polyhedralGravity/input/TetgenAdapter.h"
 #include "polyhedralGravity/model/Polyhedron.h"
 #include "polyhedralGravity/model/GravityModelData.h"
+#include "polyhedralGravity/calculation/GravityEvaluable.h"
+#include "polyhedralGravity/calculation/PolyhedronTransform.h"
 #include "polyhedralGravity/util/UtilityConstants.h"
 #include "polyhedralGravity/util/UtilityContainer.h"
 #include "spdlog/spdlog.h"
@@ -19,257 +24,219 @@
 #include "polyhedralGravity/output/Logging.h"
 #include "xsimd/xsimd.hpp"
 
-namespace polyhedralGravity {
+/**
+ * Namespace containing the methods used to evaluate the polyhedrale Gravity Model
+ * @note Naming scheme corresponds to the following:
+ * evaluate()           --> main Method for evaluating the gravity model
+ * *()                  --> Methods calculating one property for the evaluation
+ */
+namespace polyhedralGravity::GravityModel {
 
     /**
-     * Namespace containing the methods used to evaluate the polyhedrale Gravity Model
-     * @note Naming scheme corresponds to the following:
-     * evaluate()           --> main Method for evaluating the gravity model
-     * *()                  --> Methods calculating one property for the evaluation
+     * Evaluates the polyhedrale gravity model for a given constant density polyhedron at computation
+     * point P.
+     * @param polyhedron - the polyhedron consisting of vertices and triangular faces
+     * @param density - the constant density in [kg/m^3]
+     * @param computationPoint - the computation Point P
+     * @param parallel - whether to evaluate in parallel or serial
+     * @return the GravityModelResult containing the potential, the acceleration and the change of acceleration
+     * at computation Point P
      */
-    namespace GravityModel {
+    GravityModelResult evaluate(const PolyhedronOrSource &polyhedron, double density, const Array3 &computationPoint, bool parallel = true);
+
+    /**
+     * Evaluates the polyhedral gravity model for a given constant density polyhedron at multiple computation
+     * points.
+     * @param polyhedron - the polyhedron consisting of vertices and triangular faces
+     * @param density - the constant density in [kg/m^3]
+     * @param computationPoints - vector of computation points
+     * @param parallel - whether to evaluate in parallel or serial
+     * @return the GravityModelResult containing the potential, the acceleration and the change of acceleration
+     * foreach computation Point P
+     */
+    std::vector<GravityModelResult>
+    evaluate(const PolyhedronOrSource &polyhedron, double density, const std::vector<Array3> &computationPoints, bool parallel = true);
+
+    namespace detail {
 
         /**
-         * Evaluates the polyhedrale gravity model for a given constant density polyhedron at computation
-         * point P.
-         * @param polyhedron - the polyhedron consisting of vertices and triangular faces
-         * @param density - the constant density in [kg/m^3]
-         * @param computationPoint - the computation Point P
-         * @return the GravityModelResult containing the potential, the acceleration and the change of acceleration
-         * at computation Point P
+         * Computes the segment vectors G_ij for one plane of the polyhedron according to Tsoulis (18).
+         * The segment vectors G_ij represent the vector from one vertex of the face to the neighboring vertex and
+         * depict every line segment of the triangular face (A-B-C)
+         * @param vertex0 - the first vertex A
+         * @param vertex1 - the second vertex B
+         * @param vertex2 - the third vertex C
+         * @return the segment vectors for a plane
          */
-        GravityModelResult evaluate(
-                const Polyhedron &polyhedron,
-                double density,
-                const Array3 &computationPoint);
+        Array3Triplet buildVectorsOfSegments(const Array3 &vertex0, const Array3 &vertex1, const Array3 &vertex2);
 
         /**
-         * Evaluates the polyhedral gravity model for a given constant density polyhedron at multiple computation
-         * points.
-         * @param polyhedron - the polyhedron consisting of vertices and triangular faces
-         * @param density - the constant density in [kg/m^3]
-         * @param computationPoints - vector of computation points
-         * @return the GravityModelResult containing the potential, the acceleration and the change of acceleration
-         * foreach computation Point P
+         * Computes the plane unit normal N_p for one plane p of the polyhedron according to Tsoulis (19).
+         * The plane unit normal is the outward pointing normal of the face from the polyhedron.
+         * @param segmentVector1 - first edge
+         * @param segmentVector2 - second edge
+         * @return plane unit normal
          */
-        std::vector<GravityModelResult> evaluate(
-                const Polyhedron &polyhedron,
-                double density,
-                const std::vector<Array3> &computationPoints);
-
+        Array3 buildUnitNormalOfPlane(const Array3 &segmentVector1, const Array3 &segmentVector2);
 
         /**
-         * An iterator transforming the polyhedron's coordinates on demand by a given offset.
-         * This function returns a pair of transform iterators (first --> begin(), second --> end()).
-         * @param polyhedron - reference to the polyhedron
-         * @param offset - the offset to apply
-         * @return pair of transform iterators
+         * Computes the segment unit normals n_pq for one plane p of the polyhedron according to Tsoulis (20).
+         * The segment unit normal n_pq represent the normal of one line segment of a polyhedrale face.
+         * @param segmentVectors - the segment vectors of the face G_p(0-2)
+         * @param planeUnitNormal - the plane unit normal N_p
+         * @return segment unit normals n_pq for plane p with q = {0, 1, 2}
          */
-        inline auto transformPolyhedron(const Polyhedron &polyhedron, const Array3 &offset = {0.0, 0.0, 0.0}) {
-            //The offset is captured by value to ensure its lifetime!
-            const auto lambdaOffsetApplication = [&polyhedron, offset](
-                    const std::array<size_t, 3> &face) -> Array3Triplet {
-                using namespace util;
-                return {polyhedron.getVertex(face[0]) - offset,
-                        polyhedron.getVertex(face[1]) - offset,
-                        polyhedron.getVertex(face[2]) - offset};
-            };
-            auto first = thrust::make_transform_iterator(polyhedron.getFaces().begin(), lambdaOffsetApplication);
-            auto last = thrust::make_transform_iterator(polyhedron.getFaces().end(), lambdaOffsetApplication);
-            return std::make_pair(first, last);
-        }
-
-
-        namespace detail {
-
-            /**
-             * Computes the segment vectors G_ij for one plane of the polyhedron according to Tsoulis (18).
-             * The segment vectors G_ij represent the vector from one vertex of the face to the neighboring vertex and
-             * depict every line segment of the triangular face (A-B-C)
-             * @param vertex0 - the first vertex A
-             * @param vertex1 - the second vertex B
-             * @param vertex2 - the third vertex C
-             * @return the segment vectors for a plane
-             */
-            Array3Triplet
-            buildVectorsOfSegments(const Array3 &vertex0, const Array3 &vertex1, const Array3 &vertex2);
+        Array3Triplet buildUnitNormalOfSegments(const Array3Triplet &segmentVectors, const Array3 &planeUnitNormal);
 
-            /**
-             * Computes the plane unit normal N_p for one plane p of the polyhedron according to Tsoulis (19).
-             * The plane unit normal is the outward pointing normal of the face from the polyhedron.
-             * @param segmentVector1 - first edge
-             * @param segmentVector2 - second edge
-             * @return plane unit normal
-             */
-            Array3 buildUnitNormalOfPlane(const Array3 &segmentVector1, const Array3 &segmentVector2);
-
-            /**
-             * Computes the segment unit normals n_pq for one plane p of the polyhedron according to Tsoulis (20).
-             * The segment unit normal n_pq represent the normal of one line segment of a polyhedrale face.
-             * @param segmentVectors - the segment vectors of the face G_p(0-2)
-             * @param planeUnitNormal - the plane unit normal N_p
-             * @return segment unit normals n_pq for plane p with q = {0, 1, 2}
-             */
-            Array3Triplet
-            buildUnitNormalOfSegments(const Array3Triplet &segmentVectors, const Array3 &planeUnitNormal);
-
-            /**
-             * Computes the plane unit normal orientation/ direction sigma_p for one plane p of the polyhedron
-             * according to Tsoulis (21).
-             * The plane unit normal orientation values represents the relative position of computation point P
-             * with respect to the pointing direction of N_p. E. g. if N_p points to the half-space containing P, the
-             * inner product of N_p and -G_i1 will be positive, leading to a negative sigma_p.
-             * If sigma_p is zero than P and P' lie geometrically in the same plane --> P == P'.
-             * @param planeUnitNormal - the plane unit normal N_p
-             * @param vertex0 - the first vertex of the plane
-             * @return plane normal orientation
-             */
-            double computeUnitNormalOfPlaneDirection(const Array3 &planeUnitNormal, const Array3 &vertex0);
+        /**
+         * Computes the plane unit normal orientation/ direction sigma_p for one plane p of the polyhedron
+         * according to Tsoulis (21).
+         * The plane unit normal orientation values represents the relative position of computation point P
+         * with respect to the pointing direction of N_p. E. g. if N_p points to the half-space containing P, the
+         * inner product of N_p and -G_i1 will be positive, leading to a negative sigma_p.
+         * If sigma_p is zero than P and P' lie geometrically in the same plane --> P == P'.
+         * @param planeUnitNormal - the plane unit normal N_p
+         * @param vertex0 - the first vertex of the plane
+         * @return plane normal orientation
+         */
+        double computeUnitNormalOfPlaneDirection(const Array3 &planeUnitNormal, const Array3 &vertex0);
 
-            /**
-             * Calculates the Hessian Plane form spanned by three given points p, q, and r.
-             * @param p - first point on the plane
-             * @param q - second point on the plane
-             * @param r - third point on the plane
-             * @return HessianPlane
-             * @related Cross-Product method https://tutorial.math.lamar.edu/classes/calciii/eqnsofplanes.aspx
-             */
-            HessianPlane computeHessianPlane(const Array3 &p, const Array3 &q, const Array3 &r);
+        /**
+         * Calculates the Hessian Plane form spanned by three given points p, q, and r.
+         * @param p - first point on the plane
+         * @param q - second point on the plane
+         * @param r - third point on the plane
+         * @return HessianPlane
+         * @related Cross-Product method https://tutorial.math.lamar.edu/classes/calciii/eqnsofplanes.aspx
+         */
+        HessianPlane computeHessianPlane(const Array3 &p, const Array3 &q, const Array3 &r);
 
-            /**
-             * Calculates the (plane) distances h_p of computation point P to the plane S_p given in Hessian Form
-             * according to the following equation:
-             * h_p = D / sqrt(A^2+B^2+C^2)
-             * @param hessianPlane - Hessian Plane Form of S_p
-             * @return plane distance h_p
-             */
-            double distanceBetweenOriginAndPlane(const HessianPlane &hessianPlane);
+        /**
+         * Calculates the (plane) distances h_p of computation point P to the plane S_p given in Hessian Form
+         * according to the following equation:
+         * h_p = D / sqrt(A^2+B^2+C^2)
+         * @param hessianPlane - Hessian Plane Form of S_p
+         * @return plane distance h_p
+         */
+        double distanceBetweenOriginAndPlane(const HessianPlane &hessianPlane);
 
-            /**
-             * Computes P' for a given plane p according to equation (22) of Tsoulis paper.
-             * P' is the orthogonal projection of the computation point P onto the plane S_p.
-             * @param planeUnitNormal - the plane unit normal N_p
-             * @param planeDistance - the distance from P to the plane h_p
-             * @param hessianPlane - the Hessian Plane Form
-             * @return P' for this plane
-             */
-            Array3 projectPointOrthogonallyOntoPlane(
-                    const Array3 &planeUnitNormal,
-                    double planeDistance,
-                    const HessianPlane &hessianPlane);
+        /**
+         * Computes P' for a given plane p according to equation (22) of Tsoulis paper.
+         * P' is the orthogonal projection of the computation point P onto the plane S_p.
+         * @param planeUnitNormal - the plane unit normal N_p
+         * @param planeDistance - the distance from P to the plane h_p
+         * @param hessianPlane - the Hessian Plane Form
+         * @return P' for this plane
+         */
+        Array3 projectPointOrthogonallyOntoPlane(const Array3 &planeUnitNormal, double planeDistance,
+                                                 const HessianPlane &hessianPlane);
 
-            /**
-             * Computes the segment normal orientations/ directions sigma_pq for a given plane p.
-             * If sigma_pq is negative, this denotes that n_pq points to the half-plane containing P'. Nn case
-             * sigma_pq is positive, P' resides in the other half-plane and if sigma_pq is zero, then P' lies directly
-             * on the segment pq.
-             * @param vertices - the vertices of this plane
-             * @param projectionPointOnPlane - the projection point P' for this plane
-             * @param segmentUnitNormalsForPlane - the segment unit normals sigma_pq for this plane
-             * @return the segment normal orientations for the plane p
-             */
-            Array3 computeUnitNormalOfSegmentsDirections(const Array3Triplet &vertices,
-                                                         const Array3 &projectionPointOnPlane,
-                                                         const Array3Triplet &segmentUnitNormalsForPlane);
+        /**
+         * Computes the segment normal orientations/ directions sigma_pq for a given plane p.
+         * If sigma_pq is negative, this denotes that n_pq points to the half-plane containing P'. Nn case
+         * sigma_pq is positive, P' resides in the other half-plane and if sigma_pq is zero, then P' lies directly
+         * on the segment pq.
+         * @param vertices - the vertices of this plane
+         * @param projectionPointOnPlane - the projection point P' for this plane
+         * @param segmentUnitNormalsForPlane - the segment unit normals sigma_pq for this plane
+         * @return the segment normal orientations for the plane p
+         */
+        Array3
+        computeUnitNormalOfSegmentsDirections(const Array3Triplet &vertices, const Array3 &projectionPointOnPlane,
+                                              const Array3Triplet &segmentUnitNormalsForPlane);
 
-            /**
-             * Computes the orthogonal projection Points P'' foreach segment q of a given plane p.
-             * @param projectionPointOnPlane - the projection Point P'
-             * @param segmentNormalOrientations - the segment normal orientations sigma_pq for this plane p
-             * @param face - the vertices of the plane p
-             * @return the orthogonal projection points of P on the segment P'' foreach segment q of p
-             */
-            Array3Triplet projectPointOrthogonallyOntoSegments(
-                    const Array3 &projectionPointOnPlane,
-                    const Array3 &segmentNormalOrientations,
-                    const Array3Triplet &face);
+        /**
+         * Computes the orthogonal projection Points P'' foreach segment q of a given plane p.
+         * @param projectionPointOnPlane - the projection Point P'
+         * @param segmentNormalOrientations - the segment normal orientations sigma_pq for this plane p
+         * @param face - the vertices of the plane p
+         * @return the orthogonal projection points of P on the segment P'' foreach segment q of p
+         */
+        Array3Triplet projectPointOrthogonallyOntoSegments(const Array3 &projectionPointOnPlane,
+                                                           const Array3 &segmentNormalOrientations,
+                                                           const Array3Triplet &face);
 
-            /**
-             * Calculates the point P'' for a given Segment consisting of vertices v1 and v2 and the orthogonal projection
-             * point P' for the plane consisting of those vertices. Solves the three equations given in (24), (25) and (26).
-             * @param vertex1 - first endpoint of segment
-             * @param vertex2 - second endpoint of segment
-             * @param orthogonalProjectionPointOnPlane - the orthogonal projection P' of P on this plane
-             * @return P'' for this segment
-             * @note If sigma_pq is zero then P'' == P', this is not checked by this method, but has to be assured first
-             */
-            Array3 projectPointOrthogonallyOntoSegment(const Array3 &vertex1, const Array3 &vertex2,
-                                                       const Array3 &orthogonalProjectionPointOnPlane);
+        /**
+         * Calculates the point P'' for a given Segment consisting of vertices v1 and v2 and the orthogonal projection
+         * point P' for the plane consisting of those vertices. Solves the three equations given in (24), (25) and (26).
+         * @param vertex1 - first endpoint of segment
+         * @param vertex2 - second endpoint of segment
+         * @param orthogonalProjectionPointOnPlane - the orthogonal projection P' of P on this plane
+         * @return P'' for this segment
+         * @note If sigma_pq is zero then P'' == P', this is not checked by this method, but has to be assured first
+         */
+        Array3 projectPointOrthogonallyOntoSegment(const Array3 &vertex1, const Array3 &vertex2,
+                                                   const Array3 &orthogonalProjectionPointOnPlane);
 
-            /**
-             * Computes the (segment) distances h_pq between P' for a given plane p and P'' for a given segment q of plane p.
-             * @param orthogonalProjectionPointOnPlane - the orthogonal projection point P' for p
-             * @param orthogonalProjectionPointOnSegments - the orthogonal projection points P'' for each segment q of p
-             * @return distances h_pq for plane p
-             */
-            Array3 distancesBetweenProjectionPoints(
-                    const Array3 &orthogonalProjectionPointOnPlane,
-                    const Array3Triplet &orthogonalProjectionPointOnSegments);
+        /**
+         * Computes the (segment) distances h_pq between P' for a given plane p and P'' for a given segment q of plane p.
+         * @param orthogonalProjectionPointOnPlane - the orthogonal projection point P' for p
+         * @param orthogonalProjectionPointOnSegments - the orthogonal projection points P'' for each segment q of p
+         * @return distances h_pq for plane p
+         */
+        Array3 distancesBetweenProjectionPoints(const Array3 &orthogonalProjectionPointOnPlane,
+                                                const Array3Triplet &orthogonalProjectionPointOnSegments);
 
-            /**
-             * Computes the 3D distances l1_pq and l2_pq between the computation point P and the line
-             * segment endpoints of each polyhedral segment for one plane.
-             * Computes the 1D distances s1_pq and s2_pq between orthogonal projection of P on the line
-             * segment P''_pq and the line segment endpoints for each polyhedral segment for one plane
-             * @param segmentVectorsForPlane - the segment vectors G_pq for plane p
-             * @param orthogonalProjectionPointsOnSegmentForPlane - the orthogonal projection Points P'' for plane p
-             * @param face - the vertices of plane p
-             * @return distances l1_pq and l2_pq and s1_pq and s2_pq foreach segment q of plane p
-             */
-            std::array<Distance, 3> distancesToSegmentEndpoints(
-                    const Array3Triplet &segmentVectorsForPlane,
-                    const Array3Triplet &orthogonalProjectionPointsOnSegmentForPlane,
-                    const Array3Triplet &face);
+        /**
+         * Computes the 3D distances l1_pq and l2_pq between the computation point P and the line
+         * segment endpoints of each polyhedral segment for one plane.
+         * Computes the 1D distances s1_pq and s2_pq between orthogonal projection of P on the line
+         * segment P''_pq and the line segment endpoints for each polyhedral segment for one plane
+         * @param segmentVectorsForPlane - the segment vectors G_pq for plane p
+         * @param orthogonalProjectionPointsOnSegmentForPlane - the orthogonal projection Points P'' for plane p
+         * @param face - the vertices of plane p
+         * @return distances l1_pq and l2_pq and s1_pq and s2_pq foreach segment q of plane p
+         */
+        std::array<Distance, 3> distancesToSegmentEndpoints(const Array3Triplet &segmentVectorsForPlane,
+                                                            const Array3Triplet &orthogonalProjectionPointsOnSegmentForPlane,
+                                                            const Array3Triplet &face);
 
-            /**
-             * Calculates the Transcendental Expressions LN_pq and AN_pq for every line segment of the polyhedron for
-             * a given plane p.
-             * LN_pq is calculated according to (14) using the natural logarithm and AN_pq is calculated according
-             * to (15) using the arctan.
-             * @param distancesForPlane - the distances l1, l2, s1, s2 foreach segment q of plane p
-             * @param planeDistance - the plane distance h_p for plane p
-             * @param segmentDistancesForPlane - the segment distance h_pq for segment q of plane p
-             * @param segmentNormalOrientationsForPlane - the segment normal orientations n_pq for a plane p
-             * @param orthogonalProjectionPointOnPlane - the orthogonal projection point P' for plane p
-             * @param face - the vertices of plane p
-             * @return LN_pq and AN_pq foreach segment q of plane p
-             */
-            std::array<TranscendentalExpression, 3>
-            computeTranscendentalExpressions(const std::array<Distance, 3> &distancesForPlane,
-                                             double planeDistance, const Array3 &segmentDistancesForPlane,
-                                             const Array3 &segmentNormalOrientationsForPlane,
-                                             const Array3 &projectionPointVertexNorms);
+        /**
+         * Calculates the Transcendental Expressions LN_pq and AN_pq for every line segment of the polyhedron for
+         * a given plane p.
+         * LN_pq is calculated according to (14) using the natural logarithm and AN_pq is calculated according
+         * to (15) using the arctan.
+         * @param distancesForPlane - the distances l1, l2, s1, s2 foreach segment q of plane p
+         * @param planeDistance - the plane distance h_p for plane p
+         * @param segmentDistancesForPlane - the segment distance h_pq for segment q of plane p
+         * @param segmentNormalOrientationsForPlane - the segment normal orientations n_pq for a plane p
+         * @param orthogonalProjectionPointOnPlane - the orthogonal projection point P' for plane p
+         * @param face - the vertices of plane p
+         * @return LN_pq and AN_pq foreach segment q of plane p
+         */
+        std::array<TranscendentalExpression, 3>
+        computeTranscendentalExpressions(const std::array<Distance, 3> &distancesForPlane, double planeDistance,
+                                         const Array3 &segmentDistancesForPlane,
+                                         const Array3 &segmentNormalOrientationsForPlane,
+                                         const Array3 &projectionPointVertexNorms);
 
-            /**
-             * Calculates the singularities (correction) terms according to the Flow text for a given plane p.
-             * @param segmentVectorsForPlane - the segment vectors for a given plane
-             * @param segmentNormalOrientationForPlane - the segment orientation sigma_pq
-             * @param projectionPointVertexNorms - the projection point P'
-             * @param planeUnitNormal - the plane unit normal N_p
-             * @param planeDistance - the plane distance h_p
-             * @param planeNormalOrientation - the plane normal orientation sigma_p
-             * @param face - the vertices of plane p
-             * @return the singularities for a plane p
-             */
-            std::pair<double, Array3> computeSingularityTerms(const Array3Triplet &segmentVectorsForPlane,
-                                                              const Array3 &segmentNormalOrientationForPlane,
-                                                              const Array3 &projectionPointVertexNorms,
-                                                              const Array3 &planeUnitNormal,
-                                                              double planeDistance,
-                                                              double planeNormalOrientation);
+        /**
+         * Calculates the singularities (correction) terms according to the Flow text for a given plane p.
+         * @param segmentVectorsForPlane - the segment vectors for a given plane
+         * @param segmentNormalOrientationForPlane - the segment orientation sigma_pq
+         * @param projectionPointVertexNorms - the projection point P'
+         * @param planeUnitNormal - the plane unit normal N_p
+         * @param planeDistance - the plane distance h_p
+         * @param planeNormalOrientation - the plane normal orientation sigma_p
+         * @param face - the vertices of plane p
+         * @return the singularities for a plane p
+         */
+        std::pair<double, Array3> computeSingularityTerms(const Array3Triplet &segmentVectorsForPlane,
+                                                          const Array3 &segmentNormalOrientationForPlane,
+                                                          const Array3 &projectionPointVertexNorms,
+                                                          const Array3 &planeUnitNormal, double planeDistance,
+                                                          double planeNormalOrientation);
 
 
-            /**
-             * Computes the L2 norms of the orthogonal projection point P' on a plane p with each vertex of that plane p.
-             * The values are later used to determine if P' is situated at a vertex.
-             * @param orthogonalProjectionPointOnPlane - the orthogonal projection point P'
-             * @param face - the vertices of plane p
-             * @return the norms of p and each vertex
-             */
-            Array3 computeNormsOfProjectionPointAndVertices(
-                    const Array3 &orthogonalProjectionPointOnPlane,
-                    const Array3Triplet &face);
-        }
-    };
+        /**
+         * Computes the L2 norms of the orthogonal projection point P' on a plane p with each vertex of that plane p.
+         * The values are later used to determine if P' is situated at a vertex.
+         * @param orthogonalProjectionPointOnPlane - the orthogonal projection point P'
+         * @param face - the vertices of plane p
+         * @return the norms of p and each vertex
+         */
+        Array3 computeNormsOfProjectionPointAndVertices(const Array3 &orthogonalProjectionPointOnPlane,
+                                                        const Array3Triplet &face);
 
-}
\ No newline at end of file
+    }
+}
diff --git a/src/polyhedralGravity/calculation/MeshChecking.cpp b/src/polyhedralGravity/calculation/MeshChecking.cpp
index 79092fe..38ac5f3 100644
--- a/src/polyhedralGravity/calculation/MeshChecking.cpp
+++ b/src/polyhedralGravity/calculation/MeshChecking.cpp
@@ -4,7 +4,7 @@ namespace polyhedralGravity::MeshChecking {
 
     bool checkNormalsOutwardPointing(const Polyhedron &polyhedron) {
         using namespace util;
-        auto it = GravityModel::transformPolyhedron(polyhedron);
+        auto it = transformPolyhedron(polyhedron);
         // All normals have to point outwards (intersect the polyhedron even times)
         return thrust::transform_reduce(
                 thrust::device,
@@ -27,7 +27,7 @@ namespace polyhedralGravity::MeshChecking {
     }
 
     bool checkTrianglesNotDegenerated(const Polyhedron &polyhedron) {
-        auto it = GravityModel::transformPolyhedron(polyhedron);
+        auto it = transformPolyhedron(polyhedron);
         // All triangles surface area needs to be greater than zero
         return thrust::transform_reduce(
                 thrust::device,
@@ -38,7 +38,7 @@ namespace polyhedralGravity::MeshChecking {
 
     size_t
     detail::countRayPolyhedronIntersections(const Array3 &rayOrigin, const Array3 &rayVector, const Polyhedron &polyhedron) {
-        auto it = GravityModel::transformPolyhedron(polyhedron);
+        auto it = transformPolyhedron(polyhedron);
         std::set<Array3> intersections{};
         // Count every triangular face which is intersected by the ray
         std::for_each(it.first, it.second, [&rayOrigin, &rayVector, &intersections](const Array3Triplet &triangle) {
diff --git a/src/polyhedralGravity/calculation/PolyhedronTransform.h b/src/polyhedralGravity/calculation/PolyhedronTransform.h
new file mode 100644
index 0000000..c92614f
--- /dev/null
+++ b/src/polyhedralGravity/calculation/PolyhedronTransform.h
@@ -0,0 +1,29 @@
+#pragma once
+
+#include <utility>
+
+#include "polyhedralGravity/model/Polyhedron.h"
+#include "polyhedralGravity/model/GravityModelData.h"
+#include "thrust/iterator/transform_iterator.h"
+
+namespace polyhedralGravity {
+
+    /**
+     * An iterator transforming the polyhedron's coordinates on demand by a given offset.
+     * This function returns a pair of transform iterators (first --> begin(), second --> end()).
+     * @param polyhedron - reference to the polyhedron
+     * @param offset - the offset to apply
+     * @return pair of transform iterators
+     */
+    inline auto transformPolyhedron(const Polyhedron &polyhedron, const Array3 &offset = {0.0, 0.0, 0.0}) {
+        //The offset is captured by value to ensure its lifetime!
+        const auto lambdaOffsetApplication = [&polyhedron, offset](const IndexArray3 &face) -> Array3Triplet {
+            using namespace util;
+            return {polyhedron.getVertex(face[0]) - offset, polyhedron.getVertex(face[1]) - offset,
+                    polyhedron.getVertex(face[2]) - offset};
+        };
+        auto first = thrust::make_transform_iterator(polyhedron.getFaces().begin(), lambdaOffsetApplication);
+        auto last = thrust::make_transform_iterator(polyhedron.getFaces().end(), lambdaOffsetApplication);
+        return std::make_pair(first, last);
+    }
+}
\ No newline at end of file
diff --git a/src/polyhedralGravity/model/GravityModelData.h b/src/polyhedralGravity/model/GravityModelData.h
index 3ce07fc..3cb1cfb 100644
--- a/src/polyhedralGravity/model/GravityModelData.h
+++ b/src/polyhedralGravity/model/GravityModelData.h
@@ -4,23 +4,51 @@
 #include <ostream>
 #include <cmath>
 #include <algorithm>
+#include <tuple>
 #include "polyhedralGravity/util/UtilityContainer.h"
 #include "polyhedralGravity/util/UtilityConstants.h"
 
 namespace polyhedralGravity {
 
     /**
-     * Alias for an array of size 3
+     * Alias for an array of size 3 (double)
      * @example for x, y, z coordinates.
      */
     using Array3 = std::array<double, 3>;
 
+    /**
+     * Alias for an array of size 3 (size_t)
+     * @example for the vertex indices in a triangular face.
+     */
+    using IndexArray3 = std::array<size_t, 3>;
+
+    /**
+     * Alias for an array of size 6
+     * @example for xx, yy, zz, xy, xz, yz second derivatives.
+     */
+    using Array6 = std::array<double, 6>;
+
     /**
      * Alias for a triplet of arrays of size 3
      * @example for the segment of a triangular face
      */
     using Array3Triplet = std::array<Array3, 3>;
 
+    /**
+     * Contains in the order of the tuple:
+     * The gravitational potential in [m^2/s^2] <--> [J/kg] at point P.
+     * @related Equation (1) and (11) of Tsoulis Paper, here referred as V
+     *
+     * The first order derivatives of the gravitational potential in [m/s^2].
+     * The array contains the derivatives depending on the coordinates x-y-z in this order.
+     * @related Equation (2) and (12) of Tsoulis Paper, here referred as Vx, Vy, Vz
+     *
+     * The second order derivatives or also called gradiometric Tensor in [1/s^2].
+     * The array contains the second order derivatives in the following order xx, yy, zz, xy, xz, yz.
+     * @related Equation (3) and (13) of Tsoulis Paper, here referred as Vxx, Vyy, Vzz, Vxy, Vxz, Vyz
+     */
+    using GravityModelResult = std::tuple<double, Array3, Array6>;
+
     /**
      * Contains the 3D distances l1_pq and l2_pq between P and the endpoints of segment pq and
      * the 1D distances s1_pq and s2_pq between P'' and the segment endpoints.
@@ -55,10 +83,7 @@ namespace polyhedralGravity {
          * @note Just used for testing purpose
          */
         bool operator==(const Distance &rhs) const {
-            return l1 == rhs.l1 &&
-                   l2 == rhs.l2 &&
-                   s1 == rhs.s1 &&
-                   s2 == rhs.s2;
+            return l1 == rhs.l1 && l2 == rhs.l2 && s1 == rhs.s1 && s2 == rhs.s2;
         }
 
         /**
@@ -112,8 +137,7 @@ namespace polyhedralGravity {
          * @note Just used for testing purpose
          */
         bool operator==(const TranscendentalExpression &rhs) const {
-            return ln == rhs.ln &&
-                   an == rhs.an;
+            return ln == rhs.ln && an == rhs.an;
         }
 
         /**
@@ -174,10 +198,7 @@ namespace polyhedralGravity {
          * @note Just used for testing purpose
          */
         bool operator==(const HessianPlane &rhs) const {
-            return a == rhs.a &&
-                   b == rhs.b &&
-                   c == rhs.c &&
-                   d == rhs.d;
+            return a == rhs.a && b == rhs.b && c == rhs.c && d == rhs.d;
         }
 
         /**
@@ -193,81 +214,4 @@ namespace polyhedralGravity {
         }
     };
 
-    /**
-     * A data structure containing the result of the polyhedral gravity model's evaluation.
-    */
-    struct GravityModelResult {
-
-        /**
-         * The gravitational potential in [m^2/s^2] <--> [J/kg] at point P.
-         * @related Equation (1) and (11) of Tsoulis Paper, here referred as V
-         */
-        double gravitationalPotential{};
-
-        /**
-         * The first order derivatives of the gravitational potential in [m/s^2].
-         * The array contains the derivatives depending on the coordinates x-y-z in this order.
-         * @related Equation (2) and (12) of Tsoulis Paper, here referred as Vx, Vy, Vz
-         */
-        std::array<double, 3> acceleration{};
-
-        /**
-         * The second order derivatives or also called gradiometric Tensor in [1/s^2].
-         * The array contains the second order derivatives in the following order xx, yy, zz, xy, xz, yz.
-         * @related Equation (3) and (13) of Tsoulis Paper, here referred as Vxx, Vyy, Vzz, Vxy, Vxz, Vyz
-         */
-        std::array<double, 6> gradiometricTensor{};
-
-        /**
-         * Creates new empty Result
-         */
-        GravityModelResult() = default;
-
-        /**
-         * Creates a new Result with the given values
-         * @param gravitationalPotential - the potential
-         * @param gravitationalPotentialDerivative - the acceleration
-         * @param gradiometricTensor - the gradiometric tensor i.e. the second derivatives
-         */
-        GravityModelResult(double gravitationalPotential, const std::array<double, 3> &gravitationalPotentialDerivative,
-                           const std::array<double, 6> &gradiometricTensor)
-                : gravitationalPotential(gravitationalPotential),
-                  acceleration(gravitationalPotentialDerivative),
-                  gradiometricTensor(gradiometricTensor) {}
-
-        /**
-         * Pretty output of this struct on the given ostream.
-         * @param os - the ostream
-         * @param result - the GravityModelResult to output
-         * @return os
-         */
-        friend std::ostream &operator<<(std::ostream &os, const GravityModelResult &result) {
-            using util::operator<<;
-            os << "gravitationalPotential: " << result.gravitationalPotential << " acceleration: "
-               << result.acceleration << " gradiometricTensor: " << result.gradiometricTensor;
-            return os;
-        }
-
-        /**
-         * Compares the results to used epsilon and if one is smaller resets it to zero.
-         * This step is performed due to rounding errors when calculating with doubles.
-         */
-        inline void eliminateRoundingErrors(){
-            if (std::abs(gravitationalPotential) < util::EPSILON) {
-                gravitationalPotential = 0.0;
-            }
-            for (double &acc : acceleration) {
-                if (std::abs(acc) < util::EPSILON) {
-                    acc = 0.0;
-                }
-            }
-            for (double &secondDev : gradiometricTensor) {
-                if (std::abs(secondDev) < util::EPSILON) {
-                    secondDev = 0.0;
-                }
-            }
-        }
-
-    };
-
 }
diff --git a/src/polyhedralGravity/model/Polyhedron.h b/src/polyhedralGravity/model/Polyhedron.h
index cb96182..97df911 100644
--- a/src/polyhedralGravity/model/Polyhedron.h
+++ b/src/polyhedralGravity/model/Polyhedron.h
@@ -6,9 +6,28 @@
 #include <algorithm>
 #include <exception>
 #include <stdexcept>
+#include <tuple>
+#include <variant>
+#include <string>
+#include "GravityModelData.h"
 
 namespace polyhedralGravity {
 
+    /* Forward declaration of Polyhedron */
+    class Polyhedron;
+
+    /**
+     * Variant of possible polyhedron sources (composed of members, read from file), but not the polyhedron itself
+     * @example Utilized in the Python interface which does not expose the Polyhedron class
+     */
+    using PolyhedralSource = std::variant<std::tuple<std::vector<Array3>, std::vector<IndexArray3>>, std::vector<std::string>>;
+
+    /**
+     * Variant of possible polyhedral sources (direct, composed of members, read from file), including the polyhedron
+     * @example Utilized in the C++ implementation of the polyhedral-gravity model
+     */
+    using PolyhedronOrSource = std::variant<Polyhedron, std::tuple<std::vector<Array3>, std::vector<IndexArray3>>, std::vector<std::string>>;
+
     /**
      * Data structure containing the model data of one polyhedron. This includes nodes, edges (faces) and elements.
      * The index always starts with zero!
@@ -19,7 +38,7 @@ namespace polyhedralGravity {
          * A vector containing the vertices of the polyhedron.
          * Each node is an array of size three containing the xyz coordinates.
          */
-        const std::vector<std::array<double, 3>> _vertices;
+        const std::vector<Array3> _vertices;
 
         /**
          * A vector containing the faces (triangles) of the polyhedron.
@@ -28,7 +47,7 @@ namespace polyhedralGravity {
          * two nodes.
          * @example face consisting of {1, 2, 3} --> segments: {1, 2}, {2, 3}, {3, 1}
          */
-        const std::vector<std::array<size_t, 3>> _faces;
+        const std::vector<IndexArray3> _faces;
 
 
     public:
@@ -48,7 +67,7 @@ namespace polyhedralGravity {
          * @note ASSERTS PRE-CONDITION that the in the indexing in the faces vector starts with zero!
          * @throws runtime_error if no face contains the node zero indicating mathematical index
          */
-        Polyhedron(std::vector<std::array<double, 3>> nodes, std::vector<std::array<size_t, 3>> faces)
+        Polyhedron(std::vector<Array3> nodes, std::vector<IndexArray3> faces)
                 : _vertices{std::move(nodes)},
                   _faces{std::move(faces)} {
             //Checks that the node with index zero is actually used
@@ -61,6 +80,17 @@ namespace polyhedralGravity {
             }
         }
 
+        /**
+         * Creates a polyhedron from a tuple of nodes and faces.
+         * @param data - tuple of nodes and faces
+         *
+         * @note ASSERTS PRE-CONDITION that the in the indexing in the faces vector starts with zero!
+         * @throws runtime_error if no face contains the node zero indicating mathematical index
+         */
+        explicit Polyhedron(std::tuple<std::vector<Array3>, std::vector<IndexArray3>> data)
+                : Polyhedron(std::get<std::vector<Array3>>(data), std::get<std::vector<IndexArray3>>(data)) {}
+
+
         /**
          * Default destructor
          */
@@ -70,7 +100,7 @@ namespace polyhedralGravity {
          * Returns the vertices of this polyhedron
          * @return vector of cartesian coordinates
          */
-        [[nodiscard]] const std::vector<std::array<double, 3>> &getVertices() const {
+        [[nodiscard]] const std::vector<Array3> &getVertices() const {
             return _vertices;
         }
 
@@ -79,7 +109,7 @@ namespace polyhedralGravity {
          * @param index - size_t
          * @return cartesian coordinates of the vertex at index
          */
-        [[nodiscard]] const std::array<double, 3> &getVertex(size_t index) const {
+        [[nodiscard]] const Array3 &getVertex(size_t index) const {
             return _vertices[index];
         }
 
@@ -103,10 +133,19 @@ namespace polyhedralGravity {
          * Returns the triangular faces of this polyhedron
          * @return vector of triangular faces, where each element size_t references a vertex in the vertices vector
          */
-        [[nodiscard]] const std::vector<std::array<size_t, 3>> &getFaces() const {
+        [[nodiscard]] const std::vector<IndexArray3> &getFaces() const {
             return _faces;
         }
 
+        /**
+         * Returns the resolved face with its concrete cartesian coordinates at the given index.
+         * @param index - size_t
+         * @return triplet of vertices' cartesian coordinates
+         */
+        [[nodiscard]] Array3Triplet getFace(size_t index) const {
+            return {_vertices[_faces[index][0]], _vertices[_faces[index][1]], _vertices[_faces[index][2]]};
+        }
+
     };
 
 }
diff --git a/src/polyhedralGravity/output/CSVWriter.cpp b/src/polyhedralGravity/output/CSVWriter.cpp
index 3a74af8..dd8eb46 100644
--- a/src/polyhedralGravity/output/CSVWriter.cpp
+++ b/src/polyhedralGravity/output/CSVWriter.cpp
@@ -7,9 +7,7 @@ namespace polyhedralGravity {
         _logger->info("Point P,Potential [m^2/s^2],Acceleration [m/s^2],Second Derivative Gravity Tensor [1/s^2]");
         for (size_t i = 0; i < computationPoints.size() && i < gravityResults.size(); ++i) {
             const std::array<double, 3> &computationPoint = computationPoints[i];
-            const double &potential = gravityResults[i].gravitationalPotential;
-            const std::array<double, 3> &acceleration = gravityResults[i].acceleration;
-            const std::array<double, 6> &secondDerivative = gravityResults[i].gradiometricTensor;
+            const auto &[potential, acceleration, secondDerivative] = gravityResults[i];
             //Because of ADL the overload operator<< for arrays (in UtilityContainer) does not work here
             //@related https://en.cppreference.com/w/cpp/language/adl
             _logger->info("[{} {} {}],{},[{} {} {}],[{} {} {} {} {} {}]",
diff --git a/src/polyhedralGravity/util/UtilityContainer.h b/src/polyhedralGravity/util/UtilityContainer.h
index 806d510..98884d9 100644
--- a/src/polyhedralGravity/util/UtilityContainer.h
+++ b/src/polyhedralGravity/util/UtilityContainer.h
@@ -2,6 +2,7 @@
 
 #include <array>
 #include <numeric>
+#include <utility>
 #include <algorithm>
 #include <functional>
 #include <cmath>
@@ -14,8 +15,7 @@ namespace polyhedralGravity::util {
      * Alias for two-dimensional array with size M and N.
      * M is the major size.
      */
-    template<typename T, size_t M, size_t N>
-    using Matrix = std::array<std::array<T, N>, M>;
+    template<typename T, size_t M, size_t N> using Matrix = std::array<std::array<T, N>, M>;
 
     /**
      * Applies a binary function to elements of two containers piece by piece. The objects must
@@ -198,9 +198,9 @@ namespace polyhedralGravity::util {
      */
     template<typename T>
     T det(const Matrix<T, 3, 3> &matrix) {
-        return matrix[0][0] * matrix[1][1] * matrix[2][2] + matrix[0][1] * matrix[1][2] * matrix[2][0]
-               + matrix[0][2] * matrix[1][0] * matrix[2][1] - matrix[0][2] * matrix[1][1] * matrix[2][0]
-               - matrix[0][0] * matrix[1][2] * matrix[2][1] - matrix[0][1] * matrix[1][0] * matrix[2][2];
+        return matrix[0][0] * matrix[1][1] * matrix[2][2] + matrix[0][1] * matrix[1][2] * matrix[2][0] +
+               matrix[0][2] * matrix[1][0] * matrix[2][1] - matrix[0][2] * matrix[1][1] * matrix[2][0] -
+               matrix[0][0] * matrix[1][2] * matrix[2][1] - matrix[0][1] * matrix[1][0] * matrix[2][2];
     }
 
     /**
@@ -327,6 +327,35 @@ namespace polyhedralGravity::util {
         return std::abs(x - y) > maxExponentDifference;
     }
 
+    namespace detail {
+
+        /**
+         * Helper method for the tuple operator+, which expands the tuple into a parameter pack.
+         * @tparam Ts - types of the tuple
+         * @tparam Is - indices of the tuple
+         * @param t1 - first tuple
+         * @param t2 - second tuple
+         * @return a new tuple with the added values
+         */
+        template<typename... Ts, size_t... Is>
+        auto tuple_add(const std::tuple<Ts...> &t1, const std::tuple<Ts...> &t2, std::index_sequence<Is...>) {
+            return std::make_tuple(std::get<Is>(t1) + std::get<Is>(t2)...);
+        }
+
+    }
+
+    /**
+     * Adds the contents of two tuples of the same size and types with the operator +.
+     * @tparam Ts - types of the tuples
+     * @param t1 - first tuple
+     * @param t2 - second tuple
+     * @return a new tuple with the added values
+     */
+    template<typename... Ts>
+    auto operator+(const std::tuple<Ts...> &t1, const std::tuple<Ts...> &t2) {
+        return detail::tuple_add(t1, t2, std::index_sequence_for<Ts...>{});
+    }
+
     /**
      * Operator << for an array of any size.
      * @tparam T - type of the array, must have an << operator overload
diff --git a/python-binding/CMakeLists.txt b/src/polyhedralGravityPython/CMakeLists.txt
similarity index 55%
rename from python-binding/CMakeLists.txt
rename to src/polyhedralGravityPython/CMakeLists.txt
index 0ed7cf9..e48c4eb 100644
--- a/python-binding/CMakeLists.txt
+++ b/src/polyhedralGravityPython/CMakeLists.txt
@@ -1,27 +1,13 @@
-# Include pybind11
-include(pybind11)
-
-# Stores the file of the python interface in a variable
-set(pythonInterface ${CMAKE_CURRENT_SOURCE_DIR}/python-interface.cpp)
+# Exclude the YAML files from the python interface
+set(PYTHON_INTERFACE_SRC ${SRC})
 
 file(GLOB_RECURSE PYTHON_INTERFACE_SRC
-        "${PROJECT_SOURCE_DIR}/src/polyhedralGravity/*.h"
-        "${PROJECT_SOURCE_DIR}/src/polyhedralGravity/*.cpp")
-
+    "${PROJECT_SOURCE_DIR}/src/polyhedralGravity/*.cpp"
+    "${PROJECT_SOURCE_DIR}/src/polyhedralGravity/*.h"
+    )
 list(FILTER PYTHON_INTERFACE_SRC EXCLUDE REGEX ".*YAML.*")
 
-pybind11_add_module(polyhedral_gravity ${pythonInterface} ${PYTHON_INTERFACE_SRC})
-
-target_include_directories(polyhedral_gravity PUBLIC
-        "${PROJECT_SOURCE_DIR}/src"
-        )
-
-target_link_libraries(polyhedral_gravity PUBLIC
-        spdlog::spdlog
-        tetgen
-        xsimd
-        Thrust
-        )
+pybind11_add_module(polyhedral_gravity ${PYTHON_INTERFACE_SRC} PolyhedralGravityPython.cpp)
 
 # Checks if MinGW is used, in this case, the .dll of the the standard library need to be linked statically to the
 # Python interface since otherwise they won't be found
@@ -37,4 +23,4 @@ if(WIN32 AND ${CMAKE_CXX_COMPILER_ID} STREQUAL "GNU")
         -lwinpthread
         --no-whole-file
         )
-endif()
\ No newline at end of file
+endif()
diff --git a/src/polyhedralGravityPython/PolyhedralGravityPython.cpp b/src/polyhedralGravityPython/PolyhedralGravityPython.cpp
new file mode 100644
index 0000000..7160587
--- /dev/null
+++ b/src/polyhedralGravityPython/PolyhedralGravityPython.cpp
@@ -0,0 +1,133 @@
+#include <tuple>
+#include <variant>
+#include "pybind11/pybind11.h"
+#include "pybind11/stl.h"
+
+#include "polyhedralGravity/model/Polyhedron.h"
+#include "polyhedralGravity/model/GravityModelData.h"
+#include "polyhedralGravity/calculation/GravityModel.h"
+#include "polyhedralGravity/calculation/GravityEvaluable.h"
+#include "polyhedralGravity/calculation/MeshChecking.h"
+#include "polyhedralGravity/input/TetgenAdapter.h"
+
+
+namespace py = pybind11;
+
+PYBIND11_MODULE(polyhedral_gravity, m) {
+    using namespace polyhedralGravity;
+    m.doc() = "Computes the full gravity tensor for a given constant density polyhedron which consists of some "
+              "vertices and triangular faces at a given computation points";
+
+    m.def("evaluate", [](const PolyhedralSource &polyhedralSource, double density,
+                         const std::variant<Array3, std::vector<Array3>> &computationPoints,
+                         bool parallel) -> std::variant<GravityModelResult, std::vector<GravityModelResult>> {
+              if (std::holds_alternative<Array3>(computationPoints)) {
+                  return std::visit([&density, &computationPoints, parallel](const auto &source) {
+                      return GravityModel::evaluate(source, density, std::get<Array3>(computationPoints), parallel);
+                  }, polyhedralSource);
+              } else {
+                  return std::visit([&density, &computationPoints, parallel](const auto &source) {
+                      return GravityModel::evaluate(source, density, std::get<std::vector<Array3>>(computationPoints),
+                                                    parallel);
+                  }, polyhedralSource);
+              }
+          }, R"mydelimiter(
+            Evaluates the polyhedral gravity model for a given constant density polyhedron at a given computation point.
+
+            Args:
+                polyhedral_source:  The vertices & faces of the polyhedron as tuple or
+                                    the filenames of the files containing the vertices & faces as list of strings
+                density:            The constant density of the polyhedron in [kg/m^3]
+                computation_points: The computation points as tuple or list of points
+                parallel:           If true, the computation is done in parallel (default: true)
+
+            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("polyhedral_source"), py::arg("density"), py::arg("computation_points"),
+          py::arg("parallel") = true);
+
+    pybind11::class_<GravityEvaluable>(m, "GravityEvaluable", R"mydelimiter(
+                A class to evaluate the polyhedral gravity model for a given constant density polyhedron at a given computation point.
+                It provides a __call__ method to evaluate the polyhedral gravity model for computation points while
+                also caching the polyhedron data over the lifetime of the object.
+            )mydelimiter")
+            .def(pybind11::init<const PolyhedralSource &, double>(),
+                 R"mydelimiter(
+                    Creates a new GravityEvaluable for a given constant density polyhedron.
+                    It provides a __call__ method to evaluate the polyhedral gravity model for computation points while
+                    also caching the polyhedron data over the lifetime of the object.
+
+                    Args:
+                        polyhedral_source:  The vertices & faces of the polyhedron as tuple or
+                                            the filenames of the files containing the vertices & faces as list of strings
+                        density:            The constant density of the polyhedron in [kg/m^3]
+                    )mydelimiter",
+                 pybind11::arg("polyhedral_source"),
+                 pybind11::arg("density"))
+            .def("__repr__", &GravityEvaluable::toString)
+            .def("__call__", &GravityEvaluable::operator(),
+                 R"mydelimiter(
+                    Evaluates the polyhedral gravity model for a given constant density polyhedron at a given computation point.
+
+                    Args:
+                        computation_points: The computation points as tuple or list of points
+                        parallel:           If true, the computation is done in parallel (default: true)
+
+                    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);
+
+    py::module_ utility = m.def_submodule("utility",
+                                          "This submodule contains useful utility functions like parsing meshes "
+                                          "or checking if the polyhedron's mesh plane unit normals point outwards "
+                                          "like it is required by the polyhedral-gravity model.");
+
+    utility.def("read", [](const std::vector<std::string> &filenames) {
+        TetgenAdapter tetgen{filenames};
+        auto polyhedron = tetgen.getPolyhedron();
+        return std::make_tuple(polyhedron.getVertices(), polyhedron.getFaces());
+    }, R"mydelimiter(
+            Reads a polyhedron from a mesh file. The vertices and faces are read from input
+            files (either .node/.face, mesh, .ply, .off, .stl). File-Order matters in case of the first option!
+
+            Args:
+                input_files (List[str]): list of input files.
+            Returns:
+                tuple of vertices (N, 3) (floats) and faces (N, 3) (ints).
+          )mydelimiter", py::arg("input_files"));
+
+    utility.def("check_mesh", [](const std::vector<Array3> &vertices, const std::vector<IndexArray3> &faces) {
+        Polyhedron poly{vertices, faces};
+        return MeshChecking::checkTrianglesNotDegenerated(poly) && MeshChecking::checkNormalsOutwardPointing(poly);
+    }, R"mydelimiter(
+                Checks if no triangles of the polyhedral mesh are degenerated by checking that their surface area
+                is greater zero.
+                Further, Checks if all the polyhedron's plane unit normals are pointing outwards.
+
+                Args:
+                    vertices (2-D array-like): (N, 3) array of vertex coordinates (floats).
+                    faces (2-D array-like): (N, 3) array of faces, vertex-indices (ints).
+                Returns:
+                    True if no triangle is degenerate and the polyhedron's plane unit normals are all pointing outwards.
+          )mydelimiter", py::arg("vertices"), py::arg("faces"));
+
+
+    utility.def("check_mesh", [](const std::vector<std::string> &filenames) {
+        TetgenAdapter tetgen{filenames};
+        Polyhedron poly = tetgen.getPolyhedron();
+        return MeshChecking::checkTrianglesNotDegenerated(poly) && MeshChecking::checkNormalsOutwardPointing(poly);
+    }, R"mydelimiter(
+                Checks if no triangles of the polyhedral mesh are degenerated by checking that their surface area
+                is greater zero.
+                Further, Checks if all the polyhedron's plane unit normals are pointing outwards.
+                Reads a polyhedron from a mesh file. The vertices and faces are read from input
+                files (either .node/.face, mesh, .ply, .off, .stl). File-Order matters in case of the first option!
+
+                Args:
+                    input_files (List[str]): list of input files.
+                Returns:
+                    True if no triangle is degenerate and the polyhedron's plane unit normals are all pointing outwards.
+          )mydelimiter", py::arg("input_files"));
+}


Test patch:
diff --git a/.github/workflows/ctest.yml b/.github/workflows/ctest.yml
index 5825dd6..0913fd3 100644
--- a/.github/workflows/ctest.yml
+++ b/.github/workflows/ctest.yml
@@ -44,8 +44,11 @@ jobs:
       # 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 ..
diff --git a/cmake/gtest.cmake b/cmake/gtest.cmake
index da09655..c968cce 100644
--- a/cmake/gtest.cmake
+++ b/cmake/gtest.cmake
@@ -3,10 +3,10 @@ include(FetchContent)
 message(STATUS "Setting up gtest")
 
 #Adapted from https://cliutils.gitlab.io/modern-cmake/chapters/testing/googletest.html
-#Fetches the version 1.11.0 from the official github for googletest
+#Fetches the version 1.13.0 from the official github for googletest
 FetchContent_Declare(googletest
         GIT_REPOSITORY https://github.com/google/googletest.git
-        GIT_TAG release-1.11.0
+        GIT_TAG v1.13.0
         )
 
 FetchContent_MakeAvailable(googletest)
diff --git a/test/calculation/GravityModelCubeTest.cpp b/test/calculation/GravityModelCubeTest.cpp
index d6e1dd4..81ae505 100644
--- a/test/calculation/GravityModelCubeTest.cpp
+++ b/test/calculation/GravityModelCubeTest.cpp
@@ -83,10 +83,10 @@ TEST_P(GravityModelCubeTest, CubePoints) {
     double expectedPotential = std::get<1>(testData);
     std::array<double, 3> expectedAcceleration = std::get<2>(testData);
 
-    GravityModelResult actualResult = GravityModel::evaluate(_cube, _cube_density, computationPoint);
+    const auto[actualPotential, actualAcceleration, x] = GravityModel::evaluate(_cube, _cube_density, computationPoint);
 
-    ASSERT_NEAR(actualResult.gravitationalPotential, expectedPotential, LOCAL_TEST_EPSILON);
-    ASSERT_THAT(actualResult.acceleration, Pointwise(DoubleNear(LOCAL_TEST_EPSILON), expectedAcceleration));
+    ASSERT_NEAR(actualPotential, expectedPotential, LOCAL_TEST_EPSILON);
+    ASSERT_THAT(actualAcceleration, Pointwise(DoubleNear(LOCAL_TEST_EPSILON), expectedAcceleration));
 }
 
 INSTANTIATE_TEST_SUITE_P(CubeGravityModelTest, GravityModelCubeTest,