!> Input/Output routines for native hexahedral mesh files.
!!
!! This module handles the parsing of the ASCII-based native mesh format 
!! used by the solver. It reads node coordinates, cell connectivity, 
!! face metrics, and boundary patches from the directory specified 
!! in the case parameters.
module mod_mesh_io
   use mod_kinds, only : rk, path_len, fatal_error
   use mod_mesh_types, only : mesh_t, mesh_finalize, max_cell_faces
   implicit none

   private

   public :: read_native_mesh

contains

   !> Orchestrates the reading of all mesh components from a directory.
   !!
   !! @param mesh_dir Path to the directory containing `.dat` files.
   !! @param m The mesh structure to populate.
   subroutine read_native_mesh(mesh_dir, m)
      character(len=*), intent(in) :: mesh_dir
      type(mesh_t), intent(inout) :: m

      call mesh_finalize(m)
      call read_points(trim(mesh_dir)//'/points.dat', m)
      call read_cells(trim(mesh_dir)//'/cells.dat', m)
      call read_faces(trim(mesh_dir)//'/faces.dat', m)
      call read_patches(trim(mesh_dir)//'/patches.dat', m)
      call read_periodic_optional(trim(mesh_dir)//'/periodic.dat', m)
      call build_cell_faces(m)
   end subroutine read_native_mesh


   !> Reads node coordinates from `points.dat`.
   subroutine read_points(filename, m)
      character(len=*), intent(in) :: filename
      type(mesh_t), intent(inout) :: m

      integer :: unit_id, ios, i, id
      real(rk) :: x, y, z

      open(newunit=unit_id, file=trim(filename), status='old', action='read', iostat=ios)
      if (ios /= 0) call fatal_error('mesh_io', 'could not open '//trim(filename))

      read(unit_id, *, iostat=ios) m%npoints
      if (ios /= 0 .or. m%npoints <= 0) call fatal_error('mesh_io', 'invalid points header')

      allocate(m%points(3, m%npoints))
      do i = 1, m%npoints
         read(unit_id, *, iostat=ios) id, x, y, z
         if (ios /= 0) call fatal_error('mesh_io', 'failed reading point')
         if (id < 1 .or. id > m%npoints) call fatal_error('mesh_io', 'point id out of range')
         m%points(:, id) = [x, y, z]
      end do

      close(unit_id)
   end subroutine read_points


   !> Reads cell definitions from `cells.dat`.
   subroutine read_cells(filename, m)
      character(len=*), intent(in) :: filename
      type(mesh_t), intent(inout) :: m

      integer :: unit_id, ios, i, id
      integer :: nodes(8)
      real(rk) :: cx, cy, cz, volume

      open(newunit=unit_id, file=trim(filename), status='old', action='read', iostat=ios)
      if (ios /= 0) call fatal_error('mesh_io', 'could not open '//trim(filename))

      read(unit_id, *, iostat=ios) m%ncells
      if (ios /= 0 .or. m%ncells <= 0) call fatal_error('mesh_io', 'invalid cells header')

      allocate(m%cells(m%ncells))
      do i = 1, m%ncells
         read(unit_id, *, iostat=ios) id, nodes, cx, cy, cz, volume
         if (ios /= 0) call fatal_error('mesh_io', 'failed reading cell')
         if (id < 1 .or. id > m%ncells) call fatal_error('mesh_io', 'cell id out of range')
         if (volume <= 0.0_rk) call fatal_error('mesh_io', 'cell volume must be positive')
         m%cells(id)%id = id
         m%cells(id)%nodes = nodes
         m%cells(id)%center = [cx, cy, cz]
         m%cells(id)%volume = volume
      end do

      close(unit_id)
   end subroutine read_cells


   !> Reads face metrics and owner/neighbor connectivity from `faces.dat`.
   subroutine read_faces(filename, m)
      character(len=*), intent(in) :: filename
      type(mesh_t), intent(inout) :: m

      integer :: unit_id, ios, i, id
      integer :: owner, neighbor, patch
      real(rk) :: nx, ny, nz, area, cx, cy, cz

      open(newunit=unit_id, file=trim(filename), status='old', action='read', iostat=ios)
      if (ios /= 0) call fatal_error('mesh_io', 'could not open '//trim(filename))

      read(unit_id, *, iostat=ios) m%nfaces
      if (ios /= 0 .or. m%nfaces <= 0) call fatal_error('mesh_io', 'invalid faces header')

      allocate(m%faces(m%nfaces))
      do i = 1, m%nfaces
         read(unit_id, *, iostat=ios) id, owner, neighbor, patch, nx, ny, nz, area, cx, cy, cz
         if (ios /= 0) call fatal_error('mesh_io', 'failed reading face')
         if (id < 1 .or. id > m%nfaces) call fatal_error('mesh_io', 'face id out of range')
         if (owner < 1 .or. owner > m%ncells) call fatal_error('mesh_io', 'face owner out of range')
         if (neighbor < 0 .or. neighbor > m%ncells) call fatal_error('mesh_io', 'face neighbor out of range')
         if (area <= 0.0_rk) call fatal_error('mesh_io', 'face area must be positive')
         m%faces(id)%id = id
         m%faces(id)%owner = owner
         m%faces(id)%neighbor = neighbor
         m%faces(id)%patch = patch
         m%faces(id)%normal = [nx, ny, nz]
         m%faces(id)%area = area
         m%faces(id)%center = [cx, cy, cz]
      end do

      close(unit_id)
   end subroutine read_faces


   !> Reads boundary patch definitions and face lists from `patches.dat`.
   subroutine read_patches(filename, m)
      character(len=*), intent(in) :: filename
      type(mesh_t), intent(inout) :: m

      integer :: unit_id, ios, p, id, nfaces
      character(len=64) :: name

      open(newunit=unit_id, file=trim(filename), status='old', action='read', iostat=ios)
      if (ios /= 0) call fatal_error('mesh_io', 'could not open '//trim(filename))

      read(unit_id, *, iostat=ios) m%npatches
      if (ios /= 0 .or. m%npatches < 0) call fatal_error('mesh_io', 'invalid patches header')

      allocate(m%patches(m%npatches))
      do p = 1, m%npatches
         read(unit_id, *, iostat=ios) id, name, nfaces
         if (ios /= 0) call fatal_error('mesh_io', 'failed reading patch header')
         if (id < 1 .or. id > m%npatches) call fatal_error('mesh_io', 'patch id out of range')
         m%patches(id)%id = id
         m%patches(id)%name = name
         m%patches(id)%nfaces = nfaces
         allocate(m%patches(id)%face_ids(nfaces))
         if (nfaces > 0) then
            read(unit_id, *, iostat=ios) m%patches(id)%face_ids
            if (ios /= 0) call fatal_error('mesh_io', 'failed reading patch face ids')
         else
            read(unit_id, *, iostat=ios)
         end if
      end do

      close(unit_id)
   end subroutine read_patches


   !> Reads periodic link data from `periodic.dat` if it exists.
   subroutine read_periodic_optional(filename, m)
      character(len=*), intent(in) :: filename
      type(mesh_t), intent(inout) :: m

      integer :: unit_id, ios, nlinks, i
      integer :: face_id, pair_face_id, neighbor_cell

      open(newunit=unit_id, file=trim(filename), status='old', action='read', iostat=ios)
      if (ios /= 0) return

      read(unit_id, *, iostat=ios) nlinks
      if (ios /= 0 .or. nlinks < 0) call fatal_error('mesh_io', 'invalid periodic header')

      do i = 1, nlinks
         read(unit_id, *, iostat=ios) face_id, pair_face_id, neighbor_cell
         if (ios /= 0) call fatal_error('mesh_io', 'failed reading periodic link')
         if (face_id < 1 .or. face_id > m%nfaces) call fatal_error('mesh_io', 'periodic face id out of range')
         if (pair_face_id < 1 .or. pair_face_id > m%nfaces) call fatal_error('mesh_io', 'periodic pair face id out of range')
         if (neighbor_cell < 1 .or. neighbor_cell > m%ncells) call fatal_error('mesh_io', 'periodic neighbor cell out of range')
         m%faces(face_id)%periodic_face = pair_face_id
         m%faces(face_id)%periodic_neighbor = neighbor_cell
      end do

      close(unit_id)
   end subroutine read_periodic_optional


   !> Builds the cell-to-face mapping by inverting the face owner/neighbor data.
   subroutine build_cell_faces(m)
      type(mesh_t), intent(inout) :: m

      integer :: f, c

      allocate(m%ncell_faces(m%ncells))
      allocate(m%cell_faces(max_cell_faces, m%ncells))
      m%ncell_faces = 0
      m%cell_faces = 0

      do f = 1, m%nfaces
         c = m%faces(f)%owner
         call append_cell_face(m, c, f)

         c = m%faces(f)%neighbor
         if (c > 0) call append_cell_face(m, c, f)
      end do

      do c = 1, m%ncells
         if (m%ncell_faces(c) /= max_cell_faces) then
            call fatal_error('mesh_io', 'each v1 cuboid cell must have exactly six faces')
         end if
      end do
   end subroutine build_cell_faces


   !> Helper to safely add a face to a cell's local face list.
   subroutine append_cell_face(m, cell_id, face_id)
      type(mesh_t), intent(inout) :: m
      integer, intent(in) :: cell_id
      integer, intent(in) :: face_id

      if (m%ncell_faces(cell_id) >= max_cell_faces) then
         call fatal_error('mesh_io', 'cell has more than six faces')
      end if

      m%ncell_faces(cell_id) = m%ncell_faces(cell_id) + 1
      m%cell_faces(m%ncell_faces(cell_id), cell_id) = face_id
   end subroutine append_cell_face

end module mod_mesh_io

