can set a patchData w. an all_primal centering

pyphare/pyphare/core/operators.py

@@ -0,0 +1,106 @@
+import numpy as np
+from pyphare.pharesee.hierarchy import PatchLevel, Patch, FieldData
+from pyphare.pharesee.hierarchy import ScalarField, VectorField#, TensorField
+
+
+# TODO commencer par calculer z=hier*hier puis racine caree du dataset
+
+
+def sqrt(hier):
+
+    assert isinstance(hier, ScalarField)
+    domain_box = hier.domain_box
+
+    for time in list(hier.time_hier.keys()):
+        new_patch_level = {}
+
+        for ilvl, lvl in hier.levels(time).items():
+            new_patches = {}
+
+            for patch in lvl.patches:
+                layout = patch.layout
+
+                if ilvl not in new_patches:
+                    new_patches[ilvl] = []
+
+                new_pd = {}
+
+                dset = np.sqrt(np.asarray(patch.patch_datas['scalar'].dataset))
+                new_pd[ScalarField.name] = FieldData(layout, 'scalar', dset, centering=patch.patch_datas['scalar'].centerings)
+                new_patches[ilvl].append(Patch(new_pd, patch.id))
+
+            new_patch_level[ilvl] = PatchLevel(ilvl, new_patches[ilvl])
+
+    # return hier
+    return ScalarField(new_patch_level, domain_box, time=time) # TODO should this hierarchy be updated with all the times ?
+
+
+
+
+def modulus(hier):
+
+    assert isinstance(hier, VectorField)
+
+    h = hier*hier
+    return sqrt(h)
+
+
+
+
+
+def sqrt_old(hierarchy):
+
+    """
+    returns a Hierarchy of the same type of the input (which can be
+    ScalarField, VectorField or Tensor2Field) containing the square root
+    of each dataset
+
+    """
+
+    patch_levels = hierarchy.patch_levels
+    domain_box = hierarchy.domain_box
+
+    num_of_components = len(hierarchy.levels()[0].patches[0].patch_datas.keys())
+
+    if num_of_components == 1:
+        names = ['value']
+    elif num_of_components == 3:
+        names = ['x', 'y', 'z']
+    elif num_of_components == 6:
+        names = ['xx', 'xy', 'xz', 'yy', 'yz', 'zz']
+    elif num_of_components == 9:
+        names = ['xx', 'xy', 'xz', 'yx', 'yy', 'yz', 'zx', 'zy', 'zz']
+
+
+    for time in list(hierarchy.time_hier.keys()):
+        new_patch_level = {}
+
+        for ilvl, lvl in hierarchy.levels(time).items():
+            new_patches = {}
+
+            for patch in lvl.patches:
+                new_pd = {}
+                layout = patch.layout
+
+                num_of_components = len(patch.patch_datas.keys())
+
+                for (name, pd_name) in zip(names, patch.patch_datas.keys()):
+                    dset = np.sqrt(np.asarray(patch.patch_datas[pd_name].dataset))
+
+                    pd = FieldData(layout, name, dset, centering=patch.patch_datas[pd_name].centerings)
+                    new_pd[name] = pd
+
+                if ilvl not in new_patches:
+                    new_patches[ilvl] = []
+
+                new_patches[ilvl].append(Patch(new_pd, patch.id))
+
+            new_patch_level[ilvl] = PatchLevel(ilvl, new_patches[ilvl])
+
+    if num_of_components == 1:
+        return ScalarField(new_patch_level, domain_box, time=time)
+    if num_of_components == 3:
+        return VectorField(new_patch_level, domain_box, time=time)
+    if num_of_components == 6 or num_of_components == 9:
+        return Tensor2Field(new_patch_level, domain_box, time=time)
+

pyphare/pyphare/pharesee/hierarchy.py

@@ -563,6 +563,18 @@ def finest_part_data(hierarchy, time=None):
 class PatchHierarchy(object):
     """is a collection of patch levels"""
 
+    def flat_name(self, qty):
+        if qty in ['rho', 'tags']:
+            return 'scalar'
+        elif qty in ['Bx', 'Ex', 'Fx', 'Vx']:
+            return 'x'
+        elif qty in ['By', 'Ey', 'Fy', 'Vy']:
+            return 'y'
+        elif qty in ['Bz', 'Ez', 'Fz', 'Vz']:
+            return 'z'
+        else:
+            raise ValueError("{} is not a valid quantity".format(qty))
+
     def __init__(
         self,
         patch_levels,
@@ -1070,6 +1082,124 @@ def dist_plot(self, **kwargs):
         return final, dp(final, **kwargs)
 
 
+
+
+class ScalarField(PatchHierarchy):
+
+    name = 'scalar'
+
+    def __init__(self, patch_levels, domain_box, **kwargs):
+        refinement_ratio = kwargs.get("refinement_ratio", 2)
+        time = kwargs.get("time", .0)
+        data_files = kwargs.get("data_files", None)
+
+        self.name = ['scalar'] # TODO should we use the static name ?
+
+        super().__init__(patch_levels, domain_box, refinement_ratio, time, data_files)
+
+
+
+
+class VectorField(PatchHierarchy):
+
+    def __init__(self, patch_levels, domain_box, **kwargs):
+        refinement_ratio = kwargs.get("refinement_ratio", 2)
+        time = kwargs.get("time", .0)
+        data_files = kwargs.get("data_files", None)
+
+        self.names = ['x', 'y', 'z']
+
+        super().__init__(patch_levels, domain_box, refinement_ratio, time, data_files)
+
+
+    def __mul__(self, operand):
+
+        # scalar product between 2 vectors
+        if isinstance(operand, VectorField):
+
+            times_self = list(self.time_hier.keys())
+            times_operand = list(operand.time_hier.keys())
+
+            patch_levels = self.patch_levels
+            domain_box = self.domain_box
+
+            # TODO verifier qu'on a le meme centering et la meme size pur les dataset !
+            assert(all(l == r for l, r in zip(times_self, times_operand)))
+            for (time_self, time_operand) in zip(list(self.time_hier.keys()), list(operand.time_hier.keys())):
+                new_patch_level = {}
+
+                assert(all (l == r for l, r in zip(self.levels(time_self).keys(), operand.levels(time_operand).keys())))
+                for (ilvl_self, lvl_self), (ilvl_operand, lvl_operand) in zip(self.levels(time_self).items(), operand.levels(time_operand).items()):
+                    new_patches = {}
+
+                    assert(all (l.id == r.id for l, r in zip(lvl_self.patches, lvl_operand.patches)))
+                    for (patch_self, patch_operand) in zip(lvl_self.patches, lvl_operand.patches):
+                        layout = patch_self.layout
+
+                        if ilvl_self not in new_patches:
+                            new_patches[ilvl_self] = []
+
+                        new_pd = {}
+                        dset = 0.
+
+                        for (name, name_self, name_operand) in zip(self.names, patch_self.patch_datas.keys(), patch_operand.patch_datas.keys()):
+                            pd_self = patch_self.patch_datas[name_self]
+                            pd_operand = patch_operand.patch_datas[name_operand]
+
+                            dset = dset+np.asarray(pd_self.dataset)*np.asarray(pd_operand.dataset)
+
+                        # TODO not very nice to keep the last centering from the previous for loop !
+                        pd = FieldData(layout, ScalarField.name, dset, centering=pd_operand.centerings)
+                        new_pd[ScalarField.name] = pd
+
+                        new_patches[ilvl_self].append(Patch(new_pd, patch_self.id))
+
+                    new_patch_level[ilvl_self] = PatchLevel(ilvl_self, new_patches[ilvl_self])
+
+            return ScalarField(new_patch_level, domain_box, time=time_self)
+
+        # multiplication by a scalar
+        elif isinstance(operand, (int, float)):
+
+            patch_levels = self.patch_levels
+            domain_box = self.domain_box
+
+            names = ['x', 'y', 'z']
+
+            for time in list(self.time_hier.keys()):
+                new_patch_level = {}
+
+                for ilvl, lvl in self.levels(time).items():
+                    new_patches = {}
+
+                    for patch in lvl.patches:
+                        new_pd = {}
+                        layout = patch.layout
+
+                        for (name, pd_name) in zip(names, patch.patch_datas.keys()):
+                            dset = np.asarray(patch.patch_datas[pd_name].dataset)*operand
+                            pd = FieldData(layout, name, dset, centering=patch.patch_datas[pd_name].centerings)
+                            new_pd[name] = pd
+
+                        if ilvl not in new_patches:
+                            new_patches[ilvl] = []
+
+                        new_patches[ilvl].append(Patch(new_pd, patch.id))
+
+                    new_patch_level[ilvl] = PatchLevel(ilvl, new_patches[ilvl])
+
+            return VectorField(new_patch_level, domain_box, time=time)
+
+        else:
+            raise ValueError("Unsupported type of operand : {}".format(type(operand)))
+
+
+    def __rmul__(self, operand):
+        return self.__mul__(operand)
+
+
+
+
 def amr_grid(hierarchy, time):
     """returns a non-uniform contiguous primal grid
     associated to the given hierarchy