!---------------------------------------------------------------------------
! pmlib_mesh_map.f
!   version: 2D single-core
!---------------------------------------------------------------------------
!                  __    __
!    ___  _ _     / /__ / /
!   | _ \/ ` \ _ / // // _ \
!   | ,_/_,_,_\  \_/\_/\___/
!   |_|                      Mads M. Hejlesen* and Jens H. Walther*^
!                            * Technical University of Denmark
!                            ^ ETH-Zurich
!
!---------------------------------------------------------------------------
! Copyright (c) 2014 Technical University of Denmark, 
!                    Department of Mechanical Engineering
!
! This program is free software; you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation; either version 2 of the License, or
! (at your option) any later version.
! 
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with this program; if not, see <http://www.gnu.org/licenses/>
!
! Contact: mmhej@mek.dtu.dk or jhwa@mek.dtu.dk
!---------------------------------------------------------------------------
SUBROUTINE pmlib_mesh_map(topo_send,topo_recieve,nvar_send,ierr)

!---------------------------------------------------------------------------
! Sets up the inter topology mesh communication
!  - Finds the overlapping mesh indices
!  - Finds the communication sequence
!  - Allocates the send and recieve buffers
!
! The communication is calculated on processor 0 only
! after which it is broadcasted to the remaining processors.
!---------------------------------------------------------------------------

USE pmlib_mod_communication
USE pmlib_mod_topology

IMPLICIT NONE
!---------------------------------------------------------------------------
! Arguments
!---------------------------------------------------------------------------
  TYPE(class_topology)                      :: topo_send
  TYPE(class_topology)                      :: topo_recieve
  INTEGER, INTENT(IN)                       :: nvar_send
  INTEGER, INTENT(OUT)                      :: ierr

!---------------------------------------------------------------------------
! Local variables
!---------------------------------------------------------------------------
  CHARACTER(LEN=14)                     :: caller = 'pmlib_mesh_map'
  INTEGER                               :: ncomm, icomm
  INTEGER                               :: iproc, jproc
  INTEGER                               :: i,j,k

  REAL(MK), DIMENSION(pmlib_ndim)       :: xmin_send, xmax_send
  REAL(MK), DIMENSION(pmlib_ndim)       :: xmin_recieve, xmax_recieve
  REAL(MK), DIMENSION(pmlib_ndim)       :: dx
  INTEGER,DIMENSION(2*pmlib_ndim)       :: nghost_recieve, nghost_send
  INTEGER,DIMENSION(2*pmlib_ndim)       :: bc_send, bc_recieve
  INTEGER                               :: nghost
  INTEGER, DIMENSION(pmlib_ndim)        :: ncell_send, ncell_recieve


!---------------------------------------------------------------------------
! Initiate subroutine
!---------------------------------------------------------------------------
  ierr = 0

!---------------------------------------------------------------------------
! Allocate communication lists
!---------------------------------------------------------------------------
  ncomm               = 1
  communication%ncomm = 1
  communication%ctype = 1

  ALLOCATE( communication%imin_send(pmlib_ndim,ncomm) )
  ALLOCATE( communication%imax_send(pmlib_ndim,ncomm) )
  ALLOCATE( communication%imin_recieve(pmlib_ndim,ncomm) )
  ALLOCATE( communication%imax_recieve(pmlib_ndim,ncomm) )

  ALLOCATE( communication%ncell_send(pmlib_ndim,ncomm) )
  ALLOCATE( communication%ncell_recieve(pmlib_ndim,ncomm) )

!---------------------------------------------------------------------------
! Subdomain parameters
!---------------------------------------------------------------------------
  xmin_send      = topo_send%xmin
  xmax_send      = topo_send%xmax

  xmin_recieve   = topo_recieve%xmin
  xmax_recieve   = topo_recieve%xmax

  ncell_send     = topo_send%ncell
  ncell_recieve  = topo_recieve%ncell

  nghost_send    = topo_send%nghost
  nghost_recieve = topo_recieve%nghost

  bc_send        = topo_send%bound_cond
  bc_recieve     = topo_recieve%bound_cond

  dx = topo_recieve%dx

!---------------------------------------------------------------------------
! Check overlap
!---------------------------------------------------------------------------
  IF( NINT( (xmax_recieve(1) - xmin_send(1))/dx(1) ) .GT. 0 .AND. &
    & NINT( (xmax_recieve(2) - xmin_send(2))/dx(2) ) .GT. 0 .AND. &
    & NINT( (xmax_send(1) - xmin_recieve(1))/dx(1) ) .GT. 0 .AND. &
    & NINT( (xmax_send(2) - xmin_recieve(2))/dx(2) ) .GT. 0 )THEN

!---------------------------------------------------------------------------
! Find indices of intersection
!---------------------------------------------------------------------------
    icomm = 1
    DO i = 1,pmlib_ndim

! min index
      IF( NINT( (xmin_recieve(i) - xmin_send(i))/dx(i) ) .EQ. 0 )THEN

        communication%imin_send(i,icomm)    = 1
        communication%imin_recieve(i,icomm) = 1

      ELSEIF( NINT( (xmin_send(i) - xmin_recieve(i))/dx(i) ) .GT. 0 )THEN

        communication%imin_send(i,icomm)    = 1
        communication%imin_recieve(i,icomm) = & 
              & NINT((xmin_send(i)-xmin_recieve(i))/dx(i)) + 1

      ELSEIF( NINT( (xmin_recieve(i) - xmin_send(i))/dx(i) ) .GT. 0 .AND. &
            & bc_recieve(2*i-1) .NE. 0 )THEN

        communication%imin_send(i,icomm)    = & 
              & NINT((xmin_recieve(i)-xmin_send(i))/dx(i)) + 1
        communication%imin_recieve(i,icomm) = 1

      ELSEIF( NINT( (xmin_recieve(i) - xmin_send(i))/dx(i) ) .GT. 0 .AND. &
            & bc_recieve(2*i-1) .EQ. 0 )THEN

        nghost = MIN( nghost_recieve(2*i-1), &
               & NINT((xmin_recieve(i)-xmin_send(i))/dx(i)) &
               & + nghost_send(2*i-1) )

        communication%imin_send(i,icomm)    = & 
              & NINT((xmin_recieve(i)-xmin_send(i))/dx(i)) + 1 - nghost
        communication%imin_recieve(i,icomm) = 1 - nghost

      ELSE
        ierr = -1
        CALL pmlib_write(caller,'Failed to identify intersection indexes')
        GOTO 9999
      ENDIF


! max index
      IF( NINT( (xmax_recieve(i) - xmax_send(i))/dx(i) ) .EQ. 0 )THEN

        communication%imax_send(i,icomm)    = ncell_send(i)
        communication%imax_recieve(i,icomm) = ncell_recieve(i)

      ELSEIF( NINT( (xmax_recieve(i) - xmax_send(i))/dx(i) ) .GT. 0 )THEN

        communication%imax_send(i,icomm)    = ncell_send(i)
        communication%imax_recieve(i,icomm) = & 
              & NINT((xmax_send(i)-xmin_recieve(i))/dx(i))

      ELSEIF( NINT( (xmax_send(i) - xmax_recieve(i))/dx(i) ) .GT. 0 .AND. &
            & bc_recieve(2*i) .NE. 0)THEN

        communication%imax_send(i,icomm) = & 
              & NINT((xmax_recieve(i)-xmin_send(i))/dx(i))
        communication%imax_recieve(i,icomm) = ncell_recieve(i)

      ELSEIF( NINT( (xmax_send(i) - xmax_recieve(i))/dx(i) ) .GT. 0 .AND. &
            & bc_recieve(2*i) .EQ. 0 )THEN

        nghost = MIN( nghost_recieve(2*i), &
               & NINT((xmax_recieve(i)-xmin_send(i))/dx(i)) &
               & + nghost_send(2*i) )

        communication%imax_send(i,icomm) = & 
              & NINT((xmax_recieve(i)-xmin_send(i))/dx(i)) + nghost
        communication%imax_recieve(i,icomm) = ncell_recieve(i) + nghost

      ELSE
        ierr = -1
        CALL pmlib_write(caller,'Failed to identify intersection indexes')
        GOTO 9999
      ENDIF

      communication%ncell_send(i,icomm) = &
                   & communication%imax_send(i,icomm) - &
                   &  communication%imin_send(i,icomm) + 1
      
      communication%ncell_recieve(i,icomm) = &
                   & communication%imax_recieve(i,icomm) - &
                   &  communication%imin_recieve(i,icomm) + 1

! check if the number of sended and recieved cells coinside
      IF(  communication%ncell_send(i,icomm) .NE. &
        &  communication%ncell_recieve(i,icomm))THEN
        ierr = -1
        CALL pmlib_write(caller, &
               & 'Number of sended and recieved cells does not correspond' )
        GOTO 9999
      END IF

      IF(communication%imin_send(i,icomm) .LT. 1-nghost_send(2*i-1))THEN
        CALL pmlib_write(caller,'Sending index is too low')
      END IF
      IF(communication%imax_send(i,icomm) .GT. ncell_send(i)+nghost_send(2*i))THEN
        CALL pmlib_write(caller,'Sending index is too high')
      END IF
      IF(communication%imin_recieve(i,icomm) .LT. 1-nghost_recieve(2*i-1))THEN
        CALL pmlib_write(caller,'Recieving index is too low')
      END IF
      IF(communication%imax_recieve(i,icomm) .GT. ncell_recieve(i)+nghost_recieve(2*i))THEN
        CALL pmlib_write(caller,'Recieving index is too high')
      END IF

    END DO ! direction
  END IF ! intersecting

!---------------------------------------------------------------------------
! Allocate communication buffer
!---------------------------------------------------------------------------
  IF( ASSOCIATED(comm_buffer) ) THEN
    DEALLOCATE(comm_buffer)
    IF (ierr .NE. 0) THEN
      CALL pmlib_write(caller,'Failed to deallocate array.')
      GOTO 9999
    ENDIF
  END IF 
  ALLOCATE( comm_buffer(ncomm) )

  DO j = 1,ncomm

    ncell_send(1) = communication%ncell_send(1,j)  
    ncell_send(2) = communication%ncell_send(2,j)
 
! Allocate send array for communication j
    ALLOCATE( comm_buffer(j)%mesh_send(nvar_send, &
            & ncell_send(1),ncell_send(2)) )
    comm_buffer(j)%nsend = nvar_send*ncell_send(1)*ncell_send(2)
    comm_buffer(j)%ncell_send = ncell_send

! Allocate send array for communication j
    ALLOCATE( comm_buffer(j)%mesh_recieve(nvar_send, &
            & ncell_send(1),ncell_send(2)) )
    comm_buffer(j)%nrecieve = nvar_send*ncell_send(1)*ncell_send(2)
    comm_buffer(j)%ncell_recieve = ncell_send
  END DO

!---------------------------------------------------------------------------
! set counter for number of packed variables
!---------------------------------------------------------------------------
  communication%npacked = 0

!---------------------------------------------------------------------------
! Return
!---------------------------------------------------------------------------
 9999 CONTINUE
RETURN
END SUBROUTINE pmlib_mesh_map

