Getting exact function values at guass points (query points)

dolfinx
1

Dear FEniCSx community,

I have solutions over Lagrange functionspace (nodal dofs). However, I want to query the solutions at guass points.

Using boxbeam tree- gives a constant value for all points lying in an element. (one value per element).

Can you suggest functionspace which uses guass points as dofs (like how lagrange uses nodes) or way to obtain interpolated value at query point (guass point for my case)?



bb_tree = dolfinx.geometry.bb_tree(mesh, mesh.topology.dim)

points = np.array(points, dtype=mesh.geometry.x.dtype)
potential_colliding_cells = dolfinx.geometry.compute_collisions_points(
    bb_tree, points
)
colliding_cells = dolfinx.geometry.compute_colliding_cells(
    mesh, potential_colliding_cells, points
)
points_on_proc = []
cells = []
for i, point in enumerate(points):
    if len(colliding_cells.links(i)) > 0:
        points_on_proc.append(point)
        cells.append(colliding_cells.links(i)[0])
points_on_proc = np.array(points_on_proc, dtype=np.float64).reshape(-1, 3)
cells = np.array(cells, dtype=np.int32)

I remember in legacy dolfin, we can directly give query to solved ufl-function with any coordinate point- to obtain function value.


u=Function(V)  # solved ufl function
u.vector(2,3,0.5) # extracting function value at (x,y,z) coordinate--> (2,3,0.5).
2

This is surely not right. However, as you haven’t provided a reproducible example, it is not really easy to give you guidance.

You can use a “quadrature space” for this, as illustrated below:

from mpi4py import MPI
import dolfinx
import basix.ufl

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 2, 2)
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 2))
u = dolfinx.fem.Function(V)
u.interpolate(lambda x: x[0] * x[1])

q_el = basix.ufl.quadrature_element(mesh.basix_cell(), degree=2, scheme="default")
Q = dolfinx.fem.functionspace(mesh, q_el)

q = dolfinx.fem.Function(Q)
u_expr = dolfinx.fem.Expression(u, Q.element.interpolation_points)
q.interpolate(u_expr)

q_coords = Q.tabulate_dof_coordinates()
for coord, val in zip(q_coords, q.x.array):
    print(f"q({coord}) = {val}")
1 Like