Open Boundary Problems

Multiple effective-field terms, namely the demagnetization field (DemagField) and the Oersted field (OerstedField) involve the solution of so called open-boundary problems of the form

\[\begin{split}\Delta u &= \nabla \cdot \vec{f} \quad \text{in} \quad \omega \\ \Delta u &= 0 \quad \text{in} \quad \mathbb{R}^3 \setminus \omega \\ [u] &= 0 \quad \text{on} \quad \partial\omega \\ \left[\frac{\partial u}{\partial \vec{n}}\right] &= - \vec{f} \cdot \vec{n} \quad \text{on} \quad \partial\omega \\ u(\vec{x}) &= \mathcal{O}(1/|\vec{x}|) \quad \text{if} \quad |\vec{x}| \rightarrow \infty\end{split}\]

different algorithms are available for the solution of these problems.

Shell-Transformation Method

The shell-transformation method for the solution of open-boundary problems with the finite-element method was introduced in [Brunotte1992]. The region of interest \(\omega_\text{notrans}\) is surrounded by a finite shell \(\omega_\text{trans}\). A bijective transformation is constructed that maps the shell region \(\omega_\text{trans}\) onto the complete exterior region \(\mathbb{R}^3 \setminus \omega_\text{notrans}\). The open-boundary conditions can thus be implemented by requiring homogeneous Dirichlet boundary condition of the outer boundary of the shell \(\omega_\text{trans}\), resulting in the following weak formulation

\[\int_{\omega_\text{notrans}} \nabla u \cdot \nabla v \dx + \int_{\omega_\text{trans}} (\nabla u)^T \mat{g} \nabla v \dx = \int_{\omega} \vec{f} \cdot \nabla v \dx\]

with the metric tensor \(\mat{g}\) defined as

\[\mat{g} = (\mat{J}^{-1})^T | \det \mat{J} | \mat{J}^{-1}\]

where \(\mat{J}\) denotes the Jacobian of the transformation. The untransformed region \(\omega_\text{notrans}\) does not have to coincide with the sample region \(\omega\) of the original problem, but it holds \(\omega \subseteq \omega_\text{notrans}\).

magnum.fe uses cuboid shells, i.e. \(\omega_\text{notrans}\) has to be of cuboidal shape. The applied method is described in detail in [Abert2013a].

Hybrid FEM-BEM Method (Fredkin and Koehler)

A hybrid FEM-BEM approach for the solution of open-boundary problems was proposed by Fredkin and Koehler [Fredkin1990]. Consider the following splitting of the solution \(u\):

\[u = u_1 + u_2\]

\(u_1\) is defined by

\[\begin{split}\Delta u_1 &= \nabla \cdot \vec{f} \quad \text{in} \quad \omega \\ \frac{\partial u_1}{\partial \vec{n}} &= - \vec{f} \cdot \vec{n} \quad \text{on} \quad \partial\omega \\ u_1 &= 0 \quad \text{in} \quad \mathbb{R}^3 \setminus \omega.\end{split}\]

This Neuman problem within \(\omega\) is solved with the finite-element method. While \(u_1\) solves for the right-hand side \(\nabla \cdot \vec{f}\) and fullfills the jump condition of the normal derivative \(- \vec{f} \cdot \vec{n}\) it is not continuous across \(\partial \omega\). This jump is compensated by \(u_2\) which is defined as

\[\begin{split}\Delta u_2 &= 0 \quad \text{in} \quad \omega \\ [u_2] &= - [u_1] \quad \text{on} \quad \partial \omega \\ \left[\frac{\partial u_2}{\partial \vec{n}}\right] &= 0 \quad \text{on} \quad \partial\omega \\ u_2(\vec{x}) &= \mathcal{O}(1/|\vec{x}|) \quad \text{if} \quad |\vec{x}| \rightarrow \infty\end{split}\]

This system is solved by the double-layer potential

\[u_2 = \int_{\partial \omega} u_1 \frac{\partial}{\partial \vec{n}} \frac{1}{|\vec{x} - \vec{x}'|} \dx\]

The double-layer potential is computed on the boundary \(\partial \omega\) with the boundary-element method. This can be done either by numeric integration of a Galerkin-BEM approach (bem_method = "BEM++") or by analytic integration of the Collocation-BEM operators (bem_method = "H2LIB") [Lindholm1984]. The collocation method uses \(H^2\)-matrices provided by H2Lib, whereas the Galerkin method uses \(H\)-matrices provided by BEM++.

These values are used as Dirichlet boundary conditions to solve \(u_2\) within \(\omega\) with the finite-element method. For a detailed description of the method, see [Fredkin1990].

Direct FEM-BEM Method (Johnson and Nedelec)

The solution of (magnetostatic) Maxwell equations requires to handle more complicated jump conditions than the ones assumed for the methods above.

A general open-boundary problem with arbitrary jump-conditions may be solved by the method proposed by Johnson and Nedelec [Johnson1980]. Compared with the methods presented above it is able to handle arbitrary (scaled) boundary jump conditions for the potential \(u\) as well as for the normal derivative \(\phi = \frac{\partial u}{\partial \vec{n}}\) The current implementation assumes non-scaled boundary jump conditions for the potential \(u\)) whereas the jump of the normal derivative may be scaled by a function \(\mu\):

\[\begin{split}L u &= f \quad \text{in} \quad \omega \\ L u &= 0 \quad \text{in} \quad \mathbb{R}^3 \setminus \omega \\ [u] &= u_0 \quad \text{on} \quad \partial\omega \\ \left[\mu \, \phi \right] &= \phi_0 \quad \text{on} \quad \partial\omega \\ u(\vec{x}) &= \mathcal{O}(1/|\vec{x}|) \quad \text{if} \quad |\vec{x}| \rightarrow \infty\end{split}\]

The problem is solved by combining FEM and BEM equations and solve them simultaneously by means of an iterative method like GMRES. The BEM part is based on the external representation formula derived from Green’s second identity, whereas the FEM part directly follows from a Galerkin FEM approach, where the Neumann boundary condition are left as unknowns. The combined Galerkin formulation of the upper problem results in the following system of equations

\[\begin{split}a(u, v) - \left\langle \phi, v \right\rangle_\Gamma &= \left\langle f, v\right\rangle_\Omega + \left\langle\phi_0, v\right\rangle_\Gamma \\ \left\langle \left(\frac{1}{2} - K \right) u, \psi \right\rangle_\Gamma + \left\langle V \, \phi, \psi \right\rangle_\Gamma &= \left\langle \left(\frac{1}{2} - K \right) u_0, \psi \right\rangle_\Gamma\end{split}\]

with the weak form \(a(u,v)\) of the differential operator \(L\), the unknown potential \(u\) and (outer) normal derivative \(\phi\). The scalar product \(\langle u , v \rangle_D\) stands for the integral of the function \(u\) weighted by the test function \(v\) over a domain \(D\). The occuring single-layer operator \(V\) as well as the double-layer operator \(K\) are defined as

\[\begin{split}V \phi(x) &= \int_\Gamma G(x,x') \, \phi(x') \, d\Gamma(x') \\ K u(x) &= \int_\Gamma \partial_{n(x')} G(x,x') \, u(x') \, d\Gamma(x')\end{split}\]

References:

[Brunotte1992]

Brunotte, X., Meunier, G., & Imhoff, J. F. (1992). Finite element modeling of unbounded problems using transformations: a rigorous, powerful and easy solution. IEEE Transactions on Magnetics, 28(2), 1663-1666.

[Abert2013a]

Abert, C., Exl, L., Selke, G., Drews, A., & Schrefl, T. (2013). Numerical methods for the stray-field calculation: A comparison of recently developed algorithms. Journal of Magnetism and Magnetic Materials, 326, 176-185.

[Fredkin1990]

(1, 2) Fredkin, D. R., & Koehler, T. R. (1990). Hybrid method for computing demagnetizing fields. IEEE Transactions on Magnetics, 26(2), 415-417.

[Lindholm1984]

Lindholm, D. A. (1984). Three-Dimensional Magnetostatic Fields from Point-Matched Integral Equations with Linearly Varying Scalar Sources. IEEE Transactions on Magnetics, 20(5), 2025-2032.

[Johnson1980]

Johnson, C., & Nedelec, J. C. (1980). On the Coupling of Boundary Integral and Finite Element Methods, Mathematics of Computation, 35(152).