Interpolating UFL expressions into quadrature functions

dolfinx
1

Hi all,

I am looking into interpolating ufl expressions into quadrature functions. Consider for instance the example borrowed from

However, the following MWE

from dolfinx.mesh import create_unit_cube
from dolfinx import fem
from mpi4py import MPI
from ufl import TensorElement, SpatialCoordinate, grad, dx
import basix, numpy as np
from petsc4py import PETSc

metadata = {"quadrature_degree":5}
mesh = create_unit_cube(MPI.COMM_WORLD, 5, 5, 5)
Wq = TensorElement("Quadrature", mesh.ufl_cell(), degree=metadata["quadrature_degree"], symmetry=True, quad_scheme="default")
Vq = fem.FunctionSpace(mesh, Wq)

quadrature_degree = metadata["quadrature_degree"]
quadrature_points, wts = basix.make_quadrature(basix.CellType.triangle, quadrature_degree)


expr = grad(SpatialCoordinate(mesh))
e_expr = fem.Expression(expr, quadrature_points)
map_c = mesh.topology.index_map(mesh.topology.dim)
num_cells = map_c.size_local + map_c.num_ghosts
cells = np.arange(0, num_cells, dtype=np.int32)

e_eval = e_expr.eval(cells)
# print(e_eval)
e_Q = fem.Function(Vq)
with e_Q.vector.localForm() as e_Q_local:
    e_Q_local.setBlockSize(e_Q.function_space.dofmap.bs)
    e_Q_local.setValuesBlocked(Vq.dofmap.list.array, e_eval, addv=PETSc.InsertMode.INSERT)

throws an error when trying to set the corresponding values

    e_Q_local.setValuesBlocked(Vq.dofmap.list.array, e_eval, addv=PETSc.InsertMode.INSERT)
  File "PETSc/Vec.pyx", line 1136, in petsc4py.PETSc.Vec.setValuesBlocked
  File "PETSc/petscvec.pxi", line 360, in petsc4py.PETSc.vecsetvalues
ValueError: incompatible array sizes: ni=11250, nv=47250, bs=6

Would anyone have a clue as to what Iā€™m doing incorrectly?

1 Like
2

Hi @bhaveshshrimali , you must create quadrature rule consistently with the Quadrature function space (the same cell type),

quadrature_points, wts = basix.make_quadrature(basix.CellType.tetrahedron, quadrature_degree)

or use

basix.cell.string_to_type(mesh.ufl_cell().cellname())

as cell type for a general approach.

3 Likes
3

Indeed, Thanks @michalhabera !