Merge pull request #225 from arcaneframework/dev/mab-fourier-3D

fourier 3D with tests
arcaneframework · Feb 27, 2025 · 81ff754 · 81ff754
2 parents 671c720 + bb218a8
commit 81ff754
Show file tree

Hide file tree

Showing 5 changed files with 270 additions and 21 deletions.
diff --git a/femutils/ArcaneFemFunctions.h b/femutils/ArcaneFemFunctions.h
@@ -241,6 +241,29 @@ class ArcaneFemFunctions
       return { Center_x, Center_y, 0 };
     }
 
+    /*---------------------------------------------------------------------------*/
+    /**
+     * @brief Computes the barycenter (centroid) of a tetrahedron.
+     *
+     * This method calculates the barycenter of a tetrahedron defined by its nodes.
+     * The barycenter is computed as the average of the vertices' coordinates.
+     */
+    /*---------------------------------------------------------------------------*/
+
+    static inline Real3 computeBaryCenterTetra4(Cell cell, const VariableNodeReal3& node_coord)
+    {
+      Real3 vertex0 = node_coord[cell.nodeId(0)];
+      Real3 vertex1 = node_coord[cell.nodeId(1)];
+      Real3 vertex2 = node_coord[cell.nodeId(2)];
+      Real3 vertex3 = node_coord[cell.nodeId(3)];
+
+      Real Center_x = (vertex0.x + vertex1.x + vertex2.x + vertex3.x) / 4.;
+      Real Center_y = (vertex0.y + vertex1.y + vertex2.y + vertex3.x) / 4.;
+      Real Center_z = (vertex0.z + vertex1.z + vertex2.z + vertex3.z) / 4.;
+
+      return { Center_x, Center_y, Center_z };
+    }
+
     /*---------------------------------------------------------------------------*/
     /**
      * @brief Computes the length of the edge defined by a given face.
@@ -864,6 +887,36 @@ class ArcaneFemFunctions
       }
     }
 
+    /*---------------------------------------------------------------------------*/
+    /**
+     * @brief Applies a manufactured source term to the RHS vector.
+     *
+     * This method adds a manufactured source term to the RHS vector for each
+     * node in the mesh. The contribution to each node is weighted by the area of
+     * the cell and evenly distributed among the nodes of the cell.
+     *
+     * @param [IN]  qdot       : The constant source term.
+     * @param [IN]  mesh       : The mesh containing all cells.
+     * @param [IN]  node_dof   : DOF connectivity view.
+     * @param [IN]  node_coord : The coordinates of the nodes.
+     * @param [OUT] rhs_values : The RHS values to update.
+     */
+    /*---------------------------------------------------------------------------*/
+
+    static inline void applyManufacturedSourceToRhs(IBinaryMathFunctor<Real, Real3, Real>* manufactured_source, IMesh* mesh, const IndexedNodeDoFConnectivityView& node_dof, const VariableNodeReal3& node_coord, VariableDoFReal& rhs_values)
+    {
+      ENUMERATE_ (Cell, icell, mesh->allCells()) {
+        Cell cell = *icell;
+        Real volume = ArcaneFemFunctions::MeshOperation::computeVolumeTetra4(cell, node_coord);
+        Real3 bcenter = ArcaneFemFunctions::MeshOperation::computeBaryCenterTetra4(cell, node_coord);
+
+        for (Node node : cell.nodes()) {
+          if (node.isOwn())
+            rhs_values[node_dof.dofId(node, 0)] += manufactured_source->apply(volume / cell.nbNode(), bcenter);
+        }
+      }
+    }
+
     /*---------------------------------------------------------------------------*/
     /**
      * @brief Applies Dirichlet boundary conditions to RHS and LHS.
@@ -993,6 +1046,40 @@ class ArcaneFemFunctions
         }
       }
     }
+
+    /*---------------------------------------------------------------------------*/
+    /**
+     * @brief Applies Manufactured Dirichlet boundary conditions to RHS and LHS.
+     *
+     * Updates the LHS matrix and RHS vector to enforce the Dirichlet.
+     *
+     * - For LHS matrix `A`, the diagonal term for the Dirichlet DOF is set to `P`.
+     * - For RHS vector `b`, the Dirichlet DOF term is scaled by `P`.
+     *
+     * @param [IN]  manufactured_dirichlet   : External function for Dirichlet.
+     * @param [IN]  group           : Group of all external faces.
+     * @param [IN]  bs              : Boundary condition values.
+     * @param [IN]  node_dof        : DOF connectivity view.
+     * @param [IN]  node_coord      : Node coordinates.
+     * @param [OUT] m_linear_system : Linear system for LHS.
+     * @param [OUT] rhs_values RHS  : RHS values to update.
+     */
+    /*---------------------------------------------------------------------------*/
+    static inline void applyManufacturedDirichletToLhsAndRhs(IBinaryMathFunctor<Real, Real3, Real>* manufactured_dirichlet, Real /*lambda*/, const FaceGroup& group, BC::IManufacturedSolution* bs, const IndexedNodeDoFConnectivityView& node_dof, const VariableNodeReal3& node_coord, DoFLinearSystem& m_linear_system, VariableDoFReal& rhs_values)
+    {
+      Real Penalty = bs->getPenalty();
+
+      ENUMERATE_ (Face, iface, group) {
+        for (Node node : iface->nodes()) {
+          if (node.isOwn()) {
+            m_linear_system.matrixSetValue(node_dof.dofId(node, 0), node_dof.dofId(node, 0), Penalty);
+            double tt = 1.;
+            Real u_g = Penalty * manufactured_dirichlet->apply(tt, node_coord[node]);
+            rhs_values[node_dof.dofId(node, 0)] = u_g;
+          }
+        }
+      }
+    }
   };
 
   /*---------------------------------------------------------------------------*/

diff --git a/modules/fourier/CMakeLists.txt b/modules/fourier/CMakeLists.txt
@@ -34,6 +34,7 @@ set(MESH_FILES
   multi-material.msh
   plancher.quad4.msh
   square_-2pi_to_2pi.msh
+  bar_dynamic_3D.msh
 )
 foreach(MESH_FILE IN LISTS MESH_FILES)
     file(COPY ${MSH_DIR}/${MESH_FILE} DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/meshes)
@@ -73,6 +74,13 @@ if(FEMUTILS_HAS_SOLVER_BACKEND_PETSC)
   add_test(NAME [Fourier]manufacture_solution COMMAND Fourier
     -A,UsingDotNet=1
     inputs/manufacture.solution.arc)
+  add_test(NAME [Fourier]conduction_3D COMMAND Fourier
+    inputs/conduction.3D.arc)
+  add_test(NAME [Fourier]manufacture_solution_3D COMMAND Fourier
+    -A,UsingDotNet=1
+    -A,//meshes/mesh/filename=meshes/bar_dynamic_3D.msh 
+    -A,//fem/mesh-type=TETRA4
+    inputs/manufacture.solution.arc)
 
   if(FEMUTILS_HAS_PARALLEL_SOLVER AND MPIEXEC_EXECUTABLE)
     add_test(NAME [Fourier]conduction_2p COMMAND ${MPIEXEC_EXECUTABLE} -n 2 ./Fourier
@@ -82,6 +90,11 @@ if(FEMUTILS_HAS_SOLVER_BACKEND_PETSC)
     add_test(NAME [Fourier]manufacture_solution_2p COMMAND ${MPIEXEC_EXECUTABLE} -n 2 ./Fourier
       -A,UsingDotNet=1
       inputs/manufacture.solution.arc)
+    add_test(NAME [Fourier]manufacture_solution_3D_2p COMMAND ${MPIEXEC_EXECUTABLE} -n 2 ./Fourier
+      -A,UsingDotNet=1
+      -A,//meshes/mesh/filename=meshes/bar_dynamic_3D.msh 
+      -A,//fem/mesh-type=TETRA4
+      inputs/manufacture.solution.arc)
     if(FEMTEST_HAS_GMSH_TEST)
       add_test(NAME [Fourier]conduction_10k_2p COMMAND ${MPIEXEC_EXECUTABLE} -n 2 ./Fourier
         -A,//meshes/mesh/filename=meshes/plancher.10k.msh

diff --git a/modules/fourier/FemModule.cc b/modules/fourier/FemModule.cc
@@ -202,31 +202,39 @@ _assembleLinearOperator()
 
   auto node_dof(m_dofs_on_nodes.nodeDoFConnectivityView());
 
-  if (options()->qdot.isPresent())
-    ArcaneFemFunctions::BoundaryConditions2D::applyConstantSourceToRhs(qdot, mesh(), node_dof, m_node_coord, rhs_values);
-
   BC::IArcaneFemBC* bc = options()->boundaryConditions();
 
   if (bc) {
-    for (BC::INeumannBoundaryCondition* bs : bc->neumannBoundaryConditions())
-      ArcaneFemFunctions::BoundaryConditions2D::applyNeumannToRhs(bs, node_dof, m_node_coord, rhs_values);
-
-    for (BC::IDirichletBoundaryCondition* bs : bc->dirichletBoundaryConditions())
-      ArcaneFemFunctions::BoundaryConditions2D::applyDirichletToLhsAndRhs(bs, node_dof, m_node_coord, m_linear_system, rhs_values);
-
-    for (BC::IManufacturedSolution* bs : bc->manufacturedSolutions()) {
-      if (bs->getManufacturedSource()) {
-        ARCANE_CHECK_POINTER(m_manufactured_source);
-        info() << "Apply manufactured Source condition to all cells";
-        ArcaneFemFunctions::BoundaryConditions2D::applyManufacturedSourceToRhs(m_manufactured_source, mesh(), node_dof, m_node_coord, rhs_values);
-      }
-      if (bs->getManufacturedDirichlet()) {
-        ARCANE_CHECK_POINTER(m_manufactured_dirichlet);
-        info() << "Apply manufactured dirichlet condition to all borders";
-        FaceGroup group = mesh()->outerFaces();
-        ArcaneFemFunctions::BoundaryConditions2D::applyManufacturedDirichletToLhsAndRhs(m_manufactured_dirichlet, lambda, group, bs, node_dof, m_node_coord, m_linear_system, rhs_values);
+    auto applyBoundaryConditions = [&](auto BCFunctions) {
+      if (options()->qdot.isPresent())
+        BCFunctions.applyConstantSourceToRhs(qdot, mesh(), node_dof, m_node_coord, rhs_values);
+
+      BC::IArcaneFemBC* bc = options()->boundaryConditions();
+      for (BC::INeumannBoundaryCondition* bs : bc->neumannBoundaryConditions())
+        BCFunctions.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);
+
+      for (BC::IManufacturedSolution* bs : bc->manufacturedSolutions()) {
+        if (bs->getManufacturedSource()) {
+          ARCANE_CHECK_POINTER(m_manufactured_source);
+          info() << "Apply manufactured Source condition to all cells";
+          BCFunctions.applyManufacturedSourceToRhs(m_manufactured_source, mesh(), node_dof, m_node_coord, rhs_values);
+        }
+        if (bs->getManufacturedDirichlet()) {
+          ARCANE_CHECK_POINTER(m_manufactured_dirichlet);
+          info() << "Apply manufactured dirichlet condition to all borders";
+          FaceGroup group = mesh()->outerFaces();
+          BCFunctions.applyManufacturedDirichletToLhsAndRhs(m_manufactured_dirichlet, lambda, group, bs, node_dof, m_node_coord, m_linear_system, rhs_values);
+        }
       }
-    }
+    };
+
+    if (mesh()->dimension() == 3)
+      applyBoundaryConditions(ArcaneFemFunctions::BoundaryConditions3D());
+    else
+      applyBoundaryConditions(ArcaneFemFunctions::BoundaryConditions2D());
   }
 
   elapsedTime = platform::getRealTime() - elapsedTime;

diff --git a/modules/fourier/check/test_conduction_3D.txt b/modules/fourier/check/test_conduction_3D.txt
@@ -0,0 +1,64 @@
+1 12
+2 55
+3 55
+4 55
+5 55
+6 12
+7 12
+8 12
+9 50.7216511903551
+10 46.4387425490711
+11 42.1509493047928
+12 37.8582773154181
+13 33.5608291716414
+14 29.2581845749637
+15 24.950861190741
+16 20.6388656181528
+17 16.3220581180551
+18 16.3219078995032
+19 20.6387256401124
+20 24.9509397115797
+21 29.2582772718981
+22 33.5608291616598
+23 37.8581842921888
+24 42.1508586952028
+25 46.4386526503503
+26 50.7215780969766
+27 50.7216511903551
+28 46.4387425490711
+29 42.1509493047928
+30 37.858277315418
+31 33.5608291716414
+32 29.2581845749637
+33 24.9508714701844
+34 20.638744473218
+35 16.3220014135257
+36 16.321953690758
+37 20.6388647258524
+38 24.950956569867
+39 29.2582772718981
+40 33.5608291616598
+41 37.8581842921888
+42 42.1508586952028
+43 46.4386526503503
+44 50.7215780969766
+45 31.4101765280069
+46 48.5809477492477
+47 35.7102611938306
+48 40.0053634820613
+49 44.2955976300534
+50 52.8613623960177
+51 18.480906619199
+52 14.1613710092729
+53 27.1053100470222
+54 22.7955731592032
+55 35.7101764872987
+56 14.1614036258102
+57 18.4809848361747
+58 31.4102611531306
+59 27.1053624167412
+60 22.7955960075257
+61 48.5808986275508
+62 40.0053089816822
+63 44.2955443989129
+64 52.8614231112239
diff --git a/modules/fourier/inputs/conduction.3D.arc b/modules/fourier/inputs/conduction.3D.arc
@@ -0,0 +1,77 @@
+<?xml version="1.0"?>
+<!--
+  Case configuration for a Fourier analysis simulation.
+  This file includes settings for:
+    - General simulation parameters
+    - Mesh configuration details
+    - Finite Element Method (FEM) configurations
+    - Post-processing options
+-->
+<case codename="Fourier" xml:lang="en" codeversion="1.0">
+
+  <!--
+    Arcane-specific settings:
+      - title: Descriptive name for the simulation case.
+      - timeloop: Specifies the time-stepping loop used in this Fourier simulation.
+  -->
+  <arcane>
+    <title>Fouriers equation FEM code</title>
+    <timeloop>FourierLoop</timeloop>
+  </arcane>
+
+  <!--
+    Mesh configuration:
+      - filename: The path to the mesh file used in the simulation.
+  -->
+  <meshes>
+    <mesh>
+      <filename>meshes/bar_dynamic_3D.msh</filename>
+    </mesh>
+  </meshes>
+
+  <!--
+    FEM (Finite Element Method) settings:
+      - lambda: Thermal conductivity or diffusivity coefficient.
+      - qdot: Heat source term or volumetric heat generation.
+      - result-file: File where simulation results will be saved.
+      - boundary-conditions: Defines the boundary conditions for the simulation.
+        - dirichlet: Fixed value boundary condition with penalty enforcement for specified surfaces.
+        - neumann: Flux or gradient boundary condition for specified surfaces.
+  -->
+  <fem>
+    <lambda>0.023</lambda>
+    <qdot>1.123e-2</qdot>
+    <mesh-type>TETRA4</mesh-type>
+    <result-file>check/test_conduction_3D.txt</result-file>
+    <boundary-conditions>
+      <dirichlet>
+        <enforce-Dirichlet-method>Penalty</enforce-Dirichlet-method>
+        <surface>surfaceleft</surface>
+        <value>55.0</value>
+      </dirichlet>
+      <dirichlet>
+        <enforce-Dirichlet-method>Penalty</enforce-Dirichlet-method>
+        <surface>surfaceright</surface>
+        <value>12.0</value>
+      </dirichlet>
+      <neumann>
+        <surface>sidesurfaces</surface>
+        <value>0.0</value>
+      </neumann>
+    </boundary-conditions>
+  </fem>
+
+  <!--
+    Post-processing settings:
+      - output-period: Defines how often (in simulation steps) the output is generated.
+      - format: Specifies the post-processing format, here using VtkHdfV2.
+      - output: Defines the variables to be included in the post-processing output.
+  -->
+  <arcane-post-processing>
+   <output-period>1</output-period>
+   <output>
+     <variable>U</variable>
+   </output>
+  </arcane-post-processing>
+
+</case>