#include <malloc/_malloc.h>
#include <petscsystypes.h>
static char help[] = "Tests DMSwarm with DMShell\n\n";

#include <fstream>
#include <iostream>
#include <petsc/private/dmimpl.h>
#include <petscdm.h>
#include <petscdmda.h>
#include <petscdmshell.h>
#include <petscdmswarm.h>
#include <petscsf.h>
#include <petscvec.h>
#include <string>

#ifndef PETSC_SUCCESS
#define PETSC_SUCCESS 0
#endif

#define maxneigh 250
#define maxneigh10 1000
#define NumberDimensions 3
#define BufferLenght 4
#define occu_H_zero 5e-2
#define r_cutoff_ADP 6.2934034034 //! Angular dependant potential cutoff

static double dmax_arg1, dmax_arg2;
#define DMAX(a, b)                                                             \
  (dmax_arg1 = (a), dmax_arg2 = (b),                                           \
   (dmax_arg1) > (dmax_arg2) ? (dmax_arg1) : (dmax_arg2))

double box_x_min = 0.0;
double box_x_max = 0.0;
double box_y_min = 0.0;
double box_y_max = 0.0;
double box_z_min = 0.0;
double box_z_max = 0.0;

PetscInt ndiv_mesh_X;
PetscInt ndiv_mesh_Y;
PetscInt ndiv_mesh_Z;

double box_origin_0[3];

typedef struct DMD {

  /*! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
   @brief System information
    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !*/
  //! @param n_atoms Number of atoms
  int n_atoms;

  //! @param Variable: with the atomistic data
  DM atomistic_data;

  //! @param Variable: with the background mesh
  DM background_mesh;

  //! @param Variable: with the bounding cell
  DM bounding_cell;

  /*! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
   @brief Interatomic potential information
    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !*/
  //! @param r2_cutoff: Cutoff radius of the interatomic potential
  double r2_cutoff;

} DMD;

/*******************************************************/

typedef struct dump_file {

  /** @param n_atoms Number of atoms */
  int n_atoms;

  /** @param specie: Integer which defines the atomic specie (H = 1 and Mg = 0)
   */
  int *specie;

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    System Lagrange multipliers
    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  /** @param beta: Thermodynamic Lagrange multiplier */
  double *beta;

  /** @param gamma: Chemical Lagrange multiplier */
  double *gamma;

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    Position-related variables
    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  /** @param mean_q: Mean value of each atomic position */
  double *mean_q;

  /** @param stdv_q: Standard desviation of each atomic position */
  double *stdv_q;

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    Chemical variables
    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  /** @param xi: Molar fraction (mean occupancy) */
  double *xi;

} dump_file;

/*******************************************************/

dump_file read_dump_information(const char *SimulationFile) {

  unsigned int dim = NumberDimensions;

  dump_file Simulation;

  //! Variables to read atomistic data
  FILE *simulation_data = fopen(SimulationFile, "r");
  char aux_line[10000];
  int error;

  //! @brief Read the time step of the file
  int time_step;
  fgets(aux_line, sizeof(aux_line), simulation_data);
  error = fscanf(simulation_data, "%d\n", &time_step);

  //! @brief Read number of atomic position
  int n_atoms;
  fgets(aux_line, sizeof(aux_line), simulation_data);
  error = fscanf(simulation_data, "%d\n", &n_atoms);
  Simulation.n_atoms = n_atoms;

  //! @brief Read box bounds
  fgets(aux_line, sizeof(aux_line), simulation_data);
  error = fscanf(simulation_data, "%lf %lf\n", &box_x_min, &box_x_max);
  error = fscanf(simulation_data, "%lf %lf\n", &box_y_min, &box_y_max);
  error = fscanf(simulation_data, "%lf %lf\n", &box_z_min, &box_z_max);

  //! @brief Compute box-parameters
  box_origin_0[0] = box_x_min;
  box_origin_0[1] = box_y_min;
  box_origin_0[2] = box_z_min;

  //! @brief Create arrays
  Simulation.specie = (int *)calloc(n_atoms, sizeof(int));

  Simulation.beta = (double *)calloc(n_atoms, sizeof(double));

  Simulation.gamma = (double *)calloc(n_atoms, sizeof(double));

  Simulation.mean_q = (double *)calloc(n_atoms * dim, sizeof(double));

  Simulation.stdv_q = (double *)calloc(n_atoms, sizeof(double));

  Simulation.xi = (double *)calloc(n_atoms, sizeof(double));

  //! @brief Read properties from dump file
  fgets(aux_line, sizeof(aux_line), simulation_data);
  unsigned int Particle_Type;

  for (unsigned int i_site = 0; i_site < n_atoms; i_site++) {

    // clang-format off
    error = fscanf(simulation_data, "%i %lf %lf %lf %lf %lf %lf %lf\n",
                   &(Particle_Type), 
                   &(Simulation.mean_q[3 * i_site + 0]),
                   &(Simulation.mean_q[3 * i_site + 1]),
                   &(Simulation.mean_q[3 * i_site + 2]),
                   &(Simulation.stdv_q[i_site]), 
                   &(Simulation.xi[i_site]), 
                   &(Simulation.gamma[i_site]),
                   &(Simulation.beta[i_site]));
    // clang-format on

    // Initialize particle type.
    if (Particle_Type == 1) {
      Simulation.specie[i_site] = 0; //! Magnesium
    } else if (Particle_Type == 2) {
      Simulation.specie[i_site] = 1; //! Hydrogen
    } else {
      std::cout << "FAIL Particle_Type:" << i_site << std::endl;
    }
  }

  //  Close file
  fclose(simulation_data);

  // Check Hx molar fraction
  double xi_i;
  double xi_min_Hx = occu_H_zero;
  double xi_max_Hx = 1.0 - occu_H_zero;
  for (unsigned int i_site = 0; i_site < n_atoms; i_site++) {
    if (Simulation.specie[i_site] == 1) { // Hydrogen
      xi_i = Simulation.xi[i_site];

      Simulation.xi[i_site] = xi_i < xi_min_Hx ? xi_min_Hx : xi_i;

      Simulation.xi[i_site] = xi_i > xi_max_Hx ? xi_max_Hx : xi_i;
    }
  }

  return Simulation;
}

/*******************************************************/

PetscErrorCode _DMLocatePoints_DMDARegular_IS(DM dm, Vec pos, IS *iscell) {
  PetscInt p, n, bs, npoints, si, sj, sk, milocal, mjlocal, mklocal, mx, my, mz;
  DM background_mesh;
  PetscInt *cellidx;
  const PetscScalar *coor;
  PetscMPIInt rank;

  PetscFunctionBegin;
  PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
  PetscCall(VecGetLocalSize(pos, &n));
  PetscCall(VecGetBlockSize(pos, &bs));
  npoints = n / bs;

  PetscCall(PetscMalloc1(npoints, &cellidx));
  PetscCall(DMGetApplicationContext(dm, &background_mesh));
  PetscCall(DMDAGetCorners(background_mesh, &si, &sj, &sk, &milocal, &mjlocal,
                           &mklocal));
  PetscCall(DMDAGetInfo(background_mesh, NULL, &mx, &my, &mz, NULL, NULL, NULL,
                        NULL, NULL, NULL, NULL, NULL, NULL));

  double dx = (box_x_max - box_x_min) / ((PetscReal)mx - 1);
  double dy = (box_y_max - box_y_min) / ((PetscReal)my - 1);
  double dz = (box_z_max - box_z_min) / ((PetscReal)mz - 1);

  char name[PETSC_MAX_PATH_LEN];
  PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1,
                          "outputs/atoms-rank%d.csv", rank));

  std::ofstream file_f;
  file_f.open(name, std::ios::trunc);

  file_f << "rank: " << rank << std::endl;
  file_f << "si: " << si << " , "
         << "milocal: " << milocal << std::endl;
  file_f << "sj: " << sj << " , "
         << "mjlocal: " << mjlocal << std::endl;
  file_f << "sk: " << sk << " , "
         << "mklocal: " << mklocal << std::endl;

  file_f << "npoints: " << npoints << std::endl;

  file_f << "p"
         << " , "
         << "Coord.x"
         << " , "
         << "Coord.y"
         << " , "
         << "Coord.z"
         << " , "
         << "cell_idx"
         << " , "
         << "mi"
         << " , "
         << "mj"
         << " , "
         << "mk" << std::endl;

  PetscCall(VecGetArrayRead(pos, &coor));
  for (p = 0; p < npoints; p++) {
    PetscReal coorx, coory, coorz;
    PetscInt mi, mj, mk;

    coorx = PetscRealPart(coor[3 * p + 0]);
    coory = PetscRealPart(coor[3 * p + 1]);
    coorz = PetscRealPart(coor[3 * p + 2]);

    mi = (PetscInt)round((coorx - (box_x_min)) / dx);
    mj = (PetscInt)round((coory - (box_y_min)) / dy);
    mk = (PetscInt)round((coorz - (box_z_min)) / dz);

    cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND;

    if ((mk >= sk) && (mk < sk + mklocal)) {
      if ((mj >= sj) && (mj < sj + mjlocal)) {
        if ((mi >= si) && (mi < si + milocal))
          cellidx[p] =
              (mi - si) + (mj - sj) * milocal + (mk - sk) * milocal * mklocal;
      }
    }

    //  if (coorx < box_x_min)
    //    cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND;
    //  if (coorx > box_x_max)
    //    cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND;
    //  if (coory < box_y_min)
    //    cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND;
    //  if (coory > box_y_max)
    //    cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND;
    //  if (coorz < box_z_min)
    //    cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND;
    //  if (coorz > box_z_max)
    //    cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND;

    file_f << p << " , " << coorx << " , " << coory << " , " << coorz << " , "
           << cellidx[p] << " , " << mi << " , " << mj << " , " << mk
           << std::endl;
  }

  file_f.close();

  return 1;

  PetscCall(VecRestoreArrayRead(pos, &coor));
  PetscCall(ISCreateGeneral(PETSC_COMM_SELF, npoints, cellidx,
                            PETSC_OWN_POINTER, iscell));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*******************************************************/

PetscErrorCode DMLocatePoints_DMDARegular(DM dm, Vec pos,
                                          DMPointLocationType ltype,
                                          PetscSF cellSF) {
  IS iscell;
  PetscSFNode *cells;
  PetscInt p, bs, npoints, nfound;
  const PetscInt *boxCells;

  PetscFunctionBegin;
  PetscCall(_DMLocatePoints_DMDARegular_IS(dm, pos, &iscell));
  PetscCall(VecGetLocalSize(pos, &npoints));
  PetscCall(VecGetBlockSize(pos, &bs));
  npoints = npoints / bs;

  PetscCall(PetscMalloc1(npoints, &cells));
  PetscCall(ISGetIndices(iscell, &boxCells));

  for (p = 0; p < npoints; p++) {
    cells[p].rank = 0;
    cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
    cells[p].index = boxCells[p];
    std::cout << "site " << p << " is located in cell: " << boxCells[p]
              << std::endl;
  }

  PetscCall(ISRestoreIndices(iscell, &boxCells));
  PetscCall(ISDestroy(&iscell));
  nfound = npoints;
  PetscCall(PetscSFSetGraph(cellSF, npoints, nfound, NULL, PETSC_OWN_POINTER,
                            cells, PETSC_OWN_POINTER));
  PetscCall(ISDestroy(&iscell));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*******************************************************/

PetscErrorCode DMGetNeighbors_DMDARegular(DM dm, PetscInt *nneighbors,
                                          const PetscMPIInt **neighbors) {
  DM background_mesh;

  PetscFunctionBegin;
  PetscCall(DMGetApplicationContext(dm, &background_mesh));
  PetscCall(DMGetNeighbors(background_mesh, nneighbors, neighbors));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*******************************************************/

PetscErrorCode ParaviewViewGP(DM atomistic_data, const char prefix[]) {

  char name_vtk[PETSC_MAX_PATH_LEN];

  PetscFunctionBegin;

  PetscCall(PetscSNPrintf(name_vtk, PETSC_MAX_PATH_LEN - 1, "%s.xmf", prefix));

  PetscCall(DMSwarmViewXDMF(atomistic_data, name_vtk));

  PetscFunctionReturn(PETSC_SUCCESS);
}

/*******************************************************/

/*
 Create a DMShell and attach a regularly spaced DMDA for point location
 Override methods for point location
*/
PetscErrorCode init_DMD_simulation(DMD *Simulation, dump_file Simulation_file) {

  unsigned int dim = NumberDimensions;

  DM atomistic_data, bounding_cell, background_mesh;
  PetscMPIInt rank;
  PetscInt blocksize;

  //! Topological variables
  AO dump2petsc_mapping;

  //!
  PetscInt n_atoms_local;
  PetscInt n_atoms_global;
  PetscInt n_atoms = Simulation_file.n_atoms;

  //!
  IS beta_bcc;
  IS gamma_bcc;
  IS interstitial_sites;

  PetscFunctionBegin;
  PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));

  /* Create a regularly spaced DMDA */
  PetscInt overlap = 1;
  PetscInt dof = 1;
  PetscCall(DMDACreate3d(PETSC_COMM_WORLD,                                  //!
                         DM_BOUNDARY_PERIODIC,                              //!
                         DM_BOUNDARY_PERIODIC,                              //!
                         DM_BOUNDARY_PERIODIC,                              //!
                         DMDA_STENCIL_BOX,                                  //!
                         ndiv_mesh_X, ndiv_mesh_Y, ndiv_mesh_Z,             //!
                         PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE, //!
                         dof, overlap, NULL, NULL, NULL, &background_mesh));

  PetscCall(DMSetFromOptions(background_mesh));

  PetscCall(DMSetUp(background_mesh));

  PetscCall(DMDASetUniformCoordinates(background_mesh,
                                      box_x_min, //!
                                      box_x_max, //!
                                      box_y_min, //!
                                      box_y_max, //!
                                      box_z_min, //!
                                      box_z_max));

  /* Create a DMShell for point location purposes */
  PetscCall(DMShellCreate(PETSC_COMM_WORLD, &bounding_cell));
  PetscCall(DMSetApplicationContext(bounding_cell, background_mesh));
  bounding_cell->ops->locatepoints = DMLocatePoints_DMDARegular;
  bounding_cell->ops->getneighbors = DMGetNeighbors_DMDARegular;

  /* Create the swarm */
  PetscCall(DMCreate(PETSC_COMM_WORLD, &atomistic_data));
  PetscCall(DMSetType(atomistic_data, DMSWARM));
  PetscCall(DMSetDimension(atomistic_data, dim));
  PetscCall(PetscObjectSetName((PetscObject)atomistic_data, "Atoms"));

  PetscCall(DMSwarmSetType(atomistic_data, DMSWARM_PIC));
  PetscCall(DMSwarmSetCellDM(atomistic_data, bounding_cell));

  //! Site idx
  PetscCall(
      DMSwarmRegisterPetscDatatypeField(atomistic_data, "idx", 1, PETSC_INT));

  //! MPI rank
  PetscCall(DMSwarmRegisterPetscDatatypeField(atomistic_data, "MPI-rank", 1,
                                              PETSC_INT));

  //! Molar fraction: Xi := <n>_0
  PetscCall(DMSwarmRegisterPetscDatatypeField(atomistic_data, "molar-fraction",
                                              1, PETSC_REAL));

  //! Energy density: mf_rho := <rho>_0
  PetscCall(DMSwarmRegisterPetscDatatypeField(atomistic_data, "mf-rho", 1,
                                              PETSC_REAL));

  //! Thermal multiplier: Beta := 1/(kb * T)
  PetscCall(
      DMSwarmRegisterPetscDatatypeField(atomistic_data, "beta", 1, PETSC_REAL));

  //! Chemical multiplier: Gamma
  PetscCall(DMSwarmRegisterPetscDatatypeField(atomistic_data, "gamma", 1,
                                              PETSC_REAL));

  //! Standard desviation of the position: stdv-q
  PetscCall(DMSwarmRegisterPetscDatatypeField(atomistic_data, "stdv-q", 1,
                                              PETSC_REAL));

  //! First Piola-Kirchoff stress tensor: FPK
  PetscCall(
      DMSwarmRegisterPetscDatatypeField(atomistic_data, "FPK", 9, PETSC_REAL));

  //! Boundary condition idex
  PetscCall(DMSwarmRegisterPetscDatatypeField(atomistic_data, "idx-bcc-beta", 1,
                                              PETSC_INT));
  PetscCall(DMSwarmRegisterPetscDatatypeField(atomistic_data, "idx-bcc-gamma",
                                              1, PETSC_INT));
  PetscCall(DMSwarmRegisterPetscDatatypeField(atomistic_data, "idx-bcc-fix", 1,
                                              PETSC_INT));

  PetscCall(DMSwarmFinalizeFieldRegister(atomistic_data));
  PetscCall(DMSwarmSetLocalSizes(atomistic_data, n_atoms, BufferLenght));

  //! @brief Get the bounding box of the background mesh
  PetscReal local_domain_ll[3];
  PetscReal local_domain_ur[3];
  PetscCall(
      DMGetLocalBoundingBox(background_mesh, local_domain_ll, local_domain_ur));

  //! Lower left coordinates of the brick
  PetscReal xI_lw, yI_lw, zI_lw;
  xI_lw = local_domain_ll[0];
  yI_lw = local_domain_ll[1];
  zI_lw = local_domain_ll[2];

  //! Upper right coordinates of the brick
  PetscReal xI_up, yI_up, zI_up;
  xI_up = local_domain_ur[0];
  yI_up = local_domain_ur[1];
  zI_up = local_domain_ur[2];

  //!
  PetscReal local_dx, local_dy, local_dz;
  local_dx = xI_up - xI_lw;
  local_dy = yI_up - yI_lw;
  local_dz = zI_up - zI_lw;

  //! @brief mean_q: Mean value of each atomic position
  double *mean_q_ptr;
  PetscCall(DMSwarmGetField(atomistic_data, DMSwarmPICField_coor, &blocksize,
                            NULL, (void **)&mean_q_ptr));

  //! @brief idx_ptr: Global index of each atomic position
  int *idx_ptr;
  PetscCall(
      DMSwarmGetField(atomistic_data, "idx", NULL, NULL, (void **)&idx_ptr));

  //! @brief Loop in the elemnts of the background mesh and find atoms in the
  //! mesh
  n_atoms_local = 0;
  for (PetscInt site_i = 0; site_i < n_atoms; site_i++) {

    if ((Simulation_file.mean_q[site_i * dim + 0] >= xI_lw) && //!
        (Simulation_file.mean_q[site_i * dim + 0] < xI_up) &&  //!
        (Simulation_file.mean_q[site_i * dim + 1] >= yI_lw) && //!
        (Simulation_file.mean_q[site_i * dim + 1] < yI_up) &&  //!
        (Simulation_file.mean_q[site_i * dim + 2] >= zI_lw) && //!
        (Simulation_file.mean_q[site_i * dim + 2] < zI_up)) {

      idx_ptr[n_atoms_local] = site_i;

      mean_q_ptr[n_atoms_local * dim + 0] =
          Simulation_file.mean_q[site_i * dim + 0];
      mean_q_ptr[n_atoms_local * dim + 1] =
          Simulation_file.mean_q[site_i * dim + 1];
      mean_q_ptr[n_atoms_local * dim + 2] =
          Simulation_file.mean_q[site_i * dim + 2];

      n_atoms_local++;
    }
  }

  PetscCall(DMSwarmRestoreField(atomistic_data, DMSwarmPICField_coor,
                                &blocksize, NULL, (void **)&mean_q_ptr));
  PetscCall(DMSwarmRestoreField(atomistic_data, "idx", NULL, NULL,
                                (void **)&idx_ptr));

  //! @brief Set local sizes
  PetscCall(DMSwarmSetLocalSizes(atomistic_data, n_atoms_local, BufferLenght));

  //! @brief idx_ptr: index of the site
  PetscCall(DMSwarmGetField(atomistic_data, "idx", &blocksize, NULL,
                            (void **)&idx_ptr));

  //! @brief specie: Integer which defines MPI rank location
  PetscInt *site_mpi_rank;
  PetscCall(DMSwarmGetField(atomistic_data, "MPI-rank", &blocksize, NULL,
                            (void **)&site_mpi_rank));

  //! @brief stdv_q: Standard desviation of each atomic position
  double *stdv_q;
  PetscCall(DMSwarmGetField(atomistic_data, "stdv-q", &blocksize, NULL,
                            (void **)&stdv_q));

  //! @brief mf-rho: Meanfield energy density
  double *mf_rho;
  PetscCall(DMSwarmGetField(atomistic_data, "mf-rho", &blocksize, NULL,
                            (void **)&mf_rho));

  //! @brief xi: Molar fraction (mean occupancy)
  double *xi;
  PetscCall(DMSwarmGetField(atomistic_data, "molar-fraction", &blocksize, NULL,
                            (void **)&xi));

  //! @brief beta: Thermal Lagrange multiplier
  double *beta;
  PetscCall(DMSwarmGetField(atomistic_data, "beta", &blocksize, NULL,
                            (void **)&beta));

  //! @brief gamma: Chemical Lagrange multiplier
  double *gamma;
  PetscCall(DMSwarmGetField(atomistic_data, "gamma", &blocksize, NULL,
                            (void **)&gamma));

  //! @brief Pointers with the information of the boundary condition
  PetscInt *beta_bcc_ptr;
  PetscCall(DMSwarmGetField(atomistic_data, "idx-bcc-beta", &blocksize, NULL,
                            (void **)&beta_bcc_ptr));

  PetscInt *gamma_bcc_ptr;
  PetscCall(DMSwarmGetField(atomistic_data, "idx-bcc-gamma", &blocksize, NULL,
                            (void **)&gamma_bcc_ptr));

  //! @brief Pointer with the list of interstitial sites
  PetscInt *interstitial_ptr =
      (PetscInt *)calloc(n_atoms_local, sizeof(PetscInt));

  //! Set information from the .dump file
  for (PetscInt local_site_i = 0; local_site_i < n_atoms_local;
       local_site_i++) {
    PetscInt site_i = idx_ptr[local_site_i];
    site_mpi_rank[local_site_i] = (PetscInt)rank;
    stdv_q[local_site_i] = Simulation_file.stdv_q[site_i];
    xi[local_site_i] = Simulation_file.xi[site_i];
    mf_rho[local_site_i] = 0.0;
    beta[local_site_i] = Simulation_file.beta[site_i];
    gamma[local_site_i] = Simulation_file.gamma[site_i];
  }

  //! Create a context to relate petsc numbering and dump numbering
  IS dump_ordering;
  PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, n_atoms_local, idx_ptr,
                            PETSC_COPY_VALUES, &dump_ordering));
  PetscCall(AOCreateBasicIS(dump_ordering, NULL, &dump2petsc_mapping));
  PetscCall(ISDestroy(&dump_ordering));

  //! Tunrn dump numbering of "idx", and
  //! "interstitial" variables into petsc numbering
  PetscCall(AOApplicationToPetsc(dump2petsc_mapping, n_atoms_local, idx_ptr));
  PetscCall(AOApplicationToPetsc(dump2petsc_mapping, n_atoms_local,
                                 interstitial_ptr));

  //! Structures with the list of interstitial sites
  PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, n_atoms_local, interstitial_ptr,
                            PETSC_COPY_VALUES, &interstitial_sites));
  PetscCall(PetscObjectSetName((PetscObject)interstitial_sites,
                               "interstitial-sites"));

  //! Remove negative index
  PetscCall(ISGeneralFilter(interstitial_sites, 0, n_atoms));

  //! Restore particle fields
  PetscCall(DMSwarmRestoreField(atomistic_data, "idx", &blocksize, NULL,
                                (void **)&idx_ptr));
  PetscCall(DMSwarmRestoreField(atomistic_data, "MPI-rank", &blocksize, NULL,
                                (void **)&site_mpi_rank));
  PetscCall(DMSwarmRestoreField(atomistic_data, "stdv-q", &blocksize, NULL,
                                (void **)&stdv_q));
  PetscCall(DMSwarmRestoreField(atomistic_data, "mf-rho", &blocksize, NULL,
                                (void **)&mf_rho));
  PetscCall(DMSwarmRestoreField(atomistic_data, "molar-fraction", &blocksize,
                                NULL, (void **)&xi));
  PetscCall(DMSwarmRestoreField(atomistic_data, "beta", &blocksize, NULL,
                                (void **)&beta));
  PetscCall(DMSwarmRestoreField(atomistic_data, "gamma", &blocksize, NULL,
                                (void **)&gamma));
  PetscCall(DMSwarmRestoreField(atomistic_data, "idx-bcc-beta", &blocksize,
                                NULL, (void **)&beta_bcc_ptr));
  PetscCall(DMSwarmRestoreField(atomistic_data, "idx-bcc-gamma", &blocksize,
                                NULL, (void **)&gamma_bcc_ptr));

  free(interstitial_ptr);

  //! Check global size
  DMSwarmGetSize(atomistic_data, &n_atoms_global);
  if (n_atoms_global != Simulation_file.n_atoms) {
    PetscCall(PetscError(
        PETSC_COMM_WORLD, __LINE__, "init_DMD_simulation", __FILE__,
        PETSC_ERR_RETURN, PETSC_ERROR_INITIAL,
        "The global size (%i) of the DMSwarm does not match the input (%i) !",
        n_atoms_global, Simulation_file.n_atoms));
    PetscFunctionReturn(PETSC_ERR_RETURN);
  }

  //! Generate output structure
  Simulation->atomistic_data = atomistic_data;
  Simulation->background_mesh = background_mesh;
  Simulation->bounding_cell = bounding_cell;

  PetscFunctionReturn(PETSC_SUCCESS);
}

/*******************************************************/

PetscErrorCode destroy_DMD_simulation(DMD *Simulation) {

  //! @brief Atomistic variables
  DMDestroy(&Simulation->atomistic_data);
  DMDestroy(&Simulation->background_mesh);
  DMDestroy(&Simulation->bounding_cell);

  PetscFunctionReturn(PETSC_SUCCESS);
}

/*******************************************************/

int main(int argc, char **argv) {
  PetscFunctionBeginUser;
  PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));

  //  const char SimulationFile[10000] = "Mg-hcp-cube-x3-x3-x3.dump";
  //  const char SimulationFile[10000] = "Mg-hcp-cube-x5-x5-x5.dump";
  const char SimulationFile[10000] = "Mg-hcp-cube-x17-x10-x10.dump";
  //  const char SimulationFile[10000] = "Mg-sphere-radius-20.dump";
  //  const char SimulationFile[10000] = "MgHx-hcp-hex-nanowire-D-40-A.dump";

  dump_file Simulation_file = read_dump_information(SimulationFile);

  //!

  ndiv_mesh_X = 3;
  ndiv_mesh_Y = 3;
  ndiv_mesh_Z = 3;

  DMD Simulation;
  PetscCall(init_DMD_simulation(&Simulation, Simulation_file));

  PetscCall(DMView(Simulation.atomistic_data, PETSC_VIEWER_STDOUT_WORLD));

  /* Get outputs */
  //  PetscCall(DMView(Simulation.atomistic_data, PETSC_VIEWER_STDOUT_WORLD));
  //  PetscCall(DMView(Simulation.background_mesh, PETSC_VIEWER_STDOUT_WORLD));
  // PetscCall(DMView(Simulation.bounding_cell, PETSC_VIEWER_STDOUT_WORLD));
  PetscCall(ParaviewViewGP(Simulation.atomistic_data, "outputs/atoms"));

  // Export the mesh to VTK format (VTU)
  Vec mesh_coors;
  DMGetCoordinates(Simulation.background_mesh, &mesh_coors);
  PetscViewer viewer;
  PetscViewerVTKOpen(PETSC_COMM_WORLD, "outputs/mesh.vts", FILE_MODE_WRITE,
                     &viewer);
  VecView(mesh_coors, viewer);
  PetscViewerDestroy(&viewer);

  /* Free data */
  PetscCall(destroy_DMD_simulation(&Simulation));

  PetscCall(PetscFinalize());
  return 0;
}

/*TEST

   build:
      requires: double !complex

   test:
      filter: grep -v atomic
      filter_output: grep -v atomic
TEST*/