Merge branch 'sans-master-diffGratTopOpt' into 'master'

Topology optimization on diffraction gratings

See merge request !8

DiffractionGratingsTopOpt/README.md

+# Directory in git repository : ```src_blaze/PhD/Maxwell_solvers/meshGrating_topopt_conical```
+
+File : ```README.md```\
+Author : ©Simon Ans\
+Last update : 05/03/2024\
+To get a pdf file from this ```.md``` file, run in a shell terminal :
+
+> ```pandoc README.md --pdf-engine=xelatex -o README.pdf```
+
+## About this directory
+
+---
+
+This directory gathers all the code files necessary to compute the mesh-based topology optimization of a periodic grating. Each element of the Finite Element triangulation can have a different permittivity in order to get the best reflection efficiency possible on a particular order for a given wavelength range.
+
+**One should read the paper "Topology optimization of blazed gratings under conical incidence" available on <u>link</u> before using this script.**
+
+## Quick Start
+
+---
+
+#### **Software Requirements**
+
+Make sure you have got all the following softwares and that they are updated :
+
+* ```gmsh``` 4.11.1 and more (**The location is given in the header of ```config_topopt_conical_data.py```**),
+
+* ```getdp``` 3.5.0 and more (**The location is given in the header of ```config_topopt_conical_data.py```**),
+
+* ```python``` 3.10 and more, with the libraries ```numpy```, ```scipy```, ```matplotlib.pyplot```, ```joblib```, ```os```, ```multiprocessing```, ```time```, ```meshio``` and most importantly ```nlopt``` (https://github.com/stevengj/nlopt). The code has been written using a machine with ```Miniconda3``` installed, hence the ```python``` command is used instead of ```python3```.
+
+#### **Lunch the optimization**
+
+While being in this directory, run in the shell terminal
+
+> ```$ python topopt_conical_main.py```
+
+The two main files of interest are ```rho_opt.pos``` and ```rho_opt.txt```. They define the optimal density distribution in two ways. The first one can be directly displayed with ```Gmsh``` and the second one is automatically displayed into the PDF file ```rho_opt.pdf```.
+
+Once the plots and the results are analysed and/or moved, run
+
+> ```$ sh topopt_conical_clean.sh```
+
+## More detailed description of the directory
+
+---
+
+#### General organization
+
+* ```config_topopt_conical_data.py``` : a Python file that gathers all the data for the simulation,
+
+* ```topopt_conical_main.py``` : the main Python file that handles all the resolution and plots the results,
+
+* ```resolution_functions.py``` : the Python file with all the functions that enable the resolution of the Finite Element problem and the plot of the results,
+
+* ```optimization_functions.py``` : the Python file with the definition of the target and the computation of the Jacobian. There is also the test function for the adjoint method,
+
+* ```tabs_material.py``` : a function that enables the main code to handle any material with its $\varepsilon_r$. See the website https://refractiveindex.info to add the materials,
+
+* ```topopt_conical.geo``` : the geometrical description of the pattern and the definition of the mesh using ```Gmsh```,
+
+* ```topopt_conical.pro``` : the computation file, that solves a given problem on the geometry defined by ```topopt_conical.geo``` with the Finite Element Method using GetDP,
+
+* ```topopt_conical_clean.sh``` : a Shell script that cleans the directory after using.
+  
+#### **Parameters and results obtained**
+
+The output is the local optimal density map ```rho_opt``` (depending on the initial density), saved in the files ```rho_opt.txt``` and ```rho_opt.pos```.
+
+The main parameters are as follows - they are all defined in the ```config``` Python file :
+
+| **Parameter**             | **Name in the scripts**    |  **Default value**   | **Description**|
+|---------------------------|----------------------------|----------------------|----------------|
+| $d$                       | ```period```               | 3300nm               | Period of the pattern |
+| $h_{\text{PML}^+}$        | ```h_PML_sup```            | 500nm                | Superior PML height |
+| $h_{\mathcal{S}^+}$       | ```h_super```              | 1500nm               | Superstrate height |
+| $h_{\mathcal{S}^-}$       | ```h_subs```               | 200nm                | Substrate height |
+| $h_{\text{PML}^-}$        | ```h_PML_inf```            | 200nm                | Inferior PML height |
+| $h_d$                     | ```h_design```             | 650nm                | Height of the design region |
+| $\Lambda_t$               | ```target_wavelength```    | [700nm]              | Target wavelengths for the optimization (Python list) |
+| $A_e$                     | ```amplitude```            | 1                    | Amplitude of the incident electric field |
+| $\theta_i$                | ```theta_i```              | 5°                   | Incidence angle on the (Oxy) plane |
+| $\varphi_i$               | ```phi_i```                | -66°                 | Incidence angle on the (Oyz) plane |
+| $\psi_i$                  | ```psi_i```                | 90°                  | Minimum polarization angle |
+| $N$                       | ```nb_orders```            | 1                    | Number of diffraction orders computed in one direction (total is $[-N,N]$) |
+|                           | ```wanted_order```         | -1                   | Order on which the optimization is proceeded |
+|                           | ```material_superstrate``` | Air                  | Material for the superstrate |
+|                           | ```material_min```         | Air                  | Material with minimal permittivity |
+|                           | ```material_max```         | Fused silica         | Material with maximal permittivity  |
+|                           | ```material_substrate```   | Silver (widerange)   | Material for the substrate |
+|                           | ```num_threads_multilambda```| ```len(lambda_target)```                   | Number of threads to run the multi-wavelength optimization in parallel |
+|                           | ```num_threads_resolution```| 1                   | Number of threads to run the Finite Element Method in parallel |
+|                           | ```tol```                  | 1e-5                 | Tolerance of oscillation of the target function |
+|                           | ```maxiter```              | 100                  | Maximum number of optimization iterations in a single binarization iteration |
+|                           | ```filter_radius```        | 200nm                | Radius of diffusion for the connectedness |
+
+#### Customization
+
+##### **A flexible script**
+
+There are other parameters in the ```config``` file, such as the initial configuration given in the ```analytic_init_density``` function (/!\ every length must be multiplied by ```nm```) and the refinement of the mesh with ```paramaille```, which is typically $\min(\Lambda_t)/h_{\max}$ ($h$ being the size of a mesh element). This size of the mesh can be defined differently in each region.
+
+There are also several flags in order to decide if one wants binarization, connectedness, a metal, etc.
+
+##### **Some other features**
+
+During the optimization, the directory ```res_conical``` is created and filled with a lot of data. The most interesting are the efficiency files ```efficiency_*.txt``` on the targeted diffraction order. The current efficiency is always updated but there are also the ```efficiency_r_*_list.txt``` files that keep the evolutions of each diffraction efficiency in reflection.
+
+Furthermore the total electric field is saved and updated. It can be directly displayed using ```Gmsh```. The adjoint fields are also available, although they are commented by default since they are not necessary for the optimization. The total absorption is also computed.
+
+Finally, at each optimization step, an image of the current density distribution is saved with the name ```myimage_*.png```. It enables to directly make a film of the optimization process.
+
+##### **Check the adjoint method**
+
+In order to be sure that the adjoint method provides the same derivatives as the finite difference method (with the step ```h_0``` in the ```config``` file), there is the function ```optimization_functions.test_Target_n```. One can run some tests. In order to launch a run quickly, it is better to take a small period (max 100nm) on the order 0. However the pattern can be random.
+
+The lines to launch this test is in the main function, and one can launch it by changing ```Flag_test``` to 1 in the ```config``` file. The relative difference between both Jacobians should be around $10^{-2}$ for a simple case.
\ No newline at end of file

DiffractionGratingsTopOpt/config_topopt_conical_data.py

+import numpy as np
+import matplotlib.pyplot as plt
+import meshio as mio
+import random
+
+'''
+© Simon Ans
+File : config_topopt_conical_main.py
+
+Gathers all the parameters of the simulation.
+        ____________________________________________________
+       | /!\ The reference unity for distances is nanometer |
+       ------------------------------------------------------
+'''
+
+## Paths to the softwares and documents -----------------------------------------------------------------------------------------------
+
+res_folder_name = "res_conical/"
+gmsh_path = ""
+GetDP_path = "~/programs/getdp/build/"
+mesh_file = "design_region.msh"
+
+## Global constants -------------------------------------------------------------------------------------------------------------------
+
+nm           = 1000      # resizing factor for GetDP
+ep0          = 8.854187817e-3 * nm
+mu0          = 400 * np.pi * nm
+cel          = 1 / (np.sqrt(ep0 * mu0))
+deg2rad      = np.pi/180
+
+## Parameters of the simulation -------------------------------------------------------------------------------------------------------
+
+## Dimensions ##
+period    = 3300
+h_PML_sup = 700          # superior PML thickness
+h_super   = 1000         # superstrate thickness
+h_subs    = 200          # substrate thickness
+h_PML_inf = 200          # inferior PML thickness
+
+h_design = 650           # thickness of the design space
+
+## Incident field ##
+target_wavelength = [700]
+#target_wavelength = [400,425,450,475,500,525,550,575,600,650,700,750,800,850,900,950,1000,1050,1100,1150,1200,1300,1400,1500]
+nb_wavelengths    = len(target_wavelength)
+amplitude         = 1
+# Incidence angles (in degrees)
+theta_i           = 5
+phi_i             = -66
+psi_i             = 90
+# phi = -90 and psi = 90  : 2D case with E// polarization
+#               psi = 180 : 2D case with H// polarization
+
+nb_orders    = 1
+wanted_order = -1                       # (integer) diffraction order on which the optimization is carried out
+target_order = nb_orders + wanted_order # (>= 0 integer) shifting of the wanted order for numerical purpose
+
+## Materials ##
+material_superstrate = 'air'        #'silicon_dioxide'  #'niobium_pentoxide'  #'PMMA'  #'silver_IEF' #silver_widerange  #'silver_Johnson_Christy'
+                                    #'ITO_Refractive'  #'BK7_schott' #'ITO_GT'  #'gold_ief_V2'  #'titanium'  #'titanium_dioxyde'
+                                    #'fused_silica'  #'silicon_nitride' #'silicon_internet2'  #'silver_palik'  #'copper_palik'
+                                    #'alu_palik'  #'gold_johnson' #'gold_olmon_evap'  #'chromium_BBM'  #'solar_etr'  #'solar_global_tilt'
+                                    #'solar_direct_circumsolar' #'silicon_salzberg'  #'silicon_odp0'  #'silicon'
+                                    #'silicon_webber'  #'silicon_internet'
+material_min    = 'air'           # minimal epsilon
+material_max    = 'fused_silica'  # maximal epsilon
+                                  # /!\ depending on the type of material, change the Flag_design_metal below !
+material_substrate = 'silver_widerange'
+
+## Numerical parameters ##
+paramaille  = 15         # minimum number of mesh elements per wavelength
+# Refinement coefficients on each region of the geometry (fraction of paramaille)
+paramaille_pmlSup      = 1
+paramaille_superstrate = 1
+paramaille_design      = 1.5
+paramaille_substrate   = 2
+paramaille_pmlInf      = 1
+
+## Parallel computation
+num_threads_multilambda = nb_wavelengths
+num_threads_resolution  = 1
+# /!\ num_threads_multilambda * num_threads_resolution must be < number of available threads on the machine
+
+## Initial density field ##
+def analytic_init_density(x,y):
+    return -(x - period*nm/2)/(period*nm)
+    # return random.random()
+
+## Optimization parameter ##
+# Boundaries of the densities for the simulation
+low_bnd = 0.
+up_bnd  = 1.
+
+# Stopping criteria
+tol = 1e-5
+maxiter = 100
+
+# Filters parameter
+filter_radius = 200   # diffusion
+nu = 1/2              # binarization
+nb_binarizations = 7
+
+# (For the finite differences only) Step of the finite differences
+h_0 = 1e-2
+
+## Choices for the simulation ##
+# Flags
+Flag_diffusion    = 1      # Impose the connectedness (1) or not (0)
+Flag_binarize     = 1      # Binarize the densities (1) or not (0)
+Flag_design_metal = 0      # Is the nanostructure made of metal (1) or not (0)
+Flag_pml_inf      = 0      # Do we need an inferior PML (1) or not (0)?
+Flag_test         = 0      # Launch an optimization (0) or a test on the Jacobian (1)
+Flag_o2_geom      = 0      # Type of finite elements (0 : tetrahedron ; 1 : curved)
+Flag_o2_interp    = 1      # Type of interpolation method (0 : P^2 ; 1 : more)
+
+## Identification of the surfaces of integration --------------------------------------------------------------------------------------
+#  (normally one does not have to change these values)
+PML_SUP     = 2000
+SUPERSTRATE = 2100
+DESIGN      = 2200
+SUBSTRATE   = 2300
+PML_INF     = 2400
+
+SURF_LEFT      = 100
+SURF_RIGHT     = 101
+SURF_UP        = 102
+SURF_DOWN      = 103
+SURF_INTEG_SUP = 120
+SURF_INTEG_SUB = 121
+SURF_PLOT      = 140
+PRINT_POINT    = 150
+
+## Optimization feedbacks -------------------------------------------------------------------------------------------------------------
+
+# Tool for the feedback
+optim_iter = 0
+np.savetxt("optim_iter.txt", np.array([optim_iter]))
+
+def plot_msh_density(current_msh_density, fig_msh_density_filename):
+    ''' Plot and save the density field written in a Python simple array '''
+
+    mesh = mio.read("design_region.msh")
+    fig  = plt.figure(figsize=(6,6))
+    im   = plt.tripcolor(mesh.points[:, 0],mesh.points[:, 1],mesh.cells_dict['triangle'], facecolors=current_msh_density, cmap='Blues', vmin=0, vmax=1)
+    plt.axis('scaled')
+    plt.axis('off')
+    plt.colorbar(im, fraction=0.046, pad=0.04, shrink=0.82)
+    plt.savefig(fig_msh_density_filename, bbox_inches='tight')
+    plt.close(fig)
+
+    return(0)
+
+def display_density(current_rho):
+    ''' Update the iteration of the optimization process and store the current density '''
+
+    global optim_iter
+    optim_iter += 1
+    print("###   Solving problem (direct and inverse) - iteration %03d...   ###"%optim_iter)
+    plot_msh_density(current_rho, "myimage_%03d.png"%optim_iter)
+    np.savetxt("optim_iter.txt", np.array([optim_iter]))
+
+    return(0)
\ No newline at end of file

DiffractionGratingsTopOpt/optimization_functions.py

+import numpy as np
+import matplotlib.pyplot as plt
+import meshio as mio
+import os
+import re
+import time
+from joblib import Parallel, delayed
+from sklearn.metrics.pairwise import euclidean_distances
+from config_topopt_conical_data import *
+import resolution_functions
+import tabs_material
+
+'''
+© Simon Ans
+File : optimization_functions.py
+'''
+
+##########################################################
+##                   DENSITY FUNCTIONS                  ##
+##   Settle and change the density field on the mesh    ##
+##########################################################
+
+def init_density(templatefile, outfile, initfunc):
+    ''' Initialize a density field on the mesh with an analytic function '''
+    
+    print("###### Initializing template file... ######")
+
+    # Solve a P^0 problem in order to spot the triangles of the mesh
+    os.system(GetDP_path + 'getdp topopt_conical.pro -pre res_init_densities -msh topopt_conical.msh -cal -v 1')
+
+    # Prepare the files that read and store the positions of the triangles
+    f_in  = open(templatefile, 'r')
+    f_out = open(outfile, 'w')
+    lines_in = f_in.readlines()
+    f_out.write(lines_in[0])
+    vals = []
+
+    # Store the value of the density inside each triangle of the mesh (detected through its barycentre)
+    for line in lines_in[1:-1]:
+        bary = np.mean(np.array(re.search(r'\((.*?)\)',line).group(1).split(','),float).reshape(3,3),axis=0)
+        val = initfunc(bary[0],bary[1])
+        f_out.write(re.sub('{.*?}','{%.10f,%.10f,%.10f}'%(val,val,val),line))
+        vals.append(val)
+
+    f_out.write(lines_in[-1])
+    f_in.close()
+    f_out.close()
+
+    print('######  Template file initialized.   ######\n')
+
+    return(np.array(vals, dtype=float))
+
+## -------------------------------------------------------------------------------------------------------------
+
+def get_barycenters():
+    ''' Get the barycenter matrix of the triangles of the mesh topopt_conical.msh '''
+
+    msh = mio.read(mesh_file)
+    res=[]
+    for k in range(len(msh.cells_dict['triangle'])):
+        t = msh.cells_dict['triangle'][k]
+        a,b,c = msh.points[t[0]],msh.points[t[1]],msh.points[t[2]]
+        res.append(np.mean(np.array([a,b,c]),axis=0).tolist())
+    
+    return np.array(res)
+
+## -------------------------------------------------------------------------------------------------------------
+
+def update_density(templatefile, outfile, tab):
+    ''' Update the density field with the current density in the optimization process '''
+
+    # Same as init_density
+    f_in  = open(templatefile, 'r')
+    f_out = open(outfile, 'w')
+    lines_in = f_in.readlines()
+    f_out.write(lines_in[0])
+        
+    # Update (read the table of the template mesh and put the values of the current densities)
+    counter = 0
+    for line in lines_in[1:-1]:
+        template_density = tab[counter]
+        f_out.write(re.sub('{.*?}','{%.10f,%.10f,%.10f}'%(template_density,template_density,template_density),line))
+        counter += 1
+    f_out.write(lines_in[-1])
+    f_in.close()
+    f_out.close()
+
+    return(0)
+
+## -------------------------------------------------------------------------------------------------------------
+
+def read_real_posfile(posfile):
+    ''' Read a complex .pos file to make an array '''
+
+    # Same as init_density
+    f_in  = open(posfile, 'r')
+    lines_in = f_in.readlines()
+    real_table = []
+
+    # Read the accurate inputs of the .pos file
+    for line in lines_in[1:-1]:
+        current_line = np.array(re.search(r'\{(.*?)\}',line).group(1).split(','),float)
+        real_table.append(current_line[0])
+    f_in.close()
+
+    real_table = np.array(np.real(real_table), dtype=float)
+    
+    return(real_table)
+
+## -------------------------------------------------------------------------------------------------------------
+
+def read_complex_posfile(posfile):
+    ''' Read a complex .pos file to make an array '''
+
+    # Same as init_density
+    f_in  = open(posfile, 'r')
+    lines_in = f_in.readlines()
+    complex_table = []
+
+    # Read the accurate inputs of the .pos file
+    for line in lines_in[1:-2]:
+        current_line = np.array(re.search(r'\{(.*?)\}',line).group(1).split(','),float)
+        complex_table.append(current_line[0]+1j*current_line[3])
+    f_in.close()
+
+    complex_table = np.array(np.real(complex_table), dtype=float) + 1j*np.array(np.imag(complex_table), dtype=float)
+    
+    return(complex_table)
+
+##########################################################
+##                CONSTRAINTS FUNCTIONS                 ##
+##   Apply the diffusion and/or binarization filters    ##
+##########################################################
+
+def diffusion_filter(density):
+    ''' Filter of the mesh to avoid the isolated points of density '''
+
+    midpoint = get_barycenters()
+    distance_temp = filter_radius*nm - euclidean_distances(midpoint, midpoint)
+    distance_temp[distance_temp < 0] = 0
+    distance_mat = distance_temp.copy()
+    distance_sum = distance_mat.sum(1)
+    rho_diffusion_Filtered = np.divide(distance_mat @ (density.flatten()), distance_sum)
+    return rho_diffusion_Filtered
+
+## -------------------------------------------------------------------------------------------------------------
+
+def diffusion_filter_jacobian():
+    ''' Jacobian of the diffusion filter '''
+
+    midpoint = get_barycenters()
+    distance_temp = filter_radius*nm - euclidean_distances(midpoint, midpoint)
+    distance_temp[distance_temp < 0] = 0
+    distance_mat = distance_temp.copy()
+    distance_sum = distance_mat.sum(1)
+    diffusion_jacobian = np.divide(distance_mat, distance_sum)
+    return diffusion_jacobian
+
+## -------------------------------------------------------------------------------------------------------------
+
+def binarize_filter(density, iter_bin):
+    ''' Binarize the mesh with a degree given by iter_bin '''
+
+    beta_f = 2**iter_bin
+    return ( (np.tanh(beta_f*nu) + np.tanh(beta_f*(density - nu))) / (np.tanh(beta_f*nu) + np.tanh(beta_f*(1 - nu))) )
+
+## -------------------------------------------------------------------------------------------------------------
+
+def binarize_filter_jacobian(density, iter_bin):
+    ''' Jacobian of the binarization filter '''
+
+    beta_f = 2**iter_bin
+    return np.diag( (beta_f * (1 - np.tanh(beta_f*(density-nu))**2)) / (np.tanh(beta_f*nu) + np.tanh(beta_f*(1 - nu))) )
+
+## -------------------------------------------------------------------------------------------------------------
+
+def composition_filters(density, iter_bin):
+    ''' Compose the filters previously defined depending on the flags '''
+
+    raw_density = density.copy()
+    if ((Flag_diffusion == 0) and (Flag_binarize == 0)): return raw_density
+    if ((Flag_diffusion == 1) and (Flag_binarize == 0)): return diffusion_filter(raw_density)
+    if ((Flag_diffusion == 0) and (Flag_binarize == 1)): return binarize_filter(raw_density, iter_bin)
+    if ((Flag_diffusion == 1) and (Flag_binarize == 1)): return binarize_filter(diffusion_filter(raw_density), iter_bin)
+
+## -------------------------------------------------------------------------------------------------------------
+
+def composition_filters_jacobian(density, iter_bin):
+    ''' Compute the Jacobian matrix of the composition of the filters depending on the flags '''
+
+    raw_density = density.copy()
+    problem_size = np.size(raw_density)
+    if ((Flag_diffusion == 0) and (Flag_binarize == 0)): return np.identity(problem_size)
+    if ((Flag_diffusion == 1) and (Flag_binarize == 0)): return diffusion_filter_jacobian()
+    if ((Flag_diffusion == 0) and (Flag_binarize == 1)): return binarize_filter_jacobian(raw_density, iter_bin)
+    if ((Flag_diffusion == 1) and (Flag_binarize == 1)): return np.matmul( diffusion_filter_jacobian(), binarize_filter_jacobian(diffusion_filter(raw_density), iter_bin) )
+
+## -------------------------------------------------------------------------------------------------------------
+
+def final_binarization(templatefile, nonbin_posfile):
+    ''' After the optimization, a final binarization to counter the last diffusion filter applied '''
+
+    rho_opt_nonbin = read_real_posfile(nonbin_posfile)
+    rho_opt_bin = binarize_filter(rho_opt_nonbin, nb_binarizations)
+    np.savetxt("rho_opt.txt", rho_opt_bin)
+    update_density(templatefile, "rho_opt.pos", rho_opt_bin)
+    last_iter = (nb_binarizations+1)*maxiter+1
+    plot_msh_density(rho_opt_bin, "myimage_%03d.png"%last_iter)
+
+    return(0)
+
+##########################################################
+##                OPTIMIZATION FUNCTIONS                ##
+##     Main functions for the topology optimization     ##
+##########################################################
+
+def compute_target_and_jacobian(i_wavelength):
+    ''' Calculate the Jacobian of the target using the complex amplitudes and their derivatives '''
+
+    wavelength = target_wavelength[i_wavelength]
+
+    # Components of the wave vectors (same reasoning as in the .pro file)
+    lambda0 = wavelength*nm
+    k0      = 2*np.pi / lambda0
+    n_super, k_super, epsilon_r_super = np.real(tabs_material.epsilon(material_superstrate, wavelength*1e-9))
+    deg2rad = np.pi/180
+
+    alpha      = -k0*n_super*np.sin(theta_i*deg2rad)*np.sin(phi_i*deg2rad)
+    beta_super = -k0*n_super*np.cos(theta_i*deg2rad)
+    gamma      = -k0*n_super*np.sin(theta_i*deg2rad)*np.cos(phi_i*deg2rad)
+
+    alpha_n      = alpha + 2*np.pi/(period*nm) * wanted_order
+    beta_n_super = -np.sqrt(epsilon_r_super*k0**2 - alpha_n**2 - gamma**2)
+
+    # Read the result files in order to get the complex amplitudes and their derivatives
+    comp_amp_x_tab = np.loadtxt(res_folder_name + "comp_amp_x_r_%d_%d.txt"%(wanted_order,i_wavelength))
+    comp_amp_z_tab = np.loadtxt(res_folder_name + "comp_amp_z_r_%d_%d.txt"%(wanted_order,i_wavelength))
+    comp_amp_x     = comp_amp_x_tab[1] + 1j*comp_amp_x_tab[2]
+    comp_amp_z     = comp_amp_z_tab[1] + 1j*comp_amp_z_tab[2]
+
+    comp_amp_x_derivative = read_complex_posfile(res_folder_name + "comp_amp_deriv_x_r_project_%d.pos"%i_wavelength)
+    comp_amp_z_derivative = read_complex_posfile(res_folder_name + "comp_amp_deriv_z_r_project_%d.pos"%i_wavelength)
+
+    target_jacobian = -2/(beta_n_super*beta_super*amplitude**2) * np.real((beta_n_super**2 + alpha_n**2) * comp_amp_x_derivative*np.conjugate(comp_amp_x)
+                                                                        + (beta_n_super**2 + gamma**2)   * comp_amp_z_derivative*np.conjugate(comp_amp_z)
+                                                                        + alpha_n*gamma * (comp_amp_x_derivative*np.conjugate(comp_amp_z) + comp_amp_x*np.conjugate(comp_amp_z_derivative)))
+
+    return(np.array(target_jacobian, dtype=float))
+
+## -------------------------------------------------------------------------------------------------------------
+
+def get_sum_targetJacobian(target_jacobian):
+    ''' Solve the direct and adjoint problems for each wavelength and add them together
+        in order to get the wanted efficiencies and their Jacobian '''
+
+    # Enable each wavelength to access the current density
+    for i_wavelength in range(nb_wavelengths):
+        os.system("cp -r rho_opt.pos topopt_conical_wavelength" + str(i_wavelength))
+
+    # Run the direct and adjoint problem to compute the target and its Jacobian 
+    temp_target_jacobian = np.zeros(np.size(target_jacobian))
+    temp_target   = 0
+    Parallel(n_jobs=num_threads_multilambda)(delayed(resolution_functions.solve)(i_wavelength) for i_wavelength in range(nb_wavelengths))
+    for i_wavelength in range(nb_wavelengths) :
+        temp_target_jacobian += compute_target_and_jacobian(i_wavelength)
+        temp_target   += np.loadtxt(res_folder_name + "efficiency_r_%d_%d.txt"%(wanted_order,i_wavelength))[1]
+
+    return temp_target, temp_target_jacobian
+
+## -------------------------------------------------------------------------------------------------------------
+
+def target(input_density, target_jacobian):
+    ''' Definition of the target function and its jacobian for the optimization problem
+        The fact that 'target_jacobian' is an input is a requirement of the NLopt class '''
+
+    iter_bin = int(np.loadtxt("iter_bin.txt"))
+
+    # Filter the densities depending on the filtering flags
+    updated_density  = composition_filters(input_density, iter_bin)
+    # Compute the Jacobian of this filter
+    filters_jacobian = composition_filters_jacobian(input_density, iter_bin)
+
+    # Store the filtered density file for the solver
+    update_density(res_folder_name + "template_density_table.pos", "rho_opt.pos", updated_density)
+
+    if (Flag_test == 0):
+        # Display the updated densities on the pattern
+        display_density(updated_density)
+    # Access the efficiencies and their Jacobian for each wavelength, and add them together
+    temp_target, temp_target_jacobian = get_sum_targetJacobian(target_jacobian)
+
+    # Update the Jacobian
+    if (np.size(target_jacobian) > 0):
+        target_jacobian[:] = 1/nb_wavelengths * np.matmul(filters_jacobian, temp_target_jacobian)
+        np.savetxt("jacobian.txt", target_jacobian)
+
+    # Update the target
+    target = 1 - temp_target/nb_wavelengths
+
+    return(target)
+
+## -------------------------------------------------------------------------------------------------------------
+
+def feedback_result(min_value, successful_termination):
+    ''' Displays the termination of the optimization with sentences '''
+
+    if (successful_termination > 0):
+        print("######              Optimization ended successfully.          ######\n")
+        if (successful_termination == 3):
+            print("The tolerance on the target function is reached.")
+        elif (successful_termination == 5):
+            print("The maximal number of iterations is reached.")
+        print("Reflection efficiency : %.6f"%(1-min_value))
+    else :
+        print("######          Issue encountered during optimization.        ######\n")
+        print(("Value of the termination : %d (see the nlopt reference https://nlopt.readthedocs.io/en/latest/ for details).\n"%successful_termination))
+    
+    return(0)
+
+##########################################################
+##                   TEST FUNCTIONS                     ##
+##    Comparison adjoint / finite element derivatives   ##
+##########################################################
+
+def test_jacobian(rho_init, iter_bin):
+    ''' Compare the derivatives computed with the adjoint problem and with the finite differences for the target function
+        * The first one is given by the Jacobian computed by target(rho_init, jacobian)
+        * The second one is the result of (target(rho_init,...) - target(rho_init + h_0,...)) / h_0 on each triangle of the mesh '''
+
+    print("######  Test of the adjoint method   ######\n")
+    plot_msh_density(rho_init, "rho_init.pdf")
+    np.savetxt("iter_bin.txt", np.array([iter_bin]))
+    problem_size = np.size(rho_init)
+
+    jacobian = np.zeros(problem_size)
+    target_rho_plus_h0 = np.zeros_like(jacobian)
+
+    # Estimation of the test duration ---------------------------------------------
+    t1 = time.time()
+
+    # Solve the direct and adjoint problems in order to access the target with the initial densities
+    # and the Jacobian computed by the adjoint method
+    target_rho = target(rho_init, jacobian)
+    numerical_jacobian_adjoint = np.loadtxt("jacobian.txt")
+
+    t2 = time.time()
+    print("## Finite differences computation : beginning the loop on %d triangles."%problem_size)
+    print("## Estimated remaining time : %dmin%02ds"%((t2-t1)*problem_size//60, (t2-t1)*problem_size%60))
+    # -----------------------------------------------------------------------------
+    
+    # Loop on the triangles to compute the finite differences for each one    
+    for i_triangle in range(problem_size):
+        rho_plus_h0 = rho_init.copy()
+        rho_plus_h0[i_triangle] += h_0
+
+        # Compute the target with the slightly different densities
+        target_rho_plus_h0[i_triangle] = target(rho_plus_h0, jacobian)
+
+    # Final computation of the finite differences
+    numerical_jacobian_finite_diff = (target_rho_plus_h0 - target_rho) / h_0
+    # Relative difference between both methods
+    diff_rel = (numerical_jacobian_finite_diff - numerical_jacobian_adjoint) / np.mean(np.abs(numerical_jacobian_finite_diff))
+
+    # Save the results in text files
+    np.savetxt("adjoint.txt", numerical_jacobian_adjoint)
+    np.savetxt("finite_diff.txt", numerical_jacobian_finite_diff)
+    np.savetxt("diff_rel.txt", diff_rel)
+
+    # Plot the results on the pattern
+    plot_msh_density(numerical_jacobian_adjoint, "adjoint_jacobian_Target_n.pdf")
+    plot_msh_density(numerical_jacobian_finite_diff, "finite_diff_jacobian_Target_n.pdf")
+    plot_msh_density(np.abs(diff_rel), "diff_jacobian_Target_n.pdf")
+
+    print("######         End of the test       ######\n")
+
+    return(0)
+
+## -------------------------------------------------------------------------------------------------------------
+
+def test_jacobian_soloPoint(rho_init, iter_bin, x_point, y_point):
+    ''' Compare the derivatives computed with the adjoint problem and with the finite differences for the target function,
+        but only for one specific point '''
+
+    print("######  Test of the adjoint method   ######\n")
+    plot_msh_density(rho_init, "rho_init.pdf")
+    np.savetxt("iter_bin.txt", np.array([iter_bin]))
+    problem_size = np.size(rho_init)
+    jacobian = np.zeros(problem_size)
+
+    # Spot the closest triangle to the point (x_point, y_point)
+    midpoint = get_barycenters()/nm
+    closest_point = np.argmin( np.sqrt((x_point-midpoint[:,0])**2 + (y_point-midpoint[:,1])**2) )
+
+    # Solve the direct and adjoint problems in order to access the target with the initial densities
+    # and the Jacobian computed by the adjoint method
+    target_rho = target(rho_init, jacobian)
+    numerical_jacobian_adjoint = np.loadtxt("jacobian.txt")[closest_point]
+
+    # Compute the target on the triangle closest_point
+    rho_plus_h0 = rho_init.copy()
+    rho_plus_h0[closest_point] += h_0
+    target_rho_plus_h0 = target(rho_plus_h0, jacobian)
+
+    # Final computation of the finite differences
+    numerical_jacobian_finite_diff = (target_rho_plus_h0 - target_rho) / h_0
+    # Relative difference between both methods
+    diff_rel = (numerical_jacobian_finite_diff - numerical_jacobian_adjoint) / np.abs(numerical_jacobian_finite_diff)
+    print("######         End of the test       ######\n")
+    print("Jacobian computed : * with the adjoint method     on the point (%.0f,%.0f) : %.5f"%(x_point,y_point,numerical_jacobian_adjoint))
+    print("                    * with the finite differences on the point (%.0f,%.0f) : %.5f"%(x_point,y_point,numerical_jacobian_finite_diff))
+    print("Relative differences between both Jacobian methods : %.5f"%diff_rel)
+    
+    return(0)
\ No newline at end of file

DiffractionGratingsTopOpt/resolution_functions.py

+import numpy as np
+import matplotlib.pyplot as plt
+import os
+import multiprocessing
+import tabs_material
+from config_topopt_conical_data import *
+
+##########################################################
+##             FUNCTIONS USING THE TERMINAL             ##
+##                  Use Gmsh and GetDP                  ##
+##########################################################
+
+def generate_mesh():
+    ''' Once the file param.dat is created, run this line to generate the corresponding mesh '''
+    print("######            Meshing...         ######")
+    os.system(gmsh_path + "gmsh topopt_conical.geo -v 0 -")
+    print("######       Meshing completed.      ######\n")
+
+    return(0)
+
+## -------------------------------------------------------------------------------------------------------------
+
+def make_computation_directories():
+    ''' Prepare the directories to run the resolution for every wavelength
+        Adapt the param.dat file'''
+
+    # Initial wavelength
+    initial_wavelength = target_wavelength[0]
+    # Initial permittivities and permeabilities
+    initial_n_super, initial_k_super, initial_epsilon_r_super = tabs_material.epsilon(material_superstrate, target_wavelength[0]*1e-9)
+    initial_n_min  , initial_k_min  , initial_epsilon_r_min   = tabs_material.epsilon(material_min, target_wavelength[0]*1e-9)
+    initial_n_max  , initial_k_max  , initial_epsilon_r_max   = tabs_material.epsilon(material_max, target_wavelength[0]*1e-9)
+    initial_n_subs , initial_k_subs , initial_epsilon_r_subs  = tabs_material.epsilon(material_substrate, target_wavelength[0]*1e-9)
+
+    for i_wavelength in range(nb_wavelengths):
+        ## Make directory and copy the useful files
+        os.system("mkdir topopt_conical_wavelength" + str(i_wavelength))
+        os.system("cp -r param.dat topopt_conical.geo topopt_conical.msh topopt_conical.pro tabs_material.py topopt_conical_wavelength" + str(i_wavelength))
+    
+        ## Adapt the param.dat file
+        os.chdir("topopt_conical_wavelength" + str(i_wavelength))
+
+        # Wavelength and iteration
+        current_wavelength  = target_wavelength[i_wavelength]
+        os.system("sed -i 's/lambda_wavelength = %.15e;/lambda_wavelength = %.15e;/g' param.dat"%(initial_wavelength, current_wavelength))
+        os.system("sed -i 's/i_lambda = 0;/i_lambda = %d;/g' param.dat"%i_wavelength)
+
+        # Permittivities and permeabilities for every wavelength
+        current_n_super, current_k_super, current_epsilon_r_super = tabs_material.epsilon(material_superstrate, target_wavelength[i_wavelength]*1e-9)
+        current_n_min  , current_k_min  , current_epsilon_r_min   = tabs_material.epsilon(material_min, target_wavelength[i_wavelength]*1e-9)
+        current_n_max  , current_k_max  , current_epsilon_r_max   = tabs_material.epsilon(material_max, target_wavelength[i_wavelength]*1e-9)
+        current_n_subs , current_k_subs , current_epsilon_r_subs  = tabs_material.epsilon(material_substrate, target_wavelength[i_wavelength]*1e-9)
+
+        os.system("sed -i 's/epsilon_r_super_re = %.15e;/epsilon_r_super_re = %.15e;/g' param.dat"%(np.real(initial_epsilon_r_super), np.real(current_epsilon_r_super)))
+        os.system("sed -i 's/epsilon_r_super_im = %.15e;/epsilon_r_super_im = %.15e;/g' param.dat"%(np.imag(initial_epsilon_r_super), np.imag(current_epsilon_r_super)))
+        os.system("sed -i 's/n_super = %.15e;/n_super = %.15e;/g' param.dat"%(initial_n_super, current_n_super))
+        os.system("sed -i 's/k_super = %.15e;/k_super = %.15e;/g' param.dat"%(initial_k_super, current_k_super))
+        os.system("sed -i 's/epsilon_r_min_re = %.15e;/epsilon_r_min_re = %.15e;/g' param.dat"%(np.real(initial_epsilon_r_min), np.real(current_epsilon_r_min)))
+        os.system("sed -i 's/epsilon_r_min_im = %.15e;/epsilon_r_min_im = %.15e;/g' param.dat"%(np.imag(initial_epsilon_r_min), np.imag(current_epsilon_r_min)))
+        os.system("sed -i 's/n_min = %.15e;/n_min = %.15e;/g' param.dat"%(initial_n_min, current_n_min))
+        os.system("sed -i 's/k_min = %.15e;/k_min = %.15e;/g' param.dat"%(initial_k_min, current_k_min))
+        os.system("sed -i 's/epsilon_r_max_re = %.15e;/epsilon_r_max_re = %.15e;/g' param.dat"%(np.real(initial_epsilon_r_max), np.real(current_epsilon_r_max)))
+        os.system("sed -i 's/epsilon_r_max_im = %.15e;/epsilon_r_max_im = %.15e;/g' param.dat"%(np.imag(initial_epsilon_r_max), np.imag(current_epsilon_r_max)))
+        os.system("sed -i 's/n_max = %.15e;/n_max = %.15e;/g' param.dat"%(initial_n_max, current_n_max))
+        os.system("sed -i 's/k_max = %.15e;/k_max = %.15e;/g' param.dat"%(initial_k_max, current_k_max))
+        os.system("sed -i 's/epsilon_r_subs_re = %.15e;/epsilon_r_subs_re = %.15e;/g' param.dat"%(np.real(initial_epsilon_r_subs), np.real(current_epsilon_r_subs)))
+        os.system("sed -i 's/epsilon_r_subs_im = %.15e;/epsilon_r_subs_im = %.15e;/g' param.dat"%(np.imag(initial_epsilon_r_subs), np.imag(current_epsilon_r_subs)))
+        os.system("sed -i 's/n_subs = %.15e;/n_subs = %.15e;/g' param.dat"%(initial_n_subs, current_n_subs))
+        os.system("sed -i 's/k_subs = %.15e;/k_subs = %.15e;/g' param.dat"%(initial_k_subs, current_k_subs))
+        os.chdir("..")
+
+    return(0)
+
+## -------------------------------------------------------------------------------------------------------------
+
+def solve(i_wavelength):
+    ''' Launch GetDP to solve the problem on the mesh generated with the current data file '''
+
+    current_wavelength = target_wavelength[i_wavelength]
+    # Move into the corresponding directory
+    os.chdir("topopt_conical_wavelength" + str(i_wavelength))
+
+    os.system(GetDP_path + "getdp topopt_conical.pro -pre helmholtz_conical_direct_inverse -msh topopt_conical.msh -cal -petsc_prealloc 200 -setnumber lambda0 " + str(current_wavelength * nm) + " -v 1")
+
+    # Move again into the main directory
+    os.chdir("..")
+
+    return(0)
+
+##########################################################
+##                 TOOLS FOR THE DATA                   ##
+##            Functions to reduce set_data              ##
+##########################################################
+
+def check_quantites():
+    ''' In the set_data program, check if the quantities are coherent '''
+
+    if h_design > h_super:
+        raise ValueError("h_design > h_super, please increase h_super.")
+    if ((theta_i > 90) or (theta_i < 0)):
+        raise ValueError("The angle of incidence out of the plane is not in the correct interval (0° < theta_i < 90°)")
+    if ((phi_i > 180) or (phi_i < -180)):
+        raise ValueError("The angle of incidence in the plane is not in the correct interval (-180° < phi_i < 180°)")
+    if ((psi_i > 180) or (psi_i < -180)):
+        raise ValueError("The angle of polarization is not in the correct interval (-180° < psi_i < 180°)")
+    if multiprocessing.cpu_count() < num_threads_multilambda * num_threads_resolution :
+        raise ValueError("The number of CPUs of the current machine is too low regarding the number of threads allocated (for the multi-wavelength and the Finite Element resolutions).")
+    if paramaille < 10:
+        print("Warning : The number of mesh elements per wavelength is quite low (minimum should be 10). The result might not be precise enough.")
+        wait = input("Do you want to continue anyway ? (Press Enter if yes)")
+
+    return(0)
+
+## -------------------------------------------------------------------------------------------------------------
+
+def write_param(integer_parameters_name, integer_parameters_value, float_parameters_name, float_parameters_value):
+    ''' Write in a file called param.dat all the data contained in parameters '''
+
+    param = open("param.dat","w")
+
+    length_integers = len(integer_parameters_name)
+    length_float    = len(float_parameters_name)
+
+    for i_integers in range (length_integers):
+        param.write("%s = %d;\n"%(integer_parameters_name[i_integers], integer_parameters_value[i_integers]))
+    
+    for i_float in range (length_float):
+        param.write("%s = %.15e;\n"%(float_parameters_name[i_float], float_parameters_value[i_float]))
+
+    param.close()
+
+    return(0)
+
+##########################################################
+##                   MAIN FUNCTIONS                     ##
+##       Parameters and resolution of the problem       ##
+##########################################################
+
+def set_data():
+    ''' Input the data into a file to create the periodic 2D geometry with its physical parameters
+        It is read by the .geo and .pro files '''    
+
+    ## Check if the quantities are coherent  -------------------------------------------------------------------
+
+    check_quantites()
+
+    ## Permittivities ------------------------------------------------------------------------------------------
+    n_super, k_super, epsilon_r_super = tabs_material.epsilon(material_superstrate, target_wavelength[0]*1e-9)
+    n_min, k_min, epsilon_r_min = tabs_material.epsilon(material_min, target_wavelength[0]*1e-9)
+    n_max, k_max, epsilon_r_max = tabs_material.epsilon(material_max, target_wavelength[0]*1e-9)
+    n_subs, k_subs, epsilon_r_subs    = tabs_material.epsilon(material_substrate, target_wavelength[0]*1e-9)
+
+    epsilon_r_super_re = epsilon_r_super.real
+    epsilon_r_super_im = epsilon_r_super.imag
+    mu_r_super_re = 1
+    mu_r_super_im = 0
+    epsilon_r_min_re = epsilon_r_min.real
+    epsilon_r_min_im = epsilon_r_min.imag
+    epsilon_r_max_re = epsilon_r_max.real
+    epsilon_r_max_im = epsilon_r_max.imag
+    mu_r_voxel_re = 1
+    mu_r_voxel_im = 0
+    epsilon_r_subs_re = epsilon_r_subs.real
+    epsilon_r_subs_im = epsilon_r_subs.imag
+    mu_r_subs_re = 1
+    mu_r_subs_im = 0    
+
+    ## Resize the quantities -----------------------------------------------------------------------------------
+
+    Period  = nm * period
+    hpmlsup = nm * h_PML_sup
+    hsuper  = nm * h_super
+    hdesign = nm * h_design
+    hsubs   = nm * h_subs
+    hpmlinf = nm * h_PML_inf
+
+    ## Write in the file ---------------------------------------------------------------------------------------
+
+    # list of the names of the integer parameters to write the param.dat file
+    integer_parameters_name = ["i_lambda", "target_order", "nb_orders", "Flag_o2_geom", "Flag_o2_interp", "Flag_pml_inf", "Flag_design_metal", "PML_SUP", "SUPERSTRATE", "SUBSTRATE", "PML_INF", "SURF_LEFT", "SURF_RIGHT", "SURF_UP", "SURF_DOWN", "SURF_INTEG_SUP", "SURF_INTEG_SUB", "SURF_PLOT", "PRINT_POINT", "DESIGN"]
+
+    # list of the values of the integers parameters to write the param.dat file
+    integer_parameters_value = [0, target_order, nb_orders, Flag_o2_geom, Flag_o2_interp, Flag_pml_inf, Flag_design_metal, PML_SUP,  SUPERSTRATE, SUBSTRATE, PML_INF, SURF_LEFT, SURF_RIGHT, SURF_UP, SURF_DOWN, SURF_INTEG_SUP, SURF_INTEG_SUB, SURF_PLOT, PRINT_POINT, DESIGN]
+
+    # list of the names of the float parameters to write the param.dat file
+    float_parameters_name = ["lambda_wavelength", "nm", "ep0", "mu0", "cel", "deg2rad", "period", "h_PML_sup", "h_super", "h_design", "h_subs", "h_PML_inf", "paramaille", "paramaille_pmlSup", "paramaille_superstrate", "paramaille_design", "paramaille_substrate", "paramaille_pmlInf", "theta_incidence", "phi_incidence", "psi_incidence", "A"]
+    float_parameters_name += ["epsilon_r_super_re", "epsilon_r_super_im", "n_super", "k_super", "mu_r_super_re", "mu_r_super_im", "epsilon_r_min_re", "epsilon_r_min_im", "n_min", "k_min", "epsilon_r_max_re", "epsilon_r_max_im", "n_max", "k_max", "mu_r_voxel_re", "mu_r_voxel_im", "epsilon_r_subs_re", "epsilon_r_subs_im", "n_subs", "k_subs", "mu_r_subs_re", "mu_r_subs_im"]
+    
+    # list of the values of the float parameters to write the param.dat file
+    float_parameters_value = [target_wavelength[0], nm, ep0, mu0, cel, deg2rad, Period, hpmlsup, hsuper, hdesign, hsubs, hpmlinf, paramaille, paramaille_pmlSup, paramaille_superstrate, paramaille_design, paramaille_substrate, paramaille_pmlInf, theta_i, phi_i, psi_i, amplitude]
+    float_parameters_value += [epsilon_r_super_re, epsilon_r_super_im, n_super, k_super, mu_r_super_re, mu_r_super_im, epsilon_r_min_re, epsilon_r_min_im, n_min, k_min, epsilon_r_max_re, epsilon_r_max_im, n_max, k_max, mu_r_voxel_re, mu_r_voxel_im, epsilon_r_subs_re, epsilon_r_subs_im, n_subs, k_subs, mu_r_subs_re, mu_r_subs_im]
+
+    write_param(integer_parameters_name, integer_parameters_value, float_parameters_name, float_parameters_value)
+
+    return(0)
+
+## -------------------------------------------------------------------------------------------------------------
+
+def print_sum_energies():
+    ''' Check the sum of the efficiencies and losses (must be 1 is nb_orders is high) '''
+
+    sum_efficiency_r = np.zeros(nb_wavelengths)
+    sum_efficiency_t = np.zeros_like(sum_efficiency_r)
+    losses_tot       = np.zeros_like(sum_efficiency_r)
+    sum_energies     = np.zeros_like(sum_efficiency_r)
+    for i_wavelength in range (nb_wavelengths):
+        for i_order in range (2*nb_orders+1):
+            sum_efficiency_r[i_wavelength] += np.loadtxt(res_folder_name + "efficiency_r_%d_%d.txt"%(i_order-nb_orders,i_wavelength))[1]
+            sum_efficiency_t[i_wavelength] += np.loadtxt(res_folder_name + "efficiency_t_%d_%d.txt"%(i_order-nb_orders,i_wavelength))[1]
+        losses_tot[i_wavelength] = np.loadtxt(res_folder_name + "absorption-Q_tot_%d.txt"%i_wavelength)[1]
+        sum_energies[i_wavelength] = sum_efficiency_r[i_wavelength] + sum_efficiency_t[i_wavelength] + losses_tot[i_wavelength]
+
+    print("Sum of the efficiencies and losses for the wavelengths of target_wavelength :")
+    print(sum_energies)
+
+    return(0)
\ No newline at end of file

DiffractionGratingsTopOpt/topopt_conical.geo

+/*
+© Simon Ans
+File : topopt_conical.geo
+
+Construction of the geometry of the periodical grating made with sticks using
+the data in param.dat. The geometry is detailed in resolution_functions.set_data.
+*/
+
+Include "param.dat";
+SetFactory("OpenCASCADE");
+
+ref_paramaille = 400 * nm / paramaille;
+lc_pmlSup      = ref_paramaille / paramaille_pmlSup;
+lc_superstrate = ref_paramaille / paramaille_superstrate;
+lc_slide       = ref_paramaille / paramaille_design;
+lc_substrate   = ref_paramaille / paramaille_substrate;
+lc_pmlInf      = ref_paramaille / paramaille_pmlInf;
+
+geometry_sup = h_super + h_PML_sup;
+geometry_inf = h_subs  + h_PML_inf;
+
+// Box
+Rectangle(1) = {-period/2, h_super      , 0, period, h_PML_sup};
+Rectangle(2) = {-period/2, h_design     , 0, period, h_super - h_design};
+Rectangle(3) = {-period/2, 0            , 0, period, h_design};
+Rectangle(4) = {-period/2, -h_subs      , 0, period, h_subs};
+Rectangle(5) = {-period/2, -geometry_inf, 0, period, h_PML_inf};
+
+Coherence;
+
+Lines_Left()    = Line In BoundingBox{-period/2 - 1, -geometry_inf - 1, -1, -period/2 + 1,  geometry_sup + 1, 1};
+Lines_Right()   = Line In BoundingBox{ period/2 - 1, -geometry_inf - 1, -1,  period/2 + 1,  geometry_sup + 1, 1};
+Lines_Up()      = Line In BoundingBox{-period/2 - 1,  geometry_sup - 1, -1,  period/2 + 1,  geometry_sup + 1, 1};
+Lines_Down()    = Line In BoundingBox{-period/2 - 1, -geometry_inf - 1, -1,  period/2 + 1, -geometry_inf + 1, 1};
+Lines_Super()   = Line In BoundingBox{-period/2 - 1,  h_super      - 1, -1,  period/2 + 1,  h_super      + 1, 1};
+Lines_Sub()     = Line In BoundingBox{-period/2 - 1, -h_subs       - 1, -1,  period/2 + 1, -h_subs       + 1, 1};
+
+phys_plot_bnd()  = Lines_Super();
+phys_plot_bnd() += Line In BoundingBox{-period/2 - 1, h_design - 1, -1,  period/2 + 1, h_design + 1, 1};
+phys_plot_bnd() += Line In BoundingBox{-period/2 - 1,         - 1, -1,  period/2 + 1,         + 1, 1};
+phys_plot_bnd() += Lines_Sub() ;
+
+Physical Line(SURF_LEFT)  = Lines_Left();
+Physical Line(SURF_RIGHT) = Lines_Right();
+Physical Line(SURF_UP)    = Lines_Up();
+Physical Line(SURF_DOWN)  = Lines_Down();
+
+Physical Line(SURF_INTEG_SUP) = Lines_Super(); // superstrate/pml cut
+Physical Line(SURF_INTEG_SUB) = Lines_Sub();   // substrate/pml cut
+Physical Line(SURF_PLOT) = phys_plot_bnd[];    // final plot
+
+Physical Surface(PML_SUP)     = {1};
+Physical Surface(SUPERSTRATE) = {2};
+Physical Surface(DESIGN)      = {3};
+Physical Surface(SUBSTRATE)   = {4};
+Physical Surface(PML_INF)     = {5};
+
+Characteristic Length{PointsOf{Physical Surface{PML_SUP};}}     = lc_pmlSup;
+Characteristic Length{PointsOf{Physical Surface{PML_INF};}}     = lc_pmlInf;
+Characteristic Length{PointsOf{Physical Surface{SUPERSTRATE};}} = lc_superstrate;
+Characteristic Length{PointsOf{Physical Surface{SUBSTRATE};}}   = lc_substrate;
+Characteristic Length{PointsOf{Physical Surface{DESIGN};}}      = lc_slide;
+
+Periodic Line {Lines_Left()} = {Lines_Right()} Translate {period,0,0};
+
+Mesh 2;
+Save "topopt_conical.msh";
+
+Delete Physicals;
+Physical Surface(DESIGN) = {3};
+Save "design_region.msh"; 
\ No newline at end of file

DiffractionGratingsTopOpt/topopt_conical_clean.sh

+#!/bin/bash
+# File : topopt_conical_clean.sh
+# Clean the folder by deleting all the files created by the computation
+
+rm *.msh 2> /dev/null
+rm *.txt 2> /dev/null
+rm *.npz 2> /dev/null
+rm *.dat 2> /dev/null
+rm *.pre 2> /dev/null
+rm *.out 2> /dev/null
+rm *.pdf 2> /dev/null
+rm *.pos 2> /dev/null
+rm *.png 2> /dev/null
+rm *.jpg 2> /dev/null
+rm *.mp4 2> /dev/null
+rm -rf topopt_conical_wavelength* 2> /dev/null
+rm -rf __pycache__/ 2> /dev/null
+rm -rf res_conical/ 2> /dev/null

DiffractionGratingsTopOpt/topopt_conical_main.py

+import numpy as np
+import nlopt
+import resolution_functions
+import optimization_functions
+import matplotlib.pyplot as plt
+import os
+import meshio as mio
+from config_topopt_conical_data import *
+
+'''
+© Simon Ans
+File : topopt_conical_main.py
+
+Optimize the permittivities on the mesh to maximize the reflection efficiencies
+for different wavelengths on a given order.
+'''
+
+if __name__=="__main__":
+    
+    ## Prepare the computation ------------------------------------------------------------------------------------------------------------
+
+    # Parallelization of the Finite Element Method
+    os.system("export OPENBLAS_NUM_THREADS=%d"%num_threads_resolution)
+
+    # Create the paramater file and the mesh
+    resolution_functions.set_data()
+    resolution_functions.generate_mesh()
+    resolution_functions.make_computation_directories()
+    
+    ## Density field ----------------------------------------------------------------------------------------------------------------------
+
+    msh_density_init = optimization_functions.init_density(res_folder_name + "template_density_table.pos", "rho_opt.pos", analytic_init_density)
+    #msh_density_init = np.loadtxt("rho_opt.txt")
+
+    if (Flag_test == 1):
+        ## Test of the adjoint method -----------------------------------------------------------------------------------------------------
+        if (Flag_design_metal == 1):
+            optimization_functions.test_jacobian(msh_density_init, 7)
+            #optimization_functions.test_jacobian_soloPoint(msh_density_init, 7, 40, 300)
+        else :
+            optimization_functions.test_jacobian(msh_density_init, 3)
+            #optimization_functions.test_jacobian_soloPoint(msh_density_init, 3, 40, 350)
+        resolution_functions.print_sum_energies()
+
+    else :
+        ## Optimization process -----------------------------------------------------------------------------------------------------------
+        
+        problem_size = np.size(msh_density_init)
+        # Choose the MMA optimization (actually GCMMA from the 2002 paper of Svanberg)
+        opt = nlopt.opt(nlopt.LD_MMA, problem_size)
+
+        # Define the target
+        opt.set_min_objective(optimization_functions.target)
+
+        # Boundaries
+        lb_array = low_bnd*np.ones(problem_size)
+        ub_array = up_bnd *np.ones(problem_size)
+        opt.set_lower_bounds(lb_array)
+        opt.set_upper_bounds(ub_array)
+
+        # Stopping criteria
+        opt.set_ftol_abs(tol)    
+        opt.set_maxeval(maxiter)
+
+        if (Flag_binarize == 1):
+            # Optimization with binarization
+            for iter_bin in range(nb_binarizations+1):
+                np.savetxt("iter_bin.txt", np.array([iter_bin]))
+                print("######                Binarization iteration %d...            ######"%iter_bin)
+                print("----------------------------------------------------------------------------------------------\n")
+                rho_opt = opt.optimize(msh_density_init)
+                msh_density_init = rho_opt.copy()
+            optimization_functions.final_binarization(res_folder_name + "template_density_table.pos", "rho_opt.pos")
+        else :
+            np.savetxt("iter_bin.txt", np.array([0]))
+            rho_opt = opt.optimize(msh_density_init)
+
+        plot_msh_density(rho_opt, "rho_opt.pdf")
+        plot_msh_density(np.loadtxt("jacobian.txt"), "jac_opt.pdf")
+        
+        # Print the optimization outcome
+        sum_energies = resolution_functions.print_sum_energies()
+        min_value = opt.last_optimum_value()
+        successful_termination = opt.last_optimize_result()
+        optimization_functions.feedback_result(min_value, successful_termination)
\ No newline at end of file