How to extract vector components and magnitude at a specific point in FEniCSx

General
1

I have defined a simple vector function in a domain. I would like to know how to extract the components of the vector function at a specific point, as well as calculate its magnitude (Euclidean norm) at that point.
Here is MWE:

import dolfinx
from mpi4py import MPI
import ufl
from ufl import grad
from dolfinx import fem

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 2, 2)

element = ufl.VectorElement("Lagrange", mesh.ufl_cell(), 1)
V = fem.FunctionSpace(mesh, element)

u = fem.Function(V)
x = ufl.SpatialCoordinate(mesh)
u.interpolate(fem.Expression(grad(x[0]**2+2*x[1]**2),
                                                  V.element.interpolation_points()))
print(u.x.array)

The output is:

[1. 0. 2. 0. 2. 2. 1. 2. 0. 0. 2. 4. 0. 2. 1. 4. 0. 4.]

It seems that the degrees of freedom mixed the two components together. What should be done to extract its components, even its norm? Moreover, how to obtain the value of any point？

Thank you for your help!

2

What I can think of right now is to compute the norm of a vector, one possible approach is to define a function for the norm and then interpolate it.
i.e.

import dolfinx
from mpi4py import MPI
import ufl
from ufl import grad, sqrt, dot
from dolfinx import fem

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 2, 2)

element = ufl.VectorElement("Lagrange", mesh.ufl_cell(), 1)
V = fem.FunctionSpace(mesh, element)

u = fem.Function(V)
x = ufl.SpatialCoordinate(mesh)
u.interpolate(fem.Expression(grad(x[0]**2+2*x[1]**2),
                                                  V.element.interpolation_points()))
print(u.x.array)

N = fem.functionspace(mesh, ("Lagrange", 1))
u_norm = fem.Function(N)
u_norm.interpolate(fem.Expression(sqrt(dot(u, u)), N.element.interpolation_points()))

print(u_norm.x.array)

However, I’m uncertain if this is the most efficient way to compute the vector norm in terms of both memory and computational cost. Additionally, extracting individual vector components is still unclear to me.

3

Is the point going to align with a degree of freedom, or can it be at an arbitrary position in your domain?

You can consider something along the lines of:


from mpi4py import MPI
import dolfinx
import ufl
import numpy as np
comm = MPI.COMM_WORLD
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)


def determine_ownership_and_pull_back(mesh, points):
    tol = float(1e2 * np.finfo(points.dtype).eps)
    if dolfinx.__version__ == "0.8.0":
        _, _, recv_points, cells = dolfinx.cpp.geometry.determine_point_ownership(
            mesh._cpp_object, points, tol
        )
        recv_points = np.array(points).reshape(-1, 3)
    elif dolfinx.__version__ == "0.9.0.0":
        collision_data = dolfinx.cpp.geometry.determine_point_ownership(
            mesh._cpp_object, points, tol
        )
        recv_points = collision_data.dest_points
        cells = collision_data.dest_cells
    else:
        raise RuntimeError("Unsupported version of dolfinx")
    mesh_nodes = mesh.geometry.x
    cmap = mesh.geometry.cmap

    ref_x = np.zeros((len(cells), mesh.topology.dim), dtype=mesh.geometry.x.dtype)
    for i, (point, cell) in enumerate(zip(recv_points, cells)):
        geom_dofs = mesh.geometry.dofmap[cell]
        ref_x[i] = cmap.pull_back(point.reshape(-1, 3), mesh_nodes[geom_dofs])
    return cells, ref_x


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

if MPI.COMM_WORLD.rank == 0:
    points = np.array([[0.5, 0.5,0], [0.3, 0.2,0]], dtype=np.float64)
else:
    points = np.empty((0, 3), dtype=np.float64)
cells, ref_x = determine_ownership_and_pull_back(mesh, points)

grad_u = ufl.grad(u)
mag_grad = ufl.sqrt(ufl.inner(grad_u, grad_u))
u_comp = grad_u[0]


for ufl_expr in [grad_u, mag_grad, u_comp]:
    for cell, point in zip(cells, ref_x):

        expr = dolfinx.fem.Expression(ufl_expr,point, comm=MPI.COMM_SELF)
        print(str(ufl_expr), cell, point, expr.eval(mesh, np.array([cell], dtype=np.int32)))

if you just want to do it at a small subset of points.

To get individual components you can do u.sub(i).collapse() to get that single component.