Questions regarding the FEM modeling of the outer layer

Malte Bernhard Höltershinken, Takfarinas Medani

Two questions:

First, the question asks how the outer layer is determined.

Second, the question asks how boundary conditions are assigned.

(a) How is the outer layer determined? In the FEM code, the volume conductor is determined from

(1): A list of nodes, more concretely, a list of 3d positions.

(2): A list of elements. Here, an element is given by the nodes it consists of.

(3): A conductivity for each element

Now, each element has a set of associated faces. If we, for example, have a tetrahedral mesh, each element consists of 4 nodes and has 4 triangular faces. If an element now consists of the nodes (n_1, n_2, n_3, n_4), we know that its 4 triangular faces are given by the node-triples (n_1, n_2, n_3), (n_1, n_2, n_4), (n_1, n_3, n_4), and (n_2, n_3, n_4). Now, a face is a boundary face if, and only if, it is contained in exactly one element. The boundary is then given as the union over all boundary faces.

Note that the determination of the boundary is independent of the conductivity assigned to the elements, and thus does not rely on the notion of an "outer layer." You do not need to assume a layered structure of the head (you could, for example, have skull holes or you could assign a different conductivity tensor to each white matter element).

(b) How are boundary conditions assigned?

In FEM, one typically differentiates between so-called "essential" and "natural" boundary conditions. The first step towards a finite element discretization consists of going from the so-called "strong formulation" to the so-called "weak formulation." The strong formulation is the partial differential equation you want to solve, e.g.

• div(\sigma \nabla u) = div(j^P) inside the head

  <\sigma \nabla u, \eta> = 0 on the boundary of the head

A weak formulation now consists of rewriting this in the form Find u in U such that a(u, v) = l(v) for all v in V, where a is some (oftentimes bilinear) form, l is a linear functional, and U and V are suitable spaces of functions. A boundary condition is now called "essential," if it is explicitly enforced in the space U. A boundary condition is furthermore called "natural," if it is implicitly present in the functional l and fulfilling a(u, v) = l(v) for all v in V will enforce its validity.

Typical examples of essential boundary conditions are so-called Dirichlet boundary conditions, where values for the solution on the boundary are prescribed. If we, for example, prescribe the value 0 on the boundary, the standard approach would be to choose as U a space of functions that vanish on the boundary. For FEM approaches, this means that we would have to explicitly identify the degrees of freedom corresponding to the boundary values and set their values.

Typical examples of natural boundary conditions are so-called Neumann boundary conditions, where values for the normal derivative of the solution on the boundary are prescribed. Here, the space U also contains functions that do not fulfill the boundary condition. But, at least in the continuous formulation, fulfilling a(u, v) = l(v) for all v in V enforces that the solution u fulfills the boundary condition. This can typically be derived by using the divergence theorem to link volume and surface integrals. Importantly, in FEM approaches one does not need to explicitly identify the boundary degrees of freedom and assign prescribed values. Note that the two paragraphs above refer to the "standard" case. There are, for example, FEM approaches that do not treat Dirichlet boundary conditions as essential boundary conditions.