~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 supermesh_difference(vtu1_filename, vtu1_filename_len, vtu2_filename, vtu2_filename_len, output_filename, output_filename_len)

  use fields

  use fldebug

  use interpolation_module

  use intersection_finder_module

  use linked_lists

  use reference_counting

  use state_module

  use supermesh_construction

  use unify_meshes_module

  use vtk_interfaces

  implicit none

  integer, intent(in) :: vtu1_filename_len

  integer, intent(in) :: vtu2_filename_len

  integer, intent(in) :: output_filename_len

  character(len = vtu1_filename_len), intent(in) :: vtu1_filename

  character(len = vtu2_filename_len), intent(in) :: vtu2_filename

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

  integer :: dim, ele_1, ele_2, ele_3, i, j, nintersections, nintersections_max

  integer, dimension(:), allocatable :: ele_map_13, ele_map_23, node_map_13, &

    & node_map_23

  integer, dimension(:, :), allocatable :: intersection_parents

  logical :: p0

  type(element_type), pointer :: shape

  type(ilist), dimension(:), allocatable :: ele_map_21

  type(inode), pointer :: node

  type(mesh_type) :: pwc_mesh_3

  type(mesh_type), pointer :: mesh_3

  type(scalar_field) :: s_field_3

  type(scalar_field), pointer :: s_field, s_field_13, s_field_23

  type(state_type) :: state_1, state_2, state_3

  type(state_type), dimension(2) :: state_13, state_1_split, state_23, &

    & state_2_split, state_3_split

  type(tensor_field) :: t_field_3

  type(tensor_field), pointer :: t_field, t_field_13, t_field_23

  type(vector_field), pointer :: positions_1, positions_2, v_field, &

    & v_field_13, v_field_23

  type(vector_field) :: intersection, v_field_3

  type(vector_field), target :: positions_3

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

  ewrite(1, *) "In supermesh_difference"

  call vtk_read_state(vtu1_filename, state_1)

  call vtk_read_state(vtu2_filename, state_2)

  p0 = has_mesh(state_1, "P0Mesh")

  positions_1 => extract_vector_field(state_1, "Coordinate")

  positions_2 => extract_vector_field(state_2, "Coordinate")

  dim = positions_1%dim

  if(positions_2%dim /= dim) then

    FLExit("Input vtu dimensions do not match")

  end if

  call intersector_set_dimension(dim)

  call intersector_set_exactness(.false.)

  allocate(ele_map_21(ele_count(positions_1)))

  ! Use the rtree to avoid continuity assumptions

  ele_map_21 = rtree_intersection_finder(positions_1, positions_2)

  nintersections_max = sum(ele_map_21%length)

  ewrite(2, "(a,i0)") "Maximum number of intersections: ", nintersections_max

  allocate(intersections(nintersections_max))

  allocate(intersection_parents(nintersections_max, 2))

  nintersections = 0

  do ele_1 = 1, ele_count(positions_1)

    node => ele_map_21(ele_1)%firstnode

    do while(associated(node))

      ele_2 = node%value

      ! TODO: Integrate the tet intersector

      intersection = intersect_elements(positions_1, ele_1, ele_val(positions_2, ele_2), ele_shape(positions_1, ele_1))

      if(ele_count(intersection) > 0) then

        nintersections = nintersections + 1

        intersections(nintersections) = intersection

        intersection_parents(nintersections, :) = (/ele_1, ele_2/)

      else

        call deallocate(intersection)

      end if

      node => node%next

    end do

  end do  

  call deallocate(ele_map_21)

  deallocate(ele_map_21)

  ewrite(2, "(a,i0)") "Number of intersections: ", nintersections

  positions_3 = unify_meshes(intersections(:nintersections))

  positions_3%name = "Coordinate"

  allocate(node_map_13(node_count(positions_3)))

  allocate(node_map_23(node_count(positions_3)))

  if(p0) then

    allocate(ele_map_13(ele_count(positions_3)))

    allocate(ele_map_23(ele_count(positions_3)))

  end if

  ele_3 = 0

  do i = 1, nintersections

    do j = 1, ele_count(intersections(i))

      ele_3 = ele_3 + 1

      node_map_13(ele_nodes(positions_3, ele_3)) = intersection_parents(i, 1)

      node_map_23(ele_nodes(positions_3, ele_3)) = intersection_parents(i, 2)

      if(p0) then

        ele_map_13(ele_3) = intersection_parents(i, 1)

        ele_map_23(ele_3) = intersection_parents(i, 2)

      end if

    end do

    call deallocate(intersections(i))

  end do

  deallocate(intersections)

  deallocate(intersection_parents)

  mesh_3 => positions_3%mesh

  pwc_mesh_3 = piecewise_constant_mesh(mesh_3, "PiecewiseConstantMesh")

  call insert(state_3, positions_3, "Coordinate")

  call insert(state_1_split(1), positions_1, "Coordinate")

  call insert(state_2_split(1), positions_2, "Coordinate")

  call insert(state_13(1), positions_3, "Coordinate")

  call insert(state_23(1), positions_3, "Coordinate")

  call insert(state_3_split(1), positions_3, "Coordinate")

  call insert(state_1_split(1), positions_1%mesh, "Mesh")

  call insert(state_13(1), mesh_3, "Mesh")

  call insert(state_2_split(1), positions_2%mesh, "Mesh")

  call insert(state_23(1), mesh_3, "Mesh")

  call insert(state_3, mesh_3, "Mesh")

  call insert(state_3_split(1), mesh_3, "Mesh")

  if(p0) then

    call insert(state_1_split(2), positions_1, "Coordinate")

    call insert(state_2_split(2), positions_2, "Coordinate")

    call insert(state_13(2), positions_3, "Coordinate")

    call insert(state_23(2), positions_3, "Coordinate")

    call insert(state_3_split(2), positions_3, "Coordinate")

    call insert(state_1_split(2), extract_mesh(state_1, "P0Mesh"), "Mesh")

    call insert(state_2_split(2), extract_mesh(state_2, "P0Mesh"), "Mesh")

    call insert(state_13(2), pwc_mesh_3, "Mesh")

    call insert(state_23(2), pwc_mesh_3, "Mesh")

    call insert(state_3_split(2), pwc_mesh_3, "Mesh")

  end if

  call deallocate(positions_3)

  do i = 1, scalar_field_count(state_1)

    s_field => extract_scalar_field(state_1, i)

    shape => ele_shape(s_field, 1)

    if(shape%degree == 0) then    

      assert(p0)      

      call insert(state_1_split(2), s_field, s_field%name)

      call insert(state_2_split(2), extract_scalar_field(state_2, s_field%name), s_field%name)

      call allocate(s_field_3, pwc_mesh_3, s_field%name)

      call insert(state_13(2), s_field_3, s_field_3%name)

    else

      call insert(state_1_split(1), s_field, s_field%name)

      call insert(state_2_split(1), extract_scalar_field(state_2, s_field%name), s_field%name)

      call allocate(s_field_3, mesh_3, s_field%name)

      call insert(state_13(1), s_field_3, s_field_3%name)

    end if

    call deallocate(s_field_3)

    if(shape%degree == 0) then

      call allocate(s_field_3, pwc_mesh_3, s_field%name)

      call insert(state_23(2), s_field_3, s_field_3%name)

    else

      call allocate(s_field_3, mesh_3, s_field%name)

      call insert(state_23(1), s_field_3, s_field_3%name)

    end if

    call deallocate(s_field_3)

    if(shape%degree == 0) then

      call allocate(s_field_3, pwc_mesh_3, s_field%name)

      call insert(state_3_split(2), s_field_3, s_field_3%name)

    else

      call allocate(s_field_3, mesh_3, s_field%name)

      call insert(state_3_split(1), s_field_3, s_field_3%name)

    end if

    call insert(state_3, s_field_3, s_field_3%name)

    call deallocate(s_field_3)    

  end do    

  do i = 1, vector_field_count(state_1)

    v_field => extract_vector_field(state_1, i)

    if(v_field%name == "Coordinate") cycle

    shape => ele_shape(v_field, 1)

    if(shape%degree == 0) then    

      assert(p0)      

      call insert(state_1_split(2), v_field, v_field%name)

      call insert(state_2_split(2), extract_scalar_field(state_2, v_field%name), v_field%name)

      call allocate(v_field_3, dim, pwc_mesh_3, v_field%name)

      call insert(state_13(2), v_field_3, v_field_3%name)

    else

      call insert(state_1_split(1), v_field, v_field%name)

      call insert(state_2_split(1), extract_scalar_field(state_2, v_field%name), v_field%name)

      call allocate(v_field_3, dim, mesh_3, v_field%name)

      call insert(state_13(1), v_field_3, v_field_3%name)

    end if

    call deallocate(v_field_3)

    if(shape%degree == 0) then

      call allocate(v_field_3, dim, pwc_mesh_3, v_field%name)

      call insert(state_23(2), v_field_3, v_field_3%name)

    else

      call allocate(v_field_3, dim, mesh_3, v_field%name)

      call insert(state_23(1), v_field_3, v_field_3%name)

    end if

    call deallocate(v_field_3)

    if(shape%degree == 0) then

      call allocate(v_field_3, dim, pwc_mesh_3, v_field%name)

      call insert(state_3_split(2), v_field_3, v_field_3%name)

    else

      call allocate(v_field_3, dim, mesh_3, v_field%name)

      call insert(state_3_split(1), v_field_3, v_field_3%name)

    end if

    call insert(state_3, v_field_3, v_field_3%name)

    call deallocate(v_field_3)    

  end do

  do i = 1, tensor_field_count(state_1)

    t_field => extract_tensor_field(state_1, i)

    shape => ele_shape(t_field, 1)

    if(shape%degree == 0) then    

      assert(p0)      

      call insert(state_1_split(2), t_field, t_field%name)

      call insert(state_2_split(2), extract_scalar_field(state_2, t_field%name), t_field%name)

      call allocate(t_field_3, pwc_mesh_3, t_field%name)

      call insert(state_13(2), t_field_3, t_field_3%name)

    else

      call insert(state_1_split(1), t_field, t_field%name)

      call insert(state_2_split(1), extract_scalar_field(state_2, t_field%name), t_field%name)

      call allocate(t_field_3, mesh_3, t_field%name)

      call insert(state_13(1), t_field_3, t_field_3%name)

    end if

    call deallocate(t_field_3)

    if(shape%degree == 0) then

      call allocate(t_field_3, pwc_mesh_3, t_field%name)

      call insert(state_23(2), t_field_3, t_field_3%name)

    else

      call allocate(t_field_3, mesh_3, t_field%name)

      call insert(state_23(1), t_field_3, t_field_3%name)

    end if

    call deallocate(t_field_3)

    if(shape%degree == 0) then

      call allocate(t_field_3, pwc_mesh_3, t_field%name)

      call insert(state_3_split(2), t_field_3, t_field_3%name)

    else

      call allocate(t_field_3, mesh_3, t_field%name)

      call insert(state_3_split(1), t_field_3, t_field_3%name)

    end if

    call insert(state_3, t_field_3, t_field_3%name)

    call deallocate(t_field_3)    

  end do  

  call deallocate(pwc_mesh_3)

  call linear_interpolation(state_1_split(1), state_13(1), map = node_map_13)

  deallocate(node_map_13)  

  call linear_interpolation(state_2_split(1), state_23(1), map = node_map_23)

  deallocate(node_map_23)   

  if(p0) then 

    call linear_interpolation(state_1_split(2), state_13(2), map = ele_map_13)

    deallocate(ele_map_13)

    call linear_interpolation(state_2_split(2), state_23(2), map = ele_map_23)

    deallocate(ele_map_23)

  end if

  call deallocate(state_1_split)

  call deallocate(state_1)

  call deallocate(state_2_split)

  call deallocate(state_2)

  do i = 1, scalar_field_count(state_13(1))

    s_field_13 => extract_scalar_field(state_13(1), i)

    s_field_23 => extract_scalar_field(state_23(1), s_field_13%name)

    s_field_3 = extract_scalar_field(state_3_split(1), s_field_13%name)

    s_field_3%val = s_field_13%val - s_field_23%val

  end do

  do i = 1, vector_field_count(state_13(1))

    v_field_13 => extract_vector_field(state_13(1), i)

    if(v_field_13%name == "Coordinate") cycle

    v_field_23 => extract_vector_field(state_23(1), v_field_13%name)

    v_field_3 = extract_vector_field(state_3_split(1), v_field_13%name)

    do j = 1, dim

      v_field_3%val(i,:) = v_field_13%val(i,:) - v_field_23%val(i,:)

    end do

  end do

  do i = 1, tensor_field_count(state_13(1))

    t_field_13 => extract_tensor_field(state_13(1), i)

    t_field_23 => extract_tensor_field(state_23(1), t_field_13%name)

    t_field_3 = extract_tensor_field(state_3_split(1), t_field_13%name)    

    t_field_3%val = t_field_13%val - t_field_23%val

  end do

  if(p0) then

    do i = 1, scalar_field_count(state_13(2))

      s_field_13 => extract_scalar_field(state_13(2), i)

      s_field_23 => extract_scalar_field(state_23(2), s_field_13%name)

      s_field_3 = extract_scalar_field(state_3_split(2), s_field_13%name)

      s_field_3%val = s_field_13%val - s_field_23%val

    end do

    do i = 1, vector_field_count(state_13(2))

      v_field_13 => extract_vector_field(state_13(2), i)

      if(v_field_13%name == "Coordinate") cycle

      v_field_23 => extract_vector_field(state_23(2), v_field_13%name)

      v_field_3 = extract_vector_field(state_3_split(2), v_field_13%name)

      do j = 1, dim

        v_field_3%val(i,:) = v_field_13%val(i,:) - v_field_23%val(i,:)

      end do

    end do

    do i = 1, tensor_field_count(state_13(2))

      t_field_13 => extract_tensor_field(state_13(2), i)

      t_field_23 => extract_tensor_field(state_23(2), t_field_13%name)

      t_field_3 = extract_tensor_field(state_3_split(2), t_field_13%name)    

      t_field_3%val = t_field_13%val - t_field_23%val

    end do

  end if

  call deallocate(state_13)

  call deallocate(state_23)

  call deallocate(state_3_split)

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

  call deallocate(state_3)

  call print_references(0)

  ewrite(1, *) "Exiting supermesh_difference"

end subroutine supermesh_difference