tensorprod/m_tensorprod_quadrature.f90
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | !> @brief Module for precomputed quadrature for tensor product B-splines | ||
| 2 | !> | ||
| 3 | !> A TensorProdQuadrature object is provided wherein the values of the B-splines at the quadrature nodes are precomputed and stored | ||
| 4 | module m_tensorprod_quadrature | ||
| 5 | use m_common, only: wp | ||
| 6 | use m_bspline, only: BSplineSpace | ||
| 7 | use m_bspline_quadrature, only: BSplineQuadrature, DEFAULT_EXTRA_QUADRATURE_POINTS | ||
| 8 | |||
| 9 | implicit none | ||
| 10 | |||
| 11 | private | ||
| 12 | |||
| 13 | public :: TensorProdQuadrature, quadrature, quadrature_product, evaluate_at_nodes, DEFAULT_EXTRA_QUADRATURE_POINTS | ||
| 14 | |||
| 15 | !> @brief A tensor product of BSplineQuadrature objects | ||
| 16 | type TensorProdQuadrature | ||
| 17 | !> The BSplineQuadrature objects for each direction | ||
| 18 | type(BSplineQuadrature) :: bspline_quad(3) | ||
| 19 | contains | ||
| 20 | procedure :: destroy => destroy_tensorprod_quadrature | ||
| 21 | end type TensorProdQuadrature | ||
| 22 | |||
| 23 | !> @brief Construct a TensorProdQuadrature object | ||
| 24 | interface TensorProdQuadrature | ||
| 25 | module procedure init_bspline_quadrature_from_quads, init_bspline_quadrature_from_tp | ||
| 26 | end interface TensorProdQuadrature | ||
| 27 | |||
| 28 | !> @brief Evaluate the quadrature for a tensor product B-spline space by making use of the precomputed quadrature | ||
| 29 | interface quadrature | ||
| 30 | procedure quadrature_nof, quadrature_3x1d, quadrature_3d, quadrature_3d_interval | ||
| 31 | end interface quadrature | ||
| 32 | |||
| 33 | !> @brief Evaluate the quadrature for a product tensor product B-spline spaces by making use of the precomputed quadratures | ||
| 34 | interface quadrature_product | ||
| 35 | procedure quadrature_nof_product, quadrature_3x1d_product, quadrature_3d_product, quadrature_3d_product_interval | ||
| 36 | end interface quadrature_product | ||
| 37 | |||
| 38 | ! TODO much of the functionality of this module is in fact unused, and only used in testing; maybe remove it or move it | ||
| 39 | ! TODO specifically to the testing module? | ||
| 40 | contains | ||
| 41 | |||
| 42 | !> @brief Initialize a TensorProdQuadrature object from three BSplineQuadrature objects | ||
| 43 | !> | ||
| 44 | !> @param[in] bspline1_quad The BSplineQuadrature object for the x direction | ||
| 45 | !> @param[in] bspline2_quad The BSplineQuadrature object for the y direction | ||
| 46 | !> @param[in] bspline3_quad The BSplineQuadrature object for the z direction | ||
| 47 | !> | ||
| 48 | !> @return The initialized TensorProdQuadrature object | ||
| 49 | 7886 | function init_bspline_quadrature_from_quads(bspline1_quad, bspline2_quad, bspline3_quad) result(bsplinetp_quad) | |
| 50 | implicit none | ||
| 51 | |||
| 52 | type(BSplineQuadrature), intent(in) :: bspline1_quad, bspline2_quad, bspline3_quad | ||
| 53 | |||
| 54 | type(TensorProdQuadrature) :: bsplinetp_quad | ||
| 55 | |||
| 56 | 7886 | bsplinetp_quad%bspline_quad(1) = bspline1_quad | |
| 57 | 7886 | bsplinetp_quad%bspline_quad(2) = bspline2_quad | |
| 58 | 7886 | bsplinetp_quad%bspline_quad(3) = bspline3_quad | |
| 59 |
2/2✓ Branch 3 → 4 taken 23658 times.
✓ Branch 3 → 5 taken 7886 times.
|
31544 | end function init_bspline_quadrature_from_quads |
| 60 | |||
| 61 | !> @brief Initialize a TensorProdQuadrature object from a TensorProdSpace object | ||
| 62 | !> | ||
| 63 | !> @param[in] tp_space The TensorProdSpace object | ||
| 64 | !> @param[in] _(optional)_ n_quad The number of quadrature points to use in each direction (default is sufficient for exact | ||
| 65 | !> integration of the product of the B-spline spaces) | ||
| 66 | !> | ||
| 67 | !> @return The initialized TensorProdQuadrature object | ||
| 68 | 32 | function init_bspline_quadrature_from_tp(tp_space, n_quad) result(bsplinetp_quad) | |
| 69 | use m_bspline_functions, only: bspline_eval | ||
| 70 | use m_tensorprod_basis, only: TensorProdSpace | ||
| 71 | use m_quadrature_data | ||
| 72 | |||
| 73 | implicit none | ||
| 74 | |||
| 75 | type(TensorProdSpace), intent(in) :: tp_space | ||
| 76 | type(TensorProdQuadrature) :: bsplinetp_quad | ||
| 77 | integer, optional, intent(in) :: n_quad | ||
| 78 | |||
| 79 |
3/6✓ Branch 7 → 8 taken 32 times.
✗ Branch 7 → 9 not taken.
✓ Branch 9 → 10 taken 32 times.
✗ Branch 9 → 11 not taken.
✓ Branch 11 → 12 taken 32 times.
✗ Branch 11 → 13 not taken.
|
32 | bsplinetp_quad%bspline_quad(1) = BSplineQuadrature(tp_space%spaces(1), n_quad) |
| 80 |
3/6✓ Branch 15 → 16 taken 32 times.
✗ Branch 15 → 17 not taken.
✓ Branch 17 → 18 taken 32 times.
✗ Branch 17 → 19 not taken.
✓ Branch 19 → 20 taken 32 times.
✗ Branch 19 → 21 not taken.
|
32 | bsplinetp_quad%bspline_quad(2) = BSplineQuadrature(tp_space%spaces(2), n_quad) |
| 81 |
3/6✓ Branch 23 → 24 taken 32 times.
✗ Branch 23 → 25 not taken.
✓ Branch 25 → 26 taken 32 times.
✗ Branch 25 → 27 not taken.
✓ Branch 27 → 28 taken 32 times.
✗ Branch 27 → 29 not taken.
|
32 | bsplinetp_quad%bspline_quad(3) = BSplineQuadrature(tp_space%spaces(3), n_quad) |
| 82 |
2/2✓ Branch 3 → 4 taken 96 times.
✓ Branch 3 → 5 taken 32 times.
|
128 | end function init_bspline_quadrature_from_tp |
| 83 | |||
| 84 | !> @brief Destroy a TensorProdQuadrature object | ||
| 85 | !> | ||
| 86 | !> @param[inout] this The TensorProdQuadrature object to destroy | ||
| 87 | 6824 | subroutine destroy_tensorprod_quadrature(this) | |
| 88 | implicit none | ||
| 89 | |||
| 90 | class(TensorProdQuadrature), intent(inout) :: this | ||
| 91 | |||
| 92 | 6824 | call this%bspline_quad(1)%destroy() | |
| 93 | 6824 | call this%bspline_quad(2)%destroy() | |
| 94 | 6824 | call this%bspline_quad(3)%destroy() | |
| 95 | 6824 | end subroutine destroy_tensorprod_quadrature | |
| 96 | |||
| 97 | !> @brief Integrate the (i,j,k)-th B-spline from tensor product B-spline space by making use of the precomputed quadrature | ||
| 98 | !> | ||
| 99 | !> @param[in] this The TensorProdQuadrature object | ||
| 100 | !> @param[in] i The index of the B-spline in the x direction | ||
| 101 | !> @param[in] j The index of the B-spline in the y direction | ||
| 102 | !> @param[in] k The index of the B-spline in the z direction | ||
| 103 | !> | ||
| 104 | !> @return The integral | ||
| 105 | 24653 | pure real(wp) function quadrature_nof(this, i, j, k) result(ans) | |
| 106 | use m_bspline_quadrature, only: quadrature | ||
| 107 | |||
| 108 | implicit none | ||
| 109 | |||
| 110 | type(TensorProdQuadrature), intent(in) :: this | ||
| 111 | integer, intent(in) :: i, j, k | ||
| 112 | |||
| 113 | ans = quadrature(this%bspline_quad(1), i) * & | ||
| 114 | & quadrature(this%bspline_quad(2), j) * & | ||
| 115 | 24653 | & quadrature(this%bspline_quad(3), k) | |
| 116 | 24653 | end function quadrature_nof | |
| 117 | |||
| 118 | !> @brief Integrate the product of the (i,j,k)-th B-spline from tensor product B-spline space with a function | ||
| 119 | !> f(x,y,z)=fx(x)*fy(y)*fz(z) by making use of the precomputed quadrature | ||
| 120 | !> | ||
| 121 | !> @param[in] this The TensorProdQuadrature object | ||
| 122 | !> @param[in] i The index of the B-spline in the x direction | ||
| 123 | !> @param[in] j The index of the B-spline in the y direction | ||
| 124 | !> @param[in] k The index of the B-spline in the z direction | ||
| 125 | !> @param[in] fx The function to integrate in the x direction | ||
| 126 | !> @param[in] fy The function to integrate in the y direction | ||
| 127 | !> @param[in] fz The function to integrate in the z direction | ||
| 128 | !> | ||
| 129 | !> @return The integral | ||
| 130 | 24653 | pure real(wp) function quadrature_3x1d(this, i, j, k, fx, fy, fz) result(ans) | |
| 131 | use m_bspline_quadrature, only: quadrature | ||
| 132 | use m_common, only: user_function_1d_interface | ||
| 133 | |||
| 134 | implicit none | ||
| 135 | |||
| 136 | type(TensorProdQuadrature), intent(in) :: this | ||
| 137 | integer, intent(in) :: i, j, k | ||
| 138 | procedure(user_function_1d_interface) :: fx, fy, fz | ||
| 139 | |||
| 140 | ans = quadrature(this%bspline_quad(1), i, userfun=fx) * & | ||
| 141 | & quadrature(this%bspline_quad(2), j, userfun=fy) * & | ||
| 142 | 24653 | & quadrature(this%bspline_quad(3), k, userfun=fz) | |
| 143 | 24653 | end function quadrature_3x1d | |
| 144 | |||
| 145 | !> @brief Integrate the product of the (i,j,k)-th B-spline from tensor product B-spline space with a function f(x,y,z) by making | ||
| 146 | !> use of the precomputed quadrature | ||
| 147 | !> | ||
| 148 | !> @param[in] this The TensorProdQuadrature object | ||
| 149 | !> @param[in] i The index of the B-spline in the x direction | ||
| 150 | !> @param[in] j The index of the B-spline in the y direction | ||
| 151 | !> @param[in] k The index of the B-spline in the z direction | ||
| 152 | !> @param[in] f The function f(x,y,z) to integrate | ||
| 153 | !> | ||
| 154 | !> @return The integral | ||
| 155 | 73959 | pure real(wp) function quadrature_3d(this, i, j, k, f) result(ans) | |
| 156 | use m_bspline, only: get_index_ranges | ||
| 157 | use m_common, only: user_function_3d_interface | ||
| 158 | |||
| 159 | implicit none | ||
| 160 | |||
| 161 | type(TensorProdQuadrature), intent(in) :: this | ||
| 162 | integer, intent(in) :: i, j, k | ||
| 163 | procedure(user_function_3d_interface) :: f | ||
| 164 | |||
| 165 | integer :: ii, jj, kk ! Interval indices | ||
| 166 | integer :: ii_min0, ii_max0, jj_min0, jj_max0, kk_min0, kk_max0 | ||
| 167 | integer :: ii_min1, ii_max1, jj_min1, jj_max1, kk_min1, kk_max1 | ||
| 168 | |||
| 169 | ! In each direction we get up to two index ranges due to the presence of the boundary | ||
| 170 | 24653 | call get_index_ranges(this%bspline_quad(1)%bspline, i, ii_min0, ii_max0, ii_min1, ii_max1) | |
| 171 | 24653 | call get_index_ranges(this%bspline_quad(2)%bspline, j, jj_min0, jj_max0, jj_min1, jj_max1) | |
| 172 | 24653 | call get_index_ranges(this%bspline_quad(3)%bspline, k, kk_min0, kk_max0, kk_min1, kk_max1) | |
| 173 | |||
| 174 | ! Hence, in 3D, we get up to 8 contributions to the integral | ||
| 175 | |||
| 176 | ans = 0._wp | ||
| 177 |
2/2✓ Branch 6 → 7 taken 106362 times.
✓ Branch 6 → 32 taken 24653 times.
|
131015 | do ii = ii_min0, ii_max0 |
| 178 |
2/2✓ Branch 8 → 9 taken 518190 times.
✓ Branch 8 → 18 taken 106362 times.
|
624552 | do jj = jj_min0, jj_max0 |
| 179 |
2/2✓ Branch 10 → 11 taken 2644208 times.
✓ Branch 10 → 12 taken 518190 times.
|
3162398 | do kk = kk_min0, kk_max0 |
| 180 | 3162398 | ans = ans + quadrature_3d_interval(this, i, j, k, ii, jj, kk, f) | |
| 181 | end do | ||
| 182 |
2/2✓ Branch 14 → 15 taken 617418 times.
✓ Branch 14 → 16 taken 518190 times.
|
1241970 | do kk = kk_min1, kk_max1 |
| 183 | 1135608 | ans = ans + quadrature_3d_interval(this, i, j, k, ii, jj, kk, f) | |
| 184 | end do | ||
| 185 | end do | ||
| 186 |
2/2✓ Branch 20 → 21 taken 105462 times.
✓ Branch 20 → 30 taken 106362 times.
|
236477 | do jj = jj_min1, jj_max1 |
| 187 |
2/2✓ Branch 22 → 23 taken 617418 times.
✓ Branch 22 → 24 taken 105462 times.
|
722880 | do kk = kk_min0, kk_max0 |
| 188 | 722880 | ans = ans + quadrature_3d_interval(this, i, j, k, ii, jj, kk, f) | |
| 189 | end do | ||
| 190 |
2/2✓ Branch 26 → 27 taken 227700 times.
✓ Branch 26 → 28 taken 105462 times.
|
439524 | do kk = kk_min1, kk_max1 |
| 191 | 333162 | ans = ans + quadrature_3d_interval(this, i, j, k, ii, jj, kk, f) | |
| 192 | end do | ||
| 193 | end do | ||
| 194 | end do | ||
| 195 |
2/2✓ Branch 34 → 35 taken 7770 times.
✓ Branch 34 → 60 taken 24653 times.
|
32423 | do ii = ii_min1, ii_max1 |
| 196 |
2/2✓ Branch 36 → 37 taken 44646 times.
✓ Branch 36 → 46 taken 7770 times.
|
52416 | do jj = jj_min0, jj_max0 |
| 197 |
2/2✓ Branch 38 → 39 taken 260742 times.
✓ Branch 38 → 40 taken 44646 times.
|
305388 | do kk = kk_min0, kk_max0 |
| 198 | 305388 | ans = ans + quadrature_3d_interval(this, i, j, k, ii, jj, kk, f) | |
| 199 | end do | ||
| 200 |
2/2✓ Branch 42 → 43 taken 95934 times.
✓ Branch 42 → 44 taken 44646 times.
|
148350 | do kk = kk_min1, kk_max1 |
| 201 | 140580 | ans = ans + quadrature_3d_interval(this, i, j, k, ii, jj, kk, f) | |
| 202 | end do | ||
| 203 | end do | ||
| 204 |
2/2✓ Branch 48 → 49 taken 16170 times.
✓ Branch 48 → 58 taken 7770 times.
|
48593 | do jj = jj_min1, jj_max1 |
| 205 |
2/2✓ Branch 50 → 51 taken 95934 times.
✓ Branch 50 → 52 taken 16170 times.
|
112104 | do kk = kk_min0, kk_max0 |
| 206 | 112104 | ans = ans + quadrature_3d_interval(this, i, j, k, ii, jj, kk, f) | |
| 207 | end do | ||
| 208 |
2/2✓ Branch 54 → 55 taken 35832 times.
✓ Branch 54 → 56 taken 16170 times.
|
59772 | do kk = kk_min1, kk_max1 |
| 209 | 52002 | ans = ans + quadrature_3d_interval(this, i, j, k, ii, jj, kk, f) | |
| 210 | end do | ||
| 211 | end do | ||
| 212 | end do | ||
| 213 | 24653 | end function quadrature_3d | |
| 214 | |||
| 215 | !> @brief Integrate the product of the (i,j,k)-th B-spline from tensor product B-spline space with a function f(x,y,z) over a | ||
| 216 | !> specific interval (ii,jj,kk) by making use of the precomputed quadrature | ||
| 217 | !> | ||
| 218 | !> @param[in] this The TensorProdQuadrature object | ||
| 219 | !> @param[in] i The index of the B-spline in the x direction | ||
| 220 | !> @param[in] j The index of the B-spline in the y direction | ||
| 221 | !> @param[in] k The index of the B-spline in the z direction | ||
| 222 | !> @param[in] ii The index of the interval in the x direction | ||
| 223 | !> @param[in] jj The index of the interval in the y direction | ||
| 224 | !> @param[in] kk The index of the interval in the z direction | ||
| 225 | !> @param[in] f The function f(x,y,z) to integrate | ||
| 226 | !> | ||
| 227 | !> @return The integral over the specified interval | ||
| 228 | 4595186 | pure real(wp) function quadrature_3d_interval(this, i, j, k, ii, jj, kk, f) result(ans) | |
| 229 | use m_bspline, only: get_j_minus_i_quad | ||
| 230 | use m_bspline_quadrature, only: evaluate | ||
| 231 | use m_common, only: user_function_3d_interface | ||
| 232 | |||
| 233 | implicit none | ||
| 234 | |||
| 235 | type(TensorProdQuadrature), intent(in) :: this | ||
| 236 | integer, intent(in) :: i, j, k, ii, jj, kk | ||
| 237 | procedure(user_function_3d_interface) :: f | ||
| 238 | |||
| 239 | integer :: i_minus_ii, j_minus_jj, k_minus_kk, ndx, ndy, ndz | ||
| 240 | real(wp) :: x, y, z, val | ||
| 241 | |||
| 242 | 4595186 | i_minus_ii = get_j_minus_i_quad(this%bspline_quad(1)%bspline, i, ii) | |
| 243 | 4595186 | j_minus_jj = get_j_minus_i_quad(this%bspline_quad(2)%bspline, j, jj) | |
| 244 | 4595186 | k_minus_kk = get_j_minus_i_quad(this%bspline_quad(3)%bspline, k, kk) | |
| 245 | |||
| 246 | ans = 0._wp | ||
| 247 | if (0 <= i_minus_ii .and. i_minus_ii <= this%bspline_quad(1)%bspline%degree .and. & | ||
| 248 | 0 <= j_minus_jj .and. j_minus_jj <= this%bspline_quad(2)%bspline%degree .and. & | ||
| 249 |
9/12✓ Branch 2 → 3 taken 4433906 times.
✓ Branch 2 → 21 taken 161280 times.
✓ Branch 3 → 4 taken 4433906 times.
✗ Branch 3 → 21 not taken.
✓ Branch 4 → 5 taken 4125710 times.
✓ Branch 4 → 21 taken 308196 times.
✓ Branch 5 → 6 taken 4125710 times.
✗ Branch 5 → 21 not taken.
✓ Branch 6 → 7 taken 3844946 times.
✓ Branch 6 → 21 taken 280764 times.
✓ Branch 7 → 8 taken 3844946 times.
✗ Branch 7 → 21 not taken.
|
4595186 | 0 <= k_minus_kk .and. k_minus_kk <= this%bspline_quad(3)%bspline%degree) then |
| 250 | |||
| 251 |
2/2✓ Branch 9 → 10 taken 17091980 times.
✓ Branch 9 → 20 taken 3844946 times.
|
20936926 | do ndx = 1, this%bspline_quad(1)%n_quad |
| 252 |
2/2✓ Branch 11 → 12 taken 77627064 times.
✓ Branch 11 → 18 taken 17091980 times.
|
98563990 | do ndy = 1, this%bspline_quad(2)%n_quad |
| 253 |
2/2✓ Branch 13 → 14 taken 358736306 times.
✓ Branch 13 → 16 taken 77627064 times.
|
453455350 | do ndz = 1, this%bspline_quad(3)%n_quad |
| 254 | val = evaluate(this%bspline_quad(1), i, i_minus_ii, ndx) * & | ||
| 255 | & evaluate(this%bspline_quad(2), j, j_minus_jj, ndy) * & | ||
| 256 | 358736306 | evaluate(this%bspline_quad(3), k, k_minus_kk, ndz) | |
| 257 | 358736306 | x = (ii + .5_wp + this%bspline_quad(1)%nodes(ndx) / 2) * this%bspline_quad(1)%bspline%h | |
| 258 | 358736306 | y = (jj + .5_wp + this%bspline_quad(2)%nodes(ndy) / 2) * this%bspline_quad(2)%bspline%h | |
| 259 | 358736306 | z = (kk + .5_wp + this%bspline_quad(3)%nodes(ndz) / 2) * this%bspline_quad(3)%bspline%h | |
| 260 | 358736306 | val = val * f(x, y, z) | |
| 261 | 358736306 | val = val * this%bspline_quad(1)%weights(ndx) * this%bspline_quad(2)%weights(ndy) * this%bspline_quad(3)%weights(ndz) | |
| 262 | 436363370 | ans = ans + val | |
| 263 | end do | ||
| 264 | end do | ||
| 265 | end do | ||
| 266 | end if | ||
| 267 | 4595186 | end function quadrature_3d_interval | |
| 268 | |||
| 269 | !> @brief Integrate the product of two tensor product B-splines by making use of the precomputed quadratures | ||
| 270 | !> | ||
| 271 | !> @param[in] this1 The first TensorProdQuadrature object | ||
| 272 | !> @param[in] this2 The second TensorProdQuadrature object | ||
| 273 | !> @param[in] i1 The index of the B-spline in the x direction for the first object | ||
| 274 | !> @param[in] j1 The index of the B-spline in the y direction for the first object | ||
| 275 | !> @param[in] k1 The index of the B-spline in the z direction for the first object | ||
| 276 | !> @param[in] i2 The index of the B-spline in the x direction for the second object | ||
| 277 | !> @param[in] j2 The index of the B-spline in the y direction for the second object | ||
| 278 | !> @param[in] k2 The index of the B-spline in the z direction for the second object | ||
| 279 | !> | ||
| 280 | !> @return The integral of the product of the two B-splines | ||
| 281 | 665445 | pure real(wp) function quadrature_nof_product(this1, this2, i1, j1, k1, i2, j2, k2) result(ans) | |
| 282 | use m_bspline_quadrature, only: quadrature_product | ||
| 283 | |||
| 284 | implicit none | ||
| 285 | |||
| 286 | type(TensorProdQuadrature), intent(in) :: this1, this2 | ||
| 287 | integer, intent(in) :: i1, j1, k1, i2, j2, k2 | ||
| 288 | |||
| 289 | ans = quadrature_product(this1%bspline_quad(1), this2%bspline_quad(1), i1, i2) * & | ||
| 290 | & quadrature_product(this1%bspline_quad(2), this2%bspline_quad(2), j1, j2) * & | ||
| 291 | 665445 | & quadrature_product(this1%bspline_quad(3), this2%bspline_quad(3), k1, k2) | |
| 292 | 665445 | end function quadrature_nof_product | |
| 293 | |||
| 294 | !> @brief Integrate the product of two tensor product B-splines with the function f(x,y,z)=fx(x)*fy(y)*fz(z) by making use of the | ||
| 295 | !> precomputed quadratures | ||
| 296 | !> | ||
| 297 | !> @param[in] this1 The first TensorProdQuadrature object | ||
| 298 | !> @param[in] this2 The second TensorProdQuadrature object | ||
| 299 | !> @param[in] i1 The index of the B-spline in the x direction for the first object | ||
| 300 | !> @param[in] j1 The index of the B-spline in the y direction for the first object | ||
| 301 | !> @param[in] k1 The index of the B-spline in the z direction for the first object | ||
| 302 | !> @param[in] i2 The index of the B-spline in the x direction for the second object | ||
| 303 | !> @param[in] j2 The index of the B-spline in the y direction for the second object | ||
| 304 | !> @param[in] k2 The index of the B-spline in the z direction for the second object | ||
| 305 | !> @param[in] fx The function to integrate in the x direction | ||
| 306 | !> @param[in] fy The function to integrate in the y direction | ||
| 307 | !> @param[in] fz The function to integrate in the z direction | ||
| 308 | !> | ||
| 309 | !> @return The integral of the product of the two B-splines with the function | ||
| 310 | 665445 | pure real(wp) function quadrature_3x1d_product(this1, this2, i1, j1, k1, i2, j2, k2, fx, fy, fz) result(ans) | |
| 311 | use m_bspline_quadrature, only: quadrature_product | ||
| 312 | use m_common, only: user_function_1d_interface | ||
| 313 | |||
| 314 | implicit none | ||
| 315 | |||
| 316 | type(TensorProdQuadrature), intent(in) :: this1, this2 | ||
| 317 | integer, intent(in) :: i1, j1, k1, i2, j2, k2 | ||
| 318 | procedure(user_function_1d_interface) :: fx, fy, fz | ||
| 319 | |||
| 320 | ans = quadrature_product(this1%bspline_quad(1), this2%bspline_quad(1), i1, i2, userfun=fx) * & | ||
| 321 | & quadrature_product(this1%bspline_quad(2), this2%bspline_quad(2), j1, j2, userfun=fy) * & | ||
| 322 | 665445 | & quadrature_product(this1%bspline_quad(3), this2%bspline_quad(3), k1, k2, userfun=fz) | |
| 323 | 665445 | end function quadrature_3x1d_product | |
| 324 | |||
| 325 | !> @brief Integrate the product of two tensor product B-splines with the function f(x,y,z) by making use of the precomputed | ||
| 326 | !> quadratures | ||
| 327 | !> | ||
| 328 | !> @param[in] this1 The first TensorProdQuadrature object | ||
| 329 | !> @param[in] this2 The second TensorProdQuadrature object | ||
| 330 | !> @param[in] i1 The index of the B-spline in the x direction for the first object | ||
| 331 | !> @param[in] j1 The index of the B-spline in the y direction for the first object | ||
| 332 | !> @param[in] k1 The index of the B-spline in the z direction for the first object | ||
| 333 | !> @param[in] i2 The index of the B-spline in the x direction for the second object | ||
| 334 | !> @param[in] j2 The index of the B-spline in the y direction for the second object | ||
| 335 | !> @param[in] k2 The index of the B-spline in the z direction for the second object | ||
| 336 | !> @param[in] f The function f(x,y,z) to integrate | ||
| 337 | !> | ||
| 338 | !> @return The integral of the product of the two B-splines with the function | ||
| 339 | 1996335 | pure real(wp) function quadrature_3d_product(this1, this2, i1, j1, k1, i2, j2, k2, f) result(ans) | |
| 340 | use m_bspline, only: get_index_ranges | ||
| 341 | use m_common, only: user_function_3d_interface | ||
| 342 | |||
| 343 | implicit none | ||
| 344 | |||
| 345 | type(TensorProdQuadrature), intent(in) :: this1, this2 | ||
| 346 | integer, intent(in) :: i1, j1, k1, i2, j2, k2 | ||
| 347 | procedure(user_function_3d_interface) :: f | ||
| 348 | |||
| 349 | integer :: ii, jj, kk ! Interval indices | ||
| 350 | integer :: ii_min0, ii_max0, jj_min0, jj_max0, kk_min0, kk_max0 | ||
| 351 | integer :: ii_min1, ii_max1, jj_min1, jj_max1, kk_min1, kk_max1 | ||
| 352 | |||
| 353 | ! In each direction we get up to two index ranges due to the presence of the boundary | ||
| 354 | 665445 | call get_index_ranges(this1%bspline_quad(1)%bspline, this2%bspline_quad(1)%bspline, i1, i2, ii_min0, ii_max0, ii_min1, ii_max1) | |
| 355 | 665445 | call get_index_ranges(this1%bspline_quad(2)%bspline, this2%bspline_quad(2)%bspline, j1, j2, jj_min0, jj_max0, jj_min1, jj_max1) | |
| 356 | 665445 | call get_index_ranges(this1%bspline_quad(3)%bspline, this2%bspline_quad(3)%bspline, k1, k2, kk_min0, kk_max0, kk_min1, kk_max1) | |
| 357 | |||
| 358 | ! Hence, in 3D, we get up to 8 contributions to the integral | ||
| 359 | |||
| 360 | ans = 0._wp | ||
| 361 |
2/2✓ Branch 6 → 7 taken 912088 times.
✓ Branch 6 → 32 taken 665445 times.
|
1577533 | do ii = ii_min0, ii_max0 |
| 362 |
2/2✓ Branch 8 → 9 taken 1261462 times.
✓ Branch 8 → 18 taken 912088 times.
|
2173550 | do jj = jj_min0, jj_max0 |
| 363 |
2/2✓ Branch 10 → 11 taken 1756862 times.
✓ Branch 10 → 12 taken 1261462 times.
|
3018324 | do kk = kk_min0, kk_max0 |
| 364 | 3018324 | ans = ans + quadrature_3d_product_interval(this1, this2, i1, j1, k1, i2, j2, k2, ii, jj, kk, f) | |
| 365 | end do | ||
| 366 |
1/2✗ Branch 14 → 15 not taken.
✓ Branch 14 → 16 taken 1261462 times.
|
2173550 | do kk = kk_min1, kk_max1 |
| 367 | 1261462 | ans = ans + quadrature_3d_product_interval(this1, this2, i1, j1, k1, i2, j2, k2, ii, jj, kk, f) | |
| 368 | end do | ||
| 369 | end do | ||
| 370 |
1/2✗ Branch 20 → 21 not taken.
✓ Branch 20 → 30 taken 912088 times.
|
1577533 | do jj = jj_min1, jj_max1 |
| 371 | ✗ | do kk = kk_min0, kk_max0 | |
| 372 | ✗ | ans = ans + quadrature_3d_product_interval(this1, this2, i1, j1, k1, i2, j2, k2, ii, jj, kk, f) | |
| 373 | end do | ||
| 374 |
0/2✗ Branch 26 → 27 not taken.
✗ Branch 26 → 28 not taken.
|
912088 | do kk = kk_min1, kk_max1 |
| 375 | ✗ | ans = ans + quadrature_3d_product_interval(this1, this2, i1, j1, k1, i2, j2, k2, ii, jj, kk, f) | |
| 376 | end do | ||
| 377 | end do | ||
| 378 | end do | ||
| 379 |
1/2✗ Branch 34 → 35 not taken.
✓ Branch 34 → 60 taken 665445 times.
|
665445 | do ii = ii_min1, ii_max1 |
| 380 | ✗ | do jj = jj_min0, jj_max0 | |
| 381 | ✗ | do kk = kk_min0, kk_max0 | |
| 382 | ✗ | ans = ans + quadrature_3d_product_interval(this1, this2, i1, j1, k1, i2, j2, k2, ii, jj, kk, f) | |
| 383 | end do | ||
| 384 | ✗ | do kk = kk_min1, kk_max1 | |
| 385 | ✗ | ans = ans + quadrature_3d_product_interval(this1, this2, i1, j1, k1, i2, j2, k2, ii, jj, kk, f) | |
| 386 | end do | ||
| 387 | end do | ||
| 388 |
0/2✗ Branch 48 → 49 not taken.
✗ Branch 48 → 58 not taken.
|
665445 | do jj = jj_min1, jj_max1 |
| 389 | ✗ | do kk = kk_min0, kk_max0 | |
| 390 | ✗ | ans = ans + quadrature_3d_product_interval(this1, this2, i1, j1, k1, i2, j2, k2, ii, jj, kk, f) | |
| 391 | end do | ||
| 392 | ✗ | do kk = kk_min1, kk_max1 | |
| 393 | ✗ | ans = ans + quadrature_3d_product_interval(this1, this2, i1, j1, k1, i2, j2, k2, ii, jj, kk, f) | |
| 394 | end do | ||
| 395 | end do | ||
| 396 | end do | ||
| 397 | 665445 | end function quadrature_3d_product | |
| 398 | |||
| 399 | !> @brief Integrate the product of two tensor product B-splines with the function f(x,y,z) over a specific interval (ii,jj,kk) by | ||
| 400 | !> making use of the precomputed quadratures | ||
| 401 | !> | ||
| 402 | !> @param[in] this1 The first TensorProdQuadrature object | ||
| 403 | !> @param[in] this2 The second TensorProdQuadrature object | ||
| 404 | !> @param[in] i1 The index of the B-spline in the x direction for the first object | ||
| 405 | !> @param[in] j1 The index of the B-spline in the y direction for the first object | ||
| 406 | !> @param[in] k1 The index of the B-spline in the z direction for the first object | ||
| 407 | !> @param[in] i2 The index of the B-spline in the x direction for the second object | ||
| 408 | !> @param[in] j2 The index of the B-spline in the y direction for the second object | ||
| 409 | !> @param[in] k2 The index of the B-spline in the z direction for the second object | ||
| 410 | !> @param[in] ii The index of the interval in the x direction | ||
| 411 | !> @param[in] jj The index of the interval in the y direction | ||
| 412 | !> @param[in] kk The index of the interval in the z direction | ||
| 413 | !> @param[in] f The function f(x,y,z) to integrate | ||
| 414 | !> | ||
| 415 | !> @return The integral of the product of the two B-splines with the function over the specified interval | ||
| 416 | 1756862 | pure real(wp) function quadrature_3d_product_interval(this1, this2, i1, j1, k1, i2, j2, k2, ii, jj, kk, f) result(ans) | |
| 417 | use m_bspline, only: get_j_minus_i_quad | ||
| 418 | use m_bspline_quadrature, only: evaluate | ||
| 419 | use m_common, only: user_function_3d_interface | ||
| 420 | |||
| 421 | implicit none | ||
| 422 | |||
| 423 | type(TensorProdQuadrature), intent(in) :: this1, this2 | ||
| 424 | integer, intent(in) :: i1, j1, k1, i2, j2, k2, ii, jj, kk | ||
| 425 | procedure(user_function_3d_interface) :: f | ||
| 426 | |||
| 427 | integer :: i1_minus_ii, j1_minus_jj, k1_minus_kk, i2_minus_ii, j2_minus_jj, k2_minus_kk, ndx, ndy, ndz | ||
| 428 | real(wp) :: x, y, z, val | ||
| 429 | |||
| 430 | 1756862 | i1_minus_ii = get_j_minus_i_quad(this1%bspline_quad(1)%bspline, i1, ii) | |
| 431 | 1756862 | j1_minus_jj = get_j_minus_i_quad(this1%bspline_quad(2)%bspline, j1, jj) | |
| 432 | 1756862 | k1_minus_kk = get_j_minus_i_quad(this1%bspline_quad(3)%bspline, k1, kk) | |
| 433 | |||
| 434 | 1756862 | i2_minus_ii = get_j_minus_i_quad(this2%bspline_quad(1)%bspline, i2, ii) | |
| 435 | 1756862 | j2_minus_jj = get_j_minus_i_quad(this2%bspline_quad(2)%bspline, j2, jj) | |
| 436 | 1756862 | k2_minus_kk = get_j_minus_i_quad(this2%bspline_quad(3)%bspline, k2, kk) | |
| 437 | |||
| 438 | ans = 0._wp | ||
| 439 | if (0 <= i1_minus_ii .and. i1_minus_ii <= this1%bspline_quad(1)%bspline%degree .and. & | ||
| 440 | 0 <= j1_minus_jj .and. j1_minus_jj <= this1%bspline_quad(2)%bspline%degree .and. & | ||
| 441 | 0 <= k1_minus_kk .and. k1_minus_kk <= this1%bspline_quad(3)%bspline%degree .and. & | ||
| 442 | 0 <= i2_minus_ii .and. i2_minus_ii <= this2%bspline_quad(1)%bspline%degree .and. & | ||
| 443 | 0 <= j2_minus_jj .and. j2_minus_jj <= this2%bspline_quad(2)%bspline%degree .and. & | ||
| 444 |
12/24✓ Branch 2 → 3 taken 1756862 times.
✗ Branch 2 → 27 not taken.
✓ Branch 3 → 4 taken 1756862 times.
✗ Branch 3 → 27 not taken.
✓ Branch 4 → 5 taken 1756862 times.
✗ Branch 4 → 27 not taken.
✓ Branch 5 → 6 taken 1756862 times.
✗ Branch 5 → 27 not taken.
✓ Branch 6 → 7 taken 1756862 times.
✗ Branch 6 → 27 not taken.
✓ Branch 7 → 8 taken 1756862 times.
✗ Branch 7 → 27 not taken.
✓ Branch 8 → 9 taken 1756862 times.
✗ Branch 8 → 27 not taken.
✓ Branch 9 → 10 taken 1756862 times.
✗ Branch 9 → 27 not taken.
✓ Branch 10 → 11 taken 1756862 times.
✗ Branch 10 → 27 not taken.
✓ Branch 11 → 12 taken 1756862 times.
✗ Branch 11 → 27 not taken.
✓ Branch 12 → 13 taken 1756862 times.
✗ Branch 12 → 27 not taken.
✓ Branch 13 → 14 taken 1756862 times.
✗ Branch 13 → 27 not taken.
|
1756862 | 0 <= k2_minus_kk .and. k2_minus_kk <= this2%bspline_quad(3)%bspline%degree) then |
| 445 | |||
| 446 |
2/2✓ Branch 15 → 16 taken 8577446 times.
✓ Branch 15 → 26 taken 1756862 times.
|
10334308 | do ndx = 1, this1%bspline_quad(1)%n_quad |
| 447 |
2/2✓ Branch 17 → 18 taken 42081820 times.
✓ Branch 17 → 24 taken 8577446 times.
|
52416128 | do ndy = 1, this1%bspline_quad(2)%n_quad |
| 448 |
2/2✓ Branch 19 → 20 taken 207252848 times.
✓ Branch 19 → 22 taken 42081820 times.
|
257912114 | do ndz = 1, this1%bspline_quad(3)%n_quad |
| 449 | val = evaluate(this1%bspline_quad(1), i1, i1_minus_ii, ndx) * & | ||
| 450 | & evaluate(this1%bspline_quad(2), j1, j1_minus_jj, ndy) * & | ||
| 451 | & evaluate(this1%bspline_quad(3), k1, k1_minus_kk, ndz) * & | ||
| 452 | & evaluate(this2%bspline_quad(1), i2, i2_minus_ii, ndx) * & | ||
| 453 | & evaluate(this2%bspline_quad(2), j2, j2_minus_jj, ndy) * & | ||
| 454 | 207252848 | evaluate(this2%bspline_quad(3), k2, k2_minus_kk, ndz) | |
| 455 | 207252848 | x = (ii + .5_wp + this1%bspline_quad(1)%nodes(ndx) / 2) * this1%bspline_quad(1)%bspline%h | |
| 456 | 207252848 | y = (jj + .5_wp + this1%bspline_quad(2)%nodes(ndy) / 2) * this1%bspline_quad(2)%bspline%h | |
| 457 | 207252848 | z = (kk + .5_wp + this1%bspline_quad(3)%nodes(ndz) / 2) * this1%bspline_quad(3)%bspline%h | |
| 458 | 207252848 | val = val * f(x, y, z) | |
| 459 | 207252848 | val = val * this1%bspline_quad(1)%weights(ndx) * this1%bspline_quad(2)%weights(ndy) * this1%bspline_quad(3)%weights(ndz) | |
| 460 | 249334668 | ans = ans + val | |
| 461 | end do | ||
| 462 | end do | ||
| 463 | end do | ||
| 464 | end if | ||
| 465 | 1756862 | end function quadrature_3d_product_interval | |
| 466 | |||
| 467 | !> @brief Evaluate the function f at the quadrature nodes on the (ii,jj,kk)-th interval | ||
| 468 | !> | ||
| 469 | !> @param[in] this The TensorProdQuadrature object | ||
| 470 | !> @param[in] f The function to evaluate | ||
| 471 | !> @param[in] ii The index of the interval in the x direction | ||
| 472 | !> @param[in] jj The index of the interval in the y direction | ||
| 473 | !> @param[in] kk The index of the interval in the z direction | ||
| 474 | !> @param[out] f_vals The values of the function at the quadrature nodes, shape (n_quad_z, n_quad_y, n_quad_x) | ||
| 475 | !> | ||
| 476 | !> @note The output array f_vals is intentionally in reverse order, as this facilitates its use with the current loop-order in the | ||
| 477 | !> mass matrix assembly | ||
| 478 |
1/2✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 3433700 times.
|
3433700 | pure subroutine evaluate_at_nodes(this, f, ii, jj, kk, f_vals) |
| 479 | use m_common, only: user_function_3d_interface | ||
| 480 | implicit none | ||
| 481 | |||
| 482 | type(TensorProdQuadrature), intent(in) :: this | ||
| 483 | procedure(user_function_3d_interface) :: f | ||
| 484 | integer, intent(in) :: ii, jj, kk | ||
| 485 | real(wp), intent(out) :: f_vals(1:, 1:, 1:) | ||
| 486 | |||
| 487 | integer :: ndx, ndy, ndz | ||
| 488 | real(wp) :: x, y, z | ||
| 489 | |||
| 490 |
2/2✓ Branch 5 → 6 taken 19829155 times.
✓ Branch 5 → 16 taken 3433700 times.
|
23262855 | do ndx = 1, this%bspline_quad(1)%n_quad |
| 491 |
2/2✓ Branch 7 → 8 taken 117290620 times.
✓ Branch 7 → 14 taken 19829155 times.
|
140553475 | do ndy = 1, this%bspline_quad(2)%n_quad |
| 492 |
2/2✓ Branch 9 → 10 taken 534854538 times.
✓ Branch 9 → 12 taken 117290620 times.
|
671974313 | do ndz = 1, this%bspline_quad(3)%n_quad |
| 493 | 534854538 | x = (ii + .5_wp + this%bspline_quad(1)%nodes(ndx) / 2) * this%bspline_quad(1)%bspline%h | |
| 494 | 534854538 | y = (jj + .5_wp + this%bspline_quad(2)%nodes(ndy) / 2) * this%bspline_quad(2)%bspline%h | |
| 495 | 534854538 | z = (kk + .5_wp + this%bspline_quad(3)%nodes(ndz) / 2) * this%bspline_quad(3)%bspline%h | |
| 496 | 652145158 | f_vals(ndz, ndy, ndx) = f(x, y, z) | |
| 497 | end do | ||
| 498 | end do | ||
| 499 | end do | ||
| 500 | 3433700 | end subroutine evaluate_at_nodes | |
| 501 | |||
| 502 | ✗ | end module m_tensorprod_quadrature | |
| 503 |