--- /dev/null
+version: 2
+
+jobs:
+
+ build:
+ working_directory: ~/pygalmesh
+ docker:
+ - image: ubuntu:18.04
+ steps:
+ - run: apt-get update
+ - run: apt-get install -y libcgal-dev libeigen3-dev python3-pip clang++-6.0
+ - run: pip3 install -U pytest pytest-cov black flake8 meshio pybind11
+ - checkout
+ # install
+ # Use clang++ for its smaller memory footprint.
+ - run: CC="clang++-6.0" pip3 install .
+ # lint
+ - run: LC_ALL=C.UTF-8 black --check setup.py pygalmesh/ test/*.py
+ - run: flake8 setup.py pygalmesh/ test/*.py
+ # The actual test
+ - run: cd test/ && pytest
--- /dev/null
+[flake8]
+ignore = E203, E266, E501, W503
+max-line-length = 80
+max-complexity = 18
+select = B,C,E,F,W,T4,B9
--- /dev/null
+*.off
+*.vtu
+*.mesh
+.cache/
+build/
+dist/
+MANIFEST
+README.rst
+do-configure.sh
+.pytest_cache/
+*.egg-info/
--- /dev/null
+[MESSAGES CONTROL]
+
+disable=
+ fixme,
+ invalid-name,
+ missing-docstring,
+ no-member,
+ too-few-public-methods
--- /dev/null
+dist: trusty
+
+language: python
+
+python:
+ - '2.7'
+ - '3.6'
+
+addons:
+ apt:
+ sources:
+ - sourceline: 'ppa:nschloe/cgal-backports'
+ - sourceline: 'ppa:nschloe/eigen-backports'
+ - sourceline: 'ppa:nschloe/swig-backports'
+ packages:
+ - libcgal-dev
+ - libcgal-qt5-11 # bug in cgal, should be installed automatically
+ - libeigen3-dev
+ - pandoc
+ - python-numpy
+ - swig
+
+# test dependencies
+before_install:
+ - pip install meshio
+
+install:
+ - pip install pytest
+ - pip install .
+
+script:
+ - cd test
+ - curl https://raw.githubusercontent.com/libigl/libigl-unit-tests/master/data/elephant.off > elephant.off
+ - pytest
--- /dev/null
+cmake_minimum_required(VERSION 3.0)
+
+project(pygalmesh CXX)
+
+add_subdirectory(src)
--- /dev/null
+The MIT License (MIT)
+
+Copyright (c) 2016-2018 Nico Schlömer
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
--- /dev/null
+include src/domain.hpp
+include src/generate_from_off.hpp
+include src/generate.hpp
+include src/generate_surface_mesh.hpp
+include src/polygon2d.hpp
+include src/primitives.hpp
--- /dev/null
+VERSION=$(shell python3 -c "import pygalmesh; print(pygalmesh.__version__)")
+
+default:
+ @echo "\"make publish\"?"
+
+tag:
+ @if [ "$(shell git rev-parse --abbrev-ref HEAD)" != "master" ]; then exit 1; fi
+ @echo "Tagging v$(VERSION)..."
+ git tag v$(VERSION)
+ git push --tags
+
+upload: setup.py
+ # Make sure we're on the master branch
+ @if [ "$(shell git rev-parse --abbrev-ref HEAD)" != "master" ]; then exit 1; fi
+ rm -rf dist/*
+ python3 setup.py 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 --universal
+ # twine upload dist/*.whl
+
+publish: tag upload
+
+clean:
+ @find . | grep -E "(__pycache__|\.pyc|\.pyo$\)" | xargs rm -rf
+
+black:
+ black setup.py pygalmesh/ test/*.py
+
+lint:
+ black --check setup.py pygalmesh/ test/*.py
+ flake8 setup.py pygalmesh/ test/*.py
--- /dev/null
+# pygalmesh
+
+A Python frontend to [CGAL](https://www.cgal.org/)'s [3D mesh generation
+capabilities](https://doc.cgal.org/latest/Mesh_3/index.html).
+
+[](https://circleci.com/gh/nschloe/pygalmesh/tree/master)
+[](https://github.com/ambv/black)
+[](https://pypi.org/project/pygalmesh)
+[](https://github.com/nschloe/pygalmesh)
+
+pygalmesh makes it easy to create high-quality 3D volume 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.
+
+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.
+
+Other Python mesh generators are [pygmsh](https://github.com/nschloe/pygmsh) (a
+frontend to [gmsh](http://gmsh.info/)) and
+[MeshPy](https://github.com/inducer/meshpy).
+[meshzoo](https://github.com/nschloe/meshzoo) provides some basic canonical
+meshes.
+
+### Examples
+
+#### A simple ball
+<img src="https://nschloe.github.io/pygalmesh/ball.png" width="30%">
+
+```python
+import pygalmesh
+
+s = pygalmesh.Ball([0, 0, 0], 1.0)
+pygalmesh.generate_mesh(s, "out.mesh", cell_size=0.2)
+```
+CGAL's mesh generator returns Medit-files, which can be processed by, e.g.,
+[meshio](https://github.com/nschloe/meshio).
+```python
+import meshio
+vertices, cells, _, _, _ = meshio.read("out.mesh")
+```
+The mesh generation comes with many more options, described
+[here](https://doc.cgal.org/latest/Mesh_3/). Try, for example,
+```python
+pygalmesh.generate_mesh(
+ s,
+ "out.mesh",
+ cell_size=0.2,
+ edge_size=0.1,
+ odt=True,
+ lloyd=True,
+ verbose=False
+)
+```
+
+#### Other primitive shapes
+<img src="https://nschloe.github.io/pygalmesh/tetra.png" width="30%">
+
+pygalmesh provides out-of-the-box support for balls, cuboids, ellipsoids, tori,
+cones, cylinders, and tetrahedra. Try for example
+```python
+import pygalmesh
+
+s0 = pygalmesh.Tetrahedron(
+ [0.0, 0.0, 0.0],
+ [1.0, 0.0, 0.0],
+ [0.0, 1.0, 0.0],
+ [0.0, 0.0, 1.0]
+)
+pygalmesh.generate_mesh(s0, "out.mesh", cell_size=0.1, edge_size=0.1)
+```
+
+#### Domain combinations
+<img src="https://nschloe.github.io/pygalmesh/ball-difference.png" width="30%">
+
+Supported are unions, intersections, and differences of all domains. As
+mentioned above, however, the sharp intersections between two domains are not
+automatically handled. Try for example
+```python
+import pygalmesh
+
+radius = 1.0
+displacement = 0.5
+s0 = pygalmesh.Ball([displacement, 0, 0], radius)
+s1 = pygalmesh.Ball([-displacement, 0, 0], radius)
+u = pygalmesh.Difference(s0, s1)
+```
+To sharpen the intersection circle, add it as a feature edge polygon line,
+e.g.,
+```python
+a = numpy.sqrt(radius**2 - displacement**2)
+edge_size = 0.15
+n = int(2*numpy.pi*a / edge_size)
+circ = [
+ [
+ 0.0,
+ a * numpy.cos(i * 2*numpy.pi / n),
+ a * numpy.sin(i * 2*numpy.pi / n)
+ ] for i in range(n)
+ ]
+circ.append(circ[0])
+
+pygalmesh.generate_mesh(
+ u,
+ "out.mesh",
+ feature_edges=[circ],
+ cell_size=0.15,
+ edge_size=edge_size,
+ facet_angle=25,
+ facet_size=0.15,
+ cell_radius_edge_ratio=2.0
+)
+```
+Note that the length of the polygon legs are kept in sync with the `edge_size`
+of the mesh generation. This makes sure that it fits in nicely with the rest of
+the mesh.
+
+#### Domain deformations
+<img src="https://nschloe.github.io/pygalmesh/egg.png" width="30%">
+
+You can of course translate, rotate, scale, and stretch any domain. Try, for
+example,
+```python
+import pygalmesh
+
+s = pygalmesh.Stretch(
+ pygalmesh.Ball([0, 0, 0], 1.0),
+ [1.0, 2.0, 0.0]
+)
+
+pygalmesh.generate_mesh(s, "out.mesh", cell_size=0.1)
+```
+
+#### Extrusion of 2D polygons
+<img src="https://nschloe.github.io/pygalmesh/triangle-rotated.png" width="30%">
+
+pygalmesh lets you extrude any polygon into a 3D body. It even supports
+rotation alongside!
+```python
+import pygalmesh
+
+p = pygalmesh.Polygon2D([[-0.5, -0.3], [0.5, -0.3], [0.0, 0.5]])
+edge_size = 0.1
+domain = pygalmesh.Extrude(
+ p,
+ [0.0, 0.0, 1.0],
+ 0.5 * 3.14159265359,
+ edge_size
+)
+pygalmesh.generate_mesh(
+ domain,
+ "out.mesh",
+ cell_size=0.1,
+ edge_size=edge_size,
+ verbose=False
+)
+```
+Feature edges are automatically preserved here, which is why an edge length
+needs to be given to `pygalmesh.Extrude`.
+
+#### Rotation bodies
+<img src="https://nschloe.github.io/pygalmesh/circle-rotate-extr.png" width="30%">
+
+Polygons in the x-z-plane can also be rotated around the z-axis to yield a
+rotation body.
+```python
+import pygalmesh
+
+p = pygalmesh.Polygon2D([[0.5, -0.3], [1.5, -0.3], [1.0, 0.5]])
+edge_size = 0.1
+domain = pygalmesh.RingExtrude(p, edge_size)
+pygalmesh.generate_mesh(
+ domain,
+ "out.mesh",
+ cell_size=0.1,
+ edge_size=edge_size,
+ verbose=False
+)
+```
+
+#### Your own custom level set function
+<img src="https://nschloe.github.io/pygalmesh/heart.png" width="30%">
+
+If all of the variety is not enough for you, you can define your own custom
+level set function. You simply need to subclass `pygalmesh.DomainBase` and
+specify a function, e.g.,
+```python
+import pygalmesh
+class Heart(pygalmesh.DomainBase):
+ def __init__(self):
+ super(Heart, self).__init__()
+ return
+
+ def eval(self, x):
+ return (x[0]**2 + 9.0/4.0 * x[1]**2 + x[2]**2 - 1)**3 \
+ - x[0]**2 * x[2]**3 - 9.0/80.0 * x[1]**2 * x[2]**3
+
+ def get_bounding_sphere_squared_radius(self):
+ return 10.0
+
+d = Heart()
+pygalmesh.generate_mesh(d, "out.mesh", cell_size=0.1)
+```
+Note that you need to specify the square of a bounding sphere radius, used as
+an input to CGAL's mesh generator.
+
+#### Surface meshes
+
+If you're only after the surface of a body, pygalmesh has
+`generate_surface_mesh` for you. It offers fewer options (obviously,
+`cell_size` is gone), but otherwise works the same way:
+```python
+import pygalmesh
+
+s = pygalmesh.Ball([0, 0, 0], 1.0)
+pygalmesh.generate_surface_mesh(
+ s,
+ "out.off",
+ angle_bound=30,
+ radius_bound=0.1,
+ distance_bound=0.1
+)
+```
+The output format is
+[OFF](http://segeval.cs.princeton.edu/public/off_format.html) which again is
+handled by [meshio](https://github.com/nschloe/meshio).
+
+Refer to [CGAL's
+documention](https://doc.cgal.org/latest/Surface_mesher/index.html) for the
+options.
+
+#### Meshes from OFF files
+<img src="https://nschloe.github.io/pygalmesh/elephant.png" width="30%">
+
+If you have an OFF file at hand (like
+[elephant.off](https://raw.githubusercontent.com/CGAL/cgal-swig-bindings/master/examples/data/elephant.off)
+or [these](https://github.com/CGAL/cgal/tree/master/Surface_mesher/demo/Surface_mesher/inputs)),
+pygalmesh generates the mesh via
+```python
+import pygalmesh
+
+pygalmesh.generate_from_off(
+ "elephant.off",
+ "out.mesh",
+ facet_angle=25.0,
+ facet_size=0.15,
+ facet_distance=0.008,
+ cell_radius_edge_ratio=3.0,
+ verbose=False
+)
+```
+
+### Installation
+
+For installation, pygalmesh needs [CGAL](https://www.cgal.org/) and
+[Eigen](http://eigen.tuxfamily.org/index.php?title=Main_Page) installed on your
+system. They are typically available on your Linux distribution, e.g., on
+Ubuntu
+```
+sudo apt install libcgal-dev libeigen3-dev
+```
+After that, pygalmesh can be [installed from the Python Package
+Index](https://pypi.org/project/pygalmesh/), so with
+```
+pip install -U pygalmesh
+```
+you can install/upgrade.
+
+[meshio](https://github.com/nschloe/meshio) (`sudo -H pip install meshio`)
+can be helpful in processing the meshes.
+
+#### Manual installation
+
+For manual installation (if you're a developer or just really keen on getting
+the bleeding edge version of pygalmesh), there are two possibilities:
+
+ * Get the sources, type `sudo python setup.py install`. This does the trick
+ most the time.
+ * As a fallback, there's a CMake-based installation. Simply go `cmake
+ /path/to/sources/` and `make`.
+
+### Testing
+
+To run the pygalmesh unit tests, check out this repository and type
+```
+pytest
+```
+
+### Distribution
+
+To create a new release
+
+1. bump the `__version__` number (in `setup.py` _and_ `src/pygalmesh.i`)
+
+2. publish to PyPi and GitHub:
+ ```
+ make publish
+ ```
+
+### License
+
+pygalmesh is published under the [MIT license](https://en.wikipedia.org/wiki/MIT_License).
--- /dev/null
+# -*- coding: utf-8 -*-
+#
+
+__author__ = u"Nico Schlömer"
+__author_email__ = "nico.schloemer@gmail.com"
+__copyright__ = u"Copyright (c) 2016-2018, {} <{}>".format(__author__, __author_email__)
+__license__ = "License :: OSI Approved :: MIT License"
+__version__ = "0.2.6"
+__maintainer__ = u"Nico Schlömer"
+__status__ = "Development Status :: 4 - Beta"
+__url__ = "https://github.com/nschloe/pygalmesh"
--- /dev/null
+# -*- coding: utf-8 -*-
+#
+# https://github.com/pybind/pybind11/issues/1004
+from _pygalmesh import (
+ DomainBase,
+ Translate,
+ Rotate,
+ Scale,
+ Stretch,
+ Intersection,
+ Union,
+ Difference,
+ Extrude,
+ Ball,
+ Cuboid,
+ Ellipsoid,
+ Tetrahedron,
+ Cone,
+ Cylinder,
+ Torus,
+ HalfSpace,
+ Polygon2D,
+ RingExtrude,
+ generate_mesh,
+ generate_surface_mesh,
+ generate_from_off,
+)
+
+from .__about__ import (
+ __author__,
+ __author_email__,
+ __copyright__,
+ __license__,
+ __version__,
+ __maintainer__,
+ __status__,
+)
+
+__all__ = [
+ "__author__",
+ "__author_email__",
+ "__copyright__",
+ "__license__",
+ "__version__",
+ "__maintainer__",
+ "__status__",
+ #
+ "DomainBase",
+ "Translate",
+ "Rotate",
+ "Scale",
+ "Stretch",
+ "Intersection",
+ "Union",
+ "Difference",
+ "Extrude",
+ "Ball",
+ "Cuboid",
+ "Ellipsoid",
+ "Tetrahedron",
+ "Cone",
+ "Cylinder",
+ "Torus",
+ "HalfSpace",
+ "Polygon2D",
+ "RingExtrude",
+ "generate_mesh",
+ "generate_surface_mesh",
+ "generate_from_off",
+]
--- /dev/null
+# -*- coding: utf-8 -*-
+#
+import os
+
+import codecs
+from setuptools import setup, Extension, find_packages
+
+
+# https://packaging.python.org/single_source_version/
+base_dir = os.path.abspath(os.path.dirname(__file__))
+about = {}
+with open(os.path.join(base_dir, "pygalmesh", "__about__.py"), "rb") as handle:
+ # pylint: disable=exec-used
+ exec(handle.read(), about)
+
+
+class get_pybind_include(object):
+ """Helper class to determine the pybind11 include path
+ The purpose of this class is to postpone importing pybind11
+ until it is actually installed, so that the ``get_include()``
+ method can be invoked.
+ """
+
+ def __init__(self, user=False):
+ self.user = user
+
+ def __str__(self):
+ import pybind11
+
+ return pybind11.get_include(self.user)
+
+
+def read(fname):
+ return codecs.open(os.path.join(base_dir, fname), encoding="utf-8").read()
+
+
+ext_modules = [
+ Extension(
+ "_pygalmesh",
+ [
+ "src/generate.cpp",
+ "src/generate_from_off.cpp",
+ "src/generate_surface_mesh.cpp",
+ "src/pybind11.cpp",
+ ],
+ language="c++",
+ include_dirs=[
+ "/usr/include/eigen3/",
+ # Path to pybind11 headers
+ get_pybind_include(),
+ get_pybind_include(user=True),
+ ],
+ libraries=["CGAL", "gmp", "mpfr"],
+ # extra_compile_args=['-std=c++11']
+ )
+]
+
+setup(
+ name="pygalmesh",
+ packages=find_packages(),
+ # cmdclass={'build_ext': BuildExt},
+ ext_modules=ext_modules,
+ version=about["__version__"],
+ url=about["__url__"],
+ author=about["__author__"],
+ author_email=about["__author_email__"],
+ install_requires=["numpy", "pybind11 >= 2.2", "pipdate"],
+ description="Python frontend to CGAL's 3D mesh generation capabilities",
+ long_description=read("README.md"),
+ long_description_content_type="text/markdown",
+ license=about["__license__"],
+ classifiers=[
+ about["__status__"],
+ about["__license__"],
+ "Operating System :: OS Independent",
+ "Programming Language :: Python",
+ "Programming Language :: Python :: 2",
+ "Programming Language :: Python :: 3",
+ "Topic :: Scientific/Engineering",
+ "Topic :: Scientific/Engineering :: Mathematics",
+ "Topic :: Scientific/Engineering :: Physics",
+ "Topic :: Scientific/Engineering :: Visualization",
+ ],
+)
--- /dev/null
+FIND_PACKAGE(pybind11 REQUIRED)
+# include_directories(${PYBIND11_INCLUDE_DIR})
+
+FIND_PACKAGE(Eigen3 REQUIRED)
+include_directories(${EIGEN3_INCLUDE_DIR})
+
+FILE(GLOB pygalmesh_SRCS "${CMAKE_CURRENT_SOURCE_DIR}/*.cpp")
+# FILE(GLOB pygalmesh_HEADERS "${CMAKE_CURRENT_SOURCE_DIR}/*.hpp")
+
+pybind11_add_module(pygalmesh ${pygalmesh_SRCS})
+
+# Call CGAL as late as possible. It messes around with some
+# compiler variables. See bugs
+# <https://github.com/CGAL/cgal/issues/2849> and
+# <https://github.com/CGAL/cgal/issues/2850>.
+FIND_PACKAGE(CGAL REQUIRED)
+include(${CGAL_USE_FILE})
+
+# ADD_LIBRARY(pygalmesh ${pygalmesh_SRCS})
+target_link_libraries(pygalmesh PRIVATE ${CGAL_LIBRARIES})
+
+# execute_process(
+# COMMAND python -c "from distutils.sysconfig import get_python_lib; print get_python_lib()"
+# OUTPUT_VARIABLE PYTHON_SITE_PACKAGES
+# OUTPUT_STRIP_TRAILING_WHITESPACE
+# )
+# install(TARGETS _pygalmesh DESTINATION ${PYTHON_SITE_PACKAGES})
+# install(
+# FILES ${CMAKE_BINARY_DIR}/src/pygalmesh.py
+# DESTINATION ${PYTHON_SITE_PACKAGES}
+# )
--- /dev/null
+#ifndef DOMAIN_HPP
+#define DOMAIN_HPP
+
+#include <Eigen/Dense>
+#include <array>
+#include <limits>
+#include <memory>
+#include <vector>
+
+namespace pygalmesh {
+
+class DomainBase
+{
+ public:
+
+ virtual ~DomainBase() = default;
+
+ virtual
+ double
+ eval(const std::array<double, 3> & x) const = 0;
+
+ virtual
+ double
+ get_bounding_sphere_squared_radius() const = 0;
+
+ virtual
+ std::vector<std::vector<std::array<double, 3>>>
+ get_features() const
+ {
+ return {};
+ };
+};
+
+class Translate: public pygalmesh::DomainBase
+{
+ public:
+ Translate(
+ const std::shared_ptr<const pygalmesh::DomainBase> & domain,
+ const std::array<double, 3> & direction
+ ):
+ domain_(domain),
+ direction_(Eigen::Vector3d(direction.data())),
+ translated_features_(translate_features(domain->get_features(), direction_))
+ {
+ }
+
+ virtual ~Translate() = default;
+
+ std::vector<std::vector<std::array<double, 3>>>
+ translate_features(
+ const std::vector<std::vector<std::array<double, 3>>> & features,
+ const Eigen::Vector3d & direction
+ ) const
+ {
+ std::vector<std::vector<std::array<double, 3>>> translated_features;
+ for (const auto & feature: features) {
+ std::vector<std::array<double, 3>> translated_feature;
+ for (const auto & point: feature) {
+ const std::array<double, 3> translated_point = {
+ point[0] + direction[0],
+ point[1] + direction[1],
+ point[2] + direction[2]
+ };
+ translated_feature.push_back(translated_point);
+ }
+ translated_features.push_back(translated_feature);
+ }
+ return translated_features;
+ }
+
+ virtual
+ double
+ eval(const std::array<double, 3> & x) const
+ {
+ const std::array<double, 3> d = {
+ x[0] - direction_[0],
+ x[1] - direction_[1],
+ x[2] - direction_[2]
+ };
+ return domain_->eval(d);
+ }
+
+ virtual
+ double
+ get_bounding_sphere_squared_radius() const
+ {
+ const double radius = sqrt(domain_->get_bounding_sphere_squared_radius());
+ const double dir_norm = direction_.norm();
+ return (radius + dir_norm)*(radius + dir_norm);
+ }
+
+ virtual
+ std::vector<std::vector<std::array<double, 3>>>
+ get_features() const
+ {
+ return translated_features_;
+ };
+
+ private:
+ const std::shared_ptr<const pygalmesh::DomainBase> domain_;
+ const Eigen::Vector3d direction_;
+ const std::vector<std::vector<std::array<double, 3>>> translated_features_;
+};
+
+class Rotate: public pygalmesh::DomainBase
+{
+ public:
+ Rotate(
+ const std::shared_ptr<const pygalmesh::DomainBase> & domain,
+ const std::array<double, 3> & axis,
+ const double angle
+ ):
+ domain_(domain),
+ normalized_axis_(Eigen::Vector3d(axis.data()).normalized()),
+ sinAngle_(sin(angle)),
+ cosAngle_(cos(angle)),
+ rotated_features_(rotate_features(domain_->get_features()))
+ {
+ }
+
+ virtual ~Rotate() = default;
+
+ Eigen::Vector3d
+ rotate(
+ const Eigen::Vector3d & vec,
+ const Eigen::Vector3d & axis,
+ const double sinAngle,
+ const double cosAngle
+ ) const
+ {
+ // Rotate a vector `v` by the angle `theta` in the plane perpendicular
+ // to the axis given by `u`.
+ // Refer to
+ // http://en.wikipedia.org/wiki/Rotation_matrix#Rotation_matrix_from_axis_and_angle
+ //
+ // cos(theta) * I * v
+ // + sin(theta) u\cross v
+ // + (1-cos(theta)) (u*u^T) * v
+ return cosAngle * vec
+ + sinAngle * axis.cross(vec)
+ + (1.0-cosAngle) * axis.dot(vec) * axis;
+ }
+
+ virtual
+ double
+ eval(const std::array<double, 3> & x) const
+ {
+ // rotate with negative angle
+ const auto p2 = rotate(
+ Eigen::Vector3d(x.data()),
+ normalized_axis_,
+ -sinAngle_,
+ cosAngle_
+ );
+ return domain_->eval({p2[0], p2[1], p2[2]});
+ }
+
+ std::vector<std::vector<std::array<double, 3>>>
+ rotate_features(
+ const std::vector<std::vector<std::array<double, 3>>> & features
+ ) const
+ {
+ std::vector<std::vector<std::array<double, 3>>> rotated_features;
+ for (const auto & feature: features) {
+ std::vector<std::array<double, 3>> rotated_feature;
+ for (const auto & point: feature) {
+ const auto p2 = rotate(
+ Eigen::Vector3d(point.data()),
+ normalized_axis_,
+ sinAngle_,
+ cosAngle_
+ );
+ rotated_feature.push_back({p2[0], p2[1], p2[2]});
+ }
+ rotated_features.push_back(rotated_feature);
+ }
+ return rotated_features;
+ }
+
+ virtual
+ double
+ get_bounding_sphere_squared_radius() const
+ {
+ return domain_->get_bounding_sphere_squared_radius();
+ }
+
+ virtual
+ std::vector<std::vector<std::array<double, 3>>>
+ get_features() const
+ {
+ return rotated_features_;
+ };
+
+ private:
+ const std::shared_ptr<const pygalmesh::DomainBase> domain_;
+ const Eigen::Vector3d normalized_axis_;
+ const double sinAngle_;
+ const double cosAngle_;
+ const std::vector<std::vector<std::array<double, 3>>> rotated_features_;
+};
+
+class Scale: public pygalmesh::DomainBase
+{
+ public:
+ Scale(
+ std::shared_ptr<const pygalmesh::DomainBase> & domain,
+ const double alpha
+ ):
+ domain_(domain),
+ alpha_(alpha),
+ scaled_features_(scale_features(domain_->get_features()))
+ {
+ assert(alpha_ > 0.0);
+ }
+
+ virtual ~Scale() = default;
+
+ virtual
+ double
+ eval(const std::array<double, 3> & x) const
+ {
+ return domain_->eval({x[0]/alpha_, x[1]/alpha_, x[2]/alpha_});
+ }
+
+ virtual
+ double
+ get_bounding_sphere_squared_radius() const
+ {
+ return alpha_*alpha_ * domain_->get_bounding_sphere_squared_radius();
+ }
+
+ std::vector<std::vector<std::array<double, 3>>>
+ scale_features(
+ const std::vector<std::vector<std::array<double, 3>>> & features
+ ) const
+ {
+ std::vector<std::vector<std::array<double, 3>>> scaled_features;
+ for (const auto & feature: features) {
+ std::vector<std::array<double, 3>> scaled_feature;
+ for (const auto & point: feature) {
+ scaled_feature.push_back({
+ alpha_ * point[0],
+ alpha_ * point[1],
+ alpha_ * point[2]
+ });
+ }
+ scaled_features.push_back(scaled_feature);
+ }
+ return scaled_features;
+ }
+
+ virtual
+ std::vector<std::vector<std::array<double, 3>>>
+ get_features() const
+ {
+ return scaled_features_;
+ };
+
+ private:
+ std::shared_ptr<const pygalmesh::DomainBase> domain_;
+ const double alpha_;
+ const std::vector<std::vector<std::array<double, 3>>> scaled_features_;
+};
+
+class Stretch: public pygalmesh::DomainBase
+{
+ public:
+ Stretch(
+ std::shared_ptr<const pygalmesh::DomainBase> & domain,
+ const std::array<double, 3> & direction
+ ):
+ domain_(domain),
+ normalized_direction_(Eigen::Vector3d(direction.data()).normalized()),
+ alpha_(Eigen::Vector3d(direction.data()).norm()),
+ stretched_features_(stretch_features(domain_->get_features()))
+ {
+ assert(alpha_ > 0.0);
+ }
+
+ virtual ~Stretch() = default;
+
+ virtual
+ double
+ eval(const std::array<double, 3> & x) const
+ {
+ const Eigen::Vector3d v(x.data());
+ const double beta = normalized_direction_.dot(v);
+ // scale the component of normalized_direction_ by 1/alpha_
+ const auto v2 = beta/alpha_ * normalized_direction_
+ + (v - beta * normalized_direction_);
+ return domain_->eval({v2[0], v2[1], v2[2]});
+ }
+
+ virtual
+ double
+ get_bounding_sphere_squared_radius() const
+ {
+ return alpha_*alpha_ * domain_->get_bounding_sphere_squared_radius();
+ }
+
+ std::vector<std::vector<std::array<double, 3>>>
+ stretch_features(
+ const std::vector<std::vector<std::array<double, 3>>> & features
+ ) const
+ {
+ std::vector<std::vector<std::array<double, 3>>> stretched_features;
+ for (const auto & feature: features) {
+ std::vector<std::array<double, 3>> stretched_feature;
+ for (const auto & point: feature) {
+ // scale the component of normalized_direction_ by alpha_
+ const Eigen::Vector3d v(point.data());
+ const double beta = normalized_direction_.dot(v);
+ const auto v2 = beta * alpha_ * normalized_direction_
+ + (v - beta * normalized_direction_);
+ stretched_feature.push_back({v2[0], v2[1], v2[2]});
+ }
+ stretched_features.push_back(stretched_feature);
+ }
+ return stretched_features;
+ }
+
+ virtual
+ std::vector<std::vector<std::array<double, 3>>>
+ get_features() const
+ {
+ return stretched_features_;
+ };
+
+ private:
+ std::shared_ptr<const pygalmesh::DomainBase> domain_;
+ const Eigen::Vector3d normalized_direction_;
+ const double alpha_;
+ const std::vector<std::vector<std::array<double, 3>>> stretched_features_;
+};
+
+class Intersection: public pygalmesh::DomainBase
+{
+ public:
+ explicit Intersection(
+ std::vector<std::shared_ptr<const pygalmesh::DomainBase>> & domains
+ ):
+ domains_(domains)
+ {
+ }
+
+ virtual ~Intersection() = default;
+
+ virtual
+ double
+ eval(const std::array<double, 3> & x) const
+ {
+ // TODO find a differentiable expression
+ double maxval = std::numeric_limits<double>::lowest();
+ for (const auto & domain: domains_) {
+ maxval = std::max(maxval, domain->eval(x));
+ }
+ return maxval;
+ }
+
+ virtual
+ double
+ get_bounding_sphere_squared_radius() const
+ {
+ double min = std::numeric_limits<double>::max();
+ for (const auto & domain: domains_) {
+ min = std::min(min, domain->get_bounding_sphere_squared_radius());
+ }
+ return min;
+ }
+
+ virtual
+ std::vector<std::vector<std::array<double, 3>>>
+ get_features() const
+ {
+ std::vector<std::vector<std::array<double, 3>>> features;
+ for (const auto & domain: domains_) {
+ const auto f = domain->get_features();
+ features.insert(std::end(features), std::begin(f), std::end(f));
+ }
+ return features;
+ };
+
+ private:
+ std::vector<std::shared_ptr<const pygalmesh::DomainBase>> domains_;
+};
+
+class Union: public pygalmesh::DomainBase
+{
+ public:
+ explicit Union(
+ std::vector<std::shared_ptr<const pygalmesh::DomainBase>> & domains
+ ):
+ domains_(domains)
+ {
+ }
+
+ virtual ~Union() = default;
+
+ virtual
+ double
+ eval(const std::array<double, 3> & x) const
+ {
+ // TODO find a differentiable expression
+ double minval = std::numeric_limits<double>::max();
+ for (const auto & domain: domains_) {
+ minval = std::min(minval, domain->eval(x));
+ }
+ return minval;
+ }
+
+ virtual
+ double
+ get_bounding_sphere_squared_radius() const
+ {
+ double max = 0.0;
+ for (const auto & domain: domains_) {
+ max = std::max(max, domain->get_bounding_sphere_squared_radius());
+ }
+ return max;
+ }
+
+ virtual
+ std::vector<std::vector<std::array<double, 3>>>
+ get_features() const
+ {
+ std::vector<std::vector<std::array<double, 3>>> features;
+ for (const auto & domain: domains_) {
+ const auto f = domain->get_features();
+ features.insert(std::end(features), std::begin(f), std::end(f));
+ }
+ return features;
+ };
+
+ private:
+ std::vector<std::shared_ptr<const pygalmesh::DomainBase>> domains_;
+};
+
+class Difference: public pygalmesh::DomainBase
+{
+ public:
+ Difference(
+ std::shared_ptr<const pygalmesh::DomainBase> & domain0,
+ std::shared_ptr<const pygalmesh::DomainBase> & domain1
+ ):
+ domain0_(domain0),
+ domain1_(domain1)
+ {
+ }
+
+ virtual ~Difference() = default;
+
+ virtual
+ double
+ eval(const std::array<double, 3> & x) const
+ {
+ // TODO find a continuous (perhaps even differentiable) expression
+ const double val0 = domain0_->eval(x);
+ const double val1 = domain1_->eval(x);
+ return (val0 < 0.0 && val1 >= 0.0) ? val0 : std::max(val0, -val1);
+ }
+
+ virtual
+ double
+ get_bounding_sphere_squared_radius() const
+ {
+ return domain0_->get_bounding_sphere_squared_radius();
+ }
+
+ virtual
+ std::vector<std::vector<std::array<double, 3>>>
+ get_features() const
+ {
+ std::vector<std::vector<std::array<double, 3>>> features;
+
+ const auto f0 = domain0_->get_features();
+ features.insert(std::end(features), std::begin(f0), std::end(f0));
+
+ const auto f1 = domain1_->get_features();
+ features.insert(std::end(features), std::begin(f1), std::end(f1));
+
+ return features;
+ };
+
+ private:
+ std::shared_ptr<const pygalmesh::DomainBase> domain0_;
+ std::shared_ptr<const pygalmesh::DomainBase> domain1_;
+};
+
+} // namespace pygalmesh
+#endif // DOMAIN_HPP
--- /dev/null
+#define CGAL_MESH_3_VERBOSE 1
+
+#include "generate.hpp"
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <CGAL/Mesh_triangulation_3.h>
+#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
+#include <CGAL/Mesh_criteria_3.h>
+
+#include <CGAL/Implicit_mesh_domain_3.h>
+#include <CGAL/Mesh_domain_with_polyline_features_3.h>
+#include <CGAL/make_mesh_3.h>
+
+namespace pygalmesh {
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+
+// Wrapper for DomainBase for translating to K::Point_3/FT.
+class CgalDomainWrapper
+{
+ public:
+ explicit CgalDomainWrapper(const std::shared_ptr<DomainBase> & domain):
+ domain_(domain)
+ {
+ }
+
+ virtual ~CgalDomainWrapper() = default;
+
+ virtual
+ K::FT
+ operator()(K::Point_3 p) const
+ {
+ return domain_->eval({p.x(), p.y(), p.z()});
+ }
+
+ private:
+ const std::shared_ptr<DomainBase> domain_;
+};
+
+typedef CGAL::Mesh_domain_with_polyline_features_3<
+ CGAL::Implicit_mesh_domain_3<CgalDomainWrapper,K>
+ >
+ Mesh_domain;
+
+// Triangulation
+typedef CGAL::Mesh_triangulation_3<Mesh_domain>::type Tr;
+typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
+
+// Mesh Criteria
+typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
+typedef Mesh_criteria::Facet_criteria Facet_criteria;
+typedef Mesh_criteria::Cell_criteria Cell_criteria;
+
+// translate vector<vector<array<double, 3>> to list<vector<Point_3>>
+std::list<std::vector<K::Point_3>>
+translate_feature_edges(
+ const std::vector<std::vector<std::array<double, 3>>> & feature_edges
+ )
+{
+ std::list<std::vector<K::Point_3>> polylines;
+ for (const auto & feature_edge: feature_edges) {
+ std::vector<K::Point_3> polyline;
+ for (const auto & point: feature_edge) {
+ polyline.push_back(K::Point_3(point[0], point[1], point[2]));
+ }
+ polylines.push_back(polyline);
+ }
+ return polylines;
+}
+
+void
+generate_mesh(
+ 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 double boundary_precision,
+ 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 double cell_size,
+ const bool verbose
+ )
+{
+ 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();
+
+ const auto d = CgalDomainWrapper(domain);
+ Mesh_domain cgal_domain(
+ d,
+ K::Sphere_3(CGAL::ORIGIN, bounding_sphere_radius2),
+ boundary_precision
+ );
+
+ 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
+ if (!verbose) {
+ // suppress output
+ std::cerr.setstate(std::ios_base::failbit);
+ }
+ 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
--- /dev/null
+#ifndef GENERATE_HPP
+#define GENERATE_HPP
+
+#include "domain.hpp"
+
+#include <memory>
+#include <string>
+#include <vector>
+
+namespace pygalmesh {
+
+void generate_mesh(
+ 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 = 0.0,
+ const double boundary_precision = 1.0e-4,
+ const bool lloyd = false,
+ const bool odt = false,
+ const bool perturb = true,
+ const bool exude = true,
+ const double edge_size = 0.0, // std::numeric_limits<double>::max(),
+ const double facet_angle = 0.0,
+ const double facet_size = 0.0,
+ const double facet_distance = 0.0,
+ const double cell_radius_edge_ratio = 0.0,
+ const double cell_size = 0.0,
+ const bool verbose = true
+ );
+
+} // namespace pygalmesh
+
+#endif // GENERATE_HPP
--- /dev/null
+#include "generate_from_off.hpp"
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Mesh_triangulation_3.h>
+#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
+#include <CGAL/Mesh_criteria_3.h>
+#include <CGAL/Polyhedral_mesh_domain_3.h>
+#include <CGAL/make_mesh_3.h>
+#include <CGAL/refine_mesh_3.h>
+
+// IO
+#include <CGAL/IO/Polyhedron_iostream.h>
+
+namespace pygalmesh {
+
+// Domain
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+typedef CGAL::Polyhedron_3<K> Polyhedron;
+typedef CGAL::Polyhedral_mesh_domain_3<Polyhedron, K> Mesh_domain;
+
+// Triangulation
+typedef CGAL::Mesh_triangulation_3<Mesh_domain>::type Tr;
+
+typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
+
+// Criteria
+typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
+
+// To avoid verbose function and named parameters call
+using namespace CGAL::parameters;
+
+void
+generate_from_off(
+ const std::string & infile,
+ const std::string & outfile,
+ 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 double cell_size,
+ const bool verbose
+ )
+{
+ // Create input polyhedron
+ Polyhedron polyhedron;
+ std::ifstream input(infile);
+ input >> polyhedron;
+ if (!input.good()) {
+ std::stringstream msg;
+ msg << "Cannot read input file \"" << infile << "\"" << std::endl;
+ throw std::runtime_error(msg.str());
+ }
+ input.close();
+
+ // Create domain
+ Mesh_domain cgal_domain(polyhedron);
+
+ // Mesh criteria
+ 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
+ if (!verbose) {
+ // suppress output
+ std::cerr.setstate(std::ios_base::failbit);
+ }
+ // TODO make_mesh_3 segfaults on travis. hm.
+ 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
--- /dev/null
+#ifndef GENERATE_FROM_OFF_HPP
+#define GENERATE_FROM_OFF_HPP
+
+#include <string>
+#include <vector>
+
+namespace pygalmesh {
+
+void
+generate_from_off(
+ const std::string & infile,
+ const std::string & outfile,
+ const bool lloyd = false,
+ const bool odt = false,
+ const bool perturb = true,
+ const bool exude = true,
+ const double edge_size = 0.0, // std::numeric_limits<double>::max(),
+ const double facet_angle = 0.0,
+ const double facet_size = 0.0,
+ const double facet_distance = 0.0,
+ const double cell_radius_edge_ratio = 0.0,
+ const double cell_size = 0.0,
+ const bool verbose = true
+ );
+
+} // namespace pygalmesh
+
+#endif // GENERATE_FROM_OFF_HPP
--- /dev/null
+#define CGAL_SURFACE_MESHER_VERBOSE 1
+
+#include "generate_surface_mesh.hpp"
+
+#include <CGAL/Surface_mesh_default_triangulation_3.h>
+#include <CGAL/Complex_2_in_triangulation_3.h>
+#include <CGAL/make_surface_mesh.h>
+#include <fstream>
+#include <CGAL/IO/Complex_2_in_triangulation_3_file_writer.h>
+#include <CGAL/Implicit_surface_3.h>
+
+namespace pygalmesh {
+
+// default triangulation for Surface_mesher
+typedef CGAL::Surface_mesh_default_triangulation_3 Tr;
+// c2t3
+typedef CGAL::Complex_2_in_triangulation_3<Tr> C2t3;
+typedef Tr::Geom_traits GT;
+
+// Wrapper for DomainBase for translating to GT.
+class CgalDomainWrapper
+{
+ public:
+ explicit CgalDomainWrapper(const std::shared_ptr<DomainBase> & domain):
+ domain_(domain)
+ {
+ }
+
+ virtual ~CgalDomainWrapper() = default;
+
+ virtual
+ GT::FT
+ operator()(GT::Point_3 p) const
+ {
+ return domain_->eval({p.x(), p.y(), p.z()});
+ }
+
+ private:
+ const std::shared_ptr<DomainBase> domain_;
+};
+
+typedef CGAL::Implicit_surface_3<GT, CgalDomainWrapper> Surface_3;
+
+void
+generate_surface_mesh(
+ const std::shared_ptr<pygalmesh::DomainBase> & domain,
+ const std::string & outfile,
+ const double bounding_sphere_radius,
+ const double angle_bound,
+ const double radius_bound,
+ const double distance_bound,
+ const bool verbose
+ )
+{
+ const double bounding_sphere_radius2 = bounding_sphere_radius > 0 ?
+ bounding_sphere_radius*bounding_sphere_radius :
+ // add a little wiggle room
+ 1.01 * domain->get_bounding_sphere_squared_radius();
+
+ Tr tr; // 3D-Delaunay triangulation
+ C2t3 c2t3 (tr); // 2D-complex in 3D-Delaunay triangulation
+
+ const auto d = CgalDomainWrapper(domain);
+ Surface_3 surface(
+ d,
+ GT::Sphere_3(CGAL::ORIGIN, bounding_sphere_radius2)
+ );
+
+ CGAL::Surface_mesh_default_criteria_3<Tr> criteria(
+ angle_bound,
+ radius_bound,
+ distance_bound
+ );
+
+ if (!verbose) {
+ // suppress output
+ std::cout.setstate(std::ios_base::failbit);
+ std::cerr.setstate(std::ios_base::failbit);
+ }
+ CGAL::make_surface_mesh(
+ c2t3,
+ surface,
+ criteria,
+ CGAL::Non_manifold_tag()
+ );
+ if (!verbose) {
+ std::cout.clear();
+ std::cerr.clear();
+ }
+
+ // Output
+ std::ofstream off_file(outfile);
+ CGAL::output_surface_facets_to_off(off_file, c2t3);
+ off_file.close();
+
+ return;
+}
+
+} // namespace pygalmesh
--- /dev/null
+#ifndef GENERATE_SURFACE_MESH_HPP
+#define GENERATE_SURFACE_MESH_HPP
+
+#include "domain.hpp"
+
+#include <memory>
+#include <string>
+
+namespace pygalmesh {
+
+void generate_surface_mesh(
+ const std::shared_ptr<pygalmesh::DomainBase> & domain,
+ const std::string & outfile,
+ const double bounding_sphere_radius = 0.0,
+ const double angle_bound = 0.0,
+ const double radius_bound = 0.0,
+ const double distance_bound = 0.0,
+ const bool verbose = true
+ );
+
+} // namespace pygalmesh
+
+#endif // GENERATE_SURFACE_MESH_HPP
--- /dev/null
+#ifndef POLYGON2D_HPP
+#define POLYGON2D_HPP
+
+#include "domain.hpp"
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Polygon_2_algorithms.h>
+#include <array>
+#include <memory>
+#include <vector>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+
+namespace pygalmesh {
+
+class Polygon2D {
+ public:
+ explicit Polygon2D(const std::vector<std::array<double, 2>> & _points):
+ points(vector_to_cgal_points(_points))
+ {
+ }
+
+ virtual ~Polygon2D() = default;
+
+ std::vector<K::Point_2>
+ vector_to_cgal_points(const std::vector<std::array<double, 2>> & _points) const
+ {
+ std::vector<K::Point_2> points2(_points.size());
+ for (size_t i = 0; i < _points.size(); i++) {
+ assert(_points[i].size() == 2);
+ points2[i] = K::Point_2(_points[i][0], _points[i][1]);
+ }
+ return points2;
+ }
+
+ bool
+ is_inside(const std::array<double, 2> & point)
+ {
+ K::Point_2 pt(point[0], point[1]);
+ switch(CGAL::bounded_side_2(this->points.begin(), this->points.end(), pt, K())) {
+ case CGAL::ON_BOUNDED_SIDE:
+ return true;
+ case CGAL::ON_BOUNDARY:
+ return true;
+ case CGAL::ON_UNBOUNDED_SIDE:
+ return false;
+ default:
+ return false;
+ }
+ return false;
+ }
+
+ public:
+ const std::vector<K::Point_2> points;
+};
+
+
+class Extrude: public pygalmesh::DomainBase {
+ public:
+ Extrude(
+ const std::shared_ptr<pygalmesh::Polygon2D> & poly,
+ const std::array<double, 3> & direction,
+ const double alpha = 0.0,
+ const double edge_size = 0.0
+ ):
+ poly_(poly),
+ direction_(direction),
+ alpha_(alpha),
+ edge_size_(edge_size)
+ {
+ }
+
+ virtual ~Extrude() = default;
+
+ virtual
+ double
+ eval(const std::array<double, 3> & x) const
+ {
+ if (x[2] < 0.0 || x[2] > direction_[2]) {
+ return 1.0;
+ }
+
+ const double beta = x[2] / direction_[2];
+
+ std::array<double, 2> x2 = {
+ x[0] - beta * direction_[0],
+ x[1] - beta * direction_[1]
+ };
+
+ if (alpha_ != 0.0) {
+ std::array<double, 2> x3;
+ // turn by -beta*alpha
+ const double sinAlpha = sin(beta*alpha_);
+ const double cosAlpha = cos(beta*alpha_);
+ x3[0] = cosAlpha * x2[0] + sinAlpha * x2[1];
+ x3[1] = -sinAlpha * x2[0] + cosAlpha * x2[1];
+ x2 = x3;
+ }
+
+ return poly_->is_inside(x2) ? -1.0 : 1.0;
+ }
+
+ virtual
+ double
+ get_bounding_sphere_squared_radius() const
+ {
+ double max = 0.0;
+ for (const auto & pt: poly_->points) {
+ // bottom polygon
+ const double nrm0 = pt.x()*pt.x() + pt.y()*pt.y();
+ if (nrm0 > max) {
+ max = nrm0;
+ }
+
+ // TODO rotation
+
+ // top polygon
+ const double x = pt.x() + direction_[0];
+ const double y = pt.y() + direction_[1];
+ const double z = direction_[2];
+ const double nrm1 = x*x + y*y + z*z;
+ if (nrm1 > max) {
+ max = nrm1;
+ }
+ }
+ return max;
+ }
+
+ virtual
+ std::vector<std::vector<std::array<double, 3>>>
+ get_features() const
+ {
+ std::vector<std::vector<std::array<double, 3>>> features = {};
+
+ size_t n;
+
+ // bottom polygon
+ n = poly_->points.size();
+ for (size_t i=0; i < n-1; i++) {
+ features.push_back({
+ {poly_->points[i].x(), poly_->points[i].y(), 0.0},
+ {poly_->points[i+1].x(), poly_->points[i+1].y(), 0.0}
+ });
+ }
+ features.push_back({
+ {poly_->points[n-1].x(), poly_->points[n-1].y(), 0.0},
+ {poly_->points[0].x(), poly_->points[0].y(), 0.0}
+ });
+
+ // top polygon, R*x + d
+ n = poly_->points.size();
+ const double sinAlpha = sin(alpha_);
+ const double cosAlpha = cos(alpha_);
+ for (size_t i=0; i < n-1; i++) {
+ features.push_back({
+ {
+ cosAlpha * poly_->points[i].x() - sinAlpha * poly_->points[i].y() + direction_[0],
+ sinAlpha * poly_->points[i].x() + cosAlpha * poly_->points[i].y() + direction_[1],
+ direction_[2]
+ },
+ {
+ cosAlpha * poly_->points[i+1].x() - sinAlpha * poly_->points[i+1].y() + direction_[0],
+ sinAlpha * poly_->points[i+1].x() + cosAlpha * poly_->points[i+1].y() + direction_[1],
+ direction_[2]
+ }
+ });
+ }
+ features.push_back({
+ {
+ cosAlpha * poly_->points[n-1].x() - sinAlpha * poly_->points[n-1].y() + direction_[0],
+ sinAlpha * poly_->points[n-1].x() + cosAlpha * poly_->points[n-1].y() + direction_[1],
+ direction_[2]
+ },
+ {
+ cosAlpha * poly_->points[0].x() - sinAlpha * poly_->points[0].y() + direction_[0],
+ sinAlpha * poly_->points[0].x() + cosAlpha * poly_->points[0].y() + direction_[1],
+ direction_[2]
+ }
+ });
+
+ // features connecting the top and bottom
+ if (alpha_ == 0) {
+ for (const auto & pt: poly_->points) {
+ std::vector<std::array<double, 3>> line = {
+ {pt.x(), pt.y(), 0.0},
+ {pt.x() + direction_[0], pt.y() + direction_[1], direction_[2]}
+ };
+ features.push_back(line);
+ }
+ } else {
+ // Alright, we need to chop the lines on which the polygon corners are
+ // sitting into pieces. How long? About edge_size. For the starting point
+ // (x0, y0, z0) height h and angle alpha, the lines are given by
+ //
+ // f(beta) = (
+ // cos(alpha*beta) x0 - sin(alpha*beta) y0,
+ // sin(alpha*beta) x0 + cos(alpha*beta) y0,
+ // z0 + beta * h
+ // )
+ //
+ // with beta in [0, 1]. The length from beta0 till beta1 is then
+ //
+ // l = sqrt(alpha^2 (x0^2 + y0^2) + h^2) * (beta1 - beta0).
+ //
+ const double height = direction_[2];
+ for (const auto & pt: poly_->points) {
+ const double l = sqrt(alpha_*alpha_ * (pt.x()*pt.x() + pt.y()*pt.y()) + height*height);
+ assert(edge_size_ > 0.0);
+ const size_t n = int(l / edge_size_ - 0.5) + 1;
+ std::vector<std::array<double, 3>> line = {
+ {pt.x(), pt.y(), 0.0},
+ };
+ for (size_t i=0; i < n; i++) {
+ const double beta = double(i+1) / n;
+ const double sinAB = sin(alpha_*beta);
+ const double cosAB = cos(alpha_*beta);
+ line.push_back({
+ cosAB * pt.x() - sinAB * pt.y(),
+ sinAB * pt.x() + cosAB * pt.y(),
+ beta * height
+ });
+ }
+ features.push_back(line);
+ }
+ }
+
+ return features;
+ };
+
+ private:
+ const std::shared_ptr<pygalmesh::Polygon2D> poly_;
+ const std::array<double, 3> direction_;
+ const double alpha_;
+ const double edge_size_;
+};
+
+
+class ring_extrude: public pygalmesh::DomainBase {
+ public:
+ ring_extrude(
+ const std::shared_ptr<pygalmesh::Polygon2D> & poly,
+ const double edge_size
+ ):
+ poly_(poly),
+ edge_size_(edge_size)
+ {
+ assert(edge_size > 0.0);
+ }
+
+ virtual ~ring_extrude() = default;
+
+ virtual
+ double
+ eval(const std::array<double, 3> & x) const
+ {
+ const double r = sqrt(x[0]*x[0] + x[1]*x[1]);
+ const double z = x[2];
+
+ return poly_->is_inside({r, z}) ? -1.0 : 1.0;
+ }
+
+ virtual
+ double
+ get_bounding_sphere_squared_radius() const
+ {
+ double max = 0.0;
+ for (const auto & pt: poly_->points) {
+ const double nrm1 = pt.x()*pt.x() + pt.y()*pt.y();
+ if (nrm1 > max) {
+ max = nrm1;
+ }
+ }
+ return max;
+ }
+
+ virtual
+ std::vector<std::vector<std::array<double, 3>>>
+ get_features() const
+ {
+ std::vector<std::vector<std::array<double, 3>>> features = {};
+
+ for (const auto & pt: poly_->points) {
+ const double r = pt.x();
+ const double circ = 2 * 3.14159265359 * r;
+ const size_t n = int(circ / edge_size_ - 0.5) + 1;
+ std::vector<std::array<double, 3>> line;
+ for (size_t i=0; i < n; i++) {
+ const double alpha = (2 * 3.14159265359 * i) / n;
+ line.push_back({
+ r * cos(alpha),
+ r * sin(alpha),
+ pt.y()
+ });
+ }
+ line.push_back(line.front());
+ features.push_back(line);
+ }
+ return features;
+ }
+
+ private:
+ const std::shared_ptr<pygalmesh::Polygon2D> poly_;
+ const double edge_size_;
+};
+
+} // namespace pygalmesh
+
+#endif // POLYGON2D_HPP
--- /dev/null
+#ifndef PRIMITIVES_HPP
+#define PRIMITIVES_HPP
+
+#include "domain.hpp"
+
+#include <memory>
+#include <vector>
+
+namespace pygalmesh {
+
+class Ball: public pygalmesh::DomainBase
+{
+ public:
+ Ball(
+ const std::array<double, 3> & x0,
+ const double radius
+ ):
+ x0_(x0),
+ radius_(radius)
+ {
+ assert(x0_.size() == 3);
+ }
+
+ virtual ~Ball() = default;
+
+ virtual
+ double
+ eval(const std::array<double, 3> & x) const
+ {
+ const double xx0 = x[0] - x0_[0];
+ const double yy0 = x[1] - x0_[1];
+ const double zz0 = x[2] - x0_[2];
+ return xx0*xx0 + yy0*yy0 + zz0*zz0 - radius_*radius_;
+ }
+
+ virtual
+ double
+ get_bounding_sphere_squared_radius() const
+ {
+ const double x0_nrm = sqrt(x0_[0]*x0_[0] + x0_[1]*x0_[1] + x0_[2]*x0_[2]);
+ return (x0_nrm + radius_) * (x0_nrm + radius_);
+ }
+
+ private:
+ const std::array<double, 3> x0_;
+ const double radius_;
+};
+
+
+class Cuboid: public pygalmesh::DomainBase
+{
+ public:
+ Cuboid(
+ const std::array<double, 3> & x0,
+ const std::array<double, 3> & x1
+ ):
+ x0_(x0),
+ x1_(x1)
+ {
+ }
+
+ virtual ~Cuboid() = default;
+
+ virtual
+ double
+ eval(const std::array<double, 3> & x) const
+ {
+ // TODO differentiable expression?
+ return std::max(std::max(
+ (x[0] - x0_[0]) * (x[0] - x1_[0]),
+ (x[1] - x0_[1]) * (x[1] - x1_[1])
+ ),
+ (x[2] - x0_[2]) * (x[2] - x1_[2])
+ );
+ }
+
+ virtual
+ double
+ get_bounding_sphere_squared_radius() const
+ {
+ const double x0_nrm2 = x0_[0]*x0_[0] + x0_[1]*x0_[1] + x0_[2]*x0_[2];
+ const double x1_nrm2 = x1_[0]*x1_[0] + x1_[1]*x1_[1] + x1_[2]*x1_[2];
+ return std::max({x0_nrm2, x1_nrm2});
+ }
+
+ virtual
+ std::vector<std::vector<std::array<double, 3>>>
+ get_features() const
+ {
+ std::vector<std::array<double, 3>> corners = {
+ {x0_[0], x0_[1], x0_[2]},
+ {x1_[0], x0_[1], x0_[2]},
+ {x0_[0], x1_[1], x0_[2]},
+ {x0_[0], x0_[1], x1_[2]},
+ {x1_[0], x1_[1], x0_[2]},
+ {x1_[0], x0_[1], x1_[2]},
+ {x0_[0], x1_[1], x1_[2]},
+ {x1_[0], x1_[1], x1_[2]}
+ };
+ return {
+ {corners[0], corners[1]},
+ {corners[0], corners[2]},
+ {corners[0], corners[3]},
+ {corners[1], corners[4]},
+ {corners[1], corners[5]},
+ {corners[2], corners[4]},
+ {corners[2], corners[6]},
+ {corners[3], corners[5]},
+ {corners[3], corners[6]},
+ {corners[4], corners[7]},
+ {corners[5], corners[7]},
+ {corners[6], corners[7]}
+ };
+ };
+
+ private:
+ const std::array<double, 3> x0_;
+ const std::array<double, 3> x1_;
+};
+
+
+class Ellipsoid: public pygalmesh::DomainBase
+{
+ public:
+ Ellipsoid(
+ const std::array<double, 3> & x0,
+ const double a0,
+ const double a1,
+ const double a2
+ ):
+ x0_(x0),
+ a0_2_(a0*a0),
+ a1_2_(a1*a1),
+ a2_2_(a2*a1)
+ {
+ }
+
+ virtual ~Ellipsoid() = default;
+
+ virtual
+ double
+ eval(const std::array<double, 3> & x) const
+ {
+ const double xx0 = x[0] - x0_[0];
+ const double yy0 = x[1] - x0_[1];
+ const double zz0 = x[2] - x0_[2];
+ return xx0*xx0/a0_2_ + yy0*yy0/a1_2_ + zz0*zz0/a2_2_ - 1.0;
+ }
+
+ virtual
+ double
+ get_bounding_sphere_squared_radius() const
+ {
+ return std::max({a0_2_, a1_2_, a2_2_});
+ }
+
+ private:
+ const std::array<double, 3> x0_;
+ const double a0_2_;
+ const double a1_2_;
+ const double a2_2_;
+};
+
+
+class Cylinder: public pygalmesh::DomainBase
+{
+ public:
+ Cylinder(
+ const double z0,
+ const double z1,
+ const double radius,
+ const double feature_edge_length
+ ):
+ z0_(z0),
+ z1_(z1),
+ radius_(radius),
+ feature_edge_length_(feature_edge_length)
+ {
+ assert(z1_ > z0_);
+ }
+
+ virtual ~Cylinder() = default;
+
+ virtual
+ double
+ eval(const std::array<double, 3> & x) const
+ {
+ return (z0_ < x[2] && x[2] < z1_) ?
+ x[0]*x[0] + x[1]*x[1] - radius_*radius_ :
+ 1.0;
+ }
+
+ virtual
+ double
+ get_bounding_sphere_squared_radius() const
+ {
+ const double zmax = std::max({abs(z0_), abs(z1_)});
+ return zmax*zmax + radius_*radius_;
+ }
+
+ virtual
+ std::vector<std::vector<std::array<double, 3>>>
+ get_features() const
+ {
+ const double pi = 3.1415926535897932384;
+ const size_t n = 2 * pi * radius_ / feature_edge_length_;
+ std::vector<std::array<double, 3>> circ0(n+1);
+ std::vector<std::array<double, 3>> circ1(n+1);
+ for (size_t i=0; i < n; i++) {
+ const double c = radius_ * cos((2*pi * i) / n);
+ const double s = radius_ * sin((2*pi * i) / n);
+ circ0[i] = {c, s, z0_};
+ circ1[i] = {c, s, z1_};
+ }
+ // close the circles
+ circ0[n] = circ0[0];
+ circ1[n] = circ1[0];
+ return {circ0, circ1};
+ };
+
+ private:
+ const double z0_;
+ const double z1_;
+ const double radius_;
+ const double feature_edge_length_;
+};
+
+
+class Cone: public pygalmesh::DomainBase
+{
+ public:
+ Cone(
+ const double radius,
+ const double height,
+ const double feature_edge_length
+ ):
+ radius_(radius),
+ height_(height),
+ feature_edge_length_(feature_edge_length)
+ {
+ assert(radius_ > 0.0);
+ assert(height_ > 0.0);
+ }
+
+ virtual ~Cone() = default;
+
+ virtual
+ double
+ eval(const std::array<double, 3> & x) const
+ {
+ const double rad = radius_ * (1.0 - x[2] / height_);
+
+ return (0.0 < x[2] && x[2] < height_) ?
+ x[0]*x[0] + x[1]*x[1] - rad*rad :
+ 1.0;
+ }
+
+ virtual
+ double
+ get_bounding_sphere_squared_radius() const
+ {
+ const double max = std::max({radius_, height_});
+ return max*max;
+ }
+
+ virtual
+ std::vector<std::vector<std::array<double, 3>>>
+ get_features() const
+ {
+ const double pi = 3.1415926535897932384;
+ const size_t n = 2 * pi * radius_ / feature_edge_length_;
+ std::vector<std::array<double, 3>> circ0(n+1);
+ for (size_t i=0; i < n; i++) {
+ const double c = radius_ * cos((2*pi * i) / n);
+ const double s = radius_ * sin((2*pi * i) / n);
+ circ0[i] = {c, s, 0.0};
+ }
+ circ0[n] = circ0[0];
+ return {circ0};
+ };
+
+ private:
+ const double radius_;
+ const double height_;
+ const double feature_edge_length_;
+};
+
+
+class Tetrahedron: public pygalmesh::DomainBase
+{
+ public:
+ Tetrahedron(
+ const std::array<double, 3> & x0,
+ const std::array<double, 3> & x1,
+ 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]))
+ {
+ }
+
+ 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)
+ );
+ }
+
+ 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;
+ }
+
+ virtual
+ double
+ get_bounding_sphere_squared_radius() const
+ {
+ return std::max({
+ x0_.dot(x0_),
+ x1_.dot(x1_),
+ x2_.dot(x2_),
+ x3_.dot(x3_)
+ });
+ }
+
+ virtual
+ std::vector<std::vector<std::array<double, 3>>>
+ get_features() const
+ {
+ std::vector<std::array<double, 3>> pts = {
+ {x0_[0], x0_[1], x0_[2]},
+ {x1_[0], x1_[1], x1_[2]},
+ {x2_[0], x2_[1], x2_[2]},
+ {x3_[0], x3_[1], x3_[2]}
+ };
+ return {
+ {pts[0], pts[1]},
+ {pts[0], pts[2]},
+ {pts[0], pts[3]},
+ {pts[1], pts[2]},
+ {pts[1], pts[3]},
+ {pts[2], pts[3]}
+ };
+ };
+
+ private:
+ const Eigen::Vector3d x0_;
+ const Eigen::Vector3d x1_;
+ const Eigen::Vector3d x2_;
+ const Eigen::Vector3d x3_;
+};
+
+
+class Torus: public pygalmesh::DomainBase
+{
+ public:
+ Torus(
+ const double major_radius,
+ const double minor_radius
+ ):
+ major_radius_(major_radius),
+ minor_radius_(minor_radius)
+ {
+ }
+
+ virtual ~Torus() = default;
+
+ virtual
+ double
+ eval(const std::array<double, 3> & x) const
+ {
+ const double r = sqrt(x[0]*x[0] + x[1]*x[1]);
+ return (
+ (r - major_radius_)*(r - major_radius_) + x[2]*x[2]
+ - minor_radius_*minor_radius_
+ );
+ }
+
+ virtual
+ double
+ get_bounding_sphere_squared_radius() const
+ {
+ return (major_radius_ + minor_radius_)*(major_radius_ + minor_radius_);
+ }
+
+ private:
+ const double major_radius_;
+ const double minor_radius_;
+};
+
+
+class HalfSpace: public pygalmesh::DomainBase
+{
+ public:
+ HalfSpace(
+ const std::array<double, 3> & n,
+ const double alpha,
+ const double bounding_sphere_squared_radius
+ ):
+ n_(n),
+ alpha_(alpha),
+ bounding_sphere_squared_radius_(bounding_sphere_squared_radius)
+ {
+ }
+
+ virtual ~HalfSpace() = default;
+
+ virtual
+ double
+ eval(const std::array<double, 3> & x) const
+ {
+ return n_[0]*x[0] + n_[1]*x[1] + n_[2]*x[2] - alpha_;
+ }
+
+ virtual
+ double
+ get_bounding_sphere_squared_radius() const
+ {
+ return bounding_sphere_squared_radius_;
+ }
+
+ private:
+ const std::array<double, 3> n_;
+ const double alpha_;
+ const double bounding_sphere_squared_radius_;
+};
+
+} // namespace pygalmesh
+
+#endif // PRIMITIVES_HPP
--- /dev/null
+#include "domain.hpp"
+#include "generate.hpp"
+#include "generate_from_off.hpp"
+#include "generate_surface_mesh.hpp"
+#include "polygon2d.hpp"
+#include "primitives.hpp"
+
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+
+namespace py = pybind11;
+using namespace pygalmesh;
+
+
+// https://pybind11.readthedocs.io/en/stable/advanced/classes.html#overriding-virtual-functions-in-python
+class PyDomainBase: public DomainBase {
+public:
+ using DomainBase::DomainBase;
+
+ double
+ eval(const std::array<double, 3> & x) const override {
+ PYBIND11_OVERLOAD_PURE(double, DomainBase, eval, x);
+ }
+
+ double
+ get_bounding_sphere_squared_radius() const override {
+ PYBIND11_OVERLOAD_PURE(double, DomainBase, get_bounding_sphere_squared_radius);
+ }
+
+ // std::vector<std::vector<std::array<double, 3>>>
+ // get_features() const override {
+ // PYBIND11_OVERLOAD(
+ // std::vector<std::vector<std::array<double, 3>>>,
+ // DomainBase,
+ // get_features,
+ // 0.0
+ // );
+ // }
+};
+
+
+PYBIND11_MODULE(_pygalmesh, m) {
+ // m.doc() = "documentation string";
+
+ // Domain base.
+ // shared_ptr b/c of
+ // <https://github.com/pybind/pybind11/issues/956#issuecomment-317022720>
+ py::class_<DomainBase, PyDomainBase, std::shared_ptr<DomainBase>>(m, "DomainBase")
+ .def(py::init<>())
+ .def("eval", &DomainBase::eval)
+ .def("get_bounding_sphere_squared_radius", &DomainBase::get_bounding_sphere_squared_radius)
+ .def("get_features", &DomainBase::get_features);
+
+ // Domain transformations
+ py::class_<Translate, DomainBase, std::shared_ptr<Translate>>(m, "Translate")
+ .def(py::init<
+ const std::shared_ptr<const DomainBase> &,
+ const std::array<double, 3> &
+ >())
+ .def("eval", &Translate::eval)
+ .def("translate_features", &Translate::translate_features)
+ .def("get_bounding_sphere_squared_radius", &Translate::get_bounding_sphere_squared_radius)
+ .def("get_features", &Translate::get_features);
+
+ py::class_<Rotate, DomainBase, std::shared_ptr<Rotate>>(m, "Rotate")
+ .def(py::init<
+ const std::shared_ptr<const pygalmesh::DomainBase> &,
+ const std::array<double, 3> &,
+ const double
+ >())
+ .def("eval", &Rotate::eval)
+ .def("rotate", &Rotate::rotate)
+ .def("rotate_features", &Rotate::rotate_features)
+ .def("get_bounding_sphere_squared_radius", &Rotate::get_bounding_sphere_squared_radius)
+ .def("get_features", &Rotate::get_features);
+
+ py::class_<Scale, DomainBase, std::shared_ptr<Scale>>(m, "Scale")
+ .def(py::init<
+ std::shared_ptr<const pygalmesh::DomainBase> &,
+ const double
+ >())
+ .def("eval", &Scale::eval)
+ .def("scale_features", &Scale::scale_features)
+ .def("get_bounding_sphere_squared_radius", &Scale::get_bounding_sphere_squared_radius)
+ .def("get_features", &Scale::get_features);
+
+ py::class_<Stretch, DomainBase, std::shared_ptr<Stretch>>(m, "Stretch")
+ .def(py::init<
+ std::shared_ptr<const pygalmesh::DomainBase> &,
+ const std::array<double, 3> &
+ >())
+ .def("eval", &Stretch::eval)
+ .def("stretch_features", &Stretch::stretch_features)
+ .def("get_bounding_sphere_squared_radius", &Stretch::get_bounding_sphere_squared_radius)
+ .def("get_features", &Stretch::get_features);
+
+ py::class_<Intersection, DomainBase, std::shared_ptr<Intersection>>(m, "Intersection")
+ .def(py::init<
+ std::vector<std::shared_ptr<const pygalmesh::DomainBase>> &
+ >())
+ .def("eval", &Intersection::eval)
+ .def("get_bounding_sphere_squared_radius", &Intersection::get_bounding_sphere_squared_radius)
+ .def("get_features", &Intersection::get_features);
+
+ py::class_<Union, DomainBase, std::shared_ptr<Union>>(m, "Union")
+ .def(py::init<
+ std::vector<std::shared_ptr<const pygalmesh::DomainBase>> &
+ >())
+ .def("eval", &Union::eval)
+ .def("get_bounding_sphere_squared_radius", &Union::get_bounding_sphere_squared_radius)
+ .def("get_features", &Union::get_features);
+
+ py::class_<Difference, DomainBase, std::shared_ptr<Difference>>(m, "Difference")
+ .def(py::init<
+ std::shared_ptr<const pygalmesh::DomainBase> &,
+ std::shared_ptr<const pygalmesh::DomainBase> &
+ >())
+ .def("eval", &Difference::eval)
+ .def("get_bounding_sphere_squared_radius", &Difference::get_bounding_sphere_squared_radius)
+ .def("get_features", &Difference::get_features);
+
+ // Primitives
+ py::class_<Ball, DomainBase, std::shared_ptr<Ball>>(m, "Ball")
+ .def(py::init<
+ const std::array<double, 3> &,
+ const double
+ >())
+ .def("eval", &Ball::eval)
+ .def("get_bounding_sphere_squared_radius", &Ball::get_bounding_sphere_squared_radius);
+
+ py::class_<Cuboid, DomainBase, std::shared_ptr<Cuboid>>(m, "Cuboid")
+ .def(py::init<
+ const std::array<double, 3> &,
+ const std::array<double, 3> &
+ >())
+ .def("eval", &Cuboid::eval)
+ .def("get_bounding_sphere_squared_radius", &Cuboid::get_bounding_sphere_squared_radius)
+ .def("get_features", &Cuboid::get_features);
+
+ py::class_<Ellipsoid, DomainBase, std::shared_ptr<Ellipsoid>>(m, "Ellipsoid")
+ .def(py::init<
+ const std::array<double, 3> &,
+ const double,
+ const double,
+ const double
+ >())
+ .def("eval", &Ellipsoid::eval)
+ .def("get_bounding_sphere_squared_radius", &Ellipsoid::get_bounding_sphere_squared_radius)
+ .def("get_features", &Ellipsoid::get_features);
+
+ py::class_<Cylinder, DomainBase, std::shared_ptr<Cylinder>>(m, "Cylinder")
+ .def(py::init<
+ const double,
+ const double,
+ const double,
+ const double
+ >())
+ .def("eval", &Cylinder::eval)
+ .def("get_bounding_sphere_squared_radius", &Cylinder::get_bounding_sphere_squared_radius)
+ .def("get_features", &Cylinder::get_features);
+
+ py::class_<Cone, DomainBase, std::shared_ptr<Cone>>(m, "Cone")
+ .def(py::init<
+ const double,
+ const double,
+ const double
+ >())
+ .def("eval", &Cone::eval)
+ .def("get_bounding_sphere_squared_radius", &Cone::get_bounding_sphere_squared_radius)
+ .def("get_features", &Cone::get_features);
+
+ py::class_<Tetrahedron, DomainBase, std::shared_ptr<Tetrahedron>>(m, "Tetrahedron")
+ .def(py::init<
+ const std::array<double, 3> &,
+ const std::array<double, 3> &,
+ const std::array<double, 3> &,
+ const std::array<double, 3> &
+ >())
+ .def("eval", &Tetrahedron::eval)
+ .def("get_bounding_sphere_squared_radius", &Tetrahedron::get_bounding_sphere_squared_radius)
+ .def("get_features", &Tetrahedron::get_features);
+
+ py::class_<Torus, DomainBase, std::shared_ptr<Torus>>(m, "Torus")
+ .def(py::init<
+ const double,
+ const double
+ >())
+ .def("eval", &Torus::eval)
+ .def("get_bounding_sphere_squared_radius", &Torus::get_bounding_sphere_squared_radius)
+ .def("get_features", &Torus::get_features);
+
+ py::class_<HalfSpace, DomainBase, std::shared_ptr<HalfSpace>>(m, "HalfSpace")
+ .def(py::init<
+ const std::array<double, 3> &,
+ const double,
+ const double
+ >())
+ .def("eval", &HalfSpace::eval)
+ .def("get_bounding_sphere_squared_radius", &HalfSpace::get_bounding_sphere_squared_radius);
+
+ // polygon2d
+ py::class_<Polygon2D, std::shared_ptr<Polygon2D>>(m, "Polygon2D")
+ .def(py::init<
+ const std::vector<std::array<double, 2>> &
+ >())
+ .def("vector_to_cgal_points", &Polygon2D::vector_to_cgal_points)
+ .def("is_inside", &Polygon2D::is_inside);
+
+ py::class_<Extrude, DomainBase, std::shared_ptr<Extrude>>(m, "Extrude")
+ .def(py::init<
+ const std::shared_ptr<pygalmesh::Polygon2D> &,
+ const std::array<double, 3> &,
+ const double,
+ const double
+ >(),
+ py::arg("poly"),
+ py::arg("direction"),
+ py::arg("alpha") = 0.0,
+ py::arg("edge_size") = 0.0
+ )
+ .def("eval", &Extrude::eval)
+ .def("get_bounding_sphere_squared_radius", &Extrude::get_bounding_sphere_squared_radius)
+ .def("get_features", &Extrude::get_features);
+
+ py::class_<ring_extrude, DomainBase, std::shared_ptr<ring_extrude>>(m, "RingExtrude")
+ .def(py::init<
+ const std::shared_ptr<pygalmesh::Polygon2D> &,
+ const double
+ >())
+ .def("eval", &ring_extrude::eval)
+ .def("get_bounding_sphere_squared_radius", &ring_extrude::get_bounding_sphere_squared_radius)
+ .def("get_features", &ring_extrude::get_features);
+
+ // functions
+ m.def("generate_from_off", &generate_from_off);
+ m.def(
+ "generate_mesh", &generate_mesh,
+ py::arg("domain"),
+ py::arg("outfile"),
+ py::arg("feature_edges") = std::vector<std::vector<std::array<double, 3>>>(),
+ py::arg("bounding_sphere_radius") = 0.0,
+ py::arg("boundary_precision") = 1.0e-4,
+ py::arg("lloyd") = false,
+ py::arg("odt") = false,
+ py::arg("perturb") = true,
+ py::arg("exude") = true,
+ py::arg("edge_size") = 0.0,
+ py::arg("facet_angle") = 0.0,
+ py::arg("facet_size") = 0.0,
+ py::arg("facet_distance") = 0.0,
+ py::arg("cell_radius_edge_ratio") = 0.0,
+ py::arg("cell_size") = 0.0,
+ py::arg("verbose") = true
+ );
+ m.def(
+ "generate_surface_mesh", &generate_surface_mesh,
+ py::arg("domain"),
+ py::arg("outfile"),
+ py::arg("bounding_sphere_radius") = 0.0,
+ py::arg("angle_bound") = 0.0,
+ py::arg("radius_bound") = 0.0,
+ py::arg("distance_bound") = 0.0,
+ py::arg("verbose") = true
+ );
+}
--- /dev/null
+# -*- coding: utf-8 -*-
+#
+import numpy
+import meshio
+
+# pylint: disable=import-error
+import pygalmesh
+
+
+def _row_dot(a, b):
+ # http://stackoverflow.com/a/26168677/353337
+ return numpy.einsum("ij, ij->i", a, b)
+
+
+def compute_volumes(vertices, tets):
+ cell_coords = vertices[tets]
+
+ a = cell_coords[:, 1, :] - cell_coords[:, 0, :]
+ b = cell_coords[:, 2, :] - cell_coords[:, 0, :]
+ c = cell_coords[:, 3, :] - cell_coords[:, 0, :]
+
+ # omega = <a, b x c>
+ omega = _row_dot(a, numpy.cross(b, c))
+
+ # https://en.wikipedia.org/wiki/Tetrahedron#Volume
+ return abs(omega) / 6.0
+
+
+def compute_triangle_areas(vertices, triangles):
+ e0 = vertices[triangles[:, 1]] - vertices[triangles[:, 0]]
+ e1 = vertices[triangles[:, 2]] - vertices[triangles[:, 1]]
+
+ e0_cross_e1 = numpy.cross(e0, e1)
+ areas = 0.5 * numpy.sqrt(_row_dot(e0_cross_e1, e0_cross_e1))
+
+ return areas
+
+
+def test_ball():
+ s = pygalmesh.Ball([0.0, 0.0, 0.0], 1.0)
+ pygalmesh.generate_mesh(s, "out.mesh", cell_size=0.2, verbose=False)
+
+ mesh = meshio.read("out.mesh")
+
+ assert abs(max(mesh.points[:, 0]) - 1.0) < 0.02
+ assert abs(min(mesh.points[:, 0]) + 1.0) < 0.02
+ assert abs(max(mesh.points[:, 1]) - 1.0) < 0.02
+ assert abs(min(mesh.points[:, 1]) + 1.0) < 0.02
+ assert abs(max(mesh.points[:, 2]) - 1.0) < 0.02
+ assert abs(min(mesh.points[:, 2]) + 1.0) < 0.02
+
+ vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"]))
+ assert abs(vol - 4.0 / 3.0 * numpy.pi) < 0.15
+ return
+
+
+def test_balls_union():
+ radius = 1.0
+ displacement = 0.5
+ s0 = pygalmesh.Ball([displacement, 0, 0], radius)
+ s1 = pygalmesh.Ball([-displacement, 0, 0], radius)
+ u = pygalmesh.Union([s0, s1])
+
+ a = numpy.sqrt(radius ** 2 - displacement ** 2)
+ edge_size = 0.1
+ n = int(2 * numpy.pi * a / edge_size)
+ circ = [
+ [0.0, a * numpy.cos(i * 2 * numpy.pi / n), a * numpy.sin(i * 2 * numpy.pi / n)]
+ for i in range(n)
+ ]
+ circ.append(circ[0])
+
+ pygalmesh.generate_mesh(
+ u,
+ "out.mesh",
+ feature_edges=[circ],
+ cell_size=0.15,
+ edge_size=edge_size,
+ verbose=False,
+ )
+
+ mesh = meshio.read("out.mesh")
+
+ assert abs(max(mesh.points[:, 0]) - (radius + displacement)) < 0.02
+ assert abs(min(mesh.points[:, 0]) + (radius + displacement)) < 0.02
+ assert abs(max(mesh.points[:, 1]) - radius) < 0.02
+ assert abs(min(mesh.points[:, 1]) + radius) < 0.02
+ assert abs(max(mesh.points[:, 2]) - radius) < 0.02
+ assert abs(min(mesh.points[:, 2]) + radius) < 0.02
+
+ vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"]))
+ h = radius - displacement
+ ref_vol = 2 * (
+ 4.0 / 3.0 * numpy.pi * radius ** 3 - h * numpy.pi / 6.0 * (3 * a ** 2 + h ** 2)
+ )
+
+ assert abs(vol - ref_vol) < 0.1
+
+ return
+
+
+def test_balls_intersection():
+ radius = 1.0
+ displacement = 0.5
+ s0 = pygalmesh.Ball([displacement, 0, 0], radius)
+ s1 = pygalmesh.Ball([-displacement, 0, 0], radius)
+ u = pygalmesh.Intersection([s0, s1])
+
+ a = numpy.sqrt(radius ** 2 - displacement ** 2)
+ edge_size = 0.1
+ n = int(2 * numpy.pi * a / edge_size)
+ circ = [
+ [0.0, a * numpy.cos(i * 2 * numpy.pi / n), a * numpy.sin(i * 2 * numpy.pi / n)]
+ for i in range(n)
+ ]
+ circ.append(circ[0])
+
+ pygalmesh.generate_mesh(
+ u,
+ "out.mesh",
+ feature_edges=[circ],
+ cell_size=0.15,
+ edge_size=edge_size,
+ verbose=False,
+ )
+
+ mesh = meshio.read("out.mesh")
+
+ assert abs(max(mesh.points[:, 0]) - (radius - displacement)) < 0.02
+ assert abs(min(mesh.points[:, 0]) + (radius - displacement)) < 0.02
+ assert abs(max(mesh.points[:, 1]) - a) < 0.02
+ assert abs(min(mesh.points[:, 1]) + a) < 0.02
+ assert abs(max(mesh.points[:, 2]) - a) < 0.02
+ assert abs(min(mesh.points[:, 2]) + a) < 0.02
+
+ vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"]))
+ h = radius - displacement
+ ref_vol = 2 * (h * numpy.pi / 6.0 * (3 * a ** 2 + h ** 2))
+
+ assert abs(vol - ref_vol) < 0.1
+
+ return
+
+
+# pylint: disable=too-many-locals
+def test_balls_difference():
+ radius = 1.0
+ displacement = 0.5
+ s0 = pygalmesh.Ball([displacement, 0, 0], radius)
+ s1 = pygalmesh.Ball([-displacement, 0, 0], radius)
+ u = pygalmesh.Difference(s0, s1)
+
+ a = numpy.sqrt(radius ** 2 - displacement ** 2)
+ edge_size = 0.15
+ n = int(2 * numpy.pi * a / edge_size)
+ circ = [
+ [0.0, a * numpy.cos(i * 2 * numpy.pi / n), a * numpy.sin(i * 2 * numpy.pi / n)]
+ for i in range(n)
+ ]
+ circ.append(circ[0])
+
+ pygalmesh.generate_mesh(
+ u,
+ "out.mesh",
+ feature_edges=[circ],
+ cell_size=0.15,
+ edge_size=edge_size,
+ facet_angle=25,
+ facet_size=0.15,
+ cell_radius_edge_ratio=2.0,
+ verbose=False,
+ )
+
+ mesh = meshio.read("out.mesh")
+
+ tol = 0.02
+ assert abs(max(mesh.points[:, 0]) - (radius + displacement)) < tol
+ assert abs(min(mesh.points[:, 0]) - 0.0) < tol
+ assert abs(max(mesh.points[:, 1]) - radius) < tol
+ assert abs(min(mesh.points[:, 1]) + radius) < tol
+ assert abs(max(mesh.points[:, 2]) - radius) < tol
+ assert abs(min(mesh.points[:, 2]) + radius) < tol
+
+ vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"]))
+ h = radius - displacement
+ ref_vol = 4.0 / 3.0 * numpy.pi * radius ** 3 - 2 * h * numpy.pi / 6.0 * (
+ 3 * a ** 2 + h ** 2
+ )
+
+ assert abs(vol - ref_vol) < 0.05
+
+ return
+
+
+def test_cuboids_intersection():
+ c0 = pygalmesh.Cuboid([0, 0, -0.5], [3, 3, 0.5])
+ c1 = pygalmesh.Cuboid([1, 1, -2], [2, 2, 2])
+ u = pygalmesh.Intersection([c0, c1])
+
+ # In CGAL, feature edges must not intersect, and that's a problem here: The
+ # intersection edges of the cuboids share eight points with the edges of
+ # the tall and skinny cuboid.
+ # eps = 1.0e-2
+ # extra_features = [
+ # [[1.0, 1.0 + eps, 0.5], [1.0, 2.0 - eps, 0.5]],
+ # [[1.0 + eps, 2.0, 0.5], [2.0 - eps, 2.0, 0.5]],
+ # [[2.0, 2.0 - eps, 0.5], [2.0, 1.0 + eps, 0.5]],
+ # [[2.0 - eps, 1.0, 0.5], [1.0 + eps, 1.0, 0.5]],
+ # ]
+
+ pygalmesh.generate_mesh(u, "out.mesh", cell_size=0.1, edge_size=0.1, verbose=False)
+
+ mesh = meshio.read("out.mesh")
+
+ # filter the vertices that belong to cells
+ verts = mesh.points[numpy.unique(mesh.cells["tetra"])]
+
+ tol = 1.0e-2
+ assert abs(max(verts[:, 0]) - 2.0) < tol
+ assert abs(min(verts[:, 0]) - 1.0) < tol
+ assert abs(max(verts[:, 1]) - 2.0) < tol
+ assert abs(min(verts[:, 1]) - 1.0) < tol
+ assert abs(max(verts[:, 2]) - 0.5) < 0.05
+ assert abs(min(verts[:, 2]) + 0.5) < 0.05
+
+ vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"]))
+ assert abs(vol - 1.0) < 0.05
+
+ return
+
+
+def test_cuboids_union():
+ c0 = pygalmesh.Cuboid([0, 0, -0.5], [3, 3, 0.5])
+ c1 = pygalmesh.Cuboid([1, 1, -2], [2, 2, 2])
+ u = pygalmesh.Union([c0, c1])
+
+ pygalmesh.generate_mesh(u, "out.mesh", cell_size=0.2, edge_size=0.2, verbose=False)
+
+ mesh = meshio.read("out.mesh")
+
+ # filter the vertices that belong to cells
+ verts = mesh.points[numpy.unique(mesh.cells["tetra"])]
+
+ tol = 1.0e-2
+ assert abs(max(verts[:, 0]) - 3.0) < tol
+ assert abs(min(verts[:, 0]) - 0.0) < tol
+ assert abs(max(verts[:, 1]) - 3.0) < tol
+ assert abs(min(verts[:, 1]) - 0.0) < tol
+ assert abs(max(verts[:, 2]) - 2.0) < tol
+ assert abs(min(verts[:, 2]) + 2.0) < tol
+
+ vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"]))
+ assert abs(vol - 12.0) < 0.1
+ return
+
+
+def test_cuboid():
+ s0 = pygalmesh.Cuboid([0, 0, 0], [1, 2, 3])
+ pygalmesh.generate_mesh(s0, "out.mesh", edge_size=0.1, verbose=False)
+
+ mesh = meshio.read("out.mesh")
+
+ tol = 1.0e-3
+ assert abs(max(mesh.points[:, 0]) - 1.0) < tol
+ assert abs(min(mesh.points[:, 0]) + 0.0) < tol
+ assert abs(max(mesh.points[:, 1]) - 2.0) < tol
+ assert abs(min(mesh.points[:, 1]) + 0.0) < tol
+ assert abs(max(mesh.points[:, 2]) - 3.0) < tol
+ assert abs(min(mesh.points[:, 2]) + 0.0) < tol
+
+ vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"]))
+ assert abs(vol - 6.0) < tol
+ return
+
+
+def test_cone():
+ base_radius = 1.0
+ height = 2.0
+ edge_size = 0.1
+ s0 = pygalmesh.Cone(base_radius, height, edge_size)
+ pygalmesh.generate_mesh(
+ s0, "out.mesh", cell_size=0.1, edge_size=edge_size, verbose=False
+ )
+
+ mesh = meshio.read("out.mesh")
+
+ tol = 2.0e-1
+ assert abs(max(mesh.points[:, 0]) - base_radius) < tol
+ assert abs(min(mesh.points[:, 0]) + base_radius) < tol
+ assert abs(max(mesh.points[:, 1]) - base_radius) < tol
+ assert abs(min(mesh.points[:, 1]) + base_radius) < tol
+ assert abs(max(mesh.points[:, 2]) - height) < tol
+ assert abs(min(mesh.points[:, 2]) + 0.0) < tol
+
+ vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"]))
+ ref_vol = numpy.pi * base_radius * base_radius / 3.0 * height
+ assert abs(vol - ref_vol) < tol
+ return
+
+
+def test_cylinder():
+ radius = 1.0
+ z0 = 0.0
+ z1 = 1.0
+ edge_length = 0.1
+ s0 = pygalmesh.Cylinder(z0, z1, radius, edge_length)
+ pygalmesh.generate_mesh(
+ s0, "out.mesh", cell_size=0.1, edge_size=edge_length, verbose=False
+ )
+
+ mesh = meshio.read("out.mesh")
+
+ tol = 1.0e-1
+ assert abs(max(mesh.points[:, 0]) - radius) < tol
+ assert abs(min(mesh.points[:, 0]) + radius) < tol
+ assert abs(max(mesh.points[:, 1]) - radius) < tol
+ assert abs(min(mesh.points[:, 1]) + radius) < tol
+ assert abs(max(mesh.points[:, 2]) - z1) < tol
+ assert abs(min(mesh.points[:, 2]) + z0) < tol
+
+ vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"]))
+ ref_vol = numpy.pi * radius * radius * (z1 - z0)
+ assert abs(vol - ref_vol) < tol
+ return
+
+
+def test_tetrahedron():
+ s0 = pygalmesh.Tetrahedron(
+ [0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]
+ )
+ pygalmesh.generate_mesh(s0, "out.mesh", cell_size=0.1, edge_size=0.1, verbose=False)
+
+ mesh = meshio.read("out.mesh")
+
+ tol = 1.0e-4
+ assert abs(max(mesh.points[:, 0]) - 1.0) < tol
+ assert abs(min(mesh.points[:, 0]) + 0.0) < tol
+ assert abs(max(mesh.points[:, 1]) - 1.0) < tol
+ assert abs(min(mesh.points[:, 1]) + 0.0) < tol
+ assert abs(max(mesh.points[:, 2]) - 1.0) < tol
+ assert abs(min(mesh.points[:, 2]) + 0.0) < tol
+
+ vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"]))
+ assert abs(vol - 1.0 / 6.0) < tol
+ return
+
+
+def test_torus():
+ major_radius = 1.0
+ minor_radius = 0.5
+ s0 = pygalmesh.Torus(major_radius, minor_radius)
+ pygalmesh.generate_mesh(s0, "out.mesh", cell_size=0.1, verbose=False)
+
+ mesh = meshio.read("out.mesh")
+
+ tol = 1.0e-2
+ radii_sum = major_radius + minor_radius
+ assert abs(max(mesh.points[:, 0]) - radii_sum) < tol
+ assert abs(min(mesh.points[:, 0]) + radii_sum) < tol
+ assert abs(max(mesh.points[:, 1]) - radii_sum) < tol
+ assert abs(min(mesh.points[:, 1]) + radii_sum) < tol
+ assert abs(max(mesh.points[:, 2]) - minor_radius) < tol
+ assert abs(min(mesh.points[:, 2]) + minor_radius) < tol
+
+ vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"]))
+ ref_vol = (numpy.pi * minor_radius * minor_radius) * (2 * numpy.pi * major_radius)
+ assert abs(vol - ref_vol) < 1.0e-1
+ return
+
+
+def test_custom_function():
+ class Hyperboloid(pygalmesh.DomainBase):
+ def __init__(self, edge_size):
+ super(Hyperboloid, self).__init__()
+ self.z0 = -1.0
+ self.z1 = 1.0
+ self.waist_radius = 0.5
+ self.edge_size = edge_size
+ return
+
+ def eval(self, x):
+ if self.z0 < x[2] and x[2] < self.z1:
+ return x[0] ** 2 + x[1] ** 2 - (x[2] ** 2 + self.waist_radius) ** 2
+ return 1.0
+
+ def get_bounding_sphere_squared_radius(self):
+ z_max = max(abs(self.z0), abs(self.z1))
+ r_max = z_max ** 2 + self.waist_radius
+ return r_max * r_max + z_max * z_max
+
+ def get_features(self):
+ radius0 = self.z0 ** 2 + self.waist_radius
+ n0 = int(2 * numpy.pi * radius0 / self.edge_size)
+ circ0 = [
+ [
+ radius0 * numpy.cos((2 * numpy.pi * k) / n0),
+ radius0 * numpy.sin((2 * numpy.pi * k) / n0),
+ self.z0,
+ ]
+ for k in range(n0)
+ ]
+ circ0.append(circ0[0])
+
+ radius1 = self.z1 ** 2 + self.waist_radius
+ n1 = int(2 * numpy.pi * radius1 / self.edge_size)
+ circ1 = [
+ [
+ radius1 * numpy.cos((2 * numpy.pi * k) / n1),
+ radius1 * numpy.sin((2 * numpy.pi * k) / n1),
+ self.z1,
+ ]
+ for k in range(n1)
+ ]
+ circ1.append(circ1[0])
+ return [circ0, circ1]
+
+ edge_size = 0.12
+ d = Hyperboloid(edge_size)
+
+ pygalmesh.generate_mesh(
+ d, "out.mesh", cell_size=0.1, edge_size=edge_size, verbose=False
+ )
+
+ mesh = meshio.read("out.mesh")
+
+ # TODO check the reference values
+ tol = 1.0e-1
+ assert abs(max(mesh.points[:, 0]) - 1.4) < tol
+ assert abs(min(mesh.points[:, 0]) + 1.4) < tol
+ assert abs(max(mesh.points[:, 1]) - 1.4) < tol
+ assert abs(min(mesh.points[:, 1]) + 1.4) < tol
+ assert abs(max(mesh.points[:, 2]) - 1.0) < tol
+ assert abs(min(mesh.points[:, 2]) + 1.0) < tol
+
+ vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"]))
+ assert abs(vol - 2 * numpy.pi * 47.0 / 60.0) < 0.15
+ return
+
+
+def test_scaling():
+ alpha = 1.3
+ s = pygalmesh.Scale(pygalmesh.Cuboid([0, 0, 0], [1, 2, 3]), alpha)
+ pygalmesh.generate_mesh(s, "out.mesh", cell_size=0.2, edge_size=0.1, verbose=False)
+
+ mesh = meshio.read("out.mesh")
+
+ tol = 1.0e-3
+ assert abs(max(mesh.points[:, 0]) - 1 * alpha) < tol
+ assert abs(min(mesh.points[:, 0]) + 0.0) < tol
+ assert abs(max(mesh.points[:, 1]) - 2 * alpha) < tol
+ assert abs(min(mesh.points[:, 1]) + 0.0) < tol
+ assert abs(max(mesh.points[:, 2]) - 3 * alpha) < tol
+ assert abs(min(mesh.points[:, 2]) + 0.0) < tol
+
+ vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"]))
+ assert abs(vol - 6.0 * alpha ** 3) < tol
+ return
+
+
+def test_stretch():
+ alpha = 2.0
+ s = pygalmesh.Stretch(pygalmesh.Cuboid([0, 0, 0], [1, 2, 3]), [alpha, 0.0, 0.0])
+ pygalmesh.generate_mesh(s, "out.mesh", cell_size=0.2, edge_size=0.2, verbose=False)
+
+ mesh = meshio.read("out.mesh")
+
+ tol = 1.0e-3
+ assert abs(max(mesh.points[:, 0]) - alpha) < tol
+ assert abs(min(mesh.points[:, 0]) + 0.0) < tol
+ assert abs(max(mesh.points[:, 1]) - 2.0) < tol
+ assert abs(min(mesh.points[:, 1]) + 0.0) < tol
+ assert abs(max(mesh.points[:, 2]) - 3.0) < tol
+ assert abs(min(mesh.points[:, 2]) + 0.0) < tol
+
+ vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"]))
+ assert abs(vol - 12.0) < tol
+
+ return
+
+
+def test_rotation():
+ s0 = pygalmesh.Rotate(
+ pygalmesh.Cuboid([0, 0, 0], [1, 2, 3]), [1.0, 0.0, 0.0], numpy.pi / 12.0
+ )
+ pygalmesh.generate_mesh(s0, "out.mesh", cell_size=0.1, edge_size=0.1, verbose=False)
+
+ mesh = meshio.read("out.mesh")
+
+ tol = 1.0e-3
+ vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"]))
+ assert abs(vol - 6.0) < tol
+ return
+
+
+def test_translation():
+ s0 = pygalmesh.Translate(pygalmesh.Cuboid([0, 0, 0], [1, 2, 3]), [1.0, 0.0, 0.0])
+ pygalmesh.generate_mesh(s0, "out.mesh", cell_size=0.1, edge_size=0.1, verbose=False)
+
+ mesh = meshio.read("out.mesh")
+
+ tol = 1.0e-3
+ assert abs(max(mesh.points[:, 0]) - 2.0) < tol
+ assert abs(min(mesh.points[:, 0]) - 1.0) < tol
+ assert abs(max(mesh.points[:, 1]) - 2.0) < tol
+ assert abs(min(mesh.points[:, 1]) + 0.0) < tol
+ assert abs(max(mesh.points[:, 2]) - 3.0) < tol
+ assert abs(min(mesh.points[:, 2]) + 0.0) < tol
+ vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"]))
+ assert abs(vol - 6.0) < tol
+ return
+
+
+# # segfaults on travis, works locally
+# def test_off():
+# pygalmesh.generate_from_off(
+# 'elephant.off',
+# 'out.mesh',
+# facet_angle=25.0,
+# facet_size=0.15,
+# facet_distance=0.008,
+# cell_radius_edge_ratio=3.0,
+# verbose=False
+# )
+#
+# vertices, cells, _, _, _ = meshio.read('out.mesh')
+#
+# tol = 1.0e-3
+# assert abs(max(vertices[:, 0]) - 0.357612477657) < tol
+# assert abs(min(vertices[:, 0]) + 0.358747130015) < tol
+# assert abs(max(vertices[:, 1]) - 0.496137874959) < tol
+# assert abs(min(vertices[:, 1]) + 0.495301051456) < tol
+# assert abs(max(vertices[:, 2]) - 0.298780230629) < tol
+# assert abs(min(vertices[:, 2]) + 0.300472866512) < tol
+#
+# vol = sum(compute_volumes(vertices, cells['tetra']))
+# assert abs(vol - 0.044164693065) < tol
+# return
+
+
+def test_extrude():
+ p = pygalmesh.Polygon2D([[-0.5, -0.3], [0.5, -0.3], [0.0, 0.5]])
+ domain = pygalmesh.Extrude(p, [0.0, 0.3, 1.0])
+ pygalmesh.generate_mesh(
+ domain, "out.mesh", cell_size=0.1, edge_size=0.1, verbose=False
+ )
+
+ mesh = meshio.read("out.mesh")
+
+ tol = 1.0e-3
+ assert abs(max(mesh.points[:, 0]) - 0.5) < tol
+ assert abs(min(mesh.points[:, 0]) + 0.5) < tol
+ assert abs(max(mesh.points[:, 1]) - 0.8) < tol
+ assert abs(min(mesh.points[:, 1]) + 0.3) < tol
+ assert abs(max(mesh.points[:, 2]) - 1.0) < tol
+ assert abs(min(mesh.points[:, 2]) + 0.0) < tol
+
+ vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"]))
+ assert abs(vol - 0.4) < tol
+ return
+
+
+def test_extrude_rotate():
+ p = pygalmesh.Polygon2D([[-0.5, -0.3], [0.5, -0.3], [0.0, 0.5]])
+ edge_size = 0.1
+ domain = pygalmesh.Extrude(p, [0.0, 0.0, 1.0], 0.5 * 3.14159265359, edge_size)
+ pygalmesh.generate_mesh(
+ domain, "out.mesh", cell_size=0.1, edge_size=edge_size, verbose=False
+ )
+
+ mesh = meshio.read("out.mesh")
+
+ tol = 1.0e-3
+ assert abs(max(mesh.points[:, 0]) - 0.583012701892) < tol
+ assert abs(min(mesh.points[:, 0]) + 0.5) < tol
+ assert abs(max(mesh.points[:, 1]) - 0.5) < tol
+ assert abs(min(mesh.points[:, 1]) + 0.583012701892) < tol
+ assert abs(max(mesh.points[:, 2]) - 1.0) < tol
+ assert abs(min(mesh.points[:, 2]) + 0.0) < tol
+
+ vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"]))
+ assert abs(vol - 0.4) < 0.05
+ return
+
+
+def test_ring_extrude():
+ p = pygalmesh.Polygon2D([[0.5, -0.3], [1.5, -0.3], [1.0, 0.5]])
+ edge_size = 0.1
+ domain = pygalmesh.RingExtrude(p, edge_size)
+ pygalmesh.generate_mesh(
+ domain, "out.mesh", cell_size=0.1, edge_size=edge_size, verbose=False
+ )
+
+ mesh = meshio.read("out.mesh")
+
+ tol = 1.0e-3
+ assert abs(max(mesh.points[:, 0]) - 1.5) < tol
+ assert abs(min(mesh.points[:, 0]) + 1.5) < tol
+ assert abs(max(mesh.points[:, 1]) - 1.5) < tol
+ assert abs(min(mesh.points[:, 1]) + 1.5) < tol
+ assert abs(max(mesh.points[:, 2]) - 0.5) < tol
+ assert abs(min(mesh.points[:, 2]) + 0.3) < tol
+
+ vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"]))
+ assert abs(vol - 2 * numpy.pi * 0.4) < 0.05
+ return
+
+
+# def test_heart():
+# class Heart(pygalmesh.DomainBase):
+# def __init__(self, edge_size):
+# super(Heart, self).__init__()
+# return
+#
+# def eval(self, x):
+# return (x[0]**2 + 9.0/4.0 * x[1]**2 + x[2]**2 - 1)**3 \
+# - x[0]**2 * x[2]**3 - 9.0/80.0 * x[1]**2 * x[2]**3
+#
+# def get_bounding_sphere_squared_radius(self):
+# return 10.0
+#
+# edge_size = 0.1
+# d = Heart(edge_size)
+#
+# pygalmesh.generate_mesh(
+# d,
+# 'out.mesh',
+# cell_size=0.1,
+# edge_size=edge_size,
+# # odt=True,
+# # lloyd=True,
+# # verbose=True
+# )
+#
+# return
+
+
+def test_sphere():
+ radius = 1.0
+ s = pygalmesh.Ball([0.0, 0.0, 0.0], radius)
+ pygalmesh.generate_surface_mesh(
+ s,
+ "out.off",
+ angle_bound=30,
+ radius_bound=0.1,
+ distance_bound=0.1,
+ verbose=False,
+ )
+
+ mesh = meshio.read("out.off")
+
+ tol = 1.0e-2
+ assert abs(max(mesh.points[:, 0]) - radius) < tol
+ assert abs(min(mesh.points[:, 0]) + radius) < tol
+ assert abs(max(mesh.points[:, 1]) - radius) < tol
+ assert abs(min(mesh.points[:, 1]) + radius) < tol
+ assert abs(max(mesh.points[:, 2]) - radius) < tol
+ assert abs(min(mesh.points[:, 2]) + radius) < tol
+
+ areas = compute_triangle_areas(mesh.points, mesh.cells["triangle"])
+ surface_area = sum(areas)
+ assert abs(surface_area - 4 * numpy.pi * radius ** 2) < 0.1
+ return
+
+
+def test_halfspace():
+ c = pygalmesh.Cuboid([0, 0, 0], [1, 1, 1])
+ s = pygalmesh.HalfSpace([1.0, 2.0, 3.0], 1.0, 2.0)
+ u = pygalmesh.Intersection([c, s])
+
+ pygalmesh.generate_mesh(u, "out.mesh", cell_size=0.2, edge_size=0.2, verbose=False)
+
+ mesh = meshio.read("out.mesh")
+
+ vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"]))
+ assert abs(vol - 1 / 750) < 1.0e-3
+ return
+
+
+if __name__ == "__main__":
+ test_sphere()
projects / pygalmesh.git / commitdiff