This is my header of the msh file:
$MeshFormat
2.2 0 8
$EndMeshFormat
$Nodes
2191
I’m quite new to GMSH and its formats, so i’m sorry if its somethin one should know.
But by reading your link, i think this part is irrelevant for me?
Thanks for your replies anyways!
From the header I can tell that you do not have any physical groups marked in your mesh.
Besides the fact that you only have tetrahedra (3D mesh), the cell_data dictionary
will be empty. Even if you used ‘tetra’ instead of ‘triangle’.
Actually i changed "triangle" to "tetra" and i did not get any error messages. But since this is now working i’m quite happy anyways.
But does anyone of you know how to impose boundary conditions on a complex geometry, especially if i want to implement only Neumann and Robin boundary conditions on different parts of my mesh? Are there any examples where i can maybe look it up?
I’m thinking of something, if i f.e. split my mesh in four parts, therefore i can write
boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
ds = Measure("ds", domain=mesh, subdomain_data=boundary_markers)
and evaluate each boundary condition on
ds(1), ds(2), ds(3), ds(4)
Thank you anyway for your help, i highly appreciate it!
That is how it works.
You can use for example gmsh or other CAD software to mark the boundaries and subdomains and then convert the mesh. This would give you 2 mesh functions - 1 for the subdomains, 1 for the boundaries.
Hello, I am facing some issues accessing subdomains when using the XDMF format. When reading a 3D mesh from GMSH with two distinct subdomains, I followed the steps mentioned earlier in the post:
However,
subdomains = MeshFunction("size_t", mesh, mesh.topology.dim, 0)
gives me only one single domain (i.e. all the elements have tag 0) when I query:
print(subdomains.where_equal(0))
(perhaps I am reseting the tag created by GMsh…)
When I open the xdmf mesh in paraview, I can see that there are two distinct subdomains. What’s the correct way of accessing the subdomains when using xdmf meshes (I did not have trouble when using the legacy xml format). For eg I run into errors if I write
a = inner(u, v)*dx(1) + inner(u,v)*dx(2)
Any help is highly appreciated!!
You Need to load the cell-function into dolfin, using
mvc = MeshValueCollection("size_t", mesh, 3)
with XDMFFile("cf.xdmf") as infile:
infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
dx_cf = Measure(“dx”, subdomain_data=mf)```Sorry forgot to mention before that I am using the complex development version for fenics i.e., dolfin-x. So apparently I need to use read_mvc_int function to read the subdomain data but the syntax isn’t clear to me. As you suggested, when I try
mvc = MeshValueCollection("size_t", mesh, 3)
with XDMFFile(MPI.comm_world,"cf.xdmf") as infile:
infile.read_mvc_int(mvc, "name_to_read")
mf = dolfin.cpp.mesh.MeshFunctionSizet(mesh, mvc, 0)
dx_cf = Measure("dx", subdomain_data=mf)
I get:
File “test.py”, line 139, in
infile.read_mvc_int(mvc, “name_to_read”)
File “/usr/local/lib/python3.6/dist-packages/dolfin/io.py”, line 176, in read_mvc_int
return self._cpp_object.read_mvc_int(mesh, name)
TypeError: read_mvc_int(): incompatible function arguments. The following argument types are supported:
1. (self: dolfin.cpp.io.XDMFFile, mesh: dolfin.cpp.mesh.Mesh, name: str=’’) -> dolfin.cpp.mesh.MeshValueCollection_int
Invoked with: <dolfin.cpp.io.XDMFFile object at 0x7f82e822c340>, <dolfin.cpp.mesh.MeshValueCollection_sizet object at 0x7f82e822c378>, ‘name_to_read’
It is not clear to me what should be second argument i.e. name: str=’’ for dolfin to be able to read the subdomain labels. If needed i can post a MWE if that will be of any help. Thanks in advance!!
What I want to finally achieve is to be able to write the weak form by integrating over subdomains, i.e.
a = inner(u,v)*dx(1) + inner(u,v)*dx(2)
Here is the MWE. The “test.msh” could be any 3D mesh created in GMSh with 2 physical volumes.
## dolfin imports
import dolfin
from dolfin import (MPI, CellType, Expression, DirichletBC, Function, FunctionSpace,
RectangleMesh, TestFunction, TrialFunction, solve, MeshFunction,
FacetNormal, dx, ds, dot, grad, inner, CellVolume, MeshValueCollection,
Measure)
from dolfin.fem.assembling import assemble_system
from dolfin.la import PETScKrylovSolver, PETScOptions, PETScVector
from dolfin.io import XDMFFile
from dolfin.cpp.mesh import Mesh, GhostMode, CellType
# other imports
import numpy as np
import os
import sys
import subprocess
import meshio
# read the 3D GMsh file:
msh = meshio.read('mesh.msh')
# write the gmsh file in xdmf format
meshio.write("mesh.xdmf", meshio.Mesh(
points=msh.points,
cells={"tetra": msh.cells["tetra"]}))
# write surface data:
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"]}}))
# write volume data:
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(
MPI.comm_world, CellType.Type.tetrahedron, msh.points[:, 0:3],
msh.cells["tetra"], [], GhostMode.none)
mvc = MeshValueCollection("size_t", mesh, 3)
with XDMFFile(MPI.comm_world,"cf.xdmf") as infile:
infile.read_mvc_int(mesh, "name_to_read") #!!! goes through but does not return subdomain labels in mvc
#infile.read_mvc_int(mvc, "name_to_read") #!!! doesn't work, ends up in incompatible function arguments.
mf = dolfin.cpp.mesh.MeshFunctionSizet(mesh, mvc, 0)
dx_cf = Measure("dx", subdomain_data=mf)
## and finally we want integrate over individual subdomains:
cmap = dolfin.fem.create_coordinate_map(mesh.ufl_domain())
mesh.geometry.coord_mapping = cmap
boundaries = dolfin.MeshFunction("size_t", mesh, mesh.topology.dim - 1, 0)
# define fs:
V = FunctionSpace(mesh, ("Lagrange", 1))
u = TrialFunction(V)
v = TestFunction(V)
a = inner(u,v)*dx ##!!! works
#a = inner(u,v)*dx(1) + inner(u,v)*dx(2) ##!!! does not work
g = 1.0 ## let's compute l2-projection of constant
L = inner(g, v)*dx ##!!! works
#L = inner(g, v)*dx(1) + inner(g, v)*dx(2) ##!!! does not work
A, b = assemble_system(a, L, [])
solver = PETScKrylovSolver(mesh.mpi_comm())
PETScOptions.set("ksp_type", "preonly")
PETScOptions.set("pc_type", "lu")
solver.set_from_options()
solver.set_operator(A)
x = PETScVector()
n_iter = solver.solve(x, b)
u = Function(V)
u.vector().vec().array = x.vec().array
with XDMFFile(MPI.comm_world, "test.xdmf") as xdmf:
xdmf.write(u)Im currently not using dolfinx, but I’m guessing that the command you should use is:
cf = infile.read_mvc_int(mesh, "name_to_read"), as you can see from https://github.com/FEniCS/dolfinx/blob/master/python/test/unit/io/test_XDMF.py Line 80, this command returns a mesh function.
Then, using the syntax:
dx = ufl.Measure('dx', subdomain_data=cf, domain=mesh)
you should be able to do dx(1) and dx(2).
I took slightly different route (perhaps very ad-hoc):
cfvec = infile.read_mf_int(mesh, "name_to_read")
and then
# get element indices for a given domain id:
domain1_ids = cfvec.where_equal(1)
domain2_ids = cfvec.where_equal(2)
# first create a marker and flush it with zeros everywhere:
marker = dolfin.MeshFunction("size_t", mesh, mesh.topology.dim, 0)
# now create an integer array of the size equal to number of elements in mesh
indarray = np.zeros(np.size(marker.array()),'int64')
indarray[domain1_ids] = 1
indarray[domain2_ids] = 2
marker.set_values(indarray)
dx = Measure('dx', subdomain_data=marker, domain=mesh)
But i think there is no way at the moment one could use subdomains in dolfinx. So after defining dx, when i do:
cmap = dolfin.fem.create_coordinate_map(mesh.ufl_domain())
mesh.geometry.coord_mapping = cmap
##### define fs:
V = FunctionSpace(mesh, ("Lagrange", 1))
u = TrialFunction(V)
v = TestFunction(V)
a = inner(u, v) * dx(1) + inner(u, v) * dx(2)
u = Function(V)
g = 1.0 ## let's compute l2-projection of constant
L = inner(g, v) * dx(1) + inner(g, v) * dx(2)
solve(a==L,u,[])
I get:
Traceback (most recent call last):
File “test0.py”, line 62, in
solve(a==L,u,)
File “/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py”, line 113, in solve
return _solve_varproblem(*args, **kwargs)
File “/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py”, line 126, in _solve_varproblem
a = fem.Form(eq.lhs, form_compiler_parameters=form_compiler_parameters)
File “/usr/local/lib/python3.6/dist-packages/dolfin/fem/form.py”, line 64, in init
self._cpp_object = cpp.fem.Form(ufc_form, function_spaces)
RuntimeError: Cell integral subdomain not supported. Under development.
Perhaps it works with a branch other than the master.
I was using the wrong branch i guess.
Using the image from:
docker pull quay.io/fenicsproject/dolfinx:complex
fixed my problem.
Hi I was able to use meshio to read in my files from Gmsh thanks to your instructions here. Currently, my code works fine using the following:
#!/usr/bin/env python3
from fenics import *
import numpy as np
import os
import time
import meshio
#Folder and Directory Names
output_folder_name = "./Kolachalama/"
os.system("rm "+output_folder_name+"*.xdmf")
os.system("rm "+output_folder_name+"*.h5")
output_vtk_folder = "Solution/"
output_vtk_filename = "Kola.pvd"
os.system("rm "+output_folder_name+output_vtk_folder+"*.vtu")
os.system("rm "+output_folder_name+output_vtk_folder+"*.pvd")
os.system("rm "+output_folder_name+output_vtk_folder+"*.pvtu")
vtkfile = File(output_folder_name+output_vtk_folder+output_vtk_filename)
#Define Mesh and Function Space
msh = meshio.read(output_folder_name+"Vessel.msh")
meshio.write(output_folder_name+"Kola.xdmf", meshio.Mesh(points=msh.points, cells{"tetra":msh.cells["tetra"]}))
meshio.write(output_folder_name+"mfKola.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(output_folder_name+"cfKola.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()
with XDMFFile(output_folder_name+"Kola.xdmf") as infile:
infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 2) # Creates a container for mesh data - ie facet data
with XDMFFile(output_folder_name+"mfKola.xdmf") as infile:
infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh,mvc)
boundaries = MeshFunction("size_t", mesh, mvc)
V = FunctionSpace(mesh, 'P', 1)
However, I run into problems when I try to run the code on parallel using:
mpirun -np 6 python3 Kolachalama.py
I get the following error:
Traceback (most recent call last):
File "Kolachalama.py", line 25, in <module>
meshio.write(output_folder_name+"Kola.xdmf", meshio.Mesh(points=msh.points, cells={"tetra":msh.cells["tetra"]}))
File "/home/karthic/.local/lib/python3.6/site-packages/meshio/_helpers.py", line 245, in write
return interface.write(filename, mesh, *args, **_kwargs)
File "/home/karthic/.local/lib/python3.6/site-packages/meshio/_xdmf/main.py", line 501, in write
XdmfWriter(*args, **kwargs)
File "/home/karthic/.local/lib/python3.6/site-packages/meshio/_xdmf/main.py", line 297, in __init__
self.h5_file = h5py.File(self.h5_filename, "w")
File "/home/karthic/.local/lib/python3.6/site-packages/h5py/_hl/files.py", line 394, in __init__
swmr=swmr)
File "/home/karthic/.local/lib/python3.6/site-packages/h5py/_hl/files.py", line 176, in make_fid
fid = h5f.create(name, h5f.ACC_TRUNC, fapl=fapl, fcpl=fcpl)
File "h5py/_objects.pyx", line 54, in h5py._objects.with_phil.wrapper
File "h5py/_objects.pyx", line 55, in h5py._objects.with_phil.wrapper
File "h5py/h5f.pyx", line 105, in h5py.h5f.create
OSError: Unable to create file (unable to flock file, errno = 11, error message = 'Resource temporarily unavailable')
Traceback (most recent call last):
File "Kolachalama.py", line 25, in <module>
meshio.write(output_folder_name+"Kola.xdmf", meshio.Mesh(points=msh.points, cells={"tetra":msh.cells["tetra"]}))
File "/home/karthic/.local/lib/python3.6/site-packages/meshio/_helpers.py", line 245, in write
return interface.write(filename, mesh, *args, **_kwargs)
File "/home/karthic/.local/lib/python3.6/site-packages/meshio/_xdmf/main.py", line 501, in write
XdmfWriter(*args, **kwargs)
File "/home/karthic/.local/lib/python3.6/site-packages/meshio/_xdmf/main.py", line 297, in __init__
self.h5_file = h5py.File(self.h5_filename, "w")
File "/home/karthic/.local/lib/python3.6/site-packages/h5py/_hl/files.py", line 394, in __init__
swmr=swmr)
File "/home/karthic/.local/lib/python3.6/site-packages/h5py/_hl/files.py", line 176, in make_fid
fid = h5f.create(name, h5f.ACC_TRUNC, fapl=fapl, fcpl=fcpl)
File "h5py/_objects.pyx", line 54, in h5py._objects.with_phil.wrapper
File "h5py/_objects.pyx", line 55, in h5py._objects.with_phil.wrapper
File "h5py/h5f.pyx", line 105, in h5py.h5f.create
OSError: Unable to create file (unable to flock file, errno = 11, error message = 'Resource temporarily unavailable')
I’m assuming the error has to do with each process attempting to write the same xdmf file? I’m not too sure how to resolve this. Any assistance would be greatly appreciated.
Split your program in two, run meshio in serial (saving tve xdmf), then run dolfin and loading the mesh in parallel .
@dokken
Hi, I adapted your snippet to read a 2D rectangular channel mesh (triangle elements) produced in gmsh, with named and numbered boundaries (inlet:1, outlet:2, bottom:3, top:4). It’s written in 2.2 format and all boundaries are correctly identified in the msh file. I then import the mesh using the following adaptation of your snippet so the inlet and outlet boundaries are imported separately:
import meshio
msh = meshio.read("channel2Dtri.msh")
meshio.write("meshtriangles.xdmf",
meshio.Mesh(points=msh.points,
cells={"triangle": msh.cells["triangle"]}))
meshio.write("meshlines-inlet.xdmf",
meshio.Mesh(points=msh.points,
cells={"line": msh.cells["line"]},
cell_data={"line": {"inlet": msh.cell_data["line"]["gmsh:physical"]}}))
meshio.write("meshlines-outlet.xdmf",
meshio.Mesh(points=msh.points,
cells={"line": msh.cells["line"]},
cell_data={"line": {"outlet": msh.cell_data["line"]["gmsh:physical"]}}))
from dolfin import *
mesh = Mesh()
with XDMFFile("meshtriangles.xdmf") as infile:
infile.read(mesh)
mvc_i = MeshValueCollection("size_t", mesh, 2)
mvc_o = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("meshlines-inlet.xdmf") as infile:
infile.read(mvc_i, "inlet")
with XDMFFile("meshlines-outlet.xdmf") as infile:
infile.read(mvc_o, "outlet")
mf_i = cpp.mesh.MeshFunctionSizet(mesh, mvc_i)
mf_o = cpp.mesh.MeshFunctionSizet(mesh, mvc_o)
File("dolfinmesh.pvd").write(mesh)
File("dolfinmesh_lines_inlet.pvd").write(mf_i)
File("dolfinmesh_lines_outlet.pvd").write(mf_o)
The code executes fine without any errors, but when I check the xdmf files in paraview meshlines-inlet.xdmf and meshlines-outlet.xdmf show the entire mesh lines, colored by a function with a value of 1 on boundary lines (all ds) and 1.8e+19 on internal lines (all dS), shown in the figure below. Clearly I’m not doing it right, or is it possible at all with meshio and XDMFFILE to import different boundaries separately?
First, you do not need to save one file for each boundary. The gmsh:physical will contain all boundaries.
meshio.write("boundary_markers.xdmf",
meshio.Mesh(points=msh.points,
cells={"line": msh.cells["line"]},
cell_data={"line": {"inlet": msh.cell_data["line"]["gmsh:physical"]}}))
Secondly, the MeshValueCollection should be initialized as mvc_i = MeshValueCollection("size_t", mesh,1)
Thirdly, try scaling the plots in Paraview from 1 to 4. The interior boundaries are not contained in the mesh value collection, and is therefore not initialized, causing the output to be arbitrary large for interior cells.
Thanks for clarifying. I corrected my code and ended up with pretty much your original code snippet. The boundaries are indeed identified with the correct numbers in paraview after scaling 1-4.
When implementing Dirichlet BCs in FEniCS, should I just use the boundary markers 1-4 directly or are there additional steps with the MeshValueCollection?
Since you have loaded the meshvaluecollection into a meshfunction, it can be used as a meshfunction is usually used .
Yes, just found out that DirichletBC can take the meshfunction as an input which allows using the boundary markets directly. I wans’t familiar with that form.
Thanks…
Actually, the topology dimension argument is ignored.
Assume, the xdmf-file is created as follows (2D):
outboundary = meshio.Mesh(points=outpoints,
cells={'line': inmsh.cells['line']},
cell_data={'line': {'Boundary': inmsh.cell_data['line']['gmsh:physical']}},
field_data=inmsh.field_data)
meshio.write(filename=outfile_boundary, mesh=outboundary, file_format=outfile_format)
and afterwards read back in to dolfin as:
mvc = MeshValueCollection('size_t', mesh) # mvc.dim() == 18446744073709551615 (-1 signed int)
with XDMFFile(infile_boundary) as infile:
infile.read(mvc, 'Boundary') # mvc.dim() == mesh.topology().dim()
boundaries = cpp.mesh.MeshFunctionSizet(mesh, mvc)
Indeed, XDMFFile.read clears the MeshValueCollection and reinitializes it with the given mesh.
So, after the line infile.read(mvc, 'Boundary'), the topological dimension of the MeshValueCollection is that of the mesh written by meshio in the first step, no matter what you write in the constructor.
So is that why the internal lines (dS) are getting a boundary marker of 1.8e19?
BTW using infile.read(mvc, 'Boundary') throws an error:
*** Error: Unable to recognise cell type.
*** Reason: Unknown value “”.
I’m using DOLFIN version: 2018.1.0