1!! Copyright (C) 2008 X. Andrade
2!!
3!! This program is free software; you can redistribute it and/or modify
4!! it under the terms of the GNU General Public License as published by
5!! the Free Software Foundation; either version 2, or (at your option)
6!! any later version.
7!!
8!! This program is distributed in the hope that it will be useful,
9!! but WITHOUT ANY WARRANTY; without even the implied warranty of
10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
11!! GNU General Public License for more details.
12!!
13!! You should have received a copy of the GNU General Public License
14!! along with this program; if not, write to the Free Software
15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
16!! 02110-1301, USA.
17!!
18
19#include "global.h"
20
21module stencil_oct_m
22  use global_oct_m
23  use messages_oct_m
24  use profiling_oct_m
25
26  implicit none
27
28  private
29  public ::                        &
30    stencil_t,                     &
31    stencil_allocate,              &
32    stencil_copy,                  &
33    stencil_end,                   &
34    stencil_init_center,           &
35    stencil_union
36
37  type stargeneral_arms_t
38    ! Components are public by default
39    integer          :: arms(1:3,1:3)
40    integer          :: narms
41  end type stargeneral_arms_t
42
43
44  type stencil_t
45    ! Components are public by default
46    integer          :: center
47    integer          :: size
48    integer          :: npoly
49    integer, pointer :: points(:, :)
50
51    ! The stargeneral arms
52    type(stargeneral_arms_t) :: stargeneral
53  end type stencil_t
54
55contains
56
57  !-------------------------------------------------------
58  subroutine stencil_allocate(this, size)
59    type(stencil_t), intent(out) :: this
60    integer,         intent(in)  :: size
61
62    PUSH_SUB(stencil_allocate)
63
64    this%size = size
65    this%npoly = size
66
67    SAFE_ALLOCATE(this%points(1:MAX_DIM, 1:size))
68
69    this%points = 0
70
71    POP_SUB(stencil_allocate)
72  end subroutine stencil_allocate
73
74  !-------------------------------------------------------
75  subroutine stencil_copy(input, output)
76    type(stencil_t), intent(in)  :: input
77    type(stencil_t), intent(out) :: output
78
79    PUSH_SUB(stencil_copy)
80
81    call stencil_allocate(output, input%size)
82    output%points(1:MAX_DIM, 1:output%size) = input%points(1:MAX_DIM, 1:output%size)
83    output%center = input%center
84    output%npoly = input%npoly
85
86    POP_SUB(stencil_copy)
87  end subroutine stencil_copy
88
89
90  !-------------------------------------------------------
91  subroutine stencil_end(this)
92    type(stencil_t), intent(inout) :: this
93
94    PUSH_SUB(stencil_end)
95
96    SAFE_DEALLOCATE_P(this%points)
97
98    POP_SUB(stencil_end)
99  end subroutine stencil_end
100
101
102  !-------------------------------------------------------
103  subroutine stencil_union(dim, st1, st2, stu)
104    integer,         intent(in)    :: dim
105    type(stencil_t), intent(inout) :: st1
106    type(stencil_t), intent(inout) :: st2
107    type(stencil_t), intent(inout) :: stu
108
109    integer :: idir, ii, jj, nstu
110    logical :: not_in_st1
111
112    PUSH_SUB(stencil_union)
113
114    call stencil_allocate(stu, st1%size + st2%size)
115
116    ! copy the first stencil
117    do ii = 1, st1%size
118      do idir = 1, dim
119        stu%points(idir, ii) = st1%points(idir, ii)
120      end do
121    end do
122
123    nstu = st1%size
124
125    do ii = 1, st2%size
126
127      not_in_st1 = .true.
128
129      ! check whether that point was already in the stencil
130      do jj = 1, st1%size
131        if(all(st1%points(1:dim, jj) == st2%points(1:dim, ii))) then
132          not_in_st1 = .false.
133          exit
134        end if
135      end do
136
137      if(not_in_st1) then !add it
138        nstu = nstu + 1
139        stu%points(1:dim, nstu) = st2%points(1:dim, ii)
140      end if
141
142    end do
143
144    stu%points(dim + 1:MAX_DIM, 1:nstu) = 0
145
146    stu%size = nstu
147
148    ! this is not defined for a union, which could be any combination. The
149    ! weights should already be known here
150    stu%npoly = -1
151
152    call stencil_init_center(stu)
153
154    POP_SUB(stencil_union)
155  end subroutine stencil_union
156
157
158  !-------------------------------------------------------
159  subroutine stencil_init_center(this)
160    type(stencil_t), intent(inout) :: this
161
162    integer :: ii
163
164    PUSH_SUB(stencil_init_center)
165
166    this%center = -1
167
168    do ii = 1, this%size
169      if(all(this%points(1:MAX_DIM, ii) == 0)) this%center = ii
170    end do
171
172    POP_SUB(stencil_init_center)
173  end subroutine stencil_init_center
174
175end module stencil_oct_m
176
177!! Local Variables:
178!! mode: f90
179!! coding: utf-8
180!! End:
181