# How to create a Dolfinx expression from a cross product?

**URL:** <https://fenicsproject.discourse.group/t/how-to-create-a-dolfinx-expression-from-a-cross-product/18132>\
**Category:** dolfinx\
**Tags:** dolfinx\
**Created:** [August 20, 2025, 9:51pm UTC](https://fenicsproject.discourse.group/t/how-to-create-a-dolfinx-expression-from-a-cross-product/18132 "2025-08-20T21:51:42Z")\
**Posts on this page:** 1\
**Showing post:** 4

<div class="post-metadata">

**Author:** ![fea](https://avatars.discourse-cdn.com/v4/letter/f/e79b87/32.png) [@fea](https://fenicsproject.discourse.group/u/fea)\
**Post date:** [August 23, 2025, 12:23am UTC](https://fenicsproject.discourse.group/t/how-to-create-a-dolfinx-expression-from-a-cross-product/18132/4 "2025-08-23T00:23:58Z")

</div>

[Looks like the cross product between the normal vector to the mesh, and a 2D vector is handled with ufl.perp()](https://fenicsproject.discourse.group/t/implementing-coriolis-term-in-shallow-water-equations/15191/2)

```auto
from mpi4py import MPI
from petsc4py import PETSc

import numpy as np

from basix.ufl import element, mixed_element
from dolfinx import fem, io, mesh, plot
from dolfinx.fem.petsc import LinearProblem
from ufl import (Measure, SpatialCoordinate, TestFunctions, TrialFunctions,
                 div, exp, inner, as_vector, cross, CellNormal, FacetNormal)
import ufl

import pyvista as pv

domain = mesh.create_unit_square(MPI.COMM_WORLD, 32, 32, mesh.CellType.quadrilateral)

k = 1
Q_el = element("BDMCF", domain.basix_cell(), k)
P_el = element("DG", domain.basix_cell(), k - 1)
V_el = mixed_element([Q_el, P_el])
V = fem.functionspace(domain, V_el)

Ve = fem.functionspace(domain, ("Lagrange", k + 1))

(sigma, u) = TrialFunctions(V)
(tau, v) = TestFunctions(V)

x = SpatialCoordinate(domain)
f = 10.0 * exp(-((x[0] - 0.5) * (x[0] - 0.5) + (x[1] - 0.5) * (x[1] - 0.5)) / 0.02)

dx = Measure("dx", domain)
a = inner(sigma, tau) * dx + inner(u, div(tau)) * dx + inner(div(sigma), v) * dx
L = -inner(f, v) * dx

fdim = domain.topology.dim - 1
facets_top = mesh.locate_entities_boundary(domain, fdim, lambda x: np.isclose(x[1], 1.0))
Q, _ = V.sub(0).collapse()
dofs_top = fem.locate_dofs_topological((V.sub(0), Q), fdim, facets_top)

def f1(x):
    values = np.zeros((2, x.shape[1]))
    values[1, :] = np.sin(5 * x[0])
    return values

f_h1 = fem.Function(Q)
f_h1.interpolate(f1)
bc_top = fem.dirichletbc(f_h1, dofs_top, V.sub(0))

facets_bottom = mesh.locate_entities_boundary(domain, fdim, lambda x: np.isclose(x[1], 0.0))
dofs_bottom = fem.locate_dofs_topological((V.sub(0), Q), fdim, facets_bottom)

def f2(x):
    values = np.zeros((2, x.shape[1]))
    values[1, :] = -np.sin(5 * x[0])
    return values

f_h2 = fem.Function(Q)
f_h2.interpolate(f2)
bc_bottom = fem.dirichletbc(f_h2, dofs_bottom, V.sub(0))

bcs = [bc_top, bc_bottom]

problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu",
                                                      'petsc_options_prefix': 'mixed_poisson',
                                                      "pc_factor_mat_solver_type": "mumps"},
                                                       petsc_options_prefix="mixed_poisson_",
                                                        kind="mpi")

w_h = problem.solve()

sigma_h, u_h = w_h.split()

sigma_v = as_vector(sigma_h)
sigma_cross = ufl.perp(sigma_v) #FIXED: Use ufl.perp() instead

W = fem.functionspace(domain, ("DG", 0, (domain.geometry.dim, )))
fem.Expression(sigma_cross, W.element.interpolation_points)

W = fem.functionspace(domain, ("DG", 0, (domain.geometry.dim, )))

sigma_v = ufl.as_vector(sigma_h)
sigma_d = fem.Function(W)
expr_d = fem.Expression(sigma_cross,W.element.interpolation_points)
sigma_d.interpolate(expr_d)

u_d = fem.Function(W)
expr_d = fem.Expression(sigma_h,W.element.interpolation_points)
u_d.interpolate(expr_d)

num_dofs = W.dofmap.index_map.size_local
num_verts = W.mesh.topology.index_map(0).size_local

values_d = np.zeros((num_dofs, 3), dtype=np.float64)
values_d[:, :domain.geometry.dim] = sigma_d.x.array.real.reshape(num_dofs, W.dofmap.index_map_bs)
values = np.zeros((num_dofs, 3), dtype=np.float64)
values[:, :domain.geometry.dim] = u_d.x.array.real.reshape(num_dofs, W.dofmap.index_map_bs)

grid = pv.UnstructuredGrid(*plot.vtk_mesh(Ve))
grid["sigma"] = values_d
glyphs = grid.glyph(orient = "sigma", scale = "sigma", factor = 0.1)# tolerance=0.1)
grid.set_active_vectors("sigma")
grid["u"] = values
grid.set_active_scalars("u")

u_plotter = pv.Plotter()
u_plotter.add_mesh(grid,scalars="u",show_edges= False)
u_plotter.add_mesh(glyphs, cmap = 'bwr')#, clim = [np.min(grid["sigma"]),np.max(grid["sigma"])])
u_plotter.view_xy()
u_plotter.show()

```

 ![cross](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/1/1504c23a7d22ba302627d425903dbc4bf0dce6b6.png)

Next, I need to determine if this works when setting up a Linear Problem in KSP or PETSc. Once I have that, I can feel more confident in how this is handled in the 2D case. The Shallow Water Equation in 2D like in the link might be a good one to test.

---

_[View the full topic](https://fenicsproject.discourse.group/t/how-to-create-a-dolfinx-expression-from-a-cross-product/18132)._