diff --git a/tfields/__about__.py b/tfields/__about__.py
index 31f8969f0ff3c5802cbd5d7d4d3faf0a4b9ec79b..564163345ccf3a8a84b823034a90dd16f3350da0 100644
--- a/tfields/__about__.py
+++ b/tfields/__about__.py
@@ -9,7 +9,7 @@ __all__ = [
     "__license__",
     "__copyright__",
 ]
-__version__ = '0.1.0.dev9'
+__version__ = '0.1.0.dev10'
 __title__ = 'tfields'
 __summary__ = "numpy + sympy implementation of tensor fields with attached coordinate systems"
 __keywords__ = "tensors coordinate system trafo sympy numpy"
diff --git a/tfields/core.py b/tfields/core.py
index 6ab30325a286dec75a91ff96e14d92f63ff44a64..14d7daaace3956bdeaf27ec08e9cf709911a7097 100644
--- a/tfields/core.py
+++ b/tfields/core.py
@@ -126,12 +126,12 @@ class AbstractNdarray(np.ndarray):
             >>> from tfields import Tensors, TensorFields
             >>> scalars = Tensors([0, 1, 2])
             >>> vectors = Tensors([[0, 0, 0], [0, 0, 1], [0, -1, 0]])
-            >>> scalarField = TensorFields(vectors, scalars, coord_sys='cylinder')
+            >>> scalar_field = TensorFields(vectors, scalars, coord_sys='cylinder')
 
             Save it and restore it
             >>> out_file = NamedTemporaryFile(suffix='.pickle')
 
-            >>> pickle.dump(scalarField,
+            >>> pickle.dump(scalar_field,
             ...             out_file)
             >>> _ = out_file.seek(0)
 
@@ -301,7 +301,7 @@ class AbstractNdarray(np.ndarray):
         Recursively walk trough all __slots__ and describe all elements
         """
         d = {}
-        d['bulk'] = np.array(self)
+        d['bulk'] = self.bulk
         d['bulk_type'] = self.__class__.__name__
         for attr in self._iter_slots():
             value = getattr(self, attr)
@@ -517,6 +517,21 @@ class Tensors(AbstractNdarray):
 
         return obj
 
+    def __iter__(self):
+        """
+        Forwarding iterations to the bulk array. Otherwise __getitem__ would
+        kick in and slow down imensely.
+        Examples:
+            >>> import tfields
+            >>> vectors = tfields.Tensors([[0, 0, 0], [0, 0, 1], [0, -1, 0]])
+            >>> scalar_field = tfields.TensorFields(vectors, [42, 21, 10.5], [1, 2, 3])
+            >>> [(point.rank, point.dim) for point in scalar_field]
+            [(0, 1), (0, 1), (0, 1)]
+
+        """
+        for index in range(len(self)):
+            yield super(Tensors, self).__getitem__(index).view(Tensors)
+
     @classmethod
     def merged(cls, *objects, **kwargs):
         """
@@ -664,6 +679,14 @@ class Tensors(AbstractNdarray):
         inst = cls.__new__(cls, tfields.lib.grid.igrid(*base_vectors, **kwargs))
         return inst
 
+    @property
+    def bulk(self):
+        """
+        The pure ndarray version of the actual state
+            -> nothing attached
+        """
+        return np.array(self)
+
     @property
     def rank(self):
         """
@@ -1203,10 +1226,10 @@ class TensorFields(Tensors):
         >>> from tfields import Tensors, TensorFields
         >>> scalars = Tensors([0, 1, 2])
         >>> vectors = Tensors([[0, 0, 0], [0, 0, 1], [0, -1, 0]])
-        >>> scalarField = TensorFields(vectors, scalars)
-        >>> scalarField.rank
+        >>> scalar_field = TensorFields(vectors, scalars)
+        >>> scalar_field.rank
         1
-        >>> scalarField.fields[0].rank
+        >>> scalar_field.fields[0].rank
         0
         >>> vectorField = TensorFields(vectors, vectors)
         >>> vectorField.fields[0].rank
@@ -1288,26 +1311,26 @@ class TensorFields(Tensors):
             >>> import tfields
             >>> import numpy as np
             >>> vectors = tfields.Tensors([[0, 0, 0], [0, 0, 1], [0, -1, 0]])
-            >>> scalarField = tfields.TensorFields(vectors, [42, 21, 10.5], [1, 2, 3])
+            >>> scalar_field = tfields.TensorFields(vectors, [42, 21, 10.5], [1, 2, 3])
 
             Slicing
-            >>> sliced = scalarField[2:]
+            >>> sliced = scalar_field[2:]
             >>> assert isinstance(sliced, tfields.TensorFields)
             >>> assert isinstance(sliced.fields[0], tfields.Tensors)
             >>> assert sliced.fields[0].equal([10.5])
 
             Picking
-            >>> picked = scalarField[1]
+            >>> picked = scalar_field[1]
             >>> assert np.array_equal(picked, [0, 0, 1])
 
             Masking
-            >>> masked = scalarField[[True, False, True]]
+            >>> masked = scalar_field[[True, False, True]]
             >>> assert masked.equal([[0, 0, 0], [0, -1, 0]])
             >>> assert masked.fields[0].equal([42, 10.5])
             >>> assert masked.fields[1].equal([1, 3])
 
             Iteration
-            >>> _ = [point for point in scalarField]
+            >>> _ = [point for point in scalar_field]
 
         """
         item = super(TensorFields, self).__getitem__(index)
diff --git a/tfields/mesh3D.py b/tfields/mesh3D.py
index c3e7c0dc6ecc8023db79d4e2378460b2a95bb51f..0237f03d46362703287fd7b837089d6fe837f107 100644
--- a/tfields/mesh3D.py
+++ b/tfields/mesh3D.py
@@ -590,6 +590,8 @@ class Mesh3D(tfields.TensorMaps):
             ...                       [[-42, -21], [42, 21]])
                                    
         """
+        # Possible Extension (small todo): check: len(field(s)) == len(self/maps)
+
         # Redirect fields
         fields = []
         if template.fields:
@@ -610,12 +612,14 @@ class Mesh3D(tfields.TensorMaps):
                     projected_field[nan_mask] = np.nan  # correction for nan
                     fields.append(projected_field)
 
-        # Redirect maps fields
+        # Redirect maps and their fields
         maps = []
         for mp, template_mp in zip(self.maps, template.maps):
             mp_fields = []
-            if len(template_mp.fields) > 0:
-                for field in mp.fields:
+            for field in mp.fields:
+                if len(template_mp) == 0 and len(template_mp.fields) == 0:
+                    mp_fields.append(field[0:0])  # np.empty
+                else:
                     mp_fields.append(field[template_mp.fields[0].astype(int)])
             new_mp = tfields.TensorFields(tfields.Tensors(template_mp),
                                           *mp_fields)
diff --git a/tfields/plotting/mpl.py b/tfields/plotting/mpl.py
index 0c40c822cc752ccc05770784dfedbfbb02014ab6..8588685b3020ca964fc102c101c6718ca23606ba 100644
--- a/tfields/plotting/mpl.py
+++ b/tfields/plotting/mpl.py
@@ -194,8 +194,6 @@ def plot_mesh(vertices, faces, **kwargs):
         vmin
         vmax
     """
-    vertices = np.array(vertices)
-    faces = np.array(faces)
     if faces.shape[0] == 0:
         warnings.warn("No faces to plot")
         return None
diff --git a/tfields/triangles3D.py b/tfields/triangles3D.py
index 54fbba700ecd2651331a26e721312859a3d2bf51..d2a67afb2056c001870fbcf657b65fab54c89a3e 100644
--- a/tfields/triangles3D.py
+++ b/tfields/triangles3D.py
@@ -215,18 +215,21 @@ class Triangles3D(tfields.TensorFields):
 
             if not np.array_equal(transform, np.eye(3)):
                 a = np.linalg.norm(np.linalg.solve(transform.T,
-                                                   (self[aIndices, :] - self[bIndices, :]).T),
+                                                   (self.bulk[aIndices, :] -
+                                                    self.bulk[bIndices, :]).T),
                                    axis=0)
                 b = np.linalg.norm(np.linalg.solve(transform.T,
-                                                   (self[aIndices, :] - self[cIndices, :]).T),
+                                                   (self.bulk[aIndices, :] -
+                                                    self.bulk[cIndices, :]).T),
                                    axis=0)
                 c = np.linalg.norm(np.linalg.solve(transform.T,
-                                                   (self[bIndices, :] - self[cIndices, :]).T),
+                                                   (self.bulk[bIndices, :] -
+                                                    self.bulk[cIndices, :]).T),
                                    axis=0)
             else:
-                a = np.linalg.norm(self[aIndices, :] - self[bIndices, :], axis=1)
-                b = np.linalg.norm(self[aIndices, :] - self[cIndices, :], axis=1)
-                c = np.linalg.norm(self[bIndices, :] - self[cIndices, :], axis=1)
+                a = np.linalg.norm(self.bulk[aIndices, :] - self.bulk[bIndices, :], axis=1)
+                b = np.linalg.norm(self.bulk[aIndices, :] - self.bulk[cIndices, :], axis=1)
+                c = np.linalg.norm(self.bulk[bIndices, :] - self.bulk[cIndices, :], axis=1)
 
             # sort by length for numerical stability
             lengths = np.concatenate((a.reshape(-1, 1), b.reshape(-1, 1), c.reshape(-1, 1)), axis=1)
@@ -245,9 +248,9 @@ class Triangles3D(tfields.TensorFields):
         bIndices = [i * 3 + 1 for i in indices]
         cIndices = [i * 3 + 2 for i in indices]
 
-        a = self[aIndices, :]
-        b = self[bIndices, :]
-        c = self[cIndices, :]
+        a = self.bulk[aIndices, :]
+        b = self.bulk[bIndices, :]
+        c = self.bulk[cIndices, :]
         return a, b, c
 
     def circumcenters(self):
@@ -500,6 +503,7 @@ class Triangles3D(tfields.TensorFields):
             ValueError: point must be castable to shape 3 but is of shape (2, 3)
 
         """
+
         if self.ntriangles() == 0:
             return np.array([], dtype=bool)
 
@@ -537,9 +541,9 @@ class Triangles3D(tfields.TensorFields):
 
         Returns:
             np.ndarray: 2-d mask which is True, where a point is in a triangle
-                              triangle index ->
+                              triangle index (axis = 1)->
                 points index
-                                    1        0       0
+                (axis = 0)          1        0       0
                       |             1        0       0
                       V             0        0       1
 
@@ -548,7 +552,7 @@ class Triangles3D(tfields.TensorFields):
                 triangle indices can have multiple true values
 
                 For Example, if you want to know the number of points in one
-                face, just do:
+                face, just sum over all points:
                 >> tris.in_triangles(poits).sum(axis=0)[face_index]
 
         """
@@ -557,9 +561,10 @@ class Triangles3D(tfields.TensorFields):
 
         masks = np.zeros((len(tensors), self.ntriangles()), dtype=bool)
         with tensors.tmp_transform(self.coord_sys):
-            for i, point in enumerate(tensors):
+            for i, point in enumerate(iter(tensors)):
                 # print("Status: {i} points processed. Yet {nP} are found "
-                #       "to be in triangles.".format(nP=masks.sum() - 1, **locals()))
+                #       "to be in triangles (before handling assign_multiple)."
+                #       .format(nP=masks.sum(), **locals()))  # SLOW!
                 masks[i] = self._in_triangles(point, delta)
         # print("Status: All points processed. Yet {nP} are found "
         #       "to be in triangles.".format(nP=masks.sum() - 1, **locals()))