Changeset 2624 in lmdz_wrf for trunk

Timestamp:

Jun 22, 2019, 9:56:59 PM (5 years ago)

Author:

lfita

Message:

Working on piling polygons

File:

• trunk/tools/geometry_tools.py (modified) (4 diffs)

trunk/tools/geometry_tools.py

@@ -292 +292 @@
     1.57079632679
     >>> angle_vectors2D([0.,1.], [1.,0.])
-    1.57079632679
+    -1.57079632679
     """
     fname = 'angle_vectors2D'
@@ -310 +310 @@
 
     return theta
-
-print angle_vectors2D([1.,0.], [0.,1.])
-print angle_vectors2D([0.,1.], [1.,0.])
-
-quit()
 
 def max_coords_poly(polygon):
@@ -1577 +1572 @@
     """
     fname = 'pile_polygons'
-
-    Npolys = len(polygons)
+    pilepolygons = dict(polygons)
+
+    Npolys = len(polyns)
     ipoly = polygons[polyns[Npolys-1]]
     iNpts = ipoly.shape[0]
-    pilepolygons = ipoly
 
     pilesortpolys = polyns[::-1]
     pilesortpolys.remove(polyns[Npolys-1])
+    newipoly = []
     for polyn in pilesortpolys:
         poly = polygons[polyn]
@@ -1591 +1587 @@
           nvertexa=iNpts, nvertexb=Npts, nvertexab=iNpts*Npts, polya=ipoly,          \
           polyb=poly)
-        iip=0
-        for ip in Nint:
-            newipoly = ipoly[iip:inti]
-        
+        Npi=0
+        for ip in range(Nint):
+            # Determine side of the point of the polygon to cut
+            #   Here are assumed that all polygons are clockwise created, thus any
+            #   cutting point will lay to the left side of the face which cuts
+            #   but we need to deterime if it is the starting point or the ending point
+            #   the one which lays to the left of the over pile polygon
+            iip = inti[ip]
+            iip1 = iip + 1
+            if iip1 > iNpts-1: iip1=0
+            pip = intp[ip]
+            pip1 = pip + 1
+            if pip1 > Npts-1: pip1=0
+
+            iface = ipoly[iip,:]-poly[pip,:]
+            pface = poly[pip1,:]-poly[pip,:]
+            ifangle=angle_vectors2D(pface,iface)
+
+            iface = ipoly[iip1,:]-poly[pip,:]
+            pface = poly[pip1,:]-poly[pip,:]
+            efangle=angle_vectors2D(pface,iface)
+
+            if ifangle < 0.: 
+                # Initial point of the cuted face to the left of the overlying polygon
+                for iv in range(Npi,Npi+iip): newipoly.append(ipoly[iv,:])
+                Npi = Npi + iip
+                newipoly.append([gen.fillvalueF, gen.fillvalueF])
+
+            elif efangle < 0.:
+                # Ending point of the cuted face to the left of the overlying polygon
+                newipoly.append([gen.fillvalueF, gen.fillvalueF])
+                newipoly.append(ipoly[iip1,:])
+                Npi = iip1
+            else:
+                print errormsg
+                print '  ' + fname + ': No left point found !!'
+                print '    for pile face'
+    
+        pilepolygons[polyns[Npolys-1]] = newipoly
        
     return pilepolygons