How to create a Dolfinx expression from a cross product?

dolfinx
4

Looks like the cross product between the normal vector to the mesh, and a 2D vector is handled with ufl.perp()

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
cross1024×768 183 KB

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.