cr_spline_interpol.F

Go to the documentation of this file.

Copyright>        OpenRadioss
Copyright>        Copyright (C) 1986-2025 Altair Engineering Inc.
Copyright>
Copyright>        This program is free software: you can redistribute it and/or modify
Copyright>        it under the terms of the GNU Affero General Public License as published by
Copyright>        the Free Software Foundation, either version 3 of the License, or
Copyright>        (at your option) any later version.
Copyright>
Copyright>        This program is distributed in the hope that it will be useful,
Copyright>        but WITHOUT ANY WARRANTY; without even the implied warranty of
Copyright>        MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
Copyright>        GNU Affero General Public License for more details.
Copyright>
Copyright>        You should have received a copy of the GNU Affero General Public License
Copyright>        along with this program.  If not, see <https://www.gnu.org/licenses/>.
Copyright>
Copyright>
Copyright>        Commercial Alternative: Altair Radioss Software
Copyright>
Copyright>        As an alternative to this open-source version, Altair also offers Altair Radioss
Copyright>        software under a commercial license.  Contact Altair to discuss further if the
Copyright>        commercial version may interest you: https://www.altair.com/radioss/.
!||====================================================================
!||    cr_spline_interpol     ../common_source/tools/interpolation/cr_spline_interpol.F
!||--- called by ------------------------------------------------------
!||    cr_spline_length       ../common_source/tools/interpolation/cr_spline_length.F
!||    cr_spline_point_proj   ../common_source/tools/interpolation/cr_spline_point_proj.F
!||====================================================================
      SUBROUTINE cr_spline_interpol( PTS , KNOTS, T, C, C_D, C_DD)
C-----------------------------------------------
        !FOR A GIVEN POSITION T in [0,1] which is position on [P1,P2], GET VALUE ON SPLINE P0-P1-P2-P3
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C   D e s c r i p t i o n
C-----------------------------------------------
        ! INPUT  - PTS(0:3,:) are 4 control points
        ! INPUT  - T is in [0,1]  : T=0 => VAL = P(2,:)
        ! OUTPUT - C is interpolated value.
        ! OUTPUT - C_D is 1st derivative
        ! OUTPUT - C_DD is 2nd derivatice
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
        my_real, INTENT(IN)    :: pts(0:3,3),knots(0:3),t
        my_real, INTENT(INOUT) :: c(3),c_d(3),c_dd(3)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
        my_real :: a1(3),a2(3),a3(3)
        my_real :: b1(3),b2(3)
        my_real :: a1p(3),a2p(3),a3p(3), b1p(3),b2p(3),b1pp(3),b2pp(3)
        my_real :: tt  !in [T1,T2]=[KNOT(1),KNOT(2)]
C-----------------------------------------------
C   S o u r c e   L i n e s
C-----------------------------------------------
        tt=knots(1)+t*(knots(2)-knots(1))
        a1(1:3)   = ((knots(1)-tt)*pts(0,1:3)+(tt-knots(0))*pts(1,1:3)) /(knots(1)-knots(0))
        a2(1:3)   = ((knots(2)-tt)*pts(1,1:3)+(tt-knots(1))*pts(2,1:3)) /(knots(2)-knots(1))
        a3(1:3)   = ((knots(3)-tt)*pts(2,1:3)+(tt-knots(2))*pts(3,1:3)) /(knots(3)-knots(2))
        b1(1:3)   = ((knots(2)-tt)*a1(1:3) + (tt-knots(0))*a2(1:3)) / (knots(2)-knots(0))
        b2(1:3)   = ((knots(3)-tt)*a2(1:3) + (tt-knots(1))*a3(1:3)) / (knots(3)-knots(1))
        c(1:3)    = ((knots(2)-tt)*b1(1:3) + (tt-knots(1))*b2(1:3)) / (knots(2)-knots(1))
        a1p(1:3)  = (pts(1,1:3)-pts(0,1:3))/(knots(1)-knots(0))
        a2p(1:3)  = (pts(2,1:3)-pts(1,1:3))/(knots(2)-knots(1))
        a3p(1:3)  = (pts(3,1:3)-pts(2,1:3))/(knots(3)-knots(2))
        b1p(1:3)  =  (a2(1:3)-a1(1:3))/(knots(2)-knots(0))
     .            + (knots(2)-tt)/(knots(2)-knots(0))*a1p(1:3)
     .            + (tt-knots(0))/(knots(2)-knots(0))*a2p(1:3)
        b2p(1:3)  =  (a3(1:3)-a2(1:3))/(knots(3)-knots(1))
     .            + (knots(3)-tt)/(knots(3)-knots(1))*a2p(1:3)
     .            + (tt-knots(1))/(knots(3)-knots(1))*a3p(1:3)
        c_d(1:3)  =  (b2(1:3)-b1(1:3))/(knots(2)-knots(1))
     .            + (knots(2)-tt)/(knots(2)-knots(1))*b1p(1:3)
     .            + (tt-knots(1))/(knots(2)-knots(1))*b2p(1:3)
        b1pp(1:3) =  (2.*a2p(1:3)-2.*a1p(1:3))/(knots(2)-knots(0))
        b2pp(1:3) =  (2.*a3p(1:3)-2.*a2p(1:3))/(knots(3)-knots(1))
        c_dd(1:3) = (2.*b2p(1:3)-2.*b1p(1:3)+(knots(2)-tt)*b1pp(1:3)+(tt-knots(1))*b2pp(1:3))/(knots(2)-knots(1))
        RETURN
      END SUBROUTINE cr_spline_interpol