Merge pull request #116 from nschloe/new-generate

better generate

Author: nschloe

README.md

@@ -36,7 +36,7 @@ import pygalmesh
 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)
+mesh = pygalmesh.generate_2d(points, constraints, edge_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
@@ -213,24 +213,19 @@ to CGAL's mesh generator.
 #### Local refinement
 <img src="https://nschloe.github.io/pygalmesh/ball-local-refinement.png" width="30%">
 
-Use `generate_mesh` with a `SizingFieldBase` object as `cell_size`.
+Use `generate_mesh` with a function (regular or lambda) as `cell_size`. The same goes
+for `edge_size`, `facet_size`, and `facet_distance`.
 ```python
 import numpy
 import pygalmesh
 
-# define a cell_size function
-class Field(pygalmesh.SizingFieldBase):
-    def eval(self, x):
-        return abs(numpy.sqrt(numpy.dot(x, x)) - 0.5) / 5 + 0.025
-
-
 mesh = pygalmesh.generate_mesh(
     pygalmesh.Ball([0.0, 0.0, 0.0], 1.0),
     facet_angle=30,
     facet_size=0.1,
     facet_distance=0.025,
     cell_radius_edge_ratio=2,
-    cell_size=Field(),
+    cell_size=lambda x: abs(numpy.sqrt(numpy.dot(x, x)) - 0.5) / 5 + 0.025,
 )
 ```
 

pygalmesh/__init__.py

@@ -14,7 +14,6 @@
     RingExtrude,
     Rotate,
     Scale,
-    SizingFieldBase,
     Stretch,
     Tetrahedron,
     Torus,
@@ -42,7 +41,6 @@
     "_cli",
     #
     "DomainBase",
-    "SizingFieldBase",
     "Translate",
     "Rotate",
     "Scale",

pygalmesh/main.py

@@ -13,11 +13,19 @@
     _generate_mesh,
     _generate_periodic_mesh,
     _generate_surface_mesh,
-    _generate_with_sizing_field,
     _remesh_surface,
 )
 
 
+class Wrapper(SizingFieldBase):
+    def __init__(self, f):
+        self.f = f
+        super().__init__()
+
+    def eval(self, x):
+        return self.f(x)
+
+
 def generate_mesh(
     domain,
     feature_edges=None,
@@ -40,13 +48,18 @@ def generate_mesh(
     fh, outfile = tempfile.mkstemp(suffix=".mesh")
     os.close(fh)
 
-    if isinstance(cell_size, float):
-        genfun = _generate_mesh
-    else:
-        assert isinstance(cell_size, SizingFieldBase)
-        genfun = _generate_with_sizing_field
+    def _select(obj):
+        if isinstance(obj, float):
+            return obj, None
+        assert callable(obj)
+        return -1.0, Wrapper(obj)
 
-    genfun(
+    edge_size_value, edge_size_field = _select(edge_size)
+    cell_size_value, cell_size_field = _select(cell_size)
+    facet_size_value, facet_size_field = _select(facet_size)
+    facet_distance_value, facet_distance_field = _select(facet_distance)
+
+    _generate_mesh(
         domain,
         outfile,
         feature_edges=feature_edges,
@@ -55,12 +68,16 @@ def generate_mesh(
         odt=odt,
         perturb=perturb,
         exude=exude,
-        edge_size=edge_size,
+        edge_size_value=edge_size_value,
+        edge_size_field=edge_size_field,
         facet_angle=facet_angle,
-        facet_size=facet_size,
-        facet_distance=facet_distance,
+        facet_size_value=facet_size_value,
+        facet_size_field=facet_size_field,
+        facet_distance_value=facet_distance_value,
+        facet_distance_field=facet_distance_field,
         cell_radius_edge_ratio=cell_radius_edge_ratio,
-        cell_size=cell_size,
+        cell_size_value=cell_size_value,
+        cell_size_field=cell_size_field,
         verbose=verbose,
         seed=seed,
     )
@@ -70,7 +87,7 @@ def generate_mesh(
     return mesh
 
 
-def generate_2d(points, constraints, B=math.sqrt(2), cell_size=0.0, num_lloyd_steps=0):
+def generate_2d(points, constraints, B=math.sqrt(2), edge_size=0.0, num_lloyd_steps=0):
     # some sanity checks
     points = numpy.asarray(points)
     constraints = numpy.asarray(constraints)
@@ -82,7 +99,7 @@ def generate_2d(points, constraints, B=math.sqrt(2), cell_size=0.0, num_lloyd_st
     if numpy.any(length2 < 1.0e-15):
         raise RuntimeError("Constraint of (near)-zero length.")
 
-    points, cells = _generate_2d(points, constraints, B, cell_size, num_lloyd_steps)
+    points, cells = _generate_2d(points, constraints, B, edge_size, num_lloyd_steps)
     return meshio.Mesh(numpy.array(points), {"triangle": numpy.array(cells)})
 
 

setup.cfg

@@ -1,6 +1,6 @@
 [metadata]
 name = pygalmesh
-version = 0.7.2
+version = 0.8.0
 url = https://github.com/nschloe/pygalmesh
 author = Nico Schlömer
 author_email = nico.schloemer@gmail.com
@@ -21,6 +21,13 @@ classifiers =
     Topic :: Scientific/Engineering :: Mathematics
     Topic :: Scientific/Engineering :: Physics
     Topic :: Scientific/Engineering :: Visualization
+keywords =
+    mathematics
+    physics
+    engineering
+    cgal
+    mesh
+    mesh generation
 
 [options]
 packages = find:

src/generate.cpp

@@ -24,6 +24,7 @@ typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
 
 // Mesh Criteria
 typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
+typedef Mesh_criteria::Edge_criteria Edge_criteria;
 typedef Mesh_criteria::Facet_criteria Facet_criteria;
 typedef Mesh_criteria::Cell_criteria Cell_criteria;
 
@@ -54,12 +55,23 @@ generate_mesh(
     const bool odt,
     const bool perturb,
     const bool exude,
-    const double edge_size,
+    //
+    const double edge_size_value,
+    const std::shared_ptr<pygalmesh::SizingFieldBase> & edge_size_field,
+    //
     const double facet_angle,
-    const double facet_size,
-    const double facet_distance,
+    //
+    const double facet_size_value,
+    const std::shared_ptr<pygalmesh::SizingFieldBase> & facet_size_field,
+    //
+    const double facet_distance_value,
+    const std::shared_ptr<pygalmesh::SizingFieldBase> & facet_distance_field,
+    //
     const double cell_radius_edge_ratio,
-    const double cell_size,
+    //
+    const double cell_size_value,
+    const std::shared_ptr<pygalmesh::SizingFieldBase> & cell_size_field,
+    //
     const bool verbose,
     const int seed
     )
@@ -81,26 +93,73 @@ generate_mesh(
        K::Sphere_3(CGAL::ORIGIN, bounding_sphere_radius2)
        );
 
+  // cgal_domain.detect_features();
+
   const auto native_features = translate_feature_edges(domain->get_features());
   cgal_domain.add_features(native_features.begin(), native_features.end());
 
   const auto polylines = translate_feature_edges(feature_edges);
   cgal_domain.add_features(polylines.begin(), polylines.end());
 
-  Mesh_criteria criteria(
-      CGAL::parameters::edge_size=edge_size,
-      CGAL::parameters::facet_angle=facet_angle,
-      CGAL::parameters::facet_size=facet_size,
-      CGAL::parameters::facet_distance=facet_distance,
-      CGAL::parameters::cell_radius_edge_ratio=cell_radius_edge_ratio,
-      CGAL::parameters::cell_size=cell_size
-      );
-
-  // Mesh generation
+  // perhaps there's a more elegant solution here
+  // see <https://github.com/CGAL/cgal/issues/4145>
   if (!verbose) {
     // suppress output
     std::cerr.setstate(std::ios_base::failbit);
   }
+
+  // Build the float/field values according to
+  // <https://github.com/CGAL/cgal/issues/5044#issuecomment-705526982>.
+
+  // nested ternary operator
+  const auto facet_criteria = facet_size_field ? (
+      facet_distance_field ?
+      Facet_criteria(
+        facet_angle,
+         [&](K::Point_3 p, const int, const Mesh_domain::Index&) {
+           return facet_size_field->eval({p.x(), p.y(), p.z()});
+         },
+         [&](K::Point_3 p, const int, const Mesh_domain::Index&) {
+           return facet_distance_field->eval({p.x(), p.y(), p.z()});
+         }
+      ) : Facet_criteria(
+        facet_angle,
+         [&](K::Point_3 p, const int, const Mesh_domain::Index&) {
+           return facet_size_field->eval({p.x(), p.y(), p.z()});
+         },
+         facet_distance_value
+      )
+    ) : (
+      facet_distance_field ?
+      Facet_criteria(
+        facet_angle,
+        facet_size_value,
+         [&](K::Point_3 p, const int, const Mesh_domain::Index&) {
+           return facet_distance_field->eval({p.x(), p.y(), p.z()});
+         }
+      ) : Facet_criteria(
+        facet_angle,
+        facet_size_value,
+        facet_distance_value
+      )
+    );
+
+  const auto edge_criteria = edge_size_field ?
+     Edge_criteria(
+         [&](K::Point_3 p, const int, const Mesh_domain::Index&) {
+           return edge_size_field->eval({p.x(), p.y(), p.z()});
+          }) : Edge_criteria(edge_size_value);
+
+  const auto cell_criteria = cell_size_field ?
+     Cell_criteria(
+         cell_radius_edge_ratio,
+         [&](K::Point_3 p, const int, const Mesh_domain::Index&) {
+           return cell_size_field->eval({p.x(), p.y(), p.z()});
+          }) : Cell_criteria(cell_radius_edge_ratio, cell_size_value);
+
+  const auto criteria = Mesh_criteria(edge_criteria, facet_criteria, cell_criteria);
+
+  // Mesh generation
   C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(
       cgal_domain,
       criteria,
@@ -121,115 +180,4 @@ generate_mesh(
   return;
 }
 
-void
-generate_with_sizing_field(
-    const std::shared_ptr<pygalmesh::DomainBase> & domain,
-    const std::string & outfile,
-    const std::vector<std::vector<std::array<double, 3>>> & feature_edges,
-    const double bounding_sphere_radius,
-    const bool lloyd,
-    const bool odt,
-    const bool perturb,
-    const bool exude,
-    const double edge_size,
-    const double facet_angle,
-    const double facet_size,
-    const double facet_distance,
-    const double cell_radius_edge_ratio,
-    const std::shared_ptr<pygalmesh::SizingFieldBase> & cell_size,
-    const bool verbose,
-    const int seed
-    )
-{
-  CGAL::get_default_random() = CGAL::Random(seed);
-
-  const double bounding_sphere_radius2 = bounding_sphere_radius > 0 ?
-    bounding_sphere_radius*bounding_sphere_radius :
-    // some wiggle room
-    1.01 * domain->get_bounding_sphere_squared_radius();
-
-  // wrap domain
-  const auto d = [&](K::Point_3 p) {
-    return domain->eval({p.x(), p.y(), p.z()});
-  };
-
-  Mesh_domain cgal_domain = Mesh_domain::create_implicit_mesh_domain(
-       d,
-       K::Sphere_3(CGAL::ORIGIN, bounding_sphere_radius2)
-       );
-
-  // cgal_domain.detect_features();
-
-  const auto native_features = translate_feature_edges(domain->get_features());
-  cgal_domain.add_features(native_features.begin(), native_features.end());
-
-  const auto polylines = translate_feature_edges(feature_edges);
-  cgal_domain.add_features(polylines.begin(), polylines.end());
-
-  // perhaps there's a more elegant solution here
-  // see <https://github.com/CGAL/cgal/issues/4145>
-  if (!verbose) {
-    // suppress output
-    std::cerr.setstate(std::ios_base::failbit);
-  }
-  if (cell_size) {
-    const auto size = [&](K::Point_3 p, const int, const Mesh_domain::Index&) {
-        return cell_size->eval({p.x(), p.y(), p.z()});
-    };
-    auto criteria = Mesh_criteria(
-        CGAL::parameters::edge_size=edge_size,
-        CGAL::parameters::facet_angle=facet_angle,
-        CGAL::parameters::facet_size=facet_size,
-        CGAL::parameters::facet_distance=facet_distance,
-        CGAL::parameters::cell_radius_edge_ratio=cell_radius_edge_ratio,
-        CGAL::parameters::cell_size=size
-        );
-
-    // Mesh generation
-    C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(
-        cgal_domain,
-        criteria,
-        lloyd ? CGAL::parameters::lloyd() : CGAL::parameters::no_lloyd(),
-        odt ? CGAL::parameters::odt() : CGAL::parameters::no_odt(),
-        perturb ? CGAL::parameters::perturb() : CGAL::parameters::no_perturb(),
-        exude ? CGAL::parameters::exude() : CGAL::parameters::no_exude()
-        );
-    if (!verbose) {
-      std::cerr.clear();
-    }
-
-    // Output
-    std::ofstream medit_file(outfile);
-    c3t3.output_to_medit(medit_file);
-    medit_file.close();
-  } else {
-    auto criteria = Mesh_criteria(
-        CGAL::parameters::edge_size=edge_size,
-        CGAL::parameters::facet_angle=facet_angle,
-        CGAL::parameters::facet_size=facet_size,
-        CGAL::parameters::facet_distance=facet_distance,
-        CGAL::parameters::cell_radius_edge_ratio=cell_radius_edge_ratio
-        );
-
-    // Mesh generation
-    C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(
-        cgal_domain,
-        criteria,
-        lloyd ? CGAL::parameters::lloyd() : CGAL::parameters::no_lloyd(),
-        odt ? CGAL::parameters::odt() : CGAL::parameters::no_odt(),
-        perturb ? CGAL::parameters::perturb() : CGAL::parameters::no_perturb(),
-        exude ? CGAL::parameters::exude() : CGAL::parameters::no_exude()
-        );
-    if (!verbose) {
-      std::cerr.clear();
-    }
-
-    // Output
-    std::ofstream medit_file(outfile);
-    c3t3.output_to_medit(medit_file);
-    medit_file.close();
-  }
-  return;
-}
-
 } // namespace pygalmesh