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)