wip (broken)

vlm/backends/cpu/include/vlm_backend_cpu.hpp

@@ -7,17 +7,7 @@ namespace vlm {
 
 class BackendCPU : public Backend {
     public:
-        std::vector<f32> lhs; // influence matrix (LHS)
-        SoA_3D_t<f32> rollup_vertices; // store displaced vertices at rollup
-        std::vector<f32> rhs; // inpenetrability condition (RHS)
-        SoA_3D_t<f32> local_velocities; // local velocities at collocation points
-        std::vector<i32> ipiv; // LU solver row pivots
-        std::vector<f32> gamma; // vortex strengths of the wing (current timestep) and wake (shed gammas)
-        std::vector<f32> gamma_prev; // vortex strengths of the wing (previous timestep)
-        std::vector<f32> delta_gamma; // chordwise gamma difference
-        std::vector<f32> trefftz_buffer;
-
-        BackendCPU(Mesh& mesh);
+        BackendCPU(const MeshGeom* mesh_geom, int timesteps);
         ~BackendCPU();
         void reset() override;
         void lhs_assemble() override;

vlm/include/vlm_allocator.hpp

@@ -0,0 +1,13 @@
+#pragma once
+
+namespace vlm {
+
+struct Allocator {
+    Allocator() = default;
+    virtual ~Allocator() = default;
+
+    virtual void* malloc(unsigned long long size) const = 0;
+    virtual void free(void* ptr) const = 0;
+};
+
+} // namespace vlm

vlm/include/vlm_backend.hpp

@@ -3,15 +3,21 @@
 #include "vlm_fwd.hpp"
 #include "vlm_types.hpp"
 #include "vlm_mesh.hpp"
+#include "vlm_data.hpp"
 
 namespace vlm {
 
 class Backend {
     public:
-        Mesh& mesh;
+        Allocator* h_allocator = nullptr;
+        Allocator* d_allocator = nullptr;
+        MeshProxy* h_mesh; 
+        MeshProxy* d_mesh; 
+        Data* h_data; // host (CPU) data
+        Data* d_data; // device (accelerator) data
         f32 sigma_vatistas = 0.0f;
 
-        Backend(Mesh& mesh) : mesh(mesh) {};
+        Backend(const MeshGeom* mesh_geom, int timesteps) {};
         virtual void reset() = 0;
         virtual void lhs_assemble() = 0;
         virtual void compute_rhs() = 0;
@@ -33,7 +39,7 @@ class Backend {
         virtual ~Backend() = default;
 };
 
-std::unique_ptr<Backend> create_backend(const std::string& backend_name, Mesh& mesh);
+std::unique_ptr<Backend> create_backend(const std::string& backend_name, const MeshGeom* mesh, int timesteps);
 std::vector<std::string> get_available_backends();
 
 } // namespace vlm

vlm/include/vlm_data.hpp

@@ -6,6 +6,22 @@
 
 namespace vlm {
 
+struct Data {
+    f32* lhs = nullptr;
+    f32* rhs = nullptr;
+    f32* gamma = nullptr;
+    f32* delta_gamma = nullptr;
+    f32* rollup_vertices = nullptr;
+    f32* local_velocities = nullptr;
+    f32* trefftz_buffer = nullptr;
+
+    i32* solver_info = nullptr;
+    i32* solver_ipiv = nullptr;
+    f32* solver_buffer = nullptr;
+};
+
+void data_alloc(const Allocator* allocator, Data* data, u64 n, u64 m);
+
 // Flow characteristics
 class FlowData {
     public:

vlm/include/vlm_fwd.hpp

@@ -3,6 +3,10 @@
 namespace vlm {
 
 class FlowData;
+struct MeshView;
+struct MeshGeom;
+struct Data;
+struct Allocator;
 class Mesh;
 class Backend;
 template<typename Interpolator> class Database2D;

vlm/include/vlm_mesh.hpp

@@ -1,73 +1,55 @@
 #pragma once
 
 #include "vlm_types.hpp"
+#include "vlm_allocator.hpp"
+
 #include "tinyfwd.hpp"
 #include "tinyview.hpp"
 
 namespace vlm {
 
-// class SMesh {
-//     public:
-
-//     tiny::View<f32, 3> vertices_wing();
-//     tiny::View<f32, 3> vertices_wake();
-//     tiny::View<f32, 3> normals();
-//     tiny::View<f32, 3> colloc();
-//     tiny::View<f32, 2> areas();
-
-//     virtual void alloc(u64 nc, u64 ns, u64 nw) = 0;
-//     virtual void dealloc() = 0;
-
-//     private:
-//     f32* _vertices = nullptr;
-//     f32* _normals = nullptr;
-//     f32* _colloc = nullptr;
-//     f32* _area = nullptr;
-
-//     u64 _nc = 0; // number of wing panels chordwise
-//     u64 _ns = 0; // number of wing panels spanwise
-//     u64 _nw = 0; // number of wake panels chordwise
-//     u64 _naw = 0; // number of active wake panels chordwise
-
-//     f32 _s_ref = 0.0f; // reference area
-//     f32 _c_ref = 0.0f; // reference chord
-//     linalg::alias::float4x4 frame = linalg::identity; // mesh orientation and position in inertial frame
-// };
-
-
-/*
-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;
+// C interface
+struct MeshProxy {
+    f32* vertices = nullptr; // (nc+nw+1)*(ns+1)*3
+    // TODO: evaluate if we really need those values or we can just compute on the fly in the kernels
+    f32* normals = nullptr; // nc*ns*3
+    f32* colloc = nullptr; // nc*ns*3
+    f32* area = nullptr; // nc*ns
+
+    u64 nc = 0; // number of wing panels chordwise
+    u64 ns = 0; // number of wing panels spanwise
+    u64 nw = 0; // number of wake panels chordwise
+    u64 nwa = 0; // number of wake panels active chordwise
+
+    f32 s_ref = 0.0f; // reference area
+    f32 c_ref = 0.0f; // reference chord
+    f32 frame[16] = {}; // 4x4 homogeneous matrix in col major order
+
+    char* name = nullptr; // mesh name
+    bool lifting = true;
 };
 
-class VLMSurface : public StructuredSurfaceMesh {
-    public:
-    SoA_3D_t<f32> colloc; // collocation points
-
-    const f32 s_ref; // reference area
-    const f32 c_ref; // reference chord
+// Contains only the raw geometry 
+struct MeshGeom {
+    f32* vertices = nullptr;
+    u64 nc = 0; // number of wing panels chordwise
+    u64 ns = 0;
+};
 
-    VLMSurface(const u64 ni, const u64 nj);
-    ~VLMSurface() = default;
+struct System {
+    MeshProxy* meshes = nullptr; // meshes
+    // Data data; // unified (linearized) data for the whole system
+    u32 n_meshes = 0;
 };
 
-class LiftingBody {
-    public:
-    VLMSurface body;
-    VLMSurface wake;
+void mesh_read_file(const std::string& filename, MeshGeom* mesh_geom);
+void mesh_alloc(const Allocator* allocator, MeshProxy* mesh, u64 nc, u64 ns, u64 nw, u64 nwa);
+void mesh_free(const Allocator* allocator, MeshProxy* mesh);
+
+struct AllocatorCPU : Allocator {
+    void* malloc(u64 size) const override { return std::malloc(size); }
+    void free(void* ptr) const override { std::free(ptr); }
 };
-*/
 
 
 // === STRUCTURED MESH ===

vlm/include/vlm_types.hpp

@@ -67,38 +67,16 @@ struct SoA_3D_t {
     }
 };
 
-template<typename T, size_t N>
-class StackArray {
-    T data[N];
-};
-
 template<typename T>
-class HeapArray {
-    T* first;
-    T* last;
-    T* end;
-};
-
-template<typename T, size_t N>
-class View {
-    private:
-        const T* data;
-        const StackArray<T, N> dims;
-        const StackArray<T, N> strides;
-    public:
-        T* operator()() {
-            return data;
-        }
-        template<typename... Idx>
-        const T& operator()(Idx... idx) const {
-            static_assert(sizeof...(idx) == N, "The number of indices must match the dimension N.");
-            const size_t indices[N] = {idx...};
-            size_t index = 0;
-            for (size_t i = 0; i < N; ++i) {
-                index += indices[i] * strides.data[i];
-            }
-            return data[index];
-        }
+struct Status {
+    bool ok = false; // true if no error
+    union {
+        struct err_t {
+            int args = 0; // bitfield for position of wrong arguments
+            char* msg = nullptr; // error message
+        } err;
+        T data; // data to return
+    };
 };
 
 } // namespace vlm

vlm/src/vlm_mesh.cpp

@@ -14,6 +14,100 @@ using namespace vlm;
 
 constexpr f32 EPS = std::numeric_limits<f32>::epsilon();
 
+void vlm::mesh_alloc(const Allocator* allocator, MeshProxy* mesh, u64 nc, u64 ns, u64 nw, u64 nwa) {
+    mesh->vertices = (f32*)allocator->malloc((nc+nw+1)*(ns+1)*3*sizeof(f32));
+    mesh->normals = (f32*)allocator->malloc(nc*ns*3*sizeof(f32));
+    mesh->colloc = (f32*)allocator->malloc(nc*ns*3*sizeof(f32));
+    mesh->area = (f32*)allocator->malloc(nc*ns*sizeof(f32));
+    mesh->ns = nc;
+    mesh->ns = ns;
+    mesh->nw = nw;
+    mesh->nwa = nwa;
+}
+
+void vlm::mesh_free(const Allocator* allocator, MeshProxy* mesh) {
+    allocator->free(mesh->vertices);
+    allocator->free(mesh->normals);
+    allocator->free(mesh->colloc);
+    allocator->free(mesh->area);
+}
+
+// Reads the Plot3D Structured file, allocates vertex buffer and fills it with the file data
+void mesh_io_read_file_plot3d(std::ifstream& f, MeshGeom* mesh_geom) {
+    std::cout << "reading plot3d mesh" << std::endl;
+    u64 ni, nj, nk, blocks, nb_panels;
+    f32 x, y, z;
+    f >> blocks;
+    if (blocks != 1) {
+        throw std::runtime_error("Only single block plot3d mesh is supported");
+    }
+    f >> ni >> nj >> nk;
+    if (nk != 1) {
+        throw std::runtime_error("Only 2D plot3d mesh is supported");
+    }
+
+    mesh_geom->ns = nj - 1;
+    mesh_geom->nc = ni - 1;
+    mesh_geom->vertices = new f32[ni*nj*3];
+    nb_panels = mesh_geom->ns * mesh_geom->nc;
+
+    std::cout << "number of panels: " << nb_panels << std::endl;
+    std::cout << "ns: " << mesh_geom->ns << std::endl;
+    std::cout << "nc: " << mesh_geom->nc << std::endl;
+
+    for (u64 j = 0; j < nj; j++) {
+        for (u64 i = 0; i < ni; i++) {
+            f >> x;
+            mesh_geom->vertices[nb_panels * 0 + nj*i + j] = x;
+        }
+    }
+    for (u64 j = 0; j < nj; j++) {
+        for (u64 i = 0; i < ni; i++) {
+            f >> y;
+            mesh_geom->vertices[nb_panels * 1 + nj*i + j] = y;
+        }
+    }
+    for (u64 j = 0; j < nj; j++) {
+        for (u64 i = 0; i < ni; i++) {
+            f >> z;
+            mesh_geom->vertices[nb_panels * 2 + nj*i + j] = z;
+        }
+    }
+
+    // TODO: this validation is not robust enough when eventually we will be using multiple bodies
+    // Validation
+    // const f32 eps = std::numeric_limits<f32>::epsilon();
+    // if (std::abs(m.v.x[0]) != eps || std::abs(m.v.y[0]) != eps) {
+    //     throw std::runtime_error("First vertex of plot3d mesh must be at origin");
+    // }
+    const f32 first_y = mesh_geom->vertices[nb_panels * 1 + 0];
+    for (u64 i = 1; i < ni; i++) {
+        if ( mesh_geom->vertices[nb_panels * 1 + i*nj] != first_y) {
+            throw std::runtime_error("Mesh vertices should be ordered in chordwise direction");
+        }
+    }
+}
+
+void mesh_io_read_file(const std::string& filename, MeshGeom* mesh_geom) {
+    const std::filesystem::path path(filename);
+    if (!std::filesystem::exists(path)) {
+        throw std::runtime_error("Mesh file not found"); // TODO: consider not using exceptions anymore
+    }
+
+    std::ifstream f(path);
+    if (f.is_open()) {
+        if (path.extension() == ".x") {
+            mesh_io_read_file_plot3d(f, mesh_geom);
+        } else {
+            throw std::runtime_error("Only structured gridpro mesh format is supported");
+        }
+        f.close();
+    } else {
+        throw std::runtime_error("Failed to open mesh file");
+    }
+}
+
+// TODO: deprecate everything below this line 
 Mesh::Mesh(const tiny::Config& cfg) :
     Mesh(
         cfg().section("files").get<std::string>("mesh"),