! Copyright (c) 2015-2016, 2019, 2022 Harald Klimach <harald.klimach@dlr.de>
! Copyright (c) 2015 Samuel Ziegenberg
! Copyright (c) 2015 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2016 Neda Ebrahimi Pour <neda.epour@uni-siegen.de>
! Copyright (c) 2016 Tobias Girresser <tobias.girresser@student.uni-siegen.de>
! Copyright (c) 2016 Peter Vitt <peter.vitt2@uni-siegen.de>
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice, this
! list of conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice,
! this list of conditions and the following disclaimer in the documentation
! and/or other materials provided with the distribution.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
! DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
! FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!******************************************************************************!
!> Seeder Harvesting Tool
!! Visualization of meshes, generated by seeder.
!! (c) 2015 University of Siegen
!!
!! For a documentation, run ./waf gendoxy and find the documentation at
!! ./Documentation/html/index.html
program sdr_harvesting
  use iso_c_binding,            only: c_loc

  ! treelm modules
  use env_module,               only: pathLen

  use treelmesh_module,         only: treelmesh_type

  use tem_bc_prop_module,       only: tem_bc_prop_type
  use tem_color_prop_module,    only: tem_color_prop_type
  use tem_general_module,       only: tem_general_type, tem_start, tem_finalize
  use tem_logging_module,       only: logunit
  use tem_meshInfo_module,      only: tem_varSys_append_meshInfoVar
  use tem_stencil_module,       only: tem_stencilHeader_type, tem_create_stencil
  use tem_time_module,          only: tem_time_type, tem_time_reset
  use tem_tracking_module,      only: tem_tracking_getData,  &
    &                                 tem_init_tracker_subTree
  use tem_varsys_module,        only: tem_varsys_type, tem_varSys_init, &
    &                                 tem_varSys_proc_element
  use tem_varMap_module,        only: tem_varMap_type, tem_create_varMap
  use tem_restart_module,       only: tem_restart_type, tem_load_restart, &
    &                                 tem_init_restart
  use ply_dof_module,           only: q_space
  use ply_sampled_tracking_module, only: ply_sampled_track_init, &
    &                                    ply_sampled_track_output

  ! libharvesting
  use hvs_output_module,        only: hvs_output_file_type, hvs_output_open, &
    &                                 hvs_output_close, hvs_output_write,    &
    &                                 hvs_output_init, hvs_output_finalize

  ! sdr_harvesting
  use sdr_hvs_config_module,    only: sdr_hvs_config_type, sdr_hvs_config_load
  use sdr_hvs_props_module,     only: sdr_hvs_props_type,     &
    &                                 sdr_hvs_props_import_dofs

  ! seeder
  use sdr_protoData_module,     only: sdr_append_protoVar, access_state, &
                                      sdr_readRestart, sdr_protoData_type

  implicit none

  !-----------------------------------------------------------------------------
  type(tem_restart_type)           :: restart
  type(sdr_protoData_type), target :: protoData

  type(tem_general_type)       :: general
  type(treelmesh_type)         :: mesh
  type(sdr_hvs_config_type)    :: config
  type(sdr_hvs_props_type)     :: property
  type(hvs_output_file_type)   :: out_file
  type(tem_varsys_type)        :: varsys
  type(tem_varMap_type)        :: varMap
  type(tem_stencilHeader_type) :: stencil
  type(tem_time_type)          :: time
  procedure(tem_varSys_proc_element), pointer :: get_element
  integer              :: iTrack
  integer              :: nDofs
  integer              :: nSimpleVars
  integer, allocatable :: var_degree(:)
  integer, allocatable :: lvl_degree(:)
  integer, allocatable :: var_space(:)
  ! basename for output if tracking table is not defined
  ! trim(config%prefix)//trim(general%solver%simName)//trim(varSys%systemName)
  character(len=pathLen) :: basename
  !-----------------------------------------------------------------------------

  ! Initialize environment.
  call tem_start( codeName = 'sdr_harvesting', &
    &             general  = general           )


  ! Definition of the mesh variables.
  call tem_varSys_init( me         = varsys,        &
    &                   systemName = 'seeder_vars', &
    &                   length     = 8              )

  ! Is there a mesh, else use the restart files
  call tem_varSys_append_meshInfoVar(varsys)

  nSimpleVars = varsys%varname%nVals

  ! load config file
  call sdr_hvs_config_load( me         = config,   &
    &                       mesh       = mesh,     &
    &                       property   = property, &
    &                       varsys     = varsys,   &
    &                       general    = general,  &
    &                       restart    = restart,  &
    &                       time       = time      )


  if (restart%controller%readRestart) then
    get_element => access_state
    call sdr_append_protoVar(protoVar    = restart%header%varSys%varname   &
      &                                            %val(1:restart%header   &
      &                                            %varSys%varname%nVals), &
      &                      varSys      = varSys,                         &
      &                      method_data = c_loc(protoData),               &
      &                      get_element = get_element                     )

    ! create varMap with protoData variables
    ! VarMap is required to create MPI_type to read restart file
    call tem_create_varMap( varName = restart%header%varSys%varname  &
      &                                      %val(1:restart%header   &
      &                                      %varSys%varname%nVals), &
      &                     varSys  = varSys,                        &
      &                     varMap  = varMap                         )

    ! This routine creates mpi data_type to load data from restart file
    ! using MPI I/O. So, this must be initilized before readRestart
    call tem_init_restart( me       = restart,        &
      &                    solver   = general%solver, &
      &                    varMap   = varMap,         &
      &                    tree     = mesh            )

    ! read restart and fill protoData
    call sdr_readRestart(restart, mesh, protoData)
  end if

  allocate(var_degree(varsys%varname%nVals))
  allocate(lvl_degree(mesh%global%maxlevel))
  allocate(var_space(varsys%varname%nVals))
  var_degree = 0
  var_space = q_space

  if (config%do_subsampling) then
    call sdr_hvs_props_import_dofs( me        = property,                   &
      &                             mesh      = mesh,                       &
      &                             proc      = general%proc,               &
      &                             maxdegree = property%subres%polydegree, &
      &                             ndims     = 3                           )
    var_degree(nSimpleVars+1:varsys%varname%nVals) = property%subres%polydegree
  else
    ! No subsampling: load only integral mean values from the subresolution.
    write(logunit(2),*) 'Importing integral means from subresolution data'
    call sdr_hvs_props_import_dofs( me        = property,     &
      &                             mesh      = mesh,         &
      &                             proc      = general%proc, &
      &                             maxdegree = 0,            &
      &                             ndims     = 3             )
  end if
  lvl_degree = maxval(var_degree)

  if (config%ply_sample_track%tracking%control%active) then

    ! The stencil is needed for the identification of shapes, that are based
    ! on the boundary conditions.
    ! HK: For now we just use all 26 directions, that are available in the
    !     mesh right now. Probably, this should be configurable.
    call tem_create_stencil(stencil, 'd3q27')

    ! Mesh infos currently have just one degree of freedom per element.
    nDofs = 1

    ! No time associated with the mesh.
    call tem_time_reset(time)

    ! Initialize tracker subTree and remote empty trackers.
    ! It also adds solver%simName and varSys%systemName as prefix to
    ! tracking label
    call ply_sampled_track_init( me      = config%ply_sample_track, &
      &                          mesh    = mesh,                    &
      &                          solver  = general%solver,          &
      &                          varSys  = varSys,                  &
      &                          bc      = property%bc,             &
      &                          stencil = stencil,                 &
      &                          proc    = general%proc,            &
      &                          ndofs   = 1,                       &
      &                          ndims   = 3                        )

    call ply_sampled_track_output( me         = config%ply_sample_track, &
      &                            mesh       = mesh,                    &
      &                            bc         = property%bc,             &
      &                            solver     = general%solver,          &
      &                            proc       = general%proc,            &
      &                            varSys     = varSys,                  &
      &                            lvl_degree = lvl_degree,              &
      &                            var_degree = var_degree,              &
      &                            var_space  = var_space,               &
      &                            time       = time                     )

    do iTrack=1,config%ply_sample_track%tracking%control%nActive

      ! Finialize output
      call hvs_output_finalize( out_file = config%ply_sample_track          &
        &                                        %tracking%instance(iTrack) &
        &                                        %output_file               )

    end do

  else

    ! No tracking defined, just dump the plain original mesh.
    ! Open the output files, this also generates the vertices for the mesh,
    ! and writes the mesh data to disk.

    !!@todo HK: no subsampling if we dump all data (no tracking)?
    !!          Probably should do subsampling also here?

    ! filename for output
    basename = trim(config%prefix) // trim(general%solver%simName)

    ! Deactivate spatial reduction
    out_file%ascii%isReduce = .false.

    if (restart%controller%readRestart) then
      call hvs_output_init(out_file    = out_file,                   &
        &                  out_config  = config%output,              &
        &                  tree        = mesh,                       &
        &                  varSys      = varSys,                     &
        &                  varPos      = restart%varMap%varPos       &
        &                                %val(:varMap%varPos%nVals), &
        &                  basename    = trim(basename),             &
        &                  globProc    = general%proc,               &
        &                  solver      = general%solver              )
    else
      call hvs_output_init(out_file    = out_file,       &
        &                  out_config  = config%output,  &
        &                  tree        = mesh,           &
        &                  varSys      = varSys,         &
        &                  basename    = trim(basename), &
        &                  globProc    = general%proc,   &
        &                  solver      = general%solver  )
    end if

    ! Open output file handle
    call hvs_output_open( out_file   = out_file, &
      &                   use_iter   = .false.,  &
      &                   mesh       = mesh,     &
      &                   varsys     = varsys    )

    ! Fill output files with data.
    call hvs_output_write( out_file = out_file, &
      &                    varsys   = varsys,   &
      &                    mesh     = mesh      )

    ! Close output files again.
    call hvs_output_close( out_file = out_file, &
      &                    varsys   = varsys,   &
      &                    mesh     = mesh      )

    ! Finialize output
    call hvs_output_finalize( out_file = out_file )
  end if

  ! Finalize environment.
  call tem_finalize(general)

end program sdr_harvesting
