FEniCSx transition

[federicobetti-ethz]

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

[KristianValen-Sendstad]

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) <#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 <#113>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AA4X54NPSRJB5ZZW5DH67FL3YDPIDAVCNFSM6AAAAACJPFMB5OVHI2DSMVQWIX3LMV43ASLTON2WKOZTGUZDMMBQGM2DCNQ> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[federicobetti-ethz]

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 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.

[jorgensd]

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 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.

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.

[federicobetti-ethz]

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,
    )

[KristianValen-Sendstad]

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) <#113 (comment)> 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_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, ) — Reply to this email directly, view it on GitHub <#113 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AA4X54NKRFMAJHL6B5FG2HD3YYNVJAVCNFSM6AAAAACJPFMB5OVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZTIMRWGIZTIMRYGI> . You are receiving this because you commented.Message ID: ***@***.***>

[federicobetti-ethz]

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?

[KristianValen-Sendstad]

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) <#113 (comment)> 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 <#113 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AA4X54J3ULKWNHOTMQJ7P7T3ZIC35AVCNFSM6AAAAACJPFMB5OVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZTINBSGU3TMNJVGU> . You are receiving this because you commented.Message ID: ***@***.***>

[federicobetti-ethz]

But the rest should be correct right? I am not sure why it happens, maybe I just need to add a small tolerance.

[jorgensd]

Could someone point me to the initial mesh file and tags that you are basing this on?
The question above is really:

1. How to read in a specific mesh, with existing cell and facet markers, to DOLFINx?

[federicobetti-ethz]

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).

[jorgensd]

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:

?

[federicobetti-ethz]

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)

[jorgensd]

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)

I would do: ```python 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, 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)

[federicobetti-ethz]

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.