meshfree_interpolation.h

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2019 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
 
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
 
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
 
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
 
 
#ifndef LIBMESH_MESHFREE_INTERPOLATION_H
#define LIBMESH_MESHFREE_INTERPOLATION_H
 
// Local includes
#include "libmesh/libmesh_config.h"
#include "libmesh/libmesh_common.h"
#include "libmesh/point.h"
#include "libmesh/parallel_object.h"
#ifdef LIBMESH_HAVE_NANOFLANN
#  include "libmesh/ignore_warnings.h"
#  include "libmesh/nanoflann.hpp"
#  include "libmesh/restore_warnings.h"
#endif
 
// C++ includes
#include <string>
#include <vector>
#include <memory>
 
 
namespace libMesh
{
 
class MeshfreeInterpolation : public ParallelObject
{
public:
 
  enum ParallelizationStrategy {SYNC_SOURCES     = 0,
                                INVALID_STRATEGY};
  MeshfreeInterpolation (const libMesh::Parallel::Communicator & comm_in) :
    ParallelObject(comm_in),
    _parallelization_strategy (SYNC_SOURCES)
  {}
 
  void print_info (std::ostream & os=libMesh::out) const;
 
  friend std::ostream & operator << (std::ostream & os,
                                     const MeshfreeInterpolation & mfi);
 
  virtual void clear();
 
  unsigned int n_field_variables () const
  { return cast_int<unsigned int>(_names.size()); }
 
  void set_field_variables (const std::vector<std::string> & names)
  { _names = names; }
 
  const std::vector<std::string> & field_variables() const
  { return _names; }
 
  std::vector<Point> & get_source_points ()
  { return _src_pts; }
 
  std::vector<Number> & get_source_vals ()
  { return _src_vals; }
 
  virtual void add_field_data (const std::vector<std::string> & field_names,
                               const std::vector<Point>  & pts,
                               const std::vector<Number> & vals);
 
  virtual void prepare_for_use ();
 
  virtual void interpolate_field_data (const std::vector<std::string> & field_names,
                                       const std::vector<Point>  & tgt_pts,
                                       std::vector<Number> & tgt_vals) const = 0;
 
protected:
 
  virtual void gather_remote_data ();
 
  ParallelizationStrategy  _parallelization_strategy;
  std::vector<std::string> _names;
  std::vector<Point>       _src_pts;
  std::vector<Number>      _src_vals;
};
 
 
 
template <unsigned int KDDim>
class InverseDistanceInterpolation : public MeshfreeInterpolation
{
protected:
 
#ifdef LIBMESH_HAVE_NANOFLANN
 
  template <unsigned int PLDim>
  class PointListAdaptor
  {
  private:
    const std::vector<Point> & _pts;
 
  public:
    PointListAdaptor (const std::vector<Point> & pts) :
      _pts(pts)
    {}
 
    typedef Real coord_t;
 
    inline size_t kdtree_get_point_count() const { return _pts.size(); }
 
    inline coord_t kdtree_distance(const coord_t * p1, const size_t idx_p2, size_t size) const
    {
      libmesh_assert_equal_to (size, PLDim);
      libmesh_assert_less (idx_p2, _pts.size());
 
      const Point & p2(_pts[idx_p2]);
 
      switch (size)
        {
        case 3:
          {
            const coord_t d0=p1[0] - p2(0);
            const coord_t d1=p1[1] - p2(1);
            const coord_t d2=p1[2] - p2(2);
 
            return d0*d0 + d1*d1 + d2*d2;
          }
 
        case 2:
          {
            const coord_t d0=p1[0] - p2(0);
            const coord_t d1=p1[1] - p2(1);
 
            return d0*d0 + d1*d1;
          }
 
        case 1:
          {
            const coord_t d0=p1[0] - p2(0);
 
            return d0*d0;
          }
 
        default:
          libmesh_error_msg("ERROR: unknown size " << size);
        }
 
      return -1.;
    }
 
    inline coord_t kdtree_get_pt(const size_t idx, int dim) const
    {
      libmesh_assert_less (dim, PLDim);
      libmesh_assert_less (idx, _pts.size());
      libmesh_assert_less (dim, 3);
 
      const Point & p(_pts[idx]);
 
      if (dim==0) return p(0);
      if (dim==1) return p(1);
      return p(2);
    }
 
    template <class BBOX>
    bool kdtree_get_bbox(BBOX & /* bb */) const { return false; }
  };
 
  PointListAdaptor<KDDim> _point_list_adaptor;
 
  // template <int KDDIM>
  // class KDTree : public KDTreeSingleIndexAdaptor<L2_Simple_Adaptor<num_t, PointListAdaptor >,
  //  PointListAdaptor,
  //  KDDIM>
  // {
  // };
 
  typedef nanoflann::KDTreeSingleIndexAdaptor<nanoflann::L2_Simple_Adaptor<Real, PointListAdaptor<KDDim>>,
                                              PointListAdaptor<KDDim>, KDDim> kd_tree_t;
 
  mutable std::unique_ptr<kd_tree_t> _kd_tree;
 
#endif // LIBMESH_HAVE_NANOFLANN
 
  virtual void construct_kd_tree ();
 
  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;
 
  const Real         _half_power;
  const unsigned int _n_interp_pts;
 
  mutable std::vector<Number> _vals;
 
public:
 
  InverseDistanceInterpolation (const libMesh::Parallel::Communicator & comm_in,
                                const unsigned int n_interp_pts = 8,
                                const Real  power               = 2) :
    MeshfreeInterpolation(comm_in),
#if LIBMESH_HAVE_NANOFLANN
    _point_list_adaptor(_src_pts),
#endif
    _half_power(power/2.0),
    _n_interp_pts(n_interp_pts)
  {}
 
  virtual void clear() override;
 
  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;
};
 
} // namespace libMesh
 
 
#endif // #define LIBMESH_MESHFREE_INTERPOLATION_H