Setting boundary conditions for meshes generated by gmsh

Question

tldr: while converting msh file to xml dolfin-convert does not produces [filename]_facet_region.xml for square meshes.

I have the following gmsh geo file with square mesh:

    size = 1; //m
    meshThickness = size / 3; 
    gridsize = size / 5;

    // All numbering counterclockwise from bottom-left corner
    Point(1) = {0, 0, 0, gridsize};
    Point(2) = {0, size, 0, gridsize};
    Point(3) = {size, size, 0, gridsize};
    Point(4) = {size, 0, 0, gridsize};
    Line(1) = {1, 2};               
    Line(2) = {2, 3};           
    Line(3) = {3, 4};       
    Line(4) = {4, 1};       
    Line Loop(4) = {1, 2, 3, 4};
    Plane Surface(8) = {4};
    Physical Line(12) = {1};
    Physical Line(23) = {2};
    Physical Line(34) = {3};
    Physical Line(41) = {4};

    //Physical Surface(1) = {8};

    Transfinite Surface{8};
    Recombine Surface{8};

Mesh.CharacteristicLengthFromCurvature = 1;
Mesh.Algorithm = 8; // delannay=1

I want to set Dirichlet boundary conditions for physical lines 12, 23, 34, and 41.
I convert geo to msh file with

gmsh -3 square.geo

and later use

dolfin-convert quad.msh quad.xml 

It was suggested http://fenicsproject.org/qa/2986/how-to-define-boundary-condition-for-mesh-generated-by-gmsh to use

mesh = Mesh("yourmeshfile.xml")
subdomains = MeshFunction("size_t", mesh, "yourmeshfile_physical_region.xml")
boundaries = MeshFunction("size_t", mesh, "yourmeshfile_facet_region.xml")
bcs = [DirichletBC(V, 5.0, boundaries, 1),# of course with your boundary
    DirichletBC(V, 0.0, boundaries, 0)]

But I see no facet_region file! It exists in case of triangular one, though. How do I set up boundary condition in this case?

Answer 1

It works just fine for me (Version 1.6) without:

Transfinite Surface{8};
Recombine Surface{8};

Otherwise it seems to be a quadrilateral mesh which, according to
http://fenicsproject.org/qa/7496/can-dolfin-convert-work-on-quadrilaterals-from-gmsh,
is not supported yet (version 1.6 and especially by dolfin-convert, at least not the standard one).

I think firedrake has some support for quadrilaterals (i never used it) so maybe try that out (or some dev version).

Answer 2

I had the same problem for triangle meshes. It turns out that dolfin-convert ignores physical regions if you only define one. Split your geometry in two arbitrary regions, and define them as

Physical Surface(1) = ...
Physical Surface(2) = ...

and it might work.