Partial refactoring of mesh representation

I now evaluate the nodes in the poloidal plane eagerly (but the 3D
field traces will still be done lazily). I have refactored mesh.py,
hypnotoad_interface.py, nektar_writer.py, meshio_writer.py, and
wall.py to use this new approach. Everything else has just had some
renaming done and won't actually run. None of the tests have been
refactored yet.

neso_fame/coordinates.py

@@ -15,6 +15,7 @@
 from decimal import ROUND_HALF_UP, Context, Decimal
 from enum import Enum
 from functools import wraps
+from types import EllipsisType
 from typing import (
     Callable,
     Concatenate,
@@ -117,7 +118,13 @@ def approx_eq(
 
 
 CoordIndex = int | tuple[int, ...]
-IndexSlice = int | slice | list[int | tuple[int, ...]] | tuple[int | slice, ...]
+IndexSlice = (
+    int
+    | slice
+    | EllipsisType
+    | list[int | tuple[int, ...]]
+    | tuple[int | slice | EllipsisType, ...]
+)
 
 
 @dataclass
@@ -163,12 +170,28 @@ def __getitem__(self, idx: CoordIndex) -> SliceCoord:
         x1, x2 = np.broadcast_arrays(self.x1, self.x2)
         return SliceCoord(float(x1[idx]), float(x2[idx]), self.system)
 
+    @property
+    def get(self) -> _CoordWrapper[SliceCoords]:
+        """A view of this data which can be sliced.
+
+        :method:`~neso_fame.coordinates.SliceCoords.__getitem__` and
+        returns scalar :class:`~neso_fame.coordinates.SliceCoord`
+        objects and will raise an error if the index does not
+        correspond to a scalar value. Indexing the object returned by
+        this property will return a
+        :class:`~neso_fame.coordinates.SliceCoords` object instead,
+        potentially containing an array of coordinates.
+
+        """
+        return _CoordWrapper(self)
+
+    def _get(self, index: IndexSlice) -> SliceCoords:
+        x1, x2 = np.broadcast_arrays(self.x1, self.x2)
+        return SliceCoords(x1[index], x2[index], self.system)
+
     def get_set(self, index: IndexSlice) -> FrozenCoordSet[SliceCoord]:
         """Get a set of individual point objects from the collection."""
-        x1, x2 = np.broadcast_arrays(self.x1, self.x2)
-        return FrozenCoordSet(
-            SliceCoords(x1[index], x2[index], self.system).iter_points()
-        )
+        return FrozenCoordSet(self._get(index).iter_points())
 
     def round_to(self, figures: int = 8) -> SliceCoords:
         """Round coordinate values to the desired number of significant figures."""
@@ -318,12 +341,28 @@ def __getitem__(self, idx: CoordIndex) -> Coord:
         x1, x2, x3 = np.broadcast_arrays(self.x1, self.x2, self.x3)
         return Coord(float(x1[idx]), float(x2[idx]), float(x3[idx]), self.system)
 
+    @property
+    def get(self) -> _CoordWrapper[Coords]:
+        """A view of this data which can be sliced.
+
+        :method:`~neso_fame.coordinates.Coords.__getitem__` and
+        returns scalar :class:`~neso_fame.coordinates.Coord`
+        objects and will raise an error if the index does not
+        correspond to a scalar value. Indexing the object returned by
+        this property will return a
+        :class:`~neso_fame.coordinates.Coords` object instead,
+        potentially containing an array of coordinates.
+
+        """
+        return _CoordWrapper(self)
+
+    def _get(self, index: IndexSlice) -> Coords:
+        x1, x2, x3 = np.broadcast_arrays(self.x1, self.x2, self.x3)
+        return Coords(x1[index], x2[index], x3[index], self.system)
+
     def get_set(self, index: IndexSlice) -> FrozenCoordSet[Coord]:
         """Get a set of individual point objects from the collection."""
-        x1, x2, x3 = np.broadcast_arrays(self.x1, self.x2, self.x3)
-        return FrozenCoordSet(
-            Coords(x1[index], x2[index], x3[index], self.system).iter_points()
-        )
+        return FrozenCoordSet(self._get(index).iter_points())
 
     def to_coord(self) -> Coord:
         """Convert the object to a `Coord` object.
@@ -357,9 +396,35 @@ def to_slice_coords(self) -> SliceCoords:
 
 
 C = TypeVar("C", Coord, SliceCoord)
+Cs = TypeVar("Cs", Coords, SliceCoords)
 T = TypeVar("T")
 
 
+@dataclass(frozen=True)
+class _CoordWrapper(Generic[Cs]):
+    """Simple class to allow slicing coordinate objects.
+
+    :method:`~neso_fame.coordinates.SliceCoords.__getitem__` and
+    :method:`~neso_fame.coordinates.Coords.__getitem__` are designed
+    to return scalar :class:`~neso_fame.coordinates.SliceCoord` and
+    :class:`~neso_fame.coordinates.SliceCoord` objects,
+    respectively. This wrapper class will return
+    :class:`~neso_fame.coordinates.SliceCoords` and
+    :class:`~neso_fame.coordinates.SliceCoords` instead.
+
+    Group
+    -----
+    coordinates
+
+    """
+
+    data: Cs
+
+    def __getitem__(self, index: IndexSlice) -> Cs:
+        """Get a slice of the contained data."""
+        return self.data._get(index)
+
+
 class _CoordContainer(Generic[C]):
     """Base type for approximate coordinate comparison.
 

neso_fame/element_builder.py

@@ -33,7 +33,7 @@
     FieldTracer,
     Prism,
     Quad,
-    StraightLineAcrossField,
+    straight_line_across_field,
 )
 
 
@@ -483,7 +483,7 @@ def make_quad_for_prism(
         # aligned but then there would be no guarantee that the centre
         # would remain within the walls.
         q = Quad(
-            StraightLineAcrossField(north, south),
+            straight_line_across_field(north, south),
             self._tracer,
             self._dx3,
             north_start_weight=self._vertex_start_weights.get(north, 1.0),

neso_fame/generators.py

@@ -39,8 +39,8 @@
     PrismMesh,
     Quad,
     QuadMesh,
-    StraightLineAcrossField,
     control_points,
+    straight_line_across_field,
 )
 from neso_fame.wall import (
     Connections,
@@ -98,7 +98,7 @@ def _get_vec(
 
 def _constrain_to_plain(
     field: FieldTrace,
-    shape: StraightLineAcrossField,
+    shape: straight_line_across_field,
     north_bound: BoundType,
     south_bound: BoundType,
 ) -> FieldTrace:
@@ -231,7 +231,7 @@ def field_aligned_2d(
         connectivity = _ordered_connectivity(num_nodes)
 
     def make_quad(node1: int, node2: int) -> Quad:
-        shape = StraightLineAcrossField(lower_dim_mesh[node1], lower_dim_mesh[node2])
+        shape = straight_line_across_field(lower_dim_mesh[node1], lower_dim_mesh[node2])
         north_bound = node1 in (0, num_nodes - 1) and conform_to_bounds
         south_bound = node2 in (0, num_nodes - 1) and conform_to_bounds
         return Quad(
@@ -399,7 +399,7 @@ def bound_type(edges: list[NodePair]) -> BoundType:
 
     @cache
     def make_quad(node1: Index, node2: Index) -> Quad:
-        shape = StraightLineAcrossField(lower_dim_mesh[node1], lower_dim_mesh[node2])
+        shape = straight_line_across_field(lower_dim_mesh[node1], lower_dim_mesh[node2])
         if conform_to_bounds:
             north_bound = boundary_nodes.get(node1, False)
             south_bound = boundary_nodes.get(node2, False)
@@ -834,6 +834,7 @@ def _validate_wall_elements(
     changed to prevent self-intersecting elements.
 
     """
+    # FIXME: My refactor of how I represent meshes will mean I pretty much need to rewrite this.
     # TODO: Should I validate internal elements too? If I flatten one
     # then that could end up makign a further element invalid, which
     # sounds unpleasant to have to deal with...
@@ -852,6 +853,8 @@ def _validate_wall_elements(
         # Try merging with adjacent triangles (which haven't already
         # been merged with another element, which would remove them
         # from new_elements)
+
+        # FIXME: Will need to change how I map faces to elements; just use pairs of points?
         merge_candidates = frozenset(
             item
             for item in itertools.chain.from_iterable(
@@ -1056,6 +1059,7 @@ def whole_line_in_tokamak(start: SliceCoord, wall: Sequence[WallSegment]) -> boo
         )
         if corners_within_vessel(corners)
     ]
+    # Probably more efficient just to iterate over inner regions
     if mesh_to_core:
         core_elements = list(
             itertools.starmap(
@@ -1069,7 +1073,10 @@ def whole_line_in_tokamak(start: SliceCoord, wall: Sequence[WallSegment]) -> boo
         inner_bounds = frozenset(factory.innermost_quads())
     if mesh_to_wall:
         # FIXME: Not capturing the curves of the outermost hypnotoad quads now, for some reason.
+
+        # FIXME: Assemble coordinate pairs and mapping between these pairs and the list of Coords defining the curve
         plasma_points = [tuple(p) for p in factory.outermost_vertices()]
+        # FIXME: Assemble list of Coords (one for each wall segment) and also coordinate pairs?
         if wall_resolution is not None:
             target = _average_poloidal_spacing(hypnotoad_poloidal_mesh)
             wall: list[Point2D] = adjust_wall_resolution(
@@ -1095,17 +1102,6 @@ def whole_line_in_tokamak(start: SliceCoord, wall: Sequence[WallSegment]) -> boo
             + list(periodic_pairwise(iter(range(n, n + len(plasma_points)))))
         )
         info.set_holes([tuple(hypnotoad_poloidal_mesh.equilibrium.o_point)])
-        # It may be worth adding the option to start merging quads
-        # with a really weird aspect ratio (i.e, leading away from the
-        # x-point towards the centre). When aspect ratio of two of
-        # them gets too high, create a triangle. Would need to be
-        # careful about direction though. Would this also be useful
-        # around seperatrix? Want some increase resolution there, but
-        # not necessarily too much. Might a mapping between
-        # hypntoad-generate points and FAME elements be useful here?
-        # Or maybe I can extend _element_corners to be able to do
-        # this? Might also be good to reduce perpendicular resolution
-        # in PFR. Try tracing out from x-points? Quite hard to coordinate.
         wall_mesh = triangle.build(
             info, allow_volume_steiner=True, allow_boundary_steiner=False
         )
@@ -1115,6 +1111,7 @@ def whole_line_in_tokamak(start: SliceCoord, wall: Sequence[WallSegment]) -> boo
             wall_mesh_points[:, 0], wall_mesh_points[:, 1], system
         )
         initial: tuple[list[Prism], frozenset[Quad]] = ([], frozenset())
+        # FIXME: Take coordinate pairs, check if either of them correspond to curves and use those or else just use the pair. If any of the pairs are from the wall, create a boundary item as well.
         initial_wall_elements, initial_outer_bounds = reduce(
             lambda left, right: (left[0] + [right[0]], left[1] | right[1]),
             (