Hele-Shaw flow past a circle in dolfinx

I was trying to learn to simulate in fenicsx the Hele-Shaw flow past a circle based on a tutorial — Fig 1. Hele-Shaw flow past a circle · An album of computational fluid motion. But the tutorial uses older package dolfin, I am trying to recreate in dolfinx (fenicx) but getting error, and unable to understand the reason of error.

from dolfinx import fem, io, mesh, plot, cpp
from mpi4py import MPI
import numpy as np
import ufl
import gmsh
import dolfinx
from mpi4py import MPI
from dolfinx.io import (VTXWriter, distribute_entity_data, gmshio)

# Physical parameters
D = 1.0  # Circle diameter
L = 6.0  # Channel length
Ht = 3.0  # Channel height
Hb = 3.0  # Channel height
U_in = 1.0  # Inlet velocity

# Numerical parameters
resolution = 0.05  # Target characteristic length factor


EPS = np.finfo(float).eps  # Machine epsilon for float

mesh_comm = MPI.COMM_WORLD
model_rank = 0

mesh_dim=2

# Create gmsh geometry
gmsh.initialize()
# gmsh.model.occ.add("Channel with circle")

# Create channel
rectangle = gmsh.model.occ.addRectangle(-L/2, -Hb, 0, L, Ht + Hb)

# Create circle
circle = gmsh.model.occ.addDisk(0, 0, 0, D/2,D/2)


gmsh.model.occ.cut([(2, rectangle)], [(2, circle)])
gmsh.model.occ.synchronize()

# Set mesh resolution
gmsh.option.setNumber("Mesh.CharacteristicLengthFactor", resolution)
gmsh.option.setNumber("Mesh.Algorithm", 6)
gmsh.model.mesh.generate(mesh_dim)


# Convert gmsh mesh to dolfinx mesh
mesh_obj, a, b = dolfinx.io.gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=2)

# mesh_obj = dolfinx.mesh.create_mesh(gmsh.model.getMesh(False, False, dolfinx.io.gmsh_codim(gmsh.model.getDimension())))
gmsh.finalize()
gmsh.write("mesh2.msh")


# Create finite element spaces
BDM = ufl.VectorElement("Lagrange", mesh_obj.ufl_cell(), 1)
DG = ufl.FiniteElement("Discontinuous Lagrange", mesh_obj.ufl_cell(), 0)
V = fem.FunctionSpace(mesh_obj, BDM * DG)

# Define test, trial and solution function
u, p = ufl.TrialFunctions(V)
v, q = ufl.TestFunctions(V)
u_sol = fem.Function(V)

I am getting error with mesh_obj, a, b = dolfinx.io.gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=2). Can anyone please help?

Next time it would be useful if you post the whole error message:

Traceback (most recent call last):
  File "/root/shared/mwe123.py", line 49, in <module>
    mesh_obj, a, b = dolfinx.io.gmshio.model_to_mesh(
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/io/gmshio.py", line 232, in model_to_mesh
    cell_id = cell_information[perm_sort[-1]]["id"]
IndexError: index -1 is out of bounds for axis 0 with size 0

You can search for the IndexError on the forum to find

which all points you to the fact that you need to add physical groups to your mesh in GMSH to be able to read it with gmshio.model_to_mesh
https://docs.fenicsproject.org/dolfinx/v0.7.3/python/generated/dolfinx.io.gmshio.html?highlight=model_to_mesh#dolfinx.io.gmshio.model_to_mesh

Thanks. I will try your suggestions, and update in the forum. So once a mesh is created,

gmsh.model.occ.cut([(2, rectangle)], [(2, circle)])
gmsh.model.occ.synchronize()

I need to add a physical group to it before calling

dolfinx.io.gmshio.model_to_mesh

Yes (and do not call syncronize after adding the groups, as it tends to mess them up).

Thanks, @dokken . Adding the physics group actually worked. I am stuck at the last few steps, which I tried to resolve but could not understand how to proceed.

from dolfinx import fem, io, mesh, plot, cpp
from mpi4py import MPI
import numpy as np
import ufl
import gmsh
import dolfinx
from mpi4py import MPI
from dolfinx.io import (VTXWriter, distribute_entity_data, gmshio)

# Physical parameters
D = 1
Ht = 1.62 # height till top
Hb = 1.9 # height till bottom

L = 6
H = Ht+Hb
c_x = L/2
c_y = Hb
r = D/2

# Numerical parameters
resolution = 0.05  # Target characteristic length factor


EPS = np.finfo(float).eps  # Machine epsilon for float

mesh_comm = MPI.COMM_WORLD
model_rank = 0

mesh_dim=2

# Create gmsh geometry
gmsh.initialize()
# gmsh.model.occ.add("Channel with circle")

# Create channel
rectangle = gmsh.model.occ.addRectangle(0, 0, 0, L, H)

# Create circle
circle = gmsh.model.occ.addDisk(c_x, c_y, 0, r, r)


gmsh.model.occ.cut([(2, rectangle)], [(2, circle)])
# gmsh.model.occ.synchronize()



fluid_marker = 1
volumes = gmsh.model.getEntities(dim=mesh_dim)
gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], fluid_marker)
gmsh.model.setPhysicalName(volumes[0][0], fluid_marker, "Fluid")


# Set mesh resolution
gmsh.option.setNumber("Mesh.CharacteristicLengthFactor", resolution)
gmsh.option.setNumber("Mesh.Algorithm", 6)
gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2)
gmsh.option.setNumber("Mesh.RecombineAll", 1)
gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)
gmsh.model.mesh.generate(2)
gmsh.model.mesh.setOrder(2)
gmsh.model.mesh.optimize("Netgen")

gmsh.write("mesh2.msh")
# gmsh.finalize()

# Convert gmsh mesh to dolfinx mesh
mesh_obj, a, b = dolfinx.io.gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=2)


# Create finite element spaces
BDM = ufl.FiniteElement("Lagrange", mesh_obj.ufl_cell(), 1)
DG  = ufl.FiniteElement("Discontinuous Lagrange", mesh_obj.ufl_cell(), 0)

V = fem.FunctionSpace(mesh_obj, BDM*DG)

# Define test, trial and solution function
u, p = ufl.TrialFunctions(V)
v, q = ufl.TestFunctions(V)

u_sol = fem.Function(V)

# Define boundary conditions
def boundary_noflow(x):
    return x[1] < -Hb+EPS or x[1] > Ht - EPS \
        or x[0]**2+x[1]**2 < (D/2)**2+EPS
def boundary_inflow(x):
    return x[0] < -L/2+EPS

bcs = [fem.dirichletbc(V.sub(0), fem.Constant((0.0, 0.0)), boundary_noflow), \
       fem.dirichletbc(V.sub(0), fem.Constant((1.0, 0.0)), boundary_inflow),]

# Weak formulation
a = (ufl.dot(u, v) + ufl.div(v) * p + ufl.div(u) * q) * ufl.dx
Lrhs = - fem.Constant(0) * q * ufl.dx


# Solve and export
solve(a == Lrhs, u_sol, bcs)
(u,p) = u_sol.split()
u.rename('u','u')
File("sol.pvd") << u

This is the code. From

bcs = [fem.dirichletbc(V.sub(0), fem.Constant((0.0, 0.0)), boundary_noflow), \
       fem.dirichletbc(V.sub(0), fem.Constant((1.0, 0.0)), boundary_inflow),]

# Weak formulation
a = (ufl.dot(u, v) + ufl.div(v) * p + ufl.div(u) * q) * ufl.dx
Lrhs = - fem.Constant(0) * q * ufl.dx


# Solve and export
solve(a == Lrhs, u_sol, bcs)
(u,p) = u_sol.split()
u.rename('u','u')
File("sol.pvd") << u

This section of code I have issue with. I am following old code as provided -
Fig 1. Hele-Shaw flow past a circle · An album of computational fluid motion. I am getting error bcs and Lrhs. The solve function too I don’t undertand.

  Cell In[7], line 1
    bcs = [fem.dirichletbc(V.sub(0), fem.Constant((0.0, 0.0)), boundary_noflow), \

TypeError: Constant.__init__() missing 1 required positional argument: 'c'

< minor highjack of thread >
Funny I should see this. Do you think you could update the code on the website once you have the fenicsx version running? That would be valuable for future readers.
Maybe the best approach would be to add your codeblock underneath mine. The website is opensource (using typical git workflow of forking, editing, merge-requesting). You can find the .md file on: content/chapters/01-creeping/Fig1/index.md · main · Stein Stoter / album-of-computational-fluid-motion · GitLab.
< /minor highjack of thread >

Yes, Definitely. That’s a good learning exercise for me too.

1 Like

I would advice you to have a look at the API of DOLFINx:
dolfinx.fem — DOLFINx 0.7.3 documentation which states that you first need to send in the mesh, ie mesh_obj.

Secondly, for the solve command, see any of the demos or the DOLFINx tutorial for the updated API

from dolfinx import fem, io, mesh, plot, cpp
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import ufl
import gmsh
import dolfinx
from mpi4py import MPI
from dolfinx.io import (VTXWriter, distribute_entity_data, gmshio)
from dolfinx.fem.petsc import LinearProblem


# Physical parameters
D = 1
Ht = 1.62 # height till top
Hb = 1.9 # height till bottom

L = 6
H = Ht+Hb
c_x = 0
c_y = 0
r = D/2

# Numerical parameters
resolution = 0.05  # Target characteristic length factor


EPS = np.finfo(float).eps  # Machine epsilon for float

mesh_comm = MPI.COMM_WORLD
model_rank = 0

mesh_dim=2

# Create gmsh geometry
gmsh.initialize()

# Create channel
rectangle = gmsh.model.occ.addRectangle(-L/2, -Hb, 0, L, H)

# Create circle
circle = gmsh.model.occ.addDisk(c_x, c_y, 0, r, r)


gmsh.model.occ.cut([(2, rectangle)], [(2, circle)])
gmsh.model.occ.synchronize()



fluid_marker = 1
volumes = gmsh.model.getEntities(dim=mesh_dim)
gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], fluid_marker)
gmsh.model.setPhysicalName(volumes[0][0], fluid_marker, "Fluid")


# Set mesh resolution
gmsh.option.setNumber("Mesh.CharacteristicLengthFactor", resolution)
gmsh.option.setNumber("Mesh.Algorithm", 8)
# gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 1)
# gmsh.option.setNumber("Mesh.RecombineAll", 1)
# gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)
gmsh.model.mesh.generate(2)
# gmsh.model.mesh.setOrder(2)
gmsh.model.mesh.optimize("Netgen")

# gmsh.write("mesh2.msh")

# Convert gmsh mesh to dolfinx mesh
mesh_obj, _, _ = dolfinx.io.gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=2)
gmsh.finalize()

# Create finite element spaces
BDM = ufl.FiniteElement("BDM", mesh_obj.ufl_cell(), 1)
DG  = ufl.FiniteElement("DG", mesh_obj.ufl_cell(), 0)

V = fem.FunctionSpace(mesh_obj, BDM*DG)

# V = fem.FunctionSpace(mesh_obj, BDM)
# Q = fem.FunctionSpace(mesh_obj, DG)

# Define test, trial and solution function
u, p = ufl.TrialFunctions(V)
v, q = ufl.TestFunctions(V)

u_sol = fem.Function(V)

# Define boundary conditions
def boundary_noflow(x):
    return np.logical_or(x[1] < -Hb+EPS, x[1] > Ht-EPS, x[0]**2 + x[1]**2 < (D/2)**2 + EPS)
 
def boundary_inflow(x):
    return np.isclose(x[0], -L/2)

fdim = mesh_obj.topology.dim - 1

V1, _ = V.sub(0).collapse()
u_D = fem.Function(V1)
u_D.x.array[:] =1


facets_noflow = mesh.locate_entities_boundary(mesh_obj, fdim, boundary_noflow)
bc_noflow = fem.dirichletbc(u_D, fem.locate_dofs_topological((V.sub(0), V1), 1, facets_noflow), V.sub(0))

facets_inflow = mesh.locate_entities_boundary(mesh_obj, 1, boundary_inflow)
bc_inflow = fem.dirichletbc(u_D, fem.locate_dofs_topological((V.sub(0), V1), 1, facets_inflow), V.sub(0))


# inflow_dofs = fem.locate_dofs_geometrical(V, boundary_inflow)
# bc_inflow = fem.dirichletbc(PETSc.ScalarType((1, 0)), inflow_dofs)


# Weak formulation
a = (ufl.dot(u, v) + ufl.div(v) * p + ufl.div(u) * q) * ufl.dx
Lrhs = - fem.Constant(mesh_obj, PETSc.ScalarType(0)) * q * ufl.dx

# Solve and export
problem_loc = LinearProblem(a, Lrhs,  [bc_noflow, bc_inflow],u_sol)
# uh = problem_loc.solve()
# u_sol.interpolate(uh)
(u,p) = u_sol.split()


with io.XDMFFile(mesh_obj.comm, "Out_Hele_Shaw_flow/Hele-Shaw-flow.xdmf", "w") as file:
    file.write_mesh(mesh_obj)
    file.write_function(u, 0)

This is the final code, but some error is coming, which I am unable to correct.

    file.write_function(u, 0)

  File /usr/lib/petsc/lib/python3/dist-packages/dolfinx/io/utils.py:235 in write_function
    super().write_function(getattr(u, "_cpp_object", u), t, mesh_xpath)

RuntimeError: Only Lagrange functions are supported. Interpolate Functions before output.

See Interpolating Functions on higher order elements - #2 by dokken

The problem is mostly solved. I am able to export uh_p, but I am having issue exporting the velocity field, uh_u. It would be greatly appreciated if anyone could help with export.

from dolfinx import fem, io, mesh, plot, cpp
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import ufl
import gmsh
import dolfinx
from mpi4py import MPI
from dolfinx.io import (VTXWriter, distribute_entity_data, gmshio)
from dolfinx.fem.petsc import LinearProblem
from basix.ufl import element, mixed_element


# Physical parameters
D = 1
Ht = 1.62 # height till top
Hb = 1.9 # height till bottom

L = 6
H = Ht+Hb
c_x = 0
c_y = 0
r = D/2

# Numerical parameters
resolution = 0.05  # Target characteristic length factor


EPS = np.finfo(float).eps  # Machine epsilon for float

mesh_comm = MPI.COMM_WORLD
model_rank = 0

mesh_dim=2

# Create gmsh geometry
gmsh.initialize()

# Create channel
rectangle = gmsh.model.occ.addRectangle(-L/2, -Hb, 0, L, H)

# Create circle
circle = gmsh.model.occ.addDisk(c_x, c_y, 0, r, r)


gmsh.model.occ.cut([(2, rectangle)], [(2, circle)])
gmsh.model.occ.synchronize()



fluid_marker = 1
volumes = gmsh.model.getEntities(dim=mesh_dim)
gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], fluid_marker)
gmsh.model.setPhysicalName(volumes[0][0], fluid_marker, "Fluid")


# Set mesh resolution
gmsh.option.setNumber("Mesh.CharacteristicLengthFactor", resolution)
gmsh.option.setNumber("Mesh.Algorithm", 8)
# gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 1)
# gmsh.option.setNumber("Mesh.RecombineAll", 1)
# gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)
gmsh.model.mesh.generate(2)
# gmsh.model.mesh.setOrder(2)
gmsh.model.mesh.optimize("Netgen")

# gmsh.write("mesh2.msh")

# Convert gmsh mesh to dolfinx mesh
mesh_obj, _, _ = dolfinx.io.gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=2)
gmsh.finalize()


# with io.XDMFFile(MPI.COMM_WORLD, "problem_interval_mesh.xdmf", "w") as xdmf:
#       xdmf.write_mesh(mesh_obj)

# Create finite element spaces
BDM = ufl.FiniteElement("BDM", mesh_obj.ufl_cell(), 1)
DG  = ufl.FiniteElement("DG", mesh_obj.ufl_cell(), 0)
V = fem.FunctionSpace(mesh_obj, BDM*DG)

# V = fem.FunctionSpace(mesh_obj, BDM)
# Q = fem.FunctionSpace(mesh_obj, DG)

# Define test, trial and solution functshoes that fits feet designion
(u, p) = ufl.TrialFunctions(V)
(v, q) = ufl.TestFunctions(V)

u_sol = fem.Function(V)


# Define boundary conditions
def boundary_noflow(x):
    return np.logical_or(np.isclose(x[1], -Hb+EPS), np.isclose(x[1], Ht-EPS), np.isclose(x[0]**2 + x[1]**2, r**2 + EPS))
 
def boundary_inflow(x):
    return np.isclose(x[0], -L/2)

fdim = mesh_obj.topology.dim - 1

V0 = V.sub(0)
Q, _ = V0.collapse()

# Noflow
facets_noflow = mesh.locate_entities_boundary(mesh_obj, fdim, boundary_noflow)
dofs_noflow = fem.locate_dofs_topological((V0, Q), fdim, facets_noflow)

def f1(x):
    values = np.zeros((2, x.shape[1]))
    return values

f_h1 = fem.Function(Q)
f_h1.interpolate(f1)

bc_noflow = fem.dirichletbc(f_h1, dofs_noflow, V0)

#Inflow
facets_inflow = mesh.locate_entities_boundary(mesh_obj, fdim, boundary_inflow)
dofs_inflow = fem.locate_dofs_topological((V0, Q), fdim, facets_inflow)

def f2(x):
    values = np.zeros((2, x.shape[1]))
    values[0, :] = 1
    return values

f_h2 = fem.Function(Q)
f_h2.interpolate(f2)

bc_inflow = fem.dirichletbc(f_h2, dofs_inflow, V0)

bcs = [bc_noflow, bc_inflow]

# inflow_dofs = fem.locate_dofs_geometrical(V, boundary_inflow)
# bc_inflow = fem.dirichletbc(PETSc.ScalarType((1, 0)), inflow_dofs)


# Weak formulation
a = (ufl.inner(u, v) + ufl.div(v) * p + ufl.div(u) * q) * ufl.dx
# Lrhs = - fem.Constant(mesh_obj, PETSc.ScalarType(0)) * q * ufl.dx

Lrhs = - fem.Constant(mesh_obj, PETSc.ScalarType(0)) * q * ufl.dx

# Solve and export
# problem_loc = LinearProblem(a, Lrhs,  bcs)
# uh = problem_loc.solve()
# uh_u, uh_p = uh.split()
# u, p = u_sol.split()
# u.interpolate(uh_u)
# p.interpolate(uh_p)




# problem = LinearProblem(
#     a,
#     Lrhs,
#     bcs=bcs,
#     petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"},
# )
# try:
#     uh = problem.solve()
# except PETSc.Error as e:  # type: ignore
#     if e.ierr == 92:
#         print("The required PETSc solver/preconditioner is not available. Exiting.")
#         print(e)
#         exit(0)
#     else:
#         raise e


problem_loc = LinearProblem(a, Lrhs,  bcs)

uh = problem_loc.solve()
uh_u, uh_p = uh.split()



# W0_DG = fem.FunctionSpace(mesh_obj, ("BDM", 1))
# sigma_DG = fem.Function(W0_DG)
# sigma_DG.interpolate(uh_u)


# uh_u_expr = fem.Expression(uh_u, V1.element.interpolation_points())
# u.interpolate(uh_u_expr)
# u.x.scatter_forward()

# from pathlib import Path
# name = "DFG_navier_stokes" 
# folder = Path("results")

# with io.VTXWriter(mesh_obj.comm, folder / f"{name}_velocity_final.bp", uh_p,  engine="BP4") as f:
#     f.write(0)

with io.XDMFFile(mesh_obj.comm, "out_mixed_poisson/u.xdmf", "w") as file:
    file.write_mesh(mesh_obj)
    file.write_function(uh_p)
    

# with dolfinx.io.VTXWriter(mesh_obj.comm, "moisture_solution1D.bp", uh_u, engine="BP4") as vtx:
#     vtx.write(0.0)
    

Reference:

[1] Mixed formulation for the Poisson equation - DOLFINx Tutorial

You are not mentioning what issue you are having.
You need to interpolate uh_u in to a suitable Lagrange (or DG) space for visualization in paraview with VTXWriter.

Suppose If I run for uh_u to XDMFFILE

with io.XDMFFile(mesh_obj.comm, "out_mixed_poisson/u.xdmf", "w") as file:
    file.write_mesh(mesh_obj)
    file.write_function(uh_u)

I am getting this error

 Cell In[2], line 3
    file.write_function(uh_u)

  File /usr/lib/petsc/lib/python3/dist-packages/dolfinx/io/utils.py:235 in write_function
    super().write_function(getattr(u, "_cpp_object", u), t, mesh_xpath)

RuntimeError: Only Lagrange functions are supported. Interpolate Functions before output.

As I said, you need to interpolate uh_u into a suitable (discontinuous) Lagrange space.
For XDMFFile this means a first order Lagrange space if the geometry is linear.
For VTXWriter you can use an arbitrary order (discontinuous) Lagrange space for visualization.

1 Like

Thanks @dokken. @Stein

I suppose you are asking for this

gdim = mesh_obj.geometry.dim
V1 = fem.functionspace(mesh_obj, ("Discontinuous Lagrange", 1, (gdim,)))
u1 = fem.Function(V1, dtype=default_scalar_type)
u1.interpolate(uh_u)

Complete code:

from dolfinx import fem, io, mesh, plot, cpp
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import ufl
import gmsh
import dolfinx
from mpi4py import MPI
from dolfinx.io import (VTXWriter, distribute_entity_data, gmshio)
from dolfinx.fem.petsc import LinearProblem
from basix.ufl import element, mixed_element
from dolfinx import default_scalar_type, plot


# Physical parameters
D = 1
Ht = 1.62 # height till top
Hb = 1.9 # height till bottom

L = 6
H = Ht+Hb
c_x = 0
c_y = 0
r = D/2

# Numerical parameters
resolution = 0.05  # Target characteristic length factor


EPS = np.finfo(float).eps  # Machine epsilon for float

mesh_comm = MPI.COMM_WORLD
model_rank = 0

mesh_dim=2

# Create gmsh geometry
gmsh.initialize()

# Create channel
rectangle = gmsh.model.occ.addRectangle(-L/2, -Hb, 0, L, H)

# Create circle
circle = gmsh.model.occ.addDisk(c_x, c_y, 0, r, r)


gmsh.model.occ.cut([(2, rectangle)], [(2, circle)])
gmsh.model.occ.synchronize()



fluid_marker = 1
volumes = gmsh.model.getEntities(dim=mesh_dim)
gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], fluid_marker)
gmsh.model.setPhysicalName(volumes[0][0], fluid_marker, "Fluid")


# Set mesh resolution
gmsh.option.setNumber("Mesh.CharacteristicLengthFactor", resolution)
gmsh.option.setNumber("Mesh.Algorithm", 8)
# gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 1)
# gmsh.option.setNumber("Mesh.RecombineAll", 1)
# gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)
gmsh.model.mesh.generate(2)
# gmsh.model.mesh.setOrder(2)
gmsh.model.mesh.optimize("Netgen")

# gmsh.write("mesh2.msh")

# Convert gmsh mesh to dolfinx mesh
mesh_obj, _, _ = dolfinx.io.gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=2)
gmsh.finalize()


# with io.XDMFFile(MPI.COMM_WORLD, "problem_interval_mesh.xdmf", "w") as xdmf:
#       xdmf.write_mesh(mesh_obj)

# Create finite element spaces
BDM = ufl.FiniteElement("BDM", mesh_obj.ufl_cell(), 1)
DG  = ufl.FiniteElement("DG", mesh_obj.ufl_cell(), 0)
V = fem.FunctionSpace(mesh_obj, BDM*DG)

# V = fem.FunctionSpace(mesh_obj, BDM)
# Q = fem.FunctionSpace(mesh_obj, DG)

# Define test, trial and solution functshoes that fits feet designion
(u, p) = ufl.TrialFunctions(V)
(v, q) = ufl.TestFunctions(V)

u_sol = fem.Function(V)


# Define boundary conditions
def boundary_noflow(x):
    return np.logical_or(np.isclose(x[1], -Hb+EPS), np.isclose(x[1], Ht-EPS), np.isclose(x[0]**2 + x[1]**2, r**2 + EPS))
 
def boundary_inflow(x):
    return np.isclose(x[0], -L/2)

fdim = mesh_obj.topology.dim - 1

V0 = V.sub(0)
Q, _ = V0.collapse()

# Noflow
facets_noflow = mesh.locate_entities_boundary(mesh_obj, fdim, boundary_noflow)
dofs_noflow = fem.locate_dofs_topological((V0, Q), fdim, facets_noflow)

def f1(x):
    values = np.zeros((2, x.shape[1]))
    return values

f_h1 = fem.Function(Q)
f_h1.interpolate(f1)

bc_noflow = fem.dirichletbc(f_h1, dofs_noflow, V0)

#Inflow
facets_inflow = mesh.locate_entities_boundary(mesh_obj, fdim, boundary_inflow)
dofs_inflow = fem.locate_dofs_topological((V0, Q), fdim, facets_inflow)

def f2(x):
    values = np.zeros((2, x.shape[1]))
    values[0, :] = 1
    return values

f_h2 = fem.Function(Q)
f_h2.interpolate(f2)

bc_inflow = fem.dirichletbc(f_h2, dofs_inflow, V0)

bcs = [bc_noflow, bc_inflow]

# inflow_dofs = fem.locate_dofs_geometrical(V, boundary_inflow)
# bc_inflow = fem.dirichletbc(PETSc.ScalarType((1, 0)), inflow_dofs)


# Weak formulation
a = (ufl.inner(u, v) + ufl.div(v) * p + ufl.div(u) * q) * ufl.dx
# Lrhs = - fem.Constant(mesh_obj, PETSc.ScalarType(0)) * q * ufl.dx

Lrhs = - fem.Constant(mesh_obj, PETSc.ScalarType(0)) * q * ufl.dx

# Solve and export
# problem_loc = LinearProblem(a, Lrhs,  bcs)
# uh = problem_loc.solve()
# uh_u, uh_p = uh.split()
# u, p = u_sol.split()
# u.interpolate(uh_u)
# p.interpolate(uh_p)




# problem = LinearProblem(
#     a,
#     Lrhs,
#     bcs=bcs,
#     petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"},
# )
# try:
#     uh = problem.solve()
# except PETSc.Error as e:  # type: ignore
#     if e.ierr == 92:
#         print("The required PETSc solver/preconditioner is not available. Exiting.")
#         print(e)
#         exit(0)
#     else:
#         raise e


problem_loc = LinearProblem(a, Lrhs,  bcs)

uh = problem_loc.solve()
uh_u, uh_p = uh.split()



# W0_DG = fem.FunctionSpace(mesh_obj, ("BDM", 1))
# sigma_DG = fem.Function(W0_DG)
# uh_u.interpolate(sigma_DG)



gdim = mesh_obj.geometry.dim
V1 = fem.functionspace(mesh_obj, ("Discontinuous Lagrange", 1, (gdim,)))
u1 = fem.Function(V1, dtype=default_scalar_type)
u1.interpolate(uh_u)


# uh_u_expr = fem.Expression(uh_u, V1.element.interpolation_points())
# u.interpolate(uh_u_expr)
# u.x.scatter_forward()

# from pathlib import Path
# name = "DFG_navier_stokes" 
# folder = Path("results")

# with io.VTXWriter(mesh_obj.comm, folder / f"{name}_velocity_final.bp", uh_p,  engine="BP4") as f:
#     f.write(0)

# with io.XDMFFile(mesh_obj.comm, "out_mixed_poisson/u.xdmf", "w") as file:
#     file.write_mesh(mesh_obj)
#     file.write_function(uh_u)
    

with dolfinx.io.VTXWriter(mesh_obj.comm, "moisture_solution1D.bp", u1, engine="BP4") as vtx:
    vtx.write(0.0)

I have to check the values though for the correctness.