tensorprod/m_tensorprod_indices.F90
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | !> @brief Module for handling tensor product indices and lists of indices | ||
| 2 | module m_tensorprod_indices | ||
| 3 | use m_common, only: wp | ||
| 4 | implicit none | ||
| 5 | |||
| 6 | private | ||
| 7 | public :: TensorProdIndices, TensorProdIndexList | ||
| 8 | public :: Box | ||
| 9 | public :: size, intersect, size_2d | ||
| 10 | |||
| 11 | !> @brief Type to hold tensor product indices | ||
| 12 | !> | ||
| 13 | !> This type holds the indices for a tensor product space, where the indices themselves are also a tensor product of 1D indices | ||
| 14 | type :: TensorProdIndices | ||
| 15 | !> Start index in x-direction | ||
| 16 | integer :: i0 | ||
| 17 | !> End index in x-direction | ||
| 18 | integer :: i1 | ||
| 19 | !> Start index in y-direction | ||
| 20 | integer :: j0 | ||
| 21 | !> End index in y-direction | ||
| 22 | integer :: j1 | ||
| 23 | !> Start index in z-direction | ||
| 24 | integer :: k0 | ||
| 25 | !> End index in z-direction | ||
| 26 | integer :: k1 | ||
| 27 | !> Total number in x-direction | ||
| 28 | integer :: total_nr_x | ||
| 29 | !> Total number in y-direction | ||
| 30 | integer :: total_nr_y | ||
| 31 | !> Total number in z-direction | ||
| 32 | integer :: total_nr_z | ||
| 33 | !> If true, then the indices represent a complement of the specified indices | ||
| 34 | !> (e.g. {0, ..., i0-1} U {i1+1, ..., total_nr_x-1} for x-direction) | ||
| 35 | logical :: complement | ||
| 36 | contains | ||
| 37 | procedure :: init => init_indices | ||
| 38 | procedure :: empty => empty_indices | ||
| 39 | end type | ||
| 40 | |||
| 41 | !> @brief Type to hold lists of tensor product indices | ||
| 42 | !> | ||
| 43 | !> This type holds a list of indices for a tensor product space, where the indices themselves are not necessarily a tensor product | ||
| 44 | !> of 1D indices | ||
| 45 | type TensorProdIndexList | ||
| 46 | |||
| 47 | !> Linear index | ||
| 48 | integer, allocatable :: lin(:) | ||
| 49 | |||
| 50 | !> Cartesian index i | ||
| 51 | integer, allocatable :: i(:) | ||
| 52 | !> Cartesian index j | ||
| 53 | integer, allocatable :: j(:) | ||
| 54 | !> Cartesian index k | ||
| 55 | integer, allocatable :: k(:) | ||
| 56 | |||
| 57 | !> The corresponding TensorProdIndices | ||
| 58 | type(TensorProdIndices) :: tp_inds | ||
| 59 | contains | ||
| 60 | procedure :: init => init_lists | ||
| 61 | procedure :: destroy => destroy_indices | ||
| 62 | procedure, private :: get_local_index_2d | ||
| 63 | procedure, private :: get_local_index_3d | ||
| 64 | generic :: get_local_index => get_local_index_2d, get_local_index_3d | ||
| 65 | procedure :: get_bounds => get_bounds_list | ||
| 66 | end type | ||
| 67 | |||
| 68 | !> @brief Type to represent a 3D box | ||
| 69 | type Box | ||
| 70 | !> Left most x-coordinate | ||
| 71 | real(wp) :: x0 = 0._wp | ||
| 72 | !> Right most x-coordinate | ||
| 73 | real(wp) :: x1 = -1._wp | ||
| 74 | !> Left most y-coordinate | ||
| 75 | real(wp) :: y0 = 0._wp | ||
| 76 | !> Right most y-coordinate | ||
| 77 | real(wp) :: y1 = -1._wp | ||
| 78 | !> Left most z-coordinate | ||
| 79 | real(wp) :: z0 = 0._wp | ||
| 80 | !> Right most z-coordinate | ||
| 81 | real(wp) :: z1 = -1._wp | ||
| 82 | contains | ||
| 83 | procedure :: is_inside | ||
| 84 | end type Box | ||
| 85 | |||
| 86 | !> @brief Get the size of the tensor product indices | ||
| 87 | interface size | ||
| 88 | module procedure size_tpinds, size_index_list | ||
| 89 | end interface | ||
| 90 | |||
| 91 | interface size_2d | ||
| 92 | module procedure size_tp_inds_2d, size_index_list_2d | ||
| 93 | end interface | ||
| 94 | |||
| 95 | contains | ||
| 96 | !> @brief Initialize the TensorProdIndices | ||
| 97 | !> | ||
| 98 | !> @param this[inout] The TensorProdIndices to initialize | ||
| 99 | !> @param i0[in] Start index in x-direction | ||
| 100 | !> @param i1[in] End index in x-direction | ||
| 101 | !> @param j0[in] Start index in y-direction | ||
| 102 | !> @param j1[in] End index in y-direction | ||
| 103 | !> @param k0[in] Start index in z-direction | ||
| 104 | !> @param k1[in] End index in z-direction | ||
| 105 | !> @param bsplines[in] The B-spline spaces for each direction | ||
| 106 | !> @param complement[in] _(optional)_ If true, the indices represent a complement of the specified indices (default is false) | ||
| 107 | !> @param intervals[in] _(optional)_ If true, the indices represent intervals instead of B-splines (default is false) | ||
| 108 | 313108 | pure subroutine init_indices(this, i0, i1, j0, j1, k0, k1, bsplines, complement, intervals) | |
| 109 | use m_bspline_basis, only: BSplineSpace | ||
| 110 | implicit none | ||
| 111 | |||
| 112 | class(TensorProdIndices), intent(inout) :: this | ||
| 113 | integer, intent(in) :: i0, i1, j0, j1, k0, k1 | ||
| 114 | type(BSplineSpace), intent(in) :: bsplines(3) | ||
| 115 | logical, intent(in), optional :: complement | ||
| 116 | logical, intent(in), optional :: intervals | ||
| 117 | |||
| 118 | logical :: intervals_ | ||
| 119 | |||
| 120 | 313108 | this%complement = .false. | |
| 121 |
1/2✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 313108 times.
|
313108 | if (present(complement)) this%complement = complement |
| 122 |
2/2✓ Branch 4 → 5 taken 312891 times.
✓ Branch 4 → 6 taken 217 times.
|
313108 | if (i0 <= i1) then |
| 123 | 312891 | this%i0 = i0 | |
| 124 | 312891 | this%i1 = i1 | |
| 125 | else | ||
| 126 | 217 | this%i0 = 0 | |
| 127 | 217 | this%i1 = -1 | |
| 128 | end if | ||
| 129 |
2/2✓ Branch 7 → 8 taken 311277 times.
✓ Branch 7 → 9 taken 1831 times.
|
313108 | if (j0 <= j1) then |
| 130 | 311277 | this%j0 = j0 | |
| 131 | 311277 | this%j1 = j1 | |
| 132 | else | ||
| 133 | 1831 | this%j0 = 0 | |
| 134 | 1831 | this%j1 = -1 | |
| 135 | end if | ||
| 136 |
2/2✓ Branch 10 → 11 taken 297538 times.
✓ Branch 10 → 12 taken 15570 times.
|
313108 | if (k0 <= k1) then |
| 137 | 297538 | this%k0 = k0 | |
| 138 | 297538 | this%k1 = k1 | |
| 139 | else | ||
| 140 | 15570 | this%k0 = 0 | |
| 141 | 15570 | this%k1 = -1 | |
| 142 | end if | ||
| 143 | |||
| 144 | intervals_ = .false. | ||
| 145 |
2/2✓ Branch 13 → 14 taken 122837 times.
✓ Branch 13 → 16 taken 190271 times.
|
313108 | if (present(intervals)) intervals_ = intervals |
| 146 |
1/2✓ Branch 14 → 15 taken 122837 times.
✗ Branch 14 → 16 not taken.
|
122837 | if (intervals_) then |
| 147 | 122837 | this%total_nr_x = bsplines(1)%nr_intervals | |
| 148 | 122837 | this%total_nr_y = bsplines(2)%nr_intervals | |
| 149 | 122837 | this%total_nr_z = bsplines(3)%nr_intervals | |
| 150 | else | ||
| 151 | 190271 | this%total_nr_x = bsplines(1)%nr_bsplines | |
| 152 | 190271 | this%total_nr_y = bsplines(2)%nr_bsplines | |
| 153 | 190271 | this%total_nr_z = bsplines(3)%nr_bsplines | |
| 154 | end if | ||
| 155 | |||
| 156 | 313108 | end subroutine init_indices | |
| 157 | |||
| 158 | !> @brief Initialize the TensorProdIndices to empty | ||
| 159 | !> | ||
| 160 | !> @param this[inout] The constructed TensorProdIndices | ||
| 161 | 10452110 | subroutine empty_indices(this) | |
| 162 | class(TensorProdIndices), intent(inout) :: this | ||
| 163 | |||
| 164 | 10452110 | this%i0 = 0 | |
| 165 | 10452110 | this%i1 = -1 | |
| 166 | 10452110 | this%j0 = 0 | |
| 167 | 10452110 | this%j1 = -1 | |
| 168 | 10452110 | this%k0 = 0 | |
| 169 | 10452110 | this%k1 = -1 | |
| 170 | 10452110 | this%total_nr_x = 0 | |
| 171 | 10452110 | this%total_nr_y = 0 | |
| 172 | 10452110 | this%total_nr_z = 0 | |
| 173 | 10452110 | this%complement = .false. | |
| 174 | |||
| 175 | 10452110 | end subroutine empty_indices | |
| 176 | |||
| 177 | !> @brief Initialize the TensorProdIndexList from the given TensorProdIndices | ||
| 178 | !> | ||
| 179 | !> @param this[inout] The TensorProdIndexList to initialize | ||
| 180 | !> @param inds[in] The TensorProdIndices to use for initialization | ||
| 181 | !> @param my_inds[in] The TensorProdIndices of the current rank (used for determining the global linear index, as well as the | ||
| 182 | !> complement when applicable) | ||
| 183 | !> @param l0[in] The starting linear index for the current rank (used for determining the global linear index) | ||
| 184 | !> | ||
| 185 | !> @note If inds%complement is true, the indices will be initialized as the complement of the specified indices within my_inds | ||
| 186 | 2625 | subroutine init_lists(this, inds, my_inds, l0) | |
| 187 | implicit none | ||
| 188 | |||
| 189 | class(TensorProdIndexList), intent(inout) :: this | ||
| 190 | type(TensorProdIndices), intent(in) :: inds, my_inds | ||
| 191 | integer, intent(in) :: l0 | ||
| 192 | |||
| 193 | integer :: sze, i, j, k, row_global, row_local, my_sze_x, my_sze_y, my_sze_z | ||
| 194 | |||
| 195 | ! NOTE: complement is understood as the complement of the indices within the responsible B-spline indices of this rank | ||
| 196 | 2625 | my_sze_x = my_inds%i1 - my_inds%i0 + 1 | |
| 197 | 2625 | my_sze_y = my_inds%j1 - my_inds%j0 + 1 | |
| 198 | 2625 | my_sze_z = my_inds%k1 - my_inds%k0 + 1 | |
| 199 | |||
| 200 | 2625 | sze = (inds%i1 - inds%i0 + 1) * (inds%j1 - inds%j0 + 1) * (inds%k1 - inds%k0 + 1) | |
| 201 |
2/2✓ Branch 2 → 3 taken 385 times.
✓ Branch 2 → 4 taken 2240 times.
|
2625 | if (inds%complement) then |
| 202 | 385 | sze = my_sze_x * my_sze_y * my_sze_z - sze | |
| 203 | end if | ||
| 204 | |||
| 205 |
8/14✓ Branch 4 → 5 taken 1908 times.
✓ Branch 4 → 12 taken 717 times.
✓ Branch 5 → 6 taken 1908 times.
✗ Branch 5 → 11 not taken.
✓ Branch 6 → 7 taken 1908 times.
✗ Branch 6 → 11 not taken.
✓ Branch 7 → 8 taken 1908 times.
✗ Branch 7 → 11 not taken.
✓ Branch 8 → 9 taken 1908 times.
✗ Branch 8 → 11 not taken.
✓ Branch 9 → 10 taken 1908 times.
✗ Branch 9 → 11 not taken.
✗ Branch 10 → 11 not taken.
✓ Branch 10 → 12 taken 1908 times.
|
2625 | if (sze > 0 .and. (inds%i0 < my_inds%i0 .or. & |
| 206 | inds%i1 > my_inds%i1 .or. & | ||
| 207 | inds%j0 < my_inds%j0 .or. & | ||
| 208 | inds%j1 > my_inds%j1 .or. & | ||
| 209 | inds%k0 < my_inds%k0 .or. & | ||
| 210 | inds%k1 > my_inds%k1)) then | ||
| 211 | ✗ | error stop "TensorProdIndexList::init_lists provided indices must be a subset of responsible B-spline indices of this rank." | |
| 212 | end if | ||
| 213 | |||
| 214 |
1/2✗ Branch 12 → 13 not taken.
✓ Branch 12 → 14 taken 2625 times.
|
2625 | if (allocated(this%lin)) deallocate (this%lin) |
| 215 |
1/2✗ Branch 14 → 15 not taken.
✓ Branch 14 → 16 taken 2625 times.
|
2625 | if (allocated(this%i)) deallocate (this%i) |
| 216 |
1/2✗ Branch 16 → 17 not taken.
✓ Branch 16 → 18 taken 2625 times.
|
2625 | if (allocated(this%j)) deallocate (this%j) |
| 217 |
1/2✗ Branch 18 → 19 not taken.
✓ Branch 18 → 20 taken 2625 times.
|
2625 | if (allocated(this%k)) deallocate (this%k) |
| 218 | |||
| 219 |
9/14✓ Branch 20 → 21 taken 717 times.
✓ Branch 20 → 22 taken 1908 times.
✓ Branch 22 → 21 taken 1908 times.
✗ Branch 22 → 23 not taken.
✓ Branch 23 → 24 taken 2625 times.
✗ Branch 23 → 25 not taken.
✓ Branch 25 → 26 taken 1908 times.
✓ Branch 25 → 27 taken 717 times.
✗ Branch 27 → 28 not taken.
✓ Branch 27 → 29 taken 2625 times.
✗ Branch 29 → 30 not taken.
✓ Branch 29 → 31 taken 2625 times.
✗ Branch 31 → 32 not taken.
✓ Branch 31 → 33 taken 2625 times.
|
7875 | allocate (this%lin(0:sze - 1)) |
| 220 |
6/10✓ Branch 33 → 34 taken 717 times.
✓ Branch 33 → 35 taken 1908 times.
✓ Branch 35 → 34 taken 1908 times.
✗ Branch 35 → 36 not taken.
✗ Branch 36 → 37 not taken.
✓ Branch 36 → 38 taken 2625 times.
✗ Branch 38 → 39 not taken.
✓ Branch 38 → 40 taken 2625 times.
✗ Branch 40 → 41 not taken.
✓ Branch 40 → 42 taken 2625 times.
|
5250 | allocate (this%i(0:sze - 1)) |
| 221 |
6/10✓ Branch 42 → 43 taken 717 times.
✓ Branch 42 → 44 taken 1908 times.
✓ Branch 44 → 43 taken 1908 times.
✗ Branch 44 → 45 not taken.
✗ Branch 45 → 46 not taken.
✓ Branch 45 → 47 taken 2625 times.
✗ Branch 47 → 48 not taken.
✓ Branch 47 → 49 taken 2625 times.
✗ Branch 49 → 50 not taken.
✓ Branch 49 → 51 taken 2625 times.
|
5250 | allocate (this%j(0:sze - 1)) |
| 222 |
6/10✓ Branch 51 → 52 taken 717 times.
✓ Branch 51 → 53 taken 1908 times.
✓ Branch 53 → 52 taken 1908 times.
✗ Branch 53 → 54 not taken.
✗ Branch 54 → 55 not taken.
✓ Branch 54 → 56 taken 2625 times.
✗ Branch 56 → 57 not taken.
✓ Branch 56 → 58 taken 2625 times.
✗ Branch 58 → 59 not taken.
✓ Branch 58 → 60 taken 2625 times.
|
5250 | allocate (this%k(0:sze - 1)) |
| 223 | |||
| 224 | row_local = 0 | ||
| 225 |
2/2✓ Branch 60 → 61 taken 385 times.
✓ Branch 60 → 62 taken 2240 times.
|
2625 | if (.not. inds%complement) then |
| 226 |
2/2✓ Branch 63 → 64 taken 6243 times.
✓ Branch 63 → 92 taken 2240 times.
|
8483 | do k = inds%k0, inds%k1 |
| 227 |
2/2✓ Branch 65 → 66 taken 157790 times.
✓ Branch 65 → 71 taken 6243 times.
|
166273 | do j = inds%j0, inds%j1 |
| 228 |
2/2✓ Branch 67 → 68 taken 433712 times.
✓ Branch 67 → 69 taken 157790 times.
|
597745 | do i = inds%i0, inds%i1 |
| 229 | 433712 | row_global = l0 + (i - my_inds%i0) + (j - my_inds%j0) * my_sze_x + (k - my_inds%k0) * my_sze_x * my_sze_y | |
| 230 | |||
| 231 | 433712 | this%lin(row_local) = row_global | |
| 232 | 433712 | this%i(row_local) = i | |
| 233 | 433712 | this%j(row_local) = j | |
| 234 | 433712 | this%k(row_local) = k | |
| 235 | |||
| 236 | 591502 | row_local = row_local + 1 | |
| 237 | end do | ||
| 238 | end do | ||
| 239 | end do | ||
| 240 | else | ||
| 241 |
2/2✓ Branch 73 → 74 taken 385 times.
✓ Branch 73 → 75 taken 2079 times.
|
2464 | do k = my_inds%k0, my_inds%k1 |
| 242 |
2/2✓ Branch 76 → 77 taken 34628 times.
✓ Branch 76 → 90 taken 2079 times.
|
37092 | do j = my_inds%j0, my_inds%j1 |
| 243 |
2/2✓ Branch 78 → 79 taken 575786 times.
✓ Branch 78 → 88 taken 34628 times.
|
612493 | do i = my_inds%i0, my_inds%i1 |
| 244 |
12/12✓ Branch 79 → 80 taken 565507 times.
✓ Branch 79 → 86 taken 10279 times.
✓ Branch 80 → 81 taken 555228 times.
✓ Branch 80 → 86 taken 10279 times.
✓ Branch 81 → 82 taken 539578 times.
✓ Branch 81 → 86 taken 15650 times.
✓ Branch 82 → 83 taken 523928 times.
✓ Branch 82 → 86 taken 15650 times.
✓ Branch 83 → 84 taken 520060 times.
✓ Branch 83 → 86 taken 3868 times.
✓ Branch 84 → 85 taken 497508 times.
✓ Branch 84 → 86 taken 22552 times.
|
575786 | if (k >= inds%k0 .and. k <= inds%k1 .and. j >= inds%j0 .and. j <= inds%j1 .and. i >= inds%i0 .and. i <= inds%i1) cycle |
| 245 | 78278 | row_global = l0 + (i - my_inds%i0) + (j - my_inds%j0) * my_sze_x + (k - my_inds%k0) * my_sze_x * my_sze_y | |
| 246 | |||
| 247 | 78278 | this%lin(row_local) = row_global | |
| 248 | 78278 | this%i(row_local) = i | |
| 249 | 78278 | this%j(row_local) = j | |
| 250 | 78278 | this%k(row_local) = k | |
| 251 | |||
| 252 | 610414 | row_local = row_local + 1 | |
| 253 | end do | ||
| 254 | end do | ||
| 255 | end do | ||
| 256 | end if | ||
| 257 | |||
| 258 | ! Set the corresponding TensorProdIndices | ||
| 259 | 2625 | this%tp_inds = inds | |
| 260 | 2625 | end subroutine init_lists | |
| 261 | |||
| 262 | !> @brief Destroy the TensorProdIndexList | ||
| 263 | !> | ||
| 264 | !> @param this[inout] The TensorProdIndexList to destroy | ||
| 265 | 1595 | subroutine destroy_indices(this) | |
| 266 | implicit none | ||
| 267 | |||
| 268 | class(TensorProdIndexList), intent(inout) :: this | ||
| 269 | |||
| 270 |
1/2✓ Branch 2 → 3 taken 1595 times.
✗ Branch 2 → 4 not taken.
|
1595 | if (allocated(this%lin)) deallocate (this%lin) |
| 271 |
1/2✓ Branch 4 → 5 taken 1595 times.
✗ Branch 4 → 6 not taken.
|
1595 | if (allocated(this%i)) deallocate (this%i) |
| 272 |
1/2✓ Branch 6 → 7 taken 1595 times.
✗ Branch 6 → 8 not taken.
|
1595 | if (allocated(this%j)) deallocate (this%j) |
| 273 |
1/2✓ Branch 8 → 9 taken 1595 times.
✗ Branch 8 → 10 not taken.
|
1595 | if (allocated(this%k)) deallocate (this%k) |
| 274 | |||
| 275 | 1595 | end subroutine destroy_indices | |
| 276 | |||
| 277 | !> @brief Get the size of the TensorProdIndices | ||
| 278 | !> | ||
| 279 | !> @param a The TensorProdIndices for which to get the size | ||
| 280 | !> @param dir _(optional)_ The direction for which to get the size (1 for x, 2 for y, 3 for z). If not provided, the total size is | ||
| 281 | !> returned. | ||
| 282 | 716039 | pure integer function size_tpinds(a, dir) | |
| 283 | implicit none | ||
| 284 | |||
| 285 | type(TensorProdIndices), intent(in) :: a | ||
| 286 | integer, intent(in), optional :: dir | ||
| 287 | |||
| 288 |
2/2✓ Branch 2 → 3 taken 712202 times.
✓ Branch 2 → 4 taken 3837 times.
|
716039 | if (.not. present(dir)) then |
| 289 | 712202 | size_tpinds = (a%i1 - a%i0 + 1) * (a%j1 - a%j0 + 1) * (a%k1 - a%k0 + 1) | |
| 290 | else | ||
| 291 | 4413 | select case (dir) | |
| 292 | case (1) ! x-direction | ||
| 293 | 576 | size_tpinds = a%i1 - a%i0 + 1 | |
| 294 | case (2) ! y-direction | ||
| 295 | 576 | size_tpinds = a%j1 - a%j0 + 1 | |
| 296 | case (3) ! z-direction | ||
| 297 | 2685 | size_tpinds = a%k1 - a%k0 + 1 | |
| 298 | case default | ||
| 299 |
3/4✓ Branch 4 → 5 taken 576 times.
✓ Branch 4 → 6 taken 576 times.
✓ Branch 4 → 7 taken 2685 times.
✗ Branch 4 → 8 not taken.
|
3837 | error stop "Invalid direction for TensorProdIndices::size_tpinds" |
| 300 | end select | ||
| 301 | end if | ||
| 302 | |||
| 303 |
1/2✗ Branch 9 → 10 not taken.
✓ Branch 9 → 11 taken 716039 times.
|
716039 | if (a%complement) then |
| 304 | ✗ | size_tpinds = a%total_nr_x * a%total_nr_y * a%total_nr_z - size_tpinds | |
| 305 | end if | ||
| 306 | 716039 | end function size_tpinds | |
| 307 | |||
| 308 | !> @brief Get the size of the TensorProdIndices in 2D | ||
| 309 | !> | ||
| 310 | !> @param a The TensorProdIndices for which to get the size | ||
| 311 | !> @param dir _(optional)_ The direction for which to get the size (1 for x, 2 for y). If not provided, the total size is returned. | ||
| 312 | 1662 | pure integer function size_tp_inds_2d(a, dir) | |
| 313 | implicit none | ||
| 314 | |||
| 315 | type(TensorProdIndices), intent(in) :: a | ||
| 316 | integer, intent(in), optional :: dir | ||
| 317 | |||
| 318 |
1/2✓ Branch 2 → 3 taken 1662 times.
✗ Branch 2 → 4 not taken.
|
1662 | if (.not. present(dir)) then |
| 319 | 1662 | size_tp_inds_2d = (a%i1 - a%i0 + 1) * (a%j1 - a%j0 + 1) | |
| 320 | else | ||
| 321 | ✗ | select case (dir) | |
| 322 | case (1) ! x-direction | ||
| 323 | ✗ | size_tp_inds_2d = a%i1 - a%i0 + 1 | |
| 324 | case (2) ! y-direction | ||
| 325 | ✗ | size_tp_inds_2d = a%j1 - a%j0 + 1 | |
| 326 | case default | ||
| 327 | ✗ | error stop "Invalid direction for TensorProdIndices::size_tp_inds_2d" | |
| 328 | end select | ||
| 329 | end if | ||
| 330 | |||
| 331 | 1662 | end function size_tp_inds_2d | |
| 332 | |||
| 333 | !> @brief Get the size of the TensorProdIndexList | ||
| 334 | !> | ||
| 335 | !> @param a The TensorProdIndexList for which to get the size | ||
| 336 | !> @param dir _(optional)_ The direction for which to get the size (1 for x, 2 for y, 3 for z). If not provided, the total size is | ||
| 337 | !> returned. | ||
| 338 | 7838 | pure integer function size_index_list(a, dir) | |
| 339 | implicit none | ||
| 340 | type(TensorProdIndexList), intent(in) :: a | ||
| 341 | integer, intent(in), optional :: dir | ||
| 342 | 7838 | size_index_list = size(a%tp_inds, dir) | |
| 343 | 7838 | end function size_index_list | |
| 344 | |||
| 345 | !> @brief Get the size of the TensorProdIndexList in 2D | ||
| 346 | !> | ||
| 347 | !> @param a The TensorProdIndexList for which to get the size | ||
| 348 | !> @param dir _(optional)_ The direction for which to get the size (1 for x, 2 for y). If not provided, the total size is returned. | ||
| 349 | 1662 | pure integer function size_index_list_2d(a, dir) | |
| 350 | implicit none | ||
| 351 | type(TensorProdIndexList), intent(in) :: a | ||
| 352 | integer, intent(in), optional :: dir | ||
| 353 | |||
| 354 |
1/2✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 1662 times.
|
1662 | if (a%tp_inds%complement) then |
| 355 | ✗ | error stop "TensorProdIndexList::size_index_list_2d not implemented for complement indices" | |
| 356 | end if | ||
| 357 | |||
| 358 | 1662 | size_index_list_2d = size_tp_inds_2d(a%tp_inds, dir) | |
| 359 | 1662 | end function size_index_list_2d | |
| 360 | |||
| 361 | !> @brief Compute the intersection of two TensorProdIndices | ||
| 362 | !> | ||
| 363 | !> @param a The first TensorProdIndices | ||
| 364 | !> @param b The second TensorProdIndices | ||
| 365 | !> | ||
| 366 | !> @return The intersection of the two TensorProdIndices | ||
| 367 | 385 | pure type(TensorProdIndices) function intersect(a, b) | |
| 368 | implicit none | ||
| 369 | |||
| 370 | type(TensorProdIndices), intent(in) :: a, b | ||
| 371 | |||
| 372 | ! Compute the intersection of two TensorProdIndices | ||
| 373 | |||
| 374 |
2/4✓ Branch 2 → 3 taken 385 times.
✗ Branch 2 → 5 not taken.
✓ Branch 3 → 4 taken 385 times.
✗ Branch 3 → 5 not taken.
|
385 | if (.not. a%complement .and. .not. b%complement) then |
| 375 | 385 | intersect%i0 = max(a%i0, b%i0) | |
| 376 | 385 | intersect%i1 = min(a%i1, b%i1) | |
| 377 | 385 | intersect%j0 = max(a%j0, b%j0) | |
| 378 | 385 | intersect%j1 = min(a%j1, b%j1) | |
| 379 | 385 | intersect%k0 = max(a%k0, b%k0) | |
| 380 | 385 | intersect%k1 = min(a%k1, b%k1) | |
| 381 | |||
| 382 | intersect%complement = .false. | ||
| 383 | |||
| 384 | 385 | intersect%total_nr_x = a%total_nr_x | |
| 385 | 385 | intersect%total_nr_y = a%total_nr_y | |
| 386 | 385 | intersect%total_nr_z = a%total_nr_z | |
| 387 | else | ||
| 388 | ✗ | error stop "Intersecting TensorProdIndices with complement is not implemented" | |
| 389 | end if | ||
| 390 | |||
| 391 | 385 | end function intersect | |
| 392 | |||
| 393 | !> @brief Check if a point is inside the box | ||
| 394 | !> | ||
| 395 | !> @param this The box | ||
| 396 | !> @param x The x-coordinate of the point | ||
| 397 | !> @param y The y-coordinate of the point | ||
| 398 | !> @param z The z-coordinate of the point | ||
| 399 | !> | ||
| 400 | !> @return True if the point is inside the box, false otherwise | ||
| 401 | ✗ | pure logical function is_inside(this, x, y, z) | |
| 402 | implicit none | ||
| 403 | |||
| 404 | class(Box), intent(in) :: this | ||
| 405 | real(wp), intent(in) :: x, y, z | ||
| 406 | |||
| 407 | ! TODO: check what happens on the boundary, and periodicity? | ||
| 408 | is_inside = (x >= this%x0 .and. x <= this%x1 .and. & | ||
| 409 | y >= this%y0 .and. y <= this%y1 .and. & | ||
| 410 | ✗ | z >= this%z0 .and. z <= this%z1) | |
| 411 | |||
| 412 | ✗ | end function is_inside | |
| 413 | |||
| 414 | !> @brief Get the local index in the TensorProdIndexList for a given global Cartesian index | ||
| 415 | !> | ||
| 416 | !> @param this The TensorProdIndexList | ||
| 417 | !> @param i The global x-index | ||
| 418 | !> @param j The global y-index | ||
| 419 | !> @param k The global z-index | ||
| 420 | !> | ||
| 421 | !> @return The local index in the TensorProdIndexList | ||
| 422 | 379024 | pure integer function get_local_index_3d(this, i, j, k) | |
| 423 | implicit none | ||
| 424 | |||
| 425 | class(TensorProdIndexList), intent(in) :: this | ||
| 426 | integer, intent(in) :: i, j, k | ||
| 427 | |||
| 428 |
1/2✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 379024 times.
|
379024 | if (this%tp_inds%complement) then |
| 429 | ✗ | error stop "TensorProdIndexList::get_local_index_3d not implemented for complement indices" | |
| 430 | end if | ||
| 431 | |||
| 432 | get_local_index_3d = (i - this%tp_inds%i0) + (j - this%tp_inds%j0) * (1 + this%tp_inds%i1 - this%tp_inds%i0) + & | ||
| 433 | 379024 | (k - this%tp_inds%k0) * (1 + this%tp_inds%i1 - this%tp_inds%i0) * (1 + this%tp_inds%j1 - this%tp_inds%j0) | |
| 434 | |||
| 435 | 379024 | end function get_local_index_3d | |
| 436 | |||
| 437 | !> @brief Get the local index in the TensorProdIndexList for a given global Cartesian index (2D version) | ||
| 438 | !> | ||
| 439 | !> @param this The TensorProdIndexList | ||
| 440 | !> @param i The global x-index | ||
| 441 | !> @param j The global y-index | ||
| 442 | !> | ||
| 443 | !> @return The local index in the TensorProdIndexList | ||
| 444 | ✗ | pure integer function get_local_index_2d(this, i, j) | |
| 445 | implicit none | ||
| 446 | class(TensorProdIndexList), intent(in) :: this | ||
| 447 | integer, intent(in) :: i, j | ||
| 448 | ✗ | if (this%tp_inds%complement) then | |
| 449 | ✗ | error stop "TensorProdIndexList::get_local_index_2d not implemented for complement indices" | |
| 450 | end if | ||
| 451 | ✗ | get_local_index_2d = (i - this%tp_inds%i0) + (j - this%tp_inds%j0) * (1 + this%tp_inds%i1 - this%tp_inds%i0) | |
| 452 | ✗ | end function get_local_index_2d | |
| 453 | |||
| 454 | !> @brief Get the index bounds from a single TensorProdIndexList | ||
| 455 | !> | ||
| 456 | !> @param[in] index_list The TensorProdIndexList | ||
| 457 | !> @param[out] i_min The minimum i-index | ||
| 458 | !> @param[out] i_max The maximum i-index | ||
| 459 | !> @param[out] j_min The minimum j-index | ||
| 460 | !> @param[out] j_max The maximum j-index | ||
| 461 | !> @param[out] _(optional)_ k_min The minimum k-index | ||
| 462 | !> @param[out] _(optional)_ k_max The maximum k-index | ||
| 463 | 2603 | pure subroutine get_bounds_list(this, i_min, i_max, j_min, j_max, k_min, k_max) | |
| 464 | implicit none | ||
| 465 | |||
| 466 | class(TensorProdIndexList), intent(in) :: this | ||
| 467 | integer, intent(out) :: i_min, i_max, j_min, j_max | ||
| 468 | integer, optional, intent(out) :: k_min, k_max | ||
| 469 | |||
| 470 | 2603 | i_min = this%tp_inds%i0 | |
| 471 | 2603 | i_max = this%tp_inds%i1 | |
| 472 | 2603 | j_min = this%tp_inds%j0 | |
| 473 | 2603 | j_max = this%tp_inds%j1 | |
| 474 |
2/2✓ Branch 2 → 3 taken 2175 times.
✓ Branch 2 → 4 taken 428 times.
|
2603 | if (present(k_min)) k_min = this%tp_inds%k0 |
| 475 |
2/2✓ Branch 4 → 5 taken 2175 times.
✓ Branch 4 → 6 taken 428 times.
|
2603 | if (present(k_max)) k_max = this%tp_inds%k1 |
| 476 | 2603 | end subroutine get_bounds_list | |
| 477 | ✗ | end module m_tensorprod_indices | |
| 478 |