3d-ish

pyphare/pyphare/core/box.py

@@ -67,6 +67,9 @@ def __eq__(self, other):
         return isinstance(other, Box) and (self.lower == other.lower).all() and (self.upper == other.upper).all()
 
     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)
 
@@ -183,5 +186,26 @@ 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)])]

pyphare/pyphare/core/gridlayout.py

@@ -275,6 +275,22 @@ 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):
         """
         from a qty and a direction, returns a 1d array containing

pyphare/pyphare/pharein/maxwellian_fluid_model.py

@@ -142,104 +142,54 @@ def to_dict(self):
 #------------------------------------------------------------------------------
 
     def validate(self, sim):
-        import itertools
         import numpy as np
         from ..core.gridlayout import GridLayout, directions as all_directions
         from ..core.box import Box
         from pyphare.pharein import py_fn_wrapper
 
-        dl = np.array(sim.dl)
         directions = all_directions[:sim.ndim]
-
         domain_box = Box([0] * sim.ndim, np.asarray(sim.cells) - 1)
         assert len(sim.origin) == domain_box.ndim
+        ndim = domain_box.ndim
 
         def _primal_directions(qty):
             return [1 if yee_element_is_primal(qty, direction) else 0 for direction in directions]
 
-        def _nbrGhosts(qty, layout, primal_directions):
+        def _nbrGhosts(layout, qty, primal_directions):
             return [
-                layout.nbrGhosts(sim.interp_order, "primal" if primal_directions[dim] else "dual") for dim in range(sim.ndim)
+                layout.nbrGhosts(sim.interp_order, "primal" if primal_directions[dim] else "dual") for dim in range(ndim)
             ]
 
-        def qty_shifts(qty, nbrGhosts):
-          shifts = []
-          for dim in range(sim.ndim):
-              if yee_element_is_primal(qty, directions[dim]):
-                  shifts += [np.arange(-nbrGhosts[dim], nbrGhosts[dim] + 1) * dl[dim]]
-              else:
-                  shifts += [np.arange(-nbrGhosts[dim], nbrGhosts[dim]) *  dl[dim] + (.5 * dl[dim])]
-          return shifts
+        def split_to_columns(ar):
+            ar = ar.reshape(-1, ndim) # drop outer arrays if any
+            return [ar[:,dim] for dim in range(ndim)]
 
-        def _1d_check(layout, nbrGhosts, primal_directions, qty, fn):
-            fnX = py_fn_wrapper(fn)(layout.yeeCoordsFor(qty, "x"))
-            include_border = primal_directions[0]
-            np.testing.assert_allclose(fnX[:nbrGhosts * 2 + include_border] , fnX[-(nbrGhosts * 2) - include_border:], atol=1e-15)
-            return (
-                np.allclose(fnX[:nbrGhosts * 2 + include_border] , fnX[-(nbrGhosts * 2) - include_border:], atol=1e-15)
-            )
+        def _nd_check(layout, qty, nbrGhosts, fn):
+            from pyphare.pharesee.geometry import hierarchy_overlaps
+            from pyphare.pharesee.hierarchy import hierarchy_from_box
+            from ..core import box as boxm
+
+            hier = hierarchy_from_box(domain_box, nbrGhosts)
+            coords = layout.meshCoords(qty)
+
+            for ilvl, overlaps in hierarchy_overlaps(hier).items():
+                for overlap in overlaps:
+                    (pd1, pd2), box, offsets = overlap["pdatas"], overlap["box"], overlap["offset"]
+                    slice1 = boxm.select(coords, boxm.amr_to_local(box, boxm.shift(pd1.ghost_box, offsets[0])))
+                    slice2 = boxm.select(coords, boxm.amr_to_local(box, boxm.shift(pd2.ghost_box, offsets[1])))
+                    if not np.allclose(fn(*split_to_columns(slice1)) , fn(*split_to_columns(slice2)), atol=1e-15):
+                        return False
+            return True
 
 
-        def split_to_columns(ar):
-            ar = ar.reshape(-1, sim.ndim) # drop outer arrays if any
-            return [ar[:,dim] for dim in range(sim.ndim)]
-
-
-        def _2d_check(nbrGhosts, primal_directions, qty, fn, shift):
-            layout = GridLayout(domain_box, sim.origin + (shift * sim.dl), sim.dl, sim.interp_order)
-            x = layout.yeeCoordsFor(qty, "x")[nbrGhosts[0]:-(nbrGhosts[0]-1)-(primal_directions[0])]
-            y = layout.yeeCoordsFor(qty, "y")[nbrGhosts[1]:-(nbrGhosts[1]-1)-(primal_directions[1])]
-            X,Y = np.meshgrid(x,y,indexing="ij")
-            coords = np.array([X.flatten(),Y.flatten()]).T
-            top = coords[:len(x)]
-            bot = coords[-len(x):]
-            left  = coords[0 :: len(x)]
-            right = coords[len(x) - 1 :: len(x)]
-            return (
-                np.allclose(fn(*split_to_columns(left)) , fn(*split_to_columns(right)), atol=1e-15) and
-                np.allclose(fn(*split_to_columns(top))  , fn(*split_to_columns(bot)), atol=1e-15)
-            )
-
-        def _3d_check(nbrGhosts, primal_directions, qty, fn, shift):
-            def get_3d_faces(M):
-                def get_face(M, dim, front_side):
-                    return M[tuple((0 if front_side else -1) if i == dim else slice(None) for i in range(M.ndim))]
-                return [get_face(M, dim, front_side) for dim in range(sim.ndim) for front_side in (True, False)]
-
-            layout = GridLayout(domain_box, sim.origin + (shift * sim.dl), sim.dl, sim.interp_order)
-            x = layout.yeeCoordsFor(qty, "x")[nbrGhosts[0]:-(nbrGhosts[0]-1)-(primal_directions[0])]
-            y = layout.yeeCoordsFor(qty, "y")[nbrGhosts[1]:-(nbrGhosts[1]-1)-(primal_directions[1])]
-            z = layout.yeeCoordsFor(qty, "z")[nbrGhosts[2]:-(nbrGhosts[2]-1)-(primal_directions[2])]
-            coords = np.array([X.flatten(), Y.flatten(), Z.flatten()]).T.reshape((len(x), len(y), len(z), sim.ndim))
-            left, right, top, bot, front, back = get_3d_faces(coords)
-            return (
-                np.allclose(fn(split_to_columns(left)) , fn(split_to_columns(right)), atol=1e-15) and
-                np.allclose(fn(split_to_columns(top))  , fn(split_to_columns(bot)), atol=1e-15) and
-                np.allclose(fn(split_to_columns(front)), fn(split_to_columns(back)), atol=1e-15)
-            )
 
         def periodic_function_check(vec_field, dic):
             valid = True
             for xyz in ["x", "y", "z"]:
                 qty = vec_field + xyz
-                fn = dic[qty]
-
                 layout = GridLayout(domain_box, sim.origin, sim.dl, sim.interp_order)
-                primal_directions = _primal_directions(qty)
-                nbrGhosts = _nbrGhosts(qty, layout, primal_directions)
-
-                if sim.ndim == 1:
-                    valid &= _1d_check(layout, nbrGhosts[0], primal_directions, qty, fn)
-
-                if sim.ndim == 2:
-                    permutations = itertools.product(*qty_shifts(qty, nbrGhosts))
-                    valid &= all([_2d_check(nbrGhosts, primal_directions, qty, fn, np.asarray(perm)) for perm in permutations])
-
-                if sim.ndim == 3:
-                    # untested block - add in during 3d/splitting PR https://github.com/PHAREHUB/PHARE/pull/314
-                    permutations = itertools.product(*qty_shifts(qty, nbrGhosts))
-                    valid &= all([_3d_check(nbrGhosts, primal_directions, qty, fn, np.asarray(perm)) for perm in permutations])
-
+                nbrGhosts = _nbrGhosts(layout, qty, _primal_directions(qty))
+                valid &= _nd_check(layout, qty, nbrGhosts, fn=dic[qty])
             return valid
 
         if sim.boundary_types[0] == "periodic":

pyphare/pyphare/pharein/simulation.py

@@ -194,9 +194,13 @@ def check_boundaries(ndim, **kwargs):
     3: [4, 5, 8, 9, 25]
   },
   3: {
-    1: [6, 7, 8, 9, 12, 13, 14, 15, 18, 19, 20, 21, 26, 27],
-    2: [6, 7, 8, 9, 12, 13, 14, 15, 18, 19, 20, 21, 26, 27, 64],
-    3: [6, 7, 8, 9, 12, 13, 14, 15, 18, 19, 20, 21, 26, 27, 125]
+    1: [6, 12],
+    2: [6, 12],
+    3: [6, 12]
+    # full below
+    # 1: [6, 7, 8, 9, 12, 13, 14, 15, 18, 19, 20, 21, 26, 27],
+    # 2: [6, 7, 8, 9, 12, 13, 14, 15, 18, 19, 20, 21, 26, 27, 64],
+    # 3: [6, 7, 8, 9, 12, 13, 14, 15, 18, 19, 20, 21, 26, 27, 125]
   }
 } # Default refined_particle_nbr per dim/interp is considered index 0 of list
 def check_refined_particle_nbr(ndim, **kwargs):

pyphare/pyphare/pharesee/geometry.py

@@ -49,13 +49,25 @@ 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)),
+        }
+
+    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)),
         }
 
-    raise ValueError("Unhandeled dimension")
+
+
 
 def touch_domain_border(box, domain_box, border):
     if border == "upper":
@@ -70,13 +82,18 @@ def periodicity_shifts(domain_box):
 
     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)],
@@ -104,13 +121,33 @@ def periodicity_shifts(domain_box):
           "topleftright" : [
               *shifts["leftright"], *shifts["top"],
               shifts["topleft"][-1], shifts["topright"][-1]],
-          "bottomtopleftright" : [ # one patch covers domain
+          "leftrightbottomtop" : [ # one patch covers domain
               *shifts["bottomleft"], *shifts["topright"],
               shifts["bottomright"][-1], shifts["topleft"][-1]]
         })
 
     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"],
+            ]
+        })
+
+    shifts = {"".join(sorted(k)) : l for k,l in shifts.items()}
 
     return shifts
 
@@ -206,7 +243,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():
@@ -249,6 +286,7 @@ def hierarchy_overlaps(hierarchy, time=0):
 
 
 
+
 def get_periodic_list(patches, domain_box, n_ghosts):
     """
     given a list of patches and a domain box the function
@@ -278,38 +316,50 @@ def get_periodic_list(patches, domain_box, n_ghosts):
             shift_patch(first_patch, domain_box.shape)
             sorted_patches.append(first_patch)
 
-    if dim == 2:
+        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:
+    def borders_per(box):
+        return "".join([key for key, side in sides.items() if box * side is not None])
 
-            in_sides = borders_per(boxm.grow(patch.box, n_ghosts))
+    for patch in patches:
 
-            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)
+        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
 
 
 
 
+
+
 def ghost_area_boxes(hierarchy, quantities, levelNbrs=[], time=0):
     """
     this function returns boxes representing ghost cell boxes for all levels
@@ -412,18 +462,7 @@ def level_ghost_boxes(hierarchy, quantities, levelNbrs=[], time=None):
 
                 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] = {}