From b71b76b1129bcf57148be84c076404a9715b106d Mon Sep 17 00:00:00 2001
From: Paul RASCLE <paul.rascle@openfields.fr>
Date: Fri, 10 Jul 2020 12:15:45 +0200
Subject: [PATCH] shape adjustement on free border mesh points

---
 src/HYDROTools/shapesGroups.py | 77 +++++++++++++++++++++++++++++++---
 1 file changed, 71 insertions(+), 6 deletions(-)

diff --git a/src/HYDROTools/shapesGroups.py b/src/HYDROTools/shapesGroups.py
index 14f6e8ee..85b8a28b 100644
--- a/src/HYDROTools/shapesGroups.py
+++ b/src/HYDROTools/shapesGroups.py
@@ -1,7 +1,4 @@
 
-
-# -----------------------------------------------------------------------------
-
 import salome
 salome.salome_init()
 
@@ -13,6 +10,9 @@ import MEDLoader as ml
 import medcoupling as mc
 
 import shapefile
+import math
+import os
+
 
 def freeBordersGroup(meshFile):
     print(" === freeBordersGroup", meshFile)
@@ -35,6 +35,7 @@ def freeBordersGroup(meshFile):
         print('ExportMED() failed. Invalid file name?')
     return newMeshName
 
+
 def explodeGroup(grp, grpName):
     print(" === explodeGroup", grpName)
     nbCells=grp.getNumberOfCells()
@@ -106,6 +107,7 @@ def explodeGroup(grp, grpName):
         print(chainDesc[1:])
     return nodeChains
 
+
 def writeShapeLines(mcMesh, grpName, nodeChains, offsetX=0., offsetY=0.):
     print(" === writeShapeLines", grpName)
     coords = mcMesh.getCoords()
@@ -124,6 +126,7 @@ def writeShapeLines(mcMesh, grpName, nodeChains, offsetX=0., offsetY=0.):
         w.record(chainName)
     w.close()
 
+
 def writeShapePoints(mcMesh, grpName, nodeChains, offsetX=0., offsetY=0.):
     print(" === writeShapePoints", grpName)
     coords = mcMesh.getCoords()
@@ -142,6 +145,7 @@ def writeShapePoints(mcMesh, grpName, nodeChains, offsetX=0., offsetY=0.):
         w.record(chainName)
     w.close()
 
+
 def exploreEdgeGroups(meshFile, offsetX=0., offsetY=0.):
     print(" === exploreEdgeGroups", meshFile)
     mcMesh = ml.MEDFileMesh.New(meshFile)
@@ -157,7 +161,68 @@ def exploreEdgeGroups(meshFile, offsetX=0., offsetY=0.):
         writeShapeLines(mcMesh, grpName, nodeChains, offsetX, offsetY)
         writeShapePoints(mcMesh, grpName, nodeChains, offsetX, offsetY)
 
-# ---
 
-#meshFile = freeBordersGroup('/home/paul/projets/hydro95/V9_5_BR/tests/Maill_SB1610_ALB1501_barr_eff.med')
-#exploreEdgeGroups(meshFile, 1000000., 6000000.)
+def fitShapePointsToMesh(freeBorderShapefile, shapefileToAdjust):
+    print(" === fitShapePointsToMesh", freeBorderShapefile, shapefileToAdjust)
+    fb = shapefile.Reader(freeBorderShapefile)
+    fbShapes = fb.shapes()
+    maxbbox=[1.e30, 1.e30, -1.e30, -1.e30]
+    outerBboxIndex = -1
+    for i,s in enumerate(fbShapes):
+        bb = s.bbox
+        if (bb[0] < maxbbox[0] and bb[1] < maxbbox[1] and bb[2] > maxbbox[2] and bb[3] > maxbbox[3]):
+            maxbbox = bb
+            outerBboxIndex = i
+    fbs = fbShapes[outerBboxIndex]
+
+    sf = shapefile.Reader(shapefileToAdjust)
+    shapes = sf.shapes()
+    sfta = sf.shape(0)
+    pdist =[]
+    for i,p in enumerate(sfta.points):
+        dmin = 1.e30
+        jmin = -1
+        for j,pfb in enumerate(fbs.points):
+            d = math.sqrt((pfb[0]-p[0])*(pfb[0]-p[0]) + (pfb[1]-p[1])*(pfb[1]-p[1]))
+            if d < dmin:
+                dmin = d
+                jmin = j
+        pdist.append((dmin, jmin, i)) # distance, index in freeBorder, index in shapeToAdjust
+    pdist.sort() # the 2 closest points must be set on the mesh freeBorder
+    print(pdist)
+    i1 = min(pdist[0][2], pdist[1][2]) # index of first adjusted point in shapeToAdjust
+    i2 = max(pdist[0][2], pdist[1][2]) # index of second adjusted point in shapeToAdjust
+    if i1 == pdist[0][2]:
+        ifb1 = pdist[0][1]             # index of first adjusted point in freeBorder
+        ifb2 = pdist[1][1]             # index of second adjusted point in freeBorder
+    else:
+        ifb1 = pdist[1][1]
+        ifb2 = pdist[0][1]
+
+    a = os.path.splitext(shapefileToAdjust)
+    shapefileAdjusted = a[0] + '_adj' + a[1]
+    chainName = os.path.basename(a[0])
+
+    w = shapefile.Writer(shapefileAdjusted)
+    w.shapeType = 3
+    w.field('name', 'C')
+
+    chaincoords = []
+    chaincoords.append(fbs.points[ifb1])
+    for i in range(i1+1, i2):
+        chaincoords.append(sfta.points[i])
+    chaincoords.append(fbs.points[ifb2])
+    w.line([chaincoords])
+    w.record(chainName + '_0')
+
+    chaincoords = []
+    chaincoords.append(fbs.points[ifb2])
+    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])
+    w.line([chaincoords])
+    w.record(chainName + '_1')
+    w.close()
+
-- 
2.39.2