# What is wrong with my mesh as related to quadrature-space projections?

Dear everyone,

I have been using different versions of a mesh, and it was working well. Now, I am having issues when projecting. All I get are `np.inf`. If anybody has the time to look at this, I would appreciate it. Thanks in advance, and sorry if it is too simple.

## MWE

``````import ufl
import basix
import dolfinx
from dolfinx import fem
from dolfinx import mesh
import dolfinx.fem.petsc
from mpi4py import MPI
from dolfinx import default_scalar_type as num_type

# dolfinx_cyclic_vM_project.py
def project(v, target_func, bc=[]) -> fem.petsc.LinearProblem:
"""Project UFL expression.

---------

This function was modified from dolfiny (MIT license)
Michal Habera, michal.habera@uni.lu,
Andreas Zilian, andreas.zilian@uni.lu.

--------
https://hplgit.github.io/INF5620/doc/pub/H14/fem/html/._main_fem003.html#___sec9
https://fenicsproject.discourse.group/t/how-to-compute-the-strain-stress-matrices/2756/7
"""
petsc = fem.petsc
# Ensure we have a mesh and attach to measure
V = target_func.function_space
dx = ufl.dx(V.mesh)

# Define variational problem for projection
w = ufl.TestFunction(V)
Pv = ufl.TrialFunction(V)
a = ufl.inner(Pv, w) * dx
L = ufl.inner(v, w) * dx

problem = petsc.LinearProblem(
a, L, bcs=bc, u=target_func,
petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
problem.solve()
return problem

"my_mesh.msh", comm=MPI.COMM_WORLD)[0]
# # This works ##############################
# domain = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
# #########################################

# 589
msh_cell_type = domain.basix_cell()

# 599
msh_cell_type, value_shape=(), degree=2)

# 605
p.x.array[:] = range(len(p.x.array))

# 746
Ve_dg_scal = basix.ufl.element(
"Discontinuous Lagrange", msh_cell_type, degree=2,
shape=())
V_dg_scal = fem.functionspace(domain, Ve_dg_scal)
p_dg_scal = fem.Function(V_dg_scal, dtype=num_type)
p_dg_scal_project = project(p, p_dg_scal)
p_dg_scal.x.array
# array([inf, inf, inf, ..., inf, inf, inf])
``````

## My system

``````- FEniCSx software
dolfinx: 0.8.0_r27835.5a20e2b-1
basix: 0.8.0_r1067.0ed41b0-1
ufl: 2024.2.0.dev0_r3610.72a1bfa-1
ffcx: 0.8.0_r7150.38260c0-1

- Dependencies
blas-openblas: 0.3.26-3
boost: 1.83.0-7
chrpath: 0.17-1
cmake: 3.29.2-1
cython: 3.0.10-3
eigen: 3.4.0-2
glew: 2.2.0-6
gcc-fortran: 13.2.1-6
hdf5-openmpi: 1.14.3-3
ninja: 1.12.0-2
openmpi: 5.0.3-1
parallel: 20240322-1
pybind11: 2.12.0-4
python-anyio: 4.3.0-3
python-build: 1.2.1-3
python-cattrs: 23.2.3-3
python-cffi: 1.16.0-2
python-docutils:
python-fsspec:
python-hatch-vcs: 0.4.0-3
python-hatchling: 1.24.2-1
python-installer: 0.7.0-8
python-mako: 1.3.3-2
python-mpi4py: 3.1.5-5
python-numpy: 1.26.4-2
python-pathspec: 0.12.1-2
python-pygments: 2.17.2-3
python-pytest: 1:8.1.2-1
python-pytest-asyncio: 0.23.6-2
python-pytest-rerunfailures: 14.0-2
python-pyvista: 0.43.4.r0.g3d6e327a7-1
python-rich: 13.7.1-2
python-scikit-build: 0.17.6-3
python-scipy: 1.13.0-2
python-setuptools-scm: 8.0.4-3
python-sphinx:
python-sympy: 1.12.1rc1-1
python-virtualenv: 20.25.0-3
python-wheel: 0.43.0-4
pugixml: 1.14-1
sowing: 1.1.26.9-1
superlu: 6.0.1-1
nanobind: 1.9.2-1
petsc: 3.21.0-1
python-cppimport: 22.08.02-2
python-imageio: 2.34.1-2
python-imageio-ffmpeg: 0.4.9-1
python-meshio: 5.3.5-2
python-scikit-build-core: 0.8.2-1
python-setuptools: 1:69.0.3-6
python-tifffile: 2024.5.3-1
python-wheel: 0.43.0-4
superlu_dist: 8.1.2-1

- Insight-toolkit & dependencies
castxml: 0.6.4-1
eigen: 3.4.0-2
expat: 2.6.2-1
gcc-libs: 13.2.1-6
gdcm: 3.0.12-1
glibc: 2.39-4
hdf5: 1.14.3-3
insight-toolkit: 5.4.0.rc04-1
libjpeg-turbo: 3.0.2-2
libpng: 1.6.43-1
libtiff: 4.6.0-4
meson-python: 0.16.0-3
openjpeg2: 2.5.2-1
openmp: 17.0.6-2
python-flit-core: 3.9.0-4
python-itk-elastix: 5.4.0.rc04-1
python-networkx: 3.3-1
python-pythran: 0.16.0-1
swig: 4.2.1-3
zlib: 1:1.3.1-2

- Others
python-scooby: 0.10.0-1
python-scikit-image: 0.23.2-1
paraview: 5.12.1-8
dcmtk: 3.6.8-1
python-tifffile: 2024.5.3-1
mumps: 5.7.1-1
hypre: 2.31.0-1
metis: 5.2.1-2
scalapack: 2.2.0-1
superlu_dist: 8.1.2-1
gklib: 5.1.1-4
mfront: 20240324.r2770.1e94d01fd-1
libfyaml: 0.9-1

- Operating system
Arch Linux: 6.9.0-rc6-gnu-rc6
CPU: 12th Gen Intel(R) Core(TM) i7-1260P
``````

I could not upload the mesh file (the link will expire on the 18th of July 2024):

Looks like the mesh might not correctly connected. Have you tried running it through @dokkenâ€™s checker script Mesh topology verification, testing if there are cells that are only connected through edges or vertices Â· GitHub ?

1 Like

Hi, @nate! Thank you very much. I donâ€™t really know how to use the script. It seems to be made for triangles and tetrahedra, and expects some data structure that I donâ€™t understand. I was trying to get the relevant bits (geometry and topology) out of the `.msh`, but I donâ€™t know how. For example, `topology` in `mesh_checker.py` seems to be a Numpy array (`topology.shape`), but I got a dictionary:

``````import dolfinx
import gmsh

gmsh.initialize()
msh = gmsh.open("my_mesh.msh")

# https://fenicsproject.discourse.group/t/8741
# help(dolfinx.io.gmshio)
geometry = dolfinx.io.gmshio.extract_geometry(gmsh.model)
topology = dolfinx.io.gmshio.extract_topology_and_markers(gmsh.model)

print(topology)
``````

[Contents of topology (in case that it helps)]

``````{1: {'topology': array([[ 15,  87],
[ 87,  88],
[ 88,  89],
[ 89,  14],
[ 15, 104],
[104, 105],
[105, 106],
[106,  20],
[ 21, 113],
[113, 114],
[114, 115],
[115,  20],
[ 21, 129],
[129, 130],
[130, 131],
[131,  24],
[  7,  52],
[ 52,  53],
[ 53,  54],
[ 54,   8],
[  7, 135],
[135, 136],
[136, 137],
[137,  24]], dtype=uint64), 'cell_data': array([1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3,
3, 3], dtype=uint64)}, 3: {'topology': array([[  4,  46, 156,  40],
[ 40, 156, 157,  41],
[ 41, 157, 158,  42],
...,
[270, 271,  26,  25],
[271, 272,  27,  26],
[272, 123,   1,  27]], dtype=uint64), 'cell_data': array([4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5], dtype=uint64)}}
``````

I donâ€™t really think that the following would make a difference, but Gmsh reports a 2D geometry, while DOLFINx says that itâ€™s 3D.

``````import dolfinx
from mpi4py import MPI
import gmsh

gmsh.initialize()
msh = gmsh.open("my_mesh.msh")
# https://gmsh.info/doc/texinfo/gmsh.html#x1
print('Gmsh model ' + gmsh.model.getCurrent() + ' (' +
str(gmsh.model.getDimension()) + 'D)')

"my_mesh.msh", comm=MPI.COMM_WORLD)[0]
print("DOLFINx domain " + domain.name + f" ({domain.geometry.dim}D)")
``````

[Results of previous block: Gmsh shows 2D; DOLFINx shows 3D]

``````Gmsh model my_mesh (2D)
Dolfinx domain mesh (3D)
``````

What may be relevant is that other types of data can be interpolated

``````import ufl
import basix
import dolfinx
from dolfinx import fem
from mpi4py import MPI
from dolfinx import default_scalar_type as num_type

"my_mesh.msh", comm=MPI.COMM_WORLD)[0]

# 589
msh_cell_type = domain.basix_cell()

# 590
Ve_cg_vec = basix.ufl.element(
'Lagrange', msh_cell_type, degree=2, shape=(domain.geometry.dim,))
V_cg_vec = fem.functionspace(domain, Ve_cg_vec)

# 606
u = fem.Function(V_cg_vec, name="Total_displacement")

# 793
V_dg_vec = fem.functionspace(
domain, ("Discontinuous Lagrange", 2, (domain.geometry.dim,)))
u_dg_vec = fem.Function(
V_dg_vec, dtype=num_type, name="displacement")

u.x.array[:] = range(len(u.x.array))
u_dg_vec.interpolate(u)
print(u_dg_vec.x.array)
``````
``````[-1.95676808e-15  1.00000000e+00  2.00000000e+00 ...  3.16800000e+03
3.16900000e+03  3.17000000e+03]
``````

However, if I try to use `project(â€¦)` with that data, I get the same `np.inf`.

``````import ufl
import basix
import dolfinx
from dolfinx import fem
from dolfinx import mesh
import dolfinx.fem.petsc
from mpi4py import MPI
from dolfinx import default_scalar_type as num_type

# dolfinx_cyclic_vM_project.py
def project(v, target_func, bc=[]) -> fem.petsc.LinearProblem:
"""Project UFL expression.

---------

This function was modified from dolfiny (MIT license)
Michal Habera, michal.habera@uni.lu,
Andreas Zilian, andreas.zilian@uni.lu.

--------
https://hplgit.github.io/INF5620/doc/pub/H14/fem/html/._main_fem003.html#___sec9
https://fenicsproject.discourse.group/t/how-to-compute-the-strain-stress-matrices/2756/7
"""
petsc = fem.petsc
# Ensure we have a mesh and attach to measure
V = target_func.function_space
dx = ufl.dx(V.mesh)

# Define variational problem for projection
w = ufl.TestFunction(V)
Pv = ufl.TrialFunction(V)
a = ufl.inner(Pv, w) * dx
L = ufl.inner(v, w) * dx

problem = petsc.LinearProblem(
a, L, bcs=bc, u=target_func,
petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
problem.solve()
return problem

"my_mesh.msh", comm=MPI.COMM_WORLD)[0]

# 589
msh_cell_type = domain.basix_cell()

# 590
Ve_cg_vec = basix.ufl.element(
'Lagrange', msh_cell_type, degree=2, shape=(domain.geometry.dim,))
V_cg_vec = fem.functionspace(domain, Ve_cg_vec)

# 606
u = fem.Function(V_cg_vec, name="Total_displacement")

# 793
V_dg_vec = fem.functionspace(
domain, ("Discontinuous Lagrange", 2, (domain.geometry.dim,)))
u_dg_vec = fem.Function(
V_dg_vec, dtype=num_type, name="displacement")

u.x.array[:] = range(len(u.x.array))
project(u, u_dg_vec)
print(u_dg_vec.x.array)
``````

[All `inf` when trying to project a variable which can be interpolated]

``````[inf inf inf ... inf inf inf]
``````

Thanks again for your attention and time.

There is clearly something fishy with your mesh, as trying to read it with meshio yields:

``````import meshio
``````
``````Traceback (most recent call last):
File "/root/shared/mwe.py", line 5, in <module>
File "/usr/lib/python3/dist-packages/meshio/_helpers.py", line 71, in read
File "/usr/lib/python3/dist-packages/meshio/_helpers.py", line 103, in _read_file
File "/usr/lib/python3/dist-packages/meshio/gmsh/main.py", line 19, in read
File "/usr/lib/python3/dist-packages/meshio/gmsh/main.py", line 48, in read_buffer
File "/usr/lib/python3/dist-packages/meshio/gmsh/_gmsh41.py", line 104, in read_buffer
return Mesh(
File "/usr/lib/python3/dist-packages/meshio/_mesh.py", line 174, in __init__
raise ValueError(
ValueError: Incompatible cell data 'gmsh:physical'. 89 cell blocks, but 'gmsh:physical' has 26 blocks.
``````

This seem to indicate that there are elements in your mesh that are either duplicated or marked not marked with a physical marker.

Could you elaborate on how you made the mesh?

You need to send gdim=2 into read_from_msh to make dolfinx realize it is a 2D mesh.

Finally, if you change the lu backend to `mumps` the code produces non-inf values:

``````        petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"})
``````

something to always check is if `ksp.getConvergedReason()` is bigger than 0, i.e.

``````    if problem.solver.getConvergedReason() <= 0:
problem.solver.view()
print(problem.solver.getConvergedReason())
assert problem.solver.getConvergedReason()>0, "Linear solver did not converge"
``````

which will yield -11 without setting mumps as the backend, meaning that the lu solver failed:

``````  KSP_DIVERGED_PC_FAILED                 = -11,
KSP_DIVERGED_PCSETUP_FAILED_DEPRECATED = -11,
``````
1 Like

Dear JĂ¸rgen @dokken, thank you very much! May your light shine for us for many years to come.

I want to confirm that I get the same thing.

What I found is that when exporting, I should not select â€śSave all elementsâ€ť. If I do that, `meshio` can read the mesh (NINJAFILES - my_new_mesh.msh):

``````  import meshio
print(msh)
``````

#+results:

``````  <meshio mesh object>
Number of points: 273
Number of cells:
line: 4
line: 4
line: 4
line: 4
line: 4
line: 4
Cell sets: left, bottom, right, Polymer, Fibre, gmsh:bounding_entities
Point data: gmsh:dim_tags
Cell data: gmsh:physical, gmsh:geometrical
Field data: left, bottom, right, Polymer, Fibre
``````

Mostly irrelevant now, but I am attaching some pictures from the FLTK interface of Gmsh which shows the structure (for completeness). Most of the geometrical elements are copies of rotations from surfaces.

Also, thanks (really THANKS) to both of you, this is the new version for anyone who finds it useful.

``````  # Projection of
# 1. a function with quadrature elements onto a function
#    with discontinuous Lagrange elements
# 2. a function with Lagrange elements onto a function with
#    discontinuous Lagrange elements
#
#   (GNU Public License version 3.0)

import ufl
import basix
import dolfinx
from dolfinx import fem
from dolfinx import mesh
import dolfinx.fem.petsc
from mpi4py import MPI
from dolfinx import default_scalar_type as num_type

# dolfinx_cyclic_vM_project.py
def project(v, target_func, bc=[]) -> fem.petsc.LinearProblem:
"""Project UFL expression.

---------

This function was modified from dolfiny (MIT license)
Michal Habera, michal.habera@uni.lu,
Andreas Zilian, andreas.zilian@uni.lu.

--------
https://hplgit.github.io/INF5620/doc/pub/H14/fem/html/._main_fem003.html#___sec9
https://fenicsproject.discourse.group/t/how-to-compute-the-strain-stress-matrices/2756/7
"""
petsc = fem.petsc
# Ensure we have a mesh and attach to measure
V = target_func.function_space
dx = ufl.dx(V.mesh)

# Define variational problem for projection
w = ufl.TestFunction(V)
Pv = ufl.TrialFunction(V)
a = ufl.inner(Pv, w) * dx
L = ufl.inner(v, w) * dx

problem = petsc.LinearProblem(
a, L, bcs=bc, u=target_func,
petsc_options={"ksp_type": "preonly", "pc_type": "lu",
"pc_factor_mat_solver_type": "mumps"})
problem.solve()

if problem.solver.getConvergedReason() <= 0:
problem.solver.view()
print("Projection convergence criterion: ",
problem.solver.getConvergedReason())
assert problem.solver.getConvergedReason()>0, \
"Linear solver of projection did not converge"

return problem

# "/tmp/my_mesh.msh", comm=MPI.COMM_WORLD, gdim=2)[0]
"/tmp/my_new_mesh.msh", comm=MPI.COMM_WORLD, gdim=2)[0]

# 589
msh_cell_type = domain.basix_cell()

# * Create and project function with Lagrange elements
# 590
Ve_cg_vec = basix.ufl.element(
'Lagrange', msh_cell_type, degree=2, shape=(domain.geometry.dim,))
V_cg_vec = fem.functionspace(domain, Ve_cg_vec)

# 606
u = fem.Function(V_cg_vec)

# 793
V_dg_vec = fem.functionspace(
domain, ("Discontinuous Lagrange", 2, (domain.geometry.dim,)))
u_dg_vec = fem.Function(
V_dg_vec, dtype=num_type)

u.x.array[:] = range(len(u.x.array))
project(u, u_dg_vec)
print(f"{u_dg_vec.x.array = }")

# * Create and project function with quadrature elements
# 599
msh_cell_type, value_shape=(), degree=2)

# 605
p.x.array[:] = range(len(p.x.array))

# 746
Ve_dg_scal = basix.ufl.element(
"Discontinuous Lagrange", msh_cell_type, degree=2,
shape=())
V_dg_scal = fem.functionspace(domain, Ve_dg_scal)
p_dg_scal = fem.Function(V_dg_scal, dtype=num_type)
p_dg_scal_project = project(p, p_dg_scal)
print(f"{p_dg_scal.x.array = }")
``````

#+results:

``````  Info    : Reading '/tmp/my_mesh.msh'...
Info    : 89 entities
Info    : 273 nodes
Info    : 441 elements