from gmsh geo file to fenics mesh with subdomain and boundary facet exported

from dolfin import *
import os.path
from subprocess import Popen, PIPE, check_output

# https://fenicsproject.org/qa/9983/setting-boundary-conditions-for-meshes-generated-by-gmsh
# generete mesh by: gmsh test.geo -3
# dolfin-convert -i gmsh test.msh test.xml

geo_file_content = """
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};  // do not comment out this line

//Transfinite Surface{8};  // fenics does not support recomine
//Recombine Surface{8};


// boundary layer
Field[4] = BoundaryLayer;
Field[4].EdgesList = {1,3};  // FacesList  for 3D
Field[4].hfar = 10;
Field[4].hwall_n = 0.015;
Field[4].hwall_t = 0.01;
Field[4].thickness = 0.02;
Field[4].ratio = 1.0;

BoundaryLayer Field = 4;

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

fname="test"
with open(fname+".geo", 'w') as f:
    f.write(geo_file_content)

cmd = 'gmsh {}.geo -3'.format(fname)
print check_output([cmd], shell=True)  # run in shell mode in case you are not run in terminal

cmd = 'dolfin-convert -i gmsh {}.msh {}.xml'.format(fname, fname)
Popen([cmd], stdout=PIPE, shell=True).communicate()

mesh = Mesh(fname+".xml")
if os.path.exists( fname+"_physical_region.xml"):
    subdomains = MeshFunction("size_t", mesh, fname+"_physical_region.xml")
    plot(subdomains)
if os.path.exists( fname+"_facet_region.xml"):
    boundaries = MeshFunction("size_t", mesh, fname+"_facet_region.xml")
    plot(boundaries)

plot(mesh)
interactive()