forcad_nurbs_volume.F90 Source File


This file depends on

sourcefile~~forcad_nurbs_volume.f90~~EfferentGraph sourcefile~forcad_nurbs_volume.f90 forcad_nurbs_volume.F90 sourcefile~forcad_interface.f90 forcad_interface.F90 sourcefile~forcad_nurbs_volume.f90->sourcefile~forcad_interface.f90 sourcefile~forcad_kinds.f90 forcad_kinds.F90 sourcefile~forcad_nurbs_volume.f90->sourcefile~forcad_kinds.f90 sourcefile~forcad_utils.f90 forcad_utils.F90 sourcefile~forcad_nurbs_volume.f90->sourcefile~forcad_utils.f90 sourcefile~forcad_interface.f90->sourcefile~forcad_utils.f90 sourcefile~forcad_utils.f90->sourcefile~forcad_kinds.f90

Files dependent on this one

sourcefile~~forcad_nurbs_volume.f90~~AfferentGraph sourcefile~forcad_nurbs_volume.f90 forcad_nurbs_volume.F90 sourcefile~forcad.f90 forcad.f90 sourcefile~forcad.f90->sourcefile~forcad_nurbs_volume.f90 sourcefile~cmp_area.f90 cmp_area.f90 sourcefile~cmp_area.f90->sourcefile~forcad.f90 sourcefile~cmp_length.f90 cmp_length.f90 sourcefile~cmp_length.f90->sourcefile~forcad.f90 sourcefile~cmp_volume.f90 cmp_volume.f90 sourcefile~cmp_volume.f90->sourcefile~forcad.f90 sourcefile~demo_curve.f90 demo_curve.f90 sourcefile~demo_curve.f90->sourcefile~forcad.f90 sourcefile~demo_surface.f90 demo_surface.f90 sourcefile~demo_surface.f90->sourcefile~forcad.f90 sourcefile~demo_volume.f90 demo_volume.f90 sourcefile~demo_volume.f90->sourcefile~forcad.f90 sourcefile~example_bend_pipe.f90 example_bend_pipe.f90 sourcefile~example_bend_pipe.f90->sourcefile~forcad.f90 sourcefile~example_curve_1.f90 example_curve_1.f90 sourcefile~example_curve_1.f90->sourcefile~forcad.f90 sourcefile~example_helix_pipe.f90 example_helix_pipe.f90 sourcefile~example_helix_pipe.f90->sourcefile~forcad.f90 sourcefile~example_ppm1.f90 example_ppm1.f90 sourcefile~example_ppm1.f90->sourcefile~forcad.f90 sourcefile~example_ppm2.f90 example_ppm2.f90 sourcefile~example_ppm2.f90->sourcefile~forcad.f90 sourcefile~example_ppm3.f90 example_ppm3.f90 sourcefile~example_ppm3.f90->sourcefile~forcad.f90 sourcefile~example_surface_1.f90 example_surface_1.f90 sourcefile~example_surface_1.f90->sourcefile~forcad.f90 sourcefile~example_toroidal_pipe.f90 example_toroidal_pipe.f90 sourcefile~example_toroidal_pipe.f90->sourcefile~forcad.f90 sourcefile~example_twist_taper.f90 example_twist_taper.f90 sourcefile~example_twist_taper.f90->sourcefile~forcad.f90 sourcefile~example_volume_1.f90 example_volume_1.f90 sourcefile~example_volume_1.f90->sourcefile~forcad.f90 sourcefile~fdm_curve.f90 fdm_curve.f90 sourcefile~fdm_curve.f90->sourcefile~forcad.f90 sourcefile~fdm_elevate_and_insert_1d.f90 fdm_elevate_and_insert_1d.f90 sourcefile~fdm_elevate_and_insert_1d.f90->sourcefile~forcad.f90 sourcefile~fdm_elevate_and_insert_2d.f90 fdm_elevate_and_insert_2d.f90 sourcefile~fdm_elevate_and_insert_2d.f90->sourcefile~forcad.f90 sourcefile~fdm_elevate_and_insert_3d.f90 fdm_elevate_and_insert_3d.f90 sourcefile~fdm_elevate_and_insert_3d.f90->sourcefile~forcad.f90 sourcefile~fdm_surface.f90 fdm_surface.f90 sourcefile~fdm_surface.f90->sourcefile~forcad.f90 sourcefile~fdm_volume.f90 fdm_volume.f90 sourcefile~fdm_volume.f90->sourcefile~forcad.f90 sourcefile~lsq_fit_bspline_1d.f90 lsq_fit_bspline_1d.f90 sourcefile~lsq_fit_bspline_1d.f90->sourcefile~forcad.f90 sourcefile~lsq_fit_bspline_2d.f90 lsq_fit_bspline_2d.f90 sourcefile~lsq_fit_bspline_2d.f90->sourcefile~forcad.f90 sourcefile~lsq_fit_bspline_3d.f90 lsq_fit_bspline_3d.f90 sourcefile~lsq_fit_bspline_3d.f90->sourcefile~forcad.f90 sourcefile~nearest_point_1d.f90 nearest_point_1d.f90 sourcefile~nearest_point_1d.f90->sourcefile~forcad.f90 sourcefile~nearest_point_2d.f90 nearest_point_2d.f90 sourcefile~nearest_point_2d.f90->sourcefile~forcad.f90 sourcefile~nearest_point_2d_bench.f90 nearest_point_2d_bench.f90 sourcefile~nearest_point_2d_bench.f90->sourcefile~forcad.f90 sourcefile~nearest_point_3d.f90 nearest_point_3d.f90 sourcefile~nearest_point_3d.f90->sourcefile~forcad.f90 sourcefile~poisson_iga_solver_2d.f90 poisson_iga_solver_2d.f90 sourcefile~poisson_iga_solver_2d.f90->sourcefile~forcad.f90 sourcefile~poisson_iga_solver_3d.f90 poisson_iga_solver_3d.f90 sourcefile~poisson_iga_solver_3d.f90->sourcefile~forcad.f90 sourcefile~put_to_nurbs.f90 put_to_nurbs.f90 sourcefile~put_to_nurbs.f90->sourcefile~forcad.f90 sourcefile~shape_c_1d.f90 shape_C_1d.f90 sourcefile~shape_c_1d.f90->sourcefile~forcad.f90 sourcefile~shape_c_2d.f90 shape_C_2d.f90 sourcefile~shape_c_2d.f90->sourcefile~forcad.f90 sourcefile~shape_c_3d.f90 shape_C_3d.f90 sourcefile~shape_c_3d.f90->sourcefile~forcad.f90 sourcefile~shape_circle.f90 shape_circle.f90 sourcefile~shape_circle.f90->sourcefile~forcad.f90 sourcefile~shape_half_circle.f90 shape_half_circle.f90 sourcefile~shape_half_circle.f90->sourcefile~forcad.f90 sourcefile~shape_half_ring_2d.f90 shape_half_ring_2d.f90 sourcefile~shape_half_ring_2d.f90->sourcefile~forcad.f90 sourcefile~shape_half_ring_3d.f90 shape_half_ring_3d.f90 sourcefile~shape_half_ring_3d.f90->sourcefile~forcad.f90 sourcefile~shape_hexahedron.f90 shape_hexahedron.f90 sourcefile~shape_hexahedron.f90->sourcefile~forcad.f90 sourcefile~shape_ring_2d.f90 shape_ring_2d.f90 sourcefile~shape_ring_2d.f90->sourcefile~forcad.f90 sourcefile~shape_ring_3d.f90 shape_ring_3d.f90 sourcefile~shape_ring_3d.f90->sourcefile~forcad.f90 sourcefile~shape_tetragon.f90 shape_tetragon.f90 sourcefile~shape_tetragon.f90->sourcefile~forcad.f90 sourcefile~test_nurbs_surface.f90 test_nurbs_surface.f90 sourcefile~test_nurbs_surface.f90->sourcefile~forcad.f90 sourcefile~test_nurbs_volume.f90 test_nurbs_volume.f90 sourcefile~test_nurbs_volume.f90->sourcefile~forcad.f90

Source Code

!> author: Seyed Ali Ghasemi
!> license: BSD 3-Clause
!> This module defines the 'nurbs_volume' type for representing a Non-Uniform Rational B-Spline (NURBS) volume.
module forcad_nurbs_volume

    use forcad_kinds, only: rk
    use forcad_utils, only: 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, remove_knots_A_5_8, &
        elemConn_Cn, unique, rotation, det, inv, gauss_leg, export_vtk_legacy, basis_bspline_2der
    use fordebug, only: debug

    implicit none

    private
    public nurbs_volume, compute_Tgc, compute_dTgc

    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    type nurbs_volume
        real(rk), allocatable, private :: Xc(:,:)  !! Control points (2D array: [nc(1)*nc(2)*nc(3), dim])
        real(rk), allocatable, private :: Xg(:,:)  !! Geometry points (2D array: [ng(1)*ng(2)*ng(3), dim])
        real(rk), allocatable, private :: Wc(:)    !! Weights for the control points (1D array: [nc(1)*nc(2)*nc(3)])
        real(rk), allocatable, private :: Xt1(:)   !! Evaluation parameter values in the first direction (1D array: [ng(1)])
        real(rk), allocatable, private :: Xt2(:)   !! Evaluation parameter values in the second direction (1D array: [ng(2)])
        real(rk), allocatable, private :: Xt3(:)   !! Evaluation parameter values in the third direction (1D array: [ng(3)])
        real(rk), allocatable, private :: Xt(:,:)  !! Evaluation parameter values (2D array: [ng(1)*ng(2)*ng(3), dim]
        real(rk), allocatable, private :: knot1(:) !! Knot vector in the first direction (1D array)
        real(rk), allocatable, private :: knot2(:) !! Knot vector in the second direction (1D array)
        real(rk), allocatable, private :: knot3(:) !! Knot vector in the third direction (1D array)
        integer, private :: degree(3)              !! Degree (order) of the volume
        integer, private :: nc(3)                  !! Number of control points in each direction
        integer, private :: ng(3)                  !! Number of geometry points in each direction
        integer, allocatable, private :: elemConn_Xc_vis(:,:) !! Connectivity for visualization of control points
        integer, allocatable, private :: elemConn_Xg_vis(:,:) !! Connectivity for visualization of geometry points
        integer, allocatable, private :: elemConn(:,:)        !! IGA element connectivity

        type(debug) :: err !! 101: size mismatch (weights vs control points), 102: missing control points, 103: missing knot vector, 104: missing geometry points, 105: missing weights, 106: lsq fit underdetermined
    contains
        procedure, private :: set1                   !!> Set knot vectors, control points and weights for the NURBS volume object
        procedure, private :: set2                   !!> Set NURBS volume using nodes of parameter space, degree, continuity, control points and weights
        procedure, private :: set3                   !!> Set Bezier or Rational Bezier volume using control points and weights
        procedure, private :: set4                   !!> Set NURBS volume using degree, number of control points, control points and weights
        generic :: set => set1, set2, set3, set4  !!> 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
        generic :: get_Xc => get_Xc_all, get_Xci, get_Xcid !!> Get control points
        procedure, private :: get_Xg_all   !!> Get all geometry points
        procedure, private :: get_Xgi      !!> Get i-th geometry point
        procedure, private :: get_Xgid     !!> Get i-th geometry point in a specific direction
        generic :: get_Xg => get_Xg_all, get_Xgi, get_Xgid !!> Get geometry points
        procedure, private :: get_Wc_all   !!> Get all weights
        procedure, private :: get_Wci      !!> Get i-th weight
        generic :: get_Wc => get_Wc_all, get_Wci !!> Get weights
        procedure :: get_Xt                 !!> Get parameter values
        procedure, private :: get_knot_all  !!> Get all knot vectors
        procedure, private :: get_knoti     !!> Get i-th knot value
        generic :: get_knot => get_knoti, get_knot_all !!> Get knot vector
        procedure :: get_ng                 !!> Get number of geometry points
        procedure, private :: get_nc_dir             !!> Get number of control points in a specific direction
        procedure, private :: get_nc_all             !!> Get number of control points in all directions
        generic :: get_nc => get_nc_all, get_nc_dir !!> Get number of control points
        procedure :: cmp_degree             !!> Compute degree of the NURBS volume
        procedure, private :: get_degree_all!!> Get degree of the NURBS volume in all directions
        procedure, private :: get_degree_dir!!> Get degree of the NURBS volume in a specific direction
        generic :: get_degree => get_degree_all, get_degree_dir !!> Get degree of the NURBS volume
        procedure :: finalize               !!> Finalize the NURBS volume object
        procedure :: cmp_elem_Xc_vis        !!> Generate connectivity for control points
        procedure :: cmp_elem_Xg_vis        !!> Generate connectivity for geometry points
        procedure :: cmp_elem_Xth           !!> Generate connectivity for parameter points
        procedure :: cmp_elem               !!> Generate IGA element connectivity
        procedure :: get_elem_Xc_vis        !!> Get connectivity for control points
        procedure :: get_elem_Xg_vis        !!> Get connectivity for geometry points
        procedure :: get_elem               !!> Get IGA element connectivity
        procedure :: set_elem_Xc_vis        !!> Set connectivity for control points
        procedure :: set_elem_Xg_vis        !!> Set connectivity for geometry points
        procedure :: set_elem               !!> Set IGA element connectivity
        procedure :: export_Xc              !!> Export control points to VTK file
        procedure :: export_Xg              !!> Export geometry points to VTK file
        procedure :: export_Xth             !!> Export parameter space to VTK file
        procedure :: export_Xth_in_Xg       !!> Export parameter space in geometry points to VTK file
        procedure :: modify_Xc              !!> Modify control points
        procedure :: modify_Wc              !!> Modify weights
        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, 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
        procedure :: put_to_nurbs           !!> Put a shape to a NURBS volume
        procedure :: remove_knots           !!> Remove knots from the knot vector
        procedure :: rotate_Xc              !!> Rotate control points
        procedure :: rotate_Xg              !!> Rotate geometry points
        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 (Approximation)
        procedure :: nearest_point2         !!> Find the nearest point on the NURBS volume (Minimization - Newtons method)
        procedure :: ansatz                 !!> Compute the shape functions, derivative of shape functions and dV
        procedure :: cmp_volume             !!> Compute the volume of the NURBS volume
        procedure :: lsq_fit_bspline        !!> Fit B-spline volume to structured data points using least squares

        ! Faces
        procedure :: cmp_elemFace_Xc_vis   !!> Compute faces of the control points
        procedure :: cmp_elemFace_Xg_vis   !!> Compute faces of the geometry points
        procedure :: cmp_elemFace          !!> Compute faces of the IGA elements
        procedure :: cmp_degreeFace        !!> Compute degrees of the faces

        ! Shapes
        procedure :: set_hexahedron         !!> Set a hexahedron
        procedure :: set_ring               !!> Set a ring
        procedure :: set_half_ring          !!> Set a half ring
        procedure :: set_C                  !!> Set a C-shape
    end type
    !===============================================================================

    interface compute_Xg
        module procedure compute_Xg_nurbs_3d
        module procedure compute_Xg_bspline_3d
        module procedure compute_Xg_nurbs_3d_1point
        module procedure compute_Xg_bspline_3d_1point
    end interface

    interface compute_Tgc
        module procedure compute_Tgc_nurbs_3d_vector
        module procedure compute_Tgc_bspline_3d_vector
        module procedure compute_Tgc_nurbs_3d_scalar
        module procedure compute_Tgc_bspline_3d_scalar
    end interface

    interface compute_dTgc
        module procedure compute_dTgc_nurbs_3d_vector
        module procedure compute_dTgc_bspline_3d_vector
        module procedure compute_dTgc_nurbs_3d_scalar
        module procedure compute_dTgc_bspline_3d_scalar
    end interface

    interface compute_d2Tgc
        module procedure compute_d2Tgc_nurbs_3d_vector
        module procedure compute_d2Tgc_bspline_3d_vector
        module procedure compute_d2Tgc_nurbs_3d_scalar
        module procedure compute_d2Tgc_bspline_3d_scalar
    end interface

contains

    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    !> Set control points and weights for the NURBS volume object.
    pure subroutine set1(this, knot1, knot2, knot3, Xc, Wc)
        class(nurbs_volume), intent(inout) :: this
        real(rk), intent(in), contiguous :: knot1(:), knot2(:), knot3(:)
        real(rk), intent(in), contiguous :: Xc(:,:)
        real(rk), intent(in), contiguous, optional :: Wc(:)

        if (.not. this%err%ok) return

        if (allocated(this%knot1)) then
            if (size(this%knot1) /= size(knot1)) deallocate(this%knot1)
        end if
        if (allocated(this%knot2)) then
            if (size(this%knot2) /= size(knot2)) deallocate(this%knot2)
        end if
        if (allocated(this%knot3)) then
            if (size(this%knot3) /= size(knot3)) deallocate(this%knot3)
        end if
        if (allocated(this%Xc)) then
            if (size(this%Xc,1) /= size(Xc,1) .or. size(this%Xc,2) /= size(Xc,2)) deallocate(this%Xc)
        end if

        this%knot1 = knot1
        this%knot2 = knot2
        this%knot3 = knot3
        call this%cmp_degree()
        call this%cmp_nc()
        this%Xc = Xc
        if (present(Wc)) then
            if (size(Wc) /= this%nc(1)*this%nc(2)*this%nc(3)) then
                call this%err%set(&
                    code       = 101,&
                    severity   = 1,&
                    category   = 'forcad_nurbs_volume',&
                    message    = 'Weights length mismatch: size(Wc) must equal number of control points.',&
                    location   = 'set1',&
                    suggestion = 'Provide Wc with size(Wc) == nc(1)*nc(2)*nc(3).')
                return
            else
                if (allocated(this%Wc)) then
                    if (size(this%Wc) /= size(Wc)) deallocate(this%Wc)
                end if
                this%Wc = Wc
            end if
        end if
    end subroutine
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    !> Set control points and weights for the NURBS volume object.
    pure subroutine set2(this, Xth_dir1, Xth_dir2, Xth_dir3, degree, continuity1, continuity2, continuity3, Xc, Wc)
        class(nurbs_volume), intent(inout) :: this
        real(rk), intent(in), contiguous :: Xth_dir1(:), Xth_dir2(:), Xth_dir3(:)
        integer, intent(in), contiguous :: degree(:)
        integer, intent(in), contiguous :: continuity1(:), continuity2(:), continuity3(:)
        real(rk), intent(in), contiguous, optional :: Xc(:,:)
        real(rk), intent(in), contiguous, optional :: Wc(:)

        if (.not. this%err%ok) return

        if (allocated(this%knot1)) deallocate(this%knot1)
        if (allocated(this%knot2)) deallocate(this%knot2)
        if (allocated(this%knot3)) deallocate(this%knot3)
        this%knot1 = compute_knot_vector(Xth_dir1, degree(1), continuity1)
        this%knot2 = compute_knot_vector(Xth_dir2, degree(2), continuity2)
        this%knot3 = compute_knot_vector(Xth_dir3, degree(3), continuity3)
        this%degree(1) = degree(1)
        this%degree(2) = degree(2)
        this%degree(3) = degree(3)
        call this%cmp_nc()

        if (present(Xc)) then
            if (size(Xc,1) /= this%nc(1)*this%nc(2)*this%nc(3)) then
                call this%err%set(&
                    code       = 101,&
                    severity   = 1,&
                    category   = 'forcad_nurbs_volume', &
                    message    = 'Control points size mismatch in set2',&
                    location   = 'set2', &
                    suggestion = 'size(Xc,1) must equal nc(1)*nc(2)*nc(3).')
                return
            end if
        end if
        if (present(Wc)) then
            if (size(Wc) /= this%nc(1)*this%nc(2)*this%nc(3)) then
                call this%err%set(&
                    code       = 101,&
                    severity   = 1,&
                    category   = 'forcad_nurbs_volume', &
                    message    = 'Weights size mismatch in set2',&
                    location   = 'set2', &
                    suggestion = 'size(Wc) must equal nc(1)*nc(2)*nc(3).')
                return
            end if
        end if

        if (present(Xc)) then
            if (allocated(this%Xc)) then
                if (size(this%Xc,1) /= size(Xc,1) .or. size(this%Xc,2) /= size(Xc,2)) deallocate(this%Xc)
            end if
            this%Xc = Xc
        end if
        if (present(Wc)) then
            if (allocated(this%Wc)) then
                if (size(this%Wc) /= size(Wc)) deallocate(this%Wc)
            end if
            this%Wc = Wc
        end if
    end subroutine
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    !> Set Bezier or Rational Bezier volume using control points and weights.
    pure subroutine set3(this, nc, Xc, Wc)
        class(nurbs_volume), intent(inout) :: this
        integer, intent(in), contiguous :: nc(:)
        real(rk), intent(in), contiguous :: Xc(:,:)
        real(rk), intent(in), contiguous, optional :: Wc(:)

        if (.not. this%err%ok) return

        if (allocated(this%Xc)) then
            if (size(this%Xc,1) /= nc(1)*nc(2)*nc(3) .or. size(this%Xc,2) /= size(Xc,2)) deallocate(this%Xc)
        end if

        this%Xc = Xc
        this%nc = nc

        if (allocated(this%knot1)) then
            if (size(this%knot1) /= 2*nc(1)) then
                deallocate(this%knot1)
                allocate(this%knot1(2*this%nc(1)))
            end if
        else
            allocate(this%knot1(2*this%nc(1)))
        end if
        this%knot1(1:this%nc(1)) = 0.0_rk
        this%knot1(this%nc(1)+1:2*this%nc(1)) = 1.0_rk

        if (allocated(this%knot2)) then
            if (size(this%knot2) /= 2*nc(2)) then
                deallocate(this%knot2)
                allocate(this%knot2(2*this%nc(2)))
            end if
        else
            allocate(this%knot2(2*this%nc(2)))
        end if
        this%knot2(1:this%nc(2)) = 0.0_rk
        this%knot2(this%nc(2)+1:2*this%nc(2)) = 1.0_rk

        if (allocated(this%knot3)) then
            if (size(this%knot3) /= 2*nc(3)) then
                deallocate(this%knot3)
                allocate(this%knot3(2*this%nc(3)))
            end if
        else
            allocate(this%knot3(2*this%nc(3)))
        end if
        this%knot3(1:this%nc(3)) = 0.0_rk
        this%knot3(this%nc(3)+1:2*this%nc(3)) = 1.0_rk

        call this%cmp_degree()
        if (present(Wc)) then
            if (size(Wc) /= this%nc(1)*this%nc(2)*this%nc(3)) then
                call this%err%set(&
                    code       = 101,&
                    severity   = 1,&
                    category   = 'forcad_nurbs_volume',&
                    message    = 'Weights length mismatch: size(Wc) must equal number of control points.',&
                    location   = 'set3',&
                    suggestion = 'Provide Wc with size(Wc) == nc(1)*nc(2)*nc(3).')
                return
            else
                if (allocated(this%Wc)) then
                    if (size(this%Wc) /= size(Wc)) deallocate(this%Wc)
                end if
                this%Wc = Wc
            end if
        end if
    end subroutine
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure subroutine set4(this, degree, nc, Xc, Wc)
        class(nurbs_volume), intent(inout) :: this
        integer, intent(in), contiguous :: degree(:)
        integer, intent(in), contiguous :: nc(:)
        real(rk), intent(in), contiguous :: Xc(:,:)
        real(rk), intent(in), contiguous, optional :: Wc(:)
        integer :: m(3), i

        if (.not. this%err%ok) return

        if (allocated(this%Xc)) then
            if (size(this%Xc,1) /= nc(1)*nc(2)*nc(3) .or. size(this%Xc,2) /= size(Xc,2)) deallocate(this%Xc)
        end if

        this%Xc = Xc
        this%nc = nc
        this%degree = degree

        ! Size of knot vectors
        m = nc + degree + 1

        if (allocated(this%knot1)) then
            if (size(this%knot1) /= m(1)) then
                deallocate(this%knot1)
                allocate(this%knot1(m(1)))
            end if
        else
            allocate(this%knot1(m(1)))
        end if
        this%knot1(1:degree(1)+1) = 0.0_rk
        this%knot1(degree(1)+2:m(1)-degree(1)-1) = [(real(i, rk)/(m(1)-2*degree(1)-1), i=1, m(1)-2*degree(1)-2)]
        this%knot1(m(1)-degree(1):m(1)) = 1.0_rk

        if (allocated(this%knot2)) then
            if (size(this%knot2) /= m(2)) then
                deallocate(this%knot2)
                allocate(this%knot2(m(2)))
            end if
        else
            allocate(this%knot2(m(2)))
        end if
        this%knot2(1:degree(2)+1) = 0.0_rk
        this%knot2(degree(2)+2:m(2)-degree(2)-1) = [(real(i, rk)/(m(2)-2*degree(2)-1), i=1, m(2)-2*degree(2)-2)]
        this%knot2(m(2)-degree(2):m(2)) = 1.0_rk

        if (allocated(this%knot3)) then
            if (size(this%knot3) /= m(3)) then
                deallocate(this%knot3)
                allocate(this%knot3(m(3)))
            end if
        else
            allocate(this%knot3(m(3)))
        end if
        this%knot3(1:degree(3)+1) = 0.0_rk
        this%knot3(degree(3)+2:m(3)-degree(3)-1) = [(real(i, rk)/(m(3)-2*degree(3)-1), i=1, m(3)-2*degree(3)-2)]
        this%knot3(m(3)-degree(3):m(3)) = 1.0_rk

        if (present(Wc)) then
            if (size(Wc) /= nc(1)*nc(2)*nc(3)) then
                call this%err%set(&
                    code       = 101,&
                    severity   = 1,&
                    category   = 'forcad_nurbs_volume',&
                    message    = 'Weights length mismatch: size(Wc) must equal number of control points.',&
                    location   = 'set4',&
                    suggestion = 'Provide Wc with size(Wc) == nc(1)*nc(2)*nc(3).')
                return
            else
                if (allocated(this%Wc)) then
                    if (size(this%Wc) /= size(Wc)) deallocate(this%Wc)
                end if
                this%Wc = Wc
            end if
        end if
    end subroutine
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure subroutine create(this, res1, res2, res3, Xt1, Xt2, Xt3, Xt)
        class(nurbs_volume), intent(inout) :: this
        integer, intent(in), optional :: res1, res2, res3
        real(rk), intent(in), contiguous, optional :: Xt1(:), Xt2(:), Xt3(:)
        real(rk), intent(in), contiguous, optional :: Xt(:,:)
        integer :: i

        if (.not. this%err%ok) return

        ! check
        if (.not.allocated(this%Xc)) then
            call this%err%set(&
                code       = 102,&
                severity   = 1,&
                category   = 'forcad_nurbs_volume',&
                message    = 'Control points are not set.',&
                location   = 'create',&
                suggestion = 'Call set(...) first before create().')
            return
        end if

        if (.not.allocated(this%knot1) .or. .not.allocated(this%knot2) .or. .not.allocated(this%knot3)) then
            call this%err%set(&
                code       = 103,&
                severity   = 1,&
                category   = 'forcad_nurbs_volume',&
                message    = 'Knot vector is not set.',&
                location   = 'create',&
                suggestion = 'Call set(...) first before create().')
            return
        end if

        ! Set parameter values
        if (present(Xt1)) then
            if (allocated(this%Xt1)) then
                if (size(this%Xt1) /= size(Xt1)) deallocate(this%Xt1)
            end if
            this%Xt1 = Xt1
        elseif (present(res1)) then
            if (allocated(this%Xt1)) then
                if (size(this%Xt1) /= res1) then
                    deallocate(this%Xt1)
                    allocate(this%Xt1(res1))
                end if
            else
                allocate(this%Xt1(res1))
            end if
            this%Xt1 = [(this%knot1(1)+(this%knot1(size(this%knot1))-this%knot1(1))*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)) then
                if (size(this%Xt2) /= size(Xt2)) deallocate(this%Xt2)
            end if
            this%Xt2 = Xt2
        elseif (present(res2)) then
            if (allocated(this%Xt2)) then
                if (size(this%Xt2) /= res2) then
                    deallocate(this%Xt2)
                    allocate(this%Xt2(res2))
                end if
            else
                allocate(this%Xt2(res2))
            end if
            this%Xt2 = [(this%knot2(1)+(this%knot2(size(this%knot2))-this%knot2(1))*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)) then
                if (size(this%Xt3) /= size(Xt3)) deallocate(this%Xt3)
            end if
            this%Xt3 = Xt3
        elseif (present(res3)) then
            if (allocated(this%Xt3)) then
                if (size(this%Xt3) /= res3) then
                    deallocate(this%Xt3)
                    allocate(this%Xt3(res3))
                end if
            else
                allocate(this%Xt3(res3))
            end if
            this%Xt3 = [(this%knot3(1)+(this%knot3(size(this%knot3))-this%knot3(1))*real(i-1,rk)/real(res3-1,rk), i=1, res3)]
            ! else
            ! this%Xt3 = this%Xt3
        end if

        if (present(Xt)) then
            this%Xt = Xt
        else

            ! Set number of geometry points
            this%ng(1) = size(this%Xt1,1)
            this%ng(2) = size(this%Xt2,1)
            this%ng(3) = size(this%Xt3,1)

            call ndgrid(this%Xt1, this%Xt2, this%Xt3, this%Xt)
        end if

        if (allocated(this%Xg)) then
            if (size(this%Xg,1) /= product(this%ng) .or. size(this%Xg,2) /= size(this%Xc,2)) then
                deallocate(this%Xg)
                allocate(this%Xg(product(this%ng), size(this%Xc,2)))
            end if
        else
            allocate(this%Xg(product(this%ng), size(this%Xc,2)))
        end if

        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(:)

        if (.not. this%err%ok) return

        ! 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
    pure function get_Xc_all(this) result(Xc)
        class(nurbs_volume), intent(in) :: this
        real(rk), allocatable :: Xc(:,:)

        if (.not. this%err%ok) return

        if (allocated(this%Xc)) then
            Xc = this%Xc
        else
            error stop 'Control points are not set.'
        end if
    end function
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure function get_Xci(this, n) result(Xc)
        class(nurbs_volume), intent(in) :: this
        integer, intent(in) :: n
        real(rk), allocatable :: Xc(:)

        if (.not. this%err%ok) return

        if (allocated(this%Xc)) then
            if (n<lbound(this%Xc,1) .or. n>ubound(this%Xc,1)) then
                error stop 'Invalid index for control points.'
            end if
            Xc = this%Xc(n,:)
        else
            error stop 'Control points are not set.'
        end if
    end function
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure function get_Xcid(this, n, dir) result(Xc)
        class(nurbs_volume), intent(in) :: this
        integer, intent(in) :: n
        integer, intent(in) :: dir
        real(rk) :: Xc

        if (.not. this%err%ok) return

        if (allocated(this%Xc)) then
            if (n<lbound(this%Xc,1) .or. n>ubound(this%Xc,1)) then
                error stop 'Invalid index for control points.'
            end if
            if (dir<lbound(this%Xc,2) .or. dir>ubound(this%Xc,2)) then
                error stop 'Invalid direction for control points.'
            end if
            Xc = this%Xc(n, dir)
        else
            error stop 'Control points are not set.'
        end if
    end function
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure function get_Xg_all(this) result(Xg)
        class(nurbs_volume), intent(in) :: this
        real(rk), allocatable :: Xg(:,:)

        if (.not. this%err%ok) return

        if (allocated(this%Xg)) then
            Xg = this%Xg
        else
            error stop 'Geometry points are not set.'
        end if
    end function
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure function get_Xgi(this, n) result(Xg)
        class(nurbs_volume), intent(in) :: this
        integer, intent(in) :: n
        real(rk), allocatable :: Xg(:)

        if (.not. this%err%ok) return

        if (allocated(this%Xg)) then
            if (n<lbound(this%Xg,1) .or. n>ubound(this%Xg,1)) then
                error stop 'Invalid index for geometry points.'
            end if
            Xg = this%Xg(n,:)
        else
            error stop 'Control points are not set.'
        end if
    end function
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure function get_Xgid(this, n, dir) result(Xg)
        class(nurbs_volume), intent(in) :: this
        integer, intent(in) :: n
        integer, intent(in) :: dir
        real(rk) :: Xg

        if (.not. this%err%ok) return

        if (allocated(this%Xg)) then
            if (n<lbound(this%Xg,1) .or. n>ubound(this%Xg,1)) then
                error stop 'Invalid index for geometry points.'
            end if
            if (dir<lbound(this%Xg,2) .or. dir>ubound(this%Xg,2)) then
                error stop 'Invalid direction for geometry points.'
            end if
            Xg = this%Xg(n, dir)
        else
            error stop 'Control points are not set.'
        end if
    end function
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure function get_Wc_all(this) result(Wc)
        class(nurbs_volume), intent(in) :: this
        real(rk), allocatable :: Wc(:)

        if (.not. this%err%ok) return

        if (allocated(this%Wc)) then
            Wc = this%Wc
        else
            error stop 'The NURBS volume is not rational or weights are not set.'
        end if
    end function
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure function get_Wci(this, n) result(Wc)
        class(nurbs_volume), intent(in) :: this
        integer, intent(in) :: n
        real(rk) :: Wc

        if (.not. this%err%ok) return

        if (allocated(this%Wc)) then
            if (n<lbound(this%Wc,1) .or. n>ubound(this%Wc,1)) then
                error stop 'Invalid index for weights.'
            end if
            Wc = this%Wc(n)
        else
            error stop 'The NURBS volume is not rational or weights are not set.'
        end if
    end function
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure function get_Xt(this, dir) result(Xt)
        class(nurbs_volume), intent(in) :: this
        integer, intent(in) :: dir
        real(rk), allocatable :: Xt(:)

        if (.not. this%err%ok) return

        if (dir == 1) then
            if (allocated(this%Xt1)) then
                Xt = this%Xt1
            else
                error stop 'Parameter values are not set.'
            end if
        elseif (dir == 2) then
            if (allocated(this%Xt2)) then
                Xt = this%Xt2
            else
                error stop 'Parameter values are not set.'
            end if
        elseif (dir == 3) then
            if (allocated(this%Xt3)) then
                Xt = this%Xt3
            else
                error stop 'Parameter values are not set.'
            end if
        else
            error stop 'Invalid direction for parameter values.'
        end if

    end function
    !===============================================================================


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

        if (.not. this%err%ok) return

        ng = this%ng
    end function
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure subroutine cmp_degree(this, dir)
        class(nurbs_volume), intent(inout) :: this
        integer, intent(in), optional :: dir
        integer, allocatable :: m1(:), m2(:), m3(:)

        if (.not. this%err%ok) return

        if (present(dir)) then
            if (dir == 1) then
                m1 = this%get_multiplicity(1)
                this%degree(1) = m1(1) - 1
            else if (dir == 2) then
                m2 = this%get_multiplicity(2)
                this%degree(2) = m2(1) - 1
            else if (dir == 3) then
                m3 = this%get_multiplicity(3)
                this%degree(3) = m3(1) - 1
            else
                call this%err%set(&
                    code       = 100,&
                    severity   = 1,&
                    category   = 'forcad_nurbs_volume',&
                    message    = 'Invalid direction for degree.',&
                    location   = 'cmp_degree',&
                    suggestion = 'Check the direction argument.')
                return
            end if
        else
            m1 = this%get_multiplicity(1)
            this%degree(1) = m1(1) - 1

            m2 = this%get_multiplicity(2)
            this%degree(2) = m2(1) - 1

            m3 = this%get_multiplicity(3)
            this%degree(3) = m3(1) - 1
        end if

    end subroutine
    !===============================================================================


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

        if (.not. this%err%ok) return

        degree(1) = this%degree(1)
        degree(2) = this%degree(2)
        degree(3) = this%degree(3)
    end function
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure function get_degree_dir(this, dir) result(degree)
        class(nurbs_volume), intent(in) :: this
        integer, intent(in) :: dir
        integer :: degree

        if (.not. this%err%ok) return

        if (dir == 1) then
            degree = this%degree(1)
        else if (dir == 2) then
            degree = this%degree(2)
        else if (dir == 3) then
            degree = this%degree(3)
        else
            error stop 'Invalid direction for degree.'
        end if
    end function
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure function get_knot_all(this, dir) result(knot)
        class(nurbs_volume), intent(in) :: this
        integer, intent(in) :: dir
        real(rk), allocatable :: knot(:)

        if (.not. this%err%ok) return

        if (dir == 1) then
            if (allocated(this%knot1)) then
                knot = this%knot1
            else
                error stop 'Knot vector is not set.'
            end if
        elseif (dir == 2) then
            if (allocated(this%knot2)) then
                knot = this%knot2
            else
                error stop 'Knot vector is not set.'
            end if
        elseif (dir == 3) then
            if (allocated(this%knot3)) then
                knot = this%knot3
            else
                error stop 'Knot vector is not set.'
            end if
        else
            error stop 'Invalid direction for knot vector.'
        end if

    end function
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure function get_knoti(this, dir, i) result(knot)
        class(nurbs_volume), intent(in) :: this
        integer, intent(in) :: dir
        integer, intent(in) :: i
        real(rk) :: knot

        if (.not. this%err%ok) return

        if (dir == 1) then
            if (allocated(this%knot1)) then
                if (i < 1 .or. i > size(this%knot1)) then
                    error stop 'Invalid index for knot vector.'
                else
                    knot = this%knot1(i)
                end if
            else
                error stop 'Knot vector is not set.'
            end if
        elseif (dir == 2) then
            if (allocated(this%knot2)) then
                if (i < 1 .or. i > size(this%knot2)) then
                    error stop 'Invalid index for knot vector.'
                else
                    knot = this%knot2(i)
                end if
            else
                error stop 'Knot vector is not set.'
            end if
        elseif (dir == 3) then
            if (allocated(this%knot3)) then
                if (i < 1 .or. i > size(this%knot3)) then
                    error stop 'Invalid index for knot vector.'
                else
                    knot = this%knot3(i)
                end if
            else
                error stop 'Knot vector is not set.'
            end if
        else
            error stop 'Invalid direction for knot vector.'
        end if

    end function
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure subroutine finalize(this)
        class(nurbs_volume), intent(inout) :: this
        if (allocated(this%Xc)) deallocate(this%Xc)
        if (allocated(this%Xg)) deallocate(this%Xg)
        if (allocated(this%Wc)) deallocate(this%Wc)
        if (allocated(this%Xt1)) deallocate(this%Xt1)
        if (allocated(this%Xt2)) deallocate(this%Xt2)
        if (allocated(this%Xt3)) deallocate(this%Xt3)
        if (allocated(this%Xt)) deallocate(this%Xt)
        if (allocated(this%knot1)) deallocate(this%knot1)
        if (allocated(this%knot2)) deallocate(this%knot2)
        if (allocated(this%knot3)) deallocate(this%knot3)
        if (allocated(this%elemConn_Xc_vis)) deallocate(this%elemConn_Xc_vis)
        if (allocated(this%elemConn_Xg_vis)) deallocate(this%elemConn_Xg_vis)
        if (allocated(this%elemConn)) deallocate(this%elemConn)
    end subroutine
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure function cmp_elem_Xc_vis(this, p) result(elemConn)
        class(nurbs_volume), intent(in) :: this
        integer, allocatable :: elemConn(:,:)
        integer, intent(in), contiguous, optional :: p(:)

        if (.not. this%err%ok) return

        if (present(p)) then
            elemConn = elemConn_C0(this%nc(1), this%nc(2), this%nc(3), p(1), p(2), p(3))
        else
            elemConn = elemConn_C0(this%nc(1), this%nc(2), this%nc(3), 1, 1, 1)
        end if
    end function
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure function cmp_elem_Xg_vis(this, p) result(elemConn)
        class(nurbs_volume), intent(in) :: this
        integer, allocatable :: elemConn(:,:)
        integer, intent(in), contiguous, optional :: p(:)

        if (.not. this%err%ok) return

        if (present(p)) then
            elemConn = elemConn_C0(this%ng(1), this%ng(2), this%ng(3), p(1), p(2), p(3))
        else
            elemConn = elemConn_C0(this%ng(1), this%ng(2), this%ng(3), 1, 1, 1)
        end if
    end function
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure function cmp_elem_Xth(this, p) result(elemConn)
        class(nurbs_volume), intent(in) :: this
        integer, allocatable :: elemConn(:,:)
        integer, intent(in), contiguous, optional :: p(:)

        if (.not. this%err%ok) return

        if (present(p)) then
            elemConn = elemConn_C0(size(unique(this%knot1)), size(unique(this%knot2)), size(unique(this%knot3)), p(1), p(2), p(3))
        else
            elemConn = elemConn_C0(size(unique(this%knot1)), size(unique(this%knot2)), size(unique(this%knot3)), 1, 1, 1)
        end if
    end function
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    impure subroutine export_Xc(this, filename, point_data, field_names, encoding)
        class(nurbs_volume), intent(inout) :: this
        character(len=*), intent(in) :: filename
        real(rk), intent(in), contiguous, optional :: point_data(:,:)
        character(len=*), intent(in), contiguous, optional :: field_names(:)
        character(len=*), intent(in), optional :: encoding
        integer, allocatable :: elemConn(:,:)

        if (.not. this%err%ok) return

        ! check
        if (.not.allocated(this%Xc)) then
            call this%err%set(&
                code       = 102,&
                severity   = 1,&
                category   = 'forcad_nurbs_volume',&
                message    = 'Control points are not set.',&
                location   = 'export_Xc',&
                suggestion = 'Call set(...) first before exporting.')
            return
        end if

        if (.not.allocated(this%elemConn_Xc_vis)) then
            elemConn = this%cmp_elem_Xc_vis()
        else
            elemConn = this%elemConn_Xc_vis
        end if

        call export_vtk_legacy(filename=filename, points=this%Xc, elemConn=elemConn, vtkCellType=12, &
                               point_data=point_data, field_names=field_names, encoding=encoding)
    end subroutine
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    impure subroutine export_Xg(this, filename, point_data, field_names, encoding)
        class(nurbs_volume), intent(inout) :: this
        character(len=*), intent(in) :: filename
        real(rk), intent(in), contiguous, optional :: point_data(:,:)
        character(len=*), intent(in), contiguous, optional :: field_names(:)
        character(len=*), intent(in), optional :: encoding
        integer, allocatable :: elemConn(:,:)

        if (.not. this%err%ok) return

        ! check
        if (.not.allocated(this%Xg)) then
            call this%err%set(&
                code       = 104,&
                severity   = 1,&
                category   = 'forcad_nurbs_volume',&
                message    = 'Geometry points are not set.',&
                location   = 'export_Xg',&
                suggestion = 'Generate Xg by calling create(...) before exporting.')
            return
        end if

        if (.not.allocated(this%elemConn_Xg_vis)) then
            elemConn = this%cmp_elem_Xg_vis()
        else
            elemConn = this%elemConn_Xg_vis
        end if

        call export_vtk_legacy(filename=filename, points=this%Xg, elemConn=elemConn, vtkCellType=12, &
                               point_data=point_data, field_names=field_names, encoding=encoding)
    end subroutine
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    impure subroutine export_Xth(this, filename, point_data, field_names, encoding)
        class(nurbs_volume), intent(in) :: this
        character(len=*), intent(in) :: filename
        real(rk), intent(in), contiguous, optional :: point_data(:,:)
        character(len=*), intent(in), contiguous, optional :: field_names(:)
        character(len=*), intent(in), optional :: encoding
        integer, allocatable :: elemConn(:,:)
        real(rk), allocatable :: Xth(:,:), Xth1(:), Xth2(:), Xth3(:)

        if (.not. this%err%ok) return

        elemConn = this%cmp_elem_Xth()
        Xth1 = unique(this%knot1)
        Xth2 = unique(this%knot2)
        Xth3 = unique(this%knot3)
        call ndgrid(Xth1, Xth2, Xth3, Xth)

        call export_vtk_legacy(filename=filename, points=Xth, elemConn=elemConn, vtkCellType=12, &
                               point_data=point_data, field_names=field_names, encoding=encoding)
    end subroutine
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    impure subroutine export_Xth_in_Xg(this, filename, res, encoding)
        class(nurbs_volume), intent(inout) :: this
        character(len=*),   intent(in) :: filename
        integer, intent(in), optional  :: res   ! min points per span (>=2)
        character(len=*), intent(in), optional :: encoding

        integer :: a, b, t, s, o, r, res_min, dim, i, j, k, ne_total, np, line_idx, m, N1sp, N2sp, N3sp, L, N, res1, res2, res3, offsetP

        real(rk), allocatable :: U1(:), U2(:), U3(:)        ! unique knots
        real(rk), allocatable :: U1r(:), U2r(:), U3r(:)     ! refined per dir (length N)
        real(rk), allocatable :: Xt_all(:,:), Xg_all(:,:)   ! batched params & geometry
        integer,  allocatable :: elemConn(:,:)              ! [ne_total, N]

        if (.not. this%err%ok) return

        if (.not. allocated(this%Xc)) then
            call this%err%set(&
                code       = 102,&
                severity   = 1,&
                category   = 'forcad_nurbs_volume',&
                message    = 'Control points are not set.',&
                location   = 'export_Xth_in_Xg',&
                suggestion = 'Call set(...) first before exporting.')
            return
        end if

        if (.not. allocated(this%knot1) .or. .not. allocated(this%knot2) .or. .not. allocated(this%knot3)) then
            call this%err%set(&
                code       = 103,&
                severity   = 1,&
                category   = 'forcad_nurbs_volume',&
                message    = 'Knot vector is not set.',&
                location   = 'export_Xth_in_Xg',&
                suggestion = 'Call set(...) first before exporting.')
            return
        end if

        res_min = 10
        if (present(res)) res_min = max(2, res)

        U1 = unique(this%knot1)
        if (size(U1) < 2) then
            call this%err%set(&
                code       = 100,&
                severity   = 1,&
                category   = 'forcad_nurbs_volume',&
                message    = 'knot1 needs >= 2 unique values.',&
                location   = 'export_Xth_in_Xg',&
                suggestion = 'Check the knot vector for sufficient unique values.')
            return
        end if
        U2 = unique(this%knot2)
        if (size(U2) < 2) then
            call this%err%set(&
                code       = 100,&
                severity   = 1,&
                category   = 'forcad_nurbs_volume',&
                message    = 'knot2 needs >= 2 unique values.',&
                location   = 'export_Xth_in_Xg',&
                suggestion = 'Check the knot vector for sufficient unique values.')
            return
        end if
        U3 = unique(this%knot3)
        if (size(U3) < 2) then
            call this%err%set(&
                code       = 100,&
                severity   = 1,&
                category   = 'forcad_nurbs_volume',&
                message    = 'knot3 needs >= 2 unique values.',&
                location   = 'export_Xth_in_Xg',&
                suggestion = 'Check the knot vector for sufficient unique values.')
            return
        end if

        ! spans per direction
        N1sp = size(U1)-1
        N2sp = size(U2)-1
        N3sp = size(U3)-1

        L = N1sp
        a = L
        b = N2sp
        do
            t = mod(a,b)
            if (t==0) exit
            a = b
            b = t
        end do
        L = (L / b) * N2sp

        a = L
        b = N3sp
        do
            t = mod(a,b)
            if (t==0) exit
            a = b
            b = t
        end do
        L = (L / b) * N3sp

        L = L * max(1, res_min-1)

        N = L + 1
        res1 = L / N1sp + 1
        res2 = L / N2sp + 1
        res3 = L / N3sp + 1

        dim = size(this%Xc,2)
        if (dim < 2 .or. dim > 3) then
            call this%err%set(&
                code       = 100,&
                severity   = 1,&
                category   = 'forcad_nurbs_volume',&
                message    = 'Invalid geometry dimension.',&
                location   = 'export_Xth_in_Xg',&
                suggestion = 'Check the geometry dimension before exporting the NURBS volume.')
            return
        end if

        ! Allocate refined knot vectors
        allocate(U1r( (size(U1)-1)*(res1-1) + 1 ))
        allocate(U2r( (size(U2)-1)*(res2-1) + 1 ))
        allocate(U3r( (size(U3)-1)*(res3-1) + 1 ))

        do s = 1, size(U1)-1
            o = (s-1)*(res1-1)
            do r = 1, res1
                U1r(o+r) = U1(s) + (U1(s+1)-U1(s)) * real(r-1, rk) / real(res1-1, rk)
            end do
        end do
        do s = 1, size(U2)-1
            o = (s-1)*(res2-1)
            do r = 1, res2
                U2r(o+r) = U2(s) + (U2(s+1)-U2(s)) * real(r-1, rk) / real(res2-1, rk)
            end do
        end do
        do s = 1, size(U3)-1
            o = (s-1)*(res3-1)
            do r = 1, res3
                U3r(o+r) = U3(s) + (U3(s+1)-U3(s)) * real(r-1, rk) / real(res3-1, rk)
            end do
        end do

        if (size(U1r) /= N .or. size(U2r) /= N .or. size(U3r) /= N) then
            call this%err%set(&
                code       = 100,&
                severity   = 1,&
                category   = 'forcad_nurbs_volume',&
                message    = 'Refinement size mismatch.',&
                location   = 'export_Xth_in_Xg',&
                suggestion = 'Check the refinement process for consistency.')
            return
        end if

        ! total element count and node count
        ne_total = size(U2)*size(U3) + size(U1)*size(U3) + size(U1)*size(U2)
        np = ne_total * N

        ! Allocate global arrays
        allocate(Xt_all(np,3), Xg_all(np,dim), elemConn(ne_total, N))

        ! build all parametric points
        line_idx = 0
        offsetP  = 0

        ! dir-1: u varies (v=U2(j), w=U3(k))
        do k = 1, size(U3)
            do j = 1, size(U2)
                line_idx = line_idx + 1
                Xt_all(offsetP+1:offsetP+N,1) = U1r
                Xt_all(offsetP+1:offsetP+N,2) = U2(j)
                Xt_all(offsetP+1:offsetP+N,3) = U3(k)
                offsetP = offsetP + N
            end do
        end do

        ! dir-2: v varies (u=U1(i), w=U3(k))
        do k = 1, size(U3)
            do i = 1, size(U1)
                line_idx = line_idx + 1
                Xt_all(offsetP+1:offsetP+N,1) = U1(i)
                Xt_all(offsetP+1:offsetP+N,2) = U2r
                Xt_all(offsetP+1:offsetP+N,3) = U3(k)
                offsetP = offsetP + N
            end do
        end do

        ! dir-3: w varies (u=U1(i), v=U2(j))
        do j = 1, size(U2)
            do i = 1, size(U1)
                line_idx = line_idx + 1
                Xt_all(offsetP+1:offsetP+N,1) = U1(i)
                Xt_all(offsetP+1:offsetP+N,2) = U2(j)
                Xt_all(offsetP+1:offsetP+N,3) = U3r
                offsetP = offsetP + N
            end do
        end do

        ! compute global points
        if (this%is_rational()) then
            Xg_all = compute_Xg(Xt_all, this%knot1, this%knot2, this%knot3, this%degree, this%nc, &
                                [np,1,1], this%Xc, this%Wc)
        else
            Xg_all = compute_Xg(Xt_all, this%knot1, this%knot2, this%knot3, this%degree, this%nc, &
                                [np,1,1], this%Xc)
        end if

        ! connectivity
        do concurrent (l = 1:ne_total, m = 1:N)
            elemConn(l, m) = (l-1)*N + m
        end do

        ! write VTK file
        call export_vtk_legacy(filename=filename, points=Xg_all, elemConn=elemConn, vtkCellType=4, encoding=encoding)
    end subroutine
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure subroutine modify_Xc(this,X,num,dir)
        class(nurbs_volume), intent(inout) :: this
        real(rk), intent(in) :: X
        integer, intent(in) :: num
        integer, intent(in) :: dir

        if (.not. this%err%ok) return

        if (allocated(this%Xc)) then
            this%Xc(num,dir) = X
            if (allocated(this%Wc)) then
                call this%set(knot1=this%get_knot(1), knot2=this%get_knot(2), knot3=this%get_knot(3),&
                    Xc=this%get_Xc(), Wc=this%get_Wc())
            else
                call this%set(knot1=this%get_knot(1), knot2=this%get_knot(2), knot3=this%get_knot(3),&
                    Xc=this%get_Xc())
            end if
        else
            call this%err%set(&
                code       = 102,&
                severity   = 1,&
                category   = 'forcad_nurbs_volume',&
                message    = 'Control points are not set.',&
                location   = 'modify_Xc',&
                suggestion = 'Call set(...) before modifying it.')
            return
        end if
    end subroutine
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure subroutine modify_Wc(this,W,num)
        class(nurbs_volume), intent(inout) :: this
        real(rk), intent(in) :: W
        integer, intent(in) :: num

        if (.not. this%err%ok) return

        if (allocated(this%Wc)) then
            this%Wc(num) = W
            if (allocated(this%knot1) .and. allocated(this%knot2) .and. allocated(this%knot3)) then
                call this%set(knot1=this%get_knot(1), knot2=this%get_knot(2), knot3=this%get_knot(3),&
                    Xc=this%get_Xc(), Wc=this%get_Wc())
            else
                call this%set(nc=this%nc, Xc=this%get_Xc(), Wc=this%get_Wc())
            end if
        else
            call this%err%set(&
                code       = 105,&
                severity   = 1,&
                category   = 'forcad_nurbs_volume',&
                message    = 'Weights are not set.',&
                location   = 'modify_Wc',&
                suggestion = 'Pass Wc when calling set(...), before modifying weights.')
            return
        end if
    end subroutine
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure function get_multiplicity(this, dir) result(m)
        class(nurbs_volume), intent(in) :: this
        integer, intent(in) :: dir
        integer, allocatable :: m(:)

        if (.not. this%err%ok) return

        if (dir == 1) then

            ! check
            if (.not.allocated(this%knot1)) then
                error stop 'Knot vector is not set.'
            else
                m = compute_multiplicity(this%knot1)
            end if

        elseif (dir == 2) then

            ! check
            if (.not.allocated(this%knot2)) then
                error stop 'Knot vector is not set.'
            else
                m = compute_multiplicity(this%knot2)
            end if

        elseif (dir == 3) then

            ! check
            if (.not.allocated(this%knot3)) then
                error stop 'Knot vector is not set.'
            else
                m = compute_multiplicity(this%knot3)
            end if

        else
            error stop 'Invalid direction.'
        end if

    end function
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure function get_continuity(this, dir) result(c)
        class(nurbs_volume), intent(in) :: this
        integer, intent(in) :: dir
        integer, allocatable :: c(:)

        if (.not. this%err%ok) return

        if (dir == 1) then

            ! check
            if (.not.allocated(this%knot1)) then
                error stop 'Knot vector is not set.'
            else
                c = this%degree(1) - compute_multiplicity(this%knot1)
            end if

        elseif (dir == 2) then

            ! check
            if (.not.allocated(this%knot2)) then
                error stop 'Knot vector is not set.'
            else
                c = this%degree(2) - compute_multiplicity(this%knot2)
            end if

        elseif (dir == 3) then

            ! check
            if (.not.allocated(this%knot3)) then
                error stop 'Knot vector is not set.'
            else
                c = this%degree(3) - compute_multiplicity(this%knot3)
            end if

        else
            error stop 'Invalid direction.'
        end if

    end function
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure subroutine cmp_nc(this, dir)
        class(nurbs_volume), intent(inout) :: this
        integer, intent(in), optional :: dir

        if (.not. this%err%ok) return

        if (present(dir)) then

            if (dir == 1) then

                ! check
                if (.not.allocated(this%knot1)) then
                    call this%err%set(&
                        code       = 103,&
                        severity   = 1,&
                        category   = 'forcad_nurbs_volume',&
                        message    = 'Knot vector is not set.',&
                        location   = 'cmp_nc',&
                        suggestion = 'Call set(...) first before computing nc.' )
                    return
                else
                    this%nc(1) = sum(compute_multiplicity(this%knot1)) - this%degree(1) - 1
                end if

            elseif (dir == 2) then

                ! check
                if (.not.allocated(this%knot2)) then
                    call this%err%set(&
                        code       = 103,&
                        severity   = 1,&
                        category   = 'forcad_nurbs_volume',&
                        message    = 'Knot vector is not set.',&
                        location   = 'cmp_nc',&
                        suggestion = 'Call set(...) first before computing nc.' )
                    return
                else
                    this%nc(2) = sum(compute_multiplicity(this%knot2)) - this%degree(2) - 1
                end if

            elseif (dir == 3) then

                ! check
                if (.not.allocated(this%knot3)) then
                    call this%err%set(&
                        code       = 103,&
                        severity   = 1,&
                        category   = 'forcad_nurbs_volume',&
                        message    = 'Knot vector is not set.',&
                        location   = 'cmp_nc',&
                        suggestion = 'Call set(...) first before computing nc.' )
                    return
                else
                    this%nc(3) = sum(compute_multiplicity(this%knot3)) - this%degree(3) - 1
                end if

            else
                call this%err%set(&
                    code       = 103,&
                    severity   = 1,&
                    category   = 'forcad_nurbs_volume',&
                    message    = 'Invalid direction for computing number of control points.',&
                    location   = 'cmp_nc',&
                    suggestion = 'Use dir=1 or dir=2 to specify the direction.')
                return
            end if

        else

            ! check
            if (.not.allocated(this%knot1)) then
                call this%err%set(&
                    code       = 103,&
                    severity   = 1,&
                    category   = 'forcad_nurbs_volume',&
                    message    = 'Knot vector is not set.',&
                    location   = 'cmp_nc',&
                    suggestion = 'Call set(...) first before computing nc.')
                return
            else
                this%nc(1) = sum(compute_multiplicity(this%knot1)) - this%degree(1) - 1
            end if

            ! check
            if (.not.allocated(this%knot2)) then
                call this%err%set(&
                    code       = 103,&
                    severity   = 1,&
                    category   = 'forcad_nurbs_volume',&
                    message    = 'Knot vector is not set.',&
                    location   = 'cmp_nc',&
                    suggestion = 'Call set(...) first before computing nc.')
                return
            else
                this%nc(2) = sum(compute_multiplicity(this%knot2)) - this%degree(2) - 1
            end if

            ! check
            if (.not.allocated(this%knot3)) then
                call this%err%set(&
                    code       = 103,&
                    severity   = 1,&
                    category   = 'forcad_nurbs_volume',&
                    message    = 'Knot vector is not set.',&
                    location   = 'cmp_nc',&
                    suggestion = 'Call set(...) first before computing nc.')
                return
            else
                this%nc(3) = sum(compute_multiplicity(this%knot3)) - this%degree(3) - 1
            end if

        end if

    end subroutine
    !===============================================================================


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

        if (.not. this%err%ok) return

        nc = this%nc
    end function
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure function get_nc_dir(this, dir) result(nc)
        class(nurbs_volume), intent(in) :: this
        integer, intent(in) :: dir
        integer :: nc

        if (.not. this%err%ok) return

        if (dir == 1) then

            ! check
            if (.not.allocated(this%knot1)) then
                error stop 'Knot vector is not set.'
            else
                nc = sum(compute_multiplicity(this%knot1)) - this%degree(1) - 1
            end if

        elseif (dir == 2) then

            ! check
            if (.not.allocated(this%knot2)) then
                error stop 'Knot vector is not set.'
            else
                nc = sum(compute_multiplicity(this%knot2)) - this%degree(2) - 1
            end if

        elseif (dir == 3) then

            ! check
            if (.not.allocated(this%knot3)) then
                error stop 'Knot vector is not set.'
            else
                nc = sum(compute_multiplicity(this%knot3)) - this%degree(3) - 1
            end if

        else
            error stop 'Invalid direction.'
        end if

    end function
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    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), optional :: Tgc(:,:)
        integer :: i
        real(rk), allocatable :: Xt(:,:)

        if (.not. this%err%ok) return

        ! Set parameter values
        if (present(Xt1)) then
            if (allocated(this%Xt1)) then
                if (size(Xt1,1) /= size(this%Xt1,1)) deallocate(this%Xt1)
            end if
            this%Xt1 = Xt1
        elseif (present(res1)) then
            if (allocated(this%Xt1)) then
                if (size(this%Xt1,1) /= res1) then
                    deallocate(this%Xt1)
                    allocate(this%Xt1(res1))
                end if
            else
                allocate(this%Xt1(res1))
            end if
            this%Xt1 = [(this%knot1(1)+(this%knot1(size(this%knot1))-this%knot1(1))*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)) then
                if (size(Xt2,1) /= size(this%Xt2,1)) deallocate(this%Xt2)
            end if
            this%Xt2 = Xt2
        elseif (present(res2)) then
            if (allocated(this%Xt2)) then
                if (size(this%Xt2,1) /= res2) then
                    deallocate(this%Xt2)
                    allocate(this%Xt2(res2))
                end if
            else
                allocate(this%Xt2(res2))
            end if
            this%Xt2 = [(this%knot2(1)+(this%knot2(size(this%knot2))-this%knot2(1))*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)) then
                if (size(Xt3,1) /= size(this%Xt3,1)) deallocate(this%Xt3)
            end if
            this%Xt3 = Xt3
        elseif (present(res3)) then
            if (allocated(this%Xt3)) then
                if (size(this%Xt3,1) /= res3) then
                    deallocate(this%Xt3)
                    allocate(this%Xt3(res3))
                end if
            else
                allocate(this%Xt3(res3))
            end if
            this%Xt3 = [(this%knot3(1)+(this%knot3(size(this%knot3))-this%knot3(1))*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, elem)
        class(nurbs_volume), intent(inout) :: this
        real(rk), intent(in), contiguous :: Xt(:)
        integer, intent(in), contiguous, optional :: elem(:)
        real(rk), allocatable, intent(out) :: dTgc(:,:)
        real(rk), allocatable, intent(out), optional :: Tgc(:)

        if (.not. this%err%ok) return

        if (this%is_rational()) then ! NURBS
            if (present(elem)) then
                associate(Wce => this%Wc(elem))
                    call compute_dTgc(Xt, this%knot1, this%knot2, this%knot3, this%degree, this%nc, Wce, dTgc, Tgc, elem)
                end associate
            else
                call compute_dTgc(Xt, this%knot1, this%knot2, this%knot3, this%degree, this%nc, this%Wc, dTgc, Tgc)
            end if
        else
            call compute_dTgc(Xt, this%knot1, this%knot2, this%knot3, this%degree, this%nc, dTgc, Tgc, elem)
        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(:,:)

        if (.not. this%err%ok) return

        ! Set parameter values
        if (present(Xt1)) then
            if (allocated(this%Xt1)) then
                if (size(Xt1,1) /= size(this%Xt1,1)) deallocate(this%Xt1)
            end if
            this%Xt1 = Xt1
        elseif (present(res1)) then
            if (allocated(this%Xt1)) then
                if (size(this%Xt1,1) /= res1) then
                    deallocate(this%Xt1)
                    allocate(this%Xt1(res1))
                end if
            else
                allocate(this%Xt1(res1))
            end if
            this%Xt1 = [(this%knot1(1)+(this%knot1(size(this%knot1))-this%knot1(1))*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)) then
                if (size(Xt2,1) /= size(this%Xt2,1)) deallocate(this%Xt2)
            end if
            this%Xt2 = Xt2
        elseif (present(res2)) then
            if (allocated(this%Xt2)) then
                if (size(this%Xt2,1) /= res2) then
                    deallocate(this%Xt2)
                    allocate(this%Xt2(res2))
                end if
            else
                allocate(this%Xt2(res2))
            end if
            this%Xt2 = [(this%knot2(1)+(this%knot2(size(this%knot2))-this%knot2(1))*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)) then
                if (size(Xt3,1) /= size(this%Xt3,1)) deallocate(this%Xt3)
            end if
            this%Xt3 = Xt3
        elseif (present(res3)) then
            if (allocated(this%Xt3)) then
                if (size(this%Xt3,1) /= res3) then
                    deallocate(this%Xt3)
                    allocate(this%Xt3(res3))
                end if
            else
                allocate(this%Xt3(res3))
            end if
            this%Xt3 = [(this%knot3(1)+(this%knot3(size(this%knot3))-this%knot3(1))*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_d2Tgc(Xt, this%knot1, this%knot2, this%knot3, this%degree, this%nc, this%ng, this%Wc, d2Tgc, dTgc, Tgc)
        else
            call compute_d2Tgc(Xt, this%knot1, this%knot2, this%knot3, 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_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 (.not. this%err%ok) return

        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(:)
        real(rk), allocatable, intent(out) :: Tgc(:,:)
        integer :: i
        real(rk), allocatable :: Xt(:,:)

        if (.not. this%err%ok) return

        ! Set parameter values
        if (present(Xt1)) then
            if (allocated(this%Xt1)) then
                if (size(Xt1,1) /= size(this%Xt1,1)) deallocate(this%Xt1)
            end if
            this%Xt1 = Xt1
        elseif (present(res1)) then
            if (allocated(this%Xt1)) then
                if (size(this%Xt1,1) /= res1) then
                    deallocate(this%Xt1)
                    allocate(this%Xt1(res1))
                end if
            else
                allocate(this%Xt1(res1))
            end if
            this%Xt1 = [(this%knot1(1)+(this%knot1(size(this%knot1))-this%knot1(1))*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)) then
                if (size(Xt2,1) /= size(this%Xt2,1)) deallocate(this%Xt2)
            end if
            this%Xt2 = Xt2
        elseif (present(res2)) then
            if (allocated(this%Xt2)) then
                if (size(this%Xt2,1) /= res2) then
                    deallocate(this%Xt2)
                    allocate(this%Xt2(res2))
                end if
            else
                allocate(this%Xt2(res2))
            end if
            this%Xt2 = [(this%knot2(1)+(this%knot2(size(this%knot2))-this%knot2(1))*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)) then
                if (size(Xt3,1) /= size(this%Xt3,1)) deallocate(this%Xt3)
            end if
            this%Xt3 = Xt3
        elseif (present(res3)) then
            if (allocated(this%Xt3)) then
                if (size(this%Xt3,1) /= res3) then
                    deallocate(this%Xt3)
                    allocate(this%Xt3(res3))
                end if
            else
                allocate(this%Xt3(res3))
            end if
            this%Xt3 = [(this%knot3(1)+(this%knot3(size(this%knot3))-this%knot3(1))*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
            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 (.not. this%err%ok) return

        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(Xt, this%knot1, this%knot2, this%knot3, this%degree, this%nc)
        end if
    end subroutine
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure subroutine insert_knots(this, dir ,Xth, r, B, Bs)
        class(nurbs_volume),            intent(inout) :: this
        integer,                        intent(in)    :: dir
        real(rk), contiguous,           intent(in)    :: Xth(:)
        integer,  contiguous,           intent(in)    :: r(:)
        real(rk), allocatable, optional,intent(out)   :: B(:,:)
        real(rk), allocatable, optional,intent(out)   :: Bs(:,:)

        integer :: k, i, s, j, n_new
        real(rk), allocatable :: Xc(:,:), Xcw(:,:), Xcw_new(:,:), knot_new(:)
        real(rk), allocatable :: Xc4(:,:,:,:), H(:,:)
        integer :: nc_old(3), dim, ncp_old, n1_old
        real(rk), allocatable :: Wc_old(:)
        real(rk), allocatable :: A1(:,:), A_re_loc(:,:)

        if (.not. this%err%ok) return

        dim     = size(this%Xc,2)

        if (present(B) .or. present(Bs)) then
            nc_old  = this%nc
            ncp_old = size(this%Xc,1)
            if (this%is_rational()) then
                allocate(Wc_old(ncp_old))
                Wc_old = this%Wc
            end if
            select case(dir)
            case(1)
                n1_old = nc_old(1)
            case(2)
                n1_old = nc_old(2)
            case(3)
                n1_old = nc_old(3)
            case default
                call this%err%set(&
                    code       = 100,&
                    severity   = 1,&
                    category   = 'forcad_nurbs_volume',&
                    message    = 'Invalid direction for inserting knots.',&
                    location   = 'insert_knots',&
                    suggestion = 'Use dir=1 or dir=2 or dir=3 to specify the direction.')
                return
            end select
            allocate(A1(n1_old,n1_old), source=0.0_rk)
            do concurrent (j = 1:n1_old)
                A1(j,j) = 1.0_rk
            end do
        end if

        if (dir == 1) then

            if (this%is_rational()) then
                allocate(Xcw(size(this%Xc,1), dim+1))
                do concurrent (j = 1:size(this%Xc,1))
                    Xcw(j,1:dim) = this%Xc(j,1:dim) * this%Wc(j)
                end do
                Xcw(:,dim+1) = this%Wc(:)
                Xcw = reshape(Xcw, [this%nc(1), this%nc(2)*this%nc(3)*(dim+1)])

                do i = 1, size(Xth)
                    k = findspan(this%nc(1)-1,this%degree(1),Xth(i),this%knot1)
                    if (abs(this%knot1(k+1) - Xth(i)) < 2.0_rk*epsilon(0.0_rk)) then
                        s = compute_multiplicity(this%knot1, Xth(i))
                    else
                        s = 0
                    end if

                    if (present(B) .or. present(Bs)) then
                        call insert_knot_A_5_1(this%degree(1), this%knot1, Xcw, Xth(i), k, s, r(i), n_new, knot_new, Xcw_new, A_re_loc)
                        A1 = matmul(A_re_loc, A1)
                    else
                        call insert_knot_A_5_1(this%degree(1), this%knot1, Xcw, Xth(i), k, s, r(i), n_new, knot_new, Xcw_new)
                    end if

                    call move_alloc(Xcw_new, Xcw)

                    H = reshape(Xcw, [ (n_new+1)*this%nc(2)*this%nc(3), dim+1 ])

                    associate (C => H(:,1:dim), W => H(:,dim+1))
                        do j=1,dim
                            C(:,j) = C(:,j) / H(:,dim+1)
                        end do
                        call this%set(knot1=knot_new, knot2=this%get_knot(2), knot3=this%get_knot(3), Xc=C, Wc=W)
                    end associate
                    deallocate(H)
                end do

            else
                Xc = reshape(this%Xc, [this%nc(1), this%nc(2)*this%nc(3)*dim])

                do i = 1, size(Xth)
                    k = findspan(this%nc(1)-1,this%degree(1),Xth(i),this%knot1)
                    if (abs(this%knot1(k+1) - Xth(i)) < 2.0_rk*epsilon(0.0_rk)) then
                        s = compute_multiplicity(this%knot1, Xth(i))
                    else
                        s = 0
                    end if

                    if (present(B) .or. present(Bs)) then
                        call insert_knot_A_5_1(this%degree(1), this%knot1, Xc, Xth(i), k, s, r(i), n_new, knot_new, Xcw_new, A_re_loc)
                        A1 = matmul(A_re_loc, A1)
                    else
                        call insert_knot_A_5_1(this%degree(1), this%knot1, Xc, Xth(i), k, s, r(i), n_new, knot_new, Xcw_new)
                    end if

                    call move_alloc(Xcw_new, Xc)

                    H = reshape(Xc, [(n_new+1)*this%nc(2)*this%nc(3), dim])
                    call this%set(knot1=knot_new, knot2=this%get_knot(2), knot3=this%get_knot(3), Xc=H)
                    deallocate(H)
                end do
            end if

        elseif (dir == 2) then

            if (this%is_rational()) then
                allocate(Xcw(size(this%Xc,1), dim+1))
                do concurrent (j = 1:size(this%Xc,1))
                    Xcw(j,1:dim) = this%Xc(j,1:dim) * this%Wc(j)
                end do
                Xcw(:,dim+1) = this%Wc(:)
                Xc4 = reshape(Xcw, [this%nc(1), this%nc(2),this%nc(3), dim+1])
                Xc4 = reshape(Xc4, [this%nc(2), this%nc(1),this%nc(3), dim+1], order=[2,1,3,4])
                Xcw = reshape(Xc4, [this%nc(2), this%nc(1)*this%nc(3)*(dim+1)])

                do i = 1, size(Xth)
                    k = findspan(this%nc(2)-1,this%degree(2),Xth(i),this%knot2)
                    if (abs(this%knot2(k+1) - Xth(i)) < 2.0_rk*epsilon(0.0_rk)) then
                        s = compute_multiplicity(this%knot2, Xth(i))
                    else
                        s = 0
                    end if

                    if (present(B) .or. present(Bs)) then
                        call insert_knot_A_5_1(this%degree(2), this%knot2, Xcw, Xth(i), k, s, r(i), n_new, knot_new, Xcw_new, A_re_loc)
                        A1 = matmul(A_re_loc, A1)
                    else
                        call insert_knot_A_5_1(this%degree(2), this%knot2, Xcw, Xth(i), k, s, r(i), n_new, knot_new, Xcw_new)
                    end if

                    call move_alloc(Xcw_new, Xcw)

                    Xc4 = reshape(Xcw, [n_new+1, this%nc(1), this%nc(3), dim+1])
                    Xc4 = reshape(Xc4, [this%nc(1), n_new+1, this%nc(3), dim+1], order=[2,1,3,4])
                    H   = reshape(Xc4, [this%nc(1)*(n_new+1)*this%nc(3), dim+1])

                    associate (C => H(:,1:dim), W => H(:,dim+1))
                        do j=1,dim
                            C(:,j) = C(:,j) / W(:)
                        end do
                        call this%set(knot1=this%get_knot(1), knot2=knot_new, knot3=this%get_knot(3), Xc=C, Wc=W)
                    end associate
                    deallocate(H)
                end do

            else
                Xc4 = reshape(this%Xc, [this%nc(1),this%nc(2),this%nc(3), dim])
                Xc4 = reshape(Xc4,     [this%nc(2),this%nc(1),this%nc(3), dim], order=[2,1,3,4])
                Xc  = reshape(Xc4,     [this%nc(2), this%nc(1)*this%nc(3)*dim])

                do i = 1, size(Xth)
                    k = findspan(this%nc(2)-1,this%degree(2),Xth(i),this%knot2)
                    if (abs(this%knot2(k+1) - Xth(i)) < 2.0_rk*epsilon(0.0_rk)) then
                        s = compute_multiplicity(this%knot2, Xth(i))
                    else
                        s = 0
                    end if

                    if (present(B) .or. present(Bs)) then
                        call insert_knot_A_5_1(this%degree(2), this%knot2, Xc, Xth(i), k, s, r(i), n_new, knot_new, Xcw_new, A_re_loc)
                        A1 = matmul(A_re_loc, A1)
                    else
                        call insert_knot_A_5_1(this%degree(2), this%knot2, Xc, Xth(i), k, s, r(i), n_new, knot_new, Xcw_new)
                    end if

                    call move_alloc(Xcw_new, Xc)

                    Xc4 = reshape(Xc,  [n_new+1, this%nc(1), this%nc(3), dim])
                    Xc4 = reshape(Xc4, [this%nc(1), n_new+1, this%nc(3), dim], order=[2,1,3,4])
                    H   = reshape(Xc4, [this%nc(1)*(n_new+1)*this%nc(3), dim])
                    call this%set(knot1=this%get_knot(1), knot2=knot_new, knot3=this%get_knot(3), Xc=H)
                    deallocate(H)
                end do
            end if

        elseif (dir == 3) then

            if (this%is_rational()) then
                allocate(Xcw(size(this%Xc,1), dim+1))
                do concurrent (j = 1:size(this%Xc,1))
                    Xcw(j,1:dim) = this%Xc(j,1:dim) * this%Wc(j)
                end do
                Xcw(:,dim+1) = this%Wc(:)
                Xc4 = reshape(Xcw, [this%nc(1), this%nc(2),this%nc(3), dim+1])
                Xc4 = reshape(Xc4, [this%nc(3), this%nc(2),this%nc(1), dim+1], order=[3,2,1,4])
                Xcw = reshape(Xc4, [this%nc(3), this%nc(2)*this%nc(1)*(dim+1)])

                do i = 1, size(Xth)
                    k = findspan(this%nc(3)-1,this%degree(3),Xth(i),this%knot3)
                    if (abs(this%knot3(k+1) - Xth(i)) < 2.0_rk*epsilon(0.0_rk)) then
                        s = compute_multiplicity(this%knot3, Xth(i))
                    else
                        s = 0
                    end if

                    if (present(B) .or. present(Bs)) then
                        call insert_knot_A_5_1(this%degree(3), this%knot3, Xcw, Xth(i), k, s, r(i), n_new, knot_new, Xcw_new, A_re_loc)
                        A1 = matmul(A_re_loc, A1)
                    else
                        call insert_knot_A_5_1(this%degree(3), this%knot3, Xcw, Xth(i), k, s, r(i), n_new, knot_new, Xcw_new)
                    end if

                    call move_alloc(Xcw_new, Xcw)

                    Xc4 = reshape(Xcw, [n_new+1, this%nc(2), this%nc(1), dim+1])
                    Xc4 = reshape(Xc4, [this%nc(1), this%nc(2), n_new+1, dim+1], order=[3,2,1,4])
                    H   = reshape(Xc4, [this%nc(1)*this%nc(2)*(n_new+1), dim+1])

                    associate (C => H(:,1:dim), W => H(:,dim+1))
                        do j=1,dim
                            C(:,j) = C(:,j) / W(:)
                        end do
                        call this%set(knot1=this%get_knot(1), knot2=this%get_knot(2), knot3=knot_new, Xc=C, Wc=W)
                    end associate
                    deallocate(H)
                end do

            else
                Xc4 = reshape(this%Xc, [this%nc(1),this%nc(2),this%nc(3), dim])
                Xc4 = reshape(Xc4,     [this%nc(3),this%nc(2),this%nc(1), dim], order=[3,2,1,4])
                Xc  = reshape(Xc4,     [this%nc(3), this%nc(2)*this%nc(1)*dim])

                do i = 1, size(Xth)
                    k = findspan(this%nc(3)-1,this%degree(3),Xth(i),this%knot3)
                    if (abs(this%knot3(k+1) - Xth(i)) < 2.0_rk*epsilon(0.0_rk)) then
                        s = compute_multiplicity(this%knot3, Xth(i))
                    else
                        s = 0
                    end if

                    if (present(B) .or. present(Bs)) then
                        call insert_knot_A_5_1(this%degree(3), this%knot3, Xc, Xth(i), k, s, r(i), n_new, knot_new, Xcw_new, A_re_loc)
                        A1 = matmul(A_re_loc, A1)
                    else
                        call insert_knot_A_5_1(this%degree(3), this%knot3, Xc, Xth(i), k, s, r(i), n_new, knot_new, Xcw_new)
                    end if

                    call move_alloc(Xcw_new, Xc)

                    Xc4 = reshape(Xc,  [n_new+1, this%nc(2), this%nc(1), dim])
                    Xc4 = reshape(Xc4, [this%nc(1), this%nc(2), n_new+1, dim], order=[3,2,1,4])
                    H   = reshape(Xc4, [this%nc(1)*this%nc(2)*(n_new+1), dim])
                    call this%set(knot1=this%get_knot(1), knot2=this%get_knot(2), knot3=knot_new, Xc=H)
                    deallocate(H)
                end do
            end if

        else
            call this%err%set(&
                code       = 100,&
                severity   = 1,&
                category   = 'forcad_nurbs_volume',&
                message    = 'Invalid direction for inserting knots.',&
                location   = 'insert_knots',&
                suggestion = 'Use dir=1 or dir=2 or dir=3 to specify the direction.')
            return
        end if

        if (present(B) .or. present(Bs)) then
            block
                real(rk), allocatable :: S_loc(:,:)
                integer :: nc1, nc2, nc3, n1_new
                integer :: i1, j1, i2, i3, i2_old, i2_new, ii, i3_old, i3_new
                integer :: mS, nS, c

                nc1 = this%nc(1)
                nc2 = this%nc(2)
                nc3 = this%nc(3)

                select case(dir)
                case (1)
                    n1_new = this%nc(1)
                    mS = n1_new*nc2*nc3
                    nS = nc_old(1)*nc2*nc3
                    allocate(S_loc(mS,nS), source=0.0_rk)
                    if (this%is_rational()) then
                        do concurrent (i3=0:nc3-1, i2=0:nc2-1, j1=1:nc_old(1), i1=1:n1_new, A1(i1,j1) /= 0.0_rk)
                            S_loc((i3*nc2 + i2)*n1_new    + i1,(i3*nc2 + i2)*nc_old(1) + j1) = A1(i1,j1) * Wc_old((i3*nc2 + i2)*nc_old(1) + j1) / this%Wc((i3*nc2 + i2)*n1_new    + i1)
                        end do
                    else
                        do concurrent (i3=0:nc3-1, i2=0:nc2-1, j1=1:nc_old(1), i1=1:n1_new, A1(i1,j1) /= 0.0_rk)
                            S_loc((i3*nc2 + i2)*n1_new    + i1,(i3*nc2 + i2)*nc_old(1) + j1) = A1(i1,j1)
                        end do
                    end if

                case (2)
                    n1_new = this%nc(2)
                    mS = n1_new*nc1*nc3
                    nS = nc_old(2)*nc1*nc3
                    allocate(S_loc(mS,nS), source=0.0_rk)
                    if (this%is_rational()) then
                        do concurrent (i3=0:nc3-1, i2_old=1:nc_old(2), i2_new=1:n1_new, ii=0:nc1-1, A1(i2_new,i2_old) /= 0.0_rk)
                            S_loc(ii + 1 + (i2_new-1)*nc1 + i3*nc1*n1_new,ii + 1 + (i2_old-1)*nc1 + i3*nc1*nc_old(2)) = A1(i2_new,i2_old) * Wc_old(ii + 1 + (i2_old-1)*nc1 + i3*nc1*nc_old(2)) / this%Wc(ii + 1 + (i2_new-1)*nc1 + i3*nc1*n1_new)
                        end do
                    else
                        do concurrent (i3=0:nc3-1, i2_old=1:nc_old(2), i2_new=1:n1_new, ii=0:nc1-1, A1(i2_new,i2_old) /= 0.0_rk)
                            S_loc(ii + 1 + (i2_new-1)*nc1 + i3*nc1*n1_new,ii + 1 + (i2_old-1)*nc1 + i3*nc1*nc_old(2)) = A1(i2_new,i2_old)
                        end do
                    end if

                case (3)
                    n1_new = this%nc(3)
                    mS = n1_new*nc1*nc2
                    nS = nc_old(3)*nc1*nc2
                    allocate(S_loc(mS,nS), source=0.0_rk)
                    if (this%is_rational()) then
                        do concurrent (i3_old=1:nc_old(3), i3_new=1:n1_new, i2=0:nc2-1, ii=0:nc1-1, A1(i3_new,i3_old) /= 0.0_rk)
                            S_loc(ii + 1 + i2*nc1 + (i3_new-1)*nc1*nc2,ii + 1 + i2*nc1 + (i3_old-1)*nc1*nc2) = A1(i3_new,i3_old) * Wc_old(ii + 1 + i2*nc1 + (i3_old-1)*nc1*nc2) / this%Wc(ii + 1 + i2*nc1 + (i3_new-1)*nc1*nc2)
                        end do
                    else
                        do concurrent (i3_old=1:nc_old(3), i3_new=1:n1_new, i2=0:nc2-1, ii=0:nc1-1, A1(i3_new,i3_old) /= 0.0_rk)
                            S_loc(ii + 1 + i2*nc1 + (i3_new-1)*nc1*nc2,ii + 1 + i2*nc1 + (i3_old-1)*nc1*nc2) = A1(i3_new,i3_old)
                        end do
                    end if
                end select

                if (present(B)) then
                    allocate(B(mS*dim, nS*dim), source=0.0_rk)
                    do c = 1, dim
                        B(c:mS*dim:dim, c:nS*dim:dim) = S_loc
                    end do
                end if

                if (present(Bs)) then
                    call move_alloc(S_loc, Bs)
                else
                    deallocate(S_loc)
                end if
            end block
        end if
    end subroutine
    !===============================================================================


    !===============================================================================
    !> author:  Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure subroutine elevate_degree(this, dir, t, B, Bs)
        class(nurbs_volume),             intent(inout) :: this
        integer,                         intent(in)    :: dir
        integer,                         intent(in)    :: t
        real(rk), allocatable, optional, intent(out)   :: B(:,:)
        real(rk), allocatable, optional, intent(out)   :: Bs(:,:)

        integer :: n1_new
        real(rk), allocatable :: knot_new(:)
        real(rk), allocatable :: Xc(:,:), Xcw(:,:), Xcw_new(:,:), H(:,:), Xc4(:,:,:,:)
        real(rk), allocatable :: Tdir(:,:)
        integer :: nc_old(3), dim, ncp_old, mS, nS, c
        real(rk), allocatable :: Wc_old(:), S_loc(:,:)
        integer :: i1, j1, i2, i3, i2_old, i2_new, ii, i3_old, i3_new

        if (.not. this%err%ok) return

        dim = size(this%Xc,2)

        if (present(B) .or. present(Bs)) then
            nc_old  = this%nc
            ncp_old = size(this%Xc,1)
            if (this%is_rational()) then
                allocate(Wc_old(ncp_old)); Wc_old = this%Wc
            end if
        end if


        select case(dir)
        case (1)
            if (this%is_rational()) then
                allocate(Xcw(size(this%Xc,1), dim+1))
                do concurrent (i1=1:size(this%Xc,1))
                    Xcw(i1,1:dim) = this%Xc(i1,1:dim) * this%Wc(i1)
                end do
                Xcw(:,dim+1) = this%Wc(:)
                Xcw = reshape(Xcw, [ this%nc(1), this%nc(2)*this%nc(3)*(dim+1) ])

                if (present(B) .or. present(Bs)) then
                    call elevate_degree_A_5_9(t, this%knot1, this%degree(1), Xcw, n1_new, knot_new, Xcw_new, Tdir)
                else
                    call elevate_degree_A_5_9(t, this%knot1, this%degree(1), Xcw, n1_new, knot_new, Xcw_new)
                end if

                H = reshape(Xcw_new, [ n1_new*this%nc(2)*this%nc(3), dim+1 ])
                associate (C => H(:,1:dim), W => H(:,dim+1))
                    do i1=1,dim
                        C(:,i1) = C(:,i1) / W(:)
                    end do
                    call this%set(knot1=knot_new, knot2=this%get_knot(2), knot3=this%get_knot(3), Xc=C, Wc=W)
                end associate
                deallocate(H, Xcw, Xcw_new)

            else
                Xc = reshape(this%Xc, [ this%nc(1), this%nc(2)*this%nc(3)*dim ])
                if (present(B) .or. present(Bs)) then
                    call elevate_degree_A_5_9(t, this%knot1, this%degree(1), Xc, n1_new, knot_new, Xcw_new, Tdir)
                else
                    call elevate_degree_A_5_9(t, this%knot1, this%degree(1), Xc, n1_new, knot_new, Xcw_new)
                end if
                H = reshape(Xcw_new, [ n1_new*this%nc(2)*this%nc(3), dim ])
                call this%set(knot1=knot_new, knot2=this%get_knot(2), knot3=this%get_knot(3), Xc=H)
                deallocate(H, Xc, Xcw_new)
            end if

        case (2)
            if (this%is_rational()) then
                allocate(Xcw(size(this%Xc,1), dim+1))
                do concurrent (i1=1:size(this%Xc,1))
                    Xcw(i1,1:dim) = this%Xc(i1,1:dim) * this%Wc(i1)
                end do
                Xcw(:,dim+1) = this%Wc(:)

                Xc4 = reshape(Xcw, [ this%nc(1), this%nc(2), this%nc(3), dim+1 ])
                Xc4 = reshape(Xc4, [ this%nc(2), this%nc(1), this%nc(3), dim+1 ], order=[2,1,3,4])
                Xcw = reshape(Xc4, [ this%nc(2), this%nc(1)*this%nc(3)*(dim+1) ])

                if (present(B) .or. present(Bs)) then
                    call elevate_degree_A_5_9(t, this%knot2, this%degree(2), Xcw, n1_new, knot_new, Xcw_new, Tdir)
                else
                    call elevate_degree_A_5_9(t, this%knot2, this%degree(2), Xcw, n1_new, knot_new, Xcw_new)
                end if

                Xc4 = reshape(Xcw_new, [ n1_new, this%nc(1), this%nc(3), dim+1 ])
                Xc4 = reshape(Xc4,     [ this%nc(1), n1_new, this%nc(3), dim+1 ], order=[2,1,3,4])
                H   = reshape(Xc4,     [ this%nc(1)*n1_new*this%nc(3), dim+1 ])
                associate (C => H(:,1:dim), W => H(:,dim+1))
                    do i1=1,dim
                        C(:,i1) = C(:,i1) / W(:)
                    end do
                    call this%set(knot1=this%get_knot(1), knot2=knot_new, knot3=this%get_knot(3), Xc=C, Wc=W)
                end associate
                deallocate(H, Xcw, Xcw_new, Xc4)

            else
                Xc4 = reshape(this%Xc, [ this%nc(1), this%nc(2), this%nc(3), dim ])
                Xc4 = reshape(Xc4,     [ this%nc(2), this%nc(1), this%nc(3), dim ], order=[2,1,3,4])
                Xc  = reshape(Xc4,     [ this%nc(2), this%nc(1)*this%nc(3)*dim ])

                if (present(B) .or. present(Bs)) then
                    call elevate_degree_A_5_9(t, this%knot2, this%degree(2), Xc, n1_new, knot_new, Xcw_new, Tdir)
                else
                    call elevate_degree_A_5_9(t, this%knot2, this%degree(2), Xc, n1_new, knot_new, Xcw_new)
                end if

                Xc4 = reshape(Xcw_new, [ n1_new, this%nc(1), this%nc(3), dim ])
                Xc4 = reshape(Xc4,     [ this%nc(1), n1_new, this%nc(3), dim ], order=[2,1,3,4])
                H   = reshape(Xc4,     [ this%nc(1)*n1_new*this%nc(3), dim ])
                call this%set(knot1=this%get_knot(1), knot2=knot_new, knot3=this%get_knot(3), Xc=H)
                deallocate(H, Xc, Xc4, Xcw_new)
            end if

        case (3)
            if (this%is_rational()) then
                allocate(Xcw(size(this%Xc,1), dim+1))
                do concurrent (i1=1:size(this%Xc,1))
                    Xcw(i1,1:dim) = this%Xc(i1,1:dim) * this%Wc(i1)
                end do
                Xcw(:,dim+1) = this%Wc(:)

                Xc4 = reshape(Xcw, [ this%nc(1), this%nc(2), this%nc(3), dim+1 ])
                Xc4 = reshape(Xc4, [ this%nc(3), this%nc(2), this%nc(1), dim+1 ], order=[3,2,1,4])
                Xcw = reshape(Xc4, [ this%nc(3), this%nc(2)*this%nc(1)*(dim+1) ])

                if (present(B) .or. present(Bs)) then
                    call elevate_degree_A_5_9(t, this%knot3, this%degree(3), Xcw, n1_new, knot_new, Xcw_new, Tdir)
                else
                    call elevate_degree_A_5_9(t, this%knot3, this%degree(3), Xcw, n1_new, knot_new, Xcw_new)
                end if

                Xc4 = reshape(Xcw_new, [ n1_new, this%nc(2), this%nc(1), dim+1 ])
                Xc4 = reshape(Xc4,     [ this%nc(1), this%nc(2), n1_new, dim+1 ], order=[3,2,1,4])
                H   = reshape(Xc4,     [ this%nc(1)*this%nc(2)*n1_new, dim+1 ])
                associate (C => H(:,1:dim), W => H(:,dim+1))
                    do i1=1,dim
                        C(:,i1) = C(:,i1) / W(:)
                    end do
                    call this%set(knot1=this%get_knot(1), knot2=this%get_knot(2), knot3=knot_new, Xc=C, Wc=W)
                end associate
                deallocate(H, Xcw, Xcw_new, Xc4)

            else
                Xc4 = reshape(this%Xc, [ this%nc(1), this%nc(2), this%nc(3), dim ])
                Xc4 = reshape(Xc4,     [ this%nc(3), this%nc(2), this%nc(1), dim ], order=[3,2,1,4])
                Xc  = reshape(Xc4,     [ this%nc(3), this%nc(2)*this%nc(1)*dim ])

                if (present(B) .or. present(Bs)) then
                    call elevate_degree_A_5_9(t, this%knot3, this%degree(3), Xc, n1_new, knot_new, Xcw_new, Tdir)
                else
                    call elevate_degree_A_5_9(t, this%knot3, this%degree(3), Xc, n1_new, knot_new, Xcw_new)
                end if

                Xc4 = reshape(Xcw_new, [ n1_new, this%nc(2), this%nc(1), dim ])
                Xc4 = reshape(Xc4,     [ this%nc(1), this%nc(2), n1_new, dim ], order=[3,2,1,4])
                H   = reshape(Xc4,     [ this%nc(1)*this%nc(2)*n1_new, dim ])
                call this%set(knot1=this%get_knot(1), knot2=this%get_knot(2), knot3=knot_new, Xc=H)
                deallocate(H, Xc, Xc4, Xcw_new)
            end if

        case default
            call this%err%set(&
                code       = 100,&
                severity   = 1,&
                category   = 'forcad_nurbs_volume',&
                message    = 'Invalid direction for elevating degree.',&
                location   = 'elevate_degree',&
                suggestion = 'Use dir=1 or dir=2 or dir=3 to specify the direction.')
            return
        end select

        if (present(B) .or. present(Bs)) then
            select case(dir)
            case(1)
                mS = this%nc(1)*this%nc(2)*this%nc(3)
                nS = nc_old(1) *this%nc(2)*this%nc(3)
                allocate(S_loc(mS,nS), source=0.0_rk)

                if (this%is_rational()) then
                    do concurrent (i3=0:this%nc(3)-1, i2=0:this%nc(2)-1, j1=1:nc_old(1), i1=1:this%nc(1), Tdir(i1,j1) /= 0.0_rk)
                        S_loc((i3*this%nc(2) + i2)*this%nc(1)    + i1,(i3*this%nc(2) + i2)*nc_old(1)     + j1) = Tdir(i1,j1) * Wc_old((i3*this%nc(2) + i2)*nc_old(1)     + j1) / this%Wc((i3*this%nc(2) + i2)*this%nc(1)    + i1)
                    end do
                else
                    do concurrent (i3=0:this%nc(3)-1, i2=0:this%nc(2)-1, j1=1:nc_old(1), i1=1:this%nc(1), Tdir(i1,j1) /= 0.0_rk)
                        S_loc((i3*this%nc(2) + i2)*this%nc(1)    + i1,(i3*this%nc(2) + i2)*nc_old(1)     + j1) = Tdir(i1,j1)
                    end do
                end if

            case(2)
                mS = this%nc(2)*this%nc(1)*this%nc(3)
                nS = nc_old(2) *this%nc(1)*this%nc(3)
                allocate(S_loc(mS,nS), source=0.0_rk)

                if (this%is_rational()) then
                    do concurrent (i3=0:this%nc(3)-1, i2_old=1:nc_old(2), i2_new=1:this%nc(2), ii=0:this%nc(1)-1, Tdir(i2_new,i2_old) /= 0.0_rk)
                        S_loc(ii + 1 + (i2_new-1)*this%nc(1) + i3*this%nc(1)*this%nc(2),ii + 1 + (i2_old-1)*this%nc(1) + i3*this%nc(1)*nc_old(2)) = Tdir(i2_new,i2_old) * Wc_old(ii + 1 + (i2_old-1)*this%nc(1) + i3*this%nc(1)*nc_old(2)) / this%Wc(ii + 1 + (i2_new-1)*this%nc(1) + i3*this%nc(1)*this%nc(2))
                    end do
                else
                    do concurrent (i3=0:this%nc(3)-1, i2_old=1:nc_old(2), i2_new=1:this%nc(2), ii=0:this%nc(1)-1, Tdir(i2_new,i2_old) /= 0.0_rk)
                        S_loc(ii + 1 + (i2_new-1)*this%nc(1) + i3*this%nc(1)*this%nc(2),ii + 1 + (i2_old-1)*this%nc(1) + i3*this%nc(1)*nc_old(2)) = Tdir(i2_new,i2_old)
                    end do
                end if

            case(3)
                mS = this%nc(3)*this%nc(1)*this%nc(2)
                nS = nc_old(3) *this%nc(1)*this%nc(2)
                allocate(S_loc(mS,nS), source=0.0_rk)

                if (this%is_rational()) then
                    do concurrent (i3_old=1:nc_old(3), i3_new=1:this%nc(3), i2=0:this%nc(2)-1, ii=0:this%nc(1)-1, Tdir(i3_new,i3_old) /= 0.0_rk)
                        S_loc(ii + 1 + i2*this%nc(1) + (i3_new-1)*this%nc(1)*this%nc(2),ii + 1 + i2*this%nc(1) + (i3_old-1)*this%nc(1)*this%nc(2)) = Tdir(i3_new,i3_old) * Wc_old(ii + 1 + i2*this%nc(1) + (i3_old-1)*this%nc(1)*this%nc(2)) / this%Wc(ii + 1 + i2*this%nc(1) + (i3_new-1)*this%nc(1)*this%nc(2))
                    end do
                else
                    do concurrent (i3_old=1:nc_old(3), i3_new=1:this%nc(3), i2=0:this%nc(2)-1, ii=0:this%nc(1)-1, Tdir(i3_new,i3_old) /= 0.0_rk)
                        S_loc(ii + 1 + i2*this%nc(1) + (i3_new-1)*this%nc(1)*this%nc(2),ii + 1 + i2*this%nc(1) + (i3_old-1)*this%nc(1)*this%nc(2)) = Tdir(i3_new,i3_old)
                    end do
                end if
            end select

            if (present(Bs)) then
                call move_alloc(S_loc, Bs)
            end if

            if (present(B)) then
                if (.not. present(Bs)) then
                    allocate(B(mS*dim, nS*dim), source=0.0_rk)
                    do c = 1, dim
                        B(c:mS*dim:dim, c:nS*dim:dim) = S_loc
                    end do
                    deallocate(S_loc)
                else
                    allocate(B(mS*dim, nS*dim), source=0.0_rk)
                    do c = 1, dim
                        B(c:mS*dim:dim, c:nS*dim:dim) = Bs
                    end do
                end if
            end if

            if (allocated(Wc_old)) deallocate(Wc_old)
        end if
    end subroutine
    !===============================================================================


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

        if (.not. this%err%ok) return

        r = .false.
        if (allocated(this%Wc)) then
            ! if (any(this%Wc /= this%Wc(1))) then
            if (any(abs(this%Wc - this%Wc(1)) > 2.0_rk*epsilon(0.0_rk))) then
                r = .true.
            end if
        end if
    end function
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure subroutine set_elem_Xc_vis(this, elemConn)
        class(nurbs_volume), intent(inout) :: this
        integer, intent(in), contiguous :: elemConn(:,:)

        if (.not. this%err%ok) return

        if (allocated(this%elemConn_Xc_vis)) then
            if (size(this%elemConn_Xc_vis,1) /= size(elemConn,1) .or. size(this%elemConn_Xc_vis,2) /= size(elemConn,2)) then
                deallocate(this%elemConn_Xc_vis)
            end if
        end if
        this%elemConn_Xc_vis = elemConn
    end subroutine
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure subroutine set_elem_Xg_vis(this, elemConn)
        class(nurbs_volume), intent(inout) :: this
        integer, intent(in), contiguous :: elemConn(:,:)

        if (.not. this%err%ok) return

        if (allocated(this%elemConn_Xg_vis)) then
            if (size(this%elemConn_Xg_vis,1) /= size(elemConn,1) .or. size(this%elemConn_Xg_vis,2) /= size(elemConn,2)) then
                deallocate(this%elemConn_Xg_vis)
            end if
        end if
        this%elemConn_Xg_vis = elemConn
    end subroutine
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure subroutine set_elem(this, elemConn)
        class(nurbs_volume), intent(inout) :: this
        integer, intent(in), contiguous :: elemConn(:,:)

        if (.not. this%err%ok) return

        if (allocated(this%elemConn)) then
            if (size(this%elemConn,1) /= size(elemConn,1) .or. size(this%elemConn,2) /= size(elemConn,2)) then
                deallocate(this%elemConn)
            end if
        end if
        this%elemConn = elemConn
    end subroutine
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure function get_elem_Xc_vis(this) result(elemConn)
        class(nurbs_volume), intent(in) :: this
        integer, allocatable :: elemConn(:,:)

        if (.not. this%err%ok) return

        elemConn = this%elemConn_Xc_vis
    end function
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure function get_elem_Xg_vis(this) result(elemConn)
        class(nurbs_volume), intent(in) :: this
        integer, allocatable :: elemConn(:,:)

        if (.not. this%err%ok) return

        elemConn = this%elemConn_Xg_vis
    end function
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure function get_elem(this) result(elemConn)
        class(nurbs_volume), intent(in) :: this
        integer, allocatable :: elemConn(:,:)

        if (.not. this%err%ok) return

        elemConn = this%elemConn
    end function
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure subroutine set_hexahedron(this, L, nc, Wc)
        class(nurbs_volume), intent(inout) :: this
        real(rk), intent(in), contiguous :: L(:)
        integer, intent(in), contiguous :: nc(:)
        real(rk), intent(in), contiguous, optional :: Wc(:)

        if (.not. this%err%ok) return

        if (present(Wc)) then
            call this%set(nc, hexahedron_Xc(L, nc), Wc)
        else
            call this%set(nc, hexahedron_Xc(L, nc))
        end if
    end subroutine
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure subroutine put_to_nurbs(this, X, elemConn)
        class(nurbs_volume), intent(inout) :: this
        real(rk), intent(in), contiguous :: X(:,:)
        integer, intent(in), contiguous, optional :: elemConn(:,:)
        real(rk), allocatable :: Xt(:,:)
        real(rk) :: min_X1, max_X1, min_X2, max_X2, min_X3, max_X3

        if (.not. this%err%ok) return

        ! Normalize the X coordinates to the range of knot vectors
        allocate(Xt(size(X,1), size(X,2)))
        min_X1 = minval(X(:,1))
        max_X1 = maxval(X(:,1))

        min_X2 = minval(X(:,2))
        max_X2 = maxval(X(:,2))

        min_X3 = minval(X(:,3))
        max_X3 = maxval(X(:,3))

        Xt(:,1) = (X(:,1) - min_X1) / (max_X1 - min_X1) * (this%knot1(size(this%knot1)) - this%knot1(1)) + this%knot1(1)
        Xt(:,2) = (X(:,2) - min_X2) / (max_X2 - min_X2) * (this%knot2(size(this%knot2)) - this%knot2(1)) + this%knot2(1)
        Xt(:,3) = (X(:,3) - min_X3) / (max_X3 - min_X3) * (this%knot3(size(this%knot3)) - this%knot3(1)) + this%knot3(1)

        if (this%is_rational()) then
            this%Xg = compute_Xg(Xt=Xt, knot1=this%knot1, knot2=this%knot2, knot3=this%knot3, degree=this%degree, nc=this%nc, Xc=this%Xc, Wc=this%Wc)
        else
            this%Xg = compute_Xg(Xt=Xt, knot1=this%knot1, knot2=this%knot2, knot3=this%knot3, degree=this%degree, nc=this%nc, Xc=this%Xc)
        end if

        if (present(elemConn)) call this%set_elem_Xg_vis(elemConn)
    end subroutine
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure subroutine remove_knots(this, dir ,Xth,r)
        class(nurbs_volume), intent(inout) :: this
        integer, intent(in) :: dir
        real(rk), intent(in), contiguous :: Xth(:)
        integer, intent(in), contiguous :: r(:)
        integer :: k, i, s, d, j, nc_new, t
        real(rk), allocatable :: Xc(:,:), Xcw(:,:), Xcw_new(:,:), Xc_new(:,:), Wc_new(:), knot_new(:)
        real(rk), allocatable :: Xc4(:,:,:,:)

        if (.not. this%err%ok) return

        if (dir == 1) then ! direction 1

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

                do i = 1, size(Xth)
                    k = findspan(this%nc(1)-1,this%degree(1),Xth(i),this%knot1)
                    ! if (this%knot1(k+1) == Xth(i)) then
                    if (abs(this%knot1(k+1) - Xth(i)) < 2.0_rk*epsilon(0.0_rk)) then
                        s = compute_multiplicity(this%knot1, Xth(i))
                    else
                        s = 0
                    end if
                    k = k + 1

                    d = size(this%Xc,2)
                    allocate(Xcw(size(this%Xc,1),d+1))
                    do j = 1, size(this%Xc,1)
                        Xcw(j,1:d) = this%Xc(j,1:d)*this%Wc(j)
                    end do
                    Xcw(:,d+1) = this%Wc(:)

                    Xcw = reshape(Xcw,[this%nc(1),this%nc(2)*this%nc(3)*(d+1)],order=[1,2])

                    call remove_knots_A_5_8(&
                        this%degree(1),&
                        this%knot1,&
                        Xcw,&
                        Xth(i),&
                        k,&
                        s,&
                        r(i),&
                        t,&
                        knot_new,&
                        Xcw_new)

                    if (allocated(Xcw)) deallocate(Xcw)

                    if (t == 0) then
                        ! no change
                    else
                        nc_new = size(Xcw_new,1)
                        Xcw_new = reshape(Xcw_new,[(nc_new)*this%nc(2)*this%nc(3),d+1],order=[1,2])

                        allocate(Xc_new(1:(nc_new)*this%nc(2)*this%nc(3),1:d))
                        allocate(Wc_new(1:(nc_new)*this%nc(2)*this%nc(3)))
                        do j = 1, (nc_new)*this%nc(2)*this%nc(3)
                            Xc_new(j,1:d) = Xcw_new(j,1:d)/Xcw_new(j,d+1)
                        end do
                        Wc_new(:) = Xcw_new(:,d+1)

                        call this%set(knot1=knot_new, knot2=this%get_knot(2), knot3=this%get_knot(3), Xc=Xc_new, Wc=Wc_new)
                        deallocate(Xcw_new, Xc_new, Wc_new)
                    end if
                end do


            else ! B-Spline

                do i = 1, size(Xth)
                    k = findspan(this%nc(1)-1,this%degree(1),Xth(i),this%knot1)
                    ! if (this%knot1(k+1) == Xth(i)) then
                    if (abs(this%knot1(k+1) - Xth(i)) < 2.0_rk*epsilon(0.0_rk)) then
                        s = compute_multiplicity(this%knot1, Xth(i))
                    else
                        s = 0
                    end if
                    k = k + 1

                    d = size(this%Xc,2)

                    Xc = reshape(this%Xc,[this%nc(1),this%nc(2)*this%nc(3)*d])

                    call remove_knots_A_5_8(&
                        this%degree(1),&
                        this%knot1,&
                        Xc,&
                        Xth(i),&
                        k,&
                        s,&
                        r(i),&
                        t,&
                        knot_new,&
                        Xc_new)

                    if (allocated(Xc)) deallocate(Xc)

                    if (t == 0) then
                        ! no change
                    else
                        nc_new = size(Xcw_new,1)
                        Xc_new = reshape(Xc_new,[(nc_new)*this%nc(2)*this%nc(3),d])

                        call this%set(knot1=knot_new, knot2=this%get_knot(2), knot3=this%get_knot(3), Xc=Xc_new)
                    end if
                end do

            end if

        elseif (dir == 2) then! direction 2

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

                do i = 1, size(Xth)
                    k = findspan(this%nc(2)-1,this%degree(2),Xth(i),this%knot2)
                    ! if (this%knot2(k+1) == Xth(i)) then
                    if (abs(this%knot2(k+1) - Xth(i)) < 2.0_rk*epsilon(0.0_rk)) then
                        s = compute_multiplicity(this%knot2, Xth(i))
                    else
                        s = 0
                    end if
                    k = k + 1

                    d = size(this%Xc,2)
                    allocate(Xcw(size(this%Xc,1),d+1))
                    do j = 1, size(this%Xc,1)
                        Xcw(j,1:d) = this%Xc(j,1:d)*this%Wc(j)
                    end do

                    Xcw(:,d+1) = this%Wc(:)

                    Xc4 = reshape(Xcw, [this%nc(1),this%nc(2),this%nc(3),d+1])
                    Xc4 = reshape(Xc4, [this%nc(2),this%nc(1),this%nc(3),d+1], order=[2,1,3,4])
                    Xcw = reshape(Xc4,[this%nc(2),this%nc(1)*this%nc(3)*(d+1)])

                    call remove_knots_A_5_8(&
                        this%degree(2),&
                        this%knot2,&
                        Xcw,&
                        Xth(i),&
                        k,&
                        s,&
                        r(i),&
                        t,&
                        knot_new,&
                        Xcw_new)

                    if (allocated(Xcw)) deallocate(Xcw)

                    if (t == 0) then
                        ! no change
                    else
                        nc_new = size(Xcw_new,1)

                        Xc4 = reshape(Xcw_new, [nc_new,this%nc(1),this%nc(3),d+1])
                        Xc4 = reshape(Xc4, [this%nc(1),nc_new,this%nc(3),d+1], order=[2,1,3,4])
                        Xcw_new = reshape(Xc4,[this%nc(1)*(nc_new)*this%nc(3),d+1])

                        allocate(Xc_new(1:this%nc(1)*(nc_new)*this%nc(3),1:d))
                        allocate(Wc_new(1:this%nc(1)*(nc_new)*this%nc(3)))
                        do j = 1, this%nc(1)*(nc_new)*this%nc(3)
                            Xc_new(j,1:d) = Xcw_new(j,1:d)/Xcw_new(j,d+1)
                        end do

                        Wc_new(:) = Xcw_new(:,d+1)

                        call this%set(knot1=this%get_knot(1), knot2=knot_new, knot3=this%get_knot(3), Xc=Xc_new, Wc=Wc_new)
                        deallocate(Xcw_new, Xc_new, Wc_new)
                    end if
                end do

            else ! B-Spline

                do i = 1, size(Xth)
                    k = findspan(this%nc(2)-1,this%degree(2),Xth(i),this%knot2)
                    ! if (this%knot2(k+1) == Xth(i)) then
                    if (abs(this%knot2(k+1) - Xth(i)) < 2.0_rk*epsilon(0.0_rk)) then
                        s = compute_multiplicity(this%knot2, Xth(i))
                    else
                        s = 0
                    end if
                    k = k + 1

                    d = size(this%Xc,2)

                    Xc4 = reshape(this%Xc, [this%nc(1),this%nc(2),this%nc(3),d])
                    Xc4 = reshape(Xc4, [this%nc(2),this%nc(1),this%nc(3),d], order=[2,1,3,4])
                    Xc = reshape(Xc4, [this%nc(2),this%nc(1)*this%nc(3)*d])

                    call remove_knots_A_5_8(&
                        this%degree(2),&
                        this%knot2,&
                        Xc,&
                        Xth(i),&
                        k,&
                        s,&
                        r(i),&
                        t,&
                        knot_new,&
                        Xc_new)

                    if (allocated(Xc)) deallocate(Xc)

                    if (t == 0) then
                        ! no change
                    else
                        nc_new = size(Xcw_new,1)

                        Xc4 = reshape(Xc_new, [nc_new,this%nc(1),this%nc(3),d])
                        Xc4 = reshape(Xc4, [this%nc(1),nc_new,this%nc(3),d], order=[2,1,3,4])
                        Xc_new = reshape(Xc4, [this%nc(1)*(nc_new)*this%nc(3),d])

                        call this%set(knot1=this%get_knot(1), knot2=knot_new, knot3=this%get_knot(3), Xc=Xc_new)
                    end if
                end do

            end if

        elseif (dir == 3) then! direction 3

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

                do i = 1, size(Xth)
                    k = findspan(this%nc(3)-1,this%degree(3),Xth(i),this%knot3)
                    ! if (this%knot3(k+1) == Xth(i)) then
                    if (abs(this%knot3(k+1) - Xth(i)) < 2.0_rk*epsilon(0.0_rk)) then
                        s = compute_multiplicity(this%knot3, Xth(i))
                    else
                        s = 0
                    end if
                    k = k + 1

                    d = size(this%Xc,2)
                    allocate(Xcw(size(this%Xc,1),d+1))
                    do j = 1, size(this%Xc,1)
                        Xcw(j,1:d) = this%Xc(j,1:d)*this%Wc(j)
                    end do
                    Xcw(:,d+1) = this%Wc(:)

                    Xc4 = reshape(Xcw, [this%nc(1),this%nc(2),this%nc(3),d+1])
                    Xc4 = reshape(Xc4, [this%nc(3),this%nc(2),this%nc(1),d+1], order=[3,2,1,4])
                    Xcw = reshape(Xc4, [this%nc(3),this%nc(2)*this%nc(1)*(d+1)])

                    call remove_knots_A_5_8(&
                        this%degree(3),&
                        this%knot3,&
                        Xcw,&
                        Xth(i),&
                        k,&
                        s,&
                        r(i),&
                        t,&
                        knot_new,&
                        Xcw_new)

                    if (allocated(Xcw)) deallocate(Xcw)

                    if (t == 0) then
                        ! no change
                    else
                        nc_new = size(Xcw_new,1)

                        Xc4 = reshape(Xcw_new, [nc_new,this%nc(2),this%nc(1),d+1])
                        Xc4 = reshape(Xc4, [this%nc(1),this%nc(2),nc_new,d+1], order=[3,2,1,4])
                        Xcw_new = reshape(Xc4, [this%nc(1)*this%nc(2)*(nc_new),d+1])

                        allocate(Xc_new(1:this%nc(1)*this%nc(2)*(nc_new),1:d))
                        allocate(Wc_new(1:this%nc(1)*this%nc(2)*(nc_new)))
                        do j = 1, this%nc(1)*this%nc(2)*(nc_new)
                            Xc_new(j,1:d) = Xcw_new(j,1:d)/Xcw_new(j,d+1)
                        end do
                        Wc_new(:) = Xcw_new(:,d+1)

                        call this%set(knot1=this%get_knot(1), knot2=this%get_knot(2), knot3=knot_new, Xc=Xc_new, Wc=Wc_new)
                        deallocate(Xcw_new, Xc_new, Wc_new)
                    end if
                end do

            else ! B-Spline

                do i = 1, size(Xth)
                    k = findspan(this%nc(3)-1,this%degree(3),Xth(i),this%knot3)
                    ! if (this%knot3(k+1) == Xth(i)) then
                    if (abs(this%knot3(k+1) - Xth(i)) < 2.0_rk*epsilon(0.0_rk)) then
                        s = compute_multiplicity(this%knot3, Xth(i))
                    else
                        s = 0
                    end if
                    k = k + 1

                    d = size(this%Xc,2)

                    Xc4 = reshape(this%Xc, [this%nc(1),this%nc(2),this%nc(3),d])
                    Xc4 = reshape(Xc4, [this%nc(3),this%nc(2),this%nc(1),d], order=[3,2,1,4])
                    Xc = reshape(Xc4, [this%nc(3),this%nc(2)*this%nc(1)*d])

                    call remove_knots_A_5_8(&
                        this%degree(3),&
                        this%knot3,&
                        Xc,&
                        Xth(i),&
                        k,&
                        s,&
                        r(i),&
                        t,&
                        knot_new,&
                        Xc_new)

                    if (allocated(Xc)) deallocate(Xc)

                    if (t == 0) then
                        ! no change
                    else
                        nc_new = size(Xcw_new,1)

                        Xc4 = reshape(Xc_new, [nc_new,this%nc(2),this%nc(1),d])
                        Xc4 = reshape(Xc4, [this%nc(1),this%nc(2),nc_new,d], order=[3,2,1,4])
                        Xc_new = reshape(Xc4, [this%nc(1)*this%nc(2)*(nc_new),d])

                        call this%set(knot1=this%get_knot(1), knot2=this%get_knot(2), knot3=knot_new, Xc=Xc_new)
                    end if
                end do


            end if

        else
            call this%err%set(&
                code       = 100,&
                severity   = 1,&
                category   = 'forcad_nurbs_volume',&
                message    = 'Invalid direction for removing knots.',&
                location   = 'remove_knots',&
                suggestion = 'Use dir=1 or dir=2 or dir=3 to specify the direction.')
            return
        end if

    end subroutine
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure function cmp_elem(this) result(elemConn)
        class(nurbs_volume), intent(in) :: this
        integer, allocatable :: elemConn(:,:)

        if (.not. this%err%ok) return

        call elemConn_Cn(this%nc(1), this%nc(2), this%nc(3),&
            this%degree(1),this%degree(2),this%degree(3),&
            unique(this%knot1),unique(this%knot2),unique(this%knot3),&
            this%get_multiplicity(1),this%get_multiplicity(2),this%get_multiplicity(3),&
            elemConn)
    end function
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure subroutine rotate_Xc(this, alpha, beta, theta)
        class(nurbs_volume), intent(inout) :: this
        real(rk), intent(in) :: alpha, beta, theta
        integer :: i

        if (.not. this%err%ok) return

        do i = 1, this%nc(1)*this%nc(2)*this%nc(3)
            this%Xc(i, :) = matmul(rotation(alpha,beta,theta), this%Xc(i, :))
        end do
    end subroutine
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure subroutine rotate_Xg(this, alpha, beta, theta)
        class(nurbs_volume), intent(inout) :: this
        real(rk), intent(in) :: alpha, beta, theta
        integer :: i

        if (.not. this%err%ok) return

        do i = 1, this%ng(1)*this%ng(2)*this%ng(3)
            this%Xg(i, :) = matmul(rotation(alpha,beta,theta), this%Xg(i, :))
        end do
    end subroutine
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure subroutine translate_Xc(this, vec)
        class(nurbs_volume), intent(inout) :: this
        real(rk), intent(in), contiguous :: vec(:)
        integer :: i

        if (.not. this%err%ok) return

        do i = 1, this%nc(1)*this%nc(2)*this%nc(3)
            this%Xc(i, :) = this%Xc(i, :) + vec
        end do
    end subroutine
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure subroutine translate_Xg(this, vec)
        class(nurbs_volume), intent(inout) :: this
        real(rk), intent(in), contiguous :: vec(:)
        integer :: i

        if (.not. this%err%ok) return

        do i = 1, this%ng(1)*this%ng(2)*this%ng(3)
            this%Xg(i, :) = this%Xg(i, :) + vec
        end do
    end subroutine
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    impure subroutine show(this, vtkfile_Xc, vtkfile_Xg, vtkfile_Xth_in_Xg)
        class(nurbs_volume), intent(inout) :: this
        character(len=*), intent(in) :: vtkfile_Xc, vtkfile_Xg
        character(len=*), intent(in), optional :: vtkfile_Xth_in_Xg
#ifndef NOSHOW_PYVISTA

        block
        character(len=3000) :: pyvista_script

        if (.not. this%err%ok) return

        pyvista_script = &
            "import pyvista as pv"//achar(10)//&
            "pv.global_theme.color = 'white'"//achar(10)//&
            "Xc = pv.read('"//trim(vtkfile_Xc)//"')"//achar(10)//&
            "Xg = pv.read('"//trim(vtkfile_Xg)//"')"//achar(10)//&
            "Xg = Xg.clean(tolerance=1e-12)"//achar(10)//&
            "p = pv.Plotter(lighting='light kit')"//achar(10)//&
            "actor_Xcp = p.add_mesh("//achar(10)//&
            "    Xc,"//achar(10)//&
            "    style='points',"//achar(10)//&
            "    point_size=10,"//achar(10)//&
            "    color='red',"//achar(10)//&
            "    render_points_as_spheres=True,"//achar(10)//&
            "    opacity=0.5,"//achar(10)//&
            ")"//achar(10)//&
            "actor_Xcw = p.add_mesh("//achar(10)//&
            "    Xc,"//achar(10)//&
            "    show_edges=True,"//achar(10)//&
            "    color='yellow',"//achar(10)//&
            "    line_width=3,"//achar(10)//&
            "    style='wireframe',"//achar(10)//&
            "    opacity=0.2"//achar(10)//&
            ")"//achar(10)//&
            "actor_Xg = p.add_mesh("//achar(10)//&
            "    Xg,"//achar(10)//&
            "    show_edges=False,"//achar(10)//&
            "    color='cyan',"//achar(10)//&
            "    line_width=1,"//achar(10)//&
            "    metallic=0.6,"//achar(10)//&
            "    pbr=True,"//achar(10)//&
            "    smooth_shading=True,"//achar(10)//&
            "    split_sharp_edges=True,"//achar(10)//&
            ")"//achar(10)//&
            "p.add_axes(interactive=False)"//achar(10)//&
            "def point_picker_callback(point):"//achar(10)//&
            "    mesh = Xc"//achar(10)//&
            "    point_id = mesh.find_closest_point(point)"//achar(10)//&
            "    point_coords = mesh.points[point_id]"//achar(10)//&
            "    label = f'ID: {point_id + 1}\n({point_coords[0]:.3f}, {point_coords[1]:.3f}, {point_coords[2]:.3f})'"//achar(10)//&
            "    p.add_point_labels([point_coords],[label],font_size=14,text_color='black',show_points=False,fill_shape=False,shape=None)"//achar(10)//&
            "picker = p.enable_point_picking(callback=point_picker_callback, show_message=False)"//achar(10)//&
            "window_size = p.window_size"//achar(10)//&
            "y_pos = window_size[1]"//achar(10)//&
            "def Xcp_toggle_vis(flag): actor_Xcp.SetVisibility(flag)"//achar(10)//&
            "def Xcw_toggle_vis(flag): actor_Xcw.SetVisibility(flag)"//achar(10)//&
            "def Xg_toggle_vis(flag):  actor_Xg.SetVisibility(flag)"//achar(10)

        if (present(vtkfile_Xth_in_Xg)) then
            pyvista_script = trim(adjustl(pyvista_script))//achar(10)//&
                "Xth = pv.read('"//trim(vtkfile_Xth_in_Xg)//"')"//achar(10)//&
                "actor_Xth = p.add_mesh("//achar(10)//&
                "    Xth,"//achar(10)//&
                "    style='wireframe',"//achar(10)//&
                "    color='magenta',"//achar(10)//&
                "    line_width=2,"//achar(10)//&
                "    opacity=0.7"//achar(10)//&
                ")"//achar(10)//&
                "def Xth_toggle_vis(flag): actor_Xth.SetVisibility(flag)"
        end if

        pyvista_script = trim(adjustl(pyvista_script))//achar(10)//&
            "p.add_checkbox_button_widget(Xcp_toggle_vis, value=True, color_on='red',    size=25, position=(0, y_pos - 1*25))"//achar(10)//&
            "p.add_checkbox_button_widget(Xcw_toggle_vis, value=True, color_on='yellow', size=25, position=(0, y_pos - 2*25))"//achar(10)//&
            "p.add_checkbox_button_widget(Xg_toggle_vis,  value=True, color_on='cyan',   size=25, position=(0, y_pos - 3*25))"//achar(10)//&
            "p.add_text('Xc (Points)',             position=(28, y_pos - 1*25), font_size=8, color='black', font='times')"//achar(10)//&
            "p.add_text('Xc (Control geometry)',   position=(28, y_pos - 2*25), font_size=8, color='black', font='times')"//achar(10)//&
            "p.add_text('Xg (Geometry)',           position=(28, y_pos - 3*25), font_size=8, color='black', font='times')"

        if (present(vtkfile_Xth_in_Xg)) then
            pyvista_script = trim(adjustl(pyvista_script))//achar(10)//&
                "p.add_checkbox_button_widget(Xth_toggle_vis, value=True, color_on='magenta', size=25, position=(0, y_pos - 4*25))"//achar(10)//&
                "p.add_text('Xth (Parameter)', position=(28, y_pos - 4*25), font_size=8, color='black', font='times')"
        end if

        pyvista_script = trim(adjustl(pyvista_script))//achar(10)//&
            "p.add_text('ForCAD', position=(0.0, 10.0), font_size=14, color='black', font='times')"//achar(10)//&
            "p.add_text('https://github.com/gha3mi/forcad', position=(0.0, 0.0), font_size=7, color='blue', font='times')"//achar(10)//&
            "p.show(title='ForCAD', interactive=True)"//achar(10)//&
            "p.deep_clean()"//achar(10)//&
            "del p"

        call execute_command_line('python -c "'//trim(adjustl(pyvista_script))//'"')
        end block
#endif
    end subroutine
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure subroutine set_ring(this, center, radius1, radius2, length)
        class(nurbs_volume), intent(inout) :: this
        real(rk), intent(in), contiguous :: center(:)
        real(rk), intent(in) :: radius1, radius2, length
        real(rk), allocatable :: Xc(:,:), Wc(:), knot1(:), knot2(:), knot3(:)
        integer :: i

        if (.not. this%err%ok) return

        ! Define control points for ring
        allocate(Xc(28, 3))
        Xc(1,:) = [ 1.0_rk,  0.0_rk,              0.0_rk]
        Xc(2,:) = [ 1.0_rk,  sqrt(3.0_rk),        0.0_rk]
        Xc(3,:) = [-0.5_rk,  sqrt(3.0_rk)/2.0_rk, 0.0_rk]
        Xc(4,:) = [-2.0_rk,  0.0_rk,              0.0_rk]
        Xc(5,:) = [-0.5_rk, -sqrt(3.0_rk)/2.0_rk, 0.0_rk]
        Xc(6,:) = [ 1.0_rk, -sqrt(3.0_rk),        0.0_rk]
        Xc(7,:) = [ 1.0_rk,  0.0_rk,              0.0_rk]

        Xc(1:7,1:2) = Xc(1:7,1:2) * radius1

        Xc(8,:) = [ 1.0_rk,  0.0_rk,              0.0_rk]
        Xc(9,:) = [ 1.0_rk,  sqrt(3.0_rk),        0.0_rk]
        Xc(10,:)= [-0.5_rk,  sqrt(3.0_rk)/2.0_rk, 0.0_rk]
        Xc(11,:)= [-2.0_rk,  0.0_rk,              0.0_rk]
        Xc(12,:)= [-0.5_rk, -sqrt(3.0_rk)/2.0_rk, 0.0_rk]
        Xc(13,:)= [ 1.0_rk, -sqrt(3.0_rk),        0.0_rk]
        Xc(14,:)= [ 1.0_rk,  0.0_rk,              0.0_rk]

        Xc(8:14,1:2) = Xc(8:14,1:2) * radius2

        Xc(15,:)= [ 1.0_rk,  0.0_rk,              length]
        Xc(16,:)= [ 1.0_rk,  sqrt(3.0_rk),        length]
        Xc(17,:)= [-0.5_rk,  sqrt(3.0_rk)/2.0_rk, length]
        Xc(18,:)= [-2.0_rk,  0.0_rk,              length]
        Xc(19,:)= [-0.5_rk, -sqrt(3.0_rk)/2.0_rk, length]
        Xc(20,:)= [ 1.0_rk, -sqrt(3.0_rk),        length]
        Xc(21,:)= [ 1.0_rk,  0.0_rk,              length]

        Xc(15:21,1:2) = Xc(15:21,1:2) * radius1

        Xc(22,:)= [ 1.0_rk,  0.0_rk,              length]
        Xc(23,:)= [ 1.0_rk,  sqrt(3.0_rk),        length]
        Xc(24,:)= [-0.5_rk,  sqrt(3.0_rk)/2.0_rk, length]
        Xc(25,:)= [-2.0_rk,  0.0_rk,              length]
        Xc(26,:)= [-0.5_rk, -sqrt(3.0_rk)/2.0_rk, length]
        Xc(27,:)= [ 1.0_rk, -sqrt(3.0_rk),        length]
        Xc(28,:)= [ 1.0_rk,  0.0_rk,              length]

        Xc(22:28,1:2) = Xc(22:28,1:2) * radius2

        ! Translate the control points
        do i = 1, size(Xc, 1)
            Xc(i,:) = center + Xc(i,:)
        end do

        ! Define weights for the control points
        Wc = [1.0_rk, 0.5_rk, 1.0_rk, 0.5_rk, 1.0_rk, 0.5_rk, 1.0_rk,&
            1.0_rk, 0.5_rk, 1.0_rk, 0.5_rk, 1.0_rk, 0.5_rk, 1.0_rk,&
            1.0_rk, 0.5_rk, 1.0_rk, 0.5_rk, 1.0_rk, 0.5_rk, 1.0_rk,&
            1.0_rk, 0.5_rk, 1.0_rk, 0.5_rk, 1.0_rk, 0.5_rk, 1.0_rk]

        ! Define knot vector
        knot1 = [0.0_rk, 0.0_rk, 0.0_rk, 1.0_rk/3.0_rk, 1.0_rk/3.0_rk, 2.0_rk/3.0_rk, 2.0_rk/3.0_rk, 1.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]

        ! Set knot vector, control points, and weights
        call this%set(knot1, knot2, knot3, Xc, Wc)
    end subroutine
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure subroutine set_C(this, center, radius1, radius2, length)
        class(nurbs_volume), intent(inout) :: this
        real(rk), intent(in), contiguous :: center(:)
        real(rk), intent(in) :: radius1, radius2, length
        real(rk), allocatable :: Xc(:,:), Wc(:), knot1(:), knot2(:), knot3(:)
        integer :: i

        if (.not. this%err%ok) return

        ! Define control points for C-shape
        allocate(Xc(20, 3))
        Xc(1,:)= [ 1.0_rk,  0.0_rk,              0.0_rk]
        Xc(2,:)= [ 1.0_rk,  sqrt(3.0_rk),        0.0_rk]
        Xc(3,:)= [-0.5_rk,  sqrt(3.0_rk)/2.0_rk, 0.0_rk]
        Xc(4,:)= [-2.0_rk,  0.0_rk,              0.0_rk]
        Xc(5,:)= [-0.5_rk, -sqrt(3.0_rk)/2.0_rk, 0.0_rk]

        Xc(1:5,1:2) = Xc(1:5,1:2) * radius1

        Xc(6,:)= [ 1.0_rk,  0.0_rk,              0.0_rk]
        Xc(7,:)= [ 1.0_rk,  sqrt(3.0_rk),        0.0_rk]
        Xc(8,:)= [-0.5_rk,  sqrt(3.0_rk)/2.0_rk, 0.0_rk]
        Xc(9,:)= [-2.0_rk,  0.0_rk,              0.0_rk]
        Xc(10,:)=[-0.5_rk, -sqrt(3.0_rk)/2.0_rk, 0.0_rk]

        Xc(6:10,1:2) = Xc(6:10,1:2) * radius2

        Xc(11,:)= [ 1.0_rk,  0.0_rk,              length]
        Xc(12,:)= [ 1.0_rk,  sqrt(3.0_rk),        length]
        Xc(13,:)= [-0.5_rk,  sqrt(3.0_rk)/2.0_rk, length]
        Xc(14,:)= [-2.0_rk,  0.0_rk,              length]
        Xc(15,:)= [-0.5_rk, -sqrt(3.0_rk)/2.0_rk, length]

        Xc(11:15,1:2) = Xc(11:15,1:2) * radius1

        Xc(16,:)= [ 1.0_rk,  0.0_rk,              length]
        Xc(17,:)= [ 1.0_rk,  sqrt(3.0_rk),        length]
        Xc(18,:)= [-0.5_rk,  sqrt(3.0_rk)/2.0_rk, length]
        Xc(19,:)= [-2.0_rk,  0.0_rk,              length]
        Xc(20,:)= [-0.5_rk, -sqrt(3.0_rk)/2.0_rk, length]

        Xc(16:20,1:2) = Xc(16:20,1:2) * radius2

        ! Translate the control points
        do i = 1, size(Xc, 1)
            Xc(i,:) = center + Xc(i,:)
        end do

        ! Define weights for the control points
        Wc = [1.0_rk, 0.5_rk, 1.0_rk, 0.5_rk, 1.0_rk,&
            1.0_rk, 0.5_rk, 1.0_rk, 0.5_rk, 1.0_rk,&
            1.0_rk, 0.5_rk, 1.0_rk, 0.5_rk, 1.0_rk,&
            1.0_rk, 0.5_rk, 1.0_rk, 0.5_rk, 1.0_rk]

        ! Define knot vector
        knot1 = [0.0_rk, 0.0_rk, 0.0_rk, 1.0_rk/2.0_rk, 1.0_rk/2.0_rk, 1.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]

        ! Set knot vector, control points, and weights
        call this%set(knot1, knot2, knot3, Xc, Wc)
    end subroutine
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure subroutine set_half_ring(this, center, radius1, radius2, length)
        class(nurbs_volume), intent(inout) :: this
        real(rk), intent(in), contiguous :: center(:)
        real(rk), intent(in) :: radius1, radius2, length
        real(rk), allocatable :: Xc(:,:), Wc(:), knot1(:), knot2(:), knot3(:)
        integer :: i

        if (.not. this%err%ok) return

        ! Define control points for half ring
        allocate(Xc(20, 3))
        Xc(1,:)  = [ 0.5_rk, 0.0_rk, 0.0_rk]
        Xc(2,:)  = [ 0.5_rk, 0.5_rk, 0.0_rk]
        Xc(3,:)  = [ 0.0_rk, 0.5_rk, 0.0_rk]
        Xc(4,:)  = [-0.5_rk, 0.5_rk, 0.0_rk]
        Xc(5,:)  = [-0.5_rk, 0.0_rk, 0.0_rk]

        Xc(1:5,1:2) = Xc(1:5,1:2) * radius1

        Xc(6,:)  = [ 0.5_rk, 0.0_rk, 0.0_rk]
        Xc(7,:)  = [ 0.5_rk, 0.5_rk, 0.0_rk]
        Xc(8,:)  = [ 0.0_rk, 0.5_rk, 0.0_rk]
        Xc(9,:)  = [-0.5_rk, 0.5_rk, 0.0_rk]
        Xc(10,:) = [-0.5_rk, 0.0_rk, 0.0_rk]

        Xc(6:10,1:2) = Xc(6:10,1:2) * radius2

        Xc(11,:) = [ 0.5_rk, 0.0_rk, length]
        Xc(12,:) = [ 0.5_rk, 0.5_rk, length]
        Xc(13,:) = [ 0.0_rk, 0.5_rk, length]
        Xc(14,:) = [-0.5_rk, 0.5_rk, length]
        Xc(15,:) = [-0.5_rk, 0.0_rk, length]

        Xc(11:15,1:2) = Xc(11:15,1:2) * radius1

        Xc(16,:) = [ 0.5_rk, 0.0_rk, length]
        Xc(17,:) = [ 0.5_rk, 0.5_rk, length]
        Xc(18,:) = [ 0.0_rk, 0.5_rk, length]
        Xc(19,:) = [-0.5_rk, 0.5_rk, length]
        Xc(20,:) = [-0.5_rk, 0.0_rk, length]

        Xc(16:20,1:2) = Xc(16:20,1:2) * radius2

        ! Translate the control points
        do i = 1, size(Xc, 1)
            Xc(i,:) = center + Xc(i,:)
        end do

        ! Define weights for the control points
        Wc = [1.0_rk, 1.0_rk/sqrt(2.0_rk), 1.0_rk, 1.0_rk/sqrt(2.0_rk), 1.0_rk,&
            1.0_rk, 1.0_rk/sqrt(2.0_rk), 1.0_rk, 1.0_rk/sqrt(2.0_rk), 1.0_rk,&
            1.0_rk, 1.0_rk/sqrt(2.0_rk), 1.0_rk, 1.0_rk/sqrt(2.0_rk), 1.0_rk,&
            1.0_rk, 1.0_rk/sqrt(2.0_rk), 1.0_rk, 1.0_rk/sqrt(2.0_rk), 1.0_rk]

        ! Define knot vector
        knot1 = [0.0_rk, 0.0_rk, 0.0_rk, 1.0_rk/2.0_rk, &
            1.0_rk/2.0_rk, 1.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]

        ! Set knot vector, control points, and weights
        call this%set(knot1, knot2, knot3, Xc, Wc)
    end subroutine
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure subroutine nearest_point(this, point_Xg, nearest_Xg, nearest_Xt, id)
        class(nurbs_volume), intent(in) :: this
        real(rk), intent(in), contiguous :: point_Xg(:)
        real(rk), intent(out), optional :: nearest_Xg(size(point_Xg))
        real(rk), intent(out), optional :: nearest_Xt(3)
        integer, intent(out), optional :: id
        integer :: id_, i
        real(rk), allocatable :: distances(:)

        if (.not. this%err%ok) return

        allocate(distances(this%ng(1)*this%ng(2)*this%ng(3)))

#if defined(__NVCOMPILER)
        do i = 1, this%ng(1)*this%ng(2)*this%ng(3)
#else
        do concurrent (i = 1: this%ng(1)*this%ng(2)*this%ng(3))
#endif
            distances(i) = norm2(this%Xg(i,:) - point_Xg)
        end do

        ! replaced minloc due to NVFortran bug
#if defined(__NVCOMPILER)
        id_ = 1
        do i = 2, size(distances)
            if (distances(i) < distances(id_)) id_ = i
        end do
#else
        id_ = minloc(distances, dim=1)
#endif
        if (present(id)) id = id_
        if (present(nearest_Xg)) nearest_Xg = this%Xg(id_,:)
        if (present(nearest_Xt)) nearest_Xt = this%Xt(id_,:)
    end subroutine
    !===============================================================================


    !===============================================================================
    !> 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), contiguous :: point_Xg(:)
        real(rk), intent(in) :: tol
        integer, intent(in) :: maxit
        real(rk), intent(out) :: nearest_Xt(3)
        real(rk), intent(out), optional :: nearest_Xg(size(this%Xc,2))

        real(rk) :: obj, obj_trial, grad(3), hess(3,3), dk(3)
        real(rk) :: alphak, alpha_max, alpha_i, tau, beta, eps
        real(rk) :: lower_bounds(3), upper_bounds(3), xt(3)
        real(rk), allocatable ::Tgc(:), dTgc(:,:), d2Tgc(:,:)
        real(rk) :: Xg(size(this%Xc,2)), xk(3), xkn(3)
        integer :: k, l, i
        logical :: convergenz
        type(nurbs_volume) :: copy_this

        if (.not. this%err%ok) return

        alphak = 0.0_rk
        dk     = 0.0_rk
        k      = 0
        eps    = 10.0_rk*tiny(1.0_rk)

        ! bounds
        lower_bounds = [minval(this%knot1), minval(this%knot2), minval(this%knot3)]
        upper_bounds = [maxval(this%knot1), maxval(this%knot2), maxval(this%knot3)]

        ! initial guess (coarse search)
        copy_this = this
        call copy_this%create(50,50,50)
        call copy_this%nearest_point(point_Xg=point_Xg, nearest_Xt=xk)
        call copy_this%finalize()

        ! clamp initial guess to bounds
        xk = max(min(xk, upper_bounds), lower_bounds)

        xkn = xk
        convergenz = .false.

        do while (.not. convergenz .and. k < maxit)

            ! objective, gradient, hessian
            Xg = this%cmp_Xg(xk)
            call this%derivative2(Xt=xk, d2Tgc=d2Tgc, dTgc=dTgc, Tgc=Tgc)  ! Tgc unused

            obj = norm2(Xg - point_Xg) + 0.001_rk  ! small epsilon to avoid divide-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 .or. (k>0 .and. norm2(xk-xkn) <= tol)) then
                convergenz  = .true.
                nearest_Xt  = xk
                if (present(nearest_Xg)) nearest_Xg = this%cmp_Xg(nearest_Xt)
            else
                ! Newton step
                dk = - matmul(inv(hess), grad)

                ! Backtracking-Armijo with feasibility (box constraints)
                tau  = 0.5_rk
                beta = 1.0e-4_rk

                ! compute maximum feasible step so xk + alpha*dk stays in [lower_bounds, upper_bounds]
                alpha_max = 1.0_rk
                do i = 1, 3
                    if (dk(i) > 0.0_rk) then
                        if (upper_bounds(i) > xk(i)) then
                            alpha_i = (upper_bounds(i) - xk(i)) / dk(i)
                            alpha_max = min(alpha_max, max(0.0_rk, alpha_i))
                        else
                            alpha_max = 0.0_rk
                        end if
                    else if (dk(i) < 0.0_rk) then
                        if (lower_bounds(i) < xk(i)) then
                            alpha_i = (lower_bounds(i) - xk(i)) / dk(i)
                            alpha_max = min(alpha_max, max(0.0_rk, alpha_i))
                        else
                            alpha_max = 0.0_rk
                        end if
                    end if
                end do

                if (alpha_max <= eps) then
                    convergenz = .true.
                    nearest_Xt = xk
                    if (present(nearest_Xg)) nearest_Xg = this%cmp_Xg(nearest_Xt)
                    exit
                end if

                alphak = min(1.0_rk, alpha_max)
                l = 0
                do
                    if (alphak <= eps .or. l >= 50) exit
                    xt = xk + alphak*dk        ! feasible since alphak ≤ alpha_max
                    obj_trial = norm2(this%cmp_Xg(xt) - point_Xg) + 0.001_rk
                    if (obj_trial <= obj + alphak*beta*dot_product(grad, dk)) exit
                    alphak = min(tau*alphak, alpha_max)  ! shrink but stay feasible
                    l = l + 1
                end do

                xkn = xk
                if (alphak > eps) then
                    xk = xk + alphak*dk
                end if

                ! clamp updated iterate
                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
    pure function cmp_elemFace(this, elem, face) result(elemConn)
        class(nurbs_volume), intent(in) :: this
        integer, intent(in) :: elem
        integer, intent(in) :: face
        integer, allocatable :: elemConn(:)
        integer :: n(3), ii, jj, k

        if (.not. this%err%ok) return

        !> number of nodes in each direction
        n = this%degree + 1

        select case(face)
          case(1)
            allocate(elemConn(n(1)*n(2)))
            elemConn = this%elemConn(elem, 1 : n(1)*n(2))
          case(2)
            allocate(elemConn(n(1)*n(2)))
            elemConn = this%elemConn(elem, n(1)*n(2)*n(3)-n(1)*n(2)+1 : n(1)*n(2)*n(3))
          case(3)
            allocate(elemConn(n(1)*n(3)))
            k = 0
            do jj = 1, n(3)
                do ii = 1, n(1)
                    k = k+1
                    elemConn(k) = this%elemConn(elem, n(1)*n(2)*jj + ii - n(1)*n(2))
                end do
            end do
          case(4)
            allocate(elemConn(n(1)*n(3)))
            k = 0
            do jj = 1, n(3)
                do ii = 1, n(1)
                    k = k+1
                    elemConn(k) = this%elemConn(elem, n(1)*n(2)*jj + ii - n(1))
                end do
            end do
          case(5)
            allocate(elemConn(n(2)*n(3)))
            do ii = 1, n(2)*n(3)
                elemConn(ii) = this%elemConn(elem, n(1)*ii - n(1) + 1)
            end do
          case(6)
            allocate(elemConn(n(2)*n(3)))
            do ii = 1, n(2)*n(3)
                elemConn(ii) = this%elemConn(elem, n(1)*ii)
            end do
        end select
    end function
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure function cmp_elemFace_Xc_vis(this, elem, face) result(elemConn)
        class(nurbs_volume), intent(in) :: this
        integer, intent(in) :: elem
        integer, intent(in) :: face
        integer, allocatable :: elemConn(:)
        integer :: n(3), ii, jj, k

        if (.not. this%err%ok) return

        !> number of nodes in each direction
        n = [2,2,2]

        select case(face)
          case(1)
            allocate(elemConn(n(1)*n(2)))
            elemConn = this%elemConn_Xc_vis(elem, 1 : n(1)*n(2))
          case(2)
            allocate(elemConn(n(1)*n(2)))
            elemConn = this%elemConn_Xc_vis(elem, n(1)*n(2)*n(3)-n(1)*n(2)+1 : n(1)*n(2)*n(3))
          case(3)
            allocate(elemConn(n(1)*n(3)))
            k = 0
            do jj = 1, n(3)
                do ii = 1, n(1)
                    k = k+1
                    elemConn(k) = this%elemConn_Xc_vis(elem, n(1)*n(2)*jj + ii - n(1)*n(2))
                end do
            end do
          case(4)
            allocate(elemConn(n(1)*n(3)))
            k = 0
            do jj = 1, n(3)
                do ii = 1, n(1)
                    k = k+1
                    elemConn(k) = this%elemConn_Xc_vis(elem, n(1)*n(2)*jj + ii - n(1))
                end do
            end do
          case(5)
            allocate(elemConn(n(2)*n(3)))
            do ii = 1, n(2)*n(3)
                elemConn(ii) = this%elemConn_Xc_vis(elem, n(1)*ii - n(1) + 1)
            end do
          case(6)
            allocate(elemConn(n(2)*n(3)))
            do ii = 1, n(2)*n(3)
                elemConn(ii) = this%elemConn_Xc_vis(elem, n(1)*ii)
            end do
        end select
    end function
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure function cmp_elemFace_Xg_vis(this, elem, face) result(elemConn)
        class(nurbs_volume), intent(in) :: this
        integer, intent(in) :: elem
        integer, intent(in) :: face
        integer, allocatable :: elemConn(:)
        integer :: n(3), ii, jj, k

        if (.not. this%err%ok) return

        !> number of nodes in each direction
        n = [2,2,2]

        select case(face)
          case(1)
            allocate(elemConn(n(1)*n(2)))
            elemConn = this%elemConn_Xg_vis(elem, 1 : n(1)*n(2))
          case(2)
            allocate(elemConn(n(1)*n(2)))
            elemConn = this%elemConn_Xg_vis(elem, n(1)*n(2)*n(3)-n(1)*n(2)+1 : n(1)*n(2)*n(3))
          case(3)
            allocate(elemConn(n(1)*n(3)))
            k = 0
            do jj = 1, n(3)
                do ii = 1, n(1)
                    k = k+1
                    elemConn(k) = this%elemConn_Xg_vis(elem, n(1)*n(2)*jj + ii - n(1)*n(2))
                end do
            end do
          case(4)
            allocate(elemConn(n(1)*n(3)))
            k = 0
            do jj = 1, n(3)
                do ii = 1, n(1)
                    k = k+1
                    elemConn(k) = this%elemConn_Xg_vis(elem, n(1)*n(2)*jj + ii - n(1))
                end do
            end do
          case(5)
            allocate(elemConn(n(2)*n(3)))
            do ii = 1, n(2)*n(3)
                elemConn(ii) = this%elemConn_Xg_vis(elem, n(1)*ii - n(1) + 1)
            end do
          case(6)
            allocate(elemConn(n(2)*n(3)))
            do ii = 1, n(2)*n(3)
                elemConn(ii) = this%elemConn_Xg_vis(elem, n(1)*ii)
            end do
        end select
    end function
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure function cmp_degreeFace(this, face) result(degree)
        class(nurbs_volume), intent(in) :: this
        integer, intent(in) :: face
        integer :: degree(3)

        if (.not. this%err%ok) return

        select case (face)
          case(1)
            degree = [this%degree(1), this%degree(2), 0]
          case(2)
            degree = [this%degree(1), this%degree(2), 0]
          case(3)
            degree = [this%degree(1), 0, this%degree(3)]
          case(4)
            degree = [this%degree(1), 0, this%degree(3)]
          case(5)
            degree = [0, this%degree(2), this%degree(3)]
          case(6)
            degree = [0, this%degree(2), this%degree(3)]
          case default
            error stop 'Invalid face number'
        end select
    end function
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure subroutine ansatz(this, ie, ig, Tgc, dTgc_dXg, dV, ngauss)
        class(nurbs_volume), intent(inout) :: this
        integer, intent(in) :: ie, ig
        real(rk), intent(out) :: dV
        real(rk), allocatable, intent(out) :: Tgc(:), dTgc_dXg(:,:)
        integer, intent(in), optional :: ngauss(3)
        real(rk), allocatable :: Xth(:,:), Xth_e(:,:), Xth_eT(:,:), Xc_eT(:,:), Xth1(:), Xth2(:), Xth3(:), Xksi(:,:), Wksi(:)
        integer, allocatable :: elem_th(:,:), elem_c(:,:), elem_ce(:)
        type(nurbs_volume) :: th, th_e
        real(rk), allocatable :: dTtth_dXksi(:,:), Ttth(:), dTgc_dXt(:,:), Xt(:), dXt_dXksi(:,:), dXg_dXt(:,:)
        real(rk), allocatable :: dXg_dXksi(:,:) !! Jacobian matrix
        real(rk) :: det_dXg_dXksi !! Determinant of the Jacobian matrix
        real(rk) :: Xksii(3)

        if (.not. this%err%ok) return

        if (present(ngauss)) then
            call gauss_leg([0.0_rk, 1.0_rk], [0.0_rk, 1.0_rk], [0.0_rk, 1.0_rk], ngauss-1, Xksi, Wksi)
        else
            call gauss_leg([0.0_rk, 1.0_rk], [0.0_rk, 1.0_rk], [0.0_rk, 1.0_rk], this%degree, Xksi, Wksi)
        end if

        Xth1 = unique(this%knot1)
        Xth2 = unique(this%knot2)
        Xth3 = unique(this%knot3)
        call ndgrid(Xth1, Xth2, Xth3, Xth)
        call th%set([0.0_rk,Xth1,1.0_rk], [0.0_rk,Xth2,1.0_rk], [0.0_rk,Xth3,1.0_rk], Xth)
        elem_th = th%cmp_elem()

        elem_c = this%cmp_elem()

        Xth_e = Xth(elem_th(ie,:),:)
        call th_e%set([0.0_rk,0.0_rk,1.0_rk,1.0_rk], [0.0_rk,0.0_rk,1.0_rk,1.0_rk], [0.0_rk,0.0_rk,1.0_rk,1.0_rk], Xth_e)

        Xth_eT = transpose(Xth_e)
        elem_ce = elem_c(ie,:)
        Xc_eT = transpose(this%Xc(elem_ce,:))

        Xksii = Xksi(ig,:)
        call th_e%derivative(Xksii, dTtth_dXksi, Ttth)
        Xt = matmul(Xth_eT, Ttth)
        dXt_dXksi = matmul(Xth_eT, dTtth_dXksi)

        call this%derivative(Xt, dTgc_dXt, Tgc, elem_ce)
        dXg_dXt = matmul(Xc_eT, dTgc_dXt)

        dTgc_dXg = matmul(dTgc_dXt, inv(dXg_dXt))

        dXg_dXksi = matmul(dXg_dXt, dXt_dXksi)
        det_dXg_dXksi = det(dXg_dXksi)

        dV = det_dXg_dXksi*Wksi(ig)
    end subroutine
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure subroutine cmp_volume(this, volume, ngauss)
        class(nurbs_volume), intent(inout) :: this
        real(rk), intent(out) :: volume
        integer, intent(in), optional :: ngauss(3)
        real(rk), allocatable :: Tgc(:), dTgc_dXg(:,:)
        integer :: ie, ig
        integer :: ngauss_(3)
        real(rk) :: dV, dV_ig

        if (.not. this%err%ok) return

        if (present(ngauss)) then
            ngauss_ = ngauss
        else
            ngauss_ = this%degree + 1
        end if

        volume = 0.0_rk
#if defined(__NVCOMPILER)
        do ie = 1, size(this%cmp_elem(),1)
#else
        do concurrent (ie = 1:size(this%cmp_elem(),1)) reduce(+:volume)
#endif
            dV = 0.0_rk
            do ig = 1, product(ngauss_)
                call this%ansatz(ie, ig, Tgc, dTgc_dXg, dV_ig, ngauss_)
                dV = dV + dV_ig
            end do
            volume = volume + dV
        end do
    end subroutine
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure function cmp_Tgc_3d(Xti, knot1, knot2, knot3, nc, degree, Wc) result(Tgc)
        real(rk), intent(in), contiguous :: Xti(:)
        real(rk), intent(in), contiguous :: knot1(:)
        real(rk), intent(in), contiguous :: knot2(:)
        real(rk), intent(in), contiguous :: knot3(:)
        integer, intent(in) :: degree(3), nc(3)
        real(rk), intent(in), contiguous :: Wc(:)
        real(rk) :: Tgc(nc(1)*nc(2)*nc(3))
        real(rk) :: tmp
        integer :: i

        Tgc = kron(basis_bspline(Xti(3), knot3, nc(3), degree(3)),&
            kron(&
            basis_bspline(Xti(2), knot2, nc(2), degree(2)),&
            basis_bspline(Xti(1), knot1, nc(1), degree(1))))
        tmp = dot_product(Tgc, Wc)
        do concurrent (i = 1: nc(1)*nc(2)*nc(3))
           Tgc(i) = (Tgc(i) * Wc(i)) / tmp
        end do
    end function
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure function compute_Xg_nurbs_3d(Xt, knot1, knot2, knot3, degree, nc, ng, Xc, Wc) result(Xg)
        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), optional :: ng(3)
        real(rk), intent(in), contiguous :: Xc(:,:)
        real(rk), intent(in), contiguous :: Wc(:)
        real(rk), allocatable :: Xg(:,:)
        real(rk) :: Xti(size(Xt,2))
        integer :: i, ng_

        if (present(ng)) then
            ng_ = product(ng)
        else
            ng_ = size(Xt, 1)
        end if

        allocate(Xg(ng_, size(Xc,2)))
#if defined(__NVCOMPILER) || defined(__GFORTRAN__)
        do i = 1, ng_
#else
        do concurrent (i = 1: ng_) local(Xti)
#endif
            Xti = Xt(i,:)
            Xg(i,:) = matmul(cmp_Tgc_3d(Xti, knot1, knot2, knot3, nc, degree, Wc), Xc)
        end do
    end function
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure function compute_Xg_nurbs_3d_1point(Xt, knot1, knot2, knot3, degree, nc, Xc, Wc) result(Xg)
        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) :: Xg(size(Xc,2))
        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)))
        Xg = matmul(Tgc, Xc)
    end function
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure function compute_Xg_bspline_3d(Xt, knot1, knot2, knot3, degree, nc, ng, Xc) result(Xg)
        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), optional :: ng(3)
        real(rk), intent(in), contiguous :: Xc(:,:)
        real(rk), allocatable :: Xg(:,:)
        integer :: i, ng_

        if (present(ng)) then
            ng_ = product(ng)
        else
            ng_ = size(Xt, 1)
        end if

        allocate(Xg(ng_, size(Xc,2)))
#if defined(__NVCOMPILER)
        do i = 1, ng_
#else
        do concurrent (i = 1: ng_)
#endif
            Xg(i,:) = matmul(kron(basis_bspline(Xt(i,3), knot3, nc(3), degree(3)), kron(&
                basis_bspline(Xt(i,2), knot2, nc(2), degree(2)),&
                basis_bspline(Xt(i,1), knot1, nc(1), degree(1)))),&
                Xc)
        end do
    end function
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure function compute_Xg_bspline_3d_1point(Xt, knot1, knot2, knot3, degree, nc, Xc) result(Xg)
        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) :: 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
    pure subroutine compute_dTgc_nurbs_3d_vector(Xt, knot1, knot2, knot3, degree, nc, ng, Wc, dTgc, Tgc)
        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), optional :: ng(3)
        real(rk), intent(in), contiguous :: Wc(:)
        real(rk), allocatable, intent(out) :: dTgc(:,:,:)
        real(rk), allocatable, intent(out) :: Tgc(:,:)
        real(rk) :: dB1(nc(1)), dB2(nc(2)), dB3(nc(3))
        real(rk) :: B1(nc(1)), B2(nc(2)), B3(nc(3))
        real(rk) :: dBi(nc(1)*nc(2)*nc(3), 3), Bi(nc(1)*nc(2)*nc(3))
        integer :: i, ng_

        if (present(ng)) then
            ng_ = product(ng)
        else
            ng_ = size(Xt, 1)
        end if

        allocate(dTgc(ng_, nc(1)*nc(2)*nc(3), 3))
        allocate(Tgc(ng_, nc(1)*nc(2)*nc(3)))
#if defined(__NVCOMPILER) || defined(__GFORTRAN__)
        do i = 1, ng_
#else
        do concurrent (i = 1: ng_) local(B1, B2, B3, dB1, dB2, dB3, Bi, dBi)
#endif
            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
    end subroutine
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    !> If `elem` is not present: `Wc` refers to the full weight vector.
    !> If `elem` is present:     `Wc` refers to the element-local weight vector (`Wce`).
    pure subroutine compute_dTgc_nurbs_3d_scalar(Xt, knot1, knot2, knot3, degree, nc, Wc, dTgc, Tgc, elem)
        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(:)
        integer, intent(in), contiguous, optional :: elem(:)
        real(rk), allocatable, intent(out) :: dTgc(:,:)
        real(rk), allocatable, intent(out) :: Tgc(:)
        real(rk) :: dB1(nc(1)), dB2(nc(2)), dB3(nc(3))
        real(rk) :: B1(nc(1)), B2(nc(2)), B3(nc(3))
        real(rk), allocatable :: dBi(:,:), Bi(:)

        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)

        if (.not. present(elem)) then
            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)))

            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)
        else
            allocate(dTgc(size(elem), 3))
            allocate(Tgc(size(elem)))
            allocate(dBi(size(elem), 3), Bi(size(elem)))

            associate(Biall => kron( B3, kron( B2,  B1)))
                Bi = Biall(elem)
                Tgc = Bi*(Wc/(dot_product(Bi,Wc)))
            end associate

            associate(dB1all => kron(kron(B3,B2),dB1), dB2all => kron(kron(B3,dB2),B1), dB3all => kron(kron(dB3,B2),B1))
                dBi(:,1) = dB1all(elem)
                dBi(:,2) = dB2all(elem)
                dBi(:,3) = dB3all(elem)
            end associate

            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 if
    end subroutine
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure subroutine compute_dTgc_bspline_3d_vector(Xt, knot1, knot2, knot3, degree, nc, ng, dTgc, Tgc)
        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), optional :: ng(3)
        real(rk), allocatable, intent(out) :: dTgc(:,:,:)
        real(rk), allocatable, intent(out) :: Tgc(:,:)
        real(rk) :: dB1(nc(1)), dB2(nc(2)), dB3(nc(3))
        real(rk) :: B1(nc(1)), B2(nc(2)), B3(nc(3))
        integer :: i, ng_

        if (present(ng)) then
            ng_ = product(ng)
        else
            ng_ = size(Xt, 1)
        end if

        allocate(dTgc(ng_, nc(1)*nc(2)*nc(3), 3))
        allocate(Tgc(ng_, nc(1)*nc(2)*nc(3)))

#if defined(__NVCOMPILER) || defined(__GFORTRAN__)
        do i = 1, ng_
#else
        do concurrent (i = 1: ng_) local(B1, B2, B3, dB1, dB2, dB3)
#endif
            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

    end subroutine
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure subroutine compute_dTgc_bspline_3d_scalar(Xt, knot1, knot2, knot3, degree, nc, dTgc, Tgc, elem)
        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), contiguous, optional :: elem(:)
        real(rk), allocatable, intent(out) :: dTgc(:,:)
        real(rk), allocatable, intent(out) :: Tgc(:)
        real(rk) :: dB1(nc(1)), dB2(nc(2)), dB3(nc(3))
        real(rk) :: B1(nc(1)), B2(nc(2)), B3(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)

        if (.not. present(elem)) then
            allocate(dTgc(nc(1)*nc(2)*nc(3), 3))
            allocate(Tgc(nc(1)*nc(2)*nc(3)))

            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)
        else
            allocate(dTgc(size(elem), 3))
            allocate(Tgc(size(elem)))

            associate(B => kron(B3, kron(B2, B1)))
                Tgc = B(elem)
            end associate

            associate(dB1 => kron(kron(B3,B2),dB1), dB2 => kron(kron(B3,dB2),B1), dB3 => kron(kron(dB3,B2),B1))
                dTgc(:,1) = dB1(elem)
                dTgc(:,2) = dB2(elem)
                dTgc(:,3) = dB3(elem)
            end associate
        end if
    end subroutine
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure subroutine compute_d2Tgc_nurbs_3d_vector(Xt, knot1, knot2, knot3, degree, nc, ng, Wc, d2Tgc, dTgc, Tgc)
        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), optional :: 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) :: d2B1(nc(1)), d2B2(nc(2)), d2B3(nc(3))
        real(rk) :: dB1(nc(1)), dB2(nc(2)), dB3(nc(3))
        real(rk) :: B1(nc(1)), B2(nc(2)), B3(nc(3))
        real(rk), allocatable :: Tgci(:), dTgci(:)
        real(rk) :: d2Bi(3*nc(1)*nc(2)*nc(3), 3), dBi(nc(1)*nc(2)*nc(3), 3), Bi(nc(1)*nc(2)*nc(3))
        integer :: i, ng_

        if (present(ng)) then
            ng_ = product(ng)
        else
            ng_ = size(Xt, 1)
        end if

        allocate(Tgci(nc(1)*nc(2)*nc(3)), dTgci(nc(1)*nc(2)*nc(3)))
        allocate(d2Tgc(ng_, 3*nc(1)*nc(2)*nc(3), 3))
        allocate(dTgc(ng_, nc(1)*nc(2)*nc(3), 3))
        allocate(Tgc(ng_, nc(1)*nc(2)*nc(3)))
#if defined(__NVCOMPILER) || defined(__GFORTRAN__)
        do i = 1, ng_
#else
        do concurrent (i = 1: ng_) local(B1, B2, B3, dB1, dB2, dB3, d2B1, d2B2, d2B3, Bi, dBi, d2Bi)
#endif
            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
    pure subroutine compute_d2Tgc_nurbs_3d_scalar(Xt, knot1, knot2, knot3, degree, nc, Wc, d2Tgc, dTgc, Tgc)
        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) :: d2B1(nc(1)), d2B2(nc(2)), d2B3(nc(3))
        real(rk) :: dB1(nc(1)), dB2(nc(2)), dB3(nc(3))
        real(rk) :: B1(nc(1)), B2(nc(2)), B3(nc(3))
        real(rk), allocatable :: d2Bi(:,:), dBi(:,:), Bi(:)


        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
    pure subroutine compute_d2Tgc_bspline_3d_vector(Xt, knot1, knot2, knot3, degree, nc, ng, d2Tgc, dTgc, Tgc)
        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), optional :: ng(3)
        real(rk), allocatable, intent(out) :: d2Tgc(:,:,:)
        real(rk), allocatable, intent(out) :: dTgc(:,:,:)
        real(rk), allocatable, intent(out) :: Tgc(:,:)
        integer :: i, ng_
        real(rk) :: d2B1(nc(1)), d2B2(nc(2)), d2B3(nc(3))
        real(rk) :: dB1(nc(1)), dB2(nc(2)), dB3(nc(3))
        real(rk) :: B1(nc(1)), B2(nc(2)), B3(nc(3))

        if (present(ng)) then
            ng_ = product(ng)
        else
            ng_ = size(Xt, 1)
        end if

        allocate(d2Tgc(ng_, 3*nc(1)*nc(2)*nc(3), 3))
        allocate(dTgc(ng_, nc(1)*nc(2)*nc(3), 3))
        allocate(Tgc(ng_, nc(1)*nc(2)*nc(3)))
#if defined(__NVCOMPILER) || defined(__GFORTRAN__)
        do i = 1, ng_
#else
        do concurrent (i = 1: ng_) local(B1, B2, B3, dB1, dB2, dB3, d2B1, d2B2, d2B3)
#endif
            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
    pure subroutine compute_d2Tgc_bspline_3d_scalar(Xt, knot1, knot2, knot3, degree, nc, d2Tgc, dTgc, Tgc)
        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) :: d2B1(nc(1)), d2B2(nc(2)), d2B3(nc(3))
        real(rk) :: dB1(nc(1)), dB2(nc(2)), dB3(nc(3))
        real(rk) :: B1(nc(1)), B2(nc(2)), B3(nc(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)

        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
    pure function compute_Tgc_nurbs_3d_vector(Xt, knot1, knot2, knot3, degree, nc, ng, Wc) result(Tgc)
        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), optional :: ng(3)
        real(rk), intent(in), contiguous :: Wc(:)
        real(rk), allocatable :: Tgc(:,:)
        real(rk) :: Tgci(nc(1)*nc(2)*nc(3))
        integer :: i, ng_

        if (present(ng)) then
            ng_ = product(ng)
        else
            ng_ = size(Xt, 1)
        end if

        allocate(Tgc(ng_, nc(1)*nc(2)*nc(3)))
#if defined(__NVCOMPILER) || defined(__GFORTRAN__)
        do i = 1, ng_
#else
        do concurrent (i = 1: ng_) local(Tgci)
#endif
            Tgci = kron(basis_bspline(Xt(i,3), knot3, nc(3), degree(3)), kron(&
                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
    end function
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure function compute_Tgc_nurbs_3d_scalar(Xt, knot1, knot2, knot3, degree, nc, Wc) result(Tgc)
        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
    pure function compute_Tgc_bspline_3d_vector(Xt, knot1, knot2, knot3, degree, nc, ng) result(Tgc)
        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), optional :: ng(3)
        real(rk), allocatable :: Tgc(:,:)
        integer :: i, ng_

        if (present(ng)) then
            ng_ = product(ng)
        else
            ng_ = size(Xt, 1)
        end if

        allocate(Tgc(ng_, nc(1)*nc(2)*nc(3)))
#if defined(__NVCOMPILER)
        do i = 1, ng_
#else
        do concurrent (i = 1: ng_)
#endif
            Tgc(i,:) = kron(basis_bspline(Xt(i,3), knot3, nc(3), degree(3)), kron(&
                basis_bspline(Xt(i,2), knot2, nc(2), degree(2)),&
                basis_bspline(Xt(i,1), knot1, nc(1), degree(1))))
        end do
    end function
    !===============================================================================


    !===============================================================================
    !> author: Seyed Ali Ghasemi
    !> license: BSD 3-Clause
    pure function compute_Tgc_bspline_3d_scalar(Xt, knot1, knot2, knot3, degree, nc) result(Tgc)
        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
    pure subroutine lsq_fit_bspline(this, Xt, Xdata, ndata)
        use forcad_interface, only: solve
        class(nurbs_volume), intent(inout) :: this
        real(rk), intent(in), contiguous :: Xt(:,:), Xdata(:,:)
        integer, intent(in) :: ndata(3)
        real(rk), allocatable :: T(:,:), Tt(:,:), TtT(:,:), TtX(:,:)
        integer :: i, n

        if (.not. this%err%ok) return

        if (this%nc(1) > ndata(1)) then
            call this%err%set(&
                code       = 106,&
                severity   = 1,&
                category   = 'forcad_nurbs_volume',&
                message    = 'Invalid number of control points in the first direction.',&
                location   = 'lsq_fit_bspline',&
                suggestion = 'Ensure that the number of control points does not exceed the number of data points.')
            return
        end if
        if (this%nc(2) > ndata(2)) then
            call this%err%set(&
                code       = 106,&
                severity   = 1,&
                category   = 'forcad_nurbs_volume',&
                message    = 'Invalid number of control points in the second direction.',&
                location   = 'lsq_fit_bspline',&
                suggestion = 'Ensure that the number of control points does not exceed the number of data points.')
            return
        end if
        if (this%nc(3) > ndata(3)) then
            call this%err%set(&
                code       = 106,&
                severity   = 1,&
                category   = 'forcad_nurbs_volume',&
                message    = 'Invalid number of control points in the third direction.',&
                location   = 'lsq_fit_bspline',&
                suggestion = 'Ensure that the number of control points does not exceed the number of data points.')
            return
        end if

        n = ndata(1)*ndata(2)*ndata(3)

        allocate(T(n, this%nc(1)*this%nc(2)*this%nc(3)))
#if defined(__NVCOMPILER) || (defined(__GFORTRAN__) && (__GNUC__ < 15 || (__GNUC__ == 15 && __GNUC_MINOR__ < 1)))
        do i = 1, n
#else
        do concurrent (i = 1: n)
#endif
            T(i,:) = kron(&
            basis_bspline(Xt(i,3), this%knot3, this%nc(3), this%degree(3)), kron(&
            basis_bspline(Xt(i,2), this%knot2, this%nc(2), this%degree(2)),&
            basis_bspline(Xt(i,1), this%knot1, this%nc(1), this%degree(1))))
        end do
        Tt = transpose(T)
        TtT = matmul(Tt, T)
        TtX = matmul(Tt, Xdata)
        this%Xc = solve(TtT, TtX)
    end subroutine
    !===============================================================================

end module forcad_nurbs_volume