3D mesh with subdomains, xdmf not working?


272
views
0
5 months ago by
Hi!

I need 3D meshes with subdomains and mshr does not have such a thing yet. In search for a (simple) way of creating meshes I looked into using gmsh through pygmsh, which seems manageable to use. However, I run into problems when writing the mesh with meshio and reading it again with fenics (2017.2.0). The test I used is this:
from dolfin import *
import pygmsh as pg
import meshio

geom = pg.built_in.Geometry()
poly = geom.add_polygon([
    [ 0.0,  0.5, 0.0],
    [-0.1,  0.1, 0.0],
    [-0.5,  0.0, 0.0],
    [-0.1, -0.1, 0.0],
    [ 0.0, -0.5, 0.0],
    [ 0.1, -0.1, 0.0],
    [ 0.5,  0.0, 0.0],
    [ 0.1,  0.1, 0.0]
    ],
    lcar=0.05
    )

axis = [0, 0, 1]

geom.extrude(
    poly,
    translation_axis=axis,
    rotation_axis=axis,
    point_on_axis=[0, 0, 0],
    angle=2.0 / 6.0 * 3.1415926
    )

points, cells, point_data, cell_data, field_data =    pg.generate_mesh(geom,verbose=False, num_lloyd_steps=0)
meshio.write("mesh.xdmf", points, cells, cell_data=cell_data,point_data=point_data,field_data=field_data)

mesh=Mesh()
f = XDMFFile(mpi_comm_world(),"mesh.xdmf")
f.read(mesh)
subdomains = MeshFunction("size_t", mesh,2)
dx = dx(domain=mesh, subdomain_data=subdomains)
print(assemble(1*dx))
​

but I get the error
Traceback (most recent call last):
  File "test_mesh.py", line 35, in <module>
    f.read(mesh)
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to recognise cell type.
*** Reason:  Unknown value "".
*** Where:   This error was encountered inside XDMFFile.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2017.2.0
*** Git changeset:  0baf73825079a581e43ab1705370043040aa213d
*** -------------------------------------------------------------------------​


Can someone here see the mistake in the above code?
Is XDMF not the fileformat to use for meshes with fenics?
Is there a better way of creating simple 3D meshes with subdomains that can easily be parametrized?

It is too bad that mshr does not support 3D meshes with subdomains as it is really a very nice and simple interface to meshing :(

Best Regards,
Søren

Community: FEniCS Project

1 Answer


0
5 months ago by
In case anyone is interested: I gave up on xdmf and wrote in gmsh format instead with meshio.
points, cells, point_data, cell_data, field_data = pg.generate_mesh(geom,verbose=False)
if MPI.rank(mpi_comm_world())==0:
   meshio.write(meshfile.split(".")[0]+".msh", points, cells, cell_data=cell_data,field_data=field_data,file_format='gmsh-ascii')
MPI.barrier(mpi_comm_world())
​

Using dolfin-convert, I convert to xml which I read in and save in h5 format. I use this function for converting from xml to h5:
def convertMesh2h5(meshfile):
   from dolfin_utils.meshconvert import meshconvert
   meshconvert.convert2xml(meshfile.split(".")[0]+".msh", meshfile.split(".")[0]+".xml", iformat='gmsh')
   
   if MPI.rank(mpi_comm_world())==0:
      mesh = Mesh(mpi_comm_self(), meshfile.split(".")[0]+".xml")
      mesh_file = HDF5File(mpi_comm_self(), meshfile.split(".")[0]+".h5", 'w')
      mesh_file.write(mesh, '/mesh')
      subdomains = MeshFunction("size_t", mesh, meshfile.split(".")[0]+"_physical_region.xml")
      boundaries = MeshFunction("size_t", mesh, meshfile.split(".")[0]+"_facet_region.xml")
      mesh_file.write(subdomains, "/subdomains")
      mesh_file.write(boundaries, "/boundaries")
   MPI.barrier(mpi_comm_world())
​

And to read the mesh again

         self.mesh=Mesh()
         f = HDF5File(mpi_comm_world(), meshfile, 'r')
         f.read(self.mesh,"/mesh",False)
         self.subdomains = MeshFunction("size_t", self.mesh)
         self.boundaries = MeshFunction("size_t", self.mesh)
         f.read(self.subdomains, "/subdomains")
         f.read(self.boundaries, "/boundaries")


Seems like a big hazzle to generate a mesh for parallel use, but it seems to allow the flexibility of pygmsh.

Please login to add an answer/comment or follow this question.

Similar posts:
Search »