libMesh
meshfree_interpolation.C
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 // C++ includes
21 #include <iomanip>
22 
23 // Local includes
24 #include "libmesh/point.h"
25 #include "libmesh/meshfree_interpolation.h"
26 #include "libmesh/libmesh_logging.h"
27 #include "libmesh/parallel.h"
28 #include "libmesh/parallel_algebra.h"
29 
30 
31 namespace libMesh
32 {
33 
34 //--------------------------------------------------------------------------------
35 // MeshfreeInterpolation methods
36 void MeshfreeInterpolation::print_info (std::ostream & os) const
37 {
38  os << "MeshfreeInterpolation"
39  << "\n n_source_points()=" << _src_pts.size()
40  << "\n n_field_variables()=" << this->n_field_variables()
41  << "\n";
42 
43  if (this->n_field_variables())
44  {
45  os << " variables = ";
46  for (unsigned int v=0; v<this->n_field_variables(); v++)
47  os << _names[v] << " ";
48  os << std::endl;
49  }
50 
51 }
52 
53 
54 
55 std::ostream & operator << (std::ostream & os, const MeshfreeInterpolation & mfi)
56 {
57  mfi.print_info(os);
58  return os;
59 }
60 
61 
62 
64 {
65  _names.clear();
66  _src_pts.clear();
67  _src_vals.clear();
68 }
69 
70 
71 
72 void MeshfreeInterpolation::add_field_data (const std::vector<std::string> & field_names,
73  const std::vector<Point> & pts,
74  const std::vector<Number> & vals)
75 {
76  libmesh_experimental();
77  libmesh_assert_equal_to (field_names.size()*pts.size(), vals.size());
78 
79  // If we already have field variables, we assume we are appending.
80  // that means the names and ordering better be identical!
81  if (!_names.empty())
82  {
83  if (_names.size() != field_names.size())
84  libmesh_error_msg("ERROR: when adding field data to an existing list the \nvariable list must be the same!");
85 
86  for (std::size_t v=0; v<_names.size(); v++)
87  if (_names[v] != field_names[v])
88  libmesh_error_msg("ERROR: when adding field data to an existing list the \nvariable list must be the same!");
89  }
90 
91  // otherwise copy the names
92  else
93  _names = field_names;
94 
95  // append the data
96  _src_pts.insert (_src_pts.end(),
97  pts.begin(),
98  pts.end());
99 
100  _src_vals.insert (_src_vals.end(),
101  vals.begin(),
102  vals.end());
103 
104  libmesh_assert_equal_to (_src_vals.size(),
105  _src_pts.size()*this->n_field_variables());
106 }
107 
108 
109 
111 {
113  {
114  case SYNC_SOURCES:
115  this->gather_remote_data();
116  break;
117 
118  case INVALID_STRATEGY:
119  libmesh_error_msg("Invalid _parallelization_strategy = " << _parallelization_strategy);
120 
121  default:
122  // no preparation required for other strategies
123  break;
124  }
125 }
126 
127 
128 
130 {
131 #ifndef LIBMESH_HAVE_MPI
132 
133  // no MPI -- no-op
134  return;
135 
136 #else
137 
138  // This function must be run on all processors at once
139  parallel_object_only();
140 
141  LOG_SCOPE ("gather_remote_data()", "MeshfreeInterpolation");
142 
143  this->comm().allgather(_src_pts);
144  this->comm().allgather(_src_vals);
145 
146 #endif // LIBMESH_HAVE_MPI
147 }
148 
149 
150 
151 //--------------------------------------------------------------------------------
152 // InverseDistanceInterpolation methods
153 template <unsigned int KDDim>
155 {
156 #ifdef LIBMESH_HAVE_NANOFLANN
157 
158  LOG_SCOPE ("construct_kd_tree()", "InverseDistanceInterpolation<>");
159 
160  // Initialize underlying KD tree
161  if (_kd_tree.get() == libmesh_nullptr)
162  _kd_tree.reset (new kd_tree_t (KDDim,
163  _point_list_adaptor,
164  nanoflann::KDTreeSingleIndexAdaptorParams(10 /* max leaf */)));
165 
166  libmesh_assert (_kd_tree.get() != libmesh_nullptr);
167 
168  _kd_tree->buildIndex();
169 #endif
170 }
171 
172 
173 
174 template <unsigned int KDDim>
176 {
177 #ifdef LIBMESH_HAVE_NANOFLANN
178  // Delete the KD Tree and start fresh
179  if (_kd_tree.get())
180  _kd_tree.reset (libmesh_nullptr);
181 #endif
182 
183  // Call base class clear method
185 }
186 
187 
188 
189 template <unsigned int KDDim>
190 void InverseDistanceInterpolation<KDDim>::interpolate_field_data (const std::vector<std::string> & field_names,
191  const std::vector<Point> & tgt_pts,
192  std::vector<Number> & tgt_vals) const
193 {
194  libmesh_experimental();
195 
196  // forcibly initialize, if needed
197 #ifdef LIBMESH_HAVE_NANOFLANN
198  if (_kd_tree.get() == libmesh_nullptr)
200 #endif
201 
202  LOG_SCOPE ("interpolate_field_data()", "InverseDistanceInterpolation<>");
203 
204  libmesh_assert_equal_to (field_names.size(), this->n_field_variables());
205 
206  // If we already have field variables, we assume we are appending.
207  // that means the names and ordering better be identical!
208  if (_names.size() != field_names.size())
209  libmesh_error_msg("ERROR: when adding field data to an existing list the \nvariable list must be the same!");
210 
211  for (std::size_t v=0; v<_names.size(); v++)
212  if (_names[v] != field_names[v])
213  libmesh_error_msg("ERROR: when adding field data to an existing list the \nvariable list must be the same!");
214 
215  tgt_vals.resize (tgt_pts.size()*this->n_field_variables());
216 
217 #ifdef LIBMESH_HAVE_NANOFLANN
218  {
219  std::vector<Number>::iterator out_it = tgt_vals.begin();
220 
221  const size_t num_results = std::min((size_t) _n_interp_pts, _src_pts.size());
222 
223  std::vector<size_t> ret_index(num_results);
224  std::vector<Real> ret_dist_sqr(num_results);
225 
226  for (std::vector<Point>::const_iterator tgt_it=tgt_pts.begin();
227  tgt_it != tgt_pts.end(); ++tgt_it)
228  {
229  const Point & tgt(*tgt_it);
230  const Real query_pt[] = { tgt(0), tgt(1), tgt(2) };
231 
232  _kd_tree->knnSearch(&query_pt[0], num_results, &ret_index[0], &ret_dist_sqr[0]);
233 
234  this->interpolate (tgt, ret_index, ret_dist_sqr, out_it);
235 
236  // libMesh::out << "knnSearch(): num_results=" << num_results << "\n";
237  // for (size_t i=0;i<num_results;i++)
238  // libMesh::out << "idx[" << i << "]="
239  // << std::setw(6) << ret_index[i]
240  // << "\t dist["<< i << "]=" << ret_dist_sqr[i]
241  // << "\t val[" << std::setw(6) << ret_index[i] << "]=" << _src_vals[ret_index[i]]
242  // << std::endl;
243  // libMesh::out << "\n";
244 
245  // libMesh::out << "ival=" << _vals[0] << '\n';
246  }
247  }
248 #else
249 
250  libmesh_error_msg("ERROR: This functionality requires the library to be configured with nanoflann support!");
251 
252 #endif
253 }
254 
255 template <unsigned int KDDim>
257  const std::vector<size_t> & src_indices,
258  const std::vector<Real> & src_dist_sqr,
259  std::vector<Number>::iterator & out_it) const
260 {
261  // We explicitly assume that the input source points are sorted from closest to
262  // farthest. assert that assumption in DEBUG mode.
263 #ifdef DEBUG
264  if (!src_dist_sqr.empty())
265  {
266  Real min_dist = src_dist_sqr.front();
267  std::vector<Real>::const_iterator it = src_dist_sqr.begin();
268 
269  for (++it; it!= src_dist_sqr.end(); ++it)
270  {
271  if (*it < min_dist)
272  libmesh_error_msg(*it << " was less than min_dist = " << min_dist);
273 
274  min_dist = *it;
275  }
276  }
277 #endif
278 
279 
280  libmesh_assert_equal_to (src_dist_sqr.size(), src_indices.size());
281 
282 
283  // Compute the interpolation weights & interpolated value
284  const unsigned int n_fv = this->n_field_variables();
285  _vals.resize(n_fv); std::fill (_vals.begin(), _vals.end(), Number(0.));
286 
287  Real tot_weight = 0.;
288 
289  std::vector<Real>::const_iterator src_dist_sqr_it=src_dist_sqr.begin();
290  std::vector<size_t>::const_iterator src_idx_it=src_indices.begin();
291 
292  // Loop over source points
293  while ((src_dist_sqr_it != src_dist_sqr.end()) &&
294  (src_idx_it != src_indices.end()))
295  {
296  libmesh_assert_greater_equal (*src_dist_sqr_it, 0.);
297 
298  const Real
299  dist_sq = std::max(*src_dist_sqr_it, std::numeric_limits<Real>::epsilon()),
300  weight = 1./std::pow(dist_sq, _half_power);
301 
302  tot_weight += weight;
303 
304  const std::size_t src_idx = *src_idx_it;
305 
306  // loop over field variables
307  for (unsigned int v=0; v<n_fv; v++)
308  {
309  libmesh_assert_less (src_idx*n_fv+v, _src_vals.size());
310  _vals[v] += _src_vals[src_idx*n_fv+v]*weight;
311  }
312 
313  ++src_dist_sqr_it;
314  ++src_idx_it;
315  }
316 
317  // don't forget normalizing term & set the output buffer!
318  for (unsigned int v=0; v<n_fv; v++, ++out_it)
319  {
320  _vals[v] /= tot_weight;
321 
322  *out_it = _vals[v];
323  }
324 }
325 
326 
327 
328 // ------------------------------------------------------------
329 // Explicit Instantiations
330 template class InverseDistanceInterpolation<1>;
331 template class InverseDistanceInterpolation<2>;
332 template class InverseDistanceInterpolation<3>;
333 
334 } // namespace libMesh
Inverse distance interpolation.
unsigned int n_field_variables() const
The number of field variables.
const class libmesh_nullptr_t libmesh_nullptr
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.
long double max(long double a, double b)
libmesh_assert(j)
virtual void interpolate_field_data(const std::vector< std::string > &field_names, const std::vector< Point > &tgt_pts, std::vector< Number > &tgt_vals) const libmesh_override
Interpolate source data at target points.
virtual void clear()
Clears all internal data structures and restores to a pristine state.
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.
dof_id_type weight(const MeshBase &mesh, const processor_id_type pid)
Definition: mesh_tools.C:245
Base class to support various mesh-free interpolation methods.
virtual void construct_kd_tree()
Build & initialize the KD tree, if needed.
friend std::ostream & operator<<(std::ostream &os, const MeshfreeInterpolation &mfi)
Same as above, but allows you to also use stream syntax.
virtual void gather_remote_data()
Gathers source points and values that have been added on other processors.
virtual void prepare_for_use()
Prepares data structures for use.
double pow(double a, int b)
virtual void clear() libmesh_override
Clears all internal data structures and restores to a pristine state.
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.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const Parallel::Communicator & comm() const
nanoflann::KDTreeSingleIndexAdaptor< nanoflann::L2_Simple_Adaptor< Real, PointListAdaptor< KDDim > >, PointListAdaptor< KDDim >, KDDim > kd_tree_t
std::vector< std::string > _names
long double min(long double a, double b)
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
void allgather(const T &send, std::vector< T > &recv) const
Take a vector of length this->size(), and fill in recv[processor_id] = the value of send on that proc...