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