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()