bounding box mistake found. Problem was mesh cleaned. Problem with split...

bounding box mistake found. Problem was mesh cleaned. Problem with split remains (see test_mesh3D.py main TODO)

tests/test_bounding_box.py

 import unittest
 import tfields
-# import numpy as np
+import numpy as np
 
 
 class BoundingBox_Test(unittest.TestCase):
     def setUp(self):
-        self._mesh = tfields.Mesh3D.grid((5.6, 6.2, 3),
-                                         (-0.25, 0.25, 4),
-                                         (-1, 1, 10))
+        self.mesh = tfields.Mesh3D.grid((5.6, 6.2, 3), (-0.25, 0.25, 4), (-1, 1, 10))
 
-        self._cuts = {'x': [5.7, 6.1],
-                      'y': [-0.2, 0, 0.2],
-                      'z': [-0.5, 0.5]}
+        self.cuts = {"x": [5.7, 6.1], "y": [-0.2, 0, 0.2], "z": [-0.5, 0.5]}
 
     def test_tree(self):
         # test already in doctests.
-        tree = tfields.bounding_box.Node(self._mesh,
-                                         self._cuts,
-                                         at_intersection='keep')
+        tree = tfields.bounding_box.Node(self.mesh, self.cuts, at_intersection="keep")
         leaves = tree.leaves()
         leaves = tfields.bounding_box.Node.sort_leaves(leaves)
         meshes = [leaf.mesh for leaf in leaves]  # NOQA
@@ -25,21 +19,39 @@ class BoundingBox_Test(unittest.TestCase):
         special_leaf = tree.find_leaf([5.65, -0.21, 0])  # NOQA
 
 
-class Searcher_Test(unittest.TestCase):
+class Searcher_Check(object):
     def setUp(self):
-        self._mesh = tfields.Mesh3D.grid((0, 1, 2), (1, 2, 2), (2, 3, 2))
+        self.mesh = None
+        self.points = None
+        self.delta = None
+        self.n_sections = None
 
-    # not yet working again
-    # def test_tree(self):
-    #     tree = tfields.bounding_box.Searcher(self._mesh, n_sections=[5, 5, 5])
-    #     points = tfields.Tensors([[0.5, 1, 2.1],
-    #                               [0.5, 0, 0],
-    #                               [0.5, 2, 2.1],
-    #                               [0.5, 1.5, 2.5]])
-    #     box_res = tree.in_faces(points, delta=0.0001)
-    #     usual_res = self._mesh.in_faces(points, delta=0.0001)
-    #     self.assertTrue(np.array_equal(box_res, usual_res))
+    @property
+    def usual_res(self):
+        if not hasattr(self, "_usual_res"):
+            self._usual_res = self.mesh.in_faces(self.points, delta=self.delta)
+        return self._usual_res
+
+    def test_tree(self):
+        tree = tfields.bounding_box.Searcher(self.mesh, n_sections=None)
+        box_res = tree.in_faces(self.points, delta=self.delta)
+        self.assertTrue(np.array_equal(box_res, self.usual_res))
+
+    def test_tree_n_sections(self):
+        tree = tfields.bounding_box.Searcher(self.mesh, n_sections=self.n_sections)
+        box_res = tree.in_faces(self.points, delta=self.delta)
+        self.assertTrue(np.array_equal(box_res, self.usual_res))
+
+
+class Searcher_Coarse_Grid_2D_Test(Searcher_Check, unittest.TestCase):
+    def setUp(self):
+        self.mesh = tfields.Mesh3D.grid((0, 1, 2), (1, 2, 2), (2, 3, 2))
+        self.points = tfields.Tensors(
+            [[0.5, 1, 2.1], [0.5, 0, 0], [0.5, 2, 2.1], [0.5, 1.5, 2.5]]
+        )
+        self.delta = 0.0001
+        self.n_sections = [5, 5, 5]
 
 
-if __name__ == '__main__':
+if __name__ == "__main__":
     unittest.main()

tests/test_mesh3D.py

@@ -6,15 +6,17 @@ import sympy  # NOQA: F401
 import os
 import sys
 from tests.test_core import TensorMaps_Check
+
 THIS_DIR = os.path.dirname(
-    os.path.realpath(os.path.join(os.getcwd(), os.path.expanduser(__file__))))
+    os.path.realpath(os.path.join(os.getcwd(), os.path.expanduser(__file__)))
+)
 sys.path.append(os.path.normpath(os.path.join(THIS_DIR)))
 
 
 class Mesh3D_Check(TensorMaps_Check):
     def test_cut_split(self):
-        x, y, z = sympy.symbols('x y z')
-        self._inst.cut(x + 1./100*y > 0, at_intersection='split')
+        x, y, z = sympy.symbols("x y z")
+        self._inst.cut(x + 1.0 / 100 * y > 0, at_intersection="split")
 
     def test_triangle(self):
         tri = self._inst.triangles()
@@ -24,28 +26,48 @@ class Mesh3D_Check(TensorMaps_Check):
         self.assertTrue(tri.equal(tri_2))
 
     def test_save_obj(self):
-        out_file = NamedTemporaryFile(suffix='.obj')
+        out_file = NamedTemporaryFile(suffix=".obj")
         self._inst.save(out_file.name)
         _ = out_file.seek(0)  # this is only necessary in the test
         load_inst = type(self._inst).load(out_file.name)
-        print(self._inst, load_inst)
         self.demand_equal(load_inst)
 
 
-class Square_Test(Mesh3D_Check, unittest.TestCase):
+class Square_Sparse_Test(Mesh3D_Check, unittest.TestCase):
     def setUp(self):
         self._inst = tfields.Mesh3D.plane((0, 1, 2j), (0, 1, 2j), (0, 0, 1j))
 
+    def test_cut_keep(self):
+        x = sympy.symbols("x")
+        cut_inst = self._inst.cut(x < 0.5, at_intersection="keep")
+        self.assertTrue(self._inst.equal(cut_inst))
+        cut_inst = self._inst.cut(x > 0.5, at_intersection="keep")
+        self.assertTrue(self._inst.equal(cut_inst))
+
+
+class Square_Dense_Test(Mesh3D_Check, unittest.TestCase):
+    def setUp(self):
+        self._inst = tfields.Mesh3D.plane((0, 1, 5j), (0, 1, 5j), (0, 0, 1j))
+
+    def test_cut_keep(self):
+        x = sympy.symbols("x")
+        cut_inst = self._inst.cut(x < 0.9, at_intersection="keep")
+        self.assertTrue(self._inst.equal(cut_inst))
+        self.assertEqual(len(cut_inst), 25)
+        cut_inst = self._inst.cut(x > 0.9, at_intersection="keep")
+        self.assertEqual(len(cut_inst), 10)
+
 
 class Sphere_Test(Mesh3D_Check, unittest.TestCase):
     def setUp(self):
         basis_points = 4
         self._inst = tfields.Mesh3D.grid(
-                (1, 1, 1),
-                (-np.pi, np.pi, basis_points),
-                (-np.pi / 2, np.pi / 2, basis_points),
-                coord_sys='spherical')
-        self._inst.transform('cartesian')
+            (1, 1, 1),
+            (-np.pi, np.pi, basis_points),
+            (-np.pi / 2, np.pi / 2, basis_points),
+            coord_sys="spherical",
+        )
+        self._inst.transform("cartesian")
         self._inst[:, 1] += 2
         clean = self._inst.cleaned()
         # self.demand_equal(clean)
@@ -54,9 +76,17 @@ class Sphere_Test(Mesh3D_Check, unittest.TestCase):
 
 class IO_Stl_test(unittest.TestCase):  # no Mesh3D_Check for speed
     def setUp(self):
-        self._inst = tfields.Mesh3D.load(os.path.join(THIS_DIR,
-                                                      '../data/baffle.stl'))
+        self._inst = tfields.Mesh3D.load(os.path.join(THIS_DIR, "../data/baffle.stl"))
 
 
-if __name__ == '__main__':
+if __name__ == "__main__":
     unittest.main()
+    # TODO!
+    import rna
+    from sympy.abc import x
+
+    kwargs = dict(dim=3, edgecolor="k")
+    mesh = tfields.Mesh3D.grid((0, 1, 2), (1, 2, 2), (2, 3, 2))
+    cuts = mesh.cut(x < 0.2, at_intersection="split")
+    cuts.plot(**kwargs)
+    rna.plotting.show()

tfields/bounding_box.py

@@ -17,6 +17,14 @@ class Node(object):
         cut_expr: Cut that determines the seperation in left and right node
         cuts: List of cuts for the children nodes
 
+    Attrs:
+        parent (Node)
+        remaining_cuts (dict): key specifies dimension, value the cuts that
+            are still not done
+        cut_expr (dict): part of parents remaining_cuts. The dimension defines
+            what is meant by left and right
+
+
     Examples:
         >>> import tfields
         >>> mesh = tfields.Mesh3D.grid((5.6, 6.2, 3),
@@ -37,27 +45,35 @@ class Node(object):
         >>> special_leaf = tree.find_leaf([5.65, -0.21, 0])
 
     """
-    def __init__(self, mesh, cuts,
-                 coord_sys=None, at_intersection="split", delta=0.,
-                 parent=None, box=None, internal_template=None, cut_expr=None):
+
+    def __init__(
+        self,
+        mesh,
+        cuts,
+        coord_sys=None,
+        at_intersection="split",
+        delta=0.0,
+        parent=None,
+        box=None,
+        internal_template=None,
+        cut_expr=None,
+    ):
         self.parent = parent
         # initialize
         self.mesh = copy.deepcopy(mesh)
         if self.is_root():
             cuts = copy.deepcopy(cuts)  # dicts are mutable
         self.remaining_cuts = cuts
-        log = logging.getLogger('Node')
-        log.debug(cuts)
+        logging.debug(cuts)
 
         self.delta = delta
         if box is None:
             vertices = np.array(self.mesh)
-            self.box = {'x': [min(vertices[:, 0]) - delta,
-                              max(vertices[:, 0]) + delta],
-                        'y': [min(vertices[:, 1]) - delta,
-                              max(vertices[:, 1]) + delta],
-                        'z': [min(vertices[:, 2]) - delta,
-                              max(vertices[:, 2]) + delta]}
+            self.box = {
+                "x": [min(vertices[:, 0]) - delta, max(vertices[:, 0]) + delta],
+                "y": [min(vertices[:, 1]) - delta, max(vertices[:, 1]) + delta],
+                "z": [min(vertices[:, 2]) - delta, max(vertices[:, 2]) + delta],
+            }
         else:
             self.box = box
         self.left = None
@@ -100,7 +116,7 @@ class Node(object):
 
     def in_box(self, point):
         x, y, z = point
-        for key in ['x', 'y', 'z']:
+        for key in ["x", "y", "z"]:
             value = locals()[key]
             if value < self.box[key][0] or self.box[key][1] < value:
                 return False
@@ -117,23 +133,29 @@ class Node(object):
         """
         sorting the leaves first in x, then y, then z direction
         """
-        sorted_leaves = sorted(leaves_list,
-                               key=lambda x: (x.box['x'][1], x.box['y'][1], x.box['z'][1]))
+        sorted_leaves = sorted(
+            leaves_list, key=lambda x: (x.box["x"][1], x.box["y"][1], x.box["z"][1])
+        )
         return sorted_leaves
 
     def _trim_to_box(self):
         # 6 cuts to remove outer part of the box
-        x, y, z = sympy.symbols('x y z')
+        x, y, z = sympy.symbols("x y z")
         eps = 0.0000000001
-        x_cut = (float(self.box['x'][0] - eps) <= x) & (x <= float(self.box['x'][1] + eps))
-        y_cut = (float(self.box['y'][0] - eps) <= y) & (y <= float(self.box['y'][1] + eps))
-        z_cut = (float(self.box['z'][0] - eps) <= z) & (z <= float(self.box['z'][1] + eps))
+        x_cut = (float(self.box["x"][0] - eps) <= x) & (
+            x <= float(self.box["x"][1] + eps)
+        )
+        y_cut = (float(self.box["y"][0] - eps) <= y) & (
+            y <= float(self.box["y"][1] + eps)
+        )
+        z_cut = (float(self.box["z"][0] - eps) <= z) & (
+            z <= float(self.box["z"][1] + eps)
+        )
         section_cut = x_cut & y_cut & z_cut
 
         self.mesh, self._internal_template = self.mesh.cut(
-            section_cut,
-            at_intersection=self.at_intersection,
-            return_template=True)
+            section_cut, at_intersection=self.at_intersection, return_template=True
+        )
 
     def leaves(self):
         """
@@ -182,26 +204,28 @@ class Node(object):
 
     def _split(self):
         """
-        Split the node, if there is no cut_expr set and remaing cuts exist.
+        Split the node in two new nodes, if there is no cut_expr set and
+        remaing cuts exist.
         """
         if self.cut_expr is None and self.remaining_cuts is None:
-            raise RuntimeError("Cannot split the mesh without cut_expr and"
-                               "remaining_cuts")
+            raise RuntimeError(
+                "Cannot split the mesh without cut_expr and" "remaining_cuts"
+            )
         else:
             # create cut expression
-            x, y, z = sympy.symbols('x y z')
-            if 'x' in self.cut_expr:
-                left_cut_expression = x <= self.cut_expr['x']
-                right_cut_expression = x >= self.cut_expr['x']
-                key = 'x'
-            elif 'y' in self.cut_expr:
-                left_cut_expression = y <= self.cut_expr['y']
-                right_cut_expression = y >= self.cut_expr['y']
-                key = 'y'
-            elif 'z' in self.cut_expr:
-                left_cut_expression = z <= self.cut_expr['z']
-                right_cut_expression = z >= self.cut_expr['z']
-                key = 'z'
+            x, y, z = sympy.symbols("x y z")
+            if "x" in self.cut_expr:
+                left_cut_expression = x <= self.cut_expr["x"]
+                right_cut_expression = x >= self.cut_expr["x"]
+                key = "x"
+            elif "y" in self.cut_expr:
+                left_cut_expression = y <= self.cut_expr["y"]
+                right_cut_expression = y >= self.cut_expr["y"]
+                key = "y"
+            elif "z" in self.cut_expr:
+                left_cut_expression = z <= self.cut_expr["z"]
+                right_cut_expression = z >= self.cut_expr["z"]
+                key = "z"
             else:
                 raise KeyError()
 
@@ -227,38 +251,50 @@ class Node(object):
             left_mesh, self.left_template = self.mesh.cut(
                 left_cut_expression,
                 at_intersection=self.at_intersection,
-                return_template=True)
+                return_template=True,
+            )
             right_mesh, self.right_template = self.mesh.cut(
                 right_cut_expression,
                 at_intersection=self.at_intersection,
-                return_template=True)
+                return_template=True,
+            )
 
             # two new Nodes
-            self.left = Node(left_mesh,
-                             left_cuts,
-                             parent=self,
-                             internal_template=self.left_template, cut_expr=None,
-                             coord_sys=self.coord_sys,
-                             at_intersection=self.at_intersection,
-                             box=left_box)
-            self.right = Node(right_mesh,
-                              right_cuts,
-                              parent=self,
-                              internal_template=self.right_template, cut_expr=None,
-                              coord_sys=self.coord_sys,
-                              at_intersection=self.at_intersection,
-                              box=right_box)
+            self.left = Node(
+                left_mesh,
+                left_cuts,
+                parent=self,
+                internal_template=self.left_template,
+                cut_expr=None,
+                coord_sys=self.coord_sys,
+                at_intersection=self.at_intersection,
+                box=left_box,
+            )
+            self.right = Node(
+                right_mesh,
+                right_cuts,
+                parent=self,
+                internal_template=self.right_template,
+                cut_expr=None,
+                coord_sys=self.coord_sys,
+                at_intersection=self.at_intersection,
+                box=right_box,
+            )
 
     def _choose_next_cut(self):
+        """
+        Set self.cut_expr by choosing the dimension with the most remaining
+        cuts. Remove that cut from remaining cuts
+        """
         largest = 0
         for key in self.remaining_cuts:
             if len(self.remaining_cuts[key]) > largest:
                 largest = len(self.remaining_cuts[key])
                 largest_key = key
 
-        median = sorted(
-            self.remaining_cuts[largest_key]
-        )[int(0.5 * (len(self.remaining_cuts[largest_key]) - 1))]
+        median = sorted(self.remaining_cuts[largest_key])[
+            int(0.5 * (len(self.remaining_cuts[largest_key]) - 1))
+        ]
         # pop median cut from remaining cuts
         self.remaining_cuts[largest_key] = [
             x for x in self.remaining_cuts[largest_key] if x != median
@@ -323,22 +359,21 @@ class Node(object):
             if template.fields:
                 for idx in template.fields[0]:
                     template_field.append(self._convert_field_index(idx))
-                template.fields = [tfields.Tensors(template_field, dim=1,
-                                                   dtype=int)]
+                template.fields = [tfields.Tensors(template_field, dim=1, dtype=int)]
 
             template_map_field = []
             if len(template.maps[3]) > 0:
                 for idx in template.maps[3].fields[0]:
                     template_map_field.append(self._convert_map_index(idx))
-                template.maps[3].fields = [tfields.Tensors(template_map_field,
-                                                           dim=1, dtype=int)]
+                template.maps[3].fields = [
+                    tfields.Tensors(template_map_field, dim=1, dtype=int)
+                ]
             self._template = template
         return self._template
 
 
 class Searcher(Node):
-    def __init__(self, mesh, n_sections=None, delta=0.,
-                 cut_length=None):
+    def __init__(self, mesh, n_sections=None, delta=0.0, cut_length=None):
         """
         Special cutting tree root node.
         Provides a fast point in mesh search algorithm (Searcher.in_faces)
@@ -363,8 +398,9 @@ class Searcher(Node):
 
                 ab, ac = triangles.edges()
                 bc = ac - ab
-                side_lengths = np.concatenate([np.linalg.norm(side, axis=1)
-                                               for side in [ab, ac, bc]])
+                side_lengths = np.concatenate(
+                    [np.linalg.norm(side, axis=1) for side in [ab, ac, bc]]
+                )
                 # import mplTools as mpl
                 # axis= tfields.plotting.gca(2)
                 # mpl.plotHistogram(side_lengths, axis=axis)
@@ -384,18 +420,17 @@ class Searcher(Node):
         elif cut_length is not None:
             raise ValueError("cut_length not used.")
 
+        # build dictionary with cuts per dimension
         cut = {}
-        for i, key in enumerate(['x', 'y', 'z']):
+        for i, key in enumerate(["x", "y", "z"]):
             n_cuts = n_sections[i] + 1
-            if n_cuts <= 2:
-                values = []
-            else:
-                values = np.linspace(minima[i], maxima[i], n_cuts)[1:-1]
+            # [1:-1] because no need to cut at min or max
+            values = np.linspace(minima[i], maxima[i], n_cuts)[1:-1]
             cut[key] = values
 
-        return super(Searcher, self).__init__(mesh, cut,
-                                              at_intersection='keep',
-                                              delta=delta)
+        return super(Searcher, self).__init__(
+            mesh, cut, at_intersection="keep", delta=delta
+        )
 
     def in_faces(self, tensors, delta=-1, assign_multiple=False):
         """
@@ -403,25 +438,26 @@ class Searcher(Node):
             * check rare case of point+-delta outside box
 
         Examples:
-            # >>> import tfields
-            # >>> import numpy as np
-            # >>> mesh = tfields.Mesh3D.grid((0, 1, 2), (1, 2, 2), (2, 3, 2))
-            # >>> tree = tfields.bounding_box.Searcher(mesh, n_sections=[5, 5, 5])
-            # >>> points = tfields.Tensors([[0.5, 1, 2.1],
-            # ...                           [0.5, 0, 0],
-            # ...                           [0.5, 2, 2.1],
-            # ...                           [0.5, 1.5, 2.5]])
-            # >>> box_res = tree.in_faces(points, delta=0.0001)
-            # >>> usual_res = mesh.in_faces(points, delta=0.0001)
-            # >>> assert np.array_equal(box_res, usual_res)
+            >>> import tfields
+            >>> import numpy as np
+            >>> mesh = tfields.Mesh3D.grid((0, 1, 2), (1, 2, 2), (2, 3, 2))
+            >>> tree = tfields.bounding_box.Searcher(mesh)
+            >>> points = tfields.Tensors([[0.5, 1, 2.1],
+            ...                           [0.5, 0, 0],
+            ...                           [0.5, 2, 2.1],
+            ...                           [0.5, 1.5, 2.5]])
+            >>> box_res = tree.in_faces(points, delta=0.0001)
+            >>> usual_res = mesh.in_faces(points, delta=0.0001)
+            >>> assert np.array_equal(box_res, usual_res)
 
         """
-        raise ValueError("Broken feature. We are working on it!")
+        # raise ValueError("Broken feature. We are working on it!")
         if not self.is_root():
             raise ValueError("in_faces may only be called by root Node.")
-        if self.at_intersection != 'keep':
-            raise ValueError("Intersection method must be 'keep' for in_faces"
-                             "method.")
+        if self.at_intersection != "keep":
+            raise ValueError(
+                "Intersection method must be 'keep' for in_faces" "method."
+            )
 
         if self.mesh.nfaces() == 0:
             return np.empty((tensors.shape[0], 0), dtype=bool)
@@ -436,8 +472,7 @@ class Searcher(Node):
                     continue
                 if leaf.template.nfaces() == 0:
                     continue
-                leaf_mask = leaf.template.triangles()._in_triangles(point,
-                                                                    delta)
+                leaf_mask = leaf.template.triangles()._in_triangles(point, delta)
                 original_face_indices = leaf.template.maps[3].fields[0][leaf_mask]
                 if not assign_multiple and len(original_face_indices) > 0:
                     original_face_indices = original_face_indices[:1]
@@ -445,7 +480,8 @@ class Searcher(Node):
         return masks
 
 
-if __name__ == '__main__':
+if __name__ == "__main__":
     import doctest
+
     # doctest.run_docstring_examples(Searcher.in_faces, globals())
     doctest.testmod()