AMReX_EB2_IF_Polynomial.H

Go to the documentation of this file.

#ifndef AMREX_EB2_IF_POLYNOMIAL_H_
#define AMREX_EB2_IF_POLYNOMIAL_H_
#include <AMReX_Config.H>
 
#include <AMReX_EB2_IF_Base.H>
#include <AMReX_Array.H>
#include <AMReX_Vector.H>
#include <AMReX_IntVect.H>
#include <memory>
#include <cmath>
 
// For all implicit functions, >0: body; =0: boundary; <0: fluid
 
namespace amrex::EB2 {
 
/********************************************************************************
 *                                                                              *
 * Represents one term in a general polynomial                                  *
 *                                                                              *
 ********************************************************************************/
class PolyTerm
{
public:
 
    Real coef;
 
    IntVect powers;
};
 
template <unsigned int N>
class PolyIF
    : public GPUable
{
public:
 
    PolyIF (const GpuArray<PolyTerm,N> & a_polynomial, bool a_inside = true)
        : m_polynomial(a_polynomial),
          m_sign( a_inside ? 1.0_rt : -1.0_rt )
    {}
 
    AMREX_GPU_HOST_DEVICE inline
    Real operator() (AMREX_D_DECL(Real x, Real y, Real z)) const noexcept
    {
        Real retval = 0.0_rt;
        for (auto const& term : m_polynomial) {
            retval += term.coef * AMREX_D_TERM(  std::pow(x, term.powers[0]),
                                               * std::pow(y, term.powers[1]),
                                               * std::pow(z, term.powers[2]));
        }
        return m_sign*retval;
    }
 
    inline Real operator() (const RealArray& p) const noexcept {
        return this->operator()(AMREX_D_DECL(p[0],p[1],p[2]));
    }
 
protected:
    GpuArray<PolyTerm,N> m_polynomial;
    Real                 m_sign;
};
 
 
class PolynomialIF  // No GPU support
{
public:
 
    PolynomialIF (const Vector<PolyTerm> & a_polynomial, bool a_inside = true)
        : m_polynomial(a_polynomial),
          m_inside(a_inside),
          m_sign( a_inside ? 1.0_rt : -1.0_rt ),
          m_size(int(m_polynomial.size()))
    {}
 
    ~PolynomialIF () = default;
    PolynomialIF (const PolynomialIF& rhs) = default;
    PolynomialIF (PolynomialIF&& rhs) = default;
    PolynomialIF& operator= (const PolynomialIF& rhs) = delete;
    PolynomialIF& operator= (PolynomialIF&& rhs) = delete;
 
    [[nodiscard]] inline
    Real operator() (AMREX_D_DECL(Real x, Real y, Real z)) const noexcept
    {
        Real retval = 0.0_rt;
        for (auto const& term : m_polynomial) {
            retval += term.coef * AMREX_D_TERM(  std::pow(x, term.powers[0]),
                                               * std::pow(y, term.powers[1]),
                                               * std::pow(z, term.powers[2]));
        }
        return m_sign*retval;
    }
 
    [[nodiscard]] inline Real operator() (const RealArray& p) const noexcept {
        return this->operator()(AMREX_D_DECL(p[0],p[1],p[2]));
    }
 
protected:
    Vector<PolyTerm> m_polynomial;
    bool             m_inside;
    Real             m_sign;
    int              m_size;
};
 
}
 
#endif