3.1.2 Hybrid Finite elements-Boundary elements method

In this section we briefly describe a numerical method for the evaluation of the demagnetizing field which is commonly used when spatial discretization based on Finite Elements (FE) Method is performed [52]. The main advantage of FE discretization lies in the possibility of simulating processes occurring in bodies with in principle arbitrary geometry. This is crucial for the design of technological devices, whose geometry is often very far from being approximated with a 'staircase'. Nevertheless, as mentioned in section 3.1, plain FE method would require that the whole space was discretized. For this reason, many researchers have dealt with the derivation of FE modifications to take into account the open-boundary nature of the problem. In fact, in order to impose the regularity condition at infinity, Chen suggests that the FE mesh has to be extended over a large region outside the magnetic particles (at least five times the extension of the particle [54]). Various other techniques have been proposed to reduce the size of the external mesh or to avoid a discretization of the exterior space. The use of asymptotic boundary conditions [55] reduces the size of the external mesh as compared to truncation. At the external boundary, Robbin conditions, which are derived from a series expansion of the solution of the Laplace equation for outside the magnet and give the decay rate of the potential at a certain distance from the sample, are applied [56]. A similar technique that considerably reduces the size of the external mesh is the use of space transformations to evaluate the integral over the exterior space. Among the various transformations proposed to treat the open boundary problem, the parallelepipedic shell transformation [57], which maps the external space into shells enclosing the parallelepipedic interior domain, has proved to be most suitable in micromagnetic calculations. The method can be easily incorporated into standard FE programs transforming the derivatives of the nodal shape functions. This method was applied in static three-dimensional micromagnetic simulations of the magnetic properties of nanocrystalline permanent magnets (see Refs. [58] and [59]). An alternative approach was proposed by Fredkin and Koehler [70] in 1990. The main idea, due to the linearity of Poisson problem, is to split the scalar potential $\varphi$ into $\varphi_1$ and $\varphi_2$ such that:

$$\varphi=\varphi_1+\varphi_2 \quad.$$ (3.13)

The boundary value problem for the potential $\varphi_1$ can be formulated as an internal Neumann problem in the following way:

$$\left\{\begin{aligned} &\nabla^2\varphi_1=\nabla\cdot\textbf{M... ...l \Omega}=\textbf{n}\cdot\textbf{M} \end{aligned} \right. \quad.$$

By comparing the problem (3.14) with the original problem (3.3), we can derive the boundary value problem for $\varphi_2$:

$$\left\{\begin{aligned} &\nabla^2\varphi_2=0 &\text{in }&\Omega... ... at infinity }\partial\Omega_\infty \end{aligned} \right. \quad,$$

where $\partial \Omega^-$ indicates the internal layer of the boundary $\partial \Omega$. The solution $\varphi_2$ of the latter boundary value problem is the well-known double layer potential [13]:

$$\varphi_2(\textbf{r})=\frac{1}{4\pi}\int_{\partial\Omega} \varph... ...(\frac{1}{\vert\textbf{r}-\textbf{r}'\vert} \right)\cdot\textbf{n} dS \quad.$$ (3.16)

Until now, it is not yet evident the advantage of the method, since the evaluation of the potential $\varphi_2(\textbf{r})$ at each location $\textbf{r}\in \Omega$ is very expensive from computational point of view. Indeed the situation seems more complicated with respect to the direct evaluation of the original integral (3.7), since now $N\times N_b$ operations ($N_b$ is the number of boundary nodes) plus the solution of boundary value problem (3.14) are required in order to obtain $\varphi$. The nice idea is to use Eq. (3.16) to evaluate the potential $\varphi_2$ only on the boundary $\partial \Omega^-$. In fact, it is known from potential theory [66] that if $\partial \Omega$ is sufficiently smooth at location $\textbf{r}_0$, then

$$\underset{\substack{\textbf{r}\rightarrow\textbf{r}_0^- \text... ...frac{1}{\vert\textbf{r}_0-\textbf{r}'\vert} \right)\cdot\textbf{n} dS \quad.$$ (3.17)

It can be shown [67] that the discretized version of the latter equation is:

$$\underline{\varphi}_2=B\underline{\varphi}_1 \quad,$$ (3.18)

where $B$ is a suitable $N_b\times N_b$ boundary matrix and $\underline{\varphi}_1,\underline{\varphi}_2$ are the vector containing the boundary nodal values of the scalar potentials $\varphi_1,\varphi_2$. From the knowledge of $\varphi_2$ on the boundary $\partial \Omega$, which now costs $N_b^2$ operations, then the boundary value problem (3.15), now with Dirichlet boundary conditions on $\partial \Omega$, can be solved with the usual FE technique [52]. Let us summarize the costs of this hybrid technique. We assume that the mesh (typically consisted of tetrahedrons) has $N$ nodes. In most situations, one can think that the boundary nodes $N_b$ are in the order $\mathcal{O}(N^{2/3})$. In particular, this happens when characteristic dimensions of the body, along the coordinates axes, are of the same order of magnitude. We will discuss remarkable exceptions after. Therefore,

• at the preprocessing stage, the boundary matrix $B [N_b\times N_b]$ has to be computed and stored. The storage requirement is $\mathcal{O}(N_b^2)$.
• In a single evaluation of the demagnetizing field:
  1. the internal Neumann problem (3.14) has to be solved with the usual FE technique to find $\varphi_1$. This implies that it costs the inversion of an $N \times N$ linear system.
  2. The value of $\varphi_2$ on the boundary has to be computed by means of Eq. (3.18), which costs $N_b^2$ operations.
  3. The internal Dirichlet problem, given by Eq. (3.15) with boundary conditions obtained at step 2, has to be solved with FE technique, which costs another $N \times N$ linear system inversion, but with same stiffness matrix used at step 1.
  4. Finally, the demagnetizing field ${\mathbf{H}_\text{m}}=-\nabla(\varphi_1+\varphi_2)$ has to be evaluated.

We observe that this method has the advantage to manage arbitrarily complicated geometries, whereas the FFT convolution method is optimal with a structured mesh. Nevertheless, if one deals with somehow ``flat'' magnetic bodies, like for instance thin-films, the assumption $N_b\thicksim\mathcal{O}(N^{2/3})$ fails, since it happens instead $N\simeq N_b$, unless than one performs additional (and maybe useless) refinement along the ``small'' dimension. Thus, one can conclude that the hybrid FE-BE method is not optimal for thin-film geometries.

In the sequel, this method will be used in full micromagnetic simulations of damping and precessional switching processes to investigate how far they can be treated within the framework of the uniform mode approximation analyzed in chapter 2. The FFT convolution method described in section 3.1.1 will be used in chapter 4 where the solution of micromagnetic standard problem no. 4 [79] will be addressed.