meshfree_interpolation.h

Go to the documentation of this file.
00001 // The libMesh Finite Element Library.
00002 // Copyright (C) 2002-2012 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
00003 
00004 // This library is free software; you can redistribute it and/or
00005 // modify it under the terms of the GNU Lesser General Public
00006 // License as published by the Free Software Foundation; either
00007 // version 2.1 of the License, or (at your option) any later version.
00008 
00009 // This library is distributed in the hope that it will be useful,
00010 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00011 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00012 // Lesser General Public License for more details.
00013 
00014 // You should have received a copy of the GNU Lesser General Public
00015 // License along with this library; if not, write to the Free Software
00016 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00017 
00018 
00019 
00020 #ifndef LIBMESH_MESHFREE_INTERPOLATION_H
00021 #define LIBMESH_MESHFREE_INTERPOLATION_H
00022 
00023 // Local includes
00024 #include "libmesh/libmesh_config.h"
00025 #include "libmesh/libmesh_common.h"
00026 #include "libmesh/auto_ptr.h"
00027 #include "libmesh/point.h"
00028 #ifdef LIBMESH_HAVE_NANOFLANN
00029 #  include "libmesh/nanoflann.hpp"
00030 #endif
00031 
00032 // C++ includes
00033 #include <string>
00034 #include <vector>
00035 
00036 
00037 
00038 namespace libMesh
00039 {
00040 
00050 class MeshfreeInterpolation
00051 {
00052 public:
00053 
00066   enum ParallelizationStrategy { SYNC_SOURCES     = 0,
00067                                  INVALID_STRATEGY}; 
00071   MeshfreeInterpolation () :
00072     _parallelization_strategy (SYNC_SOURCES)
00073   {}
00074 
00079   void print_info (std::ostream& os=libMesh::out) const;
00080 
00084   friend std::ostream& operator << (std::ostream& os, 
00085                                     const MeshfreeInterpolation& mfi);
00086 
00091   virtual void clear();
00092 
00096   unsigned int n_field_variables () const 
00097   { return libmesh_cast_int<unsigned int>(_names.size()); }
00098   
00103   void set_field_variables (const std::vector<std::string> &names)
00104   { _names = names; }
00105 
00109   const std::vector<std::string> & field_variables() const
00110   { return _names; }
00111   
00115   std::vector<Point> & get_source_points ()
00116   { return _src_pts; }
00117   
00121   std::vector<Number> & get_source_vals ()
00122   { return _src_vals; }
00123 
00127   virtual void add_field_data (const std::vector<std::string> &field_names,
00128                                const std::vector<Point>  &pts,
00129                                const std::vector<Number> &vals);
00130 
00137   virtual void prepare_for_use ();
00138 
00143   virtual void interpolate_field_data (const std::vector<std::string> &field_names,
00144                                        const std::vector<Point>  &tgt_pts,
00145                                        std::vector<Number> &tgt_vals) const = 0;
00146 
00147 protected:
00148 
00157   virtual void gather_remote_data ();
00158   
00159   ParallelizationStrategy  _parallelization_strategy;
00160   std::vector<std::string> _names;
00161   std::vector<Point>       _src_pts;
00162   std::vector<Number>      _src_vals;
00163 };
00164 
00165 
00166 
00170 template <unsigned int KDDim>
00171 class InverseDistanceInterpolation : public MeshfreeInterpolation
00172 {
00173 private:
00174 
00175 #ifdef LIBMESH_HAVE_NANOFLANN  
00176 
00182   template <unsigned int PLDim>
00183   class PointListAdaptor
00184   {    
00185   private:
00186     const std::vector<Point> &_pts;
00187 
00188   public:
00189     PointListAdaptor (const std::vector<Point> &pts) :
00190       _pts(pts)
00191     {}
00192 
00196     typedef Real coord_t;
00197 
00201     inline size_t kdtree_get_point_count() const { return _pts.size(); }
00202     
00207     inline coord_t kdtree_distance(const coord_t *p1, const size_t idx_p2, size_t size) const
00208     {
00209       libmesh_assert_equal_to (size, PLDim);
00210       libmesh_assert_less (idx_p2, _pts.size());
00211 
00212       const Point &p2(_pts[idx_p2]);
00213 
00214       switch (size)
00215         {
00216         case 3:
00217           {
00218             const coord_t d0=p1[0] - p2(0);
00219             const coord_t d1=p1[1] - p2(1);
00220             const coord_t d2=p1[2] - p2(2);
00221 
00222             return d0*d0 + d1*d1 + d2*d2;
00223           }
00224          
00225         case 2:
00226           {
00227             const coord_t d0=p1[0] - p2(0);
00228             const coord_t d1=p1[1] - p2(1);
00229 
00230             return d0*d0 + d1*d1;
00231           }
00232          
00233         case 1:
00234           {
00235             const coord_t d0=p1[0] - p2(0);
00236 
00237             return d0*d0;
00238           }
00239          
00240         default:
00241           libMesh::err << "ERROR: unknown size " << size << std::endl;
00242           libmesh_error();
00243         }
00244 
00245       return -1.;
00246     }
00247     
00253     inline coord_t kdtree_get_pt(const size_t idx, int dim) const
00254     {
00255       libmesh_assert_less (dim, (int) PLDim);
00256       libmesh_assert_less (idx, _pts.size());
00257       libmesh_assert_less (dim, 3);
00258 
00259       const Point &p(_pts[idx]);
00260 
00261       if (dim==0) return p(0);
00262       if (dim==1) return p(1);
00263       return p(2);
00264     }
00265     
00272     template <class BBOX>
00273     bool kdtree_get_bbox(BBOX & /* bb */) const { return false; }
00274   };
00275 
00276   PointListAdaptor<KDDim> _point_list_adaptor;
00277 
00278   // template <int KDDIM>
00279   // class KDTree : public KDTreeSingleIndexAdaptor<L2_Simple_Adaptor<num_t, PointListAdaptor >,
00280   //                                             PointListAdaptor,
00281   //                                             KDDIM>
00282   // {
00283   // };
00284 
00285   typedef nanoflann::KDTreeSingleIndexAdaptor<nanoflann::L2_Simple_Adaptor<Real, PointListAdaptor<KDDim> >,
00286                                               PointListAdaptor<KDDim>, KDDim> kd_tree_t;
00287 
00288   mutable AutoPtr<kd_tree_t> _kd_tree;
00289 
00290 #endif // LIBMESH_HAVE_NANOFLANN
00291 
00295   void construct_kd_tree ();
00296 
00301   void interpolate (const Point               &pt,
00302                     const std::vector<size_t> &src_indices,
00303                     const std::vector<Real>   &src_dist_sqr,
00304                     std::vector<Number>::iterator &out_it) const;
00305   
00306   const int          _half_power;
00307   const unsigned int _n_interp_pts;
00308 
00312   mutable std::vector<Number> _vals;
00313   
00314 public:
00315 
00320   InverseDistanceInterpolation (const unsigned int n_interp_pts = 8,
00321                                 const unsigned int power        = 2) :
00322     MeshfreeInterpolation(),
00323 #if LIBMESH_HAVE_NANOFLANN
00324     _point_list_adaptor(_src_pts),
00325 #endif
00326     _half_power(power/2),
00327     _n_interp_pts(n_interp_pts)
00328   {
00329     libmesh_assert_greater_equal (power, 2);
00330   }
00331 
00336   virtual void clear();
00337 
00342   virtual void interpolate_field_data (const std::vector<std::string> &field_names,
00343                                        const std::vector<Point>  &tgt_pts,
00344                                        std::vector<Number> &tgt_vals) const;
00345 
00346 };
00347 
00348 } // namespace libMesh
00349 
00350 
00351 #endif // #define LIBMESH_MESHFREE_INTERPOLATION_H
Site Created By: libMesh Developers
Last modified: February 05 2013 19:54:47 UTC

Hosted By: