Check validity of the polyhedra before computing gravity fields?

[santisoler]

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:

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:

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.

[schuhmaj]

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

[santisoler]

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:

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?

[schuhmaj]

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)

[schuhmaj]

I will add these features in the next days, but in a new PR other than #31

[santisoler]

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.