Add Custom Grid Parameter to Geometry Class and Discretization #518

ArasuCandassamy · 2025-08-01T11:54:43Z

This pull request introduces support for custom grids in the Geometry class and updates the discretize_space function to handle such grids. It also includes extensive testing to ensure the new functionality works correctly. The most important changes are grouped as follows:

Enhancements to `Geometry` class:

Added a grid parameter to the Geometry class, allowing users to specify custom grids of breakpoints for each dimension. This parameter is stored as a private attribute _grid and exposed via a new grid property. (psydac/cad/geometry.py: [1] [2] [3] [4])

Updates to `discretize_space` function:

Modified the discretize_space function to handle custom grids provided via the domain_h object. This includes safety checks for grid consistency, such as matching dimensions, grid length, sorting, and boundary alignment with domain limits. (psydac/api/discretization.py: psydac/api/discretization.pyR463-R493)

Testing for custom grids:

Added multiple unit tests to validate the behavior of the grid parameter in various scenarios, including 1D, 2D, and 3D geometries, geometries with mappings, and non-uniform grids. Tests also verify that the grid is correctly stored and accessed. (psydac/cad/tests/test_geometry.py: [1] [2]

Testing discretization with a custom grids:

Added multiple unit tests to validate the discretization with a custom grid in various scenarios. Tests also verify the consistency between the given number of cells and the given custom grid. (psydac/api/tests/test_discretization_grid.py: [1]

Exemple of discretization using a custom grid:

import numpy as np
from sympde.topology import VectorFunctionSpace, Square
from psydac.api.discretization import discretize
from psydac.cad.geometry import Geometry

domain = Square('domain', (0, 1), (0, 1))

degree = [2,2]
ncells = [8,8]

# Manual grid with a breakpoint at y=0.8
X_grid = np.linspace(0, 1, ncells[0] + 1)  # uniform grid in x
Y_grid = np.linspace(0, 1, ncells[1]) 
Y_grid = np.sort(np.unique(np.concatenate([np.linspace(0, 1, ncells[1]),[0.8]]))) 

custom_grid = [
    X_grid, 
    Y_grid 
]

# Create a Geometry object with the custom grid
domain_h = Geometry(domain=domain, ncells={domain.name: ncells}, mappings={domain.name: None}, grid=grid)

V1 = VectorFunctionSpace('V1', domain) 
V1h = discretize(V1, domain_h, degree=degree)

Example of custom grid in 2D and 3D:

- Add grid parameter to Geometry class constructor and properties - Update discretize_space function to use custom grid from domain_h when available - Fallback to uniform grid when no custom grid is provided - Support for 1D and 2D domains with custom breakpoint grids

…en grid

- Add 15 new security validation tests for grid parameter - Test grid type validation (list/tuple required) - Test grid dimension consistency with domain - Test grid length consistency with ncells - Test monotonic ordering of grid points - Test boundary value matching with tolerance - Add parametrized tests for all error conditions - Test security in 1D, 2D, and 3D configurations - Ensure security checks don't interfere with valid functionality - Total tests: 39 (24 existing + 15 new security tests) These tests ensure robust validation of custom grids preventing: - Invalid grid types and formats - Dimension mismatches between grid and domain - Inconsistent ncells vs grid length - Non-monotonic grid points - Boundary value mismatches - Numerical precision issues All tests pass successfully, confirming grid security implementation works correctly while preserving existing functionality.

codacy-production · 2025-08-01T12:30:19Z

Coverage summary from Codacy

See diff coverage on Codacy

Coverage variation	Diff coverage
✅ +0.06%	✅ 78.86%

Coverage variation details

	Coverable lines	Covered lines	Coverage
Common ancestor commit (`42b214e`)	30762	18640	60.59%
Head commit (`6e90be7`)	61754 (+30992)	37458 (+18818)	60.66% (+0.06%)

Coverage variation is the difference between the coverage for the head and common ancestor commits of the pull request branch: <coverage of head commit> - <coverage of common ancestor commit>

Diff coverage details

	Coverable lines	Covered lines	Diff coverage
Pull request (#518)	123	97	78.86%

Diff coverage is the percentage of lines that are covered by tests out of the coverable lines that the pull request added or modified: <covered lines added or modified>/<coverable lines added or modified> * 100%

See your quality gate settings Change summary preferences

FrederikSchnack · 2025-08-04T09:25:54Z

Hi Arasu, thanks for the nice changes! I had a quick look and noticed one thing:
Most (more modern) examples in psydac that I came across were using the discretize interface to discretize a domain (in order to obtain a Geometry object) instead of calling the constructor directly. For example in

psydac/psydac/api/tests/test_2d_multipatch_poisson.py

Line 63 in e47ff5d

domain_h = discretize(domain, ncells=ncells)

It would be very nice, if you could add the grid keyword to discretize_domain in

psydac/psydac/api/discretization.py

Line 556 in e47ff5d

    
           def discretize_domain(domain, *, filename=None, ncells=None, periodic=None, comm=None, mpi_dims_mask=None):

as well (and the function it is calling).

Unrelated to this PR, I just realized that Geometry takes mapping(s) as arguments, which we also do not provide in discretize_domain. Do you know why that is @yguclu? Because I was always using the property domain_h.domain.mapping in the multipatch feec interface. Maybe we can add this here as well.

- Add grid parameter to discretize_domain() in discretization.py - Add grid parameter to Geometry.from_topological_domain() in geometry.py

ArasuCandassamy · 2025-08-04T12:37:12Z

Thank you @FrederikSchnack for your comments.

I've added the grid keyword to discretize_domain in

psydac/psydac/api/discretization.py

Line 556 in e79a11c

    
           def discretize_domain(domain, *, filename=None, ncells=None, periodic=None, comm=None, mpi_dims_mask=None, grid=None):

and to Geometry._from_topological_domain in

psydac/psydac/cad/geometry.py

Line 164 in e79a11c

    
           def from_topological_domain(cls, domain, ncells, *, periodic=None, comm=None, mpi_dims_mask=None, grid=None):

I’ve also updated the test file accordingly.
Now, to discretize a domain with a custom grid, you can proceed as follows:

import numpy as np
from sympde.topology import VectorFunctionSpace, Square
from psydac.api.discretization import discretize
from psydac.cad.geometry import Geometry

domain = Square('domain', (0, 1), (0, 1))

degree = [2,2]
ncells = [8,8]

# Manual grid with a breakpoint at y=0.8
X_grid = np.linspace(0, 1, ncells[0] + 1)  # uniform grid in x
Y_grid = np.sort(np.unique(np.concatenate([np.linspace(0, 1, ncells[1]),[0.8]]))) 

custom_grid = [
    X_grid, 
    Y_grid 
]

# Create a Geometry object with the custom grid
domain_h = discretize(domain=domain, ncells=ncells, grid=grid, comm=None)

V1 = VectorFunctionSpace('V1', domain) 
V1h = discretize(V1, domain_h, degree=degree)

yguclu · 2025-08-04T16:59:19Z

I’ve also updated the test file accordingly. Now, to discretize a domain with a custom grid, you can proceed as follows:

import numpy as np
from sympde.topology import VectorFunctionSpace, Square
from psydac.api.discretization import discretize
from psydac.cad.geometry import Geometry

domain = Square('domain', (0, 1), (0, 1))

degree = [2,2]
ncells = [8,8]

# Manual grid with a breakpoint at y=0.8
X_grid = np.linspace(0, 1, ncells[0] + 1)  # uniform grid in x
Y_grid = np.sort(np.unique(np.concatenate([np.linspace(0, 1, ncells[1]),[0.8]]))) 

custom_grid = [
    X_grid, 
    Y_grid 
]

# Create a Geometry object with the custom grid
domain_h = discretize(domain=domain, ncells=ncells, grid=grid, comm=None)

V1 = VectorFunctionSpace('V1', domain) 
V1h = discretize(V1, domain_h, degree=degree)

Do you mean grid=custom_grid?

Is there an option to only pass the grid along one axis, e.g. with grid = [None, Y_grid]?

yguclu · 2025-08-04T17:12:49Z

Unrelated to this PR, I just realized that Geometry takes mapping(s) as arguments, which we also do not provide in discretize_domain. Do you know why that is @yguclu? Because I was always using the property domain_h.domain.mapping in the multipatch feec interface. Maybe we can add this here as well.

Psydac handles spline/NURBS (discrete) mappings quite differently from analytical mappings. Discrete mappings are handled through the Geometry constructor if I am not mistaken, by reading the information from an input file. Analytical mappings are handled through the so-called "mapped domains", which are obtained by applying the mapping to a logical domain. The mapped domain is then passed to discretize_domain through the domain argument.

yguclu

Nice work @ArasuCandassamy!

I have a few comments which you will find below. I was also thinking that it would be nice to allow the user to specify the grid just for some of the directions. For example, this could be done with the call

grid_y = [0, 0.1, 0.2, 0.333, 0.444, 0.5]
domain_h = discretize(domain, ncells=[4, 5, 6], grid=[None, grid_y, None])

where grid=[None, grid_y, None] means that the grid is only specified along the axis of index 1.

yguclu · 2025-08-04T17:23:25Z

psydac/api/discretization.py

+                if hasattr(domain_h, 'grid') and domain_h.grid is not None:
+                    # Use provided grid of breakpoints
+                    grids = domain_h.grid
+
+                    # Safety checks for grid consistency
+                    if not isinstance(grids, (list, tuple)):
+                        raise TypeError(f"Grid must be a list or tuple, got {type(grids)}")
+
+                    if len(grids) != len(ncells):
+                        raise ValueError(f"Grid dimensions ({len(grids)}) must match domain dimensions ({len(ncells)})")
+
+                    for dim, (grid, nc, xmin, xmax) in enumerate(zip(grids, ncells, min_coords, max_coords)):
+                        grid = np.asarray(grid)
+
+                        # Check grid length vs ncells consistency
+                        expected_grid_length = nc + 1
+                        if len(grid) != expected_grid_length:
+                            raise ValueError(f"Dimension {dim}: grid length ({len(grid)}) must be ncells+1 ({expected_grid_length})")
+
+                        # Check if grid is sorted
+                        if not np.all(np.diff(grid) > 0):
+                            raise ValueError(f"Dimension {dim}: grid points must be strictly increasing")
+
+                        # Check domain boundaries (with tolerance for numerical precision)
+                        tol = 1e-12
+                        if abs(grid[0] - xmin) > tol:
+                            raise ValueError(f"Dimension {dim}: grid start ({grid[0]}) must match domain minimum ({xmin})")
+                        if abs(grid[-1] - xmax) > tol:
+                            raise ValueError(f"Dimension {dim}: grid end ({grid[-1]}) must match domain maximum ({xmax})")


I think that these checks should be made directly in the constructor of the Geometry class, not here. Here instead you should check that if the knots are passed, these are consistent with domain_h.ncells and domain_h.grid.

I have removed these checks from there and I have put them in geometry.py with a function in the Geometry class:

psydac/psydac/cad/geometry.py

Line 338 in 6e90be7

def _validate_grid_consistency(self, domain, ncells, grid):

yguclu · 2025-08-04T17:34:39Z

psydac/cad/geometry.py

            self._mappings = mappings
            self._cart     = None
            self._is_parallel = comm is not None
+            self._grid = grid


Before storing the input parameter grid, there are a few things that we should verify:

If filename is provided, one cannot pass ncells or grid. I would avoid more complicated solutions.

If ncells or grid is provided, one cannot pass filename.

Since ncells and grid are dictionaries, they must have for keys the names of the domain patches, i.e. domain.keys()

I've seen that the verification is done in the discretize_domain function here :

psydac/psydac/api/discretization.py

Lines 556 to 571 in e79a11c

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'")

elif filename and ncells:

raise ValueError("Cannot provide both 'filename' and 'ncells'")

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, grid=grid)

Can I add the checks on grid, filenameand ncells here ?

Yes, I know that the same checks are done in discretize_domain...😅 Nevertheless a class constructor should always check its arguments, and now is a good time to fix it!😉

yguclu · 2025-08-04T17:38:12Z

psydac/cad/geometry.py

+    grid: list of breakpoints
+        The grid of breakpoints in each direction.
+        If not given, the grid will be constructed from the ncells and periodicity.


This argument has to be a dictionary just like ncells, periodic, and mappings. I think the docstring has to be updated regarding ncells and periodic: neither a list nor a tuple is accepted.

The new docstring is:

grid: list | tuple | dict The grid of breakpoints in each direction. If not given, the grid will be constructed from the ncells and periodicity.

So, it's possible to give a list or a tuple to create a Geometry class object, but Geometry.from_topological_domain is converting tuples and lists into a dict:

psydac/psydac/cad/geometry.py

Lines 229 to 230 in 6e90be7

if isinstance(grid, (list, tuple)):

grid = {itr.name:grid for itr in interior}

yguclu · 2025-08-04T17:55:44Z

psydac/cad/geometry.py

    @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):


I am sorry that this function does not have a docstring. In practice it is a helper constructor which discretizes a multipatch domain using the same discretization parameters for all patches, i.e. the same number of cells, periodicity, and (now) grid.

I've added a docstring for this function:

psydac/psydac/cad/geometry.py

Lines 174 to 207 in 6e90be7

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.

"""

- Add consistency checks in geometry.py - Pass grid as a dictionary - Add docstrings - Add support for partially non-uniform grids - Modify and add tests for new features

ArasuCandassamy added 5 commits July 31, 2025 10:14

Bug fix

8edcd2b

Add security checks and verify consistency between ncells and the giv…

f0f73a3

…en grid

Update tests with grid feature in Geometry class

a21d841

ArasuCandassamy requested a review from yguclu August 4, 2025 08:31

Add grid parameter support to discretize_domain function

e79a11c

- Add grid parameter to discretize_domain() in discretization.py - Add grid parameter to Geometry.from_topological_domain() in geometry.py

yguclu requested changes Aug 4, 2025

View reviewed changes

ArasuCandassamy added 3 commits August 5, 2025 09:23

Docstring for Geometry.from_topological_domain

8d15fff

Pass grid as a dictionary

bd8d09a

Apply requested changes

6e90be7

- Add consistency checks in geometry.py - Pass grid as a dictionary - Add docstrings - Add support for partially non-uniform grids - Modify and add tests for new features

This was referenced Oct 15, 2025

Fix method plot_2d_decomposition of TensorFemSpace #529

Open

Make plot_2d_decomposition use Geometry instead of TensorFemSpace #532

Open

	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'")

	elif filename and ncells:
	raise ValueError("Cannot provide both 'filename' and 'ncells'")

	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, grid=grid)

	if isinstance(grid, (list, tuple)):
	grid = {itr.name:grid for itr in interior}

	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.
	"""

Add Custom Grid Parameter to Geometry Class and Discretization #518

Are you sure you want to change the base?

Add Custom Grid Parameter to Geometry Class and Discretization #518

Uh oh!

Conversation

ArasuCandassamy commented Aug 1, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Enhancements to Geometry class:

Updates to discretize_space function:

Testing for custom grids:

Testing discretization with a custom grids:

Exemple of discretization using a custom grid:

Example of custom grid in 2D and 3D:

Uh oh!

codacy-production bot commented Aug 1, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Coverage summary from Codacy

Uh oh!

FrederikSchnack commented Aug 4, 2025

Uh oh!

ArasuCandassamy commented Aug 4, 2025

Uh oh!

yguclu commented Aug 4, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

yguclu commented Aug 4, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

yguclu left a comment • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

ArasuCandassamy commented Aug 1, 2025 •

edited

Loading

Enhancements to `Geometry` class:

Updates to `discretize_space` function:

codacy-production bot commented Aug 1, 2025 •

edited

Loading

yguclu commented Aug 4, 2025 •

edited

Loading

yguclu commented Aug 4, 2025 •

edited

Loading

yguclu left a comment •

edited

Loading