# Trivial form fails on xdmf mesh

**URL:** <https://fenicsproject.discourse.group/t/trivial-form-fails-on-xdmf-mesh/8370>\
**Category:** Uncategorized\
**Created:** [May 20, 2022, 10:24am UTC](https://fenicsproject.discourse.group/t/trivial-form-fails-on-xdmf-mesh/8370 "2022-05-20T10:24:18Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![hawkspar](https://avatars.discourse-cdn.com/v4/letter/h/b487fb/32.png) [@hawkspar](https://fenicsproject.discourse.group/u/hawkspar)\
**Post date:** [May 20, 2022, 10:24am UTC](https://fenicsproject.discourse.group/t/trivial-form-fails-on-xdmf-mesh/8370/1 "2022-05-20T10:24:18Z")

</div>

Consider the short MWE :

```python
import ufl
import dolfinx as dfx
from dolfinx.io import XDMFFile
from mpi4py.MPI import COMM_WORLD

with XDMFFile(COMM_WORLD, "mesh.xdmf", "r") as file:
    mesh = file.read_mesh(name="Grid")

FE=ufl.FiniteElement("Lagrange",mesh.ufl_cell(),1)
V = dfx.FunctionSpace(mesh,FE)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

#r = ufl.SpatialCoordinate(mesh)[1]
r = dfx.Function(V)
r.interpolate(lambda x: x[1])

n = ufl.inner(u,v)*r**2*ufl.dx

N = dfx.fem.assemble_matrix(n)
N.assemble()

x, y = N.createVecs()
x.set(1)
N.mult(x,y)
t = dfx.Function(V)
t.vector[:]=y
with XDMFFile(COMM_WORLD, "sanity_check_N_dummy_test.xdmf","w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_function(t)

```

I was expecting a smooth square going from 0 at the bottom to max at the top. Instead I get :

 ![ugly](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/9/949a926c5cf75b1359793b43b52c8320bb176ac7.png)

The [mesh](https://seminaris.polytechnique.fr/share/s/PKKzkyKeBbPR2eM) is here. I strongly suspect it could be the source of the error since I couldn’t reproduce that error with a `mesh=dfx.generation.UnitSquareMesh(COMM_WORLD,100,100)`. But that gives me no hint on how to fix the original mesh…

This is a non-uniform structured mesh of quads originally generated using `gmsh`, then put through an entire ordeal, but it still looks fine to me. Besides, I’m surprised the code doesn’t fail at reading if something is really wrong with the code.

Does anyone have an idea as to the root cause of the problem ?

---

<div class="post-metadata">

**Author:** ![hawkspar](https://avatars.discourse-cdn.com/v4/letter/h/b487fb/32.png) [@hawkspar](https://fenicsproject.discourse.group/u/hawkspar)\
**Post date:** [May 20, 2022, 11:45am UTC](https://fenicsproject.discourse.group/t/trivial-form-fails-on-xdmf-mesh/8370/2 "2022-05-20T11:45:24Z")

</div>

I’m trying with the identity now. Having `n=ufl.inner(u,v)*ufl.dx` gives another ugly graph :

 ![ugly2](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/a/a58334043cf627f4b035ec453346eb03b8696dd7.png)

I don’t understand why this graph isn’t filled with red throughout and why there appears to be some sort of scar of zeros at a random position…

I’m in complex mode, running sequentially my docker container on an OpenSUSE Leap 15.2 machine.

---

<div class="post-metadata">

**Author:** ![hawkspar](https://avatars.discourse-cdn.com/v4/letter/h/b487fb/32.png) [@hawkspar](https://fenicsproject.discourse.group/u/hawkspar)\
**Post date:** [May 20, 2022, 12:07pm UTC](https://fenicsproject.discourse.group/t/trivial-form-fails-on-xdmf-mesh/8370/3 "2022-05-20T12:07:17Z")

</div>

I’m pretty sure this problem is unrelated to [this one](https://fenicsproject.discourse.group/t/assemble-matrix-in-dolfinx/7547) since in my version of dolfinx `dolfinx.fem.form` is still a module.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [May 20, 2022, 12:21pm UTC](https://fenicsproject.discourse.group/t/trivial-form-fails-on-xdmf-mesh/8370/4 "2022-05-20T12:21:00Z")

</div>

The problem is that the interior cells in your mesh are not connected. This can be illustrated by running the following:

```python
import ufl
import dolfinx as dfx
from dolfinx.io import XDMFFile
from mpi4py.MPI import COMM_WORLD
import numpy as np
from petsc4py import PETSc

with XDMFFile(COMM_WORLD, "mesh.xdmf", "r") as file:
    mesh = file.read_mesh(name="Grid")

FE=ufl.FiniteElement("Lagrange",mesh.ufl_cell(),1)
V = dfx.fem.FunctionSpace(mesh,FE)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

n = ufl.inner(u,v)*ufl.dx

N = dfx.fem.petsc.assemble_matrix(dfx.fem.form(n))
N.assemble()

y = dfx.fem.Function(V)

facets = dfx.mesh.locate_entities_boundary(mesh, mesh.topology.dim-1, lambda x: np.full(x.shape[1], True))
dofs = dfx.fem.locate_dofs_topological(V, mesh.topology.dim-1, facets)
c = dfx.fem.Constant(mesh, PETSc.ScalarType(3))
bc = dfx.fem.dirichletbc(c, dofs, V)
y.x.array[bc.dof_indices()[0]] = 2

with XDMFFile(COMM_WORLD, "sanity_check_N_dummy_test.xdmf","w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_function(y)

```

This marks all dofs that are connected to a boundary facet (a facet connected to only one cell) with the value 2.  
Looking at this in Paraview, you will observe that there are two rays through your domain, corresponding to the rays through your domain.

---

<div class="post-metadata">

**Author:** ![hawkspar](https://avatars.discourse-cdn.com/v4/letter/h/b487fb/32.png) [@hawkspar](https://fenicsproject.discourse.group/u/hawkspar)\
**Post date:** [May 20, 2022, 12:39pm UTC](https://fenicsproject.discourse.group/t/trivial-form-fails-on-xdmf-mesh/8370/5 "2022-05-20T12:39:29Z")

</div>

Many thanks for your quick response ; I can’t seem to get your code running on my machine though. I had to rollback a few of the synthax, I’m standing on :

```python
import ufl
import dolfinx as dfx
from dolfinx.io import XDMFFile
from mpi4py.MPI import COMM_WORLD
import numpy as np

with XDMFFile(COMM_WORLD, "../cases/nozzle/nozzle.xdmf", "r") as file:
    mesh = file.read_mesh(name="Grid")

FE=ufl.FiniteElement("Lagrange",mesh.ufl_cell(),1)
V = dfx.FunctionSpace(mesh,FE)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

n = ufl.inner(u,v)*ufl.dx

N = dfx.fem.assemble_matrix(n)
N.assemble()

y = dfx.Function(V)

facets = dfx.mesh.locate_entities_boundary(mesh, mesh.topology.dim-1, lambda x: np.full(x.shape[1], True))
dofs = dfx.fem.locate_dofs_topological(V, mesh.topology.dim-1, facets)
c=dfx.Function(V)
with c.vector.localForm() as loc: loc.set(3)
bc = dfx.DirichletBC(c, dofs, V)
y.x.array[bc.dof_indices()[0]] = 2

with XDMFFile(COMM_WORLD, "sanity_check_N_dummy_test.xdmf","w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_function(y)

```

Which fails at the DirichletBC with traceback :

```auto
Traceback (most recent call last):
  File "/home/shared/src/MWE2.py", line 26, in <module>
    bc = dfx.DirichletBC(c, dofs, V)
  File "/usr/local/dolfinx-complex/lib/python3.8/dist-packages/dolfinx/fem/dirichletbc.py", line 177, in __init__
    self._cpp_object = dirichletbc_obj(dtype)(_value, dofs, _V)
TypeError: __init__ (): incompatible constructor arguments. The following argument types are supported:
    1. dolfinx.cpp.fem.DirichletBC_complex128(arg0: dolfinx::fem::Function<std::complex<double> >, arg1: numpy.ndarray[numpy.int32])
    2. dolfinx.cpp.fem.DirichletBC_complex128(arg0: dolfinx::fem::Function<std::complex<double> >, arg1: List[numpy.ndarray[numpy.int32][2]], arg2: dolfinx::fem::FunctionSpace)

Invoked with: <dolfinx.cpp.fem.Function_complex128 object at 0x7f4bc023a370>, array([0, 1, 3, ..., 76138, 76139, 76140], dtype=int32), <dolfinx.cpp.fem.FunctionSpace object at 0x7f4bc0245030>

Did you forget to `#include <pybind11/stl.h>`? Or <pybind11/complex.h>,
<pybind11/functional.h>, <pybind11/chrono.h>, etc. Some automatic
conversions are optional and require extra headers to be included
when compiling your pybind11 module.

```

Do you think this mesh break up can be fixed inside dolfinx or would it be easier to start anew ?

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [May 20, 2022, 1:05pm UTC](https://fenicsproject.discourse.group/t/trivial-form-fails-on-xdmf-mesh/8370/6 "2022-05-20T13:05:05Z")

</div>

Just remove `V` from

> [@hawkspar](#):
>
> `bc = dfx.DirichletBC(c, dofs, V)`

i.e.  
`bc = dfx.DirichletBC(c, dofs)`.

> [@hawkspar](#):
>
> Do you think this mesh break up can be fixed inside dolfinx or would it be easier to start anew ?

The mesh break up should be fixed in the mesh generator, stitching together meshes in DOLFINx is non-trivial (especially in parallel). I’ve done some work on thjis with my [dolfinx\_mpc](https://github.com/jorgensd/dolfinx_mpc), but it is not the recommended way of resolving such issues.

---

<div class="post-metadata">

**Author:** ![hawkspar](https://avatars.discourse-cdn.com/v4/letter/h/b487fb/32.png) [@hawkspar](https://fenicsproject.discourse.group/u/hawkspar)\
**Post date:** [May 20, 2022, 1:28pm UTC](https://fenicsproject.discourse.group/t/trivial-form-fails-on-xdmf-mesh/8370/7 "2022-05-20T13:28:22Z")

</div>

Ok thanks I’ll try to read a fresh mesh again. Do you expect it to solve both the “scar” and the values problem ? Looking at the first figure, I can see how the problem you’re mentioning breaks specific regions of the graph, but not why my `r^2` is not looking like a square…

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [May 20, 2022, 1:51pm UTC](https://fenicsproject.discourse.group/t/trivial-form-fails-on-xdmf-mesh/8370/8 "2022-05-20T13:51:36Z")

</div>

I would suggest starting by checking that interpolating or projecting r^2 into a suitable function space gives the expected function.

In general, if something works on the unit square or a rectangle mesh, it should work on any mesh.