meshfree_interpolation.h

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2018 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/auto_ptr.h" // deprecated
#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, (int) 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