

=============================================================================
MATLAB code for the following paper:

Topic: Scattering Problems in Periodic or Locally Perturbed Periodic Domains: 
       Convergence of the Perfectly Matched Layer Approximation 
       and Efficient Numerical Solution

Authors:  Tilo Arens, Nasim Shafieeabyaneh, Ruming Zhang in June 2026
==============================================================================
Requirements:

 * MATLAB_R2023b
 * Python 3.10.14
 * Gmsh 4.11.1
 * function lgwt.m (download from https://www.mathworks.com/matlabcentral/
    fileexchange/4540-legendre-gauss-quadrature-weights-and-nodes)
 * function QuadLTri.m (download from https://de.mathworks.com/matlabcentral/
   fileexchange/104305-gauss-legendre-quadrature-triangle/files/Upload/QuadLTri.m)
=============================================================================
                   Numerical results reported in Tables 1 and 2:
==============================================================================
Run the main script MainPML.m : 

>> inputs: meshsize: mesh size in the spatial direction
                     Note that you should choose: '04', '02' or '01'
             source: location of the point source, e.g., [0,0.5]
                  K: wave number 

>> Overview:

1) GenerateMesh_periodic.py
   # Generate a triangular mesh for the bounded cell
     Inputs: it is required to enter the meshsize and the definition of 
             periodic lower boundary in the computational domain
     Output: one .msh file is produced which shows the generated mesh
 
2) Run MainPML.m which contains the following functions:

Read_Mesh                             # Read the mesh file generated by Gmsh (from the previous part) 
Split_Integration             	      # Generate a mesh for alpha domain
System                                # Compute the coefficient matrix and right-hand side 
                                        (Computed by FEM) of the linear system given in equation (39) 
LeftHandSide.m                        # Compute the left-hand side of (39)
Solve_System_Schurcomplement.m        # Solve the systems on the left and right-hand sides of (39)
QuadraturePointsTriangles.m           # Select the Gaussian quadrature nodes in each triangle
PWL_function.m                        # Compute the piecewise linear basis functions 
PML_function.m                        # Compute the PML function
CoeffMatrixPeriodic.m                 # Compute the coefficient matrix corresponding to the periodic part
GlobalStiffnessMatricesPeriodic.m     # Assemble the global stiffness matrices
LocalStiffnessMatricesPeriodic.m      # Construct the local stiffness matrices for periodic part
LocalStiffnessMatricesPerturbed.m     # Construct the local stiffness matrices for perturbed part
B_Perturbed.m                         # Compute the coefficient matrix corresponding to the perturbed part
GradientTransformation.m              # Compute the coeficients A und c 
sinc.m                                # Define the sinc function
Greensfunction.m                      # Compute exact solution only when the source is located below the scatterer

=================================================================================================
                      Numerical solutions plotted in Figures 3 and 4:
=================================================================================================
To visualize the solution in the main (where the perturbation is) and neighboring cells:
 1) GenerateMeshPeriodic.py           # Generate a triangular mesh in the bounded cell
 2) MainPML.m                         # Obtain the numerical solution only in the main cell
 3) ExtendedSolution.m                # Extend the numerical solution to the neighboring cells 
     Output: SolutionExtended01paper_k 
             ErrorExtended01paper_k    
             (Note that the file name depends on the selected mesh size and k)
    Requirements:
    transformedFieldScatteredField.m  # Use the definition of the inverse FB transform to extend
                                        the numerical solutions  
    GreensfunctionExtended.m          # Extend the exact solution to the neighboring cells 
 4) plotSolution.py                   # Plot the numerical solution and its numerical error
    On line 7, add the path where the output of ExtendSolution.m is located
    On line 8, add the path where the plots should be placed
    Note that perturbMesh.py is an additional function for this file                     

