doc

doc/conf.py

@@ -22,7 +22,7 @@ copyright = '2020, Jonathan Lambrechts'
 author = 'Jonathan Lambrechts'
 
 # The full version, including alpha/beta/rc tags
-version = "0.0.1"
+version = "0.1"
 commit_tag = os.environ.get("CI_COMMIT_TAG")
 if commit_tag and (commit_tag.startswith("v-") or commit_tag.startswith("w-")):
     version = commit_tag[2:]

doc/gen_examples.py

@@ -20,7 +20,7 @@ for filename in glob.glob("../examples/*.py") :
                     inrst = True
                     fout.write("\n")
                     continue
-                if (not lstrip and inrst) or i == 0:
+                if (not lstrip and inrst) or i == 0 or (inrst and lstrip[0] != '#'):
                     fout.write("\n.. code-block:: python\n  \n")
                     inrst = False
                     continue

doc/index.rst

@@ -5,7 +5,7 @@ Seamsh
 
 .. toctree::
    :maxdepth: 2
-   :caption: Contents:
+   :caption: Contents
 
    self
    install

examples/1-simple.py

-# seamsh - Copyright (C) <2010-2020>
-# <Universite catholique de Louvain (UCL), Belgium
-# 	
-# List of the contributors to the development of seamsh: see AUTHORS file.
-# Description and complete License: see LICENSE file.
-# 	
-# This program (seamsh) is free software: 
-# you can redistribute it and/or modify it under the terms of the GNU Lesser 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 Lesser General Public License for more details.
-# 
-# You should have received a copy of the GNU Lesser General Public License
-# along with this program (see COPYING file).  If not, 
-# see <http://www.gnu.org/licenses/>.
-
-# %%
-# my super example simple
-# =======================
+# %%
+# Simple Geometry
+# ===============
+# This example illustrates how to insert different types of curves
+# directly from numpy arrays and osr spatial references.
 
 import seamsh
+from seamsh.geometry import CurveType
 import numpy as np
 from osgeo import osr 
 
-def mesh_size(x) :
-    return (0.0025+(1+np.sin(4*np.pi*x[:,1]*np.pi*x[:,0]))*0.01)*2
+# %%
+# First, let's create a domain object and its associated projection.
+# Any projection supported by osgeo.osr can be used.
 
 domain_srs = osr.SpatialReference()
 domain_srs.ImportFromEPSG(4326)
+domain = seamsh.geometry.Domain(domain_srs)
+
+# %%
+# Boundary curves are created from a list of control points, an osr projection,
+# a physical tag and a curve type. If the projection does not match the domain's
+# projection, the points are re-projected conversion is done.
+# A STRICTPOLYLINE curve linearly interpolates control points and forces a mesh
+# on each control point.
+
+domain.add_boundary_curve([[0,0],[-0.2,0.25],[-0.2,0.75],[0,1]],
+        "wall0",domain_srs,curve_type=CurveType.STRICTPOLYLINE)
+
+# %%
+# A POLYLINE is similar to a STRICTPOLYLINE but no mesh point is forced on the
+# control points, the mesher can cut the corners.
+
+domain.add_boundary_curve([[0,1],[0.25,1.2],[0.75,1.2],[1,1]],
+        "wall1",domain_srs,curve_type=CurveType.POLYLINE)
+
+# %%
+# A SPLINE uses a cubic spline interpolation through the control points
+
+domain.add_boundary_curve([[0,0],[0.25,-0.2],[0.75,-0.2],[1,0]],
+        "wall2",domain_srs,curve_type=CurveType.SPLINE)
+
+# %%
+# A BSPLINE is smoother than a SPLINE but the control points are not on the
+# curves.
+
+domain.add_boundary_curve([[1,1],[1.2,0.75],[1.2,0.25],[1,0]],
+        "wall3",domain_srs,curve_type=CurveType.BSPLINE)
+
+# %%
+# If the last point is the same as the first point, a periodic curve is
+# created.
+
+domain.add_boundary_curve([[0.4,0.6],[0.6,0.6],[0.6,0.4],[0.4,0.4],[0.4,0.6]],
+        "island",domain_srs,curve_type=CurveType.BSPLINE)
+
+# %%
+# Interior curves are meshed on both sides, they do not define the domain
+# boundaries but forces edge alignment and mesh points. When an interior curve
+# can touch a boundary or another interior curve, a control point should be 
+# present in both curves at the intersection.
+
+domain.add_interior_curve([[-0.2,0.75],[0.1,0.6],[0.2,0.4],[0.1,0.4]],
+        "interior1",domain_srs,curve_type=CurveType.BSPLINE)
+domain.add_interior_curve([[0.1,0.4],[0.1,0.6],[0.1,0.8]],
+        "interior2",domain_srs,curve_type=CurveType.POLYLINE)
+
+# %%
+# A callback which takes an numpy array of points and a projection as argument
+# defines the mesh element size. In this case, a simple analytical function is
+# used.
+
+def mesh_size(x, projection) :
+    return (0.005+(1+np.sin(4*np.pi*x[:,1]*np.pi*x[:,0]))*0.01)
+
+# %%
+# Finally, the seamesh.gmsh.mesh function generates the mesh.
 
-domain = seamsh.Domain(domain_srs)
-domain.add_curve([[0,0],[-0.2,0.25],[-0.2,0.75],[0,1]],
-        "wall",domain_srs,curve_type=seamsh.SPLINE)
-domain.add_curve([[0,1],[0.25,1.2],[0.75,1.2],[1,1]],
-        "wall1",domain_srs,curve_type=seamsh.BSPLINE)
-domain.add_curve([[0,0],[0.25,-0.2],[0.75,-0.2],[1,0]],
-        "wall2",domain_srs,curve_type=seamsh.STRICTPOLYLINE)
-domain.add_curve([[1,1],[1.2,0.75],[1.2,0.25],[1,0]],
-        "wall3",domain_srs,curve_type=seamsh.STRICTPOLYLINE)
-domain.add_curve([[0.4,0.6],[0.6,0.6],[0.6,0.4],[0.4,0.4],[0.4,0.6]],
-        "wallin",domain_srs,curve_type=seamsh.BSPLINE)
-domain.build_topology()
-
-seamsh.mesh_gmsh(domain,"test.msh",mesh_size)
+seamsh.gmsh.mesh(domain,"basic_geometry.msh",mesh_size,version=2.0)

examples/2-gis.py

+# %%
+# GIS Data
+# ========
+# This example illustrates how to insert different types of curves
+# directly from numpy arrays and osr spatial references.
+# 
+# For the sake of this example we, start by downloading some data.
+
+import os
+import urllib.request
+import tarfile
+if not os.path.isdir("data") :
+    urllib.request.urlretrieve("ftp://braque.mema.ucl.ac.be/seamsh/data-test-1.tar.gz",
+                               "data-test-1.tar.gz")
+    f = tarfile.open("data-test-1.tar.gz","r:*")
+    f.extractall()
+
+# %%
+# The actual example starts here.
+
+import seamsh
+from seamsh.geometry import CurveType
+import seamsh.geometry
+import numpy as np
+from osgeo import osr 
+
+# %%
+# First, let's create a domain object and its associated projection.
+# Any projection supported by osgeo.osr can be used.
+
+domain_srs = osr.SpatialReference()
+domain_srs.ImportFromProj4("+proj=utm +ellps=WGS84 +zone=31")
+domain = seamsh.geometry.Domain(domain_srs)
+
+# %%
+# We load all curves from a given ESTRI shapefile as POLYLINEs.
+# In the shapefile, a field named "physical" defines the physical tag of each
+# curve. If a re-projection is required, it will be done automatically.
+
+domain.add_boundary_curves_shp("data/data_no_duplicate.shp","physical",CurveType.POLYLINE)
+
+# %%
+# Interior curves (and interior points) can be loaded in a similar way.
+# In this case we do not assign any physical tag.
+
+domain.add_interior_curves_shp("data/interior.shp",None,CurveType.STRICTPOLYLINE)
+
+# %%
+# Seamsh provides helper classes to compute the element size field.
+#
+# field.Raster allows to load geotiff files (or any raster file supported by gdal)
+
+bath_field = seamsh.field.Raster("data/medit.tiff")
+
+# %%
+# field.Distance can be used to compute the distance from given (tagged) boundaries,
+# a first field returns the distance from boundaries with physical tag "coast", or 
+# "island". The curves are sampled at regular intervals and the computed distance is 
+# actually the distance from the closest sampling point. The second argument compute
+# the length of the interval between the sampling points.
+
+dist_coast_field = seamsh.field.Distance(domain,100,["coast","island"])
+
+# %%
+# A second distance fields return the distance from the island tagged "porquerolles",
+# in the shp file.
+
+dist_porquerolles_field = seamsh.field.Distance(domain,20,["porquerolles"])
+
+# %%
+# The mesh size is defined based on those fields, smaller elements are prescribed
+# around the Porquerolles island.
+
+def mesh_size(x,projection) :
+    s_coast = np.clip((dist_coast_field(x,projection)-400)*0.5,200,5000)
+    s_porq = np.clip((dist_porquerolles_field(x,projection)-200)*0.5,50,5000)
+    s_dist = np.minimum(s_coast,s_porq)
+    return s_dist
+
+# %%
+# another option would be to take a mesh size proportional to the square root
+# of the (clipped) bathymetry :
+
+# def mesh_size(x,projection) :
+#     bath = -bath_field(x,projection)
+#     s_bath = np.sqrt(np.clip(bath,100,4000))*10
+#     return  s_bath*10
+
+# %%
+# The "intermediate_file_name" option is used to save files containing
+# intermediate meshes and mesh size fields. If this parameter takes the
+# special value "-" an interactive gmsh graphical window will pop up
+# after each meshing step.
+
+seamsh.gmsh.mesh(domain,"gis_mesh.msh",mesh_size,intermediate_file_name="debug")
+
+# %%
+# the gmsh.convert_to_gis function can be used to convert a gmsh .msh file
+# in a shape file or a geo package file.
+
+seamsh.gmsh.convert_to_gis("gis_mesh.msh",domain_srs,"gis_mesh.gpkg")

examples/3-coarsening.py

+# %%
+# Boundary Coarsening
+# ===================
+# This example revisits the previous example (GIS Data) with an extra step to
+# reduce the resolution of the input coastlines.
+#
+# This approach has two main advantages:
+#
+#   - Thanks to the fitting of the input data resolution to the final mesh
+#     element size, the quality of the elements near the boundary is improved.
+#   - The input curves do not have to form topologically valid closed loops.
+#     The only requirement is that the external boundary should be closed 
+#     but the lines defining the boundaries can intersect each others. 
+# 
+# But for now, the coarsening algorithm presents several limitations:
+#   - The physical tags are lost (this will be fixed soon).
+#   - Interior curves hare not handled, they can be present as long as they
+#     do not intersect the domain boundaries (external and islands) but they
+#     won't be coarsened.
+#
+# This example starts as the previous one.
+
+import os
+import urllib.request
+import tarfile
+if not os.path.isdir("data") :
+    urllib.request.urlretrieve("ftp://braque.mema.ucl.ac.be/seamsh/data-test-1.tar.gz",
+                               "data-test-1.tar.gz")
+    f = tarfile.open("data-test-1.tar.gz","r:*")
+    f.extractall()
+
+import seamsh
+from seamsh.geometry import CurveType
+import seamsh.geometry
+import numpy as np
+from osgeo import osr 
+
+domain_srs = osr.SpatialReference()
+domain_srs.ImportFromProj4("+proj=utm +ellps=WGS84 +zone=31")
+domain = seamsh.geometry.Domain(domain_srs)
+domain.add_boundary_curves_shp("data/data_no_duplicate.shp","physical",CurveType.POLYLINE)
+
+dist_coast_field = seamsh.field.Distance(domain,100,["coast","island"])
+dist_porquerolles_field = seamsh.field.Distance(domain,20,["porquerolles"])
+
+def mesh_size(x,projection) :
+    s_coast = np.clip((dist_coast_field(x,projection)-400)*0.5,200,5000)
+    s_porq = np.clip((dist_porquerolles_field(x,projection)-200)*0.5,50,5000)
+    s_dist = np.minimum(s_coast,s_porq)
+    return s_dist
+
+# %%
+# The geometry.coarsen_boundaries function requires in input the coordinates of one (any)
+# point which inside the domain. The last parameters should be at least two times smaller
+# than the smallest value returned by the mesh_size callback, otherwise the resulting
+# geometry will be invalid.
+
+coarse = seamsh.geometry.coarsen_boundaries(domain,(8e5,4.68e6),domain_srs,mesh_size,20)
+
+# %%
+# The resulting coarsend domain can be meshed normally.
+
+seamsh.gmsh.mesh(coarse,"coarse_boundary_mesh.msh",mesh_size)
+

examples/2-shapefiles.py → examples/4-advanced-gmsh.py

 # %%
-# my super example
-# ================
+# Advanced Gmsh / Quad Mesh
+# =========================
+# Seamsh uses gmsh internally. The whole gmsh python api is exposed through the
+# seamsh.gmsh.gmsh module (the repetition is not a typo). This example is a 
+# variation on the previous one using the gmsh api to generate a quadrangular 
+# instead of triangular mesh.
+# 
+# We start as usual.
 
 import os
 import urllib.request
 import tarfile
 if not os.path.isdir("data") :
-    urllib.request.urlretrieve("ftp://braque.mema.ucl.ac.be/seamsh/data-test-1.tar.gz", "data-test-1.tar.gz")
+    urllib.request.urlretrieve("ftp://braque.mema.ucl.ac.be/seamsh/data-test-1.tar.gz",
+                               "data-test-1.tar.gz")
     f = tarfile.open("data-test-1.tar.gz","r:*")
     f.extractall()
 
+
 import seamsh
 from seamsh.geometry import CurveType
 import seamsh.geometry
 import numpy as np
 from osgeo import osr 
 
-# %%
-# with comments in the middle
+domain_srs = osr.SpatialReference()
+domain_srs.ImportFromProj4("+proj=utm +ellps=WGS84 +zone=31")
+domain = seamsh.geometry.Domain(domain_srs)
+domain.add_boundary_curves_shp("data/data_no_duplicate.shp","physical",CurveType.POLYLINE)
+
+dist_coast_field = seamsh.field.Distance(domain,100,["coast","island"])
+dist_porquerolles_field = seamsh.field.Distance(domain,20,["porquerolles"])
 
 def mesh_size(x,projection) :
     s_coast = np.clip((dist_coast_field(x,projection)-400)*0.5,200,5000)
-    s_porquerolles = np.clip((dist_porquerolles_field(x,projection)-200)*0.5,50,5000)
-    return np.minimum(s_coast,s_porquerolles)
+    s_porq = np.clip((dist_porquerolles_field(x,projection)-200)*0.5,50,5000)
+    s_dist = np.minimum(s_coast,s_porq)
+    return s_dist
 
-domain_srs = osr.SpatialReference()
+coarse = seamsh.geometry.coarsen_boundaries(domain,(8e5,4.68e6),domain_srs,mesh_size,20)
 
-domain_srs.ImportFromEPSG(32631)
-#domain_srs.ImportFromProj4("+proj=utm +ellps=WGS84 +zone=31")
+# %%
+# The gmsh options are set. The full list of gmsh mesh options can be found in
+# `gmsh's documentation <https://gmsh.info/doc/texinfo/gmsh.html#Mesh-options-list>`_.
+
+seamsh.gmsh.gmsh.option.setNumber("Mesh.RecombineAll",1);
+seamsh.gmsh.gmsh.option.setNumber("Mesh.RecombinationAlgorithm",1);
+seamsh.gmsh.gmsh.option.setNumber("Mesh.Algorithm",8);
+
+# %%
+# Then the mesh is generated.
+
+seamsh.gmsh.mesh(coarse,"quad_mesh.msh",mesh_size)
 
-domain = seamsh.geometry.Domain(domain_srs)
-domain.add_boundary_curves_shp("data/data_no_duplicate.shp","physical",CurveType.POLYLINE)
-#domain.add_interior_curves_shp("data/interior.shp",None,CurveType.STRICTPOLYLINE)
-#domain.add_interior_points_shp("data/interior.shp")
-bath_field = seamsh.field.Raster("data/medit.tiff")
-dist_coast_field = seamsh.field.Distance(domain,100,["coast","island"])
-dist_porquerolles_field = seamsh.field.Distance(domain,20,["porquerolles"])
-#coarse = seamsh.geometry.coarsen_boundaries(domain,(8e5,4.68e6),domain_srs,mesh_size,20)
-seamsh.gmsh.mesh(domain,"test.msh",mesh_size,intermediate_file_name="log")
-seamsh.gmsh.convert_to_gis("test.msh",domain_srs,"test.gpkg")
-
-seamsh.gmsh.gmsh.model.add("test")
-seamsh.gmsh.gmsh.open("test.msh")
-tag,nodes = seamsh.gmsh.gmsh.model.mesh.getElementsByType(2)
-ntri = len(tag)
-print("ntri",ntri)
-assert(ntri>40000 and ntri<41000)

seamsh/field.py

@@ -31,7 +31,6 @@ from .gmsh import _curve_sample
 def _ensure_valid_points(x, pfrom, pto):
     x = np.asarray(x)[:, :2]
     if not pfrom.IsSame(pto):
-        logging.info("reprojecting coordinates")
         x = osr.CoordinateTransformation(pfrom, pto).TransformPoints(x)
         x = np.asarray(x)
     return x
@@ -92,7 +91,8 @@ class Raster:
         self._data = src_ds.GetRasterBand(1).ReadAsArray()
         assert(self._geo_matrix[2] == 0.)
         assert(self._geo_matrix[4] == 0.)
-        self._projection = src_ds.GetProjection()
+        self._projection = osr.SpatialReference()
+        self._projection.ImportFromWkt(src_ds.GetProjection())
 
     def __call__(self, x: np.ndarray, projection: osr.SpatialReference
                  ) -> np.ndarray:
@@ -104,7 +104,7 @@ class Raster:
         Returns:
             The field value on points x. [n]
         """
-        x = _ensure_valid_points(x, projection, self._proection)
+        x = _ensure_valid_points(x, projection, self._projection)
         gm = self._geo_matrix
         xi = (x[:, 0]-gm[3])/gm[5]
         eta = (x[:, 1]-gm[0])/gm[1]

seamsh/geometry.py

@@ -44,7 +44,6 @@ lib = c.CDLL(libpath)
 def _ensure_valid_points(x, pfrom, pto):
     x = np.asarray(x)[:, :2]
     if not pfrom.IsSame(pto):
-        logging.warning("reprojecting coordinates !!!!!!!!!!!!! ")
         x = osr.CoordinateTransformation(pfrom, pto).TransformPoints(x)
         x = np.asarray(x)[:, :2]
     return x

seamsh/gmsh.py

@@ -81,7 +81,9 @@ def mesh(domain: Domain, filename: str, mesh_size: MeshSizeCallback,
         version: msh file version (typically 2.0 or 4.0)
         intermediate_file_name: if not None, save intermediate meshes to those
             files for debugging purpose (suffixes and extensions will be
-            appended).
+            appended), if == "-", an interactive gmsh graphical window will pop up
+            after each meshing step.
+
     """
 
     domain._build_topology()
@@ -107,7 +109,7 @@ def mesh(domain: Domain, filename: str, mesh_size: MeshSizeCallback,
             curve = domain._curves[l]
             tags = _gmsh_curve_geo(curve.curve_type, curve.pointsid)
             physicals.setdefault(curve.tag, []).extend(tags)
-            ctag.extend(tags)
+            ctag.extend(tags if o else list(-t for t in reversed(tags)))
         cltag.append(gmsh.model.geo.addCurveLoop(ctag))
         ctags += ctag
     stag = gmsh.model.geo.addPlaneSurface(cltag)
@@ -125,18 +127,11 @@ def mesh(domain: Domain, filename: str, mesh_size: MeshSizeCallback,
         for i, (l, o) in enumerate(cl):
             curve = domain._curves[l]
             sizes = curve.mesh_size if o else np.flip(curve.mesh_size)
-            u0, u1 = gmsh.model.getParametrizationBounds(1, ctags[it])
+            u0, u1 = gmsh.model.getParametrizationBounds(1, abs(ctags[it]))
             u = np.linspace(u0, u1, curve.points.shape[0])
             gmsh.model.mesh.setSizeAtParametricPoints(
                 1, ctags[it], u, sizes)
             it += 1
-            """
-            elif curve.curve_type == STRICTPOLYLINE:
-                for j in range(curve.points.shape[0]-1):
-                    gmsh.model.mesh.setSizeAtParametricPoints(
-                        1, ctags[it], [0, 1], [sizes[j], sizes[j+1]])
-                    it += 1
-            """
     for name, tags in physicals.items():
         tag = gmsh.model.addPhysicalGroup(1, tags)
         gmsh.model.setPhysicalName(1, tag, name)
@@ -154,7 +149,10 @@ def mesh(domain: Domain, filename: str, mesh_size: MeshSizeCallback,
         gmsh.model.mesh.clear()
     gmsh.model.mesh.generate(1)
     if intermediate_file_name is not None :
-        gmsh.write(intermediate_file_name+"_1d.msh")
+        if intermediate_file_name == "-" :
+            gmsh.fltk.run()
+        else :
+            gmsh.write(intermediate_file_name+"_1d.msh")
     # 2D mesh ##
     _, x, u = gmsh.model.mesh.getNodes()
     x = x.reshape([-1,3])
@@ -174,14 +172,20 @@ def mesh(domain: Domain, filename: str, mesh_size: MeshSizeCallback,
         nodes_map = dict({tag: i for i, tag in enumerate(node_tags)})
         sf_view = gmsh.view.add("mesh size field")
         node_lc = mesh_size(node_x, domain._projection)
-        tri = gmsh.model.mesh.getElements(2, stag)[2][0]
-        tri = np.array(list(nodes_map[t] for t in tri)).reshape([-1, 3])
-        xnode = node_x[tri, :].swapaxes(1, 2).reshape([-1, 9])
-        data = np.column_stack([xnode, node_lc[tri]]).reshape([-1])
-        gmsh.view.addListData(sf_view, "ST", tri.shape[0], data)
+        etypes,etags,enodes = gmsh.model.mesh.getElements(2, stag)
+        for etype,enode in zip(etypes,enodes):
+            nn = 3 if etype == 2 else 4
+            enode = list(nodes_map[t] for t in enode)
+            enode = np.array(enode).reshape([-1, nn])
+            xnode = node_x[enode, :].swapaxes(1, 2).reshape([-1, nn*3])
+            data = np.column_stack([xnode, node_lc[enode]]).reshape([-1])
+            gmsh.view.addListData(sf_view, "ST"if etype == 2 else "SQ", enode.shape[0], data)
         gmsh.model.mesh.field.setNumber(bg_field, "ViewTag", sf_view)
         if intermediate_file_name is not None :
-            gmsh.view.write(sf_view,intermediate_file_name+"_2d_"+str(i)+".pos")
+            if intermediate_file_name == "-" :
+                gmsh.fltk.run()
+            else :
+                gmsh.view.write(sf_view,intermediate_file_name+"_2d_"+str(i)+".pos")
         gmsh.model.mesh.generate(2)
         gmsh.view.remove(sf_view)
     gmsh.model.mesh.field.remove(bg_field)

setup.py

@@ -25,7 +25,7 @@ import os
 with open("Readme.rst", "r") as fh:
     long_description = fh.read()
 
-version = "0.0.2"
+version = "0.1"
 commit_tag = os.environ.get("CI_COMMIT_TAG")
 if commit_tag and (commit_tag.startswith("v-") or commit_tag.startswith("w-")):
     version = commit_tag[2:]