Stokes flow (u + p)

using Rasmah

# 2D channel [0,1]×[0,1]; no-slip on the y walls, pressure-driven inlet/outlet
function channel(nx, ny)
    nv = (nx + 1) * (ny + 1)
    verts = Matrix{Float64}(undef, 3, nv)
    for j in 0:ny, i in 0:nx
        idx = i + j * (nx + 1) + 1
        verts[1, idx] = i / nx
        verts[2, idx] = j / ny
    end
    faces = Matrix{Int}(undef, 3, 2nx * ny)
    k = 0
    for j in 0:(ny - 1), i in 0:(nx - 1)
        a = i + j * (nx + 1) + 1; b = a + 1; c = a + nx + 1; d = c + 1
        k += 1; faces[:, k] = [a, b, c]
        k += 1; faces[:, k] = [b, d, c]
    end
    TriangleMesh(verts, faces)
end

m = channel(8, 8)
U = FESpace(m, ReferenceFE(2, 2); valuetype = VectorValue{2})
P = FESpace(m, ReferenceFE(2, 1))

# no-slip on y ∈ {0,1}: select the P2 wall nodes by bounding box
ufixed = vcat(dofs(U, BoundingBox([0.0, 0.0], [1.0, 0.0])),
              dofs(U, BoundingBox([0.0, 1.0], [1.0, 1.0])))

# pressure: p = 1 at inlet (x = 0), p = 0 at outlet (x = 1)
inlet  = dofs(P, BoundingBox([0.0, 0.0], [0.0, 1.0]))
outlet = dofs(P, BoundingBox([1.0, 0.0], [1.0, 1.0]))
pfixed = vcat(inlet, outlet)
pg     = vcat(fill(1.0, length(inlet)), fill(0.0, length(outlet)))

res = solve(StokesFlow(), m, water;
            bc = BoundaryConditions(FixedVelocity(ufixed), FixedPressure(pfixed; value = pg)))
println("max |u| = ", round(maximum(abs, res.u); digits = 3),
        ", p range = ", round(minimum(res.p); digits = 3), " … ",
        round(maximum(res.p); digits = 3))

vtk(m; field = res.p, fieldname = "pressure")