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