Block-Structured AMR Software Framework
 
Loading...
Searching...
No Matches
AMReX_EB2_2D_C.H
Go to the documentation of this file.
1#ifndef AMREX_EB2_2D_C_H_
2#define AMREX_EB2_2D_C_H_
3#include <AMReX_Config.H>
4
5namespace amrex::EB2 {
6
8void
9amrex_eb2_build_types (Box const& tbx, Box const& bxg2,
10 Array4<Real const> const& s,
11 Array4<EBCellFlag> const& cell,
12 Array4<Type_t> const& fx,
13 Array4<Type_t> const& fy)
14{
15 auto lo = amrex::max_lbound(tbx, bxg2);
16 auto hi = amrex::min_ubound(tbx, bxg2);
17 amrex::Loop(lo, hi,
18 [=] (int i, int j, int k) noexcept
19 {
20 if ( s(i,j ,k) < 0.0_rt && s(i+1,j ,k) < 0.0_rt
21 && s(i,j+1,k) < 0.0_rt && s(i+1,j+1,k) < 0.0_rt)
22 {
23 cell(i,j,k).setRegular();
24 }
25 else if (s(i,j ,k) >= 0.0_rt && s(i+1,j ,k) >= 0.0_rt
26 && s(i,j+1,k) >= 0.0_rt && s(i+1,j+1,k) >= 0.0_rt)
27 {
28 cell(i,j,k).setCovered();
29 }
30 else
31 {
32 cell(i,j,k).setSingleValued();
33 }
34 });
35
36 // x-face
37 const Box& xbx = amrex::surroundingNodes(bxg2,0);
38 lo = amrex::max_lbound(tbx, xbx);
39 hi = amrex::min_ubound(tbx, xbx);
40 amrex::Loop(lo, hi,
41 [=] (int i, int j, int k) noexcept
42 {
43 if (s(i,j,k) < 0.0_rt && s(i,j+1,k) < 0.0_rt) {
44 fx(i,j,k) = Type::regular;
45 } else if (s(i,j,k) >= 0.0_rt && s(i,j+1,k) >= 0.0_rt) {
46 fx(i,j,k) = Type::covered;
47 } else {
48 fx(i,j,k) = Type::irregular;
49 }
50 });
51
52 // y-face
53 const Box& ybx = amrex::surroundingNodes(bxg2,1);
54 lo = amrex::max_lbound(tbx, ybx);
55 hi = amrex::min_ubound(tbx, ybx);
56 amrex::Loop(lo, hi,
57 [=] (int i, int j, int k) noexcept
58 {
59 if (s(i,j,k) < 0.0_rt && s(i+1,j,k) < 0.0_rt) {
60 fy(i,j,k) = Type::regular;
61 } else if (s(i,j,k) >= 0.0_rt && s(i+1,j,k) >= 0.0_rt) {
62 fy(i,j,k) = Type::covered;
63 } else {
64 fy(i,j,k) = Type::irregular;
65 }
66 });
67}
68
69namespace detail {
71 int num_cuts (Real a, Real b) noexcept {
72 return (a >= 0.0_rt && b < 0.0_rt) || (b >= 0.0_rt && a < 0.0_rt);
73 }
74}
75
77int check_mvmc (int i, int j, int, Array4<Real const> const& fine)
78{
79 using detail::num_cuts;
80
81 constexpr int k = 0;
82 i *= 2;
83 j *= 2;
84 int ncuts = num_cuts(fine(i ,j ,k),fine(i+1,j ,k))
85 + num_cuts(fine(i+1,j ,k),fine(i+2,j ,k))
86 + num_cuts(fine(i ,j+2,k),fine(i+1,j+2,k))
87 + num_cuts(fine(i+1,j+2,k),fine(i+2,j+2,k))
88 + num_cuts(fine(i ,j ,k),fine(i ,j+1,k))
89 + num_cuts(fine(i ,j+1,k),fine(i ,j+2,k))
90 + num_cuts(fine(i+2,j ,k),fine(i+2,j+1,k))
91 + num_cuts(fine(i+2,j+1,k),fine(i+2,j+2,k));
92 return (ncuts != 0 && ncuts != 2);
93}
94
96int coarsen_from_fine (int i, int j, Box const& bx, int ngrow,
97 Array4<Real> const& cvol, Array4<Real> const& ccent,
98 Array4<Real> const& cba, Array4<Real> const& cbc,
99 Array4<Real> const& cbn, Array4<Real> const& capx,
100 Array4<Real> const& capy, Array4<Real> const& cfcx,
101 Array4<Real> const& cfcy, Array4<Real> const& cecx,
102 Array4<Real> const& cecy, Array4<EBCellFlag> const& cflag,
103 Array4<Real const> const& fvol, Array4<Real const> const& fcent,
104 Array4<Real const> const& fba, Array4<Real const> const& fbc,
105 Array4<Real const> const& fbn, Array4<Real const> const& fapx,
106 Array4<Real const> const& fapy, Array4<Real const> const& ffcx,
107 Array4<Real const> const& ffcy, Array4<Real const> const& fecx,
108 Array4<Real const> const& fecy, Array4<EBCellFlag const> const& fflag)
109{
110 const Box& gbx = amrex::grow(bx,ngrow);
111 const Box& xbx = amrex::surroundingNodes(bx,0);
112 const Box& ybx = amrex::surroundingNodes(bx,1);
113 const Box& xgbx = amrex::surroundingNodes(gbx,0);
114 const Box& ygbx = amrex::surroundingNodes(gbx,1);
115
116 int ierr = 0;
117 constexpr int k = 0;
118 int ii = i*2;
119 int jj = j*2;
120
121 if (bx.contains(i,j,k))
122 {
123 if (fflag(ii,jj ,k).isRegular() && fflag(ii+1,jj ,k).isRegular() &&
124 fflag(ii,jj+1,k).isRegular() && fflag(ii+1,jj+1,k).isRegular())
125 {
126 cflag(i,j,k).setRegular();
127 cvol(i,j,k) = 1.0_rt;
128 ccent(i,j,k,0) = 0.0_rt;
129 ccent(i,j,k,1) = 0.0_rt;
130 cba(i,j,k) = 0.0_rt;
131 cbc(i,j,k,0) = -1.0_rt;
132 cbc(i,j,k,1) = -1.0_rt;
133 cbn(i,j,k,0) = 0.0_rt;
134 cbn(i,j,k,1) = 0.0_rt;
135 }
136 else if (fflag(ii,jj ,k).isCovered() && fflag(ii+1,jj ,k).isCovered() &&
137 fflag(ii,jj+1,k).isCovered() && fflag(ii+1,jj+1,k).isCovered())
138 {
139 cflag(i,j,k).setCovered();
140 cvol(i,j,k) = 0.0_rt;
141 ccent(i,j,k,0) = 0.0_rt;
142 ccent(i,j,k,1) = 0.0_rt;
143 cba(i,j,k) = 0.0_rt;
144 cbc(i,j,k,0) = -1.0_rt;
145 cbc(i,j,k,1) = -1.0_rt;
146 cbn(i,j,k,0) = 0.0_rt;
147 cbn(i,j,k,1) = 0.0_rt;
148 }
149 else
150 {
151 cflag(i,j,k).setSingleValued();
152
153 cvol(i,j,k) = 0.25_rt*(fvol(ii,jj ,k) + fvol(ii+1,jj ,k) +
154 fvol(ii,jj+1,k) + fvol(ii+1,jj+1,k));
155 Real cvolinv = 1.0_rt/cvol(i,j,k);
156
157 ccent(i,j,k,0) = 0.25_rt * cvolinv *
158 (fvol(ii ,jj ,k)*(0.5_rt*fcent(ii ,jj ,k,0)-0.25_rt) +
159 fvol(ii+1,jj ,k)*(0.5_rt*fcent(ii+1,jj ,k,0)+0.25_rt) +
160 fvol(ii ,jj+1,k)*(0.5_rt*fcent(ii ,jj+1,k,0)-0.25_rt) +
161 fvol(ii+1,jj+1,k)*(0.5_rt*fcent(ii+1,jj+1,k,0)+0.25_rt));
162 ccent(i,j,k,1) = 0.25_rt * cvolinv *
163 (fvol(ii ,jj ,k)*(0.5_rt*fcent(ii ,jj ,k,1)-0.25_rt) +
164 fvol(ii+1,jj ,k)*(0.5_rt*fcent(ii+1,jj ,k,1)-0.25_rt) +
165 fvol(ii ,jj+1,k)*(0.5_rt*fcent(ii ,jj+1,k,1)+0.25_rt) +
166 fvol(ii+1,jj+1,k)*(0.5_rt*fcent(ii+1,jj+1,k,1)+0.25_rt));
167
168 cba(i,j,k) = 0.5_rt*(fba(ii,jj ,k) + fba(ii+1,jj ,k) +
169 fba(ii,jj+1,k) + fba(ii+1,jj+1,k));
170 Real cbainv = 1.0_rt/cba(i,j,k);
171
172 cbc(i,j,k,0) = 0.5_rt * cbainv *
173 (fba(ii ,jj ,k)*(0.5_rt*fbc(ii ,jj ,k,0)-0.25_rt) +
174 fba(ii+1,jj ,k)*(0.5_rt*fbc(ii+1,jj ,k,0)+0.25_rt) +
175 fba(ii ,jj+1,k)*(0.5_rt*fbc(ii ,jj+1,k,0)-0.25_rt) +
176 fba(ii+1,jj+1,k)*(0.5_rt*fbc(ii+1,jj+1,k,0)+0.25_rt));
177 cbc(i,j,k,1) = 0.5_rt * cbainv *
178 (fba(ii ,jj ,k)*(0.5_rt*fbc(ii ,jj ,k,1)-0.25_rt) +
179 fba(ii+1,jj ,k)*(0.5_rt*fbc(ii+1,jj ,k,1)-0.25_rt) +
180 fba(ii ,jj+1,k)*(0.5_rt*fbc(ii ,jj+1,k,1)+0.25_rt) +
181 fba(ii+1,jj+1,k)*(0.5_rt*fbc(ii+1,jj+1,k,1)+0.25_rt));
182
183 Real nx = fbn(ii ,jj ,k,0)*fba(ii ,jj ,k)
184 + fbn(ii+1,jj ,k,0)*fba(ii+1,jj ,k)
185 + fbn(ii ,jj+1,k,0)*fba(ii ,jj+1,k)
186 + fbn(ii+1,jj+1,k,0)*fba(ii+1,jj+1,k);
187 Real ny = fbn(ii ,jj ,k,1)*fba(ii ,jj ,k)
188 + fbn(ii+1,jj ,k,1)*fba(ii+1,jj ,k)
189 + fbn(ii ,jj+1,k,1)*fba(ii ,jj+1,k)
190 + fbn(ii+1,jj+1,k,1)*fba(ii+1,jj+1,k);
191 Real nfac = 1.0_rt/std::sqrt(nx*nx+ny*ny+1.e-30_rt);
192 cbn(i,j,k,0) = nx*nfac;
193 cbn(i,j,k,1) = ny*nfac;
194 ierr = (nx == 0.0_rt && ny == 0.0_rt)
195 // we must check the enclosing surface to make sure the coarse cell does not
196 // fully contains the geometry object.
197 || ( fapx(ii,jj ,k)==1.0_rt && fapx(ii+2,jj ,k)==1.0_rt
198 && fapx(ii,jj+1,k)==1.0_rt && fapx(ii+2,jj+1,k)==1.0_rt
199 && fapy(ii,jj ,k)==1.0_rt && fapy(ii+1,jj ,k)==1.0_rt
200 && fapy(ii,jj+2,k)==1.0_rt && fapy(ii+1,jj+2,k)==1.0_rt);
201 }
202 }
203 else if (gbx.contains(i,j,k))
204 {
205 cvol(i,j,k) = 1.0_rt;
206 ccent(i,j,k,0) = 0.0_rt;
207 ccent(i,j,k,1) = 0.0_rt;
208 cba(i,j,k) = 0.0_rt;
209 cbc(i,j,k,0) = -1.0_rt;
210 cbc(i,j,k,1) = -1.0_rt;
211 cbn(i,j,k,0) = 0.0_rt;
212 cbn(i,j,k,1) = 0.0_rt;
213 }
214
215 if (xbx.contains(i,j,k))
216 {
217 capx(i,j,k) = 0.5_rt*(fapx(ii,jj,k) + fapx(ii,jj+1,k));
218 if (capx(i,j,k) != 0.0_rt) {
219 cfcx(i,j,k) = 0.5_rt / capx(i,j,k) *
220 (fapx(ii,jj ,k)*(0.5_rt*ffcx(ii,jj ,k)-0.25_rt) +
221 fapx(ii,jj+1,k)*(0.5_rt*ffcx(ii,jj+1,k)+0.25_rt));
222 if (fecy(ii,jj,k) == Real(1.0) && fecy(ii,jj+1,k) == Real(1.0)) {
223 cecy(i,j,k) = Real(1.0);
224 } else {
225 cecy(i,j,k) = cfcx(i,j,k);
226 }
227 }
228 else {
229 cfcx(i,j,k) = 0.0_rt;
230 cecy(i,j,k) = -1.0_rt;
231 }
232 }
233 else if (xgbx.contains(i,j,k))
234 {
235 capx(i,j,k) = 1.0_rt;
236 cfcx(i,j,k) = 0.0_rt;
237 cecy(i,j,k) = 1.0_rt;
238 }
239
240 if (ybx.contains(i,j,k))
241 {
242 capy(i,j,k) = 0.5_rt*(fapy(ii,jj,k) + fapy(ii+1,jj,k));
243 if (capy(i,j,k) != 0.0_rt) {
244 cfcy(i,j,k) = 0.5_rt / capy(i,j,k) *
245 (fapy(ii ,jj,k)*(0.5_rt*ffcy(ii ,jj,k)-0.25_rt) +
246 fapy(ii+1,jj,k)*(0.5_rt*ffcy(ii+1,jj,k)+0.25_rt));
247 if (fecx(ii,jj,k) == Real(1.0) && fecx(ii+1,jj,k) == Real(1.0)) {
248 cecx(i,j,k) = Real(1.0);
249 } else {
250 cecx(i,j,k) = cfcy(i,j,k);
251 }
252 } else {
253 cfcy(i,j,k) = 0.0_rt;
254 cecx(i,j,k) = -1.0_rt;
255 }
256 }
257 else if (ygbx.contains(i,j,k))
258 {
259 capy(i,j,k) = 1.0_rt;
260 cfcy(i,j,k) = 0.0_rt;
261 cecx(i,j,k) = 1.0_rt;
262 }
263
264 return ierr;
265}
266
268void build_cellflag_from_ap (int i, int j, Array4<EBCellFlag> const& cflag,
269 Array4<Real const> const& apx, Array4<Real const> const& apy)
270{
271 constexpr int k = 0;
272
273 // By default, all neighbors are already set.
274 auto flg = cflag(i,j,k);
275
276 if (apx(i ,j ,k) == 0.0_rt) { flg.setDisconnected(-1, 0, 0); }
277 if (apx(i+1,j ,k) == 0.0_rt) { flg.setDisconnected( 1, 0, 0); }
278 if (apy(i ,j ,k) == 0.0_rt) { flg.setDisconnected( 0,-1, 0); }
279 if (apy(i ,j+1,k) == 0.0_rt) { flg.setDisconnected( 0, 1, 0); }
280
281 if ((apx(i,j ,k) == 0.0_rt || apy(i-1,j,k) == 0.0_rt) &&
282 (apx(i,j-1,k) == 0.0_rt || apy(i ,j,k) == 0.0_rt))
283 {
284 flg.setDisconnected(-1,-1,0);
285 }
286
287 if ((apx(i+1,j ,k) == 0.0_rt || apy(i+1,j,k) == 0.0_rt) &&
288 (apx(i+1,j-1,k) == 0.0_rt || apy(i ,j,k) == 0.0_rt))
289 {
290 flg.setDisconnected(1,-1,0);
291 }
292
293 if ((apx(i,j ,k) == 0.0_rt || apy(i-1,j+1,k) == 0.0_rt) &&
294 (apx(i,j+1,k) == 0.0_rt || apy(i ,j+1,k) == 0.0_rt))
295 {
296 flg.setDisconnected(-1,1,0);
297 }
298
299 if ((apx(i+1,j ,k) == 0.0_rt || apy(i+1,j+1,k) == 0.0_rt) &&
300 (apx(i+1,j+1,k) == 0.0_rt || apy(i ,j+1,k) == 0.0_rt))
301 {
302 flg.setDisconnected(1,1,0);
303 }
304
305 cflag(i,j,k) = flg;
306}
307
308}
309
310#endif
#define AMREX_FORCE_INLINE
Definition AMReX_Extension.H:119
#define AMREX_GPU_HOST_DEVICE
Definition AMReX_GpuQualifiers.H:20
Array4< Real > fine
Definition AMReX_InterpFaceRegister.cpp:90
AMREX_GPU_HOST_DEVICE bool contains(const IntVectND< dim > &p) const noexcept
Returns true if argument is contained within BoxND.
Definition AMReX_Box.H:204
static constexpr Type_t regular
Definition AMReX_EB2_Graph.H:38
static constexpr Type_t covered
Definition AMReX_EB2_Graph.H:39
static constexpr Type_t irregular
Definition AMReX_EB2_Graph.H:40
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int num_cuts(Real a, Real b) noexcept
Definition AMReX_EB2_2D_C.H:71
Definition AMReX_FabArrayBase.H:32
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int coarsen_from_fine(int i, int j, Box const &bx, int ngrow, Array4< Real > const &cvol, Array4< Real > const &ccent, Array4< Real > const &cba, Array4< Real > const &cbc, Array4< Real > const &cbn, Array4< Real > const &capx, Array4< Real > const &capy, Array4< Real > const &cfcx, Array4< Real > const &cfcy, Array4< Real > const &cecx, Array4< Real > const &cecy, Array4< EBCellFlag > const &cflag, Array4< Real const > const &fvol, Array4< Real const > const &fcent, Array4< Real const > const &fba, Array4< Real const > const &fbc, Array4< Real const > const &fbn, Array4< Real const > const &fapx, Array4< Real const > const &fapy, Array4< Real const > const &ffcx, Array4< Real const > const &ffcy, Array4< Real const > const &fecx, Array4< Real const > const &fecy, Array4< EBCellFlag const > const &fflag)
Definition AMReX_EB2_2D_C.H:96
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void build_cellflag_from_ap(int i, int j, Array4< EBCellFlag > const &cflag, Array4< Real const > const &apx, Array4< Real const > const &apy)
Definition AMReX_EB2_2D_C.H:268
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void amrex_eb2_build_types(Box const &tbx, Box const &bxg2, Array4< Real const > const &s, Array4< EBCellFlag > const &cell, Array4< Type_t > const &fx, Array4< Type_t > const &fy)
Definition AMReX_EB2_2D_C.H:9
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int check_mvmc(int i, int j, int, Array4< Real const > const &fine)
Definition AMReX_EB2_2D_C.H:77
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 min_ubound(BoxND< dim > const &b1, BoxND< dim > const &b2) noexcept
Definition AMReX_Box.H:1928
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE BoxND< dim > surroundingNodes(const BoxND< dim > &b, int dir) noexcept
Returns a BoxND with NODE based coordinates in direction dir that encloses BoxND b....
Definition AMReX_Box.H:1399
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE BoxND< dim > grow(const BoxND< dim > &b, int i) noexcept
Grow BoxND in all directions by given amount.
Definition AMReX_Box.H:1211
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 max_lbound(BoxND< dim > const &b1, BoxND< dim > const &b2) noexcept
Definition AMReX_Box.H:1909
AMREX_GPU_HOST_DEVICE AMREX_ATTRIBUTE_FLATTEN_FOR void Loop(Dim3 lo, Dim3 hi, F const &f) noexcept
Definition AMReX_Loop.H:125
Definition AMReX_FabArrayCommI.H:896
Definition AMReX_Array4.H:61