diff --git a/examples/03_variogram/README.rst b/examples/03_variogram/README.rst
index d16c16c2b..b0339f12d 100644
--- a/examples/03_variogram/README.rst
+++ b/examples/03_variogram/README.rst
@@ -1,3 +1,5 @@
+.. _tutorial_03_variogram:
+
 Tutorial 3: Variogram Estimation
 ================================
 
diff --git a/examples/08_mesh/00_mesh_intro.py b/examples/08_mesh/00_mesh_intro.py
new file mode 100644
index 000000000..b45bb123d
--- /dev/null
+++ b/examples/08_mesh/00_mesh_intro.py
@@ -0,0 +1,75 @@
+"""
+Example: Using a mesh for GSTools
+---------------------------------
+
+This example shows how external data can be analysed with GSTools.
+
+
+Pretending we have some Data
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+We will pretend that we have some external data, by creating some random data
+with GSTools, but we will delete the objects and only use the data, without the
+backend GSTools provides.
+"""
+from datetime import date
+import numpy as np
+import gstools as gs
+
+# create a circular mesh
+point_no = 10000
+rng = np.random.RandomState(20170521)
+r = 50.0 * np.sqrt(rng.uniform(size=point_no))
+theta = 2.0 * np.pi * rng.uniform(size=point_no)
+x = r * np.cos(theta)
+y = r * np.sin(theta)
+
+tmp_model = gs.Exponential(dim=2, var=1.5, len_scale=10.0)
+tmp_srf = gs.SRF(tmp_model)
+field = tmp_srf((x, y))
+tmp_srf.plot()
+
+# Now that we have our data, let's delete everything GSTools related and pretend
+# that this has never happend
+del(tmp_model)
+del(tmp_srf)
+
+
+# Creating the Mesh
+# ^^^^^^^^^^^^^^^^^
+#
+# Starting out fresh, we want to feed the mesh with our data
+mesh = gs.Mesh(2, pos=(x, y), values=field)
+
+# We can add meta data too
+mesh.set_field_data("Süderbrarup", "location")
+mesh.set_field_data(date(year=2020, month=2, day=28), "date")
+
+# This can be conviniently accessed
+print(mesh.location)
+print(mesh.date)
+
+# But the meta data is also collected as a dictionary in case you want to export
+# it
+print(mesh.field_data)
+
+
+# Estimating the Variogram
+# ^^^^^^^^^^^^^^^^^^^^^^^^
+# Now, with our mesh, which was loaded from completely external sources, we can
+# estimate the variogram of the data.
+# To speed things up, we will only use a fraction of the available data
+
+bins = np.linspace(0, 50, 50)
+bin_centre, gamma = gs.vario_estimate_unstructured(
+    (x, y), field, bins, sampling_size=2000, sampling_seed=19900408)
+
+# As we are experts, we'll do an expert guess and say, that we will most likely
+# have data that has an exponential variogram. Non-experts can have a look at
+# the "Finding the best fitting variogram model" tutorial in
+# :ref:`tutorial_03_variogram`.
+fit_model = gs.Exponential(dim=2)
+fit_model.fit_variogram(bin_centre, gamma, nugget=False)
+
+ax = fit_model.plot(x_max=max(bin_centre))
+ax.plot(bin_centre, gamma)
diff --git a/examples/08_mesh/README.rst b/examples/08_mesh/README.rst
new file mode 100644
index 000000000..2b7196277
--- /dev/null
+++ b/examples/08_mesh/README.rst
@@ -0,0 +1,12 @@
+.. _tutorial_08_mesh:
+
+Tutorial 8: Working With Meshes
+===============================
+
+In case you have external data, which does not come from the spatial random
+generation or similar methods of GSTools and you want to analyse this data, e.g.
+estimate the variogram, you can load the data into a mesh class and conveniently
+manipulate, analyse, visualise, and export the data with GSTools.
+
+Gallery
+-------
diff --git a/gstools/__init__.py b/gstools/__init__.py
index 1bafdfe4e..6d94acc6d 100644
--- a/gstools/__init__.py
+++ b/gstools/__init__.py
@@ -33,6 +33,8 @@
 
 .. autosummary::
    SRF
+   Field
+   Mesh
 
 Covariance Base-Class
 ^^^^^^^^^^^^^^^^^^^^^
@@ -100,6 +102,7 @@
 """
 
 from gstools import field, variogram, random, covmodel, tools, krige, transform
+from gstools.field.mesh import Mesh
 from gstools.field import SRF
 from gstools.tools import (
     vtk_export,
diff --git a/gstools/field/__init__.py b/gstools/field/__init__.py
index 6884e3775..a90322a4f 100644
--- a/gstools/field/__init__.py
+++ b/gstools/field/__init__.py
@@ -17,10 +17,14 @@
 
 .. autosummary::
    SRF
+   Field
+   Mesh
 
 ----
 """
 
+from gstools.field.mesh import Mesh
+from gstools.field.base import Field
 from gstools.field.srf import SRF
 
-__all__ = ["SRF"]
+__all__ = ["SRF", "Field", "Mesh"]
diff --git a/gstools/field/base.py b/gstools/field/base.py
index 75ae04679..47dd9db21 100755
--- a/gstools/field/base.py
+++ b/gstools/field/base.py
@@ -15,8 +15,8 @@
 
 import numpy as np
 
+from gstools.field import Mesh
 from gstools.covmodel.base import CovModel
-from gstools.tools.export import to_vtk, vtk_export
 from gstools.field.tools import (
     _get_select,
     check_mesh,
@@ -29,30 +29,38 @@
 __all__ = ["Field"]
 
 
-class Field:
-    """A field base class for random and kriging fields ect.
+class Field(Mesh):
+    """A field base class for random and kriging fields, etc.
 
     Parameters
     ----------
     model : :any:`CovModel`
         Covariance Model related to the field.
-    mean : :class:`float`, optional
-        Mean value of the field.
     """
 
-    def __init__(self, model, mean=0.0):
+    def __init__(
+        self,
+        model,
+        points=None,
+        name="field",
+        values=None,
+        mesh_type="unstructured",
+        mean=0.0,
+    ):
         # initialize attributes
-        self.pos = None
-        self.mesh_type = None
-        self.field = None
+        super().__init__(
+            dim=model.dim,
+            points=points,
+            name=name,
+            values=values,
+            mesh_type=mesh_type,
+        )
         # initialize private attributes
-        self._mean = None
         self._model = None
-        self.mean = mean
         self.model = model
-        self._value_type = None
+        self.mean = mean
 
-    def __call__(*args, **kwargs):
+    def __call__(self, *args, **kwargs):
         """Generate the field."""
         pass
 
@@ -74,7 +82,7 @@ def unstructured(self, *args, **kwargs):
 
     def mesh(
         self, mesh, points="centroids", direction="xyz", name="field", **kwargs
-    ):  # pragma: no cover
+    ):
         """Generate a field on a given meshio or ogs5py mesh.
 
         Parameters
@@ -119,7 +127,7 @@ def mesh(
                 pnts = mesh.centroids_flat.T[select]
             else:
                 pnts = mesh.NODES.T[select]
-            out = self.unstructured(pos=pnts, **kwargs)
+            out = self.unstructured(points=pnts, **kwargs)
         else:
             if points == "centroids":
                 # define unique order of cells
@@ -132,9 +140,9 @@ def mesh(
                     offset.append(pnts.shape[0])
                     length.append(pnt.shape[0])
                     pnts = np.vstack((pnts, pnt))
-                # generate pos for __call__
+                # generate points for __call__
                 pnts = pnts.T[select]
-                out = self.unstructured(pos=pnts, **kwargs)
+                out = self.unstructured(points=pnts, **kwargs)
                 if isinstance(out, np.ndarray):
                     field = out
                 else:
@@ -145,7 +153,7 @@ def mesh(
                     field_dict[cell] = field[offset[i] : offset[i] + length[i]]
                 mesh.cell_data[name] = field_dict
             else:
-                out = self.unstructured(pos=mesh.points.T[select], **kwargs)
+                out = self.unstructured(points=mesh.points.T[select], **kwargs)
                 if isinstance(out, np.ndarray):
                     field = out
                 else:
@@ -154,13 +162,13 @@ def mesh(
                 mesh.point_data[name] = field
         return out
 
-    def _pre_pos(self, pos, mesh_type="unstructured", make_unstruct=False):
+    def _pre_points(self, points, mesh_type="unstructured", make_unstruct=False):
         """
         Preprocessing positions and mesh_type.
 
         Parameters
         ----------
-        pos : :any:`iterable`
+        points : :any:`iterable`
             the position tuple, containing main direction and transversal
             directions
         mesh_type : :class:`str`
@@ -176,7 +184,7 @@ def _pre_pos(self, pos, mesh_type="unstructured", make_unstruct=False):
             analog to x
         z : :class:`numpy.ndarray` or None
             analog to x
-        pos : :class:`tuple` of :class:`numpy.ndarray`
+        points : :class:`tuple` of :class:`numpy.ndarray`
             the normalized position tuple
         mesh_type_gen : :class:`str`
             'structured' / 'unstructured' for the generator
@@ -185,11 +193,11 @@ def _pre_pos(self, pos, mesh_type="unstructured", make_unstruct=False):
         axis_lens : :class:`tuple` or :any:`None`
             axis lengths of the structured mesh if mesh type was changed.
         """
-        x, y, z = pos2xyz(pos, max_dim=self.model.dim)
-        pos = xyz2pos(x, y, z)
+        x, y, z = pos2xyz(points, max_dim=self.model.dim)
+        points = xyz2pos(x, y, z)
         mesh_type_gen = mesh_type
         # format the positional arguments of the mesh
-        check_mesh(self.model.dim, x, y, z, mesh_type)
+        check_mesh(self.model.dim, points, mesh_type)
         mesh_type_changed = False
         axis_lens = None
         if (
@@ -204,117 +212,7 @@ def _pre_pos(self, pos, mesh_type="unstructured", make_unstruct=False):
             x, y, z = unrotate_mesh(self.model.dim, self.model.angles, x, y, z)
         if not self.model.is_isotropic:
             y, z = make_isotropic(self.model.dim, self.model.anis, y, z)
-        return x, y, z, pos, mesh_type_gen, mesh_type_changed, axis_lens
-
-    def _to_vtk_helper(
-        self, filename=None, field_select="field", fieldname="field"
-    ):  # pragma: no cover
-        """Create a VTK/PyVista grid of the field or save it as a VTK file.
-
-        This is an internal helper that will handle saving or creating objects
-
-        Parameters
-        ----------
-        filename : :class:`str`
-            Filename of the file to be saved, including the path. Note that an
-            ending (.vtr or .vtu) will be added to the name. If ``None`` is
-            passed, a PyVista dataset of the appropriate type will be returned.
-        field_select : :class:`str`, optional
-            Field that should be stored. Can be:
-            "field", "raw_field", "krige_field", "err_field" or "krige_var".
-            Default: "field"
-        fieldname : :class:`str`, optional
-            Name of the field in the VTK file. Default: "field"
-        """
-        if self.value_type is None:
-            raise ValueError(
-                "Field value type not set! "
-                + "Specify 'scalar' or 'vector' before plotting."
-            )
-        elif self.value_type == "vector":
-            if hasattr(self, field_select):
-                field = getattr(self, field_select)
-            else:
-                field = None
-            if not (
-                self.pos is None or field is None or self.mesh_type is None
-            ):
-                suf = ["_X", "_Y", "_Z"]
-                fields = {}
-                for i in range(self.model.dim):
-                    fields[fieldname + suf[i]] = field[i]
-                if filename is None:
-                    return to_vtk(self.pos, fields, self.mesh_type)
-                else:
-                    return vtk_export(
-                        filename, self.pos, fields, self.mesh_type
-                    )
-        elif self.value_type == "scalar":
-            if hasattr(self, field_select):
-                field = getattr(self, field_select)
-            else:
-                field = None
-            if not (
-                self.pos is None or field is None or self.mesh_type is None
-            ):
-                if filename is None:
-                    return to_vtk(self.pos, {fieldname: field}, self.mesh_type)
-                else:
-                    return vtk_export(
-                        filename, self.pos, {fieldname: field}, self.mesh_type
-                    )
-            else:
-                print(
-                    "Field.to_vtk: No "
-                    + field_select
-                    + " stored in the class."
-                )
-        else:
-            raise ValueError(
-                "Unknown field value type: {}".format(self.value_type)
-            )
-
-    def to_pyvista(
-        self, field_select="field", fieldname="field"
-    ):  # pragma: no cover
-        """Create a VTK/PyVista grid of the stored field.
-
-        Parameters
-        ----------
-        field_select : :class:`str`, optional
-            Field that should be stored. Can be:
-            "field", "raw_field", "krige_field", "err_field" or "krige_var".
-            Default: "field"
-        fieldname : :class:`str`, optional
-            Name of the field in the VTK file. Default: "field"
-        """
-        grid = self._to_vtk_helper(
-            filename=None, field_select=field_select, fieldname=fieldname
-        )
-        return grid
-
-    def vtk_export(
-        self, filename, field_select="field", fieldname="field"
-    ):  # pragma: no cover
-        """Export the stored field to vtk.
-
-        Parameters
-        ----------
-        filename : :class:`str`
-            Filename of the file to be saved, including the path. Note that an
-            ending (.vtr or .vtu) will be added to the name.
-        field_select : :class:`str`, optional
-            Field that should be stored. Can be:
-            "field", "raw_field", "krige_field", "err_field" or "krige_var".
-            Default: "field"
-        fieldname : :class:`str`, optional
-            Name of the field in the VTK file. Default: "field"
-        """
-        if not isinstance(filename, str):
-            raise TypeError("Please use a string filename.")
-        return self._to_vtk_helper(
-            filename=filename, field_select=field_select, fieldname=fieldname
-        )
+        return x, y, z, points, mesh_type_gen, mesh_type_changed, axis_lens
 
     def plot(self, field="field", fig=None, ax=None):  # pragma: no cover
         """
@@ -382,18 +280,13 @@ def model(self, model):
                 "Field: 'model' is not an instance of 'gstools.CovModel'"
             )
 
-    @property
-    def value_type(self):
-        """:class:`str`: Type of the field values (scalar, vector)."""
-        return self._value_type
-
     def __str__(self):
         """Return String representation."""
         return self.__repr__()
 
     def __repr__(self):
         """Return String representation."""
-        return "Field(model={0}, mean={1})".format(self.model, self.mean)
+        return "Field(model={0})".format(self.model)
 
 
 if __name__ == "__main__":  # pragma: no cover
diff --git a/gstools/field/condition.py b/gstools/field/condition.py
index cbb0f076d..37aa857a9 100644
--- a/gstools/field/condition.py
+++ b/gstools/field/condition.py
@@ -40,7 +40,7 @@ def ordinary(srf):
     krige_ok = Ordinary(
         model=srf.model, cond_pos=srf.cond_pos, cond_val=srf.cond_val
     )
-    krige_field, krige_var = krige_ok(srf.pos, srf.mesh_type)
+    krige_field, krige_var = krige_ok(srf.points, srf.mesh_type)
 
     # evaluate the field at the conditional points
     x, y, z = pos2xyz(srf.cond_pos, max_dim=srf.model.dim)
@@ -52,7 +52,7 @@ def ordinary(srf):
     err_ok = Ordinary(
         model=srf.model, cond_pos=srf.cond_pos, cond_val=err_data
     )
-    err_field, __ = err_ok(srf.pos, srf.mesh_type)
+    err_field, __ = err_ok(srf.points, srf.mesh_type)
     cond_field = srf.raw_field + krige_field - err_field
     info = {"mean": krige_ok.mean}
     return cond_field, krige_field, err_field, krige_var, info
@@ -85,7 +85,7 @@ def simple(srf):
         cond_pos=srf.cond_pos,
         cond_val=srf.cond_val,
     )
-    krige_field, krige_var = krige_sk(srf.pos, srf.mesh_type)
+    krige_field, krige_var = krige_sk(srf.points, srf.mesh_type)
 
     # evaluate the field at the conditional points
     x, y, z = pos2xyz(srf.cond_pos, max_dim=srf.model.dim)
@@ -100,7 +100,7 @@ def simple(srf):
         cond_pos=srf.cond_pos,
         cond_val=err_data,
     )
-    err_field, __ = err_sk(srf.pos, srf.mesh_type)
+    err_field, __ = err_sk(srf.points, srf.mesh_type)
     cond_field = srf.raw_field + krige_field - err_field + srf.mean
     info = {}
     return cond_field, krige_field, err_field, krige_var, info
diff --git a/gstools/field/mesh.py b/gstools/field/mesh.py
new file mode 100644
index 000000000..e322d73a0
--- /dev/null
+++ b/gstools/field/mesh.py
@@ -0,0 +1,578 @@
+# -*- coding: utf-8 -*-
+"""
+GStools subpackage providing a base class for spatial fields.
+
+.. currentmodule:: gstools.field.base
+
+The following classes are provided
+
+.. autosummary::
+   Mesh
+"""
+# pylint: disable=C0103
+
+import numpy as np
+
+from pyevtk.hl import gridToVTK, pointsToVTK
+
+from gstools.tools.geometric import pos2xyz
+from gstools.field.tools import check_mesh, check_point_data
+
+__all__ = ["Mesh"]
+
+
+def value_type(mesh_type, shape):
+    """Determine the value type ("scalar" or "vector")"""
+    r = "scalar"
+    if mesh_type == "unstructured":
+        # TODO is this the right place for doing these checks?!
+        if len(shape) == 2 and 2 <= shape[0] <= 3:
+            r = "vector"
+    else:
+        # for very small (2-3 points) meshes, this could break
+        # a 100% solution would require dim. information.
+        if len(shape) == shape[0] + 1:
+            r = "vector"
+    return r
+
+# TODO figure out what this is all about
+def convert_points(dim, points):
+    points = np.array(points)
+    shape = points.shape
+    if dim == 1:
+        #points = points.flatten()
+        if len(points.shape) > 1:
+            raise ValueError(
+                "points shape {} does not match dim {}".format(shape, dim)
+            )
+    else:
+        pass
+
+    return points
+
+
+class Mesh:
+    """A base class encapsulating field data.
+
+    It holds a position array, which define the spatial locations of the
+    field values.
+    It can hold multiple fields in the :any:`self.point_data` dict. This assumes
+    that each field is defined on the same positions.
+    The mesh type must also be specified.
+
+    Parameters
+    ----------
+    pos : :class:`numpy.ndarray`, optional
+        positions of the field values
+    name : :any:`str`, optional
+        key of the field values
+    values : :any:`list`, optional
+        a list of the values of the fields
+    mesh_type : :any:`str`, optional
+        the type of mesh on which the field is defined on, can be
+        * unstructured
+        * structured
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> from gstools import Mesh
+    >>> pos = np.linspace(0., 100., 40)
+    >>> z = np.random.random(40)
+    >>> z2 = np.random.random(40)
+    >>> mesh = Mesh(1, (pos,))
+    >>> mesh.add_field(z, "test_field1")
+    >>> mesh.add_field(z2, "test_field2")
+    >>> mesh.set_default_field("test_field2")
+    >>> print(mesh.field)
+
+    """
+    def __init__(
+        self, dim, points=None, name="field", values=None, mesh_type="unstructured",
+    ):
+        self._dim = dim
+
+        # the pos/ points of the mesh
+        if points is not None:
+            check_mesh(dim, points, mesh_type)
+        self._points = points
+
+        # data stored at each pos/ point, the "fields"
+        if values is not None:
+            self.point_data = {name: values}
+        else:
+            self.point_data = {}
+
+        # data valid for the global field
+        self.field_data = {}
+
+        # do following line manually in order to satisfy the linters
+        # self.set_field_data(name, "default_field")
+        self.default_field = name
+        self.field_data["default_field"] = name
+
+        self.mesh_type = mesh_type
+
+    def set_field_data(self, value, name):
+        """Add an attribute to this instance and add it the `field_data`
+
+        This helper method is used to create attributes for easy access
+        while using this class, but it also adds an entry to the dictionary
+        `field_data`, which is used for exporting the data.
+        """
+        setattr(self, name, value)
+        self.field_data[name] = value
+
+    def add_field(
+        self, values, name="field", is_default_field=False,
+    ):
+        """Add a field (point data) to the mesh
+
+        .. note::
+            If no field has existed before calling this method,
+            the given field will be set to the default field.
+
+        .. warning::
+            If point data with the same `name` already exists, it will be
+            overwritten.
+
+        Parameters
+        ----------
+        values : :class:`numpy.ndarray`
+            the point data, has to be the same shape as the mesh
+        name : :class:`str`, optional
+            the name of the point data
+        is_default_field : :class:`bool`, optional
+            is this the default field of the mesh?
+
+        """
+        values = np.array(values)
+        check_point_data(
+            self.dim,
+            self.points,
+            values,
+            self.mesh_type,
+            value_type(self.mesh_type, values.shape),
+        )
+        # set the default_field to the first field added
+        if not self.point_data:
+            is_default_field = True
+        elif self.point_data.get(self.default_field, False) is None:
+            del(self.point_data[self.default_field])
+            is_default_field = True
+        if is_default_field:
+            self.set_field_data(name, "default_field")
+        self.point_data[name] = values
+
+    def del_field_data(self):
+        """Delete the field data.
+
+        Deleting the field data resets the field data dictionary and deletes
+        the attributes too.
+        """
+        for attr in self.field_data.keys():
+            if attr == "default_field" or attr == "mesh_type":
+                continue
+            try:
+                delattr(self, attr)
+            except AttributeError:
+                pass
+        self.field_data = {"default_field": "field", "mesh_type": self.mesh_type}
+
+    def __getitem__(self, key):
+        """:any:`numpy.ndarray`: The values of the field."""
+        return self.point_data[key]
+
+    def __setitem__(self, key, value):
+        self.point_data[key] = value
+
+    @property
+    def dim(self):
+        """:any:`int`: The dimension of the mesh."""
+        return self._dim
+
+    @property
+    def points(self):
+        """:any:`numpy.ndarray`: The pos. on which the field is defined."""
+        return self._points
+
+    @points.setter
+    def points(self, value):
+        """
+        Warning: setting new positions deletes all previously stored fields.
+        """
+        self.point_data = {self.default_field: None}
+        check_mesh(self.dim, value, self.mesh_type)
+        self._points = value
+
+    @property
+    def field(self):
+        """:class:`numpy.ndarray`: The point data of the default field."""
+        return self.point_data[self.default_field]
+
+    @field.setter
+    def field(self, values):
+        check_point_data(
+            self.dim,
+            self.points,
+            values,
+            self.mesh_type,
+            value_type(self.mesh_type, values.shape),
+        )
+        self.point_data[self.default_field] = values
+
+    @property
+    def value_type(self):
+        """:any:`str`: The value type of the default field."""
+        if (
+            self.default_field in self.point_data and
+            self.point_data[self.default_field] is not None
+        ):
+            r = value_type(
+                self.mesh_type, self.point_data[self.default_field].shape
+            )
+        else:
+            r = None
+        return r
+
+    @property
+    def mesh_type(self):
+        """:any:`str`: The mesh type of the fields."""
+        return self._mesh_type
+
+    @mesh_type.setter
+    def mesh_type(self, value):
+        """
+        Warning: setting a new mesh type deletes all previously stored fields.
+        """
+        self._check_mesh_type(value)
+        if value == "structured":
+            self._vtk_export_fct = self._vtk_export_structured
+            self._to_vtk_fct = self._to_vtk_structured
+        else:
+            self._vtk_export_fct = self._vtk_export_unstructured
+            self._to_vtk_fct = self._to_vtk_unstructured
+        self.point_data = {}
+        self._mesh_type = value
+        self.field_data["mesh_type"] = value
+
+    def _vtk_naming_helper(
+        self, field_select="field", fieldname="field"
+    ):  # pragma: no cover
+        """Prepare the field names for export.
+
+        Parameters
+        ----------
+        field_select : :class:`str`, optional
+            Field that should be stored. Can be:
+            "field", "raw_field", "krige_field", "err_field" or "krige_var".
+            Default: "field"
+        fieldname : :class:`str`, optional
+            Name of the field in the VTK file. Default: "field"
+
+        Returns
+        -------
+        fields : :class:`dict`
+            a dictionary containing the fields to be exported
+        """
+        if self.value_type is None:
+            raise ValueError(
+                "Field value type not set! "
+                + "Specify 'scalar' or 'vector' before plotting."
+            )
+        elif self.value_type == "vector":
+            if hasattr(self, field_select):
+                field = getattr(self, field_select)
+            else:
+                field = None
+            if not (
+                self.points is None or field is None or self.mesh_type is None
+            ):
+                suf = ["_X", "_Y", "_Z"]
+                fields = {}
+                for i in range(self.dim):
+                    fields[fieldname + suf[i]] = field[i]
+                return fields
+        elif self.value_type == "scalar":
+            if hasattr(self, field_select):
+                field = getattr(self, field_select)
+            else:
+                field = None
+            if not (
+                self.points is None or field is None or self.mesh_type is None
+            ):
+                return {fieldname: field}
+            else:
+                print(
+                    "Field.to_vtk: No "
+                    + field_select
+                    + " stored in the class."
+                )
+        else:
+            raise ValueError(
+                "Unknown field value type: {}".format(self.value_type)
+            )
+
+    def convert(
+        self, data_format="pyvista", field_select="field", fieldname="field"
+    ):  # pragma: no cover
+        """Create a VTK/PyVista grid of the stored field.
+
+        Parameters
+        ----------
+        data_format : :class:`str`
+            the format in which to convert the data, possible choices:
+            * 'pyvista'
+        field_select : :class:`str`, optional
+            Field that should be stored. Can be:
+            "field", "raw_field", "krige_field", "err_field" or "krige_var".
+            Default: "field"
+        fieldname : :class:`str`, optional
+            Name of the field in the VTK file. Default: "field"
+        """
+        if data_format.lower() != "pyvista":
+            raise NotImplementedError(
+                "Only pyvista format is supported at the moment."
+            )
+        field_names = self._vtk_naming_helper(
+            field_select=field_select, fieldname=fieldname
+        )
+
+        grid = self._to_vtk_fct(field_names)
+        return grid
+
+    def export(
+        self, filename, data_format="vtk", field_select="field", fieldname="field"
+    ):  # pragma: no cover
+        """Export the stored field to vtk.
+
+        Parameters
+        ----------
+        filename : :class:`str`
+            Filename of the file to be saved, including the path. Note that an
+            ending (.vtr or .vtu) will be added to the name.
+        data_format : :class:`str`
+            the format in which to export the data, possible choices:
+            * 'vtk'
+        field_select : :class:`str`, optional
+            Field that should be stored. Can be:
+            "field", "raw_field", "krige_field", "err_field" or "krige_var".
+            Default: "field"
+        fieldname : :class:`str`, optional
+            Name of the field in the VTK file. Default: "field"
+        """
+        if data_format.lower() != "vtk":
+            raise NotImplementedError(
+                "Only VTK format is supported at the moment."
+            )
+        if not isinstance(filename, str):
+            raise TypeError("Please use a string as a filename.")
+        field_names = self._vtk_naming_helper(
+            field_select=field_select, fieldname=fieldname
+        )
+        return self._vtk_export_fct(filename, field_names)
+
+    def _vtk_export_unstructured(self, filename, fields):
+        """Export a field to vtk unstructured grid file.
+
+        Parameters
+        ----------
+        filename : :class:`str`
+            Filename of the file to be saved, including the path. Note that an
+            ending (.vtu) will be added to the name.
+        fields : :class:`dict` or :class:`numpy.ndarray`
+            Unstructured fields to be saved.
+            Either a single numpy array as returned by SRF,
+            or a dictionary of fields with theirs names as keys.
+        """
+        x, y, z, fields = self._vtk_unstructured_reshape(fields)
+        return pointsToVTK(filename, x, y, z, data=fields)
+
+    def _vtk_export_structured(self, filename, fields):
+        """Export a field to vtk structured rectilinear grid file.
+
+        Parameters
+        ----------
+        filename : :class:`str`
+            Filename of the file to be saved, including the path. Note that an
+            ending (.vtr) will be added to the name.
+        fields : :class:`dict` or :class:`numpy.ndarray`
+            Structured fields to be saved.
+            Either a single numpy array as returned by SRF,
+            or a dictionary of fields with theirs names as keys.
+        """
+        x, y, z, fields = self._vtk_structured_reshape(fields=fields)
+        return gridToVTK(filename, x, y, z, pointData=fields)
+
+    def _to_vtk_unstructured(self, fields):  # pragma: no cover
+        """Export a field to vtk structured rectilinear grid file.
+
+        Parameters
+        ----------
+        fields : :class:`dict` or :class:`numpy.ndarray`
+            Unstructured fields to be saved.
+            Either a single numpy array as returned by SRF,
+            or a dictionary of fields with theirs names as keys.
+
+        Returns
+        -------
+        :class:`pyvista.UnstructuredGrid`
+            A PyVista unstructured grid of the unstructured field data. Data arrays
+            live on the point data of this PyVista dataset. This is essentially
+            a point cloud with no topology.
+        """
+        x, y, z, fields = self._vtk_unstructured_reshape(fields)
+        try:
+            import pyvista as pv
+
+            grid = pv.PolyData(np.c_[x, y, z]).cast_to_unstructured_grid()
+            grid.point_arrays.update(fields)
+        except ImportError:
+            raise ImportError("Please install PyVista to create VTK datasets.")
+        return grid
+
+    def _to_vtk_structured(self, fields):  # pragma: no cover
+        """Create a vtk structured rectilinear grid from a field.
+
+        Parameters
+        ----------
+        fields : :class:`dict` or :class:`numpy.ndarray`
+            Structured fields to be saved.
+            Either a single numpy array as returned by SRF,
+            or a dictionary of fields with theirs names as keys.
+
+        Returns
+        -------
+        :class:`pyvista.RectilinearGrid`
+            A PyVista rectilinear grid of the structured field data. Data arrays
+            live on the point data of this PyVista dataset.
+        """
+        x, y, z, fields = self._vtk_structured_reshape(fields=fields)
+        try:
+            import pyvista as pv
+
+            grid = pv.RectilinearGrid(x, y, z)
+            grid.point_arrays.update(fields)
+        except ImportError:
+            raise ImportError("Please install PyVista to create VTK datasets.")
+        return grid
+
+    def _vtk_unstructured_reshape(self, fields):
+        if not isinstance(fields, dict):
+            fields = {"field": fields}
+        x, y, z = pos2xyz(self.points)
+        if y is None:
+            y = np.zeros_like(x)
+        if z is None:
+            z = np.zeros_like(x)
+        for field in fields:
+            fields[field] = fields[field].reshape(-1)
+            if (
+                len(fields[field]) != len(x)
+                or len(fields[field]) != len(y)
+                or len(fields[field]) != len(z)
+            ):
+                raise ValueError(
+                    "gstools.vtk_export_unstructured: "
+                    "field shape doesn't match the given mesh"
+                )
+        return x, y, z, fields
+
+    def _vtk_structured_reshape(self, fields):
+        """An internal helper to extract what is needed for the vtk rectilinear grid
+        """
+        if not isinstance(fields, dict):
+            fields = {"field": fields}
+        x, y, z = pos2xyz(self.points)
+        if y is None:
+            y = np.array([0])
+        if z is None:
+            z = np.array([0])
+        # need fortran order in VTK
+        for field in fields:
+            fields[field] = fields[field].reshape(-1, order="F")
+            if len(fields[field]) != len(x) * len(y) * len(z):
+                raise ValueError(
+                    "gstools.vtk_export_structured: "
+                    "field shape doesn't match the given mesh"
+                )
+        return x, y, z, fields
+
+    def _check_mesh_type(self, mesh_type):
+        if mesh_type != "unstructured" and mesh_type != "structured":
+            raise ValueError("Unknown 'mesh_type': {}".format(mesh_type))
+
+
+#class Mesh(MeshBase):
+#    def __init__(self):
+#        super().__init__()
+#
+#    def _vtk_reshape(self, fields):
+#        if not isinstance(fields, dict):
+#            fields = {"field": fields}
+#        x, y, z = pos2xyz(self.pos)
+#        if y is None:
+#            y = np.zeros_like(x)
+#        if z is None:
+#            z = np.zeros_like(x)
+#        for field in fields:
+#            fields[field] = fields[field].reshape(-1)
+#            if (
+#                len(fields[field]) != len(x)
+#                or len(fields[field]) != len(y)
+#                or len(fields[field]) != len(z)
+#            ):
+#                raise ValueError(
+#                    "gstools.vtk_export_unstructured: "
+#                    "field shape doesn't match the given mesh"
+#                )
+#        return x, y, z, fields
+#
+#    def _vtk_export(self, filename, fields):
+#        """Export a field to vtk unstructured grid file.
+#
+#        Parameters
+#        ----------
+#        filename : :class:`str`
+#            Filename of the file to be saved, including the path. Note that an
+#            ending (.vtu) will be added to the name.
+#        fields : :class:`dict` or :class:`numpy.ndarray`
+#            Unstructured fields to be saved.
+#            Either a single numpy array as returned by SRF,
+#            or a dictionary of fields with theirs names as keys.
+#        """
+#        x, y, z, fields = self._vtk_reshape(fields=fields)
+#        return pointsToVTK(filename, x, y, z, data=fields)
+#
+#    def _to_vtk(self, fields):
+#        """Export a field to vtk structured rectilinear grid file.
+#
+#        Parameters
+#        ----------
+#        fields : :class:`dict` or :class:`numpy.ndarray`
+#            Unstructured fields to be saved.
+#            Either a single numpy array as returned by SRF,
+#            or a dictionary of fields with theirs names as keys.
+#
+#        Returns
+#        -------
+#        :class:`pyvista.UnstructuredGrid`
+#            A PyVista unstructured grid of the unstructured field data. Data arrays
+#            live on the point data of this PyVista dataset. This is essentially
+#            a point cloud with no topology.
+#        """
+#        x, y, z, fields = self._vtk_unstructured_reshape(fields=fields)
+#        try:
+#            import pyvista as pv
+#
+#            grid = pv.PolyData(np.c_[x, y, z]).cast_to_unstructured_grid()
+#            grid.point_arrays.update(fields)
+#        except ImportError:
+#            raise ImportError("Please install PyVista to create VTK datasets.")
+#        return grid
+
+
+if __name__ == "__main__":  # pragma: no cover
+    import doctest
+
+    doctest.testmod()
diff --git a/gstools/field/srf.py b/gstools/field/srf.py
index 74b6ea4c6..c5fda834f 100644
--- a/gstools/field/srf.py
+++ b/gstools/field/srf.py
@@ -11,6 +11,7 @@
 """
 # pylint: disable=C0103
 
+from warnings import warn
 import numpy as np
 from gstools.field.generator import RandMeth, IncomprRandMeth
 from gstools.field.tools import reshape_field_from_unstruct_to_struct
@@ -77,25 +78,23 @@ class SRF(Field):
     def __init__(
         self,
         model,
+        points=None,
+        mesh_type="unstructured",
         mean=0.0,
         upscaling="no_scaling",
         generator="RandMeth",
         **generator_kwargs
     ):
-        super().__init__(model, mean)
+        super().__init__(model, points=points, mesh_type=mesh_type, mean=mean)
         # initialize private attributes
         self._generator = None
         self._upscaling = None
         self._upscaling_func = None
         # condition related
+        # TODO rename to self._cond_points
         self._cond_pos = None
         self._cond_val = None
         self._krige_type = None
-        # initialize attributes
-        self.raw_field = None
-        self.krige_field = None
-        self.err_field = None
-        self.krige_var = None
         self.set_generator(generator, **generator_kwargs)
         self.upscaling = upscaling
         if self._value_type is None:
@@ -105,7 +104,12 @@ def __init__(
             )
 
     def __call__(
-        self, pos, seed=np.nan, point_volumes=0.0, mesh_type="unstructured"
+        self,
+        points=None,
+        name="field",
+        seed=np.nan,
+        point_volumes=0.0,
+        mesh_type=None,
     ):
         """Generate the spatial random field.
 
@@ -113,7 +117,7 @@ def __call__(
 
         Parameters
         ----------
-        pos : :class:`list`
+        points : :class:`list`, optional
             the position tuple, containing main direction and transversal
             directions
         seed : :class:`int`, optional
@@ -132,20 +136,33 @@ def __call__(
         field : :class:`numpy.ndarray`
             the SRF
         """
-        self.mesh_type = mesh_type
+        if mesh_type is not None:
+            warn(
+                "'mesh_type' will be ignored in the srf.__call__, use " +
+                "__init__ or setter instead.",
+                DeprecationWarning,
+            )
+            self.mesh_type = mesh_type
+        if points is not None:
+            warn(
+                "'points' will be ignored in the srf.__call__, use "
+                "__init__ or setter instead.",
+                DeprecationWarning,
+            )
+            self.points = points
         # update the model/seed in the generator if any changes were made
         self.generator.update(self.model, seed)
         # internal conversation
-        x, y, z, self.pos, mt_gen, mt_changed, axis_lens = self._pre_pos(
-            pos, mesh_type
+        x, y, z, points, mt_gen, mt_changed, axis_lens = self._pre_points(
+            self.points, self.mesh_type
         )
         # generate the field
         self.raw_field = self.generator.__call__(x, y, z, mt_gen)
         # reshape field if we got an unstructured mesh
         if mt_changed:
             self.raw_field = reshape_field_from_unstruct_to_struct(
-                self.model.dim, self.raw_field, axis_lens
-            )
+                    self.model.dim, self.raw_field, axis_lens
+                )
         # apply given conditions to the field
         if self.condition:
             (
@@ -156,21 +173,23 @@ def __call__(
                 info,
             ) = self.cond_func(self)
             # store everything in the class
-            self.field = cond_field
-            self.krige_field = krige_field
-            self.err_field = err_field
-            self.krige_var = krigevar
+            self.add_field(cond_field, name)
+            self.add_field(krige_field, "krige_field")
+            self.add_field(err_field, "err_field")
+            self.set_field_data(krigevar, "krige_var")
             if "mean" in info:  # ordinary krging estimates mean
-                self.mean = info["mean"]
+                self.set_field_data(info["mean"], "mean")
         else:
-            self.field = self.raw_field + self.mean
+            self.add_field(self.raw_field + self.mean, name)
         # upscaled variance
         if not np.isscalar(point_volumes) or not np.isclose(point_volumes, 0):
             scaled_var = self.upscaling_func(self.model, point_volumes)
-            self.field -= self.mean
-            self.field *= np.sqrt(scaled_var / self.model.sill)
-            self.field += self.mean
-        return self.field
+            self[name] -= self.mean
+            self[name] *= np.sqrt(scaled_var / self.model.sill)
+            self[name] += self.mean
+        self.add_field(self.raw_field, "raw_field")
+        del(self.raw_field)
+        return self[name]
 
     def set_condition(
         self, cond_pos=None, cond_val=None, krige_type="ordinary"
diff --git a/gstools/field/tools.py b/gstools/field/tools.py
index e9f289a23..4de6a7f14 100644
--- a/gstools/field/tools.py
+++ b/gstools/field/tools.py
@@ -74,34 +74,34 @@ def reshape_input_axis_from_struct(x, y=None, z=None):
 # SRF helpers #################################################################
 
 
-def check_mesh(dim, x, y, z, mesh_type):
+def check_mesh(dim, points, mesh_type):
     """Do a basic check of the shapes of the input arrays."""
     if dim >= 2:
-        if y is None:
+        if len(points) < 2:
             raise ValueError(
                 "The y-component is missing for " "{0} dimensions".format(dim)
             )
     if dim == 3:
-        if z is None:
+        if len(points) < 3:
             raise ValueError(
                 "The z-component is missing for " "{0} dimensions".format(dim)
             )
     if mesh_type == "unstructured":
         if dim >= 2:
             try:
-                if len(x) != len(y):
+                if len(points[0]) != len(points[1]):
                     raise ValueError(
                         "len(x) = {0} != len(y) = {1} "
-                        "for unstructured grids".format(len(x), len(y))
+                        "for unstructured grids".format(len(points[0]), len(points[1]))
                     )
             except TypeError:
                 pass
             if dim == 3:
                 try:
-                    if len(x) != len(z):
+                    if len(points[0]) != len(points[2]):
                         raise ValueError(
                             "len(x) = {0} != len(z) = {1} "
-                            "for unstructured grids".format(len(x), len(z))
+                            "for unstructured grids".format(len(points[0]), len(points[2]))
                         )
                 except TypeError:
                     pass
@@ -110,6 +110,40 @@ def check_mesh(dim, x, y, z, mesh_type):
     else:
         raise ValueError("Unknown mesh type {0}".format(mesh_type))
 
+def check_point_data(dim, points, data, mesh_type, value_type='scalar'):
+    """Do a basic check and compare shapes of the field data to mesh field."""
+    offset = 0
+    if value_type == 'vector':
+        if dim != data.shape[0]:
+            raise ValueError(
+                "Field is {} dim., but mesh is {}".
+                format(data.shape[0], dim)
+            )
+        offset = 1
+    if mesh_type == 'unstructured':
+        for d in range(dim):
+            if points[d].shape[0] != data.shape[offset]:
+                raise ValueError(
+                    "Field shape {} does not match existing mesh shape {}".
+                    format(data.shape[offset], points[d].shape[0])
+                )
+    elif mesh_type == 'structured':
+        try:
+            if dim != len(data.squeeze().shape)-offset:
+                raise ValueError(
+                    "Field is {} dim., but mesh {}".
+                    format(len(data.squeeze().shape), dim)
+                )
+        except AttributeError:
+            pass
+        for d in range(dim):
+            if points[d].shape[0] != data.shape[offset+d]:
+                raise ValueError(
+                    "Field shape {} does not match existing mesh shape {}".
+                    format(data.shape[offset+d], points[d].shape[0])
+                )
+    else:
+        raise ValueError("Unknown mesh type {0}".format(mesh_type))
 
 def make_isotropic(dim, anis, y, z):
     """Stretch given axes in order to implement anisotropy."""
diff --git a/gstools/krige/base.py b/gstools/krige/base.py
index 366d02cd3..32e9d5996 100755
--- a/gstools/krige/base.py
+++ b/gstools/krige/base.py
@@ -69,7 +69,7 @@ def __init__(
         drift_functions=None,
         trend_function=None,
     ):
-        super().__init__(model, mean)
+        super().__init__(model, mean=mean)
         self.krige_var = None
         # initialize private attributes
         self._unbiased = True
@@ -88,7 +88,7 @@ def __init__(
         self.set_condition(cond_pos, cond_val, ext_drift)
 
     def __call__(
-        self, pos, mesh_type="unstructured", ext_drift=None, chunk_size=None
+        self, points, mesh_type="unstructured", ext_drift=None, chunk_size=None
     ):
         """
         Generate the kriging field.
@@ -98,7 +98,7 @@ def __call__(
 
         Parameters
         ----------
-        pos : :class:`list`
+        points : :class:`list`
             the position tuple, containing main direction and transversal
             directions (x, [y, z])
         mesh_type : :class:`str`, optional
@@ -119,8 +119,8 @@ def __call__(
         """
         self.mesh_type = mesh_type
         # internal conversation
-        x, y, z, self.pos, __, mt_changed, axis_lens = self._pre_pos(
-            pos, mesh_type, make_unstruct=True
+        x, y, z, self.points, __, mt_changed, axis_lens = self._pre_points(
+            points, mesh_type, make_unstruct=True
         )
         point_no = len(x)
         # set chunk size
@@ -215,7 +215,7 @@ def _post_field(self, field, krige_var):
             self.field = field
         else:
             self.field = field + eval_func(
-                self.trend_function, self.pos, self.mesh_type
+                self.trend_function, self.points, self.mesh_type
             )
         self.krige_var = krige_var
 
@@ -309,7 +309,7 @@ def set_drift_functions(self, drift_functions=None):
 
     def update(self):
         """Update the kriging settings."""
-        x, y, z, __, __, __, __ = self._pre_pos(self.cond_pos)
+        x, y, z, __, __, __, __ = self._pre_points(self.cond_pos)
         # krige pos are the unrotated and isotropic condition positions
         self._krige_pos = (x, y, z)[: self.model.dim]
         self._krige_mat = self._get_krige_mat()
diff --git a/gstools/krige/methods.py b/gstools/krige/methods.py
index bad1b059d..2be07b64b 100644
--- a/gstools/krige/methods.py
+++ b/gstools/krige/methods.py
@@ -83,7 +83,7 @@ def _post_field(self, field, krige_var):
             self.field = (
                 field
                 + self.mean
-                + eval_func(self.trend_function, self.pos, self.mesh_type)
+                + eval_func(self.trend_function, self.points, self.mesh_type)
             )
         # add the given mean
         self.krige_var = self.model.sill - krige_var
diff --git a/gstools/tools/export.py b/gstools/tools/export.py
index 9174830b9..3d69e5d3f 100644
--- a/gstools/tools/export.py
+++ b/gstools/tools/export.py
@@ -16,6 +16,7 @@
 """
 # pylint: disable=C0103, E1101
 
+import warnings
 import numpy as np
 from pyevtk.hl import gridToVTK, pointsToVTK
 
@@ -39,7 +40,7 @@
 # export routines #############################################################
 
 
-def _vtk_structured_helper(pos, fields):
+def _vtk_structured_reshape(pos, fields):
     """An internal helper to extract what is needed for the vtk rectilinear grid
     """
     if not isinstance(fields, dict):
@@ -79,7 +80,11 @@ def to_vtk_structured(pos, fields):  # pragma: no cover
         A PyVista rectilinear grid of the structured field data. Data arrays
         live on the point data of this PyVista dataset.
     """
-    x, y, z, fields = _vtk_structured_helper(pos=pos, fields=fields)
+    warnings.warn(
+        "Export routines are deprecated, use Mesh export methods instead",
+        DeprecationWarning
+    )
+    x, y, z, fields = _vtk_structured_reshape(pos=pos, fields=fields)
     try:
         import pyvista as pv
 
@@ -106,11 +111,15 @@ def vtk_export_structured(filename, pos, fields):  # pragma: no cover
         Either a single numpy array as returned by SRF,
         or a dictionary of fields with theirs names as keys.
     """
-    x, y, z, fields = _vtk_structured_helper(pos=pos, fields=fields)
+    warnings.warn(
+        "Export routines are deprecated, use Mesh export methods instead",
+        DeprecationWarning
+    )
+    x, y, z, fields = _vtk_structured_reshape(pos=pos, fields=fields)
     return gridToVTK(filename, x, y, z, pointData=fields)
 
 
-def _vtk_unstructured_helper(pos, fields):
+def _vtk_unstructured_reshape(pos, fields):
     if not isinstance(fields, dict):
         fields = {"field": fields}
     x, y, z = pos2xyz(pos)
@@ -152,7 +161,11 @@ def to_vtk_unstructured(pos, fields):  # pragma: no cover
         live on the point data of this PyVista dataset. This is essentially
         a point cloud with no topology.
     """
-    x, y, z, fields = _vtk_unstructured_helper(pos=pos, fields=fields)
+    warnings.warn(
+        "Export routines are deprecated, use Mesh export methods instead",
+        DeprecationWarning
+    )
+    x, y, z, fields = _vtk_unstructured_reshape(pos=pos, fields=fields)
     try:
         import pyvista as pv
 
@@ -179,7 +192,11 @@ def vtk_export_unstructured(filename, pos, fields):  # pragma: no cover
         Either a single numpy array as returned by SRF,
         or a dictionary of fields with theirs names as keys.
     """
-    x, y, z, fields = _vtk_unstructured_helper(pos=pos, fields=fields)
+    warnings.warn(
+        "Export routines are deprecated, use Mesh export methods instead",
+        DeprecationWarning
+    )
+    x, y, z, fields = _vtk_unstructured_reshape(pos=pos, fields=fields)
     return pointsToVTK(filename, x, y, z, data=fields)
 
 
diff --git a/gstools/transform/field.py b/gstools/transform/field.py
index eb0a1954e..52c68ba50 100644
--- a/gstools/transform/field.py
+++ b/gstools/transform/field.py
@@ -335,12 +335,12 @@ def _zinnharvey(field, conn="high", mean=None, var=None):
         mean = np.mean(field)
     if var is None:
         var = np.var(field)
-    field = np.abs((field - mean) / np.sqrt(var))
+    field = np.abs((field - mean) / var)
     field = 2 * erf(field / np.sqrt(2)) - 1
     field = np.sqrt(2) * erfinv(field)
     if conn == "high":
         field = -field
-    return field * np.sqrt(var) + mean
+    return field * var + mean
 
 
 def _normal_force_moments(field, mean=0, var=1):
@@ -508,5 +508,5 @@ def _uniform_to_uquad(field, a=0, b=1):
     y_raw = 3 * field / al + ga
     out = np.zeros_like(y_raw)
     out[y_raw > 0] = y_raw[y_raw > 0] ** (1 / 3)
-    out[y_raw < 0] = -((-y_raw[y_raw < 0]) ** (1 / 3))
+    out[y_raw < 0] = -(-y_raw[y_raw < 0]) ** (1 / 3)
     return out + be
diff --git a/tests/test_condition.py b/tests/test_condition.py
index e7c588791..5544366b3 100644
--- a/tests/test_condition.py
+++ b/tests/test_condition.py
@@ -44,10 +44,11 @@ def test_simple(self):
             model = Model(
                 dim=1, var=0.5, len_scale=2, anis=[0.1, 1], angles=[0.5, 0, 0]
             )
+            # import pudb; pu.db
             srf = SRF(model, self.mean, seed=19970221)
             srf.set_condition(self.cond_pos[0], self.cond_val, "simple")
-            field_1 = srf.unstructured(self.pos[0])
-            field_2 = srf.structured(self.pos[0])
+            field_1 = srf.unstructured((self.pos[0],))
+            field_2 = srf.structured((self.pos[0],))
             for i, val in enumerate(self.cond_val):
                 self.assertAlmostEqual(val, field_1[i], places=2)
                 self.assertAlmostEqual(val, field_2[(i,)], places=2)
@@ -60,7 +61,7 @@ def test_simple(self):
                     anis=[0.1, 1],
                     angles=[0.5, 0, 0],
                 )
-                srf = SRF(model, self.mean, seed=19970221)
+                srf = SRF(model, mean=self.mean, seed=19970221)
                 srf.set_condition(self.cond_pos[:dim], self.cond_val, "simple")
                 field_1 = srf.unstructured(self.pos[:dim])
                 field_2 = srf.structured(self.pos[:dim])
@@ -73,10 +74,10 @@ def test_ordinary(self):
             model = Model(
                 dim=1, var=0.5, len_scale=2, anis=[0.1, 1], angles=[0.5, 0, 0]
             )
-            srf = SRF(model, seed=19970221)
+            srf = SRF(model, points=(self.pos[0],), seed=19970221)
             srf.set_condition(self.cond_pos[0], self.cond_val, "ordinary")
-            field_1 = srf.unstructured(self.pos[0])
-            field_2 = srf.structured(self.pos[0])
+            field_1 = srf.unstructured()
+            field_2 = srf.structured()
             for i, val in enumerate(self.cond_val):
                 self.assertAlmostEqual(val, field_1[i], places=2)
                 self.assertAlmostEqual(val, field_2[(i,)], places=2)
diff --git a/tests/test_export.py b/tests/test_export.py
index 8c483d092..c9f25231d 100644
--- a/tests/test_export.py
+++ b/tests/test_export.py
@@ -6,8 +6,7 @@
 import tempfile
 import shutil
 
-from gstools import SRF, Gaussian, Exponential
-from gstools.random import MasterRNG
+import gstools as gs
 
 HAS_PYVISTA = False
 try:
@@ -21,42 +20,78 @@
 class TestExport(unittest.TestCase):
     def setUp(self):
         self.test_dir = tempfile.mkdtemp()
-        # structured field with a size 100x100x100 and a grid-size of 1x1x1
-        x = y = z = range(50)
-        model = Gaussian(dim=3, var=0.6, len_scale=20)
-        self.srf_structured = SRF(model)
-        self.srf_structured((x, y, z), mesh_type="structured")
-        # unstrucutred field
-        seed = MasterRNG(19970221)
+        self.x_grid = self.y_grid = self.z_grid = np.arange(25)
+        model = gs.Gaussian(dim=3, var=0.6, len_scale=20)
+        self.srf_structured = gs.SRF(model)
+        self.srf_structured(
+            (self.x_grid, self.y_grid, self.z_grid), mesh_type="structured"
+        )
+        # unstructured field
+        seed = gs.random.MasterRNG(19970221)
         rng = np.random.RandomState(seed())
-        x = rng.randint(0, 100, size=1000)
-        y = rng.randint(0, 100, size=1000)
-        model = Exponential(
+        self.x_tuple = rng.randint(0, 100, size=100)
+        self.y_tuple = rng.randint(0, 100, size=100)
+        model = gs.Exponential(
             dim=2, var=1, len_scale=[12.0, 3.0], angles=np.pi / 8.0
         )
-        self.srf_unstructured = SRF(model, seed=20170519)
-        self.srf_unstructured([x, y])
+        self.srf_unstructured = gs.SRF(model, seed=20170519)
+        self.srf_unstructured([self.x_tuple, self.y_tuple])
 
     def tearDown(self):
         # Remove the test data directory after the test
         shutil.rmtree(self.test_dir)
 
     @unittest.skipIf(not HAS_PYVISTA, "PyVista is not installed.")
-    def test_pyvista(self):
-        mesh = self.srf_structured.to_pyvista()
+    def test_pyvista_struct(self):
+        mesh = self.srf_structured.convert()
+        self.assertIsInstance(mesh, pv.RectilinearGrid)
+
+    @unittest.skipIf(not HAS_PYVISTA, "PyVista is not installed.")
+    def test_pyvista_unstruct(self):
+        mesh = self.srf_unstructured.convert()
+        self.assertIsInstance(mesh, pv.UnstructuredGrid)
+
+    def test_pyevtk_export_struct(self):
+        filename = os.path.join(self.test_dir, "structured")
+        self.srf_structured.export(filename, "vtk")
+        self.assertTrue(os.path.isfile(filename + ".vtr"))
+
+    def test_pyevtk_export_unstruct(self):
+        filename = os.path.join(self.test_dir, "unstructured")
+        self.srf_unstructured.export(filename, "vtk")
+        self.assertTrue(os.path.isfile(filename + ".vtu"))
+
+    @unittest.skipIf(not HAS_PYVISTA, "PyVista is not installed.")
+    def test_pyvista_vector_struct(self):
+        model = gs.Gaussian(dim=2, var=1, len_scale=10)
+        srf = gs.SRF(model, generator="VectorField")
+        srf((self.x_grid, self.y_grid), mesh_type="structured", seed=19841203)
+        mesh = srf.convert()
         self.assertIsInstance(mesh, pv.RectilinearGrid)
-        mesh = self.srf_unstructured.to_pyvista()
+
+    @unittest.skipIf(not HAS_PYVISTA, "PyVista is not installed.")
+    def test_pyvista_vector_unstruct(self):
+        model = gs.Gaussian(dim=2, var=1, len_scale=10)
+        srf = gs.SRF(model, generator="VectorField")
+        srf((self.x_tuple, self.y_tuple), mesh_type="unstructured", seed=19841203)
+        mesh = srf.convert()
         self.assertIsInstance(mesh, pv.UnstructuredGrid)
 
-    def test_pyevtk_export(self):
-        # Structured
-        sfilename = os.path.join(self.test_dir, "structured")
-        self.srf_structured.vtk_export(sfilename)
-        self.assertTrue(os.path.isfile(sfilename + ".vtr"))
-        # Unstructured
-        ufilename = os.path.join(self.test_dir, "unstructured")
-        self.srf_unstructured.vtk_export(ufilename)
-        self.assertTrue(os.path.isfile(ufilename + ".vtu"))
+    def test_pyevtk_vector_export_struct(self):
+        filename = os.path.join(self.test_dir, "vector")
+        model = gs.Gaussian(dim=2, var=1, len_scale=10)
+        srf = gs.SRF(model, generator="VectorField")
+        srf((self.x_grid, self.y_grid), mesh_type="structured", seed=19841203)
+        srf.export(filename, "vtk")
+        self.assertTrue(os.path.isfile(filename + ".vtr"))
+
+    def test_pyevtk_vector_export_unstruct(self):
+        filename = os.path.join(self.test_dir, "vector")
+        model = gs.Gaussian(dim=2, var=1, len_scale=10)
+        srf = gs.SRF(model, generator="VectorField")
+        srf((self.x_tuple, self.y_tuple), mesh_type="unstructured", seed=19841203)
+        srf.export(filename, "vtk")
+        self.assertTrue(os.path.isfile(filename + ".vtu"))
 
 
 if __name__ == "__main__":
diff --git a/tests/test_krige.py b/tests/test_krige.py
index 3afee2b0a..75886c90b 100644
--- a/tests/test_krige.py
+++ b/tests/test_krige.py
@@ -156,8 +156,12 @@ def test_extdrift(self):
                 anis=[0.9, 0.8],
                 angles=[2, 1, 0.5],
             )
-            srf = SRF(model)
-            field = srf(grid)
+            # TODO should meshgrids be accepted?!
+            if dim > 1:
+                grid = [g.flatten() for g in grid]
+            # TODO the pre_points shit has to be applied here too
+            srf = SRF(model, grid)
+            field = srf()
             ext_drift.append(field)
             field = field.reshape(self.grid_shape[:dim])
             cond_drift.append(field[self.data_idx[:dim]])
diff --git a/tests/test_mesh.py b/tests/test_mesh.py
new file mode 100644
index 000000000..284097455
--- /dev/null
+++ b/tests/test_mesh.py
@@ -0,0 +1,225 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+"""
+This is the unittest of Mesh class.
+"""
+import unittest
+import numpy as np
+
+from gstools import Mesh
+from gstools.field.mesh import convert_points
+
+
+class TestMesh(unittest.TestCase):
+    def setUp(self):
+        self.x_grid = np.linspace(0.0, 12.0, 48)
+        self.y_grid = np.linspace(0.0, 10.0, 46)
+        self.z_grid = np.linspace(0.0, 10.0, 40)
+
+        self.f1_grid = self.x_grid
+        self.f2_grid = self.x_grid.reshape(
+            (len(self.x_grid), 1)
+        ) * self.y_grid.reshape((1, len(self.y_grid)))
+        self.f3_grid = (
+            self.x_grid.reshape((len(self.x_grid), 1, 1))
+            * self.y_grid.reshape((1, len(self.y_grid), 1))
+            * self.z_grid.reshape((1, 1, len(self.z_grid)))
+        )
+
+        self.rng = np.random.RandomState(123018)
+        self.x_tuple = self.rng.uniform(0.0, 10, 100)
+        self.y_tuple = self.rng.uniform(0.0, 10, 100)
+        self.z_tuple = self.rng.uniform(0.0, 10, 100)
+
+        self.f1_tuple = self.x_tuple
+        self.f2_tuple = self.x_tuple * self.y_tuple
+        self.f3_tuple = self.x_tuple * self.y_tuple * self.z_tuple
+
+        self.m1_grid = Mesh(1, (self.x_grid,), mesh_type="structured")
+        self.m2_grid = Mesh(
+            2, (self.x_grid, self.y_grid), mesh_type="structured"
+        )
+        self.m3_grid = Mesh(
+            3, (self.x_grid, self.y_grid, self.z_grid), mesh_type="structured"
+        )
+        self.m1_tuple = Mesh(1, (self.x_tuple,))
+        self.m2_tuple = Mesh(2, (self.x_tuple, self.y_tuple))
+        self.m3_tuple = Mesh(3, (self.x_tuple, self.y_tuple, self.z_tuple))
+
+    def test_convert_points(self):
+        pass
+        # TODO figure out what this is all about
+        #1d
+        #p1_target = np.array((1, 2, 3, 4)).reshape(4, 1)
+        #print(p1_target.shape)
+        #p1_is = convert_points(1, p1_target)
+        #np.testing.assert_array_equal(p1_is, p1_target)
+
+        #p2 = [1, 2, 3, 4, 5]
+        #p2_target = np.array(p2)
+        #p2_is = convert_points(1, p2)
+        #np.testing.assert_array_equal(p2_is, p2_target)
+
+        #p3 = [[1], [2], [3], [4], [5]]
+        #p3_target = np.array(p3).reshape(-1)
+        #p3_is = convert_points(1, p3)
+        #np.testing.assert_array_equal(p3_is, p3_target)
+
+        ##2d
+        #p4 = [[1, 2], [3, 4], [5, 6]]
+        #p4_target = np.array(p4)
+        #print(p4_target.shape)
+
+
+    def test_item_getter_setter(self):
+        self.m3_grid.add_field(256.0 * self.f3_grid, name="2nd")
+        self.m3_grid.add_field(512.0 * self.f3_grid, name="3rd")
+        self.assertEqual(
+            self.m3_grid["2nd"].all(), (256.0 * self.f3_grid).all()
+        )
+        self.assertEqual(
+            self.m3_grid["3rd"].all(), (512.0 * self.f3_grid).all()
+        )
+
+        self.m3_tuple["tmp_field"] = 2.0 * self.f3_tuple
+        self.assertEqual(
+            self.m3_tuple["tmp_field"].all(), (2.0 * self.f3_tuple).all()
+        )
+
+    def test_points_getter(self):
+        self.assertEqual(self.m1_grid.points, (self.x_grid,))
+        self.assertEqual(self.m1_tuple.points, (self.x_tuple,))
+        self.assertEqual(self.m2_grid.points, (self.x_grid, self.y_grid))
+        self.assertEqual(
+            self.m3_grid.points, (self.x_grid, self.y_grid, self.z_grid)
+        )
+        self.assertEqual(
+            self.m3_tuple.points, (self.x_tuple, self.y_tuple, self.z_tuple)
+        )
+
+    def test_default_value_type(self):
+        # value_type getter with no field set
+        self.assertEqual(self.m1_tuple.value_type, None)
+        self.assertEqual(self.m2_tuple.value_type, None)
+        self.assertEqual(self.m3_tuple.value_type, None)
+        self.assertEqual(self.m1_grid.value_type, None)
+        self.assertEqual(self.m2_grid.value_type, None)
+        self.assertEqual(self.m3_grid.value_type, None)
+
+    def test_field_data_setter(self):
+        # attribute creation by adding field_data
+        self.m2_tuple.set_field_data(3.14, "mean")
+        self.assertEqual(self.m2_tuple.field_data["mean"], 3.14)
+        self.assertEqual(self.m2_tuple.mean, 3.14)
+
+    def test_new_points(self):
+        # set new points (which causes reset)
+        x_tuple2 = self.rng.uniform(0.0, 10, 100)
+        y_tuple2 = self.rng.uniform(0.0, 10, 100)
+        self.m2_tuple.add_field(self.f2_tuple)
+
+        self.m2_tuple.points = (x_tuple2, y_tuple2)
+
+        self.assertEqual(self.m2_tuple.points, (x_tuple2, y_tuple2))
+
+        # previous field has to be deleted
+        self.assertEqual(self.m2_tuple.field, None)
+
+    def test_add_field(self):
+        # structured
+        self.m1_grid.add_field(self.f1_grid)
+        self.assertEqual(self.m1_grid.field.all(), self.f1_grid.all())
+        self.m2_grid.add_field(self.f2_grid)
+        self.assertEqual(self.m2_grid.field.all(), self.f2_grid.all())
+        self.m3_grid.add_field(self.f3_grid)
+        self.assertEqual(self.m3_grid.field.all(), self.f3_grid.all())
+
+        # unstructured
+        self.m1_tuple.add_field(self.f1_tuple)
+        self.assertEqual(self.m1_tuple.field.all(), self.f1_tuple.all())
+        self.m2_tuple.add_field(self.f2_tuple)
+        self.assertEqual(self.m2_tuple.field.all(), self.f2_tuple.all())
+        self.m3_tuple.add_field(self.f3_tuple)
+        self.assertEqual(self.m3_tuple.field.all(), self.f3_tuple.all())
+
+        # multiple fields
+        new_field = 10.0 * self.f1_grid
+        self.m1_grid.add_field(new_field, name="2nd")
+        # default field
+        self.assertEqual(self.m1_grid.field.all(), self.f1_grid.all())
+        self.assertEqual(self.m1_grid["2nd"].all(), new_field.all())
+        # overwrite default field
+        newer_field = 100.0 * self.f1_grid
+        self.m1_grid.add_field(newer_field, name="3rd", is_default_field=True)
+        self.assertEqual(self.m1_grid.field.all(), newer_field.all())
+
+    def test_default_field(self):
+        # first field added, should be default_field
+        f2_grid = 5.0 * self.f1_grid
+        self.m1_grid.add_field(self.f1_grid, name="test_field1")
+        self.m1_grid.add_field(f2_grid, name="test_field2")
+        self.assertEqual(self.m1_grid.default_field, "test_field1")
+        self.m1_grid.default_field = "test_field2"
+        self.assertEqual(self.m1_grid.default_field, "test_field2")
+        self.assertEqual(self.m1_grid.field[5], self.m1_grid["test_field2"][5])
+
+    def test_reset(self):
+        self.m3_tuple.set_field_data("TestLoc", "location")
+        self.m3_tuple.del_field_data()
+        self.assertEqual(self.m3_tuple.field_data["default_field"], "field")
+        self.assertEqual(self.m3_tuple.field_data["mesh_type"], "unstructured")
+        self.assertEqual(self.m3_tuple.default_field, "field")
+        self.assertEqual(self.m3_tuple.mesh_type, "unstructured")
+
+    def test_point_data_check(self):
+        self.assertRaises(ValueError, self.m1_tuple.add_field, self.f1_grid)
+        self.assertRaises(ValueError, self.m1_tuple.add_field, self.f2_grid)
+        self.assertRaises(ValueError, self.m1_tuple.add_field, self.f3_grid)
+        self.assertRaises(ValueError, self.m2_tuple.add_field, self.f1_grid)
+        self.assertRaises(ValueError, self.m2_tuple.add_field, self.f2_grid)
+        self.assertRaises(ValueError, self.m2_tuple.add_field, self.f3_grid)
+        self.assertRaises(ValueError, self.m3_tuple.add_field, self.f1_grid)
+        self.assertRaises(ValueError, self.m3_tuple.add_field, self.f2_grid)
+        self.assertRaises(ValueError, self.m3_tuple.add_field, self.f3_grid)
+
+        self.assertRaises(ValueError, self.m1_grid.add_field, self.f2_grid)
+        self.assertRaises(ValueError, self.m1_grid.add_field, self.f3_grid)
+        self.assertRaises(ValueError, self.m2_grid.add_field, self.f1_grid)
+        self.assertRaises(ValueError, self.m2_grid.add_field, self.f3_grid)
+        self.assertRaises(ValueError, self.m3_grid.add_field, self.f1_grid)
+        self.assertRaises(ValueError, self.m3_grid.add_field, self.f2_grid)
+        self.assertRaises(ValueError, self.m1_grid.add_field, self.f1_tuple)
+        self.assertRaises(ValueError, self.m1_grid.add_field, self.f2_tuple)
+        self.assertRaises(ValueError, self.m1_grid.add_field, self.f3_tuple)
+        self.assertRaises(ValueError, self.m2_grid.add_field, self.f1_tuple)
+        self.assertRaises(ValueError, self.m2_grid.add_field, self.f2_tuple)
+        self.assertRaises(ValueError, self.m2_grid.add_field, self.f3_tuple)
+        self.assertRaises(ValueError, self.m3_grid.add_field, self.f1_tuple)
+        self.assertRaises(ValueError, self.m3_grid.add_field, self.f2_tuple)
+        self.assertRaises(ValueError, self.m3_grid.add_field, self.f3_tuple)
+
+        x_tuple2 = self.rng.uniform(0.0, 10, 100)
+        y_tuple2 = self.rng.uniform(0.0, 10, 100)
+        f_tuple2 = np.vstack((x_tuple2, y_tuple2))
+        x_tuple3 = self.rng.uniform(0.0, 10, (3, 100))
+        y_tuple3 = self.rng.uniform(0.0, 10, (3, 100))
+        z_tuple3 = self.rng.uniform(0.0, 10, (3, 100))
+        f_tuple3 = np.vstack((x_tuple2, y_tuple2, z_tuple3))
+
+        m2_tuple = Mesh(2, (x_tuple2, y_tuple2))
+        m3_tuple = Mesh(3, (x_tuple3, y_tuple3, z_tuple3))
+
+        self.assertRaises(ValueError, m2_tuple.add_field, f_tuple3)
+        self.assertRaises(ValueError, m3_tuple.add_field, f_tuple2)
+
+        f_grid2 = np.zeros((2, len(self.x_grid), len(self.y_grid)))
+        f_grid3 = np.zeros(
+            (3, len(self.x_grid), len(self.y_grid), len(self.z_grid))
+        )
+
+        self.assertRaises(ValueError, self.m2_grid.add_field, f_grid3)
+        self.assertRaises(ValueError, self.m3_grid.add_field, f_grid2)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/tests/test_srf.py b/tests/test_srf.py
index 082a7ec38..53d4e2c02 100644
--- a/tests/test_srf.py
+++ b/tests/test_srf.py
@@ -32,8 +32,14 @@ def setUp(self):
 
     def test_shape_1d(self):
         self.cov_model.dim = 1
-        srf = SRF(self.cov_model, mean=self.mean, mode_no=self.mode_no)
-        field_str = srf([self.x_grid], seed=self.seed, mesh_type="structured")
+        srf = SRF(
+            self.cov_model,
+            (self.x_grid,),
+            mesh_type="unstructured",
+            mean=self.mean,
+            mode_no=self.mode_no,
+        )
+        field_str = srf(seed=self.seed)
         field_unstr = srf(
             [self.x_tuple], seed=self.seed, mesh_type="unstructured"
         )
@@ -73,6 +79,21 @@ def test_shape_3d(self):
         )
         self.assertEqual(field_unstr.shape, (len(self.x_tuple),))
 
+    def test_ens(self):
+        srf = SRF(
+            self.cov_model,
+            (self.x_tuple, self.y_tuple, self.z_tuple),
+            mean=self.mean,
+            mode_no=self.mode_no,
+        )
+        for i in range(3):
+            name = "ens_{}".format(i+1)
+            srf(name=name)
+            self.assertTrue(name in srf.point_data.keys())
+        self.assertEqual(srf.default_field, "ens_1")
+        # 3 ens + 1 raw_field
+        self.assertEqual(len(srf.point_data), 4)
+
     def test_anisotropy_2d(self):
         self.cov_model.dim = 2
         srf = SRF(self.cov_model, mean=self.mean, mode_no=self.mode_no)
@@ -125,8 +146,13 @@ def test_rotation_unstruct_2d(self):
         field_str = np.reshape(field, (y_len, x_len))
 
         self.cov_model.angles = -np.pi / 2.0
-        srf = SRF(self.cov_model, mean=self.mean, mode_no=self.mode_no)
-        field_rot = srf((x_u, y_u), seed=self.seed, mesh_type="unstructured")
+        srf = SRF(
+            self.cov_model,
+            (x_u, y_u),
+            mean=self.mean,
+            mode_no=self.mode_no
+        )
+        field_rot = srf(seed=self.seed)
         field_rot_str = np.reshape(field_rot, (y_len, x_len))
 
         self.assertAlmostEqual(field_str[0, 0], field_rot_str[-1, 0])