PLaSK library

Loading...
Searching...
No Matches
interpolation.hpp
Go to the documentation of this file.
1/*
2 * This file is part of PLaSK (https://plask.app) by Photonics Group at TUL
3 * Copyright (c) 2022 Lodz University of Technology
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, version 3.
8 *
9 * This program 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
12 * GNU General Public License for more details.
13 */
14#ifndef PLASK__INTERPOLATION_H
15#define PLASK__INTERPOLATION_H
16
85#include <typeinfo> // for 'typeid'
86
87#include "mesh.hpp"
88#include "../exceptions.hpp"
89#include "../memory.hpp"
90#include "../lazydata.hpp"
91#include "../log/log.hpp"
92#include <plask/geometry/space.hpp>
93
94 namespace plask {
95
100 enum InterpolationMethod: unsigned {
101 INTERPOLATION_DEFAULT = 0,
102 INTERPOLATION_NEAREST,
103 INTERPOLATION_LINEAR,
104 INTERPOLATION_SPLINE,
105 INTERPOLATION_SMOOTH_SPLINE,
106 INTERPOLATION_PERIODIC_SPLINE,
107 INTERPOLATION_FOURIER,
108 // ...add new interpolation algorithms here...
109# ifndef DOXYGEN
110 __ILLEGAL_INTERPOLATION_METHOD__ // necessary for metaprogram loop and automatic Python enums
111# endif // DOXYGEN
112};
113
114 PLASK_API extern const char* interpolationMethodNames[__ILLEGAL_INTERPOLATION_METHOD__+1];
115
119 template <InterpolationMethod default_method>
120 inline InterpolationMethod getInterpolationMethod(InterpolationMethod method) {
121 return (method == INTERPOLATION_DEFAULT)? default_method : method;
122}
123
127 class InterpolationFlags {
128 unsigned char sym[3];
129 unsigned char per;
130 double lo[3], hi[3];
131
132 public:
133
134 enum class Symmetry: unsigned char {
135 NO = 0,
136 POSITIVE = 1,
137 PP = 1,
138 NP = 3,
139 PN = 5,
140 NN = 7,
141 PPP = 1,
142 NPP = 3,
143 PNP = 5,
144 NNP = 7,
145 PPN = 9,
146 NPN = 11,
147 PNN = 13,
148 NNN = 15,
149 NEGATIVE = 15,
150 };
151
152 InterpolationFlags(): sym{0,0,0}, per(0), lo{0.,0.,0.} , hi{0., 0., 0.} {}
153
154 InterpolationFlags(shared_ptr<const GeometryD<2>> geometry, Symmetry sym0, Symmetry sym1):
155 sym{geometry->isSymmetric(Geometry::DIRECTION_TRAN)? (unsigned char)sym0 : (unsigned char)0, geometry->isSymmetric(Geometry::DIRECTION_VERT)? (unsigned char)sym1 : (unsigned char)0},
156 per((unsigned char)((geometry->isPeriodic(Geometry::DIRECTION_TRAN)? 1 : 0) + (geometry->isPeriodic(Geometry::DIRECTION_VERT)? 2 : 0))),
157 lo{geometry->getChildBoundingBox().left(), geometry->getChildBoundingBox().bottom(), 0.},
158 hi{geometry->getChildBoundingBox().right(), geometry->getChildBoundingBox().top(), 0.}
159 {
160 if (geometry->isSymmetric(Geometry::DIRECTION_TRAN)) {
161 if (lo[0] < 0. && hi[0] > 0.)
162 throw Exception("interpolation: Symmetric geometry spans at both sides of transverse axis");
163 if (!sym[0]) {
164 hi[0] = max(-lo[0], hi[0]); lo[0] = -hi[0];
165 }
166 }
167 if (geometry->isSymmetric(Geometry::DIRECTION_VERT)) {
168 if (lo[1] < 0. && hi[1] > 0.)
169 throw Exception("interpolation: Symmetric geometry spans at both sides of vertical axis");
170 if (!sym[1]) {
171 hi[1] = max(-lo[1], hi[1]); lo[1] = -hi[1];
172 }
173 }
174 }
175
176 InterpolationFlags(shared_ptr<const Geometry2DCartesian> geometry): InterpolationFlags(geometry, Symmetry::POSITIVE, Symmetry::POSITIVE) {}
177 InterpolationFlags(shared_ptr<Geometry2DCartesian> geometry): InterpolationFlags(geometry, Symmetry::POSITIVE, Symmetry::POSITIVE) {}
178 InterpolationFlags(shared_ptr<const Geometry2DCylindrical> geometry): InterpolationFlags(geometry, Symmetry::POSITIVE, Symmetry::POSITIVE) {}
179 InterpolationFlags(shared_ptr<Geometry2DCylindrical> geometry): InterpolationFlags(geometry, Symmetry::POSITIVE, Symmetry::POSITIVE) {}
180
181 InterpolationFlags(shared_ptr<const GeometryD<3>> geometry, Symmetry sym0, Symmetry sym1, Symmetry sym2):
182 sym{geometry->isSymmetric(Geometry::DIRECTION_LONG)? (unsigned char)sym0 : (unsigned char)0, geometry->isSymmetric(Geometry::DIRECTION_TRAN)? (unsigned char)sym1 : (unsigned char)0, geometry->isSymmetric(Geometry::DIRECTION_VERT)? (unsigned char)sym2 : (unsigned char)0},
183 per((unsigned char)((geometry->isPeriodic(Geometry::DIRECTION_LONG)? 1 : 0) + (geometry->isPeriodic(Geometry::DIRECTION_TRAN)? 2 : 0) + (geometry->isPeriodic(Geometry::DIRECTION_VERT)? 4 : 0))),
184 lo{geometry->getChildBoundingBox().back(), geometry->getChildBoundingBox().left(), geometry->getChildBoundingBox().bottom()},
185 hi{geometry->getChildBoundingBox().front(), geometry->getChildBoundingBox().right(), geometry->getChildBoundingBox().top()}
186 {
187 if (geometry->isSymmetric(Geometry::DIRECTION_LONG)) {
188 if (lo[0] < 0. && hi[0] > 0.)
189 throw Exception("interpolation: Symmetric geometry spans at both sides of longitudinal axis");
190 if (!sym[0]) {
191 hi[0] = max(-lo[0], hi[0]); lo[0] = -hi[0];
192 }
193 }
194 if (geometry->isSymmetric(Geometry::DIRECTION_TRAN)) {
195 if (lo[1] < 0. && hi[1] > 0.)
196 throw Exception("interpolation: Symmetric geometry spans at both sides of transverse axis");
197 if (!sym[1]) {
198 hi[1] = max(-lo[1], hi[1]); lo[1] = -hi[1];
199 }
200 }
201 if (geometry->isSymmetric(Geometry::DIRECTION_VERT)) {
202 if (lo[2] < 0. && hi[2] > 0.)
203 throw Exception("interpolation: Symmetric geometry spans at both sides of vertical axis");
204 if (!sym[2]) {
205 hi[2] = max(-lo[2], hi[2]); lo[2] = -hi[2];
206 }
207 }
208 }
209
210 InterpolationFlags(shared_ptr<const Geometry3D> geometry): InterpolationFlags(geometry, Symmetry::POSITIVE, Symmetry::POSITIVE, Symmetry::POSITIVE) {}
211 InterpolationFlags(shared_ptr<Geometry3D> geometry): InterpolationFlags(geometry, Symmetry::POSITIVE, Symmetry::POSITIVE, Symmetry::POSITIVE) {}
212
213 unsigned char symmetric(int axis) const { return sym[axis]; }
214
215 bool periodic(int axis) const { return (per & (1 << axis)) != 0; }
216
217 double low(int axis) const {
218 if (sym[axis]) return min(lo[axis], -hi[axis]);
219 else return lo[axis];
220 }
221
222 double high(int axis) const {
223 if (sym[axis]) return max(lo[axis], hi[axis]);
224 else return hi[axis];
225 }
226
232 double wrap(int ax, double x) const {
233 double d = hi[ax] - lo[ax];
234 if (periodic(ax)) {
235 if (sym[ax]) {
236 x = std::fmod(abs(x), 2.*d);
237 if (x == 0.) x = 1e-12;
238 if (x > d) x = - (x - 2.*d);
239 if (hi[ax] < 0) x = - x;
240 } else {
241 x = std::fmod(x-lo[ax], d);
242 x += (x >= 0)? lo[ax] : hi[ax];
243 }
244 } else if (sym[ax]) {
245 if (x == 0.) x = 1e-12;
246 if (lo[ax] >= 0) x = abs(x);
247 else x = - abs(x);
248 }
249 return x;
250 }
251
252 template <int dim>
253 Vec<dim> wrap(Vec<dim> pos) const {
254 for (int i = 0; i != dim; ++i) {
255 pos[i] = wrap(i, pos[i]);
256 }
257 return pos;
258 }
259
260 template <int dim, typename DataT>
261 Vec<dim,DataT> reflect(int ax, Vec<dim,DataT> vec) const {
262 for (int i = 0; i != dim; ++i) if (sym[ax] & (2 << i)) vec[i] = -vec[i];
263 return vec;
264 }
265
266 template <typename DataT>
267 DataT reflect(int ax, DataT val) const {
268 if (sym[ax] & 14) return -val;
269 else return val;
270 }
271
272 template <int dim, typename DataT>
273 DataT postprocess(Vec<dim> pos, DataT data) const {
274 for (int i = 0; i != dim; ++i) {
275 if (sym[i]) {
276 if (periodic(i)) {
277 double d = hi[i] - lo[i];
278 pos[i] = std::fmod(pos[i], 2.*d);
279 if (pos[i] > d || (pos[i] < 0. && pos[i] > -d)) data = reflect(i, data);
280 } else {
281 if (lo[i] >= 0.) { if (pos[i] < 0.) data = reflect(i, data); }
282 else { if (pos[i] > 0.) data = reflect(i, data); }
283 }
284 }
285 }
286 return data;
287 }
288};
289
295 template <typename SrcMeshT, typename SrcT, typename DstT, InterpolationMethod method>
296 struct InterpolationAlgorithm
297{
298 static LazyData<DstT> interpolate(const shared_ptr<const SrcMeshT>& src_mesh, const DataVector<const SrcT>&,
299 const shared_ptr<const MeshD<SrcMeshT::DIM>>&, const InterpolationFlags&) {
300 std::string msg = "interpolate (source mesh type: ";
301 msg += typeid(*src_mesh.get()).name();
302 msg += ", interpolation method: ";
303 msg += interpolationMethodNames[method];
304 msg += ")";
305 throw NotImplemented(msg);
306 }
307};
308
312 template <typename SrcMeshT, typename SrcT, typename DstT>
313 struct InterpolationAlgorithm<SrcMeshT, SrcT, DstT, INTERPOLATION_DEFAULT>
314{
315 static LazyData<DstT> interpolate(const shared_ptr<const SrcMeshT>&, const DataVector<const SrcT>&,
316 const shared_ptr<const MeshD<SrcMeshT::DIM>>&, const InterpolationFlags& PLASK_UNUSED(flags)) {
317 throw CriticalException("interpolate(...) called for INTERPOLATION_DEFAULT method. Contact solver author to fix this issue."
318#ifndef NDEBUG
319 "\n\nINFO FOR SOLVER AUTHOR: To avoid this error use 'getInterpolationMethod<YOUR_DEFAULT_METHOD>(interpolation_method) in C++ code of the provider in your solver.\n"
320#endif
321 );
322 }
323};
324
325#ifndef DOXYGEN
326// The following structures are solely used for metaprogramming
327 template <typename SrcMeshT, typename SrcT, typename DstT, int iter>
328 struct __InterpolateMeta__
329{
330 inline static LazyData<typename std::remove_const<DstT>::type> interpolate(
331 const shared_ptr<const SrcMeshT>& src_mesh, const DataVector<const SrcT>& src_vec,
332 const shared_ptr<const MeshD<SrcMeshT::DIM>>& dst_mesh, InterpolationMethod method, const InterpolationFlags& flags) {
333 if (int(method) == iter)
334 return InterpolationAlgorithm<SrcMeshT, SrcT, typename std::remove_const<DstT>::type, (InterpolationMethod)iter>
335 ::interpolate(src_mesh, DataVector<const SrcT>(src_vec), dst_mesh, flags);
336 else
337 return __InterpolateMeta__<SrcMeshT, SrcT, DstT, iter+1>::interpolate(src_mesh, src_vec, dst_mesh, method, flags);
338 }
339};
340 template <typename SrcMeshT, typename SrcT, typename DstT>
341 struct __InterpolateMeta__<SrcMeshT, SrcT, DstT, __ILLEGAL_INTERPOLATION_METHOD__>
342{
343 inline static LazyData<typename std::remove_const<DstT>::type> interpolate(const shared_ptr<const SrcMeshT>&, const DataVector<const SrcT>&,
344 const shared_ptr<const MeshD<SrcMeshT::DIM>>&, InterpolationMethod, const InterpolationFlags&) {
345 throw CriticalException("no such interpolation method");
346 }
347};
348#endif // DOXYGEN
349
350
365 template <typename SrcMeshT, typename SrcT, typename DstT=SrcT>
366 LazyData<typename std::remove_const<DstT>::type> interpolate(
367 shared_ptr<const SrcMeshT> src_mesh,
368 DataVector<const SrcT> src_vec,
369 shared_ptr<const MeshD<SrcMeshT::DIM>> dst_mesh,
370 InterpolationMethod method=INTERPOLATION_DEFAULT,
371 const InterpolationFlags& flags=InterpolationFlags(),
372 bool verbose=true)
373{
374 if (src_mesh->size() != src_vec.size())
375 throw BadMesh("interpolate", "mesh size ({1}) and values size ({0}) do not match", src_vec.size(), src_mesh->size());
376 if (src_mesh == dst_mesh) return new LazyDataFromVectorImpl<typename std::remove_const<DstT>::type>(src_vec); // meshes are identical, so just return src_vec
377 if (verbose && method < __ILLEGAL_INTERPOLATION_METHOD__) writelog(LOG_DEBUG, "interpolate: Running {0} interpolation", interpolationMethodNames[method]);
378 return __InterpolateMeta__<SrcMeshT, SrcT, DstT, 0>::interpolate(src_mesh, src_vec, dst_mesh, method, flags);
379}
380
381 template <typename SrcMeshT, typename SrcT, typename DstT=SrcT, typename DstMeshT>
382 LazyData<typename std::remove_const<DstT>::type> interpolate(
383 shared_ptr<SrcMeshT> src_mesh,
384 DataVector<SrcT> src_vec,
385 shared_ptr<DstMeshT> dst_mesh,
386 InterpolationMethod method=INTERPOLATION_DEFAULT,
387 const InterpolationFlags& flags=InterpolationFlags(),
388 bool verbose=true)
389{
390 return interpolate(shared_ptr<const SrcMeshT>(src_mesh), DataVector<const SrcT>(src_vec),
391 shared_ptr<const MeshD<SrcMeshT::DIM>>(dst_mesh), method, flags, verbose);
392}
393
394
395/*template <typename SrcMeshT, typename SrcT, typename DstT=SrcT>
396 LazyData<typename std::remove_const<DstT>::type> interpolate(shared_ptr<SrcMeshT> src_mesh, DataVector<const SrcT> src_vec,
397 shared_ptr<MeshD<SrcMeshT::DIM>> dst_mesh,
398 InterpolationMethod method=INTERPOLATION_DEFAULT, bool verbose=true)
399{
400 return interpolate(shared_ptr<const SrcMeshT>(src_mesh), src_vec, shared_ptr<const MeshD<SrcMeshT::DIM>>(dst_mesh), method, verbose);
401}
402
403 template <typename SrcMeshT, typename SrcT, typename DstT=SrcT>
404 LazyData<typename std::remove_const<DstT>::type> interpolate(shared_ptr<const SrcMeshT> src_mesh, DataVector<const SrcT> src_vec,
405 shared_ptr<MeshD<SrcMeshT::DIM>> dst_mesh,
406 InterpolationMethod method=INTERPOLATION_DEFAULT, bool verbose=true)
407{
408 return interpolate(src_mesh, src_vec, shared_ptr<const MeshD<SrcMeshT::DIM>>(dst_mesh), method, verbose);
409}*/
410
416 template <typename DstT, typename SrcMeshType, typename SrcT = DstT>
417 struct InterpolatedLazyDataImpl: public LazyDataImpl<DstT> {
418
419 shared_ptr<const SrcMeshType> src_mesh;
420 shared_ptr<const MeshD<SrcMeshType::DIM>> dst_mesh;
421 DataVector<const SrcT> src_vec;
422 InterpolationFlags flags;
423
424 InterpolatedLazyDataImpl(const shared_ptr<const SrcMeshType>& src_mesh, const DataVector<const SrcT>& src_vec,
425 const shared_ptr<const MeshD<SrcMeshType::DIM>>& dst_mesh, const InterpolationFlags& flags)
426 : src_mesh(src_mesh), dst_mesh(dst_mesh), src_vec(src_vec), flags(flags) {}
427
428 std::size_t size() const override { return dst_mesh->size(); }
429
430};
431
437 template <typename DstT, typename SrcMeshType, typename SrcT = DstT>
438 struct LinearInterpolatedLazyDataImpl: public InterpolatedLazyDataImpl<DstT, SrcMeshType, SrcT> {
439
440 LinearInterpolatedLazyDataImpl(shared_ptr<const SrcMeshType> src_mesh, const DataVector<const SrcT>& src_vec,
441 shared_ptr<const MeshD<SrcMeshType::DIM>> dst_mesh, const InterpolationFlags& flags):
442 InterpolatedLazyDataImpl<DstT, SrcMeshType>(src_mesh, src_vec, dst_mesh, flags) {}
443
444 DstT at(std::size_t index) const override {
445 return this->src_mesh->interpolateLinear(this->src_vec, this->dst_mesh->at(index), this->flags);
446 }
447
448};
449
455 template <typename DstT, typename SrcMeshType, typename SrcT=DstT>
456 struct NearestNeighborInterpolatedLazyDataImpl: public InterpolatedLazyDataImpl<DstT, SrcMeshType, SrcT> {
457
458 NearestNeighborInterpolatedLazyDataImpl(shared_ptr<const SrcMeshType> src_mesh, const DataVector<const SrcT>& src_vec,
459 shared_ptr<const MeshD<SrcMeshType::DIM>> dst_mesh, const InterpolationFlags& flags):
460 InterpolatedLazyDataImpl<DstT, SrcMeshType, DstT>(src_mesh, src_vec, dst_mesh, flags) {}
461
462 SrcT at(std::size_t index) const override {
463 return this->src_mesh->interpolateNearestNeighbor(this->src_vec, this->dst_mesh->at(index), this->flags);
464 }
465
466};
467
468
469
470} // namespace plask
471
472
473 namespace boost {
474
475 template <>
476 inline plask::InterpolationMethod lexical_cast<plask::InterpolationMethod, std::string>(const std::string& arg) {
477 std::string val = arg; to_upper(val);
478 for (unsigned i = plask::INTERPOLATION_DEFAULT+1; i != plask::__ILLEGAL_INTERPOLATION_METHOD__; ++i) {
479 if (val == plask::interpolationMethodNames[i]) return (plask::InterpolationMethod)i;
480 }
481 throw bad_lexical_cast(typeid(std::string), typeid(plask::InterpolationMethod));
482 }
483
484}
485
486
487
488#endif //PLASK__INTERPOLATION_H