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

dolfinx
, ,
1

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.

    Copyright
    ---------
    SPDX-License-Identifier: GPL-3.0-only
    (GNU Public License version 3.0)

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

    See also
    --------
    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

domain = dolfinx.io.gmshio.read_from_msh(
    "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
Ve_quad_scal = basix.ufl.quadrature_element(
    msh_cell_type, value_shape=(), degree=2)
V_quad_scal = fem.functionspace(domain, Ve_quad_scal)

# 605
p = fem.Function(V_quad_scal)
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
  adios2: 2.10.0-2
  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-fancy-pypi-readme: 24.1.0-3
  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-pyproject-metadata: 0.8.0-1
  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-lazy-loader: 0.4-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):

2

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
3

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)')

domain = dolfinx.io.gmshio.read_from_msh(
    "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


domain = dolfinx.io.gmshio.read_from_msh(
    "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.

    Copyright
    ---------
    SPDX-License-Identifier: GPL-3.0-only
    (GNU Public License version 3.0)

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

    See also
    --------
    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

domain = dolfinx.io.gmshio.read_from_msh(
    "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.

4

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

import meshio
mesh = meshio.read("my_mesh.msh")

Traceback (most recent call last):
  File "/root/shared/mwe.py", line 5, in <module>
    mesh = meshio.read("my_mesh.msh")
  File "/usr/lib/python3/dist-packages/meshio/_helpers.py", line 71, in read
    return _read_file(Path(filename), file_format)
  File "/usr/lib/python3/dist-packages/meshio/_helpers.py", line 103, in _read_file
    return reader_map[file_format](str(path))
  File "/usr/lib/python3/dist-packages/meshio/gmsh/main.py", line 19, in read
    mesh = read_buffer(f)
  File "/usr/lib/python3/dist-packages/meshio/gmsh/main.py", line 48, in read_buffer
    return reader.read_buffer(f, is_ascii, data_size)
  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,

from KSPConvergedReason — PETSc 3.21.2 documentation

1 Like
5

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
  msh = meshio.read("/tmp/my_new_mesh.msh")
  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
      quad: 16
      quad: 8
      quad: 8
      quad: 16
      quad: 16
      quad: 16
      quad: 8
      quad: 8
      quad: 16
      quad: 16
      quad: 8
      quad: 8
      quad: 16
      quad: 16
      quad: 8
      quad: 8
      quad: 16
      quad: 16
      quad: 16
      quad: 16
    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.

mesh_structure
mesh_structure6517Ă—2739 980 KB

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
  #
  # Copyright 2024 edgar (fenicsproject.discourse.group)
  # License:
  #   SPDX-License-Identifier: GPL-3.0-only
  #   (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.

      Copyright
      ---------
      SPDX-License-Identifier: GPL-3.0-only
      (GNU Public License version 3.0)

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

      See also
      --------
      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
      https://fenicsproject.discourse.group/t/what-is-wrong-with-my-mesh-as-related-to-quadrature-space-projections/15028
      """
      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

  # * Read mesh file or load Gmsh model
  domain = dolfinx.io.gmshio.read_from_msh(
      # "/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
  Ve_quad_scal = basix.ufl.quadrature_element(
      msh_cell_type, value_shape=(), degree=2)
  V_quad_scal = fem.functionspace(domain, Ve_quad_scal)

  # 605
  p = fem.Function(V_quad_scal)
  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
  Info    : Done reading '/tmp/my_mesh.msh'
  Projection convergence criterion:  4
  u_dg_vec.x.array = array([-1.16662321e-14,  1.00000000e+00,  2.00000000e+00, ...,
          2.08900000e+03,  2.11200000e+03,  2.11300000e+03])
  Projection convergence criterion:  4
  p_dg_scal.x.array = array([-1.28287321e+00,  2.19328011e+00,  4.87108255e-01, ...,
          1.02329427e+03,  1.02244613e+03,  1.02135152e+03])
6

Yes, do not select that as duplicate elements (and unused nodes will be stored).