DiskDrag.h

#pragma once
 
#include "DiskBase.h"
 
#include <cmath>
 
namespace disk_drag
{
struct DragFactors
{
   double normal;     
   double tangential; 
};
 
inline DragFactors Factors(double rho, double cdNormal, double cdTangential, double diskD, double diskT, double cableD)
{
   return {DiskBase::DragForceNormal(rho, cdNormal, diskD, diskT, cableD, 1.0, 1.0),
         DiskBase::DragForceTangential(rho, cdTangential, diskD, diskT, cableD, 1.0, 1.0)};
}
 
inline void Decompose(const double n[3], const double v[3], double vNormal[3], double vTang[3])
{
   const double along = n[0] * v[0] + n[1] * v[1] + n[2] * v[2];
   for (int i = 0; i < 3; i++) {
      vNormal[i] = along * n[i];
      vTang[i]   = v[i] - vNormal[i];
   }
}
 
inline double Norm(const double w[3])
{
   return std::sqrt(w[0] * w[0] + w[1] * w[1] + w[2] * w[2]);
}
 
inline void DragForce(const double n[3], const double relVel[3], const DragFactors& k, double force[3])
{
   double vNormal[3], vTang[3];
   Decompose(n, relVel, vNormal, vTang);
   const double speedNormal = Norm(vNormal);
   const double speedTang   = Norm(vTang);
   for (int i = 0; i < 3; i++)
      force[i] = -k.normal * speedNormal * vNormal[i] - k.tangential * speedTang * vTang[i];
}
 
inline void DragVelocityJacobian(const double n[3], const double relVel[3], const DragFactors& k, double dForceDv[9])
{
   double vNormal[3], vTang[3];
   Decompose(n, relVel, vNormal, vTang);
   const double speedNormal = Norm(vNormal);
   const double speedTang   = Norm(vTang);
   for (int i = 0; i < 3; i++) {
      for (int j = 0; j < 3; j++) {
         const double identity   = (i == j) ? 1.0 : 0.0;
         const double projNormal = n[i] * n[j];
         const double projTang   = identity - projNormal;
         double       dNormal    = speedNormal * projNormal;
         if (speedNormal > 0.0) dNormal += vNormal[i] * vNormal[j] / speedNormal;
         double dTang = speedTang * projTang;
         if (speedTang > 0.0) dTang += vTang[i] * vTang[j] / speedTang;
         dForceDv[3 * i + j] = -k.normal * dNormal - k.tangential * dTang;
      }
   }
}
} // namespace disk_drag