diff --git a/tfields/core.py b/tfields/core.py
index 2e6669768b76b0823d8cd0e4fdd59051e657e734..71b6dde0ef1230d174075643d7c1ca38ed6e8a10 100644
--- a/tfields/core.py
+++ b/tfields/core.py
@@ -1186,17 +1186,17 @@ class TensorFields(Tensors):
             >>> assert masked.fields[0].equal([42, 10.5])
             >>> assert masked.fields[1].equal([1, 3])
 
+            Iteration
+            >>> _ = [point for point in scalarField]
+
         """
         item = super(TensorFields, self).__getitem__(index)
         try:
             if issubclass(type(item), TensorFields):
                 if isinstance(index, tuple):
                     index = index[0]
-                if isinstance(index, slice):
+                if item.fields:
                     item.fields = [field.__getitem__(index) for field in item.fields]
-                else:
-                    if item.fields:
-                        item.fields = [field.__getitem__(index) for field in item.fields]
         except IndexError as err:
             warnings.warn("Index error occured for field.__getitem__. Error "
                           "message: {err}".format(**locals()))
@@ -1330,6 +1330,86 @@ class TensorMaps(TensorFields):
         obj.maps = maps
         return obj
 
+    def __getitem__(self, index):
+        """
+        In addition to the usual, also slice fields
+
+        Examples:
+            >>> import tfields
+            >>> import numpy as np
+            >>> vectors = tfields.Tensors([[0, 0, 0], [0, 0, 1], [0, -1, 0],
+            ...                            [1, 1, 1], [-1, -1, -1]])
+            >>> maps=[tfields.TensorFields([[0, 1, 2], [0, 1, 3], [2, 3, 4]],
+            ...                            [[1, 2], [3, 4], [5, 6]]),
+            ...       tfields.TensorFields([[0], [1], [2], [3], [4]])]
+            >>> mesh = tfields.TensorMaps(vectors,
+            ...                           [42, 21, 10.5, 1, 1],
+            ...                           [1, 2, 3, 3, 3],
+            ...                           maps=maps)
+
+            Slicing
+            >>> sliced = mesh[2:]
+            >>> assert isinstance(sliced, tfields.TensorMaps)
+            >>> assert isinstance(sliced.fields[0], tfields.Tensors)
+            >>> assert isinstance(sliced.maps[0], tfields.TensorFields)
+            >>> assert sliced.fields[0].equal([10.5, 1, 1])
+            >>> assert sliced.maps[0].equal([[0, 1, 2]])
+            >>> assert sliced.maps[0].fields[0].equal([[5, 6]])
+
+            Picking
+            >>> picked = mesh[1]
+            >>> assert np.array_equal(picked, [0, 0, 1])
+            >>> assert np.array_equal(picked.maps[0], np.empty((0, 3)))
+            >>> assert np.array_equal(picked.maps[1], [[0]])
+
+            Masking
+            >>> masked = mesh[[True, False, True, True, True]]
+            >>> assert masked.equal([[0, 0, 0], [0, -1, 0], [1, 1, 1], [-1, -1, -1]])
+            >>> assert masked.fields[0].equal([42, 10.5, 1, 1])
+            >>> assert masked.fields[1].equal([1, 3, 3, 3])
+            >>> assert masked.maps[0].equal([[1, 2, 3]])
+            >>> assert masked.maps[1].equal([[0], [1], [2], [3]])
+
+            Iteration
+            >>> _ = [vertex for vertex in mesh]
+
+        """
+        item = super(TensorMaps, self).__getitem__(index)
+        try:
+            if issubclass(type(item), TensorMaps):
+                item.maps = [mp.copy() for mp in item.maps]
+                if isinstance(index, tuple):
+                    index = index[0]
+                if item.maps:
+                    indices = np.array(range(len(self)))
+                    keep_indices = indices[index]
+                    if isinstance(keep_indices, int):
+                        keep_indices = [keep_indices]
+                    delete_indices = set(indices).difference(set(keep_indices))
+
+                    # correct all maps that contain deleted indices
+                    for mp_idx in range(len(self.maps)):
+                        # build mask, where the map should be deleted
+                        map_delete_mask = np.full((len(self.maps[mp_idx]),), False, dtype=bool)
+                        for i, mp in enumerate(self.maps[mp_idx]):
+                            for index in mp:
+                                if index in delete_indices:
+                                    map_delete_mask[i] = True
+                                    break
+                        map_mask = ~map_delete_mask
+
+                        # build the correction counters
+                        move_up_counter = np.zeros(self.maps[mp_idx].shape, dtype=int)
+                        for p in delete_indices:
+                            move_up_counter[self.maps[mp_idx] > p] -= 1
+
+                        item.maps[mp_idx] = (self.maps[mp_idx] + move_up_counter)[map_mask]
+        except IndexError as err:
+            warnings.warn("Index error occured for field.__getitem__. Error "
+                          "message: {err}".format(**locals()))
+
+        return item
+
     @classmethod
     def merged(cls, *objects, **kwargs):
         if not all([isinstance(o, cls) for o in objects]):
@@ -1345,12 +1425,12 @@ class TensorMaps(TensorFields):
             for map_field in o.maps:
                 map_field = map_field + cum_tensor_lengths[i]
                 try:
-                    map_index = dims.index(map_field.dim)
+                    mp_idx = dims.index(map_field.dim)
                 except ValueError:
                     maps.append(map_field)
                     dims.append(map_field.dim)
                 else:
-                    maps[map_index] = TensorFields.merged(maps[map_index], map_field)
+                    maps[mp_idx] = TensorFields.merged(maps[mp_idx], map_field)
         kwargs['maps'] = maps
 
         obj = super(TensorMaps, cls).merged(*objects, **kwargs) 
@@ -1424,16 +1504,14 @@ class TensorMaps(TensorFields):
             ...                          [3,3,3], [4,4,4], [5,6,7]],
             ...                         maps=[mp1, mp2])
 
-            # >>> c = tm.cleaned()
-            # >>> assert c.equal([[0., 0., 0.],
-            # ...                 [1., 1., 1.],
-            # ...                 [2., 2., 2.],
-            # ...                 [3., 3., 3.],
-            # ...                 [4., 4., 4.]])
-            # >>> c.maps
-            # [array([[0, 1, 2],
-            #         [0, 3, 4]]),
-            #  array([[0], [0]])]
+            >>> c = tm.cleaned()
+            >>> assert c.equal([[0., 0., 0.],
+            ...                 [1., 1., 1.],
+            ...                 [2., 2., 2.],
+            ...                 [3., 3., 3.],
+            ...                 [4., 4., 4.]])
+            >>> assert np.array_equal(c.maps[0], [[0, 1, 2], [0, 3, 4]])
+            >>> assert np.array_equal(c.maps[1], [[0], [0]])
 
 
         Returns:
@@ -1453,11 +1531,11 @@ class TensorMaps(TensorFields):
                 if duplicate_index != tensor_index:
                     stale_mask[tensor_index] = True
                     # redirect maps
-                    for map_index in range(len(self.maps)):
-                        for f in range(self.maps[map_index].shape[0]):
-                            if tensor_index in self.maps[f]:
-                                index = tfields.index(self.maps[map_index][f], tensor_index)
-                                inst.maps[map_index][f][index] = duplicate_index
+                    for mp_idx in range(len(self.maps)):
+                        for f in range(len(self.maps[mp_idx])):
+                            if tensor_index in self.maps[mp_idx][f]:
+                                index = tfields.index(self.maps[mp_idx][f], tensor_index)
+                                inst.maps[mp_idx][f][index] = duplicate_index
 
         return inst.removed(stale_mask)
 
@@ -1475,7 +1553,7 @@ class TensorMaps(TensorFields):
             ...                                                    [3, 4, 6]],
             ...                                                   [1, 3, 5, 7, 9],
             ...                                                   [2, 4, 6, 8, 0])])
-            >>> c = m.removed([True, True, True, False, False, False, False])
+            >>> c = m.keep([False, False, False, True, True, True, True])
             >>> c.equal([[0, 0, 0],
             ...          [3, 3, 3],
             ...          [4, 4, 4],
@@ -1489,20 +1567,51 @@ class TensorMaps(TensorFields):
 
         """
         remove_condition = np.array(remove_condition)
-        # built instance that only contains the vaild points
-        inst = self[~remove_condition].copy()
-        delete_indices = np.arange(self.shape[0])[remove_condition]
-        face_keep_masks = self.to_maps_masks(~remove_condition)
+        # # built instance that only contains the vaild points
+        # inst = self[~remove_condition].copy()
+        # delete_indices = np.arange(self.shape[0])[remove_condition]
+        # face_keep_masks = self.to_maps_masks(~remove_condition)
 
-        for map_index, face_keep_mask in enumerate(face_keep_masks):
-            moveUpCounter = np.zeros(self.maps[map_index].shape, dtype=int)
+        # for mp_idx, face_keep_mask in enumerate(face_keep_masks):
+        #     move_up_counter = np.zeros(self.maps[mp_idx].shape, dtype=int)
 
-            # correct map:
-            for p in delete_indices:
-                moveUpCounter[self.maps[map_index] > p] -= 1
+        #     # correct map:
+        #     for p in delete_indices:
+        #         move_up_counter[self.maps[mp_idx] > p] -= 1
 
-            inst.maps[map_index] = (self.maps[map_index] + moveUpCounter)[face_keep_mask]
-        return inst
+        #     inst.maps[mp_idx] = (self.maps[mp_idx] + move_up_counter)[face_keep_mask]
+        # return inst
+        return self[~remove_condition]
+
+    def keep(self, keep_condition):
+        """
+        Return copy of self with vertices where keep_condition is True
+        Copy because self is immutable
+
+        Examples:
+            >>> import tfields
+            >>> m = tfields.TensorMaps([[0,0,0], [1,1,1], [2,2,2], [0,0,0],
+            ...                         [3,3,3], [4,4,4], [5,5,5]],
+            ...                        maps=[tfields.TensorFields([[0, 1, 2], [0, 1, 3],
+            ...                                                    [3, 4, 5], [3, 4, 1],
+            ...                                                    [3, 4, 6]],
+            ...                                                   [1, 3, 5, 7, 9],
+            ...                                                   [2, 4, 6, 8, 0])])
+            >>> c = m.removed([True, True, True, False, False, False, False])
+            >>> c.equal([[0, 0, 0],
+            ...          [3, 3, 3],
+            ...          [4, 4, 4],
+            ...          [5, 5, 5]])
+            True
+            >>> c.maps[0]
+            TensorFields([[0, 1, 2],
+                          [0, 1, 3]])
+            >>> assert c.maps[0].fields[0].equal([5, 9])
+            >>> assert c.maps[0].fields[1].equal([6, 0])
+
+        """
+        keep_condition = np.array(keep_condition)
+        return self[keep_condition]
 
     def to_maps_masks(self, mask):
         """
@@ -1510,24 +1619,30 @@ class TensorMaps(TensorFields):
             >>> from tfields import TensorMaps
             >>> import numpy as np
             >>> m = TensorMaps([[1,2,3], [3,3,3], [0,0,0], [5,6,7]],
-            ...                maps=[[[0, 1, 2], [1, 2, 3]]])
+            ...                maps=[[[0, 1, 2], [1, 2, 3]],
+            ...                      [[0], [3]]])
             >>> from sympy.abc import x,y,z
             >>> vertexMask = m.evalf(z < 6)
             >>> faceMask = m.to_maps_masks(vertexMask)
-            >>> assert np.array_equal(faceMask, [[True, False]])
+            >>> assert np.array_equal(faceMask, [[True, False], [True, False]])
+            >>> index_face_mask = m.to_maps_masks(0)
+            >>> assert np.array_equal(index_face_mask, [[False, False], [True, False]])
 
         Returns:
-            mask of faces with all vertices in mask
+            masks of maps with all vertices in mask
         """
-
-        indices = np.array(range(len(mask)))
-        delete_indices = indices[~mask]
-        delete_indices = set(delete_indices)  # set speeds up everything enormously
+        indices = np.array(range(len(self)))
+        keep_indices = indices[mask]
+        if isinstance(keep_indices, int):
+            keep_indices = [keep_indices]
+        delete_indices = set(indices).difference(set(keep_indices))
+        # delete_indices = indices[~mask]
+        # delete_indices = set(delete_indices)  # set speeds up everything enormously
 
         masks = []
-        for map_index in range(len(self.maps)):
-            map_delete_mask = np.full((len(self.maps[map_index]),), False, dtype=bool)
-            for i, mp in enumerate(self.maps[map_index]):
+        for mp_idx in range(len(self.maps)):
+            map_delete_mask = np.full((len(self.maps[mp_idx]),), False, dtype=bool)
+            for i, mp in enumerate(self.maps[mp_idx]):
                 for index in mp:
                     if index in delete_indices:
                         map_delete_mask[i] = True
@@ -1539,6 +1654,10 @@ class TensorMaps(TensorFields):
 if __name__ == '__main__':  # pragma: no cover
     import doctest
     doctest.testmod()
+    # doctest.run_docstring_examples(TensorFields.__getitem__, globals())
+    # doctest.run_docstring_examples(TensorMaps.__getitem__, globals())
+    # doctest.run_docstring_examples(TensorMaps.keep, globals())
+    # doctest.run_docstring_examples(TensorMaps.to_maps_masks, globals())
     # doctest.run_docstring_examples(TensorMaps, globals())
     # doctest.run_docstring_examples(TensorFields, globals())
     # doctest.run_docstring_examples(AbstractNdarray.copy, globals())
diff --git a/tfields/mesh3D.py b/tfields/mesh3D.py
index 7ccb8442d38ef21b81b6fe36652545f443057d71..9b9ea64a0a0b7caffb1e2584e7da13b2c52d7d1b 100644
--- a/tfields/mesh3D.py
+++ b/tfields/mesh3D.py
@@ -119,7 +119,7 @@ class Mesh3D(tfields.TensorMaps):
             raise ValueError("Conflicting options maps and faces")
         if maps is not None:
             kwargs['maps'] = maps
-        if faceScalars:
+        if len(faceScalars) > 0:
             map_fields = scalars_to_fields(faceScalars)
         else:
             map_fields = []
@@ -133,6 +133,61 @@ class Mesh3D(tfields.TensorMaps):
             raise ValueError("Face dimension should be 3")
         return obj
 
+    @classmethod
+    def plane(cls, *base_vectors, **kwargs):
+        vertices = tfields.Tensors.grid(*base_vectors)
+        fix_coord = None
+        for coord in range(3):
+            if len(base_vectors[coord]) > 1:
+                continue
+            if len(base_vectors[coord]) == 0:
+                continue
+            fix_coord = coord
+
+        var_coords = list(range(3))
+        var_coords.pop(var_coords.index(fix_coord))
+
+        faces = []
+        base0, base1 = base_vectors[var_coords[0]], base_vectors[var_coords[1]]
+        for i1 in range(len(base1) - 1):
+            for i0 in range(len(base0) - 1):
+                idx_top_left = len(base1) * (i0 + 0) + (i1 + 0)
+                idx_top_right = len(base1) * (i0 + 0) + (i1 + 1)
+                idx_bot_left = len(base1) * (i0 + 1) + (i1 + 0)
+                idx_bot_right = len(base1) * (i0 + 1) + (i1 + 1)
+                faces.append([idx_top_left, idx_top_right, idx_bot_left])
+                faces.append([idx_top_right, idx_bot_left, idx_bot_right])
+        inst = cls.__new__(cls, vertices, faces=faces, **kwargs)
+        return inst
+
+    @classmethod
+    def grid(cls, *base_vectors, **kwargs):
+        if not len(base_vectors) == 3:
+            raise AttributeError("3 base_vectors vectors required")
+
+        indices = [0, -1]
+        coords = range(3)
+        baseLengthsAbove1 = [len(b) > 1 for b in base_vectors]
+        # if one plane is given: rearrange indices and coords
+        if not all(baseLengthsAbove1):
+            indices = [0]
+            for i, b in enumerate(baseLengthsAbove1):
+                if not b:
+                    coords = [i]
+                    break
+
+        base_vectors = list(base_vectors)
+        planes = []
+        for ind in indices:
+            for coord in coords:
+                basePart = base_vectors[:]
+                basePart[coord] = np.array([base_vectors[coord][ind]],
+                                           dtype=float)
+
+                planes.append(cls.plane(*basePart))
+        inst = cls.merged(*planes, **kwargs)
+        return inst
+
     @property
     def faces(self):
         return maps_to_faces(self.maps)
@@ -330,108 +385,6 @@ class Mesh3D(tfields.TensorMaps):
 
         return ~faceDeleteMask
 
-    def getRemovedVertices(self, vertexDeleteMask):
-        """
-        Return copy of self without vertices where vertexDeleteMask is True
-        Copy because self is immutable
-
-        Examples:
-            >>> m = tfields.Mesh3D([[0,0,0], [1,1,1], [2,2,2], [0,0,0],
-            ...                     [3,3,3], [4,4,4], [5,5,5]],
-            ...                    [[0, 1, 2], [0, 1, 3], [3, 4, 5], [3, 4, 1],
-            ...                     [3, 4, 6]],
-            ...                    faceScalars=[[1,2], [3,4], [5,6], [7,8], [9,0]])
-            >>> c = m.getRemovedVertices([True, True, True, False, False,
-            ...                           False, False])
-            >>> c
-            Mesh3D([[ 0.,  0.,  0.],
-                    [ 3.,  3.,  3.],
-                    [ 4.,  4.,  4.],
-                    [ 5.,  5.,  5.]])
-            >>> c.faces
-            array([[0, 1, 2],
-                   [0, 1, 3]])
-            >>> c.faceScalars
-            array([[ 5.,  6.],
-                   [ 9.,  0.]])
-        
-        """
-        log = logger.new()
-        vertexDeleteMask = np.array(vertexDeleteMask)
-        log.verbose("Remove {0} vertices.".format(vertexDeleteMask.sum()))
-        # built instance that only contains the vaild points
-        inst = self[~vertexDeleteMask].copy()
-
-        moveUpCounter = np.zeros(self.faces.shape, dtype=int)
-
-        # correct faces:
-        deleteIndices = np.arange(len(self))[vertexDeleteMask]
-        for p in deleteIndices:
-            moveUpCounter[self.faces > p] -= 1
-
-        faceKeepMask = self.getFaceMask(~vertexDeleteMask)
-        inst.faces = (self.faces + moveUpCounter)[faceKeepMask]
-        inst.faceScalars = self.faceScalars[faceKeepMask]
-        return inst
-
-    def cleaned(self, stale=True, duplicates=True):
-        """
-        Args:
-            stale (bool): remove stale vertices
-            duplicates (bool): replace duplicate vertices by originals
-        Examples:
-            >>> m = tfields.Mesh3D([[0,0,0], [1,1,1], [2,2,2], [0,0,0],
-            ...                     [3,3,3], [4,4,4], [5,6,7]],
-            ...                    [[0, 1, 2], [3, 4, 5]],
-            ...                    faceScalars=[[1,2,3,4,5], [6,7,8,9,0]])
-            >>> c = m.clean()
-            >>> c
-            Mesh3D([[ 0.,  0.,  0.],
-                    [ 1.,  1.,  1.],
-                    [ 2.,  2.,  2.],
-                    [ 3.,  3.,  3.],
-                    [ 4.,  4.,  4.]])
-            >>> c.faces
-            array([[0, 1, 2],
-                   [0, 3, 4]])
-
-        Returns:
-            copy of self without stale vertices and duplicat points
-        """
-        log = logger.new()
-        log.verbose("Cleaning up.")
-        # remove stale vertices
-        if stale:
-            vertexDeleteMask = self.staleVertices()
-        else:
-            vertexDeleteMask = np.full(self.shape[0], False, dtype=bool)
-        # remove duplicates in order to not have any artificial separations
-        inst = self
-        if duplicates:
-            inst = self.copy()
-            log.verbose("Finding Duplicates")
-            dups = tfields.duplicates(self, axis=0)
-            for i, dupi in zip(range(self.shape[0]), dups):
-                if dupi != i:
-                    log.verbose("Found Duplicate point @ index {0}".format(i))
-                    vertexDeleteMask[i] = True
-                    # redirect faces
-                    log.verbose("Run trough all faces to let it point to the"
-                                "original")
-                    for f in range(self.nfaces()):
-                        if i in self.faces[f]:
-                            index = tfields.index(self.faces[f], i)
-                            inst.faces[f][index] = dupi
-
-        return inst.getRemovedVertices(vertexDeleteMask)
-
-    def clean(self, *args, **kwargs):
-        """
-        Deprecated
-        """
-        warnings.warn("Name clean is deprecated. take cleaned instead", DeprecationWarning)
-        return self.cleaned(*args, **kwargs)
-
     def getScalarMap(self, mask):
         """
         Return copy of self without vertices where mask is True
@@ -568,7 +521,7 @@ class Mesh3D(tfields.TensorMaps):
                               else False
                               for i in faceIndices]
             mesh.removeFaces(faceDeleteMask)
-            mesh = mesh.getRemovedVertices(mesh.staleVertices())
+            mesh = mesh.cleaned(duplicates=False)
             parts.append(mesh)
         return parts
 
@@ -628,7 +581,7 @@ class Mesh3D(tfields.TensorMaps):
             >>> base = [np.linspace(0, 1, 10),
             ...         np.linspace(0, 1, 10),
             ...         np.linspace(0, 1, 10)]
-            >>> mesh = tfields.Mesh3D.createMeshGrid(*base).cleaned()
+            >>> mesh = tfields.Mesh3D.grid(*base).cleaned()
 
             Some small mistake occured in the test. Check that.
             # Select the first face as a seedFace
@@ -688,7 +641,7 @@ class Mesh3D(tfields.TensorMaps):
             >>> base = [np.linspace(0, 1, 2),
             ...         np.linspace(0, 1, 4),
             ...         np.linspace(0, 1, 4)]
-            >>> mesh = tfields.Mesh3D.createMeshGrid(*base).cleaned()
+            >>> mesh = tfields.Mesh3D.grid(*base).cleaned()
 
             Select the first face as a seedFace
             >>> faceGroups = mesh.getSides([[1,0,0],[0,1,0],[0,0,1]])
@@ -903,9 +856,9 @@ class Mesh3D(tfields.TensorMaps):
             return self.copy(), self.copy()
 
         # mask for points that do not fulfill the cut expression
-        mask = self.evalf(expression, coordSys=coordSys)
+        mask = self.evalf(expression)
         # remove the points
-        inst = self.getRemovedVertices(~mask)
+        inst = self[mask].copy()
         scalarMap = self.getScalarMap(mask)
         scalarMap = scalarMap.reshape((-1, 1)).astype(float)
 
@@ -1057,9 +1010,10 @@ class Mesh3D(tfields.TensorMaps):
             >>> from sympy.abc import x,y,z
             >>> cutExpr = x > 1.5
 
-            >>> m = tfields.Mesh3D.createMeshGrid(np.linspace(0,3,4),
-            ...                           np.linspace(0,3,4),
-            ...                           np.linspace(0, 0, 1), faceScalars=np.linspace(4, 8, 18))
+            >>> m = tfields.Mesh3D.grid((0,3,4),
+            ...                         (0,3,4),
+            ...                         (0, 0, 1))
+            >>> m.faceScalars = np.linspace(4, 8, 18)
             >>> mNew = m.cut(cutExpr)
             >>> len(mNew)
             8
@@ -1085,7 +1039,7 @@ class Mesh3D(tfields.TensorMaps):
             copy of cut mesh
 
         """
-        with self.tmp_transform(coordSys):
+        with self.tmp_transform(coordSys or self.coordSys):
             if isinstance(expression, Mesh3D):
                 obj = self._cut_template(expression)
             else: