stencil.F90

Go to the documentation of this file.

!! Copyright (C) 2008 X. Andrade
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program 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 General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
 
module stencil_oct_m
  use debug_oct_m
  use global_oct_m
  use messages_oct_m
  use profiling_oct_m
 
  implicit none
 
  private
  public ::                        &
    stencil_t,                     &
    stencil_allocate,              &
    stencil_copy,                  &
    stencil_end,                   &
    stencil_init_center,           &
    stencil_union
 
  type stargeneral_arms_t
    ! Components are public by default
    integer          :: arms(1:3,1:3) = reshape([0,0,0,0,0,0,0,0,0], (/3, 3/))
    integer          :: narms = 0
  end type stargeneral_arms_t
 
 
  type stencil_t
    ! Components are public by default
    integer, private     :: dim
    integer              :: center
    !                                     !! A value of -1 indicates that the central point is not in the stencil.
    integer              :: size
    integer, allocatable :: points(:, :)
 
    ! The stargeneral arms
    type(stargeneral_arms_t) :: stargeneral
  end type stencil_t
 
contains
 
  !-------------------------------------------------------
  subroutine stencil_allocate(this, dim, size)
    type(stencil_t), intent(inout) :: this
    integer,         intent(in)    :: dim
    integer,         intent(in)    :: size
 
    push_sub(stencil_allocate)
 
    this%dim = dim
    this%size = size
 
    safe_allocate(this%points(1:this%dim, 1:size))
 
    this%points = 0
 
    pop_sub(stencil_allocate)
  end subroutine stencil_allocate
 
  !-------------------------------------------------------
  subroutine stencil_copy(input, output)
    type(stencil_t), intent(in)  :: input
    type(stencil_t), intent(out) :: output
 
    push_sub(stencil_copy)
 
    call stencil_allocate(output, input%dim, input%size)
    output%points = input%points
    output%center = input%center
 
    output%stargeneral%narms = input%stargeneral%narms
    output%stargeneral%arms = input%stargeneral%arms
 
    pop_sub(stencil_copy)
  end subroutine stencil_copy
 
 
  !-------------------------------------------------------
  subroutine stencil_end(this)
    type(stencil_t), intent(inout) :: this
 
    push_sub(stencil_end)
 
    safe_deallocate_a(this%points)
 
    pop_sub(stencil_end)
  end subroutine stencil_end
 
 
  !-------------------------------------------------------
  subroutine stencil_union(st1, st2, stu)
    type(stencil_t), intent(in)    :: st1
    type(stencil_t), intent(in)    :: st2
    type(stencil_t), intent(inout) :: stu
 
    integer :: ii, jj, nstu
    logical :: not_in_st1
 
    push_sub(stencil_union)
 
    assert(st1%dim == st2%dim)
 
    call stencil_allocate(stu, st1%dim, st1%size + st2%size)
 
    ! copy the first stencil
    do ii = 1, st1%size
      stu%points(:, ii) = st1%points(:, ii)
    end do
 
    nstu = st1%size
 
    do ii = 1, st2%size
 
      not_in_st1 = .true.
 
      ! check whether that point was already in the stencil
      do jj = 1, nstu
        if (all(stu%points(:, jj) == st2%points(:, ii))) then
          not_in_st1 = .false.
          exit
        end if
      end do
 
      if (not_in_st1) then !add it
        nstu = nstu + 1
        stu%points(:, nstu) = st2%points(:, ii)
      end if
 
    end do
 
    stu%size = nstu
 
    call stencil_init_center(stu)
 
    pop_sub(stencil_union)
  end subroutine stencil_union
 
 
  !-------------------------------------------------------
  subroutine stencil_init_center(this)
    type(stencil_t), intent(inout) :: this
 
    integer :: ii
 
    push_sub(stencil_init_center)
 
    this%center = -1
 
    do ii = 1, this%size
      if (all(this%points(:, ii) == 0)) this%center = ii
    end do
 
    pop_sub(stencil_init_center)
  end subroutine stencil_init_center
 
end module stencil_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: