bspline/m_bspline_basis.f90
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | !> The B-spline basis module | ||
| 2 | !> | ||
| 3 | !> This module provides | ||
| 4 | !> - BSplineSpace: the space of B-spline basis functions | ||
| 5 | !> - BSplineFun: the B-spline function which is represented by a linear combination of the basis functions from a space | ||
| 6 | !> | ||
| 7 | !> The space is defined by the number of _uniform_ intervals, the B-spline degree, as well as the boundary type: | ||
| 8 | !> - The B-spline space can be _clamped_ meaning that only the first B-spline is nonzero at the left boundary and only the last | ||
| 9 | !> B-spline is nonzero at the right boundary. | ||
| 10 | !> - Or the B-spline space can be _periodic_ meaning that the first 'degree' B-splines are combined with the last 'degree' | ||
| 11 | !> B-splines, thereby enforcing periodicity. | ||
| 12 | module m_bspline_basis | ||
| 13 | use m_common, only: wp | ||
| 14 | use m_bspline_functions, only: MAX_BSPLINE_DEGREE | ||
| 15 | #include "petsc.fi" | ||
| 16 | |||
| 17 | implicit none | ||
| 18 | |||
| 19 | private | ||
| 20 | |||
| 21 | public :: BSplineSpace, BSplineFun | ||
| 22 | public :: MAX_BSPLINE_DEGREE | ||
| 23 | public :: derivative, anti_derivative | ||
| 24 | public :: evaluate, evaluate_support, quadrature, quadrature_product, evaluate_single | ||
| 25 | public :: size, get_index_ranges, get_j_minus_i_quad, get_internal_and_actual_inds, get_j_internal, get_petsc | ||
| 26 | public :: get_interval_index | ||
| 27 | public :: find_root_brent, find_root_newton | ||
| 28 | |||
| 29 | !> @brief The B-spline space type | ||
| 30 | !> | ||
| 31 | !> The B-spline space represents the space of B-spline basis functions. | ||
| 32 | type :: BSplineSpace | ||
| 33 | !> The number of B-splines | ||
| 34 | integer :: nr_bsplines | ||
| 35 | |||
| 36 | !> The number of intervals | ||
| 37 | integer :: nr_intervals | ||
| 38 | |||
| 39 | !> The degree of the B-splines | ||
| 40 | integer :: degree | ||
| 41 | |||
| 42 | !> The mesh-width | ||
| 43 | real(wp) :: h | ||
| 44 | |||
| 45 | !> Flag indicating whether the B-splines are periodic | ||
| 46 | logical :: is_periodic | ||
| 47 | |||
| 48 | !> Flag indicating whether the B-splines are scaled (such that the derivative of an unscaled degree d B-spline is simply given | ||
| 49 | !> by the difference of scaled degree d-1 B-splines) | ||
| 50 | logical :: is_scaled | ||
| 51 | contains | ||
| 52 | |||
| 53 | end type BSplineSpace | ||
| 54 | |||
| 55 | !> @brief The B-spline function type | ||
| 56 | !> | ||
| 57 | !> The B-spline function represents a linear combination of B-spline basis functions from a B-spline space. | ||
| 58 | type BSplineFun | ||
| 59 | !> The corresponding B-spline space | ||
| 60 | type(BSplineSpace) :: bspline | ||
| 61 | |||
| 62 | !> The coefficients of the B-spline function | ||
| 63 | real(wp), allocatable :: data(:) | ||
| 64 | contains | ||
| 65 | !> Initialize the B-spline function | ||
| 66 | procedure :: init => init_bspline_fun | ||
| 67 | procedure :: destroy => destroy_bspline_fun | ||
| 68 | procedure :: periodicity => bspline_fun_periodicity | ||
| 69 | ! final :: destroy_bspline_fun_final | ||
| 70 | ! TODO implement the assignment(=) operator | ||
| 71 | end type | ||
| 72 | |||
| 73 | !> @brief Compute quadrature of a B-spline | ||
| 74 | interface quadrature | ||
| 75 | procedure quadrature_eval, quadrature_interval | ||
| 76 | end interface | ||
| 77 | |||
| 78 | !> @brief Compute quadrature of the product of two B-splines | ||
| 79 | interface quadrature_product | ||
| 80 | procedure quadrature_product_eval, quadrature_product_interval | ||
| 81 | end interface | ||
| 82 | |||
| 83 | !> @brief Evaluate a B-spline (function) | ||
| 84 | interface evaluate | ||
| 85 | procedure evaluate_1d, evaluate_single | ||
| 86 | end interface | ||
| 87 | |||
| 88 | interface evaluate_support | ||
| 89 | procedure evaluate_support_1d | ||
| 90 | end interface | ||
| 91 | |||
| 92 | !> @brief Get the size of the B-spline space (i.e., the number of linearly independent B-splines) | ||
| 93 | interface size | ||
| 94 | procedure size_1d | ||
| 95 | end interface | ||
| 96 | |||
| 97 | !> @brief Get the index ranges for integrating a B-spline or the product of two B-splines | ||
| 98 | interface get_index_ranges | ||
| 99 | procedure get_index_ranges_single, get_index_ranges_product | ||
| 100 | end interface | ||
| 101 | |||
| 102 | !> @brief Create a B-spline space | ||
| 103 | interface BSplineSpace | ||
| 104 | module procedure init_BSplineSpace | ||
| 105 | end interface | ||
| 106 | |||
| 107 | !> @brief Convert a B-spline function to a PETSc vector | ||
| 108 | interface get_petsc | ||
| 109 | module procedure get_petsc_1dvec | ||
| 110 | end interface | ||
| 111 | contains | ||
| 112 | |||
| 113 | !> @brief Create a B-spline space | ||
| 114 | !> | ||
| 115 | !> @param[in] nr_intervals The number of intervals | ||
| 116 | !> @param[in] degree The degree of the B-splines | ||
| 117 | !> @param[in] is_scaled _(optional)_ If true, the B-splines are scaled (i.e., this space can be used to express the derivative of | ||
| 118 | !> an unscaled B-spline of one degree higher) | ||
| 119 | !> @param[in] is_periodic _(optional)_ If true, the B-splines are periodic, otherwise they are clamped | ||
| 120 | !> | ||
| 121 | !> @return The B-spline space | ||
| 122 | !> | ||
| 123 | !> @note This requires the number of intervals to be larger than the degree | ||
| 124 | 101374 | pure function init_BSplineSpace(nr_intervals, degree, is_scaled, is_periodic) result(ans) | |
| 125 | |||
| 126 | implicit none | ||
| 127 | |||
| 128 | integer, intent(in) :: nr_intervals | ||
| 129 | integer, intent(in) :: degree | ||
| 130 | logical, intent(in), optional :: is_scaled, is_periodic | ||
| 131 | |||
| 132 | type(BSplineSpace) :: ans | ||
| 133 | character(len=256) :: error_msg | ||
| 134 | |||
| 135 |
2/2✓ Branch 2 → 3 taken 99504 times.
✓ Branch 2 → 4 taken 1870 times.
|
101374 | if (present(is_periodic)) then |
| 136 | 99504 | ans%is_periodic = is_periodic | |
| 137 | else | ||
| 138 | ans%is_periodic = .false. | ||
| 139 | end if | ||
| 140 | |||
| 141 |
2/2✓ Branch 4 → 5 taken 77106 times.
✓ Branch 4 → 6 taken 24268 times.
|
101374 | if (present(is_scaled)) then |
| 142 | 77106 | ans%is_scaled = is_scaled | |
| 143 | else | ||
| 144 | ans%is_scaled = .false. | ||
| 145 | end if | ||
| 146 | |||
| 147 |
5/8✓ Branch 6 → 7 taken 4645 times.
✓ Branch 6 → 15 taken 96729 times.
✓ Branch 7 → 8 taken 4645 times.
✗ Branch 7 → 10 not taken.
✓ Branch 8 → 9 taken 4645 times.
✗ Branch 8 → 10 not taken.
✗ Branch 9 → 10 not taken.
✓ Branch 9 → 15 taken 4645 times.
|
101374 | if (degree < 0 .and. .not. (degree == -1 .and. nr_intervals == 1 .and. ans%is_scaled)) then |
| 148 | ✗ | write (error_msg, '(A,I0)') "ERROR in BSplineSpace: degree must be non-negative, got degree = ", degree | |
| 149 | ✗ | error stop error_msg | |
| 150 | end if | ||
| 151 |
1/2✗ Branch 15 → 16 not taken.
✓ Branch 15 → 23 taken 101374 times.
|
101374 | if (degree > MAX_BSPLINE_DEGREE) then |
| 152 | ✗ | write (error_msg, '(A,I0,A,I0)') "ERROR in BSplineSpace: degree ", degree, & | |
| 153 | ✗ | " exceeds MAX_BSPLINE_DEGREE = ", MAX_BSPLINE_DEGREE | |
| 154 | ✗ | error stop error_msg | |
| 155 | end if | ||
| 156 |
1/2✗ Branch 23 → 24 not taken.
✓ Branch 23 → 31 taken 101374 times.
|
101374 | if (nr_intervals <= degree) then |
| 157 | write (error_msg, '(A,I0,A,I0)') "ERROR in BSplineSpace: periodic B-splines require nr_intervals > degree to avoid "// & | ||
| 158 | ✗ | & " multi-valuedness. Got nr_intervals = ", nr_intervals, " and degree = ", degree | |
| 159 | ✗ | error stop error_msg | |
| 160 | end if | ||
| 161 | |||
| 162 |
2/2✓ Branch 31 → 32 taken 40051 times.
✓ Branch 31 → 33 taken 61323 times.
|
101374 | if (ans%is_periodic) then |
| 163 | ans%nr_bsplines = nr_intervals ! Continuity is imposed | ||
| 164 | else | ||
| 165 | 40051 | ans%nr_bsplines = nr_intervals + degree | |
| 166 | end if | ||
| 167 | ans%nr_intervals = nr_intervals | ||
| 168 | ans%degree = degree | ||
| 169 | 101374 | ans%h = 1._wp / nr_intervals | |
| 170 | |||
| 171 | 101374 | end function | |
| 172 | |||
| 173 | !> @brief Returns the B-spline space containing the derivatives of the given B-spline space | ||
| 174 | !> | ||
| 175 | !> @param[in] bspline The B-spline space | ||
| 176 | !> @param[in] allow_periodic_degree0 Whether to allow periodic degree 0 B-splines in the 0-form space (default: false) | ||
| 177 | !> | ||
| 178 | !> @return The B-spline space containing the derivatives of the given B-spline space | ||
| 179 | !> | ||
| 180 | !> @note This requires the B-spline space to be unscaled | ||
| 181 | 76683 | pure type(BSplineSpace) function derivative(bspline, allow_periodic_degree0) result(dbspline) | |
| 182 | implicit none | ||
| 183 | |||
| 184 | type(BSplineSpace), intent(in) :: bspline | ||
| 185 | logical, intent(in), optional :: allow_periodic_degree0 | ||
| 186 | character(len=256) :: error_msg | ||
| 187 | |||
| 188 | logical :: allow_periodic_degree0_ | ||
| 189 | integer :: derivative_degree | ||
| 190 | |||
| 191 | allow_periodic_degree0_ = .false. | ||
| 192 |
2/2✓ Branch 2 → 3 taken 1770 times.
✓ Branch 2 → 4 taken 74913 times.
|
76683 | if (present(allow_periodic_degree0)) allow_periodic_degree0_ = allow_periodic_degree0 |
| 193 | |||
| 194 |
1/2✗ Branch 4 → 5 not taken.
✓ Branch 4 → 9 taken 76683 times.
|
76683 | if (bspline%is_scaled) then |
| 195 | write (error_msg, '(A)') "ERROR in BSplineSpace::derivative: the derivative of scaled B-spline space is not implemented. " & | ||
| 196 | ✗ | //"B-spline space must be unscaled." | |
| 197 | ✗ | error stop error_msg | |
| 198 | end if | ||
| 199 | |||
| 200 | 76683 | derivative_degree = bspline%degree - 1 | |
| 201 |
4/4✓ Branch 9 → 10 taken 1770 times.
✓ Branch 9 → 12 taken 74913 times.
✓ Branch 10 → 11 taken 1135 times.
✓ Branch 10 → 12 taken 635 times.
|
76683 | if (allow_periodic_degree0_ .and. bspline%is_periodic) derivative_degree = max(derivative_degree, 0) |
| 202 | 76683 | dbspline = BSplineSpace(bspline%nr_intervals, derivative_degree, is_scaled=.true., is_periodic=bspline%is_periodic) | |
| 203 | 76683 | end function derivative | |
| 204 | |||
| 205 | !> @brief Returns the B-spline space containing the anti-derivatives of the given B-spline space | ||
| 206 | !> | ||
| 207 | !> @param[in] bspline The B-spline space | ||
| 208 | !> | ||
| 209 | !> @return The B-spline space containing the anti-derivatives of the given B-spline space | ||
| 210 | !> | ||
| 211 | !> @note This requires the B-spline space to be scaled | ||
| 212 | 375 | pure type(BSplineSpace) function anti_derivative(bspline) result(dbspline) | |
| 213 | implicit none | ||
| 214 | |||
| 215 | type(BSplineSpace), intent(in) :: bspline | ||
| 216 | character(len=256) :: error_msg | ||
| 217 | |||
| 218 |
1/2✗ Branch 2 → 3 not taken.
✓ Branch 2 → 7 taken 375 times.
|
375 | if (.not. bspline%is_scaled) then |
| 219 | write (error_msg, '(A)') "ERROR in BSplineSpace::anti_derivative: the anti-derivative of unscaled B-spline space is not"& | ||
| 220 | ✗ | & //" implemented. B-spline space must be scaled." | |
| 221 | ✗ | error stop error_msg | |
| 222 | end if | ||
| 223 | |||
| 224 | 375 | dbspline = BSplineSpace(bspline%nr_intervals, bspline%degree + 1, is_scaled=.false., is_periodic=bspline%is_periodic) | |
| 225 | 375 | end function anti_derivative | |
| 226 | |||
| 227 | !> @brief Determine the number of B-splines in the space | ||
| 228 | !> | ||
| 229 | !> @param[in] bspline The B-spline space | ||
| 230 | !> | ||
| 231 | !> @return The number of B-splines | ||
| 232 | 1700121 | pure integer function size_1d(bspline) | |
| 233 | implicit none | ||
| 234 | type(BSplineSpace), intent(in) :: bspline | ||
| 235 | 1715769 | size_1d = bspline%nr_bsplines | |
| 236 | 1700121 | end function | |
| 237 | |||
| 238 | !> @brief Evaluate the j-th B-spline of the B-spline space 'this' at x | ||
| 239 | !> | ||
| 240 | !> @param[in] this The B-spline space | ||
| 241 | !> @param[in] x The evaluation point | ||
| 242 | !> @param[in] j The B-spline index | ||
| 243 | !> @param[in] derivative _(optional)_ If true, the derivative of the B-spline is evaluated | ||
| 244 | !> | ||
| 245 | !> @return The value of the B-spline at x | ||
| 246 | 182387748 | pure real(wp) function evaluate_single(this, x, j, derivative) result(ans) | |
| 247 | use m_bspline_functions | ||
| 248 | |||
| 249 | implicit none | ||
| 250 | |||
| 251 | type(BSplineSpace), intent(in) :: this | ||
| 252 | real(wp), intent(in) :: x | ||
| 253 | integer, intent(in) :: j | ||
| 254 | logical, intent(in), optional :: derivative | ||
| 255 | |||
| 256 | integer :: j_minus_i, j_actual, minus_x | ||
| 257 | real(wp) :: x_eval | ||
| 258 | logical :: derivative_ | ||
| 259 | |||
| 260 |
2/2✓ Branch 2 → 3 taken 1241240 times.
✓ Branch 2 → 4 taken 89952634 times.
|
91193874 | if (present(derivative)) then |
| 261 | 1241240 | derivative_ = derivative | |
| 262 | else | ||
| 263 | derivative_ = .false. | ||
| 264 | end if | ||
| 265 | |||
| 266 | 91193874 | call determine_evaluation_index(this, x, j, x_eval, j_minus_i, j_actual, minus_x=minus_x) | |
| 267 | |||
| 268 |
4/4✓ Branch 5 → 6 taken 86312762 times.
✓ Branch 5 → 12 taken 4881112 times.
✓ Branch 6 → 7 taken 59580230 times.
✓ Branch 6 → 12 taken 26732532 times.
|
91193874 | if (this%degree < j_minus_i .or. j_minus_i < 0) then |
| 269 | ans = 0._wp | ||
| 270 | else | ||
| 271 |
2/2✓ Branch 7 → 8 taken 364000 times.
✓ Branch 7 → 9 taken 59216230 times.
|
59580230 | if (derivative_) then |
| 272 | 364000 | ans = bspline_eval_derivative(x_eval, j_actual, j_minus_i, this%degree, this%is_scaled) | |
| 273 | 364000 | ans = ans * minus_x * this%nr_intervals | |
| 274 | else | ||
| 275 | 59216230 | ans = bspline_eval(x_eval, j_actual, j_minus_i, this%degree, this%is_scaled) | |
| 276 | end if | ||
| 277 |
2/2✓ Branch 10 → 11 taken 343200 times.
✓ Branch 10 → 12 taken 59237030 times.
|
59580230 | if (this%is_scaled) then |
| 278 | 343200 | ans = ans * this%nr_intervals | |
| 279 | end if | ||
| 280 | end if | ||
| 281 | 91193874 | end function | |
| 282 | |||
| 283 | !> @brief Find a root of the B-spline function 'this' minus rhs using Newton's method | ||
| 284 | !> | ||
| 285 | !> @param[in] this The B-spline function | ||
| 286 | !> @param[in] x0 _(optional)_ The initial guess for the root (default: 0.5) | ||
| 287 | !> @param[in] rhs _(optional)_ The right-hand side value to find the root of (default: 0) | ||
| 288 | !> @param[in] max_iter _(optional)_ The maximum number of iterations for the root-finding method (default: 20) | ||
| 289 | !> @param[in] tol _(optional)_ The tolerance for convergence of the root-finding method (default: 1e-8) | ||
| 290 | !> @param[out] success _(optional)_ A logical flag indicating whether the root-finding method converged successfully | ||
| 291 | !> | ||
| 292 | !> @return The root of the B-spline function 'this' minus rhs | ||
| 293 | 285 | real(wp) function find_root_newton(this, x0, rhs, max_iter, tol, success) result(ans) | |
| 294 | implicit none | ||
| 295 | |||
| 296 | type(BSplineFun), intent(in) :: this | ||
| 297 | real(wp), optional, intent(in) :: x0 | ||
| 298 | real(wp), optional, intent(in) :: rhs | ||
| 299 | integer, optional, intent(in) :: max_iter | ||
| 300 | real(wp), optional, intent(in) :: tol | ||
| 301 | logical, optional, intent(out) :: success | ||
| 302 | |||
| 303 | real(wp) :: rhs_, x0_, tol_, error_est, f, df, update | ||
| 304 | integer :: max_iter_, iter | ||
| 305 | logical :: success_ | ||
| 306 | |||
| 307 | rhs_ = 0.0_wp | ||
| 308 |
1/2✓ Branch 2 → 3 taken 285 times.
✗ Branch 2 → 4 not taken.
|
285 | if (present(rhs)) rhs_ = rhs |
| 309 | |||
| 310 | x0_ = 0.5_wp | ||
| 311 |
1/2✗ Branch 4 → 5 not taken.
✓ Branch 4 → 6 taken 285 times.
|
285 | if (present(x0)) x0_ = x0 |
| 312 | |||
| 313 | max_iter_ = 20 | ||
| 314 |
1/2✗ Branch 6 → 7 not taken.
✓ Branch 6 → 8 taken 285 times.
|
285 | if (present(max_iter)) max_iter_ = max_iter |
| 315 | |||
| 316 | tol_ = 1e-6_wp | ||
| 317 |
1/2✓ Branch 8 → 9 taken 285 times.
✗ Branch 8 → 10 not taken.
|
285 | if (present(tol)) tol_ = tol |
| 318 | |||
| 319 |
2/2✓ Branch 10 → 11 taken 1 time.
✓ Branch 10 → 12 taken 284 times.
|
285 | if (this%bspline%degree == 0) max_iter_ = 0 |
| 320 | |||
| 321 | 285 | ans = x0_ | |
| 322 | |||
| 323 | success_ = .false. | ||
| 324 |
2/2✓ Branch 13 → 14 taken 1233 times.
✓ Branch 13 → 30 taken 2 times.
|
1235 | do iter = 1, max_iter_ |
| 325 | 1233 | f = evaluate(this, ans) - rhs_ | |
| 326 |
2/2✓ Branch 14 → 15 taken 240 times.
✓ Branch 14 → 16 taken 993 times.
|
1233 | if (abs(f) < tol_) then |
| 327 | success_ = .true. | ||
| 328 | exit | ||
| 329 | end if | ||
| 330 | |||
| 331 | 993 | df = evaluate(this, ans, derivative=.true.) | |
| 332 |
2/2✓ Branch 16 → 17 taken 43 times.
✓ Branch 16 → 18 taken 950 times.
|
993 | if (abs(df) < 1000 * epsilon(1._wp)) exit |
| 333 | 950 | update = f / df | |
| 334 | 950 | ans = ans - update | |
| 335 | |||
| 336 |
3/4✓ Branch 18 → 19 taken 950 times.
✗ Branch 18 → 20 not taken.
✓ Branch 19 → 20 taken 70 times.
✓ Branch 19 → 28 taken 880 times.
|
950 | if (ans < 0._wp .or. ans > 1._wp) then |
| 337 |
1/2✗ Branch 20 → 21 not taken.
✓ Branch 20 → 25 taken 70 times.
|
70 | if (this%bspline%is_periodic) then |
| 338 | ✗ | ans = modulo(ans, 1._wp) | |
| 339 |
1/2✗ Branch 25 → 26 not taken.
✓ Branch 25 → 27 taken 70 times.
|
70 | elseif (ans < 0._wp) then |
| 340 | ✗ | ans = 0._wp | |
| 341 | else | ||
| 342 | 70 | ans = 1._wp | |
| 343 | end if | ||
| 344 | end if | ||
| 345 | |||
| 346 |
1/2✗ Branch 28 → 15 not taken.
✓ Branch 28 → 29 taken 950 times.
|
995 | if (abs(update) < tol_) then |
| 347 | success_ = .true. | ||
| 348 | exit | ||
| 349 | end if | ||
| 350 | end do | ||
| 351 | |||
| 352 |
1/2✓ Branch 31 → 32 taken 285 times.
✗ Branch 31 → 33 not taken.
|
285 | if (present(success)) success = success_ |
| 353 | |||
| 354 | 285 | end function | |
| 355 | |||
| 356 | !> @brief Find a root of the B-spline function 'this' minus rhs using Brent's method | ||
| 357 | !> | ||
| 358 | !> @param[in] this The B-spline function | ||
| 359 | !> @param[in] x0 The left-hand side of the bracketed interval for the root | ||
| 360 | !> @param[in] x1 The right-hand side of the bracketed interval for the root | ||
| 361 | !> @param[in] rhs _(optional)_ The right-hand side value to find the root of (default: 0) | ||
| 362 | !> @param[in] max_iter _(optional)_ The maximum number of iterations for the root-finding method (default: 20) | ||
| 363 | !> @param[in] tol _(optional)_ The tolerance for convergence of the root-finding method (default: 1e-6) | ||
| 364 | !> @param[out] success _(optional)_ A logical flag indicating whether the root-finding method converged successfully | ||
| 365 | !> @param[in] res0 _(optional)_ The value of the function minus rhs at x0 (if cheaply available) | ||
| 366 | !> @param[in] res1 _(optional)_ The value of the function minus rhs at x1 (if cheaply available) | ||
| 367 | 285 | real(wp) function find_root_brent(this, x0, x1, rhs, max_iter, tol, success, res0, res1) result(ans) | |
| 368 | use m_common, only: brent | ||
| 369 | implicit none | ||
| 370 | |||
| 371 | type(BSplineFun), intent(in) :: this | ||
| 372 | real(wp), intent(in) :: x0, x1 | ||
| 373 | real(wp), optional, intent(in) :: rhs | ||
| 374 | integer, optional, intent(in) :: max_iter | ||
| 375 | real(wp), optional, intent(in) :: tol | ||
| 376 | logical, optional, intent(out) :: success | ||
| 377 | real(wp), optional, intent(in) :: res0, res1 | ||
| 378 | |||
| 379 | real(wp) :: rhs_, x0_, tol_, error_est, f, df, update | ||
| 380 | integer :: max_iter_, iter | ||
| 381 | |||
| 382 | 285 | rhs_ = 0.0_wp | |
| 383 |
1/2✓ Branch 2 → 3 taken 285 times.
✗ Branch 2 → 4 not taken.
|
285 | if (present(rhs)) rhs_ = rhs |
| 384 | |||
| 385 | 285 | max_iter_ = 20 | |
| 386 |
1/2✗ Branch 4 → 5 not taken.
✓ Branch 4 → 6 taken 285 times.
|
285 | if (present(max_iter)) max_iter_ = max_iter |
| 387 | |||
| 388 | 285 | tol_ = 1e-6_wp | |
| 389 |
1/2✓ Branch 6 → 7 taken 285 times.
✗ Branch 6 → 8 not taken.
|
285 | if (present(tol)) tol_ = tol |
| 390 | |||
| 391 |
2/2✓ Branch 8 → 9 taken 1 time.
✓ Branch 8 → 10 taken 284 times.
|
285 | if (this%bspline%degree == 0) max_iter_ = 0 |
| 392 | |||
| 393 | 285 | ans = brent(f_brent, x0, x1, tol_, max_iter_, fa_in=res0, fb_in=res1, success=success) | |
| 394 | contains | ||
| 395 | 2297 | pure real(wp) function f_brent(x) result(fx) | |
| 396 | implicit none | ||
| 397 | real(wp), intent(in) :: x | ||
| 398 | 2297 | fx = evaluate(this, x) - rhs_ | |
| 399 | 2297 | end function | |
| 400 | end function | ||
| 401 | |||
| 402 | !> @brief Determine the B-spline index and the interval index | ||
| 403 | !> | ||
| 404 | !> @param[in] this The B-spline space | ||
| 405 | !> @param[in] x The evaluation point | ||
| 406 | !> @param[in] j The B-spline index (0:nr_bsplines-1) | ||
| 407 | !> @param[out] x_eval The evaluation point in the interval [0, 1) | ||
| 408 | !> @param[out] j_minus_i The difference between the B-spline index and the interval index | ||
| 409 | !> @param[out] j_actual The internal B-spline index (0:degree) | ||
| 410 | !> @param[out] minus_x _(optional)_ Equals -1 if clamped and the j-th B-spline is on the right boundary, equals 1 otherwise | ||
| 411 | 91193874 | pure subroutine determine_evaluation_index(this, x, j, x_eval, j_minus_i, j_actual, minus_x) | |
| 412 | implicit none | ||
| 413 | |||
| 414 | type(BSplineSpace), intent(in) :: this | ||
| 415 | real(wp), intent(in) :: x | ||
| 416 | integer, intent(in) :: j | ||
| 417 | |||
| 418 | real(wp), intent(out) :: x_eval | ||
| 419 | integer, intent(out) :: j_minus_i, j_actual | ||
| 420 | integer, intent(out), optional :: minus_x | ||
| 421 | |||
| 422 | integer :: i_x | ||
| 423 | |||
| 424 | 91193874 | i_x = int(x * this%nr_intervals) | |
| 425 | |||
| 426 |
1/2✓ Branch 2 → 3 taken 91193874 times.
✗ Branch 2 → 4 not taken.
|
91193874 | if (present(minus_x)) then |
| 427 | 91193874 | minus_x = 1 | |
| 428 | end if | ||
| 429 | |||
| 430 |
2/2✓ Branch 4 → 5 taken 8082776 times.
✓ Branch 4 → 9 taken 83111098 times.
|
91193874 | if (this%is_periodic) then |
| 431 | 8082776 | x_eval = x * this%nr_intervals - i_x | |
| 432 | 8082776 | i_x = mod(i_x, this%nr_intervals) | |
| 433 | |||
| 434 | ! All of the B-splines are the same | ||
| 435 | 8082776 | j_actual = this%degree | |
| 436 |
4/4✓ Branch 5 → 6 taken 2912906 times.
✓ Branch 5 → 8 taken 5169870 times.
✓ Branch 6 → 7 taken 1415884 times.
✓ Branch 6 → 8 taken 1497022 times.
|
8082776 | if (i_x > j .and. j < this%degree) then |
| 437 | ! Impose continuity/periodicity if x is on the right side of the domain, | ||
| 438 | ! whereas j refers to a left boundary B-spline, which due to periodicity overlaps | ||
| 439 | ! with a right boundary B-spline | ||
| 440 | 1415884 | j_minus_i = this%nr_bsplines + j - i_x | |
| 441 | else | ||
| 442 | 6666892 | j_minus_i = j - i_x | |
| 443 | end if | ||
| 444 | |||
| 445 | else | ||
| 446 | 83111098 | i_x = max(0, min(i_x, this%nr_intervals - 1)) | |
| 447 | 83111098 | j_minus_i = j - i_x | |
| 448 | |||
| 449 |
2/2✓ Branch 9 → 10 taken 28102174 times.
✓ Branch 9 → 11 taken 55008924 times.
|
83111098 | if (j < this%degree) then |
| 450 | 28102174 | j_actual = j | |
| 451 | 28102174 | x_eval = x * this%nr_intervals - i_x | |
| 452 |
2/2✓ Branch 11 → 12 taken 18975800 times.
✓ Branch 11 → 13 taken 36033124 times.
|
55008924 | else if (j < this%nr_intervals) then |
| 453 | ! The {1+degree, ..., nr_intervals} B-splines are all the same | ||
| 454 | 18975800 | j_actual = this%degree | |
| 455 | 18975800 | x_eval = x * this%nr_intervals - i_x | |
| 456 | else | ||
| 457 | 36033124 | j_actual = this%nr_intervals - 1 + this%degree - j | |
| 458 | 36033124 | x_eval = i_x + 1 - x * this%nr_intervals | |
| 459 | 36033124 | j_minus_i = this%degree - j_minus_i | |
| 460 | |||
| 461 |
1/2✓ Branch 13 → 14 taken 36033124 times.
✗ Branch 13 → 15 not taken.
|
36033124 | if (present(minus_x)) then |
| 462 | 36033124 | minus_x = -1 | |
| 463 | end if | ||
| 464 | end if | ||
| 465 | end if | ||
| 466 | |||
| 467 | 91193874 | end subroutine | |
| 468 | |||
| 469 | !> @brief Integrate the j-th B-spline of the B-spline space 'this' using a quadrature rule | ||
| 470 | !> | ||
| 471 | !> @param[in] this The B-spline space | ||
| 472 | !> @param[in] j The B-spline index | ||
| 473 | !> @param[in] userfun _(optional)_ The user-defined function to multiply the B-spline with | ||
| 474 | !> @param[in] n_quad _(optional)_ The number of quadrature points | ||
| 475 | !> | ||
| 476 | !> @return The integral of the B-spline multiplied by the optional user function | ||
| 477 | 208964 | pure real(wp) function quadrature_eval(this, j, userfun, n_quad) result(ans) | |
| 478 | use m_common, only: user_function_1d_interface | ||
| 479 | |||
| 480 | implicit none | ||
| 481 | |||
| 482 | type(BSplineSpace), intent(in) :: this | ||
| 483 | integer, intent(in) :: j | ||
| 484 | procedure(user_function_1d_interface), optional :: userfun | ||
| 485 | integer, intent(in), optional :: n_quad | ||
| 486 | |||
| 487 | integer :: i, n_quad_, i0_min, i0_max, i1_min, i1_max | ||
| 488 |
2/2✓ Branch 2 → 3 taken 73959 times.
✓ Branch 2 → 4 taken 30523 times.
|
104482 | if (present(n_quad)) then |
| 489 | 73959 | n_quad_ = n_quad | |
| 490 | else | ||
| 491 | 30523 | n_quad_ = 1 + int(this%degree / 2) | |
| 492 | end if | ||
| 493 | |||
| 494 | 104482 | call get_index_ranges_single(this, j, i0_min, i0_max, i1_min, i1_max) | |
| 495 | |||
| 496 | ans = 0._wp | ||
| 497 |
2/2✓ Branch 7 → 8 taken 435055 times.
✓ Branch 7 → 12 taken 104482 times.
|
539537 | do i = i0_min, i0_max |
| 498 |
2/2✓ Branch 8 → 9 taken 42690 times.
✓ Branch 8 → 10 taken 392365 times.
|
539537 | if (present(userfun)) then |
| 499 | 42690 | ans = ans + quadrature_interval(this, j, i, n_quad=n_quad_, userfun=userfun) | |
| 500 | else | ||
| 501 | 392365 | ans = ans + quadrature_interval(this, j, i, n_quad=n_quad_) | |
| 502 | end if | ||
| 503 | end do | ||
| 504 |
2/2✓ Branch 14 → 15 taken 55881 times.
✓ Branch 14 → 19 taken 104482 times.
|
160363 | do i = i1_min, i1_max |
| 505 |
2/2✓ Branch 15 → 16 taken 3090 times.
✓ Branch 15 → 17 taken 52791 times.
|
160363 | if (present(userfun)) then |
| 506 | 3090 | ans = ans + quadrature_interval(this, j, i, n_quad=n_quad_, userfun=userfun) | |
| 507 | else | ||
| 508 | 52791 | ans = ans + quadrature_interval(this, j, i, n_quad=n_quad_) | |
| 509 | end if | ||
| 510 | end do | ||
| 511 | |||
| 512 | 104482 | end function | |
| 513 | |||
| 514 | !> @brief Determine j_actual - i, where the input j may exceed the number of B-splines | ||
| 515 | !> | ||
| 516 | !> @param[in] this The B-spline basis | ||
| 517 | !> @param[in] j The B-spline index | ||
| 518 | !> @param[in] i The interval index | ||
| 519 | !> | ||
| 520 | !> @return The j_actual - i | ||
| 521 | 44403521 | pure integer function get_j_minus_i_quad(this, j, i) result(j_minus_i) | |
| 522 | implicit none | ||
| 523 | |||
| 524 | type(BSplineSpace), intent(in) :: this | ||
| 525 | integer, intent(in) :: j, i | ||
| 526 | |||
| 527 | 44403521 | j_minus_i = j - i | |
| 528 |
2/2✓ Branch 2 → 3 taken 11075104 times.
✓ Branch 2 → 8 taken 33328417 times.
|
44403521 | if (this%is_periodic) then |
| 529 |
2/2✓ Branch 3 → 4 taken 860508 times.
✓ Branch 3 → 5 taken 10214596 times.
|
11075104 | if (j >= this%nr_bsplines) then |
| 530 | 860508 | j_minus_i = j_minus_i - this%nr_bsplines | |
| 531 |
4/4✓ Branch 5 → 6 taken 3085735 times.
✓ Branch 5 → 8 taken 7128861 times.
✓ Branch 6 → 7 taken 2689933 times.
✓ Branch 6 → 8 taken 395802 times.
|
10214596 | else if (i > j .and. j < this%degree) then |
| 532 | ! Impose continuity/periodicity if x is on the right side of the domain, | ||
| 533 | ! whereas j refers to a left boundary B-spline, which due to periodicity overlaps | ||
| 534 | ! with a right boundary B-spline | ||
| 535 | 2689933 | j_minus_i = j_minus_i + this%nr_bsplines | |
| 536 | end if | ||
| 537 | end if | ||
| 538 | |||
| 539 | 44403521 | end function | |
| 540 | |||
| 541 | !> @brief Integrate the j-th B-spline of the B-spline space 'this' over the i-th interval using a quadrature rule | ||
| 542 | !> | ||
| 543 | !> @param[in] this The B-spline space | ||
| 544 | !> @param[in] j The B-spline index | ||
| 545 | !> @param[in] i The interval index | ||
| 546 | !> @param[in] userfun _(optional)_ The user-defined function to multiply the B-spline with | ||
| 547 | !> @param[in] n_quad _(optional)_ The number of quadrature points | ||
| 548 | !> | ||
| 549 | !> @return The integral of the B-spline multiplied by the optional userfun over the interval | ||
| 550 | 1390733 | pure real(wp) function quadrature_interval(this, j, i, userfun, n_quad) result(ans) | |
| 551 | use m_common, only: user_function_1d_interface | ||
| 552 | use m_quadrature, only: quadrature | ||
| 553 | |||
| 554 | implicit none | ||
| 555 | |||
| 556 | type(BSplineSpace), intent(in) :: this | ||
| 557 | integer, intent(in) :: j, i | ||
| 558 | procedure(user_function_1d_interface), optional :: userfun | ||
| 559 | integer, intent(in), optional :: n_quad | ||
| 560 | |||
| 561 | integer :: j_minus_i, n_quad_ | ||
| 562 | |||
| 563 |
2/2✓ Branch 2 → 3 taken 490936 times.
✓ Branch 2 → 4 taken 899797 times.
|
1390733 | if (present(n_quad)) then |
| 564 | 490936 | n_quad_ = n_quad | |
| 565 | else | ||
| 566 | 899797 | n_quad_ = 1 + int(this%degree / 2) | |
| 567 | end if | ||
| 568 | |||
| 569 | 1390733 | j_minus_i = get_j_minus_i_quad(this, j, i) | |
| 570 | |||
| 571 |
4/4✓ Branch 5 → 6 taken 962834 times.
✓ Branch 5 → 8 taken 427899 times.
✓ Branch 6 → 7 taken 533828 times.
✓ Branch 6 → 8 taken 429006 times.
|
1390733 | if (j_minus_i > this%degree .or. j_minus_i < 0) then |
| 572 | ans = 0._wp | ||
| 573 | else | ||
| 574 | 533828 | ans = quadrature(n_quad_, evaluate_jth_bspline, i * this%h, (i + 1) * this%h) | |
| 575 | end if | ||
| 576 | contains | ||
| 577 | 2043906 | pure real(wp) function evaluate_jth_bspline(x) result(f) | |
| 578 | implicit none | ||
| 579 | |||
| 580 | real(wp), intent(in) :: x | ||
| 581 | |||
| 582 | ! TODO efficiency | ||
| 583 | 2043906 | f = evaluate(this, x, j) | |
| 584 |
2/2✓ Branch 2 → 3 taken 312992 times.
✓ Branch 2 → 5 taken 1730914 times.
|
2043906 | if (present(userfun)) then |
| 585 | 312992 | f = f * userfun(x) | |
| 586 | end if | ||
| 587 | 2043906 | end function | |
| 588 | end function | ||
| 589 | |||
| 590 | !> @brief Get the interval index ranges (i.e., the support) for integrating the j-th B-spline | ||
| 591 | !> | ||
| 592 | !> @param[in] this The B-spline space | ||
| 593 | !> @param[in] j The B-spline index | ||
| 594 | !> @param[out] i0_min The minimum interval index for the first range | ||
| 595 | !> @param[out] i0_max The maximum interval index for the first range | ||
| 596 | !> @param[out] i1_min The minimum interval index for the second range | ||
| 597 | !> @param[out] i1_max The maximum interval index for the second range | ||
| 598 | 336571 | pure subroutine get_index_ranges_single(this, j, i0_min, i0_max, i1_min, i1_max) | |
| 599 | implicit none | ||
| 600 | |||
| 601 | type(BSplineSpace), intent(in) :: this | ||
| 602 | |||
| 603 | integer, intent(in) :: j | ||
| 604 | integer, intent(out) :: i0_min, i0_max, i1_min, i1_max | ||
| 605 | |||
| 606 |
4/4✓ Branch 2 → 3 taken 134427 times.
✓ Branch 2 → 4 taken 202144 times.
✓ Branch 3 → 4 taken 58385 times.
✓ Branch 3 → 5 taken 76042 times.
|
336571 | if (.not. this%is_periodic .or. j > this%degree) then |
| 607 | 260529 | i0_min = max(0, j - this%degree) | |
| 608 | 260529 | i0_max = min(this%nr_intervals - 1, j) | |
| 609 | 260529 | i1_min = 1 | |
| 610 | 260529 | i1_max = -1 | |
| 611 | else ! if (j <= this%degree) then | ||
| 612 | 76042 | i0_min = 0 | |
| 613 | 76042 | i0_max = j | |
| 614 | 76042 | i1_min = this%nr_intervals - this%degree + j | |
| 615 | 76042 | i1_max = this%nr_intervals - 1 | |
| 616 | end if | ||
| 617 | 336571 | end subroutine get_index_ranges_single | |
| 618 | |||
| 619 | !> @brief Get the interval index ranges (i.e., the support) for integrating the product of the j1-th and j2-th B-splines | ||
| 620 | !> | ||
| 621 | !> @param[in] this1 The first B-spline space | ||
| 622 | !> @param[in] this2 The second B-spline space | ||
| 623 | !> @param[in] j1 The first B-spline index | ||
| 624 | !> @param[in] j2 The second B-spline index | ||
| 625 | !> @param[out] i0_min The minimum interval index for the first range | ||
| 626 | !> @param[out] i0_max The maximum interval index for the first range | ||
| 627 | !> @param[out] i1_min The minimum interval index for the second range | ||
| 628 | !> @param[out] i1_max The maximum interval index for the second range | ||
| 629 | 9197893 | pure subroutine get_index_ranges_product(this1, this2, j1, j2, i0_min, i0_max, i1_min, i1_max) | |
| 630 | implicit none | ||
| 631 | |||
| 632 | type(BSplineSpace), intent(in) :: this1, this2 | ||
| 633 | |||
| 634 | integer, intent(in) :: j1, j2 | ||
| 635 | integer, intent(out) :: i0_min, i0_max, i1_min, i1_max | ||
| 636 | |||
| 637 |
6/6✓ Branch 2 → 3 taken 779110 times.
✓ Branch 2 → 5 taken 8418783 times.
✓ Branch 3 → 4 taken 750758 times.
✓ Branch 3 → 6 taken 28352 times.
✓ Branch 4 → 5 taken 732783 times.
✓ Branch 4 → 6 taken 17975 times.
|
9197893 | if (.not. this1%is_periodic .or. (j1 >= this1%degree .and. j2 >= this2%degree)) then |
| 638 | 9151566 | i0_min = max(0, j1 - this1%degree, j2 - this2%degree) | |
| 639 | 9151566 | i0_max = min(this1%nr_intervals - 1, j1, j2) | |
| 640 | 9151566 | i1_min = 1 | |
| 641 | 9151566 | i1_max = -1 | |
| 642 |
3/4✓ Branch 6 → 7 taken 17975 times.
✓ Branch 6 → 9 taken 28352 times.
✓ Branch 7 → 8 taken 17975 times.
✗ Branch 7 → 9 not taken.
|
46327 | else if (j1 >= this1%degree .and. j2 < this2%degree) then |
| 643 | ! j1: j1-degree,j1 | ||
| 644 | ! j2: 1,j2 | ||
| 645 | 17975 | i0_min = max(j1 - this1%degree, 0) | |
| 646 | 17975 | i0_max = min(j1, j2) | |
| 647 | ! j1: j1-degree,j1 | ||
| 648 | ! j2: nr_intervals-(degree-j2),nr_intervals | ||
| 649 | 17975 | i1_min = max(j1 - this1%degree, this1%nr_intervals - this2%degree + j2) | |
| 650 | 17975 | i1_max = min(j1, this1%nr_intervals - 1) | |
| 651 |
3/4✓ Branch 9 → 10 taken 28352 times.
✗ Branch 9 → 12 not taken.
✓ Branch 10 → 11 taken 23022 times.
✓ Branch 10 → 12 taken 5330 times.
|
28352 | else if (j1 < this1%degree .and. j2 >= this2%degree) then |
| 652 | ! j1: 1,j1 | ||
| 653 | ! j2: j2-degree,j2 | ||
| 654 | 23022 | i0_min = max(j2 - this2%degree, 0) | |
| 655 | 23022 | i0_max = min(j1, j2) | |
| 656 | ! j1: nr_intervals-(degree-j1),nr_intervals | ||
| 657 | ! j2: j2-degree,j2 | ||
| 658 | 23022 | i1_min = max(j2 - this2%degree, this1%nr_intervals - this1%degree + j1) | |
| 659 | 23022 | i1_max = min(j2, this1%nr_intervals - 1) | |
| 660 | else ! if (j1 < this%degree .and. j2 < this%degree) then | ||
| 661 | ! j1: 1,j1 | ||
| 662 | ! j2: 1,j2 | ||
| 663 | 5330 | i0_min = 0 | |
| 664 | 5330 | i0_max = max(j1, j2) | |
| 665 | ! j1: nr_intervals-(degree-j1),nr_intervals | ||
| 666 | ! j2: nr_intervals-(degree-j2),nr_intervals | ||
| 667 | 5330 | i1_min = max(this1%nr_intervals - this1%degree + j1, this1%nr_intervals - this2%degree + j2) | |
| 668 | 5330 | i1_max = this1%nr_intervals - 1 | |
| 669 | end if | ||
| 670 | 9197893 | end subroutine get_index_ranges_product | |
| 671 | |||
| 672 | !> @brief Integrate the product of the j1-th and j2-th B-splines of the B-spline spaces 'this1' and 'this2' using a quadrature | ||
| 673 | !> rule | ||
| 674 | !> | ||
| 675 | !> @param[in] this1 The first B-spline basis | ||
| 676 | !> @param[in] this2 The second B-spline basis | ||
| 677 | !> @param[in] j1 The first B-spline index | ||
| 678 | !> @param[in] j2 The second B-spline index | ||
| 679 | !> @param[in] userfun _(optional)_ The user-defined function to multiply the B-splines with | ||
| 680 | !> @param[in] n_quad _(optional)_ The number of quadrature points | ||
| 681 | !> | ||
| 682 | !> @return The integral of the product of the B-splines multiplied by the optional userfun | ||
| 683 | !> | ||
| 684 | !> @note It is assumed that the two B-spline spaces have the same number of intervals as well as periodicity. The degree of the | ||
| 685 | !> B-splines may differ. | ||
| 686 | 6129440 | pure real(wp) function quadrature_product_eval(this1, this2, j1, j2, userfun, n_quad) result(ans) | |
| 687 | use m_common, only: user_function_1d_interface | ||
| 688 | |||
| 689 | implicit none | ||
| 690 | |||
| 691 | type(BSplineSpace), intent(in) :: this1, this2 | ||
| 692 | integer, intent(in) :: j1, j2 | ||
| 693 | procedure(user_function_1d_interface), optional :: userfun | ||
| 694 | integer, intent(in), optional :: n_quad | ||
| 695 | |||
| 696 | integer :: i, n_quad_, i0_min, i0_max, i1_min, i1_max | ||
| 697 | |||
| 698 |
2/2✓ Branch 2 → 3 taken 1996335 times.
✓ Branch 2 → 4 taken 1068385 times.
|
3064720 | if (present(n_quad)) then |
| 699 | 1996335 | n_quad_ = n_quad | |
| 700 | else | ||
| 701 | 1068385 | n_quad_ = 1 + max(this1%degree, this2%degree) | |
| 702 | end if | ||
| 703 | |||
| 704 | 3064720 | call get_index_ranges_product(this1, this2, j1, j2, i0_min, i0_max, i1_min, i1_max) | |
| 705 | |||
| 706 | ans = 0._wp | ||
| 707 |
2/2✓ Branch 7 → 8 taken 3178102 times.
✓ Branch 7 → 12 taken 3064720 times.
|
6242822 | do i = i0_min, i0_max |
| 708 |
1/2✗ Branch 8 → 9 not taken.
✓ Branch 8 → 10 taken 3178102 times.
|
6242822 | if (present(userfun)) then |
| 709 | ✗ | ans = ans + quadrature_product_interval(this1, this2, j1, j2, i, n_quad=n_quad_, userfun=userfun) | |
| 710 | else | ||
| 711 | 3178102 | ans = ans + quadrature_product_interval(this1, this2, j1, j2, i, n_quad=n_quad_) | |
| 712 | end if | ||
| 713 | end do | ||
| 714 |
2/2✓ Branch 14 → 15 taken 29906 times.
✓ Branch 14 → 19 taken 3064720 times.
|
3094626 | do i = i1_min, i1_max |
| 715 |
1/2✗ Branch 15 → 16 not taken.
✓ Branch 15 → 17 taken 29906 times.
|
3094626 | if (present(userfun)) then |
| 716 | ✗ | ans = ans + quadrature_product_interval(this1, this2, j1, j2, i, n_quad=n_quad_, userfun=userfun) | |
| 717 | else | ||
| 718 | 29906 | ans = ans + quadrature_product_interval(this1, this2, j1, j2, i, n_quad=n_quad_) | |
| 719 | end if | ||
| 720 | end do | ||
| 721 | 3064720 | end function | |
| 722 | |||
| 723 | !> @brief Integrate the product of the j1-th and j2-th B-splines of the B-spline spaces 'this1' and 'this2' over the i-th interval | ||
| 724 | !> using a quadrature rule | ||
| 725 | !> | ||
| 726 | !> @param[in] this1 The first B-spline basis | ||
| 727 | !> @param[in] this2 The second B-spline basis | ||
| 728 | !> @param[in] j1 The first B-spline index | ||
| 729 | !> @param[in] j2 The second B-spline index | ||
| 730 | !> @param[in] i The interval index | ||
| 731 | !> @param[in] userfun _(optional)_ The user-defined function to multiply the B-splines with | ||
| 732 | !> @param[in] n_quad _(optional)_ The number of quadrature points | ||
| 733 | !> | ||
| 734 | !> @return The integral of the product of the B-splines multiplied by the optional userfun over the interval | ||
| 735 | !> | ||
| 736 | !> @note It is assumed that the two B-spline spaces have the same number of intervals as well as periodicity. The degree of the | ||
| 737 | !> B-splines may differ. | ||
| 738 | 3208008 | pure real(wp) function quadrature_product_interval(this1, this2, j1, j2, i, userfun, n_quad) result(ans) | |
| 739 | use m_common, only: user_function_1d_interface | ||
| 740 | use m_quadrature, only: quad => quadrature | ||
| 741 | |||
| 742 | implicit none | ||
| 743 | |||
| 744 | type(BSplineSpace), intent(in) :: this1, this2 | ||
| 745 | integer, intent(in) :: j1, j2, i | ||
| 746 | procedure(user_function_1d_interface), optional :: userfun | ||
| 747 | integer, intent(in), optional :: n_quad | ||
| 748 | |||
| 749 | integer :: j1_minus_i, j2_minus_i, n_quad_ | ||
| 750 | |||
| 751 |
1/2✓ Branch 2 → 3 taken 3208008 times.
✗ Branch 2 → 4 not taken.
|
3208008 | if (present(n_quad)) then |
| 752 | 3208008 | n_quad_ = n_quad | |
| 753 | else | ||
| 754 | ✗ | n_quad_ = 1 + max(this1%degree, this2%degree) | |
| 755 | end if | ||
| 756 | |||
| 757 | 3208008 | j1_minus_i = get_j_minus_i_quad(this1, j1, i) | |
| 758 | 3208008 | j2_minus_i = get_j_minus_i_quad(this2, j2, i) | |
| 759 | |||
| 760 |
6/8✓ Branch 5 → 6 taken 3203178 times.
✓ Branch 5 → 10 taken 4830 times.
✓ Branch 6 → 7 taken 3203178 times.
✗ Branch 6 → 10 not taken.
✓ Branch 7 → 8 taken 3198348 times.
✓ Branch 7 → 10 taken 4830 times.
✓ Branch 8 → 9 taken 3198348 times.
✗ Branch 8 → 10 not taken.
|
3208008 | if (j1_minus_i > this1%degree .or. j1_minus_i < 0 .or. j2_minus_i > this2%degree .or. j2_minus_i < 0) then |
| 761 | ans = 0._wp | ||
| 762 | else | ||
| 763 | 3198348 | ans = quad(n_quad_, evaluate_bspline_product, i * this1%h, (i + 1) * this1%h) | |
| 764 | end if | ||
| 765 | contains | ||
| 766 | 16295178 | pure real(wp) function evaluate_bspline_product(x) result(f) | |
| 767 | implicit none | ||
| 768 | |||
| 769 | real(wp), intent(in) :: x | ||
| 770 | |||
| 771 | ! TODO efficiency: internal evaluation function with internal indices and no bounds checking | ||
| 772 | 16295178 | f = evaluate(this1, x, j1) * evaluate(this2, x, j2) | |
| 773 |
1/2✗ Branch 2 → 3 not taken.
✓ Branch 2 → 5 taken 16295178 times.
|
16295178 | if (present(userfun)) then |
| 774 | ✗ | f = f * userfun(x) | |
| 775 | end if | ||
| 776 | 16295178 | end function | |
| 777 | end function | ||
| 778 | |||
| 779 | !> @brief Get the internal B-spline index and whether the B-spline is mirrored | ||
| 780 | !> | ||
| 781 | !> @param[in] this The B-spline space | ||
| 782 | !> @param[in] j The B-spline index | ||
| 783 | !> @param[out] j_internal The internal B-spline index (0:degree) | ||
| 784 | !> @param[out] bspline_mirrored _(optional)_ Whether the B-spline is mirrored (i.e., the B-spline is on the right boundary of the | ||
| 785 | !> domain) | ||
| 786 | ✗ | pure subroutine get_j_internal(this, j, j_internal, bspline_mirrored) | |
| 787 | !$acc routine seq | ||
| 788 | implicit none | ||
| 789 | |||
| 790 | type(BSplineSpace), intent(in) :: this | ||
| 791 | integer, intent(in) :: j | ||
| 792 | integer, intent(out) :: j_internal | ||
| 793 | logical, intent(out) :: bspline_mirrored | ||
| 794 | |||
| 795 | ✗ | bspline_mirrored = .false. | |
| 796 | ✗ | if (this%is_periodic) then | |
| 797 | ! All of the B-splines are the same | ||
| 798 | ✗ | j_internal = this%degree | |
| 799 | else | ||
| 800 |
2/2✓ Branch 4 → 5 taken 886412608 times.
✓ Branch 4 → 6 taken 2146419613 times.
|
3032832221 | if (j < this%degree) then |
| 801 | 886412608 | j_internal = j | |
| 802 |
2/2✓ Branch 6 → 7 taken 1162252703 times.
✓ Branch 6 → 8 taken 984166910 times.
|
2146419613 | else if (j < this%nr_intervals) then |
| 803 | ! The {1+degree, ..., nr_intervals} B-splines are all the same | ||
| 804 | 1162252703 | j_internal = this%degree | |
| 805 | else | ||
| 806 | 984166910 | j_internal = this%nr_intervals - 1 + this%degree - j | |
| 807 | 984166910 | bspline_mirrored = .true. | |
| 808 | end if | ||
| 809 | end if | ||
| 810 | ✗ | end subroutine get_j_internal | |
| 811 | |||
| 812 | !> @brief Get the internal B-spline index and the actual B-spline index | ||
| 813 | !> | ||
| 814 | !> @param[in] bspline The B-spline space | ||
| 815 | !> @param[in] i The interval index | ||
| 816 | !> @param[in] j_minus_i The difference between the B-spline index and the interval index | ||
| 817 | !> @param[out] j The actual B-spline index (0:nr_bsplines-1) | ||
| 818 | !> @param[out] j_internal _(optional)_ The internal B-spline index (0:degree) | ||
| 819 | !> @param[out] bspline_mirrored _(optional)_ Whether the B-spline is mirrored (i.e., the B-spline is on the right boundary of the | ||
| 820 | !> domain) | ||
| 821 | ✗ | pure subroutine get_internal_and_actual_inds(bspline, i, j_minus_i, j, j_internal, bspline_mirrored) | |
| 822 | !$acc routine seq | ||
| 823 | implicit none | ||
| 824 | |||
| 825 | type(BSplineSpace), intent(in) :: bspline | ||
| 826 | integer, intent(in) :: i, j_minus_i | ||
| 827 | |||
| 828 | integer, intent(out) :: j | ||
| 829 | integer, intent(out), optional :: j_internal | ||
| 830 | logical, intent(out), optional :: bspline_mirrored | ||
| 831 | |||
| 832 | ✗ | j = j_minus_i + i | |
| 833 | ✗ | if (bspline%is_periodic .and. j >= bspline%nr_intervals) then | |
| 834 | 291046106 | j = j - bspline%nr_intervals | |
| 835 | end if | ||
| 836 | |||
| 837 | ✗ | if (present(j_internal) .and. present(bspline_mirrored)) then | |
| 838 | ✗ | call get_j_internal(bspline, j, j_internal, bspline_mirrored) | |
| 839 | end if | ||
| 840 | ✗ | end subroutine get_internal_and_actual_inds | |
| 841 | |||
| 842 | !> @brief Initialize the B-spline function | ||
| 843 | !> | ||
| 844 | !> @param[out] bfun The B-spline function | ||
| 845 | !> @param[in] bspline The B-spline space | ||
| 846 |
1/2✓ Branch 2 → 3 taken 18551 times.
✗ Branch 2 → 5 not taken.
|
18551 | subroutine init_bspline_fun(bfun, bspline, vec) |
| 847 | implicit none | ||
| 848 | |||
| 849 | class(BSplineFun), intent(out) :: bfun | ||
| 850 | type(BSplineSpace), intent(in) :: bspline | ||
| 851 | Vec, intent(in), optional :: vec | ||
| 852 | |||
| 853 | integer :: ierr, i | ||
| 854 | integer, allocatable :: inds(:) | ||
| 855 | |||
| 856 | 18551 | bfun%bspline = bspline | |
| 857 |
1/2✗ Branch 5 → 6 not taken.
✓ Branch 5 → 7 taken 18551 times.
|
18551 | if (allocated(bfun%data)) deallocate (bfun%data) |
| 858 |
7/14✗ Branch 7 → 8 not taken.
✓ Branch 7 → 9 taken 18551 times.
✓ Branch 9 → 8 taken 18551 times.
✗ Branch 9 → 10 not taken.
✓ Branch 10 → 11 taken 18551 times.
✗ Branch 10 → 12 not taken.
✓ Branch 12 → 13 taken 18551 times.
✗ Branch 12 → 14 not taken.
✗ Branch 14 → 15 not taken.
✓ Branch 14 → 16 taken 18551 times.
✗ Branch 16 → 17 not taken.
✓ Branch 16 → 18 taken 18551 times.
✗ Branch 18 → 19 not taken.
✓ Branch 18 → 20 taken 18551 times.
|
55653 | allocate (bfun%data(0:bspline%nr_intervals - 1 + bspline%degree)) |
| 859 | |||
| 860 |
2/2✓ Branch 20 → 21 taken 7824 times.
✓ Branch 20 → 43 taken 10727 times.
|
18551 | if (present(vec)) then |
| 861 |
6/12✗ Branch 21 → 22 not taken.
✓ Branch 21 → 23 taken 7824 times.
✓ Branch 23 → 22 taken 7824 times.
✗ Branch 23 → 24 not taken.
✓ Branch 24 → 25 taken 7824 times.
✗ Branch 24 → 26 not taken.
✓ Branch 26 → 27 taken 7824 times.
✗ Branch 26 → 28 not taken.
✗ Branch 28 → 29 not taken.
✓ Branch 28 → 30 taken 7824 times.
✗ Branch 30 → 31 not taken.
✓ Branch 30 → 32 taken 7824 times.
|
23472 | allocate (inds(size(bspline))) |
| 862 |
2/2✓ Branch 32 → 33 taken 193143 times.
✓ Branch 32 → 34 taken 7824 times.
|
200967 | do i = 1, size(bspline) |
| 863 | 200967 | inds(i) = i - 1 | |
| 864 | end do | ||
| 865 | |||
| 866 |
1/2✗ Branch 36 → 37 not taken.
✓ Branch 36 → 40 taken 7824 times.
|
7824 | PetscCall(VecGetValues(vec, size(bspline), inds, bfun%data, ierr)) |
| 867 |
1/2✗ Branch 40 → 41 not taken.
✓ Branch 40 → 42 taken 7824 times.
|
7824 | deallocate (inds) |
| 868 | |||
| 869 | 7824 | call bfun%periodicity() | |
| 870 | else | ||
| 871 |
2/2✓ Branch 43 → 44 taken 287578 times.
✓ Branch 43 → 45 taken 10727 times.
|
298305 | bfun%data = 0._wp |
| 872 | end if | ||
| 873 | ✗ | end subroutine init_bspline_fun | |
| 874 | |||
| 875 | !> @brief Enforce periodicity on the B-spline function | ||
| 876 | !> | ||
| 877 | !> @param[inout] bfun The B-spline function | ||
| 878 | 7824 | subroutine bspline_fun_periodicity(bfun) | |
| 879 | implicit none | ||
| 880 | |||
| 881 | class(BSplineFun), intent(inout) :: bfun | ||
| 882 | integer :: i | ||
| 883 | |||
| 884 |
2/2✓ Branch 2 → 3 taken 4225 times.
✓ Branch 2 → 7 taken 3599 times.
|
7824 | if (bfun%bspline%is_periodic) then |
| 885 |
2/2✓ Branch 4 → 5 taken 15219 times.
✓ Branch 4 → 6 taken 4225 times.
|
19444 | do i = 0, bfun%bspline%degree - 1 |
| 886 | 19444 | bfun%data(i + bfun%bspline%nr_intervals) = bfun%data(i) | |
| 887 | end do | ||
| 888 | end if | ||
| 889 | 7824 | end subroutine bspline_fun_periodicity | |
| 890 | |||
| 891 | !> @brief Destroy the B-spline function | ||
| 892 | !> | ||
| 893 | !> @param[inout] bfun The B-spline function | ||
| 894 | 4876 | subroutine destroy_bspline_fun(bfun) | |
| 895 | implicit none | ||
| 896 | |||
| 897 | class(BSplineFun), intent(inout) :: bfun | ||
| 898 | |||
| 899 |
1/2✓ Branch 2 → 3 taken 4876 times.
✗ Branch 2 → 4 not taken.
|
4876 | if (allocated(bfun%data)) deallocate (bfun%data) |
| 900 | 4876 | end subroutine destroy_bspline_fun | |
| 901 | |||
| 902 | ! TODO the 'final' functionality is not reliably supported by all (any?) compilers, so we do not use it here | ||
| 903 | |||
| 904 | ! !> @brief Destroy the B-spline function when it goes out of scope (final) | ||
| 905 | ! !> | ||
| 906 | ! !> @param[inout] bfun The B-spline function | ||
| 907 | ! subroutine destroy_bspline_fun_final(bfun) | ||
| 908 | ! implicit none | ||
| 909 | |||
| 910 | ! type(BSplineFun), intent(inout) :: bfun | ||
| 911 | |||
| 912 | ! call bfun%destroy() | ||
| 913 | ! end subroutine destroy_bspline_fun_final | ||
| 914 | |||
| 915 | !> @brief Get the interval index for a point in the B-spline space | ||
| 916 | !> | ||
| 917 | !> @param[in] bspline The B-spline space | ||
| 918 | !> @param[in] x The point | ||
| 919 | !> | ||
| 920 | !> @return The interval index for the point in the B-spline space (0:nr_intervals-1) | ||
| 921 | 670124591 | pure integer function get_interval_index(bspline, x) result(i_x) | |
| 922 | !$acc routine seq | ||
| 923 | implicit none | ||
| 924 | |||
| 925 | type(BSplineSpace), intent(in) :: bspline | ||
| 926 | real(wp), intent(in) :: x | ||
| 927 | |||
| 928 | 670124591 | i_x = int(x * bspline%nr_intervals) | |
| 929 |
1/2✗ Branch 2 → 3 not taken.
✓ Branch 2 → 5 taken 670124591 times.
|
670124591 | if (i_x < 0) then |
| 930 | ✗ | if (bspline%is_periodic) then | |
| 931 | ✗ | i_x = modulo(i_x, bspline%nr_intervals) | |
| 932 | else | ||
| 933 | i_x = 0 | ||
| 934 | end if | ||
| 935 |
2/2✓ Branch 5 → 6 taken 423 times.
✓ Branch 5 → 9 taken 670124168 times.
|
670124591 | else if (i_x >= bspline%nr_intervals) then |
| 936 |
1/2✗ Branch 6 → 7 not taken.
✓ Branch 6 → 8 taken 423 times.
|
423 | if (bspline%is_periodic) then |
| 937 | ✗ | i_x = modulo(i_x, bspline%nr_intervals) | |
| 938 | else | ||
| 939 | 423 | i_x = bspline%nr_intervals - 1 | |
| 940 | end if | ||
| 941 | end if | ||
| 942 | 670124591 | end function get_interval_index | |
| 943 | |||
| 944 | !> @brief Evaluate the 1D B-spline function at a point | ||
| 945 | !> | ||
| 946 | !> @param[in] bfun The B-spline function | ||
| 947 | !> @param[in] x The point | ||
| 948 | !> @param[in] derivative _(optional)_ Whether to evaluate the derivative of the B-spline function | ||
| 949 | !> | ||
| 950 | !> @return The value of the B-spline function at the point | ||
| 951 | 670124591 | pure real(wp) function evaluate_1d(bfun, x, derivative) result(val) | |
| 952 | !$acc routine seq | ||
| 953 | use m_bspline_functions, only: bspline_eval | ||
| 954 | |||
| 955 | implicit none | ||
| 956 | |||
| 957 | type(BSplineFun), intent(in) :: bfun | ||
| 958 | real(wp), intent(in) :: x | ||
| 959 | logical, intent(in), optional :: derivative | ||
| 960 | |||
| 961 | integer :: i_x, j_minus_i | ||
| 962 | real(wp) :: bspline_vals(0:MAX_BSPLINE_DEGREE) | ||
| 963 | logical :: derivative_ | ||
| 964 | |||
| 965 |
2/2✓ Branch 2 → 3 taken 265327753 times.
✓ Branch 2 → 4 taken 404796838 times.
|
670124591 | if (present(derivative)) then |
| 966 | 265327753 | derivative_ = derivative | |
| 967 | else | ||
| 968 | 404796838 | derivative_ = .false. | |
| 969 | end if | ||
| 970 | |||
| 971 | val = 0._wp | ||
| 972 |
1/2✓ Branch 5 → 6 taken 670124591 times.
✗ Branch 5 → 11 not taken.
|
670124591 | if (bfun%bspline%nr_bsplines == 0) return |
| 973 | |||
| 974 | 670124591 | i_x = get_interval_index(bfun%bspline, x) | |
| 975 | 670124591 | call evaluate_support_1d(bfun%bspline, i_x, x, derivative_, bspline_vals) | |
| 976 |
2/2✓ Branch 8 → 9 taken 2543614481 times.
✓ Branch 8 → 10 taken 670124591 times.
|
3213739072 | do j_minus_i = 0, bfun%bspline%degree |
| 977 | 3213739072 | val = val + bfun%data(j_minus_i + i_x) * bspline_vals(j_minus_i) | |
| 978 | end do | ||
| 979 | end function evaluate_1d | ||
| 980 | |||
| 981 | !> @brief Evaluate the B-spline basis functions that have support on the interval i at a point | ||
| 982 | !> | ||
| 983 | !> @param[in] bspline The B-spline space | ||
| 984 | !> @param[in] i The interval index | ||
| 985 | !> @param[in] x The point | ||
| 986 | !> @param[in] derivative Whether to evaluate the derivative of the B-spline basis functions | ||
| 987 | !> @param[out] bspline_vals The values of the B-spline basis functions at the point | ||
| 988 | 1313628731 | pure subroutine evaluate_support_1d(bspline, i, x, derivative, bspline_vals) | |
| 989 | !$acc routine seq | ||
| 990 | use m_bspline_functions, only: bspline_eval, bspline_eval_derivative | ||
| 991 | |||
| 992 | implicit none | ||
| 993 | |||
| 994 | type(BSplineSpace), intent(in) :: bspline | ||
| 995 | integer, intent(in) :: i | ||
| 996 | real(wp), intent(in) :: x | ||
| 997 | logical, intent(in) :: derivative | ||
| 998 | real(wp), intent(out) :: bspline_vals(0:bspline%degree) | ||
| 999 | |||
| 1000 | integer :: j_minus_i, j_internal, j | ||
| 1001 | logical :: bspline_mirrored | ||
| 1002 | real(wp) :: x_eval, bspline_val | ||
| 1003 | |||
| 1004 | ✗ | do j_minus_i = 0, bspline%degree | |
| 1005 | ✗ | call get_internal_and_actual_inds(bspline, i, j_minus_i, j, j_internal, bspline_mirrored) | |
| 1006 | |||
| 1007 | ✗ | if (bspline_mirrored) then | |
| 1008 | 244255941 | x_eval = i + 1 - x * bspline%nr_intervals | |
| 1009 |
2/2✓ Branch 6 → 7 taken 8709 times.
✓ Branch 6 → 8 taken 244247232 times.
|
244255941 | if (derivative) then |
| 1010 | bspline_val = -bspline%nr_intervals * bspline_eval_derivative(x_eval, j_internal, bspline%degree - j_minus_i, & | ||
| 1011 | 8709 | bspline%degree, .false.) | |
| 1012 | else | ||
| 1013 | 244247232 | bspline_val = bspline_eval(x_eval, j_internal, bspline%degree - j_minus_i, bspline%degree, bspline%is_scaled) | |
| 1014 | end if | ||
| 1015 | else | ||
| 1016 | ✗ | x_eval = x * bspline%nr_intervals - i | |
| 1017 | ✗ | if (derivative) then | |
| 1018 | 981581394 | bspline_val = bspline%nr_intervals * bspline_eval_derivative(x_eval, j_internal, j_minus_i, bspline%degree, .false.) | |
| 1019 | else | ||
| 1020 | 3848145730 | bspline_val = bspline_eval(x_eval, j_internal, j_minus_i, bspline%degree, bspline%is_scaled) | |
| 1021 | end if | ||
| 1022 | end if | ||
| 1023 | ✗ | if (bspline%is_scaled) then | |
| 1024 | 1105868480 | bspline_val = bspline_val * bspline%nr_intervals | |
| 1025 | end if | ||
| 1026 | ✗ | bspline_vals(j_minus_i) = bspline_val | |
| 1027 | end do | ||
| 1028 | |||
| 1029 | 1313628731 | end subroutine | |
| 1030 | |||
| 1031 | !> @brief Get a PETSc vector containing the B-spline function data | ||
| 1032 | !> | ||
| 1033 | !> @param[out] vec The PETSc vector | ||
| 1034 | !> @param[in] bfun The B-spline function | ||
| 1035 | 7824 | subroutine get_petsc_1dvec(vec, bfun) | |
| 1036 | implicit none | ||
| 1037 | |||
| 1038 | Vec, intent(out) :: vec | ||
| 1039 | type(BSplineFun), intent(in) :: bfun | ||
| 1040 | |||
| 1041 | integer :: ierr, i | ||
| 1042 | integer, allocatable :: inds(:) | ||
| 1043 | |||
| 1044 |
6/12✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 7824 times.
✓ Branch 4 → 3 taken 7824 times.
✗ Branch 4 → 5 not taken.
✓ Branch 5 → 6 taken 7824 times.
✗ Branch 5 → 7 not taken.
✓ Branch 7 → 8 taken 7824 times.
✗ Branch 7 → 9 not taken.
✗ Branch 9 → 10 not taken.
✓ Branch 9 → 11 taken 7824 times.
✗ Branch 11 → 12 not taken.
✓ Branch 11 → 13 taken 7824 times.
|
23472 | allocate (inds(size(bfun%bspline))) |
| 1045 |
2/2✓ Branch 13 → 14 taken 193143 times.
✓ Branch 13 → 15 taken 7824 times.
|
200967 | do i = 1, size(bfun%bspline) |
| 1046 | 200967 | inds(i) = i - 1 | |
| 1047 | end do | ||
| 1048 | |||
| 1049 |
1/2✗ Branch 17 → 18 not taken.
✓ Branch 17 → 19 taken 7824 times.
|
7824 | PetscCall(VecCreate(PETSC_COMM_SELF, vec, ierr)) |
| 1050 |
1/2✗ Branch 20 → 21 not taken.
✓ Branch 20 → 22 taken 7824 times.
|
7824 | PetscCall(VecSetSizes(vec, size(bfun%bspline), size(bfun%bspline), ierr)) |
| 1051 |
1/2✗ Branch 23 → 24 not taken.
✓ Branch 23 → 25 taken 7824 times.
|
7824 | PetscCall(VecSetFromOptions(vec, ierr)) |
| 1052 |
1/2✗ Branch 26 → 27 not taken.
✓ Branch 26 → 28 taken 7824 times.
|
7824 | PetscCall(VecSetValues(vec, size(bfun%bspline), inds, bfun%data, INSERT_VALUES, ierr)) |
| 1053 |
1/2✗ Branch 29 → 30 not taken.
✓ Branch 29 → 31 taken 7824 times.
|
7824 | PetscCall(VecAssemblyBegin(vec, ierr)) |
| 1054 |
1/2✗ Branch 32 → 33 not taken.
✓ Branch 32 → 34 taken 7824 times.
|
7824 | PetscCall(VecAssemblyEnd(vec, ierr)) |
| 1055 | |||
| 1056 |
1/2✗ Branch 34 → 35 not taken.
✓ Branch 34 → 36 taken 7824 times.
|
7824 | deallocate (inds) |
| 1057 | ✗ | end subroutine get_petsc_1dvec | |
| 1058 |
7/16__m_bspline_basis_MOD___copy_m_bspline_basis_Bsplinefun:
✗ Branch 2 → 3 not taken.
✗ Branch 2 → 6 not taken.
✗ Branch 3 → 4 not taken.
✗ Branch 3 → 5 not taken.
__m_bspline_basis_MOD___final_m_bspline_basis_Bsplinefun:
✗ Branch 3 → 4 not taken.
✓ Branch 3 → 7 taken 18551 times.
✗ Branch 4 → 5 not taken.
✗ Branch 4 → 6 not taken.
✓ Branch 8 → 9 taken 18551 times.
✓ Branch 8 → 17 taken 18551 times.
✗ Branch 10 → 11 not taken.
✓ Branch 10 → 12 taken 18551 times.
✓ Branch 12 → 13 taken 18551 times.
✗ Branch 12 → 16 not taken.
✓ Branch 13 → 14 taken 328 times.
✓ Branch 13 → 15 taken 18223 times.
|
37102 | end module m_bspline_basis |
| 1059 |