~reducedmodelling/fluidity/ReducedModel

!    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 differentiate_vtu(input_filename, input_filename_len, output_filename, output_filename_len, input_fieldname, input_fieldname_len)

! these 5 need to be on top and in this order, so as not to confuse silly old intel compiler 

  use quadrature

  use elements

  use sparse_tools

  use fields

  use state_module

!

  use field_derivatives

  use fldebug

  use state_module

  use vtk_interfaces

  implicit none

  integer, intent(in) :: input_filename_len

  integer, intent(in) :: output_filename_len

  integer, intent(in) :: input_fieldname_len

  character(len = input_filename_len), intent(in) :: input_filename

  character(len = output_filename_len), intent(in) :: output_filename

  character(len = input_fieldname_len), intent(in) :: input_fieldname

  integer :: dim, i, j, nfields

  logical :: allocated

  type(mesh_type), pointer :: mesh

  type(scalar_field) :: masslump

  type(scalar_field), dimension(:), allocatable :: s_fields

  type(scalar_field), pointer :: s_field

  type(state_type) :: collapsed_state, state

  type(vector_field) :: field_grad

  type(vector_field), dimension(:), allocatable :: field_grads

  type(vector_field), pointer :: positions

  ewrite(1, *) "In differentiate_vtu"

  call vtk_read_state(trim(input_filename), state)

  positions => extract_vector_field(state, "Coordinate")

  dim = positions%dim

  mesh => extract_mesh(state, "Mesh")

  if(len_trim(input_fieldname) == 0) then

    call collapse_fields_in_state(state, collapsed_state)

    nfields = scalar_field_count(collapsed_state)

    allocate(s_fields(nfields))

    do i = 1, scalar_field_count(collapsed_state)

      s_fields(i) = extract_scalar_field(collapsed_state, i)

    end do

    allocate(field_grads(nfields))

    do i = 1, nfields

      s_field => extract_scalar_field(collapsed_state, i)

      call allocate(field_grads(i), dim, mesh, trim(s_field%name) // "Gradient")

    end do

    call deallocate(collapsed_state)

    select case(continuity(positions))

      case(-1)

        do i = 1, ele_count(mesh)

          call solve_grad_ele(i, positions, s_fields, field_grads)

        end do

      case(0)

        call allocate(masslump, mesh, "LumpedMass")

        call zero(masslump)

        do i = 1, nfields

          call zero(field_grads(i))

        end do

        do i = 1, ele_count(mesh)

          call assemble_grad_ele(i, positions, s_fields, masslump, field_grads)

        end do

        do i = 1, nfields

          do j = 1, dim

            field_grads(i)%val(j,:) = field_grads(i)%val(j,:) / masslump%val

          end do

        end do

        call deallocate(masslump)

      case default

        ewrite(-1, *) "For continuity ", continuity(positions)

        FLAbort("Unrecognised continuity")

    end select

    do i = 1, nfields

      call insert(state, field_grads(i), field_grads(i)%name)

      call deallocate(field_grads(i))

    end do

    deallocate(field_grads)

    deallocate(s_fields)

  else

    s_field => extract_scalar_field(state, trim(input_fieldname), allocated = allocated)

    call allocate(field_grad, dim, mesh, trim(s_field%name) // "Gradient")

    call grad(s_field, positions, field_grad)

    call insert(state, field_grad, field_grad%name) 

    call deallocate(field_grad)   

    if(allocated) deallocate(s_field)

  end if

  call vtk_write_state(output_filename, state = (/state/))

  call deallocate(state)

  call print_references(0)

  ewrite(1, *) "Exiting differentate_vtu"

contains

  subroutine solve_grad_ele(ele, positions, fields, field_grads)

    integer, intent(in) :: ele

    type(vector_field), intent(in) :: positions

    type(scalar_field), dimension(:), intent(in) :: fields

    type(vector_field), dimension(size(fields)), intent(inout) :: field_grads

    integer :: i

    integer, dimension(:), pointer :: nodes

    real, dimension(ele_ngi(positions, ele)) :: detwei

    real, dimension(ele_loc(positions, ele), ele_loc(positions, ele)) :: little_mass

    real, dimension(ele_loc(positions, ele), size(fields) * positions%dim) :: little_rhs

    real, dimension(ele_loc(positions, ele), ele_ngi(positions, ele), positions%dim) :: dshape

    type(element_type), pointer :: shape

    shape => ele_shape(positions, ele)

    call transform_to_physical(positions, ele, shape, &

      & detwei = detwei, dshape = dshape)

    little_mass = shape_shape(shape, shape, detwei)

    do i = 1, size(fields)

      little_rhs(:, (i - 1) * positions%dim + 1:i * positions%dim) = transpose(shape_vector_rhs(shape, ele_grad_at_quad(fields(i), ele, dshape), detwei))

    end do

    call solve(little_mass, little_rhs)

    nodes => ele_nodes(positions, ele)

    do i = 1, size(field_grads)

      call set(field_grads(i), nodes, little_rhs(:, (i - 1) * positions%dim + 1:i * positions%dim))

    end do

  end subroutine solve_grad_ele

  subroutine assemble_grad_ele(ele, positions, fields, masslump, rhs)

    integer, intent(in) :: ele

    type(vector_field), intent(in) :: positions

    type(scalar_field), dimension(:), intent(in) :: fields

    type(scalar_field), intent(inout) :: masslump

    type(vector_field), dimension(size(fields)), intent(inout) :: rhs

    integer :: i

    integer, dimension(:), pointer :: nodes

    real, dimension(ele_ngi(positions, ele)) :: detwei

    real, dimension(ele_loc(positions, ele), ele_ngi(positions, ele), positions%dim) :: dshape

    type(element_type), pointer :: shape

    shape => ele_shape(positions, ele)

    call transform_to_physical(positions, ele, shape, &

      & detwei = detwei, dshape = dshape)

    nodes => ele_nodes(positions, ele)

    call addto(masslump, nodes, sum(shape_shape(shape, shape, detwei), 2))

    do i = 1, size(rhs)

      call addto(rhs(i), nodes, shape_vector_rhs(shape, ele_grad_at_quad(fields(i), ele, dshape), detwei))

    end do

  end subroutine assemble_grad_ele

end subroutine differentiate_vtu