Hello there, it’s been a while I have been using FEniCS(x) from time to time, and I have some old codes and therefore ways of doing things, which I have no idea whether they are well adapted to today’s code and way of doing things.
My goal for now is to solve a simple equation (Ohm’s law) in a simple 3D shape (a bar). I have created a file.msh using the Gmsh GUI, this produced the following file.geo:
//+
SetFactory("OpenCASCADE");
Box(1) = {0, 0, 0, 1e-2, 1e-2, 1};
//+
Physical Volume("the_wire", 13) = {1};
//+
Physical Surface("left_end", 14) = {5};
//+
Physical Surface("right_end", 15) = {6};
and a corresponding file.msh. As far as I understand, nowadays, it is preferred to use pygmsh, but I’m not sure how to reproduce the very same file.msh.
Now, the next step is to use FEniCSx on the file.msh. My attempt is the following:
import meshio
# 1: sample volume
# surfaces
# 5 left end
# 6 right end
msh = meshio.read("file.msh")
So far so good, “msh” yields:
<meshio mesh object>
Number of points: 47
Number of cells:
triangle: 4
triangle: 4
tetra: 71
Cell sets: left_end, right_end, the_wire, gmsh:bounding_entities
Point data: gmsh:dim_tags
Cell data: gmsh:physical, gmsh:geometrical
Field data: left_end, right_end, the_wire
I then proceed with:
meshio.write("mesh.xdmf", meshio.Mesh(points=msh.points, cells={"tetra": msh.cells["tetra"]}))
meshio.write("mf.xdmf", meshio.Mesh(points=msh.points, cells={"triangle": msh.cells["triangle"]},
cell_data={"triangle": {"name_to_read": msh.cell_data["triangle"]["gmsh:physical"]}}))
meshio.write("cf.xdmf", meshio.Mesh(
points=msh.points, cells={"tetra": msh.cells["tetra"]},
cell_data={"tetra": {"name_to_read":
msh.cell_data["tetra"]["gmsh:physical"]}}))
mesh = Mesh()
but this yielded:
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Input In [5], in <cell line: 1>()
----> 1 meshio.write("mesh.xdmf", meshio.Mesh(points=msh.points, cells={"tetra": msh.cells["tetra"]}))
2 meshio.write("mf.xdmf", meshio.Mesh(points=msh.points, cells={"triangle": msh.cells["triangle"]},
3 cell_data={"triangle": {"name_to_read": msh.cell_data["triangle"]["gmsh:physical"]}}))
4 meshio.write("cf.xdmf", meshio.Mesh(
5 points=msh.points, cells={"tetra": msh.cells["tetra"]},
6 cell_data={"tetra": {"name_to_read":
7 msh.cell_data["tetra"]["gmsh:physical"]}}))
TypeError: list indices must be integers or slices, not str
Googling the error yielded me to try:
msh = meshio.read("file.msh")
for key in msh.cell_data_dict["gmsh:physical"].keys():
if key == "triangle":
triangle_data = msh.cell_data_dict["gmsh:physical"][key]
elif key == "tetra":
tetra_data = msh.cell_data_dict["gmsh:physical"][key]
for cell in msh.cells:
if cell.type == "tetra":
tetra_cells = cell.data
elif cell.type == "triangle":
triangle_cells = cell.data
tetra_mesh = meshio.Mesh(points=msh.points, cells={"tetra": tetra_cells},
cell_data={"name_to_read":[tetra_data]})
triangle_mesh =meshio.Mesh(points=msh.points,
cells=[("triangle", triangle_cells)],
cell_data={"name_to_read":[triangle_data]})
meshio.write("plate.xdmf", tetra_mesh)
meshio.write("mf.xdmf", triangle_mesh)
but this yields:
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Input In [9], in <cell line: 14>()
11 triangle_cells = cell.data
12 tetra_mesh = meshio.Mesh(points=msh.points, cells={"tetra": tetra_cells},
13 cell_data={"name_to_read":[tetra_data]})
---> 14 triangle_mesh =meshio.Mesh(points=msh.points,
15 cells=[("triangle", triangle_cells)],
16 cell_data={"name_to_read":[triangle_data]})
17 meshio.write("plate.xdmf", tetra_mesh)
18 meshio.write("mf.xdmf", triangle_mesh)
File /usr/local/lib/python3.9/dist-packages/meshio/_mesh.py:182, in Mesh.__init__(self, points, cells, point_data, cell_data, field_data, point_sets, cell_sets, gmsh_periodic, info)
180 data[k] = np.asarray(data[k])
181 if len(data[k]) != len(self.cells[k]):
--> 182 raise ValueError(
183 "Incompatible cell data. "
184 + f"Cell block {k} ('{self.cells[k].type}') "
185 + f"has length {len(self.cells[k])}, but "
186 + f"corresponding cell data item has length {len(data[k])}."
187 )
ValueError: Incompatible cell data. Cell block 0 ('triangle') has length 4, but corresponding cell data item has length 8.
Googling the error yielded me to try:
import meshio
def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
out_mesh = meshio.Mesh(points=mesh.points, cells={
cell_type: cells}, cell_data={"name_to_read": [cell_data]})
if prune_z:
out_mesh.prune_z_0()
return out_mesh
msh = meshio.read("file.msh")
facet_mesh = create_mesh(msh, "triangle", prune_z=True)
meshio.write("facet_mesh.xdmf", facet_mesh)
triangle_mesh = create_mesh(msh, "tetra", prune_z=True)
meshio.write("mesh.xdmf", triangle_mesh)
but this yielded the error:
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
Input In [15], in <cell line: 15>()
11 return out_mesh
14 msh = meshio.read("file.msh")
---> 15 facet_mesh = create_mesh(msh, "triangle", prune_z=True)
16 meshio.write("facet_mesh.xdmf", facet_mesh)
18 triangle_mesh = create_mesh(msh, "tetra", prune_z=True)
Input In [15], in create_mesh(mesh, cell_type, prune_z)
7 out_mesh = meshio.Mesh(points=mesh.points, cells={
8 cell_type: cells}, cell_data={"name_to_read": [cell_data]})
9 if prune_z:
---> 10 out_mesh.prune_z_0()
11 return out_mesh
AttributeError: 'Mesh' object has no attribute 'prune_z_0'
I ran out of ideas on how to achieve my goal. Any pointer is welcome, and if there’s a quicker/neater way of doing it with pygmsh, then let me know, I’d rather learn good habits than using an archaic way of doing things.
