Is there a mapping of (cell number, local entity) to global edge index?

Question

I have a Mesh and a MeshValueCollection. I would like to end up with a MeshFunction. I could set the values in the MeshFunction if I could use the (cell number, local entity number) as an index into the edges.

# loading data
mesh = Mesh("test.xml")
mvc = MeshValueCollection("size_t", mesh, "sides.xml")

# empty mesh function
meshf = MeshFunction("size_t", mesh, 1, 0)

# I could iterate over the mesh value collection like this:
for cell_num, entity_num in mvc.values():
  # ? but I can't figure out how I could use these to index the mesh function ?

Answer 1

Have you tried:

meshf = MeshFunction("size_t", mesh, mvc)

Comments

Hi - thanks for looking at this question - I am using the snapshot Mac OS X 10.7 download from today, and here's what happens when I try that:

>>> mvc = MeshValueCollection("size_t", mesh, "test_facets.xml")
>>> meshf = MeshFunction("size_t", mesh, mvc)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Applications/FEniCS.app/Contents/Resources/lib/python2.7/site-packages/dolfin/cpp/mesh.py", line 6612, in __new__
    return MeshFunctionSizet(*args)
  File "/Applications/FEniCS.app/Contents/Resources/lib/python2.7/site-packages/dolfin/cpp/mesh.py", line 4955, in __init__
    _mesh.MeshFunctionSizet_swiginit(self,_mesh.new_MeshFunctionSizet(*args))
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics@fenicsproject.org
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to assign mesh value collection to mesh function.
*** Reason:  Mesh value collection does not contain all values for all entities.
*** Where:   This error was encountered inside MeshFunction.h.
*** Process: 0
*** DOLFIN version: 1.2.0+
*** Git changeset:  
*** -------------------------------------------------------------------------

---

It is raising an error because the MeshValueCollection is sparse and not assigning all values of the MeshFunction. I am not sure why it is raising an error. We could adapt the behaviour we have for creating a MeshFunction from MeshDomains, which also comes with a sparse data structure.

I will ask the fenics list what others think.

---

Ok, thanks for taking a look at it.

---

I have now pushed a fix for this in the master branch. It will soon be part of the DOLFIN 1.3 release.

---

Ok - would that appear in the snapshot release in the next few days or so?
At any rate, thanks very much for providing the sparse MeshValueCollection functionality ...

---

Hi,
I wrote an exporter myself, faced the same problem and I just could resolve it using your recent fix, Johan - thanks for that! However I am still wondering why the location is organized like that. Why can the location of a DirichletBC defined quite simple, however for the location of a NeumannBC I have to set up a dedicated MeshFunction. The BC demo is a good example: I just put the domain indexes into the mesh file such as

  <domains>
   <mesh_value_collection type="uint" dim="2" size="3226">
    <value cell_index="2" local_entity="3" value="0" />
    ...

And then I use in my code:

 DirichletBC bcl(V, fixed, 0);

without touching the subdomain definition. However to set up a Neumann BC I need a MeshFunction, for that I have to read in a dedicated MeshFunction file. Maybe I am missing something, but wouldn't it be easier to make the domains defined in a mesh file available for the assignment of NeumannBCs ? To make the question even more concrete, how would I assign the domains in the aneursm.xml mesh of the BC demo to a Neumann boundary condition? If that is already possible, it would be great to get a hint how to do that.

Thanks in advance and all the best,

David

---

I was under the impression that the project was moving away from having all of the domain information in the mesh file (I think that "dolfin-convert" calls this the xml-old format), I could very easily be wrong though.

It turns out that (I think) Johan answered a question similar to yours on SciComp (his answer is presently the one with the most votes). I tried it out in the Mac OS X snapshot dated 17-November (haven't bothered to update it, since everything I wanted to work presently works).

To summarize his answer (in Python):

mesh = Mesh("aneurysm.xml.gz")
facet_function = MeshFunction("size_t", mesh, 2, mesh.domains())
cell_function = MeshFunction("size_t", mesh, 3, mesh.domains())

---

xml-old format is something else. Previously we stored the values of MeshFunctions directly in the xml file, that did not work in parallel. The solution since 1.0.0 release is to use MeshValueCollections which uses cells and local cell entities as indices inside the XML file. The information stored in the MeshValueCollection is read into and stored in a MeshDomains class, available by Mesh::domains.
Then we have changed how one can generate a MeshFunction from the stored Mesh::domains. This is what I answered on StackExchange. Here I answered how one can generate a MeshFunction from a MeshValueCollection XML-file.