Added derivative2, cmp_Xg, nearest_point2

Signed-off-by: Seyed Ali Ghasemi <info@gha3mi.com>
gha3mi · Jun 22, 2024 · dd7f080 · dd7f080
1 parent cf4908e
commit dd7f080
Show file tree

Hide file tree

Showing 12 changed files with 3,025 additions and 414 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -28,6 +28,12 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
 - Added `cmp_elemFace()` to extract element connectivity of faces.
 - Added `cmp_degreeFace()` to extract the degrees of faces.
 - Updated `example_volume_1` to use `cmp_elemFace()` and `cmp_degreeFace()`.
+- Added `cmp_Xg()` to evaluate the geometry points.
+- Added generic method `derivative2()` to compute the second derivative of a NURBS objects.
+- Added `nearest_point2()` to compute the nearest point on a NURBS object using optimization.
+- Added Interfaces.
+- Added new tests: fdm_curve.f90, fdm_surface.f90, fdm_volume.f90.
+- Updated nearest_point_* examples to use the new `nearest_point2()` method.
 
 ### Changed
 
@@ -46,6 +52,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
 - Used `do concurrent` within `kron` and `basis_bspline` functions.
 - Converted the `basis_bspline_der` function to a subroutine.
 - Fixed NURBS derivative calculations.
+- Made `basis` and `derivative` generic methods.
 
 ### Removed
 

diff --git a/example/nearest_point_1d.f90 b/example/nearest_point_1d.f90
@@ -4,24 +4,40 @@ program nearest_point_1d
 
     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
+    type(nurbs_curve) :: shape              !! Declare a NURBS curve object
+    real(rk), allocatable :: Xc(:,:), Wc(:) !! Arrays for control points and weights
+    real(rk) :: knot(6)                     !! Array for knot vector
+    real(rk), allocatable :: nearest_Xg(:)  !! Array for the nearest point on the curve
+    real(rk) :: nearest_Xt                  !! Array for the parametric coordinates of the nearest point
+    integer :: id                           !! Variable for the id of the nearest point
 
     !-----------------------------------------------------------------------------
-    ! Setting up the NURBS circle
+    ! Setting up the NURBS curve
     !-----------------------------------------------------------------------------
 
-    !> 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)
+    !> Define control points for the NURBS curve
+    allocate(Xc(3, 3))
+    Xc(1,:) = [0.0_rk, 0.0_rk, 0.0_rk]
+    Xc(2,:) = [0.0_rk, 5.0_rk, 0.0_rk]
+    Xc(3,:) = [5.0_rk, 5.0_rk, 0.0_rk]
+
+    !> Define weights for the control points (optional)
+    allocate(Wc(3))
+    Wc = [1.0_rk, 1.1_rk, 1.0_rk]
+
+    !> Define knot vector
+    knot = [0.0_rk, 0.0_rk, 0.0_rk, 1.0_rk, 1.0_rk, 1.0_rk]
+
+    !> Set knot vector, control points, and weights for the NURBS curve object.
+    !> Wc is optional
+    call shape%set(knot, Xc, Wc)
 
     !-----------------------------------------------------------------------------
-    ! Creating circle
+    ! Creating the NURBS curve
     !-----------------------------------------------------------------------------
 
-    !> Generate the NURBS circle with a resolution of 100
-    call shape%create(res = 100)
+    !> Generate the NURBS curve with a resolution of 20
+    call shape%create(20)
 
     !-----------------------------------------------------------------------------
     ! Nearest point on the curve
@@ -31,15 +47,27 @@ program nearest_point_1d
     ! 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
+    call shape%nearest_point([4.5_rk, 4.5_rk, 5.0_rk], nearest_Xg, nearest_Xt, id)
+    print '(a,1x,g0,2x,g0,2x,g0,a,2x,g0,2x,a,1x,g0)','Nearest point on the curve:', nearest_Xg, ' with parametric coordinates:', nearest_Xt, ' and id:', id
+
+    !-----------------------------------------------------------------------------
+    ! Nearest point on the curve (Optimization)
+    !-----------------------------------------------------------------------------
+
+    !> Find the nearest point on the curve to a given point
+    !> The optimization method is used to find the nearest point
+    !> The optimization method is based on the Newton-Raphson method
+    ! nearest_Xt: Corresponding parametric coordinates of the nearest point
+    ! nearest_Xg: Coordinates of the nearest point on the curve (optional)
+    call shape%nearest_point2([4.5_rk, 4.5_rk, 5.0_rk], 1.0e-11_rk, 30, nearest_Xt, nearest_Xg)
+    print '(a,1x,g0,2x,g0,a,2x,g0,2x,g0)', 'Nearest point on the curve:', nearest_Xg, ' with parametric coordinates:', nearest_Xt
 
     !-----------------------------------------------------------------------------
     ! Finalizing
     !-----------------------------------------------------------------------------
 
     !> Finalize the NURBS curve object
     call shape%finalize()
-    deallocate(nearest_Xg)
+    deallocate(nearest_Xg, Xc, Wc)
 
 end program
diff --git a/example/nearest_point_2d.f90 b/example/nearest_point_2d.f90
@@ -8,14 +8,23 @@ program nearest_point_2d
     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
+    real(rk) :: Xc(4,3)                    !! Control points
+    real(rk) :: Wc(4)                      !! Weights of the control points
 
     !-----------------------------------------------------------------------------
     ! Setting up the NURBS tetrangon
     !-----------------------------------------------------------------------------
 
-    !> Set a tetragon with lengths of 2.0 and 3.0 and 3 and 4 control points in each direction
+    !> Set a surface with 4 control points
+    Xc(1,:) = [0.0_rk, 0.0_rk, 0.0_rk]
+    Xc(2,:) = [2.0_rk, 0.0_rk, 0.0_rk]
+    Xc(3,:) = [0.0_rk, 2.0_rk, 0.0_rk]
+    Xc(4,:) = [2.0_rk, 2.0_rk, 0.0_rk]
+
     !> The weights of the control points (Wc) are optional.
-    call shape%set_tetragon(L=[2.0_rk, 3.0_rk], nc=[3,4])
+    Wc = [1.0_rk, 1.1_rk, 0.7_rk, 1.0_rk]
+
+    call shape%set(knot1=[0.0_rk, 0.0_rk, 1.0_rk, 1.0_rk], knot2=[0.0_rk, 0.0_rk, 1.0_rk, 1.0_rk], Xc=Xc, Wc=Wc)
 
     !-----------------------------------------------------------------------------
     ! Creating the NURBS surface
@@ -25,22 +34,34 @@ program nearest_point_2d
     call shape%create(30, 30)
 
     !-----------------------------------------------------------------------------
-    ! Nearest point on the surface
+    ! Nearest point on the surface (Approximation)
     !-----------------------------------------------------------------------------
 
     !> Find the nearest point on the surface to a given point
     ! 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
+    call shape%nearest_point([1.3_rk, 1.0_rk, 1.999999999_rk], nearest_Xg, nearest_Xt, id)
+    print '(a,1x,g0,2x,g0,2x,g0,a,2x,g0,2x,g0,2x,a,1x,g0)','Nearest point on the surface:', nearest_Xg, ' with parametric coordinates:', nearest_Xt, ' and id:', id
+
+    !-----------------------------------------------------------------------------
+    ! Nearest point on the surface (Optimization)
+    !-----------------------------------------------------------------------------
+
+    !> Find the nearest point on the surface to a given point
+    !> The optimization method is used to find the nearest point
+    !> The optimization method is based on the Newton-Raphson method
+    ! nearest_Xt: Corresponding parametric coordinates of the nearest point
+    ! nearest_Xg: Coordinates of the nearest point on the surface (optional)
+    call shape%nearest_point2([1.3_rk, 1.0_rk, 1.999999999_rk], 1.0e-11_rk, 30, nearest_Xt, nearest_Xg)
+    print '(a,1x,g0,2x,g0,2x,g0,a,2x,g0,2x,g0)', 'Nearest point on the surface:', nearest_Xg, ' with parametric coordinates:', nearest_Xt
 
     !-----------------------------------------------------------------------------
     ! Finalizing
     !-----------------------------------------------------------------------------
 
     !> Finalize the NURBS surface object
     call shape%finalize()
-    deallocate(nearest_Xg, nearest_Xt)
+    ! deallocate(nearest_Xg, nearest_Xt)
 
 end program
diff --git a/example/nearest_point_3d.f90 b/example/nearest_point_3d.f90
@@ -8,21 +8,33 @@ program nearest_point_3d
     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
+    real(rk) :: Xc(8,3)                    !! Control points
+    real(rk) :: Wc(8)                      !! Weights of the control points
 
     !-----------------------------------------------------------------------------
     ! 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].
+    Xc(1,:) = [0.0_rk, 0.0_rk, 0.0_rk]
+    Xc(2,:) = [2.0_rk, 0.0_rk, 0.0_rk]
+    Xc(3,:) = [0.0_rk, 4.0_rk, 0.0_rk]
+    Xc(4,:) = [2.0_rk, 4.0_rk, 0.0_rk]
+    Xc(5,:) = [0.0_rk, 0.0_rk, 2.0_rk]
+    Xc(6,:) = [2.0_rk, 0.0_rk, 2.0_rk]
+    Xc(7,:) = [0.0_rk, 4.0_rk, 2.0_rk]
+    Xc(8,:) = [2.0_rk, 4.0_rk, 2.0_rk]
+
     !> 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])
+    Wc = [1.0_rk, 1.1_rk, 1.11_rk, 1.0_rk, 0.5_rk, 0.5_rk, 1.2_rk, 1.0_rk]
+
+    call shape%set(knot1=[0.0_rk, 0.0_rk, 1.0_rk, 1.0_rk], knot2=[0.0_rk, 0.0_rk, 1.0_rk, 1.0_rk], knot3=[0.0_rk, 0.0_rk, 1.0_rk, 1.0_rk], Xc=Xc, Wc=Wc)
 
     !-----------------------------------------------------------------------------
     ! Creating the NURBS volume
     !-----------------------------------------------------------------------------
 
-    !> Generate the NURBS volume with resolutions of 8, 16 and 32
-    call shape%create(8, 16, 32)
+    !> Generate the NURBS volume with resolutions of 20, 20, 20
+    call shape%create(30, 30, 30)
 
     !-----------------------------------------------------------------------------
     ! Nearest point on the volume
@@ -32,15 +44,27 @@ program nearest_point_3d
     ! 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
+    call shape%nearest_point([1.5_rk, 3.5_rk, 1.1_rk], nearest_Xg, nearest_Xt, id)
+    print '(a,1x,g0,2x,g0,2x,g0,a,2x,g0,2x,g0,2x,g0,2x,a,1x,g0)','Nearest point on the volume:', nearest_Xg, ' with parametric coordinates:', nearest_Xt, ' and id:', id
+
+    !-----------------------------------------------------------------------------
+    ! Nearest point on the volume (Optimization)
+    !-----------------------------------------------------------------------------
+
+    !> Find the nearest point on the volume to a given point
+    !> The optimization method is used to find the nearest point
+    !> The optimization method is based on the Newton-Raphson method
+    ! nearest_Xt: Corresponding parametric coordinates of the nearest point
+    ! nearest_Xg: Coordinates of the nearest point on the volume (optional)
+    call shape%nearest_point2([1.5_rk, 3.5_rk, 1.1_rk], 1.0e-11_rk, 500, nearest_Xt, nearest_Xg)
+    print '(a,1x,g0,2x,g0,2x,g0,a,2x,g0,2x,g0,2x,g0)', 'Nearest point on the volume:', nearest_Xg, ' with parametric coordinates:', nearest_Xt
 
     !-----------------------------------------------------------------------------
     ! Finalizing
     !-----------------------------------------------------------------------------
 
     !> Finalize the NURBS volume object
     call shape%finalize()
-    deallocate(nearest_Xg, nearest_Xt)
+    ! deallocate(nearest_Xg, nearest_Xt)
 
 end program