Repository: gerlero/foamlib

Issue:
Performance problems when reading large non-dictionary OpenFoam files
Hi again @gerlero,

The other you change you made is working great now. I've come across one other thing that I was hoping you might be able to look at.

I am also trying to parse files for the mesh which are essentially just a large array. For these I am having performance problems once the the number of rows start to appraoch 1M. Interestingly, this happens not on parsing but on extracting the values.

What is interesting is that these files are way smaller than some of the field files that are being processed easily. My guess it's to do with the values consisting of one large list.

I am using foamlib 0.8.10 and python 3.11

Thanks!

**FoamFile (real version has all numbers):**

> /*--------------------------------*- C++ -*----------------------------------*\
> |       o          |                                                          |
> |    o     o       | HELYX(R) Core: Open-source CFD for Enterprise            |
> |   o   O   o      | Version: 4.1.1                                           |
> |    o     o       | ENGYS Ltd. <http://engys.com/>                           |
> |       o          |                                                          |
> \*---------------------------------------------------------------------------*/
> FoamFile
> {
>     version     2.0;
>     format      binary;
>     class       regIOobject;
>     arch        "LSB;label=32;scalar=64";
>     location    "constant/polyMesh";
>     object      cellZones;
> }
> // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
> 
> 1000000
> (
>     1
>     2
>     3
>     4
>     ...
>     999997
>     999998
>     999999
>     1000000    
> )
> // ************************************************************************* //
> 

```
from foamlib import FoamFile

path = 'path_to_foam_file...'
U = FoamFile(path)
print(U[None])
```

Hints from issue comments:
@gerlero: > My guess it's to do with the values consisting of one large list.

@cjbraley Yeah, only `nonuniform` fields currently go through the `re`-based "fast path", while this is just a plain list as far as foamlib is concerned and as such goes through the normal `pyparsing`-based parser which is slower. 

It should definitely be possible to also "fast-path" a simple list of integers like the one you show. I guess it would be interesting to see which other types would really benefit from being on the fast path (I'm thinking the one in the `points` file, etc.) and come up with a solution that handles all of those.@cjbraley: @gerlero ,

Thanks for your response. That makes sense and is what I suspected was probably happening. For the work I'm doing, it is simple lists of integers, like the one I shared, but also vectors (see below).

Thankfully you have already done the hard work parsing the more complicated files, while it would be nice to use foamlib for everything, if it's difficult for you to change, I can also write a simple parser for these files.

**faces**
```
446448146
(
3(904040 904479 924424)
3(904479 904039 924424)
3(904481 904479 924424)
3(904040 904481 924424)
3(55880317 904479 904039)
3(904039 924424 55880317)
3(904479 55880317 924424)
...
3(924424 55850535 904479)
3(924424 904481 55850535)
3(904481 904479 55850535)
3(65097588 904481 924424)
3(65097588 924424 904040)
)
```

**points**
```
69229432
(
(1.9651 3.571886497513121e-18 0.468999999999996)
(1.964838616647891 -0.006931990539489501 0.4689999999999958)
(1.964037752895968 -0.01391940630148043 0.4689999999999963)
(1.962548278381021 -0.02147948475297981 0.4689999999999959)
(1.960348384467494 -0.02913366520348652 0.4689999999999965)
...
(1.092260254127397 -0.3151451624161979 0.05406851296852164)
(2.461613398668586 -0.5905097235548232 0.0415838960382229)
(0.2032821554166777 -0.1886504205189976 0.3134021318602374)
(0.767235839694649 -0.5245151548925305 0.03938924639997934)
)
```


@gerlero: @cjbraley Then we should ideally fast-path:

- [x] list[int]
- [x] list[3-list[int]]
- [x] list[3-list[float]] (~ `vector` field)
- [x] also on the fast path now:  list[float] | list[4-list[int|float]]

I'll see what I can come up with.@gerlero: BTW, fast-pathing these would also imply returning them as NumPy arrays (treating these as Python (nested) lists would have a performance impact and take just too much memory).@cjbraley: NumPy arrays would be perfect. We are trying to end up with dataframes as an end result.@gerlero: Implemented in `0.9.0`. I haven't really done any performance testing though; please re-open if it doesn't work as expected.@cjbraley: Thanks again for looking at this @gerlero. Very impressed at your speed!

I have been testing the latest release today. I am now able to parse scalar (list[int]) files up to ~200M lines and vector (list[3-list[float]] ) to around ~25M, anything larger crashes. It goes from finishing quickly to crashing very quickly as you increase the number of rows. This is frustratingly close to the about 500M and 70M we need.

However, **I can only parse the ASCII formatted files**, the equivalent binary formatted files fail with the following error. It seems like it is not recognising when a file is binary?
```
  File "/Users/Christopher.Braley/Documents/projects/apqx/openfoam/foamlib_constant.py", line 27, in parse_file
    print(U[None])
          ~^^^^^^
  File "/Users/Christopher.Braley/Documents/projects/apqx/openfoam/venv/lib/python3.11/site-packages/foamlib/_files/_files.py", line 246, in __getitem__
    parsed = self._get_parsed()
             ^^^^^^^^^^^^^^^^^^
  File "/Users/Christopher.Braley/Documents/projects/apqx/openfoam/venv/lib/python3.11/site-packages/foamlib/_files/_io.py", line 66, in _get_parsed
    self.__parsed = Parsed(contents)
                    ^^^^^^^^^^^^^^^^
  File "/Users/Christopher.Braley/Documents/projects/apqx/openfoam/venv/lib/python3.11/site-packages/foamlib/_files/_parsing.py", line 439, in __init__
    for parse_result in _FILE.parse_string(
                        ^^^^^^^^^^^^^^^^^^^
  File "/Users/Christopher.Braley/Documents/projects/apqx/openfoam/venv/lib/python3.11/site-packages/pyparsing/core.py", line 1219, in parse_string
    raise exc.with_traceback(None)
pyparsing.exceptions.ParseException: Expected end of text, found '('  (at char 841), (line:21, col:1)
```

Another thing I was thinking about is whether it might make sense to have standard parsing for mesh files because their format is predictable. Here is another module, which takes the approach of a MeshParser, similar to how you have the field parser. https://github.com/ApolloLV/openfoamparser/blob/master/openfoamparser/mesh_parser.py . This one is problematic for other reasons in that it tries to parse all mesh files at once and uses standard lists which kills performance. It also writes the functions based on file names but I think doing it by type would make more sense (e.g. parse_scalar_content).@gerlero: > I have been testing the latest release today. I am now able to parse scalar (list[int]) files up to ~200M lines and vector (list[3-list[float]] ) to around ~25M, anything larger crashes. It goes from finishing quickly to crashing very quickly as you increase the number of rows. This is frustratingly close to the about 500M and 70M we need.

Yeah, the parsing engine reads the entire file into memory; plus the conversions between `str`<->`bytes` (since we care about the bytes, but `pyparsing` only works with `str`s) may make it so that multiple copies of the entire file are loaded into memory. I'm not saying there are no possible optimizations here, but I don't clearly see any low-hanging fruit right now.

> However, **I can only parse the ASCII formatted files**, the equivalent binary formatted files fail with the following error. It seems like it is not recognising when a file is binary?

Yes. I didn't enable support for non-field binary-encoded lists, given I wasn't sure if it would work properly. I haven't given it much thought since, but now I believe it might work (the issue is, to work within the current constraints of foamlib's parser, whether we can predict the right shape of the stored array). 

If you're interested, please open a new issue for this specifically. It should be simple to get an implementation ready to try out.

> Another thing I was thinking about is whether it might make sense to have standard parsing for mesh files because their format is predictable. Here is another module, which takes the approach of a MeshParser, similar to how you have the field parser. https://github.com/ApolloLV/openfoamparser/blob/master/openfoamparser/mesh_parser.py . This one is problematic for other reasons in that it tries to parse all mesh files at once and uses standard lists which kills performance. It also writes the functions based on file names but I think doing it by type would make more sense (e.g. parse_scalar_content).

The current philosophy of foamlib's parser is that it uses the same generic parsing grammar for every file—no matter if it's a field, a plain dictionary, a `points` file in the mesh... you get the point. Every OpenFOAM file can be interacted with using the `FoamFile` class: even subclasses like `FoamFieldFile` only add syntactic sugar (i.e. extra properties), but are not really required. Keeping with this philosophy, I'd prefer to have the mesh files also be parseable via `FoamFile` (even if we added other convenience classes in the future for dealing with the mesh). I'd rather not change this part of foamlib if at all possible.@cjbraley: @gerlero. Very reasonable response. I completely understand wanting to keep the parsing generic, it's what makes the library very user friendly. I agree that it's probably not worth it.

I think that if you are able to get the non-field binary-encoded lists working that will hopefully get us over the line anyway. I will make a separate issue for that.

I did also have more of a question than an issue that I was hoping you would be able to tell me if I'm wasting my time on. As you can imagine, with our file sizes, the data is actually stored in S3. My ideal workflow would actually be, downloading the file into a buffer and passing it straight into foamlib for conversion, then saving the processed file back to s3, without ever writing to disk. From the testing I've done, it seems foamlib can only accept paths (str) and is dependent on the file being on disk - do I have that right? 

Thanks again@gerlero: > I did also have more of a question than an issue that I was hoping you would be able to tell me if I'm wasting my time on. As you can imagine, with our file sizes, the data is actually stored in S3. My ideal workflow would actually be, downloading the file into a buffer and passing it straight into foamlib for conversion, then saving the processed file back to s3, without ever writing to disk. From the testing I've done, it seems foamlib can only accept paths (str) and is dependent on the file being on disk - do I have that right?

Actually, while the `FoamFile` class is designed to work with actual files, there is no reason why foamlib couldn't offer a (e.g. `json`-like) `loads()` standalone function that parses directly from memory. Feel free to open an issue if a function like that make your workflow easier.

Solution patch:
diff --git a/foamlib/_files/_parsing.py b/foamlib/_files/_parsing.py
index 2f6cb80..09fa28d 100644
--- a/foamlib/_files/_parsing.py
+++ b/foamlib/_files/_parsing.py
@@ -25,6 +25,7 @@ from pyparsing import (
     LineEnd,
     Literal,
     Located,
+    NoMatch,
     Opt,
     ParserElement,
     ParseResults,
@@ -37,131 +38,194 @@ from pyparsing import (
     printables,
 )
 
-from ._types import Data, Dimensioned, DimensionSet, File, TensorKind
+from ._types import Data, Dimensioned, DimensionSet, File
 
 
-def _tensor(tensor_kind: TensorKind, *, ignore: Regex | None = None) -> Regex:
+def _numeric_list(
+    *, nested: int | None = None, ignore: Regex | None = None, force_float: bool = False
+) -> ParserElement:
+    if not force_float:
+        int_pattern = r"(?:-?\d+)"
     float_pattern = r"(?i:[+-]?(?:(?:\d+\.?\d*(?:e[+-]?\d+)?)|nan|inf(?:inity)?))"
+    spacing_pattern = (
+        rf"(?:(?:\s|{ignore.re.pattern})+)" if ignore is not None else r"\s+"
+    )
 
-    if tensor_kind == TensorKind.SCALAR:
-        ret = Regex(float_pattern)
-        ret.add_parse_action(lambda tks: [float(tks[0])])
-        return ret
+    if nested is None:
+        if not force_float:
+            int_element_pattern = int_pattern
+            int_element = common.integer
+        float_element_pattern = float_pattern
+        float_element = common.ieee_float
+    else:
+        if not force_float:
+            int_element_pattern = rf"(?:(?:{nested})?{spacing_pattern}?\({spacing_pattern}?(?:{int_pattern}{spacing_pattern}){{{nested - 1}}}{int_pattern}{spacing_pattern}?\))"
+            int_element = (
+                Opt(Literal(str(nested))).suppress()
+                + Literal("(").suppress()
+                + Group(common.integer[nested])
+                + Literal(")").suppress()
+            )
+        float_element_pattern = rf"(?:(?:{nested})?{spacing_pattern}?\({spacing_pattern}?(?:{float_pattern}{spacing_pattern}){{{nested - 1}}}{float_pattern}{spacing_pattern}?\))"
+        float_element = (
+            Opt(Literal(str(nested))).suppress()
+            + Literal("(").suppress()
+            + Group(common.ieee_float[nested])
+            + Literal(")").suppress()
+        )
 
-    ignore_pattern = rf"(?:\s|{ignore.re.pattern})+" if ignore is not None else r"\s+"
+    if not force_float:
+        int_list = Forward()
+    float_list = Forward()
 
-    ret = Regex(
-        rf"\((?:{ignore_pattern})?(?:{float_pattern}{ignore_pattern}){{{tensor_kind.size - 1}}}{float_pattern}(?:{ignore_pattern})?\)"
-    )
-    ret.add_parse_action(
-        lambda tks: np.fromstring(
-            re.sub(ignore.re, " ", tks[0][1:-1])
-            if ignore is not None
-            else tks[0][1:-1],
-            sep=" ",
-        )
-    )
-    return ret
+    def process_count(tks: ParseResults) -> None:
+        nonlocal int_list, float_list
 
+        if not tks:
+            count = None
+        else:
+            (count,) = tks
+            assert isinstance(count, int)
 
-def _list_of(entry: ParserElement) -> ParserElement:
-    return Opt(
-        Literal("List") + Literal("<") + _IDENTIFIER + Literal(">")
-    ).suppress() + (
-        (
-            counted_array(entry, common.integer + Literal("(").suppress())
-            + Literal(")").suppress()
-        ).set_parse_action(lambda tks: [tks.as_list()])
-        | (
-            Literal("(").suppress()
-            + Group((entry)[...], aslist=True)
-            + Literal(")").suppress()
+        if count is None:
+            if not force_float:
+                int_list_pattern = rf"\({spacing_pattern}?(?:{int_element_pattern}{spacing_pattern})*{int_element_pattern}{spacing_pattern}?\)"
+                float_list_pattern = rf"\({spacing_pattern}?(?:{float_element_pattern}{spacing_pattern})*{float_element_pattern}{spacing_pattern}?\)"
+            else:
+                float_list_pattern = rf"\({spacing_pattern}?(?:{float_element_pattern}{spacing_pattern})*{float_element_pattern}?{spacing_pattern}?\)"
+
+        elif count == 0:
+            if not force_float:
+                int_list <<= NoMatch()
+                float_list <<= NoMatch()
+            else:
+                float_list <<= (Literal("(") + Literal(")")).add_parse_action(
+                    lambda: np.empty((0, nested) if nested else 0, dtype=float)
+                )
+            return
+
+        else:
+            if not force_float:
+                int_list_pattern = rf"\({spacing_pattern}?(?:{int_element_pattern}{spacing_pattern}){{{count - 1}}}{int_element_pattern}{spacing_pattern}?\)"
+            float_list_pattern = rf"\({spacing_pattern}?(?:{float_element_pattern}{spacing_pattern}){{{count - 1}}}{float_element_pattern}{spacing_pattern}?\)"
+
+        if not force_float:
+            int_list <<= Regex(int_list_pattern).add_parse_action(
+                lambda tks: to_array(tks, dtype=int)
+            )
+        float_list <<= Regex(float_list_pattern).add_parse_action(
+            lambda tks: to_array(tks, dtype=float)
         )
-        | (
-            common.integer + Literal("{").suppress() + entry + Literal("}").suppress()
-        ).set_parse_action(lambda tks: [[tks[1]] * tks[0]])
-    )
 
+    def to_array(
+        tks: ParseResults, *, dtype: type
+    ) -> np.ndarray[tuple[int] | tuple[int, int], np.dtype[np.int64 | np.float64]]:
+        (s,) = tks
+        s = s.replace("(", "").replace(")", "")
 
-def _parse_ascii_field(
-    s: str, tensor_kind: TensorKind, *, ignore: Regex | None
-) -> np.ndarray[tuple[int] | tuple[int, int], np.dtype[np.float64]]:
-    if ignore is not None:
-        s = re.sub(ignore.re, " ", s)
-    s = s.replace("(", " ").replace(")", " ")
+        if ignore is not None:
+            s = re.sub(ignore.re, " ", s)
 
-    return np.fromstring(s, sep=" ").reshape(-1, *tensor_kind.shape)  # type: ignore [return-value]
+        ret: np.ndarray[
+            tuple[int] | tuple[int, int], np.dtype[np.int64 | np.float64]
+        ] = np.fromstring(s, sep=" ", dtype=dtype)  # type: ignore[assignment]
 
+        if nested is not None:
+            ret = ret.reshape(-1, nested)
 
-def _unpack_binary_field(
-    b: bytes, tensor_kind: TensorKind, *, length: int
-) -> np.ndarray[tuple[int] | tuple[int, int], np.dtype[np.float64 | np.float32]]:
-    float_size = len(b) / tensor_kind.size / length
-    assert float_size in (4, 8)
+        return ret
 
-    dtype = np.float32 if float_size == 4 else float
-    return np.frombuffer(b, dtype=dtype).reshape(-1, *tensor_kind.shape)  # type: ignore [return-value]
+    def to_full_array(
+        tks: ParseResults, *, dtype: type
+    ) -> np.ndarray[tuple[int] | tuple[int, int], np.dtype[np.int64 | np.float64]]:
+        count, lst = tks
+        assert isinstance(count, int)
 
+        if nested is None:
+            return np.full(count, lst, dtype=dtype)
 
-def _tensor_list(
-    tensor_kind: TensorKind, *, ignore: Regex | None = None
-) -> ParserElement:
-    tensor = _tensor(tensor_kind, ignore=ignore)
-    ignore_pattern = rf"(?:\s|{ignore.re.pattern})+" if ignore is not None else r"\s+"
-
-    list_ = Forward()
-
-    list_ <<= Regex(
-        rf"\((?:{ignore_pattern})?(?:{tensor.re.pattern}{ignore_pattern})*{tensor.re.pattern}(?:{ignore_pattern})?\)"
-    ).add_parse_action(
-        lambda tks: [_parse_ascii_field(tks[0], tensor_kind, ignore=ignore)]
-    ) | (
-        (Literal("(") + Literal(")"))
-        .suppress()
-        .add_parse_action(lambda: np.array([]).reshape(0, *tensor_kind.shape))
+        return np.full((count, nested), lst, dtype=dtype)  # type: ignore[return-value]
+
+    count = Opt(common.integer).add_parse_action(process_count)
+
+    ret = count.suppress() + (
+        (int_list | float_list) if not force_float else float_list
     )
 
-    def count_parse_action(tks: ParseResults) -> None:
-        nonlocal list_
-        length = tks[0]
-        assert isinstance(length, int)
+    if not force_float:
+        ret |= (
+            common.integer
+            + Literal("{").suppress()
+            + int_element
+            + Literal("}").suppress()
+        ).add_parse_action(lambda tks: to_full_array(tks, dtype=int))
+    ret |= (
+        common.integer
+        + Literal("{").suppress()
+        + float_element
+        + Literal("}").suppress()
+    ).add_parse_action(lambda tks: to_full_array(tks, dtype=float))
 
-        if not length:
-            list_ <<= (
-                ((Literal("(") + Literal(")")) | (Literal("{") + Literal("}")))
-                .suppress()
-                .add_parse_action(lambda: np.array([]).reshape(0, *tensor_kind.shape))
-            )
-            return
+    if ignore is not None:
+        ret.ignore(ignore)
 
-        list_ <<= (
-            Regex(
-                rf"\((?:{ignore_pattern})?(?:{tensor.re.pattern}{ignore_pattern}){{{length - 1}}}{tensor.re.pattern}(?:{ignore_pattern})?\)"
-            ).add_parse_action(
-                lambda tks: [_parse_ascii_field(tks[0], tensor_kind, ignore=ignore)]
-            )
-            | Regex(
-                rf"\((?s:({'.' * 8 * tensor_kind.size}|{'.' * 4 * tensor_kind.size}){{{length}}})\)"
-            ).add_parse_action(
-                lambda tks: [
-                    _unpack_binary_field(
-                        tks[0][1:-1].encode("latin-1"), tensor_kind, length=length
-                    )
-                ]
-            )
-            | (
-                Literal("{").suppress() + tensor + Literal("}").suppress()
-            ).add_parse_action(
-                lambda tks: [np.full((length, *tensor_kind.shape), tks[0], dtype=float)]
-            )
+    return ret
+
+
+def _binary_field(*, nested: int | None = None) -> ParserElement:
+    elsize = nested if nested is not None else 1
+
+    binary_field = Forward()
+
+    def process_count(tks: ParseResults) -> None:
+        nonlocal binary_field
+        (size,) = tks
+        assert isinstance(size, int)
+
+        binary_field <<= Regex(
+            rf"\((?s:({'.' * 8 * elsize}|{'.' * 4 * elsize}){{{size}}})\)"
         )
 
-    count = common.integer.copy().add_parse_action(count_parse_action)
+    def to_array(
+        tks: ParseResults,
+    ) -> np.ndarray[tuple[int] | tuple[int, int], np.dtype[np.float64 | np.float32]]:
+        size, s = tks
+        assert isinstance(size, int)
+        assert isinstance(s, str)
+        assert s[0] == "("
+        assert s[-1] == ")"
+        s = s[1:-1]
+
+        float_size = len(s) / elsize / size
+        assert float_size in (4, 8)
+
+        dtype = np.float32 if float_size == 4 else float
+        ret = np.frombuffer(s.encode("latin-1"), dtype=dtype)
 
+        if nested is not None:
+            ret = ret.reshape(-1, nested)
+
+        return ret  # type: ignore[return-value]
+
+    count = common.integer.copy().add_parse_action(process_count)
+
+    return (count + binary_field).add_parse_action(to_array)
+
+
+def _list_of(entry: ParserElement) -> ParserElement:
     return (
-        Opt(Literal("List") + Literal("<") + str(tensor_kind) + Literal(">")).suppress()
-        + Opt(count).suppress()
-        + list_
+        (
+            counted_array(entry, common.integer + Literal("(").suppress())
+            + Literal(")").suppress()
+        ).set_parse_action(lambda tks: [tks.as_list()])
+        | (
+            Literal("(").suppress()
+            + Group((entry)[...], aslist=True)
+            + Literal(")").suppress()
+        )
+        | (
+            common.integer + Literal("{").suppress() + entry + Literal("}").suppress()
+        ).set_parse_action(lambda tks: [[tks[1]] * tks[0]])
     )
 
 
@@ -250,12 +314,11 @@ _SWITCH = (
 _DIMENSIONS = (
     Literal("[").suppress() + common.number[0, 7] + Literal("]").suppress()
 ).set_parse_action(lambda tks: DimensionSet(*tks))
-_TENSOR = (
-    _tensor(TensorKind.SCALAR)
-    | _tensor(TensorKind.VECTOR)
-    | _tensor(TensorKind.SYMM_TENSOR)
-    | _tensor(TensorKind.TENSOR)
-)
+_TENSOR = common.ieee_float | (
+    Literal("(").suppress()
+    + Group(common.ieee_float[3] | common.ieee_float[6] | common.ieee_float[9])
+    + Literal(")").suppress()
+).add_parse_action(lambda tks: np.array(tks[0], dtype=float))
 _PARENTHESIZED = Forward()
 _IDENTIFIER = Combine(Word(_IDENTCHARS, _IDENTBODYCHARS) + Opt(_PARENTHESIZED))
 _PARENTHESIZED <<= Combine(
@@ -270,12 +333,46 @@ _DIMENSIONED = (Opt(_IDENTIFIER) + _DIMENSIONS + _TENSOR).set_parse_action(
 _FIELD = (Keyword("uniform", _IDENTBODYCHARS).suppress() + _TENSOR) | (
     Keyword("nonuniform", _IDENTBODYCHARS).suppress()
     + (
-        _tensor_list(TensorKind.SCALAR, ignore=_COMMENT)
-        | _tensor_list(TensorKind.VECTOR, ignore=_COMMENT)
-        | _tensor_list(TensorKind.SYMM_TENSOR, ignore=_COMMENT)
-        | _tensor_list(TensorKind.TENSOR, ignore=_COMMENT)
+        (
+            Opt(
+                Literal("List") + Literal("<") + Literal("scalar") + Literal(">")
+            ).suppress()
+            + (_numeric_list(force_float=True, ignore=_COMMENT) | _binary_field())
+        )
+        | (
+            Opt(
+                Literal("List") + Literal("<") + Literal("vector") + Literal(">")
+            ).suppress()
+            + (
+                _numeric_list(nested=3, force_float=True, ignore=_COMMENT)
+                | _binary_field(nested=3)
+            )
+        )
+        | (
+            Opt(
+                Literal("List") + Literal("<") + Literal("symmTensor") + Literal(">")
+            ).suppress()
+            + (
+                _numeric_list(
+                    nested=6,
+                    force_float=True,
+                    ignore=_COMMENT,
+                )
+                | _binary_field(nested=6)
+            )
+        )
+        | (
+            Opt(
+                Literal("List") + Literal("<") + Literal("tensor") + Literal(">")
+            ).suppress()
+            + (
+                _numeric_list(nested=9, force_float=True, ignore=_COMMENT)
+                | _binary_field(nested=9)
+            )
+        )
     )
 )
+
 _DIRECTIVE = Word("#", _IDENTBODYCHARS)
 _TOKEN = dbl_quoted_string | _DIRECTIVE | _IDENTIFIER
 _DATA = Forward()
@@ -283,7 +380,12 @@ _KEYWORD_ENTRY = _keyword_entry_of(_TOKEN | _list_of(_IDENTIFIER), _DATA)
 _DICT = _dict_of(_TOKEN, _DATA)
 _DATA_ENTRY = Forward()
 _LIST_ENTRY = _DICT | _KEYWORD_ENTRY | _DATA_ENTRY
-_LIST = _list_of(_LIST_ENTRY)
+_LIST = (
+    _numeric_list(ignore=_COMMENT)
+    | _numeric_list(nested=3, ignore=_COMMENT)
+    | _numeric_list(nested=4, ignore=_COMMENT)
+    | _list_of(_LIST_ENTRY)
+)
 _NUMBER = (
     common.number
     | CaselessKeyword("nan").set_parse_action(lambda: np.nan)


Test patch:
diff --git a/tests/test_files/test_files.py b/tests/test_files/test_files.py
index 1f8e405..8c9c19a 100644
--- a/tests/test_files/test_files.py
+++ b/tests/test_files/test_files.py
@@ -3,9 +3,9 @@ import sys
 from pathlib import Path
 
 if sys.version_info >= (3, 9):
-    from collections.abc import Generator, Sequence
+    from collections.abc import Generator
 else:
-    from typing import Generator, Sequence
+    from typing import Generator
 
 import numpy as np
 import pytest
@@ -153,16 +153,16 @@ def test_mesh(cavity: FoamCase) -> None:
 
     points = file[None]
 
-    assert isinstance(points, Sequence)
-    assert isinstance(points[0], Sequence)
-    assert len(points[0]) == 3
+    assert isinstance(points, np.ndarray)
+    assert points.ndim == 2
+    assert points.shape[-1] == 3
 
 
 def test_internal_field(cavity: FoamCase) -> None:
     blocks = cavity.block_mesh_dict["blocks"]
     assert isinstance(blocks, list)
     sizes = blocks[2]
-    assert isinstance(sizes, list)
+    assert isinstance(sizes, np.ndarray)
     size = np.prod(sizes)
 
     p_arr = np.zeros(size)
diff --git a/tests/test_files/test_parsing/test_advanced.py b/tests/test_files/test_parsing/test_advanced.py
index 500026e..d9fdc65 100644
--- a/tests/test_files/test_parsing/test_advanced.py
+++ b/tests/test_files/test_parsing/test_advanced.py
@@ -1,5 +1,6 @@
 # Based on https://foss.heptapod.net/fluiddyn/fluidsimfoam/-/blob/branch/default/tests/test_parser_advanced.py
 
+import numpy as np
 import pytest
 from foamlib._files._parsing import Parsed
 
@@ -276,9 +277,10 @@ def test_assignment_strange_name() -> None:
     """
     )
     assert parsed[("equations", '"(U|e|k).*"')] == 0.7
-    assert parsed[("equations", '"(U|e|k|epsilon).*"')] == (
-        "table",
-        [[0, 0.4], [0.5, 0.7]],
+    assert isinstance(parsed[("equations", '"(U|e|k|epsilon).*"')], tuple)
+    assert parsed[("equations", '"(U|e|k|epsilon).*"')][0] == "table"
+    assert np.array_equal(
+        parsed[("equations", '"(U|e|k|epsilon).*"')][1], [[0, 0.4], [0.5, 0.7]]
     )
 
 
@@ -331,17 +333,23 @@ def test_list_u() -> None:
         )
         """
     )
-    assert parsed[()] == [
-        [4.50773, 1.79963, 0.0],
-        [6.06208, 2.40831, 0.0],
-        [6.874, 2.72079, 0.0],
-        [7.42929, 2.931, 0.0],
-        [7.85095, 3.08805, 0.0],
-        [8.19202, 3.21306, 0.0],
-        [17.5, 1.92559e-09, 0.0],
-        [17.5, 6.81045e-12, 0.0],
-        [17.5, 6.81045e-12, 0.0],
-    ]
+    data = parsed[()]
+    assert isinstance(data, np.ndarray)
+    assert data.shape == (9, 3)
+    assert np.array_equal(
+        data,
+        [
+            [4.507730000e00, 1.799630000e00, 0.000000000e00],
+            [6.062080000e00, 2.408310000e00, 0.000000000e00],
+            [6.874000000e00, 2.720790000e00, 0.000000000e00],
+            [7.429290000e00, 2.931000000e00, 0.000000000e00],
+            [7.850950000e00, 3.088050000e00, 0.000000000e00],
+            [8.192020000e00, 3.213060000e00, 0.000000000e00],
+            [1.750000000e01, 1.925590000e-09, 0.000000000e00],
+            [1.750000000e01, 6.810450000e-12, 0.000000000e00],
+            [1.750000000e01, 6.810450000e-12, 0.000000000e00],
+        ],
+    )
 
 
 def test_list_as_write_cell_centers() -> None:
diff --git a/tests/test_files/test_parsing/test_basic.py b/tests/test_files/test_parsing/test_basic.py
index 3155608..52302e9 100644
--- a/tests/test_files/test_parsing/test_basic.py
+++ b/tests/test_files/test_parsing/test_basic.py
@@ -16,7 +16,16 @@ def test_parse_value() -> None:
     assert Parsed(b"uniform 1")[()] == 1
     assert Parsed(b"uniform 1.0")[()] == 1.0
     assert Parsed(b"uniform 1.0e-3")[()] == 1.0e-3
-    assert Parsed(b"(1.0 2.0 3.0)")[()] == [1.0, 2.0, 3.0]
+    assert Parsed(b"(word word)")[()] == ["word", "word"]
+    lst = Parsed(b"(1 2 3)")[()]
+    assert isinstance(lst, np.ndarray)
+    assert lst.dtype == np.int64
+    assert np.array_equal(lst, [1, 2, 3])
+    lst = Parsed(b"(1.0 2 3)")[()]
+    assert isinstance(lst, np.ndarray)
+    assert lst.dtype == np.float64
+    assert np.array_equal(lst, [1.0, 2.0, 3.0])
+    assert Parsed(b"()")[()] == []
     field = Parsed(b"uniform (1 2 3)")[()]
     assert isinstance(field, np.ndarray)
     assert np.array_equal(field, [1, 2, 3])
@@ -32,15 +41,28 @@ def test_parse_value() -> None:
     field = Parsed(b"nonuniform List<tensor> ()")[()]
     assert isinstance(field, np.ndarray)
     assert field.shape == (0, 9)
-    assert Parsed(b"3(1 2 3)")[()] == [1, 2, 3]
-    assert Parsed(b"2((1 2 3) (4 5 6))")[()] == [
-        [1, 2, 3],
-        [4, 5, 6],
-    ]
-    assert Parsed(b"2{(1 2 3)}")[()] == [
-        [1, 2, 3],
-        [1, 2, 3],
-    ]
+    lst = Parsed(b"3(1 2 3)")[()]
+    assert isinstance(lst, np.ndarray)
+    assert np.array_equal(lst, [1, 2, 3])
+    lst = Parsed(b"2((1 2 3) (4 5 6))")[()]
+    assert isinstance(lst, np.ndarray)
+    assert np.array_equal(
+        lst,
+        [
+            [1, 2, 3],
+            [4, 5, 6],
+        ],
+    )
+    lst = Parsed(b"2{(1 2 3)}")[()]
+    assert isinstance(lst, np.ndarray)
+    assert np.array_equal(
+        lst,
+        [
+            [1, 2, 3],
+            [1, 2, 3],
+        ],
+    )
+    assert Parsed(b"0()")[()] == []
     field = Parsed(b"nonuniform List<vector> 2((1 2 3) (4 5 6))")[()]
     assert isinstance(field, np.ndarray)
     assert np.array_equal(
@@ -83,19 +105,26 @@ def test_parse_value() -> None:
     assert Parsed(b"[1 1 -2 0 0 0 0] 9.81")[()] == FoamFile.Dimensioned(
         dimensions=FoamFile.DimensionSet(mass=1, length=1, time=-2), value=9.81
     )
-    assert Parsed(b"hex (0 1 2 3 4 5 6 7) (1 1 1) simpleGrading (1 1 1)")[()] == (
-        "hex",
-        [0, 1, 2, 3, 4, 5, 6, 7],
-        [1, 1, 1],
-        "simpleGrading",
-        [1, 1, 1],
-    )
+    tpl = Parsed(b"hex (0 1 2 3 4 5 6 7) (1 1 1) simpleGrading (1 1 1)")[()]
+    assert isinstance(tpl, tuple)
+    assert len(tpl) == 5
+    assert tpl[0] == "hex"
+    assert isinstance(tpl[1], np.ndarray)
+    assert tpl[1].dtype == np.int64
+    assert np.array_equal(tpl[1], [0, 1, 2, 3, 4, 5, 6, 7])
+    assert isinstance(tpl[2], np.ndarray)
+    assert tpl[2].dtype == np.int64
+    assert np.array_equal(tpl[2], [1, 1, 1])
+    assert tpl[3] == "simpleGrading"
+    assert isinstance(tpl[4], np.ndarray)
+    assert tpl[4].dtype == np.int64
+    assert np.array_equal(tpl[4], [1, 1, 1])
     assert Parsed(b"(a b; c d;)")[()] == [("a", "b"), ("c", "d")]
     assert Parsed(b"(a {b c;} d {e g;})")[()] == [
         ("a", {"b": "c"}),
         ("d", {"e": "g"}),
     ]
-    assert Parsed(b"(a (0 1 2); b {})")[()] == [("a", [0, 1, 2]), ("b", {})]
+    assert Parsed(b"(a (b c d); e {})")[()] == [("a", ["b", "c", "d"]), ("e", {})]
     assert Parsed(b"({a b; c d;} {e g;})")[()] == [{"a": "b", "c": "d"}, {"e": "g"}]
     assert Parsed(b"(water oil mercury air)")[()] == ["water", "oil", "mercury", "air"]
     assert Parsed(b"div(phi,U)")[()] == "div(phi,U)"
diff --git a/tests/test_files/test_parsing/test_decompose_par.py b/tests/test_files/test_parsing/test_decompose_par.py
index e3515b0..07497a8 100644
--- a/tests/test_files/test_parsing/test_decompose_par.py
+++ b/tests/test_files/test_parsing/test_decompose_par.py
@@ -3,6 +3,7 @@
 from pathlib import Path
 from textwrap import dedent
 
+import numpy as np
 from foamlib import FoamFile
 
 
@@ -43,7 +44,8 @@ def test_simple(tmp_path: Path) -> None:
 
     assert decompose_par_dict["numberOfSubdomains"] == 8
     assert decompose_par_dict["method"] == "simple"
-    assert decompose_par_dict["coeffs", "n"] == [4, 2, 1]
+    assert isinstance(decompose_par_dict["coeffs", "n"], np.ndarray)
+    assert np.array_equal(decompose_par_dict["coeffs", "n"], [4, 2, 1])
     assert decompose_par_dict["coeffs", "order"] == "xyz"
     assert decompose_par_dict["coeffs", "delta"] == 0.001
 
diff --git a/tests/test_files/test_parsing/test_intermediate.py b/tests/test_files/test_parsing/test_intermediate.py
index 6fd7819..3057621 100644
--- a/tests/test_files/test_parsing/test_intermediate.py
+++ b/tests/test_files/test_parsing/test_intermediate.py
@@ -1,5 +1,6 @@
 # Based on https://foss.heptapod.net/fluiddyn/fluidsimfoam/-/blob/branch/default/tests/test_parser.py
 
+import numpy as np
 import pytest
 from foamlib import FoamFile
 from foamlib._files._parsing import Parsed
@@ -61,20 +62,20 @@ def test_strange_names() -> None:
 
 
 def test_list_simple() -> None:
-    assert Parsed(b"""
+    faces = Parsed(b"""
         faces
         (
             (1 5 4 0)
             (2 3 4 5)
         );
-    """)[("faces",)] == [
-        [1, 5, 4, 0],
-        [2, 3, 4, 5],
-    ]
+    """)[("faces",)]
+
+    assert isinstance(faces, np.ndarray)
+    assert np.array_equal(faces, [[1, 5, 4, 0], [2, 3, 4, 5]])
 
 
 def test_list_assignment() -> None:
-    assert Parsed(b"""
+    faces = Parsed(b"""
         faces
         (
             1
@@ -82,7 +83,9 @@ def test_list_assignment() -> None:
             4
             0
         );
-    """)[("faces",)] == [1, 5, 4, 0]
+    """)[("faces",)]
+    assert isinstance(faces, np.ndarray)
+    assert np.array_equal(faces, [1, 5, 4, 0])
 
 
 def test_dict_simple() -> None:
@@ -135,11 +138,12 @@ def test_dict_with_list() -> None:
             pRefValue                   0;
         }
     """)
-    assert parsed[("PISO", "pRefPoint")] == [0, 0, 0]
+    assert isinstance(parsed[("PISO", "pRefPoint")], np.ndarray)
+    assert np.array_equal(parsed[("PISO", "pRefPoint")], [0, 0, 0])
 
 
 def test_list_with_dict() -> None:
-    assert Parsed(b"""
+    boundary = Parsed(b"""
         boundary
         (
             upperBoundary
@@ -152,33 +156,36 @@ def test_list_with_dict() -> None:
                 );
             }
         );
-    """)[("boundary",)] == [
-        (
-            "upperBoundary",
-            {
-                "type": "cyclic",
-                "neighbourPatch": "lowerBoundary",
-                "faces": [
-                    [3, 7, 6, 2],
-                ],
-            },
-        ),
-    ]
+    """)[("boundary",)]
+    assert isinstance(boundary, list)
+    assert len(boundary) == 1
+
+    upper_boundary = boundary[0]
+    assert upper_boundary[0] == "upperBoundary"
+    assert upper_boundary[1]["type"] == "cyclic"
+    assert upper_boundary[1]["neighbourPatch"] == "lowerBoundary"
+    assert isinstance(upper_boundary[1]["faces"], np.ndarray)
+    assert np.array_equal(upper_boundary[1]["faces"], [[3, 7, 6, 2]])
 
 
 def test_list_with_str() -> None:
-    assert Parsed(b"""
+    blocks = Parsed(b"""
         blocks
         (
             hex (0 1 2 3 4 5 6 7) (40 40 40) simpleGrading (1 1 1)
         );
-    """)[("blocks",)] == [
-        "hex",
-        [0, 1, 2, 3, 4, 5, 6, 7],
-        [40, 40, 40],
-        "simpleGrading",
-        [1, 1, 1],
-    ]
+    """)[("blocks",)]
+    assert isinstance(blocks, list)
+    assert len(blocks) == 5
+
+    assert blocks[0] == "hex"
+    assert isinstance(blocks[1], np.ndarray)
+    assert np.array_equal(blocks[1], [0, 1, 2, 3, 4, 5, 6, 7])
+    assert isinstance(blocks[2], np.ndarray)
+    assert np.array_equal(blocks[2], [40, 40, 40])
+    assert blocks[3] == "simpleGrading"
+    assert isinstance(blocks[4], np.ndarray)
+    assert np.array_equal(blocks[4], [1, 1, 1])
 
 
 def test_file_simple() -> None:
@@ -230,10 +237,9 @@ def test_file() -> None:
     """)
     assert parsed[("a",)] == 1
     assert parsed[("b",)] == 2
-    assert parsed[("faces",)] == [
-        [1, 5, 4, 0],
-        [2, 3, 4, 5],
-    ]
+    faces = parsed[("faces",)]
+    assert isinstance(faces, np.ndarray)
+    assert np.array_equal(faces, [[1, 5, 4, 0], [2, 3, 4, 5]])
     assert parsed[("my_dict", "a")] == 1
 
 
@@ -455,12 +461,15 @@ def test_for_blockmesh() -> None:
         );
     """)
     assert parsed[("negHalfWidth",)] == ("#neg", "$halfWidth")
-    assert parsed[("blocks",)] == [
-        "hex",
-        [4, 6, 14, 12, 0, 2, 10, 8],
-        [1, "$upstreamCells", "$cylinderBoxCells"],
-        "$expandBlock",
-    ]
+    blocks = parsed[("blocks",)]
+    assert isinstance(blocks, list)
+    assert len(blocks) == 4
+
+    assert blocks[0] == "hex"
+    assert isinstance(blocks[1], np.ndarray)
+    assert np.array_equal(blocks[1], [4, 6, 14, 12, 0, 2, 10, 8])
+    assert blocks[2] == [1, "$upstreamCells", "$cylinderBoxCells"]
+    assert blocks[3] == "$expandBlock"
 
 
 def test_for_u() -> None:
@@ -483,30 +492,44 @@ def test_blocks() -> None:
             hex (16 17 18 19 20 21 22 23) (96 1 72) simpleGrading (1 1 1)
         );
     """)
-    assert parsed[("blocks",)] == [
-        "hex",
-        [0, 1, 2, 3, 4, 5, 6, 7],
-        "inletChannel",
-        [40, 1, 64],
-        "simpleGrading",
-        [1, 1, 1],
-        "hex",
-        [4, 5, 6, 7, 8, 9, 10, 11, 12],
-        "inletChannel",
-        [40, 1, 16],
-        "simpleGrading",
-        [1, 1, 1],
-        "hex",
-        [12, 13, 14, 15, 16, 17, 18, 19],
-        [96, 1, 8],
-        "simpleGrading",
-        [1, 1, 1],
-        "hex",
-        [16, 17, 18, 19, 20, 21, 22, 23],
-        [96, 1, 72],
-        "simpleGrading",
-        [1, 1, 1],
-    ]
+    blocks = parsed[("blocks",)]
+    assert isinstance(blocks, list)
+    assert len(blocks) == 22
+
+    assert blocks[0] == "hex"
+    assert isinstance(blocks[1], np.ndarray)
+    assert np.array_equal(blocks[1], [0, 1, 2, 3, 4, 5, 6, 7])
+    assert blocks[2] == "inletChannel"
+    assert isinstance(blocks[3], np.ndarray)
+    assert np.array_equal(blocks[3], [40, 1, 64])
+    assert blocks[4] == "simpleGrading"
+    assert isinstance(blocks[5], np.ndarray)
+    assert np.array_equal(blocks[5], [1, 1, 1])
+    assert blocks[6] == "hex"
+    assert isinstance(blocks[7], np.ndarray)
+    assert np.array_equal(blocks[7], [4, 5, 6, 7, 8, 9, 10, 11, 12])
+    assert blocks[8] == "inletChannel"
+    assert isinstance(blocks[9], np.ndarray)
+    assert np.array_equal(blocks[9], [40, 1, 16])
+    assert blocks[10] == "simpleGrading"
+    assert isinstance(blocks[11], np.ndarray)
+    assert np.array_equal(blocks[11], [1, 1, 1])
+    assert blocks[12] == "hex"
+    assert isinstance(blocks[13], np.ndarray)
+    assert np.array_equal(blocks[13], [12, 13, 14, 15, 16, 17, 18, 19])
+    assert isinstance(blocks[14], np.ndarray)
+    assert np.array_equal(blocks[14], [96, 1, 8])
+    assert blocks[15] == "simpleGrading"
+    assert isinstance(blocks[16], np.ndarray)
+    assert np.array_equal(blocks[16], [1, 1, 1])
+    assert blocks[17] == "hex"
+    assert isinstance(blocks[18], np.ndarray)
+    assert np.array_equal(blocks[18], [16, 17, 18, 19, 20, 21, 22, 23])
+    assert isinstance(blocks[19], np.ndarray)
+    assert np.array_equal(blocks[19], [96, 1, 72])
+    assert blocks[20] == "simpleGrading"
+    assert isinstance(blocks[21], np.ndarray)
+    assert np.array_equal(blocks[21], [1, 1, 1])
 
 
 @pytest.mark.xfail(reason="Not currently supported")
@@ -586,10 +609,16 @@ def test_list_edges() -> None:
             spline 6 5 ((0.6 0.0124 0.05) (0.7 0.0395 0.05) (0.8 0.0724 0.05) (0.9 0.132 0.05) (1 0.172 0.05) (1.1 0.132 0.05) (1.2 0.0724 0.05) (1.3 0.0395 0.05) (1.4 0.0124 0.05))
         );
     """)
-    assert parsed[("edges",)] == [
-        "spline",
-        1,
-        2,
+    edges = parsed[("edges",)]
+    assert isinstance(edges, list)
+    assert len(edges) == 8
+
+    assert edges[0] == "spline"
+    assert edges[1] == 1
+    assert edges[2] == 2
+    assert isinstance(edges[3], np.ndarray)
+    assert np.array_equal(
+        edges[3],
         [
             [0.6, 0.0124, 0.0],
             [0.7, 0.0395, 0.0],
@@ -601,9 +630,13 @@ def test_list_edges() -> None:
             [1.3, 0.0395, 0.0],
             [1.4, 0.0124, 0.0],
         ],
-        "spline",
-        6,
-        5,
+    )
+    assert edges[4] == "spline"
+    assert edges[5] == 6
+    assert edges[6] == 5
+    assert isinstance(edges[7], np.ndarray)
+    assert np.array_equal(
+        edges[7],
         [
             [0.6, 0.0124, 0.05],
             [0.7, 0.0395, 0.05],
@@ -615,7 +648,7 @@ def test_list_edges() -> None:
             [1.3, 0.0395, 0.05],
             [1.4, 0.0124, 0.05],
         ],
-    ]
+    )
 
 
 def test_list_edges_arcs() -> None:
@@ -626,18 +659,21 @@ def test_list_edges_arcs() -> None:
             arc 5 10 origin (0 0 0)
         );
     """)
-    assert parsed[("edges",)] == [
-        "arc",
-        0,
-        5,
-        "origin",
-        [0, 0, 0],
-        "arc",
-        5,
-        10,
-        "origin",
-        [0, 0, 0],
-    ]
+    edges = parsed[("edges",)]
+    assert isinstance(edges, list)
+    assert len(edges) == 10
+    assert edges[0] == "arc"
+    assert edges[1] == 0
+    assert edges[2] == 5
+    assert edges[3] == "origin"
+    assert isinstance(edges[4], np.ndarray)
+    assert np.array_equal(edges[4], [0, 0, 0])
+    assert edges[5] == "arc"
+    assert edges[6] == 5
+    assert edges[7] == 10
+    assert edges[8] == "origin"
+    assert isinstance(edges[9], np.ndarray)
+    assert np.array_equal(edges[9], [0, 0, 0])
 
 
 def test_list_blocks() -> None:
@@ -649,23 +685,47 @@ def test_list_blocks() -> None:
             hex (2 3 11 10 5 4 12 13) (225 100 1) simpleGrading (1 ((0.1 0.25 41.9) (0.9 0.75 1)) 1)
         );
     """)
-    assert parsed[("blocks",)] == [
-        "hex",
-        [0, 1, 9, 8, 7, 6, 14, 15],
-        [50, 100, 1],
-        "simpleGrading",
-        [1, [[0.1, 0.25, 41.9], [0.9, 0.75, 1]], 1],
-        "hex",
-        [1, 2, 10, 9, 6, 5, 13, 14],
-        [50, 100, 1],
-        "simpleGrading",
-        [1, [[0.1, 0.25, 41.9], [0.9, 0.75, 1]], 1],
-        "hex",
-        [2, 3, 11, 10, 5, 4, 12, 13],
-        [225, 100, 1],
-        "simpleGrading",
-        [1, [[0.1, 0.25, 41.9], [0.9, 0.75, 1]], 1],
-    ]
+    blocks = parsed[("blocks",)]
+
+    assert isinstance(blocks, list)
+    assert len(blocks) == 15
+
+    assert blocks[0] == "hex"
+    assert isinstance(blocks[1], np.ndarray)
+    assert np.array_equal(blocks[1], [0, 1, 9, 8, 7, 6, 14, 15])
+    assert isinstance(blocks[2], np.ndarray)
+    assert np.array_equal(blocks[2], [50, 100, 1])
+    assert blocks[3] == "simpleGrading"
+    assert isinstance(blocks[4], list)
+    assert len(blocks[4]) == 3
+    assert blocks[4][0] == 1
+    assert isinstance(blocks[4][1], np.ndarray)
+    assert np.array_equal(blocks[4][1], [[0.1, 0.25, 41.9], [0.9, 0.75, 1]])
+    assert blocks[4][2] == 1
+    assert blocks[5] == "hex"
+    assert isinstance(blocks[6], np.ndarray)
+    assert np.array_equal(blocks[6], [1, 2, 10, 9, 6, 5, 13, 14])
+    assert isinstance(blocks[7], np.ndarray)
+    assert np.array_equal(blocks[7], [50, 100, 1])
+    assert blocks[8] == "simpleGrading"
+    assert isinstance(blocks[9], list)
+    assert len(blocks[9]) == 3
+    assert blocks[9][0] == 1
+    assert isinstance(blocks[9][1], np.ndarray)
+    assert np.array_equal(blocks[9][1], [[0.1, 0.25, 41.9], [0.9, 0.75, 1]])
+    assert blocks[9][2] == 1
+    assert blocks[10] == "hex"
+    assert isinstance(blocks[11], np.ndarray)
+    assert np.array_equal(blocks[11], [2, 3, 11, 10, 5, 4, 12, 13])
+    assert isinstance(blocks[12], np.ndarray)
+    assert np.array_equal(blocks[12], [225, 100, 1])
+    assert blocks[13] == "simpleGrading"
+    assert isinstance(blocks[14], list)
+    assert len(blocks[14]) == 3
+    assert blocks[14][0] == 1
+    assert isinstance(blocks[14][1], np.ndarray)
+    assert np.array_equal(blocks[14][1], [[0.1, 0.25, 41.9], [0.9, 0.75, 1]])
+    assert blocks[14][2] == 1
 
 
 @pytest.mark.xfail(reason="Not currently supported")
diff --git a/tests/test_files/test_parsing/test_poly_mesh.py b/tests/test_files/test_parsing/test_poly_mesh.py
index 7c32ff2..f2b4016 100644
--- a/tests/test_files/test_parsing/test_poly_mesh.py
+++ b/tests/test_files/test_parsing/test_poly_mesh.py
@@ -2,6 +2,7 @@
 
 from pathlib import Path
 
+import numpy as np
 from foamlib import FoamFile
 
 contents = r"""
@@ -49,18 +50,18 @@ def test_get_cells_coords(tmp_path: Path) -> None:
     file = FoamFile(path)
 
     points = file[None]
-    assert isinstance(points, list)
+    assert isinstance(points, np.ndarray)
 
-    assert points[0] == [0, 0, 0]
-    assert points[1] == [0.15707963268, 0, 0]
-    assert points[2] == [0.314159265359, 0, 0]
-    assert points[3] == [0.471238898038, 0, 0]
-    assert points[4] == [0.628318530718, 0, 0]
-    assert points[5] == [0, 0, 0]
-    assert points[6] == [0.15707963268, 0, 0]
-    assert points[7] == [0.314159265359, 0, 0]
-    assert points[8] == [0.471238898038, 0, 0]
-    assert points[9] == [0.628318530718, 0, 0]
+    assert np.array_equal(points[0], [0, 0, 0])
+    assert np.array_equal(points[1], [0.15707963268, 0, 0])
+    assert np.array_equal(points[2], [0.314159265359, 0, 0])
+    assert np.array_equal(points[3], [0.471238898038, 0, 0])
+    assert np.array_equal(points[4], [0.628318530718, 0, 0])
+    assert np.array_equal(points[5], [0, 0, 0])
+    assert np.array_equal(points[6], [0.15707963268, 0, 0])
+    assert np.array_equal(points[7], [0.314159265359, 0, 0])
+    assert np.array_equal(points[8], [0.471238898038, 0, 0])
+    assert np.array_equal(points[9], [0.628318530718, 0, 0])
 
     assert len(points) == 10