diff --git a/CHANGELOG.md b/CHANGELOG.md
index d8fc460e8..48039f532 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,7 +5,8 @@ All notable changes to this project will be documented in this file.
 ## Unreleased
 
 ### Added
-
+-   #596 : Add `topetsc` method to `psydac.linalg.basic.ZeroOperator`
+-   #596 : Add `out` argument to `array_to_psydac` and `vec_topetsc`
 -   #576 : Add module `psydac.feec.polar.conga_projections` with broken FEEC polar projections in 2D
 -   #576 : Add 2D Poisson and TE Maxwell examples on polar mapped domains in `psydac/feec/polar/examples`
 -   #576 : Add module `psydac.utilities.parallel_utils` for parallel execution and gathering variable-length arrays
@@ -15,6 +16,8 @@ All notable changes to this project will be documented in this file.
 
 ### Fixed
 
+-   #596 : Allow export to VTK with a prescribed grid containing exact domain bounds
+-   #596 : Fix global geometric projectors with complex `dtype`
 -   #576 : Fix bug in `TensorFemSpace.eval_field` caused by round-off at MPI subdomain boundaries
 -   #576 : Require `sympde==0.19.3` which fixes a bug in the linearity checks
 -   #579 : Require `h5py>=3.16` which installs correctly with `setuptools>=81.0`
diff --git a/psydac/core/bsplines_kernels.py b/psydac/core/bsplines_kernels.py
index caaf75a95..105bf1a95 100644
--- a/psydac/core/bsplines_kernels.py
+++ b/psydac/core/bsplines_kernels.py
@@ -1095,7 +1095,7 @@ def cell_index_p(breaks: 'float[:]', i_grid: 'float[:]', tol: float, out: 'int[:
         i_cell = (low + high)//2
         it = 0
         max_it = 2*(high-low) # this number of iterations should not be reached
-        while it < max_it and (x < breaks[i_cell] - tol or x >= breaks[i_cell + 1] + tol):
+        while it < max_it and (x < breaks[i_cell] - tol or x > breaks[i_cell + 1] + tol):
             if x < breaks[i_cell]:
                 high = i_cell
             else:
diff --git a/psydac/feec/tests/test_global_projectors.py b/psydac/feec/tests/test_global_projectors.py
index f4a549e38..246179f36 100644
--- a/psydac/feec/tests/test_global_projectors.py
+++ b/psydac/feec/tests/test_global_projectors.py
@@ -25,8 +25,9 @@
 @pytest.mark.parametrize('degree', [1, 2, 3, 4, 5, 6, 7])
 @pytest.mark.parametrize('periodic', [False, True])
 @pytest.mark.parametrize('multiplicity', [1, 2])
+@pytest.mark.parametrize('dtype', [float, complex])
 
-def test_H1_projector_1d(domain, ncells, degree, periodic, multiplicity, verbose=False):
+def test_H1_projector_1d(domain, ncells, degree, periodic, multiplicity, dtype, verbose=False):
     if verbose:
         print('\n1d_h1')
 
@@ -39,13 +40,16 @@ def test_H1_projector_1d(domain, ncells, degree, periodic, multiplicity, verbose
 
     # H1 space (0-forms)
     N  = SplineSpace(degree=degree, knots=knots, periodic=periodic, basis='B')
-    V0 = TensorFemSpace(domain_decomposition, N)
+    V0 = TensorFemSpace(domain_decomposition, N, dtype=dtype)
 
     # Projector onto H1 space (1D interpolation)
     P0 = GlobalGeometricProjectorH1(V0)
 
     # Function to project
-    f  = lambda xi1 : np.sin( xi1 + 0.5 )
+    if dtype == complex:
+        f  = lambda xi1 : np.sin( xi1 + 0.5 -1j )
+    else:
+        f  = lambda xi1 : np.sin( xi1 + 0.5)
 
     # Compute the projection
     u0 = P0(f)
@@ -69,8 +73,9 @@ def test_H1_projector_1d(domain, ncells, degree, periodic, multiplicity, verbose
 @pytest.mark.parametrize('periodic', [True, False])
 @pytest.mark.parametrize('nquads', [100, 120, 140, 160])
 @pytest.mark.parametrize('multiplicity', [1, 2])
+@pytest.mark.parametrize('dtype', [float, complex])
 
-def test_L2_projector_1d(domain, ncells, degree, periodic, nquads, multiplicity, verbose=False):
+def test_L2_projector_1d(domain, ncells, degree, periodic, nquads, multiplicity, dtype, verbose=False):
     if verbose:
         print('\n1d_l2')
 
@@ -86,7 +91,7 @@ def test_L2_projector_1d(domain, ncells, degree, periodic, nquads, multiplicity,
     multiplicity = min(multiplicity, degree)
     # H1 space (0-forms)
     N  = SplineSpace(degree=degree, knots=knots, periodic=periodic, basis='B')
-    V0 = TensorFemSpace(domain_decomposition, N)
+    V0 = TensorFemSpace(domain_decomposition, N, dtype=dtype)
 
     # L2 space (1-forms)
     V1 = V0.reduce_degree(axes=[0], basis='M')
@@ -95,7 +100,10 @@ def test_L2_projector_1d(domain, ncells, degree, periodic, nquads, multiplicity,
     P1 = GlobalGeometricProjectorL2(V1, nquads=[nquads])
 
     # Function to project
-    f  = lambda xi1 : np.sin( xi1 + 0.5 )
+    if dtype == complex:
+        f  = lambda xi1 : np.sin( xi1 + 0.5j )
+    else:
+        f  = lambda xi1 : np.sin( xi1 + 0.5 )
 
     # Compute the projection
     u1 = P1(f)
@@ -117,8 +125,9 @@ def test_L2_projector_1d(domain, ncells, degree, periodic, nquads, multiplicity,
 @pytest.mark.parametrize('degree', [[2,2], [2,3], [3,3]])
 @pytest.mark.parametrize('periodic', [[False, False], [True, True]])
 @pytest.mark.parametrize('multiplicity', [(1, 1), (2, 2)])
+@pytest.mark.parametrize('dtype', [float, complex])
 
-def test_derham_projector_2d_hdiv(ncells, degree, periodic, multiplicity, verbose=False):
+def test_derham_projector_2d_hdiv(ncells, degree, periodic, multiplicity, dtype, verbose=False):
     if verbose:
         print('\n2d_hdiv')
 
@@ -129,12 +138,20 @@ def test_derham_projector_2d_hdiv(ncells, degree, periodic, multiplicity, verbos
     multiplicity = [min(m, p) for p, m in zip(degree, multiplicity)]
 
     derham   = Derham(domain, ["H1", "Hdiv", "L2"])
+    if dtype == complex:
+        derham.V0.codomain_type='complex'
+        derham.V1.codomain_type='complex' 
+        derham.V2.codomain_type='complex'
+
     derham_h   = discretize(derham, domain_h, degree=degree, get_H1vec_space = True, multiplicity=multiplicity)
     P0, P1, P2, PX = derham_h.projectors(nquads=[2*p+1 for p in degree])
 
     # Function to project
     f1  = lambda xi1, xi2 : np.sin( xi1 + 0.5 ) * np.cos( xi2 + 0.3 )
-    f2  = lambda xi1, xi2 : np.cos( xi1 + 0.5 ) * np.sin( xi2 - 0.2 )
+    if dtype == complex:
+        f2  = lambda xi1, xi2 : np.cos( xi1 + 0.5 ) * np.sin( xi2 - 0.2j )
+    else:
+        f2  = lambda xi1, xi2 : np.cos( xi1 + 0.5 ) * np.sin( xi2 - 0.2 )
 
     # Compute the projection
     u0 = P0(f1)
@@ -181,8 +198,9 @@ def test_derham_projector_2d_hdiv(ncells, degree, periodic, multiplicity, verbos
 @pytest.mark.parametrize('degree', [[2,2], [2,3], [3,3]])
 @pytest.mark.parametrize('periodic', [[False, False], [True, True]])
 @pytest.mark.parametrize('multiplicity', [(1, 1), (2, 2)])
+@pytest.mark.parametrize('dtype', [float, complex])
 
-def test_derham_projector_2d_hdiv_2(ncells, degree, periodic, multiplicity, verbose=False):
+def test_derham_projector_2d_hdiv_2(ncells, degree, periodic, multiplicity, dtype, verbose=False):
     if verbose:
         print('\n2d_hdiv_2')
 
@@ -193,11 +211,18 @@ def test_derham_projector_2d_hdiv_2(ncells, degree, periodic, multiplicity, verb
     multiplicity = [min(m, p) for p, m in zip(degree, multiplicity)]
     
     derham   = Derham(domain, ["H1", "Hdiv", "L2"])
+    if dtype == complex:
+        derham.V0.codomain_type='complex'
+        derham.V1.codomain_type='complex' 
+        derham.V2.codomain_type='complex'
     derham_h   = discretize(derham, domain_h, degree=degree, get_H1vec_space = True, multiplicity=multiplicity)
     P0, P1, P2, PX = derham_h.projectors()
 
     # Function to project
-    f1  = lambda xi1, xi2 : 20 * xi1**2*(xi1-1.)**2 
+    if dtype == complex:
+        f1  = lambda xi1, xi2 : (20+2j) * xi1**2*(xi1-1.)**2
+    else:
+        f1  = lambda xi1, xi2 : 20 * xi1**2*(xi1-1.)**2 
     f2  = lambda xi1, xi2 : 10 * xi2**2*(xi2-1.)**2
 
     # Compute the projection
@@ -253,8 +278,9 @@ def test_derham_projector_2d_hdiv_2(ncells, degree, periodic, multiplicity, verb
 @pytest.mark.parametrize('degree', [[2,2], [2,3], [3,3]])
 @pytest.mark.parametrize('periodic', [[False, False], [True, False] ,[True, True]])
 @pytest.mark.parametrize('multiplicity', [[1,1],[2,2]])
+@pytest.mark.parametrize('dtype', [float, complex])
 
-def test_derham_projector_2d_hcurl(ncells, degree, periodic, multiplicity, verbose=False):
+def test_derham_projector_2d_hcurl(ncells, degree, periodic, multiplicity, dtype, verbose=False):
     if verbose:
         print('\n2d_hcurl')
 
@@ -265,12 +291,20 @@ def test_derham_projector_2d_hcurl(ncells, degree, periodic, multiplicity, verbo
     multiplicity = [min(m, p) for p, m in zip (degree, multiplicity)]
     
     derham   = Derham(domain, ["H1", "Hcurl", "L2"])
+    if dtype == complex:
+        derham.V0.codomain_type='complex'
+        derham.V1.codomain_type='complex' 
+        derham.V2.codomain_type='complex'
     derham_h   = discretize(derham, domain_h, degree=degree, get_H1vec_space = True, multiplicity=multiplicity)
     P0, P1, P2, PX = derham_h.projectors()
 
     # Function to project
-    f1  = lambda xi1, xi2 : np.sin( xi1 + 0.5 ) * np.cos( xi2 + 0.3 )
-    f2  = lambda xi1, xi2 : np.cos( xi1 + 0.5 ) * np.sin( xi2 - 0.2 )
+    if dtype == complex:
+        f1  = lambda xi1, xi2 : np.sin( xi1 + 0.5 ) * np.cos( xi2 + 0.3j )
+        f2  = lambda xi1, xi2 : 1j * np.cos( xi1 + 0.5 ) * np.sin( xi2 - 0.2 )
+    else:
+        f1  = lambda xi1, xi2 : np.sin( xi1 + 0.5 ) * np.cos( xi2 + 0.3 )
+        f2  = lambda xi1, xi2 : np.cos( xi1 + 0.5 ) * np.sin( xi2 - 0.2 )
 
     # Compute the projection
     u0 = P0(f1)
@@ -311,7 +345,7 @@ def test_derham_projector_2d_hcurl(ncells, degree, periodic, multiplicity, verbo
     if verbose:
         print(ncells, maxnorm_error / error_estim)
     assert maxnorm_error <= max(10 * error_estim, 1e-13)
-    
+
 #==============================================================================
 @pytest.mark.parametrize('ncells', [[10,9,12]])
 @pytest.mark.parametrize('degree', [[2,2,2], [2,3,2], [3,3,3]])
@@ -324,6 +358,10 @@ def test_derham_projector_3d(ncells, degree, periodic, multiplicity, verbose=Fal
     domain_h = discretize(domain, ncells=ncells, periodic=periodic)
     
     derham   = Derham(domain)
+    derham.V0.codomain_type='complex'
+    derham.V1.codomain_type='complex' 
+    derham.V2.codomain_type='complex'
+    derham.V3.codomain_type='complex'
     #change multiplicity if higher than degree to avoid problems (case p<m doesn't work)
     multiplicity = [min(m, p) for p, m in zip(degree, multiplicity)]
 
@@ -333,7 +371,7 @@ def test_derham_projector_3d(ncells, degree, periodic, multiplicity, verbose=Fal
     # Function to project
     f1 = lambda xi1, xi2, xi3 : np.sin( xi1 + 0.51 ) * np.cos( xi2 + 0.32 ) * np.sin( xi3 - 0.43)
     f2 = lambda xi1, xi2, xi3 : np.cos( xi1 + 0.27 ) * np.sin( xi2 - 0.29 ) * np.cos( xi3 + 0.67)
-    f3 = lambda xi1, xi2, xi3 : np.cos( xi1 + 0.72 ) * np.sin( xi2 - 0.73 ) * np.cos( xi3 - 0.14)
+    f3 = lambda xi1, xi2, xi3 : 1j * np.cos( xi1 + 0.72 ) * np.sin( xi2 - 0.73 ) * np.cos( xi3 - 0.14)
 
     # Compute the projection
     u0 = P0(f1)
@@ -407,23 +445,25 @@ def manual_convergence_tests(dim):
     if dim == 1:
         degrees = [3,4,5,6,7]
         ncells = [10, 20, 40, 80, 160, 320, 640]
+        dtype = float
         for deg in degrees:
             print(f"degree = {deg}")
             for nc in ncells:
-                test_H1_projector_1d(domain, nc, deg, periodic, multiplicity=1, verbose=True)
-                test_H1_projector_1d(domain, nc, deg, periodic, multiplicity=2, verbose=True)
-                test_L2_projector_1d(domain, nc, deg, periodic, nquads=deg, multiplicity=1, verbose=True)
+                test_H1_projector_1d(domain, nc, deg, periodic, multiplicity=1, dtype=dtype, verbose=True)
+                test_H1_projector_1d(domain, nc, deg, periodic, multiplicity=2, dtype=dtype, verbose=True)
+                test_L2_projector_1d(domain, nc, deg, periodic, nquads=deg, multiplicity=1, dtype=dtype, verbose=True)
         
     elif dim == 2:
         degrees = [[3,2], [3, 3], [3,4]]
         ncells = [10, 20, 40, 80, 160]
+        dtype = float
         for deg in degrees:
             print(f"degree = {deg}")
             for nc in ncells:
-                test_derham_projector_2d_hcurl ([nc, nc], deg, [periodic, periodic], multiplicity=[1,1], verbose=True)
-                test_derham_projector_2d_hdiv  ([nc, nc], deg, [periodic, periodic], multiplicity=[1,1], verbose=True)
-                test_derham_projector_2d_hdiv  ([nc, nc], deg, [periodic, periodic], multiplicity=[2,2], verbose=True)
-                test_derham_projector_2d_hdiv_2([nc, nc], deg, [periodic, periodic], multiplicity=[1,1], verbose=True)
+                test_derham_projector_2d_hcurl ([nc, nc], deg, [periodic, periodic], multiplicity=[1,1], dtype=dtype, verbose=True)
+                test_derham_projector_2d_hdiv  ([nc, nc], deg, [periodic, periodic], multiplicity=[1,1], dtype=dtype, verbose=True)
+                test_derham_projector_2d_hdiv  ([nc, nc], deg, [periodic, periodic], multiplicity=[2,2], dtype=dtype, verbose=True)
+                test_derham_projector_2d_hdiv_2([nc, nc], deg, [periodic, periodic], multiplicity=[1,1], dtype=dtype, verbose=True)
     
     elif dim == 3:
         print("\n3d")
diff --git a/psydac/fem/splines.py b/psydac/fem/splines.py
index f36f54015..7c6b9e394 100644
--- a/psydac/fem/splines.py
+++ b/psydac/fem/splines.py
@@ -185,7 +185,7 @@ def init_interpolation( self, dtype=float ):
 
         if self.periodic:
             # Convert to CSC format and compute sparse LU decomposition
-            self._interpolator = SparseSolver( csc_matrix( imat ) )
+            self._interpolator = SparseSolver( csc_matrix(imat, dtype=dtype) )
         else:
             # Convert to LAPACK banded format (see DGBTRF function)
             dmat = dia_matrix( imat )
@@ -221,7 +221,7 @@ def init_histopolation( self, dtype=float):
         self.hmat= imat
         if self.periodic:
             # Convert to CSC format and compute sparse LU decomposition
-            self._histopolator = SparseSolver( csc_matrix( imat ) )
+            self._histopolator = SparseSolver( csc_matrix(imat, dtype=dtype) )
         else:
             # Convert to LAPACK banded format (see DGBTRF function)
             dmat = dia_matrix( imat )
diff --git a/psydac/linalg/basic.py b/psydac/linalg/basic.py
index a6bc99b66..402afc67f 100644
--- a/psydac/linalg/basic.py
+++ b/psydac/linalg/basic.py
@@ -553,12 +553,16 @@ def copy(self):
         return ZeroOperator(self.domain, self.codomain)
 
     def toarray(self):
-        return np.zeros(self.shape, dtype=self.dtype) 
+        return np.zeros(self.shape, dtype=self.dtype)
 
     def tosparse(self):
         from scipy.sparse import csr_matrix
         return csr_matrix(self.shape, dtype=self.dtype)
 
+    def topetsc(self):
+        from psydac.linalg.topetsc import mat_topetsc
+        return mat_topetsc(self)
+
     def transpose(self, conjugate=False):
         return ZeroOperator(domain=self.codomain, codomain=self.domain)
 
diff --git a/psydac/linalg/block.py b/psydac/linalg/block.py
index 846535d75..dc3f73388 100644
--- a/psydac/linalg/block.py
+++ b/psydac/linalg/block.py
@@ -482,11 +482,11 @@ def toarray_local(self, order='C'):
         return np.block([blocks])[0]
 
     # ...
-    def topetsc(self):
+    def topetsc(self, out=None):
         """ Convert to petsc data structure.
         """
         from psydac.linalg.topetsc import vec_topetsc
-        vec = vec_topetsc( self )
+        vec = vec_topetsc( self, out=out )
         return vec
 
 #===============================================================================
@@ -886,17 +886,20 @@ def nonzero_block_indices(self):
     # ...
     def update_ghost_regions(self):
         for Lij in self._blocks.values():
-            Lij.update_ghost_regions()
+            if hasattr(Lij, 'update_ghost_regions'):
+                Lij.update_ghost_regions()
 
     # ...
     def exchange_assembly_data(self):
         for Lij in self._blocks.values():
-            Lij.exchange_assembly_data()
+            if hasattr(Lij, 'exchange_assembly_data'):
+                Lij.exchange_assembly_data()
 
     # ...
-    def remove_spurious_entries(self ):
+    def remove_spurious_entries(self):
         for Lij in self._blocks.values():
-            Lij.remove_spurious_entries()
+            if hasattr(Lij, 'remove_spurious_entries'):
+                Lij.remove_spurious_entries()
 
     @property
     def ghost_regions_in_sync(self):
diff --git a/psydac/linalg/stencil.py b/psydac/linalg/stencil.py
index fe3288be0..a95bdd927 100644
--- a/psydac/linalg/stencil.py
+++ b/psydac/linalg/stencil.py
@@ -698,11 +698,11 @@ def _toarray_parallel_with_pads(self, order='C'):
         return out.flatten( order=order)
 
     #...
-    def topetsc(self):
+    def topetsc(self, out=None):
         """ Convert to petsc data structure.
         """
         from psydac.linalg.topetsc import vec_topetsc
-        vec = vec_topetsc( self )
+        vec = vec_topetsc( self, out=out )
         return vec
 
     # ...
diff --git a/psydac/linalg/tests/test_block.py b/psydac/linalg/tests/test_block.py
index 69159297b..9760396b6 100644
--- a/psydac/linalg/tests/test_block.py
+++ b/psydac/linalg/tests/test_block.py
@@ -904,6 +904,8 @@ def test_block_vector_2d_serial_topetsc( dtype, n1, n2, p1, p2, P1, P2 ):
     W = BlockVectorSpace(V1, V2)
 
     x = BlockVector(W)
+    y = BlockVector(W)
+    z = BlockVector(W)
 
     # Initialize random number generator (set seed for reproducibility)
     rng = np.random.default_rng(seed=n1*n2*p1*p2)
@@ -917,10 +919,16 @@ def test_block_vector_2d_serial_topetsc( dtype, n1, n2, p1, p2, P1, P2 ):
     x.update_ghost_regions()
 
     v = x.topetsc()
-    v = petsc_to_psydac(v, W)
+    u = v.copy()
+    x.topetsc(out=u)
 
-    # The vectors can only be compared in the serial case
-    assert np.allclose( x.toarray() , v.toarray() )
+    # Convert PETSc.Vec to BlockVector of W
+    y = petsc_to_psydac(v, W)
+    petsc_to_psydac(u, W, out=z)
+
+    assert np.array_equal(x.toarray(), y.toarray())
+    assert np.array_equal(x.toarray(), z.toarray())
+    assert np.array_equal(v.array, u.array) 
 
 #===============================================================================
 @pytest.mark.parametrize( 'dtype', [float, complex] )
@@ -1343,6 +1351,8 @@ def test_block_vector_2d_parallel_topetsc( dtype, n1, n2, p1, p2, P1, P2 ):
     #TODO: implement conversion to PETSc recursively to treat case of block of blocks 
 
     x = BlockVector(W)
+    y = BlockVector(W)
+    z = BlockVector(W)
 
     # Initialize random number generator (set seed for reproducibility)
     rng = np.random.default_rng(seed=n1*n2*p1*p2)
@@ -1355,9 +1365,18 @@ def test_block_vector_2d_parallel_topetsc( dtype, n1, n2, p1, p2, P1, P2 ):
 
     x.update_ghost_regions()
 
-    v = petsc_to_psydac(x.topetsc(), W)
+    v = x.topetsc()
+    u = v.copy()
+    x.topetsc(out=u)
+
+    # Convert PETSc.Vec to BlockVector of W
+    y = petsc_to_psydac(v, W)
+    petsc_to_psydac(u, W, out=z)
+
+    assert np.array_equal(x.toarray(), y.toarray())
+    assert np.array_equal(x.toarray(), z.toarray())
+    assert np.array_equal(v.array, u.array)
 
-    assert np.allclose( x.toarray() , v.toarray(), rtol=1e-12, atol=1e-12 )
 
 #=============================================================================== 
 @pytest.mark.parametrize( 'dtype', [float, complex] )
diff --git a/psydac/linalg/tests/test_matrix_free.py b/psydac/linalg/tests/test_matrix_free.py
index 1bf1a9fa7..49a657904 100644
--- a/psydac/linalg/tests/test_matrix_free.py
+++ b/psydac/linalg/tests/test_matrix_free.py
@@ -128,6 +128,31 @@ def test_solvers_matrix_free(solver, use_jacobi_pc):
     error = np.linalg.norm( (b - y).toarray())
     assert np.linalg.norm( (b - y).toarray() ) < tol
 
+@pytest.mark.parametrize('n1', [3, 5])
+@pytest.mark.parametrize('n2', [4, 7])
+@pytest.mark.parametrize('p1', [2, 6])
+@pytest.mark.parametrize('p2', [3, 9])
+@pytest.mark.petsc
+def test_ZeroOperator_topetsc(n1, n2, p1, p2):
+    P1 = False
+    P2 = False
+    m1 = (n2+n1)//2
+    m2 = n1+1
+    q1 = p1
+    q2 = p2 
+    V1 = get_StencilVectorSpace([n1, n2], [p1, p2], [P1, P2])
+    V2 = get_StencilVectorSpace([m1, m2], [q1, q2], [P1, P2])
+    O = ZeroOperator(domain=V1, codomain=V2)
+    Op = O.topetsc()
+
+    v1 = get_random_StencilVector(V1)
+    v1p = v1.topetsc()
+    v2p = Op.createVecLeft()
+
+    Op.mult(v1p, v2p)
+
+    assert np.linalg.norm( v2p.array ) == 0.
+
 #===============================================================================
 # SCRIPT FUNCTIONALITY
 #===============================================================================
diff --git a/psydac/linalg/tests/test_stencil_vector.py b/psydac/linalg/tests/test_stencil_vector.py
index 0d0f0c7ca..432f21eee 100644
--- a/psydac/linalg/tests/test_stencil_vector.py
+++ b/psydac/linalg/tests/test_stencil_vector.py
@@ -432,6 +432,8 @@ def test_stencil_vector_2d_serial_topetsc(dtype, n1, n2, p1, p2, s1, s2, P1, P2)
     # Create vector space and stencil vector
     V = StencilVectorSpace(C, dtype=dtype)
     x = StencilVector(V)
+    y = StencilVector(V)
+    z = StencilVector(V)
 
     # Fill the vector with data
     if dtype == complex:
@@ -444,19 +446,24 @@ def test_stencil_vector_2d_serial_topetsc(dtype, n1, n2, p1, p2, s1, s2, P1, P2)
 
     # Convert vector to PETSc.Vec
     v = x.topetsc()
+    w = v.copy()
+    x.topetsc(out=w)
 
     # Convert PETSc.Vec to StencilVector of V
-    v = petsc_to_psydac(v, V)
+    z = petsc_to_psydac(v, V)
+    petsc_to_psydac(w, V, out=y)
 
     # Test properties of v and data contained
-    assert v.space is V
-    assert v.dtype == dtype
-    assert v.starts == (0, 0)
-    assert v.ends == (n1 - 1, n2 - 1)
-    assert v.pads == (p1, p2)
-    assert v._data.shape == (n1 + 2 * p1 * s1, n2 + 2 * p2 * s2)
-    assert v._data.dtype == dtype
-    assert np.array_equal(x.toarray(), v.toarray())
+    assert z.space is V
+    assert z.dtype == dtype
+    assert z.starts == (0, 0)
+    assert z.ends == (n1 - 1, n2 - 1)
+    assert z.pads == (p1, p2)
+    assert z._data.shape == (n1 + 2 * p1 * s1, n2 + 2 * p2 * s2)
+    assert z._data.dtype == dtype
+    assert np.array_equal(x.toarray(), z.toarray())
+    assert np.array_equal(x.toarray(), y.toarray())
+    assert np.array_equal(v.array, w.array) 
 
 # ===============================================================================
 @pytest.mark.parametrize('dtype', [float, complex])
@@ -619,6 +626,8 @@ def test_stencil_vector_2d_parallel_topetsc(dtype, n1, n2, p1, p2, s1, s2, P1, P
     # Create vector space and stencil vector
     V = StencilVectorSpace(C, dtype=dtype)
     x = StencilVector(V)
+    y = StencilVector(V)
+    z = StencilVector(V)
 
     # Fill the vector with data
     if dtype == complex:
@@ -633,11 +642,16 @@ def test_stencil_vector_2d_parallel_topetsc(dtype, n1, n2, p1, p2, s1, s2, P1, P
 
     # Convert vector to PETSc.Vec
     v = x.topetsc()
+    u = v.copy()
+    x.topetsc(out=u)
 
     # Convert PETSc.Vec to StencilVector of V
-    v = petsc_to_psydac(v, V)
+    y = petsc_to_psydac(v, V)
+    petsc_to_psydac(u, V, out=z)
 
-    assert np.array_equal(x.toarray(), v.toarray())
+    assert np.array_equal(x.toarray(), y.toarray())
+    assert np.array_equal(x.toarray(), z.toarray())
+    assert np.array_equal(v.array, u.array)    
     
 # ===============================================================================
 @pytest.mark.parametrize('dtype', [float, complex])
@@ -662,6 +676,8 @@ def test_stencil_vector_1d_parallel_topetsc(dtype, n1, p1, s1, P1):
     # Create vector space and stencil vector
     V = StencilVectorSpace(C, dtype=dtype)
     x = StencilVector(V)
+    y = StencilVector(V)
+    z = StencilVector(V)
 
     # Fill the vector with data
     if dtype == complex:
@@ -675,11 +691,16 @@ def test_stencil_vector_1d_parallel_topetsc(dtype, n1, p1, s1, P1):
 
     # Convert vector to PETSc.Vec
     v = x.topetsc()
+    u = v.copy()
+    x.topetsc(out=u)    
 
     # Convert PETSc.Vec to StencilVector of V
-    v = petsc_to_psydac(v, V)
+    y = petsc_to_psydac(v, V)
+    petsc_to_psydac(u, V, out=z)
 
-    assert np.array_equal(x.toarray(), v.toarray())
+    assert np.array_equal(x.toarray(), y.toarray())
+    assert np.array_equal(x.toarray(), z.toarray())
+    assert np.array_equal(v.array, u.array) 
 
 # ===============================================================================
 @pytest.mark.parametrize('dtype', [float, complex])
@@ -713,6 +734,8 @@ def test_stencil_vector_3d_parallel_topetsc(dtype, n1, n2, n3, p1, p2, p3, s1, s
     # Create vector space and stencil vector
     V = StencilVectorSpace(C, dtype=dtype)
     x = StencilVector(V)
+    y = StencilVector(V)
+    z = StencilVector(V)
 
     # Fill the vector with data
 
@@ -729,11 +752,16 @@ def test_stencil_vector_3d_parallel_topetsc(dtype, n1, n2, n3, p1, p2, p3, s1, s
 
     # Convert vector to PETSc.Vec
     v = x.topetsc()
+    u = v.copy()
+    x.topetsc(out=u)
 
     # Convert PETSc.Vec to StencilVector of V
-    v = petsc_to_psydac(v, V)
+    y = petsc_to_psydac(v, V)
+    petsc_to_psydac(u, V, out=z)
 
-    assert np.array_equal(x.toarray(), v.toarray())
+    assert np.array_equal(x.toarray(), y.toarray())
+    assert np.array_equal(x.toarray(), z.toarray())
+    assert np.array_equal(v.array, u.array) 
 
 # ===============================================================================
 @pytest.mark.parametrize('dtype', [float, complex])
diff --git a/psydac/linalg/topetsc.py b/psydac/linalg/topetsc.py
index c6365a21a..44a9afa59 100644
--- a/psydac/linalg/topetsc.py
+++ b/psydac/linalg/topetsc.py
@@ -7,7 +7,7 @@
 
 import numpy as np
 
-from psydac.linalg.basic   import VectorSpace
+from psydac.linalg.basic   import VectorSpace, ZeroOperator
 from psydac.linalg.block   import BlockVectorSpace, BlockVector, BlockLinearOperator
 from psydac.linalg.stencil import StencilVectorSpace, StencilVector, StencilMatrix
 from psydac.linalg.kernels.stencil2IJV_kernels import stencil2IJV_1d_C, stencil2IJV_2d_C, stencil2IJV_3d_C
@@ -336,14 +336,16 @@ def get_npts_per_block(V : VectorSpace) -> list:
     return npts_local_per_block
 
 
-def vec_topetsc(vec):
+def vec_topetsc(vec, out=None):
     """ Convert vector from PSYDAC format to a PETSc.Vec object.
 
     Parameters
     ----------
     vec : psydac.linalg.stencil.StencilVector | psydac.linalg.block.BlockVector
       PSYDAC StencilVector or BlockVector. In the case of a BlockVector, only the case where the blocks are StencilVector is implemented.
-
+    out : PETSc.Vec, optional
+      PETSc vector where to store the returned value. In case of `None`, a new vector is created.
+     
     Returns
     -------
     gvec : PETSc.Vec
@@ -374,13 +376,19 @@ def vec_topetsc(vec):
     # Sum over the blocks to get the total local size
     localsize = np.sum(np.prod(npts_local, axis=1))
 
-    gvec  = PETSc.Vec().create(comm=carts[0].global_comm)    
+    if out is not None:
+        assert isinstance(out, PETSc.Vec)
+        assert out.getSizes() == (localsize, globalsize)
+        out.zeroEntries()
+        gvec = out
+    else:
+        gvec  = PETSc.Vec().create(comm=carts[0].global_comm)    
 
-    # Set global and local size:
-    gvec.setSizes(size=(localsize, globalsize))
+        # Set global and local size:
+        gvec.setSizes(size=(localsize, globalsize))
 
-    gvec.setFromOptions()
-    gvec.setUp()
+        gvec.setFromOptions()
+        gvec.setUp()
 
     petsc_indices = []
     petsc_data = []
@@ -452,7 +460,7 @@ def mat_topetsc(mat):
 
     from petsc4py import PETSc
 
-    assert isinstance(mat, StencilMatrix) or isinstance(mat, BlockLinearOperator), 'Conversion only implemented for StencilMatrix and BlockLinearOperator.'
+    assert isinstance(mat, StencilMatrix) or isinstance(mat, BlockLinearOperator) or isinstance(mat, ZeroOperator), 'Conversion only implemented for StencilMatrix, BlockLinearOperator and ZeroOperator.'
 
     if (isinstance(mat.domain, BlockVectorSpace) and any([isinstance(mat.domain.spaces[b], BlockVectorSpace) for b in range(len(mat.domain.spaces))]))\
         or (isinstance(mat.codomain, BlockVectorSpace) and any([isinstance(mat.codomain.spaces[b], BlockVectorSpace) for b in range(len(mat.codomain.spaces))])):
@@ -463,23 +471,11 @@ def mat_topetsc(mat):
     elif isinstance(mat.domain, BlockVectorSpace):
         comm = mat.domain.spaces[0].cart.global_comm
 
-    nonzero_block_indices = ((0,0),) if isinstance(mat, StencilMatrix) else mat.nonzero_block_indices
-
-    mat.update_ghost_regions()
-    mat.remove_spurious_entries()
 
     # Get the number of points local to the current process:
     dnpts_local = get_npts_local(mat.domain) # indexed [block, dimension]. Different for each process.
     cnpts_local = get_npts_local(mat.codomain) # indexed [block, dimension]. Different for each process. 
 
-    # Get the number of points per block, per process and per dimension:
-    dnpts_per_block_per_process = np.array(get_npts_per_block(mat.domain)) # global variable, indexed as [block, process, dimension]
-    cnpts_per_block_per_process = np.array(get_npts_per_block(mat.codomain)) # global variable, indexed as [block, process, dimension]
-
-    # Get the index shift for each block and each process:
-    dindex_shift = get_index_shift_per_block_per_process(mat.domain) # global variable, indexed as [block, process, dimension]
-    cindex_shift = get_index_shift_per_block_per_process(mat.codomain) # global variable, indexed as [block, process, dimension]
-
     globalsize = mat.shape
 
     # Sum over the blocks to get the total local size
@@ -500,6 +496,23 @@ def mat_topetsc(mat):
     gmat.setFromOptions()
     gmat.setUp()
 
+    if isinstance(mat, ZeroOperator):
+        gmat.assemble()
+        return gmat
+
+    nonzero_block_indices = ((0,0),) if isinstance(mat, StencilMatrix) else mat.nonzero_block_indices
+
+    mat.update_ghost_regions()
+    mat.remove_spurious_entries()
+
+    # Get the number of points per block, per process and per dimension (global variable) indexed as [block, process, dimension]
+    dnpts_per_block_per_process = np.array(get_npts_per_block(mat.domain))
+    cnpts_per_block_per_process = np.array(get_npts_per_block(mat.codomain))
+
+    # Get the index shift per block, per process and per dimension (global variable) indexed as [block, process, dimension]
+    dindex_shift = get_index_shift_per_block_per_process(mat.domain)
+    cindex_shift = get_index_shift_per_block_per_process(mat.codomain)
+
     I = [0] # Row pointers
     J = [] # Column indices
     V = [] # Values
@@ -510,6 +523,10 @@ def mat_topetsc(mat):
     for bc, bd in nonzero_block_indices:
         if isinstance(mat, BlockLinearOperator):
             mat_block = mat.blocks[bc][bd]
+
+        if isinstance(mat_block, ZeroOperator):
+            continue
+
         dnpts_block = dnpts_per_block_per_process[bd]
         cnpts_block = cnpts_per_block_per_process[bc]
         dshift_block = dindex_shift[bd]
diff --git a/psydac/linalg/utilities.py b/psydac/linalg/utilities.py
index db86e075d..4b8b84640 100644
--- a/psydac/linalg/utilities.py
+++ b/psydac/linalg/utilities.py
@@ -19,7 +19,7 @@
 )
 
 #==============================================================================
-def array_to_psydac(x, V):
+def array_to_psydac(x, V, out=None):
     """ 
     Convert a NumPy array to a Vector of the space V. This function is designed
     to be the inverse of the method .toarray() of the class Vector.
@@ -32,6 +32,9 @@ def array_to_psydac(x, V):
     V : psydac.linalg.stencil.StencilVectorSpace or psydac.linalg.block.BlockVectorSpace
         Space of the final PSYDAC Vector.
 
+    out : psydac.linalg.stencil.StencilVector | psydac.linalg.block.BlockVector, optional
+      The Psydac vector where to store the result.
+
     Returns
     -------
     u : psydac.linalg.stencil.StencilVector or psydac.linalg.block.BlockVector
@@ -47,7 +50,13 @@ def array_to_psydac(x, V):
         assert V.dtype==complex, 'Complex array cannot be converted to a real StencilVector'
     assert x.size == V.dimension, 'Array must have the same global size as the space.'
 
-    u = V.zeros()
+    if out is None:
+        u = V.zeros()
+    else:
+        assert isinstance(out, BlockVector)
+        assert out.space is V
+        u = out
+
     _array_to_psydac_recursive(x, u)
     u.update_ghost_regions()