Trivial form fails on xdmf mesh

1

Consider the short MWE :

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
ugly1356×311 29 KB

The mesh 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 ?

2

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

ugly2
ugly21345×308 35.4 KB

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.

3

I’m pretty sure this problem is unrelated to this one since in my version of dolfinx dolfinx.fem.form is still a module.

4

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

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.

5

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 :

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 :

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 ?

6

Just remove V from

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

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, but it is not the recommended way of resolving such issues.

7

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…

8

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.