How to assign values to the array of DOF based on mesh markers

Table of Contents

  1. Context
  2. Warning
  3. My system
  4. Code


I was going to ask a question to know how to assign values to the array of DOF based on mesh markers. However, writing the question helped me to solve it by myself. Nonetheless, I want to share how I did it, in case that someone finds it useful. Thank you!


This helps if you have subdomains, and–for whatever reason–you want to access the values of the polynomial coefficients of each quadrature point (DOF) as an array.

When the markers or tags are already known (e.g. as physical groups from Gmsh), DOLFINx usually helps to do away with the complexity of dealing with the arrays. Some times, however, you may want to access each element of the array which holds all the DOF. When there are no subdomains, this can be done by “unrolling” (flattening and reshaping) the whole array ( & However, with subdomains, the shape of the array for each subdomain will depend on the polynomial degree and the number of quadrature (or integration) points. Even when the indices for each element is previously known (from markers), the indices may not be sequential, and the DOF will be indexed depending on the type of element.

Since DOLFINx already provides a way to go from elements’ index to index of DOF, the only thing which is presented here is how to do it when there are subdomains.

It may be evident to others, but it wasn’t for me. I hope that it serves someone


The calculations by themselves do not represent anything meaningful. DO NOT USE to estimate any physical variable.

My system

- FEniCSx software
  dolfinx: 0.8.0.dev0_r27744.59bcec3-1
  basix: 0.8.0.dev0_r1019.c9ca4e5-1
  ufl: 2023.3.0.dev0_r3597.0f2cb21-1
  ffcx: 0.8.0.dev0_r7125.88871c5-1

- Dependencies
  adios2: 2.9.2-2
  blas-openblas: 0.3.26-3
  boost: 1.83.0-5
  chrpath: 0.17-1
  cmake: 3.29.0-1
  cython: 3.0.9-1
  eigen: 3.4.0-2
  gcc-fortran: 13.2.1-5
  hdf5-openmpi: 1.14.3-2
  ninja: 1.11.1-3
  openmpi: 5.0.2-7
  parallel: 20240122-1
  pybind11: 2.11.1-2
  python-build: 1.1.1-2
  python-cattrs: 23.2.3-1
  python-cffi: 1.16.0-1
  python-docutils: 1:0.20.1-2
  python-fsspec: 2023.9.2-1
  python-hatch-vcs: 0.4.0-1
  python-hatchling: 1.21.1-2
  python-installer: 0.7.0-4
  python-mako: 1.3.2-1
  python-mpi4py: 3.1.5-2
  python-numpy: 1.26.4-1
  python-pathspec: 0.12.1-1
  python-pygments: 2.17.2-1
  python-pyproject-metadata: 0.7.1-2
  python-pytest: 1:8.1.1-1
  python-rich: 13.7.1-1
  python-scikit-build: 0.17.6-2
  python-scipy: 1.12.0-1
  python-setuptools-scm: 8.0.4-1
  python-sphinx: 7.2.6-1
  python-sympy: 1.12-2
  python-wheel: 0.43.0-1
  pugixml: 1.14-1
  sowing: 1.1.26.p8-1
  superlu: 6.0.1-1
  nanobind: 1.8.0-1
  python-cppimport: 22.08.02.r6.g0849d17-1
  python-imageio: 2.33.0-2
  python-meshio: 5.3.4-2
  python-scikit-build-core: 0.5.0-1
  python-setuptools: 1:69.0.3-4
  python-tifffile: 2024.2.12-1
  python-wheel: 0.43.0-1
  superlu_dist: 8.1.2-1

- Operating system
  Arch Linux: 6.8.1-arch1-1
  CPU: 11th Gen Intel(R) Core(TM) i7-1185G7 @ 3.00GHz


    # Copyright (C) 2024 eDgar (
    #   SPDX-License-Identifier: GPL-3.0-only (GNU Public License version 3.0)
    import dolfinx as dfx
    import basix
    import ufl
    from dolfinx import default_scalar_type as num_type
    from dolfinx import mesh
    from dolfinx import fem
    from mpi4py import MPI
    import numpy as np
    # * Create a square mesh
    msh = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
    # * Retrieve the cell tags
    #   Assume that this is previously known (as if physical
    #   groups from Gmsh)
    bottom_cells = mesh.locate_entities(
        msh, msh.topology.dim, lambda x: x[1] < 0.5)
    top_cells = mesh.locate_entities(
        msh, msh.topology.dim, lambda x: x[1] >= - np.infty)
    # * Get the meshtags
    # Based on
    bottom_markers = np.arange(len(bottom_cells), dtype=np.int32)
    top_markers = np.arange(len(top_cells), dtype=np.int32)
    bottom_ct = dfx.mesh.meshtags(
        msh, msh.topology.dim, bottom_cells, bottom_markers)
    top_ct = dfx.mesh.meshtags(
        msh, msh.topology.dim, top_cells, top_markers)
    # * Set the polynomial degree
    deg = 2
    q_deg = 4
    msh_cell_type = msh.basix_cell()
    # * Create the variable space
    # ** Elements for displacement and equivalent strain
    Ed = basix.ufl.element(
        'Lagrange', msh_cell_type, degree=deg,
    Ee = basix.ufl.element(
        "Lagrange", msh.basix_cell(), degree=deg)
    # ** Space with mixed elements
    mixed_el = basix.ufl.mixed_element([Ed, Ee])
    V = dfx.fem.functionspace(msh, mixed_el)
    u = dfx.fem.Function(V)
    d, e = ufl.split(u)
    # ** Quadrature points
    basix_celltype = getattr(
    q_points, wts = basix.make_quadrature(
                basix_celltype, q_deg)
    # # * This is not the right way to access the array by markers
    # len(u.x.array)
    # len([u.x.array[bottom_ct.find(val)] for val in bottom_markers])
    # len([u.x.array[top_ct.find(val)] for val in top_markers])
    # print(f"{len(u.x.array) = }")
    e = ufl.grad(d)
    strain_expr = dfx.fem.Expression(
        ufl.as_vector([e[0, 0], e[1, 1], e[0, 1]]), q_points)
    strain_bottom = strain_expr.eval(msh, bottom_cells)
    m = 1
    l = 1
    strain_bottom_flat = strain_bottom.flatten()
    C = np.array([[m, l, 0], [l, m, 0], [0, 0, m]])
    res = strain_bottom_flat.reshape(-1, 3) @ C

Final note

Thank you for FEniCSx, the time, the answers, the effort and this forum. If somebody has suggestions on how to improve the code, please, let me know.