# Stokes Operator is Singular After Applying Dirichlet Boundary Conditions

**URL:** <https://discourse.julialang.org/t/stokes-operator-is-singular-after-applying-dirichlet-boundary-conditions/135278>\
**Category:** Modelling & Simulations\
**Tags:** fem, finite-element, fluid-flow\
**Created:** [January 26, 2026, 10:14pm UTC](https://discourse.julialang.org/t/stokes-operator-is-singular-after-applying-dirichlet-boundary-conditions/135278 "2026-01-26T22:14:56Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![noetheriankoala](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/noetheriankoala/32/218462_2.png) [@noetheriankoala](https://discourse.julialang.org/u/noetheriankoala)\
**Post date:** [January 26, 2026, 10:14pm UTC](https://discourse.julialang.org/t/stokes-operator-is-singular-after-applying-dirichlet-boundary-conditions/135278/1 "2026-01-26T22:14:56Z")

</div>

Hello! I am trying to implement a simple incompressible Navier-Stokes simulation. I have a mesh consisting of a cylinder with the Stanford bunny excised from its center. When I apply Drichlet boundary conditions, the Stokes operator becomes singular, and I cannot figure out why. I did not include it in the below example, but the problem persists even if I fix one of the pressure DOFs. You can find the mesh I am using [here](https://drive.proton.me/urls/BPTTR2M2P4#dveEwlN3Rb5T).

In case it helps you help me, I am well versed in linear algebra but not in how it applies to finite element analysis. In particular I am not familiar with all the lingo.

I have attached a minimal Jupyter notebook demonstrating my issue as well as the mesh.

```julia-auto
using Ferrite
using LinearAlgebra
using SparseArrays, BlockArrays
using Ferrite
using FerriteGmsh
using FerriteGmsh: Gmsh

# load mesh
grid = togrid("/workspaces/BasisSubset/examples/advection-diffusion/bunny-cylinder.unv")

# specify velocity shape functions
ipu = Lagrange{RefTetrahedron, 2}()^3
qr = QuadratureRule{RefTetrahedron}(2)
cvu = CellValues(qr, ipu);

# specify pressure shape functions
ipp = Lagrange{RefTetrahedron, 1}() 
cvp = CellValues(qr, ipp);

# DOF handler for Navier-Stokes
dh_ns = DofHandler(grid)
add!(dh_ns, :u, ipu)
add!(dh_ns, :p, ipp)
close!(dh_ns);

# constraint handler for Navier-Stokes
ch_ns = ConstraintHandler(dh_ns);

# delete existing facesets
empty!(grid.facetsets)
empty!(grid.nodesets)

# determines if boundary is on wall of cylinder
cyl_wall_pred = x -> (abs(sqrt((x[1]+2.0)^2 + (x[2]+2.0)^2) - 50.0) < 1e-12 )

topo = ExclusiveTopology(grid)

# add facetset for everything on boundary
addboundaryfacetset!(
    grid, topo, "boundary_all",
    x -> true
)

# add facetset for everything on the wall of the cylinder
addboundaryfacetset!(
    grid, topo, "cylinder_wall", 
    x -> cyl_wall_pred(x)
)

# add facetset for everything that is not on the wall of the cylinder
addboundaryfacetset!(
    grid, topo, "boundary_other",
    x -> !cyl_wall_pred(x)
)

# WHAT I WANT
# add!(ch_ns, Dirichlet(
# :u, getfacetset(grid, "cylinder_wall"), v -> begin
# return (-v[2], v[1], 0)
# end
# ))
# add!(ch_ns, Dirichlet(
# :u, getfacetset(grid, "boundary_other"), v -> begin
# return (0.0, 0.0, 0.0)
# end
# ))

# BUT THIS DOES NOT WORK EITHER
add!(ch_ns, Dirichlet(
    :u, getfacetset(grid, "boundary_all"), v -> begin
        return (0.0, 0.0, 0.0)
    end
))

close!(ch_ns)
;

# maybe not relevant, but I noticed that the
# number of dofs in the facetsets "cylinder_wall"
# and "boundary_other" do not sum to that of 
# "boundary_all"
println(length(getfacetset(grid, "cylinder_wall")))
println(length(getfacetset(grid, "boundary_other")))
println(length(getfacetset(grid, "boundary_all")))

function assemble_stokes_operator!(K::SparseMatrixCSC, dh::DofHandler, cvu::CellValues, cvp::CellValues)
    n_basefuncs_u = getnbasefunctions(cvu)
    n_basefuncs_p = getnbasefunctions(cvp)
    n_basefuncs = n_basefuncs_u + n_basefuncs_p
    
    # will store element assembelies
    Ke = BlockedArray(
        Matrix{Float64}(undef, (n_basefuncs, n_basefuncs)), 
        [n_basefuncs_u, n_basefuncs_p], 
        [n_basefuncs_u, n_basefuncs_p]
    )

    # create assembler
    assembler_K = start_assemble(K)

    # loop over all cells
    for cell in CellIterator(dh)
        # update shape functions
        reinit!(cvu, cell)
        reinit!(cvp, cell)

        # clear element assembelies
        fill!(Ke, 0)

        # loop over quadrature points
        for q_point in 1:getnquadpoints(cvu)
            # get shape function coordinate transform
            dΩ = getdetJdV(cvu, q_point)

            # assemble velocity-velocity blocks
            for i in 1:n_basefuncs_u
                ∇a = shape_gradient(cvu, q_point, i)
                for j in 1:n_basefuncs_u
                    ∇u = shape_gradient(cvu, q_point, j)
                    Ke[BlockIndex((1,1), (i, j))] -= ∇a ⊡ ∇u * dΩ
                end
            end
            
            # assemble velocity-pressure blocks
            for j in 1:n_basefuncs_p
                p = shape_value(cvp, q_point, j)
                for i in 1:n_basefuncs_u
                    divu = shape_divergence(cvu, q_point, i)
                    Ke[BlockIndex((1,2), (i, j))] += (divu * p) * dΩ
                    Ke[BlockIndex((2,1), (j, i))] += (p * divu) * dΩ
                end
            end
        end

        # Assemble Ke into K
        assemble!(assembler_K, celldofs(cell), Ke)
    end
end

# allocate sparse matrices
K = allocate_matrix(dh_ns, ch_ns)

# allocate right hand side of solve
rhs = zeros(ndofs(dh_ns))

# asseble Stokes operator
assemble_stokes_operator!(K, dh_ns, cvu, cvp)

println("logdet before `apply!`: ", logdet(K))

# solve for velocity field
apply!(K, rhs, ch_ns)

println("logdet after `apply!`: ", logdet(K))

x0 = K \ rhs # this lines fails because K is singular
apply!(x0, ch_ns)

```

---

<div class="post-metadata">

**Author:** ![KnutAM](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/knutam/32/37720_2.png) [@KnutAM](https://discourse.julialang.org/u/KnutAM)\
**Post date:** [January 27, 2026, 2:27pm UTC](https://discourse.julialang.org/t/stokes-operator-is-singular-after-applying-dirichlet-boundary-conditions/135278/2 "2026-01-27T14:27:09Z")

</div>

> maybe not relevant, but I noticed that the  
> number of dofs in the facetsets “cylinder\_wall”  
> and “boundary\_other” do not sum to that of  
> “boundary\_all”

Answered here: [`addboundaryfacetset!` adding inconsistant number of facets · Issue #1273 · Ferrite-FEM/Ferrite.jl · GitHub](https://github.com/Ferrite-FEM/Ferrite.jl/issues/1273), but for others seeing this:  
probably need to use `all = false` in `addboundaryfacetset!` for `"boundary_other"`

---

<div class="post-metadata">

**Author:** ![ziolai](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ziolai/32/23422_2.png) [@ziolai](https://discourse.julialang.org/u/ziolai)\
**Post date:** [January 27, 2026, 3:01pm UTC](https://discourse.julialang.org/t/stokes-operator-is-singular-after-applying-dirichlet-boundary-conditions/135278/3 "2026-01-27T15:01:34Z")

</div>

Not sure.

Is inadvertedly a boundary condition on the pressure field missing, causing the pressure block to be singular?

---

<div class="post-metadata">

**Author:** ![noetheriankoala](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/noetheriankoala/32/218462_2.png) [@noetheriankoala](https://discourse.julialang.org/u/noetheriankoala)\
**Post date:** [January 27, 2026, 3:50pm UTC](https://discourse.julialang.org/t/stokes-operator-is-singular-after-applying-dirichlet-boundary-conditions/135278/4 "2026-01-27T15:50:29Z")

</div>

Even if I fix the pressure field, the operator remains singular. Since the pressure is defined up to addition by a constant, I have fixed the pressure field using

```julia-auto
addnodeset!(
    grid, "p_anchor", 
    x -> x ≈ grid.nodes[1].x
)

add!(ch_ns, Dirichlet(
    :p, getnodeset(grid, "p_anchor"), x -> 1.0
))

```

---

<div class="post-metadata">

**Author:** ![ziolai](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ziolai/32/23422_2.png) [@ziolai](https://discourse.julialang.org/u/ziolai)\
**Post date:** [January 27, 2026, 3:59pm UTC](https://discourse.julialang.org/t/stokes-operator-is-singular-after-applying-dirichlet-boundary-conditions/135278/5 "2026-01-27T15:59:04Z")

</div>

Right.

Again, not sure.

Is this error caused by nuisances in the mesh?

Should one try a 2D rectangular or 3D mesh cube mesh using the built-in Ferrite mesh generator? See e.g. [finite\_element\_electrical\_engineering/project-based-assignment/metal-hydride-storage/notebooks/stokes\_channel.ipynb at main · ziolai/finite\_element\_electrical\_engineering · GitHub](https://github.com/ziolai/finite_element_electrical_engineering/blob/main/project-based-assignment/metal-hydride-storage/notebooks/stokes_channel.ipynb)

Good luck.

---

<div class="post-metadata">

**Author:** ![noetheriankoala](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/noetheriankoala/32/218462_2.png) [@noetheriankoala](https://discourse.julialang.org/u/noetheriankoala)\
**Post date:** [January 27, 2026, 4:12pm UTC](https://discourse.julialang.org/t/stokes-operator-is-singular-after-applying-dirichlet-boundary-conditions/135278/6 "2026-01-27T16:12:17Z")

</div>

Yes, you are correct. Something is wrong with my mesh. I created a cylinder with a sphere removed from its center as below, adjusted `cyl_wall_pred`, and now everything is working as expected. Any suggestions for how I can diagnose/fix the mesh? I have the original model in FreeCAD, and used gmsh to create the `.unv` file.

Here is how I generated the test grid:

```julia-auto
function create_grid()
    gmsh.initialize()
    gmsh.clear()

    R_cyl = 1.0
    H_cyl = 5.0

    R_sph = 0.5

    lc = 0.2 # target mesh size

    cyl = gmsh.model.occ.addCylinder(
        0.0, 0.0, 0.0, # base center
        0.0, 0.0, H_cyl, # axis direction
        R_cyl
    )

    sph = gmsh.model.occ.addSphere(
        0.0, 0.0, H_cyl/2,
        R_sph
    )

    cut_entities, _ = gmsh.model.occ.cut(
        [(3, cyl)],
        [(3, sph)]
    )

    gmsh.model.occ.synchronize()

    gmsh.model.mesh.generate(3)
    
    grid = togrid()

    gmsh.finalize()

    return grid
end

grid = create_grid()

```

---

<div class="post-metadata">

**Author:** ![ziolai](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ziolai/32/23422_2.png) [@ziolai](https://discourse.julialang.org/u/ziolai)\
**Post date:** [January 27, 2026, 4:20pm UTC](https://discourse.julialang.org/t/stokes-operator-is-singular-after-applying-dirichlet-boundary-conditions/135278/7 "2026-01-27T16:20:49Z")

</div>

Nice!

Is it beneficial here to elaborate on how you define the boundary patches?

Is it wise here to perform tests in 2D first?

---

<div class="post-metadata">

**Author:** ![noetheriankoala](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/noetheriankoala/32/218462_2.png) [@noetheriankoala](https://discourse.julialang.org/u/noetheriankoala)\
**Post date:** [January 27, 2026, 5:53pm UTC](https://discourse.julialang.org/t/stokes-operator-is-singular-after-applying-dirichlet-boundary-conditions/135278/8 "2026-01-27T17:53:11Z")

</div>

I have recreated the mesh using a different Stanford bunny model and everything is now working as expected! I do not know what the issue was with the original mesh. For anyone interested, you can find the new mesh [here](https://drive.proton.me/urls/N15SHRB4HR#4DTfKrOs5K0d).

Thanks for your help!