Dear Fenicsx community,
While porting a project from dolfin to dolfinx, I’m still struggling to build an explicit matrix to map from Nedelec function spaces to points inside cells. I need this operator for electromagentic geophysical inversion purposes.
I managed to build such an operator in dolfinx for Lagrange spaces by using the “tabulate” method (working example attached). However, this is not yet working to map from Nedelec spaces due to the missing transformation.
From my current understanding, the “evaluate_basis_all” functionality in legacy dolfin included the transformations (Jacobian), but the “tabulate” method from dolfinx does not (as mentioned, e.g., here: Evaluate_basis_all in dolfinx - #12 by reckx_sanchez).
-
To be sure: Is combining the transformation matrix (Jacobian) with the outcome of the “tabulate” method the correct way to evaluate, e.g., Nedelec basis functions, in dolfinx?
-
The “eval” method should require the same functionalities. Are the transformaion matrices build internally just-in-time or are they stored somethere and accessible using the Python interface?
-
If it is not available: I already found this guide for 2D Advanced finite elements — FEniCS Workshop to manually build the Jacobian. Is there also an example available how to build the Jacobian for 3D Nedelec spaces using the Python interface?
I’m grateful for any advice,
Nico Skibbe
Working code for Lagrange spaces:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
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, 1, 1, 1)
# random point
point = np.array([[0.1, 0.2, 0.3]])
# build functionspaces
n = 1
ele = element("Lagrange", mesh.basix_cell(), n, shape=(3,))
V = dfx.fem.functionspace(mesh, ele)
# define function
def func(x):
return np.array([3*x[0], x[1] + 1, x[1] + x[2]])
u_ref = dfx.fem.Function(V)
u_ref.interpolate(func)
# Step 1/4: Find cells whose bounding-box collide with the the points
# aka from dolfin as mesh.bounding_box_tree().compute_first_entity_collision(p)
bb_tree = dfx.geometry.bb_tree(mesh, mesh.topology.dim)
cells = []
points_on_proc = []
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/4: 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/4: Tabulate basis functions on the reference element
tabulated = ele.tabulate(0, x_ref)
# Q is the matrix we are interested in
Q = tabulated.reshape(1, 1, 3, -1)
# Step 4/4: test
# figure out the values to multiply with Q and do the multiplication
dofs_to_eval = V.dofmap.cell_dofs(cells[0])
u_ref_vals = u_ref.x.array.reshape(-1, 3)
x = np.dot(Q[0, 0, 0, ::3], u_ref_vals[dofs_to_eval, 0])
y = np.dot(Q[0, 0, 1, 1::3], u_ref_vals[dofs_to_eval, 1])
z = np.dot(Q[0, 0, 2, 2::3], u_ref_vals[dofs_to_eval, 2])
# calc reference solution
reference = u_ref.eval(point, cells)
print('reference: ', reference)
# x, y, z supposed to be equal to reference
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 (reference == dot product): {assertion}')
Modified code for Nedelec spaces:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
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, 1, 1, 1)
# random point
point = np.array([[0.1, 0.2, 0.3]])
# build functionspaces
n = 1
ele = element("N1curl", mesh.basix_cell(), n)
V = dfx.fem.functionspace(mesh, ele)
# define function
def func(x):
return np.array([3*x[0], x[1] + 1, x[1] + x[2]])
u_ref = dfx.fem.Function(V)
u_ref.interpolate(func)
# Step 1/4: Find cells whose bounding-box collide with the the points
# aka from dolfin as mesh.bounding_box_tree().compute_first_entity_collision(p)
bb_tree = dfx.geometry.bb_tree(mesh, mesh.topology.dim)
cells = []
points_on_proc = []
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/4: 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/4: Tabulate basis functions on the reference element
tabulated = ele.tabulate(0, x_ref)
# !!! Missing part to get from reference basis function to the physical
# basis function as described in the tutorial link above
# Q is the matrix we are interested in
Q = tabulated.reshape(1, 1, 3, -1)
# Attention:
# Whats missing here, according to the advanced element workshop code, is the
# Jacobian of the mapping from the reference element to the physical element
# at a single point in the reference element.
# But still if continuing the reference solution using "eval" should produce
# the correct output, no?
# Step 4/4: test
# figure out the values to multiply with Q and do the multiplication
dofs_to_eval = V.dofmap.cell_dofs(cells[0])
u_ref_vals = u_ref.x.array
x = np.dot(Q[0, 0, 0], u_ref_vals[dofs_to_eval])
y = np.dot(Q[0, 0, 1], u_ref_vals[dofs_to_eval])
z = np.dot(Q[0, 0, 2], u_ref_vals[dofs_to_eval])
# calc reference solution
reference = u_ref.eval(point, cells)
print('reference: ', reference)
print('reference should be: [0.3 1.2 0.5]')
# x, y, z supposed to be equal to reference
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 (reference == dot product): {assertion}')