This should work, it is Arnold's example. The only things that have changed since then are:
- You have to explicitly set the Expression's degree
- You must create Mixed Spaces with a MixedElement (at least that's how it works for me)
from dolfin import *
from numpy import log
# forcing function, high degree for better approximation
f = Expression("exp(-10.*(pow(x[0], 2) + pow(x[1]-1/sqrt(2.), 2) + pow(x[2]-1/sqrt(2.), 2)))",
degree = 4)
mesh = Mesh("hemisphere_mesh.xml")
global_normal = Expression(("x[0]", "x[1]", "x[2]"), degree = 1)
mesh.init_cell_orientations(global_normal)
# finite element mesh and spaces
V0 = FunctionSpace(mesh, "RT", 3)
V1 = FunctionSpace(mesh, "DG", 2)
R = FunctionSpace(mesh, "Real", 0)
Mel = MixedElement([V0.ufl_element(), V1.ufl_element(), R.ufl_element()])
V = FunctionSpace(mesh, Mel)
# essential boundary condition
def bdry(x, on_boundary):
return on_boundary
bc = DirichletBC(V.sub(0), Constant((0., 0., 0.0)), bdry)
# trial and test functions
(u, p, c0) = TrialFunctions(V)
(v, q, c1) = TestFunctions(V)
# bilinear form, right hand side
a = (dot(u, v) - p*div(v) - div(u)*q + c0*q + p*c1) * dx
L = -f*q * dx
# assemble, solve
up = Function(V)
solve(a == L, up, bc)
(u, p, c0) = up.split()
I didn't know about setting the normal, so thanks for the question. Let me know in case of anything. Best regards!