Repository: meshpro/dmsh (language: python)

Issue:
potentially better force function
* I've been experimenting with this force function from the paper attached Eq. 3.9 
* It indeed converges faster and produces higher mesh qualities than the default DistMesh one. <img width="770" alt="Screen Shot 2021-09-05 at 12 36 56 PM" src="https://user-images.githubusercontent.com/18619644/132134456-d4ca40d6-c469-4de9-a2b5-f639e9a84ed4.png">. 
* Perhaps it would be useful to implement here?

[imp_distmesh2d.pdf](https://github.com/nschloe/dmsh/files/7112006/imp_distmesh2d.pdf)





Discussion / hints:
Interesting! I'll check it out later.

Fix patch:
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 7e71406..f40c088 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -3,27 +3,21 @@ name: ci
 on:
   push:
     branches:
-    - main
+      - main
   pull_request:
     branches:
-    - main
+      - main
 
 jobs:
   lint:
     runs-on: ubuntu-latest
     steps:
-    - uses: actions/setup-python@v2
-      with:
-        python-version: "3.x"
-    - uses: actions/checkout@v2
-    - name: Lint with flake8
-      run: |
-        pip install flake8
-        flake8 .
-    - name: Lint with black
-      run: |
-        pip install black
-        black --check .
+      - name: Check out repo
+        uses: actions/checkout@v2
+      - name: Set up Python
+        uses: actions/setup-python@v2
+      - name: Run pre-commit
+        uses: pre-commit/action@v2.0.3
 
   build:
     runs-on: ubuntu-latest
@@ -31,13 +25,13 @@ jobs:
       matrix:
         python-version: [3.7, 3.8, 3.9]
     steps:
-    - uses: actions/setup-python@v2
-      with:
-        python-version: ${{ matrix.python-version }}
-    - uses: actions/checkout@v2
-    - name: Test with tox
-      run: |
-        pip install tox
-        tox -- --cov dmsh --cov-report xml --cov-report term
-    - uses: codecov/codecov-action@v1
-      if: ${{ matrix.python-version == '3.9' }}
+      - uses: actions/setup-python@v2
+        with:
+          python-version: ${{ matrix.python-version }}
+      - uses: actions/checkout@v2
+      - name: Test with tox
+        run: |
+          pip install tox
+          tox -- --cov dmsh --cov-report xml --cov-report term
+      - uses: codecov/codecov-action@v1
+        if: ${{ matrix.python-version == '3.9' }}
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
new file mode 100644
index 0000000..57b6939
--- /dev/null
+++ b/.pre-commit-config.yaml
@@ -0,0 +1,16 @@
+repos:
+  - repo: https://github.com/PyCQA/isort
+    rev: 5.9.3
+    hooks:
+      - id: isort
+
+  - repo: https://github.com/psf/black
+    rev: 21.8b0
+    hooks:
+      - id: black
+        language_version: python3
+
+  - repo: https://github.com/PyCQA/flake8
+    rev: 3.9.2
+    hooks:
+      - id: flake8
diff --git a/README.md b/README.md
index 47ec875..eb524b9 100644
--- a/README.md
+++ b/README.md
@@ -22,9 +22,9 @@ requires a lot of memory, and isn't terribly robust either.
 
 On the plus side,
 
- * it's got a user-friendly interface,
- * is pure Python (and hence easily installable on any system), and
- * it produces pretty high-quality meshes.
+- it's got a user-friendly interface,
+- is pure Python (and hence easily installable on any system), and
+- it produces pretty high-quality meshes.
 
 Combined with [optimesh](https://github.com/nschloe/optimesh), dmsh produces the
 highest-quality 2D meshes in the west.
@@ -33,8 +33,8 @@ highest-quality 2D meshes in the west.
 
 #### Primitives
 
-<img alt="circle" src="https://nschloe.github.io/dmsh/circle.svg" width="100%"> | <img alt="circle" src="https://nschloe.github.io/dmsh/rectangle.svg" width="100%"> | <img alt="circle" src="https://nschloe.github.io/dmsh/polygon.svg" width="100%">
-|:---:|:---:|:---:|
+| <img alt="circle" src="https://nschloe.github.io/dmsh/circle.svg" width="100%"> | <img alt="circle" src="https://nschloe.github.io/dmsh/rectangle.svg" width="100%"> | <img alt="circle" src="https://nschloe.github.io/dmsh/polygon.svg" width="100%"> |
+| :-----------------------------------------------------------------------------: | :--------------------------------------------------------------------------------: | :------------------------------------------------------------------------------: |
 
 ```python
 import dmsh
@@ -82,8 +82,8 @@ X, cells = dmsh.generate(geo, 0.1)
 
 ##### Difference
 
-<img src="https://nschloe.github.io/dmsh/moon.svg" width="100%"> | <img src="https://nschloe.github.io/dmsh/pacman.svg" width="100%"> | <img src="https://nschloe.github.io/dmsh/rectangle-hole-refinement.svg" width="100%">
-:-------------------:|:------------------:|:----:|
+| <img src="https://nschloe.github.io/dmsh/moon.svg" width="100%"> | <img src="https://nschloe.github.io/dmsh/pacman.svg" width="100%"> | <img src="https://nschloe.github.io/dmsh/rectangle-hole-refinement.svg" width="100%"> |
+| :--------------------------------------------------------------: | :----------------------------------------------------------------: | :-----------------------------------------------------------------------------------: |
 
 ```python
 import dmsh
@@ -91,6 +91,7 @@ import dmsh
 geo = dmsh.Circle([-0.5, 0.0], 1.0) - dmsh.Circle([+0.5, 0.0], 1.0)
 X, cells = dmsh.generate(geo, 0.1)
 ```
+
 ```python
 import dmsh
 
@@ -100,6 +101,7 @@ X, cells = dmsh.generate(geo, 0.1, tol=1.0e-10)
 
 The following example uses a nonconstant edge length; it depends on the distance to the
 circle `c`.
+
 ```python
 import dmsh
 import numpy as np
@@ -113,8 +115,8 @@ X, cells = dmsh.generate(geo, lambda pts: np.abs(c.dist(pts)) / 5 + 0.05, tol=1.
 
 ##### Union
 
-<img src="https://nschloe.github.io/dmsh/union-circles.svg" width="100%"> | <img src="https://nschloe.github.io/dmsh/union-rectangles.svg" width="100%"> | <img src="https://nschloe.github.io/dmsh/union-three-circles.svg" width="100%">
-:-------------------:|:------------------:|:----:|
+| <img src="https://nschloe.github.io/dmsh/union-circles.svg" width="100%"> | <img src="https://nschloe.github.io/dmsh/union-rectangles.svg" width="100%"> | <img src="https://nschloe.github.io/dmsh/union-three-circles.svg" width="100%"> |
+| :-----------------------------------------------------------------------: | :--------------------------------------------------------------------------: | :-----------------------------------------------------------------------------: |
 
 ```python
 import dmsh
@@ -122,12 +124,14 @@ import dmsh
 geo = dmsh.Circle([-0.5, 0.0], 1.0) + dmsh.Circle([+0.5, 0.0], 1.0)
 X, cells = dmsh.generate(geo, 0.15)
 ```
+
 ```python
 import dmsh
 
 geo = dmsh.Rectangle(-1.0, +0.5, -1.0, +0.5) + dmsh.Rectangle(-0.5, +1.0, -0.5, +1.0)
 X, cells = dmsh.generate(geo, 0.15)
 ```
+
 ```python
 import dmsh
 import numpy as np
@@ -145,8 +149,8 @@ X, cells = dmsh.generate(geo, 0.15)
 
 #### Intersection
 
-<img src="https://nschloe.github.io/dmsh/intersection-circles.svg" width="100%"> | <img src="https://nschloe.github.io/dmsh/intersection-three-circles.svg" width="100%"> | <img src="https://nschloe.github.io/dmsh/intersection-circle-halfspace.svg" width="100%">
-:-------------------:|:------------------:|:----:|
+| <img src="https://nschloe.github.io/dmsh/intersection-circles.svg" width="100%"> | <img src="https://nschloe.github.io/dmsh/intersection-three-circles.svg" width="100%"> | <img src="https://nschloe.github.io/dmsh/intersection-circle-halfspace.svg" width="100%"> |
+| :------------------------------------------------------------------------------: | :------------------------------------------------------------------------------------: | :---------------------------------------------------------------------------------------: |
 
 ```python
 import dmsh
@@ -171,6 +175,7 @@ X, cells = dmsh.generate(geo, 0.1, tol=1.0e-10)
 ```
 
 The following uses the `HalfSpace` primtive for cutting off a circle.
+
 ```python
 import dmsh
 
@@ -180,8 +185,8 @@ X, cells = dmsh.generate(geo, 0.1)
 
 ### Rotation, translation, scaling
 
-<img src="https://nschloe.github.io/dmsh/rotation.svg" width="100%"> | <img src="https://nschloe.github.io/dmsh/scaling.svg" width="100%">
-|:----:|:----:|
+| <img src="https://nschloe.github.io/dmsh/rotation.svg" width="100%"> | <img src="https://nschloe.github.io/dmsh/scaling.svg" width="100%"> |
+| :------------------------------------------------------------------: | :-----------------------------------------------------------------: |
 
 ```python
 import dmsh
@@ -190,12 +195,14 @@ import numpy as np
 geo = dmsh.Rotation(dmsh.Rectangle(-1.0, +2.0, -1.0, +1.0), 0.1 * np.pi)
 X, cells = dmsh.generate(geo, 0.1, tol=1.0e-10)
 ```
+
 ```python
 import dmsh
 
 geo = dmsh.Rectangle(-1.0, +2.0, -1.0, +1.0) + [1.0, 1.0]
 X, cells = dmsh.generate(geo, 0.1)
 ```
+
 ```python
 import dmsh
 
@@ -209,6 +216,7 @@ X, cells = dmsh.generate(geo, 0.1, tol=1.0e-5)
 
 All objects can be used to refine the mesh according to the distance to the object;
 e.g. a `Path`:
+
 ```python
 import dmsh
 
@@ -224,8 +232,8 @@ def edge_size(x):
 X, cells = dmsh.generate(geo, edge_size, tol=1.0e-10)
 ```
 
-
 ### Custom shapes
+
 It is also possible to define your own geometry. Simply create a class derived from
 `dmsh.Geometry` that contains a `dist` method and a method to project points onto the
 boundary.
@@ -261,33 +269,39 @@ X, cells = dmsh.generate(geo, 0.1)
 
 ### Debugging
 
-![level-set-poly](https://nschloe.github.io/dmsh/levelset-polygon.png) | ![level-set-rect-hole](https://nschloe.github.io/dmsh/levelset-rect-hole.png)
-|:----:|:----:|
+| ![level-set-poly](https://nschloe.github.io/dmsh/levelset-polygon.png) | ![level-set-rect-hole](https://nschloe.github.io/dmsh/levelset-rect-hole.png) |
+| :--------------------------------------------------------------------: | :---------------------------------------------------------------------------: |
 
 dmsh is rather fragile, but sometimes the break-downs are due to an incorrectly defined
 geometry. Use
+
 ```
 geo.show()
 ```
+
 to inspect the level set function of your domain. (It must be negative inside the
 domain and positive outside. The 0-level set forms the domain boundary.)
 
-
 ### Installation
 
 dmsh is [available from the Python Package
 Index](https://pypi.org/project/dmsh/), so simply type
+
 ```
 pip install dmsh
 ```
+
 to install.
 
 ### Testing
+
 To run the dmsh unit tests, check out this repository and type
+
 ```
 tox
 ```
 
 ### License
+
 This software is published under the [MIT
 license](https://en.wikipedia.org/wiki/MIT_License).
diff --git a/pyproject.toml b/pyproject.toml
index 8fe2f47..3c2aae9 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,3 +1,6 @@
 [build-system]
 requires = ["setuptools>=42", "wheel"]
 build-backend = "setuptools.build_meta"
+
+[tool.isort]
+profile = "black"
diff --git a/setup.cfg b/setup.cfg
index 0fe0714..fffa58c 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -32,6 +32,7 @@ packages = find:
 install_requires =
     importlib_metadata;python_version<"3.8"
     meshplex >= 0.16.0, < 0.17.0
+    npx
     numpy
     scipy
 python_requires = >=3.7
diff --git a/src/dmsh/geometry/pypathlib/closed_path.py b/src/dmsh/geometry/pypathlib/closed_path.py
index 82695ce..fd5aa8e 100644
--- a/src/dmsh/geometry/pypathlib/closed_path.py
+++ b/src/dmsh/geometry/pypathlib/closed_path.py
@@ -21,7 +21,7 @@ def __init__(self, points):
         return
 
     def signed_squared_distance(self, x):
-        """Negative inside the polgon."""
+        """Negative inside the polygon."""
         x = np.array(x)
         assert x.shape[1] == 2
         t, dist2, idx = self._all_distances(x)
@@ -30,7 +30,7 @@ def signed_squared_distance(self, x):
         return dist2
 
     def signed_distance(self, x):
-        """Negative inside the polgon."""
+        """Negative inside the polygon."""
         x = np.array(x)
         assert x.shape[1] == 2
         t, dist2, idx = self._all_distances(x)
diff --git a/src/dmsh/geometry/pypathlib/path.py b/src/dmsh/geometry/pypathlib/path.py
index 23dcdbd..2848767 100644
--- a/src/dmsh/geometry/pypathlib/path.py
+++ b/src/dmsh/geometry/pypathlib/path.py
@@ -39,7 +39,7 @@ def _all_distances(self, x):
         #    (<e, e> <x-x0, x-x0> - <x-x0, e>**2) / <e, e>
         #
         # but this expression is numerically disadvantageous. (For example, the
-        # expresison can become negative due to round-off.) Simply compute the
+        # expression can become negative due to round-off.) Simply compute the
         # projection and the dot product.
         x_min_proj = diff - t[:, :, None] * self.edges[None, :, :]
         dist2_sides = np.einsum("ijk,ijk->ij", x_min_proj, x_min_proj)
diff --git a/src/dmsh/helpers.py b/src/dmsh/helpers.py
index 7f5035c..8d4e0b5 100644
--- a/src/dmsh/helpers.py
+++ b/src/dmsh/helpers.py
@@ -95,7 +95,7 @@ def find_feature_points(geometries, num_steps=10):
 
 
 def _find_feature_points_between_two_paths(path0, path1, num_steps, nx=11, ny=11):
-    """Given two geometries with their parameterization, this methods finds feature
+    """Given two geometries with their parametrization, this methods finds feature
     points, i.e., points where the boundaries meet. This is done by casting a net over
     the parameter space and performing `num_steps` Newton steps. Found solutions are
     checked for uniqueness.
diff --git a/src/dmsh/main.py b/src/dmsh/main.py
index fad2c5f..c3527e3 100644
--- a/src/dmsh/main.py
+++ b/src/dmsh/main.py
@@ -2,6 +2,7 @@
 from typing import Callable, Union
 
 import meshplex
+import npx
 import numpy as np
 import scipy.spatial
 
@@ -256,35 +257,45 @@ def distmesh_smoothing(
 
         edges = mesh.edges["points"]
 
-        edges_vec = mesh.points[edges[:, 1]] - mesh.points[edges[:, 0]]
-        edge_lengths = np.sqrt(np.einsum("ij,ij->i", edges_vec, edges_vec))
-        edges_vec /= edge_lengths[..., None]
+        edges_vec_normalized = mesh.points[edges[:, 1]] - mesh.points[edges[:, 0]]
+        edge_lengths = np.sqrt(
+            np.einsum("ij,ij->i", edges_vec_normalized, edges_vec_normalized)
+        )
+        edges_vec_normalized /= edge_lengths[..., None]
 
         # Evaluate element sizes at edge midpoints
         edge_midpoints = (mesh.points[edges[:, 1]] + mesh.points[edges[:, 0]]) / 2
         p = target_edge_size_function(edge_midpoints.T)
-        desired_lengths = (
+        target_lengths = (
             f_scale * p * np.sqrt(np.dot(edge_lengths, edge_lengths) / np.dot(p, p))
         )
 
-        force_abs = desired_lengths - edge_lengths
+        force_abs = target_lengths - edge_lengths
         # only consider repulsive forces
         force_abs[force_abs < 0.0] = 0.0
 
+        # In <https://github.com/nschloe/dmsh/issues/85>, there's a suggestion for a
+        # better forcing function. The below doesn't seem to work too well though.
+        #
+        # Need to set delta_t to 1.0e-2 or smaller to accommodate for the missing factor
+        # `target_lengths`.
+        # force_type = "persson"
+        # relative_length = edge_lengths / target_lengths
+        # if force_type.lower() == "persson":
+        #     force_abs = 1.0 - relative_length
+        #     # only consider repulsive forces
+        #     force_abs[relative_length > 1.0] = 0.0
+        # else:
+        #     assert force_type.lower() == "bossens"
+        #     force_abs = (1 - relative_length ** 4) * np.exp(-(relative_length ** 4))
+
         # force vectors
-        force = edges_vec * force_abs[..., None]
+        force = edges_vec_normalized * force_abs[..., None]
 
-        # bincount replacement for the slow np.add.at
-        # more speed-up can be achieved if the weights were contiguous in memory, i.e.,
-        # if force[k] was used
         n = mesh.points.shape[0]
-        force_per_point = np.array(
-            [
-                np.bincount(edges[:, 0], weights=-force[:, k], minlength=n)
-                + np.bincount(edges[:, 1], weights=+force[:, k], minlength=n)
-                for k in range(force.shape[1])
-            ]
-        ).T
+        force_per_point = npx.sum_at(-force, edges[:, 0], minlength=n) + npx.sum_at(
+            +force, edges[:, 1], minlength=n
+        )
 
         update = delta_t * force_per_point
 
@@ -322,8 +333,6 @@ def distmesh_smoothing(
         if np.all(move2 < tol ** 2):
             break
 
-    # print("num steps:  ", k)
-
     # The cell removal steps in _recell_and_boundary_step() might create points which
     # aren't part of any cell (dangling points). Remove them now.
     mesh.remove_dangling_points()
diff --git a/tox.ini b/tox.ini
index fc2e2c2..95b9534 100644
--- a/tox.ini
+++ b/tox.ini
@@ -9,6 +9,7 @@ deps =
   pytest
   pytest-codeblocks
   pytest-cov
+  pytest-randomly
 extras = all
 commands =
     pytest {posargs} --codeblocks


Test patch:
diff --git a/tests/test_polygon.py b/tests/test_polygon.py
index f538c7b..d5ccaf5 100644
--- a/tests/test_polygon.py
+++ b/tests/test_polygon.py
@@ -24,7 +24,7 @@ def test(show=False):
     return X, cells
 
 
-def test_boundary_step2():
+def test_boundary_step2(plot=False):
     geo = dmsh.Polygon(
         [
             [0.0, 0.0],
@@ -39,11 +39,12 @@ def test_boundary_step2():
     np.random.seed(0)
     pts = np.random.uniform(-2.0, 2.0, (2, 100))
     pts = geo.boundary_step(pts)
-    # geo.plot()
-    # import matplotlib.pyplot as plt
-    # plt.plot(pts[0], pts[1], "xk")
-    # plt.show()
-    # print(geo.dist(pts).shape)
+    if plot:
+        geo.plot()
+        import matplotlib.pyplot as plt
+
+        plt.plot(pts[0], pts[1], "xk")
+        plt.show()
     dist = geo.dist(pts)
     assert np.all(np.abs(dist) < 1.0e-12)
 
@@ -52,4 +53,4 @@ def test_boundary_step2():
     # from helpers import save
     # X, cells = test(show=False)
     # save("polygon.svg", X, cells)
-    test_boundary_step2()
+    test_boundary_step2(plot=True)