Update to 3D

.github/workflows/cmake_macos.yml

@@ -78,7 +78,7 @@ jobs:
               -DCMAKE_BUILD_TYPE=RelWithDebInfo \
               -DENABLE_SAMRAI_TESTS=OFF -DCMAKE_C_COMPILER_LAUNCHER=ccache \
               -DCMAKE_CXX_COMPILER_LAUNCHER=ccache -DlowResourceTests=ON \
-              -DCMAKE_CXX_FLAGS="-DPHARE_DIAG_DOUBLES=1 "
+              -DCMAKE_CXX_FLAGS="-DPHARE_DIAG_DOUBLES=1 -DPHARE_SIMULATORS=2 "
 
     - name: Build
       working-directory: ${{runner.workspace}}/build

pyphare/pyphare/core/box.py

@@ -62,6 +62,9 @@ def __eq__(self, other):
         )
 
     def __sub__(self, other):
+        if isinstance(other, (list, tuple)):
+            assert all([isinstance(item, Box) for item in other])
+            return remove_all(self, other)
         assert isinstance(other, Box)
         return remove(self, other)
 
@@ -179,5 +182,42 @@ def copy(arr, replace):
     return list(boxes.values())
 
 
+def remove_all(box, to_remove):
+    if len(to_remove) > 0:
+        remaining = box - to_remove[0]
+        for to_rm in to_remove[1:]:
+            tmp, remove = [], []
+            for i, rem in enumerate(remaining):
+                if rem * to_rm is not None:
+                    remove.append(i)
+                    tmp += rem - to_rm
+            for rm in reversed(remove):
+                del remaining[rm]
+            remaining += tmp
+        return remaining
+    return box
+
+
+
 def amr_to_local(box, ref_box):
     return Box(box.lower - ref_box.lower, box.upper - ref_box.lower)
+
+
+def select(data, box):
+    return data[tuple([slice(l, u + 1) for l,u in zip(box.lower, box.upper)])]
+
+class DataSelector:
+    """
+        can't assign return to function unless []
+        usage
+        DataSelector(data)[box] = val
+    """
+    def __init__(self, data):
+        self.data = data
+    def __getitem__(self, box_or_slice):
+        if isinstance(box_or_slice, Box):
+            return select(self.data, box_or_slice)
+        return self.data[box_or_slice]
+
+    def __setitem__(self, box_or_slice, val):
+        self.__getitem__(box_or_slice)[:] = val

pyphare/pyphare/core/gridlayout.py

@@ -366,6 +366,20 @@ def yeeCoords(self, knode, iStart, centering, direction, ds, origin, derivOrder)
 
         return x
 
+    def meshCoords(self, qty):
+        ndim = self.ndim
+        assert ndim > 0 and ndim < 4
+        x = self.yeeCoordsFor(qty, "x")
+        if ndim == 1:
+            return x
+        y = self.yeeCoordsFor(qty, "y")
+        if ndim == 2:
+            X,Y = np.meshgrid(x,y,indexing="ij")
+            return np.array([X.flatten(),Y.flatten()]).T.reshape((len(x), len(y), ndim))
+        z = self.yeeCoordsFor(qty, "z")
+        X ,Y, Z = np.meshgrid(x,y,z,indexing="ij")
+        return np.array([X.flatten(), Y.flatten(), Z.flatten()]).T.reshape((len(x), len(y), len(z), ndim))
+
     def yeeCoordsFor(self, qty, direction, withGhosts=True, **kwargs):
         """
         from a qty and a direction, returns a 1d array containing

pyphare/pyphare/pharesee/geometry.py

@@ -54,19 +54,24 @@ def domain_border_ghost_boxes(domain_box, patches):
     elif domain_box.ndim == 2:
         upper_x, upper_y = domain_box.upper
         return {
-            "bottom": Box(
-                (
-                    0,
-                    0,
-                ),
-                (upper_x, ghost_box_width),
-            ),
-            "top": Box((0, upper_y - ghost_box_width), (upper_x, upper_y)),
             "left": Box((0, 0), (ghost_box_width, upper_y)),
             "right": Box((upper_x - ghost_box_width, 0), (upper_x, upper_y)),
+            "bottom": Box((0, 0), (upper_x, ghost_box_width)),
+            "top": Box((0, upper_y - ghost_box_width), (upper_x, upper_y)),
         }
 
-    raise ValueError("Unhandeled dimension")
+    else:
+        upper_x, upper_y, upper_z = domain_box.upper
+        return {
+            "left": Box((0, 0, 0), (ghost_box_width, upper_y, upper_z)),
+            "right": Box(
+                (upper_x - ghost_box_width, 0, 0), (upper_x, upper_y, upper_z)
+            ),
+            "bottom": Box((0, 0, 0), (upper_x, ghost_box_width, upper_z)),
+            "top": Box((0, upper_y - ghost_box_width, 0), (upper_x, upper_y, upper_z)),
+            "front": Box((0, 0, 0), (upper_x, upper_y, ghost_box_width)),
+            "back": Box((0, 0, upper_z - ghost_box_width), (upper_x, upper_y, upper_z)),
+        }
 
 
 def touch_domain_border(box, domain_box, border):
@@ -79,21 +84,29 @@ def touch_domain_border(box, domain_box, border):
 
 
 def periodicity_shifts(domain_box):
+    shifts = {}
+
     if domain_box.ndim == 1:
         shape_x = domain_box.shape
-        return {
+        shifts = {
             "left": shape_x,
             "right": -shape_x,
         }
+        shifts.update({"leftright": [shifts["left"], shifts["right"]]})
 
     if domain_box.ndim == 2:
         shape_x, shape_y = domain_box.shape
+    if domain_box.ndim == 3:
+        shape_x, shape_y, shape_z = domain_box.shape
+
+    if domain_box.ndim > 1:
         shifts = {
             "left": [(shape_x, 0)],
             "right": [(-shape_x, 0)],
             "bottom": [(0, shape_y)],
             "top": [(0, -shape_y)],
         }
+
         shifts.update(
             {
                 "bottomleft": [*shifts["left"], *shifts["bottom"], (shape_x, shape_y)],
@@ -134,7 +147,7 @@ def periodicity_shifts(domain_box):
                     shifts["topleft"][-1],
                     shifts["topright"][-1],
                 ],
-                "bottomtopleftright": [  # one patch covers domain
+                "leftrightbottomtop": [  # one patch covers domain
                     *shifts["bottomleft"],
                     *shifts["topright"],
                     shifts["bottomright"][-1],
@@ -144,12 +157,122 @@ def periodicity_shifts(domain_box):
         )
 
     if domain_box.ndim == 3:
-        raise ValueError("Unhandeled dimension")
+        front = {
+            f"{k}front": [(v[0], v[1], shape_z) for v in l] for k, l in shifts.items()
+        }
+        back = {
+            f"{k}back": [([v[0], v[1], -shape_z]) for v in l] for k, l in shifts.items()
+        }
+
+        shifts = {k: [([v[0], v[1], 0]) for v in l] for k, l in shifts.items()}
+
+        shifts.update(front)
+        shifts.update(back)
+        shifts.update(
+            {
+                "back": [(0, 0, -shape_z)],
+                "front": [(0, 0, shape_z)],
+                "leftrightbottomtopfrontback": [
+                    *shifts["bottomleftfront"],
+                    *shifts["bottomrightback"],
+                    *shifts["topleftfront"],
+                    *shifts["toprightback"],
+                ],
+            }
+        )
+
+    assert len(list(shifts.keys())) == len(
+        set(["".join(sorted(k)) for k in list(shifts.keys())])
+    )
+
+    shifts = {"".join(sorted(k)): l for k, l in shifts.items()}
 
     return shifts
 
 
+def possible_periodic_shifts(box, domain_box):
+    from pyphare.core.box import shift
+
+    boxes = []
+    dim = domain_box.ndim
+
+    if dim == 1:
+        for ishift in (-1, 0, 1):
+            ioffset = ishift * domain_box.shape[0]
+            offset = ioffset
+
+            shifted = shift(box, offset)
+            boxes += [(offset, shifted)]
+    if dim == 2:
+        for ishift in (-1, 0, 1):
+            ioffset = ishift * domain_box.shape[0]
+            for jshift in (-1, 0, 1):
+                joffset = jshift * domain_box.shape[1]
+                offset = [ioffset, joffset]
+                shifted = shift(box, offset)
+                boxes += [(offset, shifted)]
+    if dim == 3:
+        for ishift in (-1, 0, 1):
+            ioffset = ishift * domain_box.shape[0]
+            for jshift in (-1, 0, 1):
+                joffset = jshift * domain_box.shape[1]
+                for kshift in (-1, 0, 1):
+                    koffset = kshift * domain_box.shape[2]
+
+                    offset = [ioffset, joffset, koffset]
+                    shifted = shift(box, offset)
+                    boxes += [(offset, shifted)]
+    return boxes
+
+
 def compute_overlaps(patches, domain_box):
+    from pyphare.core.phare_utilities import np_array_ify
+
+    overlaps = []
+    zero_offset = [0] * domain_box.ndim if domain_box.ndim > 1 else 0
+
+    for ip, refPatch in enumerate(patches):
+        for cmpPatch in patches[ip:]:
+
+            for ref_pdname, ref_pd in refPatch.patch_datas.items():
+                cmp_pd = cmpPatch.patch_datas[ref_pdname]
+
+                gb_ref = ref_pd.ghost_box
+                gb_cmp = cmp_pd.ghost_box
+
+                for offset, shifted_cmp in possible_periodic_shifts(gb_cmp, domain_box):
+                    overlap = gb_ref * shifted_cmp
+                    if overlap is not None and not np.all(
+                        overlap.shape == gb_ref.shape
+                    ):
+                        if ref_pd.quantity == "field":
+                            overlap = toFieldBox(overlap, ref_pd)
+
+                        overlaps.append(
+                            {
+                                "pdatas": (ref_pd, cmp_pd),
+                                "patches": (refPatch, cmpPatch),
+                                "box": overlap,
+                                "offset": (zero_offset, np_array_ify(offset)),
+                            }
+                        )
+                        if offset != zero_offset:
+                            other_overlap = (
+                                boxm.shift(gb_ref, -np_array_ify(offset)) * gb_cmp
+                            )
+                            overlaps.append(
+                                {
+                                    "pdatas": (ref_pd, cmp_pd),
+                                    "patches": (refPatch, cmpPatch),
+                                    "box": other_overlap,
+                                    "offset": (-np_array_ify(offset), zero_offset),
+                                }
+                            )
+
+    return overlaps
+
+
+def compute_overlaps_(patches, domain_box):
     """
     returns a list of overlaps for all patch datas in given patches
     and for a domain box. An overlap is defined as an intersection of
@@ -246,7 +369,7 @@ def borders_per(patch):
     ]
 
     for patch_i, ref_patch in enumerate(border_patches):
-        in_sides = borders_per_patch[ref_patch]
+        in_sides = "".join(sorted(borders_per_patch[ref_patch]))
         assert in_sides in shifts
 
         for ref_pdname, ref_pd in ref_patch.patch_datas.items():
@@ -336,36 +459,41 @@ def get_periodic_list(patches, domain_box, n_ghosts):
             shift_patch(first_patch, domain_box.shape)
             sorted_patches.append(first_patch)
 
+        return sorted_patches
+
+    dbu = domain_box.upper
+
     if dim == 2:
         sides = {
-            "bottom": Box([0, 0], [domain_box.upper[0], 0]),
-            "top": Box(
-                [0, domain_box.upper[1]], [domain_box.upper[0], domain_box.upper[1]]
-            ),
-            "left": Box([0, 0], [0, domain_box.upper[1]]),
-            "right": Box(
-                [domain_box.upper[0], 0], [domain_box.upper[0], domain_box.upper[1]]
-            ),
+            "left": Box([0, 0], [0, dbu[1]]),
+            "right": Box([dbu[0], 0], [dbu[0], dbu[1]]),
+            "bottom": Box([0, 0], [dbu[0], 0]),
+            "top": Box([0, dbu[1]], [dbu[0], dbu[1]]),
         }
 
-        shifts = periodicity_shifts(domain_box)
+    else:
+        sides = {
+            "left": Box([0, 0, 0], [0, dbu[1], dbu[2]]),
+            "right": Box([dbu[0], 0, 0], [dbu[0], dbu[1], dbu[2]]),
+            "bottom": Box([0, 0, 0], [dbu[0], 0, dbu[2]]),
+            "top": Box([0, dbu[1], 0], [dbu[0], dbu[1], dbu[2]]),
+            "front": Box([0, 0, 0], [dbu[0], dbu[1], 0]),
+            "back": Box([0, 0, dbu[2]], [dbu[0], dbu[1], dbu[2]]),
+        }
 
-        def borders_per(box):
-            return "".join(
-                [key for key, side in sides.items() if box * side is not None]
-            )
+    shifts = periodicity_shifts(domain_box)
 
-        for patch in patches:
-            in_sides = borders_per(boxm.grow(patch.box, n_ghosts))
+    def borders_per(box):
+        return "".join([key for key, side in sides.items() if box * side is not None])
 
-            if in_sides in shifts:  # in_sides might be empty, so no borders
-                for shift in shifts[in_sides]:
-                    patch_copy = copy(patch)
-                    shift_patch(patch_copy, shift)
-                    sorted_patches.append(patch_copy)
+    for patch in patches:
+        in_sides = "".join(sorted(borders_per(boxm.grow(patch.box, n_ghosts))))
 
-    if dim == 3:
-        raise ValueError("not yet implemented")
+        if in_sides in shifts:  # in_sides might be empty, so no borders
+            for shift in shifts[in_sides]:
+                patch_copy = copy(patch)
+                shift_patch(patch_copy, shift)
+                sorted_patches.append(patch_copy)
 
     return sorted_patches
 
@@ -470,18 +598,7 @@ def level_ghost_boxes(hierarchy, quantities, levelNbrs=[], time=None):
                     check_patches = patches
 
                 for gabox in ghostAreaBoxes:
-                    remaining = gabox - check_patches[0].box
-
-                    for patch in check_patches[1:]:
-                        tmp = []
-                        remove = []
-                        for i, rem in enumerate(remaining):
-                            if rem * patch.box is not None:
-                                remove.append(i)
-                                tmp += rem - patch.box
-                        for rm in reversed(remove):
-                            del remaining[rm]
-                        remaining += tmp
+                    remaining = gabox - [p.box for p in check_patches]
 
                     if ilvl not in lvl_gaboxes:
                         lvl_gaboxes[ilvl] = {}

pyphare/pyphare/pharesee/hierarchy/hierarchy.py

@@ -561,6 +561,8 @@ def plot(self, **kwargs):
             return self.plot1d(**kwargs)
         elif self.ndim == 2:
             return self.plot2d(**kwargs)
+        elif self.ndim == 3:
+            raise NotImplementedError("cannot make 3D plots")
 
     def dist_plot(self, **kwargs):
         """