From 3adf21c3293ce273c3f5a9ace1c60462ca207f3e Mon Sep 17 00:00:00 2001 From: Seyed Ali Ghasemi Date: Tue, 16 Apr 2024 14:25:36 +0200 Subject: [PATCH] Improve set1() method --- src/forcad_nurbs_surface.f90 | 52 +++++++++++++++------ src/forcad_nurbs_volume.f90 | 89 ++++++++++++++++++++++++++++-------- 2 files changed, 109 insertions(+), 32 deletions(-) diff --git a/src/forcad_nurbs_surface.f90 b/src/forcad_nurbs_surface.f90 index 6aa280829..a67ae3ac1 100644 --- a/src/forcad_nurbs_surface.f90 +++ b/src/forcad_nurbs_surface.f90 @@ -93,13 +93,24 @@ pure subroutine set1(this, knot1, knot2, Xc, Wc) real(rk), intent(in), contiguous :: Xc(:,:) real(rk), intent(in), contiguous, optional :: Wc(:) + if (allocated(this%knot1)) deallocate(this%knot1) + if (allocated(this%knot2)) deallocate(this%knot2) + if (allocated(this%Xc)) deallocate(this%Xc) + this%knot1 = knot1 this%knot2 = knot2 this%degree = this%get_degree() this%nc(1) = this%get_nc(1) this%nc(2) = this%get_nc(2) this%Xc = Xc - if (present(Wc)) this%Wc = Wc + if (present(Wc)) then + if (size(Wc) /= this%nc(1)*this%nc(2)) then + error stop 'Number of weights does not match the number of control points.' + else + if (allocated(this%Wc)) deallocate(this%Wc) + this%Wc = Wc + end if + end if end subroutine !=============================================================================== @@ -880,6 +891,7 @@ pure subroutine insert_knots(this, dir ,Xth,r) integer :: k, i, s, dim, j, n_new real(rk), allocatable :: Xc(:,:), Xcw(:,:), Xcw_new(:,:), Xc_new(:,:), Wc_new(:), knot_new(:) real(rk), allocatable:: Xc3(:,:,:) + real(rk), allocatable :: knot1(:), knot2(:) if (dir == 1) then ! direction 1 @@ -924,7 +936,8 @@ pure subroutine insert_knots(this, dir ,Xth,r) end do Wc_new(:) = Xcw_new(:,dim+1) - call this%set(knot1=knot_new, knot2=this%knot2, Xc=Xc_new, Wc=Wc_new) + knot2 = this%knot2 + call this%set(knot1=knot_new, knot2=knot2, Xc=Xc_new, Wc=Wc_new) deallocate(Xcw, Xcw_new, Xc_new, Wc_new) end do @@ -957,7 +970,8 @@ pure subroutine insert_knots(this, dir ,Xth,r) Xc_new = reshape(Xc_new,[(this%nc(2))*(n_new+1),dim]) - call this%set(knot1=knot_new, knot2=this%knot2, Xc=Xc_new) + knot2 = this%knot2 + call this%set(knot1=knot_new, knot2=knot2, Xc=Xc_new) end do end if @@ -1009,7 +1023,8 @@ pure subroutine insert_knots(this, dir ,Xth,r) end do Wc_new(:) = Xcw_new(:,dim+1) - call this%set(knot2=knot_new, knot1=this%knot1, Xc=Xc_new, Wc=Wc_new) + knot1 = this%knot1 + call this%set(knot2=knot_new, knot1=knot1, Xc=Xc_new, Wc=Wc_new) deallocate(Xcw, Xcw_new, Xc_new, Wc_new) end do @@ -1045,7 +1060,8 @@ pure subroutine insert_knots(this, dir ,Xth,r) Xc3 = reshape(Xc3, [this%nc(1),n_new+1,dim], order=[2,1,3]) Xc_new = reshape(Xc3,[(this%nc(1))*(n_new+1),dim]) - call this%set(knot2=knot_new, knot1=this%knot1, Xc=Xc_new) + knot1 = this%knot1 + call this%set(knot2=knot_new, knot1=knot1, Xc=Xc_new) end do @@ -1069,6 +1085,7 @@ pure subroutine elevate_degree(this, dir, t) real(rk), allocatable :: Xc(:,:), Xcw(:,:), Xcw_new(:,:), knot_new(:), Xc_new(:,:), Wc_new(:) integer :: dim, j, nc_new real(rk), allocatable:: Xc3(:,:,:) + real(rk), allocatable :: knot1(:), knot2(:) if (dir == 1) then ! direction 1 @@ -1096,7 +1113,8 @@ pure subroutine elevate_degree(this, dir, t) Wc_new(:) = Xcw_new(:,dim+1) - call this%set(knot1=knot_new, knot2=this%knot2, Xc=Xc_new, Wc=Wc_new) + knot2 = this%knot2 + call this%set(knot1=knot_new, knot2=knot2, Xc=Xc_new, Wc=Wc_new) deallocate(Xcw, Xcw_new, Xc_new, Wc_new) else ! B-Spline @@ -1108,7 +1126,8 @@ pure subroutine elevate_degree(this, dir, t) Xc_new = reshape(Xc_new,[this%nc(2)*nc_new,dim],order=[1,2]) - call this%set(knot1=knot_new, knot2=this%knot2, Xc=Xc_new) + knot2 = this%knot2 + call this%set(knot1=knot_new, knot2=knot2, Xc=Xc_new) deallocate(Xc, Xc_new) end if @@ -1143,7 +1162,8 @@ pure subroutine elevate_degree(this, dir, t) Wc_new(:) = Xcw_new(:,dim+1) - call this%set(knot2=knot_new, knot1=this%knot1, Xc=Xc_new, Wc=Wc_new) + knot1 = this%knot1 + call this%set(knot2=knot_new, knot1=knot1, Xc=Xc_new, Wc=Wc_new) deallocate(Xcw, Xcw_new, Xc_new, Wc_new) else ! B-Spline @@ -1160,7 +1180,8 @@ pure subroutine elevate_degree(this, dir, t) Xc3 = reshape(Xc3, [this%nc(1),nc_new,dim], order=[2,1,3]) Xc_new = reshape(Xc3,[(this%nc(1))*nc_new,dim]) - call this%set(knot2=knot_new, knot1=this%knot1, Xc=Xc_new) + knot1 = this%knot1 + call this%set(knot2=knot_new, knot1=knot1, Xc=Xc_new) end if @@ -1275,6 +1296,7 @@ pure subroutine remove_knots(this, dir ,Xth,r) integer :: k, i, s, dim, j, nc_new, t real(rk), allocatable :: Xc(:,:), Xcw(:,:), Xcw_new(:,:), Xc_new(:,:), Wc_new(:), knot_new(:) real(rk), allocatable:: Xc3(:,:,:) + real(rk), allocatable :: knot1(:), knot2(:) if (dir == 1) then ! direction 1 @@ -1327,7 +1349,8 @@ pure subroutine remove_knots(this, dir ,Xth,r) Wc_new(:) = Xcw_new(:,dim+1) - call this%set(knot1=knot_new, knot2=this%knot2, Xc=Xc_new, Wc=Wc_new) + knot2 = this%knot2 + call this%set(knot1=knot_new, knot2=knot2, Xc=Xc_new, Wc=Wc_new) deallocate(Xcw_new, Xc_new, Wc_new) end if end do @@ -1368,7 +1391,8 @@ pure subroutine remove_knots(this, dir ,Xth,r) nc_new = size(Xc_new,1) Xc_new = reshape(Xc_new,[(this%nc(2))*(nc_new),dim],order=[1,2]) - call this%set(knot1=knot_new, knot2=this%knot2, Xc=Xc_new) + knot2 = this%knot2 + call this%set(knot1=knot_new, knot2=knot2, Xc=Xc_new) end if end do @@ -1431,7 +1455,8 @@ pure subroutine remove_knots(this, dir ,Xth,r) Wc_new(:) = Xcw_new(:,dim+1) - call this%set(knot2=knot_new, knot1=this%knot1, Xc=Xc_new, Wc=Wc_new) + knot1 = this%knot1 + call this%set(knot2=knot_new, knot1=knot1, Xc=Xc_new, Wc=Wc_new) deallocate(Xcw_new, Xc_new, Wc_new) end if @@ -1477,7 +1502,8 @@ pure subroutine remove_knots(this, dir ,Xth,r) Xc3 = reshape(Xc3, [this%nc(1),nc_new,dim], order=[2,1,3]) Xc_new = reshape(Xc3,[(this%nc(1))*(nc_new),dim]) - call this%set(knot2=knot_new, knot1=this%knot1, Xc=Xc_new) + knot1 = this%knot1 + call this%set(knot2=knot_new, knot1=knot1, Xc=Xc_new) end if end do diff --git a/src/forcad_nurbs_volume.f90 b/src/forcad_nurbs_volume.f90 index 1cddfa4c3..50fc9eb6b 100644 --- a/src/forcad_nurbs_volume.f90 +++ b/src/forcad_nurbs_volume.f90 @@ -95,6 +95,11 @@ pure subroutine set1(this, knot1, knot2, knot3, Xc, Wc) real(rk), intent(in), contiguous :: Xc(:,:) real(rk), intent(in), contiguous, optional :: Wc(:) + if (allocated(this%knot1)) deallocate(this%knot1) + if (allocated(this%knot2)) deallocate(this%knot2) + if (allocated(this%knot3)) deallocate(this%knot3) + if (allocated(this%Xc)) deallocate(this%Xc) + this%knot1 = knot1 this%knot2 = knot2 this%knot3 = knot3 @@ -103,7 +108,14 @@ pure subroutine set1(this, knot1, knot2, knot3, Xc, Wc) this%nc(2) = this%get_nc(2) this%nc(3) = this%get_nc(3) this%Xc = Xc - if (present(Wc)) this%Wc = Wc + if (present(Wc)) then + if (size(Wc) /= this%nc(1)*this%nc(2)*this%nc(3)) then + error stop 'Number of weights does not match the number of control points.' + else + if (allocated(this%Wc)) deallocate(this%Wc) + this%Wc = Wc + end if + end if end subroutine !=============================================================================== @@ -995,6 +1007,7 @@ pure subroutine insert_knots(this, dir ,Xth,r) integer :: k, i, s, dim, j, n_new real(rk), allocatable :: Xc(:,:), Xcw(:,:), Xcw_new(:,:), Xc_new(:,:), Wc_new(:), knot_new(:) real(rk), allocatable :: Xc4(:,:,:,:) + real(rk), allocatable :: knot1(:), knot2(:), knot3(:) if (dir == 1) then ! direction 1 @@ -1039,7 +1052,9 @@ pure subroutine insert_knots(this, dir ,Xth,r) end do Wc_new(:) = Xcw_new(:,dim+1) - call this%set(knot1=knot_new, knot2=this%knot2, knot3=this%knot3, Xc=Xc_new, Wc=Wc_new) + knot2 = this%knot2 + knot3 = this%knot3 + call this%set(knot1=knot_new, knot2=knot2, knot3=knot3, Xc=Xc_new, Wc=Wc_new) deallocate(Xcw, Xcw_new, Xc_new, Wc_new) end do @@ -1072,7 +1087,9 @@ pure subroutine insert_knots(this, dir ,Xth,r) Xc_new = reshape(Xc_new,[(n_new+1)*this%nc(2)*this%nc(3),dim]) - call this%set(knot1=knot_new, knot2=this%knot2, knot3=this%knot3, Xc=Xc_new) + knot2 = this%knot2 + knot3 = this%knot3 + call this%set(knot1=knot_new, knot2=knot2, knot3=knot3, Xc=Xc_new) end do end if @@ -1123,7 +1140,9 @@ pure subroutine insert_knots(this, dir ,Xth,r) end do Wc_new(:) = Xcw_new(:,dim+1) - call this%set(knot1=this%knot1, knot2=knot_new, knot3=this%knot3, Xc=Xc_new, Wc=Wc_new) + knot1 = this%knot1 + knot3 = this%knot3 + call this%set(knot1=knot1, knot2=knot_new, knot3=knot3, Xc=Xc_new, Wc=Wc_new) deallocate(Xcw, Xcw_new, Xc_new, Wc_new) end do @@ -1159,7 +1178,9 @@ pure subroutine insert_knots(this, dir ,Xth,r) Xc4 = reshape(Xc4, [this%nc(1),n_new+1,this%nc(3),dim], order=[2,1,3,4]) Xc_new = reshape(Xc4,[this%nc(1)*(n_new+1)*this%nc(3),dim]) - call this%set(knot1=this%knot1, knot2=knot_new, knot3= this%knot3, Xc=Xc_new) + knot1 = this%knot1 + knot3 = this%knot3 + call this%set(knot1=knot1, knot2=knot_new, knot3=knot3, Xc=Xc_new) end do @@ -1211,7 +1232,9 @@ pure subroutine insert_knots(this, dir ,Xth,r) end do Wc_new(:) = Xcw_new(:,dim+1) - call this%set(knot1=this%knot1, knot2=this%knot2, knot3=knot_new, Xc=Xc_new, Wc=Wc_new) + knot1 = this%knot1 + knot2 = this%knot2 + call this%set(knot1=knot1, knot2=knot2, knot3=knot_new, Xc=Xc_new, Wc=Wc_new) deallocate(Xcw, Xcw_new, Xc_new, Wc_new) end do @@ -1247,7 +1270,9 @@ pure subroutine insert_knots(this, dir ,Xth,r) Xc4 = reshape(Xc4, [this%nc(1),this%nc(2),n_new+1,dim], order=[3,2,1,4]) Xc_new = reshape(Xc4, [this%nc(1)*this%nc(2)*(n_new+1),dim]) - call this%set(knot1=this%knot1, knot2=this%knot2, knot3=knot_new, Xc=Xc_new) + knot1 = this%knot1 + knot2 = this%knot2 + call this%set(knot1=knot1, knot2=knot2, knot3=knot_new, Xc=Xc_new) end do end if @@ -1270,6 +1295,7 @@ pure subroutine elevate_degree(this, dir, t) real(rk), allocatable :: Xc(:,:), Xcw(:,:), Xcw_new(:,:), Xc_new(:,:), Wc_new(:), knot_new(:) integer :: nc_new, dim, j real(rk), allocatable:: Xc4(:,:,:,:) + real(rk), allocatable :: knot1(:), knot2(:), knot3(:) if (dir == 1) then ! direction 1 @@ -1296,7 +1322,9 @@ pure subroutine elevate_degree(this, dir, t) end do Wc_new(:) = Xcw_new(:,dim+1) - call this%set(knot1=knot_new, knot2=this%knot2, knot3=this%knot3, Xc=Xc_new, Wc=Wc_new) + knot2 = this%knot2 + knot3 = this%knot3 + call this%set(knot1=knot_new, knot2=knot2, knot3=knot3, Xc=Xc_new, Wc=Wc_new) deallocate(Xcw, Xcw_new, Xc_new, Wc_new) else ! B-Spline @@ -1309,7 +1337,9 @@ pure subroutine elevate_degree(this, dir, t) Xc_new = reshape(Xc_new,[nc_new*this%nc(2)*this%nc(3),dim],order=[1,2]) - call this%set(knot1=knot_new, knot2=this%knot2, knot3=this%knot3, Xc=Xc_new) + knot2 = this%knot2 + knot3 = this%knot3 + call this%set(knot1=knot_new, knot2=knot2, knot3=knot3, Xc=Xc_new) end if @@ -1345,7 +1375,9 @@ pure subroutine elevate_degree(this, dir, t) Wc_new(:) = Xcw_new(:,dim+1) - call this%set(knot1=this%knot1, knot2=knot_new, knot3=this%knot3, Xc=Xc_new, Wc=Wc_new) + knot1 = this%knot1 + knot3 = this%knot3 + call this%set(knot1=knot1, knot2=knot_new, knot3=knot3, Xc=Xc_new, Wc=Wc_new) deallocate(Xcw, Xcw_new, Xc_new, Wc_new) else ! B-Spline @@ -1362,7 +1394,9 @@ pure subroutine elevate_degree(this, dir, t) Xc4 = reshape(Xc4, [this%nc(1),nc_new,this%nc(3),dim], order=[2,1,3,4]) Xc_new = reshape(Xc4,[this%nc(1)*nc_new*this%nc(3),dim]) - call this%set(knot1=this%knot1, knot2=knot_new, knot3= this%knot3, Xc=Xc_new) + knot1 = this%knot1 + knot3 = this%knot3 + call this%set(knot1=knot1, knot2=knot_new, knot3=knot3, Xc=Xc_new) end if @@ -1396,7 +1430,9 @@ pure subroutine elevate_degree(this, dir, t) Wc_new(:) = Xcw_new(:,dim+1) - call this%set(knot1=this%knot1, knot2=this%knot2, knot3=knot_new, Xc=Xc_new, Wc=Wc_new) + knot1 = this%knot1 + knot2 = this%knot2 + call this%set(knot1=knot1, knot2=knot2, knot3=knot_new, Xc=Xc_new, Wc=Wc_new) deallocate(Xcw, Xcw_new, Xc_new, Wc_new) else ! B-Spline @@ -1413,7 +1449,9 @@ pure subroutine elevate_degree(this, dir, t) Xc4 = reshape(Xc4, [this%nc(1),this%nc(2),nc_new,dim], order=[3,2,1,4]) Xc_new = reshape(Xc4,[this%nc(1)*this%nc(2)*nc_new,dim]) - call this%set(knot1=this%knot1, knot2=this%knot2, knot3=knot_new, Xc=Xc_new) + knot1 = this%knot1 + knot2 = this%knot2 + call this%set(knot1=knot1, knot2=knot2, knot3=knot_new, Xc=Xc_new) end if @@ -1605,6 +1643,7 @@ pure subroutine remove_knots(this, dir ,Xth,r) integer :: k, i, s, dim, j, nc_new, t real(rk), allocatable :: Xc(:,:), Xcw(:,:), Xcw_new(:,:), Xc_new(:,:), Wc_new(:), knot_new(:) real(rk), allocatable :: Xc4(:,:,:,:) + real(rk), allocatable :: knot1(:), knot2(:), knot3(:) if (dir == 1) then ! direction 1 @@ -1656,7 +1695,9 @@ pure subroutine remove_knots(this, dir ,Xth,r) end do Wc_new(:) = Xcw_new(:,dim+1) - call this%set(knot1=knot_new, knot2=this%knot2, knot3=this%knot3, Xc=Xc_new, Wc=Wc_new) + knot2 = this%knot2 + knot3 = this%knot3 + call this%set(knot1=knot_new, knot2=knot2, knot3=knot3, Xc=Xc_new, Wc=Wc_new) deallocate(Xcw_new, Xc_new, Wc_new) end if end do @@ -1697,7 +1738,9 @@ pure subroutine remove_knots(this, dir ,Xth,r) nc_new = size(Xcw_new,1) Xc_new = reshape(Xc_new,[(nc_new)*this%nc(2)*this%nc(3),dim]) - call this%set(knot1=knot_new, knot2=this%knot2, knot3=this%knot3, Xc=Xc_new) + knot2 = this%knot2 + knot3 = this%knot3 + call this%set(knot1=knot_new, knot2=knot2, knot3=knot3, Xc=Xc_new) end if end do @@ -1759,7 +1802,9 @@ pure subroutine remove_knots(this, dir ,Xth,r) Wc_new(:) = Xcw_new(:,dim+1) - call this%set(knot1=this%knot1, knot2=knot_new, knot3=this%knot3, Xc=Xc_new, Wc=Wc_new) + knot1 = this%knot1 + knot3 = this%knot3 + call this%set(knot1=knot1, knot2=knot_new, knot3=knot3, Xc=Xc_new, Wc=Wc_new) deallocate(Xcw_new, Xc_new, Wc_new) end if end do @@ -1804,7 +1849,9 @@ pure subroutine remove_knots(this, dir ,Xth,r) Xc4 = reshape(Xc4, [this%nc(1),nc_new,this%nc(3),dim], order=[2,1,3,4]) Xc_new = reshape(Xc4, [this%nc(1)*(nc_new)*this%nc(3),dim]) - call this%set(knot1=this%knot1, knot2=knot_new, knot3= this%knot3, Xc=Xc_new) + knot1 = this%knot1 + knot3 = this%knot3 + call this%set(knot1=knot1, knot2=knot_new, knot3=knot3, Xc=Xc_new) end if end do @@ -1864,7 +1911,9 @@ pure subroutine remove_knots(this, dir ,Xth,r) end do Wc_new(:) = Xcw_new(:,dim+1) - call this%set(knot1=this%knot1, knot2=this%knot2, knot3=knot_new, Xc=Xc_new, Wc=Wc_new) + knot1 = this%knot1 + knot2 = this%knot2 + call this%set(knot1=knot1, knot2=knot2, knot3=knot_new, Xc=Xc_new, Wc=Wc_new) deallocate(Xcw_new, Xc_new, Wc_new) end if end do @@ -1909,7 +1958,9 @@ pure subroutine remove_knots(this, dir ,Xth,r) Xc4 = reshape(Xc4, [this%nc(1),this%nc(2),nc_new,dim], order=[3,2,1,4]) Xc_new = reshape(Xc4, [this%nc(1)*this%nc(2)*(nc_new),dim]) - call this%set(knot1=this%knot1, knot2=this%knot2, knot3=knot_new, Xc=Xc_new) + knot1 = this%knot1 + knot2 = this%knot2 + call this%set(knot1=knot1, knot2=knot2, knot3=knot_new, Xc=Xc_new) end if end do