Geometry: Add meshing action with mage internal meshing functions.

src/Meshing/AMeshingTool.py

@@ -0,0 +1,27 @@
+# AMeshingTool.py -- Pamhyr
+# Copyright (C) 2023  INRAE
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <https://www.gnu.org/licenses/>.
+
+# -*- coding: utf-8 -*-
+
+from Model.Except import NotImplementedMethodeError
+
+
+class AMeshingTool(object):
+    def __init__(self):
+        super(AMeshingTool, self).__init__()
+
+    def meshing(self, reach, **kwargs):
+        raise NotImplementedMethodeError(self, self.meshing)

src/Meshing/Mage.py

@@ -0,0 +1,232 @@
+# Mage.py -- Pamhyr
+# Copyright (C) 2023  INRAE
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <https://www.gnu.org/licenses/>.
+
+# -*- coding: utf-8 -*-
+
+import os
+import logging
+import tempfile
+
+from ctypes import (
+    cdll,
+    byref, Structure,
+    c_char_p, c_wchar_p,
+    create_string_buffer,
+    POINTER, c_void_p,
+    c_int, c_double, c_bool
+)
+
+from Meshing.AMeshingTool import AMeshingTool
+
+logger = logging.getLogger()
+
+
+class MeshingWithMage(AMeshingTool):
+    def __init__(self):
+        super(MeshingWithMage, self).__init__()
+
+        self._init_bief_lib()
+
+    def _init_bief_lib(self):
+        self._bief_lib = cdll.LoadLibrary(self._lib_path())
+
+        self._init_c_init_bief_from_geo_file()
+        self._init_c_get_nb_sections()
+        self._init_c_get_nb_points_section()
+        self._init_c_set_bief_name()
+        self._init_c_st_to_m_compl()
+        self._init_c_interpolate_profils_pas_transversal()
+        self._init_c_output_bief()
+
+    def _lib_path(self):
+        return os.path.abspath(
+            os.path.join(
+                os.path.dirname(__file__),
+                "..", "..", "mage", "libbief.so"
+            )
+        )
+
+    def _init_c_init_bief_from_geo_file(self):
+        self._c_init_bief_from_geo_file = getattr(
+            self._bief_lib, 'c_init_bief_from_geo_file'
+        )
+        self._c_init_bief_from_geo_file.argtypes = [
+            c_char_p, POINTER(c_int), POINTER(c_int)
+        ]
+        self._c_init_bief_from_geo_file.restype = None
+
+    def _init_c_get_nb_sections(self):
+        self._c_get_nb_sections = getattr(self._bief_lib, 'c_get_nb_sections')
+        self._c_get_nb_sections.argtypes = [POINTER(c_int)]
+        self._c_get_nb_sections.restype = None
+
+    def _init_c_get_nb_points_section(self):
+        self._c_get_nb_points_section = getattr(
+            self._bief_lib, 'c_get_nb_points_section'
+        )
+        self._c_get_nb_points_section.argtypes = [
+            POINTER(c_int), POINTER(c_int)
+        ]
+        self._c_get_nb_points_section.restype = None
+
+    def _init_c_set_bief_name(self):
+        self._c_set_bief_name = getattr(self._bief_lib, 'c_set_bief_name')
+        self._c_set_bief_name.argtypes = [c_char_p]
+        self._c_set_bief_name.restype = None
+
+    def _init_c_st_to_m_compl(self):
+        self._c_st_to_m_compl = getattr(self._bief_lib, 'c_st_to_m_compl')
+        self._c_st_to_m_compl.argtypes = [POINTER(c_int), c_char_p, c_char_p]
+        self._c_st_to_m_compl.restype = None
+
+    def _init_c_interpolate_profils_pas_transversal(self):
+        self._c_interpolate_profils_pas_transversal = getattr(
+            self._bief_lib, 'c_interpolate_profils_pas_transversal'
+        )
+        self._c_interpolate_profils_pas_transversal.argtypes = [
+            POINTER(c_int), POINTER(c_int),
+            c_char_p, c_char_p,
+            POINTER(c_double), POINTER(c_bool),
+            POINTER(c_int), POINTER(c_bool)
+        ]
+        self._c_interpolate_profils_pas_transversal.restype = None
+
+    def _init_c_output_bief(self):
+        self._c_output_bief = getattr(self._bief_lib, 'c_output_bief')
+        self._c_output_bief.argtypes = [c_char_p]
+        self._c_output_bief.restype = None
+
+    #####################
+    # Binding functions #
+    #####################
+
+    def init_bief_from_geo_file(self, name, with_charriage, with_water):
+        cname = create_string_buffer(name.encode())
+        self._c_init_bief_from_geo_file(
+            cname,
+            byref(c_int(with_charriage)),
+            byref(c_int(with_water))
+        )
+
+    def get_nb_sections(self):
+        nb_sections = c_int(0)
+        self._c_get_nb_sections(byref(nb_sections))
+        return nb_sections.value
+
+    def get_nb_points_section(self, section):
+        nb_points = c_int(0)
+        self._c_get_nb_points_section(
+            byref(c_int(section)), byref(nb_points)
+        )
+        return nb_points.value
+
+    def set_bief_name(self, name):
+        cname = create_string_buffer(name.encode())
+        self._c_set_bief_name(cname)
+
+    def st_to_m_compl(self, npoints, tag1='   ', tag2='   '):
+        cnpoints = c_int(npoints)
+        ctag1 = create_string_buffer(tag1.encode())
+        ctag2 = create_string_buffer(tag2.encode())
+
+        self._c_st_to_m_compl(byref(cnpoints), ctag1, ctag2)
+
+    def interpolate_profils_pas_transversal(self, limite1, limite2,
+                                            directrice1, directrice2,
+                                            pas, lplan=False,
+                                            lm=3, lineaire=False):
+        climite1 = c_int(limite1)
+        climite2 = c_int(limite2)
+        cpas = c_double(pas)
+        clplan = c_bool(lplan)
+        clm = c_int(lm)
+        clineaire = c_bool(lineaire)
+        cdirectrice1 = create_string_buffer(directrice1.encode())
+        cdirectrice2 = create_string_buffer(directrice2.encode())
+
+        self._c_interpolate_profils_pas_transversal(
+            byref(climite1), byref(climite2),
+            cdirectrice1, cdirectrice2,
+            byref(cpas), byref(clplan),
+            byref(clm), byref(clineaire)
+        )
+
+    def output_bief(self, name):
+        cname = create_string_buffer(name.encode())
+        self._c_output_bief(cname)
+
+    ###########
+    # Meshing #
+    ###########
+
+    def meshing(self, reach, step: float = 50):
+        with tempfile.TemporaryDirectory() as tmp:
+            st_file = self.export_reach_to_st(reach, tmp)
+            m_file = st_file.rsplit(".ST", 1)[0] + ".M"
+
+            self.load_st_file(st_file)
+
+            ns, npts_max = self.get_reach_stat()
+
+            self.complete_cross_section()
+            self.interpolate_cross_section(ns, step)
+
+            self.export_to_m(m_file)
+
+            self.import_m_file(reach, m_file)
+            return reach
+
+    def export_reach_to_st(self, reach, tmp):
+        tmp_st = os.path.join(tmp, "meshing.ST")
+
+        logger.debug(f"meshing: Export ST to {tmp_st}")
+
+        reach.export_reach(tmp_st)
+        return tmp_st
+
+    def load_st_file(self, st):
+        self.init_bief_from_geo_file(st, 0, 0)
+        self.set_bief_name("tmp")
+
+    def get_reach_stat(self):
+        ns = self.get_nb_sections()
+        npts_max = max(
+            map(
+                lambda i: self.get_nb_points_section(i),
+                range(1, ns)
+            )
+        )
+
+        return ns, npts_max
+
+    def complete_cross_section(self):
+        self.st_to_m_compl(0, ' ', ' ')
+
+    def interpolate_cross_section(self, ns, step: float):
+        self.interpolate_profils_pas_transversal(
+            1, ns,
+            'un', 'np',
+            step
+        )
+
+    def export_to_m(self, m):
+        self.output_bief(m)
+
+    def import_m_file(self, reach, m):
+        reach.purge()
+
+        logger.debug(f"meshing: Import geometry from {m}")
+        reach.import_geometry(m)

src/Model/Geometry/Reach.py

@@ -247,6 +247,11 @@ class Reach(SQLSubModel):
         self._update_profile_numbers()
         self._status.modified()
 
+    def purge(self):
+        self._profiles = []
+        self._update_profile_numbers()
+        self._status.modified()
+
     def move_up_profile(self, index: int):
         if index < len(self.profiles):
             next = index - 1
@@ -525,6 +530,8 @@ class Reach(SQLSubModel):
         list_header = []
         stop_code = "999.999"
 
+        logger.debug(f"Import geometry from {filename}")
+
         with open(filename, encoding="utf-8") as file_st:
             for line in file_st:
                 if not (line.startswith("#") or line.startswith("*")):
@@ -563,21 +570,31 @@ class Reach(SQLSubModel):
 
     # TODO: Move this function to model reach
     def export_reach(self, filename):
-        with open(f"{filename}", "w") as file_st:
+        with open(f"{filename}", "w") as st:
-            for index in range(len(self.profiles)):
+            cnt = 0
-                profile = self.profiles[index]
+            for profile in self.profiles:
+                num = f"{cnt:>6}"
+                c1 = f"{profile.code1:>6}"
+                c2 = f"{profile.code2:>6}"
+                t = f"{len(profile.points):>6}"
+                kp = f"{profile.kp:>12f}"[0:12]
+                pname = profile.name
+                if profile.name == "":
+                    pname = f"p{profile.id:>3}".replace(" ", "0")
 
-                for v in profile.header:
+                st.write(f"{num}{c1}{c2}{t} {kp} {pname}\n")
-                    file_st.write(f"{v:>6}")
-                file_st.write("\n")
 
-                for point in self._data.profile[index_pro].points:
+                for point in profile.points:
-                    for i in [point.x, point.y, point.z, point.name]:
+                    x = f"{point.x:<12.4f}"[0:12]
-                        file_st.write(f"{i:>13.4f}")
+                    y = f"{point.y:<12.4f}"[0:12]
-                        file_st.write("\n")
+                    z = f"{point.z:<12.4f}"[0:12]
+                    n = f"{point.name:<3}"
 
-                file_st.write("    999.9990     999.9990     999.9990")
+                    st.write(f"{x} {y} {z} {n}\n")
-                file_st.write("\n")
+
+                st.write("    999.9990     999.9990     999.9990")
+                st.write("\n")
+                cnt += 1
 
     def get_incline(self):
         first = self.profile(0)

src/Scripts/MageMesh.py

@@ -139,6 +139,7 @@ def init_c_output_bief(bief_lib):
 # Binding functions #
 #####################
 
+
 def init_bief_from_geo_file(name, with_charriage, with_water):
     logger.info("! call init_bief_from_geo_file:")
     cname = create_string_buffer(name.encode())
@@ -148,6 +149,7 @@ def init_bief_from_geo_file(name, with_charriage, with_water):
         byref(c_int(with_water))
     )
 
+
 def get_nb_sections():
     nb_sections = c_int(0)
     c_get_nb_sections(byref(nb_sections))

src/View/Geometry/Window.py

@@ -43,6 +43,8 @@ from View.Tools.PamhyrWindow import PamhyrWindow
 from View.Tools.Plot.PamhyrToolbar import PamhyrPlotToolbar
 from View.Tools.Plot.PamhyrCanvas import MplCanvas
 
+from Meshing.Mage import MeshingWithMage
+
 from View.Geometry.Table import GeometryReachTableModel
 from View.Geometry.PlotXY import PlotXY
 from View.Geometry.PlotAC import PlotAC
@@ -168,6 +170,7 @@ class GeometryWindow(PamhyrWindow):
             "action_add": self.add,
             "action_delete": self.delete,
             "action_edit": self.edit_profile,
+            "action_meshing": self.edit_meshing,
         }
 
         for action in actions:
@@ -243,6 +246,20 @@ class GeometryWindow(PamhyrWindow):
 
         self.tableView.model().blockSignals(False)
 
+    def edit_meshing(self):
+        self.tableView.model().blockSignals(True)
+
+        mesher = MeshingWithMage()
+        mesher.meshing(self._reach)
+
+        self.tableView.model().blockSignals(False)
+
+        self.update_profile_windows()
+        self.plot_xy()
+        self.plot_kpc()
+        self.plot_ac()
+
+
     pyqtSlot(bool)
 
     def changed_profile_slot(self, status):

src/View/ui/GeometryReach.ui

@@ -126,6 +126,7 @@
    <addaction name="action_up"/>
    <addaction name="action_down"/>
    <addaction name="action_export"/>
+   <addaction name="action_meshing"/>
   </widget>
   <action name="action_import">
    <property name="text">
@@ -227,6 +228,11 @@
     <string>Move down selected cross-section(s)</string>
    </property>
   </action>
+  <action name="action_meshing">
+   <property name="text">
+    <string>meshing</string>
+   </property>
+  </action>
  </widget>
  <resources/>
  <connections/>