3d C++ / python

.github/workflows/cmake_macos.yml

@@ -87,7 +87,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,7 +157,35 @@ 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
 
@@ -246,7 +287,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 +377,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
 
@@ -471,18 +517,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

@@ -556,11 +556,44 @@ def plot2d(self, **kwargs):
 
         return fig, ax
 
+    def plot3d(self, **kwargs):
+        """!HAX!"""
+        time = kwargs.get("time", self._default_time())
+        usr_lvls = kwargs.get("levels", self.levelNbrs(time))
+        default_qty = None
+        if len(self.quantities()) == 1:
+            default_qty = self.quantities()[0]
+        qty = kwargs.get("qty", default_qty)
+        for lvl_nbr, lvl in self.levels(time).items():
+            if lvl_nbr not in usr_lvls:
+                continue
+            for patch in self.level(lvl_nbr, time).patches:
+                pdat = patch.patch_datas[qty]
+                primals = pdat.primal_directions()
+                if primals[0]:
+                    pdat._x = pdat.x[:-1]
+                if primals[1]:
+                    pdat._y = pdat.y[:-1]
+                pdat.dataset = pdat.dataset[:, :, int(pdat.ghost_box.shape[2] / 2)]
+                patch.box.lower = patch.box.lower[:-1]
+                patch.box.upper = patch.box.upper[:-1]
+                patch.box.ndim = 2
+
+                pdat.ghost_box.lower = pdat.ghost_box.lower[:-1]
+                pdat.ghost_box.upper = pdat.ghost_box.upper[:-1]
+                pdat.ghost_box.ndim = 2
+                pdat.size = np.copy(pdat.ghost_box.shape)
+                pdat.layout.dl = pdat.layout.dl[:-1]
+
+        return self.plot2d(**kwargs)  # ¯\_(ツ)_/¯
+
     def plot(self, **kwargs):
         if self.ndim == 1:
             return self.plot1d(**kwargs)
         elif self.ndim == 2:
             return self.plot2d(**kwargs)
+        elif self.ndim == 3:
+            return self.plot3d(**kwargs)
 
     def dist_plot(self, **kwargs):
         """

pyphare/pyphare/pharesee/run/utils.py

@@ -359,7 +359,7 @@ def _get_rank(patchdatas, patch_id, **kwargs):
     layout = reference_pd.layout
     centering = "dual"
     nbrGhosts = layout.nbrGhosts(layout.interp_order, centering)
-    shape = grow(reference_pd.box, [nbrGhosts] * 2).shape
+    shape = grow(reference_pd.box, [nbrGhosts] * ndim).shape
 
     if ndim == 1:
         pass

res/amr/splitting.yml

@@ -99,3 +99,43 @@ dimension_2:
       delta:  1 2
       weight: .140625 .09375 .0625 .0234375 .015625 .00390625
 
+
+dimension_3:
+  interp_1:
+    N_particles_6:
+      delta:  .966431
+      weight: .166666
+
+    N_particles_12:
+      delta:  .74823
+      weight: .083333
+
+    N_particles_27:
+      delta:  1
+      weight: .125 .0625
+
+  interp_2:
+    N_particles_6:
+      delta:  1.149658
+      weight: .166666
+
+    N_particles_12:
+      delta:  .888184
+      weight: .083333
+
+    N_particles_27:
+      delta:  1.111333
+      weight: .099995 .055301
+
+  interp_3:
+    N_particles_6:
+      delta:  1.312622
+      weight: .166666
+
+    N_particles_12:
+      delta:  1.012756
+      weight: .083333
+
+    N_particles_27:
+      delta:  1.276815
+      weight: .104047 .05564