For the variational form of the (inhomogeneous) shallow water equations, I need to find
u^\bot = \textbf{k}\times \textbf{u}
where \textbf{k} is the unit vector pointing normal to the surface. I tried to use ufl.FacetNormal taking inspiration from this example to create the unit vector normal to the mesh k = ufl.FacetNormal(mesh), and from this, I thought simply running u_n = ufl.cross(k,u) would suffice. Running the following compiled
from mpi4py import MPI
import meshio
import numpy as np
import ufl
import basix.ufl
import dolfinx as df
import scifem
from petsc4py import PETSc
from dolfinx.fem.petsc import apply_lifting, assemble_matrix, assemble_vector, set_bc
from Projector import Projector
order = 1
exact_case = 2
tmp_data = meshio.dolfin.read("supplement/examples/mixed-poisson/hdiv-l2/meshes/sphere_ico4.xml")
mesh = df.mesh.create_mesh(
comm=MPI.COMM_WORLD,
cells=tmp_data.cells_dict["triangle"],
x=tmp_data.points,
e=ufl.Mesh(basix.ufl.element("Lagrange", "triangle", 1, shape=(3,))),
)
def create_function_spaces(mesh, family, order):
S = df.fem.functionspace(mesh, (family, order))
V = df.fem.functionspace(mesh, ("DG", order-1))
PlottingMesh = df.fem.functionspace(mesh, ("CG", order))
return (S, V, PlottingMesh)
def initial_conditions(S, V):
u0 = df.fem.Function(S)
u0.interpolate(lambda x: np.zeros(x[:S.mesh.topology.dim].shape))
f = ufl.exp(-((-x[1]-1)/0.1)**2)
petsc_options = {"ksp_type": "preonly", "pc_type": "lu"}
V_projector = Projector(V, petsc_options=petsc_options)
D0 = V_projector.project(f)
return (u0, D0)
# Define global normal
x = ufl.SpatialCoordinate(mesh)
global_normal = x # Not sure exactly what this is used for
family = "RT"
order = 1
# Define mixed mesh
(S,V,PlottingMesh) = create_function_spaces(mesh,family,order)
W = ufl.MixedFunctionSpace(*(S, V))
# Constants
dt = df.default_scalar_type(0.05)
H = 1.0
g = 1.0
# forcing term
f = x[2]
# Initial values of solutions u and D
u_, D_ = initial_conditions(S,V)
# normal unit vector
k = ufl.FacetNormal(mesh)
# u_n = k x u
u_n = ufl.cross(k,u)
F = (ufl.inner(u - u_, w) - dt*ufl.div(w)*g*D_mid + dt*f*ufl.inner(u_n, w)
+ ufl.inner(D - D_, phi) + dt*H*ufl.div(u_mid)*phi)*ufl.dx
(a, L) = ufl.system(F)
where the mesh is a 2D spherical manifold of radius 1 and the specific mesh used is found in the following file (make sure to extract the .gz file containing the mesh), and the code for the projector can be found here.
However, the code breaks when attempting to setup the problem
petsc_options = {
"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_mat_solver_type": "mumps",
"ksp_error_if_not_converged": True,
}
problem = df.fem.petsc.LinearProblem(
ufl.extract_blocks(a),
L,
bcs=[],
petsc_options=petsc_options,
petsc_options_prefix="mixed_poisson_",
kind="mpi",
)
with error
ValueError: Integral of type cell cannot contain a ReferenceNormal.
which I presume to mean that the cross product with k being in the variational form is causing the error. Are there things I can do to get around this, or should ufl.FacetNormal generally be avoided for what I’m trying to do?