Repository: ubermag/discretisedfield

Issue:
Bug in computing topological charge
When trying to compute the topological charge density in a plane (y,z) using the following code:
```python
import discretisedfield as df
import discretisedfield.tools as dft
p1 = (0, 0, 0)
p2 = (10, 10, 10)
cell = (2, 2, 2)
mesh = df.Mesh(p1=p1, p2=p2, cell=cell)
f = df.Field(mesh, nvdim=3, value=(1, 1, 1))
dft.topological_charge_density(f.sel('x'))
```
we get the following error
```python
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[26], line 8
      6 mesh = df.Mesh(p1=p1, p2=p2, cell=cell)
      7 f = df.Field(mesh, nvdim=3, value=(1, 1, 1))
----> 8 dft.topological_charge_density(f.sel('x'))
      9 #dft.topological_charge_density(f.sel('z'), method='berg-luescher')

File ~/ubermag_devtools/repos/discretisedfield/discretisedfield/tools/tools.py:124, in topological_charge_density(field, method)
    122 axis1 = field.vdim_mapping[field.vdims[0]]
    123 axis2 = field.vdim_mapping[field.vdims[1]]
--> 124 axis1_idx = field.mesh.region._dim2index(axis1)
    125 axis2_idx = field.mesh.region._dim2index(axis2)
    127 if method == "continuous":

File ~/ubermag_devtools/repos/discretisedfield/discretisedfield/region.py:284, in Region._dim2index(self, dim)
    282     return self.dims.index(dim)
    283 except ValueError:
--> 284     raise ValueError(f"'{dim}' not in region.dims={self.dims}.") from None

ValueError: 'x' not in region.dims=('y', 'z').
```

Discussion:

Patch:
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 30feaa81..325dfff0 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -29,7 +29,7 @@ repos:
     additional_dependencies: [flake8-rst-docstrings]  #, flake8-docstrings]
 
 - repo: https://github.com/psf/black
-  rev: 23.1.0
+  rev: 22.12.0
   hooks:
   - id: black-jupyter
 
diff --git a/discretisedfield/tools/tools.py b/discretisedfield/tools/tools.py
index bb68609c..27ba98c4 100644
--- a/discretisedfield/tools/tools.py
+++ b/discretisedfield/tools/tools.py
@@ -102,8 +102,9 @@ def topological_charge_density(field, /, method="continuous"):
 
     """
     if field.nvdim != 3:
-        msg = f"Cannot compute topological charge density for {field.nvdim=} field."
-        raise ValueError(msg)
+        raise ValueError(
+            f"Cannot compute topological charge density for {field.nvdim=} field."
+        )
 
     if field.mesh.region.ndim != 2:
         raise ValueError(
@@ -111,52 +112,27 @@ def topological_charge_density(field, /, method="continuous"):
             f" spatial dimensions, not {field.mesh.region.ndim=}."
         )
 
-    if method not in ["continuous", "berg-luescher"]:
-        msg = "Method can be either continuous or berg-luescher"
-        raise ValueError(msg)
-
     of = field.orientation  # unit field - orientation field
 
-    # Spatial axes associated with the first and second vector
-    # dimensions.
-    axis1 = field.vdim_mapping[field.vdims[0]]
-    axis2 = field.vdim_mapping[field.vdims[1]]
-    axis1_idx = field.mesh.region._dim2index(axis1)
-    axis2_idx = field.mesh.region._dim2index(axis2)
-
     if method == "continuous":
+        axis1 = field.mesh.region.dims[0]
+        axis2 = field.mesh.region.dims[1]
         return 1 / (4 * np.pi) * of.dot(of.diff(axis1).cross(of.diff(axis2)))
 
     elif method == "berg-luescher":
         q = df.Field(field.mesh, nvdim=1)
 
         # Area of a single triangle
-        area = 0.5 * field.mesh.cell[axis1_idx] * field.mesh.cell[axis2_idx]
+        area = 0.5 * field.mesh.cell[0] * field.mesh.cell[1]
 
-        for i, j in itertools.product(
-            range(of.mesh.n[axis1_idx]), range(of.mesh.n[axis2_idx])
-        ):
-            index = dfu.assemble_index(0, 2, {axis1_idx: i, axis2_idx: j})
-            v0 = of.array[index]
+        for i, j in itertools.product(range(of.mesh.n[0]), range(of.mesh.n[1])):
+            v0 = of.array[i, j]
 
             # Extract 4 neighbouring vectors (if they exist)
-            v1 = v2 = v3 = v4 = None
-            if i + 1 < of.mesh.n[axis1_idx]:
-                v1 = of.array[
-                    dfu.assemble_index(0, 2, {axis1_idx: i + 1, axis2_idx: j})
-                ]
-            if j + 1 < of.mesh.n[axis2_idx]:
-                v2 = of.array[
-                    dfu.assemble_index(0, 2, {axis1_idx: i, axis2_idx: j + 1})
-                ]
-            if i - 1 >= 0:
-                v3 = of.array[
-                    dfu.assemble_index(0, 2, {axis1_idx: i - 1, axis2_idx: j})
-                ]
-            if j - 1 >= 0:
-                v4 = of.array[
-                    dfu.assemble_index(0, 2, {axis1_idx: i, axis2_idx: j - 1})
-                ]
+            v1 = of.array[i + 1, j] if i + 1 < of.mesh.n[0] else None
+            v2 = of.array[i, j + 1] if j + 1 < of.mesh.n[1] else None
+            v3 = of.array[i - 1, j] if i - 1 >= 0 else None
+            v4 = of.array[i, j - 1] if j - 1 >= 0 else None
 
             charge = 0
             triangle_count = 0
@@ -178,13 +154,17 @@ def topological_charge_density(field, /, method="continuous"):
                 charge += dfu.bergluescher_angle(v0, v4, v1)
 
             if triangle_count > 0:
-                q.array[index] = charge / (area * triangle_count)
-            else:
-                # If the cell has no neighbouring cells
-                q.array[index] = 0
+                q.array[i, j] = charge / (area * triangle_count)
+            else:  # the cell has no neighbouring cells
+                q.array[i, j] = 0
 
         return q
 
+    else:
+        raise ValueError(
+            f"'method' can be either 'continuous' or 'berg-luescher', not '{method}'."
+        )
+
 
 def topological_charge(field, /, method="continuous", absolute=False):
     """Topological charge.
@@ -277,8 +257,7 @@ def topological_charge(field, /, method="continuous", absolute=False):
 
     """
     if field.nvdim != 3:
-        msg = f"Cannot compute topological charge for {field.nvdim=} field."
-        raise ValueError(msg)
+        raise ValueError(f"Cannot compute topological charge for {field.nvdim=} field.")
 
     if field.mesh.region.ndim != 2:
         raise ValueError(
@@ -286,10 +265,6 @@ def topological_charge(field, /, method="continuous", absolute=False):
             f" spatial dimensions, not {field.mesh.region.ndim=}."
         )
 
-    if method not in ["continuous", "berg-luescher"]:
-        msg = "Method can be either continuous or berg-luescher"
-        raise ValueError(msg)
-
     if method == "continuous":
         q = topological_charge_density(field, method=method)
         if absolute:
@@ -298,24 +273,14 @@ def topological_charge(field, /, method="continuous", absolute=False):
             return float(q.integrate())
 
     elif method == "berg-luescher":
-        # Spatial axes associated with the first and second vector
-        # dimensions.
-        axis1 = field.vdim_mapping[field.vdims[0]]
-        axis2 = field.vdim_mapping[field.vdims[1]]
-        axis1_idx = field.mesh.region._dim2index(axis1)
-        axis2_idx = field.mesh.region._dim2index(axis2)
         of = field.orientation
 
         topological_charge = 0
-        for i, j in itertools.product(
-            range(of.mesh.n[axis1_idx] - 1), range(of.mesh.n[axis2_idx] - 1)
-        ):
-            v1 = of.array[dfu.assemble_index(0, 2, {axis1_idx: i, axis2_idx: j})]
-            v2 = of.array[dfu.assemble_index(0, 2, {axis1_idx: i + 1, axis2_idx: j})]
-            v3 = of.array[
-                dfu.assemble_index(0, 2, {axis1_idx: i + 1, axis2_idx: j + 1})
-            ]
-            v4 = of.array[dfu.assemble_index(0, 2, {axis1_idx: i, axis2_idx: j + 1})]
+        for i, j in itertools.product(range(of.mesh.n[0] - 1), range(of.mesh.n[1] - 1)):
+            v1 = of.array[i, j]
+            v2 = of.array[i + 1, j]
+            v3 = of.array[i + 1, j + 1]
+            v4 = of.array[i, j + 1]
 
             triangle1 = dfu.bergluescher_angle(v1, v2, v4)
             triangle2 = dfu.bergluescher_angle(v2, v3, v4)
@@ -328,6 +293,11 @@ def topological_charge(field, /, method="continuous", absolute=False):
 
         return topological_charge
 
+    else:
+        raise ValueError(
+            f"'method' can be either 'continuous' or 'berg-luescher', not '{method}'."
+        )
+
 
 def emergent_magnetic_field(field):
     r"""Emergent magnetic field.


Test patch:
diff --git a/discretisedfield/tests/test_field.py b/discretisedfield/tests/test_field.py
index dd7bb49c..844b8219 100644
--- a/discretisedfield/tests/test_field.py
+++ b/discretisedfield/tests/test_field.py
@@ -128,6 +128,7 @@ def test_init_scalar_valid_args(valid_mesh, value, dtype):
 
 @pytest.mark.parametrize("value, dtype", iters + vfuncs)
 def test_init_vector_valid_args(valid_mesh, value, dtype):
+
     f = df.Field(valid_mesh, nvdim=3, value=value, dtype=dtype)
     assert isinstance(f, df.Field)
 
diff --git a/discretisedfield/tests/test_mesh.py b/discretisedfield/tests/test_mesh.py
index a3608724..3f08ba2d 100644
--- a/discretisedfield/tests/test_mesh.py
+++ b/discretisedfield/tests/test_mesh.py
@@ -1726,6 +1726,7 @@ def test_sel_subregions():
 
 
 def test_fftn_mesh():
+
     # 3D mesh with even number of cells
     p1 = (-10, -10, -5)
     p2 = (10, 10, 5)
diff --git a/discretisedfield/tests/test_tools.py b/discretisedfield/tests/test_tools.py
index 932552e6..0478c2be 100644
--- a/discretisedfield/tests/test_tools.py
+++ b/discretisedfield/tests/test_tools.py
@@ -18,12 +18,18 @@ def test_topological_charge(method):
     # -> Q(f) = 0
     f = df.Field(mesh, nvdim=3, value=(0, 0, 0))
 
-    q = dft.topological_charge_density(f.sel("z"), method=method)
-
-    assert q.nvdim == 1
-    assert q.mean() == 0
-    for absolute in [True, False]:
-        assert dft.topological_charge(f.sel("z"), method=method, absolute=absolute) == 0
+    for normal_direction in "xyz":
+        q = dft.topological_charge_density(f.sel(normal_direction), method=method)
+        assert q.nvdim == 1
+        assert np.allclose(q.mean(), 0)
+
+        for absolute in [True, False]:
+            assert np.allclose(
+                dft.topological_charge(
+                    f.sel(normal_direction), method=method, absolute=absolute
+                ),
+                0,
+            )
 
     # Skyrmion (with PBC) from a file
     test_filename = os.path.join(