Adds a Delaunay triangulation capability to Quest by kennyweiss · Pull Request #742 · llnl/axom

kennyweiss · 2021-12-07T06:48:41Z

Summary

This PR is a feature
It adds a new Delaunay triangulation capability to quest in 2D and 3D.
The triangulation is built incrementally by inserting successive points and updating to mesh to ensure it satisfies the Delaunay condition: No vertices are in the circumsphere of any triangle/tetrahedron in the mesh. For each inserted point, we remove the triangles/tetrahedra in its "Delaunay cavity" and fill this hole with a new set of triangles/tetrahedra containing the new point (the "Delaunay ball").
The underlying triangulation is encoded using the Indexed mesh data structure with Adjacencies (IA), which is implemented using slam's dynamic sets, relations and maps.
This PR contains several other optimizations and improvements throughout axom in support of this new capability, primarily in core, slam and primal.
- The main new additions are quest::Delaunay.hpp and slam's IAMesh class as well as optimizations and improvements to some geometric primitives in primal, such as Triangle, Tetrahedron and Sphere.
- There is also a driver example delaunay_triangulation which generates a triangle or tetrahedral mesh from a specified number of random points within a rectilinear bounding box.

Credits: This capability was initially implemented by @raineyeh during Summer 2017 and has been sitting in a branch.

Update: Here's some rough performance data on the insertion rate (in vertices per second) while refactoring and improving the code:

Delaunay	2D				3D
	1,000	10,000	100,000	1,000,000	1,000	10,000	100,000	1,000,000
Initial sample (after converting driver to CLI on 12/10/2021)	8K	1.5K	.15K	-	0.4K	.2K	-	-
<lots of optimizations that I didn't measure>
Find element using most recently introduced element (~12/24)	160K	94K	43K	11K	25K	23K	21K	13.5K
Find element using small regular grid of vertex indices (~12/26)	130K	158K	152K	115K	29K	29K	27.5 K	23.5K
Use StackArray in fixVertexNeighborhood() instead of std::vector (~12/27)	150K	177K	168K	126K	30K	32K	30K	25K
Use vector instead of map to fix star neighborhood (~12/28)	165K	190K	180K	134K	33K	31K	32K	27.5K
Refactor fixVertexNeighborhood and findCavity() to use slam relation data without extra copies	193K	255K	230K	165K	50K	48K	49K	38.5K
Final (for this PR)	180K	250K	233K	170K	50K	53K	47K	38.5K

I.e. the performance improved by ~200x.
There are still a few opportunities for optimization:

For finding the initial triangle/tetrahedron: I'm currently using a small grid for the start vertex. Depending on the parameters, this could still involve a lot of mesh traversal depending on the ratio of vertices per grid cell. An octree on the vertices might improve performance at very little additional overhead
For the point insertion: I'm inserting the points in random order. Spatially sorting the vertices could improve the coherence and reduce the average number of invalidated elements per vertex insertion

Update: The following shows the insertion of 100 random points within the unit square:

src/axom/core/tests/core_bit_utilities.hpp

src/axom/core/tests/utils_utilities.hpp

src/axom/core/utilities/BitUtilities.hpp

src/axom/primal/tests/primal_in_sphere.cpp

kennyweiss · 2022-01-07T18:37:25Z

Note: This PR got pretty big as I was optimizing the Delaunay triangulation, which depended on large parts of core, slam, primal and spin. I'll add some comments to highlight the key parts.

kennyweiss · 2022-01-07T18:41:03Z

src/axom/core/Array.hpp

+  Array(ArrayOptions::Uninitialized,
+        IndexType num_elements,
+        IndexType capacity = 0,
+        int allocator_id = axom::detail::getAllocatorID<SPACE>());


Note: axom::Array supports two features that can only be enabled through its constructors -- uninitialized data (via the ArrayOptions::Uninitialized tag) and a custom allocator (via the allocatorID argument) -- so I needed to add some extra overloads to support this.

I also improved the testing of this feature and ensured that when you copy and array with these features, the copy gets them as well.

So the default behavior is to initialize the data and a user can choose to leave it uninitialized via an arg to the ctor? I think it should be the other way around.

We did run into some issues with axom::Arrays of non-trivial type, e.g. axom::Array<BitSet> in ImplicitGrid, from the array data being uninitialized - e.g. assigning a new bitset like array[i] = BitsetType{...} would call the destructor on junk data, which led to some weird segfaults.

I think making uninitialized data opt-in also keeps axom::Array more drop-in compatible with std::vector, since vector's behavior is to default-init all elements.

@kanye-quest I agree 100%

@kanye-quest -- is there a test for this?

I added this because, Delaunay::isValid() computation is using large bitsets, e.g. w/ millions of bit ( one per tetrahedron) and in ImplicitGrid::findCandidates(), the cost of initializing these bits to immediately get overwritten was substantial. I saw a substantial speedup when I disabled the initialization.

kennyweiss · 2022-01-07T18:42:39Z

src/axom/core/tests/core_bit_utilities.hpp

+//
+// SPDX-License-Identifier: (BSD-3-Clause)
+
+#include "gtest/gtest.h"


We didn't have any tests for the bit utility functions: popCount(), trailingZeros() and leadingZeros(), so I added new tests before making changes.

kennyweiss · 2022-01-07T18:43:29Z

src/axom/core/tests/utils_utilities.hpp

 // UNIT TESTS
 //------------------------------------------------------------------------------
-TEST(utils_Utilities, log2)
+TEST(utils_utilities, log2)


These tests weren't being run because the name didn't match the file name and we were filtering them out.

kennyweiss · 2022-01-07T18:46:28Z

src/axom/core/utilities/BitUtilities.hpp

+#include "axom/config.hpp"
+#include "axom/core/Macros.hpp"
 #include "axom/core/Types.hpp"



I added compiler intrinsics for our CPU-based compilers, which are used extensively in slam::BitSet
In some runs of the Delaunay triangulation, axom::popCount() was taking >50% of the runtime (in a release build).

Note: the underlying optimized assembly is only enabled when the code is built in an architecture-aware manner, e.g. via the -march=native flag. I'm planning to add a ticket for us to enable this, e.g. via a new config option

kennyweiss · 2022-01-07T18:49:28Z

src/axom/primal/examples/primal_introduction.cpp

 // almost all our examples are in 3D
 constexpr int in3D = 3;

 // primitives represented by doubles in 3D


While I was editing the primal tests, I modernized them a bit, e.g. type aliases via using instead of typedef and using the initializer-list constructors for primal's primitives instead of make_point, make_vector, ....

kennyweiss · 2022-01-07T18:52:32Z

src/axom/primal/geometry/Sphere.hpp

@@ -6,14 +6,16 @@
 #ifndef AXOM_PRIMAL_SPHERE_HPP_
 #define AXOM_PRIMAL_SPHERE_HPP_



The changes to primal::Sphere change its API and implementation to be based on primal primitives -- e.g. it is based on primal::Point (which always has the right dimensions) instead of a pointer to an array.

kennyweiss · 2022-01-07T18:54:37Z

src/axom/primal/geometry/Vector.hpp

 AXOM_HOST_DEVICE Point<T, NDIMS> operator+(const Vector<T, NDIMS>& V,
                                           const Point<T, NDIMS>& P);

+/*!


These changes are related to #708
I think this completes the set.

kennyweiss · 2022-01-07T18:56:40Z

src/axom/primal/operators/in_sphere.hpp

                      const Point<T, 3>& p2,
-                      const Point<T, 3>& p3)
+                      const Point<T, 3>& p3,
+                      double EPS = 1e-8)


This change significantly improves the efficiency of in_sphere() in 3D by reformulating the 5x5 determinant as a 4x4 determinant.

We have an explicit implementation of 4x4 determinants, but 5x5 determinants are based on a general LU decomposition.

kennyweiss · 2022-01-07T18:58:44Z

src/axom/primal/tests/primal_in_sphere.cpp

+      EXPECT_TRUE(sphere.getOrientation(tet[i]) == primal::ON_BOUNDARY);
+    }
+
+    // check that the orientation matches expectations


These tests ensure that we get the same result for in_sphere() as we do when we compute the circumsphere of a triangle/tet and then call getOrientation()

… size of the BitSet Previously, there was an assertion to check that the index was in range.

This makes it easier to loop over the cells and access the indexed data.

The previous implementation was using an ImplicitGrid. For each of vertex in a mesh with V vertices and E elements, it was getting a Bitset of O(E) bits, even though there are, on average, only O(1) elements that we actually need to check. When V and E are large, e.g. when V is O(E^6), it was taking 10-20x longer to validate the Delaunay complex than to generate it incrementally. The current approach is to use an ImplicitGrid to bootstrap the index and then copy over the precomputing intersections to a UniformGrid. This avoids excessive allocations associated with insertions into a UniformGrid.

…om a Point Other operators were added in #708

This allows us to use the function (which only works when `NDIMS==2`) without specifying the template parameter.

In anticipation of refactoring those functions.

Also adds some tests that check the the function formally satisfies the definition of orientation in terms of normals and dot products, and that `primal::orientation()` is consistent with `Plane::getOrientation()`

And adds some tests.

Also: * slight optimization for 2D area computation using this function. * optimizes some determinant-based computations in Triangle -- `ppedVolume()` and `physToBarycentric()`

Compares results to `Sphere::getOrientation()`.

Uses initializer lists for Points and Vectors as well as some more modern C++ idioms.

…y checks This is more efficient and also uses a different algorithm to validate the Delauanay complex than to construct it.

…` wants the number of zeros to its left

…y validation tests

Adds newlines at the end of files. Co-authored-by: Chris White <white238@llnl.gov>

* Revert change in default initialization for members of `axom::Array` * Use `constexpr` for some compile-time `const` variables

Thanks Arlie and Chris!

dtaller · 2022-01-13T17:26:22Z

src/axom/quest/Delaunay.hpp

  /**
   * \brief Default constructor
-   * \note User need to call initializeBoundary(BoundingBox) before adding points.
+   * \note User must to call initializeBoundary(BoundingBox) before adding points.


I know that I'm not on the axom team and this is a small thing, but you have a typo here "User must to call". Probably you just want to say "User must call", no 'to'.

Thanks @dtaller !

Thanks Danny! I'll take care of it.

kennyweiss added Quest Issues related to Axom's 'quest' component Slam Issues related to Axom's 'slam' component labels Dec 7, 2021

kennyweiss self-assigned this Dec 7, 2021

kennyweiss force-pushed the feature/kweiss/delaunay-revamp branch 2 times, most recently from c50c502 to 1d98f07 Compare December 13, 2021 07:17

kennyweiss force-pushed the feature/kweiss/delaunay-revamp branch 3 times, most recently from db0ad93 to 5573c19 Compare December 24, 2021 06:44

kennyweiss force-pushed the feature/kweiss/delaunay-revamp branch from ac1d500 to cf7e47a Compare January 6, 2022 00:49

white238 reviewed Jan 7, 2022

View reviewed changes

src/axom/core/tests/core_bit_utilities.hpp Outdated Show resolved Hide resolved

white238 reviewed Jan 7, 2022

View reviewed changes

src/axom/core/tests/utils_utilities.hpp Outdated Show resolved Hide resolved

white238 reviewed Jan 7, 2022

View reviewed changes

src/axom/core/utilities/BitUtilities.hpp Show resolved Hide resolved

white238 reviewed Jan 7, 2022

View reviewed changes

src/axom/primal/tests/primal_in_sphere.cpp Show resolved Hide resolved

kennyweiss marked this pull request as ready for review January 7, 2022 18:19

kennyweiss requested review from agcapps, bmhan12, cyrush, publixsubfan, rhornung67 and white238 January 7, 2022 18:19

kennyweiss commented Jan 7, 2022

View reviewed changes

kennyweiss and others added 25 commits January 12, 2022 17:50

It is now ok to call slam::Bitset::test() with an index outside the…

95070cb

… size of the BitSet Previously, there was an assertion to check that the index was in range.

Adds a new overload to ImplicitGrid::getCandidates() from a grid cell

db67485

This makes it easier to loop over the cells and access the indexed data.

Adds missing geometric overload in primal for subtracting a Vector fr…

7cd33f4

…om a Point Other operators were added in #708

Adds default template parameter to primal::Segment::normal()

e9408fc

This allows us to use the function (which only works when `NDIMS==2`) without specifying the template parameter.

Improves testing of primal::orientation

e2ae042

In anticipation of refactoring those functions.

Refactors primal::orientation() and fixes its doxygen comments

184b955

Also adds some tests that check the the function formally satisfies the definition of orientation in terms of normals and dot products, and that `primal::orientation()` is consistent with `Plane::getOrientation()`

Exposes tolerance parameter EPS in primal::orientation()

d8ddf7a

And adds some tests.

Minor cleanup to primal::Plane and associated tests

6d91772

Adds signedArea() function to 2D Triangles

28d481c

Also: * slight optimization for 2D area computation using this function. * optimizes some determinant-based computations in Triangle -- `ppedVolume()` and `physToBarycentric()`

Improves Tetrahedron's signedVolume() implementation and testing

24bff98

Improves testing of primal::in_sphere()

e3c983b

Compares results to `Sphere::getOrientation()`.

Cleans up some primal tests and examples

0986ef4

Uses initializer lists for Points and Vectors as well as some more modern C++ idioms.

Use cached spheres instead of in_sphere operator for Delaunay validit…

1ce15ab

…y checks This is more efficient and also uses a different algorithm to validate the Delauanay complex than to construct it.

Refactors primal::Sphere based on primal::Point and primal::Vector

e6992ba

Bugfix: _BitScanReverse finds the index of highest bit; `leadingZeros…

64c8115

…` wants the number of zeros to its left

Bugfix: Use same index type for Implicit and Uniform Grids in Delauna…

be0416a

…y validation tests

Fix build for nvcc

2d64279

Apply suggestions from code review

271dd36

Adds newlines at the end of files. Co-authored-by: Chris White <white238@llnl.gov>

Updates RELEASE-NOTES

2eabc4e

Applying some PR suggestions

80d21e4

* Revert change in default initialization for members of `axom::Array` * Use `constexpr` for some compile-time `const` variables

Clarification: Adds a comment about which compilers define __GNUC__

368b2c8

Applies PR suggestions from Max

98e9f19

Applies suggestions and corrections from code review

d6ecb6a

Thanks Arlie and Chris!

kennyweiss force-pushed the feature/kweiss/delaunay-revamp branch from 479b0cd to d6ecb6a Compare January 13, 2022 03:10

kennyweiss merged commit aabfcf8 into develop Jan 13, 2022

kennyweiss deleted the feature/kweiss/delaunay-revamp branch January 13, 2022 04:40

dtaller reviewed Jan 13, 2022

View reviewed changes

kennyweiss mentioned this pull request Jun 30, 2022

Reorders input points in Scattered Interpolation query #859

Merged

2 tasks

		@@ -6,14 +6,16 @@
		#ifndef AXOM_PRIMAL_SPHERE_HPP_
		#define AXOM_PRIMAL_SPHERE_HPP_

Conversation

kennyweiss commented Dec 7, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Summary

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

kennyweiss commented Jan 7, 2022

Uh oh!

Choose a reason for hiding this comment

Uh oh!

rhornung67 Jan 7, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

kennyweiss Jan 7, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

kennyweiss Jan 7, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

7 participants

kennyweiss commented Dec 7, 2021 •

edited

Loading

rhornung67 Jan 7, 2022 •

edited

Loading

kennyweiss Jan 7, 2022 •

edited

Loading

kennyweiss Jan 7, 2022 •

edited

Loading