diff --git a/.github/workflows/fpm.yml b/.github/workflows/fpm.yml index abcb974a4..a72b02f65 100644 --- a/.github/workflows/fpm.yml +++ b/.github/workflows/fpm.yml @@ -13,15 +13,12 @@ jobs: toolchain: - {compiler: gcc} - {compiler: intel} - - {compiler: intel-classic} - {compiler: nvidia-hpc} exclude: - os: windows-latest toolchain: {compiler: nvidia-hpc} - os: windows-latest toolchain: {compiler: intel} - - os: windows-latest - toolchain: {compiler: intel-classic} - os: macos-latest toolchain: {compiler: nvidia-hpc} diff --git a/CHANGELOG.md b/CHANGELOG.md index 3915ea6fc..85b35f637 100644 --- 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 @@ -44,7 +50,11 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), - Improved `elemConn_C0` to support higher order elements and check for the number of elements. - Renamed test files to match the module name. - 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 -- Excluded macOS from CI due to a problem with fpm. \ No newline at end of file +- Excluded macOS from CI due to a problem with fpm. +- Removed `intel-classic` from CI (`ifort` is deprecated). \ No newline at end of file diff --git a/example/nearest_point_1d.f90 b/example/nearest_point_1d.f90 index 4bfcfd5c3..962bce256 100644 --- 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,8 +47,22 @@ 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 @@ -40,6 +70,6 @@ program nearest_point_1d !> 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 index 1aaa6fd8d..5de270744 100644 --- 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,15 +34,29 @@ 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 @@ -41,6 +64,6 @@ program nearest_point_2d !> 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 index b2252ab37..41682e922 100644 --- a/example/nearest_point_3d.f90 +++ b/example/nearest_point_3d.f90 @@ -8,21 +8,37 @@ 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,8 +48,22 @@ 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 @@ -41,6 +71,6 @@ program nearest_point_3d !> Finalize the NURBS volume object call shape%finalize() - deallocate(nearest_Xg, nearest_Xt) + ! deallocate(nearest_Xg, nearest_Xt) end program diff --git a/src/forcad_nurbs_curve.f90 b/src/forcad_nurbs_curve.f90 index c32154b87..f59c4c8f0 100644 --- a/src/forcad_nurbs_curve.f90 +++ b/src/forcad_nurbs_curve.f90 @@ -33,6 +33,7 @@ module forcad_nurbs_curve procedure :: set3 !!> Set Bezier or Rational Bezier curve using control points and weights generic :: set => set1, set2, set3 !!> Set NURBS curve procedure :: create !!> Generate geometry points + procedure :: cmp_Xg !!> Compute geometry points procedure, private :: get_Xc_all !!> Get all control points procedure, private :: get_Xci !!> Get i-th control point procedure, private :: get_Xcid !!> Get i-th control point in a specific direction @@ -71,8 +72,15 @@ module forcad_nurbs_curve procedure :: get_nc !!> Get number of control points procedure :: insert_knots !!> Insert knots into the knot vector 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, private :: basis_vector !!> Compute the basis functions of the NURBS curve + procedure, private :: basis_scalar !!> Compute the basis functions of the NURBS curve + generic :: basis => basis_vector, basis_scalar !!> Compute the basis functions of the NURBS curve + procedure, private :: derivative_vector !!> Compute the derivative of the NURBS curve + procedure, private :: derivative_scalar !!> Compute the derivative of the NURBS curve + generic :: derivative => derivative_vector, derivative_scalar !!> Compute the derivative of the NURBS curve + procedure, private :: derivative2_vector !!> Compute the second derivative of the NURBS curve + procedure, private :: derivative2_scalar !!> Compute the second derivative of the NURBS curve + generic :: derivative2 => derivative2_vector, derivative2_scalar !!> Compute the second derivative of the NURBS curve procedure :: is_rational !!> Check if the NURBS curve is rational procedure :: remove_knots !!> Remove knots from the knot vector procedure :: rotate_Xc !!> Rotate control points @@ -80,7 +88,8 @@ 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 + procedure :: nearest_point !!> Find the nearest point on the NURBS curve (Approximation) + procedure :: nearest_point2 !!> Find the nearest point on the NURBS curve (Minimization - Newton's method) ! Shapes procedure :: set_circle !!> Set a circle @@ -89,6 +98,201 @@ module forcad_nurbs_curve end type !=============================================================================== + + interface compute_Xg + pure function compute_Xg_nurbs_1d(f_Xt, f_knot, f_degree, f_nc, f_ng, f_Xc, f_Wc) result(f_Xg) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:) + real(rk), intent(in), contiguous :: f_knot(:) + integer, intent(in) :: f_degree + integer, intent(in) :: f_nc + integer, intent(in) :: f_ng + real(rk), intent(in), contiguous :: f_Xc(:,:) + real(rk), intent(in), contiguous :: f_Wc(:) + real(rk), allocatable :: f_Xg(:,:) + end function + + pure function compute_Xg_bspline_1d(f_Xt, f_knot, f_degree, f_nc, f_ng, f_Xc) result(f_Xg) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:) + real(rk), intent(in), contiguous :: f_knot(:) + integer, intent(in) :: f_degree + integer, intent(in) :: f_nc + integer, intent(in) :: f_ng + real(rk), intent(in), contiguous :: f_Xc(:,:) + real(rk), allocatable :: f_Xg(:,:) + end function + + pure function compute_Xg_nurbs_1d_1point(f_Xt, f_knot, f_degree, f_nc, f_Xc, f_Wc) result(f_Xg) + import :: rk + real(rk), intent(in) :: f_Xt + real(rk), intent(in), contiguous :: f_knot(:) + integer, intent(in) :: f_degree + integer, intent(in) :: f_nc + real(rk), intent(in), contiguous :: f_Xc(:,:) + real(rk), intent(in), contiguous :: f_Wc(:) + real(rk), allocatable :: f_Xg(:) + end function + + pure function compute_Xg_bspline_1d_1point(f_Xt, f_knot, f_degree, f_nc, f_Xc) result(f_Xg) + import :: rk + real(rk), intent(in) :: f_Xt + real(rk), intent(in), contiguous :: f_knot(:) + integer, intent(in) :: f_degree + integer, intent(in) :: f_nc + real(rk), intent(in), contiguous :: f_Xc(:,:) + real(rk), allocatable :: f_Xg(:) + end function + end interface + + interface compute_d2Tgc + pure subroutine compute_d2Tgc_nurbs_1d_vector(f_Xt, f_knot, f_degree, f_nc, f_ng, f_Wc, f_d2Tgc, f_dTgc, f_Tgc) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:) + real(rk), intent(in), contiguous :: f_knot(:) + integer, intent(in) :: f_degree + integer, intent(in) :: f_nc + integer, intent(in) :: f_ng + real(rk), intent(in), contiguous :: f_Wc(:) + real(rk), allocatable, intent(out) :: f_d2Tgc(:,:) + real(rk), allocatable, intent(out) :: f_dTgc(:,:) + real(rk), allocatable, intent(out) :: f_Tgc(:,:) + end subroutine + + pure subroutine compute_d2Tgc_bspline_1d_vector(f_Xt, f_knot, f_degree, f_nc, f_ng, f_d2Tgc, f_dTgc, f_Tgc) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:) + real(rk), intent(in), contiguous :: f_knot(:) + integer, intent(in) :: f_degree + integer, intent(in) :: f_nc + integer, intent(in) :: f_ng + real(rk), allocatable, intent(out) :: f_d2Tgc(:,:) + real(rk), allocatable, intent(out) :: f_dTgc(:,:) + real(rk), allocatable, intent(out) :: f_Tgc(:,:) + end subroutine + + pure subroutine compute_d2Tgc_nurbs_1d_scalar(f_Xt, f_knot, f_degree, f_nc, f_Wc, f_d2Tgc, f_dTgc, f_Tgc) + import :: rk + real(rk), intent(in) :: f_Xt + real(rk), intent(in), contiguous :: f_knot(:) + integer, intent(in) :: f_degree + integer, intent(in) :: f_nc + real(rk), intent(in), contiguous :: f_Wc(:) + real(rk), allocatable, intent(out) :: f_d2Tgc(:) + real(rk), allocatable, intent(out) :: f_dTgc(:) + real(rk), allocatable, intent(out) :: f_Tgc(:) + end subroutine + + pure subroutine compute_d2Tgc_bspline_1d_scalar(f_Xt, f_knot, f_degree, f_nc, f_d2Tgc, f_dTgc, f_Tgc) + import :: rk + real(rk), intent(in) :: f_Xt + real(rk), intent(in), contiguous :: f_knot(:) + integer, intent(in) :: f_degree + integer, intent(in) :: f_nc + real(rk), allocatable, intent(out) :: f_d2Tgc(:) + real(rk), allocatable, intent(out) :: f_dTgc(:) + real(rk), allocatable, intent(out) :: f_Tgc(:) + end subroutine + end interface + + interface compute_dTgc + pure subroutine compute_dTgc_nurbs_1d_vector(f_Xt, f_knot, f_degree, f_nc, f_ng, f_Wc, f_dTgc, f_Tgc) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:) + real(rk), intent(in), contiguous :: f_knot(:) + integer, intent(in) :: f_degree + integer, intent(in) :: f_nc + integer, intent(in) :: f_ng + real(rk), intent(in), contiguous :: f_Wc(:) + real(rk), allocatable, intent(out) :: f_dTgc(:,:) + real(rk), allocatable, intent(out) :: f_Tgc(:,:) + end subroutine + + pure subroutine compute_dTgc_bspline_1d_vector(f_Xt, f_knot, f_degree, f_nc, f_ng, f_dTgc, f_Tgc) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:) + real(rk), intent(in), contiguous :: f_knot(:) + integer, intent(in) :: f_degree + integer, intent(in) :: f_nc + integer, intent(in) :: f_ng + real(rk), allocatable, intent(out) :: f_dTgc(:,:) + real(rk), allocatable, intent(out) :: f_Tgc(:,:) + end subroutine + + pure subroutine compute_dTgc_nurbs_1d_scalar(f_Xt, f_knot, f_degree, f_nc, f_Wc, f_dTgc, f_Tgc) + import :: rk + real(rk), intent(in) :: f_Xt + real(rk), intent(in), contiguous :: f_knot(:) + integer, intent(in) :: f_degree + integer, intent(in) :: f_nc + real(rk), intent(in), contiguous :: f_Wc(:) + real(rk), allocatable, intent(out) :: f_dTgc(:) + real(rk), allocatable, intent(out) :: f_Tgc(:) + end subroutine + + pure subroutine compute_dTgc_bspline_1d_scalar(f_Xt, f_knot, f_degree, f_nc, f_dTgc, f_Tgc) + import :: rk + real(rk), intent(in) :: f_Xt + real(rk), intent(in), contiguous :: f_knot(:) + integer, intent(in) :: f_degree + integer, intent(in) :: f_nc + real(rk), allocatable, intent(out) :: f_dTgc(:) + real(rk), allocatable, intent(out) :: f_Tgc(:) + end subroutine + end interface + + interface compute_Tgc + pure function compute_Tgc_nurbs_1d_vector(f_Xt, f_knot, f_degree, f_nc, f_ng, f_Wc) result(f_Tgc) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:) + real(rk), intent(in), contiguous :: f_knot(:) + integer, intent(in) :: f_degree + integer, intent(in) :: f_nc + integer, intent(in) :: f_ng + real(rk), intent(in), contiguous :: f_Wc(:) + real(rk), allocatable :: f_Tgc(:,:) + end function + + pure function compute_Tgc_bspline_1d_vector(f_Xt, f_knot, f_degree, f_nc, f_ng) result(f_Tgc) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:) + real(rk), intent(in), contiguous :: f_knot(:) + integer, intent(in) :: f_degree + integer, intent(in) :: f_nc + integer, intent(in) :: f_ng + real(rk), allocatable :: f_Tgc(:,:) + end function + + pure function compute_Tgc_nurbs_1d_scalar(f_Xt, f_knot, f_degree, f_nc, f_Wc) result(f_Tgc) + import :: rk + real(rk), intent(in) :: f_Xt + real(rk), intent(in), contiguous :: f_knot(:) + integer, intent(in) :: f_degree + integer, intent(in) :: f_nc + real(rk), intent(in), contiguous :: f_Wc(:) + real(rk), allocatable :: f_Tgc(:) + end function + + pure function compute_Tgc_bspline_1d_scalar(f_Xt, f_knot, f_degree, f_nc) result(f_Tgc) + import :: rk + real(rk), intent(in) :: f_Xt + real(rk), intent(in), contiguous :: f_knot(:) + integer, intent(in) :: f_degree + integer, intent(in) :: f_nc + real(rk), allocatable :: f_Tgc(:) + end function + end interface + + interface + pure function nearest_point_help_1d(f_ng, f_Xg, f_point_Xg) result(f_distances) + import :: rk + integer, intent(in) :: f_ng + real(rk), intent(in), contiguous :: f_Xg(:,:) + real(rk), intent(in), contiguous :: f_point_Xg(:) + real(rk), allocatable :: f_distances(:) + end function + end interface + contains !=============================================================================== @@ -193,33 +397,6 @@ pure subroutine create(this, res, Xt) real(rk), intent(in), contiguous, optional :: Xt(:) integer :: i - interface - pure function compute_Xg_nurbs_1d(f_Xt, f_knot, f_degree, f_nc, f_ng, f_Xc, f_Wc) result(f_Xg) - import :: rk - real(rk), intent(in), contiguous :: f_Xt(:) - real(rk), intent(in), contiguous :: f_knot(:) - integer, intent(in) :: f_degree - integer, intent(in) :: f_nc - integer, intent(in) :: f_ng - real(rk), intent(in), contiguous :: f_Xc(:,:) - real(rk), intent(in), contiguous :: f_Wc(:) - real(rk), allocatable :: f_Xg(:,:) - end function - end interface - - interface - pure function compute_Xg_bspline_1d(f_Xt, f_knot, f_degree, f_nc, f_ng, f_Xc) result(f_Xg) - import :: rk - real(rk), intent(in), contiguous :: f_Xt(:) - real(rk), intent(in), contiguous :: f_knot(:) - integer, intent(in) :: f_degree - integer, intent(in) :: f_nc - integer, intent(in) :: f_ng - real(rk), intent(in), contiguous :: f_Xc(:,:) - real(rk), allocatable :: f_Xg(:,:) - end function - end interface - ! check if (.not.allocated(this%Xc)) then error stop 'Control points are not set.' @@ -248,16 +425,42 @@ 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(& + 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(& + this%Xt, this%knot, this%degree, this%nc, this%ng, this%Xc) end if end subroutine !=============================================================================== + !=============================================================================== + !> author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + pure function cmp_Xg(this, Xt) result(Xg) + class(nurbs_curve), intent(in) :: this + real(rk), intent(in) :: Xt + real(rk), allocatable :: Xg(:) + + ! check + if (.not.allocated(this%Xc)) then + error stop 'Control points are not set.' + end if + + if (.not.allocated(this%knot)) then + error stop 'Knot vector is not set.' + end if + + if (this%is_rational()) then ! NURBS + Xg = compute_Xg(Xt, this%knot, this%degree, this%nc, this%Xc, this%Wc) + else ! B-Spline + Xg = compute_Xg(Xt, this%knot, this%degree, this%nc, this%Xc) + end if + end function + !=============================================================================== + + !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause @@ -883,37 +1086,67 @@ pure subroutine elevate_degree(this, t) !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause - pure subroutine derivative(this, res, Xt, dTgc) + pure subroutine derivative_vector(this, res, Xt, dTgc, Tgc) class(nurbs_curve), intent(inout) :: this integer, intent(in), optional :: res real(rk), intent(in), contiguous, optional :: Xt(:) real(rk), allocatable, intent(out) :: dTgc(:,:) + real(rk), allocatable, intent(out), optional :: Tgc(:,:) integer :: i - interface - pure function compute_dTgc_nurbs_1d(f_Xt, f_knot, f_degree, f_nc, f_ng, f_Wc) result(f_dTgc) - import :: rk - real(rk), intent(in), contiguous :: f_Xt(:) - real(rk), intent(in), contiguous :: f_knot(:) - integer, intent(in) :: f_degree - integer, intent(in) :: f_nc - integer, intent(in) :: f_ng - real(rk), intent(in), contiguous :: f_Wc(:) - real(rk), allocatable :: f_dTgc(:,:) - end function - end interface - - interface - pure function compute_dTgc_bspline_1d(f_Xt, f_knot, f_degree, f_nc, f_ng) result(f_dTgc) - import :: rk - real(rk), intent(in), contiguous :: f_Xt(:) - real(rk), intent(in), contiguous :: f_knot(:) - integer, intent(in) :: f_degree - integer, intent(in) :: f_nc - integer, intent(in) :: f_ng - real(rk), allocatable :: f_dTgc(:,:) - end function - end interface + ! Set parameter values + if (present(Xt)) then + if (allocated(this%Xt)) deallocate(this%Xt) + this%Xt = Xt + elseif (present(res)) then + if (allocated(this%Xt)) deallocate(this%Xt) + allocate(this%Xt(res)) + this%Xt = [(real(i-1, rk) / real(res-1, rk), i=1, res)] + ! else + ! this%Xt = this%Xt + end if + + ! Set number of geometry points + this%ng = size(this%Xt,1) + + if (this%is_rational()) then ! NURBS + call compute_dTgc(this%Xt, this%knot, this%degree, this%nc, this%ng, this%Wc, dTgc, Tgc) + else ! B-Spline + call compute_dTgc(this%Xt, this%knot, this%degree, this%nc, this%ng, dTgc, Tgc) + end if + end subroutine + !=============================================================================== + + + !=============================================================================== + !> author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + pure subroutine derivative_scalar(this, Xt, dTgc, Tgc) + class(nurbs_curve), intent(inout) :: this + real(rk), intent(in) :: Xt + real(rk), allocatable, intent(out) :: dTgc(:) + real(rk), allocatable, intent(out), optional :: Tgc(:) + + if (this%is_rational()) then ! NURBS + call compute_dTgc(Xt, this%knot, this%degree, this%nc, this%Wc, dTgc, Tgc) + else ! B-Spline + call compute_dTgc(Xt, this%knot, this%degree, this%nc, dTgc, Tgc) + end if + end subroutine + !=============================================================================== + + + !=============================================================================== + !> author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + pure subroutine derivative2_vector(this, res, Xt, d2Tgc, dTgc, Tgc) + class(nurbs_curve), intent(inout) :: this + integer, intent(in), optional :: res + real(rk), intent(in), contiguous, optional :: Xt(:) + real(rk), allocatable, intent(out) :: d2Tgc(:,:) + real(rk), allocatable, intent(out), optional :: dTgc(:,:) + real(rk), allocatable, intent(out), optional :: Tgc(:,:) + integer :: i ! Set parameter values if (present(Xt)) then @@ -931,9 +1164,28 @@ pure function compute_dTgc_bspline_1d(f_Xt, f_knot, f_degree, f_nc, f_ng) result this%ng = size(this%Xt,1) if (this%is_rational()) then ! NURBS - dTgc = compute_dTgc_nurbs_1d(this%Xt, this%knot, this%degree, this%nc, this%ng, this%Wc) + call compute_d2Tgc(this%Xt, this%knot, this%degree, this%nc, this%ng, this%Wc, d2Tgc, dTgc, Tgc) + else ! B-Spline + call compute_d2Tgc(this%Xt, this%knot, this%degree, this%nc, this%ng, d2Tgc, dTgc, Tgc) + end if + end subroutine + !=============================================================================== + + + !=============================================================================== + !> author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + pure subroutine derivative2_scalar(this, Xt, d2Tgc, dTgc, Tgc) + class(nurbs_curve), intent(inout) :: this + real(rk), intent(in) :: Xt + real(rk), allocatable, intent(out) :: d2Tgc(:) + real(rk), allocatable, intent(out), optional :: dTgc(:) + real(rk), allocatable, intent(out), optional :: Tgc(:) + + if (this%is_rational()) then ! NURBS + call compute_d2Tgc(Xt, this%knot, this%degree, this%nc, this%Wc, d2Tgc, dTgc, Tgc) else ! B-Spline - dTgc = compute_dTgc_bspline_1d(this%Xt, this%knot, this%degree, this%nc, this%ng) + call compute_d2Tgc(Xt, this%knot, this%degree, this%nc, d2Tgc, dTgc, Tgc) end if end subroutine !=============================================================================== @@ -942,38 +1194,13 @@ pure function compute_dTgc_bspline_1d(f_Xt, f_knot, f_degree, f_nc, f_ng) result !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause - pure subroutine basis(this, res, Xt, Tgc) + pure subroutine basis_vector(this, res, Xt, Tgc) class(nurbs_curve), intent(inout) :: this integer, intent(in), optional :: res real(rk), intent(in), contiguous, optional :: Xt(:) real(rk), allocatable, intent(out) :: Tgc(:,:) integer :: i - interface - pure function compute_Tgc_nurbs_1d(f_Xt, f_knot, f_degree, f_nc, f_ng, f_Wc) result(f_Tgc) - import :: rk - real(rk), intent(in), contiguous :: f_Xt(:) - real(rk), intent(in), contiguous :: f_knot(:) - integer, intent(in) :: f_degree - integer, intent(in) :: f_nc - integer, intent(in) :: f_ng - real(rk), intent(in), contiguous :: f_Wc(:) - real(rk), allocatable :: f_Tgc(:,:) - end function - end interface - - interface - pure function compute_Tgc_bspline_1d(f_Xt, f_knot, f_degree, f_nc, f_ng) result(f_Tgc) - import :: rk - real(rk), intent(in), contiguous :: f_Xt(:) - real(rk), intent(in), contiguous :: f_knot(:) - integer, intent(in) :: f_degree - integer, intent(in) :: f_nc - integer, intent(in) :: f_ng - real(rk), allocatable :: f_Tgc(:,:) - end function - end interface - ! Set parameter values if (present(Xt)) then if (allocated(this%Xt)) deallocate(this%Xt) @@ -990,9 +1217,26 @@ pure function compute_Tgc_bspline_1d(f_Xt, f_knot, f_degree, f_nc, f_ng) result( this%ng = size(this%Xt,1) if (this%is_rational()) then ! NURBS - Tgc = compute_Tgc_nurbs_1d(this%Xt, this%knot, this%degree, this%nc, this%ng, this%Wc) + Tgc = compute_Tgc(this%Xt, this%knot, this%degree, this%nc, this%ng, this%Wc) else ! B-Spline - Tgc = compute_Tgc_bspline_1d(this%Xt, this%knot, this%degree, this%nc, this%ng) + Tgc = compute_Tgc(this%Xt, this%knot, this%degree, this%nc, this%ng) + end if + end subroutine + !=============================================================================== + + + !=============================================================================== + !> author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + pure subroutine basis_scalar(this, Xt, Tgc) + class(nurbs_curve), intent(inout) :: this + real(rk), intent(in) :: Xt + real(rk), allocatable, intent(out) :: Tgc(:) + + if (this%is_rational()) then ! NURBS + Tgc = compute_Tgc(Xt, this%knot, this%degree, this%nc, this%Wc) + else ! B-Spline + Tgc = compute_Tgc(Xt, this%knot, this%degree, this%nc) end if end subroutine !=============================================================================== @@ -1499,19 +1743,9 @@ pure subroutine nearest_point(this, point_Xg, nearest_Xg, nearest_Xt, id) integer :: id_ real(rk), allocatable :: distances(:) - interface - pure function nearest_point_help_1d(f_ng, f_Xg, f_point_Xg) result(f_distances) - import :: rk - integer, intent(in) :: f_ng - real(rk), intent(in), contiguous :: f_Xg(:,:) - real(rk), intent(in), contiguous :: f_point_Xg(:) - real(rk), allocatable :: f_distances(:) - end function - end interface - allocate(distances(this%ng)) distances = nearest_point_help_1d(this%ng, this%Xg, point_Xg) - + id_ = minloc(distances, dim=1) if (present(id)) id = id_ if (present(nearest_Xg)) nearest_Xg = this%Xg(id_,:) @@ -1519,6 +1753,92 @@ pure function nearest_point_help_1d(f_ng, f_Xg, f_point_Xg) result(f_distances) end subroutine !=============================================================================== + + !=============================================================================== + !> author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + impure subroutine nearest_point2(this, point_Xg, tol, maxit, nearest_Xt, nearest_Xg) + + class(nurbs_curve), intent(inout) :: this + real(rk), intent(in) :: point_Xg(:) + real(rk), intent(in) :: tol + integer, intent(in) :: maxit + real(rk), intent(out) :: nearest_Xt + real(rk), allocatable, intent(out), optional :: nearest_Xg(:) + real(rk):: xk, obj, grad, hess, dk, alphak, tau, beta, lower_bounds, upper_bounds + real(rk), allocatable :: Xg(:), Tgc(:), dTgc(:), d2Tgc(:), distances(:) + integer :: k, l + logical :: convergenz + type(nurbs_curve) :: copy_this + + k = 0 + + ! lower and upper bounds + lower_bounds = minval(this%knot) + upper_bounds = maxval(this%knot) + + ! guess initial point + copy_this = this + call this%create(10) + allocate(distances(copy_this%ng)) + distances = nearest_point_help_1d(copy_this%ng, copy_this%Xg, point_Xg) + xk = copy_this%Xt(minloc(distances, dim=1)) + call copy_this%finalize() + + ! Check if xk is within the knot vector range + if (xk < minval(this%knot)) then + xk = minval(this%knot) + else if (xk > maxval(this%knot)) then + xk = maxval(this%knot) + end if + + convergenz = .false. + + allocate(Xg(size(this%Xc,2))) + ! allocate(dTgc(size(this%Xc,1))) + ! allocate(d2Tgc(size(this%Xc,1))) + + do while (.not. convergenz .and. k < maxit) + + ! objection, gradient and hessian + Xg = this%cmp_Xg(xk) + call this%derivative2(Xt=xk, d2Tgc=d2Tgc, dTgc=dTgc, Tgc=Tgc) ! Tgc is not needed + + obj = norm2(Xg - point_Xg) + 0.001_rk ! add a small number to avoid division by zero + grad = dot_product((Xg-point_Xg)/obj, matmul(dTgc,this%Xc)) + hess = dot_product(matmul(dTgc,this%Xc) - (Xg-point_Xg)/obj*grad, matmul(dTgc,this%Xc))/obj& + + dot_product((Xg-point_Xg)/obj, matmul(d2Tgc,this%Xc)) + + ! debug + print '(i3,1x,e20.10,1x,e20.10)', k, xk, abs(grad) + + if (abs(grad) <= tol) then + convergenz = .true. + nearest_Xt = xk + if (present(nearest_Xg)) nearest_Xg = this%cmp_Xg(nearest_Xt) + else + dk = - grad / hess + + ! Backtracking-Armijo Line Search + alphak = 1.0_rk + tau = 0.5_rk ! 0 < tau < 1 + beta = 1.0e-4_rk ! 0 < beta < 1 + l = 0 + do while (.not. norm2(this%cmp_Xg(xk + alphak*dk) - point_Xg) <= obj + alphak*beta*grad*dk .and. l<50) + alphak = tau * alphak + l = l + 1 + end do + + xk = xk + alphak*dk + ! Check if xk is within the knot vector range + xk = max(min(xk, upper_bounds), lower_bounds) + k = k + 1 + end if + end do + + end subroutine + !=============================================================================== + end module forcad_nurbs_curve !=============================================================================== @@ -1552,6 +1872,31 @@ impure function compute_Xg_nurbs_1d(Xt, knot, degree, nc, ng, Xc, Wc) result(Xg) !=============================================================================== +!=============================================================================== +!> author: Seyed Ali Ghasemi +!> license: BSD 3-Clause +impure function compute_Xg_nurbs_1d_1point(Xt, knot, degree, nc, Xc, Wc) result(Xg) + use forcad_utils, only: rk, basis_bspline + + implicit none + real(rk), intent(in) :: Xt + real(rk), intent(in), contiguous :: knot(:) + integer, intent(in) :: degree + integer, intent(in) :: nc + real(rk), intent(in), contiguous :: Xc(:,:) + real(rk), intent(in), contiguous :: Wc(:) + real(rk), allocatable :: Xg(:) + real(rk), allocatable :: Tgc(:) + + allocate(Xg(size(Xc))) + allocate(Tgc(nc)) + Tgc = basis_bspline(Xt, knot, nc, degree) + Tgc = Tgc*(Wc/(dot_product(Tgc,Wc))) + Xg = matmul(Tgc, Xc) +end function +!=============================================================================== + + !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause @@ -1581,7 +1926,27 @@ impure function compute_Xg_bspline_1d(Xt, knot, degree, nc, ng, Xc) result(Xg) !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause -impure function compute_dTgc_nurbs_1d(Xt, knot, degree, nc, ng, Wc) result(dTgc) +impure function compute_Xg_bspline_1d_1point(Xt, knot, degree, nc, Xc) result(Xg) + use forcad_utils, only: rk, basis_bspline + + implicit none + real(rk), intent(in) :: Xt + real(rk), intent(in), contiguous :: knot(:) + integer, intent(in) :: degree + integer, intent(in) :: nc + real(rk), intent(in), contiguous :: Xc(:,:) + real(rk), allocatable :: Xg(:) + + allocate(Xg(size(Xc))) + Xg = matmul(basis_bspline(Xt, knot, nc, degree), Xc) +end function +!=============================================================================== + + +!=============================================================================== +!> author: Seyed Ali Ghasemi +!> license: BSD 3-Clause +impure subroutine compute_dTgc_nurbs_1d_vector(Xt, knot, degree, nc, ng, Wc, dTgc, Tgc) use forcad_utils, only: rk, basis_bspline_der implicit none @@ -1591,26 +1956,49 @@ impure function compute_dTgc_nurbs_1d(Xt, knot, degree, nc, ng, Wc) result(dTgc) integer, intent(in) :: nc integer, intent(in) :: ng real(rk), intent(in), contiguous :: Wc(:) - real(rk), allocatable :: dTgc(:,:) - real(rk), allocatable :: dTgci(:) + real(rk), allocatable, intent(out) :: dTgc(:,:) + real(rk), allocatable, intent(out) :: Tgc(:,:) + real(rk), allocatable :: dBi(:), Bi(:) integer :: i - allocate(dTgc(ng, nc)) - allocate(dTgci(nc)) - !$OMP PARALLEL DO PRIVATE(dTgci) + allocate(dTgc(ng, nc), Tgc(ng, nc), dBi(nc), Bi(nc)) do i = 1, size(Xt) - dTgci = basis_bspline_der(Xt(i), knot, nc, degree) - dTgc(i,:) = dTgci*(Wc/(dot_product(dTgci,Wc))) + call basis_bspline_der(Xt(i), knot, nc, degree, dBi, Bi) + Tgc(i,:) = Bi*(Wc/(dot_product(Bi,Wc))) + dTgc(i,:) = ( dBi*Wc - Tgc(i,:)*dot_product(dBi,Wc) ) / dot_product(Bi,Wc) end do - !$OMP END PARALLEL DO -end function +end subroutine !=============================================================================== !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause -impure function compute_dTgc_bspline_1d(Xt, knot, degree, nc, ng) result(dTgc) +impure subroutine compute_dTgc_nurbs_1d_scalar(Xt, knot, degree, nc, Wc, dTgc, Tgc) + use forcad_utils, only: rk, basis_bspline_der + + implicit none + real(rk), intent(in) :: Xt + real(rk), intent(in), contiguous :: knot(:) + integer, intent(in) :: degree + integer, intent(in) :: nc + real(rk), intent(in), contiguous :: Wc(:) + real(rk), allocatable, intent(out) :: dTgc(:) + real(rk), allocatable, intent(out) :: Tgc(:) + real(rk), allocatable :: dBi(:), Bi(:) + + allocate(dTgc(nc), Tgc(nc)) + call basis_bspline_der(Xt, knot, nc, degree, dBi, Bi) + Tgc = Bi*(Wc/(dot_product(Bi,Wc))) + dTgc = ( dBi*Wc - Tgc*dot_product(dBi,Wc) ) / dot_product(Bi,Wc) +end subroutine +!=============================================================================== + + +!=============================================================================== +!> author: Seyed Ali Ghasemi +!> license: BSD 3-Clause +impure subroutine compute_dTgc_bspline_1d_vector(Xt, knot, degree, nc, ng, dTgc, Tgc) use forcad_utils, only: rk, basis_bspline_der implicit none @@ -1619,23 +2007,157 @@ impure function compute_dTgc_bspline_1d(Xt, knot, degree, nc, ng) result(dTgc) integer, intent(in) :: degree integer, intent(in) :: nc integer, intent(in) :: ng - real(rk), allocatable :: dTgc(:,:) + real(rk), allocatable, intent(out) :: dTgc(:,:) + real(rk), allocatable, intent(out) :: Tgc(:,:) + real(rk), allocatable :: dBi(:), Bi(:) integer :: i - allocate(dTgc(ng, nc)) - !$OMP PARALLEL DO + allocate(dTgc(ng, nc), Tgc(ng, nc), dBi(nc), Bi(nc)) do i = 1, size(Xt) - dTgc(i,:) = basis_bspline_der(Xt(i), knot, nc, degree) + call basis_bspline_der(Xt(i), knot, nc, degree, dBi, Bi) + Tgc(i,:) = Bi + dTgc(i,:) = dBi end do - !$OMP END PARALLEL DO -end function +end subroutine +!=============================================================================== + + +!=============================================================================== +!> author: Seyed Ali Ghasemi +!> license: BSD 3-Clause +impure subroutine compute_dTgc_bspline_1d_scalar(Xt, knot, degree, nc, dTgc, Tgc) + use forcad_utils, only: rk, basis_bspline_der + + implicit none + real(rk), intent(in) :: Xt + real(rk), intent(in), contiguous :: knot(:) + integer, intent(in) :: degree + integer, intent(in) :: nc + real(rk), allocatable, intent(out) :: dTgc(:) + real(rk), allocatable, intent(out) :: Tgc(:) + + allocate(dTgc(nc), Tgc(nc)) + call basis_bspline_der(Xt, knot, nc, degree, dTgc, Tgc) +end subroutine !=============================================================================== !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause -impure function compute_Tgc_nurbs_1d(Xt, knot, degree, nc, ng, Wc) result(Tgc) +impure subroutine compute_d2Tgc_nurbs_1d_vector(Xt, knot, degree, nc, ng, Wc, d2Tgc, dTgc, Tgc) + use forcad_utils, only: rk, basis_bspline_2der + + implicit none + real(rk), intent(in), contiguous :: Xt(:) + real(rk), intent(in), contiguous :: knot(:) + integer, intent(in) :: degree + integer, intent(in) :: nc + integer, intent(in) :: ng + real(rk), intent(in), contiguous :: Wc(:) + real(rk), allocatable, intent(out) :: d2Tgc(:,:) + real(rk), allocatable, intent(out) :: dTgc(:,:) + real(rk), allocatable, intent(out) :: Tgc(:,:) + real(rk), allocatable :: d2Bi(:), dBi(:), Tgci(:), dTgci(:), Bi(:) + integer :: i + + allocate(d2Tgc(ng, nc), dTgc(ng, nc), Tgc(ng, nc), d2Bi(nc), dTgci(nc), dBi(nc), Tgci(nc), Bi(nc)) + + do i = 1, size(Xt) + call basis_bspline_2der(Xt(i), knot, nc, degree, d2Bi, dBi, Bi) + Tgci = Bi*(Wc/(dot_product(Bi,Wc))) + Tgc(i,:) = Tgci + + dTgci = ( dBi*Wc - Tgci*dot_product(dBi,Wc) ) / dot_product(Bi,Wc) + dTgc(i,:) = dTgci + + d2Tgc(i,:) = (d2Bi*Wc - 2.0_rk*dTgci*dot_product(dBi,Wc) - Tgci*dot_product(d2Bi,Wc)) / dot_product(Bi,Wc) + end do +end subroutine +!=============================================================================== + + +!=============================================================================== +!> author: Seyed Ali Ghasemi +!> license: BSD 3-Clause +impure subroutine compute_d2Tgc_nurbs_1d_scalar(Xt, knot, degree, nc, Wc, d2Tgc, dTgc, Tgc) + use forcad_utils, only: rk, basis_bspline_2der + + implicit none + real(rk), intent(in) :: Xt + real(rk), intent(in), contiguous :: knot(:) + integer, intent(in) :: degree + integer, intent(in) :: nc + real(rk), intent(in), contiguous :: Wc(:) + real(rk), allocatable, intent(out) :: d2Tgc(:) + real(rk), allocatable, intent(out) :: dTgc(:) + real(rk), allocatable, intent(out) :: Tgc(:) + real(rk), allocatable :: d2Bi(:), dBi(:), Bi(:) + + allocate(d2Tgc(nc), dTgc(nc), Tgc(nc), d2Bi(nc), dBi(nc), Bi(nc)) + + call basis_bspline_2der(Xt, knot, nc, degree, d2Bi, dBi, Bi) + Tgc = Bi*(Wc/(dot_product(Bi,Wc))) + dTgc = ( dBi*Wc - Tgc*dot_product(dBi,Wc) ) / dot_product(Bi,Wc) + d2Tgc = (d2Bi*Wc - 2.0_rk*dTgc*dot_product(dBi,Wc) - Tgc*dot_product(d2Bi,Wc)) / dot_product(Bi,Wc) +end subroutine +!=============================================================================== + + +!=============================================================================== +!> author: Seyed Ali Ghasemi +!> license: BSD 3-Clause +impure subroutine compute_d2Tgc_bspline_1d_vector(Xt, knot, degree, nc, ng, d2Tgc, dTgc, Tgc) + use forcad_utils, only: rk, basis_bspline_2der + + implicit none + real(rk), intent(in), contiguous :: Xt(:) + real(rk), intent(in), contiguous :: knot(:) + integer, intent(in) :: degree + integer, intent(in) :: nc + integer, intent(in) :: ng + real(rk), allocatable, intent(out) :: d2Tgc(:,:) + real(rk), allocatable, intent(out) :: dTgc(:,:) + real(rk), allocatable, intent(out) :: Tgc(:,:) + real(rk), allocatable :: d2Tgci(:), dTgci(:), Tgci(:) + integer :: i + + allocate(d2Tgc(ng, nc), dTgc(ng, nc), Tgc(ng, nc), dTgci(nc), Tgci(nc), d2Tgci(nc)) + do i = 1, size(Xt) + call basis_bspline_2der(Xt(i), knot, nc, degree, d2Tgci, dTgci, Tgci) + Tgc(i,:) = Tgci + dTgc(i,:) = dTgci + d2Tgc(i,:) = d2Tgci + end do +end subroutine +!=============================================================================== + + +!=============================================================================== +!> author: Seyed Ali Ghasemi +!> license: BSD 3-Clause +impure subroutine compute_d2Tgc_bspline_1d_scalar(Xt, knot, degree, nc, d2Tgc, dTgc, Tgc) + use forcad_utils, only: rk, basis_bspline_2der + + implicit none + real(rk), intent(in) :: Xt + real(rk), intent(in), contiguous :: knot(:) + integer, intent(in) :: degree + integer, intent(in) :: nc + real(rk), allocatable, intent(out) :: d2Tgc(:) + real(rk), allocatable, intent(out) :: dTgc(:) + real(rk), allocatable, intent(out) :: Tgc(:) + + allocate(d2Tgc(nc), dTgc(nc), Tgc(nc)) + call basis_bspline_2der(Xt, knot, nc, degree, d2Tgc, dTgc, Tgc) +end subroutine +!=============================================================================== + + +!=============================================================================== +!> author: Seyed Ali Ghasemi +!> license: BSD 3-Clause +impure function compute_Tgc_nurbs_1d_vector(Xt, knot, degree, nc, ng, Wc) result(Tgc) use forcad_utils, only: rk, basis_bspline implicit none @@ -1649,8 +2171,7 @@ impure function compute_Tgc_nurbs_1d(Xt, knot, degree, nc, ng, Wc) result(Tgc) real(rk), allocatable :: Tgci(:) integer :: i - allocate(Tgc(ng, nc)) - allocate(Tgci(nc)) + allocate(Tgc(ng, nc), Tgci(nc)) !$OMP PARALLEL DO PRIVATE(Tgci) do i = 1, size(Xt,1) Tgci = basis_bspline(Xt(i), knot, nc, degree) @@ -1664,7 +2185,28 @@ impure function compute_Tgc_nurbs_1d(Xt, knot, degree, nc, ng, Wc) result(Tgc) !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause -impure function compute_Tgc_bspline_1d(Xt, knot, degree, nc, ng) result(Tgc) +impure function compute_Tgc_nurbs_1d_scalar(Xt, knot, degree, nc, Wc) result(Tgc) + use forcad_utils, only: rk, basis_bspline + + implicit none + real(rk), intent(in) :: Xt + real(rk), intent(in), contiguous :: knot(:) + integer, intent(in) :: degree + integer, intent(in) :: nc + real(rk), intent(in), contiguous :: Wc(:) + real(rk), allocatable :: Tgc(:) + + allocate(Tgc(nc)) + Tgc = basis_bspline(Xt, knot, nc, degree) + Tgc = Tgc*(Wc/(dot_product(Tgc,Wc))) +end function +!=============================================================================== + + +!=============================================================================== +!> author: Seyed Ali Ghasemi +!> license: BSD 3-Clause +impure function compute_Tgc_bspline_1d_vector(Xt, knot, degree, nc, ng) result(Tgc) use forcad_utils, only: rk, basis_bspline implicit none @@ -1686,6 +2228,25 @@ impure function compute_Tgc_bspline_1d(Xt, knot, degree, nc, ng) result(Tgc) !=============================================================================== +!=============================================================================== +!> author: Seyed Ali Ghasemi +!> license: BSD 3-Clause +impure function compute_Tgc_bspline_1d_scalar(Xt, knot, degree, nc) result(Tgc) + use forcad_utils, only: rk, basis_bspline + + implicit none + real(rk), intent(in) :: Xt + real(rk), intent(in), contiguous :: knot(:) + integer, intent(in) :: degree + integer, intent(in) :: nc + real(rk), allocatable :: Tgc(:) + + allocate(Tgc(nc)) + Tgc = basis_bspline(Xt, knot, nc, degree) +end function +!=============================================================================== + + !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause @@ -1707,4 +2268,4 @@ impure function nearest_point_help_1d(ng, Xg, point_Xg) result(distances) !$OMP END PARALLEL DO end function -!=============================================================================== \ No newline at end of file +!=============================================================================== diff --git a/src/forcad_nurbs_surface.f90 b/src/forcad_nurbs_surface.f90 index 1a7d776cb..259f4e524 100644 --- a/src/forcad_nurbs_surface.f90 +++ b/src/forcad_nurbs_surface.f90 @@ -36,6 +36,7 @@ module forcad_nurbs_surface procedure :: set3 !!> Set Bezier or Rational Bezier surface using control points and weights generic :: set => set1, set2, set3 !!> Set NURBS surface procedure :: create !!> Generate geometry points + procedure :: cmp_Xg !!> Compute geometry points procedure, private :: get_Xc_all !!> Get all control points procedure, private :: get_Xci !!> Get i-th control point procedure, private :: get_Xcid !!> Get i-th control point in a specific direction @@ -74,8 +75,15 @@ module forcad_nurbs_surface procedure :: get_continuity !!> Compute and return the continuity of the NURBS surface procedure :: cmp_nc !!> Compute number of required control points procedure :: get_nc !!> Get number of control points - procedure :: derivative !!> Compute the derivative of the NURBS surface - procedure :: basis !!> Compute the basis functions of the NURBS surface + procedure, private :: basis_vector !!> Compute the basis functions of the NURBS surface + procedure, private :: basis_scalar !!> Compute the basis functions of the NURBS surface + generic :: basis => basis_vector, basis_scalar !!> Compute the basis functions of the NURBS surface + procedure, private :: derivative_vector !!> Compute the derivative of the NURBS surface + procedure, private :: derivative_scalar !!> Compute the derivative of the NURBS surface + generic :: derivative => derivative_vector, derivative_scalar !!> Compute the derivative of the NURBS surface + procedure, private :: derivative2_vector !!> Compute the second derivative of the NURBS surface + procedure, private :: derivative2_scalar !!> Compute the second derivative of the NURBS surface + generic :: derivative2 => derivative2_vector, derivative2_scalar !!> Compute the second derivative 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 @@ -85,7 +93,8 @@ module forcad_nurbs_surface 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 surface + procedure :: nearest_point !!> Find the nearest point on the NURBS surface (Approximation) + procedure :: nearest_point2 !!> Find the nearest point on the NURBS surface (Minimization - Newton's method) ! Shapes procedure :: set_tetragon !!> Set a tetragon @@ -95,6 +104,202 @@ module forcad_nurbs_surface end type !=============================================================================== + interface compute_Xg + 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) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:,:) + real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:) + integer, intent(in) :: f_degree(2) + integer, intent(in) :: f_nc(2) + integer, intent(in) :: f_ng(2) + real(rk), intent(in), contiguous :: f_Xc(:,:) + real(rk), intent(in), contiguous :: f_Wc(:) + real(rk), allocatable :: f_Xg(:,:) + end function + + pure function compute_Xg_bspline_2d(f_Xt, f_knot1, f_knot2, f_degree, f_nc, f_ng, f_Xc) result(f_Xg) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:,:) + real(rk), intent(in), contiguous :: f_knot1(:) + real(rk), intent(in), contiguous :: f_knot2(:) + integer, intent(in) :: f_degree(2) + integer, intent(in) :: f_nc(2) + integer, intent(in) :: f_ng(2) + real(rk), intent(in), contiguous :: f_Xc(:,:) + real(rk), allocatable :: f_Xg(:,:) + end function + + pure function compute_Xg_nurbs_2d_1point(f_Xt, f_knot1, f_knot2, f_degree, f_nc, f_Xc, f_Wc) result(f_Xg) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:) + real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:) + integer, intent(in) :: f_degree(2) + integer, intent(in) :: f_nc(2) + real(rk), intent(in), contiguous :: f_Xc(:,:) + real(rk), intent(in), contiguous :: f_Wc(:) + real(rk), allocatable :: f_Xg(:) + end function + + pure function compute_Xg_bspline_2d_1point(f_Xt, f_knot1, f_knot2, f_degree, f_nc, f_Xc) result(f_Xg) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:) + real(rk), intent(in), contiguous :: f_knot1(:) + real(rk), intent(in), contiguous :: f_knot2(:) + integer, intent(in) :: f_degree(2) + integer, intent(in) :: f_nc(2) + real(rk), intent(in), contiguous :: f_Xc(:,:) + real(rk), allocatable :: f_Xg(:) + end function + end interface + + interface compute_dTgc + pure subroutine compute_dTgc_nurbs_2d_vector(f_Xt, f_knot1, f_knot2, f_degree, f_nc, f_ng, f_Wc, f_dTgc, f_Tgc) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:,:) + real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:) + integer, intent(in) :: f_degree(2) + integer, intent(in) :: f_nc(2) + integer, intent(in) :: f_ng(2) + real(rk), intent(in), contiguous :: f_Wc(:) + real(rk), allocatable, intent(out) :: f_dTgc(:,:,:) + real(rk), allocatable, intent(out) :: f_Tgc(:,:) + end subroutine + + pure subroutine compute_dTgc_bspline_2d_vector(f_Xt, f_knot1, f_knot2, f_degree, nc, f_ng, f_dTgc, f_Tgc) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:,:) + real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:) + integer, intent(in) :: f_degree(2) + integer, intent(in) :: nc(2) + integer, intent(in) :: f_ng(2) + real(rk), allocatable, intent(out) :: f_dTgc(:,:,:) + real(rk), allocatable, intent(out) :: f_Tgc(:,:) + end subroutine + + pure subroutine compute_dTgc_nurbs_2d_scalar(f_Xt, f_knot1, f_knot2, f_degree, f_nc, f_Wc, f_dTgc, f_Tgc) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:) + real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:) + integer, intent(in) :: f_degree(2) + integer, intent(in) :: f_nc(2) + real(rk), intent(in), contiguous :: f_Wc(:) + real(rk), allocatable, intent(out) :: f_dTgc(:,:) + real(rk), allocatable, intent(out) :: f_Tgc(:) + end subroutine + + pure subroutine compute_dTgc_bspline_2d_scalar(f_Xt, f_knot1, f_knot2, f_degree, nc, f_dTgc, f_Tgc) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:) + real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:) + integer, intent(in) :: f_degree(2) + integer, intent(in) :: nc(2) + real(rk), allocatable, intent(out) :: f_dTgc(:,:) + real(rk), allocatable, intent(out) :: f_Tgc(:) + end subroutine + end interface + + interface compute_d2Tgc + pure subroutine compute_d2Tgc_nurbs_2d_vector(f_Xt, f_knot1, f_knot2, f_degree, f_nc, f_ng, f_Wc, f_d2Tgc, f_dTgc, f_Tgc) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:,:) + real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:) + integer, intent(in) :: f_degree(2) + integer, intent(in) :: f_nc(2) + integer, intent(in) :: f_ng(2) + real(rk), intent(in), contiguous :: f_Wc(:) + real(rk), allocatable, intent(out) :: f_d2Tgc(:,:,:) + real(rk), allocatable, intent(out) :: f_dTgc(:,:,:) + real(rk), allocatable, intent(out) :: f_Tgc(:,:) + end subroutine + + pure subroutine compute_d2Tgc_bspline_2d_vector(f_Xt, f_knot1, f_knot2, f_degree, nc, f_ng, f_d2Tgc, f_dTgc, f_Tgc) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:,:) + real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:) + integer, intent(in) :: f_degree(2) + integer, intent(in) :: nc(2) + integer, intent(in) :: f_ng(2) + real(rk), allocatable, intent(out) :: f_d2Tgc(:,:,:) + real(rk), allocatable, intent(out) :: f_dTgc(:,:,:) + real(rk), allocatable, intent(out) :: f_Tgc(:,:) + end subroutine + + pure subroutine compute_d2Tgc_nurbs_2d_scalar(f_Xt, f_knot1, f_knot2, f_degree, f_nc, f_Wc, f_d2Tgc, f_dTgc, f_Tgc) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:) + real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:) + integer, intent(in) :: f_degree(2) + integer, intent(in) :: f_nc(2) + real(rk), intent(in), contiguous :: f_Wc(:) + real(rk), allocatable, intent(out) :: f_d2Tgc(:,:) + real(rk), allocatable, intent(out) :: f_dTgc(:,:) + real(rk), allocatable, intent(out) :: f_Tgc(:) + end subroutine + + pure subroutine compute_d2Tgc_bspline_2d_scalar(f_Xt, f_knot1, f_knot2, f_degree, nc, f_d2Tgc, f_dTgc, f_Tgc) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:) + real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:) + integer, intent(in) :: f_degree(2) + integer, intent(in) :: nc(2) + real(rk), allocatable, intent(out) :: f_d2Tgc(:,:) + real(rk), allocatable, intent(out) :: f_dTgc(:,:) + real(rk), allocatable, intent(out) :: f_Tgc(:) + end subroutine + end interface + + interface compute_Tgc + pure function compute_Tgc_nurbs_2d_vector(f_Xt, f_knot1, f_knot2, f_degree, f_nc, f_ng, f_Wc) result(f_Tgc) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:,:) + real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:) + integer, intent(in) :: f_degree(2) + integer, intent(in) :: f_nc(2) + integer, intent(in) :: f_ng(2) + real(rk), intent(in), contiguous :: f_Wc(:) + real(rk), allocatable :: f_Tgc(:,:) + end function + + pure function compute_Tgc_bspline_2d_vector(f_Xt, f_knot1, f_knot2, f_degree, f_nc, f_ng) result(f_Tgc) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:,:) + real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:) + integer, intent(in) :: f_degree(2) + integer, intent(in) :: f_nc(2) + integer, intent(in) :: f_ng(2) + real(rk), allocatable :: f_Tgc(:,:) + end function + + pure function compute_Tgc_nurbs_2d_scalar(f_Xt, f_knot1, f_knot2, f_degree, f_nc, f_Wc) result(f_Tgc) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:) + real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:) + integer, intent(in) :: f_degree(2) + integer, intent(in) :: f_nc(2) + real(rk), intent(in), contiguous :: f_Wc(:) + real(rk), allocatable :: f_Tgc(:) + end function + + pure function compute_Tgc_bspline_2d_scalar(f_Xt, f_knot1, f_knot2, f_degree, f_nc) result(f_Tgc) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:) + real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:) + integer, intent(in) :: f_degree(2) + integer, intent(in) :: f_nc(2) + real(rk), allocatable :: f_Tgc(:) + end function + end interface + + interface + pure function nearest_point_help_2d(f_ng, f_Xg, f_point_Xg) result(f_distances) + import :: rk + integer, intent(in) :: f_ng(2) + real(rk), intent(in), contiguous :: f_Xg(:,:) + real(rk), intent(in), contiguous :: f_point_Xg(:) + real(rk), allocatable :: f_distances(:) + end function + end interface + contains !=============================================================================== @@ -200,34 +405,6 @@ pure subroutine create(this, res1, res2, Xt1, Xt2, Xt) real(rk), contiguous, intent(in), optional :: Xt(:,:) integer :: i - 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) - import :: rk - real(rk), intent(in), contiguous :: f_Xt(:,:) - real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:) - integer, intent(in) :: f_degree(2) - integer, intent(in) :: f_nc(2) - integer, intent(in) :: f_ng(2) - real(rk), intent(in), contiguous :: f_Xc(:,:) - real(rk), intent(in), contiguous :: f_Wc(:) - real(rk), allocatable :: f_Xg(:,:) - end function - end interface - - interface - pure function compute_Xg_bspline_2d(f_Xt, f_knot1, f_knot2, f_degree, f_nc, f_ng, f_Xc) result(f_Xg) - import :: rk - real(rk), intent(in), contiguous :: f_Xt(:,:) - real(rk), intent(in), contiguous :: f_knot1(:) - real(rk), intent(in), contiguous :: f_knot2(:) - integer, intent(in) :: f_degree(2) - integer, intent(in) :: f_nc(2) - integer, intent(in) :: f_ng(2) - real(rk), intent(in), contiguous :: f_Xc(:,:) - real(rk), allocatable :: f_Xg(:,:) - end function - end interface - ! check if (.not.allocated(this%Xc)) then error stop 'Control points are not set.' @@ -276,16 +453,42 @@ pure function compute_Xg_bspline_2d(f_Xt, f_knot1, f_knot2, f_degree, f_nc, f_ng if (allocated(this%Xg)) deallocate(this%Xg) if (this%is_rational()) then ! NURBS - this%Xg = compute_Xg_nurbs_2d(& - this%Xt, this%knot1, this%knot2, this%degree, this%nc, this%ng, this%Xc, this%Wc) + this%Xg = compute_Xg(& + 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(& - this%Xt, this%knot1, this%knot2, this%degree, this%nc, this%ng, this%Xc) + this%Xg = compute_Xg(& + this%Xt, this%knot1, this%knot2, this%degree, this%nc, this%ng, this%Xc) end if end subroutine !=============================================================================== + !=============================================================================== + !> author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + pure function cmp_Xg(this, Xt) result(Xg) + class(nurbs_surface), intent(in) :: this + real(rk), contiguous, intent(in) :: Xt(:) + real(rk), allocatable :: Xg(:) + + ! check + if (.not.allocated(this%Xc)) then + error stop 'Control points are not set.' + end if + + if (.not.allocated(this%knot1) .or. .not.allocated(this%knot2)) then + error stop 'Knot vector(s) is/are not set.' + end if + + if (this%is_rational()) then ! NURBS + Xg = compute_Xg(Xt, this%knot1, this%knot2, this%degree, this%nc, this%Xc, this%Wc) + else ! B-Spline + Xg = compute_Xg(Xt, this%knot1, this%knot2, this%degree, this%nc, this%Xc) + end if + end function + !=============================================================================== + + !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause @@ -960,38 +1163,84 @@ pure function get_nc(this, dir) result(nc) !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause - pure subroutine derivative(this, res1, res2, Xt1, Xt2, dTgc) + pure subroutine derivative_vector(this, res1, res2, Xt1, Xt2, dTgc, Tgc) class(nurbs_surface), intent(inout) :: this integer, intent(in), optional :: res1, res2 real(rk), intent(in), contiguous, optional :: Xt1(:), Xt2(:) - real(rk), allocatable, intent(out) :: dTgc(:,:) + real(rk), allocatable, intent(out) :: dTgc(:,:,:) + real(rk), allocatable, intent(out), optional :: Tgc(:,:) integer :: i real(rk), allocatable :: Xt(:,:) - interface - pure function compute_dTgc_nurbs_2d(f_Xt, f_knot1, f_knot2, f_degree, f_nc, f_ng, f_Wc) result(f_dTgc) - import :: rk - real(rk), intent(in), contiguous :: f_Xt(:,:) - real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:) - integer, intent(in) :: f_degree(2) - integer, intent(in) :: f_nc(2) - integer, intent(in) :: f_ng(2) - real(rk), intent(in), contiguous :: f_Wc(:) - real(rk), allocatable :: f_dTgc(:,:) - end function - end interface - - interface - pure function compute_dTgc_bspline_2d(f_Xt, f_knot1, f_knot2, f_degree, nc, f_ng) result(f_dTgc) - import :: rk - real(rk), intent(in), contiguous :: f_Xt(:,:) - real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:) - integer, intent(in) :: f_degree(2) - integer, intent(in) :: nc(2) - integer, intent(in) :: f_ng(2) - real(rk), allocatable :: f_dTgc(:,:) - end function - end interface + ! Set parameter values + if (present(Xt1)) then + if (allocated(this%Xt1)) deallocate(this%Xt1) + this%Xt1 = Xt1 + elseif (present(res1)) then + if (allocated(this%Xt1)) deallocate(this%Xt1) + allocate(this%Xt1(res1)) + this%Xt1 = [(real(i-1, rk) / real(res1-1, rk), i=1, res1)] + ! else + ! this%Xt1 = this%Xt1 + end if + + ! Set parameter values + if (present(Xt2)) then + if (allocated(this%Xt2)) deallocate(this%Xt2) + this%Xt2 = Xt2 + elseif (present(res2)) then + if (allocated(this%Xt2)) deallocate(this%Xt2) + allocate(this%Xt2(res2)) + this%Xt2 = [(real(i-1, rk) / real(res2-1, rk), i=1, res2)] + ! else + ! this%Xt2 = this%Xt2 + end if + + ! 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) + + if (this%is_rational()) then ! NURBS + call compute_dTgc(Xt, this%knot1, this%knot2, this%degree, this%nc, this%ng, this%Wc, dTgc, Tgc) + else ! B-Spline + call compute_dTgc(Xt, this%knot1, this%knot2, this%degree, this%nc, this%ng, dTgc, Tgc) + end if + end subroutine + !=============================================================================== + + + !=============================================================================== + !> author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + pure subroutine derivative_scalar(this, Xt, dTgc, Tgc) + class(nurbs_surface), intent(inout) :: this + real(rk), intent(in), contiguous :: Xt(:) + real(rk), allocatable, intent(out) :: dTgc(:,:) + real(rk), allocatable, intent(out), optional :: Tgc(:) + + if (this%is_rational()) then ! NURBS + call compute_dTgc(Xt, this%knot1, this%knot2, this%degree, this%nc, this%Wc, dTgc, Tgc) + else ! B-Spline + call compute_dTgc(Xt, this%knot1, this%knot2, this%degree, this%nc, dTgc, Tgc) + end if + end subroutine + !=============================================================================== + + + !=============================================================================== + !> author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + pure subroutine derivative2_vector(this, res1, res2, Xt1, Xt2, d2Tgc, dTgc, Tgc) + class(nurbs_surface), intent(inout) :: this + integer, intent(in), optional :: res1, res2 + real(rk), intent(in), contiguous, optional :: Xt1(:), Xt2(:) + real(rk), allocatable, intent(out) :: d2Tgc(:,:,:) + real(rk), allocatable, intent(out), optional :: dTgc(:,:,:) + real(rk), allocatable, intent(out), optional :: Tgc(:,:) + integer :: i + real(rk), allocatable :: Xt(:,:) ! Set parameter values if (present(Xt1)) then @@ -1024,9 +1273,28 @@ pure function compute_dTgc_bspline_2d(f_Xt, f_knot1, f_knot2, f_degree, nc, f_ng call ndgrid(this%Xt1, this%Xt2, Xt) if (this%is_rational()) then ! NURBS - dTgc = compute_dTgc_nurbs_2d(Xt, this%knot1, this%knot2, this%degree, this%nc, this%ng, this%Wc) + call compute_d2Tgc(Xt, this%knot1, this%knot2, this%degree, this%nc, this%ng, this%Wc, d2Tgc, dTgc, Tgc) + else ! B-Spline + call compute_d2Tgc(Xt, this%knot1, this%knot2, this%degree, this%nc, this%ng, d2Tgc, dTgc, Tgc) + end if + end subroutine + !=============================================================================== + + + !=============================================================================== + !> author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + pure subroutine derivative2_scalar(this, Xt, d2Tgc, dTgc, Tgc) + class(nurbs_surface), intent(inout) :: this + real(rk), intent(in), contiguous :: Xt(:) + real(rk), allocatable, intent(out) :: d2Tgc(:,:) + real(rk), allocatable, intent(out), optional :: dTgc(:,:) + real(rk), allocatable, intent(out), optional :: Tgc(:) + + if (this%is_rational()) then ! NURBS + call compute_d2Tgc(Xt, this%knot1, this%knot2, this%degree, this%nc, this%Wc, d2Tgc, dTgc, Tgc) else ! B-Spline - dTgc = compute_dTgc_bspline_2d(Xt, this%knot1, this%knot2, this%degree, this%nc, this%ng) + call compute_d2Tgc(Xt, this%knot1, this%knot2, this%degree, this%nc, d2Tgc, dTgc, Tgc) end if end subroutine !=============================================================================== @@ -1035,7 +1303,7 @@ pure function compute_dTgc_bspline_2d(f_Xt, f_knot1, f_knot2, f_degree, nc, f_ng !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause - pure subroutine basis(this, res1, res2, Xt1, Xt2, Tgc) + pure subroutine basis_vector(this, res1, res2, Xt1, Xt2, Tgc) class(nurbs_surface), intent(inout) :: this integer, intent(in), optional :: res1, res2 real(rk), intent(in), contiguous, optional :: Xt1(:), Xt2(:) @@ -1043,31 +1311,6 @@ pure subroutine basis(this, res1, res2, Xt1, Xt2, Tgc) integer :: i real(rk), allocatable :: Xt(:,:) - interface - pure function compute_Tgc_nurbs_2d(f_Xt, f_knot1, f_knot2, f_degree, f_nc, f_ng, f_Wc) result(f_Tgc) - import :: rk - real(rk), intent(in), contiguous :: f_Xt(:,:) - real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:) - integer, intent(in) :: f_degree(2) - integer, intent(in) :: f_nc(2) - integer, intent(in) :: f_ng(2) - real(rk), intent(in), contiguous :: f_Wc(:) - real(rk), allocatable :: f_Tgc(:,:) - end function - end interface - - interface - pure function compute_Tgc_bspline_2d(f_Xt, f_knot1, f_knot2, f_degree, f_nc, f_ng) result(f_Tgc) - import :: rk - real(rk), intent(in), contiguous :: f_Xt(:,:) - real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:) - integer, intent(in) :: f_degree(2) - integer, intent(in) :: f_nc(2) - integer, intent(in) :: f_ng(2) - real(rk), allocatable :: f_Tgc(:,:) - end function - end interface - ! Set parameter values if (present(Xt1)) then if (allocated(this%Xt1)) deallocate(this%Xt1) @@ -1099,9 +1342,26 @@ pure function compute_Tgc_bspline_2d(f_Xt, f_knot1, f_knot2, f_degree, f_nc, f_n call ndgrid(this%Xt1, this%Xt2, Xt) if (this%is_rational()) then ! NURBS - Tgc = compute_Tgc_nurbs_2d(Xt, this%knot1, this%knot2, this%degree, this%nc, this%ng, this%Wc) + Tgc = compute_Tgc(Xt, this%knot1, this%knot2, this%degree, this%nc, this%ng, this%Wc) + else ! B-Spline + Tgc = compute_Tgc(Xt, this%knot1, this%knot2, this%degree, this%nc, this%ng) + end if + end subroutine + !=============================================================================== + + + !=============================================================================== + !> author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + pure subroutine basis_scalar(this, Xt, Tgc) + class(nurbs_surface), intent(inout) :: this + real(rk), intent(in) :: Xt(:) + real(rk), allocatable, intent(out) :: Tgc(:) + + if (this%is_rational()) then ! NURBS + Tgc = compute_Tgc(Xt, this%knot1, this%knot2, this%degree, this%nc, this%Wc) else ! B-Spline - Tgc = compute_Tgc_bspline_2d(Xt, this%knot1, this%knot2, this%degree, this%nc, this%ng) + Tgc = compute_Tgc(Xt, this%knot1, this%knot2, this%degree, this%nc) end if end subroutine !=============================================================================== @@ -2099,19 +2359,9 @@ pure subroutine nearest_point(this, point_Xg, nearest_Xg, nearest_Xt, id) integer :: id_ real(rk), allocatable :: distances(:) - interface - pure function nearest_point_help_2d(f_ng, f_Xg, f_point_Xg) result(f_distances) - import :: rk - integer, intent(in) :: f_ng(2) - real(rk), intent(in), contiguous :: f_Xg(:,:) - real(rk), intent(in), contiguous :: f_point_Xg(:) - real(rk), allocatable :: f_distances(:) - end function - end interface - allocate(distances(this%ng(1)*this%ng(2))) distances = nearest_point_help_2d(this%ng, this%Xg, point_Xg) - + id_ = minloc(distances, dim=1) if (present(id)) id = id_ if (present(nearest_Xg)) nearest_Xg = this%Xg(id_,:) @@ -2119,6 +2369,120 @@ pure function nearest_point_help_2d(f_ng, f_Xg, f_point_Xg) result(f_distances) end subroutine !=============================================================================== + + !=============================================================================== + !> author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + impure subroutine nearest_point2(this, point_Xg, tol, maxit, nearest_Xt, nearest_Xg) + + class(nurbs_surface), intent(inout) :: this + real(rk), intent(in) :: point_Xg(:) + real(rk), intent(in) :: tol + integer, intent(in) :: maxit + real(rk), intent(out) :: nearest_Xt(2) + real(rk), allocatable, intent(out), optional :: nearest_Xg(:) + real(rk):: xk(2), obj, grad(2), hess(2,2), dk(2), alphak, tau, beta, det_inv, Ainv(2,2), lower_bounds(2), upper_bounds(2) + real(rk), allocatable :: Xg(:), Tgc(:), dTgc(:,:), d2Tgc(:,:), distances(:) + integer :: k, l + logical :: convergenz + type(nurbs_surface) :: copy_this + + k = 0 + + ! lower and upper bounds + lower_bounds = [minval(this%knot1), minval(this%knot2)] + upper_bounds = [maxval(this%knot1), maxval(this%knot2)] + + ! guess initial point + copy_this = this + call this%create(10,10) + allocate(distances(copy_this%ng(1)*copy_this%ng(2))) + distances = nearest_point_help_2d(copy_this%ng, copy_this%Xg, point_Xg) + xk = copy_this%Xt(minloc(distances, dim=1),:) + call copy_this%finalize() + + ! Check if xk is within the knot vector range + if (xk(1) < minval(this%knot1)) then + xk(1) = minval(this%knot1) + else if (xk(1) > maxval(this%knot1)) then + xk(1) = maxval(this%knot1) + end if + + if (xk(2) < minval(this%knot2)) then + xk(2) = minval(this%knot2) + else if (xk(2) > maxval(this%knot2)) then + xk(2) = maxval(this%knot2) + end if + + convergenz = .false. + + allocate(Xg(size(this%Xc,2))) + ! allocate(dTgc(size(this%Xc,1), 2)) + ! allocate(d2Tgc(size(this%Xc,1), 2)) + + do while (.not. convergenz .and. k < maxit) + + ! objection, gradient and hessian + Xg = this%cmp_Xg(xk) + call this%derivative2(Xt=xk, d2Tgc=d2Tgc, dTgc=dTgc, Tgc=Tgc) ! Tgc is not needed + + obj = norm2(Xg - point_Xg) + 0.001_rk ! add a small number to avoid division by zero + + grad(1) = dot_product((Xg-point_Xg)/obj, matmul(dTgc(:,1),this%Xc)) + grad(2) = dot_product((Xg-point_Xg)/obj, matmul(dTgc(:,2),this%Xc)) + + hess(1,1) = ( dot_product(matmul(dTgc(:,1),this%Xc), matmul(dTgc(:,1),this%Xc)) + & + dot_product((Xg-point_Xg), matmul(d2Tgc(1:this%nc(1)*this%nc(2) ,1),this%Xc)) ) /obj & + - ( dot_product(Xg-point_Xg, matmul(dTgc(:,1), this%Xc))*grad(1) ) / obj**2 + hess(2,1) = ( dot_product(matmul(dTgc(:,1),this%Xc), matmul(dTgc(:,2),this%Xc)) + & + dot_product((Xg-point_Xg), matmul(d2Tgc(this%nc(1)*this%nc(2)+1:2*this%nc(1)*this%nc(2),1),this%Xc)) ) /obj & + - ( dot_product(Xg-point_Xg, matmul(dTgc(:,2), this%Xc))*grad(1) ) / obj**2 + hess(1,2) = ( dot_product(matmul(dTgc(:,2),this%Xc), matmul(dTgc(:,1),this%Xc)) + & + dot_product((Xg-point_Xg), matmul(d2Tgc(1:this%nc(1)*this%nc(2) ,2),this%Xc)) ) /obj & + - ( dot_product(Xg-point_Xg, matmul(dTgc(:,1), this%Xc))*grad(2) ) / obj**2 + hess(2,2) = ( dot_product(matmul(dTgc(:,2),this%Xc), matmul(dTgc(:,2),this%Xc)) + & + dot_product((Xg-point_Xg), matmul(d2Tgc(this%nc(1)*this%nc(2)+1:2*this%nc(1)*this%nc(2),2),this%Xc)) ) /obj & + - ( dot_product(Xg-point_Xg, matmul(dTgc(:,2), this%Xc))*grad(2) ) / obj**2 + + ! debug + print '(i3,1x,2e20.10,1x,e20.10)', k, xk, norm2(grad) + + if (norm2(grad) <= tol) then + convergenz = .true. + nearest_Xt = xk + if (present(nearest_Xg)) nearest_Xg = this%cmp_Xg(nearest_Xt) + else + + ! Inverse of Hessian + det_inv = 1.0_rk/(hess(1,1)*hess(2,2) - hess(1,2)*hess(2,1)) + Ainv(1,1) = hess(2,2) + Ainv(2,1) = -hess(2,1) + Ainv(1,2) = -hess(1,2) + Ainv(2,2) = hess(1,1) + Ainv = Ainv * det_inv + + dk = - matmul(Ainv, grad) + + ! Backtracking-Armijo Line Search + alphak = 1.0_rk + tau = 0.5_rk ! 0 < tau < 1 + beta = 1.0e-4_rk ! 0 < beta < 1 + l = 0 + do while (.not. norm2(this%cmp_Xg(xk + alphak*dk) - point_Xg) <= obj + alphak*beta*dot_product(grad,dk) .and. l<50) + alphak = tau * alphak + l = l + 1 + end do + + xk = xk + alphak*dk + ! Check if xk is within the knot vector range + xk = max(min(xk, upper_bounds), lower_bounds) + k = k + 1 + end if + end do + + end subroutine + !=============================================================================== + end module forcad_nurbs_surface !=============================================================================== @@ -2155,6 +2519,34 @@ impure function compute_Xg_nurbs_2d(Xt, knot1, knot2, degree, nc, ng, Xc, Wc) re !=============================================================================== +!=============================================================================== +!> author: Seyed Ali Ghasemi +!> license: BSD 3-Clause +impure function compute_Xg_nurbs_2d_1point(Xt, knot1, knot2, degree, nc, Xc, Wc) result(Xg) + use forcad_utils, only: rk, basis_bspline, kron + + implicit none + real(rk), intent(in), contiguous :: Xt(:) + real(rk), intent(in), contiguous :: knot1(:), knot2(:) + integer, intent(in) :: degree(2) + integer, intent(in) :: nc(2) + real(rk), intent(in), contiguous :: Xc(:,:) + real(rk), intent(in), contiguous :: Wc(:) + real(rk), allocatable :: Xg(:) + real(rk), allocatable :: Tgc(:) + + allocate(Xg(size(Xc,2))) + allocate(Tgc(nc(1)*nc(2))) + + Tgc = kron(& + basis_bspline(Xt(2), knot2, nc(2), degree(2)),& + basis_bspline(Xt(1), knot1, nc(1), degree(1))) + Tgc = Tgc*(Wc/(dot_product(Tgc,Wc))) + Xg= matmul(Tgc, Xc) +end function +!=============================================================================== + + !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause @@ -2187,7 +2579,30 @@ impure function compute_Xg_bspline_2d(Xt, knot1, knot2, degree, nc, ng, Xc) resu !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause -impure function compute_dTgc_nurbs_2d(Xt, knot1, knot2, degree, nc, ng, Wc) result(dTgc) +impure function compute_Xg_bspline_2d_1point(Xt, knot1, knot2, degree, nc, Xc) result(Xg) + use forcad_utils, only: rk, basis_bspline, kron + + implicit none + real(rk), intent(in), contiguous :: Xt(:) + real(rk), intent(in), contiguous :: knot1(:), knot2(:) + integer, intent(in) :: degree(2) + integer, intent(in) :: nc(2) + real(rk), intent(in), contiguous :: Xc(:,:) + real(rk), allocatable :: Xg(:) + + allocate(Xg(size(Xc,2))) + Xg= matmul(kron(& + basis_bspline(Xt(2), knot2, nc(2), degree(2)),& + basis_bspline(Xt(1), knot1, nc(1), degree(1))),& + Xc) +end function +!=============================================================================== + + +!=============================================================================== +!> author: Seyed Ali Ghasemi +!> license: BSD 3-Clause +impure subroutine compute_dTgc_nurbs_2d_vector(Xt, knot1, knot2, degree, nc, ng, Wc, dTgc, Tgc) use forcad_utils, only: rk, basis_bspline_der, kron implicit none @@ -2197,28 +2612,71 @@ impure function compute_dTgc_nurbs_2d(Xt, knot1, knot2, degree, nc, ng, Wc) resu integer, intent(in) :: nc(2) integer, intent(in) :: ng(2) real(rk), intent(in), contiguous :: Wc(:) - real(rk), allocatable :: dTgc(:,:) - real(rk), allocatable :: dTgci(:) + real(rk), allocatable, intent(out) :: dTgc(:,:,:) + real(rk), allocatable, intent(out) :: Tgc(:,:) + real(rk), allocatable :: dBi(:,:), dB1(:), dB2(:) + real(rk), allocatable :: Bi(:), B1(:), B2(:) integer :: i - allocate(dTgc(ng(1)*ng(2), nc(1)*nc(2))) - allocate(dTgci(nc(1)*nc(2))) - !$OMP PARALLEL DO PRIVATE(dTgci) + allocate(dTgc(ng(1)*ng(2), nc(1)*nc(2), 2), Tgc(ng(1)*ng(2), nc(1)*nc(2))) + allocate(Bi(nc(1)*nc(2)), dBi(nc(1)*nc(2), 2)) do i = 1, size(Xt, 1) - dTgci = kron(& - basis_bspline_der(Xt(i,2), knot2, nc(2), degree(2)),& - basis_bspline_der(Xt(i,1), knot1, nc(1), degree(1))) - dTgc(i,:) = dTgci*(Wc/(dot_product(dTgci,Wc))) + call basis_bspline_der(Xt(i,1), knot1, nc(1), degree(1), dB1, B1) + call basis_bspline_der(Xt(i,2), knot2, nc(2), degree(2), dB2, B2) + + Bi = kron(B2, B1) + + Tgc(i,:) = Bi*(Wc/(dot_product(Bi,Wc))) + + dBi(:,1) = kron(B2, dB1) + dBi(:,2) = kron(dB2, B1) + + dTgc(i,:,1) = ( dBi(:,1)*Wc - Tgc(i,:)*dot_product(dBi(:,1),Wc) ) / dot_product(Bi,Wc) + dTgc(i,:,2) = ( dBi(:,2)*Wc - Tgc(i,:)*dot_product(dBi(:,2),Wc) ) / dot_product(Bi,Wc) end do - !$OMP END PARALLEL DO -end function +end subroutine +!=============================================================================== + + +!=============================================================================== +!> author: Seyed Ali Ghasemi +!> license: BSD 3-Clause +impure subroutine compute_dTgc_nurbs_2d_scalar(Xt, knot1, knot2, degree, nc, Wc, dTgc, Tgc) + use forcad_utils, only: rk, basis_bspline_der, kron + + implicit none + real(rk), intent(in), contiguous :: Xt(:) + real(rk), intent(in), contiguous :: knot1(:), knot2(:) + integer, intent(in) :: degree(2) + integer, intent(in) :: nc(2) + real(rk), intent(in), contiguous :: Wc(:) + real(rk), allocatable, intent(out) :: dTgc(:,:) + real(rk), allocatable, intent(out) :: Tgc(:) + real(rk), allocatable :: dB1(:), dB2(:), dBi(:,:) + real(rk), allocatable :: B1(:), B2(:), Bi(:) + + allocate(dTgc(nc(1)*nc(2), 2), Tgc(nc(1)*nc(2))) + allocate(dBi(nc(1)*nc(2), 2), Bi(nc(1)*nc(2))) + + call basis_bspline_der(Xt(1), knot1, nc(1), degree(1), dB1, B1) + call basis_bspline_der(Xt(2), knot2, nc(2), degree(2), dB2, B2) + + Bi = kron(B2, B1) + Tgc = Bi*(Wc/(dot_product(Bi,Wc))) + + dBi(:,1) = kron(B2, dB1) + dBi(:,2) = kron(dB2, B1) + + dTgc(:,1) = ( dBi(:,1)*Wc - Tgc*dot_product(dBi(:,1),Wc) ) / dot_product(Bi,Wc) + dTgc(:,2) = ( dBi(:,2)*Wc - Tgc*dot_product(dBi(:,2),Wc) ) / dot_product(Bi,Wc) +end subroutine !=============================================================================== !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause -impure function compute_dTgc_bspline_2d(Xt, knot1, knot2, degree, nc, ng) result(dTgc) +impure subroutine compute_dTgc_bspline_2d_vector(Xt, knot1, knot2, degree, nc, ng, dTgc, Tgc) use forcad_utils, only: rk, basis_bspline_der, kron implicit none @@ -2227,25 +2685,263 @@ impure function compute_dTgc_bspline_2d(Xt, knot1, knot2, degree, nc, ng) result integer, intent(in) :: degree(2) integer, intent(in) :: nc(2) integer, intent(in) :: ng(2) - real(rk), allocatable :: dTgc(:,:) + real(rk), allocatable, intent(out) :: dTgc(:,:,:) + real(rk), allocatable, intent(out) :: Tgc(:,:) + real(rk), allocatable :: dB1(:), dB2(:) + real(rk), allocatable :: B1(:), B2(:) integer :: i - allocate(dTgc(ng(1)*ng(2), nc(1)*nc(2))) - !$OMP PARALLEL DO + allocate(dTgc(ng(1)*ng(2), nc(1)*nc(2), 2)) + !$OMP PARALLEL DO PRIVATE(dB1, dB2) do i = 1, size(Xt, 1) - dTgc(i,:) = kron(& - basis_bspline_der(Xt(i,2), knot2, nc(2), degree(2)),& - basis_bspline_der(Xt(i,1), knot1, nc(1), degree(1))) + call basis_bspline_der(Xt(i,1), knot1, nc(1), degree(1), dB1, B1) + call basis_bspline_der(Xt(i,2), knot2, nc(2), degree(2), dB2, B2) + Tgc(i,:) = kron(B2, B1) + + dTgc(i,:,1) = kron(B2, dB1) + dTgc(i,:,2) = kron(dB2, B1) end do !$OMP END PARALLEL DO -end function +end subroutine +!=============================================================================== + + +!=============================================================================== +!> author: Seyed Ali Ghasemi +!> license: BSD 3-Clause +impure subroutine compute_dTgc_bspline_2d_scalar(Xt, knot1, knot2, degree, nc, dTgc, Tgc) + use forcad_utils, only: rk, basis_bspline_der, kron + + implicit none + real(rk), intent(in), contiguous :: Xt(:) + real(rk), intent(in), contiguous :: knot1(:), knot2(:) + integer, intent(in) :: degree(2) + integer, intent(in) :: nc(2) + real(rk), allocatable, intent(out) :: dTgc(:,:) + real(rk), allocatable, intent(out) :: Tgc(:) + real(rk), allocatable :: dTgc1(:), dTgc2(:) + real(rk), allocatable :: Tgc1(:), Tgc2(:) + + allocate(dTgc(nc(1)*nc(2), 2)) + call basis_bspline_der(Xt(1), knot1, nc(1), degree(1), dTgc1, Tgc1) + call basis_bspline_der(Xt(2), knot2, nc(2), degree(2), dTgc2, Tgc2) + Tgc = kron(Tgc2, Tgc1) + + dTgc(:,1) = kron(Tgc2, dTgc1) + dTgc(:,2) = kron(dTgc2, Tgc1) +end subroutine !=============================================================================== !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause -impure function compute_Tgc_nurbs_2d(Xt, knot1, knot2, degree, nc, ng, Wc) result(Tgc) +impure subroutine compute_d2Tgc_nurbs_2d_vector(Xt, knot1, knot2, degree, nc, ng, Wc, d2Tgc, dTgc, Tgc) + use forcad_utils, only: rk, basis_bspline_2der, kron + + implicit none + real(rk), intent(in), contiguous :: Xt(:,:) + real(rk), intent(in), contiguous :: knot1(:), knot2(:) + integer, intent(in) :: degree(2) + integer, intent(in) :: nc(2) + integer, intent(in) :: ng(2) + real(rk), intent(in), contiguous :: Wc(:) + real(rk), allocatable, intent(out) :: d2Tgc(:,:,:) + real(rk), allocatable, intent(out) :: dTgc(:,:,:) + real(rk), allocatable, intent(out) :: Tgc(:,:) + real(rk), allocatable :: d2Bi(:,:), d2B1(:), d2B2(:) + real(rk), allocatable :: dBi(:,:), dB1(:), dB2(:) + real(rk), allocatable :: Bi(:), B1(:), B2(:) + real(rk), allocatable :: Tgci(:), dTgci(:) + + integer :: i + + allocate(d2Tgc(ng(1)*ng(2), 2*nc(1)*nc(2), 2)) + + allocate(B1(nc(1)), B2(nc(2))) + allocate(dB1(nc(1)), dB2(nc(2))) + allocate(d2B1(nc(1)), d2B2(nc(2))) + allocate(Bi(nc(1)*nc(2)), dBi(nc(1)*nc(2), 2), d2Bi(2*nc(1)*nc(2), 2)) + + allocate(Tgci(nc(1)*nc(2)), dTgci(nc(1)*nc(2))) + allocate(Tgc(ng(1)*ng(2), nc(1)*nc(2)), dTgc(ng(1)*ng(2), nc(1)*nc(2), 2)) + do i = 1, size(Xt, 1) + call basis_bspline_2der(Xt(i,1), knot1, nc(1), degree(1), d2B1, dB1, B1) + call basis_bspline_2der(Xt(i,2), knot2, nc(2), degree(2), d2B2, dB2, B2) + + Bi = kron(B2, B1) + + Tgc(i,:) = Bi*(Wc/(dot_product(Bi,Wc))) + + dBi(:,1) = kron(B2, dB1) + dBi(:,2) = kron(dB2, B1) + + dTgc(i,:,1) = ( dBi(:,1)*Wc - Tgc(i,:)*dot_product(dBi(:,1),Wc) ) / dot_product(Bi,Wc) + dTgc(i,:,2) = ( dBi(:,2)*Wc - Tgc(i,:)*dot_product(dBi(:,2),Wc) ) / dot_product(Bi,Wc) + + d2Bi(1:nc(1)*nc(2) ,1) = kron(B2, d2B1) + d2Bi(nc(1)*nc(2)+1:2*nc(1)*nc(2),1) = kron(dB2, dB1) + d2Bi(1:nc(1)*nc(2) ,2) = kron(dB2, dB1) + d2Bi(nc(1)*nc(2)+1:2*nc(1)*nc(2),2) = kron(d2B2, B1) + + d2Tgc(i,1:nc(1)*nc(2) ,1) = & + (d2Bi(1:nc(1)*nc(2) ,1)*Wc - 2.0_rk*dTgc(i,:,1)*dot_product(dBi(:,1),Wc) & + - Tgc(i,:)*dot_product(d2Bi(1:nc(1)*nc(2) ,1),Wc)) / dot_product(Bi,Wc) + d2Tgc(i,nc(1)*nc(2)+1:2*nc(1)*nc(2),1) = & + (d2Bi(nc(1)*nc(2)+1:2*nc(1)*nc(2),1)*Wc - dTgc(i,:,1)*dot_product(dBi(:,2),Wc) - dTgc(i,:,2)*dot_product(dBi(:,1),Wc) & + - Tgc(i,:)*dot_product(d2Bi(nc(1)*nc(2)+1:2*nc(1)*nc(2),1),Wc)) / dot_product(Bi,Wc) + d2Tgc(i,1:nc(1)*nc(2) ,2) = & + (d2Bi(1:nc(1)*nc(2) ,2)*Wc - dTgc(i,:,1)*dot_product(dBi(:,2),Wc) - dTgc(i,:,2)*dot_product(dBi(:,1),Wc) & + - Tgc(i,:)*dot_product(d2Bi(1:nc(1)*nc(2) ,2),Wc)) / dot_product(Bi,Wc) + d2Tgc(i,nc(1)*nc(2)+1:2*nc(1)*nc(2),2) = & + (d2Bi(nc(1)*nc(2)+1:2*nc(1)*nc(2),2)*Wc - 2.0_rk*dTgc(i,:,2)*dot_product(dBi(:,2),Wc) & + - Tgc(i,:)*dot_product(d2Bi(nc(1)*nc(2)+1:2*nc(1)*nc(2),2),Wc)) / dot_product(Bi,Wc) + end do +end subroutine +!=============================================================================== + + +!=============================================================================== +!> author: Seyed Ali Ghasemi +!> license: BSD 3-Clause +impure subroutine compute_d2Tgc_nurbs_2d_scalar(Xt, knot1, knot2, degree, nc, Wc, d2Tgc, dTgc, Tgc) + use forcad_utils, only: rk, basis_bspline_2der, kron + + implicit none + real(rk), intent(in), contiguous :: Xt(:) + real(rk), intent(in), contiguous :: knot1(:), knot2(:) + integer, intent(in) :: degree(2) + integer, intent(in) :: nc(2) + real(rk), intent(in), contiguous :: Wc(:) + real(rk), allocatable, intent(out) :: d2Tgc(:,:) + real(rk), allocatable, intent(out) :: dTgc(:,:) + real(rk), allocatable, intent(out) :: Tgc(:) + real(rk), allocatable :: d2Bi(:,:), d2B1(:), d2B2(:) + real(rk), allocatable :: dBi(:,:), dB1(:), dB2(:) + real(rk), allocatable :: Bi(:), B1(:), B2(:) + + allocate(B1(nc(1)), B2(nc(2))) + allocate(dB1(nc(1)), dB2(nc(2))) + allocate(d2B1(nc(1)), d2B2(nc(2))) + allocate(Bi(nc(1)*nc(2)), dBi(nc(1)*nc(2), 2), d2Bi(2*nc(1)*nc(2), 2)) + + allocate(Tgc(nc(1)*nc(2)), dTgc(nc(1)*nc(2), 2), d2Tgc(2*nc(1)*nc(2), 2)) + + call basis_bspline_2der(Xt(1), knot1, nc(1), degree(1), d2B1, dB1, B1) + call basis_bspline_2der(Xt(2), knot2, nc(2), degree(2), d2B2, dB2, B2) + + Bi = kron(B2, B1) + + Tgc = Bi*(Wc/(dot_product(Bi,Wc))) + + dBi(:,1) = kron(B2, dB1) + dBi(:,2) = kron(dB2, B1) + + dTgc(:,1) = ( dBi(:,1)*Wc - Tgc*dot_product(dBi(:,1),Wc) ) / dot_product(Bi,Wc) + dTgc(:,2) = ( dBi(:,2)*Wc - Tgc*dot_product(dBi(:,2),Wc) ) / dot_product(Bi,Wc) + + d2Bi(1:nc(1)*nc(2) ,1) = kron(B2, d2B1) + d2Bi(nc(1)*nc(2)+1:2*nc(1)*nc(2),1) = kron(dB2, dB1) + d2Bi(1:nc(1)*nc(2) ,2) = kron(dB2, dB1) + d2Bi(nc(1)*nc(2)+1:2*nc(1)*nc(2),2) = kron(d2B2, B1) + + d2Tgc(1:nc(1)*nc(2) ,1) = & + (d2Bi(1:nc(1)*nc(2) ,1)*Wc - 2.0_rk*dTgc(:,1)*dot_product(dBi(:,1),Wc) & + - Tgc*dot_product(d2Bi(1:nc(1)*nc(2) ,1),Wc)) / dot_product(Bi,Wc) + d2Tgc(nc(1)*nc(2)+1:2*nc(1)*nc(2),1) = & + (d2Bi(nc(1)*nc(2)+1:2*nc(1)*nc(2),1)*Wc - dTgc(:,1)*dot_product(dBi(:,2),Wc) - dTgc(:,2)*dot_product(dBi(:,1),Wc) & + - Tgc*dot_product(d2Bi(nc(1)*nc(2)+1:2*nc(1)*nc(2),1),Wc)) / dot_product(Bi,Wc) + d2Tgc(1:nc(1)*nc(2) ,2) = & + (d2Bi(1:nc(1)*nc(2) ,2)*Wc - dTgc(:,1)*dot_product(dBi(:,2),Wc) - dTgc(:,2)*dot_product(dBi(:,1),Wc) & + - Tgc*dot_product(d2Bi(1:nc(1)*nc(2) ,2),Wc)) / dot_product(Bi,Wc) + d2Tgc(nc(1)*nc(2)+1:2*nc(1)*nc(2),2) = & + (d2Bi(nc(1)*nc(2)+1:2*nc(1)*nc(2),2)*Wc - 2.0_rk*dTgc(:,2)*dot_product(dBi(:,2),Wc) & + - Tgc*dot_product(d2Bi(nc(1)*nc(2)+1:2*nc(1)*nc(2),2),Wc)) / dot_product(Bi,Wc) +end subroutine +!=============================================================================== + + +!=============================================================================== +!> author: Seyed Ali Ghasemi +!> license: BSD 3-Clause +impure subroutine compute_d2Tgc_bspline_2d_vector(Xt, knot1, knot2, degree, nc, ng, d2Tgc, dTgc, Tgc) + use forcad_utils, only: rk, basis_bspline_2der, kron + + implicit none + real(rk), intent(in), contiguous :: Xt(:,:) + real(rk), intent(in), contiguous :: knot1(:), knot2(:) + integer, intent(in) :: degree(2) + integer, intent(in) :: nc(2) + integer, intent(in) :: ng(2) + real(rk), allocatable, intent(out) :: d2Tgc(:,:,:) + real(rk), allocatable, intent(out) :: dTgc(:,:,:) + real(rk), allocatable, intent(out) :: Tgc(:,:) + real(rk), allocatable :: d2B1(:), d2B2(:) + real(rk), allocatable :: dB1(:), dB2(:) + real(rk), allocatable :: B1(:), B2(:) + integer :: i + + allocate(d2Tgc(ng(1)*ng(2), 2*nc(1)*nc(2), 2)) + allocate(dTgc(ng(1)*ng(2), nc(1)*nc(2), 2)) + allocate(Tgc(ng(1)*ng(2), nc(1)*nc(2))) + do i = 1, size(Xt, 1) + call basis_bspline_2der(Xt(i,1), knot1, nc(1), degree(1), d2B1, dB1, B1) + call basis_bspline_2der(Xt(i,2), knot2, nc(2), degree(2), d2B2, dB2, B2) + + Tgc(i,:) = kron(B2, B1) + + dTgc(i,:,1) = kron(B2, dB1) + dTgc(i,:,2) = kron(dB2, B1) + + d2Tgc(i,1:nc(1)*nc(2) ,1) = kron(B2, d2B1) + d2Tgc(i,nc(1)*nc(2)+1:2*nc(1)*nc(2),1) = kron(dB2, dB1) + d2Tgc(i,1:nc(1)*nc(2) ,2) = kron(dB2, dB1) + d2Tgc(i,nc(1)*nc(2)+1:2*nc(1)*nc(2),2) = kron(d2B2, B1) + end do +end subroutine +!=============================================================================== + + +!=============================================================================== +!> author: Seyed Ali Ghasemi +!> license: BSD 3-Clause +impure subroutine compute_d2Tgc_bspline_2d_scalar(Xt, knot1, knot2, degree, nc, d2Tgc, dTgc, Tgc) + use forcad_utils, only: rk, basis_bspline_2der, kron + + implicit none + real(rk), intent(in), contiguous :: Xt(:) + real(rk), intent(in), contiguous :: knot1(:), knot2(:) + integer, intent(in) :: degree(2) + integer, intent(in) :: nc(2) + real(rk), allocatable, intent(out) :: d2Tgc(:,:) + real(rk), allocatable, intent(out) :: dTgc(:,:) + real(rk), allocatable, intent(out) :: Tgc(:) + real(rk), allocatable :: d2B1(:), d2B2(:) + real(rk), allocatable :: dB1(:), dB2(:) + real(rk), allocatable :: B1(:), B2(:) + + allocate(d2Tgc(2*nc(1)*nc(2), 2)) + allocate(dTgc(nc(1)*nc(2), 2)) + allocate(Tgc(nc(1)*nc(2))) + call basis_bspline_2der(Xt(1), knot1, nc(1), degree(1), d2B1, dB1, B1) + call basis_bspline_2der(Xt(2), knot2, nc(2), degree(2), d2B2, dB2, B2) + Tgc = kron(B2, B1) + + dTgc(:,1) = kron(B2, dB1) + dTgc(:,2) = kron(dB2, B1) + + d2Tgc(1:nc(1)*nc(2) ,1) = kron(B2, d2B1) + d2Tgc(nc(1)*nc(2)+1:2*nc(1)*nc(2),1) = kron(dB2, dB1) + d2Tgc(1:nc(1)*nc(2) ,2) = kron(dB2, dB1) + d2Tgc(nc(1)*nc(2)+1:2*nc(1)*nc(2),2) = kron(d2B2, B1) +end subroutine +!=============================================================================== + + +!=============================================================================== +!> author: Seyed Ali Ghasemi +!> license: BSD 3-Clause +impure function compute_Tgc_nurbs_2d_vector(Xt, knot1, knot2, degree, nc, ng, Wc) result(Tgc) use forcad_utils, only: rk, basis_bspline, kron implicit none @@ -2264,8 +2960,8 @@ impure function compute_Tgc_nurbs_2d(Xt, knot1, knot2, degree, nc, ng, Wc) resul !$OMP PARALLEL DO PRIVATE(Tgci) do i = 1, size(Xt, 1) Tgci = kron(& - basis_bspline(Xt(i,2), knot2, nc(2), degree(2)),& - basis_bspline(Xt(i,1), knot1, nc(1), degree(1))) + basis_bspline(Xt(i,2), knot2, nc(2), degree(2)),& + basis_bspline(Xt(i,1), knot1, nc(1), degree(1))) Tgc(i,:) = Tgci*(Wc/(dot_product(Tgci,Wc))) end do !$OMP END PARALLEL DO @@ -2276,7 +2972,30 @@ impure function compute_Tgc_nurbs_2d(Xt, knot1, knot2, degree, nc, ng, Wc) resul !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause -impure function compute_Tgc_bspline_2d(Xt, knot1, knot2, degree, nc, ng) result(Tgc) +impure function compute_Tgc_nurbs_2d_scalar(Xt, knot1, knot2, degree, nc, Wc) result(Tgc) + use forcad_utils, only: rk, basis_bspline, kron + + implicit none + real(rk), intent(in), contiguous :: Xt(:) + real(rk), intent(in), contiguous :: knot1(:), knot2(:) + integer, intent(in) :: degree(2) + integer, intent(in) :: nc(2) + real(rk), intent(in), contiguous :: Wc(:) + real(rk), allocatable :: Tgc(:) + + allocate(Tgc(nc(1)*nc(2))) + Tgc = kron(& + basis_bspline(Xt(2), knot2, nc(2), degree(2)),& + basis_bspline(Xt(1), knot1, nc(1), degree(1))) + Tgc = Tgc*(Wc/(dot_product(Tgc,Wc))) +end function +!=============================================================================== + + +!=============================================================================== +!> author: Seyed Ali Ghasemi +!> license: BSD 3-Clause +impure function compute_Tgc_bspline_2d_vector(Xt, knot1, knot2, degree, nc, ng) result(Tgc) use forcad_utils, only: rk, basis_bspline, kron implicit none @@ -2292,14 +3011,35 @@ impure function compute_Tgc_bspline_2d(Xt, knot1, knot2, degree, nc, ng) result( !$OMP PARALLEL DO do i = 1, size(Xt, 1) Tgc(i,:) = kron(& - basis_bspline(Xt(i,2), knot2, nc(2), degree(2)),& - basis_bspline(Xt(i,1), knot1, nc(1), degree(1))) + basis_bspline(Xt(i,2), knot2, nc(2), degree(2)),& + basis_bspline(Xt(i,1), knot1, nc(1), degree(1))) end do !$OMP END PARALLEL DO end function !=============================================================================== +!=============================================================================== +!> author: Seyed Ali Ghasemi +!> license: BSD 3-Clause +impure function compute_Tgc_bspline_2d_scalar(Xt, knot1, knot2, degree, nc) result(Tgc) + use forcad_utils, only: rk, basis_bspline, kron + + implicit none + real(rk), intent(in), contiguous :: Xt(:) + real(rk), intent(in), contiguous :: knot1(:), knot2(:) + integer, intent(in) :: degree(2) + integer, intent(in) :: nc(2) + real(rk), allocatable :: Tgc(:) + + allocate(Tgc(nc(1)*nc(2))) + Tgc = kron(& + basis_bspline(Xt(2), knot2, nc(2), degree(2)),& + basis_bspline(Xt(1), knot1, nc(1), degree(1))) +end function +!=============================================================================== + + !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause diff --git a/src/forcad_nurbs_volume.f90 b/src/forcad_nurbs_volume.f90 index 174b2bb47..31bce847f 100644 --- a/src/forcad_nurbs_volume.f90 +++ b/src/forcad_nurbs_volume.f90 @@ -38,6 +38,7 @@ module forcad_nurbs_volume procedure :: set3 !!> Set Bezier or Rational Bezier volume using control points and weights generic :: set => set1, set2, set3 !!> Set NURBS volume procedure :: create !!> Generate geometry points + procedure :: cmp_Xg !!> Compute geometry points procedure, private :: get_Xc_all !!> Get all control points procedure, private :: get_Xci !!> Get i-th control point procedure, private :: get_Xcid !!> Get i-th control point in a specific direction @@ -75,8 +76,15 @@ module forcad_nurbs_volume procedure :: get_multiplicity !!> Compute and return the multiplicity of the knots procedure :: get_continuity !!> Compute and return the continuity of the NURBS volume procedure :: cmp_nc !!> Compute number of required control points - procedure :: derivative !!> Compute the derivative of the NURBS volume - procedure :: basis !!> Compute the basis functions of the NURBS volume + procedure, private :: basis_vector !!> Compute the basis functions of the NURBS volume + procedure, private :: basis_scalar !!> Compute the basis functions of the NURBS volume + generic :: basis => basis_vector, basis_scalar !!> Compute the basis functions of the NURBS volume + procedure, private :: derivative_vector !!> Compute the derivative of the NURBS volume + procedure, private :: derivative_scalar !!> Compute the derivative of the NURBS volume + generic :: derivative => derivative_vector, derivative_scalar !!> Compute the derivative of the NURBS volume + procedure, private :: derivative2_vector !!> Compute the second derivative of the NURBS volume + procedure, private :: derivative2_scalar !!> Compute the second derivative of the NURBS volume + generic :: derivative2 => derivative2_vector, derivative2_scalar !!> Compute the second derivative 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 @@ -87,7 +95,8 @@ 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 + procedure :: nearest_point !!> Find the nearest point on the NURBS volume (Approximation) + procedure :: nearest_point2 !!> Find the nearest point on the NURBS volume (Minimization - Newton's method) ! Faces procedure :: cmp_elemFace_Xc_vis !!> Compute faces of the control points @@ -103,6 +112,202 @@ module forcad_nurbs_volume end type !=============================================================================== + interface compute_Xg + 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) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:,:) + real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:), f_knot3(:) + integer, intent(in) :: f_degree(3) + integer, intent(in) :: f_nc(3) + integer, intent(in) :: f_ng(3) + real(rk), intent(in), contiguous :: f_Xc(:,:) + real(rk), intent(in), contiguous :: f_Wc(:) + real(rk), allocatable :: f_Xg(:,:) + end function + + pure function compute_Xg_bspline_3d(f_Xt, f_knot1, f_knot2, f_knot3, f_degree, f_nc, f_ng, f_Xc) result(f_Xg) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:,:) + real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:), f_knot3(:) + integer, intent(in) :: f_degree(3) + integer, intent(in) :: f_nc(3) + integer, intent(in) :: f_ng(3) + real(rk), intent(in), contiguous :: f_Xc(:,:) + real(rk), allocatable :: f_Xg(:,:) + end function + + pure function compute_Xg_nurbs_3d_1point(f_Xt, f_knot1, f_knot2, f_knot3, f_degree, f_nc, f_Xc, f_Wc) result(f_Xg) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:) + real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:), f_knot3(:) + integer, intent(in) :: f_degree(3) + integer, intent(in) :: f_nc(3) + real(rk), intent(in), contiguous :: f_Xc(:,:) + real(rk), intent(in), contiguous :: f_Wc(:) + real(rk), allocatable :: f_Xg(:) + end function + + pure function compute_Xg_bspline_3d_1point(f_Xt, f_knot1, f_knot2, f_knot3, f_degree, f_nc, f_Xc) result(f_Xg) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:) + real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:), f_knot3(:) + integer, intent(in) :: f_degree(3) + integer, intent(in) :: f_nc(3) + real(rk), intent(in), contiguous :: f_Xc(:,:) + real(rk), allocatable :: f_Xg(:) + end function + end interface + + interface compute_dTgc + pure subroutine compute_dTgc_nurbs_3d_vector(f_Xt, f_knot1, f_knot2, f_knot3, f_degree, f_nc, f_ng, f_Wc, f_dTgc, f_Tgc) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:,:) + real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:), f_knot3(:) + integer, intent(in) :: f_degree(3) + integer, intent(in) :: f_nc(3) + integer, intent(in) :: f_ng(3) + real(rk), intent(in), contiguous :: f_Wc(:) + real(rk), allocatable, intent(out) :: f_dTgc(:,:,:) + real(rk), allocatable, intent(out) :: f_Tgc(:,:) + end subroutine + + pure subroutine compute_dTgc_bspline_3d_vector(f_Xt, f_knot1, f_knot2, f_knot3, f_degree, f_nc, f_ng, f_dTgc, f_Tgc) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:,:) + real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:), f_knot3(:) + integer, intent(in) :: f_degree(3) + integer, intent(in) :: f_nc(3) + integer, intent(in) :: f_ng(3) + real(rk), allocatable, intent(out) :: f_dTgc(:,:,:) + real(rk), allocatable, intent(out) :: f_Tgc(:,:) + end subroutine + + pure subroutine compute_dTgc_nurbs_3d_scalar(f_Xt, f_knot1, f_knot2, f_knot3, f_degree, f_nc, f_Wc, f_dTgc, f_Tgc) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:) + real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:), f_knot3(:) + integer, intent(in) :: f_degree(3) + integer, intent(in) :: f_nc(3) + real(rk), intent(in), contiguous :: f_Wc(:) + real(rk), allocatable, intent(out) :: f_dTgc(:,:) + real(rk), allocatable, intent(out) :: f_Tgc(:) + end subroutine + + pure subroutine compute_dTgc_bspline_3d_scalar(f_Xt, f_knot1, f_knot2, f_knot3, f_degree, f_nc, f_dTgc, f_Tgc) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:) + real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:), f_knot3(:) + integer, intent(in) :: f_degree(3) + integer, intent(in) :: f_nc(3) + real(rk), allocatable, intent(out) :: f_dTgc(:,:) + real(rk), allocatable, intent(out) :: f_Tgc(:) + end subroutine + end interface + + interface compute_d2Tgc + pure subroutine compute_d2Tgc_nurbs_3d_vector(& + f_Xt, f_knot1, f_knot2, f_knot3, f_degree, f_nc, f_ng, f_Wc, f_d2Tgc, f_dTgc, f_Tgc) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:,:) + real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:), f_knot3(:) + integer, intent(in) :: f_degree(3) + integer, intent(in) :: f_nc(3) + integer, intent(in) :: f_ng(3) + real(rk), intent(in), contiguous :: f_Wc(:) + real(rk), allocatable, intent(out) :: f_d2Tgc(:,:,:) + real(rk), allocatable, intent(out) :: f_dTgc(:,:,:) + real(rk), allocatable, intent(out) :: f_Tgc(:,:) + end subroutine + + pure subroutine compute_d2Tgc_bspline_3d_vector(& + f_Xt, f_knot1, f_knot2, f_knot3, f_degree, f_nc, f_ng, f_d2Tgc, f_dTgc, f_Tgc) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:,:) + real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:), f_knot3(:) + integer, intent(in) :: f_degree(3) + integer, intent(in) :: f_nc(3) + integer, intent(in) :: f_ng(3) + real(rk), allocatable, intent(out) :: f_d2Tgc(:,:,:) + real(rk), allocatable, intent(out) :: f_dTgc(:,:,:) + real(rk), allocatable, intent(out) :: f_Tgc(:,:) + end subroutine + + pure subroutine compute_d2Tgc_nurbs_3d_scalar(f_Xt, f_knot1, f_knot2, f_knot3, f_degree, f_nc, f_Wc, f_d2Tgc, f_dTgc, f_Tgc) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:) + real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:), f_knot3(:) + integer, intent(in) :: f_degree(3) + integer, intent(in) :: f_nc(3) + real(rk), intent(in), contiguous :: f_Wc(:) + real(rk), allocatable, intent(out) :: f_d2Tgc(:,:) + real(rk), allocatable, intent(out) :: f_dTgc(:,:) + real(rk), allocatable, intent(out) :: f_Tgc(:) + end subroutine + + pure subroutine compute_d2Tgc_bspline_3d_scalar(f_Xt, f_knot1, f_knot2, f_knot3, f_degree, f_nc, f_d2Tgc, f_dTgc, f_Tgc) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:) + real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:), f_knot3(:) + integer, intent(in) :: f_degree(3) + integer, intent(in) :: f_nc(3) + real(rk), allocatable, intent(out) :: f_d2Tgc(:,:) + real(rk), allocatable, intent(out) :: f_dTgc(:,:) + real(rk), allocatable, intent(out) :: f_Tgc(:) + end subroutine + end interface + + interface compute_Tgc + pure function compute_Tgc_nurbs_3d_vector(f_Xt, f_knot1, f_knot2, f_knot3, f_degree, f_nc, f_ng, f_Wc) result(f_Tgc) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:,:) + real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:), f_knot3(:) + integer, intent(in) :: f_degree(3) + integer, intent(in) :: f_nc(3) + integer, intent(in) :: f_ng(3) + real(rk), intent(in), contiguous :: f_Wc(:) + real(rk), allocatable :: f_Tgc(:,:) + end function + + pure function compute_Tgc_bspline_3d_vector(f_Xt, f_knot1, f_knot2, f_knot3, f_degree, f_nc, f_ng) result(f_Tgc) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:,:) + real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:), f_knot3(:) + integer, intent(in) :: f_degree(3) + integer, intent(in) :: f_nc(3) + integer, intent(in) :: f_ng(3) + real(rk), allocatable :: f_Tgc(:,:) + end function + + pure function compute_Tgc_nurbs_3d_scalar(f_Xt, f_knot1, f_knot2, f_knot3, f_degree, f_nc, f_Wc) result(f_Tgc) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:) + real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:), f_knot3(:) + integer, intent(in) :: f_degree(3) + integer, intent(in) :: f_nc(3) + real(rk), intent(in), contiguous :: f_Wc(:) + real(rk), allocatable :: f_Tgc(:) + end function + + pure function compute_Tgc_bspline_3d_scalar(f_Xt, f_knot1, f_knot2, f_knot3, f_degree, f_nc) result(f_Tgc) + import :: rk + real(rk), intent(in), contiguous :: f_Xt(:) + real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:), f_knot3(:) + integer, intent(in) :: f_degree(3) + integer, intent(in) :: f_nc(3) + real(rk), allocatable :: f_Tgc(:) + end function + end interface + + interface + pure function nearest_point_help_3d(f_ng, f_Xg, f_point_Xg) result(f_distances) + import :: rk + integer, intent(in) :: f_ng(3) + real(rk), intent(in), contiguous :: f_Xg(:,:) + real(rk), intent(in), contiguous :: f_point_Xg(:) + real(rk), allocatable :: f_distances(:) + end function + end interface + contains !=============================================================================== @@ -216,33 +421,6 @@ pure subroutine create(this, res1, res2, res3, Xt1, Xt2, Xt3, Xt) real(rk), intent(in), contiguous, optional :: Xt(:,:) integer :: i - 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) - import :: rk - real(rk), intent(in), contiguous :: f_Xt(:,:) - real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:), f_knot3(:) - integer, intent(in) :: f_degree(3) - integer, intent(in) :: f_nc(3) - integer, intent(in) :: f_ng(3) - real(rk), intent(in), contiguous :: f_Xc(:,:) - real(rk), intent(in), contiguous :: f_Wc(:) - real(rk), allocatable :: f_Xg(:,:) - end function - end interface - - interface - pure function compute_Xg_bspline_3d(f_Xt, f_knot1, f_knot2, f_knot3, f_degree, f_nc, f_ng, f_Xc) result(f_Xg) - import :: rk - real(rk), intent(in), contiguous :: f_Xt(:,:) - real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:), f_knot3(:) - integer, intent(in) :: f_degree(3) - integer, intent(in) :: f_nc(3) - integer, intent(in) :: f_ng(3) - real(rk), intent(in), contiguous :: f_Xc(:,:) - real(rk), allocatable :: f_Xg(:,:) - end function - end interface - ! check if (.not.allocated(this%Xc)) then error stop 'Control points are not set.' @@ -303,17 +481,43 @@ pure function compute_Xg_bspline_3d(f_Xt, f_knot1, f_knot2, f_knot3, f_degree, f 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(& - 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(& - this%Xt, this%knot1, this%knot2, this%knot3, this%degree, this%nc, this%ng, this%Xc) + if (this%is_rational()) then ! NURBS + this%Xg = compute_Xg(& + this%Xt, this%knot1, this%knot2, this%knot3, this%degree, this%nc, this%ng, this%Xc, this%Wc) + else ! B-Spline + this%Xg = compute_Xg(& + this%Xt, this%knot1, this%knot2, this%knot3, this%degree, this%nc, this%ng, this%Xc) end if end subroutine !=============================================================================== + !=============================================================================== + !> author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + pure function cmp_Xg(this, Xt) result(Xg) + class(nurbs_volume), intent(in) :: this + real(rk), intent(in), contiguous :: Xt(:) + real(rk), allocatable :: Xg(:) + + ! check + if (.not.allocated(this%Xc)) then + error stop 'Control points are not set.' + end if + + if (.not.allocated(this%knot1) .or. .not.allocated(this%knot2) .or. .not.allocated(this%knot3)) then + error stop 'Knot vector(s) is/are not set.' + end if + + if (this%is_rational()) then ! NURBS + Xg = compute_Xg(Xt, this%knot1, this%knot2, this%knot3, this%degree, this%nc, this%Xc, this%Wc) + else ! B-Spline + Xg = compute_Xg(Xt, this%knot1, this%knot2, this%knot3, this%degree, this%nc, this%Xc) + end if + end function + !=============================================================================== + + !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause @@ -1071,38 +1275,97 @@ pure function get_nc(this, dir) result(nc) !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause - pure subroutine derivative(this, res1, res2, res3, Xt1, Xt2, Xt3, dTgc) + pure subroutine derivative_vector(this, res1, res2, res3, Xt1, Xt2, Xt3, dTgc, Tgc) class(nurbs_volume), intent(inout) :: this integer, intent(in), optional :: res1, res2, res3 real(rk), intent(in), contiguous, optional :: Xt1(:), Xt2(:), Xt3(:) - real(rk), allocatable, intent(out) :: dTgc(:,:) + real(rk), allocatable, intent(out) :: dTgc(:,:,:) + real(rk), allocatable, intent(out), optional :: Tgc(:,:) integer :: i real(rk), allocatable :: Xt(:,:) - interface - pure function compute_dTgc_nurbs_3d(f_Xt, f_knot1, f_knot2, f_knot3, f_degree, f_nc, f_ng, f_Wc) result(f_dTgc) - import :: rk - real(rk), intent(in), contiguous :: f_Xt(:,:) - real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:), f_knot3(:) - integer, intent(in) :: f_degree(3) - integer, intent(in) :: f_nc(3) - integer, intent(in) :: f_ng(3) - real(rk), intent(in), contiguous :: f_Wc(:) - real(rk), allocatable :: f_dTgc(:,:) - end function - end interface - - interface - pure function compute_dTgc_bspline_3d(f_Xt, f_knot1, f_knot2, f_knot3, f_degree, f_nc, f_ng) result(f_dTgc) - import :: rk - real(rk), intent(in), contiguous :: f_Xt(:,:) - real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:), f_knot3(:) - integer, intent(in) :: f_degree(3) - integer, intent(in) :: f_nc(3) - integer, intent(in) :: f_ng(3) - real(rk), allocatable :: f_dTgc(:,:) - end function - end interface + ! Set parameter values + if (present(Xt1)) then + if (allocated(this%Xt1)) deallocate(this%Xt1) + this%Xt1 = Xt1 + elseif (present(res1)) then + if (allocated(this%Xt1)) deallocate(this%Xt1) + allocate(this%Xt1(res1)) + this%Xt1 = [(real(i-1, rk) / real(res1-1, rk), i=1, res1)] + ! else + ! this%Xt1 = this%Xt1 + end if + + ! Set parameter values + if (present(Xt2)) then + if (allocated(this%Xt2)) deallocate(this%Xt2) + this%Xt2 = Xt2 + elseif (present(res2)) then + if (allocated(this%Xt2)) deallocate(this%Xt2) + allocate(this%Xt2(res2)) + this%Xt2 = [(real(i-1, rk) / real(res2-1, rk), i=1, res2)] + ! else + ! this%Xt2 = this%Xt2 + end if + + ! Set parameter values + if (present(Xt3)) then + if (allocated(this%Xt3)) deallocate(this%Xt3) + this%Xt3 = Xt3 + elseif (present(res3)) then + if (allocated(this%Xt3)) deallocate(this%Xt3) + allocate(this%Xt3(res3)) + this%Xt3 = [(real(i-1, rk) / real(res3-1, rk), i=1, res3)] + ! else + ! this%Xt3 = this%Xt3 + end if + + ! 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) + + if (this%is_rational()) then ! NURBS + call compute_dTgc(Xt, this%knot1, this%knot2, this%knot3, this%degree, this%nc, this%ng, this%Wc, dTgc, Tgc) + else + call compute_dTgc(Xt, this%knot1, this%knot2, this%knot3, this%degree, this%nc, this%ng, dTgc, Tgc) + end if + end subroutine + !=============================================================================== + + + !=============================================================================== + !> author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + pure subroutine derivative_scalar(this, Xt, dTgc, Tgc) + class(nurbs_volume), intent(inout) :: this + real(rk), intent(in), contiguous :: Xt(:) + real(rk), allocatable, intent(out) :: dTgc(:,:) + real(rk), allocatable, intent(out), optional :: Tgc(:) + + if (this%is_rational()) then ! NURBS + call compute_dTgc(Xt, this%knot1, this%knot2, this%knot3, this%degree, this%nc, this%Wc, dTgc, Tgc) + else + call compute_dTgc(Xt, this%knot1, this%knot2, this%knot3, this%degree, this%nc, dTgc, Tgc) + end if + end subroutine + !=============================================================================== + + + !=============================================================================== + !> author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + pure subroutine derivative2_vector(this, res1, res2, res3, Xt1, Xt2, Xt3, d2Tgc, dTgc, Tgc) + class(nurbs_volume), intent(inout) :: this + integer, intent(in), optional :: res1, res2, res3 + real(rk), intent(in), contiguous, optional :: Xt1(:), Xt2(:), Xt3(:) + real(rk), allocatable, intent(out) :: d2Tgc(:,:,:) + real(rk), allocatable, intent(out), optional :: dTgc(:,:,:) + real(rk), allocatable, intent(out), optional :: Tgc(:,:) + integer :: i + real(rk), allocatable :: Xt(:,:) ! Set parameter values if (present(Xt1)) then @@ -1148,9 +1411,9 @@ pure function compute_dTgc_bspline_3d(f_Xt, f_knot1, f_knot2, f_knot3, f_degree, call ndgrid(this%Xt1, this%Xt2, this%Xt3, Xt) if (this%is_rational()) then ! NURBS - dTgc = compute_dTgc_nurbs_3d(Xt, this%knot1, this%knot2, this%knot3, this%degree, this%nc, this%ng, this%Wc) + call compute_d2Tgc(Xt, this%knot1, this%knot2, this%knot3, this%degree, this%nc, this%ng, this%Wc, d2Tgc, dTgc, Tgc) else - dTgc = compute_dTgc_bspline_3d(Xt, this%knot1, this%knot2, this%knot3, this%degree, this%nc, this%ng) + call compute_d2Tgc(Xt, this%knot1, this%knot2, this%knot3, this%degree, this%nc, this%ng, d2Tgc, dTgc, Tgc) end if end subroutine !=============================================================================== @@ -1159,7 +1422,27 @@ pure function compute_dTgc_bspline_3d(f_Xt, f_knot1, f_knot2, f_knot3, f_degree, !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause - pure subroutine basis(this, res1, res2, res3, Xt1, Xt2, Xt3, Tgc) + pure subroutine derivative2_scalar(this, Xt, d2Tgc, dTgc, Tgc) + class(nurbs_volume), intent(inout) :: this + real(rk), intent(in), contiguous :: Xt(:) + real(rk), allocatable, intent(out) :: d2Tgc(:,:) + real(rk), allocatable, intent(out), optional :: dTgc(:,:) + real(rk), allocatable, intent(out), optional :: Tgc(:) + + if (this%is_rational()) then ! NURBS + call compute_d2Tgc(Xt, this%knot1, this%knot2, this%knot3, this%degree, this%nc, this%Wc, d2Tgc, dTgc, Tgc) + else + call compute_d2Tgc(Xt, this%knot1, this%knot2, this%knot3, this%degree, this%nc, d2Tgc, dTgc, Tgc) + end if + end subroutine + !=============================================================================== + + + + !=============================================================================== + !> author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + pure subroutine basis_vector(this, res1, res2, res3, Xt1, Xt2, Xt3, Tgc) class(nurbs_volume), intent(inout) :: this integer, intent(in), optional :: res1, res2, res3 real(rk), intent(in), contiguous, optional :: Xt1(:), Xt2(:), Xt3(:) @@ -1167,31 +1450,6 @@ pure subroutine basis(this, res1, res2, res3, Xt1, Xt2, Xt3, Tgc) integer :: i real(rk), allocatable :: Xt(:,:) - interface - pure function compute_Tgc_nurbs_3d(f_Xt, f_knot1, f_knot2, f_knot3, f_degree, f_nc, f_ng, f_Wc) result(f_Tgc) - import :: rk - real(rk), intent(in), contiguous :: f_Xt(:,:) - real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:), f_knot3(:) - integer, intent(in) :: f_degree(3) - integer, intent(in) :: f_nc(3) - integer, intent(in) :: f_ng(3) - real(rk), intent(in), contiguous :: f_Wc(:) - real(rk), allocatable :: f_Tgc(:,:) - end function - end interface - - interface - pure function compute_Tgc_bspline_3d(f_Xt, f_knot1, f_knot2, f_knot3, f_degree, f_nc, f_ng) result(f_Tgc) - import :: rk - real(rk), intent(in), contiguous :: f_Xt(:,:) - real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:), f_knot3(:) - integer, intent(in) :: f_degree(3) - integer, intent(in) :: f_nc(3) - integer, intent(in) :: f_ng(3) - real(rk), allocatable :: f_Tgc(:,:) - end function - end interface - ! Set parameter values if (present(Xt1)) then if (allocated(this%Xt1)) deallocate(this%Xt1) @@ -1235,10 +1493,27 @@ pure function compute_Tgc_bspline_3d(f_Xt, f_knot1, f_knot2, f_knot3, f_degree, call ndgrid(this%Xt1, this%Xt2, this%Xt3, Xt) - if (allocated(this%Wc)) then ! NURBS volume - Tgc = compute_Tgc_nurbs_3d(Xt, this%knot1, this%knot2, this%knot3, this%degree, this%nc, this%ng, this%Wc) + if (this%is_rational()) then ! NURBS + Tgc = compute_Tgc(Xt, this%knot1, this%knot2, this%knot3, this%degree, this%nc, this%ng, this%Wc) + else + Tgc = compute_Tgc(Xt, this%knot1, this%knot2, this%knot3, this%degree, this%nc, this%ng) + end if + end subroutine + !=============================================================================== + + + !=============================================================================== + !> author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + pure subroutine basis_scalar(this, Xt, Tgc) + class(nurbs_volume), intent(inout) :: this + real(rk), intent(in), contiguous :: Xt(:) + real(rk), allocatable, intent(out) :: Tgc(:) + + if (this%is_rational()) then ! NURBS + Tgc = compute_Tgc(Xt, this%knot1, this%knot2, this%knot3, this%degree, this%nc, this%Wc) else - Tgc = compute_Tgc_bspline_3d(Xt, this%knot1, this%knot2, this%knot3, this%degree, this%nc, this%ng) + Tgc = compute_Tgc(Xt, this%knot1, this%knot2, this%knot3, this%degree, this%nc) end if end subroutine !=============================================================================== @@ -2596,16 +2871,6 @@ pure subroutine nearest_point(this, point_Xg, nearest_Xg, nearest_Xt, id) integer :: id_ real(rk), allocatable :: distances(:) - interface - pure function nearest_point_help_3d(f_ng, f_Xg, f_point_Xg) result(f_distances) - import :: rk - integer, intent(in) :: f_ng(3) - real(rk), intent(in), contiguous :: f_Xg(:,:) - real(rk), intent(in), contiguous :: f_point_Xg(:) - real(rk), allocatable :: f_distances(:) - end function - end interface - allocate(distances(this%ng(1)*this%ng(2)*this%ng(3))) distances = nearest_point_help_3d(this%ng, this%Xg, point_Xg) @@ -2617,6 +2882,152 @@ pure function nearest_point_help_3d(f_ng, f_Xg, f_point_Xg) result(f_distances) !=============================================================================== + + !=============================================================================== + !> author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + impure subroutine nearest_point2(this, point_Xg, tol, maxit, nearest_Xt, nearest_Xg) + + class(nurbs_volume), intent(inout) :: this + real(rk), intent(in) :: point_Xg(:) + real(rk), intent(in) :: tol + integer, intent(in) :: maxit + real(rk), intent(out) :: nearest_Xt(3) + real(rk), allocatable, intent(out), optional :: nearest_Xg(:) + real(rk):: xk(3), obj, grad(3), hess(3,3), dk(3), alphak, tau, beta, det_inv, Ainv(3,3), lower_bounds(3), upper_bounds(3) + real(rk), allocatable :: Xg(:), Tgc(:), dTgc(:,:), d2Tgc(:,:), distances(:) + integer :: k, l + logical :: convergenz + type(nurbs_volume) :: copy_this + + k = 0 + + ! lower and upper bounds + lower_bounds = [minval(this%knot1), minval(this%knot2), minval(this%knot3)] + upper_bounds = [maxval(this%knot1), maxval(this%knot2), maxval(this%knot3)] + + ! guess initial point + copy_this = this + call this%create(50, 50, 50) + allocate(distances(copy_this%ng(1)*copy_this%ng(2)*copy_this%ng(3))) + distances = nearest_point_help_3d(copy_this%ng, copy_this%Xg, point_Xg) + xk = copy_this%Xt(minloc(distances, dim=1),:) + call copy_this%finalize() + + ! Check if xk is within the knot vector range + if (xk(1) < minval(this%knot1)) then + xk(1) = minval(this%knot1) + else if (xk(1) > maxval(this%knot1)) then + xk(1) = maxval(this%knot1) + end if + + if (xk(2) < minval(this%knot2)) then + xk(2) = minval(this%knot2) + else if (xk(2) > maxval(this%knot2)) then + xk(2) = maxval(this%knot2) + end if + + if (xk(3) < minval(this%knot3)) then + xk(3) = minval(this%knot3) + else if (xk(3) > maxval(this%knot3)) then + xk(3) = maxval(this%knot3) + end if + + convergenz = .false. + + allocate(Xg(size(this%Xc,2))) + ! allocate(dTgc(size(this%Xc,1), 2)) + ! allocate(d2Tgc(size(this%Xc,1), 2)) + + do while (.not. convergenz .and. k < maxit) + + ! objection, gradient and hessian + Xg = this%cmp_Xg(xk) + call this%derivative2(Xt=xk, d2Tgc=d2Tgc, dTgc=dTgc, Tgc=Tgc) ! Tgc is not needed + + obj = norm2(Xg - point_Xg) + 0.001_rk ! add a small number to avoid division by zero + + grad(1) = dot_product((Xg-point_Xg)/obj, matmul(dTgc(:,1),this%Xc)) + grad(2) = dot_product((Xg-point_Xg)/obj, matmul(dTgc(:,2),this%Xc)) + grad(3) = dot_product((Xg-point_Xg)/obj, matmul(dTgc(:,3),this%Xc)) + + hess(1,1) = ( dot_product(matmul(dTgc(:,1),this%Xc), matmul(dTgc(:,1),this%Xc)) + dot_product((Xg-point_Xg), & + matmul(d2Tgc(1:this%nc(1)*this%nc(2)*this%nc(3) ,1),this%Xc)) ) /obj & + - ( dot_product(Xg-point_Xg, matmul(dTgc(:,1), this%Xc))*grad(1) ) / obj**2 + hess(2,1) = ( dot_product(matmul(dTgc(:,1),this%Xc), matmul(dTgc(:,2),this%Xc)) + dot_product((Xg-point_Xg), & + matmul(d2Tgc(this%nc(1)*this%nc(2)*this%nc(3)+1:2*this%nc(1)*this%nc(2)*this%nc(3) ,1),this%Xc)) ) /obj & + - ( dot_product(Xg-point_Xg, matmul(dTgc(:,2), this%Xc))*grad(1) ) / obj**2 + hess(3,1) = ( dot_product(matmul(dTgc(:,1),this%Xc), matmul(dTgc(:,3),this%Xc)) + dot_product((Xg-point_Xg), & + matmul(d2Tgc(2*this%nc(1)*this%nc(2)*this%nc(3)+1:3*this%nc(1)*this%nc(2)*this%nc(3) ,1),this%Xc)) ) /obj & + - ( dot_product(Xg-point_Xg, matmul(dTgc(:,3), this%Xc))*grad(1) ) / obj**2 + hess(1,2) = ( dot_product(matmul(dTgc(:,2),this%Xc), matmul(dTgc(:,1),this%Xc)) + dot_product((Xg-point_Xg), & + matmul(d2Tgc(1:this%nc(1)*this%nc(2)*this%nc(3) ,2),this%Xc)) ) /obj & + - ( dot_product(Xg-point_Xg, matmul(dTgc(:,1), this%Xc))*grad(2) ) / obj**2 + hess(2,2) = ( dot_product(matmul(dTgc(:,2),this%Xc), matmul(dTgc(:,2),this%Xc)) + dot_product((Xg-point_Xg), & + matmul(d2Tgc(this%nc(1)*this%nc(2)*this%nc(3)+1:2*this%nc(1)*this%nc(2)*this%nc(3) ,2),this%Xc)) ) /obj & + - ( dot_product(Xg-point_Xg, matmul(dTgc(:,2), this%Xc))*grad(2) ) / obj**2 + hess(3,2) = ( dot_product(matmul(dTgc(:,2),this%Xc), matmul(dTgc(:,3),this%Xc)) + dot_product((Xg-point_Xg), & + matmul(d2Tgc(2*this%nc(1)*this%nc(2)*this%nc(3)+1:3*this%nc(1)*this%nc(2)*this%nc(3) ,2),this%Xc)) ) /obj & + - ( dot_product(Xg-point_Xg, matmul(dTgc(:,3), this%Xc))*grad(2) ) / obj**2 + hess(1,3) = ( dot_product(matmul(dTgc(:,3),this%Xc), matmul(dTgc(:,1),this%Xc)) + dot_product((Xg-point_Xg), & + matmul(d2Tgc(1:this%nc(1)*this%nc(2)*this%nc(3) ,3),this%Xc)) ) /obj & + - ( dot_product(Xg-point_Xg, matmul(dTgc(:,1), this%Xc))*grad(3) ) / obj**2 + hess(2,3) = ( dot_product(matmul(dTgc(:,3),this%Xc), matmul(dTgc(:,2),this%Xc)) + dot_product((Xg-point_Xg), & + matmul(d2Tgc(this%nc(1)*this%nc(2)*this%nc(3)+1:2*this%nc(1)*this%nc(2)*this%nc(3) ,3),this%Xc)) ) /obj & + - ( dot_product(Xg-point_Xg, matmul(dTgc(:,2), this%Xc))*grad(3) ) / obj**2 + hess(3,3) = ( dot_product(matmul(dTgc(:,3),this%Xc), matmul(dTgc(:,3),this%Xc)) + dot_product((Xg-point_Xg), & + matmul(d2Tgc(2*this%nc(1)*this%nc(2)*this%nc(3)+1:3*this%nc(1)*this%nc(2)*this%nc(3) ,3),this%Xc)) ) /obj & + - ( dot_product(Xg-point_Xg, matmul(dTgc(:,3), this%Xc))*grad(3) ) / obj**2 + + ! debug + print '(i3,1x,3e20.10,1x,e20.10)', k, xk, norm2(grad) + + if (norm2(grad) <= tol) then + convergenz = .true. + nearest_Xt = xk + if (present(nearest_Xg)) nearest_Xg = this%cmp_Xg(nearest_Xt) + else + + ! Inverse of Hessian + det_inv = 1.0_rk/(& + + hess(1,1)*hess(2,2)*hess(3,3) - hess(1,1)*hess(2,3)*hess(3,2)& + - hess(1,2)*hess(2,1)*hess(3,3) + hess(1,2)*hess(2,3)*hess(3,1)& + + hess(1,3)*hess(2,1)*hess(3,2) - hess(1,3)*hess(2,2)*hess(3,1)) + Ainv(1,1) = +(hess(2,2)*hess(3,3) - hess(2,3)*hess(3,2)) + Ainv(2,1) = -(hess(2,1)*hess(3,3) - hess(2,3)*hess(3,1)) + Ainv(3,1) = +(hess(2,1)*hess(3,2) - hess(2,2)*hess(3,1)) + Ainv(1,2) = -(hess(1,2)*hess(3,3) - hess(1,3)*hess(3,2)) + Ainv(2,2) = +(hess(1,1)*hess(3,3) - hess(1,3)*hess(3,1)) + Ainv(3,2) = -(hess(1,1)*hess(3,2) - hess(1,2)*hess(3,1)) + Ainv(1,3) = +(hess(1,2)*hess(2,3) - hess(1,3)*hess(2,2)) + Ainv(2,3) = -(hess(1,1)*hess(2,3) - hess(1,3)*hess(2,1)) + Ainv(3,3) = +(hess(1,1)*hess(2,2) - hess(1,2)*hess(2,1)) + + Ainv = det_inv*Ainv + + dk = - matmul(Ainv, grad) + + ! Backtracking-Armijo Line Search + alphak = 1.0_rk + tau = 0.5_rk ! 0 < tau < 1 + beta = 1.0e-4_rk ! 0 < beta < 1 + l = 0 + do while (.not. norm2(this%cmp_Xg(xk + alphak*dk) - point_Xg) <= obj + alphak*beta*dot_product(grad,dk) .and. l<50) + alphak = tau * alphak + l = l + 1 + end do + + xk = xk + alphak*dk + ! Check if xk is within the knot vector range + xk = max(min(xk, upper_bounds), lower_bounds) + k = k + 1 + end if + end do + + end subroutine + !=============================================================================== + + !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause @@ -2840,6 +3251,35 @@ impure function compute_Xg_nurbs_3d(Xt, knot1, knot2, knot3, degree, nc, ng, Xc, !=============================================================================== +!=============================================================================== +!> author: Seyed Ali Ghasemi +!> license: BSD 3-Clause +impure function compute_Xg_nurbs_3d_1point(Xt, knot1, knot2, knot3, degree, nc, Xc, Wc) result(Xg) + use forcad_utils, only: rk, basis_bspline, kron + + implicit none + real(rk), intent(in), contiguous :: Xt(:) + real(rk), intent(in), contiguous :: knot1(:), knot2(:), knot3(:) + integer, intent(in) :: degree(3) + integer, intent(in) :: nc(3) + real(rk), intent(in), contiguous :: Xc(:,:) + real(rk), intent(in), contiguous :: Wc(:) + real(rk), allocatable :: Xg(:) + real(rk), allocatable :: Tgc(:) + + allocate(Xg(size(Xc,2))) + allocate(Tgc(nc(1)*nc(2)*nc(3))) + + Tgc = kron(basis_bspline(Xt(3), knot3, nc(3), degree(3)),& + kron(& + basis_bspline(Xt(2), knot2, nc(2), degree(2)),& + basis_bspline(Xt(1), knot1, nc(1), degree(1)))) + Tgc = Tgc*(Wc/(dot_product(Tgc,Wc))) + Xg = matmul(Tgc, Xc) +end function +!=============================================================================== + + !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause @@ -2872,7 +3312,30 @@ impure function compute_Xg_bspline_3d(Xt, knot1, knot2, knot3, degree, nc, ng, X !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause -impure function compute_dTgc_nurbs_3d(Xt, knot1, knot2, knot3, degree, nc, ng, Wc) result(dTgc) +impure function compute_Xg_bspline_3d_1point(Xt, knot1, knot2, knot3, degree, nc, Xc) result(Xg) + use forcad_utils, only: rk, basis_bspline, kron + + implicit none + real(rk), intent(in), contiguous :: Xt(:) + real(rk), intent(in), contiguous :: knot1(:), knot2(:), knot3(:) + integer, intent(in) :: degree(3) + integer, intent(in) :: nc(3) + real(rk), intent(in), contiguous :: Xc(:,:) + real(rk), allocatable :: Xg(:) + + allocate(Xg(size(Xc,2))) + Xg = matmul(kron(basis_bspline(Xt(3), knot3, nc(3), degree(3)), kron(& + basis_bspline(Xt(2), knot2, nc(2), degree(2)),& + basis_bspline(Xt(1), knot1, nc(1), degree(1)))),& + Xc) +end function +!=============================================================================== + + +!=============================================================================== +!> author: Seyed Ali Ghasemi +!> license: BSD 3-Clause +impure subroutine compute_dTgc_nurbs_3d_vector(Xt, knot1, knot2, knot3, degree, nc, ng, Wc, dTgc, Tgc) use forcad_utils, only: rk, basis_bspline_der, kron implicit none @@ -2882,28 +3345,78 @@ impure function compute_dTgc_nurbs_3d(Xt, knot1, knot2, knot3, degree, nc, ng, W integer, intent(in) :: nc(3) integer, intent(in) :: ng(3) real(rk), intent(in), contiguous :: Wc(:) - real(rk), allocatable :: dTgc(:,:) - real(rk), allocatable :: dTgci(:) + real(rk), allocatable, intent(out) :: dTgc(:,:,:) + real(rk), allocatable, intent(out) :: Tgc(:,:) + real(rk), allocatable :: dBi(:,:), dB1(:), dB2(:), dB3(:) + real(rk), allocatable :: Bi(:), B1(:), B2(:), B3(:) integer :: i - allocate(dTgc(ng(1)*ng(2)*ng(3), nc(1)*nc(2)*nc(3))) - allocate(dTgci(nc(1)*nc(2)*nc(3))) - !$OMP PARALLEL DO PRIVATE(dTgci) + allocate(dTgc(ng(1)*ng(2)*ng(3), nc(1)*nc(2)*nc(3), 3)) + allocate(Tgc(ng(1)*ng(2)*ng(3), nc(1)*nc(2)*nc(3))) + allocate(Bi(nc(1)*nc(2)*nc(3)), dBi(nc(1)*nc(2)*nc(3), 3)) do i = 1, size(Xt, 1) - dTgci = kron(basis_bspline_der(Xt(i,3), knot3, nc(3), degree(3)), kron(& - basis_bspline_der(Xt(i,2), knot2, nc(2), degree(2)),& - basis_bspline_der(Xt(i,1), knot1, nc(1), degree(1)))) - dTgc(i,:) = dTgci*(Wc/(dot_product(dTgci,Wc))) + call basis_bspline_der(Xt(i,1), knot1, nc(1), degree(1), dB1, B1) + call basis_bspline_der(Xt(i,2), knot2, nc(2), degree(2), dB2, B2) + call basis_bspline_der(Xt(i,3), knot3, nc(3), degree(3), dB3, B3) + + Bi = kron( B3, kron( B2, B1)) + Tgc(i,:) = Bi*(Wc/(dot_product(Bi,Wc))) + + dBi(:,1) = kron(kron(B3,B2),dB1) + dBi(:,2) = kron(kron(B3,dB2),B1) + dBi(:,3) = kron(kron(dB3,B2),B1) + + dTgc(i,:,1) = ( dBi(:,1)*Wc - Tgc(i,:)*dot_product(dBi(:,1),Wc) ) / dot_product(Bi,Wc) + dTgc(i,:,2) = ( dBi(:,2)*Wc - Tgc(i,:)*dot_product(dBi(:,2),Wc) ) / dot_product(Bi,Wc) + dTgc(i,:,3) = ( dBi(:,3)*Wc - Tgc(i,:)*dot_product(dBi(:,3),Wc) ) / dot_product(Bi,Wc) end do - !$OMP END PARALLEL DO -end function +end subroutine +!=============================================================================== + + +!=============================================================================== +!> author: Seyed Ali Ghasemi +!> license: BSD 3-Clause +impure subroutine compute_dTgc_nurbs_3d_scalar(Xt, knot1, knot2, knot3, degree, nc, Wc, dTgc, Tgc) + use forcad_utils, only: rk, basis_bspline_der, kron + + implicit none + real(rk), intent(in), contiguous :: Xt(:) + real(rk), intent(in), contiguous :: knot1(:), knot2(:), knot3(:) + integer, intent(in) :: degree(3) + integer, intent(in) :: nc(3) + real(rk), intent(in), contiguous :: Wc(:) + real(rk), allocatable, intent(out) :: dTgc(:,:) + real(rk), allocatable, intent(out) :: Tgc(:) + real(rk), allocatable :: dB1(:), dB2(:), dB3(:), dBi(:,:) + real(rk), allocatable :: B1(:), B2(:), B3(:), Bi(:) + + allocate(dTgc(nc(1)*nc(2)*nc(3), 3)) + allocate(Tgc(nc(1)*nc(2)*nc(3))) + allocate(dBi(nc(1)*nc(2)*nc(3), 3), Bi(nc(1)*nc(2)*nc(3))) + + call basis_bspline_der(Xt(1), knot1, nc(1), degree(1), dB1, B1) + call basis_bspline_der(Xt(2), knot2, nc(2), degree(2), dB2, B2) + call basis_bspline_der(Xt(3), knot3, nc(3), degree(3), dB3, B3) + + Bi = kron( B3, kron( B2, B1)) + Tgc = Bi*(Wc/(dot_product(Bi,Wc))) + + dBi(:,1) = kron(kron(B3,B2),dB1) + dBi(:,2) = kron(kron(B3,dB2),B1) + dBi(:,3) = kron(kron(dB3,B2),B1) + + dTgc(:,1) = ( dBi(:,1)*Wc - Tgc*dot_product(dBi(:,1),Wc) ) / dot_product(Bi,Wc) + dTgc(:,2) = ( dBi(:,2)*Wc - Tgc*dot_product(dBi(:,2),Wc) ) / dot_product(Bi,Wc) + dTgc(:,3) = ( dBi(:,3)*Wc - Tgc*dot_product(dBi(:,3),Wc) ) / dot_product(Bi,Wc) +end subroutine !=============================================================================== !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause -impure function compute_dTgc_bspline_3d(Xt, knot1, knot2, knot3, degree, nc, ng) result(dTgc) +impure subroutine compute_dTgc_bspline_3d_vector(Xt, knot1, knot2, knot3, degree, nc, ng, dTgc, Tgc) use forcad_utils, only: rk, basis_bspline_der, kron implicit none @@ -2912,25 +3425,332 @@ impure function compute_dTgc_bspline_3d(Xt, knot1, knot2, knot3, degree, nc, ng) integer, intent(in) :: degree(3) integer, intent(in) :: nc(3) integer, intent(in) :: ng(3) - real(rk), allocatable :: dTgc(:,:) + real(rk), allocatable, intent(out) :: dTgc(:,:,:) + real(rk), allocatable, intent(out) :: Tgc(:,:) + real(rk), allocatable :: dB1(:), dB2(:), dB3(:) + real(rk), allocatable :: B1(:), B2(:), B3(:) integer :: i - allocate(dTgc(ng(1)*ng(2)*ng(3), nc(1)*nc(2)*nc(3))) - !$OMP PARALLEL DO + allocate(dTgc(ng(1)*ng(2)*ng(3), nc(1)*nc(2)*nc(3), 3)) + allocate(Tgc(ng(1)*ng(2)*ng(3), nc(1)*nc(2)*nc(3))) + !$OMP PARALLEL DO PRIVATE(dB1, dB2, dB3) do i = 1, size(Xt, 1) - dTgc(i,:) = kron(basis_bspline_der(Xt(i,3), knot3, nc(3), degree(3)), kron(& - basis_bspline_der(Xt(i,2), knot2, nc(2), degree(2)),& - basis_bspline_der(Xt(i,1), knot1, nc(1), degree(1)))) + call basis_bspline_der(Xt(i,1), knot1, nc(1), degree(1), dB1, B1) + call basis_bspline_der(Xt(i,2), knot2, nc(2), degree(2), dB2, B2) + call basis_bspline_der(Xt(i,3), knot3, nc(3), degree(3), dB3, B3) + + Tgc(i,:) = kron(B3, kron(B2, B1)) + + dTgc(i,:,1) = kron(kron(B3,B2),dB1) + dTgc(i,:,2) = kron(kron(B3,dB2),B1) + dTgc(i,:,3) = kron(kron(dB3,B2),B1) end do !$OMP END PARALLEL DO -end function +end subroutine !=============================================================================== !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause -impure function compute_Tgc_nurbs_3d(Xt, knot1, knot2, knot3, degree, nc, ng, Wc) result(Tgc) +impure subroutine compute_dTgc_bspline_3d_scalar(Xt, knot1, knot2, knot3, degree, nc, dTgc, Tgc) + use forcad_utils, only: rk, basis_bspline_der, kron + + implicit none + real(rk), intent(in), contiguous :: Xt(:) + real(rk), intent(in), contiguous :: knot1(:), knot2(:), knot3(:) + integer, intent(in) :: degree(3) + integer, intent(in) :: nc(3) + real(rk), allocatable, intent(out) :: dTgc(:,:) + real(rk), allocatable, intent(out) :: Tgc(:) + real(rk), allocatable :: dB1(:), dB2(:), dB3(:) + real(rk), allocatable :: B1(:), B2(:), B3(:) + + allocate(dTgc(nc(1)*nc(2)*nc(3), 3)) + allocate(Tgc(nc(1)*nc(2)*nc(3))) + call basis_bspline_der(Xt(1), knot1, nc(1), degree(1), dB1, B1) + call basis_bspline_der(Xt(2), knot2, nc(2), degree(2), dB2, B2) + call basis_bspline_der(Xt(3), knot3, nc(3), degree(3), dB3, B3) + + Tgc = kron(B3, kron(B2, B1)) + dTgc(:,1) = kron(kron(B3,B2),dB1) + dTgc(:,2) = kron(kron(B3,dB2),B1) + dTgc(:,3) = kron(kron(dB3,B2),B1) +end subroutine +!=============================================================================== + + +!=============================================================================== +!> author: Seyed Ali Ghasemi +!> license: BSD 3-Clause +impure subroutine compute_d2Tgc_nurbs_3d_vector(Xt, knot1, knot2, knot3, degree, nc, ng, Wc, d2Tgc, dTgc, Tgc) + use forcad_utils, only: rk, basis_bspline_2der, kron + + implicit none + real(rk), intent(in), contiguous :: Xt(:,:) + real(rk), intent(in), contiguous :: knot1(:), knot2(:), knot3(:) + integer, intent(in) :: degree(3) + integer, intent(in) :: nc(3) + integer, intent(in) :: ng(3) + real(rk), intent(in), contiguous :: Wc(:) + real(rk), allocatable, intent(out) :: d2Tgc(:,:,:) + real(rk), allocatable, intent(out) :: dTgc(:,:,:) + real(rk), allocatable, intent(out) :: Tgc(:,:) + real(rk), allocatable :: d2Bi(:,:), d2B1(:), d2B2(:), d2B3(:) + real(rk), allocatable :: dBi(:,:), dB1(:), dB2(:), dB3(:) + real(rk), allocatable :: Bi(:), B1(:), B2(:), B3(:) + real(rk), allocatable :: Tgci(:), dTgci(:) + integer :: i + + allocate(B1(nc(1)), B2(nc(2)), B3(nc(3))) + allocate(dB1(nc(1)), dB2(nc(2)), dB3(nc(3))) + allocate(d2B1(nc(1)), d2B2(nc(2)), d2B3(nc(3))) + allocate(Bi(nc(1)*nc(2)*nc(3)), dBi(nc(1)*nc(2)*nc(3), 3), d2Bi(3*nc(1)*nc(2)*nc(3), 3)) + + allocate(Tgci(nc(1)*nc(2)*nc(3)), dTgci(nc(1)*nc(2)*nc(3))) + allocate(d2Tgc(ng(1)*ng(2)*ng(3), 3*nc(1)*nc(2)*nc(3), 3)) + allocate(dTgc(ng(1)*ng(2)*ng(3), nc(1)*nc(2)*nc(3), 3)) + allocate(Tgc(ng(1)*ng(2)*ng(3), nc(1)*nc(2)*nc(3))) + do i = 1, size(Xt, 1) + call basis_bspline_2der(Xt(i,1), knot1, nc(1), degree(1), d2B1, dB1, B1) + call basis_bspline_2der(Xt(i,2), knot2, nc(2), degree(2), d2B2, dB2, B2) + call basis_bspline_2der(Xt(i,3), knot3, nc(3), degree(3), d2B3, dB3, B3) + + Bi = kron(B3, kron(B2, B1)) + + Tgc(i,:) = Bi*(Wc/(dot_product(Bi,Wc))) + + dBi(:,1) = kron(kron(B3,B2),dB1) + dBi(:,2) = kron(kron(B3,dB2),B1) + dBi(:,3) = kron(kron(dB3,B2),B1) + + dTgc(i,:,1) = ( dBi(:,1)*Wc - Tgc(i,:)*dot_product(dBi(:,1),Wc) ) / dot_product(Bi,Wc) + dTgc(i,:,2) = ( dBi(:,2)*Wc - Tgc(i,:)*dot_product(dBi(:,2),Wc) ) / dot_product(Bi,Wc) + dTgc(i,:,3) = ( dBi(:,3)*Wc - Tgc(i,:)*dot_product(dBi(:,3),Wc) ) / dot_product(Bi,Wc) + + d2Bi(1:nc(1)*nc(2)*nc(3) ,1) = kron(kron(B3,B2),d2B1) + d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,1) = kron(kron(B3,dB2),dB1) + d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,1) = kron(kron(dB3,B2),dB1) + d2Bi(1:nc(1)*nc(2)*nc(3) ,2) = kron(kron(B3,dB2),dB1) + d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,2) = kron(kron(B3,d2B2),B1) + d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,2) = kron(kron(dB3,dB2),B1) + d2Bi(1:nc(1)*nc(2)*nc(3) ,3) = kron(kron(dB3,B2),dB1) + d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,3) = kron(kron(dB3,dB2),B1) + d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,3) = kron(kron(d2B3,B2),B1) + + d2Tgc(i, 1:nc(1)*nc(2)*nc(3) ,1) = (d2Bi(1:nc(1)*nc(2)*nc(3) ,1)*Wc & + - 2.0_rk*dTgc(i, :,1)*dot_product(dBi(:,1),Wc) & + - Tgc(i,:)*dot_product(d2Bi(1:nc(1)*nc(2)*nc(3) ,1),Wc)) / dot_product(Bi,Wc) + d2Tgc(i, nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,1) = (d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,1)*Wc & + - dTgc(i, :,1)*dot_product(dBi(:,2),Wc) - dTgc(i, :,2)*dot_product(dBi(:,1),Wc) & + - Tgc(i,:)*dot_product(d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,1),Wc)) / dot_product(Bi,Wc) + d2Tgc(i, 2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,1) = (d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,1)*Wc & + - dTgc(i, :,1)*dot_product(dBi(:,3),Wc) - dTgc(i, :,3)*dot_product(dBi(:,1),Wc) & + - Tgc(i,:)*dot_product(d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,1),Wc)) / dot_product(Bi,Wc) + d2Tgc(i, 1:nc(1)*nc(2)*nc(3) ,2) = (d2Bi(1:nc(1)*nc(2)*nc(3) ,2)*Wc & + - dTgc(i, :,1)*dot_product(dBi(:,2),Wc) - dTgc(i, :,2)*dot_product(dBi(:,1),Wc) & + - Tgc(i,:)*dot_product(d2Bi(1:nc(1)*nc(2)*nc(3) ,2),Wc)) / dot_product(Bi,Wc) + d2Tgc(i, nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,2) = (d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,2)*Wc & + - 2.0_rk*dTgc(i, :,2)*dot_product(dBi(:,2),Wc) & + - Tgc(i,:)*dot_product(d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,2),Wc)) / dot_product(Bi,Wc) + d2Tgc(i, 2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,2) = (d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,2)*Wc & + - dTgc(i, :,2)*dot_product(dBi(:,3),Wc) - dTgc(i, :,3)*dot_product(dBi(:,2),Wc) & + - Tgc(i,:)*dot_product(d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,2),Wc)) / dot_product(Bi,Wc) + d2Tgc(i, 1:nc(1)*nc(2)*nc(3) ,3) = (d2Bi(1:nc(1)*nc(2)*nc(3) ,3)*Wc & + - dTgc(i, :,1)*dot_product(dBi(:,3),Wc) - dTgc(i, :,3)*dot_product(dBi(:,1),Wc) & + - Tgc(i,:)*dot_product(d2Bi(1:nc(1)*nc(2)*nc(3) ,3),Wc)) / dot_product(Bi,Wc) + d2Tgc(i, nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,3) = (d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,3)*Wc & + - dTgc(i, :,2)*dot_product(dBi(:,3),Wc) - dTgc(i, :,3)*dot_product(dBi(:,2),Wc) & + - Tgc(i,:)*dot_product(d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,3),Wc)) / dot_product(Bi,Wc) + d2Tgc(i, 2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,3) = (d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,3)*Wc & + - 2.0_rk*dTgc(i, :,3)*dot_product(dBi(:,3),Wc) & + - Tgc(i,:)*dot_product(d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,3),Wc)) / dot_product(Bi,Wc) + end do +end subroutine +!=============================================================================== + + +!=============================================================================== +!> author: Seyed Ali Ghasemi +!> license: BSD 3-Clause +impure subroutine compute_d2Tgc_nurbs_3d_scalar(Xt, knot1, knot2, knot3, degree, nc, Wc, d2Tgc, dTgc, Tgc) + use forcad_utils, only: rk, basis_bspline_2der, kron + + implicit none + real(rk), intent(in), contiguous :: Xt(:) + real(rk), intent(in), contiguous :: knot1(:), knot2(:), knot3(:) + integer, intent(in) :: degree(3) + integer, intent(in) :: nc(3) + real(rk), intent(in), contiguous :: Wc(:) + real(rk), allocatable, intent(out) :: d2Tgc(:,:) + real(rk), allocatable, intent(out) :: dTgc(:,:) + real(rk), allocatable, intent(out) :: Tgc(:) + real(rk), allocatable :: d2Bi(:,:), d2B1(:), d2B2(:), d2B3(:) + real(rk), allocatable :: dBi(:,:), dB1(:), dB2(:), dB3(:) + real(rk), allocatable :: Bi(:), B1(:), B2(:), B3(:) + + allocate(B1(nc(1)), B2(nc(2)), B3(nc(3))) + allocate(dB1(nc(1)), dB2(nc(2)), dB3(nc(3))) + allocate(d2B1(nc(1)), d2B2(nc(2)), d2B3(nc(3))) + allocate(Bi(nc(1)*nc(2)*nc(3)), dBi(nc(1)*nc(2)*nc(3), 3), d2Bi(3*nc(1)*nc(2)*nc(3), 3)) + + allocate(d2Tgc(3*nc(1)*nc(2)*nc(3), 3)) + allocate(dTgc(nc(1)*nc(2)*nc(3), 3)) + allocate(Tgc(nc(1)*nc(2)*nc(3))) + + call basis_bspline_2der(Xt(1), knot1, nc(1), degree(1), d2B1, dB1, B1) + call basis_bspline_2der(Xt(2), knot2, nc(2), degree(2), d2B2, dB2, B2) + call basis_bspline_2der(Xt(3), knot3, nc(3), degree(3), d2B3, dB3, B3) + + Bi = kron(B3, kron(B2, B1)) + + Tgc = Bi*(Wc/(dot_product(Bi,Wc))) + + dBi(:,1) = kron(kron(B3,B2),dB1) + dBi(:,2) = kron(kron(B3,dB2),B1) + dBi(:,3) = kron(kron(dB3,B2),B1) + + dTgc(:,1) = ( dBi(:,1)*Wc - Tgc*dot_product(dBi(:,1),Wc) ) / dot_product(Bi,Wc) + dTgc(:,2) = ( dBi(:,2)*Wc - Tgc*dot_product(dBi(:,2),Wc) ) / dot_product(Bi,Wc) + dTgc(:,3) = ( dBi(:,3)*Wc - Tgc*dot_product(dBi(:,3),Wc) ) / dot_product(Bi,Wc) + + d2Bi(1:nc(1)*nc(2)*nc(3) ,1) = kron(kron(B3,B2),d2B1) + d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,1) = kron(kron(B3,dB2),dB1) + d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,1) = kron(kron(dB3,B2),dB1) + d2Bi(1:nc(1)*nc(2)*nc(3) ,2) = kron(kron(B3,dB2),dB1) + d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,2) = kron(kron(B3,d2B2),B1) + d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,2) = kron(kron(dB3,dB2),B1) + d2Bi(1:nc(1)*nc(2)*nc(3) ,3) = kron(kron(dB3,B2),dB1) + d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,3) = kron(kron(dB3,dB2),B1) + d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,3) = kron(kron(d2B3,B2),B1) + + d2Tgc(1:nc(1)*nc(2)*nc(3) ,1) = (d2Bi(1:nc(1)*nc(2)*nc(3) ,1)*Wc & + - 2.0_rk*dTgc(:,1)*dot_product(dBi(:,1),Wc) & + - Tgc*dot_product(d2Bi(1:nc(1)*nc(2)*nc(3) ,1),Wc)) / dot_product(Bi,Wc) + d2Tgc(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,1) = (d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,1)*Wc & + - dTgc(:,1)*dot_product(dBi(:,2),Wc) - dTgc(:,2)*dot_product(dBi(:,1),Wc) & + - Tgc*dot_product(d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,1),Wc)) / dot_product(Bi,Wc) + d2Tgc(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,1) = (d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,1)*Wc & + - dTgc(:,1)*dot_product(dBi(:,3),Wc) - dTgc(:,3)*dot_product(dBi(:,1),Wc) & + - Tgc*dot_product(d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,1),Wc)) / dot_product(Bi,Wc) + d2Tgc(1:nc(1)*nc(2)*nc(3) ,2) = (d2Bi(1:nc(1)*nc(2)*nc(3) ,2)*Wc & + - dTgc(:,1)*dot_product(dBi(:,2),Wc) - dTgc(:,2)*dot_product(dBi(:,1),Wc) & + - Tgc*dot_product(d2Bi(1:nc(1)*nc(2)*nc(3) ,2),Wc)) / dot_product(Bi,Wc) + d2Tgc(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,2) = (d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,2)*Wc & + - 2.0_rk*dTgc(:,2)*dot_product(dBi(:,2),Wc) & + - Tgc*dot_product(d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,2),Wc)) / dot_product(Bi,Wc) + d2Tgc(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,2) = (d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,2)*Wc & + - dTgc(:,2)*dot_product(dBi(:,3),Wc) - dTgc(:,3)*dot_product(dBi(:,2),Wc) & + - Tgc*dot_product(d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,2),Wc)) / dot_product(Bi,Wc) + d2Tgc(1:nc(1)*nc(2)*nc(3) ,3) = (d2Bi(1:nc(1)*nc(2)*nc(3) ,3)*Wc & + - dTgc(:,1)*dot_product(dBi(:,3),Wc) - dTgc(:,3)*dot_product(dBi(:,1),Wc) & + - Tgc*dot_product(d2Bi(1:nc(1)*nc(2)*nc(3) ,3),Wc)) / dot_product(Bi,Wc) + d2Tgc(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,3) = (d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,3)*Wc & + - dTgc(:,2)*dot_product(dBi(:,3),Wc) - dTgc(:,3)*dot_product(dBi(:,2),Wc) & + - Tgc*dot_product(d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,3),Wc)) / dot_product(Bi,Wc) + d2Tgc(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,3) = (d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,3)*Wc & + - 2.0_rk*dTgc(:,3)*dot_product(dBi(:,3),Wc) & + - Tgc*dot_product(d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,3),Wc)) / dot_product(Bi,Wc) +end subroutine +!=============================================================================== + + +!=============================================================================== +!> author: Seyed Ali Ghasemi +!> license: BSD 3-Clause +impure subroutine compute_d2Tgc_bspline_3d_vector(Xt, knot1, knot2, knot3, degree, nc, ng, d2Tgc, dTgc, Tgc) + use forcad_utils, only: rk, basis_bspline_2der, kron + + implicit none + real(rk), intent(in), contiguous :: Xt(:,:) + real(rk), intent(in), contiguous :: knot1(:), knot2(:), knot3(:) + integer, intent(in) :: degree(3) + integer, intent(in) :: nc(3) + integer, intent(in) :: ng(3) + real(rk), allocatable, intent(out) :: d2Tgc(:,:,:) + real(rk), allocatable, intent(out) :: dTgc(:,:,:) + real(rk), allocatable, intent(out) :: Tgc(:,:) + real(rk), allocatable :: d2B1(:), d2B2(:), d2B3(:) + real(rk), allocatable :: dB1(:), dB2(:), dB3(:) + real(rk), allocatable :: B1(:), B2(:), B3(:) + integer :: i + + allocate(d2Tgc(ng(1)*ng(2)*ng(3), 3*nc(1)*nc(2)*nc(3), 3)) + allocate(dTgc(ng(1)*ng(2)*ng(3), nc(1)*nc(2)*nc(3), 3)) + allocate(Tgc(ng(1)*ng(2)*ng(3), nc(1)*nc(2)*nc(3))) + do i = 1, size(Xt, 1) + call basis_bspline_2der(Xt(i,1), knot1, nc(1), degree(1), d2B1, dB1, B1) + call basis_bspline_2der(Xt(i,2), knot2, nc(2), degree(2), d2B2, dB2, B2) + call basis_bspline_2der(Xt(i,3), knot3, nc(3), degree(3), d2B3, dB3, B3) + + Tgc(i,:) = kron(B3, kron(B2, B1)) + + dTgc(i,:,1) = kron(kron(B3,B2),dB1) + dTgc(i,:,2) = kron(kron(B3,dB2),B1) + dTgc(i,:,3) = kron(kron(dB3,B2),B1) + + d2Tgc(i,1:nc(1)*nc(2)*nc(3) ,1) = kron(kron(B3,B2),d2B1) + d2Tgc(i,nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,1) = kron(kron(B3,dB2),dB1) + d2Tgc(i,2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,1) = kron(kron(dB3,B2),dB1) + d2Tgc(i,1:nc(1)*nc(2)*nc(3) ,2) = kron(kron(B3,dB2),dB1) + d2Tgc(i,nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,2) = kron(kron(B3,d2B2),B1) + d2Tgc(i,2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,2) = kron(kron(dB3,dB2),B1) + d2Tgc(i,1:nc(1)*nc(2)*nc(3) ,3) = kron(kron(dB3,B2),dB1) + d2Tgc(i,nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,3) = kron(kron(dB3,dB2),B1) + d2Tgc(i,2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,3) = kron(kron(d2B3,B2),B1) + end do +end subroutine +!=============================================================================== + + +!=============================================================================== +!> author: Seyed Ali Ghasemi +!> license: BSD 3-Clause +impure subroutine compute_d2Tgc_bspline_3d_scalar(Xt, knot1, knot2, knot3, degree, nc, d2Tgc, dTgc, Tgc) + use forcad_utils, only: rk, basis_bspline_2der, kron + + implicit none + real(rk), intent(in), contiguous :: Xt(:) + real(rk), intent(in), contiguous :: knot1(:), knot2(:), knot3(:) + integer, intent(in) :: degree(3) + integer, intent(in) :: nc(3) + real(rk), allocatable, intent(out) :: d2Tgc(:,:) + real(rk), allocatable, intent(out) :: dTgc(:,:) + real(rk), allocatable, intent(out) :: Tgc(:) + real(rk), allocatable :: d2B1(:), d2B2(:), d2B3(:) + real(rk), allocatable :: dB1(:), dB2(:), dB3(:) + real(rk), allocatable :: B1(:), B2(:), B3(:) + + allocate(d2Tgc(3*nc(1)*nc(2)*nc(3), 3)) + allocate(dTgc(nc(1)*nc(2)*nc(3), 3)) + allocate(Tgc(nc(1)*nc(2)*nc(3))) + call basis_bspline_2der(Xt(1), knot1, nc(1), degree(1), d2B1, dB1, B1) + call basis_bspline_2der(Xt(2), knot2, nc(2), degree(2), d2B2, dB2, B2) + call basis_bspline_2der(Xt(3), knot3, nc(3), degree(3), d2B3, dB3, B3) + + Tgc = kron(B3, kron(B2, B1)) + + dTgc(:,1) = kron(kron(B3,B2),dB1) + dTgc(:,2) = kron(kron(B3,dB2),B1) + dTgc(:,3) = kron(kron(dB3,B2),B1) + + d2Tgc(1:nc(1)*nc(2)*nc(3) ,1) = kron(kron(B3,B2),d2B1) + d2Tgc(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,1) = kron(kron(B3,dB2),dB1) + d2Tgc(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,1) = kron(kron(dB3,B2),dB1) + d2Tgc(1:nc(1)*nc(2)*nc(3) ,2) = kron(kron(B3,dB2),dB1) + d2Tgc(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,2) = kron(kron(B3,d2B2),B1) + d2Tgc(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,2) = kron(kron(dB3,dB2),B1) + d2Tgc(1:nc(1)*nc(2)*nc(3) ,3) = kron(kron(dB3,B2),dB1) + d2Tgc(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,3) = kron(kron(dB3,dB2),B1) + d2Tgc(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,3) = kron(kron(d2B3,B2),B1) +end subroutine +!=============================================================================== + + +!=============================================================================== +!> author: Seyed Ali Ghasemi +!> license: BSD 3-Clause +impure function compute_Tgc_nurbs_3d_vector(Xt, knot1, knot2, knot3, degree, nc, ng, Wc) result(Tgc) use forcad_utils, only: rk, basis_bspline, kron implicit none @@ -2961,7 +3781,30 @@ impure function compute_Tgc_nurbs_3d(Xt, knot1, knot2, knot3, degree, nc, ng, Wc !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause -impure function compute_Tgc_bspline_3d(Xt, knot1, knot2, knot3, degree, nc, ng) result(Tgc) +impure function compute_Tgc_nurbs_3d_scalar(Xt, knot1, knot2, knot3, degree, nc, Wc) result(Tgc) + use forcad_utils, only: rk, basis_bspline, kron + + implicit none + real(rk), intent(in), contiguous :: Xt(:) + real(rk), intent(in), contiguous :: knot1(:), knot2(:), knot3(:) + integer, intent(in) :: degree(3) + integer, intent(in) :: nc(3) + real(rk), intent(in), contiguous :: Wc(:) + real(rk), allocatable :: Tgc(:) + + allocate(Tgc(nc(1)*nc(2)*nc(3))) + Tgc = kron(basis_bspline(Xt(3), knot3, nc(3), degree(3)), kron(& + basis_bspline(Xt(2), knot2, nc(2), degree(2)),& + basis_bspline(Xt(1), knot1, nc(1), degree(1)))) + Tgc = Tgc*(Wc/(dot_product(Tgc,Wc))) +end function +!=============================================================================== + + +!=============================================================================== +!> author: Seyed Ali Ghasemi +!> license: BSD 3-Clause +impure function compute_Tgc_bspline_3d_vector(Xt, knot1, knot2, knot3, degree, nc, ng) result(Tgc) use forcad_utils, only: rk, basis_bspline, kron implicit none @@ -2985,6 +3828,27 @@ impure function compute_Tgc_bspline_3d(Xt, knot1, knot2, knot3, degree, nc, ng) !=============================================================================== +!=============================================================================== +!> author: Seyed Ali Ghasemi +!> license: BSD 3-Clause +impure function compute_Tgc_bspline_3d_scalar(Xt, knot1, knot2, knot3, degree, nc) result(Tgc) + use forcad_utils, only: rk, basis_bspline, kron + + implicit none + real(rk), intent(in), contiguous :: Xt(:) + real(rk), intent(in), contiguous :: knot1(:), knot2(:), knot3(:) + integer, intent(in) :: degree(3) + integer, intent(in) :: nc(3) + real(rk), allocatable :: Tgc(:) + + allocate(Tgc(nc(1)*nc(2)*nc(3))) + Tgc= kron(basis_bspline(Xt(3), knot3, nc(3), degree(3)), kron(& + basis_bspline(Xt(2), knot2, nc(2), degree(2)),& + basis_bspline(Xt(1), knot1, nc(1), degree(1)))) +end function +!=============================================================================== + + !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause @@ -3006,4 +3870,4 @@ impure function nearest_point_help_3d(ng, Xg, point_Xg) result(distances) !$OMP END PARALLEL DO end function -!=============================================================================== \ No newline at end of file +!=============================================================================== diff --git a/src/forcad_utils.f90 b/src/forcad_utils.f90 index ac2bc25ff..95cac9104 100644 --- a/src/forcad_utils.f90 +++ b/src/forcad_utils.f90 @@ -8,7 +8,7 @@ module forcad_utils private public :: rk, basis_bernstein, basis_bspline, elemConn_C0, kron, ndgrid, compute_multiplicity, compute_knot_vector, & basis_bspline_der, insert_knot_A_5_1, findspan, elevate_degree_A_5_9, hexahedron_Xc, tetragon_Xc, remove_knots_A_5_8, & - elemConn_Cn, unique, rotation + elemConn_Cn, unique, rotation, basis_bspline_2der integer, parameter :: rk = kind(1.0d0) @@ -91,58 +91,102 @@ pure function basis_bspline(Xt, knot, nc, degree) result(B) !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause - pure function basis_bspline_der(Xt, knot, nc, degree) result(dB) + pure subroutine basis_bspline_der(Xt, knot, nc, degree, dB, B) integer, intent(in) :: degree real(rk), intent(in), contiguous :: knot(:) integer, intent(in) :: nc real(rk), intent(in) :: Xt - real(rk), allocatable :: dB(:) - real(rk), allocatable :: Nt(:,:), dNt_dXt(:,:) - real(rk) :: R, L, Rp, Lp, knot_i, knot_ip, knot_jk, knot_jkm, knot_end, a, b, c, d - integer :: i, k, n, m, jk + real(rk), allocatable, intent(out) :: dB(:) + real(rk), allocatable, intent(out), optional :: B(:) + real(rk), allocatable :: N(:,:), dN_dXt(:,:) + real(rk) :: temp, Xth_i, Xth_i1, Xth_ip, Xth_ip1 + integer :: i, p - k = degree + 1 - n = nc - 1 - allocate(Nt(nc+degree, degree+1)) - Nt = 0.0_rk - do i = 1, n+k - knot_i = knot(i) - knot_ip = knot(i+1) - knot_end = knot(size(knot)) - if ( abs(Xt - knot_end) > tiny(0.0_rk) ) then - if ( Xt >= knot_i .and. Xt < knot_ip ) Nt(i,1) = 1.0_rk - elseif ( abs(Xt - knot_end) < tiny(0.0_rk) ) then - if ( Xt >= knot_i .and. Xt <= knot_ip ) Nt(i,1) = 1.0_rk - end if + temp = abs(Xt - knot(size(knot))) + allocate(N(nc, 0:degree), source=0.0_rk) + allocate(dN_dXt(nc, 0:degree), source=0.0_rk) + + do p = 0, degree + do concurrent (i = 1:nc) + Xth_i = knot(i) + Xth_i1 = knot(i+1) + Xth_ip = knot(i+p) + Xth_ip1 = knot(i+p+1) + + if ( temp /= tiny(0.0_rk) .and. Xth_i <= Xt .and. Xt <= Xth_i1 ) then + N(i,0) = 1.0_rk + dN_dXt(i,0) = 0.0_rk + end if + if ( Xth_ip /= Xth_i ) then + N(i,p) = (Xt - Xth_i)/(Xth_ip - Xth_i) * N(i,p-1) + dN_dXt(i,p) = ( N(i,p-1) + (Xt - Xth_i)*dN_dXt(i,p-1) ) / (Xth_ip - Xth_i) + end if + if ( Xth_ip1 /= Xth_i1 ) then + N(i,p) = N(i,p) + (Xth_ip1 - Xt)/(Xth_ip1 - Xth_i1) * N(i+1,p-1) + dN_dXt(i,p) = dN_dXt(i,p) - ( N(i+1,p-1) - (Xth_ip1 - Xt)*dN_dXt(i+1,p-1) ) / (Xth_ip1 - Xth_i1) + end if + + end do end do - allocate(dNt_dXt(nc+degree, degree+1)) - dNt_dXt = 0.0_rk - m = 0 - do jk = 2, k - m = m + 1 - do i = 1, n + k - m - knot_i = knot(i) - knot_ip = knot(i+1) - knot_jk = knot(i+jk) - knot_jkm = knot(i+jk-1) - a = (knot_jkm - knot_i) - b = (knot_jk - Xt) - c = (knot_jk - knot_ip) - d = (Xt - knot_i) - R = d/a - if ( isnan(R) .or. isinf(R) .or. abs(R) < tiny(0.0_rk) ) R = 0.0_rk - L = b/c - if ( isnan(L) .or. isinf(L) .or. abs(L) < tiny(0.0_rk) ) L = 0.0_rk - Nt(i,jk) = R*Nt(i,jk-1) + L*Nt(i+1,jk-1) - Rp = (Nt(i,jk-1) + d*dNt_dXt(i,jk-1)) / a - if ( isnan(Rp) .or. isinf(Rp) .or. abs(Rp) < tiny(0.0_rk) ) Rp = 0.0_rk - Lp = (b*dNt_dXt(i+1,jk-1) - Nt(i+1,jk-1)) / c - if ( isnan(Lp) .or. isinf(Lp) .or. abs(Lp) < tiny(0.0_rk) ) Lp = 0.0_rk - dNt_dXt(i,jk) = Rp + Lp + + dB = dN_dXt(:,degree) + if (present(B)) B = N(:,degree) + end subroutine + !=============================================================================== + + + !=============================================================================== + !> author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + pure subroutine basis_bspline_2der(Xt, knot, nc, degree, d2B, dB, B) + integer, intent(in) :: degree + real(rk), intent(in), contiguous :: knot(:) + integer, intent(in) :: nc + real(rk), intent(in) :: Xt + real(rk), allocatable, intent(out) :: d2B(:) + real(rk), allocatable, intent(out), optional :: dB(:) + real(rk), allocatable, intent(out), optional :: B(:) + real(rk), allocatable :: N(:,:), dN_dXt(:,:), d2N_dXt2(:,:) + real(rk) :: temp, Xth_i, Xth_i1, Xth_ip, Xth_ip1 + integer :: i, p + + temp = abs(Xt - knot(size(knot))) + allocate(N(nc, 0:degree), source=0.0_rk) + allocate(dN_dXt(nc, 0:degree), source=0.0_rk) + allocate(d2N_dXt2(nc, 0:degree), source=0.0_rk) + + do p = 0, degree + do concurrent (i = 1:nc) + Xth_i = knot(i) + Xth_i1 = knot(i+1) + Xth_ip = knot(i+p) + Xth_ip1 = knot(i+p+1) + + if ( temp /= tiny(0.0_rk) .and. Xth_i <= Xt .and. Xt <= Xth_i1 ) then + N(i,0) = 1.0_rk + dN_dXt(i,0) = 0.0_rk + d2N_dXt2(i,0) = 0.0_rk + end if + if ( Xth_ip /= Xth_i) then + N(i,p) = (Xt - Xth_i)/(Xth_ip - Xth_i) * N(i,p-1) + ! dN_dXt(i,p) = p*(N(i,p-1)/(Xth_ip - Xth_i)) + dN_dXt(i,p) = ( N(i,p-1) + (Xt - Xth_i)*dN_dXt(i,p-1) ) / (Xth_ip - Xth_i) + d2N_dXt2(i,p) = (2*dN_dXt(i,p-1) + (Xt - Xth_i)*d2N_dXt2(i,p-1)) / (Xth_ip - Xth_i) + end if + if ( Xth_ip1 /= Xth_i1 ) then + N(i,p) = N(i,p) + (Xth_ip1 - Xt)/(Xth_ip1 - Xth_i1) * N(i+1,p-1) + ! dN_dXt(i,p) = dN_dXt(i,p) - p*( N(i+1,p-1)/(Xth_ip1 - Xth_i1)) + dN_dXt(i,p) = dN_dXt(i,p) - ( N(i+1,p-1) - (Xth_ip1 - Xt)*dN_dXt(i+1,p-1) ) / (Xth_ip1 - Xth_i1) + d2N_dXt2(i,p) = d2N_dXt2(i,p) - (2*dN_dXt(i+1,p-1) - (Xth_ip1 - Xt)*d2N_dXt2(i+1,p-1)) / (Xth_ip1 - Xth_i1) + end if + end do end do - dB = dNt_dXt(1:nc,k) - end function + + d2B = d2N_dXt2(:,degree) + if (present(dB)) dB = dN_dXt(:,degree) + if (present(B)) B = N(:,degree) + end subroutine !=============================================================================== diff --git a/test/fdm_curve.f90 b/test/fdm_curve.f90 new file mode 100644 index 000000000..325d8b0a8 --- /dev/null +++ b/test/fdm_curve.f90 @@ -0,0 +1,119 @@ +program fdm_test_curve + + use forcad, only: rk, nurbs_curve + + implicit none + + type(nurbs_curve) :: curve !! 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) :: Xtp, tol, Xt, Xtm + real(rk), allocatable :: Tgc(:), dTgc(:), Tgcp(:), dTgcp(:), Tgcm(:), dTgcm(:), & + CFD(:), BFD(:), FFD(:), d2Tgc(:), d2Tgcp(:), d2Tgcm(:) + + !----------------------------------------------------------------------------- + ! Setting up the NURBS curve + !----------------------------------------------------------------------------- + + !> 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 curve%set(knot, Xc, Wc) + + !----------------------------------------------------------------------------- + ! Creating the NURBS curve + !----------------------------------------------------------------------------- + + !> Generate the NURBS curve with a resolution of 20 + call curve%create(res = 20) + + !----------------------------------------------------------------------------- + ! Finite Difference Derivative Test + !----------------------------------------------------------------------------- + + tol = 1.0e-6_rk + Xt = 0.5_rk + call curve%derivative2(Xt=Xt, d2Tgc=d2Tgc, dTgc=dTgc, Tgc=Tgc) + Xtm = Xt - tol + call curve%derivative2(Xt=Xtm, d2Tgc=d2Tgcm, dTgc=dTgcm, Tgc=Tgcm) + Xtp = Xt + tol + call curve%derivative2(Xt=Xtp, d2Tgc=d2Tgcp, dTgc=dTgcp, Tgc=Tgcp) + print *, 'Tolerance:', tol + print *, 'Error BFD dTgc: ', norm2((Tgc - Tgcm)/tol - dTgc) + print *, 'Error CFD dTgc: ', norm2((Tgcp - Tgcm)/(2.0_rk*tol) - dTgc) + print *, 'Error FFD dTgc: ', norm2((Tgcp - Tgc )/tol - dTgc) + print *, 'Error BFD d2Tgc: ', norm2((dTgc - dTgcm)/tol - d2Tgc) + print *, 'Error CFD d2Tgc: ', norm2((dTgcp - dTgcm)/(2.0_rk*tol) - d2Tgc) + print *, 'Error FFD d2Tgc: ', norm2((dTgcp - dTgc )/tol - d2Tgc) + + !> Finalize the NURBS curve object + call curve%finalize() + deallocate(Xc, Wc) + + + + + + + + + + !----------------------------------------------------------------------------- + ! Setting up the NURBS curve + !----------------------------------------------------------------------------- + + !> 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 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 for the NURBS curve object. + !> Wc is optional + call curve%set(knot, Xc) + + !----------------------------------------------------------------------------- + ! Creating the NURBS curve + !----------------------------------------------------------------------------- + + !> Generate the NURBS curve with a resolution of 20 + call curve%create(res = 20) + + !----------------------------------------------------------------------------- + ! Finite Difference Derivative Test + !----------------------------------------------------------------------------- + + tol = 1.0e-6_rk + Xt = 0.5_rk + call curve%derivative2(Xt=Xt, d2Tgc=d2Tgc, dTgc=dTgc, Tgc=Tgc) + Xtm = Xt - tol + call curve%derivative2(Xt=Xtm, d2Tgc=d2Tgcm, dTgc=dTgcm, Tgc=Tgcm) + Xtp = Xt + tol + call curve%derivative2(Xt=Xtp, d2Tgc=d2Tgcp, dTgc=dTgcp, Tgc=Tgcp) + print *, 'Tolerance:', tol + print *, 'Error BFD dTgc: ', norm2((Tgc - Tgcm)/tol - dTgc) + print *, 'Error CFD dTgc: ', norm2((Tgcp - Tgcm)/(2.0_rk*tol) - dTgc) + print *, 'Error FFD dTgc: ', norm2((Tgcp - Tgc )/tol - dTgc) + print *, 'Error BFD d2Tgc: ', norm2((dTgc - dTgcm)/tol - d2Tgc) + print *, 'Error CFD d2Tgc: ', norm2((dTgcp - dTgcm)/(2.0_rk*tol) - d2Tgc) + print *, 'Error FFD d2Tgc: ', norm2((dTgcp - dTgc )/tol - d2Tgc) + + !> Finalize the NURBS curve object + call curve%finalize() + deallocate(Xc) +end program diff --git a/test/fdm_surface.f90 b/test/fdm_surface.f90 new file mode 100644 index 000000000..593328bf1 --- /dev/null +++ b/test/fdm_surface.f90 @@ -0,0 +1,144 @@ +program fdm_test_surface + + use forcad, only: rk, nurbs_surface + + implicit none + + type(nurbs_surface) :: surface !! Declare a NURBS surface object + real(rk), allocatable :: Xc(:,:), Wc(:) !! Arrays for control points and weights + real(rk) :: Xtp(2), tol, Xt(2), Xtm(2) + real(rk), allocatable :: Tgc(:), dTgc(:,:), Tgcp(:), dTgcp(:,:), Tgcm(:), dTgcm(:,:), d2Tgc(:,:), d2Tgcp(:,:), d2Tgcm(:,:) + real(rk), allocatable :: CFD(:,:), BFD(:,:), FFD(:,:), CFD2(:,:), BFD2(:,:), FFD2(:,:) + integer :: i + + !----------------------------------------------------------------------------- + ! Setting up the NURBS surface + !----------------------------------------------------------------------------- + Wc = [1.0_rk, 0.5_rk, 1.0_rk, 0.2_rk, 1.0_rk, 0.5_rk, 1.0_rk,& + 0.2_rk, 1.0_rk, 0.5_rk, 1.0_rk, 0.2_rk, 1.0_rk, 0.5_rk, 1.0_rk, 0.2_rk] + call surface%set_tetragon(L=[5.0_rk, 8.0_rk], nc=[4,4], Wc=Wc) + + !----------------------------------------------------------------------------- + ! Creating the NURBS surface + !----------------------------------------------------------------------------- + + !> Generate the NURBS surface with a resolution of 20 + call surface%create(20, 20) + + !----------------------------------------------------------------------------- + ! Finite Difference Derivative Test + !----------------------------------------------------------------------------- + + allocate(CFD(16,2), BFD(16,2), FFD(16,2)) + allocate(CFD2(2*16,2), BFD2(2*16,2), FFD2(2*16,2)) + + tol = 1.0e-6_rk + + Xt(1) = 0.5_rk + Xt(2) = 0.3_rk + call surface%derivative2(Xt=Xt, d2Tgc=d2Tgc, dTgc=dTgc, Tgc=Tgc) + + do i = 1, 2 + Xtm = Xt + Xtm(i) = Xt(i) - tol + call surface%derivative2(Xt=Xtm, d2Tgc=d2Tgcm, dTgc=dTgcm, Tgc=Tgcm) + + Xtp = Xt + Xtp(i) = Xt(i) + tol + call surface%derivative2(Xt=Xtp, d2Tgc=d2Tgcp, dTgc=dTgcp, Tgc=Tgcp) + + BFD(:,i) = (Tgc - Tgcm)/tol + CFD(:,i) = (Tgcp - Tgcm)/(2.0_rk*tol) + FFD(:,i) = (Tgcp - Tgc )/tol + + BFD2(:,i) = reshape((dTgc - dTgcm)/tol, shape=[2*16]) + CFD2(:,i) = reshape((dTgcp - dTgcm)/(2.0_rk*tol), shape=[2*16]) + FFD2(:,i) = reshape((dTgcp - dTgc )/tol, shape=[2*16]) + end do + + print *, 'Tolerance:', tol + print *, 'Error BFD dTgc: ', norm2(BFD - dTgc) + print *, 'Error CFD dTgc: ', norm2(CFD - dTgc) + print *, 'Error FFD dTgc: ', norm2(FFD - dTgc) + print *, 'Error BFD d2Tgc: ', norm2(BFD2 - d2Tgc) + print *, 'Error CFD d2Tgc: ', norm2(CFD2 - d2Tgc) + print *, 'Error FFD d2Tgc: ', norm2(FFD2 - d2Tgc) + + !> Finalize the NURBS surface object + call surface%finalize() + deallocate(CFD, BFD, FFD, CFD2, BFD2, FFD2) + + + + + + + + + + + + + + + + + + + + !----------------------------------------------------------------------------- + ! Setting up the NURBS surface + !----------------------------------------------------------------------------- + call surface%set_tetragon(L=[5.0_rk, 8.0_rk], nc=[4,4]) + + !----------------------------------------------------------------------------- + ! Creating the NURBS surface + !----------------------------------------------------------------------------- + + !> Generate the NURBS surface with a resolution of 20 + call surface%create(20, 20) + + !----------------------------------------------------------------------------- + ! Finite Difference Derivative Test + !----------------------------------------------------------------------------- + + allocate(CFD(16,2), BFD(16,2), FFD(16,2)) + allocate(CFD2(2*16,2), BFD2(2*16,2), FFD2(2*16,2)) + + tol = 1.0e-6_rk + + Xt(1) = 0.5_rk + Xt(2) = 0.3_rk + call surface%derivative2(Xt=Xt, d2Tgc=d2Tgc, dTgc=dTgc, Tgc=Tgc) + + do i = 1, 2 + Xtm = Xt + Xtm(i) = Xt(i) - tol + call surface%derivative2(Xt=Xtm, d2Tgc=d2Tgcm, dTgc=dTgcm, Tgc=Tgcm) + + Xtp = Xt + Xtp(i) = Xt(i) + tol + call surface%derivative2(Xt=Xtp, d2Tgc=d2Tgcp, dTgc=dTgcp, Tgc=Tgcp) + + BFD(:,i) = (Tgc - Tgcm)/tol + CFD(:,i) = (Tgcp - Tgcm)/(2.0_rk*tol) + FFD(:,i) = (Tgcp - Tgc )/tol + + BFD2(:,i) = reshape((dTgc - dTgcm)/tol, shape=[2*16]) + CFD2(:,i) = reshape((dTgcp - dTgcm)/(2.0_rk*tol), shape=[2*16]) + FFD2(:,i) = reshape((dTgcp - dTgc )/tol, shape=[2*16]) + end do + + print *, 'Tolerance:', tol + print *, 'Error BFD dTgc: ', norm2(BFD - dTgc) + print *, 'Error CFD dTgc: ', norm2(CFD - dTgc) + print *, 'Error FFD dTgc: ', norm2(FFD - dTgc) + print *, 'Error BFD d2Tgc: ', norm2(BFD2 - d2Tgc) + print *, 'Error CFD d2Tgc: ', norm2(CFD2 - d2Tgc) + print *, 'Error FFD d2Tgc: ', norm2(FFD2 - d2Tgc) + + !> Finalize the NURBS surface object + call surface%finalize() + deallocate(CFD, BFD, FFD, CFD2, BFD2, FFD2) + +end program diff --git a/test/fdm_volume.f90 b/test/fdm_volume.f90 new file mode 100644 index 000000000..8d890c9e1 --- /dev/null +++ b/test/fdm_volume.f90 @@ -0,0 +1,145 @@ +program fdm_test_volume + + use forcad, only: rk, nurbs_volume + + implicit none + + type(nurbs_volume) :: volume !! Declare a NURBS volume object + real(rk), allocatable :: Xc(:,:), Wc(:) !! Arrays for control points and weights + real(rk) :: Xt(3), tol, Xtm(3), Xtp(3) + real(rk), allocatable :: Tgc(:), dTgc(:,:), Tgcp(:), dTgcp(:,:), Tgcm(:), dTgcm(:,:), d2Tgc(:,:), d2Tgcp(:,:), d2Tgcm(:,:) + real(rk), allocatable :: CFD(:,:), BFD(:,:), FFD(:,:), CFD2(:,:), BFD2(:,:), FFD2(:,:) + integer :: i + + !----------------------------------------------------------------------------- + ! Setting up the NURBS volume + !----------------------------------------------------------------------------- + allocate(Wc(64)) + Wc = 1.0_rk + Wc(10) = 0.5_rk + call volume%set_hexahedron(L=[2.0_rk, 4.0_rk, 8.0_rk], nc=[4,4,4], Wc=Wc) + + !----------------------------------------------------------------------------- + ! Creating the NURBS volume + !----------------------------------------------------------------------------- + + !> Generate the NURBS volume with a resolution of 20 + call volume%create(20, 20, 20) + + !----------------------------------------------------------------------------- + ! Finite Difference Derivative Test + !----------------------------------------------------------------------------- + + allocate(CFD(64,3), BFD(64,3), FFD(64,3)) + allocate(CFD2(3*64,3), BFD2(3*64,3), FFD2(3*64,3)) + + tol = 1.0e-6_rk + + Xt(1) = 0.5_rk + Xt(2) = 0.3_rk + Xt(3) = 0.7_rk + call volume%derivative2(Xt=Xt, d2Tgc=d2Tgc, dTgc=dTgc, Tgc=Tgc) + do i = 1, 3 + Xtm = Xt + Xtm(i) = Xt(i) - tol + call volume%derivative2(Xt=Xtm, d2Tgc=d2Tgcm, dTgc=dTgcm, Tgc=Tgcm) + + Xtp = Xt + Xtp(i) = Xt(i) + tol + call volume%derivative2(Xt=Xtp, d2Tgc=d2Tgcp, dTgc=dTgcp, Tgc=Tgcp) + + BFD(:,i) = (Tgc - Tgcm)/tol + CFD(:,i) = (Tgcp - Tgcm)/(2.0_rk*tol) + FFD(:,i) = (Tgcp - Tgc )/tol + + BFD2(:,i) = reshape((dTgc - dTgcm)/tol, shape=[3*64]) + CFD2(:,i) = reshape((dTgcp - dTgcm)/(2.0_rk*tol), shape=[3*64]) + FFD2(:,i) = reshape((dTgcp - dTgc )/tol, shape=[3*64]) + end do + + print *, 'Tolerance:', tol + print *, 'Error BFD dTgc: ', norm2(BFD - dTgc) + print *, 'Error CFD dTgc: ', norm2(CFD - dTgc) + print *, 'Error FFD dTgc: ', norm2(FFD - dTgc) + print *, 'Error BFD d2Tgc: ', norm2(BFD2 - d2Tgc) + print *, 'Error CFD d2Tgc: ', norm2(CFD2 - d2Tgc) + print *, 'Error FFD d2Tgc: ', norm2(FFD2 - d2Tgc) + + !> Finalize the NURBS volume object + call volume%finalize() + deallocate(CFD, BFD, FFD, CFD2, BFD2, FFD2) + deallocate(Wc) + + + + + + + + + + + + + + + + + + + !----------------------------------------------------------------------------- + ! Setting up the NURBS volume + !----------------------------------------------------------------------------- + call volume%set_hexahedron(L=[2.0_rk, 4.0_rk, 8.0_rk], nc=[4,4,4]) + + !----------------------------------------------------------------------------- + ! Creating the NURBS volume + !----------------------------------------------------------------------------- + + !> Generate the NURBS volume with a resolution of 20 + call volume%create(20, 20, 20) + + !----------------------------------------------------------------------------- + ! Finite Difference Derivative Test + !----------------------------------------------------------------------------- + + allocate(CFD(64,3), BFD(64,3), FFD(64,3)) + allocate(CFD2(3*64,3), BFD2(3*64,3), FFD2(3*64,3)) + + tol = 1.0e-6_rk + + Xt(1) = 0.5_rk + Xt(2) = 0.3_rk + Xt(3) = 0.7_rk + call volume%derivative2(Xt=Xt, d2Tgc=d2Tgc, dTgc=dTgc, Tgc=Tgc) + do i = 1, 3 + Xtm = Xt + Xtm(i) = Xt(i) - tol + call volume%derivative2(Xt=Xtm, d2Tgc=d2Tgcm, dTgc=dTgcm, Tgc=Tgcm) + + Xtp = Xt + Xtp(i) = Xt(i) + tol + call volume%derivative2(Xt=Xtp, d2Tgc=d2Tgcp, dTgc=dTgcp, Tgc=Tgcp) + + BFD(:,i) = (Tgc - Tgcm)/tol + CFD(:,i) = (Tgcp - Tgcm)/(2.0_rk*tol) + FFD(:,i) = (Tgcp - Tgc )/tol + + BFD2(:,i) = reshape((dTgc - dTgcm)/tol, shape=[3*64]) + CFD2(:,i) = reshape((dTgcp - dTgcm)/(2.0_rk*tol), shape=[3*64]) + FFD2(:,i) = reshape((dTgcp - dTgc )/tol, shape=[3*64]) + end do + + print *, 'Tolerance:', tol + print *, 'Error BFD dTgc: ', norm2(BFD - dTgc) + print *, 'Error CFD dTgc: ', norm2(CFD - dTgc) + print *, 'Error FFD dTgc: ', norm2(FFD - dTgc) + print *, 'Error BFD d2Tgc: ', norm2(BFD2 - d2Tgc) + print *, 'Error CFD d2Tgc: ', norm2(CFD2 - d2Tgc) + print *, 'Error FFD d2Tgc: ', norm2(FFD2 - d2Tgc) + + !> Finalize the NURBS volume object + call volume%finalize() + deallocate(CFD, BFD, FFD, CFD2, BFD2, FFD2) + +end program diff --git a/test/test_nurbs_curve.f90 b/test/test_nurbs_curve.f90 index 9711c1d3d..b3fe3661b 100644 --- a/test/test_nurbs_curve.f90 +++ b/test/test_nurbs_curve.f90 @@ -167,11 +167,11 @@ program test_nurbs_curve call nurbs%basis(Xt=[(real(i-1, rk) / real(23-1, rk), i=1, 23)], Tgc=Tgc) call bsp%basis(Xt=[(real(i-1, rk) / real(23-1, rk), i=1, 23)], Tgc=Tgc) - call nurbs%derivative(res=23, dTgc=dTgc) - call bsp%derivative(res=23, dTgc=dTgc) + call nurbs%derivative(res=23, dTgc=dTgc, Tgc=Tgc) + call bsp%derivative(res=23, dTgc=dTgc, Tgc=Tgc) - call nurbs%derivative(Xt=[(real(i-1, rk) / real(23-1, rk), i=1, 23)], dTgc=dTgc) - call bsp%derivative(Xt=[(real(i-1, rk) / real(23-1, rk), i=1, 23)], dTgc=dTgc) + call nurbs%derivative(Xt=[(real(i-1, rk) / real(23-1, rk), i=1, 23)], dTgc=dTgc, Tgc=Tgc) + call bsp%derivative(Xt=[(real(i-1, rk) / real(23-1, rk), i=1, 23)], dTgc=dTgc, Tgc=Tgc) call nurbs%rotate_Xc(45.0_rk, 0.0_rk, 0.0_rk) call nurbs%rotate_Xc(-45.0_rk, 0.0_rk, 0.0_rk)