Found no facets matching domain for boundary condition

Question

Hi, everyone.

I have read most of the articles regarding the "no factet matching" listed above, but still have no clue of what went wrong with my very simple program. I have 4 different taged lines on the geometry. It seems that somehow subdomain labeling system does not work at all.
Thank you for your input.

The geometry: 2dtest4.xml

cl__1 = 1;
Point(1) = {0, 0, 0, 1};
Point(2) = {10, 0, 0, 1};
Point(3) = {20, 0, 0, 1};
Point(4) = {0, 10, 0, 1};
Point(5) = {0, 20, 0, 1};
Line(1) = {2, 3};
Line(2) = {4, 5};
Circle(3) = {4, 1, 2};
Circle(4) = {5, 1, 3};
Line Loop(6) = {2, 4, -1, -3};
Plane Surface(7) = {6};
Physical Surface(1) = {7};
Physical Line(1) = {2};
Physical Line(2) = {1};
Physical Line(3) = {3};
Physical Line(4) = {4};

====================================================

from dolfin import *

set_log_level(1)

Create mesh and define function space

mesh = Mesh("2dtest4.xml")
boundaries = MeshFunction("size_t",mesh,mesh.topology().dim()-1)
boundaries.set_all(0)

V = FunctionSpace(mesh, "CG", 1)

Define variational problem

u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0.0)
a = dot(grad(u), grad(v))dx
L = fv*dx

Define boundary condition values

u1 = Constant(100.0)
u2 = Constant(50.0)
u3 = Constant(0.0)
u4 = Constant(0.0)

Define boundary conditions

bc1 = DirichletBC(V, u1, boundaries, 0) - last argument must match with number (or label) of Physical Surface(0)

bc1 = DirichletBC(V, u1, boundaries, 1)
bc2 = DirichletBC(V, u2, boundaries, 2)
bc3 = DirichletBC(V, u3, boundaries, 3)
bc4 = DirichletBC(V, u4, boundaries, 4)

Compute solution

u = Function(V)
solve(a == L, u, [bc1, bc2, bc3, bc4])

Plot solution

plot(u, interactive=True)

Answer 1

Hi, the problem is in the way you create your boundaries . First you create an empty FacetFunction over a mesh - it knows nothing about the labels from GMSH. Even if it did you erase them by calling set_all method. There is no wonder then, that the assembler finds no facets with labels 1, 2, etc. Consider the following

from dolfin import *

mesh = Mesh('2dtest4.xml')
tdim = mesh.topology().dim()

# Create facet function over mesh, unitialized values
boundaries = MeshFunction('size_t', mesh, tdim - 1)
plot(boundaries, title='unitialized values', interactive=True)

# Set all values to 0
boundaries.set_all(0)
plot(boundaries, title='0 values', interactive=True)

# Actually use facet labels from gmsh
boundaries = MeshFunction('size_t', mesh, '2dtest4_facet_region.xml')
plot(boundaries, title='using values frm GMSH', interactive=True) 

Comments

Thank you for the answer to my question.

I tried what the answer said and got the error message saying
"Error: Unable to read data from XML file.
*** Reason: Unable to open file "[mesh]_facet_region.xml"

When doing with dolfin-convert.py, I do not see " [mesh]_facet_region.xml" or "[mesh]_physical.xml" file. Only "[mesh].xml" is there. I am using Gmsh resulting in [mesh].geo to be converted into [mesh].xml. I did "python dolfin-convert.py [mesh].msh [mesh].xml". Am I still missing something?

Thank you.

---

Okay, this is a different issue. Let's first see if the original one is solved by the answer - could you check with the mesh and boundary files? These were generated from geo file with GMSH 2.8.5 and dolfin-convert from DOLFIN 1.4.0.

---

Thank you for your time and effort to help me out with this.

I did as what you told and it worked nicely!
I got a new dolfin-convert.py form DOLFIN 1.4.0, which generates both xml files.
It took a month to figure out what went wrong ouch!

For the sake of newbies like me,

1) Make sure of using the right version of dolfin-convert.py (from 1.4.0). I had a wrong one putting me in all the troubles.

2) Upgrade to GMSH 2.8.5

Thank you.