docstrings and restructure fix #8

Makefile

 all :
 	python setup.py install --user
-	make -C doc html
 
 html:
-	make -C doc html
+	make -C doc html &> /dev/null

doc/conf.py

@@ -30,7 +30,7 @@ release = '0.1'
 # Add any Sphinx extension module names here, as strings. They can be
 # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
 # ones.
-extensions = ['sphinx.ext.autodoc','sphinx.ext.napoleon']
+extensions = ['sphinx.ext.autodoc','sphinx.ext.napoleon','sphinx_autodoc_typehints']
 
 napoleon_google_docstring = True
 napoleon_numpy_docstring = False

doc/index.rst

@@ -8,7 +8,6 @@ Welcome to msea's documentation!
    :caption: Contents:
 
 
-
 * :ref:`genindex`
 
 ..

doc/msea.Field.rst → doc/msea.field.rst

-msea.Field module
+msea.field module
 =================
 
-.. automodule:: msea.Field
+.. automodule:: msea.field
    :members:
    :undoc-members:
    :show-inheritance:

doc/msea.geometry.rst

+msea.geometry module
+====================
+
+.. automodule:: msea.geometry
+   :members:
+   :undoc-members:
+   :show-inheritance:

doc/msea.gmsh.rst

+msea.gmsh module
+=====================
+
+.. automodule:: msea.gmsh
+   :members:
+   :undoc-members:
+   :show-inheritance:

doc/msea.rst

 msea package
 ============
 
-Submodules
-----------
-
 .. toctree::
    :maxdepth: 4
-
-   msea.Field
-
-Module contents
----------------
-
-.. automodule:: msea
-   :members:
-   :undoc-members:
-   :show-inheritance:
+   
+   msea.field
+   msea.gmsh
+   msea.geometry

msea/__init__.py

-from . import Field
-from ._private import Domain, CurveType, convert_msh_gis
-__all__ = ["Field","Domain","CurveType","convert_msh_gis"]
-
+from . import field
+from . import gmsh
+from . import geometry

msea/Field.py → msea/field.py

-from ._private import Domain
+from .geometry import Domain
 import numpy as np
-from osgeo import osr
+from osgeo import osr, gdal
+from scipy.spatial import cKDTree
+import logging
+from typing import List 
+import itertools
+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
+
 
 class Distance:
     """Callable evaluating the distance to a set of discretized curves.
@@ -9,24 +24,25 @@ class Distance:
     closest point is computed.
     """
 
-    def __init__(self, domain:Domain, sampling:float, tags=None):
+    def __init__(self, domain: Domain, sampling: float,
+            tags: List[str] = None):
         """
         Args:
             domain: a Domain object containing the set of curves
             sampling: the interval between two consecutive sampling points.
             tags: List of physical tags specifying the curve from the domain.
-                if None, all curves are taken into account (default: None)
+                if None, all curves are taken into account.
         """
         points = []
-        for curve in domain._curves:
+        for curve in itertools.chain(domain._curves, domain._interior_curves):
             if (tags is None) or (curve.tag in tags):
-                points.append(curve.sample(sampling))
+                points.append(_curve_sample(curve,sampling))
         points = np.vstack(points)
         self._tree = cKDTree(points)
         self._projection = domain._projection
 
-    def __call__(self, x:np.ndarray, projection:osr.SpatialReference
-                 )->np.ndarray:
+    def __call__(self, x: np.ndarray, projection: osr.SpatialReference
+                 ) -> np.ndarray:
         """Compute the distance between each point of x and the curves.
 
         Args:
@@ -46,10 +62,10 @@ class Distance:
 class Raster:
     """Callable to evaluate a raster field loaded from a file."""
 
-    def __init__(self, filename:str):
+    def __init__(self, filename: str):
         """
         Args:
-            filename: A geotiff file or any other raster file supported by gdal.
+            filename: A geotiff file or any other raster supported by gdal.
         """
         src_ds = gdal.Open(filename)
         self._geo_matrix = src_ds.GetGeoTransform()
@@ -58,8 +74,8 @@ class Raster:
         assert(self._geo_matrix[4] == 0.)
         self._projection = src_ds.GetProjection()
 
-    def __call__(self, x:np.ndarray, projection:osr.SpatialReference
-                 )->np.ndarray:
+    def __call__(self, x: np.ndarray, projection: osr.SpatialReference
+                 ) -> np.ndarray:
         """Evaluate the field value on each point of x.
 
         Keyword arguments:
@@ -73,4 +89,3 @@ class Raster:
         xi = (x[:, 0]-gm[3])/gm[5]
         eta = (x[:, 1]-gm[0])/gm[1]
         return self._data[xi.astype(int), eta.astype(int)]
-

msea/gmsh.py

+import gmsh
+import numpy as np
+import uuid
+from typing import Dict, Any, Union, Tuple
+
+from .geometry import Domain, CurveType, MeshSizeCallback
+
+gmsh.initialize()
+gmsh.option.setNumber("General.Terminal", 1)
+GmshOptionValue = Union[float,str,Tuple[float, float, float]]
+
+def _gmsh_curve_geo(curve_type: CurveType, pointsid):
+    pts = list(i+1 for i in pointsid)
+    if curve_type == CurveType.SPLINE:
+        return [gmsh.model.geo.addSpline(pts)]
+    elif curve_type == CurveType.BSPLINE:
+        return [gmsh.model.geo.addBSpline(pts)]
+    elif curve_type == CurveType.STRICTPOLYLINE:
+        print("STRICT POLYLINE !!!")
+        pairs = zip(pts[:-1], pts[1:])
+        return list(gmsh.model.geo.addLine(p0, p1) for p0, p1 in pairs)
+    elif curve_type == CurveType.POLYLINE:
+        return [gmsh.model.geo.addPolyline(pts)]
+    else:
+        raise ValueError("unknown curve_type '%s'" % curve_type)
+
+def _curve_sample(curve, lc):
+    gmsh.model.add(str(uuid.uuid4()))
+    tags = list([gmsh.model.geo.addPoint(*x, 0) -
+                 1 for x in curve.points[:-1, :]])
+    if np.linalg.norm(curve.points[0, :]-curve.points[-1, :]) < 1e-8:
+        tags.append(tags[0])
+    else:
+        tags.append(gmsh.model.geo.addPoint(*curve.points[-1, :], 0)-1)
+    ltag = _gmsh_curve_geo(curve.curve_type, tags)
+    gmsh.model.geo.synchronize()
+    gmsh.option.setNumber("Mesh.LcIntegrationPrecision", 1e-3)
+    gmsh.option.setNumber("Mesh.CharacteristicLengthMin", lc)
+    gmsh.option.setNumber("Mesh.CharacteristicLengthMax", lc)
+    gmsh.model.mesh.generate(1)
+    nodes = gmsh.model.mesh.getNodes(1, ltag[0], includeBoundary=True)[
+        1].reshape([-1, 3])
+    gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0)
+    gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 1e22)
+    gmsh.model.remove()
+    return nodes[:, :2]
+
+
+def mesh(domain: Domain, filename: str, mesh_size: MeshSizeCallback,
+         version:float = 4.0, gmsh_options: Dict[str,GmshOptionValue] = {}):
+    """ Use gmsh to generate a mesh from a geometry and a mesh size callback
+
+    Args:
+        domain: the input geometry
+        filename: output mesh file (.msh)
+        mesh_size: callbable prescribing the mesh element size
+        version: msh file version (typically 2.0 or 4.0)
+        gmsh_options: additional options passed to gmsh
+    """
+
+    domain._build_topology()
+    nadapt = 3
+    nadapt1d = 4
+    for curve in domain._curves:
+        curve.mesh_size = mesh_size(curve.points, domain._projection)
+    gmsh.option.setNumber("Mesh.CharacteristicLengthFromPoints", 0)
+    gmsh.option.setNumber("Mesh.LcIntegrationPrecision", 1e-3)
+    gmsh.option.setNumber("Mesh.CharacteristicLengthFactor", 1)
+    gmsh.option.setNumber("Mesh.CharacteristicLengthFromCurvature", 0)
+    gmsh.option.setNumber("Mesh.CharacteristicLengthExtendFromBoundary", 0)
+    gmsh.option.setNumber("Mesh.CharacteristicLengthFromParametricPoints",
+                          1)
+    ctags = []
+    cltag = []
+    physicals = {}
+    for i, x in enumerate(domain._points):
+        gmsh.model.geo.addPoint(*x, 0, tag=i+1)
+    for cl in domain._curveloops:
+        ctag = []
+        for i, (l, o) in enumerate(cl):
+            curve = domain._curves[l]
+            tags = _gmsh_curve_geo(curve.curve_type, curve.pointsid)
+            physicals.setdefault(curve.tag, []).extend(tags)
+            ctag.extend(tags)
+        cltag.append(gmsh.model.geo.addCurveLoop(ctag))
+        ctags += ctag
+    stag = gmsh.model.geo.addPlaneSurface(cltag)
+    embeded = []
+    for curve in domain._interior_curves:
+        tags = _gmsh_curve_geo(curve.curve_type, curve.pointsid)
+        embeded.extend(tags)
+        physicals.setdefault(curve.tag, []).extend(tags)
+    gmsh.model.geo.synchronize()
+    gmsh.model.mesh.embed(1, embeded, 2, stag)
+    gmsh.model.mesh.embed(
+        0, [i+1 for i in domain._interior_points_id], 2, stag)
+    it = 0
+    for cl in domain._curveloops:
+        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])
+            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)
+    tag = gmsh.model.addPhysicalGroup(2, [stag])
+    gmsh.model.setPhysicalName(2, stag, "domain")
+    # 1D mesh
+    gmsh.model.mesh.generate(1)
+    for i in range(nadapt1d):
+        for (dim, tag) in gmsh.model.getEntities(1):
+            _, x, u = gmsh.model.mesh.getNodes(dim, tag, True)
+            size = mesh_size(x.reshape([-1, 3]), domain._projection)
+            gmsh.model.mesh.setSizeAtParametricPoints(dim, tag, u, size)
+        print("pass ", i, nadapt)
+        gmsh.model.mesh.clear()
+        gmsh.model.mesh.generate(1)
+    # gmsh.fltk.run()
+    # 2D mesh ##
+    initial_sampling = 0.1
+    gmsh.model.mesh.generate(1)
+    gmsh.option.setNumber("Mesh.CharacteristicLengthMin", initial_sampling)
+    gmsh.option.setNumber("Mesh.CharacteristicLengthMax", initial_sampling)
+    gmsh.option.setNumber("Mesh.CharacteristicLengthExtendFromBoundary", 0)
+    gmsh.option.setNumber("Mesh.CharacteristicLengthFromParametricPoints",
+                          0)
+    gmsh.model.mesh.generate(2)
+    gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0)
+    gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 1e22)
+    bg_field = gmsh.model.mesh.field.add("PostView")
+    gmsh.model.mesh.field.setAsBackgroundMesh(bg_field)
+    for i in range(nadapt):
+        node_tags, node_x, _ = gmsh.model.mesh.getNodes(-1, -1)
+        node_x = node_x.reshape([-1, 3])
+        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)
+        gmsh.model.mesh.field.setNumber(bg_field, "ViewTag", sf_view)
+        gmsh.model.mesh.generate(2)
+        gmsh.view.remove(sf_view)
+    gmsh.model.mesh.field.remove(bg_field)
+    gmsh.fltk.run()
+    gmsh.option.setNumber("Mesh.MshFileVersion", version)
+    gmsh.write(filename)
+    gmsh.finalize()
+
+
+def convet_to_gis(input_filename: str, output_filename: str):
+    """ Convert a triangular gmsh mesh into shapefile or geopackage.
+
+    Args:
+        input_filename : any mesh file readable by gmsh (typically
+            a .msh file)
+        output_filename : shape file (.shp) or geopackage (.gpkg) file
+    """
+    gmsh.model.add(str(uuid.uuid4()))
+    gmsh.open(inputfn)
+    if outputfn.endswith(".shp"):
+        shpdriver = ogr.GetDriverByName('ESRI Shapefile')
+    elif outputfn.endswith(".gpkg"):
+        shpdriver = ogr.GetDriverByName('GPKG')
+    else:
+        raise ValueError("Unknown file extension '" + outputfn+"'")
+    if os.path.exists(outputfn):
+        shpdriver.DeleteDataSource(outputfn)
+    out_data_source = shpdriver.CreateDataSource(outputfn)
+    out_layer = out_data_source.CreateLayer(outputfn, geom_type=ogr.wkbPolygon)
+    out_layer.CreateField(ogr.FieldDefn("id", ogr.OFTInteger))
+    id_field = out_layer.FindFieldIndex("id", 1)
+
+    tri_i, tri_n = gmsh.model.mesh.getElementsByType(2)
+    node_tags, x, _ = gmsh.model.mesh.getNodesByElementType(2)
+    x = x.reshape([-1, 3])
+    nodes_map = dict({tag: i for i, tag in enumerate(node_tags)})
+    tri_n = np.array(list(nodes_map[t] for t in tri_n)).reshape([-1, 3])
+
+    tri_n = np.column_stack([tri_n, tri_n[:, [0]]])
+    header = struct.pack("<biii", 1, 3, 1, 4)
+    out_data_source.StartTransaction()
+    for i, tn in zip(tri_i, tri_n):
+        wkb = header+x[tn, :2].astype('<f8').tobytes()
+        feat = ogr.Feature(out_layer.GetLayerDefn())
+        feat.SetGeometryDirectly(ogr.CreateGeometryFromWkb(wkb))
+        feat.SetField(id_field, int(i))
+        out_layer.CreateFeature(feat)
+        feat = None
+    out_data_source.CommitTransaction()
+    gmsh.model.remove()

testshp.py

-import msea
-import numpy as np
-from osgeo import osr 
-
-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)
-
-domain_srs = osr.SpatialReference()
-
-domain_srs.ImportFromEPSG(32631)
-#domain_srs.ImportFromProj4("+proj=utm +ellps=WGS84 +zone=31")
-
-domain = msea.Domain(domain_srs)
-domain.add_shapefile("test/data_no_duplicate.shp","physical",curve_type=msea.POLYLINE)
-#domain.add_shapefile("test/interior.shp",None,True,curve_type=msea.POLYLINE)
-#domain.add_interior_points("test/interior.shp")
-bath_field = msea.RasterField("medit.tiff")
-dist_coast_field = msea.DistanceField(domain,100,["coast","island"])
-dist_porquerolles_field = msea.DistanceField(domain,20,["porquerolles"])
-domain.unrefine_boundaries(8e5,4.68e6,20,mesh_size)
-msea.mesh_gmsh(domain,"test.msh",mesh_size)