example_surface_1.f90 Source File

!> This program demonstrates the usage of a NURBS (Non-Uniform Rational B-Spline) surface object to create  and finalize a NURBS surface.
!> It sets up control points, weights, and knot vectors for all three dimensions, generates the surface, and exports the control points and the surface to VTK files.

program example_surface_1

    use forcad, only: rk, nurbs_surface

    implicit none

    type(nurbs_surface) :: nurbs            !! Declare a NURBS surface object
    real(rk), allocatable :: Xc(:,:), Wc(:) !! Arrays for control points and weights
    real(rk) :: knot1(6), knot2(6)          !! Arrays for knot vectors in both dimensions

    !-----------------------------------------------------------------------------
    ! Setting up the NURBS surface
    !-----------------------------------------------------------------------------

    !> Define control points for the NURBS surface
    Xc = generate_Xc(3, 3, 1.0_rk)

    !> Define weights for the control points
    allocate(Wc(size(Xc, 1)))
    Wc = 1.0_rk
    Wc(2) = 2.0_rk

    !> Define knot vectors for both dimensions
    knot1 = [0.0_rk, 0.0_rk, 0.0_rk, 1.0_rk, 1.0_rk, 1.0_rk]
    knot2 = [0.0_rk, 0.0_rk, 0.0_rk, 1.0_rk, 1.0_rk, 1.0_rk]

    !> Set knot vectors, control points, and weights for the NURBS surface object
    call nurbs%set(knot1, knot2, Xc, Wc)

    !> Deallocate local arrays
    deallocate(Xc, Wc)

    !> Export parameter space to a VTK file
    call nurbs%export_Xth('vtk/nurbs_surface_Xth.vtk')

    !> Export the control points to a VTK file
    call nurbs%export_Xc('vtk/nurbs_surface_Xc.vtk')

    !-----------------------------------------------------------------------------
    ! Creating the NURBS surface
    !-----------------------------------------------------------------------------

    !> Generate the NURBS surface with resolutions of 30 in both dimensions
    call nurbs%create(30, 30)

    !> Export the generated surface to a VTK file
    call nurbs%export_Xg('vtk/nurbs_surface_Xg.vtk')

    !> Export the NURBS surface to an IGES file
    call nurbs%export_iges('iges/nurbs_surface.iges')

    !-----------------------------------------------------------------------------
    ! Visualization using PyVista
    ! Note: PyVista is required for visualization. Install it using `pip install pyvista`
    !-----------------------------------------------------------------------------

    !> Show the control geometry and geometry using PyVista
    call nurbs%show('vtk/nurbs_surface_Xc.vtk','vtk/nurbs_surface_Xg.vtk')

    !-----------------------------------------------------------------------------
    ! Refinements
    !-----------------------------------------------------------------------------

    !> Print size of the knot vectors
    print*, size(nurbs%get_knot(1))
    print*, size(nurbs%get_knot(2))

    !> Insert knots 0.25, twice and 0.75, once in both directions
    call nurbs%insert_knots(1, [0.25_rk, 0.75_rk], [2,1]) ! direction 1
    call nurbs%insert_knots(2, [0.25_rk, 0.75_rk], [2,1]) ! direction 2

    !> Print size of the knot vectors after inserting knots
    print*, size(nurbs%get_knot(1))
    print*, size(nurbs%get_knot(2))

    !> Print the degrees
    print*, nurbs%get_degree()

    !> Elevate degree by 2 in both directions
    call nurbs%elevate_degree(1, 2) ! direction 1
    call nurbs%elevate_degree(2, 2) ! direction 2

    !> Print the degrees after elevating
    print*, nurbs%get_degree()

    !> Print size of the knot vectors
    print*, size(nurbs%get_knot(1))
    print*, size(nurbs%get_knot(2))

    !> Remove knots 0.25, twice and 0.75, once in both directions
    call nurbs%remove_knots(1, [0.25_rk, 0.75_rk], [2,1]) ! direction 1
    call nurbs%remove_knots(2, [0.25_rk, 0.75_rk], [2,1]) ! direction 2

    !> Print size of the knot vectors after removing knots
    print*, size(nurbs%get_knot(1))
    print*, size(nurbs%get_knot(2))

    !> Generate the refined NURBS surface with resolutions of 30 in both dimensions
    call nurbs%create()

    !> Export refined parameter space to a VTK file
    call nurbs%export_Xth('vtk/nurbs_surface_Xth2.vtk')

    !> Export updated control points to a VTK file
    call nurbs%export_Xc('vtk/nurbs_surface_Xc2.vtk')

    !> Export the refined generated surface to a VTK file
    call nurbs%export_Xg('vtk/nurbs_surface_Xg2.vtk')

    !> Export the NURBS surface to an IGES file
    call nurbs%export_iges('iges/nurbs_surface2.iges')

    !-----------------------------------------------------------------------------
    ! Visualization using PyVista
    ! Note: PyVista is required for visualization. Install it using `pip install pyvista`
    !-----------------------------------------------------------------------------

    !> Show the control geometry and geometry using PyVista
    call nurbs%show('vtk/nurbs_surface_Xc2.vtk','vtk/nurbs_surface_Xg2.vtk')

    !-----------------------------------------------------------------------------
    ! Transformations
    !-----------------------------------------------------------------------------

    !> Rotate the control points
    call nurbs%rotate_Xc(alpha=45.0_rk, beta=0.0_rk, theta=90.0_rk)

    !> Rotate the generated curve
    call nurbs%rotate_Xg(alpha=-45.0_rk, beta=0.0_rk, theta=-90.0_rk)

    !> Translate the control points
    call nurbs%translate_Xc([1.0_rk, 2.0_rk, -3.0_rk])

    !> Translate the generated curve
    call nurbs%translate_Xg([-1.0_rk, -2.0_rk, 3.0_rk])

    !> Export parameter space to a VTK file
    call nurbs%export_Xth('vtk/nurbs_surface_Xth3.vtk')

    !> Export the transformed control points to a VTK file
    call nurbs%export_Xc('vtk/nurbs_surface_Xc3.vtk')

    !> Export the transformed generated volume to a VTK file
    call nurbs%export_Xg('vtk/nurbs_surface_Xg3.vtk')

    !> Export the transformed NURBS surface to an IGES file
    call nurbs%export_iges('iges/nurbs_surface3.iges')

    !-----------------------------------------------------------------------------
    ! Visualization using PyVista
    ! Note: PyVista is required for visualization. Install it using `pip install pyvista`
    !-----------------------------------------------------------------------------

    !> Show the control geometry and geometry using PyVista
    call nurbs%show('vtk/nurbs_surface_Xc3.vtk','vtk/nurbs_surface_Xg3.vtk')

    !-----------------------------------------------------------------------------
    ! Finalizing
    !-----------------------------------------------------------------------------

    !> Finalize the NURBS surface object
    call nurbs%finalize()

contains

    !-----------------------------------------------------------------------------
    pure function generate_Xc(num_rows, num_cols, peak_height) result(control_points)
        integer,  intent(in) :: num_rows, num_cols
        real(rk), intent(in) :: peak_height
        real(rk), allocatable :: control_points(:,:)

        integer  :: i, j, idx
        real(rk) :: x_spacing, y_spacing, x_offset, y_offset
        real(rk) :: x, y, z

        x_spacing = merge(0.0_rk, 1.0_rk/real(num_cols-1, rk), num_cols == 1)
        y_spacing = merge(0.0_rk, 1.0_rk/real(num_rows-1, rk), num_rows == 1)

        x_offset = -0.5_rk
        y_offset = -0.5_rk

        allocate(control_points(num_rows*num_cols, 3))

        do concurrent (i = 1:num_rows, j = 1:num_cols) local(idx, x, y, z)
            idx = (i-1)*num_cols+j

            x = x_offset+real(j-1, rk)*x_spacing
            y = y_offset+real(i-1, rk)*y_spacing
            z = peak_height*exp(-(x**2+y**2))+0.1_rk*peak_height

            control_points(idx, 1) = x
            control_points(idx, 2) = y
            control_points(idx, 3) = z
        end do
    end function
    !-----------------------------------------------------------------------------

end program example_surface_1