Table of Contents
Hello!
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!
Context
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 (gdm_constitutive.py::hooke & gradient_damagex.py::evaluate_constitutive_law). 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
Warning
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-anyio:
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-fancy-pypi-readme:
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-pytest-asyncio:
python-pytest-rerunfailures:
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-virtualenv:
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
petsc: 3.20.5.776.g55ef5d72e27-1
python-cppimport: 22.08.02.r6.g0849d17-1
python-imageio: 2.33.0-2
python-imageio-ffmpeg:
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
Code
# Copyright (C) 2024 eDgar (fenicsproject.discourse.group)
# 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
# https://fenicsproject.discourse.group/t/10423/2
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,
shape=(msh.geometry.dim,))
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(
basix.CellType, msh.topology.cell_type.name)
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.