From 0abae690a217864084861681c7c6d251214baf88 Mon Sep 17 00:00:00 2001
From: Paul RASCLE <paul.rascle@openfields.fr>
Date: Tue, 14 Jul 2020 11:33:27 +0200
Subject: [PATCH] split the free border by the exension shape

---
 src/HYDROTools/shapesGroups.py | 42 ++++++++++++++++++++++++++++++++--
 1 file changed, 40 insertions(+), 2 deletions(-)

diff --git a/src/HYDROTools/shapesGroups.py b/src/HYDROTools/shapesGroups.py
index 85b8a28b..cac3e446 100644
--- a/src/HYDROTools/shapesGroups.py
+++ b/src/HYDROTools/shapesGroups.py
@@ -164,6 +164,9 @@ def exploreEdgeGroups(meshFile, offsetX=0., offsetY=0.):
 
 def fitShapePointsToMesh(freeBorderShapefile, shapefileToAdjust):
     print(" === fitShapePointsToMesh", freeBorderShapefile, shapefileToAdjust)
+
+    # --- find domain freeBorder:  bounding box englobing all others
+
     fb = shapefile.Reader(freeBorderShapefile)
     fbShapes = fb.shapes()
     maxbbox=[1.e30, 1.e30, -1.e30, -1.e30]
@@ -175,6 +178,9 @@ def fitShapePointsToMesh(freeBorderShapefile, shapefileToAdjust):
             outerBboxIndex = i
     fbs = fbShapes[outerBboxIndex]
 
+    # --- find the intersections of the shapefile to adjust and the domain free border:
+    #     the closests points (two pairs of points)
+
     sf = shapefile.Reader(shapefileToAdjust)
     shapes = sf.shapes()
     sfta = sf.shape(0)
@@ -199,6 +205,9 @@ def fitShapePointsToMesh(freeBorderShapefile, shapefileToAdjust):
         ifb1 = pdist[1][1]
         ifb2 = pdist[0][1]
 
+    # --- write the adusted shapefile: free border closest points replaced with corresponding points
+    #     on the free border. two polylines, one inside the domain, one outside
+
     a = os.path.splitext(shapefileToAdjust)
     shapefileAdjusted = a[0] + '_adj' + a[1]
     chainName = os.path.basename(a[0])
@@ -217,8 +226,9 @@ def fitShapePointsToMesh(freeBorderShapefile, shapefileToAdjust):
 
     chaincoords = []
     chaincoords.append(fbs.points[ifb2])
-    for i in range(i2+1, len(sfta.points)):
-        chaincoords.append(sfta.points[i])
+    if i2+1 < len(sfta.points):
+        for i in range(i2+1, len(sfta.points)):
+            chaincoords.append(sfta.points[i])
     for i in range(i1):
         chaincoords.append(sfta.points[i])
     chaincoords.append(fbs.points[ifb1])
@@ -226,3 +236,31 @@ def fitShapePointsToMesh(freeBorderShapefile, shapefileToAdjust):
     w.record(chainName + '_1')
     w.close()
 
+    # write the free border splited in two polylines (cut by the adjusted shapefile)
+
+    a = os.path.splitext(freeBorderShapefile)
+    freeBorderSplit = a[0] + '_split' + a[1]
+    chainName = os.path.basename(a[0])
+
+    w = shapefile.Writer(freeBorderSplit)
+    w.shapeType = 3
+    w.field('name', 'C')
+
+    if (ifb1 > ifb2):
+        i = ifb1; ifb1 = ifb2; ifb2 = i
+
+    chaincoords = []
+    for i in range(ifb1, ifb2+1):
+        chaincoords.append(fbs.points[i])
+    w.line([chaincoords])
+    w.record(chainName + '_0')
+
+    chaincoords = []
+    for i in range(ifb2, len(fbs.points)):
+        chaincoords.append(fbs.points[i])
+    for i in range(ifb1+1):
+        chaincoords.append(fbs.points[i])
+    w.line([chaincoords])
+    w.record(chainName + '_1')
+    w.close()
+
-- 
2.39.2