Exchanging "evaluate_basis_all" with "tabulate" to map field vectors from a N1curl space to predefined output coordinates

dolfinx
1

Dear Fenicsx community,

first thank you for many many hours of saved time due to finding valuable answers in this forum, what a great community!

Im currently porting a project from the legacy dolfin to dolfinx.
One missing piece I could not yet reproduce is the calculation of a matrix I could get through the means of the “evaluate_basis_all”.

My intended uses are two in particular:

  1. To have a matrix that can map electrical field vectors from a arbitrary order (most likely 2) N1curl space to predefined output coordinates.
  2. To have a matrix that can take the electrical field vectors from a arbitrary order (most likely 2) N1curl space, curl them to magnetic fields and map to predefined output coordinates.

The calculation of these things I can do, but the matrix with the evaluted basis
functions at given physical coordinates is what I am missing.

In search for a possible solution i found another post with a general guide on how to do this in dolfinx (Evaluate_basis_all in dolfinx - #11 by dokken)

I followed the guide an came up with a minimal (not so much) working example, but I run out of ideas. Somehow exchanging “evaluate_basis_all” with a “tabulate” is not working as I thought, clearly I’m missing something.

Please find attached the MWE for the latest dolfinx as well as the MWE in dolfin (to have a working example). The content of the desired matrix Q is different, the rest should, apart from sorting and some sign coventions, be similar.

Any help that pushes me in the right direction is highly appreciated,

best regards,
Nico

dolfinx code:

(the one I try to fix, so the assertion at the end is True)

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import dolfinx.fem.petsc
import dolfinx as dfx
import numpy as np
from basix.ufl import element
from mpi4py import MPI

mesh = dfx.mesh.create_unit_cube(MPI.COMM_WORLD, 10, 10, 10)
point = np.array([[0.18, 0.22, 0.33]])

# %% build functionspaces
ele = element("N1curl", mesh.basix_cell(), 2)
V = dfx.fem.functionspace(mesh, ele)
u = dfx.fem.Function(V)

# this is exchanged with another function later, just for testing
f = lambda x: np.stack([x[0], x[1], x[2]], axis=0)
u_ref = dfx.fem.Function(V)
u_ref.interpolate(f)

bb_tree = dfx.geometry.bb_tree(mesh, mesh.topology.dim)
cells = []
points_on_proc = []

# Step 1/3: Find cells whose bounding-box collide with the the points
# aka from dolfin as mesh.bounding_box_tree().compute_first_entity_collision(p)
cell_candidates = dfx.geometry.compute_collisions_points(bb_tree, point)
colliding_cells = dfx.geometry.compute_colliding_cells(
    mesh, cell_candidates, point)
for i, pt in enumerate(point):
    if len(colliding_cells.links(i)) > 0:
        points_on_proc.append(pt)
        cells.append(colliding_cells.links(i)[0])

geom_dofs = mesh.geometry.dofmap[cells[0]]

# Step 2/3: pull back point to the reference element(s)
cmap = mesh.geometry.cmap
x_ref = cmap.pull_back(np.array(points_on_proc), mesh.geometry.x[geom_dofs])

# Step 3/3: Tabulate basis functions on the reference element
tabulated = ele.tabulate(0, x_ref)

# tabulate shape: (a, b, c * d)
# a: 1 or > 1 if asking for derivatives, 1 in our case
# b: number of points, 1 in our case
# c: number of dofs, 6 or 20 in our case, depending on the order
# d: value size (1 or > 1 for vectors, 3 in our case)

# tab_reshaped = tabulated.reshape(a, b, c, d)
Q = tabulated.reshape(1, 1, 20, 3)

dofs_to_eval = V.dofmap.cell_dofs(cells[0])
u_ref_vals = u_ref.x.array[dofs_to_eval]
# "u_ref.x.array[dofs_to_eval]" seem to be fitting to the legacy dolfin
# equivalent "u.vector()[ids]", apart from some sorting and sign convention
# dolfinx (this script):
# array([ 0.025,  0.035,  0.06 , -0.05 ,  0.075, -0.015])
# dolfin (the other, working script):
# array([0.025, 0.06 , 0.035, 0.075, 0.05 , 0.015])

x = np.dot(Q[0, 0, :, 0], u_ref_vals)
y = np.dot(Q[0, 0, :, 1], u_ref_vals)
z = np.dot(Q[0, 0, :, 2], u_ref_vals)
# calc reference solution
reference = u_ref.eval(point, cells)
print('reference: ', reference)
# x, y, z supposed to be equal to point[0]
print(f'dot product:[{x:.2f}, {y:.2f}, {z:.2f}]')
assertion = np.allclose(np.array((x, y, z)), reference, atol=1e-4)
print(f'worked: {assertion}')

legacy dolfin code:

(MWE as reference only)

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import dolfin as df
import numpy as np

mesh = df.UnitCubeMesh(df.MPI.comm_world, 10, 10, 10)
V = df.FunctionSpace(mesh, "N1curl", 2)

# func to eval (test to map on itself)
ex = df.Expression(('x[0]', 'x[1]', 'x[2]'), degree=2)
u = df.Function(V)
u.interpolate(ex)

# input
point = np.array([[0.18, 0.22, 0.33]])

coords = mesh.coordinates()
dolfin_element = V.dolfin_element()
bbt = mesh.bounding_box_tree()
dofmap = V.dofmap()

for ni in range(point.shape[0]):
    # Loop over all interpolation points
    xi = point[ni, :]
    p = df.Point(xi[0], xi[1], xi[2])
    # Find cell for the point
    cell_id = bbt.compute_first_entity_collision(p)
    # Vertex coordinates for the cell

    if cell_id < len(mesh.cells()):
        xvert = coords[mesh.cells()[cell_id, :], :]
        cell = df.Cell(mesh, cell_id)
        v = dolfin_element.evaluate_basis_all(xi, xvert,
                                              cell.orientation())
        vx = v[::3]
        vy = v[1::3]
        vz = v[2::3]
        ids = dofmap.cell_dofs(cell_id)
        Q = v.reshape(20, 3)  # (n_dofs, 3)

x = np.dot(Q[:, 0], u.vector()[ids])
y = np.dot(Q[:, 1], u.vector()[ids])
z = np.dot(Q[:, 2], u.vector()[ids])
print('reference: ', point[0])
print(f'dot product:[{x:.2f} {y:.2f} {z:.2f}]')
assertion = np.allclose(np.array((x, y, z)), point[0], atol=1e-4)
print(f'worked: {assertion}')