Merge pull request #130 from nschloe/better-tetrahedron

Better tetrahedron

.github/workflows/ci.yml

@@ -3,10 +3,10 @@ name: ci
 on:
   push:
     branches:
-    - master
+    - main
   pull_request:
     branches:
-    - master
+    - main
 
 jobs:
   lint:

Makefile

@@ -4,13 +4,13 @@ default:
 	@echo "\"make publish\"?"
 
 tag:
-	@if [ "$(shell git rev-parse --abbrev-ref HEAD)" != "master" ]; then exit 1; fi
+	@if [ "$(shell git rev-parse --abbrev-ref HEAD)" != "main" ]; then exit 1; fi
 	curl -H "Authorization: token `cat $(HOME)/.github-access-token`" -d '{"tag_name": "v$(VERSION)"}' https://github.com/repos/nschloe/pygalmesh/releases
 
 upload: clean
-	# Make sure we're on the master branch
-	@if [ "$(shell git rev-parse --abbrev-ref HEAD)" != "master" ]; then exit 1; fi
-	python3 setup.py sdist
+	# Make sure we're on the main branch
+	@if [ "$(shell git rev-parse --abbrev-ref HEAD)" != "main" ]; then exit 1; fi
+	python3 -m build --sdist .
 	twine upload dist/*.tar.gz
 	# HTTPError: 400 Client Error: Binary wheel 'pygalmesh-0.2.0-cp27-cp27mu-linux_x86_64.whl' has an unsupported platform tag 'linux_x86_64'. for url: https://upload.pypi.org/legacy/
 	# python3 setup.py bdist_wheel

pygalmesh/main.py

@@ -29,19 +29,19 @@ def eval(self, x):
 def generate_mesh(
     domain,
     extra_feature_edges=None,
-    bounding_sphere_radius=0.0,
-    lloyd=False,
-    odt=False,
-    perturb=True,
-    exude=True,
-    max_edge_size_at_feature_edges=0.0,
-    min_facet_angle=0.0,
-    max_radius_surface_delaunay_ball=0.0,
-    max_facet_distance=0.0,
-    max_circumradius_edge_ratio=0.0,
-    max_cell_circumradius=0.0,
-    verbose=True,
-    seed=0,
+    bounding_sphere_radius: float = 0.0,
+    lloyd: bool = False,
+    odt: bool = False,
+    perturb: bool = True,
+    exude: bool = True,
+    max_edge_size_at_feature_edges: float = 0.0,
+    min_facet_angle: float = 0.0,
+    max_radius_surface_delaunay_ball: float = 0.0,
+    max_facet_distance: float = 0.0,
+    max_circumradius_edge_ratio: float = 0.0,
+    max_cell_circumradius: float = 0.0,
+    verbose: bool = True,
+    seed: int = 0,
 ):
     """
     From <https://doc.cgal.org/latest/Mesh_3/classCGAL_1_1Mesh__criteria__3.html>:
@@ -130,9 +130,9 @@ def _select(obj):
 def generate_2d(
     points,
     constraints,
-    B=math.sqrt(2),
-    max_edge_size=0.0,
-    num_lloyd_steps=0,
+    B: float = math.sqrt(2),
+    max_edge_size: float = 0.0,
+    num_lloyd_steps: int = 0,
 ):
     # some sanity checks
     points = numpy.asarray(points)
@@ -158,19 +158,19 @@ def generate_2d(
 def generate_periodic_mesh(
     domain,
     bounding_cuboid,
-    lloyd=False,
-    odt=False,
-    perturb=True,
-    exude=True,
-    max_edge_size_at_feature_edges=0.0,
-    min_facet_angle=0.0,
-    max_radius_surface_delaunay_ball=0.0,
-    max_facet_distance=0.0,
-    max_circumradius_edge_ratio=0.0,
-    max_cell_circumradius=0.0,
-    number_of_copies_in_output=1,
-    verbose=True,
-    seed=0,
+    lloyd: bool = False,
+    odt: bool = False,
+    perturb: bool = True,
+    exude: bool = True,
+    max_edge_size_at_feature_edges: float = 0.0,
+    min_facet_angle: float = 0.0,
+    max_radius_surface_delaunay_ball: float = 0.0,
+    max_facet_distance: float = 0.0,
+    max_circumradius_edge_ratio: float = 0.0,
+    max_cell_circumradius: float = 0.0,
+    number_of_copies_in_output: int = 1,
+    verbose: bool = True,
+    seed: int = 0,
 ):
     fh, outfile = tempfile.mkstemp(suffix=".mesh")
     os.close(fh)
@@ -203,12 +203,12 @@ def generate_periodic_mesh(
 
 def generate_surface_mesh(
     domain,
-    bounding_sphere_radius=0.0,
-    min_facet_angle=0.0,
-    max_radius_surface_delaunay_ball=0.0,
-    max_facet_distance=0.0,
-    verbose=True,
-    seed=0,
+    bounding_sphere_radius: float = 0.0,
+    min_facet_angle: float = 0.0,
+    max_radius_surface_delaunay_ball: float = 0.0,
+    max_facet_distance: float = 0.0,
+    verbose: bool = True,
+    seed: int = 0,
 ):
     fh, outfile = tempfile.mkstemp(suffix=".off")
     os.close(fh)
@@ -230,20 +230,20 @@ def generate_surface_mesh(
 
 
 def generate_volume_mesh_from_surface_mesh(
-    filename,
-    lloyd=False,
-    odt=False,
-    perturb=True,
-    exude=True,
-    max_edge_size_at_feature_edges=0.0,
-    min_facet_angle=0.0,
-    max_radius_surface_delaunay_ball=0.0,
-    max_facet_distance=0.0,
-    max_circumradius_edge_ratio=0.0,
-    max_cell_circumradius=0.0,
-    verbose=True,
-    reorient=False,
-    seed=0,
+    filename: str,
+    lloyd: bool = False,
+    odt: bool = False,
+    perturb: bool = True,
+    exude: bool = True,
+    max_edge_size_at_feature_edges: float = 0.0,
+    min_facet_angle: float = 0.0,
+    max_radius_surface_delaunay_ball: float = 0.0,
+    max_facet_distance: float = 0.0,
+    max_circumradius_edge_ratio: float = 0.0,
+    max_cell_circumradius: float = 0.0,
+    verbose: bool = True,
+    reorient: bool = False,
+    seed: int = 0,
 ):
     mesh = meshio.read(filename)
 
@@ -279,19 +279,19 @@ def generate_volume_mesh_from_surface_mesh(
 
 
 def generate_from_inr(
-    inr_filename,
-    lloyd=False,
-    odt=False,
-    perturb=True,
-    exude=True,
-    max_edge_size_at_feature_edges=0.0,
-    min_facet_angle=0.0,
-    max_radius_surface_delaunay_ball=0.0,
-    max_facet_distance=0.0,
-    max_circumradius_edge_ratio=0.0,
-    max_cell_circumradius=0.0,
-    verbose=True,
-    seed=0,
+    inr_filename: str,
+    lloyd: bool = False,
+    odt: bool = False,
+    perturb: bool = True,
+    exude: bool = True,
+    max_edge_size_at_feature_edges: float = 0.0,
+    min_facet_angle: float = 0.0,
+    max_radius_surface_delaunay_ball: float = 0.0,
+    max_facet_distance: float = 0.0,
+    max_circumradius_edge_ratio: float = 0.0,
+    max_cell_circumradius: float = 0.0,
+    verbose: bool = True,
+    seed: int = 0,
 ):
     fh, outfile = tempfile.mkstemp(suffix=".mesh")
     os.close(fh)
@@ -348,13 +348,13 @@ def generate_from_inr(
 
 
 def remesh_surface(
-    filename,
-    max_edge_size_at_feature_edges=0.0,
-    min_facet_angle=0.0,
-    max_radius_surface_delaunay_ball=0.0,
-    max_facet_distance=0.0,
-    verbose=True,
-    seed=0,
+    filename: str,
+    max_edge_size_at_feature_edges: float = 0.0,
+    min_facet_angle: float = 0.0,
+    max_radius_surface_delaunay_ball: float = 0.0,
+    max_facet_distance: float = 0.0,
+    verbose: bool = True,
+    seed: int = 0,
 ):
     mesh = meshio.read(filename)
 
@@ -382,7 +382,7 @@ def remesh_surface(
     return mesh
 
 
-def save_inr(vol, h, fname):
+def save_inr(vol, h, fname: str):
     """
     Save a volume (described as a numpy array) to INR format.
     Code inspired by iso2mesh (http://iso2mesh.sf.net) by Q. Fang
@@ -414,18 +414,18 @@ def save_inr(vol, h, fname):
 def generate_from_array(
     vol,
     h,
-    lloyd=False,
-    odt=False,
-    perturb=True,
-    exude=True,
-    max_edge_size_at_feature_edges=0.0,
-    min_facet_angle=0.0,
-    max_radius_surface_delaunay_ball=0.0,
-    max_cell_circumradius=0.0,
-    max_facet_distance=0.0,
-    max_circumradius_edge_ratio=0.0,
-    verbose=True,
-    seed=0,
+    lloyd: bool = False,
+    odt: bool = False,
+    perturb: bool = True,
+    exude: bool = True,
+    max_edge_size_at_feature_edges: float = 0.0,
+    min_facet_angle: float = 0.0,
+    max_radius_surface_delaunay_ball: float = 0.0,
+    max_cell_circumradius: float = 0.0,
+    max_facet_distance: float = 0.0,
+    max_circumradius_edge_ratio: float = 0.0,
+    verbose: bool = True,
+    seed: int = 0,
 ):
     assert vol.dtype in ["uint8", "uint16"]
     fh, inr_filename = tempfile.mkstemp(suffix=".inr")

setup.cfg

@@ -1,6 +1,6 @@
 [metadata]
 name = pygalmesh
-version = 0.9.2
+version = 0.9.3
 url = https://github.com/nschloe/pygalmesh
 author = Nico Schlömer
 author_email = nico.schloemer@gmail.com
@@ -17,6 +17,7 @@ classifiers =
     Programming Language :: Python :: 3.6
     Programming Language :: Python :: 3.7
     Programming Language :: Python :: 3.8
+    Programming Language :: Python :: 3.9
     Topic :: Scientific/Engineering
     Topic :: Scientific/Engineering :: Mathematics
     Topic :: Scientific/Engineering :: Physics

src/primitives.hpp

@@ -306,43 +306,62 @@ class Tetrahedron: public pygalmesh::DomainBase
         const std::array<double, 3> & x2,
         const std::array<double, 3> & x3
         ):
-      x0_(Eigen::Vector3d(x0[0], x0[1], x0[2])),
-      x1_(Eigen::Vector3d(x1[0], x1[1], x1[2])),
-      x2_(Eigen::Vector3d(x2[0], x2[1], x2[2])),
-      x3_(Eigen::Vector3d(x3[0], x3[1], x3[2]))
+      x0_(Eigen::Vector3d(x0.data())),
+      x1_(Eigen::Vector3d(x1.data())),
+      x2_(Eigen::Vector3d(x2.data())),
+      x3_(Eigen::Vector3d(x3.data())),
+      A_(constructA(x0, x1, x2, x3))
     {
     }
 
+    Eigen::Matrix4d constructA(
+        const std::array<double, 3> & x0,
+        const std::array<double, 3> & x1,
+        const std::array<double, 3> & x2,
+        const std::array<double, 3> & x3
+    ) {
+      Eigen::Matrix4d A;
+      A << x0[0], x1[0], x2[0], x3[0],
+           x0[1], x1[1], x2[1], x3[1],
+           x0[2], x1[2], x2[2], x3[2],
+           1.0,   1.0,   1.0,   1.0;
+      return A;
+    }
+
     virtual ~Tetrahedron() = default;
 
-    bool isOnSameSide(
-        const Eigen::Vector3d & v0,
-        const Eigen::Vector3d & v1,
-        const Eigen::Vector3d & v2,
-        const Eigen::Vector3d & v3,
-        const Eigen::Vector3d & p
-        ) const
-    {
-      const auto normal = (v1 - v0).cross(v2 - v0);
-      const double dot_v3 = normal.dot(v3 - v0);
-      const double dot_p = normal.dot(p - v0);
-      return (
-          (dot_v3 > 0 && dot_p > 0) || (dot_v3 < 0 && dot_p < 0)
-          );
-    }
+    // bool isOnSameSide(
+    //     const Eigen::Vector3d & v0,
+    //     const Eigen::Vector3d & v1,
+    //     const Eigen::Vector3d & v2,
+    //     const Eigen::Vector3d & v3,
+    //     const Eigen::Vector3d & p
+    //     ) const
+    // {
+    //   const auto normal = (v1 - v0).cross(v2 - v0);
+    //   const double dot_v3 = normal.dot(v3 - v0);
+    //   const double dot_p = normal.dot(p - v0);
+    //   return (
+    //       (dot_v3 > 0 && dot_p > 0) || (dot_v3 < 0 && dot_p < 0)
+    //       );
+    // }
 
     virtual
     double
     eval(const std::array<double, 3> & x) const
     {
-      // TODO continuous expression
-      Eigen::Vector3d pvec(x.data());
-      const bool a =
-        isOnSameSide(x0_, x1_, x2_, x3_, pvec) &&
-        isOnSameSide(x1_, x2_, x3_, x0_, pvec) &&
-        isOnSameSide(x2_, x3_, x0_, x1_, pvec) &&
-        isOnSameSide(x3_, x0_, x1_, x2_, pvec);
-      return a ? -1.0 : 1.0;
+      Eigen::Vector4d b;
+      b << x[0], x[1], x[2], 1.0;
+      Eigen::Vector4d bary = A_.partialPivLu().solve(b);
+      return -bary.minCoeff();
+
+      // Eigen::Vector3d pvec(x.data());
+      // const bool a =
+      //   isOnSameSide(x0_, x1_, x2_, x3_, pvec) &&
+      //   isOnSameSide(x1_, x2_, x3_, x0_, pvec) &&
+      //   isOnSameSide(x2_, x3_, x0_, x1_, pvec) &&
+      //   isOnSameSide(x3_, x0_, x1_, x2_, pvec);
+      // return a ? -1.0 : 1.0;
     }
 
     virtual
@@ -382,6 +401,7 @@ class Tetrahedron: public pygalmesh::DomainBase
     const Eigen::Vector3d x1_;
     const Eigen::Vector3d x2_;
     const Eigen::Vector3d x3_;
+    const Eigen::Matrix4d A_;
 };