mesh curves in their own projection

seamsh/field.py

@@ -50,7 +50,7 @@ class Distance:
 
         for icurve, curve in enumerate(all_curves_iter):
             if (tags is None) or (curve.tag in tags):
-                points.append(_curve_sample(curve, size, None))
+                points.append(_curve_sample(curve, size, domain._projection))
                 progress.log("{} features sampled".format(icurve+1))
         for point in _tools.chain(domain._interior_points):
             if (tags is None) or (point.tag in tags):

seamsh/geometry.py

@@ -68,11 +68,12 @@ def _generate_unique_points(x):
 
 class _Curve:
 
-    def __init__(self, points, tag, curve_type):
-        self.points = points[:, :2]
+    def __init__(self, points, tag, curve_type, projection):
+        self.points = _tools.np.array(points)[:, :2]
         self.mesh_size = None
         self.tag = tag
         self.curve_type = curve_type
+        self.projection = projection
 
 
 class _Point:
@@ -133,7 +134,7 @@ class Domain:
                     nc = []
                     for i, j in zip(breaks[:-1], breaks[1:]):
                         ncurve = _Curve(
-                            curve.points[i:j+1], curve.tag, curve.curve_type)
+                            curve.points[i:j+1], curve.tag, curve.curve_type, curve.projection)
                         ncurve.pointsid = curve.pointsid[i:j+1]
                         newcurves.append(ncurve)
                         nc.append(len(newcurves)-1)
@@ -266,8 +267,7 @@ class Domain:
             projection: the points coordinate system
             curve_type: curve interpolation
         """
-        points = _tools.project_points(points, projection, self._projection)
-        curve = _Curve(points, physical_tag, curve_type)
+        curve = _Curve(points, physical_tag, curve_type, projection)
         self._curves.append(curve)
 
     def add_interior_curve(self, points: _tools.np.ndarray, physical_tag: str,
@@ -282,8 +282,7 @@ class Domain:
             projection: the points coordinate system
             curve_type: curve interpolation
         """
-        points = _tools.project_points(points, projection, self._projection)
-        curve = _Curve(points, physical_tag, curve_type)
+        curve = _Curve(points, physical_tag, curve_type, projection)
         self._interior_curves.append(curve)
 
     def add_interior_points_shp(self, filename: str,

seamsh/gmsh.py

@@ -65,20 +65,37 @@ def _curve_sample_gmsh_tag(tag, lc, projection):
         size = _tools.np.hstack([size, esize])[s]
     return x[:, :2], xi, size
 
+class _lcproj():
+    def __init__(self, lc, proj):
+        self.lc = lc
+        self.projection = proj
+    def __call__(self, x, projection):
+        xlc = _tools.project_points(x, self.projection, projection)
+        x2 = _tools.np.copy(x)
+        eps = 1e-8
+        x2[:,0] += eps
+        x2[:,1] += eps
+        xlc = _tools.project_points(x, self.projection, projection)
+        xlc2 = _tools.project_points(x2, self.projection, projection)
+        slc = self.lc(xlc, projection)
+        h = _tools.np.hypot(xlc2[:,0]-xlc[:,0], xlc2[:,1]-xlc[:,1])/_tools.np.hypot(eps,eps)
+        return slc/h
 
 def _curve_sample(curve, lc, projection):
     gmsh.model.add(str(_tools.uuid.uuid4()))
-    tags = list([gmsh.model.geo.addPoint(*x, 0) -
-                 1 for x in curve.points[:-1, :]])
+    #points = _tools.project_points(curve.points, curve.projection, projection)
+    tags = list([gmsh.model.geo.addPoint(*x, 0) - 1
+                 for x in curve.points[:-1, :]])
     if _tools.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()
-    r = _curve_sample_gmsh_tag(ltag[0], lc, projection)[0]
+    plc = _lcproj(lc, curve.projection)
+    r = _curve_sample_gmsh_tag(ltag[0], plc, projection)[0]
     gmsh.model.remove()
-    return r
+    return _tools.project_points(r, curve.projection, projection)
 
 
 def _create_gmsh_geometry(domain: _geometry.Domain):