! -*- f90 -*- subroutine IGA_Basis_BSpline(k,uu,p,d,U,B) & bind(C, name="IGA_Basis_BSpline") use PetIGA implicit none integer(kind=IGA_INTEGER_KIND), intent(in),value :: k, p, d real (kind=IGA_REAL_KIND ), intent(in),value :: uu real (kind=IGA_REAL_KIND ), intent(in) :: U(0:k+p) real (kind=IGA_REAL_KIND ), intent(out) :: B(0:d,0:p) real (kind=IGA_REAL_KIND ) :: ders(0:p,0:d) call BasisFunsDers(k,uu,p,d,U,ders) B = transpose(ders) contains pure subroutine BasisFunsDers(i,uu,p,n,U,ders) use PetIGA implicit none integer(kind=IGA_INTEGER_KIND), intent(in) :: i, p, n real (kind=IGA_REAL_KIND ), intent(in) :: uu, U(0:i+p) real (kind=IGA_REAL_KIND ), intent(out) :: ders(0:p,0:n) integer(kind=IGA_INTEGER_KIND) :: j, k, r, s1, s2, rk, pk, j1, j2 real (kind=IGA_REAL_KIND ) :: saved, temp, d real (kind=IGA_REAL_KIND ) :: left(p), right(p) real (kind=IGA_REAL_KIND ) :: ndu(0:p,0:p), a(0:1,0:p) ndu(0,0) = 1 do j = 1, p left(j) = uu - U(i+1-j) right(j) = U(i+j) - uu saved = 0 do r = 0, j-1 ndu(j,r) = right(r+1) + left(j-r) temp = ndu(r,j-1) / ndu(j,r) ndu(r,j) = saved + right(r+1) * temp saved = left(j-r) * temp end do ndu(j,j) = saved end do ders(:,0) = ndu(:,p) do r = 0, p s1 = 0; s2 = 1; a(0,0) = 1 do k = 1, n d = 0 rk = r-k; pk = p-k; if (r >= k) then a(s2,0) = a(s1,0) / ndu(pk+1,rk) d = a(s2,0) * ndu(rk,pk) end if if (rk > -1) then j1 = 1 else j1 = -rk end if if (r-1 <= pk) then j2 = k-1 else j2 = p-r end if do j = j1, j2 a(s2,j) = (a(s1,j) - a(s1,j-1)) / ndu(pk+1,rk+j) d = d + a(s2,j) * ndu(rk+j,pk) end do if (r <= pk) then a(s2,k) = - a(s1,k-1) / ndu(pk+1,r) d = d + a(s2,k) * ndu(r,pk) end if ders(r,k) = d j = s1; s1 = s2; s2 = j; end do end do r = p do k = 1, n ders(:,k) = ders(:,k) * r r = r * (p-k) end do end subroutine BasisFunsDers end subroutine IGA_Basis_BSpline subroutine IGA_Basis_Lagrange(kk,uu,p,d,U,B) & bind(C, name="IGA_Basis_Lagrange") use PetIGA implicit none integer(kind=IGA_INTEGER_KIND), intent(in),value :: kk, p, d real (kind=IGA_REAL_KIND ), intent(in),value :: uu real (kind=IGA_REAL_KIND ), intent(in) :: U(0:kk+p) real (kind=IGA_REAL_KIND ), intent(out) :: B(0:d,0:p) integer(kind=IGA_INTEGER_KIND) :: m, i, j, k, l real (kind=IGA_REAL_KIND ) :: Lp, Ls1, Ls2, Ls3 real (kind=IGA_REAL_KIND ) :: X(0:p) do m = 0, p X(m) = U(kk) + m * (U(kk+1) - U(kk)) / p end do do m = 0, p Lp = 1 do i = 0, p if (i == m) cycle Lp = Lp * (uu-X(i))/(X(m)-X(i)) end do B(0,m) = Lp end do if (d < 1) return do m = 0, p Ls1 = 0 do j = 0, p if (j == m) cycle Lp = 1 do i = 0, p if (i == m .or. i == j) cycle Lp = Lp * (uu-X(i))/(X(m)-X(i)) end do Ls1 = Ls1 + Lp/(X(m)-X(j)) end do B(1,m) = Ls1 end do if (d < 2) return do m = 0, p Ls2 = 0 do k = 0, p if (k == m) cycle Ls1 = 0 do j = 0, p if (j == m .or. j == k) cycle Lp = 1 do i = 0, p if (i == m .or. i == k .or. i == j) cycle Lp = Lp * (uu-X(i))/(X(m)-X(i)) end do Ls1 = Ls1 + Lp/(X(m)-X(j)) end do Ls2 = Ls2 + Ls1/(X(m)-X(k)) end do B(2,m) = Ls2 end do if (d < 3) return do m = 0, p Ls3 = 0 do l = 0, p if (l == m) cycle Ls2 = 0 do k = 0, p if (k == m .or. k == l) cycle Ls1 = 0 do j = 0, p if (j == m .or. j == l .or. j == k) cycle Lp = 1 do i = 0, p if (i == m .or. i == l .or. i == k .or. i == j) cycle Lp = Lp * (uu-X(i))/(X(m)-X(i)) end do Ls1 = Ls1 + Lp/(X(m)-X(j)) end do Ls2 = Ls2 + Ls1/(X(m)-X(k)) end do Ls3 = Ls3 + Ls2/(X(m)-X(l)) end do B(3,m) = Ls3 end do end subroutine IGA_Basis_Lagrange subroutine IGA_Basis_Hierarchical(kk,uu,p,d,U,B) & bind(C, name="IGA_Basis_Hierarchical") use PetIGA implicit none integer(kind=IGA_INTEGER_KIND), intent(in),value :: kk, p, d real (kind=IGA_REAL_KIND ), intent(in),value :: uu real (kind=IGA_REAL_KIND ), intent(in) :: U(0:kk+p) real (kind=IGA_REAL_KIND ), intent(out) :: B(0:d,0:p) integer(kind=IGA_INTEGER_KIND) :: i, k real (kind=IGA_REAL_KIND ) :: J, x, Lp(0:p,0:d) real (kind=IGA_REAL_KIND ), parameter :: two = 2 J = (U(kk+1)-U(kk))/2 x = (uu-U(kk))/J - 1 B(0,0) = (1-x)/2 B(0,p) = (x+1)/2 if (d > 0) then B(1,0) = -1/two B(1,p) = +1/two endif if (p > 1) then Lp(:,:) = 0 Lp(0,0) = 1 Lp(1,0) = x if (d > 0) then Lp(0,1) = 0 Lp(1,1) = 1 end if do i = 1, p-1 Lp(i+1,0) = ((2*i+1)*x*Lp(i,0) - i*Lp(i-1,0))/(i+1) B(0,i) = (-Lp(i+1,0) + Lp(i-1,0))/sqrt(two*(2*(i+1)-1)) do k = 1, d Lp(i+1,k) = (2*i+1)*Lp(i,k-1) + Lp(i-1,k) B(k,i) = (-Lp(i+1,k) + Lp(i-1,k))/sqrt(two*(2*(i+1)-1)) end do end do end if do k = 1, d B(k,:) = B(k,:)/(J**k) end do end subroutine IGA_Basis_Hierarchical