spline_interpol_2d.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/.
!||====================================================================
!||    spline_interpol_2d   ../starter/source/materials/tools/spline_interpol_2d.F
!||--- called by ------------------------------------------------------
!||    test_spline          ../starter/source/materials/tools/test_splines.F
!||====================================================================
      SUBROUTINE spline_interpol_2d(PTX, PTY ,KNOTS, T,  NX  ,NY  )
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)    :: t
      my_real, INTENT(IN)    :: ptx(4)
      my_real, INTENT(IN)    :: pty(4)
      my_real, INTENT(IN)    :: knots(4)
      my_real, INTENT(OUT)   :: nx,ny
c      my_real, INTENT(INOUT) :: C(4),C_D(4),C_DD(4)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      my_real :: dk21,dk32,dk43,dk31,dk42,kt1,kt2,kt3,kt4
      my_real :: tt  ! in [T1,T2]=[KNOT(2),KNOT(3)]
      my_real :: a1x,a1y,a2x,a2y,a3x,a3y,b1x,b1y,b2x,b2y
      my_real :: cx,cy,c_dx,c_dy,c_ddx,c_ddy
      my_real :: a1px,a1py,a2px,a2py,a3px,a3py,b1px,b1py,b2px,b2py
      my_real :: b1ppx,b1ppy,b2ppx,b2ppy
C-----------------------------------------------
C   S o u r c e   L i n e s
c=======================================================================
      dk21 = knots(2) - knots(1)
      dk32 = knots(3) - knots(2)
      dk43 = knots(4) - knots(3)
      dk31 = knots(3) - knots(1)
      dk42 = knots(4) - knots(2)
      tt   = knots(2) + t * dk32
      kt1  = knots(1) - tt
      kt2  = knots(2) - tt
      kt3  = knots(3) - tt
      kt4  = knots(4) - tt
c
      a1x  = (kt2*ptx(1) - kt1*ptx(2)) / dk21
      a1y  = (kt2*pty(1) - kt1*pty(2)) / dk21
      a2x  = (kt3*ptx(2) - kt2*ptx(3)) / dk32
      a2y  = (kt3*pty(2) - kt2*pty(3)) / dk32
      a3x  = (kt4*ptx(3) - kt3*ptx(4)) / dk43
      a3y  = (kt4*pty(3) - kt3*pty(4)) / dk43
c      
      b1x  = (kt3*a1x - kt1*a2x) / dk31
      b1y  = (kt3*a1y - kt1*a2y) / dk31
      b2x  = (kt4*a2x - kt2*a3x) / dk42
      b2y  = (kt4*a2y - kt2*a3y) / dk42
c
      nx   = (kt3*b1x - kt2*b2x) / dk32 
      ny   = (kt3*b1y - kt2*b2y) / dk32 
c
      a1px = (ptx(2) - ptx(1)) / dk21
      a1py = (pty(2) - pty(1)) / dk21
      a2px = (ptx(3) - ptx(2)) / dk32
      a2py = (pty(3) - pty(2)) / dk32
      a3px = (ptx(4) - ptx(3)) / dk43
      a3py = (pty(4) - pty(3)) / dk43
      b1px = (a2x-a1x)/ dk31 + kt3 / dk31*a1px - kt1 / dk31*a2px  
      b1py = (a2y-a1y)/ dk31 + kt3 / dk31*a1py - kt1 / dk31*a2py  
      b2px = (a3x-a2x)/ dk42 + kt4 / dk42*a2px - kt2 / dk42*a3px 
      b2py = (a3y-a2y)/ dk42 + kt4 / dk42*a2py - kt2 / dk42*a3py 
!      C_DX = (B2-B1Y)/ DK32 + KT3 / DK32*B1P - KT2 / DK32*B2P
!      C_DY = (B2-B1)/ DK32 + KT3 / DK32*B1P - KT2 / DK32*B2P
      b1ppx= two*(a2px - a1px) / dk31
      b1ppy= two*(a2py - a1py) / dk31
      b2ppx= two*(a3px - a2px) / dk42 
      b2ppy= two*(a3py - two*a2py) / dk42 
!      C_DDX= (TWO*B2P-TWO*B1P + KT3*B1PP - KT2*B2PP) / DK32
!      C_DDY= (TWO*B2P-TWO*B1P + KT3*B1PP - KT2*B2PP) / DK32
c------------
      RETURN
      END SUBROUTINE spline_interpol_2d