[geom] Support periodic boundaries with inverse order in Mesh

phi/geom/_mesh.py

@@ -202,7 +202,8 @@ def faces_to_vertices(self, values: Tensor, reduce=sum):
     @cached_property
     def _cell_deltas(self):
         bounds = bounding_box(self.vertices)
-        is_periodic = dim_mask(self.vector.item_names, self.periodic)
+        periodic = {dim[:-len('[::-1]')] if dim.endswith('[::-1]') else dim: dim.endswith('[::-1]') for dim in self.periodic}
+        is_periodic = dim_mask(self.vector.item_names, tuple(periodic))
         return pairwise_distances(self.center, format=self.cell_connectivity, periodic=is_periodic, domain=(bounds.lower, bounds.upper))
 
     @cached_property
@@ -646,8 +647,9 @@ def mesh(vertices: Geometry | Tensor,
     periodic_dims = []
     if periodic is not None:
         periodic_dims = [s.strip() for s in periodic.split(',') if s.strip()]
-        assert all(p in vertices.vector.item_names for p in periodic_dims), f"Periodic boundaries must be named after axes, e.g. {vertices.vector.item_names} but got {periodic}"
-        for base in periodic_dims:
+        periodic_base = [p[:-len('[::-1]')] if p.endswith('[::-1]') else p for p in periodic_dims]
+        assert all(p in vertices.vector.item_names for p in periodic_base), f"Periodic boundaries must be named after axes, e.g. {vertices.vector.item_names} but got {periodic}"
+        for base in periodic_base:
             assert base+'+' in boundaries and base+'-' in boundaries, f"Missing boundaries for periodicity '{base}'. Make sure '{base}+' and '{base}-' are keys in boundaries dict, got {tuple(boundaries)}"
     return Mesh(vertices, elements, element_rank, boundaries, periodic_dims, face_format=face_format, max_cell_walk=max_cell_walk)
 
@@ -675,9 +677,10 @@ def build_faces(vertices: Tensor,  # (vertices:i, vector)
     n_e = instance(elements).size
     # --- Periodic: map vertices of boundary+ to the corresponding vertex in boundary- ---
     vertex_id = np.arange(instance(vertices).size)
-    for dim in periodic:
-        vertex_id[np.concatenate(boundaries[dim+'+'])] = np.concatenate(boundaries[dim+'-'])
-    is_periodic = dim_mask(vertices.vector.item_names, periodic)
+    periodic = {dim[:-len('[::-1]')] if dim.endswith('[::-1]') else dim: dim.endswith('[::-1]') for dim in periodic}
+    for dim, is_flipped in periodic.items():
+        vertex_id[np.concatenate(boundaries[dim+'+'])] = np.concatenate(boundaries[dim+'-'])[::-1] if is_flipped else np.concatenate(boundaries[dim+'-'])
+    is_periodic = dim_mask(vertices.vector.item_names, tuple(periodic))
     # --- element-facet and facet-vertex matrix. A facet describes a single oriented face of an element, i.e. shared faces get two entries. ---
     if element_rank == 2:  # edges are the lines between neighbor vertices in the vertex lists + the edge last-to-first
         v1 = stored_indices(elements).index[dual(elements).name].numpy()
@@ -737,7 +740,7 @@ def build_faces(vertices: Tensor,  # (vertices:i, vector)
     f_f: csr_matrix = f_v @ f_v_b.T >= element_rank  # symmetric
     f_f.setdiag(0)
     f_f.eliminate_zeros()
-    f_f.data = np.arange(1, n_f+1)  # equal to f_f.nonzero()[0]
+    f_f.data = f_f.nonzero()[0] + 1
     e_e = e_f @ f_f @ e_f_b.T  # stores the outgoing facet_index+1 for each element pair
     shared_edge = wrap(e_e, instance(elements).without_sizes() & dual) - 1
     shared_edge = to_format(shared_edge, face_format)