some progress

headeronly/linalg.h

@@ -322,7 +322,7 @@ namespace linalg
     template<class T> struct converter<mat<T,1,1>, identity_t> { constexpr mat<T,1,1> operator() (identity_t) const { return {vec<T,1>{1}}; } };
     template<class T> struct converter<mat<T,2,2>, identity_t> { constexpr mat<T,2,2> operator() (identity_t) const { return {{1,0},{0,1}}; } };
     template<class T> struct converter<mat<T,3,3>, identity_t> { constexpr mat<T,3,3> operator() (identity_t) const { return {{1,0,0},{0,1,0},{0,0,1}}; } };
-    template<class T> struct converter<mat<T,4,4>, identity_t> { constexpr mat<T,4,4> operator() (identity_t) const { return {{1,0,0,0},{0,1,0,0},{0,0,1,0},{0,0,0,1}}; } };
+    template<class T> struct converter<mat<T,4,4>, identity_t> { constexpr mat<T,4,4> operator() (identity_t) const { return {{1,0,0,0},{0,1,0,0},{0,0,1,0},{0,0,0,1}}; } };    
     constexpr identity_t identity {1};
 
     // Produce a scalar by applying f(A,B) -> A to adjacent pairs of elements from a vec/mat in left-to-right/column-major order (matching the associativity of arithmetic and logical operators)
@@ -465,6 +465,7 @@ namespace linalg
     template<class T, int M> constexpr T dot      (const vec<T,M> & a, const vec<T,M> & b)      { return sum(a*b); }
     template<class T, int M> constexpr T length2  (const vec<T,M> & a)                          { return dot(a,a); }
     template<class T, int M> T           length   (const vec<T,M> & a)                          { return std::sqrt(length2(a)); }
+    template<class T, int M> T           norm     (const vec<T,M> & a)                          { return length(a); }
     template<class T, int M> vec<T,M>    normalize(const vec<T,M> & a)                          { return a / length(a); }
     template<class T, int M> constexpr T distance2(const vec<T,M> & a, const vec<T,M> & b)      { return length2(b-a); }
     template<class T, int M> T           distance (const vec<T,M> & a, const vec<T,M> & b)      { return length(b-a); }
@@ -549,6 +550,7 @@ namespace linalg
     template<class T> vec<T,4>   rotation_quat     (const mat<T,3,3> & m);
     template<class T> mat<T,4,4> translation_matrix(const vec<T,3> & translation)           { return {{1,0,0,0},{0,1,0,0},{0,0,1,0},{translation,1}}; }
     template<class T> mat<T,4,4> rotation_matrix   (const vec<T,4> & rotation)              { return {{qxdir(rotation),0}, {qydir(rotation),0}, {qzdir(rotation),0}, {0,0,0,1}}; }
+    template<class T> mat<T,4,4> rotation_matrix   (const vec<T,3> & axis, const vec<T,3> & p, T angle) {}
     template<class T> mat<T,4,4> scaling_matrix    (const vec<T,3> & scaling)               { return {{scaling.x,0,0,0}, {0,scaling.y,0,0}, {0,0,scaling.z,0}, {0,0,0,1}}; }
     template<class T> mat<T,4,4> pose_matrix       (const vec<T,4> & q, const vec<T,3> & p) { return {{qxdir(q),0}, {qydir(q),0}, {qzdir(q),0}, {p,1}}; }
     template<class T> mat<T,4,4> lookat_matrix     (const vec<T,3> & eye, const vec<T,3> & center, const vec<T,3> & view_y_dir, fwd_axis fwd = neg_z);

tests/test_uvlm_theodorsen.cpp

@@ -10,35 +10,75 @@
 #include "vlm_utils.hpp"
 
 using namespace vlm;
+using namespace linalg::ostream_overloads;
 
-int main(int argc, char **argv) {
+int main() {
     const std::vector<std::string> meshes = {"../../../../mesh/infinite_rectangular_5x200.x"};
     const std::vector<std::string> backends = get_available_backends();
 
     auto solvers = tiny::make_combination(meshes, backends);
 
+    // Define simulation length
+    const f32 t_final = 10.0f;
+
+    // Define wing motion
+    auto displacement_wing = [&](f32 t) -> linalg::alias::float4x4 {
+        return linalg::translation_matrix(linalg::alias::float3{
+            0.0f,
+            0.0f,
+            std::sin(0.2f * t)
+        });
+    };
+
+    // For the moment take AoA=0 and U_inf=1 
+    auto displacement_freestream = [&](f32 t) -> linalg::alias::float4x4 {
+        return linalg::translation_matrix(linalg::alias::float3{
+            -1.0f*t,
+            0.0f,
+            0.0f
+        });
+    };
+
+    auto displacement = [&](f32 t) -> linalg::alias::float4x4 {
+        return linalg::mul(displacement_freestream(t), displacement_wing(t));
+    };
+
+    auto absolute_velocity = [&](f32 t, const linalg::alias::float4& vertex) -> f32 {
+        return linalg::length((linalg::mul(displacement(t+EPS_sqrt_f), vertex)-linalg::mul(displacement(t), vertex))/EPS_sqrt_f);
+    };
+
     for (const auto& [mesh_name, backend_name] : solvers) {
+        UVLM solver{};
+        solver.mesh = create_mesh(mesh_name);
+        solver.backend = create_backend(backend_name, *solver.mesh);
+
+        // Copy initial relative pose to the body frame
+        // SoA_3D_t<f32> initial_pose;
+        // initial_pose.resize(solver.mesh->nb_vertices_wing());
+        // initial_pose.x = solver.mesh->v.x;
+        // initial_pose.y = solver.mesh->v.y;
+        // initial_pose.z = solver.mesh->v.z;
+
+        std::vector<f32> vec_dt; // pre-calculated timesteps
+
+        // Pre-calculate timesteps to determine wake size
+        for (f32 t = 0.0f; t < t_final;) {
+            const f32 dt = (0.5f * (solver.mesh->chord_length(0) + solver.mesh->chord_length(1))) / absolute_velocity(t, {0.f, 0.f, 0.f, 1.f});
+            vec_dt.push_back(dt);
+            t += dt;
+        }
 
-        auto translation = [&](f32 t) -> linalg::alias::float4x4 {
-            return linalg::translation_matrix(linalg::alias::float3{
-                -1.0f*t,
-                0.0f,
-                std::sin(0.2f * t)
-            });
-        };
-
-        // auto rotation = [&](f32 t) -> linalg::alias::float4x4 {
-        //     return linalg::translation_matrix(linalg::alias::float3{
-        //         -1.0f*t,
-        //         0.0f,
-        //         std::sin(0.2f * t)
-        //     });
-        // };
-
-        auto wing_pose = [&](f32 t) -> linalg::alias::float4x4 {
-            return translation(t);
-        };
+        // TODO: think about splitting mesh body and wake mesh
+        // solver.mesh->alloc(vec_dt.size()+1); // +1 for the initial pose
 
+        // f32 t = 0.0f;
+        // for (u64 i = 0; i < vec_dt.size(); i++) {
+        //     FlowData flow{};
+        //     solver.solve(flow);
+        //     solver.mesh->shed_wake();
+        //     solver.mesh->move(linalg::mul(displacement(t+vec_dt[i]), linalg::inverse(displacement(t))));
+        //     t += vec_dt[i];
+        // }
     }
 
     return 0;

vlm/include/vlm.hpp

@@ -47,6 +47,7 @@ class LinearVLM final: public Solver {
 
 class UVLM final: public Solver {
     public:
+
     UVLM(const tiny::Config& cfg): Solver(cfg) {}
     UVLM() = default;
     ~UVLM() = default;

vlm/include/vlm_mesh.hpp

@@ -4,6 +4,40 @@
 #include "tinyfwd.hpp"
 
 namespace vlm {
+/*
+class StructuredSurfaceMesh {
+    public:
+    
+    SoA_3D_t<f32> v;
+    SoA_3D_t<f32> normal;
+    std::vector<f32> area = {}; // size ncw*ns
+
+    StructuredSurfaceMesh(const u64 ni, const u64 nj);
+    ~StructuredSurfaceMesh() = default;
+
+    private:
+    const u64 ni;
+    const u64 nj;
+};
+
+class VLMSurface : public StructuredSurfaceMesh {
+    public:
+    SoA_3D_t<f32> colloc; // collocation points
+
+    const f32 s_ref; // reference area
+    const f32 c_ref; // reference chord
+
+    VLMSurface(const u64 ni, const u64 nj);
+    ~VLMSurface() = default;
+};
+
+class LiftingBody {
+    public:
+    VLMSurface body;
+    VLMSurface wake;
+};
+*/
+
 
 // === STRUCTURED MESH ===
 // nc : number of panels chordwise
@@ -53,6 +87,7 @@ class Mesh {
     void compute_metrics_i(u64 i);
     void transform(const linalg::alias::float4x4& transform);
     void shed_wake();
+    void move(const linalg::alias::float4x4& transform);
     u64 nb_panels_wing() const;
     u64 nb_panels_total() const;
     u64 nb_vertices_wing() const;
@@ -61,7 +96,7 @@ class Mesh {
     f32 panels_area_xy(const u64 i, const u64 j, const u64 m, const u64 n) const;
     f32 panel_width_y(const u64 i, const u64 j) const;
     f32 strip_width(const u64 j) const;
-    f32 chord_length(const u64 j) const;
+    f32 chord_length(const u64 j) const; // vertex idx
     f32 chord_mean(const u64 j, const u64 n) const;
 
     // i panel vertices coordinates

vlm/include/vlm_types.hpp

@@ -6,7 +6,6 @@
 #include <array>
 #include <string>
 #include <cmath>
-#include <limits>
 #include <memory>
 
 #include "linalg.h"
@@ -26,10 +25,14 @@ using i64 = std::int64_t;
 using f32 = float;
 using f64 = double;
 
-constexpr f64 PI_d = 3.14159265358979;
-constexpr f32 PI_f = 3.141593f;
-constexpr f64 EPS_d = std::numeric_limits<f64>::epsilon();
-constexpr f32 EPS_f = std::numeric_limits<f32>::epsilon();
+static const f64 PI_d = 3.14159265358979;
+static const f32 PI_f = 3.141593f;
+static const f64 EPS_d = 2.22045e-16;
+static const f32 EPS_f = 1.19209e-07f;
+static const f64 EPS_sqrt_d = std::sqrt(EPS_d);
+static const f32 EPS_sqrt_f = std::sqrt(EPS_f);
+static const f64 EPS_cbrt_d = std::cbrt(EPS_d);
+static const f32 EPS_cbrt_f = std::cbrt(EPS_f);
 
 // Optimized n power function that uses template folding optimisations
 // to generate a series of mul instructions
@@ -64,4 +67,12 @@ struct SoA_3D_t {
     }
 };
 
+template<typename T>
+struct SoA3D_t {
+    T* x = nullptr;
+    T* y = nullptr;
+    T* z = nullptr;
+    u64 size = 0;
+};
+
 } // namespace vlm