Submergence.h

#ifndef FHSIM_BASE_SUBMERGENCE_H
#define FHSIM_BASE_SUBMERGENCE_H
 
#include <algorithm>
#include <cmath>
#include <limits>
 
namespace submergence
{
 
template<class Environment>
double WaveElevationBound(const Environment& environment)
{
   if constexpr (requires { environment.MaxWaveElevation(); })
      return environment.MaxWaveElevation();
   else
      return std::numeric_limits<double>::infinity();
}
 
template<class Environment>
double MediumDensity(Environment& environment, double time, const double pos[3], double waterDensity)
{
   constexpr double kAirWaterThreshold = 100.0; // kg/m^3, between air and any water
   double density                      = 0.0;
   environment.GetDensity(time, pos, density);
   return density > kAirWaterThreshold ? waterDensity : density;
}
 
inline bool SurelySubmerged(double topDepth, double waveBound)
{
   return topDepth > waveBound;
}
 
inline double CircleFraction(double x)
{
   if (x >= 1.0)
      return 1.0;
   if (x <= -1.0)
      return 0.0;
   return 0.5 + (x * std::sqrt(1.0 - x * x) + std::asin(x)) / M_PI;
}
 
inline double CircleFractionIntegral(double x)
{
   if (x >= 1.0)
      return x;
   if (x <= -1.0)
      return 0.0;
   const double c = 1.0 - x * x;
   return 0.5 * x + (x * std::asin(x) + std::sqrt(c) - c * std::sqrt(c) / 3.0) / M_PI;
}
 
inline double CylinderFraction(double s1, double s2, double sectionHeight)
{
   const double low  = std::min(s1, s2);
   const double high = std::max(s1, s2);
   const double half = 0.5 * sectionHeight;
   if (low >= half)
      return 1.0;
   if (high <= -half)
      return 0.0;
   if (!(half > 0.0))
      return high / (high - low); // vertical: the wetted length fraction, here low < 0 < high
   const double x1 = low / half;
   const double x2 = high / half;
   if (x2 - x1 < 1e-6)
      return CircleFraction(0.5 * (x1 + x2)); // horizontal: one cross-section
   return (CircleFractionIntegral(x2) - CircleFractionIntegral(x1)) / (x2 - x1);
}
 
inline double CylinderFractionGradient(double s1, double s2, double sectionHeight, double& dF_ds1, double& dF_ds2, double& dF_dHeight)
{
   dF_ds1 = dF_ds2 = dF_dHeight = 0.0;
   const bool firstIsLow = s1 <= s2;
   const double low      = firstIsLow ? s1 : s2;
   const double high     = firstIsLow ? s2 : s1;
   const double half     = 0.5 * sectionHeight;
   double dLow = 0.0, dHigh = 0.0;
   double fraction;
   if (low >= half) {
      fraction = 1.0;
   } else if (high <= -half) {
      fraction = 0.0;
   } else if (!(half > 0.0)) {
      const double span = high - low;
      fraction          = high / span;
      dHigh             = -low / (span * span);
      dLow              = high / (span * span);
   } else {
      const double x1 = low / half;
      const double x2 = high / half;
      if (x2 - x1 < 1e-6) {
         const double xm    = 0.5 * (x1 + x2);
         fraction           = CircleFraction(xm);
         const double slope = (std::abs(xm) < 1.0) ? 2.0 * std::sqrt(1.0 - xm * xm) / M_PI : 0.0;
         dLow = dHigh = 0.5 * slope / half;
         dF_dHeight   = -0.5 * slope * xm / half;
      } else {
         const double span = x2 - x1;
         fraction          = (CircleFractionIntegral(x2) - CircleFractionIntegral(x1)) / span;
         const double dx1  = (fraction - CircleFraction(x1)) / span;
         const double dx2  = (CircleFraction(x2) - fraction) / span;
         dLow              = dx1 / half;
         dHigh             = dx2 / half;
         dF_dHeight        = -0.5 * (dx1 * x1 + dx2 * x2) / half;
      }
   }
   dF_ds1 = firstIsLow ? dLow : dHigh;
   dF_ds2 = firstIsLow ? dHigh : dLow;
   return fraction;
}
 
inline double CylinderHalfHeight(double axisDz, double diameter, double cosTheta)
{
   return 0.5 * (std::abs(axisDz) + diameter * std::sqrt(std::max(0.0, 1.0 - cosTheta * cosTheta)));
}
 
} // namespace submergence
 
#endif