Hello,

I’m seeking assistance with a specific section of my code. Here’s the relevant snippet:

V = VectorL2(mesh, order=order)

Qi = H1(mesh, order=order, dirichlet=“gamma1|gamma2”)

QI = Compress(Qi, active_dofs=Qi.GetDofs(mesh.Boundaries(“gammaI”)))

X = FESpace([V, Q, Qbar, QI])

(u, p, pbar, pI), (v, q, qbar, qI) = X.TnT()

dx1 = dx(definedon=mesh.Materials(“Domain1”))

ds1 = dx(element_boundary=True, definedon=mesh.Materials(“Domain1”))

dx2 = dx(definedon=mesh.Materials(“Domain2”))

ds2 = dx(element_boundary=True, definedon=mesh.Materials(“Domain2”))

dsI = ds(definedon=mesh.Boundaries(“gammaI”))

dsI2 = ds(element_boundary=True, definedon=mesh.Boundaries(“gammaI”))

dsI3 = ds(skeleton=True, definedon=mesh.Boundaries(“gammaI”))

# Bilinear form terms with pressure in the interface

…

a += pI.Deriv().Trace() * qI.Deriv().Trace() * dsI2

a += jump_u * v * n * dsI2

a += av_u * v * n * dsI2

a += jump_u * v.Other() * n * dsI2

a += -av_u * v.Other() * n * dsI3

a += pI * jump_v * dsI3

a += qI * jump_u * dsI3

…

Order=1 I am not getting convergence, and for order>1, I encounter the following error: “UmfpackInverse: Numeric factorization failed.”

I suspect the issue lies within the aforementioned terms of the bilinear form, specifically regarding the integration differences between dsI, dsI2, and dsI3. Could someone please shed some light on this problem?

I would greatly appreciate any assistance provided. Thank you!