Skip to content

Commit

Permalink
Add is_rational() method.
Browse files Browse the repository at this point in the history
  • Loading branch information
gha3mi committed Apr 7, 2024
1 parent 4d8f6a7 commit a3cb1d9
Show file tree
Hide file tree
Showing 3 changed files with 73 additions and 18 deletions.
35 changes: 27 additions & 8 deletions src/forcad_nurbs_curve.f90
Original file line number Diff line number Diff line change
Expand Up @@ -50,6 +50,7 @@ module forcad_nurbs_curve
procedure :: elevate_degree !!> Elevate the degree of the curve
procedure :: derivative !!> Compute the derivative of the NURBS curve
procedure :: basis !!> Compute the basis functions of the NURBS curve
procedure :: is_rational !!> Check if the NURBS curve is rational
end type
!===============================================================================

Expand Down Expand Up @@ -177,15 +178,15 @@ pure subroutine create(this, res, Xt)
if (allocated(this%Xg)) deallocate(this%Xg)
allocate(this%Xg(this%ng, size(this%Xc,2)))

if (allocated(this%Wc)) then
if (this%is_rational()) then ! NURBS
do i = 1, size(this%Xt, 1)
Tgc = basis_bspline(this%Xt(i), this%knot, this%nc, this%degree)
Tgc = Tgc*(this%Wc/(dot_product(Tgc,this%Wc)))
do j = 1, size(this%Xc, 2)
this%Xg(i,j) = dot_product(Tgc,this%Xc(:,j))
end do
end do
else
else ! B-Spline
do i = 1, size(this%Xt, 1)
Tgc = basis_bspline(this%Xt(i), this%knot, this%nc, this%degree)
do j = 1, size(this%Xc, 2)
Expand Down Expand Up @@ -533,7 +534,7 @@ pure subroutine insert_knots(this,Xth,r)
integer :: k, i, s, dim, j, n_new
real(rk), allocatable :: Xcw(:,:), Xcw_new(:,:), Xc_new(:,:), Wc_new(:), knot_new(:)

if (allocated(this%Wc)) then ! NURBS
if (this%is_rational()) then ! NURBS

do i = 1, size(Xth)
k = findspan(this%nc-1,this%degree,Xth(i),this%knot)
Expand Down Expand Up @@ -616,7 +617,7 @@ pure subroutine elevate_degree(this, t)
real(rk), allocatable :: Xcw(:,:), Xcw_new(:,:), knot_new(:), Xc_new(:,:), Wc_new(:)
integer :: dim, j, nc_new

if (allocated(this%Wc)) then ! NURBS
if (this%is_rational()) then ! NURBS

dim = size(this%Xc,2)
allocate(Xcw(size(this%Xc,1),dim+1))
Expand Down Expand Up @@ -679,13 +680,13 @@ pure subroutine derivative(this, res, Xt, dTgc)

allocate(dTgc(size(this%Xt, 1), this%nc))

if (allocated(this%Wc)) then
if (this%is_rational()) then ! NURBS
do i = 1, size(this%Xt, 1)
dTgci = basis_bspline_der(this%Xt(i), this%knot, this%nc, this%degree)
dTgci = dTgci*(this%Wc/(dot_product(dTgci,this%Wc)))
dTgc(i,:) = dTgci
end do
else
else ! B-Spline
do i = 1, size(this%Xt, 1)
dTgci = basis_bspline_der(this%Xt(i), this%knot, this%nc, this%degree)
dTgc(i,:) = dTgci
Expand Down Expand Up @@ -720,13 +721,13 @@ pure subroutine basis(this, res, Xt, Tgc)

allocate(Tgc(size(this%Xt, 1), this%nc))

if (allocated(this%Wc)) then
if (this%is_rational()) then ! NURBS
do i = 1, size(this%Xt, 1)
Tgci = basis_bspline(this%Xt(i), this%knot, this%nc, this%degree)
Tgci = Tgci*(this%Wc/(dot_product(Tgci,this%Wc)))
Tgc(i,:) = Tgci
end do
else
else ! B-Spline
do i = 1, size(this%Xt, 1)
Tgci = basis_bspline(this%Xt(i), this%knot, this%nc, this%degree)
Tgc(i,:) = Tgci
Expand All @@ -735,4 +736,22 @@ pure subroutine basis(this, res, Xt, Tgc)
end subroutine
!===============================================================================


!===============================================================================
!> author: Seyed Ali Ghasemi
!> license: BSD 3-Clause
pure function is_rational(this) result(r)
class(nurbs_curve), intent(in) :: this
logical :: r

r = .false.
if (allocated(this%Wc)) then
if (any(this%Wc /= this%Wc(1))) then
r = .true.
end if
end if
end function
!===============================================================================


end module forcad_nurbs_curve
38 changes: 28 additions & 10 deletions src/forcad_nurbs_surface.f90
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,7 @@ module forcad_nurbs_surface
procedure :: basis !!> Compute the basis functions of the NURBS surface
procedure :: insert_knots !!> Insert knots into the knot vector
procedure :: elevate_degree !!> Elevate degree
procedure :: is_rational !!> Check if the NURBS surface is rational
end type
!===============================================================================

Expand Down Expand Up @@ -190,7 +191,7 @@ pure subroutine create(this, res1, res2, Xt1, Xt2)
if (allocated(this%Xg)) deallocate(this%Xg)
allocate(this%Xg(this%ng(1)*this%ng(2), size(this%Xc,2)))

if (allocated(this%Wc)) then ! NURBS surface
if (this%is_rational()) then ! NURBS
do i = 1, size(Xt, 1)
Tgc1 = basis_bspline(Xt(i,1), this%knot1, this%nc(1), this%degree(1))
Tgc2 = basis_bspline(Xt(i,2), this%knot2, this%nc(2), this%degree(2))
Expand All @@ -200,7 +201,7 @@ pure subroutine create(this, res1, res2, Xt1, Xt2)
this%Xg(i,j) = dot_product(Tgc,this%Xc(:,j))
end do
end do
else
else ! B-Spline
do i = 1, size(Xt, 1)
Tgc1 = basis_bspline(Xt(i,1), this%knot1, this%nc(1), this%degree(1))
Tgc2 = basis_bspline(Xt(i,2), this%knot2, this%nc(2), this%degree(2))
Expand Down Expand Up @@ -669,15 +670,15 @@ pure subroutine derivative(this, res1, res2, Xt1, Xt2, dTgc)

allocate(dTgc(this%ng(1)*this%ng(2), this%nc(1)*this%nc(2)))

if (allocated(this%Wc)) then ! NURBS surface
if (this%is_rational()) then ! NURBS
do i = 1, size(Xt, 1)
dTgc1 = basis_bspline_der(Xt(i,1), this%knot1, this%nc(1), this%degree(1))
dTgc2 = basis_bspline_der(Xt(i,2), this%knot2, this%nc(2), this%degree(2))
dTgci = kron(dTgc2, dTgc1)
dTgci = dTgci*(this%Wc/(dot_product(dTgci,this%Wc)))
dTgc(i,:) = dTgci
end do
else
else ! B-Spline
do i = 1, size(Xt, 1)
dTgc1 = basis_bspline_der(Xt(i,1), this%knot1, this%nc(1), this%degree(1))
dTgc2 = basis_bspline_der(Xt(i,2), this%knot2, this%nc(2), this%degree(2))
Expand Down Expand Up @@ -734,15 +735,15 @@ pure subroutine basis(this, res1, res2, Xt1, Xt2, Tgc)

allocate(Tgc(this%ng(1)*this%ng(2), this%nc(1)*this%nc(2)))

if (allocated(this%Wc)) then ! NURBS surface
if (this%is_rational()) then ! NURBS
do i = 1, size(Xt, 1)
Tgc1 = basis_bspline(Xt(i,1), this%knot1, this%nc(1), this%degree(1))
Tgc2 = basis_bspline(Xt(i,2), this%knot2, this%nc(2), this%degree(2))
Tgci = kron(Tgc2, Tgc1)
Tgci = Tgci*(this%Wc/(dot_product(Tgci,this%Wc)))
Tgc(i,:) = Tgci
end do
else
else ! B-Spline
do i = 1, size(Xt, 1)
Tgc1 = basis_bspline(Xt(i,1), this%knot1, this%nc(1), this%degree(1))
Tgc2 = basis_bspline(Xt(i,2), this%knot2, this%nc(2), this%degree(2))
Expand All @@ -769,7 +770,7 @@ pure subroutine insert_knots(this, dir ,Xth,r)

if (dir == 1) then ! direction 1

if (allocated(this%Wc)) then ! NURBS
if(this%is_rational()) then ! NURBS

do i = 1, size(Xth)
k = findspan(this%nc(1)-1,this%degree(1),Xth(i),this%knot1)
Expand Down Expand Up @@ -853,7 +854,7 @@ pure subroutine insert_knots(this, dir ,Xth,r)

elseif (dir == 2) then! direction 2

if (allocated(this%Wc)) then ! NURBS
if(this%is_rational()) then ! NURBS

do i = 1, size(Xth)
k = findspan(this%nc(2)-1,this%degree(2),Xth(i),this%knot2)
Expand Down Expand Up @@ -966,7 +967,7 @@ pure subroutine elevate_degree(this, dir, t)

if (dir == 1) then ! direction 1

if (allocated(this%Wc)) then ! NURBS
if(this%is_rational()) then ! NURBS

dim = size(this%Xc,2)
allocate(Xcw(size(this%Xc,1),dim+1))
Expand Down Expand Up @@ -1010,7 +1011,7 @@ pure subroutine elevate_degree(this, dir, t)

elseif (dir == 2) then ! direction 2

if (allocated(this%Wc)) then ! NURBS
if(this%is_rational()) then ! NURBS

dim = size(this%Xc,2)
allocate(Xcw(size(this%Xc,1),dim+1))
Expand Down Expand Up @@ -1069,4 +1070,21 @@ pure subroutine elevate_degree(this, dir, t)
end subroutine
!===============================================================================


!===============================================================================
!> author: Seyed Ali Ghasemi
!> license: BSD 3-Clause
pure function is_rational(this) result(r)
class(nurbs_surface), intent(in) :: this
logical :: r

r = .false.
if(allocated(this%Wc)) then
if (any(this%Wc /= this%Wc(1))) then
r = .true.
end if
end if
end function
!===============================================================================

end module forcad_nurbs_surface
18 changes: 18 additions & 0 deletions src/forcad_nurbs_volume.f90
Original file line number Diff line number Diff line change
Expand Up @@ -54,6 +54,7 @@ module forcad_nurbs_volume
procedure :: basis !!> Compute the basis functions of the NURBS volume
procedure :: insert_knots !!> Insert knots into the knot vector
procedure :: elevate_degree !!> Elevate the degree of the NURBS volume
procedure :: is_rational !!> Check if the NURBS volume is rational
end type
!===============================================================================

Expand Down Expand Up @@ -1316,4 +1317,21 @@ pure subroutine elevate_degree(this, dir, t)
end subroutine
!===============================================================================


!===============================================================================
!> author: Seyed Ali Ghasemi
!> license: BSD 3-Clause
pure function is_rational(this) result(r)
class(nurbs_volume), intent(in) :: this
logical :: r

r = .false.
if (allocated(this%Wc)) then
if (any(this%Wc /= this%Wc(1))) then
r = .true.
end if
end if
end function
!===============================================================================

end module forcad_nurbs_volume

0 comments on commit a3cb1d9

Please sign in to comment.