Store function with meshio

[nschloe]

Coming back to FEniCS after a while, I forgot how to store mesh functions with meshio. Here’s what I tried, but the order of the points is messed up:

import numpy as np
import meshio
import meshzoo
from dolfin import FunctionSpace, Mesh, MeshEditor, Expression, project

points, cells = meshzoo.disk(6, 10)

editor = MeshEditor()
mesh = Mesh()
# topological and geometrical dimension 2
editor.open(mesh, "triangle", 2, 2, 1)
editor.init_vertices(len(points))
editor.init_cells(len(cells))
for k, point in enumerate(points):
    editor.add_vertex(k, point[:2])
for k, cell in enumerate(cells.astype(np.uintp)):
    editor.add_cell(k, cell)
editor.close()
V = FunctionSpace(mesh, "CG", 1)

f1 = Expression("x[0]", element=V.ufl_element())
e1 = project(f1, V).vector().get_local()

meshio.Mesh(
    V.mesh().coordinates(),
    {"triangle": V.mesh().cells()},
    point_data={"e1": e1},
).write("out.vtk")

If there’s a better way of creating a mesh from points and cells, let me know, too!

[dokken]

I would suggest saving your data to xdmf, as explained here: https://jsdokken.com/converted_files/tutorial_pygmsh.html#second

I’ve also made a similar version coming from cell-sets (pygmsh)

Subdomain definition based on input file

How about: import meshio import numpy as np mesh = meshio.read("mesh.inp") cells = mesh.get_cells_type("triangle") cell_sets = mesh.cell_sets_dict cell_data = np.zeros(len(cells)) for marker, set in cell_sets.items(): for type, entities in set.items(): if type == "triangle": # Create marker with int from integer from name of input cell_data[entities] = int(marker[-1]) out_mesh = meshio.Mesh(points=mesh.points, cells={ "triangle": cells…

and

Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio

I had a play with pygmsh over the weekend, and what seems to be done in the later versions is to save the physical markers to cell sets. I have made the following function for this: def create_entity_mesh(mesh, cell_type, prune_z=False, remove_unused_points=False): """ Given a meshio mesh, extract mesh and physical markers for a given entity. We assume that all unused points are at the end of the mesh.points (this happens when we use physical markers with…

[nschloe]

That doesn’t show how to store point_data, does it?

[dokken]

Sorry Nico, I was a bit to fast reading through your post. How about using the vertex_to_dof_map?

from dolfin import *
import meshio
mesh = UnitSquareMesh(5, 5)
V = FunctionSpace(mesh, "CG", 1)

f1 = Expression("x[0]", element=V.ufl_element())
e1 = project(f1, V).vector().get_local()
v_to_d = vertex_to_dof_map(V)
# vertex_to_dof_map
meshio.Mesh(
    V.mesh().coordinates(),
    {"triangle": V.mesh().cells()},
    point_data={"e1": e1[v_to_d]},
).write("out.vtk")

For higher order spaces where there is no 1-to-1 map, I would use:

dolfin.Function.compute_vertex_values()

which would interpolate to vertices of the mesh.

[nschloe]

@dokken

vertex_to_dof_map

That’s what I was looking for. Thanks!