diff --git a/psydac/api/discretization.py b/psydac/api/discretization.py index 69dc89679..8c3299356 100644 --- a/psydac/api/discretization.py +++ b/psydac/api/discretization.py @@ -460,9 +460,14 @@ def discretize_space(V, domain_h, *, degree=None, multiplicity=None, knots=None, assert len(ncells) == len(periodic) == len(degree_i) == len(multiplicity_i) == len(min_coords) == len(max_coords) if knots is None: - # Create uniform grid - grids = [np.linspace(xmin, xmax, num=ne + 1) - for xmin, xmax, ne in zip(min_coords, max_coords, ncells)] + # Check if grid is provided in domain_h + if hasattr(domain_h, 'grid') and domain_h.grid is not None and interior.name in domain_h.grid: + # Use provided grid of breakpoints + grids = domain_h.grid[interior.name] + else: + # Create uniform grid + grids = [np.linspace(xmin, xmax, num=ne + 1) + for xmin, xmax, ne in zip(min_coords, max_coords, ncells)] # Create 1D finite element spaces and precompute quadrature data spaces[i] = [SplineSpace( p, multiplicity=m, grid=grid , periodic=P) for p,m,grid,P in zip(degree_i, multiplicity_i,grids, periodic)] @@ -522,23 +527,31 @@ def discretize_space(V, domain_h, *, degree=None, multiplicity=None, knots=None, #============================================================================== -def discretize_domain(domain, *, filename=None, ncells=None, periodic=None, comm=None, mpi_dims_mask=None): +def discretize_domain(domain, *, filename=None, ncells=None, periodic=None, comm=None, mpi_dims_mask=None, grid=None): if comm is not None: # Create a copy of the communicator comm = comm.Dup() - if not (filename or ncells): - raise ValueError("Must provide either 'filename' or 'ncells'") + if not (filename or ncells or grid): + raise ValueError("Must provide either 'filename' or 'ncells' or 'grid'") elif filename and ncells: raise ValueError("Cannot provide both 'filename' and 'ncells'") + + elif filename and grid: + raise ValueError("Cannot provide both 'filename' and 'grid'") elif filename: return Geometry(filename=filename, comm=comm) elif ncells: - return Geometry.from_topological_domain(domain, ncells, periodic=periodic, comm=comm, mpi_dims_mask=mpi_dims_mask) + # Validate grid parameter if provided - basic validation only + if grid is not None: + if not isinstance(grid, (list, tuple, dict)): + raise TypeError("Grid must be a list, tuple, or dict") + + return Geometry.from_topological_domain(domain, ncells, periodic=periodic, comm=comm, mpi_dims_mask=mpi_dims_mask, grid=grid) #============================================================================== def discretize(a, *args, **kwargs): diff --git a/psydac/api/tests/test_discretization_grid.py b/psydac/api/tests/test_discretization_grid.py new file mode 100644 index 000000000..ab2a6b8a7 --- /dev/null +++ b/psydac/api/tests/test_discretization_grid.py @@ -0,0 +1,760 @@ +# -*- coding: utf-8 -*- +""" +Test module for discretization functions with custom grid parameter. + +This module tests the integration between the grid parameter in Geometry +and its usage in the discretize_space function. +""" + +import pytest +import numpy as np + +from sympde.topology import Square, Line, Cube +from sympde.topology import ScalarFunctionSpace + +from psydac.api.discretization import discretize + + +def test_discretize_space_custom_grid_1d(): + """Test discretize_space with custom grid in 1D.""" + + # Create a Line domain + domain = Line('Omega', bounds=(0.0, 1.0)) + + # Define custom grid with non-uniform spacing + custom_grid = [np.array([0.0, 0.1, 0.4, 0.7, 0.9, 1.0])] # 5 cells + + # Create geometry with custom grid + ncells = [5] + domain_h = discretize(domain, ncells=ncells, grid=custom_grid) + + # Create symbolic space + symbolic_space = ScalarFunctionSpace('V', domain) + + # Discretize the space - this should use our custom grid + degree = [2] + discrete_space = discretize(symbolic_space, domain_h=domain_h, degree=degree) + + # Verify that the breaks match our custom grid + breaks = discrete_space.spaces[0].breaks + assert np.allclose(breaks, custom_grid[0]), f"Expected {custom_grid[0]}, got {breaks}" + + # Verify that the number of cells is correct + assert discrete_space.spaces[0].ncells == 5 + + +def test_discretize_space_custom_grid_2d(): + """Test discretize_space with custom grid in 2D.""" + + # Create a Square domain + domain = Square('Omega', bounds1=(0.0, 1.0), bounds2=(0.0, 1.0)) + + # Define custom grid with non-uniform spacing + custom_grid = [ + np.array([0.0, 0.2, 0.6, 1.0]), # 3 cells in x-direction + np.array([0.0, 0.3, 0.5, 0.8, 1.0]) # 4 cells in y-direction + ] + + # Create geometry with custom grid + ncells = [3, 4] + domain_h = discretize(domain, ncells=ncells, grid=custom_grid) + + + # Create symbolic space + symbolic_space = ScalarFunctionSpace('V', domain) + + # Discretize the space - this should use our custom grid + degree = [2, 2] + discrete_space = discretize(symbolic_space, domain_h=domain_h, degree=degree) + + # Verify that the breaks match our custom grid for both directions + for i, (space, expected_grid) in enumerate(zip(discrete_space.spaces, custom_grid)): + breaks = space.breaks + assert np.allclose(breaks, expected_grid), f"Direction {i}: Expected {expected_grid}, got {breaks}" + + # Verify number of cells + expected_ncells = len(expected_grid) - 1 + assert space.ncells == expected_ncells + + +def test_discretize_space_fallback_uniform_grid(): + """Test that when grid=None, uniform grid is used.""" + + # Create a Square domain + domain = Square('Omega', bounds1=(0.0, 1.0), bounds2=(0.0, 1.0)) + + # Create geometry without custom grid (grid=None) + ncells = [4, 4] + domain_h = discretize(domain, ncells=ncells) + + # Create symbolic space + symbolic_space = ScalarFunctionSpace('V', domain) + + # Discretize the space - this should use uniform grid + degree = [2, 2] + discrete_space = discretize(symbolic_space, domain_h=domain_h, degree=degree) + + # Verify that uniform grids are used + for i, space in enumerate(discrete_space.spaces): + breaks = space.breaks + + # For uniform grid, breaks should be evenly spaced + expected_breaks = np.linspace(0.0, 1.0, ncells[i] + 1) + assert np.allclose(breaks, expected_breaks), f"Expected uniform grid {expected_breaks}, got {breaks}" + + +@pytest.mark.parametrize("degree", [[1], [2], [3]]) +def test_discretize_space_custom_grid_1d_various_degrees(degree): + """Test discretize_space with custom grid for various degrees in 1D.""" + + # Create a Line domain + domain = Line('Omega', bounds=(0.0, 1.0)) + + # Define custom grid + custom_grid = [np.array([0.0, 0.2, 0.5, 0.8, 1.0])] # 4 cells + + # Create geometry with custom grid + ncells = [4] + + domain_h = discretize(domain, ncells=ncells, grid=custom_grid) + + # Create symbolic space + symbolic_space = ScalarFunctionSpace('V', domain) + + # Discretize the space + discrete_space = discretize(symbolic_space, domain_h=domain_h, degree=degree) + + # Verify that the breaks match our custom grid + breaks = discrete_space.spaces[0].breaks + assert np.allclose(breaks, custom_grid[0]) + + # Verify degree + assert discrete_space.spaces[0].degree == degree[0] + + +@pytest.mark.parametrize("ncells,expected_cells", [ + ([3], 3), + ([5], 5), + ([8], 8) +]) +def test_discretize_space_custom_grid_1d_various_ncells(ncells, expected_cells): + """Test discretize_space with custom grid for various numbers of cells in 1D.""" + + # Create a Line domain + domain = Line('Omega', bounds=(0.0, 1.0)) + + # Define custom grid with specified number of cells + custom_grid = [np.linspace(0.0, 1.0, expected_cells + 1)] + + # Create geometry with custom grid + domain_h = discretize(domain, ncells=ncells, grid=custom_grid) + + # Create symbolic space + symbolic_space = ScalarFunctionSpace('V', domain) + + # Discretize the space + degree = [2] + discrete_space = discretize(symbolic_space, domain_h=domain_h, degree=degree) + + # Verify that the breaks match our custom grid + breaks = discrete_space.spaces[0].breaks + assert np.allclose(breaks, custom_grid[0]) + + # Verify number of cells + assert discrete_space.spaces[0].ncells == expected_cells + + +def test_discretize_space_custom_grid_2d_non_uniform(): + """Test discretize_space with highly non-uniform custom grid in 2D.""" + + # Create a Square domain + domain = Square('Omega', bounds1=(0.0, 1.0), bounds2=(0.0, 1.0)) + + # Define highly non-uniform custom grid + custom_grid = [ + np.array([0.0, 0.01, 0.02, 0.5, 0.98, 0.99, 1.0]), # 6 cells, clustered at boundaries + np.array([0.0, 0.1, 0.9, 1.0]) # 3 cells, clustered at boundaries + ] + + # Create geometry with custom grid + ncells = [6, 3] + + domain_h = discretize(domain, ncells=ncells, grid=custom_grid) + + # Create symbolic space + symbolic_space = ScalarFunctionSpace('V', domain) + + # Discretize the space + degree = [3, 2] + discrete_space = discretize(symbolic_space, domain_h=domain_h, degree=degree) + + # Verify that the breaks match our custom grid for both directions + for i, (space, expected_grid) in enumerate(zip(discrete_space.spaces, custom_grid)): + breaks = space.breaks + assert np.allclose(breaks, expected_grid), f"Direction {i}: Expected {expected_grid}, got {breaks}" + + # Verify that spacing is indeed non-uniform + spacings = np.diff(breaks) + assert not np.allclose(spacings, spacings[0]), f"Grid should be non-uniform in direction {i}" + + +def test_discretize_space_custom_grid_3d(): + """Test discretize_space with custom grid in 3D.""" + + # Create a Cube domain + domain = Cube('Omega', bounds1=(0.0, 1.0), bounds2=(0.0, 1.0), bounds3=(0.0, 1.0)) + + # Define custom grid with non-uniform spacing in 3D + custom_grid = [ + np.array([0.0, 0.3, 0.7, 1.0]), # 3 cells in x-direction + np.array([0.0, 0.2, 0.5, 0.8, 1.0]), # 4 cells in y-direction + np.array([0.0, 0.4, 1.0]) # 2 cells in z-direction + ] + + # Create geometry with custom grid + ncells = [3, 4, 2] + + domain_h = discretize(domain, ncells=ncells, grid=custom_grid) + + # Create symbolic space + symbolic_space = ScalarFunctionSpace('V', domain) + + # Discretize the space - this should use our custom grid + degree = [2, 2, 2] + discrete_space = discretize(symbolic_space, domain_h=domain_h, degree=degree) + + # Verify that the breaks match our custom grid for all directions + for i, (space, expected_grid) in enumerate(zip(discrete_space.spaces, custom_grid)): + breaks = space.breaks + assert np.allclose(breaks, expected_grid), f"Direction {i}: Expected {expected_grid}, got {breaks}" + + # Verify number of cells + expected_ncells = len(expected_grid) - 1 + assert space.ncells == expected_ncells + + +def test_discretize_space_custom_grid_3d_non_uniform(): + """Test discretize_space with highly non-uniform custom grid in 3D.""" + + # Create a Cube domain + domain = Cube('Omega', bounds1=(0.0, 1.0), bounds2=(0.0, 1.0), bounds3=(0.0, 1.0)) + + # Define highly non-uniform custom grid - clustered at boundaries + custom_grid = [ + np.array([0.0, 0.01, 0.5, 0.99, 1.0]), # 4 cells in x, clustered at boundaries + np.array([0.0, 0.1, 0.9, 1.0]), # 3 cells in y, clustered at boundaries + np.array([0.0, 0.02, 0.03, 0.97, 0.98, 1.0]) # 5 cells in z, very clustered + ] + + # Create geometry with custom grid + ncells = [4, 3, 5] + + domain_h = discretize(domain, ncells=ncells, grid=custom_grid) + + # Create symbolic space + symbolic_space = ScalarFunctionSpace('V', domain) + + # Discretize the space + degree = [2, 3, 2] + discrete_space = discretize(symbolic_space, domain_h=domain_h, degree=degree) + + # Verify that the breaks match our custom grid for all directions + for i, (space, expected_grid) in enumerate(zip(discrete_space.spaces, custom_grid)): + breaks = space.breaks + assert np.allclose(breaks, expected_grid), f"Direction {i}: Expected {expected_grid}, got {breaks}" + + # Verify that spacing is indeed non-uniform + spacings = np.diff(breaks) + assert not np.allclose(spacings, spacings[0]), f"Grid should be non-uniform in direction {i}" + + +@pytest.mark.parametrize("degree", [[1, 1, 1], [2, 2, 2], [3, 2, 1], [1, 3, 2]]) +def test_discretize_space_custom_grid_3d_various_degrees(degree): + """Test discretize_space with custom grid for various degrees in 3D.""" + + # Create a Cube domain + domain = Cube('Omega', bounds1=(0.0, 1.0), bounds2=(0.0, 1.0), bounds3=(0.0, 1.0)) + + # Define custom grid + custom_grid = [ + np.array([0.0, 0.25, 0.5, 0.75, 1.0]), # 4 cells in x + np.array([0.0, 0.33, 0.67, 1.0]), # 3 cells in y + np.array([0.0, 0.5, 1.0]) # 2 cells in z + ] + + # Create geometry with custom grid + ncells = [4, 3, 2] + + domain_h = discretize(domain, ncells=ncells, grid=custom_grid) + + # Create symbolic space + symbolic_space = ScalarFunctionSpace('V', domain) + + # Discretize the space + discrete_space = discretize(symbolic_space, domain_h=domain_h, degree=degree) + + # Verify that the breaks match our custom grid for all directions + for i, (space, expected_grid) in enumerate(zip(discrete_space.spaces, custom_grid)): + breaks = space.breaks + assert np.allclose(breaks, expected_grid) + + # Verify degree + assert space.degree == degree[i] + + +def test_discretize_space_custom_grid_3d_preserves_boundary_values(): + """Test that custom grid preserves boundary values in 3D.""" + + # Create a Cube domain with specific bounds + domain = Cube('Omega', bounds1=(-1.0, 2.0), bounds2=(0.5, 3.5), bounds3=(-0.5, 1.5)) + + # Define custom grid that respects the domain bounds + custom_grid = [ + np.array([-1.0, 0.0, 1.0, 2.0]), # x: 3 cells from -1 to 2 + np.array([0.5, 1.5, 2.5, 3.5]), # y: 3 cells from 0.5 to 3.5 + np.array([-0.5, 0.0, 0.5, 1.0, 1.5]) # z: 4 cells from -0.5 to 1.5 + ] + + # Create geometry with custom grid + ncells = [3, 3, 4] + + domain_h = discretize(domain, ncells=ncells, grid=custom_grid) + + # Create symbolic space + symbolic_space = ScalarFunctionSpace('V', domain) + + # Discretize the space + degree = [2, 2, 2] + discrete_space = discretize(symbolic_space, domain_h=domain_h, degree=degree) + + # Verify that the breaks match our custom grid for all directions + for i, (space, expected_grid) in enumerate(zip(discrete_space.spaces, custom_grid)): + breaks = space.breaks + assert np.allclose(breaks, expected_grid) + + # Verify boundary values are preserved + assert breaks[0] == expected_grid[0] # First boundary + assert breaks[-1] == expected_grid[-1] # Last boundary + + # Verify specific boundary values + assert discrete_space.spaces[0].breaks[0] == -1.0 # x min + assert discrete_space.spaces[0].breaks[-1] == 2.0 # x max + assert discrete_space.spaces[1].breaks[0] == 0.5 # y min + assert discrete_space.spaces[1].breaks[-1] == 3.5 # y max + assert discrete_space.spaces[2].breaks[0] == -0.5 # z min + assert discrete_space.spaces[2].breaks[-1] == 1.5 # z max + + +def test_discretize_space_fallback_uniform_grid_3d(): + """Test that when grid=None, uniform grid is used in 3D.""" + + # Create a Cube domain + domain = Cube('Omega', bounds1=(0.0, 1.0), bounds2=(0.0, 1.0), bounds3=(0.0, 1.0)) + + # Create geometry without custom grid (grid=None) + ncells = [3, 4, 2] + + domain_h = discretize(domain, ncells=ncells, grid=None) + + # Create symbolic space + symbolic_space = ScalarFunctionSpace('V', domain) + + # Discretize the space - this should use uniform grid + degree = [2, 2, 2] + discrete_space = discretize(symbolic_space, domain_h=domain_h, degree=degree) + + # Verify that uniform grids are used in all directions + for i, space in enumerate(discrete_space.spaces): + breaks = space.breaks + + # For uniform grid, breaks should be evenly spaced + expected_breaks = np.linspace(0.0, 1.0, ncells[i] + 1) + assert np.allclose(breaks, expected_breaks), f"Expected uniform grid {expected_breaks}, got {breaks}" + + +# ============================================================================== +# GRID SECURITY TESTS +# ============================================================================== + +def test_grid_security_type_validation(): + """Test that grid type validation works correctly.""" + + domain = Line('Omega', bounds=(0.0, 1.0)) + symbolic_space = ScalarFunctionSpace('V', domain) + + # Test invalid grid type + ncells = [4] + + + with pytest.raises(TypeError, match="Grid must be a list, tuple, or dict"): + domain_h = discretize(domain, ncells=ncells, grid="invalid") + discretize(symbolic_space, domain_h=domain_h, degree=[2]) + + +def test_grid_security_dimension_mismatch(): + """Test that grid dimension validation works correctly.""" + + domain = Line('Omega', bounds=(0.0, 1.0)) + symbolic_space = ScalarFunctionSpace('V', domain) + + # Test 2D grid for 1D domain + ncells = [4] + + + with pytest.raises(ValueError, match="Grid dimensions .* must match domain dimensions"): + invalid_grid = [[0, 0.25, 0.5, 0.75, 1], [0, 0.5, 1]] # 2D grid for 1D domain + domain_h = discretize(domain, ncells=ncells, grid=invalid_grid) + discretize(symbolic_space, domain_h=domain_h, degree=[2]) + + +def test_grid_security_length_consistency(): + """Test that grid length vs ncells consistency is enforced.""" + + domain = Line('Omega', bounds=(0.0, 1.0)) + symbolic_space = ScalarFunctionSpace('V', domain) + + # Test grid length inconsistent with ncells + ncells = [4] + + + with pytest.raises(ValueError, match="grid length .* must be ncells\\+1"): + invalid_grid = [[0, 0.33, 0.66, 1]] # 3 cells but ncells=[4] + domain_h = discretize(domain, ncells=ncells, grid=invalid_grid) + discretize(symbolic_space, domain_h=domain_h, degree=[2]) + + +def test_grid_security_monotonic_order(): + """Test that grid points must be in strictly increasing order.""" + + domain = Line('Omega', bounds=(0.0, 1.0)) + symbolic_space = ScalarFunctionSpace('V', domain) + + # Test non-monotonic grid + ncells = [4] + + + with pytest.raises(ValueError, match="grid points must be strictly increasing"): + invalid_grid = [[0, 0.5, 0.25, 0.75, 1]] # Not monotonic + domain_h = discretize(domain, ncells=ncells, grid=invalid_grid) + discretize(symbolic_space, domain_h=domain_h, degree=[2]) + + +def test_grid_security_boundary_mismatch(): + """Test that grid boundaries must match domain boundaries.""" + + domain = Line('Omega', bounds=(0.0, 1.0)) + symbolic_space = ScalarFunctionSpace('V', domain) + + # Test start boundary mismatch + ncells = [4] + + + with pytest.raises(ValueError, match="grid start .* must match domain minimum"): + invalid_grid = [[0.1, 0.3, 0.6, 0.8, 1]] # Start doesn't match domain [0,1] + domain_h = discretize(domain, ncells=ncells, grid=invalid_grid) + discretize(symbolic_space, domain_h=domain_h, degree=[2]) + + # Test end boundary mismatch + with pytest.raises(ValueError, match="grid end .* must match domain maximum"): + invalid_grid = [[0, 0.2, 0.4, 0.6, 0.9]] # End doesn't match domain [0,1] + domain_h = discretize(domain, ncells=ncells, grid=invalid_grid) + discretize(symbolic_space, domain_h=domain_h, degree=[2]) + + +def test_grid_security_2d_validation(): + """Test grid security checks work correctly in 2D.""" + + domain = Square('Omega', bounds1=(0.0, 1.0), bounds2=(0.0, 1.0)) + symbolic_space = ScalarFunctionSpace('V', domain) + + # Test inconsistent 2D grid dimensions + ncells = [3, 2] + + + with pytest.raises(ValueError, match="Dimension 1: grid length .* must be ncells\\+1"): + invalid_grid = [ + [0, 0.33, 0.66, 1], # 3 cells ✓ + [0, 0.5, 1, 1.5] # 4 points but should be 3 ✗ + ] + domain_h = discretize(domain, ncells=ncells, grid=invalid_grid) + discretize(symbolic_space, domain_h=domain_h, degree=[2, 2]) + + +def test_grid_security_boundary_tolerance(): + """Test that small numerical differences in boundaries are tolerated.""" + + domain = Line('Omega', bounds=(0.0, 1.0)) + symbolic_space = ScalarFunctionSpace('V', domain) + + # Test grid within tolerance (should work) + ncells = [2] + + + # Small boundary differences (within tolerance) + tolerance_grid = [[1e-15, 0.5, 1.0 - 1e-15]] + domain_h = discretize(domain, ncells=ncells, grid=tolerance_grid) + + # This should work without raising an exception + discrete_space = discretize(symbolic_space, domain_h=domain_h, degree=[2]) + assert discrete_space.spaces[0].ncells == 2 + + +def test_grid_security_boundary_tolerance_exceeded(): + """Test that boundary differences outside tolerance are rejected.""" + + domain = Line('Omega', bounds=(0.0, 1.0)) + symbolic_space = ScalarFunctionSpace('V', domain) + + # Test boundaries outside tolerance + ncells = [2] + + + with pytest.raises(ValueError, match="grid start .* must match domain minimum"): + # Boundaries outside tolerance + invalid_grid = [[1e-10, 0.5, 1.0 - 1e-10]] + domain_h = discretize(domain, ncells=ncells, grid=invalid_grid) + discretize(symbolic_space, domain_h=domain_h, degree=[2]) + + +def test_grid_security_3d_validation(): + """Test grid security checks work correctly in 3D.""" + + domain = Cube('Omega', bounds1=(0.0, 1.0), bounds2=(0.0, 1.0), bounds3=(0.0, 1.0)) + symbolic_space = ScalarFunctionSpace('V', domain) + + # Test valid 3D grid (should work) + ncells = [2, 3, 2] + + + valid_grid = [ + [0, 0.5, 1.0], # 2 cells in x + [0, 0.33, 0.67, 1.0], # 3 cells in y + [0, 0.4, 1.0] # 2 cells in z + ] + domain_h = discretize(domain, ncells=ncells, grid=valid_grid) + discrete_space = discretize(symbolic_space, domain_h=domain_h, degree=[2, 2, 2]) + + # Verify all dimensions work correctly + for i, expected_ncells in enumerate([2, 3, 2]): + assert discrete_space.spaces[i].ncells == expected_ncells + + # Test invalid 3D grid (dimension 2 has wrong length) + with pytest.raises(ValueError, match="Dimension 2: grid length .* must be ncells\\+1"): + invalid_grid = [ + [0, 0.5, 1.0], # 2 cells ✓ + [0, 0.33, 0.67, 1.0], # 3 cells ✓ + [0, 0.25, 0.5, 0.75, 1.0] # 4 cells instead of 2 ✗ + ] + domain_h = discretize(domain, ncells=ncells, grid=invalid_grid) + discretize(symbolic_space, domain_h=domain_h, degree=[2, 2, 2]) + + +@pytest.mark.parametrize("error_type,grid_config,error_pattern", [ + ("type", "invalid_string", "Grid must be a list, tuple, or dict"), + ("length", [[0, 0.5, 1]], "grid length .* must be ncells\\+1"), + ("monotonic", [[0, 0.8, 0.5, 0.9, 1]], "grid points must be strictly increasing"), + ("boundary_start", [[0.1, 0.3, 0.6, 0.8, 1]], "grid start .* must match domain minimum"), + ("boundary_end", [[0, 0.2, 0.4, 0.6, 0.9]], "grid end .* must match domain maximum"), +]) +def test_grid_security_parametrized_errors(error_type, grid_config, error_pattern): + """Parametrized test for various grid security error conditions.""" + + domain = Line('Omega', bounds=(0.0, 1.0)) + symbolic_space = ScalarFunctionSpace('V', domain) + ncells = [4] + + + if error_type == "type": + with pytest.raises(TypeError, match=error_pattern): + domain_h = discretize(domain, ncells=ncells, grid=grid_config) + discretize(symbolic_space, domain_h=domain_h, degree=[2]) + else: + with pytest.raises(ValueError, match=error_pattern): + geo = discretize(domain, ncells=ncells, grid=grid_config) + discretize(symbolic_space, domain_h=geo, degree=[2]) + + +def test_grid_security_preserves_functionality(): + """Test that security checks don't interfere with normal functionality.""" + + # Test that valid grids still work perfectly + domain = Square('Omega', bounds1=(0.0, 1.0), bounds2=(0.0, 1.0)) + symbolic_space = ScalarFunctionSpace('V', domain) + + # Define valid but complex custom grid + ncells = [4, 3] + + + valid_grid = [ + [0, 0.1, 0.3, 0.7, 1.0], # 4 cells, non-uniform + [0, 0.2, 0.8, 1.0] # 3 cells, non-uniform + ] + + domain_h = discretize(domain, ncells=ncells, grid=valid_grid) + discrete_space = discretize(symbolic_space, domain_h=domain_h, degree=[3, 2]) + + # Verify that everything works as expected + for i, (space, expected_grid) in enumerate(zip(discrete_space.spaces, valid_grid)): + assert np.allclose(space.breaks, expected_grid) + assert space.ncells == len(expected_grid) - 1 + assert space.degree == [3, 2][i] + + # Verify that the grid is indeed non-uniform + for i, space in enumerate(discrete_space.spaces): + spacings = np.diff(space.breaks) + assert not np.allclose(spacings, spacings[0]), f"Grid should be non-uniform in direction {i}" + +def test_grid_with_none_values_2d(): + """Test grid with None values in 2D - mixed None and custom grids.""" + + domain = Square('Omega', bounds1=(0.0, 1.0), bounds2=(0.0, 1.0)) + ncells = [3, 4] + + # None for dimension 0, custom for dimension 1 + grid_mixed = [ + None, # Uniform grid (will be auto-generated) + [0.0, 0.2, 0.6, 0.8, 1.0] # Custom grid + ] + + domain_h = discretize(domain, ncells=ncells, grid=grid_mixed) + + # Verify the grid was processed correctly + actual_grid = domain_h.grid + assert 'Omega' in actual_grid + + # Check dimension 0 (should be uniform) + expected_dim0 = np.linspace(0.0, 1.0, 4) # 3 cells + 1 = 4 points + assert np.allclose(actual_grid['Omega'][0], expected_dim0) + + # Check dimension 1 (should be preserved custom) + expected_dim1 = [0.0, 0.2, 0.6, 0.8, 1.0] + assert np.allclose(actual_grid['Omega'][1], expected_dim1) + +def test_grid_with_all_none_values_2d(): + """Test grid with all None values in 2D - should generate uniform grids.""" + + domain = Square('Omega', bounds1=(0.0, 1.0), bounds2=(0.0, 1.0)) + ncells = [2, 3] + + # All None values + grid_all_none = [None, None] + + domain_h = discretize(domain, ncells=ncells, grid=grid_all_none) + + # Verify uniform grids were generated + actual_grid = domain_h.grid + assert 'Omega' in actual_grid + + # Both dimensions should be uniform + expected_dim0 = np.linspace(0.0, 1.0, 3) # 2 cells + 1 = 3 points + expected_dim1 = np.linspace(0.0, 1.0, 4) # 3 cells + 1 = 4 points + + assert np.allclose(actual_grid['Omega'][0], expected_dim0) + assert np.allclose(actual_grid['Omega'][1], expected_dim1) + +def test_grid_with_none_values_3d(): + """Test grid with None values in 3D - mixed None and custom grids.""" + + domain = Cube('Omega', bounds1=(0.0, 1.0), bounds2=(0.0, 1.0), bounds3=(0.0, 1.0)) + ncells = [2, 3, 2] + + # Mixed None and custom grids + grid_mixed_3d = [ + None, # Uniform for dim 0 + [0.0, 0.33, 0.67, 1.0], # Custom for dim 1 + None # Uniform for dim 2 + ] + + domain_h = discretize(domain, ncells=ncells, grid=grid_mixed_3d) + + # Verify the grid was processed correctly + actual_grid = domain_h.grid + assert 'Omega' in actual_grid + + # Check dimension 0 (should be uniform) + expected_dim0 = np.linspace(0.0, 1.0, 3) # 2 cells + 1 = 3 points + assert np.allclose(actual_grid['Omega'][0], expected_dim0) + + # Check dimension 1 (should be preserved custom) + expected_dim1 = [0.0, 0.33, 0.67, 1.0] + assert np.allclose(actual_grid['Omega'][1], expected_dim1) + + # Check dimension 2 (should be uniform) + expected_dim2 = np.linspace(0.0, 1.0, 3) # 2 cells + 1 = 3 points + assert np.allclose(actual_grid['Omega'][2], expected_dim2) + +def test_grid_with_none_values_1d(): + """Test grid with None values in 1D.""" + + domain = Line('Omega', bounds=(0.0, 1.0)) + ncells = [4] + + # None for 1D + grid_none_1d = [None] + + domain_h = discretize(domain, ncells=ncells, grid=grid_none_1d) + + # Verify uniform grid was generated + actual_grid = domain_h.grid + assert 'Omega' in actual_grid + + expected_grid = np.linspace(0.0, 1.0, 5) # 4 cells + 1 = 5 points + assert np.allclose(actual_grid['Omega'], expected_grid) + +def test_grid_none_values_respects_domain_bounds(): + """Test that None values generate uniform grids respecting domain bounds.""" + + # Test with non-standard domain bounds + domain = Square('Omega', bounds1=(2.0, 5.0), bounds2=(-1.0, 3.0)) + ncells = [3, 2] + + grid_all_none = [None, None] + + domain_h = discretize(domain, ncells=ncells, grid=grid_all_none) + + actual_grid = domain_h.grid + assert 'Omega' in actual_grid + + # Check that bounds are respected + expected_dim0 = np.linspace(2.0, 5.0, 4) # bounds1 with 3 cells + 1 = 4 points + expected_dim1 = np.linspace(-1.0, 3.0, 3) # bounds2 with 2 cells + 1 = 3 points + + assert np.allclose(actual_grid['Omega'][0], expected_dim0) + assert np.allclose(actual_grid['Omega'][1], expected_dim1) + +def test_grid_none_mixed_with_discretize_space(): + """Test that grid with None values works with discretize_space.""" + + domain = Square('Omega', bounds1=(0.0, 1.0), bounds2=(0.0, 1.0)) + symbolic_space = ScalarFunctionSpace('V', domain) + + ncells = [2, 3] + grid_mixed = [ + [0.0, 0.3, 1.0], # Custom for dim 0 (2 cells) + None # Uniform for dim 1 + ] + + domain_h = discretize(domain, ncells=ncells, grid=grid_mixed) + discrete_space = discretize(symbolic_space, domain_h=domain_h, degree=[2, 2]) + + # Verify it works + assert discrete_space.spaces[0].ncells == 2 + assert discrete_space.spaces[1].ncells == 3 + + # Verify the grid + actual_grid = domain_h.grid + expected_dim0 = [0.0, 0.3, 1.0] + expected_dim1 = np.linspace(0.0, 1.0, 4) # 3 cells + 1 = 4 points + + assert np.allclose(actual_grid['Omega'][0], expected_dim0) + assert np.allclose(actual_grid['Omega'][1], expected_dim1) + + +# CLEAN UP SYMPY NAMESPACE +def teardown_module(): + from sympy.core import cache + cache.clear_cache() + +def teardown_function(): + from sympy.core import cache + cache.clear_cache() diff --git a/psydac/cad/geometry.py b/psydac/cad/geometry.py index 066dc8231..6f6f7e20f 100644 --- a/psydac/cad/geometry.py +++ b/psydac/cad/geometry.py @@ -61,6 +61,9 @@ class Geometry( object ): True if the dimension is to be used in the domain decomposition (=default for each dimension). If mpi_dims_mask[i]=False, the i-th dimension will not be decomposed. + grid: list | tuple | dict + The grid of breakpoints in each direction. + If not given, the grid will be constructed from the ncells and periodicity. """ _ldim = None _pdim = None @@ -71,7 +74,7 @@ class Geometry( object ): # Option [1]: from a (domain, mappings) or a file #-------------------------------------------------------------------------- def __init__(self, domain=None, ncells=None, periodic=None, mappings=None, - filename=None, comm=None, mpi_dims_mask=None): + filename=None, comm=None, mpi_dims_mask=None, grid=None): # ... read the geometry if the filename is given if filename is not None: @@ -81,6 +84,7 @@ def __init__(self, domain=None, ncells=None, periodic=None, mappings=None, assert isinstance(domain, Domain) assert isinstance(ncells, dict) assert isinstance(mappings, dict) + assert grid is None or isinstance(grid, dict) if periodic is not None: assert isinstance(periodic, dict) @@ -102,6 +106,14 @@ def __init__(self, domain=None, ncells=None, periodic=None, mappings=None, self._mappings = mappings self._cart = None self._is_parallel = comm is not None + # Process and validate grid parameter if provided + if grid is not None: + # Process None values first, then validate + processed_grid = self._process_grid_with_none_values(domain, ncells, grid) + self._validate_grid_consistency(domain, ncells, processed_grid) + self._grid = processed_grid + else: + self._grid = grid if len(domain) == 1: #name = domain.name @@ -110,6 +122,8 @@ def __init__(self, domain=None, ncells=None, periodic=None, mappings=None, else: ncells = [ncells[itr] for itr in interior_names] periodic = [periodic[itr] for itr in interior_names] + if grid is not None: + grid = [grid[itr] for itr in interior_names] self._ddm = MultiPatchDomainDecomposition(ncells, periodic, comm=comm) else: @@ -150,14 +164,47 @@ def from_discrete_mapping(cls, mapping, comm=None, name=None): ncells = {domain.name: mapping.space.domain_decomposition.ncells} periodic = {domain.name: mapping.space.domain_decomposition.periods} - return Geometry(domain=domain, ncells=ncells, periodic=periodic, mappings=mappings, comm=comm) + return Geometry(domain=domain, ncells=ncells, periodic=periodic, mappings=mappings, comm=comm, grid=None) #-------------------------------------------------------------------------- # Option [3]: discrete topological line/square/cube #-------------------------------------------------------------------------- @classmethod - def from_topological_domain(cls, domain, ncells, *, periodic=None, comm=None, mpi_dims_mask=None): + def from_topological_domain(cls, domain, ncells, *, periodic=None, comm=None, mpi_dims_mask=None, grid=None): + """ + Create a Geometry object from a topological domain and discretization parameters. + + Parameters + ---------- + cls : type + The class type (typically `Geometry`) to instantiate. + + domain : Sympde.topology.Domain + The symbolic domain to be discretized. + + ncells : list | tuple | dict + The number of cells of the discretized topological domain in each direction. + + periodic : list | tuple | dict + The periodicity of the topological domain in each direction. + + comm: MPI.Comm + MPI intra-communicator. + + mpi_dims_mask: list of bool + True if the dimension is to be used in the domain decomposition (=default for each dimension). + If mpi_dims_mask[i]=False, the i-th dimension will not be decomposed. + + grid: list | tuple | dict + The grid of breakpoints in each direction. + If not given, the grid will be constructed from the ncells and periodicity. + + Returns + ------- + geo : Geometry + The constructed Geometry object. + """ interior = domain.interior if not isinstance(interior, Union): interior = [interior] @@ -179,9 +226,222 @@ def from_topological_domain(cls, domain, ncells, *, periodic=None, comm=None, mp if isinstance(periodic, (list, tuple)): periodic = {itr.name:periodic for itr in interior} - geo = Geometry(domain=domain, mappings=mappings, ncells=ncells, periodic=periodic, comm=comm, mpi_dims_mask=mpi_dims_mask) + if isinstance(grid, (list, tuple)): + grid = {itr.name:grid for itr in interior} + + geo = Geometry(domain=domain, mappings=mappings, ncells=ncells, periodic=periodic, comm=comm, mpi_dims_mask=mpi_dims_mask, grid=grid) return geo + + #-------------------------------------------------------------------------- + def _generate_uniform_grid(self, domain_interior, ncells_dim, dim): + """ + Generate a uniform grid for a given dimension. + + Parameters + ---------- + domain_interior : NCubeInterior + The domain interior. + ncells_dim : int + Number of cells in this dimension. + dim : int + The dimension index. + + Returns + ------- + list + Uniform grid points for the dimension. + """ + if domain_interior is not None: + domain_min = domain_interior.min_coords[dim] + domain_max = domain_interior.max_coords[dim] + else: + # Default to [0, 1] if no domain interior available + domain_min = 0.0 + domain_max = 1.0 + + return list(np.linspace(domain_min, domain_max, ncells_dim + 1)) + + #-------------------------------------------------------------------------- + def _process_grid_with_none_values(self, domain, ncells, grid): + """ + Process grid parameter, replacing None values with uniform grids. + + Parameters + ---------- + domain : Sympde.topology.Domain + The symbolic domain. + ncells : dict + Number of cells per patch and direction. + grid : dict + Grid of breakpoints per patch and direction, may contain None values. + + Returns + ------- + dict + Processed grid with None values replaced by uniform grids. + """ + processed_grid = {} + + for patch_name, patch_grid in grid.items(): + if patch_name not in ncells: + raise ValueError(f"Grid patch '{patch_name}' not found in ncells") + + patch_ncells = ncells[patch_name] + + # Get domain interior for this patch + if len(domain) == 1: + domain_interior = domain.interior + else: + # Find the interior corresponding to this patch + domain_interior = None + for interior in domain.interior.args: + if interior.name == patch_name: + domain_interior = interior + break + + if not isinstance(patch_grid, (list, tuple, np.ndarray)): + raise TypeError(f"Grid must be a list or tuple") + + # Process the grid, replacing None values with uniform grids + processed_patch_grid = [] + + # Check if this is a 1D grid or multi-dimensional + # For 1D: either [None] or [points...] where points are not lists + # For nD: [[...], [...], ...] where each inner element is a list/array or None + if len(patch_grid) == 1 and not isinstance(patch_grid[0], (list, tuple, np.ndarray)): + # 1D case with single element that's not a list (could be None or a number) + if patch_grid[0] is None: + # Generate uniform grid for 1D + processed_patch_grid = self._generate_uniform_grid(domain_interior, patch_ncells[0], 0) + else: + # This shouldn't happen in normal 1D usage, but handle gracefully + processed_patch_grid = patch_grid + elif len(patch_grid) > 1 and not any(isinstance(x, (list, tuple, np.ndarray, type(None))) for x in patch_grid): + # 1D case with multiple points [x0, x1, x2, ...] + processed_patch_grid = patch_grid + else: + # Multi-dimensional case: [[...], [...], ...] or [None, [...], None, ...] + for dim in range(len(patch_grid)): + if patch_grid[dim] is None: + # Generate uniform grid for this dimension + uniform_grid = self._generate_uniform_grid(domain_interior, patch_ncells[dim], dim) + processed_patch_grid.append(uniform_grid) + else: + processed_patch_grid.append(patch_grid[dim]) + + processed_grid[patch_name] = processed_patch_grid + + return processed_grid + + #-------------------------------------------------------------------------- + def _validate_grid_consistency(self, domain, ncells, grid): + """ + Validate grid parameter consistency with domain and ncells. + Note: This method expects that None values have already been processed. + + Parameters + ---------- + domain : Sympde.topology.Domain + The symbolic domain. + ncells : dict + Number of cells per patch and direction. + grid : dict + Grid of breakpoints per patch and direction (already processed). + """ + if not isinstance(grid, dict): + raise TypeError("Grid must be a dictionary when domain has multiple patches") + + # Get domain dimension for validation + domain_dim = domain.dim + + # Validate each patch + for patch_name, patch_grid in grid.items(): + if patch_name not in ncells: + raise ValueError(f"Grid patch '{patch_name}' not found in ncells") + patch_ncells = ncells[patch_name] + + # Check if grid is list/tuple/array format + if not isinstance(patch_grid, (list, tuple, np.ndarray)): + raise TypeError(f"Grid must be a list, tuple, or numpy array") + + # Check grid dimensions consistency + if isinstance(patch_grid[0], (list, tuple, np.ndarray)): + # Multi-dimensional grid + grid_dim = len(patch_grid) + else: + # 1D grid + grid_dim = 1 + + if grid_dim != domain_dim: + raise ValueError(f"Grid dimensions ({grid_dim}) must match domain dimensions ({domain_dim})") + + # Get domain interior for boundary validation + if len(domain) == 1: + domain_interior = domain.interior + else: + # Find the interior corresponding to this patch + domain_interior = None + for interior in domain.interior.args: + if interior.name == patch_name: + domain_interior = interior + break + if domain_interior is None: + continue # Skip validation if interior not found + + # Validate grid length vs ncells consistency and other checks + if isinstance(patch_ncells, (list, tuple)): + if grid_dim == 1: + # For 1D case, grid is either [points] or [[points]] + if isinstance(patch_grid[0], (list, tuple, np.ndarray)): + grid_points = patch_grid[0] # [[points]] format + else: + grid_points = patch_grid # [points] format + + expected_length = patch_ncells[0] + 1 + if len(grid_points) != expected_length: + raise ValueError(f"Dimension 0: grid length ({len(grid_points)}) must be ncells+1 ({expected_length})") + + # Check if grid points are strictly increasing + grid_array = np.asarray(grid_points) + if not np.all(np.diff(grid_array) > 0): + raise ValueError(f"grid points must be strictly increasing") + + # Check if grid boundaries match domain boundaries + if domain_interior is not None: + domain_min = domain_interior.min_coords[0] + domain_max = domain_interior.max_coords[0] + tolerance = 1e-14 # Numerical tolerance for boundary comparison + + if abs(grid_array[0] - domain_min) > tolerance: + raise ValueError(f"grid start ({grid_array[0]}) must match domain minimum ({domain_min})") + + if abs(grid_array[-1] - domain_max) > tolerance: + raise ValueError(f"grid end ({grid_array[-1]}) must match domain maximum ({domain_max})") + else: + # For multi-dimensional case + for dim in range(grid_dim): + grid_points = patch_grid[dim] + expected_length = patch_ncells[dim] + 1 + if len(grid_points) != expected_length: + raise ValueError(f"Dimension {dim}: grid length ({len(grid_points)}) must be ncells+1 ({expected_length})") + + # Check if grid points are strictly increasing + grid_array = np.asarray(grid_points) + if not np.all(np.diff(grid_array) > 0): + raise ValueError(f"grid points must be strictly increasing") + + # Check if grid boundaries match domain boundaries + if domain_interior is not None: + domain_min = domain_interior.min_coords[dim] + domain_max = domain_interior.max_coords[dim] + tolerance = 1e-14 # Numerical tolerance for boundary comparison + + if abs(grid_array[0] - domain_min) > tolerance: + raise ValueError(f"grid start ({grid_array[0]}) must match domain minimum ({domain_min})") + + if abs(grid_array[-1] - domain_max) > tolerance: + raise ValueError(f"grid end ({grid_array[-1]}) must match domain maximum ({domain_max})") #-------------------------------------------------------------------------- @property @@ -220,6 +480,10 @@ def is_parallel(self): def mappings(self): return self._mappings + @property + def grid(self): + return self._grid + def __len__(self): return len(self.domain) diff --git a/psydac/cad/tests/test_geometry.py b/psydac/cad/tests/test_geometry.py index 0f0f1a93d..61b6d8b9b 100644 --- a/psydac/cad/tests/test_geometry.py +++ b/psydac/cad/tests/test_geometry.py @@ -270,6 +270,275 @@ def test_geometry_1(): square = Geometry.as_square(ncells=[10, 10]) cube = Geometry.as_cube(ncells=[10, 10, 10]) +#============================================================================== +def test_geometry_grid_parameter_1d(): + """Test that the grid parameter is correctly stored and accessed in 1D geometry.""" + + + # Create a Line domain + domain = Line('Omega', bounds=(0.0, 1.0)) + + # Create a 1D geometry with custom grid + custom_grid = {domain.name: np.array([0.0, 0.2, 0.5, 0.8, 1.0])} + + # Define ncells and mappings + ncells = {domain.name: [4]} # 4 cells from the 5 breakpoints + mappings = {domain.name: None} # No mapping for topological domain + + # Create geometry with custom grid + geo = Geometry(domain=domain, ncells=ncells, mappings=mappings, grid=custom_grid) + + # Test that grid is correctly stored + assert geo.grid is not None + assert len(geo.grid) == 1 + assert domain.name in geo.grid + assert np.allclose(geo.grid[domain.name], custom_grid[domain.name]) + +#============================================================================== +def test_geometry_grid_parameter_2d(): + """Test that the grid parameter is correctly stored and accessed in 2D geometry.""" + + # Create a 2D geometry with custom grid + custom_grid = [ + np.array([0.0, 0.3, 0.7, 1.0]), # 3 cells in first direction + np.array([0.0, 0.25, 0.5, 0.75, 1.0]) # 4 cells in second direction + ] + + # Create a Square domain + domain = Square('Omega', bounds1=(0.0, 1.0), bounds2=(0.0, 1.0)) + + # Define ncells + ncells = [3, 4] + + # Create geometry with custom grid using from_topological_domain + geo = Geometry.from_topological_domain(domain, ncells=ncells, grid=custom_grid) + + # Test that grid is correctly stored + assert geo.grid is not None + assert len(geo.grid) == 1 # One patch named 'Omega' + assert domain.name in geo.grid + assert len(geo.grid[domain.name]) == 2 # Two dimensions + assert np.allclose(geo.grid[domain.name][0], custom_grid[0]) + assert np.allclose(geo.grid[domain.name][1], custom_grid[1]) + +#============================================================================== +def test_geometry_grid_parameter_3d(): + """Test that the grid parameter is correctly stored and accessed in 3D geometry.""" + + # Create a 3D geometry with custom grid + custom_grid_list = [ + np.array([0.0, 0.3, 0.7, 1.0]), # 3 cells in first direction + np.array([0.0, 0.2, 0.5, 0.8, 1.0]), # 4 cells in second direction + np.array([0.0, 0.4, 1.0]) # 2 cells in third direction + ] + + # Create a Cube domain + domain = Cube('Omega', bounds1=(0.0, 1.0), bounds2=(0.0, 1.0), bounds3=(0.0, 1.0)) + + # Define ncells and convert grid to dictionary format + ncells = {domain.name: [3, 4, 2]} + mappings = {domain.name: None} + custom_grid = {domain.name: custom_grid_list} + + # Create geometry with custom grid + geo = Geometry(domain=domain, ncells=ncells, mappings=mappings, grid=custom_grid) + + # Test that grid is correctly stored + assert geo.grid is not None + assert len(geo.grid) == 1 # One patch named 'Omega' + assert domain.name in geo.grid + assert len(geo.grid[domain.name]) == 3 # Three dimensions + assert np.allclose(geo.grid[domain.name][0], custom_grid_list[0]) + assert np.allclose(geo.grid[domain.name][1], custom_grid_list[1]) + assert np.allclose(geo.grid[domain.name][2], custom_grid_list[2]) + +#============================================================================== +def test_geometry_grid_parameter_none(): + """Test that when grid is None, the geometry still works correctly.""" + + # Create a Square domain + domain = Square('Omega', bounds1=(0.0, 1.0), bounds2=(0.0, 1.0)) + + # Define ncells and mappings + ncells = {domain.name: [4, 4]} + mappings = {domain.name: None} + + # Create geometry without custom grid (default None) + geo = Geometry(domain=domain, ncells=ncells, mappings=mappings, grid=None) + + # Test that grid is None + assert geo.grid is None + +#============================================================================== +def test_geometry_grid_parameter_with_mapping(): + """Test grid parameter with an actual mapping.""" + + # Create a simple identity mapping + ncells = [3, 3] + degree = [2, 2] + mapping = discrete_mapping('identity', ncells=ncells, degree=degree) + + # Create custom grid matching the ncells + custom_grid_list = [ + np.array([0.0, 0.4, 0.7, 1.0]), # 3 cells + np.array([0.0, 0.3, 0.6, 1.0]) # 3 cells + ] + + # Create a topological domain + F = Mapping('F', dim=2) + domain = F(Square(name='Omega')) + + # Associate the mapping to the topological domain + mappings = {domain.name: mapping} + ncells_dict = {domain.name: ncells} + custom_grid = {domain.name: custom_grid_list} + + # Create geometry with both mapping and custom grid + geo = Geometry(domain=domain, ncells=ncells_dict, mappings=mappings, grid=custom_grid) + + # Test that grid is correctly stored + assert geo.grid is not None + assert len(geo.grid) == 1 # One patch + assert domain.name in geo.grid + assert len(geo.grid[domain.name]) == 2 # Two dimensions + assert np.allclose(geo.grid[domain.name][0], custom_grid_list[0]) + assert np.allclose(geo.grid[domain.name][1], custom_grid_list[1]) + + # Test that mapping is also correctly stored + assert geo.mappings is not None + assert domain.name in geo.mappings + assert geo.mappings[domain.name] is mapping + +#============================================================================== +def test_geometry_grid_parameter_3d_with_mapping(): + """Test grid parameter with topological cube domain in 3D.""" + + # Create custom grid for 3D + custom_grid_list = [ + np.array([0.0, 0.5, 1.0]), # 2 cells + np.array([0.0, 0.3, 0.6, 1.0]), # 3 cells + np.array([0.0, 0.7, 1.0]) # 2 cells + ] + + # Create a topological cube domain + domain = Cube('Omega', bounds1=(0.0, 1.0), bounds2=(0.0, 1.0), bounds3=(0.0, 1.0)) + + # Define ncells and mappings + ncells = {domain.name: [2, 3, 2]} + mappings = {domain.name: None} # No mapping for topological domain + custom_grid = {domain.name: custom_grid_list} + + # Create geometry with custom grid + geo = Geometry(domain=domain, ncells=ncells, mappings=mappings, grid=custom_grid) + + # Test that grid is correctly stored + assert geo.grid is not None + assert len(geo.grid) == 1 # One patch named 'Omega' + assert domain.name in geo.grid + assert len(geo.grid[domain.name]) == 3 # Three dimensions + assert np.allclose(geo.grid[domain.name][0], custom_grid_list[0]) + assert np.allclose(geo.grid[domain.name][1], custom_grid_list[1]) + assert np.allclose(geo.grid[domain.name][2], custom_grid_list[2]) + +#============================================================================== +def test_geometry_grid_parameter_3d_non_uniform(): + """Test grid parameter with non-uniform spacing in 3D.""" + + # Create a 3D geometry with highly non-uniform grid + custom_grid_list = [ + np.array([0.0, 0.1, 0.15, 0.2, 0.8, 0.9, 1.0]), # 6 cells, clustered + np.array([0.0, 0.6, 0.65, 1.0]), # 3 cells, mostly at the end + np.array([0.0, 0.05, 0.1, 0.85, 0.95, 1.0]) # 5 cells, clustered at both ends + ] + + # Create a Cube domain + domain = Cube('Omega', bounds1=(0.0, 1.0), bounds2=(0.0, 1.0), bounds3=(0.0, 1.0)) + + # Define ncells and mappings + ncells = {domain.name: [6, 3, 5]} + mappings = {domain.name: None} + custom_grid = {domain.name: custom_grid_list} + + # Create geometry with custom grid + geo = Geometry(domain=domain, ncells=ncells, mappings=mappings, grid=custom_grid) + + # Test that grid is correctly stored + assert geo.grid is not None + assert len(geo.grid) == 1 # One patch named 'Omega' + assert domain.name in geo.grid + assert len(geo.grid[domain.name]) == 3 # Three dimensions + assert np.allclose(geo.grid[domain.name][0], custom_grid_list[0]) + assert np.allclose(geo.grid[domain.name][1], custom_grid_list[1]) + assert np.allclose(geo.grid[domain.name][2], custom_grid_list[2]) + +#============================================================================== +def test_geometry_grid_parameter_topological(): + """Test grid parameter with topological domain using regular constructor.""" + + # Create a square domain + domain = Square('Omega', bounds1=(0.0, 1.0), bounds2=(0.0, 1.0)) + ncells = [4, 5] + + # Create custom grid + custom_grid_list = [ + np.array([0.0, 0.2, 0.5, 0.8, 1.0]), # 4 cells + np.array([0.0, 0.1, 0.3, 0.6, 0.9, 1.0]) # 5 cells + ] + + # Create geometry from topological domain with custom grid + ncells_dict = {domain.name: ncells} + mappings = {domain.name: None} + custom_grid = {domain.name: custom_grid_list} + + geo = Geometry(domain=domain, ncells=ncells_dict, mappings=mappings, grid=custom_grid) + + # Test that grid is correctly stored + assert geo.grid is not None + assert len(geo.grid) == 1 # One patch named 'Omega' + assert domain.name in geo.grid + assert len(geo.grid[domain.name]) == 2 # Two dimensions + assert np.allclose(geo.grid[domain.name][0], custom_grid_list[0]) + assert np.allclose(geo.grid[domain.name][1], custom_grid_list[1]) + +#============================================================================== +def test_geometry_grid_parameter_export_import(): + """Test that grid parameter is preserved through export/import cycle.""" + + # Create a geometry with custom grid and identity mapping + ncells = [3, 3] + degree = [2, 2] + mapping = discrete_mapping('identity', ncells=ncells, degree=degree) + + custom_grid_list = [ + np.array([0.0, 0.3, 0.7, 1.0]), + np.array([0.0, 0.4, 0.8, 1.0]) + ] + + # Create topological domain + F = Mapping('F', dim=2) + domain = F(Square(name='Omega')) + + mappings = {domain.name: mapping} + ncells_dict = {domain.name: ncells} + custom_grid = {domain.name: custom_grid_list} + + # Create original geometry with grid + geo_original = Geometry(domain=domain, ncells=ncells_dict, mappings=mappings, grid=custom_grid) + + # Export geometry + filename = 'test_grid_geometry.h5' + geo_original.export(filename) + + # Import geometry + geo_imported = Geometry(filename=filename) + + # Note: The current implementation might not preserve the grid parameter in HDF5 files + # This test documents the current behavior and can be updated when grid export/import is implemented + + # Clean up + if os.path.exists(filename): + os.remove(filename) + #============================================================================== # CLEAN UP SYMPY NAMESPACE #============================================================================== @@ -289,7 +558,8 @@ def teardown_module(): 'quart_circle_1.h5', 'circle.h5', 'pipe.h5', - 'L_shaped.h5' + 'L_shaped.h5', + 'test_grid_geometry.h5' ] for fname in filenames: if os.path.exists(fname):