Gmsh + Meshio fails with mixed elements

1

Hello all,

I am a little stuck with the following code and I haven’t found any useful examples.
I made a 2D Gmsh mesh with some physical groups.
I used the following code to retrieve the facets and the mesh itself:

Wri_path = "./Gmsh_meshes/"
msh = meshio.read(Wri_path+'SinglePBR.msh')

meshio.write(Wri_path+"md_.xdmf",
             meshio.Mesh(points=msh.points,
                         cells={"triangle": msh.cells["triangle"]}))

meshio.write(Wri_path+"mf_.xdmf",
             meshio.Mesh(points=msh.points,
                         cells={"line": msh.cells["line"]},
                         cell_data={"line": {"name_to_read": msh.cell_data["line"]["gmsh:physical"]}}))

# Reading mesh data stored in .xdmf files.
mesh = Mesh()
with XDMFFile(Wri_path+"md_.xdmf") as infile:
    infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 1)

with XDMFFile(Wri_path+"mf_.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

However, the following lines make it fails:

# Define function spaces for PDEs variational formulation.
P1 = FiniteElement('P', triangle, 1)  # Lagrange 1-order polynomials family
element = MixedElement([P1, P1, P1, P1])
V = FunctionSpace(mesh, element)  # Test functions
function_space = FunctionSpace(mesh, P1)

It seems to be that it works just fine when there is only one finite element function, but I haven’t been able to make it runs when there are multiple PDEs. I always obtain the following error message:

File “Reading_Gmsh.py”, line 76, in
V = FunctionSpace(mesh, element)
File “/usr/lib/python3/dist-packages/dolfin/function/functionspace.py”, line 31, in init
self._init_from_ufl(*args, **kwargs)
File “/usr/lib/python3/dist-packages/dolfin/function/functionspace.py”, line 39, in _init_from_ufl
ufl.FunctionSpace.init(self, mesh.ufl_domain(), element)
File “/usr/lib/python3/dist-packages/ufl/functionspace.py”, line 58, in init
error(“Non-matching cell of finite element and domain.”)
File “/usr/lib/python3/dist-packages/ufl/log.py”, line 172, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Non-matching cell of finite element and domain.

I really appreciate your help on this issue or any feedback.

Regards,
Santiago

2

Hello,
try replacing triangle with mesh.ufl_cell().

3

Thank you bleyerj for your fast reply. Your solution solved the issue.

Still, I got this error message:

Dimension mismatch in dot product.
Traceback (most recent call last):
  File "Pht_Anh_PBR_Robin.py", line 212, in <module>
    parameter*u_A*u_B*np.power(u_T, 2*v_A*dx + \
  File "/usr/lib/python3/dist-packages/ufl/operators.py", line 187, in dot
    return Dot(a, b)
  File "/usr/lib/python3/dist-packages/ufl/tensoralgebra.py", line 208, in __new__
    error("Dimension mismatch in dot product.")
  File "/usr/lib/python3/dist-packages/ufl/log.py", line 172, in error
    raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Dimension mismatch in dot product.

I am not sure if it has to be with the way how I am defining the Robin-type boundary conditions in the variational formulation, e.g.:

ds = Measure("ds", domain=mesh, subdomain_data=mf)
v_A, v_B, v_C, v_T = TestFunctions(V)
u = Function(V)
u_n = Function(V)
u_A, u_B, u_C, u_T = split(u)

F_A = ((u_A - u_An)/dt)*v_A*dx + dot(w, grad(u_A))*v_A*dx + \
           D*dot(grad(u_A), grad(v_A))*dx - \
           parameter*u_A*u_B*np.power(u_T, 2*v_A*dx + \
           Vz*(u_A - CAin)*v_A*ds(1)

...

Is it the appropriate way? or do you have any idea?

Thank you once again.

Regards,
Santiago

4

I would guess from the error message that you have non-matching entities in the dot product. In your case, most likely, it is in dot(w, grad(u_A)). In any case, it is difficult to pin down the error without a MWE that can be run to reproduce your error ? Also, I would use ufl’s pow function instead of numpy’s np.power

5

Hello @bhaveshshrimali,

Thank you for your reply. The code was running without a problem before I incorporated the Gmsh mesh with physical entities and meshio abstractions. So despite the error message, I guess that it has to be with the meshio-FEniCS communication or my definition of boundary conditions.

Unfortunately, I couldn’t reproduce it with a short example, so I just created the following repository. It contains the source code and the Gmsh mesh, so you (or anybody) may check it out and reproduce the error as it is.

Nevertheless, I believe the code has good readability, and it is not extense. I would be so grateful if you take some minutes to take a look. Hopefully, it is just a silly mistake and I’m stuck with it.

Best regards,
Santiago

6

With the actual mesh replaced by a template mesh, say UnitSquareMesh(10, 10), your code runs fine.
What version of meshio are you using? I presume it is a slightly older release (guessing from the syntax).

In any case, can you commit your xdmf and h5 files too (Also please try the coarsest possible mesh, I believe the file is too large right now to be rendered properly on git), that makes fetching those directly easier. I could take a look.

Also, to logically join paths, try something like

from dolfin import *
import os
mesh = Mesh()
XDMFFile(os.path.join(Wri_path,"md_.xdmf")).read(mesh) # or a `with`-`as` statement
7

@bhaveshshrimali,

Sorry for the delay, I had some annoying troubles with h5py while updating meshio.
You were right, I was using version 3.1.1 so I updated it to 4.0.11

I know I could use mshr, mark the domain as traditionally, and not mess with meshio. Nevertheless, the idea of using meshio and Gmsh is to be able to mark the domain in the future with much more complex meshes, so this is the first approach.

I committed the following:

  1. New syntax considering new meshio version 4.0.11. I am not sure if I did it correctly, but it seems to be ok, maybe you can tell. I am also using os to logically join paths as suggested.
  2. A much more coarse mesh (msh extension).
  3. The respective xdmf and h5. But remember that the meshio commands call for the msh file to retrieve the physical and geometric entities.

In the end I got stuck again with the same error message.

Thank you once again Bhavesh for following up.

Public repository: