Alamo
Stencil.H
Go to the documentation of this file.
1 #ifndef NUMERIC_STENCIL_H_
2 #define NUMERIC_STENCIL_H_
3 
4 #include <AMReX.H>
5 #include <AMReX_MultiFab.H>
6 #include "Set/Set.H"
7 
8 /// \brief This namespace contains some numerical tools
9 namespace Numeric
10 {
11 
12 enum StencilType { Lo, Hi, Central };
13 static std::array<StencilType, AMREX_SPACEDIM>
14 DefaultType = { AMREX_D_DECL(StencilType::Central, StencilType::Central, StencilType::Central) };
15 
16 [[nodiscard]]
17 static
18 AMREX_FORCE_INLINE
19 std::array<StencilType, AMREX_SPACEDIM>
20 GetStencil(const int i, const int j, const int k, const amrex::Box domain)
21 {
22  (void)i; (void)j; (void)k; // Suppress "unused variable" warnings.
23  std::array<StencilType, AMREX_SPACEDIM> sten;
24  const amrex::Dim3 lo = amrex::lbound(domain), hi = amrex::ubound(domain);
25  AMREX_D_TERM(sten[0] = (i == lo.x ? Numeric::StencilType::Hi :
26  i == hi.x ? Numeric::StencilType::Lo :
27  Numeric::StencilType::Central);,
28  sten[1] = (j == lo.y ? Numeric::StencilType::Hi :
29  j == hi.y ? Numeric::StencilType::Lo :
30  Numeric::StencilType::Central);,
31  sten[2] = (k == lo.z ? Numeric::StencilType::Hi :
32  k == hi.z ? Numeric::StencilType::Lo :
33  Numeric::StencilType::Central););
34  return sten;
35 }
36 
37 template<class T, int x, int y, int z>
38 struct Stencil
39 {};
40 
41 //
42 // FIRST order derivatives
43 //
44 
45 template<class T>
46 struct Stencil<T, 1, 0, 0>
47 {
48  [[nodiscard]] AMREX_FORCE_INLINE
49  static T D(const amrex::Array4<const T>& f,
50  const int& i, const int& j, const int& k, const int& m,
51  const Set::Scalar dx[AMREX_SPACEDIM],
52  std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
53  {
54  if (stencil[0] == StencilType::Lo)
55  return (f(i, j, k, m) - f(i - 1, j, k, m)) / dx[0]; // 1st order stencil
56  else if (stencil[0] == StencilType::Hi)
57  return (f(i + 1, j, k, m) - f(i, j, k, m)) / dx[0]; // 1st order stencil
58  else
59  return (f(i + 1, j, k, m) - f(i - 1, j, k, m)) * 0.5 / dx[0];
60  };
61  [[nodiscard]] AMREX_FORCE_INLINE
62  static std::pair<Set::Scalar,T>
63  Dsplit( const amrex::Array4<const T>& f,
64  const int& i, const int& j, const int& k, const int& m,
65  const Set::Scalar dx[AMREX_SPACEDIM],
66  std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
67  {
68  if (stencil[0] == StencilType::Lo)
69  return std::make_pair(
70  1.0 / dx[0],
71  -f(i - 1, j, k, m) / dx[0]
72  ); // 1st order
73  else if (stencil[0] == StencilType::Hi)
74  return std::make_pair(
75  -1.0 / dx[0],
76  f(i + 1, j, k, m) / dx[0]
77  ); // 1st order
78  else
79  return std::make_pair(
80  0.0,
81  (f(i + 1, j, k, m) - f(i - 1, j, k, m)) * 0.5 / dx[0]
82  );
83  };
84 };
85 
86 template<class T>
87 struct Stencil<T, 0, 1, 0>
88 {
89  [[nodiscard]] AMREX_FORCE_INLINE
90  static T D(const amrex::Array4<const T>& f,
91  const int& i, const int& j, const int& k, const int& m,
92  const Set::Scalar dx[AMREX_SPACEDIM],
93  std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
94  {
95  if (stencil[1] == StencilType::Lo)
96  return (f(i, j, k, m) - f(i, j - 1, k, m)) / dx[1];
97  else if (stencil[1] == StencilType::Hi)
98  return (f(i, j + 1, k, m) - f(i, j, k, m)) / dx[1];
99  else
100  return (f(i, j + 1, k, m) - f(i, j - 1, k, m)) * 0.5 / dx[1];
101  };
102  [[nodiscard]] AMREX_FORCE_INLINE
103  static std::pair<Set::Scalar,T>
104  Dsplit(const amrex::Array4<const T>& f,
105  const int& i, const int& j, const int& k, const int& m,
106  const Set::Scalar dx[AMREX_SPACEDIM],
107  std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
108  {
109  if (stencil[1] == StencilType::Lo)
110  return std::make_pair(
111  1.0 / dx[1],
112  - f(i, j - 1, k, m) / dx[1]
113  );
114  else if (stencil[1] == StencilType::Hi)
115  return std::make_pair(
116  - 1.0 / dx[1],
117  f(i, j + 1, k, m) / dx[1]
118  );
119  else
120  return std::make_pair(
121  0.0,
122  (f(i, j + 1, k, m) - f(i, j - 1, k, m)) * 0.5 / dx[1]
123  );
124  };
125 };
126 
127 template<class T>
128 struct Stencil<T, 0, 0, 1>
129 {
130  [[nodiscard]] AMREX_FORCE_INLINE
131  static T D(const amrex::Array4<const T>& f,
132  const int& i, const int& j, const int& k, const int& m,
133  const Set::Scalar dx[AMREX_SPACEDIM],
134  std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
135  {
136  if (stencil[2] == StencilType::Lo)
137  return (f(i, j, k, m) - f(i, j, k - 1, m)) / dx[2];
138  else if (stencil[2] == StencilType::Hi)
139  return (f(i, j, k + 1, m) - f(i, j, k, m)) / dx[2];
140  else
141  return (f(i, j, k + 1, m) - f(i, j, k - 1, m)) * 0.5 / dx[2];
142  };
143  [[nodiscard]] AMREX_FORCE_INLINE
144  static std::pair<Set::Scalar, T>
145  Dsplit(const amrex::Array4<const T>& f,
146  const int& i, const int& j, const int& k, const int& m,
147  const Set::Scalar dx[AMREX_SPACEDIM],
148  std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
149  {
150  if (stencil[2] == StencilType::Lo)
151  return std::make_pair(
152  1.0 / dx[2],
153  -f(i, j, k - 1, m) / dx[2]
154  );
155  else if (stencil[2] == StencilType::Hi)
156  return std::make_pair(
157  -1.0 / dx[2],
158  f(i, j, k + 1, m) / dx[2]
159  );
160  else
161  return std::make_pair(
162  0.0,
163  (f(i, j, k + 1, m) - f(i, j, k - 1, m)) * 0.5 / dx[2]
164  );
165  };
166 };
167 
168 //
169 // SECOND order derivatives
170 //
171 
172 template<class T>
173 struct Stencil<T, 2, 0, 0>
174 {
175  [[nodiscard]] AMREX_FORCE_INLINE
176  static T D(const amrex::Array4<const T>& f,
177  const int& i, const int& j, const int& k, const int& m,
178  const Set::Scalar dx[AMREX_SPACEDIM],
179  std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
180  {
181  if (stencil[0] == StencilType::Central)
182  return (f(i + 1, j, k, m) - 2.0 * f(i, j, k, m) + f(i - 1, j, k, m)) / dx[0] / dx[0];
183  else
184  return 0.0 * f(i, j, k, m); // TODO this is not a great way to do this
185  };
186 };
187 
188 template<class T>
189 struct Stencil<T, 0, 2, 0>
190 {
191  [[nodiscard]] AMREX_FORCE_INLINE
192  static T D(const amrex::Array4<const T>& f,
193  const int& i, const int& j, const int& k, const int& m,
194  const Set::Scalar dx[AMREX_SPACEDIM],
195  std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
196  {
197  if (stencil[1] == StencilType::Central)
198  return (f(i, j + 1, k, m) - 2.0 * f(i, j, k, m) + f(i, j - 1, k, m)) / dx[1] / dx[1];
199  else
200  return 0.0 * f(i, j, k, m); // TODO this is not a great way to do this
201  };
202 };
203 
204 template<class T>
205 struct Stencil<T, 0, 0, 2>
206 {
207  [[nodiscard]] AMREX_FORCE_INLINE
208  static T D(const amrex::Array4<const T>& f,
209  const int& i, const int& j, const int& k, const int& m,
210  const Set::Scalar dx[AMREX_SPACEDIM],
211  std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
212  {
213  if (stencil[2] == StencilType::Central)
214  return (f(i, j, k + 1, m) - 2.0 * f(i, j, k, m) + f(i, j, k - 1, m)) / dx[2] / dx[2];
215  else
216  return 0.0 * f(i, j, k, m); // TODO this is not a great way to do this
217  };
218 };
219 
220 template<class T>
221 struct Stencil<T, 1, 1, 0>
222 {
223  [[nodiscard]] AMREX_FORCE_INLINE
224  static T D(const amrex::Array4<const T>& f,
225  const int& i, const int& j, const int& k, const int& m,
226  const Set::Scalar dx[AMREX_SPACEDIM],
227  std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
228  {
229  int ihi = 1, ilo = 1, jhi = 1, jlo = 1;
230  Set::Scalar ifac = 0.5, jfac = 0.5;
231  if (stencil[0] == StencilType::Hi) { ilo = 0; ifac = 1.0; }
232  if (stencil[0] == StencilType::Lo) { ihi = 0; ifac = 1.0; }
233  if (stencil[1] == StencilType::Hi) { jlo = 0; jfac = 1.0; }
234  if (stencil[1] == StencilType::Lo) { jhi = 0; jfac = 1.0; }
235 
236  return ifac * jfac * (f(i + ihi, j + jhi, k, m) + f(i - ilo, j - jlo, k, m) - f(i + ihi, j - jlo, k, m) - f(i - ilo, j + jhi, k, m)) / (dx[0] * dx[1]);
237  };
238 };
239 template<class T>
240 struct Stencil<T, 1, 0, 1>
241 {
242  [[nodiscard]] AMREX_FORCE_INLINE
243  static T D(const amrex::Array4<const T>& f,
244  const int& i, const int& j, const int& k, const int& m,
245  const Set::Scalar dx[AMREX_SPACEDIM],
246  std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
247  {
248  int khi = 1, klo = 1, ihi = 1, ilo = 1;
249  Set::Scalar kfac = 0.5, ifac = 0.5;
250  if (stencil[0] == StencilType::Hi) { ilo = 0; ifac = 1.0; }
251  if (stencil[0] == StencilType::Lo) { ihi = 0; ifac = 1.0; }
252  if (stencil[2] == StencilType::Hi) { klo = 0; kfac = 1.0; }
253  if (stencil[2] == StencilType::Lo) { khi = 0; kfac = 1.0; }
254 
255  return kfac * ifac * (f(i + ihi, j, k + khi, m) + f(i - ilo, j, k - klo, m) - f(i + ihi, j, k - klo, m) - f(i - ilo, j, k + khi, m)) / (dx[0] * dx[2]);
256  };
257 };
258 template<class T>
259 struct Stencil<T, 0, 1, 1>
260 {
261  [[nodiscard]] AMREX_FORCE_INLINE
262  static T D(const amrex::Array4<const T>& f,
263  const int& i, const int& j, const int& k, const int& m,
264  const Set::Scalar dx[AMREX_SPACEDIM],
265  std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
266  {
267  int jhi = 1, jlo = 1, khi = 1, klo = 1;
268  Set::Scalar jfac = 0.5, kfac = 0.5;
269  if (stencil[1] == StencilType::Hi) { jlo = 0; jfac = 1.0; }
270  if (stencil[1] == StencilType::Lo) { jhi = 0; jfac = 1.0; }
271  if (stencil[2] == StencilType::Hi) { klo = 0; kfac = 1.0; }
272  if (stencil[2] == StencilType::Lo) { khi = 0; kfac = 1.0; }
273 
274  return jfac * kfac * (f(i, j + jhi, k + khi, m) + f(i, j - jlo, k - klo, m) - f(i, j + jhi, k - klo, m) - f(i, j - jlo, k + khi, m)) / (dx[1] * dx[2]);
275  };
276 };
277 
278 //
279 // FOURTH order derivatives
280 //
281 
282 template<class T>
283 struct Stencil<T, 4, 0, 0>
284 {
285  [[nodiscard]] AMREX_FORCE_INLINE
286  static T D(const amrex::Array4<const T>& f,
287  const int& i, const int& j, const int& k, const int& m,
288  const Set::Scalar dx[AMREX_SPACEDIM])
289  {
290  return ((f(i + 2, j, k, m)) - 4. * (f(i + 1, j, k, m)) + 6. * (f(i, j, k, m)) - 4. * (f(i - 1, j, k, m)) + (f(i - 2, j, k, m))) /
291  (dx[0] * dx[0] * dx[0] * dx[0]);
292  };
293 };
294 template<class T>
295 struct Stencil<T, 0, 4, 0>
296 {
297  [[nodiscard]] AMREX_FORCE_INLINE
298  static T D(const amrex::Array4<const T>& f,
299  const int& i, const int& j, const int& k, const int& m,
300  const Set::Scalar dx[AMREX_SPACEDIM])
301  {
302  return ((f(i, j + 2, k, m)) - 4. * (f(i, j + 1, k, m)) + 6. * (f(i, j, k, m)) - 4. * (f(i, j - 1, k, m)) + (f(i, j - 2, k, m))) /
303  (dx[1] * dx[1] * dx[1] * dx[1]);
304  };
305 };
306 template<class T>
307 struct Stencil<T, 0, 0, 4>
308 {
309  [[nodiscard]] AMREX_FORCE_INLINE
310  static T D(const amrex::Array4<const T>& f,
311  const int& i, const int& j, const int& k, const int& m,
312  const Set::Scalar dx[AMREX_SPACEDIM])
313  {
314  return ((f(i, j, k + 2, m)) - 4. * (f(i, j, k + 1, m)) + 6. * (f(i, j, k, m)) - 4. * (f(i, j, k - 1, m)) + (f(i, j, k - 2, m))) /
315  (dx[2] * dx[2] * dx[2] * dx[2]);
316  };
317 };
318 
319 /// Compute
320 /// \f[ \frac{\partial^4}{\partial^3 x_1\partial x_2\f]
321 template<class T>
322 struct Stencil<T, 3, 1, 0>
323 {
324  [[nodiscard]] AMREX_FORCE_INLINE
325  static T D(const amrex::Array4<const T>& f,
326  const int& i, const int& j, const int& k, const int& m,
327  const Set::Scalar dx[AMREX_SPACEDIM])
328  {
329  return ((-f(i + 2, j + 2, k, m) + 8.0 * f(i + 2, j + 1, k, m) - 8.0 * f(i + 2, j - 1, k, m) + f(i + 2, j - 2, k, m))
330  - 2 * (-f(i + 1, j + 2, k, m) + 8.0 * f(i + 1, j + 1, k, m) - 8.0 * f(i + 1, j - 1, k, m) + f(i + 1, j - 2, k, m))
331  + 2 * (-f(i - 1, j + 2, k, m) + 8.0 * f(i - 1, j + 1, k, m) - 8.0 * f(i - 1, j - 1, k, m) + f(i - 1, j - 2, k, m))
332  - (-f(i - 2, j + 2, k, m) + 8.0 * f(i - 2, j + 1, k, m) - 8.0 * f(i - 2, j - 1, k, m) + f(i - 2, j - 2, k, m))) /
333  (24.0 * dx[0] * dx[0] * dx[0] * dx[1]);
334  };
335 };
336 
337 /// Compute
338 /// \f[ \frac{\partial^4}{\partial^3 x_2\partial x_1\f]
339 template<class T>
340 struct Stencil<T, 1, 3, 0>
341 {
342  [[nodiscard]] AMREX_FORCE_INLINE
343  static T D(const amrex::Array4<const T>& f,
344  const int& i, const int& j, const int& k, const int& m,
345  const Set::Scalar dx[AMREX_SPACEDIM])
346  {
347  return ((-f(i + 2, j + 2, k, m) + 8.0 * f(i + 1, j + 2, k, m) - 8.0 * f(i - 1, j + 2, k, m) + f(i - 2, j + 2, k, m))
348  - 2 * (-f(i + 2, j + 1, k, m) + 8.0 * f(i + 1, j + 1, k, m) - 8.0 * f(i - 1, j + 1, k, m) + f(i - 2, j + 1, k, m))
349  + 2 * (-f(i + 2, j - 1, k, m) + 8.0 * f(i + 1, j - 1, k, m) - 8.0 * f(i - 1, j - 1, k, m) + f(i - 2, j - 1, k, m))
350  - (-f(i + 2, j - 2, k, m) + 8.0 * f(i + 1, j - 2, k, m) - 8.0 * f(i - 1, j - 2, k, m) + f(i - 2, j - 2, k, m))) /
351  (24.0 * dx[0] * dx[1] * dx[1] * dx[1]);
352  };
353 };
354 
355 /// Compute
356 /// \f[ \frac{\partial^4}{\partial^3 x_2\partial x_3\f]
357 template<class T>
358 struct Stencil<T, 0, 3, 1>
359 {
360  [[nodiscard]] AMREX_FORCE_INLINE
361  static T D(const amrex::Array4<const T>& f,
362  const int& i, const int& j, const int& k, const int& m,
363  const Set::Scalar dx[AMREX_SPACEDIM])
364  {
365  return ((-f(i, j + 2, k + 2, m) + 8.0 * f(i, j + 2, k + 1, m) - 8.0 * f(i, j + 2, k - 1, m) + f(i, j + 2, k - 2, m))
366  - 2 * (-f(i, j + 1, k + 2, m) + 8.0 * f(i, j + 1, k + 1, m) - 8.0 * f(i, j + 1, k - 1, m) + f(i, j + 1, k - 2, m))
367  + 2 * (-f(i, j - 1, k + 2, m) + 8.0 * f(i, j - 1, k + 1, m) - 8.0 * f(i, j - 1, k - 1, m) + f(i, j - 1, k - 2, m))
368  - (-f(i, j - 2, k + 2, m) + 8.0 * f(i, j - 2, k + 1, m) - 8.0 * f(i, j - 2, k - 1, m) + f(i, j - 2, k - 2, m))) /
369  (24.0 * dx[1] * dx[1] * dx[1] * dx[2]);
370  };
371 };
372 /// \brief Compute \f[ \frac{\partial^4}{\partial^3 x_3\partial x_2}\f]
373 template<class T>
374 struct Stencil<T, 0, 1, 3>
375 {
376  [[nodiscard]] AMREX_FORCE_INLINE
377  static T D(const amrex::Array4<const T>& f,
378  const int& i, const int& j, const int& k, const int& m,
379  const Set::Scalar dx[AMREX_SPACEDIM])
380  {
381  return ((-f(i, j + 2, k + 2, m) + 8.0 * f(i, j + 1, k + 2, m) - 8.0 * f(i, j - 1, k + 2, m) + f(i, j - 2, k + 2, m))
382  - 2 * (-f(i, j + 2, k + 1, m) + 8.0 * f(i, j + 1, k + 1, m) - 8.0 * f(i, j - 1, k + 1, m) + f(i, j - 2, k + 1, m))
383  + 2 * (-f(i, j + 2, k - 1, m) + 8.0 * f(i, j + 1, k - 1, m) - 8.0 * f(i, j - 1, k - 1, m) + f(i, j - 2, k - 1, m))
384  - (-f(i, j + 2, k - 2, m) + 8.0 * f(i, j + 1, k - 2, m) - 8.0 * f(i, j - 1, k - 2, m) + f(i, j - 2, k - 2, m))) /
385  (24.0 * dx[1] * dx[2] * dx[2] * dx[2]);
386  };
387 };
388 /// \brief Compute \f[ \frac{\partial^4}{\partial^3 x_3\partial x_1}\f]
389 template<class T>
390 struct Stencil<T, 1, 0, 3>
391 {
392  [[nodiscard]] AMREX_FORCE_INLINE
393  static T D(const amrex::Array4<const T>& f,
394  const int& i, const int& j, const int& k, const int& m,
395  const Set::Scalar dx[AMREX_SPACEDIM])
396  {
397  return ((-f(i + 2, j, k + 2, m) + 8.0 * f(i + 1, j, k + 2, m) - 8.0 * f(i - 1, j, k + 2, m) + f(i - 2, j, k + 2, m))
398  - 2 * (-f(i + 2, j, k + 1, m) + 8.0 * f(i + 1, j, k + 1, m) - 8.0 * f(i - 1, j, k + 1, m) + f(i - 2, j, k + 1, m))
399  + 2 * (-f(i + 2, j, k - 1, m) + 8.0 * f(i + 1, j, k - 1, m) - 8.0 * f(i - 1, j, k - 1, m) + f(i - 2, j, k - 1, m))
400  - (-f(i + 2, j, k - 2, m) + 8.0 * f(i + 1, j, k - 2, m) - 8.0 * f(i - 1, j, k - 2, m) + f(i - 2, j, k - 2, m))) /
401  (24.0 * dx[0] * dx[2] * dx[2] * dx[2]);
402 
403  };
404 };
405 /// \brief Compute \f[ \frac{\partial^4}{\partial^3 x_1\partial x_3}\f]
406 template<class T>
407 struct Stencil<T, 3, 0, 1>
408 {
409  [[nodiscard]] AMREX_FORCE_INLINE
410  static T D(const amrex::Array4<const T>& f,
411  const int& i, const int& j, const int& k, const int& m,
412  const Set::Scalar dx[AMREX_SPACEDIM])
413  {
414  return ((-f(i + 2, j, k + 2, m) + 8.0 * f(i + 2, j, k + 1, m) - 8.0 * f(i + 2, j, k - 1, m) + f(i + 2, j, k - 2, m))
415  - 2 * (-f(i + 1, j, k + 2, m) + 8.0 * f(i + 1, j, k + 1, m) - 8.0 * f(i + 1, j, k - 1, m) + f(i + 1, j, k - 2, m))
416  + 2 * (-f(i - 1, j, k + 2, m) + 8.0 * f(i - 1, j, k + 1, m) - 8.0 * f(i - 1, j, k - 1, m) + f(i - 1, j, k - 2, m))
417  - (-f(i - 2, j, k + 2, m) + 8.0 * f(i - 2, j, k + 1, m) - 8.0 * f(i - 2, j, k - 1, m) + f(i - 2, j, k - 2, m))) /
418  (24.0 * dx[0] * dx[0] * dx[0] * dx[2]);
419 
420  };
421 };
422 /// \brief Compute \f[ \frac{\partial^4}{\partial^2 x_1\partial^2 x_2}\f]
423 template<class T>
424 struct Stencil<T, 2, 2, 0>
425 {
426  [[nodiscard]] AMREX_FORCE_INLINE
427  static T D(const amrex::Array4<const T>& f,
428  const int& i, const int& j, const int& k, const int& m,
429  const Set::Scalar dx[AMREX_SPACEDIM])
430  {
431  return (-(-f(i + 2, j + 2, k, m) + 16.0 * f(i + 1, j + 2, k, m) - 30.0 * f(i, j + 2, k, m) + 16.0 * f(i - 1, j + 2, k, m) - f(i - 2, j + 2, k, m))
432  + 16 * (-f(i + 2, j + 1, k, m) + 16.0 * f(i + 1, j + 1, k, m) - 30.0 * f(i, j + 1, k, m) + 16.0 * f(i - 1, j + 1, k, m) - f(i - 2, j + 1, k, m))
433  - 30 * (-f(i + 2, j, k, m) + 16.0 * f(i + 1, j, k, m) - 30.0 * f(i, j, k, m) + 16.0 * f(i - 1, j, k, m) - f(i - 2, j, k, m))
434  + 16 * (-f(i + 2, j - 1, k, m) + 16.0 * f(i + 1, j - 1, k, m) - 30.0 * f(i, j - 1, k, m) + 16.0 * f(i - 1, j - 1, k, m) - f(i - 2, j - 1, k, m))
435  - (-f(i + 2, j - 2, k, m) + 16.0 * f(i + 1, j - 2, k, m) - 30.0 * f(i, j - 2, k, m) + 16.0 * f(i - 1, j - 2, k, m) - f(i - 2, j - 2, k, m))) /
436  (144.0 * dx[0] * dx[0] * dx[1] * dx[1]);
437  };
438 };
439 
440 /// \brief Compute \f[ \frac{\partial^4}{\partial^2 x_2\partial^2 x_3}\f]
441 template<class T>
442 struct Stencil<T, 0, 2, 2>
443 {
444  [[nodiscard]] AMREX_FORCE_INLINE
445  static T D(const amrex::Array4<const T>& f,
446  const int& i, const int& j, const int& k, const int& m,
447  const Set::Scalar dx[AMREX_SPACEDIM])
448  {
449  return (-(-f(i, j + 2, k + 2, m) + 16.0 * f(i, j + 2, k + 1, m) - 30.0 * f(i, j + 2, k, m) + 16.0 * f(i, j + 2, k - 1, m) - f(i, j + 2, k - 2, m))
450  + 16 * (-f(i, j + 1, k + 2, m) + 16.0 * f(i, j + 1, k + 1, m) - 30.0 * f(i, j + 1, k, m) + 16.0 * f(i, j + 1, k - 1, m) - f(i, j + 1, k - 2, m))
451  - 30 * (-f(i, j, k + 2, m) + 16.0 * f(i, j, k + 1, m) - 30.0 * f(i, j, k, m) + 16.0 * f(i, j, k - 1, m) - f(i, j, k - 2, m))
452  + 16 * (-f(i, j - 1, k + 2, m) + 16.0 * f(i, j - 1, k + 1, m) - 30.0 * f(i, j - 1, k, m) + 16.0 * f(i, j - 1, k - 1, m) - f(i, j - 1, k - 2, m))
453  - (-f(i, j - 2, k + 2, m) + 16.0 * f(i, j - 2, k + 1, m) - 30.0 * f(i, j - 2, k, m) + 16.0 * f(i, j - 2, k - 1, m) - f(i, j - 2, k - 2, m))) /
454  (144.0 * dx[0] * dx[0] * dx[2] * dx[2]);
455  };
456 };
457 /// \brief Compute \f[ \frac{\partial^4}{\partial^2 x_1\partial^2 x_3}\f]
458 template<class T>
459 struct Stencil<T, 2, 0, 2>
460 {
461  [[nodiscard]] AMREX_FORCE_INLINE
462  static T D(const amrex::Array4<const T>& f,
463  const int& i, const int& j, const int& k, const int& m,
464  const Set::Scalar dx[AMREX_SPACEDIM])
465  {
466  return (-(-f(i + 2, j, k + 2, m) + 16.0 * f(i + 2, j, k + 1, m) - 30.0 * f(i + 2, j, k, m) + 16.0 * f(i + 2, j, k - 1, m) - f(i + 2, j, k - 2, m))
467  + 16 * (-f(i + 1, j, k + 2, m) + 16.0 * f(i + 1, j, k + 1, m) - 30.0 * f(i + 1, j, k, m) + 16.0 * f(i + 1, j, k - 1, m) - f(i + 1, j, k - 2, m))
468  - 30 * (-f(i, j, k + 2, m) + 16.0 * f(i, j, k + 1, m) - 30.0 * f(i, j, k, m) + 16.0 * f(i, j, k - 1, m) - f(i, j, k - 2, m))
469  + 16 * (-f(i - 1, j, k + 2, m) + 16.0 * f(i - 1, j, k + 1, m) - 30.0 * f(i - 1, j, k, m) + 16.0 * f(i - 1, j, k - 1, m) - f(i - 1, j, k - 2, m))
470  - (-f(i - 2, j, k + 2, m) + 16.0 * f(i - 2, j, k + 1, m) - 30.0 * f(i - 2, j, k, m) + 16.0 * f(i - 2, j, k - 1, m) - f(i - 2, j, k - 2, m))) /
471  (144.0 * dx[0] * dx[0] * dx[2] * dx[2]);
472  };
473 };
474 
475 
476 /// \brief Compute \f[ \frac{\partial^4}{\partial^2 x_1\partial x_2\partial x_3}\f]
477 template<class T>
478 struct Stencil<T, 2, 1, 1>
479 {
480  [[nodiscard]] AMREX_FORCE_INLINE
481  static T D(const amrex::Array4<const T>& f,
482  const int& i, const int& j, const int& k, const int& m,
483  const Set::Scalar dx[AMREX_SPACEDIM])
484  {
485  return (+(f(i + 1, j + 1, k + 1, m) - 2.0 * f(i, j + 1, k + 1, m) + f(i - 1, j + 1, k + 1, m))
486  - (f(i + 1, j - 1, k + 1, m) - 2.0 * f(i, j - 1, k + 1, m) + f(i - 1, j - 1, k + 1, m))
487  - (f(i + 1, j + 1, k - 1, m) - 2.0 * f(i, j + 1, k - 1, m) + f(i - 1, j + 1, k - 1, m))
488  + (f(i + 1, j - 1, k - 1, m) - 2.0 * f(i, j - 1, k - 1, m) + f(i - 1, j - 1, k - 1, m)))
489  / (4.0 * dx[0] * dx[0] * dx[1] * dx[2]);
490 
491  };
492 };
493 /// \brief Compute \f[ \frac{\partial^4}{\partial x_1\partial^2 x_2\partial x_3}\f]
494 template<class T>
495 struct Stencil<T, 1, 2, 1>
496 {
497  [[nodiscard]] AMREX_FORCE_INLINE
498  static T D(const amrex::Array4<const T>& f,
499  const int& i, const int& j, const int& k, const int& m,
500  const Set::Scalar dx[AMREX_SPACEDIM])
501  {
502  return (+(f(i + 1, j + 1, k + 1, m) - 2.0 * f(i + 1, j, k + 1, m) + f(i + 1, j - 1, k + 1, m))
503  - (f(i - 1, j + 1, k + 1, m) - 2.0 * f(i - 1, j, k + 1, m) + f(i - 1, j - 1, k + 1, m))
504  - (f(i + 1, j + 1, k - 1, m) - 2.0 * f(i + 1, j, k - 1, m) + f(i + 1, j - 1, k - 1, m))
505  + (f(i - 1, j + 1, k - 1, m) - 2.0 * f(i - 1, j, k - 1, m) + f(i - 1, j - 1, k - 1, m)))
506  / (4.0 * dx[0] * dx[1] * dx[1] * dx[2]);
507 
508  };
509 };
510 /// \brief Compute \f[ \frac{\partial^4}{\partial x_1\partial x_2\partial^2 x_3}\f]
511 template<class T>
512 struct Stencil<T, 1, 1, 2>
513 {
514  [[nodiscard]] AMREX_FORCE_INLINE
515  static T D(const amrex::Array4<const T>& f,
516  const int& i, const int& j, const int& k, const int& m,
517  const Set::Scalar dx[AMREX_SPACEDIM])
518  {
519  return (+(f(i + 1, j + 1, k + 1, m) - 2.0 * f(i + 1, j + 1, k, m) + f(i + 1, j + 1, k - 1, m))
520  - (f(i - 1, j + 1, k + 1, m) - 2.0 * f(i - 1, j + 1, k, m) + f(i - 1, j + 1, k - 1, m))
521  - (f(i + 1, j - 1, k + 1, m) - 2.0 * f(i + 1, j - 1, k, m) + f(i + 1, j - 1, k - 1, m))
522  + (f(i - 1, j - 1, k + 1, m) - 2.0 * f(i - 1, j - 1, k, m) + f(i - 1, j - 1, k - 1, m)))
523  / (4.0 * dx[0] * dx[1] * dx[1] * dx[2]);
524 
525  };
526 };
527 
528 [[nodiscard]] AMREX_FORCE_INLINE
529 Set::Scalar
530 Laplacian(const amrex::Array4<const Set::Scalar>& f,
531  const int& i, const int& j, const int& k, const int& m,
532  const Set::Scalar dx[AMREX_SPACEDIM])
533 {
534  Set::Scalar ret = 0.0;
535  ret += (Numeric::Stencil<Set::Scalar, 2, 0, 0>::D(f, i, j, k, m, dx));
536 #if AMREX_SPACEDIM > 1
537  ret += (Numeric::Stencil<Set::Scalar, 0, 2, 0>::D(f, i, j, k, m, dx));
538 #if AMREX_SPACEDIM > 2
539  ret += (Numeric::Stencil<Set::Scalar, 0, 0, 2>::D(f, i, j, k, m, dx));
540 #endif
541 #endif
542  return ret;
543 }
544 
545 [[nodiscard]] AMREX_FORCE_INLINE
546 Set::Vector
547 Laplacian(const amrex::Array4<const Set::Vector>& f,
548  const int& i, const int& j, const int& k,
549  const Set::Scalar dx[AMREX_SPACEDIM])
550 {
551  Set::Vector ret = Set::Vector::Zero();
552  ret += (Numeric::Stencil<Set::Vector, 2, 0, 0>::D(f, i, j, k, 0, dx));
553 #if AMREX_SPACEDIM > 1
554  ret += (Numeric::Stencil<Set::Vector, 0, 2, 0>::D(f, i, j, k, 0, dx));
555 #if AMREX_SPACEDIM > 2
556  ret += (Numeric::Stencil<Set::Vector, 0, 0, 2>::D(f, i, j, k, 0, dx));
557 #endif
558 #endif
559  return ret;
560 }
561 
562 [[nodiscard]] AMREX_FORCE_INLINE
563 Set::Vector
564 Divergence(const amrex::Array4<const Set::Matrix>& dw,
565  const int& i, const int& j, const int& k,
566  const Set::Scalar DX[AMREX_SPACEDIM],
567  std::array<StencilType, AMREX_SPACEDIM>& stencil = DefaultType)
568 {
569  Set::Vector ret = Set::Vector::Zero();
570 
571  if (stencil[0] == StencilType::Central)
572  {
573  AMREX_D_TERM(ret(0) += (dw(i + 1, j, k)(0, 0) - dw(i - 1, j, k)(0, 0)) / 2. / DX[0];,
574  ret(1) += (dw(i + 1, j, k)(1, 0) - dw(i - 1, j, k)(1, 0)) / 2. / DX[0];,
575  ret(2) += (dw(i + 1, j, k)(2, 0) - dw(i - 1, j, k)(2, 0)) / 2. / DX[0];)
576  }
577  else if (stencil[0] == StencilType::Lo)
578  {
579  AMREX_D_TERM(ret(0) += (dw(i, j, k)(0, 0) - dw(i - 1, j, k)(0, 0)) / DX[0];,
580  ret(1) += (dw(i, j, k)(1, 0) - dw(i - 1, j, k)(1, 0)) / DX[0];,
581  ret(2) += (dw(i, j, k)(2, 0) - dw(i - 1, j, k)(2, 0)) / DX[0];)
582  }
583  else if (stencil[0] == StencilType::Hi)
584  {
585  AMREX_D_TERM(ret(0) += (dw(i + 1, j, k)(0, 0) - dw(i, j, k)(0, 0)) / DX[0];,
586  ret(1) += (dw(i + 1, j, k)(1, 0) - dw(i, j, k)(1, 0)) / DX[0];,
587  ret(2) += (dw(i + 1, j, k)(2, 0) - dw(i, j, k)(2, 0)) / DX[0];)
588  }
589 
590 #if AMREX_SPACEDIM > 1
591  if (stencil[1] == StencilType::Central)
592  {
593  AMREX_D_TERM(ret(0) += (dw(i, j + 1, k)(0, 1) - dw(i, j - 1, k)(0, 1)) / 2. / DX[1];,
594  ret(1) += (dw(i, j + 1, k)(1, 1) - dw(i, j - 1, k)(1, 1)) / 2. / DX[1];,
595  ret(2) += (dw(i, j + 1, k)(2, 1) - dw(i, j - 1, k)(2, 1)) / 2. / DX[1];)
596  }
597  else if (stencil[1] == StencilType::Lo)
598  {
599  AMREX_D_TERM(ret(0) += (dw(i, j, k)(0, 1) - dw(i, j - 1, k)(0, 1)) / DX[1];,
600  ret(1) += (dw(i, j, k)(1, 1) - dw(i, j - 1, k)(1, 1)) / DX[1];,
601  ret(2) += (dw(i, j, k)(2, 1) - dw(i, j - 1, k)(2, 1)) / DX[1];)
602  }
603  else if (stencil[1] == StencilType::Hi)
604  {
605  AMREX_D_TERM(ret(0) += (dw(i, j + 1, k)(0, 1) - dw(i, j, k)(0, 1)) / DX[1];,
606  ret(1) += (dw(i, j + 1, k)(1, 1) - dw(i, j, k)(1, 1)) / DX[1];,
607  ret(2) += (dw(i, j + 1, k)(2, 1) - dw(i, j, k)(2, 1)) / DX[1];)
608  }
609 #endif
610 #if AMREX_SPACEDIM > 2
611  if (stencil[2] == StencilType::Central)
612  {
613  AMREX_D_TERM(ret(0) += (dw(i, j, k + 1)(0, 2) - dw(i, j, k - 1)(0, 2)) / 2. / DX[2];,
614  ret(1) += (dw(i, j, k + 1)(1, 2) - dw(i, j, k - 1)(1, 2)) / 2. / DX[2];,
615  ret(2) += (dw(i, j, k + 1)(2, 2) - dw(i, j, k - 1)(2, 2)) / 2. / DX[2];)
616  }
617  else if (stencil[2] == StencilType::Lo)
618  {
619  AMREX_D_TERM(ret(0) += (dw(i, j, k)(0, 2) - dw(i, j, k - 1)(0, 2)) / DX[2];,
620  ret(1) += (dw(i, j, k)(1, 2) - dw(i, j, k - 1)(1, 2)) / DX[2];,
621  ret(2) += (dw(i, j, k)(2, 2) - dw(i, j, k - 1)(2, 2)) / DX[2];)
622  }
623  else if (stencil[2] == StencilType::Hi)
624  {
625  AMREX_D_TERM(ret(0) += (dw(i, j, k + 1)(0, 2) - dw(i, j, k)(0, 2)) / DX[2];,
626  ret(1) += (dw(i, j, k + 1)(1, 2) - dw(i, j, k)(1, 2)) / DX[2];,
627  ret(2) += (dw(i, j, k + 1)(2, 2) - dw(i, j, k)(2, 2)) / DX[2];)
628  }
629 #endif
630  return ret;
631 }
632 
633 
634 [[nodiscard]] AMREX_FORCE_INLINE
635 Set::Scalar
636 Divergence(const amrex::Array4<const Set::Scalar>& f,
637  const int& i, const int& j, const int& k, const int& m,
638  const Set::Scalar dx[AMREX_SPACEDIM],
639  std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
640 {
641  Set::Scalar ret;
642  ret = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, m, dx, stencil));
643 #if AMREX_SPACEDIM > 1
644  ret += (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, m, dx, stencil));
645 #endif
646 #if AMREX_SPACEDIM > 2
647  ret += (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, m, dx, stencil));
648 #endif
649  return ret;
650 }
651 
652 
653 [[nodiscard]] AMREX_FORCE_INLINE
654 Set::Vector
655 Gradient(const amrex::Array4<const Set::Scalar>& f,
656  const int& i, const int& j, const int& k, const int& m,
657  const Set::Scalar dx[AMREX_SPACEDIM],
658  std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
659 {
660  Set::Vector ret;
661  ret(0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, m, dx, stencil));
662 #if AMREX_SPACEDIM > 1
663  ret(1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, m, dx, stencil));
664 #if AMREX_SPACEDIM > 2
665  ret(2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, m, dx, stencil));
666 #endif
667 #endif
668  return ret;
669 }
670 
671 [[nodiscard]] AMREX_FORCE_INLINE
672 Set::Vector
673 CellGradientOnNode(const amrex::Array4<const Set::Scalar>& f,
674  const int& i, const int& j, const int& k, const int& m,
675  const Set::Scalar dx[AMREX_SPACEDIM])
676 {
677  Set::Vector ret;
678 #if AMREX_SPACEDIM == 1
679  ret(0) = (f(i, j, k, m) - f(i - 1, j, k, m)) / dx[0];
680 #elif AMREX_SPACEDIM == 2
681  ret(0) = 0.5 * (f(i, j, k, m) - f(i - 1, j, k, m) + f(i, j - 1, k, m) - f(i - 1, j - 1, k, m)) / dx[0];
682  ret(1) = 0.5 * (f(i, j, k, m) - f(i, j - 1, k, m) + f(i - 1, j, k, m) - f(i - 1, j - 1, k, m)) / dx[1];
683 #elif AMREX_SPACEDIM == 3
684  ret(0) = 0.25 * (f(i, j, k, m) - f(i - 1, j, k, m) + f(i, j - 1, k, m) - f(i - 1, j - 1, k, m) + f(i, j, k - 1, m) - f(i - 1, j, k - 1, m) + f(i, j - 1, k - 1, m) - f(i - 1, j - 1, k - 1, m)) / dx[0];
685  ret(1) = 0.25 * (f(i, j, k, m) - f(i, j - 1, k, m) + f(i - 1, j, k, m) - f(i - 1, j - 1, k, m) + f(i, j, k - 1, m) - f(i, j - 1, k - 1, m) + f(i - 1, j, k - 1, m) - f(i - 1, j - 1, k - 1, m)) / dx[1];
686  ret(2) = 0.25 * (f(i, j, k, m) - f(i, j, k - 1, m) + f(i - 1, j, k, m) - f(i - 1, j, k - 1, m) + f(i, j - 1, k, m) - f(i, j - 1, k - 1, m) + f(i - 1, j - 1, k, m) - f(i - 1, j - 1, k - 1, m)) / dx[2];
687 #endif
688  return ret;
689 }
690 
691 
692 template<class T>
693 [[nodiscard]] AMREX_FORCE_INLINE
694 std::array<T, AMREX_SPACEDIM>
695 CellGradientOnNode(const amrex::Array4<const T>& f,
696  const int& i, const int& j, const int& k, const int& m,
697  const Set::Scalar dx[AMREX_SPACEDIM])
698 {
699  std::array<T, AMREX_SPACEDIM> ret;
700  Util::Abort(INFO);
701 #if AMREX_SPACEDIM == 1
702  ret[0] = (f(i, j, k, m) - f(i - 1, j, k, m)) / dx[0];
703 #elif AMREX_SPACEDIM == 2
704  ret[0] = (f(i, j, k, m) - f(i - 1, j, k, m) + f(i, j - 1, k, m) - f(i - 1, j - 1, k, m)) * 0.5 / dx[0];
705  ret[1] = (f(i, j, k, m) - f(i, j - 1, k, m) + f(i - 1, j, k, m) - f(i - 1, j - 1, k, m)) * 0.5 / dx[1];
706 #elif AMREX_SPACEDIM == 3
707  ret[0] = (f(i, j, k, m) - f(i - 1, j, k, m) + f(i, j - 1, k, m) - f(i - 1, j - 1, k, m) + f(i, j, k - 1, m) - f(i - 1, j, k - 1, m) + f(i, j - 1, k - 1, m) - f(i - 1, j - 1, k - 1, m)) * 0.25 / dx[0];
708  ret[1] = (f(i, j, k, m) - f(i, j - 1, k, m) + f(i - 1, j, k, m) - f(i - 1, j - 1, k, m) + f(i, j, k - 1, m) - f(i, j - 1, k - 1, m) + f(i - 1, j, k - 1, m) - f(i - 1, j - 1, k - 1, m)) * 0.25 / dx[1];
709  ret[2] = (f(i, j, k, m) - f(i, j, k - 1, m) + f(i - 1, j, k, m) - f(i - 1, j, k - 1, m) + f(i, j - 1, k, m) - f(i, j - 1, k - 1, m) + f(i - 1, j - 1, k, m) - f(i - 1, j - 1, k - 1, m)) * 0.25 / dx[2];
710 #endif
711  return ret;
712 }
713 
714 
715 
716 [[nodiscard]] AMREX_FORCE_INLINE
717 Set::Matrix
718 Gradient(const amrex::Array4<const Set::Scalar>& f,
719  const int& i, const int& j, const int& k,
720  const Set::Scalar dx[AMREX_SPACEDIM],
721  std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
722 {
723  Set::Matrix ret;
724  ret(0, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 0, dx, stencil));
725 #if AMREX_SPACEDIM > 1
726  ret(0, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 0, dx, stencil));
727  ret(1, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 1, dx, stencil));
728  ret(1, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 1, dx, stencil));
729 #if AMREX_SPACEDIM > 2
730  ret(0, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 0, dx, stencil));
731  ret(2, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 2, dx, stencil));
732  ret(1, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 1, dx, stencil));
733  ret(2, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 2, dx, stencil));
734  ret(2, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 2, dx, stencil));
735 #endif
736 #endif
737  return ret;
738 }
739 
740 [[nodiscard]] AMREX_FORCE_INLINE
741 Set::Matrix
742 Gradient(const amrex::Array4<const Set::Vector>& f,
743  const int& i, const int& j, const int& k,
744  const Set::Scalar dx[AMREX_SPACEDIM],
745  std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
746 {
747  Set::Matrix ret;
748 
749 #if AMREX_SPACEDIM > 0
750  ret.col(0) = Numeric::Stencil<Set::Vector, 1, 0, 0>::D(f, i, j, k, 0, dx, stencil);
751 #endif
752 #if AMREX_SPACEDIM > 1
753  ret.col(1) = Numeric::Stencil<Set::Vector, 0, 1, 0>::D(f, i, j, k, 0, dx, stencil);
754 #endif
755 #if AMREX_SPACEDIM > 2
756  ret.col(2) = Numeric::Stencil<Set::Vector, 0, 0, 1>::D(f, i, j, k, 0, dx, stencil);
757 #endif
758 
759  return ret;
760 }
761 
762 [[nodiscard]] AMREX_FORCE_INLINE
763 std::pair<Set::Vector,Set::Matrix>
764 GradientSplit( const amrex::Array4<const Set::Vector>& f,
765  const int& i, const int& j, const int& k,
766  const Set::Scalar dx[AMREX_SPACEDIM],
767  std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
768 {
769  Set::Vector diag;
770  Set::Matrix offdiag;
771 
772  std::pair<Set::Scalar,Set::Vector> ret;
773 #if AMREX_SPACEDIM > 0
774  ret = Numeric::Stencil<Set::Vector, 1, 0, 0>::Dsplit(f, i, j, k, 0, dx, stencil);
775  diag(0) = ret.first;
776  offdiag.col(0) = ret.second;
777 #endif
778 #if AMREX_SPACEDIM > 1
779  ret = Numeric::Stencil<Set::Vector, 0, 1, 0>::Dsplit(f, i, j, k, 0, dx, stencil);
780  diag(1) = ret.first;
781  offdiag.col(1) = ret.second;
782 #endif
783 #if AMREX_SPACEDIM > 2
784  ret = Numeric::Stencil<Set::Vector, 0, 0, 1>::Dsplit(f, i, j, k, 0, dx, stencil);
785  diag(2) = ret.first;
786  offdiag.col(2) = ret.second;
787 #endif
788  return std::make_pair(diag,offdiag);
789 }
790 
791 
792 [[nodiscard]] AMREX_FORCE_INLINE
793 Set::Matrix
794 NodeGradientOnCell(const amrex::Array4<const Set::Vector>& f,
795  const int& i, const int& j, const int& k,
796  const Set::Scalar dx[AMREX_SPACEDIM])
797 {
798  Set::Matrix ret;
799 #if AMREX_SPACEDIM == 1
800  ret.col(0) = (f(i + 1, j, k) - f(i, j, k)) / dx[0];
801 #elif AMREX_SPACEDIM == 2
802  ret.col(0) = (f(i + 1, j, k) - f(i, j, k) + f(i + 1, j + 1, k) - f(i, j + 1, k)) * 0.5 / dx[0];
803  ret.col(1) = (f(i, j + 1, k) - f(i, j, k) + f(i + 1, j + 1, k) - f(i + 1, j, k)) * 0.5 / dx[1];
804 #elif AMREX_SPACEDIM == 3
805  ret.col(0) = (f(i + 1, j, k) - f(i, j, k) + f(i + 1, j + 1, k) - f(i, j + 1, k) + f(i + 1, j, k + 1) - f(i, j, k + 1) + f(i + 1, j + 1, k + 1) - f(i, j + 1, k + 1)) * 0.25 / dx[0];
806  ret.col(1) = (f(i, j + 1, k) - f(i, j, k) + f(i + 1, j + 1, k) - f(i + 1, j, k) + f(i, j + 1, k + 1) - f(i, j, k + 1) + f(i + 1, j + 1, k + 1) - f(i + 1, j, k + 1)) * 0.25 / dx[1];
807  ret.col(2) = (f(i, j, k + 1) - f(i, j, k) + f(i, j + 1, k + 1) - f(i, j + 1, k) + f(i + 1, j, k + 1) - f(i + 1, j, k) + f(i + 1, j + 1, k + 1) - f(i + 1, j + 1, k)) * 0.25 / dx[2];
808 #endif
809  return ret;
810 }
811 
812 [[nodiscard]] AMREX_FORCE_INLINE
813 Set::Matrix3
814 Gradient(const amrex::Array4<const Set::Matrix>& f,
815  const int& i, const int& j, const int& k,
816  const Set::Scalar dx[AMREX_SPACEDIM],
817  std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
818 {
819  Set::Matrix3 ret;
820 
821 #if AMREX_SPACEDIM > 0
822  ret[0] = Numeric::Stencil<Set::Matrix, 1, 0, 0>::D(f, i, j, k, 0, dx, stencil);
823 #endif
824 #if AMREX_SPACEDIM > 1
825  ret[1] = Numeric::Stencil<Set::Matrix, 0, 1, 0>::D(f, i, j, k, 0, dx, stencil);
826 #endif
827 #if AMREX_SPACEDIM > 2
828  ret[2] = Numeric::Stencil<Set::Matrix, 0, 0, 1>::D(f, i, j, k, 0, dx, stencil);
829 #endif
830 
831  return ret;
832 }
833 [[nodiscard]] AMREX_FORCE_INLINE
834 Set::Matrix3
835 NodeGradientOnCell(const amrex::Array4<const Set::Matrix>& f,
836  const int& i, const int& j, const int& k,
837  const Set::Scalar dx[AMREX_SPACEDIM])
838 {
839  Set::Matrix3 ret;
840 
841 #if AMREX_SPACEDIM > 0
842  ret[0] = (f(i+1,j,k) - f(i,j,k)) / dx[0];
843 #endif
844 #if AMREX_SPACEDIM > 1
845  ret[1] = (f(i,j+1,k) - f(i,j,k)) / dx[1];
846 #endif
847 #if AMREX_SPACEDIM > 2
848  ret[2] = (f(i,j,k+1) - f(i,j,k)) / dx[2];
849 #endif
850 
851  return ret;
852 }
853 
854 
855 [[nodiscard]] AMREX_FORCE_INLINE
856 Set::Matrix3
857 MatrixGradient(const amrex::Array4<const Set::Scalar>& f,
858  const int& i, const int& j, const int& k,
859  const Set::Scalar dx[AMREX_SPACEDIM],
860  std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
861 {
862  Set::Matrix3 ret;
863 #if AMREX_SPACEDIM == 1
864  ret[0](0, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 0, dx, stencil));
865 #elif AMREX_SPACEDIM == 2
866  ret[0](0, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 0, dx, stencil));
867  ret[0](0, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 0, dx, stencil));
868  ret[0](1, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 1, dx, stencil));
869  ret[0](1, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 1, dx, stencil));
870  ret[1](0, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 2, dx, stencil));
871  ret[1](0, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 2, dx, stencil));
872  ret[1](1, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 3, dx, stencil));
873  ret[1](1, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 3, dx, stencil));
874 #elif AMREX_SPACEDIM == 3
875  ret[0](0, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 0, dx, stencil));
876  ret[0](1, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 1, dx, stencil));
877  ret[0](2, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 2, dx, stencil));
878  ret[1](0, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 3, dx, stencil));
879  ret[1](1, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 4, dx, stencil));
880  ret[1](2, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 5, dx, stencil));
881  ret[2](0, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 6, dx, stencil));
882  ret[2](1, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 7, dx, stencil));
883  ret[2](2, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 8, dx, stencil));
884  ret[0](0, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 0, dx, stencil));
885  ret[0](1, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 1, dx, stencil));
886  ret[0](2, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 2, dx, stencil));
887  ret[1](0, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 3, dx, stencil));
888  ret[1](1, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 4, dx, stencil));
889  ret[1](2, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 5, dx, stencil));
890  ret[2](0, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 6, dx, stencil));
891  ret[2](1, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 7, dx, stencil));
892  ret[2](2, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 8, dx, stencil));
893  ret[0](0, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 0, dx, stencil));
894  ret[0](1, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 1, dx, stencil));
895  ret[0](2, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 2, dx, stencil));
896  ret[1](0, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 3, dx, stencil));
897  ret[1](1, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 4, dx, stencil));
898  ret[1](2, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 5, dx, stencil));
899  ret[2](0, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 6, dx, stencil));
900  ret[2](1, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 7, dx, stencil));
901  ret[2](2, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 8, dx, stencil));
902 #endif
903  return ret;
904 }
905 
906 
907 [[nodiscard]] AMREX_FORCE_INLINE
908 Set::Matrix
909 Hessian(const amrex::Array4<const Set::Scalar>& f,
910  const int& i, const int& j, const int& k, const int& m,
911  const Set::Scalar dx[AMREX_SPACEDIM],
912  std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType
913 )
914 {
915  Set::Matrix ret;
916  ret(0, 0) = (Numeric::Stencil<Set::Scalar, 2, 0, 0>::D(f, i, j, k, m, dx, stencil));
917 #if AMREX_SPACEDIM > 1
918  ret(1, 1) = (Numeric::Stencil<Set::Scalar, 0, 2, 0>::D(f, i, j, k, m, dx, stencil));
919  ret(0, 1) = (Numeric::Stencil<Set::Scalar, 1, 1, 0>::D(f, i, j, k, m, dx, stencil));
920  ret(1, 0) = ret(0, 1);
921 #if AMREX_SPACEDIM > 2
922  ret(2, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 2>::D(f, i, j, k, m, dx, stencil));
923  ret(1, 2) = (Numeric::Stencil<Set::Scalar, 0, 1, 1>::D(f, i, j, k, m, dx, stencil));
924  ret(2, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 1>::D(f, i, j, k, m, dx, stencil));
925  ret(2, 1) = ret(1, 2);
926  ret(0, 2) = ret(2, 0);
927 #endif
928 #endif
929  return ret;
930 }
931 
932 [[nodiscard]] AMREX_FORCE_INLINE
933 Set::Matrix3
934 Hessian(const amrex::Array4<const Set::Scalar>& f,
935  const int& i, const int& j, const int& k,
936  const Set::Scalar DX[AMREX_SPACEDIM],
937  std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
938 {
939  Set::Matrix3 ret;
940  // 1D
941 #if AMREX_SPACEDIM>0
942  ret(0, 0, 0) = (Numeric::Stencil<Set::Scalar, 2, 0, 0>::D(f, i, j, k, 0, DX, stencil));
943 #endif
944  // 2D
945 #if AMREX_SPACEDIM>1
946  ret(0, 0, 1) = (Numeric::Stencil<Set::Scalar, 1, 1, 0>::D(f, i, j, k, 0, DX, stencil));
947  ret(0, 1, 0) = ret(0, 0, 1);
948  ret(0, 1, 1) = (Numeric::Stencil<Set::Scalar, 0, 2, 0>::D(f, i, j, k, 0, DX, stencil));
949  ret(1, 0, 0) = (Numeric::Stencil<Set::Scalar, 2, 0, 0>::D(f, i, j, k, 1, DX, stencil));
950  ret(1, 0, 1) = (Numeric::Stencil<Set::Scalar, 1, 1, 0>::D(f, i, j, k, 1, DX, stencil));
951  ret(1, 1, 0) = ret(1, 0, 1);
952  ret(1, 1, 1) = (Numeric::Stencil<Set::Scalar, 0, 2, 0>::D(f, i, j, k, 1, DX, stencil));
953 #endif
954  // 3D
955 #if AMREX_SPACEDIM>2
956  ret(0, 2, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 2>::D(f, i, j, k, 0, DX, stencil));
957  ret(0, 0, 2) = (Numeric::Stencil<Set::Scalar, 1, 0, 1>::D(f, i, j, k, 0, DX, stencil));
958  ret(0, 1, 2) = (Numeric::Stencil<Set::Scalar, 0, 1, 1>::D(f, i, j, k, 0, DX, stencil));
959  ret(0, 2, 0) = ret(0, 0, 2);
960  ret(0, 2, 1) = ret(0, 1, 2);;
961  ret(1, 2, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 2>::D(f, i, j, k, 1, DX, stencil));
962  ret(1, 0, 2) = (Numeric::Stencil<Set::Scalar, 1, 0, 1>::D(f, i, j, k, 1, DX, stencil));
963  ret(1, 1, 2) = (Numeric::Stencil<Set::Scalar, 0, 1, 1>::D(f, i, j, k, 1, DX, stencil));
964  ret(1, 2, 0) = ret(1, 0, 2);
965  ret(1, 2, 1) = ret(1, 1, 2);;
966  ret(2, 0, 0) = (Numeric::Stencil<Set::Scalar, 2, 0, 0>::D(f, i, j, k, 2, DX, stencil));
967  ret(2, 1, 1) = (Numeric::Stencil<Set::Scalar, 0, 2, 0>::D(f, i, j, k, 2, DX, stencil));
968  ret(2, 2, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 2>::D(f, i, j, k, 2, DX, stencil));
969  ret(2, 0, 1) = (Numeric::Stencil<Set::Scalar, 1, 1, 0>::D(f, i, j, k, 2, DX, stencil));
970  ret(2, 1, 0) = ret(2, 0, 1);
971  ret(2, 0, 2) = (Numeric::Stencil<Set::Scalar, 1, 0, 1>::D(f, i, j, k, 2, DX, stencil));
972  ret(2, 1, 2) = (Numeric::Stencil<Set::Scalar, 0, 1, 1>::D(f, i, j, k, 2, DX, stencil));
973  ret(2, 2, 0) = ret(2, 0, 2);
974  ret(2, 2, 1) = ret(2, 1, 2);;
975 #endif
976  return ret;
977 }
978 
979 
980 // Returns Hessian of a vector field.
981 // Return value: ret[i](j,k) = ret_{i,jk}
982 [[nodiscard]] AMREX_FORCE_INLINE
983 Set::Matrix3
984 Hessian(const amrex::Array4<const Set::Vector>& f,
985  const int& i, const int& j, const int& k,
986  const Set::Scalar dx[AMREX_SPACEDIM],
987  std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
988 {
989  Set::Matrix3 ret;
990 
991 #if AMREX_SPACEDIM>0
992  Set::Vector f_11 = Numeric::Stencil<Set::Vector, 2, 0, 0>::D(f, i, j, k, 0, dx, stencil);
993  ret[0](0, 0) = f_11(0);
994 #endif
995 #if AMREX_SPACEDIM>1
996  ret[1](0, 0) = f_11(1);
997  Set::Vector f_12 = Numeric::Stencil<Set::Vector, 1, 1, 0>::D(f, i, j, k, 0, dx, stencil);
998  ret[0](1, 0) = ret[0](0, 1) = f_12(0);
999  ret[1](1, 0) = ret[1](0, 1) = f_12(1);
1000  Set::Vector f_22 = Numeric::Stencil<Set::Vector, 0, 2, 0>::D(f, i, j, k, 0, dx, stencil);
1001  ret[0](1, 1) = f_22(0);
1002  ret[1](1, 1) = f_22(1);
1003 #endif
1004 #if AMREX_SPACEDIM>2
1005  ret[2](0, 0) = f_11(2);
1006  ret[2](1, 0) = ret[2](0, 1) = f_12(2);
1007  Set::Vector f_13 = Numeric::Stencil<Set::Vector, 1, 0, 1>::D(f, i, j, k, 0, dx, stencil);
1008  ret[0](2, 0) = ret[0](0, 2) = f_13(0);
1009  ret[1](2, 0) = ret[1](0, 2) = f_13(1);
1010  ret[2](2, 0) = ret[2](0, 2) = f_13(2);
1011  ret[2](1, 1) = f_22(2);
1012  Set::Vector f_23 = Numeric::Stencil<Set::Vector, 0, 1, 1>::D(f, i, j, k, 0, dx, stencil);
1013  ret[0](1, 2) = ret[0](2, 1) = f_23(0);
1014  ret[1](1, 2) = ret[1](2, 1) = f_23(1);
1015  ret[2](1, 2) = ret[2](2, 1) = f_23(2);
1016  Set::Vector f_33 = Numeric::Stencil<Set::Vector, 0, 0, 2>::D(f, i, j, k, 0, dx, stencil);
1017  ret[0](2, 2) = f_33(0);
1018  ret[1](2, 2) = f_33(1);
1019  ret[2](2, 2) = f_33(2);
1020 #endif
1021  return ret;
1022 }
1023 
1024 
1025 [[nodiscard]] AMREX_FORCE_INLINE
1026 Set::Matrix
1027 FieldToMatrix(const amrex::Array4<const Set::Scalar>& f,
1028  const int& i, const int& j, const int& k)
1029 {
1030  Set::Matrix ret;
1031 #if AMREX_SPACEDIM == 1
1032  ret(0, 0) = f(i, j, k, 0);
1033 
1034 #elif AMREX_SPACEDIM == 2
1035  ret(0, 0) = f(i, j, k, 0); ret(0, 1) = f(i, j, k, 1);
1036  ret(1, 0) = f(i, j, k, 2); ret(1, 1) = f(i, j, k, 3);
1037 
1038 #elif AMREX_SPACEDIM == 3
1039  ret(0, 0) = f(i, j, k, 0); ret(0, 1) = f(i, j, k, 1); ret(0, 2) = f(i, j, k, 2);
1040  ret(1, 0) = f(i, j, k, 3); ret(1, 1) = f(i, j, k, 4); ret(1, 2) = f(i, j, k, 5);
1041  ret(2, 0) = f(i, j, k, 6); ret(2, 1) = f(i, j, k, 7); ret(2, 2) = f(i, j, k, 8);
1042 #endif
1043 
1044  return ret;
1045 }
1046 
1047 [[nodiscard]] AMREX_FORCE_INLINE
1048 Set::Matrix
1049 FieldToMatrix(const amrex::Array4<Set::Scalar>& f,
1050  const int& i, const int& j, const int& k)
1051 {
1052  Set::Matrix ret;
1053 #if AMREX_SPACEDIM == 1
1054  ret(0, 0) = f(i, j, k, 0);
1055 
1056 #elif AMREX_SPACEDIM == 2
1057  ret(0, 0) = f(i, j, k, 0); ret(0, 1) = f(i, j, k, 1);
1058  ret(1, 0) = f(i, j, k, 2); ret(1, 1) = f(i, j, k, 3);
1059 
1060 #elif AMREX_SPACEDIM == 3
1061  ret(0, 0) = f(i, j, k, 0); ret(0, 1) = f(i, j, k, 1); ret(0, 2) = f(i, j, k, 2);
1062  ret(1, 0) = f(i, j, k, 3); ret(1, 1) = f(i, j, k, 4); ret(1, 2) = f(i, j, k, 5);
1063  ret(2, 0) = f(i, j, k, 6); ret(2, 1) = f(i, j, k, 7); ret(2, 2) = f(i, j, k, 8);
1064 #endif
1065 
1066  return ret;
1067 }
1068 
1069 [[nodiscard]] AMREX_FORCE_INLINE
1070 Set::Vector
1071 FieldToVector(const amrex::Array4<const Set::Scalar>& f,
1072  const int& i, const int& j, const int& k)
1073 {
1074  Set::Vector ret;
1075  ret(0) = f(i, j, k, 0);
1076 #if AMREX_SPACEDIM > 1
1077  ret(1) = f(i, j, k, 1);
1078 #if AMREX_SPACEDIM > 2
1079  ret(2) = f(i, j, k, 2);
1080 #endif
1081 #endif
1082  return ret;
1083 }
1084 
1085 [[nodiscard]] AMREX_FORCE_INLINE
1086 Set::Vector
1087 FieldToVector(const amrex::Array4<Set::Scalar>& f,
1088  const int& i, const int& j, const int& k)
1089 {
1090  Set::Vector ret;
1091  ret(0) = f(i, j, k, 0);
1092 #if AMREX_SPACEDIM > 1
1093  ret(1) = f(i, j, k, 1);
1094 #if AMREX_SPACEDIM > 2
1095  ret(2) = f(i, j, k, 2);
1096 #endif
1097 #endif
1098  return ret;
1099 }
1100 
1101 AMREX_FORCE_INLINE
1102 void
1103 MatrixToField(const amrex::Array4<Set::Scalar>& f,
1104  const int& i, const int& j, const int& k,
1105  Set::Matrix matrix)
1106 {
1107 #if AMREX_SPACEDIM == 1
1108  f(i, j, k, 0) = matrix(0, 0);
1109 #elif AMREX_SPACEDIM == 2
1110  f(i, j, k, 0) = matrix(0, 0); f(i, j, k, 1) = matrix(0, 1);
1111  f(i, j, k, 2) = matrix(1, 0); f(i, j, k, 3) = matrix(1, 1);
1112 #elif AMREX_SPACEDIM == 3
1113  f(i, j, k, 0) = matrix(0, 0); f(i, j, k, 1) = matrix(0, 1); f(i, j, k, 2) = matrix(0, 2);
1114  f(i, j, k, 3) = matrix(1, 0); f(i, j, k, 4) = matrix(1, 1); f(i, j, k, 5) = matrix(1, 2);
1115  f(i, j, k, 6) = matrix(2, 0); f(i, j, k, 7) = matrix(2, 1); f(i, j, k, 8) = matrix(2, 2);
1116 #endif
1117 }
1118 
1119 AMREX_FORCE_INLINE
1120 void
1121 VectorToField(const amrex::Array4<Set::Scalar>& f,
1122  const int& i, const int& j, const int& k,
1123  Set::Vector vector)
1124 {
1125  f(i, j, k, 0) = vector(0);
1126 #if AMREX_SPACEDIM > 1
1127  f(i, j, k, 1) = vector(1);
1128 #if AMREX_SPACEDIM > 2
1129  f(i, j, k, 2) = vector(2);
1130 #endif
1131 #endif
1132 }
1133 
1134 template<int index, int SYM>
1135 [[nodiscard]]
1136 Set::Matrix3
1137 Divergence(const amrex::Array4<const Set::Matrix4<AMREX_SPACEDIM, SYM>>&,
1138  const int, const int, const int, const Set::Scalar[AMREX_SPACEDIM],
1139  std::array<StencilType, AMREX_SPACEDIM> /*stecil*/ = DefaultType)
1140 {
1141  Util::Abort(INFO, "Not implemented yet"); return Set::Matrix3::Zero();
1142 }
1143 
1144 template<>
1145 [[nodiscard]] AMREX_FORCE_INLINE
1146 Set::Matrix3 Divergence<2, Set::Sym::Isotropic>(const amrex::Array4<const Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Isotropic>>& C,
1147  const int i, const int j, const int k, const Set::Scalar dx[AMREX_SPACEDIM],
1148  std::array<StencilType, AMREX_SPACEDIM> stencil)
1149 {
1150  Set::Matrix3 ret;
1151 
1152  Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Isotropic> gradCx =
1153  Numeric::Stencil<Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Isotropic>, 1, 0, 0>::D(C, i, j, k, 0, dx, stencil);
1154 #if AMREX_SPACEDIM>1
1155  Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Isotropic> gradCy =
1156  Numeric::Stencil<Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Isotropic>, 0, 1, 0>::D(C, i, j, k, 0, dx, stencil);
1157 #endif
1158 #if AMREX_SPACEDIM>2
1159  Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Isotropic> gradCz =
1160  Numeric::Stencil<Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Isotropic>, 0, 0, 1>::D(C, i, j, k, 0, dx, stencil);
1161 #endif
1162 
1163  for (int i = 0; i < AMREX_SPACEDIM; i++)
1164  for (int k = 0; k < AMREX_SPACEDIM; k++)
1165  for (int l = 0; l < AMREX_SPACEDIM; l++)
1166  {
1167  ret(i, k, l) = 0.0;
1168  ret(i, k, l) = gradCx(i, 0, k, l);
1169 #if AMREX_SPACEDIM>1
1170  ret(i, k, l) = gradCy(i, 1, k, l);
1171 #endif
1172 #if AMREX_SPACEDIM>2
1173  ret(i, k, l) = gradCz(i, 2, k, l);
1174 #endif
1175  }
1176  return ret;
1177 }
1178 
1179 
1180 
1181 
1182 template<int dim>
1183 [[nodiscard]] AMREX_FORCE_INLINE
1184 Set::Matrix4<dim, Set::Sym::Full>
1185 DoubleHessian(const amrex::Array4<const Set::Scalar>& f,
1186  const int& i, const int& j, const int& k, const int& m,
1187  const Set::Scalar dx[AMREX_SPACEDIM]);
1188 
1189 template<>
1190 [[nodiscard]] AMREX_FORCE_INLINE
1191 Set::Matrix4<2, Set::Sym::Full>
1192 DoubleHessian<2>(const amrex::Array4<const Set::Scalar>& f,
1193  const int& i, const int& j, const int& k, const int& m,
1194  const Set::Scalar dx[AMREX_SPACEDIM])
1195 {
1196  Set::Matrix4<2, Set::Sym::Full> ret;
1197  // [0,0,0,0]
1198  ret(0, 0, 0, 0) = Stencil<Set::Scalar, 4, 0, 0>::D(f, i, j, k, m, dx);
1199  // [0, 0, 0, 1]
1200  ret(0, 0, 0, 1) = Stencil<Set::Scalar, 3, 1, 0>::D(f, i, j, k, m, dx);
1201  // [0, 0, 1, 1]
1202  ret(0, 0, 1, 1) = Stencil<Set::Scalar, 2, 2, 0>::D(f, i, j, k, m, dx);
1203  // [0, 1, 1, 1]
1204  ret(0, 1, 1, 1) = Stencil<Set::Scalar, 1, 3, 0>::D(f, i, j, k, m, dx);
1205  // [1, 1, 1, 1]
1206  ret(1, 1, 1, 1) = Stencil<Set::Scalar, 0, 4, 0>::D(f, i, j, k, m, dx);
1207  return ret;
1208 }
1209 
1210 template<>
1211 [[nodiscard]] AMREX_FORCE_INLINE
1212 Set::Matrix4<3, Set::Sym::Full>
1213 DoubleHessian<3>(const amrex::Array4<const Set::Scalar>& f,
1214  const int& i, const int& j, const int& k, const int& m,
1215  const Set::Scalar dx[AMREX_SPACEDIM])
1216 {
1217  Set::Matrix4<3, Set::Sym::Full> ret;
1218  // [0,0,0,0]
1219  ret(0, 0, 0, 0) = Stencil<Set::Scalar, 4, 0, 0>::D(f, i, j, k, m, dx);
1220  // [0, 0, 0, 1]
1221  ret(0, 0, 0, 1) = Stencil<Set::Scalar, 3, 1, 0>::D(f, i, j, k, m, dx);
1222  // [0, 0, 0, 2]
1223  ret(0, 0, 0, 2) = Stencil<Set::Scalar, 3, 0, 1>::D(f, i, j, k, m, dx);
1224  // [0, 0, 1, 1]
1225  ret(0, 0, 1, 1) = Stencil<Set::Scalar, 2, 2, 0>::D(f, i, j, k, m, dx);
1226  // [0, 0, 1, 2]
1227  ret(0, 0, 1, 2) = Stencil<Set::Scalar, 2, 1, 1>::D(f, i, j, k, m, dx);
1228  // [0, 0, 2, 2]
1229  ret(0, 0, 2, 2) = Stencil<Set::Scalar, 2, 0, 2>::D(f, i, j, k, m, dx);
1230  // [0, 1, 1, 1]
1231  ret(0, 1, 1, 1) = Stencil<Set::Scalar, 1, 3, 0>::D(f, i, j, k, m, dx);
1232  // [0, 1, 1, 2]
1233  ret(0, 1, 1, 2) = Stencil<Set::Scalar, 1, 2, 1>::D(f, i, j, k, m, dx);
1234  // [0, 1, 2, 2]
1235  ret(0, 1, 2, 2) = Stencil<Set::Scalar, 1, 1, 2>::D(f, i, j, k, m, dx);
1236  // [0, 2, 2, 2]
1237  ret(0, 2, 2, 2) = Stencil<Set::Scalar, 1, 0, 3>::D(f, i, j, k, m, dx);
1238  // [1, 1, 1, 1]
1239  ret(1, 1, 1, 1) = Stencil<Set::Scalar, 0, 4, 0>::D(f, i, j, k, m, dx);
1240  // [1, 1, 1, 2]
1241  ret(1, 1, 1, 2) = Stencil<Set::Scalar, 0, 3, 1>::D(f, i, j, k, m, dx);
1242  // [1, 1, 2, 2]
1243  ret(1, 1, 2, 2) = Stencil<Set::Scalar, 0, 2, 2>::D(f, i, j, k, m, dx);
1244  // [1, 2, 2, 2]
1245  ret(1, 2, 2, 2) = Stencil<Set::Scalar, 0, 1, 3>::D(f, i, j, k, m, dx);
1246  // [2, 2, 2, 2]
1247  ret(2, 2, 2, 2) = Stencil<Set::Scalar, 0, 0, 4>::D(f, i, j, k, m, dx);
1248  return ret;
1249 }
1250 
1251 struct Interpolate
1252 {
1253 public:
1254  template<class T>
1255  [[nodiscard]] AMREX_FORCE_INLINE
1256  static T CellToNodeAverage(const amrex::Array4<const T>& f,
1257  const int& i, const int& j, const int& k, const int& m,
1258  std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
1259  {
1260  int AMREX_D_DECL(
1261  ilo = (stencil[0] == Numeric::StencilType::Lo ? 0 : 1),
1262  jlo = (stencil[1] == Numeric::StencilType::Lo ? 0 : 1),
1263  klo = (stencil[2] == Numeric::StencilType::Lo ? 0 : 1));
1264 
1265  return (AMREX_D_TERM(f(i, j, k, m) + f(i - ilo, j, k, m)
1266  ,
1267  +f(i, j - jlo, k, m) + f(i - ilo, j - jlo, k, m)
1268  ,
1269  +f(i, j, k - klo, m) + f(i - ilo, j, k - klo, m)
1270  + f(i, j - jlo, k - klo, m) + f(i - ilo, j - jlo, k - klo, m)
1271  )) * fac;
1272  }
1273  template<class T>
1274  [[nodiscard]] AMREX_FORCE_INLINE
1275  static T CellToNodeAverage(const amrex::Array4<T>& f,
1276  const int& i, const int& j, const int& k, const int& m,
1277  std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
1278  {
1279  int AMREX_D_DECL(
1280  ilo = (stencil[0] == Numeric::StencilType::Lo ? 0 : 1),
1281  jlo = (stencil[1] == Numeric::StencilType::Lo ? 0 : 1),
1282  klo = (stencil[2] == Numeric::StencilType::Lo ? 0 : 1));
1283 
1284  return (AMREX_D_TERM(f(i, j, k, m) + f(i - ilo, j, k, m)
1285  ,
1286  +f(i, j - jlo, k, m) + f(i - ilo, j - jlo, k, m)
1287  ,
1288  +f(i, j, k - klo, m) + f(i - ilo, j, k - klo, m)
1289  + f(i, j - jlo, k - klo, m) + f(i - ilo, j - jlo, k - klo, m)
1290  )) * fac;
1291  }
1292  template<class T>
1293  [[nodiscard]] AMREX_FORCE_INLINE
1294  static T NodeToCellAverage(const amrex::Array4<const T>& f,
1295  const int& i, const int& j, const int& k, const int& m)
1296  {
1297  return (AMREX_D_TERM(f(i, j, k, m) + f(i + 1, j, k, m)
1298  ,
1299  +f(i, j + 1, k, m) + f(i + 1, j + 1, k, m)
1300  ,
1301  +f(i, j, k + 1, m) + f(i + 1, j, k + 1, m)
1302  + f(i, j + 1, k + 1, m) + f(i + 1, j + 1, k + 1, m)
1303  )) * fac;
1304  }
1305  template<class T>
1306  [[nodiscard]] AMREX_FORCE_INLINE
1307  static T NodeToCellAverage(const amrex::Array4<T>& f,
1308  const int& i, const int& j, const int& k, const int& m)
1309  {
1310  return (AMREX_D_TERM(f(i, j, k, m) + f(i + 1, j, k, m)
1311  ,
1312  +f(i, j + 1, k, m) + f(i + 1, j + 1, k, m)
1313  ,
1314  +f(i, j, k + 1, m) + f(i + 1, j, k + 1, m)
1315  + f(i, j + 1, k + 1, m) + f(i + 1, j + 1, k + 1, m)
1316  )) * fac;
1317  }
1318  constexpr static Set::Scalar fac = AMREX_D_PICK(0.5, 0.25, 0.125);
1319 };
1320 
1321 }
1322 #endif
Numeric::CellGradientOnNode
AMREX_FORCE_INLINE Set::Vector CellGradientOnNode(const amrex::Array4< const Set::Scalar > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM])
Definition: Stencil.H:673
Numeric::Interpolate::CellToNodeAverage
static AMREX_FORCE_INLINE T CellToNodeAverage(const amrex::Array4< const T > &f, const int &i, const int &j, const int &k, const int &m, std::array< StencilType, AMREX_SPACEDIM > stencil=DefaultType)
Definition: Stencil.H:1256
Numeric::GradientSplit
AMREX_FORCE_INLINE std::pair< Set::Vector, Set::Matrix > GradientSplit(const amrex::Array4< const Set::Vector > &f, const int &i, const int &j, const int &k, const Set::Scalar dx[AMREX_SPACEDIM], std::array< StencilType, AMREX_SPACEDIM > stencil=DefaultType)
Definition: Stencil.H:764
Numeric::MatrixToField
AMREX_FORCE_INLINE void MatrixToField(const amrex::Array4< Set::Scalar > &f, const int &i, const int &j, const int &k, Set::Matrix matrix)
Definition: Stencil.H:1103
Numeric::FieldToVector
AMREX_FORCE_INLINE Set::Vector FieldToVector(const amrex::Array4< const Set::Scalar > &f, const int &i, const int &j, const int &k)
Definition: Stencil.H:1071
Numeric::Interpolate::NodeToCellAverage
static AMREX_FORCE_INLINE T NodeToCellAverage(const amrex::Array4< T > &f, const int &i, const int &j, const int &k, const int &m)
Definition: Stencil.H:1307
Numeric::Stencil< T, 0, 3, 1 >::D
static AMREX_FORCE_INLINE T D(const amrex::Array4< const T > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM])
Definition: Stencil.H:361
Numeric::Stencil< T, 1, 1, 0 >::D
static AMREX_FORCE_INLINE T D(const amrex::Array4< const T > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM], std::array< StencilType, AMREX_SPACEDIM > stencil=DefaultType)
Definition: Stencil.H:224
BC::AMREX_D_DECL
@ AMREX_D_DECL
Definition: BC.H:34
Numeric::Stencil< T, 2, 0, 0 >::D
static AMREX_FORCE_INLINE T D(const amrex::Array4< const T > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM], std::array< StencilType, AMREX_SPACEDIM > stencil=DefaultType)
Definition: Stencil.H:176
Numeric::VectorToField
AMREX_FORCE_INLINE void VectorToField(const amrex::Array4< Set::Scalar > &f, const int &i, const int &j, const int &k, Set::Vector vector)
Definition: Stencil.H:1121
Numeric::Stencil< T, 0, 1, 0 >::Dsplit
static AMREX_FORCE_INLINE std::pair< Set::Scalar, T > Dsplit(const amrex::Array4< const T > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM], std::array< StencilType, AMREX_SPACEDIM > stencil=DefaultType)
Definition: Stencil.H:104
Numeric::Stencil< T, 3, 0, 1 >::D
static AMREX_FORCE_INLINE T D(const amrex::Array4< const T > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM])
Definition: Stencil.H:410
Numeric::Stencil< T, 0, 1, 3 >::D
static AMREX_FORCE_INLINE T D(const amrex::Array4< const T > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM])
Definition: Stencil.H:377
Numeric::Interpolate::CellToNodeAverage
static AMREX_FORCE_INLINE T CellToNodeAverage(const amrex::Array4< T > &f, const int &i, const int &j, const int &k, const int &m, std::array< StencilType, AMREX_SPACEDIM > stencil=DefaultType)
Definition: Stencil.H:1275
Numeric::Interpolate::fac
constexpr static Set::Scalar fac
Definition: Stencil.H:1318
Numeric::Stencil< T, 0, 2, 2 >::D
static AMREX_FORCE_INLINE T D(const amrex::Array4< const T > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM])
Definition: Stencil.H:445
Numeric::Stencil< T, 0, 0, 2 >::D
static AMREX_FORCE_INLINE T D(const amrex::Array4< const T > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM], std::array< StencilType, AMREX_SPACEDIM > stencil=DefaultType)
Definition: Stencil.H:208
Numeric::Stencil< T, 1, 0, 0 >::Dsplit
static AMREX_FORCE_INLINE std::pair< Set::Scalar, T > Dsplit(const amrex::Array4< const T > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM], std::array< StencilType, AMREX_SPACEDIM > stencil=DefaultType)
Definition: Stencil.H:63
Set::Vector
Eigen::Matrix< amrex::Real, AMREX_SPACEDIM, 1 > Vector
Definition: Base.H:20
Numeric::Gradient
AMREX_FORCE_INLINE Set::Vector Gradient(const amrex::Array4< const Set::Scalar > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM], std::array< StencilType, AMREX_SPACEDIM > stencil=DefaultType)
Definition: Stencil.H:655
Numeric::Stencil< T, 2, 1, 1 >::D
static AMREX_FORCE_INLINE T D(const amrex::Array4< const T > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM])
Definition: Stencil.H:481
Numeric::Divergence
AMREX_FORCE_INLINE Set::Vector Divergence(const amrex::Array4< const Set::Matrix > &dw, const int &i, const int &j, const int &k, const Set::Scalar DX[AMREX_SPACEDIM], std::array< StencilType, AMREX_SPACEDIM > &stencil=DefaultType)
Definition: Stencil.H:564
Numeric::Stencil< T, 1, 1, 2 >::D
static AMREX_FORCE_INLINE T D(const amrex::Array4< const T > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM])
Definition: Stencil.H:515
Numeric::Hi
@ Hi
Definition: Stencil.H:12
Numeric::Interpolate::NodeToCellAverage
static AMREX_FORCE_INLINE T NodeToCellAverage(const amrex::Array4< const T > &f, const int &i, const int &j, const int &k, const int &m)
Definition: Stencil.H:1294
Set::Scalar
amrex::Real Scalar
Definition: Base.H:19
Numeric::Central
@ Central
Definition: Stencil.H:12
Set::Matrix3
Definition: Matrix3.H:8
Numeric::Stencil< T, 1, 0, 3 >::D
static AMREX_FORCE_INLINE T D(const amrex::Array4< const T > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM])
Definition: Stencil.H:393
Set::Matrix
Eigen::Matrix< amrex::Real, AMREX_SPACEDIM, AMREX_SPACEDIM > Matrix
Definition: Base.H:23
Set::Matrix3::Zero
static Matrix3 Zero()
Definition: Matrix3.H:42
Numeric::MatrixGradient
AMREX_FORCE_INLINE Set::Matrix3 MatrixGradient(const amrex::Array4< const Set::Scalar > &f, const int &i, const int &j, const int &k, const Set::Scalar dx[AMREX_SPACEDIM], std::array< StencilType, AMREX_SPACEDIM > stencil=DefaultType)
Definition: Stencil.H:857
Numeric::Laplacian
AMREX_FORCE_INLINE Set::Scalar Laplacian(const amrex::Array4< const Set::Scalar > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM])
Definition: Stencil.H:530
Numeric::DoubleHessian< 2 >
AMREX_FORCE_INLINE Set::Matrix4< 2, Set::Sym::Full > DoubleHessian< 2 >(const amrex::Array4< const Set::Scalar > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM])
Definition: Stencil.H:1192
Numeric::Stencil< T, 4, 0, 0 >::D
static AMREX_FORCE_INLINE T D(const amrex::Array4< const T > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM])
Definition: Stencil.H:286
Numeric::Stencil< T, 0, 0, 1 >::Dsplit
static AMREX_FORCE_INLINE std::pair< Set::Scalar, T > Dsplit(const amrex::Array4< const T > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM], std::array< StencilType, AMREX_SPACEDIM > stencil=DefaultType)
Definition: Stencil.H:145
Numeric::StencilType
StencilType
Definition: Stencil.H:12
Numeric::Stencil< T, 3, 1, 0 >::D
static AMREX_FORCE_INLINE T D(const amrex::Array4< const T > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM])
Definition: Stencil.H:325
Util::Abort
void Abort(const char *msg)
Definition: Util.cpp:165
Numeric::Stencil< T, 0, 0, 1 >::D
static AMREX_FORCE_INLINE T D(const amrex::Array4< const T > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM], std::array< StencilType, AMREX_SPACEDIM > stencil=DefaultType)
Definition: Stencil.H:131
Numeric::NodeGradientOnCell
AMREX_FORCE_INLINE Set::Matrix NodeGradientOnCell(const amrex::Array4< const Set::Vector > &f, const int &i, const int &j, const int &k, const Set::Scalar dx[AMREX_SPACEDIM])
Definition: Stencil.H:794
Set::Matrix4< AMREX_SPACEDIM, SYM >
Numeric::Interpolate
Definition: Stencil.H:1251
Numeric::Stencil< T, 0, 4, 0 >::D
static AMREX_FORCE_INLINE T D(const amrex::Array4< const T > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM])
Definition: Stencil.H:298
Numeric::Stencil< T, 1, 2, 1 >::D
static AMREX_FORCE_INLINE T D(const amrex::Array4< const T > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM])
Definition: Stencil.H:498
Numeric::DoubleHessian< 3 >
AMREX_FORCE_INLINE Set::Matrix4< 3, Set::Sym::Full > DoubleHessian< 3 >(const amrex::Array4< const Set::Scalar > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM])
Definition: Stencil.H:1213
Set.H
Numeric::Stencil< T, 2, 2, 0 >::D
static AMREX_FORCE_INLINE T D(const amrex::Array4< const T > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM])
Definition: Stencil.H:427
Numeric::Hessian
AMREX_FORCE_INLINE Set::Matrix Hessian(const amrex::Array4< const Set::Scalar > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM], std::array< StencilType, AMREX_SPACEDIM > stencil=DefaultType)
Definition: Stencil.H:909
Numeric::FieldToMatrix
AMREX_FORCE_INLINE Set::Matrix FieldToMatrix(const amrex::Array4< const Set::Scalar > &f, const int &i, const int &j, const int &k)
Definition: Stencil.H:1027
Numeric::Stencil< T, 0, 1, 0 >::D
static AMREX_FORCE_INLINE T D(const amrex::Array4< const T > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM], std::array< StencilType, AMREX_SPACEDIM > stencil=DefaultType)
Definition: Stencil.H:90
Numeric::Stencil< T, 0, 0, 4 >::D
static AMREX_FORCE_INLINE T D(const amrex::Array4< const T > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM])
Definition: Stencil.H:310
Numeric::Stencil< T, 1, 3, 0 >::D
static AMREX_FORCE_INLINE T D(const amrex::Array4< const T > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM])
Definition: Stencil.H:343
Numeric::DoubleHessian
AMREX_FORCE_INLINE Set::Matrix4< dim, Set::Sym::Full > DoubleHessian(const amrex::Array4< const Set::Scalar > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM])
Numeric::Stencil< T, 2, 0, 2 >::D
static AMREX_FORCE_INLINE T D(const amrex::Array4< const T > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM])
Definition: Stencil.H:462
Numeric::Stencil
Definition: Stencil.H:38
Numeric::Lo
@ Lo
Definition: Stencil.H:12
Numeric::Stencil< T, 0, 1, 1 >::D
static AMREX_FORCE_INLINE T D(const amrex::Array4< const T > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM], std::array< StencilType, AMREX_SPACEDIM > stencil=DefaultType)
Definition: Stencil.H:262
Numeric::DefaultType
static std::array< StencilType, AMREX_SPACEDIM > DefaultType
Definition: Stencil.H:14
Numeric
This namespace contains some numerical tools.
Definition: Function.H:5
INFO
#define INFO
Definition: Util.H:20
Numeric::Stencil< T, 1, 0, 1 >::D
static AMREX_FORCE_INLINE T D(const amrex::Array4< const T > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM], std::array< StencilType, AMREX_SPACEDIM > stencil=DefaultType)
Definition: Stencil.H:243
Numeric::GetStencil
static AMREX_FORCE_INLINE std::array< StencilType, AMREX_SPACEDIM > GetStencil(const int i, const int j, const int k, const amrex::Box domain)
Definition: Stencil.H:20
Numeric::Stencil< T, 1, 0, 0 >::D
static AMREX_FORCE_INLINE T D(const amrex::Array4< const T > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM], std::array< StencilType, AMREX_SPACEDIM > stencil=DefaultType)
Definition: Stencil.H:49
Numeric::Stencil< T, 0, 2, 0 >::D
static AMREX_FORCE_INLINE T D(const amrex::Array4< const T > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM], std::array< StencilType, AMREX_SPACEDIM > stencil=DefaultType)
Definition: Stencil.H:192
Generated on Fri Aug 16 2024 16:30:29 for Alamo by   doxygen 1.8.17