Advanced Multi-Physics (AMP)
On-Line Documentation
DelaunayHelpers.h
Go to the documentation of this file.
1#ifndef included_AMP_DelaunayHelpers
2#define included_AMP_DelaunayHelpers
3
4#include <array>
5#include <stdexcept>
6#include <stdlib.h>
7
8#include "AMP/utils/Array.h"
9#include "AMP/utils/UtilityMacros.h"
10#include "AMP/utils/extended_int.h"
11
12
13namespace AMP::DelaunayHelpers {
14
15
16// clang-format off
17template<int NDIM, class TYPE> struct getETYPE;
18template<> struct getETYPE<1,int> { typedef int ETYPE; };
19template<> struct getETYPE<2,int> { typedef int64_t ETYPE; };
20#ifdef __SIZEOF_INT128__
21DISABLE_WARNINGS
22template<> struct getETYPE<3,int> { typedef __int128 ETYPE; };
23template<> struct getETYPE<4,int> { typedef __int128 ETYPE; };
24ENABLE_WARNINGS
25#else
26template<> struct getETYPE<3,int> { typedef AMP::extended::int128_t ETYPE; };
27template<> struct getETYPE<4,int> { typedef AMP::extended::int128_t ETYPE; };
28#endif
29template<> struct getETYPE<1,double> { typedef long double ETYPE; };
30template<> struct getETYPE<2,double> { typedef long double ETYPE; };
31template<> struct getETYPE<3,double> { typedef long double ETYPE; };
32template<> struct getETYPE<4,double> { typedef long double ETYPE; };
33// clang-format on
34
35
36/********************************************************************
37 * Create triangles neighbors from the triangles *
38 ********************************************************************/
39template<size_t NG>
40std::vector<std::array<int, NG + 1>>
41create_tri_neighbors( const std::vector<std::array<int, NG + 1>> &tri );
42
43
44/********************************************************************
45 * Sort the triangles by index *
46 ********************************************************************/
47template<size_t NG>
48std::vector<int> sortTri( std::vector<std::array<int, NG + 1>> &tri );
49template<size_t NG>
50std::vector<std::array<int, NG + 1>> uniqueTri( const std::vector<std::array<int, NG + 1>> &tri );
51
52
53/********************************************************************
54 * Compute the determinant of a matrix *
55 * Note: The matrix is stored in column-major order *
56 * Note: For the determinant to be exact, we must support a signed *
57 * range of +- N^D, where N is the largest input value and D is *
58 * the size of the matrix. *
59 ********************************************************************/
60template<class TYPE, std::size_t NDIM>
61inline TYPE det( const TYPE *M )
62{
63 if constexpr ( NDIM == 1 ) {
64 return M[0];
65 } else if constexpr ( NDIM == 2 ) {
66 return M[0] * M[3] - M[1] * M[2];
67 } else if constexpr ( NDIM == 3 ) {
68 TYPE det( 0 );
69 det += M[0] * ( M[4] * M[8] - M[7] * M[5] );
70 det -= M[3] * ( M[1] * M[8] - M[7] * M[2] );
71 det += M[6] * ( M[1] * M[5] - M[4] * M[2] );
72 return det;
73 } else if constexpr ( NDIM == 4 ) {
74 TYPE tmp[6];
75 tmp[0] = M[2] * M[7] - M[6] * M[3];
76 tmp[1] = M[2] * M[11] - M[10] * M[3];
77 tmp[2] = M[2] * M[15] - M[14] * M[3];
78 tmp[3] = M[6] * M[11] - M[10] * M[7];
79 tmp[4] = M[6] * M[15] - M[14] * M[7];
80 tmp[5] = M[10] * M[15] - M[14] * M[11];
81 TYPE det( 0 );
82 det += M[0] * ( M[5] * tmp[5] - M[9] * tmp[4] + M[13] * tmp[3] );
83 det -= M[4] * ( M[1] * tmp[5] - M[9] * tmp[2] + M[13] * tmp[1] );
84 det += M[8] * ( M[1] * tmp[4] - M[5] * tmp[2] + M[13] * tmp[0] );
85 det -= M[12] * ( M[1] * tmp[3] - M[5] * tmp[1] + M[9] * tmp[0] );
86 return det;
87 } else {
88 throw std::logic_error( "Not programmed" );
89 }
90}
91template<std::size_t NDIM>
92inline typename getETYPE<NDIM, int>::ETYPE det2( const int *M );
93template<>
94inline int det2<1>( const int *M )
95{
96 return M[0];
97}
98template<>
99inline int64_t det2<2>( const int *M )
100{
101 return int64_t( M[0] ) * int64_t( M[3] ) - int64_t( M[1] ) * int64_t( M[2] );
102}
103template<>
104inline typename getETYPE<3, int>::ETYPE det2<3>( const int *M )
105{
106 using int64 = int64_t;
107 using int128 = typename getETYPE<3, int>::ETYPE;
108 int128 t1 = int128( int64( M[4] ) * int64( M[8] ) - int64( M[7] ) * int64( M[5] ) );
109 int128 t2 = int128( int64( M[1] ) * int64( M[8] ) - int64( M[7] ) * int64( M[2] ) );
110 int128 t3 = int128( int64( M[1] ) * int64( M[5] ) - int64( M[4] ) * int64( M[2] ) );
111 auto det = int128( M[0] ) * t1 + int128( -M[3] ) * t2 + int128( M[6] ) * t3;
112 return det;
113}
114template<std::size_t NDIM>
115inline long double det2( const double *M )
116{
117 long double M2[NDIM * NDIM];
118 for ( size_t i = 0; i < NDIM * NDIM; i++ )
119 M2[i] = M[i];
120 return det<long double, NDIM>( M2 );
121}
122
123
124/********************************************************************
125 * Solve a small dense system *
126 * Note: The matrix is stored in column-major order *
127 * Note: The solve functions do not normalize the inverse: *
128 * x = det * M^-1 * b *
129 * Instead they return the normalization constant det so the user *
130 * can perform the normalization if necessary. This helps to *
131 * preserve accuracy. *
132 ********************************************************************/
133template<class TYPE, std::size_t NDIM>
134inline void solve( const TYPE *M, const TYPE *b, TYPE *x, TYPE &det_M )
135{
136 if constexpr ( NDIM == 1 ) {
137 det_M = M[0];
138 x[0] = b[0];
139 } else if constexpr ( NDIM == 2 ) {
140 det_M = M[0] * M[3] - M[1] * M[2];
141 x[0] = M[3] * b[0] - M[2] * b[1];
142 x[1] = M[0] * b[1] - M[1] * b[0];
143 } else if constexpr ( NDIM == 3 ) {
144 TYPE inv[9];
145 inv[0] = M[4] * M[8] - M[7] * M[5];
146 inv[1] = M[7] * M[2] - M[1] * M[8];
147 inv[2] = M[1] * M[5] - M[4] * M[2];
148 inv[3] = M[6] * M[5] - M[3] * M[8];
149 inv[4] = M[0] * M[8] - M[6] * M[2];
150 inv[5] = M[3] * M[2] - M[0] * M[5];
151 inv[6] = M[3] * M[7] - M[6] * M[4];
152 inv[7] = M[6] * M[1] - M[0] * M[7];
153 inv[8] = M[0] * M[4] - M[3] * M[1];
154 det_M = M[0] * inv[0] + M[3] * inv[1] + M[6] * inv[2];
155 x[0] = inv[0] * b[0] + inv[3] * b[1] + inv[6] * b[2];
156 x[1] = inv[1] * b[0] + inv[4] * b[1] + inv[7] * b[2];
157 x[2] = inv[2] * b[0] + inv[5] * b[1] + inv[8] * b[2];
158 } else if constexpr ( NDIM == 4 ) {
159 TYPE inv[16];
160 TYPE m00 = M[0], m01 = M[4], m02 = M[8], m03 = M[12];
161 TYPE m10 = M[1], m11 = M[5], m12 = M[9], m13 = M[13];
162 TYPE m20 = M[2], m21 = M[6], m22 = M[10], m23 = M[14];
163 TYPE m30 = M[3], m31 = M[7], m32 = M[11], m33 = M[15];
164 det_M = det<TYPE, NDIM>( M );
165 inv[0] = m12 * m23 * m31 - m13 * m22 * m31 + m13 * m21 * m32 - m11 * m23 * m32 -
166 m12 * m21 * m33 + m11 * m22 * m33;
167 inv[1] = m03 * m22 * m31 - m02 * m23 * m31 - m03 * m21 * m32 + m01 * m23 * m32 +
168 m02 * m21 * m33 - m01 * m22 * m33;
169 inv[2] = m02 * m13 * m31 - m03 * m12 * m31 + m03 * m11 * m32 - m01 * m13 * m32 -
170 m02 * m11 * m33 + m01 * m12 * m33;
171 inv[3] = m03 * m12 * m21 - m02 * m13 * m21 - m03 * m11 * m22 + m01 * m13 * m22 +
172 m02 * m11 * m23 - m01 * m12 * m23;
173 inv[4] = m13 * m22 * m30 - m12 * m23 * m30 - m13 * m20 * m32 + m10 * m23 * m32 +
174 m12 * m20 * m33 - m10 * m22 * m33;
175 inv[5] = m02 * m23 * m30 - m03 * m22 * m30 + m03 * m20 * m32 - m00 * m23 * m32 -
176 m02 * m20 * m33 + m00 * m22 * m33;
177 inv[6] = m03 * m12 * m30 - m02 * m13 * m30 - m03 * m10 * m32 + m00 * m13 * m32 +
178 m02 * m10 * m33 - m00 * m12 * m33;
179 inv[7] = m02 * m13 * m20 - m03 * m12 * m20 + m03 * m10 * m22 - m00 * m13 * m22 -
180 m02 * m10 * m23 + m00 * m12 * m23;
181 inv[8] = m11 * m23 * m30 - m13 * m21 * m30 + m13 * m20 * m31 - m10 * m23 * m31 -
182 m11 * m20 * m33 + m10 * m21 * m33;
183 inv[9] = m03 * m21 * m30 - m01 * m23 * m30 - m03 * m20 * m31 + m00 * m23 * m31 +
184 m01 * m20 * m33 - m00 * m21 * m33;
185 inv[10] = m01 * m13 * m30 - m03 * m11 * m30 + m03 * m10 * m31 - m00 * m13 * m31 -
186 m01 * m10 * m33 + m00 * m11 * m33;
187 inv[11] = m03 * m11 * m20 - m01 * m13 * m20 - m03 * m10 * m21 + m00 * m13 * m21 +
188 m01 * m10 * m23 - m00 * m11 * m23;
189 inv[12] = m12 * m21 * m30 - m11 * m22 * m30 - m12 * m20 * m31 + m10 * m22 * m31 +
190 m11 * m20 * m32 - m10 * m21 * m32;
191 inv[13] = m01 * m22 * m30 - m02 * m21 * m30 + m02 * m20 * m31 - m00 * m22 * m31 -
192 m01 * m20 * m32 + m00 * m21 * m32;
193 inv[14] = m02 * m11 * m30 - m01 * m12 * m30 - m02 * m10 * m31 + m00 * m12 * m31 +
194 m01 * m10 * m32 - m00 * m11 * m32;
195 inv[15] = m01 * m12 * m20 - m02 * m11 * m20 + m02 * m10 * m21 - m00 * m12 * m21 -
196 m01 * m10 * m22 + m00 * m11 * m22;
197 x[0] = ( inv[0] * b[0] + inv[1] * b[1] + inv[2] * b[2] + inv[3] * b[3] );
198 x[1] = ( inv[4] * b[0] + inv[5] * b[1] + inv[6] * b[2] + inv[7] * b[3] );
199 x[2] = ( inv[8] * b[0] + inv[9] * b[1] + inv[10] * b[2] + inv[11] * b[3] );
200 x[3] = ( inv[12] * b[0] + inv[13] * b[1] + inv[14] * b[2] + inv[15] * b[3] );
201 }
202}
203
204
205/********************************************************************
206 * Solve a small dense system *
207 * Note: The matrix is stored in column-major order *
208 ********************************************************************/
209void solve( int NDIM, const double *M, const double *b, double *x );
210
211
212/****************************************************************
213 * Function to calculate the inverse of a matrix *
214 ****************************************************************/
215void inverse( int NDIM, const double *M, double *M_inv );
216
217
218/****************************************************************
219 * Function to compute the Barycentric coordinates *
220 * Note: we use exact math until we perform the normalization *
221 * The exact math component requires N^D precision *
222 ****************************************************************/
223template<int NDIM, class TYPE>
224std::array<double, NDIM + 1> computeBarycentric( const std::array<TYPE, NDIM> *x,
225 const std::array<TYPE, NDIM> &xi )
226{
227 // Compute the barycentric coordinates T*L=r-r0
228 // http://en.wikipedia.org/wiki/Barycentric_coordinate_system_(mathematics)
229 using ETYPE = typename getETYPE<NDIM, TYPE>::ETYPE;
230 ETYPE T[NDIM * NDIM];
231 for ( int i = 0; i < NDIM; i++ ) {
232 for ( int j = 0; j < NDIM; j++ )
233 T[j + i * NDIM] = ETYPE( x[i][j] - x[NDIM][j] );
234 }
235 ETYPE r[NDIM];
236 for ( int i = 0; i < NDIM; i++ )
237 r[i] = ETYPE( xi[i] - x[NDIM][i] );
238 ETYPE L2[NDIM + 1], det( 0 );
239 DelaunayHelpers::solve<ETYPE, NDIM>( T, r, L2, det );
240 L2[NDIM] = det;
241 for ( int i = 0; i < NDIM; i++ )
242 L2[NDIM] -= L2[i];
243 // Perform the normalization (will require inexact math)
244 double scale = 1.0 / static_cast<double>( det );
245 std::array<double, NDIM + 1> L = { 0 };
246 for ( int i = 0; i < NDIM + 1; i++ )
247 L[i] = static_cast<double>( L2[i] ) * scale;
248 return L;
249}
250
251
252/********************************************************************
253 * This function calculates the volume of a N-dimensional simplex *
254 * 1 | x1-x4 x2-x4 x3-x4 | *
255 * V = -- | y1-y4 y2-y4 y3-y4 | (3D) *
256 * n! | z1-z4 z2-z4 z3-z4 | *
257 * Note: the sign of the volume depends on the order of the points. *
258 * It will be positive for points stored in a clockwise mannor *
259 * This version is optimized for performance. *
260 * Note: If the volume is zero, then the simplex is invalid *
261 * Eg. a line in 2D or a plane in 3D. *
262 * Note: exact math requires N^D precision *
263 ********************************************************************/
264static constexpr double inv_factorial( int N )
265{
266 double x = 1;
267 for ( int i = 2; i <= N; i++ )
268 x *= i;
269 return 1.0 / x;
270}
271template<int NDIM, class TYPE>
272double calcVolume( const std::array<TYPE, NDIM> *x )
273{
274 if constexpr ( NDIM == 1 )
275 return static_cast<double>( x[1][0] - x[0][0] );
276 TYPE M[NDIM * NDIM];
277 for ( int d = 0; d < NDIM; d++ ) {
278 auto tmp = x[NDIM][d];
279 for ( int j = 0; j < NDIM; j++ )
280 M[d + j * NDIM] = x[j][d] - tmp;
281 }
282 constexpr double C = inv_factorial( NDIM );
283 return C * static_cast<double>( DelaunayHelpers::det2<NDIM>( M ) );
284}
285
286
287/********************************************************************
288 * Helper interfaces to convert arrays *
289 ********************************************************************/
290template<class TYPE, std::size_t NDIM>
291AMP::Array<TYPE> convert( const std::vector<std::array<TYPE, NDIM>> &x )
292{
293 AMP::Array<TYPE> y( NDIM, x.size() );
294 for ( size_t i = 0; i < x.size(); i++ ) {
295 for ( size_t d = 0; d < NDIM; d++ )
296 y( d, i ) = x[i][d];
297 }
298 return y;
299}
300template<class TYPE, std::size_t NDIM>
301std::vector<std::array<TYPE, NDIM>> convert( const AMP::Array<TYPE> &x )
302{
303 AMP_ASSERT( x.size( 0 ) == NDIM );
304 std::vector<std::array<TYPE, NDIM>> y( x.size( 1 ) );
305 for ( size_t i = 0; i < y.size(); i++ ) {
306 for ( size_t d = 0; d < NDIM; d++ )
307 y[i][d] = x( d, i );
308 }
309 return y;
310}
311
312
313} // namespace AMP::DelaunayHelpers
314
315#endif
AMP::Array::size
ARRAY_INLINE const ArraySize & size() const
Return the size of the Array.
Definition Array.h:381
DISABLE_WARNINGS
#define DISABLE_WARNINGS
Re-enable warnings.
Definition UtilityMacros.h:247
ENABLE_WARNINGS
#define ENABLE_WARNINGS
Suppress all warnings.
Definition UtilityMacros.h:248
AMP_ASSERT
#define AMP_ASSERT(EXP)
Assert error.
Definition UtilityMacros.h:114
AMP::DelaunayHelpers::create_tri_neighbors
std::vector< std::array< int, NG+1 > > create_tri_neighbors(const std::vector< std::array< int, NG+1 > > &tri)
AMP::DelaunayHelpers::inverse
void inverse(int NDIM, const double *M, double *M_inv)
AMP::DelaunayHelpers::det2< 1 >
int det2< 1 >(const int *M)
Definition DelaunayHelpers.h:94
AMP::DelaunayHelpers::calcVolume
double calcVolume(const std::array< TYPE, NDIM > *x)
Definition DelaunayHelpers.h:272
AMP::DelaunayHelpers::det
TYPE det(const TYPE *M)
Definition DelaunayHelpers.h:61
AMP::DelaunayHelpers::det2< 2 >
int64_t det2< 2 >(const int *M)
Definition DelaunayHelpers.h:99
AMP::DelaunayHelpers::solve
void solve(const TYPE *M, const TYPE *b, TYPE *x, TYPE &det_M)
Definition DelaunayHelpers.h:134
AMP::DelaunayHelpers::det2
getETYPE< NDIM, int >::ETYPE det2(const int *M)
AMP::DelaunayHelpers::inv_factorial
static constexpr double inv_factorial(int N)
Definition DelaunayHelpers.h:264
AMP::DelaunayHelpers::uniqueTri
std::vector< std::array< int, NG+1 > > uniqueTri(const std::vector< std::array< int, NG+1 > > &tri)
AMP::DelaunayHelpers::convert
AMP::Array< TYPE > convert(const std::vector< std::array< TYPE, NDIM > > &x)
Definition DelaunayHelpers.h:291
AMP::DelaunayHelpers::sortTri
std::vector< int > sortTri(std::vector< std::array< int, NG+1 > > &tri)
AMP::DelaunayHelpers::computeBarycentric
std::array< double, NDIM+1 > computeBarycentric(const std::array< TYPE, NDIM > *x, const std::array< TYPE, NDIM > &xi)
Definition DelaunayHelpers.h:224
AMP::DelaunayHelpers::det2< 3 >
getETYPE< 3, int >::ETYPE det2< 3 >(const int *M)
Definition DelaunayHelpers.h:104

Advanced Multi-Physics (AMP)
Oak Ridge National Laboratory
Idaho National Laboratory
Los Alamos National Laboratory 		This page automatically produced from the
source code by doxygen
Last updated: Tue Mar 10 2026 13:06:41.
Comments on this page