Normal component of a vector function on interior facets

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

1 Answer

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)

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.

Excellent. Thanks for your tips.
written 12 weeks ago by Nate  
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

written 11 weeks ago by Chris Richardson  
Tormod did some plotting of DGT0, and you can see that his logic is the same,  see

But it won't work for higher order spaces, there are more dofs per facet.
written 11 weeks ago by Michal Habera  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »