RBF coarsening strategies: first selected points are deterministic (#337)

When increasing the number of cores, the selected points are always the same.

Author: David Blom

src/RBFMeshMotionSolver/AdaptiveCoarsening.C

@@ -82,8 +82,38 @@ namespace rbf
     {
         selectedPositions.clear();
 
-        selectedPositions.push_back( 0 );
-        selectedPositions.push_back( 1 );
+        // An initial selection is needed before the greedy algorithm starts
+        // adding points to the selection.
+        // The first point is the point with the largest disp/value
+        // The second point is the point with the largest distance from the
+        // first point.
+
+        {
+            // Find first point: largest value
+            El::DistMatrix<double, El::MC, El::STAR> norms;
+            El::RowTwoNorms( *values, norms );
+            El::Entry<double> locMax = El::MaxAbsLoc( norms );
+            selectedPositions.push_back( locMax.i );
+
+            // Find second point: largest distance from the first point
+            El::DistMatrix<double> distance = *positions;
+            El::DistMatrix<double> tmp;
+            El::Ones( tmp, distance.Height(), distance.Width() );
+
+            for ( int iColumn = 0; iColumn < tmp.Width(); iColumn++ )
+            {
+                El::DistMatrix<double> view;
+                El::View( view, tmp, 0, iColumn, tmp.Height(), 1 );
+                El::Scale( positions->Get( locMax.i, iColumn ), view );
+            }
+
+            El::Axpy( -1, tmp, distance );
+
+            El::RowTwoNorms( distance, norms );
+            locMax = El::MaxAbsLoc( norms );
+            selectedPositions.push_back( locMax.i );
+        }
+
         int maxPoints = std::min( this->maxPoints, positions->Height() );
         int minPoints = std::min( this->minPoints, positions->Height() );
         double error = 0;

src/RBFMeshMotionSolver/RBFCoarsening.C

@@ -261,9 +261,53 @@ namespace rbf
             rbf::vector errorList( positions.rows() );
             selectedPositions.resize( 2 );
 
-            for ( int i = 0; i < selectedPositions.rows(); i++ )
-                selectedPositions( i ) = i;
+            if ( livePointSelection )
+            {
+                // Select the point with the largest displacment
+                int maxDisplacementIndex = -1;
+                ( values.rowwise().norm() ).maxCoeff( &maxDisplacementIndex );
+
+                // Add first point
+                selectedPositions( 0 ) = maxDisplacementIndex;
+            }
+            else
+            {
+                // With unit displacement, first point is point with largest radius from origin and displacement > 0
+                vector rad = positions.rowwise().norm();
+
+                int maxRadiusFromOriginIndex = -1;
+                double maxRadius = -1;
+
+                for ( int i = 0; i < rad.rows(); i++ )
+                {
+                    if ( rad( i ) > maxRadius && values.row( i ).norm() > SMALL )
+                    {
+                        maxRadius = rad( i );
+                        maxRadiusFromOriginIndex = i;
+                    }
+                }
+
+                selectedPositions( 0 ) = maxRadiusFromOriginIndex;
+            }
+
+            // Find point with largest distance from first point
+            vector rad = ( positions - ( matrix::Constant( positions.rows(), 1, 1.0 ) * positions.row( selectedPositions( 0 ) ) ) ).rowwise().norm();
+            int maxRadiusIndex = -1;
+            double maxRadius = -1;
+
+            for ( int i = 0; i < rad.rows(); i++ )
+            {
+                if ( rad( i ) > maxRadius && (rad( i ) < 1.0 - SMALL || rad( i ) > 1.0 + SMALL) )
+                {
+                    maxRadius = rad( i );
+                    maxRadiusIndex = i;
+                }
+            }
+
+            // Add second point
+            selectedPositions( 1 ) = maxRadiusIndex;
 
+            assert( selectedPositions( 0 ) != selectedPositions( 1 ) );
             assert( positions.rows() >= selectedPositions.rows() );
 
             rbf::matrix positionsInterpolationCoarse = positions;

src/RBFMeshMotionSolver/UnitCoarsening.C

@@ -37,8 +37,38 @@ namespace rbf
 
         selectedPositions.clear();
 
-        selectedPositions.push_back( 0 );
-        selectedPositions.push_back( 1 );
+        // An initial selection is needed before the greedy algorithm starts
+        // adding points to the selection.
+        // The first point is the point with the largest radius from the origin.
+        // The second point is the point with the largest distance from the
+        // first point.
+
+        {
+            // Find first point: largest radius from origin
+            El::DistMatrix<double, El::MC, El::STAR> norms;
+            El::RowTwoNorms( *positions, norms );
+            El::Entry<double> locMax = El::MaxAbsLoc( norms );
+            selectedPositions.push_back( locMax.i );
+
+            // Find second point: largest distance from the first point
+            El::DistMatrix<double> distance = *positions;
+            El::DistMatrix<double> tmp;
+            El::Ones( tmp, distance.Height(), distance.Width() );
+
+            for ( int iColumn = 0; iColumn < tmp.Width(); iColumn++ )
+            {
+                El::DistMatrix<double> view;
+                El::View( view, tmp, 0, iColumn, tmp.Height(), 1 );
+                El::Scale( positions->Get( locMax.i, iColumn ), view );
+            }
+
+            El::Axpy( -1, tmp, distance );
+
+            El::RowTwoNorms( distance, norms );
+            locMax = El::MaxAbsLoc( norms );
+            selectedPositions.push_back( locMax.i );
+        }
+
         int maxPoints = std::min( this->maxPoints, positions->Height() );
         int minPoints = std::min( this->minPoints, positions->Height() );