i22clip_tools.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/.
!||====================================================================
!||    crossprod2d            ../engine/source/interfaces/int22/i22clip_tools.F
!||--- called by ------------------------------------------------------
!||    intersectp             ../engine/source/interfaces/int22/i22clip_tools.F
!||    is_on_1st_half_plane   ../engine/source/interfaces/int22/i22clip_tools.F
!||====================================================================
        function crossprod2d( v1, v2)
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-----------------------------------------------
          my_real         , intent(inout) :: v1(2),v2(2)
          my_real                      :: crossprod2d
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------  
C-----------------------------------------------
C   S o u r c e   L i n e s
C-----------------------------------------------         
          crossprod2d = v1(1)*v2(2) - v1(2)*v2(1)
          return
        end function crossprod2d
C----------------------------------------------- 
 
 
 
!||====================================================================
!||    polygonalclipping   ../engine/source/interfaces/int22/i22clip_tools.F
!||--- called by ------------------------------------------------------
!||    i22wetsurf          ../engine/source/interfaces/int22/i22wetsurf.F
!||--- calls      -----------------------------------------------------
!||    clipedge            ../engine/source/interfaces/int22/i22clip_tools.F
!||====================================================================
        subroutine polygonalclipping( basePolygon, clipPolygon, resultPolygon )
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"        
          ! based on Sutherland Hodgman  algorithm
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
          type polygon
            my_real          :: node(16,2)
            integer :: NumNodes          
          end type polygon
          type(polygon)    :: basePolygon, clipPolygon, resultPolygon
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
          type(polygon)    :: tmpPolygon
          my_real          :: edge_p1(2),edge_p2(2)    ! vertices of edge to clipPolygon tmpPolygon
          integer :: i  
          tmppolygon%NumNodes = clippolygon%NumNodes
          IF(tmppolygon%NumNodes > 0)THEN
            tmppolygon%node(1:tmppolygon%NumNodes,:) = clippolygon%node(1:tmppolygon%NumNodes,:)
          END IF
          
          do i=1,basepolygon%NumNodes-1 ! for each edge i of the polygon basePolygon
            edge_p1(:) = basepolygon%node(i,:)   !  node 1 of edge i
            edge_p2(:) = basepolygon%node(i+1,:) !  node 2 of edge i
            ! clipPolygon the work polygon by edge i
            call clipedge( tmppolygon, edge_p1, edge_p2, resultpolygon)
            ! tmpPolygon <= resultPolygon
            tmppolygon%NumNodes = resultpolygon%NumNodes
            tmppolygon%node(1:tmppolygon%NumNodes,:) = resultpolygon%node(1:tmppolygon%NumNodes,:)
          end do
        end subroutine polygonalclipping
C----------------------------------------------- 
 
 
 
!||====================================================================
!||    clipedge               ../engine/source/interfaces/int22/i22clip_tools.F
!||--- called by ------------------------------------------------------
!||    polygonalclipping      ../engine/source/interfaces/int22/i22clip_tools.F
!||--- calls      -----------------------------------------------------
!||    intersectp             ../engine/source/interfaces/int22/i22clip_tools.F
!||    is_on_1st_half_plane   ../engine/source/interfaces/int22/i22clip_tools.F
!||====================================================================
        subroutine clipedge( polyg, p1, p2, resultPolyg )
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-----------------------------------------------
          type polygon
            my_real          :: node(16,2)
            integer :: NumNodes          
          end type polygon
          type(polygon) :: polyg, resultPolyg
          my_real         , intent(inout) :: p1(2), p2(2) 
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------         
          my_real                      :: x1(2), x2(2), intersecpoint(2)
          integer  ::  i, k
          INTERFACE
            function intersectp(a,b,c,d,iflg)
              my_real         , intent(inout) :: a(2),b(2),c(2),d(2)
              integer, intent(in)             :: iflg
              my_real          intersectp(2)
            end function
            function is_on_1st_half_plane(a,b,c)
              my_real         , intent(inout) :: a(2),b(2),c(2)
              logical                         :: IS_ON_1ST_HALF_PLANE
            end function
          END INTERFACE
C-----------------------------------------------
C   S o u r c e   L i n e s
C-----------------------------------------------           
          k = 0
          do i=1,polyg%NumNodes-1 ! for each edge i of poly
            x1(:) = polyg%node(i,:)   ! node 1 of edge i
            x2(:) = polyg%node(i+1,:) ! node 2 of edge i
 
            if ( is_on_1st_half_plane(x1, p1, p2) ) then 
              if ( is_on_1st_half_plane(x2, p1, p2) ) then
                ! add the node 2 to the output polygon
                k = k+1
                resultpolyg%node(k,:) = x2(:)
 
              else ! node i+1 is outside
                intersecpoint = intersectp(x1, x2, p1,p2  ,1)
                k = k+1
                resultpolyg%node(k,:) = intersecpoint(:)
              end if
            else ! node i is outside
              if ( is_on_1st_half_plane(x2, p1, p2) ) then
                intersecpoint = intersectp(x1, x2, p1,p2  ,1)
                k = k+1
                resultpolyg%node(k,:) = intersecpoint(:)
 
                k = k+1
                resultpolyg%node(k,:) = x2(:)
              end if
            end if
          end do
          if (k > 0) then
            ! if the last vertice is not equal to the first one
            if ( (resultpolyg%node(1,1) /= resultpolyg%node(k,1)) .or.  (resultpolyg%node(1,2) /= resultpolyg%node(k,2))) then
              k=k+1
              resultpolyg%node(k,:) = resultpolyg%node(1,:)
            end if
          end if
          ! set the size of the resultPolyggon
          resultpolyg%NumNodes = k
        end subroutine clipedge
C-----------------------------------------------         
 
 
 
!||====================================================================
!||    intersectp             ../engine/source/interfaces/int22/i22clip_tools.F
!||--- called by ------------------------------------------------------
!||    clipedge               ../engine/source/interfaces/int22/i22clip_tools.F
!||    getprojectedfacetype   ../engine/source/interfaces/int22/i22wetsurf.F
!||--- calls      -----------------------------------------------------
!||    crossprod2d            ../engine/source/interfaces/int22/i22clip_tools.f
!||====================================================================
        function intersectp( Seg_P1, Seg_P2, Line_P1, Line_P2, IFLG)
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-----------------------------------------------
          my_real         ,intent(inout) :: seg_p1(2), seg_p2(2)
          my_real         ,intent(inout) :: line_p1(2), line_p2(2)
          INTEGER         ,intent(in)    :: iflg
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
          my_real                     :: intersectp(2), v1(2), v2(2), seg_p1_line_p1(2)
          my_real          :: alpha, tmp,tmp2, l
          my_real          :: crossprod2d, tol
          external crossprod2d
C-----------------------------------------------
C   S o u r c e   L i n e s
C-----------------------------------------------
          tol = em06
 
          intersectp(1:2) = ep30
          v1(:) = seg_p2(:) - seg_p1(:)
          v2(:) = line_p2(:) - line_p1(:)
          l=max(sum(v1(:)*v1(:)),sum(v2(:)*v2(:)))
          tmp =  crossprod2d(v1,v2)
          if (abs(tmp) <= tol*l) then
!          if (tmp == 0.0D00) then
            !colinear vectors
            seg_p1_line_p1(:) = line_p1(:) - seg_p1(:)
            ! if the the segment [Seg_P1Seg_P2] is included in the line (Line_P1Line_P2)
            tmp =  crossprod2d(seg_p1_line_p1,v1)
!            if ( tmp == 0.0d00) then
            if(abs(tmp) <= tol*l)then
              ! the intersection is the last point of the segment
              intersectp(:) = seg_p2(:)
            end if
          else
            ! non colinear vectors
            seg_p1_line_p1(:) = line_p1(:) - seg_p1(:)
            ! parametric coordinates
            tmp2 = crossprod2d(seg_p1_line_p1,v2)
            alpha = tmp2
            tmp2 = crossprod2d(v1,v2)
            alpha = alpha/tmp2
            ! if a is not in [0;1]
            if ( (alpha >= -em03) .and. (alpha <= one+em03)) then
              alpha=max(zero,min(alpha,one))
              intersectp(:) = seg_p1(:) + alpha*v1(:)
            end if
          end if
 
        end function intersectp
C-----------------------------------------------
 
 
 
!||====================================================================
!||    is_on_1st_half_plane   ../engine/source/interfaces/int22/i22clip_tools.F
!||--- called by ------------------------------------------------------
!||    clipedge               ../engine/source/interfaces/int22/i22clip_tools.F
!||--- calls      -----------------------------------------------------
!||    crossprod2d            ../engine/source/interfaces/int22/i22clip_tools.F
!||====================================================================
        function is_on_1st_half_plane( point, p1, p2)
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-----------------------------------------------
          my_real         , intent(inout) :: point(2), p1(2), p2(2)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
          my_real                         :: v1(2), v2 (2)
          logical :: is_on_1st_half_plane
          my_real          :: crossp, l
          my_real          :: crossprod2d
          EXTERNAL         :: crossprod2d
C-----------------------------------------------
C   S o u r c e   L i n e s
C-----------------------------------------------
          v1(:) = p2(:) -  p1(:)
          v2(:) = point(:)  -  p1(:)
          l=max(sum(v1(:)*v1(:)),sum(v2(:)*v2(:)))
          crossp = crossprod2d(v1,v2)
          !DONT BE TOO MUCH TOLERANT OTHERWISE ALGO WILL ADD AN EMPTY INTERSECTION POINT WHICH WOULD RECQUIRE ADDITIONAL CODE TO BE SKIPPED.
          !EM06 : issue
          !EM10 : OK
          if (  crossp >= -em10*l) then
            is_on_1st_half_plane = .true.
          else
            is_on_1st_half_plane = .false.
          end if
        end function is_on_1st_half_plane
C-----------------------------------------------
 
 
 
!||====================================================================
!||    setclockwisepolyg   ../engine/source/interfaces/int22/i22clip_tools.F
!||====================================================================
        subroutine setclockwisepolyg( Polyg, Area )
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-----------------------------------------------
          type polygon
            my_real          :: node(16,2)
            integer :: NumNodes          
          end type polygon
          type(polygon)   ,intent(inout)  :: Polyg
          my_real         , intent(inout) :: area 
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------         
          my_real :: x(16),y(16),total
          integer ::  i,n
C-----------------------------------------------
C   S o u r c e   L i n e s
C-----------------------------------------------           
 
          total = zero
          n     = polyg%NumNodes
          
          DO i=1,n
            x(i) = polyg%node(i,1)            
            y(i) = polyg%node(i,2)                        
          ENDDO
          
          DO i=1,n-1
            total = total + (x(i+1)-x(i))*(y(i+1)+y(i))
          ENDDO
          
          area = half * abs(total)
          
          !setclockwise
          IF (total < zero) THEN
            DO i=1,n
              polyg%node(i,1) = x(n-i+1)           
              polyg%node(i,2) = y(n-i+1)                       
            ENDDO 
          ENDIF   
          
      END SUBROUTINE
C----------------------------------------------- 
 
 
 
!||====================================================================
!||    setcounterclockwisepolyg   ../engine/source/interfaces/int22/i22clip_tools.F
!||--- called by ------------------------------------------------------
!||    i22wetsurf                 ../engine/source/interfaces/int22/i22wetsurf.F
!||====================================================================
        subroutine setcounterclockwisepolyg( Polyg, Area )
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-----------------------------------------------
          type polygon
            my_real          :: node(16,2)
            integer :: NumNodes          
          end type polygon
          type(polygon)   ,intent(inout)  :: Polyg
          my_real         , intent(inout) :: area 
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------         
          my_real                      :: x(16),y(16),total
          integer  ::  i,n
C-----------------------------------------------
C   S o u r c e   L i n e s
C-----------------------------------------------           
 
          total = zero
          n     = polyg%NumNodes
          
          DO i=1,n
            x(i) = polyg%node(i,1)            
            y(i) = polyg%node(i,2)                        
          ENDDO
          
          DO i=1,n-1
            total = total + (x(i+1)-x(i))*(y(i+1)+y(i))
          ENDDO
          
          area = half * abs(total)
          
          !setclockwise
          IF (total > zero) THEN
            DO i=1,n
              polyg%node(i,1) = x(n-i+1)           
              polyg%node(i,2) = y(n-i+1)                       
            ENDDO 
          ENDIF   
          
      END SUBROUTINE
C-----------------------------------------------