ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
ERF_EBCutCell.H
Go to the documentation of this file.
1 #ifndef ERF_EB_CUT_CELL_H_
2 #define ERF_EB_CUT_CELL_H_
3 
4 #include <AMReX_REAL.H>
5 #include <AMReX_Array.H>
6 #include <AMReX_RealBox.H>
7 #include <AMReX_RealVect.H>
8 #include <AMReX_Algorithm.H>
9 #include <AMReX_EBCellFlag.H>
10 
11 #include "ERF_EBPolygon.H"
12 #include "ERF_EBUtils.H"
13 
14 class eb_cut_cell_ {
15 
16  public:
17 
18  AMREX_GPU_HOST_DEVICE
19  eb_cut_cell_ ( amrex::EBCellFlag const& a_flag,
20  amrex::RealBox const& a_rbox,
21  amrex::RealVect const& a_point,
22  amrex::RealVect const& a_normal );
23 
24  [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
25  bool isCovered () const noexcept { return m_flag.isCovered(); }
26 
27  [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
28  bool isRegular () const noexcept { return m_flag.isRegular(); }
29 
30  [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
31  bool isSingleValued () const noexcept { return m_flag.isSingleValued(); }
32 
33  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
34  void set_covered () {
35  m_flag.setCovered();
36  }
37 
38  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
39  void set_regular () {
40 
41  m_flag.setRegular();
42 
43  m_F1.set_area( m_rbox.length(0)*m_rbox.length(1) );
44  m_F3.set_area( m_rbox.length(0)*m_rbox.length(1) );
45 
46  m_F2.set_area( m_rbox.length(1)*m_rbox.length(2) );
47  m_F4.set_area( m_rbox.length(1)*m_rbox.length(2) );
48 
49  m_F5.set_area( m_rbox.length(0)*m_rbox.length(2) );
50  m_F6.set_area( m_rbox.length(0)*m_rbox.length(2) );
51 
52  m_F7.set_area(0.);
53 
54  }
55 
56  [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
57  amrex::Real volume () const {
58 
59  if (m_flag.isCovered() ) { return 0.; }
60  if (m_flag.isRegular() ) { return m_rbox.volume(); }
61 
62  amrex::Real volume(0.);
63 
64  amrex::Real const* lo = m_rbox.lo();
65  amrex::RealVect v0(lo[0], lo[1], lo[2]);
66 
67  if( m_F2.ok() ) { volume += m_F2.area() * m_F2.distance(v0); }
68  if (m_F3.ok() ) { volume += m_F3.area() * m_F3.distance(v0); }
69  if (m_F6.ok() ) { volume += m_F6.area() * m_F6.distance(v0); }
70  if (m_F7.ok() ) { volume += m_F7.area() * m_F7.distance(v0); }
71 
72  volume /= 3.;
73 
74  return m_invert*volume + (1.-m_invert)*(m_rbox.volume()-volume);
75  }
76 
77  [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
78  amrex::Real areaLo ( int const a_idim ) const noexcept {
79  AMREX_ASSERT( a_idim >=0 && a_idim < AMREX_SPACEDIM );
80  if (m_flag.isCovered() ) { return 0.; }
81  if (m_flag.isRegular() ) { return m_lo_faces[a_idim]->area(); }
82  amrex::Real const area(m_lo_faces[a_idim]->area());
83  return m_invert*area + (1.-m_invert)*(m_rbox_area[a_idim] - area);
84  }
85 
86  [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
87  amrex::Real areaHi ( int const a_idim ) const noexcept {
88  AMREX_ASSERT( a_idim >= 0 && a_idim < AMREX_SPACEDIM );
89  if (m_flag.isCovered() ) { return 0.; }
90  if (m_flag.isRegular() ) { return m_hi_faces[a_idim]->area(); }
91  amrex::Real const area(m_hi_faces[a_idim]->area());
92  return m_invert*area + (1.-m_invert)*(m_rbox_area[a_idim] - area);
93  }
94 
95  [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
96  amrex::Real areaBoun () const noexcept {
97  return m_F7.area();
98  }
99 
100  [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
101  amrex::RealVect centLo ( int const a_idim ) const noexcept {
102  AMREX_ASSERT( a_idim >=0 && a_idim < AMREX_SPACEDIM );
103  amrex::RealVect cent = m_cellface_cent[ m_lo_faces_id[a_idim] ];
104  if (m_flag.isCovered() || m_flag.isRegular() ) {
105  // Default cent
106  } else {
107  amrex::Real const area_R = m_lo_faces[a_idim]->area();
108  if (amrex::almostEqual(area_R,0.0) || amrex::almostEqual(area_R, m_rbox_area[a_idim]) ){
109  // Default cent
110  } else {
111  amrex::RealVect cent_O = cent;
112  amrex::RealVect cent_R = m_lo_faces[a_idim]->get_centroid();
113  amrex::Real const area_C = m_rbox_area[a_idim] - area_R;
114  cent = m_invert * cent_R + (1.-m_invert) * ((1.+area_R/area_C)*cent_O - area_R/area_C*cent_R);
115  }
116  }
117  return cent;
118  }
119 
120  [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
121  amrex::RealVect centHi ( int const a_idim ) const noexcept {
122  AMREX_ASSERT( a_idim >=0 && a_idim < AMREX_SPACEDIM );
123  amrex::RealVect cent = m_cellface_cent[ m_hi_faces_id[a_idim] ];
124  if (m_flag.isCovered() || m_flag.isRegular() ) {
125  // Default cent
126  } else {
127  amrex::Real const area_R = m_hi_faces[a_idim]->area();
128  if (amrex::almostEqual(area_R,0.0) || amrex::almostEqual(area_R, m_rbox_area[a_idim]) ){
129  // Default cent
130  } else {
131  amrex::RealVect cent_O = cent;
132  amrex::RealVect cent_R = m_hi_faces[a_idim]->get_centroid();
133  amrex::Real const area_C = m_rbox_area[a_idim] - area_R;
134  cent = m_invert * cent_R + (1.-m_invert) * ((1.+area_R/area_C)*cent_O - area_R/area_C*cent_R);
135  }
136  }
137  return cent;
138  }
139 
140  [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
141  amrex::RealVect centBoun () const noexcept {
142  amrex::RealVect cent{0.,0.,0.};
143  if (m_flag.isSingleValued()) {
144  cent = m_F7.get_centroid();
145  }
146  return cent;
147  }
148 
149  [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
150  amrex::RealVect normBoun () const noexcept {
151  amrex::RealVect normal{0.,0.,0.};
152  if (m_flag.isSingleValued()) {
153  normal = m_eb_normal;
154  }
155  return normal;
156  }
157 
158  [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
159  amrex::RealVect centVol () const noexcept {
160  amrex::RealVect vcent = m_cell_cent;
161  if (m_flag.isSingleValued()) {
162  amrex::Real xm = m_rbox.lo(0);
163  amrex::Real xp = m_rbox.hi(0);
164  amrex::Real ym = m_rbox.lo(1);
165  amrex::Real yp = m_rbox.hi(1);
166  amrex::Real zm = m_rbox.lo(2);
167  amrex::Real zp = m_rbox.hi(2);
168 
169  amrex::Real axm = areaLo(0);
170  amrex::Real axp = areaHi(0);
171  amrex::Real aym = areaLo(1);
172  amrex::Real ayp = areaHi(1);
173  amrex::Real azm = areaLo(2);
174  amrex::Real azp = areaHi(2);
175 
176  amrex::Real barea = m_F7.area();
177  amrex::RealVect bcent = m_F7.get_centroid();
178 
179  amrex::Real vol = volume();
180 
181  vcent[0] = 0.5 * ( - xm * xm * axm + xp * xp * axp + bcent[0] * bcent[0] * m_eb_normal[0] * barea ) / vol;
182  vcent[1] = 0.5 * ( - ym * ym * aym + yp * yp * ayp + bcent[1] * bcent[1] * m_eb_normal[1] * barea ) / vol;
183  vcent[2] = 0.5 * ( - zm * zm * azm + zp * zp * azp + bcent[2] * bcent[2] * m_eb_normal[2] * barea ) / vol;
184  }
185  return vcent;
186  }
187 
188  void debug ( int const a_face = -1 );
189 
190  private:
191 
192  amrex::RealBox const m_rbox;
193  amrex::RealVect const m_eb_point;
194  amrex::RealVect const m_eb_normal;
195 
196  amrex::Real m_invert;
197 
198  amrex::RealVect m_rbox_area;
199 
200  amrex::EBCellFlag m_flag;
201 
202  // Cell faces
203 
204  polygon_ m_F1;
205  polygon_ m_F2;
206  polygon_ m_F3;
207  polygon_ m_F4;
208  polygon_ m_F5;
209  polygon_ m_F6;
210 
211  static int constexpr m_max_faces = 6;
212  amrex::Array<amrex::RealVect,m_max_faces> m_cellface_cent; // Centroids of cell faces
213  amrex::RealVect m_cell_cent; // Centroid of cell
214  amrex::RealVect m_dx;
215 
216  // cut face
217  polygon_ m_F7;
218 
219  amrex::Array<polygon_ const* const,3> m_lo_faces {&m_F4, &m_F5, &m_F1};
220  amrex::Array<polygon_ const* const,3> m_hi_faces {&m_F2, &m_F6, &m_F3};
221 
222  amrex::Array<int const,3> m_lo_faces_id {3,4,0};
223  amrex::Array<int const,3> m_hi_faces_id {1,5,2};
224 
225  AMREX_GPU_HOST_DEVICE
226  void calc_edge_intersections ( int const a_dry_run = 0 );
227 
228  AMREX_GPU_HOST_DEVICE
229  void set_covered_regular_cell_vertices ();
230 
231 };
232 
233 AMREX_GPU_HOST_DEVICE
234 AMREX_FORCE_INLINE
235 eb_cut_cell_::
236 eb_cut_cell_ ( amrex::EBCellFlag const& a_flag,
237  amrex::RealBox const& a_rbox,
238  amrex::RealVect const& a_point,
239  amrex::RealVect const& a_normal )
240  : m_rbox(a_rbox)
241  , m_eb_point(a_point)
242  , m_eb_normal(a_normal)
243  , m_invert(0.0)
244  , m_F1(a_point, a_normal)
245  , m_F2(a_point, a_normal)
246  , m_F3(a_point, a_normal)
247  , m_F4(a_point, a_normal)
248  , m_F5(a_point, a_normal)
249  , m_F6(a_point, a_normal)
250  , m_cellface_cent({amrex::RealVect(0.), amrex::RealVect(0.), amrex::RealVect(0.),
251  amrex::RealVect(0.), amrex::RealVect(0.), amrex::RealVect(0.)})
252 {
253  using namespace amrex;
254 
255  m_rbox_area[0] = m_rbox.length(1)*m_rbox.length(2);
256  m_rbox_area[1] = m_rbox.length(0)*m_rbox.length(2);
257  m_rbox_area[2] = m_rbox.length(0)*m_rbox.length(1);
258 
259  RealVect v0(m_rbox.lo(0), m_rbox.lo(1), m_rbox.lo(2));
260  RealVect v1(m_rbox.hi(0), m_rbox.lo(1), m_rbox.lo(2));
261  RealVect v2(m_rbox.lo(0), m_rbox.hi(1), m_rbox.lo(2));
262  RealVect v3(m_rbox.lo(0), m_rbox.lo(1), m_rbox.hi(2));
263  RealVect v4(m_rbox.hi(0), m_rbox.lo(1), m_rbox.hi(2));
264  RealVect v5(m_rbox.hi(0), m_rbox.hi(1), m_rbox.lo(2));
265  RealVect v6(m_rbox.lo(0), m_rbox.hi(1), m_rbox.hi(2));
266  RealVect v7(m_rbox.hi(0), m_rbox.hi(1), m_rbox.hi(2));
267 
268  // Centoids of cell faces
269 
270  m_cellface_cent[0] = 0.25 * ( v0 + v2 + v5 + v1 ); // F1
271  m_cellface_cent[1] = 0.25 * ( v1 + v5 + v7 + v4 ); // F2
272  m_cellface_cent[2] = 0.25 * ( v3 + v4 + v7 + v6 ); // F3
273  m_cellface_cent[3] = 0.25 * ( v0 + v3 + v6 + v2 ); // F4
274  m_cellface_cent[4] = 0.25 * ( v0 + v1 + v4 + v3 ); // F5
275  m_cellface_cent[5] = 0.25 * ( v2 + v5 + v7 + v6 ); // F6
276 
277  // Cell centroid
278 
279  m_cell_cent[0] = 0.5 * ( m_rbox.lo(0) + m_rbox.hi(0) );
280  m_cell_cent[1] = 0.5 * ( m_rbox.lo(1) + m_rbox.hi(1) );
281  m_cell_cent[2] = 0.5 * ( m_rbox.lo(2) + m_rbox.hi(2) );
282 
283  // Cell size
284 
285  m_dx[0] = m_rbox.hi(0) - m_rbox.lo(0);
286  m_dx[1] = m_rbox.hi(1) - m_rbox.lo(1);
287  m_dx[2] = m_rbox.hi(2) - m_rbox.lo(2);
288 
289  if (a_flag.isCovered() ) {
290 
291  set_covered();
292 
293  } else if (a_flag.isRegular() ) {
294 
295  set_regular();
296 
297  } else { // Check that the box and plane intersect.
298 
299  RealVect c = 0.5*(v0 + v7);
300  RealVect e = v7 - c;
301 
302  Real r = e[0]*amrex::Math::abs(a_normal[0]) +
303  e[1]*amrex::Math::abs(a_normal[1]) +
304  e[2]*amrex::Math::abs(a_normal[2]);
305 
306  Real s = amrex::Math::abs(c.dotProduct(a_normal)
307  - a_point.dotProduct(a_normal));
308 
309  if (s > r) {
310  if (a_normal.dotProduct(v0 - a_point) > 0.)
311  { set_covered(); } else { set_regular(); }
312  } else { m_flag.setSingleValued(); }
313  }
314 
315  if ( m_flag.isSingleValued() ) {
316 
317  m_invert = ((m_eb_normal.dotProduct(v0 - m_eb_point)) > 0.) ? 0.0 : 1.0;
318 
319  calc_edge_intersections();
320 
321  } // end singleValued
322 
323  m_F1.define();
324  m_F2.define();
325  m_F3.define();
326  m_F4.define();
327  m_F5.define();
328  m_F6.define();
329  m_F7.define();
330 
331  // For covered and regular cut cells, add vertices to utilize e.g., get_centroid.
332  if ( !m_flag.isSingleValued() ) {
333  set_covered_regular_cell_vertices();
334  }
335 
336 }
337 
338 AMREX_GPU_HOST_DEVICE
339 AMREX_FORCE_INLINE
340 void
341 eb_cut_cell_::
342 calc_edge_intersections ( int const a_dry_run )
343 {
344  using namespace amrex;
345 
346  RealVect v0(m_rbox.lo(0), m_rbox.lo(1), m_rbox.lo(2));
347  RealVect v1(m_rbox.hi(0), m_rbox.lo(1), m_rbox.lo(2));
348  RealVect v2(m_rbox.lo(0), m_rbox.hi(1), m_rbox.lo(2));
349  RealVect v3(m_rbox.lo(0), m_rbox.lo(1), m_rbox.hi(2));
350  RealVect v4(m_rbox.hi(0), m_rbox.lo(1), m_rbox.hi(2));
351  RealVect v5(m_rbox.hi(0), m_rbox.hi(1), m_rbox.lo(2));
352  RealVect v6(m_rbox.lo(0), m_rbox.hi(1), m_rbox.hi(2));
353  RealVect v7(m_rbox.hi(0), m_rbox.hi(1), m_rbox.hi(2));
354 
355  if (!a_dry_run) {
356  m_F1.add_vertex(v0);
357  m_F4.add_vertex(v0);
358  m_F5.add_vertex(v0);
359  }
360 
361  int add_v7(1);
362 
363  RealVect vIP;
364  Real distIP;
365 
366  // Path 1
367  { int cuts(0);
368 
369  if (utils::intersect_plane_edge(m_eb_point, m_eb_normal, v0, v1, vIP, distIP)) {
370 
371  ++cuts;
372 #ifndef AMREX_USE_GPU
373  if (a_dry_run) { Print() << "P1: v0--v1: add vIP to F1, F5 and F7 :: " << vIP[0] << "\n"; }
374  else
375 #endif
376  {
377  m_F1.add_vertex(vIP);
378  m_F5.add_vertex(vIP);
379  m_F7.add_vertex(vIP);
380  }
381 
382  if ( almostEqual(distIP, 0.0) ) {
383 #ifndef AMREX_USE_GPU
384  if (a_dry_run) { Print() << "P1: v0--v1: intersection ~ 0 :: add vIP to F4\n"; }
385  else
386 #endif
387  { m_F4.add_vertex(vIP); }
388 
389  } else if ( almostEqual(distIP, 1.0) ) {
390 #ifndef AMREX_USE_GPU
391  if (a_dry_run) { Print() << "P1: v0--v1: intersection ~ 1 :: add vIP to F2\n"; }
392  else
393 #endif
394  { m_F2.add_vertex(vIP); }
395  }
396 
397  }
398 
399  if (cuts%2 == 0) {
400 
401 #ifndef AMREX_USE_GPU
402  if (a_dry_run) { Print() << "P1: Add v1 to F1, F2 and F5\n"; }
403  else
404 #endif
405  {
406  m_F1.add_vertex(v1);
407  m_F2.add_vertex(v1);
408  m_F5.add_vertex(v1);
409  }
410  }
411 
412  if (utils::intersect_plane_edge(m_eb_point, m_eb_normal, v1, v4, vIP, distIP)) {
413 
414  ++cuts;
415 
416 #ifndef AMREX_USE_GPU
417  if (a_dry_run) { Print() << "P1: v1--v4: add vIP to F2, F5 and F7 :: " << vIP[2] << "\n"; }
418  else
419 #endif
420  {
421  m_F2.add_vertex(vIP);
422  m_F5.add_vertex(vIP);
423  m_F7.add_vertex(vIP);
424  }
425 
426  if ( almostEqual(distIP, 0.0) ) {
427 #ifndef AMREX_USE_GPU
428  if (a_dry_run) { Print() << "P1: v1--v4: intersection ~ 0 :: add vIP to F1\n"; }
429  else
430 #endif
431  {m_F1.add_vertex(vIP); }
432 
433  } else if ( almostEqual(distIP, 1.0) ) {
434 #ifndef AMREX_USE_GPU
435  if (a_dry_run) { Print() << "P1: v1--v4: intersection ~ 1 :: add vIP to F3\n"; }
436  else
437 #endif
438  { m_F3.add_vertex(vIP); }
439  }
440 
441  }
442 
443  if (cuts%2 == 0) {
444 
445 #ifndef AMREX_USE_GPU
446  if (a_dry_run) { Print() << "P1: Add v4 to F2, F3 and F5\n"; }
447  else
448 #endif
449  {
450  m_F2.add_vertex(v4);
451  m_F3.add_vertex(v4);
452  m_F5.add_vertex(v4);
453  }
454  }
455 
456  if (utils::intersect_plane_edge(m_eb_point, m_eb_normal, v4, v7, vIP, distIP)) {
457 
458  ++cuts;
459 
460 #ifndef AMREX_USE_GPU
461  if (a_dry_run) { Print() << "P1: v4--v7: add vIP to F2, F3 and F7 :: " << vIP[1] << "\n"; }
462  else
463 #endif
464  {
465  m_F2.add_vertex(vIP);
466  m_F3.add_vertex(vIP);
467  m_F7.add_vertex(vIP);
468  }
469 
470  if ( almostEqual(distIP, 0.0)) {
471 #ifndef AMREX_USE_GPU
472  if (a_dry_run) { Print() << "P1: v4--v7: intersection ~ 0 :: add vIP to F5\n"; }
473  else
474 #endif
475  { m_F5.add_vertex(vIP); }
476 
477  } else if ( almostEqual(distIP, 1.0)) {
478 #ifndef AMREX_USE_GPU
479  if (a_dry_run) { Print() << "P1: v4--v7: intersection ~ 1 :: add vIP to F6\n"; }
480  else
481 #endif
482  { m_F6.add_vertex(vIP); }
483  }
484  }
485 
486  if (cuts == 2 && add_v7) {
487 
488 #ifndef AMREX_USE_GPU
489  if (a_dry_run) { Print() << "P1: Add v7 to F2, F3 and F6\n"; }
490  else
491 #endif
492  {
493  m_F2.add_vertex(v7);
494  m_F3.add_vertex(v7);
495  m_F6.add_vertex(v7);
496  }
497 
498  add_v7 = 0;
499  }
500  } // end Path 1
501 
502  // Path 4
503  if (utils::intersect_plane_edge(m_eb_point, m_eb_normal, v1, v5, vIP, distIP)) {
504 
505 #ifndef AMREX_USE_GPU
506  if (a_dry_run) { Print() << "P4: v1--v5: add vIP to F1, F2 and F7 :: " << vIP[1] << "\n"; }
507  else
508 #endif
509  {
510  m_F1.add_vertex(vIP);
511  m_F2.add_vertex(vIP);
512  m_F7.add_vertex(vIP);
513  }
514 
515  if ( almostEqual(distIP, 0.0)) {
516 #ifndef AMREX_USE_GPU
517  if (a_dry_run) { Print() << "P1: v1--v5: intersection ~ 0 :: add vIP to F5\n"; }
518  else
519 #endif
520  { m_F5.add_vertex(vIP); }
521 
522  } else if ( almostEqual(distIP, 1.0)) {
523 #ifndef AMREX_USE_GPU
524  if (a_dry_run) { Print() << "P1: v1--v5: intersection ~ 1 :: add vIP to F6\n"; }
525  else
526 #endif
527  { m_F6.add_vertex(vIP); }
528  }
529  }
530 
531 
532  // Path 2
533  { int cuts(0);
534 
535  if (utils::intersect_plane_edge(m_eb_point, m_eb_normal, v0, v2, vIP, distIP)) {
536 
537  ++cuts;
538 
539 #ifndef AMREX_USE_GPU
540  if (a_dry_run) { Print() << "P2: v0--v2: add vIP to F1, F4 and F7 :: " << vIP[1] << "\n"; }
541  else
542 #endif
543  {
544  m_F1.add_vertex(vIP);
545  m_F4.add_vertex(vIP);
546  m_F7.add_vertex(vIP);
547  }
548 
549  if ( almostEqual(distIP, 0.0)) {
550 #ifndef AMREX_USE_GPU
551  if (a_dry_run) { Print() << "P2: v0--v2: intersection ~ 0 :: add vIP to F5\n"; }
552  else
553 #endif
554  { m_F5.add_vertex(vIP); }
555 
556  } else if ( almostEqual(distIP, 1.0)) {
557 #ifndef AMREX_USE_GPU
558  if (a_dry_run) { Print() << "P2: v0--v2: intersection ~ 1 :: add vIP to F6\n"; }
559  else
560 #endif
561  { m_F6.add_vertex(vIP); }
562  }
563 
564  }
565 
566  if (cuts%2 == 0) {
567 
568 #ifndef AMREX_USE_GPU
569  if (a_dry_run) { Print() << "P2: Add v2 to F1, F4, F6\n"; }
570  else
571 #endif
572  {
573  m_F1.add_vertex(v2);
574  m_F4.add_vertex(v2);
575  m_F6.add_vertex(v2);
576  }
577  }
578 
579  if (utils::intersect_plane_edge(m_eb_point, m_eb_normal, v2, v5, vIP, distIP)) {
580 
581  ++cuts;
582 
583 #ifndef AMREX_USE_GPU
584  if (a_dry_run) { Print() << "P2: v2--v5: add vIP to F1, F6 and F7 :: " << vIP[0] << "\n"; }
585  else
586 #endif
587  {
588  m_F1.add_vertex(vIP);
589  m_F6.add_vertex(vIP);
590  m_F7.add_vertex(vIP);
591  }
592 
593  if ( almostEqual(distIP, 0.0)) {
594 #ifndef AMREX_USE_GPU
595  if (a_dry_run) { Print() << "P2: v2--v5: intersection ~ 0 :: add vIP to F4\n"; }
596  else
597 #endif
598  { m_F4.add_vertex(vIP); }
599 
600  } else if ( almostEqual(distIP, 1.0)) {
601 #ifndef AMREX_USE_GPU
602  if (a_dry_run) { Print() << "P2: v2--v5: intersection ~ 1 :: add vIP to F2\n"; }
603  else
604 #endif
605  { m_F2.add_vertex(vIP); }
606  }
607  }
608 
609  if (cuts%2 == 0) {
610 
611 #ifndef AMREX_USE_GPU
612  if (a_dry_run) { Print() << "P2: Add v5 to F1, F2 and F6\n"; }
613  else
614 #endif
615  {
616  m_F1.add_vertex(v5);
617  m_F2.add_vertex(v5);
618  m_F6.add_vertex(v5);
619  }
620  }
621 
622  if (utils::intersect_plane_edge(m_eb_point, m_eb_normal, v5, v7, vIP, distIP)) {
623 
624  ++cuts;
625 
626 #ifndef AMREX_USE_GPU
627  if (a_dry_run) { Print() << "P2: v5--v7: add vIP to F2, F6 and F7 :: " << vIP[2] << "\n"; }
628  else
629 #endif
630  {
631  m_F2.add_vertex(vIP);
632  m_F6.add_vertex(vIP);
633  m_F7.add_vertex(vIP);
634  }
635 
636  if ( almostEqual(distIP, 0.0)) {
637 #ifndef AMREX_USE_GPU
638  if (a_dry_run) { Print() << "P2: v5--v7: intersection ~ 0 :: add vIP to F1\n"; }
639  else
640 #endif
641  { m_F1.add_vertex(vIP); }
642 
643  } else if ( almostEqual(distIP, 1.0)) {
644 #ifndef AMREX_USE_GPU
645  if (a_dry_run) { Print() << "P2: v5--v7: intersection ~ 1 :: add vIP to F3\n"; }
646  else
647 #endif
648  { m_F3.add_vertex(vIP); }
649  }
650  }
651 
652  if (cuts == 2 && add_v7) {
653 
654 #ifndef AMREX_USE_GPU
655  if (a_dry_run) { Print() << "P2: Add v7 to F2, F3 and F6\n"; }
656  else
657 #endif
658  {
659  m_F2.add_vertex(v7);
660  m_F3.add_vertex(v7);
661  m_F6.add_vertex(v7);
662  }
663 
664  add_v7 = 0;
665  }
666 
667  } // end Path 2
668 
669  // Path 5
670  if (utils::intersect_plane_edge(m_eb_point, m_eb_normal, v2, v6, vIP, distIP)) {
671 
672 #ifndef AMREX_USE_GPU
673  if (a_dry_run) { Print() << "P5: v2--v6: add vIP to F4, F6 and F7 :: " << vIP[2] << "\n"; }
674  else
675 #endif
676  {
677  m_F4.add_vertex(vIP);
678  m_F6.add_vertex(vIP);
679  m_F7.add_vertex(vIP);
680  }
681 
682  if ( almostEqual(distIP, 0.0)) {
683 #ifndef AMREX_USE_GPU
684  if (a_dry_run) { Print() << "P5: v2--v6: intersection ~ 0 :: add vIP to F1\n"; }
685  else
686 #endif
687  { m_F1.add_vertex(vIP); }
688 
689  } else if ( almostEqual(distIP, 1.0)) {
690 #ifndef AMREX_USE_GPU
691  if (a_dry_run) { Print() << "P5: v2--v6: intersection ~ 1 :: add vIP to F3\n"; }
692  else
693 #endif
694  { m_F3.add_vertex(vIP); }
695  }
696  }
697 
698 
699  // Path 3
700  { int cuts(0);
701 
702  if (utils::intersect_plane_edge(m_eb_point, m_eb_normal, v0, v3, vIP, distIP)) {
703 
704  ++cuts;
705 
706 #ifndef AMREX_USE_GPU
707  if (a_dry_run) { Print() << "P3: v0--v3: add vIP to F4, F5 and F7 :: " << vIP[2] << "\n"; }
708  else
709 #endif
710  {
711  m_F4.add_vertex(vIP);
712  m_F5.add_vertex(vIP);
713  m_F7.add_vertex(vIP);
714  }
715 
716  if ( almostEqual(distIP, 0.0)) {
717 #ifndef AMREX_USE_GPU
718  if (a_dry_run) { Print() << "P3: v0--v3: intersection ~ 0 :: add vIP to F1\n"; }
719  else
720 #endif
721  { m_F1.add_vertex(vIP); }
722 
723  } else if ( almostEqual(distIP, 1.0)) {
724 #ifndef AMREX_USE_GPU
725  if (a_dry_run) { Print() << "P3: v0--v3: intersection ~ 1 :: add vIP to F3\n"; }
726  else
727 #endif
728  { m_F3.add_vertex(vIP); }
729  }
730 
731  }
732 
733  if (cuts%2 == 0) {
734 
735 #ifndef AMREX_USE_GPU
736  if (a_dry_run) { Print() << "P3: Add v3 to F3, F4 and F5\n"; }
737  else
738 #endif
739  {
740  m_F3.add_vertex(v3);
741  m_F4.add_vertex(v3);
742  m_F5.add_vertex(v3);
743  }
744  }
745 
746  if (utils::intersect_plane_edge(m_eb_point, m_eb_normal, v3, v6, vIP, distIP)) {
747 
748  ++cuts;
749 
750 #ifndef AMREX_USE_GPU
751  if (a_dry_run) { Print() << "P3: v3--v6: add vIP to F3, F4 and F7 :: " << vIP[1] << "\n"; }
752  else
753 #endif
754  {
755  m_F3.add_vertex(vIP);
756  m_F4.add_vertex(vIP);
757  m_F7.add_vertex(vIP);
758  }
759 
760  if ( almostEqual(distIP, 0.0)) {
761 #ifndef AMREX_USE_GPU
762  if (a_dry_run) { Print() << "P3: v3--v6: intersection ~ 0 :: add vIP to F5\n"; }
763  else
764 #endif
765  { m_F5.add_vertex(vIP); }
766 
767  } else if ( almostEqual(distIP, 1.0)) {
768 #ifndef AMREX_USE_GPU
769  if (a_dry_run) { Print() << "P3: v3--v6: intersection ~ 1 :: add vIP to F6\n"; }
770  else
771 #endif
772  { m_F6.add_vertex(vIP); }
773  }
774 
775  }
776 
777  if (cuts%2 == 0) {
778 
779 #ifndef AMREX_USE_GPU
780  if (a_dry_run) { Print() << "P3: Add v6 to F3, F4 and F6\n"; }
781  else
782 #endif
783  {
784  m_F3.add_vertex(v6);
785  m_F4.add_vertex(v6);
786  m_F6.add_vertex(v6);
787  }
788  }
789 
790  if (utils::intersect_plane_edge(m_eb_point, m_eb_normal, v6, v7, vIP, distIP)) {
791 
792  ++cuts;
793 
794 #ifndef AMREX_USE_GPU
795  if (a_dry_run) { Print() << "P3: v6--v7: add vIP to F3, F6 and F7 :: " << vIP[0] << "\n"; }
796  else
797 #endif
798  {
799  m_F3.add_vertex(vIP);
800  m_F6.add_vertex(vIP);
801  m_F7.add_vertex(vIP);
802  }
803 
804  if ( almostEqual(distIP, 0.0)) {
805 #ifndef AMREX_USE_GPU
806  if (a_dry_run) { Print() << "P3: v6--v7: intersection ~ 0 :: add vIP to F4\n"; }
807  else
808 #endif
809  { m_F4.add_vertex(vIP); }
810 
811  } else if ( almostEqual(distIP, 1.0)) {
812 #ifndef AMREX_USE_GPU
813  if (a_dry_run) { Print() << "P3: v6--v7: intersection ~ 1 :: add vIP to F2\n"; }
814  else
815 #endif
816  { m_F2.add_vertex(vIP); }
817  }
818  }
819 
820  if (cuts == 2 && add_v7) {
821 
822 #ifndef AMREX_USE_GPU
823  if (a_dry_run) { Print() << "P1: Add v7 to F2, F3 and F6\n"; }
824  else
825 #endif
826  {
827  m_F2.add_vertex(v7);
828  m_F3.add_vertex(v7);
829  m_F6.add_vertex(v7);
830  }
831  add_v7 = 0;
832  }
833 
834  } // end Path 3
835 
836  if (utils::intersect_plane_edge(m_eb_point, m_eb_normal, v3, v4, vIP, distIP)) {
837 
838 #ifndef AMREX_USE_GPU
839  if (a_dry_run) { Print() << "P6: v3--v4: add vIP to F3, F5 and F7 :: " << vIP[0] << "\n"; }
840  else
841 #endif
842  {
843  m_F3.add_vertex(vIP);
844  m_F5.add_vertex(vIP);
845  m_F7.add_vertex(vIP);
846  }
847 
848  if ( almostEqual(distIP, 0.0)) {
849 #ifndef AMREX_USE_GPU
850  if (a_dry_run) { Print() << "P6: v3--v4: intersection ~ 0 :: add vIP to F4\n"; }
851  else
852 #endif
853  { m_F4.add_vertex(vIP); }
854 
855  } else if ( almostEqual(distIP, 1.0)) {
856 #ifndef AMREX_USE_GPU
857  if (a_dry_run) { Print() << "P6: v3--v4: intersection ~ 1 :: add vIP to F2\n"; }
858  else
859 #endif
860  { m_F2.add_vertex(vIP); }
861  }
862  }
863 }
864 
865 
866 
867 AMREX_GPU_HOST_DEVICE
868 AMREX_FORCE_INLINE
869 void
870 eb_cut_cell_::
871 set_covered_regular_cell_vertices ()
872 {
873  using namespace amrex;
874 
875  RealVect v0(m_rbox.lo(0), m_rbox.lo(1), m_rbox.lo(2));
876  RealVect v1(m_rbox.hi(0), m_rbox.lo(1), m_rbox.lo(2));
877  RealVect v2(m_rbox.lo(0), m_rbox.hi(1), m_rbox.lo(2));
878  RealVect v3(m_rbox.lo(0), m_rbox.lo(1), m_rbox.hi(2));
879  RealVect v4(m_rbox.hi(0), m_rbox.lo(1), m_rbox.hi(2));
880  RealVect v5(m_rbox.hi(0), m_rbox.hi(1), m_rbox.lo(2));
881  RealVect v6(m_rbox.lo(0), m_rbox.hi(1), m_rbox.hi(2));
882  RealVect v7(m_rbox.hi(0), m_rbox.hi(1), m_rbox.hi(2));
883 
884  // Add vertices in the order of outward normal vector
885 
886  m_F1.add_vertex(v0);
887  m_F1.add_vertex(v2);
888  m_F1.add_vertex(v5);
889  m_F1.add_vertex(v1);
890 
891  m_F2.add_vertex(v1);
892  m_F2.add_vertex(v5);
893  m_F2.add_vertex(v7);
894  m_F2.add_vertex(v4);
895 
896  m_F3.add_vertex(v3);
897  m_F3.add_vertex(v4);
898  m_F3.add_vertex(v7);
899  m_F3.add_vertex(v6);
900 
901  m_F4.add_vertex(v0);
902  m_F4.add_vertex(v3);
903  m_F4.add_vertex(v6);
904  m_F4.add_vertex(v2);
905 
906  m_F5.add_vertex(v0);
907  m_F5.add_vertex(v1);
908  m_F5.add_vertex(v4);
909  m_F5.add_vertex(v3);
910 
911  m_F6.add_vertex(v2);
912  m_F6.add_vertex(v5);
913  m_F6.add_vertex(v7);
914  m_F6.add_vertex(v6);
915 
916 }
917 
918 
919 
920 #endif
ERF_EBPolygon.H
ERF_EBUtils.H
eb_cut_cell_
Definition: ERF_EBCutCell.H:14
eb_cut_cell_::volume
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real volume() const
Definition: ERF_EBCutCell.H:57
eb_cut_cell_::m_F6
polygon_ m_F6
Definition: ERF_EBCutCell.H:209
eb_cut_cell_::m_lo_faces
amrex::Array< polygon_ const *const, 3 > m_lo_faces
Definition: ERF_EBCutCell.H:219
eb_cut_cell_::m_rbox_area
amrex::RealVect m_rbox_area
Definition: ERF_EBCutCell.H:198
eb_cut_cell_::isRegular
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool isRegular() const noexcept
Definition: ERF_EBCutCell.H:28
eb_cut_cell_::m_cell_cent
amrex::RealVect m_cell_cent
Definition: ERF_EBCutCell.H:213
eb_cut_cell_::m_dx
amrex::RealVect m_dx
Definition: ERF_EBCutCell.H:214
eb_cut_cell_::areaHi
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real areaHi(int const a_idim) const noexcept
Definition: ERF_EBCutCell.H:87
eb_cut_cell_::areaLo
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real areaLo(int const a_idim) const noexcept
Definition: ERF_EBCutCell.H:78
eb_cut_cell_::m_eb_point
amrex::RealVect const m_eb_point
Definition: ERF_EBCutCell.H:193
eb_cut_cell_::m_invert
amrex::Real m_invert
Definition: ERF_EBCutCell.H:196
eb_cut_cell_::m_rbox
amrex::RealBox const m_rbox
Definition: ERF_EBCutCell.H:192
eb_cut_cell_::centBoun
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::RealVect centBoun() const noexcept
Definition: ERF_EBCutCell.H:141
eb_cut_cell_::set_covered_regular_cell_vertices
AMREX_GPU_HOST_DEVICE void set_covered_regular_cell_vertices()
Definition: ERF_EBCutCell.H:871
eb_cut_cell_::set_regular
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void set_regular()
Definition: ERF_EBCutCell.H:39
eb_cut_cell_::m_F1
polygon_ m_F1
Definition: ERF_EBCutCell.H:204
eb_cut_cell_::areaBoun
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real areaBoun() const noexcept
Definition: ERF_EBCutCell.H:96
eb_cut_cell_::m_F3
polygon_ m_F3
Definition: ERF_EBCutCell.H:206
eb_cut_cell_::m_flag
amrex::EBCellFlag m_flag
Definition: ERF_EBCutCell.H:200
eb_cut_cell_::calc_edge_intersections
AMREX_GPU_HOST_DEVICE void calc_edge_intersections(int const a_dry_run=0)
Definition: ERF_EBCutCell.H:342
eb_cut_cell_::m_F5
polygon_ m_F5
Definition: ERF_EBCutCell.H:208
eb_cut_cell_::centLo
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::RealVect centLo(int const a_idim) const noexcept
Definition: ERF_EBCutCell.H:101
eb_cut_cell_::centVol
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::RealVect centVol() const noexcept
Definition: ERF_EBCutCell.H:159
eb_cut_cell_::set_covered
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void set_covered()
Definition: ERF_EBCutCell.H:34
eb_cut_cell_::m_F7
polygon_ m_F7
Definition: ERF_EBCutCell.H:217
eb_cut_cell_::normBoun
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::RealVect normBoun() const noexcept
Definition: ERF_EBCutCell.H:150
eb_cut_cell_::m_hi_faces_id
amrex::Array< int const, 3 > m_hi_faces_id
Definition: ERF_EBCutCell.H:223
eb_cut_cell_::eb_cut_cell_
AMREX_GPU_HOST_DEVICE eb_cut_cell_(amrex::EBCellFlag const &a_flag, amrex::RealBox const &a_rbox, amrex::RealVect const &a_point, amrex::RealVect const &a_normal)
Definition: ERF_EBCutCell.H:236
eb_cut_cell_::isSingleValued
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool isSingleValued() const noexcept
Definition: ERF_EBCutCell.H:31
eb_cut_cell_::isCovered
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool isCovered() const noexcept
Definition: ERF_EBCutCell.H:25
eb_cut_cell_::debug
void debug(int const a_face=-1)
Definition: ERF_EBCutCell.cpp:9
eb_cut_cell_::centHi
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::RealVect centHi(int const a_idim) const noexcept
Definition: ERF_EBCutCell.H:121
eb_cut_cell_::m_max_faces
static constexpr int m_max_faces
Definition: ERF_EBCutCell.H:211
eb_cut_cell_::m_eb_normal
amrex::RealVect const m_eb_normal
Definition: ERF_EBCutCell.H:194
eb_cut_cell_::m_lo_faces_id
amrex::Array< int const, 3 > m_lo_faces_id
Definition: ERF_EBCutCell.H:222
eb_cut_cell_::m_F2
polygon_ m_F2
Definition: ERF_EBCutCell.H:205
eb_cut_cell_::m_hi_faces
amrex::Array< polygon_ const *const, 3 > m_hi_faces
Definition: ERF_EBCutCell.H:220
eb_cut_cell_::m_F4
polygon_ m_F4
Definition: ERF_EBCutCell.H:207
eb_cut_cell_::m_cellface_cent
amrex::Array< amrex::RealVect, m_max_faces > m_cellface_cent
Definition: ERF_EBCutCell.H:212
polygon_
Definition: ERF_EBPolygon.H:9
polygon_::get_centroid
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::RealVect get_centroid() const noexcept
Definition: ERF_EBPolygon.H:138
polygon_::ok
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int ok() const noexcept
Definition: ERF_EBPolygon.H:119
polygon_::add_vertex
AMREX_GPU_HOST_DEVICE void add_vertex(amrex::RealVect const &a_v)
Definition: ERF_EBPolygon.H:43
polygon_::distance
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real distance(amrex::RealVect const &a_point) const noexcept
Definition: ERF_EBPolygon.H:130
polygon_::area
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real area() const noexcept
Definition: ERF_EBPolygon.H:123
polygon_::set_area
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void set_area(amrex::Real const &a_area)
Definition: ERF_EBPolygon.H:50
amrex
Definition: ERF_ConsoleIO.cpp:12
utils::intersect_plane_edge
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int intersect_plane_edge(amrex::RealVect const &a_plane_point, amrex::RealVect const &a_plane_normal, amrex::RealVect const &a_edge_point0, amrex::RealVect const &a_edge_point1, amrex::RealVect &a_intersection_point, amrex::Real &a_intersection_dist)
Definition: ERF_EBUtils.H:30