Changeset 2617 in lmdz_wrf

Timestamp:

Jun 19, 2019, 3:53:19 AM (5 years ago)

Author:

lfita

Message:

Adding:

• `crossingpoints_polys': Subroutine to provide the crossing points between two polygons
• starting... `pile_polygons'

Location:

trunk/tools

Files:

• geometry_tools.py (modified) (2 diffs)
• module_scientific.f90 (modified) (2 diffs)

trunk/tools/geometry_tools.py

@@ -17 +17 @@
 import generic_tools as gen
 import numpy.ma as ma
-import module_ForSci as sci
+import module_ForSci as fsci
 
 errormsg = 'ERROR -- error -- ERROR -- error'
@@ -1514 +1514 @@
     return Npts, finalcutpolygon
 
-def cut_between_xpolygon_old(polygon, xval1, xval2, Nadd=20):
-    """ Function to cut a polygon between 2 given value of the x-axis
-      polygon: polygon to cut
-      xval1: first value to use to cut the polygon
-      xval2: first value to use to cut the polygon
-      Nadd: additional points to add to draw the line (20, default)
-    """
-    fname = 'cut_betwen_xpolygon'
-
-    N = polygon.shape[0]
-
-    ipt = None
-    ept = None
-
-    dx = np.zeros((2), dtype=np.float)
-    dy = np.zeros((2), dtype=np.float)
-    icut = np.zeros((2), dtype=int)
-    ecut = np.zeros((2), dtype=int)
-    ipt = np.zeros((2,2), dtype=np.float)
-    ept = np.zeros((2,2), dtype=np.float)
-
-    if xval1 > xval2:
-        print errormsg
-        print '  ' + fname + ': wrong between cut values !!'
-        print '     it is expected xval1 < xval2'
-        print '     values provided xval1: (', xval1, ')> xval2 (', xval2, ')'
-        quit(-1)
-
-    xvals = [xval1, xval2]
-
-    for ic in range(2):
-        xval = xvals[ic]
-        if type(polygon) == type(gen.mamat):
-            # Assuming clockwise polygons
-            for ip in range(N-1):
-                if not polygon.mask[ip,0]:
-                    eep = ip + 1
-                    if eep == N: eep = 0
-      
-                    if (polygon[ip,1] <= xval and polygon[eep,1] > xval) or          \
-                      (polygon[ip,1] < xval and polygon[eep,1] >= xval):
-                        icut[ic] = ip
-                        dx[ic] = polygon[eep,1] - polygon[ip,1]
-                        dy[ic] = polygon[eep,0] - polygon[ip,0]
-                        dd = xval - polygon[ip,1]
-                        ipt[ic,:] = [polygon[ip,0]+dy[ic]*dd/dx[ic], xval]
-
-                    if (polygon[ip,1] >= yval and polygon[eep,1] < xval) or          \
-                      (polygon[ip,1] > yval and polygon[eep,1] <= xval):
-                        ecut[ic] = ip
-                        dx[ic] = polygon[eep,1] - polygon[ip,1]
-                        dy[ic] = polygon[eep,0] - polygon[ip,0]
-                        dd = xval - polygon[ip,1]
-                        ept[ic,:] = [polygon[ip,0]+dy[ic]*dd/dx[ic], xval]
-        else:
-            # Assuming clockwise polygons
-            for ip in range(N-1):
-                eep = ip + 1
-                if eep == N: eep = 0
-      
-                if (polygon[ip,1] <= xval and polygon[eep,1] > xval) or              \
-                  (polygon[ip,1] < xval and polygon[eep,1] >= xval):
-                    icut[ic] = ip
-                    dx[ic] = polygon[eep,1] - polygon[ip,1]
-                    dy[ic] = polygon[eep,0] - polygon[ip,0]
-                    dd = xval - polygon[ip,1]
-                    print 'Lluis ip', ip, 'poly:', polygon[ip,:], 'xval:', xval, 'ip+1', polygon[eep,:]
-                    print '  dx:', dx, 'dy:', dy, 'dd', dd
-                    ipt[ic,:] = [polygon[ip,0]+dy[ic]*dd/dx[ic], xval]
-
-                if (polygon[ip,1] >= xval and polygon[eep,1] < xval) or              \
-                  (polygon[ip,1] > xval and polygon[eep,1] <= xval):
-                    ecut[ic] = ip
-                    dx[ic] = polygon[eep,1] - polygon[ip,1]
-                    dy[ic] = polygon[eep,0] - polygon[ip,0]
-                    dd = xval - polygon[ip,1]
-                    if dx[ic] == 0.:
-                        ept[ic,:] = [polygon[eep,0], xval]
-                    else:
-                        ept[ic,:] = [polygon[ip,0]+dy[ic]*dd/dx[ic], xval]
-
-        if ipt is None or ept is None:
-            print errormsg
-            print '  ' + fname + ': no cutting for polygon at x=', xval, '!!'
-
-    Npts = icut[1] - icut[0] + Nadd + ecut[0] - ecut[1]
-    cutpolygon = np.zeros((Npts,2), dtype=np.float)
-    cutpolygon[0,:] = ipt[0,:]
-    cutpolygon[1:icut[1]-icut[0]+1,:] = polygon[icut[0]+1:icut[1]+1,:]
-    iip = icut[1]-icut[0]
-
-    # cutting lines
-    Nadd2 = int(Nadd/2)
-    cutlines = np.zeros((2,Nadd2,2), dtype=np.float)
-
-    for ic in range(2):
-        print ic, 'Lluis ipt:', ipt[ic,:], 'ept:', ept[ic,:]
-        dx = (ept[ic,1] - ipt[ic,1])/(Nadd2-2)
-        dy = (ept[ic,0] - ipt[ic,0])/(Nadd2-2)
-        print '    dx:', dx, 'dy', dy
-        cutlines[ic,0,:] = ipt[ic,:]
-        for ip in range(1,Nadd2-1):
-            cutlines[ic,ip,:] = ipt[ic,:] + np.array([dy*ip,dx*ip])
-        cutlines[ic,Nadd2-1,:] = ept[ic,:]
-
-    cutpolygon[iip:iip+Nadd2,:] = cutlines[1,:,:]
-    iip = iip + Nadd2
-    cutpolygon[iip:iip+(ecut[0]-ecut[1]),:] = polygon[ecut[1]+1:ecut[0]+1,:]
-    iip = iip + ecut[0]-ecut[1]
-    cutpolygon[iip:iip+Nadd2,:] = cutlines[0,::-1,:]
-
-    cutpolygon = ma.masked_equal(cutpolygon, gen.fillValueF)
-
-    return Npts, cutpolygon
+def pile_polygons(polyns, polygons):
+    """ Function to pile polygons one over the following one
+      polyns: ordered list of polygons. First over all. last below all
+      polygons: dictionary with the polygons
+    >>> pns = ['sqra', 'sqrb']
+    >>> polya = np.array([[-0.5, -0.75], [0.5, -0.75], [0.5, 0.75], [-0.5, 0.75]])
+    >>> polyb = np.array([[-0.75, -0.5], [0.75, -0.5], [0.75, 0.5], [-0.75, 0.5]])
+    >>> plgs = {'sqra': polya, 'sqrb': polyb}
+    >>> pile_polygons(pns, plgs)
+
+    """
+    fname = 'pile_polygons'
+
+    Npolys = len(polygons)
+    ipoly = polygons[polyns[Npolys-1]]
+    iNpts = ipoly.shape[0]
+    pilepolygons = ipoly
+
+    pilesortpolys = polyns[::-1]
+    pilesortpolys.remove(polyns[Npolys-1])
+    for polyn in pilesortpolys:
+        poly = polygons[polyn]
+        Npts = poly.shape[0]
+        print ' LLuis shapes ipoly', ipoly.shape, 'poly', poly.shape
+        Nint, inti, intp, pts = fsci.module_scientific.crossingpoints_polys(         \
+          nvertexa=iNpts, nvertexb=Npts, nvertexab=iNpts*Npts, polya=ipoly,          \
+          polyb=poly)
+        print Nint
+        for ip in range(Nint):
+            print ip, ':', inti[ip], intp[ip], '=', pts[ip]
+       
+    return pilepolygons
+
+#pns = ['sqra', 'sqrb']
+#polya = np.array([[-0.5, -0.75], [0.5, -0.75], [0.5, 0.75], [-0.5, 0.75]])
+#polyb = np.array([[-0.75, -0.5], [0.75, -0.5], [0.75, 0.5], [-0.75, 0.5]])
+#plgs = {'sqra': polya, 'sqrb': polyb}
+#pile_polygons(pns, plgs)
+
+
 
 ####### ###### ##### #### ### ## #

trunk/tools/module_scientific.f90

@@ -17 +17 @@
 ! clean_polygons: Subroutine to clean polygons from non-used paths, polygons only left as path since they are inner path of a hole
 ! coincident_polygon: Subroutine to provide the intersection polygon between two polygons
+! crossingpoints_polys: Subroutine to provide the crossing points between two polygons
 ! distanceRK: Function to provide the distance between two points
 ! FindMinimumR_K*: Function returns the location of the minimum in the section between Start and End.
@@ -7507 +7508 @@
   END SUBROUTINE multi_spaceweightcount_in3DRK3_slc3v3
 
+  SUBROUTINE crossingpoints_polys(NvertexA, NvertexB, NvertexAB, polyA, polyB, Ncross, crossA, crossB,&
+    crosspoints)
+  ! Subroutine to provide the crossing points between two polygons
+
+    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)                                 :: Ncross
+    INTEGER, DIMENSION(NvertexAB), INTENT(out)           :: crossA, crossB
+    REAL(r_k), DIMENSION(NvertexAB,2), INTENT(out)       :: crosspoints
+
+! Local
+    INTEGER                                              :: iA, iB, icoin, ii
+    REAL(r_k), DIMENSION(2,2)                            :: face1, face2
+    INTEGER                                              :: intersect
+    REAL(r_k), DIMENSION(2)                              :: intersectpt
+
+
+!!!!!!! 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)
+! Ncross: number of crossing vertices
+! crossA: vertex of polygon A where there is a crossing
+! crossB: vertex of polygon B where there is a crossing
+! crosspoints: coordinates of the crossing points
+
+    fname = 'crossingpoints_polys'
+
+    Ncross = 0
+    crossA = 0
+    crossB = 0
+    crosspoints = 0.
+
+    ! Looping
+    DO iA=1, NvertexA-1
+      face1(1,:) = polyA(iA,:)
+      face1(2,:) = polyA(iA+1,:)
+      DO iB=1, NvertexB-1
+        face2(1,:) = polyB(iB,:)
+        face2(2,:) = polyB(iB+1,:)
+        CALL intersectfaces(face1, face2, intersect, intersectpt)
+        IF (intersect /= 0) THEN
+          Ncross = Ncross + 1
+          crossA(Ncross) = iA
+          crossB(Ncross) = iB
+          crosspoints(Ncross,:) = intersectpt
+        END IF
+      END DO
+      ! Last face B
+      iB = NvertexB
+      face2(1,:) = polyB(iB,:)
+      face2(2,:) = polyB(1,:)
+      CALL intersectfaces(face1, face2, intersect, intersectpt)
+      IF (intersect /= 0) THEN
+        Ncross = Ncross + 1
+        crossA(Ncross) = iA
+        crossB(Ncross) = iB
+        crosspoints(Ncross,:) = intersectpt
+      END IF
+    END DO
+
+    ! Last face A
+    iA = NvertexA
+    face1(1,:) = polyA(iA,:)
+    face1(2,:) = polyA(1,:)
+    DO iB=1, NvertexB-1
+      face2(1,:) = polyB(iB,:)
+      face2(2,:) = polyB(iB+1,:)
+      CALL intersectfaces(face1, face2, intersect, intersectpt)
+      IF (intersect /= 0) THEN
+        Ncross = Ncross + 1
+        crossA(Ncross) = iA
+        crossB(Ncross) = iB
+        crosspoints(Ncross,:) = intersectpt
+      END IF
+    END DO
+    ! Last face B
+    iB = NvertexB
+    face2(1,:) = polyB(iB,:)
+    face2(2,:) = polyB(1,:)
+    CALL intersectfaces(face1, face2, intersect, intersectpt)
+    IF (intersect /= 0) THEN
+      Ncross = Ncross + 1
+      crossA(Ncross) = iA
+      crossB(Ncross) = iB
+      crosspoints(Ncross,:) = intersectpt
+    END IF
+
+  END SUBROUTINE crossingpoints_polys
+
 END MODULE module_scientific