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_cube_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_cube_size_lapl,        &
31    stencil_cube_get_lapl,         &
32    stencil_cube_polynomials_lapl, &
33    stencil_cube_size_grad,        &
34    stencil_cube_get_grad,         &
35    stencil_cube_polynomials_grad
36
37
38contains
39
40  ! ---------------------------------------------------------
41  integer function stencil_cube_size_lapl(dim, order)
42    integer, intent(in) :: dim
43    integer, intent(in) :: order
44
45    PUSH_SUB(stencil_cube_size_lapl)
46
47    stencil_cube_size_lapl = (2*order+1)**dim
48
49    POP_SUB(stencil_cube_size_lapl)
50  end function stencil_cube_size_lapl
51
52
53  ! ---------------------------------------------------------
54  subroutine stencil_cube_get_lapl(this, dim, order)
55    type(stencil_t), intent(out) :: this
56    integer,         intent(in)  :: dim
57    integer,         intent(in)  :: order
58
59    integer :: i, j, k, n
60
61    PUSH_SUB(stencil_cube_get_lapl)
62
63    call stencil_allocate(this, stencil_cube_size_lapl(dim, order))
64
65    n = 1
66    select case(dim)
67    case(1)
68      do i = -order, order
69        this%points(1, n) = i
70        n = n + 1
71      end do
72    case(2)
73      do i = -order, order
74        do j = -order, order
75          this%points(1, n) = i
76          this%points(2, n) = j
77          n = n + 1
78        end do
79      end do
80    case(3)
81      do i = -order, order
82        do j = -order, order
83          do k = -order, order
84            this%points(1, n) = i
85            this%points(2, n) = j
86            this%points(3, n) = k
87            n = n + 1
88          end do
89        end do
90      end do
91    end select
92
93    call stencil_init_center(this)
94
95    POP_SUB(stencil_cube_get_lapl)
96  end subroutine stencil_cube_get_lapl
97
98
99  ! ---------------------------------------------------------
100  subroutine stencil_cube_polynomials_lapl(dim, order, pol)
101    integer, intent(in)  :: dim
102    integer, intent(in)  :: order
103    integer, intent(out) :: pol(:,:) !< pol(dim, order)
104
105    integer :: i, j, k, n
106
107    PUSH_SUB(stencil_cube_polynomials_lapl)
108
109    n = 1
110    select case(dim)
111    case(1)
112      do i = 0, 2*order
113        pol(1, n) = i
114        n = n + 1
115      end do
116    case(2)
117      do i = 0, 2*order
118        do j = 0, 2*order
119          pol(1, n) = i
120          pol(2, n) = j
121          n = n + 1
122        end do
123      end do
124    case(3)
125      do i = 0, 2*order
126        do j = 0, 2*order
127          do k = 0, 2*order
128            pol(1, n) = i
129            pol(2, n) = j
130            pol(3, n) = k
131            n = n + 1
132          end do
133        end do
134      end do
135    end select
136
137    POP_SUB(stencil_cube_polynomials_lapl)
138  end subroutine stencil_cube_polynomials_lapl
139
140
141  !> Now come the gradient routines. As this stencil is the same for
142  !! the laplacian and the gradient, these routines just call the
143  !! corresponding ones for the laplacian
144
145  ! ---------------------------------------------------------
146  integer function stencil_cube_size_grad(dim, order)
147    integer, intent(in) :: dim
148    integer, intent(in) :: order
149
150    PUSH_SUB(stencil_cube_size_grad)
151
152    stencil_cube_size_grad = stencil_cube_size_lapl(dim, order)
153
154    POP_SUB(stencil_cube_size_grad)
155  end function stencil_cube_size_grad
156
157
158  ! ---------------------------------------------------------
159  subroutine stencil_cube_get_grad(this, dim, order)
160    type(stencil_t), intent(out) :: this
161    integer,         intent(in)  :: dim
162    integer,         intent(in)  :: order
163
164    PUSH_SUB(stencil_cube_get_grad)
165
166    call stencil_cube_get_lapl(this, dim, order)
167
168    POP_SUB(stencil_cube_get_grad)
169  end subroutine stencil_cube_get_grad
170
171
172  ! ---------------------------------------------------------
173  subroutine stencil_cube_polynomials_grad(dim, order, pol)
174    integer, intent(in)  :: dim
175    integer, intent(in)  :: order
176    integer, intent(out) :: pol(:,:) !< pol(sb%dim, order)
177
178    PUSH_SUB(stencil_cube_polynomials_grad)
179
180    call stencil_cube_polynomials_lapl(dim, order, pol)
181
182    POP_SUB(stencil_cube_polynomials_grad)
183  end subroutine stencil_cube_polynomials_grad
184
185end module stencil_cube_oct_m
186
187!! Local Variables:
188!! mode: f90
189!! coding: utf-8
190!! End:
191