The fractal interpolation techniques are powerful alternatives to classical interpolation methods in case of complex irregular data such as financial time series, seismic data and biological signals. The main goal of this study is to apply fractal interpolation techniques over two-dimensional irregular meshes. We partition the polygonal interpolation domain into triangular regions and consider affine contraction mappings from triangular regions to triangular patches of the mesh. Calculation of the fractal interpolation coefficients is reduced to solving a linear system of equations. Freely chosen scaling variables of the iterated function system transformations provide flexibility. Additionally, the constructed iterated function system coefficients for the interpolation can be used directly to evaluate two-dimensional numerical integrals over the domain. We provide basic error results via numerical simulations.