How do I integrate including the facet normal with quadrilateral elements


233
views
1
6 months ago by
Nate  

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 *

for ctype in [CellType.Type.triangle, CellType.Type.quadrilateral]:

	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)
	j = assemble(dot(grad(u), n)*ds(1))
	k = assemble(dot(grad(u), Constant((0.0, 1.0)))*ds(1))

	print("CellType.%s, grad(u).n: %.3e, grad(u).(0, 1): %.3e" % (ctype, j, k))
Community: FEniCS Project
1
Thanks, an issue filed: https://bitbucket.org/fenics-project/ufl/issues/101
written 6 months ago by Jan Blechta  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »