Cell mean of a Function used in a UFL form.

Question

Hello,

I would like to form an expression of the element mean over the space $P_1(\Omega)$, or the $L^2$-projection of a function $p \in P_1(\Omega)$ onto the discontinuous function space $ P_0(\Omega)$.

So, for each vertex $i$ of a 3D mesh located on element $e$ containing 3 other vertices $j$, $k$, and $l$, I need a function that returns

$$p_{avg}^i = \frac{p_i^e + p_j^e + p_k^e + p_l^e}{4},$$

to be used in a form with $p$ a test or trial function.

The reason is that I'm going to try and implement the pressure-projection stabilization method described here to stabilize the Stokes equations.

Answer 1

It seems that CellAvg in ufl.restriction does what you want.

Example:

from ufl.restriction import CellAvg

mesh = UnitSquareMesh(4, 4)
V = FunctionSpace(mesh, "Lagrange", 1)

u, v = TrialFunction(V), TestFunction(V)
a = CellAvg(u) * CellAvg(v) * dx

A = assemble(a)

Comments

Cool, I have to try this. I always thought there is no other way in FEniCS than the method withe the additional multiplier, that KristianE mentioned. Is this a new addition to UFL?

Answer 2

Is this example useful to you?

from dolfin import *
from mshr import *

L1 = 4.; L2 = 7.; w = 4.
noslip = Constant((0.,0.))
inflow = Expression(("1.5*(1.-pow(x[1]*2./w,2.))","0."),w=w)
r1 = Rectangle(Point(-L1,-w/2.),Point(L2,w/2.))
c1 = Circle(Point(0.,0.),1.,8)
domain = r1-c1
mesh = generate_mesh(domain, 45)

def upper_lower(x, on_boundary):
  return x[1] < -w/2.+1e3*DOLFIN_EPS or w/2.-1e3*DOLFIN_EPS < x[1]
def left(x, on_boundary):
  return x[0] < -L1+1e3*DOLFIN_EPS
def right(x, on_boundary):
  return L2 - 1e3*DOLFIN_EPS < x[0]
def cylbnd(x, on_boundary):
  return x[0]*x[0]+x[1]*x[1] < 1.5 and on_boundary

class right2(SubDomain):
  def inside(self, x, on_boundary):
    return L2 - 1e3*DOLFIN_EPS < x[0]

Q = VectorFunctionSpace(mesh, "CG", 1)
V = FunctionSpace(mesh,'CG',1)
P = FunctionSpace(mesh,'DG',0)
W = MixedFunctionSpace([Q, V, P])
U = Function(W)

#boundaries 
bc0 = DirichletBC(W.sub(0), noslip, cylbnd)
bc1 = DirichletBC(W.sub(0), noslip, upper_lower)
bc2 = DirichletBC(W.sub(1), Constant(0.), right)
bc3 = DirichletBC(W.sub(0), inflow, left)
bcs = [bc0, bc1, bc2, bc3]

(u, p, pDG) = TrialFunctions(W)
(u_test, p_test, pDG_test) = TestFunctions(W)

a = (inner(sym(grad(u)), sym(grad(u_test))) - div(u_test)*pDG \
+ p_test*div(u)  + (p-pDG)*pDG_test)*dx
L = (inner(noslip, u_test))*dx 
solve(a==L, U, bcs)
v, p, pDG = U.split()

plot(v[0], interactive=True)
plot(p, interactive=True)

Note that there is a problem related to mshr for this code, i.e. it will fail if if you change the resolution (45) or increase the number of circle segments (8).