diff --git a/MANIFEST.in b/MANIFEST.in index 4022a96..9c1d980 100644 --- a/MANIFEST.in +++ b/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 diff --git a/README.md b/README.md index bb3084e..58541a8 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@

pygalmesh -

Create high-quality 3D meshes with ease.

+

Create high-quality meshes with ease.

[![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 + -* 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 @@ -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). diff --git a/pygalmesh/__init__.py b/pygalmesh/__init__.py index 697dff0..f5bd6dc 100644 --- a/pygalmesh/__init__.py +++ b/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", diff --git a/pygalmesh/main.py b/pygalmesh/main.py index c597be8..cafe383 100644 --- a/pygalmesh/main.py +++ b/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, diff --git a/setup.cfg b/setup.cfg index 02dbc42..210bae3 100644 --- a/setup.cfg +++ b/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 = diff --git a/setup.py b/setup.py index 3b20d21..0cd468d 100644 --- a/setup.py +++ b/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", diff --git a/src/generate_2d.cpp b/src/generate_2d.cpp new file mode 100644 index 0000000..b7596ce --- /dev/null +++ b/src/generate_2d.cpp @@ -0,0 +1,90 @@ +#define CGAL_MESH_3_VERBOSE 1 + +#include "generate_2d.hpp" + +#include +#include +#include +#include +#include +#include +#include + +namespace pygalmesh { + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Delaunay_mesh_vertex_base_2 Vb; +typedef CGAL::Delaunay_mesh_face_base_2 Fb; +typedef CGAL::Triangulation_data_structure_2 Tds; +typedef CGAL::Constrained_Delaunay_triangulation_2 CDT; +typedef CGAL::Delaunay_mesh_size_criteria_2 Criteria; +typedef CDT::Vertex_handle Vertex_handle; +typedef CDT::Point Point; + +std::tuple>, std::vector>> +generate_2d( + const std::vector> & points, + const std::vector> & 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 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_index; + std::vector> 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> 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 diff --git a/src/generate_2d.hpp b/src/generate_2d.hpp new file mode 100644 index 0000000..25bdcce --- /dev/null +++ b/src/generate_2d.hpp @@ -0,0 +1,20 @@ +#ifndef GENERATE_2D_HPP +#define GENERATE_2D_HPP + +#include +#include + +namespace pygalmesh { + +std::tuple>, std::vector>> +generate_2d( + const std::vector> & points, + const std::vector> & constraints, + const double max_circumradius_shortest_edge_ratio, + const double cell_size, + const int num_lloyd_steps +); + +} // namespace pygalmesh + +#endif // GENERATE_2D_HPP diff --git a/src/pybind11.cpp b/src/pybind11.cpp index 73df5df..769de6c 100644 --- a/src/pybind11.cpp +++ b/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"), diff --git a/test/test_2d.py b/test/test_2d.py new file mode 100644 index 0000000..8c0f922 --- /dev/null +++ b/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()