From 77bf64f917af295f221074842ba6bc194d301fe5 Mon Sep 17 00:00:00 2001 From: ScSteffen Date: Tue, 15 Oct 2024 14:29:30 -0400 Subject: [PATCH] changed tesslation implementation from recursive to direct --- include/quadratures/qsphericaltessalation.hpp | 5 +- src/quadratures/qsphericaltessalation.cpp | 132 +++++++++--------- src/solvers/snsolver_hpc.cpp | 1 + 3 files changed, 67 insertions(+), 71 deletions(-) diff --git a/include/quadratures/qsphericaltessalation.hpp b/include/quadratures/qsphericaltessalation.hpp index c7721d3e..02f4aaaa 100644 --- a/include/quadratures/qsphericaltessalation.hpp +++ b/include/quadratures/qsphericaltessalation.hpp @@ -20,9 +20,6 @@ class QSphericalTessalation : public QuadratureBase void SetConnectivity() override; private: - std::vector, 3>> subdivide_triangle( const std::array, 3>& triangle ); - std::vector, 3>> recursive_subdivide( const std::vector, 3>>& triangles, - int order ); std::array compute_centroid( const std::array, 3>& triangle ); std::array map_to_unit_sphere( const std::array& point ); double dot_product( const std::array& v1, const std::array& v2 ); @@ -33,6 +30,8 @@ class QSphericalTessalation : public QuadratureBase const std::vector& weights, std::vector>& full_points, std::vector& full_weights ); + + std::vector, 3>> generate_tessellation( const std::array, 3>& triangle, int order ); }; #endif // QSPHERICALTRIANGLE_H diff --git a/src/quadratures/qsphericaltessalation.cpp b/src/quadratures/qsphericaltessalation.cpp index 678c7438..a454086b 100644 --- a/src/quadratures/qsphericaltessalation.cpp +++ b/src/quadratures/qsphericaltessalation.cpp @@ -14,33 +14,72 @@ QSphericalTessalation::QSphericalTessalation( Config* settings ) : QuadratureBas _supportedDimensions = { 2 }; // Extension to 3d is straightforward } +std::vector, 3>> QSphericalTessalation::generate_tessellation( const std::array, 3>& triangle, + int order ) { + std::vector, 3>> triangles; + + // Vertices of the triangle + std::array A = triangle[0]; + std::array B = triangle[1]; + std::array C = triangle[2]; + + // Generate grid points along edges and inside the triangle + std::vector> grid_points; + for( int i = 0; i <= order; ++i ) { + for( int j = 0; j <= order - i; ++j ) { + double alpha = static_cast( i ) / order; + double beta = static_cast( j ) / order; + double gamma = 1.0 - alpha - beta; + std::array point = { + alpha * A[0] + beta * B[0] + gamma * C[0], alpha * A[1] + beta * B[1] + gamma * C[1], alpha * A[2] + beta * B[2] + gamma * C[2] }; + grid_points.push_back( point ); + } + } + + // Create triangles from the grid points + auto idx = [&]( int i, int j ) -> int { return i * ( order + 1 ) - ( i * ( i - 1 ) ) / 2 + j; }; + + for( int i = 0; i < order; ++i ) { + for( int j = 0; j < order - i; ++j ) { + std::array p1 = grid_points[idx( i, j )]; + std::array p2 = grid_points[idx( i + 1, j )]; + std::array p3 = grid_points[idx( i, j + 1 )]; + triangles.push_back( { p1, p2, p3 } ); + + if( j < order - i - 1 ) { + std::array p4 = grid_points[idx( i + 1, j + 1 )]; + triangles.push_back( { p2, p4, p3 } ); + } + } + } + + return triangles; +} + void QSphericalTessalation::SetPointsAndWeights() { // Define basal triangle vertices - std::array A = { 1.0, 0.0, 0.0 }; - std::array B = { 0.0, 1.0, 0.0 }; - std::array C = { 0.0, 0.0, 1.0 }; - std::array, 3> basal_triangle = { A, B, C }; + std::array A = { 1.0, 0.0, 0.0 }; + std::array B = { 0.0, 1.0, 0.0 }; + std::array C = { 0.0, 0.0, 1.0 }; - // Perform recursive subdivision - std::vector, 3>> final_triangles = recursive_subdivide( { basal_triangle }, _order ); + // Generate the tessellation for the given order + std::vector, 3>> final_triangles = generate_tessellation( { A, B, C }, _order ); - // Compute centroids and map to unit sphere + // Compute centroids of triangles and map them to the unit sphere std::vector> centroids; std::vector> mapped_points; std::vector weights; for( const auto& tri : final_triangles ) { + // Compute the centroid auto centroid = compute_centroid( tri ); auto mapped_centroid = map_to_unit_sphere( centroid ); - centroids.push_back( centroid ); - mapped_points.push_back( mapped_centroid ); + centroids.push_back( mapped_centroid ); - // Map triangle vertices to unit sphere - std::array a = map_to_unit_sphere( tri[0] ); - std::array b = map_to_unit_sphere( tri[1] ); - std::array c = map_to_unit_sphere( tri[2] ); - - // Calculate area using spherical excess + // Compute the area of the spherical triangle + auto a = map_to_unit_sphere( tri[0] ); + auto b = map_to_unit_sphere( tri[1] ); + auto c = map_to_unit_sphere( tri[2] ); double area = spherical_triangle_area( a, b, c ); weights.push_back( area ); } @@ -49,40 +88,36 @@ void QSphericalTessalation::SetPointsAndWeights() { std::vector> full_points; std::vector full_weights; - // Perform reflection and permutation - reflect_and_permute( mapped_points, weights, full_points, full_weights ); + // Perform reflection and permutation (as done in the original code) + reflect_and_permute( centroids, weights, full_points, full_weights ); + // Set the points and weights for the quadrature _nq = full_points.size(); _pointsKarth.resize( _nq ); _pointsSphere.resize( _nq ); _weights.resize( _nq ); + for( size_t i = 0; i < _nq; ++i ) { _pointsKarth[i] = { full_points[i][0], full_points[i][1], full_points[i][2] }; _weights[i] = full_weights[i]; } + double w = 0.0; - // Transform _points to _pointsSphere ==>transform (x,y,z) into (my,phi) for( unsigned idx = 0; idx < _nq; idx++ ) { - _pointsSphere[idx].resize( 3 ); // (my,phi) - _pointsSphere[idx][0] = _pointsKarth[idx][2]; // my = z - _pointsSphere[idx][1] = atan2( _pointsKarth[idx][1], _pointsKarth[idx][0] ); // phi in [-pi,pi] - _pointsSphere[idx][2] = 1.0; // radius r + _pointsSphere[idx].resize( 3 ); + _pointsSphere[idx][0] = _pointsKarth[idx][2]; + _pointsSphere[idx][1] = atan2( _pointsKarth[idx][1], _pointsKarth[idx][0] ); + _pointsSphere[idx][2] = 1.0; - // adapt intervall s.t. phi in [0,2pi] if( _pointsSphere[idx][1] < 0 ) { _pointsSphere[idx][1] = 2 * M_PI + _pointsSphere[idx][1]; } - // std::cout << _pointsKarth[idx][0] << " " << _pointsKarth[idx][1] << " " << _pointsKarth[idx][2] << std::endl; - // std::cout << _pointsSphere[idx][0] << " " << _pointsSphere[idx][1] << " " << _pointsSphere[idx][2] << std::endl; - // std::cout << _weights[idx] << std::endl; w += _weights[idx]; } - // std::cout << w << std::endl; - // exit( 1 ); } -void QSphericalTessalation::SetNq() { _nq = pow( GetOrder(), 4 ); } +void QSphericalTessalation::SetNq() { _nq = 4 * pow( GetOrder(), 2 ); } void QSphericalTessalation::SetConnectivity() { // TODO // Not initialized for this quadrature. @@ -92,45 +127,6 @@ void QSphericalTessalation::SetConnectivity() { // TODO bool QSphericalTessalation::CheckOrder() { return true; } -// Helper function to compute midpoint of two vectors -std::array QSphericalTessalation::midpoint( const std::array& p1, const std::array& p2 ) { - return { ( p1[0] + p2[0] ) / 2, ( p1[1] + p2[1] ) / 2, ( p1[2] + p2[2] ) / 2 }; -} - -// Subdivide a triangle into four smaller triangles -std::vector, 3>> QSphericalTessalation::subdivide_triangle( const std::array, 3>& triangle ) { - std::array A = triangle[0]; - std::array B = triangle[1]; - std::array C = triangle[2]; - - // Compute midpoints of edges - std::array AB_mid = midpoint( A, B ); - std::array BC_mid = midpoint( B, C ); - std::array CA_mid = midpoint( C, A ); - - // Define four new triangles - std::vector, 3>> new_triangles = { - { A, AB_mid, CA_mid }, { AB_mid, B, BC_mid }, { CA_mid, BC_mid, C }, { AB_mid, BC_mid, CA_mid } }; - - return new_triangles; -} - -// Recursively subdivide triangles to the desired order -std::vector, 3>> -QSphericalTessalation::recursive_subdivide( const std::vector, 3>>& triangles, int order ) { - if( order == 1 ) { - return triangles; - } - else { - std::vector, 3>> new_triangles; - for( const auto& triangle : triangles ) { - std::vector, 3>> subdivided = subdivide_triangle( triangle ); - new_triangles.insert( new_triangles.end(), subdivided.begin(), subdivided.end() ); - } - return recursive_subdivide( new_triangles, order - 1 ); - } -} - // Compute the centroid of a triangle std::array QSphericalTessalation::compute_centroid( const std::array, 3>& triangle ) { std::array A = triangle[0]; diff --git a/src/solvers/snsolver_hpc.cpp b/src/solvers/snsolver_hpc.cpp index 89b69347..c361bb51 100644 --- a/src/solvers/snsolver_hpc.cpp +++ b/src/solvers/snsolver_hpc.cpp @@ -1295,6 +1295,7 @@ void SNSolverHPC::PrintVolumeOutput( int idx_iter ) { } } } + void SNSolverHPC::PrepareVolumeOutput() { unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();