bspline/m_bspline_quadrature.F90
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | !> The B-spline quadrature module | ||
| 2 | !> | ||
| 3 | !> This module provides functionality for numerical integration (i.e., quadrature rules) over B-spline spaces. | ||
| 4 | module m_bspline_quadrature | ||
| 5 | use m_common, only: wp | ||
| 6 | use m_bspline_basis, only: BSplineSpace | ||
| 7 | use m_quadrature, only: MAX_QUADRATURE_NODES | ||
| 8 | |||
| 9 | implicit none | ||
| 10 | |||
| 11 | !> The default number of extra quadrature points to use for B-spline quadrature | ||
| 12 | integer, parameter :: DEFAULT_EXTRA_QUADRATURE_POINTS = 2 | ||
| 13 | |||
| 14 | private | ||
| 15 | |||
| 16 | public :: DEFAULT_EXTRA_QUADRATURE_POINTS | ||
| 17 | public :: BSplineQuadrature, quadrature, quadrature_product, evaluate, n_quad_exact | ||
| 18 | |||
| 19 | !> @brief A type to hold the B-spline values at quadrature nodes | ||
| 20 | !> | ||
| 21 | !> This type merely wraps an array such that we can have an array of arrays | ||
| 22 | type QuadArray | ||
| 23 | real(wp) :: data(MAX_QUADRATURE_NODES) | ||
| 24 | end type QuadArray | ||
| 25 | |||
| 26 | !> @brief A type to hold the B-spline quadrature data: the values of all of the B-splines at the quadrature nodes | ||
| 27 | type BSplineQuadrature | ||
| 28 | !> The number of quadrature nodes | ||
| 29 | integer :: n_quad | ||
| 30 | |||
| 31 | !> The B-spline space for which the quadrature is defined | ||
| 32 | type(BSplineSpace) :: bspline | ||
| 33 | |||
| 34 | !> The quadrature weights (scaled by the interval length 'h') | ||
| 35 | real(wp), allocatable :: weights(:) | ||
| 36 | |||
| 37 | !> The quadrature nodes (on the interval (-1, 1)) | ||
| 38 | real(wp), allocatable :: nodes(:) | ||
| 39 | |||
| 40 | !> The B-spline values at the quadrature nodes | ||
| 41 | type(QuadArray), allocatable :: bspline_at_nodes(:, :) ! j, j-1 | ||
| 42 | contains | ||
| 43 | procedure :: destroy => destroy_bspline_quadrature | ||
| 44 | ! final :: destroy_bspline_quadrature_final | ||
| 45 | procedure, private, pass(to) :: copy => copy_bspline_quadrature | ||
| 46 | generic, public :: assignment(=) => copy | ||
| 47 | end type BSplineQuadrature | ||
| 48 | |||
| 49 | !> @brief Create a B-spline quadrature object | ||
| 50 | interface BSplineQuadrature | ||
| 51 | module procedure init_bspline_quadrature | ||
| 52 | end interface BSplineQuadrature | ||
| 53 | |||
| 54 | !> @brief Quadrature rules on a B-spline space | ||
| 55 | interface quadrature | ||
| 56 | procedure quadrature_1d, quadrature_1d_interval | ||
| 57 | end interface quadrature | ||
| 58 | |||
| 59 | !> @brief Quadrature rules on a product of two B-spline spaces | ||
| 60 | interface quadrature_product | ||
| 61 | procedure quadrature_product_1d, quadrature_product_1d_interval | ||
| 62 | end interface quadrature_product | ||
| 63 | |||
| 64 | !> @brief Determine the number of quadrature points needed for exact integration of a (product of) B-spline space(s) | ||
| 65 | interface n_quad_exact | ||
| 66 | procedure n_quad_exact_single, n_quad_exact_double | ||
| 67 | end interface n_quad_exact | ||
| 68 | contains | ||
| 69 | |||
| 70 | !> @brief Initialize a B-spline quadrature object ! TODO use TBP init | ||
| 71 | !> | ||
| 72 | !> @param[in] bspline The B-spline space for which the quadrature is defined | ||
| 73 | !> @param[in] n_quad _(optional)_ The number of quadrature points to use (default is the number needed for exact integration of | ||
| 74 | !> the product of the B-spline space with itself) | ||
| 75 | !> | ||
| 76 | !> @return The B-spline quadrature object | ||
| 77 | 35406 | function init_bspline_quadrature(bspline, n_quad) result(bspline_quadrature) | |
| 78 | use m_bspline_functions, only: bspline_eval | ||
| 79 | use m_quadrature_data | ||
| 80 | |||
| 81 | type(BSplineSpace), intent(in) :: bspline | ||
| 82 | integer, optional, intent(in) :: n_quad | ||
| 83 | |||
| 84 | type(BSplineQuadrature) :: bspline_quadrature | ||
| 85 | |||
| 86 | integer :: j, j_minus_i, ndx, j0, j1, n_quad_ | ||
| 87 | real(wp) :: x | ||
| 88 | |||
| 89 |
1/2✓ Branch 2 → 3 taken 35406 times.
✗ Branch 2 → 4 not taken.
|
35406 | if (present(n_quad)) then |
| 90 | 35406 | n_quad_ = n_quad | |
| 91 | else | ||
| 92 | ! By default we ensure that we can integrate products of the B-spline space exactly | ||
| 93 | n_quad_ = n_quad_exact_double(bspline, bspline) | ||
| 94 | end if | ||
| 95 | 35406 | n_quad_ = max(n_quad_, 1) ! Support 2D problems having degree = -1 in one direction | |
| 96 | |||
| 97 | 35406 | bspline_quadrature%bspline = bspline | |
| 98 | bspline_quadrature%n_quad = n_quad_ | ||
| 99 | |||
| 100 | if (allocated(bspline_quadrature%weights)) deallocate (bspline_quadrature%weights) | ||
| 101 |
5/10✗ Branch 5 → 6 not taken.
✓ Branch 5 → 7 taken 35406 times.
✓ Branch 7 → 6 taken 35406 times.
✗ Branch 7 → 8 not taken.
✓ Branch 8 → 9 taken 35406 times.
✗ Branch 8 → 10 not taken.
✗ Branch 10 → 11 not taken.
✓ Branch 10 → 12 taken 35406 times.
✗ Branch 12 → 13 not taken.
✓ Branch 12 → 14 taken 35406 times.
|
106218 | allocate (bspline_quadrature%weights(n_quad_)) |
| 102 | if (allocated(bspline_quadrature%nodes)) deallocate (bspline_quadrature%nodes) | ||
| 103 |
4/8✗ Branch 14 → 15 not taken.
✓ Branch 14 → 16 taken 35406 times.
✓ Branch 16 → 15 taken 35406 times.
✗ Branch 16 → 17 not taken.
✗ Branch 17 → 18 not taken.
✓ Branch 17 → 19 taken 35406 times.
✗ Branch 19 → 20 not taken.
✓ Branch 19 → 21 taken 35406 times.
|
70812 | allocate (bspline_quadrature%nodes(n_quad_)) |
| 104 | |||
| 105 | 2611 | select case (n_quad_) | |
| 106 | case (1) | ||
| 107 |
3/6✗ Branch 22 → 23 not taken.
✓ Branch 22 → 26 taken 2611 times.
✗ Branch 23 → 24 not taken.
✗ Branch 23 → 25 not taken.
✓ Branch 27 → 28 taken 2611 times.
✓ Branch 27 → 29 taken 2611 times.
|
7833 | bspline_quadrature%weights = gl_weights_1 * bspline%h / 2 |
| 108 |
3/6✗ Branch 29 → 30 not taken.
✓ Branch 29 → 33 taken 2611 times.
✗ Branch 30 → 31 not taken.
✗ Branch 30 → 32 not taken.
✓ Branch 34 → 35 taken 2611 times.
✓ Branch 34 → 247 taken 2611 times.
|
7833 | bspline_quadrature%nodes = gl_nodes_1 |
| 109 | case (2) | ||
| 110 |
3/6✗ Branch 36 → 37 not taken.
✓ Branch 36 → 40 taken 2541 times.
✗ Branch 37 → 38 not taken.
✗ Branch 37 → 39 not taken.
✓ Branch 41 → 42 taken 5082 times.
✓ Branch 41 → 43 taken 2541 times.
|
10164 | bspline_quadrature%weights = gl_weights_2 * bspline%h / 2 |
| 111 |
3/6✗ Branch 43 → 44 not taken.
✓ Branch 43 → 47 taken 2541 times.
✗ Branch 44 → 45 not taken.
✗ Branch 44 → 46 not taken.
✓ Branch 48 → 49 taken 5082 times.
✓ Branch 48 → 247 taken 2541 times.
|
10164 | bspline_quadrature%nodes = gl_nodes_2 |
| 112 | case (3) | ||
| 113 |
3/6✗ Branch 50 → 51 not taken.
✓ Branch 50 → 54 taken 5556 times.
✗ Branch 51 → 52 not taken.
✗ Branch 51 → 53 not taken.
✓ Branch 55 → 56 taken 16668 times.
✓ Branch 55 → 57 taken 5556 times.
|
27780 | bspline_quadrature%weights = gl_weights_3 * bspline%h / 2 |
| 114 |
3/6✗ Branch 57 → 58 not taken.
✓ Branch 57 → 61 taken 5556 times.
✗ Branch 58 → 59 not taken.
✗ Branch 58 → 60 not taken.
✓ Branch 62 → 63 taken 16668 times.
✓ Branch 62 → 247 taken 5556 times.
|
27780 | bspline_quadrature%nodes = gl_nodes_3 |
| 115 | case (4) | ||
| 116 |
3/6✗ Branch 64 → 65 not taken.
✓ Branch 64 → 68 taken 6605 times.
✗ Branch 65 → 66 not taken.
✗ Branch 65 → 67 not taken.
✓ Branch 69 → 70 taken 26420 times.
✓ Branch 69 → 71 taken 6605 times.
|
39630 | bspline_quadrature%weights = gl_weights_4 * bspline%h / 2 |
| 117 |
3/6✗ Branch 71 → 72 not taken.
✓ Branch 71 → 75 taken 6605 times.
✗ Branch 72 → 73 not taken.
✗ Branch 72 → 74 not taken.
✓ Branch 76 → 77 taken 26420 times.
✓ Branch 76 → 247 taken 6605 times.
|
39630 | bspline_quadrature%nodes = gl_nodes_4 |
| 118 | case (5) | ||
| 119 |
3/6✗ Branch 78 → 79 not taken.
✓ Branch 78 → 82 taken 10348 times.
✗ Branch 79 → 80 not taken.
✗ Branch 79 → 81 not taken.
✓ Branch 83 → 84 taken 51740 times.
✓ Branch 83 → 85 taken 10348 times.
|
72436 | bspline_quadrature%weights = gl_weights_5 * bspline%h / 2 |
| 120 |
3/6✗ Branch 85 → 86 not taken.
✓ Branch 85 → 89 taken 10348 times.
✗ Branch 86 → 87 not taken.
✗ Branch 86 → 88 not taken.
✓ Branch 90 → 91 taken 51740 times.
✓ Branch 90 → 247 taken 10348 times.
|
72436 | bspline_quadrature%nodes = gl_nodes_5 |
| 121 | case (6) | ||
| 122 |
3/6✗ Branch 92 → 93 not taken.
✓ Branch 92 → 96 taken 4400 times.
✗ Branch 93 → 94 not taken.
✗ Branch 93 → 95 not taken.
✓ Branch 97 → 98 taken 26400 times.
✓ Branch 97 → 99 taken 4400 times.
|
35200 | bspline_quadrature%weights = gl_weights_6 * bspline%h / 2 |
| 123 |
3/6✗ Branch 99 → 100 not taken.
✓ Branch 99 → 103 taken 4400 times.
✗ Branch 100 → 101 not taken.
✗ Branch 100 → 102 not taken.
✓ Branch 104 → 105 taken 26400 times.
✓ Branch 104 → 247 taken 4400 times.
|
35200 | bspline_quadrature%nodes = gl_nodes_6 |
| 124 | case (7) | ||
| 125 |
3/6✗ Branch 106 → 107 not taken.
✓ Branch 106 → 110 taken 2341 times.
✗ Branch 107 → 108 not taken.
✗ Branch 107 → 109 not taken.
✓ Branch 111 → 112 taken 16387 times.
✓ Branch 111 → 113 taken 2341 times.
|
21069 | bspline_quadrature%weights = gl_weights_7 * bspline%h / 2 |
| 126 |
3/6✗ Branch 113 → 114 not taken.
✓ Branch 113 → 117 taken 2341 times.
✗ Branch 114 → 115 not taken.
✗ Branch 114 → 116 not taken.
✓ Branch 118 → 119 taken 16387 times.
✓ Branch 118 → 247 taken 2341 times.
|
21069 | bspline_quadrature%nodes = gl_nodes_7 |
| 127 | case (8) | ||
| 128 |
3/6✗ Branch 120 → 121 not taken.
✓ Branch 120 → 124 taken 922 times.
✗ Branch 121 → 122 not taken.
✗ Branch 121 → 123 not taken.
✓ Branch 125 → 126 taken 7376 times.
✓ Branch 125 → 127 taken 922 times.
|
9220 | bspline_quadrature%weights = gl_weights_8 * bspline%h / 2 |
| 129 |
3/6✗ Branch 127 → 128 not taken.
✓ Branch 127 → 131 taken 922 times.
✗ Branch 128 → 129 not taken.
✗ Branch 128 → 130 not taken.
✓ Branch 132 → 133 taken 7376 times.
✓ Branch 132 → 247 taken 922 times.
|
9220 | bspline_quadrature%nodes = gl_nodes_8 |
| 130 | case (9) | ||
| 131 |
3/6✗ Branch 134 → 135 not taken.
✓ Branch 134 → 138 taken 60 times.
✗ Branch 135 → 136 not taken.
✗ Branch 135 → 137 not taken.
✓ Branch 139 → 140 taken 540 times.
✓ Branch 139 → 141 taken 60 times.
|
660 | bspline_quadrature%weights = gl_weights_9 * bspline%h / 2 |
| 132 |
3/6✗ Branch 141 → 142 not taken.
✓ Branch 141 → 145 taken 60 times.
✗ Branch 142 → 143 not taken.
✗ Branch 142 → 144 not taken.
✓ Branch 146 → 147 taken 540 times.
✓ Branch 146 → 247 taken 60 times.
|
660 | bspline_quadrature%nodes = gl_nodes_9 |
| 133 | case (10) | ||
| 134 |
3/6✗ Branch 148 → 149 not taken.
✓ Branch 148 → 152 taken 22 times.
✗ Branch 149 → 150 not taken.
✗ Branch 149 → 151 not taken.
✓ Branch 153 → 154 taken 220 times.
✓ Branch 153 → 155 taken 22 times.
|
264 | bspline_quadrature%weights = gl_weights_10 * bspline%h / 2 |
| 135 |
3/6✗ Branch 155 → 156 not taken.
✓ Branch 155 → 159 taken 22 times.
✗ Branch 156 → 157 not taken.
✗ Branch 156 → 158 not taken.
✓ Branch 160 → 161 taken 220 times.
✓ Branch 160 → 247 taken 22 times.
|
264 | bspline_quadrature%nodes = gl_nodes_10 |
| 136 | case (11) | ||
| 137 | ✗ | bspline_quadrature%weights = gl_weights_11 * bspline%h / 2 | |
| 138 | ✗ | bspline_quadrature%nodes = gl_nodes_11 | |
| 139 | case (12) | ||
| 140 | ✗ | bspline_quadrature%weights = gl_weights_12 * bspline%h / 2 | |
| 141 | ✗ | bspline_quadrature%nodes = gl_nodes_12 | |
| 142 | case (13) | ||
| 143 | ✗ | bspline_quadrature%weights = gl_weights_13 * bspline%h / 2 | |
| 144 | ✗ | bspline_quadrature%nodes = gl_nodes_13 | |
| 145 | case (14) | ||
| 146 | ✗ | bspline_quadrature%weights = gl_weights_14 * bspline%h / 2 | |
| 147 | ✗ | bspline_quadrature%nodes = gl_nodes_14 | |
| 148 | case (15) | ||
| 149 | ✗ | bspline_quadrature%weights = gl_weights_15 * bspline%h / 2 | |
| 150 | ✗ | bspline_quadrature%nodes = gl_nodes_15 | |
| 151 | case (16) | ||
| 152 | ✗ | bspline_quadrature%weights = gl_weights_16 * bspline%h / 2 | |
| 153 | ✗ | bspline_quadrature%nodes = gl_nodes_16 | |
| 154 | case default | ||
| 155 |
10/17✓ Branch 21 → 22 taken 2611 times.
✓ Branch 21 → 36 taken 2541 times.
✓ Branch 21 → 50 taken 5556 times.
✓ Branch 21 → 64 taken 6605 times.
✓ Branch 21 → 78 taken 10348 times.
✓ Branch 21 → 92 taken 4400 times.
✓ Branch 21 → 106 taken 2341 times.
✓ Branch 21 → 120 taken 922 times.
✓ Branch 21 → 134 taken 60 times.
✓ Branch 21 → 148 taken 22 times.
✗ Branch 21 → 162 not taken.
✗ Branch 21 → 176 not taken.
✗ Branch 21 → 190 not taken.
✗ Branch 21 → 204 not taken.
✗ Branch 21 → 218 not taken.
✗ Branch 21 → 232 not taken.
✗ Branch 21 → 246 not taken.
|
35406 | error stop "Number of quadrature points not supported" |
| 156 | end select | ||
| 157 | |||
| 158 | ! allocate(bspline_quadrature%bspline_at_nodes(1:n_quad_, 1:bspline%degree+1, 0:bspline%degree)) | ||
| 159 | if (allocated(bspline_quadrature%bspline_at_nodes)) deallocate (bspline_quadrature%bspline_at_nodes) | ||
| 160 |
11/16✓ Branch 247 → 248 taken 1168 times.
✓ Branch 247 → 249 taken 34238 times.
✓ Branch 249 → 248 taken 34238 times.
✗ Branch 249 → 250 not taken.
✓ Branch 250 → 251 taken 1168 times.
✓ Branch 250 → 252 taken 34238 times.
✓ Branch 252 → 251 taken 34238 times.
✗ Branch 252 → 253 not taken.
✓ Branch 253 → 254 taken 35406 times.
✗ Branch 253 → 255 not taken.
✓ Branch 255 → 256 taken 34238 times.
✓ Branch 255 → 257 taken 1168 times.
✗ Branch 257 → 258 not taken.
✓ Branch 257 → 259 taken 35406 times.
✗ Branch 259 → 260 not taken.
✓ Branch 259 → 261 taken 35406 times.
|
141624 | allocate (bspline_quadrature%bspline_at_nodes(0:bspline%degree, 0:bspline%degree)) |
| 161 | |||
| 162 |
2/2✓ Branch 261 → 262 taken 18389 times.
✓ Branch 261 → 263 taken 17017 times.
|
35406 | if (bspline%is_periodic) then |
| 163 | ! All of the B-splines are the same, and 1 + degree is the first one without any boundary involvement | ||
| 164 | j0 = bspline%degree | ||
| 165 | j1 = bspline%degree | ||
| 166 | else | ||
| 167 | ! The first degree B-splines have some boundary involvement (as do the final degree B-splines) | ||
| 168 | ! and the 1 + degree B-spline is the first one without any boundary involvement | ||
| 169 | ! and thus represents all the 'interior' B-splines | ||
| 170 | j0 = 0 | ||
| 171 | j1 = bspline%degree | ||
| 172 | end if | ||
| 173 |
2/2✓ Branch 264 → 265 taken 79018 times.
✓ Branch 264 → 275 taken 35406 times.
|
114424 | do j = j0, j1 |
| 174 |
2/2✓ Branch 266 → 267 taken 333383 times.
✓ Branch 266 → 273 taken 79018 times.
|
447807 | do j_minus_i = 0, bspline%degree |
| 175 |
2/2✓ Branch 267 → 268 taken 1778368 times.
✓ Branch 267 → 271 taken 333383 times.
|
2190769 | do ndx = 1, n_quad_ |
| 176 | 1778368 | x = (1 + bspline_quadrature%nodes(ndx)) / 2 | |
| 177 | bspline_quadrature%bspline_at_nodes(j, j_minus_i)%data(ndx) = bspline_eval(x, j, j_minus_i, bspline%degree, & | ||
| 178 | 1778368 | & bspline%is_scaled) | |
| 179 |
2/2✓ Branch 268 → 269 taken 284272 times.
✓ Branch 268 → 270 taken 1494096 times.
|
2111751 | if (bspline%is_scaled) then |
| 180 | bspline_quadrature%bspline_at_nodes(j, j_minus_i)%data(ndx) = & | ||
| 181 | 284272 | & bspline_quadrature%bspline_at_nodes(j, j_minus_i)%data(ndx) * bspline%nr_intervals | |
| 182 | end if | ||
| 183 | end do | ||
| 184 | end do | ||
| 185 | end do | ||
| 186 | |||
| 187 | 35406 | end function init_bspline_quadrature | |
| 188 | |||
| 189 | !> @brief Destroy a B-spline quadrature object | ||
| 190 | !> | ||
| 191 | !> @param[inout] this The B-spline quadrature object to destroy | ||
| 192 | 16947 | subroutine destroy_bspline_quadrature(this) | |
| 193 | implicit none | ||
| 194 | |||
| 195 | class(BSplineQuadrature), intent(inout) :: this | ||
| 196 | |||
| 197 |
1/2✓ Branch 2 → 3 taken 16947 times.
✗ Branch 2 → 4 not taken.
|
16947 | if (allocated(this%weights)) deallocate (this%weights) |
| 198 |
1/2✓ Branch 4 → 5 taken 16947 times.
✗ Branch 4 → 6 not taken.
|
16947 | if (allocated(this%nodes)) deallocate (this%nodes) |
| 199 |
1/2✓ Branch 6 → 7 taken 16947 times.
✗ Branch 6 → 8 not taken.
|
16947 | if (allocated(this%bspline_at_nodes)) deallocate (this%bspline_at_nodes) |
| 200 | 16947 | end subroutine destroy_bspline_quadrature | |
| 201 | |||
| 202 | ! TODO the 'final' functionality is not reliably supported by all (any?) compilers, so we do not use it here | ||
| 203 | |||
| 204 | ! !> @brief Destroy a B-spline quadrature object when it goes out of scope (final) | ||
| 205 | ! !> | ||
| 206 | ! !> @param[inout] this The B-spline quadrature object to destroy | ||
| 207 | ! subroutine destroy_bspline_quadrature_final(this) | ||
| 208 | ! implicit none | ||
| 209 | |||
| 210 | ! type(BSplineQuadrature), intent(inout) :: this | ||
| 211 | |||
| 212 | ! call this%destroy() | ||
| 213 | ! end subroutine destroy_bspline_quadrature_final | ||
| 214 | |||
| 215 | !> @brief Copy a B-spline quadrature object | ||
| 216 | !> | ||
| 217 | !> @param[out] to The B-spline quadrature object to copy to | ||
| 218 | !> @param[in] from The B-spline quadrature object to copy | ||
| 219 |
1/2✓ Branch 2 → 3 taken 35406 times.
✗ Branch 2 → 5 not taken.
|
35406 | subroutine copy_bspline_quadrature(to, from) |
| 220 | implicit none | ||
| 221 | |||
| 222 | class(BSplineQuadrature), intent(out) :: to | ||
| 223 | class(BSplineQuadrature), intent(in) :: from | ||
| 224 | |||
| 225 | integer :: j, j_minus_i, j0, j1 | ||
| 226 | |||
| 227 | 35406 | to%bspline = from%bspline | |
| 228 | 35406 | to%n_quad = from%n_quad | |
| 229 | |||
| 230 |
1/2✗ Branch 5 → 6 not taken.
✓ Branch 5 → 7 taken 35406 times.
|
35406 | if (allocated(to%weights)) deallocate (to%weights) |
| 231 |
7/14✗ Branch 7 → 8 not taken.
✓ Branch 7 → 9 taken 35406 times.
✓ Branch 9 → 8 taken 35406 times.
✗ Branch 9 → 10 not taken.
✓ Branch 10 → 11 taken 35406 times.
✗ Branch 10 → 12 not taken.
✓ Branch 12 → 13 taken 35406 times.
✗ Branch 12 → 14 not taken.
✗ Branch 14 → 15 not taken.
✓ Branch 14 → 16 taken 35406 times.
✗ Branch 16 → 17 not taken.
✓ Branch 16 → 18 taken 35406 times.
✗ Branch 18 → 19 not taken.
✓ Branch 18 → 20 taken 35406 times.
|
106218 | allocate (to%weights(from%n_quad)) |
| 232 |
4/12✓ Branch 20 → 21 taken 35406 times.
✗ Branch 20 → 22 not taken.
✗ Branch 21 → 22 not taken.
✓ Branch 21 → 30 taken 35406 times.
✗ Branch 22 → 23 not taken.
✗ Branch 22 → 24 not taken.
✗ Branch 24 → 25 not taken.
✗ Branch 24 → 26 not taken.
✗ Branch 26 → 27 not taken.
✗ Branch 26 → 28 not taken.
✓ Branch 31 → 32 taken 153444 times.
✓ Branch 31 → 33 taken 35406 times.
|
224256 | to%weights = from%weights |
| 233 | |||
| 234 |
1/2✗ Branch 33 → 34 not taken.
✓ Branch 33 → 35 taken 35406 times.
|
35406 | if (allocated(to%nodes)) deallocate (to%nodes) |
| 235 |
7/14✗ Branch 35 → 36 not taken.
✓ Branch 35 → 37 taken 35406 times.
✓ Branch 37 → 36 taken 35406 times.
✗ Branch 37 → 38 not taken.
✓ Branch 38 → 39 taken 35406 times.
✗ Branch 38 → 40 not taken.
✓ Branch 40 → 41 taken 35406 times.
✗ Branch 40 → 42 not taken.
✗ Branch 42 → 43 not taken.
✓ Branch 42 → 44 taken 35406 times.
✗ Branch 44 → 45 not taken.
✓ Branch 44 → 46 taken 35406 times.
✗ Branch 46 → 47 not taken.
✓ Branch 46 → 48 taken 35406 times.
|
106218 | allocate (to%nodes(from%n_quad)) |
| 236 |
4/12✓ Branch 48 → 49 taken 35406 times.
✗ Branch 48 → 50 not taken.
✗ Branch 49 → 50 not taken.
✓ Branch 49 → 58 taken 35406 times.
✗ Branch 50 → 51 not taken.
✗ Branch 50 → 52 not taken.
✗ Branch 52 → 53 not taken.
✗ Branch 52 → 54 not taken.
✗ Branch 54 → 55 not taken.
✗ Branch 54 → 56 not taken.
✓ Branch 59 → 60 taken 153444 times.
✓ Branch 59 → 61 taken 35406 times.
|
224256 | to%nodes = from%nodes |
| 237 | |||
| 238 |
1/2✗ Branch 61 → 62 not taken.
✓ Branch 61 → 63 taken 35406 times.
|
35406 | if (allocated(to%bspline_at_nodes)) deallocate (to%bspline_at_nodes) |
| 239 |
12/18✓ Branch 63 → 64 taken 1168 times.
✓ Branch 63 → 65 taken 34238 times.
✓ Branch 65 → 64 taken 34238 times.
✗ Branch 65 → 66 not taken.
✓ Branch 66 → 67 taken 1168 times.
✓ Branch 66 → 68 taken 34238 times.
✓ Branch 68 → 67 taken 34238 times.
✗ Branch 68 → 69 not taken.
✓ Branch 69 → 70 taken 35406 times.
✗ Branch 69 → 71 not taken.
✓ Branch 71 → 72 taken 34238 times.
✓ Branch 71 → 73 taken 1168 times.
✗ Branch 73 → 74 not taken.
✓ Branch 73 → 75 taken 35406 times.
✗ Branch 75 → 76 not taken.
✓ Branch 75 → 77 taken 35406 times.
✗ Branch 77 → 78 not taken.
✓ Branch 77 → 79 taken 35406 times.
|
141624 | allocate (to%bspline_at_nodes(0:from%bspline%degree, 0:from%bspline%degree)) |
| 240 | |||
| 241 |
2/2✓ Branch 79 → 80 taken 18389 times.
✓ Branch 79 → 81 taken 17017 times.
|
35406 | if (from%bspline%is_periodic) then |
| 242 | ! All of the B-splines are the same, and 1 + degree is the first one without any boundary involvement | ||
| 243 | j0 = from%bspline%degree | ||
| 244 | j1 = from%bspline%degree | ||
| 245 | else | ||
| 246 | ! The first degree B-splines have some boundary involvement (as do the final degree B-splines) | ||
| 247 | ! and the 1 + degree B-spline is the first one without any boundary involvement | ||
| 248 | ! and thus represents all the 'interior' B-splines | ||
| 249 | j0 = 0 | ||
| 250 | j1 = from%bspline%degree | ||
| 251 | end if | ||
| 252 |
2/2✓ Branch 82 → 83 taken 79018 times.
✓ Branch 82 → 91 taken 35406 times.
|
114424 | do j = j0, j1 |
| 253 |
2/2✓ Branch 84 → 85 taken 333383 times.
✓ Branch 84 → 89 taken 79018 times.
|
447807 | do j_minus_i = 0, from%bspline%degree |
| 254 |
2/2✓ Branch 86 → 87 taken 1778368 times.
✓ Branch 86 → 88 taken 333383 times.
|
2190769 | to%bspline_at_nodes(j, j_minus_i)%data(1:from%n_quad) = from%bspline_at_nodes(j, j_minus_i)%data(1:from%n_quad) |
| 255 | end do | ||
| 256 | end do | ||
| 257 | |||
| 258 | 35406 | end subroutine copy_bspline_quadrature | |
| 259 | |||
| 260 | !> @brief Determine the number of quadrature points needed for exact integration of a single B-spline space | ||
| 261 | !> | ||
| 262 | !> @param[in] bspline The B-spline space for which the quadrature is defined | ||
| 263 | !> | ||
| 264 | !> @return The number of quadrature points needed for exact integration | ||
| 265 | ✗ | pure integer function n_quad_exact_single(bspline) result(ans) | |
| 266 | implicit none | ||
| 267 | |||
| 268 | type(BSplineSpace), intent(in) :: bspline | ||
| 269 | |||
| 270 | ✗ | ans = 1 + int(bspline%degree / 2) | |
| 271 | ✗ | end function n_quad_exact_single | |
| 272 | |||
| 273 | !> @brief Determine the number of quadrature points needed for exact integration of a product of two B-spline spaces | ||
| 274 | !> | ||
| 275 | !> @param[in] bspline1 The first B-spline space for which the quadrature is defined | ||
| 276 | !> @param[in] bspline2 The second B-spline space for which the quadrature is defined | ||
| 277 | !> | ||
| 278 | !> @return The number of quadrature points needed for exact integration | ||
| 279 | 25831 | pure integer function n_quad_exact_double(bspline1, bspline2) result(ans) | |
| 280 | implicit none | ||
| 281 | |||
| 282 | type(BSplineSpace), intent(in) :: bspline1, bspline2 | ||
| 283 | |||
| 284 | ✗ | ans = 1 + max(bspline1%degree, bspline2%degree) | |
| 285 | 25831 | end function n_quad_exact_double | |
| 286 | |||
| 287 | !> @brief Evaluate the B-spline at a quadrature node | ||
| 288 | !> | ||
| 289 | !> @param[in] this The B-spline quadrature object | ||
| 290 | !> @param[in] j The index of the B-spline | ||
| 291 | !> @param[in] j_minus_i The index of the B-spline minus the index of the interval | ||
| 292 | !> @param[in] ndx The index of the quadrature node | ||
| 293 | !> | ||
| 294 | !> @return The value of the B-spline at the quadrature node | ||
| 295 | 1319723844 | pure real(wp) function evaluate(this, j, j_minus_i, ndx) result(ans) | |
| 296 | use m_bspline_basis, only: get_j_internal | ||
| 297 | implicit none | ||
| 298 | |||
| 299 | type(BSplineQuadrature), intent(in) :: this | ||
| 300 | integer, intent(in) :: j, j_minus_i, ndx | ||
| 301 | |||
| 302 | integer :: j_internal, ndx_actual, j_internal_minus_i | ||
| 303 | logical :: bspline_mirrored | ||
| 304 | |||
| 305 | 1319723844 | call get_j_internal(this%bspline, j, j_internal, bspline_mirrored) | |
| 306 |
2/2✓ Branch 3 → 4 taken 560995730 times.
✓ Branch 3 → 5 taken 758728114 times.
|
1319723844 | if (bspline_mirrored) then |
| 307 | 560995730 | ndx_actual = this%n_quad - ndx + 1 | |
| 308 | 560995730 | j_internal_minus_i = this%bspline%degree - j_minus_i | |
| 309 | else | ||
| 310 | 758728114 | ndx_actual = ndx | |
| 311 | 758728114 | j_internal_minus_i = j_minus_i | |
| 312 | end if | ||
| 313 | |||
| 314 | 1319723844 | ans = this%bspline_at_nodes(j_internal, j_internal_minus_i)%data(ndx_actual) | |
| 315 | 1319723844 | end function evaluate | |
| 316 | |||
| 317 | !> @brief Quadrature over a B-spline | ||
| 318 | !> | ||
| 319 | !> @param[in] this The B-spline quadrature object | ||
| 320 | !> @param[in] j The index of the B-spline | ||
| 321 | !> @param[in] userfun _(optional)_ A user-defined function to be multiplied with the B-spline | ||
| 322 | !> | ||
| 323 | !> @return The result of the quadrature rule | ||
| 324 | 42055 | pure real(wp) function quadrature_1d(this, j, userfun) result(ans) | |
| 325 | use m_bspline_basis, only: get_index_ranges | ||
| 326 | use m_common, only: user_function_1d_interface | ||
| 327 | use m_quadrature_data | ||
| 328 | |||
| 329 | type(BSplineQuadrature), intent(in) :: this | ||
| 330 | integer, intent(in) :: j | ||
| 331 | procedure(user_function_1d_interface), optional :: userfun | ||
| 332 | |||
| 333 | integer :: i, i0_min, i0_max, i1_min, i1_max | ||
| 334 | |||
| 335 | 42055 | call get_index_ranges(this%bspline, j, i0_min, i0_max, i1_min, i1_max) | |
| 336 | |||
| 337 | ans = 0._wp | ||
| 338 |
2/2✓ Branch 4 → 5 taken 136294 times.
✓ Branch 4 → 9 taken 42055 times.
|
178349 | do i = i0_min, i0_max |
| 339 |
2/2✓ Branch 5 → 6 taken 78570 times.
✓ Branch 5 → 7 taken 57724 times.
|
178349 | if (present(userfun)) then |
| 340 | 78570 | ans = ans + quadrature_1d_interval(this, j, i, userfun) | |
| 341 | else | ||
| 342 | 57724 | ans = ans + quadrature_1d_interval(this, j, i) | |
| 343 | end if | ||
| 344 | end do | ||
| 345 |
2/2✓ Branch 11 → 12 taken 14178 times.
✓ Branch 11 → 16 taken 42055 times.
|
56233 | do i = i1_min, i1_max |
| 346 |
2/2✓ Branch 12 → 13 taken 7579 times.
✓ Branch 12 → 14 taken 6599 times.
|
56233 | if (present(userfun)) then |
| 347 | 7579 | ans = ans + quadrature_1d_interval(this, j, i, userfun) | |
| 348 | else | ||
| 349 | 6599 | ans = ans + quadrature_1d_interval(this, j, i) | |
| 350 | end if | ||
| 351 | end do | ||
| 352 | |||
| 353 | 42055 | end function quadrature_1d | |
| 354 | |||
| 355 | !> @brief Quadrature over a B-spline on a single interval | ||
| 356 | !> | ||
| 357 | !> @param[in] this The B-spline quadrature object | ||
| 358 | !> @param[in] j The index of the B-spline | ||
| 359 | !> @param[in] i The index of the interval | ||
| 360 | !> @param[in] userfun _(optional)_ A user-defined function to be multiplied with the B-spline | ||
| 361 | !> | ||
| 362 | !> @return The result of the quadrature rule on the interval | ||
| 363 | 295322 | pure real(wp) function quadrature_1d_interval(this, j, i, userfun) result(ans) | |
| 364 | use m_bspline_basis, only: get_j_minus_i_quad | ||
| 365 | use m_common, only: user_function_1d_interface | ||
| 366 | |||
| 367 | implicit none | ||
| 368 | |||
| 369 | type(BSplineQuadrature), intent(in) :: this | ||
| 370 | integer, intent(in) :: j, i | ||
| 371 | procedure(user_function_1d_interface), optional :: userfun | ||
| 372 | |||
| 373 | integer :: j_minus_i, ndx | ||
| 374 | real(wp) :: val, x | ||
| 375 | |||
| 376 | 295322 | j_minus_i = get_j_minus_i_quad(this%bspline, j, i) | |
| 377 | |||
| 378 | ans = 0._wp | ||
| 379 |
4/4✓ Branch 2 → 3 taken 228260 times.
✓ Branch 2 → 11 taken 67062 times.
✓ Branch 3 → 4 taken 163838 times.
✓ Branch 3 → 11 taken 64422 times.
|
295322 | if (0 <= j_minus_i .and. j_minus_i <= this%bspline%degree) then |
| 380 |
2/2✓ Branch 5 → 6 taken 440084 times.
✓ Branch 5 → 10 taken 163838 times.
|
603922 | do ndx = 1, this%n_quad |
| 381 | 440084 | val = evaluate(this, j, j_minus_i, ndx) | |
| 382 |
2/2✓ Branch 6 → 7 taken 271866 times.
✓ Branch 6 → 9 taken 168218 times.
|
440084 | if (present(userfun)) then |
| 383 | 271866 | x = (i + .5_wp + this%nodes(ndx) / 2) * this%bspline%h | |
| 384 | 271866 | val = val * userfun(x) | |
| 385 | end if | ||
| 386 | 603922 | ans = ans + val * this%weights(ndx) | |
| 387 | end do | ||
| 388 | end if | ||
| 389 | |||
| 390 | 295322 | end function quadrature_1d_interval | |
| 391 | |||
| 392 | !> @brief Quadrature over a product of two B-splines | ||
| 393 | !> | ||
| 394 | !> @param[in] this1 The first B-spline quadrature object | ||
| 395 | !> @param[in] this2 The second B-spline quadrature object | ||
| 396 | !> @param[in] j1 The index of the first B-spline | ||
| 397 | !> @param[in] j2 The index of the second B-spline | ||
| 398 | !> @param[in] userfun _(optional)_ A user-defined function to be multiplied with the product of the B-splines | ||
| 399 | !> | ||
| 400 | !> @return The result of the quadrature rule on the product of the B-splines | ||
| 401 | 4078300 | pure real(wp) function quadrature_product_1d(this1, this2, j1, j2, userfun) result(ans) | |
| 402 | use m_bspline_basis, only: get_index_ranges | ||
| 403 | use m_common, only: user_function_1d_interface | ||
| 404 | use m_quadrature_data | ||
| 405 | |||
| 406 | type(BSplineQuadrature), intent(in) :: this1, this2 | ||
| 407 | integer, intent(in) :: j1, j2 | ||
| 408 | procedure(user_function_1d_interface), optional :: userfun | ||
| 409 | |||
| 410 | integer :: i, i0_min, i0_max, i1_min, i1_max | ||
| 411 | |||
| 412 | 4078300 | call get_index_ranges(this1%bspline, this2%bspline, j1, j2, i0_min, i0_max, i1_min, i1_max) | |
| 413 | |||
| 414 | ans = 0._wp | ||
| 415 |
2/2✓ Branch 4 → 5 taken 5519659 times.
✓ Branch 4 → 9 taken 4078300 times.
|
9597959 | do i = i0_min, i0_max |
| 416 |
2/2✓ Branch 5 → 6 taken 2736264 times.
✓ Branch 5 → 7 taken 2783395 times.
|
9597959 | if (present(userfun)) then |
| 417 | 2736264 | ans = ans + quadrature_product_1d_interval(this1, this2, j1, j2, i, userfun) | |
| 418 | else | ||
| 419 | 2783395 | ans = ans + quadrature_product_1d_interval(this1, this2, j1, j2, i) | |
| 420 | end if | ||
| 421 | end do | ||
| 422 |
1/2✗ Branch 11 → 12 not taken.
✓ Branch 11 → 16 taken 4078300 times.
|
4078300 | do i = i1_min, i1_max |
| 423 |
0/2✗ Branch 12 → 13 not taken.
✗ Branch 12 → 14 not taken.
|
4078300 | if (present(userfun)) then |
| 424 | ✗ | ans = ans + quadrature_product_1d_interval(this1, this2, j1, j2, i, userfun) | |
| 425 | else | ||
| 426 | ✗ | ans = ans + quadrature_product_1d_interval(this1, this2, j1, j2, i) | |
| 427 | end if | ||
| 428 | end do | ||
| 429 | |||
| 430 | 4078300 | end function quadrature_product_1d | |
| 431 | |||
| 432 | !> @brief Quadrature over a product of two B-splines on a single interval | ||
| 433 | !> | ||
| 434 | !> @param[in] this1 The first B-spline quadrature object | ||
| 435 | !> @param[in] this2 The second B-spline quadrature object | ||
| 436 | !> @param[in] j1 The index of the first B-spline | ||
| 437 | !> @param[in] j2 The index of the second B-spline | ||
| 438 | !> @param[in] i The index of the interval | ||
| 439 | !> @param[in] userfun _(optional)_ A user-defined function to be multiplied with the product of the B-splines | ||
| 440 | !> | ||
| 441 | !> @return The result of the quadrature rule on the product of the B-splines on the interval | ||
| 442 | 5519659 | pure real(wp) function quadrature_product_1d_interval(this1, this2, j1, j2, i, userfun) result(ans) | |
| 443 | use m_bspline_basis, only: get_j_minus_i_quad | ||
| 444 | use m_common, only: user_function_1d_interface | ||
| 445 | |||
| 446 | implicit none | ||
| 447 | |||
| 448 | type(BSplineQuadrature), intent(in) :: this1, this2 | ||
| 449 | integer, intent(in) :: j1, j2, i | ||
| 450 | procedure(user_function_1d_interface), optional :: userfun | ||
| 451 | |||
| 452 | integer :: j1_minus_i, j2_minus_i, ndx | ||
| 453 | real(wp) :: val, x | ||
| 454 | |||
| 455 | 5519659 | j1_minus_i = get_j_minus_i_quad(this1%bspline, j1, i) | |
| 456 | 5519659 | j2_minus_i = get_j_minus_i_quad(this2%bspline, j2, i) | |
| 457 | |||
| 458 | ans = 0._wp | ||
| 459 | if (0 <= j1_minus_i .and. j1_minus_i <= this1%bspline%degree .and. & | ||
| 460 |
4/8✓ Branch 2 → 3 taken 5519659 times.
✗ Branch 2 → 13 not taken.
✓ Branch 3 → 4 taken 5519659 times.
✗ Branch 3 → 13 not taken.
✓ Branch 4 → 5 taken 5519659 times.
✗ Branch 4 → 13 not taken.
✓ Branch 5 → 6 taken 5519659 times.
✗ Branch 5 → 13 not taken.
|
5519659 | 0 <= j2_minus_i .and. j2_minus_i <= this2%bspline%degree) then |
| 461 |
2/2✓ Branch 7 → 8 taken 26597093 times.
✓ Branch 7 → 12 taken 5519659 times.
|
32116752 | do ndx = 1, this1%n_quad |
| 462 | 26597093 | val = evaluate(this1, j1, j1_minus_i, ndx) * evaluate(this2, j2, j2_minus_i, ndx) | |
| 463 |
2/2✓ Branch 8 → 9 taken 13184604 times.
✓ Branch 8 → 11 taken 13412489 times.
|
26597093 | if (present(userfun)) then |
| 464 | 13184604 | x = (i - .5_wp + this1%nodes(ndx) / 2) * this1%bspline%h | |
| 465 | 13184604 | val = val * userfun(x) | |
| 466 | end if | ||
| 467 | 32116752 | ans = ans + val * this1%weights(ndx) | |
| 468 | end do | ||
| 469 | end if | ||
| 470 | |||
| 471 | 5519659 | end function quadrature_product_1d_interval | |
| 472 | |||
| 473 |
12/26__m_bspline_quadrature_MOD___copy_m_bspline_quadrature_Bsplinequadrature:
✗ Branch 2 → 3 not taken.
✗ Branch 2 → 12 not taken.
✗ Branch 3 → 4 not taken.
✗ Branch 3 → 5 not taken.
✗ Branch 6 → 7 not taken.
✗ Branch 6 → 8 not taken.
✗ Branch 9 → 10 not taken.
✗ Branch 9 → 11 not taken.
__m_bspline_quadrature_MOD___final_m_bspline_quadrature_Bsplinequadrature:
✗ Branch 3 → 4 not taken.
✓ Branch 3 → 7 taken 35406 times.
✗ Branch 4 → 5 not taken.
✗ Branch 4 → 6 not taken.
✓ Branch 8 → 9 taken 35406 times.
✓ Branch 8 → 23 taken 35406 times.
✗ Branch 10 → 11 not taken.
✓ Branch 10 → 12 taken 35406 times.
✓ Branch 12 → 13 taken 35406 times.
✗ Branch 12 → 18 not taken.
✓ Branch 13 → 14 taken 529 times.
✓ Branch 13 → 15 taken 34877 times.
✓ Branch 15 → 16 taken 529 times.
✓ Branch 15 → 17 taken 34877 times.
✓ Branch 18 → 19 taken 35406 times.
✗ Branch 18 → 22 not taken.
✓ Branch 19 → 20 taken 529 times.
✓ Branch 19 → 21 taken 34877 times.
|
106218 | end module m_bspline_quadrature |
| 474 |