How to use MeshValueCollection?

mesh
1

Hi,

It’s been several days that i am stuck on this problem : “Reason: Error reading MeshValueCollection.”

Please do you know how to fix it ?

Context: i have just imported my .geo file on fenics (google colabs), then i used this program :

!add-apt-repository -y ppa:fenics-packages/fenics

!apt-get update -qq

!apt-get install -y -qq software-properties-common python3-software-properties module-init-tools

!apt install -y -qq --no-install-recommends fenics

!pip3 install --upgrade pip

!pip3 install -U --no-binary=h5py h5py meshio

from dolfin import *

!apt-get install gmsh

!gmsh --version

!gmsh -2 wear.geo -o wear.msh

!pip install -U meshio

!apt-get install python-lxml

import meshio

in_mesh = meshio.read("wear.msh")

import numpy

def create_mesh(mesh, cell_type, prune_z=False):
    cells = numpy.vstack([cell.data for cell in mesh.cells if cell.type==cell_type])
    cell_data = numpy.hstack([mesh.cell_data_dict["gmsh:physical"][key]
                         for key in mesh.cell_data_dict["gmsh:physical"].keys() if key==cell_type])
    # Remove z-coordinates from mesh if we have a 2D cell and all points have the same third coordinate
    points= mesh.points
    if prune_z:
        points = points[:,:2]
    mesh_new = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    return mesh_new
line_mesh = create_mesh(in_mesh, "line", prune_z=True)
meshio.write("facet_mesh.xdmf", line_mesh)

triangle_mesh = create_mesh(in_mesh, "triangle", prune_z=True)
meshio.write("mesh.xdmf", triangle_mesh)

mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 1)
with XDMFFile("facet_mesh.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

and this is the error message :

---------------------------------------------------------------------------

RuntimeError Traceback (most recent call last)

[<ipython-input-10-1c9f148725b4>](https://localhost:8080/#) in <module>() 4 mvc = MeshValueCollection("size_t", mesh, 35) 5 with XDMFFile("facet_mesh.xdmf") as infile: ----> 6 infile.read(mvc, "name_to_read") 7 mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

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-support@googlegroups.com](mailto:fenics-support@googlegroups.com) *** *** Remember to include the error message listed below and, if possible, *** include a *minimal* running example to reproduce the error. *** *** ------------------------------------------------------------------------- *** Error: Unable to find entity in map. *** Reason: Error reading MeshValueCollection. *** Where: This error was encountered inside HDF5File.cpp. *** Process: 0 *** *** DOLFIN version: 2019.2.0.dev0 *** Git changeset: unknown *** -------------------------------------------------------------------------

SEARCH STACK OVERFLOW
2

I don’t know, but do you think it comes from the definition of my physical groups in the .geo file ?

Physical Line(1) = {1} ;
Physical Line(2) = {2} ;
Physical Line(3) = {20} ;
Physical Line(4) = {3} ;
Physical Line(5) = {6} ;
Physical Line(6) = {19} ;
Physical Line(7) = {7} ;
Physical Line(8) = {10} ;
Physical Line(9) = {11} ;
Physical Line(10) = {12} ;
Physical Line(11) = {18} ;
Physical Line(12) = {13} ;
Physical Line(13) = {9} ;
Physical Line(14) = {17} ;
Physical Line(15) = {8} ;
Physical Line(16) = {4} ;
Physical Line(17) = {55} ;
Physical Line(18) = {144} ;

Physical Surface("Surface") = {29};
Physical Surface("Surface") += {30};
Physical Surface("Surface") += {31};
Physical Surface("Surface") += {32};
Physical Surface("Surface") += {33};
Physical Surface("Surface") += {34};
3

i found my mistake by chance :slight_smile:

4

Hi @rkenko,
I am facing the same error as you.

Can you please point out the mistake you realized?
Thanks in Advance.

5

Hi,

1 - You have to follow this procedure :

2 - If you still face errors, try to see if your mesh is well defined. Typically, with your mesh you have to be able to run this example without no errors :

Good luck and please tell me if it solved your problem :slight_smile:

2 Likes
6

Thank you so much!
This worked finally.

I uploaded my .msh file into the colab link, and everything went by smoothly, except a minor replacement of “gmsh:physical” with “gmsh:geometrical” (which looks like arising from version diffs).

2 Likes
7

Hi,

You are talking about these lines ?

cell_data = numpy.hstack([mesh.cell_data_dict["gmsh:physical"][key]
                         for key in mesh.cell_data_dict["gmsh:physical"].keys() if key==cell_type])

You changed the both ´gmsh:physical’ ? Or which one did you changed ? The first or the second ?

It is strange because it was supposed to work without any changes ^^

8

Yes, I changed both of them
Is it though? It was throwing an error that “gmsh:physical” not found
When I looked into the cell_data_dict, it was named “gmsh:geometrical”, so I changed it.

Anyways, I think there might be some updates because of version diffs.

9

The keys in the cell_data_dict is dependent on how the gmsh mesh has been generated. It is normal to use Physical Line, Physical Surface, Physical Volume to mark domains in gmsh, and they will be found in gmsh:physical. If you do not use Physical Markers gmsh sometimes supplies geometrical markers, which I think is based on how you created each line/surface/volume

10

Hi @dokken,

Did you find the solution to this error :

Unable to find entity in map

I saw many discussions, and one for example with @nschloe , that you had dealing with this issue but no solution on that topics…

12

In your gmsh model, you should add Physical tags, such as Physical Volume, Physical Surface, Physical Curve. With the GMSH python API, this is done with: gmsh.model.addPhysicalGroup
as shown in for instance: Using the GMSH Python API to generate complex meshes | Jørgen S. Dokken

Also, please post your code as a minimal example in the forum, rather than linking to colab. Also please remove all code not required to reproduce your error.

14

Please reduce your problem to the following:

  • mesh generation script
  • script for converting mesh to xdmf
    -script reading in the mesh and plotting it (no solving of a PDE required). Just call plot(mesh)