Why doesn't a simple pressure expression using FacetNormal give the expected analytical results?

variational formulation
1

Hi everyone,

I am trying to implement a pressure field p, which has the following form (simplified here for my minimal working example): p = \underline{\lambda} \cdot \underline{n}; where \underline{\lambda} is a known vector (equal for example here to [0.,1.,0.]). My pressure field is only defined on the boundary of my geometry.

I would like to enforce this pressure value via a variational formulation (I need to use a variational formulation because I cannot have directly an analytical expression of p in my real problem, which is a little bit more complex):
\int_\gamma (p-\underline{\lambda} \cdot \underline{n})p^* dS = 0 \, \forall p^*

Here is my minimal working code:

import dolfin

mesh = dolfin.UnitCubeMesh(10, 10, 10)

n = dolfin.FacetNormal(mesh)

p_elts = dolfin.FiniteElement("CG", mesh.ufl_cell(), 1)

W = dolfin.FunctionSpace(mesh, p_elts)

boundaries = dolfin.MeshFunction("size_t", mesh, mesh.topology().dim()-1)

internal = dolfin.CompiledSubDomain(" 1.0 > x[0] && x[0]>0. && 1.0 > x[1] && x[1]> 0. && 1.0 > x[2] && x[2]> 0.")

bc_p = dolfin.DirichletBC(W, 0., internal, method="pointwise")

p = dolfin.TrialFunction(W)

p_test = dolfin.TestFunction(W)

lbda = dolfin.Constant([0., 1., 0.])

dV = dolfin.Measure("dx", domain=mesh, subdomain_data=boundaries)

dS = dolfin.Measure( "exterior_facet", domain=mesh, subdomain_data=boundaries)

F = dolfin.Constant(0.)*dolfin.inner(p, p_test)*dV + (p-dolfin.dot(lbda, n)) * p_test * dS

a = dolfin.lhs(F)

L = dolfin.rhs(F)

w = dolfin.Function(W)

dolfin.solve(a==L, w, bc_p)

p = w

file_pvd = dolfin.File('test_pressure_field'+".pvd")

file_pvd << (p)

I tested the formulation on a cube. Given the form of my pressure field p, I would expect to have p=1 at the face y=1, p=-1 at the face y=-1 and 0 on all other faces. However, it is not the case. From the discussion in Output solution across surface of mesh with element coordinates to a .txt file - #10 by dokken, I understood that the function FacetNormal is interpolated at the nodes, which (at least, this is my guess) may cause the differences I observe with the expected solution.

Is my guess correct? If it is the case, does anyone have an idea of how I can get the expected analytical results (i.e. p=+/-1 at the faces y=+/-1 and 0 on all the other faces)?

Thank you in advance for your help!