Understanding dS and friends

1

I am trying to understand, how to work with discontinuous elements. Specifically, I am interested in how dS is defined.

So here is what I imagine is the semantics of assemble(f*dS):

When performing an integral f*dS, each interior face is visited once. Fenics then does the following:

  • Chooses n, which is a normal unit vector of the face
  • The two cells adjacent to the face are labeled “+”, “-” such that n points towards + and away from -
  • The choice of n should be considered an implementation detail, that cannot be relied upon

Is that correct?

2

Here is an example that puzzles me:

from fenics import *
mesh = UnitIntervalMesh(4)
V = FunctionSpace(mesh, "DG", 0)
n = FacetNormal(mesh)
direction = Constant(("1",))

assemble(dot(n, direction) * dS)

Which throws:

Discontinuous type ReferenceNormal must be restricted.

If I choose “+” or “-” fenics accepts it and gives some results:

assemble(dot(n, direction)("+") * dS) # == 3
assemble(dot(n, direction)("-") * dS) # == -3

I am confused about this.? What does dot(n, direction)("+") mean? How is it defined on cells and why does it make a difference when coming from "+" or “-”?

3

Hi,
See this thread Integrating over an interior surface for how to specify how to restrict arguments to a specific side of the facet.

3 Likes
4

How f*dS is computed

Each interior facet is visited once and then the following happens:

  • There are two cells adjecent to the facet. One is called positive and labeled “+”, the other is called negative and labeled “-”. Which one is labeled positive should be thought of as an implementation detai. (It is decided by a cell function g. The cell on which g has the bigger value is labeled positive.)

  • Now if f involves a discontinuous quantity q, we need to decide if we want to evalutate it on the positive or on the negative cell. This is done by q("+") or q("-").

  • Similar, if we want to use a unit normal vector to the facet, we need to decide which one.
    FacetNormal(mesh)("+") gives the unit normal of the positive cell. FacetNormal(mesh)("-") gives the unit normal of the negative cell.

1 Like