# Construct mapping operator (jacobian) to evaluate 3D Nedelec basis functions

**URL:** <https://fenicsproject.discourse.group/t/construct-mapping-operator-jacobian-to-evaluate-3d-nedelec-basis-functions/17414>\
**Category:** dolfinx\
**Tags:** dolfinx\
**Created:** [April 10, 2025, 9:54am UTC](https://fenicsproject.discourse.group/t/construct-mapping-operator-jacobian-to-evaluate-3d-nedelec-basis-functions/17414 "2025-04-10T09:54:55Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![NicoSkibbe](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nicoskibbe/32/8246_2.png) [@NicoSkibbe](https://fenicsproject.discourse.group/u/NicoSkibbe)\
**Post date:** [April 10, 2025, 9:54am UTC](https://fenicsproject.discourse.group/t/construct-mapping-operator-jacobian-to-evaluate-3d-nedelec-basis-functions/17414/1 "2025-04-10T09:54:55Z")

</div>

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](https://fenicsproject.discourse.group/t/evaluate-basis-all-in-dolfinx/3968/12)).

1. 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?

2. 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?

3. If it is not available: I already found this guide for 2D [Advanced finite elements — FEniCS Workshop](https://jsdokken.com/FEniCS-workshop/src/finite_element_method/advanced_elements.html) 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:

* * *

```Python
#!/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:

* * *

```Python
#!/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}')

```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [April 10, 2025, 12:26pm UTC](https://fenicsproject.discourse.group/t/construct-mapping-operator-jacobian-to-evaluate-3d-nedelec-basis-functions/17414/2 "2025-04-10T12:26:03Z")

</div>

> [@NicoSkibbe](#):
>
> 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?

Yes. You would use the correct Piola transform to do this evaluation in “physical” space.

> [@NicoSkibbe](#):
>
> 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?

They are built just-in-time and discarded after computation.

> [@NicoSkibbe](#):
>
> If it is not available: I already found this guide for 2D [Advanced finite elements — FEniCS Workshop](https://jsdokken.com/FEniCS-workshop/src/finite_element_method/advanced_elements.html) to manually build the Jacobian. Is there also an example available how to build the Jacobian for 3D Nedelec spaces using the Python interface?

You can use `dolfinx.fem.Expression` to evaluate your function in any point in physical space.  
Expression takes in:

- Your expression: Say the test or trial function (or coefficient) you would like to evaluate
- Points in reference cell

Say:

```bash
v = ufl.TestFunction(V)
expr = dolfinx.fem.Expression(v, my_points_in_reference_element, comm=MPI.COMM_SELF)
expr.eval(mesh, np.array([my_cell_index, dtype=np.int32]))

```

This is for instance how `PointSource` is implemented in scifem:

> <https://github.com/scientificcomputing/scifem/blob/98dc80f776c3d6bfbb452abf67f62d59043bfacb/src/scifem/point_source.py#L22-L177>

which is based on my original post at:

> [@Point Sources (Redux)](https://fenicsproject.discourse.group/t/point-sources-redux/13496/6):
>
> Following my post above, with a minor adaptations, a point source can be implemented as follows: # Create a point source for Poisson problem # Author: Jørgen S. Dokken # SPDX-License-Identifier: MIT from mpi4py import MPI from petsc4py import PETSc import dolfinx import dolfinx.fem.petsc import numpy as np import ufl def compute\_cell\_contributions(V, points): # Determine what process owns a point and what cells it lies within mesh = V.mesh \_, \_, owning\_points, cells = dolfinx.cp…

---

<div class="post-metadata">

**Author:** ![NicoSkibbe](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nicoskibbe/32/8246_2.png) [@NicoSkibbe](https://fenicsproject.discourse.group/u/NicoSkibbe)\
**Post date:** [April 14, 2025, 3:42pm UTC](https://fenicsproject.discourse.group/t/construct-mapping-operator-jacobian-to-evaluate-3d-nedelec-basis-functions/17414/3 "2025-04-14T15:42:22Z")

</div>

Dear Dokken, dear community,

thank you for the fast answer!

I adapted your function to compute the Nedelec basis function contributions in 3D. Please consider the attached MWE where this is shown for 2 different points. I simplified the problem (serial code, small meshes) and found that my implemented algorithm to build the interpolation Operator Q (Qx, Qy, Qz) work for all points within the cells 0 or 1 of the mesh (tested varying N and also in 2D). However, it does not work for all the other cells.

I assume some kind of sorting/indexing/cell-orientation problem that I oversee? Since only one process is used, it should not be related to local and global dof ordering (as can be confirmed by viewing the local-to-global dof lists)

I’m grateful for any further advice to identifiy the issue,  
Nico Skibbe

MWE:

```Python
from mpi4py import MPI
import dolfinx
import dolfinx.fem.petsc
import numpy as np
import ufl

def compute_cell_contributions(V, points):
    """ Modified code, Original at https://fenicsproject.discourse.group/t/point-sources-redux/13496/2 """
    # Determine what process owns a point and what cells it lies within
    mesh = V.mesh
    _, _, own_points, cells = dolfinx.cpp.geometry.determine_point_ownership(
        mesh._cpp_object, points, 1e-6)
    own_points = np.asarray(own_points).reshape(-1, 3)

    # Pull owning points back to reference cell
    mesh_nodes = mesh.geometry.x
    cmap = mesh.geometry.cmap
    ref_x = np.zeros((len(cells), mesh.geometry.dim),
                     dtype=mesh.geometry.x.dtype)

    # Create expression evaluating a trial function (i.e. the basis function)
    u = ufl.TrialFunction(V)
    basis_values = []

    for i, (point, cell) in enumerate(zip(own_points, cells)):
        geom_dofs = mesh.geometry.dofmap[cell]
        ref_x[i] = cmap.pull_back(point.reshape(-1, 3), mesh_nodes[geom_dofs])
        expr = dolfinx.fem.Expression(u, ref_x[i], comm=MPI.COMM_SELF)
        values = expr.eval(mesh, np.asarray([cell], dtype=np.int32))

        # Strip out basis function values per cell
        basis_values.append(values)

    return cells, basis_values

# %% set up test problem
p = 2 # 1
nde = 20 # 6
N = 1
domain = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, N, N, N)
domain.name = "mesh"
V = dolfinx.fem.functionspace(domain, ("N1curl", p))

points = np.array([[0.7, 0.5, 0.71],
                   [0.9, 0.14, 0.2]],
                  dtype=domain.geometry.x.dtype)

# %% compute basis values
cells, basis_values = compute_cell_contributions(V, points)

# %% calc reference solution
# this is exchanged with another function later, just for testing
target_eq = lambda x: np.stack([1*x[0], 2*x[1], 3*x[2]], axis=0)
target_function = dolfinx.fem.Function(V)
target_function.interpolate(target_eq)

evaluated = target_function.eval(points, cells)

# %% use basis values to define interpolation matrix and test it
x, y, z = [], [], []
for n in range(2):
    Qx = basis_values[n][0][0:nde]
    Qy = basis_values[n][0][nde:2*nde]
    Qz = basis_values[n][0][2*nde:3*nde]

    dofs_to_eval = V.dofmap.cell_dofs(cells[n])
    u_ref_vals = target_function.x.array[dofs_to_eval]

    x.append(np.dot(Qx, u_ref_vals))
    y.append(np.dot(Qy, u_ref_vals))
    z.append(np.dot(Qz, u_ref_vals))

    validation = target_eq(points[n])
    # x, y, z supposed to be equal to point[n]
    print('######################################')
    print('Point: {} -> Cell No: {}'.format(points[n], cells[n]))
    print('dofs: ', dofs_to_eval)
    print('f(point): ', validation)
    print('eval(point): ', evaluated[n])
    print(f'dot product:[{x[n]:.4f}, {y[n]:.4f}, {z[n]:.4f}]')
    assertion = np.allclose(np.array((x[n], y[n], z[n])), validation, atol=1e-4)
    print(f'worked: {assertion}')

```

---

<div class="post-metadata">

**Author:** ![NicoSkibbe](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nicoskibbe/32/8246_2.png) [@NicoSkibbe](https://fenicsproject.discourse.group/u/NicoSkibbe)\
**Post date:** [April 17, 2025, 12:36pm UTC](https://fenicsproject.discourse.group/t/construct-mapping-operator-jacobian-to-evaluate-3d-nedelec-basis-functions/17414/4 "2025-04-17T12:36:45Z")

</div>

A quick update: Comparing the dofs and basis values against the legacy dolfin code, the signs of the ref values and the basis values seem to be different for particular entries. Thus the summation in the dot product returns the wrong value.

E.g., consider the difference between legacy dolfin and dolfinx for the same physical point (0.7, 0.3, 0.8) in a cell of a unit cube mesh:

dolfin, X direction

```auto

ref values x: [0.50 1.00 0.50 1.50 1.00 0.50]

basis values x: [-0.30 0.30 0.50 0.00 0.20 -0.20]

dot product: 0.50

```

dolfinx, X direction

```auto

ref values x: [-1.00 -0.50 -0.50 0.50 1.50 1.00]

basis values x: [-0.30 -0.50 -0.30 -0.20 0.00 0.20]

dot product: 0.80

```

dolfin, Y direction

```auto

ref values y: [0.50 1.00 0.50 1.50 1.00 0.50]

basis values y: [0.70 0.10 -0.10 0.20 -0.20 -0.00]

dot product: 0.50

```

dolfinx, Y direction

```auto

ref values y: [-1.00 -0.50 -0.50 0.50 1.50 1.00]

basis values y: [-0.10 0.10 0.70 -0.00 0.20 -0.20]

dot product: -0.20

```

dolfin, Z direction

```auto

ref values z: [0.50 1.00 0.50 1.50 1.00 0.50]

basis values z: [0.00 -0.30 -0.40 0.30 0.40 0.30]

dot product: 0.50

```

dolfinx, Z direction

```auto

ref values z: [-1.00 -0.50 -0.50 0.50 1.50 1.00]

basis values z: [0.30 0.40 0.00 0.30 0.30 0.40]

dot product: 0.50

```

Independent of the different orderings between dolfin and dolfinx, both contain identical values but sometimes with different signs. For points in cells 0 and 1, all signs are consistent. But for this specific cell (No 3 in dolfinx mesh, No 2 in dolfin mesh), the first entry for dolfin or the related thrid entry for dolfinx lead to a different sign of the resulting product. According to, e.g., [https://dl.acm.org/doi/fullHtml/10.1145/3524456](https://dl.acm.org/doi/fullHtml/10.1145/3524456), I assume the missing part is to adujst the local ordering of basis functions in each cell.

Is the local sign convention of the Nedelec basis accessible for each cell in dolfinx? I appreciate any advice!

Best regards,  
Nico

Here is our python code for the legacy dolfin:

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

N = 1
mesh = df.UnitCubeMesh(df.MPI.comm_world, N, N, N)
V = df.FunctionSpace(mesh, "N1curl", 1)

# 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.7, 0.3, 0.8]])

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(6, 3) # (n_dofs, 3)

# basis values to map from nedelec space to x-component)
print('ref values x: [{:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f}]'.format(*u.vector()[ids])) 
print('basis values x: [{:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f}]'.format(*Q[:, 0])) 
print('dot product: {:.2f}'.format(np.dot(Q[:, 0], u.vector()[ids])))

print('ref values y: [{:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f}]'.format(*u.vector()[ids])) 
print('basis values y: [{:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f}]'.format(*Q[:, 1])) 
print('dot product: {:.2f}'.format(np.dot(Q[:, 1], u.vector()[ids])))

print('ref values z: [{:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f}]'.format(*u.vector()[ids])) 
print('basis values z: [{:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f}]'.format(*Q[:, 2])) 
print('dot product: {:.2f}'.format(np.dot(Q[:, 2], u.vector()[ids])))

```

and the code for the dolfinx variant:

```python
#!/usr/bin/env python3
# -*- coding: utf-8 -*-

from mpi4py import MPI
import dolfinx
import dolfinx.fem.petsc
import numpy as np
import ufl

def compute_cell_contributions(V, points):
    """ Modified code, Original at https://fenicsproject.discourse.group/t/point-sources-redux/13496/2 """
    # Determine what process owns a point and what cells it lies within
    mesh = V.mesh
    _, _, own_points, cells = dolfinx.cpp.geometry.determine_point_ownership(
        mesh._cpp_object, points, 1e-6)
    own_points = np.asarray(own_points).reshape(-1, 3)

    # Pull owning points back to reference cell
    mesh_nodes = mesh.geometry.x
    cmap = mesh.geometry.cmap
    ref_x = np.zeros((len(cells), mesh.geometry.dim),
                     dtype=mesh.geometry.x.dtype)

    # Create expression evaluating a trial function (i.e. the basis function)
    u = ufl.TrialFunction(V)
    basis_values = []

    for i, (point, cell) in enumerate(zip(own_points, cells)):
        geom_dofs = mesh.geometry.dofmap[cell]
        ref_x[i] = cmap.pull_back(point.reshape(-1, 3), mesh_nodes[geom_dofs])
        expr = dolfinx.fem.Expression(u, ref_x[i], comm=MPI.COMM_SELF)
        values = expr.eval(mesh, np.asarray([cell], dtype=np.int32))

        # Strip out basis function values per cell
        basis_values.append(values)

    return cells, basis_values

# %% set up test problem
p = 1 # 1
nde = 6 # 6
N = 1
domain = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, N, N, N)
domain.name = "mesh"
V = dolfinx.fem.functionspace(domain, ("N1curl", p))

points = np.array([[0.7, 0.3, 0.8]],
                  dtype=domain.geometry.x.dtype)

# %% compute basis values
cells, basis_values = compute_cell_contributions(V, points)

# %% calc reference solution
# this is exchanged with another function later, just for testing
target_eq = lambda x: np.stack([1*x[0], 1*x[1], 1*x[2]], axis=0)
target_function = dolfinx.fem.Function(V)
target_function.interpolate(target_eq)

evaluated = target_function.eval(points, cells)

Q = np.zeros((nde, 3))

# %% use basis values to define interpolation matrix and test it
x, y, z = [], [], []
for n in range(1):
    Q[:, 0] = basis_values[n][0][0:nde]
    Q[:, 1] = basis_values[n][0][nde:2*nde]
    Q[:, 2] = basis_values[n][0][2*nde:3*nde]

    dofs_to_eval = V.dofmap.cell_dofs(cells[n])
    u_ref_vals = target_function.x.array[dofs_to_eval]

print('ref values x: [{:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f}]'.format(*u_ref_vals)) 
print('basis values x: [{:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f}]'.format(*Q[:, 0])) 
print('dot product: {:.2f}'.format(np.dot(Q[:, 0], u_ref_vals)))

print('ref values y: [{:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f}]'.format(*u_ref_vals)) 
print('basis values y: [{:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f}]'.format(*Q[:, 1])) 
print('dot product: {:.2f}'.format(np.dot(Q[:, 1], u_ref_vals)))

print('ref values z: [{:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f}]'.format(*u_ref_vals)) 
print('basis values z: [{:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f}]'.format(*Q[:, 2])) 
print('dot product: {:.2f}'.format(np.dot(Q[:, 2], u_ref_vals)))

```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [April 17, 2025, 12:58pm UTC](https://fenicsproject.discourse.group/t/construct-mapping-operator-jacobian-to-evaluate-3d-nedelec-basis-functions/17414/5 "2025-04-17T12:58:32Z")

</div>

Legacy dolfin orders the mesh to ensure a globally consistent orientation.  
This is not done in DOLFINx, as explained in  
[https://dl.acm.org/doi/full/10.1145/3524456](https://dl.acm.org/doi/full/10.1145/3524456)  
and to some extent in

> **[DOLFINx: The next generation FEniCS problem solving environment](https://zenodo.org/records/10447666)**
>
> DOLFINx is the next generation problem solving environment from the FEniCS Project; it provides an expressive and performant environment for solving partial differential equations using the finite element method. We present the modern...

To map values from the reference to the physical cell, one applies the push forward, as well as the corresponding dof transformation.

You can see this in

> <https://github.com/FEniCS/dolfinx/blob/2cc4dd164a6cb13fd20bd248b664d330f1be6d64/cpp/dolfinx/fem/Function.h#L655>

---

<div class="post-metadata">

**Author:** ![NicoSkibbe](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nicoskibbe/32/8246_2.png) [@NicoSkibbe](https://fenicsproject.discourse.group/u/NicoSkibbe)\
**Post date:** [April 17, 2025, 3:27pm UTC](https://fenicsproject.discourse.group/t/construct-mapping-operator-jacobian-to-evaluate-3d-nedelec-basis-functions/17414/6 "2025-04-17T15:27:43Z")

</div>

Hi,

thank you for the fast feedback and your time!

So basically we are missing the dof transformation, I see. I understand this as an identity matrix or vector somewhere that contains either 1 or -1 depending on the direction of a local dof compared to the global one.

Is this by any chance available through the basix python interface? Or any means to calculate it?

Best regards,  
Nico

---

<div class="post-metadata">

**Author:** ![daaffy](https://avatars.discourse-cdn.com/v4/letter/d/aca169/32.png) [@daaffy](https://fenicsproject.discourse.group/u/daaffy)\
**Post date:** [November 3, 2025, 10:20pm UTC](https://fenicsproject.discourse.group/t/construct-mapping-operator-jacobian-to-evaluate-3d-nedelec-basis-functions/17414/7 "2025-11-03T22:20:51Z")

</div>

Hi,

I’m also wondering if we’ve figured out how to access the dof transformation through the python interface.

I’m modelling a solid-fluid interaction problem and I need to be able to evaluate global BDM basis functions at a boundary in order to define the interface constraints. I’m trying to do this by evaluating the local basis functions through tabulation followed by a Piola transform, however, as we’ve seen above, there appears to be a missing dof transformation step. It would be great if this could be accessed in python.

Regards,  
Jack

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [November 5, 2025, 7:50pm UTC](https://fenicsproject.discourse.group/t/construct-mapping-operator-jacobian-to-evaluate-3d-nedelec-basis-functions/17414/8 "2025-11-05T19:50:41Z")

</div>

It is possible through the Python interface.  
I for instance use it in my scifem interpolation matrix: [API reference — scifem](https://scientificcomputing.github.io/scifem/docs/api.html#scifem.petsc_interpolation_matrix)  
(see: [scifem/src/scifem/interpolation.py at 06d1089e8b52d0e86445af221f1d53df1f29a9dd · scientificcomputing/scifem · GitHub](https://github.com/scientificcomputing/scifem/blob/06d1089e8b52d0e86445af221f1d53df1f29a9dd/src/scifem/interpolation.py#L110-L121))

It is also used in another package I make, namely adios4dolfinx, for reading Nedelec, BDM etc from file: [adios4dolfinx/src/adios4dolfinx/checkpointing.py at eef7270977d5378346333116b42add46338c01ec · jorgensd/adios4dolfinx · GitHub](https://github.com/jorgensd/adios4dolfinx/blob/eef7270977d5378346333116b42add46338c01ec/src/adios4dolfinx/checkpointing.py#L474-L482)