Dear Dokken,
Thank you for your reply. It seems the first step has been done successfully since I haven’t received any error. But my question is, is this code works for 3D domains where the elements are tetrahedrons. I am asking because I opened the XDMF files (i.e. facet_mesh.xdmf and mesh.xdmf) in Paraview. It shows me a 2D geometry (i.e. only a cross-section of my geometry), while my geometry is a 3D meshed body.
Thanks
To use it with a 3D geometry, you have to extract 3D data. The code above extracts line segments:
And triangular elements:
Thus you should call:
tetra_mesh = create_mesh(mesh_from_file, "tetra")
Dear Dokken,
Now, I am using the following code:
import meshio
mesh_from_file = meshio.read(“mesh.msh”)
def create_mesh(mesh, cell_type, prune_z=False):
cells = numpy.vstack([cell.data for cell in mesh.cells if cell.type==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})
return mesh_new
line_mesh = create_mesh(mesh_from_file, “line”, prune_z=True)
meshio.write(“facet_mesh.xdmf”, line_mesh)
triangle_mesh = create_mesh(mesh_from_file, “triangle”, prune_z=True)
meshio.write(“mesh.xdmf”, triangle_mesh)
tetra_mesh = create_mesh(mesh_from_file, “tetra”)
meshio.write(“mesh_3d.xdmf”, tetra_mesh)
print(mesh_from_file.cell_data_dict) #To ensure the mesh extraction is done successfully.
However, it does not create the mesh_3d.xdmf file.
Please simplify the example by removing the line_mesh and traingle_mesh, and only consider the
tetra_mesh. If you do print tetra_mesh, what does it show you?
import meshio
mesh_from_file = meshio.read(“mesh.msh”)
def create_mesh(mesh, cell_type, prune_z=False):
cells = numpy.vstack([cell.data for cell in mesh.cells if cell.type==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})
return mesh_new
tetra_mesh = create_mesh(mesh_from_file, “tetra”)
print(tetra_mesh)
meshio.write(“mesh_3d.xdmf”, tetra_mesh)
Here is the message that I am getting:
Traceback (most recent call last):
File “code.py”, line 1, in
from fenics import *
File “/usr/lib/python3/dist-packages/fenics/init.py”, line 7, in
from dolfin import *
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/init.py”, line 138, in
from . import parameter
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/parameter/init.py”, line 11, in
from ffc import default_jit_parameters
File “/usr/lib/python3/dist-packages/ffc/init.py”, line 24, in
from ffc.compiler import compile_form, compile_element
File “/usr/lib/python3/dist-packages/ffc/compiler.py”, line 127, in
from ffc.representation import compute_ir
File “/usr/lib/python3/dist-packages/ffc/representation.py”, line 44, in
from ffc.fiatinterface import create_element, reference_cell
File “/usr/lib/python3/dist-packages/ffc/fiatinterface.py”, line 31, in
import FIAT
File “/usr/lib/python3/dist-packages/FIAT/init.py”, line 8, in
from FIAT.finite_element import FiniteElement, CiarletElement # noqa: F401
File “/usr/lib/python3/dist-packages/FIAT/finite_element.py”, line 13, in
from FIAT.polynomial_set import PolynomialSet
File “/usr/lib/python3/dist-packages/FIAT/polynomial_set.py”, line 19, in
from FIAT import expansions
File “/usr/lib/python3/dist-packages/FIAT/expansions.py”, line 12, in
import sympy
File “/usr/lib/python3/dist-packages/sympy/init.py”, line 57, in
from .core import *
File “/usr/lib/python3/dist-packages/sympy/core/init.py”, line 8, in
from .expr import Expr, AtomicExpr, UnevaluatedExpr
File “/usr/lib/python3/dist-packages/sympy/core/expr.py”, line 6, in
from .evalf import EvalfMixin, pure_complex
File “/usr/lib/python3/dist-packages/sympy/core/evalf.py”, line 28, in
from sympy.utilities.iterables import is_sequence
File “/usr/lib/python3/dist-packages/sympy/utilities/init.py”, line 13, in
from .lambdify import lambdify
File “/usr/lib/python3/dist-packages/sympy/utilities/lambdify.py”, line 19, in
from sympy.utilities.decorator import doctest_depends_on
File “/usr/lib/python3/dist-packages/sympy/utilities/decorator.py”, line 11, in
from sympy.utilities.runtests import DependencyError, SymPyDocTests, PyTestReporter
File “/usr/lib/python3/dist-packages/sympy/utilities/runtests.py”, line 22, in
import pdb
File “/usr/lib/python3.8/pdb.py”, line 77, in
import code
File “/mnt/e/research/fenics_codes/meshing/code.py”, line 23, in
tetra_mesh = create_mesh(mesh_from_file, “tetra”)
File “/mnt/e/research/fenics_codes/meshing/code.py”, line 14, in create_mesh
cells = numpy.vstack([cell.data for cell in mesh.cells if cell.type==cell_type])
File “<array_function internals>”, line 5, in vstack
File “/usr/lib/python3/dist-packages/numpy/core/shape_base.py”, line 282, in vstack
return _nx.concatenate(arrs, 0)
File “<array_function internals>”, line 5, in concatenate
ValueError: need at least one array to concatenate
That means there are no tetrahedral elements in your mesh.
What does
print(mesh_from_file)
show?
Dear Dokken,
Since I am using Ubuntu as the Windows10 subsystem, I am not able to see the print contents.
Please use colab as I showed above:
Then you can see print contents.
I have given you as much assistance as I am able to. If you are still having issue, please make a minimal failing example.
Thank you so much, Dear Dokken, for all your efforts!
It seems, all that is discussed above, implies that the mesh is present in msh22 format, not in msh41 (default in newer gmsh versions).
I’m just noting this here because if it was mentioned I missed it. And it caused me quite some cryptic errors and head scratching…
Create a msh22 formatted mesh using e.g.
gmsh -format msh22 thefile.geo
or with pygmsh/meshio:
# ...
mesh.write("thefile.msh", file_format="gmsh22")
I disagree with your claim. All the code I have posted above is using msh 4.1.
The following illustrates a use case with the latest version of GMSH. Starting from a clean ubuntu 20.04 environment using docker:
docker run -ti -v $PWD:/home/shared -w /home/shared --rm ubuntu:20.04
Installation instructions in docker
apt-get update -qq
apt-get install -y -qq software-properties-common python3-pip libglu1 libxrender1 libcursor1 libxft2 libxinerama1
add-apt-repository -y ppa:fenics-packages/fenics
apt install -y --no-install-recommends fenics
pip3 -q install --upgrade sympy
pip3 install gmsh
export HDF5_MPI="ON"
export CC=mpicc
export HDF5_DIR="/usr/lib/x86_64-linux-gnu/hdf5/openmpi/"
pip3 install --no-binary=h5py h5py meshio
# Download tutorial_t1.geo from https://gitlab.onelab.info/gmsh/gmsh/-/blob/master/tutorial/t1.geo
/usr/local/lib/python3.8/site-packages/gmsh-4.8.0-Linux64-sdk/bin/gmsh -2 tutorial_t1.geo
The mesh, as shown above has the following msh format:
$MeshFormat
4.1 0 8
$EndMeshFormat
$PhysicalNames
1
2 6 "My surface"
$EndPhysicalNames
$Entities
# ....
The following code runs without error:
import numpy
import meshio
mesh_from_file = meshio.read("tutorial_t1.msh")
def create_mesh(mesh, cell_type, prune_z=False):
cells = numpy.vstack(
[cell.data for cell in mesh.cells if cell.type == 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})
return mesh_new
triangle_mesh = create_mesh(mesh_from_file, "triangle")
meshio.write("mesh_3d.xdmf", triangle_mesh)
You are right, your example works.
The problem I encountered seems to be related to pygmsh. I think the MWE* would be a bit off-topic here, though.
*) Minimum failing example, that should be 
Okay, so checking your code and my problem, it seems it is related to the cell data.
Your function create_mesh does not seem to read any cell data, i.e. the physical labels, right?
If I wanted to do that, I’d need a slightly different function (that I think you posted earlier), like this one:
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])
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
But this one only seems to work with msh22 because of the way cell data is stored there. And this is what fails with msh41 because, at least for me, that key is gone (there are others like gmsh:geometrical).
I think I’ll have to understand the changes in the 4.1 mesh format better. And maybe this is even related with some recent changes in meshio. But I am not sure.
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 pygmsh)
"""
cells = mesh.get_cells_type(cell_type)
try:
# If mesh created with gmsh API it is simple to extract entity data
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
except KeyError:
# If mehs created with pygmsh, we need to parse through cell sets and sort the data
cell_entities = []
cell_data = []
cell_sets = mesh.cell_sets_dict
for marker, set in cell_sets.items():
for type, entities in set.items():
if type == cell_type:
cell_entities.append(entities)
cell_data.append(np.full(len(entities), int(marker)))
cell_entities = np.hstack(cell_entities)
sorted = np.argsort(cell_entities)
cell_data = np.hstack(cell_data)[sorted]
if remove_unused_points:
num_vertices = len(np.unique(cells.reshape(-1)))
# We assume that the mesh has been created with physical tags,
# then unused points are at the end of the array
points = mesh.points[:num_vertices]
else:
points = mesh.points
# Create output mesh
out_mesh = meshio.Mesh(points=points, cells={cell_type: cells},
cell_data={"name_to_read": [cell_data]})
if prune_z:
out_mesh.prune_z_0()
return out_mesh
Thank you! I think this is going to be very useful.
I tried something similar but stopped when I found that the cell entities and the cells seemed to have a different order, resulting in mixed up assignments. That’s when I (prematurely) concluded that it’s to do with gmsh 4.1 (and/or meshio, pygmsh).
I think this kind of still holds and I would like to give a short account here in case it helps anyone:
I tried to use pygmsh to generate my mesh with physical labels (“1” and “2” for 2d subdomains, “10” and “11” for boundaries). I’ll post an MWE below to illustrate.
First, I use mesh.write("gmsh_mwe.msh", file_format="gmsh22"). It works and produces a .msh file.
However, converting this to xdmf results in field data name_to_read filled with 1e-38 values, i.e. something is either uninitialised or type-wise lost in translation. I don’t know where to search here and I suspect it’s to do with pygmsh and meshio.
Second, I use mesh.write("gmsh_mwe.msh", file_format="gmsh") (i.e. gmsh 4.1). This results in a WriteError: Specify entity information to deal with more than one cell type which is certainly pygmsh/meshio related. Probably @nschloe can tell me that I’m making a really stupid mistake.
The way I chose now is to unroll the geo file using the gmsh API, then manually call gmsh -format msh22 gmsh_mwe.geo_unrolled and go on from there. Funny enough, you code @dokken works like a breeze with msh files created with gmsh on the command line.
This workflow is a bit of a bummer however for what I wanted to do (use pygmsh to update the geometry programmatically) :-/
MWE: gmsh_mwe.py
#!/usr/bin/python3
import gmsh
import pygmsh as msh
w1, h1 = 10., 10.
w2, h2 = 5., 5.
meshsize = 2.
with msh.geo.Geometry() as geom:
points = [
[0., 0.],
[-w1, 0.],
[-w1, h1],
[0., h1],
[w2, 0.],
[w2, h2],
[0., h2],
]
p = [geom.add_point(P[:2], mesh_size=meshsize) for P in points]
n_pts = len(points)
lines = [
[0, 1],
[1, 2],
[2, 3],
[3, 6],
[6, 5],
[5, 4],
[4, 0],
[6, 0],
]
l = [geom.add_line(p[p1], p[p2]) for p1, p2 in lines]
# Rect 1
print("== Rect 1 "+"="*50)
ls = [l[0], l[1], l[2], l[3], l[7]]
rect1_cl = geom.add_curve_loop(ls)
rect1_s = geom.add_plane_surface(rect1_cl)
geom.add_physical(rect1_s, "1")
# Rect 2
print("== Rect 2 "+"="*50)
ls = [-l[7], l[4], l[5], l[6]]
rect2_cl = geom.add_curve_loop(ls)
rect2_s = geom.add_plane_surface(rect2_cl)
geom.add_physical(rect2_s, "2")
# Make the two boundaries
left = rect1_cl.curves[:-1]
left = geom.add_physical(left, "10")
right = rect2_cl.curves[1:]
right = geom.add_physical(right, "11")
mesh = geom.generate_mesh() # (verbose=True)
name_root = "gmsh_mwe"
# Variant 1: gmsh22
mesh.write(name_root + ".msh", file_format="gmsh22")
# Variant 2: gmsh 4.1
# Results in WriteError: Specify entity information to deal with more than one cell type
mesh.write(name_root + ".msh", file_format="gmsh")
# Variant 3: unroll geo file
gmsh.write(name_root + ".geo_unrolled")
For Variant 3, run:
gmsh -format msh22 gmsh_mwe.geo_unrolled
The mesh using variants 1 or 2 (check legend, 1e-38 ish):
The mesh using variant 3 (manual gmsh):
Probably nschloe can tell me that I’m making a really stupid mistake.
Nah, this is pygmsh related. Best post an MWE at the issue tracker: Issues · nschloe/pygmsh · GitHub.
Hi, I was trying to apply the above code to my own mesh/problem. The mesh I want to work with is test3.msh, which is giving me the error:
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
File "/usr/local/lib/python3.8/dist-packages/meshio/_mesh.py", line 213, in get_cell_data
return np.concatenate(
File "<__array_function__ internals>", line 5, in concatenate
ValueError: need at least one array to concatenate
Since I couldn’t figure out, what the problem is, I started working with a simpler mesh, test4.msh.
I can read this mesh without any problems. However, applying the DirichletBC gives the following answer:
*** Warning: Found no facets matching domain for boundary condition.
Here I also tried renaming “top” to “0” in the mesh file.
So my questions are:
Thank You.
The code and .msh files are listed below. Both are generated in gmsh and exported as .msh-files
from dolfin import cpp
from fenics import *
from vtkplotter.dolfin import plot as dolfin_plot
import meshio
from ufl import nabla_div
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("test4.msh")
print(msh)
triangle_mesh = create_mesh(msh, "triangle", True)
line_mesh = create_mesh(msh, "line", True)
meshio.write("mesh.xdmf", triangle_mesh)
meshio.write("mf.xdmf", line_mesh)
from dolfin import *
mesh = Mesh()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
infile.read(mesh)
infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim() - 1)
with XDMFFile("mf.xdmf") as infile:
infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
#ds = Measure('ds', domain=mesh, subdomain_data=mf)
#ds_custom = Measure("ds", domain=mesh, subdomain_data=mf)
#print(assemble(1*ds_custom(1)), assemble(1*ds_custom(2)))
#dolfin_plot(mf)
#boundaries = MeshFunction("size_t", mesh, dim=0)
#dolfin_plot(boundaries)
# Define function space for system of PDEs
degree = 2
lambda_ = 1
mu = 1
V = VectorFunctionSpace(mesh, 'P', degree)
# Define boundary conditions
tol = 1E-14
def clamped_boundary(x, on_boundary): # beam is only fixed on the left end
return on_boundary and near(x[0], 0, tol)#x[0] < tol
def force_on_rhs_boundary(x, on_boundary):
return on_boundary and near(x[0], 1.0, tol)
# https://fenicsproject.discourse.group/t/transitioning-from-mesh-xml-to-mesh-xdmf-from-dolfin-convert-to-meshio/412/93
# Boundary conditions
bc1 = DirichletBC(V, Constant((0, 1)), mf, 1)
#bc3 = DirichletBC(V, Constant((0, 0)), clamped_boundary)
#bc2 = DirichletBC(V, Constant((1, 0)), force_on_rhs_boundary)
# Combine dirichlet boundary conditions
bcs = [bc1] #, bc2, bc3]
# Define strain and stress
def epsilon(u): # strain
return 0.5*(nabla_grad(u) + nabla_grad(u).T)
# return sym(nabla_grad(u))
def sigma(u): # stress
return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u)
# Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension() # space dimension
v = TestFunction(V)
f = Constant((0, 0))
T = Constant((0, 0))
a = inner(sigma(u), epsilon(v)) * dx
L = dot(f, v) * dx + dot(T, v) * ds
# Compute solution
u = Function(V)
solve(a == L, u, bcs)
# Plot solution
print('sol u')
dolfin_plot(u, title='Displacement', mode='displacement')
test4.msh:
$MeshFormat
4.1 0 8
$EndMeshFormat
$PhysicalNames
1
1 8 "top"
$EndPhysicalNames
$Entities
6 7 2 0
1 0 0 0 0
2 1 0 0 0
3 0 1 0 0
4 1 1 0 0
5 0 0.5 0 0
6 1 0.5 0 0
1 0 0.5 0 0 1 0 1 8 2 3 -5
2 0 1 0 1 1 0 2 3 8 2 3 -4
3 1 0.5 0 1 1 0 2 4 8 2 4 -6
4 0 0.5 0 1 0.5 0 1 8 2 5 -6
5 0 0 0 0 0.5 0 0 2 5 -1
6 1 0 0 1 0.5 0 1 4 2 6 -2
7 0 0 0 1 0 0 0 2 2 -1
1 0 0.5 0 1 1 0 1 1 4 1 4 -3 -2
2 0 0 0 1 0.5 0 1 2 4 7 -5 4 6
$EndEntities
$Nodes
13 8 1 8
0 1 0 1
1
0 0 0
0 2 0 1
2
1 0 0
0 3 0 1
3
0 1 0
0 4 0 1
4
1 1 0
0 5 0 1
5
0 0.5 0
0 6 0 1
6
1 0.5 0
1 1 0 0
1 2 0 0
1 3 0 0
1 4 0 0
1 6 0 0
2 1 0 1
7
0.4999999999999999 0.75 0
2 2 0 1
8
0.4999999999999999 0.25 0
$EndNodes
$Elements
7 13 1 13
1 1 1 1
1 3 5
1 2 1 1
2 3 4
1 3 1 1
3 4 6
1 4 1 1
4 5 6
1 6 1 1
5 6 2
2 1 2 4
6 6 4 7
7 3 5 7
8 4 3 7
9 5 6 7
2 2 2 4
10 6 2 8
11 1 5 8
12 2 1 8
13 5 6 8
$EndElements
test3.msh was too big, I will post it in another answer or give a link.
Edit: test3.msh can be found here:
https://raw.githubusercontent.com/lukasschulth/EBZ/master/coding/Elasticity%20Problems/test3.msh?token=AJJXEY7LSDWXIBCB6AI7JNTARJ4IA
Hi,
I recently started using gmsh, so my answer could be totally wrong, but maybe worth trying.
In your test4.msh file,
There is a line
1 8 ‘’top”
This means that “top” has tag=8, so you should use
bc1 = DirichletBC(V, Constant((0, 1)), mf, 8)
instead of
bc1 = DirichletBC(V, Constant((0, 1)), mf, 1)
The first number indicates the dimension and the second number indicates the tag if I understand correctly.
I do not know why test3h.msh does not work, sorry
Hey, your answer regarding test4.msh is absolutely correct.
Thanks alot!
Hey,
I solved the import problem by exporting the mesh as Version 2 ASCII with the option ‘save all elements’. Before I was always exporting as Version 4 ASCII.
The import into fenics works as I want to, but now I can’t choose the boundaries by their tags anymore as you explained in your post. I’m getting the following error again:
*** Warning: Found no facets matching domain for boundary condition.
What am I missing here?
How can I both import the mesh file and access the boundaries?
The first few lines of the mesh file are:
$MeshFormat
2.2 0 8
$EndMeshFormat
$PhysicalNames
4
1 21 "right"
1 22 "top"
1 23 "left"
1 24 "bot"
$EndPhysicalNames