/**
 * @file matrix3d_boundary_traversal.cpp
 *
 * @author Johannes Sjolund (j.sjolund@gmail.com)
 *
 * @brief Traverse the boundary elements on the faces of a 3D matrix, by
 * iterating sequentially over integers 0 to the sum of the numbers of
 * elements on each face (2 * (nx * ny + nx * nz + ny * nz)) where matrix
 * dimensions nx >= ny >= nz.
 *
 * @copyright Copyright 2019 Johannes Sjolund. All rights reserved.
 * @license This project is released under the GNU Public License.
 *
 */
#include <iomanip>
#include <iostream>

//! Defines indexing for flattened 3D matrix
#define I3D(x, y, z, nx, ny, nz) ((x) + (y) * (nx) + (z) * (nx) * (ny))

/**
 * @brief Calculate the coordinate of a boundary element on the faces of a 3D
 * matrix from an index between 0 and the sum of the numbers of elements on each
 * face. Matrix dimensions must be nx >= ny >= nz. Note that corners and edges
 * are repeated.
 *
 * @param i Index in range [0, 2 * (nx * ny + nx * nz + ny * nz) )
 * @param nx Matrix length on X-axis, range [0, ?]
 * @param ny Matrix length on Y-axis, range [0, nx]
 * @param nz Matrix length on Z-axis, range [0, ny]
 * @param x Matrix boundary face coordinate on X-axis
 * @param y Matrix boundary face coordinate on Y-axis
 * @param z Matrix boundary face coordinate on Z-axis
 */
void getBoundaryElement(int i, int nx, int ny, int nz, int *x, int *y, int *z) {
  int n = nx * ny + nx * nz + ny * nz;
  int j = i % n;
  int xy = (2 * nx * ny) / (j + nx * ny + 1);
  int xy_xz = (2 * nx * ny + 2 * nx * nz) / (j + nx * ny + nx * nz + 1);
  int xy_xz_yz = (2 * n) / (j + n + 1);
  int xz = xy_xz - xy;
  int yz = xy_xz_yz - xy_xz;
  int xy_x = j % nx;
  int xy_y = j / nx;
  int xy_z = (nz - 1) * (i / n);
  int xz_x = j % nx;
  int xz_y = (ny - 1) * (i / n);
  int xz_z = (j % (nx * ny)) / nx;
  int yz_x = (nx - 1) * (i / n);
  int yz_y = (j % (nx * ny + nx * nz)) % ny;
  int yz_z = (j % (nx * ny + nx * nz)) / ny;
  *x = xy * xy_x + xz * xz_x + yz * yz_x;
  *y = xy * xy_y + xz * xz_y + yz * yz_y;
  *z = xy * xy_z + xz * xz_z + yz * yz_z;
}

/**
 * @brief Fill a matrix with numbers > 0
 */
void fillMatrix(int *matrix, int nx, int ny, int nz) {
  for (int x = 0; x < nx; x++)
    for (int y = 0; y < ny; y++)
      for (int z = 0; z < nz; z++) {
        int i = I3D(x, y, z, nx, ny, nz);
        matrix[i] = i + 1;
      }
}

/**
 * @brief Check that all boundary face elements of the matrix are equal to zero,
 * while all others are not equal to zero
 */
bool verifyMatrix(int *matrix, int nx, int ny, int nz) {
  for (int x = 0; x < nx; x++)
    for (int y = 0; y < ny; y++)
      for (int z = 0; z < nz; z++) {
        int val = matrix[I3D(x, y, z, nx, ny, nz)];
        if (x == 0 || x == nx - 1 || y == 0 || y == ny - 1 || z == 0 ||
            z == nz - 1) {
          if (val != 0) return false;
        } else {
          if (val == 0) return false;
        }
      }
  return true;
}

/**
 * @brief Print the contents of a matrix to standard output
 */
void printMatrix(int *matrix, int nx, int ny, int nz) {
  for (int z = 0; z < nz; z++) {
    for (int y = 0; y < ny; y++) {
      for (int x = 0; x < nx; x++) {
        std::cout << std::setfill('0') << std::setw(3)
                  << matrix[I3D(x, y, z, nx, ny, nz)] << " ";
      }
      std::cout << std::endl;
    }
    std::cout << std::endl;
  }
}

/**
 * @brief Test the matrix boundary element iteration with different size
 * matrices and check the result.
 */
int main(int argc, char **argv) {
  int n_min = 1;
  int n_max = 50;

  // Matrix dimensions must be nx >= ny >= nz
  for (int nx = n_min; nx <= n_max; nx++)
    for (int ny = n_min; ny <= nx; ny++)
      for (int nz = n_min; nz <= ny; nz++) {
        // 3D matrix flattened into 1D array
        int matrix[nx * ny * nz];
        // Fill matrix with numbers except 0
        fillMatrix(matrix, nx, ny, nz);

        // Set all matrix boundary elements to zero
        const int n = 2 * (nx * ny + nx * nz + ny * nz);
        for (int i = 0; i < n; i++) {
          int x, y, z;
          getBoundaryElement(i, nx, ny, nz, &x, &y, &z);
          matrix[I3D(x, y, z, nx, ny, nz)] = 0;
        }

        // Verify that the matrix has the correct content
        if (!verifyMatrix(matrix, nx, ny, nz)) {
          std::cout << "FAIL: "
                    << "nx=" << nx << ", ny=" << ny << ", nz=" << nz
                    << std::endl;
          printMatrix(matrix, nx, ny, nz);
          return 1;
        } else {
          std::cout << "OK: "
                    << "nx=" << nx << ", ny=" << ny << ", nz=" << nz
                    << std::endl;
        }
      }
  return 0;
}