Create mesh in fenicsx from point and element array

12

Question: How to construct meshtags from the pair (facets, values) and not (facet_indices, values)?

I am trying to read mesh data from here and saving it as .xdmf file with meshtag information.


import h5py
import numpy as np
import dolfinx
import ufl
from mpi4py import MPI

filename = "colin27_coarse_boundaries.h5"
file = h5py.File(filename, 'r')

coordinates = np.array(file["mesh/coordinates"])
topology = np.array(file["mesh/topology"])
boundary_coordinates = np.array(file["boundaries/coordinates"])
boundary_facets = np.array(file["boundaries/topology"])
tags = np.array(file["boundaries/values"], dtype=np.int64)

gdim = 3
shape = "tetrahedron"
degree = 1

cell = ufl.Cell(shape, geometric_dimension=gdim)
ufl_tetra = ufl.Mesh(ufl.VectorElement("Lagrange", cell, degree))

mesh = dolfinx.mesh.create_mesh(MPI.COMM_WORLD, topology, coordinates, ufl_tetra)

tdim = mesh.topology.dim
fdim = tdim - 1

facet_tag = dolfinx.mesh.meshtags_from_entities(mesh, fdim, boundary_facets, tags)

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "brain.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(facet_tag)

Thanks