## 2 Materials and Methods  
A new approach to identifying an arbitrary number of inclusions, their geometry and their location in 2D and 3D structures was based on the influence of temperature and heat flux on the structure. This solves the inverse heat conduction problem.  
Consider an isotropic flat cuboid body (Figure 1), inside which there are inclusions occupying the area \(V=\{x\in[0;a],\;y\in[0;b],\;z\in[0;c]\}\;V_{1},\;V_{2},\ldots\;V_{m},\;m=1,2 \ldots.\)  
Figure 1: The cuboid isotropic flat body with selected areas \(V_{1},\;V_{2},\;\ldots V_{m}\).  
The body was under the action of both the temperature at the boundary \(B_{T}\) (the Dirichlet (or first type) boundary condition) and the heat flux at the boundary \(B_{q}\) (the Neumann (or second type) boundary condition). The temperature field \(T\) inside \(V\) satisfies the heat conduction equation [40]  
\[k\Delta T=0,\text{ in }V(x,y,z), \tag{1}\]  
and boundary conditions  
\[T|_{B_{T}}=T_{0}, \tag{2}\]  
\[k\frac{dT}{dn}\Big{|}_{B_{q}}=q_{0}, \tag{3}\]  
where \(\Delta(\cdot)=\
abla^{2}(\cdot)=\partial^{2}(\cdot)/\partial x^{2}+\partial^ {2}(\cdot)/\partial y^{2}+\partial^{2}(\cdot)/\partial z^{2}\), \(T_{0}\) is the set temperature at \(B_{T}\); \(q_{0}\) is the heat flux determined on \(B_{q}\); \(n\)--is the unit normal vector directed outward of the area. In the area of \(V\) defined subdomains \(V_{1}\), \(V_{2}\),..., \(V_{m}\) with the coefficient of thermal conductivity \(k_{1}\), \(k_{2}\),..., \(k_{m}\).  
To form the topological optimization procedure, the area \(V\) was divided into \(N\) finite elements. Then, each element was assigned a variable density, \(\rho_{n}(\mathbf{x})(n=1,N)\) which depends on the vector of design variables \(\mathbf{x}\)\((x,\)\(y,\)\(z)\). The thermal conductivity coefficient depends on the artificially introduced density \(\rho_{n}(\mathbf{x})\), which will be the control variable when using the Solid Isotropic Material with Penalty method (SIMP) [45].  
The thermal conductivity coefficient was introduced as follows:  
\[k(\mathbf{x})=k_{1}\rho_{n}^{p}(\mathbf{x})+k_{2}(1-\rho_{n}^{p}(\mathbf{x})) \tag{4}\]  
where \(k_{1},k_{2}\)--are thermal conductivity coefficients in volumes \(V_{m}\), respectively, and \(p\geq 1\) is a penalty factor for the SIMP method. The goal function was defined as follows [40]:  
\[f(\mathbf{x})=\chi\frac{h_{0}h_{\max}}{mes\left(V_{m}\right)}\int\limits_{V_{m }}|\
abla\rho(\mathbf{x})|^{2}dV_{m}+(1-\chi)\int\limits_{V_{m}}\Big{[}(T- \widetilde{T})^{2}+(q-\widetilde{q})^{2}\Big{]}dV_{m} \tag{5}\]  
where the first term is the penalty function to eliminate the \"checkboard\" effect; \(0<\chi<1\) is the coefficient for matching the target function and the penalty function; \(h_{0}\) represents the initial grid size; \(h_{\max}\) is the current grid size; \(mes\left(V_{m}\right)\) is the area of the inclusions; \(\rho\)\((\mathbf{x})\) is the physical density; \(T\) is the temperature on a flat area and \(\widetilde{T}\) is calculated temperature on a flat area; \(q\) is a given flow on an area and \(\widetilde{q}\) is the calculated flow on a flat area.  
Limitations for physical density \(\rho\)\((\mathbf{x})\) are chosen as the following:  
\[0\leq\int\limits_{V_{m}}\rho(\mathbf{x})\ d\ V_{m}\leq\gamma\cdot mes\ (V) \tag{6}\]  
where \(\gamma\) denotes the fraction of material in the inclusions. Identification consists of minimizing the difference in temperature distribution and heat fluxes in the original design, and obtaining them in this iteration step.