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