Defining boundary conditions for mesh generated in gmsh

Praveen:

If you want to use e.g. a tet mesh, and then define surface boundaries using tags, it is unfortunately not possible yet. I recently added some functionality to the dolfin mesh convert utility to import integer tags for regions (i.e. if you have a tet mesh, you can assign number to different volumetric regions). The integer tags are then written to a separate xml file as a mesh function.

To define 2D boundary regions on a 3D mesh you'll have to extend the functionality of dolfin-convert. At the moment it looks for the highest dimentional entity in a mesh file, and ignores the rest. This is done since, by default, gmsh writes out all line, surface and volume entities unless you define physical regions. You'll have to modify the behaviour, perhaps to check for the presence of tagged entities, and then write out an additional dim-1 mesh function for the boundary properties.

Modifying dolfin-convert is fairly easy (quite high level Python code); I can see one major difficulty though: mapping between the mesh numbering of gmsh and dolfin. Dolfin will number all the faces in the mesh, wheras gmsh will probably only number those that you define as boundaries.

Thanks for the answer. I always set physical types so that gmsh will export only the entities that I need, volume cells and boundary edges/faces. My python skill at present is very basic. So I may not have much success trying to modify dolfin-convert.

I would like to add the functionality as described by Neilen.

Can anybody give me information regarding the numbering of edges (in 2D) or facets (in 3D) when constructing the mesh?

Thanks

On Fri, Jun 17, 2011 at 05:41:16AM -0000, Evan Lezar wrote:
> Question #161172 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/161172
>
> Status: Open => Needs information
>
> Evan Lezar requested more information:
> I would like to add the functionality as described by Neilen.
>
> Can anybody give me information regarding the numbering of edges (in 2D)
> or facets (in 3D) when constructing the mesh?

Do you mean local or global numbering? (Relative to local cell or
global edge/face numbers.)

Local numbers follow the UFC numbering convention (same number as
opposite vertex, described in the FEniCS book) and global numbering is
undefined. It comes out as whatever the algorithm implemented in
TopologyComputation.cpp computes.

--
Anders

I think that I mean global numbering. I want to be able to construct a mesh-function that has the same values as physical regions specified in a gmsh mesh file. Specifically on boundary facets.

How is the order of the facets in the SubDomain.mark ( meshfunction, value) determined?

Ok, I see in SubDomain::mark_meshfunction, that the mesh is first initialised, before checking to see if the vertices of an entity fall in a subdomain.

I assume that the numbering of the facets of a given mesh are constant, otherwise it would not be possible to load a saved mesh function?

On Thursday June 16 2011 23:41:24 Evan Lezar wrote:
> Question #161172 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/161172
>
> Evan Lezar requested more information:
> I think that I mean global numbering. I want to be able to construct a
> mesh-function that has the same values as physical regions specified in
> a gmsh mesh file. Specifically on boundary facets.
>
> How is the order of the facets in the SubDomain.mark ( meshfunction,
> value) determined?

If you make your mesh converter script depend on DOLFIN you can use
dolfin::MeshEditor to create the mesh with the correct order. Then use the
connectivities to map boundary faces to a FacetFunction and save.

I do this in TriTetMesh,

  https://launchpad.net/tritetmesh

where I map bondary faces generated from Tetgen to a FacetFunction. The logic
is within the file:

  TetMesh::save_dolfin

in the file

  tritetmesh/tetmesh/tetmesh.cpp

Something similare have to be possible to do in Python.

Johan

On Fri, Jun 17, 2011 at 06:45:49AM -0000, Evan Lezar wrote:
> Question #161172 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/161172
>
> Evan Lezar posted a new comment:
> Ok, I see in SubDomain::mark_meshfunction, that the mesh is first
> initialised, before checking to see if the vertices of an entity fall in
> a subdomain.
>
> I assume that the numbering of the facets of a given mesh are constant,
> otherwise it would not be possible to load a saved mesh function?

In general, one can't rely on the numbering of facets in DOLFIN. The
numbering depends on the specific algorithm used to compute the vertex
numbers, and the numbering may change if that algorithm is optimized.

For a robust numbering, you need to relate the facets to which cell
they belong to, and the local number of the facet relative to that
cell as we do for storing boundary facet numbers as part of the DOLFIN
XML format. Look at the example aneurysm.xml.gz in data/meshes.

--
Anders

On Friday June 17 2011 10:41:23 Anders Logg wrote:
> Question #161172 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/161172
>
> Anders Logg proposed the following answer:
>
> On Fri, Jun 17, 2011 at 06:45:49AM -0000, Evan Lezar wrote:
> > Question #161172 on DOLFIN changed:
> > https://answers.launchpad.net/dolfin/+question/161172
> >
> > Evan Lezar posted a new comment:
> > Ok, I see in SubDomain::mark_meshfunction, that the mesh is first
> > initialised, before checking to see if the vertices of an entity fall in
> > a subdomain.
> >
> > I assume that the numbering of the facets of a given mesh are constant,
> > otherwise it would not be possible to load a saved mesh function?
>
> In general, one can't rely on the numbering of facets in DOLFIN. The
> numbering depends on the specific algorithm used to compute the vertex
> numbers, and the numbering may change if that algorithm is optimized.
>
> For a robust numbering, you need to relate the facets to which cell
> they belong to, and the local number of the facet relative to that
> cell as we do for storing boundary facet numbers as part of the DOLFIN
> XML format. Look at the example aneurysm.xml.gz in data/meshes.

That is true, but if you link to dolfin you will always get what ever mapping
dolfin is using there and then. I think it is much easier to work with a
global facetfunction, as you can use it for a lot of other boundary specific
things.

Johan

Johan

> --
> Anders

On Fri, Jun 17, 2011 at 06:01:41PM -0000, Johan Hake wrote:
> Question #161172 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/161172
>
> Johan Hake proposed the following answer:
> On Friday June 17 2011 10:41:23 Anders Logg wrote:
> > Question #161172 on DOLFIN changed:
> > https://answers.launchpad.net/dolfin/+question/161172
> >
> > Anders Logg proposed the following answer:
> >
> > On Fri, Jun 17, 2011 at 06:45:49AM -0000, Evan Lezar wrote:
> > > Question #161172 on DOLFIN changed:
> > > https://answers.launchpad.net/dolfin/+question/161172
> > >
> > > Evan Lezar posted a new comment:
> > > Ok, I see in SubDomain::mark_meshfunction, that the mesh is first
> > > initialised, before checking to see if the vertices of an entity fall in
> > > a subdomain.
> > >
> > > I assume that the numbering of the facets of a given mesh are constant,
> > > otherwise it would not be possible to load a saved mesh function?
> >
> > In general, one can't rely on the numbering of facets in DOLFIN. The
> > numbering depends on the specific algorithm used to compute the vertex
> > numbers, and the numbering may change if that algorithm is optimized.
> >
> > For a robust numbering, you need to relate the facets to which cell
> > they belong to, and the local number of the facet relative to that
> > cell as we do for storing boundary facet numbers as part of the DOLFIN
> > XML format. Look at the example aneurysm.xml.gz in data/meshes.
>
> That is true, but if you link to dolfin you will always get what ever mapping
> dolfin is using there and then. I think it is much easier to work with a
> global facetfunction, as you can use it for a lot of other boundary specific
> things.

Sure, but it can't be used for storing data to a file for later use.
(Yes it can, but it can break anytime.)

--
Anders

On Friday June 17 2011 11:11:29 Anders Logg wrote:
> Question #161172 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/161172
>
> Anders Logg proposed the following answer:
>
> On Fri, Jun 17, 2011 at 06:01:41PM -0000, Johan Hake wrote:
> > Question #161172 on DOLFIN changed:
> > https://answers.launchpad.net/dolfin/+question/161172
> >
> > Johan Hake proposed the following answer:
> >
> > On Friday June 17 2011 10:41:23 Anders Logg wrote:
> > > Question #161172 on DOLFIN changed:
> > > https://answers.launchpad.net/dolfin/+question/161172
> > >
> > > Anders Logg proposed the following answer:
> > >
> > > On Fri, Jun 17, 2011 at 06:45:49AM -0000, Evan Lezar wrote:
> > > > Question #161172 on DOLFIN changed:
> > > > https://answers.launchpad.net/dolfin/+question/161172
> > > >
> > > > Evan Lezar posted a new comment:
> > > > Ok, I see in SubDomain::mark_meshfunction, that the mesh is first
> > > > initialised, before checking to see if the vertices of an entity fall
> > > > in a subdomain.
> > > >
> > > > I assume that the numbering of the facets of a given mesh are
> > > > constant, otherwise it would not be possible to load a saved mesh
> > > > function?
> > >
> > > In general, one can't rely on the numbering of facets in DOLFIN. The
> > > numbering depends on the specific algorithm used to compute the vertex
> > > numbers, and the numbering may change if that algorithm is optimized.
> > >
> > > For a robust numbering, you need to relate the facets to which cell
> > > they belong to, and the local number of the facet relative to that
> > > cell as we do for storing boundary facet numbers as part of the DOLFIN
> > > XML format. Look at the example aneurysm.xml.gz in data/meshes.
> >
> > That is true, but if you link to dolfin you will always get what ever
> > mapping dolfin is using there and then. I think it is much easier to
> > work with a global facetfunction, as you can use it for a lot of other
> > boundary specific things.
>
> Sure, but it can't be used for storing data to a file for later use.
> (Yes it can, but it can break anytime.)

Doing it all the time :)

Johan

>
> --
> Anders

Thanks for all the input. I will see what I can come up with.

Will let you know how it goes.

Hi,

I have sent a patch of my initial attempt through to the mailing list.

Would like to hear if it works as expected!

Regards
Evan

I don't know if this is related, but I would like to use both Physical Volume and Physical Surface in the same mesh, but then dolfin-convert seems to get stuck after printing "Found all cells"?

Example gmsh file 'markedSquare.geo':
---
cl1 = .1;
Point(3) = {1, 1, 1, cl1};
Point(4) = {1, 1, 1, cl1};
Point(5) = {-1, 1, 1, cl1};
Point(6) = {-1, -1, 1, cl1};
Point(7) = {-1, -1, -1, cl1};
Point(8) = {1, -1, -1, cl1};
Point(9) = {1, 1, -1, cl1};
Point(10) = {1, -1, 1, cl1};
Point(11) = {-1, 1, -1, cl1};
Point(12) = {0, 1, 1, cl1};
Point(13) = {0, -1, 1, cl1};
Point(14) = {0, -1, -1, cl1};
Point(15) = {0, 1, -1, cl1};
Line(1) = {3, 9};
Line(3) = {11, 5};
Line(5) = {7, 6};
Line(7) = {10, 8};
Line(9) = {7, 11};
Line(10) = {6, 5};
Line(11) = {10, 3};
Line(12) = {8, 9};
Line(13) = {15, 12};
Line(14) = {12, 13};
Line(15) = {13, 14};
Line(16) = {14, 15};
Line(35) = {7, 14};
Line(36) = {14, 8};
Line(37) = {6, 13};
Line(38) = {13, 10};
Line(39) = {5, 12};
Line(40) = {12, 3};
Line(41) = {11, 15};
Line(42) = {15, 9};
Line Loop(44) = {37, 15, -35, 5};
Plane Surface(44) = {44};
Line Loop(46) = {15, 36, -7, -38};
Plane Surface(46) = {46};
Line Loop(47) = {1, -12, -7, 11};
Plane Surface(47) = {47};
Line Loop(48) = {13, 14, 15, 16};
Plane Surface(48) = {48};
Line Loop(50) = {16, 42, -12, -36};
Plane Surface(50) = {50};
Line Loop(52) = {9, 41, -16, -35};
Plane Surface(52) = {52};
Line Loop(53) = {5, 10, -3, -9};
Plane Surface(53) = {53};
Line Loop(55) = {41, 13, -39, -3};
Plane Surface(55) = {55};
Line Loop(57) = {40, 1, -42, 13};
Plane Surface(57) = {57};
Line Loop(59) = {11, -40, 14, 38};
Plane Surface(59) = {59};
Line Loop(61) = {39, 14, -37, 10};
Plane Surface(61) = {61};

Surface Loop(62) = {57, 59, 47, 50, 46, 48};
Surface Loop(64) = {53, 44, 61, 55, 52, 48};
Surface Loop(66) = {61, 55, 52, 53, 44, 47, 57, 59, 46, 50};

Volume(67) = {66};
//Volume(68) = {62};
//Volume(69) = {64};
//Physical Volume(71) = {69};

Physical Volume(222) = {67};
Physical Surface(111) = {61};
---
This gets stuck when I run it with
gmsh markedSquare.geo -3
dolfin-convert markedSquare.msh markedSquare.xml

if I comment out the last Physical Volume at the end, then
--
mesh = Mesh("markedSquare.xml")
plot(mesh)
interactive()
--
plots a two-dimensional mesh?