From: Drew Parsons
Date: Thu, 15 Oct 2020 13:05:40 +0000 (+0800)
Subject: Import pygalmesh_0.9.1.orig.tar.xz
X-Git-Tag: archive/raspbian/0.10.6-5+rpi1~1^2^2~3
X-Git-Url: https://dgit.raspbian.org/?a=commitdiff_plain;h=9ddbacb8f14e34d037d5ce9c124607ff3e6737c3;p=pygalmesh.git
Import pygalmesh_0.9.1.orig.tar.xz
[dgit import orig pygalmesh_0.9.1.orig.tar.xz]
---
9ddbacb8f14e34d037d5ce9c124607ff3e6737c3
diff --git a/.codecov.yml b/.codecov.yml
new file mode 100644
index 0000000..a052f98
--- /dev/null
+++ b/.codecov.yml
@@ -0,0 +1 @@
+comment: no
diff --git a/.flake8 b/.flake8
new file mode 100644
index 0000000..c321e71
--- /dev/null
+++ b/.flake8
@@ -0,0 +1,5 @@
+[flake8]
+ignore = E203, E266, E501, W503
+max-line-length = 80
+max-complexity = 18
+select = B,C,E,F,W,T4,B9
diff --git a/.gitattributes b/.gitattributes
new file mode 100644
index 0000000..2c8f9bb
--- /dev/null
+++ b/.gitattributes
@@ -0,0 +1,3 @@
+*.inr filter=lfs diff=lfs merge=lfs -text
+*.off filter=lfs diff=lfs merge=lfs -text
+*.vtu filter=lfs diff=lfs merge=lfs -text
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
new file mode 100644
index 0000000..a9460bc
--- /dev/null
+++ b/.github/workflows/ci.yml
@@ -0,0 +1,51 @@
+name: ci
+
+on:
+ push:
+ branches:
+ - master
+ pull_request:
+ branches:
+ - master
+
+jobs:
+ lint:
+ runs-on: ubuntu-latest
+ steps:
+ - uses: actions/setup-python@v2
+ with:
+ python-version: "3.x"
+ - uses: actions/checkout@v2
+ - name: Lint with flake8
+ run: |
+ pip install flake8
+ flake8 .
+ - name: Lint with black
+ run: |
+ pip install black
+ black --check .
+
+ build:
+ runs-on: ubuntu-20.04
+ steps:
+ - uses: actions/setup-python@v2
+ with:
+ python-version: "3.x"
+ - uses: actions/checkout@v2
+ with:
+ lfs: true
+ - name: Install CGAL 5 from PPA
+ run: |
+ sudo apt-get install -y software-properties-common
+ sudo apt-add-repository -y ppa:nschloe/cgal-backports
+ sudo apt update
+ sudo apt install -y libcgal-dev
+ - name: Install other dependencies
+ run: |
+ sudo apt-get install -y libeigen3-dev python3-pip clang
+ - name: Test with tox
+ run: |
+ pip install tox
+ CC="clang" tox
+ - name: Submit to codecov
+ run: bash <(curl -s https://codecov.io/bash)
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..696f5a3
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,9 @@
+*.mesh
+.cache/
+build/
+dist/
+MANIFEST
+README.rst
+do-configure.sh
+.pytest_cache/
+*.egg-info/
diff --git a/.pylintrc b/.pylintrc
new file mode 100644
index 0000000..31f545d
--- /dev/null
+++ b/.pylintrc
@@ -0,0 +1,8 @@
+[MESSAGES CONTROL]
+
+disable=
+ fixme,
+ invalid-name,
+ missing-docstring,
+ no-member,
+ too-few-public-methods
diff --git a/CMakeLists.txt b/CMakeLists.txt
new file mode 100644
index 0000000..56b6e55
--- /dev/null
+++ b/CMakeLists.txt
@@ -0,0 +1,7 @@
+# CMake is used for debugging in pygalmesh. Like every other Python package, the
+# production build system is setuptools.
+cmake_minimum_required(VERSION 3.0)
+
+project(pygalmesh CXX)
+
+add_subdirectory(src)
diff --git a/LICENSE b/LICENSE
new file mode 100644
index 0000000..94a9ed0
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,674 @@
+ GNU GENERAL PUBLIC LICENSE
+ Version 3, 29 June 2007
+
+ Copyright (C) 2007 Free Software Foundation, Inc.
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+ Preamble
+
+ The GNU General Public License is a free, copyleft license for
+software and other kinds of works.
+
+ The licenses for most software and other practical works are designed
+to take away your freedom to share and change the works. By contrast,
+the GNU General Public License is intended to guarantee your freedom to
+share and change all versions of a program--to make sure it remains free
+software for all its users. We, the Free Software Foundation, use the
+GNU General Public License for most of our software; it applies also to
+any other work released this way by its authors. You can apply it to
+your programs, too.
+
+ When we speak of free software, we are referring to freedom, not
+price. Our General Public Licenses are designed to make sure that you
+have the freedom to distribute copies of free software (and charge for
+them if you wish), that you receive source code or can get it if you
+want it, that you can change the software or use pieces of it in new
+free programs, and that you know you can do these things.
+
+ To protect your rights, we need to prevent others from denying you
+these rights or asking you to surrender the rights. Therefore, you have
+certain responsibilities if you distribute copies of the software, or if
+you modify it: responsibilities to respect the freedom of others.
+
+ For example, if you distribute copies of such a program, whether
+gratis or for a fee, you must pass on to the recipients the same
+freedoms that you received. You must make sure that they, too, receive
+or can get the source code. And you must show them these terms so they
+know their rights.
+
+ Developers that use the GNU GPL protect your rights with two steps:
+(1) assert copyright on the software, and (2) offer you this License
+giving you legal permission to copy, distribute and/or modify it.
+
+ For the developers' and authors' protection, the GPL clearly explains
+that there is no warranty for this free software. For both users' and
+authors' sake, the GPL requires that modified versions be marked as
+changed, so that their problems will not be attributed erroneously to
+authors of previous versions.
+
+ Some devices are designed to deny users access to install or run
+modified versions of the software inside them, although the manufacturer
+can do so. This is fundamentally incompatible with the aim of
+protecting users' freedom to change the software. The systematic
+pattern of such abuse occurs in the area of products for individuals to
+use, which is precisely where it is most unacceptable. Therefore, we
+have designed this version of the GPL to prohibit the practice for those
+products. If such problems arise substantially in other domains, we
+stand ready to extend this provision to those domains in future versions
+of the GPL, as needed to protect the freedom of users.
+
+ Finally, every program is threatened constantly by software patents.
+States should not allow patents to restrict development and use of
+software on general-purpose computers, but in those that do, we wish to
+avoid the special danger that patents applied to a free program could
+make it effectively proprietary. To prevent this, the GPL assures that
+patents cannot be used to render the program non-free.
+
+ The precise terms and conditions for copying, distribution and
+modification follow.
+
+ TERMS AND CONDITIONS
+
+ 0. Definitions.
+
+ "This License" refers to version 3 of the GNU General Public License.
+
+ "Copyright" also means copyright-like laws that apply to other kinds of
+works, such as semiconductor masks.
+
+ "The Program" refers to any copyrightable work licensed under this
+License. Each licensee is addressed as "you". "Licensees" and
+"recipients" may be individuals or organizations.
+
+ To "modify" a work means to copy from or adapt all or part of the work
+in a fashion requiring copyright permission, other than the making of an
+exact copy. The resulting work is called a "modified version" of the
+earlier work or a work "based on" the earlier work.
+
+ A "covered work" means either the unmodified Program or a work based
+on the Program.
+
+ To "propagate" a work means to do anything with it that, without
+permission, would make you directly or secondarily liable for
+infringement under applicable copyright law, except executing it on a
+computer or modifying a private copy. Propagation includes copying,
+distribution (with or without modification), making available to the
+public, and in some countries other activities as well.
+
+ To "convey" a work means any kind of propagation that enables other
+parties to make or receive copies. Mere interaction with a user through
+a computer network, with no transfer of a copy, is not conveying.
+
+ An interactive user interface displays "Appropriate Legal Notices"
+to the extent that it includes a convenient and prominently visible
+feature that (1) displays an appropriate copyright notice, and (2)
+tells the user that there is no warranty for the work (except to the
+extent that warranties are provided), that licensees may convey the
+work under this License, and how to view a copy of this License. If
+the interface presents a list of user commands or options, such as a
+menu, a prominent item in the list meets this criterion.
+
+ 1. Source Code.
+
+ The "source code" for a work means the preferred form of the work
+for making modifications to it. "Object code" means any non-source
+form of a work.
+
+ A "Standard Interface" means an interface that either is an official
+standard defined by a recognized standards body, or, in the case of
+interfaces specified for a particular programming language, one that
+is widely used among developers working in that language.
+
+ The "System Libraries" of an executable work include anything, other
+than the work as a whole, that (a) is included in the normal form of
+packaging a Major Component, but which is not part of that Major
+Component, and (b) serves only to enable use of the work with that
+Major Component, or to implement a Standard Interface for which an
+implementation is available to the public in source code form. A
+"Major Component", in this context, means a major essential component
+(kernel, window system, and so on) of the specific operating system
+(if any) on which the executable work runs, or a compiler used to
+produce the work, or an object code interpreter used to run it.
+
+ The "Corresponding Source" for a work in object code form means all
+the source code needed to generate, install, and (for an executable
+work) run the object code and to modify the work, including scripts to
+control those activities. However, it does not include the work's
+System Libraries, or general-purpose tools or generally available free
+programs which are used unmodified in performing those activities but
+which are not part of the work. For example, Corresponding Source
+includes interface definition files associated with source files for
+the work, and the source code for shared libraries and dynamically
+linked subprograms that the work is specifically designed to require,
+such as by intimate data communication or control flow between those
+subprograms and other parts of the work.
+
+ The Corresponding Source need not include anything that users
+can regenerate automatically from other parts of the Corresponding
+Source.
+
+ The Corresponding Source for a work in source code form is that
+same work.
+
+ 2. Basic Permissions.
+
+ All rights granted under this License are granted for the term of
+copyright on the Program, and are irrevocable provided the stated
+conditions are met. This License explicitly affirms your unlimited
+permission to run the unmodified Program. The output from running a
+covered work is covered by this License only if the output, given its
+content, constitutes a covered work. This License acknowledges your
+rights of fair use or other equivalent, as provided by copyright law.
+
+ You may make, run and propagate covered works that you do not
+convey, without conditions so long as your license otherwise remains
+in force. You may convey covered works to others for the sole purpose
+of having them make modifications exclusively for you, or provide you
+with facilities for running those works, provided that you comply with
+the terms of this License in conveying all material for which you do
+not control copyright. Those thus making or running the covered works
+for you must do so exclusively on your behalf, under your direction
+and control, on terms that prohibit them from making any copies of
+your copyrighted material outside their relationship with you.
+
+ Conveying under any other circumstances is permitted solely under
+the conditions stated below. Sublicensing is not allowed; section 10
+makes it unnecessary.
+
+ 3. Protecting Users' Legal Rights From Anti-Circumvention Law.
+
+ No covered work shall be deemed part of an effective technological
+measure under any applicable law fulfilling obligations under article
+11 of the WIPO copyright treaty adopted on 20 December 1996, or
+similar laws prohibiting or restricting circumvention of such
+measures.
+
+ When you convey a covered work, you waive any legal power to forbid
+circumvention of technological measures to the extent such circumvention
+is effected by exercising rights under this License with respect to
+the covered work, and you disclaim any intention to limit operation or
+modification of the work as a means of enforcing, against the work's
+users, your or third parties' legal rights to forbid circumvention of
+technological measures.
+
+ 4. Conveying Verbatim Copies.
+
+ You may convey verbatim copies of the Program's source code as you
+receive it, in any medium, provided that you conspicuously and
+appropriately publish on each copy an appropriate copyright notice;
+keep intact all notices stating that this License and any
+non-permissive terms added in accord with section 7 apply to the code;
+keep intact all notices of the absence of any warranty; and give all
+recipients a copy of this License along with the Program.
+
+ You may charge any price or no price for each copy that you convey,
+and you may offer support or warranty protection for a fee.
+
+ 5. Conveying Modified Source Versions.
+
+ You may convey a work based on the Program, or the modifications to
+produce it from the Program, in the form of source code under the
+terms of section 4, provided that you also meet all of these conditions:
+
+ a) The work must carry prominent notices stating that you modified
+ it, and giving a relevant date.
+
+ b) The work must carry prominent notices stating that it is
+ released under this License and any conditions added under section
+ 7. This requirement modifies the requirement in section 4 to
+ "keep intact all notices".
+
+ c) You must license the entire work, as a whole, under this
+ License to anyone who comes into possession of a copy. This
+ License will therefore apply, along with any applicable section 7
+ additional terms, to the whole of the work, and all its parts,
+ regardless of how they are packaged. This License gives no
+ permission to license the work in any other way, but it does not
+ invalidate such permission if you have separately received it.
+
+ d) If the work has interactive user interfaces, each must display
+ Appropriate Legal Notices; however, if the Program has interactive
+ interfaces that do not display Appropriate Legal Notices, your
+ work need not make them do so.
+
+ A compilation of a covered work with other separate and independent
+works, which are not by their nature extensions of the covered work,
+and which are not combined with it such as to form a larger program,
+in or on a volume of a storage or distribution medium, is called an
+"aggregate" if the compilation and its resulting copyright are not
+used to limit the access or legal rights of the compilation's users
+beyond what the individual works permit. Inclusion of a covered work
+in an aggregate does not cause this License to apply to the other
+parts of the aggregate.
+
+ 6. Conveying Non-Source Forms.
+
+ You may convey a covered work in object code form under the terms
+of sections 4 and 5, provided that you also convey the
+machine-readable Corresponding Source under the terms of this License,
+in one of these ways:
+
+ a) Convey the object code in, or embodied in, a physical product
+ (including a physical distribution medium), accompanied by the
+ Corresponding Source fixed on a durable physical medium
+ customarily used for software interchange.
+
+ b) Convey the object code in, or embodied in, a physical product
+ (including a physical distribution medium), accompanied by a
+ written offer, valid for at least three years and valid for as
+ long as you offer spare parts or customer support for that product
+ model, to give anyone who possesses the object code either (1) a
+ copy of the Corresponding Source for all the software in the
+ product that is covered by this License, on a durable physical
+ medium customarily used for software interchange, for a price no
+ more than your reasonable cost of physically performing this
+ conveying of source, or (2) access to copy the
+ Corresponding Source from a network server at no charge.
+
+ c) Convey individual copies of the object code with a copy of the
+ written offer to provide the Corresponding Source. This
+ alternative is allowed only occasionally and noncommercially, and
+ only if you received the object code with such an offer, in accord
+ with subsection 6b.
+
+ d) Convey the object code by offering access from a designated
+ place (gratis or for a charge), and offer equivalent access to the
+ Corresponding Source in the same way through the same place at no
+ further charge. You need not require recipients to copy the
+ Corresponding Source along with the object code. If the place to
+ copy the object code is a network server, the Corresponding Source
+ may be on a different server (operated by you or a third party)
+ that supports equivalent copying facilities, provided you maintain
+ clear directions next to the object code saying where to find the
+ Corresponding Source. Regardless of what server hosts the
+ Corresponding Source, you remain obligated to ensure that it is
+ available for as long as needed to satisfy these requirements.
+
+ e) Convey the object code using peer-to-peer transmission, provided
+ you inform other peers where the object code and Corresponding
+ Source of the work are being offered to the general public at no
+ charge under subsection 6d.
+
+ A separable portion of the object code, whose source code is excluded
+from the Corresponding Source as a System Library, need not be
+included in conveying the object code work.
+
+ A "User Product" is either (1) a "consumer product", which means any
+tangible personal property which is normally used for personal, family,
+or household purposes, or (2) anything designed or sold for incorporation
+into a dwelling. In determining whether a product is a consumer product,
+doubtful cases shall be resolved in favor of coverage. For a particular
+product received by a particular user, "normally used" refers to a
+typical or common use of that class of product, regardless of the status
+of the particular user or of the way in which the particular user
+actually uses, or expects or is expected to use, the product. A product
+is a consumer product regardless of whether the product has substantial
+commercial, industrial or non-consumer uses, unless such uses represent
+the only significant mode of use of the product.
+
+ "Installation Information" for a User Product means any methods,
+procedures, authorization keys, or other information required to install
+and execute modified versions of a covered work in that User Product from
+a modified version of its Corresponding Source. The information must
+suffice to ensure that the continued functioning of the modified object
+code is in no case prevented or interfered with solely because
+modification has been made.
+
+ If you convey an object code work under this section in, or with, or
+specifically for use in, a User Product, and the conveying occurs as
+part of a transaction in which the right of possession and use of the
+User Product is transferred to the recipient in perpetuity or for a
+fixed term (regardless of how the transaction is characterized), the
+Corresponding Source conveyed under this section must be accompanied
+by the Installation Information. But this requirement does not apply
+if neither you nor any third party retains the ability to install
+modified object code on the User Product (for example, the work has
+been installed in ROM).
+
+ The requirement to provide Installation Information does not include a
+requirement to continue to provide support service, warranty, or updates
+for a work that has been modified or installed by the recipient, or for
+the User Product in which it has been modified or installed. Access to a
+network may be denied when the modification itself materially and
+adversely affects the operation of the network or violates the rules and
+protocols for communication across the network.
+
+ Corresponding Source conveyed, and Installation Information provided,
+in accord with this section must be in a format that is publicly
+documented (and with an implementation available to the public in
+source code form), and must require no special password or key for
+unpacking, reading or copying.
+
+ 7. Additional Terms.
+
+ "Additional permissions" are terms that supplement the terms of this
+License by making exceptions from one or more of its conditions.
+Additional permissions that are applicable to the entire Program shall
+be treated as though they were included in this License, to the extent
+that they are valid under applicable law. If additional permissions
+apply only to part of the Program, that part may be used separately
+under those permissions, but the entire Program remains governed by
+this License without regard to the additional permissions.
+
+ When you convey a copy of a covered work, you may at your option
+remove any additional permissions from that copy, or from any part of
+it. (Additional permissions may be written to require their own
+removal in certain cases when you modify the work.) You may place
+additional permissions on material, added by you to a covered work,
+for which you have or can give appropriate copyright permission.
+
+ Notwithstanding any other provision of this License, for material you
+add to a covered work, you may (if authorized by the copyright holders of
+that material) supplement the terms of this License with terms:
+
+ a) Disclaiming warranty or limiting liability differently from the
+ terms of sections 15 and 16 of this License; or
+
+ b) Requiring preservation of specified reasonable legal notices or
+ author attributions in that material or in the Appropriate Legal
+ Notices displayed by works containing it; or
+
+ c) Prohibiting misrepresentation of the origin of that material, or
+ requiring that modified versions of such material be marked in
+ reasonable ways as different from the original version; or
+
+ d) Limiting the use for publicity purposes of names of licensors or
+ authors of the material; or
+
+ e) Declining to grant rights under trademark law for use of some
+ trade names, trademarks, or service marks; or
+
+ f) Requiring indemnification of licensors and authors of that
+ material by anyone who conveys the material (or modified versions of
+ it) with contractual assumptions of liability to the recipient, for
+ any liability that these contractual assumptions directly impose on
+ those licensors and authors.
+
+ All other non-permissive additional terms are considered "further
+restrictions" within the meaning of section 10. If the Program as you
+received it, or any part of it, contains a notice stating that it is
+governed by this License along with a term that is a further
+restriction, you may remove that term. If a license document contains
+a further restriction but permits relicensing or conveying under this
+License, you may add to a covered work material governed by the terms
+of that license document, provided that the further restriction does
+not survive such relicensing or conveying.
+
+ If you add terms to a covered work in accord with this section, you
+must place, in the relevant source files, a statement of the
+additional terms that apply to those files, or a notice indicating
+where to find the applicable terms.
+
+ Additional terms, permissive or non-permissive, may be stated in the
+form of a separately written license, or stated as exceptions;
+the above requirements apply either way.
+
+ 8. Termination.
+
+ You may not propagate or modify a covered work except as expressly
+provided under this License. Any attempt otherwise to propagate or
+modify it is void, and will automatically terminate your rights under
+this License (including any patent licenses granted under the third
+paragraph of section 11).
+
+ However, if you cease all violation of this License, then your
+license from a particular copyright holder is reinstated (a)
+provisionally, unless and until the copyright holder explicitly and
+finally terminates your license, and (b) permanently, if the copyright
+holder fails to notify you of the violation by some reasonable means
+prior to 60 days after the cessation.
+
+ Moreover, your license from a particular copyright holder is
+reinstated permanently if the copyright holder notifies you of the
+violation by some reasonable means, this is the first time you have
+received notice of violation of this License (for any work) from that
+copyright holder, and you cure the violation prior to 30 days after
+your receipt of the notice.
+
+ Termination of your rights under this section does not terminate the
+licenses of parties who have received copies or rights from you under
+this License. If your rights have been terminated and not permanently
+reinstated, you do not qualify to receive new licenses for the same
+material under section 10.
+
+ 9. Acceptance Not Required for Having Copies.
+
+ You are not required to accept this License in order to receive or
+run a copy of the Program. Ancillary propagation of a covered work
+occurring solely as a consequence of using peer-to-peer transmission
+to receive a copy likewise does not require acceptance. However,
+nothing other than this License grants you permission to propagate or
+modify any covered work. These actions infringe copyright if you do
+not accept this License. Therefore, by modifying or propagating a
+covered work, you indicate your acceptance of this License to do so.
+
+ 10. Automatic Licensing of Downstream Recipients.
+
+ Each time you convey a covered work, the recipient automatically
+receives a license from the original licensors, to run, modify and
+propagate that work, subject to this License. You are not responsible
+for enforcing compliance by third parties with this License.
+
+ An "entity transaction" is a transaction transferring control of an
+organization, or substantially all assets of one, or subdividing an
+organization, or merging organizations. If propagation of a covered
+work results from an entity transaction, each party to that
+transaction who receives a copy of the work also receives whatever
+licenses to the work the party's predecessor in interest had or could
+give under the previous paragraph, plus a right to possession of the
+Corresponding Source of the work from the predecessor in interest, if
+the predecessor has it or can get it with reasonable efforts.
+
+ You may not impose any further restrictions on the exercise of the
+rights granted or affirmed under this License. For example, you may
+not impose a license fee, royalty, or other charge for exercise of
+rights granted under this License, and you may not initiate litigation
+(including a cross-claim or counterclaim in a lawsuit) alleging that
+any patent claim is infringed by making, using, selling, offering for
+sale, or importing the Program or any portion of it.
+
+ 11. Patents.
+
+ A "contributor" is a copyright holder who authorizes use under this
+License of the Program or a work on which the Program is based. The
+work thus licensed is called the contributor's "contributor version".
+
+ A contributor's "essential patent claims" are all patent claims
+owned or controlled by the contributor, whether already acquired or
+hereafter acquired, that would be infringed by some manner, permitted
+by this License, of making, using, or selling its contributor version,
+but do not include claims that would be infringed only as a
+consequence of further modification of the contributor version. For
+purposes of this definition, "control" includes the right to grant
+patent sublicenses in a manner consistent with the requirements of
+this License.
+
+ Each contributor grants you a non-exclusive, worldwide, royalty-free
+patent license under the contributor's essential patent claims, to
+make, use, sell, offer for sale, import and otherwise run, modify and
+propagate the contents of its contributor version.
+
+ In the following three paragraphs, a "patent license" is any express
+agreement or commitment, however denominated, not to enforce a patent
+(such as an express permission to practice a patent or covenant not to
+sue for patent infringement). To "grant" such a patent license to a
+party means to make such an agreement or commitment not to enforce a
+patent against the party.
+
+ If you convey a covered work, knowingly relying on a patent license,
+and the Corresponding Source of the work is not available for anyone
+to copy, free of charge and under the terms of this License, through a
+publicly available network server or other readily accessible means,
+then you must either (1) cause the Corresponding Source to be so
+available, or (2) arrange to deprive yourself of the benefit of the
+patent license for this particular work, or (3) arrange, in a manner
+consistent with the requirements of this License, to extend the patent
+license to downstream recipients. "Knowingly relying" means you have
+actual knowledge that, but for the patent license, your conveying the
+covered work in a country, or your recipient's use of the covered work
+in a country, would infringe one or more identifiable patents in that
+country that you have reason to believe are valid.
+
+ If, pursuant to or in connection with a single transaction or
+arrangement, you convey, or propagate by procuring conveyance of, a
+covered work, and grant a patent license to some of the parties
+receiving the covered work authorizing them to use, propagate, modify
+or convey a specific copy of the covered work, then the patent license
+you grant is automatically extended to all recipients of the covered
+work and works based on it.
+
+ A patent license is "discriminatory" if it does not include within
+the scope of its coverage, prohibits the exercise of, or is
+conditioned on the non-exercise of one or more of the rights that are
+specifically granted under this License. You may not convey a covered
+work if you are a party to an arrangement with a third party that is
+in the business of distributing software, under which you make payment
+to the third party based on the extent of your activity of conveying
+the work, and under which the third party grants, to any of the
+parties who would receive the covered work from you, a discriminatory
+patent license (a) in connection with copies of the covered work
+conveyed by you (or copies made from those copies), or (b) primarily
+for and in connection with specific products or compilations that
+contain the covered work, unless you entered into that arrangement,
+or that patent license was granted, prior to 28 March 2007.
+
+ Nothing in this License shall be construed as excluding or limiting
+any implied license or other defenses to infringement that may
+otherwise be available to you under applicable patent law.
+
+ 12. No Surrender of Others' Freedom.
+
+ If conditions are imposed on you (whether by court order, agreement or
+otherwise) that contradict the conditions of this License, they do not
+excuse you from the conditions of this License. If you cannot convey a
+covered work so as to satisfy simultaneously your obligations under this
+License and any other pertinent obligations, then as a consequence you may
+not convey it at all. For example, if you agree to terms that obligate you
+to collect a royalty for further conveying from those to whom you convey
+the Program, the only way you could satisfy both those terms and this
+License would be to refrain entirely from conveying the Program.
+
+ 13. Use with the GNU Affero General Public License.
+
+ Notwithstanding any other provision of this License, you have
+permission to link or combine any covered work with a work licensed
+under version 3 of the GNU Affero General Public License into a single
+combined work, and to convey the resulting work. The terms of this
+License will continue to apply to the part which is the covered work,
+but the special requirements of the GNU Affero General Public License,
+section 13, concerning interaction through a network will apply to the
+combination as such.
+
+ 14. Revised Versions of this License.
+
+ The Free Software Foundation may publish revised and/or new versions of
+the GNU General Public License from time to time. Such new versions will
+be similar in spirit to the present version, but may differ in detail to
+address new problems or concerns.
+
+ Each version is given a distinguishing version number. If the
+Program specifies that a certain numbered version of the GNU General
+Public License "or any later version" applies to it, you have the
+option of following the terms and conditions either of that numbered
+version or of any later version published by the Free Software
+Foundation. If the Program does not specify a version number of the
+GNU General Public License, you may choose any version ever published
+by the Free Software Foundation.
+
+ If the Program specifies that a proxy can decide which future
+versions of the GNU General Public License can be used, that proxy's
+public statement of acceptance of a version permanently authorizes you
+to choose that version for the Program.
+
+ Later license versions may give you additional or different
+permissions. However, no additional obligations are imposed on any
+author or copyright holder as a result of your choosing to follow a
+later version.
+
+ 15. Disclaimer of Warranty.
+
+ THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY
+APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT
+HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY
+OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO,
+THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM
+IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF
+ALL NECESSARY SERVICING, REPAIR OR CORRECTION.
+
+ 16. Limitation of Liability.
+
+ IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS
+THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY
+GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE
+USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF
+DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD
+PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS),
+EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF
+SUCH DAMAGES.
+
+ 17. Interpretation of Sections 15 and 16.
+
+ If the disclaimer of warranty and limitation of liability provided
+above cannot be given local legal effect according to their terms,
+reviewing courts shall apply local law that most closely approximates
+an absolute waiver of all civil liability in connection with the
+Program, unless a warranty or assumption of liability accompanies a
+copy of the Program in return for a fee.
+
+ END OF TERMS AND CONDITIONS
+
+ How to Apply These Terms to Your New Programs
+
+ If you develop a new program, and you want it to be of the greatest
+possible use to the public, the best way to achieve this is to make it
+free software which everyone can redistribute and change under these terms.
+
+ To do so, attach the following notices to the program. It is safest
+to attach them to the start of each source file to most effectively
+state the exclusion of warranty; and each file should have at least
+the "copyright" line and a pointer to where the full notice is found.
+
+
+ Copyright (C)
+
+ This program is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see .
+
+Also add information on how to contact you by electronic and paper mail.
+
+ If the program does terminal interaction, make it output a short
+notice like this when it starts in an interactive mode:
+
+ Copyright (C)
+ This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
+ This is free software, and you are welcome to redistribute it
+ under certain conditions; type `show c' for details.
+
+The hypothetical commands `show w' and `show c' should show the appropriate
+parts of the General Public License. Of course, your program's commands
+might be different; for a GUI interface, you would use an "about box".
+
+ You should also get your employer (if you work as a programmer) or school,
+if any, to sign a "copyright disclaimer" for the program, if necessary.
+For more information on this, and how to apply and follow the GNU GPL, see
+.
+
+ The GNU General Public License does not permit incorporating your program
+into proprietary programs. If your program is a subroutine library, you
+may consider it more useful to permit linking proprietary applications with
+the library. If this is what you want to do, use the GNU Lesser General
+Public License instead of this License. But first, please read
+.
diff --git a/MANIFEST.in b/MANIFEST.in
new file mode 100644
index 0000000..0163cb5
--- /dev/null
+++ b/MANIFEST.in
@@ -0,0 +1,14 @@
+include src/domain.hpp
+include src/generate_from_inr.hpp
+include src/generate_from_off.hpp
+include src/generate_periodic.hpp
+include src/generate.hpp
+include src/generate_2d.hpp
+include src/generate_surface_mesh.hpp
+include src/polygon2d.hpp
+include src/primitives.hpp
+include src/remesh_surface.hpp
+include src/sizing_field.hpp
+
+include test/helpers.py
+recursive-include test/meshes *
diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..79d6e0d
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,34 @@
+VERSION=$(shell python3 -c "from configparser import ConfigParser; p = ConfigParser(); p.read('setup.cfg'); print(p['metadata']['version'])")
+
+default:
+ @echo "\"make publish\"?"
+
+tag:
+ @if [ "$(shell git rev-parse --abbrev-ref HEAD)" != "master" ]; then exit 1; fi
+ curl -H "Authorization: token `cat $(HOME)/.github-access-token`" -d '{"tag_name": "v$(VERSION)"}' https://api.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
+ 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
+ # twine upload dist/*.whl
+
+publish: tag upload
+
+clean:
+ @find . | grep -E "(__pycache__|\.pyc|\.pyo$\)" | xargs rm -rf
+ @rm -rf build/*
+ @rm -rf pygalmesh.egg-info/
+ @rm -rf dist/
+
+format:
+ isort .
+ black .
+ blacken-docs README.md
+
+lint:
+ black --check .
+ flake8 setup.py pygalmesh/ test/*.py
diff --git a/README.md b/README.md
new file mode 100644
index 0000000..18bf50c
--- /dev/null
+++ b/README.md
@@ -0,0 +1,495 @@
+
+
+
Create high-quality meshes with ease.
+
+
+[](https://pypi.org/project/pygalmesh)
+[](https://anaconda.org/conda-forge/pygalmesh/)
+[](https://repology.org/project/pygalmesh/versions)
+[](https://pypi.org/pypi/pygalmesh/)
+[](https://github.com/nschloe/pygalmesh)
+[](https://pypistats.org/packages/pygalmesh)
+
+[](https://join.slack.com/t/nschloe/shared_invite/zt-cofhrwm8-BgdrXAtVkOjnDmADROKD7A
+)
+
+[](https://github.com/nschloe/pygalmesh/actions?query=workflow%3Aci)
+[](https://codecov.io/gh/nschloe/pygalmesh)
+[](https://lgtm.com/projects/g/nschloe/pygalmesh)
+[](https://github.com/psf/black)
+
+pygalmesh is a Python frontend to [CGAL](https://www.cgal.org/)'s
+[2D](https://doc.cgal.org/latest/Mesh_2/index.html) and [3D mesh generation
+capabilities](https://doc.cgal.org/latest/Mesh_3/index.html). pygalmesh makes it easy
+to create high-quality 2D, 3D volume meshes, periodic volume meshes, and surface meshes.
+
+### Examples
+
+#### 2D meshes
+
+
+CGAL generates 2D meshes from linear contraints.
+```python
+import numpy
+import pygalmesh
+
+points = numpy.array([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])
+constraints = [[0, 1], [1, 2], [2, 3], [3, 0]]
+
+mesh = pygalmesh.generate_2d(
+ points,
+ constraints,
+ max_edge_size=1.0e-1,
+ num_lloyd_steps=10,
+)
+# mesh.points, mesh.cells
+```
+The quality of the mesh isn't very good, but can be improved with
+[optimesh](https://github.com/nschloe/optimesh).
+
+#### A simple ball
+
+
+```python
+import pygalmesh
+
+s = pygalmesh.Ball([0, 0, 0], 1.0)
+mesh = pygalmesh.generate_mesh(s, max_cell_circumradius=0.2)
+
+# mesh.points, mesh.cells, ...
+```
+You can write the mesh with
+
+```python
+mesh.write("out.vtk")
+```
+You can use any format supported by [meshio](https://github.com/nschloe/meshio).
+
+The mesh generation comes with many more options, described
+[here](https://doc.cgal.org/latest/Mesh_3/). Try, for example,
+
+```python
+mesh = pygalmesh.generate_mesh(
+ s, max_cell_circumradius=0.2, odt=True, lloyd=True, verbose=False
+)
+```
+
+#### Other primitive shapes
+
+
+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]
+)
+mesh = pygalmesh.generate_mesh(
+ s0,
+ max_cell_circumradius=0.1,
+ max_edge_size_at_feature_edges=0.1,
+)
+```
+
+#### Domain combinations
+
+
+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
+import numpy
+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)
+
+# add circle
+a = numpy.sqrt(radius ** 2 - displacement ** 2)
+max_edge_size_at_feature_edges = 0.15
+n = int(2 * numpy.pi * a / max_edge_size_at_feature_edges)
+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])
+
+mesh = pygalmesh.generate_mesh(
+ u,
+ extra_feature_edges=[circ],
+ max_cell_circumradius=0.15,
+ max_edge_size_at_feature_edges=max_edge_size_at_feature_edges,
+ min_facet_angle=25,
+ max_radius_surface_delaunay_ball=0.15,
+ max_circumradius_edge_ratio=2.0,
+)
+```
+Note that the length of the polygon legs are kept in sync with
+`max_edge_size_at_feature_edges` of the mesh generation. This makes sure that it fits in
+nicely with the rest of the mesh.
+
+#### Domain deformations
+
+
+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])
+
+mesh = pygalmesh.generate_mesh(s, max_cell_circumradius=0.1)
+```
+
+#### Extrusion of 2D polygons
+
+
+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]])
+max_edge_size_at_feature_edges = 0.1
+domain = pygalmesh.Extrude(
+ p,
+ [0.0, 0.0, 1.0],
+ 0.5 * 3.14159265359,
+ max_edge_size_at_feature_edges,
+)
+mesh = pygalmesh.generate_mesh(
+ domain,
+ max_cell_circumradius=0.1,
+ max_edge_size_at_feature_edges=max_edge_size_at_feature_edges,
+ verbose=False,
+)
+```
+Feature edges are automatically preserved here, which is why an edge length needs to be
+given to `pygalmesh.Extrude`.
+
+#### Rotation bodies
+
+
+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]])
+max_edge_size_at_feature_edges = 0.1
+domain = pygalmesh.RingExtrude(p, max_edge_size_at_feature_edges)
+mesh = pygalmesh.generate_mesh(
+ domain,
+ max_cell_circumradius=0.1,
+ max_edge_size_at_feature_edges=max_edge_size_at_feature_edges,
+ verbose=False,
+)
+```
+
+#### Your own custom level set function
+
+
+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().__init__()
+
+ 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()
+mesh = pygalmesh.generate_mesh(d, max_cell_circumradius=0.1)
+```
+Note that you need to specify the square of a bounding sphere radius, used as an input
+to CGAL's mesh generator.
+
+
+#### Local refinement
+
+
+Use `generate_mesh` with a function (regular or lambda) as `max_cell_circumradius`. The
+same goes for `max_edge_size_at_feature_edges`, `max_radius_surface_delaunay_ball`, and
+`max_facet_distance`.
+```python
+import numpy
+import pygalmesh
+
+mesh = pygalmesh.generate_mesh(
+ pygalmesh.Ball([0.0, 0.0, 0.0], 1.0),
+ min_facet_angle=30.0,
+ max_radius_surface_delaunay_ball=0.1,
+ max_facet_distance=0.025,
+ max_circumradius_edge_ratio=2.0,
+ max_cell_circumradius=lambda x: abs(numpy.sqrt(numpy.dot(x, x)) - 0.5) / 5 + 0.025,
+)
+```
+
+#### Surface meshes
+
+If you're only after the surface of a body, pygalmesh has `generate_surface_mesh` for
+you. It offers fewer options (obviously, `max_cell_circumradius` is gone), but otherwise
+works the same way:
+```python
+import pygalmesh
+
+s = pygalmesh.Ball([0, 0, 0], 1.0)
+mesh = pygalmesh.generate_surface_mesh(
+ s,
+ min_facet_angle=30.0,
+ max_radius_surface_delaunay_ball=0.1,
+ max_facet_distance=0.1,
+)
+```
+Refer to [CGAL's
+documention](https://doc.cgal.org/latest/Surface_mesher/index.html) for the
+options.
+
+#### Periodic volume meshes
+
+
+pygalmesh also interfaces CGAL's [3D periodic
+mesh generation](https://doc.cgal.org/latest/Periodic_3_mesh_3/index.html). Besides a
+domain, one needs to specify a bounding box, and optionally the number of copies in the
+output (1, 2, 4, or 8). Example:
+```python
+import numpy
+import pygalmesh
+
+
+class Schwarz(pygalmesh.DomainBase):
+ def __init__(self):
+ super().__init__()
+
+ def eval(self, x):
+ x2 = numpy.cos(x[0] * 2 * numpy.pi)
+ y2 = numpy.cos(x[1] * 2 * numpy.pi)
+ z2 = numpy.cos(x[2] * 2 * numpy.pi)
+ return x2 + y2 + z2
+
+
+mesh = pygalmesh.generate_periodic_mesh(
+ Schwarz(),
+ [0, 0, 0, 1, 1, 1],
+ max_cell_circumradius=0.05,
+ min_facet_angle=30,
+ max_radius_surface_delaunay_ball=0.05,
+ max_facet_distance=0.025,
+ max_circumradius_edge_ratio=2.0,
+ number_of_copies_in_output=4,
+ # odt=True,
+ # lloyd=True,
+ verbose=False,
+)
+```
+
+#### Volume meshes from surface meshes
+
+
+If you have a surface mesh at hand (like
+[elephant.vtu](http://nschloe.github.io/pygalmesh/elephant.vtu)), pygalmesh generates a
+volume mesh on the command line via
+```
+pygalmesh-volume-from-surface elephant.vtu out.vtk --cell-size 1.0 --odt
+```
+(See `pygalmesh-volume-from-surface -h` for all options.)
+
+In Python, do
+
+```python
+import pygalmesh
+
+mesh = pygalmesh.generate_volume_mesh_from_surface_mesh(
+ "elephant.vtu",
+ min_facet_angle=25.0,
+ max_radius_surface_delaunay_ball=0.15,
+ max_facet_distance=0.008,
+ max_circumradius_edge_ratio=3.0,
+ verbose=False,
+)
+```
+
+#### Meshes from INR voxel files
+
+
+It is also possible to generate meshes from INR voxel files, e.g.,
+[skull_2.9.inr](https://github.com/nschloe/pygalmesh/raw/gh-pages/skull_2.9.inr)
+either on the command line
+```
+pygalmesh-from-inr skull_2.9.inr out.vtu --cell-size 5.0 --odt
+```
+(see `pygalmesh-from-inr -h` for all options) or from Python
+
+```python
+import pygalmesh
+
+mesh = pygalmesh.generate_from_inr(
+ "skull_2.9.inr",
+ max_cell_circumradius=5.0,
+ verbose=False,
+)
+```
+
+#### Meshes from numpy arrays representing 3D images
+
+
+pygalmesh can help generating unstructed meshes from 3D numpy arrays.
+
+The code below creates a mesh from the 3D breast phantom from [Lou et
+al](http://biomedicaloptics.spiedigitallibrary.org/article.aspx?articleid=2600985)
+available
+[here](https://wustl.app.box.com/s/rqivtin0xcofjwlkz43acou8jknsbfx8/file/127108205145).
+The phantom comprises four tissue types (background, fat, fibrograndular, skin, vascular
+tissues). The generated mesh conforms to tissues interfaces.
+
+```python
+import pygalmesh
+import meshio
+
+Nx = 722
+Ny = 411
+Nz = 284
+h = [0.2] * 3
+
+with open("MergedPhantom.DAT", "rb") as fid:
+ vol = np.fromfile(fid, dtype=np.uint8)
+
+vol = vol.reshape((Nx, Ny, Nz))
+
+mesh = pygalmesh.generate_from_array(
+ vol, h, max_facet_distance=0.2, max_cell_circumradius=1.0
+)
+mesh.write("breast.vtk")
+```
+
+In addition, we can specify different mesh sizes for each tissue type. The code below
+sets the mesh size to *1 mm* for the skin tissue (label `4`), *0.5 mm* for the vascular
+tissue (label `5`), and *2 mm* for all other tissues (`default`).
+
+
+```python
+mesh = pygalmesh.generate_from_array(
+ vol,
+ h,
+ max_facet_distance=0.2,
+ max_cell_circumradius={"default": 2.0, 4: 1.0, 5: 0.5},
+)
+mesh.write("breast_adapted.vtk")
+```
+
+#### Surface remeshing
+
|
+:---:|:---:|
+
+pygalmesh can help remeshing an existing surface mesh, e.g.,
+[`lion-head.off`](https://github.com/nschloe/pygalmesh/raw/gh-pages/lion-head.off). On
+the command line, use
+
+```
+pygalmesh-remesh-surface lion-head.off out.vtu -e 0.025 -a 25 -s 0.1 -d 0.001
+```
+(see `pygalmesh-remesh-surface -h` for all options) or from Python
+
+```python
+import pygalmesh
+
+mesh = pygalmesh.remesh_surface(
+ "lion-head.off",
+ max_edge_size_at_feature_edges=0.025,
+ min_facet_angle=25,
+ max_radius_surface_delaunay_ball=0.1,
+ max_facet_distance=0.001,
+ 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.
+
+#### 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 `python3 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
+```
+
+
+### Background
+
+CGAL offers two different approaches for mesh generation:
+
+1. Meshes defined implicitly by level sets of functions.
+2. Meshes defined by a set of bounding planes.
+
+pygalmesh provides a front-end to the first approach, which has the following advantages
+and disadvantages:
+
+* All boundary points are guaranteed to be in the level set within any specified
+ residual. This results in smooth curved surfaces.
+* Sharp intersections of subdomains (e.g., in unions or differences of sets) need to be
+ specified manually (via feature edges, see below), which can be tedious.
+
+On the other hand, the bounding-plane approach (realized by
+[mshr](https://bitbucket.org/fenics-project/mshr)), has the following properties:
+
+* Smooth, curved domains are approximated by a set of bounding planes, resulting in more
+ of less visible edges.
+* Intersections of domains can be computed automatically, so domain unions etc. have
+ sharp edges where they belong.
+
+See [here](https://github.com/nschloe/awesome-scientific-computing#meshing) for other
+mesh generation tools.
+
+
+### License
+
+pygalmesh is published under the [GPLv3 license](https://www.gnu.org/licenses/gpl-3.0.en.html).
diff --git a/pygalmesh/__about__.py b/pygalmesh/__about__.py
new file mode 100644
index 0000000..8c2ff90
--- /dev/null
+++ b/pygalmesh/__about__.py
@@ -0,0 +1,14 @@
+from _pygalmesh import _CGAL_VERSION_STR
+
+try:
+ # Python 3.8
+ from importlib import metadata
+except ImportError:
+ import importlib_metadata as metadata
+
+try:
+ __version__ = metadata.version("pygalmesh")
+except Exception:
+ __version__ = "unknown"
+
+__cgal_version__ = _CGAL_VERSION_STR
diff --git a/pygalmesh/__init__.py b/pygalmesh/__init__.py
new file mode 100644
index 0000000..3a73598
--- /dev/null
+++ b/pygalmesh/__init__.py
@@ -0,0 +1,73 @@
+# https://github.com/pybind/pybind11/issues/1004
+from _pygalmesh import (
+ Ball,
+ Cone,
+ Cuboid,
+ Cylinder,
+ Difference,
+ DomainBase,
+ Ellipsoid,
+ Extrude,
+ HalfSpace,
+ Intersection,
+ Polygon2D,
+ RingExtrude,
+ Rotate,
+ Scale,
+ Stretch,
+ Tetrahedron,
+ Torus,
+ Translate,
+ Union,
+)
+
+from . import _cli
+from .__about__ import __cgal_version__, __version__
+from .main import (
+ generate_2d,
+ generate_from_array,
+ generate_from_inr,
+ generate_mesh,
+ generate_periodic_mesh,
+ generate_surface_mesh,
+ generate_volume_mesh_from_surface_mesh,
+ remesh_surface,
+ save_inr,
+)
+
+__all__ = [
+ "__version__",
+ "__cgal_version__",
+ #
+ "_cli",
+ #
+ "DomainBase",
+ "Translate",
+ "Rotate",
+ "Scale",
+ "Stretch",
+ "Intersection",
+ "Union",
+ "Difference",
+ "Extrude",
+ "Ball",
+ "Cuboid",
+ "Ellipsoid",
+ "Tetrahedron",
+ "Cone",
+ "Cylinder",
+ "Torus",
+ "HalfSpace",
+ "Polygon2D",
+ "RingExtrude",
+ #
+ "generate_mesh",
+ "generate_2d",
+ "generate_periodic_mesh",
+ "generate_surface_mesh",
+ "generate_volume_mesh_from_surface_mesh",
+ "generate_from_array",
+ "generate_from_inr",
+ "remesh_surface",
+ "save_inr",
+]
diff --git a/pygalmesh/_cli/__init__.py b/pygalmesh/_cli/__init__.py
new file mode 100644
index 0000000..180a34c
--- /dev/null
+++ b/pygalmesh/_cli/__init__.py
@@ -0,0 +1,5 @@
+from ._inr import inr
+from ._remesh_surface import remesh_surface
+from ._volume_from_surface import volume_from_surface
+
+__all__ = ["inr", "volume_from_surface", "remesh_surface"]
diff --git a/pygalmesh/_cli/_inr.py b/pygalmesh/_cli/_inr.py
new file mode 100644
index 0000000..38d48a6
--- /dev/null
+++ b/pygalmesh/_cli/_inr.py
@@ -0,0 +1,131 @@
+import argparse
+
+import meshio
+
+from ..main import generate_from_inr
+from .helpers import _get_version_text
+
+
+def inr(argv=None):
+ parser = _get_inr_parser()
+ args = parser.parse_args(argv)
+
+ mesh = generate_from_inr(
+ args.infile,
+ lloyd=args.lloyd,
+ odt=args.odt,
+ perturb=args.perturb,
+ exude=args.exude,
+ max_edge_size_at_feature_edges=args.max_edge_size_at_feature_edges,
+ min_facet_angle=args.min_facet_angle,
+ max_radius_surface_delaunay_ball=args.max_radius_surface_delaunay_ball,
+ max_facet_distance=args.max_facet_distance,
+ max_circumradius_edge_ratio=args.max_circumradius_edge_ratio,
+ max_cell_circumradius=args.max_cell_circumradius,
+ verbose=not args.quiet,
+ )
+ meshio.write(args.outfile, mesh)
+
+
+def _get_inr_parser():
+ parser = argparse.ArgumentParser(
+ description="Generate volume mesh from INR voxel file.",
+ formatter_class=argparse.RawTextHelpFormatter,
+ )
+
+ parser.add_argument("infile", type=str, help="input INR file")
+
+ parser.add_argument("outfile", type=str, help="output mesh file")
+
+ parser.add_argument(
+ "--lloyd",
+ "-l",
+ type=bool,
+ default=False,
+ help="Lloyd smoothing (default: false)",
+ )
+
+ parser.add_argument(
+ "--odt",
+ "-o",
+ action="store_true",
+ default=False,
+ help="ODT smoothing (default: false)",
+ )
+
+ parser.add_argument(
+ "--perturb",
+ "-p",
+ action="store_true",
+ default=False,
+ help="perturb (default: false)",
+ )
+
+ parser.add_argument(
+ "--exude",
+ "-x",
+ action="store_true",
+ default=False,
+ help="exude (default: false)",
+ )
+
+ parser.add_argument(
+ "--max-edge-size-at-feature-edges",
+ "-e",
+ type=float,
+ default=0.0,
+ help="maximum edge size at feature edges (default: 0.0)",
+ )
+
+ parser.add_argument(
+ "--min-facet-angle",
+ "-a",
+ type=float,
+ default=0.0,
+ help="minimum facet angle (default: 0.0)",
+ )
+
+ parser.add_argument(
+ "--max-radius-surface-delaunay-ball",
+ "-s",
+ type=float,
+ default=0.0,
+ help="maximum radius of the surface facet Delaunay ball (default: 0.0)",
+ )
+
+ parser.add_argument(
+ "--max-facet-distance",
+ "-d",
+ type=float,
+ default=0.0,
+ help="maximum facet distance (default: 0.0)",
+ )
+
+ parser.add_argument(
+ "--max-circumradius-edge-ratio",
+ "-r",
+ type=float,
+ default=0.0,
+ help="cell radius/edge ratio (default: 0.0)",
+ )
+
+ parser.add_argument(
+ "--max-cell-circumradius",
+ "-c",
+ type=float,
+ default=0.0,
+ help="maximum cell circumradius (default: 0.0)",
+ )
+
+ parser.add_argument(
+ "--quiet",
+ "-q",
+ action="store_true",
+ default=False,
+ help="quiet mode (default: False)",
+ )
+
+ parser.add_argument(
+ "--version", "-v", action="version", version=_get_version_text()
+ )
+ return parser
diff --git a/pygalmesh/_cli/_remesh_surface.py b/pygalmesh/_cli/_remesh_surface.py
new file mode 100644
index 0000000..2964497
--- /dev/null
+++ b/pygalmesh/_cli/_remesh_surface.py
@@ -0,0 +1,80 @@
+import argparse
+
+import meshio
+
+from ..main import remesh_surface as rms
+from .helpers import _get_version_text
+
+
+def remesh_surface(argv=None):
+ parser = _get_parser()
+ args = parser.parse_args(argv)
+
+ mesh = rms(
+ args.infile,
+ # lloyd=args.lloyd,
+ # odt=args.odt,
+ # perturb=args.perturb,
+ # exude=args.exude,
+ max_edge_size_at_feature_edges=args.max_edge_size_at_feature_edges,
+ min_facet_angle=args.min_facet_angle,
+ max_radius_surface_delaunay_ball=args.max_radius_surface_delaunay_ball,
+ max_facet_distance=args.max_facet_distance,
+ verbose=not args.quiet,
+ )
+ meshio.write(args.outfile, mesh)
+
+
+def _get_parser():
+ parser = argparse.ArgumentParser(
+ description="Remesh surface mesh.",
+ formatter_class=argparse.RawTextHelpFormatter,
+ )
+
+ parser.add_argument("infile", type=str, help="input mesh file")
+ parser.add_argument("outfile", type=str, help="output mesh file")
+
+ parser.add_argument(
+ "--max-edge-size-at-feature-edges",
+ "-e",
+ type=float,
+ default=0.0,
+ help="maximum edge size at feature edges (default: 0.0)",
+ )
+
+ parser.add_argument(
+ "--min-facet-angle",
+ "-a",
+ type=float,
+ default=0.0,
+ help="minimum facet angle (default: 0.0)",
+ )
+
+ parser.add_argument(
+ "--max-radius-surface-delaunay-ball",
+ "-s",
+ type=float,
+ default=0.0,
+ help="maximum radius of the surface facet Delaunay ball (default: 0.0)",
+ )
+
+ parser.add_argument(
+ "--max-facet-distance",
+ "-d",
+ type=float,
+ default=0.0,
+ help="maximum facet distance (default: 0.0)",
+ )
+
+ parser.add_argument(
+ "--quiet",
+ "-q",
+ action="store_true",
+ default=False,
+ help="quiet mode (default: False)",
+ )
+
+ parser.add_argument(
+ "--version", "-v", action="version", version=_get_version_text()
+ )
+ return parser
diff --git a/pygalmesh/_cli/_volume_from_surface.py b/pygalmesh/_cli/_volume_from_surface.py
new file mode 100644
index 0000000..4216bc8
--- /dev/null
+++ b/pygalmesh/_cli/_volume_from_surface.py
@@ -0,0 +1,140 @@
+import argparse
+
+import meshio
+
+from ..main import generate_volume_mesh_from_surface_mesh
+from .helpers import _get_version_text
+
+
+def volume_from_surface(argv=None):
+ parser = _get_volume_from_surface_parser()
+ args = parser.parse_args(argv)
+
+ mesh = generate_volume_mesh_from_surface_mesh(
+ args.infile,
+ lloyd=args.lloyd,
+ odt=args.odt,
+ perturb=args.perturb,
+ exude=args.exude,
+ max_edge_size_at_feature_edges=args.max_edge_size_at_feature_edges,
+ min_facet_angle=args.min_facet_angle,
+ max_radius_surface_delaunay_ball=args.max_radius_surface_delaunay_ball,
+ max_facet_distance=args.max_facet_distance,
+ max_circumradius_edge_ratio=args.max_circumradius_edge_ratio,
+ max_cell_circumradius=args.max_cell_circumradius,
+ reorient=args.reorient,
+ verbose=not args.quiet,
+ )
+ meshio.write(args.outfile, mesh)
+
+
+def _get_volume_from_surface_parser():
+ parser = argparse.ArgumentParser(
+ description="Generate volume mesh from surface mesh.",
+ formatter_class=argparse.RawTextHelpFormatter,
+ )
+
+ parser.add_argument("infile", type=str, help="input mesh file")
+
+ parser.add_argument("outfile", type=str, help="output mesh file")
+
+ parser.add_argument(
+ "--lloyd",
+ "-l",
+ type=bool,
+ default=False,
+ help="Lloyd smoothing (default: false)",
+ )
+
+ parser.add_argument(
+ "--odt",
+ "-o",
+ action="store_true",
+ default=False,
+ help="ODT smoothing (default: false)",
+ )
+
+ parser.add_argument(
+ "--perturb",
+ "-p",
+ action="store_true",
+ default=False,
+ help="perturb (default: false)",
+ )
+
+ parser.add_argument(
+ "--exude",
+ "-x",
+ action="store_true",
+ default=False,
+ help="exude (default: false)",
+ )
+
+ parser.add_argument(
+ "--max-edge-size-at-feature-edges",
+ "-e",
+ type=float,
+ default=0.0,
+ help="maximum edge size at feature edges (default: 0.0)",
+ )
+
+ parser.add_argument(
+ "--min-facet-angle",
+ "-a",
+ type=float,
+ default=0.0,
+ help="minimum facet angle (default: 0.0)",
+ )
+
+ parser.add_argument(
+ "--max-radius-surface-delaunay-ball",
+ "-s",
+ type=float,
+ default=0.0,
+ help="maximum radius of the surface facet Delaunay ball (default: 0.0)",
+ )
+
+ parser.add_argument(
+ "--max-facet-distance",
+ "-d",
+ type=float,
+ default=0.0,
+ help="maximum facet distance (default: 0.0)",
+ )
+
+ parser.add_argument(
+ "--max-circumradius-edge-ratio",
+ "-r",
+ type=float,
+ default=0.0,
+ help="cell radius/edge ratio (default: 0.0)",
+ )
+
+ parser.add_argument(
+ "--max-cell-circumradius",
+ "-c",
+ type=float,
+ default=0.0,
+ help="maximum cell circumradius (default: 0.0)",
+ )
+
+ parser.add_argument(
+ "--reorient",
+ "-t",
+ action="store_true",
+ default=False,
+ help="automatically fix face orientation (default: False)",
+ )
+
+ parser.add_argument(
+ "--quiet",
+ "-q",
+ action="store_true",
+ default=False,
+ help="quiet mode (default: False)",
+ )
+
+ parser.add_argument(
+ "--version", "-v", action="version", version=_get_version_text()
+ )
+ return parser
diff --git a/pygalmesh/_cli/helpers.py b/pygalmesh/_cli/helpers.py
new file mode 100644
index 0000000..5ab9f20
--- /dev/null
+++ b/pygalmesh/_cli/helpers.py
@@ -0,0 +1,20 @@
+import sys
+
+from _pygalmesh import _CGAL_VERSION_STR
+
+from ..__about__ import __version__
+
+
+def _get_version_text():
+ return "\n".join(
+ [
+ "pygalmesh {} [Python {}.{}.{}, CGAL {}]".format(
+ __version__,
+ sys.version_info.major,
+ sys.version_info.minor,
+ sys.version_info.micro,
+ _CGAL_VERSION_STR,
+ ),
+ "Copyright (c) 2016-2020, Nico Schlömer ",
+ ]
+ )
diff --git a/pygalmesh/main.py b/pygalmesh/main.py
new file mode 100644
index 0000000..a870be6
--- /dev/null
+++ b/pygalmesh/main.py
@@ -0,0 +1,450 @@
+import math
+import os
+import tempfile
+
+import meshio
+import numpy
+from _pygalmesh import (
+ SizingFieldBase,
+ _generate_2d,
+ _generate_from_inr,
+ _generate_from_inr_with_subdomain_sizing,
+ _generate_from_off,
+ _generate_mesh,
+ _generate_periodic_mesh,
+ _generate_surface_mesh,
+ _remesh_surface,
+)
+
+
+class Wrapper(SizingFieldBase):
+ def __init__(self, f):
+ self.f = f
+ super().__init__()
+
+ def eval(self, x):
+ return self.f(x)
+
+
+def generate_mesh(
+ domain,
+ 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,
+):
+ """
+ From :
+
+ max_edge_size_at_feature_edges:
+ a scalar field (resp. a constant) providing a space varying (resp. a uniform)
+ upper bound for the lengths of curve edges. This parameter has to be set to a
+ positive value when 1-dimensional features protection is used.
+ min_facet_angle:
+ a lower bound for the angles (in degrees) of the surface mesh facets.
+ max_radius_surface_delaunay_ball:
+ a scalar field (resp. a constant) describing a space varying (resp. a uniform)
+ upper-bound or for the radii of the surface Delaunay balls.
+ max_facet_distance:
+ a scalar field (resp. a constant) describing a space varying (resp. a uniform)
+ upper bound for the distance between the facet circumcenter and the center of
+ its surface Delaunay ball.
+ max_circumradius_edge_ratio:
+ an upper bound for the radius-edge ratio of the mesh tetrahedra.
+ max_cell_circumradius:
+ a scalar field (resp. a constant) describing a space varying (resp. a uniform)
+ upper-bound for the circumradii of the mesh tetrahedra.
+ """
+ extra_feature_edges = [] if extra_feature_edges is None else extra_feature_edges
+
+ fh, outfile = tempfile.mkstemp(suffix=".mesh")
+ os.close(fh)
+
+ def _select(obj):
+ if isinstance(obj, float):
+ return obj, None
+ assert callable(obj)
+ return -1.0, Wrapper(obj)
+
+ (
+ max_edge_size_at_feature_edges_value,
+ max_edge_size_at_feature_edges_field,
+ ) = _select(max_edge_size_at_feature_edges)
+ max_cell_circumradius_value, max_cell_circumradius_field = _select(
+ max_cell_circumradius
+ )
+ (
+ max_radius_surface_delaunay_ball_value,
+ max_radius_surface_delaunay_ball_field,
+ ) = _select(max_radius_surface_delaunay_ball)
+ max_facet_distance_value, max_facet_distance_field = _select(max_facet_distance)
+
+ # if feature_edges:
+ # if max_edge_size_at_feature_edges == 0.0:
+ # raise ValueError(
+ # "Need a positive max_edge_size_at_feature_edges bound if feature_edges are present."
+ # )
+ # elif max_edge_size_at_feature_edges != 0.0:
+ # warnings.warn(
+ # "No feature edges. The max_edge_size_at_feature_edges argument has no effect."
+ # )
+
+ _generate_mesh(
+ domain,
+ outfile,
+ extra_feature_edges=extra_feature_edges,
+ bounding_sphere_radius=bounding_sphere_radius,
+ lloyd=lloyd,
+ odt=odt,
+ perturb=perturb,
+ exude=exude,
+ max_edge_size_at_feature_edges_value=max_edge_size_at_feature_edges_value,
+ max_edge_size_at_feature_edges_field=max_edge_size_at_feature_edges_field,
+ min_facet_angle=min_facet_angle,
+ max_radius_surface_delaunay_ball_value=max_radius_surface_delaunay_ball_value,
+ max_radius_surface_delaunay_ball_field=max_radius_surface_delaunay_ball_field,
+ max_facet_distance_value=max_facet_distance_value,
+ max_facet_distance_field=max_facet_distance_field,
+ max_circumradius_edge_ratio=max_circumradius_edge_ratio,
+ max_cell_circumradius_value=max_cell_circumradius_value,
+ max_cell_circumradius_field=max_cell_circumradius_field,
+ verbose=verbose,
+ seed=seed,
+ )
+
+ mesh = meshio.read(outfile)
+ os.remove(outfile)
+ return mesh
+
+
+def generate_2d(
+ points,
+ constraints,
+ B=math.sqrt(2),
+ max_edge_size=0.0,
+ num_lloyd_steps=0,
+):
+ # some sanity checks
+ points = numpy.asarray(points)
+ constraints = numpy.asarray(constraints)
+ assert numpy.all(constraints >= 0)
+ assert numpy.all(constraints < len(points))
+ # make sure there are no edges of 0 length
+ edges = points[constraints[:, 0]] - points[constraints[:, 1]]
+ length2 = numpy.einsum("ij,ij->i", edges, edges)
+ if numpy.any(length2 < 1.0e-15):
+ raise RuntimeError("Constraint of (near)-zero length.")
+
+ points, cells = _generate_2d(
+ points,
+ constraints,
+ B,
+ max_edge_size,
+ num_lloyd_steps,
+ )
+ return meshio.Mesh(numpy.array(points), {"triangle": numpy.array(cells)})
+
+
+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,
+):
+ fh, outfile = tempfile.mkstemp(suffix=".mesh")
+ os.close(fh)
+
+ assert number_of_copies_in_output in [1, 2, 4, 8]
+
+ _generate_periodic_mesh(
+ domain,
+ outfile,
+ bounding_cuboid,
+ lloyd=lloyd,
+ odt=odt,
+ perturb=perturb,
+ exude=exude,
+ max_edge_size_at_feature_edges=max_edge_size_at_feature_edges,
+ min_facet_angle=min_facet_angle,
+ max_radius_surface_delaunay_ball=max_radius_surface_delaunay_ball,
+ max_facet_distance=max_facet_distance,
+ max_circumradius_edge_ratio=max_circumradius_edge_ratio,
+ max_cell_circumradius=max_cell_circumradius,
+ number_of_copies_in_output=number_of_copies_in_output,
+ verbose=verbose,
+ seed=seed,
+ )
+
+ mesh = meshio.read(outfile)
+ os.remove(outfile)
+ return 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,
+):
+ fh, outfile = tempfile.mkstemp(suffix=".off")
+ os.close(fh)
+
+ _generate_surface_mesh(
+ domain,
+ outfile,
+ bounding_sphere_radius=bounding_sphere_radius,
+ min_facet_angle=min_facet_angle,
+ max_radius_surface_delaunay_ball=max_radius_surface_delaunay_ball,
+ max_facet_distance=max_facet_distance,
+ verbose=verbose,
+ seed=seed,
+ )
+
+ mesh = meshio.read(outfile)
+ os.remove(outfile)
+ return 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,
+):
+ mesh = meshio.read(filename)
+
+ fh, off_file = tempfile.mkstemp(suffix=".off")
+ os.close(fh)
+ meshio.write(off_file, mesh)
+
+ fh, outfile = tempfile.mkstemp(suffix=".mesh")
+ os.close(fh)
+
+ _generate_from_off(
+ off_file,
+ outfile,
+ lloyd=lloyd,
+ odt=odt,
+ perturb=perturb,
+ exude=exude,
+ max_edge_size_at_feature_edges=max_edge_size_at_feature_edges,
+ min_facet_angle=min_facet_angle,
+ max_radius_surface_delaunay_ball=max_radius_surface_delaunay_ball,
+ max_facet_distance=max_facet_distance,
+ max_circumradius_edge_ratio=max_circumradius_edge_ratio,
+ max_cell_circumradius=max_cell_circumradius,
+ verbose=verbose,
+ reorient=reorient,
+ seed=seed,
+ )
+
+ mesh = meshio.read(outfile)
+ os.remove(off_file)
+ os.remove(outfile)
+ return 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,
+):
+ fh, outfile = tempfile.mkstemp(suffix=".mesh")
+ os.close(fh)
+
+ if isinstance(max_cell_circumradius, float):
+ _generate_from_inr(
+ inr_filename,
+ outfile,
+ lloyd=lloyd,
+ odt=odt,
+ perturb=perturb,
+ exude=exude,
+ max_edge_size_at_feature_edges=max_edge_size_at_feature_edges,
+ min_facet_angle=min_facet_angle,
+ max_radius_surface_delaunay_ball=max_radius_surface_delaunay_ball,
+ max_facet_distance=max_facet_distance,
+ max_circumradius_edge_ratio=max_circumradius_edge_ratio,
+ max_cell_circumradius=max_cell_circumradius,
+ verbose=verbose,
+ seed=seed,
+ )
+ else:
+ assert isinstance(max_cell_circumradius, dict)
+ if "default" in max_cell_circumradius.keys():
+ default_max_cell_circumradius = max_cell_circumradius.pop("default")
+ else:
+ default_max_cell_circumradius = 0.0
+
+ max_cell_circumradiuss = list(max_cell_circumradius.values())
+ subdomain_labels = list(max_cell_circumradius.keys())
+
+ _generate_from_inr_with_subdomain_sizing(
+ inr_filename,
+ outfile,
+ default_max_cell_circumradius,
+ max_cell_circumradiuss,
+ subdomain_labels,
+ lloyd=lloyd,
+ odt=odt,
+ perturb=perturb,
+ exude=exude,
+ max_edge_size_at_feature_edges=max_edge_size_at_feature_edges,
+ min_facet_angle=min_facet_angle,
+ max_radius_surface_delaunay_ball=max_radius_surface_delaunay_ball,
+ max_facet_distance=max_facet_distance,
+ max_circumradius_edge_ratio=max_circumradius_edge_ratio,
+ verbose=verbose,
+ seed=seed,
+ )
+
+ mesh = meshio.read(outfile)
+ os.remove(outfile)
+ return mesh
+
+
+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,
+):
+ mesh = meshio.read(filename)
+
+ fh, off_file = tempfile.mkstemp(suffix=".off")
+ os.close(fh)
+ meshio.write(off_file, mesh)
+
+ fh, outfile = tempfile.mkstemp(suffix=".off")
+ os.close(fh)
+
+ _remesh_surface(
+ off_file,
+ outfile,
+ max_edge_size_at_feature_edges=max_edge_size_at_feature_edges,
+ min_facet_angle=min_facet_angle,
+ max_radius_surface_delaunay_ball=max_radius_surface_delaunay_ball,
+ max_facet_distance=max_facet_distance,
+ verbose=verbose,
+ seed=seed,
+ )
+
+ mesh = meshio.read(outfile)
+ os.remove(off_file)
+ os.remove(outfile)
+ return mesh
+
+
+def save_inr(vol, h, fname):
+ """
+ Save a volume (described as a numpy array) to INR format.
+ Code inspired by iso2mesh (http://iso2mesh.sf.net) by Q. Fang
+ INPUTS:
+ - vol: volume as numpy array
+ - h: voxel sizes as list or numpy array
+ - fname: filename for saving the inr file
+ """
+ fid = open(fname, "wb")
+
+ btype, bitlen = {
+ "uint8": ("unsigned fixed", 8),
+ "uint16": ("unsigned fixed", 16),
+ "float32": ("float", 32),
+ "float64": ("float", 64),
+ }[vol.dtype.name]
+
+ header = (
+ "#INRIMAGE-4#{8:s}\nXDIM={0:d}\nYDIM={1:d}\nZDIM={2:d}\nVDIM=1\nTYPE={3:s}\n"
+ + "PIXSIZE={4:d} bits\nCPU=decm\nVX={5:f}\nVY={6:f}\nVZ={7:f}\n"
+ ).format(*vol.shape, btype, bitlen, h[0], h[1], h[2], "{")
+
+ header = header + "\n" * (256 - 4 - len(header)) + "##}\n"
+
+ fid.write(header.encode("ascii"))
+ fid.write(vol.tobytes(order="F"))
+
+
+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,
+):
+ assert vol.dtype in ["uint8", "uint16"]
+ fh, inr_filename = tempfile.mkstemp(suffix=".inr")
+ os.close(fh)
+ save_inr(vol, h, inr_filename)
+ mesh = generate_from_inr(
+ inr_filename,
+ lloyd,
+ odt,
+ perturb,
+ exude,
+ max_edge_size_at_feature_edges,
+ min_facet_angle,
+ max_radius_surface_delaunay_ball,
+ max_facet_distance,
+ max_circumradius_edge_ratio,
+ max_cell_circumradius,
+ verbose,
+ seed,
+ )
+ os.remove(inr_filename)
+ return mesh
diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 0000000..8fe2f47
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,3 @@
+[build-system]
+requires = ["setuptools>=42", "wheel"]
+build-backend = "setuptools.build_meta"
diff --git a/setup.cfg b/setup.cfg
new file mode 100644
index 0000000..be3f277
--- /dev/null
+++ b/setup.cfg
@@ -0,0 +1,46 @@
+[metadata]
+name = pygalmesh
+version = 0.9.1
+url = https://github.com/nschloe/pygalmesh
+author = Nico Schlömer
+author_email = nico.schloemer@gmail.com
+description = Python frontend to CGAL's mesh generation capabilities
+long_description = file: README.md
+long_description_content_type = text/markdown
+license = GPL-3.0-or-later
+license_files = LICENSE
+classifiers =
+ Development Status :: 4 - Beta
+ License :: OSI Approved :: GNU General Public License v3 or later (GPLv3+)
+ Operating System :: OS Independent
+ Programming Language :: Python :: 3
+ Programming Language :: Python :: 3.6
+ Programming Language :: Python :: 3.7
+ Programming Language :: Python :: 3.8
+ Topic :: Scientific/Engineering
+ Topic :: Scientific/Engineering :: Mathematics
+ Topic :: Scientific/Engineering :: Physics
+ Topic :: Scientific/Engineering :: Visualization
+keywords =
+ mathematics
+ physics
+ engineering
+ cgal
+ mesh
+ mesh generation
+
+[options]
+packages = find:
+setup_requires = pybind11 >= 2.5
+install_requires =
+ importlib_metadata;python_version<"3.8"
+ meshio >= 4.0.0, < 5.0.0
+ numpy
+ pybind11 >= 2.5
+python_requires = >=3.6
+
+[options.entry_points]
+console_scripts =
+ pygalmesh-volume-from-surface = pygalmesh._cli:volume_from_surface
+ pygalmesh-remesh-surface = pygalmesh._cli:remesh_surface
+ pygalmesh-from-inr = pygalmesh._cli:inr
diff --git a/setup.py b/setup.py
new file mode 100644
index 0000000..86520a7
--- /dev/null
+++ b/setup.py
@@ -0,0 +1,126 @@
+import os
+import sys
+
+import setuptools
+from setuptools import Extension, setup
+from setuptools.command.build_ext import build_ext
+
+
+#
+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 __str__(self):
+ import pybind11
+
+ return pybind11.get_include()
+
+
+# cf http://bugs.python.org/issue26689
+def has_flag(compiler, flagname):
+ """Return a boolean indicating whether a flag name is supported on
+ the specified compiler.
+ """
+ import os
+ import tempfile
+
+ with tempfile.NamedTemporaryFile("w", suffix=".cpp", delete=False) as f:
+ f.write("int main (int argc, char **argv) { return 0; }")
+ fname = f.name
+ try:
+ compiler.compile([fname], extra_postargs=[flagname])
+ except setuptools.distutils.errors.CompileError:
+ return False
+ finally:
+ try:
+ os.remove(fname)
+ except OSError:
+ pass
+ return True
+
+
+def cpp_flag(compiler):
+ """Return the -std=c++[11/14/17] compiler flag.
+ The newer version is prefered over c++11 (when it is available).
+ """
+ flags = ["-std=c++17", "-std=c++14", "-std=c++11"]
+
+ for flag in flags:
+ if has_flag(compiler, flag):
+ return flag
+
+ raise RuntimeError("Unsupported compiler -- at least C++11 support " "is needed!")
+
+
+class BuildExt(build_ext):
+ """A custom build extension for adding compiler-specific options."""
+
+ c_opts = {
+ "msvc": ["/EHsc"],
+ "unix": [],
+ }
+ l_opts = {
+ "msvc": [],
+ "unix": [],
+ }
+
+ if sys.platform == "darwin":
+ darwin_opts = ["-stdlib=libc++", "-mmacosx-version-min=10.7"]
+ c_opts["unix"] += darwin_opts
+ l_opts["unix"] += darwin_opts
+
+ def build_extensions(self):
+ ct = self.compiler.compiler_type
+ opts = self.c_opts.get(ct, [])
+ link_opts = self.l_opts.get(ct, [])
+ if ct == "unix":
+ opts.append(cpp_flag(self.compiler))
+ if has_flag(self.compiler, "-fvisibility=hidden"):
+ opts.append("-fvisibility=hidden")
+
+ for ext in self.extensions:
+ ext.define_macros = [
+ ("VERSION_INFO", '"{}"'.format(self.distribution.get_version()))
+ ]
+ ext.extra_compile_args = opts
+ ext.extra_link_args = link_opts
+ build_ext.build_extensions(self)
+
+
+ext_modules = [
+ Extension(
+ "_pygalmesh",
+ # Sort input source files to ensure bit-for-bit reproducible builds
+ # (https://github.com/pybind/python_example/pull/53)
+ sorted(
+ [
+ "src/generate.cpp",
+ "src/generate_2d.cpp",
+ "src/generate_from_inr.cpp",
+ "src/generate_from_off.cpp",
+ "src/generate_periodic.cpp",
+ "src/generate_surface_mesh.cpp",
+ "src/remesh_surface.cpp",
+ "src/pybind11.cpp",
+ ]
+ ),
+ include_dirs=[
+ os.environ.get("EIGEN_INCLUDE_DIR", "/usr/include/eigen3/"),
+ # Path to pybind11 headers
+ get_pybind_include(),
+ ],
+ language="c++",
+ # no CGAL libraries necessary from CGAL 5.0 onwards
+ libraries=["gmp", "mpfr"],
+ )
+]
+
+if __name__ == "__main__":
+ setup(
+ cmdclass={"build_ext": BuildExt},
+ ext_modules=ext_modules,
+ zip_safe=False,
+ )
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
new file mode 100644
index 0000000..249ce40
--- /dev/null
+++ b/src/CMakeLists.txt
@@ -0,0 +1,27 @@
+FIND_PACKAGE(pybind11 REQUIRED)
+# include_directories(${PYBIND11_INCLUDE_DIR})
+
+FIND_PACKAGE(Eigen3 REQUIRED)
+include_directories(${EIGEN3_INCLUDE_DIR})
+
+FIND_PACKAGE(CGAL REQUIRED)
+include(${CGAL_USE_FILE})
+
+FILE(GLOB pygalmesh_SRCS "${CMAKE_CURRENT_SOURCE_DIR}/*.cpp")
+# FILE(GLOB pygalmesh_HEADERS "${CMAKE_CURRENT_SOURCE_DIR}/*.hpp")
+
+pybind11_add_module(pygalmesh ${pygalmesh_SRCS})
+
+# 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}
+# )
diff --git a/src/domain.hpp b/src/domain.hpp
new file mode 100644
index 0000000..ff8cd9c
--- /dev/null
+++ b/src/domain.hpp
@@ -0,0 +1,490 @@
+#ifndef DOMAIN_HPP
+#define DOMAIN_HPP
+
+#include
+#include
+#include
+#include
+#include
+
+namespace pygalmesh {
+
+class DomainBase
+{
+ public:
+
+ virtual ~DomainBase() = default;
+
+ virtual
+ double
+ eval(const std::array & x) const = 0;
+
+ virtual
+ double
+ get_bounding_sphere_squared_radius() const = 0;
+
+ virtual
+ std::vector>>
+ get_features() const
+ {
+ return {};
+ };
+};
+
+class Translate: public pygalmesh::DomainBase
+{
+ public:
+ Translate(
+ const std::shared_ptr & domain,
+ const std::array & direction
+ ):
+ domain_(domain),
+ direction_(Eigen::Vector3d(direction.data())),
+ translated_features_(translate_features(domain->get_features(), direction_))
+ {
+ }
+
+ virtual ~Translate() = default;
+
+ std::vector>>
+ translate_features(
+ const std::vector>> & features,
+ const Eigen::Vector3d & direction
+ ) const
+ {
+ std::vector>> translated_features;
+ for (const auto & feature: features) {
+ std::vector> translated_feature;
+ for (const auto & point: feature) {
+ const std::array 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 & x) const
+ {
+ const std::array 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>>
+ get_features() const
+ {
+ return translated_features_;
+ };
+
+ private:
+ const std::shared_ptr domain_;
+ const Eigen::Vector3d direction_;
+ const std::vector>> translated_features_;
+};
+
+class Rotate: public pygalmesh::DomainBase
+{
+ public:
+ Rotate(
+ const std::shared_ptr & domain,
+ const std::array & 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 & 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>>
+ rotate_features(
+ const std::vector>> & features
+ ) const
+ {
+ std::vector>> rotated_features;
+ for (const auto & feature: features) {
+ std::vector> 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>>
+ get_features() const
+ {
+ return rotated_features_;
+ };
+
+ private:
+ const std::shared_ptr domain_;
+ const Eigen::Vector3d normalized_axis_;
+ const double sinAngle_;
+ const double cosAngle_;
+ const std::vector>> rotated_features_;
+};
+
+class Scale: public pygalmesh::DomainBase
+{
+ public:
+ Scale(
+ std::shared_ptr & 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 & 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>>
+ scale_features(
+ const std::vector>> & features
+ ) const
+ {
+ std::vector>> scaled_features;
+ for (const auto & feature: features) {
+ std::vector> 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>>
+ get_features() const
+ {
+ return scaled_features_;
+ };
+
+ private:
+ std::shared_ptr domain_;
+ const double alpha_;
+ const std::vector>> scaled_features_;
+};
+
+class Stretch: public pygalmesh::DomainBase
+{
+ public:
+ Stretch(
+ std::shared_ptr & domain,
+ const std::array & 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 & 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>>
+ stretch_features(
+ const std::vector>> & features
+ ) const
+ {
+ std::vector>> stretched_features;
+ for (const auto & feature: features) {
+ std::vector> 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>>
+ get_features() const
+ {
+ return stretched_features_;
+ };
+
+ private:
+ std::shared_ptr domain_;
+ const Eigen::Vector3d normalized_direction_;
+ const double alpha_;
+ const std::vector>> stretched_features_;
+};
+
+class Intersection: public pygalmesh::DomainBase
+{
+ public:
+ explicit Intersection(
+ std::vector> & domains
+ ):
+ domains_(domains)
+ {
+ }
+
+ virtual ~Intersection() = default;
+
+ virtual
+ double
+ eval(const std::array & x) const
+ {
+ // TODO find a differentiable expression
+ double maxval = std::numeric_limits::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::max();
+ for (const auto & domain: domains_) {
+ min = std::min(min, domain->get_bounding_sphere_squared_radius());
+ }
+ return min;
+ }
+
+ virtual
+ std::vector>>
+ get_features() const
+ {
+ std::vector>> 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> domains_;
+};
+
+class Union: public pygalmesh::DomainBase
+{
+ public:
+ explicit Union(
+ std::vector> & domains
+ ):
+ domains_(domains)
+ {
+ }
+
+ virtual ~Union() = default;
+
+ virtual
+ double
+ eval(const std::array & x) const
+ {
+ // TODO find a differentiable expression
+ double minval = std::numeric_limits::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>>
+ get_features() const
+ {
+ std::vector>> 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> domains_;
+};
+
+class Difference: public pygalmesh::DomainBase
+{
+ public:
+ Difference(
+ std::shared_ptr & domain0,
+ std::shared_ptr & domain1
+ ):
+ domain0_(domain0),
+ domain1_(domain1)
+ {
+ }
+
+ virtual ~Difference() = default;
+
+ virtual
+ double
+ eval(const std::array & 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>>
+ get_features() const
+ {
+ std::vector>> 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 domain0_;
+ std::shared_ptr domain1_;
+};
+
+} // namespace pygalmesh
+#endif // DOMAIN_HPP
diff --git a/src/generate.cpp b/src/generate.cpp
new file mode 100644
index 0000000..d4fef73
--- /dev/null
+++ b/src/generate.cpp
@@ -0,0 +1,183 @@
+#define CGAL_MESH_3_VERBOSE 1
+
+#include "generate.hpp"
+
+#include
+
+#include
+#include
+#include
+
+#include
+#include
+#include
+
+namespace pygalmesh {
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+
+typedef CGAL::Mesh_domain_with_polyline_features_3> Mesh_domain;
+
+// Triangulation
+typedef CGAL::Mesh_triangulation_3::type Tr;
+typedef CGAL::Mesh_complex_3_in_triangulation_3 C3t3;
+
+// Mesh Criteria
+typedef CGAL::Mesh_criteria_3
Mesh_criteria;
+typedef Mesh_criteria::Edge_criteria Edge_criteria;
+typedef Mesh_criteria::Facet_criteria Facet_criteria;
+typedef Mesh_criteria::Cell_criteria Cell_criteria;
+
+// translate vector> to list>
+std::list>
+translate_feature_edges(
+ const std::vector>> & feature_edges
+ )
+{
+ std::list> polylines;
+ for (const auto & feature_edge: feature_edges) {
+ std::vector 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 & domain,
+ const std::string & outfile,
+ const std::vector>> & extra_feature_edges,
+ const double bounding_sphere_radius,
+ const bool lloyd,
+ const bool odt,
+ const bool perturb,
+ const bool exude,
+ //
+ const double max_edge_size_at_feature_edges_value,
+ const std::shared_ptr & max_edge_size_at_feature_edges_field,
+ //
+ const double min_facet_angle,
+ //
+ const double max_radius_surface_delaunay_ball_value,
+ const std::shared_ptr & max_radius_surface_delaunay_ball_field,
+ //
+ const double max_facet_distance_value,
+ const std::shared_ptr & max_facet_distance_field,
+ //
+ const double max_circumradius_edge_ratio,
+ //
+ const double max_cell_circumradius_value,
+ const std::shared_ptr & max_cell_circumradius_field,
+ //
+ const bool verbose,
+ const int seed
+ )
+{
+ CGAL::get_default_random() = CGAL::Random(seed);
+
+ const double bounding_sphere_radius2 = bounding_sphere_radius > 0 ?
+ bounding_sphere_radius*bounding_sphere_radius :
+ // some wiggle room
+ 1.01 * domain->get_bounding_sphere_squared_radius();
+
+ // wrap domain
+ const auto d = [&](K::Point_3 p) {
+ return domain->eval({p.x(), p.y(), p.z()});
+ };
+
+ Mesh_domain cgal_domain = Mesh_domain::create_implicit_mesh_domain(
+ d,
+ K::Sphere_3(CGAL::ORIGIN, bounding_sphere_radius2)
+ );
+
+ // cgal_domain.detect_features();
+
+ const auto native_features = translate_feature_edges(domain->get_features());
+ cgal_domain.add_features(native_features.begin(), native_features.end());
+
+ const auto polylines = translate_feature_edges(extra_feature_edges);
+ cgal_domain.add_features(polylines.begin(), polylines.end());
+
+ // perhaps there's a more elegant solution here
+ // see
+ if (!verbose) {
+ // suppress output
+ std::cerr.setstate(std::ios_base::failbit);
+ }
+
+ // Build the float/field values according to
+ // .
+
+ // nested ternary operator
+ const auto facet_criteria = max_radius_surface_delaunay_ball_field ? (
+ max_facet_distance_field ?
+ Facet_criteria(
+ min_facet_angle,
+ [&](K::Point_3 p, const int, const Mesh_domain::Index&) {
+ return max_radius_surface_delaunay_ball_field->eval({p.x(), p.y(), p.z()});
+ },
+ [&](K::Point_3 p, const int, const Mesh_domain::Index&) {
+ return max_facet_distance_field->eval({p.x(), p.y(), p.z()});
+ }
+ ) : Facet_criteria(
+ min_facet_angle,
+ [&](K::Point_3 p, const int, const Mesh_domain::Index&) {
+ return max_radius_surface_delaunay_ball_field->eval({p.x(), p.y(), p.z()});
+ },
+ max_facet_distance_value
+ )
+ ) : (
+ max_facet_distance_field ?
+ Facet_criteria(
+ min_facet_angle,
+ max_radius_surface_delaunay_ball_value,
+ [&](K::Point_3 p, const int, const Mesh_domain::Index&) {
+ return max_facet_distance_field->eval({p.x(), p.y(), p.z()});
+ }
+ ) : Facet_criteria(
+ min_facet_angle,
+ max_radius_surface_delaunay_ball_value,
+ max_facet_distance_value
+ )
+ );
+
+ const auto edge_criteria = max_edge_size_at_feature_edges_field ?
+ Edge_criteria(
+ [&](K::Point_3 p, const int, const Mesh_domain::Index&) {
+ return max_edge_size_at_feature_edges_field->eval({p.x(), p.y(), p.z()});
+ }) : Edge_criteria(max_edge_size_at_feature_edges_value);
+
+ const auto cell_criteria = max_cell_circumradius_field ?
+ Cell_criteria(
+ max_circumradius_edge_ratio,
+ [&](K::Point_3 p, const int, const Mesh_domain::Index&) {
+ return max_cell_circumradius_field->eval({p.x(), p.y(), p.z()});
+ }) : Cell_criteria(max_circumradius_edge_ratio, max_cell_circumradius_value);
+
+ const auto criteria = Mesh_criteria(edge_criteria, facet_criteria, cell_criteria);
+
+ // Mesh generation
+ C3t3 c3t3 = CGAL::make_mesh_3(
+ 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
diff --git a/src/generate.hpp b/src/generate.hpp
new file mode 100644
index 0000000..f1d42e4
--- /dev/null
+++ b/src/generate.hpp
@@ -0,0 +1,46 @@
+#ifndef GENERATE_HPP
+#define GENERATE_HPP
+
+#include "domain.hpp"
+#include "sizing_field.hpp"
+
+#include
+#include
+#include
+#include
+
+namespace pygalmesh {
+
+void generate_mesh(
+ const std::shared_ptr & domain,
+ const std::string & outfile,
+ const std::vector>> & extra_feature_edges = {},
+ const double bounding_sphere_radius = 0.0,
+ const bool lloyd = false,
+ const bool odt = false,
+ const bool perturb = true,
+ const bool exude = true,
+ //
+ const double max_edge_size_at_feature_edges_value = 0.0, // std::numeric_limits::max(),
+ const std::shared_ptr & max_edge_size_at_feature_edges_field = nullptr,
+ //
+ const double min_facet_angle = 0.0,
+ //
+ const double max_radius_surface_delaunay_ball_value = 0.0,
+ const std::shared_ptr & max_radius_surface_delaunay_ball_field = nullptr,
+ //
+ const double max_facet_distance_value = 0.0,
+ const std::shared_ptr & max_facet_distance_field = nullptr,
+ //
+ const double max_circumradius_edge_ratio = 0.0,
+ //
+ const double max_cell_circumradius_value = 0.0,
+ const std::shared_ptr & max_cell_circumradius_field = nullptr,
+ //
+ const bool verbose = true,
+ const int seed = 0
+ );
+
+} // namespace pygalmesh
+
+#endif // GENERATE_HPP
diff --git a/src/generate_2d.cpp b/src/generate_2d.cpp
new file mode 100644
index 0000000..5b573ce
--- /dev/null
+++ b/src/generate_2d.cpp
@@ -0,0 +1,96 @@
+#define CGAL_MESH_3_VERBOSE 1
+
+#include "generate_2d.hpp"
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+namespace pygalmesh {
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+typedef CGAL::Delaunay_mesh_vertex_base_2 Vb;
+typedef CGAL::Delaunay_mesh_face_base_2 Fb;
+typedef CGAL::Triangulation_data_structure_2 Tds;
+typedef CGAL::Constrained_Delaunay_triangulation_2 CDT;
+typedef CGAL::Delaunay_mesh_size_criteria_2 Criteria;
+typedef CDT::Vertex_handle Vertex_handle;
+typedef CDT::Point Point;
+
+std::tuple>, std::vector>>
+generate_2d(
+ const std::vector> & points,
+ const std::vector> & constraints,
+ // See
+ // https://doc.cgal.org/latest/Mesh_2/classCGAL_1_1Delaunay__mesh__size__criteria__2.html#a58b0186eae407ba76b8f4a3d0aa85a1a
+ // for what the bounds mean. Spoiler:
+ // B = circumradius / shortest_edge,
+ // relates to the smallest angle via sin(alpha_min) = 1 / (2B)
+ // edge_size S: "all segments of all triangles must be shorter than a bound S."
+ const double max_circumradius_shortest_edge_ratio,
+ const double max_edge_size,
+ const int num_lloyd_steps
+)
+{
+ CDT cdt;
+ // construct a constrained triangulation
+ std::vector vertices(points.size());
+ int k = 0;
+ for (auto pt: points) {
+ vertices[k] = cdt.insert(Point(pt[0], pt[1]));
+ k++;
+ }
+ for (auto c: constraints) {
+ cdt.insert_constraint(vertices[c[0]], vertices[c[1]]);
+ }
+
+ // create proper mesh
+ CGAL::refine_Delaunay_mesh_2(
+ cdt,
+ Criteria(
+ 0.25 / (max_circumradius_shortest_edge_ratio * max_circumradius_shortest_edge_ratio),
+ max_edge_size
+ )
+ );
+
+ if (num_lloyd_steps > 0) {
+ CGAL::lloyd_optimize_mesh_2(
+ cdt,
+ CGAL::parameters::max_iteration_number = num_lloyd_steps
+ );
+ }
+
+ // convert points to vector of arrays
+ std::map vertex_index;
+ std::vector> out_points(cdt.number_of_vertices());
+ k = 0;
+ for (auto vit: cdt.finite_vertex_handles()) {
+ out_points[k][0] = vit->point()[0];
+ out_points[k][1] = vit->point()[1];
+ vertex_index[vit] = k;
+ k++;
+ }
+
+ // https://github.com/CGAL/cgal/issues/5068#issuecomment-706213606
+ auto nb_faces = 0u;
+ for (auto fit: cdt.finite_face_handles()) {
+ if(fit->is_in_domain()) ++nb_faces;
+ }
+ std::vector> out_cells(nb_faces);
+ k = 0;
+ for (auto fit: cdt.finite_face_handles()) {
+ if(!fit->is_in_domain()) continue;
+ out_cells[k][0] = vertex_index[fit->vertex(0)];
+ out_cells[k][1] = vertex_index[fit->vertex(1)];
+ out_cells[k][2] = vertex_index[fit->vertex(2)];
+ k++;
+ }
+
+ return std::make_tuple(out_points, out_cells);
+}
+
+} // namespace pygalmesh
diff --git a/src/generate_2d.hpp b/src/generate_2d.hpp
new file mode 100644
index 0000000..9b55073
--- /dev/null
+++ b/src/generate_2d.hpp
@@ -0,0 +1,25 @@
+#ifndef GENERATE_2D_HPP
+#define GENERATE_2D_HPP
+
+#include
+#include
+#include
+
+namespace pygalmesh {
+
+std::tuple>, std::vector>>
+generate_2d(
+ const std::vector> & points,
+ const std::vector> & constraints,
+ const double max_circumradius_shortest_edge_ratio,
+ // https://github.com/CGAL/cgal/issues/5061#issuecomment-705520984
+ // the "default" size criterion for a triangle in the 2D mesh generator refers to its
+ // edge lengths. In the output mesh, all segments of all triangles must be shorter
+ // than the given bound.
+ const double max_edge_size,
+ const int num_lloyd_steps
+);
+
+} // namespace pygalmesh
+
+#endif // GENERATE_2D_HPP
diff --git a/src/generate_from_inr.cpp b/src/generate_from_inr.cpp
new file mode 100644
index 0000000..830f33f
--- /dev/null
+++ b/src/generate_from_inr.cpp
@@ -0,0 +1,164 @@
+#define CGAL_MESH_3_VERBOSE 1
+
+#include "generate_from_inr.hpp"
+
+#include
+
+#include
+#include
+
+#include
+#include
+#include
+
+#include
+#include
+#include
+
+namespace pygalmesh {
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+
+typedef CGAL::Labeled_mesh_domain_3 Mesh_domain;
+
+// Triangulation
+typedef CGAL::Mesh_triangulation_3::type Tr;
+typedef CGAL::Mesh_complex_3_in_triangulation_3 C3t3;
+
+// Mesh Criteria
+typedef CGAL::Mesh_criteria_3
Mesh_criteria;
+typedef Mesh_criteria::Facet_criteria Facet_criteria;
+typedef Mesh_criteria::Cell_criteria Cell_criteria;
+
+typedef CGAL::Mesh_constant_domain_field_3 Sizing_field_cell;
+
+void
+generate_from_inr(
+ const std::string & inr_filename,
+ const std::string & outfile,
+ const bool lloyd,
+ const bool odt,
+ const bool perturb,
+ const bool exude,
+ const double max_edge_size_at_feature_edges,
+ const double min_facet_angle,
+ const double max_radius_surface_delaunay_ball,
+ const double max_facet_distance,
+ const double max_circumradius_edge_ratio,
+ const double max_cell_circumradius,
+ const bool verbose,
+ const int seed
+ )
+{
+ CGAL::get_default_random() = CGAL::Random(seed);
+
+ CGAL::Image_3 image;
+ const bool success = image.read(inr_filename.c_str());
+ if (!success) {
+ throw "Could not read image file";
+ }
+ Mesh_domain cgal_domain = Mesh_domain::create_labeled_image_mesh_domain(image);
+
+ Mesh_criteria criteria(
+ CGAL::parameters::edge_size=max_edge_size_at_feature_edges,
+ CGAL::parameters::facet_angle=min_facet_angle,
+ CGAL::parameters::facet_size=max_radius_surface_delaunay_ball,
+ CGAL::parameters::facet_distance=max_facet_distance,
+ CGAL::parameters::cell_radius_edge_ratio=max_circumradius_edge_ratio,
+ CGAL::parameters::cell_size=max_cell_circumradius
+ );
+
+ // Mesh generation
+ if (!verbose) {
+ // suppress output
+ std::cerr.setstate(std::ios_base::failbit);
+ }
+ C3t3 c3t3 = CGAL::make_mesh_3(
+ 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;
+}
+
+
+void
+generate_from_inr_with_subdomain_sizing(
+ const std::string & inr_filename,
+ const std::string & outfile,
+ const double default_max_cell_circumradius,
+ const std::vector & max_cell_circumradiuss,
+ const std::vector & cell_labels,
+ const bool lloyd,
+ const bool odt,
+ const bool perturb,
+ const bool exude,
+ const double max_edge_size_at_feature_edges,
+ const double min_facet_angle,
+ const double max_radius_surface_delaunay_ball,
+ const double max_facet_distance,
+ const double max_circumradius_edge_ratio,
+ const bool verbose,
+ const int seed
+ )
+{
+ CGAL::get_default_random() = CGAL::Random(seed);
+
+ CGAL::Image_3 image;
+ const bool success = image.read(inr_filename.c_str());
+ if (!success) {
+ throw "Could not read image file";
+ }
+ Mesh_domain cgal_domain = Mesh_domain::create_labeled_image_mesh_domain(image);
+
+ Sizing_field_cell max_cell_circumradius(default_max_cell_circumradius);
+ const int ndimensions = 3;
+ for(std::vector::size_type i(0); i < max_cell_circumradiuss.size(); ++i)
+ max_cell_circumradius.set_size(max_cell_circumradiuss[i], ndimensions, cgal_domain.index_from_subdomain_index(cell_labels[i]));
+
+ Mesh_criteria criteria(
+ CGAL::parameters::edge_size=max_edge_size_at_feature_edges,
+ CGAL::parameters::facet_angle=min_facet_angle,
+ CGAL::parameters::facet_size=max_radius_surface_delaunay_ball,
+ CGAL::parameters::facet_distance=max_facet_distance,
+ CGAL::parameters::cell_radius_edge_ratio=max_circumradius_edge_ratio,
+ CGAL::parameters::cell_size=max_cell_circumradius
+ );
+
+ // Mesh generation
+ if (!verbose) {
+ // suppress output
+ std::cerr.setstate(std::ios_base::failbit);
+ }
+ C3t3 c3t3 = CGAL::make_mesh_3(
+ 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
diff --git a/src/generate_from_inr.hpp b/src/generate_from_inr.hpp
new file mode 100644
index 0000000..ddd60c7
--- /dev/null
+++ b/src/generate_from_inr.hpp
@@ -0,0 +1,48 @@
+#ifndef GENERATE_FROM_INR_HPP
+#define GENERATE_FROM_INR_HPP
+
+#include
+#include
+
+namespace pygalmesh {
+
+void generate_from_inr(
+ const std::string & inr_filename,
+ const std::string & outfile,
+ const bool lloyd = false,
+ const bool odt = false,
+ const bool perturb = true,
+ const bool exude = true,
+ const double max_edge_size_at_feature_edges = 0.0, // std::numeric_limits::max(),
+ const double min_facet_angle = 0.0,
+ const double max_radius_surface_delaunay_ball = 0.0,
+ const double max_facet_distance = 0.0,
+ const double max_circumradius_edge_ratio = 0.0,
+ const double max_cell_circumradius = 0.0,
+ const bool verbose = true,
+ const int seed = 0
+ );
+
+void
+generate_from_inr_with_subdomain_sizing(
+ const std::string & inr_filename,
+ const std::string & outfile,
+ const double default_max_cell_circumradius,
+ const std::vector & max_cell_circumradiuss,
+ const std::vector & cell_labels,
+ const bool lloyd = false,
+ const bool odt = false,
+ const bool perturb = true,
+ const bool exude = true,
+ const double max_edge_size_at_feature_edges = 0.0,
+ const double min_facet_angle = 0.0,
+ const double max_radius_surface_delaunay_ball = 0.0,
+ const double max_facet_distance = 0.0,
+ const double max_circumradius_edge_ratio = 0.0,
+ const bool verbose = true,
+ const int seed = 0
+ );
+
+} // namespace pygalmesh
+
+#endif // GENERATE_FROM_INR_HPP
diff --git a/src/generate_from_off.cpp b/src/generate_from_off.cpp
new file mode 100644
index 0000000..2696075
--- /dev/null
+++ b/src/generate_from_off.cpp
@@ -0,0 +1,138 @@
+#include "generate_from_off.hpp"
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+// IO
+#include
+
+// for re-orientation
+#include
+#include
+#include
+#include
+
+// for sharp features
+//#include
+
+namespace pygalmesh {
+
+// Domain
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+typedef CGAL::Polyhedron_3 Polyhedron;
+typedef CGAL::Polyhedral_mesh_domain_3 Mesh_domain;
+// for sharp features
+//typedef CGAL::Polyhedral_mesh_domain_with_features_3 Mesh_domain;
+//typedef CGAL::Mesh_polyhedron_3::type Polyhedron;
+
+// Triangulation
+typedef CGAL::Mesh_triangulation_3::type Tr;
+
+typedef CGAL::Mesh_complex_3_in_triangulation_3 C3t3;
+
+// Criteria
+typedef CGAL::Mesh_criteria_3
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 max_edge_size_at_feature_edges,
+ const double min_facet_angle,
+ const double max_radius_surface_delaunay_ball,
+ const double max_facet_distance,
+ const double max_circumradius_edge_ratio,
+ const double max_cell_circumradius,
+ const bool verbose,
+ const bool reorient,
+ const int seed
+) {
+ CGAL::get_default_random() = CGAL::Random(seed);
+
+ std::ifstream input(infile);
+ Polyhedron polyhedron;
+ // fix the orientation of the faces of the input file
+ if (reorient) {
+ std::stringstream msg;
+ msg << "fixing face orientation for \"" << infile <<"\""<< std::endl;
+ std::vector points;
+ std::vector > polygons;
+
+ if(!input || !CGAL::read_OFF(input, points, polygons) || points.empty())
+ {
+ std::stringstream msg;
+ msg << "Cannot read .off file \"" << infile <<"\""<< std::endl;
+ throw std::runtime_error(msg.str());
+ }
+ // orient the polygons
+ CGAL::Polygon_mesh_processing::orient_polygon_soup(points, polygons);
+ // create the polyhedron
+ CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh(points, polygons, polyhedron);
+
+ } else {
+
+ // Create input polyhedron
+ input >> polyhedron;
+ if (!input.good()) {
+ // Even if the mesh exists, it may not be valid, see
+ //
+ std::stringstream msg;
+ msg << "Invalid input file \"" << infile << "\"" << std::endl;
+ msg << "If this is due to wrong face orientation, retry with the option --reorient \"" << std::endl;
+ throw std::runtime_error(msg.str());
+ }
+ }
+
+ input.close();
+
+ // Create domain
+ Mesh_domain cgal_domain(polyhedron);
+
+ // Get sharp features
+ // cgal_domain.detect_features();
+
+
+ // Mesh criteria
+ Mesh_criteria criteria(
+ CGAL::parameters::edge_size = max_edge_size_at_feature_edges,
+ CGAL::parameters::facet_angle = min_facet_angle,
+ CGAL::parameters::facet_size = max_radius_surface_delaunay_ball,
+ CGAL::parameters::facet_distance = max_facet_distance,
+ CGAL::parameters::cell_radius_edge_ratio = max_circumradius_edge_ratio,
+ CGAL::parameters::cell_size = max_cell_circumradius);
+
+ // Mesh generation
+ if (!verbose) {
+ // suppress output
+ std::cerr.setstate(std::ios_base::failbit);
+ }
+ C3t3 c3t3 = CGAL::make_mesh_3(
+ 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
diff --git a/src/generate_from_off.hpp b/src/generate_from_off.hpp
new file mode 100644
index 0000000..f3e5a21
--- /dev/null
+++ b/src/generate_from_off.hpp
@@ -0,0 +1,30 @@
+#ifndef GENERATE_FROM_OFF_HPP
+#define GENERATE_FROM_OFF_HPP
+
+#include
+#include
+
+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 max_edge_size_at_feature_edges = 0.0, // std::numeric_limits::max(),
+ const double min_facet_angle = 0.0,
+ const double max_radius_surface_delaunay_ball = 0.0,
+ const double max_facet_distance = 0.0,
+ const double max_circumradius_edge_ratio = 0.0,
+ const double max_cell_circumradius = 0.0,
+ const bool verbose = true,
+ const bool reorient = false,
+ const int seed = 0
+ );
+
+} // namespace pygalmesh
+
+#endif // GENERATE_FROM_OFF_HPP
diff --git a/src/generate_periodic.cpp b/src/generate_periodic.cpp
new file mode 100644
index 0000000..2a4e8ce
--- /dev/null
+++ b/src/generate_periodic.cpp
@@ -0,0 +1,123 @@
+#define CGAL_MESH_3_VERBOSE 1
+
+#include "generate_periodic.hpp"
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include // CGAL_PI
+#include
+#include
+#include
+
+
+namespace pygalmesh {
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+
+typedef CGAL::Labeled_mesh_domain_3 Periodic_mesh_domain;
+
+// Triangulation
+typedef CGAL::Periodic_3_mesh_triangulation_3::type Tr;
+typedef CGAL::Mesh_complex_3_in_triangulation_3 C3t3;
+
+// Mesh Criteria
+typedef CGAL::Mesh_criteria_3
Mesh_criteria;
+typedef Mesh_criteria::Facet_criteria Facet_criteria;
+typedef Mesh_criteria::Cell_criteria Cell_criteria;
+
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+typedef K::FT FT;
+typedef K::Point_3 Point;
+typedef K::Iso_cuboid_3 Iso_cuboid;
+// Domain
+typedef FT (Function)(const Point&);
+typedef CGAL::Labeled_mesh_domain_3 Periodic_mesh_domain;
+// Triangulation
+typedef CGAL::Periodic_3_mesh_triangulation_3::type Tr;
+typedef CGAL::Mesh_complex_3_in_triangulation_3 C3t3;
+// Criteria
+typedef CGAL::Mesh_criteria_3
Periodic_mesh_criteria;
+// To avoid verbose function and named parameters call
+using namespace CGAL::parameters;
+
+void
+generate_periodic_mesh(
+ const std::shared_ptr & domain,
+ const std::string & outfile,
+ const std::array bounding_cuboid,
+ const bool lloyd,
+ const bool odt,
+ const bool perturb,
+ const bool exude,
+ const double max_edge_size_at_feature_edges,
+ const double min_facet_angle,
+ const double max_radius_surface_delaunay_ball,
+ const double max_facet_distance,
+ const double max_circumradius_edge_ratio,
+ const double max_cell_circumradius,
+ const int number_of_copies_in_output,
+ const bool verbose,
+ const int seed
+ )
+{
+ CGAL::get_default_random() = CGAL::Random(seed);
+
+ K::Iso_cuboid_3 cuboid(
+ bounding_cuboid[0],
+ bounding_cuboid[1],
+ bounding_cuboid[2],
+ bounding_cuboid[3],
+ bounding_cuboid[4],
+ bounding_cuboid[5]
+ );
+
+ // wrap domain
+ const auto d = [&](K::Point_3 p) {
+ return domain->eval({p.x(), p.y(), p.z()});
+ };
+ Periodic_mesh_domain cgal_domain =
+ Periodic_mesh_domain::create_implicit_mesh_domain(d, cuboid);
+
+ Mesh_criteria criteria(
+ CGAL::parameters::edge_size=max_edge_size_at_feature_edges,
+ CGAL::parameters::facet_angle=min_facet_angle,
+ CGAL::parameters::facet_size=max_radius_surface_delaunay_ball,
+ CGAL::parameters::facet_distance=max_facet_distance,
+ CGAL::parameters::cell_radius_edge_ratio=max_circumradius_edge_ratio,
+ CGAL::parameters::cell_size=max_cell_circumradius
+ );
+
+ // Mesh generation
+ if (!verbose) {
+ // suppress output
+ std::cerr.setstate(std::ios_base::failbit);
+ }
+ C3t3 c3t3 = CGAL::make_periodic_3_mesh_3(
+ 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);
+ CGAL::output_periodic_mesh_to_medit(medit_file, c3t3, number_of_copies_in_output);
+ medit_file.close();
+
+ return;
+}
+
+} // namespace pygalmesh
diff --git a/src/generate_periodic.hpp b/src/generate_periodic.hpp
new file mode 100644
index 0000000..25ad5ca
--- /dev/null
+++ b/src/generate_periodic.hpp
@@ -0,0 +1,33 @@
+#ifndef GENERATE_PERIODIC_HPP
+#define GENERATE_PERIODIC_HPP
+
+#include "domain.hpp"
+
+#include
+#include
+#include
+
+namespace pygalmesh {
+
+void generate_periodic_mesh(
+ const std::shared_ptr & domain,
+ const std::string & outfile,
+ const std::array bounding_cuboid,
+ const bool lloyd = false,
+ const bool odt = false,
+ const bool perturb = true,
+ const bool exude = true,
+ const double max_edge_size_at_feature_edges = 0.0, // std::numeric_limits::max(),
+ const double min_facet_angle = 0.0,
+ const double max_radius_surface_delaunay_ball = 0.0,
+ const double max_facet_distance = 0.0,
+ const double max_circumradius_edge_ratio = 0.0,
+ const double max_cell_circumradius = 0.0,
+ const int number_of_copies_in_output = 1,
+ const bool verbose = true,
+ const int seed = 0
+ );
+
+} // namespace pygalmesh
+
+#endif // GENERATE_PERIODIC_HPP
diff --git a/src/generate_surface_mesh.cpp b/src/generate_surface_mesh.cpp
new file mode 100644
index 0000000..4d3f367
--- /dev/null
+++ b/src/generate_surface_mesh.cpp
@@ -0,0 +1,102 @@
+#define CGAL_SURFACE_MESHER_VERBOSE 1
+
+#include "generate_surface_mesh.hpp"
+
+#include
+#include
+#include
+#include
+#include
+#include
+
+namespace pygalmesh {
+
+// default triangulation for Surface_mesher
+typedef CGAL::Surface_mesh_default_triangulation_3 Tr;
+// c2t3
+typedef CGAL::Complex_2_in_triangulation_3 C2t3;
+typedef Tr::Geom_traits GT;
+
+// Wrapper for DomainBase for translating to GT.
+class CgalDomainWrapper
+{
+ public:
+ explicit CgalDomainWrapper(const std::shared_ptr & 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 domain_;
+};
+
+typedef CGAL::Implicit_surface_3 Surface_3;
+
+void
+generate_surface_mesh(
+ const std::shared_ptr & domain,
+ const std::string & outfile,
+ const double bounding_sphere_radius,
+ const double min_facet_angle,
+ const double max_radius_surface_delaunay_ball,
+ const double max_facet_distance,
+ const bool verbose,
+ const int seed
+ )
+{
+ CGAL::get_default_random() = CGAL::Random(seed);
+
+ const double bounding_sphere_radius2 = bounding_sphere_radius > 0 ?
+ bounding_sphere_radius*bounding_sphere_radius :
+ // 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 criteria(
+ min_facet_angle,
+ max_radius_surface_delaunay_ball,
+ max_facet_distance
+ );
+
+ 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
diff --git a/src/generate_surface_mesh.hpp b/src/generate_surface_mesh.hpp
new file mode 100644
index 0000000..df88da4
--- /dev/null
+++ b/src/generate_surface_mesh.hpp
@@ -0,0 +1,24 @@
+#ifndef GENERATE_SURFACE_MESH_HPP
+#define GENERATE_SURFACE_MESH_HPP
+
+#include "domain.hpp"
+
+#include
+#include
+
+namespace pygalmesh {
+
+void generate_surface_mesh(
+ const std::shared_ptr & domain,
+ const std::string & outfile,
+ const double bounding_sphere_radius = 0.0,
+ const double min_facet_angle = 0.0,
+ const double max_radius_surface_delaunay_ball = 0.0,
+ const double max_facet_distance = 0.0,
+ const bool verbose = true,
+ const int seed = 0
+ );
+
+} // namespace pygalmesh
+
+#endif // GENERATE_SURFACE_MESH_HPP
diff --git a/src/polygon2d.hpp b/src/polygon2d.hpp
new file mode 100644
index 0000000..40144b1
--- /dev/null
+++ b/src/polygon2d.hpp
@@ -0,0 +1,308 @@
+#ifndef POLYGON2D_HPP
+#define POLYGON2D_HPP
+
+#include "domain.hpp"
+
+#include
+#include
+#include
+#include
+#include
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+
+namespace pygalmesh {
+
+class Polygon2D {
+ public:
+ explicit Polygon2D(const std::vector> & _points):
+ points(vector_to_cgal_points(_points))
+ {
+ }
+
+ virtual ~Polygon2D() = default;
+
+ std::vector
+ vector_to_cgal_points(const std::vector> & _points) const
+ {
+ std::vector 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 & 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 points;
+};
+
+
+class Extrude: public pygalmesh::DomainBase {
+ public:
+ Extrude(
+ const std::shared_ptr & poly,
+ const std::array & direction,
+ const double alpha = 0.0,
+ const double max_edge_size_at_feature_edges = 0.0
+ ):
+ poly_(poly),
+ direction_(direction),
+ alpha_(alpha),
+ max_edge_size_at_feature_edges_(max_edge_size_at_feature_edges)
+ {
+ }
+
+ virtual ~Extrude() = default;
+
+ virtual
+ double
+ eval(const std::array & x) const
+ {
+ if (x[2] < 0.0 || x[2] > direction_[2]) {
+ return 1.0;
+ }
+
+ const double beta = x[2] / direction_[2];
+
+ std::array x2 = {
+ x[0] - beta * direction_[0],
+ x[1] - beta * direction_[1]
+ };
+
+ if (alpha_ != 0.0) {
+ std::array 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>>
+ get_features() const
+ {
+ std::vector>> 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> 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 max_edge_size_at_feature_edges. 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(max_edge_size_at_feature_edges_ > 0.0);
+ const size_t n = int(l / max_edge_size_at_feature_edges_ - 0.5) + 1;
+ std::vector> 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 poly_;
+ const std::array direction_;
+ const double alpha_;
+ const double max_edge_size_at_feature_edges_;
+};
+
+
+class ring_extrude: public pygalmesh::DomainBase {
+ public:
+ ring_extrude(
+ const std::shared_ptr & poly,
+ const double max_edge_size_at_feature_edges
+ ):
+ poly_(poly),
+ max_edge_size_at_feature_edges_(max_edge_size_at_feature_edges)
+ {
+ assert(max_edge_size_at_feature_edges > 0.0);
+ }
+
+ virtual ~ring_extrude() = default;
+
+ virtual
+ double
+ eval(const std::array & 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>>
+ get_features() const
+ {
+ std::vector>> 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 / max_edge_size_at_feature_edges_ - 0.5) + 1;
+ std::vector> 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