Changeset 2263 in lmdz_wrf

Timestamp:

Dec 20, 2018, 1:54:00 PM (6 years ago)

Author:

lfita

Message:

Adding final functions for the coincident polygons !!

File:

• trunk/tools/module_scientific.f90 (modified) (10 diffs)

trunk/tools/module_scientific.f90

@@ -30 +30 @@
 ! path_properties: Subroutine to determine the properties of a path
 ! percentiles_R_K*D: Subroutine to compute the percentiles of a *D R_K array along given set of axis
+! point_in_face: Function to determine if a given point is on a face of a polygon
 ! point_inside: Function to determine if a given point is inside a polygon providing its sorted vertices
 ! poly_has_point: Function to determine if a polygon has already a given point as one of its vertex
@@ -52 +53 @@
 ! runmean_F1D: Subroutine fo computing the running mean of a given set of float 1D values
 ! shoelace_area_polygon: Computing the area of a polygon using sholace formula
-! sort_polygon: Subroutine to sort a polygon using its center as average of the coordinates
+! sort_polygon: Subroutine to sort a polygon using its center as average of the coordinates and remove duplicates
 ! SortR_K*: Subroutine receives an array x() r_K and sorts it into ascending order.
 ! StatsR_K: Subroutine to provide the minmum, maximum, mean, the quadratic mean, and the standard deviation of a 
@@ -4413 +4414 @@
     ! First, include that vertex which do not lay within any polygon
     DO iA=1, NvertexA
-      PRINT *, '  iA:', iA, ':', polyA(iA,:), point_inside(polyA(iA,:), NvertexB, polyB)
+      !PRINT *, '  iA:', iA, ':', polyA(iA,:), point_inside(polyA(iA,:), NvertexB, polyB)
       IF (.NOT. point_inside(polyA(iA,:), NvertexB, polyB)) THEN
         icoin = icoin + 1
@@ -4421 +4422 @@
 
     DO iB=1, NvertexB
-      PRINT *, '  iB:', iB, ':', polyB(iB,:), point_inside(polyB(iB,:), NvertexA, polyA)
+      !PRINT *, '  iB:', iB, ':', polyB(iB,:), point_inside(polyB(iB,:), NvertexA, polyA)
       IF (.NOT. point_inside(polyB(iB,:), NvertexA, polyA)) THEN
         icoin = icoin + 1
@@ -4448 +4449 @@
         ! Compute areas of the four possible triangles. Introduce the coincident vertexs not included
         CALL intersectfaces(face1, face2, intersct, ptintersct)
-        PRINT *,iA,':',face1(1,:),';',face1(2,:), '=', iB, face2(1,:),';',face2(2,:), '<>', intersct,':', ptintersct
+        !PRINT *,iA,':',face1(1,:),';',face1(2,:), '=', iB, face2(1,:),';',face2(2,:), '<>', intersct,':', ptintersct
         IF (intersct == 1) THEN
           IF (.NOT.poly_has_point(icoin,coinpoly(1:icoin,:),ptintersct) ) THEN
@@ -4467 +4468 @@
   END SUBROUTINE join_polygon
 
-  SUBROUTINE sort_polygon(Nvertex, polygon, sortpoly)
-  ! Subroutine to sort a polygon using its center as average of the coordinates
+  SUBROUTINE sort_polygon(Nvertex, polygon, sense, Nnewvertex, newpoly)
+  ! Subroutine to sort a polygon using its center as average of the coordinates and remove duplicates
   !    Should be used the centroid instead, but by now let do it simple
   !    https://en.wikipedia.org/wiki/Centroid
@@ -4475 +4476 @@
     IMPLICIT NONE
 
-    INTEGER, INTENT(in)                                  :: Nvertex
+    INTEGER, INTENT(in)                                  :: Nvertex, sense
     REAL(r_k), DIMENSION(Nvertex,2), INTENT(in)          :: polygon
-    REAL(r_k), DIMENSION(Nvertex,2), INTENT(out)         :: sortpoly
+    INTEGER, INTENT(out)                                 :: Nnewvertex
+    REAL(r_k), DIMENSION(Nvertex,2), INTENT(out)         :: newpoly
 
 ! Local
@@ -4484 +4486 @@
     REAL(r_k), DIMENSION(2)                              :: center
     REAL(r_k), DIMENSION(Nvertex)                        :: angles
+    REAL(r_k), DIMENSION(Nvertex,2)                      :: sortpoly
 
 !!!!!!! Variables
-! Nvertex: number of vetexs
+! Nvertex: number of vertices
 ! polygon: coordinates of the vertices of the polygon
-! sortpoly: sorted polygon (clockwise)
+! sense: sens of sorting thepolygon (1: clockwise, -1: anti-clockwise)
+! sortpoly: sorted polygon
+! Nnewvertex: number of vertices new polygon
+! newpoly: sorted and duplicate removed polygon
 
     fname = 'sort_polygon'
@@ -4506 +4512 @@
         vang = ATAN2(polygon(j,2)-center(2), polygon(j,1)-center(1))
         IF (angles(iv) == vang) THEN
-          sortpoly(iv,:) = polygon(j,:)
+          IF (sense == -1) THEN
+            sortpoly(iv,:) = polygon(j,:)
+          ELSE
+            sortpoly(Nvertex-iv+1,:) = polygon(j,:)
+          END IF
           EXIT
         END IF
       END DO
     END DO
+
+    newpoly(1,:) = sortpoly(1,:)
+    j = 1
+    DO iv=2, Nvertex
+      IF (.NOT.poly_has_point(j,newpoly(1:j,:),sortpoly(iv,:)) ) THEN
+        j = j+1
+        newpoly(j,:) = sortpoly(iv,:)
+      END IF
+    END DO
+    Nnewvertex = j
 
   END SUBROUTINE sort_polygon
@@ -4569 +4589 @@
   END FUNCTION point_inside
 
+  LOGICAL FUNCTION point_in_face(pt, Nvertex, poly)
+  ! Function to determine if a given point is on a face of a polygon
+
+    IMPLICIT NONE
+
+    REAL(r_k), DIMENSION(2), INTENT(in)                  :: pt
+    INTEGER, INTENT(in)                                  :: Nvertex
+    REAL(r_k), DIMENSION(Nvertex,2), INTENT(in)          :: poly
+! Local
+    INTEGER                                              :: iv
+    REAL(r_k)                                            :: ix, ex, iy, ey, tmpv
+    REAL(r_k)                                            :: dx, dy, A, B
+
+!!!!!!! Variables
+! pt: point to look for
+! Nvertex: Number of vertices of the polygon
+! poly: polygon
+    fname = 'point_in_face'
+
+    point_in_face = .FALSE.
+    DO iv=1, Nvertex
+      IF (iv < Nvertex) THEN
+        ix = poly(iv,1)
+        ex = poly(iv+1,1)
+        iy = poly(iv,2)
+        ey = poly(iv+1,2)
+      ELSE
+        ix = poly(iv,1)
+        ex = poly(1,1)
+        iy = poly(iv,2)
+        ey = poly(1,2)
+      END IF    
+      dx = ex - ix
+      dy = ey - iy
+
+      IF (dx == zeroRK) THEN
+        IF (pt(1) == ix) THEN
+          IF ( (iy < ey) .AND. (pt(2) >= iy) .AND. pt(2) <= ey) THEN
+            point_in_face = .TRUE.
+            EXIT
+          ELSE IF ( (iy > ey) .AND. (pt(2) >= ey) .AND. pt(2) <= iy) THEN
+            point_in_face = .TRUE.
+            EXIT
+          END IF
+        END IF
+      ELSE
+        IF (dy == zeroRK) THEN
+          IF (pt(2) == iy) THEN
+            IF ((ix < ex) .AND. (pt(1) >= ix) .AND. pt(1) <= ex) THEN
+              point_in_face = .TRUE.
+              EXIT            
+            ELSE IF ((ix > ex) .AND. (pt(1) >= ex) .AND. pt(1) <= ix) THEN
+              point_in_face = .TRUE.
+              EXIT
+            END IF
+          END IF
+        ELSE
+          A = iy
+          B = (ey-iy)/(ex-ix)
+          IF (A+B*(pt(1)-ix) == pt(2)) THEN
+            point_in_face = .TRUE.
+            EXIT
+          END IF
+        END IF
+      END IF
+    END DO
+
+  END FUNCTION point_in_face
+
+  SUBROUTINE coincident_polygon(NvertexA, NvertexB, NvertexAB, polyA, polyB, Ncoinvertex, coinpoly)
+  ! Subroutine to provide the intersection polygon between two polygons
+  ! AFTER: http://www.cap-lore.com/MathPhys/IP/ and http://www.cap-lore.com/MathPhys/IP/IL.html
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(in)                                  :: NvertexA, NvertexB, NvertexAB
+    REAL(r_k), DIMENSION(NvertexA,2), INTENT(in)         :: polyA
+    REAL(r_k), DIMENSION(NvertexB,2), INTENT(in)         :: polyB
+    INTEGER, INTENT(out)                                 :: Ncoinvertex
+    REAL(r_k), DIMENSION(NvertexAB,2), INTENT(out)       :: coinpoly
+
+! Local
+    INTEGER                                              :: iA, iB, icoin, ii
+    REAL(r_k), DIMENSION(2,2)                            :: face1, face2
+    INTEGER                                              :: intersct
+    REAL(r_k), DIMENSION(2)                              :: ptintersct
+
+!!!!!!! variables
+! NvertexA: number of vertexs polygon A
+! NvertexB: number of vertexs polygon B
+! polyA: pairs of coordinates for the polygon A (clockwise)
+! polyB: pairs of coordinates for the polygon B (clockwise)
+! Ncoinvertex: number of vertexes for the coincident polygon
+! coinpoly: pairs of coordinates for the coincident polygon (clockwise)
+
+    fname = 'coincident_polygon'
+
+    icoin = 0
+    coinpoly = 0.
+    ! First, include that vertex which lay within any polygon
+    DO iA=1, NvertexA
+      IF (point_inside(polyA(iA,:), NvertexB, polyB)) THEN
+        icoin = icoin + 1
+        coinpoly(icoin,:) = polyA(iA,:)
+      END IF
+      IF (point_in_face(polyA(iA,:), NvertexB, polyB)) THEN
+        icoin = icoin + 1
+        coinpoly(icoin,:) = polyA(iA,:)
+      END IF
+    END DO
+
+    DO iB=1, NvertexB
+      IF (point_inside(polyB(iB,:), NvertexA, polyA)) THEN
+        icoin = icoin + 1
+        coinpoly(icoin,:) = polyB(iB,:)
+      END IF
+      IF (point_in_face(polyB(iB,:), NvertexA, polyA)) THEN
+        icoin = icoin + 1
+        coinpoly(icoin,:) = polyB(iB,:)
+      END IF
+    END DO
+
+    ! Look interesections
+    DO iA=1, NvertexA
+      ! Getting couple of vertexs from polyA and polyB
+      IF (iA /= NvertexA) THEN
+        face1(1,:) = polyA(iA,:)
+        face1(2,:) = polyA(iA+1,:)
+      ELSE
+        face1(1,:) = polyA(iA,:)
+        face1(2,:) = polyA(1,:)
+      END IF
+      DO iB=1, NvertexB
+        IF (iB /= NvertexB) THEN
+          face2(1,:) = polyB(iB,:)
+          face2(2,:) = polyB(iB+1,:)
+        ELSE
+          face2(1,:) = polyB(iB,:)
+          face2(2,:) = polyB(1,:)
+        END IF
+        
+        ! Compute areas of the four possible triangles. Introduce the coincident vertexs not included
+        CALL intersectfaces(face1, face2, intersct, ptintersct)
+        !PRINT *,iA,':',face1(1,:),';',face1(2,:), '=', iB, face2(1,:),';',face2(2,:), '<>', intersct,':', ptintersct
+        IF ((intersct /= 0) .AND. (.NOT.poly_has_point(icoin,coinpoly(1:icoin,:),ptintersct)) ) THEN
+          icoin = icoin + 1
+          coinpoly(icoin,:) = ptintersct
+        END IF
+
+      END DO
+    END DO
+    Ncoinvertex = icoin
+
+  END SUBROUTINE coincident_polygon
+
 END MODULE module_scientific