# Topology Optimization for Magnetic Circuits with Adjoint Method in 3D

Zakaria Houta, Frédéric Messine, Thomas Huguet

<figure><img src="image_1.png" /></figure>

To cite this version:

Zakaria Houta, Frédéric Messine, Thomas Huguet. Topology Optimization for Magnetic Circuits with Adjoint Method in 3D. 2023. hal-04095069v2

HAL Id: hal-04095069

https://hal.science/hal-04095069v2

Preprint submitted on 3 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Topology Optimization for Magnetic Circuits with Continuous Adjoint Method in 3D

Zakaria Houta, Frédéric Messine, Thomas Huguet

► To cite this version:

Zakaria Houta, Frédéric Messine, Thomas Huguet. Topology Optimization for Magnetic Circuits with Continuous Adjoint Method in 3D. COMPEL: The International Journal for Computation and Mathematics in Electrical and Electronic Engineering, 2023, 43 (2), pp.370-389. 10.1108/COMPEL-12-2023-0644. hal-04095069v3

HAL Id: hal-04095069

https://hal.science/hal-04095069v3

Submitted on 12 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Topology Optimization for Magnetic Circuits with Adjoint Method in 3D

Zakaria Houta,* houta@laplace.univ-tlse.fr

Frédéric Messine,† messine@laplace.univ-tlse.fr

Thomas Huguet,‡ huguet@laplace.univ-tlse.fr

October 3, 2023

## Abstract

In this paper a 3D-topology optimization method in electromagnetism is introduced to design 3D-magnetic circuits. Some new theoretical results are provided to deal with 3D design cases in order to compute the derivatives using the adjoint method. Moreover, an homothopy algorithmic optimization approach is developed based on SIMP to enforce the convergence to discrete solutions close to be manufacturable. This new approach is applied and validated to design a 3D-magnetic circuit with close to 20000 variables.

**Keywords:** Topology Optimization, Adjoint Method, Sensitivity Analysis, Magnetostatic inverse problem

## 1 Introduction

Topology optimization was introduced since 2000 to design 2D mechanical structures, [2, 5, 17, 13]; as for example, to design the inside mechanical structures of airplane wings, [1]. Among different topology optimization approaches [18], the main interesting ones to solve large problems are based on continuous relaxations of the binary variables which are used to show if a cell is air or material. Thus, local descent based algorithms can be efficiently used to solve these large problems. During the optimization process, the intermediate points have to be penalized in order to enforce local optimization solvers to converge to binary solutions representing for example air or materials. Bendsoe et al. have introduced different approaches to perform SIMP method [5, 4, 15, 9]. Note SIMP method is largely used in professional solvers, as for example OptiStruct of Altair company. Because local search solvers are based on gradient descent steps, the computations of the derivative has to be very efficient. Thus adjoint methods were developed to provide the gradient of a function which depends on expensive simulation computations, [2, 5, 17].

More recently, new codes were implemented to solve 3D design problems in mechanic mainly based on voxels (boxes which are extensions of the 2D-pixel). The most famous result was obtained by Sigmud et al. in [1] who designed a 3D airplane wing for Boeing with a billion of voxels. This result was obtained in 5 days using a parallel code based on SIMP, adjoint method and local search algorithm.

<footnote>* LAPLACE-CNRS, Toulouse-INP, ENSEEIHT, University of Toulouse</footnote>

<footnote>† LAPLACE-CNRS, Toulouse-INP, ENSEEIHT, University of Toulouse and Deeper Pulse, Toulouse</footnote>

<footnote>‡ LAPLACE-CNRS, Toulouse-INP, ENSEEIHT, University of Toulouse</footnote>
# Topology Optimization for Magnetic Circuits with
Continuous Adjoint Method in 3D

November 12, 2024

## Abstract

* Purpose: The aim of this article is to present a new approach to optimizing the design of 3D magnetic circuits. This approach is based on topology optimization, where derivative calculations are performed using the continuous adjoint method. Thus, the continuous adjoint method for magnetostatics has to be developed in 3D and has to be combined with penalization, filtering and homotopy approaches in order to provide an efficient optimization code.

* Design/methodology/approach: In order to provide this new topology optimization code, we start from 2D-magnetostatic results to perform the sensitivity analysis, and this approach is extended to 3D. From this sensitivity analysis, the continuous adjoint method is derived to compute the gradient of an objective function of a 3D topological optimization design problem. From this result, this design problem is discretized and can then be solved by finite element software. Thus, by adding the Solid Isotropic Material with Penalization (SIMP) penalization approach and developing a homotopy-based optimization algorithm, we provide an interesting means for designing 3D magnetic circuits.

* Findings: In this paper, we present the 3D continuous adjoint method for magnetostatic problems involving an objective least-squares function. Based on 2D results, we provide new theoretical results for developing sensitivity analysis in 3D taking into account different parameters including the ferromagnetic material, the current density and the magnetization. Then, by discretizing, filtering and penalizing using SIMP approaches, a topology optimization code has been derived to address only the ferromagnetic material parameters. Based on this efficient gradient computation method, we develop a homotopy-based optimization algorithm for solving large-scale 3D design problems.

* Originality/value: In this paper, we propose an approach based on topology optimization to solve 3D magnetostatic design problems when an objective least-squares function is involved. This approach is based on the continuous adjoint method we have derived for 3D magnetostatic design problems. The effectiveness of this topology optimization code is demonstrated by solving the design of a 3D magnetic circuit with up to 100,000 design variables.

**Keywords:** Topology Optimization, Adjoint Method, Sensitivity Analysis, Magnetostatic Inverse Problem.

## 1 Introduction

Topology optimization was first introduced to design 2D mechanical structures, (Allaire 2007, Bendsøe & Sigmund 2003, Sigmund 2007); as for example, to design the inside mechanical
From 2010, some studies were done to solve topology optimization design problems in elec-
tromagnetism, dealing with Maxwell equations. To solve these kinds of problems, the adjoint
method has also been developed, [16, 21, 9]. In [14], magnetic circuits were designed in such a
way. More recently the nonlinearity of ferromagnetic material has been taken into account in [20]
and [9]. Moreover, some new material schemes also called penalization schemes were introduced
for this SIMP method dedicated to problems in electromagnetism, [15, 9].

As far as we know, topology optimization using adjoint methods in electromagnetism has
only been developed for 2D design problems. In this paper, we extend the topology optimization
with adjoint method for 3D design problems in electromagnetism with the assumption that the
ferromagnetic materials are linear. Some extensions to the nonlinear case could be quite directly
done from [20, 9]. Moreover, topology optimization problems in electromagnetism using multi-
materials are addressed in [9, 6].

In Section 2, topology optimization with adjoint method is recalled for 2D design problems
in electromagnetism. Section 3 is dedicated to the extension to 3D design problems with a new
theoretical formulation of the gradient. In Section 4, we detail our algorithmic approach to solve
topology optimization problems using SIMP techniques. Section 5 is dedicated to the gradient
validation (by comparing it to the finite difference technique) and dedicated to numerical tests to
design a 3D magnetic circuit. Section 6 concludes.

**2 Problem description**

In this work, Ω represents the whole studied domain. Within Ω, Ω<sub>V</sub> represents the variable domain in which the matter can be added or removed. This matter can be permanent magnet in Ω<sub>V</sub>(M), coils with current in Ω<sub>V</sub>(J) or ferromagnetic material in Ω<sub>V</sub>(ν<sub>r</sub>). Ω<sub>T</sub> represents the target domain in which a particular magnetic flux density B<sub>0</sub> has to be imposed playing with the matter in Ω<sub>V</sub>. Fig. 1 shows these different domains. In this work Ω<sub>V</sub>(M), Ω<sub>V</sub>(J) and Ω<sub>V</sub>(ν<sub>r</sub>) have to be disconnected.

<figure><img src="image_3.png" /><figcaption>Figure 1: Subdivision of the domain for topology optimization: Ω<sub>T</sub> is the target area and Ω<sub>V</sub>(ν<sub>r</sub>), Ω<sub>V</sub>(M), Ω<sub>V</sub>(J) are the design variable areas.</figcaption></figure>

The topology optimization problem is defined as follows:

* The objective function is the squared error in a target zone Ω_T between an imposed magnetic flux density B_0 and the magnetic flux density B obtained with respect to the magnetic sources and the material distributions in all the domain Ω. Note that B has to be computed by solving magnetostatic Maxwell equations.

* The design variables are magnetic sources (active material) such as ferromagnetic material, coil with current density and permanent magnets. The associated variables for all the variable domains Ω_V are the distribution of relative reluctivity ν_r, current density J (in the 3 directions (x, y, z)) and magnetization M (in the 3 directions (x, y, z)). Hence, we define a
structures of airplane wings, (Aage, Andreassen, Lazarov & Sigmund 2017). Among differ-
ent topology optimization approaches, the main interesting ones for solving large problems, are
based on continuous relaxations of the binary variables which are used to show whether a cell
is air or material. Thus, local descent based algorithms can be efficiently used to solve these
large-scale problems. During the optimization process, the intermediate points have to be pe-
nalized in order to enforce local optimization solvers to converge to binary solutions represent-
ing for example air or material. Bendsoe et al. have introduced the well-known Solid Isotropic
Material with Penalization (SIMP) approach and various penalization schemes can be found in
(Bendsøe & Sigmund 2003, Bendsøe & Sigmund 1999, Sanogo & Messine 2018, Gauthey, Gangl
& Hassan 2021). It is worth noting that SIMP-based methods are widely used in professional
solvers, as for example OptiStruct of Altair company. Because local search solvers are based
on gradient descent steps, the computations of the derivative has to be very efficient. Thus adjoint
methods were developed to provide the gradient of a function which depends on expensive com-
putations based numerical simulations, (Allaire 2007, Bendsøe & Sigmund 2003, Sigmund 2007).

More recently, new codes were implemented to solve 3D design problems in mechanics mainly based on voxels (boxes which are extensions of the 2D-pixel). The most famous result was obtained by Aage et al. in (Aage et al. 2017) who designed a 3D airplane wing for Boeing with a billion of voxels. This result was obtained in 5 days using a parallel code based on SIMP, adjoint method and a local search algorithm.

From approaches in mechanics, some studies were done to solve topology optimization design problems in electromagnetism dealing with Maxwell equations. Specially in magnetostatics, to solve these kinds of problems, the adjoint method has also been developed, (Sanogo, Messine, Henaux & Vilamot 2014, Youness & Messine 2019b, Gauthey et al. 2021, Gauthey, Hassan, Messine & Gillon 2023). In (Sanogo 2016), magnetic circuits were designed in such a way. More recently the nonlinearity of ferromagnetic material has been taken into account in (Youness & Messine 2019a) and (Gauthey et al. 2021). Moreover, some new material schemes also called penalization schemes were introduced for this SIMP method dedicated to problems in magnetostatics, (Sanogo & Messine 2018, Gauthey et al. 2021). A recent survey paper present and discuss about all the conventional approaches developed in order to solve topology optimization design problems in electromagnetism, (Lucchini, Torchio, Cirimele, Alotto & Paolo 2022). This paper emphasizes on the intrinsic advantages and limits of all these approaches, (Lucchini et al. 2022).

As far as we know, topology optimization using continuous adjoint methods in magneto-
statics has only been developed for 2D design problems. Nevertheless, note that in (Seebacher,
Kaltenbacher, Wein & Lehmann 2021), a first recent study to design 3D-magnetic circuit is done
but using a discrete adjoint method associated with SIMP approach. This adjoint method needs to
first discretize the problem by a finite element methods and then, the adjoint problem can be quite
directly derived. In the continuous approach that we developed in this work, the adjoint method is
directly applied to the continuous formulation of the problem and therefore, it is not directly linked
to the used discretization numerical method. Some other works deal also with electromagnetic de-
sign in 3D but they are based on topological derivative methods which are derived from shape op-
timization, (Masmoudi, Pommier & Samet 2005, Gangl & Sturm 2021). In this paper, we extend
the topology optimization with adjoint method for 3D design problems in magnetostatics with the
assumption that the ferromagnetic materials are linear. Note that, this could not be a big issue
to extend the 3D-adjoint method including the non-linearity of ferromagnetic materials as it is
shown in recent works for the 2D case, (Youness & Messine 2019a, Gauthey et al. 2021). Indeed,
instead to solve a linear system, we will have to solve a nonlinear one and as shown in (Gauthey
et al. 2021), that will multiply by about 7 the total CPU-times. Moreover, another limitation of
this paper is that we focus on topology optimization to address design problems of a unique ma-
terial. This also could be extended following recent works in 2D, (Gauthey et al. 2021, Cherrière,
Laurent, Hlioui, Louf, Duysinx, Geuzaine, Ben Ahmed, Gabsi & Fernández 2022).

In Section 2, topology optimization with continuous adjoint method is recalled for 2D design problems in magnetostatics. Section 3 is dedicated to the extension to 3D design problems with a new theoretical formulation of the gradient based on the continuous adjoint method. In Section 4,
global variable **p** as:

$$ \mathbf{p} : \left| \begin{array}{l} \Omega_V \quad \rightarrow \quad \mathbb{R}^7 \\ (x, y, z)^T \quad \mapsto \quad \mathbf{p}(x, y, z) = (p_1, \dots, p_7)^T \end{array} \right. $$

Where $\mathbf{p} = (\nu_r, J_x, J_y, J_z, M_x, M_y, M_z)^T$ and we define the three following functions: $\nu_r(\mathbf{p}) = \nu_r, \mathbf{J}(\mathbf{p}) = (J_x, J_y, J_z)^T$ and $\mathbf{M}(\mathbf{p}) = (M_x, M_y, M_z)^T$.

Note that, **p** is defined in the whole domain $\Omega$, **p** is fixed on $\Omega \setminus \Omega_V$ and variable on $\Omega_V$.

Therefore, the optimization problem we addressed is the following:

$$ \min_{\mathbf{p}} F_o(\mathbf{p}) = \int_{\Omega_T} \| \mathbf{B}(\mathbf{A}(\mathbf{p})) - \mathbf{B}_0 \|^2 d\Omega, \quad (1) $$

where **A** corresponds to the magnetic potential vector obtained by solving the Maxwell equations.

In this work, the constitutive law of the magnetic materials is supposed to be linear in order to simplify the computations in 3D; some extensions to the non-linear case can be done from [9].

Hence, considering the linear case, we have:

$$ \mathbf{H} = \nu \mathbf{B} - \mathbf{M} \quad (2) $$

Where $\nu = \nu_0\nu_r$ is the magnetic reluctivity of the material, $\nu_0$ is the air reluctivity and $\nu_r$ is the relative reluctivity of the material.

In different points of domain $\Omega$, we have for **p**:

* For air: $\nu_r = 1, \mathbf{M} = \mathbf{0}$ and $\mathbf{J} = \mathbf{0}$ in the 3 directions.

* For ferromagnetic material: $\nu_r$ is a constant depending on the material (in this paper, $\nu_r = \frac{1}{1000}$), $\mathbf{M} = \mathbf{0}$ and $\mathbf{J} = \mathbf{0}$ in the 3 directions.

* For permanent magnets: $\nu_r = 1, \mathbf{M}$ is a constant bounded vector in 3D and $\mathbf{J} = \mathbf{0}$ in the 3 directions.

* For coils: $\nu_r = 1, \mathbf{M} = \mathbf{0}$ and $\mathbf{J}$ is a constant bounded vector in 3D.

Because **p** is extended on the whole domain $\Omega$, in $\Omega \setminus \Omega_V$ (with $\Omega_V = \Omega_V(\nu_r) \cup \Omega_V(\mathbf{M}) \cup \Omega_V(\mathbf{J})$), **p** is fixed and cannot change during the optimization process. Thus, in $\Omega \setminus \Omega_V$, **p** can be air, a ferromagnetic material, a permanent magnet or a coil. In the variable domain $\Omega_V$, **p** is variable and can be air or a ferromagnetic material if $(x, y, z)^T$ belongs to $\Omega_V(\nu_r)$, air or a permanent magnet if $(x, y, z)^T$ belongs to $\Omega_V(\mathbf{M})$ and air or a coil if $(x, y, z)^T$ belongs to $\Omega_V(\mathbf{J})$.

The Maxwell-Ampere equation $\nabla \times \mathbf{H} = \mathbf{J}$ becomes, using the constitutive law of the material described in equation (2) and the existence of the potential vector **A** such that $\mathbf{B} = \nabla \times \mathbf{A}$:

$$ -\nabla \times [\nu_0 \nu_r(\mathbf{p}) \mathbf{B}(\mathbf{A}(\mathbf{p})) - \mathbf{M}(\mathbf{p})] + \mathbf{J}(\mathbf{p}) = 0, \text{ in } \Omega. \quad (3) $$

Equation (3) is called, in this paper, the primal problem. Note that $\mathbf{B}(\cdot)$ is the rotational operator $\nabla \times (\cdot)$ which is also denoted by $\mathbf{B}(\mathbf{A}) = \text{curl}(\mathbf{A})$. The Dirichlet condition $\mathbf{A} \times \mathbf{n} = \mathbf{0}$ is imposed on the domain boundary $\partial\Omega$, where $\mathbf{n}$ is the normal vector (this hypothesis is relevant because the border of the domain $\Omega$ is far from the zone of interest $\Omega_T$), as shown in Fig. 1.

The partial derivative $\frac{\partial F_o}{\partial p_i}$ is calculated using the adjoint method. Hence, from equations (1) and (3), the optimization problem (1) can be equivalently formulated as:

$$ (\mathcal{P}) \left\{ \begin{array}{ll} \min_{\mathbf{p},\mathbf{A}} & F(\mathbf{p}, \mathbf{A}) = \int_{\Omega_T} \| \mathbf{B}(\mathbf{A}) - \mathbf{B}_0 \|^2 d\Omega \\ \text{uc} & \\ & -\nabla \times [\nu_0 \nu_r(\mathbf{p}) \mathbf{B}(\mathbf{A}) - \mathbf{M}(\mathbf{p})] + \mathbf{J}(\mathbf{p}) = 0 \text{ in } \Omega, \\ \text{with} & \mathbf{A} \times \mathbf{n} = \mathbf{0} \text{ in } \partial\Omega \\ & \mathbf{p} \in \mathcal{P} := \{\mathbf{p} \in (L^\infty(\Omega))^7 : p_i^{\min} \le p_i \le p_i^{\max}, i = 1, \dots, 7\} \end{array} \right. \quad (4) $$

Where $\mathbf{p}^{\min}$ and $\mathbf{p}^{\max}$ are the lower and upper bounds of the variable **p**. Note that in $(\mathcal{P})$, the variables **p** and **A** are separated and linked by the constraint which is the Maxwell equation.
we detail our algorithmic approach to solve topology optimization problems using SIMP techniques. Section 5 is dedicated to the gradient verification (by comparing it to the finite difference technique) and to numerical tests to design a 3D magnetic circuit. Section 6 concludes.

## 2 Problem description

In this work, $\Omega$ represents the whole studied domain. Within $\Omega$, $\Omega_V$ represents the variable domain in which the matter can be added or removed. This matter can be permanent magnet in $\Omega_V(\mathbf{M})$, coils with current in $\Omega_V(\mathbf{J})$ or ferromagnetic material in $\Omega_V(\nu_r)$. $\Omega_T$ represents the target domain in which a particular magnetic flux density $\mathbf{B}_0$ has to be imposed playing with the matter in $\Omega_V$. In this work $\Omega_V(\mathbf{M})$, $\Omega_V(\mathbf{J})$ and $\Omega_V(\nu_r)$ have to be disconnected.

The topology optimization problem is defined as follows:

* The objective function is the squared error in a target zone $\Omega_T$ between an imposed magnetic flux density $\mathbf{B}_0$ and the magnetic flux density $\mathbf{B}$ obtained with respect to the magnetic sources and the material distributions in all the domain $\Omega$. Note that $\mathbf{B}$ has to be computed by solving Maxwell magnetostatic equations.

* The design variables are magnetic sources (active material) such as ferromagnetic material, coil with current density and permanent magnets. The associated variables for all the variable domains $\Omega_V$ are the distribution of relative reluctivity $\nu_r$, current density $\mathbf{J}$ (in the 3 directions $(x, y, z)$) and magnetization $\mathbf{M}$ (in the 3 directions $(x, y, z)$). Hence, we define a global variable $\mathbf{p}$ as:

$$ \mathbf{p}: \left| \begin{array}{l} \Omega_V \quad \rightarrow \quad \mathbb{R}^7 \\ (x, y, z)^T \quad \longmapsto \quad \mathbf{p}(x, y, z) = (p_1, \dots, p_7)^T \end{array} \right. $$

Where $\mathbf{p} = (\nu_r, J_x, J_y, J_z, M_x, M_y, M_z)^T$ and we define the three following functions: $\nu_r(\mathbf{p}) = \nu_r, \mathbf{J}(\mathbf{p}) = (J_x, J_y, J_z)^T$ and $\mathbf{M}(\mathbf{p}) = (M_x, M_y, M_z)^T$.

Note that, $\mathbf{p}$ is defined in the whole domain $\Omega$, $\mathbf{p}$ is fixed on $\Omega \setminus \Omega_V$ and is variable on $\Omega_V$.

Therefore, the optimization problem that we addressed, is the following:

$$ \min_{\mathbf{p}} F_o(\mathbf{p}) = \int_{\Omega_T} \| \mathbf{B}(\mathbf{A}(\mathbf{p})) - \mathbf{B}_0 \|^2 d\Omega, \quad (1) $$

where $\mathbf{A}$ corresponds to the magnetic potential vector obtained by solving the Maxwell equations.

In this work, the constitutive law of the magnetic materials is supposed to be linear in order to simplify the computations in 3D; some extensions to the non-linear case can be done from (Gauthey et al. 2021) and (Youness & Messine 2019a).

Hence, considering the linear case, we have:

$$ \mathbf{H} = \nu \mathbf{B} - \mathbf{M} \quad (2) $$

Where $\nu = \nu_0\nu_r$ is the magnetic reluctivity of the material, $\nu_0$ is the air reluctivity and $\nu_r$ is the relative reluctivity of the material.

In different points of domain $\Omega$, we have for $\mathbf{p}$:

* For air: $\nu_r = 1$, $\mathbf{M} = \mathbf{0}$ and $\mathbf{J} = \mathbf{0}$ in the 3 directions.

* For ferromagnetic material: $\nu_r$ is a constant depending on the material (in this paper, $\nu_r = \frac{1}{1000}$), $\mathbf{M} = \mathbf{0}$ and $\mathbf{J} = \mathbf{0}$ in the 3 directions.

* For permanent magnets: $\nu_r = 1$, $\mathbf{M}$ is a constant bounded vector in 3D and $\mathbf{J} = \mathbf{0}$ in the 3 directions.

* For coils: $\nu_r = 1$, $\mathbf{M} = \mathbf{0}$ and $\mathbf{J}$ is a constant bounded vector in 3D.
Moreover if **p** satisfies the Maxwell constraint, then $F_o(\mathbf{p}) = F(\mathbf{p}, \mathbf{A})$ and therefore, for such a point **p**:

$$ \frac{\partial F_o}{\partial p_i}(\mathbf{p}) = \frac{\partial F}{\partial p_i}(\mathbf{p}, \mathbf{A}) $$

# 3 Sensibility analysis from 2D to 3D

Efficient optimization codes are based on computations of the derivatives. These derivatives can be obtained with evaluations of the objective function considering small variations of the variable **p**. This defines the well-known finite difference method. However, all these evaluations require to solve Maxwell equations using a numerical solver. Hence, depending directly to the number of variables, the finite difference technique is extremely expensive.

Hence, adjoint methods has been developed in order to compute the derivatives with more efficiency. It was first developed in mechanics [2, 5], then extended to 2D magnetostatic [14, 16] and to 2D electrical machines [9, 6].

Thus, in the first subsection 3.1, the 2D magnetosatic adjoint method is recalled. In the following subsection 3.2, the adjoint method is extended to the 3D-case.

## 3.1 Sensibility analysis in 2D

Generally, due to symmetries and invariances along one direction, 3D electromagnetism problems can be simplified in 2D problems.

Considering the Cartesian coordinates in 3D: $(x, y, z)$ and an invariance following $z$, the magnetic flux density can be reduced to $\mathbf{B} = (B_x, B_y)^T$. Furthermore:

* The potential vector is a scalar:

$$ \mathbf{A} = (0, 0, A(x, y)) \equiv A(x, y) $$

* The magnetic density flux is a 2D vector:

$$ \mathbf{B}(A) = \mathrm{curl}(\mathbf{A}) = (\partial_y A, -\partial_x A, 0) \equiv (\partial_y A, -\partial_x A) $$

* The boundary condition is reduced to $A = 0$ on the boundary domain $\partial\Omega$.

In order to solve problem ($\mathcal{P}$), it is necessary to compute the derivative of the function $F$ with respect to the design variable **p**. The sensitivity calculation in 2D is summarized below (for more details see [14] and [16]):

* For problem ($\mathcal{P}$), the Lagrangian function $\mathcal{L}$ is introduced as follows:

$$ \mathcal{L}(\mathbf{p}, A, \lambda) = F(\mathbf{p}, A) + \int_{\Omega} \lambda \mu_0 [-\nabla \times (\nu_0 \nu_r(\mathbf{p})\mathbf{B}(A) - \mathbf{M}(\mathbf{p})) + J(\mathbf{p})] d\Omega, $$

where the equality constraint is multiplied by the Lagrangian multiplier $\lambda$ which is a scalar function. $\lambda$ belongs to $H_0^1(\Omega)$ (a standard Sobolev space) with the same properties as $A$. Note that, the equality constraint of ($\mathcal{P}$) is also multiplied by $\mu_0$ (the void permeability = $\frac{1}{\nu_0}$) in order to ensure the homogeneity considering the physical units in the summed terms. Note that this Lagrangian formulation is in 2D, with $A$ a scalar function, $\mathbf{M}$ a vector (in 2D), $J$ is also a scalar (in direction $z$) and the curl operator in 2D is scalar $\nabla \times (v_x, v_y)^T = \partial_x v_y - \partial_y v_x$.

* Using the Green formula in 2D, we obtain:

$$ \frac{\partial \mathbf{B}}{\partial A}(A) = \mathbf{B}(.), $$
Because **p** is extended on the whole domain Ω, in Ω \ Ω<sub>V</sub> (with Ω<sub>V</sub> = Ω<sub>V</sub>(ν<sub>r</sub>) ∪ Ω<sub>V</sub>(**M**) ∪ Ω<sub>V</sub>(**J**)), **p** is fixed and cannot change during the optimization process. Thus, in Ω \ Ω<sub>V</sub>, **p** can be air, a ferromagnetic material, a permanent magnet or a coil. In the variable domain Ω<sub>V</sub>, **p** is variable and can be air or a ferromagnetic material if (x, y, z)<sup>T</sup> belongs to Ω<sub>V</sub>(ν<sub>r</sub>), air or a permanent magnet if (x, y, z)<sup>T</sup> belongs to Ω<sub>V</sub>(**M**) and air or a coil if (x, y, z)<sup>T</sup> belongs to Ω<sub>V</sub>(**J**).

The Maxwell-Ampere equation ∇ × **H** = **J** becomes, using the constitutive law of the material described in equation (2) and the existence of the potential vector **A** such that **B** = ∇ × **A**:

$$ -\nabla \times [\nu_0 \nu_r (\mathbf{p}) \mathbf{B}(\mathbf{A}(\mathbf{p})) - \mathbf{M}(\mathbf{p})] + \mathbf{J}(\mathbf{p}) = 0, \text{ in } \Omega. \quad (3) $$

Equation (3) is called, in this paper, the primal problem. Note that **B**(.) is the curl operator ∇ × (.) which is also denoted by **B**(**A**) = **curl**(**A**). The Dirichlet condition **A** × **n** = **0** is imposed on the domain boundary ∂Ω, where **n** is the normal vector (this hypothesis is relevant because the border of the domain Ω is far from the zone of interest Ω<sub>T</sub>).

The partial derivative ∂F<sub>o</sub>/∂p<sub>i</sub> is calculated using the adjoint method. Hence, from equations (1) and (3), the optimization problem (1) can be equivalently formulated as:

$$ (\mathcal{P}) \quad \left\{ \begin{array}{ll} \min_{\mathbf{p},\mathbf{A}} & F(\mathbf{p}, \mathbf{A}) = \int_{\Omega_T} || \mathbf{B}(\mathbf{A}) - \mathbf{B}_0 ||^2 d\Omega \\ \text{uc} & -\nabla \times [\nu_0 \nu_r(\mathbf{p})\mathbf{B}(\mathbf{A}) - \mathbf{M}(\mathbf{p})] + \mathbf{J}(\mathbf{p}) = 0 \text{ in } \Omega, \\ & \text{with } \mathbf{A} \times \mathbf{n} = \mathbf{0} \text{ in } \partial\Omega \\ & \mathbf{p} \in \mathcal{P} := \{\mathbf{p} \in (L^\infty(\Omega))^7 : p_i^{\min} \le p_i \le p_i^{\max}, i = 1, \dots, 7\} \end{array} \right. \quad (4) $$

Where **p**<sup>min</sup> and **p**<sup>max</sup> are the lower and upper bounds of the variable **p**. Note that in (P), the variables **p** and **A** are separated and linked by the constraint which is the Maxwell equation. Moreover if **p** satisfies the Maxwell constraint, then F<sub>o</sub>(**p**) = F(**p**, **A**) and therefore, for such a point **p**:

$$ \frac{\partial F_o}{\partial p_i}(\mathbf{p}) = \frac{\partial F}{\partial p_i}(\mathbf{p}, \mathbf{A}). $$

# 3 Sensitivity analysis from 2D to 3D

Efficient optimization codes are based on computations of the derivatives. These derivatives can be obtained with evaluations of the objective function considering small variations of the variable **p**. This defines the well-known finite difference method. However, all these evaluations require to solve Maxwell equations using a numerical solver. Hence, because it depends directly to the number of variables, the finite difference technique is extremely expensive.

Thus, adjoint methods has been developed in order to compute the derivatives with more efficiency. It was first developed in mechanics (Allaire 2007, Bendsøe & Sigmund 2003), then extended to 2D magnetostatics (Sanogo 2016, Sanogo et al. 2014) and to 2D electrical machines (Gauthey et al. 2021, Cherrière et al. 2022).

Thus, in the first subsection 3.1, the 2D magnetostatic continuous adjoint method is recalled. In the following subsection 3.2, the continuous adjoint method is extended to the 3D-case.

## 3.1 Sensitivity analysis in 2D

Generally, due to symmetries and invariances along one direction, 3D magnetostatic problems can be simplified in 2D problems.

Considering the Cartesian coordinates in 3D: (x, y, z) and an invariance following z, the magnetic flux density can be reduced to **B** = (B<sub>x</sub>, B<sub>y</sub>)<sup>T</sup>. Furthermore:

* The potential vector is a scalar:

$$ \mathbf{A} = (0, 0, A(x, y)) \equiv A(x, y) $$
and the adjoint problem in 2D is:

$$ \nabla \times [\nu \mathbf{B}(\lambda) - \mathbf{M}_a \mathbf{1}_{\Omega_T}] = 0, \text{ in } \Omega, \quad (5) $$

where $\lambda$ is named the adjoint variable and $\mathbf{M}_a = \frac{2}{\mu_0}(B_x - B_{0x}, B_x - B_{0x})^T$ can be understood as a permanent magnet in $\Omega_T$. The adjoint problem is therefore equivalent to a magnetostatic problem with only one source $\mathbf{M}_a$.

Note that, the computation of $\mathbf{M}_a$ depends on the value of the field $\mathbf{B} = (B_x, B_y)^T$ in $\Omega_T$ which must be previously calculated by solving the primal problem defined by equation (3).

* From the the equality $\frac{\partial F_o(\mathbf{p})}{\partial \mathbf{p}_i} = \frac{\partial \mathcal{L}}{\partial \mathbf{p}_i}(\mathbf{p}, A(\mathbf{p}), \lambda)$ which is verified for any couple $(\mathbf{p}, A)$ satisfying the constraint (3), the sensitivity formula is obtained by deriving the Lagrangian $\mathcal{L}$ with respect to the design variable $\mathbf{p}$:

$$ \frac{\partial F_o(\mathbf{p})}{\partial \mathbf{p}_i} = \int_{\Omega} \mu_0 \left[ \lambda \frac{\partial J(\mathbf{p})}{\partial \mathbf{p}_i} + \mathbf{B}(\lambda)^T \frac{\partial \mathbf{M}(\mathbf{p})}{\partial \mathbf{p}_i} - \nu_0 \frac{\partial \nu_r(\mathbf{p})}{\partial \mathbf{p}_i} \mathbf{B}(\lambda)^T \mathbf{B}(A) \right] d\Omega, \quad (6) $$

for each component $\mathbf{p}_i$ of $\mathbf{p}$; only $p_1 = \nu_r, p_4 = J_z, p_5 = M_x$ and $p_6 = M_y$ are used in this 2D problem.

Note that in the sensitivity formula (6), $A$ and $\lambda$ are solutions of equations (3) and (5), respectively.

Equations (3) and (5) can be solved by the finite element method whose weak formulation is written as:

$$ \int_{\Omega} \nu \mathbf{B}(X) \cdot \mathbf{B}(\phi) d\Omega - \int_{\Omega} \mathbf{M} \cdot \mathbf{B}(\phi) d\Omega - \int_{\Omega} J\phi d\Omega = 0, \forall \phi \in H_0^1(\Omega). \quad (7) $$

For example, in order to solve the primal problem of designing a simple magnetic circuit with only a coil, a ferromagnetic material and no permanent magnet, we just have to take $X = A$, $\mathbf{M} = \mathbf{0}$ and the only source is the current $J$. Considering the same example, the adjoint problem is derived from equation (7) by taking $X = \lambda, J = 0$ and the only source is $\mathbf{M} = \mathbf{M}_a$.

In a 2D problem, equation (7) can be simplified, such as:

$$ \int_{\Omega} \nu \nabla X \cdot \nabla \phi \, d\Omega - \int_{\Omega} \left( \frac{-M_y}{M_x} \right) \cdot \nabla \phi \, d\Omega - \int_{\Omega} J\phi \, d\Omega = 0, \forall \phi \in H_0^1(\Omega). \quad (8) $$

## 3.2 Sensitivity analysis and adjoint method in 3D

In this subsection, we detail the steps to compute the sensitivity of the objective function $F_o$ with respect to the design variable $\mathbf{p}_i$ for a 3D problem. Note that in 3D, the potential vector and the magnetic flux density have non zero-components, and then,

* The potential vector is:

$$ \mathbf{A} = (A_x(x, y, z), A_y(x, y, z), A_z(x, y, z)) $$

* The magnetic flux density is:

$$ \mathbf{B}(\mathbf{A}) = \mathrm{curl}(\mathbf{A}) = (\partial_y A_z - \partial_z A_y, \partial_z A_x - \partial_x A_z, \partial_x A_y - \partial_y A_x)^T $$

* The boundary condition is $\mathbf{A} \times \mathbf{n} = \mathbf{0}$ on the boundary domain $\partial\Omega$.

The following theorem is dedicated to the computations in 3D of the partial derivatives of the original function $F_o$.
* The magnetic density flux is a 2D vector:

$$ \mathbf{B}(A) = \mathbf{curl}(\mathbf{A}) = (\partial_y A, -\partial_x A, 0) \equiv (\partial_y A, -\partial_x A) $$

* The boundary condition is reduced to $A = 0$ on the boundary domain $\partial\Omega$.

In order to solve problem (P), it is necessary to compute the derivative of the function $F_o$ with respect to the design variable $\mathbf{p}$. The sensitivity calculation in 2D is summarized below (for more details see (Sanogo 2016) and (Sanogo et al. 2014)):

* For problem (P), the Lagrangian function $\mathcal{L}$ is introduced as follows:

$$ \mathcal{L}(\mathbf{p}, A, \lambda) = F(\mathbf{p}, A) + \int_{\Omega} \lambda \mu_0 [-\nabla \times (\nu_0 \nu_r(\mathbf{p})\mathbf{B}(A) - \mathbf{M}(\mathbf{p})) + J(\mathbf{p})] d\Omega, $$

where the equality constraint is multiplied by the Lagrangian multiplier $\lambda$ which is a scalar function. $\lambda$ belongs to $H_0^1(\Omega)$ (a standard Sobolev space) with the same properties as $A$. Note that, the equality constraint of (P) is also multiplied by $\mu_0$ (the void permeability = $\frac{1}{\nu_0}$) in order to ensure the homogeneity considering the physical units in the summed terms. Note that this Lagrangian formulation is in 2D, with $A$ a scalar function, $\mathbf{M}$ a vector (in 2D), $J$ is also a scalar (in direction $z$) and the curl operator in 2D is scalar $\nabla \times (v_x, v_y)^T = \partial_x v_y - \partial_y v_x$.

* Using the Green formula in 2D, we obtain:

$$ \frac{\partial \mathbf{B}}{\partial A}(A) = \mathbf{B}(\cdot), $$

and the continuous adjoint problem in 2D is:

$$ \nabla \times [\nu \mathbf{B}(\lambda) - \mathbf{M}_a \mathbf{1}_{\Omega_T}] = 0, \text{ in } \Omega, \quad (5) $$

with $\mathbf{1}_{\Omega_T}$ is the indicator function, defined by:

$$ \mathbf{1}_{\Omega_T}(x) = \begin{cases} 1 & \text{if } x \in \Omega_T \\ 0 & \text{if } x \in \Omega \setminus \Omega_T \end{cases} $$

where $\lambda$ is named the adjoint variable and $\mathbf{M}_a = \frac{2}{\mu_0}(B_x - B_{0x}, B_y - B_{0y})^T$ can be understood as a permanent magnet in $\Omega_T$. This adjoint problem is therefore equivalent to a magnetostatic problem with only one source $\mathbf{M}_a$.

Note that, the computation of $\mathbf{M}_a$ depends on the value of the field $\mathbf{B} = (B_x, B_y)^T$ in $\Omega_T$ which must be previously calculated by solving the primal problem defined by equation (3).

* From the equality $\frac{\partial F_o(\mathbf{p})}{\partial \mathbf{p}_i} = \frac{\partial \mathcal{L}}{\partial \mathbf{p}_i}(\mathbf{p}, A(\mathbf{p}), \lambda)$ which is verified for any couple $(\mathbf{p}, A)$ satisfying the constraint (3), the sensitivity formula is obtained by deriving the Lagrangian $\mathcal{L}$ with respect to the design variable $\mathbf{p}$:

$$ \frac{\partial F_o(\mathbf{p})}{\partial \mathbf{p}_i} = \int_{\Omega} \mu_0 [\lambda \frac{\partial J(\mathbf{p})}{\partial \mathbf{p}_i} + \mathbf{B}(\lambda) \cdot \frac{\partial \mathbf{M}(\mathbf{p})}{\partial \mathbf{p}_i} - \nu_0 \frac{\partial \nu_r(\mathbf{p})}{\partial \mathbf{p}_i} \mathbf{B}(\lambda) \cdot \mathbf{B}(A)] d\Omega, \quad (6) $$

for each component $\mathbf{p}_i$ of $\mathbf{p}$; only $p_1 = \nu_r, p_4 = J_z, p_5 = M_x$ and $p_6 = M_y$ are used in this 2D problem.

Note that in the sensitivity formula (6), $A$ and $\lambda$ are solutions of equations (3) and (5), respectively.

Equations (3) and (5) can be solved by the finite element method whose weak formulation is written as:

$$ \int_{\Omega} \nu \mathbf{B}(X) \cdot \mathbf{B}(\phi) d\Omega - \int_{\Omega} \mathbf{M} \cdot \mathbf{B}(\phi) d\Omega - \int_{\Omega} J \phi d\Omega = 0, \quad (7) $$
**Theorem 1.** The expression in 3D of the sensitivity of $F_o$ with respect to each component $p_i$ of $\mathbf{p}$ is:

$$ \frac{\partial F_o(\mathbf{p})}{\partial p_i} = \int_{\Omega} \mu_0 \left( \lambda^T \frac{\partial \mathbf{J}(\mathbf{p})}{\partial p_i} + \mathbf{B}(\lambda)^T \frac{\partial \mathbf{M}(\mathbf{p})}{\partial p_i} - \nu_0 \frac{\partial \nu_r(\mathbf{p})}{\partial p_i} \mathbf{B}(\lambda)^T \mathbf{B}(\mathbf{A}) \right) d\Omega, \quad (9) $$

where $\mathbf{B}(\lambda)$ is obtained by solving the following adjoint problem:

$$ \nabla \times [\mathbf{M}_a \mathbf{1}_{\Omega_T} - \nu \mathbf{B}(\lambda)] = 0, \quad (10) $$

with $\mathbf{M}_a$ which is the adjoint magnetization defined as follows:

$$ \mathbf{M}_a = \frac{2}{\mu_0} \times \begin{pmatrix} B_x - B_{0x} \\ B_y - B_{0y} \\ B_z - B_{0z} \end{pmatrix}, $$

with $\mathbf{B} = (B_x, B_y, B_z)^T$ which is the magnetic flux density.

Note that equation (10) is the adjoint problem in 3D. It is equivalent to a magnetostatic problem where $\mathbf{M}_a$ is the only source. This problem has to be solved numerically in order to provide $\lambda$ and $\mathbf{B}(\lambda)$.

*Proof.* The proof of Theorem 1 is based on the following 3 steps:

1. Construction of the Lagrangian function by considering the formulation of problem (P).

2. Determination of the adjoint problem from the Karush-Khun-Tucker theorem.

3. Computation of the topological sensitivity.

Note that, in 3D, we work in the following Sobolev spaces

$$ \mathbf{H}(\mathrm{curl}, \Omega) = \{ \mathbf{w} \in \mathbf{L}^2(\Omega), \mathrm{curl}(\mathbf{w}) \in \mathbf{L}^2(\Omega) \}, $$

$$ \mathbf{H}_0(\mathrm{curl}, \Omega) = \{ \mathbf{w} \in \mathbf{H}(\mathrm{curl}, \Omega), \text{ with } \mathbf{w} \times \mathbf{n} = \mathbf{0} \text{ in } \partial\Omega \}, $$

because we will need to integrate the vectors $\mathbf{A}$, $\lambda$, $\mathrm{curl}(\mathbf{A})$ and $\mathrm{curl}(\lambda)$ in $\Omega$. $\mathbf{L}^\infty(\Omega)$ is the space of bounded vector functions on $\Omega$.

1. Construction of the Lagrangian function:

The Lagrangian $\mathcal{L}$ is defined as the addition of the objective function and the sum of the equal constraints multiplied by a Lagrange multiplier $\lambda$:

$$ \mathcal{L}(\mathbf{p}, \mathbf{A}, \lambda) = F(\mathbf{p}, \mathbf{A}) + \int_{\Omega} \lambda^T \mu_0 [-\nabla \times (\nu_0 \nu_r(\mathbf{p})\mathbf{B}(\mathbf{A}) - \mathbf{M}(\mathbf{p})) + \mathbf{J}(\mathbf{p})] d\Omega, $$

where $(\mathbf{p}, \mathbf{A}, \lambda) \in \mathbf{L}^\infty(\Omega) \times \mathbf{H}_0(\mathrm{curl}, \Omega) \times \mathbf{H}_0(\mathrm{curl}, \Omega)$ and $\mu_0 = \frac{1}{\nu_0}$ the void permeability. The equality constraint of $(\mathcal{P})$ is multiplied by $\mu_0$ (the void permeability $= \frac{1}{\nu_0}$) in order to ensure homogeneity considering the physical units in the summed terms.

Using Green's formula in 3D, we obtain:

$$ \int_{\Omega} \lambda^T \mu_0 [-\nabla \times (\nu \mathbf{B}(\mathbf{A}) - \mathbf{M}) + \mathbf{J}] d\Omega = \int_{\Omega} \mu_0 \lambda^T \mathbf{J} d\Omega \\ - \int_{\Omega} \mu_0 \mathrm{curl}(\lambda)^T (\nu \mathbf{B}(\mathbf{A}) - \mathbf{M}) d\Omega \\ + \int_{\partial\Omega} \mu_0 (\lambda \times \mathbf{n})^T (\nu \mathbf{B}(\mathbf{A}) - \mathbf{M}) d\partial\Omega. $$
where $\phi$ is a test function of $H_0^1(\Omega)$.

For example, in order to solve the primal problem of designing a simple magnetic circuit with only a coil, a ferromagnetic material and no permanent magnet, we just have to take $X = A$, $M = 0$ and the only source is the current $J$. Considering the same example, the adjoint problem is derived from equation (7) by taking $X = \lambda$, $J = 0$ and the only source is $M = M_a$.

In a 2D problem, equation (7) can be simplified, such as:

$$ \int_{\Omega} \nu \nabla X \cdot \nabla \phi \, d\Omega - \int_{\Omega} \left( \frac{-M_y}{M_x} \right) \cdot \nabla \phi \, d\Omega - \int_{\Omega} J \phi \, d\Omega = 0, \forall \phi \in H_0^1(\Omega). \quad (8) $$

## 3.2 Sensitivity analysis and continuous adjoint method in 3D

In this subsection, we detail the steps to compute the sensitivity of the objective function $F_o$ with respect to the design variable $p_i$ for a 3D problem. Note that in 3D, the potential vector and the magnetic flux density have non zero-components, and then,

* The potential vector is:

$$ \mathbf{A} = (A_x(x, y, z), A_y(x, y, z), A_z(x, y, z)) $$

* The magnetic flux density is:

$$ \mathbf{B}(\mathbf{A}) = \mathrm{curl}(\mathbf{A}) = (\partial_y A_z - \partial_z A_y, \partial_z A_x - \partial_x A_z, \partial_x A_y - \partial_y A_x)^T $$

* The boundary condition is $\mathbf{A} \times \mathbf{n} = 0$ on the boundary domain $\partial\Omega$.

The following theorem is dedicated to the computations in 3D of the partial derivatives of the original function $F_o$.

**Theorem 1.** The expression in 3D of the sensitivity of $F_o$ with respect to each component $p_i$ of $\mathbf{p}$ is:

$$ \frac{\partial F_o(\mathbf{p})}{\partial p_i} = \int_{\Omega} \mu_0 \left( \lambda \cdot \frac{\partial \mathbf{J}(\mathbf{p})}{\partial p_i} + \mathbf{B}(\lambda) \cdot \frac{\partial \mathbf{M}(\mathbf{p})}{\partial p_i} - \nu_0 \frac{\partial \nu_r(\mathbf{p})}{\partial p_i} \mathbf{B}(\lambda) \cdot \mathbf{B}(\mathbf{A}) \right) d\Omega, \quad (9) $$

where $\mathbf{B}(\lambda)$ is obtained by solving the following continuous adjoint problem:

$$ \nabla \times [\mathbf{M}_a \mathbf{1}_{\Omega_T} - \nu \mathbf{B}(\lambda)] = 0, \quad (10) $$

with $\mathbf{M}_a$ which is the adjoint magnetization defined as follows:

$$ \mathbf{M}_a = \frac{2}{\mu_0} \begin{pmatrix} B_x - B_{0x} \\ B_y - B_{0y} \\ B_z - B_{0z} \end{pmatrix}, $$

with $\mathbf{B} = (B_x, B_y, B_z)^T$ which is the magnetic flux density.

Note that equation (10) is the adjoint problem in 3D. It is equivalent to a magnetostatic problem where $\mathbf{M}_a$ is the only source. This problem has to be solved numerically in order to provide $\lambda$ and $\mathbf{B}(\lambda)$.

*Proof.* The proof of Theorem 1 is based on the following 3 steps:

1. Construction of the Lagrangian function by considering the formulation of problem (P).

2. Determination of the continuous adjoint problem from the Karush-Khun-Tucker theorem.

3. Computation of the topological sensitivity.
Because $\lambda \in H_0(\text{curl}, \Omega)$, $\lambda \times \mathbf{n} = \mathbf{0}$ in $\partial\Omega$, hence we have:

$$ 
\begin{aligned}
\int_{\Omega} \lambda^T \mu_0 [-\nabla \times (\nu \mathbf{B}(\mathbf{A}) - \mathbf{M}) + \mathbf{J}] d\Omega &= \int_{\Omega} \mu_0 \lambda^T \mathbf{J} d\Omega - \int_{\Omega} \mu_0 \mathbf{B}(\lambda)^T (\nu \mathbf{B}(\mathbf{A}) - \mathbf{M}) d\Omega \\
&= \int_{\Omega} \mu_0 [\lambda^T \mathbf{J} - \nu \mathbf{B}(\lambda)^T \mathbf{B}(\mathbf{A}) + \mathbf{B}(\lambda)^T \mathbf{M}] d\Omega.
\end{aligned}
$$

Finally, the Lagrangian function is:

$$ 
\begin{aligned}
\mathcal{L}(\mathbf{p}, \mathbf{A}, \lambda) = & \int_{\Omega_T} \| \mathbf{B}(\mathbf{A}) - \mathbf{B}_0 \|^2 d\Omega \\
& + \int_{\Omega} \mu_0 [\lambda^T \mathbf{J}(\mathbf{p}) - \nu(\mathbf{p})\mathbf{B}(\lambda)^T \mathbf{B}(\mathbf{A}) + \mathbf{B}(\lambda)^T \mathbf{M}(\mathbf{p})] d\Omega.
\end{aligned}
\quad (11) $$

## 2. Determination of the adjoint problem:

The directional derivative of a function $f: E \to F$ in $u$ in the direction $h$ with $(u, h) \in E^2$ is defined by:

$$ D_h f(u) = \lim_{t \to 0} \frac{f(u + th) - f(u)}{t}. $$

If $E$ is an euclidean vector space, and $f$ is a real-valued differentiable application, then we can use the gradient of $f$ to express the directional derivative:

$$ D_h f(u) = \langle \nabla f(u) | h \rangle. $$

In order to determine the adjoint system, we use the necessary optimality conditions at the first order of the Karush-Kuhn-Tucker (KKT) theorem. Thus, the directional derivative of $\mathcal{L}$ in $\mathbf{A}$ in any direction $\varphi$ must satisfy:

$$ \langle \frac{\partial \mathcal{L}(\mathbf{p}, \mathbf{A}, \lambda)}{\partial \mathbf{A}} | \varphi \rangle = 0, \forall \varphi \in H_0(\text{rot}, \Omega). \quad (12) $$

By considering a function $f$ defined by:

$$ \mathcal{L}(\mathbf{p}, \mathbf{A}, \lambda) = \int_{\Omega} f(\mathbf{p}, \mathbf{A}, \lambda) d\Omega, \quad (13) $$

where,

$$ f(\mathbf{p}, \mathbf{A}, \lambda) = \| \mathbf{B}(\mathbf{A}) - \mathbf{B}_0 \|^2 \mathbb{1}_{\Omega_T} + \mu_0 [\lambda^T \mathbf{J}(\mathbf{p}) - \nu(\mathbf{p})\mathbf{B}(\lambda)^T \mathbf{B}(\mathbf{A}) + \mathbf{B}(\lambda)^T \mathbf{M}(\mathbf{p})]. $$

By developing the expression of the Lagrangian (13), and using the linearity property of integrals, we obtain:

$$ 
\begin{aligned}
\langle \frac{\partial \mathcal{L}(\mathbf{p}, \mathbf{A}, \lambda)}{\partial \mathbf{A}} | \varphi \rangle &= \lim_{t \to 0} \frac{\mathcal{L}(\mathbf{p}, \mathbf{A} + t\varphi, \lambda) - \mathcal{L}(\mathbf{p}, \mathbf{A}, \lambda)}{t} \\
&= \lim_{t \to 0} \frac{\int_{\Omega} f(\mathbf{p}, \mathbf{A} + t\varphi, \lambda) d\Omega - \int_{\Omega} f(\mathbf{p}, \mathbf{A}, \lambda) d\Omega}{t} \\
&= \int_{\Omega} \lim_{t \to 0} \frac{f(\mathbf{p}, \mathbf{A} + t\varphi, \lambda) - f(\mathbf{p}, \mathbf{A}, \lambda)}{t} d\Omega \\
&= \int_{\Omega} D_{\varphi} f(\mathbf{A}) d\Omega \\
&= \int_{\Omega} \langle \nabla_{\mathbf{A}} f | \varphi \rangle d\Omega \\
&= \int_{\Omega} \frac{\partial f}{\partial A_x} \cdot \varphi_x + \frac{\partial f}{\partial A_y} \cdot \varphi_y + \frac{\partial f}{\partial A_z} \cdot \varphi_z d\Omega.
\end{aligned}
$$
Note that, in 3D, we work in the following Sobolev spaces:

$$ \mathbf{H}(\mathbf{curl}, \Omega) = \{\mathbf{w} \in \mathbf{L}^2(\Omega), \mathbf{curl}(\mathbf{w}) \in \mathbf{L}^2(\Omega)\}, $$

$$ \mathbf{H}_0(\mathbf{curl}, \Omega) = \{\mathbf{w} \in \mathbf{H}(\mathbf{curl}, \Omega), \text{ with } \mathbf{w} \times \mathbf{n} = \mathbf{0} \text{ in } \partial\Omega\}, $$

because we will need to integrate the vectors **A**, **λ**, **curl(A)** and **curl(λ)** in Ω. L∞(Ω) is the space of bounded vector functions on Ω.

### 1. Construction of the Lagrangian function:

The Lagrangian function L is defined as the addition of the objective function and the sum of the equal constraints multiplied by a Lagrange multiplier λ:

$$ \mathcal{L}(\mathbf{p}, \mathbf{A}, \mathbf{\lambda}) = F(\mathbf{p}, \mathbf{A}) + \int_{\Omega} \mathbf{\lambda} \cdot \mu_0 [-\nabla \times (\nu_0 \nu_r (\mathbf{p}) \mathbf{B}(\mathbf{A}) - \mathbf{M}(\mathbf{p})) + \mathbf{J}(\mathbf{p})] d\Omega, $$

where (p, A, λ) ∈ L∞(Ω) × H₀(curl, Ω) × H₀(curl, Ω) and μ₀ = 1/ν₀ the void permeability. The equality constraint of (P) is multiplied by μ₀ (the void permeability = 1/ν₀) in order to ensure homogeneity considering the physical units in the summed terms. By using Green formula in 3D:

$$ - \int_{\Omega} \mathbf{v} \cdot (\nabla \times \mathbf{u}) d\Omega = - \int_{\Omega} (\nabla \times \mathbf{v}) \cdot \mathbf{u} d\Omega + \int_{\partial\Omega} \mathbf{v} \times \mathbf{n} \cdot \mathbf{u} d\partial\Omega $$

and by taking: v = μ₀λ and u = νB(A) - M, we obtain:

$$ \int_{\Omega} \mathbf{\lambda} \cdot \mu_0 [-\nabla \times (\nu \mathbf{B}(\mathbf{A}) - \mathbf{M}) + \mathbf{J}] d\Omega = \int_{\Omega} \mu_0 \mathbf{\lambda} \cdot \mathbf{J} d\Omega - \int_{\Omega} \mu_0 \mathbf{curl}(\mathbf{\lambda}) \cdot (\nu \mathbf{B}(\mathbf{A}) - \mathbf{M}) d\Omega + \int_{\partial\Omega} \mu_0 (\mathbf{\lambda} \times \mathbf{n}) \cdot (\nu \mathbf{B}(\mathbf{A}) - \mathbf{M}) d\partial\Omega. $$

Because λ ∈ H₀(curl, Ω), λ × n = 0 in ∂Ω, hence we have:

$$ \int_{\Omega} \mathbf{\lambda} \cdot \mu_0 [-\nabla \times (\nu \mathbf{B}(\mathbf{A}) - \mathbf{M}) + \mathbf{J}] d\Omega = \int_{\Omega} \mu_0 \mathbf{\lambda} \cdot \mathbf{J} d\Omega - \int_{\Omega} \mu_0 \mathbf{B}(\mathbf{\lambda}) \cdot (\nu \mathbf{B}(\mathbf{A}) - \mathbf{M}) d\Omega \\ = \int_{\Omega} \mu_0 [\mathbf{\lambda} \cdot \mathbf{J} - \nu \mathbf{B}(\mathbf{\lambda}) \cdot \mathbf{B}(\mathbf{A}) + \mathbf{B}(\mathbf{\lambda}) \cdot \mathbf{M}(\mathbf{p})] d\Omega. $$

Finally, the Lagrangian function is:

$$ \mathcal{L}(\mathbf{p}, \mathbf{A}, \mathbf{\lambda}) = \int_{\Omega_T} \| \mathbf{B}(\mathbf{A}) - \mathbf{B}_0 \|^2 d\Omega + \int_{\Omega} \mu_0 [\mathbf{\lambda} \cdot \mathbf{J}(\mathbf{p}) - \nu(\mathbf{p})\mathbf{B}(\mathbf{\lambda}) \cdot \mathbf{B}(\mathbf{A}) + \mathbf{B}(\mathbf{\lambda}) \cdot \mathbf{M}(\mathbf{p})] d\Omega. \quad (11) $$

### 2. Determination of the continuous adjoint problem:

The directional derivative of a function f : E → F in u in the direction h with (u, h) ∈ E² is defined by:

$$ D_h f(u) = \lim_{t \to 0} \frac{f(u + th) - f(u)}{t}. $$

If E is an euclidean vector space, and f is a real-valued differentiable application, then we can use the gradient of f to express the directional derivative:

$$ D_h f(u) = \langle \nabla f(u) | h \rangle. $$
Using the derivation of compound functions, for each of the partial derivatives, where $A$, denotes $A_x$, $A_y$ or $A_z$, then we have:

$$
\begin{align*}
\frac{\partial f}{\partial A_{\cdot}} &= \frac{\partial}{\partial A_{\cdot}} \left[ \| \mathbf{B}(\mathbf{A}) - \mathbf{B}_0 \|^2 \mathbf{1}_{\Omega_T} + \mu_0 \left[ \boldsymbol{\lambda}^T \cdot \mathbf{J} - \nu \mathbf{B}(\boldsymbol{\lambda})^T \cdot \mathbf{B}(\mathbf{A}) + \mathbf{B}(\boldsymbol{\lambda})^T \cdot \mathbf{M} \right] \right] \\
&= \frac{\partial}{\partial \mathbf{B}(\mathbf{A})^T} (\| \mathbf{B}(\mathbf{A}) - \mathbf{B}_0 \|^2 \mathbf{1}_{\Omega_T}) \cdot \frac{\partial \mathbf{B}(\mathbf{A})}{\partial A_{\cdot}} - \mu_0 \nu \mathbf{B}(\boldsymbol{\lambda})^T \cdot \frac{\partial \mathbf{B}(\mathbf{A})}{\partial A_{\cdot}} \\
&= \mathbf{B}_a^T \cdot \frac{\partial \mathbf{B}(\mathbf{A})}{\partial A_{\cdot}} \mathbf{1}_{\Omega_T} - \mu_0 \nu \mathbf{B}(\boldsymbol{\lambda})^T \cdot \frac{\partial \mathbf{B}(\mathbf{A})}{\partial A_{\cdot}} \\
&= \left[ \mathbf{B}_a \mathbf{1}_{\Omega_T} - \mu_0 \nu \mathbf{B}(\boldsymbol{\lambda}) \right]^T \cdot \frac{\partial \mathbf{B}(\mathbf{A})}{\partial A_{\cdot}},
\end{align*}
$$

with:

$$
\mathbf{B}_a = 2 \times \begin{pmatrix} B_x - B_{0x} \\ B_y - B_{0y} \\ B_z - B_{0z} \end{pmatrix}.
$$

Thus, we obtain:

$$
< \frac{\partial \mathcal{L}(\mathbf{p}, \mathbf{A}, \boldsymbol{\lambda})}{\partial \mathbf{A}} | \boldsymbol{\varphi} > = \int_{\Omega} [\mathbf{B}_a \mathbf{1}_{\Omega_T} - \mu_0 \nu \mathbf{B}(\boldsymbol{\lambda})]^T \cdot \left[ \varphi_x \frac{\partial \mathbf{B}(\mathbf{A})}{\partial A_x} + \varphi_y \frac{\partial \mathbf{B}(\mathbf{A})}{\partial A_y} + \varphi_z \frac{\partial \mathbf{B}(\mathbf{A})}{\partial A_z} \right] d\Omega. \quad (14)
$$

Expanding the term with partial derivatives according to the components of **A**, we obtain:

$$
\begin{align*}
\left[ \varphi_x \frac{\partial \mathbf{B}(\mathbf{A})}{\partial A_x} + \varphi_y \frac{\partial \mathbf{B}(\mathbf{A})}{\partial A_y} + \varphi_z \frac{\partial \mathbf{B}(\mathbf{A})}{\partial A_z} \right] &= J_{\mathbf{B}}(\mathbf{A}) \cdot \boldsymbol{\varphi} \\
&= D_{\boldsymbol{\varphi}} \mathbf{B}(\mathbf{A}).
\end{align*}
$$

Where $J_{\mathbf{B}}(\mathbf{A})$ is the Jacobian matrix of the function $\mathbf{B}$ with respect to $\mathbf{A}$. For each element $\mathbf{A} = (A_x, A_y, A_z)^T \in \mathbb{R}^3$ it is associated $\mathbf{B}(\mathbf{A}) = (B_x(\mathbf{A}), B_y(\mathbf{A}), B_z(\mathbf{A}))^T \in \mathbb{R}^3$.

Using the definition of the directional derivative and the linearity of the **curl** operator, we have:

$$
D_{\varphi} \mathbf{B}(\mathbf{A}) = \mathbf{B}(\varphi). \tag{15}
$$

From equation (14) and using equation (15), we have:

$$
\int_{\Omega} [\mathbf{B}_a \mathbf{1}_{\Omega_T} - \mu_0 \nu \mathbf{B}(\boldsymbol{\lambda})]^T \cdot \mathbf{B}(\boldsymbol{\varphi}) d\Omega = 0.
$$

By applying the Green formula in 3D, we obtain:

$$
\int_{\Omega} \varphi^T \nabla \times [\mathbf{B}_a \mathbf{1}_{\Omega_T} - \mu_0 \nu \mathbf{B}(\boldsymbol{\lambda})] d\Omega + \int_{\partial \Omega} (\boldsymbol{\varphi} \times \mathbf{n})^T [\mathbf{B}_a \mathbf{1}_{\Omega_T} - \mu_0 \nu \mathbf{B}(\boldsymbol{\lambda})] d\partial \Omega = 0. \quad (16)
$$

Because $\boldsymbol{\varphi} \in \mathbf{H}_0(\mathrm{curl}, \Omega)$, $\boldsymbol{\varphi} \times \mathbf{n} = \mathbf{0}$ on $\partial\Omega$, we have:

$$
\int_{\Omega} \varphi^T \nabla \times [\mathbf{B}_a \mathbf{1}_{\Omega_T} - \mu_0 \nu \mathbf{B}(\boldsymbol{\lambda})] d\Omega = 0, \quad \forall \boldsymbol{\varphi} \in \mathbf{H}_0(\mathrm{rot}, \Omega).
$$

Using Corollary 2.2.2 in [3] on elements of $\mathbf{L}^2(\Omega)$:

$$
\nabla \times [\mathbf{B}_a \mathbf{1}_{\Omega_T} - \mu_0 \nu \mathbf{B}(\boldsymbol{\lambda})] = 0.
$$

By multiplying the above equation by $\frac{1}{\mu_0}$, we find the following expression of the adjoint problem:

$$
\nabla \times [\mathbf{M}_a \mathbf{1}_{\Omega_T} - \nu \mathbf{B}(\boldsymbol{\lambda})] = 0,
$$

where $\mathbf{M}_a$ represents the adjoint magnetization as defined in the Theorem 1.
In order to determine the adjoint system, we use the necessary optimality conditions at the first order of the Karush-Kuhn-Tucker (KKT) theorem. Thus, the directional derivative of $\mathcal{L}$ in $\mathbf{A}$ in any direction $\varphi$ must satisfy:

$$ < \frac{\partial \mathcal{L}(\mathbf{p}, \mathbf{A}, \boldsymbol{\lambda})}{\partial \mathbf{A}} | \boldsymbol{\varphi} > = 0, \forall \boldsymbol{\varphi} \in \mathbf{H}_0(\mathbf{r}\mathbf{o}\mathbf{t}, \Omega). \quad (12) $$

By considering a function $f$ defined by:

$$ \mathcal{L}(\mathbf{p}, \mathbf{A}, \boldsymbol{\lambda}) = \int_{\Omega} f(\mathbf{p}, \mathbf{A}, \boldsymbol{\lambda}) d\Omega, \quad (13) $$

where,

$$ f(\mathbf{p}, \mathbf{A}, \boldsymbol{\lambda}) = \| \mathbf{B}(\mathbf{A}) - \mathbf{B}_0 \|^2 \mathbf{1}_{\Omega_T} + \mu_0 [\boldsymbol{\lambda}.\mathbf{J}(\mathbf{p}) - \nu(\mathbf{p})\mathbf{B}(\boldsymbol{\lambda}).\mathbf{B}(\mathbf{A}) + \mathbf{B}(\boldsymbol{\lambda}).\mathbf{M}(\mathbf{p})]. $$

By developing the expression of the Lagrangian (13), and using the linearity property of integrals, we obtain:

$$ \begin{aligned} < \frac{\partial \mathcal{L}(\mathbf{p}, \mathbf{A}, \boldsymbol{\lambda})}{\partial \mathbf{A}} | \boldsymbol{\varphi} > &= \lim_{t \to 0} \frac{\mathcal{L}(\mathbf{p}, \mathbf{A} + t\boldsymbol{\varphi}, \boldsymbol{\lambda}) - \mathcal{L}(\mathbf{p}, \mathbf{A}, \boldsymbol{\lambda})}{t} \\ &= \lim_{t \to 0} \frac{\int_{\Omega} f(\mathbf{p}, \mathbf{A} + t\boldsymbol{\varphi}, \boldsymbol{\lambda}) d\Omega - \int_{\Omega} f(\mathbf{p}, \mathbf{A}, \boldsymbol{\lambda}) d\Omega}{t} \\ &= \int_{\Omega} \lim_{t \to 0} \frac{f(\mathbf{p}, \mathbf{A} + t\boldsymbol{\varphi}, \boldsymbol{\lambda}) - f(\mathbf{p}, \mathbf{A}, \boldsymbol{\lambda})}{t} d\Omega \\ &= \int_{\Omega} D_{\boldsymbol{\varphi}} f(\mathbf{A}) d\Omega \\ &= \int_{\Omega} \langle \nabla_{\mathbf{A}} f | \boldsymbol{\varphi} > d\Omega \\ &= \int_{\Omega} \frac{\partial f}{\partial A_x} \boldsymbol{\varphi}_x + \frac{\partial f}{\partial A_y} \boldsymbol{\varphi}_y + \frac{\partial f}{\partial A_z} \boldsymbol{\varphi}_z d\Omega. \end{aligned} $$

Using the derivation of compound functions, for each of the partial derivatives, where $A_x$, $A_y$ or $A_z$, then we have:

$$ \begin{aligned} \frac{\partial f}{\partial A_x} &= \frac{\partial}{\partial A_x} [ \| \mathbf{B}(\mathbf{A}) - \mathbf{B}_0 \|^2 \mathbf{1}_{\Omega_T} + \mu_0 [\boldsymbol{\lambda}.\mathbf{J} - \nu\mathbf{B}(\boldsymbol{\lambda}).\mathbf{B}(\mathbf{A}) + \mathbf{B}(\boldsymbol{\lambda}).\mathbf{M}] ] \\ &= \frac{\partial}{\partial \mathbf{B}(\mathbf{A})} (\| \mathbf{B}(\mathbf{A}) - \mathbf{B}_0 \|^2 \mathbf{1}_{\Omega_T}) \cdot \frac{\partial \mathbf{B}(\mathbf{A})}{\partial A_x} - \mu_0 \nu \mathbf{B}(\boldsymbol{\lambda}) \cdot \frac{\partial \mathbf{B}(\mathbf{A})}{\partial A_x} \\ &= \mathbf{B}_a \cdot \frac{\partial \mathbf{B}(\mathbf{A})}{\partial A_x} \mathbf{1}_{\Omega_T} - \mu_0 \nu \mathbf{B}(\boldsymbol{\lambda}) \cdot \frac{\partial \mathbf{B}(\mathbf{A})}{\partial A_x} \\ &= [\mathbf{B}_a \mathbf{1}_{\Omega_T} - \mu_0 \nu \mathbf{B}(\boldsymbol{\lambda})] \cdot \frac{\partial \mathbf{B}(\mathbf{A})}{\partial A_x}, \end{aligned} $$

with:

$$ \mathbf{B}_a = 2 \begin{pmatrix} B_x - B_{0x} \\ B_y - B_{0y} \\ B_z - B_{0z} \end{pmatrix}. $$

Thus, we obtain:

$$ < \frac{\partial \mathcal{L}(\mathbf{p}, \mathbf{A}, \boldsymbol{\lambda})}{\partial \mathbf{A}} | \boldsymbol{\varphi} > = \int_{\Omega} [\mathbf{B}_a \mathbf{1}_{\Omega_T} - \mu_0 \nu \mathbf{B}(\boldsymbol{\lambda})] \cdot \left[ \varphi_x \frac{\partial \mathbf{B}(\mathbf{A})}{\partial A_x} + \varphi_y \frac{\partial \mathbf{B}(\mathbf{A})}{\partial A_y} + \varphi_z \frac{\partial \mathbf{B}(\mathbf{A})}{\partial A_z} \right] d\Omega. \quad (14) $$
### 3. Computation of the sensitivity of the problem:

From the independence of the variables **p**, **A** and **λ**, in the definition of the Lagrangian function $\mathcal{L}$, and using the property of derivation under the integral, we obtain:

$$ \frac{\partial \mathcal{L}(\mathbf{p}, \mathbf{A}, \mathbf{\lambda})}{\partial \mathbf{p}_i} = \int_{\Omega} \mu_0 [\mathbf{\lambda}^T \frac{\partial \mathbf{J}(\mathbf{p})}{\partial \mathbf{p}_i} + \mathbf{B}(\mathbf{\lambda})^T \frac{\partial \mathbf{M}(\mathbf{p})}{\partial \mathbf{p}_i} - \nu_0 \frac{\partial \nu_r(\mathbf{p})}{\partial \mathbf{p}_i} \mathbf{B}(\mathbf{\lambda})^T \mathbf{B}(\mathbf{A})] d\Omega. \quad (17) $$

For each couple ($\mathbf{p}$,$\mathbf{A}$) satisfying constraint (3), we have $F_o(\mathbf{p}) = \mathcal{L}(\mathbf{p}, \mathbf{A}(\mathbf{p}), \mathbf{\lambda})$. Thus, using the compound derivative theorem, we obtain:

$$ \frac{\partial F_o(\mathbf{p})}{\partial \mathbf{p}_i} = \frac{\partial \mathcal{L}(\mathbf{p}, \mathbf{A}(\mathbf{p}), \mathbf{\lambda})}{\partial \mathbf{p}_i} + \frac{\partial \mathcal{L}(\mathbf{p}, \mathbf{A}(\mathbf{p}), \mathbf{\lambda})^T}{\partial \mathbf{A}} \cdot \frac{\partial \mathbf{A}(\mathbf{p})}{\partial \mathbf{p}_i}. \quad (18) $$

Note that $\frac{\partial \mathcal{L}(\mathbf{p},\mathbf{A}(\mathbf{p}),\mathbf{\lambda})}{\partial \mathbf{p}_i} = 0$ because $\mathbf{p}$ verify the adjoint system which is equivalent to the stationarity of the Lagrangian function in $\mathbf{A}$, hence, we have:

$$ \frac{\partial F_o(\mathbf{p})}{\partial \mathbf{p}_i} = \int_{\Omega} \mu_0 [\mathbf{\lambda}^T \frac{\partial \mathbf{J}(\mathbf{p})}{\partial \mathbf{p}_i} + \mathbf{B}(\mathbf{\lambda})^T \frac{\partial \mathbf{M}(\mathbf{p})}{\partial \mathbf{p}_i} - \nu_0 \frac{\partial \nu_r(\mathbf{p})}{\partial \mathbf{p}_i} \mathbf{B}(\mathbf{\lambda})^T \mathbf{B}(\mathbf{A})] d\Omega $$

for each component $\mathbf{p}_i$ of $\mathbf{p}$.

□

From Theorem 1, we obtain directly (using the derivative with respect to $p_1 = \nu_r$) the following corollary:

**Corollary 2.** The expression of the sensitivity with respect to $\nu$ is in 3D given by:

$$ \frac{\partial F_o(\nu_r)}{\partial \nu_r} = - \int_{\Omega} \mathbf{B}(\mathbf{\lambda})^T \mathbf{B}(\mathbf{A}) d\Omega. \quad (19) $$

Where $\mathbf{B}(\mathbf{\lambda})$ is obtained by solving the adjoint problem provided by equation (10) in Theorem 1.

The weak formulation for those kind of 3D magnetostatic problem is:

$$ \int_{\Omega} \nu \mathbf{B}(\mathbf{X}) \cdot \mathbf{B}(\varphi) d\Omega - \int_{\Omega} \mathbf{M} \cdot \mathbf{B}(\varphi) d\Omega - \int_{\Omega} \mathbf{J} \cdot \varphi d\Omega = 0, \forall \varphi \in \mathbf{H}_0(\text{rot}, \Omega) \quad (20) $$

Hence, by taking a simple magnetic circuit with only one coil, one ferromagnetic material part and no permanent magnet, we have to take $X = A$, $\mathbf{M} = \mathbf{0}$ and the only source is the current $\mathbf{J}$. Moreover, the adjoint problem is provided by taking $X = \mathbf{\lambda}$, $J = 0$ and the only source is now $\mathbf{M} = \mathbf{M}_a$.

It should be noticed that in the 3D case, it is necessary to add a gauge condition, in order to ensure the uniqueness of the solution using finite element method [8].

## 4 Topology optimization algorithm

In Section 3, a way to compute efficiently the derivative of the function $F_o$ is derived in Theorem 1 for 3D problems, these formulas are provided in continuous spaces. In this work, those formulas are solved numerically using finite element method. This involves discretization steps where the entire domain is discretized into small tetrahedron elements. The variables of the associated topology problem could be directly the tetrahedron elements in the variable domain $\Omega_V$ but they also can be defined by other cells, as for example hyperrectangles (see Fig.4 in the next section).

It remains a last difficulty. In some cases, the variables are discrete. Indeed, $p_1 = \nu_r$ represents the relative reluctivity in a point $(x, y, z)^T \in \Omega_V$. For a specific ferromagnetic material, $p_1$ can only takes two discrete values: $\nu_r = \nu_r^{\max} = \nu_0$ if we put air at $(x, y, z)^T$ or $\nu_r = \nu_r^{\min}$ if we put material at $(x, y, z)^T$. Note that, the other six parameters are continuous and can be directly used in a standard optimization code based on gradient descent steps.

In the next sub-section, the main method to deal with this issue is recalled. This is the well-known SIMP method[5, 4, 15, 9]. In sub-section 4.2, our optimization code based on the SIMP method is provided and explained.
Expanding the term with partial derivatives according to the components of **A**, we obtain:

$$ 
\begin{aligned}
\left[ \varphi_x \frac{\partial \mathbf{B}(\mathbf{A})}{\partial A_x} + \varphi_y \frac{\partial \mathbf{B}(\mathbf{A})}{\partial A_y} + \varphi_z \frac{\partial \mathbf{B}(\mathbf{A})}{\partial A_z} \right] &= J_{\mathbf{B}}(\mathbf{A})\varphi \\
&= D_{\varphi}\mathbf{B}(\mathbf{A}).
\end{aligned}
$$

Where $J_{\mathbf{B}}(\mathbf{A})$ is the Jacobian matrix of the function **B** with respect to **A**. For each element **A** = ($A_x, A_y, A_z$)$^T \in \mathbb{R}^3$ it is associated $\mathbf{B}(\mathbf{A}) = (B_x(\mathbf{A}), B_y(\mathbf{A}), B_z(\mathbf{A}))^T \in \mathbb{R}^3$.

Using the definition of the directional derivative and the linearity of the **curl** operator, we have:

$$ D_{\varphi}\mathbf{B}(\mathbf{A}) = \mathbf{B}(\varphi). \qquad (15) $$

From equation (14) and using equation (15), we have:

$$ \int_{\Omega} [\mathbf{B}_a \mathbf{1}_{\Omega_T} - \mu_0 \nu \mathbf{B}(\lambda)] \cdot \mathbf{B}(\varphi) d\Omega = 0. $$

By applying the Green formula in 3D, we obtain:

$$ \int_{\Omega} \varphi \cdot \nabla \times [\mathbf{B}_a \mathbf{1}_{\Omega_T} - \mu_0 \nu \mathbf{B}(\lambda)] d\Omega + \int_{\partial\Omega} (\varphi \times \mathbf{n}) \cdot [\mathbf{B}_a \mathbf{1}_{\Omega_T} - \mu_0 \nu \mathbf{B}(\lambda)] d\partial\Omega = 0. \quad (16) $$

Because $\varphi \in \mathbf{H}_0(\text{curl}, \Omega)$, $\varphi \times \mathbf{n} = \mathbf{0}$ on $\partial\Omega$, we have:

$$ \int_{\Omega} \varphi \cdot \nabla \times [\mathbf{B}_a \mathbf{1}_{\Omega_T} - \mu_0 \nu \mathbf{B}(\lambda)] d\Omega = 0, \quad \forall \varphi \in \mathbf{H}_0(\text{rot}, \Omega). $$

Using Corollary 2.2.2 in (Allaire & Alouges 2015) on elements of $\mathbf{L}^2(\Omega)$:

$$ \nabla \times [\mathbf{B}_a \mathbf{1}_{\Omega_T} - \mu_0 \nu \mathbf{B}(\lambda)] = 0. $$

By multiplying the above equation by $\frac{1}{\mu_0}$, we find the following expression of the adjoint problem:

$$ \nabla \times [\mathbf{M}_a \mathbf{1}_{\Omega_T} - \nu \mathbf{B}(\lambda)] = 0, $$

where $\mathbf{M}_a$ represents the adjoint magnetization as defined in the Theorem 1.

### 3. Computation of the sensitivity of the problem:

From the independence of the variables **p**, **A** and **λ**, in the definition of the Lagrangian function $\mathcal{L}$, and using the property of derivation under the integral, we obtain:

$$ \frac{\partial \mathcal{L}(\mathbf{p}, \mathbf{A}, \mathbf{\lambda})}{\partial \mathbf{p}_i} = \int_{\Omega} \mu_0 \left[ \mathbf{\lambda} \cdot \frac{\partial \mathbf{J}(\mathbf{p})}{\partial \mathbf{p}_i} + \mathbf{B}(\mathbf{\lambda}) \cdot \frac{\partial \mathbf{M}(\mathbf{p})}{\partial \mathbf{p}_i} - \nu_0 \frac{\partial \nu_r(\mathbf{p})}{\partial \mathbf{p}_i} \mathbf{B}(\mathbf{\lambda}) \cdot \mathbf{B}(\mathbf{A}) \right] d\Omega. \quad (17) $$

For each couple $(\mathbf{p}, \mathbf{A})$ satisfying constraint (3), we have $F_o(\mathbf{p}) = \mathcal{L}(\mathbf{p}, \mathbf{A}(\mathbf{p}), \mathbf{\lambda})$. Thus, using the chain rule theorem, we obtain:

$$ \frac{\partial F_o(\mathbf{p})}{\partial \mathbf{p}_i} = \frac{\partial \mathcal{L}(\mathbf{p}, \mathbf{A}(\mathbf{p}), \mathbf{\lambda})}{\partial \mathbf{p}_i} + \frac{\partial \mathcal{L}(\mathbf{p}, \mathbf{A}(\mathbf{p}), \mathbf{\lambda})}{\partial \mathbf{A}} \cdot \frac{\partial \mathbf{A}(\mathbf{p})}{\partial \mathbf{p}_i}. \quad (18) $$

Note that $\frac{\partial \mathcal{L}(\mathbf{p}, \mathbf{A}(\mathbf{p}), \mathbf{\lambda})}{\partial \mathbf{A}} = \mathbf{0}$ because of the stationarity of the Lagrangian function in **A**, hence, we have:

$$ \frac{\partial F_o(\mathbf{p})}{\partial \mathbf{p}_i} = \int_{\Omega} \mu_0 \left[ \mathbf{\lambda} \cdot \frac{\partial \mathbf{J}(\mathbf{p})}{\partial \mathbf{p}_i} + \mathbf{B}(\mathbf{\lambda}) \cdot \frac{\partial \mathbf{M}(\mathbf{p})}{\partial \mathbf{p}_i} - \nu_0 \frac{\partial \nu_r(\mathbf{p})}{\partial \mathbf{p}_i} \mathbf{B}(\mathbf{\lambda}) \cdot \mathbf{B}(\mathbf{A}) \right] d\Omega $$

for each component $\mathbf{p}_i$ of **p**.
## 4.1 Penalization with material density method - SIMP

Focusing on the parameter $p_1$ which is discrete, we apply the SIMP method that is recalled in the following; for details see, [15, 9]. This SIMP method could be applied in the same way on other parameters $p_i, i = 2, \dots, 7$ if necessary. Here, the other parameters $p_i, i = 2, \dots, 7$ are fixed.

The discrete parameter $p_1$ is artificially relaxed and becomes continuous between its two previous discrete values $p_1^{\min}$ and $p_1^{\max}$. This introduces composite materials between $p_1^{\min} = \nu_r^{\min}$ and $p_1^{\max} = \nu_r^{\max}$. This involves a magnetic structure which is generally not manufacturable. Hence, to prevent the appearance of these intermediate values, we penalize them introducing a density variable $\rho$ belonging to [0,1] and linked to the design variable $p_1$ via a material interpolation scheme $g$ such that, at a point $x$ of $\Omega_V$, one obtains: $g(\rho(x, y, z)) = p_1(x, y, z)$ with $g(0) = p_1^{\min} = \nu_r^{\min}$ and $g(1) = p_1^{\max} = \nu_r^{\max}$.

In the literature, there are several choices of interpolation functions $g$, that influences the final result of topology optimization. In this work, we have used two penalization functions, as shown in Fig. 2:

* Polynomial interpolation function [4]:

$$p_1 = g(\rho) = p_1^{\min} + (p_1^{\max} - p_1^{\min})\rho^n, \quad (21)$$

where $n > 0$ is a penalty parameter. When the penalization parameter $n$ is high we enforce the solution to be discrete close to $p_1 \in \{p_1^{\min}, p_1^{\max}\}$. The main problem with this scheme is not symmetrical and hence it favors materials with low values of $\rho$.

* Arctangent interpolation function (introduced by D. Lukàš in [11]):

$$p_1 = g(\rho) = p_1^{\min} + (p_1^{\max} - p_1^{\min}) \left[ \frac{1}{2 \arctan(d)} \arctan(d(2\rho - 1)) + \frac{1}{2} \right], \quad (22)$$

where $d > 0$ is the penalty parameter of this scheme. This interpolation avoids the symmetry problems of polynomial penalization. High values of $d$ allows penalizing intermediate materials but this can also lead to a poor convergence of the optimization algorithm.

In Fig. 2, the two used interpolation schemes are represented.

<figure><img src="image_4.png" /><figcaption>Figure 2: Polynomial and arctangent interpolation schemes for several penalization degrees</figcaption></figure>

Therefore, the new variable of problem ($\mathcal{P}$) is now the material density $\rho$ and then, using chain rule and from Theorem 1 in the case where $p_1 = \nu_r$, the new sensitivity with respect to $\rho$ becomes:

$$\frac{\partial F_o}{\partial \rho} = - \int_{\Omega} \frac{\partial g}{\partial \rho} \mathbf{B}(\lambda)^T \mathbf{B}(\mathbf{A}) d\Omega \quad (23)$$
<figure><img src="image_2.png" /></figure>

From Theorem 1, we obtain directly (using the derivative with respect to $p_1 = v_r$) the following corollary:

**Corollary 2.** The expression of the sensitivity with respect to $v$ is in 3D given by:

$$ \frac{\partial F_o(v_r)}{\partial v_r} = - \int_{\Omega} \mathbf{B}(\lambda) \cdot \mathbf{B}(\mathbf{A}) d\Omega. \qquad (19) $$

Where $\mathbf{B}(\lambda)$ is obtained by solving the continuous adjoint problem provided by equation (10) in Theorem 1.

The weak formulation for those kind of 3D magnetostatic problem is:

$$ \int_{\Omega} \nu \mathbf{B}(\mathbf{X}) \cdot \mathbf{B}(\varphi) d\Omega - \int_{\Omega} \mathbf{M} \cdot \mathbf{B}(\varphi) d\Omega - \int_{\Omega} \mathbf{J} \cdot \varphi d\Omega = 0, \qquad (20) $$

where $\varphi$ is a test function of $\mathbf{H}_0(\text{rot}, \Omega)$ to define the weak formulation in the finite element method.

Hence, by taking a simple magnetic circuit with only one coil, one ferromagnetic material part and no permanent magnet, we have to take $\mathbf{X} = \mathbf{A}$, $\mathbf{M} = \mathbf{0}$ and the only source is the current $\mathbf{J}$. Moreover, the adjoint problem is provided by taking $\mathbf{X} = \lambda$, $\mathbf{J} = \mathbf{0}$ and the only source is now $\mathbf{M} = \mathbf{M}_a$.

It should be noticed that in the 3D case, it is necessary to add a gauge condition, in order to ensure the uniqueness of the solution using finite element method; this has to be done directly inside the 3D-finite element method solver.

# 4 Topology optimization algorithm

In Section 3, a way to compute efficiently the derivative of the function $F_o$ is derived in Theorem 1 for 3D problems, these formulas are provided in continuous spaces. In this work, those formulas are solved numerically using finite element method. This involves discretization steps where the entire domain is discretized into small tetrahedral elements. The variables of the associated topology problem could be directly the tetrahedral elements in the variable domain $\Omega_V$ but they can also be defined by other cells, as for example hyperrectangles.

It remains a last difficulty. In some cases, the variables are discrete. Indeed, $p_1 = v_r$ represents the relative reluctivity in a point $(x, y, z)^T \in \Omega_V$. For a specific ferromagnetic material, $p_1$ can only takes two discrete values: $v_r = v_r^{\max} = v_0$ if we put air at $(x, y, z)^T$ or $v_r = v_r^{\min}$ if we put material at $(x, y, z)^T$. Note that, the other six parameters are continuous and can be directly used in a standard optimization code based on gradient descent steps.

In the next sub-section, the main method to deal with this issue is recalled. This is the well-known SIMP method (Bendsøe & Sigmund 2003, Bendsøe & Sigmund 1999, Sanogo & Messine 2018, Gauthey et al. 2021). In sub-section 4.2, our optimization code based on the SIMP method is provided and explained.

## 4.1 Penalization with material density method - SIMP

Focusing on the parameter $p_1$ which is discrete, we apply the SIMP method that is recalled in the following; for details see, (Sanogo & Messine 2018, Gauthey et al. 2021). This SIMP method could be applied in the same way on other parameters $p_i, i = 2, \dots, 7$ if necessary. Here, the other parameters $p_i, i = 2, \dots, 7$ are fixed; note, that those remaining parameters are continuous and can be taken into account directly (without the use of SIMP) and hence, the main difficulty concerns the parameter $p_1$ which is discrete.
## 4.2 Optimization algorithm

In order to solve equation 23, the domain $\Omega_V$ is discretized into $N$ cells $\Omega_{V_i}$, each of which contains matter or not according to the value of the density $\rho_i \in [0, 1], i = 1, \dots, N$ in $\Omega_{V_i}$. Then, we pose: $\rho = (\rho_1, \dots, \rho_N)^T$.

We can write the density variable $\rho$ in all space as follows:

$$ \rho = \sum_{i=1}^{N} \rho_i \mathbb{1}_{\Omega_{V_i}} + \rho_{\Omega \setminus \Omega_V} \mathbb{1}_{\Omega \setminus \Omega_V} \quad (24) $$

with $\rho_{\Omega \setminus \Omega_V}$ equals 1 in areas with ferromagnetic material and 0 in areas without ferromagnetic material outside the variable domain $\Omega_V$.

Hence, using again chain rule and Theorem 1 in the case where $p_1 = \nu_r$, the sensitivity with respect to each element $\rho_i$ of $\rho$ becomes:

$$ \frac{\partial F_o}{\partial \rho_i} = - \int_{\Omega_{V_i}} \frac{\partial g}{\partial \rho} \mathbf{B}(\lambda)^T \cdot \mathbf{B}(\mathbf{A}) d\Omega. \quad (25) $$

In order to solve problem (P), the optimization algorithm needs the following information:

* The objective function and sensitivity calculation (see equation (25)).

* The constraints which are only the lower and upper bounds of the material density: $\rho \in [0, 1]^N$.

* The constraint imposed by the Maxwell equations is indirectly taken into account using the finite element method.

Our algorithm can be summarized in Fig. 3. Like any local optimization algorithm based on gradient descent, this optimization module starts from a starting point denoted by $\rho^{[0]}$, build a sequence of points $\rho^{[k]}$ and ends to $\rho^{opt}$. Then, this point $\rho^{opt}$ is analyzed because this solution is not necessarily discrete (does not belong necessarily to the set $\{0, 1\}^N$). Hence, a test is added to check if $\rho^{opt}$ is close to one of its two limits 0 or 1 with an imposed tolerance tol. If the tolerance tol is satisfied, the solution provided is $\rho^{opt}$, otherwise the new starting point is updated by doing $\rho^{[0]} = \rho^{opt}$) and by increasing the penalty parameter. This makes it possible to converge slowly but strongly to efficient design solutions. Indeed, by directly penalizing with a big parameter $n$ or $d$, the optimization module converges to a local minimum which is very far from those obtained using our homotopy algorithm.

<figure><img src="image_5.png" /><figcaption>Figure 3: SIMP based topology optimization algorithm</figcaption></figure>
The discrete parameter $p_1$ is relaxed and becomes continuous between its two previous discrete values $p_1^{\min}$ and $p_1^{\max}$. This introduces composite materials between $p_1^{\min} = \nu_r^{\min}$ and $p_1^{\max} = \nu_r^{\max}$. This involves a magnetic structure which is generally not manufacturable. Hence, to prevent the appearance of these intermediate values, we penalize them introducing a density variable $\rho$ belonging to $[0,1]$ and linked to the design variable $p_1$ via a material interpolation scheme $g$ such that, at a point $x$ of $\Omega_V$, one obtains: $g(\rho(x, y, z)) = p_1(x, y, z)$ with $g(0) = p_1^{\min} = \nu_r^{\min}$ and $g(1) = p_1^{\max} = \nu_r^{\max}$.

In the literature, there are several choices of interpolation functions $g$, that influences the final result of topology optimization, (Bendsøe & Sigmund 1999, Bendsøe & Sigmund 2003, Sanogo & Messine 2018, Gauthey et al. 2021, Lucchini et al. 2022). In this work, plenty penalization approaches were tested, however we present only results coming from the use of the polynomial interpolation scheme, as defined in (Bendsøe & Sigmund 1999):

$$p_1 = g(\rho) = p_1^{\min} + (p_1^{\max} - p_1^{\min})\rho^n, \qquad (21)$$

where $n > 0$ is a penalty parameter. When the penalization parameter $n$ is high, we enforce the solution to be discrete close to $p_1 \in \{p_1^{\min}, p_1^{\max}\}$. The main difficulty with this scheme is that it is not symmetrical and hence it favors materials with low values of $\rho$.

Therefore, the new variable of problem ($\mathcal{P}$) is now the material density $\rho$ and then, using the chain rule theorem and from Theorem 1 in the case where $p_1 = \nu_r$, the new sensitivity with respect to $\rho$ becomes:

$$\frac{\partial F_o}{\partial \rho} = - \int_{\Omega} \frac{\partial g}{\partial \rho} \mathbf{B}(\lambda) \cdot \mathbf{B}(\mathbf{A}) d\Omega \qquad (22)$$

## 4.2 Optimization algorithm

In order to solve equation 22, the domain $\Omega_V$ is discretized into $N$ cells $\Omega_{V_i}$, each of which contains matter or not, according to the value of the density $\rho_i \in [0, 1]$, $i = 1, \dots, N$ in $\Omega_{V_i}$. Then, by denoting $\rho = (\rho_1, \dots, \rho_N)^T$, we can write the density variable $\rho$ in all space as follows:

$$\rho = \sum_{i=1}^{N} \rho_i \mathbb{1}_{\Omega_{V_i}} + \rho_{\Omega \setminus \Omega_V} \mathbb{1}_{\Omega \setminus \Omega_V} \qquad (23)$$

with $\rho_{\Omega \setminus \Omega_V}$ equals 1 in areas with ferromagnetic material and 0 in areas without ferromagnetic material outside the variable domain $\Omega_V$.

Hence, using again the chain rule theorem and Theorem 1 in the case where $p_1 = \nu_r$, the sensitivity with respect to each element $\rho_i$ of $\rho$ becomes:

$$\frac{\partial F_o}{\partial \rho_i} = - \int_{\Omega_{V_i}} \frac{\partial g}{\partial \rho} \mathbf{B}(\lambda) \cdot \mathbf{B}(\mathbf{A}) d\Omega. \qquad (24)$$

In order to solve problem ($\mathcal{P}$), the optimization algorithm needs the following information:

* The objective function and sensitivity calculation (see equation (24)).

* The constraints which are only the lower and upper bounds of the material density: $\rho \in [0, 1]^N$.

* The constraint imposed by the Maxwell equations is indirectly taken into account using the finite element method.

Our algorithm can be summarized in Fig. 1. Like any local optimization algorithm based on gradient descent, this optimization module starts from a starting point denoted by $\rho^{[0]}$, build a sequence of points $\rho^{[k]}$ and ends to $\rho^{opt}$. Then, this point $\rho^{opt}$ is analyzed because this solution is not necessarily discrete (does not belong necessarily to the set $\{0, 1\}^N$). Hence, a test is added to check if $\rho^{opt}$ is close to one of its two limits 0 or 1 with an imposed tolerance $tol$. If the tolerance $tol$ is satisfied, the solution provided is $\rho^{opt}$, otherwise the new starting point is updated by doing
## 4.3 Filtering technique

If we consider the expression of the derivative of $F_o$ with respect to each element of $\rho$, shown in equation (25), the final solution has generally checkerboard patterns, [17]. In order to make the solution more manufacturable, a regularization method must be added. Thus, in this work, a sensitivity filtering [17] method has been used. The latter consists in modifying the values of the sensitivity of the cost function at each iteration through a spatial filter whose goal is to weight the value of the sensitivity in a target element by the value of its neighboring elements. In this work, several sensitivity filters have been tested. For our problem the bi-lateral sensitivity filter provided the most efficient discrete solutions, [17], and so, the expression of the sensitivity with respect to an element $e$ is replaced by:

$$ \frac{\partial F_o}{\partial \rho_e} = \frac{\sum_{i \in N_e} w_{ie} \overline{w} \left( \frac{\partial F_o}{\partial \rho_i} \right) v_i \frac{\partial F_o}{\partial \rho_i}}{\sum_{i \in N_e} w_{ie} \overline{w} \left( \frac{\partial F_o}{\partial \rho_i} \right) v_i} \quad (26) $$

the sensibility weights are computed as follows:

$$ \overline{w} \left( \frac{\partial F_o}{\partial \rho_i} \right) = \exp \left[ -\frac{1}{2} \left( \frac{\frac{\partial F_o}{\partial \rho_i} - \frac{\partial F_o}{\partial \rho_e}}{\sigma_r} \right)^2 \right] \quad \text{and} \quad \sigma_r \in [0, 1] $$

Moreover, $\sigma_r$ is fixed at 0.5, and the distance weighting $w_{ie}$ is equal to $\max(0, r_{min_e} - \text{dist}(i, e))$. The following parameters are defined in this work such as:

* $v_i$ is the volume of the element $i$.

* $\text{dist}(i, e)$ is the distance between the center of the element $i$ and the center of the element $e$.

* $r_{min_e} = \delta h_e$; with $\delta \in \mathbb{R}_+$ and $h_e$ is the greatest distance between the points constituting the element $e$ (Euclidean norm of the difference between the coordinates of the most distant nodes of the element $e$).

* The set of neighbors of the element $e$ is defined as follows: $N_e = \{i | \text{dist}(i, e) \le r_{min_e}\}$.

Thus, in this section an efficient topology optimization algorithm based on adjoint methods and SIMP penalization technique is provided. Moreover, our methodology is able to provide more regularized solution close to a manufacturable one.

# 5 Numerical validations on a 3D magnetic circuit

The aim of this section is to validate our methodology on a U-shaped magnetic circuits. In a first step, a simple design problem with 16 variables is addressed. This allows us to verify that the computation of the gradient based on our theoretical results (mainly Theorem 1 and Corollary 2) is close to the one obtained using the finite difference method. Then, in a second step, we compare two penalization functions: the polynomial and the arctangent interpolation functions. Finally, in a last sub-section, we expand our study by addressing a design problem with a large number of variables.

## 5.1 Description of the design problems and numerical simulations

We consider here a simple magnetic circuit shown in Fig. 4 with one coil and two distinct variable zones which can be air or ferromagnetic materials. For this magnetic circuit, the design variable $p_1$ is the relative reluctivity $\nu_r$ which varies between the minimum value $p_1^{\min} = \nu_r^{\min} = 1/1000$ (corresponding to ferromagnetic material) and the maximum value $p_1^{\max} = \nu_r^{\max} = 1$ (corresponding to air).

For this example, the magnetic circuit can be divided into 4 parts, as follows:
$\rho^{[0]} = \rho^{opt}$) and by increasing the penalty parameter. This makes it possible to converge slowly but strongly to efficient design solutions. Indeed, by directly penalizing with a big parameter $n$, the optimization module converges to a local minimum which is very far from those obtained using our homotopy algorithm.

<figure><img src="image_3.png" /><figcaption>Figure 1: SIMP based topology optimization algorithm. (Source: figure created by authors)</figcaption></figure>

## 4.3 Filtering technique

If we consider the expression of the derivative of $F_o$ with respect to each element of $\rho$, shown in equation (24), the final solution has generally checkerboard patterns, (Sigmund 2007). In order to make the solution more manufacturable, a regularization method must be added. Thus, in this work, a sensitivity filtering method has been used, (Sigmund 2007). The latter consists in modifying the values of the sensitivity of the cost function at each iteration through a spatial filter whose goal is to weight the value of the sensitivity in a target element by the value of its neighboring elements. In this work, several sensitivity filters have been tested, (Sigmund 2007, Sigmund & Maute 2013, Lucchini et al. 2022). For our problem the bi-lateral sensitivity filter, defined in (Sigmund 2007), provided the most efficient discrete solutions and so, the expression of the sensitivity with respect to an element $e$ is replaced by:

$$ \frac{\partial F_o}{\partial \rho_e} = \frac{\sum_{i \in N_e} w_{ie} \bar{w} \left( \frac{\partial F_o}{\partial \rho_i} \right) v_i \frac{\partial F_o}{\partial \rho_i}}{\sum_{i \in N_e} w_{ie} \bar{w} \left( \frac{\partial F_o}{\partial \rho_i} \right) v_i} \quad (25) $$

the sensibility weights are computed as follows:

$$ \bar{w} \left( \frac{\partial F_o}{\partial \rho_i} \right) = \exp \left[ -\frac{1}{2} \left( \frac{\frac{\partial F_o}{\partial \rho_i} - \frac{\partial F_o}{\partial \rho_e}}{\sigma_r} \right)^2 \right] \quad \text{and} \quad \sigma_r \in [0, 1]. $$

Moreover, $\sigma_r$ is fixed at 0.5, and the distance weighting $w_{ie}$ is equal to $\max(0, r_{min_e} - \text{dist}(i, e))$. The following parameters are defined in this work such as:

* $v_i$ is the volume of the element $i$.
* A U-shaped ferromagnetic part which defines the fixed part of the magnetic circuit. It is represented in white on Fig. 4. Its dimensions are 105mm × 40mm × 15mm, the thickness of the bottom hyperrectangle containing the coil is 10mm and the thickness of the both sides of the U-shaped circuit is 15mm.

* A coil which is represented in yellow on Fig. 4. It is made up from 2 bottom and top bars of 15mm × 0.4mm × 5mm and 2 lateral bars of 10mm × 0.4mm × 5mm which are linked by 4 curved bars. This coil is crossed by a current of 0.4A and has 200 windings.

* Two design variable zones which are hyperrectangles of size 30mm × 15mm × 15mm which are represented in grey on Fig. 4. Each of these two parts is discretized into 8 variable cells that contain air or ferromagnetic material.

* One target zone composed by small cubes of size 1mm and represented in black on Fig. 4. These cubes are centered in the air gap between the two variable hyperrectangles. The target zone $\Omega_T$ is composed by $27 = 3 \times 3 \times 3$ cubes, where the magnetic flux density is evaluated.

In order to solve the primal problem (3) and the adjoint problem (17), we use a finite element method solver which is GetDP software, [7]. To use GetDP, it is necessary to specify some elements of our problems: the Sobolev space, the weak formulation of the problem with the boundary conditions, the Jacobian and the type of numerical integration.

GetDP is based on the mesher Gmsh, [10, 12] and the mesh of a domain owns some great importance because this involves the accuracy obtained on the solution versus the calculation time; so a previous study has been done to obtain efficiently a good accuracy on the solutions.

In this first example, we use two distinct meshes, as follows:

* For the calculation of magnetic flux density: we use a fine tetrahedral mesh to ensure the accuracy of numerical values. This tetrahedral mesh is made on the whole study area $\Omega$ using Gmsh.

* For the optimization cell variables: we discretize the variable area $\Omega_V$ (grey area on Fig. 4) into small cells. Each cell represents a variable. For this first example, the variable zone $\Omega_V$ is made up by 16 cubes (2 × 2 × 2 in each side).

<figure><img src="image_6.png" /><figcaption>Figure 4: The studied 3D magnetic circuit</figcaption></figure>

For this first example, the given magnetic flux density $B_0$ in the target zone $\Omega_T$ is first computed using a symmetric ferromagnetic shape as it is represented in Fig. 5.

The corresponding magnetic flux density is plotted in Fig. 6.

## 5.2 Validation of the sensitivity calculations

In order to validate the accuracy of the adjoint based approach as it was defined in Theorem 1 (without taking into account the SIMP penalization technique), the first example presented in the previous subsection is considered but by dividing again by two all the optimization variable cells.
* $\text{dist}(i, e)$ is the distance between the center of the element $i$ and the center of the element $e$.

* $r_{\text{min}_e} = \delta h_e$; with $\delta \in \mathbb{R}_+^*$ and $h_e$ is the greatest distance between the points constituting the element $e$ (Euclidean norm of the difference between the coordinates of the most distant nodes of the element $e$).

* The set of neighbors of the element $e$ is defined as follows: $N_e = \{i \mid \text{dist}(i, e) \le r_{\text{min}_e}\}$.

Thus, in this section an efficient topology optimization algorithm based on adjoint methods and SIMP penalization technique is provided. Moreover, our methodology tends to provide more regularized solution close to manufacturable ones.

# 5 Numerical verifications on a 3D magnetic circuit

The aim of this section is to validate our methodology on a U-shaped magnetic circuit. This allows us to verify that the computation of the gradient based on our theoretical results (mainly Theorem 1 and Corollary 2) is close to the one obtained using the finite difference method. In a last sub-section, we address a design problem with a large number of variables (over 100,000 design variables).

## 5.1 Description of the design problems and numerical simulations

We consider here a simple magnetic circuit shown in Fig. 2 with one coil and two distinct variable zones which can be air or ferromagnetic materials. For this magnetic circuit, the design variable $p_1$ is the relative reluctivity $\nu_r$ which varies between the minimum value $p_1^{\min} = \nu_r^{\min} = 1/1000$ (corresponding to ferromagnetic material) and the maximum value $p_1^{\max} = \nu_r^{\max} = 1$ (corresponding to air).

<figure><img src="image_4.png" /><figcaption>Figure 2: 3D magnetic circuit design problem. (Source: figure created by authors)</figcaption></figure>

For this example, the magnetic circuit can be divided into 4 parts, as follows:

* A U-shaped ferromagnetic part which defines the fixed part of the magnetic circuit. It is represented in white on Fig. 2. Its dimensions are $105mm \times 40mm \times 15mm$, the thickness of the bottom hyperrectangle containing the coil is $10mm$ and the thickness of the both sides of the U-shaped circuit is $15mm$.

* A coil which is represented in yellow on Fig. 2. It is made up from 2 bottom and top bars of $15mm \times 0.4mm \times 5mm$ and 2 lateral bars of $10mm \times 0.4mm \times 5mm$ which are linked by 4 curved bars. This coil is crossed by a current density of $40A/mm^2$. This current density
<figure><img src="image_7.png" /><figcaption>Figure 5: Material distribution providing the magnetic flux density **B**₀. Black areas contain material with relative reluctivity *ν*ₜ = 1/1000 (a ferromagnetic material with magnetic permeability *μ*ₜ = 1000) and the white areas represent air (*ν*ₜ = 1).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6: Shape of the magnetic flux density **B**₀</figcaption></figure>

This yields $4^3 = 64$ variables for the two distinct zones $\Omega_V$ and thus, the objective function $F_o$ depends on 128 variables. Hence, the 128 partial derivatives based on our adjoint method were compared to those computed using a standard finite difference method. The relative percentage differences between the two gradient of $F_o$ are shown in Fig. 7. The maximum value of the relative percentage difference is less than 0.03%. The maximal errors are mainly concentrated in cell variables close to the corners of the hyperrectangles because the evaluation of the magnetic flux density is more difficult in these zones (even for fine meshes).

Therefore, according to those numerical results, the adjoint method is validated and thus, it can be used to provide efficient information about the gradient that will be used in the optimization algorithms.

## 5.3 Validation of our topology optimization algorithm based on the two penalization techniques

To solve the topology design problem ($\mathcal{P}$), a gradient based algorithm is used. At each iteration of the optimizaiton algorithm, the gradient is computed using the adjoint method and hence, a descent step is so-obtained. In this paper, we choose to use the solver MMA (Method of Moving Asymptots), [19]; MMA is a nonlinear optimization code with a MatLab implemantation in open-source. MMA is a local search based algorithm and hence, a starting point has to be given. At each iteration of MMA, a convex approximation of the problem ($\mathcal{P}$) is computed. These approximations are mainly based on the gradients at the current iteration but also on some information of previous iterations. The approximated convex optimization problem is then solved and yields a unique optimal solution which becomes the new point at this iteration. Note that, MMA can take into account constraints. Here, problem ($\mathcal{P}$) has only bound constraints.

For this first tests, the starting point is $\rho^{[0]} = 0.1$ (ferromagnetic material) if it corresponds to a variable cell in the left pole of $\Omega_V$ and $\rho^{[0]} = 1$ (air) for the right pole, as shown in Fig. 8.

For the penalization, both arctangent and polynomial interpolation functions are tested. Each
is high for a real-life coil but the main idea here is to create a magnetic field in the target zone in order to validate our method whatever the origin of this field.

* Two design variable zones which are hyperrectangles of size $30mm \times 15mm \times 15mm$ which are represented in grey on Fig. 2. Each of these two parts is discretized into small variable cells that contain air or ferromagnetic material.

* One target zone which is represented in black on Fig. 2. It is bisected into small cubes of size $1mm$ into 3 lines, which are centered in the air gap between the two variable regions.

In order to solve the primal problem (3) and the adjoint problem (17), we use a finite element method solver which is GetDP software, (Dular & Geuzaine n.d.). To use GetDP, it is necessary to specify some elements of our problems: the Sobolev space, the weak formulation of the problem with the boundary conditions, the Jacobian and the type of numerical integration. GetDP is based on the mesher Gmsh, (Geuzaine & Remacle 2009, Remacle & Geuzaine n.d.) and the mesh of a domain owns some great importance because this involves the accuracy obtained on the solution versus the calculation time; so a previous study has been done to obtain efficiently a good accuracy on the solutions.

To solve instances based on this U-shape circuit, we can use two distinct meshes, as follows:

* For the calculation of magnetic flux density: we use a tetrahedral mesh in $\Omega$. In order to ensure the accuracy of numerical values, this tetrahedral mesh is made on the whole study area $\Omega$ using Gmsh but is finer in $\Omega_T$. Note that, a sufficiently accurate mesh has to be taken in all $\Omega$ in order to obtain a convergence on the computations of the flux density in $\Omega_T$.

* For the optimization cell variables: we can discretize the variable area $\Omega_V$ (grey area on Fig. 2) into small hyperrectangle cells. Each cell represents a variable. Another possibility is to take the tetrahedral elements inside $\Omega_V$, obtained from a finite element mesh, directly as variables.

**Remark 1.** In order to validate the derivative, $\Omega_V$ is divided into 128 hyperrectangles ($4 \times 4 \times 4$ in each side). This involves that the objective function $F_o$ depends on 128 variables. Hence, the 128 partial derivatives based on our adjoint method were compared to those computed using a standard finite difference method for a material-filled configuration. The maximum value of the relative percentage difference between the two methods is less than 0.03%. We remark that the maximal errors are mainly concentrated in cell variables close to the corners of the hyperrectangles because the evaluation of the magnetic flux density is more difficult in these zones (even for fine meshes). Therefore, according to those numerical results, the continuous adjoint method as it is defined in Theorem 1 (without taking into account the SIMP penalization technique), is validated and thus, it can be used to provide efficient information about the gradient that will be used in topology optimization algorithms.

**Remark 2.** The adjoint method provides an adjoint problem which has the same difficulty to be solved than the primal problem (here the PDE Maxwell equation in magnetostatics). Thus, to perform the gradient of a function depending on $N$ design variables by using a finite difference method, we need to solve $N + 1$ times the original PDE problem vs the resolution of 2 PDE problems for the adjoint method (the primal and the adjoint PDE problems). Note that using the adjoint method, the computation of the gradient is independent of the number of design variables.

## 5.2 Design problem with a large number of variables

In this subsection, the goal is to solve problem (P) on the magnetic circuit described in subsection 5.1 by considering a large number of design variables. In order to achieve that, the tetrahedrons of the finite element mesh in the variable zone $\Omega_V$ are taken as cell variables of the optimization design. This yields 106,801 design variables; note that the simulation problem owns 192,855 meshes, and remark that the variable region owns thinner meshes.
<figure><img src="image_9.png" /><figcaption>Figure 7: Validation of the adjoint method by comparison with finite differences</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 8: Structure of the starting point $\rho^{[0]}$</figcaption></figure>

of these functions has a particular penalty parameter ($d$ for arctangent and $n$ for polynomial). A common parameter have to be introduced in order to compare the penalization factor of those two functions. Thus, we have chosen to use the concept of penality bandwidth $bd$ at $\epsilon$ (in %) as illustrated in Fig. 9. In the following, the method used to find $n$ and $d$ for a given bandwidth $bd$ at $\epsilon$ is provided:

* Give $\epsilon$

* Compute the tolerance: $tol := \epsilon \times (p_1^{\max} - p_1^{\min})$. This tolerance is used to define the zone where the variable cells are well defined (air or ferromagnetic material within a tolerance $tol$).

* Give a bandwidth $bd$

* Determine $n$ or $d$ by solving the nonlinear equation:

$$bd = \rho_{\max} - \rho_{\min} = g^{-1}(p_{1}^{\max} - tol) - g^{-1}(p_{1}^{\min} + tol)$$

where $g^{-1}$ is the inverse function of the penalty function defined in equation (21) and (22) with $n$ or $d$. To solve this non linear equation, we used `vpasolve` of MatLab.

* Return $n$ or $d$ corresponding to the wished bandwidth $bd$ at $\epsilon$.

This algorithm provide an efficient way to compare rigorously penalization schemes. If the penalty bandwidth $bd$ decreases then the penalization factor $d$ or $n$ (respectively for arctangent or polynomial penalization functions) increases. For example, for $bd \in \{0.75, 0.50, 0.25\}$, we obtain $n \in \{2.0, 4.2, 10\}$ and $d \in \{2.3, 6.7, 19\}$. We can now compare these two penalization functions by considering the same bandwidth (i.e, the same range width of $p$ to describe the intermediate solutions $[p_{\min} + tol, p_{\max} - tol]$).

First we solve the design problem without penalization $[p_1 = p_1^{\min} + (p_1^{\max} - p_1^{\min}).\rho]$ with a starting point $\rho^{[0]}$ and then, we decrease the bandwidth at 5% for both polynomial and arctangent
In this design example, the target magnetic flux density $\mathbf{B}_0$ is a field carried only along the $z$ direction with a Gaussian shape according to the coordinate $x$:

$$ \mathbf{B}_0 = \begin{pmatrix} 0 \\ 0 \\ B_{0z}(x) \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ Ae^{-\frac{(x-x_0)^2}{2\sigma^2}} \end{pmatrix}, $$

where $\sigma = 10^{-2}$ is the standard deviation of the Gaussian, $A = 4 \times 10^{-3}$ is a constant equal to the peak value of the Gaussian, and $x_0$ is the place in space where the peak is reached. Note that these kinds of imposed magnetic flux densities can be found in many applications, as for example for designing Hall effect thrusters, (Sanogo 2016, Sanogo & Messine 2018, Sanogo et al. 2014, Youness & Messine 2019b, Youness & Messine 2019a).

The choice of the penalty parameter $n$ is very difficult because, if $n$ is very large at the beginning, this leads to converge to a local minimum which will be not so interesting. Thus, we carry out a continuation process consisting in progressively increasing the penalization degree from $n = 1$ to $n = 100$. In this work, we take $n \in \{1, 2, 3, 4, 10, 20, 100\}$. Note that, for a fixed value of $n$ the optimization module stops either if the maximum number of iterations 200 (400 for $n = 1$) is reached or if the stopping criteria is $||\rho^{[k+1]} - \rho^{[k]}||_{\infty} \le 10^{-4}$ is satisfied.

In this paper, we choose to use the solver MMA (Method of Moving Asymptots), (Svanberg 1987); MMA is a nonlinear optimization code with a Matlab implementation in opensource. MMA is a local search based algorithm and hence, a starting point has to be given. At each iteration of MMA, a convex approximation of the problem ($\mathcal{P}$) is computed. These approximations are mainly based on the gradients at the current iteration but also on some information of previous iterations. The approximated convex optimization problem is then solved and yields a unique optimal solution which becomes the new point at this iteration. Note that MMA can take into account constraints. Here, problem ($\mathcal{P}$) has only bound constraints.

Increasing the penalty degree $n$, makes it possible to enforce our method to converge slowly but with more trust to a very efficient (local) optimal design which is close to be manufacturable. In Table 1, we provide all the intermediate numerical results obtained at each distinct steps of our homotopy algorithm (see Fig. 1). Note that the termination test is always the maximal number of iterations, providing so, optimized solution but not local optimal ones. We can remark that even if the values of the objective function decrease during the iterations, it is not the same for the number of intermediate values. Indeed, this number decreases until $n = 10$ and then increases a bit, see the fourth column of Table 1. Moreover, note that for the three optimizations, the final number of iterations is 1600 corresponding to a CPU-time of about 1 day of computations on a server (DELL server with 64GB of memory and with an Intel(R) Xeon(R) Gold 5218R CPU of 2.1GHz and 40 processors). We also note that the three optimizations carried out with or without the filter technique and using two distinct starting points lead to the same design solution with only a few pixels difference. The optimization without the filter brings the checkboarder effect into play, (Lucchini et al. 2022). Moreover, we also remark that the value of the objective function on the non-filtered solution is greater than the value of the filtered one. Hence, we think that this is due to the fact that using a filtering technique, the realizable domain is reduced and thus, in our design case, the optimization solver converges to a more interesting local minimum; i.e., it seems that without the use of a filtering technique, the algorithm stops too prematurely to a local minimum less interesting than the filtered one. In addition, we remark that on our design problem the starting point has not a big impact on the so obtained local solutions.

Thus, the best optimized design of the magnetic circuit, that we obtained so far, is presented in Fig. 3. This optimized design corresponds to the solution obtained at iteration 200 with the penalization parameter $n = 100$ when the starting point is full of air. Note that this solution is not close to be entirely discrete and several variables do not respect a tolerance about 5% of $p_1^{\min}$ and $p_1^{\max}$. Indeed, 11566 variables from 106801 are in $]p_1^{\min}+5\%,p_1^{\max} - 5\%[$ which represents 10.8% of all the cell variables. Remark that those intermediate cell variables which have different
<figure><img src="image_11.png" /><figcaption>Figure 9: Exemple of penalty bandwidth of 0.5 at 5% for both polynomial and arctangent penalization functions.</figcaption></figure>

<table><caption>Table 1: Circuit optimization results using polynomial and arctangent interpolation functions</caption><thead><tr><th>Penalization</th><th>Polynomial</th><th>Arctangent</th></tr></thead><tbody><tr><td>F<sub>o</sub>(ρ<sup>opt</sup>)</td><td>0.042390</td><td>0.042391</td></tr><tr><td>Number of iterations</td><td>378 (178+200)</td><td>778 (178+200+200+200)</td></tr><tr><td>CPU time (PC-8GB-2.1GHz)</td><td>17.74h</td><td>35.61h</td></tr><tr><td>Glob. Sol. (Fig. 5)</td><td>true</td><td>true</td></tr></tbody></table>

functions to the following values: $bd \in \{0.75, 0.50, 0.25\}$. For each new bandwidth value, the starting point is set to the solution of the previous step.

The results are shown in Table 1. Both penalizations lead to the same global minimum (shown in Fig. 5). Note that, the limit about the number of iteration that MMA can takes is fixed to 200 iterations; the MMA stopping criteria is $||ρ^{[k+1]} - ρ^{[k]}||_{\infty} \le 1.1 \times 10^{-5}$ (note that MMA does not converge in 200 iterations if we take $1 \times 10^{-5}$). For the unpenalized first step, the tolerance stopping criterion of our algorithm in Fig. 3, is not reached and then one supplementary iteration is done by taking $n = 2.0$ and $d = 2.3$ for respectively the polynomial and arctangent penalization schemes. For this second step, with the polynomial scheme, MMA stopped in 200 iterations, providing the correct global solution (see Fig. 5) within a tolerance $tol = 5%(p_1^{max} - p_1^{min})$ which stops our homothopy topology optimization algorithm (see Fig. 3). Considering the arctangent scheme, MMA stopped also in 200 iterations which is the given limit, but without reaching the tolerance $tol$. Then, two other iterations are performed corresponding to the values 6.7 and 19 of the penalization parameter $d$ but the tolerance $tol$ is still not reached. Nevertheless, the obtained solution is also the global one, see Fig. 5. These results are summarized in Table 1.

To conclude, both algorithms based on polynomial and arctangent penalization schemes provide the same global solution: the polynomial scheme used one step of penalization $n = 2.0$ and the arctangent scheme used all the values proposed $d \in \{2.3, 6.7, 19\}$ without reaching the tolerance of our homothopy algorithm. However, polynomial penalization allows us to improve the convergence of our topological optimization code (see Fig. 3) more efficiently. This shows that on this example and using MMA code, the polynomial penalization scheme based method surpasses the one based on arctangent.
<table><caption>Table 1: Numerical results using algorithm Fig. 1 to solve a 3D-design problem with over 100,000 variables, for three optimizations: i) in blue, filter technique with a starting point corresponding to air everywhere; ii) in red, filter technique with a starting from a point corresponding to full material distribution; iii) in black, without using the filter technique with a starting point corresponding to full material distribution.</caption><thead><tr><th rowspan="3">n</th><th rowspan="3">#its (=[k])</th><th colspan="3">CPU-time</th><th colspan="3">% intermediate densities<br/>ν<sub>r</sub> ∈ ]p<sub>1</sub><sup>min</sup>+5%,p<sub>1</sub><sup>max</sup> - 5%[</th><th colspan="3">F<sub>o</sub>(ρ)<br/>(×10<sup>-2</sup>)</th></tr><tr><th rowspan="2">Filt.<br/>air</th><th rowspan="2">Filt.<br/>full</th><th rowspan="2">No Filt.<br/>full</th><th rowspan="2">Filt.<br/>air</th><th rowspan="2">Filt.<br/>full</th><th rowspan="2">No Filt.<br/>full</th><th rowspan="2">Filt.<br/>air</th><th rowspan="2">Filt.<br/>full</th><th rowspan="2">No Filt.<br/>full</th></tr><tr></tr></thead><tbody><tr><td>1</td><td>400</td><td>6h06</td><td>9h06</td><td>6h13</td><td>48.3%</td><td>32.9%</td><td>65.3%</td><td>19.6</td><td>4.59</td><td>4.37</td></tr><tr><td>2</td><td>200</td><td>3h03</td><td>4h06</td><td>3h00</td><td>30.3%</td><td>15.1%</td><td>18.6%</td><td>2.22</td><td>3.78</td><td>2.92</td></tr><tr><td>3</td><td>200</td><td>3h03</td><td>4h07</td><td>3h00</td><td>14.7%</td><td>10.3%</td><td>31.4%</td><td>1.62</td><td>2.81</td><td>2.34</td></tr><tr><td>4</td><td>200</td><td>3h03</td><td>3h35</td><td>3h00</td><td>7.34%</td><td>7.69%</td><td>7.80%</td><td>1.47</td><td>2.39</td><td>2.08</td></tr><tr><td>10</td><td>200</td><td>3h05</td><td>3h04</td><td>3h00</td><td>3.68%</td><td>3.34%</td><td>1.13%</td><td>1.36</td><td>2.08</td><td>2.10</td></tr><tr><td>20</td><td>200</td><td>3h05</td><td>3h04</td><td>3h00</td><td>4.63%</td><td>6.20%</td><td>1.01%</td><td>1.32</td><td>1.85</td><td>2.05</td></tr><tr><td>100</td><td>200</td><td>3h05</td><td>3h04</td><td>3h05</td><td>10.8%</td><td>12.4%</td><td>3.85%</td><td>1.28</td><td>1.58</td><td>2.19</td></tr></tbody></table>

colors between blue (ferromagnetic material) and red (air) in Fig. 3 are in the border between the ferromagnetic material and air in the variable zone $\Omega_V$. These difficulties are intrinsic in the use of SIMP approaches even if the penalization degree is high.

<figure><img src="image_5.png" /><figcaption>Figure 3: Topological optimization solution before post processing. (Source: figure created by authors)</figcaption></figure>

**Remark 3.** In order to obtain a circuit closer to be manufacturable, we could enforce the intermediate values of the reluctivity to be ferromagnetic material or air as follows: for $\nu_r > \frac{\nu_r^{\max}+\nu_r^{\min}}{2}$ we take $\nu_r := \nu_r^{\max}$, and for $\nu_r < \frac{\nu_r^{\max}+\nu_r^{\min}}{2}$ we take $\nu_r := \nu_r^{\min}$.

Finally in Fig. 4, in $\Omega_T$ the reference magnetic flux density $\mathbf{B}_0$ and the magnetic flux density $\mathbf{B}$, computed from the best optimized design, are plotted. Note that, $\mathbf{B}_1, \mathbf{B}_2$ and $\mathbf{B}_3$ correspond respectively to the magnetic flux density $\mathbf{B}$ at the three bars of the target area $\Omega_T$ from the bottom to the top, see Fig. 2.
## 5.4 Design problem with a large number of variables

In the above subsection, we have shown that polynomial penalization scheme improves the efficiency of our topology optimization code in 3D. Hence, the polynomial penalization is used in the following.

In this subsection, the goal is to solve problem (P) on the same magnetic circuit described in subsection 5.1, but by considering a large number of design variables. In order to achieve that, the tetrahedrons of the finite element mesh in the variable zone $\Omega_V$ are taken as cell variables of the optimization design. This yields 18893 design variables.

In this design example, the target magnetic flux density $\mathbf{B}_0$ is a field carried only along the $z$ direction with a Gaussian shape according to the coordinate $x$:

$$ \mathbf{B}_0 = \begin{pmatrix} 0 \\ 0 \\ B_{0z}(x) \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ Ae^{-\frac{(x-x_0)^2}{2\sigma^2}} \end{pmatrix}, $$

where $\sigma = 10^{-2}$ is the standard deviation of the Gaussian, $A = 4.4 \times 10^{-2}$ is a constant equal to the peak value of the Gaussian, and $x_0$ is the place in space where the peak is reached. The target zone $\Omega_T$ in this case is the black area in the Fig.10.

<figure><img src="image_12.png" /><figcaption>Figure 10: The new 3D magnetic circuit</figcaption></figure>

The choice of the penalty parameter $n$ is very difficult because, if $n$ is very large at the beginning, this leads to converge to a local minimum which will be not so interesting. Thus, we carry out a continuation process consisting in progressively increasing the penalization degree from $n = 1$ to $n = 100$. In this work, we take $n \in \{1, 2, 3, 4, 10, 20, 100\}$. Note that, for a fixed value of $n$ the optimization module stops either if the maximum number of iterations 600 is reached or if the stopping criteria is $\|\rho^{[k+1]} - \rho^{[k]}\|_{\infty} \le 10^{-4}$ is satisfied.

Increasing the penalty degree $n$, makes it possible to enforce our method to converge slowly but with more trust to a very efficient (local) optimal design which is close to be manufacturable. In Fig. 11, we plot all the iterations of our topology optimization method. The small peaks in Fig. 11, correspond to the increase of $n$. Note that the value of $F_o$ does not decrease a lot after 300 iterations but the solution becomes more and more discrete by increasing $n$; hence, the final number of iterations is 3509 corresponding to a CPU-time about 6 days on a PC-laptop with 8GB and with an AMD-Ryzen-5 processor of 2.1GHz.

Thus, the optimal design of the magnetic circuit is presented in Fig. 12. Note that this solution is not close to be entirely discrete and several variables do not respect a tolerance about 5% of $p_1^{\min}$ and $p_1^{\max}$. Indeed, 1167 variables from 18893 are in $]p_1^{\min}+5\%, p_1^{\max}-5\%[$ which represents 6.17% of all the cell variables. Note that those intermediate cell variables which have different colors between blue (ferromagnetic material) and red (air) in Fig. 12 are in the border between the
<figure><img src="image_6.png" /><figcaption>Figure 4: Comparison in $\Omega_T$ between the components of the imposed flux density $\mathbf{B}_0$ and the calculated flux density $\mathbf{B}$ in the three bars of $\Omega_T$ yielding the three density flux $\mathbf{B}_1$, $\mathbf{B}_2$ and $\mathbf{B}_3$, and a carthography of the flux density $\mathbf{B}$ of the solution design. (Source: figure created by authors)</figcaption></figure>

In Fig. 4, it is shown that the components $x$ and $y$ of the flux density $\mathbf{B}$ are almost null (less than $10^{-5}$) as asked. Moreover, note that the component $z$ of $\mathbf{B}$ in the three bars are very close to the desired Gaussian function represented by the blue curve, see Fig. 4. These results shows that we are very close to obtain the asked flux density $\mathbf{B}_0$ and thus, that we are also close to a global optimal design solution with our best optimized design obtained so far.

# 6 Conclusion

In this paper, we have extended 2D-topology optimization methods in magnetostatics based on continuous adjoint and SIMP approaches to the 3D case. Hence, some new theoretical results concerning efficient ways to compute derivative based on the adjoint method are provided for magnetostatic design problems in 3D. From this theoretical results, a quite new optimization algorithm is developed based on homotopy approaches to penalize more and more the continuous relaxation of the original discrete problem in order to enforce the convergence to a discrete solution (air or ferromagnetic material).

In a last part of the paper, numerical validations show that the gap between the gradients computed with adjoint method or with the finite difference technique is less than 0.03% for our design 3D problem in magnetostatics. Finally, our topology optimisation solver is efficiently used to solve a 3D-magnetic circuit design problem which owns up to 100,000 cell-variables.

It should be noticed that in 3D, simulation times and numerical efforts are much greater than for 2D resolutions. This limits the possibilities for numerical explorations in 3D. However, some extensions of this work could focus on developing this method to take into account the non-linearity of materials, the consideration of several physics (mechanical and magnetostatic) and other variable regions such as magnets or coils. Thus, the theory presented in this paper could be directly used to make a code capable of solving 3D multi-material design problems. Based on theoretical results in 2D, the non-linearity of materials and the consideration of mechanical constraints could be envisaged.
<figure><img src="image_13.png" /><figcaption>Figure 11: Evolution of the value of the objective function during the iterations</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 12: Topological optimization solution before post processing</figcaption></figure>
References

Aage, N., Andreassen, E., Lazarov, B. S. & Sigmund, O. (2017), ‘Giga-voxel computational
morphogenesis for structural design’, *Nature* **550**(7674), 84–86.

Allaire, G. (2007), *Conception optimale de structures*, Springer, Berlin.

Allaire, G. & Alouges, F. (2015), *Polycopié du cours MAP 431: Analyse variationnelles des équations aux dérivées partielles*, École Polytechnique, Paris.

Bendsøe, M. P. & Sigmund, O. (1999), ‘Material interpolation schemes in topology optimization’, *Archive of applied mechanics* **69**, 635–654.

Bendsøe, M. P. & Sigmund, O. (2003), *Topology optimization: theory, methods, and applications*, Springer, Berlin.

Cherrière, T., Laurent, L., Hlioui, S., Louf, F., Duysinx, P., Geuzaine, C., Ben Ahmed, H., Gabsi, M. & Fernández, E. (2022), ‘Multi-material topology optimization using wachspress interpolations for designing a 3-phase electrical machine stator’, *Structural and Multidisciplinary Optimization* **65**(12), 352.

Dular, P. & Geuzaine, C. (n.d.), GetDP reference manual.
URL: https://getdp.info/doc/texinfo/getdp.pdf

Gangl, P. & Sturm, K. (2021), ‘Asymptotic analysis and topological derivative for 3d quasi-linear magnetostatics’, *ESAIM: Mathematical Modelling and Numerical Analysis* **55**, S853–S875.

Gauthey, T., Gangl, P. & Hassan, M. H. (2021), ‘Multi-material topology optimization with continuous magnetization direction for permanent magnet synchronous reluctance motors’, *arXiv preprint arXiv:2107.04825*.

Gauthey, T., Hassan, M. H., Messine, F. & Gillon, F. (2023), ‘Topology optimization of the harmonic content for torque ripple minimization’, *IEEE Transactions on Magnetics* pp. 1–4.

Geuzaine, C. & Remacle, J.-F. (2009), ‘Gmsh: A 3-d finite element mesh generator with built-in pre-and post-processing facilities’, *International journal for numerical methods in engineering* **79**(11), 1309–1331.

Lucchini, F., Torchio, R., Cirimele, V., Alotto, P. & Paolo, B. (2022), ‘Topology optimization for electromagnetics: A survey’, *IEEE Access* **10**, 98593–98611.

Masmoudi, M., Pommier, J. & Samet, B. (2005), ‘The topological asymptotic expansion for the maxwell equations and some applications’, *Inverse Problems* **21**(2), 547.

Remacle, J.-F. & Geuzaine, C. (n.d.), Gmsh reference manual.
URL: https://gmsh.info/doc/texinfo/gmsh.pdf

Sanogo, S. (2016), Conception optimale de circuits magnétiques dédiés à la propulsion spatiale électrique par des méthodes d'optimisation topologique, PhD thesis, Université Paul Sabatier-Toulouse III.

Sanogo, S. & Messine, F. (2018), ‘Topology optimization in electromagnetism using simp method: Issues of material interpolation schemes’, *COMPEL-The international journal for computation and mathematics in electrical and electronic engineering* **37**(6), 2138–2157.

Sanogo, S., Messine, F., Henaux, C. & Vilamot, R. (2014), ‘Topology optimization for magnetic circuits dedicated to electric propulsion’, *IEEE transactions on magnetics* **50**(12), 1–13.

Seebacher, P., Kaltenbacher, M., Wein, F. & Lehmann, H. (2021), ‘A pseudo density topology optimization approach in nonlinear electromagnetism applied to a 3d actuator’, *International Journal of Applied Electromagnetics and Mechanics* **65**(3), 545–559.

Sigmund, O. (2007), ‘Morphology-based black and white filters for topology optimization’, *Structural and Multidisciplinary Optimization* **33**, 401–424.

Sigmund, O. & Maute, K. (2013), ‘Topology optimization approaches: A comparative review’, *Structural and Multidisciplinary Optimization* **48**(6), 1031–1055.
ferromagnetic material and air in the variable zone $\Omega_V$. These difficulties are intrinsic in the use of SIMP approaches even if the penalization degree is high.

Hence, to construct a manufacturable circuit, we enforce the intermediate values of the reluc-
tivity to be ferromagnetic material or air as follows: for $\nu_r > \frac{\nu_r^{\max}+\nu_r^{\min}}{2}$ we take $\nu_r := \nu_r^{\max}$, and
for $\nu_r < \frac{\nu_r^{\max}+\nu_r^{\min}}{2}$ we take $\nu_r := \nu_r^{\min}$.

In Fig. 13, we display the manufacturable optimal design circuit and in Fig. 14, the manufac-
turable design solution is also plotted in 3 distinct 3D-views.

<figure><img src="image_15.png" /><figcaption>Figure 13: Topological optimization solution after post processing</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 14: The manufacturable optimal design in several views</figcaption></figure>

Finally in Fig. 15, the difference in $\Omega_T$ between $\mathbf{B}_0$ and the magnetic flux density for the optimal manufacturable design is plotted. $\mathbf{B}_1$, $\mathbf{B}_2$ and $\mathbf{B}_3$ correspond respectively to the magnetic flux density at the three bars of the target area $\Omega_T$ from bottom to top, see Fig. 10.
Svanberg, K. (1987), ‘The method of moving asymptotes—a new method for structural optimization’, *International journal for numerical methods in engineering* **24**(2), 359–373.

Youness, R. & Messine, F. (2019a), Adjoint based topology optimization in nonlinear magnetostatics application to hall effect thrusters, in ‘2019 22nd International Conference on the Computation of Electromagnetic Fields (COMPUMAG)’, IEEE, pp. 1–4.

Youness, R. & Messine, F. (2019b), ‘An implementation of adjoint-based topology optimization in magnetostatics: Application to design hall-effect thrusters’, *COMPEL-The international journal for computation and mathematics in electrical and electronic engineering* **38**(3), 1023–1035.
<figure><img src="image_17.png" /><figcaption>Figure 15: Comparison between the components of the imposed field **B**₀ and the calculated field **B**</figcaption></figure>

Notice that the $B_x$ and $B_y$ components of the recovered field are of the order of $10^{-5}$ while the Gaussians of the three lines are very close to the desired Gaussian.

# 6 Conclusion

In this paper, we have extended 2D-topology optimization methods in electromagnetism based on adjoint and SIMP approaches to the 3D case. Hence, some new theoretical results concerning efficient ways to compute derivative based on adjoint method are provided. Then, a quite new optimization algorithm is developed based on homothopy approaches to penalize more and more the continuous relaxation of the original discrete problem in order to enforce the convergence to discrete solution (air or ferromagnetic material).

In a last part of the paper, numerical validations show that the gap between the gradients computed with adjoint method or with the finite difference technique is less than 0.03% for our design 3D problem in electromagnetism. Finally, our topology optimisation solver is efficiently used to solve two 3D-magnetic circuit design problems where the second design problems owns 18893 cell-variables.

# Acknowledgements

This work was supported by grants from the French National Research Agency (ANR-21-CE05-0001-01).

# Conflict of interest

On behalf of all authors, the corresponding author states that there is no conflict of interest.
Replication of results

Our work is based on the Matlab R2020b programming language and on the GMSH/GetDP software for solving partial differential equations by the finite element method.

All information required to replicate the results of this work is clearly detailed in the manuscript.
The models and parameters of the considered design problems are presented in Section 5. More-
over, the optimization algorithm used to obtain the results presented in Section 5, is provided
in Section 4. Interested readers may send a request to the corresponding author to obtain the
MatLab files of codes.

References

[1] Niels Aage, Erik Andreassen, Boyan S Lazarov, and Ole Sigmund. Giga-voxel computational morphogenesis for structural design. *Nature*, 550(7674):84–86, 2017.

[2] Grégoire Allaire. *Conception optimale de structures*. Springer, Berlin, 2007.

[3] Grégoire Allaire and François Alouges. *Polycopié du cours MAP 431: Analyse variationnelles des équations aux dérivées partielles*. École Polytechnique, Paris, 2015.

[4] Martin P Bendsøe and Ole Sigmund. Material interpolation schemes in topology optimization. *Archive of applied mechanics*, 69:635–654, 1999.

[5] Martin Philip Bendsøe and Ole Sigmund. *Topology optimization: theory, methods, and applications*. Springer, Berlin, 2003.

[6] Théodore Cherrière, Luc Laurent, Sami Hlioui, François Louf, Pierre Duysinx, Christophe Geuzaine, Hamid Ben Ahmed, Mohamed Gabsi, and Eduardo Fernández. Multi-material topology optimization using wachspress interpolations for designing a 3-phase electrical machine stator. *Structural and Multidisciplinary Optimization*, 65(12):352, 2022.

[7] Patrick Dular and Christophe Geuzaine. *GetDP reference manual*.

[8] Pauline Ferrouillat. *Développement de formulations éléments finis 3D en potentiel vecteur magnétique: application aux machines asynchrones en mouvement*. PhD thesis, Université Grenoble Alpes (ComUE), 2015.

[9] Thomas Gauthey, Peter Gangl, and Maya Hage Hassan. Multi-material topology optimization with continuous magnetization direction for permanent magnet synchronous reluctance motors. *arXiv preprint arXiv:2107.04825*, 2021.

[10] Christophe Geuzaine and Jean-François Remacle. Gmsh: A 3-d finite element mesh generator with built-in pre-and post-processing facilities. *International journal for numerical methods in engineering*, 79(11):1309–1331, 2009.

[11] Dalibor Lukas. An integration of optimal topology and shape design for magnetostatics. *Scientific Computing in Electrical Engineering*, 9:227–232, 2007.

[12] Jean-François Remacle and Christophe Geuzaine. *Gmsh reference manual*.

[13] George IN Rozvany. Aims, scope, methods, history and unified terminology of computer-aided topology optimization in structural mechanics. *Structural and Multidisciplinary optimization*, 21:90–108, 2001.

[14] Satafa Sanogo. *Conception optimale de circuits magnétiques dédiés à la propulsion spatiale électrique par des méthodes d'optimisation topologique*. PhD thesis, Université Paul Sabatier-Toulouse III, 2016.

[15] Satafa Sanogo and Frédéric Messine. Topology optimization in electromagnetism using simp method: Issues of material interpolation schemes. *COMPEL-The international journal for computation and mathematics in electrical and electronic engineering*, 37(6):2138–2157, 2018.
[16] Satafa Sanogo, Frederic Messine, Carole Henaux, and Raphael Vilamot. Topology optimization for magnetic circuits dedicated to electric propulsion. *IEEE transactions on magnetics*, 50(12):1–13, 2014.

[17] Ole Sigmund. Morphology-based black and white filters for topology optimization. *Structural and Multidisciplinary Optimization*, 33:401–424, 2007.

[18] Ole Sigmund and Kurt Maute. Topology optimization approaches: A comparative review. *Structural and Multidisciplinary Optimization*, 48(6):1031–1055, 2013.

[19] Krister Svanberg. The method of moving asymptotes—a new method for structural optimization. *International journal for numerical methods in engineering*, 24(2):359–373, 1987.

[20] Rtimi Youness and Frederic Messine. Adjoint based topology optimization in nonlinear magnetostatics application to hall effect thrusters. In *2019 22nd International Conference on the Computation of Electromagnetic Fields (COMPUMAG)*, pages 1–4. IEEE, 2019.

[21] Rtimi Youness and Frederic Messine. An implementation of adjoint-based topology optimization in magnetostatics: Application to design hall-effect thrusters. *COMPEL-The international journal for computation and mathematics in electrical and electronic engineering*, 38(3):1023–1035, 2019.