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 #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 ();
227 
228  AMREX_GPU_HOST_DEVICE
229  void set_covered_regular_cell_vertices ();
230 };
231 
232 AMREX_GPU_HOST_DEVICE
233 AMREX_FORCE_INLINE
234 eb_cut_cell_::
235 eb_cut_cell_ ( amrex::EBCellFlag const& a_flag,
236  amrex::RealBox const& a_rbox,
237  amrex::RealVect const& a_point,
238  amrex::RealVect const& a_normal )
239  : m_rbox(a_rbox)
240  , m_eb_point(a_point)
241  , m_eb_normal(a_normal)
242  , m_invert(0.0)
243  , m_F1(a_point, a_normal)
244  , m_F2(a_point, a_normal)
245  , m_F3(a_point, a_normal)
246  , m_F4(a_point, a_normal)
247  , m_F5(a_point, a_normal)
248  , m_F6(a_point, a_normal)
249  , m_cellface_cent({amrex::RealVect(0.), amrex::RealVect(0.), amrex::RealVect(0.),
250  amrex::RealVect(0.), amrex::RealVect(0.), amrex::RealVect(0.)})
251 {
252  using namespace amrex;
253 
254  m_rbox_area[0] = m_rbox.length(1)*m_rbox.length(2);
255  m_rbox_area[1] = m_rbox.length(0)*m_rbox.length(2);
256  m_rbox_area[2] = m_rbox.length(0)*m_rbox.length(1);
257 
258  RealVect v0(m_rbox.lo(0), m_rbox.lo(1), m_rbox.lo(2));
259  RealVect v1(m_rbox.hi(0), m_rbox.lo(1), m_rbox.lo(2));
260  RealVect v2(m_rbox.lo(0), m_rbox.hi(1), m_rbox.lo(2));
261  RealVect v3(m_rbox.lo(0), m_rbox.lo(1), m_rbox.hi(2));
262  RealVect v4(m_rbox.hi(0), m_rbox.lo(1), m_rbox.hi(2));
263  RealVect v5(m_rbox.hi(0), m_rbox.hi(1), m_rbox.lo(2));
264  RealVect v6(m_rbox.lo(0), m_rbox.hi(1), m_rbox.hi(2));
265  RealVect v7(m_rbox.hi(0), m_rbox.hi(1), m_rbox.hi(2));
266 
267  // Centoids of cell faces
268 
269  m_cellface_cent[0] = 0.25 * ( v0 + v2 + v5 + v1 ); // F1
270  m_cellface_cent[1] = 0.25 * ( v1 + v5 + v7 + v4 ); // F2
271  m_cellface_cent[2] = 0.25 * ( v3 + v4 + v7 + v6 ); // F3
272  m_cellface_cent[3] = 0.25 * ( v0 + v3 + v6 + v2 ); // F4
273  m_cellface_cent[4] = 0.25 * ( v0 + v1 + v4 + v3 ); // F5
274  m_cellface_cent[5] = 0.25 * ( v2 + v5 + v7 + v6 ); // F6
275 
276  // Cell centroid
277 
278  m_cell_cent[0] = 0.5 * ( m_rbox.lo(0) + m_rbox.hi(0) );
279  m_cell_cent[1] = 0.5 * ( m_rbox.lo(1) + m_rbox.hi(1) );
280  m_cell_cent[2] = 0.5 * ( m_rbox.lo(2) + m_rbox.hi(2) );
281 
282  // Cell size
283 
284  m_dx[0] = m_rbox.hi(0) - m_rbox.lo(0);
285  m_dx[1] = m_rbox.hi(1) - m_rbox.lo(1);
286  m_dx[2] = m_rbox.hi(2) - m_rbox.lo(2);
287 
288  if (a_flag.isCovered() ) {
289 
290  set_covered();
291 
292  } else if (a_flag.isRegular() ) {
293 
294  set_regular();
295 
296  } else { // Check that the box and plane intersect.
297 
298  RealVect c = 0.5*(v0 + v7);
299  RealVect e = v7 - c;
300 
301  Real r = e[0]*amrex::Math::abs(a_normal[0]) +
302  e[1]*amrex::Math::abs(a_normal[1]) +
303  e[2]*amrex::Math::abs(a_normal[2]);
304 
305  Real s = amrex::Math::abs(c.dotProduct(a_normal)
306  - a_point.dotProduct(a_normal));
307 
308  if (s > r) {
309  if ((a_normal.dotProduct(v0) - a_normal.dotProduct(a_point)) > 0.)
310  { set_covered(); } else { set_regular(); }
311  } else { m_flag.setSingleValued(); }
312  }
313 
314  if ( m_flag.isSingleValued() ) {
315 
316  m_invert = ((m_eb_normal.dotProduct(v0) - m_eb_normal.dotProduct(m_eb_point)) > 0.) ? 0.0 : 1.0;
317 
318  calc_edge_intersections();
319 
320  } // end singleValued
321 
322  m_F1.define();
323  m_F2.define();
324  m_F3.define();
325  m_F4.define();
326  m_F5.define();
327  m_F6.define();
328  m_F7.define();
329 
330  // For covered and regular cut cells, add vertices to utilize e.g., get_centroid.
331  if ( !m_flag.isSingleValued() ) {
332  set_covered_regular_cell_vertices();
333  }
334 
335 }
336 
337 AMREX_GPU_HOST_DEVICE
338 AMREX_FORCE_INLINE
339 void
340 eb_cut_cell_::
341 calc_edge_intersections ()
342 {
343  using namespace amrex;
344 
345  RealVect v0(m_rbox.lo(0), m_rbox.lo(1), m_rbox.lo(2));
346  RealVect v1(m_rbox.hi(0), m_rbox.lo(1), m_rbox.lo(2));
347  RealVect v2(m_rbox.lo(0), m_rbox.hi(1), m_rbox.lo(2));
348  RealVect v3(m_rbox.lo(0), m_rbox.lo(1), m_rbox.hi(2));
349  RealVect v4(m_rbox.hi(0), m_rbox.lo(1), m_rbox.hi(2));
350  RealVect v5(m_rbox.hi(0), m_rbox.hi(1), m_rbox.lo(2));
351  RealVect v6(m_rbox.lo(0), m_rbox.hi(1), m_rbox.hi(2));
352  RealVect v7(m_rbox.hi(0), m_rbox.hi(1), m_rbox.hi(2));
353 
354 #ifndef AMREX_USE_GPU
355  constexpr bool print_F1 = false;
356  constexpr bool print_F2 = false;
357  constexpr bool print_F3 = false;
358  constexpr bool print_F4 = false;
359  constexpr bool print_F5 = false;
360  constexpr bool print_F6 = false;
361  constexpr bool print_F7 = false;
362 #endif
363 
364 #ifdef AMREX_USE_GPU
365  m_F1.add_vertex(v0);
366  m_F4.add_vertex(v0);
367  m_F5.add_vertex(v0);
368 #else
369  m_F1.add_vertex(v0,"Initial: add v0 -> F1",print_F1);
370  m_F4.add_vertex(v0,"Initial: add v0 -> F4",print_F4);
371  m_F5.add_vertex(v0,"Initial: add v0 -> F5",print_F5);
372 #endif
373 
374  int add_v7(1);
375 
376  RealVect vIP;
377  Real distIP;
378  constexpr Real tol = 1.e-15;
379 
380  amrex::Array<bool, 8> vertex_intersected{}; // true if vertex has been intersected within a path.
381 
382  // Path 1
383  {
384  int cuts(0);
385  for (int i = 0; i < 8; ++i) vertex_intersected[i] = false;
386 
387  if (utils::intersect_plane_edge(m_eb_point, m_eb_normal, v0, v1, vIP, distIP)) {
388 
389  amrex::Real length = (v0-v1).vectorLength();
390 
391  if ( Math::abs(distIP) < tol ) {
392 #ifdef AMREX_USE_GPU
393  m_F1.add_vertex(v0);
394  m_F4.add_vertex(v0);
395  m_F5.add_vertex(v0);
396  m_F7.add_vertex(v0);
397 #else
398  m_F1.add_vertex(v0,"Path 1: v0--v1: add v0 -> F1",print_F1);
399  m_F4.add_vertex(v0,"Path 1: v0--v1: add v0 -> F4",print_F4);
400  m_F5.add_vertex(v0,"Path 1: v0--v1: add v0 -> F5",print_F5);
401  m_F7.add_vertex(v0,"Path 1: v0--v1: add v0 -> F7",print_F7);
402 #endif
403  if (!vertex_intersected[0]) {
404  vertex_intersected[0] = true;
405  ++cuts;
406  }
407  } else if ( Math::abs(distIP - length) < tol ) {
408 #ifdef AMREX_USE_GPU
409  m_F1.add_vertex(v1);
410  m_F2.add_vertex(v1);
411  m_F5.add_vertex(v1);
412  m_F7.add_vertex(v1);
413 #else
414  m_F1.add_vertex(v1,"Path 1: v0--v1: add v1 -> F1",print_F1);
415  m_F2.add_vertex(v1,"Path 1: v0--v1: add v1 -> F2",print_F2);
416  m_F5.add_vertex(v1,"Path 1: v0--v1: add v1 -> F5",print_F5);
417  m_F7.add_vertex(v1,"Path 1: v0--v1: add v1 -> F7",print_F7);
418 #endif
419  if (!vertex_intersected[1]) {
420  vertex_intersected[1] = true;
421  ++cuts;
422  }
423  } else {
424 #ifdef AMREX_USE_GPU
425  m_F1.add_vertex(vIP);
426  m_F5.add_vertex(vIP);
427  m_F7.add_vertex(vIP);
428 #else
429  m_F1.add_vertex(vIP,"Path 1: v0--v1: add vIP -> F1",print_F1);
430  m_F5.add_vertex(vIP,"Path 1: v0--v1: add vIP -> F5",print_F5);
431  m_F7.add_vertex(vIP,"Path 1: v0--v1: add vIP -> F7",print_F7);
432 #endif
433  ++cuts;
434  }
435 
436  } // Path 1: v0--v1
437 
438  if (cuts%2 == 0) {
439 #ifdef AMREX_USE_GPU
440  m_F1.add_vertex(v1);
441  m_F2.add_vertex(v1);
442  m_F5.add_vertex(v1);
443 #else
444  m_F1.add_vertex(v1,"Path 1: after v0--v1, cuts%2 == 0: add v1 -> F1",print_F1);
445  m_F2.add_vertex(v1,"Path 1: after v0--v1, cuts%2 == 0: add v1 -> F2",print_F2);
446  m_F5.add_vertex(v1,"Path 1: after v0--v1, cuts%2 == 0: add v1 -> F5",print_F5);
447 #endif
448  }
449 
450  if (utils::intersect_plane_edge(m_eb_point, m_eb_normal, v1, v4, vIP, distIP)) {
451 
452  amrex::Real length = (v1-v4).vectorLength();
453 
454  if ( Math::abs(distIP) < tol ) {
455 #ifdef AMREX_USE_GPU
456  m_F1.add_vertex(v1);
457  m_F2.add_vertex(v1);
458  m_F5.add_vertex(v1);
459  m_F7.add_vertex(v1);
460 #else
461  m_F1.add_vertex(v1,"Path 1: v1--v4: add v1 -> F1",print_F1);
462  m_F2.add_vertex(v1,"Path 1: v1--v4: add v1 -> F2",print_F2);
463  m_F5.add_vertex(v1,"Path 1: v1--v4: add v1 -> F5",print_F5);
464  m_F7.add_vertex(v1,"Path 1: v1--v4: add v1 -> F7",print_F7);
465 #endif
466  if (!vertex_intersected[1]) {
467  vertex_intersected[1] = true;
468  ++cuts;
469  }
470  } else if ( Math::abs(distIP - length) < tol ) {
471 #ifdef AMREX_USE_GPU
472  m_F2.add_vertex(v4);
473  m_F3.add_vertex(v4);
474  m_F5.add_vertex(v4);
475  m_F7.add_vertex(v4);
476 #else
477  m_F2.add_vertex(v4,"Path 1: v1--v4: add v4 -> F2",print_F2);
478  m_F3.add_vertex(v4,"Path 1: v1--v4: add v4 -> F3",print_F3);
479  m_F5.add_vertex(v4,"Path 1: v1--v4: add v4 -> F5",print_F5);
480  m_F7.add_vertex(v4,"Path 1: v1--v4: add v4 -> F7",print_F7);
481 #endif
482  if (!vertex_intersected[4]) {
483  vertex_intersected[4] = true;
484  ++cuts;
485  }
486  } else {
487 #ifdef AMREX_USE_GPU
488  m_F2.add_vertex(vIP);
489  m_F5.add_vertex(vIP);
490  m_F7.add_vertex(vIP);
491 #else
492  m_F2.add_vertex(vIP,"Path 1: v1--v4: add vIP -> F2",print_F2);
493  m_F5.add_vertex(vIP,"Path 1: v1--v4: add vIP -> F5",print_F5);
494  m_F7.add_vertex(vIP,"Path 1: v1--v4: add vIP -> F7",print_F7);
495 #endif
496  ++cuts;
497  }
498 
499  } // P1: v1--v4
500 
501  if (cuts%2 == 0) {
502 #ifdef AMREX_USE_GPU
503  m_F2.add_vertex(v4);
504  m_F3.add_vertex(v4);
505  m_F5.add_vertex(v4);
506 #else
507  m_F2.add_vertex(v4,"Path 1: after v1--v4, cuts%2 == 0: add v4 -> F2",print_F2);
508  m_F3.add_vertex(v4,"Path 1: after v1--v4, cuts%2 == 0: add v4 -> F3",print_F3);
509  m_F5.add_vertex(v4,"Path 1: after v1--v4, cuts%2 == 0: add v4 -> F5",print_F5);
510 #endif
511  }
512 
513  if (utils::intersect_plane_edge(m_eb_point, m_eb_normal, v4, v7, vIP, distIP)) {
514 
515  amrex::Real length = (v4-v7).vectorLength();
516 
517  if ( Math::abs(distIP) < tol) {
518 #ifdef AMREX_USE_GPU
519  m_F2.add_vertex(v4);
520  m_F3.add_vertex(v4);
521  m_F5.add_vertex(v4);
522  m_F7.add_vertex(v4);
523 #else
524  m_F2.add_vertex(v4,"Path 1: v4--v7: add v4 -> F2",print_F2);
525  m_F3.add_vertex(v4,"Path 1: v4--v7: add v4 -> F3",print_F3);
526  m_F5.add_vertex(v4,"Path 1: v4--v7: add v4 -> F5",print_F5);
527  m_F7.add_vertex(v4,"Path 1: v4--v7: add v4 -> F7",print_F7);
528 #endif
529  if (!vertex_intersected[4]) {
530  vertex_intersected[4] = true;
531  ++cuts;
532  }
533  } else if ( Math::abs(distIP - length) < tol) {
534 #ifdef AMREX_USE_GPU
535  m_F2.add_vertex(v7);
536  m_F3.add_vertex(v7);
537  m_F6.add_vertex(v7);
538  m_F7.add_vertex(v7);
539 #else
540  m_F2.add_vertex(v7,"Path 1: v4--v7: add v7 -> F2",print_F2);
541  m_F3.add_vertex(v7,"Path 1: v4--v7: add v7 -> F3",print_F3);
542  m_F6.add_vertex(v7,"Path 1: v4--v7: add v7 -> F6",print_F6);
543  m_F7.add_vertex(v7,"Path 1: v4--v7: add v7 -> F7",print_F7);
544 #endif
545  if (!vertex_intersected[7]) {
546  vertex_intersected[7] = true;
547  ++cuts;
548  }
549  } else {
550 #ifdef AMREX_USE_GPU
551  m_F2.add_vertex(vIP);
552  m_F3.add_vertex(vIP);
553  m_F7.add_vertex(vIP);
554 #else
555  m_F2.add_vertex(vIP,"Path 1: v4--v7: add vIP -> F2",print_F2);
556  m_F3.add_vertex(vIP,"Path 1: v4--v7: add vIP -> F3",print_F3);
557  m_F7.add_vertex(vIP,"Path 1: v4--v7: add vIP -> F7",print_F7);
558 #endif
559  ++cuts;
560  }
561 
562  } // P1: v4--v7
563 
564  if (cuts == 2 && add_v7) {
565  if (!utils::intersect_plane_edge(m_eb_point, m_eb_normal, v0, v7, vIP, distIP)) {
566 #ifdef AMREX_USE_GPU
567  m_F2.add_vertex(v7);
568  m_F3.add_vertex(v7);
569  m_F6.add_vertex(v7);
570 #else
571  m_F2.add_vertex(v7,"Path 1: after v4--v7, cuts == 2 && add_v7: add v7 -> F2",print_F2);
572  m_F3.add_vertex(v7,"Path 1: after v4--v7, cuts == 2 && add_v7: add v7 -> F3",print_F3);
573  m_F6.add_vertex(v7,"Path 1: after v4--v7, cuts == 2 && add_v7: add v7 -> F6",print_F6);
574 #endif
575  add_v7 = 0;
576  }
577  }
578  } // end Path 1
579 
580  // Path 4
581  if (utils::intersect_plane_edge(m_eb_point, m_eb_normal, v1, v5, vIP, distIP)) {
582 
583  amrex::Real length = (v1-v5).vectorLength();
584 
585  if ( Math::abs(distIP) < tol) {
586 #ifdef AMREX_USE_GPU
587  m_F1.add_vertex(v1);
588  m_F2.add_vertex(v1);
589  m_F5.add_vertex(v1);
590  m_F7.add_vertex(v1);
591 #else
592  m_F1.add_vertex(v1,"Path 4: v1--v5: add v1 -> F1",print_F1);
593  m_F2.add_vertex(v1,"Path 4: v1--v5: add v1 -> F2",print_F2);
594  m_F5.add_vertex(v1,"Path 4: v1--v5: add v1 -> F5",print_F5);
595  m_F7.add_vertex(v1,"Path 4: v1--v5: add v1 -> F7",print_F7);
596 #endif
597  } else if ( Math::abs(distIP - length) < tol) {
598 #ifdef AMREX_USE_GPU
599  m_F1.add_vertex(v5);
600  m_F2.add_vertex(v5);
601  m_F6.add_vertex(v5);
602  m_F7.add_vertex(v5);
603 #else
604  m_F1.add_vertex(v5,"Path 4: v1--v5: add v5 -> F1",print_F1);
605  m_F2.add_vertex(v5,"Path 4: v1--v5: add v5 -> F2",print_F2);
606  m_F6.add_vertex(v5,"Path 4: v1--v5: add v5 -> F6",print_F6);
607  m_F7.add_vertex(v5,"Path 4: v1--v5: add v5 -> F7",print_F7);
608 #endif
609  } else {
610 #ifdef AMREX_USE_GPU
611  m_F1.add_vertex(vIP);
612  m_F2.add_vertex(vIP);
613  m_F7.add_vertex(vIP);
614 #else
615  m_F1.add_vertex(vIP,"Path 4: v1--v5: add vIP -> F1",print_F1);
616  m_F2.add_vertex(vIP,"Path 4: v1--v5: add vIP -> F2",print_F2);
617  m_F7.add_vertex(vIP,"Path 4: v1--v5: add vIP -> F7",print_F7);
618 #endif
619  }
620  }
621 
622  // Path 2
623  { int cuts(0);
624  for (int i = 0; i < 8; ++i) vertex_intersected[i] = false;
625 
626  if (utils::intersect_plane_edge(m_eb_point, m_eb_normal, v0, v2, vIP, distIP)) {
627 
628  amrex::Real length = (v0-v2).vectorLength();
629 
630  if ( Math::abs(distIP) < tol ) {
631 #ifdef AMREX_USE_GPU
632  m_F1.add_vertex(v0);
633  m_F4.add_vertex(v0);
634  m_F5.add_vertex(v0);
635  m_F7.add_vertex(v0);
636 #else
637  m_F1.add_vertex(v0,"Path 2: v0--v2: add v0 -> F1",print_F1);
638  m_F4.add_vertex(v0,"Path 2: v0--v2: add v0 -> F4",print_F4);
639  m_F5.add_vertex(v0,"Path 2: v0--v2: add v0 -> F5",print_F5);
640  m_F7.add_vertex(v0,"Path 2: v0--v2: add v0 -> F7",print_F7);
641 #endif
642  if (!vertex_intersected[0]) {
643  vertex_intersected[0] = true;
644  ++cuts;
645  }
646  } else if ( Math::abs(distIP - length) < tol ) {
647 #ifdef AMREX_USE_GPU
648  m_F1.add_vertex(v2);
649  m_F4.add_vertex(v2);
650  m_F6.add_vertex(v2);
651  m_F7.add_vertex(v2);
652 #else
653  m_F1.add_vertex(v2,"Path 2: v0--v2: add v2 -> F1",print_F1);
654  m_F4.add_vertex(v2,"Path 2: v0--v2: add v2 -> F4",print_F4);
655  m_F6.add_vertex(v2,"Path 2: v0--v2: add v2 -> F6",print_F6);
656  m_F7.add_vertex(v2,"Path 2: v0--v2: add v2 -> F7",print_F7);
657 #endif
658  if (!vertex_intersected[2]) {
659  vertex_intersected[2] = true;
660  ++cuts;
661  }
662  } else {
663 #ifdef AMREX_USE_GPU
664  m_F1.add_vertex(vIP);
665  m_F4.add_vertex(vIP);
666  m_F7.add_vertex(vIP);
667 #else
668  m_F1.add_vertex(vIP,"Path 2: v0--v2: add vIP -> F1",print_F1);
669  m_F4.add_vertex(vIP,"Path 2: v0--v2: add vIP -> F4",print_F4);
670  m_F7.add_vertex(vIP,"Path 2: v0--v2: add vIP -> F7",print_F7);
671 #endif
672  ++cuts;
673  }
674 
675  } // P2: v0--v2
676 
677  if (cuts%2 == 0) {
678 #ifdef AMREX_USE_GPU
679  m_F1.add_vertex(v2);
680  m_F4.add_vertex(v2);
681  m_F6.add_vertex(v2);
682 #else
683  m_F1.add_vertex(v2,"Path 2: after v0--v2, cuts%2 == 0: add v2 -> F1",print_F1);
684  m_F4.add_vertex(v2,"Path 2: after v0--v2, cuts%2 == 0: add v2 -> F4",print_F4);
685  m_F6.add_vertex(v2,"Path 2: after v0--v2, cuts%2 == 0: add v2 -> F6",print_F6);
686 #endif
687  }
688 
689  if (utils::intersect_plane_edge(m_eb_point, m_eb_normal, v2, v5, vIP, distIP)) {
690 
691  amrex::Real length = (v2-v5).vectorLength();
692 
693  if ( Math::abs(distIP) < tol) {
694 #ifdef AMREX_USE_GPU
695  m_F1.add_vertex(v2);
696  m_F4.add_vertex(v2);
697  m_F6.add_vertex(v2);
698  m_F7.add_vertex(v2);
699 #else
700  m_F1.add_vertex(v2,"Path 2: v2--v5: add v2 -> F1",print_F1);
701  m_F4.add_vertex(v2,"Path 2: v2--v5: add v2 -> F4",print_F4);
702  m_F6.add_vertex(v2,"Path 2: v2--v5: add v2 -> F6",print_F6);
703  m_F7.add_vertex(v2,"Path 2: v2--v5: add v2 -> F7",print_F7);
704 #endif
705  if (!vertex_intersected[2]) {
706  vertex_intersected[2] = true;
707  ++cuts;
708  }
709  } else if ( Math::abs(distIP - length) < tol) {
710 #ifdef AMREX_USE_GPU
711  m_F1.add_vertex(v5);
712  m_F2.add_vertex(v5);
713  m_F6.add_vertex(v5);
714  m_F7.add_vertex(v5);
715 #else
716  m_F1.add_vertex(v5,"Path 2: v2--v5: add v5 -> F1",print_F1);
717  m_F2.add_vertex(v5,"Path 2: v2--v5: add v5 -> F2",print_F2);
718  m_F6.add_vertex(v5,"Path 2: v2--v5: add v5 -> F6",print_F6);
719  m_F7.add_vertex(v5,"Path 2: v2--v5: add v5 -> F7",print_F7);
720 #endif
721  if (!vertex_intersected[5]) {
722  vertex_intersected[5] = true;
723  ++cuts;
724  }
725  } else {
726 #ifdef AMREX_USE_GPU
727  m_F1.add_vertex(vIP);
728  m_F6.add_vertex(vIP);
729  m_F7.add_vertex(vIP);
730 #else
731  m_F1.add_vertex(vIP,"Path 2: v2--v5: add vIP -> F1",print_F1);
732  m_F6.add_vertex(vIP,"Path 2: v2--v5: add vIP -> F6",print_F6);
733  m_F7.add_vertex(vIP,"Path 2: v2--v5: add vIP -> F7",print_F7);
734 #endif
735  ++cuts;
736  }
737  }
738 
739  if (cuts%2 == 0) {
740 #ifdef AMREX_USE_GPU
741  m_F1.add_vertex(v5);
742  m_F2.add_vertex(v5);
743  m_F6.add_vertex(v5);
744 #else
745  m_F1.add_vertex(v5,"Path 2: after v2--v5, cuts%2 == 0: add v5 -> F1",print_F1);
746  m_F2.add_vertex(v5,"Path 2: after v2--v5, cuts%2 == 0: add v5 -> F2",print_F2);
747  m_F6.add_vertex(v5,"Path 2: after v2--v5, cuts%2 == 0: add v5 -> F6",print_F6);
748 #endif
749  }
750 
751  if (utils::intersect_plane_edge(m_eb_point, m_eb_normal, v5, v7, vIP, distIP)) {
752 
753  amrex::Real length = (v5-v7).vectorLength();
754 
755  if ( Math::abs(distIP) < tol) {
756 #ifdef AMREX_USE_GPU
757  m_F1.add_vertex(v5);
758  m_F2.add_vertex(v5);
759  m_F6.add_vertex(v5);
760  m_F7.add_vertex(v5);
761 #else
762  m_F1.add_vertex(v5,"Path 2: v5--v7: add v5 -> F1",print_F1);
763  m_F2.add_vertex(v5,"Path 2: v5--v7: add v5 -> F2",print_F2);
764  m_F6.add_vertex(v5,"Path 2: v5--v7: add v5 -> F6",print_F6);
765  m_F7.add_vertex(v5,"Path 2: v5--v7: add v5 -> F7",print_F7);
766 #endif
767  if (!vertex_intersected[5]) {
768  vertex_intersected[5] = true;
769  ++cuts;
770  }
771  } else if ( Math::abs(distIP - length) < tol) {
772 #ifdef AMREX_USE_GPU
773  m_F2.add_vertex(v7);
774  m_F3.add_vertex(v7);
775  m_F6.add_vertex(v7);
776  m_F7.add_vertex(v7);
777 #else
778  m_F2.add_vertex(v7,"Path 2: v5--v7: add v7 -> F2",print_F2);
779  m_F3.add_vertex(v7,"Path 2: v5--v7: add v7 -> F3",print_F3);
780  m_F6.add_vertex(v7,"Path 2: v5--v7: add v7 -> F6",print_F6);
781  m_F7.add_vertex(v7,"Path 2: v5--v7: add v7 -> F7",print_F7);
782 #endif
783  if (!vertex_intersected[7]) {
784  vertex_intersected[7] = true;
785  ++cuts;
786  }
787  } else {
788 #ifdef AMREX_USE_GPU
789  m_F2.add_vertex(vIP);
790  m_F6.add_vertex(vIP);
791  m_F7.add_vertex(vIP);
792 #else
793  m_F2.add_vertex(vIP,"Path 2: v5--v7: add vIP -> F2",print_F2);
794  m_F6.add_vertex(vIP,"Path 2: v5--v7: add vIP -> F6",print_F6);
795  m_F7.add_vertex(vIP,"Path 2: v5--v7: add vIP -> F7",print_F7);
796 #endif
797  ++cuts;
798  }
799  } // Path 2: v5--v7
800 
801  if (cuts == 2 && add_v7) {
802  if (!utils::intersect_plane_edge(m_eb_point, m_eb_normal, v0, v7, vIP, distIP)) {
803 #ifdef AMREX_USE_GPU
804  m_F2.add_vertex(v7);
805  m_F3.add_vertex(v7);
806  m_F6.add_vertex(v7);
807 #else
808  m_F2.add_vertex(v7,"Path 2: after v5--v7, cuts == 2 && add_v7: add v7 -> F2",print_F2);
809  m_F3.add_vertex(v7,"Path 2: after v5--v7, cuts == 2 && add_v7: add v7 -> F3",print_F3);
810  m_F6.add_vertex(v7,"Path 2: after v5--v7, cuts == 2 && add_v7: add v7 -> F6",print_F6);
811 #endif
812  add_v7 = 0;
813  }
814  }
815 
816  } // end Path 2
817 
818  // Path 5
819  if (utils::intersect_plane_edge(m_eb_point, m_eb_normal, v2, v6, vIP, distIP)) {
820 
821  amrex::Real length = (v2-v6).vectorLength();
822 
823  if ( Math::abs(distIP) < tol) {
824 #ifdef AMREX_USE_GPU
825  m_F1.add_vertex(v2);
826  m_F4.add_vertex(v2);
827  m_F6.add_vertex(v2);
828  m_F7.add_vertex(v2);
829 #else
830  m_F1.add_vertex(v2,"Path 5: v2--v6: add v2 -> F1",print_F1);
831  m_F4.add_vertex(v2,"Path 5: v2--v6: add v2 -> F4",print_F4);
832  m_F6.add_vertex(v2,"Path 5: v2--v6: add v2 -> F6",print_F6);
833  m_F7.add_vertex(v2,"Path 5: v2--v6: add v2 -> F7",print_F7);
834 #endif
835  } else if ( Math::abs(distIP - length) < tol) {
836 #ifdef AMREX_USE_GPU
837  m_F3.add_vertex(v6);
838  m_F4.add_vertex(v6);
839  m_F6.add_vertex(v6);
840  m_F7.add_vertex(v6);
841 #else
842  m_F3.add_vertex(v6,"Path 5: v2--v6: add v6 -> F3",print_F3);
843  m_F4.add_vertex(v6,"Path 5: v2--v6: add v6 -> F4",print_F4);
844  m_F6.add_vertex(v6,"Path 5: v2--v6: add v6 -> F6",print_F6);
845  m_F7.add_vertex(v6,"Path 5: v2--v6: add v6 -> F7",print_F7);
846 #endif
847  } else {
848 #ifdef AMREX_USE_GPU
849  m_F4.add_vertex(vIP);
850  m_F6.add_vertex(vIP);
851  m_F7.add_vertex(vIP);
852 #else
853  m_F4.add_vertex(vIP,"Path 5: v2--v6: add vIP -> F4",print_F4);
854  m_F6.add_vertex(vIP,"Path 5: v2--v6: add vIP -> F6",print_F6);
855  m_F7.add_vertex(vIP,"Path 5: v2--v6: add vIP -> F7",print_F7);
856 #endif
857  }
858  } // end Path 5
859 
860 
861  // Path 3
862  { int cuts(0);
863  for (int i = 0; i < 8; ++i) vertex_intersected[i] = false;
864 
865  if (utils::intersect_plane_edge(m_eb_point, m_eb_normal, v0, v3, vIP, distIP)) {
866 
867  amrex::Real length = (v0-v3).vectorLength();
868 
869  if ( Math::abs(distIP) < tol) {
870 #ifdef AMREX_USE_GPU
871  m_F1.add_vertex(v0);
872  m_F4.add_vertex(v0);
873  m_F5.add_vertex(v0);
874  m_F7.add_vertex(v0);
875 #else
876  m_F1.add_vertex(v0,"Path 3: v0--v3: add v0 -> F1",print_F1);
877  m_F4.add_vertex(v0,"Path 3: v0--v3: add v0 -> F4",print_F4);
878  m_F5.add_vertex(v0,"Path 3: v0--v3: add v0 -> F5",print_F5);
879  m_F7.add_vertex(v0,"Path 3: v0--v3: add v0 -> F7",print_F7);
880 #endif
881  if (!vertex_intersected[0]) {
882  vertex_intersected[0] = true;
883  ++cuts;
884  }
885  } else if ( Math::abs(distIP - length) < tol) {
886 #ifdef AMREX_USE_GPU
887  m_F3.add_vertex(v3);
888  m_F4.add_vertex(v3);
889  m_F5.add_vertex(v3);
890  m_F7.add_vertex(v3);
891 #else
892  m_F3.add_vertex(v3,"Path 3: v0--v3: add v3 -> F3",print_F3);
893  m_F4.add_vertex(v3,"Path 3: v0--v3: add v3 -> F4",print_F4);
894  m_F5.add_vertex(v3,"Path 3: v0--v3: add v3 -> F5",print_F5);
895  m_F7.add_vertex(v3,"Path 3: v0--v3: add v3 -> F7",print_F7);
896 #endif
897  if (!vertex_intersected[3]) {
898  vertex_intersected[3] = true;
899  ++cuts;
900  }
901  } else {
902 #ifdef AMREX_USE_GPU
903  m_F4.add_vertex(vIP);
904  m_F5.add_vertex(vIP);
905  m_F7.add_vertex(vIP);
906 #else
907  m_F4.add_vertex(vIP,"Path 3: v0--v3: add vIP -> F4",print_F4);
908  m_F5.add_vertex(vIP,"Path 3: v0--v3: add vIP -> F5",print_F5);
909  m_F7.add_vertex(vIP,"Path 3: v0--v3: add vIP -> F7",print_F7);
910 #endif
911  ++cuts;
912  }
913 
914  } // P3: v0--v3
915 
916  if (cuts%2 == 0) {
917 #ifdef AMREX_USE_GPU
918  m_F3.add_vertex(v3);
919  m_F4.add_vertex(v3);
920  m_F5.add_vertex(v3);
921 #else
922  m_F3.add_vertex(v3,"Path 3: after v0--v3, cuts%2 == 0: add v3 -> F3",print_F3);
923  m_F4.add_vertex(v3,"Path 3: after v0--v3, cuts%2 == 0: add v3 -> F4",print_F4);
924  m_F5.add_vertex(v3,"Path 3: after v0--v3, cuts%2 == 0: add v3 -> F5",print_F5);
925 #endif
926  }
927 
928  if (utils::intersect_plane_edge(m_eb_point, m_eb_normal, v3, v6, vIP, distIP)) {
929 
930  amrex::Real length = (v3-v6).vectorLength();
931 
932  if ( Math::abs(distIP) < tol) {
933 #ifdef AMREX_USE_GPU
934  m_F3.add_vertex(v3);
935  m_F4.add_vertex(v3);
936  m_F5.add_vertex(v3);
937  m_F7.add_vertex(v3);
938 #else
939  m_F3.add_vertex(v3,"Path 3: v3--v6: add v3 -> F3",print_F3);
940  m_F4.add_vertex(v3,"Path 3: v3--v6: add v3 -> F4",print_F4);
941  m_F5.add_vertex(v3,"Path 3: v3--v6: add v3 -> F5",print_F5);
942  m_F7.add_vertex(v3,"Path 3: v3--v6: add v3 -> F7",print_F7);
943 #endif
944  if (!vertex_intersected[3]) {
945  vertex_intersected[3] = true;
946  ++cuts;
947  }
948  } else if ( Math::abs(distIP - length) < tol) {
949 #ifdef AMREX_USE_GPU
950  m_F3.add_vertex(v6);
951  m_F4.add_vertex(v6);
952  m_F6.add_vertex(v6);
953  m_F7.add_vertex(v6);
954 #else
955  m_F3.add_vertex(v6,"Path 3: v3--v6: add v6 -> F3",print_F3);
956  m_F4.add_vertex(v6,"Path 3: v3--v6: add v6 -> F4",print_F4);
957  m_F6.add_vertex(v6,"Path 3: v3--v6: add v6 -> F6",print_F6);
958  m_F7.add_vertex(v6,"Path 3: v3--v6: add v6 -> F7",print_F7);
959 #endif
960  if (!vertex_intersected[6]) {
961  vertex_intersected[6] = true;
962  ++cuts;
963  }
964  } else {
965 #ifdef AMREX_USE_GPU
966  m_F3.add_vertex(vIP);
967  m_F4.add_vertex(vIP);
968  m_F7.add_vertex(vIP);
969 #else
970  m_F3.add_vertex(vIP,"Path 3: v3--v6: add vIP -> F3",print_F3);
971  m_F4.add_vertex(vIP,"Path 3: v3--v6: add vIP -> F4",print_F4);
972  m_F7.add_vertex(vIP,"Path 3: v3--v6: add vIP -> F7",print_F7);
973 #endif
974  ++cuts;
975  }
976 
977  } // P3: v3--v6
978 
979  if (cuts%2 == 0) {
980 #ifdef AMREX_USE_GPU
981  m_F3.add_vertex(v6);
982  m_F4.add_vertex(v6);
983  m_F6.add_vertex(v6);
984 #else
985  m_F3.add_vertex(v6,"Path 3: after v3--v6, cuts%2 == 0: add v6 -> F3",print_F3);
986  m_F4.add_vertex(v6,"Path 3: after v3--v6, cuts%2 == 0: add v6 -> F4",print_F4);
987  m_F6.add_vertex(v6,"Path 3: after v3--v6, cuts%2 == 0: add v6 -> F6",print_F6);
988 #endif
989  }
990 
991  if (utils::intersect_plane_edge(m_eb_point, m_eb_normal, v6, v7, vIP, distIP)) {
992 
993  amrex::Real length = (v6-v7).vectorLength();
994 
995  if ( Math::abs(distIP) < tol) {
996 #ifdef AMREX_USE_GPU
997  m_F3.add_vertex(v6);
998  m_F4.add_vertex(v6);
999  m_F6.add_vertex(v6);
1000  m_F7.add_vertex(v6);
1001 #else
1002  m_F3.add_vertex(v6,"Path 3: v6--v7: add v6 -> F3",print_F3);
1003  m_F4.add_vertex(v6,"Path 3: v6--v7: add v6 -> F4",print_F4);
1004  m_F6.add_vertex(v6,"Path 3: v6--v7: add v6 -> F6",print_F6);
1005  m_F7.add_vertex(v6,"Path 3: v6--v7: add v6 -> F7",print_F7);
1006 #endif
1007  if (!vertex_intersected[6]) {
1008  vertex_intersected[6] = true;
1009  ++cuts;
1010  }
1011  } else if ( Math::abs(distIP - length) < tol) {
1012 #ifdef AMREX_USE_GPU
1013  m_F2.add_vertex(v7);
1014  m_F3.add_vertex(v7);
1015  m_F6.add_vertex(v7);
1016  m_F7.add_vertex(v7);
1017 #else
1018  m_F2.add_vertex(v7,"Path 3: v6--v7: add v7 -> F2",print_F2);
1019  m_F3.add_vertex(v7,"Path 3: v6--v7: add v7 -> F3",print_F3);
1020  m_F6.add_vertex(v7,"Path 3: v6--v7: add v7 -> F6",print_F6);
1021  m_F7.add_vertex(v7,"Path 3: v6--v7: add v7 -> F7",print_F7);
1022 #endif
1023  if (!vertex_intersected[7]) {
1024  vertex_intersected[7] = true;
1025  ++cuts;
1026  }
1027  } else {
1028 #ifdef AMREX_USE_GPU
1029  m_F3.add_vertex(vIP);
1030  m_F6.add_vertex(vIP);
1031  m_F7.add_vertex(vIP);
1032 #else
1033  m_F3.add_vertex(vIP,"Path 3: v6--v7: add vIP -> F3",print_F3);
1034  m_F6.add_vertex(vIP,"Path 3: v6--v7: add vIP -> F6",print_F6);
1035  m_F7.add_vertex(vIP,"Path 3: v6--v7: add vIP -> F7",print_F7);
1036 #endif
1037  ++cuts;
1038  }
1039  } // P3: v6--v7
1040 
1041  if (cuts == 2 && add_v7) {
1042  if (!utils::intersect_plane_edge(m_eb_point, m_eb_normal, v0, v7, vIP, distIP)) {
1043 #ifdef AMREX_USE_GPU
1044  m_F2.add_vertex(v7);
1045  m_F3.add_vertex(v7);
1046  m_F6.add_vertex(v7);
1047 #else
1048  m_F2.add_vertex(v7,"Path 3: after v6--v7, cuts == 2 && add_v7: add v7 -> F2",print_F2);
1049  m_F3.add_vertex(v7,"Path 3: after v6--v7, cuts == 2 && add_v7: add v7 -> F3",print_F3);
1050  m_F6.add_vertex(v7,"Path 3: after v6--v7, cuts == 2 && add_v7: add v7 -> F6",print_F6);
1051 #endif
1052  add_v7 = 0;
1053  }
1054  }
1055 
1056  } // end Path 3
1057 
1058  // Path 6
1059  if (utils::intersect_plane_edge(m_eb_point, m_eb_normal, v3, v4, vIP, distIP)) {
1060 
1061  amrex::Real length = (v3-v4).vectorLength();
1062 
1063  if ( Math::abs(distIP) < tol) {
1064 #ifdef AMREX_USE_GPU
1065  m_F3.add_vertex(v3);
1066  m_F4.add_vertex(v3);
1067  m_F5.add_vertex(v3);
1068  m_F7.add_vertex(v3);
1069 #else
1070  m_F3.add_vertex(v3,"Path 6: v3--v4: add v3 -> F3",print_F3);
1071  m_F4.add_vertex(v3,"Path 6: v3--v4: add v3 -> F4",print_F4);
1072  m_F5.add_vertex(v3,"Path 6: v3--v4: add v3 -> F5",print_F5);
1073  m_F7.add_vertex(v3,"Path 6: v3--v4: add v3 -> F7",print_F7);
1074 #endif
1075  } else if ( Math::abs(distIP - length) < tol) {
1076 #ifdef AMREX_USE_GPU
1077  m_F2.add_vertex(v4);
1078  m_F3.add_vertex(v4);
1079  m_F5.add_vertex(v4);
1080  m_F7.add_vertex(v4);
1081 #else
1082  m_F2.add_vertex(v4,"Path 6: v3--v4: add v4 -> F2",print_F2);
1083  m_F3.add_vertex(v4,"Path 6: v3--v4: add v4 -> F3",print_F3);
1084  m_F5.add_vertex(v4,"Path 6: v3--v4: add v4 -> F5",print_F5);
1085  m_F7.add_vertex(v4,"Path 6: v3--v4: add v4 -> F7",print_F7);
1086 #endif
1087  } else {
1088 #ifdef AMREX_USE_GPU
1089  m_F3.add_vertex(vIP);
1090  m_F5.add_vertex(vIP);
1091  m_F7.add_vertex(vIP);
1092 #else
1093  m_F3.add_vertex(vIP,"Path 6: v3--v4: add vIP -> F3",print_F3);
1094  m_F5.add_vertex(vIP,"Path 6: v3--v4: add vIP -> F5",print_F5);
1095  m_F7.add_vertex(vIP,"Path 6: v3--v4: add vIP -> F7",print_F7);
1096 #endif
1097  }
1098  } // end Path 6
1099 
1100 #ifndef AMREX_USE_GPU
1101  if (print_F1 || print_F2 || print_F3 || print_F4 || print_F5 || print_F6 || print_F7){
1102  amrex::Print()<<"\n";
1103  }
1104 #endif
1105 
1106 }
1107 
1108 
1109 
1110 AMREX_GPU_HOST_DEVICE
1111 AMREX_FORCE_INLINE
1112 void
1113 eb_cut_cell_::
1114 set_covered_regular_cell_vertices ()
1115 {
1116  using namespace amrex;
1117 
1118  RealVect v0(m_rbox.lo(0), m_rbox.lo(1), m_rbox.lo(2));
1119  RealVect v1(m_rbox.hi(0), m_rbox.lo(1), m_rbox.lo(2));
1120  RealVect v2(m_rbox.lo(0), m_rbox.hi(1), m_rbox.lo(2));
1121  RealVect v3(m_rbox.lo(0), m_rbox.lo(1), m_rbox.hi(2));
1122  RealVect v4(m_rbox.hi(0), m_rbox.lo(1), m_rbox.hi(2));
1123  RealVect v5(m_rbox.hi(0), m_rbox.hi(1), m_rbox.lo(2));
1124  RealVect v6(m_rbox.lo(0), m_rbox.hi(1), m_rbox.hi(2));
1125  RealVect v7(m_rbox.hi(0), m_rbox.hi(1), m_rbox.hi(2));
1126 
1127  // Add vertices in the order of outward normal vector
1128 
1129  m_F1.add_vertex(v0);
1130  m_F1.add_vertex(v2);
1131  m_F1.add_vertex(v5);
1132  m_F1.add_vertex(v1);
1133 
1134  m_F2.add_vertex(v1);
1135  m_F2.add_vertex(v5);
1136  m_F2.add_vertex(v7);
1137  m_F2.add_vertex(v4);
1138 
1139  m_F3.add_vertex(v3);
1140  m_F3.add_vertex(v4);
1141  m_F3.add_vertex(v7);
1142  m_F3.add_vertex(v6);
1143 
1144  m_F4.add_vertex(v0);
1145  m_F4.add_vertex(v3);
1146  m_F4.add_vertex(v6);
1147  m_F4.add_vertex(v2);
1148 
1149  m_F5.add_vertex(v0);
1150  m_F5.add_vertex(v1);
1151  m_F5.add_vertex(v4);
1152  m_F5.add_vertex(v3);
1153 
1154  m_F6.add_vertex(v2);
1155  m_F6.add_vertex(v5);
1156  m_F6.add_vertex(v7);
1157  m_F6.add_vertex(v6);
1158 
1159 }
1160 
1161 #endif
ERF_EBPolygon.H
ERF_EBUtils.H
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
eb_cut_cell_
Definition: ERF_EBCutCell.H:14
eb_cut_cell_::calc_edge_intersections
AMREX_GPU_HOST_DEVICE void calc_edge_intersections()
Definition: ERF_EBCutCell.H:341
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:1114
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_::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:235
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:179
polygon_::ok
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int ok() const noexcept
Definition: ERF_EBPolygon.H:160
polygon_::distance
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real distance(amrex::RealVect const &a_point) const noexcept
Definition: ERF_EBPolygon.H:171
polygon_::area
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real area() const noexcept
Definition: ERF_EBPolygon.H:164
polygon_::set_area
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void set_area(amrex::Real const &a_area)
Definition: ERF_EBPolygon.H:72
polygon_::add_vertex
AMREX_GPU_HOST void add_vertex(amrex::RealVect const &a_v, const std::string &msg, bool print_msg=false)
Definition: ERF_EBPolygon.H:43
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