DeepAI
Log In Sign Up

Compact 200 line MATLAB code for inverse design in photonics by topology optimization: tutorial

08/23/2020
by   Rasmus E. Christiansen, et al.
DTU
0

We provide a compact 200 line MATLAB code demonstrating how topology optimization (TopOpt) as an inverse design tool may be used in photonics, targeting the design of two-dimensional dielectric metalenses and a metallic reflector as examples. The physics model is solved using the finite element method, and the code utilizes MATLAB's fmincon algorithm to solve the optimization problem. In addition to presenting the code itself, we briefly discuss a number of extensions and provide the code required to implement some of these. Finally, we demonstrate the superiority of using a gradient-based method compared to a genetic-algorithm-based method (using MATLAB's ga algorithm) for solving inverse design problems in photonics. The MATLAB software is freely available in the paper and may be downloaded from https://www.topopt.mek.dtu.dk.

READ FULL TEXT VIEW PDF
01/14/2019

Electrical Impedance Tomography based on Genetic Algorithm

In this paper, we applies GA algorithm into Electrical Impedance Tomogra...
01/25/2022

HoneyTop90: A 90-line MATLAB code for topology optimization using honeycomb tessellation

This paper provides a simple, compact and efficient 90-line pedagogical ...
08/10/2021

Topology optimization of thermal problems in a nonsmooth variational setting: closed-form optimality criteria

This paper extends the nonsmooth Relaxed Variational Approach (RVA) to t...
09/20/2020

Topology Optimization through Differentiable Finite Element Solver

In this paper, a topology optimization framework utilizing automatic dif...
01/08/2021

Topology Optimization with linearized buckling criteria in 250 lines of Matlab

We present a 250 line Matlab code for topology optimization for lineariz...
05/05/2020

A new generation 99 line Matlab code for compliance Topology Optimization and its extension to 3D

Compact and efficient Matlab implementations of compliance Topology Opti...

I Introduction

This paper details a 200 line MATLAB code, which demonstrates how density-based Topology Optimization (TopOpt) can be applied to photonics design. The code is written for scientists and students with a basic knowledge of programming, numerical modelling and photonics, who desire to start using inverse design in their research. We briefly detail the model of the physics, followed by the discretized TopOpt design problem (Sec. II). The MATLAB code is then explained in detail (Sec. III), followed by two application examples providing the reader with targets for reproduction (Sec. IV). Then, a number of possible extensions are discussed and code-snippets for easy implementation are provided along with design examples (Sec. V). Finally, we demonstrate the superiority of using gradient-based TopOpt compared to a genetic algorithm when solving a photonic design problem (Sec. VI).

Topology optimization BOOK_TOPOPT_BENDSOE as an inverse design tool was first developed in the context of solid mechanics in the late 1980s BENDSOE_KIKUCHI_1988. Since its inception the method has developed rapidly and expanded across most areas of physics ALEXANDERSEN_ANDREASEN_2020; DILGEN_ET_AL_2019; LUNDGAARD_2019; JENSEN_SIGMUND_2011. Over the last two decades the interest in applying TopOpt for photonics has increased rapidly MOLESKY_2018 with applications within cavity design LIANG_JOHNSON_2013; WANG_2018, photonic demultiplexers PIGGOTT_2015, metasurfaces ZIN_2019; CHUNG_MILLER_2020 and topological insulators CHRISTIANSEN_NP_2019 to name a few. While the interest in TopOpt within the photonics community has grown markedly in recent years, significant barriers hinder newcomers to the field from adopting the tool in their work. These are: the required knowledge of numerical modelling, of advanced mathematical concepts and scientific programming experience. This paper seeks to lower these barriers by providing the reader with a simple 2D finite element based MATLAB implementation of TopOpt for photonics, which is straightforwardly extendable to a range of other design problems. Within the field of structural optimization in mechanics similar simple MATLAB codes SIGMUND_2001; ANDREASSEN_2011; FERRARI_2020 have proven themselves highly successful in raising the awareness of TopOpt and serving as a basic platform for further expansion of the method, thus broadening its application as a design tool and successfully driving the field forward.

For readers who are less interested in the programming and method development aspects of TopOpt as an inverse design tool, we have authored a parallel tutorial paper on TopOpt for photonics applications, utilizing the GUI based commercial finite element software COMSOL Multiphysics as the numerical tool to model the physics and solve the optimization problem CHRISTIANSEN_SIGMUND_COMSOL_2020.

Ii The Physics and the Discretized Optimization Problem

We model the physics in the rectangular domain , with the boundary , (see Fig. 1) using Maxwell’s equations, assuming time-harmonic temporal behaviour. We define a subset of as the design domain and denote this region . We assume out-of-plane (

-direction) material invariance and that all involved materials are linear, static, homogeneous, isotropic, non-dispersive, non-magnetic and without inherent polarization. Finally, we assume out-of-plane polarization of the electric field (TE polarization). From these assumptions we derive a two-dimensional Helmholtz-type partial differential equation for the out-of-plane component of the electric field

in ,

(1)

where denotes the z-component of the electric field, is the wavenumber with lambda in the top200EM interface being the wavelength, denotes the relative electric permittivity and denotes a forcing term used to introduce an incident plane wave from the bottom boundary of . We apply first order absorbing boundary conditions on all four exterior boundaries,

(2)

where n denotes the surface normal and the imaginary unit. Note that first order boundary conditions are not as accurate as certain other boundary conditions, e.g. perfectly matched layers BERENGER_1994, however they are conceptually simpler and simpler to implement. Next, we introduce a design field to control the material distribution in

through the interpolation

function,

(3)

where it is assumed that the background has the value (e.g. air) and where eps_r denotes the relative permittivity of the material used for the structure under design and is a problem dependent scaling factor. The non-physical imaginary term discourages intermediate values of in the design for the focussing problem at hand, by introducing attenuation JENSEN_SIGMUND_2005.

As the baseline example, we consider the design of a monochromatic focussing metalens situated in . To this end, we select the magnitude of at a point in space targetXY as the figure of merit (FOM) denoted , i.e.

(4)

where denotes the complex conjugate.

Figure 1: Model domain, , of height and width with a designable region, , of height and width on top of a substrate, , of height .

The model equation, boundary conditions, material-interpolation function and figure of merit are all discretized using the finite element method (FEM) BOOM_FEM_JIN using nElX nElY bi-linear quadratic elements. The discretized model uses

nodal degrees of freedom (DOFs) for

and as well as elementwise constant DOFs for , with of the elements situated in . The following constrained continuous optimization problem is formulated for the discretized problem,

where and F

are vectors containing the nodal DOFs for the electric field and forcing term and

dVs and are vectors of element DOFs for the design field and the physical filtered and thresholded field (eqs. (6)-(7)), respectively.

The diagonal selection matrix P weights the -DOFs that enter FOM. In the baseline example it has four entries for the selected element’s nodal DOFs. I.e. it selects the field intensity at the focal point, which for simplicity is taken to be at the center of a single finite element. Finally, denotes the conjugate transpose.

To ameliorate numerical issues, such as pixel-by-pixel design variations, and to introduce a weak sense of geometric length scale, a filter and threshold scheme is applied to , before using it to interpolate the material parameters GUEST_ET_AL_2004; WANG_ET_AL_2011; CHRISTIANSEN_SMO_2015. First, the following convolution-based filter operation is applied,

(6)

Here denotes the area of the ’th element and fR the desired spatial filtering radius, finally denotes the ’th set of finite elements who’s center point is within of the ’th element. Second, a smoothed approximation of the Heaviside function is applied to the filtered design variables as,

(7)

Here and control the threshold sharpness and value, respectively.

Adjoint sensitivity analysis TORTORELLI_ET_AL_1994; JENSEN_SIGMUND_2011

is carried out to compute the design sensitivities, utilizing the chain rule for the filter and threshold steps as

Duhring_2008; CHRISTIANSEN_SMO_2015,

(8)

where denotes the real part, the transpose and a vector obtained by solving,

(9)

where denotes the imaginary part. The ’th entry of the right-hand side in eq. (9) is given by,

(10)

The derivations of the expression for in eq. (8) and for eq. (9) are given in Appendix A.

The fundamental advantage of using the adjoint approach to compute the sensitivities of the FOM with respect to the design variables is that only one single additional system of equations, (namely eq. (9)), must be solved. Hence, the sensitivity information is obtained at approximately the same computational cost as the one required to compute the field information itself. In fact, for the examples treated here, it is possible to reuse the LU-factorization used to compute the field information (second line of eq. (II)), making the computational cost associated with computing the sensitivity of the FOM almost ignorable.

Iii The MATLAB Code

The design problem stated in eq. (II) is implemented in top200EM (see the full code in Appendix B), which has the interface:

function [dVs,FOM] = ...
top200EM(targetXY,dVElmIdx,nElX,nElY,dVini,eps_r,lambda,fR,maxItr);

The function takes the input parameters:

  • targetXY: two-values 1D-array with the x- and y-position of the finite element containing the focal point.

  • dVElmIdx: 1D-array of indices of the finite elements which are designable, i.e. .

  • nElX: Number of finite elements in the x-direction.

  • nElY: Number of finite elements in the y-direction.

  • dVini: Initial guess for discretized -field. Accepts a scalar for all elements or a 1D-array of identical length to dVElmIdx.

  • eps_r: Relative permittivity for the material constituting the structure under design.

  • lambda: The targeted wavelength, , measured in number of finite elements.

  • fR: Filter radius, , measured in number of finite elements.

  • maxItr: Maximum number of iterations allowed by fmincon for solving the optimization problem in eq. (II).

And returns the output parameters:

  • dVs: 1D-array of the optimized discretized -field in the design domain, .

  • FOM: Value of the figure of merit.

During execution, the data related to the physics, the discretization and the filter and threshold operations are stored in the structures phy, dis and filThr, respectively.

The spatial scaling, threshold sharpness and threshold level are hard coded in top200EM as,

  • % SETUP OF PHYSICS PARAMETERS
    phy.scale = 1e-9; % Scaling finite element side length to nanometers
    
    % SETUP FILTER AND THRESHOLDING PARAMETERS
    filThr.beta = 5; % Thresholding sharpness
    filThr.eta = 0.5; % Thresholding level
    

Note: For simplicity the code uses the unit of nanometers to measure length and the finite elements are taken to have a side length of 1 nm. This may be changed by changing the scaling parameter phy.scale.

The algorithm used to solve the design problem is MATLABs fmincon,

  • [dVs,~] = fmincon(FOM,dVs(:),[],[],[],[],LBdVs,UBdVs,[],options);
    

with the design-variable bounds and options set up as,

  • LBdVs = zeros(length(dVs),1); % Lower bound on design variables
    UBdVs = ones(length(dVs),1); % Upper bound on design variables
    options = optimoptions(’fmincon’,’Algorithm’,’interior-point’,...
    ’SpecifyObjectiveGradient’,true,’HessianApproximation’,’lbfgs’,...
    ’Display’,’off’,’MaxIterations’,maxItr,’MaxFunctionEvaluations’,maxItr);
    

The FOM and sensitivites provided to fmincon are computed using the inline function:

  • FOM = @(dVs)OBJECTIVE_GRAD(dVs,dis,phy,filThr);
    

The discretized design field is distributed in the model domain with a hard coded background of air in the top 90% of the domain and solid material in the bottom 10% as,

  • % DISTRIBUTE MATERIAL IN MODEL DOMAIN BASED ON DESIGN FIELD
    dFP(1:dis.nElY,1:dis.nElX) = 0; % Design field in physics, 0: air
    dFP(dis.nElY:-1:ceil(dis.nElY*9/10),1:dis.nElX) = 1; % 1: material
    dFP(dis.dVElmIdx(:)) = dVs; % Design variables inserted in design field
    

Followed by the application of the filter and threshold operations and the material interpolation,

  • % COMPUTE MATERIAL FIELD FROM DESIGN FIELD
    dFPS = DENSITY_FILTER(filThr.filKer,filThr.filSca,dFP,ones(dis.nElY,dis.nElX));
    dFPST = THRESHOLD( dFPS, filThr.beta, filThr.eta);
    [A,dAdx] = MATERIAL_INTERPOLATION(phy.eps_r,dFPST,1.0); % Material field
    

The system matrix for the state equation is constructed,

  • % CONSTRUCT SYSTEM MATRIX
    [dis,F] = BOUNDARY_CONDITIONS_RHS(phy.k,dis,phy.scale);
    dis.vS = reshape(dis.LEM(:)-phy.k^2*dis.MEM(:)*(A(:).’),16*dis.nElX*dis.nElY,1);
    SysMat = sparse([dis.iS(:);dis.iBC(:)],[dis.jS(:);dis.jBC(:)],[dis.vS(:);dis.vBC(:)]);
    

The state system is solved using LU-factorization,

  • % SOLVING STATE SYSTEM: SysMat * Ez = F
    [L,U,Q1,Q2] = lu(SysMat); % LU - factorization
    Ez = Q2 * (U\(L\(Q1 * F)));  Ez = full(Ez); % Solving
    

The FOM is computed,

  • % FIGURE OF MERIT
    P = sparse(dis.edofMat(dis.tElmIdx,:),dis.edofMat(dis.tElmIdx,:),1/4,...
    (dis.nElX+1)*(dis.nElY+1),(dis.nElX+1)*(dis.nElY+1)); % Weighting matrix
    FOM = Ez’ * P * Ez; % Solution in target element
    

The adjoint system is solved by reusing the LU-factorization from the state problem,

  • % ADJOINT RIGHT HAND SIDE (0th-order quadrature)
    AdjRHS = P*(2*real(Ez) - 1i*2*imag(Ez));
    
    % SOLVING THE ADJOING SYSTEM: S.’ * AdjLambda = AdjRHS
    AdjLambda = (Q1.’) * ((L.’)\((U.’)\((Q2.’) * (-1/2*AdjRHS)))); % Solving
    

The sensitivities in are computed and filtered, and the values in extracted,

  • % SENSITIVITIES
    dis.vDS = reshape(-phy.k^2*dis.MEM(:)*(dAdx(:).’),16*dis.nElX*dis.nElY,1);
    DSdx = sparse(dis.iElFull,dis.jElFull,dis.vDS); % Constructing dS/dx
    DSdxMulV = DSdx * Ez(dis.idxDSdx); % Computing dS/dx * Field values
    DsdxMulV = sparse(dis.iElSens,dis.jElSens,DSdxMulV);
    sens = 2*real(AdjF(dis.idxDSdx).’ * DsdxMulV); % Computing sensitivites
    sens = full(reshape(sens,dis.nElY,dis.nElX));
    
    % FILTERING SENSITIVITIES
    DdFSTDFS = DERIVATIVE_OF_THRESHOLD( dFPS, filThr.beta, filThr.eta);
    sensFOM = DENSITY_FILTER(filThr.filKer,filThr.filSca,sens,DdFSTDFS);
    
    % EXTRACTING SENSITIVITIES FOR DESIGN REGION
    sensFOM = sensFOM(dis.dVElmIdx);
    

Finally, and are plotted and the current FOM-value printed,

  • % PLOTTING AND PRINTING
    figure(1); % Field intensity, |Ez|^2
    imagesc((reshape(Ez.*conj(Ez),dis.nElY+1,dis.nElX+1))); colorbar; axis equal;
    figure(2); % Physical design field
    imagesc(1-dFPST); colormap(gray); axis equal; drawnow;
    disp([’FOM: ’ num2str(-FOM)]); % Display FOM value
    

After the TopOpt procedure is finished, a thresholded version of the final design is evaluated and the resulting field and design field are plotted.

  • % FINAL BINARIZED DESIGN EVALUATION filThr.beta = 1000; disp(’Black/white design evaluation:’) [obj_2,dFPST_2,F_2] = OBJECTIVE_GRAD(DVini(:),dis,phy,filThr);

Iv Using the Code

Next, we demonstrate how to use top200EM by designing a focusing metalens as follows.

iv.1 Designing a Metalens

First, we define the domain size in terms of the number of finite elements in each spatial direction and the element indices for the design domain.

  • % DESIGN FIELD INDICES
    DomainElementsX = 400;
    DomainElementsY = 200;
    DesignThicknessElements = 15;
    DDIdx = repmat([1:DomainElementsY:DomainElementsX*DomainElementsY],...
    DesignThicknessElements,1);
    DDIdx = DDIdx+repmat([165:165+DesignThicknessElements-1]’,1,DomainElementsX);
    

Second, the optimization problem is solved by executing the command,

  • [DVs,obj]=top200EM([200,80],DDIdx,DomainElementsX,DomainElementsY,...
                       0.5,3.0,35,6.0,200);
    
Figure 2: (a) Max-normalized -field in . (b) Metalens design, (black) and (white).

The final binarized design is shown in Fig. 2b with black(white) representing solid(air). Figure 2a shows the -field resulting from exciting the metalens in Fig. 2b for the targeted incident field, demonstrating the focussing effect at the targeted focal spot. The numerical aperture of the metalens is NA and the transmission efficiency is computed as the power propagating through the lens relative to the power incident on the lens.

iv.2 Solving the same design problem with different resolutions

For various reasons, such as performing a mesh-convergence study, it may be necessary to solve the same physical model problem using different mesh resolutions. This may be done with top200EM by multiplying the following inputs by an integer scaling factor: The number of finite elements in each spatial direction, nElX and nElY, the wavelength lambda and the filter radius fR and dividing the hard-coded value of phys.scale in the code by the same factor.

iv.3 The Effect of Filtering

Next, we demonstrate the effect of applying the filtering step LAZAROV_ET_AL_2016 by changing the filter radius and designing four metalenses (See Fig. 3) using top200EM. Again, the model-domain size and indices for the design domain are defined first,

  • % DESIGN FIELD INDICES
    DomainElementsX = 400;
    DomainElementsY = 200;
    DesignThicknessElements = 15;
    DDIdx = repmat([1:DomainElementsY:DomainElementsX*DomainElementsY],...
    DesignThicknessElements,1);
    DDIdx = DDIdx+repmat([165:165+DesignThicknessElements-1]’,1,DomainElementsX);
    

Then, the optimization problem is solved with the four filtering radii fR,

  • [DVs,obj]=top200EM([200,80],DDIdx,DomainElementsX,DomainElementsY,...
    0.5,3.0,35,fR,200);
    
Figure 3: Metalens designs obtained (a) without filtering , (b) using a filter radius of , (c) using a filter radius of , (d) using a filter radius of . These results illustrate the effect of applying the cone-shaped filter to the design field as part of the optimization process.

Looking at the four final binarized designs in Fig. 3, it is clearly observed that as is increased the features in the designs grow. When no filtering is applied (Fig. 3a) single-pixel sized features are observed. Such features may be problematic from a numerical-modelling point of view as well as being detrimental to fabrication.

It is noted that without filtering TopOpt is more prone to identify a local minimum, which performs worse than ones identified with filtering. In other words, filtering tends to have a convexifying effect, as long as the filter size is not too big.

V Modifying the Code

There exists a vast amount of auxiliary tools developed to extend the applicability of density-based Topology Optimization across a wide range of different problems and physics. The following sections provide examples of how simple some of these tools are to implement in top200EM.

v.1 Plasmonics

In recent work CHRISTIANSEN_VESTER_2019 it was demonstrated that using a refractive index and extinction cross section based non-linear material interpolation yielded significantly improved results when designing Au, Ag and Cu nano-particles for localized field enhancement, with recent applications to enhanced upconversion CHRISTIANSEN_SEMSC_2020. thermal emission KUDYSHEV_ET_AL_2020 and Raman scattering CHRISTIANSEN_OE_2020. This interpolation function avoids artificial resonances in connection with the transition from positive to negative values and reads,

Here and denote the refractive index and extinction cross section, respectively. The subscripts and denote the two materials being interpolated.

The interpolation scheme is straight forward to implement in the code as follows.

First the following lines of code,

  • function [A,dAdx] = MATERIAL_INTERPOLATION(eps_r,x,1.0)
    A = 1 + x*(eps_r-1) - 1i * alpha_i * x .* (1 - x); % Interpolation
    dAdx = (eps_r-1)*(1+0*x) - 1i * alpha_i * (1 - 2*x); % Derivative of interpolation
    end
    

are replaced with,

  • function [A,dAdx] = MATERIAL_INTERPOLATION(n_r,k_r,x)
    n_eff = 1 + x*(n_r-1);
    k_eff = 0 + x*(k_r-0);
    A = (n_eff.^2 - k_eff.^2) - 1i*(2.*n_eff.*k_eff);
    dAdx = 2*n_eff*(n_r-1)-2*k_eff*(k_r-1)-1i*(2*(n_r-1)*k_eff+(2*n_eff*(k_r-1)));
    end
    

where for simplicity it is assumed that is air, i.e. and .

Second, the scalar input parameter eps_r is changed to a two-valued 1D-array nk_r.

Third, the line,

  • phy.eps_r = eps_r; % Relative permittivity
    

is replaced with,

  • phy.nk_r = nk_r; % Refracitve index and extrinction coefficient
    

and fourth, the call to the material interpolation function is changed from,

  • [A,dAdx] = MATERIAL_INTERPOLATION(phy.eps_r,dFPST,1.0);
    Ψ
    

to,

  • [A,dAdx] = MATERIAL_INTERPOLATION(phy.nk_r(1),phy.nk_r(2),dFPST);
    Ψ
    

v.2 Excitation

The excitation considered in the baseline problem may be changed straightforwardly. As an example, the boundary at which the incident field enters the domain may be changed from the bottom to the top boundary as follows.

First, the index set controlling where the boundary condition is imposed in the vector, F(=F), is changed by moving dis.iRHS = TMP; from line 157 to above line 151.

Second, the values stored in F are changed to account for the propagation direction of the wave by replacing,

  • F(dis.iRHS(1,:)) = F(dis.iRHS(1,:))+1i*waveVector;
    F(dis.iRHS(2,:)) = F(dis.iRHS(2,:))+1i*waveVector;
    

with,

  • F(dis.iRHS(1,:)) = F(dis.iRHS(1,:))-1i*waveVector;
    F(dis.iRHS(2,:)) = F(dis.iRHS(2,:))-1i*waveVector;
    

v.3 Designing a Metallic Reflector

By introducing the changes presented in Sec. V.1 and Sec. V.2, one may design a metallic reflector using topEM200.

Figure 4: (a) Max-normalized -field in . (b) Metallic reflector design (black) in air background (white).

First, we set the domain dimensions and design thickness as in the previous examples,

  • % DESIGN FIELD INDICES
    DomainElementsX = 400;
    DomainElementsY = 200;
    DesignThicknessElements = 15;
    DDIdx = repmat([1:DomainElementsY:DomainElementsX*DomainElementsY],...
    DesignThicknessElements,1);
    DDIdx = DDIdx+repmat([165:165+DesignThicknessElements-1]’,1,DomainElementsX);
    Ψ
    

Second, we solve the optimization problem by executing the command,

  • [DVs,obj]=top200EM([200,100],DDIdx,DomainElementsX,DomainElementsY,...
    0.5,[1.9,1.5],35,3.0,200);
    Ψ
    

Note: For simplicity we selected the following values for the refractive index, (=nk_r(1)), and extinction cross section, (=nk_r(2)).111These values corresponds to values for gold at nm. I.e. one may think of this choice of material parameters as an unspoken rescaling of space by a factor of 10 i.e. changing the element size to pixels of 10 nm by 10 nm rather than of 1 nm by 1 nm.

Considering the final binarized reflector design in Fig. 4b, one can interpret the design as the well-known parabolic reflector broken into pieces to fit the spatially-limited design domain. From the max-normalized -field presented in Fig. 4a the focussing effect of the reflector is clearly observed.

v.4 Linking Design Variables

Certain fabrication techniques limit the allowable geometric variations in a design. For example, optical-projection lithography and electron-beam lithography restricts variations in design geometries to two-dimensional patterns, which can then be extruded in the out-of-plane direction. It is straight forward to introduce such a geometric restriction using TopOpt by linking design variables and sensitivities across elements. In top200EM the design field may be restricted to only exhibit in-plane (x-direction) variations as follows.

First, the code-line representing the values of the design variables,

  • dVs(length(dis.dVElmIdx(:))) = dVini; % Design variables
    

is replaced with,

  • dVs(1:nElX) = dVini; % Design variables for 1D design
    

Second, the code-line transferring the design variables to the elements in the physics model,

  • dFP(dis.dVElmIdx(:)) = dVs; % Design variables inserted in design field
    

is replaced with,

  • nRows=length(dis.dVElmIdx(:))/dis.nElX; % Number of rows in the 1D design
    dFP(dis.dVElmIdx) = repmat(dVs,1,nRows)’; % Design variables inserted in design field
    

and finally the code-line representing individual element sensitivities,

  • sensFOM = sensFOM(dis.dVElmIdx);
    

is replaced with,

  • sensFOM = sensFOM(dis.dVElmIdx);
    sensFOM = sum(sensFOM,1); % Correcting sensitivities for 1D design
    

which sums the sensitivity contributions from the linked elements.

v.5 Designing a 1D metalens

By introducing the changes listed in Sec. V.4 topEM200 can be used to design metasurfaces of fixed height with in-plane variations as follows. Again we define the domain dimensions as,

  • % DESIGN FIELD INDICES
    DomainElementsX = 400;
    DomainElementsY = 200;
    DesignThicknessElements = 15;
    DDIdx = repmat([1:DomainElementsY:DomainElementsX*DomainElementsY],...
    DesignThicknessElements,1);
    DDIdx = DDIdx+repmat([165:165+DesignThicknessElements-1]’,1,DomainElementsX);
    Ψ
    

Followed by the execution of the command,

  • [DVs,obj]=top200EM([200,80],DDIdx,DomainElementsX,DomainElementsY,...
    0.5,3.0,35,3.0,200);
    Ψ
    
Figure 5: (a) Max-normalized -field in . (b) Metalens design restricted to one dimensional variations, (black) and (white).

The final binarized design, resulting from solving the design problem, is shown in Fig. 5b. Here it is clear to see that the design is now restricted to vary only in the x-direction. It is worth noting that the design is still filtered in both spatial directions, hence the corners of the design appear rounded in Fig. 5b. The max-normalized -field presented in Fig. 4a demonstrates the focussing effect of the lens at the targeted point in the modelling domain, which, due to the reduced design freedom, is not as high as in the un-restricted case (Sec. IV.1).

v.6 Continuation of Threshold Sharpness

For some design problems in photonics and plasmonics, e.g. CHRISTIANSEN_VESTER_2019; CHRISTIANSEN_OE_2020; WANG_ET_AL_2011, intermediate values may be present in the final optimized design, i.e. , despite the applied penalization scheme, as they prove beneficial to optimizing the FOM. However, (in most cases) intermediate values hold no physical meaning and it is therefore not possible to realize designs containing such intermediate values experimentally. A way to promote that (almost) no design variables take intermediate values in the final design, is by using a continuation scheme for the threshold sharpness, gradually increasing it until an (almost) pure 0/1-design is achieved, ensuring that the optimized designs are physically realizable.

To implement the continuation scheme in the code, replace the following lines,

  • Ψ% SOLVE DESIGN PROBLEM USING MATLAB BUILT-IN OPTIMIZER: FMINCON
    ΨFOM = @(dVs)OBJECTIVE_GRAD(dVs,dis,phy,filThr);
    Ψ[dVs,~] = fmincon(FOM,dVs(:),[],[],[],[],LBdVs,UBdVs,[],options);
    Ψ
    

with,

  • Ψwhile filThr.beta<betaMax % Thresholding sharpness bound
    ΨΨ% SOLVE DESIGN PROBLEM USING MATLAB BUILT-IN OPTIMIZER: FMINCON
    ΨΨFOM = @(dVs)OBJECTIVE_GRAD(dVs,dis,phy,filThr);
    ΨΨ[dVs,~] = fmincon(FOM,dVs(:),[],[],[],[],LBdVs,UBdVs,[],options);
    ΨΨfilThr.beta = betaInc * filThr.beta; % Increasing thresholding sharpness
    Ψend
    Ψ
    

and select suitable values for betaMax and betaInc. These values are problem dependent and some experimentation may be required to identify the best values for a given problem. For the metalens design example in Sec. IV.1, when considering high material contrast, i.e. large values of eps_r, the values betaMax = 20.0 and betaInc = 1.5 have been found to work well.

Vi What about genetic algorithms?

For a wide range of inverse design problems, where gradient-based optimization methods are applicable222Inverse design problems where it is possible to compute the sensitivities of the FOM (and any constraints), e.g. using adjoint sensitivity analysis, a pre-requisite for using gradient-based methods., they have been found to severely outperform non-gradient-based methods, such as genetic algorithms (GA) GOLDBERG_1989. This, both in terms of the computational effort required to identify a local optimum for the FOM and, by direct extension, in terms of the number of design degrees of freedom it is feasible to consider (a difference of many orders of magnitude AAGE_ET_AL2017). Further, in many cases, gradient-based methods are able to identify better local optima for the FOM SIGMUND_2011.

In the following, we provide a demonstration of these claims, by comparing the solution of a metalens design problem obtained using the gradient-based top200EM to the solution obtained using MATLABs built in genetic algorithm ga.

Readers may perform their own comparisons by rewriting top200EM to utilize ga instead of fmincon as follows.

First, a (near perfect) 0/1-design is ensured by changing the thresholding strength from,

  • ΨfilThr.beta = 5; % Thresholding sharpness
    Ψ
    

to,

  • ΨfilThr.beta = 1e5; % Thresholding sharpness
    Ψ
    

Second, ga is used instead of fmincon by replacing the following lines of codes,

  • Ψoptions = optimoptions(’fmincon’,’Algorithm’,’interior-point’,...
    Ψ’SpecifyObjectiveGradient’,true,’HessianApproximation’,’lbfgs’,...
    Ψ’Display’,’off’,’MaxIterations’,maxItr,’MaxFunctionEvaluations’,maxItr);
    
    Ψ% SOLVE DESIGN PROBLEM USING MATLAB BUILT-IN OPTIMIZER: FMINCON
    ΨFOM = @(dVs)OBJECTIVE_GRAD(dVs,dis,phy,filThr);
    Ψ[dVs,~] = fmincon(FOM,dVs(:),[],[],[],[],LBdVs,UBdVs,[],options);
    Ψ
    

with,

  • Ψoptions = optimoptions(’ga’,’MaxGenerations’,maxItr,’Display’,’iter’);
    
    Ψ% SOLVE DESIGN PROBLEM USING MATLAB BUILT-IN OPTIMIZER: GA
    ΨFOM = @(dVs)OBJECTIVE(dVs,dis,phy,filThr);
    Ψrng(1,’twister’); % SETTING RANDOM SEED TO ENSURE REPRODUCTION OF RESULT
    Ψ[dVs,~,~,~] = ga(FOM,length(dVs),[],[],[],[],LBdVs,UBdVs,[],[],options);
    Ψ
    

where the command rng(1,’twister’) fixes the random seed for reproducibility.

Third, the computation of the gradient of the FOM is removed by changing,

  • Ψfunction [FOM,sensFOM] = OBJECTIVE_GRAD(dVs,dis,phy,filThr)
    Ψ
    

to,

  • Ψfunction [FOM] = OBJECTIVE(dVs,dis,phy,filThr)
    Ψ
    

Fourth, deleting the following lines in OBJECTIVE,

  • Ψ% ADJOINT RIGHT HAND SIDE
    ΨAdjRHS = P*(2*real(Ez) - 1i*2*imag(Ez));
    
    Ψ% SOLVING THE ADJOING SYSTEM: S.’ * AdjLambda = AdjRHS
    ΨAdjLambda = (Q1.’) * ((L.’)\((U.’)\((Q2.’) * (-1/2*AdjRHS)))); % Solving
    
    Ψ% COMPUTING SENSITIVITIES
    Ψdis.vDS = reshape(-phy.k^2*dis.MEM(:)*(dAdx(:).’),16*dis.nElX*dis.nElY,1);
    ΨDSdx = sparse(dis.iElFull,dis.jElFull,dis.vDS); % Constructing dS/dx
    ΨDSdxMulV = DSdx * Ez(dis.idxDSdx); % Computing dS/dx * Field values
    ΨDsdxMulV = sparse(dis.iElSens,dis.jElSens,DSdxMulV);
    Ψsens = 2*real(AdjLambda(dis.idxDSdx).’ * DsdxMulV); % Computing sensitivites
    Ψsens = full(reshape(sens,dis.nElY,dis.nElX));
    
    Ψ% FILTERING SENSITIVITIES
    ΨDdFSTDFS = DERIVATIVE_OF_THRESHOLD( dFPS, filThr.beta, filThr.eta);
    ΨsensFOM = DENSITY_FILTER(filThr.filKer,filThr.filSca,sens,DdFSTDFS);
    
    Ψ% EXTRACTING SENSITIVITIES FOR DESIGNABLE REGION
    ΨsensFOM = sensFOM(dis.dVElmIdx);
    
    Ψ% FMINCON DOES MINIMIZATION
    ΨsensFOM = -sensFOM(:);
    
    Ψ
    

and finally changing the code evaluating the final design,

  • Ψ% FINAL BINARIZED DESIGN EVALUATION
    ΨfilThr.beta = 1000;
    Ψdisp(’Black/white design evaluation:’);
    Ψ[FOM,~] = OBJECTIVE(dVs(:),dis,phy,filThr);
    Ψ
    

to,

  • Ψ% FINAL BINARIZED DESIGN EVALUATION
    ΨfilThr.beta = 1e5;
    Ψdisp(’Black/white design evaluation:’);
    Ψ[FOM,~] = OBJECTIVE(dVs(:),dis,phy,filThr);
    Ψ
    

Note: The code-lines plotting the design in OBJECTIVE can be removed to reduce wallclock time.

In order to reduce the computational effort required to reproduce the following example, we consider a problem in a smaller spatial domain, considering a shorter wavelength and fewer degrees of freedom than in the previous examples. This is done by defining a model problem considering 1000 design degrees of freedom as,

  • Ψ% DESIGN FIELD INDICES
    ΨDomainElementsX = 100;
    ΨDomainElementsY = 50;
    ΨDesignThicknessElements = 10;
    ΨDDIdx = repmat([1:DomainElementsY:DomainElementsX*DomainElementsY],...
    ΨDesignThicknessElements,1);
    ΨDDIdx = DDIdx+repmat([35:35+DesignThicknessElements-1]’,1,DomainElementsX);
    Ψ
    

Followed by the solution of the design problem using the modified GA-based code,

  • Ψ[dVs_GA,FOM_GA]=top200EMGA([50,10],DDIdx,DomainElementsX,DomainElementsY,...
    Ψ0.5,3.0,20,3.0,500);
    Ψ
    

and then using the original gradient-based code,

  • Ψ[dVs,FOM]=top200EM([50,10],DDIdx,DomainElementsX,DomainElementsY,...
    Ψ0.5,3.0,20,3.0,500);
    Ψ
    

The binarized designs obtained using top200EMGA and top200EM are presented in Fig. 6a and Fig. 6b, respectively. The evolution of the FOM value as function of the number of FOM evaluations is plotted in Fig. 6c for both methods.

Figure 6: (a) Design obtained using the GA-based method. (b) Design obtained using the gradient-based method. (black) and (white). (c) Convergence graph. Remark that each design iteration for the GA-based code requires 200 solutions of the physics model equation, while it requires 1 solution of the physics model equation and 1 solution of the adjoint equation for the gradient-based code!

The GA-based method uses 137 design iterations (27,600 FOM evaluations) to identify a solution. In contrast the gradient-based method uses 208 iterations (208 FOM evaluation) to identify a superior solution. Evaluating the final binarized designs the GA-based solution has a FOM value of while the gradient-based solution has an better FOM value of .

Crucially, the GA-based method uses 200 FOM evaluations per design iteration to drive the optimization process. This results in a total of 27.600 FOM evaluations for solving the design problem using the GA-based method, which corresponds to 27.600 solutions of the physics model equation, i.e. eq. (1). In contrast, when solving the design problem using the gradient-based method a total of 208 FOM evaulations are performed, corresponding to 208 solution of the physics model equation plus 208 solutions of the adjoint system (eq. (9)), which for the baseline example is almost free due to the reuse of the LU-factorization.

Hence, in our example the identification of a local optimum for the FOM using the gradient-based method required of the computational effort spend by the GA-based method. A ratio, which only becomes more pronounced as the size of the design space grows.

This example is based on standard settings for MATLAB’s ga optimizer and can surely be improved and refined. However, this will not change the basic conclusion that non-gradient-based approaches are unsuitable for large-scale TopOpt problems. A more thorough discussion of non-gradient-based methods in TopOpt can be found in Ref. SIGMUND_2011. Readers are invited to form their own conclusions by extending the example and codes provided here.

Vii Conclusion

We have presented a simple finite element based MATLAB code (downloadable from https://www.topopt.dtu.dk) for TopOpt-based inverse design of photonic structures. We have provided examples of how to use the code, as well as a set of suggestions for code extensions enabling it to handle metallic structures, different model excitations, linked design variables and introducing a continuation scheme for the threshold sharpness designed to promote 0/1-designs. Finally, we have demonstrated the superiority of a gradient-based method over a genetic-algorithm-based method when performing inverse design in photonics.

The code can be used for educational purposes as is, and is otherwise meant to serve as a starting point for the reader to develop software to handle their more advanced research applications within photonics. For simplicity and computational speed, the code treats problems in two spatial dimensions, however it is directly extendable to three spatial dimensions by modifying the finite element matrices, boundary conditions and index sets appropriately.

Appendix A Adjoint Sensitivity Analysis

The expression for in eq. (8) and the expression in eq. (9) may be derived as follows. First, zero is added to twice,

(11)

where and is a vector of nodal complex Lagrange multipliers, also called the adjoint variables. Second, one takes the derivative of with respect to and exploits that, for the optimization problem in eq. (II), does not depend explicitly on and neither nor F depend on at all, yielding,

(12)

Collecting terms including and and reducing the remaining terms yields,

(13)

To eliminate the first two terms in eq. (13) containing the derivates and , the two parenthesis must equal zero,

(14)

multiplying the second equation by i, subtracting it from the first and transposing it yields,

(15)

Using eq. (15) the expression in eq. (13) reduces to the expression for in eq. (8) and the second equation in eq. (15) is equal to eq. (9).

Appendix B MATLAB code

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% A 200 LINE TOPOLOGY OPTIMIZATION CODE FOR ELECTROMAGNETISM %%%%%%%%
% --------------------------- EXAMPLE GOAL ------------------------------ %
% Designs a 2D metalens with relative permittivity eps_r capable of       %
% monocromatic focusing of TE-polarized light at a point in space.        %
% --------------------- FIGURE OF MERIT MAXIMIZED ----------------------- %
%      Phi = |Ez|^2 in a "point" (in the center of a finite element)      %
% ------------------------- EQUATION SOLVED ----------------------------- %
%      \nabla * (\nabla Ez) + k^2 A Ez = F                                %
%      With first order absorping boundary condition:                     %
%      n * \nabla Ez = - i k Ez on boundaries                             %
%      and an incident plane wave propagating from bottom to top          %
% ---------------------- DOMAIN AND DISCRETIZATION ---------------------- %
% The equation is solved in a rectangular domain, discretized using       %
% quadrilateral bi-linear finite elements                                 %
% ----------------------------------------------------------------------- %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% Author: Rasmus E. Christansen, June 2020 %%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Disclaimer:                                                             %
% The authors reserves all rights but does not guaranty that the code is  %
% free from errors. Furthermore, we shall not be liable in any event      %
% caused by the use of the program.                                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [dVs,FOM]=top200EM(targetXY,dVElmIdx,nElX,nElY,dVini,eps_r,lambda,fR,maxItr)
% SETUP OF PHYSICS PARAMETERS
phy.scale = 1e-9; % Scaling finite element side length to nanometers
phy.eps_r = eps_r; % Relative permittivity
phy.k = 2*pi/(lambda*phy.scale); % Free-space wavenumber

% SETUP OF ALL INDEX SETS, ELEMENT MATRICES AND RELATED QUANTITIES
dis.nElX = nElX; % number of elements in x direction
dis.nElY = nElY; % number of elements in y direction
dis.tElmIdx = (targetXY(1)-1)*nElY+targetXY(2); % target index
dis.dVElmIdx = dVElmIdx; % design field element indices in model of physics
[dis.LEM,dis.MEM] = ELEMENT_MATRICES(phy.scale);
[dis]=INDEX_SETS_SPARSE(dis); % Index sets for discretized model

% SETUP FILTER AND THRESHOLDING PARAMETERS
filThr.beta = 5; % Thresholding sharpness
filThr.eta = 0.5; % Thresholding level
[filThr.filKer, filThr.filSca] = DENSITY_FILTER_SETUP( fR, nElX, nElY);

% INITIALIZE DESIGN VARIABLES, BOUNDS AND OPTIMIZER OPTIONS
dVs(length(dis.dVElmIdx(:))) = dVini; % Design variables
LBdVs = zeros(length(dVs),1); % Lower bound on design variables
UBdVs = ones(length(dVs),1); % Upper bound on design variables
options = optimoptions(’fmincon’,’Algorithm’,’interior-point’,...
’SpecifyObjectiveGradient’,true,’HessianApproximation’,’lbfgs’,...
’Display’,’off’,’MaxIterations’,maxItr,’MaxFunctionEvaluations’,maxItr);

% SOLVE DESIGN PROBLEM USING MATLAB BUILT-IN OPTIMIZER: FMINCON
FOM = @(dVs)OBJECTIVE_GRAD(dVs,dis,phy,filThr);
[dVs,~] = fmincon(FOM,dVs(:),[],[],[],[],LBdVs,UBdVs,[],options);

% FINAL BINARIZED DESIGN EVALUATION
filThr.beta = 1000;
disp(’Black/white design evaluation:’)
[FOM,~] = OBJECTIVE_GRAD(dVs(:),dis,phy,filThr);
end

%%%%%%%%%%%%%%% OBJECTIVE FUNCTION AND GRADIENT EVALUATION %%%%%%%%%%%%%%%%
function [FOM,sensFOM] = OBJECTIVE_GRAD(dVs,dis,phy,filThr)
% DISTRIBUTE MATERIAL IN MODEL DOMAIN BASED ON DESIGN FIELD
dFP(1:dis.nElY,1:dis.nElX) = 0; % Design field in physics, 0: air
dFP(dis.nElY:-1:ceil(dis.nElY*9/10),1:dis.nElX) = 1; % 1: material
dFP(dis.dVElmIdx(:)) = dVs; % Design variables inserted in design field

% FILTERING THE DESIGN FIELD AND COMPUTE THE MATERIAL FIELD
dFPS = DENSITY_FILTER(filThr.filKer,filThr.filSca,dFP,ones(dis.nElY,dis.nElX));
dFPST = THRESHOLD( dFPS, filThr.beta, filThr.eta);
[A,dAdx] = MATERIAL_INTERPOLATION(phy.eps_r,dFPST,1.0); % Material field

% CONSTRUCT THE SYSTEM MATRIX
[dis,F] = BOUNDARY_CONDITIONS_RHS(phy.k,dis,phy.scale);
dis.vS = reshape(dis.LEM(:)-phy.k^2*dis.MEM(:)*(A(:).’),16*dis.nElX*dis.nElY,1);
S = sparse([dis.iS(:);dis.iBC(:)],[dis.jS(:);dis.jBC(:)],[dis.vS(:);dis.vBC(:)]);

tic;
% SOLVING THE STATE SYSTEM: S * Ez = F
[L,U,Q1,Q2] = lu(S); % LU - factorization
Ez = Q2 * (U\(L\(Q1 * F)));  Ez = full(Ez); % Solving
toc;

% FIGURE OF MERIT
P = sparse(dis.edofMat(dis.tElmIdx,:),dis.edofMat(dis.tElmIdx,:),1/4,...
(dis.nElX+1)*(dis.nElY+1),(dis.nElX+1)*(dis.nElY+1)); % Weighting matrix
FOM = Ez’ * P * Ez; % Solution in target element

% ADJOINT RIGHT HAND SIDE
AdjRHS = P*(2*real(Ez) - 1i*2*imag(Ez));

% SOLVING THE ADJOING SYSTEM: S.’ * AdjLambda = AdjRHS
AdjLambda = (Q1.’) * ((L.’)\((U.’)\((Q2.’) * (-1/2*AdjRHS)))); % Solving

% COMPUTING SENSITIVITIES
dis.vDS = reshape(-phy.k^2*dis.MEM(:)*(dAdx(:).’),16*dis.nElX*dis.nElY,1);
DSdx = sparse(dis.iElFull(:),dis.jElFull(:),dis.vDS(:)); % Constructing dS/dx
DSdxMulV = DSdx * Ez(dis.idxDSdx); % Computing dS/dx * Field values
DsdxMulV = sparse(dis.iElSens,dis.jElSens,DSdxMulV);
sens = 2*real(AdjLambda(dis.idxDSdx).’ * DsdxMulV); % Computing sensitivites
sens = full(reshape(sens,dis.nElY,dis.nElX));

% FILTERING SENSITIVITIES
DdFSTDFS = DERIVATIVE_OF_THRESHOLD( dFPS, filThr.beta, filThr.eta);
sensFOM = DENSITY_FILTER(filThr.filKer,filThr.filSca,sens,DdFSTDFS);

% EXTRACTING SENSITIVITIES FOR DESIGNABLE REGION
sensFOM = sensFOM(dis.dVElmIdx);

% FMINCON DOES MINIMIZATION
FOM = -FOM; sensFOM = -sensFOM(:);

% PLOTTING AND PRINTING
figure(1); % Field intensity, |Ez|^2
imagesc((reshape(Ez.*conj(Ez),dis.nElY+1,dis.nElX+1))); colorbar; axis equal;
figure(2); % Physical design field
imagesc(1-dFPST); colormap(gray); axis equal; drawnow;
disp([’FOM: ’ num2str(-FOM)]); % Display FOM value
end

%%%%%%%%%%%%%%%%%%%%%%%%%% AUXILIARY FUNCTIONS %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% ABSORBING BOUNDARY CONDITIONS AND RIGHT HAND SIDE %%%%%%%%%%%%
function [dis,F] = BOUNDARY_CONDITIONS_RHS(waveVector,dis,scaling)
AbsBCMatEdgeValues = 1i*waveVector*scaling*[1/6 ; 1/6 ; 1/3 ; 1/3];
% ALL BOUNDARIES HAVE ABSORBING BOUNDARY CONDITIONS
dis.iBC = [dis.iB1(:);dis.iB2(:);dis.iB3(:);dis.iB4(:)];
dis.jBC = [dis.jB1(:);dis.jB2(:);dis.jB3(:);dis.jB4(:)];
dis.vBC = repmat(AbsBCMatEdgeValues,2*(dis.nElX+dis.nElY),1);
% BOTTOM BOUNDARY HAS INCIDENT PLANE WAVE
F = zeros((dis.nElX+1)*(dis.nElY+1),1); % System right hand side
F(dis.iRHS(1,:)) = F(dis.iRHS(1,:))+1i*waveVector;
F(dis.iRHS(2,:)) = F(dis.iRHS(2,:))+1i*waveVector;
F = scaling*F;
end
%%%%%%%%%%%%%%%%%%%%% CONNECTIVITY AND INDEX SETS %%%%%%%%%%%%%%%%%%%%%%%%%
function [dis]=INDEX_SETS_SPARSE(dis)
% INDEX SETS FOR SYSTEM MATRIX
nEX = dis.nElX; nEY = dis.nElY; % Extracting number of elements
nodenrs = reshape(1:(1+nEX)*(1+nEY),1+nEY,1+nEX); % Node numbering
edofVec = reshape(nodenrs(1:end-1,1:end-1)+1,nEX*nEY,1); % First DOF in element
dis.edofMat = repmat(edofVec,1,4)+repmat([0 nEY+[1 0] -1],nEX*nEY,1);
dis.iS = reshape(kron(dis.edofMat,ones(4,1))’,16*nEX*nEY,1);
dis.jS = reshape(kron(dis.edofMat,ones(1,4))’,16*nEX*nEY,1);
dis.idxDSdx = reshape(dis.edofMat’,1,4*nEX*nEY);
% INDEX SETS FOR BOUNDARY CONDITIONS
TMP = repmat([[1:nEY];[2:nEY+1]],2,1);
dis.iB1 = reshape(TMP,4*nEY,1); % Row indices
dis.jB1 = reshape([TMP(2,:);TMP(1,:);TMP(3,:);TMP(4,:)],4*nEY,1); % Column indices
TMP = repmat([1:(nEY+1):(nEY+1)*nEX;(nEY+1)+1:(nEY+1):(nEY+1)*nEX+1],2,1);
dis.iB2 = reshape(TMP,4*nEX,1);
dis.jB2 = reshape([TMP(2,:);TMP(1,:);TMP(3,:);TMP(4,:)],4*nEX,1);
TMP = repmat([(nEY+1)*(nEX)+1:(nEY+1)*(nEX+1)-1;(nEY+1)*(nEX)+2:(nEY+1)*(nEX+1)],2,1);
dis.iB3 = reshape(TMP,4*nEY,1);
dis.jB3 = reshape([TMP(2,:);TMP(1,:);TMP(3,:);TMP(4,:)],4*nEY,1);
TMP = repmat([2*(nEY+1):nEY+1:(nEY+1)*(nEX+1);(nEY+1):nEY+1:(nEY+1)*(nEX)],2,1);
dis.iB4 = reshape(TMP,4*nEX,1);
dis.jB4 = reshape([TMP(2,:);TMP(1,:);TMP(3,:);TMP(4,:)],4*nEX,1);
dis.iRHS = TMP;
% INDEX SETS FOR INTEGRATION OF ALL ELEMENTS
ima0 = repmat([1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4],1,nEX*nEY).’;
jma0 = repmat([1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4],1,nEX*nEY).’;
addTMP = repmat(4*[0:nEX*nEY-1],16,1);
addTMP = addTMP(:);
dis.iElFull = ima0+addTMP;
dis.jElFull = jma0+addTMP;
% INDEX SETS FOR SENSITIVITY COMPUTATIONS
dis.iElSens = [1:4*nEX*nEY]’;
jElSens = repmat([1:nEX*nEY],4,1);
dis.jElSens = jElSens(:);
end
%%%%%%%%%%%%%%%%%% MATERIAL PARAMETER INTERPOLATION %%%%%%%%%%%%%%%%%%%%%%%
function [A,dAdx] = MATERIAL_INTERPOLATION(eps_r,x,alpha_i)
A = 1 + x*(eps_r-1) - 1i * alpha_i * x .* (1 - x); % Interpolation
dAdx = (eps_r-1)*(1+0*x) - 1i * alpha_i * (1 - 2*x); % Derivative of interpolation
end
%%%%%%%%%%%%%%%%%%%%%%%%%%% DENSITY FILTER  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [xS]=DENSITY_FILTER(FilterKernel,FilterScaling,x,func)
xS = conv2((x .* func)./FilterScaling,FilterKernel,’same’);
end
function [ Kernel, Scaling ] = DENSITY_FILTER_SETUP( fR, nElX, nElY )
[dy,dx] = meshgrid(-ceil(fR)+1:ceil(fR)-1,-ceil(fR)+1:ceil(fR)-1);
Kernel = max(0,fR-sqrt(dx.^2+dy.^2)); % Cone filter kernel
Scaling = conv2(ones(nElY,nElX),Kernel,’same’); % Filter scaling
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%% THRESHOLDING %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ xOut ] = THRESHOLD( xIn, beta, eta)
xOut = (tanh(beta*eta)+tanh(beta*(xIn-eta)))./(tanh(beta*eta)+tanh(beta*(1-eta)));
end
function [ xOut ] = DERIVATIVE_OF_THRESHOLD( xIn, beta, eta)
xOut = (1-tanh(beta*(xIn-eta)).^2)*beta./(tanh(beta*eta)+tanh(beta*(1-eta)));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%% ELEMENT MATRICES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [LaplaceElementMatrix,MassElementMatrix] = ELEMENT_MATRICES(scaling)
% FIRST ORDER QUADRILATERAL FINITE ELEMENTS
aa=scaling/2; bb=scaling/2; % Element size scaling
k1=(aa^2+bb^2)/(aa*bb); k2=(aa^2-2*bb^2)/(aa*bb); k3=(bb^2-2*aa^2)/(aa*bb);
LaplaceElementMatrix = [k1/3 k2/6 -k1/6 k3/6 ; k2/6 k1/3 k3/6 -k1/6; ...
-k1/6  k3/6  k1/3  k2/6; k3/6 -k1/6  k2/6  k1/3];
MassElementMatrix = aa*bb*[4/9 2/9 1/9 2/9 ; 2/9 4/9 2/9 1/9 ; ...
1/9 2/9 4/9 2/9; 2/9 1/9 2/9 4/9];
end

Funding

This work was supported in by Villum Fonden through the NATEC (NAnophotonics for TErabit Communications) Centre (grant no. 8692) and by the Danish National Research Foundation through NanoPhoton Center for Nanophotonics (grant no. DNRF147).

Disclosures

The authors declare that there are no conflicts of interest related to this article.

References