Mapped quadrature points, weights and solution at these quadrature points

General
1

Hi,

I want to evaluate the mapped quadrature points, mapped quadrature weights i.e. quadrature weight X det(jacobian), and the value of the solution at each of the quadrature points in every element. I need this for postprocessing.

Currently, first I obtain the quadrature points and weight on the reference triangle and then perform the required affine transformation to get the mapped points and multiply the weights with det(Jacobian) to obtain the mapped weights. But is there any internal functionality within fenicsx where I can loop over elements and obtain the quadrature weights, points, and values of my approximate solution at these quadrature points?

regards
Ankit

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
3

Hello everyone,

I get an error when I execute the code, it’s about the inside of the print part

TypeError Traceback (most recent call last)
Cell In[34], line 19
16 detJ = dolfinx.fem.Expression(ufl.JacobianDeterminant(mesh), quadrature_points)
18 for i in range(mesh.topology.index_map(mesh.topology.dim).size_local):
—> 19 print(f"Cell {i}, quadrature points {x_expr.eval([i])}, detJ {detJ.eval([i])}")

I get the error TypeError: Expression.eval() missing 1 required positional argument: ‘cells’.

Isn’t this code compatible with Fenicsx and Dolfinx? If it is not compatible, how can I perform the mapping of quadrature points, weights, and solution asked in this question in Fenicsx and Dolfinx?

4

Please note that the answer to this question is over a year old, using an older version of DOLFINx, it is quite straightforward to make it compatible with 0.7.x, i.e.

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(mesh, [i])}, detJ {detJ.eval(mesh, [i])}")

gives

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]]
1 Like
5

Hi! Thank you @dokken for your answer. As I am trying to use it, I wonder when it is convenient and not convenient to add the ghost cells here:

for i in range(mesh.topology.index_map(mesh.topology.dim).size_local):

like so:

map_c = mesh.topology.index_map(mesh.topology.dim)
for i in range(map_c.size_local + map_c.num_ghosts):

I am also going to add this little note in case that it helps others (and myself) when looking for info:

Getting the coordinates of the quadrature points by element

6

This depends on what you want to use the quadrature points for, if you want to integrate over non-owned cells on every process or not.