libMesh
meshfree_interpolation.h
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 #ifndef LIBMESH_MESHFREE_INTERPOLATION_H
21 #define LIBMESH_MESHFREE_INTERPOLATION_H
22 
23 // Local includes
24 #include "libmesh/libmesh_config.h"
25 #include "libmesh/libmesh_common.h"
26 #include "libmesh/auto_ptr.h"
27 #include "libmesh/point.h"
28 #include "libmesh/parallel_object.h"
29 #ifdef LIBMESH_HAVE_NANOFLANN
30 # include "libmesh/nanoflann.hpp"
31 #endif
32 
33 // C++ includes
34 #include <string>
35 #include <vector>
36 
37 
38 
39 namespace libMesh
40 {
41 
56 {
57 public:
58 
77  LIBMESH_CAN_DEFAULT_TO_COMMWORLD) :
78  ParallelObject(comm_in),
80  {}
81 
86  void print_info (std::ostream & os=libMesh::out) const;
87 
91  friend std::ostream & operator << (std::ostream & os,
92  const MeshfreeInterpolation & mfi);
93 
98  virtual void clear();
99 
103  unsigned int n_field_variables () const
104  { return cast_int<unsigned int>(_names.size()); }
105 
110  void set_field_variables (const std::vector<std::string> & names)
111  { _names = names; }
112 
116  const std::vector<std::string> & field_variables() const
117  { return _names; }
118 
122  std::vector<Point> & get_source_points ()
123  { return _src_pts; }
124 
128  std::vector<Number> & get_source_vals ()
129  { return _src_vals; }
130 
134  virtual void add_field_data (const std::vector<std::string> & field_names,
135  const std::vector<Point> & pts,
136  const std::vector<Number> & vals);
137 
144  virtual void prepare_for_use ();
145 
150  virtual void interpolate_field_data (const std::vector<std::string> & field_names,
151  const std::vector<Point> & tgt_pts,
152  std::vector<Number> & tgt_vals) const = 0;
153 
154 protected:
155 
165  virtual void gather_remote_data ();
166 
168  std::vector<std::string> _names;
169  std::vector<Point> _src_pts;
170  std::vector<Number> _src_vals;
171 };
172 
173 
174 
178 template <unsigned int KDDim>
180 {
181 protected:
182 
183 #ifdef LIBMESH_HAVE_NANOFLANN
184 
190  template <unsigned int PLDim>
192  {
193  private:
194  const std::vector<Point> & _pts;
195 
196  public:
197  PointListAdaptor (const std::vector<Point> & pts) :
198  _pts(pts)
199  {}
200 
204  typedef Real coord_t;
205 
209  inline size_t kdtree_get_point_count() const { return _pts.size(); }
210 
215  inline coord_t kdtree_distance(const coord_t * p1, const size_t idx_p2, size_t size) const
216  {
217  libmesh_assert_equal_to (size, PLDim);
218  libmesh_assert_less (idx_p2, _pts.size());
219 
220  const Point & p2(_pts[idx_p2]);
221 
222  switch (size)
223  {
224  case 3:
225  {
226  const coord_t d0=p1[0] - p2(0);
227  const coord_t d1=p1[1] - p2(1);
228  const coord_t d2=p1[2] - p2(2);
229 
230  return d0*d0 + d1*d1 + d2*d2;
231  }
232 
233  case 2:
234  {
235  const coord_t d0=p1[0] - p2(0);
236  const coord_t d1=p1[1] - p2(1);
237 
238  return d0*d0 + d1*d1;
239  }
240 
241  case 1:
242  {
243  const coord_t d0=p1[0] - p2(0);
244 
245  return d0*d0;
246  }
247 
248  default:
249  libmesh_error_msg("ERROR: unknown size " << size);
250  }
251 
252  return -1.;
253  }
254 
260  inline coord_t kdtree_get_pt(const size_t idx, int dim) const
261  {
262  libmesh_assert_less (dim, (int) PLDim);
263  libmesh_assert_less (idx, _pts.size());
264  libmesh_assert_less (dim, 3);
265 
266  const Point & p(_pts[idx]);
267 
268  if (dim==0) return p(0);
269  if (dim==1) return p(1);
270  return p(2);
271  }
272 
283  template <class BBOX>
284  bool kdtree_get_bbox(BBOX & /* bb */) const { return false; }
285  };
286 
288 
289  // template <int KDDIM>
290  // class KDTree : public KDTreeSingleIndexAdaptor<L2_Simple_Adaptor<num_t, PointListAdaptor >,
291  // PointListAdaptor,
292  // KDDIM>
293  // {
294  // };
295 
296  typedef nanoflann::KDTreeSingleIndexAdaptor<nanoflann::L2_Simple_Adaptor<Real, PointListAdaptor<KDDim>>,
298 
300 
301 #endif // LIBMESH_HAVE_NANOFLANN
302 
306  virtual void construct_kd_tree ();
307 
312  virtual void interpolate (const Point & pt,
313  const std::vector<size_t> & src_indices,
314  const std::vector<Real> & src_dist_sqr,
315  std::vector<Number>::iterator & out_it) const;
316 
318  const unsigned int _n_interp_pts;
319 
323  mutable std::vector<Number> _vals;
324 
325 public:
326 
332  const unsigned int n_interp_pts = 8,
333  const Real power = 2) :
334  MeshfreeInterpolation(comm_in),
335 #if LIBMESH_HAVE_NANOFLANN
336  _point_list_adaptor(_src_pts),
337 #endif
338  _half_power(power/2.0),
339  _n_interp_pts(n_interp_pts)
340  {}
341 
346  virtual void clear() libmesh_override;
347 
352  virtual void interpolate_field_data (const std::vector<std::string> & field_names,
353  const std::vector<Point> & tgt_pts,
354  std::vector<Number> & tgt_vals) const libmesh_override;
355 };
356 
357 } // namespace libMesh
358 
359 
360 #endif // #define LIBMESH_MESHFREE_INTERPOLATION_H
Inverse distance interpolation.
Encapsulates the MPI_Comm object.
Definition: parallel.h:657
MeshfreeInterpolation(const libMesh::Parallel::Communicator &comm_in LIBMESH_CAN_DEFAULT_TO_COMMWORLD)
Constructor.
std::vector< Number > _vals
Temporary work array.
unsigned int n_field_variables() const
The number of field variables.
const std::vector< std::string > & field_variables() const
unsigned int dim
coord_t kdtree_distance(const coord_t *p1, const size_t idx_p2, size_t size) const
void print_info(std::ostream &os=libMesh::out) const
Prints information about this object, by default to libMesh::out.
The libMesh namespace provides an interface to certain functionality in the library.
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
std::vector< Point > & get_source_points()
virtual void clear()
Clears all internal data structures and restores to a pristine state.
coord_t kdtree_get_pt(const size_t idx, int dim) const
ParallelizationStrategy _parallelization_strategy
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.
Base class to support various mesh-free interpolation methods.
friend std::ostream & operator<<(std::ostream &os, const MeshfreeInterpolation &mfi)
Same as above, but allows you to also use stream syntax.
std::vector< Number > & get_source_vals()
This class adapts list of libMesh Point types for use in a nanoflann KD-Tree.
virtual void gather_remote_data()
Gathers source points and values that have been added on other processors.
bool kdtree_get_bbox(BBOX &) const
Optional bounding-box computation.
This class forms the base class for all other classes that are expected to be implemented in parallel...
virtual void prepare_for_use()
Prepares data structures for use.
InverseDistanceInterpolation(const libMesh::Parallel::Communicator &comm_in, const unsigned int n_interp_pts=8, const Real power=2)
Constructor.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
OStreamProxy out
nanoflann::KDTreeSingleIndexAdaptor< nanoflann::L2_Simple_Adaptor< Real, PointListAdaptor< KDDim > >, PointListAdaptor< KDDim >, KDDim > kd_tree_t
void set_field_variables(const std::vector< std::string > &names)
Defines the field variable(s) we are responsible for, and importantly their assumed ordering...
std::vector< std::string > _names
ParallelizationStrategy
"ParallelizationStrategy" to employ.
size_t kdtree_get_point_count() const
Must return the number of data points.
if(!eq) SETERRQ2(((PetscObject) dm) -> comm, PETSC_ERR_ARG_WRONG,"DM of type %s, not of type %s",((PetscObject) dm) ->type, DMLIBMESH)
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
Interpolate source data at target points.
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
unsigned int idx(const ElemType type, const unsigned int nx, const unsigned int i, const unsigned int j)
A useful inline function which replaces the macros used previously.