PLaSK library
Loading...
Searching...
No Matches
rectangular_spline.cpp
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#include "rectangular_spline.hpp"
15#include "hyman.hpp"
16#include "../exceptions.hpp"
17
18namespace plask {
19
20
21template <typename DstT, typename SrcT>
22SplineRect2DLazyDataImpl<DstT, SrcT>::SplineRect2DLazyDataImpl(const shared_ptr<const RectangularMesh2D>& src_mesh,
23 const DataVector<const SrcT>& src_vec,
24 const shared_ptr<const MeshD<2>>& dst_mesh,
25 const InterpolationFlags& flags):
26 InterpolatedLazyDataImpl<DstT, RectangularMesh2D, const SrcT>(src_mesh, src_vec, dst_mesh, flags),
27 diff0(src_mesh->size()), diff1(src_mesh->size()) {}
28
29template <typename DstT, typename SrcT>
30DstT SplineRect2DLazyDataImpl<DstT, SrcT>::at(std::size_t index) const
31{
32 Vec<2> p = this->flags.wrap(this->dst_mesh->at(index));
33
34 size_t i0_lo, i0_hi;
35 double left, right;
36 bool invert_left, invert_right;
37 prepareInterpolationForAxis(*this->src_mesh->axis[0], this->flags, p.c0, 0, i0_lo, i0_hi, left, right, invert_left, invert_right);
38
39 size_t i1_lo, i1_hi;
40 double bottom, top;
41 bool invert_bottom, invert_top;
42 prepareInterpolationForAxis(*this->src_mesh->axis[1], this->flags, p.c1, 1, i1_lo, i1_hi, bottom, top, invert_bottom, invert_top);
43
44 double d0 = right - left,
45 d1 = top - bottom;
46 double x0 = (p.c0 - left) / d0,
47 x1 = (p.c1 - bottom) / d1;
48
49 // Hermite 3rd order spline polynomials (in Horner form)
50 double hl = ( 2.*x0 - 3.) * x0*x0 + 1.,
51 hr = (-2.*x0 + 3.) * x0*x0,
52 gl = ((x0 - 2.) * x0 + 1.) * x0 * d0,
53 gr = (x0 - 1.) * x0 * x0 * d0,
54 hb = ( 2.*x1 - 3.) * x1*x1 + 1.,
55 ht = (-2.*x1 + 3.) * x1*x1,
56 gb = ((x1 - 2.) * x1 + 1.) * x1 * d1,
57 gt = (x1 - 1.) * x1 * x1 * d1;
58
59 std::size_t ilb = this->src_mesh->index(i0_lo, i1_lo),
60 ilt = this->src_mesh->index(i0_lo, i1_hi),
61 irb = this->src_mesh->index(i0_hi, i1_lo),
62 irt = this->src_mesh->index(i0_hi, i1_hi);
63
64 SrcT diff0lb = diff0[ilb],
65 diff0lt = diff0[ilt],
66 diff0rb = diff0[irb],
67 diff0rt = diff0[irt],
68 diff1lb = diff1[ilb],
69 diff1lt = diff1[ilt],
70 diff1rb = diff1[irb],
71 diff1rt = diff1[irt];
72
73 if (invert_left) { diff0lb = -this->flags.reflect(0, diff0lb); diff0lt = -this->flags.reflect(0, diff0lt); };
74 if (invert_right) { diff0rb = -this->flags.reflect(0, diff0rb); diff0rt = -this->flags.reflect(0, diff0rt); };
75 if (invert_top) { diff1lt = -this->flags.reflect(1, diff1lt); diff1rt = -this->flags.reflect(1, diff1rt); };
76 if (invert_bottom) { diff1lb = -this->flags.reflect(1, diff1lb); diff1rb = -this->flags.reflect(1, diff1rb); };
77
78 SrcT data_lb = this->src_vec[ilb],
79 data_lt = this->src_vec[ilt],
80 data_rb = this->src_vec[irb],
81 data_rt = this->src_vec[irt],
82 diff_l = gb * diff1lb + gt * diff1lt,
83 diff_r = gb * diff1rb + gt * diff1rt,
84 diff_b = gl * diff0lb + gr * diff0rb,
85 diff_t = gl * diff0lt + gr * diff0rt;
86
87 if (invert_left) { data_lb = this->flags.reflect(0, data_lb); data_lt = this->flags.reflect(0, data_lt); diff_l = this->flags.reflect(0, diff_l); }
88 if (invert_right) { data_rb = this->flags.reflect(0, data_rb); data_rt = this->flags.reflect(0, data_rt); diff_r = this->flags.reflect(0, diff_r); }
89 if (invert_top) { data_lt = this->flags.reflect(1, data_lt); data_rt = this->flags.reflect(1, data_rt); diff_t = this->flags.reflect(1, diff_t); }
90 if (invert_bottom) { data_lb = this->flags.reflect(1, data_lb); data_rb = this->flags.reflect(1, data_rb); diff_b = this->flags.reflect(1, diff_b); }
91
92 return this->flags.postprocess(this->dst_mesh->at(index),
93 hl * (hb * data_lb + ht * data_lt) + hr * (hb * data_rb + ht * data_rt) +
94 hb * diff_b + ht * diff_t + hl * diff_l + hr * diff_r
95 );
96}
97
98
99template <typename DstT, typename SrcT>
100SplineRect3DLazyDataImpl<DstT, SrcT>::SplineRect3DLazyDataImpl(const shared_ptr<const RectilinearMesh3D>& src_mesh,
101 const DataVector<const SrcT>& src_vec,
102 const shared_ptr<const MeshD<3>>& dst_mesh,
103 const InterpolationFlags& flags):
104 InterpolatedLazyDataImpl<DstT, RectilinearMesh3D, const SrcT>(src_mesh, src_vec, dst_mesh, flags),
105 diff0(src_mesh->size()), diff1(src_mesh->size()), diff2(src_mesh->size()) {}
106
107template <typename DstT, typename SrcT>
108DstT SplineRect3DLazyDataImpl<DstT, SrcT>::at(std::size_t index) const
109{
110 Vec<3> p = this->flags.wrap(this->dst_mesh->at(index));
111
112 size_t i0_lo, i0_hi;
113 double back, front;
114 bool invert_back, invert_front;
115 prepareInterpolationForAxis(*this->src_mesh->axis[0], this->flags, p.c0, 0, i0_lo, i0_hi, back, front, invert_back, invert_front);
116
117 size_t i1_lo, i1_hi;
118 double left, right;
119 bool invert_left, invert_right;
120 prepareInterpolationForAxis(*this->src_mesh->axis[1], this->flags, p.c1, 1, i1_lo, i1_hi, left, right, invert_left, invert_right);
121
122 size_t i2_lo, i2_hi;
123 double bottom, top;
124 bool invert_bottom, invert_top;
125 prepareInterpolationForAxis(*this->src_mesh->axis[2], this->flags, p.c2, 2, i2_lo, i2_hi, bottom, top, invert_bottom, invert_top);
126
127 std::size_t illl = this->src_mesh->index(i0_lo, i1_lo, i2_lo),
128 illh = this->src_mesh->index(i0_lo, i1_lo, i2_hi),
129 ilhl = this->src_mesh->index(i0_lo, i1_hi, i2_lo),
130 ilhh = this->src_mesh->index(i0_lo, i1_hi, i2_hi),
131 ihll = this->src_mesh->index(i0_hi, i1_lo, i2_lo),
132 ihlh = this->src_mesh->index(i0_hi, i1_lo, i2_hi),
133 ihhl = this->src_mesh->index(i0_hi, i1_hi, i2_lo),
134 ihhh = this->src_mesh->index(i0_hi, i1_hi, i2_hi);
135
136 double d0 = front - back,
137 d1 = right - left,
138 d2 = top - bottom;
139 double x0 = (p.c0 - back) / d0,
140 x1 = (p.c1 - left) / d1,
141 x2 = (p.c2 - bottom) / d2;
142
143 // Hermite 3rd order spline polynomials (in Horner form)
144 double h0l = ( 2.*x0 - 3.) * x0*x0 + 1.,
145 h0h = (-2.*x0 + 3.) * x0*x0,
146 g0l = ((x0 - 2.) * x0 + 1.) * x0 * d0,
147 g0h = (x0 - 1.) * x0 * x0 * d0,
148 h1l = ( 2.*x1 - 3.) * x1*x1 + 1.,
149 h1h = (-2.*x1 + 3.) * x1*x1,
150 g1l = ((x1 - 2.) * x1 + 1.) * x1 * d1,
151 g1h = (x1 - 1.) * x1 * x1 * d1,
152 h2l = ( 2.*x2 - 3.) * x2*x2 + 1.,
153 h2h = (-2.*x2 + 3.) * x2*x2,
154 g2l = ((x2 - 2.) * x2 + 1.) * x2 * d2,
155 g2h = (x2 - 1.) * x2 * x2 * d2;
156
157 SrcT diff0lll = diff0[illl],
158 diff0llh = diff0[illh],
159 diff0lhl = diff0[ilhl],
160 diff0lhh = diff0[ilhh],
161 diff0hll = diff0[ihll],
162 diff0hlh = diff0[ihlh],
163 diff0hhl = diff0[ihhl],
164 diff0hhh = diff0[ihhh],
165 diff1lll = diff1[illl],
166 diff1llh = diff1[illh],
167 diff1lhl = diff1[ilhl],
168 diff1lhh = diff1[ilhh],
169 diff1hll = diff1[ihll],
170 diff1hlh = diff1[ihlh],
171 diff1hhl = diff1[ihhl],
172 diff1hhh = diff1[ihhh],
173 diff2lll = diff2[illl],
174 diff2llh = diff2[illh],
175 diff2lhl = diff2[ilhl],
176 diff2lhh = diff2[ilhh],
177 diff2hll = diff2[ihll],
178 diff2hlh = diff2[ihlh],
179 diff2hhl = diff2[ihhl],
180 diff2hhh = diff2[ihhh];
181
182 if (invert_back) { diff0lll = -this->flags.reflect(0, diff0lll); diff0llh = -this->flags.reflect(0, diff0llh);
183 diff0lhl = -this->flags.reflect(0, diff0lhl); diff0lhh = -this->flags.reflect(0, diff0lhh); };
184 if (invert_front) { diff0hll = -this->flags.reflect(0, diff0hll); diff0hlh = -this->flags.reflect(0, diff0hlh);
185 diff0hhl = -this->flags.reflect(0, diff0hhl); diff0hhh = -this->flags.reflect(0, diff0hhh); };
186 if (invert_left) { diff1lll = -this->flags.reflect(1, diff1lll); diff0llh = -this->flags.reflect(1, diff1llh);
187 diff1hll = -this->flags.reflect(1, diff1hll); diff0hlh = -this->flags.reflect(1, diff1hlh); };
188 if (invert_right) { diff1lhl = -this->flags.reflect(1, diff1lhl); diff0lhh = -this->flags.reflect(1, diff1lhh);
189 diff1hhl = -this->flags.reflect(1, diff1hhl); diff0hhh = -this->flags.reflect(1, diff1hhh); };
190 if (invert_top) { diff2lll = -this->flags.reflect(2, diff2lll); diff0lhl = -this->flags.reflect(2, diff2lhl);
191 diff2hll = -this->flags.reflect(2, diff2hll); diff0hhl = -this->flags.reflect(2, diff2hhl); };
192 if (invert_bottom) { diff2llh = -this->flags.reflect(2, diff2llh); diff0lhh = -this->flags.reflect(2, diff2lhh);
193 diff2hlh = -this->flags.reflect(2, diff2hlh); diff0hhh = -this->flags.reflect(2, diff2hhh); };
194
195
196 SrcT data_lll = this->src_vec[illl],
197 data_llh = this->src_vec[illh],
198 data_lhl = this->src_vec[ilhl],
199 data_lhh = this->src_vec[ilhh],
200 data_hll = this->src_vec[ihll],
201 data_hlh = this->src_vec[ihlh],
202 data_hhl = this->src_vec[ihhl],
203 data_hhh = this->src_vec[ihhh],
204 D_ll = g0l * diff0lll + g0h * diff0hll, Dl_l = g1l * diff1lll + g1h * diff1lhl, Dll_ = g2l * diff2lll + g2h * diff2llh,
205 D_lh = g0l * diff0llh + g0h * diff0hlh, Dl_h = g1l * diff1llh + g1h * diff1lhh, Dlh_ = g2l * diff2lhl + g2h * diff2lhh,
206 D_hl = g0l * diff0lhl + g0h * diff0hhl, Dh_l = g1l * diff1hll + g1h * diff1hhl, Dhl_ = g2l * diff2hll + g2h * diff2hlh,
207 D_hh = g0l * diff0lhh + g0h * diff0hhh, Dh_h = g1l * diff1hlh + g1h * diff1hhh, Dhh_ = g2l * diff2hhl + g2h * diff2hhh;
208
209
210 if (invert_back) { data_lll = this->flags.reflect(0, data_lll); data_llh = this->flags.reflect(0, data_llh);
211 data_lhl = this->flags.reflect(0, data_lhl); data_lhh = this->flags.reflect(0, data_lhh);
212 Dl_l = this->flags.reflect(0, Dl_l); Dl_h = this->flags.reflect(0, Dl_h);
213 Dll_ = this->flags.reflect(0, Dll_); Dlh_ = this->flags.reflect(0, Dlh_); }
214 if (invert_front) { data_hll = this->flags.reflect(0, data_hll); data_llh = this->flags.reflect(0, data_hlh);
215 data_lhl = this->flags.reflect(0, data_hhl); data_lhh = this->flags.reflect(0, data_hhh);
216 Dh_l = this->flags.reflect(0, Dh_l); Dh_h = this->flags.reflect(0, Dh_h);
217 Dhl_ = this->flags.reflect(0, Dhl_); Dhh_ = this->flags.reflect(0, Dhh_); }
218 if (invert_left) { data_lll = this->flags.reflect(1, data_lll); data_llh = this->flags.reflect(1, data_llh);
219 data_hll = this->flags.reflect(1, data_hll); data_hlh = this->flags.reflect(1, data_hlh);
220 Dll_ = this->flags.reflect(1, Dll_); Dhl_ = this->flags.reflect(1, Dhl_);
221 D_ll = this->flags.reflect(1, D_ll); D_lh = this->flags.reflect(1, D_lh); }
222 if (invert_right) { data_lhl = this->flags.reflect(1, data_lhl); data_llh = this->flags.reflect(1, data_lhh);
223 data_hll = this->flags.reflect(1, data_hhl); data_hlh = this->flags.reflect(1, data_hhh);
224 Dlh_ = this->flags.reflect(1, Dlh_); Dhh_ = this->flags.reflect(1, Dhh_);
225 D_hl = this->flags.reflect(1, D_hl); D_hh = this->flags.reflect(1, D_hh); }
226 if (invert_bottom) { data_lll = this->flags.reflect(2, data_lll); data_lhl = this->flags.reflect(2, data_lhl);
227 data_hll = this->flags.reflect(2, data_hll); data_hhl = this->flags.reflect(2, data_hhl);
228 D_ll = this->flags.reflect(2, D_ll); D_hl = this->flags.reflect(2, D_hl);
229 Dl_l = this->flags.reflect(2, Dl_l); Dh_l = this->flags.reflect(2, Dh_l); }
230 if (invert_top) { data_llh = this->flags.reflect(2, data_llh); data_lhl = this->flags.reflect(2, data_lhh);
231 data_hll = this->flags.reflect(2, data_hlh); data_hhl = this->flags.reflect(2, data_hhh);
232 D_lh = this->flags.reflect(2, D_lh); D_hh = this->flags.reflect(2, D_hh);
233 Dl_h = this->flags.reflect(2, Dl_h); Dh_h = this->flags.reflect(2, Dh_h); }
234
235 return this->flags.postprocess(this->dst_mesh->at(index),
236 h0l * h1l * h2l * data_lll +
237 h0l * h1l * h2h * data_llh +
238 h0l * h1h * h2l * data_lhl +
239 h0l * h1h * h2h * data_lhh +
240 h0h * h1l * h2l * data_hll +
241 h0h * h1l * h2h * data_hlh +
242 h0h * h1h * h2l * data_hhl +
243 h0h * h1h * h2h * data_hhh +
244 h1l * h2l * D_ll + h0l * h2l * Dl_l + h0l * h1l * Dll_ +
245 h1l * h2h * D_lh + h0l * h2h * Dl_h + h0l * h1h * Dlh_ +
246 h1h * h2l * D_hl + h0h * h2l * Dh_l + h0h * h1l * Dhl_ +
247 h1h * h2h * D_hh + h0h * h2h * Dh_h + h0h * h1h * Dhh_
248 );
249}
250
251
252
253namespace hyman {
254 template <typename DataT>
255 static void computeDiffs(DataT* diffs, int ax, const shared_ptr<MeshAxis>& axis,
256 const DataT* data, std::size_t stride, const InterpolationFlags& flags)
257 {
258 const std::size_t n1 = axis->size() - 1;
259
260 for (std::size_t i = 1; i != n1; ++i) {
261 const std::size_t idx = stride * i;
262 const double da = axis->at(i) - axis->at(i-1),
263 db = axis->at(i+1) - axis->at(i);
264 const DataT sa = (data[idx] - data[idx-stride]) / da,
265 sb = (data[idx+stride] - data[idx]) / db;
266 // Use parabolic estimation of the derivative
267 diffs[idx] = (da * sb + db * sa) / (da + db);
268 // Hyman filter
269 Hyman<DataT>::filter(diffs[idx], sa, sb);
270 }
271
272 const size_t in0 = stride * n1;
273 double da0, db0, dan, dbn;
274 DataT sa0, sb0, san, sbn;
275
276 if (flags.symmetric(ax)) {
277 da0 = axis->at(0);
278 db0 = axis->at(1) - axis->at(0);
279 sb0 = (data[1] - data[0]) / db0;
280 if (da0 < 0. && flags.periodic(ax)) {
281 da0 += flags.high(ax) - flags.low(ax);
282 }
283 if (da0 == 0.)
284 sa0 = (data[1] - flags.reflect(ax, data[1])) / (2.*db0);
285 else if (da0 > 0.)
286 sa0 = (data[0] - flags.reflect(ax, data[0])) / (2.*da0);
287 else {
288 da0 = db0 = 0.5;
289 sa0 = sb0 = Zero<DataT>();
290 }
291 dan = axis->at(n1) - axis->at(n1-1);
292 san = (data[in0] - data[in0-stride]) / dan;
293 dbn = - axis->at(n1);
294 if (dbn < 0. && flags.periodic(ax)) {
295 dbn += flags.high(ax) - flags.low(ax);
296 }
297 if (dbn == 0.)
298 sbn = (data[in0-stride] - flags.reflect(ax, data[in0-stride])) / (2.*dan);
299 else if (dbn > 0.)
300 sbn = (data[in0] - flags.reflect(ax, data[in0])) / (2.*dbn);
301 else {
302 dan = dbn = 0.5;
303 san = sbn = Zero<DataT>();
304 }
305 } else if (flags.periodic(ax)) {
306 da0 = axis->at(0) - axis->at(n1) + flags.high(ax) - flags.low(ax);
307 db0 = axis->at(1) - axis->at(0);
308 dan = axis->at(n1) - axis->at(n1-1);
309 dbn = da0;
310 sb0 = (data[1] - data[0]) / db0;
311 san = (data[in0] - data[in0-stride]) / dan;
312 if (da0 == 0.) { da0 = dan; sa0 = san; }
313 else sa0 = (data[0] - data[in0]) / da0;
314 if (dbn == 0.) { dbn = db0; sbn = sb0; }
315 else sbn = (data[0] - data[in0]) / dbn;
316 } else {
317 da0 = db0 = dan = dbn = 0.5;
318 sa0 = sb0 = san = sbn = Zero<DataT>();
319 }
320
321 // Use parabolic estimation of the derivative
322 diffs[0] = (da0 * sb0 + db0 * sa0) / (da0 + db0);
323 diffs[in0] = (dan * sbn + dbn * san) / (dan + dbn);
324 // Hyman filter
325 Hyman<DataT>::filter(diffs[0], sa0, sb0);
326 Hyman<DataT>::filter(diffs[in0], san, sbn);
327 }
328}
329
330
331template <typename DstT, typename SrcT>
332HymanSplineRect2DLazyDataImpl<DstT, SrcT>::HymanSplineRect2DLazyDataImpl(const shared_ptr<const RectangularMesh2D>& src_mesh,
333 const DataVector<const SrcT>& src_vec,
334 const shared_ptr<const MeshD<2>>& dst_mesh,
335 const InterpolationFlags& flags):
336 SplineRect2DLazyDataImpl<DstT, SrcT>(src_mesh, src_vec, dst_mesh, flags)
337{
338 const int n0 = int(src_mesh->axis[0]->size()), n1 = int(src_mesh->axis[1]->size());
339
340 if (n0 == 0 || n1 == 0)
341 throw BadMesh("interpolate", "source mesh empty");
342
343 size_t stride0 = src_mesh->index(1, 0),
344 stride1 = src_mesh->index(0, 1);
345
346 if (n0 > 1)
347 for (size_t i1 = 0, i = 0; i1 < src_mesh->axis[1]->size(); ++i1, i += stride1)
348 hyman::computeDiffs<SrcT>(this->diff0.data()+i, 0, src_mesh->axis[0], src_vec.data()+i, stride0, flags);
349 else
350 std::fill(this->diff0.begin(), this->diff0.end(), Zero<SrcT>());
351 if (n1 > 1)
352 for (size_t i0 = 0, i = 0; i0 < src_mesh->axis[0]->size(); ++i0, i += stride0)
353 hyman::computeDiffs<SrcT>(this->diff1.data()+i, 1, src_mesh->axis[1], src_vec.data()+i, stride1, flags);
354 else
355 std::fill(this->diff1.begin(), this->diff1.end(), Zero<SrcT>());
356}
357
358
359
360template <typename DstT, typename SrcT>
361HymanSplineRect3DLazyDataImpl<DstT, SrcT>::HymanSplineRect3DLazyDataImpl(const shared_ptr<const RectilinearMesh3D>& src_mesh,
362 const DataVector<const SrcT>& src_vec,
363 const shared_ptr<const MeshD<3>>& dst_mesh,
364 const InterpolationFlags& flags):
365 SplineRect3DLazyDataImpl<DstT, SrcT>(src_mesh, src_vec, dst_mesh, flags)
366{
367 const int n0 = int(src_mesh->axis[0]->size()), n1 = int(src_mesh->axis[1]->size()), n2 = int(src_mesh->axis[2]->size());
368
369 if (n0 == 0 || n1 == 0 || n2 == 0)
370 throw BadMesh("interpolate", "source mesh empty");
371
372 if (n0 > 1) {
373 size_t stride0 = src_mesh->index(1, 0, 0);
374 for (size_t i2 = 0; i2 < src_mesh->axis[2]->size(); ++i2) {
375 for (size_t i1 = 0; i1 < src_mesh->axis[1]->size(); ++i1) {
376 size_t offset = src_mesh->index(0, i1, i2);
377 hyman::computeDiffs<SrcT>(this->diff0.data()+offset, 0, src_mesh->axis[0], src_vec.data()+offset, stride0, flags);
378 }
379 }
380 } else
381 std::fill(this->diff0.begin(), this->diff0.end(), Zero<SrcT>());
382
383 if (n1 > 1) {
384 size_t stride1 = src_mesh->index(0, 1, 0);
385 for (size_t i2 = 0; i2 < src_mesh->axis[2]->size(); ++i2) {
386 for (size_t i0 = 0; i0 < src_mesh->axis[0]->size(); ++i0) {
387 size_t offset = src_mesh->index(i0, 0, i2);
388 hyman::computeDiffs<SrcT>(this->diff1.data()+offset, 1, src_mesh->axis[1], src_vec.data()+offset, stride1, flags);
389 }
390 }
391 } else
392 std::fill(this->diff1.begin(), this->diff1.end(), Zero<SrcT>());
393
394 if (n2 > 1) {
395 size_t stride2 = src_mesh->index(0, 0, 1);
396 for (size_t i1 = 0; i1 < src_mesh->axis[1]->size(); ++i1) {
397 for (size_t i0 = 0; i0 < src_mesh->axis[0]->size(); ++i0) {
398 size_t offset = src_mesh->index(i0, i1, 0);
399 hyman::computeDiffs<SrcT>(this->diff2.data()+offset, 2, src_mesh->axis[2], src_vec.data()+offset, stride2, flags);
400 }
401 }
402 } else
403 std::fill(this->diff2.begin(), this->diff2.end(), Zero<SrcT>());
404
405}
406
407
408template struct PLASK_API HymanSplineRect2DLazyDataImpl<double, double>;
409template struct PLASK_API HymanSplineRect2DLazyDataImpl<dcomplex, dcomplex>;
410
411template struct PLASK_API HymanSplineRect2DLazyDataImpl<Vec<2,double>, Vec<2,double>>;
412template struct PLASK_API HymanSplineRect2DLazyDataImpl<Vec<2,dcomplex>, Vec<2,dcomplex>>;
413
414template struct PLASK_API HymanSplineRect2DLazyDataImpl<Vec<3,double>, Vec<3,double>>;
415template struct PLASK_API HymanSplineRect2DLazyDataImpl<Vec<3,dcomplex>, Vec<3,dcomplex>>;
416
417template struct PLASK_API HymanSplineRect2DLazyDataImpl<Tensor2<double>, Tensor2<double>>;
418template struct PLASK_API HymanSplineRect2DLazyDataImpl<Tensor2<dcomplex>, Tensor2<dcomplex>>;
419
420template struct PLASK_API HymanSplineRect2DLazyDataImpl<Tensor3<double>, Tensor3<double>>;
421template struct PLASK_API HymanSplineRect2DLazyDataImpl<Tensor3<dcomplex>, Tensor3<dcomplex>>;
422
423template struct PLASK_API HymanSplineRect3DLazyDataImpl<double, double>;
424template struct PLASK_API HymanSplineRect3DLazyDataImpl<dcomplex, dcomplex>;
425
426template struct PLASK_API HymanSplineRect3DLazyDataImpl<Vec<2,double>, Vec<2,double>>;
427template struct PLASK_API HymanSplineRect3DLazyDataImpl<Vec<2,dcomplex>, Vec<2,dcomplex>>;
428
429template struct PLASK_API HymanSplineRect3DLazyDataImpl<Vec<3,double>, Vec<3,double>>;
430template struct PLASK_API HymanSplineRect3DLazyDataImpl<Vec<3,dcomplex>, Vec<3,dcomplex>>;
431
432template struct PLASK_API HymanSplineRect3DLazyDataImpl<Tensor2<double>, Tensor2<double>>;
433template struct PLASK_API HymanSplineRect3DLazyDataImpl<Tensor2<dcomplex>, Tensor2<dcomplex>>;
434
435template struct PLASK_API HymanSplineRect3DLazyDataImpl<Tensor3<double>, Tensor3<double>>;
436template struct PLASK_API HymanSplineRect3DLazyDataImpl<Tensor3<dcomplex>, Tensor3<dcomplex>>;
437
438
439
440namespace spline {
441 template <typename DataT>
442 static void computeDiffs(DataT* data, size_t stride, size_t stride1, size_t size1, size_t stride2, size_t size2,
443 int ax, const shared_ptr<MeshAxis>& axis, const InterpolationFlags& flags)
444 {
445 const size_t n0 = axis->size();
446 const size_t n1 = n0 - 1;
447
448 std::unique_ptr<double[]>
449 dl(new double[n1]),
450 dd(new double[n1 + 1]),
451 du(new double[n1]);
452
453 std::unique_ptr<DataT[]> lastdata(new DataT[size1*size2]);
454 {
455 double left = (axis->at(0) >= 0.)? 0. : flags.low(ax);
456 const double da = 2. * (axis->at(0) - left), db = axis->at(1) - axis->at(0);
457 const double dab = da/db, dba = db/da;
458 for (size_t c1 = 0, s1 = 0; c1 != size1; ++c1, s1 += stride1) {
459 size_t cc = c1 * size2;
460 for (size_t c2 = 0, s2 = 0; c2 != size2; ++c2, s2 += stride2) {
461 lastdata[cc+c2] = data[s1+s2];
462 if (da == 0.) {
463 data[s1+s2] = 0.5 * (data[s1+s2+stride] - flags.reflect(ax, data[s1+s2+stride])) / db;
464 } else {
465 data[s1+s2] = 3. * (dab * data[s1+s2+stride] + (dba-dab) * lastdata[cc+c2] - dba * flags.reflect(ax, lastdata[cc+c2]));
466 }
467 }
468 }
469 }
470
471 for (size_t i = 1, si = stride; i != n1; ++i, si += stride) {
472 const double da = axis->at(i) - axis->at(i-1), // d[i-1]
473 db = axis->at(i+1) - axis->at(i); // d[i]
474 const double dab = da/db, dba = db/da;
475 dl[i-1] = db;
476 dd[i] = 2. * (da + db);
477 du[i] = da;
478 for (size_t c1 = 0, s1 = 0; c1 != size1; ++c1, s1 += stride1) {
479 size_t cc = c1 * size2;
480 for (size_t c2 = 0, s2 = 0; c2 != size2; ++c2, s2 += stride2) {
481 DataT current = data[s1+s2+si];
482 data[s1+s2+si] = 3. * (dab * data[s1+s2+si+stride] + (dba-dab) * current - dba * lastdata[cc+c2]);
483 lastdata[cc+c2] = current;
484 }
485 }
486 }
487
488 if (!flags.periodic(ax) || flags.symmetric(ax)) {
489
490 if (flags.symmetric(ax)) {
491 if (axis->at(0) == 0. || (flags.periodic(ax) && is_zero(axis->at(0) - flags.low(ax)))) {
492 dd[0] = 1.;
493 du[0] = 0.;
494 } else if (axis->at(0) > 0. || flags.periodic(ax)) {
495 const double left = (axis->at(0) >= 0.)? 0. : flags.low(ax);
496 // dd[0] = 3.*axis->at(0) + axis->at(1) - 4.*left;
497 dd[0] = 2.*axis->at(0) + 2.*axis->at(1) - 4.*left; // hack, but should work for asymmetric as well
498 du[0] = 2. * (axis->at(0) - left);
499 } else {
500 dd[0] = 1.;
501 du[0] = 0.;
502 for (size_t c1 = 0, s1 = 0; c1 != size1; ++c1, s1 += stride1)
503 for (size_t c2 = 0, s2 = 0; c2 != size2; ++c2, s2 += stride2)
504 data[s1+s2] = Zero<DataT>();
505 }
506 if (axis->at(n1) == 0. || (flags.periodic(ax) && is_zero(flags.high(ax) - axis->at(n1)))) {
507 dd[n1] = 1.;
508 dl[n1-1] = 0.;
509 size_t ns = n1 * stride;
510 double ih = 0.5 / (axis->at(n1) - axis->at(n1-1));
511 for (size_t c1 = 0, s1 = 0; c1 != size1; ++c1, s1 += stride1) {
512 size_t cc = c1 * size2;
513 for (size_t c2 = 0, s2 = 0; c2 != size2; ++c2, s2 += stride2) {
514 data[s1+s2+ns] = (flags.reflect(ax, lastdata[cc+c2]) - lastdata[cc+c2]) * ih;
515 }
516 }
517 } else if (axis->at(n1) < 0. || flags.periodic(ax)) {
518 const double right = (axis->at(n1) <= 0.)? 0. : flags.high(ax);
519 const double da = axis->at(n1) - axis->at(n1-1), db = 2. * (right - axis->at(n1));
520 const double dab = da/db, dba = db/da;
521 dl[n1-1] = db;
522 // dd[n1] = 4.*right - 3.*axis->at(n1) - axis->at(n1-1);
523 dd[n1] = 4.*right - 2.*axis->at(n1) - 2.*axis->at(n1-1); // hack
524 size_t ns = n1 * stride;
525 for (size_t c1 = 0, s1 = 0; c1 != size1; ++c1, s1 += stride1) {
526 size_t cc = c1 * size2;
527 for (size_t c2 = 0, s2 = 0; c2 != size2; ++c2, s2 += stride2)
528 data[s1+s2+ns] = 3. * (dab * flags.reflect(ax, data[s1+s2+ns]) + (dba-dab) * data[s1+s2+ns] - dba * lastdata[cc+c2]);
529 }
530 } else {
531 dl[n1-1] = 0.;
532 dd[n1] = 1.;
533 size_t ns = n1 * stride;
534 for (size_t c1 = 0, s1 = 0; c1 != size1; ++c1, s1 += stride1) {
535 for (size_t c2 = 0, s2 = 0; c2 != size2; ++c2, s2 += stride2)
536 data[s1+s2+ns] = Zero<DataT>();
537 }
538 }
539 } else {
540 du[0] = dl[n1-1] = 0.;
541 dd[0] = dd[n1] = 1.;
542 size_t ns = n1 * stride;
543 for (size_t c1 = 0, s1 = 0; c1 != size1; ++c1, s1 += stride1) {
544 for (size_t c2 = 0, s2 = 0; c2 != size2; ++c2, s2 += stride2)
545 data[s1+s2] = data[s1+s2+ns] = Zero<DataT>();
546 }
547 }
548
549 // Thomas algorithm
550 double id0 = 1. / dd[0];
551 du[0] *= id0;
552 for (size_t c1 = 0, s1 = 0; c1 != size1; ++c1, s1 += stride1) {
553 for (size_t c2 = 0, s2 = 0; c2 != size2; ++c2, s2 += stride2)
554 data[s1+s2] *= id0;
555 }
556
557 /* loop from 1 to X - 1 inclusive, performing the forward sweep */
558 for (size_t i = 1, si = stride; i < n0; i++, si += stride) {
559 const double m = 1. / (dd[i] - dl[i-1] * du[i-1]);
560 du[i] *= m;
561 for (size_t c1 = 0, s1 = 0; c1 != size1; ++c1, s1 += stride1) {
562 for (size_t c2 = 0, s2 = 0; c2 != size2; ++c2, s2 += stride2)
563 data[s1+s2+si] = (data[s1+s2+si] - dl[i-1] * data[s1+s2+si-stride]) * m;
564 }
565 }
566
567 /* loop from X - 2 to 0 inclusive (safely testing loop condition for an unsignedeger), to perform the back substitution */
568 for (size_t i = n1, si = n1*stride; i-- > 0; si -= stride) {
569 for (size_t c1 = 0, s1 = 0; c1 != size1; ++c1, s1 += stride1) {
570 for (size_t c2 = 0, s2 = 0; c2 != size2; ++c2, s2 += stride2)
571 data[s1+s2+si-stride] -= du[i] * data[s1+s2+si];
572 }
573 }
574 } else {
575 throw NotImplemented("smooth spline for periodic, non-symmetric geometry");
576 }
577 }
578
579}
580
581
582template <typename DstT, typename SrcT>
583SmoothSplineRect2DLazyDataImpl<DstT, SrcT>::SmoothSplineRect2DLazyDataImpl(const shared_ptr<const RectangularMesh2D>& src_mesh,
584 const DataVector<const SrcT>& src_vec,
585 const shared_ptr<const MeshD<2>>& dst_mesh,
586 const InterpolationFlags& flags):
587 SplineRect2DLazyDataImpl<DstT, SrcT>(src_mesh, src_vec, dst_mesh, flags)
588{
589 const size_t n0 = int(src_mesh->axis[0]->size()), n1 = int(src_mesh->axis[1]->size());
590
591 if (n0 == 0 || n1 == 0)
592 throw BadMesh("interpolate", "source mesh empty");
593
594 size_t stride0 = src_mesh->index(1, 0),
595 stride1 = src_mesh->index(0, 1);
596
597 DataVector<double> data;
598
599 if (n0 > 1) {
600 std::copy(src_vec.begin(), src_vec.end(), this->diff0.begin());
601 spline::computeDiffs<SrcT>(this->diff0.data(), stride0, stride1, src_mesh->axis[1]->size(), 0, 1, 0, src_mesh->axis[0], flags);
602 } else {
603 std::fill(this->diff0.begin(), this->diff0.end(), Zero<SrcT>());
604 }
605 if (n1 > 1) {
606 std::copy(src_vec.begin(), src_vec.end(), this->diff1.begin());
607 spline::computeDiffs<SrcT>(this->diff1.data(), stride1, stride0, src_mesh->axis[0]->size(), 0, 1, 1, src_mesh->axis[1], flags);
608 } else {
609 std::fill(this->diff1.begin(), this->diff1.end(), Zero<SrcT>());
610 }
611}
612
613
614template <typename DstT, typename SrcT>
615SmoothSplineRect3DLazyDataImpl<DstT, SrcT>::SmoothSplineRect3DLazyDataImpl(const shared_ptr<const RectilinearMesh3D>& src_mesh,
616 const DataVector<const SrcT>& src_vec,
617 const shared_ptr<const MeshD<3>>& dst_mesh,
618 const InterpolationFlags& flags):
619 SplineRect3DLazyDataImpl<DstT, SrcT>(src_mesh, src_vec, dst_mesh, flags)
620{
621 const size_t n0 = int(src_mesh->axis[0]->size()), n1 = int(src_mesh->axis[1]->size()), n2 = int(src_mesh->axis[2]->size());
622
623 if (n0 == 0 || n1 == 0)
624 throw BadMesh("interpolate", "source mesh empty");
625
626 size_t stride0 = src_mesh->index(1, 0, 0),
627 stride1 = src_mesh->index(0, 1, 0),
628 stride2 = src_mesh->index(0, 0, 1);
629
630 DataVector<double> data;
631
632 if (n0 > 1) {
633 std::copy(src_vec.begin(), src_vec.end(), this->diff0.begin());
634 spline::computeDiffs<SrcT>(this->diff0.data(), stride0,
635 stride1, src_mesh->axis[1]->size(), stride2, src_mesh->axis[2]->size(),
636 0, src_mesh->axis[0], flags);
637 } else {
638 std::fill(this->diff0.begin(), this->diff0.end(), Zero<SrcT>());
639 }
640 if (n1 > 1) {
641 std::copy(src_vec.begin(), src_vec.end(), this->diff1.begin());
642 spline::computeDiffs<SrcT>(this->diff1.data(), stride1,
643 stride0, src_mesh->axis[0]->size(), stride2, src_mesh->axis[2]->size(),
644 1, src_mesh->axis[1], flags);
645 } else {
646 std::fill(this->diff1.begin(), this->diff1.end(), Zero<SrcT>());
647 }
648 if (n2 > 1) {
649 std::copy(src_vec.begin(), src_vec.end(), this->diff2.begin());
650 spline::computeDiffs<SrcT>(this->diff2.data(), stride2,
651 stride0, src_mesh->axis[0]->size(), stride1, src_mesh->axis[1]->size(),
652 2, src_mesh->axis[2], flags);
653 } else {
654 std::fill(this->diff2.begin(), this->diff2.end(), Zero<SrcT>());
655 }
656}
657
658
659template struct PLASK_API SmoothSplineRect2DLazyDataImpl<double, double>;
660template struct PLASK_API SmoothSplineRect2DLazyDataImpl<dcomplex, dcomplex>;
661
662template struct PLASK_API SmoothSplineRect2DLazyDataImpl<Vec<2,double>, Vec<2,double>>;
663template struct PLASK_API SmoothSplineRect2DLazyDataImpl<Vec<2,dcomplex>, Vec<2,dcomplex>>;
664
665template struct PLASK_API SmoothSplineRect2DLazyDataImpl<Vec<3,double>, Vec<3,double>>;
666template struct PLASK_API SmoothSplineRect2DLazyDataImpl<Vec<3,dcomplex>, Vec<3,dcomplex>>;
667
668template struct PLASK_API SmoothSplineRect2DLazyDataImpl<Tensor2<double>, Tensor2<double>>;
669template struct PLASK_API SmoothSplineRect2DLazyDataImpl<Tensor2<dcomplex>, Tensor2<dcomplex>>;
670
671template struct PLASK_API SmoothSplineRect2DLazyDataImpl<Tensor3<double>, Tensor3<double>>;
672template struct PLASK_API SmoothSplineRect2DLazyDataImpl<Tensor3<dcomplex>, Tensor3<dcomplex>>;
673
674
675template struct PLASK_API SmoothSplineRect3DLazyDataImpl<double, double>;
676template struct PLASK_API SmoothSplineRect3DLazyDataImpl<dcomplex, dcomplex>;
677
678template struct PLASK_API SmoothSplineRect3DLazyDataImpl<Vec<2,double>, Vec<2,double>>;
679template struct PLASK_API SmoothSplineRect3DLazyDataImpl<Vec<2,dcomplex>, Vec<2,dcomplex>>;
680
681template struct PLASK_API SmoothSplineRect3DLazyDataImpl<Vec<3,double>, Vec<3,double>>;
682template struct PLASK_API SmoothSplineRect3DLazyDataImpl<Vec<3,dcomplex>, Vec<3,dcomplex>>;
683
684template struct PLASK_API SmoothSplineRect3DLazyDataImpl<Tensor2<double>, Tensor2<double>>;
685template struct PLASK_API SmoothSplineRect3DLazyDataImpl<Tensor2<dcomplex>, Tensor2<dcomplex>>;
686
687template struct PLASK_API SmoothSplineRect3DLazyDataImpl<Tensor3<double>, Tensor3<double>>;
688template struct PLASK_API SmoothSplineRect3DLazyDataImpl<Tensor3<dcomplex>, Tensor3<dcomplex>>;
689
690
691
692} // namespace plask