vw_smooth.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/.
!||====================================================================
!||    vw_smooth     ../starter/source/tools/curve/vw_smooth.F
!||--- called by ------------------------------------------------------
!||    law70_table   ../starter/source/materials/mat/mat070/law70_table.F
!||    law70_upd     ../starter/source/materials/mat/mat070/law70_upd.F
!||====================================================================
      SUBROUTINE vw_smooth(NPT,NTARGET,X,Y)
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
      INTEGER                     ,INTENT(INOUT)  :: NPT      ! number of input function points
      INTEGER                     ,INTENT(IN)     :: NTARGET  ! target length of abscissa vector
      my_real ,DIMENSION(NPT)     ,INTENT(INOUT)  :: x,y
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER :: I,I0,I1,I2,IDX,IMIN,IPREV,INEXT,ITER,NITER,IZERO
      my_real :: amin,aa,s,t
      INTEGER ,DIMENSION(NPT) :: PREV,NEXT
      my_real ,DIMENSION(NPT) :: area
c-----------------------------------------------
c     smooth input curve using Visvalingam-Whyatt algorithm 
c     reduce number of points from NPT to NTARGET
c     first and last points are never eliminated, as well as (0,0)
c=======================================================================
      niter = npt - ntarget
      IF (niter < 1) RETURN
      
      prev(1)   = 0
      next(1)   = 2
      area(1)   = infinity / ten
      prev(npt) = npt - 1
      next(npt) = 10000000
      area(npt) = infinity / ten 
c
      DO i = 2,npt-1
        i1 = i-1
        i2 = i+1
        prev(i) = i1
        next(i) = i2
        IF (x(i) == zero .and. y(i) == zero) THEN
          area(i) = infinity / ten
        ELSE
          area(i) = x(i1)*y(i)  + x(i)*y(i2) + x(i2)*y(i1)
     .            - x(i1)*y(i2) - x(i)*y(i1) - x(i2)*y(i)
          area(i)  = abs(area(i))
        END IF
      END DO
c----------------------------------------------------      
c     eliminate points with min area until NPT = target number of points
c----------------------------------------------------      
      imin = 1
      DO iter = 1,niter                           
        amin = infinity
        idx  = 1
        DO i=1,npt
          IF (area(i) < amin) THEN
            amin = area(i)
            imin = i
          ENDIF
        END DO      
c     
        iprev = prev(imin)
        inext = next(imin)      
        next(iprev) = inext
        prev(inext) = iprev     
        area(imin)  = infinity
c
        ! recalculate PREV POINT AREA
        i0 = iprev
        i1 = prev(iprev)
        i2 = inext
        IF (i1 > 0) THEN
          aa = x(i0)*y(i0) + x(i0)*y(i2) + x(i2)*y(i1) 
     .       - x(i0)*y(i2) - x(i0)*y(i1) - x(i2)*y(i0)
          area(i0) = abs(aa)
        END IF
c
        ! recalculate NEXT POINT AREA
        i0 = inext                                                                 
        i1 = iprev                                                                 
        i2 = next(inext)                                               
        IF (i2 < npt) THEN                                                   
          aa = x(i0)*y(i0) + x(i0)*y(i2) + x(i2)*y(i1) 
     .       - x(i0)*y(i2) - x(i0)*y(i1) - x(i2)*y(i0)
          area(i0) = abs(aa)
        END IF
      END DO   ! ITER
c------------   
      idx = 0
      DO i=1,npt
        IF (area(i) < infinity) THEN
          idx = idx + 1
          x(idx) = x(i)
          y(idx) = y(i)
        END IF
      END DO
      npt = 0
      izero = 0
      DO i = 1,idx
        s = x(i)
        t = y(i)
        IF (s <  zero .and. t >  zero .or.
     .      s >  zero .and. t <  zero .or.
     .      s == zero .and. t /= zero) THEN
          IF (izero == 1) CONTINUE
          s = zero
          t = zero
          izero = 1
        ELSE IF (s < zero .and. t < zero .and.
     .           x(i+1) > zero .and. y(i+1) > zero) THEN ! .or. 
!     .           S > ZERO .and. T > ZERO .and.
!     .           X(I-1) < ZERO .and. Y(I-1) < ZERO .and. I > 1) THEN
          IF (izero == 1) CONTINUE
          s = zero
          t = zero
          izero = 1
        END IF
        npt = npt + 1
        x(npt) = s
        y(npt) = t
      END DO
c-----------
      RETURN
      END