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 > 0.);
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) < 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 (0. <= 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(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 0.; }
93  if (m_flag.isRegular() ) { return m_rbox.volume(); }
94 
95  amrex::Real volume(0.);
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 /= 3.;
106 
107  return m_invert*volume + (1.-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 0.; }
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 + (1.-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 0.; }
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 + (1.-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,0.0) || 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 + (1.-m_invert) * ((1.+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,0.0) || 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 + (1.-m_invert) * ((1.+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{0.,0.,0.};
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{0.,0.,0.};
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] = 0.5 * ( - xm * xm * axm + xp * xp * axp + bcent[0] * bcent[0] * m_eb_normal[0] * barea ) / vol;
215  vcent[1] = 0.5 * ( - ym * ym * aym + yp * yp * ayp + bcent[1] * bcent[1] * m_eb_normal[1] * barea ) / vol;
216  vcent[2] = 0.5 * ( - 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(0.0)
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(0.), amrex::RealVect(0.), amrex::RealVect(0.),
285  amrex::RealVect(0.), amrex::RealVect(0.), amrex::RealVect(0.)})
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] = 0.25 * ( v0 + v2 + v5 + v1 ); // F1
305  m_cellface_cent[1] = 0.25 * ( v1 + v5 + v7 + v4 ); // F2
306  m_cellface_cent[2] = 0.25 * ( v3 + v4 + v7 + v6 ); // F3
307  m_cellface_cent[3] = 0.25 * ( v0 + v3 + v6 + v2 ); // F4
308  m_cellface_cent[4] = 0.25 * ( v0 + v1 + v4 + v3 ); // F5
309  m_cellface_cent[5] = 0.25 * ( v2 + v5 + v7 + v6 ); // F6
310 
311  // Cell centroid
312 
313  m_cell_cent[0] = 0.5 * ( m_rbox.lo(0) + m_rbox.hi(0) );
314  m_cell_cent[1] = 0.5 * ( m_rbox.lo(1) + m_rbox.hi(1) );
315  m_cell_cent[2] = 0.5 * ( 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 = 0.5*(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)) > 0.)
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)) > 0.) ? 0.0 : 1.0;
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 = 1.e-15;
393  amrex::RealVect v_intersect;
394  amrex::Real dist_intersect;
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  // check if intersection is at vertices
402  intersected_start = (amrex::Math::abs(dist_intersect) < tol);
403  intersected_end = (amrex::Math::abs(edge_length - dist_intersect) < tol);
404  }
405 };
406 
407 AMREX_GPU_HOST_DEVICE
408 AMREX_FORCE_INLINE
409 void
410 eb_cut_cell_::
411 calc_edge_intersections ()
412 {
413  using namespace amrex;
414 
415  amrex::RealVect v0(m_rbox.lo(0), m_rbox.lo(1), m_rbox.lo(2));
416  amrex::RealVect v1(m_rbox.hi(0), m_rbox.lo(1), m_rbox.lo(2));
417  amrex::RealVect v2(m_rbox.lo(0), m_rbox.hi(1), m_rbox.lo(2));
418  amrex::RealVect v3(m_rbox.lo(0), m_rbox.lo(1), m_rbox.hi(2));
419  amrex::RealVect v4(m_rbox.hi(0), m_rbox.lo(1), m_rbox.hi(2));
420  amrex::RealVect v5(m_rbox.hi(0), m_rbox.hi(1), m_rbox.lo(2));
421  amrex::RealVect v6(m_rbox.lo(0), m_rbox.hi(1), m_rbox.hi(2));
422  amrex::RealVect v7(m_rbox.hi(0), m_rbox.hi(1), m_rbox.hi(2));
423 
424 // #ifndef AMREX_USE_GPU
425 
426  bool print_initial = ( Math::abs(m_eb_point[0]-1.666667e-01)<1.e-4 &&
427  Math::abs(m_eb_point[1]+4.194018e-01)<1.e-4 &&
428  Math::abs(m_eb_point[2]+1.666667e-01)<1.e-4 );
429  const bool print_F1 = print_initial && false;
430  const bool print_F2 = print_initial && false;
431  const bool print_F3 = print_initial && false;
432  const bool print_F4 = print_initial && false;
433  const bool print_F5 = print_initial && false;
434  const bool print_F6 = print_initial && false;
435  const bool print_F7 = print_initial && false;
436 // #endif
437 
438  m_F1.add_vertex(v0); if(print_F1) AMREX_DEVICE_PRINTF("%s \n","Initial: add v0 -> F1");
439  m_F4.add_vertex(v0); if(print_F4) AMREX_DEVICE_PRINTF("%s \n","Initial: add v0 -> F4");
440  m_F5.add_vertex(v0); if(print_F5) AMREX_DEVICE_PRINTF("%s \n","Initial: add v0 -> F5");
441 
442  int add_v7(1);
443 
444  amrex::RealVect vIP;
445  amrex::Real distIP;
446 
447  amrex::Array<bool,8> vertex_intersected{}; // true if vertex has been intersected within a path.
448 
449  //------------------------------------------------------------
450  // STEP 1: Predefine edge-plane intersections along 6 paths.
451  // Check how the edges intersect the plane and whether the configuration is valid.
452  // If the configuration is not valid, adjust it.
453  //------------------------------------------------------------
454 
455  amrex::Array<path_data,3> p1, p2, p3;
456  path_data p4, p5, p6;
457 
458  // Path 1: v0->v1->v4->v7
459  p1[0].set(m_eb_point, m_eb_normal, v0, v1);
460  p1[1].set(m_eb_point, m_eb_normal, v1, v4);
461  p1[2].set(m_eb_point, m_eb_normal, v4, v7);
462 
463  // Path 2: v0->v2->v5->v7
464  p2[0].set(m_eb_point, m_eb_normal, v0, v2);
465  p2[1].set(m_eb_point, m_eb_normal, v2, v5);
466  p2[2].set(m_eb_point, m_eb_normal, v5, v7);
467 
468  // Path 3: v0->v3->v6->v7
469  p3[0].set(m_eb_point, m_eb_normal, v0, v3);
470  p3[1].set(m_eb_point, m_eb_normal, v3, v6);
471  p3[2].set(m_eb_point, m_eb_normal, v6, v7);
472 
473  // Path 4: v1->v5
474  p4.set(m_eb_point, m_eb_normal, v1, v5);
475 
476  // Path 5: v2->v6
477  p5.set(m_eb_point, m_eb_normal, v2, v6);
478 
479  // Path 6: v3->v4
480  p6.set(m_eb_point, m_eb_normal, v3, v4);
481 
482  //------------------------------------------------------------
483  // STEP 2: Add vertices to faces.
484  //------------------------------------------------------------
485 
486  // Path 1
487  {
488  int cuts(0);
489  for (int i = 0; i < 8; ++i) vertex_intersected[i] = false;
490 
491  if (p1[0].intersected) {
492  if (p1[0].intersected_start) {
493  m_F1.add_vertex(v0); if(print_F1) AMREX_DEVICE_PRINTF("%s \n","Path 1: v0--v1: add v0 -> F1");
494  m_F4.add_vertex(v0); if(print_F4) AMREX_DEVICE_PRINTF("%s \n","Path 1: v0--v1: add v0 -> F4");
495  m_F5.add_vertex(v0); if(print_F5) AMREX_DEVICE_PRINTF("%s \n","Path 1: v0--v1: add v0 -> F5");
496  m_F7.add_vertex(v0); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 1: v0--v1: add v0 -> F7");
497  if (!vertex_intersected[0]) {
498  vertex_intersected[0] = true;
499  ++cuts;
500  }
501  } else if (p1[0].intersected_end) {
502  m_F1.add_vertex(v1); if(print_F1) AMREX_DEVICE_PRINTF("%s \n","Path 1: v0--v1: add v1 -> F1");
503  m_F2.add_vertex(v1); if(print_F2) AMREX_DEVICE_PRINTF("%s \n","Path 1: v0--v1: add v1 -> F2");
504  m_F5.add_vertex(v1); if(print_F5) AMREX_DEVICE_PRINTF("%s \n","Path 1: v0--v1: add v1 -> F5");
505  m_F7.add_vertex(v1); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 1: v0--v1: add v1 -> F7");
506  if (!vertex_intersected[1]) {
507  vertex_intersected[1] = true;
508  ++cuts;
509  }
510  } else {
511  m_F1.add_vertex(p1[0].vIP); if(print_F1) AMREX_DEVICE_PRINTF("%s \n","Path 1: v0--v1: add vIP -> F1");
512  m_F5.add_vertex(p1[0].vIP); if(print_F5) AMREX_DEVICE_PRINTF("%s \n","Path 1: v0--v1: add vIP -> F5");
513  m_F7.add_vertex(p1[0].vIP); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 1: v0--v1: add vIP -> F7");
514  ++cuts;
515  }
516 
517  } // Path 1: v0--v1
518 
519  if (cuts%2 == 0) {
520  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");
521  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");
522  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");
523  }
524 
525  if (p1[1].intersected) {
526  if (p1[1].intersected_start) {
527  m_F1.add_vertex(v1); if(print_F1) AMREX_DEVICE_PRINTF("%s \n","Path 1: v1--v4: add v1 -> F1");
528  m_F2.add_vertex(v1); if(print_F2) AMREX_DEVICE_PRINTF("%s \n","Path 1: v1--v4: add v1 -> F2");
529  m_F5.add_vertex(v1); if(print_F5) AMREX_DEVICE_PRINTF("%s \n","Path 1: v1--v4: add v1 -> F5");
530  m_F7.add_vertex(v1); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 1: v1--v4: add v1 -> F7");
531  if (!vertex_intersected[1]) {
532  vertex_intersected[1] = true;
533  ++cuts;
534  }
535  } else if (p1[1].intersected_end) {
536  m_F2.add_vertex(v4); if(print_F2) AMREX_DEVICE_PRINTF("%s \n","Path 1: v1--v4: add v4 -> F2");
537  m_F3.add_vertex(v4); if(print_F3) AMREX_DEVICE_PRINTF("%s \n","Path 1: v1--v4: add v4 -> F3");
538  m_F5.add_vertex(v4); if(print_F5) AMREX_DEVICE_PRINTF("%s \n","Path 1: v1--v4: add v4 -> F5");
539  m_F7.add_vertex(v4); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 1: v1--v4: add v4 -> F7");
540  if (!vertex_intersected[4]) {
541  vertex_intersected[4] = true;
542  ++cuts;
543  }
544  } else {
545  m_F2.add_vertex(p1[1].vIP); if(print_F2) AMREX_DEVICE_PRINTF("%s \n","Path 1: v1--v4: add vIP -> F2");
546  m_F5.add_vertex(p1[1].vIP); if(print_F5) AMREX_DEVICE_PRINTF("%s \n","Path 1: v1--v4: add vIP -> F5");
547  m_F7.add_vertex(p1[1].vIP); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 1: v1--v4: add vIP -> F7");
548  ++cuts;
549  }
550 
551  } // P1: v1--v4
552 
553  if (cuts%2 == 0) {
554  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");
555  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");
556  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");
557  }
558 
559  if (p1[2].intersected) {
560  if (p1[2].intersected_start) {
561  m_F2.add_vertex(v4); if(print_F2) AMREX_DEVICE_PRINTF("%s \n","Path 1: v4--v7: add v4 -> F2");
562  m_F3.add_vertex(v4); if(print_F3) AMREX_DEVICE_PRINTF("%s \n","Path 1: v4--v7: add v4 -> F3");
563  m_F5.add_vertex(v4); if(print_F5) AMREX_DEVICE_PRINTF("%s \n","Path 1: v4--v7: add v4 -> F5");
564  m_F7.add_vertex(v4); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 1: v4--v7: add v4 -> F7");
565  if (!vertex_intersected[4]) {
566  vertex_intersected[4] = true;
567  ++cuts;
568  }
569  } else if (p1[2].intersected_end) {
570  m_F2.add_vertex(v7); if(print_F2) AMREX_DEVICE_PRINTF("%s \n","Path 1: v4--v7: add v7 -> F2");
571  m_F3.add_vertex(v7); if(print_F3) AMREX_DEVICE_PRINTF("%s \n","Path 1: v4--v7: add v7 -> F3");
572  m_F6.add_vertex(v7); if(print_F6) AMREX_DEVICE_PRINTF("%s \n","Path 1: v4--v7: add v7 -> F6");
573  m_F7.add_vertex(v7); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 1: v4--v7: add v7 -> F7");
574  if (!vertex_intersected[7]) {
575  vertex_intersected[7] = true;
576  ++cuts;
577  }
578  } else {
579  m_F2.add_vertex(p1[2].vIP); if(print_F2) AMREX_DEVICE_PRINTF("%s \n","Path 1: v4--v7: add vIP -> F2");
580  m_F3.add_vertex(p1[2].vIP); if(print_F3) AMREX_DEVICE_PRINTF("%s \n","Path 1: v4--v7: add vIP -> F3");
581  m_F7.add_vertex(p1[2].vIP); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 1: v4--v7: add vIP -> F7");
582  ++cuts;
583  }
584 
585  } // P1: v4--v7
586 
587  if (cuts == 2 && add_v7) {
588  if (!intersect_plane_edge(m_eb_point, m_eb_normal, v0, v7, vIP, distIP)) {
589  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");
590  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");
591  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");
592  add_v7 = 0;
593  }
594  }
595  } // end Path 1
596 
597  // Path 4
598  if (p4.intersected) {
599  if (p4.intersected_start) {
600  m_F1.add_vertex(v1); if(print_F1) AMREX_DEVICE_PRINTF("%s \n","Path 4: v1--v5: add v1 -> F1");
601  m_F2.add_vertex(v1); if(print_F2) AMREX_DEVICE_PRINTF("%s \n","Path 4: v1--v5: add v1 -> F2");
602  m_F5.add_vertex(v1); if(print_F5) AMREX_DEVICE_PRINTF("%s \n","Path 4: v1--v5: add v1 -> F5");
603  m_F7.add_vertex(v1); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 4: v1--v5: add v1 -> F7");
604  } else if (p4.intersected_end) {
605  m_F1.add_vertex(v5); if(print_F1) AMREX_DEVICE_PRINTF("%s \n","Path 4: v1--v5: add v5 -> F1");
606  m_F2.add_vertex(v5); if(print_F2) AMREX_DEVICE_PRINTF("%s \n","Path 4: v1--v5: add v5 -> F2");
607  m_F6.add_vertex(v5); if(print_F6) AMREX_DEVICE_PRINTF("%s \n","Path 4: v1--v5: add v5 -> F6");
608  m_F7.add_vertex(v5); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 4: v1--v5: add v5 -> F7");
609  } else {
610  m_F1.add_vertex(p4.vIP); if(print_F1) AMREX_DEVICE_PRINTF("%s \n","Path 4: v1--v5: add vIP -> F1");
611  m_F2.add_vertex(p4.vIP); if(print_F2) AMREX_DEVICE_PRINTF("%s \n","Path 4: v1--v5: add vIP -> F2");
612  m_F7.add_vertex(p4.vIP); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 4: v1--v5: add vIP -> F7");
613  }
614  }
615 
616  // Path 2
617  { int cuts(0);
618  for (int i = 0; i < 8; ++i) vertex_intersected[i] = false;
619 
620  if (p2[0].intersected) {
621  if (p2[0].intersected_start) {
622  m_F1.add_vertex(v0); if(print_F1) AMREX_DEVICE_PRINTF("%s \n","Path 2: v0--v2: add v0 -> F1");
623  m_F4.add_vertex(v0); if(print_F4) AMREX_DEVICE_PRINTF("%s \n","Path 2: v0--v2: add v0 -> F4");
624  m_F5.add_vertex(v0); if(print_F5) AMREX_DEVICE_PRINTF("%s \n","Path 2: v0--v2: add v0 -> F5");
625  m_F7.add_vertex(v0); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 2: v0--v2: add v0 -> F7");
626  if (!vertex_intersected[0]) {
627  vertex_intersected[0] = true;
628  ++cuts;
629  }
630  } else if (p2[0].intersected_end) {
631  m_F1.add_vertex(v2); if(print_F1) AMREX_DEVICE_PRINTF("%s \n","Path 2: v0--v2: add v2 -> F1");
632  m_F4.add_vertex(v2); if(print_F4) AMREX_DEVICE_PRINTF("%s \n","Path 2: v0--v2: add v2 -> F4");
633  m_F6.add_vertex(v2); if(print_F6) AMREX_DEVICE_PRINTF("%s \n","Path 2: v0--v2: add v2 -> F6");
634  m_F7.add_vertex(v2); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 2: v0--v2: add v2 -> F7");
635  if (!vertex_intersected[2]) {
636  vertex_intersected[2] = true;
637  ++cuts;
638  }
639  } else {
640  m_F1.add_vertex(p2[0].vIP); if(print_F1) AMREX_DEVICE_PRINTF("%s \n","Path 2: v0--v2: add vIP -> F1");
641  m_F4.add_vertex(p2[0].vIP); if(print_F4) AMREX_DEVICE_PRINTF("%s \n","Path 2: v0--v2: add vIP -> F4");
642  m_F7.add_vertex(p2[0].vIP); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 2: v0--v2: add vIP -> F7");
643  ++cuts;
644  }
645 
646  } // P2: v0--v2
647 
648  if (cuts%2 == 0) {
649  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");
650  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");
651  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");
652  }
653 
654  if (p2[1].intersected) {
655  if (p2[1].intersected_start) {
656  m_F1.add_vertex(v2); if(print_F1) AMREX_DEVICE_PRINTF("%s \n","Path 2: v2--v5: add v2 -> F1");
657  m_F4.add_vertex(v2); if(print_F4) AMREX_DEVICE_PRINTF("%s \n","Path 2: v2--v5: add v2 -> F4");
658  m_F6.add_vertex(v2); if(print_F6) AMREX_DEVICE_PRINTF("%s \n","Path 2: v2--v5: add v2 -> F6");
659  m_F7.add_vertex(v2); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 2: v2--v5: add v2 -> F7");
660  if (!vertex_intersected[2]) {
661  vertex_intersected[2] = true;
662  ++cuts;
663  }
664  } else if (p2[1].intersected_end) {
665  m_F1.add_vertex(v5); if(print_F1) AMREX_DEVICE_PRINTF("%s \n","Path 2: v2--v5: add v5 -> F1");
666  m_F2.add_vertex(v5); if(print_F2) AMREX_DEVICE_PRINTF("%s \n","Path 2: v2--v5: add v5 -> F2");
667  m_F6.add_vertex(v5); if(print_F6) AMREX_DEVICE_PRINTF("%s \n","Path 2: v2--v5: add v5 -> F6");
668  m_F7.add_vertex(v5); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 2: v2--v5: add v5 -> F7");
669  if (!vertex_intersected[5]) {
670  vertex_intersected[5] = true;
671  ++cuts;
672  }
673  } else {
674  m_F1.add_vertex(p2[1].vIP); if(print_F1) AMREX_DEVICE_PRINTF("%s \n","Path 2: v2--v5: add vIP -> F1");
675  m_F6.add_vertex(p2[1].vIP); if(print_F6) AMREX_DEVICE_PRINTF("%s \n","Path 2: v2--v5: add vIP -> F6");
676  m_F7.add_vertex(p2[1].vIP); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 2: v2--v5: add vIP -> F7");
677  ++cuts;
678  }
679  }
680 
681  if (cuts%2 == 0) {
682  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");
683  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");
684  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");
685  }
686 
687  if (p2[2].intersected) {
688  if (p2[2].intersected_start) {
689  m_F1.add_vertex(v5); if(print_F1) AMREX_DEVICE_PRINTF("%s \n","Path 2: v5--v7: add v5 -> F1");
690  m_F2.add_vertex(v5); if(print_F2) AMREX_DEVICE_PRINTF("%s \n","Path 2: v5--v7: add v5 -> F2");
691  m_F6.add_vertex(v5); if(print_F6) AMREX_DEVICE_PRINTF("%s \n","Path 2: v5--v7: add v5 -> F6");
692  m_F7.add_vertex(v5); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 2: v5--v7: add v5 -> F7");
693  if (!vertex_intersected[5]) {
694  vertex_intersected[5] = true;
695  ++cuts;
696  }
697  } else if (p2[2].intersected_end) {
698  m_F2.add_vertex(v7); if(print_F2) AMREX_DEVICE_PRINTF("%s \n","Path 2: v5--v7: add v7 -> F2");
699  m_F3.add_vertex(v7); if(print_F3) AMREX_DEVICE_PRINTF("%s \n","Path 2: v5--v7: add v7 -> F3");
700  m_F6.add_vertex(v7); if(print_F6) AMREX_DEVICE_PRINTF("%s \n","Path 2: v5--v7: add v7 -> F6");
701  m_F7.add_vertex(v7); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 2: v5--v7: add v7 -> F7");
702  if (!vertex_intersected[7]) {
703  vertex_intersected[7] = true;
704  ++cuts;
705  }
706  } else {
707  m_F2.add_vertex(p2[2].vIP); if(print_F2) AMREX_DEVICE_PRINTF("%s \n","Path 2: v5--v7: add vIP -> F2");
708  m_F6.add_vertex(p2[2].vIP); if(print_F6) AMREX_DEVICE_PRINTF("%s \n","Path 2: v5--v7: add vIP -> F6");
709  m_F7.add_vertex(p2[2].vIP); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 2: v5--v7: add vIP -> F7");
710  ++cuts;
711  }
712  } // Path 2: v5--v7
713 
714  if (cuts == 2 && add_v7) {
715  if (!intersect_plane_edge(m_eb_point, m_eb_normal, v0, v7, vIP, distIP)) {
716  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");
717  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");
718  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");
719  add_v7 = 0;
720  }
721  }
722 
723  } // end Path 2
724 
725  // Path 5
726  if (p5.intersected) {
727  if (p5.intersected_start) {
728  m_F1.add_vertex(v2); if(print_F1) AMREX_DEVICE_PRINTF("%s \n","Path 5: v2--v6: add v2 -> F1");
729  m_F4.add_vertex(v2); if(print_F4) AMREX_DEVICE_PRINTF("%s \n","Path 5: v2--v6: add v2 -> F4");
730  m_F6.add_vertex(v2); if(print_F6) AMREX_DEVICE_PRINTF("%s \n","Path 5: v2--v6: add v2 -> F6");
731  m_F7.add_vertex(v2); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 5: v2--v6: add v2 -> F7");
732  } else if (p5.intersected_end) {
733  m_F3.add_vertex(v6); if(print_F3) AMREX_DEVICE_PRINTF("%s \n","Path 5: v2--v6: add v6 -> F3");
734  m_F4.add_vertex(v6); if(print_F4) AMREX_DEVICE_PRINTF("%s \n","Path 5: v2--v6: add v6 -> F4");
735  m_F6.add_vertex(v6); if(print_F6) AMREX_DEVICE_PRINTF("%s \n","Path 5: v2--v6: add v6 -> F6");
736  m_F7.add_vertex(v6); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 5: v2--v6: add v6 -> F7");
737  } else {
738  m_F4.add_vertex(p5.vIP); if(print_F4) AMREX_DEVICE_PRINTF("%s \n","Path 5: v2--v6: add vIP -> F4");
739  m_F6.add_vertex(p5.vIP); if(print_F6) AMREX_DEVICE_PRINTF("%s \n","Path 5: v2--v6: add vIP -> F6");
740  m_F7.add_vertex(p5.vIP); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 5: v2--v6: add vIP -> F7");
741  }
742  } // end Path 5
743 
744 
745  // Path 3
746  { int cuts(0);
747  for (int i = 0; i < 8; ++i) vertex_intersected[i] = false;
748 
749  if (p3[0].intersected) {
750  if (p3[0].intersected_start) {
751  m_F1.add_vertex(v0); if(print_F1) AMREX_DEVICE_PRINTF("%s \n","Path 3: v0--v3: add v0 -> F1");
752  m_F4.add_vertex(v0); if(print_F4) AMREX_DEVICE_PRINTF("%s \n","Path 3: v0--v3: add v0 -> F4");
753  m_F5.add_vertex(v0); if(print_F5) AMREX_DEVICE_PRINTF("%s \n","Path 3: v0--v3: add v0 -> F5");
754  m_F7.add_vertex(v0); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 3: v0--v3: add v0 -> F7");
755  if (!vertex_intersected[0]) {
756  vertex_intersected[0] = true;
757  ++cuts;
758  }
759  } else if (p3[0].intersected_end) {
760  m_F3.add_vertex(v3); if(print_F3) AMREX_DEVICE_PRINTF("%s \n","Path 3: v0--v3: add v3 -> F3");
761  m_F4.add_vertex(v3); if(print_F4) AMREX_DEVICE_PRINTF("%s \n","Path 3: v0--v3: add v3 -> F4");
762  m_F5.add_vertex(v3); if(print_F5) AMREX_DEVICE_PRINTF("%s \n","Path 3: v0--v3: add v3 -> F5");
763  m_F7.add_vertex(v3); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 3: v0--v3: add v3 -> F7");
764  if (!vertex_intersected[3]) {
765  vertex_intersected[3] = true;
766  ++cuts;
767  }
768  } else {
769  m_F4.add_vertex(p3[0].vIP); if(print_F4) AMREX_DEVICE_PRINTF("%s \n","Path 3: v0--v3: add vIP -> F4");
770  m_F5.add_vertex(p3[0].vIP); if(print_F5) AMREX_DEVICE_PRINTF("%s \n","Path 3: v0--v3: add vIP -> F5");
771  m_F7.add_vertex(p3[0].vIP); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 3: v0--v3: add vIP -> F7");
772  ++cuts;
773  }
774 
775  } // P3: v0--v3
776 
777  if (cuts%2 == 0) {
778  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");
779  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");
780  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");
781  }
782 
783  if (p3[1].intersected) {
784  if (p3[1].intersected_start) {
785  m_F3.add_vertex(v3); if(print_F3) AMREX_DEVICE_PRINTF("%s \n","Path 3: v3--v6: add v3 -> F3");
786  m_F4.add_vertex(v3); if(print_F4) AMREX_DEVICE_PRINTF("%s \n","Path 3: v3--v6: add v3 -> F4");
787  m_F5.add_vertex(v3); if(print_F5) AMREX_DEVICE_PRINTF("%s \n","Path 3: v3--v6: add v3 -> F5");
788  m_F7.add_vertex(v3); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 3: v3--v6: add v3 -> F7");
789  if (!vertex_intersected[3]) {
790  vertex_intersected[3] = true;
791  ++cuts;
792  }
793  } else if (p3[1].intersected_end) {
794  m_F3.add_vertex(v6); if(print_F3) AMREX_DEVICE_PRINTF("%s \n","Path 3: v3--v6: add v6 -> F3");
795  m_F4.add_vertex(v6); if(print_F4) AMREX_DEVICE_PRINTF("%s \n","Path 3: v3--v6: add v6 -> F4");
796  m_F6.add_vertex(v6); if(print_F6) AMREX_DEVICE_PRINTF("%s \n","Path 3: v3--v6: add v6 -> F6");
797  m_F7.add_vertex(v6); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 3: v3--v6: add v6 -> F7");
798  if (!vertex_intersected[6]) {
799  vertex_intersected[6] = true;
800  ++cuts;
801  }
802  } else {
803  m_F3.add_vertex(p3[1].vIP); if(print_F3) AMREX_DEVICE_PRINTF("%s \n","Path 3: v3--v6: add vIP -> F3");
804  m_F4.add_vertex(p3[1].vIP); if(print_F4) AMREX_DEVICE_PRINTF("%s \n","Path 3: v3--v6: add vIP -> F4");
805  m_F7.add_vertex(p3[1].vIP); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 3: v3--v6: add vIP -> F7");
806  ++cuts;
807  }
808 
809  } // P3: v3--v6
810 
811  if (cuts%2 == 0) {
812  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");
813  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");
814  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");
815  }
816 
817  if (p3[2].intersected) {
818  if (p3[2].intersected_start) {
819  m_F3.add_vertex(v6); if(print_F3) AMREX_DEVICE_PRINTF("%s \n","Path 3: v6--v7: add v6 -> F3");
820  m_F4.add_vertex(v6); if(print_F4) AMREX_DEVICE_PRINTF("%s \n","Path 3: v6--v7: add v6 -> F4");
821  m_F6.add_vertex(v6); if(print_F6) AMREX_DEVICE_PRINTF("%s \n","Path 3: v6--v7: add v6 -> F6");
822  m_F7.add_vertex(v6); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 3: v6--v7: add v6 -> F7");
823  if (!vertex_intersected[6]) {
824  vertex_intersected[6] = true;
825  ++cuts;
826  }
827  } else if (p3[2].intersected_end) {
828  m_F2.add_vertex(v7); if(print_F2) AMREX_DEVICE_PRINTF("%s \n","Path 3: v6--v7: add v7 -> F2");
829  m_F3.add_vertex(v7); if(print_F3) AMREX_DEVICE_PRINTF("%s \n","Path 3: v6--v7: add v7 -> F3");
830  m_F6.add_vertex(v7); if(print_F6) AMREX_DEVICE_PRINTF("%s \n","Path 3: v6--v7: add v7 -> F6");
831  m_F7.add_vertex(v7); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 3: v6--v7: add v7 -> F7");
832  if (!vertex_intersected[7]) {
833  vertex_intersected[7] = true;
834  ++cuts;
835  }
836  } else {
837  m_F3.add_vertex(p3[2].vIP); if(print_F3) AMREX_DEVICE_PRINTF("%s \n","Path 3: v6--v7: add vIP -> F3");
838  m_F6.add_vertex(p3[2].vIP); if(print_F6) AMREX_DEVICE_PRINTF("%s \n","Path 3: v6--v7: add vIP -> F6");
839  m_F7.add_vertex(p3[2].vIP); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 3: v6--v7: add vIP -> F7");
840  ++cuts;
841  }
842  } // P3: v6--v7
843 
844  if (cuts == 2 && add_v7) {
845  if (!intersect_plane_edge(m_eb_point, m_eb_normal, v0, v7, vIP, distIP)) {
846  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");
847  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");
848  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");
849  add_v7 = 0;
850  }
851  }
852 
853  } // end Path 3
854 
855  // Path 6
856  if (p6.intersected) {
857  if (p6.intersected_start) {
858  m_F3.add_vertex(v3); if(print_F3) AMREX_DEVICE_PRINTF("%s \n","Path 6: v3--v4: add v3 -> F3");
859  m_F4.add_vertex(v3); if(print_F4) AMREX_DEVICE_PRINTF("%s \n","Path 6: v3--v4: add v3 -> F4");
860  m_F5.add_vertex(v3); if(print_F5) AMREX_DEVICE_PRINTF("%s \n","Path 6: v3--v4: add v3 -> F5");
861  m_F7.add_vertex(v3); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 6: v3--v4: add v3 -> F7");
862  } else if (p6.intersected_end) {
863  m_F2.add_vertex(v4); if(print_F2) AMREX_DEVICE_PRINTF("%s \n","Path 6: v3--v4: add v4 -> F2");
864  m_F3.add_vertex(v4); if(print_F3) AMREX_DEVICE_PRINTF("%s \n","Path 6: v3--v4: add v4 -> F3\n");
865  m_F5.add_vertex(v4); if(print_F5) AMREX_DEVICE_PRINTF("%s \n","Path 6: v3--v4: add v4 -> F5");
866  m_F7.add_vertex(v4); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 6: v3--v4: add v4 -> F7");
867  } else {
868  m_F3.add_vertex(p6.vIP); if(print_F3) AMREX_DEVICE_PRINTF("%s \n","Path 6: v3--v4: add vIP -> F3");
869  m_F5.add_vertex(p6.vIP); if(print_F5) AMREX_DEVICE_PRINTF("%s \n","Path 6: v3--v4: add vIP -> F5");
870  m_F7.add_vertex(p6.vIP); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 6: v3--v4: add vIP -> F7");
871  }
872  } // end Path 6
873 
874  if (print_F1 || print_F2 || print_F3 || print_F4 || print_F5 || print_F6 || print_F7) {
875  AMREX_DEVICE_PRINTF("%s \n"," ");
876  }
877 
878 }
879 
880 AMREX_GPU_HOST_DEVICE
881 AMREX_FORCE_INLINE
882 void
883 eb_cut_cell_::
884 set_covered_regular_cell_vertices ()
885 {
886  using namespace amrex;
887 
888  RealVect v0(m_rbox.lo(0), m_rbox.lo(1), m_rbox.lo(2));
889  RealVect v1(m_rbox.hi(0), m_rbox.lo(1), m_rbox.lo(2));
890  RealVect v2(m_rbox.lo(0), m_rbox.hi(1), m_rbox.lo(2));
891  RealVect v3(m_rbox.lo(0), m_rbox.lo(1), m_rbox.hi(2));
892  RealVect v4(m_rbox.hi(0), m_rbox.lo(1), m_rbox.hi(2));
893  RealVect v5(m_rbox.hi(0), m_rbox.hi(1), m_rbox.lo(2));
894  RealVect v6(m_rbox.lo(0), m_rbox.hi(1), m_rbox.hi(2));
895  RealVect v7(m_rbox.hi(0), m_rbox.hi(1), m_rbox.hi(2));
896 
897  // Add vertices in the order of outward normal vector
898 
899  m_F1.add_vertex(v0);
900  m_F1.add_vertex(v2);
901  m_F1.add_vertex(v5);
902  m_F1.add_vertex(v1);
903 
904  m_F2.add_vertex(v1);
905  m_F2.add_vertex(v5);
906  m_F2.add_vertex(v7);
907  m_F2.add_vertex(v4);
908 
909  m_F3.add_vertex(v3);
910  m_F3.add_vertex(v4);
911  m_F3.add_vertex(v7);
912  m_F3.add_vertex(v6);
913 
914  m_F4.add_vertex(v0);
915  m_F4.add_vertex(v3);
916  m_F4.add_vertex(v6);
917  m_F4.add_vertex(v2);
918 
919  m_F5.add_vertex(v0);
920  m_F5.add_vertex(v1);
921  m_F5.add_vertex(v4);
922  m_F5.add_vertex(v3);
923 
924  m_F6.add_vertex(v2);
925  m_F6.add_vertex(v5);
926  m_F6.add_vertex(v7);
927  m_F6.add_vertex(v6);
928 
929 }
930 
931 #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:16
ERF_EBPolygon.H
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:411
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:884
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