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