ERF_EBPolygon.H

Go to the documentation of this file.

/**
 * \file ERF_EBPolygon.H
 * \brief Defines polygon geometry utilities for EB cut-cell reconstruction.
 */
#ifndef ERF_EB_POLYGON_H_
#define ERF_EB_POLYGON_H_
 
#include <AMReX_REAL.H>
#include <AMReX_RealVect.H>
 
#include <cmath>
#include "ERF_NumericalConstants.H"
 
/**
 * \brief Polygon helper used to accumulate cut-cell face geometry.
 *
 * The polygon stores up to six vertices, sorts them in a local plane, and
 * computes area, centroid, distance, and normal data used by eb_cut_cell_.
 */
class polygon_ {
 
  public:
 
    /**
     * \brief Construct a cell-face polygon clipped by an EB plane.
     * \param a_point Point on the EB plane.
     * \param a_normal Unit normal for the EB plane.
     */
    AMREX_GPU_HOST_DEVICE
    polygon_ ( amrex::RealVect a_point,
               amrex::RealVect a_normal )
      : m_cell_face(1)
      , m_eb_point(a_point)
      , m_eb_normal(a_normal)
      , m_defined(0)
      , m_num_vertices(0)
      , m_area(zero)
      , m_sorted(0)
      , m_vertices({amrex::RealVect(zero), amrex::RealVect(zero), amrex::RealVect(zero),
                  amrex::RealVect(zero), amrex::RealVect(zero), amrex::RealVect(zero)})
      , m_zdir(zero)
    {}
 
    //! Construct an empty EB-boundary polygon.
    AMREX_GPU_HOST_DEVICE
    polygon_ ( )
      : m_cell_face(0)
      , m_eb_point(zero)
      , m_eb_normal(zero)
      , m_defined(0)
      , m_num_vertices(0)
      , m_area(zero)
      , m_sorted(0)
      , m_vertices({amrex::RealVect(zero), amrex::RealVect(zero), amrex::RealVect(zero),
                  amrex::RealVect(zero), amrex::RealVect(zero), amrex::RealVect(zero)})
      , m_zdir(zero)
    {}
 
    /**
     * \brief Add a unique vertex to the polygon.
     * \param a_v Vertex coordinates.
     */
    AMREX_GPU_HOST_DEVICE
    void add_vertex ( amrex::RealVect const& a_v ) {
      for ( int i(0); i<m_num_vertices; ++i) {
        if ( amrex::almostEqual(m_vertices[i][0],a_v[0]) &&
             amrex::almostEqual(m_vertices[i][1],a_v[1]) &&
             amrex::almostEqual(m_vertices[i][2],a_v[2]) ) {
          return;
        }
      }
      AMREX_ASSERT( m_num_vertices < m_max_vertices );
      m_vertices[m_num_vertices] = a_v;
      m_num_vertices++;
    }
 
    //! Return the number of stored vertices.
    AMREX_GPU_HOST_DEVICE
    int get_num_vertices ( ) {
      return m_num_vertices;
    }
 
    /**
     * \brief Set the polygon area directly.
     * \param a_area Area value assigned to the polygon.
     */
    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    void set_area ( amrex::Real const& a_area ) { m_area = a_area; }
 
    /**
     * \brief Finalize polygon ordering and compute its area.
     *
     * Vertices are sorted counter-clockwise in the polygon plane before area is
     * computed by a triangle fan.
     */
    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    void define () {
 
      AMREX_ALWAYS_ASSERT( m_defined == 0 );  // TODO ---------------------------- remove ALWAYS
 
      m_defined = 1;
 
      // We need at least 3 vertices.
      if (m_num_vertices < 3) { return; }
 
      // Check to see if the vertices of the face are inside or outside of the plane.
      // If they are outside, this face doesn't belong to the volume.
      if ( m_cell_face ) {
        for ( int i(0); i<m_num_vertices; ++i) {
          if ((m_eb_normal.dotProduct(m_vertices[i] - m_eb_point)) >= zero) {
          } else {
          }
        }
      }
 
      // Calculate the centroid of the polygon.
      amrex::RealVect centroid = get_centroid();
 
      // Shift the vertices relative to the centroid
      amrex::Array<amrex::RealVect,m_max_vertices> vertex_cent;
      for ( int i(0); i<m_num_vertices; ++i) {
        vertex_cent[i] = m_vertices[i] - centroid;
      }
 
      // Compute the normal vector m_zdir by cross product of two vectors in vertex_cent.
      // Choose the vectors with the largest cross-product magnitude.
      amrex::RealVect v_normal;
      amrex::Real max_norm2 = -one;
      const amrex::RealVect& vertex_cent_0 = vertex_cent[0];
 
      for (int i = 1; i < m_num_vertices; ++i)
      {
        amrex::RealVect vi = vertex_cent[i];
        amrex::RealVect v_cross = vertex_cent_0.crossProduct(vi);
        amrex::Real n2 = v_cross.radSquared();
        if (n2 > max_norm2)
        {
          max_norm2 = n2;
          v_normal = v_cross;
        }
      }
      if (!amrex::almostEqual(max_norm2, zero))
      {
          v_normal /= std::sqrt(max_norm2);
      }
      m_zdir = v_normal;
 
      //
 
      m_theta[0] = zero;
      for ( int i(1); i<m_num_vertices; ++i) {
 
        m_theta[i] = std::atan2(m_zdir.dotProduct( vertex_cent[0].crossProduct(vertex_cent[i]) ),
                              vertex_cent[0].dotProduct(vertex_cent[i]));
 
        m_theta[i] += ((m_theta[i] >= zero) ? zero : two*PI);
      }
 
      // Sort counter clockwise based on theta.
      for (int i(0); i<m_num_vertices; ++i) {
        for (int j(0); j < m_num_vertices-i-1; ++j) {
          if ( m_theta[j] > m_theta[j+1] ) {
            amrex::Swap(m_theta[j], m_theta[j+1]);
            amrex::Swap(m_vertices[j], m_vertices[j+1]);
            amrex::Swap(vertex_cent[j], vertex_cent[j+1]);
          }
        } // j-loop
      } // i-loop
      m_sorted = 1;
 
      // Compute areas of triangles
 
      for (int i(0); i<m_num_vertices; ++i) {
        int const j( (i+1 == m_num_vertices) ? 0 : i+1 );
        amrex::RealVect vi_cross_vj = vertex_cent[i].crossProduct(vertex_cent[j]);
        m_area += myhalf*vi_cross_vj.vectorLength();
      }
      AMREX_ALWAYS_ASSERT( m_area > zero ); // <------------------------------- TODO remove ALWAYS
    } // void define
 
    //! Return whether the polygon has been defined with a valid area state.
    [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    int ok ( ) const noexcept
    { return ((m_area > zero || (m_area == zero && m_defined == 1)) ? 1 : 0); }
 
    //! Return the polygon area.
    [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    amrex::Real area ( ) const noexcept {
      AMREX_ALWAYS_ASSERT( ok() ); // <-------------------------------------- TODO remove ALWAYS
      return m_area;
    }
 
    /**
     * \brief Return the perpendicular distance from this polygon to a point.
     * \param a_point Point used in the distance calculation.
     */
    [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    amrex::Real distance ( amrex::RealVect const& a_point ) const noexcept {
      AMREX_ALWAYS_ASSERT( m_defined == 1 );  // <--------------------------- TODO remove ALWAYS
      amrex::RealVect x0 = a_point - m_vertices[0];
      return amrex::Math::abs(x0.dotProduct(m_zdir));
    }
 
    //! Return the polygon centroid from its sub-triangulation.
    [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    amrex::RealVect get_centroid ( ) const noexcept {
      amrex::RealVect cent(zero);
      if (m_num_vertices==3) {
        cent = (m_vertices[0] + m_vertices[1] + m_vertices[2]) / three;
      } else {
        // Loop over sub-triangles
        amrex::Real area(zero);
        for ( int i(0); i<m_num_vertices-2; ++i) {
          amrex::RealVect const v0 (m_vertices[i+1] - m_vertices[0]);
          amrex::RealVect const v1 (m_vertices[i+2] - m_vertices[0]);
          amrex::RealVect v0_cross_v1 = v0.crossProduct(v1);
          amrex::Real area_tri = myhalf * v0_cross_v1.vectorLength();
          amrex::RealVect cent_tri = (m_vertices[0] + m_vertices[i+1] + m_vertices[i+2]) / three;
          area += area_tri;
          cent += area_tri * cent_tri;
        }
        cent = cent / area;
      }
      return cent;
    }
 
    //! Return the unit normal vector of the polygon plane.
    [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    amrex::RealVect normal () const noexcept {
      return m_zdir;
    }
 
    /**
     * \brief Print a diagnostic report for the polygon on CPU builds.
     * \param a_id Face identifier used in the diagnostic output.
     * \param a_v0 Reference point used to compute the reported distance.
     */
#ifndef AMREX_USE_GPU
    void report ( int const a_id, amrex::RealVect a_v0 )
#else
    void report ( int const /*a_id*/, amrex::RealVect /*a_v0*/ )
#endif
    {
 
#ifndef AMREX_USE_GPU
      amrex::RealVect centroid = get_centroid();
      amrex::Print() << "Face " << a_id
        << " -------------------------------------------\n"
        << "\nok? " << (ok() ? "yes" : "no") << "\n\n";
      for (int i(0); i<m_num_vertices; ++i) {
        amrex::Print() << "v" << i << ":  " << m_vertices[i] << '\n';
      }
      amrex::Print() << "\nvc: " << "  " << centroid << "\n";
 
      amrex::Real const dist = distance( a_v0 );
      // amrex::Real const vol  = ok() ? m_area * dist : zero;
 
      amrex::Print() << "\narea:     " << m_area
                     << "\ndistance: " << dist
                     << "\nvolume:   " << dist*m_area
      << "\n==================================================\n\n";
#endif
    }
 
    /**
     * \brief Print polygon vertices with a face identifier.
     * \param a_id Face identifier used in the diagnostic output.
     */
    void debug( int const a_id ) {
      amrex::Print() << "EBPolygon: id = " << a_id << ", m_num_vertices = " << m_num_vertices << '\n';
      for (int i(0); i<m_num_vertices; ++i) {
        amrex::Print() << "EBPolygon: v" << i << ":  " << m_vertices[i] << '\n';
      }
    }
    //! Print polygon vertices without a face identifier.
    void debug() {
      debug(-1);
    }
 
  private:
 
    static int constexpr m_max_vertices = 6;
 
    int const m_cell_face;
 
    amrex::RealVect const m_eb_point;
    amrex::RealVect const m_eb_normal;
 
    int m_defined;
 
    int m_num_vertices;
 
    amrex::Real m_area;
 
    int m_sorted;
 
    amrex::Array<amrex::RealVect,m_max_vertices> m_vertices;
 
    amrex::GpuArray<amrex::Real,m_max_vertices> m_theta;
 
    amrex::RealVect m_zdir; // normal to polygon
};
 
 
#endif