Advice for meshing an internal boundary

1

I have an internal boundary that separates two domains, very similar to this: ArityMismatch | Acoustic-Elastic interaction problem - #7 by Thomas2

My code works great with a structured mesh.

When I switch from a structured mesh to an unstructured, there aren’t any vertices “near” (e.g., near(x,x0)) the interface between the two domains, and so I get this error:

Warning: Found no facets matching domain for boundary condition.

I tried increasing the resolution, but this is not going to be the solution: While some vertices become close enough to be found, not all do and the result is not well behaved. Apparently an unfeasible amount of refinement would be required.

I’m next thinking to manually find coordinates near the interface, move them to the interface, and then refine the mesh. I’m not sure if this will work, and so I thought to ask for advice in this venue.

FWIW, I am using the old mshr to make my structured mesh. I know it isn’t maintained any more but has generally worked very well for me.

2

Mshr has some limited functionality for boundary marking. However, I would strongly suggest moving to Gmsh, as you can either use the GUI or the Python interface to get the boundaries without having to use geometrical conditions.

3

After a few hours of tinkering, I have installed pygmsh and made an MWE. The MWE below creates a mesh that looks fine in Paraview (take that statement with a grain of salt as I didn’t examine it in any sophisticated way). But when I run the function create_mesh here (also reproduced on this form in different topics) I get the error:

    line_mesh = create_mesh(msh, "line", prune_z=True)
  File "convert-my-mesh.py", line 10, in create_mesh
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
  File "/home/fenics/.local/lib/python3.6/site-packages/meshio/_mesh.py", line 228, in get_cell_data
    [d for c, d in zip(self.cells, self.cell_data[name]) if c.type == cell_type]
ValueError: need at least one array to concatenate

MWE for the mesh generation:

import pygmsh

with pygmsh.occ.Geometry() as geom:
    geom.characteristic_length_max = 0.1
    rectangle1 = geom.add_rectangle([0.0, 0.0, 0.0], 1.0, 1.0)
    rectangle2 = geom.add_rectangle([1.0, -1.0, 0.0], 1.0, 2.00)
    domain=geom.boolean_fragments(rectangle1,rectangle2)
    for i,entity in enumerate(domain):
                  geom.add_physical(entity,f'{i}')
    mesh = geom.generate_mesh()
4

This is because:

does not tag surfaces (i.e. lines), and therefore

 line_mesh = create_mesh(mesh_from_file, "line", prune_z=True)

will not work.

As you can observe, the following command runs with no error:

import meshio
import numpy
import pygmsh
import gmsh


def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    points = mesh.points[:, :2] if prune_z else mesh.points
    out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={
                           "name_to_read": [cell_data]})
    return out_mesh


with pygmsh.occ.Geometry() as geom:
    geom.characteristic_length_max = 0.1
    rectangle1 = geom.add_rectangle([0.0, 0.0, 0.0], 1.0, 1.0)
    rectangle2 = geom.add_rectangle([1.0, -1.0, 0.0], 1.0, 2.00)
    domain = geom.boolean_fragments(rectangle1, rectangle2)
    for i, entity in enumerate(domain):
        print(i, entity)
        geom.add_physical(entity, f'{i}')
    mesh = geom.generate_mesh()

    gmsh.write("mesh.msh")

mesh_from_file = meshio.read("mesh.msh")
triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True)
meshio.write("mesh.xdmf", triangle_mesh)
5

Thank you so much for your help. The code that you provided works as desired. I’m now able to read the mesh and run my original fenics code with the new mesh. It isn’t quite working yet, however, because I still have my original issue about locating the interface between the two rectangles and apply a boundary condition there. I attempted the following.

I attempt to label the interface by adding two lines to the with-block:

with pygmsh.occ.Geometry() as geom:
    ...
    interface = geom.add_surface_loop([domain[1]])
    geom.add_physical(interface,'3')
    ...

However, I don’t think this is correct because the following code does not appear the correctly identify domains/boundaries (solution is zero when it shouldn’t be).

'''
Read mesh and mesh function
'''
mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)

mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)


'''
Mark domains and boundaries
'''
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1,0)
domains = MeshFunction("size_t", mesh, mesh.topology().dim())
dS = Measure('dS', domain=mesh, subdomain_data=mf)
dx = Measure('dx', domain=mesh, subdomain_data=mf)
ds = Measure('ds', domain=mesh, subdomain_data=mf)

I’ve tried incorporating information from this thread: How to extract and use physical labels from GMSH surfaces on FEniCS (2D mesh) - #5 by dokken but I can’t get things quite right.

6

For posterity, I’ll answer my own question. Based on another one of @dokken’s responses somewhere else, I decided to abandon pygmsh and use the python gmsh interface instead. I then followed these two examples and got what I needed: