Radial Basis Function interpolation. More...

#include <radial_basis_interpolation.h>

Inheritance diagram for libMesh::RadialBasisInterpolation< KDDim, RBF >:

Public Types
enum	ParallelizationStrategy { SYNC_SOURCES = 0, INVALID_STRATEGY }
	"ParallelizationStrategy" to employ. More...

Public Member Functions
	RadialBasisInterpolation (const libMesh::Parallel::Communicator &comm_in, Real radius=-1)
	Constructor. More...

virtual void	clear () override
	Clears all internal data structures and restores to a pristine state. More...

virtual void	prepare_for_use () override
	Prepares data structures for use. More...

virtual void	interpolate_field_data (const std::vector< std::string > &field_names, const std::vector< Point > &tgt_pts, std::vector< Number > &tgt_vals) const override
	Interpolate source data at target points. More...

void	print_info (std::ostream &os=libMesh::out) const
	Prints information about this object, by default to libMesh::out. More...

unsigned int	n_field_variables () const
	The number of field variables. More...

void	set_field_variables (const std::vector< std::string > &names)
	Defines the field variable(s) we are responsible for, and importantly their assumed ordering. More...

const std::vector< std::string > &	field_variables () const

std::vector< Point > &	get_source_points ()

std::vector< Number > &	get_source_vals ()

virtual void	add_field_data (const std::vector< std::string > &field_names, const std::vector< Point > &pts, const std::vector< Number > &vals)
	Sets source data at specified points. More...

const Parallel::Communicator &	comm () const

processor_id_type	n_processors () const

processor_id_type	processor_id () const

Protected Types
typedef nanoflann::KDTreeSingleIndexAdaptor< nanoflann::L2_Simple_Adaptor< Real, PointListAdaptor< KDDim > >, PointListAdaptor< KDDim >, KDDim >	kd_tree_t

Protected Member Functions
virtual void	construct_kd_tree ()
	Build & initialize the KD tree, if needed. More...

virtual void	interpolate (const Point &pt, const std::vector< size_t > &src_indices, const std::vector< Real > &src_dist_sqr, std::vector< Number >::iterator &out_it) const
	Performs inverse distance interpolation at the input point from the specified points. More...

virtual void	gather_remote_data ()
	Gathers source points and values that have been added on other processors. More...

Protected Attributes
BoundingBox	_src_bbox
	Bounding box for our source points. More...

std::vector< Number >	_weights
	basis coefficients. More...

Real	_r_bbox
	Diagonal of the bounding box. More...

Real	_r_override
	Diagonal override. More...

PointListAdaptor< KDDim >	_point_list_adaptor

std::unique_ptr< kd_tree_t >	_kd_tree

const Real	_half_power

const unsigned int	_n_interp_pts

std::vector< Number >	_vals
	Temporary work array. More...

ParallelizationStrategy	_parallelization_strategy

std::vector< std::string >	_names

std::vector< Point >	_src_pts

std::vector< Number >	_src_vals

const Parallel::Communicator &	_communicator

Detailed Description

template<unsigned int KDDim, class RBF = WendlandRBF<KDDim, 2>>
class libMesh::RadialBasisInterpolation< KDDim, RBF >

Radial Basis Function interpolation.

Author: Benjamin S. Kirk

Date: 2013

Does radial basis function interpolation using Nanoflann.

Definition at line 43 of file radial_basis_interpolation.h.

Member Typedef Documentation

◆ kd_tree_t

template<unsigned int KDDim>

typedef nanoflann::KDTreeSingleIndexAdaptor<nanoflann::L2_Simple_Adaptor<Real, PointListAdaptor<KDDim> >, PointListAdaptor<KDDim>, KDDim> libMesh::InverseDistanceInterpolation< KDDim >::kd_tree_t

protectedinherited

Definition at line 297 of file meshfree_interpolation.h.

Member Enumeration Documentation

◆ ParallelizationStrategy

enum libMesh::MeshfreeInterpolation::ParallelizationStrategy

inherited

"ParallelizationStrategy" to employ.

SYNC_SOURCES assumes that the data added on each processor are independent and relatively small. Calling the prepare_for_use() method with this ParallelizationStrategy will copy remote data from other processors, so all interpolation can be performed locally.

Other ParallelizationStrategy techniques will be implemented as needed.

Enumerator
SYNC_SOURCES
INVALID_STRATEGY

Definition at line 72 of file meshfree_interpolation.h.

72 {SYNC_SOURCES = 0,

73 INVALID_STRATEGY};

Constructor & Destructor Documentation

◆ RadialBasisInterpolation()

template<unsigned int KDDim, class RBF = WendlandRBF<KDDim, 2>>

libMesh::RadialBasisInterpolation< KDDim, RBF >::RadialBasisInterpolation	(	const libMesh::Parallel::Communicator &	comm_in,
		Real	radius = `-1`
	)

inline

Constructor.

Definition at line 79 of file radial_basis_interpolation.h.

                                             :
     InverseDistanceInterpolation<KDDim> (comm_in,8,2),
     _r_bbox(0.),
     _r_override(radius)
   { libmesh_experimental(); }

Member Function Documentation

◆ add_field_data()

void libMesh::MeshfreeInterpolation::add_field_data	(	const std::vector< std::string > &	field_names,
		const std::vector< Point > &	pts,
		const std::vector< Number > &	vals
	)

virtualinherited

Sets source data at specified points.

Definition at line 72 of file meshfree_interpolation.C.

 {
   libmesh_experimental();
   libmesh_assert_equal_to (field_names.size()*pts.size(), vals.size());
  
   // If we already have field variables, we assume we are appending.
   // that means the names and ordering better be identical!
   if (!_names.empty())
     {
       if (_names.size() != field_names.size())
         libmesh_error_msg("ERROR:  when adding field data to an existing list the \nvariable list must be the same!");
  
       for (std::size_t v=0; v<_names.size(); v++)
         if (_names[v] != field_names[v])
           libmesh_error_msg("ERROR:  when adding field data to an existing list the \nvariable list must be the same!");
     }
  
   // otherwise copy the names
   else
     _names = field_names;
  
   // append the data
   _src_pts.insert (_src_pts.end(),
                    pts.begin(),
                    pts.end());
  
   _src_vals.insert (_src_vals.end(),
                     vals.begin(),
                     vals.end());
  
   libmesh_assert_equal_to (_src_vals.size(),
                            _src_pts.size()*this->n_field_variables());
 }

References libMesh::MeshfreeInterpolation::_names, libMesh::MeshfreeInterpolation::_src_pts, libMesh::MeshfreeInterpolation::_src_vals, and libMesh::MeshfreeInterpolation::n_field_variables().

◆ clear()

template<unsigned int KDDim, class RBF >

void libMesh::RadialBasisInterpolation< KDDim, RBF >::clear ( )

overridevirtual

Clears all internal data structures and restores to a pristine state.

Reimplemented from libMesh::InverseDistanceInterpolation< KDDim >.

Definition at line 42 of file radial_basis_interpolation.C.

 {
   // Call base class clear method
   InverseDistanceInterpolation<KDDim>::clear();
 }

References libMesh::InverseDistanceInterpolation< KDDim >::clear().

◆ comm()

const Parallel::Communicator& libMesh::ParallelObject::comm ( ) const

inlineinherited

Returns: A reference to the Parallel::Communicator object used by this mesh.

Definition at line 94 of file parallel_object.h.

95 { return _communicator; }

References libMesh::ParallelObject::_communicator.

◆ construct_kd_tree()

template<unsigned int KDDim>

void libMesh::InverseDistanceInterpolation< KDDim >::construct_kd_tree ( )

protectedvirtualinherited

Build & initialize the KD tree, if needed.

Definition at line 154 of file meshfree_interpolation.C.

 {
 #ifdef LIBMESH_HAVE_NANOFLANN
  
   LOG_SCOPE ("construct_kd_tree()", "InverseDistanceInterpolation<>");
  
   // Initialize underlying KD tree
   if (_kd_tree.get() == nullptr)
     _kd_tree = libmesh_make_unique<kd_tree_t>
       (KDDim,
        _point_list_adaptor,
        nanoflann::KDTreeSingleIndexAdaptorParams(10 /* max leaf */));
  
   libmesh_assert (_kd_tree.get() != nullptr);
  
   _kd_tree->buildIndex();
 #endif
 }

References libMesh::libmesh_assert().

Referenced by libMesh::InverseDistanceInterpolation< KDDim >::interpolate_field_data(), and libMesh::RadialBasisInterpolation< KDDim, RBF >::prepare_for_use().

◆ field_variables()

const std::vector<std::string>& libMesh::MeshfreeInterpolation::field_variables ( ) const

inlineinherited

Returns: The field variables as a read-only reference.

Definition at line 116 of file meshfree_interpolation.h.

117 { return _names; }

References libMesh::MeshfreeInterpolation::_names.

Referenced by libMesh::MeshlessInterpolationFunction::operator()().

◆ gather_remote_data()

void libMesh::MeshfreeInterpolation::gather_remote_data ( )

protectedvirtualinherited

Gathers source points and values that have been added on other processors.

Note: The user is responsible for adding points only once per processor if this method is called. No attempt is made to identify duplicate points.

This method is virtual so that it can be overwritten or extended as required in derived classes.

Definition at line 129 of file meshfree_interpolation.C.

 {
 #ifndef LIBMESH_HAVE_MPI
  
   // no MPI -- no-op
   return;
  
 #else
  
   // This function must be run on all processors at once
   parallel_object_only();
  
   LOG_SCOPE ("gather_remote_data()", "MeshfreeInterpolation");
  
   this->comm().allgather(_src_pts);
   this->comm().allgather(_src_vals);
  
 #endif // LIBMESH_HAVE_MPI
 }

References libMesh::MeshfreeInterpolation::_src_pts, libMesh::MeshfreeInterpolation::_src_vals, and libMesh::ParallelObject::comm().

Referenced by libMesh::MeshfreeInterpolation::prepare_for_use().

◆ get_source_points()

std::vector<Point>& libMesh::MeshfreeInterpolation::get_source_points ( )

inlineinherited

Returns: A writable reference to the point list.

Definition at line 122 of file meshfree_interpolation.h.

123 { return _src_pts; }

References libMesh::MeshfreeInterpolation::_src_pts.

◆ get_source_vals()

std::vector<Number>& libMesh::MeshfreeInterpolation::get_source_vals ( )

inlineinherited

Returns: A writable reference to the point list.

Definition at line 128 of file meshfree_interpolation.h.

129 { return _src_vals; }

References libMesh::MeshfreeInterpolation::_src_vals.

◆ interpolate()

template<unsigned int KDDim>

void libMesh::InverseDistanceInterpolation< KDDim >::interpolate	(	const Point &	pt,
		const std::vector< size_t > &	src_indices,
		const std::vector< Real > &	src_dist_sqr,
		std::vector< Number >::iterator &	out_it
	)		const

protectedvirtualinherited

Performs inverse distance interpolation at the input point from the specified points.

Definition at line 255 of file meshfree_interpolation.C.

 {
   // We explicitly assume that the input source points are sorted from closest to
   // farthest.  assert that assumption in DEBUG mode.
 #ifdef DEBUG
   if (!src_dist_sqr.empty())
     {
       Real min_dist = src_dist_sqr.front();
  
       for (auto i : src_dist_sqr)
         {
           if (i < min_dist)
             libmesh_error_msg(i << " was less than min_dist = " << min_dist);
  
           min_dist = i;
         }
     }
 #endif
  
  
   libmesh_assert_equal_to (src_dist_sqr.size(), src_indices.size());
  
  
   // Compute the interpolation weights & interpolated value
   const unsigned int n_fv = this->n_field_variables();
   _vals.resize(n_fv);  std::fill (_vals.begin(), _vals.end(), Number(0.));
  
   Real tot_weight = 0.;
  
   std::vector<Real>::const_iterator src_dist_sqr_it=src_dist_sqr.begin();
   std::vector<size_t>::const_iterator src_idx_it=src_indices.begin();
  
   // Loop over source points
   while ((src_dist_sqr_it != src_dist_sqr.end()) &&
          (src_idx_it      != src_indices.end()))
     {
       libmesh_assert_greater_equal (*src_dist_sqr_it, 0.);
  
       const Real
         dist_sq = std::max(*src_dist_sqr_it, std::numeric_limits<Real>::epsilon()),
         weight = 1./std::pow(dist_sq, _half_power);
  
       tot_weight += weight;
  
       const std::size_t src_idx = *src_idx_it;
  
       // loop over field variables
       for (unsigned int v=0; v<n_fv; v++)
         {
           libmesh_assert_less (src_idx*n_fv+v, _src_vals.size());
           _vals[v] += _src_vals[src_idx*n_fv+v]*weight;
         }
  
       ++src_dist_sqr_it;
       ++src_idx_it;
     }
  
   // don't forget normalizing term & set the output buffer!
   for (unsigned int v=0; v<n_fv; v++, ++out_it)
     {
       _vals[v] /= tot_weight;
  
       *out_it = _vals[v];
     }
 }

References std::pow(), libMesh::Real, and libMesh::MeshTools::weight().

◆ interpolate_field_data()

template<unsigned int KDDim, class RBF >

void libMesh::RadialBasisInterpolation< KDDim, RBF >::interpolate_field_data	(	const std::vector< std::string > &	field_names,
		const std::vector< Point > &	tgt_pts,
		std::vector< Number > &	tgt_vals
	)		const

overridevirtual

Interpolate source data at target points.

Pure virtual, must be overridden in derived classes.

Reimplemented from libMesh::InverseDistanceInterpolation< KDDim >.

Definition at line 154 of file radial_basis_interpolation.C.

 {
   LOG_SCOPE ("interpolate_field_data()", "RadialBasisInterpolation<>");
  
   libmesh_experimental();
  
   const unsigned int
     n_vars    = this->n_field_variables();
  
   const std::size_t
     n_src_pts = this->_src_pts.size(),
     n_tgt_pts = tgt_pts.size();
  
   libmesh_assert_equal_to (_weights.size(),    this->_src_vals.size());
   libmesh_assert_equal_to (field_names.size(), this->n_field_variables());
  
   // If we already have field variables, we assume we are appending.
   // that means the names and ordering better be identical!
   if (this->_names.size() != field_names.size())
     libmesh_error_msg("ERROR:  when adding field data to an existing list the \nvariable list must be the same!");
  
   for (std::size_t v=0; v<this->_names.size(); v++)
     if (_names[v] != field_names[v])
       libmesh_error_msg("ERROR:  when adding field data to an existing list the \nvariable list must be the same!");
  
  
   RBF rbf(_r_bbox);
  
   tgt_vals.resize (n_tgt_pts*n_vars);  std::fill (tgt_vals.begin(), tgt_vals.end(), Number(0.));
  
   for (std::size_t tgt=0; tgt<n_tgt_pts; tgt++)
     {
       const Point & p (tgt_pts[tgt]);
  
       for (std::size_t i=0; i<n_src_pts; i++)
         {
           const Point & x_i(_src_pts[i]);
           const Real
             r_i   = (p - x_i).norm(),
             phi_i = rbf(r_i);
  
           for (unsigned int var=0; var<n_vars; var++)
             tgt_vals[tgt*n_vars + var] += _weights[i*n_vars + var]*phi_i;
         }
     }
 }

References n_vars, std::norm(), and libMesh::Real.

◆ n_field_variables()

unsigned int libMesh::MeshfreeInterpolation::n_field_variables ( ) const

inlineinherited

The number of field variables.

Definition at line 103 of file meshfree_interpolation.h.

104 { return cast_int<unsigned int>(_names.size()); }

References libMesh::MeshfreeInterpolation::_names.

Referenced by libMesh::MeshfreeInterpolation::add_field_data(), and libMesh::MeshfreeInterpolation::print_info().

◆ n_processors()

processor_id_type libMesh::ParallelObject::n_processors ( ) const

inlineinherited

Returns: The number of processors in the group.

Definition at line 100 of file parallel_object.h.

101 { return cast_int<processor_id_type>(_communicator.size()); }

References libMesh::ParallelObject::_communicator.

◆ prepare_for_use()

template<unsigned int KDDim, class RBF >

void libMesh::RadialBasisInterpolation< KDDim, RBF >::prepare_for_use ( )

overridevirtual

Prepares data structures for use.

Reimplemented from libMesh::MeshfreeInterpolation.

Definition at line 51 of file radial_basis_interpolation.C.

 {
   // Call base class methods for prep
   InverseDistanceInterpolation<KDDim>::prepare_for_use();
   InverseDistanceInterpolation<KDDim>::construct_kd_tree();
  
 #ifndef LIBMESH_HAVE_EIGEN
  
   libmesh_error_msg("ERROR: this functionality presently requires Eigen!");
  
 #else
   LOG_SCOPE ("prepare_for_use()", "RadialBasisInterpolation<>");
  
   // Construct a bounding box for our source points
   _src_bbox.invalidate();
  
   const std::size_t  n_src_pts = this->_src_pts.size();
   const unsigned int n_vars    = this->n_field_variables();
   libmesh_assert_equal_to (this->_src_vals.size(), n_src_pts*this->n_field_variables());
  
   {
     Point
       &p_min(_src_bbox.min()),
       &p_max(_src_bbox.max());
  
     for (std::size_t p=0; p<n_src_pts; p++)
       {
         const Point & p_src(_src_pts[p]);
  
         for (unsigned int d=0; d<LIBMESH_DIM; d++)
           {
             p_min(d) = std::min(p_min(d), p_src(d));
             p_max(d) = std::max(p_max(d), p_src(d));
           }
       }
   }
  
   libMesh::out << "bounding box is \n"
                << _src_bbox.min() << '\n'
                << _src_bbox.max() << std::endl;
  
  
   // Construct the Radial Basis Function, giving it the size of the domain
   if (_r_override < 0)
     _r_bbox = (_src_bbox.max() - _src_bbox.min()).norm();
   else
     _r_bbox = _r_override;
  
   RBF rbf(_r_bbox);
  
   libMesh::out << "bounding box is \n"
                << _src_bbox.min() << '\n'
                << _src_bbox.max() << '\n'
                << "r_bbox = " << _r_bbox << '\n'
                << "rbf(r_bbox/2) = " << rbf(_r_bbox/2) << std::endl;
  
  
   // Construct the projection Matrix
   typedef Eigen::Matrix<Number, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> DynamicMatrix;
   //typedef Eigen::Matrix<Number, Eigen::Dynamic,              1, Eigen::ColMajor> DynamicVector;
  
   DynamicMatrix A(n_src_pts, n_src_pts), x(n_src_pts,n_vars), b(n_src_pts,n_vars);
  
   for (std::size_t i=0; i<n_src_pts; i++)
     {
       const Point & x_i (_src_pts[i]);
  
       // Diagonal
       A(i,i) = rbf(0.);
  
       for (std::size_t j=i+1; j<n_src_pts; j++)
         {
           const Point & x_j (_src_pts[j]);
  
           const Real r_ij = (x_j - x_i).norm();
  
           A(i,j) = A(j,i) = rbf(r_ij);
         }
  
       // set source data
       for (unsigned int var=0; var<n_vars; var++)
         b(i,var) = _src_vals[i*n_vars + var];
     }
  
  
   // Solve the linear system
   x = A.ldlt().solve(b);
   //x = A.fullPivLu().solve(b);
  
   // save  the weights for each variable
   _weights.resize (this->_src_vals.size());
  
   for (std::size_t i=0; i<n_src_pts; i++)
     for (unsigned int var=0; var<n_vars; var++)
       _weights[i*n_vars + var] = x(i,var);
  
 #endif
  
 }

References A, libMesh::InverseDistanceInterpolation< KDDim >::construct_kd_tree(), n_vars, std::norm(), libMesh::out, libMesh::MeshfreeInterpolation::prepare_for_use(), and libMesh::Real.

◆ print_info()

void libMesh::MeshfreeInterpolation::print_info ( std::ostream & os = libMesh::out ) const

inherited

Prints information about this object, by default to libMesh::out.

Definition at line 36 of file meshfree_interpolation.C.

 {
   os << "MeshfreeInterpolation"
      << "\n n_source_points()=" << _src_pts.size()
      << "\n n_field_variables()=" << this->n_field_variables()
      <<  "\n";
  
   if (this->n_field_variables())
     {
       os << "  variables = ";
       for (unsigned int v=0; v<this->n_field_variables(); v++)
         os << _names[v] << " ";
       os << std::endl;
     }
  
 }