Keep the physical dimension of mappings in Domain.export / Domain.from_file

AUTHORS

@@ -15,3 +15,4 @@ Contributors
 * Valentin Carlier
 * Paul Rigor
 * Frederik Schnack
+* Syed Ahnaf Raza

CHANGELOG.md

@@ -18,6 +18,7 @@ All notable changes to this project will be documented in this file.
 
 ### Changed
 
+- #200: Store `ldim` and `pdim` instead of `dim` in the topology of domain files. This is a breaking change: `Domain.from_file()` rejects files with `dim`.
 - #190: Convert the README to Markdown and expand the installation, testing, and documentation-build instructions.
 - #194: Use explicit internal imports and load optional dependencies only when needed.
 - #194: Remove the leading underscore from `partial_derivatives`, `logical_partial_derivatives` (`sympde.topology.derivatives`) and `split_test_function` (`sympde.expr.evaluation`), and add them to `sympde.api`.
@@ -37,6 +38,7 @@ All notable changes to this project will be documented in this file.
 
 ### Fixed
 
+- #200: Keep the physical dimension of the mappings in `Domain.export()` and `Domain.from_file()`.
 - #140: Fix the rendering of the `Domain.join()` API documentation.
 - #166: Restore warning-free documentation builds, including bibliography and mathematical notation support.
 - #106: Preserve symbol assumptions and free-symbol information with modern SymPy.

sympde/topology/domain.py

@@ -31,7 +31,6 @@
 # TODO fix circular dependency between domain and mapping
 if TYPE_CHECKING:
     from sympde.topology.mapping import Mapping
-# TODO add pdim
 
 iterable_types = (tuple, list, Tuple, Union)
 
@@ -336,7 +335,11 @@ def get_interior(self, name):
 
     def todict(self):
         name         = str(self.name)
-        dim          = str(self.dim)
+        ldim         = str(self.dim)
+        pdims        = {m.pdim for m in self.mappings.values() if m is not None}
+        if len(pdims) > 1:
+            raise ValueError(f'all mappings must have the same physical dimension, got {sorted(pdims)}')
+        pdim         = str(pdims.pop()) if pdims else ldim
         interior     = self.interior.todict()
         boundary     = self.boundary.todict()
         connectivity = self.connectivity.todict()
@@ -346,7 +349,8 @@ def todict(self):
             dtype = 'None'
 
         d = {'name':         name,
-             'dim':          dim,
+             'ldim':         ldim,
+             'pdim':         pdim,
              'dtype':        dtype,
              'interior':     interior,
              'boundary':     boundary,
@@ -404,8 +408,13 @@ def from_file(cls, filename):
         h5  = h5py.File( filename, mode='r' )
         yml = yaml.load( h5['topology.yml'][()], Loader=yaml.SafeLoader )
 
+        if 'dim' in yml:
+            raise ValueError(f"'{filename}' was written by an older version of SymPDE:"
+                             " replace its 'dim' key with the keys 'ldim' and 'pdim'")
+
         domain_name    = yml['name']
-        dim            = int(yml['dim'])
+        ldim           = int(yml['ldim'])
+        pdim           = int(yml['pdim'])
         dtype          = yml['dtype']
         d_interior     = yml['interior']
         d_boundary     = yml['boundary']
@@ -422,7 +431,7 @@ def from_file(cls, filename):
 
         constructors = [globals()[dt['type']] for dt in dtype]
         interiors    = [cs(i['name'], **dt['parameters']) for cs,i,dt in zip(constructors, d_interior, dtype)]
-        mappings     = [Mapping(I['mapping'], dim=dim) if I.get('mapping', "None") != "None" else None for I in d_interior]
+        mappings     = [Mapping(I['mapping'], ldim=ldim, pdim=pdim) if I.get('mapping', "None") != "None" else None for I in d_interior]
         domains      = [mapping(i) if mapping else i for i,mapping in zip(interiors, mappings)]
         patch_index  = {I.name:ind for ind,I in enumerate(interiors)}
 

sympde/topology/tests/test_topology.py

@@ -1,5 +1,9 @@
 import os
 
+import h5py
+import numpy as np
+import pytest
+import yaml
 from sympy import preorder_traversal
 from sympy.core import Basic
 
@@ -114,6 +118,40 @@ def test_topology_1(tmp_path):
 
     assert D.todict() == Omega.todict()
 
+#==============================================================================
+def test_export_surface_mapping(tmp_path):
+    domain = Mapping('F', ldim=2, pdim=3)(Square('Omega'))
+
+    filename = str(tmp_path / 'surface.h5')
+    domain.export(filename)
+    D = Domain.from_file(filename)
+
+    assert D.mapping.ldim == 2
+    assert D.mapping.pdim == 3
+    assert D.todict() == domain.todict()
+
+#==============================================================================
+def test_from_file_with_old_dim_key(tmp_path):
+    yml = Mapping('F', dim=2)(Square('Omega')).todict()
+    yml['dim'] = yml.pop('ldim')
+    del yml['pdim']
+
+    filename = str(tmp_path / 'old.h5')
+    with h5py.File(filename, mode='w') as h5:
+        h5['topology.yml'] = np.array(yaml.dump(data=yml, sort_keys=None), dtype='S')
+
+    with pytest.raises(ValueError, match="older version of SymPDE"):
+        Domain.from_file(filename)
+
+#==============================================================================
+def test_export_mappings_with_different_pdim(tmp_path):
+    A = Mapping('FA', ldim=2, pdim=3)(Square('A'))
+    B = Mapping('FB', ldim=2, pdim=2)(Square('B'))
+
+    with pytest.raises(ValueError):
+        domain = Domain.join([A, B], [((0, 0, 1), (1, 0, -1), 1)], name='Omega')
+        domain.export(str(tmp_path / 'omega.h5'))
+
 #==============================================================================
 def test_domain_1():
     Omega_1 = InteriorDomain('Omega_1', dim=2)