[JOSS Review] Minimal Python example doesn't run as expected

[santisoler]

The issue

While trying to run the minimal Python example showcased in the README.md, it runs into an error:

polyhedral-gravity-model/README.md

Lines 76 to 106 in ebe47ec

	The following example shows how to use the python interface to compute the gravity
	around a cube:
	```python
	import numpy as np
	import polyhedral_gravity
	# We define the cube as a polyhedron with 8 vertices and 12 triangular faces
	# 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]]
	)
	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]])
	```
	The simplest way to compute the gravity is to use the `evaluate` function:
	```python
	potential, acceleration, tensor = polyhedral_gravity.evaluate(
	polyhedral_source=(cube_vertices, cube_faces),
	density=cube_density,
	computation_points=computation_point
	parallel=True
	)
	```

Here I paste a copy of it that includes the fix in #20:

import numpy as np
import polyhedral_gravity

# We define the cube as a polyhedron with 8 vertices and 12 triangular faces
# 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],
    ]
)
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]])

potential, acceleration, tensor = polyhedral_gravity.evaluate(
    polyhedral_source=(cube_vertices, cube_faces),
    density=cube_density,
    computation_points=computation_point,
    parallel=True,
)

Traceback (most recent call last):
  File ".../example_01.py", line 37, in <module>
    potential, acceleration, tensor = polyhedral_gravity.evaluate(
    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
ValueError: not enough values to unpack (expected 3, got 1)

The output of polyhedral_gravity.evaluate() is a list with a single element, so it cannot be unpacked into three variables.

The solution

I would recommend to fix the example so it runs properly, unless this bug is highlighting some issue in the source code of evaluate() that should be tackled.

Update 2024-03-07

I noticed that the example using GravityEvaluable falls in the same error as the one that uses the evaluate() function.

Moreover, the same examples are replicated in the docs/quick_start_python.rst and they won't run as expected either.

I add more points to the previous recommendation:

• Fix the GravityEvaluable example in README.md and in docs/quick_start_python.rst.
• Fix the evaluate() example in docs/quick_start_python.rst.

---

This issue is part of the JOSS review: openjournals/joss-reviews#6384

[santisoler]

I just noticed by comparing the example with the Jupyter Notebook that if we define the computation_point as a 1d array (computation_point = np.array([0, 0, 0])) the example works fine.

This version of it doesn't fail on my end:

import numpy as np
import polyhedral_gravity

# We define the cube as a polyhedron with 8 vertices and 12 triangular faces
# 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],
    ]
)
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])

potential, acceleration, tensor = polyhedral_gravity.evaluate(
    polyhedral_source=(cube_vertices, cube_faces),
    density=cube_density,
    computation_points=computation_point,
    parallel=True,
)