Merge pull request #11 from Fuad-HH/assemble_load

Assemble Load Vector. It just atomic adds.

Merging without review. Tests are passing.

Author: Fuad-HH

src/StiffnessMatrix.cpp

@@ -182,3 +182,24 @@ void StiffnessMatrix::printDenseMatrix() const {
     printf("\n");
   }
 }
+
+Kokkos::View<double *> assembleLoadVector(Kokkos::View<double *> elementLoad,
+                                          Mesh mesh) {
+  Kokkos::Profiling::pushRegion("Assemble Load Vector");
+
+  Kokkos::View<double *> loadVector("loadVector", mesh.GetNumVertices());
+  auto mesh_data = mesh.GetData();
+  int node_per_elem = mesh.GetMeshType();
+
+  Kokkos::parallel_for(
+      "Fill Load Vector", mesh.GetNumElements(),
+      KOKKOS_LAMBDA(const int elem_id) {
+        for (int node = 0; node < node_per_elem; ++node) {
+          int node_id = int(mesh_data(elem_id, node, 0));
+          Kokkos::atomic_add(&loadVector(node_id), elementLoad(elem_id));
+        }
+      });
+
+  Kokkos::Profiling::popRegion();
+  return loadVector;
+}

src/StiffnessMatrix.h

@@ -9,6 +9,11 @@
 
 #include "Mesh.h"
 
+// This function assembles the element load vector into the global load vector
+[[nodiscard]]
+Kokkos::View<double *> assembleLoadVector(Kokkos::View<double *> elementLoad,
+                                          Mesh mesh);
+
 // Contains the global indices of the element stiffness matrix
 // This is paired so that we can sort
 // data using Kokkos bitonic sorting algorithm

src/main.cpp

@@ -54,10 +54,21 @@ int main(int argc, char** argv) {
     stiffnessMatrix.sortDataByRowCol(element_stiffness);
     stiffnessMatrix.assemble(element_stiffness);
 
+    auto load_vector =
+        assembleLoadVector(el_stiffness_and_load.allElementLoadVector, mesh);
+
 #ifndef NDEBUG
     stiffnessMatrix.printStiffnessMatrix();
     printf("=>----------- Dense Matrix -----------<=\n");
     stiffnessMatrix.printDenseMatrix();
+
+    printf("=>----------- Load Vector -----------<=\n");
+    auto load_vector_host = Kokkos::create_mirror_view(load_vector);
+    Kokkos::deep_copy(load_vector_host, load_vector);
+    printf("Load vector:\n");
+    for (int i = 0; i < load_vector_host.size(); i++) {
+      printf("%f, \n", load_vector_host(i));
+    }
 #endif
   }
   Kokkos::finalize();

tests/CMakeLists.txt

@@ -46,3 +46,9 @@ add_executable(test_element_stiffness test_element_stiffness.cpp ${CMAKE_SOURCE_
 target_link_libraries(test_element_stiffness PRIVATE Catch2::Catch2WithMain Kokkos::kokkos)
 target_include_directories(test_element_stiffness PRIVATE ${CMAKE_SOURCE_DIR}/src)
 catch_discover_tests(test_element_stiffness WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR})
+
+# Test Load Vector
+add_executable(test_load_vector test_load_vector_assembly.cpp ${CMAKE_SOURCE_DIR}/src/Mesh.cpp ${CMAKE_SOURCE_DIR}/src/StiffnessMatrix.cpp ${CMAKE_SOURCE_DIR}/src/CalculateStiffnessMatrixAndLoadVector.cpp)
+target_link_libraries(test_load_vector PRIVATE Catch2::Catch2WithMain Kokkos::kokkos)
+target_include_directories(test_load_vector PRIVATE ${CMAKE_SOURCE_DIR}/src)
+catch_discover_tests(test_load_vector WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR})

tests/test_load_vector_assembly.cpp

@@ -0,0 +1,37 @@
+//
+// Created by Fuad Hasan on 4/17/25.
+//
+
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include "CalculateStiffnessMatrixAndLoadVector.hpp"
+#include "Mesh.h"
+#include "StiffnessMatrix.h"
+
+TEST_CASE("Test Load Assembly") {
+  Kokkos::initialize();
+  {
+    std::string mesh_filename = "assets/TriMesh.msh";
+    Mesh mesh(mesh_filename);
+
+    auto el_stiffness_load =
+        calculateAllElementStiffnessMatrixAndLoadVector(mesh, 1.0);
+    auto load_vector =
+        assembleLoadVector(el_stiffness_load.allElementLoadVector, mesh);
+
+    REQUIRE(load_vector.size() == mesh.GetNumVertices());
+
+    auto load_vector_host = Kokkos::create_mirror_view(load_vector);
+    Kokkos::deep_copy(load_vector_host, load_vector);
+
+    std::vector<double> expected_load_vector = {0.166667, 0.166667, 0.166667,
+                                                0.166667, 0.333333};
+
+    for (size_t i = 0; i < expected_load_vector.size(); ++i) {
+      REQUIRE(load_vector_host(i) ==
+              Catch::Approx(expected_load_vector[i]).epsilon(0.0001));
+    }
+  }
+  Kokkos::finalize();
+}

tests/test_matvec_mult.cpp

@@ -4,7 +4,6 @@
 
 #include <Kokkos_Core.hpp>
 #include <catch2/catch_test_macros.hpp>
-#include <iostream>
 
 #include "MatVecMult.h"
 #include "catch2/matchers/catch_matchers_floating_point.hpp"
@@ -85,7 +84,7 @@ TEST_CASE("Test matrix-vector multiplication") {
 
     // Expected result
     std::vector<double> expected_y = {5.0, 9.0, 18.0};
-    printf("y size = %d\n", y_h.size());
+    printf("y size = %zu\n", y_h.size());
     REQUIRE(y_h.size() == N);
     for (int i = 0; i < N; ++i) {
       printf("Expected y[%d] = %f, Actual y[%d] = %f\n", i, expected_y[i], i,