KVSlab
FEniCSx transition
Describe new/missing feature
Hello, thank you very much for the effort in maintaining this nice repository to do FSI. I wonder whether a transition to dolfinx, thus moving away from legacy DOLFIN, is going to happen soon or is it at least planned?
Thank you very much.
Suggestion user interface
Ciao Federico,
Thanks a lot for your interest in turtleFSI. We truly appreciate you reaching out. You're right to ask about the transition; it's something we're excited about and really want to do as soon as possible. While we are currently seeking the necessary funding, we'd welcome your thoughts on the subject. If you have any innovative ideas or suggestions, we'd love to hear them.
Best, Kristian
On Fri, 17 Oct 2025 at 14:37, Federico Betti @.***> wrote:
federicobetti-ethz created an issue (KVSlab/turtleFSI#113) https://github.com/KVSlab/turtleFSI/issues/113 Describe new/missing feature
Hello, thank you very much for the effort in maintaining this nice repository to do FSI. I wonder whether a transition to dolfinx, thus moving away from legacy DOLFIN, is going to happen soon or is it at least planned?
Thank you very much. Suggestion user interface
— Reply to this email directly, view it on GitHub https://github.com/KVSlab/turtleFSI/issues/113, or unsubscribe https://github.com/notifications/unsubscribe-auth/AA4X54NPSRJB5ZZW5DH67FL3YDPIDAVCNFSM6AAAAACJPFMB5OVHI2DSMVQWIX3LMV43ASLTON2WKOZTGUZDMMBQGM2DCNQ . You are receiving this because you are subscribed to this thread.Message ID: @.***>
Hello Kristian, thank you very much for your kind reply. I am currently trying to translate the Turek Hron benchmark to dolfinx. The code I wrote runs just fine (and it was a reasonable effort to be honest) but I get incorrect values of drag and lift, so I wonder if I missed something that is a bit hidden. But, the actual transition to dolfinx was smooth I would say.
The main things are:
- changing all SubDomains to mesh tags given by dolfinx
- changing the newton solver and the application of boundary conditions
- changing the loading of the mesh
Again, perhaps I am still missing something as I get different values for drag and lift than the ones mentioned in the original benchmark paper. What are your thoughts? Would be great to hear from you.
Hello Kristian, thank you very much for your kind reply. I am currently trying to translate the Turek Hron benchmark to dolfinx. The code I wrote runs just fine (and it was a reasonable effort to be honest) but I get incorrect values of drag and lift, so I wonder if I missed something that is a bit hidden. But, the actual transition to dolfinx was smooth I would say.
The main things are:
1. changing all SubDomains to mesh tags given by dolfinx 2. changing the newton solver and the application of boundary conditions 3. changing the loading of the meshAgain, perhaps I am still missing something as I get different values for drag and lift than the ones mentioned in the original benchmark paper. What are your thoughts? Would be great to hear from you.
Without the actually translated code, it is a bit hard to pinpoint what your particular issue is.
That being said; doing a 1 to 1 translation of turtleFSI shouldn't be to tiresome to do. However, there are many design choices that should be re-visited properly (for instance the whole API). This would require some sketches of what functionality one would like, etc.
Hi @jorgensd, thank you for your reply. I attach below the two functions I wrote for translating in dolfinx the creation of the meshtags and bcs for the Turek Hron benchmark.
def setup_mesh_and_tags(R, H, L, f_L, f_H, c_x, c_y):
"""Setup mesh and create boundary/domain tags for FSI problem.
Args:
R: Radius of circle
H: Channel height
L: Channel length
f_L: Flag length
f_H: Flag height
c_x: Circle center x-coordinate
c_y: Circle center y-coordinate
Returns:
mesh: The computational mesh
domain_tag: Tags for fluid/solid domains
facet_tag: Tags for boundaries
"""
filepath = path.join(
path.dirname(path.dirname(path.abspath(__file__))), "mesh", "mesh.xdmf"
)
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, filepath, "r") as xdmf:
mesh = xdmf.read_mesh(name="Grid")
def inlet_boundary(x):
return np.isclose(x[0], 0)
def outlet_boundary(x):
return np.isclose(x[0], L)
def wall_boundary(x):
return np.logical_or(np.isclose(x[1], H), np.isclose(x[1], 0))
def bar_boundary(x):
return np.logical_or(
np.logical_or(
np.isclose(x[1], c_y + f_H / 2), np.isclose(x[1], c_y - f_H / 2)
),
np.isclose(x[0], c_x + R + f_L),
)
def circle_boundary(x):
return (x[0] - c_x) ** 2 + (x[1] - c_y) ** 2 < R**2 + DOLFIN_EPS * 1e5
def barwall_boundary(x):
belongs_to_circle = circle_boundary(x)
belongs_to_bar = np.logical_and(x[1] <= c_y + f_H / 2, x[1] >= c_y - f_H / 2)
x_larger_than_bar = x[0] > c_x
return np.logical_and(
belongs_to_circle, np.logical_and(belongs_to_bar, x_larger_than_bar)
)
boundaries = [
(2, wall_boundary),
(3, inlet_boundary),
(4, outlet_boundary),
(5, bar_boundary),
(6, circle_boundary),
(7, barwall_boundary),
]
facet_indices, facet_markers = [], []
facetdim = mesh.topology.dim - 1
for marker, locator in boundaries:
boundary_facets = dolfinx.mesh.locate_entities(mesh, facetdim, locator)
facet_indices.append(boundary_facets)
facet_markers.append(np.full_like(boundary_facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = dolfinx.mesh.meshtags(
mesh, facetdim, facet_indices[sorted_facets], facet_markers[sorted_facets]
)
def solid_interior(x):
flag_y_region = np.logical_and(x[1] <= c_y + f_H / 2, x[1] >= c_y - f_H / 2)
flag_x_region = np.logical_and(x[0] >= c_x, x[0] <= c_x + R + f_L)
flag_region = np.logical_and(flag_y_region, flag_x_region)
return flag_region
interiors = [
(2, solid_interior),
]
domain_indices, domain_markers = [], []
vdim = mesh.topology.dim
all_cells = np.arange(mesh.topology.index_map(vdim).size_local)
for marker, locator in interiors:
interior_cells = dolfinx.mesh.locate_entities(mesh, vdim, locator)
remaining_cells = np.setdiff1d(all_cells, interior_cells)
domain_indices.append(interior_cells)
domain_markers.append(np.full_like(interior_cells, marker, dtype=np.int32))
domain_indices.append(remaining_cells)
domain_markers.append(np.full_like(remaining_cells, 1, dtype=np.int32))
domain_indices = np.hstack(domain_indices).astype(np.int32)
domain_markers = np.hstack(domain_markers).astype(np.int32)
sorted_domains = np.argsort(domain_indices)
domain_tag = dolfinx.mesh.meshtags(
mesh, vdim, domain_indices[sorted_domains], domain_markers[sorted_domains]
)
mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim - 1)
mesh_dir = path.dirname(path.dirname(path.abspath(__file__)))
facet_tag_file = path.join(mesh_dir, "mesh", "facet_tags.xdmf")
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, facet_tag_file, "w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_meshtags(facet_tag)
domain_tag_file = path.join(mesh_dir, "mesh", "domain_tags.xdmf")
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, domain_tag_file, "w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_meshtags(domain_tag)
return mesh, domain_tag, facet_tag
def setup_function_spaces(mesh, domain_tag, facet_tag, Um, H, dt):
"""Setup function spaces, measures, and boundary conditions.
Args:
mesh: The computational mesh
domain_tag: Tags for fluid/solid domains
facet_tag: Tags for boundaries
Um: Maximum inlet velocity
H: Channel height
dt: Time step
Returns:
Tuple containing:
- dx, ds, dS: Measures for domain, boundary, interior facets
- n: Facet normal
- bcs: List of boundary conditions
- DVP: Mixed function space
- d_, v_, p_, w_: Field dictionaries
- dvp_: Solution function dictionary
- k: Time step constant
- phi, psi, gamma, beta: Test functions
"""
ds = ufl.Measure("ds", domain=mesh, subdomain_data=facet_tag)
dS = ufl.Measure("dS", domain=mesh, subdomain_data=facet_tag)
dx = ufl.Measure("dx", domain=mesh, subdomain_data=domain_tag)
de = ufl.VectorElement("CG", mesh.ufl_cell(), 2)
ve = ufl.VectorElement("CG", mesh.ufl_cell(), 2)
pe = ufl.FiniteElement("CG", mesh.ufl_cell(), 1)
k = dolfinx.fem.Constant(mesh, c=dt)
DVP_elem = ufl.MixedElement(de, ve, pe, de)
DVP = dolfinx.fem.FunctionSpace(mesh, DVP_elem)
de_space = dolfinx.fem.FunctionSpace(mesh, de)
ve_space = dolfinx.fem.FunctionSpace(mesh, ve)
pe_space = dolfinx.fem.FunctionSpace(mesh, pe)
d_ = {}
v_ = {}
p_ = {}
w_ = {}
dvp_ = {}
times = ["n-2", "n-1", "n"]
for time_ in times:
dvp = dolfinx.fem.Function(DVP)
dvp_[time_] = dvp
dvp_list = ufl.split(dvp)
d_[time_] = dvp_list[0]
v_[time_] = dvp_list[1]
p_[time_] = dvp_list[2]
w_[time_] = dvp_list[3]
phi, psi, gamma, beta = ufl.TestFunctions(DVP)
time_constant = dolfinx.fem.Constant(mesh, c=0.5)
def inlet_profile(x):
u_x = 1.5 * Um * x[1] * (H - x[1]) / (H / 2.0) ** 2
u_y = np.zeros_like(x[1])
return time_constant.value * np.stack([u_x, u_y])
inlet_profile_func = dolfinx.fem.Function(ve_space)
inlet_profile_func.interpolate(inlet_profile)
bcs = []
inlet_facets = facet_tag.find(3)
inlet_dofs = dolfinx.fem.locate_dofs_topological(
(DVP.sub(1), ve_space), mesh.topology.dim - 1, inlet_facets
)
u_inlet = dolfinx.fem.dirichletbc(inlet_profile_func, inlet_dofs, DVP.sub(1))
wall_facets = facet_tag.find(2)
wall_dofs = dolfinx.fem.locate_dofs_topological(
(DVP.sub(1), ve_space), mesh.topology.dim - 1, wall_facets
)
zero_v_func = dolfinx.fem.Function(ve_space)
u_wall = dolfinx.fem.dirichletbc(zero_v_func, wall_dofs, DVP.sub(1))
circle_facets = facet_tag.find(6)
circle_dofs = dolfinx.fem.locate_dofs_topological(
(DVP.sub(1), ve_space), mesh.topology.dim - 1, circle_facets
)
u_circle = dolfinx.fem.dirichletbc(zero_v_func, circle_dofs, DVP.sub(1))
barwall_facets = facet_tag.find(7)
barwall_dofs = dolfinx.fem.locate_dofs_topological(
(DVP.sub(1), ve_space), mesh.topology.dim - 1, barwall_facets
)
u_barwall = dolfinx.fem.dirichletbc(zero_v_func, barwall_dofs, DVP.sub(1))
outlet_facets = facet_tag.find(4)
outlet_dofs = dolfinx.fem.locate_dofs_topological(
(DVP.sub(2), pe_space), mesh.topology.dim - 1, outlet_facets
)
zero_p_func = dolfinx.fem.Function(pe_space)
p_out = dolfinx.fem.dirichletbc(zero_p_func, outlet_dofs, DVP.sub(2))
bcs = [u_wall, u_inlet, u_circle, u_barwall, p_out]
zero_d_func = dolfinx.fem.Function(de_space)
wall_dofs = dolfinx.fem.locate_dofs_topological(
(DVP.sub(0), de_space), mesh.topology.dim - 1, wall_facets
)
d_wall = dolfinx.fem.dirichletbc(zero_d_func, wall_dofs, DVP.sub(0))
inlet_dofs = dolfinx.fem.locate_dofs_topological(
(DVP.sub(0), de_space), mesh.topology.dim - 1, inlet_facets
)
d_inlet = dolfinx.fem.dirichletbc(zero_d_func, inlet_dofs, DVP.sub(0))
outlet_dofs = dolfinx.fem.locate_dofs_topological(
(DVP.sub(0), de_space), mesh.topology.dim - 1, outlet_facets
)
d_outlet = dolfinx.fem.dirichletbc(zero_d_func, outlet_dofs, DVP.sub(0))
circle_dofs = dolfinx.fem.locate_dofs_topological(
(DVP.sub(0), de_space), mesh.topology.dim - 1, circle_facets
)
d_circle = dolfinx.fem.dirichletbc(zero_d_func, circle_dofs, DVP.sub(0))
barwall_dofs = dolfinx.fem.locate_dofs_topological(
(DVP.sub(0), de_space), mesh.topology.dim - 1, barwall_facets
)
d_barwall = dolfinx.fem.dirichletbc(zero_d_func, barwall_dofs, DVP.sub(0))
zero_w_func = dolfinx.fem.Function(de_space)
wall_dofs = dolfinx.fem.locate_dofs_topological(
(DVP.sub(3), de_space), mesh.topology.dim - 1, wall_facets
)
w_wall = dolfinx.fem.dirichletbc(zero_w_func, wall_dofs, DVP.sub(3))
inlet_dofs = dolfinx.fem.locate_dofs_topological(
(DVP.sub(3), de_space), mesh.topology.dim - 1, inlet_facets
)
w_inlet = dolfinx.fem.dirichletbc(zero_w_func, inlet_dofs, DVP.sub(3))
outlet_dofs = dolfinx.fem.locate_dofs_topological(
(DVP.sub(3), de_space), mesh.topology.dim - 1, outlet_facets
)
w_outlet = dolfinx.fem.dirichletbc(zero_w_func, outlet_dofs, DVP.sub(3))
circle_dofs = dolfinx.fem.locate_dofs_topological(
(DVP.sub(3), de_space), mesh.topology.dim - 1, circle_facets
)
w_circle = dolfinx.fem.dirichletbc(zero_w_func, circle_dofs, DVP.sub(3))
barwall_dofs = dolfinx.fem.locate_dofs_topological(
(DVP.sub(3), de_space), mesh.topology.dim - 1, barwall_facets
)
w_barwall = dolfinx.fem.dirichletbc(zero_w_func, barwall_dofs, DVP.sub(3))
bcs.extend(
[
d_wall,
d_inlet,
d_outlet,
d_circle,
d_barwall,
w_wall,
w_inlet,
w_outlet,
w_circle,
w_barwall,
]
)
n = ufl.FacetNormal(mesh)
return (
dx,
ds,
dS,
n,
bcs,
DVP,
d_,
v_,
p_,
w_,
dvp_,
k,
phi,
psi,
gamma,
beta,
time_constant,
)
Dear Federico,
I believe Dokken will get back to you soon, either here or on The FEniCS Discourse site.
Best, Kristian
On Tue, 21 Oct 2025 at 13:58, Federico Betti @.***> wrote:
federicobetti-ethz left a comment (KVSlab/turtleFSI#113) https://github.com/KVSlab/turtleFSI/issues/113#issuecomment-3426234282
Hi @jorgensd https://github.com/jorgensd, thank you for your reply. I attach below the two functions I wrote for translating in dolfinx the creation of the meshtags and bcs for the Turek Hron benchmark.
def setup_mesh_and_tags(R, H, L, f_L, f_H, c_x, c_y): """Setup mesh and create boundary/domain tags for FSI problem.
Args: R: Radius of circle H: Channel height L: Channel length f_L: Flag length f_H: Flag height c_x: Circle center x-coordinate c_y: Circle center y-coordinate Returns: mesh: The computational mesh domain_tag: Tags for fluid/solid domains facet_tag: Tags for boundaries """ filepath = path.join( path.dirname(path.dirname(path.abspath(__file__))), "mesh", "mesh.xdmf" ) with dolfinx.io.XDMFFile(MPI.COMM_WORLD, filepath, "r") as xdmf: mesh = xdmf.read_mesh(name="Grid") def inlet_boundary(x): return np.isclose(x[0], 0) def outlet_boundary(x): return np.isclose(x[0], L) def wall_boundary(x): return np.logical_or(np.isclose(x[1], H), np.isclose(x[1], 0)) def bar_boundary(x): return np.logical_or( np.logical_or( np.isclose(x[1], c_y + f_H / 2), np.isclose(x[1], c_y - f_H / 2) ), np.isclose(x[0], c_x + R + f_L), ) def circle_boundary(x): return (x[0] - c_x) ** 2 + (x[1] - c_y) ** 2 < R**2 + DOLFIN_EPS * 1e5 def barwall_boundary(x): belongs_to_circle = circle_boundary(x) belongs_to_bar = np.logical_and(x[1] <= c_y + f_H / 2, x[1] >= c_y - f_H / 2) x_larger_than_bar = x[0] > c_x return np.logical_and( belongs_to_circle, np.logical_and(belongs_to_bar, x_larger_than_bar) ) boundaries = [ (2, wall_boundary), (3, inlet_boundary), (4, outlet_boundary), (5, bar_boundary), (6, circle_boundary), (7, barwall_boundary), ] facet_indices, facet_markers = [], [] facetdim = mesh.topology.dim - 1 for marker, locator in boundaries: boundary_facets = dolfinx.mesh.locate_entities(mesh, facetdim, locator) facet_indices.append(boundary_facets) facet_markers.append(np.full_like(boundary_facets, marker)) facet_indices = np.hstack(facet_indices).astype(np.int32) facet_markers = np.hstack(facet_markers).astype(np.int32) sorted_facets = np.argsort(facet_indices) facet_tag = dolfinx.mesh.meshtags( mesh, facetdim, facet_indices[sorted_facets], facet_markers[sorted_facets] ) def solid_interior(x): flag_y_region = np.logical_and(x[1] <= c_y + f_H / 2, x[1] >= c_y - f_H / 2) flag_x_region = np.logical_and(x[0] >= c_x, x[0] <= c_x + R + f_L) flag_region = np.logical_and(flag_y_region, flag_x_region) return flag_region interiors = [ (2, solid_interior), ] domain_indices, domain_markers = [], [] vdim = mesh.topology.dim all_cells = np.arange(mesh.topology.index_map(vdim).size_local) for marker, locator in interiors: interior_cells = dolfinx.mesh.locate_entities(mesh, vdim, locator) remaining_cells = np.setdiff1d(all_cells, interior_cells) domain_indices.append(interior_cells) domain_markers.append(np.full_like(interior_cells, marker, dtype=np.int32)) domain_indices.append(remaining_cells) domain_markers.append(np.full_like(remaining_cells, 1, dtype=np.int32)) domain_indices = np.hstack(domain_indices).astype(np.int32) domain_markers = np.hstack(domain_markers).astype(np.int32) sorted_domains = np.argsort(domain_indices) domain_tag = dolfinx.mesh.meshtags( mesh, vdim, domain_indices[sorted_domains], domain_markers[sorted_domains] ) mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim) mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim - 1) mesh_dir = path.dirname(path.dirname(path.abspath(__file__))) facet_tag_file = path.join(mesh_dir, "mesh", "facet_tags.xdmf") with dolfinx.io.XDMFFile(MPI.COMM_WORLD, facet_tag_file, "w") as xdmf: xdmf.write_mesh(mesh) xdmf.write_meshtags(facet_tag) domain_tag_file = path.join(mesh_dir, "mesh", "domain_tags.xdmf") with dolfinx.io.XDMFFile(MPI.COMM_WORLD, domain_tag_file, "w") as xdmf: xdmf.write_mesh(mesh) xdmf.write_meshtags(domain_tag) return mesh, domain_tag, facet_tagdef setup_function_spaces(mesh, domain_tag, facet_tag, Um, H, dt): """Setup function spaces, measures, and boundary conditions.
Args: mesh: The computational mesh domain_tag: Tags for fluid/solid domains facet_tag: Tags for boundaries Um: Maximum inlet velocity H: Channel height dt: Time step Returns: Tuple containing: - dx, ds, dS: Measures for domain, boundary, interior facets - n: Facet normal - bcs: List of boundary conditions - DVP: Mixed function space - d_, v_, p_, w_: Field dictionaries - dvp_: Solution function dictionary - k: Time step constant - phi, psi, gamma, beta: Test functions """ ds = ufl.Measure("ds", domain=mesh, subdomain_data=facet_tag) dS = ufl.Measure("dS", domain=mesh, subdomain_data=facet_tag) dx = ufl.Measure("dx", domain=mesh, subdomain_data=domain_tag) de = ufl.VectorElement("CG", mesh.ufl_cell(), 2) ve = ufl.VectorElement("CG", mesh.ufl_cell(), 2) pe = ufl.FiniteElement("CG", mesh.ufl_cell(), 1) k = dolfinx.fem.Constant(mesh, c=dt) DVP_elem = ufl.MixedElement(de, ve, pe, de) DVP = dolfinx.fem.FunctionSpace(mesh, DVP_elem) de_space = dolfinx.fem.FunctionSpace(mesh, de) ve_space = dolfinx.fem.FunctionSpace(mesh, ve) pe_space = dolfinx.fem.FunctionSpace(mesh, pe) d_ = {} v_ = {} p_ = {} w_ = {} dvp_ = {} times = ["n-2", "n-1", "n"] for time_ in times: dvp = dolfinx.fem.Function(DVP) dvp_[time_] = dvp dvp_list = ufl.split(dvp) d_[time_] = dvp_list[0] v_[time_] = dvp_list[1] p_[time_] = dvp_list[2] w_[time_] = dvp_list[3] phi, psi, gamma, beta = ufl.TestFunctions(DVP) time_constant = dolfinx.fem.Constant(mesh, c=0.5) def inlet_profile(x): u_x = 1.5 * Um * x[1] * (H - x[1]) / (H / 2.0) ** 2 u_y = np.zeros_like(x[1]) return time_constant.value * np.stack([u_x, u_y]) inlet_profile_func = dolfinx.fem.Function(ve_space) inlet_profile_func.interpolate(inlet_profile) bcs = [] inlet_facets = facet_tag.find(3) inlet_dofs = dolfinx.fem.locate_dofs_topological( (DVP.sub(1), ve_space), mesh.topology.dim - 1, inlet_facets ) u_inlet = dolfinx.fem.dirichletbc(inlet_profile_func, inlet_dofs, DVP.sub(1)) wall_facets = facet_tag.find(2) wall_dofs = dolfinx.fem.locate_dofs_topological( (DVP.sub(1), ve_space), mesh.topology.dim - 1, wall_facets ) zero_v_func = dolfinx.fem.Function(ve_space) u_wall = dolfinx.fem.dirichletbc(zero_v_func, wall_dofs, DVP.sub(1)) circle_facets = facet_tag.find(6) circle_dofs = dolfinx.fem.locate_dofs_topological( (DVP.sub(1), ve_space), mesh.topology.dim - 1, circle_facets ) u_circle = dolfinx.fem.dirichletbc(zero_v_func, circle_dofs, DVP.sub(1)) barwall_facets = facet_tag.find(7) barwall_dofs = dolfinx.fem.locate_dofs_topological( (DVP.sub(1), ve_space), mesh.topology.dim - 1, barwall_facets ) u_barwall = dolfinx.fem.dirichletbc(zero_v_func, barwall_dofs, DVP.sub(1)) outlet_facets = facet_tag.find(4) outlet_dofs = dolfinx.fem.locate_dofs_topological( (DVP.sub(2), pe_space), mesh.topology.dim - 1, outlet_facets ) zero_p_func = dolfinx.fem.Function(pe_space) p_out = dolfinx.fem.dirichletbc(zero_p_func, outlet_dofs, DVP.sub(2)) bcs = [u_wall, u_inlet, u_circle, u_barwall, p_out] zero_d_func = dolfinx.fem.Function(de_space) wall_dofs = dolfinx.fem.locate_dofs_topological( (DVP.sub(0), de_space), mesh.topology.dim - 1, wall_facets ) d_wall = dolfinx.fem.dirichletbc(zero_d_func, wall_dofs, DVP.sub(0)) inlet_dofs = dolfinx.fem.locate_dofs_topological( (DVP.sub(0), de_space), mesh.topology.dim - 1, inlet_facets ) d_inlet = dolfinx.fem.dirichletbc(zero_d_func, inlet_dofs, DVP.sub(0)) outlet_dofs = dolfinx.fem.locate_dofs_topological( (DVP.sub(0), de_space), mesh.topology.dim - 1, outlet_facets ) d_outlet = dolfinx.fem.dirichletbc(zero_d_func, outlet_dofs, DVP.sub(0)) circle_dofs = dolfinx.fem.locate_dofs_topological( (DVP.sub(0), de_space), mesh.topology.dim - 1, circle_facets ) d_circle = dolfinx.fem.dirichletbc(zero_d_func, circle_dofs, DVP.sub(0)) barwall_dofs = dolfinx.fem.locate_dofs_topological( (DVP.sub(0), de_space), mesh.topology.dim - 1, barwall_facets ) d_barwall = dolfinx.fem.dirichletbc(zero_d_func, barwall_dofs, DVP.sub(0)) zero_w_func = dolfinx.fem.Function(de_space) wall_dofs = dolfinx.fem.locate_dofs_topological( (DVP.sub(3), de_space), mesh.topology.dim - 1, wall_facets ) w_wall = dolfinx.fem.dirichletbc(zero_w_func, wall_dofs, DVP.sub(3)) inlet_dofs = dolfinx.fem.locate_dofs_topological( (DVP.sub(3), de_space), mesh.topology.dim - 1, inlet_facets ) w_inlet = dolfinx.fem.dirichletbc(zero_w_func, inlet_dofs, DVP.sub(3)) outlet_dofs = dolfinx.fem.locate_dofs_topological( (DVP.sub(3), de_space), mesh.topology.dim - 1, outlet_facets ) w_outlet = dolfinx.fem.dirichletbc(zero_w_func, outlet_dofs, DVP.sub(3)) circle_dofs = dolfinx.fem.locate_dofs_topological( (DVP.sub(3), de_space), mesh.topology.dim - 1, circle_facets ) w_circle = dolfinx.fem.dirichletbc(zero_w_func, circle_dofs, DVP.sub(3)) barwall_dofs = dolfinx.fem.locate_dofs_topological( (DVP.sub(3), de_space), mesh.topology.dim - 1, barwall_facets ) w_barwall = dolfinx.fem.dirichletbc(zero_w_func, barwall_dofs, DVP.sub(3)) bcs.extend( [ d_wall, d_inlet, d_outlet, d_circle, d_barwall, w_wall, w_inlet, w_outlet, w_circle, w_barwall, ] ) n = ufl.FacetNormal(mesh) return ( dx, ds, dS, n, bcs, DVP, d_, v_, p_, w_, dvp_, k, phi, psi, gamma, beta, time_constant, )— Reply to this email directly, view it on GitHub https://github.com/KVSlab/turtleFSI/issues/113#issuecomment-3426234282, or unsubscribe https://github.com/notifications/unsubscribe-auth/AA4X54NKRFMAJHL6B5FG2HD3YYNVJAVCNFSM6AAAAACJPFMB5OVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZTIMRWGIZTIMRYGI . You are receiving this because you commented.Message ID: @.***>
Dear Kristian, thank you very much for your reply. With the code above, I get the following. I guess the boundary of the flag (with ID 5) is tagged incorrectly right now, right?
Dear Federico,
Yes, I agree unfortunately. That doesn't seem correct to me.
Best, Kristian
On Fri, 24 Oct 2025 at 13:15, Federico Betti @.***> wrote:
federicobetti-ethz left a comment (KVSlab/turtleFSI#113) https://github.com/KVSlab/turtleFSI/issues/113#issuecomment-3442576555
Dear Kristian, thank you very much for your reply. With the code above, I get the following. I guess the boundary of the flag is tagged incorrectly right now, right? image.png (view on web) https://github.com/user-attachments/assets/c4abf9e7-79fb-410e-a750-6857e0933992
— Reply to this email directly, view it on GitHub https://github.com/KVSlab/turtleFSI/issues/113#issuecomment-3442576555, or unsubscribe https://github.com/notifications/unsubscribe-auth/AA4X54J3ULKWNHOTMQJ7P7T3ZIC35AVCNFSM6AAAAACJPFMB5OVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZTINBSGU3TMNJVGU . You are receiving this because you commented.Message ID: @.***>
But the rest should be correct right? I am not sure why it happens, maybe I just need to add a small tolerance.
Could someone point me to the initial mesh file and tags that you are basing this on? The question above is really:
- How to read in a specific mesh, with existing cell and facet markers, to DOLFINx?
The cell and facet markers are not existing already, they are created. The mesh file is the one in the turtleFSI/mesh/TF_fsi.xml.gz file.
The picture above corresponds to the facet and cell tags that are created by my dolfinx code after importing the above mesh and using the code I sent a couple messages ago to create the meshtags. And there is weirdly enough an error in the boundary of the flag (yellow tag in the Paraview screenshot).
Could you give me a code that can reproduce the example. The code above is missing imports and run details.
Is the issue you are referring to the tagging of the little square here:
Hi @jorgensd , yes the tagging of the little square is the issue. Please find attached the MFE below.
I had initially converted the mesh file in src/mesh/TF_FSI.xml.gz (after extracting) to a xdmf that was compatible with dolfinx, with the command
meshio convert TF_FSI.xml mesh.xdmf
Upon doing this and having the file mesh.xdmf locally you should be able to reproduce the issue.
Let me know if you encounter any problem.
import dolfinx
from mpi4py import MPI
import numpy as np
import os
import ufl
DOLFIN_EPS = 1e-12
R = 0.05
H = 0.41
L = 2.5
f_L = 0.35
f_H = 0.02
c_x = 0.2
c_y = 0.2
filepath = os.path.relpath("mesh.xdmf")
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, filepath, "r") as xdmf:
mesh = xdmf.read_mesh(name="Grid")
def inlet_boundary(x):
return np.isclose(x[0], 0)
def outlet_boundary(x):
return np.isclose(x[0], L)
def wall_boundary(x):
return np.logical_or(np.isclose(x[1], H), np.isclose(x[1], 0))
def bar_boundary(x):
top_bottom = np.logical_and(
np.logical_or(np.isclose(x[1], c_y + f_H / 2), np.isclose(x[1], c_y - f_H / 2)),
np.logical_and(x[0] <= c_x + R + f_L, x[0] >= c_x)
)
right_edge = np.logical_and(
np.isclose(x[0], c_x + R + f_L),
np.logical_and(x[1] <= c_y + f_H / 2, x[1] >= c_y - f_H / 2)
)
return np.logical_or(top_bottom, right_edge)
def circle_boundary(x):
return (x[0] - c_x) ** 2 + (x[1] - c_y) ** 2 < R**2 + DOLFIN_EPS * 1e5
def barwall_boundary(x):
belongs_to_circle = circle_boundary(x)
belongs_to_bar = np.logical_and(x[1] <= c_y + f_H / 2, x[1] >= c_y - f_H / 2)
x_larger_than_bar = x[0] > c_x
return np.logical_and(
belongs_to_circle, np.logical_and(belongs_to_bar, x_larger_than_bar)
)
boundaries = [
(2, wall_boundary),
(3, inlet_boundary),
(4, outlet_boundary),
(5, bar_boundary),
(6, circle_boundary),
(7, barwall_boundary),
]
facet_indices, facet_markers = [], []
facetdim = mesh.topology.dim - 1
for marker, locator in boundaries:
boundary_facets = dolfinx.mesh.locate_entities(mesh, facetdim, locator)
facet_indices.append(boundary_facets)
facet_markers.append(np.full_like(boundary_facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = dolfinx.mesh.meshtags(
mesh, facetdim, facet_indices[sorted_facets], facet_markers[sorted_facets]
)
def solid_interior(x):
flag_y_region = np.logical_and(x[1] <= c_y + f_H / 2, x[1] >= c_y - f_H / 2)
flag_x_region = np.logical_and(x[0] >= c_x, x[0] <= c_x + R + f_L)
flag_region = np.logical_and(flag_y_region, flag_x_region)
return flag_region
interiors = [
(2, solid_interior),
]
domain_indices, domain_markers = [], []
vdim = mesh.topology.dim
all_cells = np.arange(mesh.topology.index_map(vdim).size_local)
for marker, locator in interiors:
interior_cells = dolfinx.mesh.locate_entities(mesh, vdim, locator)
remaining_cells = np.setdiff1d(all_cells, interior_cells)
domain_indices.append(interior_cells)
domain_markers.append(np.full_like(interior_cells, marker, dtype=np.int32))
domain_indices.append(remaining_cells)
domain_markers.append(np.full_like(remaining_cells, 1, dtype=np.int32))
domain_indices = np.hstack(domain_indices).astype(np.int32)
domain_markers = np.hstack(domain_markers).astype(np.int32)
sorted_domains = np.argsort(domain_indices)
domain_tag = dolfinx.mesh.meshtags(
mesh, vdim, domain_indices[sorted_domains], domain_markers[sorted_domains]
)
mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim - 1)
facet_tag_file = os.path.relpath("facet_tags.xdmf")
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, facet_tag_file, "w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_meshtags(facet_tag)
domain_tag_file = os.path.relpath("domain_tags.xdmf")
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, domain_tag_file, "w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_meshtags(domain_tag)
This facet is marked because both it's vertices satisfies the constraint. You can observe this if you look at the mesh
If you change the logic to also check midpoints, you get:
mesh.topology.create_connectivity(facetdim, mesh.topology.dim)
for marker, locator in boundaries:
boundary_facets = dolfinx.mesh.locate_entities(mesh, facetdim, locator)
midpoints = dolfinx.mesh.compute_midpoints(mesh, facetdim, boundary_facets)
boundary_facets = boundary_facets[locator(midpoints.T)]
facet_indices.append(boundary_facets)
facet_markers.append(np.full_like(boundary_facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
def circle_boundary(x, tol=1e-10): return (x[0] - c_x) ** 2 + (x[1] - c_y) ** 2 < R**2 + tol
def barwall_boundary(x, tol=1e-10): belongs_to_circle = circle_boundary(x, tol) belongs_to_bar = np.logical_and( x[1] <= c_y + f_H / 2 + tol, x[1] >= c_y - f_H / 2 - tol ) x_larger_than_bar = x[0] > c_x return np.logical_and( belongs_to_circle, np.logical_and(belongs_to_bar, x_larger_than_bar) )
boundaries = [ (2, wall_boundary), (3, inlet_boundary), (4, outlet_boundary), (5, bar_boundary), (6, circle_boundary), (7, barwall_boundary), ]
facet_indices, facet_markers = [], [] facetdim = mesh.topology.dim - 1 mesh.topology.create_connectivity(facetdim, mesh.topology.dim) num_facets_local = ( mesh.topology.index_map(facetdim).size_local + mesh.topology.index_map(facetdim).num_ghosts ) facet_markers = np.full(num_facets_local, -1, dtype=np.int32) for marker, locator in boundaries: boundary_facets = dolfinx.mesh.locate_entities(mesh, facetdim, locator) midpoints = dolfinx.mesh.compute_midpoints(mesh, facetdim, boundary_facets) boundary_facets = boundary_facets[locator(midpoints.T)] facet_markers[boundary_facets] = marker
is_marked = facet_markers >= 0 facets = np.flatnonzero(is_marked) values = facet_markers[facets] facet_tag = dolfinx.mesh.meshtags(mesh, facetdim, facets, values)
Hi @jorgensd , thank you so much for the input. The facet tag seems now to be working fine. Thank you again for all of your help.
@KristianValen-Sendstad in the past few weeks I was able to work a bit on this and I now have a working dolfinx solver both for the CFD and FSI cases.
The only thing left is that I am having some trouble reproducing the results of the fsi2 case in Turek-Hron (so with U_bar=1 and rho_s=1e4), both with my code and with yours. Is it something you experienced as well?
On the other hand, I get results converging to the reference values in the publication for all the other five cases (cfd1,2,3 and fsi1, fsi3).