!    Copyright (C) 2006 Imperial College London and others.
!    
!    Please see the AUTHORS file in the main source directory for a full list
!    of copyright holders.
!
!    Prof. C Pain
!    Applied Modelling and Computation Group
!    Department of Earth Science and Engineeringp
!    Imperial College London
!
!    amcgsoftware@imperial.ac.uk
!    
!    This library is free software; you can redistribute it and/or
!    modify it under the terms of the GNU Lesser General Public
!    License as published by the Free Software Foundation,
!    version 2.1 of the License.
!
!    This library 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
!    Lesser General Public License for more details.
!
!    You should have received a copy of the GNU Lesser General Public
!    License along with this library; if not, write to the Free Software
!    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
!    USA

#include "fdebug.h"

subroutine project_to_continuous_parallel(vtuname, vtuname_len, meshname,&
     & meshname_len) 
  !!< Given a vtu file containing fields on a discontinuous mesh, and the
  !!< mesh files for the corresponding continuous mesh, produce a vtu
  !!< with its fields projected onto the continuous mesh.
  
  use fefields
  use fields
  use halos
  use global_parameters, only : current_debug_level
  use parallel_tools
  use mesh_files
  use solvers
  use sparsity_patterns
  use sparse_tools
  use state_module
  use vtk_interfaces
  
  implicit none
  
  integer, intent(in):: vtuname_len, meshname_len
  character(len=vtuname_len), intent(in):: vtuname
  character(len=meshname_len), intent(in):: meshname

  character(len = *), parameter :: solver_option_path = "/dummy"
  integer :: i, nfields
  ! VTK can store at most P2 fields. Hence degree 4 allows complete quadrature
  ! in the Galerkin projection.
  integer, parameter :: quad_degree = 4
  type(csr_matrix) :: mass
  type(csr_sparsity) :: sparsity
  type(mesh_type), pointer :: cg_mesh
  type(scalar_field) :: cg_s_field
  type(scalar_field), dimension(:), allocatable :: rhs
  type(scalar_field), dimension(:), pointer :: cg_state_fields, dg_state_fields
  type(scalar_field), pointer :: dg_s_field
  type(state_type) :: cg_state, dg_state
  type(tensor_field) :: cg_t_field
  type(tensor_field), pointer :: dg_t_field
  type(vector_field) :: cg_v_field
  type(vector_field), pointer :: dg_v_field
  type(vector_field), target :: cg_coordinate
  
  !current_debug_level = 2

  if(isparallel()) then
    call vtk_read_state(parallel_filename(vtuname, ".vtu"), dg_state, quad_degree = quad_degree)
  else
    call vtk_read_state(trim(vtuname) //".vtu", dg_state, quad_degree = quad_degree)
  end if
  
  cg_coordinate = read_mesh_files(meshname, quad_degree = quad_degree)
  if(isparallel()) call read_halos(meshname, cg_coordinate)
  cg_mesh => cg_coordinate%mesh

  do i = 1, scalar_field_count(dg_state)
    dg_s_field => extract_scalar_field(dg_state, i)
    call allocate(cg_s_field, cg_mesh, name = dg_s_field%name)
    call insert(cg_state, cg_s_field, cg_s_field%name)
    call deallocate(cg_s_field)
  end do
  
  do i = 1, vector_field_count(dg_state)
    dg_v_field => extract_vector_field(dg_state, i)
    call allocate(cg_v_field, dg_v_field%dim, cg_mesh, name = dg_v_field%name)
    call insert(cg_state, cg_v_field, cg_v_field%name)
    call deallocate(cg_v_field)
  end do
  
  do i = 1, tensor_field_count(dg_state)
    dg_t_field => extract_tensor_field(dg_state, i)
    call allocate(cg_t_field, cg_mesh, name = dg_t_field%name)
    call insert(cg_state, cg_t_field, cg_t_field%name)
    call deallocate(cg_t_field)
  end do

  sparsity = make_sparsity(cg_mesh, cg_mesh, name = "MassSparsity")
  call allocate(mass, sparsity, name = "MassMatrix")
  call deallocate(sparsity)
  call compute_mass(cg_coordinate, cg_mesh, mass)

  call collapse_state(dg_state, dg_state_fields)
  nfields = size(dg_state_fields)
  
  allocate(rhs(nfields))
  do i = 1, nfields
    call allocate(rhs(i), cg_mesh, name = trim(dg_state_fields(i)%name) // "Rhs")
    call zero(rhs(i))
  end do

  do i = 1, ele_count(cg_mesh)
    call assemble_gp_rhs_ele(i, cg_mesh, cg_coordinate, dg_state_fields, rhs)
  end do
  deallocate(dg_state_fields)
  call deallocate(dg_state)
  
  call collapse_state(cg_state, cg_state_fields)
  call set_solver_options(solver_option_path, &
    & ksptype = "cg", pctype = "sor", atol = epsilon(0.0), rtol = 0.0, max_its = 2000, start_from_zero = .true.)
  call petsc_solve(cg_state_fields, mass, rhs, option_path = solver_option_path)
  
  deallocate(cg_state_fields)
  call deallocate(mass)
  do i = 1, nfields
    call deallocate(rhs(i))
  end do
  deallocate(rhs)
  
  call insert(cg_state, cg_coordinate, "Coordinate")
  call insert(cg_state, cg_mesh, "CoordinateMesh")
  call deallocate(cg_coordinate)
  
  if(isparallel()) then
    call vtk_write_state(trim(vtuname) //"_continuous.pvtu", state = (/cg_state/))
  else
    call vtk_write_state(trim(vtuname) //"_continuous.vtu", state = (/cg_state/))
  end if
  call deallocate(cg_state)
  
  call print_references(0)
       
contains

  subroutine assemble_gp_rhs_ele(ele, mesh, positions, fields, rhs)
    integer, intent(in) :: ele
    type(mesh_type), intent(in) :: mesh
    type(vector_field), intent(in) :: positions
    type(scalar_field), dimension(:), intent(in) :: fields
    type(scalar_field), dimension(size(fields)), intent(inout) :: rhs
       
    integer :: i
    integer, dimension(:), pointer :: element_nodes
    real, dimension(ele_ngi(mesh, ele)) :: detwei
    type(element_type), pointer :: shape
    
    call transform_to_physical(positions, ele, detwei = detwei)
    
    shape => ele_shape(mesh, ele)
    element_nodes => ele_nodes(mesh, ele)

    do i = 1, size(fields)
      call addto(rhs(i), element_nodes, shape_rhs(shape, detwei * ele_val_at_quad(fields(i), ele)))
    end do
    
  end subroutine assemble_gp_rhs_ele
       
end subroutine project_to_continuous_parallel
