From 269b6256fb89d2785ae27c831dc9e2c37528606b Mon Sep 17 00:00:00 2001
From: dboe <dboe@ipp.mpg.de>
Date: Fri, 20 Aug 2021 12:17:41 +0200
Subject: [PATCH] added PHYSICAL_CYLINDER, transform_field and lib.io

---
 tests/test_grid.py          |  12 ++
 tfields/bases/__init__.py   |  23 +++-
 tfields/bases/bases.py      |  51 +++++----
 tfields/bases/manifold_3.py | 214 ++++++++++++++++++++++--------------
 tfields/core.py             | 125 ++++++++++++++++++---
 tfields/lib/io.py           |  94 ++++++++++++++++
 tfields/points_3d.py        |   2 +-
 tfields/tensor_grid.py      |  16 ++-
 8 files changed, 415 insertions(+), 122 deletions(-)
 create mode 100644 tfields/lib/io.py

diff --git a/tests/test_grid.py b/tests/test_grid.py
index 6340f7f..f0e1acf 100644
--- a/tests/test_grid.py
+++ b/tests/test_grid.py
@@ -120,6 +120,18 @@ class TensorGrid_Test_Permutation1(TensorGrid_Test):
 
 
 class TensorGrid_Test_IO_Change(unittest.TestCase):
+    def test_copy_constructor(self):
+        tgt = TensorGrid_Test()
+        tgt.setUp()
+        tg_orig = tgt._inst
+
+        tg_copy = tfields.TensorGrid(tg_orig)
+        tgt.check_filled(tg_copy)
+
+        tg_copy_new_num = tfields.TensorGrid(tg_orig, num=[1, 1, 1])
+        self.assertListEqual(list(tg_copy_new_num.num), [1, 1, 1])
+        self.assertEqual(len(tg_copy_new_num), 1)
+
     def test_change_iter_order(self):
         t1 = TensorGrid_Test()
         t1.setUp()
diff --git a/tfields/bases/__init__.py b/tfields/bases/__init__.py
index c23cec7..33390e0 100644
--- a/tfields/bases/__init__.py
+++ b/tfields/bases/__init__.py
@@ -7,7 +7,24 @@ Mail:       daniel.boeckenhoff@ipp.mpg.de
 part of tfields library
 Tools for sympy coordinate transformation
 """
-from tfields.bases.bases import get_coord_system, get_coord_system_name, lambdified_trafo, transform  # NOQA
+from tfields.bases.bases import (
+    get_coord_system,
+    get_coord_system_name,
+    lambdified_trafo,
+    transform,
+)  # NOQA
 from tfields.bases import manifold_3  # NOQA
-from tfields.bases.manifold_3 import CARTESIAN, CYLINDER, SPHERICAL  # NOQA
-from tfields.bases.manifold_3 import cartesian, cylinder, spherical  # NOQA
+from tfields.bases.manifold_3 import (
+    CARTESIAN,
+    CYLINDER,
+    SPHERICAL,
+    NATURAL_CARTESIAN,
+    PHYSICAL_CYLINDER,
+)  # NOQA
+from tfields.bases.manifold_3 import (
+    cartesian,
+    cylinder,
+    spherical,
+    natural_cartesian,
+    physical_cylinder,
+)  # NOQA
diff --git a/tfields/bases/bases.py b/tfields/bases/bases.py
index 04d7f30..6bd1e58 100644
--- a/tfields/bases/bases.py
+++ b/tfields/bases/bases.py
@@ -22,15 +22,17 @@ def get_coord_system(base):
     Return:
         sympy.diffgeom.get_coord_system
     """
-    if isinstance(base, string_types) or (isinstance(base, np.ndarray)
-                                          and base.dtype.kind in {'U', 'S'}):
+    if isinstance(base, string_types) or (
+        isinstance(base, np.ndarray) and base.dtype.kind in {"U", "S"}
+    ):
         base = getattr(tfields.bases, str(base))
     if not isinstance(base, sympy.diffgeom.CoordSystem):
         bse_tpe = type(base)
         expctd_tpe = type(sympy.diffgeom.CoordSystem)
-        raise TypeError("Wrong type of coord_system base {bse_tpe}. "
-                        "Expected {expctd_tpe}"
-                        .format(**locals()))
+        raise TypeError(
+            "Wrong type of coord_system base {bse_tpe}. "
+            "Expected {expctd_tpe}".format(**locals())
+        )
     return base
 
 
@@ -42,7 +44,7 @@ def get_coord_system_name(base):
         str: name of base
     """
     if isinstance(base, sympy.diffgeom.CoordSystem):
-        base = getattr(base, 'name')
+        base = getattr(base, "name")
     # if not (isinstance(base, string_types) or base is None):
     #     baseType = type(base)
     #     raise ValueError("Coordinate system must be string_type."
@@ -76,14 +78,15 @@ def lambdified_trafo(base_old, base_new):
 
     """
     coords = tuple(base_old.coord_function(i) for i in range(base_old.dim))
-    f = sympy.lambdify(coords,
-                       base_old.coord_tuple_transform_to(base_new,
-                                                         list(coords)),
-                       modules='numpy')
+    f = sympy.lambdify(
+        coords,
+        base_old.coord_tuple_transform_to(base_new, list(coords)),
+        modules="numpy",
+    )
     return f
 
 
-def transform(array, base_old, base_new):
+def transform(array, base_old, base_new, **kwargs):
     """
     Transform the input array in place
     Args:
@@ -119,31 +122,35 @@ def transform(array, base_old, base_new):
     """
     base_old = get_coord_system(base_old)
     base_new = get_coord_system(base_new)
-    if base_new not in base_old.transforms:
-        for baseTmp in base_new.transforms:
-            if baseTmp in base_old.transforms:
-                transform(array, base_old, baseTmp)
-                transform(array, baseTmp, base_new)
-                return
-        raise ValueError("Trafo not found.")
 
     # very fast trafos in numpy only
     short_trafo = None
     try:
-        short_trafo = getattr(base_old, 'to_{base_new.name}'.format(**locals()))
+        short_trafo = getattr(base_old, f"to_{base_new.name}".format(**locals()))
     except AttributeError:
         pass
     if short_trafo:
-        short_trafo(array)
+        short_trafo(array, **kwargs)
         return
 
+    if base_new not in base_old.transforms:
+        for baseTmp in base_new.transforms:
+            if baseTmp in base_old.transforms:
+                transform(array, base_old, baseTmp, **kwargs)
+                transform(array, baseTmp, base_new, **kwargs)
+                return
+        raise ValueError(f"Transformation {base_old} -> {base_new} not found.")
+
     # trafo via lambdified sympy expressions
     trafo = tfields.bases.lambdified_trafo(base_old, base_new)
     with warnings.catch_warnings():
-        warnings.filterwarnings('ignore', message="invalid value encountered in double_scalars")
+        warnings.filterwarnings(
+            "ignore", message="invalid value encountered in double_scalars"
+        )
         array[:] = np.concatenate([trafo(*coords).T for coords in array])
 
 
-if __name__ == '__main__':  # pragma: no cover
+if __name__ == "__main__":  # pragma: no cover
     import doctest
+
     doctest.testmod()
diff --git a/tfields/bases/manifold_3.py b/tfields/bases/manifold_3.py
index 203a87f..f3d2544 100644
--- a/tfields/bases/manifold_3.py
+++ b/tfields/bases/manifold_3.py
@@ -4,80 +4,96 @@ import sympy.diffgeom
 import numpy as np
 import warnings
 
-CARTESIAN = 'cartesian'
-CYLINDER = 'cylinder'
-SPHERICAL = 'spherical'
-
-m = sympy.diffgeom.Manifold('M', 3)
-patch = sympy.diffgeom.Patch('P', m)
-
-# cartesian
-x, y, z = sympy.symbols('x, y, z')
-cartesian = sympy.diffgeom.CoordSystem(CARTESIAN, patch, ['x', 'y', 'z'])
-
-# cylinder
-
-R, Phi, Z = sympy.symbols('R, Phi, Z')
-cylinder = sympy.diffgeom.CoordSystem(CYLINDER, patch, ['R', 'Phi', 'Z'])
-cylinder.connect_to(cartesian,
-                    [R, Phi, Z],
-                    [R * sympy.cos(Phi),
-                     R * sympy.sin(Phi),
-                     Z],
-                    inverse=False,
-                    fill_in_gaps=False)
+CARTESIAN = "cartesian"
+CYLINDER = "cylinder"
+SPHERICAL = "spherical"
+
+m = sympy.diffgeom.Manifold("M", 3)
+patch = sympy.diffgeom.Patch("P", m)
+
+#############
+# cartesian #
+#############
+x, y, z = sympy.symbols("x, y, z")
+cartesian = sympy.diffgeom.CoordSystem(CARTESIAN, patch, ["x", "y", "z"])
+
+############
+# cylinder #
+############
+
+R, Phi, Z = sympy.symbols("R, Phi, Z")
+cylinder = sympy.diffgeom.CoordSystem(CYLINDER, patch, ["R", "Phi", "Z"])
+cylinder.connect_to(
+    cartesian,
+    [R, Phi, Z],
+    [R * sympy.cos(Phi), R * sympy.sin(Phi), Z],
+    inverse=False,
+    fill_in_gaps=False,
+)
 
 
 def cylinder_to_cartesian(array):
-    rPoints = np.copy(array[:, 0])
-    phiPoints = np.copy(array[:, 1])
-    array[:, 0] = rPoints * np.cos(phiPoints)
-    array[:, 1] = rPoints * np.sin(phiPoints)
+    r_points = np.copy(array[:, 0])
+    phi_points = np.copy(array[:, 1])
+    array[:, 0] = r_points * np.cos(phi_points)
+    array[:, 1] = r_points * np.sin(phi_points)
 
 
 cylinder.to_cartesian = cylinder_to_cartesian
 
 
-cartesian.connect_to(cylinder,
-                     [x, y, z],
-                     [sympy.sqrt(x**2 + y**2),
-                      sympy.Piecewise(
-                          (sympy.pi, ((x < 0) & sympy.Eq(y, 0))),
-                          (0., sympy.Eq(sympy.sqrt(x**2 + y**2), 0)),
-                          (sympy.sign(y) * sympy.acos(x / sympy.sqrt(x**2 + y**2)), True)),
-                      z],
-                     inverse=False,
-                     fill_in_gaps=False)
+cartesian.connect_to(
+    cylinder,
+    [x, y, z],
+    [
+        sympy.sqrt(x ** 2 + y ** 2),
+        sympy.Piecewise(
+            (sympy.pi, ((x < 0) & sympy.Eq(y, 0))),
+            (0.0, sympy.Eq(sympy.sqrt(x ** 2 + y ** 2), 0)),
+            (sympy.sign(y) * sympy.acos(x / sympy.sqrt(x ** 2 + y ** 2)), True),
+        ),
+        z,
+    ],
+    inverse=False,
+    fill_in_gaps=False,
+)
 
 
 def cartesian_to_cylinder(array):
     x_vals = np.copy(array[:, 0])
     y_vals = np.copy(array[:, 1])
     problem_phi_indices = np.where((x_vals < 0) & (y_vals == 0))
-    array[:, 0] = np.sqrt(x_vals**2 + y_vals**2)  # r
-    np.seterr(divide='ignore', invalid='ignore')
+    array[:, 0] = np.sqrt(x_vals ** 2 + y_vals ** 2)  # r
+    np.seterr(divide="ignore", invalid="ignore")
     # phi
     array[:, 1] = np.sign(y_vals) * np.arccos(x_vals / array[:, 0])
-    np.seterr(divide='warn', invalid='warn')
+    np.seterr(divide="warn", invalid="warn")
     array[:, 1][problem_phi_indices] = np.pi
-    tfields.lib.util.convert_nan(array, 0.)  # for cases like cartesian 0, 0, 1
+    tfields.lib.util.convert_nan(array, 0.0)  # for cases like cartesian 0, 0, 1
 
 
 cartesian.to_cylinder = cartesian_to_cylinder
 
-# spherical
-r, phi, theta = sympy.symbols('r, phi, theta')
-spherical = sympy.diffgeom.CoordSystem(SPHERICAL, patch, ['r', 'phi', 'theta'])
-spherical.connect_to(cartesian,
-                     [r, phi, theta],
-                     [r * sympy.cos(phi - np.pi) * sympy.sin(theta - np.pi / 2),
-                      r * sympy.sin(phi - np.pi) * sympy.sin(theta - np.pi / 2),
-                      r * sympy.cos(theta - np.pi / 2)],
-                     # [r * sympy.cos(phi) * sympy.sin(theta),
-                     #  r * sympy.sin(phi) * sympy.sin(theta),
-                     #  r * sympy.cos(theta)],
-                     inverse=False,
-                     fill_in_gaps=False)
+#############
+# spherical #
+#############
+
+r, phi, theta = sympy.symbols("r, phi, theta")
+spherical = sympy.diffgeom.CoordSystem(SPHERICAL, patch, ["r", "phi", "theta"])
+spherical.connect_to(
+    cartesian,
+    [r, phi, theta],
+    [
+        r * sympy.cos(phi - np.pi) * sympy.sin(theta - np.pi / 2),
+        r * sympy.sin(phi - np.pi) * sympy.sin(theta - np.pi / 2),
+        r * sympy.cos(theta - np.pi / 2),
+    ],
+    # [r * sympy.cos(phi) * sympy.sin(theta),
+    #  r * sympy.sin(phi) * sympy.sin(theta),
+    #  r * sympy.cos(theta)],
+    inverse=False,
+    fill_in_gaps=False,
+)
 
 
 def spherical_to_cartesian(array):
@@ -97,26 +113,31 @@ def spherical_to_cartesian(array):
 spherical.to_cartesian = spherical_to_cartesian
 
 
-cartesian.connect_to(spherical,
-                     [x, y, z],
-                     [sympy.sqrt(x**2 + y**2 + z**2),
-                      sympy.Piecewise(
-                          (sympy.pi, ((x < 0) & sympy.Eq(y, 0))),
-                          (0., sympy.Eq(sympy.sqrt(x**2 + y**2), 0)),
-                          (sympy.sign(y) * sympy.acos(x / sympy.sqrt(x**2 + y**2)), True)),
-                      sympy.Piecewise(
-                          (0., sympy.Eq(x**2 + y**2 + z**2, 0)),
-                          (sympy.atan2(z, sympy.sqrt(x**2 + y**2)), True)),
-                      # sympy.Piecewise(
-                      #     (0., (sympy.Eq(x, 0) & sympy.Eq(y, 0))),
-                      #     (sympy.atan(y / x), True)),
-                      # sympy.Piecewise(
-                      #     (0., sympy.Eq(x**2 + y**2 + z**2, 0)),
-                      #     # (0., (sympy.Eq(x, 0) & sympy.Eq(y, 0) & sympy.Eq(z, 0))),
-                      #     (sympy.acos(z / sympy.sqrt(x**2 + y**2 + z**2)), True))
-                      ],
-                     inverse=False,
-                     fill_in_gaps=False)
+cartesian.connect_to(
+    spherical,
+    [x, y, z],
+    [
+        sympy.sqrt(x ** 2 + y ** 2 + z ** 2),
+        sympy.Piecewise(
+            (sympy.pi, ((x < 0) & sympy.Eq(y, 0))),
+            (0.0, sympy.Eq(sympy.sqrt(x ** 2 + y ** 2), 0)),
+            (sympy.sign(y) * sympy.acos(x / sympy.sqrt(x ** 2 + y ** 2)), True),
+        ),
+        sympy.Piecewise(
+            (0.0, sympy.Eq(x ** 2 + y ** 2 + z ** 2, 0)),
+            (sympy.atan2(z, sympy.sqrt(x ** 2 + y ** 2)), True),
+        ),
+        # sympy.Piecewise(
+        #     (0., (sympy.Eq(x, 0) & sympy.Eq(y, 0))),
+        #     (sympy.atan(y / x), True)),
+        # sympy.Piecewise(
+        #     (0., sympy.Eq(x**2 + y**2 + z**2, 0)),
+        #     # (0., (sympy.Eq(x, 0) & sympy.Eq(y, 0) & sympy.Eq(z, 0))),
+        #     (sympy.acos(z / sympy.sqrt(x**2 + y**2 + z**2)), True))
+    ],
+    inverse=False,
+    fill_in_gaps=False,
+)
 
 
 def cartesian_to_spherical(array):
@@ -126,25 +147,25 @@ def cartesian_to_spherical(array):
     phi in [-pi, pi]
     theta element [0, pi]: elevation angle defined from - Z-axis up
     """
-    xy = array[:, 0]**2 + array[:, 1]**2
+    xy = array[:, 0] ** 2 + array[:, 1] ** 2
 
     # phi for phi between -pi, pi
     problemPhiIndices = np.where((array[:, 0] < 0) & (array[:, 1] == 0))
     with warnings.catch_warnings():
         # python2.7
-        warnings.filterwarnings('ignore',
-                                message="invalid value encountered in divide")
+        warnings.filterwarnings("ignore", message="invalid value encountered in divide")
         # python3.x
-        warnings.filterwarnings('ignore',
-                                message="invalid value encountered in true_divide")
+        warnings.filterwarnings(
+            "ignore", message="invalid value encountered in true_divide"
+        )
         array[:, 1] = np.sign(array[:, 1]) * np.arccos(array[:, 0] / np.sqrt(xy))
     array[:, 1][problemPhiIndices] = np.pi
 
-    tfields.lib.util.convert_nan(array, 0.)  # for cases like cartesian 0, 0, 1
+    tfields.lib.util.convert_nan(array, 0.0)  # for cases like cartesian 0, 0, 1
     # array[:,1] = np.arctan2(array[:,1], array[:,0])  # phi for phi between 0, 2pi
 
     # r
-    array[:, 0] = np.sqrt(xy + array[:, 2]**2)
+    array[:, 0] = np.sqrt(xy + array[:, 2] ** 2)
 
     # theta
     # for elevation angle defined from - Z-axis up:
@@ -154,3 +175,36 @@ def cartesian_to_spherical(array):
 
 
 cartesian.to_spherical = cartesian_to_spherical
+
+#################################################################################
+# TODO: the part below is a temporary hack. I do not know yet the maths properly.
+#       will have to learn that and generalize.
+#################################################################################
+
+PHYSICAL_CYLINDER = "physical_cylinder"
+NATURAL_CARTESIAN = "natural_cartesian"
+
+natural_cartesian = sympy.diffgeom.CoordSystem(
+    NATURAL_CARTESIAN, patch, ["x", "y", "z"]
+)
+physical_cylinder = sympy.diffgeom.CoordSystem(
+    PHYSICAL_CYLINDER, patch, ["R", "Phi", "Z"]
+)
+
+
+def physical_cylinder_to_natural_cartesian(vector, position=None):
+    """
+    We need to project from PHYSICAL cylinder coordinates to NATURAL cartesian coordinates
+    """
+    with position.tmp_transform(CYLINDER):
+        sin_phi = np.copy(np.sin(position[:, 1]))
+        cos_phi = np.copy(np.cos(position[:, 1]))
+
+    vec_r = np.copy(vector[:, 0])
+    vec_phi = np.copy(vector[:, 1])
+
+    vector[:, 0] = vec_r * cos_phi - vec_phi * sin_phi
+    vector[:, 1] = vec_phi * cos_phi - vec_r * sin_phi
+
+
+physical_cylinder.to_natural_cartesian = physical_cylinder_to_natural_cartesian
diff --git a/tfields/core.py b/tfields/core.py
index 560874e..f62f0cd 100644
--- a/tfields/core.py
+++ b/tfields/core.py
@@ -17,6 +17,7 @@ TODO:
 """
 # builtin
 import warnings
+import inspect
 from contextlib import contextmanager
 from collections import Counter
 from copy import deepcopy
@@ -287,7 +288,7 @@ class AbstractNdarray(np.ndarray, AbstractObject):
             defaults for the attributes in __slots__ will be None
             other values will be treaded as defaults to the corresponding
             arg at the same position in the __slots__ list.
-        __slot_dtype__ (List(dtypes)): for the conversion of the
+        __slot_dtypes__ (List(dtypes)): for the conversion of the
             args in __slots__ to numpy arrays. None values mean no
             conversion.
         __slot_setters__ (List(callable)): Because __slots__ and properties are
@@ -704,6 +705,89 @@ class Tensors(AbstractNdarray):  # pylint: disable=too-many-public-methods
         for index in range(len(self)):
             yield super(Tensors, self).__getitem__(index).view(Tensors)
 
+    @classmethod
+    def _load_txt(cls, path, set_header: bool = False, **load_kwargs):
+        """
+        Factory method
+        Given a path to a txt file, construct the object
+
+        Args:
+            path: see :meth:`Tensors.load`
+            set_header: if Treu, set the header list to the `name` attribute
+            **loadkwargs: see :meth:`np.loadtxt`
+        """
+        load_txt_keys = [
+            k
+            for k, v in inspect.signature(np.loadtxt).parameters.items()
+            if v.default is not inspect._empty  # pylint:disable=protected-access
+        ]
+        load_txt_kwargs = {}
+        for key in load_txt_keys:
+            if key in load_kwargs:
+                load_txt_kwargs[key] = load_kwargs.pop(key)
+
+        # use the same defaults as np.savetxt
+        load_txt_kwargs.setdefault("comments", "# ")  # default is "#"
+        comments = load_txt_kwargs.get("comments")
+        load_txt_kwargs.setdefault("delimiter", " ")
+
+        skiprows = 0
+        header = []
+        with open(rna.path.resolve(path)) as file_:
+            while True:
+                line = file_.readline()
+                if line.startswith(comments):
+                    header.append(line.lstrip(comments).rstrip("\n"))
+                else:
+                    file_.seek(0)
+                    break
+                skiprows += 1
+        load_txt_kwargs["skiprows"] = max(skiprows, load_txt_kwargs.pop("skiprows", 0))
+
+        arr = np.loadtxt(path, **load_txt_kwargs)
+        obj = cls(arr, **load_kwargs)
+        i = 0
+        for key, line in zip(cls._iter_slots(), header):
+            slot, rest = line.split(": ")
+            tpe, val_str = rest.split(" = ")
+            if tpe == "NoneType":
+                val = None
+            else:
+                tpe = __builtins__[tpe]
+                val = tpe(val_str)
+            setattr(obj, slot, val)
+            i += 1
+        header = header[i:]
+
+        if set_header:
+            obj.name = (
+                obj.name,
+                header,
+            )  # pylint:disable=attribute-defined-outside-init
+        return obj
+
+    def _save_txt(self, path, **kwargs):
+        """
+        Save as text file.
+
+        Args:
+            **kwargs passed to np.savetxt.
+        """
+        header = kwargs.get("header", [])
+        if isinstance(header, dict):
+            header = [
+                f"{key}: {type(value).__name__} = {value}"
+                for key, value in header.items()
+            ]
+        if isinstance(header, list):
+            # statictyping like attribute saving
+            header = [
+                f"{key}: {type(getattr(self, key)).__name__} = {getattr(self, key)}"
+                for key in self._iter_slots()
+            ] + header
+            kwargs["header"] = "\n".join(header)
+        np.savetxt(path, self, **kwargs)
+
     @classmethod
     def merged(cls, *objects, **kwargs):
         """
@@ -925,7 +1009,7 @@ class Tensors(AbstractNdarray):  # pylint: disable=too-many-public-methods
         """
         return dim(self)
 
-    def transform(self, coord_sys):
+    def transform(self, coord_sys, **kwargs):
         """
         Args:
             coord_sys (str)
@@ -1000,7 +1084,7 @@ class Tensors(AbstractNdarray):  # pylint: disable=too-many-public-methods
             # already correct
             return
 
-        tfields.bases.transform(self, self.coord_sys, coord_sys)
+        tfields.bases.transform(self, self.coord_sys, coord_sys, **kwargs)
         # self[:] = tfields.bases.transform(self, self.coord_sys, coord_sys)
         self.coord_sys = coord_sys  # pylint: disable=attribute-defined-outside-init
 
@@ -1695,6 +1779,8 @@ class Tensors(AbstractNdarray):  # pylint: disable=too-many-public-methods
 
     def plot(self, *args, **kwargs):
         """
+        Generic plotting method of Tensors.
+
         Forwarding to rna.plotting.plot_tensor
         """
         artist = rna.plotting.plot_tensor(
@@ -1967,6 +2053,15 @@ class TensorFields(Tensors):
         else:
             return inst
 
+    def transform_field(self, coord_sys, field_index=0):
+        """
+        Transform the field to the coordinate system of choice.
+
+        NOTE: This is not yet any way generic!!! Have a look at Einsteinpy and actual status of
+            sympy for further implementation
+        """
+        self.fields[field_index].transform(coord_sys, position=self)
+
     @property
     def names(self):
         """
@@ -2002,8 +2097,7 @@ class TensorFields(Tensors):
 
         Args:
             other (iterable)
-            optional:
-                see Tensors.equal
+            **kwargs: see Tensors.equal
         """
         if not issubclass(type(other), Tensors):
             return super(TensorFields, self).equal(other, **kwargs)
@@ -2032,11 +2126,11 @@ class TensorFields(Tensors):
 
     def plot(self, *args, **kwargs):
         """
-        Plotting the tensor field.
+        Generic plotting method of TensorFields.
 
         Args:
             field_index: index of the field to plot (as quiver by default)
-            normalize: If true, normalize the field vectors to show only the direction
+            normalize: If True, normalize the field vectors to show only the direction
             color: additional str argument 'norm' added. If color="norm", color with the norm.
         """
         field_index = kwargs.pop("field_index", None)
@@ -2053,13 +2147,14 @@ class TensorFields(Tensors):
             color = kwargs.get("color", None)
 
             field = self.fields[field_index].copy()
-            if self.dim == field.dim:
-                field.transform(self.coord_sys)
-            else:
-                logging.debug(
-                    "Careful: Plotting tensors with field of"
-                    "different dimension. No coord_sys check performed."
-                )
+            # this tranfomration is compmlicated. Transforming the field is not yet understood
+            # if self.dim == field.dim:
+            #     field.transform(self.coord_sys)
+            # else:
+            #     logging.debug(
+            #         "Careful: Plotting tensors with field of"
+            #         "different dimension. No coord_sys check performed."
+            #     )
 
             if color == "norm":
                 norm = field.norm()
@@ -2928,7 +3023,7 @@ class TensorMaps(TensorFields):
 
     def plot(self, *args, **kwargs):  # pragma: no cover
         """
-        Plot the tensor map.
+        Generic plotting method of TensorMaps.
         """
         scalars_demanded = (
             "color" not in kwargs
diff --git a/tfields/lib/io.py b/tfields/lib/io.py
new file mode 100644
index 0000000..5c07717
--- /dev/null
+++ b/tfields/lib/io.py
@@ -0,0 +1,94 @@
+import typing
+import ast
+import sys
+import importlib
+import base64
+import numpy as np
+
+
+def get_module_and_name(type_) -> typing.Tuple[str, str]:
+    """
+    Examples:
+        >>> import numpy as np
+
+        This function can be used to ban your type to file as a string and (with the get_type
+        method) get it back.
+        >>> [get_type(*get_module_and_name(type_)) for type_ in (int, np.ndarray, str)]
+        [<class 'int'>, <class 'numpy.ndarray'>, <class 'str'>]
+    """
+    return sys.modules[type_.__module__].__name__, type_.__qualname__
+
+
+def get_type(module, name):
+    """
+    Inverse to :fun:`get_module_and_name`
+    """
+    importlib.import_module(module)
+    return getattr(sys.modules[module], name)
+
+
+def numpy_to_str(arr):
+    """
+    Convert an array to string representation
+
+    Examples
+        >>> import numpy as np
+        >>> arr = np.array([[1,2,3], [1,4,5]])
+        >>> enc = numpy_to_str(arr)
+        >>> str_to_numpy(enc)
+        array([[1, 2, 3],
+               [1, 4, 5]])
+    """
+    arr = np.ascontiguousarray(arr)
+    str_rep = base64.binascii.b2a_base64(arr).decode("ascii").rstrip("\n")
+    shape = arr.shape
+    dtype = arr.dtype
+    str_ = f"{str_rep}::{shape}::{dtype}"
+    return str_
+
+
+def str_to_numpy(str_):
+    """
+    Convert back from numpy_to_str
+    """
+    str_rep, shape, dtype = str_.split("::")
+    dtype = getattr(np, dtype)
+    arr = np.frombuffer(
+        base64.binascii.a2b_base64(str_rep.encode("ascii")), dtype=dtype
+    )
+    arr = arr.reshape(ast.literal_eval(shape))
+    return arr
+
+
+def numpy_to_bytes(arr: np.array) -> bytearray:
+    """
+    Convert to bytest array
+
+    Examples:
+        >>> a = np.ones((23, 23), dtype = 'int')
+        >>> a_b = numpy_to_bytes(a)
+        >>> a1 = bytes_to_numpy(a_b)
+        >>> assert np.array_equal(a, a1) and a.shape == a1.shape and a.dtype == a1.dtype
+    """
+    arr_dtype = bytearray(str(arr.dtype), "utf-8")
+    arr_shape = bytearray(",".join([str(a) for a in arr.shape]), "utf-8")
+    sep = bytearray("|", "utf-8")
+    arr_bytes = arr.ravel().tobytes()
+    to_return = arr_dtype + sep + arr_shape + sep + arr_bytes
+    return to_return
+
+
+def bytes_to_numpy(serialized_arr: bytearray) -> np.array:
+    """
+    Convert back from numpy_to_bytes
+    """
+    sep = "|".encode("utf-8")
+    i_0 = serialized_arr.find(sep)
+    i_1 = serialized_arr.find(sep, i_0 + 1)
+    arr_dtype = serialized_arr[:i_0].decode("utf-8")
+    arr_shape = tuple(
+        int(a) for a in serialized_arr[i_0 + 1 : i_1].decode("utf-8").split(",")
+    )
+    arr_str = serialized_arr[i_1 + 1 :]
+    arr = np.frombuffer(arr_str, dtype=arr_dtype).reshape(arr_shape)
+    return arr
diff --git a/tfields/points_3d.py b/tfields/points_3d.py
index e4f1ed4..4002b32 100644
--- a/tfields/points_3d.py
+++ b/tfields/points_3d.py
@@ -6,8 +6,8 @@ Mail:       daniel.boeckenhoff@ipp.mpg.de
 
 basic threedimensional tensors
 """
-import tfields
 import numpy as np
+import tfields
 
 
 class Points3D(tfields.Tensors):
diff --git a/tfields/tensor_grid.py b/tfields/tensor_grid.py
index e0f8586..3b59eb9 100644
--- a/tfields/tensor_grid.py
+++ b/tfields/tensor_grid.py
@@ -23,10 +23,11 @@ class TensorGrid(TensorFields):
     __slot_setters__ = TensorFields.__slot_setters__ + [
         None,
         None,
+        None,
     ]
 
     def __new__(cls, tensors, *fields, **kwargs):
-        if issubclass(type(tensors), TensorGrid):
+        if isinstance(tensors, TensorGrid):
             default_base_vectors = tensors.base_vectors
             default_num = tensors.num
             default_iter_order = tensors.iter_order
@@ -73,6 +74,19 @@ class TensorGrid(TensorFields):
         else:
             iter_order = np.array(iter_order, dtype=int)
         obj.iter_order = iter_order
+
+        if isinstance(tensors, TensorGrid):
+            if (obj.num != tensors.num).all() or (
+                obj.is_empty() and not obj.base_vectors.equal(tensors.base_vectors)
+            ):
+                # copy constructor with shape change
+                return obj.empty(*obj.base_num_tuples(), iter_order=iter_order)
+            if (obj.iter_order != tensors.iter_order).all():
+                # iter_order was changed in copy constructor
+                obj.iter_order = (
+                    tensors.iter_order
+                )  # set to 'default_iter_order' and change later
+                obj.change_iter_order(iter_order)
         return obj
 
     def __getitem__(self, index):
-- 
GitLab