PetIGA / src / petigabsp.f90

Lisandro Dalcin f3d7eb8 

Lisandro Dalcin ccbe5f2 
Lisandro Dalcin 87105e5 
Lisandro Dalcin a9da11f 
Lisandro Dalcin f3d7eb8 
Lisandro Dalcin ccbe5f2 
Lisandro Dalcin ca3d144 
Lisandro Dalcin ccbe5f2 

Lisandro Dalcin ca3d144 
Lisandro Dalcin ccbe5f2 

Lisandro Dalcin 87105e5 

Lisandro Dalcin a9da11f 

Lisandro Dalcin ca3d144 






Lisandro Dalcin 791c6f2 
Lisandro Dalcin f3d7eb8 


Lisandro Dalcin 791c6f2 
Lisandro Dalcin f3d7eb8 










Lisandro Dalcin 791c6f2 
Lisandro Dalcin f3d7eb8 
Lisandro Dalcin 791c6f2 
Lisandro Dalcin f3d7eb8 











Lisandro Dalcin 87105e5 
Lisandro Dalcin f3d7eb8 


















Lisandro Dalcin 87105e5 

Lisandro Dalcin ccbe5f2 













Lisandro Dalcin a11a4a6 


Lisandro Dalcin ccbe5f2 

Lisandro Dalcin 791c6f2 
Lisandro Dalcin ccbe5f2 








Lisandro Dalcin 791c6f2 
Lisandro Dalcin ccbe5f2 

Lisandro Dalcin 791c6f2 
Lisandro Dalcin ccbe5f2 










Lisandro Dalcin 791c6f2 
Lisandro Dalcin ccbe5f2 

Lisandro Dalcin 791c6f2 
Lisandro Dalcin ccbe5f2 

Lisandro Dalcin 791c6f2 
Lisandro Dalcin ccbe5f2 












Lisandro Dalcin 791c6f2 
Lisandro Dalcin ccbe5f2 

Lisandro Dalcin 791c6f2 
Lisandro Dalcin ccbe5f2 

Lisandro Dalcin 791c6f2 
Lisandro Dalcin ccbe5f2 

Lisandro Dalcin 791c6f2 
Lisandro Dalcin ccbe5f2 













Lisandro Dalcin 1df08c8 











Lisandro Dalcin 791c6f2 
Lisandro Dalcin 1df08c8 
Lisandro Dalcin 791c6f2 

Lisandro Dalcin 1df08c8 
Lisandro Dalcin 791c6f2 

Lisandro Dalcin 1df08c8 
Lisandro Dalcin 791c6f2 

Lisandro Dalcin 1df08c8 


Lisandro Dalcin 791c6f2 

Lisandro Dalcin 1df08c8 

Lisandro Dalcin 791c6f2 

Lisandro Dalcin 1df08c8 














! -*- 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
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.