! Copyright (c) INRIA ! License: CECILL-B ! !> @ingroup splines !> @brief Module for 1D splines, linear combination of B-spline functions. !> !> @details Here we define a 1D spline type as an element of the linear space !> given by the span of the given B-splines (basis functions). !> Therefore, initialization of a 1D spline object requires an existing !> B-splines object, to which a private (polymorphic) pointer is !> associated. !> The B-spline coefficients are stored in a public allocatable array; !> at initialization the array is allocated to the proper size and all !> values are set to zero. !> In most situations the B-spline coefficients are not set directly by !> the end user, but are computed by some other object (e.g., a Poisson !> solver or a spline interpolator). !> Various public methods allow the user to evaluate the 1D spline S(x) !> and its derivative ∂S(x)/∂x any position x. !> !> @author Yaman Güçlü - IPP Garching !> @author Edoardo Zoni - IPP Garching module sll_m_spline_1d !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #include "sll_assert.h" use sll_m_working_precision, only: f64 use sll_m_bsplines_base , only: sll_c_bsplines implicit none public :: & sll_t_spline_1d private !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !> Working precision integer, parameter :: wp = f64 !> 1D spline type :: sll_t_spline_1d real(wp) , allocatable :: bcoef(:) class(sll_c_bsplines), pointer, private :: bspl => null() contains procedure :: init => s_spline_1d__init procedure :: free => s_spline_1d__free procedure :: belongs_to_space => f_spline_1d__belongs_to_space procedure :: eval => f_spline_1d__eval procedure :: eval_deriv => f_spline_1d__eval_deriv procedure :: eval_array => s_spline_1d__eval_array procedure :: eval_array_deriv => s_spline_1d__eval_array_deriv end type sll_t_spline_1d !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ contains !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !----------------------------------------------------------------------------- !> @brief Initialize 1D spline object as element of span(B-splines) !> @param[out] self 1D spline: new element of 1D spline space !> @param[in] bsplines B-splines: given basis of 1D spline space !----------------------------------------------------------------------------- subroutine s_spline_1d__init( self, bsplines ) class(sll_t_spline_1d), intent( out) :: self class(sll_c_bsplines ), intent(in ), target :: bsplines ! Store pointer to B-splines self%bspl => bsplines ! Allocate array of spline coefficients ! in case of periodic BCs, a larger array of coefficients is used in order ! to avoid a loop with calls to the "mod( , )" function at evaluation. associate( n => bsplines % nbasis, & g => merge( 1 + bsplines % degree/2, 0, bsplines % periodic ) ) allocate( self%bcoef(1-g:n+g) ) end associate ! Set all coefficients to zero self%bcoef = 0.0_f64 end subroutine s_spline_1d__init !----------------------------------------------------------------------------- !> @brief Destroy 1D spline (re-initialization is possible afterwards) !> @param[inout] self 1D spline !----------------------------------------------------------------------------- subroutine s_spline_1d__free( self ) class(sll_t_spline_1d), intent(inout) :: self deallocate( self % bcoef ) nullify ( self % bspl ) end subroutine s_spline_1d__free !----------------------------------------------------------------------------- !> @brief Check if 1D spline belongs to span of given B-splines !> @param[in] self 1D spline !> @param[in] bsplines B-splines !> @returns true/false answer !----------------------------------------------------------------------------- SLL_PURE function f_spline_1d__belongs_to_space( self, bsplines ) result( in_space ) class(sll_t_spline_1d), intent(in) :: self class(sll_c_bsplines ), intent(in), target :: bsplines logical :: in_space in_space = associated( self%bspl, bsplines ) end function f_spline_1d__belongs_to_space !----------------------------------------------------------------------------- !> @brief Evaluate value of 1D spline at location x: y=S(x) !> @param[in] self 1D spline !> @param[in] x evaluation point !> @returns spline value y=S(x) !----------------------------------------------------------------------------- SLL_PURE function f_spline_1d__eval( self, x ) result( y ) class(sll_t_spline_1d), intent(in) :: self real(wp) , intent(in) :: x real(wp) :: y real(wp) :: values(1:self%bspl%degree+1) integer :: jmin, jmax call self % bspl % eval_basis( x, values, jmin ) jmax = jmin + self%bspl%degree y = dot_product( self%bcoef(jmin:jmax), values ) end function f_spline_1d__eval !----------------------------------------------------------------------------- !> @brief Evaluate derivative of 1D spline at location x: y=S'(x) !> @param[in] self 1D spline !> @param[in] x evaluation point !> @returns spline derivative y=S'(x) !----------------------------------------------------------------------------- SLL_PURE function f_spline_1d__eval_deriv( self, x ) result( y ) class(sll_t_spline_1d), intent(in) :: self real(wp) , intent(in) :: x real(wp) :: y real(wp) :: derivs(1:self%bspl%degree+1) integer :: jmin, jmax call self % bspl % eval_deriv( x, derivs, jmin ) jmax = jmin + self%bspl%degree y = dot_product( self%bcoef(jmin:jmax), derivs ) end function f_spline_1d__eval_deriv !----------------------------------------------------------------------------- !> @brief Evaluate value of 1D spline at all locations in array x !> @param[in] self 1D spline !> @param[in] x array of evaluation points x[i] !> @param[out] y array of spline values y[i]=S(x[i]) !----------------------------------------------------------------------------- SLL_PURE subroutine s_spline_1d__eval_array( self, x, y ) class(sll_t_spline_1d), intent(in ) :: self real(wp) , intent(in ) :: x(:) real(wp) , intent( out) :: y(:) integer :: i SLL_ASSERT( size(x) == size(y) ) do i = 1, size(x) y(i) = f_spline_1d__eval( self, x(i) ) end do end subroutine s_spline_1d__eval_array !----------------------------------------------------------------------------- !> @brief Evaluate derivative of 1D spline at all locations in array x !> @param[in] self 1D spline !> @param[in] x array of evaluation points x[i] !> @param[out] y array of spline derivatives y[i]=S'(x[i]) !----------------------------------------------------------------------------- SLL_PURE subroutine s_spline_1d__eval_array_deriv( self, x, y ) class(sll_t_spline_1d), intent(in ) :: self real(wp) , intent(in ) :: x(:) real(wp) , intent( out) :: y(:) integer :: i SLL_ASSERT( size(x) == size(y) ) do i = 1, size(x) y(i) = f_spline_1d__eval_deriv( self, x(i) ) end do end subroutine s_spline_1d__eval_array_deriv end module sll_m_spline_1d