Add nearest_point() method.

- Add Xt to nurbs_surface and nurbs_volume derived types.
- Add examples for 1D, 2D, and 3D.
- Update fpm.toml.
- Update CHANGELOG.md.

CHANGELOG.md

@@ -18,6 +18,9 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
 - Added screenshots to the README file.
 - Included nvidia compiler in the CI.
 - Added missing allocation for `Tgc` in the `compute_Xg_nurbs_1d` subroutine.
+- Added `nearest_point()` procedures to the `nurbs_curve`, `nurbs_surface` and `nurbs_volume`.
+- Added examples for the `nearest_point()` method.
+- Added `Xt` to the `nurbs_surface` and `nurbs_volume` derived types.
 
 ### Changed
 

example/nearest_point_1d.f90

@@ -0,0 +1,44 @@
+program nearest_point_1d
+
+    use forcad, only: rk, nurbs_curve
+
+    implicit none
+
+    type(nurbs_curve) :: shape             !! Declare a NURBS curve object
+    real(rk), allocatable :: nearest_Xg(:) !! Coordinates of the nearest point on the curve
+    real(rk) :: nearest_Xt                 !! Corresponding parametric coordinates of the nearest point
+    integer :: id                          !! id of the nearest point
+
+    !-----------------------------------------------------------------------------
+    ! Setting up the NURBS circle
+    !-----------------------------------------------------------------------------
+
+    !> Set a circle with radius 2.0 and center at [0.0, 0.0, 0.0]
+    call shape%set_circle(center = [0.0_rk, 0.0_rk, 0.0_rk], radius = 2.0_rk)
+
+    !-----------------------------------------------------------------------------
+    ! Creating circle
+    !-----------------------------------------------------------------------------
+
+    !> Generate the NURBS circle with a resolution of 100
+    call shape%create(res = 100)
+
+    !-----------------------------------------------------------------------------
+    ! Nearest point on the curve
+    !-----------------------------------------------------------------------------
+
+    !> Find the nearest point on the curve to a given point
+    ! nearest_Xg: Coordinates of the nearest point on the curve (optional)
+    ! nearest_Xt: Corresponding parametric coordinates of the nearest point (optional)
+    ! id: id of the nearest point (optional)
+    call shape%nearest_point([2.0_rk, 3.0_rk, 5.0_rk], nearest_Xg, nearest_Xt, id)
+    print *, 'Nearest point on the curve:', nearest_Xg, 'with parametric coordinates:', nearest_Xt, 'and id:', id
+
+    !-----------------------------------------------------------------------------
+    ! Finalizing
+    !-----------------------------------------------------------------------------
+
+    !> Finalize the NURBS curve object
+    call shape%finalize()
+
+end program

example/nearest_point_2d.f90

@@ -1,12 +1,13 @@
 program nearest_point_2d
 
-    use forcad
+    use forcad, only: rk, nurbs_surface
 
     implicit none
 
-    type(nurbs_surface) :: shape             !! Declare a NURBS surface object
-    real(rk), allocatable :: nearest(:)   !! Array for nearest point on the surface
-    integer :: id                         !! id of the nearest point
+    type(nurbs_surface) :: shape           !! Declare a NURBS surface object
+    real(rk), allocatable :: nearest_Xg(:) !! Coordinates of the nearest point on the surface
+    real(rk), allocatable :: nearest_Xt(:) !! Corresponding parametric coordinates of the nearest point
+    integer :: id                          !! id of the nearest point
 
     !-----------------------------------------------------------------------------
     ! Setting up the NURBS tetrangon
@@ -28,8 +29,11 @@ program nearest_point_2d
     !-----------------------------------------------------------------------------
 
     !> Find the nearest point on the surface to a given point
-    call shape%nearest_point([2.0_rk, 3.0_rk, 5.0_rk], nearest, id)
-    print *, 'Nearest point on the surface:', nearest, 'with id:', id
+    ! nearest_Xg: Coordinates of the nearest point on the surface (optional)
+    ! nearest_Xt: Corresponding parametric coordinates of the nearest point (optional)
+    ! id: id of the nearest point (optional)
+    call shape%nearest_point([2.0_rk, 3.0_rk, 5.0_rk], nearest_Xg, nearest_Xt, id)
+    print *, 'Nearest point on the surface:', nearest_Xg, 'with parametric coordinates:', nearest_Xt, 'and id:', id
 
     !-----------------------------------------------------------------------------
     ! Finalizing

example/nearest_point_3d.f90

@@ -0,0 +1,45 @@
+program nearest_point_3d
+
+    use forcad, only: rk, nurbs_volume
+
+    implicit none
+
+    type(nurbs_volume) :: shape            !! Declare a NURBS volume object
+    real(rk), allocatable :: nearest_Xg(:) !! Coordinates of the nearest point on the volume
+    real(rk), allocatable :: nearest_Xt(:) !! Corresponding parametric coordinates of the nearest point
+    integer :: id                          !! id of the nearest point
+
+    !-----------------------------------------------------------------------------
+    ! Setting up the NURBS hexahedron
+    !-----------------------------------------------------------------------------
+
+    !> Set up a hexahedron shape with dimensions L = [2.0, 4.0, 8.0] and a specified number of control points nc = [4, 6, 8].
+    !> The weights of the control points (Wc) are optional.
+    call shape%set_hexahedron(L=[2.0_rk, 4.0_rk, 8.0_rk], nc=[4,6,8])
+
+    !-----------------------------------------------------------------------------
+    ! Creating the NURBS volume
+    !-----------------------------------------------------------------------------
+
+    !> Generate the NURBS volume with resolutions of 8, 16 and 32
+    call shape%create(8, 16, 32)
+
+    !-----------------------------------------------------------------------------
+    ! Nearest point on the volume
+    !-----------------------------------------------------------------------------
+
+    !> Find the nearest point on the volume to a given point
+    ! nearest_Xg: Coordinates of the nearest point on the volume (optional)
+    ! nearest_Xt: Corresponding parametric coordinates of the nearest point (optional)
+    ! id: id of the nearest point (optional)
+    call shape%nearest_point([2.0_rk, 3.0_rk, 5.0_rk], nearest_Xg, nearest_Xt, id)
+    print *, 'Nearest point on the volume:', nearest_Xg, 'with parametric coordinates:', nearest_Xt, 'and id:', id
+
+    !-----------------------------------------------------------------------------
+    ! Finalizing
+    !-----------------------------------------------------------------------------
+
+    !> Finalize the NURBS volume object
+    call shape%finalize()
+
+end program

fpm.toml

@@ -140,7 +140,17 @@ name       = "example_ppm3"
 source-dir = "example"
 main       = "example_ppm3.f90"
 
+[[example]]
+name       = "nearest_point_1d"
+source-dir = "example"
+main       = "nearest_point_1d.f90"
+
 [[example]]
 name       = "nearest_point_2d"
 source-dir = "example"
-main       = "nearest_point_2d.f90"
+main       = "nearest_point_2d.f90"
+
+[[example]]
+name       = "nearest_point_3d"
+source-dir = "example"
+main       = "nearest_point_3d.f90"

src/forcad_nurbs_curve.f90

@@ -80,6 +80,7 @@ module forcad_nurbs_curve
         procedure :: translate_Xc          !!> Translate control points
         procedure :: translate_Xg          !!> Translate geometry points
         procedure :: show                  !!> Show the NURBS object using PyVista
+        procedure :: nearest_point         !!> Find the nearest point on the NURBS curve
 
         ! Shapes
         procedure :: set_circle            !!> Set a circle
@@ -247,9 +248,11 @@ pure function compute_Xg_bspline_1d(f_Xt, f_knot, f_degree, f_nc, f_ng, f_Xc) re
         if (allocated(this%Xg)) deallocate(this%Xg)
 
         if (this%is_rational()) then ! NURBS
-            this%Xg = compute_Xg_nurbs_1d(this%Xt, this%knot, this%degree, this%nc, this%ng, this%Xc, this%Wc)
+            this%Xg = compute_Xg_nurbs_1d(&
+            this%Xt, this%knot, this%degree, this%nc, this%ng, this%Xc, this%Wc)
         else ! B-Spline
-            this%Xg = compute_Xg_bspline_1d(this%Xt, this%knot, this%degree, this%nc, this%ng, this%Xc)
+            this%Xg = compute_Xg_bspline_1d(&
+            this%Xt, this%knot, this%degree, this%nc, this%ng, this%Xc)
         end if
     end subroutine
     !===============================================================================
@@ -1483,6 +1486,32 @@ pure subroutine set_half_circle(this, center, radius)
     end subroutine
     !===============================================================================
 
+
+    !===============================================================================
+    !> author: Seyed Ali Ghasemi
+    !> license: BSD 3-Clause
+    pure subroutine nearest_point(this, point_Xg, nearest_Xg, nearest_Xt, id)
+        class(nurbs_curve), intent(in) :: this
+        real(rk), intent(in) :: point_Xg(:)
+        real(rk), intent(out), allocatable, optional :: nearest_Xg(:)
+        real(rk), intent(out), optional :: nearest_Xt
+        integer, intent(out), optional :: id
+        integer :: id_
+        real(rk), allocatable :: distances(:)
+        integer :: i
+
+        allocate(distances(this%ng))
+        do concurrent (i = 1: this%ng)
+            distances(i) = norm2(this%Xg(i,:) - point_Xg)
+        end do
+
+        id_ = minloc(distances, dim=1)
+        if (present(id)) id = id_
+        if (present(nearest_Xg)) nearest_Xg = this%Xg(id_,:)
+        if (present(nearest_Xt)) nearest_Xt = this%Xt(id_)
+    end subroutine
+    !===============================================================================
+
 end module forcad_nurbs_curve
 
 !===============================================================================

src/forcad_nurbs_surface.f90

@@ -21,6 +21,7 @@ module forcad_nurbs_surface
         real(rk), allocatable, private :: Wc(:)    !! Weights for control points (1D array: [nc(1)*nc(2)])
         real(rk), allocatable, private :: Xt1(:)   !! Evaluation parameter values in the first direction (1D array: [ng(1)])
         real(rk), allocatable, private :: Xt2(:)   !! Evaluation parameter values in the second direction (1D array: [ng(2)])
+        real(rk), allocatable, private :: Xt(:,:)  !! Evaluation parameter values (2D array: [ng(1)*ng(2), 2])
         real(rk), allocatable, private :: knot1(:) !! Knot vector in the first direction (1D array)
         real(rk), allocatable, private :: knot2(:) !! Knot vector in the second direction (1D array)
         integer, private :: degree(2)              !! Degree (order) of the surface
@@ -198,7 +199,6 @@ pure subroutine create(this, res1, res2, Xt1, Xt2, Xt)
         real(rk), intent(in), contiguous, optional :: Xt1(:), Xt2(:)
         real(rk), contiguous, intent(in), optional :: Xt(:,:)
         integer :: i
-        real(rk), allocatable :: Xt_(:,:)
 
         interface
             pure function compute_Xg_nurbs_2d(f_Xt, f_knot1, f_knot2, f_degree, f_nc, f_ng, f_Xc, f_Wc) result(f_Xg)
@@ -263,22 +263,24 @@ pure function compute_Xg_bspline_2d(f_Xt, f_knot1, f_knot2, f_degree, f_nc, f_ng
         end if
 
         if (present(Xt)) then
-            Xt_ = Xt
+            this%Xt = Xt
         else
 
             ! Set number of geometry points
             this%ng(1) = size(this%Xt1,1)
             this%ng(2) = size(this%Xt2,1)
 
-            call ndgrid(this%Xt1, this%Xt2, Xt_)
+            call ndgrid(this%Xt1, this%Xt2, this%Xt)
         end if
 
         if (allocated(this%Xg)) deallocate(this%Xg)
 
         if (this%is_rational()) then ! NURBS
-            this%Xg = compute_Xg_nurbs_2d(Xt_, this%knot1, this%knot2, this%degree, this%nc, this%ng, this%Xc, this%Wc)
+            this%Xg = compute_Xg_nurbs_2d(&
+            this%Xt, this%knot1, this%knot2, this%degree, this%nc, this%ng, this%Xc, this%Wc)
         else ! B-Spline
-            this%Xg = compute_Xg_bspline_2d(Xt_, this%knot1, this%knot2, this%degree, this%nc, this%ng, this%Xc)
+            this%Xg = compute_Xg_bspline_2d(&
+            this%Xt, this%knot1, this%knot2, this%degree, this%nc, this%ng, this%Xc)
         end if
     end subroutine
     !===============================================================================
@@ -616,6 +618,7 @@ pure subroutine finalize(this)
         if (allocated(this%Wc)) deallocate(this%Wc)
         if (allocated(this%Xt1)) deallocate(this%Xt1)
         if (allocated(this%Xt2)) deallocate(this%Xt2)
+        if (allocated(this%Xt)) deallocate(this%Xt)
         if (allocated(this%knot1)) deallocate(this%knot1)
         if (allocated(this%knot2)) deallocate(this%knot2)
         if (allocated(this%elemConn_Xc_vis)) deallocate(this%elemConn_Xc_vis)
@@ -2087,20 +2090,25 @@ pure subroutine set_half_ring(this, center, radius1, radius2)
     !===============================================================================
     !> author: Seyed Ali Ghasemi
     !> license: BSD 3-Clause
-    pure subroutine nearest_point(this, point, nearest, id)
+    pure subroutine nearest_point(this, point_Xg, nearest_Xg, nearest_Xt, id)
         class(nurbs_surface), intent(in) :: this
-        real(rk), intent(in) :: point(:)
-        real(rk), intent(out), allocatable :: nearest(:)
-        integer, intent(out) :: id
+        real(rk), intent(in) :: point_Xg(:)
+        real(rk), intent(out), allocatable, optional :: nearest_Xg(:)
+        real(rk), intent(out), allocatable, optional :: nearest_Xt(:)
+        integer, intent(out), optional :: id
+        integer :: id_
         real(rk), allocatable :: distances(:)
         integer :: i
 
         allocate(distances(this%ng(1)*this%ng(2)))
         do concurrent (i = 1: this%ng(1)*this%ng(2))
-            distances(i) = norm2(this%Xg(i,:) - point)
+            distances(i) = norm2(this%Xg(i,:) - point_Xg)
         end do
-        id = minloc(distances, dim=1)
-        nearest = this%Xg(id,:)
+
+        id_ = minloc(distances, dim=1)
+        if (present(id)) id = id_
+        if (present(nearest_Xg)) nearest_Xg = this%Xg(id_,:)
+        if (present(nearest_Xt)) nearest_Xt = this%Xt(id_,:)
     end subroutine
     !===============================================================================
 

src/forcad_nurbs_volume.f90

@@ -22,6 +22,7 @@ module forcad_nurbs_volume
         real(rk), allocatable, private :: Xt1(:)   !! Evaluation parameter values in the first direction (1D array: [ng(1)])
         real(rk), allocatable, private :: Xt2(:)   !! Evaluation parameter values in the second direction (1D array: [ng(2)])
         real(rk), allocatable, private :: Xt3(:)   !! Evaluation parameter values in the third direction (1D array: [ng(3)])
+        real(rk), allocatable, private :: Xt(:,:)  !! Evaluation parameter values (2D array: [ng(1)*ng(2)*ng(3), dim]
         real(rk), allocatable, private :: knot1(:) !! Knot vector in the first direction (1D array)
         real(rk), allocatable, private :: knot2(:) !! Knot vector in the second direction (1D array)
         real(rk), allocatable, private :: knot3(:) !! Knot vector in the third direction (1D array)
@@ -86,6 +87,7 @@ module forcad_nurbs_volume
         procedure :: translate_Xc           !!> Translate control points
         procedure :: translate_Xg           !!> Translate geometry points
         procedure :: show                   !!> Show the NURBS object using PyVista
+        procedure :: nearest_point          !!> Find the nearest point on the NURBS volume
 
         ! Shapes
         procedure :: set_hexahedron         !!> Set a hexahedron
@@ -207,7 +209,6 @@ pure subroutine create(this, res1, res2, res3, Xt1, Xt2, Xt3, Xt)
         real(rk), intent(in), contiguous, optional :: Xt1(:), Xt2(:), Xt3(:)
         real(rk), intent(in), contiguous, optional :: Xt(:,:)
         integer :: i
-        real(rk), allocatable :: Xt_(:,:)
 
         interface
             pure function compute_Xg_nurbs_3d(f_Xt, f_knot1, f_knot2, f_knot3, f_degree, f_nc, f_ng, f_Xc, f_Wc) result(f_Xg)
@@ -282,24 +283,26 @@ pure function compute_Xg_bspline_3d(f_Xt, f_knot1, f_knot2, f_knot3, f_degree, f
         end if
 
         if (present(Xt)) then
-            Xt_ = Xt
+            this%Xt = Xt
         else
 
             ! Set number of geometry points
             this%ng(1) = size(this%Xt1,1)
             this%ng(2) = size(this%Xt2,1)
             this%ng(3) = size(this%Xt3,1)
 
-            call ndgrid(this%Xt1, this%Xt2, this%Xt3, Xt_)
+            call ndgrid(this%Xt1, this%Xt2, this%Xt3, this%Xt)
         end if
 
         if (allocated(this%Xg)) deallocate(this%Xg)
         allocate(this%Xg(this%ng(1)*this%ng(2)*this%ng(3), size(this%Xc,2)))
 
         if (allocated(this%Wc)) then ! NURBS volume
-            this%Xg = compute_Xg_nurbs_3d(Xt_, this%knot1, this%knot2, this%knot3, this%degree, this%nc, this%ng, this%Xc, this%Wc)
+            this%Xg = compute_Xg_nurbs_3d(&
+            this%Xt, this%knot1, this%knot2, this%knot3, this%degree, this%nc, this%ng, this%Xc, this%Wc)
         else
-            this%Xg = compute_Xg_bspline_3d(Xt_, this%knot1, this%knot2, this%knot3, this%degree, this%nc, this%ng, this%Xc)
+            this%Xg = compute_Xg_bspline_3d(&
+            this%Xt, this%knot1, this%knot2, this%knot3, this%degree, this%nc, this%ng, this%Xc)
         end if
     end subroutine
     !===============================================================================
@@ -670,6 +673,7 @@ pure subroutine finalize(this)
         if (allocated(this%Xt1)) deallocate(this%Xt1)
         if (allocated(this%Xt2)) deallocate(this%Xt2)
         if (allocated(this%Xt3)) deallocate(this%Xt3)
+        if (allocated(this%Xt)) deallocate(this%Xt)
         if (allocated(this%knot1)) deallocate(this%knot1)
         if (allocated(this%knot2)) deallocate(this%knot2)
         if (allocated(this%knot3)) deallocate(this%knot3)
@@ -2573,6 +2577,32 @@ pure subroutine set_half_ring(this, center, radius1, radius2, length)
     end subroutine
     !===============================================================================
 
+
+    !===============================================================================
+    !> author: Seyed Ali Ghasemi
+    !> license: BSD 3-Clause
+    pure subroutine nearest_point(this, point_Xg, nearest_Xg, nearest_Xt, id)
+        class(nurbs_volume), intent(in) :: this
+        real(rk), intent(in) :: point_Xg(:)
+        real(rk), intent(out), allocatable, optional :: nearest_Xg(:)
+        real(rk), intent(out), allocatable, optional :: nearest_Xt(:)
+        integer, intent(out), optional :: id
+        integer :: id_
+        real(rk), allocatable :: distances(:)
+        integer :: i
+
+        allocate(distances(this%ng(1)*this%ng(2)*this%ng(3)))
+        do concurrent (i = 1: this%ng(1)*this%ng(2)*this%ng(3))
+            distances(i) = norm2(this%Xg(i,:) - point_Xg)
+        end do
+
+        id_ = minloc(distances, dim=1)
+        if (present(id)) id = id_
+        if (present(nearest_Xg)) nearest_Xg = this%Xg(id_,:)
+        if (present(nearest_Xt)) nearest_Xt = this%Xt(id_,:)
+    end subroutine
+    !===============================================================================
+
 end module forcad_nurbs_volume
 
 !===============================================================================