! File: amr_driver.f90        
!
!
!
       
	Module AMRDVR_Mod
 	Use GlobalVars_mod
        Use Mesh_mod
	Use AMR_mod
	Use Util_mod
	Use AmrComm_Mod 
        USE mpi_module, only: ppinit
        implicit none
	save

        Type(Mesh) parent_mesh

	Contains

! ----------------------------------------------------------------------
! Initialize the Parallel AMR System and check status of MPI
! ----------------------------------------------------------------------
         subroutine PAMR_INIT()
           logical :: mpi_active

           call MPI_INITIALIZED( mpi_active, ierr )

           if ( .not. mpi_active ) then
                call MPI_INIT(ierr)
                if ( ierr /= MPI_SUCCESS ) then
                    call Abort_Program("PAMR_Init() failed")
                end if
           end if

           call MPI_COMM_RANK(MPI_COMM_WORLD, myPid, ierr)
           call MPI_COMM_SIZE(MPI_COMM_WORLD, totalProcs, ierr)
           call ppinit()

         end subroutine PAMR_INIT

! ----------------------------------------------------------------------
! Finalize the Parallel AMR System and check status of MPI
! ----------------------------------------------------------------------
         subroutine PAMR_FINALIZE( mpi_active )
           logical, intent(in), optional :: mpi_active

           call MPI_BARRIER( MPI_COMM_WORLD, ierr )

           if ( mpi_active ) then
                call MPI_FINALIZE( ierr )
           end if

         end subroutine PAMR_FINALIZE

        subroutine PAMR_LOGICAL_AMR(this_mesh)
	Type(Mesh) :: this_mesh
	!
	! Parallel mesh adjustment
        Call Mesh_Adjustment(this_mesh, parent_mesh)
        end subroutine PAMR_LOGICAL_AMR


        subroutine PAMR_PHYSICAL_AMR(in_mesh)
	!
	! input mesh and the mesh to be returned  
        type (Mesh) :: in_mesh
        type (Mesh), Save :: out_mesh

	! Parallel physical mesh refinement
        call Mesh_Refinement(in_mesh, out_mesh)

	! Resolve global ids
        call GlobalId_Resolve(out_mesh, in_mesh)

        call Delete_mesh(parent_mesh)
	
        ! Transfer in_mesh to parent_mesh
        ! After this call, in_mesh will be undefined
        call Mesh_transfer(parent_mesh, in_mesh)
        parent_mesh%defined = .true.; in_mesh%defined = .false.
 
        ! Transfer out_mesh to in_mesh
        ! After this call, out_mesh will be undefined
        call Mesh_transfer(in_mesh, out_mesh)
        in_mesh%defined = .true.; out_mesh%defined = .false.
	! 
        END SUBROUTINE PAMR_PHYSICAL_AMR


	Subroutine First_mesh_partition_migrate(this_mesh)
        Use meshmove_module
        Type(Mesh) :: this_mesh
	End Subroutine First_mesh_partition_migrate


	Subroutine PAMR_REPARTITION(this_mesh)
	Use GlobalVars_mod
        Use Mesh_mod
        Use meshmove_module
	Type(Mesh) :: this_mesh
        logical, save :: first_entry = .true.

          if ( first_entry ) then
               nullify( this_mesh%LocErr_List ) ; first_entry = .false.
          end if

          if (associated(this_mesh%LocErr_List)) then
              deallocate(this_mesh%LocErr_List)
          end if

          allocate(this_mesh%LocErr_List(this_mesh%loc_elements),       &
     &             stat=err_code)

          if (err_code /= 0) then
              call Abort_Program("PAMR_REPARTITION: allocate() failed")
          end if

	End Subroutine PAMR_REPARTITION


	Subroutine Compute_element_weight(this_mesh)
	Type(Mesh) :: this_mesh
        Type(Element), dimension(:), pointer :: elmt_array
        integer :: idx

        this_mesh%elements%weight = 1

	elmt_array => this_mesh%elements	  
	!
        Do idx = 1, this_mesh%loc_elements
	  if(elmt_array(idx)%refined_type /= 2) then 
	    elmt_array(idx)%weight = elmt_array(idx)%refine_ratio
	  else if(elmt_array(idx)%refined_type == 2) then
	    if(elmt_array(idx)%refine_ratio == 3) then
	      elmt_array(idx)%weight = 2
	    else if(elmt_array(idx)%refine_ratio == 4) then
	      elmt_array(idx)%weight = 3
	    endif
	  endif
	EndDo

	Do idx = 1, this_mesh%loc_elements
	  if(elmt_array(idx)%weight < 0 .or. elmt_array(idx)%weight > 4) then
	    call Abort_program("Illegal element weight!")
	  endif	
	EndDo
	End Subroutine Compute_element_weight


	! Currently, this routine is basically the same as Build_mesh_struct()
	! except that argument arrays are regular arrays instead of pointers
	Subroutine Import_mesh_struct(in_mesh, loc_elmts, &
                         elmt_node_array, elmt_edge_array, &
                         node_coord_array, elmt_neigh, &
                         my_pid, tot_procs)
        Use GlobalVars_Mod      
        Use Mesh_Mod   
        implicit none  
	! 
        Type(Mesh) in_mesh     
        integer, intent (in) :: loc_elmts, my_pid, tot_procs   
        integer, dimension(:) ::  elmt_node_array, &   
                               elmt_edge_array, elmt_neigh     
        real, dimension(:) :: node_coord_array
        integer, dimension(:), pointer :: tmp_nodes, tmp_edges 
        integer, dimension(:), pointer :: sum_array    
        integer :: idx, jdx, offset    
        logical :: found       

	in_mesh%loc_elements = loc_elmts
        allocate(tmp_nodes(ele_nodes*in_mesh%loc_elements), stat=err_code)
        allocate(tmp_edges(ele_edges*in_mesh%loc_elements), stat=err_code)
        allocate(sum_array(tot_procs), stat=err_code)

        ! Compute offset for global element id
        sum_array = 0
        if(myPid < tot_procs-1) then
          sum_array(my_pid+1:tot_procs) = in_mesh%loc_elements
        endif
        call MPI_ALLREDUCE(sum_array, sum_array, tot_procs, MPI_INTEGER, &
                           MPI_SUM, MPI_COMM_WORLD, err_code)
 
        in_mesh%loc_nodes = 0; in_mesh%loc_edges = 0
        !
        ! Build node arrays
        !
        Do idx = 1, ele_nodes*in_mesh%loc_elements
          found = .false.
          Do jdx = 1, idx-1
            if(elmt_node_array(idx) == elmt_node_array(jdx)) then
               tmp_nodes(idx) = tmp_nodes(jdx) ! local index
              found = .true.; exit
            endif
          EndDo
          if(.not. found) then
            in_mesh%loc_nodes = in_mesh%loc_nodes + 1
            tmp_nodes(idx) = in_mesh%loc_nodes ! local indx
          endif 
        EndDo

	!
        ! Build edge arrays
        !          
        Do idx = 1, ele_edges*in_mesh%loc_elements
          found = .false.
          Do jdx = 1, idx-1
            if(elmt_edge_array(idx) == elmt_edge_array(jdx)) then
              tmp_edges(idx) = tmp_edges(jdx) ! local index
              found = .true.; exit
            endif
          EndDo
          if(.not. found) then
            in_mesh%loc_edges = in_mesh%loc_edges + 1
            tmp_edges(idx) = in_mesh%loc_edges ! local indx
          endif
        EndDo
 
        !
        ! Create local mesh structure
        !
        call Mesh_init(in_mesh, in_mesh%loc_elements, &
                       in_mesh%loc_edges, in_mesh%loc_nodes)
 
        ! Assign local ids in edge and element structures
        Do idx = 1, in_mesh%loc_elements
          in_mesh%elements(idx)%node_indx(1:3) = tmp_nodes(3*idx-2:3*idx)
          in_mesh%elements(idx)%edge_indx(1:3) = tmp_edges(3*idx-2:3*idx)
        EndDo
        Do idx = 1, in_mesh%loc_elements
          in_mesh%edges(in_mesh%elements(idx)%edge_indx(1))%node_indx(1) &
            = in_mesh%elements(idx)%node_indx(1)
          in_mesh%edges(in_mesh%elements(idx)%edge_indx(1))%node_indx(2) &
            = in_mesh%elements(idx)%node_indx(2)
          in_mesh%edges(in_mesh%elements(idx)%edge_indx(2))%node_indx(1) &
            = in_mesh%elements(idx)%node_indx(2)
          in_mesh%edges(in_mesh%elements(idx)%edge_indx(2))%node_indx(2) &
            = in_mesh%elements(idx)%node_indx(3)
          in_mesh%edges(in_mesh%elements(idx)%edge_indx(3))%node_indx(1) &
            = in_mesh%elements(idx)%node_indx(3)
          in_mesh%edges(in_mesh%elements(idx)%edge_indx(3))%node_indx(2) &
            = in_mesh%elements(idx)%node_indx(1)
        EndDo

	! Assign global ids for nodes and edges
        Do idx = 1, in_mesh%loc_elements
          in_mesh%nodes(in_mesh%elements(idx)%node_indx(1))%global_id = &
                   elmt_node_array(3*idx-2)
          in_mesh%nodes(in_mesh%elements(idx)%node_indx(2))%global_id = &
                   elmt_node_array(3*idx-1)
          in_mesh%nodes(in_mesh%elements(idx)%node_indx(3))%global_id = &
                   elmt_node_array(3*idx)
          in_mesh%edges(in_mesh%elements(idx)%edge_indx(1))%global_id = &
                   elmt_edge_array(3*idx-2)
          in_mesh%edges(in_mesh%elements(idx)%edge_indx(2))%global_id = &
                   elmt_edge_array(3*idx-1)
          in_mesh%edges(in_mesh%elements(idx)%edge_indx(3))%global_id = &
                   elmt_edge_array(3*idx)
        EndDo

	! Assign global id for elements
        if(my_pid == 0) then
          offset = 0
        else
          offset = sum_array(my_pid)
        endif
        ! 
        Do idx = 1, in_mesh%loc_elements
          in_mesh%elements(idx)%global_id = offset + 1
          offset = offset + 1
        EndDo
          
        ! Build node coord
        Do idx = 1, in_mesh%loc_nodes
          offset = in_mesh%nodes(idx)%global_id
          in_mesh%nodes(idx)%coord(1:2) = node_coord_array(2*offset-1:2*offset)
        EndDo
          
        Do idx = 1, in_mesh%loc_elements
          in_mesh%elements(idx)%neigh_indx(1:3) = elmt_neigh(3*idx-2:3*idx)
        EndDo
        !
        End Subroutine Import_mesh_struct


	Subroutine Export_mesh_struct(trimesh, elmt_nodes, &
	        	  elmt_edges, node_coords, elmt_neigh)
	Type(Mesh) trimesh
	Type(Element), dimension(:), pointer :: elements
	Type(Node), dimension(:), pointer :: nodes
	integer, dimension(:) :: elmt_nodes, elmt_edges
        real, dimension(:) :: node_coords
        integer, dimension(:) :: elmt_neigh
 	integer :: idx, jdx

	elements => trimesh%elements
	jdx = 1 
	Do idx = 1, trimesh%loc_elements
	  elmt_nodes(jdx:jdx+ele_nodes-1) = &
	    elements(idx)%node_indx(1:ele_nodes) 
	  elmt_edges(jdx:jdx+ele_edges-1) = &
            elements(idx)%edge_indx(1:ele_edges)
	  !
          if(associated(trimesh%boundary_elements)) then
            elmt_neigh(jdx) = &
	     trimesh%boundary_elements(elements(idx)%neigh_indx(1))%global_id
	    elmt_neigh(jdx+1) = & 
             trimesh%boundary_elements(elements(idx)%neigh_indx(2))%global_id
	    elmt_neigh(jdx+2) = & 
             trimesh%boundary_elements(elements(idx)%neigh_indx(3))%global_id
	  endif
	  jdx = jdx + ele_nodes
	EndDo
	 
	nodes => trimesh%nodes
        jdx = 1
	Do idx = 1, trimesh%loc_nodes
	  node_coords(jdx:jdx+mesh_dim-1) = nodes(idx)%coord(1:mesh_dim)
	  jdx = jdx + mesh_dim
	EndDo
	!
	End Subroutine Export_mesh_struct


        ! Old AMR driver interface
        !
        Subroutine PAMR_Driver(in_mesh)
        implicit none
        !
        ! input-mesh info
        Type(Mesh) :: in_mesh
        !
        ! Parent and output mesh structures
        Type(Mesh), Save :: out_mesh

        ! Set neighboring elments info
        !Call Set_ElmtNeighbor_Info(in_mesh)

        ! Parallel mesh adjustment
        Call Mesh_Adjustment(in_mesh, parent_mesh)

        ! Parallel mesh refinement
        Call Mesh_Refinement(in_mesh, out_mesh)

        ! Transfer in_mesh to parent_mesh
        ! After this call, in_mesh will be undefined
        Call Mesh_transfer(parent_mesh, in_mesh)
        parent_mesh%defined = .true.; in_mesh%defined = .false.
 
        ! Transfer out_mesh to in_mesh
        ! After this call, out_mesh will be undefined
        Call Mesh_transfer(in_mesh, out_mesh)
        in_mesh%defined = .true.; out_mesh%defined = .false.
        !
        return
        End Subroutine PAMR_driver
	!
	End Module AMRDVR_Mod

        MODULE pyramid_module
          USE GlobalVars_mod
          USE Mesh_mod
          USE MPI_mod
          USE AMRDVR_Mod
          USE Util_mod
          USE mpi_module
          USE meshmove_module
        END MODULE pyramid_module
