Level set topology optimization of metamaterial-based heat manipulators using isogeometric analysis

label1]Chintan Jansari
label1,label3]Stéphane P.A<EMAIL_ADDRESS>label2]Elena Atroshchenko

[cor1]Corresponding author
[label1]Institute of Computational Engineering, Faculty of Sciences, Technology and Medicine, University of Luxembourg, Luxembourg City, Luxembourg.
[label3]Clyde Visiting Fellow, Department of Mechanical Engineering, The University of Utah, Salt Lake City, Utah, United States.
[label2]School of Civil and Environmental Engineering, University of New South Wales, Sydney, Australia.

We exploit level set topology optimization to find the optimal material distribution for metamaterial-based heat manipulators. The level set function, geometry, and solution field are parameterized using the non-uniform rational B-spline (NURBS) basis functions in order to take advantage of easy control of smoothness and continuity. In addition, NURBS approximations can produce conic geometries exactly and provide higher efficiency for higher-order elements. The values of the level set function at the control points (called expansion coefficients) are utilized as design variables. For optimization, we use an advanced mathematical programming technique, Sequential Quadratic Programming (SQP). Taking into account a large number of design variables and the small number of constraints associated with our optimization problem, the adjoint method is utilized to calculate the required sensitivities with respect to the design variables. The efficiency and robustness of the proposed method are demonstrated by solving three numerical examples. We have also shown that the current method can handle different geometries and types of objective functions. In addition, regularization techniques such as Tikhonov regularization and volume regularization have been explored to reduce unnecessary complexity and increase the manufacturability of optimized topologies.

Level set topology optimization Thermal cloak Thermal camouflage Thermal metamaterials Adjoint method Isogeometric analysis.

§ INTRODUCTION

Due to the special arrangement of the constituent materials, thermal metamaterials can have heat transfer capabilities superior to those of materials available in nature. Therefore, researchers have proposed the use of artificially created thermal metamaterials to control heat fluxes. The concept opened research opportunities in heat transfer applications. Creating devices that are thermal equivalent to resistors, capacitors, inductors, diodes, transistors, etc. is one such opportunity. These devices are called heat manipulators, and several of them have already been proposed, such as thermal cloak <cit.>, thermal concentrator <cit.>, thermal camouflage <cit.>, heat flux inverter <cit.> etc. The spatial layout of the member materials of a metamaterial-based heat manipulator has a significant impact on its performance. Because of this, structural optimization can be a useful technique for understanding the impact of the spatial arrangement of member materials and for developing superior designs. To the knowledge of the authors, work on optimization of the metamaterial-based heat manipulator is limited <cit.>. The authors already published an article on isogeometric shape optimization in combination with the gradient-free Particle Swarm Optimization (PSO) algorithm <cit.>. The present work can be considered as an extension of the previous work, as in this work, we are using the topology optimization to avoid the limitations imposed by shape optimization and explore even larger design space.
The goal of structural optimization is to determine the material distribution in the design domain that gives the best-desired performance (of a device or structure). In the early days of its inception, structural optimization primarily focused on the optimization of mechanical structures. Over a period of time, optimization transpired as a more general numerical technique that can cover a wider range of physical problems including fluids, optics, acoustics, thermodynamics, etc. Lately, topology optimization has become increasingly popular compared to other structural optimization techniques such as shape and size optimization. The reason for its popularity is its ability to allow changes in topology during optimization, hence avoiding the need for a close-to-optimal initial design. Several topology optimization approaches have been developed to date. The most common ones are density-based, level set-based, phase field-based, and evolutionary algorithm-based. Despite the fact that these approaches define diverse directions for topology optimizations, there are few conceptual differences between them. It is difficult to determine which technique is best for a given problem due to the lack of direct-relevant comparisons <cit.>. In addition, not simply the choice of a method but also factors like optimizers, filters, constraints, etc., define the overall performance. However, particularly in the level set based approach, the interface is clearly defined through out the optimization process, and therefore, we focus on level set topology optimization considering their possible advantage in applying objective function or constraints on the interface in future work <cit.>. The remaining approaches are covered in detail in review articles <cit.>.
The level set method (LSM), proposed by Osher and Sethian <cit.>, captures moving interfaces in multi-phase flow. The idea to incorporate LSM with topology optimization was suggested by Haber and Bendsoe <cit.>. Following this, several research groups began working and subsequently published the level set topology optimization method <cit.>. In Osher and Santosa <cit.>, Allaire et al. <cit.>, and Wang et al. <cit.>, the shape-sensitivity-based level set topology optimization framework was introduced. Specifically for heat manipulators, Fujii and Akimoto have explored the level set topology optimization method. They used the method to optimize a thermal cloak <cit.>, a thermal cloak-concentrator <cit.>, as well as a combined thermal-electric cloak <cit.>. In their work, the finite element method, for the boundary value problem, and a stochastic evolution strategy, for optimization, are utilized. Conversely, we use the isogeometric analysis for the boundary value problem and a gradient-based Sequential Quadratic Programming (SQP) method for optimization. In the next few paragraphs, we discuss the level set topology optimization method in addition to the other aspects of the optimization process.
In level set topology optimization, the isocontours of a level set function (LSF) implicitly define the interfaces between the material phases. Accordingly, the topology changes with the motion of these isocontours. Additionally, the mesh remains unchanged, eliminating the cost associated with creating a new mesh at each iteration.

The components of the level set topology optimization.

There are three main components of the level set topology optimization: the level set function (LSF) parameterization (design space), an efficient numerical method to solve the boundary value problem, and the optimization strategy. In this work, NURBS basis functions are employed to parameterize the LSF, geometry, and solution field (i.e., the temperature distribution). Isogeometric analysis (IGA) <cit.> is utilized to solve the thermal boundary value problem. Here, we decouple LSF parameterization from geometry and solution field by taking two different NURBS bases. However, it will be advantageous to decouple all three parameterizations using tools such as Geometry Independent Field approximaTion (GIFT) <cit.>. For optimization, we use a mathematical programming approach-Sequential Quadratic Programming (SQP) <cit.>. The motivation behind the above-mentioned choices and their respective alternatives are discussed in the following paragraphs. fig:method overview shows the overview of the method with its components and particular choices for each component.
At first, we discuss LSF parameterization. The LSF is approximated over the whole design domain from a few point values using LSF parameterization, and these values at nodes/control points work as design variables. The LSF parameterization decides the design freedom, level of detail, and nature of the optimization problem. Hence, the effort required for optimization strongly depends on the LSF parameterization. In the present work, the parameterization for the LSF is decoupled from the geometry and solution field parameterizations. Decoupling allows for securing the required accuracy of the solution without increasing the complexity of the optimization. This idea is standard in LSM for free boundary problems, as discussed in <cit.>.
The use of NURBS parameterizations and IGA provides several advantages over Lagrange parameterizations and the conventional finite element method (FEM) <cit.>: (i) easy control of smoothness and inter-element continuity, (ii) an exact representation of conic geometries, and (iii) higher efficiency for higher-order elements. In <cit.>, Hughes et al. proposed Isogeometric Analysis (IGA) based on the discretization of a Galerkin formulation using NURBS basis functions.
Later, several other variants of IGA such as isogeomteric collocation method (IGA-C) <cit.>, isogeometric boundary element method (IGABEM) <cit.>, geometry independent field approximation (GIFT) <cit.> have been proposed. IGA and its variants have been successfully implemented in shape and topology optimization framework <cit.>. In our work, we exploit the parameterized level set topology optimization incorporating IGA proposed by Wang et al. <cit.>. However, instead of the immersed boundary technique to map the geometry to a numerical model as in <cit.>, we use density-based point geometry mapping due to its relatively easy implementation. Also, density-based mapping avoids ill-conditioning issues which are common in immersed boundary techniques. In density-based mapping, no special treatment is performed for the integration. Thermal conductivity at an integration point for numerical analysis is defined directly from the LSF value. However, the material definition by this simplistic approach is not very accurate near the interface. Therefore, by using a fine enough solution mesh, it ensured that the density-based mapping does not deteriorate the solution accuracy.
Lastly, the optimization procedure is discussed. It includes two parts: (i) update information and (ii) update procedure.  Update information often consists of the objective function and constraint sensitivities with respect to design variables. The most common methods for sensitivity analysis are (a) the direct method, (b) the finite difference method, (c) the semi-analytical method, and (d) the adjoint method <cit.>. The direct method requires to solve an extra system for each design variable; therefore, the direct method is significantly costly for complex problems. Similarly, the finite difference method requires to solve ${n+1}$ boundary value problems for $n$ design variables. Furthermore, the accuracy of the finite difference method depends on the perturbation size. The semi-analytical approach is computationally efficient; however, its accuracy can be relatively unsatisfactory for special cases due to the incompatibility of the design sensitivity field with the structure <cit.>. The adjoint method requires to solve ${n+1}$ system for $n$ objectives and constraints that depend on the field solution. The adjoint system is efficient for a problem with a large number of design variables and a small number of constraints, which aligns with our case. Therefore, in this work, the adjoint method <cit.> is used to calculate the sensitivity.
The update procedure decides how to use the update information to advance the level set interfaces. For the update procedure, two classes of methods are available in the literature (a) Hamilton-Jacobi (HJ) equation-based procedures and (b) mathematical programming. In the first class, the problem is considered a quasi-temporal problem. The interface motion is calculated based on the solution of the Hamilton-Jacobi (HJ) equation in pseudo-time <cit.>. The second class, mathematical programming <cit.> often equipped with sophisticated step selection, constraint handling strategies, as well as optimized speed and efficiency. In this work, Sequential Quadratic Programming (SQP) <cit.> is chosen due to its ability to accurately solve nonlinear constrained problems. In addition, MATLAB has an inbuilt subroutine on SQP in its `fmincon' optimization tool, which provides an advantage from an implementation point of view.
The remainder of the paper is organized as follows: sec:Level set topology optimization provides the numerical formulation of the level set topology optimization that includes boundary value problem formulation, implicit interface representation, and numerical approximations. The optimization problem, the sensitivity analysis, the SQP algorithm, and the regularization techniques are explained in sec:Optimization problem. In sec:Numerical Examples, three numerical examples are demonstrated: a toy problem of an annular ring with the known solutions for the state and adjoint BVPs (sec:BenchAnncase), the thermal cloak problem (sec:Chen2015case_cloak), and the thermal camouflage problem (sec:Yang2016case_cmflg), which corroborates the efficiency and robustness of the proposed method. sec:Conclusions presents the main conclusions of the current work.

§ LEVEL SET TOPOLOGY OPTIMIZATION WITH ISOGEOMETRIC ANALYSIS

§.§ Boundary value problem description

As shown in fig:BVP domain, a heat manipulator embedded in the domain $\mathrm{\Omega} \in\mathbb{R}^{2}$ is considered. Note that all formulations given here stand for three dimensional physical space. The domain is externally bounded by $\mathrm{\Gamma}=\partial\mathrm{\Omega}=\mathrm{\Gamma}_D \cup \mathrm{\Gamma}_N$, where $\mathrm{\Gamma}_D$ and $\mathrm{\Gamma}_N$  are two parts of the boundary, where the Dirichlet and Neumann boundary conditions are applied, respectively. $\mathrm{\Gamma}_D \cap \mathrm{\Gamma}_N =\emptyset$. The embedded heat manipulator uniquely divides the domain into 3 different parts; the inside region $\mathrm{\Omega}_{\rm in}$, the heat manipulator region $\mathrm{\Omega}_{\rm design}$, and the outside region $\mathrm{\Omega}_{\rm out}$. $\mathrm{\Omega}=\mathrm{\Omega}_{\rm in} \cup \mathrm{\Omega}_{\rm design} \cup \mathrm{\Omega}_{\rm out}$. The internal boundaries between these parts are collectively denoted by $\mathrm{\Gamma}_I=\mathrm{\Gamma}_{I_{\rm in}} \cup \mathrm{\Gamma}_{I_{\rm out}}$.
In this work, we restrict our scope to metamaterials made of two isotropic member materials. Therefore, in addition to above-mentioned partition, $\mathrm{\Omega}_{\rm design}$ is divided into two parts $\mathrm{\Omega}_{l_1}$ and $\mathrm{\Omega}_{l_2}$, representing two member materials ($\mathrm{\Omega}_{\rm design}=\mathrm{\Omega}_{\ell_1}\cup\mathrm{\Omega}_{\ell_2}$). An interface $\mathrm{\Gamma}_L$ between the member materials is implicitly defined by the level set function. The description of the level set function and corresponding interface $\mathrm{\Gamma}_L$ will be explained in detail in sec:Implicit bound with LSF. We simplify the problem with the following assumptions: the temperature and normal flux are continuous along $\mathrm{\Gamma}_I$, the heat conduction is the only present form of heat transfer, and there is no internal heat generation. The steady-state thermal boundary value problem for the given temperature field $T$ can be written as,

\begin{align}
\nabla \cdot \left( \boldsymbol{\kappa}\nabla T\right) &= 0 \quad &&\textrm{in} \quad \mathrm{\Omega}, \label{eq:Laplace equation}\\
T &= T_D \quad &&\textrm{on} \quad \mathrm{\Gamma}_D, \label{eq:Boundary conditions a}\\
\nabla T\cdot \boldsymbol{n} &= Q_N \quad &&\textrm{on} \quad \mathrm{\Gamma}_N, \label{eq:Boundary conditions b}\\
\left\llbracket T\right\rrbracket &=0 \quad &&\textrm{on} \quad \mathrm{\Gamma}_I, \label{eq:Interface conditions a}\\
\boldsymbol{n}\cdot \left\llbracket \boldsymbol{\kappa}\nabla T\right\rrbracket &= 0 \quad &&\textrm{on} \quad \mathrm{\Gamma}_I, \label{eq:Interface conditions b}
\end{align}

where $\boldsymbol{\kappa}$ is the thermal conductivity matrix (for isotropic material, $\boldsymbol{\kappa}=\kappa \mathrm{\mathbf{I}_2}$ with $\mathrm{\mathbf{I}_2}$ be an identity matrix of $\mathbb{R}^{2}$), $Q_N$ is the flux applied on $\mathrm{\Gamma}_N$, $T_D$ is the prescribed temperature on $\mathrm{\Gamma}_D$,  $\boldsymbol{n}$ is the unit normal on the boundary, $\llbracket\cdot\rrbracket$ is the jump operator, and $\nabla= \small{\left({ \dfrac{\textstyle\partial}{\textstyle\partial x},\dfrac{\textstyle\partial}{\textstyle\partial y}}\right)}$. On the internal boundary $\mathrm{\Gamma}_I$, $\boldsymbol{n}=\boldsymbol{n}^1= -\boldsymbol{n}^2$  where the connected patches at $\mathrm{\Gamma}_I$ are denoted by 1 and 2.

Domain description of the boundary value problem. $\Omega_{\rm design}$ represents the region of the heat manipulator, $\Omega_{\rm in}~\&~\Omega_{\rm out}$ are, respectively, the inside and outside regions with respect to $\Omega_{\rm design}$. $\Omega$=$\Omega_{\rm in}\cup\Omega_{\rm out}\cup\Omega_{\rm design}$ . The solid black line shows an explicitly defined interface $\Gamma_I$, while the broken black & white line shows an implicitly defined interface $\Gamma_L$ (that separates $\Omega_{\rm design}$ in two parts $\Omega_{\ell_1}$ and $\Omega_{\ell_2}$). The detailed view highlights the matching of control points of connecting patches at the interface $\Gamma_I$.

To solve the boundary value problem, the strong form described in eq:Heat conduction BVP is transformed into the weak form using the standard Bubnov-Galerkin formulation. The weak formulation is given as follows: Find $
T^h \in \mathscr{T}^h \subseteq \mathscr{T} = \big\lbrace T \in \mathbb{H}^1 ({\mathrm{\Omega}}),   T= T_D \hspace{0.15cm} \textrm{on} \ {\mathrm{\Gamma}}_D \big\rbrace
$ such that $\forall S^h \in \mathscr{S}^h_0 \subseteq \mathscr{S}_0 = \left\lbrace S \in \mathbb{H}^1 ({\mathrm{\Omega}}), S=0 \hspace{0.15cm} \textrm{on} \ {\mathrm{\Gamma}}_D \right\rbrace$,

\begin{equation}
a(T^h,S^h) = \ell(S^h),
\label{eq:weak_form}
\end{equation}

\begin{equation}
a(T^h,S^h) = \int _{\mathrm{\Omega}} (\nabla  S^h)^{\rm T} \boldsymbol{\kappa}\nabla T^h d\mathrm{\Omega},
\label{eq:weak_form_a}
\end{equation}

\begin{equation}
\ell{(S^h)} = \int _{\mathrm{\Gamma}_{N}} (S^h)^{\rm T} Q_N  d\mathrm{\Gamma}.
\label{eq:weak_form_l}
\end{equation}

The interface continuity conditions described in eq:Interface conditions a-(<ref>) are applied using Nitsche's method <cit.>. Nitsche's method applies the interface conditions weakly while preserving the coercivity and consistency of the bilinear form. Nitsche's method lies between the Lagrange multiplier method and the penalty method, designed to overcome some of the limitations of these conventional methods such as over-sensitivity to the penalty parameter, inconsistency of variational form, limitations imposed by stability conditions etc. Nitsche's method modifies the bilinear form by substituting the Lagrange multipliers of the Lagrange method by their actual physical representation, i.e. normal flux.  As shown in <cit.>, the bilinear form after modification appears as follows,

\begin{multline}
a(T^h,S^h) = \int _{\mathrm{\Omega}} (\nabla  S^h)^{\rm T} \boldsymbol{\kappa}\nabla T^h d\mathrm{\Omega} - \int _{\mathrm{\Gamma}_I} \left(\boldsymbol{n}\cdot \{\boldsymbol{\kappa}\nabla S^h\}\right)^{\rm T}\llbracket T^h \rrbracket~d\mathrm{\Gamma} \\- \int _{\mathrm{\Gamma}_I} \llbracket S^h  \rrbracket^{\rm T} \left(\boldsymbol{n}\cdot\{\boldsymbol{\kappa}\nabla T^h\}\right)  d\mathrm{\Gamma} + \int _{\mathrm{\Gamma}_I} \beta~\llbracket S^h  \rrbracket^{\rm T} \llbracket T^h  \rrbracket~d\mathrm{\Gamma},
\label{eq:modified_weak_form_a}%
\end{multline}

where $\beta$ is the stabilization parameter and $\{\cdot\}$ is the averaging operator defined as $\{\theta\}=\gamma\theta^1 + (1-\gamma)\theta^2$ with $\gamma$ being the averaging parameter ($0<\gamma<1$) <cit.>. For the current work, $\beta=1\times 10^{12}$ and $\gamma=0.5$. In the literature <cit.>, it is also reported that the large stabilization parameter might cause ill-conditioning of the system, but we did not face any conditioning issue for our boundary value problem.

§.§ Implicit boundary representation with level set function

As described in sec:Boundary value problem, the  interface $\Gamma_L$, inside $\Omega_{\rm design}$ that separates the member materials, is not explicitly defined in eq:Heat conduction BVP. It is defined by a level set function, and its movement is the main act of level set topology optimization.
In the level set method in $\mathbb{R}^d$, the interface between two materials is implicitly represented by an isosurface of a scalar function  $\varPhi:\mathbb{R}^{d} \rightarrow \mathbb{R}$, called level set function (LSF). Accordingly, in our work, an isosurface of LSF $\varPhi:\mathbb{R}^{2} \rightarrow \mathbb{R}$, $\varPhi=0$, defines the interface $\Gamma_{L}$ separating $\Omega_{\ell_1}$ and $\Omega_{\ell_2}$ in $\Omega_{\rm design}$ as shown in fig:levelset representation. The overall level set representation can be given as,

\begin{align}
\varPhi(\boldsymbol{x})&>0\quad \forall\boldsymbol{x} \in \Omega_{\ell_1}\setminus \Gamma_L, \\
\varPhi(\boldsymbol{x})&=0 \quad\forall\boldsymbol{x} \in \Gamma_L,\\
\varPhi(\boldsymbol{x})&<0 \quad\forall\boldsymbol{x} \in \Omega_{\ell_2}\setminus \Gamma_L,
\end{align}

During optimization, the movement of the interface is defined by the evolution of the isosurface $\varPhi=0$, while the background Eulerian mesh remains fixed.

Level set function $\varPhi$ in 3D

Level set function $\varPhi$ in 2D

Level set representation to define the material distribution. (a) The level set function $\varPhi$ in 3D. Height represents the value of $\varPhi$. $\varPhi$ above and below the plane, $\varPhi=0$, represent two member materials (shown by two different colors). (b) The level set function $\varPhi$ in 2D. The material interface is calculated by the intersection of function $\varPhi$ with the plane, $\varPhi=0$.

In the level set method, the geometric mapping defines how the LSF information is utilized in the numerical solution of the boundary value problem. Via accuracy of the numerical solution, the geometric mapping affects the optimization results. There are three common geometric mapping approaches <cit.>: (a) conformal discretization, (b) immersed boundary techniques <cit.>, and (c) density-based mapping. In conformal discretization, the mesh conforms to the interface defined by the LSF. The method is distinct from shape optimization, as the LSF governs the changes in shape. The approach provides a crisp interface representation and the most accurate solution. However, it becomes expensive due to remeshing at each iteration. On the other hand, the immersed boundary techniques allow the nonconforming mesh. In these methods, the mesh remains fixed and the interface is captured in the numerical model using special treatment. Immersed boundary techniques also have a crisp interface representation and allow the enforcement of interface conditions directly. A specialized code for numerical integration and field approximation is needed for the elements cut by level set interfaces. Sometimes immersed boundary techniques face issues of noise and ill-conditioning due to small intersections. The last approach, density-based mapping, is the most common because of its easy implementation.
For the same advantage, we utilize density-based geometric mapping. We explore point-wise density mapping. In other words, the LSF value at an integration point defines the material density at that particular point. Hence, the thermal conductivity $\boldsymbol{\kappa}$ becomes directly a function of the LSF value. It takes the following form,

\begin{equation}
\boldsymbol{\kappa} = \begin{cases}
\boldsymbol{\kappa}^1 = \kappa^1 \mathrm{\mathbf{I}_2}  \quad\text{if } \boldsymbol{x}\in\Omega_{l_1} \quad \textit{i.e.}~\varPhi \geqslant 0,\\
\boldsymbol{\kappa}^2 = \kappa^2 \mathrm{\mathbf{I}_2} \quad\text{if } \boldsymbol{x}\in\Omega_{l_2}\setminus \Gamma_L \quad \textit{i.e.}~\varPhi < 0,
\end{cases}
\end{equation}

which can be represented using the Heaviside function $H$ as,

\begin{equation}
\boldsymbol{\kappa}(\varPhi) = \boldsymbol{\kappa}^1 H(\varPhi) + \boldsymbol{\kappa}^2 (1 - H(\varPhi)),
\label{eq:kappa based on LSF}
\end{equation}

\begin{equation}
1 \quad\text{if } \varPhi \geqslant 0,\\
0 \quad\text{if } \varPhi < 0.
\end{cases}
\label{eq:Heaviside function}
\end{equation}

For the sensitivity analysis (see sec:Sensitivity analysis), we need the derivative of $\boldsymbol{\kappa}$ with respect to $\varPhi$, and hence the derivative of the Heaviside function with respect to $\varPhi$, i.e. the Dirac delta function $\delta(\varPhi)$. The singularity of the Dirac delta function brings numerical issues while calculating the derivatives; therefore, both functions are often replaced by their approximations. There are several forms of an approximate Heaviside function available in the literature <cit.>. Here, we use a polynomial form given as,

\begin{equation}\label{eq:Heavised-approx}
\alpha \quad &\text{if} \quad \varPhi < -\Delta,\\
\dfrac{3(1-\alpha)}{4} \left(\dfrac{\varPhi}{\Delta}-\dfrac{\varPhi^{3}}{3 \Delta^{3}}\right) +\dfrac{1+\alpha}{2} \quad&\text{if} -\Delta \leqslant \varPhi < \Delta,\\
1 \quad &\text{if} \quad\varPhi  \geqslant \Delta,\\
\end{cases}
\end{equation}

where $\alpha$ is a small positive value (here we take $\alpha=0$) and $\Delta$ is the support bandwidth.
Consequently, the derivative of $\boldsymbol{\kappa}$ can be written as,

\begin{equation}
\dfrac{d \boldsymbol{\kappa} (\varPhi)}{d \varPhi} = (\boldsymbol{\kappa}^1 - \boldsymbol{\kappa}^2)\delta(\varPhi),
\label{eq:kappa der based ontakes LSF}
\end{equation}

where one-dimensional Dirac delta function $\delta(\varPhi)$ is approximated by,

\begin{equation}\label{delta-approx}
\delta(\varPhi)=
\begin{cases}
\dfrac{3(1-\alpha)}{4\Delta} \left(1-\dfrac{\varPhi^{2}}{ \Delta^{2}}\right)  \quad&\text{if} \quad |\varPhi|\leqslant \Delta,\\
0 \quad &\text{if} \quad|\varPhi| > \Delta.\\
\end{cases}
\end{equation}

The smoothed Heaviside and smoothed Dirac delta functions in the given polynomial form are shown in fig:smoothed H and delta fn.

Smoothed Heaviside function and smoothed Dirac delta function.

§.§ Solution of the boundary value problem using IGA

Following the implicit representation of the interface by the level-set function, the weak form (eq:modified_weak_form_a) is modified again to accommodate the conductivity matrix functional $\boldsymbol{\kappa}(\varPhi)$. Since the LSF is defined in domain $\Omega_{\text{design}}$, which does not have a common boundary with $\partial\Omega$, the applied boundary condition is not a function of the LSF. Therefore, the linear form (eq:weak_form_l) remains unchanged. The modified weak form transforms to,

\begin{equation}\label{eq:mmodified_weak_form}
a(T^h,S^h,\varPhi) = \ell(S^h),
\end{equation}

\begin{multline}\label{eq:mmodified_weak_form_a}
a(T^h,S^h,\varPhi) = \int _{\Omega} (\nabla  S^h)^{\rm T} \boldsymbol{\kappa}(\varPhi)\nabla T^h d\Omega - \int _{\Gamma_I} \left(\boldsymbol{n}\cdot \{\boldsymbol{\kappa}(\varPhi)\nabla S^h\}\right)^{\rm T}\llbracket T^h \rrbracket~d\Gamma \\- \int _{\Gamma_I} \llbracket S^h  \rrbracket^{\rm T} \left(\boldsymbol{n}\cdot\{\boldsymbol{\kappa}(\varPhi)\nabla T^h\}\right)  d\Gamma + \int _{\Gamma_I} \beta~\llbracket S^h  \rrbracket^{\rm T} \llbracket T^h  \rrbracket~d\Gamma,
\end{multline}

Parametrization of a point from the parametric domain to a point in the physical domain using NURBS basis functions.

The next step is to define the geometry and solution approximations to discretize the weak form. As mentioned in sec:Introduction, we use NURBS basis functions for geometry and solution field approximations. Here we use standard IGA, where both geometry and solution fields are approximated with same NURBS basis functions. Let $\boldsymbol{x}\in \Omega$, and $\boldsymbol{\xi}$ be the corresponding point in parametric domain as shown in fig:NURBS parameterization, then the domain using $n$ NURBS functions $N_{i}$ and $n$ control points $\mathbf{X}_i$ is approximated as,

\begin{equation}
\boldsymbol{x} = \sum _{i=1}^{n}\mathbf{X}_{i}N_{i}(\boldsymbol{\xi}),
\label{eq:NURBs_approx}
\end{equation}

and accordingly, the test and trial functions are also approximated with the same NURBS shape functions as,

\begin{equation}
T^h(\boldsymbol{\xi}) = \sum_{i=1}^{n}T_{i}N_{i}(\boldsymbol{\xi}), \quad \textrm{and} \quad
S^h(\boldsymbol{\xi}) = \sum_{i=1}^{n}S_{i}N_{i}(\boldsymbol{\xi}),
\label{eq:Trail test approx}
\end{equation}

Here, $T_i$ and $S_{i}$ are the temperature and the arbitrary temperature at the $i^{th}$ control point.

By substituting eq:Trail test approx in eq:mmodified_weak_form, a linear system is obtained,

\begin{equation}
\label{eq:Linear matrix system}
\mathbf{K} \mathbf{T} = \mathbf{F},
\end{equation}

where $\mathbf{T}$ is the vector of unknown temperatures at the control points. The global stiffness matrix $\mathbf{K}$ and the global flux vector $\mathbf{F}$ are written as,

\begin{equation} \label{eq:K}
\mathbf{K}= \mathbf{K}^b + \mathbf{K}^n+ (\mathbf{K}^n)^{\rm T}+\mathbf{K}^s,
\end{equation}

\begin{equation}
\mathbf{F}= \int _{\Gamma_{N}} \mathbf{N}^{\textrm{T}} Q_N~d\Gamma,
\end{equation}

where $\mathbf{K}^b$ is the bulk stiffness matrix. As $\Omega_{\rm in}$, $\Omega_{\rm out}$ and $\Omega_{\rm design}$ are considered separate patches, $\mathbf{K}^b$ is defined as follows,

\begin{equation}
\mathbf{K}^b =\sum_{k \in {\rm \{in,design,out \}}} \int _{\Omega_k} (\mathbf{B}^k)^{\textrm{T}}\boldsymbol{\kappa}^k(\varPhi)\mathbf{B}^k~d\Omega.
\label{eq:Kb}
\end{equation}

where $\mathbf{B}$ is the shape function derivative matrix (with superscript $k$ representing the patch index).
$\mathbf{K}^n$ and $\mathbf{K}^s$ are the interfacial stiffness matrices, since these matrices are used to couple adjacent patches with the conditions given in eq:Interface conditions a-(<ref>). A point to note, before defining $\mathbf{K}^n$ and $\mathbf{K}^s$, is that the connecting patches have matching control points at the interface as shown in fig:BVP domain. Following the notation
used in sec:Boundary value problem, the connecting patches at interface $\Gamma_I$ are denoted as 1 and 2 (and by superscripts 1 & 2 in the following equations). Consequently, $\mathbf{K}^n$ and $\mathbf{K}^s$ are given by the following equations,

\begin{equation}
\mathbf{K}^n =
\begin{bmatrix}
-\gamma\displaystyle\int_{\Gamma_I} (\mathbf{N}^1)^{\textrm{T}}\boldsymbol{n}\boldsymbol{\kappa}^1(\varPhi)\mathbf{B}^1~d\Gamma
& -(1-\gamma)\displaystyle\int_{\Gamma_I} (\mathbf{N}^1)^{\textrm{T}}\boldsymbol{n}\boldsymbol{\kappa}^2(\varPhi)\mathbf{B}^2~d\Gamma    \\[0.1em]
\gamma\displaystyle\int_{\Gamma_I} (\mathbf{N}^2)^{\textrm{T}}\boldsymbol{n}\boldsymbol{\kappa}^1(\varPhi)\mathbf{B}^1~d\Gamma  & (1-\gamma)\displaystyle\int_{\Gamma_I} (\mathbf{N}^2)^{\textrm{T}}\boldsymbol{n}\boldsymbol{\kappa}^2(\varPhi)\mathbf{B}^2~d\Gamma
\end{bmatrix},
\label{eq:Kn}
\end{equation}

\begin{equation}
\mathbf{K}^s =
\begin{bmatrix}
\beta\displaystyle\int_{\Gamma_I} (\mathbf{N}^1)^{\textrm{T}}\mathbf{N}^1~d\Gamma
& -\beta\displaystyle\int_{\Gamma_I} (\mathbf{N}^1)^{\textrm{T}}\mathbf{N}^2~d\Gamma  \\[0.1em]
-\beta\displaystyle\int_{\Gamma_I} (\mathbf{N}^2)^{\textrm{T}}\mathbf{N}^1~d\Gamma  & \beta\displaystyle\int_{\Gamma_I} (\mathbf{N}^2)^{\textrm{T}}\mathbf{N}^2~d\Gamma
\end{bmatrix},
\label{eq:Ks}%
\end{equation}

where $\mathbf{N}$ is the vector of shape functions. For a given patch $k$, the shape function derivative matrix $\mathbf{B}^k$ and the vector of shape functions $\mathbf{N}^k$ are given as follows,

\begin{equation}
\mathbf{B}^k =
\begin{bmatrix}
{N\strut_{1,x}^k} & {N\strut_{2,x}^k}&... & {N\strut_{I,x}^k} &...   \\[0.1em]
{N\strut_{1,y}^k} & {N\strut_{2,y}^k}&... & {N\strut_{I,y}^k} &...
\end{bmatrix}, \quad
\mathbf{N}^k =
\begin{bmatrix}
{N\strut_{1}^k} & {N\strut_{2}^k} & ... & {N\strut_{I}^k} & ...
\end{bmatrix}.
\end{equation}

In sec:Sensitivity analysis, the derivative of global stiffness matrix with respect to LSF $\varPhi$ will be needed in the sensitivity calculation. It is defined by differentiating  eq:K as follows,

\begin{equation} \label{eq:dKdphi}
\dfrac{d\mathbf{K}}{d\varPhi}= \dfrac{d\mathbf{K}^b}{d\varPhi} + \dfrac{d\mathbf{K}^n}{d\varPhi}+ \left(\dfrac{d\mathbf{K}^n}{d\varPhi}\right)^{\rm T}+\dfrac{d\mathbf{K}^s}{d\varPhi},
\end{equation}

\begin{equation}\label{eq:dKbdphi}
\dfrac{d\mathbf{K}^b}{d\varPhi}=\sum_{k \in {\rm \{in,design,out \}}}\int _{\Omega_k} (\mathbf{B}^k)^{\textrm{T}}\dfrac{d\boldsymbol{\kappa}^k(\varPhi)}{d\varPhi}\mathbf{B}^k~d\Omega,
\end{equation}

and $\dfrac{d\mathbf{K}^n}{d\varPhi}$ & $\dfrac{d\mathbf{K}^s}{d\varPhi}$ are defined similarly by differentiating eq:Kn & eq:Ks, respectively.

§.§ Level set function parameterization

The discretization strategy to decouple the LSF parametrization from the geometry and  solution field parameterizations. Two stages of refinement are provided. The first is to create the design/LSF mesh, and the second is to create the solution mesh to ensure the solution accuracy.

The objective of LSF parameterization is to use the basis functions (in our case, the NURBS basis functions) to approximate LSF on the entire domain via values at nodes/control points. Values at nodes/control points are called expansion coefficients and are utilized as design variables for optimization. The choice of LSF parameterization also decides the design freedom as well as the detailedness of the level set interface, and eventually influences the optimization results. A finer mesh for LSF parameterization often means a more complex optimization problem with higher design variables, which naturally requires more computational effort. On the other hand, the numerical analysis method requires a finer mesh to ensure adequate solution accuracy. Taking this into account, it is advantageous to decouple the LSF parameterization from the structural mesh. This way, the structural mesh can be refined without changing the design variables.
In this paper, we utilize two stages of refinement to achieve the above-mentioned decoupling. After creating a geometry on the coarsest possible NURBS mesh, we provide the first stage of refinement to get the required LSF parameterization (or the required number of design variables). After the first stage, the mesh is refined again for the second stage to ensure the required structural accuracy. The detailed refinement procedure is shown in fig:mesh independancy. The degree elevation and knot insertion algorithms are used for both stages of refinement. All refinements are provided uniformly. For future work, it will be advantageous to define a refinement criterion <cit.> based on the LSF, and employ local refinement using advanced splines <cit.>. The LSF $\varPhi$ is parameterized using $m$ NURBS basis functions $R_i$ as,

\begin{equation}\label{eq:LSF parameterization}
\varPhi(\boldsymbol{\xi})=\sum_{i=1}^m R_i(\boldsymbol{\xi})\varPhi_i
\end{equation}

where $\varPhi_i$ is the expansion coefficient corresponding to the $i^{th}$ control point.
We often need the expansion coefficients corresponding to a predefined LSF, for example, to initialize the optimization. To obtain these expansion coefficients $\varPhi_i$, a simple linear mass system is solved. The mass system can be written as,

\begin{equation}
\label{eq:Mass matrix system}
\mathbf{M}\boldsymbol{\Phi}= \mathbf{\Psi},
\end{equation}

where $\boldsymbol{\Phi}=[\varPhi_1\hspace{0.7em}\varPhi_2\hspace{0.7em}...\hspace{0.7em}\varPhi_m]^{\rm T}$, mass matrix $\mathbf{M}$ and the right-side vector $\mathbf{\Psi}$ are defined as,

\begin{equation}
\label{eq:Mass matrix}
\mathbf{M}=\int_{\Omega_{\rm design}} (\mathbf{R})^{\textrm{T}}\mathbf{R}~d\Omega,
\end{equation}

\begin{equation}
\label{eq:Rhs vector}
\mathbf{\Psi}= \int_{\Omega_{\rm design}}(\mathbf{R})^{\textrm{T}} \varPhi~d\Omega.
\end{equation}

§ OPTIMIZATION PROBLEM

§.§ Optimization problem description

In the level set topology optimization method, the expansion coefficients, as mentioned in eq:LSF parameterization, are used as design variables. The goal is to find the values of these design variables, such that the corresponding topology yields the optimal value of the function of interest (called objective function). In our case, $\boldsymbol{\Phi}=[\varPhi_1\hspace{0.7em}\varPhi_2\hspace{0.7em}...\hspace{0.7em}\varPhi_{N_{\rm var}}]^{\rm T}$ is the vector of the $N_{\mathrm{var}}$ design variables,  and $J$ is the objective function. In a general case, $N_{\text{var}}$ can be different from the number of basis functions $m$ in eq:LSF parameterization. For the most of numerical examples in the next section, $N_{\text{var}}$ are less than $m$ considering imposed $x$ and $y$ axial symmetry. The topology optimization problem for a heat manipulator can be defined in a mathematical form as,

\begin{equation}
\label{eq:optimization problem}
\min_{\boldsymbol{\Phi} \in  \mathbb{R}^{N_{\mathrm{var}}} }~J(T^h,\varPhi),
\end{equation}

J : ℝ^N_var  →ℝ,

J : Φ ↦J(T^h(Φ),(Φ)),

such that the following constraints are satisfied,

\begin{align}
&\text{Equality constraint:}\quad &a(T^h,S^h,\varPhi)= \ell(S^h), \forall S^h \in \mathscr{S}^h_0 \quad &\textrm{in} \quad \Omega, \label{eq:equality constraint_1}\\
&\text{Equality constraint:}\quad &T= T_D \quad &\textrm{on} \quad \Gamma_D,   \label{eq:equality constraint_2}\\
% &\text{Inequality constraints:}\quad &g_{j}(\boldsymbol{x})\leq~0 \quad
% % & j=1,2,...,n_g , \label{eq:Inequality constraints}\\
&\text{Box constraints:}\quad &\varPhi_{i,\rm{min}}\leq~\varPhi_{i}\leq~\varPhi_{i,\rm{max}} \quad &i=1,2,...,N_{\mathrm{var}} ,  \label{eq:Box constraints}
\end{align}

where $\varPhi_{i,\rm{min}}$ and $\varPhi_{i,\rm{max}}$ are the lower and upper bounds of the design variable $\varPhi_{i}$.
As mentioned in sec:Introduction, the given optimization problem is solved using a mathematical programming technique - Sequential Quadratic programming (SQP). In SQP, the boundary value problem described in sec:Boundary value problem is solved in each iteration. Since SQP is a gradient-based algorithm, the sensitivity of the objective functions $J$ with respect to each design variable $\varPhi_i$ needs to be calculated (using the solution of the boundary value problem) at the end of each iteration. Later, the sensitivities are fed into the algorithm to generate the new values of the design variables. The following two sections will, respectively, explain the adjoint method to calculate the sensitivity at each iteration and SQP in detail.

§.§ Sensitivity analysis

In this section, we outline the sensitivity analysis method called the adjoint method <cit.>. In the most general case, the performance objective function $J(T,\varPhi)$ can be written as a sum of two terms, corresponding to the volume and the surface integrals, respectively, i.e.

\begin{equation}
J(T,\varPhi)= \int _{\Omega_b} J_b(T,\varPhi)~d\Omega + \int _{\Gamma_s} J_s(T,\varPhi)~d\Gamma,
\label{eq:Objective fun definition}
\end{equation}

where ${\Omega_b}$ is the domain where the volume term is calculated, and ${\Gamma_s}$ is part of the boundary where the surface term is calculated. For our numerical examples, only the volume term of the objective function is considered.
Next, the Lagrangian can be obtained by augmenting the objective functional $J(T,\varPhi)$ with the weak form constraint as well as Dirichlet boundary constraint (with the help of the Lagrange multipliers $P$ and $\lambda$, respectively) that $T$ should satisfy. It is defined as,

\begin{equation}\label{eq:Langrangian 1}
\mathcal{L}:\mathbb{H}^1(\mathbb{R}^2) \times \mathbb{H}^1(\mathbb{R}^2) \times \mathbb{H}^1(\mathbb{R}^2) \times \mathbb{H}^1(\mathbb{R}^2) \rightarrow \mathbb{R},
\end{equation}

\begin{equation}\label{eq:Langrangian 2}
\mathcal{L}(T,P,\varPhi,\lambda)=J(T,\varPhi)+a(T,P,\varPhi)-\ell(P)+\int_{\mathrm{\Gamma}_D}\lambda (T-T_D)~d\mathrm{\Gamma}.
\end{equation}

The optimality conditions of the minimization problem are derived as the stationary conditions of the Lagrangian. The first stationary condition is obtained by equating the partial Fréchet derivative of $\mathcal{L}$ with respect to $P$ (in any arbitrary direction $\delta P \in \mathbb{H}^1$) to zero, i.e.

\begin{equation} \label{eq:Langrangian der wrt P 1}
\left<\dfrac{\partial \mathcal{L}(T,P,\varPhi,\lambda)}{\partial P}, \delta P \right> = \left<\dfrac{\partial a(T, P, \varPhi)}{\partial P}, \delta P \right> -  \left<\dfrac{\partial \ell(P)}{\partial P}, \delta P \right>=0.
\end{equation}

Using the bilinear form and the linear form from eq:mmodified_weak_form_a and eq:weak_form_l, it can be derived that,

\begin{multline} \label{eq:bilinear-linear form properties}
\left<\dfrac{\partial a(T, P, \varPhi)}{\partial T}, \delta T \right>= a(\delta T, P, \varPhi); \quad \left<\dfrac{\partial a(T, P, \varPhi)}{\partial P}, \delta P \right>= a(T,\delta P, \varPhi);\\ \text{and} \quad \left<\dfrac{\partial \ell(P)}{\partial P}, \delta P \right>=\ell(\delta P).
\end{multline}

Substituting the second and third equations of eq:bilinear-linear form properties in eq:Langrangian der wrt P 1,

\begin{equation} \label{eq:Langrangian der wrt P 2}
a(T, \delta P, \varPhi)= \ell(\delta P),
\end{equation}

which shows that $\forall \delta P \in \mathbb{H}^1$, given state variable $T \in \mathbb{H}^1$ satisfies the weak form.
The second stationary condition is obtained by taking the partial Fréchet derivative of $\mathcal{L}$ with respect to $\lambda$ (in any arbitrary direction $\delta \lambda \in \mathbb{H}^1$) and equating it to zero as follows,

\begin{equation} \label{eq:Langrangian der wrt lambda 1}
\left<\dfrac{\partial \mathcal{L}(T,P,\varPhi,\lambda)}{\partial \lambda}, \delta \lambda \right> = \int_{\mathrm{\Gamma}_D}\delta \lambda (T-T_D)~d\mathrm{\Gamma} =0.
\end{equation}

which readily implies that $T=T_D$ on $\mathrm{\Gamma}_D$. Together eq:Langrangian der wrt P 2 and eq:Langrangian der wrt lambda 1 solve the boundary value problem in the state variable $T$ (as described in eq:Heat conduction BVP).
The third stationary condition is obtained by taking the partial Fréchet derivative of $\mathcal{L}$ with respect to $T$ (in any arbitrary direction $\delta T \in \mathbb{H}^1$) and equating it to zero as follows,

\begin{equation} \label{eq:Langrangian der wrt T 1}
\left<\dfrac{\partial \mathcal{L}(T,P,\varPhi,\lambda)}{\partial T}, \delta T \right> =  \left<\dfrac{\partial J(T,\varPhi)}{\partial T}, \delta T \right> + \left<\dfrac{\partial a(T, P, \varPhi)}{\partial T}, \delta T \right> +\int_{\mathrm{\Gamma}_D}\lambda \delta T~d\mathrm{\Gamma}=0.
\end{equation}

Using eq:bilinear-linear form properties & the fact that $T$ is constant, i.e.$\delta T=0$, on $\mathrm{\Gamma}_D$,

\begin{equation} \label{eq:Langrangian der wrt T 2}
a(\delta T, P, \varPhi)=-\left<\dfrac{\partial J(T,\varPhi)}{\partial T}, \delta T \right>.
\end{equation}

eq:Langrangian der wrt T 2 is valid for $\forall\delta T \in \mathbb{H}^1$, therefore, with an arbitrary normal derivative ${\small \dfrac{\partial \delta T}{d\boldsymbol{n}}}$ on $\mathrm{\Gamma}_D$ and vanishing trace $\delta T=0$, it gives, $P=0$ on $\mathrm{\Gamma}_D$. Combining this fact with eq:Langrangian der wrt T 2, we can generate a well-posed adjoint problem as,

\begin{align}
\nabla \cdot \left( \boldsymbol{\kappa}\nabla P\right) &= -\dfrac{\partial J(T,\varPhi)}{\partial T} \quad &&\textrm{in} \quad \mathrm{\Omega}, \label{eq:Adjoint Laplace equation}\\
P &= 0 \quad &&\textrm{on} \quad \mathrm{\Gamma}_D, \label{eq:Adjoint Boundary conditions a}\\
\nabla P\cdot \boldsymbol{n} &= 0 \quad &&\textrm{on} \quad \mathrm{\Gamma}_N, \label{eq:Adjoint Boundary conditions b}\\
\left\llbracket P\right\rrbracket &=0 \quad &&\textrm{on} \quad \mathrm{\Gamma}_I, \label{eq:Adjoint Interface conditions a}\\
\boldsymbol{n}\cdot \left\llbracket \boldsymbol{\kappa}\nabla P\right\rrbracket &= 0 \quad &&\textrm{on} \quad \mathrm{\Gamma}_I, \label{eq:Adjoint Interface conditions b}
\end{align}

where $T$ is the primary temperature field, and $P$ is the adjoint temperature field. By employing the trial and test  function approximations, the adjoint problem is discretized into the following linear system,

\begin{equation} \label{eq:adjoint eq. matrix form}
\mathbf{K}^{\rm T} \mathbf{P}= \mathbf{F}_{\rm adj},
\end{equation}

where $\mathbf{P}$ is the vector of adjoint temperature at control points, $\mathbf{F}_{\rm adj}$ is the global adjoint flux vector defined as,

\begin{equation}
\mathbf{F}_{\rm adj}=-\int _{\Omega_b} \mathbf{N}^{\textrm{T}} \dfrac{d J_b}{d T}~d\Omega - \int _{\Gamma_s} \mathbf{N}^{\textrm{T}} \dfrac{d J_s}{d T}~d\Gamma,
\end{equation}

At last, with the fulfillment of all three stationary conditions, the sensitivity of the objective functional $\left(\dfrac{dJ}{d\varPhi}\right)$ becomes equal to the total derivative of Lagrangian $\mathcal{L}$ with respect to the LSF $\varPhi$. Therefore, it can be written as follows,

\begin{multline}\label{eq:Langrangian der wrt phi 1}
\left<\dfrac{d J(T,\varPhi)}{d \varPhi},\delta \varPhi \right>:=\left<\dfrac{d \mathcal{L}(T,P,\varPhi,\lambda)}{d \varPhi}, \delta \varPhi \right> =  \left<\dfrac{\partial J(T,\varPhi)}{\partial \varPhi}, \delta \varPhi \right>+ \left<\dfrac{\partial J(T,\varPhi)}{\partial T}, \delta T \right> \\+ \left<\dfrac{\partial a(T, P, \varPhi)}{\partial \varPhi}, \delta \varPhi \right> + \left<\dfrac{\partial a(T, P, \varPhi)}{\partial T}, \delta T \right>.
\end{multline}

By substituting eq:Langrangian der wrt T 2,

\begin{equation}\label{eq:Langrangian der wrt phi 2}
\left<\dfrac{d J(T,\varPhi)}{d \varPhi},\delta \varPhi  \right>=\left<\dfrac{\partial J(T,\varPhi)}{\partial \varPhi}, \delta \varPhi \right> + \left<\dfrac{\partial a(T, P, \varPhi)}{\partial \varPhi}, \delta \varPhi \right>,
\end{equation}

\begin{equation}\label{eq:Langrangian der wrt phi 3}
\therefore \dfrac{d J(T,\varPhi)}{d \varPhi}=\dfrac{\partial J(T,\varPhi)}{\partial \varPhi} +\dfrac{\partial a(T, P, \varPhi)}{\partial \varPhi}.
\end{equation}

By employing trial and test function approximations, eq:Langrangian der wrt phi 3 can be written as,

\begin{equation} \label{eq:dJ_dPhi matrix form}
\dfrac{d J}{d\varPhi} = \int _{\Omega_b} \dfrac{d J_b}{d\varPhi}~d\Omega + \int _{\Gamma_s} \dfrac{d J_s}{d\varPhi}~d\Gamma +(\mathbf{P})^{\rm T}~\dfrac{d \mathbf{K}}{d\varPhi}~\mathbf{T},
\end{equation}

where $\dfrac{d \mathbf{K}}{d\varPhi}$ is the derivative of the global stiffness matrix with respect to $\varPhi$ (as given in eq:dKdphi). Using the LSF parameterization (eq:LSF parameterization), the sensitivity with respect to a particular expansion coefficient $\varPhi_i$ is given as,

\begin{equation} \label{eq:dJ_dPhi matrix form_2}
\dfrac{d J}{d\varPhi_i} = \int _{\Omega_b} \dfrac{d J_b}{d\varPhi}~\dfrac{d \varPhi}{d\varPhi_i}~d\Omega + \int _{\Gamma_s} \dfrac{d J_s}{d\varPhi}~\dfrac{d \varPhi}{d\varPhi_i}~d\Gamma +(\mathbf{P})^{\rm T}~\dfrac{d \mathbf{K}}{d\varPhi}~\dfrac{d \varPhi}{d\varPhi_i}~\mathbf{T}.
\end{equation}

§.§ Update scheme - Sequential Quadratic Programming (SQP)

In the Sequential Quadratic Programming (SQP) approach, the optimization problem is approximated as a quadratic sub-problem at each iteration <cit.>.  For the SQP method, the main challenge is to form a good quadratic sub-problem that generates a suitable step for optimization. For our optimization problem, the quadratic subproblem in the $m^{th}$ iteration is modeled as,

\begin{equation}
\label{eq:Quadratic problem}
\min_{\boldsymbol{p} \in  \mathbb{R}^{N_{\mathrm{var}}} }~J_m+ \nabla J_m^{\rm T}\boldsymbol{p}+\dfrac{1}{2}\boldsymbol{p}^{\rm T}\boldsymbol{\mathcal{H}}_m\boldsymbol{p},
\end{equation}

where $J_m$, $\nabla J_m$ and $\boldsymbol{\mathcal{H}}_m$ are the objective function value, the gradient of the objective function, and the Hessian approximation of the Lagrangian $\mathcal{L}$ in the $m^{th}$ iteration, respectively.  Later, the minimizer of the subproblem  $\boldsymbol{p}=\boldsymbol{p}_m$ is used to get the next approximation $\boldsymbol{\Phi}_{m+1}$ based on a line search algorithm or a trust region algorithm. One point to note here is that equality constraints are excluded. As the temperature $T$ is evaluated from the linear system eq:Linear matrix system, the equality constraints will be satisfied explicitly at each iteration.
In this work, we use `fmincon' optimization toolbox from MATLAB with the `sqp' option. The `sqp' option solves the above-mentioned quadratic subproblem with an active set strategy. At each iteration, the user provides the objective function $J_m$ and the gradient $\nabla J_m$ (using the sensitivity analysis from the last section). For the calculation of Hessian approximation $\boldsymbol{\mathcal{H}}_m$, it uses the BFGS-update scheme. Once the solution to the subproblem $\boldsymbol{p}_m$ is found, the size of step $\alpha_m$ is found by using an inbuilt line search algorithm. Then, the design variables are updated as,

\begin{equation}
\boldsymbol{\Phi}_{m+1} = \boldsymbol{\Phi}_{m} + \alpha_m   \boldsymbol{p}_m.
\end{equation}

§.§ Regularizations

More often, the optimization problem is found to suffer from a lack of well-posedness, numerical artifacts, slower convergence, and entrapment into local minima with poor performance. To overcome these issues, regularization techniques can be used. Sometimes, regularization is utilized to control geometrical features in optimization results. Generally, a single regularization technique addresses more than one of the above-mentioned problems. For our numerical examples, we explored three regularization techniques: LSF reinitialization, Tikhonov regularization, and volume regularization.

§.§.§ Level set function reinitialization

As the interfaces are defined by the zero-level contours of LSF, only local region near these contours are uniquely described in the optimal solution. Therefore, the LSF is not unique especially in the region far from interfaces. As a consequence of this non-uniqueness property, the LSF sometimes becomes too flat or too steep during optimization. This phenomenon can deteriorate the convergence rate. To alleviate the problem, the LSF is initialized as well as maintained as a sign distance function ($||\nabla \varPhi||_{2} = 1$) during optimization. One way to maintain the LSF as a sign distance function is to reinitialize it after several iterations while maintaining the locations of zero-level contours.

$\varPhi$ before reinitialization. Interface points are evaluated using LSF parameterization.

Evaluation of sign distance from integration points to closest interface points.

Reinitialized $\varPhi$ after solving the mass matrix system using the sign distance at integration points.

Geometry-based constrained LSF reinitialization technique

Here, we utilize the geometry-based constrained reinitialization technique. fig:Geometry-based constrained LSF reinitialization outlines the steps of the technique. Using the inverse of LSF parameterization, we find several points on the interface as shown in fig:Geometry-based constrained LSF reinitialization a. For the sake of clarity, we call them the interface points. Then, the new LSF value at a point is defined as the distance to the closest interface point while keeping the sign from the original value (fig:Geometry-based constrained LSF reinitialization b). The number of interface points decides the accuracy of the reinitialization. To find the new expansion coefficients, we solve the mass matrix system as shown in eq:Mass matrix system by implementing the new LSF values at the integration/collocation points (fig:Geometry-based constrained LSF reinitialization c). However, there is one difference from the initialization mass system. In this system, we apply the constraints which enforce the LSF values at interface points to be zero to preserve the location of the interface via the penalty method. An example of the LSF before and after reinitialization is shown in fig:levelset reinitialization.

$\varPhi$ before reinitialization

$\varPhi$ after reinitialization

$\varPhi$ after reinitialization in 2D

The LSF before and after reinitialization. The green points (in (c)) represent the interface points before reinitialization, which lie exactly on the new interface.

§.§.§ Tikhonov regularization

The effect of Tikhonov regularization <cit.> is similar to perimeter regularization, where the perimeter of the level-set interface is penalized. In Tikhonov regularization, an additional term penalizing the gradient of LSF is added to the main objective function. By penalizing the gradient of LSF, the optimization is directed towards smoother LSF and eventually smoother interface. between This Tikhonov term and corresponding sensitivity are written as,

\begin{equation}
\label{eq:Tikhonov regularization}
J_{\rm Tknv}= \int_{\Omega_{\rm design}} \nabla\varPhi^{\rm T}~\nabla\varPhi~d\Omega, \quad \text{and} \quad \dfrac{dJ_{\rm Tknv}}{d\varPhi_i}= \int_{\Omega_{\rm design}} 2\nabla\varPhi^{\rm T}~\nabla R_i~d\Omega.
\end{equation}

To implement Tikhonov regularization, the primary objective function and sensitivities are augmented by corresponding Tikhonov regularization terms using a weighing parameter $\chi$ as,

\begin{equation}
J_{\rm total}= J + \chi J_{\mathrm{Tknv}}, \quad \text{and} \quad \dfrac{dJ_{\rm total}}{d\varPhi_i}=\dfrac{dJ}{d\varPhi_i}+\chi \dfrac{dJ_{\rm Tknv}}{d\varPhi_i}.
\end{equation}

§.§.§ Volume regularization

For heat manipulator optimization, we propose one more regularization, called volume regularization. In a metamaterial-based heat manipulator (made of two member materials), often one high $\kappa$ material and another low $\kappa$ material are used. The high $\kappa$ material works as the medium that allows the flow of the flux, while the low $\kappa$ material hinders the flow and guides it towards the required path. By taking into account these particular roles, the area where the flux is not flowing (and the temperature gradient is zero or negligibly small) can be filled by low $\kappa$ material. In other words, the areas, that do not contribute to the original objective, can be filled with low $\kappa$ material. By doing so, those areas are made as homogeneous as possible and free of unnecessary complex features. Mathematically, this objective is defined through a volume term as,

\begin{equation}
\label{eq:Volume regularization}
J_{\rm vol}= \int_{\Omega_{\rm design}} H(\varPhi)~d\Omega, \quad \text{and} \quad  \dfrac{dJ_{\rm vol}}{d\varPhi_i}= \int_{\Omega_{\rm design}} \delta(\varPhi)R_i~d\Omega.
\end{equation}

The term calculates the volume of high $\kappa$ material, as we assign high $\kappa$ material on the positive side of the level set interface. As our optimization problem is a minimization problem, the term guides the optimization toward a larger volume of low $\kappa$ material. When this term is added to the main objective function through a weighing parameter, the combined objective of heat manipulation and volume filling can be achieved simultaneously.
To implement volume regularization, the primary objective function and sensitivities are augmented by corresponding volume regularization terms using a weighing parameter $\rho$ as,

\begin{equation}
J_{\rm total}= J + \rho J_{\mathrm{Tknv}}, \quad \text{and} \quad \dfrac{dJ_{\rm total}}{d\varPhi_i}=\dfrac{dJ}{d\varPhi_i}+\rho\dfrac{dJ_{\rm Tknv}}{d\varPhi_i}.
\end{equation}

§ NUMERICAL EXAMPLES

In this section, we verify the proposed method, and test its efficiency and effectiveness for several numerical examples. We denote $p$ and $q$ as the order of NURBS approximation in two directions of a NURBS patch, which lies along circumferential direction and radial directions respectively for the given examples.
Since we are using `fmincon' optimization tool box in MATLAB, there are several inbuilt stopping criteria such as, `OptimalityTolerance' (tolerance value in first-order optimality measures), `StepTolerance' (tolerance value of the change in design variables' values), `ObjectiveLimit' (tolerance value of the objective function), `MaxFunctionEvaluations' (maximum number of function evaluations), `MaxIterations' (maximum number of iterations). We define `ObjectiveLimit'=$1\times 10^{-9}$, `StepTolerance'=$1\times 10^{-8}$, `OptimalityTolerance'=$1\times 10^{-6}$. If the LSF reinitialization is utilized, `MaxFunctionEvaluations' and  `MaxIterations' are employed to stop the optimization before reinitializing it. Also, the `StepTolerance' criterion often gets affected by poor LSF function. Therefore, the optimization is stopped when it fails by `StepTolerance' criterion for 4 consecutive times. Otherwise, the LSF is reinitialized and the optimization is run again. The stopping criteria, tolerance values, and applied regularizations differ slightly example-wise, and are mentioned in the descriptions of the examples. Most of the examples are symmetric along $x$ and $y$-axes. Therefore, $x$ and $y$-axes symmetry is applied in almost all numerical models, and we explicitly state when there is an exception. One more thing to note, from here onward, is that the value of the objective function for a particular test case will be mentioned in the caption.

§ CONCLUSIONS

§ APPENDIX

§ ANALYTICAL SOLUTIONS OF THE ANNULAR RING PROBLEM

The analytical solutions of the boundary value problem (eq:BenchAnn BVP) is given as follows,

\begin{align}
&T_{a}(r)=c_{a} +d_{a} \log(r),\\
&T_{b}(r)=c_{b} +d_{b} \log(r),
\end{align}

where coefficients $c_a, d_a, c_b, d_b$ are obtained from the following system of algebraic equations, resulted from satisfying boundary conditions (eq:BenchAnn BVP) for $r = R_a$, $r = R_b$ and $r = R_L$,

\begin{equation}
\begin{bmatrix}
1 & \log(R_a) & 0 & 0\\
0 & 0 & 1 & \log(R_b) \\
1 & \log(R_L) & -1 & -\log(R_L) \\
0 & \kappa_a & 0 & -\kappa_b
\end{bmatrix} \begin{bmatrix}
c_a \\
d_a \\
c_b \\
\end{bmatrix} = \begin{bmatrix}
T_a \\
T_b \\
0 \\
\end{bmatrix},
\end{equation}

which can be written as,

\begin{equation}
\begin{aligned}
c_{a}&=-\frac{-\kappa_{b}T_{b} \log(R_{a})-\kappa_{a}T_{a} \log(R_{I})+\kappa_{b}T_{a} \log(R_{I})+\kappa_{a}T_{a}\log(R_{b})}{\kappa_{b}\log(R_{a})+\kappa_{a}\log(R_{I})-\kappa_{b}\log(R_{I})-\kappa_{a}\log(R_{b})},\\
c_{b}&=-\frac{\kappa_{b}T_{b} \log(R_{a})+\kappa_{a}T_{b} \log(R_{I})-k_{b}T_{b} \log(R_{I})-\kappa_{a}T_{a}\log(R_{b})}{-\kappa_{b}\log(R_{a})-\kappa_{a}\log(R_{I})+\kappa_{b}\log(R_{I})+\kappa_{a}\log(R_{b})},\\
\end{aligned}
\end{equation}

The adjoint system (eq:BenchAnn adjoint) can be reduced to the following two ODEs,

\begin{equation}
\begin{aligned}
P_a^{''}(r) + \dfrac{1}{r} P_a^{'}(r) &= \dfrac{-2}{\kappa_a}(c_a + d_a \log(r)),\\
P_b^{''}(r) + \dfrac{1}{r} P_b^{'}(r) &= \dfrac{-2}{\kappa_b}(c_b + d_b \log(r)),
\end{aligned}
\end{equation}

whose solution takes the form,

\begin{equation}
\begin{aligned}
P_a(r) &= C_1 \log (r) + D_1 - \dfrac{(c_a - d_a) r^2}{2 \kappa_a} +\dfrac{d_a r^2 \log (r)}{2 \kappa_a}, \\
P_b(r) &= C_2 \log (r) + D_2 - \dfrac{(c_b - d_b) r^2}{2 \kappa_b} +\dfrac{d_b r^2 \log (r)}{2 \kappa_b},
\end{aligned}
\end{equation}

The coefficients $C_1, D_1, C_2, D_2$ are found by satisfying boundary conditions given in eq:BenchAnn adjoint, i.e. from the following system,

\begin{multline}
\begin{bmatrix}
1 & \log(R_a) & 0 & 0\\
0 & 0 & 1 & \log(R_b) \\
1 & \log(R_L) & -1 & -\log(R_L) \\
0 & \kappa_a & 0 & -\kappa_b
\end{bmatrix}
\begin{bmatrix}
C_1 \\
D_1 \\
C_2 \\
\end{bmatrix} \\
= \begin{bmatrix}
\dfrac{R_a^2 (c_a - d_a) + R_a^2 d_a \log (R_a)}{2 \kappa_a} \\
\dfrac{R_b^2 (c_b - d_b) + R_b^2 d_b \log (R_b)}{2 \kappa_b} \\
\dfrac{R_L^2}{2}\left(\dfrac{c_a - d_a}{\kappa_a} - \dfrac{c_b - d_b}{\kappa_b}\right) + \dfrac{R_L^2\log(R_L)}{2}\left(\dfrac{d_a}{\kappa_a} - \dfrac{d_b}{\kappa_b}\right)\\
\dfrac{R_L^2}{2} \left(2 c_a - 2 c_b - d_a + d_b\right) + R_L^2 \log(R_L) \left(d_a - d_b\right)
\end{bmatrix}.
\end{multline}