Merge pull request #7 from Fuad-HH/stiffnessMatrixLoad

Element Stiffness Matrix Calculation

Author: Fuad-HH

.gitignore

@@ -8,6 +8,7 @@ build-*
 *.o
 *.obj
 
+*.sh
 # Precompiled Headers
 *.gch
 *.pch

openmp-scorec-config.sh

@@ -3,7 +3,7 @@ cmake -B build -S . \
   -DCMAKE_INSTALL_PREFIX=build/install \
   -DCMAKE_CXX_COMPILER=`which mpicxx` \
   -DCMAKE_C_COMPILER=`which mpicc` \
-  -DKokkos_ROOT=../assignment2_deps/kokkos/build/install/ \
+  -DKokkos_ROOT=/lore/paudea/build/ADA89/kokkos/install/ \
   -DAssignment_ENABLE_CUDA=OFF \
   -DCatch2_ROOT=/lore/hasanm4/practice_projects/advanced_comp/assignment2_deps/Catch2/build/install/
 

src/CMakeLists.txt

@@ -3,12 +3,17 @@
 set(Assignment2_SOURCES
     main.cpp
     Mesh.cpp
+	CalculateStiffnessMatrixAndLoadVector.cpp	
     StiffnessMatrix.cpp
     MatVecMult.cpp
 )
 
 set (Assignment2_HEADERS
     Mesh.h
+	Element.hpp
+	TriElement.hpp
+	QuadElement.hpp
+	CalculateStiffnessMatrixAndLoadVector.hpp	
     StriffnessMatrix.h
     MatVecMult.h
 )

src/CalculateStiffnessMatrixAndLoadVector.cpp

@@ -0,0 +1,46 @@
+#include "CalculateStiffnessMatrixAndLoadVector.hpp" 
+
+
+Results calculateAllElementStiffnessMatrixAndLoadVector(const Mesh& mesh, double k) {
+
+    int numElements = mesh.GetNumElements();
+	    
+	int numNodes = mesh.GetNumNodesPerElement();
+
+    int sizePerElement = numNodes * numNodes; 
+
+	View1D allElementStiffnessMatrix("stores all element stiffness matrix", numElements * sizePerElement);
+	View1D allElementLoadVector("stores all element load vector", numElements * numNodes);
+    
+	Kokkos::parallel_for("CalculateStiffness", numElements,
+        KOKKOS_LAMBDA(const int elemIdx) {
+
+		double stiffnessMatrixPerElement[MAX_STIFFNESS_MATRIX_SIZE] = {};
+		double loadVectorPerElement[MAX_LOAD_VECTOR_SIZE] = {};
+		if (numNodes == 3) { // Triangle element
+                TriElement triElem(mesh, elemIdx);
+                triElem.setMaterialProperty(k);
+                triElem.computeElementStiffnessMatrix(stiffnessMatrixPerElement);
+				triElem.computeElementLoadVector(loadVectorPerElement);
+        } else { // Quad element
+                QuadElement quadElem(mesh, elemIdx);
+                quadElem.setMaterialProperty(k);
+                quadElem.computeElementStiffnessMatrix(stiffnessMatrixPerElement);
+				quadElem.computeElementLoadVector(loadVectorPerElement);
+        }
+			
+		int base_stiffness_idx = elemIdx * sizePerElement;
+		for (int i = 0; i < sizePerElement; ++i){
+			allElementStiffnessMatrix(base_stiffness_idx + i) = stiffnessMatrixPerElement[i]; 	
+		}
+		
+		int base_load_idx = elemIdx * numNodes; 
+
+		for (int i = 0; i < numNodes; ++i){
+			allElementLoadVector(base_load_idx + i) = loadVectorPerElement[i];	
+		}
+
+    });
+
+	return Results{allElementStiffnessMatrix, allElementLoadVector};
+}

src/CalculateStiffnessMatrixAndLoadVector.hpp

@@ -0,0 +1,20 @@
+#ifndef COMPUTING_AT_SCALE_ASSIGNMENT_CALCULATE_STIFFNESS_MATRIX_AND_LOAD_VECTOR_HPP 
+#define COMPUTING_AT_SCALE_ASSIGNMENT_CALCULATE_STIFFNESS_MATRIX_AND_LOAD_VECTOR_HPP 
+
+#include "TriElement.hpp"
+#include "QuadElement.hpp"
+
+constexpr static int MAX_STIFFNESS_MATRIX_SIZE = 16;
+constexpr static int MAX_LOAD_VECTOR_SIZE = 4;
+
+
+struct Results {
+	View1D allElementStiffnessMatrix;
+	View1D allElementLoadVector;
+
+};
+
+
+Results calculateAllElementStiffnessMatrixAndLoadVector(const Mesh& mesh, double k = 1.0);
+
+#endif 

src/Element.hpp

@@ -0,0 +1,53 @@
+#ifndef COMPUTING_AT_SCALE_ASSIGNMENT_ELEMENT_HPP
+#define COMPUTING_AT_SCALE_ASSIGNMENT_ELEMENT_HPP
+
+#include <Kokkos_Core.hpp>
+#include "Mesh.h"
+
+using View1D = Kokkos::View<double*>;
+
+class Element {
+protected:
+    const Mesh& mesh_;
+    int elemIdx_;
+	double k_;
+
+public:
+	KOKKOS_INLINE_FUNCTION
+    Element(const Mesh& mesh, int elemIdx, double k = 1.0) : mesh_(mesh), elemIdx_(elemIdx), k_(k) {}
+   
+	KOKKOS_INLINE_FUNCTION
+	virtual ~Element() {}
+
+    // Get number of nodes per element
+    KOKKOS_INLINE_FUNCTION
+    virtual int getNumNodes() const = 0;
+    
+	// Compute local basis function
+    KOKKOS_INLINE_FUNCTION
+    virtual double computeLocalBasisFunction(const int node, const double xi, const double eta) const = 0;
+
+    // Compute jacobian for an element
+    KOKKOS_INLINE_FUNCTION
+    virtual double computeJacobian(const double xi, const double eta) const = 0;
+
+    // Compute element stiffness matrix
+    KOKKOS_INLINE_FUNCTION
+    virtual void computeElementStiffnessMatrix(double* stiffness) const = 0;
+
+    // Compute element load vector
+    KOKKOS_INLINE_FUNCTION
+    virtual void computeElementLoadVector(double* load) const = 0;
+
+	// Set material property k
+	KOKKOS_INLINE_FUNCTION
+	void setMaterialProperty(double k) {k_ = k;}
+
+	// Get material property k
+	KOKKOS_INLINE_FUNCTION
+	void getMaterialProperty(double k) {k_ = k;}
+	
+    
+};
+
+#endif // COMPUTING_AT_SCALE_ASSIGNMENT_ELEMENT_HPP

src/Mesh.h

@@ -59,6 +59,10 @@ class Mesh {
   KOKKOS_INLINE_FUNCTION
   MeshType GetMeshType() const { return meshType_; }
 
+  [[nodiscard]] size_t GetNumNodesPerElement() const {
+	return static_cast<int>(GetMeshType()); 
+  }
+
  private:
   // ***************************** Private Attributes
   // ***************************** //

src/QuadElement.hpp

@@ -0,0 +1,184 @@
+#ifndef COMPUTING_AT_SCALE_ASSIGNMENT_QUAD_ELEMENT_HPP
+#define COMPUTING_AT_SCALE_ASSIGNMENT_QUAD_ELEMENT_HPP
+
+#include "Element.hpp"
+
+class QuadElement : public Element {
+private:
+    // Gauss quadrature points and weights for quads 
+    static constexpr int numQuadPoints_ = 4;
+    static constexpr int numNodes_ = 4;
+    static constexpr double gaussPoint_ = 0.57735026919; // 1/sqrt(3)
+    
+    // Get quadrature point coordinates and weight
+    KOKKOS_INLINE_FUNCTION
+    void getQuadPoint(int q, double& xi, double& eta, double& weight) const {
+        if (q == 0) {
+            xi = -gaussPoint_; eta = -gaussPoint_; weight = 1.0;
+        } else if (q == 1) {
+            xi = gaussPoint_; eta = -gaussPoint_; weight = 1.0;
+        } else if (q == 2) {
+            xi = gaussPoint_; eta = gaussPoint_; weight = 1.0;
+        } else if (q == 3) {
+            xi = -gaussPoint_; eta = gaussPoint_; weight = 1.0;
+        }
+    }
+
+public:
+	KOKKOS_INLINE_FUNCTION
+    QuadElement(const Mesh& mesh, int elemIdx) : Element(mesh, elemIdx) {}
+    
+    KOKKOS_INLINE_FUNCTION
+    int getNumNodes() const override { return numNodes_; }
+    
+    KOKKOS_INLINE_FUNCTION
+    double computeLocalBasisFunction(const int node, const double xi, const double eta) const override {
+        switch(node) {
+            case 0: return 0.25 * (1.0 - xi) * (1.0 - eta);
+            case 1: return 0.25 * (1.0 + xi) * (1.0 - eta);
+            case 2: return 0.25 * (1.0 + xi) * (1.0 + eta);
+            case 3: return 0.25 * (1.0 - xi) * (1.0 + eta);
+            default: return 0.0;
+        }
+    }
+    
+    KOKKOS_INLINE_FUNCTION
+    void computeBasisGradient(const int node, const double xi, const double eta, 
+                             double& dN_dxi, double& dN_deta) const {
+        switch(node) {
+            case 0:
+                dN_dxi = -0.25 * (1.0 - eta);
+                dN_deta = -0.25 * (1.0 - xi);
+                break;
+            case 1:
+                dN_dxi = 0.25 * (1.0 - eta);
+                dN_deta = -0.25 * (1.0 + xi);
+                break;
+            case 2:
+                dN_dxi = 0.25 * (1.0 + eta);
+                dN_deta = 0.25 * (1.0 + xi);
+                break;
+            case 3:
+                dN_dxi = -0.25 * (1.0 + eta);
+                dN_deta = 0.25 * (1.0 - xi);
+                break;
+            default:
+                dN_dxi = 0.0;
+                dN_deta = 0.0;
+        }
+    }
+    
+    KOKKOS_INLINE_FUNCTION
+    double computeJacobian(const double xi, const double eta) const override {
+        // For quads, the Jacobian varies by position
+        double x[4], y[4];
+        for (int i = 0; i < 4; i++) {
+            x[i] = mesh_.GetCoordinate(elemIdx_, i, 0);
+            y[i] = mesh_.GetCoordinate(elemIdx_, i, 1);
+        }
+        
+        // Compute derivatives of x and y w.r.t. local coordinates
+        double dxdxi = 0.0, dxdeta = 0.0, dydxi = 0.0, dydeta = 0.0;
+        
+        for (int i = 0; i < 4; i++) {
+            double dN_dxi, dN_deta;
+            computeBasisGradient(i, xi, eta, dN_dxi, dN_deta);
+            
+            dxdxi += x[i] * dN_dxi;
+            dxdeta += x[i] * dN_deta;
+            dydxi += y[i] * dN_dxi;
+            dydeta += y[i] * dN_deta;
+        }
+        
+        return dxdxi * dydeta - dxdeta * dydxi;
+    }
+    
+    KOKKOS_INLINE_FUNCTION
+    void computeElementStiffnessMatrix(double* stiffness) const override {
+        
+        // Initialize stiffness matrix
+        for (int i = 0; i < numNodes_ * numNodes_; i++) {
+            stiffness[i] = 0.0;
+        }
+        
+        // Get coordinates of quadrilateral vertices
+        double x[4], y[4];
+        for (int i = 0; i < 4; i++) {
+            x[i] = mesh_.GetCoordinate(elemIdx_, i, 0);
+            y[i] = mesh_.GetCoordinate(elemIdx_, i, 1);
+        }
+        
+        // Integrate using Gauss quadrature
+        for (int q = 0; q < numQuadPoints_; q++) {
+            double xi, eta, weight;
+            getQuadPoint(q, xi, eta, weight);
+            
+            // Compute Jacobian at this quadrature point
+            double dxdxi = 0.0, dxdeta = 0.0, dydxi = 0.0, dydeta = 0.0;
+            
+            for (int n = 0; n < numNodes_; n++) {
+                double dN_dxi, dN_deta;
+                computeBasisGradient(n, xi, eta, dN_dxi, dN_deta);
+                
+                dxdxi += x[n] * dN_dxi;
+                dxdeta += x[n] * dN_deta;
+                dydxi += y[n] * dN_dxi;
+                dydeta += y[n] * dN_deta;
+            }
+            
+            double det_J = dxdxi * dydeta - dxdeta * dydxi;
+            double abs_det_J = det_J > 0 ? det_J : -det_J;
+            
+            // compute inverse of the jacobian
+			double invJ = 1/abs_det_J;
+            
+            // Compute contribution to stiffness matrix
+            for (int i = 0; i < numNodes_; i++) {
+                double dNi_dxi, dNi_deta;
+                computeBasisGradient(i, xi, eta, dNi_dxi, dNi_deta);
+                
+                double dNi_dx = dydeta * dNi_dxi - dydxi * dNi_deta ;
+                double dNi_dy = -dxdeta * dNi_dxi +  dxdxi * dNi_deta;
+                
+                for (int j = 0; j < numNodes_; j++) {
+                    double dNj_dxi, dNj_deta;
+                    computeBasisGradient(j, xi, eta, dNj_dxi, dNj_deta);
+                    
+                    double dNj_dx = dydeta * dNj_dxi  - dydxi * dNj_deta ;
+                    double dNj_dy = dydeta * dNj_dxi  +  dxdxi * dNj_deta ;
+                    
+                    stiffness[i * numNodes_ + j] += (dNi_dx * dNj_dx + dNi_dy * dNj_dy) * invJ * weight;
+                }
+            }
+        }
+        
+    }
+    
+    KOKKOS_INLINE_FUNCTION
+    void computeElementLoadVector(double* load) const override {
+        // Create load vector (4 entries)
+        
+        // Initialize load vector
+        for (int i = 0; i < numNodes_; i++) {
+            load[i] = 0.0;
+        }
+        
+        double f = 1.0;
+        
+        // Integrate load using quadrature
+        for (int q = 0; q < numQuadPoints_; q++) {
+            double xi, eta, weight;
+            getQuadPoint(q, xi, eta, weight);
+            
+            double det_J = computeJacobian(xi, eta);
+            double abs_det_J = det_J > 0 ? det_J : -det_J;
+            
+            for (int i = 0; i < numNodes_; i++) {
+                double phi = computeLocalBasisFunction(i, xi, eta);
+                load[i] += phi * f * weight * abs_det_J;
+            }
+        }
+    }
+};
+
+#endif // COMPUTING_AT_SCALE_ASSIGNMENT_QUAD_ELEMENT_HPP

src/StiffnessMatrix.cpp

@@ -138,3 +138,38 @@ void StiffnessMatrix::printStiffnessMatrix() const {
   }
   printf("\n");
 }
+
+std::vector<std::vector<double>> StiffnessMatrix::getDenseMatrix() const {
+  std::vector<std::vector<double>> denseMatrix(nDof_,
+                                               std::vector<double>(nDof_, 0.0));
+
+  auto csr_row_host = Kokkos::create_mirror_view(csrRowIds_);
+  auto csr_col_host = Kokkos::create_mirror_view(csrColIds_);
+  auto csr_val_host = Kokkos::create_mirror_view(csrValues_);
+  Kokkos::deep_copy(csr_row_host, csrRowIds_);
+  Kokkos::deep_copy(csr_col_host, csrColIds_);
+  Kokkos::deep_copy(csr_val_host, csrValues_);
+
+  for (int row = 0; row < nDof_; ++row) {
+    auto row_start = csr_row_host(row);
+    auto row_end = csr_row_host(row + 1);
+    for (int i = row_start; i < row_end; ++i) {
+      int col = csr_col_host(i);
+      double value = csr_val_host(i);
+      denseMatrix[row][col] += value;
+    }
+  }
+
+  return denseMatrix;
+}
+
+void StiffnessMatrix::printDenseMatrix() const {
+  auto denseMatrix = getDenseMatrix();
+  printf("Dense Matrix:\n");
+  for (const auto &row : denseMatrix) {
+    for (const auto &val : row) {
+      printf("%5.1f ", val);
+    }
+    printf("\n");
+  }
+}