Merge pull request arcaneframework#260 from arcaneframework/dev/mab-acoustics-3d

mohd-afeef-badri · web-flow · commit c11dac516a28 · 2025-04-04T16:14:36.000+02:00
acoustics in 3d
diff --git a/meshes/msh/sub_3d.geo b/meshes/msh/sub_3d.geo
@@ -0,0 +1,39 @@
+//-----------------------------------------------------------------------------
+//
+// Name       : sub.geo
+// Author     : Mohd Afeef BADRI
+// Date       : 04 / April / 2025
+//
+// ----------------------------------------------------------------------------
+// Comment    : simple submarine problem mesh for ill posed acoustics
+//
+//
+// Usage      : gmsh sub_3d.geo -2 -format msh41 -clmin 1 -clmax 1
+//              gmsh sub_3d.geo -2 -format msh41 -setnumber MeshSize 1.0 -setnumber Rfactor 2.0 
+//
+//-----------------------------------------------------------------------------
+
+SetFactory("OpenCASCADE");
+
+DefineConstant[ meshSize = {1.0, Min 1e-8, Max 100, Step 1, Name "MeshSize"} ];
+DefineConstant[ Rfactor  = {1.0, Min 1, Max 100, Step 1, Name "Rfactor"} ];
+meshSize /= Rfactor;
+
+//==============================================================================
+// Outer Sphere (Water Domain)
+//==============================================================================
+R_sphere = .1;
+Sphere(1) = {0, 0, 0, R_sphere};
+
+//==============================================================================
+// Inner Torus (Cavity to be subtracted)
+//==============================================================================
+Sphere(2) = {0, 0, 0, R_sphere/10};
+
+BooleanDifference{ Volume{1}; Delete; }{ Volume{2}; Delete; }
+
+Physical Surface("outer", 10) = {3};
+Physical Surface("inner", 11) = {2};
+Physical Volume("water", 12) = {1};
+
+Mesh.MshFileVersion = 4.1;
diff --git a/meshes/msh/sub_3d.msh b/meshes/msh/sub_3d.msh
diff --git a/modules/acoustics/CMakeLists.txt b/modules/acoustics/CMakeLists.txt
@@ -21,6 +21,7 @@ file(COPY "inputs" DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
 file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/meshes)
 set(MESH_FILES
   sub.msh
+  sub_3d.msh
 )
 foreach(MESH_FILE IN LISTS MESH_FILES)
     file(COPY ${MSH_DIR}/${MESH_FILE} DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/meshes)
@@ -31,6 +32,8 @@ target_link_libraries(Acoustics PUBLIC FemUtils)
 enable_testing()
 
 add_test(NAME [Acoustics]2D_submarine COMMAND Acoustics inputs/sub.arc)
+add_test(NAME [Acoustics]3D_sphere COMMAND Acoustics inputs/3d_sub.arc)
+
 
 if(FEMUTILS_HAS_SOLVER_BACKEND_PETSC)
   add_test(NAME [Acoustics]parse_and_exit COMMAND Acoustics
@@ -42,6 +45,11 @@ endif()
 
 if(FEMUTILS_HAS_SOLVER_BACKEND_HYPRE)
   add_test(NAME [Acoustics]2D_submarine_hypre COMMAND Acoustics inputs/sub.hypre.arc)
+  add_test(NAME [Acoustics]3D_sphere_hypre COMMAND Acoustics
+    -A,//fem/linear-system/@name=AlephLinearSystem
+    -A,//fem/linear-system/solver-backend=hypre
+    -A,//fem/linear-system/solver-method=gmres
+    inputs/3d_sub.arc)
   if(FEMUTILS_HAS_PARALLEL_SOLVER AND MPIEXEC_EXECUTABLE)
     add_test(NAME [Acoustics]2D_submarine_hypre_2p COMMAND ${MPIEXEC_EXECUTABLE} -n 2 ./Acoustics inputs/sub.hypre.arc)
   endif()
diff --git a/modules/acoustics/ElementMatrix.h b/modules/acoustics/ElementMatrix.h
@@ -35,5 +35,33 @@ _computeElementMatrixTria3(Cell cell)
   Real3 dxU = ArcaneFemFunctions::FeOperation2D::computeGradientXTria3(cell, m_node_coord);
   Real3 dyU = ArcaneFemFunctions::FeOperation2D::computeGradientYTria3(cell, m_node_coord);
 
-  return -area * (dxU ^ dxU) - area * (dyU ^ dyU) + m_kc2 * area * (1/ 12.) *massMatrix(U,U);
+  return -area * (dxU ^ dxU) - area * (dyU ^ dyU) + m_kc2 * area * (1/ 12.) * massMatrix(U,U);
+}
+
+/*---------------------------------------------------------------------------*/
+/*
+  * @brief Computes the element matrix for a tetrahedral element (ℙ1 FE).
+  *
+  * This function calculates the integral of the expression:
+  *       a(𝑢,𝑣) = ∫∫∫ -(∂𝑢/∂𝑥 ∂𝑣/∂𝑥  + ∂𝑢/∂𝑦 ∂𝑣/∂𝑦 + ∂𝑢/∂𝑧 ∂𝑣/∂𝑧)dΩ + ∫∫∫ (𝑘²𝑢𝑣)dΩ
+  * 
+  * Steps involved:
+  * 1. Calculate the volume of the tetrahedron.
+  * 2. Compute the gradients of the shape functions.
+  * 3. Return a(𝑢,𝑣);
+  */
+/*---------------------------------------------------------------------------*/
+
+RealMatrix<4, 4> FemModule::
+_computeElementMatrixTetra4(Cell cell)
+{
+  Real volume = ArcaneFemFunctions::MeshOperation::computeVolumeTetra4(cell, m_node_coord);
+
+  RealVector<4> U = { 1, 1, 1, 1 };
+
+  Real4 dxU = ArcaneFemFunctions::FeOperation3D::computeGradientXTetra4(cell, m_node_coord);
+  Real4 dyU = ArcaneFemFunctions::FeOperation3D::computeGradientYTetra4(cell, m_node_coord);
+  Real4 dzU = ArcaneFemFunctions::FeOperation3D::computeGradientZTetra4(cell, m_node_coord);
+
+  return -volume * (dxU ^ dxU) - volume * (dyU ^ dyU) - volume * (dzU ^ dzU) + m_kc2 * volume * (1 / 20.) * massMatrix(U,U);
 }
diff --git a/modules/acoustics/FemModule.cc b/modules/acoustics/FemModule.cc
@@ -144,11 +144,23 @@ _assembleLinearOperator()
 
   const auto node_dof(m_dofs_on_nodes.nodeDoFConnectivityView());
 
-  BC::IArcaneFemBC* bc = options()->boundaryConditions();
-  if(bc){
-    for (BC::INeumannBoundaryCondition* bs : bc->neumannBoundaryConditions()){
-      ArcaneFemFunctions::BoundaryConditions2D::applyNeumannToRhs(bs, 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);
+      }
     }
+  };
+
+  // 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());
   }
 
   elapsedTime = platform::getRealTime() - elapsedTime;
@@ -157,11 +169,7 @@ _assembleLinearOperator()
 
 /*---------------------------------------------------------------------------*/
 /**
- * @brief Assembles the bilinear operator matrix for triangular elements.
- *
- * This method computes and assembles the global  matrix A for the FEM linear
- * system. For each cell (triangle), it calculates the local element matrix &
- * adds the contributions to the global matrix based on the nodes of the cell.
+ * @brief Calls the right function for LHS assembly
  */
 /*---------------------------------------------------------------------------*/
 
@@ -171,28 +179,51 @@ _assembleBilinearOperator()
   info() << "[ArcaneFem-Info] Started module _assembleBilinearOperator()";
   Real elapsedTime = platform::getRealTime();
 
+  if (mesh()->dimension() == 3)
+    _assembleBilinear<4>([this](const Cell& cell) {
+      return _computeElementMatrixTetra4(cell);
+    });
+  if (mesh()->dimension() == 2)
+    _assembleBilinear<3>([this](const Cell& cell) {
+      return _computeElementMatrixTria3(cell);
+    });
+
+  elapsedTime = platform::getRealTime() - elapsedTime;
+  ArcaneFemFunctions::GeneralFunctions::printArcaneFemTime(traceMng(), "lhs-matrix-assembly", elapsedTime);
+}
+
+/*---------------------------------------------------------------------------*/
+/**
+ * @brief Assembles the bilinear operator matrix for the FEM linear system.
+ *
+ * The method performs the following steps:
+ *   1. Computes element matrix using provided `compute_element_matrix` function.
+ *   2. Assembles global matrix by adding contributions from each cell's element
+ *      matrix to the corresponding entries in the global matrix.
+ */
+/*---------------------------------------------------------------------------*/
+template <int N>
+void FemModule::
+_assembleBilinear(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 matrix
+    auto K_e = compute_element_matrix(cell); // element matrix based on the provided function
     Int32 n1_index = 0;
     for (Node node1 : cell.nodes()) {
       Int32 n2_index = 0;
       for (Node node2 : cell.nodes()) {
         Real v = K_e(n1_index, n2_index);
-        if (node1.isOwn()) {
+        if (node1.isOwn())
           m_linear_system.matrixAddValue(node_dof.dofId(node1, 0), node_dof.dofId(node2, 0), v);
-        }
         ++n2_index;
       }
       ++n1_index;
     }
   }
-
-  elapsedTime = platform::getRealTime() - elapsedTime;
-  ArcaneFemFunctions::GeneralFunctions::printArcaneFemTime(traceMng(), "lhs-matrix-assembly", elapsedTime);
 }
 
 /*---------------------------------------------------------------------------*/
@@ -272,7 +303,7 @@ _validateResults()
   if (allNodes().size() < 200)
     ENUMERATE_ (Node, inode, allNodes()) {
       Node node = *inode;
-      info() << "u[" << node.localId() << "][" << node.uniqueId() << "] = " << m_u[node];
+      info() << "u[" << node.uniqueId() << "] = " << m_u[node];
     }
 
   String filename = options()->resultFile();
diff --git a/modules/acoustics/FemModule.h b/modules/acoustics/FemModule.h
@@ -92,6 +92,10 @@ class FemModule
   void _validateResults();
 
   RealMatrix<3, 3> _computeElementMatrixTria3(Cell cell);
+  RealMatrix<4, 4> _computeElementMatrixTetra4(Cell cell);
+
+  template <int N>
+  void _assembleBilinear(const std::function<RealMatrix<N, N>(const Cell&)>& compute_element_matrix);
 };
 
 #endif
diff --git a/modules/acoustics/check/sphere_3d.txt b/modules/acoustics/check/sphere_3d.txt
@@ -0,0 +1,56 @@
+100  0.000287725021892968
+101  0.000239487109431728
+102  6.51015195239284e-05
+103  0.000418360390566261
+104  0.000695100947909542
+105  0.00011841053666585
+106  0.000745230270268164
+107  0.00155214956876286
+108  -0.00218100907895198
+109  0.000502915198987293
+110  0.00206940546040083
+111  -0.000100127323893553
+112  -0.0012251385696835
+113  -0.00263987079126473
+114  -0.00236862276135281
+115  -0.000155901277895176
+116  -0.000239852535688998
+117  -0.00305089164066475
+118  0.000190996771418783
+119  -0.00545688210886177
+120  0.000579180510396156
+121  -0.000219216523123141
+122  -0.00141865600256941
+123  -0.00865160369166539
+124  -8.99509445725068e-05
+125  -0.00375564689083683
+126  -0.000125389596609557
+127  -0.00760882841914332
+128  0.000254953852276821
+129  -0.00649672558348422
+130  -0.000466848324349205
+131  -0.00011253003924276
+132  0.000371406994764287
+133  0.000105571321312606
+134  -0.000222638935562521
+135  -0.0240200497028155
+136  0.000436231362124546
+137  0.000413657809187159
+138  0.00014359469311985
+139  -0.013834859215644
+140  -0.00414893245513014
+141  -5.30925790943364e-05
+142  -0.00588048336219698
+143  0.000344793471343447
+144  -0.017452978722355
+145  -0.0109365904183506
+146  0.00178745648624318
+147  4.19894278062304e-05
+148  -0.000150093689767147
+149  4.19582047111049e-05
+150  -0.000131206554784793
+151  5.60241137246923e-05
+152  -3.55502845594731e-05
+153  0.000596956497932044
+154  3.0625429695194e-05
+155  -0.0140370936494099
diff --git a/modules/acoustics/inputs/3d_sub.arc b/modules/acoustics/inputs/3d_sub.arc
@@ -0,0 +1,65 @@
+<?xml version="1.0"?>
+<!--
+  Case configuration for an acoustics simulation. 
+  The following sections define:
+    - General simulation settings
+    - Mesh configurations details
+    - Finite Element Method (FEM) configurations
+    - Post-processing options
+-->
+<case codename="Acoustics" xml:lang="en" codeversion="1.0">
+
+  <!--
+    Arcane-specific settings:
+      - title: The name or description of the case.
+      - timeloop: Defines the time-stepping loop for the simulation.
+  -->
+  <arcane>
+    <title>Sphere in sphere Toy Case</title>
+    <timeloop>AcousticsLoop</timeloop>
+  </arcane>
+
+
+  <!--
+    Mesh configurations:
+      - filename: Path to the mesh file used in the simulation.
+  -->
+  <meshes>
+    <mesh>
+      <filename>meshes/sub_3d.msh</filename>
+    </mesh>
+  </meshes>
+
+
+  <!--
+    FEM (Finite Element Method) settings:
+      - kc2: Coefficient used in the FEM calculations.
+      - boundary-conditions: Defines neumann boundary conditions for the simulation.
+      - linear-system: Specifies the linear system solver to use.
+      - result-file: File for validation (optional)
+  -->
+  <fem>
+    <kc2>18e5</kc2>
+    <boundary-conditions>
+      <neumann>
+        <surface>inner</surface>
+        <value>11e2</value>
+      </neumann>
+    </boundary-conditions>
+    <linear-system name="SequentialBasicLinearSystem" />
+    <result-file>check/sphere_3d.txt</result-file>
+  </fem>
+
+  <!--
+    Post-processing settings:
+      - output-period: Defines how often output should be generated.
+      - output: Specifies which variables are to be outputted.
+  -->
+  <arcane-post-processing>
+   <output-period>1</output-period>
+   <output>
+     <variable>U</variable>
+   </output>
+  </arcane-post-processing>
+
+</case>
