refactor: MAke TP hierarchical cmaps work for non-uniform refinement

cgoates · cgoates · commit c687fe49ca89 · 2026-02-24T16:25:37.000-07:00
diff --git a/src/param/HierarchicalTPParametricAtlas.cpp b/src/param/HierarchicalTPParametricAtlas.cpp
@@ -48,10 +48,11 @@ ParentPoint HierarchicalTPParametricAtlas::parentPoint( const topology::Vertex&
         const std::array<topology::CellOrEndVertex, 3> fully_unflat_verts = fullyUnflattenVertex( dart_level_cmap, local_d );
         const topology::FullyUnflattenedDart fully_unflat_cell = unflattenFull( *mMap->refinementLevels().at( elem_level ), unrefined_d );
 
+        const size_t dart_component = dartComponentDirection( unflattenFull( dart_level_cmap, local_d ) );
         const size_t ratio = [&](){
             size_t out = 1;
             for( size_t level_ii = elem_level; level_ii < dart_level; level_ii++ )
-                out *= mMap->refinementRatios().at( level_ii );
+                out *= mMap->refinementRatios().at( level_ii ).at( dart_component );
             return out;
         }();
         const Vector6dMax unrefined_lengths = mRefinementLevels.at( elem_level )->parametricLengths( topology::Cell( unrefined_d, mMap->dim() ) );
diff --git a/src/topology/HierarchicalTPCombinatorialMap.cpp b/src/topology/HierarchicalTPCombinatorialMap.cpp
@@ -3,17 +3,16 @@
 #include <IndexOperations.hpp>
 #include <ranges>
 #include <numeric>
+#include <Logging.hpp>
 
 using namespace topology;
 
-namespace topology
-{
-bool checkForNoAncestor( const FullyUnflattenedDart& unflat, const size_t n_darts_per_ancestor )
+bool checkForNoAncestor( const FullyUnflattenedDart& unflat, const SmallVector<size_t, 3>& n_darts_per_ancestor )
 {
     using TPDartPos = TPCombinatorialMap::TPDartPos;
 
     const auto check_dart = [&]( const size_t idx, const bool add_one = false ) {
-        return ( unflat.unflat_darts.at( idx ).id() + ( add_one ? 1 : 0 ) ) % n_darts_per_ancestor != 0;
+        return ( unflat.unflat_darts.at( idx ).id() + ( add_one ? 1 : 0 ) ) % n_darts_per_ancestor.at( idx ) != 0;
     };
     constexpr bool add_one = true;
     if( unflat.unflat_darts.size() == 2 ) // 2D case
@@ -108,14 +107,6 @@ bool checkForNoAncestor( const FullyUnflattenedDart& unflat, const size_t n_dart
     throw std::invalid_argument( "Bad unflattened dart passed to checkForNoAncestor" );
 }
 
-bool checkForNoAncestor( const TPCombinatorialMap& tp_map, const Dart& d, const size_t n_darts_per_ancestor )
-{
-    if( n_darts_per_ancestor == 0 )
-        return true;
-    return checkForNoAncestor( unflattenFull( tp_map, d ), n_darts_per_ancestor );
-}
-}
-
 HierarchicalTPCombinatorialMap::HierarchicalTPCombinatorialMap(
     const std::vector<std::shared_ptr<const TPCombinatorialMap>>& refinement_levels,
     const std::vector<std::vector<Cell>>& leaf_elements )
@@ -136,16 +127,14 @@ HierarchicalTPCombinatorialMap::HierarchicalTPCombinatorialMap(
     for( size_t level = 1; level < refinement_levels.size(); level++ )
     {
         const util::IndexVec higher_sizes = get_tp_lengths( *refinement_levels.at( level ) );
-        mRefinementRatios.push_back( higher_sizes.front() / lower_sizes.front() );
-
-        if( higher_sizes.front() % lower_sizes.front() != 0 )
-            throw std::invalid_argument( "Refinement levels must be uniform refinements of each other." );
-
+        SmallVector<size_t, 3>& ratios = mRefinementRatios.emplace_back();
         for( size_t i = 0; i < lower_sizes.size(); i++ )
         {
-            if( higher_sizes.at( i ) / lower_sizes.at( i ) != mRefinementRatios.back() or higher_sizes.at( i ) % lower_sizes.at( i ) != 0 )
-                throw std::invalid_argument( "Refinement levels must be uniform refinements of each other." );
+            if( higher_sizes.at( i ) % lower_sizes.at( i ) != 0 )
+                throw std::invalid_argument( "Refinement levels must be nested refinements of each other." );
+            ratios.push_back( higher_sizes.at( i ) / lower_sizes.at( i ) );
         }
+
         lower_sizes = higher_sizes;
     }
 
@@ -446,56 +435,65 @@ bool HierarchicalTPCombinatorialMap::iterateDartLineage( const Dart& global_d,
 
     if( dart_level == ancestor_or_descendant_level ) return callback( global_d );
     
+    const FullyUnflattenedDart unflat = unflattenFull( *mRefinementLevels.at( dart_level ), local_d );
     if( dart_level > ancestor_or_descendant_level )
     {
         // Call back on a single ancestor dart
-        const size_t darts_per_ancestor_dart =
+        const SmallVector<size_t, 3> darts_per_ancestor_dart =
             std::reduce( std::next( mRefinementRatios.begin(), ancestor_or_descendant_level ),
                          std::next( mRefinementRatios.begin(), dart_level ),
-                         1, std::multiplies<>() );
-
-        const FullyUnflattenedDart unflat = unflattenFull( *mRefinementLevels.at( dart_level ), local_d );
+                         SmallVector<size_t, 3>( dim(), 1 ),
+                         [&]( const SmallVector<size_t, 3>& a, const SmallVector<size_t, 3>& b ) {
+                             SmallVector<size_t, 3> result( a.size() );
+                             std::transform( a.begin(), a.end(), b.begin(), result.begin(), [&]( const size_t x, const size_t y ) { return x * y; } );
+                             return result;
+                         } );
 
         if( checkForNoAncestor( unflat, darts_per_ancestor_dart ) ) return true;
 
         FullyUnflattenedDart ancestor_unflat( {}, unflat.dart_pos );
-        std::transform(
-            unflat.unflat_darts.begin(),
-            unflat.unflat_darts.end(),
-            std::back_inserter( ancestor_unflat.unflat_darts ),
-            [&darts_per_ancestor_dart]( const Dart& d ) { return Dart( d.id() / darts_per_ancestor_dart ); } );
+        std::transform( unflat.unflat_darts.begin(),
+                        unflat.unflat_darts.end(),
+                        darts_per_ancestor_dart.begin(),
+                        std::back_inserter( ancestor_unflat.unflat_darts ),
+                        []( const Dart& d, const size_t ratio ) {
+                            return Dart( d.id() / ratio );
+                        } );
 
         const Dart ancestor_dart = flattenFull( *mRefinementLevels.at( ancestor_or_descendant_level ), ancestor_unflat );
         return callback( mRanges.toGlobalDart( ancestor_or_descendant_level, ancestor_dart ) );
     }
     else
     {
         // Iterate several descendants.
-        const size_t darts_per_ancestor_dart =
+        const SmallVector<size_t, 3> darts_per_ancestor_dart =
             std::reduce( std::next( mRefinementRatios.begin(), dart_level ),
                          std::next( mRefinementRatios.begin(), ancestor_or_descendant_level ),
-                         1, std::multiplies<>() );
-
-        const FullyUnflattenedDart unflat = unflattenFull( *mRefinementLevels.at( dart_level ), local_d );
+                         SmallVector<size_t, 3>( dim(), 1 ),
+                         [&]( const SmallVector<size_t, 3>& a, const SmallVector<size_t, 3>& b ) {
+                             SmallVector<size_t, 3> result( a.size() );
+                             std::transform( a.begin(), a.end(), b.begin(), result.begin(), [&]( const size_t x, const size_t y ) { return x * y; } );
+                             return result;
+                         } );
 
         FullyUnflattenedDart start_dart( {}, unflat.dart_pos );
         std::transform(
             unflat.unflat_darts.begin(),
             unflat.unflat_darts.end(),
+            darts_per_ancestor_dart.begin(),
             std::back_inserter( start_dart.unflat_darts ),
-            [&darts_per_ancestor_dart]( const Dart& d ) { return Dart( d.id() * darts_per_ancestor_dart ); } );
+            []( const Dart& d, const size_t ratio ) { return Dart( d.id() * ratio ); } );
 
         // Add a series of TP indices to the start_dart, flatten, and call back.
-        const util::IndexVec lengths( dim(), darts_per_ancestor_dart );
-        const SmallVector<std::variant<bool, size_t>, 3> direction = tpDirectionAlongTPDartPos( unflat.dart_pos, lengths );
+        const SmallVector<std::variant<bool, size_t>, 3> direction = tpDirectionAlongTPDartPos( unflat.dart_pos, darts_per_ancestor_dart );
         const util::IndexVec order = [this]() {
             util::IndexVec order( dim() );
             std::iota( order.begin(), order.end(), 0 );
             return order;
         }();
 
         bool continue_iter = true;
-        util::iterateTensorProduct( lengths, order, direction, [&]( const util::IndexVec& iv ){
+        util::iterateTensorProduct( darts_per_ancestor_dart, order, direction, [&]( const util::IndexVec& iv ){
             if( not continue_iter ) return;
             FullyUnflattenedDart descendant_unflat = start_dart;
             std::transform( descendant_unflat.unflat_darts.begin(),
@@ -554,22 +552,21 @@ bool HierarchicalTPCombinatorialMap::iterateChildren( const Cell& local_cell,
     if( mRefinementLevels.size() <= descendant_level ) return true;
 
     // Iterate several descendants.
-    const size_t darts_per_ancestor_dart = mRefinementRatios.at( cell_level );
+    const auto darts_per_ancestor_dart = mRefinementRatios.at( cell_level );
 
     const FullyUnflattenedDart unflat = unflattenFull( *mRefinementLevels.at( cell_level ), local_cell.dart() );
 
     FullyUnflattenedDart start_dart( {}, unflat.dart_pos );
     std::transform(
         unflat.unflat_darts.begin(),
         unflat.unflat_darts.end(),
+        darts_per_ancestor_dart.begin(),
         std::back_inserter( start_dart.unflat_darts ),
-        [&darts_per_ancestor_dart]( const Dart& d ) { return Dart( d.id() * darts_per_ancestor_dart ); } );
+        []( const Dart& d, const size_t ratio ) { return Dart( d.id() * ratio ); } );
 
     // Add a series of TP indices to the start_dart, flatten, and call back.
-    const util::IndexVec lengths( dim(), darts_per_ancestor_dart );
-
     bool continue_iter = true;
-    util::iterateTensorProduct( lengths, [&]( const util::IndexVec& iv ){
+    util::iterateTensorProduct( darts_per_ancestor_dart, [&]( const util::IndexVec& iv ){
         if( not continue_iter ) return;
         FullyUnflattenedDart descendant_unflat = start_dart;
         std::transform( descendant_unflat.unflat_darts.begin(),
diff --git a/src/topology/HierarchicalTPCombinatorialMap.hpp b/src/topology/HierarchicalTPCombinatorialMap.hpp
@@ -43,7 +43,7 @@ namespace topology
 
         bool iterateLeafDescendants( const Dart& global_d, const std::function<bool( const Dart& )>& callback ) const;
 
-        const std::vector<size_t>& refinementRatios() const { return mRefinementRatios; }
+        const std::vector<SmallVector<size_t, 3>>& refinementRatios() const { return mRefinementRatios; }
 
         bool isUnrefinedLeafDart( const Dart& d ) const
         {
@@ -65,7 +65,7 @@ namespace topology
         std::vector<bool> mUnrefinedDarts;
         std::map<Dart, Dart> mPhiOnes;
         std::map<Dart, Dart> mPhiMinusOnes;
-        std::vector<size_t> mRefinementRatios;
+        std::vector<SmallVector<size_t, 3>> mRefinementRatios;
     };
 
     /// @brief A mutable version of HierarchicalTPCombinatorialMap.
@@ -95,8 +95,6 @@ namespace topology
             return mLeafDarts.at( d.id() );
         }
 
-        size_t refinementRatio( const size_t level ) const { return mRefinementRatios.at( level ); }
-
         /// Expose a protected method from the base class for initialization purposes.
         bool iterateDartLineage( const Dart& global_d,
                                  const size_t ancestor_or_descendant_level,
@@ -124,7 +122,5 @@ namespace topology
         }
     };
 
-    bool checkForNoAncestor( const TPCombinatorialMap& tp_map, const Dart& d, const size_t n_darts_per_ancestor );
-
     std::vector<std::vector<Cell>> leafElements( const HierarchicalTPCombinatorialMap& cmap );
 }; // namespace topology
diff --git a/src/topology/TPCombinatorialMap.cpp b/src/topology/TPCombinatorialMap.cpp
@@ -332,6 +332,49 @@ namespace topology
         return out;
     }
 
+    size_t dartComponentDirection( const FullyUnflattenedDart& unflat_d )
+    {
+        if( unflat_d.dart_pos.size() == 1 )
+        {
+            switch( unflat_d.dart_pos.at(0) )
+            {
+                case TPCombinatorialMap::TPDartPos::DartPos0:
+                case TPCombinatorialMap::TPDartPos::DartPos2:
+                    return 0;
+                case TPCombinatorialMap::TPDartPos::DartPos1:
+                case TPCombinatorialMap::TPDartPos::DartPos3:
+                    return 1;
+                default:
+                    throw std::runtime_error( "Invalid TP Dart Position for 2d" );
+            }
+        }
+        else if( unflat_d.dart_pos.size() == 2 )
+        {
+            switch( unflat_d.dart_pos.at(1) )
+            {
+                case TPCombinatorialMap::TPDartPos::DartPos2:
+                case TPCombinatorialMap::TPDartPos::DartPos4:
+                    return 2;
+                default:
+                    switch( unflat_d.dart_pos.at(0) )
+                    {
+                        case TPCombinatorialMap::TPDartPos::DartPos0:
+                        case TPCombinatorialMap::TPDartPos::DartPos2:
+                            return 0;
+                        case TPCombinatorialMap::TPDartPos::DartPos1:
+                        case TPCombinatorialMap::TPDartPos::DartPos3:
+                            return 1;
+                        default:
+                            throw std::runtime_error( "Invalid TP Dart Position for 2d" );
+                    }
+            }
+        }
+        else
+        {
+            throw std::runtime_error( "Invalid TP Dart Position size" );
+        }
+    }
+
     TPCombinatorialMap tensorProductCMapFromComponents( const SmallVector<std::shared_ptr<const CombinatorialMap1d>, 3>& components )
     {
         if( components.size() == 2 )
diff --git a/src/topology/TPCombinatorialMap.hpp b/src/topology/TPCombinatorialMap.hpp
@@ -104,6 +104,8 @@ namespace topology
         SmallVector<TPCombinatorialMap::TPDartPos, 2> dart_pos;
     };
 
+    size_t dartComponentDirection( const FullyUnflattenedDart& unflat_d );
+
     FullyUnflattenedDart unflattenFull( const TPCombinatorialMap& cmap, const Dart& d );
     Dart flattenFull( const TPCombinatorialMap& cmap, const FullyUnflattenedDart& unflat_d );
 
