1!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch 2!! Copyright (C) 2013 M. Oliveira 3!! 4!! This program is free software; you can redistribute it and/or modify 5!! it under the terms of the GNU General Public License as published by 6!! the Free Software Foundation; either version 2, or (at your option) 7!! any later version. 8!! 9!! This program is distributed in the hope that it will be useful, 10!! but WITHOUT ANY WARRANTY; without even the implied warranty of 11!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12!! GNU General Public License for more details. 13!! 14!! You should have received a copy of the GNU General Public License 15!! along with this program; if not, write to the Free Software 16!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 17!! 02110-1301, USA. 18!! 19 20#include "global.h" 21 22module box_oct_m 23 use global_oct_m 24 use messages_oct_m 25 use namespace_oct_m 26 use profiling_oct_m 27 28 implicit none 29 30 private 31 public :: & 32 box_t, & 33 box_create, & 34 box_end, & 35 box_inside, & 36 box_inside_vec, & 37 box_copy, & 38 box_get_center 39 40 41 integer, parameter, public :: & 42 BOX_SPHERE = 1, & 43 BOX_CYLINDER = 2, & 44 BOX_PARALLELEPIPED = 3 45 46 type box_t 47 private 48 49 integer :: dim 50 integer :: shape 51 52 FLOAT :: center(MAX_DIM) !< where is the box centered 53 54 FLOAT :: rsize !< the radius of the sphere or of the cylinder 55 FLOAT :: xsize !< the length of the cylinder in the x-direction 56 FLOAT :: lsize(MAX_DIM) !< half of the length of the parallelepiped in each direction. 57 end type box_t 58 59contains 60 61 !-------------------------------------------------------------- 62 subroutine box_create(box, shape, dim, sizes, center, namespace) 63 type(box_t), intent(out) :: box 64 integer, intent(in) :: shape 65 integer, intent(in) :: dim 66 FLOAT, intent(in) :: sizes(MAX_DIM) 67 FLOAT, intent(in) :: center(dim) 68 type(namespace_t), intent(in) :: namespace 69 70 PUSH_SUB(box_create) 71 72 box%shape = shape 73 box%dim = dim 74 box%center(1:dim) = center(1:dim) 75 76 select case (shape) 77 case (BOX_SPHERE) 78 box%rsize = sizes(1) 79 80 case (BOX_CYLINDER) 81 if (dim == 2) then 82 message(1) = "Cannot create a cylinder in 2D. Use sphere if you want a circle." 83 call messages_fatal(1, namespace=namespace) 84 end if 85 box%rsize = sizes(1) 86 box%xsize = sizes(2) 87 88 case (BOX_PARALLELEPIPED) 89 box%lsize(1:dim) = sizes(1:dim) 90 91 case default 92 message(1) = "Unknown box shape in box_create." 93 call messages_fatal(1, namespace=namespace) 94 95 end select 96 97 POP_SUB(box_create) 98 end subroutine box_create 99 100 !-------------------------------------------------------------- 101 subroutine box_end(box) 102 type(box_t), intent(inout) :: box 103 104 PUSH_SUB(box_end) 105 106 box%shape = 0 107 box%dim = 0 108 109 POP_SUB(box_end) 110 end subroutine box_end 111 112 !-------------------------------------------------------------- 113 !> Checks if a point is inside the box. 114 logical function box_inside(box, point) result(inside) 115 type(box_t), intent(in) :: box 116 FLOAT, intent(in) :: point(:) 117 118 FLOAT :: xx(1, 1:MAX_DIM) 119 logical :: inside2(1) 120 121 ! no push_sub because this function is called very frequently 122 123 xx(1, 1:box%dim) = point(1:box%dim) 124 125 call box_inside_vec(box, 1, xx, inside2) 126 inside = inside2(1) 127 128 end function box_inside 129 130 !-------------------------------------------------------------- 131 !> Checks if a vector of points are inside the box. 132 subroutine box_inside_vec(box, npoints, points, inside) 133 type(box_t), intent(in) :: box 134 integer, intent(in) :: npoints 135 FLOAT, intent(in) :: points(:, :) 136 logical, intent(out) :: inside(:) 137 138 integer :: ip 139 FLOAT, parameter :: DELTA = CNST(1e-12) 140 FLOAT :: llimit(MAX_DIM), ulimit(MAX_DIM) 141 FLOAT :: rr 142 FLOAT, allocatable :: xx(:, :) 143 144 ! no push_sub because this function is called very frequently 145 146 SAFE_ALLOCATE(xx(1:box%dim, 1:npoints)) 147 do ip = 1, npoints 148 xx(1:box%dim, ip) = points(ip, 1:box%dim) - box%center(1:box%dim) 149 end do 150 151 select case(box%shape) 152 case(BOX_SPHERE) 153 do ip = 1, npoints 154 inside(ip) = sum(xx(1:box%dim, ip)**2) <= (box%rsize + DELTA)**2 155 end do 156 157 case(BOX_CYLINDER) 158 do ip = 1, npoints 159 rr = sqrt(sum(xx(2:box%dim, ip)**2)) 160 inside(ip) = (rr <= box%rsize + DELTA .and. abs(xx(1, ip)) <= box%xsize + DELTA) 161 end do 162 163 case(BOX_PARALLELEPIPED) 164 llimit(1:box%dim) = -box%lsize(1:box%dim) - DELTA 165 ulimit(1:box%dim) = box%lsize(1:box%dim) + DELTA 166 167 do ip = 1, npoints 168 inside(ip) = all(xx(1:box%dim, ip) >= llimit(1:box%dim) .and. xx(1:box%dim, ip) <= ulimit(1:box%dim)) 169 end do 170 171 end select 172 173 SAFE_DEALLOCATE_A(xx) 174 175 end subroutine box_inside_vec 176 177 ! -------------------------------------------------------------- 178 recursive subroutine box_copy(boxout, boxin) 179 type(box_t), intent(out) :: boxout 180 type(box_t), intent(in) :: boxin 181 182 PUSH_SUB(box_copy) 183 184 boxout%shape = boxin%shape 185 boxout%dim = boxin%dim 186 boxout%center(1:boxin%dim) = boxin%center(1:boxin%dim) 187 boxout%rsize = boxin%rsize 188 boxout%xsize = boxin%xsize 189 boxout%lsize(1:boxin%dim) = boxin%lsize(1:boxin%dim) 190 191 POP_SUB(box_copy) 192 end subroutine box_copy 193 194 ! -------------------------------------------------------------- 195 pure function box_get_center(box) result(x) 196 type(box_t), intent(in) :: box 197 FLOAT, dimension(MAX_DIM) :: x 198 199 x(1:box%dim) = box%center(1:box%dim) 200 201 end function box_get_center 202end module box_oct_m 203 204 205!! Local Variables: 206!! mode: f90 207!! coding: utf-8 208!! End: 209