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, BSplineEvalData, init | ||
| 22 | public :: MAX_BSPLINE_DEGREE | ||
| 23 | public :: derivative, anti_derivative | ||
| 24 | public :: evaluate, 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 Type to store the values of the B-spline basis functions or their derivatives at a given point | ||
| 74 | type BSplineEvalData | ||
| 75 | integer :: interval_index | ||
| 76 | real(wp) :: vals(0:MAX_BSPLINE_DEGREE) | ||
| 77 | ! contains | ||
| 78 | ! procedure :: init => init_bspline_eval_data | ||
| 79 | end type | ||
| 80 | |||
| 81 | interface init | ||
| 82 | module procedure init_bspline_eval_data | ||
| 83 | end interface | ||
| 84 | |||
| 85 | !> @brief Compute quadrature of a B-spline | ||
| 86 | interface quadrature | ||
| 87 | procedure quadrature_eval, quadrature_interval | ||
| 88 | end interface | ||
| 89 | |||
| 90 | !> @brief Compute quadrature of the product of two B-splines | ||
| 91 | interface quadrature_product | ||
| 92 | procedure quadrature_product_eval, quadrature_product_interval | ||
| 93 | end interface | ||
| 94 | |||
| 95 | !> @brief Evaluate a B-spline (function) | ||
| 96 | interface evaluate | ||
| 97 | procedure evaluate_1d, evaluate_single | ||
| 98 | end interface | ||
| 99 | |||
| 100 | !> @brief Get the size of the B-spline space (i.e., the number of linearly independent B-splines) | ||
| 101 | interface size | ||
| 102 | procedure size_1d | ||
| 103 | end interface | ||
| 104 | |||
| 105 | !> @brief Get the index ranges for integrating a B-spline or the product of two B-splines | ||
| 106 | interface get_index_ranges | ||
| 107 | procedure get_index_ranges_single, get_index_ranges_product | ||
| 108 | end interface | ||
| 109 | |||
| 110 | !> @brief Create a B-spline space | ||
| 111 | interface BSplineSpace | ||
| 112 | module procedure init_BSplineSpace | ||
| 113 | end interface | ||
| 114 | |||
| 115 | !> @brief Convert a B-spline function to a PETSc vector | ||
| 116 | interface get_petsc | ||
| 117 | module procedure get_petsc_1dvec | ||
| 118 | end interface | ||
| 119 | contains | ||
| 120 | |||
| 121 | !> @brief Create a B-spline space | ||
| 122 | !> | ||
| 123 | !> @param[in] nr_intervals The number of intervals | ||
| 124 | !> @param[in] degree The degree of the B-splines | ||
| 125 | !> @param[in] is_scaled _(optional)_ If true, the B-splines are scaled (i.e., this space can be used to express the derivative of | ||
| 126 | !> an unscaled B-spline of one degree higher) | ||
| 127 | !> @param[in] is_periodic _(optional)_ If true, the B-splines are periodic, otherwise they are clamped | ||
| 128 | !> | ||
| 129 | !> @return The B-spline space | ||
| 130 | !> | ||
| 131 | !> @note This requires the number of intervals to be larger than the degree | ||
| 132 | 99237 | pure function init_BSplineSpace(nr_intervals, degree, is_scaled, is_periodic) result(ans) | |
| 133 | |||
| 134 | implicit none | ||
| 135 | |||
| 136 | integer, intent(in) :: nr_intervals | ||
| 137 | integer, intent(in) :: degree | ||
| 138 | logical, intent(in), optional :: is_scaled, is_periodic | ||
| 139 | |||
| 140 | type(BSplineSpace) :: ans | ||
| 141 | character(len=256) :: error_msg | ||
| 142 | |||
| 143 |
2/2✓ Branch 2 → 3 taken 97882 times.
✓ Branch 2 → 4 taken 1355 times.
|
99237 | if (present(is_periodic)) then |
| 144 | 97882 | ans%is_periodic = is_periodic | |
| 145 | else | ||
| 146 | ans%is_periodic = .false. | ||
| 147 | end if | ||
| 148 | |||
| 149 |
2/2✓ Branch 4 → 5 taken 76110 times.
✓ Branch 4 → 6 taken 23127 times.
|
99237 | if (present(is_scaled)) then |
| 150 | 76110 | ans%is_scaled = is_scaled | |
| 151 | else | ||
| 152 | ans%is_scaled = .false. | ||
| 153 | end if | ||
| 154 | |||
| 155 |
5/8✓ Branch 6 → 7 taken 4443 times.
✓ Branch 6 → 15 taken 94794 times.
✓ Branch 7 → 8 taken 4443 times.
✗ Branch 7 → 10 not taken.
✓ Branch 8 → 9 taken 4443 times.
✗ Branch 8 → 10 not taken.
✗ Branch 9 → 10 not taken.
✓ Branch 9 → 15 taken 4443 times.
|
99237 | if (degree < 0 .and. .not. (degree == -1 .and. nr_intervals == 1 .and. ans%is_scaled)) then |
| 156 | ✗ | write (error_msg, '(A,I0)') "ERROR in BSplineSpace: degree must be non-negative, got degree = ", degree | |
| 157 | ✗ | error stop error_msg | |
| 158 | end if | ||
| 159 |
1/2✗ Branch 15 → 16 not taken.
✓ Branch 15 → 23 taken 99237 times.
|
99237 | if (degree > MAX_BSPLINE_DEGREE) then |
| 160 | ✗ | write (error_msg, '(A,I0,A,I0)') "ERROR in BSplineSpace: degree ", degree, & | |
| 161 | ✗ | " exceeds MAX_BSPLINE_DEGREE = ", MAX_BSPLINE_DEGREE | |
| 162 | ✗ | error stop error_msg | |
| 163 | end if | ||
| 164 |
1/2✗ Branch 23 → 24 not taken.
✓ Branch 23 → 31 taken 99237 times.
|
99237 | if (nr_intervals <= degree) then |
| 165 | write (error_msg, '(A,I0,A,I0)') "ERROR in BSplineSpace: periodic B-splines require nr_intervals > degree to avoid "// & | ||
| 166 | ✗ | & " multi-valuedness. Got nr_intervals = ", nr_intervals, " and degree = ", degree | |
| 167 | ✗ | error stop error_msg | |
| 168 | end if | ||
| 169 | |||
| 170 |
2/2✓ Branch 31 → 32 taken 38980 times.
✓ Branch 31 → 33 taken 60257 times.
|
99237 | if (ans%is_periodic) then |
| 171 | ans%nr_bsplines = nr_intervals ! Continuity is imposed | ||
| 172 | else | ||
| 173 | 38980 | ans%nr_bsplines = nr_intervals + degree | |
| 174 | end if | ||
| 175 | ans%nr_intervals = nr_intervals | ||
| 176 | ans%degree = degree | ||
| 177 | 99237 | ans%h = 1._wp / nr_intervals | |
| 178 | |||
| 179 | 99237 | end function | |
| 180 | |||
| 181 | !> @brief Returns the B-spline space containing the derivatives of the given B-spline space | ||
| 182 | !> | ||
| 183 | !> @param[in] bspline The B-spline space | ||
| 184 | !> @param[in] allow_periodic_degree0 Whether to allow periodic degree 0 B-splines in the 0-form space (default: false) | ||
| 185 | !> | ||
| 186 | !> @return The B-spline space containing the derivatives of the given B-spline space | ||
| 187 | !> | ||
| 188 | !> @note This requires the B-spline space to be unscaled | ||
| 189 | 75747 | pure type(BSplineSpace) function derivative(bspline, allow_periodic_degree0) result(dbspline) | |
| 190 | implicit none | ||
| 191 | |||
| 192 | type(BSplineSpace), intent(in) :: bspline | ||
| 193 | logical, intent(in), optional :: allow_periodic_degree0 | ||
| 194 | character(len=256) :: error_msg | ||
| 195 | |||
| 196 | logical :: allow_periodic_degree0_ | ||
| 197 | integer :: derivative_degree | ||
| 198 | |||
| 199 | allow_periodic_degree0_ = .false. | ||
| 200 |
2/2✓ Branch 2 → 3 taken 1404 times.
✓ Branch 2 → 4 taken 74343 times.
|
75747 | if (present(allow_periodic_degree0)) allow_periodic_degree0_ = allow_periodic_degree0 |
| 201 | |||
| 202 |
1/2✗ Branch 4 → 5 not taken.
✓ Branch 4 → 9 taken 75747 times.
|
75747 | if (bspline%is_scaled) then |
| 203 | write (error_msg, '(A)') "ERROR in BSplineSpace::derivative: the derivative of scaled B-spline space is not implemented. " & | ||
| 204 | ✗ | //"B-spline space must be unscaled." | |
| 205 | ✗ | error stop error_msg | |
| 206 | end if | ||
| 207 | |||
| 208 | 75747 | derivative_degree = bspline%degree - 1 | |
| 209 |
4/4✓ Branch 9 → 10 taken 1404 times.
✓ Branch 9 → 12 taken 74343 times.
✓ Branch 10 → 11 taken 887 times.
✓ Branch 10 → 12 taken 517 times.
|
75747 | if (allow_periodic_degree0_ .and. bspline%is_periodic) derivative_degree = max(derivative_degree, 0) |
| 210 | 75747 | dbspline = BSplineSpace(bspline%nr_intervals, derivative_degree, is_scaled=.true., is_periodic=bspline%is_periodic) | |
| 211 | 75747 | end function derivative | |
| 212 | |||
| 213 | !> @brief Returns the B-spline space containing the anti-derivatives of the given B-spline space | ||
| 214 | !> | ||
| 215 | !> @param[in] bspline The B-spline space | ||
| 216 | !> | ||
| 217 | !> @return The B-spline space containing the anti-derivatives of the given B-spline space | ||
| 218 | !> | ||
| 219 | !> @note This requires the B-spline space to be scaled | ||
| 220 | 333 | pure type(BSplineSpace) function anti_derivative(bspline) result(dbspline) | |
| 221 | implicit none | ||
| 222 | |||
| 223 | type(BSplineSpace), intent(in) :: bspline | ||
| 224 | character(len=256) :: error_msg | ||
| 225 | |||
| 226 |
1/2✗ Branch 2 → 3 not taken.
✓ Branch 2 → 7 taken 333 times.
|
333 | if (.not. bspline%is_scaled) then |
| 227 | write (error_msg, '(A)') "ERROR in BSplineSpace::anti_derivative: the anti-derivative of unscaled B-spline space is not"& | ||
| 228 | ✗ | & //" implemented. B-spline space must be scaled." | |
| 229 | ✗ | error stop error_msg | |
| 230 | end if | ||
| 231 | |||
| 232 | 333 | dbspline = BSplineSpace(bspline%nr_intervals, bspline%degree + 1, is_scaled=.false., is_periodic=bspline%is_periodic) | |
| 233 | 333 | end function anti_derivative | |
| 234 | |||
| 235 | !> @brief Determine the number of B-splines in the space | ||
| 236 | !> | ||
| 237 | !> @param[in] bspline The B-spline space | ||
| 238 | !> | ||
| 239 | !> @return The number of B-splines | ||
| 240 | 1391911 | pure integer function size_1d(bspline) | |
| 241 | implicit none | ||
| 242 | type(BSplineSpace), intent(in) :: bspline | ||
| 243 | 1402383 | size_1d = bspline%nr_bsplines | |
| 244 | 1391911 | end function | |
| 245 | |||
| 246 | !> @brief Evaluate the j-th B-spline of the B-spline space 'this' at x | ||
| 247 | !> | ||
| 248 | !> @param[in] this The B-spline space | ||
| 249 | !> @param[in] x The evaluation point | ||
| 250 | !> @param[in] j The B-spline index | ||
| 251 | !> @param[in] derivative _(optional)_ If true, the derivative of the B-spline is evaluated | ||
| 252 | !> | ||
| 253 | !> @return The value of the B-spline at x | ||
| 254 | 92760468 | pure real(wp) function evaluate_single(this, x, j, derivative) result(ans) | |
| 255 | use m_bspline_functions | ||
| 256 | |||
| 257 | implicit none | ||
| 258 | |||
| 259 | type(BSplineSpace), intent(in) :: this | ||
| 260 | real(wp), intent(in) :: x | ||
| 261 | integer, intent(in) :: j | ||
| 262 | logical, intent(in), optional :: derivative | ||
| 263 | |||
| 264 | integer :: j_minus_i, j_actual, minus_x | ||
| 265 | real(wp) :: x_eval | ||
| 266 | logical :: derivative_ | ||
| 267 | |||
| 268 |
2/2✓ Branch 2 → 3 taken 785200 times.
✓ Branch 2 → 4 taken 45595034 times.
|
46380234 | if (present(derivative)) then |
| 269 | 785200 | derivative_ = derivative | |
| 270 | else | ||
| 271 | derivative_ = .false. | ||
| 272 | end if | ||
| 273 | |||
| 274 | 46380234 | call determine_evaluation_index(this, x, j, x_eval, j_minus_i, j_actual, minus_x=minus_x) | |
| 275 | |||
| 276 |
4/4✓ Branch 5 → 6 taken 42848977 times.
✓ Branch 5 → 12 taken 3531257 times.
✓ Branch 6 → 7 taken 34298632 times.
✓ Branch 6 → 12 taken 8550345 times.
|
46380234 | if (this%degree < j_minus_i .or. j_minus_i < 0) then |
| 277 | ans = 0._wp | ||
| 278 | else | ||
| 279 |
2/2✓ Branch 7 → 8 taken 180700 times.
✓ Branch 7 → 9 taken 34117932 times.
|
34298632 | if (derivative_) then |
| 280 | 180700 | ans = bspline_eval_derivative(x_eval, j_actual, j_minus_i, this%degree, this%is_scaled) | |
| 281 | 180700 | ans = ans * minus_x * this%nr_intervals | |
| 282 | else | ||
| 283 | 34117932 | ans = bspline_eval(x_eval, j_actual, j_minus_i, this%degree, this%is_scaled) | |
| 284 | end if | ||
| 285 |
2/2✓ Branch 10 → 11 taken 158600 times.
✓ Branch 10 → 12 taken 34140032 times.
|
34298632 | if (this%is_scaled) then |
| 286 | 158600 | ans = ans * this%nr_intervals | |
| 287 | end if | ||
| 288 | end if | ||
| 289 | 46380234 | end function | |
| 290 | |||
| 291 | !> @brief Find a root of the B-spline function 'this' minus rhs using Newton's method | ||
| 292 | !> | ||
| 293 | !> @param[in] this The B-spline function | ||
| 294 | !> @param[in] x0 _(optional)_ The initial guess for the root (default: 0.5) | ||
| 295 | !> @param[in] rhs _(optional)_ The right-hand side value to find the root of (default: 0) | ||
| 296 | !> @param[in] max_iter _(optional)_ The maximum number of iterations for the root-finding method (default: 20) | ||
| 297 | !> @param[in] tol _(optional)_ The tolerance for convergence of the root-finding method (default: 1e-8) | ||
| 298 | !> @param[out] success _(optional)_ A logical flag indicating whether the root-finding method converged successfully | ||
| 299 | !> | ||
| 300 | !> @return The root of the B-spline function 'this' minus rhs | ||
| 301 | 91 | real(wp) function find_root_newton(this, x0, rhs, max_iter, tol, success) result(ans) | |
| 302 | implicit none | ||
| 303 | |||
| 304 | type(BSplineFun), intent(in) :: this | ||
| 305 | real(wp), optional, intent(in) :: x0 | ||
| 306 | real(wp), optional, intent(in) :: rhs | ||
| 307 | integer, optional, intent(in) :: max_iter | ||
| 308 | real(wp), optional, intent(in) :: tol | ||
| 309 | logical, optional, intent(out) :: success | ||
| 310 | |||
| 311 | real(wp) :: rhs_, x0_, tol_, error_est, f, df, update | ||
| 312 | integer :: max_iter_, iter | ||
| 313 | logical :: success_ | ||
| 314 | |||
| 315 | rhs_ = 0.0_wp | ||
| 316 |
1/2✓ Branch 2 → 3 taken 91 times.
✗ Branch 2 → 4 not taken.
|
91 | if (present(rhs)) rhs_ = rhs |
| 317 | |||
| 318 | x0_ = 0.5_wp | ||
| 319 |
1/2✗ Branch 4 → 5 not taken.
✓ Branch 4 → 6 taken 91 times.
|
91 | if (present(x0)) x0_ = x0 |
| 320 | |||
| 321 | max_iter_ = 20 | ||
| 322 |
1/2✗ Branch 6 → 7 not taken.
✓ Branch 6 → 8 taken 91 times.
|
91 | if (present(max_iter)) max_iter_ = max_iter |
| 323 | |||
| 324 | tol_ = 1e-6_wp | ||
| 325 |
1/2✓ Branch 8 → 9 taken 91 times.
✗ Branch 8 → 10 not taken.
|
91 | if (present(tol)) tol_ = tol |
| 326 | |||
| 327 |
2/2✓ Branch 10 → 11 taken 1 time.
✓ Branch 10 → 12 taken 90 times.
|
91 | if (this%bspline%degree == 0) max_iter_ = 0 |
| 328 | |||
| 329 | 91 | ans = x0_ | |
| 330 | |||
| 331 | success_ = .false. | ||
| 332 |
2/2✓ Branch 13 → 14 taken 286 times.
✓ Branch 13 → 30 taken 1 time.
|
287 | do iter = 1, max_iter_ |
| 333 | 286 | f = evaluate(this, ans) - rhs_ | |
| 334 |
2/2✓ Branch 14 → 15 taken 70 times.
✓ Branch 14 → 16 taken 216 times.
|
286 | if (abs(f) < tol_) then |
| 335 | success_ = .true. | ||
| 336 | exit | ||
| 337 | end if | ||
| 338 | |||
| 339 | 216 | df = evaluate(this, ans, derivative=.true.) | |
| 340 |
2/2✓ Branch 16 → 17 taken 20 times.
✓ Branch 16 → 18 taken 196 times.
|
216 | if (abs(df) < 1000 * epsilon(1._wp)) exit |
| 341 | 196 | update = f / df | |
| 342 | 196 | ans = ans - update | |
| 343 | |||
| 344 |
2/4✓ Branch 18 → 19 taken 196 times.
✗ Branch 18 → 20 not taken.
✗ Branch 19 → 20 not taken.
✓ Branch 19 → 28 taken 196 times.
|
196 | if (ans < 0._wp .or. ans > 1._wp) then |
| 345 | ✗ | if (this%bspline%is_periodic) then | |
| 346 | ✗ | ans = modulo(ans, 1._wp) | |
| 347 | ✗ | elseif (ans < 0._wp) then | |
| 348 | ✗ | ans = 0._wp | |
| 349 | else | ||
| 350 | ✗ | ans = 1._wp | |
| 351 | end if | ||
| 352 | end if | ||
| 353 | |||
| 354 |
1/2✗ Branch 28 → 15 not taken.
✓ Branch 28 → 29 taken 196 times.
|
217 | if (abs(update) < tol_) then |
| 355 | success_ = .true. | ||
| 356 | exit | ||
| 357 | end if | ||
| 358 | end do | ||
| 359 | |||
| 360 |
1/2✓ Branch 31 → 32 taken 91 times.
✗ Branch 31 → 33 not taken.
|
91 | if (present(success)) success = success_ |
| 361 | |||
| 362 | 91 | end function | |
| 363 | |||
| 364 | !> @brief Find a root of the B-spline function 'this' minus rhs using Brent's method | ||
| 365 | !> | ||
| 366 | !> @param[in] this The B-spline function | ||
| 367 | !> @param[in] x0 The left-hand side of the bracketed interval for the root | ||
| 368 | !> @param[in] x1 The right-hand side of the bracketed interval for the root | ||
| 369 | !> @param[in] rhs _(optional)_ The right-hand side value to find the root of (default: 0) | ||
| 370 | !> @param[in] max_iter _(optional)_ The maximum number of iterations for the root-finding method (default: 20) | ||
| 371 | !> @param[in] tol _(optional)_ The tolerance for convergence of the root-finding method (default: 1e-6) | ||
| 372 | !> @param[out] success _(optional)_ A logical flag indicating whether the root-finding method converged successfully | ||
| 373 | !> @param[in] res0 _(optional)_ The value of the function minus rhs at x0 (if cheaply available) | ||
| 374 | !> @param[in] res1 _(optional)_ The value of the function minus rhs at x1 (if cheaply available) | ||
| 375 | 91 | real(wp) function find_root_brent(this, x0, x1, rhs, max_iter, tol, success, res0, res1) result(ans) | |
| 376 | use m_common, only: brent | ||
| 377 | implicit none | ||
| 378 | |||
| 379 | type(BSplineFun), intent(in) :: this | ||
| 380 | real(wp), intent(in) :: x0, x1 | ||
| 381 | real(wp), optional, intent(in) :: rhs | ||
| 382 | integer, optional, intent(in) :: max_iter | ||
| 383 | real(wp), optional, intent(in) :: tol | ||
| 384 | logical, optional, intent(out) :: success | ||
| 385 | real(wp), optional, intent(in) :: res0, res1 | ||
| 386 | |||
| 387 | real(wp) :: rhs_, x0_, tol_, error_est, f, df, update | ||
| 388 | integer :: max_iter_, iter | ||
| 389 | |||
| 390 | 91 | rhs_ = 0.0_wp | |
| 391 |
1/2✓ Branch 2 → 3 taken 91 times.
✗ Branch 2 → 4 not taken.
|
91 | if (present(rhs)) rhs_ = rhs |
| 392 | |||
| 393 | 91 | max_iter_ = 20 | |
| 394 |
1/2✗ Branch 4 → 5 not taken.
✓ Branch 4 → 6 taken 91 times.
|
91 | if (present(max_iter)) max_iter_ = max_iter |
| 395 | |||
| 396 | 91 | tol_ = 1e-6_wp | |
| 397 |
1/2✓ Branch 6 → 7 taken 91 times.
✗ Branch 6 → 8 not taken.
|
91 | if (present(tol)) tol_ = tol |
| 398 | |||
| 399 |
2/2✓ Branch 8 → 9 taken 1 time.
✓ Branch 8 → 10 taken 90 times.
|
91 | if (this%bspline%degree == 0) max_iter_ = 0 |
| 400 | |||
| 401 | 91 | ans = brent(f_brent, x0, x1, tol_, max_iter_, fa_in=res0, fb_in=res1, success=success) | |
| 402 | contains | ||
| 403 | 628 | pure real(wp) function f_brent(x) result(fx) | |
| 404 | implicit none | ||
| 405 | real(wp), intent(in) :: x | ||
| 406 | 628 | fx = evaluate(this, x) - rhs_ | |
| 407 | 628 | end function | |
| 408 | end function | ||
| 409 | |||
| 410 | !> @brief Determine the B-spline index and the interval index | ||
| 411 | !> | ||
| 412 | !> @param[in] this The B-spline space | ||
| 413 | !> @param[in] x The evaluation point | ||
| 414 | !> @param[in] j The B-spline index (0:nr_bsplines-1) | ||
| 415 | !> @param[out] x_eval The evaluation point in the interval [0, 1) | ||
| 416 | !> @param[out] j_minus_i The difference between the B-spline index and the interval index | ||
| 417 | !> @param[out] j_actual The internal B-spline index (0:degree) | ||
| 418 | !> @param[out] minus_x _(optional)_ Equals -1 if clamped and the j-th B-spline is on the right boundary, equals 1 otherwise | ||
| 419 | 46380234 | pure subroutine determine_evaluation_index(this, x, j, x_eval, j_minus_i, j_actual, minus_x) | |
| 420 | implicit none | ||
| 421 | |||
| 422 | type(BSplineSpace), intent(in) :: this | ||
| 423 | real(wp), intent(in) :: x | ||
| 424 | integer, intent(in) :: j | ||
| 425 | |||
| 426 | real(wp), intent(out) :: x_eval | ||
| 427 | integer, intent(out) :: j_minus_i, j_actual | ||
| 428 | integer, intent(out), optional :: minus_x | ||
| 429 | |||
| 430 | integer :: i_x | ||
| 431 | |||
| 432 | 46380234 | i_x = int(x * this%nr_intervals) | |
| 433 | |||
| 434 |
1/2✓ Branch 2 → 3 taken 46380234 times.
✗ Branch 2 → 4 not taken.
|
46380234 | if (present(minus_x)) then |
| 435 | 46380234 | minus_x = 1 | |
| 436 | end if | ||
| 437 | |||
| 438 |
2/2✓ Branch 4 → 5 taken 3917192 times.
✓ Branch 4 → 9 taken 42463042 times.
|
46380234 | if (this%is_periodic) then |
| 439 | 3917192 | x_eval = x * this%nr_intervals - i_x | |
| 440 | 3917192 | i_x = mod(i_x, this%nr_intervals) | |
| 441 | |||
| 442 | ! All of the B-splines are the same | ||
| 443 | 3917192 | j_actual = this%degree | |
| 444 |
4/4✓ Branch 5 → 6 taken 1563724 times.
✓ Branch 5 → 8 taken 2353468 times.
✓ Branch 6 → 7 taken 452653 times.
✓ Branch 6 → 8 taken 1111071 times.
|
3917192 | if (i_x > j .and. j < this%degree) then |
| 445 | ! Impose continuity/periodicity if x is on the right side of the domain, | ||
| 446 | ! whereas j refers to a left boundary B-spline, which due to periodicity overlaps | ||
| 447 | ! with a right boundary B-spline | ||
| 448 | 452653 | j_minus_i = this%nr_bsplines + j - i_x | |
| 449 | else | ||
| 450 | 3464539 | j_minus_i = j - i_x | |
| 451 | end if | ||
| 452 | |||
| 453 | else | ||
| 454 | 42463042 | i_x = max(0, min(i_x, this%nr_intervals - 1)) | |
| 455 | 42463042 | j_minus_i = j - i_x | |
| 456 | |||
| 457 |
2/2✓ Branch 9 → 10 taken 13074089 times.
✓ Branch 9 → 11 taken 29388953 times.
|
42463042 | if (j < this%degree) then |
| 458 | 13074089 | j_actual = j | |
| 459 | 13074089 | x_eval = x * this%nr_intervals - i_x | |
| 460 |
2/2✓ Branch 11 → 12 taken 12590658 times.
✓ Branch 11 → 13 taken 16798295 times.
|
29388953 | else if (j < this%nr_intervals) then |
| 461 | ! The {1+degree, ..., nr_intervals} B-splines are all the same | ||
| 462 | 12590658 | j_actual = this%degree | |
| 463 | 12590658 | x_eval = x * this%nr_intervals - i_x | |
| 464 | else | ||
| 465 | 16798295 | j_actual = this%nr_intervals - 1 + this%degree - j | |
| 466 | 16798295 | x_eval = i_x + 1 - x * this%nr_intervals | |
| 467 | 16798295 | j_minus_i = this%degree - j_minus_i | |
| 468 | |||
| 469 |
1/2✓ Branch 13 → 14 taken 16798295 times.
✗ Branch 13 → 15 not taken.
|
16798295 | if (present(minus_x)) then |
| 470 | 16798295 | minus_x = -1 | |
| 471 | end if | ||
| 472 | end if | ||
| 473 | end if | ||
| 474 | |||
| 475 | 46380234 | end subroutine | |
| 476 | |||
| 477 | !> @brief Integrate the j-th B-spline of the B-spline space 'this' using a quadrature rule | ||
| 478 | !> | ||
| 479 | !> @param[in] this The B-spline space | ||
| 480 | !> @param[in] j The B-spline index | ||
| 481 | !> @param[in] userfun _(optional)_ The user-defined function to multiply the B-spline with | ||
| 482 | !> @param[in] n_quad _(optional)_ The number of quadrature points | ||
| 483 | !> | ||
| 484 | !> @return The integral of the B-spline multiplied by the optional user function | ||
| 485 | 82326 | pure real(wp) function quadrature_eval(this, j, userfun, n_quad) result(ans) | |
| 486 | use m_common, only: user_function_1d_interface | ||
| 487 | |||
| 488 | implicit none | ||
| 489 | |||
| 490 | type(BSplineSpace), intent(in) :: this | ||
| 491 | integer, intent(in) :: j | ||
| 492 | procedure(user_function_1d_interface), optional :: userfun | ||
| 493 | integer, intent(in), optional :: n_quad | ||
| 494 | |||
| 495 | integer :: i, n_quad_, i0_min, i0_max, i1_min, i1_max | ||
| 496 |
2/2✓ Branch 2 → 3 taken 17682 times.
✓ Branch 2 → 4 taken 23481 times.
|
41163 | if (present(n_quad)) then |
| 497 | 17682 | n_quad_ = n_quad | |
| 498 | else | ||
| 499 | 23481 | n_quad_ = 1 + int(this%degree / 2) | |
| 500 | end if | ||
| 501 | |||
| 502 | 41163 | call get_index_ranges_single(this, j, i0_min, i0_max, i1_min, i1_max) | |
| 503 | |||
| 504 | ans = 0._wp | ||
| 505 |
2/2✓ Branch 7 → 8 taken 114465 times.
✓ Branch 7 → 12 taken 41163 times.
|
155628 | do i = i0_min, i0_max |
| 506 |
2/2✓ Branch 8 → 9 taken 20846 times.
✓ Branch 8 → 10 taken 93619 times.
|
155628 | if (present(userfun)) then |
| 507 | 20846 | ans = ans + quadrature_interval(this, j, i, n_quad=n_quad_, userfun=userfun) | |
| 508 | else | ||
| 509 | 93619 | ans = ans + quadrature_interval(this, j, i, n_quad=n_quad_) | |
| 510 | end if | ||
| 511 | end do | ||
| 512 |
2/2✓ Branch 14 → 15 taken 8658 times.
✓ Branch 14 → 19 taken 41163 times.
|
49821 | do i = i1_min, i1_max |
| 513 |
2/2✓ Branch 15 → 16 taken 980 times.
✓ Branch 15 → 17 taken 7678 times.
|
49821 | if (present(userfun)) then |
| 514 | 980 | ans = ans + quadrature_interval(this, j, i, n_quad=n_quad_, userfun=userfun) | |
| 515 | else | ||
| 516 | 7678 | ans = ans + quadrature_interval(this, j, i, n_quad=n_quad_) | |
| 517 | end if | ||
| 518 | end do | ||
| 519 | |||
| 520 | 41163 | end function | |
| 521 | |||
| 522 | !> @brief Determine j_actual - i, where the input j may exceed the number of B-splines | ||
| 523 | !> | ||
| 524 | !> @param[in] this The B-spline basis | ||
| 525 | !> @param[in] j The B-spline index | ||
| 526 | !> @param[in] i The interval index | ||
| 527 | !> | ||
| 528 | !> @return The j_actual - i | ||
| 529 | 29772937 | pure integer function get_j_minus_i_quad(this, j, i) result(j_minus_i) | |
| 530 | implicit none | ||
| 531 | |||
| 532 | type(BSplineSpace), intent(in) :: this | ||
| 533 | integer, intent(in) :: j, i | ||
| 534 | |||
| 535 | 29772937 | j_minus_i = j - i | |
| 536 |
2/2✓ Branch 2 → 3 taken 1922826 times.
✓ Branch 2 → 8 taken 27850111 times.
|
29772937 | if (this%is_periodic) then |
| 537 |
2/2✓ Branch 3 → 4 taken 101175 times.
✓ Branch 3 → 5 taken 1821651 times.
|
1922826 | if (j >= this%nr_bsplines) then |
| 538 | 101175 | j_minus_i = j_minus_i - this%nr_bsplines | |
| 539 |
4/4✓ Branch 5 → 6 taken 587609 times.
✓ Branch 5 → 8 taken 1234042 times.
✓ Branch 6 → 7 taken 207457 times.
✓ Branch 6 → 8 taken 380152 times.
|
1821651 | else if (i > j .and. j < this%degree) then |
| 540 | ! Impose continuity/periodicity if x is on the right side of the domain, | ||
| 541 | ! whereas j refers to a left boundary B-spline, which due to periodicity overlaps | ||
| 542 | ! with a right boundary B-spline | ||
| 543 | 207457 | j_minus_i = j_minus_i + this%nr_bsplines | |
| 544 | end if | ||
| 545 | end if | ||
| 546 | |||
| 547 | 29772937 | end function | |
| 548 | |||
| 549 | !> @brief Integrate the j-th B-spline of the B-spline space 'this' over the i-th interval using a quadrature rule | ||
| 550 | !> | ||
| 551 | !> @param[in] this The B-spline space | ||
| 552 | !> @param[in] j The B-spline index | ||
| 553 | !> @param[in] i The interval index | ||
| 554 | !> @param[in] userfun _(optional)_ The user-defined function to multiply the B-spline with | ||
| 555 | !> @param[in] n_quad _(optional)_ The number of quadrature points | ||
| 556 | !> | ||
| 557 | !> @return The integral of the B-spline multiplied by the optional userfun over the interval | ||
| 558 | 944672 | pure real(wp) function quadrature_interval(this, j, i, userfun, n_quad) result(ans) | |
| 559 | use m_common, only: user_function_1d_interface | ||
| 560 | use m_quadrature, only: quadrature | ||
| 561 | |||
| 562 | implicit none | ||
| 563 | |||
| 564 | type(BSplineSpace), intent(in) :: this | ||
| 565 | integer, intent(in) :: j, i | ||
| 566 | procedure(user_function_1d_interface), optional :: userfun | ||
| 567 | integer, intent(in), optional :: n_quad | ||
| 568 | |||
| 569 | integer :: j_minus_i, n_quad_ | ||
| 570 | |||
| 571 |
2/2✓ Branch 2 → 3 taken 123123 times.
✓ Branch 2 → 4 taken 821549 times.
|
944672 | if (present(n_quad)) then |
| 572 | 123123 | n_quad_ = n_quad | |
| 573 | else | ||
| 574 | 821549 | n_quad_ = 1 + int(this%degree / 2) | |
| 575 | end if | ||
| 576 | |||
| 577 | 944672 | j_minus_i = get_j_minus_i_quad(this, j, i) | |
| 578 | |||
| 579 |
4/4✓ Branch 5 → 6 taken 547049 times.
✓ Branch 5 → 8 taken 397623 times.
✓ Branch 6 → 7 taken 155867 times.
✓ Branch 6 → 8 taken 391182 times.
|
944672 | if (j_minus_i > this%degree .or. j_minus_i < 0) then |
| 580 | ans = 0._wp | ||
| 581 | else | ||
| 582 | 155867 | ans = quadrature(n_quad_, evaluate_jth_bspline, i * this%h, (i + 1) * this%h) | |
| 583 | end if | ||
| 584 | contains | ||
| 585 | 353986 | pure real(wp) function evaluate_jth_bspline(x) result(f) | |
| 586 | implicit none | ||
| 587 | |||
| 588 | real(wp), intent(in) :: x | ||
| 589 | |||
| 590 | ! TODO efficiency | ||
| 591 | 353986 | f = evaluate(this, x, j) | |
| 592 |
2/2✓ Branch 2 → 3 taken 103648 times.
✓ Branch 2 → 5 taken 250338 times.
|
353986 | if (present(userfun)) then |
| 593 | 103648 | f = f * userfun(x) | |
| 594 | end if | ||
| 595 | 353986 | end function | |
| 596 | end function | ||
| 597 | |||
| 598 | !> @brief Get the interval index ranges (i.e., the support) for integrating the j-th B-spline | ||
| 599 | !> | ||
| 600 | !> @param[in] this The B-spline space | ||
| 601 | !> @param[in] j The B-spline index | ||
| 602 | !> @param[out] i0_min The minimum interval index for the first range | ||
| 603 | !> @param[out] i0_max The maximum interval index for the first range | ||
| 604 | !> @param[out] i1_min The minimum interval index for the second range | ||
| 605 | !> @param[out] i1_max The maximum interval index for the second range | ||
| 606 | 100900 | pure subroutine get_index_ranges_single(this, j, i0_min, i0_max, i1_min, i1_max) | |
| 607 | implicit none | ||
| 608 | |||
| 609 | type(BSplineSpace), intent(in) :: this | ||
| 610 | |||
| 611 | integer, intent(in) :: j | ||
| 612 | integer, intent(out) :: i0_min, i0_max, i1_min, i1_max | ||
| 613 | |||
| 614 |
4/4✓ Branch 2 → 3 taken 51498 times.
✓ Branch 2 → 4 taken 49402 times.
✓ Branch 3 → 4 taken 32474 times.
✓ Branch 3 → 5 taken 19024 times.
|
100900 | if (.not. this%is_periodic .or. j > this%degree) then |
| 615 | 81876 | i0_min = max(0, j - this%degree) | |
| 616 | 81876 | i0_max = min(this%nr_intervals - 1, j) | |
| 617 | 81876 | i1_min = 1 | |
| 618 | 81876 | i1_max = -1 | |
| 619 | else ! if (j <= this%degree) then | ||
| 620 | 19024 | i0_min = 0 | |
| 621 | 19024 | i0_max = j | |
| 622 | 19024 | i1_min = this%nr_intervals - this%degree + j | |
| 623 | 19024 | i1_max = this%nr_intervals - 1 | |
| 624 | end if | ||
| 625 | 100900 | end subroutine get_index_ranges_single | |
| 626 | |||
| 627 | !> @brief Get the interval index ranges (i.e., the support) for integrating the product of the j1-th and j2-th B-splines | ||
| 628 | !> | ||
| 629 | !> @param[in] this1 The first B-spline space | ||
| 630 | !> @param[in] this2 The second B-spline space | ||
| 631 | !> @param[in] j1 The first B-spline index | ||
| 632 | !> @param[in] j2 The second B-spline index | ||
| 633 | !> @param[out] i0_min The minimum interval index for the first range | ||
| 634 | !> @param[out] i0_max The maximum interval index for the first range | ||
| 635 | !> @param[out] i1_min The minimum interval index for the second range | ||
| 636 | !> @param[out] i1_max The maximum interval index for the second range | ||
| 637 | 8988394 | pure subroutine get_index_ranges_product(this1, this2, j1, j2, i0_min, i0_max, i1_min, i1_max) | |
| 638 | implicit none | ||
| 639 | |||
| 640 | type(BSplineSpace), intent(in) :: this1, this2 | ||
| 641 | |||
| 642 | integer, intent(in) :: j1, j2 | ||
| 643 | integer, intent(out) :: i0_min, i0_max, i1_min, i1_max | ||
| 644 | |||
| 645 |
6/6✓ Branch 2 → 3 taken 745225 times.
✓ Branch 2 → 5 taken 8243169 times.
✓ Branch 3 → 4 taken 727351 times.
✓ Branch 3 → 6 taken 17874 times.
✓ Branch 4 → 5 taken 716157 times.
✓ Branch 4 → 6 taken 11194 times.
|
8988394 | if (.not. this1%is_periodic .or. (j1 >= this1%degree .and. j2 >= this2%degree)) then |
| 646 | 8959326 | i0_min = max(0, j1 - this1%degree, j2 - this2%degree) | |
| 647 | 8959326 | i0_max = min(this1%nr_intervals - 1, j1, j2) | |
| 648 | 8959326 | i1_min = 1 | |
| 649 | 8959326 | i1_max = -1 | |
| 650 |
3/4✓ Branch 6 → 7 taken 11194 times.
✓ Branch 6 → 9 taken 17874 times.
✓ Branch 7 → 8 taken 11194 times.
✗ Branch 7 → 9 not taken.
|
29068 | else if (j1 >= this1%degree .and. j2 < this2%degree) then |
| 651 | ! j1: j1-degree,j1 | ||
| 652 | ! j2: 1,j2 | ||
| 653 | 11194 | i0_min = max(j1 - this1%degree, 0) | |
| 654 | 11194 | i0_max = min(j1, j2) | |
| 655 | ! j1: j1-degree,j1 | ||
| 656 | ! j2: nr_intervals-(degree-j2),nr_intervals | ||
| 657 | 11194 | i1_min = max(j1 - this1%degree, this1%nr_intervals - this2%degree + j2) | |
| 658 | 11194 | i1_max = min(j1, this1%nr_intervals - 1) | |
| 659 |
3/4✓ Branch 9 → 10 taken 17874 times.
✗ Branch 9 → 12 not taken.
✓ Branch 10 → 11 taken 16241 times.
✓ Branch 10 → 12 taken 1633 times.
|
17874 | else if (j1 < this1%degree .and. j2 >= this2%degree) then |
| 660 | ! j1: 1,j1 | ||
| 661 | ! j2: j2-degree,j2 | ||
| 662 | 16241 | i0_min = max(j2 - this2%degree, 0) | |
| 663 | 16241 | i0_max = min(j1, j2) | |
| 664 | ! j1: nr_intervals-(degree-j1),nr_intervals | ||
| 665 | ! j2: j2-degree,j2 | ||
| 666 | 16241 | i1_min = max(j2 - this2%degree, this1%nr_intervals - this1%degree + j1) | |
| 667 | 16241 | i1_max = min(j2, this1%nr_intervals - 1) | |
| 668 | else ! if (j1 < this%degree .and. j2 < this%degree) then | ||
| 669 | ! j1: 1,j1 | ||
| 670 | ! j2: 1,j2 | ||
| 671 | 1633 | i0_min = 0 | |
| 672 | 1633 | i0_max = max(j1, j2) | |
| 673 | ! j1: nr_intervals-(degree-j1),nr_intervals | ||
| 674 | ! j2: nr_intervals-(degree-j2),nr_intervals | ||
| 675 | 1633 | i1_min = max(this1%nr_intervals - this1%degree + j1, this1%nr_intervals - this2%degree + j2) | |
| 676 | 1633 | i1_max = this1%nr_intervals - 1 | |
| 677 | end if | ||
| 678 | 8988394 | end subroutine get_index_ranges_product | |
| 679 | |||
| 680 | !> @brief Integrate the product of the j1-th and j2-th B-splines of the B-spline spaces 'this1' and 'this2' using a quadrature | ||
| 681 | !> rule | ||
| 682 | !> | ||
| 683 | !> @param[in] this1 The first B-spline basis | ||
| 684 | !> @param[in] this2 The second B-spline basis | ||
| 685 | !> @param[in] j1 The first B-spline index | ||
| 686 | !> @param[in] j2 The second B-spline index | ||
| 687 | !> @param[in] userfun _(optional)_ The user-defined function to multiply the B-splines with | ||
| 688 | !> @param[in] n_quad _(optional)_ The number of quadrature points | ||
| 689 | !> | ||
| 690 | !> @return The integral of the product of the B-splines multiplied by the optional userfun | ||
| 691 | !> | ||
| 692 | !> @note It is assumed that the two B-spline spaces have the same number of intervals as well as periodicity. The degree of the | ||
| 693 | !> B-splines may differ. | ||
| 694 | 5827518 | pure real(wp) function quadrature_product_eval(this1, this2, j1, j2, userfun, n_quad) result(ans) | |
| 695 | use m_common, only: user_function_1d_interface | ||
| 696 | |||
| 697 | implicit none | ||
| 698 | |||
| 699 | type(BSplineSpace), intent(in) :: this1, this2 | ||
| 700 | integer, intent(in) :: j1, j2 | ||
| 701 | procedure(user_function_1d_interface), optional :: userfun | ||
| 702 | integer, intent(in), optional :: n_quad | ||
| 703 | |||
| 704 | integer :: i, n_quad_, i0_min, i0_max, i1_min, i1_max | ||
| 705 | |||
| 706 |
2/2✓ Branch 2 → 3 taken 1996335 times.
✓ Branch 2 → 4 taken 917424 times.
|
2913759 | if (present(n_quad)) then |
| 707 | 1996335 | n_quad_ = n_quad | |
| 708 | else | ||
| 709 | 917424 | n_quad_ = 1 + max(this1%degree, this2%degree) | |
| 710 | end if | ||
| 711 | |||
| 712 | 2913759 | call get_index_ranges_product(this1, this2, j1, j2, i0_min, i0_max, i1_min, i1_max) | |
| 713 | |||
| 714 | ans = 0._wp | ||
| 715 |
2/2✓ Branch 7 → 8 taken 2902999 times.
✓ Branch 7 → 12 taken 2913759 times.
|
5816758 | do i = i0_min, i0_max |
| 716 |
1/2✗ Branch 8 → 9 not taken.
✓ Branch 8 → 10 taken 2902999 times.
|
5816758 | if (present(userfun)) then |
| 717 | ✗ | ans = ans + quadrature_product_interval(this1, this2, j1, j2, i, n_quad=n_quad_, userfun=userfun) | |
| 718 | else | ||
| 719 | 2902999 | ans = ans + quadrature_product_interval(this1, this2, j1, j2, i, n_quad=n_quad_) | |
| 720 | end if | ||
| 721 | end do | ||
| 722 |
2/2✓ Branch 14 → 15 taken 7084 times.
✓ Branch 14 → 19 taken 2913759 times.
|
2920843 | do i = i1_min, i1_max |
| 723 |
1/2✗ Branch 15 → 16 not taken.
✓ Branch 15 → 17 taken 7084 times.
|
2920843 | if (present(userfun)) then |
| 724 | ✗ | ans = ans + quadrature_product_interval(this1, this2, j1, j2, i, n_quad=n_quad_, userfun=userfun) | |
| 725 | else | ||
| 726 | 7084 | ans = ans + quadrature_product_interval(this1, this2, j1, j2, i, n_quad=n_quad_) | |
| 727 | end if | ||
| 728 | end do | ||
| 729 | 2913759 | end function | |
| 730 | |||
| 731 | !> @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 | ||
| 732 | !> using a quadrature rule | ||
| 733 | !> | ||
| 734 | !> @param[in] this1 The first B-spline basis | ||
| 735 | !> @param[in] this2 The second B-spline basis | ||
| 736 | !> @param[in] j1 The first B-spline index | ||
| 737 | !> @param[in] j2 The second B-spline index | ||
| 738 | !> @param[in] i The interval index | ||
| 739 | !> @param[in] userfun _(optional)_ The user-defined function to multiply the B-splines with | ||
| 740 | !> @param[in] n_quad _(optional)_ The number of quadrature points | ||
| 741 | !> | ||
| 742 | !> @return The integral of the product of the B-splines multiplied by the optional userfun over the interval | ||
| 743 | !> | ||
| 744 | !> @note It is assumed that the two B-spline spaces have the same number of intervals as well as periodicity. The degree of the | ||
| 745 | !> B-splines may differ. | ||
| 746 | 2910083 | pure real(wp) function quadrature_product_interval(this1, this2, j1, j2, i, userfun, n_quad) result(ans) | |
| 747 | use m_common, only: user_function_1d_interface | ||
| 748 | use m_quadrature, only: quad => quadrature | ||
| 749 | |||
| 750 | implicit none | ||
| 751 | |||
| 752 | type(BSplineSpace), intent(in) :: this1, this2 | ||
| 753 | integer, intent(in) :: j1, j2, i | ||
| 754 | procedure(user_function_1d_interface), optional :: userfun | ||
| 755 | integer, intent(in), optional :: n_quad | ||
| 756 | |||
| 757 | integer :: j1_minus_i, j2_minus_i, n_quad_ | ||
| 758 | |||
| 759 |
1/2✓ Branch 2 → 3 taken 2910083 times.
✗ Branch 2 → 4 not taken.
|
2910083 | if (present(n_quad)) then |
| 760 | 2910083 | n_quad_ = n_quad | |
| 761 | else | ||
| 762 | ✗ | n_quad_ = 1 + max(this1%degree, this2%degree) | |
| 763 | end if | ||
| 764 | |||
| 765 | 2910083 | j1_minus_i = get_j_minus_i_quad(this1, j1, i) | |
| 766 | 2910083 | j2_minus_i = get_j_minus_i_quad(this2, j2, i) | |
| 767 | |||
| 768 |
6/8✓ Branch 5 → 6 taken 2909173 times.
✓ Branch 5 → 10 taken 910 times.
✓ Branch 6 → 7 taken 2909173 times.
✗ Branch 6 → 10 not taken.
✓ Branch 7 → 8 taken 2908263 times.
✓ Branch 7 → 10 taken 910 times.
✓ Branch 8 → 9 taken 2908263 times.
✗ Branch 8 → 10 not taken.
|
2910083 | if (j1_minus_i > this1%degree .or. j1_minus_i < 0 .or. j2_minus_i > this2%degree .or. j2_minus_i < 0) then |
| 769 | ans = 0._wp | ||
| 770 | else | ||
| 771 | 2908263 | ans = quad(n_quad_, evaluate_bspline_product, i * this1%h, (i + 1) * this1%h) | |
| 772 | end if | ||
| 773 | contains | ||
| 774 | 13929471 | pure real(wp) function evaluate_bspline_product(x) result(f) | |
| 775 | implicit none | ||
| 776 | |||
| 777 | real(wp), intent(in) :: x | ||
| 778 | |||
| 779 | ! TODO efficiency: internal evaluation function with internal indices and no bounds checking | ||
| 780 | 13929471 | f = evaluate(this1, x, j1) * evaluate(this2, x, j2) | |
| 781 |
1/2✗ Branch 2 → 3 not taken.
✓ Branch 2 → 5 taken 13929471 times.
|
13929471 | if (present(userfun)) then |
| 782 | ✗ | f = f * userfun(x) | |
| 783 | end if | ||
| 784 | 13929471 | end function | |
| 785 | end function | ||
| 786 | |||
| 787 | !> @brief Get the internal B-spline index and whether the B-spline is mirrored | ||
| 788 | !> | ||
| 789 | !> @param[in] this The B-spline space | ||
| 790 | !> @param[in] j The B-spline index | ||
| 791 | !> @param[out] j_internal The internal B-spline index (0:degree) | ||
| 792 | !> @param[out] bspline_mirrored _(optional)_ Whether the B-spline is mirrored (i.e., the B-spline is on the right boundary of the | ||
| 793 | !> domain) | ||
| 794 | ✗ | pure subroutine get_j_internal(this, j, j_internal, bspline_mirrored) | |
| 795 | !$acc routine seq | ||
| 796 | implicit none | ||
| 797 | |||
| 798 | type(BSplineSpace), intent(in) :: this | ||
| 799 | integer, intent(in) :: j | ||
| 800 | integer, intent(out) :: j_internal | ||
| 801 | logical, intent(out) :: bspline_mirrored | ||
| 802 | |||
| 803 | ✗ | bspline_mirrored = .false. | |
| 804 | ✗ | if (this%is_periodic) then | |
| 805 | ! All of the B-splines are the same | ||
| 806 | 3853102511 | j_internal = this%degree | |
| 807 | else | ||
| 808 |
2/2✓ Branch 4 → 5 taken 675174795 times.
✓ Branch 4 → 6 taken 1794657114 times.
|
2469831909 | if (j < this%degree) then |
| 809 | 675174795 | j_internal = j | |
| 810 |
2/2✓ Branch 6 → 7 taken 1026807690 times.
✓ Branch 6 → 8 taken 767849424 times.
|
1794657114 | else if (j < this%nr_intervals) then |
| 811 | ! The {1+degree, ..., nr_intervals} B-splines are all the same | ||
| 812 | 1026807690 | j_internal = this%degree | |
| 813 | else | ||
| 814 | 767849424 | j_internal = this%nr_intervals - 1 + this%degree - j | |
| 815 | 767849424 | bspline_mirrored = .true. | |
| 816 | end if | ||
| 817 | end if | ||
| 818 | ✗ | end subroutine get_j_internal | |
| 819 | |||
| 820 | !> @brief Get the internal B-spline index and the actual B-spline index | ||
| 821 | !> | ||
| 822 | !> @param[in] bspline The B-spline space | ||
| 823 | !> @param[in] i The interval index | ||
| 824 | !> @param[in] j_minus_i The difference between the B-spline index and the interval index | ||
| 825 | !> @param[out] j The actual B-spline index (0:nr_bsplines-1) | ||
| 826 | !> @param[out] j_internal _(optional)_ The internal B-spline index (0:degree) | ||
| 827 | !> @param[out] bspline_mirrored _(optional)_ Whether the B-spline is mirrored (i.e., the B-spline is on the right boundary of the | ||
| 828 | !> domain) | ||
| 829 | ✗ | pure subroutine get_internal_and_actual_inds(bspline, i, j_minus_i, j, j_internal, bspline_mirrored) | |
| 830 | !$acc routine seq | ||
| 831 | implicit none | ||
| 832 | |||
| 833 | type(BSplineSpace), intent(in) :: bspline | ||
| 834 | integer, intent(in) :: i, j_minus_i | ||
| 835 | |||
| 836 | integer, intent(out) :: j | ||
| 837 | integer, intent(out), optional :: j_internal | ||
| 838 | logical, intent(out), optional :: bspline_mirrored | ||
| 839 | |||
| 840 | ✗ | j = j_minus_i + i | |
| 841 | ✗ | if (bspline%is_periodic .and. j >= bspline%nr_intervals) then | |
| 842 | 289562887 | j = j - bspline%nr_intervals | |
| 843 | end if | ||
| 844 | |||
| 845 | ✗ | if (present(j_internal) .and. present(bspline_mirrored)) then | |
| 846 | ✗ | call get_j_internal(bspline, j, j_internal, bspline_mirrored) | |
| 847 | end if | ||
| 848 | ✗ | end subroutine get_internal_and_actual_inds | |
| 849 | |||
| 850 | !> @brief Initialize the B-spline function | ||
| 851 | !> | ||
| 852 | !> @param[out] bfun The B-spline function | ||
| 853 | !> @param[in] bspline The B-spline space | ||
| 854 |
1/2✓ Branch 2 → 3 taken 12846 times.
✗ Branch 2 → 5 not taken.
|
12846 | subroutine init_bspline_fun(bfun, bspline, vec) |
| 855 | implicit none | ||
| 856 | |||
| 857 | class(BSplineFun), intent(out) :: bfun | ||
| 858 | type(BSplineSpace), intent(in) :: bspline | ||
| 859 | Vec, intent(in), optional :: vec | ||
| 860 | |||
| 861 | integer :: ierr, i | ||
| 862 | integer, allocatable :: inds(:) | ||
| 863 | |||
| 864 | 12846 | bfun%bspline = bspline | |
| 865 |
1/2✗ Branch 5 → 6 not taken.
✓ Branch 5 → 7 taken 12846 times.
|
12846 | if (allocated(bfun%data)) deallocate (bfun%data) |
| 866 |
7/14✗ Branch 7 → 8 not taken.
✓ Branch 7 → 9 taken 12846 times.
✓ Branch 9 → 8 taken 12846 times.
✗ Branch 9 → 10 not taken.
✓ Branch 10 → 11 taken 12846 times.
✗ Branch 10 → 12 not taken.
✓ Branch 12 → 13 taken 12846 times.
✗ Branch 12 → 14 not taken.
✗ Branch 14 → 15 not taken.
✓ Branch 14 → 16 taken 12846 times.
✗ Branch 16 → 17 not taken.
✓ Branch 16 → 18 taken 12846 times.
✗ Branch 18 → 19 not taken.
✓ Branch 18 → 20 taken 12846 times.
|
38538 | allocate (bfun%data(0:bspline%nr_intervals - 1 + bspline%degree)) |
| 867 | |||
| 868 |
2/2✓ Branch 20 → 21 taken 5236 times.
✓ Branch 20 → 43 taken 7610 times.
|
12846 | if (present(vec)) then |
| 869 |
6/12✗ Branch 21 → 22 not taken.
✓ Branch 21 → 23 taken 5236 times.
✓ Branch 23 → 22 taken 5236 times.
✗ Branch 23 → 24 not taken.
✓ Branch 24 → 25 taken 5236 times.
✗ Branch 24 → 26 not taken.
✓ Branch 26 → 27 taken 5236 times.
✗ Branch 26 → 28 not taken.
✗ Branch 28 → 29 not taken.
✓ Branch 28 → 30 taken 5236 times.
✗ Branch 30 → 31 not taken.
✓ Branch 30 → 32 taken 5236 times.
|
15708 | allocate (inds(size(bspline))) |
| 870 |
2/2✓ Branch 32 → 33 taken 122807 times.
✓ Branch 32 → 34 taken 5236 times.
|
128043 | do i = 1, size(bspline) |
| 871 | 128043 | inds(i) = i - 1 | |
| 872 | end do | ||
| 873 | |||
| 874 |
1/2✗ Branch 36 → 37 not taken.
✓ Branch 36 → 40 taken 5236 times.
|
5236 | PetscCall(VecGetValues(vec, size(bspline), inds, bfun%data, ierr)) |
| 875 |
1/2✗ Branch 40 → 41 not taken.
✓ Branch 40 → 42 taken 5236 times.
|
5236 | deallocate (inds) |
| 876 | |||
| 877 | 5236 | call bfun%periodicity() | |
| 878 | else | ||
| 879 |
2/2✓ Branch 43 → 44 taken 196707 times.
✓ Branch 43 → 45 taken 7610 times.
|
204317 | bfun%data = 0._wp |
| 880 | end if | ||
| 881 | ✗ | end subroutine init_bspline_fun | |
| 882 | |||
| 883 | !> @brief Enforce periodicity on the B-spline function | ||
| 884 | !> | ||
| 885 | !> @param[inout] bfun The B-spline function | ||
| 886 | 5236 | subroutine bspline_fun_periodicity(bfun) | |
| 887 | implicit none | ||
| 888 | |||
| 889 | class(BSplineFun), intent(inout) :: bfun | ||
| 890 | integer :: i | ||
| 891 | |||
| 892 |
2/2✓ Branch 2 → 3 taken 2955 times.
✓ Branch 2 → 7 taken 2281 times.
|
5236 | if (bfun%bspline%is_periodic) then |
| 893 |
2/2✓ Branch 4 → 5 taken 10235 times.
✓ Branch 4 → 6 taken 2955 times.
|
13190 | do i = 0, bfun%bspline%degree - 1 |
| 894 | 13190 | bfun%data(i + bfun%bspline%nr_intervals) = bfun%data(i) | |
| 895 | end do | ||
| 896 | end if | ||
| 897 | 5236 | end subroutine bspline_fun_periodicity | |
| 898 | |||
| 899 | !> @brief Destroy the B-spline function | ||
| 900 | !> | ||
| 901 | !> @param[inout] bfun The B-spline function | ||
| 902 | 2949 | subroutine destroy_bspline_fun(bfun) | |
| 903 | implicit none | ||
| 904 | |||
| 905 | class(BSplineFun), intent(inout) :: bfun | ||
| 906 | |||
| 907 |
1/2✓ Branch 2 → 3 taken 2949 times.
✗ Branch 2 → 4 not taken.
|
2949 | if (allocated(bfun%data)) deallocate (bfun%data) |
| 908 | 2949 | end subroutine destroy_bspline_fun | |
| 909 | |||
| 910 | ! TODO the 'final' functionality is not reliably supported by all (any?) compilers, so we do not use it here | ||
| 911 | |||
| 912 | ! !> @brief Destroy the B-spline function when it goes out of scope (final) | ||
| 913 | ! !> | ||
| 914 | ! !> @param[inout] bfun The B-spline function | ||
| 915 | ! subroutine destroy_bspline_fun_final(bfun) | ||
| 916 | ! implicit none | ||
| 917 | |||
| 918 | ! type(BSplineFun), intent(inout) :: bfun | ||
| 919 | |||
| 920 | ! call bfun%destroy() | ||
| 921 | ! end subroutine destroy_bspline_fun_final | ||
| 922 | |||
| 923 | !> @brief Get the interval index for a point in the B-spline space | ||
| 924 | !> | ||
| 925 | !> @param[in] bspline The B-spline space | ||
| 926 | !> @param[in] x The point | ||
| 927 | !> | ||
| 928 | !> @return The interval index for the point in the B-spline space (0:nr_intervals-1) | ||
| 929 | ✗ | pure integer function get_interval_index(bspline, x) result(i_x) | |
| 930 | !$acc routine seq | ||
| 931 | implicit none | ||
| 932 | |||
| 933 | type(BSplineSpace), intent(in) :: bspline | ||
| 934 | real(wp), intent(in) :: x | ||
| 935 | |||
| 936 | ✗ | i_x = int(x * bspline%nr_intervals) | |
| 937 | ✗ | if (i_x < 0) then | |
| 938 | ✗ | if (bspline%is_periodic) then | |
| 939 | ✗ | i_x = modulo(i_x, bspline%nr_intervals) | |
| 940 | else | ||
| 941 | i_x = 0 | ||
| 942 | end if | ||
| 943 | ✗ | else if (i_x >= bspline%nr_intervals) then | |
| 944 | ✗ | if (bspline%is_periodic) then | |
| 945 | ✗ | i_x = modulo(i_x, bspline%nr_intervals) | |
| 946 | else | ||
| 947 | ✗ | i_x = bspline%nr_intervals - 1 | |
| 948 | end if | ||
| 949 | end if | ||
| 950 | ✗ | end function get_interval_index | |
| 951 | |||
| 952 | !> @brief Initialize the B-spline evaluation data for evaluating the B-spline function at a point | ||
| 953 | !> | ||
| 954 | !> @param[inout] this The B-spline evaluation data | ||
| 955 | !> @param[in] bspace The B-spline space | ||
| 956 | !> @param[in] x The point at which to evaluate the B-spline function | ||
| 957 | !> @param[in] is_derivative Whether to initialize the data for evaluating the derivative of the B-spline function | ||
| 958 | 1280190857 | pure subroutine init_bspline_eval_data(this, bspace, x, is_derivative) | |
| 959 | !$acc routine seq | ||
| 960 | use m_bspline_functions, only: bspline_eval, bspline_eval_derivative | ||
| 961 | implicit none | ||
| 962 | |||
| 963 | type(BSplineEvalData), intent(inout) :: this | ||
| 964 | type(BSplineSpace), intent(in) :: bspace | ||
| 965 | real(wp), intent(in) :: x | ||
| 966 | logical, intent(in) :: is_derivative | ||
| 967 | |||
| 968 | integer :: j_minus_i, j_internal, j, i | ||
| 969 | logical :: bspline_mirrored | ||
| 970 | real(wp) :: x_eval, bspline_val | ||
| 971 | |||
| 972 | 1280190857 | i = int(x * bspace%nr_intervals) | |
| 973 |
2/2✓ Branch 2 → 3 taken 292546103 times.
✓ Branch 2 → 4 taken 987644754 times.
|
1280190857 | if (.not. bspace%is_periodic) i = max(0, min(i, bspace%nr_intervals - 1)) |
| 974 | |||
| 975 | 1280190857 | this%interval_index = i | |
| 976 | |||
| 977 | ✗ | do j_minus_i = 0, bspace%degree | |
| 978 | ✗ | call get_internal_and_actual_inds(bspace, i, j_minus_i, j, j_internal, bspline_mirrored) | |
| 979 | |||
| 980 | ✗ | if (bspline_mirrored) then | |
| 981 | 197401243 | x_eval = i + 1 - x * bspace%nr_intervals | |
| 982 |
2/2✓ Branch 8 → 9 taken 4933 times.
✓ Branch 8 → 10 taken 197396310 times.
|
197401243 | if (is_derivative) then |
| 983 | bspline_val = -bspace%nr_intervals * bspline_eval_derivative(x_eval, j_internal, bspace%degree - j_minus_i, & | ||
| 984 | 4933 | bspace%degree, .false.) | |
| 985 | else | ||
| 986 | 197396310 | bspline_val = bspline_eval(x_eval, j_internal, bspace%degree - j_minus_i, bspace%degree, bspace%is_scaled) | |
| 987 | end if | ||
| 988 | else | ||
| 989 | ✗ | x_eval = x * bspace%nr_intervals - i | |
| 990 | ✗ | if (is_derivative) then | |
| 991 | 981528761 | bspline_val = bspace%nr_intervals * bspline_eval_derivative(x_eval, j_internal, j_minus_i, bspace%degree, .false.) | |
| 992 | else | ||
| 993 | 3713116117 | bspline_val = bspline_eval(x_eval, j_internal, j_minus_i, bspace%degree, bspace%is_scaled) | |
| 994 | end if | ||
| 995 | end if | ||
| 996 | ✗ | if (bspace%is_scaled) then | |
| 997 | 1059376766 | bspline_val = bspline_val * bspace%nr_intervals | |
| 998 | end if | ||
| 999 | ✗ | this%vals(j_minus_i) = bspline_val | |
| 1000 | end do | ||
| 1001 | 1280190857 | end subroutine init_bspline_eval_data | |
| 1002 | |||
| 1003 | !> @brief Evaluate the 1D B-spline function at a point | ||
| 1004 | !> | ||
| 1005 | !> @param[in] bfun The B-spline function | ||
| 1006 | !> @param[in] x The point | ||
| 1007 | !> @param[in] derivative _(optional)_ Whether to evaluate the derivative of the B-spline function | ||
| 1008 | !> | ||
| 1009 | !> @return The value of the B-spline function at the point | ||
| 1010 | 658456499 | pure real(wp) function evaluate_1d(bfun, x, derivative) result(val) | |
| 1011 | !$acc routine seq | ||
| 1012 | |||
| 1013 | implicit none | ||
| 1014 | |||
| 1015 | type(BSplineFun), intent(in) :: bfun | ||
| 1016 | real(wp), intent(in) :: x | ||
| 1017 | logical, intent(in), optional :: derivative | ||
| 1018 | |||
| 1019 | integer :: j_minus_i | ||
| 1020 | logical :: derivative_ | ||
| 1021 | type(BSplineEvalData) :: eval_data | ||
| 1022 | |||
| 1023 |
2/2✓ Branch 2 → 3 taken 265326976 times.
✓ Branch 2 → 4 taken 393129523 times.
|
658456499 | if (present(derivative)) then |
| 1024 | 265326976 | derivative_ = derivative | |
| 1025 | else | ||
| 1026 | 393129523 | derivative_ = .false. | |
| 1027 | end if | ||
| 1028 | |||
| 1029 | val = 0._wp | ||
| 1030 |
1/2✓ Branch 5 → 6 taken 658456499 times.
✗ Branch 5 → 11 not taken.
|
658456499 | if (bfun%bspline%nr_bsplines == 0) return |
| 1031 | |||
| 1032 | 658456499 | call init(eval_data, bfun%bspline, x, derivative_) | |
| 1033 |
2/2✓ Branch 8 → 9 taken 2479291863 times.
✓ Branch 8 → 10 taken 658456499 times.
|
3137748362 | do j_minus_i = 0, bfun%bspline%degree |
| 1034 | 3137748362 | val = val + bfun%data(j_minus_i + eval_data%interval_index) * eval_data%vals(j_minus_i) | |
| 1035 | end do | ||
| 1036 | end function evaluate_1d | ||
| 1037 | |||
| 1038 | !> @brief Get a PETSc vector containing the B-spline function data | ||
| 1039 | !> | ||
| 1040 | !> @param[out] vec The PETSc vector | ||
| 1041 | !> @param[in] bfun The B-spline function | ||
| 1042 | 5236 | subroutine get_petsc_1dvec(vec, bfun) | |
| 1043 | implicit none | ||
| 1044 | |||
| 1045 | Vec, intent(out) :: vec | ||
| 1046 | type(BSplineFun), intent(in) :: bfun | ||
| 1047 | |||
| 1048 | integer :: ierr, i | ||
| 1049 | integer, allocatable :: inds(:) | ||
| 1050 | |||
| 1051 |
6/12✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 5236 times.
✓ Branch 4 → 3 taken 5236 times.
✗ Branch 4 → 5 not taken.
✓ Branch 5 → 6 taken 5236 times.
✗ Branch 5 → 7 not taken.
✓ Branch 7 → 8 taken 5236 times.
✗ Branch 7 → 9 not taken.
✗ Branch 9 → 10 not taken.
✓ Branch 9 → 11 taken 5236 times.
✗ Branch 11 → 12 not taken.
✓ Branch 11 → 13 taken 5236 times.
|
15708 | allocate (inds(size(bfun%bspline))) |
| 1052 |
2/2✓ Branch 13 → 14 taken 122807 times.
✓ Branch 13 → 15 taken 5236 times.
|
128043 | do i = 1, size(bfun%bspline) |
| 1053 | 128043 | inds(i) = i - 1 | |
| 1054 | end do | ||
| 1055 | |||
| 1056 |
1/2✗ Branch 17 → 18 not taken.
✓ Branch 17 → 19 taken 5236 times.
|
5236 | PetscCall(VecCreate(PETSC_COMM_SELF, vec, ierr)) |
| 1057 |
1/2✗ Branch 20 → 21 not taken.
✓ Branch 20 → 22 taken 5236 times.
|
5236 | PetscCall(VecSetSizes(vec, size(bfun%bspline), size(bfun%bspline), ierr)) |
| 1058 |
1/2✗ Branch 23 → 24 not taken.
✓ Branch 23 → 25 taken 5236 times.
|
5236 | PetscCall(VecSetFromOptions(vec, ierr)) |
| 1059 |
1/2✗ Branch 26 → 27 not taken.
✓ Branch 26 → 28 taken 5236 times.
|
5236 | PetscCall(VecSetValues(vec, size(bfun%bspline), inds, bfun%data, INSERT_VALUES, ierr)) |
| 1060 |
1/2✗ Branch 29 → 30 not taken.
✓ Branch 29 → 31 taken 5236 times.
|
5236 | PetscCall(VecAssemblyBegin(vec, ierr)) |
| 1061 |
1/2✗ Branch 32 → 33 not taken.
✓ Branch 32 → 34 taken 5236 times.
|
5236 | PetscCall(VecAssemblyEnd(vec, ierr)) |
| 1062 | |||
| 1063 |
1/2✗ Branch 34 → 35 not taken.
✓ Branch 34 → 36 taken 5236 times.
|
5236 | deallocate (inds) |
| 1064 | ✗ | end subroutine get_petsc_1dvec | |
| 1065 |
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 12846 times.
✗ Branch 4 → 5 not taken.
✗ Branch 4 → 6 not taken.
✓ Branch 8 → 9 taken 12846 times.
✓ Branch 8 → 17 taken 12846 times.
✗ Branch 10 → 11 not taken.
✓ Branch 10 → 12 taken 12846 times.
✓ Branch 12 → 13 taken 12846 times.
✗ Branch 12 → 16 not taken.
✓ Branch 13 → 14 taken 258 times.
✓ Branch 13 → 15 taken 12588 times.
|
25692 | end module m_bspline_basis |
| 1066 |