Hello, I have met a problem with code for stationary Stokes problem with MINI element. I would like to use non constant external force like a function - not like a constant and got error for inner product. I don`t know how to fix it. Can anybody help? Thank You in advance.
Code (Problem with “L”):
import matplotlib.pyplot as plt
from dolfin import *
# Load mesh and subdomains
mesh = Mesh("/mnt/c/docker/fenics/documented/dolfin_fine.xml.gz")
sub_domains = MeshFunction("size_t", mesh, "/mnt/c/docker/fenics/documented/dolfin_fine_subdomains.xml.gz")
plot(mesh)
plt.plot()
# Build function spaces on Mini element
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
B = FiniteElement("Bubble", mesh.ufl_cell(), mesh.topology().dim() + 1)
V = VectorElement(NodalEnrichedElement(P1, B))
Q = P1
W = FunctionSpace(mesh, V * Q)
# No-slip boundary condition for velocity
bc0 = DirichletBC(W.sub(0), (0, 0), sub_domains, 0)
# Inflow boundary condition for velocity
inflow = Expression(("0.0", "0.0"), degree=2)
bc1 = DirichletBC(W.sub(0), inflow, sub_domains, 1)
# Collect boundary conditions
bcs = [bc0, bc1]
# Define variational problem
u, p = TrialFunctions(W)
v, q = TestFunctions(W)
f = Expression("0.0, 10*exp(1-x[1]+3*pow(x[1], 2))" , degree=3)
a = (inner(grad(u), grad(v)) - div(v)*p - q*div(u))*dx
L = inner(f, v)*dx
# Compute solution
w = Function(W)
solve(a == L, w, bcs)
# Split the mixed solution using deepcopy
# (needed for further computation on coefficient vector)
u, p = w.split(True)
print("Norm of velocity coefficient vector: %.15g" % u.vector().norm("l2"))
print("Norm of pressure coefficient vector: %.15g" % p.vector().norm("l2"))
# Split the mixed solution using a shallow copy
u, p = w.split()
# Save solution in VTK format
ufile_pvd = File("velocity.pvd")
ufile_pvd << u
pfile_pvd = File("pressure.pvd")
pfile_pvd << p
# Plot solution
plt.figure()
plot(u, title="velocity")
plt.figure()
plot(p, title="pressure")
# Display plots
plt.show()
