Merge pull request #275 from arcaneframework/dev/mab-hex-quad-support-electrostatics

[feat] hexa/quad elemnt support

Author: mohd-afeef-badri

meshes/msh/box-rods.quad.msh

@@ -24,9 +24,11 @@ file(COPY "inputs" DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
 file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/meshes)
 set(MESH_FILES
   box-rods.msh
+  box-rods.quad.msh
   box-rod-circle.msh
   interdigital_capacitor.msh
   truncated_cube.msh
+  truncated_cube.hexa.msh
   spheres_in_sphere.msh
 )
 foreach(MESH_FILE IN LISTS MESH_FILES)
@@ -37,9 +39,13 @@ enable_testing()
 
 if(FEMUTILS_HAS_SOLVER_BACKEND_PETSC)
   add_test(NAME [electrostatics]box-rods COMMAND Electrostatics inputs/box-rods.arc)
+  add_test(NAME [electrostatics]box-rods_neumann COMMAND Electrostatics inputs/box-rods.neumann.arc)
+  add_test(NAME [electrostatics]box-rods_quad COMMAND Electrostatics inputs/box-rods.quad.arc)
+  add_test(NAME [electrostatics]box-rods_neumann_quad COMMAND Electrostatics inputs/box-rods.neumann.quad.arc)
   add_test(NAME [electrostatics]capacitor COMMAND Electrostatics inputs/Capacitor.arc)
   add_test(NAME [electrostatics]rod-circle COMMAND Electrostatics inputs/rod-circle.arc)
   add_test(NAME [electrostatics]spheres COMMAND Electrostatics inputs/spheres.arc)
+  add_test(NAME [electrostatics]truncated-cube_hexa COMMAND Electrostatics inputs/truncated_cube.hexa.arc)
   add_test(NAME [electrostatics]box-rods_petsc_monitor COMMAND Electrostatics
     -A,//fem/petsc-flags=-ksp_monitor
     inputs/box-rods.arc)
@@ -67,6 +73,8 @@ if(FEMUTILS_HAS_SOLVER_BACKEND_PETSC)
 
   if(FEMUTILS_HAS_PARALLEL_SOLVER AND MPIEXEC_EXECUTABLE)
     add_test(NAME [electrostatics]box-rods_2pu COMMAND ${MPIEXEC_EXECUTABLE} -n 2 ./Electrostatics inputs/box-rods.arc)
+    add_test(NAME [electrostatics]box-rods_quad_2pu COMMAND ${MPIEXEC_EXECUTABLE} -n 2 ./Electrostatics inputs/box-rods.quad.arc)
+    add_test(NAME [electrostatics]truncated-cube_hexa_2pu COMMAND ${MPIEXEC_EXECUTABLE} -n 2 ./Electrostatics inputs/truncated_cube.hexa.arc)
   endif()
 endif()
 

modules/electrostatics/CMakeLists.txt

@@ -0,0 +1,122 @@
+// -*- 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 electrostatics */
+/*---------------------------------------------------------------------------*/
+/*---------------------------------------------------------------------------*/
+
+/*---------------------------------------------------------------------------*/
+/**
+ * @brief Computes the element matrix for a quadrilateral element (QUAD4, ℙ1 FE).
+ *
+ * This function calculates the integral of:
+ *       𝑎(𝑢,𝑣) = ∫∫ (∂𝑢/∂𝑥 ∂𝑣/∂𝑥  + ∂𝑢/∂𝑦 ∂𝑣/∂𝑦)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;
+
+      // Integration weight
+      const Real integration_weight = detJ * w * w;
+
+      // stiffness matrix assembly
+      ae += (dxU ^ dxU) * integration_weight + (dyU ^ dyU) * integration_weight;
+    }
+  }
+  return ae;
+}
+
+/*---------------------------------------------------------------------------*/
+/**
+ * @brief Computes the element matrix for a hexahedral element (HEXA8, ℙ1 FE).
+ *
+ * This function calculates the integral of:
+ *       𝑎(𝑢,𝑣) = ∫∫∫ (∂𝑢/∂𝑥 ∂𝑣/∂𝑥 + ∂𝑢/∂𝑦 ∂𝑣/∂𝑦 + ∂𝑢/∂𝑧 ∂𝑣/∂𝑧)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;
+
+        // 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;
+      }
+    }
+  }
+  return ae;
+}

modules/electrostatics/ElementMatrixHexQuad.h

@@ -33,7 +33,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>
 
@@ -49,6 +49,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/electrostatics/Fem.axl

@@ -12,6 +12,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);
@@ -188,15 +190,36 @@ _assembleLinearOperatorCpu()
   auto node_dof(m_dofs_on_nodes.nodeDoFConnectivityView());
 
   if (options()->rho.isPresent()) {
-    Real qdot = -rho / epsilon;
-    ArcaneFemFunctions::BoundaryConditions2D::applyConstantSourceToRhs(qdot, mesh(), node_dof, m_node_coord, rhs_values);
+    Real f = -rho / epsilon;
+    if (mesh()->dimension() == 2) {
+      if (m_hex_quad_mesh)
+        ArcaneFemFunctions::BoundaryConditions2D::applyConstantSourceToRhsQuad4(f, mesh(), node_dof, m_node_coord, rhs_values);
+      else
+        ArcaneFemFunctions::BoundaryConditions2D::applyConstantSourceToRhs(f, mesh(), node_dof, m_node_coord, rhs_values);
+    }
+
+    if (mesh()->dimension() == 3) {
+      if (m_hex_quad_mesh)
+        ArcaneFemFunctions::BoundaryConditions3D::applyConstantSourceToRhsHexa8(f, mesh(), node_dof, m_node_coord, rhs_values);
+      else
+        ArcaneFemFunctions::BoundaryConditions3D::applyConstantSourceToRhs(f, mesh(), node_dof, m_node_coord, rhs_values);
+    }
   }
 
   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);
+        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);
 
       for (BC::IDirichletBoundaryCondition* bs : bc->dirichletBoundaryConditions())
         BCFunctions.applyDirichletToLhsAndRhs(bs, node_dof, m_node_coord, m_linear_system, rhs_values);
@@ -318,13 +341,15 @@ _assembleBilinearOperator()
 
   if (m_matrix_format == "DOK") {
     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); });
     else
-      _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;
@@ -422,15 +447,32 @@ _updateVariables()
   m_phi.synchronize();
 
   if (mesh()->dimension() == 2)
-    ENUMERATE_ (Cell, icell, allCells()) {
-      Cell cell = *icell;
-      m_E[cell] = ArcaneFemFunctions::FeOperation2D::computeGradientTria3(cell, m_node_coord, m_phi);
-    }
-  else
-    ENUMERATE_ (Cell, icell, allCells()) {
-      Cell cell = *icell;
-      m_E[cell] = ArcaneFemFunctions::FeOperation3D::computeGradientTetra4(cell, m_node_coord, m_phi);
-    }
+    if (m_hex_quad_mesh)
+      ENUMERATE_ (Cell, icell, allCells()) {
+        Cell cell = *icell;
+        Real3 grad = ArcaneFemFunctions::FeOperation2D::computeGradientQuad4(cell, m_node_coord, m_phi);
+        m_E[cell] = -grad.x * grad.x - grad.y * grad.y;
+      }
+    else
+      ENUMERATE_ (Cell, icell, allCells()) {
+        Cell cell = *icell;
+        Real3 grad = ArcaneFemFunctions::FeOperation2D::computeGradientTria3(cell, m_node_coord, m_phi);
+        m_E[cell] = -grad.x * grad.x - grad.y * grad.y;
+      }
+
+  if (mesh()->dimension() == 3)
+    if (m_hex_quad_mesh)
+      ENUMERATE_ (Cell, icell, allCells()) {
+        Cell cell = *icell;
+        Real3 grad = ArcaneFemFunctions::FeOperation3D::computeGradientHexa8(cell, m_node_coord, m_phi);
+        m_E[cell] = -grad.x * grad.x - grad.y * grad.y - grad.z * grad.z;
+      }
+    else
+      ENUMERATE_ (Cell, icell, allCells()) {
+        Cell cell = *icell;
+        Real3 grad = ArcaneFemFunctions::FeOperation3D::computeGradientTetra4(cell, m_node_coord, m_phi);
+        m_E[cell] = -grad.x * grad.x - grad.y * grad.y - grad.z * grad.z;
+      }
 
   m_E.synchronize();
 
@@ -460,11 +502,10 @@ _validateResults()
   if (allNodes().size() < 200)
     ENUMERATE_ (Node, inode, allNodes()) {
       Node node = *inode;
-      info() << "phi[" << node.localId() << "][" << node.uniqueId() << "] = " << m_phi[node];
+      info() << "phi[" << node.uniqueId() << "] = " << m_phi[node];
     }
 
   String filename = options()->resultFile();
-  info() << "ValidateResultFile filename=" << filename;
 
   if (!filename.empty())
     checkNodeResultFile(traceMng(), filename, m_phi, 1.0e-4, 1.0e-16);

modules/electrostatics/FemModule.cc

@@ -89,6 +89,7 @@ class FemModule
   bool m_assemble_linear_system = true;
   bool m_solve_linear_system = true;
   bool m_cross_validation = true;
+  bool m_hex_quad_mesh = false;
 
   void _doStationarySolve();
   void _getMaterialParameters();
@@ -100,7 +101,8 @@ class FemModule
 
   RealMatrix<3, 3> _computeElementMatrixTria3(Cell cell);
   RealMatrix<4, 4> _computeElementMatrixTetra4(Cell cell);
-
+  RealMatrix<4, 4> _computeElementMatrixQuad4(Cell cell);
+  RealMatrix<8, 8> _computeElementMatrixHexa8(Cell cell);
   template <int N>
   void _assembleBilinear(const std::function<RealMatrix<N, N>(const Cell&)>& compute_element_matrix);
 };

modules/electrostatics/FemModule.h

@@ -0,0 +1,20 @@
+260  0.20717053527226
+261  0.383610894594714
+262  0.387911285847584
+263  -0.407050243122318
+264  -0.425425596473062
+265  -0.705231441654697
+266  0.326249166788217
+267  0.349280160680684
+268  -0.364613077656847
+269  -0.0201922027407905
+270  -0.697845286522503
+271  -0.115572177922491
+272  -0.630863639596694
+273  0.302314253921705
+274  -0.162425200359167
+275  -0.114418588238129
+276  -0.2948195045131
+277  -0.375718403629807
+278  -0.77074423851435
+279  -0.17490104704671

modules/electrostatics/check/box-rods.quad.txt

@@ -0,0 +1,23 @@
+330  -0.454092321745147
+337  -0.567484524008041
+338  -0.596074907200828
+339  -0.602676785764029
+340  -0.630921162792379
+341  -0.702324315388119
+342  -0.635109745204455
+343  -0.666487677467807
+344  -0.758449272836748
+345  -0.601902087988295
+346  -0.623165072062167
+347  -0.674979546518178
+348  -0.643542875728247
+349  -0.679525088730678
+350  -0.775568041352925
+351  -0.687504081612848
+352  -0.401359711200829
+353  -0.697643007268312
+354  -0.842631723211244
+355  -0.490300612287456
+356  -0.885986586154468
+359  -0.749894190756057
+360  -0.903598044401985