Merge pull request #112 from nschloe/2d

2d mesh generation

Author: nschloe

MANIFEST.in

@@ -3,6 +3,7 @@ include src/generate_from_inr.hpp
 include src/generate_from_off.hpp
 include src/generate_periodic.hpp
 include src/generate.hpp
+include src/generate_2d.hpp
 include src/generate_surface_mesh.hpp
 include src/polygon2d.hpp
 include src/primitives.hpp

README.md

@@ -1,6 +1,6 @@
 <p align="center">
   <a href="https://github.com/nschloe/pygalmesh"><img alt="pygalmesh" src="https://nschloe.github.io/pygalmesh/pygalmesh-logo.svg" width="60%"></a>
-  <p align="center">Create high-quality 3D meshes with ease.</p>
+  <p align="center">Create high-quality meshes with ease.</p>
 </p>
 
 [![PyPi Version](https://img.shields.io/pypi/v/pygalmesh.svg?style=flat-square)](https://pypi.org/project/pygalmesh)
@@ -18,38 +18,29 @@
 [![LGTM](https://img.shields.io/lgtm/grade/python/github/nschloe/pygalmesh.svg?style=flat-square)](https://lgtm.com/projects/g/nschloe/pygalmesh)
 [![Code style: black](https://img.shields.io/badge/code%20style-black-000000.svg?style=flat-square)](https://github.com/psf/black)
 
-pygalmesh is a Python frontend to [CGAL](https://www.cgal.org/)'s [3D mesh generation
-capabilities](https://doc.cgal.org/latest/Mesh_3/index.html).
-pygalmesh makes it easy to create high-quality 3D volume meshes, periodic volume meshes,
-and surface meshes.
+pygalmesh is a Python frontend to [CGAL](https://www.cgal.org/)'s
+[2D](https://doc.cgal.org/latest/Mesh_2/index.html) and [3D mesh generation
+capabilities](https://doc.cgal.org/latest/Mesh_3/index.html).  pygalmesh makes it easy
+to create high-quality 2D, 3D volume meshes, periodic volume meshes, and surface meshes.
 
-### Background
-
-CGAL offers two different approaches for mesh generation:
-
-1. Meshes defined implicitly by level sets of functions.
-2. Meshes defined by a set of bounding planes.
-
-pygalmesh provides a front-end to the first approach, which has the following advantages
-and disadvantages:
-
-* All boundary points are guaranteed to be in the level set within any specified
-  residual. This results in smooth curved surfaces.
-* Sharp intersections of subdomains (e.g., in unions or differences of sets) need to be
-  specified manually (via feature edges, see below), which can be tedious.
+### Examples
 
-On the other hand, the bounding-plane approach (realized by
-[mshr](https://bitbucket.org/fenics-project/mshr)), has the following properties:
+#### 2D meshes
+<img src="https://nschloe.github.io/pygalmesh/rect.svg" width="30%">
 
-* Smooth, curved domains are approximated by a set of bounding planes, resulting in more
-  of less visible edges.
-* Intersections of domains can be computed automatically, so domain unions etc.  have
-  sharp edges where they belong.
+CGAL generates 2D meshes from linear contraints. 
+```python
+import numpy
+import pygalmesh
 
-See [here](https://github.com/nschloe/awesome-scientific-computing#meshing) for other
-mesh generation tools.
+points = numpy.array([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])
+constraints = [[0, 1], [1, 2], [2, 3], [3, 0]]
 
-### Examples
+mesh = pygalmesh.generate_2d(points, constraints, cell_size=1.0e-1, num_lloyd_steps=10)
+# mesh.points, mesh.cells
+```
+The quality of the mesh isn't very good, but can be improved with
+[optimesh](https://github.com/nschloe/optimesh).
 
 #### A simple ball
 <img src="https://nschloe.github.io/pygalmesh/ball.png" width="30%">
@@ -446,6 +437,34 @@ To run the pygalmesh unit tests, check out this repository and type
 pytest
 ```
 
+
+### Background
+
+CGAL offers two different approaches for mesh generation:
+
+1. Meshes defined implicitly by level sets of functions.
+2. Meshes defined by a set of bounding planes.
+
+pygalmesh provides a front-end to the first approach, which has the following advantages
+and disadvantages:
+
+* All boundary points are guaranteed to be in the level set within any specified
+  residual. This results in smooth curved surfaces.
+* Sharp intersections of subdomains (e.g., in unions or differences of sets) need to be
+  specified manually (via feature edges, see below), which can be tedious.
+
+On the other hand, the bounding-plane approach (realized by
+[mshr](https://bitbucket.org/fenics-project/mshr)), has the following properties:
+
+* Smooth, curved domains are approximated by a set of bounding planes, resulting in more
+  of less visible edges.
+* Intersections of domains can be computed automatically, so domain unions etc.  have
+  sharp edges where they belong.
+
+See [here](https://github.com/nschloe/awesome-scientific-computing#meshing) for other
+mesh generation tools.
+
+
 ### License
 
 pygalmesh is published under the [GPLv3 license](https://www.gnu.org/licenses/gpl-3.0.en.html).

pygalmesh/__init__.py

@@ -25,6 +25,7 @@
 from . import _cli
 from .__about__ import __version__
 from .main import (
+    generate_2d,
     generate_from_array,
     generate_from_inr,
     generate_mesh,
@@ -62,6 +63,7 @@
     "RingExtrude",
     #
     "generate_mesh",
+    "generate_2d",
     "generate_periodic_mesh",
     "generate_surface_mesh",
     "generate_volume_mesh_from_surface_mesh",

pygalmesh/main.py

@@ -1,9 +1,12 @@
+import math
 import os
 import tempfile
 
 import meshio
+import numpy
 from _pygalmesh import (
     SizingFieldBase,
+    _generate_2d,
     _generate_from_inr,
     _generate_from_inr_with_subdomain_sizing,
     _generate_from_off,
@@ -67,6 +70,11 @@ def generate_mesh(
     return mesh
 
 
+def generate_2d(points, constraints, B=math.sqrt(2), cell_size=0.0, num_lloyd_steps=0):
+    points, cells = _generate_2d(points, constraints, B, cell_size, num_lloyd_steps)
+    return meshio.Mesh(numpy.array(points), {"triangle": numpy.array(cells)})
+
+
 def generate_periodic_mesh(
     domain,
     bounding_cuboid,

setup.cfg

@@ -1,10 +1,10 @@
 [metadata]
 name = pygalmesh
-version = 0.7.1
+version = 0.7.2
 url = https://github.com/nschloe/pygalmesh
 author = Nico Schlömer
 author_email = nico.schloemer@gmail.com
-description = Python frontend to CGAL's 3D mesh generation capabilities
+description = Python frontend to CGAL's mesh generation capabilities
 long_description = file: README.md
 long_description_content_type = text/markdown
 license = GPL-3.0-or-later
@@ -14,7 +14,6 @@ classifiers =
     License :: OSI Approved :: GNU General Public License v3 or later (GPLv3+)
     Operating System :: OS Independent
     Programming Language :: Python :: 3
-    Programming Language :: Python :: 3.5
     Programming Language :: Python :: 3.6
     Programming Language :: Python :: 3.7
     Programming Language :: Python :: 3.8
@@ -31,7 +30,7 @@ install_requires =
     meshio >= 4.0.0, < 5.0.0
     numpy
     pybind11 >= 2.2
-python_requires = >=3.5
+python_requires = >=3.6
 
 [options.entry_points]
 console_scripts =

setup.py

@@ -23,6 +23,7 @@ def __str__(self):
         "_pygalmesh",
         [
             "src/generate.cpp",
+            "src/generate_2d.cpp",
             "src/generate_from_inr.cpp",
             "src/generate_from_off.cpp",
             "src/generate_periodic.cpp",

src/generate_2d.cpp

@@ -0,0 +1,90 @@
+#define CGAL_MESH_3_VERBOSE 1
+
+#include "generate_2d.hpp"
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Constrained_Delaunay_triangulation_2.h>
+#include <CGAL/Delaunay_mesher_2.h>
+#include <CGAL/Delaunay_mesh_face_base_2.h>
+#include <CGAL/Delaunay_mesh_vertex_base_2.h>
+#include <CGAL/Delaunay_mesh_size_criteria_2.h>
+#include <CGAL/lloyd_optimize_mesh_2.h>
+
+namespace pygalmesh {
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+typedef CGAL::Delaunay_mesh_vertex_base_2<K> Vb;
+typedef CGAL::Delaunay_mesh_face_base_2<K> Fb;
+typedef CGAL::Triangulation_data_structure_2<Vb, Fb> Tds;
+typedef CGAL::Constrained_Delaunay_triangulation_2<K, Tds> CDT;
+typedef CGAL::Delaunay_mesh_size_criteria_2<CDT> Criteria;
+typedef CDT::Vertex_handle Vertex_handle;
+typedef CDT::Point Point;
+
+std::tuple<std::vector<std::array<double, 2>>, std::vector<std::array<int, 3>>>
+generate_2d(
+  const std::vector<std::array<double, 2>> & points,
+  const std::vector<std::array<int, 2>> & constraints,
+  // See
+  // https://doc.cgal.org/latest/Mesh_2/classCGAL_1_1Delaunay__mesh__size__criteria__2.html#a58b0186eae407ba76b8f4a3d0aa85a1a
+  // for what the bounds mean. Spoiler:
+  // B = circumradius / shortest_edge,
+  //   relates to the smallest angle via sin(alpha_min) = 1 / (2B)
+  // cell_size is "size",
+  const double max_circumradius_shortest_edge_ratio,
+  const double cell_size,
+  const int num_lloyd_steps
+)
+{
+  CDT cdt;
+  // construct a constrained triangulation
+  std::vector<Vertex_handle> vertices(points.size());
+  int k = 0;
+  for (auto pt: points) {
+    vertices[k] = cdt.insert(Point(pt[0], pt[1]));
+    k++;
+  }
+  for (auto c: constraints) {
+    cdt.insert_constraint(vertices[c[0]], vertices[c[1]]);
+  }
+
+  // create proper mesh
+  CGAL::refine_Delaunay_mesh_2(
+      cdt,
+      Criteria(
+        0.25 / (max_circumradius_shortest_edge_ratio * max_circumradius_shortest_edge_ratio),
+        cell_size
+       )
+      );
+
+  if (num_lloyd_steps > 0) {
+    CGAL::lloyd_optimize_mesh_2(
+      cdt,
+      CGAL::parameters::max_iteration_number = num_lloyd_steps
+    );
+  }
+
+  // convert points to vector of arrays
+  std::map<Vertex_handle, int> vertex_index;
+  std::vector<std::array<double, 2>> out_points(cdt.number_of_vertices());
+  k = 0;
+  for (auto vit = cdt.vertices_begin(); vit!= cdt.vertices_end(); ++vit) {
+    out_points[k][0] = vit->point()[0];
+    out_points[k][1] = vit->point()[1];
+    vertex_index[vit] = k;
+    k++;
+  }
+
+  std::vector<std::array<int, 3>> out_cells(cdt.number_of_faces());
+  k = 0;
+  for (auto fit = cdt.faces_begin(); fit!= cdt.faces_end(); ++fit) {
+    out_cells[k][0] = vertex_index[fit->vertex(0)];
+    out_cells[k][1] = vertex_index[fit->vertex(1)];
+    out_cells[k][2] = vertex_index[fit->vertex(2)];
+    k++;
+  }
+
+  return std::make_tuple(out_points, out_cells);
+}
+
+} // namespace pygalmesh

src/generate_2d.hpp

@@ -0,0 +1,20 @@
+#ifndef GENERATE_2D_HPP
+#define GENERATE_2D_HPP
+
+#include <memory>
+#include <vector>
+
+namespace pygalmesh {
+
+std::tuple<std::vector<std::array<double, 2>>, std::vector<std::array<int, 3>>>
+generate_2d(
+  const std::vector<std::array<double, 2>> & points,
+  const std::vector<std::array<int, 2>> & constraints,
+  const double max_circumradius_shortest_edge_ratio,
+  const double cell_size,
+  const int num_lloyd_steps
+);
+
+} // namespace pygalmesh
+
+#endif // GENERATE_2D_HPP

src/pybind11.cpp

@@ -1,5 +1,6 @@
 #include "domain.hpp"
 #include "generate.hpp"
+#include "generate_2d.hpp"
 #include "generate_from_off.hpp"
 #include "generate_from_inr.hpp"
 #include "remesh_surface.hpp"
@@ -275,6 +276,14 @@ PYBIND11_MODULE(_pygalmesh, m) {
         py::arg("verbose") = true,
         py::arg("seed") = 0
         );
+    m.def(
+        "_generate_2d", &generate_2d,
+        py::arg("points"),
+        py::arg("constraints"),
+        py::arg("max_circumradius_shortest_edge_ratio") = 1.41421356237,
+        py::arg("cell_size") = 0.0,
+        py::arg("num_lloyd_steps") = 0
+        );
     m.def(
         "_generate_with_sizing_field", &generate_with_sizing_field,
         py::arg("domain"),

test/test_2d.py

@@ -0,0 +1,34 @@
+import numpy
+
+import pygalmesh
+
+
+def test_2d():
+    points = numpy.array([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])
+    constraints = [[0, 1], [1, 2], [2, 3], [3, 0]]
+
+    mesh = pygalmesh.generate_2d(
+        points, constraints, cell_size=1.0e-1, num_lloyd_steps=10
+    )
+
+    assert mesh.points.shape == (276, 2)
+    assert mesh.get_cells_type("triangle").shape == (486, 3)
+
+    # # show mesh
+    # import matplotlib.pyplot as plt
+    # pts = points[cells]
+    # for pt in pts:
+    #     plt.plot([pt[0][0], pt[1][0]], [pt[0][1], pt[1][1]], "-k")
+    #     plt.plot([pt[1][0], pt[2][0]], [pt[1][1], pt[2][1]], "-k")
+    #     plt.plot([pt[2][0], pt[0][0]], [pt[2][1], pt[0][1]], "-k")
+    # # for pt in points:
+    # #     plt.plot(pt[0], pt[1], "or")
+    # plt.gca().set_aspect("equal")
+    # plt.show()
+
+    # mesh.points *= 100
+    # mesh.write("rect.svg")
+
+
+if __name__ == "__main__":
+    test_2d()