"Non-square" MPC giving different sparsity pattern than square problem.

[jorgensd]

from mpi4py import MPI

import numpy as np
from dolfinx import fem, mesh
from basix.ufl import element
import dolfinx_mpc
import ufl

cell_type = mesh.CellType.hexahedron # BUG

N = 3
domain = mesh.create_unit_cube(MPI.COMM_WORLD, N, 2, N, cell_type=cell_type)

V = fem.functionspace(domain, element("Lagrange", domain.basix_cell(), 1))
Q = V.clone()

def periodic_boundary(x):
    return np.logical_or(np.isclose(x[0], 1), np.isclose(x[2], 1))

def periodic_map(x):
    out = x.copy()
    out[0][np.isclose(x[0], 1).nonzero()] -= 1
    out[2][np.isclose(x[2], 1).nonzero()] -= 1
    return out

mpc_u = dolfinx_mpc.MultiPointConstraint(V)
mpc_u.create_periodic_constraint_geometrical(V, periodic_boundary, periodic_map, [])
mpc_u.finalize()

mpc_p = dolfinx_mpc.MultiPointConstraint(Q)
mpc_p.create_periodic_constraint_geometrical(Q, periodic_boundary, periodic_map, [])
mpc_p.finalize()

# Stokes weak form
u, p = ufl.TrialFunction(V), ufl.TrialFunction(Q)
v, q = ufl.TestFunction(V), ufl.TestFunction(Q)

# Weak statementmpc_u,
a01 = v * p * ufl.dx
form_c = fem.form(a01)
mpcs = [mpc_u, mpc_p]

p0 =  dolfinx_mpc.create_sparsity_pattern(form_c, mpcs)
p0.finalize()

p1 =  dolfinx_mpc.create_sparsity_pattern(form_c, [mpc_u, mpc_u])
p1.finalize()
print(p0.num_nonzeros, p1.num_nonzeros)
breakpoint()

This in turn probably leads to an out of range error in assembly.