meshfree_interpolation.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2012 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"
27 #include "libmesh/parallel.h"
28 
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_field_variables()=" << this->n_field_variables()
40  << "\n n_source_points()=" << _src_pts.size()
41  << std::endl;
42  }
43 
44 
45 
46  std::ostream& operator << (std::ostream& os, const MeshfreeInterpolation& mfi)
47  {
48  mfi.print_info(os);
49  return os;
50  }
51 
52 
53 
55  {
56  _names.clear();
57  _src_pts.clear();
58  _src_vals.clear();
59  }
60 
61 
62 
63  void MeshfreeInterpolation::add_field_data (const std::vector<std::string> &field_names,
64  const std::vector<Point> &pts,
65  const std::vector<Number> &vals)
66  {
67  libmesh_assert_equal_to (field_names.size()*pts.size(), vals.size());
68 
69  // If we already have field variables, we assume we are appending.
70  // that means the names and ordering better be identical!
71  if (!_names.empty())
72  {
73  if (_names.size() != field_names.size())
74  {
75  libMesh::err << "ERROR: when adding field data to an existing list the\n"
76  << "varaible list must be the same!\n";
77  libmesh_error();
78  }
79  for (unsigned int v=0; v<_names.size(); v++)
80  if (_names[v] != field_names[v])
81  {
82  libMesh::err << "ERROR: when adding field data to an existing list the\n"
83  << "varaible list must be the same!\n";
84  libmesh_error();
85  }
86  }
87 
88  // otherwise copy the names
89  else
90  _names = field_names;
91 
92  // append the data
93  _src_pts.insert (_src_pts.end(),
94  pts.begin(),
95  pts.end());
96 
97  _src_vals.insert (_src_vals.end(),
98  vals.begin(),
99  vals.end());
100 
101  libmesh_assert_equal_to (_src_vals.size(),
102  _src_pts.size()*this->n_field_variables());
103  }
104 
105 
106 
108  {
110  {
111  case SYNC_SOURCES:
112  this->gather_remote_data();
113  break;
114 
115  case INVALID_STRATEGY:
116  libmesh_error();
117  break;
118 
119  default:
120  // no preparation required for other strategies
121  break;
122  }
123  }
124 
125 
126 
128  {
129 #ifndef LIBMESH_HAVE_MPI
130 
131  // no MPI -- no-op
132  return;
133 
134 #else
135 
136  // This function must be run on all processors at once
137  parallel_only();
138 
139  START_LOG ("gather_remote_data()", "MeshfreeInterpolation");
140 
141  std::vector<Real> send_buf, recv_buf;
142 
143  libmesh_assert_equal_to (_src_vals.size(),
144  _src_pts.size()*this->n_field_variables());
145 
146  send_buf.reserve (_src_pts.size()*(3 + this->n_field_variables()));
147 
148  // Everyone packs their data at the same time
149  for (unsigned int p_idx=0, v_idx=0; p_idx<_src_pts.size(); p_idx++)
150  {
151  const Point &pt(_src_pts[p_idx]);
152 
153  send_buf.push_back(pt(0));
154  send_buf.push_back(pt(1));
155  send_buf.push_back(pt(2));
156 
157  for (unsigned int var=0; var<this->n_field_variables(); var++)
158  {
159  libmesh_assert_less (v_idx, _src_vals.size());
160 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
161  send_buf.push_back (_src_vals[v_idx].real());
162  send_buf.push_back (_src_vals[v_idx].imag());
163  v_idx++;
164 
165 #else
166  send_buf.push_back (_src_vals[v_idx++]);
167 #endif
168  }
169  }
170 
171  // Send our data to everyone else. Note that MPI-1 said you could not
172  // use the same buffer in nonblocking sends, but that restriction has
173  // recently been removed.
174  std::vector<Parallel::Request> send_request(libMesh::n_processors()-1);
175 
176  // Use a tag for best practices. In debug mode parallel_only() blocks above
177  // so we can be sure there is no other shenanigarry going on, but in optimized
178  // mode there is no such guarantee - other prcoessors could be somewhere else
179  // completing some other communication, and we don't want to intercept that.
181 
182  for (unsigned int proc=0, cnt=0; proc<libMesh::n_processors(); proc++)
183  if (proc != libMesh::processor_id())
184  CommWorld.send (proc, send_buf, send_request[cnt++], tag);
185 
186  // All data has been sent. Receive remote data in any order
187  for (unsigned int comm_step=0; comm_step<(libMesh::n_processors()-1); comm_step++)
188  {
189  // blocking receive
190  CommWorld.receive (Parallel::any_source, recv_buf, tag);
191 
192  // Add their data to our list
193  Point pt;
194  Number val;
195  std::vector<Real>::const_iterator it=recv_buf.begin();
196  while (it != recv_buf.end())
197  {
198  pt(0) = *it, ++it;
199  pt(1) = *it, ++it;
200  pt(2) = *it, ++it;
201 
202  _src_pts.push_back(pt);
203 
204  for (unsigned int var=0; var<this->n_field_variables(); var++)
205  {
206 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
207  Real re = *it; ++it;
208  Real im = *it; ++it;
209 
210  val = Number(re,im);
211 #else
212  val = *it, ++it;
213 #endif
214  _src_vals.push_back(val);
215  }
216  }
217  }
218 
219  Parallel::wait (send_request);
220 
221  STOP_LOG ("gather_remote_data()", "MeshfreeInterpolation");
222 
223 #endif // LIBMESH_HAVE_MPI
224  }
225 
226 
227 
228  //--------------------------------------------------------------------------------
229  // InverseDistanceInterpolation methods
230  template <unsigned int KDDim>
232  {
233 #ifdef LIBMESH_HAVE_NANOFLANN
234 
235  START_LOG ("construct_kd_tree()", "InverseDistanceInterpolation<>");
236 
237  // Initialize underlying KD tree
238  if (_kd_tree.get() == NULL)
239  _kd_tree.reset (new kd_tree_t (KDDim,
240  _point_list_adaptor,
241  nanoflann::KDTreeSingleIndexAdaptorParams(10 /* max leaf */)));
242 
243  libmesh_assert (_kd_tree.get() != NULL);
244 
245  _kd_tree->buildIndex();
246 
247  STOP_LOG ("construct_kd_tree()", "InverseDistanceInterpolation<>");
248 #endif
249  }
250 
251 
252 
253  template <unsigned int KDDim>
254  void InverseDistanceInterpolation<KDDim>::interpolate_field_data (const std::vector<std::string> &field_names,
255  const std::vector<Point> &tgt_pts,
256  std::vector<Number> &tgt_vals) const
257  {
258  // forcibly initialize, if needed
259 #if LIBMESH_HAVE_NANOFLANN
260  if (_kd_tree.get() == NULL)
262 #endif
263 
264  START_LOG ("interpolate_field_data()", "InverseDistanceInterpolation<>");
265 
266  libmesh_assert_equal_to (field_names.size(), this->n_field_variables());
267 
268  // If we already have field variables, we assume we are appending.
269  // that means the names and ordering better be identical!
270  if (_names.size() != field_names.size())
271  {
272  libMesh::err << "ERROR: when adding field data to an existing list the\n"
273  << "varaible list must be the same!\n";
274  libmesh_error();
275  }
276  for (unsigned int v=0; v<_names.size(); v++)
277  if (_names[v] != field_names[v])
278  {
279  libMesh::err << "ERROR: when adding field data to an existing list the\n"
280  << "varaible list must be the same!\n";
281  libmesh_error();
282  }
283 
284  tgt_vals.resize (tgt_pts.size()*this->n_field_variables());
285 
286  std::vector<Number>::iterator out_it = tgt_vals.begin();
287 
288 #ifdef LIBMESH_HAVE_NANOFLANN
289  {
290  const size_t num_results = std::min((size_t) _n_interp_pts, _src_pts.size());
291 
292  std::vector<size_t> ret_index(num_results);
293  std::vector<Real> ret_dist_sqr(num_results);
294 
295  for (std::vector<Point>::const_iterator tgt_it=tgt_pts.begin();
296  tgt_it != tgt_pts.end(); ++tgt_it)
297  {
298  const Point &tgt(*tgt_it);
299  const Real query_pt[] = { tgt(0), tgt(1), tgt(2) };
300 
301  _kd_tree->knnSearch(&query_pt[0], num_results, &ret_index[0], &ret_dist_sqr[0]);
302 
303  this->interpolate (tgt, ret_index, ret_dist_sqr, out_it);
304 
305  // std::cout << "knnSearch(): num_results=" << num_results << "\n";
306  // for (size_t i=0;i<num_results;i++)
307  // std::cout << "idx[" << i << "]="
308  // << std::setw(6) << ret_index[i]
309  // << "\t dist["<< i << "]=" << ret_dist_sqr[i]
310  // << "\t val[" << std::setw(6) << ret_index[i] << "]=" << _src_vals[ret_index[i]]
311  // << std::endl;
312  // std::cout << "\n";
313 
314  // std::cout << "ival=" << _vals[0] << '\n';
315  }
316  }
317 #else
318 
319  libMesh::err << "ERROR: This functionality requires the library to be configured\n"
320  << "with nanoflann KD-Tree approximate nearest neighbor support!\n"
321  << std::endl;
322  libmesh_error();
323 
324 #endif
325 
326  STOP_LOG ("interpolate_field_data()", "InverseDistanceInterpolation<>");
327  }
328 
329  template <unsigned int KDDim>
331  const std::vector<size_t> &src_indices,
332  const std::vector<Real> &src_dist_sqr,
333  std::vector<Number>::iterator &out_it) const
334  {
335  // We explicitly assume that the input source points are sorted from closest to
336  // farthests. assert that assumption in DEBUG mode.
337 #ifdef DEBUG
338  if (!src_dist_sqr.empty())
339  {
340  Real min_dist = src_dist_sqr.front();
341  std::vector<Real>::const_iterator it = src_dist_sqr.begin();
342 
343  for (++it; it!= src_dist_sqr.end(); ++it)
344  {
345  if (*it < min_dist) libmesh_error();
346  min_dist = *it;
347  }
348  }
349 #endif
350 
351 
352  libmesh_assert_equal_to (src_dist_sqr.size(), src_indices.size());
353 
354 
355  // Compute the interpolation weights & interpolated value
356  const unsigned int n_fv = this->n_field_variables();
357  _vals.resize(n_fv); std::fill (_vals.begin(), _vals.end(), Number(0.));
358 
359  Real tot_weight = 0.;
360 
361  std::vector<Real>::const_iterator src_dist_sqr_it=src_dist_sqr.begin();
362  std::vector<size_t>::const_iterator src_idx_it=src_indices.begin();
363 
364  // Loop over source points
365  while ((src_dist_sqr_it != src_dist_sqr.end()) &&
366  (src_idx_it != src_indices.end()))
367  {
368  libmesh_assert_greater_equal (*src_dist_sqr_it, 0.);
369 
370  const Real
371  dist_sq = std::max(*src_dist_sqr_it, std::numeric_limits<Real>::epsilon()),
372  weight = 1./std::pow(dist_sq, _half_power);
373 
374  tot_weight += weight;
375 
376  const unsigned int src_idx = *src_idx_it;
377 
378  // loop over field variables
379  for (unsigned int v=0; v<n_fv; v++)
380  {
381  libmesh_assert_less (src_idx*n_fv+v, _src_vals.size());
382  _vals[v] += _src_vals[src_idx*n_fv+v]*weight;
383  }
384 
385  ++src_dist_sqr_it;
386  ++src_idx_it;
387  }
388 
389  // don't forget normalizing term & set the output buffer!
390  for (unsigned int v=0; v<n_fv; v++, ++out_it)
391  {
392  _vals[v] /= tot_weight;
393 
394  *out_it = _vals[v];
395  }
396  }
397 
398 
399 
400 // ------------------------------------------------------------
401 // Explicit Instantiations
402  template class InverseDistanceInterpolation<1>;
403  template class InverseDistanceInterpolation<2>;
404  template class InverseDistanceInterpolation<3>;
405 
406 } // namespace libMesh
407 

Site Created By: libMesh Developers
Last modified: January 12 2013 15:00:42 UTC

Hosted By:
SourceForge.net Logo