Mapped quadrature points, weights and solution at these quadrature points

General
2

Consider the following:

import dolfinx
from mpi4py import MPI
import ufl
import basix

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 2, 2)
quadrature_degree = 1
V = dolfinx.fem.FunctionSpace(mesh, ufl.FiniteElement(
    "Quadrature", mesh.ufl_cell(), quadrature_degree, quad_scheme="default"))


quadrature_points, wts = basix.make_quadrature(
    basix.cell.string_to_type(mesh.topology.cell_name()), quadrature_degree)
x = ufl.SpatialCoordinate(mesh)
x_expr = dolfinx.fem.Expression(x, quadrature_points)

detJ = dolfinx.fem.Expression(ufl.JacobianDeterminant(mesh), quadrature_points)

for i in range(mesh.topology.index_map(mesh.topology.dim).size_local):
    print(
        f"Cell {i}, quadrature points {x_expr.eval([i])}, detJ {detJ.eval([i])}")

Returning

Cell 0, quadrature points [[0.83333333 0.16666667]], detJ [[0.25]]
Cell 1, quadrature points [[0.66666667 0.33333333]], detJ [[-0.25]]
Cell 2, quadrature points [[0.33333333 0.16666667]], detJ [[0.25]]
Cell 3, quadrature points [[0.83333333 0.66666667]], detJ [[0.25]]
Cell 4, quadrature points [[0.16666667 0.33333333]], detJ [[-0.25]]
Cell 5, quadrature points [[0.66666667 0.83333333]], detJ [[-0.25]]
Cell 6, quadrature points [[0.33333333 0.66666667]], detJ [[0.25]]
Cell 7, quadrature points [[0.16666667 0.83333333]], detJ [[-0.25]]

for quadrature degree 1
and

Cell 0, quadrature points [[0.9454805  0.11596668 0.88403332 0.0545195  0.9454805  0.32951381
  0.67048619 0.0545195  0.88403332 0.32951381 0.67048619 0.11596668]], detJ [[0.25 0.25 0.25 0.25 0.25 0.25]]
Cell 1, quadrature points [[0.61596668 0.4454805  0.5545195  0.38403332 0.82951381 0.4454805
  0.5545195  0.17048619 0.82951381 0.38403332 0.61596668 0.17048619]], detJ [[-0.25 -0.25 -0.25 -0.25 -0.25 -0.25]]
Cell 2, quadrature points [[0.4454805  0.11596668 0.38403332 0.0545195  0.4454805  0.32951381
  0.17048619 0.0545195  0.38403332 0.32951381 0.17048619 0.11596668]], detJ [[0.25 0.25 0.25 0.25 0.25 0.25]]
Cell 3, quadrature points [[0.9454805  0.61596668 0.88403332 0.5545195  0.9454805  0.82951381
  0.67048619 0.5545195  0.88403332 0.82951381 0.67048619 0.61596668]], detJ [[0.25 0.25 0.25 0.25 0.25 0.25]]
Cell 4, quadrature points [[0.11596668 0.4454805  0.0545195  0.38403332 0.32951381 0.4454805
  0.0545195  0.17048619 0.32951381 0.38403332 0.11596668 0.17048619]], detJ [[-0.25 -0.25 -0.25 -0.25 -0.25 -0.25]]
Cell 5, quadrature points [[0.61596668 0.9454805  0.5545195  0.88403332 0.82951381 0.9454805
  0.5545195  0.67048619 0.82951381 0.88403332 0.61596668 0.67048619]], detJ [[-0.25 -0.25 -0.25 -0.25 -0.25 -0.25]]
Cell 6, quadrature points [[0.4454805  0.61596668 0.38403332 0.5545195  0.4454805  0.82951381
  0.17048619 0.5545195  0.38403332 0.82951381 0.17048619 0.61596668]], detJ [[0.25 0.25 0.25 0.25 0.25 0.25]]
Cell 7, quadrature points [[0.11596668 0.9454805  0.0545195  0.88403332 0.32951381 0.9454805
  0.0545195  0.67048619 0.32951381 0.88403332 0.11596668 0.67048619]], detJ [[-0.25 -0.25 -0.25 -0.25 -0.25 -0.25]]

for quadrature degree 3. Note that if you know that your domain is affine (i.e. a first order triangle or a tetrahedra, you only need to compute detJ at one point, not every quadrature point).

For further questions it would be great if you could describe in more detail what you want to achieve with your post-processing, as there might be built-in functions for what you want already

3 Likes