# Gmsh + Meshio fails with mixed elements

**URL:** <https://fenicsproject.discourse.group/t/gmsh-meshio-fails-with-mixed-elements/2989>\
**Category:** Uncategorized\
**Created:** [April 22, 2020, 7:13am UTC](https://fenicsproject.discourse.group/t/gmsh-meshio-fails-with-mixed-elements/2989 "2020-04-22T07:13:32Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![SantiagoOrtiz](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/santiagoortiz/32/201_2.png) [@SantiagoOrtiz](https://fenicsproject.discourse.group/u/SantiagoOrtiz)\
**Post date:** [April 22, 2020, 7:13am UTC](https://fenicsproject.discourse.group/t/gmsh-meshio-fails-with-mixed-elements/2989/1 "2020-04-22T07:13:32Z")

</div>

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:

```python
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:

```python
# 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

---

<div class="post-metadata">

**Author:** ![bleyerj](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bleyerj/32/6020_2.png) [@bleyerj](https://fenicsproject.discourse.group/u/bleyerj)\
**Post date:** [April 22, 2020, 7:39am UTC](https://fenicsproject.discourse.group/t/gmsh-meshio-fails-with-mixed-elements/2989/2 "2020-04-22T07:39:15Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![SantiagoOrtiz](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/santiagoortiz/32/201_2.png) [@SantiagoOrtiz](https://fenicsproject.discourse.group/u/SantiagoOrtiz)\
**Post date:** [April 22, 2020, 3:12pm UTC](https://fenicsproject.discourse.group/t/gmsh-meshio-fails-with-mixed-elements/2989/3 "2020-04-22T15:12:23Z")

</div>

Thank you [bleyerj](https://fenicsproject.discourse.group/u/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.:

```python
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

---

<div class="post-metadata">

**Author:** ![bhaveshshrimali](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bhaveshshrimali/32/124_2.png) [@bhaveshshrimali](https://fenicsproject.discourse.group/u/bhaveshshrimali)\
**Post date:** [April 22, 2020, 7:14pm UTC](https://fenicsproject.discourse.group/t/gmsh-meshio-fails-with-mixed-elements/2989/4 "2020-04-22T19:14:22Z")

</div>

> [@SantiagoOrtiz](#):
>
> ```python
> + dot(w, grad(u_A))*v_A*dx + \ 
> D*dot(grad(u_A), grad(v_A))*dx -
> 
> ```

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`

---

<div class="post-metadata">

**Author:** ![SantiagoOrtiz](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/santiagoortiz/32/201_2.png) [@SantiagoOrtiz](https://fenicsproject.discourse.group/u/SantiagoOrtiz)\
**Post date:** [April 23, 2020, 5:21am UTC](https://fenicsproject.discourse.group/t/gmsh-meshio-fails-with-mixed-elements/2989/5 "2020-04-23T05:21:16Z")

</div>

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.

> **[SantiagoOrtiz/ChemEng\_CaseStudy](https://github.com/SantiagoOrtiz/ChemEng_CaseStudy)**
>
> Selective oxidation of o-xylene to phthalic anhydride: A Convection-diffusion-reaction simulation in single packed-bed tube. - SantiagoOrtiz/ChemEng\_CaseStudy

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

---

<div class="post-metadata">

**Author:** ![bhaveshshrimali](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bhaveshshrimali/32/124_2.png) [@bhaveshshrimali](https://fenicsproject.discourse.group/u/bhaveshshrimali)\
**Post date:** [April 23, 2020, 6:16am UTC](https://fenicsproject.discourse.group/t/gmsh-meshio-fails-with-mixed-elements/2989/6 "2020-04-23T06:16:22Z")

</div>

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

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

```

---

<div class="post-metadata">

**Author:** ![SantiagoOrtiz](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/santiagoortiz/32/201_2.png) [@SantiagoOrtiz](https://fenicsproject.discourse.group/u/SantiagoOrtiz)\
**Post date:** [April 23, 2020, 7:18pm UTC](https://fenicsproject.discourse.group/t/gmsh-meshio-fails-with-mixed-elements/2989/7 "2020-04-23T19:18:25Z")

</div>

@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:

> **[SantiagoOrtiz/ChemEng\_CaseStudy](https://github.com/SantiagoOrtiz/ChemEng_CaseStudy)**
>
> Selective oxidation of o-xylene to phthalic anhydride: A Convection-diffusion-reaction simulation in single packed-bed tube. - SantiagoOrtiz/ChemEng\_CaseStudy