Hi,

I am solving for the pressure of a non-Newtonian fluid model and the boundary condition is not enforced. We specify the pressure on one side of the domain and solve, but the recovered solution is zero pretty much everywhere. I’ve tried iterative solvers with nonzero initial guesses but don’t see a difference. Maybe I’m missing something really obvious?

Thanks to anyone who gives this a look. I’ve been fiddling with this for a long time and am up for tying anything!

```
from dolfin import *
mesh = RectangleMesh(Point(-4, 0.0), Point(4, 1.0), 32, 32, "crossed")
# refine the boundary
for j in range(6):
mesh.init(1, 2)
markers = CellFunction("bool", mesh, False)
for c in cells(mesh):
for f in facets(c):
if f.exterior():
markers[c] = True
break
mesh = refine(mesh, markers)
# create function spaces
V = VectorFunctionSpace(mesh, "Lagrange", 4)
P = FunctionSpace(mesh, "Lagrange", 3)
boundary_exp_u = Expression(("6*x[1]*(1 - x[1])", "0.0"))
boundary_exp_p = Expression("-9.6*x[0]")
def p_boundary(x, on_boundary):
return on_boundary and x[0] > 4 - 3*DOLFIN_EPS
bc_u = DirichletBC(V, boundary_exp_u, "on_boundary")
bc_p = DirichletBC(P, boundary_exp_p, p_boundary)
# exact solutions
u = project(boundary_exp_u, V)
pi = project(boundary_exp_p, P)
# pressure solve
p = TrialFunction(P)
q = TestFunction(P)
a = (p*(q + 5.0*dot(u, grad(q))))*dx
L = pi*q*dx
p = Function(P)
problem = LinearVariationalProblem(a, L, p, bc_p)
solver = LinearVariationalSolver(problem)
solver.solve()
plot(p, interactive=True)
```