diff --git a/headeronly/linalg.h b/headeronly/linalg.h index 37c449c..1626393 100644 --- a/headeronly/linalg.h +++ b/headeronly/linalg.h @@ -246,23 +246,6 @@ namespace linalg }; // Small, fixed-size matrix type, consisting of exactly M rows and N columns of type T, stored in column-major order. - template struct mat - { - typedef vec V; - V x; - constexpr mat() : x() {} - constexpr mat(const V & x_) : x(x_) {} - constexpr explicit mat(const T & s) : x(s) {} - constexpr explicit mat(const T * p) : x(p+M*0) {} - template - constexpr explicit mat(const mat & m) : mat(V(m.x)) {} - constexpr vec row(int i) const { return {x[i]}; } - constexpr const V & operator[] (int j) const { return x; } - constexpr V & operator[] (int j) { return x; } - - template> constexpr mat(const U & u) : mat(converter{}(u)) {} - template> constexpr operator U () const { return converter{}(*this); } - }; template struct mat { typedef vec V; @@ -281,7 +264,7 @@ namespace linalg template> constexpr operator U () const { return converter{}(*this); } }; template struct mat - { + { typedef vec V; V x,y,z; constexpr mat() : x(), y(), z() {} @@ -292,8 +275,8 @@ namespace linalg template constexpr explicit mat(const mat & m) : mat(V(m.x), V(m.y), V(m.z)) {} constexpr vec row(int i) const { return {x[i], y[i], z[i]}; } - constexpr const V & operator[] (int j) const { return j==0?x:j==1?y:z; } - constexpr V & operator[] (int j) { return j==0?x:j==1?y:z; } + constexpr const V & col (int j) const { return j==0?x:j==1?y:z; } + constexpr V & col (int j) { return j==0?x:j==1?y:z; } template> constexpr mat(const U & u) : mat(converter{}(u)) {} template> constexpr operator U () const { return converter{}(*this); } @@ -310,8 +293,8 @@ namespace linalg template // typecast matrix ? constexpr explicit mat(const mat & m) : mat(V(m.x), V(m.y), V(m.z), V(m.w)) {} constexpr vec row(int i) const { return {x[i], y[i], z[i], w[i]}; } - constexpr const V & operator[] (int j) const { return j==0?x:j==1?y:j==2?z:w; } - constexpr V & operator[] (int j) { return j==0?x:j==1?y:j==2?z:w; } + constexpr const V & col (int j) const { return j==0?x:j==1?y:j==2?z:w; } + constexpr V & col (int j) { return j==0?x:j==1?y:j==2?z:w; } template> constexpr mat(const U & u) : mat(converter{}(u)) {} template> constexpr operator U () const { return converter{}(*this); } diff --git a/tests/test_uvlm_theodorsen.cpp b/tests/test_uvlm_theodorsen.cpp index a965315..c7d3dc6 100644 --- a/tests/test_uvlm_theodorsen.cpp +++ b/tests/test_uvlm_theodorsen.cpp @@ -30,7 +30,7 @@ int main() { }); }; - // For the moment take AoA=0 and U_inf=1 + // 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, @@ -43,6 +43,11 @@ int main() { return linalg::mul(displacement_freestream(t), displacement_wing(t)); }; + auto relative_displacement = [&](f32 t0, f32 t1) -> linalg::alias::float4x4 { + return linalg::mul(displacement(t1), linalg::inverse(displacement(t0))); + }; + + // Magnitude of the velocity of a vertex at time t calculated by finite difference 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); }; @@ -52,24 +57,18 @@ int main() { 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 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 vec_dt; // pre-calculated timesteps + std::vector vec_t; // pre-calculated timesteps // Pre-calculate timesteps to determine wake size + vec_t.push_back(0.0f); 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; + vec_t.push_back(t); } // TODO: think about splitting mesh body and wake mesh - // solver.mesh->alloc(vec_dt.size()+1); // +1 for the initial pose + solver.mesh->resize_wake(vec_t.size()-1); // +1 for the initial pose // f32 t = 0.0f; // for (u64 i = 0; i < vec_dt.size(); i++) { diff --git a/vlm/include/vlm_mesh.hpp b/vlm/include/vlm_mesh.hpp index 8743664..7bf0eb9 100644 --- a/vlm/include/vlm_mesh.hpp +++ b/vlm/include/vlm_mesh.hpp @@ -71,13 +71,15 @@ class Mesh { u64 nc = 1; // number of panels chordwise u64 ns = 1; // number of panels spanwise u64 nw; // number of wake panels chordwise (max capacity) - u64 current_nw = 0; // current number of built wake panels + u64 current_nw = 0; // current number of built wake panels (temporary) // Analytical quanities when provided // --------------------- f32 s_ref = 0.0f; // reference area (of wing) f32 c_ref = 0.0f; // reference chord (of wing) - linalg::alias::float3 ref_pt = {0.25f, 0.0f, 0.0f}; + + linalg::alias::float4x4 frame = linalg::identity; // transformation matrix + linalg::alias::float3 ref_pt = {0.25f, 0.0f, 0.0f}; // TODO: deprecate void update_wake(const linalg::alias::float3& u_inf); // update wake vertices void correction_high_aoa(f32 alpha_rad); // correct collocation points for high aoa @@ -85,9 +87,10 @@ class Mesh { void compute_metrics_wing(); // fill colloc, normal, area void compute_metrics_wake(); void compute_metrics_i(u64 i); - void transform(const linalg::alias::float4x4& transform); void shed_wake(); void move(const linalg::alias::float4x4& transform); + void resize_wake(const u64 nw); + u64 nb_panels_wing() const; u64 nb_panels_total() const; u64 nb_vertices_wing() const; diff --git a/vlm/src/vlm_mesh.cpp b/vlm/src/vlm_mesh.cpp index c690f0c..b124066 100644 --- a/vlm/src/vlm_mesh.cpp +++ b/vlm/src/vlm_mesh.cpp @@ -368,10 +368,12 @@ void Mesh::io_read(const std::string& filename) { throw std::runtime_error("Failed to open mesh file"); } } - -void Mesh::transform(const linalg::alias::float4x4& transform) { - // like a functional map but over multidimensional elements + +void Mesh::move(const linalg::alias::float4x4& transform) { for (u64 i = 0; i < nb_vertices_wing(); i++) { + linalg::alias::float4 global_vertex{v.x[i], v.y[i], v.z[i], 1.f}; + linalg::alias::float4 local_vertex = frame.col(3); + const linalg::alias::float4 transformed_pt = linalg::mul(transform, linalg::alias::float4{v.x[i], v.y[i], v.z[i], 1.f}); v.x[i] = transformed_pt.x; v.y[i] = transformed_pt.y; @@ -379,14 +381,19 @@ void Mesh::transform(const linalg::alias::float4x4& transform) { } } +void Mesh::resize_wake(const u64 nw_) { + nw = nw_; + alloc(); +} + // Duplicate the trailing edge vertices and release them in the buffer void Mesh::shed_wake() { - const u64 src_begin = nb_vertices_wing() - ns - 1; + assert(current_nw < nw); // check if we have capacity + const u64 src_begin = nb_vertices_wing() - (ns + 1); const u64 src_end = nb_vertices_wing(); - const u64 dst_begin = (nc + 1 + current_nw) * (ns + 1); + const u64 dst_begin = (nc + nw - current_nw) * (ns + 1); std::copy(v.x.data() + src_begin, v.x.data() + src_end, v.x.data() + dst_begin); std::copy(v.y.data() + src_begin, v.y.data() + src_end, v.y.data() + dst_begin); std::copy(v.z.data() + src_begin, v.z.data() + src_end, v.z.data() + dst_begin); - current_nw++; } diff --git a/xmake.lua b/xmake.lua index 9b6f440..c9f5cc6 100644 --- a/xmake.lua +++ b/xmake.lua @@ -7,7 +7,7 @@ add_rules("mode.debug", "mode.release", "mode.releasedbg", "mode.asan") -- set_toolchains("cuda") -- set_toolset("cxx", "clang") --- set_policy("build.sanitizer.address", true) +set_policy("build.sanitizer.address", true) set_policy("build.warning", true) set_policy("build.cuda.devlink", true) -- magic set_policy("run.autobuild", true)