GaussPointDrag.h

#pragma once
 
#include "WetCylinderWater.h"
 
#include <fhsim_environment/ElementFlowInteraction.h>
#include <fhsim_environment/hydrodynamics/WetCylinderLaw.h>
 
#include <Eigen/Eigen>
 
#include <algorithm>
#include <array>
#include <cmath>
#include <vector>
 
namespace wet_cylinder_water
{
using environment::wet_cylinder::EndWater; 
 
struct GaussDragLoad
{
   Eigen::Vector3d force  = Eigen::Vector3d::Zero(); 
   Eigen::Vector3d moment = Eigen::Vector3d::Zero(); 
};
 
class ElementGaussPoints
{
public:
   void Begin(double length)
   {
      m_length = length;
      m_centres.clear();
      m_axes.clear();
   }
 
   void Add(const Eigen::Vector3d& centre, const Eigen::Vector3d& axis)
   {
      m_centres.push_back(centre);
      m_axes.push_back(axis);
   }
 
   void Build()
   {
      const double offset = m_length / (2.0 * std::sqrt(3.0));
      m_positions.resize(2 * m_centres.size());
      for (size_t i = 0; i < m_centres.size(); i++) {
         m_positions[2 * i]     = m_centres[i] - offset * m_axes[i];
         m_positions[2 * i + 1] = m_centres[i] + offset * m_axes[i];
      }
   }
 
   int NumElements() const { return static_cast<int>(m_centres.size()); }
 
   int NumPoints() const { return static_cast<int>(m_positions.size()); }
 
   const Eigen::Vector3d& Position(int point) const { return m_positions[point]; }
 
   void AppendPositions(std::vector<double>& points) const
   {
      points.reserve(points.size() + 3 * m_positions.size());
      for (const Eigen::Vector3d& position : m_positions)
         points.insert(points.end(), position.data(), position.data() + 3);
   }
 
   void QueryWater(const environment::wake::ElementFlowInteraction& flow, double T, int firstPoint)
   {
      m_water.resize(m_positions.size());
      for (int point = 0; point < NumPoints(); point++) {
         m_water[point].setZero();
         flow.WaterVelocity(T, firstPoint + point, m_positions[point].data(), m_water[point].data());
      }
   }
 
   EndWater WaterOf(int element) const
   {
      const Eigen::Vector3d& atMinus = m_water[2 * element];
      const Eigen::Vector3d& atPlus  = m_water[2 * element + 1];
      const double           ratio   = std::sqrt(3.0); // (L / 2) / a
      const Eigen::Vector3d  mean    = 0.5 * (atMinus + atPlus);
      const Eigen::Vector3d  halfRise = 0.5 * ratio * (atPlus - atMinus);
      const Eigen::Vector3d  minus   = mean - halfRise;
      const Eigen::Vector3d  plus    = mean + halfRise;
      return {{minus(0), minus(1), minus(2)}, {plus(0), plus(1), plus(2)}};
   }
 
private:
   double                       m_length = 0.0;
   std::vector<Eigen::Vector3d> m_centres;   
   std::vector<Eigen::Vector3d> m_axes;      
   std::vector<Eigen::Vector3d> m_positions; 
   std::vector<Eigen::Vector3d> m_water;     
};
 
template<class DragModel = environment::wet_cylinder::ReynoldsDrag>
inline GaussDragLoad GaussPointDragLoad(const environment::wet_cylinder::WetCylinder& cylinder, double f, const Eigen::Vector3d& V, const Eigen::Vector3d& W,
      const Eigen::Vector3d& axis, const EndWater& water, const DragModel& model = DragModel {})
{
   const auto parts = environment::wet_cylinder::GaussPointDragLoad(cylinder, model, f, {V(0), V(1), V(2)}, {W(0), W(1), W(2)}, {axis(0), axis(1), axis(2)}, water);
   return {Eigen::Vector3d(parts.force.data()), Eigen::Vector3d(parts.moment.data())};
}
} // namespace wet_cylinder_water