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 : #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/point.h"
27 : #include "libmesh/parallel_object.h"
28 : #ifdef LIBMESH_HAVE_NANOFLANN
29 : # include "libmesh/ignore_warnings.h"
30 : # include "libmesh/nanoflann.hpp"
31 : # include "libmesh/restore_warnings.h"
32 : #endif
33 :
34 : // C++ includes
35 : #include <string>
36 : #include <vector>
37 : #include <memory>
38 :
39 :
40 : namespace libMesh
41 : {
42 :
43 : /**
44 : * Base class to support various mesh-free interpolation methods.
45 : * Such methods can be useful for example to pass data between
46 : * two different domains which share a physical boundary, where
47 : * that boundary may be discretized differently in each domain.
48 : * This is the case for conjugate heat transfer applications where
49 : * the common interface has overlapping but distinct boundary
50 : * discretizations.
51 : *
52 : * \author Benjamin S. Kirk
53 : * \date 2012
54 : * \brief Base class which defines the mesh-free interpolation interface.
55 : */
56 : class MeshfreeInterpolation : public ParallelObject
57 : {
58 : public:
59 :
60 : /**
61 : * "ParallelizationStrategy" to employ.
62 : *
63 : * SYNC_SOURCES assumes that the data added on each processor are
64 : * independent and relatively small. Calling the \p prepare_for_use()
65 : * method with this \p ParallelizationStrategy will copy remote data
66 : * from other processors, so all interpolation can be performed
67 : * locally.
68 : *
69 : * Other \p ParallelizationStrategy techniques will be implemented
70 : * as needed.
71 : */
72 : enum ParallelizationStrategy {SYNC_SOURCES = 0,
73 : INVALID_STRATEGY};
74 : /**
75 : * Constructor.
76 : */
77 0 : MeshfreeInterpolation (const libMesh::Parallel::Communicator & comm_in) :
78 : ParallelObject(comm_in),
79 0 : _parallelization_strategy (SYNC_SOURCES)
80 0 : {}
81 :
82 : /**
83 : * Prints information about this object, by default to
84 : * libMesh::out.
85 : */
86 : void print_info (std::ostream & os=libMesh::out) const;
87 :
88 : /**
89 : * Same as above, but allows you to also use stream syntax.
90 : */
91 : friend std::ostream & operator << (std::ostream & os,
92 : const MeshfreeInterpolation & mfi);
93 :
94 : /**
95 : * Clears all internal data structures and restores to a
96 : * pristine state.
97 : */
98 : virtual void clear();
99 :
100 : /**
101 : * The number of field variables.
102 : */
103 14962 : unsigned int n_field_variables () const
104 21070 : { return cast_int<unsigned int>(_names.size()); }
105 :
106 : /**
107 : * Defines the field variable(s) we are responsible for,
108 : * and importantly their assumed ordering.
109 : */
110 0 : void set_field_variables (std::vector<std::string> names)
111 0 : { _names = std::move(names); }
112 :
113 : /**
114 : *\returns The field variables as a read-only reference.
115 : */
116 5896 : const std::vector<std::string> & field_variables() const
117 73702 : { return _names; }
118 :
119 : /**
120 : * \returns A writable reference to the point list.
121 : */
122 0 : std::vector<Point> & get_source_points ()
123 0 : { return _src_pts; }
124 :
125 : /**
126 : * \returns A writable reference to the point list.
127 : */
128 0 : std::vector<Number> & get_source_vals ()
129 0 : { return _src_vals; }
130 :
131 : /**
132 : * Sets source data at specified points.
133 : */
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 :
138 : /**
139 : * Prepares data structures for use.
140 : *
141 : * This method is virtual so that it can be overwritten or extended as required
142 : * in derived classes.
143 : */
144 : virtual void prepare_for_use ();
145 :
146 : /**
147 : * Interpolate source data at target points.
148 : * Pure virtual, must be overridden in derived classes.
149 : */
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 :
156 : /**
157 : * Gathers source points and values that have been added on other processors.
158 : *
159 : * \note The user is responsible for adding points only once per processor if this
160 : * method is called. No attempt is made to identify duplicate points.
161 : *
162 : * This method is virtual so that it can be overwritten or extended as required
163 : * in derived classes.
164 : */
165 : virtual void gather_remote_data ();
166 :
167 : ParallelizationStrategy _parallelization_strategy;
168 : std::vector<std::string> _names;
169 : std::vector<Point> _src_pts;
170 : std::vector<Number> _src_vals;
171 : };
172 :
173 :
174 :
175 : /**
176 : * Inverse distance interpolation.
177 : */
178 : template <unsigned int KDDim>
179 : class InverseDistanceInterpolation : public MeshfreeInterpolation
180 : {
181 : protected:
182 :
183 : #ifdef LIBMESH_HAVE_NANOFLANN
184 : /**
185 : * This class adapts list of libMesh \p Point types
186 : * for use in a nanoflann KD-Tree. For more on the
187 : * basic idea see examples/pointcloud_adaptor_example.cpp
188 : * in the nanoflann source tree.
189 : */
190 : template <unsigned int PLDim>
191 : class PointListAdaptor
192 : {
193 : private:
194 : const std::vector<Point> & _pts;
195 :
196 : public:
197 0 : PointListAdaptor (const std::vector<Point> & pts) :
198 0 : _pts(pts)
199 0 : {}
200 :
201 : /**
202 : * libMesh \p Point coordinate type
203 : */
204 : typedef Real coord_t;
205 :
206 : /**
207 : * Must return the number of data points
208 : */
209 1988 : inline size_t kdtree_get_point_count() const { return _pts.size(); }
210 :
211 : /**
212 : * \returns The distance between the vector "p1[0:size-1]"
213 : * and the data point with index "idx_p2" stored in the class
214 : */
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 :
255 : /**
256 : * \returns The dim'th component of the idx'th point in the class:
257 : * Since this is inlined and the "dim" argument is typically an immediate value, the
258 : * "if's" are actually solved at compile time.
259 : */
260 1591257 : inline coord_t kdtree_get_pt(const size_t idx, int dim) const
261 : {
262 1591257 : libmesh_assert_less (dim, PLDim);
263 1591257 : libmesh_assert_less (idx, _pts.size());
264 1591257 : libmesh_assert_less (dim, 3);
265 :
266 21541607 : const Point & p(_pts[idx]);
267 :
268 51747767 : if (dim==0) return p(0);
269 35905534 : if (dim==1) return p(1);
270 18621257 : return p(2);
271 : }
272 :
273 : /**
274 : * Optional bounding-box computation.
275 : *
276 : * \returns \p true if the BBOX was already computed by the class
277 : * and returned in \p bb so we can avoid redoing it, or \p false
278 : * to default to a standard bbox computation loop.
279 : *
280 : * Look at bb.size() to find out the expected dimensionality
281 : * (e.g. 2 or 3 for point clouds)
282 : */
283 : template <class BBOX>
284 16 : bool kdtree_get_bbox(BBOX & /* bb */) const { return false; }
285 : };
286 :
287 : PointListAdaptor<KDDim> _point_list_adaptor;
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>>,
297 : PointListAdaptor<KDDim>, KDDim, std::size_t> kd_tree_t;
298 :
299 : mutable std::unique_ptr<kd_tree_t> _kd_tree;
300 :
301 : #endif // LIBMESH_HAVE_NANOFLANN
302 :
303 : /**
304 : * Build & initialize the KD tree, if needed.
305 : */
306 : virtual void construct_kd_tree ();
307 :
308 : /**
309 : * Performs inverse distance interpolation at the input point from
310 : * the specified points.
311 : */
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 :
317 : const Real _half_power;
318 : const unsigned int _n_interp_pts;
319 : const Number _background_value;
320 : const Real _background_eff_dist;
321 :
322 : /**
323 : * Temporary work array. Object level scope to avoid cache thrashing.
324 : */
325 : mutable std::vector<Number> _vals;
326 :
327 : public:
328 :
329 : /**
330 : * Constructor. Takes the inverse distance power,
331 : * which defaults to 2.
332 : *
333 : * In addition to the conventional inverse distance interpolation,
334 : * the user can specify a background value and an effective distance
335 : * for the background value. A background value is corresponding to a
336 : * virtual point that is always at a constant distance (which is set
337 : * by background_eff_dist) from the target point. This is useful when
338 : * the target point is far away from any source points and user wants
339 : * an independent value for this kind of scenarios.
340 : */
341 0 : InverseDistanceInterpolation (const libMesh::Parallel::Communicator & comm_in,
342 : const unsigned int n_interp_pts = 8,
343 : const Real power = 2,
344 : const Number background_value = 0.,
345 : const Real background_eff_dist = -1.0) :
346 : MeshfreeInterpolation(comm_in),
347 : #if LIBMESH_HAVE_NANOFLANN
348 0 : _point_list_adaptor(_src_pts),
349 : #endif
350 0 : _half_power(power/2.0),
351 0 : _n_interp_pts(n_interp_pts),
352 0 : _background_value(background_value),
353 0 : _background_eff_dist(background_eff_dist)
354 0 : {}
355 :
356 : /**
357 : * Clears all internal data structures and restores to a
358 : * pristine state.
359 : */
360 : virtual void clear() override;
361 :
362 : /**
363 : * Interpolate source data at target points.
364 : * Pure virtual, must be overridden in derived classes.
365 : */
366 : virtual void interpolate_field_data (const std::vector<std::string> & field_names,
367 : const std::vector<Point> & tgt_pts,
368 : std::vector<Number> & tgt_vals) const override;
369 : };
370 :
371 : } // namespace libMesh
372 :
373 :
374 : #endif // #define LIBMESH_MESHFREE_INTERPOLATION_H
|
