Error in boundary conditions when importing Gmsh msh file

General
1

Hi, I am facing an issue while importing Gmsh msh file which I converted to xdmf format using meshio, but when I define boundary conditions on top and bottom side of the geometry it gives an error -

*** Warning: Found no facets matching domain for boundary condition.
*** Warning: Found no facets matching domain for boundary condition.

Minimum Working Example -

import meshio

mesh_from_file = meshio.read("mesh.msh")

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

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

from dolfin import *

mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)

class top(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 1., 0.1)
class bottom(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 0., 0.1)

boundaries = MeshFunction("size_t", mesh, mesh.geometry().dim() - 1)

top().mark(boundaries, 1)
bottom().mark(boundaries, 2)

V = VectorFunctionSpace(mesh, "Lagrange", 1)
bc_top = DirichletBC(V.sub(1), Constant(2.), boundaries, 1)
bc_bottom = DirichletBC(V, Constant((0., 0.)), boundaries, 2)

u = Function(V)
bc_top.apply(u.vector())
bc_bottom.apply(u.vector())

with XDMFFile("bcs.xdmf") as xdmf:
    xdmf.write(u)

Here is my .msh file.

Thank You.

2

Something is really strange with your msh file, as it seems like dolfin cannot find a single exterior facet in your mesh, i.e. assemble(1*ds(domain=mesh)) is equal to zero.
This can also be check by adding a print of the variable on_boundary in any of your inside functions or by calling:

mesh.init(mesh.topology().dim()-1, mesh.topology().dim())
f_to_c = mesh.topology()(mesh.topology().dim()-1,mesh.topology().dim())
for facet in facets(mesh):
    print(facet.index(), f_to_c(facet.index()), len(f_to_c(facet.index())))

What is the geo file for your mesh?

3

This is my geo file -

// Gmsh project created on Mon Jun 26 15:54:04 2023
SetFactory("OpenCASCADE");
//+
Rectangle(1) = {0, 0, 0, 1, 1, 0};
//+
Curve Loop(3) = {4, 1, 2, 3};
//+
Plane Surface(3) = {3};

I ran this snippet and the output is -

0 [116 278] 2
1 [113 275] 2
2 [113 116] 2
3 [275 278] 2
4 [117 279] 2
5 [118 280] 2

and so on till 485 [318 319] 2 .
I have a doubt that isn’t cell index() supposed to give an integer value then why it is giving output like [116 278] ?
Thanks for the help.

4

Add

Physical Surface(1) = {3};

at the end of your file.

This resolves the issues as in your current code meshio reads in more data than it should

This is not the output of an index command, but

which tells you want cells are connected to a facet.

5

Thanks @dokken , now it’s working fine.

So, this gives me facet is connecting which two cells right ?

6

It gives you either one or two cells, depending on if the facet is in the interior of the mesh or is an exterior facet

7

Ok, understood. Thanks, by the way, can you please have a look at this post Unexpected results in the application of Hertzian contact tutorial ? I think I am unable to define correct boundary conditions due to which the results are not what I was expecting. Happy to learn where I am making mistake.

8

I work on this forum on my spare time, and I’ve got limited capacity as to look at longer questions.
If you can reduce the problem, i.e. pin down where things go wrong in your code (is it the BCs, the contact expression or the solver) I might be able to help you