Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 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 : // Local includes
21 : #include "libmesh/meshfree_interpolation.h"
22 :
23 : #include "libmesh/int_range.h"
24 : #include "libmesh/libmesh_logging.h"
25 : #include "libmesh/parallel.h"
26 : #include "libmesh/parallel_algebra.h"
27 : #include "libmesh/point.h"
28 :
29 : // C++ includes
30 : #include <iomanip>
31 : #include <memory>
32 :
33 : namespace libMesh
34 : {
35 :
36 : //--------------------------------------------------------------------------------
37 : // MeshfreeInterpolation methods
38 71 : void MeshfreeInterpolation::print_info (std::ostream & os) const
39 : {
40 : os << "MeshfreeInterpolation"
41 71 : << "\n n_source_points()=" << _src_pts.size()
42 69 : << "\n n_field_variables()=" << this->n_field_variables()
43 71 : << "\n";
44 :
45 71 : if (this->n_field_variables())
46 : {
47 71 : os << " variables = ";
48 213 : for (auto v : make_range(this->n_field_variables()))
49 142 : os << _names[v] << " ";
50 2 : os << std::endl;
51 : }
52 :
53 71 : }
54 :
55 :
56 :
57 71 : std::ostream & operator << (std::ostream & os, const MeshfreeInterpolation & mfi)
58 : {
59 71 : mfi.print_info(os);
60 71 : return os;
61 : }
62 :
63 :
64 :
65 0 : void MeshfreeInterpolation::clear ()
66 : {
67 0 : _names.clear();
68 0 : _src_pts.clear();
69 0 : _src_vals.clear();
70 0 : }
71 :
72 :
73 :
74 0 : void MeshfreeInterpolation::add_field_data (const std::vector<std::string> & field_names,
75 : const std::vector<Point> & pts,
76 : const std::vector<Number> & vals)
77 : {
78 : libmesh_experimental();
79 0 : libmesh_assert_equal_to (field_names.size()*pts.size(), vals.size());
80 :
81 : // If we already have field variables, we assume we are appending.
82 : // that means the names and ordering better be identical!
83 0 : if (!_names.empty())
84 : {
85 0 : libmesh_error_msg_if(_names.size() != field_names.size(),
86 : "ERROR: when adding field data to an existing list the \nvariable list must be the same!");
87 :
88 0 : for (auto v : index_range(_names))
89 0 : libmesh_error_msg_if(_names[v] != field_names[v],
90 : "ERROR: when adding field data to an existing list the \nvariable list must be the same!");
91 : }
92 :
93 : // otherwise copy the names
94 : else
95 0 : _names = field_names;
96 :
97 : // append the data
98 0 : _src_pts.insert (_src_pts.end(),
99 : pts.begin(),
100 0 : pts.end());
101 :
102 0 : _src_vals.insert (_src_vals.end(),
103 : vals.begin(),
104 0 : vals.end());
105 :
106 0 : libmesh_assert_equal_to (_src_vals.size(),
107 : _src_pts.size()*this->n_field_variables());
108 0 : }
109 :
110 :
111 :
112 284 : void MeshfreeInterpolation::prepare_for_use ()
113 : {
114 284 : switch (_parallelization_strategy)
115 : {
116 284 : case SYNC_SOURCES:
117 284 : this->gather_remote_data();
118 284 : break;
119 :
120 0 : case INVALID_STRATEGY:
121 0 : libmesh_error_msg("Invalid _parallelization_strategy = " << _parallelization_strategy);
122 :
123 0 : default:
124 : // no preparation required for other strategies
125 0 : break;
126 : }
127 284 : }
128 :
129 :
130 :
131 284 : void MeshfreeInterpolation::gather_remote_data ()
132 : {
133 : #ifndef LIBMESH_HAVE_MPI
134 :
135 : // no MPI -- no-op
136 : return;
137 :
138 : #else
139 :
140 : // This function must be run on all processors at once
141 8 : parallel_object_only();
142 :
143 16 : LOG_SCOPE ("gather_remote_data()", "MeshfreeInterpolation");
144 :
145 284 : this->comm().allgather(_src_pts);
146 284 : this->comm().allgather(_src_vals);
147 :
148 : #endif // LIBMESH_HAVE_MPI
149 284 : }
150 :
151 :
152 :
153 : //--------------------------------------------------------------------------------
154 : // InverseDistanceInterpolation methods
155 : template <unsigned int KDDim>
156 284 : void InverseDistanceInterpolation<KDDim>::construct_kd_tree ()
157 : {
158 : #ifdef LIBMESH_HAVE_NANOFLANN
159 :
160 16 : LOG_SCOPE ("construct_kd_tree()", "InverseDistanceInterpolation<>");
161 :
162 : // Initialize underlying KD tree
163 284 : if (_kd_tree.get() == nullptr)
164 292 : _kd_tree = std::make_unique<kd_tree_t>
165 : (KDDim,
166 284 : _point_list_adaptor,
167 268 : nanoflann::KDTreeSingleIndexAdaptorParams(10 /* max leaf */));
168 :
169 8 : libmesh_assert (_kd_tree.get() != nullptr);
170 :
171 284 : _kd_tree->buildIndex();
172 : #endif
173 284 : }
174 :
175 :
176 :
177 : template <unsigned int KDDim>
178 0 : void InverseDistanceInterpolation<KDDim>::clear()
179 : {
180 : #ifdef LIBMESH_HAVE_NANOFLANN
181 : // Delete the KD Tree and start fresh
182 0 : if (_kd_tree.get())
183 0 : _kd_tree.reset (nullptr);
184 : #endif
185 :
186 : // Call base class clear method
187 0 : MeshfreeInterpolation::clear();
188 0 : }
189 :
190 :
191 :
192 : template <unsigned int KDDim>
193 36922 : void InverseDistanceInterpolation<KDDim>::interpolate_field_data (const std::vector<std::string> & field_names,
194 : const std::vector<Point> & tgt_pts,
195 : std::vector<Number> & tgt_vals) const
196 : {
197 : libmesh_experimental();
198 :
199 : // forcibly initialize, if needed
200 : #ifdef LIBMESH_HAVE_NANOFLANN
201 36922 : if (_kd_tree.get() == nullptr)
202 142 : const_cast<InverseDistanceInterpolation<KDDim> *>(this)->construct_kd_tree();
203 : #endif
204 :
205 5900 : LOG_SCOPE ("interpolate_field_data()", "InverseDistanceInterpolation<>");
206 :
207 2950 : libmesh_assert_equal_to (field_names.size(), this->n_field_variables());
208 :
209 : // If we already have field variables, we assume we are appending.
210 : // that means the names and ordering better be identical!
211 42822 : libmesh_error_msg_if(_names.size() != field_names.size(),
212 : "ERROR: when adding field data to an existing list the \nvariable list must be the same!");
213 :
214 73915 : for (auto v : index_range(_names))
215 36993 : libmesh_error_msg_if(_names[v] != field_names[v],
216 : "ERROR: when adding field data to an existing list the \nvariable list must be the same!");
217 :
218 39872 : tgt_vals.resize (tgt_pts.size()*this->n_field_variables());
219 :
220 : #ifdef LIBMESH_HAVE_NANOFLANN
221 : {
222 36922 : std::vector<Number>::iterator out_it = tgt_vals.begin();
223 :
224 39872 : const size_t num_results = std::min((size_t) _n_interp_pts, _src_pts.size());
225 :
226 39872 : std::vector<size_t> ret_index(num_results);
227 39872 : std::vector<Real> ret_dist_sqr(num_results);
228 :
229 80873 : for (const auto & tgt : tgt_pts)
230 : {
231 43951 : const Real query_pt[] = { tgt(0), tgt(1), tgt(2) };
232 :
233 43951 : _kd_tree->knnSearch(query_pt, num_results, ret_index.data(), ret_dist_sqr.data());
234 :
235 43951 : this->interpolate (tgt, ret_index, ret_dist_sqr, out_it);
236 :
237 : // libMesh::out << "knnSearch(): num_results=" << num_results << "\n";
238 : // for (size_t i=0;i<num_results;i++)
239 : // libMesh::out << "idx[" << i << "]="
240 : // << std::setw(6) << ret_index[i]
241 : // << "\t dist["<< i << "]=" << ret_dist_sqr[i]
242 : // << "\t val[" << std::setw(6) << ret_index[i] << "]=" << _src_vals[ret_index[i]]
243 : // << std::endl;
244 : // libMesh::out << "\n";
245 :
246 : // libMesh::out << "ival=" << _vals[0] << '\n';
247 : }
248 : }
249 : #else
250 :
251 : libmesh_error_msg("ERROR: This functionality requires the library to be configured with nanoflann support!");
252 :
253 : #endif
254 36922 : }
255 :
256 : template <unsigned int KDDim>
257 43951 : void InverseDistanceInterpolation<KDDim>::interpolate (const Point & /* pt */,
258 : const std::vector<size_t> & src_indices,
259 : const std::vector<Real> & src_dist_sqr,
260 : std::vector<Number>::iterator & out_it) const
261 : {
262 : // We explicitly assume that the input source points are sorted from closest to
263 : // farthest. assert that assumption in DEBUG mode.
264 : #ifdef DEBUG
265 3148 : if (!src_dist_sqr.empty())
266 : {
267 3148 : Real min_dist = src_dist_sqr.front();
268 :
269 16540 : for (auto i : src_dist_sqr)
270 : {
271 13392 : libmesh_error_msg_if(i < min_dist, i << " was less than min_dist = " << min_dist);
272 :
273 13392 : min_dist = i;
274 : }
275 : }
276 : #endif
277 :
278 :
279 3148 : libmesh_assert_equal_to (src_dist_sqr.size(), src_indices.size());
280 :
281 :
282 : // Compute the interpolation weights & interpolated value
283 3148 : const unsigned int n_fv = this->n_field_variables();
284 43951 : _vals.resize(n_fv);
285 :
286 3148 : std::fill (_vals.begin(), _vals.end(), Number(0.));
287 3148 : Real tot_weight(0.);
288 : // The background value is optional
289 : // If background value option is enabled, add it to the total weight
290 : // If not, a "zero" weight is added
291 43951 : if (_background_eff_dist > 0.0)
292 : {
293 0 : const Real background_wt = _background_eff_dist * _background_eff_dist > std::numeric_limits<Real>::epsilon()
294 0 : ? 1.0 / std::pow(_background_eff_dist * _background_eff_dist, _half_power)
295 : : 0.0;
296 0 : tot_weight += background_wt;
297 0 : std::fill (_vals.begin(), _vals.end(), Number(_background_value * background_wt));
298 : }
299 :
300 6296 : std::vector<Real>::const_iterator src_dist_sqr_it=src_dist_sqr.begin();
301 6296 : std::vector<size_t>::const_iterator src_idx_it=src_indices.begin();
302 :
303 : // Loop over source points
304 261547 : while ((src_dist_sqr_it != src_dist_sqr.end()) &&
305 43324 : (src_idx_it != src_indices.end()))
306 : {
307 13392 : libmesh_assert_greater_equal (*src_dist_sqr_it, 0.);
308 :
309 : const Real
310 204204 : dist_sq = std::max(*src_dist_sqr_it, std::numeric_limits<Real>::epsilon()),
311 204204 : weight = 1./std::pow(dist_sq, _half_power);
312 :
313 204204 : tot_weight += weight;
314 :
315 204204 : const std::size_t src_idx = *src_idx_it;
316 :
317 : // loop over field variables
318 465208 : for (unsigned int v=0; v<n_fv; v++)
319 : {
320 14992 : libmesh_assert_less (src_idx*n_fv+v, _src_vals.size());
321 290988 : _vals[v] += _src_vals[src_idx*n_fv+v]*weight;
322 : }
323 :
324 13392 : ++src_dist_sqr_it;
325 13392 : ++src_idx_it;
326 : }
327 :
328 : // don't forget normalizing term & set the output buffer!
329 95002 : for (unsigned int v=0; v<n_fv; v++, ++out_it)
330 : {
331 51051 : _vals[v] /= tot_weight;
332 :
333 51051 : *out_it = _vals[v];
334 : }
335 43951 : }
336 :
337 :
338 :
339 : // ------------------------------------------------------------
340 : // Explicit Instantiations
341 : template class LIBMESH_EXPORT InverseDistanceInterpolation<1>;
342 : template class LIBMESH_EXPORT InverseDistanceInterpolation<2>;
343 : template class LIBMESH_EXPORT InverseDistanceInterpolation<3>;
344 :
345 : } // namespace libMesh
|
