How to use MeshValueCollection?

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

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};

i found my mistake by chance :slight_smile:

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

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

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

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

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 ^^

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.

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

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…

Hi @dokken,
Thank you for this. Yes I see. I had created the mesh using a different method.

One thing, in my colab notebook cell here, I used the code which @rkenko suggested.
It ran smoothly till meshing, but the solution came only for the top surface (instead of the volume of the box)

Can you help me know if there is something missing?

Thanks in Advance!

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.

Yes, I see. I had not added the PhysicalGroups. Thanks for this.!

The problem is that I can see the mesh properly, but the solution gives values only on one surface. Is this because of missing Physical Volume?
Code snippets:

  1. Read the .msh file
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:geometrical"][key]
                         for key in mesh.cell_data_dict["gmsh:geometrical"].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

msh = meshio.read(file_name + ".msh")

line_mesh = create_mesh(msh, "line", prune_z=False)
meshio.write("line_mesh.xdmf", line_mesh)

triangle_mesh = create_mesh(msh, "triangle", prune_z=False)
meshio.write("triangle_mesh.xdmf", triangle_mesh)

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

  1. Apply boundary conditions, define variational expression and iterate the solver for num_steps:
thermal_conductivity = 7.7
T_amb = 300
convec_coeff = 50
tol = 1e-14

T = 5            # final time
num_steps = 10000    # number of time steps
dt = T / num_steps # time step size

speed = 0.005
curr_x = 0
curr_y = 0

#Gaussian laser flux of diameter = 'diameter', centre = (curr_x, curr_y)
laser_flux = Expression('pow(x[0]-curr_x, 2) + pow(x[1]-curr_y, 2) < pow(diameter/2, 2) ? (8*P/(pi*pow(diameter,2)*kappa))*exp(-8*(pow(x[0]-curr_x, 2) + pow(x[1]-curr_y, 2))/pow(diameter,2)) : 0',
                       degree=2,
                       diameter=0.005,
                       curr_x = 0,
                       curr_y = 0,
                       P = 1700,
                       kappa = thermal_conductivity)

r = Expression('h*kappa',
               degree = 2,
               h = convec_coeff,
               kappa = thermal_conductivity)

s = Expression('s',
               degree = 2,
               s = T_amb)

#T_ambience
T_amb_exp = Expression('T_amb',
                       degree = 2,
                       T_amb = T_amb)

f = Constant(0)
V = FunctionSpace(mesh, 'P', 1)

boundary_conditions = {0: {'Neumann': laser_flux},
                       1: {'Robin':     (r, s)}}

class BoundaryTop(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[2], object_dim[2]/2, tol)

class BoundaryBot(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[2], -object_dim[2]/2, tol)

# Mark boundaries
boundary_markers = MeshFunction("size_t", mesh, 2)
boundary_markers.set_all(9999)
bz1 = BoundaryTop()
bz2 = BoundaryBot()

bz1.mark(boundary_markers, 0)
bz2.mark(boundary_markers, 1)

# Redefine boundary integration measure
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)

# Collect Dirichlet conditions
bcs = []
for i in boundary_conditions:
    if 'Dirichlet' in boundary_conditions[i]:
        bc = DirichletBC(V, boundary_conditions[i]['Dirichlet'],
                          boundary_markers, i)
        bcs.append(bc)


# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)

# Collect Neumann integrals
integrals_N = []
for i in boundary_conditions:
    if 'Neumann' in boundary_conditions[i]:
        if boundary_conditions[i]['Neumann'] != 0:
            g = boundary_conditions[i]['Neumann']
            integrals_N.append(g*v*dt*ds(i))

# Simpler Robin integrals
integrals_R = []
for i in boundary_conditions:
    if 'Robin' in boundary_conditions[i]:
        r, s = boundary_conditions[i]['Robin']
        integrals_R.append(r*(u - s)*v*dt*ds(i))

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)

u_old = interpolate(T_amb_exp, V)

# Simpler variational problem
F = (u - u_old)*v*dx + dot(grad(u), grad(v))*dt*dx + \
    sum(integrals_R) - f*v*dt*dx + sum(integrals_N)
a, L = lhs(F), rhs(F)

# Compute solution
u = Function(V)

t = 0
for n in range(num_steps):

    # Update current time
    t += dt
    curr_x += speed*dt
    laser_flux.curr_x = curr_x

    # Compute solution
    solve(a == L, u, bcs)

    # Plot solution
    plot(u)

    # Update previous solution
    u_old.assign(u)

File("Solution.pvd").write(u)
  1. Result of plot(u)

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)