OPAL (Object Oriented Parallel Accelerator Library)  2021.1.99
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
24 template <class Level>
25 constexpr typename AmrLagrangeInterpolater<Level>::qpattern_t
26  AmrLagrangeInterpolater<Level>::qpattern_ms;
27 
28 
29 template <class Level>
30 constexpr typename AmrLagrangeInterpolater<Level>::lpattern_t
31  AmrLagrangeInterpolater<Level>::lpattern_ms;
32 #endif
33 
34 // y_b y_t
35 template <class Level>
36 const typename AmrLagrangeInterpolater<Level>::scalar_t
37  AmrLagrangeInterpolater<Level>::lookup1a_ms[2] = {0.25, -0.25};
38 
39 template <class Level>
40 const typename AmrLagrangeInterpolater<Level>::scalar_t
41  AmrLagrangeInterpolater<Level>::lookup2a_ms[2] = {0.75, 1.25};
42 
43 template <class Level>
44 const typename AmrLagrangeInterpolater<Level>::scalar_t
45  AmrLagrangeInterpolater<Level>::lookup1b_ms[2] = {-0.25, 0.25};
46 
47 template <class Level>
48 const typename AmrLagrangeInterpolater<Level>::scalar_t
49  AmrLagrangeInterpolater<Level>::lookup2b_ms[2] = {1.25, 0.75};
50 
51 #if AMREX_SPACEDIM == 3
52 template <class Level>
53 const typename AmrLagrangeInterpolater<Level>::scalar_t
54  AmrLagrangeInterpolater<Level>::lookup3_ms[2] = {5.0 / 32.0, -3.0 / 32.0 };
55 
56 template <class Level>
57 const typename AmrLagrangeInterpolater<Level>::scalar_t
58  AmrLagrangeInterpolater<Level>::lookup3r_ms[2] = {-3.0 / 32.0, 5.0 / 32.0 };
59 
60 template <class Level>
61 const typename AmrLagrangeInterpolater<Level>::scalar_t
62  AmrLagrangeInterpolater<Level>::lookup4_ms[2] = {7.0 / 16.0, -9.0 / 16.0 };
63 
64 template <class Level>
65 const typename AmrLagrangeInterpolater<Level>::scalar_t
66  AmrLagrangeInterpolater<Level>::lookup4r_ms[2] = {-9.0 / 16.0, 7.0 / 16.0 };
67 
68 template <class Level>
69 const typename AmrLagrangeInterpolater<Level>::scalar_t
70  AmrLagrangeInterpolater<Level>::lookup5_ms[2] = {45.0 / 32.0, 21.0 / 32.0};
71 
72 template <class Level>
73 const typename AmrLagrangeInterpolater<Level>::scalar_t
74  AmrLagrangeInterpolater<Level>::lookup5r_ms[2] = {21.0 / 32.0, 45.0 / 32.0};
75 
76 template <class Level>
77 const typename AmrLagrangeInterpolater<Level>::scalar_t
78  AmrLagrangeInterpolater<Level>::lookup6_ms = 15.0 / 16.0;
79 
80 template <class Level>
81 const typename AmrLagrangeInterpolater<Level>::scalar_t
82  AmrLagrangeInterpolater<Level>::factor_ms = 8.0 / 15.0;
83 #endif
84 
85 
86 
87 
88 template <class Level>
89 AmrLagrangeInterpolater<Level>::AmrLagrangeInterpolater(Order order)
90  : AmrInterpolater<Level>( lo_t(order) + 1 )
91 { }
92 
93 
94 template <class Level>
95 void 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 
106 template <class Level>
107 void 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 
131 template <class Level>
132 void 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 
155 template <class Level>
156 void 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 
177 template <class Level>
178 void 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 
196 template <class Level>
197 void 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 
386 template <class Level>
387 void 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:735
begin
PartBunchBase< T, Dim >::ConstIterator begin(PartBunchBase< T, Dim > const &bunch)
Definition: PartBunchBase.h:730
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 Wed Aug 25 2021 16:41:04 for OPAL (Object Oriented Parallel Accelerator Library) by doxygen 1.9.1