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