Refine troublesome cell

create_cube_dolfinx.py

from mpi4py import MPI

import dolfinx
import numpy as np
from pathlib import Path

cube = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, nx=3, ny=3, nz=3)


def left_cells(x):
    return (x[0] < 0.5) & (x[1] < 0.5)


num_cells = cube.topology.index_map(cube.topology.dim).size_local
cells = dolfinx.mesh.locate_entities(cube, cube.topology.dim, left_cells)
all_cells = np.arange(num_cells, dtype=np.int32)
cell_marker = np.full(num_cells, 1, dtype=np.int32)
cell_marker[cells] = 2
cell_tags = dolfinx.mesh.meshtags(cube, cube.topology.dim, all_cells, cell_marker)
cell_tags.name = "cell_tags"

cube.topology.create_entities(cube.topology.dim - 1)
num_facets = cube.topology.index_map(cube.topology.dim - 1).size_local
all_facets = np.arange(num_facets, dtype=np.int32)
facet_marker = np.full(num_facets, 1, dtype=np.int32)
facets = dolfinx.mesh.locate_entities(cube, cube.topology.dim - 1, left_cells)
facet_marker[facets] = 2
facet_tags = dolfinx.mesh.meshtags(
    cube, cube.topology.dim - 1, all_facets, facet_marker
)
facet_tags.name = "facet_tags"
outfile = Path("mesh.xdmf")
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, outfile, "w") as file:
    file.write_mesh(cube)
    file.write_meshtags(cell_tags, cube.geometry)


outfile_facets = Path("mesh_facets.xdmf")
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, outfile_facets, "w") as file:
    file.write_mesh(cube)
    cube.topology.create_connectivity(cube.topology.dim - 1, cube.topology.dim)
    file.write_meshtags(facet_tags, cube.geometry)

identify_cell_dolfinx.py

# Mesh refinement of "overconstrained cells in the Stokes setting
#
# SPDX-License-Identifier: MIT
# Authors: Jørgen S. Dokken and Marius Causemann
import ufl
import dolfinx
import numpy as np
import meshio
import warnings
from tqdm import tqdm


def create_mesh(ct, points, cells, cell_data):
    out_mesh = meshio.Mesh(
        points=points, cells={ct: cells}, cell_data={"name_to_read": [cell_data]}
    )
    return out_mesh


def create_new_cells(current_cell, needed_vertex, new_vertex):
    """
    Given a cell `[v0, v1, v2, v3]`, a required vertex `v2`, and a new vertex `v`,
    `[v1,v2,v3,v], [v0,v2,v3,v], [v0,v1,v2,v]`
    """

    include = np.argwhere(current_cell == needed_vertex)
    new_cells = np.full((3, 4), -1, dtype=np.int64)
    assert len(include) == 1
    mask = np.full(4, False, dtype=np.bool_)
    loop = np.arange(4)
    loop = np.delete(loop, include)
    for i, insert_pos in enumerate(loop):
        mask[:] = True
        mask[insert_pos] = False
        new_cells[i] = np.hstack([current_cell[mask], [new_vertex]])
    return new_cells


def create_facet_marker_structure(facet, new_vertex):
    """Split facet of current cell in 3 and create marker structure

    Args:
        new_vertex: w
        facet: [v0,v2,v3]
    Returns [v0,v2, w], [v0,w,v3], [v2,v3,w]
    """
    mask = np.full(3, True, dtype=np.bool_)
    new_facets = np.full((3, 3), -1, dtype=np.int64)
    for i in range(3):
        mask[:] = True
        mask[i] = False
        new_facets[i] = np.hstack([facet[mask], [new_vertex]])
    return new_facets


def fix_overconstrained_cells(
    mesh: dolfinx.mesh.Mesh,
    cell_function: dolfinx.mesh.MeshTags,
    facet_function: dolfinx.mesh.MeshTags,
) -> tuple[int, dolfinx.mesh.Mesh, dolfinx.mesh.MeshTags, dolfinx.mesh.MeshTags]:
    """Refine cells that are "overconstrained", in the sense that all vertices are on the boundary.
    Additionally, remove all cells that only has exterior facets, i.e. is connected either by edges or vertices.

    Note:
        This function does not check that the resulting mesh is overconstrained and should be called iteratively.

    Note:
        This only works for 3D meshes and in serial
    Args:
        mesh: Input mesh
        cell_function: Cell markers to transfer
        facet_function: Facet markers to transfer

    Returns:
        Number of new cells, the new mesh and corresponding cell and facet markers
    """
    assert mesh.comm.size == 1

    # Create all required topological data
    tdim = mesh.topology.dim
    num_vertices = mesh.topology.index_map(0).size_local
    mesh.topology.create_connectivity(tdim, tdim - 1)
    c_to_f = mesh.topology.connectivity(tdim, tdim - 1)

    mesh.topology.create_connectivity(tdim - 1, tdim)
    f_to_c = mesh.topology.connectivity(tdim - 1, tdim)

    mesh.topology.create_connectivity(tdim - 1, 0)
    f_to_v = mesh.topology.connectivity(tdim - 1, 0)
    c_to_v = mesh.topology.connectivity(tdim, 0)

    cell_map = mesh.topology.index_map(tdim)
    num_cells = cell_map.size_local

    facet_map = mesh.topology.index_map(tdim - 1)
    num_facets = facet_map.size_local
    facet_midpoints = dolfinx.mesh.compute_midpoints(
        mesh, mesh.topology.dim - 1, np.arange(num_facets, dtype=np.int32)
    )
    new_vertex_coordinates = []
    new_cells = []

    removed_cells = []
    new_vertex_counter = 0

    mask = np.full(4, True, dtype=np.bool_)

    new_cell_marker_array = []
    new_facet_marker_array = []
    new_marked_facets = []
    removed_facets = []

    removed_cells = []
    for i in tqdm(range(num_cells)):
        if i in removed_cells:
            continue
        facets = c_to_f.links(i)
        exterior_facets = np.array(
            [facet for facet in facets if len(f_to_c.links(facet)) == 1], dtype=np.int32
        )
        interior_facets = [facet for facet in facets if len(f_to_c.links(facet)) == 2]
        if not len(interior_facets) > 0:
            assert len(exterior_facets) == 4
            print(f"Cell {i} with midpoint {dolfinx.mesh.compute_midpoints(mesh, mesh.topology.dim, np.array([i],dtype=np.int32))} has no interior facets")
            removed_cells.append(i)
            [removed_facets.append(f) for f in exterior_facets]
            continue
        exterior_vertices = dolfinx.mesh.compute_incident_entities(
            mesh.topology, exterior_facets, mesh.topology.dim - 1, 0
        )
        assert len(np.unique(exterior_vertices)) == len(exterior_vertices)
        if len(exterior_vertices) == 4:
            all_cells = f_to_c.links(interior_facets[0])
            other_cell = np.setdiff1d(all_cells, [i])[0]
            if other_cell in removed_cells:
                warnings.warn("Cell already removed, should call this function again")
                continue
            current_vertices = c_to_v.links(i)
            interior_facet_vertices = f_to_v.links(interior_facets[0])

            # Get position of new vertex on midpoint of facet
            coord = facet_midpoints[interior_facets[0]]
            new_vertex_coordinates.append(coord)
            all_needs = np.setdiff1d(current_vertices, interior_facet_vertices)

           
            # Transfer facet marker if and only if it was in original tag
            pos = np.flatnonzero(facet_function.indices==interior_facets[0])
            if len(pos)>0:
                # Get all new sub-facets
                new_marked_facets.append(
                    create_facet_marker_structure(
                        interior_facet_vertices, num_vertices + new_vertex_counter
                    )
                )
                for _ in range(3):
                    # Check if facet has been marked in previous grid
                        assert len(pos) == 1
                        new_facet_marker_array.append(facet_function.values[pos])
                removed_facets.append(interior_facets[0])

            # Split troublesome cell in 3
            assert len(all_needs) == 1
            include = np.argwhere(current_vertices == all_needs[0])
            assert len(include) == 1
            new_cells.append(
                create_new_cells(
                    current_vertices, all_needs[0], num_vertices + new_vertex_counter
                )
            )

            other_cell_connectivity = c_to_v.links(other_cell)
            other_needs = np.setdiff1d(other_cell_connectivity, interior_facet_vertices)
            new_cells.append(
                create_new_cells(
                    other_cell_connectivity,
                    other_needs[0],
                    num_vertices + new_vertex_counter,
                )
            )

            for _ in range(3):
                new_cell_marker_array.append(cell_function.values[i])
            for _ in range(3):
                new_cell_marker_array.append(cell_function.values[other_cell])

            new_vertex_counter += 1
            removed_cells.append(i)
            removed_cells.append(other_cell)
    new_marked_facets = np.array(new_marked_facets, dtype=np.int64).reshape(-1, 3)
    new_facet_marker_array = np.hstack(new_facet_marker_array).astype(np.int32)
    print("checked all cells!")
        
    if len(new_cells) == 0:
        new_cells_as_array = np.zeros((0, 4), dtype=np.int64)
    else:
        new_cells_as_array = np.vstack(new_cells)
    assert len(new_cells_as_array) == len(new_cell_marker_array)

    if len(removed_cells) == 0:
        removed_cells = np.array([], dtype=np.int64)
    else:
        removed_cells = np.unique(np.hstack(removed_cells).astype(np.int64))

    # Gather all cells
    remaining_cells = np.arange(num_cells, dtype=np.int32)
    remaining_cells = np.delete(remaining_cells, removed_cells)

    all_cells = dolfinx.mesh.entities_to_geometry(
        mesh, mesh.topology.dim, remaining_cells
    )
    all_new_cells = np.vstack([all_cells, new_cells_as_array]).astype(np.int64)

    new_vertex_numbering = np.unique(all_new_cells.flatten())
    all_to_reduced_num_vertices = np.full(num_vertices + new_vertex_counter, -1, dtype=np.int64)
    all_to_reduced_num_vertices[new_vertex_numbering] = np.arange(len(new_vertex_numbering))
    all_new_cells = all_to_reduced_num_vertices[all_new_cells]


    # Gather all coordinates
    all_coords = np.zeros((num_vertices + new_vertex_counter, 3), dtype=np.float64)
    if new_vertex_counter > 0:
        all_new_vertex_coordinates = np.vstack(new_vertex_coordinates)
    else:
        all_new_vertex_coordinates = np.zeros((0, 3), dtype=np.float64)
    all_coords[:num_vertices, :] = mesh.geometry.x
    all_coords[num_vertices:, :] = all_new_vertex_coordinates

    all_coords = all_coords[new_vertex_numbering]

    # Transfer markers
    assert len(cell_function.indices) == num_cells
    mask = np.full(num_cells, True, dtype=np.bool_)
    mask[removed_cells] = False
    new_values = cell_function.values[mask]
    all_values = np.hstack([new_values, new_cell_marker_array])
    ufl_domain = ufl.Mesh(mesh._ufl_domain.ufl_coordinate_element())  # type: ignore
    new_mesh = dolfinx.mesh.create_mesh(
        mesh.comm, all_new_cells, all_coords, ufl_domain
    )
    # Create MeshTags for cells
    print("Transferring cell markers")
    local_entities, local_values = dolfinx.io.distribute_entity_data(
        new_mesh, new_mesh.topology.dim, all_new_cells, all_values
    )
    new_mesh.topology.create_connectivity(mesh.topology.dim, 0)
    adj = dolfinx.graph.adjacencylist(local_entities)
    ct = dolfinx.mesh.meshtags_from_entities(
        new_mesh, new_mesh.topology.dim, adj, local_values.astype(np.int32, copy=False)
    )
    ct.name = cell_function.name

    # Create facet marker structure
    print("Transferring facet marker")
    facets_to_copy = facet_function.indices.copy()
    facets_to_keep = np.invert(np.isin(facets_to_copy, removed_facets))
    assert np.allclose(f_to_v.offsets[1:]-f_to_v.offsets[:-1], 3)
    conn_arr = f_to_v.array.reshape(-1, 3)
    new_facet_array = conn_arr[facets_to_copy[facets_to_keep], :]
    new_facet_values_array = facet_function.values[facets_to_keep]
    new_facet_values_array = np.hstack(new_facet_values_array, dtype=np.int32).reshape(
        -1
    )
    # Renumber the vertices of facets to align with the reduced set of vertices
    facet_connectivity = np.vstack([new_facet_array, new_marked_facets])
    facet_connectivity = all_to_reduced_num_vertices[facet_connectivity].astype(np.int64)
    assert(facet_connectivity != -1).all()
    facet_values = np.hstack([new_facet_values_array, new_facet_marker_array])

    local_entities, local_values = dolfinx.io.distribute_entity_data(
        new_mesh, new_mesh.topology.dim - 1, facet_connectivity, facet_values
    )
    new_mesh.topology.create_connectivity(mesh.topology.dim, 0)
    adj = dolfinx.graph.adjacencylist(local_entities)
    new_mesh.topology.create_connectivity(
        new_mesh.topology.dim - 1, new_mesh.topology.dim
    )
    ft = dolfinx.mesh.meshtags_from_entities(
        new_mesh,
        new_mesh.topology.dim - 1,
        adj,
        local_values.astype(np.int32, copy=False),
    )
    ft.name = facet_function.name

    return (
        new_cells_as_array.shape[0],
        new_mesh,
        ct,
        ft,
    )


import argparse
from pathlib import Path

from mpi4py import MPI

if __name__ == "__main__":
    # Create the parser
    parser = argparse.ArgumentParser(
        description="Refine cells that are overconstrained, in the sense that all vertices are on the boundary."
    )

    # Add arguments
    parser.add_argument("infile", type=Path, help="The input file to be processed")
    parser.add_argument(
        "outfile", type=Path, help="The output file to write results to"
    )
    parser.add_argument("--cell_tag", type=str, default="cell_tags", help="Name of cell tag variable in file")
    parser.add_argument("--facet_tag", type=str, default="facet_tags", help="Name of facet tag variable in file")

    # Parse the arguments
    args = parser.parse_args()

    print(f"reading in mesh: {args.infile}")
    with dolfinx.io.XDMFFile(MPI.COMM_WORLD, args.infile, "r") as xdmf:
        mesh = xdmf.read_mesh()
        cell_tags = xdmf.read_meshtags(mesh, name=args.cell_tag)
        mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
        try:
            facet_tags = xdmf.read_meshtags(mesh, name=args.facet_tag)
        except RuntimeError:
            facet_file_path = str(Path(args.infile).with_suffix("")) + "_facets.xdmf"
            print(f"reading in facet file: {facet_file_path}")
            with dolfinx.io.XDMFFile(MPI.COMM_WORLD, facet_file_path, "r") as facet_infile:
                mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
                facet_tags = facet_infile.read_meshtags(mesh, name=args.facet_tag)

    print(
        f"number of cells before refinement: {mesh.topology.index_map(mesh.topology.dim).size_global}"
    )
    (
        num_new_cells,
        new_mesh,
        new_cells,
        new_facets,
    ) = fix_overconstrained_cells(mesh, cell_tags, facet_tags)
    num_cells = new_mesh.topology.index_map(new_mesh.topology.dim).size_global
    print(f"number of cells after refinement: {num_cells}")
    print(f"writing to file: {args.outfile}")
    with dolfinx.io.XDMFFile(new_mesh.comm, args.outfile, "w") as xdmf:
        xdmf.write_mesh(new_mesh)
        xdmf.write_meshtags(new_cells, new_mesh.geometry)
    facet_outfile = (args.outfile.parent / (args.outfile.stem + "_facets")).with_suffix(
        ".xdmf"
    )
    with dolfinx.io.XDMFFile(new_mesh.comm, facet_outfile, "w") as xdmf:
        xdmf.write_mesh(new_mesh)
        xdmf.write_meshtags(new_facets, new_mesh.geometry)

identify_troublesome_cell.py

# Mesh refinement of "overconstrained cells in the Stokes setting
#
# SPDX-License-Identifier: MIT
# Authors: Jørgen S. Dokken and Marius Causemann

import dolfin
import numpy as np
import meshio
import warnings
from tqdm import tqdm

def create_mesh(ct, points, cells, cell_data):
    out_mesh = meshio.Mesh(
        points=points, cells={ct: cells}, cell_data={"name_to_read": [cell_data]}
    )
    return out_mesh


def create_new_cells(current_cell, needed_vertex, new_vertex):
    """
    Given a cell `[v0, v1, v2, v3]`, a required vertex `v2`, and a new vertex `v`,
    `[v1,v2,v3,v], [v0,v2,v3,v], [v0,v1,v2,v]`
    """

    include = np.argwhere(current_cell == needed_vertex)
    new_cells = np.full((3, 4), -1, dtype=np.int64)
    assert len(include) == 1
    mask = np.full(4, False, dtype=np.bool_)
    loop = np.arange(4)
    loop = np.delete(loop, include)
    for i, insert_pos in enumerate(loop):
        mask[:] = True
        mask[insert_pos] = False
        new_cells[i] = np.hstack([current_cell[mask], [new_vertex]])
    return new_cells


def create_facet_marker_structure(facet, new_vertex):
    """Split facet of current cell in 3 and create marker structure

    Args:
        new_vertex: w
        facet: [v0,v2,v3]
    Returns [v0,v2, w], [v0,w,v3], [v2,v3,w]
    """
    mask = np.full(3, True, dtype=np.bool_)
    new_facets = np.full((3, 3), -1, dtype=np.int64)
    for i in range(3):
        mask[:] = True
        mask[i] = False
        new_facets[i] = np.hstack([facet[mask], [new_vertex]])
    return new_facets


def fix_overconstrained_cells(
    mesh: dolfin.Mesh,
    cell_function: dolfin.MeshFunction,
    facet_function: dolfin.MeshFunction,
) -> tuple[int, meshio.Mesh, meshio.Mesh]:
    """Refine cells that are "overconstrained", in the sense that all vertices are on the boundary.

    Note:
        This function does not check that the resulting mesh is overconstrained and should be called iteratively.

    Note:
        This only works for 3D meshes and in serial
    Args:
        mesh: Input mesh
        cell_function: Cell markers to transfer
        facet_function: Facet markers to transfer

    Returns:
        _description_
    """

    # Create all required topological data
    mesh.init(mesh.topology().dim(), mesh.topology().dim() - 1)
    c_to_f = mesh.topology()(mesh.topology().dim(), mesh.topology().dim() - 1)

    mesh.init(mesh.topology().dim() - 1, mesh.topology().dim())
    f_to_c = mesh.topology()(mesh.topology().dim() - 1, mesh.topology().dim())

    assert dolfin.MPI.comm_world.size == 1

    mesh.init(mesh.topology().dim() - 1, 0)
    f_to_v = mesh.topology()(mesh.topology().dim() - 1, 0)

    c_to_v = mesh.topology()(mesh.topology().dim(), 0)

    num_vertices = mesh.num_vertices()
    new_vertex_coordinates = []
    new_cells = []

    removed_cells = []
    new_vertex_counter = 0

    mask = np.full(4, True, dtype=np.bool_)

    new_cell_marker_array = []
    new_facet_marker_array = []
    new_marked_facets = []
    removed_facets = []
    print("checking cells...")
    for cell in tqdm(dolfin.cells(mesh), total=mesh.num_cells()):
        if cell.index() in removed_cells:
            continue
        facets = c_to_f(cell.index())
        external_facets = []
        interior_facets = []
        [external_facets.append(facet) for facet in facets if len(f_to_c(facet)) == 1]
        [interior_facets.append(facet) for facet in facets if len(f_to_c(facet)) == 2]
        exterior_vertices = []
        for facet in external_facets:
            exterior_vertices.extend(f_to_v(facet))
        unique_exterior_vertices = np.unique(exterior_vertices)
        if len(unique_exterior_vertices) == 4:
            # Split cell on other side of facet in 3
            all_cells = f_to_c(interior_facets[0])
            other_cell = np.setdiff1d(all_cells, [cell.index()])[0]
            if other_cell in removed_cells:
                warnings.warn("Cell already removed, should call this function again")
                continue
            current_vertices = c_to_v(cell.index())
            interior_facet_vertices = f_to_v(interior_facets[0])
            # Create new vertex on midpoint of acet
            coord = dolfin.Facet(mesh, interior_facets[0]).midpoint().array()
            new_vertex_coordinates.append(coord)

            all_needs = np.setdiff1d(current_vertices, interior_facet_vertices)

            # Get all new sub-facets
            new_marked_facets.append(
                create_facet_marker_structure(
                    interior_facet_vertices, num_vertices + new_vertex_counter
                )
            )

            for _ in range(3):
                new_facet_marker_array.append(
                    facet_function.array()[interior_facets[0]]
                )
            removed_facets.append(interior_facets[0])

            # Split troublesome cell in 3
            assert len(all_needs) == 1
            include = np.argwhere(current_vertices == all_needs[0])
            assert len(include) == 1
            new_cells.append(
                create_new_cells(
                    current_vertices, all_needs[0], num_vertices + new_vertex_counter
                )
            )

            other_cell_connectivity = c_to_v(other_cell)
            other_needs = np.setdiff1d(other_cell_connectivity, interior_facet_vertices)
            new_cells.append(
                create_new_cells(
                    other_cell_connectivity,
                    other_needs[0],
                    num_vertices + new_vertex_counter,
                )
            )

            for _ in range(3):
                new_cell_marker_array.append(cell_function.array()[cell.index()])
            for _ in range(3):
                new_cell_marker_array.append(cell_function.array()[other_cell])

            new_vertex_counter += 1
            removed_cells.append(cell.index())
            removed_cells.append(other_cell)
    print("checked all cells!")

    if len(new_cells) == 0:
        new_cells_as_array = np.zeros((0, 4), dtype=np.int64)
    else:
        new_cells_as_array = np.vstack(new_cells)
    assert len(new_cells_as_array) == len(new_cell_marker_array)

    if len(removed_cells) == 0:
        removed_cells = np.array([], dtype=np.int64)
    else:
        removed_cells = np.unique(np.hstack(removed_cells).astype(np.int64))

    # Gather all cells
    all_cells = np.zeros((mesh.num_cells() - len(removed_cells), 4), dtype=np.int64)
    c = 0
    for cell in dolfin.cells(mesh):
        if cell.index() not in removed_cells:
            all_cells[c, :] = c_to_v(cell.index())
            c += 1

    all_new_cells = np.vstack([all_cells, new_cells_as_array])

    # Gather all coordinates
    all_coords = np.zeros(
        (mesh.num_vertices() + new_vertex_counter, 3), dtype=np.float64
    )
    if new_vertex_counter > 0:
        all_new_vertex_coordinates = np.vstack(new_vertex_coordinates)
    else:
        all_new_vertex_coordinates = np.zeros((0, 3), dtype=np.float64)
    for vertex in dolfin.vertices(mesh):
        all_coords[vertex.index(), :] = vertex.point().array()
    all_coords[mesh.num_vertices() :] = all_new_vertex_coordinates

    # Transfer markers

    mask = np.full(mesh.num_cells(), True, dtype=np.bool_)
    mask[removed_cells] = False
    new_values = cell_function.array()[mask]
    all_values = np.hstack([new_values, new_cell_marker_array])

    # New mesh
    new_mesh = create_mesh("tetra", all_coords, all_new_cells, all_values)

    # Create facet marker structure
    facet_connectivity_unrolled = []
    new_ff_array = []
    print("transfering facet marker")
    for facet in tqdm(dolfin.facets(mesh), total=mesh.num_facets()):
        if facet.index() not in removed_facets:
            if facet.index() >= mesh.num_facets():
                breakpoint()
            facet_connectivity_unrolled.append(f_to_v(facet.index()))
            new_ff_array.append(facet_function.array()[facet.index()])
    nf = np.array(new_marked_facets).reshape(-1, 3)

    facet_connectivity_unrolled = np.vstack([facet_connectivity_unrolled, nf])
    new_ff_array = np.hstack([new_ff_array, new_facet_marker_array])
    facet_mesh = create_mesh(
        "triangle", all_coords, facet_connectivity_unrolled, new_ff_array
    )

    return (
        new_cells_as_array.shape[0],
        new_mesh,
        facet_mesh,
    )




import argparse
from pathlib import Path


if __name__ == "__main__":
    # Create the parser
    parser = argparse.ArgumentParser(description="Refine cells that are overconstrained, in the sense that all vertices are on the boundary.")
    
    # Add arguments
    parser.add_argument('infile', type=str, help='The input file to be processed')
    parser.add_argument('outfile', type=str, help='The output file to write results to')
    
    # Parse the arguments
    args = parser.parse_args()

    mesh = dolfin.Mesh()

    print(f"reading in mesh: {args.infile}")
    with dolfin.XDMFFile(args.infile) as f:
        f.read(mesh)
        gdim = mesh.geometric_dimension()
        sm = dolfin.MeshFunction("size_t", mesh, gdim)
        f.read(sm)
    
    facet_file_path = str(Path(args.infile).with_suffix(""))+ "_facets.xdmf"
    print(f"reading in facet file: {facet_file_path}")
    mvc = dolfin.MeshValueCollection("size_t", mesh, gdim - 1)
    with dolfin.XDMFFile(facet_file_path) as facet_infile:
        facet_infile.read(mvc)

    print(f"number of cells before refinement: {mesh.num_cells()}")
    bm = dolfin.cpp.mesh.MeshFunctionSizet(mesh, mvc)
    (
        num_new_cells,
        new_mesh,
        new_facets,
    ) = fix_overconstrained_cells(mesh, sm, bm)
    print(f"number of cells after refinement: {new_mesh.cells[0].data.shape[0]}")
    print(f"writing to file: {args.outfile}")
    meshio.write(args.outfile, new_mesh)
    meshio.write(str(Path(args.outfile).with_suffix(""))+ "_facets.xdmf", new_facets)