Boundary Conditions Not Matching when using GMSH imported msh file

Hello, I am having difficulty while applying DirichletBC on the top and bottom surface of my domain.

class top(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 3., DOLFIN_EPS)

class bottom(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], -3., DOLFIN_EPS)

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

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

bc_top = DirichletBC(V.sub(1), Constant(-0.1), boundaries, 1)
bc_bottom = DirichletBC(V, Constant((0., 0.)), boundaries, 2)

bcs = [bc_bottom, bc_top]

The above boundary conditions work fine when I am using the below code:

from dolfin import *
from mshr import *
import matplotlib.pyplot as plt

domain = Rectangle(Point((-2, 1)), Point((2, 3))) + Rectangle(Point(-2, -3), Point( 2, -1 )) + Circle(Point(0, 0), 1) - Circle(Point(0, 0), 0.5)
mesh = generate_mesh(domain, 200)

dolfin.plot(mesh)
plt.show()

But whenever I import my GMSH created msh file, the error comes up:

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

Here is my .msh file. I appreciate any input regarding same.
Thank You

I would suggest you increase the third argument in near to a lager value.

@dokken I increased it to 0.1, but it doesn’t work. I want to ask that is there any chance that I’m importing the msh file wrongly. Could you please confirm ? Here is how I’m importing it:

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

import numpy
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("2d_mesh.xdmf", triangle_mesh)

from dolfin import * 
mesh = Mesh()
with XDMFFile("2d_mesh.xdmf") as infile:
    infile.read(mesh)
    
plot(mesh);

I can’t reproduce this issue with the following Gmsh model, converting the mesh using the code you provided, and marking the boundaries using the following minimum working example (MWE). What is the output when you run the following?

import meshio
mesh = meshio.read("mesh.msh")
print(mesh.cells_dict.keys())

Working Gmsh model

SetFactory("OpenCASCADE");
Rectangle(1) = {-2, 1, 0, 4, 2, 0};
Rectangle(2) = {-2, -3, 0, 4, 2, 0};
Disk(3) = {0, 0, 0, 1, 1};
Disk(4) = {0, 0, 0, 0.5, 0.5};
BooleanDifference{ Surface{3}; Delete; }{ Surface{4}; Delete; }
BooleanFragments{ Surface{1}; Surface{3}; Delete; }{ Surface{2}; Delete; }
Field[1] = Constant;
Field[1].SurfacesList = {1, 2, 3};
Field[1].VIn = 0.1;
Field[1].VOut = 1;
Background Field = 1;

Minimum Working Example

from dolfin import *
mesh = Mesh()
with XDMFFile("2d_mesh.xdmf") as infile:
    infile.read(mesh)

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

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(-0.1), boundaries, 1)
bc_bottom = DirichletBC(V, Constant((0., 0.)), boundaries, 2)
2 Likes

Thanks @conpierce8 for the reply, your Gmsh model file solved my problem, now I understand this is the issue of Gmsh, whenever I built a geometry and mesh it and export it to .msh file, this boundary conditions error of finding no facets occur, when I used your Gmsh model it happens to be working fine.

Now I have a question - What should I do when I have to built a geometry because Gmsh generated .msh files are error prone and one more thing is I want to import 2d or 3d geometries which I am not able to built in Gmsh, I just want to import my stl or step file and mesh it, how do I go about it ?

Kindly tell which version of Gmsh or any other alternative would solve the issue.
Thank you.

Can you share your .geo file? I don’t usually find Gmsh to be error-prone. I suspect it is an error in your use of meshio.

Can you share the output of the above commands for your mesh file?

dict_keys(['vertex', 'line', 'triangle'])

This is my .geo file -

// Gmsh project created on Wed Jun 21 07:55:41 2023
SetFactory("OpenCASCADE");
//+
Circle(1) = {0, 0, 0, 0.5, 0, 2*Pi};
//+
Circle(2) = {0, 0, 0, 1, 0, 2*Pi};
//+
Rectangle(1) = {-1, 1, 0, 2, 2, 0};
//+
Rectangle(2) = {-1, -3, 0, 2, 2, 0};
//+
Curve Loop(3) = {6, 3, 4, 5};
//+
Curve Loop(4) = {6, 3, 4, 5};
//+
Curve Loop(5) = {2};
//+
Curve Loop(6) = {2};
//+
Curve Loop(7) = {1};
//+
Curve Loop(8) = {1};
//+
Curve Loop(9) = {9, 10, 7, 8};
//+
Curve Loop(10) = {9, 10, 7, 8};
//+
//+
Curve Loop(11) = {6, 3, 4, 5};
//+
Plane Surface(3) = {11};
//+
Curve Loop(12) = {2};
//+
Curve Loop(13) = {1};
//+
Plane Surface(4) = {12, 13};
//+
Transfinite Curve {1, 2, 3, 9} = 30 Using Progression 1;
//+
Transfinite Curve {1, 2, 2, 9, 3} = 100 Using Progression 1;

When I exported .msh file from it and ran it gave me the following error -

The issue is only with the top boundary: when I commenting out bc_top, the warning disappears.

You have some redundant commands in the .geo file which lead to the top rectangle being duplicated. It appears that Gmsh meshes the two copies of the top rectangle conformally, i.e. the two rectangles are treated as being connected on their boundaries. Therefore, the mesh facets (lines) that fall on the rectangle boundaries are treated as internal facets, because each of these facets is connected to at least two cells of the mesh. Therefore, on_boundary==False for these edges, and they are not matched by the DirichletBC. With the redundant commands removed as in the following, the BCs are marked correctly.

Corrected .geo file

SetFactory("OpenCASCADE");
//+
Circle(1) = {0, 0, 0, 0.5, 0, 2*Pi};
//+
Circle(2) = {0, 0, 0, 1, 0, 2*Pi};
//+
Rectangle(1) = {-1, 1, 0, 2, 2, 0};
//+
Rectangle(2) = {-1, -3, 0, 2, 2, 0};
//+
Curve Loop(3) = {6, 3, 4, 5};
//+
//Curve Loop(4) = {6, 3, 4, 5}; // Duplicate of Curve Loop(3)
//+
Curve Loop(5) = {2};
//+
//Curve Loop(6) = {2}; // Duplicate of Curve Loop(5)
//+
Curve Loop(7) = {1};
//+
//Curve Loop(8) = {1}; // Duplicate of Curve Loop(7)
//+
Curve Loop(9) = {9, 10, 7, 8};
//+
//Curve Loop(10) = {9, 10, 7, 8}; // Duplicate of Curve Loop(9)
//+
//+
//Curve Loop(11) = {6, 3, 4, 5}; // Duplicate of Curve Loop(3)
//+
//Plane Surface(3) = {11}; // Not necessary; this surface already exists from the Rectangle(1) command
//+
Curve Loop(12) = {2};
//+
Curve Loop(13) = {1};
//+
Plane Surface(4) = {12, 13};
//+
//Transfinite Curve {1, 2, 3, 9} = 30 Using Progression 1; // Overridden by the following command
//+
Transfinite Curve {1, 2, 2, 9, 3} = 100 Using Progression 1;

Mesh conversion code

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

import numpy as np
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", True)
meshio.write("2d_mesh.xdmf", triangle_mesh)

DirichletBC verification

from dolfin import *
mesh = Mesh()
with XDMFFile("2d_mesh.xdmf") as infile:
    infile.read(mesh)

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

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((1., 2.)), boundaries, 2)

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

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

Yeah, that solved the issue, I was really oblivious to the fact that I have unknowingly duplicated some parts of mesh. Thank You very much @conpierce8 for the help.
But the two meshes (one that you provided via Gmsh Model Code and the Corrected Geo file) reacts differently when I give a displacement of -0.5 in both the cases.
Displacement Plot (when Gmsh Model Code msh file is used) -
image
Displacement Plot (when Corrected Geo file generated msh file is used) -
image
I am confused about what went wrong, if boundary conditions problem is solved it should give the same displacement plot in both the cases. Any inputs regarding the same would be appreciated.
Thank You.

I used a BooleanFragments in my .geo file, which effectively connects the three domains (two rectangles and one annulus). In the corrected .geo file, since this command is not present, the three domains are disconnected, i.e. they share no vertices or facets.

How you should proceed depends on the problem you are trying to solve. It looks like you are solving a contact problem between two half-spaces and an annulus. In this case, you should modify your weak form to enforce the no-penetration condition between the rectangle and the annulus (e.g. using a penalty formulation). If you need further help on setting up this problem, I would recommend creating a separate post rather than continuing in this thread.

2 Likes

Thanks @conpierce8 for the guidance, I will try to follow this tutorial on penalty approach and reach out to the community if got stuck anywhere.
By the way do you have more resources which I can follow to simulate 2D penalty formulation in my geometry ?
Thanks for your time.

I’ve never studied contact problems, so unfortunately I don’t know of any other resources.

Are the boxes supposed to be rigid or flexible?

Boxes are supposed to be rigid and the annulus region is to be kept soft such that boxes deform the annulus region like in the image shown below -


I want to simulate this that is as the displacement is increased, the area of contact between upper box and annulus increases.