OPAL (Object Oriented Parallel Accelerator Library) 2022.1
OPAL
AmrLagrangeInterpolater.hpp
Go to the documentation of this file.
1//
2// Class AmrLagrangeInterpolater
3// Lagrange interpolation for coarse-fine interfaces.
4//
5// Copyright (c) 2017 - 2020, Matthias Frey, Paul Scherrer Institut, Villigen PSI, Switzerland
6// All rights reserved
7//
8// Implemented as part of the PhD thesis
9// "Precise Simulations of Multibunches in High Intensity Cyclotrons"
10//
11// This file is part of OPAL.
12//
13// OPAL is free software: you can redistribute it and/or modify
14// it under the terms of the GNU General Public License as published by
15// the Free Software Foundation, either version 3 of the License, or
16// (at your option) any later version.
17//
18// You should have received a copy of the GNU General Public License
19// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
20//
21#include "Utilities/OpalException.h"
22
23#if AMREX_SPACEDIM == 3
24template <class Level>
25constexpr typename AmrLagrangeInterpolater<Level>::qpattern_t
26 AmrLagrangeInterpolater<Level>::qpattern_ms;
27
28
29template <class Level>
30constexpr typename AmrLagrangeInterpolater<Level>::lpattern_t
31 AmrLagrangeInterpolater<Level>::lpattern_ms;
32#endif
33
34// y_b y_t
35template <class Level>
36const typename AmrLagrangeInterpolater<Level>::scalar_t
37 AmrLagrangeInterpolater<Level>::lookup1a_ms[2] = {0.25, -0.25};
38
39template <class Level>
40const typename AmrLagrangeInterpolater<Level>::scalar_t
41 AmrLagrangeInterpolater<Level>::lookup2a_ms[2] = {0.75, 1.25};
42
43template <class Level>
44const typename AmrLagrangeInterpolater<Level>::scalar_t
45 AmrLagrangeInterpolater<Level>::lookup1b_ms[2] = {-0.25, 0.25};
46
47template <class Level>
48const typename AmrLagrangeInterpolater<Level>::scalar_t
49 AmrLagrangeInterpolater<Level>::lookup2b_ms[2] = {1.25, 0.75};
50
51#if AMREX_SPACEDIM == 3
52template <class Level>
53const typename AmrLagrangeInterpolater<Level>::scalar_t
54 AmrLagrangeInterpolater<Level>::lookup3_ms[2] = {5.0 / 32.0, -3.0 / 32.0 };
55
56template <class Level>
57const typename AmrLagrangeInterpolater<Level>::scalar_t
58 AmrLagrangeInterpolater<Level>::lookup3r_ms[2] = {-3.0 / 32.0, 5.0 / 32.0 };
59
60template <class Level>
61const typename AmrLagrangeInterpolater<Level>::scalar_t
62 AmrLagrangeInterpolater<Level>::lookup4_ms[2] = {7.0 / 16.0, -9.0 / 16.0 };
63
64template <class Level>
65const typename AmrLagrangeInterpolater<Level>::scalar_t
66 AmrLagrangeInterpolater<Level>::lookup4r_ms[2] = {-9.0 / 16.0, 7.0 / 16.0 };
67
68template <class Level>
69const typename AmrLagrangeInterpolater<Level>::scalar_t
70 AmrLagrangeInterpolater<Level>::lookup5_ms[2] = {45.0 / 32.0, 21.0 / 32.0};
71
72template <class Level>
73const typename AmrLagrangeInterpolater<Level>::scalar_t
74 AmrLagrangeInterpolater<Level>::lookup5r_ms[2] = {21.0 / 32.0, 45.0 / 32.0};
75
76template <class Level>
77const typename AmrLagrangeInterpolater<Level>::scalar_t
78 AmrLagrangeInterpolater<Level>::lookup6_ms = 15.0 / 16.0;
79
80template <class Level>
81const typename AmrLagrangeInterpolater<Level>::scalar_t
82 AmrLagrangeInterpolater<Level>::factor_ms = 8.0 / 15.0;
83#endif
84
85
86
87
88template <class Level>
89AmrLagrangeInterpolater<Level>::AmrLagrangeInterpolater(Order order)
90 : AmrInterpolater<Level>( lo_t(order) + 1 )
91{ }
92
93
94template <class Level>
95void AmrLagrangeInterpolater<Level>::stencil(
96 const AmrIntVect_t& /*iv*/,
97 const basefab_t& /*fab*/,
98 typename Level::umap_t& /*map*/,
99 const typename Level::scalar_t& /*scale*/,
100 Level* /*mglevel*/)
101{
102
103}
104
105
106template <class Level>
107void AmrLagrangeInterpolater<Level>::coarse(
108 const AmrIntVect_t& iv,
109 typename Level::umap_t& map,
110 const typename Level::scalar_t& scale,
111 lo_t dir, lo_t shift, const basefab_t& rfab,
112 const AmrIntVect_t& riv,
113 Level* mglevel)
114{
115 // polynomial degree = #points - 1
116 switch ( this->nPoints_m - 1 ) {
117
118 case Order::QUADRATIC:
119 this->crseQuadratic_m(iv, map, scale, dir, shift, rfab, riv, mglevel);
120 break;
121 case Order::LINEAR:
122 this->crseLinear_m(iv, map, scale, dir, shift, rfab, riv, mglevel);
123 break;
124 default:
125 throw OpalException("AmrLagrangeInterpolater::coarse()",
126 "Not implemented interpolation");
127 }
128}
129
130
131template <class Level>
132void AmrLagrangeInterpolater<Level>::fine(
133 const AmrIntVect_t& iv,
134 typename Level::umap_t& map,
135 const typename Level::scalar_t& scale,
136 lo_t dir, lo_t shift,
137 Level* mglevel)
138{
139 // polynomial degree = #points - 1
140 switch ( this->nPoints_m - 1 ) {
141
142 case Order::QUADRATIC:
143 this->fineQuadratic_m(iv, map, scale, dir, shift, mglevel);
144 break;
145 case Order::LINEAR:
146 this->fineLinear_m(iv, map, scale, dir, shift, mglevel);
147 break;
148 default:
149 throw OpalException("AmrLagrangeInterpolater::fine()",
150 "Not implemented interpolation");
151 }
152}
153
154
155template <class Level>
156void AmrLagrangeInterpolater<Level>::fineLinear_m(
157 const AmrIntVect_t& iv,
158 typename Level::umap_t& map,
159 const typename Level::scalar_t& scale,
160 lo_t dir, lo_t shift,
161 Level* mglevel)
162{
163 /*
164 * computes the ghost cell directly
165 */
166 AmrIntVect_t tmp = iv;
167 // first fine cell on refined coarse cell (closer to interface)
168 tmp[dir] += shift;
169 map[mglevel->serialize(tmp)] += 2.0 * scale;
170
171 // second fine cell on refined coarse cell (further away from interface)
172 tmp[dir] += shift;
173 map[mglevel->serialize(tmp)] -= scale;
174}
175
176
177template <class Level>
178void AmrLagrangeInterpolater<Level>::fineQuadratic_m(
179 const AmrIntVect_t& iv,
180 typename Level::umap_t& map,
181 const typename Level::scalar_t& scale,
182 lo_t dir, lo_t shift,
183 Level* mglevel)
184{
185 AmrIntVect_t tmp = iv;
186 // first fine cell on refined coarse cell (closer to interface)
187 tmp[dir] += shift;
188 map[mglevel->serialize(tmp)] += 2.0 / 3.0 * scale;
189
190 // second fine cell on refined coarse cell (further away from interface)
191 tmp[dir] += shift;
192 map[mglevel->serialize(tmp)] -= 0.2 * scale;
193}
194
195
196template <class Level>
197void AmrLagrangeInterpolater<Level>::crseLinear_m(
198 const AmrIntVect_t& iv,
199 typename Level::umap_t& map,
200 const typename Level::scalar_t& scale,
201 lo_t dir, lo_t /*shift*/, const basefab_t& rfab,
202 const AmrIntVect_t& riv,
203 Level* mglevel)
204{
205#if AMREX_SPACEDIM == 2
206 bool top = (riv[(dir+1)%AMREX_SPACEDIM] % 2 == 1);
207
208 // right / upper / back
209 AmrIntVect_t niv = iv;
210 niv[(dir+1)%AMREX_SPACEDIM ] += 1;
211
212 // left / lower / front
213 AmrIntVect_t miv = iv;
214 miv[(dir+1)%AMREX_SPACEDIM ] -= 1;
215
216 // factor for fine
217 scalar_t fac = 8.0 / 15.0 * scale;
218
219 if ( mglevel->isValid(niv) && rfab(niv) != Level::Refined::YES ) {
220 // check r / u / b --> 1: valid; 0: not valid
221 map[mglevel->serialize(iv)] += fac * lookup2a_ms[top];
222 map[mglevel->serialize(niv)] += fac * lookup1a_ms[top];
223
224 } else if ( mglevel->isValid(miv) && rfab(miv) != Level::Refined::YES ) {
225 // check l / f --> 1: valid; 0: not valid
226 map[mglevel->serialize(iv)] += fac * lookup2a_ms[top];
227 map[mglevel->serialize(miv)] += fac * lookup1a_ms[top];
228
229 } else
230 throw OpalException("AmrLagrangeInterpolater::crseLinear_m()",
231 "No valid interpolation scenario found!");
232
233#elif AMREX_SPACEDIM == 3
234
235 /* x y z
236 * ------------------
237 * dir: 0 1 2
238 * top1: 1 2 0 (i, L)
239 * top2: 2 0 1 (j, K)
240 */
241
242 /* There are 4 coefficients from Lagrange interpolation.
243 * Those are given by the product of one of
244 * L0, L1 and one of K0, K1.
245 *
246 * g(x, y) = f(x0, y0) * L0(x) * K0(y) +
247 * f(x0, y1) * L0(x) * K1(y) +
248 * f(x1, y0) * L1(x) * K0(y) +
249 * f(x1, y1) * L1(x) * K1(y) +
250 */
251
252 /*
253 * check in 3x3 area (using iv as center) if 4 cells are not covered
254 */
255 lo_t d1 = (dir+1)%AMREX_SPACEDIM;
256 lo_t d2 = (dir+2)%AMREX_SPACEDIM;
257
258 lbits_t area;
259 lo_t bit = 0;
260
261 AmrIntVect_t tmp = iv;
262 for (int i = -1; i < 2; ++i) {
263 tmp[d1] += i;
264 for (int j = -1; j < 2; ++j) {
265
266 tmp[d2] += j;
267 // make use of "short-circuit evaluation"
268 area[bit] = ( mglevel->isValid(tmp) && rfab(tmp) );
269 ++bit;
270
271 // undo
272 tmp[d2] -= j;
273 }
274 // undo
275 tmp[d1] -= i;
276 }
277
278 qpattern_t::const_iterator pit = std::begin(this->lpattern_ms);
279
280 while ( pit != std::end(this->lpattern_ms) ) {
281 if ( *pit == (area & lbits_t(*pit)).to_ulong() )
282 break;
283 ++pit;
284 }
285
286 // factor for fine
287 scalar_t fac = factor_ms * scale;
288
289 scalar_t L[2] = {0.0, 0.0};
290 lo_t top1 = riv[d1] % 2;
291
292 scalar_t K[2] = {0.0, 0.0};
293 lo_t top2 = riv[d2] % 2;
294
295 lo_t begin[2] = { 0, 0 };
296 lo_t end[2] = { 0, 0 };
297
298 switch ( *pit ) {
299 case this->lpattern_ms[0]:
300 {
301 // corner top right pattern
302 L[0] = lookup1a_ms[top1]; // L_{-1}
303 L[1] = lookup2a_ms[top1]; // L_{0}
304 begin[0] = -1;
305 end[0] = 0;
306
307 K[0] = lookup1a_ms[top2]; // K_{-1}
308 K[1] = lookup2a_ms[top2]; // K_{0}
309 begin[1] = -1;
310 end[1] = 0;
311 break;
312 }
313 case this->lpattern_ms[1]:
314 {
315 // corner bottom right pattern
316 L[0] = lookup2b_ms[top1]; // L_{0}
317 L[1] = lookup1b_ms[top1]; // L_{1}
318 begin[0] = 0;
319 end[0] = 1;
320
321 K[0] = lookup1a_ms[top2]; // K_{-1}
322 K[1] = lookup2a_ms[top2]; // K_{0}
323 begin[1] = -1;
324 end[1] = 0;
325 break;
326 }
327 case this->lpattern_ms[2]:
328 {
329 // corner bottom left pattern
330 L[0] = lookup2b_ms[top1]; // L_{0}
331 L[1] = lookup1b_ms[top1]; // L_{1}
332 begin[0] = 0;
333 end[0] = 1;
334
335 K[0] = lookup2b_ms[top2]; // K_{0}
336 K[1] = lookup1b_ms[top2]; // K_{1}
337 begin[1] = 0;
338 end[1] = 1;
339 break;
340 }
341 case this->lpattern_ms[3]:
342 {
343 // corner top left pattern
344 L[0] = lookup1a_ms[top1]; // L_{-1}
345 L[1] = lookup2a_ms[top1]; // L_{0}
346 begin[0] = -1;
347 end[0] = 0;
348
349 K[0] = lookup2b_ms[top2]; // K_{0}
350 K[1] = lookup1b_ms[top2]; // K_{1}
351 begin[1] = 0;
352 end[1] = 1;
353 break;
354 }
355 default:
356 throw OpalException("AmrLagrangeInterpolater::crseLinear_m()",
357 "No valid interpolation scenario found!");
358 }
359
360 /*
361 * if pattern is known --> add stencil
362 */
363 AmrIntVect_t niv = iv;
364 for (int i = begin[0]; i <= end[0]; ++i) {
365 niv[d1] += i;
366 for (int j = begin[1]; j <= end[1]; ++j) {
367 niv[d2] += j;
368
369 scalar_t value = fac * L[i-begin[0]] * K[j-begin[1]];
370 if ( !mglevel->applyBoundary(niv, rfab, map, value) )
371 map[mglevel->serialize(niv)] += value;
372
373 // undo
374 niv[d2] -= j;
375 }
376 // undo
377 niv[d1] -= i;
378 }
379#else
380 #error Lagrange interpolation: Only 2D and 3D are supported!
381#endif
382 // the neighbour cancels out
383}
384
385
386template <class Level>
387void AmrLagrangeInterpolater<Level>::crseQuadratic_m(
388 const AmrIntVect_t& iv,
389 typename Level::umap_t& map,
390 const typename Level::scalar_t& scale,
391 lo_t dir, lo_t shift, const basefab_t& rfab,
392 const AmrIntVect_t& riv,
393 Level* mglevel)
394{
395#if AMREX_SPACEDIM == 2
396
397 bool top = (riv[(dir+1)%AMREX_SPACEDIM] % 2 == 1);
398
399 // right / upper / back
400 AmrIntVect_t niv = iv;
401 niv[(dir+1)%AMREX_SPACEDIM ] += 1;
402
403 // left / lower / front
404 AmrIntVect_t miv = iv;
405 miv[(dir+1)%AMREX_SPACEDIM ] -= 1;
406
407 // 2nd right / upper / back
408 AmrIntVect_t n2iv = niv;
409 n2iv[(dir+1)%AMREX_SPACEDIM ] += 1;
410
411 // 2nd left / lower / front
412 AmrIntVect_t m2iv = miv;
413 m2iv[(dir+1)%AMREX_SPACEDIM ] -= 1;
414
415 /* 3 cases:
416 * --------
417 * r: right neighbour of iv (center of Laplacian)
418 * u: upper neighbour of iv
419 * b: back neighbour of iv
420 *
421 * l: lower / left neighbour of iv
422 * f: front neighbour of iv
423 *
424 * -1 --> not valid, error happend
425 * 0 --> r / u / b and 2nd r / u / b are valid
426 * 1 --> direct neighbours (r / u / b and l / f) of iv are valid (has priority)
427 * 2 --> l / f and 2nd l / f are valid
428 */
429
430 // check r / u / b --> 1: valid; 0: not valid
431 bool rub = rfab(niv);
432
433 // check l / f --> 1: valid; 0: not valid
434 bool lf = rfab(miv);
435
436 // check 2nd r / u / b
437 bool rub2 = rfab(n2iv);
438
439 // check 2nd l / f
440 bool lf2 = rfab(m2iv);
441
442 if ( rub && lf )
443 {
444 /*
445 * standard case -1, +1 are not-refined nor at physical/mesh boundary
446 */
447 // cell is not refined and not at physical boundary
448
449 // y_t or y_b
450 map[mglevel->serialize(iv)] += 0.5 * scale;
451
452 // y_t y_b
453 scalar_t value = scale * ((top) ? 1.0 / 12.0 : -0.05);
454 if ( !mglevel->applyBoundary(niv, rfab, map, value) )
455 map[mglevel->serialize(niv)] += value;
456
457 // y_t y_b
458 value = scale * ((top) ? -0.05 : 1.0 / 12.0);
459 if ( !mglevel->applyBoundary(miv, rfab, map, value) )
460 map[mglevel->serialize(miv)] += value;
461
462 } else if ( rub && rub2 ) {
463 /*
464 * corner case --> right / upper / back + 2nd right / upper / back
465 */
466 // y_t y_b
467 scalar_t value = scale * ((top) ? 7.0 / 20.0 : 0.75);
468 map[mglevel->serialize(iv)] += value;
469
470 // y_t y_b
471 value = scale * ((top) ? 7.0 / 30.0 : -0.3);
472 if ( !mglevel->applyBoundary(niv, rfab, map, value) )
473 map[mglevel->serialize(niv)] += value;
474
475 // y_t y_b
476 value = scale * ((top) ? -0.05 : 1.0 / 12.0);
477 if ( !mglevel->applyBoundary(n2iv, rfab, map, value) )
478 map[mglevel->serialize(n2iv)] += value;
479
480 } else if ( lf && lf2 ) {
481 /*
482 * corner case --> left / lower / front + 2nd left / lower / front
483 */
484 // y_t y_b
485 scalar_t value = scale * ((top) ? 0.75 : 7.0 / 20.0);
486 map[mglevel->serialize(iv)] += value;
487
488 // y_t y_b
489 value = scale * ((top) ? -0.3 : 7.0 / 30);
490 if ( !mglevel->applyBoundary(miv, rfab, map, value) )
491 map[mglevel->serialize(miv)] += value;
492
493 // y_t y_b
494 value = scale * ((top) ? 1.0 / 12.0 : -0.05);
495 if ( !mglevel->applyBoundary(m2iv, rfab, map, value) )
496 map[mglevel->serialize(m2iv)] += value;
497
498 } else {
499 /* last trial: linear Lagrange interpolation
500 * --> it throws an error if not possible
501 */
502 this->crseLinear_m(iv, map, scale, dir, shift, rfab, riv, mglevel);
503 }
504
505#elif AMREX_SPACEDIM == 3
506
507 /* x y z
508 * ------------------
509 * dir: 0 1 2
510 * top1: 1 2 0 (i, L)
511 * top2: 2 0 1 (j, K)
512 */
513
514 /* There are 9 coefficients from Lagrange interpolation.
515 * Those are given by the product of one of
516 * L0, L1, L2 and one of K0, K1, K2.
517 *
518 * g(x, y) = f(x0, y0) * L0(x) * K0(y) +
519 * f(x0, y1) * L0(x) * K1(y) +
520 * f(x0, y2) * L0(x) * K2(y) +
521 * f(x1, y0) * L1(x) * K0(y) +
522 * f(x1, y1) * L1(x) * K1(y) +
523 * f(x1, y2) * L1(x) * K2(y) +
524 * f(x2, y0) * L2(x) * K0(y) +
525 * f(x2, y1) * L2(x) * K1(y) +
526 * f(x2, y2) * L2(x) * K2(y) +
527 */
528
529
530 /*
531 * check in 5x5 area (using iv as center) if 9 cells are not covered
532 */
533 lo_t d1 = (dir+1)%AMREX_SPACEDIM;
534 lo_t d2 = (dir+2)%AMREX_SPACEDIM;
535
536 qbits_t area;
537 lo_t bit = 0;
538
539 AmrIntVect_t tmp = iv;
540 for (int i = -2; i < 3; ++i) {
541 tmp[d1] += i;
542 for (int j = -2; j < 3; ++j) {
543
544 tmp[d2] += j;
545
546 area[bit] = ( mglevel->isValid(tmp) && rfab(tmp) );
547 ++bit;
548
549 // undo
550 tmp[d2] -= j;
551 }
552 // undo
553 tmp[d1] -= i;
554 }
555
556 qpattern_t::const_iterator pit = std::begin(this->qpattern_ms);
557
558 while ( pit != std::end(this->qpattern_ms) ) {
559 if ( *pit == (area & qbits_t(*pit)).to_ulong() )
560 break;
561 ++pit;
562 }
563
564 // factor for fine
565 scalar_t fac = factor_ms * scale;
566
567 scalar_t L[3] = {0.0, 0.0, 0.0};
568 lo_t top1 = riv[d1] % 2;
569
570 scalar_t K[3] = {0.0, 0.0, 0.0};
571 lo_t top2 = riv[d2] % 2;
572
573 lo_t begin[2] = { 0, 0 };
574 lo_t end[2] = { 0, 0 };
575
576 switch ( *pit ) {
577 case this->qpattern_ms[0]:
578 {
579 // cross pattern
580 L[0] = lookup3_ms[top1]; // L_{-1}
581 L[1] = lookup6_ms; // L_{0}
582 L[2] = lookup3r_ms[top1]; // L_{1}
583 begin[0] = -1;
584 end[0] = 1;
585
586 K[0] = lookup3_ms[top2]; // K_{-1}
587 K[1] = lookup6_ms; // K_{0}
588 K[2] = lookup3r_ms[top2]; // K_{1}
589 begin[1] = -1;
590 end[1] = 1;
591 break;
592 }
593 case this->qpattern_ms[1]:
594 {
595 // T pattern
596 L[0] = lookup3r_ms[top1]; // L_{-2}
597 L[1] = lookup4_ms[top1]; // L_{-1}
598 L[2] = lookup5r_ms[top1]; // L_{0}
599 begin[0] = -2;
600 end[0] = 0;
601
602 K[0] = lookup3_ms[top2]; // K_{-1}
603 K[1] = lookup6_ms; // K_{0}
604 K[2] = lookup3r_ms[top2]; // K_{1}
605 begin[1] = -1;
606 end[1] = 1;
607 break;
608 }
609 case this->qpattern_ms[2]:
610 {
611 // right hammer pattern
612 L[0] = lookup3_ms[top1]; // L_{-1}
613 L[1] = lookup6_ms; // L_{0}
614 L[2] = lookup3r_ms[top1]; // L_{1}
615 begin[0] = -1;
616 end[0] = 1;
617
618 K[0] = lookup3r_ms[top2]; // K_{-2}
619 K[1] = lookup4_ms[top2]; // K_{-1}
620 K[2] = lookup5r_ms[top2]; // K_{0}
621 begin[1] = -2;
622 end[1] = 0;
623 break;
624 }
625 case this->qpattern_ms[3]:
626 {
627 // T on head pattern
628 L[0] = lookup5_ms[top1]; // L_{0}
629 L[1] = lookup4r_ms[top1]; // L_{1}
630 L[2] = lookup3_ms[top1]; // L_{2}
631 begin[0] = 0;
632 end[0] = 2;
633
634 K[0] = lookup3_ms[top2]; // K_{-1}
635 K[1] = lookup6_ms; // K_{0}
636 K[2] = lookup3r_ms[top2]; // K_{1}
637 begin[1] = -1;
638 end[1] = 1;
639 break;
640 }
641 case this->qpattern_ms[4]:
642 {
643 // left hammer pattern
644 L[0] = lookup3_ms[top1]; // L_{-1}
645 L[1] = lookup6_ms; // L_{0}
646 L[2] = lookup3r_ms[top1]; // L_{1}
647 begin[0] = -1;
648 end[0] = 1;
649
650 K[0] = lookup5_ms[top2]; // K_{0}
651 K[1] = lookup4r_ms[top2]; // K_{1}
652 K[2] = lookup3_ms[top2]; // K_{2}
653 begin[1] = 0;
654 end[1] = 2;
655 break;
656 }
657 case this->qpattern_ms[5]:
658 {
659 // upper left corner pattern
660 L[0] = lookup3r_ms[top1]; // L_{-2}
661 L[1] = lookup4_ms[top1]; // L_{-1}
662 L[2] = lookup5r_ms[top1]; // L_{0}
663 begin[0] = -2;
664 end[0] = 0;
665
666 K[0] = lookup5_ms[top2]; // K_{0}
667 K[1] = lookup4r_ms[top2]; // K_{1}
668 K[2] = lookup3_ms[top2]; // K_{2}
669 begin[1] = 0;
670 end[1] = 2;
671 break;
672 }
673 case this->qpattern_ms[6]:
674 {
675 // upper right corner pattern
676 L[0] = lookup3r_ms[top1]; // L_{-2}
677 L[1] = lookup4_ms[top1]; // L_{-1}
678 L[2] = lookup5r_ms[top1]; // L_{0}
679 begin[0] = -2;
680 end[0] = 0;
681
682 K[0] = lookup3r_ms[top2]; // K_{-2}
683 K[1] = lookup4_ms[top2]; // K_{-1}
684 K[2] = lookup5r_ms[top2]; // K_{0}
685 begin[1] = -2;
686 end[1] = 0;
687 break;
688 }
689 case this->qpattern_ms[7]:
690 {
691 // mirrored L pattern
692 L[0] = lookup5_ms[top1]; // L_{0}
693 L[1] = lookup4r_ms[top1]; // L_{1}
694 L[2] = lookup3_ms[top1]; // L_{2}
695 begin[0] = 0;
696 end[0] = 2;
697
698 K[0] = lookup3r_ms[top2]; // K_{-2}
699 K[1] = lookup4_ms[top2]; // K_{-1}
700 K[2] = lookup5r_ms[top2]; // K_{0}
701 begin[1] = -2;
702 end[1] = 0;
703 break;
704 }
705 case this->qpattern_ms[8]:
706 {
707 // L pattern
708 L[0] = lookup5_ms[top1]; // L_{0}
709 L[1] = lookup4r_ms[top1]; // L_{1}
710 L[2] = lookup3_ms[top1]; // L_{2}
711 begin[0] = 0;
712 end[0] = 2;
713
714 K[0] = lookup5_ms[top2]; // K_{0}
715 K[1] = lookup4r_ms[top2]; // K_{1}
716 K[2] = lookup3_ms[top2]; // K_{2}
717 begin[1] = 0;
718 end[1] = 2;
719 break;
720 }
721 default:
722 {
723 /* unknown pattern --> last trial: linear Lagrange interpolation
724 * --> it throws an error if not possible
725 */
726 this->crseLinear_m(iv, map, scale, dir, shift, rfab, riv, mglevel);
727 return;
728 }
729 }
730
731 /*
732 * if pattern is known --> add stencil
733 */
734 AmrIntVect_t tmp1 = iv;
735 for (int i = begin[0]; i <= end[0]; ++i) {
736 tmp1[d1] += i;
737 for (int j = begin[1]; j <= end[1]; ++j) {
738 tmp1[d2] += j;
739
740 scalar_t value = fac * L[i-begin[0]] * K[j-begin[1]];
741 if ( !mglevel->applyBoundary(tmp1, rfab, map, value) )
742 map[mglevel->serialize(tmp1)] += value;
743
744 // undo
745 tmp1[d2] -= j;
746 }
747 // undo
748 tmp1[d1] -= i;
749 }
750
751#else
752 #error Lagrange interpolation: Only 2D and 3D are supported!
753#endif
754}
end
PartBunchBase< T, Dim >::ConstIterator end(PartBunchBase< T, Dim > const &bunch)
Definition: PartBunchBase.h:736
begin
PartBunchBase< T, Dim >::ConstIterator begin(PartBunchBase< T, Dim > const &bunch)
Definition: PartBunchBase.h:731
QUADRATIC
@ QUADRATIC
Definition: IrregularDomain.h:39
LINEAR
@ LINEAR
Definition: IrregularDomain.h:38
amr::scalar_t
double scalar_t
Definition: AmrMultiGridDefs.h:37
AmrInterpolater::AmrIntVect_t
amr::AmrIntVect_t AmrIntVect_t
Definition: AmrInterpolater.h:35
AmrInterpolater::lo_t
Level::lo_t lo_t
Definition: AmrInterpolater.h:31
AmrInterpolater::basefab_t
Level::basefab_t basefab_t
Definition: AmrInterpolater.h:34
AmrLagrangeInterpolater::fine
void fine(const AmrIntVect_t &iv, umap_t &map, const scalar_t &scale, lo_t dir, lo_t shift, Level *mglevel)
Definition: AmrLagrangeInterpolater.hpp:132
AmrLagrangeInterpolater::crseLinear_m
void crseLinear_m(const AmrIntVect_t &iv, umap_t &map, const scalar_t &scale, lo_t dir, lo_t shift, const basefab_t &rfab, const AmrIntVect_t &riv, Level *mglevel)
Definition: AmrLagrangeInterpolater.hpp:197
AmrLagrangeInterpolater::fineLinear_m
void fineLinear_m(const AmrIntVect_t &iv, umap_t &map, const scalar_t &scale, lo_t dir, lo_t shift, Level *mglevel)
Definition: AmrLagrangeInterpolater.hpp:156
AmrLagrangeInterpolater::crseQuadratic_m
void crseQuadratic_m(const AmrIntVect_t &iv, umap_t &map, const scalar_t &scale, lo_t dir, lo_t shift, const basefab_t &rfab, const AmrIntVect_t &riv, Level *mglevel)
Definition: AmrLagrangeInterpolater.hpp:387
AmrLagrangeInterpolater::coarse
void coarse(const AmrIntVect_t &iv, umap_t &map, const scalar_t &scale, lo_t dir, lo_t shift, const basefab_t &rfab, const AmrIntVect_t &riv, Level *mglevel)
Definition: AmrLagrangeInterpolater.hpp:107
AmrLagrangeInterpolater::fineQuadratic_m
void fineQuadratic_m(const AmrIntVect_t &iv, umap_t &map, const scalar_t &scale, lo_t dir, lo_t shift, Level *mglevel)
Definition: AmrLagrangeInterpolater.hpp:178
AmrLagrangeInterpolater::stencil
void stencil(const AmrIntVect_t &iv, const basefab_t &fab, umap_t &map, const scalar_t &scale, Level *mglevel)
Definition: AmrLagrangeInterpolater.hpp:95
OpalException
The base class for all OPAL exceptions.
Definition: OpalException.h:28
Generated on Thu Oct 20 2022 17:40:11 for OPAL (Object Oriented Parallel Accelerator Library) by doxygen 1.9.3