Merge pull request #289 from arcaneframework/dev/mab-hex-quad-support-elastodynamics

Dev/mab hex quad support elastodynamics

Author: mohd-afeef-badri

modules/elastodynamics/FemModule.cc

@@ -155,7 +155,7 @@ _getParameters()
 
   //--------- material parameter ---------//
   E = options()->E(); // Youngs modulus 𝐸
-  nu = options()->nu(); // Poission ratio ν
+  nu = options()->nu(); // Poisson ratio ν
   rho = options()->rho(); // Density ρ
 
   mu = E / (2 * (1 + nu)); // lame parameter μ
@@ -392,23 +392,28 @@ _assembleBilinearOperatorTetra4Gpu()
 }
 
 /*---------------------------------------------------------------------------*/
+/**
+ * @brief Assembles the FEM bilinear operator for 2D problems.
+ *
+ * This method assembles the FEM stiffness matrix by iterating over each cell,
+ * computing the element stiffness matrix using the provided function, and
+ * populating the global stiffness matrix accordingly.
+ *
+ * @tparam N The number of nodes per element (e.g., 3 for triangles, 4 for quadrilaterals).
+ * @param compute_element_matrix A function that computes the element stiffness matrix for a given cell.
+ */
 /*---------------------------------------------------------------------------*/
 
+template <int N>
 void FemModule::
-_assembleBilinearOperatorTria3()
+_assembleBilinearOperator2d(const std::function<RealMatrix<N, N>(const Cell&)>& compute_element_matrix)
 {
   auto node_dof(m_dofs_on_nodes.nodeDoFConnectivityView());
 
   ENUMERATE_ (Cell, icell, allCells()) {
     Cell cell = *icell;
 
-    auto K_e = _computeElementMatrixTria3(cell); // element stiffness matrix
-    // assemble elementary matrix into  the global one elementary terms are
-    // positioned into  K according  to the rank of associated  node in the
-    // mesh.nodes list  and according the dof number. Here  for  each  node
-    // two dofs  exists [u1,u2]. For each TRIA3 there are 3 nodes hence the
-    // elementary stiffness matrix size is (3*2 x 3*2)=(6x6). We will  fill
-    // this below in 4 at a time.
+    auto K_e = compute_element_matrix(cell);
     Int32 n1_index = 0;
     for (Node node1 : cell.nodes()) {
       Int32 n2_index = 0;
@@ -417,13 +422,13 @@ _assembleBilinearOperatorTria3()
         Real v2 = K_e(2 * n1_index, 2 * n2_index + 1);
         Real v3 = K_e(2 * n1_index + 1, 2 * n2_index);
         Real v4 = K_e(2 * n1_index + 1, 2 * n2_index + 1);
-        // m_k_matrix(node1.localId(), node2.localId()) += v;
+
         if (node1.isOwn()) {
           DoFLocalId node1_dof1 = node_dof.dofId(node1, 0);
           DoFLocalId node1_dof2 = node_dof.dofId(node1, 1);
           DoFLocalId node2_dof1 = node_dof.dofId(node2, 0);
           DoFLocalId node2_dof2 = node_dof.dofId(node2, 1);
-          //          m_linear_system.matrixAddValue(node_dof.dofId(node1, 0), node_dof.dofId(node2, 0), v);
+
           m_linear_system.matrixAddValue(node1_dof1, node2_dof1, v1);
           m_linear_system.matrixAddValue(node1_dof1, node2_dof2, v2);
           m_linear_system.matrixAddValue(node1_dof2, node2_dof1, v3);
@@ -437,17 +442,28 @@ _assembleBilinearOperatorTria3()
 }
 
 /*---------------------------------------------------------------------------*/
+/**
+ * @brief Assembles the FEM bilinear operator for 3D problems.
+ *
+ * This method assembles the FEM stiffness matrix by iterating over each cell,
+ * computing the element stiffness matrix using the provided function, and
+ * populating the global stiffness matrix accordingly.
+ *
+ * @tparam N The number of nodes per element (e.g., 4 for tetrahedra, 8 for hexahedra).
+ * @param compute_element_matrix A function that computes the element stiffness matrix for a given cell.
+ */
 /*---------------------------------------------------------------------------*/
 
+template <int N>
 void FemModule::
-_assembleBilinearOperatorTetra4()
+_assembleBilinearOperator3d(const std::function<RealMatrix<N, N>(const Cell&)>& compute_element_matrix)
 {
   auto node_dof(m_dofs_on_nodes.nodeDoFConnectivityView());
 
   ENUMERATE_ (Cell, icell, allCells()) {
     Cell cell = *icell;
 
-    auto K_e = _computeElementMatrixTetra4(cell);
+    auto K_e = compute_element_matrix(cell);
     Int32 n1_index = 0;
     for (Node node1 : cell.nodes()) {
       Int32 n2_index = 0;
@@ -500,16 +516,25 @@ _assembleBilinearOperator()
   Real elapsedTime = platform::getRealTime();
 
   if (t <= dt) {
-    if (mesh()->dimension() == 2)
-      if(m_matrix_format == "BSR" || m_matrix_format == "AF-BSR")
+    if (m_matrix_format == "DOK") {
+      if (mesh()->dimension() == 2) {
+        _assembleBilinearOperator2d<6>([this](const Cell& cell) { return _computeElementMatrixTria3(cell); });
+      }
+      if (mesh()->dimension() == 3) {
+        _assembleBilinearOperator3d<12>([this](const Cell& cell) { return _computeElementMatrixTetra4(cell); });
+      }
+    }
+    else if (m_matrix_format == "BSR" || m_matrix_format == "AF-BSR") {
+      if (mesh()->dimension() == 2) {
         _assembleBilinearOperatorTria3Gpu();
-      else
-        _assembleBilinearOperatorTria3();
-    else
-      if(m_matrix_format == "BSR" || m_matrix_format == "AF-BSR")
+      }
+      if (mesh()->dimension() == 3) {
         _assembleBilinearOperatorTetra4Gpu();
-      else
-        _assembleBilinearOperatorTetra4();
+      }
+    }
+    else {
+      ARCANE_FATAL("Unsupported matrix type, only DOK| BSR | AF-BSR is supported.");
+    }
   }
 
   elapsedTime = platform::getRealTime() - elapsedTime;

modules/elastodynamics/FemModule.h

@@ -141,8 +141,6 @@ class FemModule
   void _validateResults();
   void _readCaseTables();
   void _assembleBilinearOperator();
-  void _assembleBilinearOperatorTria3();
-  void _assembleBilinearOperatorTetra4();
 
   inline void _applyTraction(VariableDoFReal& rhs_values, const IndexedNodeDoFConnectivityView& node_dof);
   inline void _applyDirichlet(VariableDoFReal& rhs_values, const IndexedNodeDoFConnectivityView& node_dof);
@@ -151,6 +149,11 @@ class FemModule
   RealMatrix<6, 6> _computeElementMatrixTria3(Cell cell);
   RealMatrix<12, 12> _computeElementMatrixTetra4(Cell cell);
 
+  template <int N>
+  void _assembleBilinearOperator2d(const std::function<RealMatrix<N, N>(const Cell&)>& compute_element_matrix);
+
+  template <int N>
+  void _assembleBilinearOperator3d(const std::function<RealMatrix<N, N>(const Cell&)>& compute_element_matrix);
 };
 
 #endif