### How do I integrate including the facet normal with quadrilateral elements

233
views
1
6 months ago by

Let  $\Omega$Ω be the unit square. I wish to compute the flux of a finite element function,   $u\in V^{h,p}$uVh,p  , on the top boundary,
$\int_{\partial\Omega_{\text{top}}}\nabla u\cdot n\mathrm{\mathrm{\quad d}}s$∂Ωtopu·n ds  where  $\partial\Omega_{\text{top}}:=\left[0,1\right]\times\left\{1\right\}$∂Ωtop:=[0,1]×{1}. In my understanding, this is equivalent to computing the functional  $\int_{\partial\Omega_{\text{top}}}\nabla u\cdot\left(0,1\right)^{\top}\quad\mathrm{d}s$∂Ωtopu·(0,1) ds . I am experiencing unexpected results, am I doing something wrong?

Using a triangle mesh, the results are as expected.

Using a quad mesh yields results which do not match that expected when using the facet normal.

Is this issue related to https://bitbucket.org/fenics-project/ffc/pull-requests/91/fail-on-compilation-of-dq-element/diff ?

Example code follows:

from dolfin import *

mesh = RectangleMesh.create(MPI.comm_world,
[Point(0, 0), Point(1, 1)],
[16, 16], ctype)

V = FunctionSpace(mesh, "Lagrange", 1)
u = interpolate(Expression("pow(x[0], 2) + pow(x[1], 2)", element=V.ufl_element()), V)

# Mark the top edge
ff = MeshFunction("size_t", mesh, 1, 0)
CompiledSubDomain("near(x[1], 1.0)").mark(ff, 1)
ds = Measure("ds", subdomain_data=ff)

# Compute functionals
n = FacetNormal(mesh)