1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
|
module bsplines
use types, only: dp
use lapack, only: dgesv, dgbsv
use utils, only: stop_error
implicit none
private
public bspline, bspline_der, bspline_der2
contains
pure recursive function bspline(t, i, k, r) result(B)
! Returns values of a B-spline.
! There are set of `n` B-splines of order `k`. The mesh is composed of `n+k`
! knots, stored in the t(:) array.
real(dp), intent(in) :: t(:) ! {t_i}, i = 1, 2, ..., n+k
integer, intent(in) :: i ! i = 1..n ?
integer, intent(in) :: k ! Order of the spline k = 1, 2, 3, ...
real(dp), intent(in) :: r(:) ! points to evaluate the spline at
real(dp) :: B(size(r))
if (k == 1) then
if ((t(size(t)) - t(i+1)) < tiny(1._dp) .and. &
(t(i+1) - t(i)) > tiny(1._dp)) then
! If this is the last non-zero knot span, include the right end point,
! to ensure that the last basis function goes to 1.
where (t(i) <= r .and. r <= t(i+1))
B = 1
else where
B = 0
end where
else
! Otherwise exclude the right end point
where (t(i) <= r .and. r < t(i+1))
B = 1
else where
B = 0
end where
end if
else
B = 0
if (t(i+k-1)-t(i) > tiny(1._dp)) then
B = B + (r-t(i)) / (t(i+k-1)-t(i)) * bspline(t, i, k-1, r)
end if
if (t(i+k)-t(i+1) > tiny(1._dp)) then
B = B + (t(i+k)-r) / (t(i+k)-t(i+1)) * bspline(t, i+1, k-1, r)
end if
end if
end function
pure function bspline_der(t, i, k, r) result(B)
! Returns values of a derivative of a B-spline.
! There are set of `n` B-splines of order `k`. The mesh is composed of `n+k`
! knots, stored in the t(:) array.
real(dp), intent(in) :: t(:) ! {t_i}, i = 1, 2, ..., n+k
integer, intent(in) :: i ! i = 1..n ?
integer, intent(in) :: k ! Order of the spline k = 1, 2, 3, ...
real(dp), intent(in) :: r(:) ! points to evaluate the spline at
real(dp) :: B(size(r))
if (k == 1) then
B = 0
else
B = 0
if (t(i+k-1)-t(i) > tiny(1._dp)) then
B = B + (k-1) / (t(i+k-1)-t(i)) * bspline(t, i, k-1, r)
end if
if (t(i+k)-t(i+1) > tiny(1._dp)) then
B = B - (k-1) / (t(i+k)-t(i+1)) * bspline(t, i+1, k-1, r)
end if
end if
end function
pure function bspline_der2(t, i, k, r) result(B)
! Returns values of a second derivative of a B-spline.
! There are set of `n` B-splines of order `k`. The mesh is composed of `n+k`
! knots, stored in the t(:) array.
real(dp), intent(in) :: t(:) ! {t_i}, i = 1, 2, ..., n+k
integer, intent(in) :: i ! i = 1..n ?
integer, intent(in) :: k ! Order of the spline k = 1, 2, 3, ...
real(dp), intent(in) :: r(:) ! points to evaluate the spline at
real(dp) :: B(size(r))
if (k == 1) then
B = 0
else
B = 0
if (t(i+k-1)-t(i) > tiny(1._dp)) then
B = B + (k-1) / (t(i+k-1)-t(i)) * bspline_der(t, i, k-1, r)
end if
if (t(i+k)-t(i+1) > tiny(1._dp)) then
B = B - (k-1) / (t(i+k)-t(i+1)) * bspline_der(t, i+1, k-1, r)
end if
end if
end function
end module
|
