1!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
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_starplus_oct_m
22  use global_oct_m
23  use messages_oct_m
24  use stencil_oct_m
25
26  implicit none
27
28  private
29  public ::                     &
30    stencil_starplus_size_lapl, &
31    stencil_starplus_get_lapl,  &
32    stencil_starplus_pol_lapl,  &
33    stencil_starplus_size_grad, &
34    stencil_starplus_get_grad,  &
35    stencil_starplus_pol_grad
36
37contains
38
39  ! ---------------------------------------------------------
40  integer function stencil_starplus_size_lapl(dim, order) result(n)
41    integer, intent(in) :: dim
42    integer, intent(in) :: order
43
44    PUSH_SUB(stencil_starplus_size_lapl)
45
46    n = 2*dim*order + 1
47    if(dim == 2) n = n + 12
48    if(dim == 3) n = n + 44
49
50    POP_SUB(stencil_starplus_size_lapl)
51  end function stencil_starplus_size_lapl
52
53
54  ! ---------------------------------------------------------
55  integer function stencil_starplus_size_grad(dim, order) result(n)
56    integer, intent(in) :: dim
57    integer, intent(in) :: order
58
59    PUSH_SUB(stencil_starplus_size_grad)
60
61    n = 2*order + 1
62    if(dim == 2) n = n + 2
63    if(dim == 3) n = n + 4
64
65    POP_SUB(stencil_starplus_size_grad)
66  end function stencil_starplus_size_grad
67
68
69  ! ---------------------------------------------------------
70  subroutine stencil_starplus_get_lapl(this, dim, order)
71    type(stencil_t), intent(out) :: this
72    integer,         intent(in)  :: dim
73    integer,         intent(in)  :: order
74
75    integer :: i, j, n
76
77    PUSH_SUB(stencil_starplus_get_lapl)
78
79    call stencil_allocate(this, stencil_starplus_size_lapl(dim, order))
80
81    n = 1
82    select case(dim)
83    case(1)
84      n = 1
85      do i = 1, dim
86        do j = -order, order
87          if(j == 0) cycle
88          n = n + 1
89          this%points(i, n) = j
90        end do
91      end do
92    case(2)
93      n = 1
94      do i = 1, dim
95        do j = -order, order
96          if(j == 0) cycle
97          n = n + 1
98          this%points(i, n) = j
99        end do
100      end do
101      n = n + 1; this%points(1:2, n) = (/ -2,  1 /)
102      n = n + 1; this%points(1:2, n) = (/ -2, -1 /)
103      n = n + 1; this%points(1:2, n) = (/ -1,  2 /)
104      n = n + 1; this%points(1:2, n) = (/ -1,  1 /)
105      n = n + 1; this%points(1:2, n) = (/ -1, -1 /)
106      n = n + 1; this%points(1:2, n) = (/ -1, -2 /)
107      n = n + 1; this%points(1:2, n) = (/  1,  2 /)
108      n = n + 1; this%points(1:2, n) = (/  1,  1 /)
109      n = n + 1; this%points(1:2, n) = (/  1, -1 /)
110      n = n + 1; this%points(1:2, n) = (/  1, -2 /)
111      n = n + 1; this%points(1:2, n) = (/  2,  1 /)
112      n = n + 1; this%points(1:2, n) = (/  2, -1 /)
113    case(3)
114      n = 1
115      do i = 1, dim
116        do j = -order, order
117          if(j == 0) cycle
118          n = n + 1
119          this%points(i, n) = j
120        end do
121      end do
122      n = n + 1; this%points(1:3, n) = (/ -2,  1, 0 /)
123      n = n + 1; this%points(1:3, n) = (/ -2, -1, 0 /)
124      n = n + 1; this%points(1:3, n) = (/ -1,  2, 0 /)
125      n = n + 1; this%points(1:3, n) = (/ -1,  1, 0 /)
126      n = n + 1; this%points(1:3, n) = (/ -1, -1, 0 /)
127      n = n + 1; this%points(1:3, n) = (/ -1, -2, 0 /)
128      n = n + 1; this%points(1:3, n) = (/  1,  2, 0 /)
129      n = n + 1; this%points(1:3, n) = (/  1,  1, 0 /)
130      n = n + 1; this%points(1:3, n) = (/  1, -1, 0 /)
131      n = n + 1; this%points(1:3, n) = (/  1, -2, 0 /)
132      n = n + 1; this%points(1:3, n) = (/  2,  1, 0 /)
133      n = n + 1; this%points(1:3, n) = (/  2, -1, 0 /)
134
135      n = n + 1; this%points(1:3, n) = (/ -2, 0,  1 /)
136      n = n + 1; this%points(1:3, n) = (/ -2, 0, -1 /)
137      n = n + 1; this%points(1:3, n) = (/ -1, 0,  2 /)
138      n = n + 1; this%points(1:3, n) = (/ -1, 0,  1 /)
139      n = n + 1; this%points(1:3, n) = (/ -1, 0, -1 /)
140      n = n + 1; this%points(1:3, n) = (/ -1, 0, -2 /)
141      n = n + 1; this%points(1:3, n) = (/  1, 0,  2 /)
142      n = n + 1; this%points(1:3, n) = (/  1, 0,  1 /)
143      n = n + 1; this%points(1:3, n) = (/  1, 0, -1 /)
144      n = n + 1; this%points(1:3, n) = (/  1, 0, -2 /)
145      n = n + 1; this%points(1:3, n) = (/  2, 0,  1 /)
146      n = n + 1; this%points(1:3, n) = (/  2, 0, -1 /)
147
148      n = n + 1; this%points(1:3, n) = (/ 0, -2,  1 /)
149      n = n + 1; this%points(1:3, n) = (/ 0, -2, -1 /)
150      n = n + 1; this%points(1:3, n) = (/ 0, -1,  2 /)
151      n = n + 1; this%points(1:3, n) = (/ 0, -1,  1 /)
152      n = n + 1; this%points(1:3, n) = (/ 0, -1, -1 /)
153      n = n + 1; this%points(1:3, n) = (/ 0, -1, -2 /)
154      n = n + 1; this%points(1:3, n) = (/ 0,  1,  2 /)
155      n = n + 1; this%points(1:3, n) = (/ 0,  1,  1 /)
156      n = n + 1; this%points(1:3, n) = (/ 0,  1, -1 /)
157      n = n + 1; this%points(1:3, n) = (/ 0,  1, -2 /)
158      n = n + 1; this%points(1:3, n) = (/ 0,  2,  1 /)
159      n = n + 1; this%points(1:3, n) = (/ 0,  2, -1 /)
160
161      n = n + 1; this%points(1:3, n) = (/ -1, -1, -1 /)
162      n = n + 1; this%points(1:3, n) = (/ -1, -1,  1 /)
163      n = n + 1; this%points(1:3, n) = (/ -1,  1, -1 /)
164      n = n + 1; this%points(1:3, n) = (/ -1,  1,  1 /)
165      n = n + 1; this%points(1:3, n) = (/  1, -1, -1 /)
166      n = n + 1; this%points(1:3, n) = (/  1, -1,  1 /)
167      n = n + 1; this%points(1:3, n) = (/  1,  1, -1 /)
168      n = n + 1; this%points(1:3, n) = (/  1,  1,  1 /)
169
170    end select
171
172    call stencil_init_center(this)
173
174    POP_SUB(stencil_starplus_get_lapl)
175  end subroutine stencil_starplus_get_lapl
176
177
178  ! ---------------------------------------------------------
179  subroutine stencil_starplus_get_grad(this, dim, dir, order)
180    type(stencil_t), intent(out) :: this
181    integer,         intent(in)  :: dim
182    integer,         intent(in)  :: dir
183    integer,         intent(in)  :: order
184
185    integer :: i, n, j
186
187    PUSH_SUB(stencil_starplus_get_grad)
188
189    call stencil_allocate(this, stencil_starplus_size_grad(dim, order))
190
191    n = 1
192    do i = -order, order
193      this%points(dir, n) = i
194      n = n + 1
195    end do
196    do j = 1, dim
197      if(j==dir) cycle
198      this%points(j, n) = -1
199      n = n + 1
200      this%points(j, n) =  1
201      n = n + 1
202    end do
203
204    call stencil_init_center(this)
205
206    POP_SUB(stencil_starplus_get_grad)
207  end subroutine stencil_starplus_get_grad
208
209
210  ! ---------------------------------------------------------
211  subroutine stencil_starplus_pol_lapl(dim, order, pol)
212    integer, intent(in)  :: dim
213    integer, intent(in)  :: order
214    integer, intent(out) :: pol(:,:) !< pol(dim, order)
215
216    integer :: i, j, n
217
218    PUSH_SUB(stencil_starplus_pol_lapl)
219
220    n = 1
221    select case(dim)
222    case(1)
223      n = 1
224      pol(:,:) = 0
225      do i = 1, dim
226        do j = 1, 2*order
227          n = n + 1
228          pol(i, n) = j
229        end do
230      end do
231    case(2)
232      n = 1
233      pol(:,:) = 0
234      do i = 1, dim
235        do j = 1, 2*order
236          n = n + 1
237          pol(i, n) = j
238        end do
239      end do
240      n = n + 1; pol(1:2, n) = (/ 1, 1 /)
241      n = n + 1; pol(1:2, n) = (/ 1, 2 /)
242      n = n + 1; pol(1:2, n) = (/ 1, 3 /)
243      n = n + 1; pol(1:2, n) = (/ 1, 4 /)
244      n = n + 1; pol(1:2, n) = (/ 2, 1 /)
245      n = n + 1; pol(1:2, n) = (/ 2, 2 /)
246      n = n + 1; pol(1:2, n) = (/ 2, 3 /)
247      n = n + 1; pol(1:2, n) = (/ 2, 4 /)
248      n = n + 1; pol(1:2, n) = (/ 3, 1 /)
249      n = n + 1; pol(1:2, n) = (/ 3, 2 /)
250      n = n + 1; pol(1:2, n) = (/ 4, 1 /)
251      n = n + 1; pol(1:2, n) = (/ 4, 2 /)
252    case(3)
253      n = 1
254      pol(:,:) = 0
255      do i = 1, dim
256        do j = 1, 2*order
257          n = n + 1
258          pol(i, n) = j
259        end do
260      end do
261
262      n = n + 1; pol(1:3, n) = (/ 1, 1, 0 /)
263      n = n + 1; pol(1:3, n) = (/ 1, 2, 0 /)
264      n = n + 1; pol(1:3, n) = (/ 1, 3, 0 /)
265      n = n + 1; pol(1:3, n) = (/ 1, 4, 0 /)
266      n = n + 1; pol(1:3, n) = (/ 2, 1, 0 /)
267      n = n + 1; pol(1:3, n) = (/ 2, 2, 0 /)
268      n = n + 1; pol(1:3, n) = (/ 2, 3, 0 /)
269      n = n + 1; pol(1:3, n) = (/ 2, 4, 0 /)
270      n = n + 1; pol(1:3, n) = (/ 3, 1, 0 /)
271      n = n + 1; pol(1:3, n) = (/ 3, 2, 0 /)
272      n = n + 1; pol(1:3, n) = (/ 4, 1, 0 /)
273      n = n + 1; pol(1:3, n) = (/ 4, 2, 0 /)
274
275      n = n + 1; pol(1:3, n) = (/ 1, 0, 1 /)
276      n = n + 1; pol(1:3, n) = (/ 1, 0, 2 /)
277      n = n + 1; pol(1:3, n) = (/ 1, 0, 3 /)
278      n = n + 1; pol(1:3, n) = (/ 1, 0, 4 /)
279      n = n + 1; pol(1:3, n) = (/ 2, 0, 1 /)
280      n = n + 1; pol(1:3, n) = (/ 2, 0, 2 /)
281      n = n + 1; pol(1:3, n) = (/ 2, 0, 3 /)
282      n = n + 1; pol(1:3, n) = (/ 2, 0, 4 /)
283      n = n + 1; pol(1:3, n) = (/ 3, 0, 1 /)
284      n = n + 1; pol(1:3, n) = (/ 3, 0, 2 /)
285      n = n + 1; pol(1:3, n) = (/ 4, 0, 1 /)
286      n = n + 1; pol(1:3, n) = (/ 4, 0, 2 /)
287
288      n = n + 1; pol(1:3, n) = (/ 0, 1, 1 /)
289      n = n + 1; pol(1:3, n) = (/ 0, 1, 2 /)
290      n = n + 1; pol(1:3, n) = (/ 0, 1, 3 /)
291      n = n + 1; pol(1:3, n) = (/ 0, 1, 4 /)
292      n = n + 1; pol(1:3, n) = (/ 0, 2, 1 /)
293      n = n + 1; pol(1:3, n) = (/ 0, 2, 2 /)
294      n = n + 1; pol(1:3, n) = (/ 0, 2, 3 /)
295      n = n + 1; pol(1:3, n) = (/ 0, 2, 4 /)
296      n = n + 1; pol(1:3, n) = (/ 0, 3, 1 /)
297      n = n + 1; pol(1:3, n) = (/ 0, 3, 2 /)
298      n = n + 1; pol(1:3, n) = (/ 0, 4, 1 /)
299      n = n + 1; pol(1:3, n) = (/ 0, 4, 2 /)
300
301      n = n + 1; pol(1:3, n) = (/ 1, 1, 1 /)
302      n = n + 1; pol(1:3, n) = (/ 1, 1, 2 /)
303      n = n + 1; pol(1:3, n) = (/ 1, 2, 1 /)
304      n = n + 1; pol(1:3, n) = (/ 1, 2, 2 /)
305      n = n + 1; pol(1:3, n) = (/ 2, 1, 1 /)
306      n = n + 1; pol(1:3, n) = (/ 2, 1, 2 /)
307      n = n + 1; pol(1:3, n) = (/ 2, 2, 1 /)
308      n = n + 1; pol(1:3, n) = (/ 2, 2, 2 /)
309
310    end select
311
312    POP_SUB(stencil_starplus_pol_lapl)
313  end subroutine stencil_starplus_pol_lapl
314
315
316  ! ---------------------------------------------------------
317  subroutine stencil_starplus_pol_grad(dim, dir, order, pol)
318    integer, intent(in)  :: dim
319    integer, intent(in)  :: dir
320    integer, intent(in)  :: order
321    integer, intent(out) :: pol(:,:) !< pol(dim, order)
322
323    integer :: j, n
324
325    PUSH_SUB(stencil_starplus_pol_grad)
326
327    pol(:,:) = 0
328    do j = 0, 2*order
329      pol(dir, j+1) = j
330    end do
331    n = 2*order + 1
332
333    select case(dim)
334    case(2)
335      select case(dir)
336      case(1)
337        n = n + 1; pol(1:2, n) = (/ 0, 1 /)
338        n = n + 1; pol(1:2, n) = (/ 0, 2 /)
339      case(2)
340        n = n + 1; pol(1:2, n) = (/ 1, 0 /)
341        n = n + 1; pol(1:2, n) = (/ 2, 0 /)
342      end select
343    case(3)
344      select case(dir)
345      case(1)
346        n = n + 1; pol(1:3, n) = (/ 0, 1, 0 /)
347        n = n + 1; pol(1:3, n) = (/ 0, 2, 0 /)
348        n = n + 1; pol(1:3, n) = (/ 0, 0, 1 /)
349        n = n + 1; pol(1:3, n) = (/ 0, 0, 2 /)
350      case(2)
351        n = n + 1; pol(1:3, n) = (/ 1, 0, 0 /)
352        n = n + 1; pol(1:3, n) = (/ 2, 0, 0 /)
353        n = n + 1; pol(1:3, n) = (/ 0, 0, 1 /)
354        n = n + 1; pol(1:3, n) = (/ 0, 0, 2 /)
355      case(3)
356        n = n + 1; pol(1:3, n) = (/ 1, 0, 0 /)
357        n = n + 1; pol(1:3, n) = (/ 2, 0, 0 /)
358        n = n + 1; pol(1:3, n) = (/ 0, 1, 0 /)
359        n = n + 1; pol(1:3, n) = (/ 0, 2, 0 /)
360      end select
361    end select
362
363    POP_SUB(stencil_starplus_pol_grad)
364  end subroutine stencil_starplus_pol_grad
365
366end module stencil_starplus_oct_m
367
368!! Local Variables:
369!! mode: f90
370!! coding: utf-8
371!! End:
372