## Materials and Methods
### Simulation  
The numerical models of the designed metasurface structures were developed and computed in the Comsol Multiphysics environment, based on the finite element method (FEM). The model geometry and description, as well as the computational mesh, are shown in Figure 2. The electromagnetic field excitation and detection were implemented through Port 1 and Port 2. The plane wavefront was generated by Port 1, and both ports received the scattered field (Port 1 for the \(s_{11}\) reflection coefficient and Port 2 for the \(s_{21}\) transmission coefficient determination). Because of the advances in electromagnetic field computing technology, in the case of the planned analyzes, there was no need to model the entire structure of the metasurface. In the proposed model, the Floquet-periodic boundary conditions were used in all of the external and parallel to \(z\)-direction boundary cases, to simulate an infinite 2D array of unit cell elements. This resulted in a significant reduction of computational mesh size, and thus a simplification of the numerical model. This assumption restricts the domain of the potential geometry deformations to the planar ones (in the \(x\)-\(y\) plane). Such a simplification of the model is acceptable, as it will be shown later in this paper, by comparison of the modeled and measured \(s_{11}\) reflection coefficient frequency responses. Perfectly matched layers (PMLs) are situated on the top and bottom of the computational domain in order to absorb the mode excited from the source port (Port 1), and any higher-order modes potentially generated by the periodic structure. The models were calculated in the frequency domain for 2-8 GHz frequency bands.  
\begin{table}
\begin{tabular}{c c c c} \hline \hline
**Parameter** & **Cross Unit Cell** & **Jerusalem-Cross Unit Cell** & **Split Ring Resonator Unit Cell** \\ \hline \(d\) [mm] & 6 & 6 & 6 \\ \(w\) [mm] & 0.9 & 0.9 & 0.9 \\ \(l,l_{1}\) [mm] & 19.08 & 8.75 & – \\ \(l_{2}\) [mm] & – & 6.87 & – \\ \(r\) [mm] & – & – & 3.53 \\ \hline \hline \end{tabular}
\end{table}
Table 1: The dimensions of designed unit cells.  
One of the analyzes carried out as part of the experiment concerned the impact of the surface deformation caused by the external mechanical forces. In this case, one of the side unit cell boundaries was fixed in _x,y,z_ space, while the opposite one was loaded with an external force, as it is presented in Figure 3. First, the mechanical problem was solved and the deformed geometry was obtained. Next, the high-frequency electromagnetic problem was solved, based on the previously achieved deformed geometry. The equations solved during numerical the analysis are as follows [26]:  
\[\
abla\times\mu_{\mathrm{r}}^{-1}(\
abla\times\mathbf{E})-k_{0}^{2}\bigg{(} \varepsilon_{\mathrm{r}}-\frac{j\sigma_{\mathrm{e}}}{\omega\varepsilon_{0}} \bigg{)}\mathbf{E}=0\]  
\[0=\
abla\cdot S+\mathbf{F}_{\mathrm{V}}\]  
where \(\mathbf{E}\) is the electric field, \(\mu_{\mathrm{r}}\) is the relative permeability, \(\varepsilon_{\mathrm{r}}\)-- is the relative permittivity, \(\sigma\) is the electric conductivity, \(k_{0}\) is the wave number, \(\omega\) is the angular frequency, \(S\) is the stress, and \(\mathbf{F}_{\mathrm{v}}\) is the force per unit volume. The exemplary mesh used during numerical calculations consisted of 11,250 domain elements, 2662 boundary elements, and 380 edge elements.  
Figure 3: External mechanical load application (in a \(y\)-direction).  
Figure 2: Numerical model of a metasurface: (**a**) model geometry and description and (**b**) computational mesh shown on the selected model boundaries.