diff --git a/tfields/mesh3D.py b/tfields/mesh3D.py
index 1d25dfbe4c1fec042dba5bd4aa9e92688ea25271..f05bf60a7d419e3416f13b73743be3636fb2fb3b 100644
--- a/tfields/mesh3D.py
+++ b/tfields/mesh3D.py
@@ -555,11 +555,55 @@ class Mesh3D(tfields.TensorMaps):
         """
         Check whether points lie within triangles with Barycentric Technique
         see Triangles3D.in_triangles
+        If multiple requests are done on huge meshes,
+        this can be hugely optimized by requesting the property
+        self.tree or setting it to self.tree = <saved tree> before
+        calling in_faces
         """
-        masks = self.triangles().in_triangles(points, delta,
-                                              assign_multiple=assign_multiple)
+        key = 'mesh_tree'
+        if hasattr(self, '_cache') and key in self._cache:
+            log = logging.getLogger()
+            log.info(
+                "Using cached decision tree to speed up point - face mapping.")
+            masks = self.tree.in_faces(
+                points, delta, assign_multiple=assign_multiple)
+        else:
+            masks = self.triangles().in_triangles(
+                points, delta, assign_multiple=assign_multiple)
         return masks
 
+    @property
+    def tree(self):
+        """
+        Cached property to retrieve a bounding_box Searcher. This searcher can
+        hugely optimize 'in_faces' searches
+
+        Examples:
+            >>> mesh = tfields.Mesh3D.grid((0, 1, 3), (1, 2, 3), (2, 3, 3))
+            >>> _ = mesh.tree
+            >>> assert hasattr(mesh, '_cache')
+            >>> assert 'mesh_tree' in mesh._cache
+            >>> mask = mesh.in_faces(tfields.Points3D([[0.2, 1.2, 2.0]]),
+            ...                      0.00001)
+            >>> assert mask.sum() == 1  # one point in one triangle
+        """
+        if not hasattr(self, '_cache'):
+            self._cache = {}
+        key = 'mesh_tree'
+        if key in self._cache:
+            tree = self._cache[key]
+        else:
+            tree = tfields.bounding_box.Searcher(self)
+            self._cache[key] = tree
+        return tree
+
+    @tree.setter
+    def tree(self, tree):
+        if not hasattr(self, '_cache'):
+            self._cache = {}
+        key = 'mesh_tree'
+        self._cache[key] = tree
+
     def removeFaces(self, face_delete_mask):
         """
         Remove faces where face_delete_mask is True
@@ -621,6 +665,7 @@ class Mesh3D(tfields.TensorMaps):
 
         Examples:
             >>> import tfields
+            >>> import numpy as np
             >>> mp = tfields.TensorFields([[0,1,2],[2,3,0],[3,2,5],[5,4,3]],
             ...                           [1, 2, 3, 4])
             >>> m = tfields.Mesh3D([[0,0,0], [1,0,0], [1,1,0], [0,1,0], [0,2,0], [1,2,0]],
@@ -1110,4 +1155,5 @@ if __name__ == '__main__':  # pragma: no cover
     import doctest
 
     doctest.run_docstring_examples(Mesh3D.project, globals())
+    doctest.run_docstring_examples(Mesh3D.tree, globals())
     # doctest.testmod()