[feat] quad4/hex8 support

Author: mohd-afeef-badri

femutils/ArcaneFemFunctions.h

@@ -710,6 +710,25 @@ class ArcaneFemFunctions
 
       return { grad_x, grad_y, 0.0 };
     }
+
+    /*---------------------------------------------------------------------------*/
+    /**
+     * @brief Computes the shape functions for a Quad4 element at a given point (ξ, η).
+     * 
+     * @param xi The ξ coordinate of the evaluation point (-1 to 1).
+     * @param eta The η coordinate of the evaluation point (-1 to 1).
+     * @return A RealVector<4> containing the shape functions {𝑁₁, 𝑁₂, 𝑁₃, 𝑁₄}.
+     */
+    /*---------------------------------------------------------------------------*/
+    static inline RealVector<4> computeShapeFunctionsQuad4(Real xi, Real eta)
+    {
+      RealVector<4> N;
+      N(0) = 0.25 * (1. - xi) * (1. - eta); // 𝑁₁
+      N(1) = 0.25 * (1. + xi) * (1. - eta); // 𝑁₂
+      N(2) = 0.25 * (1. + xi) * (1. + eta); // 𝑁₃
+      N(3) = 0.25 * (1. - xi) * (1. + eta); // 𝑁₄
+      return N;
+    }
   };
 
   /*---------------------------------------------------------------------------*/
@@ -908,7 +927,7 @@ class ArcaneFemFunctions
 
     /*---------------------------------------------------------------------------*/
     /**
-     * @brief Holds information for a Quad4 element at a single Gauss point.
+     * @brief Holds information for a Hexa8 element at a single Gauss point.
      *
      * This includes the gradients of the shape functions in the physical space (𝑥,𝑦,𝑧)
      * and the determinant of the Jacobian matrix.
@@ -1073,6 +1092,39 @@ class ArcaneFemFunctions
 
       return { grad_x, grad_y, grad_z };
     }
+
+
+    /*---------------------------------------------------------------------------*/
+    /**
+     * @brief Computes the shape functions for a Hexa8 element at a given point (ξ,η,ζ).
+     * 
+     * @param xi The ξ coordinate of the evaluation point (-1 to 1).
+     * @param eta The η coordinate of the evaluation point (-1 to 1).
+     * @param zeta The ζ coordinate of the evaluation point (-1 to 1).
+     * @return A RealVector<8> containing the shape functions {𝑁₁, 𝑁₂, 𝑁₃, 𝑁₄, 𝑁₅, 𝑁₆, 𝑁₇, 𝑁₈}.
+     */
+    /*---------------------------------------------------------------------------*/
+    static inline RealVector<8> computeShapeFunctionsHexa8(Real xi, Real eta, Real zeta)
+    {
+      RealVector<8> N;
+      const Real one_minus_eta = 1.0 - eta;
+      const Real one_plus_eta = 1.0 + eta;
+      const Real one_minus_xi = 1.0 - xi;
+      const Real one_plus_xi = 1.0 + xi;
+      const Real one_minus_zeta = 1.0 - zeta;
+      const Real one_plus_zeta = 1.0 + zeta;
+
+      N(0) = 0.125 * one_minus_xi * one_minus_eta * one_minus_zeta; // 𝑁₁
+      N(1) = 0.125 * one_plus_xi * one_minus_eta * one_minus_zeta; // 𝑁₂
+      N(2) = 0.125 * one_plus_xi * one_plus_eta * one_minus_zeta; // 𝑁₃
+      N(3) = 0.125 * one_minus_xi * one_plus_eta * one_minus_zeta; // 𝑁₄
+      N(4) = 0.125 * one_minus_xi * one_minus_eta * one_plus_zeta; // 𝑁₅
+      N(5) = 0.125 * one_plus_xi * one_minus_eta * one_plus_zeta; // 𝑁₆
+      N(6) = 0.125 * one_plus_xi * one_plus_eta * one_plus_zeta; // 𝑁₇
+      N(7) = 0.125 * one_minus_xi * one_plus_eta * one_plus_zeta; // 𝑁₈
+
+      return N;
+    }
   };
 
     /*---------------------------------------------------------------------------*/

meshes/med/sphere_in_sphere.py

@@ -0,0 +1,62 @@
+#!/usr/bin/env python
+
+###
+### This file is generated automatically by SALOME v9.14.0 with dump python functionality
+###
+
+import sys
+import salome
+
+salome.salome_init()
+import salome_notebook
+notebook = salome_notebook.NoteBook()
+
+###
+### SHAPER component
+###
+
+from salome.shaper import model
+
+model.begin()
+partSet = model.moduleDocument()
+
+### Create Part
+Part_1 = model.addPart(partSet)
+Part_1_doc = Part_1.document()
+
+### Create Sphere
+Sphere_1 = model.addSphere(Part_1_doc, model.selection("VERTEX", "PartSet/Origin"), 0.1)
+
+### Create Sphere
+Sphere_2 = model.addSphere(Part_1_doc, model.selection("VERTEX", "PartSet/Origin"), 0.01)
+
+### Create Cut
+Cut_1 = model.addCut(Part_1_doc, [model.selection("COMPOUND", "all-in-Sphere_1")], [model.selection("SOLID", "Sphere_2_1")], keepSubResults = True)
+
+### Create Group
+Group_1 = model.addGroup(Part_1_doc, "Faces", [model.selection("FACE", "Sphere_2_1/Face_1")])
+Group_1.setName("inner")
+Group_1.result().setName("inner")
+
+### Create Group
+Group_2 = model.addGroup(Part_1_doc, "Faces", [model.selection("FACE", "Sphere_1_1/Face_1")])
+Group_2.setName("outer")
+Group_2.result().setName("outer")
+
+### Create Group
+Group_3 = model.addGroup(Part_1_doc, "Solids", [model.selection("SOLID", "Cut_1_1")])
+Group_3.setName("vol")
+Group_3.result().setName("vol")
+
+model.end()
+
+###
+### SHAPERSTUDY component
+###
+
+model.publishToShaperStudy()
+import SHAPERSTUDY
+Cut_1_1, inner, outer, vol, = SHAPERSTUDY.shape(model.featureStringId(Cut_1))
+
+if salome.sg.hasDesktop():
+  salome.sg.updateObjBrowser()

meshes/med/truncated_cube.py

@@ -10,7 +10,6 @@
 salome.salome_init()
 import salome_notebook
 notebook = salome_notebook.NoteBook()
-sys.path.insert(0, r'/home/catA/mb258512/Install/arcanefem/arcanefem/meshes/med')
 
 ###
 ### SHAPER component

meshes/msh/sphere_in_sphere.hexa.msh

@@ -21,7 +21,9 @@ file(COPY "inputs" DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
 file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/meshes)
 set(MESH_FILES
   sub.msh
+  sub.quad.msh
   sub_3d.msh
+  sphere_in_sphere.hexa.msh
 )
 foreach(MESH_FILE IN LISTS MESH_FILES)
     file(COPY ${MSH_DIR}/${MESH_FILE} DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/meshes)
@@ -32,7 +34,9 @@ target_link_libraries(Acoustics PUBLIC FemUtils)
 enable_testing()
 
 add_test(NAME [acoustics]2D_submarine COMMAND Acoustics inputs/sub.arc)
+add_test(NAME [acoustics]2D_submarine_quad COMMAND Acoustics inputs/sub.quad.arc)
 add_test(NAME [acoustics]3D_sphere COMMAND Acoustics inputs/3d_sub.arc)
+add_test(NAME [acoustics]3D_sphere_hexa COMMAND Acoustics inputs/3d_sphere_in_sphere.hexa.arc)
 
 
 if(FEMUTILS_HAS_SOLVER_BACKEND_PETSC)

meshes/msh/sub.quad.msh

@@ -0,0 +1,129 @@
+// -*- tab-width: 2; indent-tabs-mode: nil; coding: utf-8-with-signature -*-
+//-----------------------------------------------------------------------------
+// Copyright 2000-2025 CEA (www.cea.fr) IFPEN (www.ifpenergiesnouvelles.com)
+// See the top-level COPYRIGHT file for details.
+// SPDX-License-Identifier: Apache-2.0
+//-----------------------------------------------------------------------------
+/*---------------------------------------------------------------------------*/
+/* ElementMatrixHexQuad.h                                      (C) 2022-2025 */
+/*                                                                           */
+/* Contains functions to compute the FEM element matrices for Acoustics      */
+/*---------------------------------------------------------------------------*/
+/*---------------------------------------------------------------------------*/
+
+/*---------------------------------------------------------------------------*/
+/**
+ * @brief Computes the element matrix for a quadrilateral element (QUAD4, ℙ1 FE).
+ *
+ * This function calculates the integral of:
+ *       𝑎(𝑢,𝑣) = ∫∫ (∂𝑢/∂𝑥 ∂𝑣/∂𝑥  + ∂𝑢/∂𝑦 ∂𝑣/∂𝑦)dΩ  + ∫∫ 𝑘²𝑢𝑣 dΩ
+ *
+ * Steps involved:
+ * 1. Define Gauss points (2x2) and weights.
+ * 2. Loop over Gauss points to compute the gradients in physical space
+ *    and the determinant of the Jacobian, via computeGradientsAndJacobianQuad4.
+ * 3. Compute the integration weight.
+ * 4. Assemble the element matrix using the computed gradients.
+ *
+ * @param cell The cell for which the element matrix is computed.
+ * @return The computed element matrix.
+ */
+/*---------------------------------------------------------------------------*/
+
+RealMatrix<4, 4> FemModule::_computeElementMatrixQuad4(Cell cell)
+{
+  // Gauss points and weights for 2x2 quadrature
+  constexpr Real gp[2] = { -M_SQRT1_3, M_SQRT1_3 }; // [-1/sqrt(3) , 1/sqrt(3)]
+  constexpr Real w = 1.0;
+
+  // Initialize the element matrix
+  RealMatrix<4, 4> ae;
+  ae.fill(0.0);
+
+  // Loop over Gauss points
+  for (Int8 ixi = 0; ixi < 2; ++ixi) {
+    for (Int8 ieta = 0; ieta < 2; ++ieta) {
+
+      // Get the coordinates of the Gauss point in natural coordinates (ξ,η,ζ)
+      const Real xi = gp[ixi];
+      const Real eta = gp[ieta];
+
+      // Get shape function gradients w.r.t (𝑥,𝑦) and determinant of Jacobian
+      const auto gp_info = ArcaneFemFunctions::FeOperation2D::computeGradientsAndJacobianQuad4(cell, m_node_coord, xi, eta);
+      const RealVector<4>& dxU = gp_info.dN_dx;
+      const RealVector<4>& dyU = gp_info.dN_dy;
+      const Real detJ = gp_info.det_j;
+
+      // Shape functions at the Gauss point
+      RealVector<4> N = ArcaneFemFunctions::FeOperation2D::computeShapeFunctionsQuad4(xi, eta);
+
+      // Integration weight
+      const Real integration_weight = detJ * w * w;
+
+      // stiffness matrix assembly
+      ae += (dxU ^ dxU) * integration_weight + (dyU ^ dyU) * integration_weight +(N ^ N) * m_kc2 * integration_weight;
+    }
+  }
+  return ae;
+}
+
+/*---------------------------------------------------------------------------*/
+/**
+ * @brief Computes the element matrix for a hexahedral element (HEXA8, ℙ1 FE).
+ *
+ * This function calculates the integral of:
+ *       𝑎(𝑢,𝑣) = ∫∫∫ (∂𝑢/∂𝑥 ∂𝑣/∂𝑥 + ∂𝑢/∂𝑦 ∂𝑣/∂𝑦 + ∂𝑢/∂𝑧 ∂𝑣/∂𝑧)dΩ  + ∫∫∫ 𝑘²𝑢𝑣 dΩ
+ *
+ * Steps involved:
+ * 1. Define Gauss points (2x2x2) and weights.
+ * 2. Loop over Gauss points to compute the gradients in physical space (𝑥,𝑦,𝑧)
+ *    and the determinant of the Jacobian via computeGradientsAndJacobianHexa8.
+ * 3. Compute the integration weight.
+ * 4. Assemble the element matrix using the computed gradients.
+ *
+ * @param cell The cell for which the element matrix is computed.
+ * @return The computed element matrix.
+ */
+/*---------------------------------------------------------------------------*/
+
+RealMatrix<8, 8> FemModule::_computeElementMatrixHexa8(Cell cell)
+{
+  // 2x2x2 Gauss points and weights for [-1,1]^3
+  constexpr Real gp[2] = { -M_SQRT1_3, M_SQRT1_3 }; // -1/sqrt(3), 1/sqrt(3)
+  constexpr Real w = 1.0;
+
+  // Initialize the element matrix
+  RealMatrix<8, 8> ae;
+  ae.fill(0.0);
+
+  // Loop over Gauss points
+  for (Int8 ixi = 0; ixi < 2; ++ixi) {
+    for (Int8 ieta = 0; ieta < 2; ++ieta) {
+      for (Int8 izeta = 0; izeta < 2; ++izeta) {
+
+        // Get the coordinates of Gauss points in natural coordinates (ξ,η,ζ)
+        const Real xi = gp[ixi];
+        const Real eta = gp[ieta];
+        const Real zeta = gp[izeta];
+
+        // Get shape function gradients w.r.t (𝑥,𝑦,𝑧) and determinant of Jacobian
+        const auto gp_info = ArcaneFemFunctions::FeOperation3D::computeGradientsAndJacobianHexa8(cell, m_node_coord, xi, eta, zeta);
+        const RealVector<8>& dxU = gp_info.dN_dx;
+        const RealVector<8>& dyU = gp_info.dN_dy;
+        const RealVector<8>& dzU = gp_info.dN_dz;
+        const Real detJ = gp_info.det_j;
+
+        // Shape functions at the Gauss point
+        RealVector<8> N = ArcaneFemFunctions::FeOperation3D::computeShapeFunctionsHexa8(xi, eta, zeta);
+
+        // Integration weight
+        const Real integration_weight = detJ * w * w;
+
+        // Assemble element matrix (variational form)
+        ae += (dxU ^ dxU) * integration_weight + (dyU ^ dyU) * integration_weight + (dzU ^ dzU) * integration_weight +
+        (N ^ N) * m_kc2 * integration_weight;
+      }
+    }
+  }
+  return ae;
+}

modules/acoustics/CMakeLists.txt

@@ -37,7 +37,7 @@
       <description>Which matrix format to use DOK|BSR|AF-BSR.</description>
     </simple>
 
-    <simple name="petsc-flags" type="string" default="" optional="true">
+    <simple name="petsc-flags" type="string" default="-mat_type baij" optional="true">
       <description>Flags for PETSc from commandline.</description>
     </simple>
 
@@ -53,6 +53,10 @@
       <description>Boolean to enable cross validation.</description>
     </simple>
 
+    <simple name="hex-quad-mesh" type="bool" default="false" optional="true">
+      <description>Boolean to use hexa or quad element mesh.</description>
+    </simple>
+
     <!-- Linear system service instance -->
     <service-instance name="linear-system" type="Arcane::FemUtils::IDoFLinearSystemFactory" default="AlephLinearSystem" />
 

modules/acoustics/ElementMatrixHexQuad.h

@@ -13,6 +13,7 @@
 
 #include "FemModule.h"
 #include "ElementMatrix.h"
+#include "ElementMatrixHexQuad.h"
 
 /*---------------------------------------------------------------------------*/
 /**
@@ -36,6 +37,7 @@ startInit()
   m_solve_linear_system = options()->solveLinearSystem();
   m_cross_validation = options()->crossValidation();
   m_petsc_flags = options()->petscFlags();
+  m_hex_quad_mesh = options()->hexQuadMesh();
 
   elapsedTime = platform::getRealTime() - elapsedTime;
   ArcaneFemFunctions::GeneralFunctions::printArcaneFemTime(traceMng(), "initialize", elapsedTime);
@@ -144,23 +146,19 @@ _assembleLinearOperator()
 
   const auto node_dof(m_dofs_on_nodes.nodeDoFConnectivityView());
 
-  auto applyBoundaryConditions = [&](auto BCFunctions) {
-    BC::IArcaneFemBC* bc = options()->boundaryConditions();
-    if (bc) {
-      for (BC::INeumannBoundaryCondition* bs : bc->neumannBoundaryConditions()) {
-        BCFunctions.applyNeumannToRhs(bs, node_dof, m_node_coord, rhs_values);
-      }
-    }
-  };
-
-  // Apply the correct boundary conditions based on mesh dimension
-  if (mesh()->dimension() == 3) {
-    using BCFunctions = ArcaneFemFunctions::BoundaryConditions3D;
-    applyBoundaryConditions(BCFunctions());
-  }
-  else {
-    using BCFunctions = ArcaneFemFunctions::BoundaryConditions2D;
-    applyBoundaryConditions(BCFunctions());
+  BC::IArcaneFemBC* bc = options()->boundaryConditions();
+  if (bc) {
+    for (BC::INeumannBoundaryCondition* bs : bc->neumannBoundaryConditions())
+      if (mesh()->dimension() == 2)
+        if (m_hex_quad_mesh)
+          ArcaneFemFunctions::BoundaryConditions2D::applyNeumannToRhsQuad4(bs, node_dof, m_node_coord, rhs_values);
+        else
+          ArcaneFemFunctions::BoundaryConditions2D::applyNeumannToRhs(bs, node_dof, m_node_coord, rhs_values);
+      else
+        if (m_hex_quad_mesh)
+          ArcaneFemFunctions::BoundaryConditions3D::applyNeumannToRhsHexa8(bs, node_dof, m_node_coord, rhs_values);
+        else
+          ArcaneFemFunctions::BoundaryConditions3D::applyNeumannToRhs(bs, node_dof, m_node_coord, rhs_values);
   }
 
   elapsedTime = platform::getRealTime() - elapsedTime;
@@ -180,13 +178,16 @@ _assembleBilinearOperator()
   Real elapsedTime = platform::getRealTime();
 
   if (mesh()->dimension() == 3)
-    _assembleBilinear<4>([this](const Cell& cell) {
-      return _computeElementMatrixTetra4(cell);
-    });
+    if(m_hex_quad_mesh)
+      _assembleBilinear<8>([this](const Cell& cell) { return _computeElementMatrixHexa8(cell); });
+    else
+      _assembleBilinear<4>([this](const Cell& cell) { return _computeElementMatrixTetra4(cell); });
+
   if (mesh()->dimension() == 2)
-    _assembleBilinear<3>([this](const Cell& cell) {
-      return _computeElementMatrixTria3(cell);
-    });
+    if(m_hex_quad_mesh)
+      _assembleBilinear<4>([this](const Cell& cell) { return _computeElementMatrixQuad4(cell); });
+    else
+      _assembleBilinear<3>([this](const Cell& cell) { return _computeElementMatrixTria3(cell); });
 
   elapsedTime = platform::getRealTime() - elapsedTime;
   ArcaneFemFunctions::GeneralFunctions::printArcaneFemTime(traceMng(), "lhs-matrix-assembly", elapsedTime);