Using mesher region information / tags

On Wed, Mar 23, 2011 at 10:26:57PM -0000, Neilen Marais wrote:
> New question #150226 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/150226
>
> Hi,
>
> Is there a way I can access mesh regions defined in a gmsh mesh? For instance, if I define a .geo file with labled physical volumes:
>
> ....
> Volume(25) = {24};
> Surface Loop(26) = {14, 10, 12, 17};
> Volume(27) = {26};
> Physical Volume ( 1000 ) = { 25 };
> Physical Volume ( 2000 ) = { 27 };
>
> with the resulting .msh file:
>
> $Elements
> 10
> 1 4 3 1000 25 0 1 9 6 10
> ....
> 5 4 3 1000 25 0 10 6 2 8
> 6 4 3 2000 27 0 6 9 3 7
> ....
> 10 4 3 2000 27 0 6 9 7 8
> $EndElements
>
> IOW, tets 1 through 5 have label 1000, and 6 through 10 has label 2000.
>
> After putting the .msh file through dolfin-convert the label data seems to be lost. Is this simply a limitation of dolfin-convert not supporting the gmsh feature, or does the dolfin mesh type not support labels of such a kind?
>
> Thanks
> Neilen

DOLFIN uses the MeshFunction class to store subdomain labels. If there
are no MeshFunctions defined in the <data> section of the resulting
XML file, then it's not supported by the dolfin-convert script.

If you find it's missing and know how to fix it consider submitting a
patch.

--
Anders

Anders,

I think I could cook up a patch, but I have searched (in vain) for a description of the dolfin xml format. If you could point me in the right direction I'm sure I could have a patch done in no time!

Thanks
Neilen

It looks like abacus supports writing of Meshfunctions for materials and such.

Look in:

  site-packages/dolfin/mesh/meshconverter.py

  ...
  def _abaqus(...):
    ...
    handler.start_meshfunction("material", 3, num_entities)
    # Each material is associated with a number of element sets
    for i, matname in enumerate(materials):
        try: elsetids = material2elsetids[matname]
        except KeyError:
            # No elements for this material
            continue
        # For each element set associated with this material
        elsets = []
        for eid in elsetids:
            try: elsets.append(eid2elset[eid])
            except KeyError:
                handler.error("Material '%s' is assigned to undefined element
"
                        "set '%s'" % (matname, eid))
        for elset in elsets:
            for elemid in elset:
                handler.add_entity_meshfunction(elemids.index(elemid), i)
    handler.end_meshfunction()

I guess something like this needs to be added to the

  def gmsh2xml

function.

Johan

On Wednesday March 23 2011 15:59:05 Neilen Marais wrote:
> Question #150226 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/150226
>
> Neilen Marais posted a new comment:
> Anders,
>
> I think I could cook up a patch, but I have searched (in vain) for a
> description of the dolfin xml format. If you could point me in the right
> direction I'm sure I could have a patch done in no time!
>
> Thanks
> Neilen

On Wed, Mar 23, 2011 at 10:59:05PM -0000, Neilen Marais wrote:
> Question #150226 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/150226
>
> Neilen Marais posted a new comment:
> Anders,
>
> I think I could cook up a patch, but I have searched (in vain) for a
> description of the dolfin xml format. If you could point me in the right
> direction I'm sure I could have a patch done in no time!

Just create a mesh and save it to XML to see what it looks
like. Attach some MeshFunctions to the mesh to see how they are stored:

from dolfin import *
mesh = UnitSquare(2, 2)
mesh.create_mesh_function("cell_domains", 2) # 2 = dim of cells in 2D
mesh.create_mesh_function("exterior_facet_domains", 1) # 1 = dim of facets in 2D

file = File("mesh.xml")
file << mesh

--
Anders

Note that "material" is not recognized by DOLFIN. It will be stored
but there is nothing in DOLFIN that will kick in an do anything with
this data (although it can be used from a user program).

A list of which data is recognized by DOLFIN (like "cell_domains") is
listed in MeshData.h:

/// The class MeshData is a container for auxiliary mesh data,
  /// represented either as MeshFunctions over topological mesh
  /// entities, arrays or maps. Each dataset is identified by a unique
  /// user-specified string. Only uint-valued data are currently
  /// supported.
  ///
  /// The following named mesh data are recognized by DOLFIN:
  ///
  /// Boundary indicators
  ///
  /// "boundary_facet_cells" - Array<uint> of size num_facets
  /// "boundary_facet_numbers" - Array<uint> of size num_facets
  /// "boundary_indicators" - Array<uint> of size num_facets
  /// "material_indicators" - MeshFunction<uint> of dimension D
  ///
  /// Subdomain indicators
  ///
  /// "cell_domains" - MeshFunction<uint> of dimension D
  /// "interior_facet_domains" - MeshFunction<uint> of dimension D - 1
  /// "exterior_facet_domains" - MeshFunction<uint> of dimension D - 1
  ///
  /// Facet orientation (used for assembly over interior facets)
  ///
  /// "facet orientation" - MeshFunction<uint> of dimension D - 1
  ///
  /// Boundary extraction
  ///
  /// "vertex map" - MeshFunction<uint> of dimension 0
  /// "cell map" - MeshFunction<uint> of dimension D
  ///
  /// Mesh partitioning
  ///
  /// "global entity indices %d" - MeshFunction<uint> of dimension 0, 1, ..., D
  /// "exterior facets" - MeshFunction<uint> of dimension D - 1
  /// "num global entities" - Array<uint> of size D + 1
  /// "overlap" - vector mapping
  ///
  /// Sub meshes
  ///
  /// "global vertex indices" - MeshFunction<uint> of dimension 0
  ///
  /// Mesh coloring
  ///
  /// "colors-%D-%d-%1" - MeshFunction<uint> of dimension D with colors based on connectivity %d
  /// "num colored cells" - Array<uint> listing the number of cells of each color
  /// "colored cells %d" - Array<uint> of cell indices with colors 0, 1, 2, ...

--
Anders

On Wed, Mar 23, 2011 at 04:12:40PM -0700, Johan Hake wrote:
> It looks like abacus supports writing of Meshfunctions for materials and such.
>
> Look in:
>
> site-packages/dolfin/mesh/meshconverter.py
>
> ...
> def _abaqus(...):
> ...
> handler.start_meshfunction("material", 3, num_entities)
> # Each material is associated with a number of element sets
> for i, matname in enumerate(materials):
> try: elsetids = material2elsetids[matname]
> except KeyError:
> # No elements for this material
> continue
> # For each element set associated with this material
> elsets = []
> for eid in elsetids:
> try: elsets.append(eid2elset[eid])
> except KeyError:
> handler.error("Material '%s' is assigned to undefined element
> "
> "set '%s'" % (matname, eid))
> for elset in elsets:
> for elemid in elset:
> handler.add_entity_meshfunction(elemids.index(elemid), i)
> handler.end_meshfunction()
>
>
> I guess something like this needs to be added to the
>
> def gmsh2xml
>
> function.
>
> Johan
>
> On Wednesday March 23 2011 15:59:05 Neilen Marais wrote:
> > Question #150226 on DOLFIN changed:
> > https://answers.launchpad.net/dolfin/+question/150226
> >
> > Neilen Marais posted a new comment:
> > Anders,
> >
> > I think I could cook up a patch, but I have searched (in vain) for a
> > description of the dolfin xml format. If you could point me in the right
> > direction I'm sure I could have a patch done in no time!
> >
> > Thanks
> > Neilen
>
> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~dolfin
> More help : https://help.launchpad.net/ListHelp

The abaqus plugin does not write anything to the mesh data structure. It
generates a an additional file.

I think this is what mesh converter should do. It creates the file
materials.xml, but that guy is just a MeshFunction over cells.

I have not tested this but this is my impression by just scanning the code.

Johan

On Wednesday March 23 2011 16:26:10 Anders Logg wrote:
> Note that "material" is not recognized by DOLFIN. It will be stored
> but there is nothing in DOLFIN that will kick in an do anything with
> this data (although it can be used from a user program).
>
> A list of which data is recognized by DOLFIN (like "cell_domains") is
> listed in MeshData.h:
>
> /// The class MeshData is a container for auxiliary mesh data,
> /// represented either as MeshFunctions over topological mesh
> /// entities, arrays or maps. Each dataset is identified by a unique
> /// user-specified string. Only uint-valued data are currently
> /// supported.
> ///
> /// The following named mesh data are recognized by DOLFIN:
> ///
> /// Boundary indicators
> ///
> /// "boundary_facet_cells" - Array<uint> of size num_facets
> /// "boundary_facet_numbers" - Array<uint> of size num_facets
> /// "boundary_indicators" - Array<uint> of size num_facets
> /// "material_indicators" - MeshFunction<uint> of dimension D
> ///
> /// Subdomain indicators
> ///
> /// "cell_domains" - MeshFunction<uint> of dimension D
> /// "interior_facet_domains" - MeshFunction<uint> of dimension D - 1
> /// "exterior_facet_domains" - MeshFunction<uint> of dimension D - 1
> ///
> /// Facet orientation (used for assembly over interior facets)
> ///
> /// "facet orientation" - MeshFunction<uint> of dimension D - 1
> ///
> /// Boundary extraction
> ///
> /// "vertex map" - MeshFunction<uint> of dimension 0
> /// "cell map" - MeshFunction<uint> of dimension D
> ///
> /// Mesh partitioning
> ///
> /// "global entity indices %d" - MeshFunction<uint> of dimension 0, 1,
> ..., D /// "exterior facets" - MeshFunction<uint> of dimension
> D - 1 /// "num global entities" - Array<uint> of size D + 1
> /// "overlap" - vector mapping
> ///
> /// Sub meshes
> ///
> /// "global vertex indices" - MeshFunction<uint> of dimension 0
> ///
> /// Mesh coloring
> ///
> /// "colors-%D-%d-%1" - MeshFunction<uint> of dimension D with
> colors based on connectivity %d /// "num colored cells" - Array<uint>
> listing the number of cells of each color /// "colored cells %d" -
> Array<uint> of cell indices with colors 0, 1, 2, ...
>
>
> --
> Anders
>
> On Wed, Mar 23, 2011 at 04:12:40PM -0700, Johan Hake wrote:
> > It looks like abacus supports writing of Meshfunctions for materials and
> > such.
> >
> > Look in:
> > site-packages/dolfin/mesh/meshconverter.py
> >
> > ...
> >
> > def _abaqus(...):
> > ...
> > handler.start_meshfunction("material", 3, num_entities)
> > # Each material is associated with a number of element sets
> >
> > for i, matname in enumerate(materials):
> > try: elsetids = material2elsetids[matname]
> >
> > except KeyError:
> > # No elements for this material
> > continue
> >
> > # For each element set associated with this material
> > elsets = []
> >
> > for eid in elsetids:
> > try: elsets.append(eid2elset[eid])
> >
> > except KeyError:
> > handler.error("Material '%s' is assigned to undefined
> > element
> >
> > "
> >
> > "set '%s'" % (matname, eid))
> >
> > for elset in elsets:
> > for elemid in elset:
> > handler.add_entity_meshfunction(elemids.index(elemid), i)
> >
> > handler.end_meshfunction()
> >
> > I guess something like this needs to be added to the
> >
> > def gmsh2xml
> >
> > function.
> >
> > Johan
> >
> > On Wednesday March 23 2011 15:59:05 Neilen Marais wrote:
> > > Question #150226 on DOLFIN changed:
> > > https://answers.launchpad.net/dolfin/+question/150226
> > >
> > > Neilen Marais posted a new comment:
> > > Anders,
> > >
> > > I think I could cook up a patch, but I have searched (in vain) for a
> > > description of the dolfin xml format. If you could point me in the
> > > right direction I'm sure I could have a patch done in no time!
> > >
> > > Thanks
> > > Neilen
> >
> > _______________________________________________
> > Mailing list: https://launchpad.net/~dolfin
> > Post to : <email address hidden>
> > Unsubscribe : https://launchpad.net/~dolfin
> > More help : https://help.launchpad.net/ListHelp

On Wed, Mar 23, 2011 at 04:36:31PM -0700, Johan Hake wrote:
> The abaqus plugin does not write anything to the mesh data structure. It
> generates a an additional file.
>
> I think this is what mesh converter should do. It creates the file
> materials.xml, but that guy is just a MeshFunction over cells.

It would be better to store it as part of the mesh file, at least in
the case when the data has some meaning to DOLFIN.

--
Anders

> I have not tested this but this is my impression by just scanning the code.
>
> Johan
>
> On Wednesday March 23 2011 16:26:10 Anders Logg wrote:
> > Note that "material" is not recognized by DOLFIN. It will be stored
> > but there is nothing in DOLFIN that will kick in an do anything with
> > this data (although it can be used from a user program).
> >
> > A list of which data is recognized by DOLFIN (like "cell_domains") is
> > listed in MeshData.h:
> >
> > /// The class MeshData is a container for auxiliary mesh data,
> > /// represented either as MeshFunctions over topological mesh
> > /// entities, arrays or maps. Each dataset is identified by a unique
> > /// user-specified string. Only uint-valued data are currently
> > /// supported.
> > ///
> > /// The following named mesh data are recognized by DOLFIN:
> > ///
> > /// Boundary indicators
> > ///
> > /// "boundary_facet_cells" - Array<uint> of size num_facets
> > /// "boundary_facet_numbers" - Array<uint> of size num_facets
> > /// "boundary_indicators" - Array<uint> of size num_facets
> > /// "material_indicators" - MeshFunction<uint> of dimension D
> > ///
> > /// Subdomain indicators
> > ///
> > /// "cell_domains" - MeshFunction<uint> of dimension D
> > /// "interior_facet_domains" - MeshFunction<uint> of dimension D - 1
> > /// "exterior_facet_domains" - MeshFunction<uint> of dimension D - 1
> > ///
> > /// Facet orientation (used for assembly over interior facets)
> > ///
> > /// "facet orientation" - MeshFunction<uint> of dimension D - 1
> > ///
> > /// Boundary extraction
> > ///
> > /// "vertex map" - MeshFunction<uint> of dimension 0
> > /// "cell map" - MeshFunction<uint> of dimension D
> > ///
> > /// Mesh partitioning
> > ///
> > /// "global entity indices %d" - MeshFunction<uint> of dimension 0, 1,
> > ..., D /// "exterior facets" - MeshFunction<uint> of dimension
> > D - 1 /// "num global entities" - Array<uint> of size D + 1
> > /// "overlap" - vector mapping
> > ///
> > /// Sub meshes
> > ///
> > /// "global vertex indices" - MeshFunction<uint> of dimension 0
> > ///
> > /// Mesh coloring
> > ///
> > /// "colors-%D-%d-%1" - MeshFunction<uint> of dimension D with
> > colors based on connectivity %d /// "num colored cells" - Array<uint>
> > listing the number of cells of each color /// "colored cells %d" -
> > Array<uint> of cell indices with colors 0, 1, 2, ...
> >
> >

On Wednesday March 23 2011 16:45:42 Anders Logg wrote:
> On Wed, Mar 23, 2011 at 04:36:31PM -0700, Johan Hake wrote:
> > The abaqus plugin does not write anything to the mesh data structure. It
> > generates a an additional file.
> >
> > I think this is what mesh converter should do. It creates the file
> > materials.xml, but that guy is just a MeshFunction over cells.
>
> It would be better to store it as part of the mesh file, at least in
> the case when the data has some meaning to DOLFIN.

Sure, but that is not supported by meshconverter for now. It looks like it
just write a MeshFunction to file, and that's it.

Johan

>
> --
> Anders
>
> > I have not tested this but this is my impression by just scanning the
> > code.
> >
> > Johan
> >
> > On Wednesday March 23 2011 16:26:10 Anders Logg wrote:
> > > Note that "material" is not recognized by DOLFIN. It will be stored
> > > but there is nothing in DOLFIN that will kick in an do anything with
> > > this data (although it can be used from a user program).
> > >
> > > A list of which data is recognized by DOLFIN (like "cell_domains") is
> > > listed in MeshData.h:
> > >
> > > /// The class MeshData is a container for auxiliary mesh data,
> > >
> > > /// represented either as MeshFunctions over topological mesh
> > > /// entities, arrays or maps. Each dataset is identified by a unique
> > > /// user-specified string. Only uint-valued data are currently
> > > /// supported.
> > > ///
> > > /// The following named mesh data are recognized by DOLFIN:
> > > ///
> > > /// Boundary indicators
> > > ///
> > > /// "boundary_facet_cells" - Array<uint> of size num_facets
> > > /// "boundary_facet_numbers" - Array<uint> of size num_facets
> > > /// "boundary_indicators" - Array<uint> of size num_facets
> > > /// "material_indicators" - MeshFunction<uint> of dimension D
> > > ///
> > > /// Subdomain indicators
> > > ///
> > > /// "cell_domains" - MeshFunction<uint> of dimension D
> > > /// "interior_facet_domains" - MeshFunction<uint> of dimension D -
> > > 1 /// "exterior_facet_domains" - MeshFunction<uint> of dimension D
> > > - 1 ///
> > > /// Facet orientation (used for assembly over interior facets)
> > > ///
> > > /// "facet orientation" - MeshFunction<uint> of dimension D - 1
> > > ///
> > > /// Boundary extraction
> > > ///
> > > /// "vertex map" - MeshFunction<uint> of dimension 0
> > > /// "cell map" - MeshFunction<uint> of dimension D
> > > ///
> > > /// Mesh partitioning
> > > ///
> > > /// "global entity indices %d" - MeshFunction<uint> of dimension 0,
> > > 1,
> > >
> > > ..., D /// "exterior facets" - MeshFunction<uint> of
> > > dimension D - 1 /// "num global entities" - Array<uint> of size
> > > D + 1
> > >
> > > /// "overlap" - vector mapping
> > > ///
> > > /// Sub meshes
> > > ///
> > > /// "global vertex indices" - MeshFunction<uint> of dimension 0
> > > ///
> > > /// Mesh coloring
> > > ///
> > > /// "colors-%D-%d-%1" - MeshFunction<uint> of dimension D with
> > >
> > > colors based on connectivity %d /// "num colored cells" - Array<uint>
> > > listing the number of cells of each color /// "colored cells %d" -
> > > Array<uint> of cell indices with colors 0, 1, 2, ...

Hi,

Running the example code given by Anders to generate a mesh function I get:

AttributeError: 'UnitSquare' object has no attribute 'create_mesh_function'

This is using the newest fenics from the PPA.

Another question is what to call the mesh function. The data coming from gmsh has no specific meaning, it is simply an identifier that the FEM code can do anything with. One obvious application would be the identification of materials, but it could conceivably also be used to indicate a solution method to apply, indication of a boundary condition to apply (when defined on facets) etc.

I would simply call it "group id" or similar. thoughts?

Sorry, there was a typo in the code I sent. Try this:

from dolfin import *
mesh = UnitSquare(2, 2)
mesh.data().create_mesh_function("cell_domains", 2) # 2 = dim of cells in 2D
mesh.data().create_mesh_function("exterior_facet_domains", 1) # 1 = dim of facets in 2D

file = File("mesh.xml")
file << mesh

I think "group id" sounds too generic. If it has a specific meaning (like indicators for domains or material parameters) then it might be natural to write that data as part of the mesh XML file in the <data> section. If it has no meaning, I would just write it to a separate XML file (containing a MeshFunction) and name that file something like foo_data.xml where foo is the name of the mesh file.

Anders,

Thanks for the updated code, it works.

A quick question. Once I have the properties written to a supplementary data .xml file, how do I read them in?

On Thu, Mar 24, 2011 at 12:44:44PM -0000, Neilen Marais wrote:
> Question #150226 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/150226
>
> Neilen Marais posted a new comment:
> Anders,
>
> Thanks for the updated code, it works.
>
> A quick question. Once I have the properties written to a supplementary
> data .xml file, how do I read them in?

Just do

  foo = MeshFunction("foo.xml")

--
Anders

Anders,

MeshFunction('foo.xml') doesnt seem to work for me. What does seem to work is

mf = dol.MeshFunction('uint', mesh, 'foo.xml')

where 'uint' is the mesh function type. The docstring for MeshFunction also seems to be a little confusing, since it seems to imply that you need a four parameter call. If I want to update a doc string is the best to check out the full fenics, or can I just send a patch to the ML?

On Thursday March 24 2011 08:11:35 Neilen Marais wrote:
> Question #150226 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/150226
>
> Neilen Marais posted a new comment:
> Anders,
>
> MeshFunction('foo.xml') doesnt seem to work for me.

That is because it is wrong. I guess Anders, has stopped using meshfunctions
stored to file, and only uses mesh functions stored in the mesh ;)

> What does seem to work is
>
> mf = dol.MeshFunction('uint', mesh, 'foo.xml')

Which should be the correct way to do it.

> where 'uint' is the mesh function type. The docstring for MeshFunction
> also seems to be a little confusing, since it seems to imply that you
> need a four parameter call. If I want to update a doc string is the best
> to check out the full fenics, or can I just send a patch to the ML?

I agree that is is comfusing, but it says that all the arguments but the first
string argument are optional. The rest of the arguments are just passed to the
corresponding C++ constructor. The dim argument is not used when instantiating
a MeshFunction from file, as that information is stored in the file. But it is
used when constructing an empty MeshFunction.

The docstring can made cleared with examples of the different ways you can
construct a MeshFunction. Patches to:

  dolfin/swig/mesh_post.i

are welcomed!

Johan

On Thu, Mar 24, 2011 at 03:11:35PM -0000, Neilen Marais wrote:
> Question #150226 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/150226
>
> Neilen Marais posted a new comment:
> Anders,
>
> MeshFunction('foo.xml') doesnt seem to work for me. What does seem to
> work is
>
> mf = dol.MeshFunction('uint', mesh, 'foo.xml')

Yes, sorry. Even I can't remember the interface.

> where 'uint' is the mesh function type. The docstring for MeshFunction
> also seems to be a little confusing, since it seems to imply that you
> need a four parameter call. If I want to update a doc string is the best
> to check out the full fenics, or can I just send a patch to the ML?

The best would be to check out check out the latest development
version of DOLFIN, fix the docstring, install it and verify it works,
then send it as a patch or bzr bundle.

--
Anders

Hi.

How did this work out? Is it possible to get hold of the updated dolfin-convert file?

Torbjørn Bækø Ness

I'm sorry, I just realized what was wrong....