Error generating quadratic mesh

mesh
1

Hello,

I am trying to generate a quadratic mesh in .msh format and convert it to .xdmf format.
On Gmsh, I used my .geo file which worked fine for the first order, but on the gmsh GUI I additionally clicked on “Set order 2” .
Then I also reused the code that worked fine for the linear mesh, which is the following:

msh = meshio.read("rectangle2D_order2.msh")

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)
    points = mesh.points[:,:2] if prune_z else mesh.points
    out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    return out_mesh

line_mesh = create_mesh(msh, "line", prune_z=True)
meshio.write("rectangle2D_order2_facet_mesh.xdmf", line_mesh)

triangle_mesh = create_mesh(msh, "triangle", prune_z=True)
meshio.write("rectangle2D_order2_mesh.xdmf", triangle_mesh)

However, I get the following error:


---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
/tmp/ipykernel_25068/3165531181.py in <module>
      1 # Convert mesh to XDMF
----> 2 line_mesh = create_mesh(msh, "line", prune_z=True)
      3 meshio.write("rectangle2D_order2_facet_mesh.xdmf", line_mesh)
      4 
      5 triangle_mesh = create_mesh(msh, "triangle", prune_z=True)

/tmp/ipykernel_25068/1200146592.py in create_mesh(mesh, cell_type, prune_z)
      1 # Helper function to extract the cells and physical data
      2 def create_mesh(mesh, cell_type, prune_z=False):
----> 3     cells = mesh.get_cells_type(cell_type)
      4     cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
      5     points = mesh.points[:,:2] if prune_z else mesh.points

~/anaconda3/lib/python3.8/site-packages/meshio/_mesh.py in get_cells_type(self, cell_type)
    122 
    123     def get_cells_type(self, cell_type):
--> 124         return numpy.concatenate([c.data for c in self.cells if c.type == cell_type])
    125 
    126     def get_cell_data(self, name, cell_type):

<__array_function__ internals> in concatenate(*args, **kwargs)

ValueError: need at least one array to concatenate

I am having trouble understanding why this is generating an error, do you have any suggestions on how to solve the problem?

Thank you very much for your help !

2

See: Gmsh import error/ Keyword "tetra" - #2 by dokken

You need to use “triangle6”.
You should check the cell types by printing msh.cell_data_dict to check the cell types.

1 Like
3

Thank you very much @dokken, it worked perfectly!

I have another question that follows from this. If now I want to find the cell size of my quadratic elements, can I reuse the following code (Evaluation of cell sizes)?

tdim = mesh.topology.dim
num_cells = mesh.topology.index_map(tdim).size_local
h = dolfinx.cpp.mesh.h(mesh, tdim, range(num_cells))
print(h)

Because when I run it, it gives the following error:

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Input In [14], in <cell line: 3>()
      1 tdim = mesh.topology.dim
      2 num_cells = mesh.topology.index_map(tdim).size_local
----> 3 h = dolfinx.cpp.mesh.h(mesh, tdim, range(2*num_cells))
      4 h_av = np.mean(h)

RuntimeError: Incompatible dimension of arrays, compile in DEBUG for more info

Or is there a way to convert the quadratic mesh into a linear mesh directly from fenics, without having to re-generate a mesh on gmsh for example?

Thanks again for all your help!

4

You are not using to code you are showing here, as your code says:

while your code says

You have you added a multiple by 2 here? It does not make sense as you do not have that many cells.

5

Yes indeed I copied the error from another test I had done, sorry. I can’t edit my previous post so I’m reposting the code and the error below:

tdim = mesh.topology.dim
num_cells = mesh.topology.index_map(tdim).size_local
h = dolfinx.cpp.mesh.h(mesh, tdim, range(num_cells))
h_av = np.mean(h)

Error:

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Input In [30], in <cell line: 3>()
      1 tdim = mesh.topology.dim
      2 num_cells = mesh.topology.index_map(tdim).size_local
----> 3 h = dolfinx.cpp.mesh.h(mesh, tdim, range(num_cells))
      4 h_av = np.mean(h)

RuntimeError: Incompatible dimension of arrays, compile in DEBUG for more info

The problem must be related to the fact that the mesh is quadratic but I don’t really see how to fix it.

Thank you again!

6

I’ve pushed a fix to:

2 Likes
7

Thank you for your help!
But I still have the same error, should I pull again the Docker image? I thought the latest version of fenics was automatically updated every night.

8

The code has not been merged into main yet, and thus are not in the docker image yet.