28 const IntVect vdim(IntVect::TheDimensionVector(a_idim));
30 const BoxArray& grids = amrex::convert(a_grids, vdim);
32 m_cellflags =
new FabArray<EBCellFlagFab>(grids, a_dmap, 1, a_ngrow[0], MFInfo(),
33 DefaultFabFactory<EBCellFlagFab>());
37 for (MFIter mfi(*
m_cellflags,
false); mfi.isValid(); ++mfi) {
38 auto& fab = (*m_cellflags)[mfi];
39 fab.setType(FabType::singlevalued);
42 m_volfrac =
new MultiFab(grids, a_dmap, 1, a_ngrow[1], MFInfo(), FArrayBoxFactory());
43 m_volcent =
new MultiFab(grids, a_dmap, AMREX_SPACEDIM, a_ngrow[2], MFInfo(), FArrayBoxFactory());
45 for (
int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
46 const BoxArray& faceba = amrex::convert(a_grids, IntVect::TheDimensionVector(idim));
47 m_areafrac[idim] =
new MultiFab(faceba, a_dmap, 1, a_ngrow[2], MFInfo(), FArrayBoxFactory());
48 m_facecent[idim] =
new MultiFab(faceba, a_dmap, AMREX_SPACEDIM-1, a_ngrow[2], MFInfo(), FArrayBoxFactory());
51 m_bndryarea =
new MultiFab(grids, a_dmap, 1, a_ngrow[2], MFInfo(), FArrayBoxFactory());
52 m_bndrycent =
new MultiFab(grids, a_dmap, AMREX_SPACEDIM, a_ngrow[2], MFInfo(), FArrayBoxFactory());
53 m_bndrynorm =
new MultiFab(grids, a_dmap, AMREX_SPACEDIM, a_ngrow[2], MFInfo(), FArrayBoxFactory());
58 for (
int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
66 const auto& FlagFab = a_factory->getMultiEBCellFlagFab();
68 for (MFIter mfi(*
m_cellflags,
false); mfi.isValid(); ++mfi) {
70 const Box& bx = mfi.validbox();
71 const Box domain = surroundingNodes(a_geom.Domain(), a_idim);
73 if (FlagFab[mfi].getType(bx) == FabType::singlevalued ) {
75 GpuArray<Real, AMREX_SPACEDIM> dx = a_geom.CellSizeArray();
77 Array4<EBCellFlag const>
const& flag = FlagFab.const_array(mfi);
84 Array4<Real const>
const& bnorm = a_factory->getBndryNormal()[mfi].const_array();
85 Array4<Real const>
const& bcent = a_factory->getBndryCent()[mfi].const_array();
88 Array4<EBCellFlag>
const& aux_flag =
m_cellflags->array(mfi);
89 Array4<Real>
const& aux_vfrac =
m_volfrac->array(mfi);
90 Array4<Real>
const& aux_vcent =
m_volcent->array(mfi);
92 Array4<Real>
const& aux_afrac_x =
m_areafrac[0]->array(mfi);
93 Array4<Real>
const& aux_afrac_y =
m_areafrac[1]->array(mfi);
94 Array4<Real>
const& aux_afrac_z =
m_areafrac[2]->array(mfi);
96 Array4<Real>
const& aux_fcent_x =
m_facecent[0]->array(mfi);
97 Array4<Real>
const& aux_fcent_y =
m_facecent[1]->array(mfi);
98 Array4<Real>
const& aux_fcent_z =
m_facecent[2]->array(mfi);
100 Array4<Real>
const& aux_barea =
m_bndryarea->array(mfi);
101 Array4<Real>
const& aux_bcent =
m_bndrycent->array(mfi);
102 Array4<Real>
const& aux_bnorm =
m_bndrynorm->array(mfi);
104 bool is_per = a_geom.isPeriodic(a_idim);
107 #ifndef AMREX_USE_GPU
110 dx, bx, domain, bnorm, bcent, flag,
111 aux_flag, aux_vfrac, aux_vcent,
112 aux_afrac_x, aux_afrac_y, aux_afrac_z,
113 aux_fcent_x, aux_fcent_y, aux_fcent_z,
114 aux_barea, aux_bcent, aux_bnorm,
115 vdim, idim=a_idim, is_per ]
116 AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept
121 aux_flag(i,j,k).setCovered();
122 aux_flag(i,j,k).setDisconnected();
124 aux_vfrac(i,j,k) = 0.0;
125 aux_vcent(i,j,k,0) = 0.0;
126 aux_vcent(i,j,k,1) = 0.0;
127 aux_vcent(i,j,k,2) = 0.0;
129 aux_afrac_x(i,j,k) = 0.0;
130 aux_afrac_y(i,j,k) = 0.0;
131 aux_afrac_z(i,j,k) = 0.0;
133 aux_fcent_x(i,j,k,0) = 0.0; aux_fcent_x(i,j,k,1) = 0.0;
134 aux_fcent_y(i,j,k,0) = 0.0; aux_fcent_y(i,j,k,1) = 0.0;
135 aux_fcent_z(i,j,k,0) = 0.0; aux_fcent_z(i,j,k,1) = 0.0;
137 if (i==bx.bigEnd(0)) {
138 aux_afrac_x(i+1,j,k) = 0.0;
139 aux_fcent_x(i+1,j,k,0) = 0.0;
140 aux_fcent_x(i+1,j,k,1) = 0.0;
142 if (j==bx.bigEnd(1)) {
143 aux_afrac_y(i,j+1,k) = 0.0;
144 aux_fcent_y(i,j+1,k,0) = 0.0;
145 aux_fcent_y(i,j+1,k,1) = 0.0;
147 if (k==bx.bigEnd(2)) {
148 aux_afrac_z(i,j,k+1) = 0.0;
149 aux_fcent_z(i,j,k+1,0) = 0.0;
150 aux_fcent_z(i,j,k+1,1) = 0.0;
153 aux_barea(i,j,k) = 0.0;
155 aux_bcent(i,j,k,0) = 0.0;
156 aux_bcent(i,j,k,1) = 0.0;
157 aux_bcent(i,j,k,2) = 0.0;
159 aux_bnorm(i,j,k,0) = 0.0;
160 aux_bnorm(i,j,k,1) = 0.0;
161 aux_bnorm(i,j,k,2) = 0.0;
164 IntVect iv_hi(i,j,k);
165 IntVect iv_lo(iv_hi - vdim);
166 if (!is_per && iv_hi[idim]==domain.bigEnd(idim)){
169 if (!is_per && iv_hi[idim]==domain.smallEnd(idim)){
173 if ( flag(iv_lo).isCovered() && flag(iv_hi).isCovered()) {
177 }
else if ( flag(iv_lo).isRegular() && flag(iv_hi).isRegular()) {
179 aux_flag(i,j,k).setRegular();
181 aux_vfrac(i,j,k) = 1.0;
183 aux_afrac_x(i,j,k) = 1.0;
184 aux_afrac_y(i,j,k) = 1.0;
185 aux_afrac_z(i,j,k) = 1.0;
187 aux_fcent_x(i,j,k,0) = 0.0; aux_fcent_x(i,j,k,1) = 0.0;
188 aux_fcent_y(i,j,k,0) = 0.0; aux_fcent_y(i,j,k,1) = 0.0;
189 aux_fcent_z(i,j,k,0) = 0.0; aux_fcent_z(i,j,k,1) = 0.0;
191 if (i==bx.bigEnd(0)) {
192 aux_afrac_x(i+1,j,k) = 1.0;
193 aux_fcent_x(i+1,j,k,0) = 0.0;
194 aux_fcent_x(i+1,j,k,1) = 0.0;
196 if (j==bx.bigEnd(1)) {
197 aux_afrac_y(i,j+1,k) = 1.0;
198 aux_fcent_y(i,j+1,k,0) = 0.0;
199 aux_fcent_y(i,j+1,k,1) = 0.0;
201 if (k==bx.bigEnd(2)) {
202 aux_afrac_z(i,j,k+1) = 1.0;
203 aux_fcent_z(i,j,k+1,0) = 0.0;
204 aux_fcent_z(i,j,k+1,1) = 0.0;
209 #ifndef AMREX_USE_GPU
210 if (verbose) { Print() <<
"\ncell: " << amrex::IntVect(i,j,k) <<
"\n"; }
212 Array<Real,AMREX_SPACEDIM> lo_arr = {-0.5,-0.5,-0.5};
213 Array<Real,AMREX_SPACEDIM> hi_arr = { 0.5, 0.5, 0.5};
222 RealVect lo_point (bcent(iv_lo,0), bcent(iv_lo,1), bcent(iv_lo,2));
223 RealVect lo_normal(bnorm(iv_lo,0), bnorm(iv_lo,1), bnorm(iv_lo,2));
225 if (!is_per && iv_hi[idim]==domain.smallEnd(idim)){
226 lo_point[idim] += 1.0;
229 if (flag(iv_lo).isSingleValued() ) {
231 Real bnorm_x = bnorm(iv_lo,0) * dx[0];
232 Real bnorm_y = bnorm(iv_lo,1) * dx[1];
233 Real bnorm_z = bnorm(iv_lo,2) * dx[2];
235 Real norm = sqrt( bnorm_x*bnorm_x + bnorm_y*bnorm_y + bnorm_z*bnorm_z);
237 RealVect bnorm_isoparam ( bnorm_x / norm, bnorm_y / norm, bnorm_z / norm);
241 lo_normal = bnorm_isoparam;
248 RealBox lo_rbx(lo_arr.data(), hi_arr.data());
250 eb_cut_cell_ lo_eb_cc(flag(iv_lo), lo_rbx, lo_point, lo_normal);
254 AMREX_ASSERT( !flag(iv_lo).isCovered() || lo_eb_cc.isCovered() );
255 AMREX_ASSERT( !flag(iv_lo).isRegular() || lo_eb_cc.isRegular() );
261 RealVect hi_point (bcent(iv_hi,0), bcent(iv_hi,1), bcent(iv_hi,2));
262 RealVect hi_normal(bnorm(iv_hi,0), bnorm(iv_hi,1), bnorm(iv_hi,2));
264 if (!is_per && iv_hi[idim]==domain.bigEnd(idim)){
265 lo_point[idim] += -1.0;
268 if (flag(iv_hi).isSingleValued() ) {
270 Real bnorm_x = bnorm(iv_hi,0) * dx[0];
271 Real bnorm_y = bnorm(iv_hi,1) * dx[1];
272 Real bnorm_z = bnorm(iv_hi,2) * dx[2];
274 Real norm = sqrt( bnorm_x*bnorm_x + bnorm_y*bnorm_y + bnorm_z*bnorm_z);
276 RealVect bnorm_isoparam ( bnorm_x / norm, bnorm_y / norm, bnorm_z / norm);
280 hi_normal = bnorm_isoparam;
287 RealBox hi_rbx(lo_arr.data(), hi_arr.data());
289 eb_cut_cell_ hi_eb_cc(flag(iv_hi), hi_rbx, hi_point, hi_normal);
293 AMREX_ASSERT( !flag(iv_hi).isCovered() || hi_eb_cc.isCovered() );
294 AMREX_ASSERT( !flag(iv_hi).isRegular() || hi_eb_cc.isRegular() );
297 #if defined(AMREX_DEBUG) || defined(AMREX_TESTING) || 1
307 eb_cut_cell_ hi_hi_eb_cc(flag(iv_hi), lo_rbx, hi_point, hi_normal);
311 #ifndef AMREX_USE_GPU
312 if ( !(!flag(iv_hi).isRegular() || hi_hi_eb_cc.isRegular()) ||
313 !(!flag(iv_hi).isCovered() || hi_hi_eb_cc.isCovered()) ) {
314 Print() <<
"flag(iv_hi) and hi_hi_eb_cc flags do not agree\n"
315 <<
"\n isRegular() " << flag(iv_hi).isRegular() <<
" " << hi_hi_eb_cc.isRegular()
316 <<
"\n isCovered() " << flag(iv_hi).isCovered() <<
" " << hi_hi_eb_cc.isCovered()
322 AMREX_ALWAYS_ASSERT( !flag(iv_hi).isRegular() || hi_hi_eb_cc.isRegular() );
323 AMREX_ALWAYS_ASSERT( !flag(iv_hi).isCovered() || hi_hi_eb_cc.isCovered() );
330 if ( flag(iv_hi).isSingleValued() ) {
332 Real
const adx = (idim == 0)
333 ? (hi_eb_cc.areaLo(0) - hi_hi_eb_cc.areaHi(0)) * dx[1] * dx[2]
334 : (hi_eb_cc.areaLo(0) + hi_hi_eb_cc.areaLo(0)) * dx[1] * dx[2]
335 - (hi_eb_cc.areaHi(0) + hi_hi_eb_cc.areaHi(0)) * dx[1] * dx[2];
337 Real
const ady = (idim == 1)
338 ? (hi_eb_cc.areaLo(1) - hi_hi_eb_cc.areaHi(1)) * dx[0] * dx[2]
339 : (hi_eb_cc.areaLo(1) + hi_hi_eb_cc.areaLo(1)) * dx[0] * dx[2]
340 - (hi_eb_cc.areaHi(1) + hi_hi_eb_cc.areaHi(1)) * dx[0] * dx[2];
342 Real
const adz = (idim == 2)
343 ? (hi_eb_cc.areaLo(2) - hi_hi_eb_cc.areaHi(2)) * dx[0] * dx[1]
344 : (hi_eb_cc.areaLo(2) + hi_hi_eb_cc.areaLo(2)) * dx[0] * dx[1]
345 - (hi_eb_cc.areaHi(2) + hi_hi_eb_cc.areaHi(2)) * dx[0] * dx[1];
347 Real
const apnorm = std::sqrt(adx*adx + ady*ady + adz*adz);
350 Real
const apnorminv = 1. / apnorm;
351 RealVect
const normal(adx*apnorminv, ady*apnorminv, adz*apnorminv);
352 Real
const dot_normals = normal.dotProduct(hi_normal);
354 #ifndef AMREX_USE_GPU
355 if ( !amrex::almostEqual(dot_normals, 1.0) ) {
356 Print() <<
"\nFail: check-1 dot_normals " << dot_normals
362 }
else if (verbose) {
363 Print() <<
"Pass: dot_normals = 1.0\n";
367 AMREX_ALWAYS_ASSERT( amrex::almostEqual(dot_normals, 1.0) );
372 #ifndef AMREX_USE_GPU
373 Real
const abs_err = std::abs( hi_eb_cc.areaHi(idim) - hi_hi_eb_cc.areaLo(idim) );
375 if ( abs_err >= machine_tol ) {
376 Print() <<
"\nFail: check-2 area abs_err: " << abs_err
377 <<
"\n hi_eb_cc.areaHi " << hi_eb_cc.areaHi(idim)
378 <<
"\n hi_hi_eb_cc.areaLo " << hi_hi_eb_cc.areaLo(idim)
380 }
else if (verbose) {
381 Print() <<
"Pass: hi_eb_cc.areaHi = hi_hi_eb_cc.areaLo"
382 <<
" abs_err: " << abs_err <<
"\n";
384 AMREX_ALWAYS_ASSERT( abs_err < machine_tol );
389 { Real
const abs_err = amrex::max(std::abs(lo_eb_cc.areaHi(idim) - afrac(iv_hi)),
390 std::abs(hi_eb_cc.areaLo(idim) - afrac(iv_hi)));
391 Real compare_tol = 5.0e-6;
392 #ifndef AMREX_USE_GPU
393 if ( abs_err >= compare_tol ) {
395 Print() <<
"\nFail: check-3 area abs_err " << abs_err
396 <<
"\n hi_eb_cc.areaLo(" << idim <<
") = " << hi_eb_cc.areaLo(idim)
397 <<
"\n lo_eb_cc.areaHi(" << idim <<
") = " << lo_eb_cc.areaHi(idim)
398 <<
"\n afrac" << iv_hi <<
" = " << afrac(iv_hi)
400 }
else if (verbose) {
401 Print() <<
"Pass: hi_eb_cc.areaLo = afrac = " << afrac(iv_hi)
402 <<
" abs_err: " << abs_err <<
"\n";
405 AMREX_ALWAYS_ASSERT( abs_err < compare_tol );
410 { Real
const vol = hi_eb_cc.volume() + hi_hi_eb_cc.volume();
411 Real
const abs_err = amrex::Math::abs(vfrac(iv_hi) - vol);
412 Real compare_tol = 5.0e-6;
413 #ifndef AMREX_USE_GPU
414 if ( abs_err >= compare_tol ) {
417 amrex::Print() <<
"\nFail: check-4 volume abs_err: " << abs_err
418 <<
"\n point: " << hi_point
419 <<
"\n normal: " << hi_normal
420 <<
"\n hi_eb_cc.volume() " << hi_eb_cc.volume()
421 <<
"\n hi_hi_eb_cc.volume() " << hi_hi_eb_cc.volume()
422 <<
"\n vfrac: " << vfrac(iv_hi)
424 }
else if (verbose) {
425 Print() <<
"Pass: hi_eb_cc + hi_hi_eb_cc = vfrac = " << vfrac(iv_hi)
426 <<
" abs_err: " << abs_err <<
"\n";
429 AMREX_ALWAYS_ASSERT( abs_err < compare_tol );
439 if (lo_eb_cc.isCovered() && hi_eb_cc.isCovered()) {
443 }
else if (lo_eb_cc.isRegular() && hi_eb_cc.isRegular()) {
445 aux_flag(i,j,k).setRegular();
447 aux_vfrac(i,j,k) = 1.0;
449 aux_afrac_x(i,j,k) = 1.0;
450 aux_afrac_y(i,j,k) = 1.0;
451 aux_afrac_z(i,j,k) = 1.0;
453 aux_fcent_x(i,j,k,0) = 0.0; aux_fcent_x(i,j,k,1) = 0.0;
454 aux_fcent_y(i,j,k,0) = 0.0; aux_fcent_y(i,j,k,1) = 0.0;
455 aux_fcent_z(i,j,k,0) = 0.0; aux_fcent_z(i,j,k,1) = 0.0;
457 if (i==bx.bigEnd(0)) {
458 aux_afrac_x(i+1,j,k) = 1.0;
459 aux_fcent_x(i+1,j,k,0) = 0.0; aux_fcent_x(i+1,j,k,1) = 0.0;
461 if (j==bx.bigEnd(1)) {
462 aux_afrac_y(i,j+1,k) = 1.0;
463 aux_fcent_y(i,j+1,k,0) = 0.0; aux_fcent_y(i,j+1,k,1) = 0.0;
465 if (k==bx.bigEnd(2)) {
466 aux_afrac_z(i,j,k+1) = 1.0;
467 aux_fcent_z(i,j,k+1,0) = 0.0; aux_fcent_z(i,j,k+1,1) = 0.0;
470 }
else if ( (lo_eb_cc.isRegular() && hi_eb_cc.isCovered())
471 || (lo_eb_cc.isCovered() && hi_eb_cc.isRegular()) ) {
474 #ifndef AMREX_USE_GPU
475 Print()<<
"eb_aux_ / Check: Regular and Covered cut cells are facing each other." << std::endl;
480 aux_flag(i,j,k).setSingleValued();
484 Real lo_vol {lo_eb_cc.volume()};
485 Real hi_vol {hi_eb_cc.volume()};
487 aux_vfrac(i,j,k) = lo_vol + hi_vol;
499 RealVect lo_vcent {lo_eb_cc.centVol()};
500 RealVect hi_vcent {hi_eb_cc.centVol()};
502 lo_vcent[idim] = lo_vcent[idim] - 0.5;
503 hi_vcent[idim] = hi_vcent[idim] + 0.5;
505 aux_vcent(i,j,k,0) = ( lo_vol * lo_vcent[0] + hi_vol * hi_vcent[0] ) / aux_vfrac(i,j,k);
506 aux_vcent(i,j,k,1) = ( lo_vol * lo_vcent[1] + hi_vol * hi_vcent[1] ) / aux_vfrac(i,j,k);
507 aux_vcent(i,j,k,2) = ( lo_vol * lo_vcent[2] + hi_vol * hi_vcent[2] ) / aux_vfrac(i,j,k);
511 Real lo_areaLo_x {lo_eb_cc.areaLo(0)};
512 Real lo_areaLo_y {lo_eb_cc.areaLo(1)};
513 Real lo_areaLo_z {lo_eb_cc.areaLo(2)};
515 Real hi_areaLo_x {hi_eb_cc.areaLo(0)};
516 Real hi_areaLo_y {hi_eb_cc.areaLo(1)};
517 Real hi_areaLo_z {hi_eb_cc.areaLo(2)};
519 aux_afrac_x(i,j,k) = (idim == 0) ? lo_areaLo_x : lo_areaLo_x + hi_areaLo_x;
520 aux_afrac_y(i,j,k) = (idim == 1) ? lo_areaLo_y : lo_areaLo_y + hi_areaLo_y;
521 aux_afrac_z(i,j,k) = (idim == 2) ? lo_areaLo_z : lo_areaLo_z + hi_areaLo_z;
523 if (i==bx.bigEnd(0)) {
524 Real lo_areaHi_x {lo_eb_cc.areaHi(0)};
525 Real hi_areaHi_x {hi_eb_cc.areaHi(0)};
526 aux_afrac_x(i+1,j,k) = (idim == 0) ? hi_areaHi_x : lo_areaHi_x + hi_areaHi_x;
528 if (j==bx.bigEnd(1)) {
529 Real lo_areaHi_y {lo_eb_cc.areaHi(1)};
530 Real hi_areaHi_y {hi_eb_cc.areaHi(1)};
531 aux_afrac_y(i,j+1,k) = (idim == 1) ? hi_areaHi_y : lo_areaHi_y + hi_areaHi_y;
533 if (k==bx.bigEnd(2)) {
534 Real lo_areaHi_z {lo_eb_cc.areaHi(2)};
535 Real hi_areaHi_z {hi_eb_cc.areaHi(2)};
536 aux_afrac_z(i,j,k+1) = (idim == 2) ? hi_areaHi_z : lo_areaHi_z + hi_areaHi_z;
549 RealVect lo_centLo_x {lo_eb_cc.centLo(0)};
550 RealVect lo_centLo_y {lo_eb_cc.centLo(1)};
551 RealVect lo_centLo_z {lo_eb_cc.centLo(2)};
553 RealVect hi_centLo_x {hi_eb_cc.centLo(0)};
554 RealVect hi_centLo_y {hi_eb_cc.centLo(1)};
555 RealVect hi_centLo_z {hi_eb_cc.centLo(2)};
558 aux_fcent_x(i,j,k,0) = lo_centLo_x[1];
559 aux_fcent_x(i,j,k,1) = lo_centLo_x[2];
560 aux_fcent_y(i,j,k,0) = (aux_afrac_y(i,j,k) > 0.0)
561 ? ( lo_areaLo_y * (lo_centLo_y[0] - 0.5)
562 + hi_areaLo_y * (hi_centLo_y[0] + 0.5) ) / aux_afrac_y(i,j,k)
564 aux_fcent_y(i,j,k,1) = (aux_afrac_y(i,j,k) > 0.0)
565 ? ( lo_areaLo_y * lo_centLo_y[2]
566 + hi_areaLo_y * hi_centLo_y[2] ) / aux_afrac_y(i,j,k)
568 aux_fcent_z(i,j,k,0) = (aux_afrac_z(i,j,k) > 0.0)
569 ? ( lo_areaLo_z * (lo_centLo_z[0] - 0.5)
570 + hi_areaLo_z * (hi_centLo_z[0] + 0.5) ) / aux_afrac_z(i,j,k)
572 aux_fcent_z(i,j,k,1) = (aux_afrac_z(i,j,k) > 0.0)
573 ? ( lo_areaLo_z * lo_centLo_z[1]
574 + hi_areaLo_z * hi_centLo_z[1] ) / aux_afrac_z(i,j,k)
576 }
else if (idim == 1) {
577 aux_fcent_x(i,j,k,0) = (aux_afrac_x(i,j,k) > 0.0)
578 ? ( lo_areaLo_x * (lo_centLo_x[1] - 0.5)
579 + hi_areaLo_x * (hi_centLo_x[1] + 0.5) ) / aux_afrac_x(i,j,k)
581 aux_fcent_x(i,j,k,1) = (aux_afrac_x(i,j,k) > 0.0)
582 ? ( lo_areaLo_x * lo_centLo_x[2]
583 + hi_areaLo_x * hi_centLo_x[2] ) / aux_afrac_x(i,j,k)
585 aux_fcent_y(i,j,k,0) = lo_centLo_y[0];
586 aux_fcent_y(i,j,k,1) = lo_centLo_y[2];
587 aux_fcent_z(i,j,k,0) = (aux_afrac_z(i,j,k) > 0.0)
588 ? ( lo_areaLo_z * lo_centLo_z[0]
589 + hi_areaLo_z * hi_centLo_z[0] ) / aux_afrac_z(i,j,k)
591 aux_fcent_z(i,j,k,1) = (aux_afrac_z(i,j,k) > 0.0)
592 ? ( lo_areaLo_z * (lo_centLo_z[1] - 0.5)
593 + hi_areaLo_z * (hi_centLo_z[1] + 0.5) ) / aux_afrac_z(i,j,k)
595 }
else if (idim == 2) {
596 aux_fcent_x(i,j,k,0) = (aux_afrac_x(i,j,k) > 0.0)
597 ? ( lo_areaLo_x * lo_centLo_x[1]
598 + hi_areaLo_x * hi_centLo_x[1] ) / aux_afrac_x(i,j,k)
600 aux_fcent_x(i,j,k,1) = (aux_afrac_x(i,j,k) > 0.0)
601 ? ( lo_areaLo_x * (lo_centLo_x[2] - 0.5)
602 + hi_areaLo_x * (hi_centLo_x[2] + 0.5) ) / aux_afrac_x(i,j,k)
604 aux_fcent_y(i,j,k,0) = (aux_afrac_y(i,j,k) > 0.0)
605 ? ( lo_areaLo_y * lo_centLo_y[0]
606 + hi_areaLo_y * hi_centLo_y[0] ) / aux_afrac_y(i,j,k)
608 aux_fcent_y(i,j,k,1) = (aux_afrac_y(i,j,k) > 0.0)
609 ? ( lo_areaLo_y * (lo_centLo_y[2] - 0.5)
610 + hi_areaLo_y * (hi_centLo_y[2] + 0.5) ) / aux_afrac_y(i,j,k)
612 aux_fcent_z(i,j,k,0) = lo_centLo_z[0];
613 aux_fcent_z(i,j,k,1) = lo_centLo_z[1];
616 if (i==bx.bigEnd(0)) {
617 Real lo_areaHi_x {lo_eb_cc.areaHi(0)};
618 Real hi_areaHi_x {hi_eb_cc.areaHi(0)};
619 RealVect lo_centHi_x {lo_eb_cc.centHi(0)};
620 RealVect hi_centHi_x {hi_eb_cc.centHi(0)};
622 aux_fcent_x(i+1,j,k,0) = hi_centHi_x[1];
623 aux_fcent_x(i+1,j,k,1) = hi_centHi_x[2];
624 }
else if (idim == 1) {
625 aux_fcent_x(i+1,j,k,0) = (aux_afrac_x(i+1,j,k) > 0.0)
626 ? ( lo_areaHi_x * (lo_centHi_x[1] - 0.5)
627 + hi_areaHi_x * (hi_centHi_x[1] + 0.5) ) / aux_afrac_x(i+1,j,k)
629 aux_fcent_x(i+1,j,k,1) = (aux_afrac_x(i+1,j,k) > 0.0)
630 ? ( lo_areaHi_x * lo_centHi_x[2]
631 + hi_areaHi_x * hi_centHi_x[2] ) / aux_afrac_x(i+1,j,k)
633 }
else if (idim == 2) {
634 aux_fcent_x(i+1,j,k,0) = (aux_afrac_x(i+1,j,k) > 0.0)
635 ? ( lo_areaHi_x * lo_centHi_x[1]
636 + hi_areaHi_x * hi_centHi_x[1] ) / aux_afrac_x(i+1,j,k)
638 aux_fcent_x(i+1,j,k,1) = (aux_afrac_x(i+1,j,k) > 0.0)
639 ? ( lo_areaHi_x * (lo_centHi_x[2] - 0.5)
640 + hi_areaHi_x * (hi_centHi_x[2] + 0.5) ) / aux_afrac_x(i+1,j,k)
644 if (j==bx.bigEnd(1)) {
645 Real lo_areaHi_y {lo_eb_cc.areaHi(1)};
646 Real hi_areaHi_y {hi_eb_cc.areaHi(1)};
647 RealVect lo_centHi_y {lo_eb_cc.centHi(1)};
648 RealVect hi_centHi_y {hi_eb_cc.centHi(1)};
650 aux_fcent_y(i,j+1,k,0) = (aux_afrac_y(i,j+1,k) > 0.0)
651 ? ( lo_areaHi_y * (lo_centHi_y[0] - 0.5)
652 + hi_areaHi_y * (hi_centHi_y[0] + 0.5) ) / aux_afrac_y(i,j+1,k)
654 aux_fcent_y(i,j+1,k,1) = (aux_afrac_y(i,j+1,k) > 0.0)
655 ? ( lo_areaHi_y * lo_centHi_y[2]
656 + hi_areaHi_y * hi_centHi_y[2] ) / aux_afrac_y(i,j+1,k)
658 }
else if (idim == 1) {
659 aux_fcent_y(i,j+1,k,0) = lo_centHi_y[0];
660 aux_fcent_y(i,j+1,k,1) = lo_centHi_y[2];
661 }
else if (idim == 2) {
662 aux_fcent_y(i,j+1,k,0) = (aux_afrac_y(i,j+1,k) > 0.0)
663 ? ( lo_areaHi_y * lo_centHi_y[0]
664 + hi_areaHi_y * hi_centHi_y[0] ) / aux_afrac_y(i,j+1,k)
666 aux_fcent_y(i,j+1,k,1) = (aux_afrac_y(i,j+1,k) > 0.0)
667 ? ( lo_areaHi_y * (lo_centHi_y[2] - 0.5)
668 + hi_areaHi_y * (hi_centHi_y[2] + 0.5) ) / aux_afrac_y(i,j+1,k)
672 if (k==bx.bigEnd(2)) {
673 Real lo_areaHi_z {lo_eb_cc.areaHi(2)};
674 Real hi_areaHi_z {hi_eb_cc.areaHi(2)};
675 RealVect lo_centHi_z {lo_eb_cc.centHi(2)};
676 RealVect hi_centHi_z {hi_eb_cc.centHi(2)};
678 aux_fcent_z(i,j,k+1,0) = (aux_afrac_z(i,j,k+1) > 0.0)
679 ? ( lo_areaHi_z * (lo_centHi_z[0] - 0.5)
680 + hi_areaHi_z * (hi_centHi_z[0] + 0.5) ) / aux_afrac_z(i,j,k+1)
682 aux_fcent_z(i,j,k+1,1) = (aux_afrac_z(i,j,k+1) > 0.0)
683 ? ( lo_areaHi_z * lo_centHi_z[1]
684 + hi_areaHi_z * hi_centHi_z[1] ) / aux_afrac_z(i,j,k+1)
686 }
else if (idim == 1) {
687 aux_fcent_z(i,j,k+1,0) = (aux_afrac_z(i,j,k+1) > 0.0)
688 ? ( lo_areaHi_z * lo_centHi_z[0]
689 + hi_areaHi_z * hi_centHi_z[0] ) / aux_afrac_z(i,j,k+1)
691 aux_fcent_z(i,j,k+1,1) = (aux_afrac_z(i,j,k+1) > 0.0)
692 ? ( lo_areaHi_z * (lo_centHi_z[1] - 0.5)
693 + hi_areaHi_z * (hi_centHi_z[1] + 0.5) ) / aux_afrac_z(i,j,k+1)
695 }
else if (idim == 2) {
696 aux_fcent_z(i,j,k+1,0) = lo_centHi_z[0];
697 aux_fcent_z(i,j,k+1,1) = lo_centHi_z[1];
703 Real lo_areaBoun {lo_eb_cc.areaBoun()};
704 Real hi_areaBoun {hi_eb_cc.areaBoun()};
706 aux_barea(i,j,k) = lo_areaBoun + hi_areaBoun;
710 RealVect lo_centBoun {lo_eb_cc.centBoun()};
711 RealVect hi_centBoun {hi_eb_cc.centBoun()};
714 aux_bcent(i,j,k,0) = ( lo_areaBoun * (lo_centBoun[0]-0.5) + hi_areaBoun * (hi_centBoun[0]+0.5) ) / aux_barea(i,j,k);
715 aux_bcent(i,j,k,1) = ( lo_areaBoun * lo_centBoun[1] + hi_areaBoun * hi_centBoun[1] ) / aux_barea(i,j,k);
716 aux_bcent(i,j,k,2) = ( lo_areaBoun * lo_centBoun[2] + hi_areaBoun * hi_centBoun[2] ) / aux_barea(i,j,k);
717 }
else if (idim == 1) {
718 aux_bcent(i,j,k,0) = ( lo_areaBoun * lo_centBoun[0] + hi_areaBoun * hi_centBoun[0] ) / aux_barea(i,j,k);
719 aux_bcent(i,j,k,1) = ( lo_areaBoun * (lo_centBoun[1]-0.5) + hi_areaBoun * (hi_centBoun[1]+0.5) ) / aux_barea(i,j,k);
720 aux_bcent(i,j,k,2) = ( lo_areaBoun * lo_centBoun[2] + hi_areaBoun * hi_centBoun[2] ) / aux_barea(i,j,k);
721 }
else if (idim == 2) {
722 aux_bcent(i,j,k,0) = ( lo_areaBoun * lo_centBoun[0] + hi_areaBoun * hi_centBoun[0] ) / aux_barea(i,j,k);
723 aux_bcent(i,j,k,1) = ( lo_areaBoun * lo_centBoun[1] + hi_areaBoun * hi_centBoun[1] ) / aux_barea(i,j,k);
724 aux_bcent(i,j,k,2) = ( lo_areaBoun * (lo_centBoun[2]-0.5) + hi_areaBoun * (hi_centBoun[2]+0.5) ) / aux_barea(i,j,k);
729 RealVect eb_normal = ( lo_areaBoun * lo_normal + hi_areaBoun * hi_normal )/ aux_barea(i,j,k);
731 aux_bnorm(i,j,k,0) = eb_normal[0];
732 aux_bnorm(i,j,k,1) = eb_normal[1];
733 aux_bnorm(i,j,k,2) = eb_normal[2];
742 ParallelFor(bx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept
744 EB2::build_cellflag_from_ap (i, j, k, aux_flag, aux_afrac_x, aux_afrac_y, aux_afrac_z);
753 m_volfrac->FillBoundary(a_geom.periodicity());
754 m_volcent->FillBoundary(a_geom.periodicity());
755 for (
int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
756 m_areafrac[idim]->FillBoundary(a_geom.periodicity());
757 m_facecent[idim]->FillBoundary(a_geom.periodicity());
amrex::FabArray< amrex::EBCellFlagFab > * m_cellflags
Definition: ERF_EBAux.H:43
amrex::MultiFab * m_bndrynorm
Definition: ERF_EBAux.H:48
amrex::MultiFab * m_bndryarea
Definition: ERF_EBAux.H:46
amrex::MultiFab * m_volfrac
Definition: ERF_EBAux.H:44
amrex::MultiFab * m_bndrycent
Definition: ERF_EBAux.H:47
amrex::Array< amrex::MultiFab *, AMREX_SPACEDIM > m_facecent
Definition: ERF_EBAux.H:51
amrex::MultiFab * m_volcent
Definition: ERF_EBAux.H:45
amrex::Array< amrex::MultiFab *, AMREX_SPACEDIM > m_areafrac
Definition: ERF_EBAux.H:50
Definition: ERF_EBCutCell.H:14
real(c_double), parameter epsilon
Definition: ERF_module_model_constants.F90:12