diff --git a/pyphare/pyphare/core/gridlayout.py b/pyphare/pyphare/core/gridlayout.py
index ec35311ca..9aaccb9df 100644
--- a/pyphare/pyphare/core/gridlayout.py
+++ b/pyphare/pyphare/core/gridlayout.py
@@ -44,6 +44,10 @@
         "Pyz": "primal",
         "Pzz": "primal",
         "tags": "dual",
+        "value": "primal",
+        "x": "primal",
+        "y": "primal",
+        "z": "primal",
     },
     "y": {
         "Bx": "dual",
@@ -79,6 +83,10 @@
         "Pyz": "primal",
         "Pzz": "primal",
         "tags": "dual",
+        "value": "primal",
+        "x": "primal",
+        "y": "primal",
+        "z": "primal",
     },
     "z": {
         "Bx": "dual",
@@ -114,6 +122,10 @@
         "Pyz": "primal",
         "Pzz": "primal",
         "tags": "dual",
+        "value": "primal",
+        "x": "primal",
+        "y": "primal",
+        "z": "primal",
     },
 }
 yee_centering_lower = {
diff --git a/pyphare/pyphare/core/operators.py b/pyphare/pyphare/core/operators.py
index 1bd81fb12..b8d3af86a 100644
--- a/pyphare/pyphare/core/operators.py
+++ b/pyphare/pyphare/core/operators.py
@@ -20,12 +20,12 @@ def _compute_dot_product(patch_datas, **kwargs):
 
 
 def _compute_sqrt(patch_datas, **kwargs):
-    ref_name = next(iter(patch_datas.keys()))
+    # ref_name = next(iter(patch_datas.keys())) TODO this is always "value"
 
     dset = np.sqrt(patch_datas["value"][:])
 
     return (
-        {"name": "value", "data": dset, "centering": patch_datas[ref_name].centerings},
+        {"name": "value", "data": dset, "centering": patch_datas["value"].centerings},
     )
 
 
@@ -139,3 +139,4 @@ def grad(hier, **kwargs):
     h = compute_hier_from(_compute_grad, hier, nb_ghosts=nb_ghosts)
 
     return VectorField(h)
+
diff --git a/pyphare/pyphare/core/ufuncs.py b/pyphare/pyphare/core/ufuncs.py
new file mode 100644
index 000000000..311a6dba2
--- /dev/null
+++ b/pyphare/pyphare/core/ufuncs.py
@@ -0,0 +1,170 @@
+import numpy as np
+
+from pyphare.pharesee.hierarchy import ScalarField, VectorField
+from pyphare.pharesee.hierarchy.hierarchy_utils import compute_hier_from
+
+from scipy.ndimage import gaussian_filter
+
+
+
+def _compute_gaussian_filter_on_scalarfield(patch_datas, **kwargs):
+    from scipy.ndimage import gaussian_filter
+
+    ndim = patch_datas["value"].box.ndim
+    nb_ghosts = kwargs["nb_ghosts"]
+    sigma = kwargs["sigma"]
+    ds = np.asarray(patch_datas["value"][:])
+
+    ds_ = np.full(list(ds.shape), np.nan)
+
+    gf_ = gaussian_filter(ds, sigma=sigma)
+    select = tuple([slice(nb_ghosts, -nb_ghosts) for _ in range(ndim)])
+    ds_[select] = np.asarray(gf_[select])
+
+    return (
+        {"name": "value", "data": ds_, "centering": patch_datas["value"].centerings},
+    )
+
+
+def _compute_gaussian_filter_on_vectorfield(patch_datas, **kwargs):
+    from scipy.ndimage import gaussian_filter
+
+    ref_name = next(iter(patch_datas.keys()))
+
+    ndim = patch_datas[ref_name].box.ndim
+    nb_ghosts = kwargs["nb_ghosts"]
+    sigma = kwargs["sigma"]
+    ds_x = np.asarray(patch_datas["x"][:])
+    ds_y = np.asarray(patch_datas["y"][:])
+    ds_z = np.asarray(patch_datas["z"][:])
+
+    dsx_ = np.full(list(ds_x.shape), np.nan)
+    dsy_ = np.full(list(ds_y.shape), np.nan)
+    dsz_ = np.full(list(ds_z.shape), np.nan)
+
+    gfx_ = gaussian_filter(ds_x, sigma=sigma)
+    gfy_ = gaussian_filter(ds_y, sigma=sigma)
+    gfz_ = gaussian_filter(ds_z, sigma=sigma)
+
+    select = tuple([slice(nb_ghosts, -nb_ghosts) for _ in range(ndim)])
+
+    dsx_[select] = np.asarray(gfx_[select])
+    dsy_[select] = np.asarray(gfy_[select])
+    dsz_[select] = np.asarray(gfz_[select])
+
+    return (
+        {"name": "x", "data": dsx_, "centering": patch_datas["x"].centerings},
+        {"name": "y", "data": dsy_, "centering": patch_datas["y"].centerings},
+        {"name": "z", "data": dsz_, "centering": patch_datas["z"].centerings},
+    )
+
+
+# def gFilt(hier, **kwargs):
+#     sigma = kwargs.get("sigma", 2)
+# 
+#     # time0 = list(hier.times())[0]
+#     # level0 = 0
+#     # p0 = 0
+#     # pd0 = hier.levels(time0)[level0].patches[p0].patch_datas
+#     # key0 = list(pd0.keys())[0]
+#     # nb_ghosts = np.max(pd0[key0].ghosts_nbr)
+# 
+#     nb_ghosts = np.max(list(hier.level(0).patches[0].patch_datas.values())[0].ghosts_nbr)
+# 
+#     if nb_ghosts < sigma :
+#         print("nb_ghosts ({0}) < sigma ({1}) : your gaussian filter might be dirty".format(nb_ghosts, sigma))
+# 
+#     if hier.ndim == 1:
+#         if isinstance(hier, ScalarField) :
+#             h = compute_hier_from(_compute_gaussian_filter_on_scalarfield, hier, nb_ghosts=nb_ghosts, sigma=sigma)
+#             return ScalarField(h)
+#         elif isinstance(hier, VectorField) :
+#             h = compute_hier_from(_compute_gaussian_filter_on_vectorfield, hier, nb_ghosts=nb_ghosts, sigma=sigma)
+#             return VectorField(h)
+#         else:
+#             return NotImplemented
+#     else:
+#         return NotImplemented
+
+
+def gF(hier, **kwargs):
+    sigma = kwargs.get("sigma", 4)
+    if sigma == 1:
+        raise ValueError("sigma value has to be at least 2")
+    h_ = hier.__deepcopy__(memo={})
+    ndim = hier.ndim
+    n_pad = 4*sigma+1
+    # The gaussian filter is calculated on the box extended by
+    # n_pad. Hence, the number of points is large enough so that the value
+    # at the last point of the real box is as equal as possible to the one
+    # at the first point of the next box...
+
+    for time in h_.times():
+        interp_ = hier.interpol(time)
+        for lvl in h_.levels(time).values():
+            for patch in lvl.patches:
+                names = list(patch.patch_datas.keys())
+                box = patch.box
+
+                for name in names:
+                    pdata = patch.patch_datas[name]
+                    nb_ghosts = pdata.ghosts_nbr
+                    if not n_pad > nb_ghosts:
+                        raise ValueError('sigma value is too small')
+
+                    r_ = []
+                    for i in range(ndim):
+                        s_ = np.arange(box.lower[i]-n_pad, box.upper[i]+2+n_pad)*pdata.dl[i]
+                        r_.append(s_)
+
+                    func, _ = interp_[name]
+
+                    if ndim == 1:
+                        data = func(r_[0])
+                    elif ndim == 2:
+                        data = func(r_[0], r_[1])
+                    elif ndim == 3:
+                        data = func(r_[0], r_[1], r_[2])
+                    else:
+                        raise ValueError('unvalid dimension')
+
+                    gf_ = gaussian_filter(data, sigma=sigma)
+                    select = tuple([slice(n_pad-nb_ghosts[i], -n_pad+nb_ghosts[i]) for i in range(ndim)])
+                    pdata.dataset = np.asarray(gf_[select])
+
+    return h_
+
+
+def peakIds(hier, **kwargs):
+    from scipy.signal import find_peaks
+
+    times = list(hier.times())
+    if len(times) == 1:
+        time = times[0]
+    else:
+        raise ValueError('multiple time is not yet implemented')
+
+    # pks_ = np.array([])
+    # hgs_ = np.array([])
+
+    names_ = kwargs.pop("names", None)
+    if names_ is None:
+        names_ = list(hier.levels(time)[0].patches[0].patch_datas.keys())
+
+    # ph_ = kwargs.get('peak_heights', None)
+    # if ph_ is None:
+    #     raise ValueError("the kwarg 'peak_heights' is mandatory for now...")
+
+    for lvl in hier.levels(time).values():
+        for patch in lvl.patches:
+            for name in names_:
+                pdata = patch.patch_datas[name]
+                ds = np.asarray(pdata.dataset)
+                # pks = find_peaks(ds, **kwargs)
+                # for pk, hg in zip(pks[0], pks[1]['peak_heights']):
+                #     pks_ = np.append(pks_, np.add(np.multiply(pk, patch.dl), patch.origin))
+                #     hgs_ = np.append(hgs_, hg)
+                out = find_peaks(ds, **kwargs)
+
+    # return pks_, hgs_
+    return zip(*out)
diff --git a/pyphare/pyphare/pharesee/hierarchy/hierarchy.py b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py
index 724eb6409..9b45554a2 100644
--- a/pyphare/pyphare/pharesee/hierarchy/hierarchy.py
+++ b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py
@@ -438,8 +438,9 @@ def plot1d(self, **kwargs):
                 label = "L{level}P{patch}".format(level=lvl_nbr, patch=ip)
                 marker = kwargs.get("marker", "")
                 ls = kwargs.get("ls", "--")
+                lw = kwargs.get("lw", 1)
                 color = kwargs.get("color", "k")
-                ax.plot(x, val, label=label, marker=marker, ls=ls, color=color)
+                ax.plot(x, val, label=label, marker=marker, ls=ls, lw=lw, color=color)
 
         ax.set_title(kwargs.get("title", ""))
         ax.set_xlabel(kwargs.get("xlabel", "x"))
@@ -614,6 +615,128 @@ def dist_plot(self, **kwargs):
 
         return final, dp(final, **kwargs)
 
+    def interpol(self, time, interp="nearest"):
+        from pyphare.pharesee.hierarchy.hierarchy_utils import flat_finest_field
+        from pyphare.pharesee.run.utils import build_interpolator
+
+        nbrGhosts = list(self.level(0, time).patches[0].patch_datas.values())[0].ghosts_nbr
+
+        interp_ = {}
+        for qty in self.quantities():
+            box = self.level(0, time).patches[0].patch_datas[qty].box
+            dl = self.level(0, time).patches[0].patch_datas[qty].dl
+            data, coords = flat_finest_field(self, qty, time=time)
+            interp_[qty] = build_interpolator(data, coords, interp, box, dl, qty, nbrGhosts)
+        return interp_
+
+    def __array_ufunc__(self, ufunc, method, *inputs, **kwargs):
+        print(f"__array_function__ of PatchHierarchy called for {ufunc.__name__}")
+        if method != "__call__":
+            return NotImplemented
+
+        final = []
+
+        times = inputs[0].times()
+        for x in inputs:
+            assert(times == x.times())
+            if not isinstance(x, PatchHierarchy):
+                raise TypeError("this arg should be a PatchHierarchy")
+
+
+        # for  d in inputs[0].patch_levels:
+        #     print(type(d), d)
+
+
+        ils = [key for dict_ in inputs[0].patch_levels for key in dict_]
+        h_type = type(inputs[0])
+
+        all_ = []
+        for i in range(len(times)):
+            pls = []
+            for x in inputs:
+                pls_ = []
+                for plvl in x.patch_levels[i].values():
+                    pls_.append(plvl)
+                pls.append(pls_)
+
+            out = [getattr(ufunc, method)(*pl, **kwargs)  for pl in zip(*pls)]
+
+            # out est une liste de liste de patchlevel : indice  sur le levels (pour 1 temps donne)
+            # il faut les remettre dans un dict avec ilvl
+
+            final = {}
+            for il, pl in zip(ils, out):
+                final[il] = pl
+
+            all_.append(final)
+
+        h_ = PatchHierarchy(all_,
+                           domain_box=self.domain_box,
+                           refinement_box=self.refinement_ratio,
+                           times=self.times(),
+                           data_files=self.data_files,
+                           selection_box=self.selection_box)
+
+        from .scalarfield import ScalarField
+        from .vectorfield import VectorField
+
+        if h_type is ScalarField:
+            return ScalarField(h_)
+        elif h_type is VectorField:
+            return VectorField(h_)
+
+    def __array_function__(self, func, types, args, kwargs):
+        print(f"__array_function__ of Patch {func.__name__} called for {[getattr(a, 'name', a) for a in args]}")
+
+        final = []
+        others = []
+
+        times = args[0].times()
+        for x in args:
+            if isinstance(x, PatchHierarchy):
+                assert(times == x.times())
+            else:
+                others.append(x)
+
+        ils = [key for dict_ in args[0].patch_levels for key in dict_]
+        h_type = type(args[0])
+
+
+        all_ = []
+        for i in range(len(times)):
+            pls = []
+            for x in args:
+                if isinstance(x, PatchHierarchy):
+                    pls_ = []
+                    for plvl in x.patch_levels[i].values():
+                        pls_.append(plvl)
+                    pls.append(pls_)
+
+            out = []
+            for pl in zip(*pls):
+                out.append(func(*pl, *others, **kwargs))
+
+            final = {}
+            for il, pl in zip(ils, out):
+                final[il] = pl
+
+            all_.append(final)
+
+        h_ = PatchHierarchy(all_,
+                           domain_box=self.domain_box,
+                           refinement_box=self.refinement_ratio,
+                           times=self.times(),
+                           data_files=self.data_files,
+                           selection_box=self.selection_box)
+
+        from .scalarfield import ScalarField
+        from .vectorfield import VectorField
+
+        if h_type is ScalarField:
+            return ScalarField(h_)
+        elif h_type is VectorField:
+            return VectorField(h_)
+
 
 def finest_part_data(hierarchy, time=None):
     """
diff --git a/pyphare/pyphare/pharesee/hierarchy/patch.py b/pyphare/pyphare/pharesee/hierarchy/patch.py
index 0ac18f3d8..5b857b6ca 100644
--- a/pyphare/pyphare/pharesee/hierarchy/patch.py
+++ b/pyphare/pyphare/pharesee/hierarchy/patch.py
@@ -1,4 +1,5 @@
 #
+from .patchdata import FieldData
 
 
 class Patch:
@@ -67,3 +68,55 @@ def __call__(self, qty, **kwargs):
             return pd.dataset[idx + nbrGhosts, nbrGhosts:-nbrGhosts]
         elif idim == 1:
             return pd.dataset[nbrGhosts:-nbrGhosts, idx + nbrGhosts]
+
+    def __array_ufunc__(self, ufunc, method, *inputs, **kwargs):
+        print(f"__array_ufunct__ of Patch called for {ufunc.__name__}")
+        if method != "__call__":
+            return NotImplemented
+
+        pds = []
+        for x in inputs:
+            if isinstance(x, Patch):
+                keys_ = []
+                pds_ = []
+                for k, p in x.patch_datas.items():
+                    keys_.append(k)
+                    pds_.append(p)
+                pds.append(pds_)
+            else:
+                raise TypeError("this arg should be a Patch")
+
+        out = [getattr(ufunc, method)(*pd, **kwargs) for pd in zip(*pds)]
+
+        final = {}
+        for k, pd in zip(keys_, out):  # TODO hmmmm, the output patch will keep the keys of the last patch in inputs
+            final[k] = pd
+
+        return Patch(final, patch_id=self.id, layout=self.layout, attrs=self.attrs)
+
+    def __array_function__(self, func, types, args, kwargs):
+        print(f"__array_function__ of Patch {func.__name__} called for {[getattr(a, 'name', a) for a in args]}")
+
+        pds = []
+        others = []
+        for x in args:
+            if isinstance(x, Patch):
+                keys_ = []
+                pds_ = []
+                for k, p in x.patch_datas.items():
+                    keys_.append(k)
+                    pds_.append(p)
+                pds.append(pds_)
+            else:
+                others.append(x)
+
+        out = []
+        for pd in zip(*pds):
+            out.append(func(*pd, *others, **kwargs))
+
+        final = {}
+        for k, pd in zip(keys_, out):
+            final[k] = pd
+
+        return Patch(final, patch_id=self.id, layout=self.layout, attrs=self.attrs)
+
diff --git a/pyphare/pyphare/pharesee/hierarchy/patchdata.py b/pyphare/pyphare/pharesee/hierarchy/patchdata.py
index e0727b9a2..b69b4faed 100644
--- a/pyphare/pyphare/pharesee/hierarchy/patchdata.py
+++ b/pyphare/pyphare/pharesee/hierarchy/patchdata.py
@@ -191,6 +191,32 @@ def grid():
             return tuple(g[select] for g in mesh)
         return mesh
 
+    def __array__(self, dtype=None):
+        # numpy array protocol
+        return np.asarray(self.dataset, dtype=dtype)
+
+    def __array_ufunc__(self, ufunc, method, *inputs, **kwargs):
+        print(f"__array_ufunc__ of FieldData called for {ufunc.__name__}")
+        if method != "__call__":
+            raise NotImplementedError
+
+        in_ = [i.dataset if isinstance(i, FieldData) else i for i in inputs]
+        out_ = getattr(ufunc, method)(*in_, **kwargs)
+
+        if isinstance(out_, np.ndarray):
+            return FieldData(self.layout, ufunc.__name__+"@"+self.field_name, out_, centering=self.centerings)
+
+    def __array_function__(self, func, types, args, kwargs):
+        print(f"__array_function__ of FieldData {func.__name__} called with {[getattr(a, 'name', a) for a in args]}")
+
+        in_ = [a.dataset if isinstance(a, FieldData) else a for a in args]
+        out_ = func(*in_, **kwargs)
+
+        if isinstance(out_, np.ndarray):
+            return FieldData(self.layout, func.__name__+"@"+self.field_name, out_, centering=self.centerings)
+        else:
+            return out_
+
 
 class ParticleData(PatchData):
     """
@@ -236,3 +262,12 @@ def compare(self, that, *args, **kwargs):
 
     def __eq__(self, that):
         return self.compare(that)
+
+    def __str__(self):
+        return "ParticleData: (pop_name={}, #_of_part={})".format(
+            self.pop_name, self.dataset.shape[0]  # TODO need to be tested
+        )
+
+    def __repr__(self):
+        return self.__str__()
+
diff --git a/pyphare/pyphare/pharesee/hierarchy/patchlevel.py b/pyphare/pyphare/pharesee/hierarchy/patchlevel.py
index 72fee0d5d..ded663f93 100644
--- a/pyphare/pyphare/pharesee/hierarchy/patchlevel.py
+++ b/pyphare/pyphare/pharesee/hierarchy/patchlevel.py
@@ -13,3 +13,42 @@ def level_range(self):
         return min([patch.patch_datas[name].x.min() for patch in self.patches]), max(
             [patch.patch_datas[name].x.max() for patch in self.patches]
         )
+
+    def __array_ufunc__(self, ufunc, method, *inputs, **kwargs):
+        print(f"__array_function__ of PatchLevel called for {ufunc.__name__}")
+        if method != "__call__":
+            return NotImplemented
+
+        ps = []
+        for x in inputs:
+            if isinstance(x, PatchLevel):
+                ps_ = []
+                for p in x.patches:
+                    ps_.append(p)
+                ps.append(ps_)
+            else:
+                raise TypeError("this arg should be a PatchLevel")
+
+        out = [getattr(ufunc, method)(*p, **kwargs)  for p in zip(*ps)]
+
+        return PatchLevel(self.level_number, out)
+
+    def __array_function__(self, func, types, args, kwargs):
+        print(f"__array_function__ of Patch {func.__name__} called for {[getattr(a, 'name', a) for a in args]}")
+
+        ps = []
+        others = []
+        for x in args:
+            if isinstance(x, PatchLevel):
+                ps_ = []
+                for p in x.patches:
+                    ps_.append(p)
+                ps.append(ps_)
+            else:
+                others.append(x)
+
+        out = []
+        for p in zip(*ps):
+            out.append(func(*p, *others, **kwargs))
+
+        return PatchLevel(self.level_number, out)
diff --git a/pyphare/pyphare/pharesee/hierarchy/scalarfield.py b/pyphare/pyphare/pharesee/hierarchy/scalarfield.py
index f8691098d..204e32ae8 100644
--- a/pyphare/pyphare/pharesee/hierarchy/scalarfield.py
+++ b/pyphare/pyphare/pharesee/hierarchy/scalarfield.py
@@ -210,3 +210,4 @@ def __neg__(self):
         names_self = self.quantities()
         h = compute_hier_from(_compute_neg, self, new_names=names_self)
         return ScalarField(h)
+
diff --git a/pyphare/pyphare/pharesee/run/utils.py b/pyphare/pyphare/pharesee/run/utils.py
index d6ffac24c..3cc588432 100644
--- a/pyphare/pyphare/pharesee/run/utils.py
+++ b/pyphare/pyphare/pharesee/run/utils.py
@@ -488,3 +488,50 @@ def make_interpolator(data, coords, interp, domain, dl, qty, nbrGhosts):
         raise ValueError("make_interpolator is not yet 3d")
 
     return interpolator, finest_coords
+
+
+def build_interpolator(data, coords, interp, box, dl, qty, nbrGhosts):
+    """
+    :param data: the values of the data that will be used for making
+    the interpolator, defined on coords
+    :param coords: coordinates where the data are known. they
+    can be define on an irregular grid (eg the finest)
+
+    finest_coords will be the structured coordinates defined on the
+    finest grid.
+    """
+    from pyphare.core.gridlayout import directions
+    from pyphare.core.gridlayout import yeeCoordsFor
+
+    dim = coords.ndim
+
+    nCells = box.upper-box.lower+2
+
+    fc = []
+    for i in range(dim):
+        s_ = yeeCoordsFor([0] * dim, nbrGhosts, dl, nCells, qty, directions[i])
+        fc.append(s_)
+
+    finest_coords = fc
+
+    if dim == 1:
+        from scipy.interpolate import interp1d
+
+        interpolator = interp1d(
+            coords, data, kind=interp, fill_value="extrapolate", assume_sorted=False
+        )
+
+    elif dim == 2:
+        from scipy.interpolate import LinearNDInterpolator, NearestNDInterpolator
+
+        if interp == "nearest":
+            interpolator = NearestNDInterpolator(coords, data)
+        elif interp == "bilinear":
+            interpolator = LinearNDInterpolator(coords, data)
+        else:
+            raise ValueError("the 'interp' kwarg can only be 'nearest' or 'bilinear'")
+
+    else:
+        raise ValueError("build_interpolator is not yet 3d")
+
+    return interpolator, finest_coords