diff --git a/.github/workflows/testing.yml b/.github/workflows/testing.yml
index d8ff30b1c..c3c2da8a4 100644
--- a/.github/workflows/testing.yml
+++ b/.github/workflows/testing.yml
@@ -195,7 +195,7 @@ jobs:
       - name: Run coverage tests on macOS
         if: matrix.os == 'macos-14'
         working-directory: ./pytest
-        run: >- 
+        run: >-
           python -m pytest -n auto
           --cov psydac
           --cov-config $GITHUB_WORKSPACE/pyproject.toml
diff --git a/examples/poisson_2d_mapping.py b/examples/poisson_2d_mapping.py
index 74fb620d3..79e84b843 100644
--- a/examples/poisson_2d_mapping.py
+++ b/examples/poisson_2d_mapping.py
@@ -781,7 +781,8 @@ def main(*, test_case, ncells, degree, nquads,
     ##########
 
     # Plot domain decomposition (master only)
-    V.plot_2d_decomposition(model.mapping, refine=N)
+    fig = V.plot_2d_decomposition(mapping, refine=N)
+    fig.show()
 
     # Perform other visualization using master or all processes
     if not distribute_viz:
diff --git a/psydac/fem/tensor.py b/psydac/fem/tensor.py
index 644f08cc4..9bc6678e4 100644
--- a/psydac/fem/tensor.py
+++ b/psydac/fem/tensor.py
@@ -1184,24 +1184,77 @@ def set_refined_space(self, ncells, new_space):
         self._refined_space[tuple(ncells)] = new_space
 
     # ...
-    def plot_2d_decomposition(self, mapping=None, refine=10):
+    def plot_2d_decomposition(self, mapping=None, *, refine=10, fig=None, ax=None, mpi_root=0):
+        """
+        Plot decomposition of 2D TensorFemSpace w/ mapping to 2D physical space
+
+        Plot the domain decomposition across MPI processes of a 2D
+        TensorFemSpace with a mapping between 2D logical and 2D physical spaces.
+        This function must be called collectively, and only the root process will make
+        the plot. On non-root processes the arguments `fig` and `ax` must be None.
+
+        Parameters
+        ----------
+        mapping : BasicCallableMapping
+            Mapping from (eta1, eta2) to (x1, x2).
+
+        refine : int, default=10
+            Cell refinement along the logical dimensions eta1 and eta2.
+
+        fig : plt.Figure, optional
+            Figure where the plot should be made. Must be None on non-root processes.
 
+        ax : plt.Axes, optional
+            Axes where the plot should be made. Must be None on non-root processes.
+
+        mpi_root: int, default=0
+            The rank of the MPI root process which should create the plot.
+
+        Returns
+        -------
+        plt.Figure
+            Figure where the plot was made. Coincides with `fig` if provided.
+        """
         import matplotlib.pyplot as plt
         from matplotlib.patches  import Polygon, Patch
+        from sympde.topology.mapping import BasicCallableMapping
         from psydac.utilities.utils import refine_array_1d
 
+        # Sanity check
+        assert self.ldim == 2, "Function only works in 2D"
+
+        # Check mapping
         if mapping is None:
             mapping = lambda eta: eta
         else:
-            assert mapping.ldim == self.ldim == 2
-            assert mapping.pdim == self.ldim == 2
+            assert isinstance(mapping, BasicCallableMapping)
+            assert mapping.ldim == 2, "Domain of argument `mapping` must be 2D"
+            assert mapping.pdim == 2, "Codomain of argument `mapping` must be 2D"
 
-        assert refine >= 1
-        N = int(refine)
-        V1, V2 = self.spaces
+        # Check refine argument
+        assert isinstance(refine, int), f"Argument `refine` must be int, got {type(refine)} instead"
+        assert refine >= 1, f"Argument `refine` must be >= 1, got {refine} instead"
 
+        # Extract information about MPI communicator
         mpi_comm = self.coeff_space.cart.comm
         mpi_rank = mpi_comm.rank
+        mpi_size = mpi_comm.size
+
+        # Check mpi_root argument
+        assert isinstance(mpi_root, int), f"Argument `mpi_root` must be int, got {type(mpi_root)} instead"
+        assert mpi_root >= 0, f"Argument `mpi_root` must be >= 0, got {mpi_root} instead"
+        assert mpi_root < mpi_size, f"Argument `mpi_root` must be smaller than communicator size ({mpi_size}), got {mpi_root} instead"
+
+        # Check fig and ax arguments
+        if mpi_rank == mpi_root:
+            assert isinstance(fig, plt.Figure) or fig is None, f"Argument `fig` must be matplotlib Figure, got {type(fig)} instead"
+            assert isinstance(ax, plt.Axes) or ax is None, f"Argument `ax` must be matplotlib Axes, got {type(ax)} instead"
+        else:
+            assert fig is None, f"Argument `fig` must be None on non-root process with rank {mpi_rank}"
+            assert ax is None, f"Argument `ax` must be None on non-root process with rank {mpi_rank}"
+
+        N = refine
+        V1, V2 = self.spaces
 
         # Local grid, refined
         [sk1, sk2], [ek1, ek2] = self.local_domain
@@ -1218,23 +1271,62 @@ def plot_2d_decomposition(self, mapping=None, refine=10):
         poly = Polygon(xy, edgecolor='None')
 
         # Gather polygons on master process
-        polys = mpi_comm.gather(poly)
+        polys = mpi_comm.gather(poly, root=mpi_root)
+
+        # Gather (s1, s2, e1, e2) on root
+        if mpi_rank == mpi_root:
+            s1_all = np.empty(mpi_size, dtype=int)
+            s2_all = np.empty(mpi_size, dtype=int)
+            e1_all = np.empty(mpi_size, dtype=int)
+            e2_all = np.empty(mpi_size, dtype=int)
+        else:
+            s1_all = None
+            s2_all = None
+            e1_all = None
+            e2_all = None
+
+        mpi_comm.Gather(sk1 * N, s1_all, root=mpi_root)
+        mpi_comm.Gather(sk2 * N, s2_all, root=mpi_root)
+        mpi_comm.Gather((ek1 + 1) * N, e1_all, root=mpi_root)
+        mpi_comm.Gather((ek2 + 1) * N, e2_all, root=mpi_root)
+
+        # Gather pcoords on root
+        # TODO: use Gatherv, and NumPy arrays as buffers
+        gathered_pcoords = mpi_comm.gather(pcoords, root=mpi_root)
 
         #-------------------------------
         # Non-master processes stop here
-        if mpi_rank != 0:
+        if mpi_rank != mpi_root:
             return
         #-------------------------------
 
-        # Global grid, refined
-        eta1    = refine_array_1d(V1.breaks, N)
-        eta2    = refine_array_1d(V2.breaks, N)
-        pcoords = np.array([[mapping(e1, e2) for e2 in eta2] for e1 in eta1])
-        xx      = pcoords[:, :, 0]
-        yy      = pcoords[:, :, 1]
+        # Reconstruct global grid (refined) on root process
+        global_shape   = ((V1.breaks.size - 1) * N + 1,
+                          (V2.breaks.size - 1) * N + 1,
+                          2)
+        pcoords_global = np.empty(global_shape)
+
+        for rank in range(mpi_comm.size):
+            s1 = s1_all[rank]
+            e1 = e1_all[rank]
+            s2 = s2_all[rank]
+            e2 = e2_all[rank]
+            pcoords_global[s1:e1+1, s2:e2+1, :] = gathered_pcoords[rank]
+
+        xx = pcoords_global[:, :, 0]
+        yy = pcoords_global[:, :, 1]
+
+        # If fig or ax are given, get one from the other. Otherwise create new ones
+        if fig and ax:
+            assert ax in fig.axes, "Argument `ax` must be in `fig.axes`"
+        elif fig:
+            ax = fig.gca()
+        elif ax:
+            fig = ax.figure
+        else:
+            fig, ax = plt.subplots(1, 1)
 
         # Plot decomposed domain
-        fig, ax = plt.subplots(1, 1)
         colors  = itertools.cycle(plt.rcParams['axes.prop_cycle'].by_key()['color'])
         handles = []
         for i, (poly, color) in enumerate(zip(polys, colors)):
@@ -1253,7 +1345,8 @@ def plot_2d_decomposition(self, mapping=None, refine=10):
         ax.set_aspect('equal')
         ax.legend(handles=handles, bbox_to_anchor=(1.05, 1), loc=2)
         fig.tight_layout()
-        fig.show()
+
+        return fig
 
     # ...
     def __str__(self):
diff --git a/psydac/fem/tests/data/decomp_analytical_1_procs.png b/psydac/fem/tests/data/decomp_analytical_1_procs.png
new file mode 100644
index 000000000..f4066700f
Binary files /dev/null and b/psydac/fem/tests/data/decomp_analytical_1_procs.png differ
diff --git a/psydac/fem/tests/data/decomp_analytical_4_procs.png b/psydac/fem/tests/data/decomp_analytical_4_procs.png
new file mode 100644
index 000000000..9eedf7529
Binary files /dev/null and b/psydac/fem/tests/data/decomp_analytical_4_procs.png differ
diff --git a/psydac/fem/tests/data/decomp_spline_1_procs.png b/psydac/fem/tests/data/decomp_spline_1_procs.png
new file mode 100644
index 000000000..5e4ffe257
Binary files /dev/null and b/psydac/fem/tests/data/decomp_spline_1_procs.png differ
diff --git a/psydac/fem/tests/data/decomp_spline_4_procs.png b/psydac/fem/tests/data/decomp_spline_4_procs.png
new file mode 100644
index 000000000..6160ff220
Binary files /dev/null and b/psydac/fem/tests/data/decomp_spline_4_procs.png differ
diff --git a/psydac/fem/tests/test_tensor.py b/psydac/fem/tests/test_tensor.py
new file mode 100644
index 000000000..cc0dc8200
--- /dev/null
+++ b/psydac/fem/tests/test_tensor.py
@@ -0,0 +1,220 @@
+import os
+import contextlib
+from pathlib import Path
+
+import pytest
+import numpy as np
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+from PIL import Image
+from mpi4py import MPI
+
+from sympde.topology.domain import Square
+from sympde.topology.space import ScalarFunctionSpace
+from sympde.topology.analytical_mapping import TargetMapping
+
+from psydac.mapping.discrete_gallery import discrete_mapping
+from psydac.api.discretization import discretize
+
+
+#==============================================================================
+# Machinery for comparing a PNG image with a reference on the root MPI process
+#==============================================================================
+
+def similar_images(file1, file2, tolerance=0.01):
+    """
+    Compare two PNG images and check if they are similar within tolerance.
+
+    Parameters
+    ----------
+    file1 : str
+        Path to first PNG file.
+    file2 : str
+        Path to second PNG file.
+    tolerance: float
+        Maximum allowed average difference between pixel values (0-1).
+
+    Returns
+    -------
+    bool :
+        True if images are similar enough, False otherwise.
+    """
+    # Load images and convert to numpy arrays
+    img1 = np.array(Image.open(file1)).astype(float)
+    img2 = np.array(Image.open(file2)).astype(float)
+
+    # Check dimensions match
+    if img1.shape != img2.shape:
+        return False
+
+    # Normalize pixel values to 0-1
+    img1 = img1 / 255.0
+    img2 = img2 / 255.0
+
+    # Calculate mean absolute difference
+    diff = np.mean(np.abs(img1 - img2))
+
+    return diff <= tolerance
+
+
+@contextlib.contextmanager
+def consistent_png_rendering():
+    """
+    Context manager for consistent Matplotlib rendering across platforms.
+    """
+    # Store original settings
+    orig_backend = mpl.get_backend()
+    orig_settings = {
+        'text.usetex': mpl.rcParams['text.usetex'],
+        'font.family': mpl.rcParams['font.family'],
+    }
+
+    try:
+        # Use Agg backend (pure python, no GUI)
+        mpl.use('Agg')
+        # Configure settings for consistent rendering
+        mpl.rcParams.update({
+            'text.usetex': False,
+            'font.family': 'DejaVu Sans',
+        })
+        yield  # Control returns to the with block
+    finally:
+        # Restore original settings
+        mpl.use(orig_backend)
+        mpl.rcParams.update(orig_settings)
+
+
+def compare_figure_to_reference(fig, filename, *, dpi, tol, folder, comm, root):
+    """
+    Compare a matplotlib figure to a reference PNG file on the root MPI process.
+
+    Parameters
+    ----------
+    fig : matplotlib.figure.Figure
+        The figure to compare with the reference image.
+    filename : str
+        Name of the PNG file to save and compare.
+    dpi : int
+        Dots per inch resolution for saving the figure.
+    tol : float
+        Tolerance for image comparison (between 0 and 1).
+    folder : str
+        Name of the folder containing reference images.
+    comm : mpi4py.MPI.Comm
+        MPI communicator object.
+    root : int
+        Rank of the MPI process that should perform the comparison.
+
+    Returns
+    -------
+    bool
+        True if the images are similar enough within the tolerance,
+        False otherwise. The result is broadcast to all MPI processes.
+
+    Notes
+    -----
+    The function saves the figure to a temporary file, compares it with
+    the reference image, and then removes the temporary file. Only the
+    root process performs the actual comparison, but the result is
+    broadcast to all processes.
+    """
+    if comm.rank == root:
+        test_dir = Path(__file__).parent.absolute()
+        file1 = test_dir / filename
+        file2 = test_dir / folder / filename
+        with consistent_png_rendering():
+            fig.savefig(file1, dpi=dpi)
+        close_enough = similar_images(file1, file2, tol)
+        # Clean up the temporary file
+        os.remove(file1)
+    else:
+        close_enough = None
+
+    # Broadcast the boolean result from the root process to all others
+    close_enough = comm.bcast(close_enough, root=root)
+
+    # All MPI processes return the same result
+    return close_enough
+
+#==============================================================================
+# Unit tests
+#==============================================================================
+@pytest.mark.parallel
+@pytest.mark.parametrize('root', ['first', 'last'])
+@pytest.mark.parametrize('kind', ['spline', 'analytical'])
+def test_plot_2d_decomposition(kind, root):
+
+    # MPI communicator
+    mpi_comm = MPI.COMM_WORLD
+    mpi_size = mpi_comm.size
+    mpi_rank = mpi_comm.rank
+
+    # MPI rank which should make the plot
+    if root == 'first':
+        mpi_root = 0
+    elif root == 'last':
+        mpi_root = mpi_size - 1
+    else:
+        raise ValueError(f'root argument has wrong value {root}')
+
+    # Parameters of tensor-product 2D spline space
+    ncells = (6, 9)
+    degree = (2, 2)
+
+    if kind == 'spline':
+        # 2D spline mapping and tensor FEM space (distributed)
+        F, Vh = discrete_mapping('target', ncells=ncells, degree=degree,
+                            comm=mpi_comm, return_space=True)
+    elif kind == 'analytical':
+        Omega = Square('Omega', bounds1=(0, 1), bounds2=(0, 2 * np.pi))
+        params = dict(c1=0, c2=0, k=0.3, D=0.2)
+        M = TargetMapping('M', dim=2, **params)
+        domain = M(Omega)
+        V = ScalarFunctionSpace('V', domain)
+
+        # 2D Geometry object
+        domain_h = discretize(domain, ncells=ncells, periodic=(False, True),
+                              comm=mpi_comm)
+
+        # 2D spline tensor FEM space (distributed)
+        Vh = discretize(V, domain_h, degree=degree)
+
+        # 2D callable mapping (analytical)
+        F = M.get_callable_mapping()
+    else:
+        raise ValueError(f'kind argument has wrong value {kind}')
+
+    # Name of temporary image file to be compared with reference one
+    filename = f'decomp_{kind}_{mpi_size}_procs.png'
+
+    # Relative tolerance for image comparison
+    RTOL = 0.02
+
+    # Plot 2D decomposition
+    # [1] Run without passing (fig, ax)
+    fig = Vh.plot_2d_decomposition(F, refine=5, mpi_root=mpi_root)
+    assert compare_figure_to_reference(fig, filename, folder='data', dpi=100,
+                                       tol=RTOL, comm=mpi_comm, root=mpi_root)
+
+    # [2] Run with given (fig, ax), compatible
+    fig2, ax2 = plt.subplots(1, 1) if mpi_rank == mpi_root else (None, None)
+    Vh.plot_2d_decomposition(F, refine=5, fig=fig2, ax=ax2, mpi_root=mpi_root)
+    assert compare_figure_to_reference(fig2, filename, folder='data', dpi=100,
+                                       tol=RTOL, comm=mpi_comm, root=mpi_root)
+
+    # [3] Run with given (fig, ax), incompatible
+    if mpi_rank == mpi_root:
+        fig3, ax3 = plt.subplots(1, 1)
+        with pytest.raises(AssertionError) as excinfo:
+            Vh.plot_2d_decomposition(F, refine=5, fig=fig2, ax=ax3, mpi_root=mpi_root)
+        assert "Argument `ax` must be in `fig.axes`" in str(excinfo.value)
+        plt.close(fig3)
+    else:
+        Vh.plot_2d_decomposition(F, refine=5, fig=None, ax=None, mpi_root=mpi_root)
+
+#==============================================================================
+if __name__ == '__main__':
+
+    test_plot_2d_decomposition('spline', 'first')
+    test_plot_2d_decomposition('analytical', 'last')
+    plt.show()
diff --git a/pyproject.toml b/pyproject.toml
index c20cf49d3..a8df91d59 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -60,6 +60,7 @@ test = [
     "pytest-cov >= 5.0.0",
     'pytest >= 4.5',
     'pytest-xdist >= 1.16',
+    'Pillow',  # Python Imaging Library (PIL) fork
 ]
 
 [project.urls]
@@ -76,7 +77,7 @@ exclude = ["*__psydac__*"]
 namespaces = false
 
 [tool.setuptools.package-data]
-"*" = ["*.txt"]
+"*" = ["*.txt", "**/tests/data/*.png"]
 
 [tool.coverage.run]
 branch = true