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

Issue:
Check validity of the polyhedra before computing gravity fields?
While reading some of the in #29, I noticed that the `evaluate()` function is not checking if the passed polyhedra is valid or not, and returning the gravity field values without any warning that they might not be right.

For example, grabbing the cube I used as example in #28, the following code runs with a valid cube and returns these values:

```python
import numpy as np
import polyhedral_gravity

cube_vertices = np.array(
    [
        [1, 1, 1],
        [2, 1, 1],
        [2, 2, 1],
        [1, 2, 1],
        [1, 1, 2],
        [2, 1, 2],
        [2, 2, 2],
        [1, 2, 2],
    ]
)
cube_faces = np.array(
    [
        [1, 3, 2],
        [0, 3, 1],
        [0, 1, 5],
        [0, 5, 4],
        [0, 7, 3],
        [0, 4, 7],
        [1, 2, 6],
        [1, 6, 5],
        [2, 3, 6],
        [3, 7, 6],
        [4, 5, 6],
        [4, 6, 7],
    ]
)
cube_density = 1.0
computation_point = np.array([0, 0, 0])

print("Valid cube?", polyhedral_gravity.utility.check_mesh(cube_vertices, cube_faces))

potential, acceleration, tensor = polyhedral_gravity.evaluate(
    polyhedral_source=(cube_vertices, cube_faces),
    density=cube_density,
    computation_points=computation_point,
    parallel=True,
)

print(f"potential: {potential}")
print(f"acceleration: {acceleration}")
print(f"tensor: {tensor}")
```

```
Valid cube? True
potential: 2.5695100395128033e-11
acceleration: [-5.715199089414289e-12, -5.715199089414289e-12, -5.7151990894142846e-12]
tensor: [2.2229884599767045e-26, 2.2229884599767045e-26, 2.2229884599767045e-26, 3.8189716252752374e-12, 3.818971625275311e-12, 3.8189716252752374e-12]
```

But, if I flip the order of the vertices indices in the first row of `cube_faces`, then the results look like these:

```python
import numpy as np
import polyhedral_gravity

cube_vertices = np.array(
    [
        [1, 1, 1],
        [2, 1, 1],
        [2, 2, 1],
        [1, 2, 1],
        [1, 1, 2],
        [2, 1, 2],
        [2, 2, 2],
        [1, 2, 2],
    ]
)
cube_faces = np.array(
    [
        [1, 2, 3],
        [0, 3, 1],
        [0, 1, 5],
        [0, 5, 4],
        [0, 7, 3],
        [0, 4, 7],
        [1, 2, 6],
        [1, 6, 5],
        [2, 3, 6],
        [3, 7, 6],
        [4, 5, 6],
        [4, 6, 7],
    ]
)
cube_density = 1.0
computation_point = np.array([0, 0, 0])

print("Valid cube?", polyhedral_gravity.utility.check_mesh(cube_vertices, cube_faces))

potential, acceleration, tensor = polyhedral_gravity.evaluate(
    polyhedral_source=(cube_vertices, cube_faces),
    density=cube_density,
    computation_points=computation_point,
    parallel=True,
)

print(f"potential: {potential}")
print(f"acceleration: {acceleration}")
print(f"tensor: {tensor}")
```

```
Valid cube? False
potential: 3.8680139904499555e-11
acceleration: [-5.715199089414289e-12, -5.715199089414289e-12, 2.025487992932876e-11]
tensor: [2.2229884599767045e-26, 2.2229884599767045e-26, -3.96702361734801e-12, 3.8189716252752374e-12, 3.818971625275311e-12, 3.8189716252752374e-12]
```

As you can see, the values of $g_z$ and $g_{zz}$ are significantly different between the two outputs.

I think it would be a good idea to include a validity check in the `evaluate()` function before computing the gravity fields. This would prevent users getting wrong results and therefore making wrong decisions or arriving at wrong conclusions because their results were wrong.

---

This issue is unrelated to the JOSS review. Addressing it won't be required for the JOSS paper to be published, but I think making the `evaluate()` function less error prone would be a major feature for this package.
Check validity of the polyhedra before computing gravity fields?
While reading some of the in #29, I noticed that the `evaluate()` function is not checking if the passed polyhedra is valid or not, and returning the gravity field values without any warning that they might not be right.

For example, grabbing the cube I used as example in #28, the following code runs with a valid cube and returns these values:

```python
import numpy as np
import polyhedral_gravity

cube_vertices = np.array(
    [
        [1, 1, 1],
        [2, 1, 1],
        [2, 2, 1],
        [1, 2, 1],
        [1, 1, 2],
        [2, 1, 2],
        [2, 2, 2],
        [1, 2, 2],
    ]
)
cube_faces = np.array(
    [
        [1, 3, 2],
        [0, 3, 1],
        [0, 1, 5],
        [0, 5, 4],
        [0, 7, 3],
        [0, 4, 7],
        [1, 2, 6],
        [1, 6, 5],
        [2, 3, 6],
        [3, 7, 6],
        [4, 5, 6],
        [4, 6, 7],
    ]
)
cube_density = 1.0
computation_point = np.array([0, 0, 0])

print("Valid cube?", polyhedral_gravity.utility.check_mesh(cube_vertices, cube_faces))

potential, acceleration, tensor = polyhedral_gravity.evaluate(
    polyhedral_source=(cube_vertices, cube_faces),
    density=cube_density,
    computation_points=computation_point,
    parallel=True,
)

print(f"potential: {potential}")
print(f"acceleration: {acceleration}")
print(f"tensor: {tensor}")
```

```
Valid cube? True
potential: 2.5695100395128033e-11
acceleration: [-5.715199089414289e-12, -5.715199089414289e-12, -5.7151990894142846e-12]
tensor: [2.2229884599767045e-26, 2.2229884599767045e-26, 2.2229884599767045e-26, 3.8189716252752374e-12, 3.818971625275311e-12, 3.8189716252752374e-12]
```

But, if I flip the order of the vertices indices in the first row of `cube_faces`, then the results look like these:

```python
import numpy as np
import polyhedral_gravity

cube_vertices = np.array(
    [
        [1, 1, 1],
        [2, 1, 1],
        [2, 2, 1],
        [1, 2, 1],
        [1, 1, 2],
        [2, 1, 2],
        [2, 2, 2],
        [1, 2, 2],
    ]
)
cube_faces = np.array(
    [
        [1, 2, 3],
        [0, 3, 1],
        [0, 1, 5],
        [0, 5, 4],
        [0, 7, 3],
        [0, 4, 7],
        [1, 2, 6],
        [1, 6, 5],
        [2, 3, 6],
        [3, 7, 6],
        [4, 5, 6],
        [4, 6, 7],
    ]
)
cube_density = 1.0
computation_point = np.array([0, 0, 0])

print("Valid cube?", polyhedral_gravity.utility.check_mesh(cube_vertices, cube_faces))

potential, acceleration, tensor = polyhedral_gravity.evaluate(
    polyhedral_source=(cube_vertices, cube_faces),
    density=cube_density,
    computation_points=computation_point,
    parallel=True,
)

print(f"potential: {potential}")
print(f"acceleration: {acceleration}")
print(f"tensor: {tensor}")
```

```
Valid cube? False
potential: 3.8680139904499555e-11
acceleration: [-5.715199089414289e-12, -5.715199089414289e-12, 2.025487992932876e-11]
tensor: [2.2229884599767045e-26, 2.2229884599767045e-26, -3.96702361734801e-12, 3.8189716252752374e-12, 3.818971625275311e-12, 3.8189716252752374e-12]
```

As you can see, the values of $g_z$ and $g_{zz}$ are significantly different between the two outputs.

I think it would be a good idea to include a validity check in the `evaluate()` function before computing the gravity fields. This would prevent users getting wrong results and therefore making wrong decisions or arriving at wrong conclusions because their results were wrong.

---

This issue is unrelated to the JOSS review. Addressing it won't be required for the JOSS paper to be published, but I think making the `evaluate()` function less error prone would be a major feature for this package.


Discussion / hints:
Good point.

However, this check is intentionally not included in the `evaluate(..)` function due to its $O(n^2)$ runtime.
As background, we're using the Möller–Trumbore intersection algorithm. This algorithm computes whether a ray intersects a plane. The number of intersections must be even to check that a plane unit normally points outwards. This is done for every plane normal for every other plane.
The approach works for both convex and concave polyhedrons. The disadvantage is that it is quadratic, at least in the naiv implementation of the `utility` package.

An Improvement of the `utility` would be a future TODO. But even then, $O(n)$, as it is the runtime of the gravity model, is not reachable (to my knowledge).

--------
What I can do, however, is to improve the visibility of this check in the package.

Hence, I've stated the `mesh_checking` capabilities more clearly in the `README.md` to draw the attention of the user to this requirement. I've also slightly updated the documentation in this regard.

Resolved by #31 

[1] https://stackoverflow.com/questions/45603469/how-to-calculate-the-normals-of-a-box

I understand your point. I agree with making the `utility` module more discoverable through documentation.

What if the mesh check is added to the forward modelling function but also provide an optional argument to disable it. Since I would recommend every user to run the check before computing the gravity field, I would turn it on by default. In some particular situations users might not want to run the check , so they could do so by setting the argument to `False`.

Something like:
```python
potential, acceleration, tensor = polyhedral_gravity.evaluate(
    polyhedral_source=(cube_vertices, cube_faces),
    density=cube_density,
    computation_points=computation_point,
    parallel=True,
    disable_check=True,
)
```

What do you think?

Agreed. I would even go a step further and enable the direct utilization with a "correct" mesh but reverse vertex ordering.
Further, I would add an interface to get a proper hint of where the vertex ordering is not aligned.

The algorithm requires the plane unit normals to be outwards pointing.
Hence, counterclockwise vertex ordering in a right-handed coordinate system.
I am not completely sure since this definition (clockwise/counterclockwise) also depends on the coordinate system.

Hence, I would just specify an enum "OUTWARDS" or "INWARDS" for the following feature.

# Changes

## Core Interface

- [ ] Add a parameter `polyhedral_orientation` to `evaluate(..)` --> `Enum Class NORMAL_ORIENTATION`
  - We implicitly assume always `OUTWARDS` ordering if not differently specified (Default Case)
  - The results are multiplied with `-1.0` if `INWARDS ` is specified
- [ ] Add a parameter `check_mesh` to `evaluate(..)` --> `optional<bool>`
  - The parameter is, by default, implicitly enabled and prints a warning that the mesh check is conducted requiring $O(n^2)$ (Default Case)
  - Explicitly setting the parameter to `True ` or `False ` silences the warning
  - In case the parameter is `True` a `runtime_error` is thrown if the mesh does not comply with the specified `polyhedral_orientation`
  - The method also conducts the $O(n)$ check for triangles without surface area
- [ ] Add the same behavior to `GravityEvaluable` in the constructor

## Utility Methods

- [ ] The `mesh_check` method is removed due to its integration into the core interface
- [ ] A method `check_vertex_order` is introduced, which returns
  - **Majority Ordering**
    - INWARDS --> All plane unit normals intersect the polyhedron odd times
    - OUTWARDS --> All plane unit normals intersect the polyhedron even times
  - **Indices of faces violating the order**
    - a list of face-indices violating the order/ might be empty if everything is "fine"
  The implementation determines which ordering is "more reasonable." If the majority of faces have inward-pointing normals, INWARD_POINTING is returned with "faulty" outward-pointing face indices, and vice versa.
  This way, I get the correct arguments for `polyhedral_orientation` and can abort if the list is non-empty (or exchange the vertex ordering at the returned faces)
I will add these features in the next days, but in a new PR other than #31 
I heavily support opening a new PR for this.

Just a few minor thoughts:

> Explicitly setting the parameter to True or False silences the warning

I think setting it to `False` should prevent the check to being run, not only silencing the warning. This way users can choose to save computation time by disabling the checks (maybe that's what you meant, but I wanted to clarify it).

> The mesh_check method is removed due to its integration into the core interface

I wouldn't remove it. You already have it, and it's something useful to have. Specially while using the `evaluate` method. Maybe some user wants to manually make sure the mesh is correct before running the forward.

> Add the same behavior to GravityEvaluable in the constructor

That would be great! And I think you can only run the check once, no need to run it every time the user ask for the computation of the gravity fields. Maybe the check could be run on instantiation of the class. Or maybe check before the first time the gravity fields are computed, and store a private attribute that flags that the mesh has been checked and is valid.

I really like your approach here! Thanks a lot for considering my idea.
Good point.

However, this check is intentionally not included in the `evaluate(..)` function due to its $O(n^2)$ runtime.
As background, we're using the Möller–Trumbore intersection algorithm. This algorithm computes whether a ray intersects a plane. The number of intersections must be even to check that a plane unit normally points outwards. This is done for every plane normal for every other plane.
The approach works for both convex and concave polyhedrons. The disadvantage is that it is quadratic, at least in the naiv implementation of the `utility` package.

An Improvement of the `utility` would be a future TODO. But even then, $O(n)$, as it is the runtime of the gravity model, is not reachable (to my knowledge).

--------
What I can do, however, is to improve the visibility of this check in the package.

Hence, I've stated the `mesh_checking` capabilities more clearly in the `README.md` to draw the attention of the user to this requirement. I've also slightly updated the documentation in this regard.

Resolved by #31 

[1] https://stackoverflow.com/questions/45603469/how-to-calculate-the-normals-of-a-box

I understand your point. I agree with making the `utility` module more discoverable through documentation.

What if the mesh check is added to the forward modelling function but also provide an optional argument to disable it. Since I would recommend every user to run the check before computing the gravity field, I would turn it on by default. In some particular situations users might not want to run the check , so they could do so by setting the argument to `False`.

Something like:
```python
potential, acceleration, tensor = polyhedral_gravity.evaluate(
    polyhedral_source=(cube_vertices, cube_faces),
    density=cube_density,
    computation_points=computation_point,
    parallel=True,
    disable_check=True,
)
```

What do you think?

Agreed. I would even go a step further and enable the direct utilization with a "correct" mesh but reverse vertex ordering.
Further, I would add an interface to get a proper hint of where the vertex ordering is not aligned.

The algorithm requires the plane unit normals to be outwards pointing.
Hence, counterclockwise vertex ordering in a right-handed coordinate system.
I am not completely sure since this definition (clockwise/counterclockwise) also depends on the coordinate system.

Hence, I would just specify an enum "OUTWARDS" or "INWARDS" for the following feature.

# Changes

## Core Interface

- [ ] Add a parameter `polyhedral_orientation` to `evaluate(..)` --> `Enum Class NORMAL_ORIENTATION`
  - We implicitly assume always `OUTWARDS` ordering if not differently specified (Default Case)
  - The results are multiplied with `-1.0` if `INWARDS ` is specified
- [ ] Add a parameter `check_mesh` to `evaluate(..)` --> `optional<bool>`
  - The parameter is, by default, implicitly enabled and prints a warning that the mesh check is conducted requiring $O(n^2)$ (Default Case)
  - Explicitly setting the parameter to `True ` or `False ` silences the warning
  - In case the parameter is `True` a `runtime_error` is thrown if the mesh does not comply with the specified `polyhedral_orientation`
  - The method also conducts the $O(n)$ check for triangles without surface area
- [ ] Add the same behavior to `GravityEvaluable` in the constructor

## Utility Methods

- [ ] The `mesh_check` method is removed due to its integration into the core interface
- [ ] A method `check_vertex_order` is introduced, which returns
  - **Majority Ordering**
    - INWARDS --> All plane unit normals intersect the polyhedron odd times
    - OUTWARDS --> All plane unit normals intersect the polyhedron even times
  - **Indices of faces violating the order**
    - a list of face-indices violating the order/ might be empty if everything is "fine"
  The implementation determines which ordering is "more reasonable." If the majority of faces have inward-pointing normals, INWARD_POINTING is returned with "faulty" outward-pointing face indices, and vice versa.
  This way, I get the correct arguments for `polyhedral_orientation` and can abort if the list is non-empty (or exchange the vertex ordering at the returned faces)
I will add these features in the next days, but in a new PR other than #31 
I heavily support opening a new PR for this.

Just a few minor thoughts:

> Explicitly setting the parameter to True or False silences the warning

I think setting it to `False` should prevent the check to being run, not only silencing the warning. This way users can choose to save computation time by disabling the checks (maybe that's what you meant, but I wanted to clarify it).

> The mesh_check method is removed due to its integration into the core interface

I wouldn't remove it. You already have it, and it's something useful to have. Specially while using the `evaluate` method. Maybe some user wants to manually make sure the mesh is correct before running the forward.

> Add the same behavior to GravityEvaluable in the constructor

That would be great! And I think you can only run the check once, no need to run it every time the user ask for the computation of the gravity fields. Maybe the check could be run on instantiation of the class. Or maybe check before the first time the gravity fields are computed, and store a private attribute that flags that the mesh has been checked and is valid.

I really like your approach here! Thanks a lot for considering my idea.

Fix patch:
diff --git a/.clang-format b/.clang-format
new file mode 100644
index 0000000..4ffb425
--- /dev/null
+++ b/.clang-format
@@ -0,0 +1,66 @@
+# Generated from CLion C/C++ Code Style settings
+BasedOnStyle: LLVM
+AccessModifierOffset: -4
+AlignAfterOpenBracket: Align
+AlignConsecutiveAssignments: None
+AlignOperands: Align
+AllowAllArgumentsOnNextLine: false
+AllowAllConstructorInitializersOnNextLine: false
+AllowAllParametersOfDeclarationOnNextLine: false
+AllowShortBlocksOnASingleLine: Always
+AllowShortCaseLabelsOnASingleLine: false
+AllowShortFunctionsOnASingleLine: None
+AllowShortIfStatementsOnASingleLine: Always
+AllowShortLambdasOnASingleLine: All
+AllowShortLoopsOnASingleLine: true
+AlwaysBreakAfterReturnType: None
+AlwaysBreakTemplateDeclarations: Yes
+BreakBeforeBraces: Custom
+BraceWrapping:
+  AfterCaseLabel: false
+  AfterClass: false
+  AfterControlStatement: Never
+  AfterEnum: false
+  AfterFunction: false
+  AfterNamespace: false
+  AfterUnion: false
+  BeforeCatch: false
+  BeforeElse: false
+  IndentBraces: false
+  SplitEmptyFunction: false
+  SplitEmptyRecord: true
+BreakBeforeBinaryOperators: None
+BreakBeforeTernaryOperators: true
+BreakConstructorInitializers: BeforeColon
+BreakInheritanceList: BeforeColon
+ColumnLimit: 0
+CompactNamespaces: false
+ContinuationIndentWidth: 8
+IndentCaseLabels: true
+IndentPPDirectives: None
+IndentWidth: 4
+KeepEmptyLinesAtTheStartOfBlocks: true
+MaxEmptyLinesToKeep: 2
+NamespaceIndentation: All
+ObjCSpaceAfterProperty: false
+ObjCSpaceBeforeProtocolList: true
+PointerAlignment: Right
+ReflowComments: false
+SpaceAfterCStyleCast: true
+SpaceAfterLogicalNot: false
+SpaceAfterTemplateKeyword: false
+SpaceBeforeAssignmentOperators: true
+SpaceBeforeCpp11BracedList: false
+SpaceBeforeCtorInitializerColon: true
+SpaceBeforeInheritanceColon: true
+SpaceBeforeParens: ControlStatements
+SpaceBeforeRangeBasedForLoopColon: false
+SpaceInEmptyParentheses: false
+SpacesBeforeTrailingComments: 0
+SpacesInAngles: false
+SpacesInCStyleCastParentheses: false
+SpacesInContainerLiterals: false
+SpacesInParentheses: false
+SpacesInSquareBrackets: false
+TabWidth: 4
+UseTab: Never
diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml
index 3a21df2..9eff77a 100644
--- a/.github/workflows/docs.yml
+++ b/.github/workflows/docs.yml
@@ -25,6 +25,9 @@ jobs:
       - uses: actions/checkout@v3
         with:
           lfs: true
+      - name: Install polyhedral-gravity from recent source for sphinx.ext.autodoc
+        shell: bash -l {0}
+        run: pip install . -vv --no-build-isolation
       - id: deployment
         uses: sphinx-notes/pages@v3
         with:
diff --git a/.github/workflows/wheels.yml b/.github/workflows/wheels.yml
index 39f5b87..eef7ba1 100644
--- a/.github/workflows/wheels.yml
+++ b/.github/workflows/wheels.yml
@@ -21,10 +21,11 @@ jobs:
         os: [windows-latest, ubuntu-latest, macos-latest]
     steps:
       - uses: actions/checkout@v3
+        ############################# LINUX WHEELS #############################
         # In case of Linux we need to install compiler and build tools before building the wheels
         # We further only build the manylinux wheels, but not the musllinux wheels
       - name: Build wheels (Linux)
-        uses: pypa/cibuildwheel@v2.15.0
+        uses: pypa/cibuildwheel@v2.17.0
         env:
           CIBW_BEFORE_BUILD: yum makecache && yum install -y gcc-c++ cmake && pip install ninja
           CIBW_BUILD: "*manylinux*"
@@ -34,23 +35,36 @@ jobs:
           package-dir: .
           output-dir: dist
         if: matrix.os == 'ubuntu-latest'
+      ############################# MACOS WHEELS #############################
       # Building on macOS requires an installation of gcc since the default clang compiler
       # lacks certain features required for building the package
-      - name: Build wheels (macOS)
-        uses: pypa/cibuildwheel@v2.15.0
+      - name: Build wheels (macOS ARM)
+        uses: pypa/cibuildwheel@v2.17.0
+        env:
+          CIBW_BEFORE_BUILD: brew install ninja
+          CIBW_ARCHS_MACOS: "arm64"
+          CIBW_TEST_COMMAND: 'python -c "import polyhedral_gravity"'
+        with:
+          package-dir: .
+          output-dir: dist
+        if: matrix.os == 'macos-latest'
+      - name: Build wheels (macOS x86_64)
+        uses: pypa/cibuildwheel@v2.17.0
         env:
           CIBW_BEFORE_BUILD: brew install ninja gcc@12
           CIBW_ENVIRONMENT: "CC=gcc-12 CXX=g++-12"
+          CIBW_ARCHS_MACOS: "x86_64"
           CIBW_TEST_COMMAND: 'python -c "import polyhedral_gravity"'
         with:
           package-dir: .
           output-dir: dist
         if: matrix.os == 'macos-latest'
+      ############################# WINDOWS WHEELS #############################
       # Set up the Visual Studio environment on Windows (required, so that CMake finds the compiler)
       - uses: ilammy/msvc-dev-cmd@v1
         if: matrix.os == 'windows-latest'
       - name: Build wheels (Windows)
-        uses: pypa/cibuildwheel@v2.15.0
+        uses: pypa/cibuildwheel@v2.17.0
         env:
           CIBW_BEFORE_BUILD: choco install -y ninja cmake
           CIBW_ARCHS_WINDOWS: "auto64"
@@ -94,44 +108,10 @@ jobs:
         with:
           repository-url: https://test.pypi.org/legacy/
 
-  # 3. Check if the package can be installed from testpypi
-  #    Notice, that this is more of an installation test since the
-  #    import check has already been done in the build_wheels job
-  check_testpypi:
-    needs: [upload_testpypi]
-    name: Test import on ${{ matrix.os }} with ${{ matrix.py }}
-    runs-on: ${{ matrix.os }}
-    if: always()
-    strategy:
-      fail-fast: false
-      matrix:
-        os: [windows-latest, ubuntu-latest, macos-latest]
-        py:
-          [
-            3.7.x,
-            3.8.x,
-            3.9.x,
-            3.10.x,
-            3.11.x,
-            pypy3.7,
-            pypy3.8,
-            pypy3.9,
-            pypy3.10,
-          ]
-    steps:
-      - uses: actions/checkout@v3
-      - uses: actions/setup-python@v4
-        with:
-          python-version: ${{ matrix.py }}
-      - name: Install package from testpypi
-        run: pip install --index-url https://test.pypi.org/simple/ polyhedral-gravity
-      - name: Check import
-        run: python -c "import polyhedral_gravity"
-
-  # 4. Upload the wheels to the actually Python Package Index
+  # 3. Upload the wheels to the actually Python Package Index
   #    using trusted publishing
   upload_pypi:
-    needs: [build_wheels, make_sdist, check_testpypi]
+    needs: [build_wheels, make_sdist, upload_testpypi]
     environment:
       name: pypi
       url: https://pypi.org/p/polyhedral-gravity
diff --git a/.gitignore b/.gitignore
index 3fe54c9..faa3a36 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,4 +2,5 @@
 cmake-build-*
 build
 polyhedral_gravity.egg-info
-dist
\ No newline at end of file
+dist
+docs/Doxyfile
\ No newline at end of file
diff --git a/README.md b/README.md
index ed40b38..22805c9 100644
--- a/README.md
+++ b/README.md
@@ -60,8 +60,8 @@ which is strongly based on the former implementation in FORTRAN.
 
 > [!NOTE]
 > The [GitHub Pages](https://esa.github.io/polyhedral-gravity-model) of this project
-contain the full extensive documentation.
-It also covers the content of the `polyhedral_gravity.utility` module to check the mesh sanity.
+contain the full extensive documentation of the C++ Library and Python Interface
+as well as background on the gravity model and advanced settings not detailed here.
 
 ## Input & Output (C++ and Python)
 
@@ -75,15 +75,7 @@ The evaluation of the polyhedral gravity model requires the following parameters
 | Constant Density $\rho$                                                    |
 
 The mesh and the constants density's unit must match.
-Have a look the documentation to view the [supported mesh files](https://esa.github.io/polyhedral-gravity-model/supported_input.html).
-
-> [!IMPORTANT]  
-> The plane unit normals of every face of the polyhedral mesh must point **outwards**
-of the polyhedron!
-You can check this property via [MeshChecking](https://esa.github.io/polyhedral-gravity-model/api/calculation.html#meshchecking) in C++ or
-via the [utility](https://esa.github.io/polyhedral-gravity-model/api/python.html#module-polyhedral_gravity.utility) submodule in Python.
-If the vertex order of the faces is inverted, i.e. the plane unit normals point
-inwards, then the sign of the output will be inverted.
+Have a look the documentation to view the [supported mesh files](https://esa.github.io/polyhedral-gravity-model/quickstart/supported_input.html).
 
 ### Output
 
@@ -111,31 +103,43 @@ around a cube:
 
 ```python
 import numpy as np
-import polyhedral_gravity
+from polyhedral_gravity import Polyhedron, GravityEvaluable, evaluate, PolyhedronIntegrity, NormalOrientation
 
 # We define the cube as a polyhedron with 8 vertices and 12 triangular faces
-# The polyhedron's (here: cube) normals need to point outwards (see below for checking this)
+# The polyhedron's normals point outwards (see below for checking this)
 # The density is set to 1.0
 cube_vertices = np.array(
-    [[-1, -1, -1], [1, -1, -1], [1, 1, -1], [-1, 1, -1],
-    [-1, -1, 1], [1, -1, 1], [1, 1, 1], [-1, 1, 1]]
+  [[-1, -1, -1], [1, -1, -1], [1, 1, -1], [-1, 1, -1],
+   [-1, -1, 1], [1, -1, 1], [1, 1, 1], [-1, 1, 1]]
 )
 cube_faces = np.array(
-    [[1, 3, 2], [0, 3, 1], [0, 1, 5], [0, 5, 4], [0, 7, 3], [0, 4, 7],
-    [1, 2, 6], [1, 6, 5], [2, 3, 6], [3, 7, 6], [4, 5, 6], [4, 6, 7]]
+  [[1, 3, 2], [0, 3, 1], [0, 1, 5], [0, 5, 4], [0, 7, 3], [0, 4, 7],
+   [1, 2, 6], [1, 6, 5], [2, 3, 6], [3, 7, 6], [4, 5, 6], [4, 6, 7]]
 )
 cube_density = 1.0
 computation_point = np.array([0, 0, 0])
 ```
 
-The simplest way to compute the gravity is to use the `evaluate` function:
+We first define a constant density Polyhedron from `vertices` and `faces`
 
 ```python
-potential, acceleration, tensor = polyhedral_gravity.evaluate(
-    polyhedral_source=(cube_vertices, cube_faces),
-    density=cube_density,
-    computation_points=computation_point,
-    parallel=True
+cube_polyhedron = Polyhedron(
+  polyhedral_source=(cube_vertices, cube_faces),
+  density=cube_density,
+)
+```
+
+In case you want to hand over the polyhedron via a [supported file format](https://esa.github.io/polyhedral-gravity-model/quickstart/supported_input.html),
+just replace the `polyhedral_source` argument with *a list of strings*,
+where each string is the path to a supported file format, e.g. `polyhedral_source=["eros.node","eros.face"]` or `polyhedral_source=["eros.mesh"]`.
+
+Continuing, the simplest way to compute the gravity is to use the `evaluate` function:
+
+```python
+potential, acceleration, tensor = evaluate(
+  polyhedron=cube_polyhedron,
+  computation_points=computation_point,
+  parallel=True,
 )
 ```
 
@@ -145,31 +149,36 @@ evaluations. This is especially useful if you want to compute the gravity
 for multiple computation points, but don't know the "future points" in advance.
 
 ```python
-evaluable = polyhedral_gravity.GravityEvaluable(
-    polyhedral_source=(cube_vertices, cube_faces),
-    density=cube_density
+evaluable = GravityEvaluable(polyhedron=cube_polyhedron) # stores intermediate computation steps
+potential, acceleration, tensor = evaluable(
+  computation_points=computation_point,
+  parallel=True,
 )
-potential, acceleration, tensor = evaluable(computation_point, parallel=True)
+# Any future evaluable call after this one will be faster
 ```
 
-In case you want to hand over the polyhedron via a supported file format,
-just replace the `polyhedral_source` argument with *a list of strings*,
-where each string is the path to a supported file format.
-Note that the `computation_point` could also be (N, 3)-shaped array.
+Note that the `computation_point` could also be (N, 3)-shaped array to compute multiple points at once.
 In this case, the return value of `evaluate(..)` or an `GravityEvaluable` will
 be a list of triplets comprising potential, acceleration, and tensor.
 
-The gravity model requires that the polyhedron's plane unit normals point outwards
-the polyhedron. In case you are unsure, you can check for this property by using the `utility` module beforehand.
-The method also verifies that all triangles are actually triangles with a non-zero
-surface area.
+The gravity model requires that all the polyhedron's plane unit normals consistently
+point outwards or inwards the polyhedron. You can specify this via the `normal_orientation`.
+This property is - by default - checked when constructing the `Polyhedron`! So, don't worry, it
+is impossible if not **explicitly** disabled to create an invalid `Polyhedron`.
+You can disable/ enable this setting via the optional `integrity_check` flag and can even
+automatically repair the ordering via `HEAL`.
+If you are confident that your mesh is defined correctly (e.g. checked once with the integrity check)
+you can disable this check (via `DISABLE`) to avoid the additional runtime overhead of the check.
 
 ```python
-print("Valid Mesh?", polyhedral_gravity.utility.check_mesh(cube_vertices, cube_faces))
+cube_polyhedron = Polyhedron(
+  polyhedral_source=(cube_vertices, cube_faces),
+  density=cube_density,
+  normal_orientation=NormalOrientation.INWARDS, # OUTWARDS (default) or INWARDS
+  integrity_check=PolyhedronIntegrity.VERIFY,   # VERIFY (default), DISABLE or HEAL
+)
 ```
 
-If the method returns `False`, then you need to revise the vertex-ordering.
-
 > [!TIP]
 > More examples and plots are depicted in the
 [jupyter notebook](script/polyhedral-gravity.ipynb).
@@ -184,9 +193,11 @@ It works analogously to the Python example above.
 // Defining the input like above in the Python example
 std::vector<std::array<double, 3>> vertices = ...
 std::vector<std::array<size_t, 3>> faces = ...
-// The polyhedron is defined by its vertices and faces
-Polyhedron polyhedron{vertices, faces};
 double density = 1.0;
+// The constant density polyhedron is defined by its vertices & faces
+// It also supports the hand-over of NormalOrientation and PolyhedronIntegrity as optional arguments
+// as above described for the Python Interface
+Polyhedron polyhedron{vertices, faces, density};
 std::vector<std::array<double, 3>> points = ...
 std::array<double, 3> point = points[0];
 bool parallel = true;
@@ -196,19 +207,19 @@ The C++ library provides also two ways to compute the gravity. Via
 the free function `evaluate`...
 
 ```cpp
-const auto[pot, acc, tensor] = GravityModel::evaluate(polyhedron, density, point);
+const auto[pot, acc, tensor] = GravityModel::evaluate(polyhedron, point, parallel);
 ```
 
 ... or via the `GravityEvaluable` class.
 
 ```cpp
 // Instantiation of the GravityEvaluable object
-GravityEvaluable evaluable{polyhedron, density};
+GravityEvaluable evaluable{polyhedron};
 
 // From now, we can evaluate the gravity model for any point with
 const auto[potential, acceleration, tensor] = evaluable(point, parallel);
 // or for multiple points with
-const auto results = evaluable(points);
+const auto results = evaluable(points, parallel);
 ```
 
 Similarly to Python, the C++ implementation also provides mesh checking capabilities.
@@ -304,8 +315,7 @@ The following options are available:
 During testing POLYHEDRAL_GRAVITY_PARALLELIZATION=`TBB` has been the most performant.
 It is further not recommend to change the LOGGING_LEVEL to something else than `INFO=2`.
 
-The recommended CMake command would look like this (we only need to change `PARALLELIZATION_DEVICE`, since
-the defaults of the others are already correctly set):
+The recommended CMake settings using the `TBB` backend would look like this:
 
 ```bash
 cmake .. -POLYHEDRAL_GRAVITY_PARALLELIZATION="TBB"
@@ -316,7 +326,7 @@ cmake .. -POLYHEDRAL_GRAVITY_PARALLELIZATION="TBB"
 After the build, the gravity model can be run by executing:
 
 ```bash
-    ./polyhedralGravity <YAML-Configuration-File>
+./polyhedralGravity <YAML-Configuration-File>
 ```
 
 where the YAML-Configuration-File contains the required parameters.
@@ -327,7 +337,7 @@ found in this repository in the folder `/example-config/`.
 
 The configuration should look similar to the given example below.
 It is required to specify the source-files of the polyhedron's mesh (more info
-about the supported file in the [previous paragraph](#supported-polyhedron-source-files-python-c)), the density
+about the supported file in the [documentation](https://esa.github.io/polyhedral-gravity-model/quickstart/supported_input.html)), the density
 of the polyhedron, and the wished computation points where the
 gravity tensor shall be computed.
 Further one must specify the name of the .csv output file.
@@ -342,7 +352,8 @@ gravityModel:
     density: 2670.0                             # constant density, units must match with the mesh (see section below)
     points: # Location of the computation point(s) P
       - [ 0, 0, 0 ]                             # Here it is situated at the origin
-    check_mesh: true                            # Fully optional, enables input checking (not given: false)
+    check_mesh: true                            # Fully optional, enables mesh autodetect+repair of 
+                                                # the polyhedron's vertex ordering (not given: true)
   output:
     filename: "gravity_result.csv"              # The name of the output file 
 
diff --git a/docs/api/calculation.rst b/docs/api/calculation.rst
deleted file mode 100644
index d436430..0000000
--- a/docs/api/calculation.rst
+++ /dev/null
@@ -1,42 +0,0 @@
-Calculation
-===========
-
-Overview
---------
-
-The calculation is the heart of the Polyhedral Gravity Model
-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
-------------
-
-.. doxygennamespace:: polyhedralGravity::GravityModel
-
-
-.. _mesh-checking-cpp:
-
-MeshChecking
-------------
-
-.. doxygennamespace:: polyhedralGravity::MeshChecking
-
-Other
------
-
-.. doxygenfunction:: polyhedralGravity::transformPolyhedron
diff --git a/docs/api/model.rst b/docs/api/model.rst
index 821adc9..2813775 100644
--- a/docs/api/model.rst
+++ b/docs/api/model.rst
@@ -4,18 +4,53 @@ Model
 Overview
 --------
 
-The Model contains two types of data containers. First, there
-is the class :code:`Polyhedron` which consists of a vector
-vertices and a vector of triangular faces. Second, there
-are Utility Container whose single purpose is improving
-code readability in providing named access to their
-members.
+The model is the heart of the Polyhedral Gravity Model
+since it contains the the two major classes:
+
+* :cpp:class:`polyhedralGravity::Polyhedron`
+* :cpp:class:`polyhedralGravity::GravityEvaluable`
+
+The class :cpp:class:`polyhedralGravity::Polyhedron` stores the mesh and density
+information and governs the compliance with Tsoulis et al.’s gravity model’s preconditions
+It ensures that all plane unit normals of a constructed Polyhedron are consistently
+pointing :cpp:enumerator:`polyhedralGravity::NormalOrientation::OUTWARDS` or
+:cpp:enumerator:`polyhedralGravity::NormalOrientation::INWARDS`.
+Depending on the value of :cpp:enum:`polyhedralGravity::PolyhedronIntegrity`, these
+constraints are either enforced by modifying input mesh data or by throwing
+an :cpp:class:`std::invalid_argument` exception.
+For this purpose, it uses the `Möller–Trumbore intersection algorithm <https://en.wikipedia.org/wiki/Möller–Trumbore_intersection_algorithm>`__.
+
+
+The class :cpp:class:`polyhedralGravity::GravityEvaluable` is provides 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
+caching is performed.
+It takes a :cpp:class:`polyhedralGravity::Polyhedron` and provides an
+:cpp:func:`polyhedralGravity::GravityEvaluable::operator()` to evaluate the
+model at computation point(s) :math:`P` optionally using parallelization.
+A :cpp:func:`polyhedralGravity::GravityModel::evaluate` summarizes the
+functionality :cpp:class:`polyhedralGravity::GravityEvaluable`, but does not
+provide any caching throughout multiple calls.
+
 
 Polyhedron
 ----------
 
 .. doxygenclass:: polyhedralGravity::Polyhedron
 
+.. doxygenenum:: polyhedralGravity::NormalOrientation
+
+.. doxygenenum:: polyhedralGravity::PolyhedronIntegrity
+
+
+GravityModel
+------------
+
+.. doxygenclass:: polyhedralGravity::GravityEvaluable
+
+.. doxygennamespace:: polyhedralGravity::GravityModel
+
+
 Named Tuple
 -----------
 
@@ -23,4 +58,4 @@ Named Tuple
 
 .. doxygenstruct:: polyhedralGravity::TranscendentalExpression
 
-.. doxygenstruct:: polyhedralGravity::HessianPlane
+.. doxygenstruct:: polyhedralGravity::HessianPlane
\ No newline at end of file
diff --git a/docs/api/python.rst b/docs/api/python.rst
index 640e439..043fd9f 100644
--- a/docs/api/python.rst
+++ b/docs/api/python.rst
@@ -1,125 +1,38 @@
-Python API
-==========
+Python API - polyhedral_gravity
+===============================
 
-polyhedral_gravity
-------------------
+Module Overview
+---------------
 
-Computes the full gravity tensor for a given constant density polyhedron which consists of some
-vertices and triangular faces at given computation points
+.. automodule:: polyhedral_gravity
 
-.. py:module:: polyhedral_gravity
+Polyhedron
+----------
 
-    .. py:function:: evaluate(polyhedral_source, density, computation_points, parallel=True)
+.. autoclass:: polyhedral_gravity.Polyhedron
+   :members:
+   :special-members: __init__, __getitem__, __repr__
 
-        Evaluates the polyhedral gravity model for a given constant density polyhedron at a given computation point.
+Enums to specify MeshChecks
+~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-        :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)
+.. autoclass:: polyhedral_gravity.NormalOrientation
 
-        :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
+.. autoclass:: polyhedral_gravity.PolyhedronIntegrity
 
 
-    .. py:class:: GravityEvaluable
+GravityModel
+------------
 
-        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.
+Single Function
+~~~~~~~~~~~~~~~
 
-        .. py:method:: __init__(polyhedral_source, density)
+.. autofunction:: polyhedral_gravity.evaluate
 
-            :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}`
 
-        .. py:method:: __call__(computation_points, parallel=True)
-
-            Evaluates the polyhedral gravity model for a given constant density polyhedron at a given computation point.
-
-            :param List[float[3]] computation_points:
-                The cartesian computation point
-            :param bool parallel:
-                If true, the computation is done in parallel (default: true)
-
-            :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.
-
-        .. py:method:: __repr__()
-
-            Returns a string representation of the polyhedral gravity evaluable.
-
-            :rtype: str
-            :return: A string representation of the polyhedron
-
-
-.. _mesh-checking-python:
-
-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.
-
-.. py:function:: read(input_files)
-    :noindex:
-
-    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!
-
-    :param List[str] input_files: polyhedral source files
-    :return: tuple of vertices (N, 3) (floats) and faces (N, 3) (ints)
-    :rtype: Tuple[List[List[float[3]]], List[List[int[3]]]]
-
-.. py:function:: check_mesh(vertices, faces)
-    :noindex:
-
-    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.
-
-    :param List[List[float[3]]] vertices: vertices of the polyhedron
-    :param List[List[int[3]]] faces: faces of the polyhedron
-    :return: True if no triangle is degenerate and the polyhedron's plane unit normals are all pointing outwards.
-    :rtype: Bool
-
-    .. note::
-        This method has quadratic runtime complexity :math:`O(n^2)`
-
-.. py:function:: check_mesh(input_files)
-    :noindex:
-
-    Checks if no triangles of the polyhedral mesh are degenerate 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!
-
-    :param List[str] input_files: polyhedral source files
-    :return: True if no triangle is degenerate and the polyhedron's plane unit normals are all pointing outwards.
-    :rtype: Bool
-
-    .. note::
-        This method has quadratic runtime complexity :math:`O(n^2)`.
+Cached Evaluation
+~~~~~~~~~~~~~~~~~
 
+.. autoclass:: polyhedral_gravity.GravityEvaluable
+   :members:
+   :special-members: __init__, __call__, __repr__
\ No newline at end of file
diff --git a/docs/background/evaluable_vs_evaluate.rst b/docs/background/evaluable_vs_evaluate.rst
index 100acd2..9611d18 100644
--- a/docs/background/evaluable_vs_evaluate.rst
+++ b/docs/background/evaluable_vs_evaluate.rst
@@ -7,27 +7,25 @@ The :code:`GravityEvaluable` was introduced with version :code:`v2.0`.
 Hence, the `project report <https://mediatum.ub.tum.de/doc/1695208/1695208.pdf>`__ does not contain any further details.
 This sections closes the knowledge gap and presents a runtime comparison.
 
-Below is a comparison of the performance of the standalone free function and the evaluable class.
-The benchmark was conducted with a M1 Pro 10-Core CPU (ARM64) using the TBB backend.
+Below is a comparison of the performance of the standalone free function and the evaluable class, as of :code:`v3.0`.
+The benchmark was conducted with a M1 Pro 10-Core CPU (ARM64) using the *TBB backend*.
 The calculation consisted each of 1000 points for the mesh of Eros (24235 vertices and 14744 faces).
-The results are shown in the table below (the lower the better/ faster).
+The points are either handed over one-by-one as :math:`(3)`-arrays or as one :math:`(1000, 3)`-array.
+The former user case appears when the latter locations are unknown beforehand - like in a trajectory propagation scenario.
+The results are shown in the plot below (the lower the better/ faster).
+
 Basically, as soon as you have more than one point to evaluate, the evaluable class is faster and
 thus recommended. This result comes from the fact that the polyhedral data is cached on the C++ side
 and thus does not need to be converted from Python to C++ for every call. Further, the evaluable class
 also caches the normals and the volume of the polyhedron, which is not the case for the standalone function.
+As of :code:`v3.0`, the evaluate free function is nearly as good as the :code:`GravityEvaluable` if all points are
+known beforehand, since the :code:`Polyhedron` class and its construction are now separated from the
+evaluation of the gravity model.
 
-+----------------------------------------+-------------------------------+
-| Test                                   | Time Per Point (microseconds) |
-+========================================+===============================+
-| Free Function (1000 times 1 point)     | 7765.073                      |
-+----------------------------------------+-------------------------------+
-| Free Function (1 time 1000 points)     | 275.917                       |
-+----------------------------------------+-------------------------------+
-| GravityEvaluable (1000 times 1 Point)  | 313.408                       |
-+----------------------------------------+-------------------------------+
-| GravityEvaluable (1 time 1000 Points)  | 253.031                       |
-+----------------------------------------+-------------------------------+
+.. image:: /../_static/runtime_measurements.png
+  :align: center
+  :width: 100%
+  :alt: Runtime Measurements for the mesh of (433) Eros in :code:`v3.0`
 
-As takeaway, one should always use the :code:`GravityEvaluable` due to its caching
-properties. The free function should only be used for few points/ quickly trying
-out something.
\ No newline at end of file
+If the scenario is yet to be determined, use the :code:`GravityEvaluable` due to its caching
+properties. If you know all points "beforehand", the approach does not matter!
\ No newline at end of file
diff --git a/docs/background/mesh_integrity.rst b/docs/background/mesh_integrity.rst
new file mode 100644
index 0000000..dd9dda9
--- /dev/null
+++ b/docs/background/mesh_integrity.rst
@@ -0,0 +1,76 @@
+.. _mesh-integrity-check:
+
+Mesh Integrity & Explanation
+============================
+
+The polyhedral gravity model by *Tsoulis et al.* is defined using
+a polyhedron with :code:`OUTWARDS` pointing plane unit normals.
+The results are negated if the plane unit normals are :code:`INWARDS` pointing.
+
+This property correlates to the sorting order of the vertices in the faces.
+In the context of most programs, one uses here the terminology *clockwise*
+and *anti-clockwise* ordering of vertices in the definition of the faces.
+The classification also depends on the coordinate system's definition (right- or left-handed).
+
+To be more illustrative and independent of the concrete coordinate system, we stick
+to the orientation of the normals: :code:`OUTWARDS` and :code:`INWARDS` pointing.
+
+For example, given you have the following cube:
+
+.. code-block:: python
+
+    import numpy as np
+
+    CUBE_VERTICES = np.array(
+      [[-1, -1, -1], [1, -1, -1], [1, 1, -1], [-1, 1, -1],
+       [-1, -1, 1], [1, -1, 1], [1, 1, 1], [-1, 1, 1]]
+    )
+    CUBE_FACES_OUTWARDS = np.array(
+      [[1, 3, 2], [0, 3, 1], [0, 1, 5], [0, 5, 4], [0, 7, 3], [0, 4, 7],
+       [1, 2, 6], [1, 6, 5], [2, 3, 6], [3, 7, 6], [4, 5, 6], [4, 6, 7]]
+    )
+    cube_inwards = Polyhedron(
+      polyhedral_source=(CUBE_VERTICES, CUBE_FACES_OUTWARDS),
+      normal_orientation=NormalOrientation.OUTWARDS,
+    )
+
+Its plane unit normals are pointing :code:`OUTWARDS`. However, if you simply invert the order like
+in the following, they will point :code:`INWARDS`:
+
+.. code-block:: python
+
+    CUBE_FACES_INWARDS = CUBE_FACES_OUTWARDS[:, [1, 0, 2]]
+    cube_outwards = Polyhedron(
+      polyhedral_source=(CUBE_VERTICES, CUBE_FACES_OUTWARDS),
+      normal_orientation=NormalOrientation.INWARDS,
+    )
+
+And then, there is still the possibility that one can define a polyhedron by hand
+with mixed orientations! This would lead to entirely different (wrong) results.
+
+So, what can we do?
+
+The :code:`Polyhedron` class checks the pointing direction of every single
+normal. This way, the `Polyhedron` ensures correct results even if a mistake occurs during the definition.
+
+Since we are in 3D and might have concave and convex polyhedrons,
+the viable option is the `Möller–Trumbore intersection algorithm <https://en.wikipedia.org/wiki/Möller–Trumbore_intersection_algorithm>`__.
+It checks the amount of intersections the plane unit normal has with the polyhedron.
+If its an even number, the normals is :code:`OUTWARDS` pointing, otherwise :code:`INWARDS`.
+In the *current implementation*, we implement a naiv version
+which takes :math:`O(n^2)` operations - which can get quite expensive for polyhedrons with many faces.
+
+To make this as straightforward to use as possible, we provide four options
+for the construction of a polyhedron in the form of the enum :code:`PolyhedronIntegrity`:
+
++-----------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+----------------------------+----------------+
+| :code:`PolyhedronIntegrity` | Meaning                                                                                                                                                                                                                 | Invalid Polyhedron Possible| Overhead       |
++=============================+=========================================================================================================================================================================================================================+============================+================+
+| :code:`AUTOMATIC` (default) | Throws an Exception if the :code:`NormalOrientation` is different or inconsistent than specified AND Prints a short-version of the above explanation to :code:`stdout` informing the user about :math:`O(n^2)` runtime  | NO                         | :math:`O(n^2)` |
++-----------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+----------------------------+----------------+
+| :code:`VERIFY`              | Throws an Exception if the :code:`NormalOrientation` is different or inconsistent than specified                                                                                                                        | NO                         | :math:`O(n^2)` |
++-----------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+----------------------------+----------------+
+| :code:`HEAL`                | Modifies the specified :code:`NormalOrientation` and the faces array to a consistent polyhedron if they are inconsistent                                                                                                | NO                         | :math:`O(n^2)` |
++-----------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+----------------------------+----------------+
+| :code:`DISABLE`             | Disables all checks                                                                                                                                                                                                     | YES                        | None           |
++-----------------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+----------------------------+----------------+
diff --git a/docs/background/rational.rst b/docs/background/rational.rst
index 92a3a6d..175a28b 100644
--- a/docs/background/rational.rst
+++ b/docs/background/rational.rst
@@ -6,7 +6,7 @@ The Rational
 
 .. image:: /../_static/eros_010.png
   :align: center
-  :width: 400
+  :width: 75%
   :alt: Downscaled mesh of (433) Eros to 10% of its original vertices and faces.
 
 Downscaled mesh of (433) Eros to 10% of its original vertices and faces.
diff --git a/docs/conf.py b/docs/conf.py
index cadb0cc..5772225 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -30,9 +30,14 @@ def configure_doxyfile(input_dir, output_dir):
         file.write(filedata)
 
 
-# Check if we're running on Read the Docs' servers or in the GitHub Actions Workflow
+# If you build the docs locally, but not via CMake, first set the environment variable BUILD_DOCS_CLI
 breathe_projects = {}
-if os.environ.get("READTHEDOCS", None) is not None or os.environ.get("GITHUB_PAGES_BUILD", None) is not None:
+if any([
+    os.environ.get("READTHEDOCS", None),
+    os.environ.get("GITHUB_PAGES_BUILD", None),
+    os.environ.get("BUILD_DOCS_CLI", None)]
+):
+    # These lines assume that we are in  /<repo-root>/docs/ folder
     input_dir = "../src"
     output_dir = "build"
     configure_doxyfile(input_dir, output_dir)
@@ -50,7 +55,7 @@ def configure_doxyfile(input_dir, output_dir):
 # Add any Sphinx extension module names here, as strings. They can be
 # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
 # ones.
-extensions = ["breathe"]
+extensions = ["breathe", "sphinx.ext.napoleon", "sphinx.ext.autodoc"]
 
 # Add any paths that contain templates here, relative to this directory.
 templates_path = ["_templates"]
@@ -70,7 +75,7 @@ def configure_doxyfile(input_dir, output_dir):
 # Add any paths that contain custom static files (such as style sheets) here,
 # relative to this directory. They are copied after the builtin static files,
 # so a file named "default.css" will overwrite the builtin "default.css".
-html_static_path = ["_static", "figures/eros_010.png"]
+html_static_path = ["_static", "figures/eros_010.png", "figures/runtime_measurements.png"]
 
 # Breathe Configuration
 breathe_default_project = "polyhedral-gravity-model"
diff --git a/docs/figures/runtime_measurements.png b/docs/figures/runtime_measurements.png
new file mode 100644
index 0000000..d91c365
--- /dev/null
+++ b/docs/figures/runtime_measurements.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:52ad95f6696e14aa33a1a18e41e1e52d64ddd9a6a95694d40b9fe688f95a5719
+size 82140
diff --git a/docs/index.rst b/docs/index.rst
index 70a964c..1daead0 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -53,6 +53,7 @@ And for more details, refer to the **Python API** or **C++ API**.
 
    background/rational
    background/approach
+   background/mesh_integrity
    background/evaluable_vs_evaluate
 
 .. toctree::
@@ -66,7 +67,6 @@ And for more details, refer to the **Python API** or **C++ API**.
    :maxdepth: 2
 
    api/model
-   api/calculation
    api/input
    api/output
    api/util
diff --git a/docs/quickstart/examples_cpp.rst b/docs/quickstart/examples_cpp.rst
index 8541b92..40b9135 100644
--- a/docs/quickstart/examples_cpp.rst
+++ b/docs/quickstart/examples_cpp.rst
@@ -42,7 +42,8 @@ Further one must specify the name of the .csv output file.
         density: 2670.0                             # constant density, units must match with the mesh (see section below)
         points: # Location of the computation point(s) P
           - [ 0, 0, 0 ]                             # Here it is situated at the origin
-        check_mesh: true                            # Fully optional, enables input checking (not given: false)
+        check_mesh: true                            # Fully optional, enables mesh autodetect+repair of
+                                                    # the polyhedron's vertex ordering (not given: true)
       output:
         filename: "gravity_result.csv"              # The name of the output file
 
@@ -64,26 +65,29 @@ 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).
+Further, we disable the parallelization using the optional fourth parameter (which defaults to true)
+and we check the if the plane unit normals are actually outwards pointing
+(costing :math:`O(n^2)`).
 
 .. 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 = ...
 
         // 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);
+        Polyhedron polyhedron{vertices, faces, density, PolyhedronIntegrity::VERIFY};
+        const auto[pot, acc, tensor] = GravityModel::evaluate(polyhedron, 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).
+(which defaults to true). We disable the :math:`O(n^2)` check if the
+plane unit normals are actually outwards pointing.
 
 .. code-block:: cpp
 
@@ -95,7 +99,8 @@ Further, we explicitly enable the parallelization using the optional fourth para
 
         // 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);
+        Polyhedron polyhedron{files, density, NormalOrientation::OUTWARDS, PolyhedronIntegrity::DISABLE};
+        const auto[pot, acc, tensor] = GravityModel::evaluate(polyhedron, point, true);
 
 
 **Example 3:** Evaluating the gravity model for a given configuration
@@ -105,19 +110,22 @@ from a .yaml file.
 
         // Reading the configuration from a yaml file
         std::shared_ptr<ConfigSource> config = std::make_shared<YAMLConfigReader>("config.yaml");
-        Polyhedron poly = config->getDataSource()->getPolyhedron();
+        auto polyhedralSource = config->getDataSource()->getPolyhedralSource();
         double density = config->getDensity();
+        PolyhedronIntegrity checkPolyhedralInput = config->getMeshInputCheckStatus() ? PolyhedronIntegrity::HEAL : PolyhedronIntegrity::DISABLE;
         // This time, we use multiple points
         std::vector<std::array<double, 3>> points = config->getPointsOfInterest();
 
         // 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);
+        Polyhedron polyhedron{polyhedralSource, density, NormalOrientation::OUTWARDS, checkPolyhedralInput};
+        const results = GravityModel::evaluate(polyhedron, points);
 
-**Example 4:** A guard statement checks that the plane unit
-normals are pointing outwards and no triangle is degenerated.
-Only use this statement if one needs clarification
-about the vertices' ordering due to its quadratic complexity!
+**Example 4:** If our :code:`Polyhedron` contains any inconsistencies in the definition, e.g.,
+some plane unit normals point outwards and some point inwards,
+the :code:`HEAL` option will fix the :code:`NormalOrientation` and any inconsistencies
+related to it. It won't throw an exception.
+The result will always be fine.
 
 .. code-block:: cpp
 
@@ -127,10 +135,8 @@ about the vertices' ordering due to its quadratic complexity!
         double density = config->getDensity();
         std::array<double, 3> point = config->getPointsOfInterest()[0];
 
-        // Guard statement
-        if (MeshChecking::checkTrianglesNotDegenerated(poly) && MeshChecking::checkNormalsOutwardPointing(poly)) {
-            GravityResult result = GravityModel::evaluate(poly, density, point);
-        }
+        Polyhedron polyhedron{files, density, NormalOrientation::OUTWARDS, PolyhedronIntegrity::HEAL};
+        GravityResult result = GravityModel::evaluate(poly, density, point);
 
 
 GravityEvaluable (with caching)
@@ -147,11 +153,14 @@ defined from within source code for a specific point and density.
         std::vector<std::array<size_t, 3>> faces = ...
         Polyhedron polyhedron{vertices, faces};
         double density = ...
+        Polyhedron polyhedron{vertices, faces, density, NormalOrientation::OUTWARDS, PolyhedronIntegrity::DISABLE};
+
+        // Our computation points
         std::array<double, 3> point = ...
         std::vector<std::array<double, 3>> points = ...
 
         // Instantiation of the GravityEvaluable object
-        GravityEvaluable evaluable{polyhedron, density};
+        GravityEvaluable evaluable{polyhedron};
 
         // From now, we can evaluate the gravity model for any point with
         const auto[pot, acc, tensor] = evaluable(point);
diff --git a/docs/quickstart/examples_python.rst b/docs/quickstart/examples_python.rst
index d2d738a..959c7e1 100644
--- a/docs/quickstart/examples_python.rst
+++ b/docs/quickstart/examples_python.rst
@@ -7,7 +7,7 @@ Details about mesh and input units can be found in :ref:`quick-start-io`.
 
 The use of the Python interface is pretty straight-forward since
 there is only one method: :code:`evaluate(..)` or the alternative
-class :code:`GravityEvaluable` caching the polyhedron.
+class :code:`GravityEvaluable` caching the :code:`Polyhedron` and intermediate results.
 Have a look at :ref:`evaluable-vs-eval` for further
 details about the difference.
 If you strive for maximal performance, use the the class :code:`GravityEvaluable`.
@@ -16,50 +16,19 @@ The polyhedral source can either be a tuple of vertices and faces, or
 a list of polyhedral mesh files (see :ref:`supported-polyhedron-source-files`).
 
 The method calls follow the same pattern as the C++ interface.
-
-.. 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
-    )
+You always define a :code:`Polyhedron` and then evaluate the polyhedral
+gravity model using either the :code:`GravityEvaluable` or :code:`evaluate(..)`.
 
 
 Free function
 -------------
 
-**Example 1:** Evaluating the gravity model for a given polyhedron
-defined from within source code for a specific point and density.
+**Example 1:** Evaluating the gravity model for a given constant density
+polyhedron defined from within source code for a single point.
 
 .. code-block:: python
 
-        import polyhedral_gravity as model
+        from polyhedral_gravity import Polyhedron, evaluate
 
         # Defining every input parameter in the source code
         vertices = ...          # (N-3)-array-like of type float
@@ -67,20 +36,24 @@ defined from within source code for a specific point and density.
         density = ...           # float
         computation_point = ... # (3)-array-like
 
-        # Evaluate the gravity model
+        # Create a constant density polyhedron & Evaluate the gravity model
         # Notice that the third argument could also be a list of points
         # Returns a tuple of potential, acceleration and tensor
         # If computation_point would be a (N,3)-array, the output would be list of triplets!
-        potential, acceleration, tensor = model.evaluate((vertices, faces), density, computation_point, parallel=True)
+        polyhedron = Polyhedron((vertices, faces), density)
+        potential, acceleration, tensor = evaluate(polyhedron, computation_point, parallel=True)
 
 
-**Example 2a:** Evaluating the gravity model for a given polyhedron
-in some source files for several points and density.
+**Example 2a:** Evaluating the gravity model for a given constant density polyhedron
+in some source files for several points.
 Of course, you can also use keyword arguments for the parameters.
+We also check that the polyhedron's plane unit normals are actually
+outwards pointing. This will raise a ValueError if at least on
+plane unit normal is inwards pointing.
 
 .. code-block:: python
 
-        import polyhedral_gravity as model
+    from polyhedral_gravity import Polyhedron, evaluate, PolyhedronIntegrity, NormalOrientation
 
         # Reading the vertices and files from a .node and .face file
         file_vertices = '___.node'      # str, path to file
@@ -91,20 +64,28 @@ Of course, you can also use keyword arguments for the parameters.
         # Evaluate the gravity model
         # Notice that the last argument could also be a list of points
         # Returns a list of tuple of potential, acceleration and tensor
-        results = model.evaluate(
+        polyhedron = Polyhedron(
             polyhedral_source=[file_vertices, file_nodes],
             density=density,
+            normal_orientation=NormalOrientation.OUTWARDS,
+            integrity_check=PolyhedronIntegrity.VERIFY,
+        )
+        results = evaluate(
+            polyhedron=polyhedron,
             computation_points=computation_points,
-            parallel=True
+            parallel=True,
         )
 
 
-**Example 2b:** Evaluating the gravity model for a given polyhedron
-in some source files for a specific point and density.
+**Example 2b:** Evaluating the gravity model for a given constant density polyhedron
+in some source files for a specific point.
+We also check that the polyhedron's plane unit normals are actually
+outwards pointing. We don't specify this here, as :code:`OUTWARDS` is the default.
+It will raise a ValueError if at least on plane unit normal is inwards pointing.
 
 .. code-block:: python
 
-        import polyhedral_gravity as model
+    from polyhedral_gravity import Polyhedron, evaluate, PolyhedronIntegrity
 
         # Reading the vertices and files from a single .mesh file
         file = '___.mesh'       # str, path to file
@@ -115,21 +96,24 @@ in some source files for a specific point and density.
         # Notice that the last argument could also be a list of points
         # Returns a tuple of potential, acceleration and tensor
         # If computation_point would be a (N,3)-array, the output would be list of triplets!
-        potential, acceleration, tensor = model.evaluate([mesh], density, computation_point)
+        polyhedron = Polyhedron(
+            polyhedral_source=[file],
+            density=density,
+            integrity_check=PolyhedronIntegrity.VERIFY,
+        )
+        potential, acceleration, tensor = evaluate(polyhedron, computation_point)
 
 
 For example 2a and 2b, refer to :ref:`supported-polyhedron-source-files` to view the available
 options for polyhedral input.
 
-**Example 3:** A guard statement checks that the plane unit
-normals are pointing outwards and no triangular surface is degenerated.
-Only use this statement if one needs clarification
-about the vertices' ordering due to its quadratic complexity!
+**Example 3a:** Here explicitly disable the security check.
+We **won't get an exception** if the plane unit normals are not
+oriented as specified, **but we also don't pay for the check with quadratic runtime complexity!**
 
 .. code-block:: python
 
-    import polyhedral_gravity as model
-    import polyhedral_gravity.utility as mesh_sanity
+    from polyhedral_gravity import Polyhedron, evaluate, PolyhedronIntegrity, NormalOrientation
 
     # Defining every input parameter in the source code
     vertices = ...          # (N-3)-array-like of type float
@@ -137,15 +121,44 @@ about the vertices' ordering due to its quadratic complexity!
     density = ...           # float
     computation_point = ... # (3)-array-like
 
+    # Evaluate the gravity model
+    # Returns a tuple of potential, acceleration and tensor
+    # If computation_point would be a (N,3)-array, the output would be list of triplets!
+    polyhedron = Polyhedron(
+        polyhedral_source=(vertices, faces),
+        density=density,
+        normal_orientation=NormalOrientation.OUTWARDS,
+        integrity_check=PolyhedronIntegrity.DISABLE,
+    )
+    potential, acceleration, tensor = evaluate(polyhedron, computation_point)
 
 
-    # Additional guard statement to check that the plane normals
-    # are outwards pointing
-    if mesh_sanity.check_mesh(vertices, faces):
-        # Evaluate the gravity model
-        # Returns a tuple of potential, acceleration and tensor
-        # If computation_point would be a (N,3)-array, the output would be list of triplets!
-        potential, acceleration, tensor = model.evaluate((vertices, faces), density, computation_point)
+**Example 3b:** Here we use the :code:`HEAL` option.
+This guarantees a valid polyhedron. But the ordering of the faces array and
+the normal_orientation might differ.
+And we also need to pay the additional quadratic runtime for the checking algorithmus.
+
+.. code-block:: python
+
+    from polyhedral_gravity import Polyhedron, evaluate, PolyhedronIntegrity, NormalOrientation
+
+    # Defining every input parameter in the source code
+    vertices = ...          # (N-3)-array-like of type float
+    faces = ...             # (N-3)-array-like of type int
+    density = ...           # float
+    computation_point = ... # (3)-array-like
+
+    # Actually, the normal_orientation doesn't matter! We could the argument
+    # as HEAL guarantees a valid polyhedron
+    # but the polyhedron might different properties polyhedron.faces
+    # and polyhedron.normal_orientation than specified
+    polyhedron = Polyhedron(
+        polyhedral_source=(vertices, faces),
+        density=density,
+        normal_orientation=NormalOrientation.OUTWARDS,
+        integrity_check=PolyhedronIntegrity.HEAL,
+    )
+    potential, acceleration, tensor = evaluate(polyhedron, computation_point)
 
 
 GravityEvaluable
@@ -161,7 +174,7 @@ Have a look at the example below to see how to use the :code:`GravityEvaluable`
 
 .. code-block:: python
 
-        import polyhedral_gravity as model
+        from polyhedral_gravity import Polyhedron, GravityEvaluable, evaluate, PolyhedronIntegrity
 
         # Defining every input parameter in the source code
         vertices = ...           # (N-3)-array-like of type float
@@ -169,8 +182,15 @@ Have a look at the example below to see how to use the :code:`GravityEvaluable`
         density = ...            # float
         computation_points = ... # (N,3)-array-like
 
+        # Definition of the Polyhedron in previous examples
+        polyhedron = Polyhedron(
+            polyhedral_source=(vertices, faces),
+            density=density,
+            integrity_check=PolyhedronIntegrity.HEAL,
+        )
+
         # Create the evaluable object
-        evaluable = model.GravityEvaluable(polyhedral_source, density)
+        evaluable = GravityEvaluable(polyhedron, density)
 
         for point in computation_points:
             # Evaluate the gravity model for single points (3)-array-like
diff --git a/docs/quickstart/overview.rst b/docs/quickstart/overview.rst
index 8b978b9..3a85b0d 100644
--- a/docs/quickstart/overview.rst
+++ b/docs/quickstart/overview.rst
@@ -32,12 +32,12 @@ mesh files.
 
 .. note::
 
-    The plane unit normals of every face of the polyhedral mesh must point **outwards**
-    of the polyhedron!
-    You can check this property via :ref:`mesh-checking-cpp` in C++ or
-    via the :ref:`mesh-checking-python` submodule in Python.
-    If the vertex order of the faces is inverted, i.e. the plane unit normals point
-    inwards, then the sign of the output will be inverted.
+    The plane unit normals of every face of the polyhedral mesh must point
+    consistently **outwards** or **inwards** the polyhedron!
+    This property is automatically enforced by the class :code:`Polyhedron` in
+    both C++ library and Python interface as long as not explicit set to :code:`DISABLE`.
+    Setting this to :code:`DISABLE`, is recommended for advanced users or when you "know your mesh".
+    For details, refer to the APIs or :ref:`mesh-integrity-check`.
 
 
 Gravity Model Output
diff --git a/docs/requirements.txt b/docs/requirements.txt
index 4aa9922..43f9e87 100644
--- a/docs/requirements.txt
+++ b/docs/requirements.txt
@@ -1,2 +1,2 @@
 breathe
-sphinx-book-theme
+sphinx-book-theme
\ No newline at end of file
diff --git a/script/cube_plots.py b/script/cube_plots.py
index 2a59606..6d5a848 100644
--- a/script/cube_plots.py
+++ b/script/cube_plots.py
@@ -1,5 +1,5 @@
 import numpy as np
-import polyhedral_gravity as gravity_model
+from polyhedral_gravity import evaluate, Polyhedron, PolyhedronIntegrity
 
 import mesh_plotting
 
@@ -32,6 +32,12 @@
 DENSITY = 1.0
 VALUES = np.arange(-2, 2.01, 0.01)
 
+cube = Polyhedron(
+    polyhedral_source=(cube_vertices, cube_faces),
+    density=DENSITY,
+    integrity_check=PolyhedronIntegrity.DISABLE,
+)
+
 ########################################################################################################################
 # Triangulation
 ########################################################################################################################
@@ -43,7 +49,7 @@
 # Plot of the potential and Acceleration in XY Plane (Z = 0)
 ########################################################################################################################
 computation_points = np.array(np.meshgrid(VALUES, VALUES, [0])).T.reshape(-1, 3)
-gravity_results = gravity_model.evaluate(cube_vertices, cube_faces, DENSITY, computation_points)
+gravity_results = evaluate(cube, computation_points)
 
 potentials = -1 * np.array([i[0] for i in gravity_results])
 potentials = potentials.reshape((len(VALUES), len(VALUES)))
@@ -58,7 +64,6 @@
 # We just want a slice of x, y values for z = 0
 accelerations = np.array([i[1][:] for i in gravity_results])
 acc_xy = np.delete(accelerations, 2, 1)
-acc_xy = -1 * acc_xy
 
 mesh_plotting.plot_quiver(X, Y, acc_xy, "Acceleration in $x$ and $y$ direction for $z=0$",
                           "../figures/cube/cube_field.png", plot_rectangle=True)
diff --git a/script/eros_plots.py b/script/eros_plots.py
index 7088237..b256360 100644
--- a/script/eros_plots.py
+++ b/script/eros_plots.py
@@ -1,5 +1,5 @@
 import numpy as np
-import polyhedral_gravity as gravity_model
+from polyhedral_gravity import Polyhedron, evaluate, PolyhedronIntegrity
 import mesh_plotting
 import mesh_utility
 
@@ -11,6 +11,12 @@
 DENSITY = 1.0
 VALUES = np.arange(-2, 2.01, 0.01)
 
+eros = Polyhedron(
+    polyhedral_source=(vertices, faces),
+    density=DENSITY,
+    integrity_check=PolyhedronIntegrity.DISABLE,
+)
+
 ########################################################################################################################
 # Triangulation
 ########################################################################################################################
@@ -22,7 +28,7 @@
 # Plot of the potential and Acceleration in XY Plane (Z = 0)
 ########################################################################################################################
 computation_points = np.array(np.meshgrid(VALUES, VALUES, [0])).T.reshape(-1, 3)
-gravity_results = gravity_model.evaluate(vertices, faces, DENSITY, computation_points)
+gravity_results = evaluate(eros, computation_points)
 
 potentials = -1 * np.array([i[0] for i in gravity_results])
 potentials = potentials.reshape((len(VALUES), len(VALUES)))
@@ -38,7 +44,6 @@
 # We just want a slice of x, y values for z = 0
 accelerations = np.array([i[1][:] for i in gravity_results])
 acc_xy = np.delete(accelerations, 2, 1)
-acc_xy = -1 * acc_xy
 
 mesh_plotting.plot_quiver(X, Y, acc_xy, "Acceleration in $x$ and $y$ direction for $z=0$",
                           "../figures/eros/eros_field_xy.png", vertices=vertices, coordinate=2)
@@ -47,7 +52,7 @@
 # Plot of the potential and Acceleration in XZ Plane (Y = 0)
 ########################################################################################################################
 computation_points = np.array(np.meshgrid(VALUES, [0], VALUES)).T.reshape(-1, 3)
-gravity_results = gravity_model.evaluate(vertices, faces, DENSITY, computation_points)
+gravity_results = evaluate(eros, computation_points)
 
 potentials = -1 * np.array([i[0] for i in gravity_results])
 potentials = potentials.reshape((len(VALUES), len(VALUES)))
@@ -65,7 +70,6 @@
 # We just want a slice of x, y values for z = 0
 accelerations = np.array([i[1][:] for i in gravity_results])
 acc_xy = np.delete(accelerations, 1, 1)
-acc_xy = -1 * acc_xy
 
 mesh_plotting.plot_quiver(X, Z, acc_xy, "Acceleration in $x$ and $z$ direction for $y=0$",
                           "../figures/eros/eros_field_xz.png", labels=('$x$', '$z$'),
@@ -75,7 +79,7 @@
 # Plot of the potential and Acceleration in YZ Plane (X = 0)
 ########################################################################################################################
 computation_points = np.array(np.meshgrid([0], VALUES, VALUES)).T.reshape(-1, 3)
-gravity_results = gravity_model.evaluate(vertices, faces, DENSITY, computation_points)
+gravity_results = evaluate(eros, computation_points)
 
 potentials = -1 * np.array([i[0] for i in gravity_results])
 potentials = potentials.reshape((len(VALUES), len(VALUES)))
@@ -93,7 +97,6 @@
 # We just want a slice of x, y values for z = 0
 accelerations = np.array([i[1][:] for i in gravity_results])
 acc_xy = np.delete(accelerations, 0, 1)
-acc_xy = -1 * acc_xy
 
 mesh_plotting.plot_quiver(Y, Z, acc_xy, "Acceleration in $y$ and $z$ direction for $x=0$",
                           "../figures/eros/eros_field_yz.png", labels=('$y$', '$z$'),
diff --git a/script/polyhedral-gravity.ipynb b/script/polyhedral-gravity.ipynb
index 22780f9..403169c 100644
--- a/script/polyhedral-gravity.ipynb
+++ b/script/polyhedral-gravity.ipynb
@@ -47,8 +47,7 @@
     "import matplotlib.pyplot as plt\n",
     "import matplotlib.patches as patches\n",
     "from matplotlib import cm\n",
-    "import polyhedral_gravity as gravity_model\n",
-    "import polyhedral_gravity.utility as gravity_util\n",
+    "from polyhedral_gravity import Polyhedron, evaluate, PolyhedronIntegrity, NormalOrientation, GravityEvaluable\n",
     "\n",
     "# These two scripts provide some utility for plotting and the analytical solution\n",
     "# They reside in the same directory as this notebook\n",
@@ -72,9 +71,11 @@
     "\n",
     "Next, we start with an easy example evaluating the gravity tensor for points of a cubic polyhedron consisting of $12$ triangular faces, centered at the origin $(0, 0, 0)$ and with edges of $2m$.\n",
     "Therefor, we first declare the vertices and the connectivity of this polyhedron.\n",
+    "Further, we define the density of the cube as $\\rho = 1 \\frac{kg}{m^3}$.\n",
     "\n",
-    "**Important!:** The vertices in the faces need to be sorted, so that the plane normals will point outwards of the polyhedron!\n",
-    "The Python interface provides an `utility` submodule which offers the capabilities to check for this constraint. However, one shall notice that the procedure is computationally expensive in $O(N^2)$ using the Möller–Trumbore intersection algorithm.\n"
+    "**Important!:** The vertices in the faces need to be sorted, so that the plane normals will point consistently outwards/ inwards the polyhedron!\n",
+    "By default the the `Polyhedron` class will check that and even assure the property!\n",
+    "For the cube, we use the guranteed version using `HEAL`. This might modify the order of the vertices in the faces, but gurantees that the normals point outwards/ inwards consistently.\n"
    ]
   },
   {
@@ -90,7 +91,7 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Plane normals point outwards: True\n"
+      "<polyhedral_gravity.Polyhedron, density = 1, vertices = 8, faces = 12, orientation = OUTWARDS>\n"
      ]
     }
    ],
@@ -124,9 +125,13 @@
     "        [4, 6, 7],\n",
     "    ]\n",
     ")\n",
-    "\n",
-    "outwards_pointing_normals = gravity_util.check_mesh(cube_vertices, cube_faces)\n",
-    "print(f\"Plane normals point outwards: {outwards_pointing_normals}\")"
+    "density = 1.0  # [kg/m^3]\n",
+    "cube_polyhedron = Polyhedron(\n",
+    "    polyhedral_source=(cube_vertices, cube_faces),\n",
+    "    density=density,\n",
+    "    integrity_check=PolyhedronIntegrity.HEAL,\n",
+    ")\n",
+    "print(cube_polyhedron)"
    ]
   },
   {
@@ -137,7 +142,7 @@
     }
    },
    "source": [
-    "Next we define the wished computation point $(1, 0, 0)$ and the density in $\\frac{kg}{m^3}$\n"
+    "Next we define the wished computation point $(1, 0, 0)$\n"
    ]
   },
   {
@@ -150,8 +155,7 @@
    },
    "outputs": [],
    "source": [
-    "computation_point = [1, 0, 0]  # [m]\n",
-    "density = 1.0  # [kg/m^3]"
+    "computation_point = [1, 0, 0]  # [m]"
    ]
   },
   {
@@ -184,15 +188,14 @@
      "text": [
       "=============Gravity Model Solution=============\n",
       "Potential: 4.786301362419239e-10 m^2/s^2\n",
-      "Acceleration [Vx, Vy, Vz]: [3.466493366453959e-10, 5.927969226604545e-26, 0.0] m/s^2\n",
+      "Acceleration [Vx, Vy, Vz]: [-3.466493366453959e-10, -5.927969226604545e-26, -0.0] m/s^2\n",
       "Second derivative tensor [Vxx, Vyy, Vzz, Vxy, Vxz, Vyz]: [-5.375692682923615e-11, -1.8280085506392544e-10, -1.8280085506392544e-10, 0.0, 0.0, 8.891953839906818e-26] 1/s^2\n"
      ]
     }
    ],
    "source": [
-    "potential, acceleration, tensor = gravity_model.evaluate(\n",
-    "    polyhedral_source=(cube_vertices, cube_faces),\n",
-    "    density=density,\n",
+    "potential, acceleration, tensor = evaluate(\n",
+    "    polyhedron=cube_polyhedron,\n",
     "    computation_points=computation_point,\n",
     "    parallel=True,\n",
     ")\n",
@@ -233,7 +236,7 @@
      "text": [
       "=============Analytical Solutions=============\n",
       "Analytical Potential: 4.78630136241924E-10 m^2/s^2\n",
-      "Analytical Acceleration [Vx, Vy, Vz]: [3.46649336645396e-10 0 0] m/s^2\n",
+      "Analytical Acceleration [Vx, Vy, Vz]: [-3.46649336645396e-10 0 0] m/s^2\n",
       "=================Differences==================\n",
       "Difference Potential: 2.06795153138257E-25\n",
       "Difference Acceleration: [1.55096364853693e-25 5.92796922660455e-26 0]\n"
@@ -320,9 +323,8 @@
     "Y = VALUES\n",
     "\n",
     "computation_points = np.array(np.meshgrid(X, Y, [0])).T.reshape(-1, 3)\n",
-    "gravity_results = gravity_model.evaluate(\n",
-    "    polyhedral_source=(cube_vertices, cube_faces),\n",
-    "    density=density,\n",
+    "gravity_results = evaluate(\n",
+    "    polyhedron=cube_polyhedron,\n",
     "    computation_points=computation_points,\n",
     ")"
    ]
@@ -420,15 +422,13 @@
     "Y = VALUES\n",
     "\n",
     "computation_points = np.array(np.meshgrid(X, Y, [0])).T.reshape(-1, 3)\n",
-    "gravity_results = gravity_model.evaluate(\n",
-    "    polyhedral_source=(cube_vertices, cube_faces),\n",
-    "    density=density,\n",
+    "gravity_results = evaluate(\n",
+    "    polyhedron=cube_polyhedron,\n",
     "    computation_points=computation_points,\n",
     ")\n",
     "\n",
     "accelerations = np.array([i[1][:] for i in gravity_results])\n",
     "acc_xy = np.delete(accelerations, 2, 1)\n",
-    "acc_xy = -1 * acc_xy\n",
     "\n",
     "X = computation_points[:, 0].reshape(len(VALUES), -1)\n",
     "Y = computation_points[:, 1].reshape(len(VALUES), -1)\n",
@@ -453,7 +453,8 @@
     "\n",
     "**Important Note:** The employed mesh of Eros here is unitless! Thus, units are not a concern here. Its just a demonstration of the interface.\n",
     "\n",
-    "Further, we also make use of the `Evaluable` class, which allows us to efficiently evaluate the gravity tensor for multiple computation points at once, even with interruptions in between. Certain computed properties are cached and reused.\n"
+    "Further, we also make use of the `GravityEvaluable` class, which allows us to efficiently evaluate the gravity tensor for multiple computation points at once, even with interruptions in between. Certain computed properties are cached and reused.\n",
+    "We also know this mesh and that all normals are `OUTWARDS` pointing. So we don't need the mesh check and `DISABLE` it.\n"
    ]
   },
   {
@@ -469,25 +470,31 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Potential: 1.641219381390107e-07\n",
-      "Acceleration [Vx, Vy, Vz]: [2.410422444002989e-09, 3.9101840968470284e-08, 9.039023451764005e-09]\n",
-      "Second derivative tensor [Vxx, Vyy, Vzz, Vxy, Vxz, Vyz]: [-1.224392128772979e-07, -1.1414944948899626e-06, -9.754414135835843e-07, -3.86381487022689e-08, 1.354813459901532e-08, -2.6597490639550824e-08]\n"
+      "Potential: 1.6412193813901074e-07\n",
+      "Acceleration [Vx, Vy, Vz]: [-2.4104224440030084e-09, -3.9101840968470205e-08, -9.039023451764084e-09]\n",
+      "Second derivative tensor [Vxx, Vyy, Vzz, Vxy, Vxz, Vyz]: [-1.2243921287729794e-07, -1.1414944948899624e-06, -9.75441413583584e-07, -3.8638148702268895e-08, 1.354813459901531e-08, -2.6597490639550805e-08]\n"
      ]
     }
    ],
    "source": [
     "eros_vertices_path = \"mesh/Eros.node\"\n",
     "eros_faces_path = \"mesh/Eros.face\"\n",
-    "\n",
     "density = 2670.0\n",
+    "\n",
+    "eros_polyhedron = Polyhedron(\n",
+    "    polyhedral_source=[eros_vertices_path, eros_faces_path],\n",
+    "    density=density,\n",
+    "    normal_orientation=NormalOrientation.OUTWARDS,\n",
+    "    integrity_check=PolyhedronIntegrity.DISABLE,\n",
+    ")\n",
+    "\n",
+    "\n",
     "computation_point = np.array([0, 0, 0])  # m\n",
     "\n",
     "# The differende betweeen array input & file input lies in the types\n",
     "# Array Mesh Input ==> tuple of arrays/ lists\n",
     "# File Mesh Input ==> list of string file paths\n",
-    "evaluable_eros = gravity_model.GravityEvaluable(\n",
-    "    polyhedral_source=[eros_vertices_path, eros_faces_path], density=density\n",
-    ")\n",
+    "evaluable_eros = GravityEvaluable(eros_polyhedron)\n",
     "\n",
     "# Call the evaluable object with the computation point(s) and optionally disable/ enable parallelization\n",
     "potential, acceleration, tensor = evaluable_eros(computation_points=computation_point, parallel=True)\n",
@@ -614,7 +621,6 @@
    "source": [
     "accelerations = np.array([i[1][:] for i in gravity_results])\n",
     "acc_xy = np.delete(accelerations, 2, 1)\n",
-    "acc_xy = -1 * acc_xy\n",
     "\n",
     "mesh_plotting.plot_quiver(X, Y, acc_xy, \"Acceleration in $x$ and $y$ direction for $z=0$\")"
    ]
@@ -631,7 +637,8 @@
     "\n",
     "All good thrings come in threes, so let's try another example.\n",
     "\n",
-    "Here, we use some additional utility to read in the vertices and triangles from the `.pk` file since this file format is not directly supported through the python interface. Additionally, we can plot the Torus, since we have the \"raw\" vertices and faces directly in Python.\n"
+    "Here, we use some additional utility to read in the vertices and triangles from the `.pk` file since this file format is not directly supported through the python interface. Additionally, we can plot the Torus, since we have the \"raw\" vertices and faces directly in Python.\n",
+    "Agian, we know this mesh and that all normals are `OUTWARDS` pointing. So we don't need the mesh check and `DISABLE` it.\n"
    ]
   },
   {
@@ -666,7 +673,12 @@
     "\n",
     "torus_vertices, torus_faces = mesh_utility.read_pk_file(\"mesh/torus_lp.pk\")\n",
     "density = 1.0\n",
-    "torus_evaluable = gravity_model.GravityEvaluable(polyhedral_source=(torus_vertices, torus_faces), density=density)\n",
+    "torus_polyhedron = Polyhedron(\n",
+    "    (torus_vertices, torus_faces), density, NormalOrientation.OUTWARDS, PolyhedronIntegrity.DISABLE\n",
+    ")\n",
+    "\n",
+    "\n",
+    "torus_evaluable = GravityEvaluable(torus_polyhedron)\n",
     "\n",
     "mesh_plotting.plot_triangulation(torus_vertices, torus_faces, \"Triangulation of a Torus\")"
    ]
@@ -766,8 +778,6 @@
     "accelerations = np.array([i[1][:] for i in gravity_results])\n",
     "accelerations = np.delete(accelerations, 2, 1)\n",
     "\n",
-    "accelerations = -1 * accelerations\n",
-    "\n",
     "\n",
     "fig, ax = plt.subplots()\n",
     "\n",
diff --git a/script/time.py b/script/time.py
index db56e47..492225f 100755
--- a/script/time.py
+++ b/script/time.py
@@ -1,91 +1,103 @@
 #!python3
-import polyhedral_gravity
+from polyhedral_gravity import evaluate, Polyhedron, PolyhedronIntegrity, GravityEvaluable
 import numpy as np
 import matplotlib.pyplot as plt
 import mesh_utility
 import timeit
 import pickle
-
-vertices, faces = mesh_utility.read_pk_file("/Users/schuhmaj/Programming/polyhedral-gravity-model/script/mesh/Eros.pk")
-vertices, faces = np.array(vertices), np.array(faces)
-
-# 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")
-
-
-########################################################################################################################
-# PICKLE SUPPORT
-########################################################################################################################
-
-
-with open("evaluable.pk", "wb") as f:
-    pickle.dump(evaluable, f, pickle.HIGHEST_PROTOCOL)
-    
-with open("evaluable.pk", "rb") as f:
-    evaluable2 = pickle.load(f)
-    
-
-start_time = timeit.default_timer()
-evaluable2(computation_points)
-end_time = timeit.default_timer()
-
-delta = (end_time - start_time) / N * 1e6
-# Print the time in milliseconds
-print("######## Pickle ########")
-print(f"--> 1 time {N} points with GravityEvaluable")
-print(f"--> Time taken: {delta:.3f} microseconds per point")
+from typing import Dict
+
+
+def run_time_measurements(sample_size: int = 1000) -> Dict[str, float]:
+    """Returns the RuntimeMeasurements for the four configurations as mapping.
+
+    Returns:
+        a mapping from configuration name to runtime per point in microseconds
+    """
+    results = dict()
+    vertices, faces = mesh_utility.read_pk_file("./mesh/Eros.pk")
+    vertices, faces = np.array(vertices), np.array(faces)
+    polyhedron = Polyhedron(
+        polyhedral_source=(vertices, faces),
+        density=1.0,
+        integrity_check=PolyhedronIntegrity.DISABLE,
+    )
+
+    # Generate 1000 random cartesian points
+    computation_points = np.random.uniform(-2, 2, (sample_size, 3))
+
+    start_time = timeit.default_timer()
+    for i in range(sample_size):
+        evaluate(polyhedron, computation_points[i])
+    end_time = timeit.default_timer()
+
+    delta = (end_time - start_time) / sample_size * 1e6
+    print("######## Old Single-Point ########")
+    print(f"--> {sample_size} times 1 point with old interface")
+    print(f"--> Time taken: {delta:.3f} microseconds per point")
+    results[f"evaluate \n ${sample_size} \\times 1$ point"] = delta
+
+    start_time = timeit.default_timer()
+    evaluate(polyhedron, computation_points)
+    end_time = timeit.default_timer()
+
+    delta = (end_time - start_time) / sample_size * 1e6
+    print("######## Old Multi-Point ########")
+    print(f"--> 1 time {sample_size} points with old interface")
+    print(f"--> Time taken: {delta:.3f} microseconds per point")
+    results[f"evaluate \n $1 \\times {sample_size}$ points"] = delta
+
+    evaluable = GravityEvaluable(polyhedron)
+
+    start_time = timeit.default_timer()
+    for i in range(sample_size):
+        evaluable(computation_points[i])
+    end_time = timeit.default_timer()
+
+    delta = (end_time - start_time) / sample_size * 1e6
+    print("######## GravityEvaluable (Single-Point) #########")
+    print(f"--> {sample_size} times 1 point with GravityEvaluable")
+    print(f"--> Time taken: {delta:.3f} microseconds per point")
+    results[f"GravityEvaluable \n ${sample_size} \\times 1$ point"] = delta
+
+    evaluable = GravityEvaluable(polyhedron)
+
+    start_time = timeit.default_timer()
+    evaluable(computation_points)
+    end_time = timeit.default_timer()
+
+    delta = (end_time - start_time) / sample_size * 1e6
+    print("######## GravityEvaluable (Multi-Point) ########")
+    print(f"--> 1 time {sample_size} points with GravityEvaluable")
+    print(f"--> Time taken: {delta:.3f} microseconds per point")
+    results[f"GravityEvaluable \n $1 \\times {sample_size}$ points"] = delta
+    return results
+
+
+def create_plot(runtime_results: Dict[str, float], sample_size: int = 1000) -> None:
+    """Creates a bar chart plot with given runtime results dictionary"""
+    # Legend for X-axis
+    configurations = list(runtime_results.keys())
+    # Figure of runtime measurements per call (in microseconds)
+    runtime_per_call = list(runtime_results.values())
+
+    # Create a numpy range equal to configurations list length
+    x_pos = np.arange(len(configurations))
+
+    fig, ax = plt.subplots(figsize=(7, 4))
+
+    # Now, use the cmap to pick colors for each bar
+    ax.bar(x_pos, runtime_per_call, align='center', color=["darkorange", "gold", "navy", "blue"])
+    ax.grid(True)
+    ax.set_xticks(x_pos)
+    ax.set_xticklabels(configurations)
+    ax.set_ylabel('Runtime per Point $[\mu s]$')
+    ax.set_title(f'Runtime Measurements (Sample Size = ${sample_size}$) of the Python Interface v3.0')
+
+    # Save the figure
+    plt.savefig("runtime_measurements.png", dpi=300)
+
+
+if __name__ == '__main__':
+    results = run_time_measurements()
+    create_plot(results)
\ No newline at end of file
diff --git a/script/torus_plots.py b/script/torus_plots.py
index 6917dd1..2fc8bc6 100644
--- a/script/torus_plots.py
+++ b/script/torus_plots.py
@@ -1,5 +1,5 @@
 import numpy as np
-import polyhedral_gravity as gravity_model
+from polyhedral_gravity import Polyhedron, evaluate, PolyhedronIntegrity
 import mesh_plotting
 import mesh_utility
 
@@ -11,6 +11,12 @@
 DENSITY = 1.0
 VALUES = np.arange(-2, 2.01, 0.01)
 
+torus = Polyhedron(
+    polyhedral_source=(vertices, faces),
+    density=DENSITY,
+    integrity_check=PolyhedronIntegrity.DISABLE,
+)
+
 ########################################################################################################################
 # Triangulation
 ########################################################################################################################
@@ -22,7 +28,7 @@
 # Plot of the potential and Acceleration in XY Plane (Z = 0)
 ########################################################################################################################
 computation_points = np.array(np.meshgrid(VALUES, VALUES, [0])).T.reshape(-1, 3)
-gravity_results = gravity_model.evaluate(vertices, faces, DENSITY, computation_points)
+gravity_results = evaluate(torus, computation_points)
 
 potentials = -1 * np.array([i[0] for i in gravity_results])
 potentials = potentials.reshape((len(VALUES), len(VALUES)))
@@ -39,7 +45,6 @@
 # We just want a slice of x, y values for z = 0
 accelerations = np.array([i[1][:] for i in gravity_results])
 acc_xy = np.delete(accelerations, 2, 1)
-acc_xy = -1 * acc_xy
 
 mesh_plotting.plot_quiver(X, Y, acc_xy, "Acceleration in $x$ and $y$ direction for $z=0$",
                           "../figures/torus/torus_field_xy.png", vertices=vertices, coordinate=2)
@@ -48,7 +53,7 @@
 # Plot of the potential and Acceleration in XZ Plane (Y = 0)
 ########################################################################################################################
 computation_points = np.array(np.meshgrid(VALUES, [0], VALUES)).T.reshape(-1, 3)
-gravity_results = gravity_model.evaluate(vertices, faces, DENSITY, computation_points)
+gravity_results = evaluate(torus, computation_points)
 
 potentials = -1 * np.array([i[0] for i in gravity_results])
 potentials = potentials.reshape((len(VALUES), len(VALUES)))
@@ -66,7 +71,6 @@
 # We just want a slice of x, y values for z = 0
 accelerations = np.array([i[1][:] for i in gravity_results])
 acc_xy = np.delete(accelerations, 1, 1)
-acc_xy = -1 * acc_xy
 
 mesh_plotting.plot_quiver(X, Z, acc_xy, "Acceleration in $x$ and $z$ direction for $y=0$",
                           "../figures/torus/torus_field_xz.png", labels=('$x$', '$z$'),
@@ -76,7 +80,7 @@
 # Plot of the potential and Acceleration in YZ Plane (X = 0)
 ########################################################################################################################
 computation_points = np.array(np.meshgrid([0], VALUES, VALUES)).T.reshape(-1, 3)
-gravity_results = gravity_model.evaluate(vertices, faces, DENSITY, computation_points)
+gravity_results = evaluate(torus, computation_points)
 
 potentials = -1 * np.array([i[0] for i in gravity_results])
 potentials = potentials.reshape((len(VALUES), len(VALUES)))
@@ -95,7 +99,6 @@
 # We just want a slice of x, y values for z = 0
 accelerations = np.array([i[1][:] for i in gravity_results])
 acc_xy = np.delete(accelerations, 0, 1)
-acc_xy = -1 * acc_xy
 
 mesh_plotting.plot_quiver(Y, Z, acc_xy, "Acceleration in $y$ and $z$ direction for $x=0$",
                           "../figures/torus/torus_field_yz.png", labels=('$y$', '$z$'),
diff --git a/setup.py b/setup.py
index 0bba93f..aaa5f30 100644
--- a/setup.py
+++ b/setup.py
@@ -47,7 +47,8 @@ def get_cmake_generator():
     If not, returns "Ninja" if ninja build is installed.
     Otherwise, none is returned.
     """
-    if (os_env_generator := os.environ.get("CMAKE_GENERATOR")) is not None:
+    os_env_generator = os.environ.get("CMAKE_GENERATOR", None)
+    if os_env_generator is not None:
         return os_env_generator
     elif is_ninja_installed():
         return "Ninja"
@@ -160,21 +161,27 @@ def build_extension(self, ext):
 # -----------------------------------------------------------------------------------------
 
 
+picture_in_readme = '''<p align="center">
+  <img src="paper/figures/eros_010.png" width="50%">
+  <br>
+  <em>
+    <a href="https://github.com/gomezzz/geodesyNets/blob/master/3dmeshes/eros_lp.pk">Mesh of (433) Eros</a> with 739 vertices and 1474 faces
+  </em>
+</p>'''
+
+
 # --------------------------------------------------------------------------------
 # Modify these entries to customize the package metadata of the polyhedral gravity
 # --------------------------------------------------------------------------------
 setup(
     name="polyhedral_gravity",
-    version="2.2",
+    version="3.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",
-    long_description="""
-        The package polyhedral_gravity provides a simple to use interface for the evaluation of the full gravity
-        tensor of a constant density polyhedron at given computation points according to the geodetic convention.
-        It is based on a fast, parallelized backbone in C++ capable of evaluating the gravity at thousands of
-        computation points in the fraction of a second.
-    """,
+    description="Package to compute full gravity tensor of a given constant density polyhedron for arbitrary points "
+                "according to the geodetic convention",
+    long_description=open("README.md").read().replace(picture_in_readme, ''),
+    long_description_content_type="text/markdown",
     ext_modules=[CMakeExtension("polyhedral_gravity")],
     cmdclass={"build_ext": CMakeBuild},
     license="GPLv3",
@@ -182,5 +189,23 @@ def build_extension(self, ext):
     zip_safe=False,
     python_requires=">=3.6",
     include_package_data=True,
+    project_urls={
+        "Homepage": "https://github.com/esa/polyhedral-gravity-model",
+        "Source": "https://github.com/esa/polyhedral-gravity-model",
+        "Documentation": "https://esa.github.io/polyhedral-gravity-model/",
+        "Issues": "https://github.com/esa/polyhedral-gravity-model/issues",
+        "Changelog": "https://github.com/esa/polyhedral-gravity-model/releases",
+    },
+    classifiers=[
+        "Development Status :: 5 - Production/Stable",
+        "License :: OSI Approved :: GNU General Public License v3 (GPLv3)",
+        "Programming Language :: C++",
+        "Programming Language :: Python",
+        "Operating System :: Microsoft :: Windows",
+        "Operating System :: MacOS",
+        "Operating System :: POSIX :: Linux",
+        "Intended Audience :: Science/Research",
+        "Topic :: Scientific/Engineering :: Physics",
+    ],
 )
 # --------------------------------------------------------------------------------
diff --git a/src/main.cpp b/src/main.cpp
index 42ef01a..4692e2c 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -1,8 +1,7 @@
 #include <chrono>
 #include "polyhedralGravity/input/ConfigSource.h"
 #include "polyhedralGravity/input/YAMLConfigReader.h"
-#include "polyhedralGravity/calculation/GravityModel.h"
-#include "polyhedralGravity/calculation/MeshChecking.h"
+#include "polyhedralGravity/model/GravityModel.h"
 #include "polyhedralGravity/output/Logging.h"
 #include "polyhedralGravity/output/CSVWriter.h"
 
@@ -16,32 +15,17 @@ int main(int argc, char *argv[]) {
     }
 
     try {
-
         std::shared_ptr<ConfigSource> config = std::make_shared<YAMLConfigReader>(argv[1]);
-        auto poly = config->getDataSource()->getPolyhedron();
+        auto polyhedralSource = config->getDataSource()->getPolyhedralSource();
         auto density = config->getDensity();
         auto computationPoints = config->getPointsOfInterest();
         auto outputFileName = config->getOutputFileName();
-        bool checkPolyhedralInput = config->getMeshInputCheckStatus();
-
-        // Checking that the vertices are correctly set-up in the input if activated
-        if (checkPolyhedralInput) {
-            SPDLOG_LOGGER_INFO(PolyhedralGravityLogger::DEFAULT_LOGGER.getLogger(), "Checking mesh...");
-            if (!MeshChecking::checkTrianglesNotDegenerated(poly)) {
-                throw std::runtime_error{
-                        "At least on triangle in the mesh is degenerated and its surface area equals zero!"};
-            } else if (!MeshChecking::checkNormalsOutwardPointing(poly)) {
-                throw std::runtime_error{
-                        "The plane unit normals are not pointing outwards! Please check the order "
-                        "of the vertices in the polyhedral input source!"};
-            }  else {
-                SPDLOG_LOGGER_INFO(PolyhedralGravityLogger::DEFAULT_LOGGER.getLogger(), "The mesh is fine.");
-            }
-        }
+        PolyhedronIntegrity checkPolyhedralInput = config->getMeshInputCheckStatus() ? PolyhedronIntegrity::HEAL : PolyhedronIntegrity::DISABLE;
 
         SPDLOG_LOGGER_INFO(PolyhedralGravityLogger::DEFAULT_LOGGER.getLogger(), "The calculation started...");
         auto start = std::chrono::high_resolution_clock::now();
-        auto result = GravityModel::evaluate(poly, density, computationPoints);
+        Polyhedron polyhedron{polyhedralSource, density, NormalOrientation::OUTWARDS, checkPolyhedralInput};
+        auto result = GravityModel::evaluate(polyhedron, computationPoints, true);
         auto end = std::chrono::high_resolution_clock::now();
         auto duration = end - start;
         auto ms = std::chrono::duration_cast<std::chrono::microseconds>(duration);
diff --git a/src/polyhedralGravity/calculation/GravityModel.cpp b/src/polyhedralGravity/calculation/GravityModel.cpp
deleted file mode 100644
index e95735e..0000000
--- a/src/polyhedralGravity/calculation/GravityModel.cpp
+++ /dev/null
@@ -1,24 +0,0 @@
-#include "GravityModel.h"
-
-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>
-    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));
-    }
-
-}
diff --git a/src/polyhedralGravity/calculation/MeshChecking.cpp b/src/polyhedralGravity/calculation/MeshChecking.cpp
deleted file mode 100644
index 38ac5f3..0000000
--- a/src/polyhedralGravity/calculation/MeshChecking.cpp
+++ /dev/null
@@ -1,86 +0,0 @@
-#include "MeshChecking.h"
-
-namespace polyhedralGravity::MeshChecking {
-
-    bool checkNormalsOutwardPointing(const Polyhedron &polyhedron) {
-        using namespace util;
-        auto it = transformPolyhedron(polyhedron);
-        // All normals have to point outwards (intersect the polyhedron even times)
-        return thrust::transform_reduce(
-                thrust::device,
-                it.first, it.second, [&polyhedron](const Array3Triplet &face) {
-                    // The centroid of the triangular face
-                    const Array3 centroid = (face[0] + face[1] + face[2]) / 3.0;
-
-                    // The normal of the plane calculated with two segments of the triangle
-                    const Array3 segmentVector1 = face[1] - face[0];
-                    const Array3 segmentVector2 = face[2] - face[1];
-                    const Array3 normal = util::normal(segmentVector1, segmentVector2);
-
-                    // The origin of the array has a slight offset in direction of the normal
-                    const Array3 rayOrigin = centroid + (normal * EPSILON);
-
-                    // If the ray intersects the polyhedron an even number of times then the normal points outwards
-                    const size_t intersects = detail::countRayPolyhedronIntersections(rayOrigin, normal, polyhedron);
-                    return intersects % 2 == 0;
-                }, true, thrust::logical_and<bool>());
-    }
-
-    bool checkTrianglesNotDegenerated(const Polyhedron &polyhedron) {
-        auto it = transformPolyhedron(polyhedron);
-        // All triangles surface area needs to be greater than zero
-        return thrust::transform_reduce(
-                thrust::device,
-                it.first, it.second, [](const Array3Triplet &face) {
-                    return util::surfaceArea(face) > 0.0;
-                }, true, thrust::logical_and<bool>());
-    }
-
-    size_t
-    detail::countRayPolyhedronIntersections(const Array3 &rayOrigin, const Array3 &rayVector, const Polyhedron &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) {
-            const auto intersection = rayIntersectsTriangle(rayOrigin, rayVector, triangle);
-            if (intersection != nullptr) {
-                intersections.insert(*intersection);
-            }
-        });
-        return intersections.size();
-    }
-
-    std::unique_ptr<Array3>
-    detail::rayIntersectsTriangle(const Array3 &rayOrigin, const Array3 &rayVector, const Array3Triplet &triangle) {
-        // Adapted Möller–Trumbore intersection algorithm
-        // see https://en.wikipedia.org/wiki/Möller–Trumbore_intersection_algorithm
-        using namespace util;
-        const Array3 edge1 = triangle[1] - triangle[0];
-        const Array3 edge2 = triangle[2] - triangle[0];
-        const Array3 h = cross(rayVector, edge2);
-        const double a = dot(edge1, h);
-        if (a > -EPSILON && a < EPSILON) {
-            return nullptr;
-        }
-
-        const double f = 1.0 / a;
-        const Array3 s = rayOrigin - triangle[0];
-        const double u = f * dot(s, h);
-        if (u < 0.0 || u > 1.0) {
-            return nullptr;
-        }
-
-        const Array3 q = cross(s, edge1);
-        const double v = f * dot(rayVector, q);
-        if (v < 0.0 || u + v > 1.0) {
-            return nullptr;
-        }
-
-        const double t = f * dot(edge2, q);
-        if (t > EPSILON) {
-            return std::make_unique<Array3>(rayOrigin + rayVector * t);
-        } else {
-            return nullptr;
-        }
-    }
-}
\ No newline at end of file
diff --git a/src/polyhedralGravity/calculation/MeshChecking.h b/src/polyhedralGravity/calculation/MeshChecking.h
deleted file mode 100644
index baac125..0000000
--- a/src/polyhedralGravity/calculation/MeshChecking.h
+++ /dev/null
@@ -1,62 +0,0 @@
-#pragma once
-
-#include <set>
-#include <memory>
-#include "thrust/transform_reduce.h"
-#include "thrust/execution_policy.h"
-#include "polyhedralGravity/model/Polyhedron.h"
-#include "polyhedralGravity/model/GravityModelData.h"
-#include "polyhedralGravity/util/UtilityContainer.h"
-#include "polyhedralGravity/util/UtilityConstants.h"
-#include "polyhedralGravity/calculation/PolyhedronTransform.h"
-
-namespace polyhedralGravity::MeshChecking {
-
-    /**
-     * Checks if the vertices are in such a order that the unit normals of each plane point outwards the polyhedron
-     * @param polyhedron the polyhedron consisting of vertices and triangular faces
-     * @return true if all the unit normals are pointing outwards
-     *
-     * @note This has quadratic complexity @f$O(n^2)@f$! For bigger problem sizes, K-D trees/ Octrees could improve
-     * the complexity to determine the intersections
-     * (see https://stackoverflow.com/questions/45603469/how-to-calculate-the-normals-of-a-box)
-     */
-    bool checkNormalsOutwardPointing(const Polyhedron &polyhedron);
-
-    /**
-     * Checks if no triangle is degenerated by checking the surface area being greater than zero.
-     * E.g. two points are the same or all three are collinear.
-     * @param polyhedron the polyhedron consisting of vertices and triangular faces
-     * @return true if triangles are fine and none of them is degenerate
-     */
-    bool checkTrianglesNotDegenerated(const Polyhedron &polyhedron);
-
-
-    namespace detail {
-
-        /**
-         * Calculates how often a vector starting at a specific origin intersects a polyhedron's mesh's triangles.
-         * @param rayOrigin the origin of the ray
-         * @param rayVector the vector describing the ray
-         * @param polyhedron the polyhedron consisting of vertices and triangular faces
-         * @return true if the ray intersects the triangle
-         */
-        size_t
-        countRayPolyhedronIntersections(const Array3 &rayOrigin, const Array3 &rayVector, const Polyhedron &polyhedron);
-
-        /**
-         * Calculates how often a vector starting at a specific origin intersects a triangular face.
-         * Uses the Möller–Trumbore intersection algorithm.
-         * @param rayOrigin the origin of the ray
-         * @param rayVector the vector describing the ray
-         * @param triangle a triangular face
-         * @return intersection point or null
-         *
-         * @related Adapted from https://en.wikipedia.org/wiki/Möller–Trumbore_intersection_algorithm
-         */
-        std::unique_ptr<Array3>
-        rayIntersectsTriangle(const Array3 &rayOrigin, const Array3 &rayVector, const Array3Triplet &triangle);
-
-    }
-
-}
diff --git a/src/polyhedralGravity/calculation/PolyhedronTransform.h b/src/polyhedralGravity/calculation/PolyhedronTransform.h
deleted file mode 100644
index f289272..0000000
--- a/src/polyhedralGravity/calculation/PolyhedronTransform.h
+++ /dev/null
@@ -1,29 +0,0 @@
-#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/input/DataSource.h b/src/polyhedralGravity/input/DataSource.h
index 78793ba..d01d8c0 100644
--- a/src/polyhedralGravity/input/DataSource.h
+++ b/src/polyhedralGravity/input/DataSource.h
@@ -1,5 +1,6 @@
 #pragma once
 
+#include <tuple>
 #include "polyhedralGravity/model/Polyhedron.h"
 
 namespace polyhedralGravity {
@@ -20,7 +21,7 @@ namespace polyhedralGravity {
          * Returns a Polyhedron from the underlying source.
          * @return a polyhedron
          */
-        virtual Polyhedron getPolyhedron() = 0;
+        virtual std::tuple<std::vector<Array3>, std::vector<IndexArray3>> getPolyhedralSource() = 0;
 
     };
 
diff --git a/src/polyhedralGravity/input/TetgenAdapter.cpp b/src/polyhedralGravity/input/TetgenAdapter.cpp
index 16d62db..cc8a003 100644
--- a/src/polyhedralGravity/input/TetgenAdapter.cpp
+++ b/src/polyhedralGravity/input/TetgenAdapter.cpp
@@ -2,7 +2,7 @@
 
 namespace polyhedralGravity {
 
-    Polyhedron TetgenAdapter::getPolyhedron() {
+    std::tuple<std::vector<Array3>, std::vector<IndexArray3>> TetgenAdapter::getPolyhedralSource() {
         //1. Step: Read in from files
         for (const auto &fileName: _fileNames) {
             size_t pos = fileName.find_last_of('.');
@@ -12,7 +12,7 @@ namespace polyhedralGravity {
         }
 
         //2. Convert tetgenio to Polyhedron
-        return {_vertices, _faces};
+        return std::make_tuple(_vertices, _faces);
     }
 
     void TetgenAdapter::readNode(const std::string &filename) {
diff --git a/src/polyhedralGravity/input/TetgenAdapter.h b/src/polyhedralGravity/input/TetgenAdapter.h
index 4443268..a088e61 100644
--- a/src/polyhedralGravity/input/TetgenAdapter.h
+++ b/src/polyhedralGravity/input/TetgenAdapter.h
@@ -22,7 +22,7 @@ namespace polyhedralGravity {
  *
  * The Adapter further keeps en eye on the already read in files in order to give feedback if data is in conflict.
  */
-    class TetgenAdapter : public DataSource {
+    class TetgenAdapter final : public DataSource {
 
         /**
          * The default exception message
@@ -72,7 +72,7 @@ namespace polyhedralGravity {
          * converted to a Polyhedron.
          * @return a Polyhedron
          */
-        Polyhedron getPolyhedron() override;
+        std::tuple<std::vector<Array3>, std::vector<IndexArray3>> getPolyhedralSource() override;
 
         /**
          * Reads nodes from a .node file
diff --git a/src/polyhedralGravity/input/YAMLConfigReader.cpp b/src/polyhedralGravity/input/YAMLConfigReader.cpp
index ed3b95c..827c080 100644
--- a/src/polyhedralGravity/input/YAMLConfigReader.cpp
+++ b/src/polyhedralGravity/input/YAMLConfigReader.cpp
@@ -39,7 +39,7 @@ namespace polyhedralGravity {
         if (_file[ROOT][INPUT] && _file[ROOT][INPUT][INPUT_CHECK]) {
             return _file[ROOT][INPUT][INPUT_CHECK].as<bool>();
         } else {
-            return false;
+            return true;
         }
     }
 
diff --git a/src/polyhedralGravity/input/YAMLConfigReader.h b/src/polyhedralGravity/input/YAMLConfigReader.h
index 4faa1dc..8120735 100644
--- a/src/polyhedralGravity/input/YAMLConfigReader.h
+++ b/src/polyhedralGravity/input/YAMLConfigReader.h
@@ -13,7 +13,7 @@ namespace polyhedralGravity {
      * The YAMLConfigReader serves as Interface between yaml-cpp and the Polyhedral Gravity Model and
      * reads in the input from an yaml configuration file.
      */
-    class YAMLConfigReader : public ConfigSource {
+    class YAMLConfigReader final : public ConfigSource {
 
         /*
          * The following static variables contain the names of the YAML nodes.
@@ -81,7 +81,7 @@ namespace polyhedralGravity {
 
         /**
          * Reads the enablement of the input sanity check from the yaml file.
-         * @return true if explicitly enabled, otherwise per-default false
+         * @return true or false if explicitly enabled, otherwise per-default true
          */
         bool getMeshInputCheckStatus() override;
 
diff --git a/src/polyhedralGravity/calculation/GravityEvaluable.cpp b/src/polyhedralGravity/model/GravityEvaluable.cpp
similarity index 94%
rename from src/polyhedralGravity/calculation/GravityEvaluable.cpp
rename to src/polyhedralGravity/model/GravityEvaluable.cpp
index 5b4e769..946ab7a 100644
--- a/src/polyhedralGravity/calculation/GravityEvaluable.cpp
+++ b/src/polyhedralGravity/model/GravityEvaluable.cpp
@@ -35,15 +35,13 @@ namespace polyhedralGravity {
         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);
+                            computationPoint[0], computationPoint[1], computationPoint[2], _polyhedron.getDensity());
         /*
          * 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());
+        const auto &[polyBegin, polyEnd] = _polyhedron.transformIterator(computationPoint);
+        auto zip1 = util::zip(polyBegin, _segmentVectors.begin(), _planeUnitNormals.begin(), _segmentUnitNormals.begin());
+        auto zip2 = util::zip(polyEnd, _segmentVectors.end(), _planeUnitNormals.end(), _segmentUnitNormals.end());
 
         SPDLOG_LOGGER_DEBUG(PolyhedralGravityLogger::DEFAULT_LOGGER.getLogger(),
                             "Starting to iterate over the planes...");
@@ -62,7 +60,7 @@ namespace polyhedralGravity {
                             "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;
+        const double prefix = util::GRAVITATIONAL_CONSTANT * _polyhedron.getDensity() * _polyhedron.getOrientationFactor();
 
         //10. Step: Final expressions after application of the prefix (and a division by 2 for the potential)
         potential = (potential * prefix) / 2.0;
@@ -94,6 +92,7 @@ namespace polyhedralGravity {
     }
 
     // Explicit template instantiation of the multipoint evaluate method
+
     template std::vector<GravityModelResult>
     GravityEvaluable::evaluate<true>(const std::vector<Array3> &computationPoints) const;
 
@@ -232,15 +231,15 @@ namespace polyhedralGravity {
     }
 
     std::string GravityEvaluable::toString() const {
-        std::stringstream ss;
-        ss << "<polyhedral_gravity.GravityEvaluable, density=" << _density << ", vertices= "
+        std::stringstream sstream;
+        sstream << "<polyhedral_gravity.GravityEvaluable, density=" << _polyhedron.getDensity() << ", vertices= "
            << _polyhedron.countVertices() << ", faces= " << _polyhedron.countFaces() << ">";
-        return ss.str();
+        return sstream.str();
     }
 
-    std::tuple<Polyhedron, double, std::vector<Array3Triplet>, std::vector<Array3>, std::vector<Array3Triplet>>
+    std::tuple<Polyhedron, std::vector<Array3Triplet>, std::vector<Array3>, std::vector<Array3Triplet>>
     GravityEvaluable::getState() const {
-        return std::make_tuple(_polyhedron, _density, _segmentVectors, _planeUnitNormals, _segmentUnitNormals);
+        return std::make_tuple(_polyhedron, _segmentVectors, _planeUnitNormals, _segmentUnitNormals);
     }
 
 }
diff --git a/src/polyhedralGravity/calculation/GravityEvaluable.h b/src/polyhedralGravity/model/GravityEvaluable.h
similarity index 61%
rename from src/polyhedralGravity/calculation/GravityEvaluable.h
rename to src/polyhedralGravity/model/GravityEvaluable.h
index 38335e5..b05e46e 100644
--- a/src/polyhedralGravity/calculation/GravityEvaluable.h
+++ b/src/polyhedralGravity/model/GravityEvaluable.h
@@ -3,15 +3,17 @@
 #include <tuple>
 #include <variant>
 #include <string>
+#include <optional>
 #include <sstream>
 
 #include "thrust/transform.h"
 #include "thrust/execution_policy.h"
 
-#include "polyhedralGravity/calculation/GravityModelDetail.h"
+#include "GravityModelDetail.h"
+#include "polyhedralGravity/util/UtilityContainer.h"
 #include "polyhedralGravity/input/TetgenAdapter.h"
-#include "polyhedralGravity/model/GravityModelData.h"
-#include "polyhedralGravity/model/Polyhedron.h"
+#include "GravityModelData.h"
+#include "Polyhedron.h"
 
 
 namespace polyhedralGravity {
@@ -27,45 +29,25 @@ namespace polyhedralGravity {
         /** The constant density polyhedron consisting of vertices and triangular faces */
         const Polyhedron _polyhedron;
 
-        /** The constant density of the polyhedron (the unit must match to the mesh, e.g., mesh in @f$[m]@f$ requires density in @f$[kg/m^3]@f$) */
-        const double _density;
-
         /** Cache for the segment vectors (segments between vertices of a polyhedral face) */
-        mutable std::vector<Array3Triplet> _segmentVectors;
+        mutable std::vector<Array3Triplet> _segmentVectors{};
 
         /** Cache for the plane unit normals (unit normals of the polyhedral faces) */
-        mutable std::vector<Array3> _planeUnitNormals;
+        mutable std::vector<Array3> _planeUnitNormals{};
 
         /** Cache for the segment unit normals (unit normals of each the polyhedral faces' segments) */
-        mutable std::vector<Array3Triplet> _segmentUnitNormals;
-
+        mutable std::vector<Array3Triplet> _segmentUnitNormals{};
 
     public:
-
         /**
          * Instantiates a GravityEvaluable with a given constant density polyhedron.
-         * @param polyhedron the polyhedron
-         * @param density the constant density (the unit must match to the mesh,
-         *                e.g., mesh in @f$[m]@f$ requires density in @f$[kg/m^3]@f$)
+         * In contrast to the {@link GravityModel::evaluate}, this evaluate method on the {@link GravityEvaluable}
+         * caches intermediate results and input data and subsequent evaluations will be faster.
          *
-         * @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.
+         * @param polyhedron the constant density polyhedron
          */
-        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 (the unit must match to the mesh,
-         *                e.g., mesh in @f$[m]@f$ requires density in @f$[kg/m^3]@f$)
-         */
-        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} {
+        explicit GravityEvaluable(const Polyhedron &polyhedron) :
+            _polyhedron{polyhedron} {
             this->prepare();
         }
 
@@ -73,20 +55,20 @@ namespace polyhedralGravity {
          * Instantiates a GravityEvaluable with a given constant density polyhedron and caches.
          * This is for restoring a GravityEvaluable from a previous state.
          * @param polyhedron the polyhedron
-         * @param density the constant density (the unit must match to the mesh,
-         *                e.g., mesh in @f$[m]@f$ requires density in @f$[kg/m^3]@f$)
          * @param segmentVectors the segment vectors
          * @param planeUnitNormals the plane unit normals
          * @param segmentUnitNormals the segment unit normals
          */
-        GravityEvaluable(const Polyhedron &polyhedron, double density, const std::vector<Array3Triplet> &segmentVectors,
+        GravityEvaluable(const Polyhedron &polyhedron,
+                         const std::vector<Array3Triplet> &segmentVectors,
                          const std::vector<Array3> &planeUnitNormals,
-                         const std::vector<Array3Triplet> &segmentUnitNormals) : _polyhedron{polyhedron},
-                                                                                 _density{density},
-                                                                                 _segmentVectors{segmentVectors},
-                                                                                 _planeUnitNormals{planeUnitNormals},
-                                                                                 _segmentUnitNormals{
-                                                                                         segmentUnitNormals} {}
+                         const std::vector<Array3Triplet> &segmentUnitNormals) :
+            _polyhedron{polyhedron},
+            _segmentVectors{segmentVectors},
+            _planeUnitNormals{planeUnitNormals},
+            _segmentUnitNormals{
+                    segmentUnitNormals} {
+        }
 
         /**
          * Evaluates the polyhedrale gravity model for a given constant density polyhedron at computation
@@ -117,15 +99,14 @@ namespace polyhedralGravity {
          * Returns a string representation of the GravityEvaluable.
          * @return string representation of the GravityEvaluable
          */
-        std::string toString() const;
+        [[nodiscard]] std::string toString() const;
 
         /**
          * Returns the polyhedron, the density and the internal caches.
          *
          * @return tuple of polyhedron, density, segmentVectors, planeUnitNormals and segmentUnitNormals
          */
-        std::tuple<Polyhedron, double, std::vector<Array3Triplet>, std::vector<Array3>, std::vector<Array3Triplet>>
-        getState() const;
+        std::tuple<Polyhedron, std::vector<Array3Triplet>, std::vector<Array3>, std::vector<Array3Triplet>> getState() const;
 
     private:
 
@@ -145,7 +126,7 @@ namespace polyhedralGravity {
         * at computation Point P
         */
         template<bool Parallelization = true>
-        GravityModelResult evaluate(const Array3 &computationPoint) const;
+        [[nodiscard]] GravityModelResult evaluate(const Array3 &computationPoint) const;
 
 
         /**
@@ -156,7 +137,7 @@ namespace polyhedralGravity {
          * @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;
+        [[nodiscard]] 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.
@@ -169,4 +150,4 @@ namespace polyhedralGravity {
 
     };
 
-} // namespace polyhedralGravity
\ No newline at end of file
+}// namespace polyhedralGravity
\ No newline at end of file
diff --git a/src/polyhedralGravity/model/GravityModel.cpp b/src/polyhedralGravity/model/GravityModel.cpp
new file mode 100644
index 0000000..545e585
--- /dev/null
+++ b/src/polyhedralGravity/model/GravityModel.cpp
@@ -0,0 +1,15 @@
+#include "GravityModel.h"
+
+namespace polyhedralGravity::GravityModel {
+
+    GravityModelResult evaluate(const Polyhedron &polyhedron, const Array3 &computationPoint, bool parallel) {
+        GravityEvaluable evaluable{polyhedron};
+        return std::get<GravityModelResult>(evaluable(computationPoint, parallel));
+    }
+
+    std::vector<GravityModelResult> evaluate(const Polyhedron &polyhedron, const std::vector<Array3> &computationPoints, bool parallel) {
+        GravityEvaluable evaluable{polyhedron};
+        return std::get<std::vector<GravityModelResult>>(evaluable(computationPoints, parallel));
+    }
+
+}
diff --git a/src/polyhedralGravity/calculation/GravityModel.h b/src/polyhedralGravity/model/GravityModel.h
similarity index 67%
rename from src/polyhedralGravity/calculation/GravityModel.h
rename to src/polyhedralGravity/model/GravityModel.h
index 5936380..df3fb20 100644
--- a/src/polyhedralGravity/calculation/GravityModel.h
+++ b/src/polyhedralGravity/model/GravityModel.h
@@ -1,12 +1,9 @@
 #pragma once
 
-#include <utility>
 #include <array>
 #include <vector>
-#include <algorithm>
-#include <variant>
 
-#include "polyhedralGravity/calculation/GravityEvaluable.h"
+#include "polyhedralGravity/model/GravityEvaluable.h"
 #include "polyhedralGravity/model/GravityModelData.h"
 #include "polyhedralGravity/model/Polyhedron.h"
 
@@ -22,27 +19,23 @@ namespace polyhedralGravity::GravityModel {
      * 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 (the unit must match to the mesh,
-     *                e.g., mesh in @f$[m]@f$ requires density in @f$[kg/m^3]@f$)
      * @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
      */
-    GravityModelResult evaluate(const PolyhedronOrSource &polyhedron, double density, const Array3 &computationPoint, bool parallel = true);
+    GravityModelResult evaluate(const Polyhedron &polyhedron, 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 (the unit must match to the mesh,
-     *                e.g., mesh in @f$[m]@f$ requires density in @f$[kg/m^3]@f$)
      * @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);
+    evaluate(const Polyhedron &polyhedron, const std::vector<Array3> &computationPoints, bool parallel = true);
 
-}
+}
\ No newline at end of file
diff --git a/src/polyhedralGravity/model/GravityModelData.h b/src/polyhedralGravity/model/GravityModelData.h
index 105d644..0f61aaa 100644
--- a/src/polyhedralGravity/model/GravityModelData.h
+++ b/src/polyhedralGravity/model/GravityModelData.h
@@ -49,6 +49,7 @@ namespace polyhedralGravity {
      */
     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.
diff --git a/src/polyhedralGravity/calculation/GravityModelDetail.cpp b/src/polyhedralGravity/model/GravityModelDetail.cpp
similarity index 100%
rename from src/polyhedralGravity/calculation/GravityModelDetail.cpp
rename to src/polyhedralGravity/model/GravityModelDetail.cpp
diff --git a/src/polyhedralGravity/calculation/GravityModelDetail.h b/src/polyhedralGravity/model/GravityModelDetail.h
similarity index 98%
rename from src/polyhedralGravity/calculation/GravityModelDetail.h
rename to src/polyhedralGravity/model/GravityModelDetail.h
index b7dc354..8336afd 100644
--- a/src/polyhedralGravity/calculation/GravityModelDetail.h
+++ b/src/polyhedralGravity/model/GravityModelDetail.h
@@ -15,10 +15,9 @@
 #include "thrust/execution_policy.h"
 #include "xsimd/xsimd.hpp"
 
-#include "polyhedralGravity/calculation/PolyhedronTransform.h"
+#include "Polyhedron.h"
+#include "GravityModelData.h"
 #include "polyhedralGravity/input/TetgenAdapter.h"
-#include "polyhedralGravity/model/Polyhedron.h"
-#include "polyhedralGravity/model/GravityModelData.h"
 #include "polyhedralGravity/util/UtilityConstants.h"
 #include "polyhedralGravity/util/UtilityContainer.h"
 #include "polyhedralGravity/util/UtilityThrust.h"
diff --git a/src/polyhedralGravity/model/Polyhedron.cpp b/src/polyhedralGravity/model/Polyhedron.cpp
new file mode 100644
index 0000000..a6780da
--- /dev/null
+++ b/src/polyhedralGravity/model/Polyhedron.cpp
@@ -0,0 +1,252 @@
+#include "Polyhedron.h"
+
+#include "polyhedralGravity/input/TetgenAdapter.h"
+
+
+namespace polyhedralGravity {
+
+    Polyhedron::Polyhedron(const std::vector<Array3> &vertices,
+            const std::vector<IndexArray3> &faces, double density, const NormalOrientation &orientation, const PolyhedronIntegrity &integrity)
+        : _vertices{vertices},
+          _faces{faces},
+          _density{density},
+          _orientation{orientation} {
+        //Checks that the node with index zero is actually used
+        if (_faces.end() == std::find_if(_faces.begin(), _faces.end(), [&](auto &face) {
+            return face[0] == 0 || face[1] == 0 || face[2] == 0;
+        })) {
+            throw std::invalid_argument("The node with index zero (0) was never used in any face! This is "
+                    "no valid polyhedron. Probable issue: Started numbering the vertices of "
+                    "the polyhedron at one (1).");
+        }
+        this->runIntegrityMeasures(integrity);
+    }
+
+    Polyhedron::Polyhedron(const PolyhedralSource &polyhedralSource, double density, const NormalOrientation &orientation, const PolyhedronIntegrity &integrity)
+        : Polyhedron{std::get<std::vector<Array3>>(polyhedralSource), std::get<std::vector<IndexArray3>>(polyhedralSource), density, orientation, integrity} {}
+
+    Polyhedron::Polyhedron(const PolyhedralFiles &polyhedralFiles, double density, const NormalOrientation &orientation, const PolyhedronIntegrity &integrity)
+        : Polyhedron{TetgenAdapter{polyhedralFiles}.getPolyhedralSource(), density, orientation, integrity} {}
+
+    Polyhedron::Polyhedron(const std::variant<PolyhedralSource, PolyhedralFiles> &polyhedralSource, double density, const NormalOrientation &orientation, const PolyhedronIntegrity &integrity)
+    : Polyhedron{std::holds_alternative<PolyhedralSource>(polyhedralSource) ? std::get<PolyhedralSource>(polyhedralSource) : TetgenAdapter{std::get<PolyhedralFiles>(polyhedralSource)}.getPolyhedralSource(),
+        density, orientation, integrity} {}
+
+    const std::vector<Array3> &Polyhedron::getVertices() const {
+        return _vertices;
+    }
+
+    const Array3 &Polyhedron::getVertex(size_t index) const {
+        return _vertices[index];
+    }
+
+    size_t Polyhedron::countVertices() const {
+        return _vertices.size();
+    }
+
+    const std::vector<IndexArray3> &Polyhedron::getFaces() const {
+        return _faces;
+    }
+
+    const IndexArray3 &Polyhedron::getFace(size_t index) const {
+        return _faces[index];
+    }
+
+    Array3Triplet Polyhedron::getResolvedFace(size_t index) const {
+        return {_vertices[_faces[index][0]], _vertices[_faces[index][1]], _vertices[_faces[index][2]]};
+    }
+
+    size_t Polyhedron::countFaces() const {
+        return _faces.size();
+    }
+
+    double Polyhedron::getDensity() const {
+        return _density;
+    }
+
+    void Polyhedron::setDensity(double density) {
+        _density = density;
+    }
+
+    NormalOrientation Polyhedron::getOrientation() const {
+        return _orientation;
+    }
+
+    double Polyhedron::getOrientationFactor() const {
+        return _orientation == NormalOrientation::OUTWARDS ? 1.0 : -1.0;
+    }
+
+    std::string Polyhedron::toString() const {
+        std::stringstream sstream{};
+        sstream << "<polyhedral_gravity.Polyhedron, density = " << _density << ", vertices = "
+           << countVertices() << ", faces = " << countFaces() << ", orientation = " << _orientation << ">";
+        return sstream.str();
+    }
+
+    std::tuple<std::vector<Array3>, std::vector<IndexArray3>, double, NormalOrientation> Polyhedron::getState() const {
+        return std::make_tuple(_vertices, _faces, _density, _orientation);
+    }
+
+    std::pair<NormalOrientation, std::set<size_t>> Polyhedron::checkPlaneUnitNormalOrientation() const {
+        // 1. Step: Find all indices of normals which vioate the constraint outwards pointing
+        const auto&[polyBegin, polyEnd] = this->transformIterator();
+        const size_t n = this->countFaces();
+        // Vector contains TRUE if the corrspeonding index VIOLATES the OUTWARDS cirteria
+        // Vector contains FALSE if the cooresponding index FULFILLS the OUTWARDS criteria
+        thrust::device_vector<bool> violatingBoolOutwards(n, false);
+        thrust::transform(
+            thrust::device,
+                polyBegin,
+                polyEnd,
+                violatingBoolOutwards.begin(),
+                [&](const auto &face) {
+                    // If the ray intersects the polyhedron odd number of times the normal points inwards
+                    // Hence, violating the OUTWARDS constraint
+                    const size_t intersects = this->countRayPolyhedronIntersections(face);
+                    return intersects % 2 != 0;
+                });
+        const size_t numberOfOutwardsViolations = std::count(violatingBoolOutwards.cbegin(), violatingBoolOutwards.cend(), true);
+        // 2. Step: Create a set with only the indices violating the constraint
+        std::set<size_t> violatingIndices{};
+        auto countingIterator = thrust::make_counting_iterator<size_t>(0);
+
+        if (numberOfOutwardsViolations <= n / 2) {
+            std::copy_if(countingIterator, countingIterator + n, std::inserter(violatingIndices, violatingIndices.end()), [&violatingBoolOutwards](const size_t &index) {
+                return violatingBoolOutwards[index];
+            });
+            // 3a. Step: Return the outwards pointing as major orientation
+            // and the violating faces, i.e. which have inwards pointing normals
+            return std::make_pair(NormalOrientation::OUTWARDS, violatingIndices);
+        } else {
+            std::copy_if(countingIterator, countingIterator + n, std::inserter(violatingIndices, violatingIndices.end()), [&violatingBoolOutwards](const size_t &index) {
+                return !violatingBoolOutwards[index];
+            });
+            // 3b. Step: Return the inwards pointing as major orientation and
+            // the violating faces, i.e. which have outwards pointing normals
+            return std::make_pair(NormalOrientation::INWARDS, violatingIndices);
+        }
+    }
+
+    void Polyhedron::runIntegrityMeasures(const PolyhedronIntegrity &integrity) {
+        using util::operator<<;
+        switch (integrity) {
+            case PolyhedronIntegrity::DISABLE:
+                return;
+            case PolyhedronIntegrity::AUTOMATIC:
+                SPDLOG_LOGGER_WARN(PolyhedralGravityLogger::DEFAULT_LOGGER.getLogger(),
+                                   "The mesh check is enabled and analyzes the polyhedron for degnerated faces & "
+                                   "that all plane unit normals point in the specified direction. This checks requires "
+                                   "a quadratic runtime cost which is most of the time not desirable. "
+                                   "Please explicitly set the integrity_check to either VERIFY, HEAL or DISABLE."
+                                   "You can find further details in the documentation!");
+            // NO BREAK! AUTOMATIC implies VERIFY, but with a info mesage to explcitly set the option
+            case PolyhedronIntegrity::VERIFY:
+            // NO BREAK! VERIFY terminates earlier, but does in the beginning the same as HEAL
+            case PolyhedronIntegrity::HEAL:
+                if (!this->checkTrianglesNotDegenerated()) {
+                    throw std::invalid_argument{"At least on triangle in the mesh is degenerated and its surface area equals zero!"};
+                }
+                const auto &[actualOrientation, violatingIndices] = this->checkPlaneUnitNormalOrientation();
+                if (actualOrientation != _orientation || !violatingIndices.empty()) {
+                    std::stringstream sstream{};
+                    sstream << "The plane unit normals are not all pointing in the specified direction " << _orientation << '\n';
+                    if (violatingIndices.empty()) {
+                        sstream << "Instead all plane unit normals are pointing "
+                                << actualOrientation
+                                << ". You can either reconstruct the polyhedron with the orientation set to " << actualOrientation
+                                << ". Alternativly, you can reconstruct with the inetgrity_check set to HEAL";
+                    } else {
+                        sstream << "The actual majority orientation of the polyhedron's normals is " << actualOrientation
+                                << ". You can either:\n 1) Fix the ordering of the following faces:\n"
+                                << violatingIndices << '\n'
+                                << "2) Or you reconstruct the polyhedron using the integrity_check set to HEAL.";
+                    }
+                    // In case of HEAL, don't throw but repair
+                    if (integrity != PolyhedronIntegrity::HEAL) {
+                        throw std::invalid_argument(sstream.str());
+                    } else {
+                        this->healPlaneUnitNormalOrientation(actualOrientation, violatingIndices);
+                    }
+                }
+        }
+    }
+
+    bool Polyhedron::checkTrianglesNotDegenerated() const {
+        const auto &[begin, end] = this->transformIterator();
+        // All triangles surface area needs to be greater than zero
+        return thrust::transform_reduce(
+                thrust::device,
+                begin, end, [](const Array3Triplet &face) {
+                    return util::surfaceArea(face) > 0.0;
+                }, true, thrust::logical_and<bool>());
+    }
+
+    void Polyhedron::healPlaneUnitNormalOrientation(const NormalOrientation &actualOrientation, const std::set<size_t> &violatingIndices) {
+        // Assign the majority plane unit normal orientation
+        _orientation = actualOrientation;
+        // Fix the vioalting faces by exchaning the vertex ordering (exchaning index 0 with index 1 in the face)
+        std::for_each(violatingIndices.cbegin(), violatingIndices.cend(), [this](size_t i) {
+            std::swap(this->_faces[i][0], this->_faces[i][1]);
+        });
+    }
+
+    size_t Polyhedron::countRayPolyhedronIntersections(const Array3Triplet &face) const {
+        using namespace util;
+        // The centroid of the triangular face
+        const Array3 centroid = (face[0] + face[1] + face[2]) / 3.0;
+
+        // The normal of the plane calculated with two segments of the triangle
+        // The normal is the rayVector starting at the rayOrigin
+        const Array3 segmentVector1 = face[1] - face[0];
+        const Array3 segmentVector2 = face[2] - face[1];
+        const Array3 rayVector = normal(segmentVector1, segmentVector2);
+
+        // The origin of the array has a slight offset in direction of the normal
+        const Array3 rayOrigin = centroid + (rayVector * EPSILON);
+
+        // Count every triangular face which is intersected by the ray
+        const auto &[begin, end] = this->transformIterator();
+        std::set<Array3> intersections{};
+        std::for_each(begin, end, [&rayOrigin, &rayVector, &intersections](const Array3Triplet &otherFace) {
+            const std::unique_ptr<Array3> intersection = rayIntersectsTriangle(rayOrigin, rayVector, otherFace);
+            if (intersection != nullptr) {
+                intersections.insert(*intersection);
+            }
+        });
+        return intersections.size();
+    }
+
+    std::unique_ptr<Array3> Polyhedron::rayIntersectsTriangle(const Array3 &rayOrigin, const Array3 &rayVector, const Array3Triplet &triangle) {
+        // Adapted Möller–Trumbore intersection algorithm
+        // see https://en.wikipedia.org/wiki/Möller–Trumbore_intersection_algorithm
+        using namespace util;
+        const Array3 edge1 = triangle[1] - triangle[0];
+        const Array3 edge2 = triangle[2] - triangle[0];
+        const Array3 h = cross(rayVector, edge2);
+        const double a = dot(edge1, h);
+        if (a > -EPSILON && a < EPSILON) {
+            return nullptr;
+        }
+
+        const double f = 1.0 / a;
+        const Array3 s = rayOrigin - triangle[0];
+        const double u = f * dot(s, h);
+        if (u < 0.0 || u > 1.0) {
+            return nullptr;
+        }
+
+        const Array3 q = cross(s, edge1);
+        const double v = f * dot(rayVector, q);
+        if (v < 0.0 || u + v > 1.0) {
+            return nullptr;
+        }
+
+        const double t = f * dot(edge2, q);
+        if (t > EPSILON) {
+            return std::make_unique<Array3>(rayOrigin + rayVector * t);
+        } else {
+            return nullptr;
+        }
+    }
+
+}
\ No newline at end of file
diff --git a/src/polyhedralGravity/model/Polyhedron.h b/src/polyhedralGravity/model/Polyhedron.h
index e9d21cf..c80ed30 100644
--- a/src/polyhedralGravity/model/Polyhedron.h
+++ b/src/polyhedralGravity/model/Polyhedron.h
@@ -3,30 +3,106 @@
 #include <utility>
 #include <vector>
 #include <array>
+#include <set>
 #include <algorithm>
 #include <exception>
 #include <stdexcept>
 #include <tuple>
 #include <variant>
 #include <string>
+#include <sstream>
+#include <memory>
+#include "thrust/copy.h"
+#include "thrust/device_vector.h"
+#include "polyhedralGravity/output/Logging.h"
+#include "thrust/transform_reduce.h"
+#include "thrust/execution_policy.h"
+#include "thrust/iterator/counting_iterator.h"
 #include "polyhedralGravity/model/GravityModelData.h"
+#include "thrust/iterator/transform_iterator.h"
+#include "polyhedralGravity/util/UtilityContainer.h"
+#include "polyhedralGravity/util/UtilityConstants.h"
 
 namespace polyhedralGravity {
 
     /* Forward declaration of Polyhedron */
     class Polyhedron;
 
+    /** A polyhedron defined by a set of filenames, as available in {@link TetgenAdapter} */
+    using PolyhedralFiles = std::vector<std::string>;
+
+    /** A polyhedron defined by a set of vertices and face indices */
+    using PolyhedralSource = std::tuple<std::vector<Array3>, std::vector<IndexArray3>>;
+
+    /**
+     * The orientation of the plane unit normals of the polyhedron.
+     * We use this property as the concret definition of the vertices ordering depends on the
+     * utilized cooridnate system.
+     * However, the normal alignement is independent. Tsoulis et al. equations require the
+     * normals to point outwards of the polyhedron. If the opposite hold, the result is
+     * negated.
+     */
+    enum class NormalOrientation: char {
+        /** Outwards pointing plane unit normals */
+        OUTWARDS,
+        /** Inwards pointing plane unit normals */
+        INWARDS
+    };
+
+
     /**
-     * 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
+     * Overloaded insertion operator to output the string representation of a NormalOrientation enum value.
+     *
+     * @param os The output stream to write the string representation to.
+     * @param orientation The NormalOrientation enum value to output.
+     * @return The output stream after writing the string representation.
      */
-    using PolyhedralSource = std::variant<std::tuple<std::vector<Array3>, std::vector<IndexArray3>>, std::vector<std::string>>;
+    inline std::ostream &operator<<(std::ostream &os, const NormalOrientation &orientation) {
+        switch (orientation) {
+            case NormalOrientation::OUTWARDS:
+                os << "OUTWARDS";
+                break;
+            case NormalOrientation::INWARDS:
+                os << "INWARDS";
+                break;
+            default:
+                os << "Unknown";
+                break;
+        }
+        return os;
+    }
 
     /**
-     * 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
+     * The three mode the poylhedron class takes in
+     * the constructor in order to determine what initilaization checks to conduct.
+     * This enum is exclusivly utilized in the constrcutor of a {@link Polyhedron} and its private method
+     * {@link runIntegrityMeasures}
      */
-    using PolyhedronOrSource = std::variant<Polyhedron, std::tuple<std::vector<Array3>, std::vector<IndexArray3>>, std::vector<std::string>>;
+    enum class PolyhedronIntegrity: char {
+        /**
+         * All activities regarding MeshChecking are disabled.
+         * No runtime overhead!
+         */
+        DISABLE,
+        /**
+         * Only verification of the NormalOrientation.
+         * A misalignment (e.g. specified OUTWARDS, but is not) leads to a runtime_error.
+         * Runtime Cost O(n^2)
+         */
+        VERIFY,
+        /**
+         * Like VERIFY, but also informs the user about the option in any case on the runtime costs.
+         * This is the implicit default option.
+         * Runtime Cost: O(n^2) and output to stdout in every case!
+         */
+        AUTOMATIC,
+        /**
+         * Verification and Autmatioc Healing of the NormalOrientation.
+         * A misalignemt does not lead to a runtime_error, but to an internal correction.
+         * Runtime Cost: O(n^2) and a modification of the mesh input!
+         */
+        HEAL,
+    };
 
     /**
      * Data structure containing the model data of one polyhedron. This includes nodes, edges (faces) and elements.
@@ -37,6 +113,8 @@ namespace polyhedralGravity {
         /**
          * A vector containing the vertices of the polyhedron.
          * Each node is an array of size three containing the xyz coordinates.
+         * The mesh must be scaled in the same units as the density is given
+         * (the unit must match to the mesh, e.g., mesh in @f$[m]@f$ requires density in @f$[kg/m^3]@f$)
          */
         const std::vector<Array3> _vertices;
 
@@ -47,49 +125,89 @@ namespace polyhedralGravity {
          * two nodes.
          * @example face consisting of {1, 2, 3} --> segments: {1, 2}, {2, 3}, {3, 1}
          */
-        const std::vector<IndexArray3> _faces;
+        std::vector<IndexArray3> _faces;
 
+        /** The constant density of the polyhedron (the unit must match to the mesh, e.g., mesh in @f$[m]@f$ requires density in @f$[kg/m^3]@f$) */
+        double _density;
 
-    public:
+        /**
+         * An optional variable representing the orientation of the polyhedron.
+         * This variable holds an instance of the NormalOrientation enumeration class,
+         * which indicates whether the plane unit normals are pointing outwards or inwards.
+         * To set the orientation of the polyhedron, assign a valid `NormalOrientation` value to this variable.
+         * If the value is not set, the orientation is considered to be unspecified.
+         */
+        NormalOrientation _orientation;
 
+    public:
         /**
-         * Generates an empty polyhedron.
+         * Generates a polyhedron from nodes and faces.
+         * @param vertices a vector of nodes
+         * @param faces a vector of faces containing the formation of faces off vertices
+         * @param density the density of the polyhedron (it must match the unit of the mesh, e.g., mesh in @f$[m]@f$ requires density in @f$[kg/m^3]@f$)
+         * @param orientation specify if the plane unit normals point outwards or inwards (defaults: to OUTWARDS)
+         * @param integrity specify if the mesh input is checked/ healed to fulfill the constraints of Tsoulis' algorithm (see {@link PolyhedronIntegrity})
+         *
+         * @note ASSERTS PRE-CONDITION that the in the indexing in the faces vector starts with zero!
+         * @throws std::invalid_argument if no face contains the node zero indicating mathematical index
+         * @throws std::invalid_argument dpending on the {@param integrity} flag
          */
-        Polyhedron()
-                : _vertices{},
-                  _faces{} {}
+        Polyhedron(
+                const std::vector<Array3> &vertices,
+                const std::vector<IndexArray3> &faces,
+                double density,
+                const NormalOrientation &orientation = NormalOrientation::OUTWARDS,
+                const PolyhedronIntegrity &integrity = PolyhedronIntegrity::AUTOMATIC
+                );
 
         /**
          * Generates a polyhedron from nodes and faces.
-         * @param nodes vector containing the nodes
-         * @param faces vector containing the triangle faces.
+         * @param polyhedralSource a tuple of vector containing the nodes and trianglular faces.
+         * @param density the density of the polyhedron (it must match the unit of the mesh, e.g., mesh in @f$[m]@f$ requires density in @f$[kg/m^3]@f$)
+         * @param orientation specify if the plane unit normals point outwards or inwards (defaults: to OUTWARDS)
+         * @param integrity specify if the mesh input is checked/ healed to fulfill the constraints of Tsoulis' algorithm (see {@link PolyhedronIntegrity})
          *
          * @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<Array3> nodes, std::vector<IndexArray3> faces)
-                : _vertices{std::move(nodes)},
-                  _faces{std::move(faces)} {
-            //Checks that the node with index zero is actually used
-            if (_faces.end() == std::find_if(_faces.begin(), _faces.end(), [&](auto &face) {
-                return face[0] == 0 || face[1] == 0 || face[2] == 0;
-            })) {
-                throw std::runtime_error("The node with index zero (0) was never used in any face! This is "
-                                         "no valid polyhedron. Probable issue: Started numbering the vertices of "
-                                         "the polyhedron at one (1).");
-            }
-        }
+         * @throws std::invalid_argument if no face contains the node zero indicating mathematical index
+         * @throws std::invalid_argument dpending on the {@param integrity} flag
+         */
+        Polyhedron(
+                const PolyhedralSource &polyhedralSource,
+                double density,
+                const NormalOrientation &orientation = NormalOrientation::OUTWARDS,
+                const PolyhedronIntegrity &integrity = PolyhedronIntegrity::AUTOMATIC
+                );
 
         /**
-         * Creates a polyhedron from a tuple of nodes and faces.
-         * @param data tuple of nodes and faces
+         * Generates a polyhedron from nodes and faces.
+         * @param polyhedralFiles a list of files (see {@link TetgenAdapter}
+         * @param density the density of the polyhedron (it must match the unit of the mesh, e.g., mesh in @f$[m]@f$ requires density in @f$[kg/m^3]@f$)
+         * @param orientation specify if the plane unit normals point outwards or inwards (defaults: to OUTWARDS)
+         * @param integrity specify if the mesh input is checked/ healed to fulfill the constraints of Tsoulis' algorithm (see {@link PolyhedronIntegrity})
          *
          * @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
+         * @throws std::invalid_argument if no face contains the node zero indicating mathematical index
+         * @throws std::invalid_argument dpending on the {@param integrity} flag
          */
-        explicit Polyhedron(std::tuple<std::vector<Array3>, std::vector<IndexArray3>> data)
-                : Polyhedron(std::get<std::vector<Array3>>(data), std::get<std::vector<IndexArray3>>(data)) {}
+        Polyhedron(const PolyhedralFiles &polyhedralFiles, double density,
+                   const NormalOrientation &orientation = NormalOrientation::OUTWARDS,
+                   const PolyhedronIntegrity &integrity = PolyhedronIntegrity::AUTOMATIC);
 
+        /**
+         * Generates a polyhedron from nodes and faces.
+         * This constructor using a variant is maninly utilized from the Python Interface.
+         * @param polyhedralSource a list of files (see {@link TetgenAdapter} or a tuple of vector containing the nodes and trianglular faces.
+         * @param density the density of the polyhedron (it must match the unit of the mesh, e.g., mesh in @f$[m]@f$ requires density in @f$[kg/m^3]@f$)
+         * @param orientation specify if the plane unit normals point outwards or inwards (defaults: to OUTWARDS)
+         * @param integrity specify if the mesh input is checked/ healed to fulfill the constraints of Tsoulis' algorithm (see {@link PolyhedronIntegrity})
+         *
+         * @note ASSERTS PRE-CONDITION that the in the indexing in the faces vector starts with zero!
+         * @throws std::invalid_argument if no face contains the node zero indicating mathematical index
+         * @throws std::invalid_argument dpending on the {@param integrity} flag
+         */
+        Polyhedron(const std::variant<PolyhedralSource, PolyhedralFiles> &polyhedralSource, double density,
+                   const NormalOrientation &orientation = NormalOrientation::OUTWARDS,
+                   const PolyhedronIntegrity &integrity = PolyhedronIntegrity::AUTOMATIC);
 
         /**
          * Default destructor
@@ -100,52 +218,161 @@ namespace polyhedralGravity {
          * Returns the vertices of this polyhedron
          * @return vector of cartesian coordinates
          */
-        [[nodiscard]] const std::vector<Array3> &getVertices() const {
-            return _vertices;
-        }
+        [[nodiscard]] const std::vector<Array3> &getVertices() const;
 
         /**
          * Returns the vertex at a specific index
          * @param index size_t
          * @return cartesian coordinates of the vertex at index
          */
-        [[nodiscard]] const Array3 &getVertex(size_t index) const {
-            return _vertices[index];
-        }
+        [[nodiscard]] const Array3 &getVertex(size_t index) const;
 
         /**
          * The number of points (nodes) that make up the polyhedron.
          * @return a size_t
          */
-        [[nodiscard]] size_t countVertices() const {
-            return _vertices.size();
-        }
+        [[nodiscard]] size_t countVertices() const;
 
         /**
-         * Returns the number of faces (triangles) that make up the polyhedral.
-         * @return a size_t
+         * 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]] size_t countFaces() const {
-            return _faces.size();
-        }
+        [[nodiscard]] const std::vector<IndexArray3> &getFaces() const;
 
         /**
-         * Returns the triangular faces of this polyhedron
-         * @return vector of triangular faces, where each element size_t references a vertex in the vertices vector
+         * Returns the indices of the vertices making up the face at the given index.
+         * @param index size_t
+         * @return triplet of the vertic'es indices forming the face
          */
-        [[nodiscard]] const std::vector<IndexArray3> &getFaces() const {
-            return _faces;
-        }
+        [[nodiscard]] const IndexArray3 &getFace(size_t index) const;
 
         /**
          * 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]]};
+        [[nodiscard]] Array3Triplet getResolvedFace(size_t index) const;
+
+        /**
+         * Returns the number of faces (triangles) that make up the polyhedral.
+         * @return a size_t
+         */
+        [[nodiscard]] size_t countFaces() const;
+
+        /**
+         * Returns the constant density of this polyhedron.
+         * @return the constant density a double
+         */
+        [[nodiscard]] double getDensity() const;
+
+        /**
+         * Sets the density to a new value. The density's unit must match to the scaling of the mesh.
+         * @param density the new constant density of the polyhedron
+         */
+        void setDensity(double density);
+
+        /**
+         * Returns the orientation of the plane unit normals of this polyhedron
+         * @return OUTWARDS or INWARDS
+         */
+        [[nodiscard]] NormalOrientation getOrientation() const;
+
+        /**
+         * Retruns the plane unit normal orientation factor.
+         * If the unit normals are outwards pointing, it is 1.0 as Tsoulis inteneded.
+         * If the unit normals are inwards pointing, it is -1.0 (reversed).
+         * @return 1.0 or -1.0 depending on plane unit orientation
+         */
+        [[nodiscard]] double getOrientationFactor() const;
+
+        /**
+         * Returns a string representation of the Polyhedron.
+         * Mainly used for the representation method in the Python interface.
+         *
+         * @return string representation of the Polyhedron
+         */
+        [[nodiscard]] std::string toString() const;
+
+        /**
+         * Returns the internal data strcuture of Python pickle support.
+         * @return tuple of vertices, faces, density and normal orientation
+         */
+        [[nodiscard]] std::tuple<std::vector<Array3>, std::vector<IndexArray3>, double, NormalOrientation> getState() const;
+
+        /**
+         * 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 offset the offset to apply
+         * @return pair of transform iterators
+         */
+        [[nodiscard]] inline auto transformIterator(const Array3 &offset = {0.0, 0.0, 0.0}) const {
+            //The offset must be captured by value to ensure its lifetime!
+            const auto lambdaOffsetApplication = [this, offset](const IndexArray3 &face) -> Array3Triplet {
+                using namespace util;
+                return {this->getVertex(face[0]) - offset, this->getVertex(face[1]) - offset,
+                        this->getVertex(face[2]) - offset};
+            };
+            auto first = thrust::make_transform_iterator(this->getFaces().begin(), lambdaOffsetApplication);
+            auto last = thrust::make_transform_iterator(this->getFaces().end(), lambdaOffsetApplication);
+            return std::make_pair(first, last);
         }
 
+        /**
+         * This method determines the majority vertex ordering of a polyhedron and the set of faces which
+         * violate the majority constraint and need to be adpated.
+         * Hence, if the set is empty, all faces obey to the returned ordering/ plane unit normal orientation.
+         *
+         * @return a pair consisting of majority ordering (OUTWARDS or INWARDS pointing normals)
+         *  and a set of face indices which violate the constraint
+         */
+        [[nodiscard]] std::pair<NormalOrientation, std::set<size_t>> checkPlaneUnitNormalOrientation() const;
+
+    private:
+        /**
+         * Checks the integrity of the polyhedron depending on the integrity flag.
+         *
+         * @param integrity the behavior depends on the value, see {@link PolyhedronIntegrity}
+         *
+         * @throws std::invalid_argument dpending on the {@param integrity} flag
+         */
+        void runIntegrityMeasures(const PolyhedronIntegrity &integrity);
+
+
+        /**
+         * Checks if no triangle is degenerated by checking the surface area being greater than zero.
+         * E.g. two points are the same or all three are collinear.
+         * @return true if triangles are fine and none of them is degenerate
+         */
+        [[nodiscard]] bool checkTrianglesNotDegenerated() const;
+
+        /**
+         * Fixes the orientation of the plane unit normals for a given set of violating face indices.
+         *
+         * @param actualOrientation The desired plane unit normal orientation.
+         * @param violatingIndices A set of indices representing the faces with violating orientations.
+         */
+        void healPlaneUnitNormalOrientation(const NormalOrientation &actualOrientation, const std::set<size_t> &violatingIndices);
+
+        /**
+         * Calculates how often a vector starting at a specific origin intersects a polyhedron's mesh's triangles.
+         * @param face the vector describing the ray
+         * @return true if the ray intersects the triangle
+         */
+        [[nodiscard]] size_t countRayPolyhedronIntersections(const Array3Triplet &face) const;
+
+        /**
+         * Calculates how often a vector starting at a specific origin intersects a triangular face.
+         * Uses the Möller–Trumbore intersection algorithm.
+         * @param rayOrigin the origin of the ray
+         * @param rayVector the vector describing the ray
+         * @param triangle a triangular face
+         * @return intersection point or null
+         *
+         * @related Adapted from https://en.wikipedia.org/wiki/Möller–Trumbore_intersection_algorithm
+         */
+        static std::unique_ptr<Array3> rayIntersectsTriangle(const Array3 &rayOrigin, const Array3 &rayVector, const Array3Triplet &triangle);
+
+
     };
 
 }
diff --git a/src/polyhedralGravity/util/UtilityContainer.h b/src/polyhedralGravity/util/UtilityContainer.h
index 9468dd9..8554eb3 100644
--- a/src/polyhedralGravity/util/UtilityContainer.h
+++ b/src/polyhedralGravity/util/UtilityContainer.h
@@ -1,6 +1,7 @@
 #pragma once
 
 #include <array>
+#include <set>
 #include <numeric>
 #include <utility>
 #include <algorithm>
@@ -15,7 +16,8 @@ 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
@@ -131,10 +133,10 @@ namespace polyhedralGravity::util {
     * @return a Container
      * TODO This method causes issues with the MVSC 19.31.31107.0? Although it is never used...
     */
-//    template<typename Container, typename Scalar>
-//    Container operator-(const Container &lhs, const Scalar &scalar) {
-//        return applyBinaryFunction(lhs, scalar, std::minus<>());
-//    }
+    //    template<typename Container, typename Scalar>
+    //    Container operator-(const Container &lhs, const Scalar &scalar) {
+    //        return applyBinaryFunction(lhs, scalar, std::minus<>());
+    //    }
 
     /**
     * Applies the Operation to a Container and a Scalar.
@@ -327,6 +329,34 @@ namespace polyhedralGravity::util {
         return std::abs(x - y) > maxExponentDifference;
     }
 
+    /**
+    * A utility struct that creates an overload set out of all the function objects it inherits from.
+    * It allows a lambda function to be used with std::visit in a variant.
+    * The lambda function needs to be able to handle all types contained in the variant,
+    * and this struct allows it to do so by inheriting the function-call operator from each type.
+    * @tparam Ts A template parameter pack representing all types for which the struct should be able to handle.
+    *
+    * @ref https://en.cppreference.com/w/cpp/utility/variant/visit
+    */
+    template<class... Ts>
+    struct overloaded : Ts... {
+        using Ts::operator()...;
+    };
+
+    /**
+    * This function template provides deduction guides for the overloaded struct.
+    * It deduces the template arguments for overloaded based on its constructor arguments
+    * thus saving you from explicitly mentioning them while instantiation.
+    * @tparam Ts A template parameter pack representing all types that will be passed to the overloaded struct.
+    * @param Ts Variable numbers of parameters to pass to the overloaded struct.
+    * @return This doesn't return a value, but rather assists in the creation of an overloaded object.
+    *          The type of the object will be overloaded<Ts...>.
+    *
+    * @ref https://en.cppreference.com/w/cpp/utility/variant/visit
+    */
+    template<class... Ts>
+    overloaded(Ts...) -> overloaded<Ts...>;
+
     namespace detail {
 
         /**
@@ -366,14 +396,30 @@ namespace polyhedralGravity::util {
      */
     template<typename T, size_t N>
     std::ostream &operator<<(std::ostream &os, const std::array<T, N> &array) {
-        os << "[";
-        auto it = array.begin();
-        auto end = array.end() - 1;
-        while (it != end) {
-            os << *it << " ";
-            ++it;
-        }
-        os << *it << "]";
+        os.operator<<('[');
+        os.operator<<(' ');
+        std::for_each(array.cbegin(), array.cend(), [&os](const auto& arg) {
+            os << arg << ' ';
+        });
+        os.operator<<(']');
+        return os;
+    }
+
+    /**
+     * Operator << for a set.
+     * @tparam T type of the set, must have an << operator overload
+     * @param os the ostream
+     * @param set the set itself
+     * @return ostream
+     */
+    template<typename T>
+    std::ostream &operator<<(std::ostream &os, const std::set<T> &set) {
+        os.operator<<('[');
+        os.operator<<(' ');
+        std::for_each(set.cbegin(), set.cend(), [&os](const auto& arg) {
+            os << arg << ' ';
+        });
+        os.operator<<(']');
         return os;
     }
 
diff --git a/src/polyhedralGravityPython/PolyhedralGravityPython.cpp b/src/polyhedralGravityPython/PolyhedralGravityPython.cpp
index c7c6600..75a785c 100644
--- a/src/polyhedralGravityPython/PolyhedralGravityPython.cpp
+++ b/src/polyhedralGravityPython/PolyhedralGravityPython.cpp
@@ -8,154 +8,299 @@
 
 #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"
+#include "polyhedralGravity/model/GravityModel.h"
+#include "polyhedralGravity/model/GravityEvaluable.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.doc() = R"mydelimiter(
+    The evaluation of the polyhedral gravity model requires the following parameters:
 
-    m.def("evaluate", [](const PolyhedralSource &polyhedralSource, double density,
-                         const std::variant<Array3, std::vector<Array3>> &computationPoints,
-                         bool parallel) -> std::variant<GravityModelResult, std::vector<GravityModelResult>> {
-              using namespace polyhedralGravity;
-              if (std::holds_alternative<Array3>(computationPoints)) {
-                  return std::visit([&density, &computationPoints, parallel](const auto &source) {
-                      using namespace polyhedralGravity;
-                      return GravityModel::evaluate(source, density, std::get<Array3>(computationPoints), parallel);
-                  }, polyhedralSource);
-              } else {
-                  return std::visit([&density, &computationPoints, parallel](const auto &source) {
-                      using namespace polyhedralGravity;
-                      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.
+    +------------------------------------------------------------------------------+
+    | Name                                                                         |
+    +==============================================================================+
+    | Polyhedral Mesh (either as vertices & faces or as polyhedral source files)   |
+    +------------------------------------------------------------------------------+
+    | Constant Density :math:`\rho`                                                |
+    +------------------------------------------------------------------------------+
+
+    In the Python Interface, you define these parameters as :py:class:`polyhedral_gravity.Polyhedron`
+
+    .. code-block:: python
+
+        from polyhedral_gravity import Polyhedron, GravityEvaluable, evaluate, PolyhedronIntegrity, NormalOrientation
+
+        polyhedron = Polyhedron(
+            polyhedral_source=(vertices, faces),           # (N,3) and (M,3) array-like
+            density=density,                               # Density of the Polyhedron, Unit must match to the mesh's scale
+            normal_orientation=NormalOrientation.OUTWARDS, # Possible values OUTWARDS (default) or INWARDS
+            integrity_check=PolyhedronIntegrity.VERIFY,    # Possible values AUTOMATIC (default), VERIFY, HEAL, DISABLE
+        )
+
+    .. note::
+
+        *Tsoulis et al.*'s formulation requires that the normals point :code:`OUTWARDS`.
+        The implementation **can handle both cases and also can automatically determine the property** if initiall set wrong.
+        Using :code:`AUTOMATIC` (default for first-time-user) or :code:`VERIFY` raises a :code:`ValueError` if the :py:class:`polyhedral_gravity.NormalOrientation` is wrong.
+        Using :code:`HEAL` will re-order the vertex sorting to fix errors.
+        Using :code:`DISABLE` will turn this check off and avoid :math:`O(n^2)` runtime complexcity of this check! Highly recommened, when you "know your mesh"!
+
+
+
+    The polyhedron's mesh's units must match with the constant density!
+    For example, if the mesh is in :math:`[m]`, then the constant density should be in :math:`[\frac{kg}{m^3}]`.
+
+    Afterwards one can use :py:func:`polyhedral_gravity.evaluate` the gravity at a single point *P* via:
+
+    .. code-block:: python
+
+        potential, acceleration, tensor = evaluate(
+            polyhedron=polyhedron,
+            computation_points=P,
+            parallel=True,
+        )
+
+    or via use the cached approach :py:class:`polyhedral_gravity.GravityEvaluable` (desriable for subsequent evaluations using the same :py:class:`polyhedral_gravity.Polyhedron`)
+
+    .. code-block:: python
+
+        evaluable = GravityEvaluable(polyhedron=polyhedron)
+        potential, acceleration, tensor = evaluable(
+            computation_points=P,
+            parallel=True,
+        )
+
+    .. note::
+
+        If :code:`P` would be an array of points, the return value would be a :code:`List[Tuple[potential, acceleration, tensor]]`!
+
+    The calculation outputs the following parameters for every Computation Point *P*.
+    The units of the respective output depend on the units of the input parameters (mesh and density)!
+
+    Hence, if e.g. your mesh is in :math:`km`, the density must match. Further, the output units will match the input units.
+
+    +------------------------------------------------------------------------------------------------+----------------------------------------------------------------------------+-----------------------------------------------------------------+
+    |         Name                                                                                   | If mesh :math:`[m]` and density :math:`[\frac{kg}{m^3}]`                   |                             Comment                             |
+    +================================================================================================+============================================================================+=================================================================+
+    |         :math:`V`                                                                              |  :math:`\frac{m^2}{s^2}` or :math:`\frac{J}{kg}`                           |           The potential or also called specific energy          |
+    +------------------------------------------------------------------------------------------------+----------------------------------------------------------------------------+-----------------------------------------------------------------+
+    |     :math:`V_x`, :math:`V_y`, :math:`V_z`                                                      |   :math:`\frac{m}{s^2}`                                                    |The gravitational acceleration in the three cartesian directions |
+    +------------------------------------------------------------------------------------------------+----------------------------------------------------------------------------+-----------------------------------------------------------------+
+    | :math:`V_{xx}`, :math:`V_{yy}`, :math:`V_{zz}`, :math:`V_{xy}`, :math:`V_{xz}`, :math:`V_{yz}` |   :math:`\frac{1}{s^2}`                                                    |The spatial rate of change of the gravitational acceleration     |
+    +------------------------------------------------------------------------------------------------+----------------------------------------------------------------------------+-----------------------------------------------------------------+
+
+    This model's output obeys to the geodesy and geophysics sign conventions.
+    Hence, the potential :math:`V` for a polyhedron with a mass :math:`m > 0` is defined as **positive**.
+
+    The accelerations :math:`V_x`, :math:`V_y`, :math:`V_z` are defined as
+
+    .. math::
+
+        \textbf{g} = + \nabla V = \left( \frac{\partial V}{\partial x}, \frac{\partial V}{\partial y}, \frac{\partial V}{\partial z} \right)
+
+    Accordingly, the second derivative tensor is defined as the derivative of :math:`\textbf{g}`.
+    )mydelimiter";
+
+    py::enum_<NormalOrientation>(m, "NormalOrientation", R"mydelimiter(
+        The orientation of the plane unit normals of the polyhedron.
+        *Tsoulis et al.* equations require the normals to point outwards of the polyhedron.
+        If the opposite hold, the result is negated.
+        The implementation can handle both cases.
+        )mydelimiter")
+        .value("OUTWARDS", NormalOrientation::OUTWARDS, "Outwards pointing plane unit normals")
+        .value("INWARDS", NormalOrientation::INWARDS, "Inwards pointing plane unit normals");
+
+    py::enum_<PolyhedronIntegrity>(m, "PolyhedronIntegrity", R"mydelimiter(
+        The pointing direction of the normals of a Polyhedron.
+        They can either point outwards or inwards the polyhedron.
+        )mydelimiter")
+        .value("DISABLE", PolyhedronIntegrity::DISABLE,
+               "All activities regarding MeshChecking are disabled. No runtime overhead!")
+        .value("VERIFY", PolyhedronIntegrity::VERIFY,
+               "Only verification of the NormalOrientation. "
+               "A misalignment (e.g. specified OUTWARDS, but is not) leads to a runtime_error. Runtime Cost :math:`O(n^2)`")
+        .value("AUTOMATIC", PolyhedronIntegrity::AUTOMATIC,
+               "Like :code:`VERIFY`, but also informs the user about the option in any case on the runtime costs. "
+               "This is the implicit default option. Runtime Cost: :math:`O(n^2)` and output to stdout in every case!")
+        .value("HEAL", PolyhedronIntegrity::HEAL,
+               "Verification and Autmatioc Healing of the NormalOrientation. "
+               "A misalignemt does not lead to a runtime_error, but to an internal correction of vertices ordering. Runtime Cost: :math:`O(n^2)`");
+
+    py::class_<Polyhedron>(m, "Polyhedron", R"mydelimiter(
+            A constant density Polyhedron stores the mesh data consisting of vertices and triangular faces.
+            The density and the coordinate system in which vertices and faces are defined need to have the same scale/ units.
+            Tsoulis et al.'s polyhedral gravity model requires that the plane unit normals of every face are pointing outwards
+            of the polyhedron. Otherwise the results are negated.
+            The class by default enforces this constraints and offers utility to (automatically) make the input data obey to this constraint.
+            )mydelimiter")
+            .def(py::init<const std::variant<PolyhedralSource, PolyhedralFiles> &, double, const NormalOrientation &, const PolyhedronIntegrity &>(), R"mydelimiter(
+            Creates a new Polyhedron from vertices and faces and a constant density.
+            If the integrity_check is not set to DISABLE, the mesh integrity is checked
+            (so that it fits the specification of the polyhedral model by *Tsoulis et al.*)
 
             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)
+                polyhedral_source:  The vertices (:math:`(N, 3)`-array-like) and faces (:math:`(M, 3)`-array-like) of the polyhedron as pair or
+                                    The filenames of the files containing the vertices & faces as list of strings
+                density:            The constant density of the polyhedron, it must match the mesh's units, e.g. mesh in :math:`[m]` then density in :math:`[kg/m^3]`
+                normal_orientation: The pointing direction of the mesh's plane unit normals, i.e., either :code:`OUTWARDS` or :code:`INWARDS` of the polyhedron.
+                                    One of :py:class:`polyhedral_gravity.NormalOrientation`.
+                                    (default: :code:`OUTWARDS`)
+                integrity_check:    Conducts an Integrity Check (degenerated faces/ vertex ordering) depending on the values. One of :py:class:`polyhedral_gravity.PolyhedronIntegrity`:
+
+                                        * :code:`AUTOMATIC` (Default): Prints to stdout and throws ValueError if normal_orientation is wrong/ inconsisten
+                                        * :code:`VERIFY`: Like :code:`AUTOMATIC`, but does not print to stdout
+                                        * :code:`DISABLE`: Recommened, when you know the mesh to avoid to pay :math:`O(n^2)` runtime. Disables ALL checks
+                                        * :code`HEAL`: Automatically fixes the normal_orientation and vertex ordering to the correct values
+
+            Raises:
+                ValueError: If the faces array does not contain a reference to vertex 0 indicating an index start at 1
+                ValueError: If :code:`integrity_check` is set to :code:`AUTOMATIC` or :code:`VERIFY` and the mesh is inconsistent
+
+            Note:
+                The :code:`integrity_check` is automatically enabled to avoid wrong results due to the wrong vertex ordering.
+                The check requires :math:`O(n^2)` operations. You want to turn this off, when you know you mesh!
+            )mydelimiter",
+                 py::arg("polyhedral_source"),
+                 py::arg("density"),
+                 py::arg("normal_orientation") = NormalOrientation::OUTWARDS,
+                 py::arg("integrity_check") = PolyhedronIntegrity::AUTOMATIC
+                    )
+            .def("check_normal_orientation", &Polyhedron::checkPlaneUnitNormalOrientation, R"mydelimiter(
+            Returns a tuple consisting of majority plane unit normal orientation,
+            i.e. the direction in which at least more than half of the plane unit normals point,
+            and the indices of the faces violating this orientation, i.e. the faces whose plane unit normals point in the other direction.
+            The set of incides vioalting the property is empty if the mesh has a clear ordering.
+            The set contains values if the mesh is incosistent.
+
+            Returns:
+                Tuple consisting consisting of majority plane unit normal orientation and the indices of the faces violating this orientation.
+
+            Note:
+                This utility is mainly for diagnostics and debugging purposes. If the polyhedron is constrcuted with `integrity_check`
+                set to :code:`AUTOMATIC` or :code:`VERIFY`, the construction fails anyways.
+                If set to :code:`HEAL`, this method should return an empty set (but maybe a different ordering than initially specified)
+                Only if set to code:`DISABLE`, then this method might actually return a set with faulty indices.
+                Hence, if you want to know your mesh error. Construct the polyhedron with :code:`integrity_check=DISABLE` and call this method.
+            )mydelimiter")
+            .def("__getitem__", &Polyhedron::getResolvedFace, R"mydelimiter(
+            Returns the the three cooridnates of the vertices making the face at the requested index.
+            This does not return the face as list of vertex indices, but resolved with the actual coordinates.
+
+            Args:
+                index:  The index of the face
 
             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.
+                :math:`(3, 3)`-array-like: The resolved face
+
+            Raises:
+                IndexError if face index is out-of-bounds
+            )mydelimiter", py::arg("index"))
+            .def("__repr__", &Polyhedron::toString, R"mydelimiter(
+            :py:class:`str`: A string representation of this polyhedron
+            )mydelimiter")
+            .def_property_readonly("vertices", &Polyhedron::getVertices, R"mydelimiter(
+            (N, 3)-array-like of :py:class:`float`: The vertices of the polyhedron (Read-Only)
+            )mydelimiter")
+            .def_property_readonly("faces", &Polyhedron::getFaces, R"mydelimiter(
+            (M, 3)-array-like of :py:class:`int`: The faces of the polyhedron (Read-Only)
+            )mydelimiter")
+            .def_property("density", &Polyhedron::getDensity, &Polyhedron::setDensity, R"mydelimiter(
+            :py:class:`float`: The density of the polyhedron (Read/ Write)
+            )mydelimiter")
+            .def_property_readonly("normal_orientation", &Polyhedron::getOrientation, R"mydelimiter(
+            :py:class:`polyhedral_gravity.NormalOrientation`: The orientation of the plane unit normals (Read-Only)
+            )mydelimiter")
+            .def(py::pickle(
+                    [](const Polyhedron &polyhedron) {
+                        const auto &[vertices, faces, density, orientation] = polyhedron.getState();
+                        return py::make_tuple(vertices, faces, density, orientation);
+                    },
+                    [](const py::tuple &tuple) {
+                        constexpr size_t tupleSize = 4;
+                        if (tuple.size() != tupleSize) {
+                            throw std::runtime_error("Invalid state!");
+                        }
+                        Polyhedron polyhedron{
+                                tuple[0].cast<std::vector<Array3>>(), tuple[1].cast<std::vector<IndexArray3>>(),
+                                tuple[2].cast<double>(), tuple[3].cast<NormalOrientation>(), PolyhedronIntegrity::DISABLE
+                        };
+                        return polyhedron;
+                    }
+                    ));
+
+    py::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 :py:meth:`poylhedral_gravity.GravityEvaluable.__call__` method to evaluate the polyhedral gravity model for computation points while
+             also caching the polyhedron & intermediate results over the lifetime of the object.
+             )mydelimiter")
+            .def(py::init<const Polyhedron &>(),R"mydelimiter(
+             Creates a new GravityEvaluable for a given constant density polyhedron.
+             It provides a :py:meth:`poylhedral_gravity.GravityEvaluable.__call__` method to evaluate the polyhedral gravity model for computation points while
+             also caching the polyhedron & intermediate results over the lifetime of the object.
+
+             Args:
+                 polyhedron: The polyhedron for which to evaluate the gravity model
+             )mydelimiter", py::arg("polyhedron"))
+            .def("__repr__", &GravityEvaluable::toString,R"mydelimiter(
+            :py:class:`str`: A string representation of this GravityEvaluable
             )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.
+             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)
+             Args:
+                 computation_points: The computation points as tuple or list of points
+                 parallel:           If :code:`True`, the computation is done in parallel (default: :code:`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)
+             Returns:
+                 Either a triplet 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 a list of these triplets
+             )mydelimiter", py::arg("computation_points"), py::arg("parallel") = true)
             .def(py::pickle(
                     [](const GravityEvaluable &evaluable) {
-                        const auto&[polyhedron, density, segmentVectors, planeUnitNormals, segmentUnitNormals] =
-                                evaluable.getState();
-                        return py::make_tuple(polyhedron.getVertices(), polyhedron.getFaces(), density, segmentVectors,
-                                              planeUnitNormals, segmentUnitNormals);
+                        const auto &[polyhedron, segmentVectors, planeUnitNormals, segmentUnitNormals] = evaluable.getState();
+                        return py::make_tuple(polyhedron, segmentVectors, planeUnitNormals, segmentUnitNormals);
                     },
                     [](const py::tuple &tuple) {
-                        constexpr size_t tupleSize = 6;
+                        constexpr size_t tupleSize = 4;
                         if (tuple.size() != tupleSize) {
                             throw std::runtime_error("Invalid state!");
                         }
-                        Polyhedron polyhedron {
-                            tuple[0].cast<std::vector<Array3>>(), tuple[1].cast<std::vector<IndexArray3>>()
-                        };
                         GravityEvaluable evaluable{
-                            polyhedron, tuple[2].cast<double>(), tuple[3].cast<std::vector<Array3Triplet>>(),
-                            tuple[4].cast<std::vector<Array3>>(), tuple[5].cast<std::vector<Array3Triplet>>()
+                                tuple[0].cast<Polyhedron>(), tuple[1].cast<std::vector<Array3Triplet>>(),
+                                tuple[2].cast<std::vector<Array3>>(), tuple[3].cast<std::vector<Array3Triplet>>()
                         };
                         return evaluable;
                     }
-            ));
+                    ));
+
+    m.def("evaluate", [](const Polyhedron &polyhedron,
+                         const std::variant<Array3, std::vector<Array3>> &computationPoints,
+                         bool parallel) -> std::variant<GravityModelResult, std::vector<GravityModelResult>> {
+                    return std::visit(util::overloaded{
+                            [&](const Array3 &point) {
+                                return std::variant<GravityModelResult, std::vector<GravityModelResult>>(GravityModel::evaluate(polyhedron, point, parallel));
+                            },
+                            [&](const std::vector<Array3> &points) {
+                                return std::variant<GravityModelResult, std::vector<GravityModelResult>>(GravityModel::evaluate(polyhedron, points, parallel));
+                            }
+                        }, computationPoints);
+          }, R"mydelimiter(
+             Evaluates the polyhedral gravity model for a given constant density polyhedron at a given computation point.
 
-    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.");
+             Args:
+                 polyhedron:            The polyhedron for which to evaluate the gravity model
+                 computation_points:    The computation points as tuple or list of points
+                 parallel:              If :code:`True`, the computation is done in parallel (default: :code:`True`)
 
-    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!
+             Returns:
+                 Either a triplet 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 a list of these triplets
+             )mydelimiter", py::arg("polyhedron"), py::arg("computation_points"), py::arg("parallel") = true);
 
-            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"));
-}
+}
\ No newline at end of file


Test patch:
diff --git a/.github/workflows/build-and-test.yml b/.github/workflows/build-and-test.yml
index 548883e..90f8701 100644
--- a/.github/workflows/build-and-test.yml
+++ b/.github/workflows/build-and-test.yml
@@ -10,28 +10,41 @@ on:
 
 # A workflow run is made up of one or more jobs that can run sequentially or in parallel
 jobs:
-  # Build & Test on Linux
-  build-and-test-linux:
-    runs-on: ubuntu-latest
+
+  # Build the C++ library & test that it is working using GoogleTest
+  build-and-test-cpp:
+    name: C/C++ Library on ${{ matrix.os }}
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        os: [windows-latest, ubuntu-latest, macos-latest]
     steps:
       - uses: actions/checkout@v3
+      - uses: ilammy/msvc-dev-cmd@v1
+        if: matrix.os == 'windows-latest'
       - name: Build
         run: |
           mkdir build
           cd build
           cmake -DCMAKE_BUILD_TYPE=Release ..
-          cmake --build .
-      - name: Test
+          cmake --build . --config Release
+      - name: Test (Linux & macOS)
         run: cd build/test && ./polyhedralGravity_test
+        if: matrix.os != 'windows-latest'
 
-  # Install the python interface by building from source and run pytest
-  pytest-linux:
-    runs-on: ubuntu-latest
+  # Install the python interface & test that it is working using pytest
+  build-and-test-python:
+    name: Python Interface on ${{ matrix.os }}
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        os: [windows-latest, ubuntu-latest, macos-latest]
     steps:
       - uses: actions/checkout@v3
-      - name: Install Ninja Build
-        run: |
-          sudo apt-get install ninja-build -y
+      - uses: ilammy/msvc-dev-cmd@v1
+        if: matrix.os == 'windows-latest'
       - name: Install Conda environment from environment.yml
         uses: mamba-org/provision-with-micromamba@main
         with:
@@ -44,22 +57,3 @@ jobs:
           pip install . -vv --no-build-isolation
           pytest -n 3
           
-
-  # Builds the polyhedral gravity python interface the interface
-  # No testing on Windows since it takes too long
-  # Enables the long paths feature on Windows (newer thrust versions require this)
-  build-windows:
-    runs-on: windows-latest
-    steps:
-      - uses: actions/checkout@v3
-      - uses: ilammy/msvc-dev-cmd@v1
-      - name: Build
-        run: |
-          New-ItemProperty -Path "HKLM:\SYSTEM\CurrentControlSet\Control\FileSystem" -Name "LongPathsEnabled" -Value 1 -PropertyType DWORD -Force
-          git config --system core.longpaths true
-          mkdir build
-          cd build
-          cmake -DCMAKE_BUILD_TYPE=Release -DBUILD_POLYHEDRAL_PYTHON_INTERFACE=ON ..
-          cmake --build . --config Release
-
-      
diff --git a/test/calculation/MeshCheckingTest.cpp b/test/calculation/MeshCheckingTest.cpp
deleted file mode 100644
index daa13dd..0000000
--- a/test/calculation/MeshCheckingTest.cpp
+++ /dev/null
@@ -1,221 +0,0 @@
-#include "gtest/gtest.h"
-
-#include <array>
-#include <utility>
-#include <tuple>
-#include "polyhedralGravity/calculation/MeshChecking.h"
-#include "polyhedralGravity/model/Polyhedron.h"
-#include "polyhedralGravity/input/TetgenAdapter.h"
-
-
-/**
- * Contains Tests if the sanity checks works
- */
-class MeshCheckingTest : public ::testing::Test {
-
-protected:
-
-    const polyhedralGravity::Polyhedron _correctCube{{
-                                                             {-1.0, -1.0, -1.0},
-                                                             {1.0, -1.0, -1.0},
-                                                             {1.0, 1.0, -1.0},
-                                                             {-1.0, 1.0, -1.0},
-                                                             {-1.0, -1.0, 1.0},
-                                                             {1.0, -1.0, 1.0},
-                                                             {1.0, 1.0, 1.0},
-                                                             {-1.0, 1.0, 1.0}},
-                                                     {
-                                                             {1, 3, 2},
-                                                             {0, 3, 1},
-                                                             {0, 1, 5},
-                                                             {0, 5, 4},
-                                                             {0, 7, 3},
-                                                             {0, 4, 7},
-                                                             {1, 2, 6},
-                                                             {1, 6, 5},
-                                                             {2, 3, 6},
-                                                             {3, 7, 6},
-                                                             {4, 5, 6},
-                                                             {4, 6, 7}}
-    };
-
-    const polyhedralGravity::Polyhedron _wrongCube{{
-                                                           {-1.0, -1.0, -1.0},
-                                                           {1.0, -1.0, -1.0},
-                                                           {1.0, 1.0, -1.0},
-                                                           {-1.0, 1.0, -1.0},
-                                                           {-1.0, -1.0, 1.0},
-                                                           {1.0, -1.0, 1.0},
-                                                           {1.0, 1.0, 1.0},
-                                                           {-1.0, 1.0, 1.0}},
-                                                   {
-                                                           {3, 1, 2},
-                                                           {3, 0, 1},
-                                                           {1, 0, 5},
-                                                           {5, 0, 4},
-                                                           {7, 0, 3},
-                                                           {4, 0, 7},
-                                                           {2, 1, 6},
-                                                           {6, 1, 5},
-                                                           {3, 2, 6},
-                                                           {7, 3, 6},
-                                                           {5, 4, 6},
-                                                           {6, 4, 7}}
-    };
-
-    const polyhedralGravity::Polyhedron _wrongCube2{{
-                                                           {-1.0, -1.0, -1.0},
-                                                           {1.0, -1.0, -1.0},
-                                                           {1.0, 1.0, -1.0},
-                                                           {-1.0, 1.0, -1.0},
-                                                           {-1.0, -1.0, 1.0},
-                                                           {1.0, -1.0, 1.0},
-                                                           {1.0, 1.0, 1.0},
-                                                           {-1.0, 1.0, 1.0}},
-                                                   {
-                                                           {3, 1, 2},
-                                                           {0, 3, 1},
-                                                           {0, 1, 5},
-                                                           {0, 5, 4},
-                                                           {0, 7, 3},
-                                                           {0, 4, 7},
-                                                           {1, 2, 6},
-                                                           {1, 6, 5},
-                                                           {2, 3, 6},
-                                                           {3, 7, 6},
-                                                           {4, 5, 6},
-                                                           {4, 6, 7}}
-    };
-
-    const polyhedralGravity::Polyhedron _degeneratedCube{{
-                                                                 {-1.0, -1.0, -1.0},
-                                                                 {1.0, -1.0, -1.0},
-                                                                 {1.0, 1.0, -1.0},
-                                                                 {-1.0, 1.0, -1.0},
-                                                                 {-1.0, -1.0, 1.0},
-                                                                 {1.0, -1.0, 1.0},
-                                                                 {1.0, 1.0, 1.0},
-                                                                 {-1.0, 1.0, 1.0}},
-                                                         {
-                                                                 {1,    3,    2},
-                                                                 {0,   3,    1},
-                                                                 {0,   1,   5},
-                                                                 {0,    5,   4},
-                                                                 {7,    7,    3},
-                                                                 {0,   4,    7},
-                                                                 {1,   2,   6},
-                                                                 {1,    6,   5},
-                                                                 {2, 3, 6},
-                                                                 {3, 7, 6},
-                                                                 {4, 5, 6},
-                                                                 {4, 6, 7}}
-    };
-
-
-    const polyhedralGravity::Polyhedron _correctPrism{{
-                                                              {-20.0, 0.0, 25.0},
-                                                              {0.0, 0.0, 25.0},
-                                                              {0.0, 10.0, 25.0},
-                                                              {-20.0, 10.0, 25.0},
-                                                              {-20.0, 0.0, 15.0},
-                                                              {0.0, 0.0, 15.0},
-                                                              {0.0, 10.0, 15.0},
-                                                              {-20.0, 10.0, 15.0}},
-                                                      {
-                                                              {0, 4, 5},
-                                                              {0, 5, 1},
-                                                              {0, 1, 3},
-                                                              {1, 2, 3},
-                                                              {1, 5, 6},
-                                                              {1, 6, 2},
-                                                              {0, 7, 4},
-                                                              {0, 3, 7},
-                                                              {4, 7, 5},
-                                                              {5, 7, 6},
-                                                              {2, 7, 3},
-                                                              {2, 6, 7}
-                                                      }};
-
-    const polyhedralGravity::Polyhedron _wrongPrism{{
-                                                            {-20.0, 0.0, 25.0},
-                                                            {0.0, 0.0, 25.0},
-                                                            {0.0, 10.0, 25.0},
-                                                            {-20.0, 10.0, 25.0},
-                                                            {-20.0, 0.0, 15.0},
-                                                            {0.0, 0.0, 15.0},
-                                                            {0.0, 10.0, 15.0},
-                                                            {-20.0, 10.0, 15.0}},
-                                                    {
-                                                            {4, 0, 5},
-                                                            {5, 0, 1},
-                                                            {1, 0, 3},
-                                                            {2, 1, 3},
-                                                            {5, 1, 6},
-                                                            {6, 1, 2},
-                                                            {7, 0, 4},
-                                                            {3, 0, 7},
-                                                            {7, 4, 5},
-                                                            {7, 5, 6},
-                                                            {7, 2, 3},
-                                                            {6, 2, 7}
-                                                    }};
-};
-
-TEST_F(MeshCheckingTest, CorrectCube) {
-    using namespace testing;
-    // All normals are pointing outwards
-    ASSERT_TRUE(polyhedralGravity::MeshChecking::checkNormalsOutwardPointing(_correctCube))
-    << "The vertices of the cube are correctly sorted, however the Sanity Check came to the wrong conclusion!";
-    ASSERT_TRUE(polyhedralGravity::MeshChecking::checkTrianglesNotDegenerated(_correctCube));
-}
-
-TEST_F(MeshCheckingTest, WrongCube) {
-    using namespace testing;
-    // All normals are pointing inwards (reversed order)
-    ASSERT_FALSE(polyhedralGravity::MeshChecking::checkNormalsOutwardPointing(_wrongCube))
-    << "The vertices of the cube are incorrectly sorted, however the Sanity Check failed to notice that!";
-    ASSERT_TRUE(polyhedralGravity::MeshChecking::checkTrianglesNotDegenerated(_wrongCube));
-}
-
-TEST_F(MeshCheckingTest, WrongCube2) {
-    using namespace testing;
-    // All normals are pointing outwards expect one which points inwards (and has a reversed order)
-    ASSERT_FALSE(polyhedralGravity::MeshChecking::checkNormalsOutwardPointing(_wrongCube2))
-    << "The vertices of the cube are incorrectly sorted, however the Sanity Check failed to notice that!";
-    ASSERT_TRUE(polyhedralGravity::MeshChecking::checkTrianglesNotDegenerated(_wrongCube2));
-}
-
-TEST_F(MeshCheckingTest, DegeneratedCube) {
-    using namespace testing;
-    // One surface has the same point twice as vertex
-    ASSERT_FALSE(polyhedralGravity::MeshChecking::checkTrianglesNotDegenerated(_degeneratedCube));
-}
-
-TEST_F(MeshCheckingTest, CorrectPrism) {
-    using namespace testing;
-    // All normals are pointing outwards
-    ASSERT_TRUE(polyhedralGravity::MeshChecking::checkNormalsOutwardPointing(_correctPrism))
-    << "The vertices of the prism are correctly sorted, however the Sanity Check came to the wrong conclusion!";
-    ASSERT_TRUE(polyhedralGravity::MeshChecking::checkTrianglesNotDegenerated(_correctPrism));
-}
-
-TEST_F(MeshCheckingTest, WrongPrism) {
-    using namespace testing;
-    // All normals are pointing inwards (reversed order)
-    ASSERT_FALSE(polyhedralGravity::MeshChecking::checkNormalsOutwardPointing(_wrongPrism))
-    << "The vertices of the prism are incorrectly sorted, however the Sanity Check failed to notice that!";
-    ASSERT_TRUE(polyhedralGravity::MeshChecking::checkTrianglesNotDegenerated(_wrongPrism));
-}
-
-TEST_F(MeshCheckingTest, CorrectBigPolyhedron) {
-    using namespace testing;
-    using namespace polyhedralGravity;
-    // All normals are pointing outwards, extensive Eros example
-    Polyhedron polyhedron{
-            TetgenAdapter{
-                    {"resources/GravityModelBigTest.node", "resources/GravityModelBigTest.face"}}.getPolyhedron()};
-    ASSERT_TRUE(polyhedralGravity::MeshChecking::checkNormalsOutwardPointing(polyhedron))
-    << "The vertices of the polyhedron are correctly sorted, however the Sanity Check came to the wrong conclusion!";
-    ASSERT_TRUE(polyhedralGravity::MeshChecking::checkTrianglesNotDegenerated(polyhedron));
-}
-
diff --git a/test/input/TetgenAdapterTest.cpp b/test/input/TetgenAdapterTest.cpp
index cba98e0..0459b6a 100644
--- a/test/input/TetgenAdapterTest.cpp
+++ b/test/input/TetgenAdapterTest.cpp
@@ -9,7 +9,7 @@ class TetgenAdapterTest : public ::testing::Test {
 
 protected:
 
-    std::vector<std::array<double, 3>> _expectedNodes = {
+    std::vector<std::array<double, 3>> _expectedVertices = {
             {-20, 0,  25},
             {0,   0,  25},
             {0,   10, 25},
@@ -47,9 +47,9 @@ TEST_F(TetgenAdapterTest, readSimpleNode) {
     };
 
     TetgenAdapter tetgenAdapter{simpleFiles};
-    auto actualPolyhedron = tetgenAdapter.getPolyhedron();
+    const auto&[actualVertices, actualFaces] = tetgenAdapter.getPolyhedralSource();
 
-    ASSERT_THAT(actualPolyhedron.getVertices(), ContainerEq(_expectedNodes));
+    ASSERT_THAT(actualVertices, ContainerEq(_expectedVertices));
 }
 
 
@@ -63,9 +63,9 @@ TEST_F(TetgenAdapterTest, readSimpleFace) {
     };
 
     TetgenAdapter tetgenAdapter{simpleFiles};
-    auto actualPolyhedron = tetgenAdapter.getPolyhedron();
+    const auto&[actualVertices, actualFaces] = tetgenAdapter.getPolyhedralSource();
 
-    ASSERT_THAT(actualPolyhedron.getFaces(), ContainerEq(_expectedFaces));
+    ASSERT_THAT(actualFaces, ContainerEq(_expectedFaces));
 }
 
 TEST_F(TetgenAdapterTest, readSimpleMesh) {
@@ -75,12 +75,12 @@ TEST_F(TetgenAdapterTest, readSimpleMesh) {
     std::vector<std::string> simpleFiles{"resources/TetgenAdapterTestReadSimple.mesh"};
 
     TetgenAdapter tetgenAdapter{simpleFiles};
-    auto actualPolyhedron = tetgenAdapter.getPolyhedron();
+    const auto&[actualVertices, actualFaces] = tetgenAdapter.getPolyhedralSource();
 
-    for (const auto &actualVertice: actualPolyhedron.getVertices()) {
-        ASSERT_THAT(_expectedNodes, Contains(actualVertice));
+    for (const auto &actualVertice: actualVertices) {
+        ASSERT_THAT(_expectedVertices, Contains(actualVertice));
     }
-    ASSERT_EQ(_expectedFaces.size(), actualPolyhedron.countFaces());
+    ASSERT_EQ(_expectedFaces.size(), actualFaces.size());
 }
 
 TEST_F(TetgenAdapterTest, readSimpleOff) {
@@ -90,12 +90,12 @@ TEST_F(TetgenAdapterTest, readSimpleOff) {
     std::vector<std::string> simpleFiles{"resources/TetgenAdapterTestReadSimple.off"};
 
     TetgenAdapter tetgenAdapter{simpleFiles};
-    auto actualPolyhedron = tetgenAdapter.getPolyhedron();
+    const auto&[actualVertices, actualFaces] = tetgenAdapter.getPolyhedralSource();
 
-    for (const auto &actualVertice: actualPolyhedron.getVertices()) {
-        ASSERT_THAT(_expectedNodes, Contains(actualVertice));
+    for (const auto &actualVertice: actualVertices) {
+        ASSERT_THAT(_expectedVertices, Contains(actualVertice));
     }
-    ASSERT_EQ(_expectedFaces.size(), actualPolyhedron.countFaces());
+    ASSERT_EQ(_expectedFaces.size(), actualFaces.size());
 }
 
 TEST_F(TetgenAdapterTest, readSimplePly) {
@@ -105,12 +105,12 @@ TEST_F(TetgenAdapterTest, readSimplePly) {
     std::vector<std::string> simpleFiles{"resources/TetgenAdapterTestReadSimple.ply"};
 
     TetgenAdapter tetgenAdapter{simpleFiles};
-    auto actualPolyhedron = tetgenAdapter.getPolyhedron();
+    const auto&[actualVertices, actualFaces] = tetgenAdapter.getPolyhedralSource();
 
-    for (const auto &actualVertice: actualPolyhedron.getVertices()) {
-        ASSERT_THAT(_expectedNodes, Contains(actualVertice));
+    for (const auto &actualVertice: actualVertices) {
+        ASSERT_THAT(_expectedVertices, Contains(actualVertice));
     }
-    ASSERT_EQ(_expectedFaces.size(), actualPolyhedron.countFaces());
+    ASSERT_EQ(_expectedFaces.size(), actualFaces.size());
 }
 
 TEST_F(TetgenAdapterTest, readSimpleStl) {
@@ -120,10 +120,10 @@ TEST_F(TetgenAdapterTest, readSimpleStl) {
     std::vector<std::string> simpleFiles{"resources/TetgenAdapterTestReadSimple.stl"};
 
     TetgenAdapter tetgenAdapter{simpleFiles};
-    auto actualPolyhedron = tetgenAdapter.getPolyhedron();
+    const auto&[actualVertices, actualFaces] = tetgenAdapter.getPolyhedralSource();
 
-    for (const auto &actualVertice: actualPolyhedron.getVertices()) {
-        ASSERT_THAT(_expectedNodes, Contains(actualVertice));
+    for (const auto &actualVertice: actualVertices) {
+        ASSERT_THAT(_expectedVertices, Contains(actualVertice));
     }
-    ASSERT_EQ(_expectedFaces.size(), actualPolyhedron.countFaces());
+    ASSERT_EQ(_expectedFaces.size(), actualFaces.size());
 }
\ No newline at end of file
diff --git a/test/calculation/GravityModelBigTest.cpp b/test/model/GravityModelBigTest.cpp
similarity index 96%
rename from test/calculation/GravityModelBigTest.cpp
rename to test/model/GravityModelBigTest.cpp
index e8ab66f..e80285d 100644
--- a/test/calculation/GravityModelBigTest.cpp
+++ b/test/model/GravityModelBigTest.cpp
@@ -6,9 +6,7 @@
 #include <utility>
 #include <fstream>
 #include <sstream>
-#include <string>
-#include "polyhedralGravity/input/TetgenAdapter.h"
-#include "polyhedralGravity/calculation/GravityModel.h"
+#include "polyhedralGravity/model/GravityModel.h"
 #include "polyhedralGravity/model/Polyhedron.h"
 
 #include "GravityModelVectorUtility.h"
@@ -21,6 +19,10 @@
  * The values which are used to check the C++ implementation are calculated by the
  * Tsoulis reference implementation in FORTRAN and saved in the files in test/resources.
  *
+ * This test case might fail locally using a different architecture/ compiler than the CI!
+ * This is okay - A fix would be to check container element verbosly in a loop with an epsilon.
+ * However, the ContainerEq() - as currentl used - matcher makes the test a lot more readbale.
+ *
  */
 class GravityModelBigTest : public ::testing::Test {
 
@@ -35,8 +37,11 @@ class GravityModelBigTest : public ::testing::Test {
     static constexpr size_t LOCAL_TEST_COUNT_NODES_PER_FACE = 3;
 
     polyhedralGravity::Polyhedron _polyhedron{
-            polyhedralGravity::TetgenAdapter{
-                    {"resources/GravityModelBigTest.node", "resources/GravityModelBigTest.face"}}.getPolyhedron()};
+        std::vector<std::string>{"resources/GravityModelBigTest.node", "resources/GravityModelBigTest.face"},
+        1.0,
+        polyhedralGravity::NormalOrientation::OUTWARDS,
+        polyhedralGravity::PolyhedronIntegrity::DISABLE
+    };
 
     std::array<double, 3> _computationPoint{0.0, 0.0, 0.0};
 
diff --git a/test/calculation/GravityModelCubeTest.cpp b/test/model/GravityModelCubeTest.cpp
similarity index 55%
rename from test/calculation/GravityModelCubeTest.cpp
rename to test/model/GravityModelCubeTest.cpp
index f977e91..536357f 100644
--- a/test/calculation/GravityModelCubeTest.cpp
+++ b/test/model/GravityModelCubeTest.cpp
@@ -8,8 +8,7 @@
 #include <tuple>
 #include <fstream>
 #include <sstream>
-#include <string>
-#include "polyhedralGravity/calculation/GravityModel.h"
+#include "polyhedralGravity/model/GravityModel.h"
 #include "polyhedralGravity/model/Polyhedron.h"
 
 
@@ -20,37 +19,38 @@ class GravityModelCubeTest :
         public ::testing::TestWithParam<std::tuple<std::array<double, 3>, double, double, std::array<double, 3>>> {
 
 protected:
-
     /**
      * Small epsilon since we compare to ground truth
      */
     static constexpr double LOCAL_TEST_EPSILON = 1e-20;
-    const polyhedralGravity::Polyhedron _cube{{
-                                                      {-1.0, -1.0, -1.0},
-                                                      {1.0, -1.0, -1.0},
-                                                      {1.0, 1.0, -1.0},
-                                                      {-1.0, 1.0, -1.0},
-                                                      {-1.0, -1.0, 1.0},
-                                                      {1.0, -1.0, 1.0},
-                                                      {1.0, 1.0, 1.0},
-                                                      {-1.0, 1.0, 1.0}},
-                                              {
-                                                      {1,    3,    2},
-                                                      {0,   3,    1},
-                                                      {0,   1,   5},
-                                                      {0,    5,   4},
-                                                      {0,    7,    3},
-                                                      {0,   4,    7},
-                                                      {1,   2,   6},
-                                                      {1,    6,   5},
-                                                      {2, 3, 6},
-                                                      {3, 7, 6},
-                                                      {4, 5, 6},
-                                                      {4, 6, 7}}
+    polyhedralGravity::Polyhedron _cube{std::vector<polyhedralGravity::Array3>{
+                                                {-1.0, -1.0, -1.0},
+                                                {1.0, -1.0, -1.0},
+                                                {1.0, 1.0, -1.0},
+                                                {-1.0, 1.0, -1.0},
+                                                {-1.0, -1.0, 1.0},
+                                                {1.0, -1.0, 1.0},
+                                                {1.0, 1.0, 1.0},
+                                                {-1.0, 1.0, 1.0}},
+                                        std::vector<polyhedralGravity::IndexArray3>{
+                                                {1, 3, 2},
+                                                {0, 3, 1},
+                                                {0, 1, 5},
+                                                {0, 5, 4},
+                                                {0, 7, 3},
+                                                {0, 4, 7},
+                                                {1, 2, 6},
+                                                {1, 6, 5},
+                                                {2, 3, 6},
+                                                {3, 7, 6},
+                                                {4, 5, 6},
+                                                {4, 6, 7}},
+                                        1.0,
+                                        polyhedralGravity::NormalOrientation::OUTWARDS,
+                                        polyhedralGravity::PolyhedronIntegrity::DISABLE
     };
 
 public:
-
     [[nodiscard]] static std::vector<std::tuple<std::array<double, 3>, double, double, std::array<double, 3>>>
     readCubePoints(const std::string &filename) {
         std::vector<std::tuple<std::array<double, 3>, double, double, std::array<double, 3>>> result{};
@@ -76,8 +76,6 @@ class GravityModelCubeTest :
     }
 
 
-
-
 };
 
 TEST_P(GravityModelCubeTest, CubePoints) {
@@ -90,7 +88,8 @@ TEST_P(GravityModelCubeTest, CubePoints) {
     double expectedPotential = std::get<2>(testData);
     std::array<double, 3> expectedAcceleration = std::get<3>(testData);
 
-    const auto[actualPotential, actualAcceleration, x] = GravityModel::evaluate(_cube, density, computationPoint);
+    _cube.setDensity(density);
+    const auto [actualPotential, actualAcceleration, x] = GravityModel::evaluate(_cube, computationPoint);
 
     ASSERT_NEAR(actualPotential, expectedPotential, LOCAL_TEST_EPSILON);
     ASSERT_THAT(actualAcceleration, Pointwise(DoubleNear(LOCAL_TEST_EPSILON), expectedAcceleration));
@@ -98,8 +97,8 @@ TEST_P(GravityModelCubeTest, CubePoints) {
 
 INSTANTIATE_TEST_SUITE_P(CubeGravityModelTest01, GravityModelCubeTest,
                          ::testing::ValuesIn(
-                                 GravityModelCubeTest::readCubePoints("resources/analytic_cube_solution_density1.txt")));
+                             GravityModelCubeTest::readCubePoints("resources/analytic_cube_solution_density1.txt")));
 
 INSTANTIATE_TEST_SUITE_P(CubeGravityModelTest42, GravityModelCubeTest,
                          ::testing::ValuesIn(
-                                 GravityModelCubeTest::readCubePoints("resources/analytic_cube_solution_density42.txt")));
+                             GravityModelCubeTest::readCubePoints("resources/analytic_cube_solution_density42.txt")));
diff --git a/test/calculation/GravityModelTest.cpp b/test/model/GravityModelTest.cpp
similarity index 81%
rename from test/calculation/GravityModelTest.cpp
rename to test/model/GravityModelTest.cpp
index 1b728df..418369a 100644
--- a/test/calculation/GravityModelTest.cpp
+++ b/test/model/GravityModelTest.cpp
@@ -4,7 +4,7 @@
 #include <string>
 #include <vector>
 #include <utility>
-#include "polyhedralGravity/calculation/GravityModel.h"
+#include "polyhedralGravity/model/GravityModel.h"
 #include "polyhedralGravity/model/Polyhedron.h"
 
 #include "GravityModelVectorUtility.h"
@@ -16,7 +16,6 @@
 class GravityModelTest : public ::testing::Test {
 
 protected:
-
     /**
      * Relative big epsilon due to deviations between FORTRAN implementation and C++ implementation
      */
@@ -24,31 +23,34 @@ class GravityModelTest : public ::testing::Test {
 
     //New polyhedron with given vertices and faces
     //this is the base example from Tsoulis
-    polyhedralGravity::Polyhedron _polyhedron{
-            {
-                    {-20, 0, 25},
-                    {0, 0, 25},
-                    {0, 10, 25},
-                    {-20, 10, 25},
-                    {-20, 0, 15},
-                    {0, 0, 15},
-                    {0, 10, 15},
-                    {-20, 10, 15}
-            },
-            {
-                    {0,   1, 3},
-                    {1, 2, 3},
-                    {0, 4,  5},
-                    {0,   5,  1},
-                    {0,   7, 4},
-                    {0, 3, 7},
-                    {1, 5,  6},
-                    {1,   6,  2},
-                    {3, 6, 7},
-                    {2, 6, 3},
-                    {4, 6, 5},
-                    {4, 7, 6}
-            }
+    polyhedralGravity::Polyhedron _polyhedron {
+        std::vector<polyhedralGravity::Array3>{
+                        {-20, 0, 25},
+                        {0, 0, 25},
+                        {0, 10, 25},
+                        {-20, 10, 25},
+                        {-20, 0, 15},
+                        {0, 0, 15},
+                        {0, 10, 15},
+                        {-20, 10, 15}
+                },
+                std::vector<polyhedralGravity::IndexArray3>{
+                        {0, 1, 3},
+                        {1, 2, 3},
+                        {0, 4, 5},
+                        {0, 5, 1},
+                        {0, 7, 4},
+                        {0, 3, 7},
+                        {1, 5, 6},
+                        {1, 6, 2},
+                        {3, 6, 7},
+                        {2, 6, 3},
+                        {4, 6, 5},
+                        {4, 7, 6}
+                },
+                1.0,
+                polyhedralGravity::NormalOrientation::OUTWARDS,
+                polyhedralGravity::PolyhedronIntegrity::DISABLE
     };
 
     std::array<double, 3> _computationPoint{0.0, 0.0, 0.0};
@@ -69,18 +71,18 @@ class GravityModelTest : public ::testing::Test {
     };
 
     std::vector<std::array<double, 3>> expectedPlaneUnitNormals{
-            {0.0,  -0.0, 1.0},
-            {0.0,  -0.0, 1.0},
-            {0.0,  -1.0, 0.0},
-            {0.0,  -1.0, 0.0},
+            {0.0, -0.0, 1.0},
+            {0.0, -0.0, 1.0},
+            {0.0, -1.0, 0.0},
+            {0.0, -1.0, 0.0},
             {-1.0, -0.0, -0.0},
-            {-1.0, 0.0,  0.0},
-            {1.0,  -0.0, 0.0},
-            {1.0,  -0.0, 0.0},
-            {0.0,  1.0,  0.0},
-            {0.0,  1.0,  0.0},
-            {0.0,  0.0,  -1.0},
-            {0.0,  0.0,  -1.0}
+            {-1.0, 0.0, 0.0},
+            {1.0, -0.0, 0.0},
+            {1.0, -0.0, 0.0},
+            {0.0, 1.0, 0.0},
+            {0.0, 1.0, 0.0},
+            {0.0, 0.0, -1.0},
+            {0.0, 0.0, -1.0}
     };
 
     std::vector<std::array<std::array<double, 3>, 3>> expectedSegmentUnitNormals{
@@ -113,18 +115,18 @@ class GravityModelTest : public ::testing::Test {
     std::vector<double> expectedPlaneNormalOrientations{1.0, 1.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 1.0, 1.0, -1.0, -1.0};
 
     std::vector<polyhedralGravity::HessianPlane> expectedHessianPlanes{
-            {0.0,    0.0,    200.0,  -5000.0},
-            {0.0,    -0.0,   200.0,  -5000.0},
-            {0.0,    -200.0, 0.0,    0.0},
-            {0.0,    -200.0, 0.0,    0.0},
-            {-100.0, 0.0,    0.0,    -2000.0},
-            {-100.0, 0.0,    0.0,    -2000.0},
-            {100.0,  0.0,    0.0,    -0.0},
-            {100.0,  -0.0,   0.0,    0.0},
-            {0.0,    200.0,  0.0,    -2000.0},
-            {0.0,    200.0,  0.0,    -2000.0},
-            {0.0,    0.0,    -200.0, 3000.0},
-            {0.0,    0.0,    -200.0, 3000.0}
+            {0.0, 0.0, 200.0, -5000.0},
+            {0.0, -0.0, 200.0, -5000.0},
+            {0.0, -200.0, 0.0, 0.0},
+            {0.0, -200.0, 0.0, 0.0},
+            {-100.0, 0.0, 0.0, -2000.0},
+            {-100.0, 0.0, 0.0, -2000.0},
+            {100.0, 0.0, 0.0, -0.0},
+            {100.0, -0.0, 0.0, 0.0},
+            {0.0, 200.0, 0.0, -2000.0},
+            {0.0, 200.0, 0.0, -2000.0},
+            {0.0, 0.0, -200.0, 3000.0},
+            {0.0, 0.0, -200.0, 3000.0}
     };
 
     std::vector<double> expectedPlaneDistances{25.0,
@@ -141,33 +143,33 @@ class GravityModelTest : public ::testing::Test {
                                                15.0};
 
     std::vector<std::array<double, 3>> expectedOrthogonalProjectionPointsOnPlane{
-            {0.0,   0.0,  25.0},
-            {0.0,   0.0,  25.0},
-            {0.0,   0.0,  0.0},
-            {0.0,   0.0,  0.0},
-            {-20.0, 0.0,  0.0},
-            {-20.0, 0.0,  0.0},
-            {0.0,   0.0,  0.0},
-            {0.0,   0.0,  0.0},
-            {0.0,   10.0, 0.0},
-            {0.0,   10.0, 0.0},
-            {0.0,   0.0,  15.0},
-            {0.0,   0.0,  15.0}
+            {0.0, 0.0, 25.0},
+            {0.0, 0.0, 25.0},
+            {0.0, 0.0, 0.0},
+            {0.0, 0.0, 0.0},
+            {-20.0, 0.0, 0.0},
+            {-20.0, 0.0, 0.0},
+            {0.0, 0.0, 0.0},
+            {0.0, 0.0, 0.0},
+            {0.0, 10.0, 0.0},
+            {0.0, 10.0, 0.0},
+            {0.0, 0.0, 15.0},
+            {0.0, 0.0, 15.0}
     };
 
     std::vector<std::array<double, 3>> expectedSegmentNormalOrientations{
-            {0.0,  0.0,  1.0},
-            {0.0,  1.0,  0.0},
-            {1.0,  -1.0, 1.0},
-            {-1.0, 0.0,  1.0},
-            {1.0,  -1.0, 0.0},
-            {1.0,  1.0,  -1.0},
-            {0.0,  -1.0, 1.0},
-            {-1.0, 1.0,  1.0},
-            {1.0,  -1.0, 1.0},
-            {0.0,  -1.0, 1.0},
-            {1.0,  0.0,  0.0},
-            {1.0,  1.0,  -1.0}
+            {0.0, 0.0, 1.0},
+            {0.0, 1.0, 0.0},
+            {1.0, -1.0, 1.0},
+            {-1.0, 0.0, 1.0},
+            {1.0, -1.0, 0.0},
+            {1.0, 1.0, -1.0},
+            {0.0, -1.0, 1.0},
+            {-1.0, 1.0, 1.0},
+            {1.0, -1.0, 1.0},
+            {0.0, -1.0, 1.0},
+            {1.0, 0.0, 0.0},
+            {1.0, 1.0, -1.0}
     };
 
     std::vector<std::array<std::array<double, 3>, 3>> expectedOrthogonalProjectionPointsOnSegment{
@@ -186,18 +188,18 @@ class GravityModelTest : public ::testing::Test {
     };
 
     std::vector<std::array<double, 3>> expectedSegmentDistances{
-            {0.0,                0.0,                20.0},
-            {0.0,                10.0,               0.0},
-            {20.0,               15.0,               13.416407864998739},
-            {13.416407864998739, 0.0,                25.0},
-            {17.67766952966369,  15.0,               0.0},
-            {25.0,               10.0,               17.67766952966369},
-            {0.0,                15.0,               17.67766952966369},
-            {17.67766952966369,  10.0,               25.0},
-            {13.416407864998739, 15.0,               20.0},
-            {0.0,                13.416407864998739, 25.0},
-            {8.94427190999916,   0.0,                0.0},
-            {20.0,               10.0,               8.94427190999916}
+            {0.0, 0.0, 20.0},
+            {0.0, 10.0, 0.0},
+            {20.0, 15.0, 13.416407864998739},
+            {13.416407864998739, 0.0, 25.0},
+            {17.67766952966369, 15.0, 0.0},
+            {25.0, 10.0, 17.67766952966369},
+            {0.0, 15.0, 17.67766952966369},
+            {17.67766952966369, 10.0, 25.0},
+            {13.416407864998739, 15.0, 20.0},
+            {0.0, 13.416407864998739, 25.0},
+            {8.94427190999916, 0.0, 0.0},
+            {20.0, 10.0, 8.94427190999916}
     };
 
     std::vector<std::array<std::array<double, 2>, 3>> expected3DDistancesPerSegmentEndpoint{
@@ -254,33 +256,33 @@ class GravityModelTest : public ::testing::Test {
     std::vector<std::array<polyhedralGravity::Distance, 3>> expectedDistancesPerSegmentEndpoint;
 
     std::vector<std::array<double, 3>> expectedTranscendentalLN{
-            {0.0,                0.0,                 0.30747952872839945},
-            {0.0,                0.687362255356451,   0.0},
-            {0.3544458320893136, 1.0986122886681098,  1.0345679811316213},
-            {1.034567981131622,  0.5108256237659907,  0.7326682560454109},
-            {0.4894110007366263, 0.3900353197707153,  0.3544458320893134},
+            {0.0, 0.0, 0.30747952872839945},
+            {0.0, 0.687362255356451, 0.0},
+            {0.3544458320893136, 1.0986122886681098, 1.0345679811316213},
+            {1.034567981131622, 0.5108256237659907, 0.7326682560454109},
+            {0.4894110007366263, 0.3900353197707153, 0.3544458320893134},
             {0.3074795287283993, 0.33382573681901684, 0.4894110007366262},
-            {-0.510825623765990, 0.6251451172504167,  0.6826834766703017},
-            {0.6826834766703017, 0.4524679290839864,  0.3900353197707153},
-            {0.9286653985398196, 0.9566555518497877,  0.33382573681901667},
-            {0.4524679290839866, 0.928665398539819,   0.6873622553564511},
-            {1.1518034938098078, 0.0,                 0.0},
-            {0.3900353197707153, 0.9566555518497877,  1.1518034938098078}
+            {-0.510825623765990, 0.6251451172504167, 0.6826834766703017},
+            {0.6826834766703017, 0.4524679290839864, 0.3900353197707153},
+            {0.9286653985398196, 0.9566555518497877, 0.33382573681901667},
+            {0.4524679290839866, 0.928665398539819, 0.6873622553564511},
+            {1.1518034938098078, 0.0, 0.0},
+            {0.3900353197707153, 0.9566555518497877, 1.1518034938098078}
     };
 
     std::vector<std::array<double, 3>> expectedTranscendentalAN{
-            {0.0,                 0.0,                 0.3567333885140938},
-            {0.0,                 0.9799235766494776,  0.0},
-            {0.0,                 0.0,                 0.0},
-            {0.0,                 0.0,                 0.0},
-            {0.4109023045514107,  0.45979025757734426, 0.0},
-            {0.23413936163132537, 0.1405746311094993,  0.4109023045514107},
-            {0.0,                 0.0,                 0.0},
-            {0.0,                 0.0,                 0.0},
-            {0.3029908626228055,  0.45979025757734426, 0.08507626483651975},
-            {0.0,                 0.3029908626228055,  0.23413936163132537},
-            {1.2703024256629791,  0.0,                 0.0},
-            {0.27165712367757405, 0.8393489455399783,  1.2703024256629791}
+            {0.0, 0.0, 0.3567333885140938},
+            {0.0, 0.9799235766494776, 0.0},
+            {0.0, 0.0, 0.0},
+            {0.0, 0.0, 0.0},
+            {0.4109023045514107, 0.45979025757734426, 0.0},
+            {0.23413936163132537, 0.1405746311094993, 0.4109023045514107},
+            {0.0, 0.0, 0.0},
+            {0.0, 0.0, 0.0},
+            {0.3029908626228055, 0.45979025757734426, 0.08507626483651975},
+            {0.0, 0.3029908626228055, 0.23413936163132537},
+            {1.2703024256629791, 0.0, 0.0},
+            {0.27165712367757405, 0.8393489455399783, 1.2703024256629791}
     };
 
     std::vector<std::array<polyhedralGravity::TranscendentalExpression, 3>> expectedTranscendentalExpressions;
@@ -291,10 +293,10 @@ class GravityModelTest : public ::testing::Test {
                            std::array<double, 3>{-0.0, -0.0, -0.46364760900080615}),
             std::make_pair(-27.67871794485226,
                            std::array<double, 3>{-0.0, -0.0, -1.1071487177940904}),
-            std::make_pair(0.0, std::array<double, 3>{-0.0, 0.0,  -0.0}),
-            std::make_pair(0.0, std::array<double, 3>{-0.0, 0.0,  -0.0}),
-            std::make_pair(0.0, std::array<double, 3>{0.0,  0.0,  0.0}),
-            std::make_pair(0.0, std::array<double, 3>{0.0,  -0.0, -0.0}),
+            std::make_pair(0.0, std::array<double, 3>{-0.0, 0.0, -0.0}),
+            std::make_pair(0.0, std::array<double, 3>{-0.0, 0.0, -0.0}),
+            std::make_pair(0.0, std::array<double, 3>{0.0, 0.0, 0.0}),
+            std::make_pair(0.0, std::array<double, 3>{0.0, -0.0, -0.0}),
             std::make_pair(0.0, std::array<double, 3>{-0.0, -0.0, -0.0}),
             std::make_pair(0.0, std::array<double, 3>{-0.0, -0.0, -0.0}),
             std::make_pair(0.0, std::array<double, 3>{-0.0, -0.0, -0.0}),
@@ -310,10 +312,10 @@ class GravityModelTest : public ::testing::Test {
     std::vector<std::array<double, 3>> expectedBetaSingularityTerms{
             {-0.0, -0.0, -0.46364760900080615},
             {-0.0, -0.0, -1.1071487177940904},
-            {-0.0, 0.0,  -0.0},
-            {-0.0, 0.0,  -0.0},
-            {0.0,  0.0,  0.0},
-            {0.0,  -0.0, -0.0},
+            {-0.0, 0.0, -0.0},
+            {-0.0, 0.0, -0.0},
+            {0.0, 0.0, 0.0},
+            {0.0, -0.0, -0.0},
             {-0.0, -0.0, -0.0},
             {-0.0, -0.0, -0.0},
             {-0.0, -0.0, -0.0},
@@ -323,7 +325,6 @@ class GravityModelTest : public ::testing::Test {
     };
 
 public:
-
     GravityModelTest() : ::testing::Test() {
         expectedDistancesPerSegmentEndpoint.resize(expected3DDistancesPerSegmentEndpoint.size());
         expectedTranscendentalExpressions.resize(expectedTranscendentalLN.size());
diff --git a/test/calculation/GravityModelVectorUtility.cpp b/test/model/GravityModelVectorUtility.cpp
similarity index 96%
rename from test/calculation/GravityModelVectorUtility.cpp
rename to test/model/GravityModelVectorUtility.cpp
index b48a9bf..cd9071e 100644
--- a/test/calculation/GravityModelVectorUtility.cpp
+++ b/test/model/GravityModelVectorUtility.cpp
@@ -44,7 +44,7 @@ namespace polyhedralGravity {
     GravityModel::calculatePlaneNormalOrientations(const Array3 &computationPoint, const Polyhedron &polyhedron,
                                                    const std::vector<Array3> &planeUnitNormals) {
         std::vector<double> planeNormalOrientations(planeUnitNormals.size(), 0.0);
-        auto transformedPolyhedronIt = transformPolyhedron(polyhedron, computationPoint);
+        auto transformedPolyhedronIt = polyhedron.transformIterator(computationPoint);
         //Calculate sigma_p for every plane given as input: N_p and vertex0 of every face
         std::transform(planeUnitNormals.cbegin(), planeUnitNormals.cend(), transformedPolyhedronIt.first,
                        planeNormalOrientations.begin(),
@@ -58,7 +58,7 @@ namespace polyhedralGravity {
     std::vector<HessianPlane>
     GravityModel::calculateFacesToHessianPlanes(const Array3 &computationPoint, const Polyhedron &polyhedron) {
         std::vector<HessianPlane> hessianPlanes{polyhedron.countFaces()};
-        auto transformedPolyhedronIt = transformPolyhedron(polyhedron, computationPoint);
+        auto transformedPolyhedronIt = polyhedron.transformIterator(computationPoint);
         //Calculate for each face/ plane/ triangle (here) the Hessian Plane
         std::transform(transformedPolyhedronIt.first, transformedPolyhedronIt.second, hessianPlanes.begin(),
                        [](const Array3Triplet &face) -> HessianPlane {
@@ -106,7 +106,7 @@ namespace polyhedralGravity {
             const std::vector<Array3> &orthogonalProjectionPointsOnPlane) {
         std::vector<Array3> segmentNormalOrientations{segmentUnitNormals.size()};
 
-        auto transformedPolyhedronIt = transformPolyhedron(polyhedron, computationPoint);
+        auto transformedPolyhedronIt = polyhedron.transformIterator(computationPoint);
         auto first = util::zip(transformedPolyhedronIt.first,
                                orthogonalProjectionPointsOnPlane.begin(), segmentUnitNormals.begin());
         auto last = util::zip(transformedPolyhedronIt.second,
@@ -132,7 +132,7 @@ namespace polyhedralGravity {
         std::vector<Array3Triplet> orthogonalProjectionPointsOnSegments{orthogonalProjectionPointsOnPlane.size()};
 
         //Zip the three required arguments together: P' for every plane, sigma_pq for every segment, the faces
-        auto transformedPolyhedronIt = transformPolyhedron(polyhedron, computationPoint);
+        auto transformedPolyhedronIt = polyhedron.transformIterator(computationPoint);
         auto first = util::zip(orthogonalProjectionPointsOnPlane.begin(), segmentNormalOrientation.begin(),
                                transformedPolyhedronIt.first);
         auto last = util::zip(orthogonalProjectionPointsOnPlane.end(), segmentNormalOrientation.end(),
@@ -172,7 +172,7 @@ namespace polyhedralGravity {
         std::vector<std::array<Distance, 3>> distances{segmentVectors.size()};
 
         //Zip the three required arguments together: G_ij for every segment, P'' for every segment
-        auto transformedPolyhedronIt = transformPolyhedron(polyhedron, computationPoint);
+        auto transformedPolyhedronIt = polyhedron.transformIterator(computationPoint);
         auto first = util::zip(segmentVectors.begin(), orthogonalProjectionPointsOnSegment.begin(),
                                transformedPolyhedronIt.first);
         auto last = util::zip(segmentVectors.end(), orthogonalProjectionPointsOnSegment.end(),
@@ -199,7 +199,7 @@ namespace polyhedralGravity {
         std::vector<std::array<TranscendentalExpression, 3>> transcendentalExpressions{distances.size()};
 
         //Zip iterator consisting of  3D and 1D distances l1/l2 and s1/2 | h_p | h_pq | sigma_pq | P'_p | faces
-        auto transformedPolyhedronIt = transformPolyhedron(polyhedron, computationPoint);
+        auto transformedPolyhedronIt = polyhedron.transformIterator(computationPoint);
         auto first = util::zip(distances.begin(), planeDistances.begin(), segmentDistances.begin(),
                                segmentNormalOrientation.begin(), orthogonalProjectionPointsOnPlane.begin(),
                                transformedPolyhedronIt.first);
@@ -238,7 +238,7 @@ namespace polyhedralGravity {
         std::vector<std::pair<double, Array3>> singularities{planeDistances.size()};
 
         //Zip iterator consisting of G_ij vectors | sigma_pq | faces | P' | h_p | sigma_p | N_i
-        auto transformedPolyhedronIt = transformPolyhedron(polyhedron, computationPoint);
+        auto transformedPolyhedronIt = polyhedron.transformIterator(computationPoint);
         auto first = util::zip(segmentVectors.begin(), segmentNormalOrientation.begin(),
                                orthogonalProjectionPointsOnPlane.begin(), planeUnitNormals.begin(),
                                planeDistances.begin(), planeNormalOrientations.begin(), transformedPolyhedronIt.first);
diff --git a/test/calculation/GravityModelVectorUtility.h b/test/model/GravityModelVectorUtility.h
similarity index 99%
rename from test/calculation/GravityModelVectorUtility.h
rename to test/model/GravityModelVectorUtility.h
index a2d651b..cd6a458 100644
--- a/test/calculation/GravityModelVectorUtility.h
+++ b/test/model/GravityModelVectorUtility.h
@@ -16,7 +16,7 @@
 #include "thrust/transform_reduce.h"
 #include "thrust/execution_policy.h"
 #include "polyhedralGravity/util/UtilityThrust.h"
-#include "polyhedralGravity/calculation/GravityModel.h"
+#include "polyhedralGravity/model/GravityModel.h"
 
 /**
  * Contains additional utility for working with the values of the polyhedrale Gravity Model.
diff --git a/test/model/PolyhedronBigTest.cpp b/test/model/PolyhedronBigTest.cpp
new file mode 100644
index 0000000..b6a7988
--- /dev/null
+++ b/test/model/PolyhedronBigTest.cpp
@@ -0,0 +1,101 @@
+#include "gtest/gtest.h"
+#include "gmock/gmock.h"
+
+#include <set>
+#include <chrono>
+#include <vector>
+#include <utility>
+#include <array>
+#include <algorithm>
+#include <random>
+#include <thrust/generate.h>
+#include <iostream>
+#include "polyhedralGravity/model/Polyhedron.h"
+
+/**
+ * Cheks that the detection of wrong faces works as intended.
+ */
+class PolyhedronBigTest : public ::testing::TestWithParam<std::set<size_t>> {
+
+protected:
+
+    const static inline polyhedralGravity::PolyhedralFiles FILENAMES{"resources/GravityModelBigTest.node", "resources/GravityModelBigTest.face"};
+
+    const static inline polyhedralGravity::Polyhedron CORRECT_POLYHEDRON{
+        FILENAMES, 1.0,
+        polyhedralGravity::NormalOrientation::OUTWARDS, polyhedralGravity::PolyhedronIntegrity::DISABLE
+    };
+
+    static constexpr size_t FACES_COUNT = 14744;
+    static constexpr size_t SET_SIZE = 100;
+    static constexpr size_t SET_NUMBER = 10;
+    static constexpr size_t SEED = 42;
+
+    /**
+     * Creates a polyhedron violating the OUTWARDS pointing normals constraint for the given face indices.
+     * @param violatinIndices the indices of the faces where the polyhedron's plane unit normals shall point INWARDS
+     * @return invalid polyhedron
+     */
+    polyhedralGravity::Polyhedron createViolatingPolyhedron(const std::set<size_t> &violatinIndices) {
+        using namespace polyhedralGravity;
+        std::vector<IndexArray3> violatingFaces = CORRECT_POLYHEDRON.getFaces();
+        for (const size_t index: violatinIndices) {
+            std::swap(violatingFaces[index][0], violatingFaces[index][1]);
+        }
+        return {
+            CORRECT_POLYHEDRON.getVertices(), violatingFaces,
+            1.0, NormalOrientation::OUTWARDS, PolyhedronIntegrity::DISABLE
+        };
+    }
+
+public:
+
+    // Function to generate a half sized set
+    static std::vector<std::set<size_t>> generateIndices() {
+        std::vector<std::set<size_t>> indexSets{};
+        indexSets.reserve(SET_NUMBER);
+
+        std::mt19937 engine{SEED};
+        std::uniform_int_distribution<size_t> dist{0, PolyhedronBigTest::FACES_COUNT - 1};
+
+        for (size_t i = 0; i < SET_NUMBER; ++i) {
+            std::set<size_t> generatedSet;
+            thrust::generate_n(std::inserter(generatedSet, generatedSet.end()), SET_SIZE, [&dist, &engine]() {return dist(engine);});
+            while(generatedSet.size() < SET_SIZE) {
+                generatedSet.insert(dist(engine));
+            }
+            indexSets.push_back(generatedSet);
+        }
+        return indexSets;
+    }
+
+};
+
+
+TEST_P(PolyhedronBigTest, BigPolyhedronFindWrongVertices) {
+    using namespace testing;
+    using namespace polyhedralGravity;
+
+    const std::set<size_t> expectedViolatingIndices = GetParam();
+    ASSERT_EQ(expectedViolatingIndices.size(), SET_SIZE);
+
+    Polyhedron invalidPolyhedron = createViolatingPolyhedron(expectedViolatingIndices);
+
+    auto start = std::chrono::high_resolution_clock::now();
+
+    const auto&[actualOrientation, actualViolatingIndices] = invalidPolyhedron.checkPlaneUnitNormalOrientation();
+
+    auto end = std::chrono::high_resolution_clock::now();
+    auto dur = end - start;
+    auto ms = std::chrono::duration_cast<std::chrono::milliseconds>(dur);
+    std::cout << "Measured time: " << ms.count() << " microseconds ";
+
+
+    // The orientation changes if violatingIndcies > len(faces)/2 --> So ensure that the size of the test-parameter is smaller than len(faces)/2
+    ASSERT_EQ(actualOrientation, NormalOrientation::OUTWARDS);
+    ASSERT_THAT(actualViolatingIndices, ContainerEq(expectedViolatingIndices));
+}
+
+INSTANTIATE_TEST_SUITE_P(NormalOrientationViolationTest, PolyhedronBigTest, ::testing::ValuesIn(PolyhedronBigTest::generateIndices()));
+
+
diff --git a/test/model/PolyhedronTest.cpp b/test/model/PolyhedronTest.cpp
new file mode 100644
index 0000000..496d801
--- /dev/null
+++ b/test/model/PolyhedronTest.cpp
@@ -0,0 +1,312 @@
+#include "gtest/gtest.h"
+#include "gmock/gmock.h"
+
+#include <array>
+#include "polyhedralGravity/model/Polyhedron.h"
+
+
+/**
+ * Contains Tests if the mesh sanity checks of the Polyhedron work
+ * Parameterization not possible since the patterns what and where the error exactly is varies leading
+ * to different assertions being required in the single test case.
+ */
+class PolyhedronTest : public ::testing::Test {
+
+protected:
+    const std::vector<polyhedralGravity::Array3> _cubeVertices{
+            {-1.0, -1.0, -1.0},
+            {1.0, -1.0, -1.0},
+            {1.0, 1.0, -1.0},
+            {-1.0, 1.0, -1.0},
+            {-1.0, -1.0, 1.0},
+            {1.0, -1.0, 1.0},
+            {1.0, 1.0, 1.0},
+            {-1.0, 1.0, 1.0}
+    };
+
+    const std::vector<polyhedralGravity::IndexArray3> _facesOutwards{
+            {1, 3, 2},
+            {0, 3, 1},
+            {0, 1, 5},
+            {0, 5, 4},
+            {0, 7, 3},
+            {0, 4, 7},
+            {1, 2, 6},
+            {1, 6, 5},
+            {2, 3, 6},
+            {3, 7, 6},
+            {4, 5, 6},
+            {4, 6, 7}
+    };
+
+
+    const std::vector<polyhedralGravity::IndexArray3> _facesInwards{
+            {3, 1, 2},
+            {3, 0, 1},
+            {1, 0, 5},
+            {5, 0, 4},
+            {7, 0, 3},
+            {4, 0, 7},
+            {2, 1, 6},
+            {6, 1, 5},
+            {3, 2, 6},
+            {7, 3, 6},
+            {5, 4, 6},
+            {6, 4, 7}
+    };
+
+    const std::vector<polyhedralGravity::IndexArray3> _facesOutwardsMajority{
+            {3, 1, 2},
+            {0, 3, 1},
+            {0, 1, 5},
+            {0, 5, 4},
+            {7, 0, 3},
+            {0, 4, 7},
+            {1, 2, 6},
+            {1, 6, 5},
+            {2, 3, 6},
+            {3, 7, 6},
+            {4, 5, 6},
+            {4, 6, 7}
+    };
+
+    const std::vector<polyhedralGravity::IndexArray3> _facesInwardsMajority{
+            {3, 1, 2},
+            {3, 0, 1},
+            {1, 0, 5},
+            {5, 0, 4},
+            {7, 0, 3},
+            {4, 0, 7},
+            {2, 1, 6},
+            {6, 1, 5},
+            {3, 2, 6},
+            {3, 7, 6},
+            {4, 5, 6},
+            {4, 6, 7}
+    };
+
+    const std::vector<polyhedralGravity::IndexArray3> _degeneratedFaces{
+            {1, 3, 2},
+            {0, 3, 1},
+            {0, 1, 5},
+            {0, 5, 4},
+            {7, 7, 3},
+            {0, 4, 7},
+            {1, 2, 6},
+            {1, 6, 5},
+            {2, 3, 6},
+            {3, 7, 6},
+            {4, 5, 6},
+            {4, 6, 7}
+    };
+
+
+    const std::vector<polyhedralGravity::Array3> _prismVertices{
+            {-20.0, 0.0, 25.0},
+            {0.0, 0.0, 25.0},
+            {0.0, 10.0, 25.0},
+            {-20.0, 10.0, 25.0},
+            {-20.0, 0.0, 15.0},
+            {0.0, 0.0, 15.0},
+            {0.0, 10.0, 15.0},
+            {-20.0, 10.0, 15.0}
+    };
+
+    const std::vector<polyhedralGravity::IndexArray3> _prismOutwards{
+            {0, 4, 5},
+            {0, 5, 1},
+            {0, 1, 3},
+            {1, 2, 3},
+            {1, 5, 6},
+            {1, 6, 2},
+            {0, 7, 4},
+            {0, 3, 7},
+            {4, 7, 5},
+            {5, 7, 6},
+            {2, 7, 3},
+            {2, 6, 7}
+    };
+
+    const std::vector<polyhedralGravity::IndexArray3> _prismInwards{
+
+            {4, 0, 5},
+            {5, 0, 1},
+            {1, 0, 3},
+            {2, 1, 3},
+            {5, 1, 6},
+            {6, 1, 2},
+            {7, 0, 4},
+            {3, 0, 7},
+            {7, 4, 5},
+            {7, 5, 6},
+            {7, 2, 3},
+            {6, 2, 7}
+    };
+};
+
+TEST_F(PolyhedronTest, CubeOutwardNormals) {
+    using namespace polyhedralGravity;
+    using namespace testing;
+    // Correct Set-Up, No-throw for all options
+    EXPECT_NO_THROW(Polyhedron(_cubeVertices, _facesOutwards, 1.0, NormalOrientation::OUTWARDS, PolyhedronIntegrity::DISABLE));
+    EXPECT_NO_THROW(Polyhedron(_cubeVertices, _facesOutwards, 1.0, NormalOrientation::OUTWARDS, PolyhedronIntegrity::AUTOMATIC));
+    EXPECT_NO_THROW(Polyhedron(_cubeVertices, _facesOutwards, 1.0, NormalOrientation::OUTWARDS, PolyhedronIntegrity::VERIFY));
+    EXPECT_NO_THROW(Polyhedron (_cubeVertices, _facesOutwards, 1.0, NormalOrientation::OUTWARDS, PolyhedronIntegrity::HEAL));
+
+    // Wrong Set-Up, Throws in case of AUTOMATIC and VERIFY, DISABLE and HEAL do not throw but respectivley ignore or repair
+    EXPECT_NO_THROW(Polyhedron(_cubeVertices, _facesOutwards, 1.0, NormalOrientation::INWARDS, PolyhedronIntegrity::DISABLE));
+    EXPECT_THROW(Polyhedron(_cubeVertices, _facesOutwards, 1.0, NormalOrientation::INWARDS, PolyhedronIntegrity::AUTOMATIC), std::invalid_argument);
+    EXPECT_THROW(Polyhedron(_cubeVertices, _facesOutwards, 1.0, NormalOrientation::INWARDS, PolyhedronIntegrity::VERIFY), std::invalid_argument);
+    EXPECT_NO_THROW(Polyhedron(_cubeVertices, _facesOutwards, 1.0, NormalOrientation::INWARDS, PolyhedronIntegrity::HEAL));
+
+    // Checking that the healing leads to the correct result (More than n/2 normals (actuall all) are outwards. Hence, only change of the Orientation)
+    Polyhedron healedPolyhedron(_cubeVertices, _facesOutwards, 1.0, NormalOrientation::INWARDS, PolyhedronIntegrity::HEAL);
+    EXPECT_EQ(healedPolyhedron.getOrientation(), NormalOrientation::OUTWARDS);
+    EXPECT_THAT(healedPolyhedron.getFaces(), ContainerEq(_facesOutwards));
+}
+
+TEST_F(PolyhedronTest, CubeInwardsNormals) {
+    using namespace polyhedralGravity;
+    using namespace testing;
+    // Correct Set-Up, No-throw for all options
+    EXPECT_NO_THROW(Polyhedron(_cubeVertices, _facesInwards, 1.0, NormalOrientation::INWARDS, PolyhedronIntegrity::DISABLE));
+    EXPECT_NO_THROW(Polyhedron(_cubeVertices, _facesInwards, 1.0, NormalOrientation::INWARDS, PolyhedronIntegrity::AUTOMATIC));
+    EXPECT_NO_THROW(Polyhedron(_cubeVertices, _facesInwards, 1.0, NormalOrientation::INWARDS, PolyhedronIntegrity::VERIFY));
+    EXPECT_NO_THROW(Polyhedron (_cubeVertices, _facesInwards, 1.0, NormalOrientation::INWARDS, PolyhedronIntegrity::HEAL));
+
+    // Wrong Set-Up, Throws in case of AUTOMATIC and VERIFY, DISABLE and HEAL do not throw but respectivley ignore or repair
+    EXPECT_NO_THROW(Polyhedron(_cubeVertices, _facesInwards, 1.0, NormalOrientation::OUTWARDS, PolyhedronIntegrity::DISABLE));
+    EXPECT_THROW(Polyhedron(_cubeVertices, _facesInwards, 1.0, NormalOrientation::OUTWARDS, PolyhedronIntegrity::AUTOMATIC), std::invalid_argument);
+    EXPECT_THROW(Polyhedron(_cubeVertices, _facesInwards, 1.0, NormalOrientation::OUTWARDS, PolyhedronIntegrity::VERIFY), std::invalid_argument);
+    EXPECT_NO_THROW(Polyhedron(_cubeVertices, _facesInwards, 1.0, NormalOrientation::OUTWARDS, PolyhedronIntegrity::HEAL));
+
+    // Checking that checkPlaneUnitNormalOrientation() returns the correct set of indices & majority orientation
+    Polyhedron intenionalDisabledCheckingPolyhedron(_cubeVertices, _facesOutwardsMajority, 1.0, NormalOrientation::INWARDS, PolyhedronIntegrity::DISABLE);
+    const auto &[majorityOrientation, violatingIndices] = intenionalDisabledCheckingPolyhedron.checkPlaneUnitNormalOrientation();
+    EXPECT_EQ(majorityOrientation, NormalOrientation::OUTWARDS);
+    EXPECT_THAT(violatingIndices, ContainerEq(std::set<size_t>({0, 4})));
+
+    // Checking that the healing leads to the correct result (More than n/2 normals (actuall all) are outwards. Hence, only change of the Orientation)
+    Polyhedron healedPolyhedron(_cubeVertices, _facesInwards, 1.0, NormalOrientation::OUTWARDS, PolyhedronIntegrity::HEAL);
+    EXPECT_EQ(healedPolyhedron.getOrientation(), NormalOrientation::INWARDS);
+    EXPECT_THAT(healedPolyhedron.getFaces(), ContainerEq(_facesInwards));
+}
+
+TEST_F(PolyhedronTest, CubeOutwardNormalsMajor) {
+    using namespace polyhedralGravity;
+    using namespace testing;
+    // The majority is correct. However, 2 faces have inward pointing normals (Index 0 and 4) --> Will throw but not in the healing scenario
+    EXPECT_NO_THROW(Polyhedron(_cubeVertices, _facesOutwardsMajority, 1.0, NormalOrientation::OUTWARDS, PolyhedronIntegrity::DISABLE));
+    EXPECT_THROW(Polyhedron(_cubeVertices, _facesOutwardsMajority, 1.0, NormalOrientation::OUTWARDS, PolyhedronIntegrity::AUTOMATIC), std::invalid_argument);
+    EXPECT_THROW(Polyhedron(_cubeVertices, _facesOutwardsMajority, 1.0, NormalOrientation::OUTWARDS, PolyhedronIntegrity::VERIFY), std::invalid_argument);
+    EXPECT_NO_THROW(Polyhedron (_cubeVertices, _facesOutwardsMajority, 1.0, NormalOrientation::OUTWARDS, PolyhedronIntegrity::HEAL));
+
+    // Wrong Set-Up, Throws in case of AUTOMATIC and VERIFY, DISABLE and HEAL do not throw but respectivley ignore or repair
+    EXPECT_NO_THROW(Polyhedron(_cubeVertices, _facesOutwardsMajority, 1.0, NormalOrientation::INWARDS, PolyhedronIntegrity::DISABLE));
+    EXPECT_THROW(Polyhedron(_cubeVertices, _facesOutwardsMajority, 1.0, NormalOrientation::INWARDS, PolyhedronIntegrity::AUTOMATIC), std::invalid_argument);
+    EXPECT_THROW(Polyhedron(_cubeVertices, _facesOutwardsMajority, 1.0, NormalOrientation::INWARDS, PolyhedronIntegrity::VERIFY), std::invalid_argument);
+    EXPECT_NO_THROW(Polyhedron(_cubeVertices, _facesOutwardsMajority, 1.0, NormalOrientation::INWARDS, PolyhedronIntegrity::HEAL));
+
+    // Checking that checkPlaneUnitNormalOrientation() returns the correct set of indices & majority orientation
+    Polyhedron intenionalDisabledCheckingPolyhedron(_cubeVertices, _facesOutwardsMajority, 1.0, NormalOrientation::INWARDS, PolyhedronIntegrity::DISABLE);
+    const auto &[majorityOrientation, violatingIndices] = intenionalDisabledCheckingPolyhedron.checkPlaneUnitNormalOrientation();
+    EXPECT_EQ(majorityOrientation, NormalOrientation::OUTWARDS);
+    EXPECT_THAT(violatingIndices, ContainerEq(std::set<size_t>({0, 4})));
+
+    // Checking that the healing leads to the correct result (More than n/2 normals are outwards. Here orientation & faces' vertex ordering changes)
+    Polyhedron healedPolyhedron(_cubeVertices, _facesOutwardsMajority, 1.0, NormalOrientation::INWARDS, PolyhedronIntegrity::HEAL);
+    EXPECT_EQ(healedPolyhedron.getOrientation(), NormalOrientation::OUTWARDS);
+    EXPECT_THAT(healedPolyhedron.getFaces(), ContainerEq(_facesOutwards));
+}
+
+TEST_F(PolyhedronTest, CubeOutwardInwardsMajor) {
+    using namespace polyhedralGravity;
+    using namespace testing;
+    // The majority is correct. However, 2 faces have inward pointing normals (Index 9, 10, and 11) --> Will throw but not in the healing scenario
+    EXPECT_NO_THROW(Polyhedron(_cubeVertices, _facesInwardsMajority, 1.0, NormalOrientation::INWARDS, PolyhedronIntegrity::DISABLE));
+    EXPECT_THROW(Polyhedron(_cubeVertices, _facesInwardsMajority, 1.0, NormalOrientation::INWARDS, PolyhedronIntegrity::AUTOMATIC), std::invalid_argument);
+    EXPECT_THROW(Polyhedron(_cubeVertices, _facesInwardsMajority, 1.0, NormalOrientation::INWARDS, PolyhedronIntegrity::VERIFY), std::invalid_argument);
+    EXPECT_NO_THROW(Polyhedron (_cubeVertices, _facesInwardsMajority, 1.0, NormalOrientation::INWARDS, PolyhedronIntegrity::HEAL));
+
+    // Wrong Set-Up, Throws in case of AUTOMATIC and VERIFY, DISABLE and HEAL do not throw but respectivley ignore or repair
+    EXPECT_NO_THROW(Polyhedron(_cubeVertices, _facesInwardsMajority, 1.0, NormalOrientation::OUTWARDS, PolyhedronIntegrity::DISABLE));
+    EXPECT_THROW(Polyhedron(_cubeVertices, _facesInwardsMajority, 1.0, NormalOrientation::OUTWARDS, PolyhedronIntegrity::AUTOMATIC), std::invalid_argument);
+    EXPECT_THROW(Polyhedron(_cubeVertices, _facesInwardsMajority, 1.0, NormalOrientation::OUTWARDS, PolyhedronIntegrity::VERIFY), std::invalid_argument);
+    EXPECT_NO_THROW(Polyhedron(_cubeVertices, _facesInwardsMajority, 1.0, NormalOrientation::OUTWARDS, PolyhedronIntegrity::HEAL));
+
+    // Checking that checkPlaneUnitNormalOrientation() returns the correct set of indices & majority orientation
+    Polyhedron intenionalDisabledCheckingPolyhedron(_cubeVertices, _facesInwardsMajority, 1.0, NormalOrientation::OUTWARDS, PolyhedronIntegrity::DISABLE);
+    const auto &[majorityOrientation, violatingIndices] = intenionalDisabledCheckingPolyhedron.checkPlaneUnitNormalOrientation();
+    EXPECT_EQ(majorityOrientation, NormalOrientation::INWARDS);
+    EXPECT_THAT(violatingIndices, ContainerEq(std::set<size_t>({9, 10, 11})));
+
+    // Checking that the healing leads to the correct result (More than n/2 normals are outwards. Here orientation & faces' vertex ordering changes)
+    Polyhedron healedPolyhedron(_cubeVertices, _facesInwardsMajority, 1.0, NormalOrientation::INWARDS, PolyhedronIntegrity::HEAL);
+    EXPECT_EQ(healedPolyhedron.getOrientation(), NormalOrientation::INWARDS);
+    EXPECT_THAT(healedPolyhedron.getFaces(), ContainerEq(_facesInwards));
+}
+
+TEST_F(PolyhedronTest, CubeDegenerated) {
+    using namespace polyhedralGravity;
+    using namespace testing;
+    // A degenarted mesh is not repaired, throw always expect DISABLE is checked
+    for (auto orientation: std::set({NormalOrientation::INWARDS, NormalOrientation::OUTWARDS})) {
+        EXPECT_NO_THROW(Polyhedron(_cubeVertices, _degeneratedFaces, 1.0, orientation, PolyhedronIntegrity::DISABLE));
+        EXPECT_THROW(Polyhedron(_cubeVertices, _degeneratedFaces, 1.0, orientation, PolyhedronIntegrity::AUTOMATIC), std::invalid_argument);
+        EXPECT_THROW(Polyhedron(_cubeVertices, _degeneratedFaces, 1.0, orientation, PolyhedronIntegrity::VERIFY), std::invalid_argument);
+        EXPECT_THROW(Polyhedron(_cubeVertices, _degeneratedFaces, 1.0, orientation, PolyhedronIntegrity::HEAL), std::invalid_argument);
+    }
+}
+
+TEST_F(PolyhedronTest, PrsimOutwards) {
+    using namespace polyhedralGravity;
+    using namespace testing;
+    // Correct Set-Up, No-throw for all options
+    EXPECT_NO_THROW(Polyhedron(_prismVertices, _prismOutwards, 1.0, NormalOrientation::OUTWARDS, PolyhedronIntegrity::DISABLE));
+    EXPECT_NO_THROW(Polyhedron(_prismVertices, _prismOutwards, 1.0, NormalOrientation::OUTWARDS, PolyhedronIntegrity::AUTOMATIC));
+    EXPECT_NO_THROW(Polyhedron(_prismVertices, _prismOutwards, 1.0, NormalOrientation::OUTWARDS, PolyhedronIntegrity::VERIFY));
+    EXPECT_NO_THROW(Polyhedron (_prismVertices, _prismOutwards, 1.0, NormalOrientation::OUTWARDS, PolyhedronIntegrity::HEAL));
+
+    // Wrong Set-Up, Throws in case of AUTOMATIC and VERIFY, DISABLE and HEAL do not throw but respectivley ignore or repair
+    EXPECT_NO_THROW(Polyhedron(_prismVertices, _prismOutwards, 1.0, NormalOrientation::INWARDS, PolyhedronIntegrity::DISABLE));
+    EXPECT_THROW(Polyhedron(_prismVertices, _prismOutwards, 1.0, NormalOrientation::INWARDS, PolyhedronIntegrity::AUTOMATIC), std::invalid_argument);
+    EXPECT_THROW(Polyhedron(_prismVertices, _prismOutwards, 1.0, NormalOrientation::INWARDS, PolyhedronIntegrity::VERIFY), std::invalid_argument);
+    EXPECT_NO_THROW(Polyhedron(_prismVertices, _prismOutwards, 1.0, NormalOrientation::INWARDS, PolyhedronIntegrity::HEAL));
+
+    // Checking that the healing leads to the correct result (More than n/2 normals (actuall all) are outwards. Hence, only change of the Orientation)
+    Polyhedron healedPolyhedron(_prismVertices, _prismOutwards, 1.0, NormalOrientation::INWARDS, PolyhedronIntegrity::HEAL);
+    EXPECT_EQ(healedPolyhedron.getOrientation(), NormalOrientation::OUTWARDS);
+    EXPECT_THAT(healedPolyhedron.getFaces(), ContainerEq(_prismOutwards));
+}
+
+TEST_F(PolyhedronTest, PrsimInwards) {
+    using namespace polyhedralGravity;
+    using namespace testing;
+    // Correct Set-Up, No-throw for all options
+    EXPECT_NO_THROW(Polyhedron(_prismVertices, _prismInwards, 1.0, NormalOrientation::INWARDS, PolyhedronIntegrity::DISABLE));
+    EXPECT_NO_THROW(Polyhedron(_prismVertices, _prismInwards, 1.0, NormalOrientation::INWARDS, PolyhedronIntegrity::AUTOMATIC));
+    EXPECT_NO_THROW(Polyhedron(_prismVertices, _prismInwards, 1.0, NormalOrientation::INWARDS, PolyhedronIntegrity::VERIFY));
+    EXPECT_NO_THROW(Polyhedron (_prismVertices, _prismInwards, 1.0, NormalOrientation::INWARDS, PolyhedronIntegrity::HEAL));
+
+    // Wrong Set-Up, Throws in case of AUTOMATIC and VERIFY, DISABLE and HEAL do not throw but respectivley ignore or repair
+    EXPECT_NO_THROW(Polyhedron(_prismVertices, _prismInwards, 1.0, NormalOrientation::OUTWARDS, PolyhedronIntegrity::DISABLE));
+    EXPECT_THROW(Polyhedron(_prismVertices, _prismInwards, 1.0, NormalOrientation::OUTWARDS, PolyhedronIntegrity::AUTOMATIC), std::invalid_argument);
+    EXPECT_THROW(Polyhedron(_prismVertices, _prismInwards, 1.0, NormalOrientation::OUTWARDS, PolyhedronIntegrity::VERIFY), std::invalid_argument);
+    EXPECT_NO_THROW(Polyhedron(_prismVertices, _prismInwards, 1.0, NormalOrientation::OUTWARDS, PolyhedronIntegrity::HEAL));
+
+    // Checking that the healing leads to the correct result (More than n/2 normals (actuall all) are outwards. Hence, only change of the Orientation)
+    Polyhedron healedPolyhedron(_prismVertices, _prismInwards, 1.0, NormalOrientation::OUTWARDS, PolyhedronIntegrity::HEAL);
+    EXPECT_EQ(healedPolyhedron.getOrientation(), NormalOrientation::INWARDS);
+    EXPECT_THAT(healedPolyhedron.getFaces(), ContainerEq(_prismInwards));
+}
+
+TEST_F(PolyhedronTest, CorrectBigPolyhedron) {
+    using namespace testing;
+    using namespace polyhedralGravity;
+    // All normals are pointing outwards, extensive Eros example
+    ASSERT_NO_THROW(Polyhedron(
+                std::vector<std::string>({"resources/GravityModelBigTest.node", "resources/GravityModelBigTest.face"}),
+                1.0, NormalOrientation::OUTWARDS, PolyhedronIntegrity::VERIFY)
+            );
+}
+
diff --git a/test/python/test_gravity_model.py b/test/python/test_gravity_model.py
new file mode 100644
index 0000000..18ca604
--- /dev/null
+++ b/test/python/test_gravity_model.py
@@ -0,0 +1,173 @@
+from typing import Tuple, List, Union
+from polyhedral_gravity import Polyhedron, GravityEvaluable, evaluate, PolyhedronIntegrity, NormalOrientation
+import numpy as np
+import pickle
+import pytest
+from pathlib import Path
+from functools import lru_cache
+
+CUBE_VERTICES = np.array([
+    [-1, -1, -1],
+    [1, -1, -1],
+    [1, 1, -1],
+    [-1, 1, -1],
+    [-1, -1, 1],
+    [1, -1, 1],
+    [1, 1, 1],
+    [-1, 1, 1]
+])
+
+CUBE_FACES = np.array([
+    [1, 3, 2],
+    [0, 3, 1],
+    [0, 1, 5],
+    [0, 5, 4],
+    [0, 7, 3],
+    [0, 4, 7],
+    [1, 2, 6],
+    [1, 6, 5],
+    [2, 3, 6],
+    [3, 7, 6],
+    [4, 5, 6],
+    [4, 6, 7]
+])
+
+CUBE_FACES_INVERTED = CUBE_FACES[:, [1, 0, 2]]
+
+DENSITY = 1.0
+
+CUBE_VERTICES_FILE = Path("test/resources/cube.node")
+CUBE_FACE_FILE = Path("test/resources/cube.face")
+
+
+@lru_cache(maxsize=10)
+def reference_solution(density: int) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
+    """Reads the cube reference solution from a file and returns it as a tuple"""
+    reference_file_path = Path(f"test/resources/analytic_cube_solution_density{int(density)}.txt")
+    reference = np.loadtxt(reference_file_path, skiprows=1)
+    points = reference[:, :3]
+    expected_potential = reference[:, 3].flatten()
+    expected_acceleration = reference[:, 4:]
+    return points, expected_potential, expected_acceleration
+
+
+@pytest.mark.parametrize(
+    "polyhedral_source,normal_orientation,density", [
+        ((CUBE_VERTICES, CUBE_FACES), NormalOrientation.OUTWARDS, 1.0),
+        ((CUBE_VERTICES, CUBE_FACES_INVERTED), NormalOrientation.INWARDS, 1.0),
+        ([str(CUBE_VERTICES_FILE), str(CUBE_FACE_FILE)], NormalOrientation.OUTWARDS, 1.0),
+        ((CUBE_VERTICES, CUBE_FACES), NormalOrientation.OUTWARDS, 42.0),
+        ((CUBE_VERTICES, CUBE_FACES_INVERTED), NormalOrientation.INWARDS, 42.0),
+        ([str(CUBE_VERTICES_FILE), str(CUBE_FACE_FILE)], NormalOrientation.OUTWARDS, 42.0)
+    ],
+    ids=["with_arrays_outwards01", "with_arrays_inwards01", "with_file_input01",
+         "with_arrays_outwards42", "with_arrays_inwards42", "with_file_input42"]
+)
+def test_polyhedral_gravity(
+        polyhedral_source: Union[Tuple[np.ndarray, np.ndarray], List[str]],
+        normal_orientation: NormalOrientation,
+        density: float
+) -> None:
+    """Checks that the evaluate function the correct results and
+    is callable with file/ array inputs.
+    """
+    points, expected_potential, expected_acceleration = reference_solution(density)
+    polyhedron = Polyhedron(
+        polyhedral_source=polyhedral_source,
+        density=DENSITY,
+        normal_orientation=normal_orientation,
+        integrity_check=PolyhedronIntegrity.VERIFY,
+    )
+    sol = evaluate(
+        polyhedron=polyhedron,
+        computation_points=points,
+        parallel=True,
+    )
+    potential = np.array([result[0] for result in sol])
+    acceleration = np.array([result[1] for result in sol])
+    np.testing.assert_array_almost_equal(potential, expected_potential)
+    np.testing.assert_array_almost_equal(acceleration, expected_acceleration)
+
+
+@pytest.mark.parametrize(
+    "polyhedral_source,normal_orientation,density", [
+        ((CUBE_VERTICES, CUBE_FACES), NormalOrientation.OUTWARDS, 1.0),
+        ((CUBE_VERTICES, CUBE_FACES_INVERTED), NormalOrientation.INWARDS, 1.0),
+        ([str(CUBE_VERTICES_FILE), str(CUBE_FACE_FILE)], NormalOrientation.OUTWARDS, 1.0),
+        ((CUBE_VERTICES, CUBE_FACES), NormalOrientation.OUTWARDS, 42.0),
+        ((CUBE_VERTICES, CUBE_FACES_INVERTED), NormalOrientation.INWARDS, 42.0),
+        ([str(CUBE_VERTICES_FILE), str(CUBE_FACE_FILE)], NormalOrientation.OUTWARDS, 42.0)
+    ],
+    ids=["with_arrays_outwards01", "with_arrays_inwards01", "with_file_input01",
+         "with_arrays_outwards42", "with_arrays_inwards42", "with_file_input42"]
+)
+def test_polyhedral_gravity_evaluable(
+        polyhedral_source: Union[Tuple[np.ndarray, np.ndarray], List[str]],
+        normal_orientation: NormalOrientation,
+        density: float
+) -> None:
+    """Checks that the evaluable produces the correct results and
+    is instantiable with file/ array inputs.
+    """
+    points, expected_potential, expected_acceleration = reference_solution(density)
+    polyhedron = Polyhedron(
+        polyhedral_source=polyhedral_source,
+        density=DENSITY,
+        normal_orientation=normal_orientation,
+        integrity_check=PolyhedronIntegrity.VERIFY,
+    )
+    evaluable = GravityEvaluable(polyhedron=polyhedron)
+    sol = evaluable(
+        computation_points=points,
+        parallel=True,
+    )
+    potential = np.array([result[0] for result in sol])
+    acceleration = np.array([result[1] for result in sol])
+    np.testing.assert_array_almost_equal(potential, expected_potential)
+    np.testing.assert_array_almost_equal(acceleration, expected_acceleration)
+
+
+@pytest.mark.parametrize(
+    "polyhedral_source,normal_orientation,density", [
+        ((CUBE_VERTICES, CUBE_FACES), NormalOrientation.OUTWARDS, 1.0),
+        ((CUBE_VERTICES, CUBE_FACES_INVERTED), NormalOrientation.INWARDS, 1.0),
+        ([str(CUBE_VERTICES_FILE), str(CUBE_FACE_FILE)], NormalOrientation.OUTWARDS, 1.0),
+        ((CUBE_VERTICES, CUBE_FACES), NormalOrientation.OUTWARDS, 42.0),
+        ((CUBE_VERTICES, CUBE_FACES_INVERTED), NormalOrientation.INWARDS, 42.0),
+        ([str(CUBE_VERTICES_FILE), str(CUBE_FACE_FILE)], NormalOrientation.OUTWARDS, 42.0)
+    ],
+    ids=["with_arrays_outwards01", "with_arrays_inwards01", "with_file_input01",
+         "with_arrays_outwards42", "with_arrays_inwards42", "with_file_input42"]
+)
+def test_polyhedral_evaluable_pickle(
+        polyhedral_source: Union[Tuple[np.ndarray, np.ndarray], List[str]],
+        normal_orientation: NormalOrientation,
+        density: float,
+        tmp_path: Path,) -> None:
+    """Tests that the evaluable can be pickled and unpicked and that the results
+    are still correct afterward (i.e. that the internal cache is correctly
+    pickled and unpicked as well).
+    """
+    points, expected_potential, expected_acceleration = reference_solution(density)
+    polyhedron = Polyhedron(
+        polyhedral_source=polyhedral_source,
+        density=DENSITY,
+        normal_orientation=normal_orientation,
+        integrity_check=PolyhedronIntegrity.DISABLE,
+    )
+    initial_evaluable = GravityEvaluable(polyhedron=polyhedron)
+    pickle_output = tmp_path.joinpath("evaluable.pk")
+    with open(pickle_output, "wb") as f:
+        pickle.dump(initial_evaluable, f, pickle.HIGHEST_PROTOCOL)
+
+    with open(pickle_output, "rb") as f:
+        read_evaluable = pickle.load(f)
+
+    sol = read_evaluable(
+        computation_points=points,
+        parallel=True,
+    )
+    potential = np.array([result[0] for result in sol])
+    acceleration = np.array([result[1] for result in sol])
+    np.testing.assert_array_almost_equal(potential, expected_potential)
+    np.testing.assert_array_almost_equal(acceleration, expected_acceleration)
diff --git a/test/python/test_polyhedral_gravity.py b/test/python/test_polyhedral_gravity.py
deleted file mode 100644
index bb7ec94..0000000
--- a/test/python/test_polyhedral_gravity.py
+++ /dev/null
@@ -1,130 +0,0 @@
-from typing import Tuple, List, Union
-import polyhedral_gravity
-import numpy as np
-import pickle
-import pytest
-from pathlib import Path
-
-CUBE_VERTICES = np.array([
-    [-1, -1, -1],
-    [1, -1, -1],
-    [1, 1, -1],
-    [-1, 1, -1],
-    [-1, -1, 1],
-    [1, -1, 1],
-    [1, 1, 1],
-    [-1, 1, 1]
-])
-
-CUBE_FACES = np.array([
-    [1, 3, 2],
-    [0, 3, 1],
-    [0, 1, 5],
-    [0, 5, 4],
-    [0, 7, 3],
-    [0, 4, 7],
-    [1, 2, 6],
-    [1, 6, 5],
-    [2, 3, 6],
-    [3, 7, 6],
-    [4, 5, 6],
-    [4, 6, 7]
-])
-
-DENSITY = 1.0
-
-CUBE_VERTICES_FILE = Path("test/resources/cube.node")
-CUBE_FACE_FILE = Path("test/resources/cube.face")
-
-
-@pytest.fixture(scope="session")
-def reference_solution() -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
-    """Reads the cube reference solution from a file and returns it as a tuple"""
-    reference_file_path = Path("test/resources/analytic_cube_solution_density1.txt")
-    reference = np.loadtxt(reference_file_path, skiprows=1)
-    points = reference[:, :3]
-    expected_potential = reference[:, 3].flatten()
-    expected_acceleration = reference[:, 4:]
-    return points, expected_potential, expected_acceleration
-
-
-@pytest.mark.parametrize(
-    "polyhedral_source", [(CUBE_VERTICES, CUBE_FACES), [str(CUBE_VERTICES_FILE), str(CUBE_FACE_FILE)]],
-    ids=["with_arrays", "with_file_input"]
-)
-def test_polyhedral_gravity(
-        polyhedral_source: Union[Tuple[np.ndarray, np.ndarray], List[str]],
-        reference_solution: Tuple[np.ndarray, np.ndarray, np.ndarray]) -> None:
-    """Checks that the evaluate function the correct results and
-    is callable with file/ array inputs.
-    """
-    points, expected_potential, expected_acceleration = reference_solution
-    sol = polyhedral_gravity.evaluate(
-        polyhedral_source=polyhedral_source,
-        density=DENSITY,
-        computation_points=points,
-        parallel=True
-    )
-    potential = np.array([result[0] for result in sol])
-    acceleration = np.array([result[1] for result in sol])
-    np.testing.assert_array_almost_equal(potential, expected_potential)
-    np.testing.assert_array_almost_equal(acceleration, expected_acceleration)
-
-
-@pytest.mark.parametrize(
-    "polyhedral_source", [(CUBE_VERTICES, CUBE_FACES), [str(CUBE_VERTICES_FILE), str(CUBE_FACE_FILE)]],
-    ids=["with_arrays", "with_file_input"]
-)
-def test_polyhedral_gravity_evaluable(
-        polyhedral_source: Union[Tuple[np.ndarray, np.ndarray], List[str]],
-        reference_solution: Tuple[np.ndarray, np.ndarray, np.ndarray]) -> None:
-    """Checks that the evaluable produces the correct results and
-    is instantiable with file/ array inputs.
-    """
-    points, expected_potential, expected_acceleration = reference_solution
-    evaluable = polyhedral_gravity.GravityEvaluable(
-        polyhedral_source=polyhedral_source,
-        density=DENSITY,
-    )
-    sol = evaluable(
-        computation_points=points,
-        parallel=True
-    )
-    potential = np.array([result[0] for result in sol])
-    acceleration = np.array([result[1] for result in sol])
-    np.testing.assert_array_almost_equal(potential, expected_potential)
-    np.testing.assert_array_almost_equal(acceleration, expected_acceleration)
-
-
-@pytest.mark.parametrize(
-    "polyhedral_source", [(CUBE_VERTICES, CUBE_FACES), [str(CUBE_VERTICES_FILE), str(CUBE_FACE_FILE)]],
-    ids=["with_arrays", "with_file_input"]
-)
-def test_polyhedral_evaluable_pickle(
-        polyhedral_source: Union[Tuple[np.ndarray, np.ndarray], List[str]],
-        reference_solution: Tuple[np.ndarray, np.ndarray, np.ndarray],
-        tmp_path: Path,) -> None:
-    """Tests that the evaluable can be pickled and unpicked and that the results
-    are still correct afterward (i.e. that the internal cache is correctly
-    pickled and unpicked as well).
-    """
-    points, expected_potential, expected_acceleration = reference_solution
-    initial_evaluable = polyhedral_gravity.GravityEvaluable(
-        polyhedral_source=polyhedral_source,
-        density=DENSITY
-    )
-    pickle_output = tmp_path.joinpath("evaluable.pk")
-    with open(pickle_output, "wb") as f:
-        pickle.dump(initial_evaluable, f, pickle.HIGHEST_PROTOCOL)
-
-    with open(pickle_output, "rb") as f:
-        read_evaluable = pickle.load(f)
-
-    sol = read_evaluable(
-        computation_points=points,
-        parallel=True
-    )
-    potential = np.array([result[0] for result in sol])
-    acceleration = np.array([result[1] for result in sol])
-    np.testing.assert_array_almost_equal(potential, expected_potential)
-    np.testing.assert_array_almost_equal(acceleration, expected_acceleration)
diff --git a/test/python/test_polyhedron.py b/test/python/test_polyhedron.py
new file mode 100644
index 0000000..d4c963f
--- /dev/null
+++ b/test/python/test_polyhedron.py
@@ -0,0 +1,88 @@
+from typing import Tuple, List, Union
+from polyhedral_gravity import Polyhedron, GravityEvaluable, evaluate, PolyhedronIntegrity, NormalOrientation
+import numpy as np
+import pickle
+import pytest
+from pathlib import Path
+
+CUBE_VERTICES = np.array([
+    [-1, -1, -1],
+    [1, -1, -1],
+    [1, 1, -1],
+    [-1, 1, -1],
+    [-1, -1, 1],
+    [1, -1, 1],
+    [1, 1, 1],
+    [-1, 1, 1]
+])
+
+CUBE_FACES_OUTWARDS = np.array([
+    [1, 3, 2],
+    [0, 3, 1],
+    [0, 1, 5],
+    [0, 5, 4],
+    [0, 7, 3],
+    [0, 4, 7],
+    [1, 2, 6],
+    [1, 6, 5],
+    [2, 3, 6],
+    [3, 7, 6],
+    [4, 5, 6],
+    [4, 6, 7]
+])
+
+CUBE_FACES_INWARDS = CUBE_FACES_OUTWARDS[:, [1, 0, 2]]
+CUBE_FACES_OUTWARDS_MAJOR = CUBE_FACES_OUTWARDS.copy()
+CUBE_FACES_OUTWARDS_MAJOR[0] = CUBE_FACES_OUTWARDS_MAJOR[0, [1, 0, 2]]
+
+
+CUBE_FACES_INWARDS_MAJOR = CUBE_FACES_INWARDS.copy()
+CUBE_FACES_INWARDS_MAJOR[2] = CUBE_FACES_INWARDS_MAJOR[2, [1, 0, 2]]
+CUBE_FACES_INWARDS_MAJOR[3] = CUBE_FACES_INWARDS_MAJOR[3, [1, 0, 2]]
+CUBE_FACES_INWARDS_MAJOR[5] = CUBE_FACES_INWARDS_MAJOR[5, [1, 0, 2]]
+
+DENSITY = 1.0
+
+CUBE_VERTICES_FILE = Path("test/resources/cube.node")
+CUBE_FACE_FILE = Path("test/resources/cube.face")
+
+@pytest.mark.parametrize(
+    "polyhedral_source,normal_orientation,violating_faces,expected_healed", [
+        ((CUBE_VERTICES, CUBE_FACES_OUTWARDS), NormalOrientation.OUTWARDS, [], CUBE_FACES_OUTWARDS),
+        ((CUBE_VERTICES, CUBE_FACES_INWARDS), NormalOrientation.INWARDS, [], CUBE_FACES_INWARDS),
+        ((CUBE_VERTICES, CUBE_FACES_OUTWARDS_MAJOR), NormalOrientation.OUTWARDS, [0], CUBE_FACES_OUTWARDS),
+        ((CUBE_VERTICES, CUBE_FACES_INWARDS_MAJOR), NormalOrientation.INWARDS, [2, 3, 5], CUBE_FACES_INWARDS),
+    ],
+    ids=["CubeOutwards", "CubeInwards", "CubeOutwardsMajor", "CubeInwardsMajor"]
+)
+def test_polyhedron(
+        polyhedral_source: Tuple[np.ndarray, np.ndarray],
+        normal_orientation: NormalOrientation,
+        violating_faces: List[int],
+        expected_healed: np.ndarray,
+) -> None:
+    inverse_orientation = NormalOrientation.OUTWARDS if NormalOrientation.INWARDS == normal_orientation else NormalOrientation.INWARDS
+
+    # The wrong orientation always leads to an error
+    with pytest.raises(ValueError):
+        Polyhedron(polyhedral_source, DENSITY, inverse_orientation, PolyhedronIntegrity.AUTOMATIC)
+        Polyhedron(polyhedral_source, DENSITY, inverse_orientation, PolyhedronIntegrity.VERIFY)
+    # At least one error will raise, no error won't raises
+    if len(violating_faces) == 0:
+        Polyhedron(polyhedral_source, DENSITY, normal_orientation, PolyhedronIntegrity.AUTOMATIC)
+        Polyhedron(polyhedral_source, DENSITY, normal_orientation, PolyhedronIntegrity.VERIFY)
+    else:
+        with pytest.raises(ValueError):
+            Polyhedron(polyhedral_source, DENSITY, normal_orientation, PolyhedronIntegrity.AUTOMATIC)
+            Polyhedron(polyhedral_source, DENSITY, normal_orientation, PolyhedronIntegrity.VERIFY)
+
+    # The healed version obeys to the majority ordering, Hence majority OUTWARDS --> only OUTWARDS
+    polyhedron_heal = Polyhedron(polyhedral_source, DENSITY, normal_orientation, PolyhedronIntegrity.HEAL)
+    healed_faces = np.array(polyhedron_heal.faces)
+    np.testing.assert_array_almost_equal(healed_faces, expected_healed)
+
+    # If disabled, there should not be any change - the only way to create an invalid polyhedron
+    polyhedron_non_modified = Polyhedron(polyhedral_source, DENSITY, normal_orientation, PolyhedronIntegrity.DISABLE)
+    non_modified_faces = np.array(polyhedron_non_modified.faces)
+    np.testing.assert_array_almost_equal(non_modified_faces, polyhedral_source[1])
+