ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Classes | Public Member Functions | Private Member Functions | Private Attributes | Static Private Attributes | List of all members
eb_cut_cell_ Class Reference

Reconstructs geometric moments for one EB cut cell. More...

#include <ERF_EBCutCell.H>

Collaboration diagram for eb_cut_cell_:
Collaboration graph
[legend]

Classes

struct  path_data
 Intersection metadata for one ordered edge path in the cut-cell box. More...
 

Public Member Functions

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)
 Construct cut-cell geometry from a cell flag, box, EB point, and EB normal. More...
 
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool isCovered () const noexcept
 Return true if the reconstructed cell is covered. More...
 
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool isRegular () const noexcept
 Return true if the reconstructed cell is regular. More...
 
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool isSingleValued () const noexcept
 Return true if the reconstructed cell is single-valued. More...
 
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void set_covered ()
 Mark the reconstructed cell as covered. More...
 
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void set_regular ()
 Mark the reconstructed cell as regular and fill regular face areas. More...
 
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real volume () const
 Return the physical volume of the reconstructed cell. More...
 
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real areaLo (int const idir) const noexcept
 Return the low-side face area in one coordinate direction. More...
 
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real areaHi (int const idir) const noexcept
 Return the high-side face area in one coordinate direction. More...
 
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real areaBoun () const noexcept
 Return the EB boundary area inside the cell. More...
 
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::RealVect centLo (int const idir) const noexcept
 Return the low-side face centroid in one coordinate direction. More...
 
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::RealVect centHi (int const idir) const noexcept
 Return the high-side face centroid in one coordinate direction. More...
 
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::RealVect centBoun () const noexcept
 Return the EB boundary centroid. More...
 
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::RealVect normBoun () const noexcept
 Return the EB boundary normal. More...
 
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::RealVect centVol () const noexcept
 Return the volume centroid of the reconstructed cell. More...
 
void debug (int const a_face=-1)
 Print cut-cell diagnostics on CPU builds. More...
 

Private Member Functions

AMREX_GPU_HOST_DEVICE void calc_edge_intersections ()
 Compute all edge-plane intersections and populate polygon vertices. More...
 
AMREX_GPU_HOST_DEVICE void set_covered_regular_cell_vertices ()
 Populate polygon vertices for cells that are fully covered or regular. More...
 

Private Attributes

amrex::RealBox const m_rbox
 
amrex::RealVect const m_eb_point
 
amrex::RealVect const m_eb_normal
 
amrex::Real m_invert
 
amrex::RealVect m_rbox_area
 
amrex::EBCellFlag m_flag
 
polygon_ m_F1
 
polygon_ m_F2
 
polygon_ m_F3
 
polygon_ m_F4
 
polygon_ m_F5
 
polygon_ m_F6
 
amrex::Array< amrex::RealVect, m_max_faces > m_cellface_cent
 
amrex::RealVect m_cell_cent
 
amrex::RealVect m_dx
 
polygon_ m_F7
 
amrex::Array< polygon_ const *const, 3 > m_lo_faces {&m_F4, &m_F5, &m_F1}
 
amrex::Array< polygon_ const *const, 3 > m_hi_faces {&m_F2, &m_F6, &m_F3}
 
amrex::Array< int const, 3 > m_lo_faces_id {3,4,0}
 
amrex::Array< int const, 3 > m_hi_faces_id {1,5,2}
 

Static Private Attributes

static constexpr int m_max_faces = 6
 

Detailed Description

Reconstructs geometric moments for one EB cut cell.

The class clips a rectangular cell by a planar EB approximation and exposes volume, face areas, centroids, boundary area, boundary centroid, and boundary normal data used by auxiliary face-centered EB factories.

Constructor & Destructor Documentation

◆ eb_cut_cell_()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE eb_cut_cell_::eb_cut_cell_ ( amrex::EBCellFlag const &  a_flag,
amrex::RealBox const &  a_rbox,
amrex::RealVect const &  a_point,
amrex::RealVect const &  a_normal 
)

Construct cut-cell geometry from a cell flag, box, EB point, and EB normal.

Parameters
a_flagAMReX EB cell flag for the cell.
a_rboxPhysical-space cell box.
a_pointPoint on the EB plane.
a_normalUnit normal for the EB plane.
335  : m_rbox(a_rbox)
336  , m_eb_point(a_point)
337  , m_eb_normal(a_normal)
338  , m_invert(zero)
339  , m_F1(a_point, a_normal)
340  , m_F2(a_point, a_normal)
341  , m_F3(a_point, a_normal)
342  , m_F4(a_point, a_normal)
343  , m_F5(a_point, a_normal)
344  , m_F6(a_point, a_normal)
345  , m_cellface_cent({amrex::RealVect(zero), amrex::RealVect(zero), amrex::RealVect(zero),
346  amrex::RealVect(zero), amrex::RealVect(zero), amrex::RealVect(zero)})
347 {
348  using namespace amrex;
349 
350  m_rbox_area[0] = m_rbox.length(1)*m_rbox.length(2);
351  m_rbox_area[1] = m_rbox.length(0)*m_rbox.length(2);
352  m_rbox_area[2] = m_rbox.length(0)*m_rbox.length(1);
353 
354  amrex::RealVect v0(m_rbox.lo(0), m_rbox.lo(1), m_rbox.lo(2));
355  amrex::RealVect v1(m_rbox.hi(0), m_rbox.lo(1), m_rbox.lo(2));
356  amrex::RealVect v2(m_rbox.lo(0), m_rbox.hi(1), m_rbox.lo(2));
357  amrex::RealVect v3(m_rbox.lo(0), m_rbox.lo(1), m_rbox.hi(2));
358  amrex::RealVect v4(m_rbox.hi(0), m_rbox.lo(1), m_rbox.hi(2));
359  amrex::RealVect v5(m_rbox.hi(0), m_rbox.hi(1), m_rbox.lo(2));
360  amrex::RealVect v6(m_rbox.lo(0), m_rbox.hi(1), m_rbox.hi(2));
361  amrex::RealVect v7(m_rbox.hi(0), m_rbox.hi(1), m_rbox.hi(2));
362 
363  // Centroids of cell faces
364 
365  m_cellface_cent[0] = fourth * ( v0 + v2 + v5 + v1 ); // F1
366  m_cellface_cent[1] = fourth * ( v1 + v5 + v7 + v4 ); // F2
367  m_cellface_cent[2] = fourth * ( v3 + v4 + v7 + v6 ); // F3
368  m_cellface_cent[3] = fourth * ( v0 + v3 + v6 + v2 ); // F4
369  m_cellface_cent[4] = fourth * ( v0 + v1 + v4 + v3 ); // F5
370  m_cellface_cent[5] = fourth * ( v2 + v5 + v7 + v6 ); // F6
371 
372  // Cell centroid
373 
374  m_cell_cent[0] = myhalf * ( m_rbox.lo(0) + m_rbox.hi(0) );
375  m_cell_cent[1] = myhalf * ( m_rbox.lo(1) + m_rbox.hi(1) );
376  m_cell_cent[2] = myhalf * ( m_rbox.lo(2) + m_rbox.hi(2) );
377 
378  // Cell size
379 
380  m_dx[0] = m_rbox.hi(0) - m_rbox.lo(0);
381  m_dx[1] = m_rbox.hi(1) - m_rbox.lo(1);
382  m_dx[2] = m_rbox.hi(2) - m_rbox.lo(2);
383 
384  if (a_flag.isCovered() ) {
385 
386  set_covered();
387 
388  } else if (a_flag.isRegular() ) {
389 
390  set_regular();
391 
392  } else { // Check that the box and plane intersect.
393 
394  amrex::RealVect c = myhalf*(v0 + v7);
395  amrex::RealVect e = v7 - c;
396 
397  amrex::Real r = e[0]*amrex::Math::abs(a_normal[0]) +
398  e[1]*amrex::Math::abs(a_normal[1]) +
399  e[2]*amrex::Math::abs(a_normal[2]);
400 
401  amrex::Real s = amrex::Math::abs(c.dotProduct(a_normal)
402  - a_point.dotProduct(a_normal));
403 
404  if (s > r) {
405  if ((a_normal.dotProduct(v0) - a_normal.dotProduct(a_point)) > zero)
406  { set_covered(); } else { set_regular(); }
407  } else { m_flag.setSingleValued(); }
408  }
409 
410  if ( m_flag.isSingleValued() ) {
411 
412  m_invert = ((m_eb_normal.dotProduct(v0) - m_eb_normal.dotProduct(m_eb_point)) > zero) ? zero : one;
413 
414  calc_edge_intersections();
415 
416  } // end singleValued
417 
418  m_F1.define();
419  m_F2.define();
420  m_F3.define();
421  m_F4.define();
422  m_F5.define();
423  m_F6.define();
424  m_F7.define();
425 
426  // For covered and regular cut cells, add vertices to utilize e.g., get_centroid.
427  if ( !m_flag.isSingleValued() ) {
428  set_covered_regular_cell_vertices();
429  }
430 
431 }
one
constexpr amrex::Real one
Definition: ERF_Constants.H:9
fourth
constexpr amrex::Real fourth
Definition: ERF_Constants.H:14
zero
constexpr amrex::Real zero
Definition: ERF_Constants.H:8
myhalf
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:13
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
eb_cut_cell_::calc_edge_intersections
AMREX_GPU_HOST_DEVICE void calc_edge_intersections()
Compute all edge-plane intersections and populate polygon vertices.
Definition: ERF_EBCutCell.H:495
eb_cut_cell_::m_F6
polygon_ m_F6
Definition: ERF_EBCutCell.H:300
eb_cut_cell_::m_rbox_area
amrex::RealVect m_rbox_area
Definition: ERF_EBCutCell.H:289
eb_cut_cell_::m_cell_cent
amrex::RealVect m_cell_cent
Definition: ERF_EBCutCell.H:304
eb_cut_cell_::m_dx
amrex::RealVect m_dx
Definition: ERF_EBCutCell.H:305
eb_cut_cell_::m_eb_point
amrex::RealVect const m_eb_point
Definition: ERF_EBCutCell.H:284
eb_cut_cell_::m_invert
amrex::Real m_invert
Definition: ERF_EBCutCell.H:287
eb_cut_cell_::m_rbox
amrex::RealBox const m_rbox
Definition: ERF_EBCutCell.H:283
eb_cut_cell_::set_covered_regular_cell_vertices
AMREX_GPU_HOST_DEVICE void set_covered_regular_cell_vertices()
Populate polygon vertices for cells that are fully covered or regular.
Definition: ERF_EBCutCell.H:968
eb_cut_cell_::set_regular
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void set_regular()
Mark the reconstructed cell as regular and fill regular face areas.
Definition: ERF_EBCutCell.H:105
eb_cut_cell_::m_F1
polygon_ m_F1
Definition: ERF_EBCutCell.H:295
eb_cut_cell_::m_F3
polygon_ m_F3
Definition: ERF_EBCutCell.H:297
eb_cut_cell_::m_flag
amrex::EBCellFlag m_flag
Definition: ERF_EBCutCell.H:291
eb_cut_cell_::m_F5
polygon_ m_F5
Definition: ERF_EBCutCell.H:299
eb_cut_cell_::set_covered
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void set_covered()
Mark the reconstructed cell as covered.
Definition: ERF_EBCutCell.H:99
eb_cut_cell_::m_F7
polygon_ m_F7
Definition: ERF_EBCutCell.H:308
eb_cut_cell_::m_eb_normal
amrex::RealVect const m_eb_normal
Definition: ERF_EBCutCell.H:285
eb_cut_cell_::m_F2
polygon_ m_F2
Definition: ERF_EBCutCell.H:296
eb_cut_cell_::m_F4
polygon_ m_F4
Definition: ERF_EBCutCell.H:298
eb_cut_cell_::m_cellface_cent
amrex::Array< amrex::RealVect, m_max_faces > m_cellface_cent
Definition: ERF_EBCutCell.H:303
polygon_::define
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void define()
Finalize polygon ordering and compute its area.
Definition: ERF_EBPolygon.H:96
amrex
Definition: ERF_ConsoleIO.cpp:15

Member Function Documentation

◆ areaBoun()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real eb_cut_cell_::areaBoun ( ) const
inlinenoexcept

Return the EB boundary area inside the cell.

172  {
173  return m_F7.area();
174  }
polygon_::area
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real area() const noexcept
Return the polygon area.
Definition: ERF_EBPolygon.H:187

Referenced by eb_aux_::define().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ areaHi()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real eb_cut_cell_::areaHi ( int const  idir) const
inlinenoexcept

Return the high-side face area in one coordinate direction.

Parameters
idirCoordinate direction.
162  {
163  AMREX_ASSERT( idir >= 0 && idir < AMREX_SPACEDIM );
164  if (m_flag.isCovered() ) { return zero; }
165  if (m_flag.isRegular() ) { return m_hi_faces[idir]->area(); }
166  amrex::Real const area(m_hi_faces[idir]->area());
167  return m_invert*area + (one-m_invert)*(m_rbox_area[idir] - area);
168  }
eb_cut_cell_::m_hi_faces
amrex::Array< polygon_ const *const, 3 > m_hi_faces
Definition: ERF_EBCutCell.H:311

Referenced by centVol(), and eb_aux_::define().

Here is the caller graph for this function:

◆ areaLo()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real eb_cut_cell_::areaLo ( int const  idir) const
inlinenoexcept

Return the low-side face area in one coordinate direction.

Parameters
idirCoordinate direction.
149  {
150  AMREX_ASSERT( idir >=0 && idir < AMREX_SPACEDIM );
151  if (m_flag.isCovered() ) { return zero; }
152  if (m_flag.isRegular() ) { return m_lo_faces[idir]->area(); }
153  amrex::Real const area(m_lo_faces[idir]->area());
154  return m_invert*area + (one-m_invert)*(m_rbox_area[idir] - area);
155  }
eb_cut_cell_::m_lo_faces
amrex::Array< polygon_ const *const, 3 > m_lo_faces
Definition: ERF_EBCutCell.H:310

Referenced by centVol(), and eb_aux_::define().

Here is the caller graph for this function:

◆ calc_edge_intersections()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void eb_cut_cell_::calc_edge_intersections ( )
private

Compute all edge-plane intersections and populate polygon vertices.

496 {
497  using namespace amrex;
498 
499  amrex::RealVect v0(m_rbox.lo(0), m_rbox.lo(1), m_rbox.lo(2));
500  amrex::RealVect v1(m_rbox.hi(0), m_rbox.lo(1), m_rbox.lo(2));
501  amrex::RealVect v2(m_rbox.lo(0), m_rbox.hi(1), m_rbox.lo(2));
502  amrex::RealVect v3(m_rbox.lo(0), m_rbox.lo(1), m_rbox.hi(2));
503  amrex::RealVect v4(m_rbox.hi(0), m_rbox.lo(1), m_rbox.hi(2));
504  amrex::RealVect v5(m_rbox.hi(0), m_rbox.hi(1), m_rbox.lo(2));
505  amrex::RealVect v6(m_rbox.lo(0), m_rbox.hi(1), m_rbox.hi(2));
506  amrex::RealVect v7(m_rbox.hi(0), m_rbox.hi(1), m_rbox.hi(2));
507 
508 // #ifndef AMREX_USE_GPU
509 
510  bool print_initial = ( Math::abs(m_eb_point[0]-Real(1.666667e-01))<Real(1.e-4) &&
511  Math::abs(m_eb_point[1]+Real(4.194018e-01))<Real(1.e-4) &&
512  Math::abs(m_eb_point[2]+Real(1.666667e-01))<Real(1.e-4) );
513  const bool print_F1 = print_initial && false;
514  const bool print_F2 = print_initial && false;
515  const bool print_F3 = print_initial && false;
516  const bool print_F4 = print_initial && false;
517  const bool print_F5 = print_initial && false;
518  const bool print_F6 = print_initial && false;
519  const bool print_F7 = print_initial && false;
520 // #endif
521 
522  m_F1.add_vertex(v0); if(print_F1) AMREX_DEVICE_PRINTF("%s \n","Initial: add v0 -> F1");
523  m_F4.add_vertex(v0); if(print_F4) AMREX_DEVICE_PRINTF("%s \n","Initial: add v0 -> F4");
524  m_F5.add_vertex(v0); if(print_F5) AMREX_DEVICE_PRINTF("%s \n","Initial: add v0 -> F5");
525 
526  int add_v7(1);
527 
528  amrex::RealVect vIP;
529  amrex::Real distIP;
530 
531  amrex::Array<bool,8> vertex_intersected{}; // true if vertex has been intersected within a path.
532 
533  //------------------------------------------------------------
534  // STEP 1: Predefine edge-plane intersections along 6 paths.
535  // Check how the edges intersect the plane and whether the configuration is valid.
536  // If the configuration is not valid, adjust it.
537  //------------------------------------------------------------
538 
539  amrex::Array<path_data,3> p1, p2, p3;
540  path_data p4, p5, p6;
541 
542  // Path 1: v0->v1->v4->v7
543  p1[0].set(m_eb_point, m_eb_normal, v0, v1);
544  p1[1].set(m_eb_point, m_eb_normal, v1, v4);
545  p1[2].set(m_eb_point, m_eb_normal, v4, v7);
546 
547  // Path 2: v0->v2->v5->v7
548  p2[0].set(m_eb_point, m_eb_normal, v0, v2);
549  p2[1].set(m_eb_point, m_eb_normal, v2, v5);
550  p2[2].set(m_eb_point, m_eb_normal, v5, v7);
551 
552  // Path 3: v0->v3->v6->v7
553  p3[0].set(m_eb_point, m_eb_normal, v0, v3);
554  p3[1].set(m_eb_point, m_eb_normal, v3, v6);
555  p3[2].set(m_eb_point, m_eb_normal, v6, v7);
556 
557  // Path 4: v1->v5
558  p4.set(m_eb_point, m_eb_normal, v1, v5);
559 
560  // Path 5: v2->v6
561  p5.set(m_eb_point, m_eb_normal, v2, v6);
562 
563  // Path 6: v3->v4
564  p6.set(m_eb_point, m_eb_normal, v3, v4);
565 
566  //------------------------------------------------------------
567  // STEP 2: Add vertices to faces.
568  //------------------------------------------------------------
569 
570  // Path 1
571  {
572  int cuts(0);
573  for (int i = 0; i < 8; ++i) vertex_intersected[i] = false;
574 
575  if (p1[0].intersected) {
576  if (p1[0].intersected_start) {
577  m_F1.add_vertex(v0); if(print_F1) AMREX_DEVICE_PRINTF("%s \n","Path 1: v0--v1: add v0 -> F1");
578  m_F4.add_vertex(v0); if(print_F4) AMREX_DEVICE_PRINTF("%s \n","Path 1: v0--v1: add v0 -> F4");
579  m_F5.add_vertex(v0); if(print_F5) AMREX_DEVICE_PRINTF("%s \n","Path 1: v0--v1: add v0 -> F5");
580  m_F7.add_vertex(v0); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 1: v0--v1: add v0 -> F7");
581  if (!vertex_intersected[0]) {
582  vertex_intersected[0] = true;
583  ++cuts;
584  }
585  } else if (p1[0].intersected_end) {
586  m_F1.add_vertex(v1); if(print_F1) AMREX_DEVICE_PRINTF("%s \n","Path 1: v0--v1: add v1 -> F1");
587  m_F2.add_vertex(v1); if(print_F2) AMREX_DEVICE_PRINTF("%s \n","Path 1: v0--v1: add v1 -> F2");
588  m_F5.add_vertex(v1); if(print_F5) AMREX_DEVICE_PRINTF("%s \n","Path 1: v0--v1: add v1 -> F5");
589  m_F7.add_vertex(v1); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 1: v0--v1: add v1 -> F7");
590  if (!vertex_intersected[1]) {
591  vertex_intersected[1] = true;
592  ++cuts;
593  }
594  } else {
595  m_F1.add_vertex(p1[0].vIP); if(print_F1) AMREX_DEVICE_PRINTF("%s \n","Path 1: v0--v1: add vIP -> F1");
596  m_F5.add_vertex(p1[0].vIP); if(print_F5) AMREX_DEVICE_PRINTF("%s \n","Path 1: v0--v1: add vIP -> F5");
597  m_F7.add_vertex(p1[0].vIP); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 1: v0--v1: add vIP -> F7");
598  ++cuts;
599  }
600 
601  } // Path 1: v0--v1
602 
603  if (cuts%2 == 0) {
604  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");
605  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");
606  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");
607  }
608 
609  if (p1[1].intersected) {
610  if (p1[1].intersected_start) {
611  m_F1.add_vertex(v1); if(print_F1) AMREX_DEVICE_PRINTF("%s \n","Path 1: v1--v4: add v1 -> F1");
612  m_F2.add_vertex(v1); if(print_F2) AMREX_DEVICE_PRINTF("%s \n","Path 1: v1--v4: add v1 -> F2");
613  m_F5.add_vertex(v1); if(print_F5) AMREX_DEVICE_PRINTF("%s \n","Path 1: v1--v4: add v1 -> F5");
614  m_F7.add_vertex(v1); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 1: v1--v4: add v1 -> F7");
615  if (!vertex_intersected[1]) {
616  vertex_intersected[1] = true;
617  ++cuts;
618  }
619  } else if (p1[1].intersected_end) {
620  m_F2.add_vertex(v4); if(print_F2) AMREX_DEVICE_PRINTF("%s \n","Path 1: v1--v4: add v4 -> F2");
621  m_F3.add_vertex(v4); if(print_F3) AMREX_DEVICE_PRINTF("%s \n","Path 1: v1--v4: add v4 -> F3");
622  m_F5.add_vertex(v4); if(print_F5) AMREX_DEVICE_PRINTF("%s \n","Path 1: v1--v4: add v4 -> F5");
623  m_F7.add_vertex(v4); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 1: v1--v4: add v4 -> F7");
624  if (!vertex_intersected[4]) {
625  vertex_intersected[4] = true;
626  ++cuts;
627  }
628  } else {
629  m_F2.add_vertex(p1[1].vIP); if(print_F2) AMREX_DEVICE_PRINTF("%s \n","Path 1: v1--v4: add vIP -> F2");
630  m_F5.add_vertex(p1[1].vIP); if(print_F5) AMREX_DEVICE_PRINTF("%s \n","Path 1: v1--v4: add vIP -> F5");
631  m_F7.add_vertex(p1[1].vIP); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 1: v1--v4: add vIP -> F7");
632  ++cuts;
633  }
634 
635  } // P1: v1--v4
636 
637  if (cuts%2 == 0) {
638  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");
639  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");
640  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");
641  }
642 
643  if (p1[2].intersected) {
644  if (p1[2].intersected_start) {
645  m_F2.add_vertex(v4); if(print_F2) AMREX_DEVICE_PRINTF("%s \n","Path 1: v4--v7: add v4 -> F2");
646  m_F3.add_vertex(v4); if(print_F3) AMREX_DEVICE_PRINTF("%s \n","Path 1: v4--v7: add v4 -> F3");
647  m_F5.add_vertex(v4); if(print_F5) AMREX_DEVICE_PRINTF("%s \n","Path 1: v4--v7: add v4 -> F5");
648  m_F7.add_vertex(v4); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 1: v4--v7: add v4 -> F7");
649  if (!vertex_intersected[4]) {
650  vertex_intersected[4] = true;
651  ++cuts;
652  }
653  } else if (p1[2].intersected_end) {
654  m_F2.add_vertex(v7); if(print_F2) AMREX_DEVICE_PRINTF("%s \n","Path 1: v4--v7: add v7 -> F2");
655  m_F3.add_vertex(v7); if(print_F3) AMREX_DEVICE_PRINTF("%s \n","Path 1: v4--v7: add v7 -> F3");
656  m_F6.add_vertex(v7); if(print_F6) AMREX_DEVICE_PRINTF("%s \n","Path 1: v4--v7: add v7 -> F6");
657  m_F7.add_vertex(v7); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 1: v4--v7: add v7 -> F7");
658  if (!vertex_intersected[7]) {
659  vertex_intersected[7] = true;
660  ++cuts;
661  }
662  } else {
663  m_F2.add_vertex(p1[2].vIP); if(print_F2) AMREX_DEVICE_PRINTF("%s \n","Path 1: v4--v7: add vIP -> F2");
664  m_F3.add_vertex(p1[2].vIP); if(print_F3) AMREX_DEVICE_PRINTF("%s \n","Path 1: v4--v7: add vIP -> F3");
665  m_F7.add_vertex(p1[2].vIP); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 1: v4--v7: add vIP -> F7");
666  ++cuts;
667  }
668 
669  } // P1: v4--v7
670 
671  if (cuts == 2 && add_v7) {
672  if (!intersect_plane_edge(m_eb_point, m_eb_normal, v0, v7, vIP, distIP)) {
673  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");
674  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");
675  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");
676  add_v7 = 0;
677  }
678  }
679  } // end Path 1
680 
681  // Path 4
682  if (p4.intersected) {
683  if (p4.intersected_start) {
684  m_F1.add_vertex(v1); if(print_F1) AMREX_DEVICE_PRINTF("%s \n","Path 4: v1--v5: add v1 -> F1");
685  m_F2.add_vertex(v1); if(print_F2) AMREX_DEVICE_PRINTF("%s \n","Path 4: v1--v5: add v1 -> F2");
686  m_F5.add_vertex(v1); if(print_F5) AMREX_DEVICE_PRINTF("%s \n","Path 4: v1--v5: add v1 -> F5");
687  m_F7.add_vertex(v1); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 4: v1--v5: add v1 -> F7");
688  } else if (p4.intersected_end) {
689  m_F1.add_vertex(v5); if(print_F1) AMREX_DEVICE_PRINTF("%s \n","Path 4: v1--v5: add v5 -> F1");
690  m_F2.add_vertex(v5); if(print_F2) AMREX_DEVICE_PRINTF("%s \n","Path 4: v1--v5: add v5 -> F2");
691  m_F6.add_vertex(v5); if(print_F6) AMREX_DEVICE_PRINTF("%s \n","Path 4: v1--v5: add v5 -> F6");
692  m_F7.add_vertex(v5); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 4: v1--v5: add v5 -> F7");
693  } else {
694  m_F1.add_vertex(p4.vIP); if(print_F1) AMREX_DEVICE_PRINTF("%s \n","Path 4: v1--v5: add vIP -> F1");
695  m_F2.add_vertex(p4.vIP); if(print_F2) AMREX_DEVICE_PRINTF("%s \n","Path 4: v1--v5: add vIP -> F2");
696  m_F7.add_vertex(p4.vIP); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 4: v1--v5: add vIP -> F7");
697  }
698  }
699 
700  // Path 2
701  { int cuts(0);
702  for (int i = 0; i < 8; ++i) vertex_intersected[i] = false;
703 
704  if (p2[0].intersected) {
705  if (p2[0].intersected_start) {
706  m_F1.add_vertex(v0); if(print_F1) AMREX_DEVICE_PRINTF("%s \n","Path 2: v0--v2: add v0 -> F1");
707  m_F4.add_vertex(v0); if(print_F4) AMREX_DEVICE_PRINTF("%s \n","Path 2: v0--v2: add v0 -> F4");
708  m_F5.add_vertex(v0); if(print_F5) AMREX_DEVICE_PRINTF("%s \n","Path 2: v0--v2: add v0 -> F5");
709  m_F7.add_vertex(v0); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 2: v0--v2: add v0 -> F7");
710  if (!vertex_intersected[0]) {
711  vertex_intersected[0] = true;
712  ++cuts;
713  }
714  } else if (p2[0].intersected_end) {
715  m_F1.add_vertex(v2); if(print_F1) AMREX_DEVICE_PRINTF("%s \n","Path 2: v0--v2: add v2 -> F1");
716  m_F4.add_vertex(v2); if(print_F4) AMREX_DEVICE_PRINTF("%s \n","Path 2: v0--v2: add v2 -> F4");
717  m_F6.add_vertex(v2); if(print_F6) AMREX_DEVICE_PRINTF("%s \n","Path 2: v0--v2: add v2 -> F6");
718  m_F7.add_vertex(v2); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 2: v0--v2: add v2 -> F7");
719  if (!vertex_intersected[2]) {
720  vertex_intersected[2] = true;
721  ++cuts;
722  }
723  } else {
724  m_F1.add_vertex(p2[0].vIP); if(print_F1) AMREX_DEVICE_PRINTF("%s \n","Path 2: v0--v2: add vIP -> F1");
725  m_F4.add_vertex(p2[0].vIP); if(print_F4) AMREX_DEVICE_PRINTF("%s \n","Path 2: v0--v2: add vIP -> F4");
726  m_F7.add_vertex(p2[0].vIP); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 2: v0--v2: add vIP -> F7");
727  ++cuts;
728  }
729 
730  } // P2: v0--v2
731 
732  if (cuts%2 == 0) {
733  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");
734  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");
735  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");
736  }
737 
738  if (p2[1].intersected) {
739  if (p2[1].intersected_start) {
740  m_F1.add_vertex(v2); if(print_F1) AMREX_DEVICE_PRINTF("%s \n","Path 2: v2--v5: add v2 -> F1");
741  m_F4.add_vertex(v2); if(print_F4) AMREX_DEVICE_PRINTF("%s \n","Path 2: v2--v5: add v2 -> F4");
742  m_F6.add_vertex(v2); if(print_F6) AMREX_DEVICE_PRINTF("%s \n","Path 2: v2--v5: add v2 -> F6");
743  m_F7.add_vertex(v2); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 2: v2--v5: add v2 -> F7");
744  if (!vertex_intersected[2]) {
745  vertex_intersected[2] = true;
746  ++cuts;
747  }
748  } else if (p2[1].intersected_end) {
749  m_F1.add_vertex(v5); if(print_F1) AMREX_DEVICE_PRINTF("%s \n","Path 2: v2--v5: add v5 -> F1");
750  m_F2.add_vertex(v5); if(print_F2) AMREX_DEVICE_PRINTF("%s \n","Path 2: v2--v5: add v5 -> F2");
751  m_F6.add_vertex(v5); if(print_F6) AMREX_DEVICE_PRINTF("%s \n","Path 2: v2--v5: add v5 -> F6");
752  m_F7.add_vertex(v5); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 2: v2--v5: add v5 -> F7");
753  if (!vertex_intersected[5]) {
754  vertex_intersected[5] = true;
755  ++cuts;
756  }
757  } else {
758  m_F1.add_vertex(p2[1].vIP); if(print_F1) AMREX_DEVICE_PRINTF("%s \n","Path 2: v2--v5: add vIP -> F1");
759  m_F6.add_vertex(p2[1].vIP); if(print_F6) AMREX_DEVICE_PRINTF("%s \n","Path 2: v2--v5: add vIP -> F6");
760  m_F7.add_vertex(p2[1].vIP); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 2: v2--v5: add vIP -> F7");
761  ++cuts;
762  }
763  }
764 
765  if (cuts%2 == 0) {
766  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");
767  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");
768  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");
769  }
770 
771  if (p2[2].intersected) {
772  if (p2[2].intersected_start) {
773  m_F1.add_vertex(v5); if(print_F1) AMREX_DEVICE_PRINTF("%s \n","Path 2: v5--v7: add v5 -> F1");
774  m_F2.add_vertex(v5); if(print_F2) AMREX_DEVICE_PRINTF("%s \n","Path 2: v5--v7: add v5 -> F2");
775  m_F6.add_vertex(v5); if(print_F6) AMREX_DEVICE_PRINTF("%s \n","Path 2: v5--v7: add v5 -> F6");
776  m_F7.add_vertex(v5); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 2: v5--v7: add v5 -> F7");
777  if (!vertex_intersected[5]) {
778  vertex_intersected[5] = true;
779  ++cuts;
780  }
781  } else if (p2[2].intersected_end) {
782  m_F2.add_vertex(v7); if(print_F2) AMREX_DEVICE_PRINTF("%s \n","Path 2: v5--v7: add v7 -> F2");
783  m_F3.add_vertex(v7); if(print_F3) AMREX_DEVICE_PRINTF("%s \n","Path 2: v5--v7: add v7 -> F3");
784  m_F6.add_vertex(v7); if(print_F6) AMREX_DEVICE_PRINTF("%s \n","Path 2: v5--v7: add v7 -> F6");
785  m_F7.add_vertex(v7); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 2: v5--v7: add v7 -> F7");
786  if (!vertex_intersected[7]) {
787  vertex_intersected[7] = true;
788  ++cuts;
789  }
790  } else {
791  m_F2.add_vertex(p2[2].vIP); if(print_F2) AMREX_DEVICE_PRINTF("%s \n","Path 2: v5--v7: add vIP -> F2");
792  m_F6.add_vertex(p2[2].vIP); if(print_F6) AMREX_DEVICE_PRINTF("%s \n","Path 2: v5--v7: add vIP -> F6");
793  m_F7.add_vertex(p2[2].vIP); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 2: v5--v7: add vIP -> F7");
794  ++cuts;
795  }
796  } // Path 2: v5--v7
797 
798  if (cuts == 2 && add_v7) {
799  if (!intersect_plane_edge(m_eb_point, m_eb_normal, v0, v7, vIP, distIP)) {
800  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");
801  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");
802  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");
803  add_v7 = 0;
804  }
805  }
806 
807  } // end Path 2
808 
809  // Path 5
810  if (p5.intersected) {
811  if (p5.intersected_start) {
812  m_F1.add_vertex(v2); if(print_F1) AMREX_DEVICE_PRINTF("%s \n","Path 5: v2--v6: add v2 -> F1");
813  m_F4.add_vertex(v2); if(print_F4) AMREX_DEVICE_PRINTF("%s \n","Path 5: v2--v6: add v2 -> F4");
814  m_F6.add_vertex(v2); if(print_F6) AMREX_DEVICE_PRINTF("%s \n","Path 5: v2--v6: add v2 -> F6");
815  m_F7.add_vertex(v2); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 5: v2--v6: add v2 -> F7");
816  } else if (p5.intersected_end) {
817  m_F3.add_vertex(v6); if(print_F3) AMREX_DEVICE_PRINTF("%s \n","Path 5: v2--v6: add v6 -> F3");
818  m_F4.add_vertex(v6); if(print_F4) AMREX_DEVICE_PRINTF("%s \n","Path 5: v2--v6: add v6 -> F4");
819  m_F6.add_vertex(v6); if(print_F6) AMREX_DEVICE_PRINTF("%s \n","Path 5: v2--v6: add v6 -> F6");
820  m_F7.add_vertex(v6); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 5: v2--v6: add v6 -> F7");
821  } else {
822  m_F4.add_vertex(p5.vIP); if(print_F4) AMREX_DEVICE_PRINTF("%s \n","Path 5: v2--v6: add vIP -> F4");
823  m_F6.add_vertex(p5.vIP); if(print_F6) AMREX_DEVICE_PRINTF("%s \n","Path 5: v2--v6: add vIP -> F6");
824  m_F7.add_vertex(p5.vIP); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 5: v2--v6: add vIP -> F7");
825  }
826  } // end Path 5
827 
828 
829  // Path 3
830  { int cuts(0);
831  for (int i = 0; i < 8; ++i) vertex_intersected[i] = false;
832 
833  if (p3[0].intersected) {
834  if (p3[0].intersected_start) {
835  m_F1.add_vertex(v0); if(print_F1) AMREX_DEVICE_PRINTF("%s \n","Path 3: v0--v3: add v0 -> F1");
836  m_F4.add_vertex(v0); if(print_F4) AMREX_DEVICE_PRINTF("%s \n","Path 3: v0--v3: add v0 -> F4");
837  m_F5.add_vertex(v0); if(print_F5) AMREX_DEVICE_PRINTF("%s \n","Path 3: v0--v3: add v0 -> F5");
838  m_F7.add_vertex(v0); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 3: v0--v3: add v0 -> F7");
839  if (!vertex_intersected[0]) {
840  vertex_intersected[0] = true;
841  ++cuts;
842  }
843  } else if (p3[0].intersected_end) {
844  m_F3.add_vertex(v3); if(print_F3) AMREX_DEVICE_PRINTF("%s \n","Path 3: v0--v3: add v3 -> F3");
845  m_F4.add_vertex(v3); if(print_F4) AMREX_DEVICE_PRINTF("%s \n","Path 3: v0--v3: add v3 -> F4");
846  m_F5.add_vertex(v3); if(print_F5) AMREX_DEVICE_PRINTF("%s \n","Path 3: v0--v3: add v3 -> F5");
847  m_F7.add_vertex(v3); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 3: v0--v3: add v3 -> F7");
848  if (!vertex_intersected[3]) {
849  vertex_intersected[3] = true;
850  ++cuts;
851  }
852  } else {
853  m_F4.add_vertex(p3[0].vIP); if(print_F4) AMREX_DEVICE_PRINTF("%s \n","Path 3: v0--v3: add vIP -> F4");
854  m_F5.add_vertex(p3[0].vIP); if(print_F5) AMREX_DEVICE_PRINTF("%s \n","Path 3: v0--v3: add vIP -> F5");
855  m_F7.add_vertex(p3[0].vIP); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 3: v0--v3: add vIP -> F7");
856  ++cuts;
857  }
858 
859  } // P3: v0--v3
860 
861  if (cuts%2 == 0) {
862  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");
863  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");
864  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");
865  }
866 
867  if (p3[1].intersected) {
868  if (p3[1].intersected_start) {
869  m_F3.add_vertex(v3); if(print_F3) AMREX_DEVICE_PRINTF("%s \n","Path 3: v3--v6: add v3 -> F3");
870  m_F4.add_vertex(v3); if(print_F4) AMREX_DEVICE_PRINTF("%s \n","Path 3: v3--v6: add v3 -> F4");
871  m_F5.add_vertex(v3); if(print_F5) AMREX_DEVICE_PRINTF("%s \n","Path 3: v3--v6: add v3 -> F5");
872  m_F7.add_vertex(v3); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 3: v3--v6: add v3 -> F7");
873  if (!vertex_intersected[3]) {
874  vertex_intersected[3] = true;
875  ++cuts;
876  }
877  } else if (p3[1].intersected_end) {
878  m_F3.add_vertex(v6); if(print_F3) AMREX_DEVICE_PRINTF("%s \n","Path 3: v3--v6: add v6 -> F3");
879  m_F4.add_vertex(v6); if(print_F4) AMREX_DEVICE_PRINTF("%s \n","Path 3: v3--v6: add v6 -> F4");
880  m_F6.add_vertex(v6); if(print_F6) AMREX_DEVICE_PRINTF("%s \n","Path 3: v3--v6: add v6 -> F6");
881  m_F7.add_vertex(v6); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 3: v3--v6: add v6 -> F7");
882  if (!vertex_intersected[6]) {
883  vertex_intersected[6] = true;
884  ++cuts;
885  }
886  } else {
887  m_F3.add_vertex(p3[1].vIP); if(print_F3) AMREX_DEVICE_PRINTF("%s \n","Path 3: v3--v6: add vIP -> F3");
888  m_F4.add_vertex(p3[1].vIP); if(print_F4) AMREX_DEVICE_PRINTF("%s \n","Path 3: v3--v6: add vIP -> F4");
889  m_F7.add_vertex(p3[1].vIP); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 3: v3--v6: add vIP -> F7");
890  ++cuts;
891  }
892 
893  } // P3: v3--v6
894 
895  if (cuts%2 == 0) {
896  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");
897  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");
898  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");
899  }
900 
901  if (p3[2].intersected) {
902  if (p3[2].intersected_start) {
903  m_F3.add_vertex(v6); if(print_F3) AMREX_DEVICE_PRINTF("%s \n","Path 3: v6--v7: add v6 -> F3");
904  m_F4.add_vertex(v6); if(print_F4) AMREX_DEVICE_PRINTF("%s \n","Path 3: v6--v7: add v6 -> F4");
905  m_F6.add_vertex(v6); if(print_F6) AMREX_DEVICE_PRINTF("%s \n","Path 3: v6--v7: add v6 -> F6");
906  m_F7.add_vertex(v6); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 3: v6--v7: add v6 -> F7");
907  if (!vertex_intersected[6]) {
908  vertex_intersected[6] = true;
909  ++cuts;
910  }
911  } else if (p3[2].intersected_end) {
912  m_F2.add_vertex(v7); if(print_F2) AMREX_DEVICE_PRINTF("%s \n","Path 3: v6--v7: add v7 -> F2");
913  m_F3.add_vertex(v7); if(print_F3) AMREX_DEVICE_PRINTF("%s \n","Path 3: v6--v7: add v7 -> F3");
914  m_F6.add_vertex(v7); if(print_F6) AMREX_DEVICE_PRINTF("%s \n","Path 3: v6--v7: add v7 -> F6");
915  m_F7.add_vertex(v7); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 3: v6--v7: add v7 -> F7");
916  if (!vertex_intersected[7]) {
917  vertex_intersected[7] = true;
918  ++cuts;
919  }
920  } else {
921  m_F3.add_vertex(p3[2].vIP); if(print_F3) AMREX_DEVICE_PRINTF("%s \n","Path 3: v6--v7: add vIP -> F3");
922  m_F6.add_vertex(p3[2].vIP); if(print_F6) AMREX_DEVICE_PRINTF("%s \n","Path 3: v6--v7: add vIP -> F6");
923  m_F7.add_vertex(p3[2].vIP); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 3: v6--v7: add vIP -> F7");
924  ++cuts;
925  }
926  } // P3: v6--v7
927 
928  if (cuts == 2 && add_v7) {
929  if (!intersect_plane_edge(m_eb_point, m_eb_normal, v0, v7, vIP, distIP)) {
930  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");
931  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");
932  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");
933  add_v7 = 0;
934  }
935  }
936 
937  } // end Path 3
938 
939  // Path 6
940  if (p6.intersected) {
941  if (p6.intersected_start) {
942  m_F3.add_vertex(v3); if(print_F3) AMREX_DEVICE_PRINTF("%s \n","Path 6: v3--v4: add v3 -> F3");
943  m_F4.add_vertex(v3); if(print_F4) AMREX_DEVICE_PRINTF("%s \n","Path 6: v3--v4: add v3 -> F4");
944  m_F5.add_vertex(v3); if(print_F5) AMREX_DEVICE_PRINTF("%s \n","Path 6: v3--v4: add v3 -> F5");
945  m_F7.add_vertex(v3); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 6: v3--v4: add v3 -> F7");
946  } else if (p6.intersected_end) {
947  m_F2.add_vertex(v4); if(print_F2) AMREX_DEVICE_PRINTF("%s \n","Path 6: v3--v4: add v4 -> F2");
948  m_F3.add_vertex(v4); if(print_F3) AMREX_DEVICE_PRINTF("%s \n","Path 6: v3--v4: add v4 -> F3\n");
949  m_F5.add_vertex(v4); if(print_F5) AMREX_DEVICE_PRINTF("%s \n","Path 6: v3--v4: add v4 -> F5");
950  m_F7.add_vertex(v4); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 6: v3--v4: add v4 -> F7");
951  } else {
952  m_F3.add_vertex(p6.vIP); if(print_F3) AMREX_DEVICE_PRINTF("%s \n","Path 6: v3--v4: add vIP -> F3");
953  m_F5.add_vertex(p6.vIP); if(print_F5) AMREX_DEVICE_PRINTF("%s \n","Path 6: v3--v4: add vIP -> F5");
954  m_F7.add_vertex(p6.vIP); if(print_F7) AMREX_DEVICE_PRINTF("%s \n","Path 6: v3--v4: add vIP -> F7");
955  }
956  } // end Path 6
957 
958  if (print_F1 || print_F2 || print_F3 || print_F4 || print_F5 || print_F6 || print_F7) {
959  AMREX_DEVICE_PRINTF("%s \n"," ");
960  }
961 
962 }
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)
Intersect a plane with a line segment.
Definition: ERF_EBCutCell.H:30
polygon_::add_vertex
AMREX_GPU_HOST_DEVICE void add_vertex(amrex::RealVect const &a_v)
Add a unique vertex to the polygon.
Definition: ERF_EBPolygon.H:63
Here is the call graph for this function:

◆ centBoun()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::RealVect eb_cut_cell_::centBoun ( ) const
inlinenoexcept

Return the EB boundary centroid.

226  {
227  amrex::RealVect cent{zero,zero,zero};
228  if (m_flag.isSingleValued()) {
229  cent = m_F7.get_centroid();
230  }
231  return cent;
232  }
polygon_::get_centroid
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::RealVect get_centroid() const noexcept
Return the polygon centroid from its sub-triangulation.
Definition: ERF_EBPolygon.H:205

Referenced by eb_aux_::define().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ centHi()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::RealVect eb_cut_cell_::centHi ( int const  idir) const
inlinenoexcept

Return the high-side face centroid in one coordinate direction.

Parameters
idirCoordinate direction.
205  {
206  AMREX_ASSERT( idir >=0 && idir < AMREX_SPACEDIM );
207  amrex::RealVect cent = m_cellface_cent[ m_hi_faces_id[idir] ];
208  if (m_flag.isCovered() || m_flag.isRegular() ) {
209  // Default cent
210  } else {
211  amrex::Real const area_R = m_hi_faces[idir]->area();
212  if (amrex::almostEqual(area_R,zero) || amrex::almostEqual(area_R, m_rbox_area[idir]) ){
213  // Default cent
214  } else {
215  amrex::RealVect cent_O = cent;
216  amrex::RealVect cent_R = m_hi_faces[idir]->get_centroid();
217  amrex::Real const area_C = m_rbox_area[idir] - area_R;
218  cent = m_invert * cent_R + (one-m_invert) * ((one+area_R/area_C)*cent_O - area_R/area_C*cent_R);
219  }
220  }
221  return cent;
222  }
eb_cut_cell_::m_hi_faces_id
amrex::Array< int const, 3 > m_hi_faces_id
Definition: ERF_EBCutCell.H:314

Referenced by eb_aux_::define().

Here is the caller graph for this function:

◆ centLo()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::RealVect eb_cut_cell_::centLo ( int const  idir) const
inlinenoexcept

Return the low-side face centroid in one coordinate direction.

Parameters
idirCoordinate direction.
181  {
182  AMREX_ASSERT( idir >=0 && idir < AMREX_SPACEDIM );
183  amrex::RealVect cent = m_cellface_cent[ m_lo_faces_id[idir] ];
184  if (m_flag.isCovered() || m_flag.isRegular() ) {
185  // Default cent
186  } else {
187  amrex::Real const area_R = m_lo_faces[idir]->area();
188  if (amrex::almostEqual(area_R,zero) || amrex::almostEqual(area_R, m_rbox_area[idir]) ){
189  // Default cent
190  } else {
191  amrex::RealVect cent_O = cent;
192  amrex::RealVect cent_R = m_lo_faces[idir]->get_centroid();
193  amrex::Real const area_C = m_rbox_area[idir] - area_R;
194  cent = m_invert * cent_R + (one-m_invert) * ((one+area_R/area_C)*cent_O - area_R/area_C*cent_R);
195  }
196  }
197  return cent;
198  }
eb_cut_cell_::m_lo_faces_id
amrex::Array< int const, 3 > m_lo_faces_id
Definition: ERF_EBCutCell.H:313

Referenced by eb_aux_::define().

Here is the caller graph for this function:

◆ centVol()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::RealVect eb_cut_cell_::centVol ( ) const
inlinenoexcept

Return the volume centroid of the reconstructed cell.

246  {
247  amrex::RealVect vcent = m_cell_cent;
248  if (m_flag.isSingleValued()) {
249  amrex::Real xm = m_rbox.lo(0);
250  amrex::Real xp = m_rbox.hi(0);
251  amrex::Real ym = m_rbox.lo(1);
252  amrex::Real yp = m_rbox.hi(1);
253  amrex::Real zm = m_rbox.lo(2);
254  amrex::Real zp = m_rbox.hi(2);
255 
256  amrex::Real axm = areaLo(0);
257  amrex::Real axp = areaHi(0);
258  amrex::Real aym = areaLo(1);
259  amrex::Real ayp = areaHi(1);
260  amrex::Real azm = areaLo(2);
261  amrex::Real azp = areaHi(2);
262 
263  amrex::Real barea = m_F7.area();
264  amrex::RealVect bcent = m_F7.get_centroid();
265 
266  amrex::Real vol = volume();
267 
268  vcent[0] = myhalf * ( - xm * xm * axm + xp * xp * axp + bcent[0] * bcent[0] * m_eb_normal[0] * barea ) / vol;
269  vcent[1] = myhalf * ( - ym * ym * aym + yp * yp * ayp + bcent[1] * bcent[1] * m_eb_normal[1] * barea ) / vol;
270  vcent[2] = myhalf * ( - zm * zm * azm + zp * zp * azp + bcent[2] * bcent[2] * m_eb_normal[2] * barea ) / vol;
271  }
272  return vcent;
273  }
eb_cut_cell_::volume
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real volume() const
Return the physical volume of the reconstructed cell.
Definition: ERF_EBCutCell.H:124
eb_cut_cell_::areaLo
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real areaLo(int const idir) const noexcept
Return the low-side face area in one coordinate direction.
Definition: ERF_EBCutCell.H:149
eb_cut_cell_::areaHi
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real areaHi(int const idir) const noexcept
Return the high-side face area in one coordinate direction.
Definition: ERF_EBCutCell.H:162

Referenced by eb_aux_::define().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ debug()

void eb_cut_cell_::debug ( int const  a_face = -1)

Print cut-cell diagnostics on CPU builds.

Parameters
a_faceFace identifier to report, or -1 for all faces.
13  {
14 
15  amrex::Print() << "\n\nDEBUG THIS "
16  << "--------------------------------------"
17  << "\nisCovered? " << isCovered()
18  << "\nisRegular? " << isRegular()
19  << "\nisSingleValued? " << isSingleValued()
20  << "\n";
21 
22  if ( isCovered() || isRegular() ) { return; }
23 
24  amrex::RealVect v0(m_rbox.lo(0), m_rbox.lo(1), m_rbox.lo(2));
25  amrex::RealVect v7(m_rbox.hi(0), m_rbox.hi(1), m_rbox.hi(2));
26 
27  amrex::Print() << "\n"
28  << "lo: " << v0 << '\n'
29  << "hi: " << v7 << '\n'
30  << "p: " << m_eb_point << '\n'
31  << "n: " << m_eb_normal << '\n'
32  << "invert? " << m_invert << "\n\n";
33 
34  amrex::Print() << "Edge intersections:\n";
35  calc_edge_intersections();
36  amrex::Print() << '\n';
37 
38  if ( a_face == -1 || a_face == 1 ) { m_F1.report(1, v0); }
39  if ( a_face == -1 || a_face == 2 ) { m_F2.report(2, v0); }
40  if ( a_face == -1 || a_face == 3 ) { m_F3.report(3, v0); }
41  if ( a_face == -1 || a_face == 4 ) { m_F4.report(4, v0); }
42  if ( a_face == -1 || a_face == 5 ) { m_F5.report(5, v0); }
43  if ( a_face == -1 || a_face == 6 ) { m_F6.report(6, v0); }
44  if ( a_face == -1 || a_face == 7 ) { m_F7.report(7, v0); }
45 }
eb_cut_cell_::isRegular
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool isRegular() const noexcept
Return true if the reconstructed cell is regular.
Definition: ERF_EBCutCell.H:91
eb_cut_cell_::isSingleValued
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool isSingleValued() const noexcept
Return true if the reconstructed cell is single-valued.
Definition: ERF_EBCutCell.H:95
eb_cut_cell_::isCovered
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool isCovered() const noexcept
Return true if the reconstructed cell is covered.
Definition: ERF_EBCutCell.H:87
polygon_::report
void report(int const a_id, amrex::RealVect a_v0)
Print a diagnostic report for the polygon on CPU builds.
Definition: ERF_EBPolygon.H:238

Referenced by eb_aux_::define().

Here is the caller graph for this function:

◆ isCovered()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool eb_cut_cell_::isCovered ( ) const
inlinenoexcept

Return true if the reconstructed cell is covered.

87 { return m_flag.isCovered(); }

Referenced by eb_aux_::define().

Here is the caller graph for this function:

◆ isRegular()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool eb_cut_cell_::isRegular ( ) const
inlinenoexcept

Return true if the reconstructed cell is regular.

91 { return m_flag.isRegular(); }

Referenced by eb_aux_::define().

Here is the caller graph for this function:

◆ isSingleValued()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool eb_cut_cell_::isSingleValued ( ) const
inlinenoexcept

Return true if the reconstructed cell is single-valued.

95 { return m_flag.isSingleValued(); }

◆ normBoun()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::RealVect eb_cut_cell_::normBoun ( ) const
inlinenoexcept

Return the EB boundary normal.

236  {
237  amrex::RealVect normal{zero,zero,zero};
238  if (m_flag.isSingleValued()) {
239  normal = m_eb_normal;
240  }
241  return normal;
242  }

◆ set_covered()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void eb_cut_cell_::set_covered ( )
inline

Mark the reconstructed cell as covered.

99  {
100  m_flag.setCovered();
101  }

◆ set_covered_regular_cell_vertices()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void eb_cut_cell_::set_covered_regular_cell_vertices ( )
private

Populate polygon vertices for cells that are fully covered or regular.

969 {
970  using namespace amrex;
971 
972  RealVect v0(m_rbox.lo(0), m_rbox.lo(1), m_rbox.lo(2));
973  RealVect v1(m_rbox.hi(0), m_rbox.lo(1), m_rbox.lo(2));
974  RealVect v2(m_rbox.lo(0), m_rbox.hi(1), m_rbox.lo(2));
975  RealVect v3(m_rbox.lo(0), m_rbox.lo(1), m_rbox.hi(2));
976  RealVect v4(m_rbox.hi(0), m_rbox.lo(1), m_rbox.hi(2));
977  RealVect v5(m_rbox.hi(0), m_rbox.hi(1), m_rbox.lo(2));
978  RealVect v6(m_rbox.lo(0), m_rbox.hi(1), m_rbox.hi(2));
979  RealVect v7(m_rbox.hi(0), m_rbox.hi(1), m_rbox.hi(2));
980 
981  // Add vertices in the order of outward normal vector
982 
983  m_F1.add_vertex(v0);
984  m_F1.add_vertex(v2);
985  m_F1.add_vertex(v5);
986  m_F1.add_vertex(v1);
987 
988  m_F2.add_vertex(v1);
989  m_F2.add_vertex(v5);
990  m_F2.add_vertex(v7);
991  m_F2.add_vertex(v4);
992 
993  m_F3.add_vertex(v3);
994  m_F3.add_vertex(v4);
995  m_F3.add_vertex(v7);
996  m_F3.add_vertex(v6);
997 
998  m_F4.add_vertex(v0);
999  m_F4.add_vertex(v3);
1000  m_F4.add_vertex(v6);
1001  m_F4.add_vertex(v2);
1002 
1003  m_F5.add_vertex(v0);
1004  m_F5.add_vertex(v1);
1005  m_F5.add_vertex(v4);
1006  m_F5.add_vertex(v3);
1007 
1008  m_F6.add_vertex(v2);
1009  m_F6.add_vertex(v5);
1010  m_F6.add_vertex(v7);
1011  m_F6.add_vertex(v6);
1012 
1013 }
Here is the call graph for this function:

◆ set_regular()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void eb_cut_cell_::set_regular ( )
inline

Mark the reconstructed cell as regular and fill regular face areas.

105  {
106 
107  m_flag.setRegular();
108 
109  m_F1.set_area( m_rbox.length(0)*m_rbox.length(1) );
110  m_F3.set_area( m_rbox.length(0)*m_rbox.length(1) );
111 
112  m_F2.set_area( m_rbox.length(1)*m_rbox.length(2) );
113  m_F4.set_area( m_rbox.length(1)*m_rbox.length(2) );
114 
115  m_F5.set_area( m_rbox.length(0)*m_rbox.length(2) );
116  m_F6.set_area( m_rbox.length(0)*m_rbox.length(2) );
117 
118  m_F7.set_area(amrex::Real(0));
119 
120  }
polygon_::set_area
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void set_area(amrex::Real const &a_area)
Set the polygon area directly.
Definition: ERF_EBPolygon.H:87
Here is the call graph for this function:

◆ volume()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real eb_cut_cell_::volume ( ) const
inline

Return the physical volume of the reconstructed cell.

124  {
125 
126  if (m_flag.isCovered() ) { return zero; }
127  if (m_flag.isRegular() ) { return m_rbox.volume(); }
128 
129  amrex::Real volume(zero);
130 
131  amrex::Real const* lo = m_rbox.lo();
132  amrex::RealVect v0(lo[0], lo[1], lo[2]);
133 
134  if (m_F2.ok() ) { volume += m_F2.area() * m_F2.distance(v0); }
135  if (m_F3.ok() ) { volume += m_F3.area() * m_F3.distance(v0); }
136  if (m_F6.ok() ) { volume += m_F6.area() * m_F6.distance(v0); }
137  if (m_F7.ok() ) { volume += m_F7.area() * m_F7.distance(v0); }
138 
139  volume /= three;
140 
141  return m_invert*volume + (one-m_invert)*(m_rbox.volume()-volume);
142  }
three
constexpr amrex::Real three
Definition: ERF_Constants.H:11
polygon_::ok
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int ok() const noexcept
Return whether the polygon has been defined with a valid area state.
Definition: ERF_EBPolygon.H:182
polygon_::distance
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real distance(amrex::RealVect const &a_point) const noexcept
Return the perpendicular distance from this polygon to a point.
Definition: ERF_EBPolygon.H:197

Referenced by centVol(), and eb_aux_::define().

Here is the call graph for this function:
Here is the caller graph for this function:

Member Data Documentation

◆ m_cell_cent

amrex::RealVect eb_cut_cell_::m_cell_cent
private

Referenced by centVol().

◆ m_cellface_cent

amrex::Array<amrex::RealVect,m_max_faces> eb_cut_cell_::m_cellface_cent
private

Referenced by centHi(), and centLo().

◆ m_dx

amrex::RealVect eb_cut_cell_::m_dx
private

◆ m_eb_normal

amrex::RealVect const eb_cut_cell_::m_eb_normal
private

Referenced by calc_edge_intersections(), centVol(), and normBoun().

◆ m_eb_point

amrex::RealVect const eb_cut_cell_::m_eb_point
private

Referenced by calc_edge_intersections().

◆ m_F1

polygon_ eb_cut_cell_::m_F1
private

Referenced by calc_edge_intersections(), set_covered_regular_cell_vertices(), and set_regular().

◆ m_F2

polygon_ eb_cut_cell_::m_F2
private

Referenced by calc_edge_intersections(), set_covered_regular_cell_vertices(), set_regular(), and volume().

◆ m_F3

polygon_ eb_cut_cell_::m_F3
private

Referenced by calc_edge_intersections(), set_covered_regular_cell_vertices(), set_regular(), and volume().

◆ m_F4

polygon_ eb_cut_cell_::m_F4
private

Referenced by calc_edge_intersections(), set_covered_regular_cell_vertices(), and set_regular().

◆ m_F5

polygon_ eb_cut_cell_::m_F5
private

Referenced by calc_edge_intersections(), set_covered_regular_cell_vertices(), and set_regular().

◆ m_F6

polygon_ eb_cut_cell_::m_F6
private

Referenced by calc_edge_intersections(), set_covered_regular_cell_vertices(), set_regular(), and volume().

◆ m_F7

polygon_ eb_cut_cell_::m_F7
private

Referenced by areaBoun(), calc_edge_intersections(), centBoun(), centVol(), set_regular(), and volume().

◆ m_flag

amrex::EBCellFlag eb_cut_cell_::m_flag
private

Referenced by areaHi(), areaLo(), centBoun(), centHi(), centLo(), centVol(), isCovered(), isRegular(), isSingleValued(), normBoun(), set_covered(), set_regular(), and volume().

◆ m_hi_faces

amrex::Array<polygon_ const* const,3> eb_cut_cell_::m_hi_faces {&m_F2, &m_F6, &m_F3}
private

Referenced by areaHi(), and centHi().

◆ m_hi_faces_id

amrex::Array<int const,3> eb_cut_cell_::m_hi_faces_id {1,5,2}
private

Referenced by centHi().

◆ m_invert

amrex::Real eb_cut_cell_::m_invert
private

Referenced by areaHi(), areaLo(), centHi(), centLo(), and volume().

◆ m_lo_faces

amrex::Array<polygon_ const* const,3> eb_cut_cell_::m_lo_faces {&m_F4, &m_F5, &m_F1}
private

Referenced by areaLo(), and centLo().

◆ m_lo_faces_id

amrex::Array<int const,3> eb_cut_cell_::m_lo_faces_id {3,4,0}
private

Referenced by centLo().

◆ m_max_faces

constexpr int eb_cut_cell_::m_max_faces = 6
staticconstexprprivate

◆ m_rbox

amrex::RealBox const eb_cut_cell_::m_rbox
private

Referenced by calc_edge_intersections(), centVol(), set_covered_regular_cell_vertices(), set_regular(), and volume().

◆ m_rbox_area

amrex::RealVect eb_cut_cell_::m_rbox_area
private

Referenced by areaHi(), areaLo(), centHi(), and centLo().

The documentation for this class was generated from the following files: