# Incorrect PETSc.Mat block size when using enriched\_element

**URL:** <https://fenicsproject.discourse.group/t/incorrect-petsc-mat-block-size-when-using-enriched-element/19805>\
**Category:** dolfinx\
**Created:** [August 21, 2026, 9:35am UTC](https://fenicsproject.discourse.group/t/incorrect-petsc-mat-block-size-when-using-enriched-element/19805 "2026-08-21T09:35:18Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![TimS](https://avatars.discourse-cdn.com/v4/letter/t/7ea924/32.png) [@TimS](https://fenicsproject.discourse.group/u/TimS)\
**Post date:** [August 21, 2026, 9:35am UTC](https://fenicsproject.discourse.group/t/incorrect-petsc-mat-block-size-when-using-enriched-element/19805/1 "2026-08-21T09:35:18Z")

</div>

General info (even though I do not think it matters much):

- MacOS 26.6.1
- Python 3.14.6 (conda)
- fenics-dolfinx 0.11.0 (conda)
- petsc4py 3.25.4 (conda)

In my application, I am using the MINI element (P1 Lagrange enriched with bubble) for the (Navier-)Stokes equations. As I would like to use gamg or hypre as a multigrid preconditioner (which, at the moment, is working mediocre at best), I encountered a comment about gamg in the PETSc documentation ([PCGAMG — PETSc 3.25.4 documentation](https://petsc.org/release/manualpages/PC/PCGAMG/)) that describes the requirement “to indicate the number of degrees of freedom per grid point”. I therefore investigated and found that the submatrix (0,0) from the nested PETSc.Mat matrix only had a block\_size of 1. I have broken my problem down further (the nested approach does not actually seem to matter) to pinpoint the issue (adapted from [Stokes equations using Taylor-Hood elements — DOLFINx Python 0.12.0.dev0](https://docs.fenicsproject.org/dolfinx/main/python/demos/demo_stokes.html#high-level-nested-matrix-solver) ). A MWE can be found in the following (I know it is not a fully working solver (e.g., boundary conditions are missing), but I only included the absolutely necessary parts to reproduce the issue).

Am I doing something wrong, misunderstanding something, is this information actually not required, or is there actually an issue?

```auto
from mpi4py import MPI
from petsc4py import PETSc

import numpy as np

from dolfinx.mesh import CellType, create_rectangle
from basix.ufl import element, enriched_element
from dolfinx import default_real_type
from dolfinx.fem import functionspace, Constant
from ufl import TestFunction, TrialFunction, inner, dot, grad, dx
from dolfinx.fem.petsc import LinearProblem

# Create mesh
msh = create_rectangle(
    MPI.COMM_WORLD, [np.array([0, 0]), np.array([1, 1])], (32, 32), CellType.triangle
)
gdim = msh.geometry.dim

# some definitions for the linear problem later
petsc_options_prefix = "block_size_test_"
petsc_options = {
    "ksp_type": "gmres",
    "ksp_rtol": 1e-9,
    "pc_type": "gamg",
}

# source term, zero vector for both approaches
f = Constant(msh, [PETSc.ScalarType(0.)] * gdim)

def use_P2():
    # P2 element
    P2 = element("Lagrange", msh.basix_cell(), degree=2, shape=(gdim,), dtype=default_real_type)
    V_P2 = functionspace(msh, P2)

    # velocity trial and test function
    u_P2 = TrialFunction(V_P2)
    v_P2 = TestFunction(V_P2)

    # variational problem
    a_ufl_P2 = inner(grad(u_P2), grad(v_P2)) * dx
    L_ufl_P2 = dot(f, v_P2) * dx

    # define solver
    problem_P2 = LinearProblem(
        a_ufl_P2,
        L_ufl_P2,
        # P=a_p,
        petsc_options_prefix = petsc_options_prefix + "P2_",
        petsc_options = petsc_options,
    )

    # print block size with P2 element
    print("Block size with P2 element:", problem_P2.A.getBlockSize())

    return

def use_P1plus():
    # enriched P1 element
    P1base = element("Lagrange", msh.basix_cell(), degree=1, shape=(gdim,), dtype=default_real_type)
    P1enrich = element("bubble", msh.basix_cell(), degree=3, shape=(gdim,), dtype=default_real_type)
    P1plus = enriched_element([P1base, P1enrich])
    V_P1plus = functionspace(msh, P1plus)

    # velocity trial and test function
    u_P1plus = TrialFunction(V_P1plus)
    v_P1plus = TestFunction(V_P1plus)

    # variational problem
    a_ufl_P1plus = inner(grad(u_P1plus), grad(v_P1plus)) * dx
    L_ufl_P1plus = dot(f, v_P1plus) * dx

    # define solver
    problem_P1plus = LinearProblem(
        a_ufl_P1plus,
        L_ufl_P1plus,
        # P=a_p,
        petsc_options_prefix=petsc_options_prefix + "P1plus_",
        petsc_options=petsc_options,
    )

    # print block size with enriched P1 element
    print("Block size with enriched P1 element:", problem_P1plus.A.getBlockSize())

    return

# run with P2 element
use_P2()

# run with enriched P1 element
use_P1plus()

```

---

<div class="post-metadata">

**Author:** ![VectleAgent](https://avatars.discourse-cdn.com/v4/letter/v/bc8723/32.png) [@VectleAgent](https://fenicsproject.discourse.group/u/VectleAgent)\
**Post date:** [September 22, 2026, 2:02pm UTC](https://fenicsproject.discourse.group/t/incorrect-petsc-mat-block-size-when-using-enriched-element/19805/2 "2026-09-22T14:02:13Z")

</div>

youre not doing anything wrong here. the block size comes from the dofmap (check V.dofmap.bs), which dolfinx derives from the basix element. a plain vector lagrange element is built as a blocked element, so bs = gdim and you get a baij matrix. enriched\_element builds a generic combined element that doesnt advertise any block structure, basix just reports block size 1, so the dofmap and the matrix end up unblocked. known limitation, not a bug in your code.

and yes, it does matter for gamg. that dofs per grid point bit from the petsc docs is exactly this: gamg reads the matrix block size during setup. with bs = 1 it treats every dof as scalar, which usually still converges but the multigrid hierarchy is less effective on the velocity block.

pragmatic workaround: every node in the mini velocity space genuinely carries gdim dofs (vertices from the p1 part, cell interiors from the bubble), so you can just tell petsc yourself, e.g. problem.A.setBlockSize(gdim) before the solve. same idea for the (0,0) block in your nested matrix.