"Reason: Cannot read old-style MeshFunction XML files as a MeshValueCollection." When running script using mpirun

Errors
1

Hi everyone,

Currently I have a Fenics script that imports a mesh, determines the appropriate cells and facets then uses the .xml file of the facets region to determine the surface indices of the facets that lie on a given surface. Currently, it looks like this:

from __future__ import print_function
from dolfin import *
import ufl
from mpi4py import MPI
#---------------------------------------
data = 'crucible'
mesh = Mesh()
hdf = HDF5File(mesh.mpi_comm(), data+'.h5', 'r')
hdf.read(mesh, '/mesh', False)
cells_f = MeshFunction('size_t', mesh, 3)
hdf.read(cells_f, '/cells')
facets_f = MeshFunction('size_t', mesh, 2)
hdf.read(facets_f, '/facets')
boundaries = MeshFunction('size_t', mesh, 'crucible_facet_region.xml') 
#-----------------------------------------------
# Mark the tags for 'top' and 'surface' groups
#-----------------------------------------------
tag_s = 2  #Tags for 'surface' group in .msh file
tag_t = 3
surfacefacetindices = np.where(boundaries.array() == tag_s)[0]
surfacefacetindices_2 = np.where(boundaries.array() == tag_t)[0]
surfacefacetindices_total = np.append(surfacefacetindices,surfacefacetindices_2)
surfacefacets = np.array(list(facets(mesh)))[surfacefacetindices_total]

print(surfacefacets)

However if I try running this code using mpirun, like so:

mpirun -n 20 python3 test_boundaries.py

I get the following error:

  File "test_boundaries.py", line 14, in <module>
    boundaries = MeshFunction('size_t', mesh, 'crucible_facet_region.xml') 
  File "/usr/local/fkf/Fenics_openmpi2/Python3/lib64/python3.6/site-packages/dolfin/mesh/meshfunction.py", line 18, in __new__
    return fn(mesh, dim)
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
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to read mesh function from XML file.
*** Reason:  Cannot read old-style MeshFunction XML files as a MeshValueCollection.
*** Where:   This error was encountered inside XMLMeshFunction.h.
*** Process: 0
*** 
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  b55804ecca7d010ff976967af869571b56364975
*** -------------------------------------------------------------------------

Does anyone know how I can fix this. This code runs fine without the mpirun

The required files can be found here: Crucible - Google Drive

2

In general, xml is not well suited for parallel computations. I would suggest writing the facet data to an xdmf file (running without mpirun), and then use the xdmf files to import your data when running with MPI.
Ref Bitbucket
You could also just try to rewrite your facet data to xml (as it will then use the new-style).

3

Well this code actually reads the facet data from a .h5 file so I have now adjusted the boundaries term such that

boundaries = facet_f

However, this does not provide a list of indices as desired, rather a list of dolfin.Facet objects. Do you have any further suggestions?

4

The line you where asking about reads from an xml file.

If you remove this, as you described

You get a list of dolfin.Facet. You can get the index from each of these facets

 surface_facets = [facet.index() for facet in surfacefacets]

But I dont really understand why you do:

as surfacefacetindices_total are the indices of each facet (local to process) that is tagged with either 2 or 3.

Here is a mwe:

from __future__ import print_function
from dolfin import *
import numpy as np
# ---------------------------------------
data = 'crucible'
mesh = Mesh()
hdf = HDF5File(mesh.mpi_comm(), data+'.h5', 'r')
hdf.read(mesh, '/mesh', False)
cells_f = MeshFunction('size_t', mesh, 3)
hdf.read(cells_f, '/cells')
facets_f = MeshFunction('size_t', mesh, 2)
hdf.read(facets_f, '/facets')

tag_s = 2  # Tags for 'surface' group in .msh file
tag_t = 3
surfacefacetindices = np.where(facets_f.array() == tag_s)[0]
surfacefacetindices_2 = np.where(facets_f.array() == tag_t)[0]
surfacefacetindices_total = np.append(
    surfacefacetindices, surfacefacetindices_2)
print(surfacefacetindices_total)

returning

[    0     1     3 ..., 82925 82934 82940]