### Normal component of a vector function on interior facets

194
views
0
12 weeks ago by
Given a vector function   $u\in\left[V\right]^2$u[V]2 , I want to find the projection of its trace on and perpendicular to the interior facets,  $u\cdot n_F$u·nF  , $\forall F\in\Gamma_{int}$FΓint . Rather naively I have constructed the following minimal example:

from dolfin import *

mesh = UnitSquareMesh(2, 2)
n = FacetNormal(mesh)

V = VectorFunctionSpace(mesh, "CG", 1)
u = interpolate(Constant((1, 1)), V)

Vscl = FunctionSpace(mesh, "CG", 1)
uscl, vscl = TrialFunction(Vscl), TestFunction(Vscl)
projection_eqn = (uscl*vscl*dx == (dot(u, n)*vscl)("+")*dS)

normal_pressure = Function(Vscl)
solve(projection_eqn, normal_pressure)​

Would anyone have any advice or more rigorous methods of computing such a function?

Community: FEniCS Project

6
12 weeks ago by
Discontinuous Lagrange Trace is the finite element family you are looking for. It is a function space living just on interior and exterior facets and being discontinuous piecewise polynomial.

from dolfin import *

mesh = UnitSquareMesh(1, 1)
n = FacetNormal(mesh)

V = VectorFunctionSpace(mesh, "CG", 1)
u = interpolate(Constant((1, 1)), V)

Vscl = FunctionSpace(mesh, "Discontinuous Lagrange Trace", 0)
u_n, v = TrialFunction(Vscl), TestFunction(Vscl)

normal_trace = Function(Vscl)

projection_eqn = ((u_n * v)("+") * dS + (u_n * v) * ds ==
(dot(u, n) * v)("+") * dS + (dot(u, n) * v) * ds)
solve(projection_eqn, normal_trace)
print(normal_trace.vector()[:])
​

Mathematical note: In infinite dimensional setting, the most general space with well defined normal trace is H(\text{div}, \Omega) having trace in H^{-1/2}(\partial \Omega). These H^{-1/2} are ugly, not even functions, and defined via Green's theorem. Fortunately, in finite dimension we work with H(\text{div}, \Omega) conforming subspaces and they have normal traces piecewise polynomial discontinuous.
For more see e.g.: "A Simple Introduction to the Mixed Finite Element Method" from Gatica.

Another tricky way would be to interpolate into Raviart-Thomas (degree=1) space and scale its degrees of freedom by facets area. You would get just normal component evaluations.

written 12 weeks ago by Nate
1
Thanks Michal, very neat.
When looking at the output vector, do you know if there is a simple way to associate which value belongs to which facet in the mesh.
You cannot use eval on this FunctionSpace, so I was thinking of just iterating over cells:

facet_to_dofmap = {}

for c in cells(mesh):
i = c.index()
fac = c.entities(1)
dofs = Vscl.dofmap().cell_dofs(i)

for f,d in zip(fac, dofs):
facet_to_dofmap[f] = d

print(facet_to_dofmap)
​
written 11 weeks ago by Chris Richardson
1
Tormod did some plotting of DGT0, and you can see that his logic is the same,  see https://bitbucket.org/trlandet/ocellaris/src/08c0249e8dd827e6dab763d74cf74817b7e77900/ocellaris/utils/dofmap.py

But it won't work for higher order spaces, there are more dofs per facet.
written 11 weeks ago by Michal Habera