ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Public Member Functions | Private Attributes | List of all members
eb_aux_ Class Reference

#include <ERF_EBAux.H>

Collaboration diagram for eb_aux_:
Collaboration graph
[legend]

Public Member Functions

 eb_aux_ ()
 
 ~eb_aux_ ()
 
void define (int const &a_idim, amrex::Geometry const &a_geom, amrex::BoxArray const &a_grids, amrex::DistributionMapping const &a_dmap, amrex::Vector< int > const &a_ngrow, amrex::EBFArrayBoxFactory const *a_factory)
 
void set_verbose ()
 
const amrex::MultiFab & getVolFrac () const
 
const amrex::MultiCutFab & getVolCent () const
 
const amrex::MultiCutFab & getBndryArea () const
 
const amrex::MultiCutFab & getBndryCent () const
 
const amrex::MultiCutFab & getBndryNorm () const
 
amrex::Array< const amrex::MultiCutFab *, AMREX_SPACEDIM > getAreaFrac () const
 
amrex::Array< const amrex::MultiCutFab *, AMREX_SPACEDIM > getFaceCent () const
 

Private Attributes

int m_verbose
 
amrex::FabArray< amrex::EBCellFlagFab > * m_cellflags = nullptr
 
amrex::MultiFab * m_volfrac = nullptr
 
amrex::MultiCutFab * m_volcent = nullptr
 
amrex::MultiCutFab * m_bndryarea = nullptr
 
amrex::MultiCutFab * m_bndrycent = nullptr
 
amrex::MultiCutFab * m_bndrynorm = nullptr
 
amrex::Array< amrex::MultiCutFab *, AMREX_SPACEDIM > m_areafrac {{AMREX_D_DECL(nullptr, nullptr, nullptr)}}
 
amrex::Array< amrex::MultiCutFab *, AMREX_SPACEDIM > m_facecent {{AMREX_D_DECL(nullptr, nullptr, nullptr)}}
 

Constructor & Destructor Documentation

◆ eb_aux_()

eb_aux_::eb_aux_ ( )
13  : m_verbose(0)
14 // ,m_defined(0)
15 {}
eb_aux_::m_verbose
int m_verbose
Definition: ERF_EBAux.H:38

◆ ~eb_aux_()

eb_aux_::~eb_aux_ ( )
8 {
9 }

Member Function Documentation

◆ define()

void eb_aux_::define ( int const &  a_idim,
amrex::Geometry const &  a_geom,
amrex::BoxArray const &  a_grids,
amrex::DistributionMapping const &  a_dmap,
amrex::Vector< int > const &  a_ngrow,
amrex::EBFArrayBoxFactory const *  a_factory 
)
25 {
26  // Box dbox(a_geom.Domain());
27 
28  const IntVect vdim(IntVect::TheDimensionVector(a_idim));
29 
30  const BoxArray& grids = amrex::convert(a_grids, vdim);
31 
32  m_cellflags = new FabArray<EBCellFlagFab>(grids, a_dmap, 1, a_ngrow[0], MFInfo(),
33  DefaultFabFactory<EBCellFlagFab>());
34 
35  // Set m_cellflags type to singlevalued
36  m_cellflags->setVal(EBCellFlag::TheDefaultCell());
37  for (MFIter mfi(*m_cellflags, false); mfi.isValid(); ++mfi) {
38  auto& fab = (*m_cellflags)[mfi];
39  fab.setType(FabType::singlevalued);
40  }
41 
42  m_volfrac = new MultiFab(grids, a_dmap, 1, a_ngrow[1], MFInfo(), FArrayBoxFactory());
43  m_volcent = new MultiCutFab(grids, a_dmap, AMREX_SPACEDIM, a_ngrow[2], *m_cellflags);
44 
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 MultiCutFab(faceba, a_dmap, 1, a_ngrow[2], *m_cellflags);
48  m_facecent[idim] = new MultiCutFab(faceba, a_dmap, AMREX_SPACEDIM-1, a_ngrow[2], *m_cellflags);
49  }
50 
51  m_bndryarea = new MultiCutFab(grids, a_dmap, 1, a_ngrow[2], *m_cellflags);
52  m_bndrycent = new MultiCutFab(grids, a_dmap, AMREX_SPACEDIM, a_ngrow[2], *m_cellflags);
53  m_bndrynorm = new MultiCutFab(grids, a_dmap, AMREX_SPACEDIM, a_ngrow[2], *m_cellflags);
54 
55  const auto& FlagFab = a_factory->getMultiEBCellFlagFab(); // EBFArrayBoxFactory, EBDataCollection
56 
57  for (MFIter mfi(*m_cellflags, false); mfi.isValid(); ++mfi) {
58 
59  const Box& bx = mfi.validbox();
60 
61  if (FlagFab[mfi].getType(bx) == FabType::singlevalued ) {
62 
63  GpuArray<Real, AMREX_SPACEDIM> dx = a_geom.CellSizeArray();
64 
65  Array4<EBCellFlag const> const& flag = FlagFab.const_array(mfi);
66 
67  // Array4<Real const> const& vfrac = (a_factory->getVolFrac()).const_array(mfi);
68  // Array4<Real const> const& ccent = (a_factory->getCentroid()).const_array(mfi);
69  // Array4<Real const> const& afrac = (a_factory->getAreaFrac()[a_idim])->const_array(mfi);
70 
71  // EB normal and face centroid
72  Array4<Real const> const& bnorm = a_factory->getBndryNormal()[mfi].const_array();
73  Array4<Real const> const& bcent = a_factory->getBndryCent()[mfi].const_array();
74 
75  // aux quantities
76  Array4<EBCellFlag> const& aux_flag = m_cellflags->array(mfi);
77  Array4<Real> const& aux_vfrac = m_volfrac->array(mfi);
78  Array4<Real> const& aux_vcent = m_volcent->array(mfi);
79 
80  Array4<Real> const& aux_afrac_x = m_areafrac[0]->array(mfi);
81  Array4<Real> const& aux_afrac_y = m_areafrac[1]->array(mfi);
82  Array4<Real> const& aux_afrac_z = m_areafrac[2]->array(mfi);
83 
84  Array4<Real> const& aux_fcent_x = m_facecent[0]->array(mfi);
85  Array4<Real> const& aux_fcent_y = m_facecent[1]->array(mfi);
86  Array4<Real> const& aux_fcent_z = m_facecent[2]->array(mfi);
87 
88  Array4<Real> const& aux_barea = m_bndryarea->array(mfi);
89  Array4<Real> const& aux_bcent = m_bndrycent->array(mfi);
90  Array4<Real> const& aux_bnorm = m_bndrynorm->array(mfi);
91 
92  bool is_per = a_geom.isPeriodic(a_idim);
93 
94  ParallelFor(bx, [
95 #ifndef AMREX_USE_GPU
96  verbose=m_verbose,
97 #endif
98  dx, bx, bnorm, bcent, flag,
99  aux_flag, aux_vfrac, aux_vcent,
100  aux_afrac_x, aux_afrac_y, aux_afrac_z,
101  aux_fcent_x, aux_fcent_y, aux_fcent_z,
102  aux_barea, aux_bcent, aux_bnorm,
103  vdim, idim=a_idim, is_per ]
104  AMREX_GPU_DEVICE (int i, int j, int k) noexcept
105  {
106 
107  // defaults to covered and disconnected.
108 
109  aux_flag(i,j,k).setCovered();
110  aux_flag(i,j,k).setDisconnected();
111 
112  aux_vfrac(i,j,k) = 0.0;
113  aux_vcent(i,j,k,0) = 0.0;
114  aux_vcent(i,j,k,1) = 0.0;
115  aux_vcent(i,j,k,2) = 0.0;
116 
117  aux_afrac_x(i,j,k) = 0.0;
118  aux_afrac_y(i,j,k) = 0.0;
119  aux_afrac_z(i,j,k) = 0.0;
120 
121  aux_fcent_x(i,j,k,0) = 0.0; aux_fcent_x(i,j,k,1) = 0.0;
122  aux_fcent_y(i,j,k,0) = 0.0; aux_fcent_y(i,j,k,1) = 0.0;
123  aux_fcent_z(i,j,k,0) = 0.0; aux_fcent_z(i,j,k,1) = 0.0;
124 
125  if (i==bx.bigEnd(0)) {
126  aux_afrac_x(i+1,j,k) = 0.0;
127  aux_fcent_x(i+1,j,k,0) = 0.0; aux_fcent_x(i+1,j,k,1) = 0.0;
128  }
129  if (j==bx.bigEnd(1)) {
130  aux_afrac_y(i,j+1,k) = 0.0;
131  aux_fcent_y(i,j+1,k,0) = 0.0; aux_fcent_y(i,j+1,k,1) = 0.0;
132  }
133  if (k==bx.bigEnd(2)) {
134  aux_afrac_z(i,j,k+1) = 0.0;
135  aux_fcent_z(i,j,k+1,0) = 0.0; aux_fcent_z(i,j,k+1,1) = 0.0;
136  }
137 
138  aux_barea(i,j,k) = 0.0;
139 
140  aux_bcent(i,j,k,0) = 0.0;
141  aux_bcent(i,j,k,1) = 0.0;
142  aux_bcent(i,j,k,2) = 0.0;
143 
144  aux_bnorm(i,j,k,0) = 0.0;
145  aux_bnorm(i,j,k,1) = 0.0;
146  aux_bnorm(i,j,k,2) = 0.0;
147 
148  // Index for low and hi cells
149  IntVect iv_hi(i,j,k);
150  IntVect iv_lo(iv_hi - vdim);
151  if (!is_per && iv_hi[idim]==bx.bigEnd(idim)){
152  iv_hi = iv_lo; // At the upper boundary, hi cell takes the values of the low cell.
153  }
154  if (!is_per && iv_hi[idim]==bx.smallEnd(idim)){
155  iv_lo = iv_hi; // At the lower boundary, low cell takes the values of the high cell.
156  }
157 
158  //
159 
160  if ( flag(iv_lo).isCovered() && flag(iv_hi).isCovered()) {
161 
162  // defaults to covered and disconnected.
163 
164  } else if ( flag(iv_lo).isRegular() && flag(iv_hi).isRegular()) {
165 
166  aux_flag(i,j,k).setRegular();
167  aux_flag(i,j,k).setConnected(vdim);
168 
169  aux_vfrac(i,j,k) = 1.0;
170 
171  aux_afrac_x(i,j,k) = 1.0;
172  aux_afrac_y(i,j,k) = 1.0;
173  aux_afrac_z(i,j,k) = 1.0;
174 
175  aux_fcent_x(i,j,k,0) = 0.0; aux_fcent_x(i,j,k,1) = 0.0;
176  aux_fcent_y(i,j,k,0) = 0.0; aux_fcent_y(i,j,k,1) = 0.0;
177  aux_fcent_z(i,j,k,0) = 0.0; aux_fcent_z(i,j,k,1) = 0.0;
178 
179  if (i==bx.bigEnd(0)) {
180  aux_afrac_x(i+1,j,k) = 1.0;
181  aux_fcent_x(i+1,j,k,0) = 0.0; aux_fcent_x(i+1,j,k,1) = 0.0;
182  }
183  if (j==bx.bigEnd(1)) {
184  aux_afrac_y(i,j+1,k) = 1.0;
185  aux_fcent_y(i,j+1,k,0) = 0.0; aux_fcent_y(i,j+1,k,1) = 0.0;
186  }
187  if (k==bx.bigEnd(2)) {
188  aux_afrac_z(i,j,k+1) = 1.0;
189  aux_fcent_z(i,j,k+1,0) = 0.0; aux_fcent_z(i,j,k+1,1) = 0.0;
190  }
191 
192  } else {
193 
194 #ifndef AMREX_USE_GPU
195  if (verbose) { Print() << "\ncell: " << amrex::IntVect(i,j,k) << "\n"; }
196 #endif
197  Array<Real,AMREX_SPACEDIM> lo_arr = {-0.5,-0.5,-0.5};
198  Array<Real,AMREX_SPACEDIM> hi_arr = { 0.5, 0.5, 0.5};
199 
200  //-----------------------
201  // Low EB cut cell
202  //-----------------------
203 
204  // Map bcent and bnorm to the isoparametric space for anisotropic grids.
205  // (This step is needed because bcent in AMReX is isotropically normalized.)
206 
207  RealVect lo_point (bcent(iv_lo,0), bcent(iv_lo,1), bcent(iv_lo,2));
208  RealVect lo_normal(bnorm(iv_lo,0), bnorm(iv_lo,1), bnorm(iv_lo,2));
209 
210  if (!is_per && iv_hi[idim]==bx.smallEnd(idim)){
211  lo_point[idim] += 1.0; // Move the boundary centroid upward in the idim direction.
212  }
213 
214  if (flag(iv_lo).isSingleValued() ) {
215 
216  Real bnorm_x = bnorm(iv_lo,0) * dx[0];
217  Real bnorm_y = bnorm(iv_lo,1) * dx[1];
218  Real bnorm_z = bnorm(iv_lo,2) * dx[2];
219 
220  Real norm = sqrt( bnorm_x*bnorm_x + bnorm_y*bnorm_y + bnorm_z*bnorm_z);
221 
222  RealVect bnorm_isoparam ( bnorm_x / norm, bnorm_y / norm, bnorm_z / norm);
223 
224  // plane point and normal
225  // lo_point = bcent_isoparam;
226  lo_normal = bnorm_isoparam;
227 
228  }
229 
230  // High side of low cell
231  lo_arr[idim] = 0.0;
232  hi_arr[idim] = 0.5;
233  RealBox lo_rbx(lo_arr.data(), hi_arr.data());
234 
235  eb_cut_cell_ lo_eb_cc(flag(iv_lo), lo_rbx, lo_point, lo_normal);
236 
237  // cell iv_lo covered (regular) imples lo_eb_cc is covered (regular)
238  // The inverse is not always true.
239  AMREX_ASSERT( !flag(iv_lo).isCovered() || lo_eb_cc.isCovered() );
240  AMREX_ASSERT( !flag(iv_lo).isRegular() || lo_eb_cc.isRegular() );
241 
242  //-----------------------
243  // High EB cut cell
244  //-----------------------
245 
246  RealVect hi_point (bcent(iv_hi,0), bcent(iv_hi,1), bcent(iv_hi,2));
247  RealVect hi_normal(bnorm(iv_hi,0), bnorm(iv_hi,1), bnorm(iv_hi,2));
248 
249  if (!is_per && iv_hi[idim]==bx.bigEnd(idim)){
250  lo_point[idim] += -1.0; // Move the boundary centroid downward in the idim direction.
251  }
252 
253  if (flag(iv_hi).isSingleValued() ) {
254 
255  Real bnorm_x = bnorm(iv_hi,0) * dx[0];
256  Real bnorm_y = bnorm(iv_hi,1) * dx[1];
257  Real bnorm_z = bnorm(iv_hi,2) * dx[2];
258 
259  Real norm = sqrt( bnorm_x*bnorm_x + bnorm_y*bnorm_y + bnorm_z*bnorm_z);
260 
261  RealVect bnorm_isoparam ( bnorm_x / norm, bnorm_y / norm, bnorm_z / norm);
262 
263  // plane point and normal
264  // hi_point = bcent_isoparam;
265  hi_normal = bnorm_isoparam;
266 
267  }
268 
269  // Low side of high cell
270  lo_arr[idim] = -0.5;
271  hi_arr[idim] = 0.0;
272  RealBox hi_rbx(lo_arr.data(), hi_arr.data());
273 
274  eb_cut_cell_ hi_eb_cc(flag(iv_hi), hi_rbx, hi_point, hi_normal);
275 
276  // cell iv_hi covered (regular) imples hi_eb_cc is covered (regular)
277  // The inverse is not always true.
278  AMREX_ASSERT( !flag(iv_hi).isCovered() || hi_eb_cc.isCovered() );
279  AMREX_ASSERT( !flag(iv_hi).isRegular() || hi_eb_cc.isRegular() );
280 
281 #if 0
282 #if defined(AMREX_DEBUG) || defined(AMREX_TESTING) || 1
283 
284  { /***************************** SANITY CHECK ***********************\
285  * Perform some basic sanity checks to verify that what we computed *
286  * for cell (i,j,k) compares to what we know to be true. *
287  \******************************************************************/
288 
289  // Compute the cut-cell for the high side of the high cell. This is
290  // only needed for sanity checks.
291 
292  eb_cut_cell_ hi_hi_eb_cc(flag(iv_hi), lo_rbx, hi_point, hi_normal);
293 
294  // cell iv_hi covered (regular) imples hi_hi_eb_cc is covered (regular)
295  // The inverse is not always true.
296 #ifndef AMREX_USE_GPU
297  if ( !(!flag(iv_hi).isRegular() || hi_hi_eb_cc.isRegular()) ||
298  !(!flag(iv_hi).isCovered() || hi_hi_eb_cc.isCovered()) ) {
299  Print() << "flag(iv_hi) and hi_hi_eb_cc flags do not agree\n"
300  << "\n isRegular() " << flag(iv_hi).isRegular() << " " << hi_hi_eb_cc.isRegular()
301  << "\n isCovered() " << flag(iv_hi).isCovered() << " " << hi_hi_eb_cc.isCovered()
302  << "\n";
303  }
304 #endif
305  // If cell iv_hi is regular or covered, then hi_hi_eb_cc must also
306  // be regular or covered. The inverse is not true.
307  AMREX_ALWAYS_ASSERT( !flag(iv_hi).isRegular() || hi_hi_eb_cc.isRegular() );
308  AMREX_ALWAYS_ASSERT( !flag(iv_hi).isCovered() || hi_hi_eb_cc.isCovered() );
309 
310  // The area and volume fractions that are computed for the scalar grid
311  // are slightly different than those we compute from the geometric
312  // reconstruction using the EB point and normal. However, we expect
313  // that the area fractions computed here will give back the same
314  // normal we used to compute them.
315  if ( flag(iv_hi).isSingleValued() ) {
316 
317  Real const adx = (idim == 0)
318  ? (hi_eb_cc.areaLo(0) - hi_hi_eb_cc.areaHi(0)) * dx[1] * dx[2]
319  : (hi_eb_cc.areaLo(0) + hi_hi_eb_cc.areaLo(0)) * dx[1] * dx[2]
320  - (hi_eb_cc.areaHi(0) + hi_hi_eb_cc.areaHi(0)) * dx[1] * dx[2];
321 
322  Real const ady = (idim == 1)
323  ? (hi_eb_cc.areaLo(1) - hi_hi_eb_cc.areaHi(1)) * dx[0] * dx[2]
324  : (hi_eb_cc.areaLo(1) + hi_hi_eb_cc.areaLo(1)) * dx[0] * dx[2]
325  - (hi_eb_cc.areaHi(1) + hi_hi_eb_cc.areaHi(1)) * dx[0] * dx[2];
326 
327  Real const adz = (idim == 2)
328  ? (hi_eb_cc.areaLo(2) - hi_hi_eb_cc.areaHi(2)) * dx[0] * dx[1]
329  : (hi_eb_cc.areaLo(2) + hi_hi_eb_cc.areaLo(2)) * dx[0] * dx[1]
330  - (hi_eb_cc.areaHi(2) + hi_hi_eb_cc.areaHi(2)) * dx[0] * dx[1];
331 
332  Real const apnorm = std::sqrt(adx*adx + ady*ady + adz*adz);
333 
334  // EB normal
335  Real const apnorminv = 1. / apnorm;
336  RealVect const normal(adx*apnorminv, ady*apnorminv, adz*apnorminv);
337  Real const dot_normals = normal.dotProduct(hi_normal);
338 
339 #ifndef AMREX_USE_GPU
340  if ( !amrex::almostEqual(dot_normals, 1.0) ) {
341  Print() << "\nFail: check-1 dot_normals " << dot_normals
342  << '\n';
343 
344  hi_eb_cc.debug();
345  hi_hi_eb_cc.debug();
346 
347  } else if (verbose) {
348  Print() << "Pass: dot_normals = 1.0\n";
349 
350  }
351 #endif
352  AMREX_ALWAYS_ASSERT( amrex::almostEqual(dot_normals, 1.0) );
353  }
354 
355  // The idim area of hi_eb_cc.areaHi() should equal hi_hi_eb_cc.areaLo()
356  {
357 #ifndef AMREX_USE_GPU
358  Real const abs_err = std::abs( hi_eb_cc.areaHi(idim) - hi_hi_eb_cc.areaLo(idim) );
359  Real machine_tol = 10.0*std::numeric_limits<amrex::Real>::epsilon();
360  if ( abs_err >= machine_tol ) {
361  Print() << "\nFail: check-2 area abs_err: " << abs_err
362  << "\n hi_eb_cc.areaHi " << hi_eb_cc.areaHi(idim)
363  << "\n hi_hi_eb_cc.areaLo " << hi_hi_eb_cc.areaLo(idim)
364  << '\n';
365  } else if (verbose) {
366  Print() << "Pass: hi_eb_cc.areaHi = hi_hi_eb_cc.areaLo"
367  << " abs_err: " << abs_err << "\n";
368  }
369  AMREX_ALWAYS_ASSERT( abs_err < machine_tol );
370 #endif
371  }
372 
373  // The low-side area of hi_eb_cc should equal idim afrac.
374  { Real const abs_err = amrex::max(std::abs(lo_eb_cc.areaHi(idim) - afrac(iv_hi)),
375  std::abs(hi_eb_cc.areaLo(idim) - afrac(iv_hi)));
376  Real compare_tol = 5.0e-6;
377 #ifndef AMREX_USE_GPU
378  if ( abs_err >= compare_tol ) {
379  //hi_eb_cc.debug();
380  Print() << "\nFail: check-3 area abs_err " << abs_err
381  << "\n hi_eb_cc.areaLo(" << idim << ") = " << hi_eb_cc.areaLo(idim)
382  << "\n lo_eb_cc.areaHi(" << idim << ") = " << lo_eb_cc.areaHi(idim)
383  << "\n afrac" << iv_hi << " = " << afrac(iv_hi)
384  << '\n';
385  } else if (verbose) {
386  Print() << "Pass: hi_eb_cc.areaLo = afrac = " << afrac(iv_hi)
387  << " abs_err: " << abs_err << "\n";
388  }
389 #endif
390  AMREX_ALWAYS_ASSERT( abs_err < compare_tol );
391  }
392 
393  // The combined volumes of hi_eb_cc.areaHi() and hi_hi_eb_cc should
394  // equal vfrac(iv_hi).
395  { Real const vol = hi_eb_cc.volume() + hi_hi_eb_cc.volume();
396  Real const abs_err = amrex::Math::abs(vfrac(iv_hi) - vol);
397  Real compare_tol = 5.0e-6;
398 #ifndef AMREX_USE_GPU
399  if ( abs_err >= compare_tol ) {
400  hi_eb_cc.debug();
401  hi_hi_eb_cc.debug();
402  amrex::Print() << "\nFail: check-4 volume abs_err: " << abs_err
403  << "\n point: " << hi_point
404  << "\n normal: " << hi_normal
405  << "\n hi_eb_cc.volume() " << hi_eb_cc.volume()
406  << "\n hi_hi_eb_cc.volume() " << hi_hi_eb_cc.volume()
407  << "\n vfrac: " << vfrac(iv_hi)
408  << '\n';
409  } else if (verbose) {
410  Print() << "Pass: hi_eb_cc + hi_hi_eb_cc = vfrac = " << vfrac(iv_hi)
411  << " abs_err: " << abs_err << "\n";
412  }
413 #endif
414  AMREX_ALWAYS_ASSERT( abs_err < compare_tol );
415  }
416  } //
417 #endif
418 #endif // 0
419 
420  //-----------------------
421  // Fill out aux_ arrays
422  //-----------------------
423 
424  if (lo_eb_cc.isCovered() && hi_eb_cc.isCovered()) {
425 
426  // defaults to covered and disconnected.
427 
428  } else if (lo_eb_cc.isRegular() && hi_eb_cc.isRegular()) {
429 
430  aux_flag(i,j,k).setRegular();
431  aux_flag(i,j,k).setConnected(vdim);
432 
433  aux_vfrac(i,j,k) = 1.0;
434 
435  aux_afrac_x(i,j,k) = 1.0;
436  aux_afrac_y(i,j,k) = 1.0;
437  aux_afrac_z(i,j,k) = 1.0;
438 
439  aux_fcent_x(i,j,k,0) = 0.0; aux_fcent_x(i,j,k,1) = 0.0;
440  aux_fcent_y(i,j,k,0) = 0.0; aux_fcent_y(i,j,k,1) = 0.0;
441  aux_fcent_z(i,j,k,0) = 0.0; aux_fcent_z(i,j,k,1) = 0.0;
442 
443  if (i==bx.bigEnd(0)) {
444  aux_afrac_x(i+1,j,k) = 1.0;
445  aux_fcent_x(i+1,j,k,0) = 0.0; aux_fcent_x(i+1,j,k,1) = 0.0;
446  }
447  if (j==bx.bigEnd(1)) {
448  aux_afrac_y(i,j+1,k) = 1.0;
449  aux_fcent_y(i,j+1,k,0) = 0.0; aux_fcent_y(i,j+1,k,1) = 0.0;
450  }
451  if (k==bx.bigEnd(2)) {
452  aux_afrac_z(i,j,k+1) = 1.0;
453  aux_fcent_z(i,j,k+1,0) = 0.0; aux_fcent_z(i,j,k+1,1) = 0.0;
454  }
455 
456  } else if ( (lo_eb_cc.isRegular() && hi_eb_cc.isCovered())
457  || (lo_eb_cc.isCovered() && hi_eb_cc.isRegular()) ) {
458 
459  // This is a problematic situation.
460 #ifndef AMREX_USE_GPU
461  Print()<< "eb_aux_ / Check: Regular and Covered cut cells are facing each other." << std::endl;
462 #endif
463 
464  } else {
465 
466  aux_flag(i,j,k).setSingleValued();
467  aux_flag(i,j,k).setConnected(vdim);
468 
469  Real lo_vol {lo_eb_cc.volume()};
470  Real hi_vol {hi_eb_cc.volume()};
471 
472  aux_vfrac(i,j,k) = lo_vol + hi_vol;
473 
474  /* centVol() returns the coordinates based on m_rbx.
475  The coordinates in the idim direction are in [0.0,0.5] for the low cell and in [-0.5,0.0] for the hi cell.
476  Therefore, they need to be mapped to the eb_aux space, by shifting:
477  x' = x - 0.5 (low cell), x + 0.5 (hi cell) if idim = 0
478  y' = y - 0.5 (low cell), y + 0.5 (hi cell) if idim = 1
479  z' = z - 0.5 (low cell), z + 0.5 (hi cell) if idim = 2
480  */
481 
482  RealVect lo_vcent {lo_eb_cc.centVol()};
483  RealVect hi_vcent {hi_eb_cc.centVol()};
484 
485  lo_vcent[idim] = lo_vcent[idim] - 0.5;
486  hi_vcent[idim] = hi_vcent[idim] + 0.5;
487 
488  aux_vcent(i,j,k,0) = ( lo_vol * lo_vcent[0] + hi_vol * hi_vcent[0] ) / aux_vfrac(i,j,k);
489  aux_vcent(i,j,k,1) = ( lo_vol * lo_vcent[1] + hi_vol * hi_vcent[1] ) / aux_vfrac(i,j,k);
490  aux_vcent(i,j,k,2) = ( lo_vol * lo_vcent[2] + hi_vol * hi_vcent[2] ) / aux_vfrac(i,j,k);
491 
492  Real lo_areaLo_x {lo_eb_cc.areaLo(0)};
493  Real lo_areaLo_y {lo_eb_cc.areaLo(1)};
494  Real lo_areaLo_z {lo_eb_cc.areaLo(2)};
495 
496  Real hi_areaLo_x {hi_eb_cc.areaLo(0)};
497  Real hi_areaLo_y {hi_eb_cc.areaLo(1)};
498  Real hi_areaLo_z {hi_eb_cc.areaLo(2)};
499 
500  aux_afrac_x(i,j,k) = (idim == 0) ? lo_areaLo_x : lo_areaLo_x + hi_areaLo_x;
501  aux_afrac_y(i,j,k) = (idim == 1) ? lo_areaLo_y : lo_areaLo_y + hi_areaLo_y;
502  aux_afrac_z(i,j,k) = (idim == 2) ? lo_areaLo_z : lo_areaLo_z + hi_areaLo_z;
503 
504  /* fcentLo returns the coordinates based on m_rbx.
505  The coordinates in the idim direction are in [0.0,0.5] for the low cell and in [-0.5,0.0] for the hi cell.
506  Therefore, they need to be mapped to the eb_aux space, by shifting:
507  x' = x - 0.5 (low cell), x + 0.5 (hi cell) if idim = 0
508  y' = y - 0.5 (low cell), y + 0.5 (hi cell) if idim = 1
509  z' = z - 0.5 (low cell), z + 0.5 (hi cell) if idim = 2
510  */
511 
512  RealVect lo_centLo_x {lo_eb_cc.centLo(0)};
513  RealVect lo_centLo_y {lo_eb_cc.centLo(1)};
514  RealVect lo_centLo_z {lo_eb_cc.centLo(2)};
515 
516  RealVect hi_centLo_x {hi_eb_cc.centLo(0)};
517  RealVect hi_centLo_y {hi_eb_cc.centLo(1)};
518  RealVect hi_centLo_z {hi_eb_cc.centLo(2)};
519 
520  if (idim == 0) {
521  aux_fcent_x(i,j,k,0) = lo_centLo_x[1]; // y
522  aux_fcent_x(i,j,k,1) = lo_centLo_x[2]; // z
523  aux_fcent_y(i,j,k,0) = (aux_afrac_y(i,j,k) > 0.0) // x (mapped)
524  ? ( lo_areaLo_y * (lo_centLo_y[0] - 0.5)
525  + hi_areaLo_y * (hi_centLo_y[0] + 0.5) ) / aux_afrac_y(i,j,k)
526  : 0.0;
527  aux_fcent_y(i,j,k,1) = (aux_afrac_y(i,j,k) > 0.0) // z
528  ? ( lo_areaLo_y * lo_centLo_y[2]
529  + hi_areaLo_y * hi_centLo_y[2] ) / aux_afrac_y(i,j,k)
530  : 0.0;
531  aux_fcent_z(i,j,k,0) = (aux_afrac_z(i,j,k) > 0.0) // x (mapped)
532  ? ( lo_areaLo_z * (lo_centLo_z[0] - 0.5)
533  + hi_areaLo_z * (hi_centLo_z[0] + 0.5) ) / aux_afrac_z(i,j,k)
534  : 0.0;
535  aux_fcent_z(i,j,k,1) = (aux_afrac_z(i,j,k) > 0.0) // y
536  ? ( lo_areaLo_z * lo_centLo_z[1]
537  + hi_areaLo_z * hi_centLo_z[1] ) / aux_afrac_z(i,j,k)
538  : 0.0;
539  } else if (idim == 1) {
540  aux_fcent_x(i,j,k,0) = (aux_afrac_x(i,j,k) > 0.0) // y (mapped)
541  ? ( lo_areaLo_x * (lo_centLo_x[1] - 0.5)
542  + hi_areaLo_x * (hi_centLo_x[1] + 0.5) ) / aux_afrac_x(i,j,k)
543  : 0.0;
544  aux_fcent_x(i,j,k,1) = (aux_afrac_x(i,j,k) > 0.0) // z
545  ? ( lo_areaLo_x * lo_centLo_x[2]
546  + hi_areaLo_x * hi_centLo_x[2] ) / aux_afrac_x(i,j,k)
547  : 0.0;
548  aux_fcent_y(i,j,k,0) = lo_centLo_y[0]; // x
549  aux_fcent_y(i,j,k,1) = lo_centLo_y[2]; // z
550  aux_fcent_z(i,j,k,0) = (aux_afrac_z(i,j,k) > 0.0) // x
551  ? ( lo_areaLo_z * lo_centLo_z[0]
552  + hi_areaLo_z * hi_centLo_z[0] ) / aux_afrac_z(i,j,k)
553  : 0.0;
554  aux_fcent_z(i,j,k,1) = (aux_afrac_z(i,j,k) > 0.0) // y (mapped)
555  ? ( lo_areaLo_z * (lo_centLo_z[1] - 0.5)
556  + hi_areaLo_z * (hi_centLo_z[1] + 0.5) ) / aux_afrac_z(i,j,k)
557  : 0.0;
558  } else if (idim == 2) {
559  aux_fcent_x(i,j,k,0) = (aux_afrac_x(i,j,k) > 0.0) // y
560  ? ( lo_areaLo_x * lo_centLo_x[1]
561  + hi_areaLo_x * hi_centLo_x[1] ) / aux_afrac_x(i,j,k)
562  : 0.0;
563  aux_fcent_x(i,j,k,1) = (aux_afrac_x(i,j,k) > 0.0) // z (mapped)
564  ? ( lo_areaLo_x * (lo_centLo_x[2] - 0.5)
565  + hi_areaLo_x * (hi_centLo_x[2] + 0.5) ) / aux_afrac_x(i,j,k)
566  : 0.0;
567  aux_fcent_y(i,j,k,0) = (aux_afrac_y(i,j,k) > 0.0) // x
568  ? ( lo_areaLo_y * lo_centLo_y[0]
569  + hi_areaLo_y * hi_centLo_y[0] ) / aux_afrac_y(i,j,k)
570  : 0.0;
571  aux_fcent_y(i,j,k,1) = (aux_afrac_y(i,j,k) > 0.0) // z (mapped)
572  ? ( lo_areaLo_y * (lo_centLo_y[2] - 0.5)
573  + hi_areaLo_y * (hi_centLo_y[2] + 0.5) ) / aux_afrac_y(i,j,k)
574  : 0.0;
575  aux_fcent_z(i,j,k,0) = lo_centLo_z[0]; // x
576  aux_fcent_z(i,j,k,1) = lo_centLo_z[1]; // y
577  }
578 
579  // Need to fill the nodes the big ends?
580 
581  Real lo_areaBoun {lo_eb_cc.areaBoun()};
582  Real hi_areaBoun {hi_eb_cc.areaBoun()};
583 
584  aux_barea(i,j,k) = lo_areaBoun + hi_areaBoun;
585 
586  RealVect lo_centBoun {lo_eb_cc.centBoun()};
587  RealVect hi_centBoun {hi_eb_cc.centBoun()};
588 
589  if (idim == 0) {
590  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); // x (mapped)
591  aux_bcent(i,j,k,1) = ( lo_areaBoun * lo_centBoun[1] + hi_areaBoun * hi_centBoun[1] ) / aux_barea(i,j,k); // y
592  aux_bcent(i,j,k,2) = ( lo_areaBoun * lo_centBoun[2] + hi_areaBoun * hi_centBoun[2] ) / aux_barea(i,j,k); // z
593  } else if (idim == 1) {
594  aux_bcent(i,j,k,0) = ( lo_areaBoun * lo_centBoun[0] + hi_areaBoun * hi_centBoun[0] ) / aux_barea(i,j,k); // x
595  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); // y (mapped)
596  aux_bcent(i,j,k,2) = ( lo_areaBoun * lo_centBoun[2] + hi_areaBoun * hi_centBoun[2] ) / aux_barea(i,j,k); // z
597  } else if (idim == 2) {
598  aux_bcent(i,j,k,0) = ( lo_areaBoun * lo_centBoun[0] + hi_areaBoun * hi_centBoun[0] ) / aux_barea(i,j,k); // x
599  aux_bcent(i,j,k,1) = ( lo_areaBoun * lo_centBoun[1] + hi_areaBoun * hi_centBoun[1] ) / aux_barea(i,j,k); // y
600  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); // z (mapped)
601  }
602 
603  RealVect eb_normal = ( lo_areaBoun * lo_normal + hi_areaBoun * hi_normal )/ aux_barea(i,j,k);
604 
605  aux_bnorm(i,j,k,0) = eb_normal[0];
606  aux_bnorm(i,j,k,1) = eb_normal[1];
607  aux_bnorm(i,j,k,2) = eb_normal[2];
608  }
609 
610  } // flag(iv_lo) and flag(iv_hi)
611 
612  });
613 
614  }
615 
616  }
617 }
eb_aux_::m_bndrycent
amrex::MultiCutFab * m_bndrycent
Definition: ERF_EBAux.H:48
eb_aux_::m_bndryarea
amrex::MultiCutFab * m_bndryarea
Definition: ERF_EBAux.H:47
eb_aux_::m_cellflags
amrex::FabArray< amrex::EBCellFlagFab > * m_cellflags
Definition: ERF_EBAux.H:42
eb_aux_::m_volfrac
amrex::MultiFab * m_volfrac
Definition: ERF_EBAux.H:44
eb_aux_::m_bndrynorm
amrex::MultiCutFab * m_bndrynorm
Definition: ERF_EBAux.H:49
eb_aux_::m_areafrac
amrex::Array< amrex::MultiCutFab *, AMREX_SPACEDIM > m_areafrac
Definition: ERF_EBAux.H:51
eb_aux_::m_volcent
amrex::MultiCutFab * m_volcent
Definition: ERF_EBAux.H:46
eb_aux_::m_facecent
amrex::Array< amrex::MultiCutFab *, AMREX_SPACEDIM > m_facecent
Definition: ERF_EBAux.H:52
eb_cut_cell_
Definition: ERF_EBCutCell.H:14

Referenced by eb_::make_factory().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ getAreaFrac()

Array< const MultiCutFab *, AMREX_SPACEDIM > eb_aux_::getAreaFrac ( ) const
656 {
657  AMREX_ASSERT(m_areafrac[0] != nullptr);
658  return {AMREX_D_DECL(m_areafrac[0], m_areafrac[1], m_areafrac[2])};
659 }

◆ getBndryArea()

const MultiCutFab & eb_aux_::getBndryArea ( ) const
635 {
636  AMREX_ASSERT(m_bndryarea != nullptr);
637  return *m_bndryarea;
638 }

◆ getBndryCent()

const MultiCutFab & eb_aux_::getBndryCent ( ) const
642 {
643  AMREX_ASSERT(m_bndrycent != nullptr);
644  return *m_bndrycent;
645 }

◆ getBndryNorm()

const MultiCutFab & eb_aux_::getBndryNorm ( ) const
649 {
650  AMREX_ASSERT(m_bndrynorm != nullptr);
651  return *m_bndrynorm;
652 }

◆ getFaceCent()

Array< const MultiCutFab *, AMREX_SPACEDIM > eb_aux_::getFaceCent ( ) const
663 {
664  AMREX_ASSERT(m_facecent[0] != nullptr);
665  return {AMREX_D_DECL(m_facecent[0], m_facecent[1], m_facecent[2])};
666 }

◆ getVolCent()

const MultiCutFab & eb_aux_::getVolCent ( ) const
628 {
629  AMREX_ASSERT(m_volcent != nullptr);
630  return *m_volcent;
631 }

◆ getVolFrac()

const MultiFab & eb_aux_::getVolFrac ( ) const
621 {
622  AMREX_ASSERT(m_volfrac != nullptr);
623  return *m_volfrac;
624 }

◆ set_verbose()

void eb_aux_::set_verbose ( )
inline
25 { m_verbose = 1; }

Member Data Documentation

◆ m_areafrac

amrex::Array<amrex::MultiCutFab*,AMREX_SPACEDIM> eb_aux_::m_areafrac {{AMREX_D_DECL(nullptr, nullptr, nullptr)}}
private

Referenced by define(), and getAreaFrac().

◆ m_bndryarea

amrex::MultiCutFab* eb_aux_::m_bndryarea = nullptr
private

Referenced by define(), and getBndryArea().

◆ m_bndrycent

amrex::MultiCutFab* eb_aux_::m_bndrycent = nullptr
private

Referenced by define(), and getBndryCent().

◆ m_bndrynorm

amrex::MultiCutFab* eb_aux_::m_bndrynorm = nullptr
private

Referenced by define(), and getBndryNorm().

◆ m_cellflags

amrex::FabArray<amrex::EBCellFlagFab>* eb_aux_::m_cellflags = nullptr
private

Referenced by define().

◆ m_facecent

amrex::Array<amrex::MultiCutFab*,AMREX_SPACEDIM> eb_aux_::m_facecent {{AMREX_D_DECL(nullptr, nullptr, nullptr)}}
private

Referenced by define(), and getFaceCent().

◆ m_verbose

int eb_aux_::m_verbose
private

Referenced by define(), and set_verbose().

◆ m_volcent

amrex::MultiCutFab* eb_aux_::m_volcent = nullptr
private

Referenced by define(), and getVolCent().

◆ m_volfrac

amrex::MultiFab* eb_aux_::m_volfrac = nullptr
private

Referenced by define(), and getVolFrac().

The documentation for this class was generated from the following files: