test_splines.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/.
!||====================================================================
!||    test_spline          ../starter/source/materials/tools/test_splines.F
!||--- calls      -----------------------------------------------------
!||    spline_interpol_2d   ../starter/source/materials/tools/spline_interpol_2d.F
!||====================================================================
       SUBROUTINE test_spline(NPT0,NSUB,XF,YF,XX,YY)
C-----------------------------------------------
C   D e s c r i p t i o n
C-----------------------------------------------
       !compute SPLINE length with third order method (SIMPSON)
       !  INPUT  - LOCAL_PT : 4 control points P0,P1,P2,P3
       !  INPUT  - ALPHA : CCR SPLINE PARAMETER
       !  INPUT  - T : position [0,1] on SPLINE [T1,T2].  1.0 means full spline length
       !  OUTPUT - LEN : length of the curve parametrised with t in [0,T]  T<=1.0
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER  ,INTENT(IN)    :: NPT0   ! number of points of input function
      INTEGER  ,INTENT(IN)    :: NSUB   ! 
      my_real  ,DIMENSION(NPT0)   ,INTENT(IN)  :: xf,yf   ! initial curve coordinates
      my_real  ,DIMENSION((NPT0-1)*NSUB+1) ,INTENT(OUT) :: xx,yy   ! curve coordinates build with splines
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER :: I,J,K,NPTS,NSEG
      my_real :: tt,dx,dy,nx,ny
      my_real ,DIMENSION(:,:)  ,ALLOCATABLE :: spline_knots
      my_real ,DIMENSION(:)    ,ALLOCATABLE :: ctrl_ptx,ctrl_pty
      my_real ,DIMENSION(4)  :: ptx,pty,knots
      my_real ,PARAMETER  :: alpha = 0.5
C-----------------------------------------------
C   S o u r c e   L i n e s
c=======================================================================
      npts = npt0 + 2
      nseg = npt0 - 1
 
      ALLOCATE (ctrl_ptx(npts))
      ALLOCATE (ctrl_pty(npts))
      ALLOCATE (spline_knots(nseg,4))
 
c     calculate spline control points
      ctrl_ptx(2:npts-1) = xf(1:npt0)
      ctrl_pty(2:npts-1) = yf(1:npt0)
      ! Add start point - minimum of bending energy            
      ctrl_ptx(1) = (half*ctrl_ptx(2) - four*ctrl_ptx(3) + ctrl_ptx(4)) * half                  
      ctrl_pty(1) = (half*ctrl_pty(2) - four*ctrl_pty(3) + ctrl_pty(4)) * half                  
      ! Add end point - minimum of bending energy               
      ctrl_ptx(npts) = (ctrl_ptx(npts-3) - four*ctrl_ptx(npts-2) + five*ctrl_ptx(npts-1)) * half                  
      ctrl_pty(npts) = (ctrl_pty(npts-3) - four*ctrl_pty(npts-2) + five*ctrl_pty(npts-1)) * half                  
c
      k = 0
      DO i = 1,nseg                         
        ptx(1) = ctrl_ptx(i)                  
        ptx(2) = ctrl_ptx(i+1)                  
        ptx(3) = ctrl_ptx(i+2)                  
        ptx(4) = ctrl_ptx(i+3)
        pty(1) = ctrl_pty(i)                  
        pty(2) = ctrl_pty(i+1)                  
        pty(3) = ctrl_pty(i+2)                  
        pty(4) = ctrl_pty(i+3)
c
        knots(1) = zero
        dx = ptx(2) - ptx(1)
        dy = pty(2) - pty(1)
        knots(2) = spline_knots(i,1) + exp(alpha*log(sqrt(dx**2 + dy**2)))
        dx = ptx(3) - ptx(2)
        dy = pty(3) - pty(2)
        knots(3) = knots(2) + exp(alpha*log(sqrt(dx**2 + dy**2)))
        dx = ptx(4) - ptx(3)
        dy = pty(4) - pty(3)
        knots(4) = knots(3) + exp(alpha*log(sqrt(dx**2 + dy**2)))
        spline_knots(i,1:4) = knots(1:4) 
c
        DO j = 1,nsub                                      
          tt = (j-one) / nsub
          k = k + 1
          CALL spline_interpol_2d(ptx, pty ,knots, tt,  nx  ,ny  )
          xx(k) = nx
          yy(k) = ny
        ENDDO                                            
        
      ENDDO 
      ! last point                                           
      tt = one
      k = k + 1
      xx(k) = xf(npt0)
      yy(k) = yf(npt0)
 
c      CALL SPLINE_INTERPOL_2D(PTX, PTY ,KNOTS, TT,  NX  ,NY  )
c      XX(K) = NX
c      YY(K) = NY
c
      DEALLOCATE (spline_knots)
      DEALLOCATE (ctrl_pty)
      DEALLOCATE (ctrl_ptx)
c-----------
      RETURN
      END SUBROUTINE