AmrLagrangeInterpolater.hpp 22.4 KB
Newer Older
1 2 3
//
// Class AmrLagrangeInterpolater
//   Lagrange interpolation for coarse-fine interfaces.
4 5 6 7
//
// Copyright (c) 2017 - 2020, Matthias Frey, Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved
//
8 9 10 11 12 13 14 15 16 17 18 19 20
// Implemented as part of the PhD thesis
// "Precise Simulations of Multibunches in High Intensity Cyclotrons"
//
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
21 22
#include "Utilities/OpalException.h"

23
#if AMREX_SPACEDIM == 3
frey_m's avatar
frey_m committed
24 25 26
template <class Level>
constexpr typename AmrLagrangeInterpolater<Level>::qpattern_t
    AmrLagrangeInterpolater<Level>::qpattern_ms;
27 28


frey_m's avatar
frey_m committed
29 30 31
template <class Level>
constexpr typename AmrLagrangeInterpolater<Level>::lpattern_t
    AmrLagrangeInterpolater<Level>::lpattern_ms;
32 33
#endif

34
//                                                      y_b   y_t
frey_m's avatar
frey_m committed
35 36
template <class Level>
const typename AmrLagrangeInterpolater<Level>::scalar_t
37
    AmrLagrangeInterpolater<Level>::lookup1a_ms[2] = {0.25, -0.25};
38

frey_m's avatar
frey_m committed
39 40
template <class Level>
const typename AmrLagrangeInterpolater<Level>::scalar_t
41 42 43 44 45 46 47 48 49
    AmrLagrangeInterpolater<Level>::lookup2a_ms[2] = {0.75, 1.25};

template <class Level>
const typename AmrLagrangeInterpolater<Level>::scalar_t
    AmrLagrangeInterpolater<Level>::lookup1b_ms[2] = {-0.25, 0.25};

template <class Level>
const typename AmrLagrangeInterpolater<Level>::scalar_t
    AmrLagrangeInterpolater<Level>::lookup2b_ms[2] = {1.25, 0.75};
50 51

#if AMREX_SPACEDIM == 3
frey_m's avatar
frey_m committed
52 53 54
template <class Level>
const typename AmrLagrangeInterpolater<Level>::scalar_t
    AmrLagrangeInterpolater<Level>::lookup3_ms[2] = {5.0 / 32.0, -3.0 / 32.0 };
55

frey_m's avatar
frey_m committed
56 57 58
template <class Level>
const typename AmrLagrangeInterpolater<Level>::scalar_t
    AmrLagrangeInterpolater<Level>::lookup3r_ms[2] = {-3.0 / 32.0, 5.0 / 32.0 };
59

frey_m's avatar
frey_m committed
60 61 62
template <class Level>
const typename AmrLagrangeInterpolater<Level>::scalar_t
    AmrLagrangeInterpolater<Level>::lookup4_ms[2] = {7.0 / 16.0, -9.0 / 16.0 };
63

frey_m's avatar
frey_m committed
64 65 66
template <class Level>
const typename AmrLagrangeInterpolater<Level>::scalar_t
    AmrLagrangeInterpolater<Level>::lookup4r_ms[2] = {-9.0 / 16.0, 7.0 / 16.0 };
67

frey_m's avatar
frey_m committed
68 69 70
template <class Level>
const typename AmrLagrangeInterpolater<Level>::scalar_t
    AmrLagrangeInterpolater<Level>::lookup5_ms[2] = {45.0 / 32.0, 21.0 / 32.0};
71

frey_m's avatar
frey_m committed
72 73 74
template <class Level>
const typename AmrLagrangeInterpolater<Level>::scalar_t
    AmrLagrangeInterpolater<Level>::lookup5r_ms[2] = {21.0 / 32.0, 45.0 / 32.0};
75

frey_m's avatar
frey_m committed
76 77 78
template <class Level>
const typename AmrLagrangeInterpolater<Level>::scalar_t
    AmrLagrangeInterpolater<Level>::lookup6_ms = 15.0 / 16.0;
79

frey_m's avatar
frey_m committed
80 81 82
template <class Level>
const typename AmrLagrangeInterpolater<Level>::scalar_t
    AmrLagrangeInterpolater<Level>::factor_ms = 8.0 / 15.0;
83 84 85 86 87
#endif




frey_m's avatar
frey_m committed
88 89 90
template <class Level>
AmrLagrangeInterpolater<Level>::AmrLagrangeInterpolater(Order order)
    : AmrInterpolater<Level>( lo_t(order) + 1 )
91 92 93
{ }


frey_m's avatar
frey_m committed
94 95
template <class Level>
void AmrLagrangeInterpolater<Level>::stencil(
96
    const AmrIntVect_t& iv,
97
    const basefab_t& fab,
frey_m's avatar
frey_m committed
98 99 100
    typename Level::umap_t& map,
    const typename Level::scalar_t& scale,
    Level* mglevel)
101 102 103 104 105
{
    
}


frey_m's avatar
frey_m committed
106 107
template <class Level>
void AmrLagrangeInterpolater<Level>::coarse(
108
    const AmrIntVect_t& iv,
frey_m's avatar
frey_m committed
109 110
    typename Level::umap_t& map,
    const typename Level::scalar_t& scale,
111
    lo_t dir, lo_t shift, const basefab_t& rfab,
112
    const AmrIntVect_t& riv,
frey_m's avatar
frey_m committed
113
    Level* mglevel)
114 115 116 117 118
{
    // polynomial degree = #points - 1
    switch ( this->nPoints_m - 1 ) {
        
        case Order::QUADRATIC:
119 120 121 122
            this->crseQuadratic_m(iv, map, scale, dir, shift, rfab, riv, mglevel);
            break;
        case Order::LINEAR:
            this->crseLinear_m(iv, map, scale, dir, shift, rfab, riv, mglevel);
123 124
            break;
        default:
125 126
            throw OpalException("AmrLagrangeInterpolater::coarse()",
                                "Not implemented interpolation");
127 128 129 130
    }
}


frey_m's avatar
frey_m committed
131 132
template <class Level>
void AmrLagrangeInterpolater<Level>::fine(
133
    const AmrIntVect_t& iv,
frey_m's avatar
frey_m committed
134 135
    typename Level::umap_t& map,
    const typename Level::scalar_t& scale,
136
    lo_t dir, lo_t shift,
frey_m's avatar
frey_m committed
137
    Level* mglevel)
138 139 140 141 142
{
    // polynomial degree = #points - 1
    switch ( this->nPoints_m - 1 ) {
        
        case Order::QUADRATIC:
143 144 145 146
            this->fineQuadratic_m(iv, map, scale, dir, shift, mglevel);
            break;
        case Order::LINEAR:
            this->fineLinear_m(iv, map, scale, dir, shift, mglevel);
147 148
            break;
        default:
149 150
            throw OpalException("AmrLagrangeInterpolater::fine()",
                                "Not implemented interpolation");
151 152 153 154
    }
}


frey_m's avatar
frey_m committed
155 156
template <class Level>
void AmrLagrangeInterpolater<Level>::fineLinear_m(
157
    const AmrIntVect_t& iv,
frey_m's avatar
frey_m committed
158 159
    typename Level::umap_t& map,
    const typename Level::scalar_t& scale,
160
    lo_t dir, lo_t shift,
frey_m's avatar
frey_m committed
161
    Level* mglevel)
162 163 164 165 166 167 168
{
    /*
     * computes the ghost cell directly
     */
    AmrIntVect_t tmp = iv;
    // first fine cell on refined coarse cell (closer to interface)
    tmp[dir] += shift;
169
    map[mglevel->serialize(tmp)] += 2.0 * scale;
170 171 172
    
    // second fine cell on refined coarse cell (further away from interface)
    tmp[dir] += shift;
173
    map[mglevel->serialize(tmp)] -= scale;
174 175 176
}


frey_m's avatar
frey_m committed
177 178
template <class Level>
void AmrLagrangeInterpolater<Level>::fineQuadratic_m(
179
    const AmrIntVect_t& iv,
frey_m's avatar
frey_m committed
180 181
    typename Level::umap_t& map,
    const typename Level::scalar_t& scale,
182
    lo_t dir, lo_t shift,
frey_m's avatar
frey_m committed
183
    Level* mglevel)
184 185 186 187
{
    AmrIntVect_t tmp = iv;
    // first fine cell on refined coarse cell (closer to interface)
    tmp[dir] += shift;
188
    map[mglevel->serialize(tmp)] += 2.0 / 3.0 * scale;
189 190 191
                        
    // second fine cell on refined coarse cell (further away from interface)
    tmp[dir] += shift;
192
    map[mglevel->serialize(tmp)] -= 0.2 * scale;
193 194 195
}


frey_m's avatar
frey_m committed
196 197
template <class Level>
void AmrLagrangeInterpolater<Level>::crseLinear_m(
198
    const AmrIntVect_t& iv,
frey_m's avatar
frey_m committed
199 200
    typename Level::umap_t& map,
    const typename Level::scalar_t& scale,
201
    lo_t dir, lo_t shift, const basefab_t& rfab,
202
    const AmrIntVect_t& riv,
frey_m's avatar
frey_m committed
203
    Level* mglevel)
204 205 206 207 208 209 210 211 212 213 214 215 216
{
#if AMREX_SPACEDIM == 2
    bool top = (riv[(dir+1)%AMREX_SPACEDIM] % 2 == 1);
    
    // right / upper / back
    AmrIntVect_t niv = iv;
    niv[(dir+1)%AMREX_SPACEDIM ] += 1;
    
    // left / lower / front
    AmrIntVect_t miv = iv;
    miv[(dir+1)%AMREX_SPACEDIM ] -= 1;
    
    // 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 222
        map[mglevel->serialize(iv)] += fac * lookup2a_ms[top];
        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 227
        map[mglevel->serialize(iv)] += fac * lookup2a_ms[top];
        map[mglevel->serialize(miv)] += fac * lookup1a_ms[top];
228 229
        
    } else
230 231
        throw OpalException("AmrLagrangeInterpolater::crseLinear_m()",
                            "No valid interpolation scenario found!");
232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254
    
#elif AMREX_SPACEDIM == 3
    
    /*          x   y   z
     * ------------------
     * dir:     0   1   2
     * top1:    1   2   0   (i, L)
     * top2:    2   0   1   (j, K)
     */
    
    /* There are 4 coefficients from Lagrange interpolation.
     * Those are given by the product of one of
     * L0, L1 and one of K0, K1.
     * 
     * g(x, y) = f(x0, y0) * L0(x) * K0(y) +
     *           f(x0, y1) * L0(x) * K1(y) +
     *           f(x1, y0) * L1(x) * K0(y) +
     *           f(x1, y1) * L1(x) * K1(y) +
     */
    
    /*
     * check in 3x3 area (using iv as center) if 4 cells are not covered
     */
255 256
    lo_t d1 = (dir+1)%AMREX_SPACEDIM;
    lo_t d2 = (dir+2)%AMREX_SPACEDIM;
257 258
    
    lbits_t area;
259
    lo_t bit = 0;
260 261 262 263 264 265 266
    
    AmrIntVect_t tmp = iv;
    for (int i = -1; i < 2; ++i) {
        tmp[d1] += i;
        for (int j = -1; j < 2; ++j) {
            
            tmp[d2] += j;
267 268
            // make use of "short-circuit evaluation"
            area[bit] = ( mglevel->isValid(tmp) && rfab(tmp) );
269 270 271 272 273 274 275 276 277 278 279
            ++bit;
            
            // undo
            tmp[d2] -= j;
        }
        // undo
        tmp[d1] -= i;
    }
    
    qpattern_t::const_iterator pit = std::begin(this->lpattern_ms);
    
280
    while ( pit != std::end(this->lpattern_ms) ) {
281 282 283 284 285 286
        if ( *pit == (area & lbits_t(*pit)).to_ulong() )
            break;
        ++pit;
    }
    
    // factor for fine
287
    scalar_t fac = factor_ms * scale;
288
    
289 290
    scalar_t L[2] = {0.0, 0.0};
    lo_t top1 = riv[d1] % 2;
291
    
292 293
    scalar_t K[2] = {0.0, 0.0};
    lo_t top2 = riv[d2] % 2;
294
    
295 296
    lo_t begin[2] = { 0, 0 };
    lo_t end[2]   = { 0, 0 };
297 298 299 300 301
    
    switch ( *pit ) {
        case this->lpattern_ms[0]:
        {
            // corner top right pattern
302 303
            L[0] = lookup1a_ms[top1]; // L_{-1}
            L[1] = lookup2a_ms[top1]; // L_{0}
304 305 306
            begin[0] = -1;
            end[0]   =  0;
            
307 308
            K[0] = lookup1a_ms[top2]; // K_{-1}
            K[1] = lookup2a_ms[top2]; // K_{0}
309 310 311 312 313 314 315
            begin[1] = -1;
            end[1]   =  0;
            break;
        }
        case this->lpattern_ms[1]:
        {
            // corner bottom right pattern
316 317
            L[0] = lookup2b_ms[top1]; // L_{0}
            L[1] = lookup1b_ms[top1]; // L_{1}
318 319 320
            begin[0] = 0;
            end[0]   = 1;
            
321 322
            K[0] = lookup1a_ms[top2]; // K_{-1}
            K[1] = lookup2a_ms[top2]; // K_{0}
323 324 325 326 327 328 329
            begin[1] = -1;
            end[1]   =  0;
            break;
        }
        case this->lpattern_ms[2]:
        {
            // corner bottom left pattern
330 331
            L[0] = lookup2b_ms[top1]; // L_{0}
            L[1] = lookup1b_ms[top1]; // L_{1}
332 333 334
            begin[0] = 0;
            end[0]   = 1;
            
335 336
            K[0] = lookup2b_ms[top2]; // K_{0}
            K[1] = lookup1b_ms[top2]; // K_{1}
337 338 339 340 341 342 343
            begin[1] = 0;
            end[1]   = 1;
            break;
        }
        case this->lpattern_ms[3]:
        {
            // corner top left pattern
344 345
            L[0] = lookup1a_ms[top1]; // L_{-1}
            L[1] = lookup2a_ms[top1]; // L_{0}
346 347 348
            begin[0] = -1;
            end[0]   =  0;
            
349 350
            K[0] = lookup2b_ms[top2]; // K_{0}
            K[1] = lookup1b_ms[top2]; // K_{1}
351 352 353 354 355
            begin[1] = 0;
            end[1]   = 1;
            break;
        }
        default:
356 357
            throw OpalException("AmrLagrangeInterpolater::crseLinear_m()",
                                "No valid interpolation scenario found!");
358 359 360 361 362 363 364 365 366 367 368
    }
    
    /*
     * if pattern is known --> add stencil
     */
    AmrIntVect_t niv = iv;
    for (int i = begin[0]; i <= end[0]; ++i) {
        niv[d1] += i;
        for (int j = begin[1]; j <= end[1]; ++j) {
            niv[d2] += j;
            
369 370 371
            scalar_t value = fac * L[i-begin[0]] * K[j-begin[1]];
            if ( !mglevel->applyBoundary(niv, rfab, map, value) )
                map[mglevel->serialize(niv)] += value;
372 373 374 375 376 377 378 379 380 381 382 383 384 385
            
            // undo
            niv[d2] -= j;
        }
        // undo
        niv[d1] -= i;
    }
#else
    #error Lagrange interpolation: Only 2D and 3D are supported!
#endif
    // the neighbour cancels out
}


frey_m's avatar
frey_m committed
386 387
template <class Level>
void AmrLagrangeInterpolater<Level>::crseQuadratic_m(
388
    const AmrIntVect_t& iv,
frey_m's avatar
frey_m committed
389 390
    typename Level::umap_t& map,
    const typename Level::scalar_t& scale,
391
    lo_t dir, lo_t shift, const basefab_t& rfab,
392
    const AmrIntVect_t& riv,
frey_m's avatar
frey_m committed
393
    Level* mglevel)
394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430
{
#if AMREX_SPACEDIM == 2
    
    bool top = (riv[(dir+1)%AMREX_SPACEDIM] % 2 == 1);
    
    // right / upper / back
    AmrIntVect_t niv = iv;
    niv[(dir+1)%AMREX_SPACEDIM ] += 1;
    
    // left / lower / front
    AmrIntVect_t miv = iv;
    miv[(dir+1)%AMREX_SPACEDIM ] -= 1;
    
    // 2nd right / upper / back
    AmrIntVect_t n2iv = niv;
    n2iv[(dir+1)%AMREX_SPACEDIM ] += 1;
    
    // 2nd left / lower / front
    AmrIntVect_t m2iv = miv;
    m2iv[(dir+1)%AMREX_SPACEDIM ] -= 1;
        
    /* 3 cases:
     * --------
     * r: right neighbour of iv (center of Laplacian)
     * u: upper neighbour of iv
     * b: back neighbour of iv
     * 
     * l: lower / left neighbour of iv
     * f: front neighbour of iv
     * 
     * -1 --> not valid, error happend
     * 0 --> r / u / b and 2nd r / u / b are valid
     * 1 --> direct neighbours (r / u / b and l / f) of iv are valid (has priority)
     * 2 --> l / f and 2nd l / f are valid
     */
    
    // 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 443 444 445 446 447 448 449
    
    if ( rub && lf )
    {
        /*
         * standard case -1, +1 are not-refined nor at physical/mesh boundary
         */            
        // cell is not refined and not at physical boundary
        
        // 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 455
        if ( !mglevel->applyBoundary(niv, rfab, map, value) )
            map[mglevel->serialize(niv)] += value;
456
        
457
        //                      y_t     y_b
458
        value = scale * ((top) ? -0.05 : 1.0 / 12.0);
459 460
        if ( !mglevel->applyBoundary(miv, rfab, map, value) )
            map[mglevel->serialize(miv)] += value;
461 462 463 464 465 466
        
    } else if ( rub && rub2 ) {
        /*
         * corner case --> right / upper / back + 2nd right / upper / back
         */
        //                     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 473
        if ( !mglevel->applyBoundary(niv, rfab, map, value) )
            map[mglevel->serialize(niv)] += value;
474
        
475
        //                      y_t     y_b
476
        value = scale * ((top) ? -0.05 : 1.0 / 12.0);
477 478
        if ( !mglevel->applyBoundary(n2iv, rfab, map, value) )
            map[mglevel->serialize(n2iv)] += value;
479 480 481 482 483
        
    } else if ( lf && lf2 ) {
        /*
         * corner case --> left / lower / front + 2nd left / lower / front
         */
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 491
        if ( !mglevel->applyBoundary(miv, rfab, map, value) )
            map[mglevel->serialize(miv)] += value;
492
        
493
        //                      y_t          y_b
494
        value = scale * ((top) ? 1.0 / 12.0 : -0.05);
495 496
        if ( !mglevel->applyBoundary(m2iv, rfab, map, value) )
            map[mglevel->serialize(m2iv)] += value;
497 498 499 500 501
        
    } else {
        /* last trial: linear Lagrange interpolation
         * --> it throws an error if not possible
         */
502
        this->crseLinear_m(iv, map, scale, dir, shift, rfab, riv, mglevel);
503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532
    }
    
#elif AMREX_SPACEDIM == 3
    
    /*          x   y   z
     * ------------------
     * dir:     0   1   2
     * top1:    1   2   0   (i, L)
     * top2:    2   0   1   (j, K)
     */
    
    /* There are 9 coefficients from Lagrange interpolation.
     * Those are given by the product of one of
     * L0, L1, L2 and one of K0, K1, K2.
     * 
     * g(x, y) = f(x0, y0) * L0(x) * K0(y) +
     *           f(x0, y1) * L0(x) * K1(y) +
     *           f(x0, y2) * L0(x) * K2(y) +
     *           f(x1, y0) * L1(x) * K0(y) +
     *           f(x1, y1) * L1(x) * K1(y) +
     *           f(x1, y2) * L1(x) * K2(y) +
     *           f(x2, y0) * L2(x) * K0(y) +
     *           f(x2, y1) * L2(x) * K1(y) +
     *           f(x2, y2) * L2(x) * K2(y) +
     */
    
    
    /*
     * check in 5x5 area (using iv as center) if 9 cells are not covered
     */
533 534
    lo_t d1 = (dir+1)%AMREX_SPACEDIM;
    lo_t d2 = (dir+2)%AMREX_SPACEDIM;
535 536
    
    qbits_t area;
537
    lo_t bit = 0;
538 539 540 541 542 543 544 545
    
    AmrIntVect_t tmp = iv;
    for (int i = -2; i < 3; ++i) {
        tmp[d1] += i;
        for (int j = -2; j < 3; ++j) {
            
            tmp[d2] += j;
            
546
            area[bit] = ( mglevel->isValid(tmp) && rfab(tmp) );
547 548 549 550 551 552 553 554 555 556
            ++bit;
            
            // undo
            tmp[d2] -= j;
        }
        // undo
        tmp[d1] -= i;
    }
    
    qpattern_t::const_iterator pit = std::begin(this->qpattern_ms);
557 558

    while ( pit != std::end(this->qpattern_ms) ) {
559 560 561 562 563 564
        if ( *pit == (area & qbits_t(*pit)).to_ulong() )
            break;
        ++pit;
    }
    
    // factor for fine
565
    scalar_t fac = factor_ms * scale;
566
    
567 568
    scalar_t L[3] = {0.0, 0.0, 0.0};
    lo_t top1 = riv[d1] % 2;
569
    
570 571
    scalar_t K[3] = {0.0, 0.0, 0.0};
    lo_t top2 = riv[d2] % 2;
572
    
573 574
    lo_t begin[2] = { 0, 0 };
    lo_t end[2]   = { 0, 0 };
575 576 577 578 579
    
    switch ( *pit ) {
        case this->qpattern_ms[0]:
        {
            // cross pattern
580 581 582
            L[0] = lookup3_ms[top1];  // L_{-1}
            L[1] = lookup6_ms;        // L_{0}
            L[2] = lookup3r_ms[top1]; // L_{1}
583 584 585
            begin[0] = -1;
            end[0]   =  1;
            
586 587 588
            K[0] = lookup3_ms[top2];  // K_{-1}
            K[1] = lookup6_ms;        // K_{0}
            K[2] = lookup3r_ms[top2]; // K_{1}
589 590 591 592 593 594 595
            begin[1] = -1;
            end[1]   =  1;
            break;
        }
        case this->qpattern_ms[1]:
        {
            // T pattern
596 597 598
            L[0] = lookup3r_ms[top1]; // L_{-2}
            L[1] = lookup4_ms[top1];  // L_{-1}
            L[2] = lookup5r_ms[top1]; // L_{0}
599 600 601
            begin[0] = -2;
            end[0]   =  0;
            
602 603 604
            K[0] = lookup3_ms[top2];  // K_{-1}
            K[1] = lookup6_ms;        // K_{0}
            K[2] = lookup3r_ms[top2]; // K_{1}
605 606 607 608 609 610 611
            begin[1] = -1;
            end[1]   =  1;
            break;
        }
        case this->qpattern_ms[2]:
        {
            // right hammer pattern
612 613 614
            L[0] = lookup3_ms[top1];  // L_{-1}
            L[1] = lookup6_ms;        // L_{0}
            L[2] = lookup3r_ms[top1]; // L_{1}
615 616 617
            begin[0] = -1;
            end[0] = 1;
            
618 619 620
            K[0] = lookup3r_ms[top2]; // K_{-2}
            K[1] = lookup4_ms[top2];  // K_{-1}
            K[2] = lookup5r_ms[top2]; // K_{0}
621 622 623 624 625 626 627
            begin[1] = -2;
            end[1]   = 0;
            break;
        }
        case this->qpattern_ms[3]:
        {
            // T on head pattern
628 629 630
            L[0] = lookup5_ms[top1];  // L_{0}
            L[1] = lookup4r_ms[top1]; // L_{1}
            L[2] = lookup3_ms[top1];  // L_{2}
631 632 633
            begin[0] = 0;
            end[0]   = 2;
            
634 635 636
            K[0] = lookup3_ms[top2];  // K_{-1}
            K[1] = lookup6_ms;        // K_{0}
            K[2] = lookup3r_ms[top2]; // K_{1}
637 638 639 640 641 642 643
            begin[1] = -1;
            end[1]   =  1;
            break;
        }
        case this->qpattern_ms[4]:
        {
            // left hammer pattern
644 645 646
            L[0] = lookup3_ms[top1];  // L_{-1}
            L[1] = lookup6_ms;        // L_{0}
            L[2] = lookup3r_ms[top1]; // L_{1}
647 648 649
            begin[0] = -1;
            end[0] = 1;
            
650 651 652
            K[0] = lookup5_ms[top2];  // K_{0}
            K[1] = lookup4r_ms[top2]; // K_{1}
            K[2] = lookup3_ms[top2];  // K_{2}
653 654 655 656 657 658 659
            begin[1] = 0;
            end[1]   = 2;
            break;
        }
        case this->qpattern_ms[5]:
        {
            // upper left corner pattern
660 661 662
            L[0] = lookup3r_ms[top1]; // L_{-2}
            L[1] = lookup4_ms[top1];  // L_{-1}
            L[2] = lookup5r_ms[top1]; // L_{0}
663 664 665
            begin[0] = -2;
            end[0]   =  0;
            
666 667 668
            K[0] = lookup5_ms[top2];  // K_{0}
            K[1] = lookup4r_ms[top2]; // K_{1}
            K[2] = lookup3_ms[top2];  // K_{2}
669 670 671 672 673 674 675
            begin[1] = 0;
            end[1]   = 2;
            break;
        }
        case this->qpattern_ms[6]:
        {
            // upper right corner pattern
676 677 678
            L[0] = lookup3r_ms[top1]; // L_{-2}
            L[1] = lookup4_ms[top1];  // L_{-1}
            L[2] = lookup5r_ms[top1]; // L_{0}
679 680 681
            begin[0] = -2;
            end[0]   =  0;
            
682 683 684
            K[0] = lookup3r_ms[top2]; // K_{-2}
            K[1] = lookup4_ms[top2];  // K_{-1}
            K[2] = lookup5r_ms[top2]; // K_{0}
685 686 687 688 689 690 691
            begin[1] = -2;
            end[1]   =  0;
            break;
        }
        case this->qpattern_ms[7]:
        {
            // mirrored L pattern
692 693 694
            L[0] = lookup5_ms[top1];  // L_{0}
            L[1] = lookup4r_ms[top1]; // L_{1}
            L[2] = lookup3_ms[top1];  // L_{2}
695 696 697
            begin[0] = 0;
            end[0]   = 2;
            
698 699 700
            K[0] = lookup3r_ms[top2]; // K_{-2}
            K[1] = lookup4_ms[top2];  // K_{-1}
            K[2] = lookup5r_ms[top2]; // K_{0}
701 702 703 704 705 706 707
            begin[1] = -2;
            end[1]   =  0;
            break;
        }
        case this->qpattern_ms[8]:
        {
            // L pattern
708 709 710
            L[0] = lookup5_ms[top1];  // L_{0}
            L[1] = lookup4r_ms[top1]; // L_{1}
            L[2] = lookup3_ms[top1];  // L_{2}
711 712 713
            begin[0] = 0;
            end[0]   = 2;
            
714 715 716
            K[0] = lookup5_ms[top2];  // K_{0}
            K[1] = lookup4r_ms[top2]; // K_{1}
            K[2] = lookup3_ms[top2];  // K_{2}
717 718 719 720 721 722 723 724 725
            begin[1] = 0;
            end[1]   = 2;
            break;
        }
        default:
        {
            /* unknown pattern --> last trial: linear Lagrange interpolation
             * --> it throws an error if not possible
             */
726 727
            this->crseLinear_m(iv, map, scale, dir, shift, rfab, riv, mglevel);
            return;
728 729 730
        }
    }
    
731 732 733 734 735 736 737 738 739 740 741 742 743
    /*
     * if pattern is known --> add stencil
     */
    AmrIntVect_t tmp1 = iv;
    for (int i = begin[0]; i <= end[0]; ++i) {
        tmp1[d1] += i;
        for (int j = begin[1]; j <= end[1]; ++j) {
            tmp1[d2] += j;
            
            scalar_t value = fac * L[i-begin[0]] * K[j-begin[1]];
            if ( !mglevel->applyBoundary(tmp1, rfab, map, value) )
                map[mglevel->serialize(tmp1)] += value;
            
744
            // undo
745
            tmp1[d2] -= j;
746
        }
747 748
        // undo
        tmp1[d1] -= i;
749 750 751 752 753 754
    }
    
#else
    #error Lagrange interpolation: Only 2D and 3D are supported!
#endif
}