4.1 Finite element micromagnetics

The solution of Poisson's equation (), which is required to calculate the demagnetizing field $\mathbf{H}_\mathrm{ms}$, has to be solved for a given magnetization distribution $\mathbf{M}(\mathbf{r})$. We write Poisson's equation in a more general form

$$\Delta u(\mathbf{r}) = f(\mathbf{r}) \quad.$$ (4.1)

In order to apply the finite element method, we have to find a so-called ``weak'' or variational formulation.

But first we have to define the vector spaces, in which we are searching for a solution [33]. Given a bounded domain $\Omega$, we denote by $L^2(\Omega)$ the space of quadratically integrable functions defined on $\Omega$. The usual inner product

$$(u,v):=\int_\Omega u(\mathbf{r}) v(\mathbf{r}) \,d{\mathbf{r}}\,$$

induces the norm $\Vert u \Vert=\sqrt{(u,u)}$ and $L^2(\Omega)$ becomes a Hilbert space. The space $H^1(\Omega)$ consists of those functions in $L^2(\Omega)$, whose (weak) derivative of order one also lie in $L^2(\Omega)$. The $H^1$-inner product is defined as

$$(u,v)_1:= \int_\Omega \left( uv+\frac{\partial u}{\parti... ...\right) \,d{\mathbf{r}}\,= (u,v)+(\nabla u, \nabla v) \quad.$$

The space $H^1(\Omega)$ belongs to a family of function spaces known as Sobolev spaces. They can be physically interpreted as the space of functions of finite energy with respect to the problem under consideration. In this sense, it is the correct space in which to seek solutions of the weak formulation.

Further, we define the trial space $S$

$$S=\{ v \in H^1(\Omega): v(\mathbf{r})=g_1(\mathbf{r}) \mathrm{~on~} \Gamma\}$$

and the test space or weighting space $V$

$$V=\{ v \in H^1(\Omega): v(\mathbf{r})=0 \mathrm{~on~} \Gamma\} \quad.$$

The functions in $V$ are called test functions.

To derive the weak formulation we multiply Poisson's equation () with a test function $v(\mathbf{r})$ and integrate over the solution domain

$$\int_\Omega \Delta u(\mathbf{r}) v(\mathbf{r}) \,d{\mathbf{... ...t_\Omega f(\mathbf{r}) v(\mathbf{r}) \,d{\mathbf{r}}\, \quad.$$

Integration by parts gives

$$-\int_\Omega \nabla u(\mathbf{r}) \nabla v(\mathbf{r}) \,d{... ...nt_\Omega f(\mathbf{r}) v(\mathbf{r}) \,d{\mathbf{r}}\, \quad,$$

where $\mathbf{r}_n$ denotes the surface normal on the boundary $\Gamma$. If appropriate boundary conditions define the values of $u$ (Dirichlet boundary conditions) or of its derivatives $\nabla u =: g$ (Neumann boundary conditions) on the boundary, we can simplify (since $v$ vanishes, where Dirichlet boundary conditions apply)

$$-\int_\Omega \nabla u(\mathbf{r}) \nabla v(\mathbf{r}) \,d{... ...nt_\Omega f(\mathbf{r}) v(\mathbf{r}) \,d{\mathbf{r}}\, \quad.$$ (4.2)

The variational formulation of () can then be stated in the following general form: Find $u \in S$ such, that

$$a(u,v)=F(v) \quad \forall v \in V \quad,$$ (4.3)

where the bilinear form $a(u,v)$ and the linear functional $F(v)$ are given by

$$a(u,v)= -\int_\Omega \nabla u(\mathbf{r}) \nabla v(\mathbf... ...+ \int_\Omega f(\mathbf{r}) v(\mathbf{r}) \,d{\mathbf{r}}\, .$$ (4.4)

We can homogenize the Dirichlet boundary condition by which the trial space, in which the solution $u$ is sought, becomes equivalent to the test space $V$. Without loss of generality, we can thus assume that we seek the solution in the space $V$. It is noted, that the weak form is a generalization of the classical formulation. Therefore, the solution of the weak formulation need not be a classical solution at the same time.

Not only Poisson's equation, but a large number of boundary value problems lead to symmetric and positive definite bilinear forms. Also for () we find

$$a(u,v)=a(v,u) \quad \forall u, v$$

and

$$a(u,u)>0 \quad \forall u \not = 0 \quad.$$

In this case the weak formulation is equivalent to a minimization problem. $u$ is the sought solution, if it minimizes the functional

$$J(u)=\frac{1}{2}a(u,u)-F(u) \quad.$$

$J(u)$ can often be interpreted as an energy functional. In the context of Poisson's equation for the demagnetizing field (), this functional gives the energy of the magnetization distribution $\mathbf{M}$ in the magnetic field $\mathbf{H}$

$$J(\mathbf{H}) \geq -E(\mathbf{H}) = \mu_0 \int \mathbf{H} ... ...,d{^3r}\, +\frac{\mu_0}{2}\int \mathbf{H}^2 \,d{^3r}\, \quad.$$

Minimization of this functional with respect to $\mathbf{H}$ reduces $E$ to the stray field energy $E_\mathrm{ms}$ and makes $\mathbf{H}$ equal to the self demagnetizing field $\mathbf{H}_\mathrm{ms}$ [34].

In general, the trial and test space $V$ is too large and complex to deal with numerically. Thus, the Galerkin discretization seeks an approximation of the solution $u_h \in V_h$ by restricting it to a finite dimensional subspace $V_h$. We rewrite the weak formulation () as follows: Find $u_h \in V_h$ such, that

$$a(u_h,v)=F(v) \quad \forall v \in V_h \quad.$$ (4.5)

The exact solution $u(\mathbf{r})$ shall be approximated by a linear combination of trial functions $\phi_i(\mathbf{r})$ from a finite dimensional subspace $V_h$ of $V$

$$u_h(\mathbf{r})=\sum_{i=0}^n u_i \phi_i(\mathbf{r}) \quad.$$

If we insert this expansion in (), we obtain

$$a(u_h,v)= a \left( \sum_{i=1}^n u_i \phi_i, v \right)= ... ...}^n u_i a(\phi_i, v)= F(v_h) \quad \forall v \in V_h \quad.$$ (4.6)

Since each basis function $\phi_i$ lies in $V$, relation () holds trivially for $v=\phi_i, i \in \{1, \ldots, N\}$. Conversely, if relation () holds for each basis function $\phi_i$, then it also holds for all $v \in V_h$. Hence it is sufficient to determine the coefficients $u_i$ of $u_h$ such, that

$$\sum_{i=1}^n u_i a(\phi_i, \phi_j)=F(\phi_j), \quad j=1, \ldots, N \quad.$$ (4.7)

Therefore we arrive at a linear system of algebraic equations, which can be solved with any standard method, such as the Gauß method, by Cholesky decomposition or iterative schemes like the conjugate gradient method.

The finite element method is a particular Galerkin method [35], which uses piecewise polynomial functions to construct the finite dimensional subspace $V_h$. The solution domain is divided into many small subdomains, referred to as elements. In two space dimensions these elements are usually triangles (fig. ) or convex quadrilaterals, while in three dimensions tetrahedra, prisms and hexahedra are commonly employed. This subdivision process is usually called triangulation. The collection of all elements is referred to as the finite element mesh or grid.

Figure: Triangulation of a 2D domain

[Polygonal domain] [Triangulation]

In the finite element method the basis functions $\phi_i \in V_h$ (fig. ) are chosen in such a way, that

1. The support $\mathrm{supp}(\phi_i):=\overline{\{\mathbf{r} \in \Omega: \phi_i(\mathbf{r}) \not = 0\}}$ of each basis function, i.e. the closure of the set, where $\phi_i$ is nonzero, is small in the sense, that it consists of only a few (connected) elements. (fig. )
2. Globally, each function $v \in V_h$ has a simple description in terms of $N$ so called degrees of freedom which uniquely characterize $v$. Each basis function is characterized by possessing exactly one nonvanishing degree of freedom.

Figure: Nodal basis functions

[Nodal basis function $\phi_i$] [Nodal basis function $\phi_j$]

As the support is restricted to a very small local area, the integrals occurring in () need only be computed over the small support of each basis function. In fact, most of the integrals are zero, and so the matrix of the linear system of algebraic equations is very sparse.

Figure: Common support of two basis functions

[Support of basis function $\phi_i$] [Support of basis function $\phi_j$] [Common support of $\phi_i$ and $\phi_j$]