diff --git a/electrostatics/Fem.axl b/electrostatics/Fem.axl
index d7acb5b9..b78f2ab8 100644
--- a/electrostatics/Fem.axl
+++ b/electrostatics/Fem.axl
@@ -9,9 +9,6 @@
FEM variable phi on nodes for electrostatic potential
-
- Boolean which is true if Dirichlet node
-
electric field vector on cells
@@ -32,16 +29,6 @@
Type of mesh provided to the solver
-
-
- Method via which Dirichlet boundary condition is imposed
-
-
-
-
- Penalty value for enforcing Dirichlet condition
-
-
+
+
+ Method via which Dirichlet boundary condition is imposed
+
+
+
+
+ Penalty value for enforcing Dirichlet condition
+
+
+
+
+
+
+
+ Dirichlet boundary condition
+
+
+
+ Function for Dirichlet boundary condition
+
+
+
+
+ Method via which Dirichlet boundary condition is imposed
+
+
+
+
+ Penalty value for enforcing Dirichlet condition
+
+
+
+
+ Function for manufactured source term condition
+
+
@@ -115,6 +143,16 @@
Value of the point Dirichlet condition
+
+
+ Method via which Dirichlet boundary condition is imposed
+
+
+
+
+ Penalty value for enforcing Dirichlet condition
+
+
-#include
-#include
-
-#include
-#include
-#include
-#include
-#include
-
-#include "IDoFLinearSystemFactory.h"
-#include "Fem_axl.h"
-#include "FemUtils.h"
-#include "DoFLinearSystem.h"
-#include "FemDoFsOnNodes.h"
+#include "FemModule.h"
/*---------------------------------------------------------------------------*/
+/**
+ * @brief Initializes the FemModule at the start of the simulation.
+ *
+ * - initializes degrees of freedom (DoFs) on nodes.
+ * - builds support for manufactured test case (optional).
+ */
/*---------------------------------------------------------------------------*/
-using namespace Arcane;
-using namespace Arcane::FemUtils;
-
-/*---------------------------------------------------------------------------*/
-/*---------------------------------------------------------------------------*/
-/*!
- * \brief Module Fem.
- */
-class FemModule
-: public ArcaneFemObject
+void FemModule::
+startInit()
{
- public:
-
- explicit FemModule(const ModuleBuildInfo& mbi)
- : ArcaneFemObject(mbi)
- , m_dofs_on_nodes(mbi.subDomain()->traceMng())
- {
- ICaseMng* cm = mbi.subDomain()->caseMng();
- cm->setTreatWarningAsError(true);
- cm->setAllowUnkownRootElelement(false);
- }
-
- public:
-
- //! Method called at each iteration
- void compute() override;
-
- //! Method called at the beginning of the simulation
- void startInit() override;
-
- VersionInfo versionInfo() const override
- {
- return VersionInfo(1, 0, 0);
- }
-
- private:
-
- Real rho;
- Real epsilon;
- Real ElementNodes;
-
- DoFLinearSystem m_linear_system;
- IItemFamily* m_dof_family = nullptr;
- FemDoFsOnNodes m_dofs_on_nodes;
-
- private:
-
- void _doStationarySolve();
- void _getMaterialParameters();
- void _assembleBilinearOperatorTRIA3();
- void _solve();
- void _getE();
- void _initBoundaryconditions();
- void _assembleLinearOperator();
- void _applyDirichletBoundaryConditions();
- void _checkResultFile();
- FixedMatrix<3, 3> _computeElementMatrixTRIA3(Cell cell);
- Real _computeAreaTriangle3(Cell cell);
- Real _computeEdgeLength2(Face face);
- Real2 _computeEdgeNormal2(Face face);
- Real2 _computeDxDyOfRealTRIA3(Cell cell);
+ info() << "Module Fem INIT";
-};
+ m_dofs_on_nodes.initialize(mesh(), 1);
+ m_dof_family = m_dofs_on_nodes.dofFamily();
+}
/*---------------------------------------------------------------------------*/
+/**
+ * @brief Performs the main computation for the FemModule.
+ *
+ * - Stops the time loop after 1 iteration since the equation is steady state.
+ * - Resets, configures, and initializes the linear system.
+ * - Executes the stationary solve.
+ */
/*---------------------------------------------------------------------------*/
void FemModule::
@@ -113,8 +56,7 @@ compute()
// This is only used for the first call.
{
StringList string_list;
- string_list.add("-trmalloc");
- string_list.add("-log_trace");
+ string_list.add("-ksp_monitor");
CommandLineArguments args(string_list);
m_linear_system.setSolverCommandLineArguments(args);
}
@@ -124,42 +66,37 @@ compute()
}
/*---------------------------------------------------------------------------*/
-/*---------------------------------------------------------------------------*/
-
-void FemModule::
-startInit()
-{
- info() << "Module Fem INIT";
-
- m_dofs_on_nodes.initialize(mesh(), 1);
- m_dof_family = m_dofs_on_nodes.dofFamily();
-
- _initBoundaryconditions();
-}
-
-/*---------------------------------------------------------------------------*/
+/**
+ * @brief Performs a stationary solve for the FEM system.
+ *
+ * This method follows a sequence of steps to solve FEM system:
+ * 1. _getMaterialParameters() Retrieves material parameters via
+ * 2. _assembleBilinearOperator() Assembles the FEM matrix A
+ * 3. _assembleLinearOperator() Assembles the FEM RHS vector b
+ * 4. _solve() Solves for solution vector u = A^-1*b
+ * 5. _validateResults() Regression test
+ */
/*---------------------------------------------------------------------------*/
void FemModule::
_doStationarySolve()
{
- // # get material parameters
_getMaterialParameters();
-
- // Assemble the FEM bilinear operator (LHS - matrix A)
- _assembleBilinearOperatorTRIA3();
-
- // Assemble the FEM linear operator (RHS - vector b)
+ _assembleBilinearOperator();
_assembleLinearOperator();
-
- // # T=linalg.solve(K,RHS)
_solve();
-
- // Check results
- _checkResultFile();
+ _validateResults();
}
/*---------------------------------------------------------------------------*/
+/**
+ * @brief Retrieves and sets the material parameters for the simulation.
+ *
+ * This method initializes:
+ * - material properties:
+ * # charge density (`rho`)
+ * # freespace permittivity (`epsilon`)
+ */
/*---------------------------------------------------------------------------*/
void FemModule::
@@ -168,56 +105,6 @@ _getMaterialParameters()
info() << "Get material parameters...";
rho = options()->rho(); // charge density
epsilon = options()->epsilon(); // freespace permittivity
- ElementNodes = 3.; // 3 nodes of triangle
-}
-
-/*---------------------------------------------------------------------------*/
-/*---------------------------------------------------------------------------*/
-
-void FemModule::
-_initBoundaryconditions()
-{
- info() << "Init boundary conditions...";
-
- info() << "Apply boundary conditions";
- _applyDirichletBoundaryConditions();
-}
-
-/*---------------------------------------------------------------------------*/
-/*---------------------------------------------------------------------------*/
-
-void FemModule::
-_applyDirichletBoundaryConditions()
-{
- // Handle all the Dirichlet boundary conditions.
- // In the 'arc' file, there are in the following format:
- //
- // Haut
- // 21.0
- //
-
- for (const auto& bs : options()->dirichletBoundaryCondition()) {
- FaceGroup group = bs->surface();
- Real value = bs->value();
- info() << "Apply Dirichlet boundary condition surface=" << group.name() << " v=" << value;
- ENUMERATE_ (Face, iface, group) {
- for (Node node : iface->nodes()) {
- m_phi[node] = value;
- m_phi_dirichlet[node] = true;
- }
- }
- }
-
- for (const auto& bs : options()->dirichletPointCondition()) {
- NodeGroup group = bs->node();
- Real value = bs->value();
- info() << "Apply Dirichlet point condition node=" << group.name() << " v=" << value;
- ENUMERATE_ (Node, inode, group) {
- Node node = *inode;
- m_phi[node] = value;
- m_phi_dirichlet[node] = true;
- }
- }
}
/*---------------------------------------------------------------------------*/
@@ -231,203 +118,22 @@ void FemModule::
_assembleLinearOperator()
{
info() << "Assembly of FEM linear operator ";
- info() << "Applying Dirichlet boundary condition via penalty method ";
- // Temporary variable to keep values for the RHS part of the linear system
VariableDoFReal& rhs_values(m_linear_system.rhsVariable());
rhs_values.fill(0.0);
auto node_dof(m_dofs_on_nodes.nodeDoFConnectivityView());
- if (options()->enforceDirichletMethod() == "Penalty") {
-
- //----------------------------------------------
- // penalty method to enforce Dirichlet BC
- //----------------------------------------------
- // Let 'P' be the penalty term and let 'i' be the set of DOF for which
- // Dirichlet condition needs to be applied
- //
- // - For LHS matrix A the diag term corresponding to the Dirichlet DOF
- // a_{i,i} = 1. * P
- //
- // - For RHS vector b the term that corresponds to the Dirichlet DOF
- // b_{i} = b_{i} * P
- //----------------------------------------------
-
- info() << "Applying Dirichlet boundary condition via "
- << options()->enforceDirichletMethod() << " method ";
-
- Real Penalty = options()->penalty(); // 1.0e30 is the default
-
- ENUMERATE_ (Node, inode, ownNodes()) {
- NodeLocalId node_id = *inode;
- if (m_phi_dirichlet[node_id]) {
- DoFLocalId dof_id = node_dof.dofId(*inode, 0);
- m_linear_system.matrixSetValue(dof_id, dof_id, Penalty);
- Real u_g = Penalty * m_phi[node_id];
- rhs_values[dof_id] = u_g;
- }
- }
- }else if (options()->enforceDirichletMethod() == "WeakPenalty") {
-
- //----------------------------------------------
- // weak penalty method to enforce Dirichlet BC
- //----------------------------------------------
- // Let 'P' be the penalty term and let 'i' be the set of DOF for which
- // Dirichlet condition needs to be applied
- //
- // - For LHS matrix A the diag term corresponding to the Dirichlet DOF
- // a_{i,i} = a_{i,i} + P
- //
- // - For RHS vector b the term that corresponds to the Dirichlet DOF
- // b_{i} = b_{i} * P
- //----------------------------------------------
-
- info() << "Applying Dirichlet boundary condition via "
- << options()->enforceDirichletMethod() << " method ";
-
- Real Penalty = options()->penalty(); // 1.0e30 is the default
-
- ENUMERATE_ (Node, inode, ownNodes()) {
- NodeLocalId node_id = *inode;
- if (m_phi_dirichlet[node_id]) {
- DoFLocalId dof_id = node_dof.dofId(*inode, 0);
- m_linear_system.matrixAddValue(dof_id, dof_id, Penalty);
- Real u_g = Penalty * m_phi[node_id];
- rhs_values[dof_id] = u_g;
- }
- }
- }else if (options()->enforceDirichletMethod() == "RowElimination") {
-
- //----------------------------------------------
- // Row elimination method to enforce Dirichlet BC
- //----------------------------------------------
- // Let 'I' be the set of DOF for which Dirichlet condition needs to be applied
- //
- // to apply the Dirichlet on 'i'th DOF
- // - For LHS matrix A the row terms corresponding to the Dirichlet DOF
- // a_{i,j} = 0. : i!=j
- // a_{i,j} = 1. : i==j
- //----------------------------------------------
-
- info() << "Applying Dirichlet boundary condition via "
- << options()->enforceDirichletMethod() << " method ";
-
- }else if (options()->enforceDirichletMethod() == "RowColumnElimination") {
-
- //----------------------------------------------
- // Row elimination method to enforce Dirichlet BC
- //----------------------------------------------
- // Let 'I' be the set of DOF for which Dirichlet condition needs to be applied
- //
- // to apply the Dirichlet on 'i'th DOF
- // - For LHS matrix A the row terms corresponding to the Dirichlet DOF
- // a_{i,j} = 0. : i!=j for all j
- // a_{i,j} = 1. : i==j
- // also the column terms corresponding to the Dirichlet DOF
- // a_{i,j} = 0. : i!=j for all i
- //----------------------------------------------
-
- info() << "Applying Dirichlet boundary condition via "
- << options()->enforceDirichletMethod() << " method ";
-
- }else {
-
- info() << "Applying Dirichlet boundary condition via "
- << options()->enforceDirichletMethod() << " is not supported \n"
- << "enforce-Dirichlet-method only supports:\n"
- << " - Penalty\n"
- << " - WeakPenalty\n"
- << " - RowElimination\n"
- << " - RowColumnElimination\n";
- }
-
-
- //----------------------------------------------
- // Constant source term assembly
- //----------------------------------------------
- //
- // $int_{Omega}((-rho/epsilon)*v^h)$
- // only for noded that are non-Dirichlet
- //----------------------------------------------
- ENUMERATE_ (Cell, icell, allCells()) {
- Cell cell = *icell;
- Real area = _computeAreaTriangle3(cell);
- for (Node node : cell.nodes()) {
- if (!(m_phi_dirichlet[node]) && node.isOwn())
- rhs_values[node_dof.dofId(node, 0)] += (- rho/epsilon) * area / ElementNodes;
- }
+ if (options()->rho.isPresent()){
+ Real qdot = - rho /epsilon;
+ ArcaneFemFunctions::BoundaryConditions2D::applyConstantSourceToRhs(qdot, mesh(), node_dof, m_node_coord, rhs_values);
}
- //----------------------------------------------
- // Constant flux term assembly
- //----------------------------------------------
- //
- // only for noded that are non-Dirichlet
- // $int_{dOmega_N}((q.n)*v^h)$
- // or
- // $int_{dOmega_N}((n_x*q_x + n_y*q_y)*v^h)$
- //----------------------------------------------
- for (const auto& bs : options()->neumannBoundaryCondition()) {
- FaceGroup group = bs->surface();
-
- if(bs->value.isPresent()) {
- Real value = bs->value();
- ENUMERATE_ (Face, iface, group) {
- Face face = *iface;
- Real length = _computeEdgeLength2(face);
- for (Node node : iface->nodes()) {
- if (!(m_phi_dirichlet[node]) && node.isOwn())
- rhs_values[node_dof.dofId(node, 0)] += value * length / 2.;
- }
- }
- continue;
- }
-
- if(bs->valueX.isPresent() && bs->valueY.isPresent()) {
- Real valueX = bs->valueX();
- Real valueY = bs->valueY();
- ENUMERATE_ (Face, iface, group) {
- Face face = *iface;
- Real length = _computeEdgeLength2(face);
- Real2 Normal = _computeEdgeNormal2(face);
- for (Node node : iface->nodes()) {
- if (!(m_phi_dirichlet[node]) && node.isOwn())
- rhs_values[node_dof.dofId(node, 0)] += (Normal.x*valueX + Normal.y*valueY) * length / 2.;
- }
- }
- continue;
- }
+ for (const auto& bs : options()->neumannBoundaryCondition())
+ ArcaneFemFunctions::BoundaryConditions2D::applyNeumannToRhs(bs, node_dof, m_node_coord, rhs_values);
- if(bs->valueX.isPresent()) {
- Real valueX = bs->valueX();
- ENUMERATE_ (Face, iface, group) {
- Face face = *iface;
- Real length = _computeEdgeLength2(face);
- Real2 Normal = _computeEdgeNormal2(face);
- for (Node node : iface->nodes()) {
- if (!(m_phi_dirichlet[node]) && node.isOwn())
- rhs_values[node_dof.dofId(node, 0)] += (Normal.x*valueX) * length / 2.;
- }
- }
- continue;
- }
-
- if(bs->valueY.isPresent()) {
- Real valueY = bs->valueY();
- ENUMERATE_ (Face, iface, group) {
- Face face = *iface;
- Real length = _computeEdgeLength2(face);
- Real2 Normal = _computeEdgeNormal2(face);
- for (Node node : iface->nodes()) {
- if (!(m_phi_dirichlet[node]) && node.isOwn())
- rhs_values[node_dof.dofId(node, 0)] += (Normal.y*valueY) * length / 2.;
- }
- }
- continue;
- }
-
- }
+ for (const auto& bs : options()->dirichletBoundaryCondition())
+ ArcaneFemFunctions::BoundaryConditions2D::applyDirichletToLhsAndRhs(bs, node_dof, m_node_coord, m_linear_system, rhs_values);
}
/*---------------------------------------------------------------------------*/
@@ -440,7 +146,7 @@ _getE()
ENUMERATE_ (Cell, icell, allCells()) {
Cell cell = *icell;
- Real2 DX = _computeDxDyOfRealTRIA3(cell);
+ Real2 DX = _computeDxDyOfRealTria3(cell);
m_E[cell].x = -DX.x ;
m_E[cell].y = -DX.y ;
m_E[cell].z = 0.0 ;
@@ -452,48 +158,8 @@ _getE()
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
-Real FemModule::
-_computeAreaTriangle3(Cell cell)
-{
- Real3 m0 = m_node_coord[cell.nodeId(0)];
- Real3 m1 = m_node_coord[cell.nodeId(1)];
- Real3 m2 = m_node_coord[cell.nodeId(2)];
- return 0.5 * ((m1.x - m0.x) * (m2.y - m0.y) - (m2.x - m0.x) * (m1.y - m0.y));
-}
-
-/*---------------------------------------------------------------------------*/
-/*---------------------------------------------------------------------------*/
-
-Real FemModule::
-_computeEdgeLength2(Face face)
-{
- Real3 m0 = m_node_coord[face.nodeId(0)];
- Real3 m1 = m_node_coord[face.nodeId(1)];
- return math::sqrt((m1.x-m0.x)*(m1.x-m0.x) + (m1.y-m0.y)*(m1.y - m0.y));
-}
-
-/*---------------------------------------------------------------------------*/
-/*---------------------------------------------------------------------------*/
-
-Real2 FemModule::
-_computeEdgeNormal2(Face face)
-{
- Real3 m0 = m_node_coord[face.nodeId(0)];
- Real3 m1 = m_node_coord[face.nodeId(1)];
- if (!face.isSubDomainBoundaryOutside())
- std::swap(m0,m1);
- Real2 N;
- Real norm_N = math::sqrt( (m1.y - m0.y)*(m1.y - m0.y) + (m1.x - m0.x)*(m1.x - m0.x) ); // for normalizing
- N.x = (m1.y - m0.y)/ norm_N;
- N.y = (m0.x - m1.x)/ norm_N;
- return N;
-}
-
-/*---------------------------------------------------------------------------*/
-/*---------------------------------------------------------------------------*/
-
Real2 FemModule::
-_computeDxDyOfRealTRIA3(Cell cell)
+_computeDxDyOfRealTria3(Cell cell)
{
Real3 m0 = m_node_coord[cell.nodeId(0)];
Real3 m1 = m_node_coord[cell.nodeId(1)];
@@ -513,80 +179,78 @@ _computeDxDyOfRealTRIA3(Cell cell)
}
/*---------------------------------------------------------------------------*/
+/**
+ * @brief Computes the element matrix for a triangular element (P1 FE).
+ *
+ * This function calculates the integral of the expression:
+ * (u.dx * v.dx + u.dy * v.dy)
+ *
+ * Steps involved:
+ * 1. Calculate the area of the triangle.
+ * 2. Compute the gradients of the shape functions.
+ * 3. Return (u.dx * v.dx + u.dy * v.dy);
+ */
/*---------------------------------------------------------------------------*/
FixedMatrix<3, 3> FemModule::
-_computeElementMatrixTRIA3(Cell cell)
+_computeElementMatrixTria3(Cell cell)
{
- // Get coordiantes of the triangle element TRI3
- //------------------------------------------------
- // 0 o
- // . .
- // . .
- // . .
- // 1 o . . . o 2
- //------------------------------------------------
- Real3 m0 = m_node_coord[cell.nodeId(0)];
- Real3 m1 = m_node_coord[cell.nodeId(1)];
- Real3 m2 = m_node_coord[cell.nodeId(2)];
-
- Real area = _computeAreaTriangle3(cell); // calculate area
-
- Real2 dPhi0(m1.y - m2.y, m2.x - m1.x);
- Real2 dPhi1(m2.y - m0.y, m0.x - m2.x);
- Real2 dPhi2(m0.y - m1.y, m1.x - m0.x);
+ Real area = ArcaneFemFunctions::MeshOperation::computeAreaTria3(cell, m_node_coord);
- FixedMatrix<2, 3> b_matrix;
- b_matrix(0, 0) = dPhi0.x;
- b_matrix(0, 1) = dPhi1.x;
- b_matrix(0, 2) = dPhi2.x;
+ Real3 dxU = ArcaneFemFunctions::FeOperation2D::computeGradientXTria3(cell, m_node_coord);
+ Real3 dyU = ArcaneFemFunctions::FeOperation2D::computeGradientYTria3(cell, m_node_coord);
- b_matrix(1, 0) = dPhi0.y;
- b_matrix(1, 1) = dPhi1.y;
- b_matrix(1, 2) = dPhi2.y;
-
- b_matrix.multInPlace(1.0 / (2.0 * area));
+ return area * (dxU ^ dxU) + area * (dyU ^ dyU);
+}
- FixedMatrix<3, 3> int_cdPi_dPj = matrixMultiplication(matrixTranspose(b_matrix), b_matrix);
- int_cdPi_dPj.multInPlace(area);
- //info() << "Cell=" << cell.localId();
- //std::cout << " int_cdPi_dPj=";
- //int_cdPi_dPj.dump(std::cout);
- //std::cout << "\n";
+/*---------------------------------------------------------------------------*/
+/**
+ * @brief Calls the right function for LHS assembly given as mesh type.
+ */
+/*---------------------------------------------------------------------------*/
- return int_cdPi_dPj;
+void FemModule::
+_assembleBilinearOperator()
+{
+ if (options()->meshType == "QUAD4")
+ ARCANE_FATAL("Non supported meshType");
+ else if (options()->meshType == "TRIA3")
+ _assembleBilinear<3>([this](const Cell& cell) {
+ return _computeElementMatrixTria3(cell);
+ });
+ else
+ ARCANE_FATAL("Non supported meshType");
}
+
/*---------------------------------------------------------------------------*/
+/**
+ * @brief Assembles the bilinear operator matrix for the FEM linear system.
+ *
+ * The method performs the following steps:
+ * 1. For each cell, retrieves the cell-specific constant `lambda`.
+ * 2. Computes element matrix using provided `compute_element_matrix` function.
+ * 3. Assembles global matrix by adding contributions from each cell's element
+ * matrix to the corresponding entries in the global matrix.
+ */
/*---------------------------------------------------------------------------*/
+template
void FemModule::
-_assembleBilinearOperatorTRIA3()
+_assembleBilinear(const std::function(const Cell&)>& compute_element_matrix)
{
auto node_dof(m_dofs_on_nodes.nodeDoFConnectivityView());
ENUMERATE_ (Cell, icell, allCells()) {
Cell cell = *icell;
- if (cell.type() != IT_Triangle3)
- ARCANE_FATAL("Only Triangle3 cell type is supported");
-
- auto K_e = _computeElementMatrixTRIA3(cell); // element stifness matrix
- // # assemble elementary matrix into the global one
- // # elementary terms are positionned into K according
- // # to the rank of associated node in the mesh.nodes list
- // for node1 in elem.nodes:
- // inode1=elem.nodes.index(node1) # get position of node1 in nodes list
- // for node2 in elem.nodes:
- // inode2=elem.nodes.index(node2)
- // K[node1.rank,node2.rank]=K[node1.rank,node2.rank]+K_e[inode1,inode2]
+
+ 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()) {
- // K[node1.rank,node2.rank]=K[node1.rank,node2.rank]+K_e[inode1,inode2]
Real v = K_e(n1_index, n2_index);
- // m_k_matrix(node1.localId(), node2.localId()) += v;
if (node1.isOwn()) {
m_linear_system.matrixAddValue(node_dof.dofId(node1, 0), node_dof.dofId(node2, 0), v);
}
@@ -598,6 +262,14 @@ _assembleBilinearOperatorTRIA3()
}
/*---------------------------------------------------------------------------*/
+/**
+ * @brief Solves the linear system and updates the solution vector.
+ *
+ * This method performs the following actions:
+ * 1. Solves the linear system to compute the solution.
+ * 2. Copies the computed solution from the DoF to the node values.
+ * 3. Synchronizes the updated node values.
+ */
/*---------------------------------------------------------------------------*/
void FemModule::
@@ -605,10 +277,6 @@ _solve()
{
m_linear_system.solve();
- // Re-Apply boundary conditions because the solver has modified the value
- // of u on all nodes
- _applyDirichletBoundaryConditions();
-
{
VariableDoFReal& dof_u(m_linear_system.solutionVariable());
// Copy RHS DoF to Node u
@@ -621,31 +289,35 @@ _solve()
}
m_phi.synchronize();
-
- const bool do_print = (allNodes().size() < 200);
- if (do_print) {
- ENUMERATE_ (Node, inode, allNodes()) {
- Node node = *inode;
- info() << "Phi[" << node.localId() << "][" << node.uniqueId() << "] = "
- << m_phi[node];
- //info() << "Phi[]" << node.uniqueId() << " "
- // << m_phi[node];
- }
- }
}
/*---------------------------------------------------------------------------*/
+/**
+ * @brief Validates and prints the results of the FEM computation.
+ *
+ * This method performs the following actions:
+ * 1. If number of nodes < 200, prints the computed values for each node.
+ * 2. Retrieves the filename for the result file from options.
+ * 3. If a filename is provided, checks the computed results against result file.
+ *
+ * @note The result comparison uses a tolerance of 1.0e-4.
+ */
/*---------------------------------------------------------------------------*/
void FemModule::
-_checkResultFile()
+_validateResults()
{
+ if (allNodes().size() < 200)
+ ENUMERATE_ (Node, inode, allNodes()) {
+ Node node = *inode;
+ info() << "phi[" << node.localId() << "][" << node.uniqueId() << "] = " << m_phi[node];
+ }
+
String filename = options()->resultFile();
- info() << "CheckResultFile filename=" << filename;
- if (filename.empty())
- return;
- const double epsilon = 1.0e-4;
- checkNodeResultFile(traceMng(), filename, m_phi, epsilon);
+ info() << "ValidateResultFile filename=" << filename;
+
+ if (!filename.empty())
+ checkNodeResultFile(traceMng(), filename, m_phi, 1.0e-4);
}
/*---------------------------------------------------------------------------*/
diff --git a/electrostatics/FemModule.h b/electrostatics/FemModule.h
new file mode 100644
index 00000000..09da1d0d
--- /dev/null
+++ b/electrostatics/FemModule.h
@@ -0,0 +1,99 @@
+// -*- tab-width: 2; indent-tabs-mode: nil; coding: utf-8-with-signature -*-
+//-----------------------------------------------------------------------------
+// Copyright 2000-2024 CEA (www.cea.fr) IFPEN (www.ifpenergiesnouvelles.com)
+// See the top-level COPYRIGHT file for details.
+// SPDX-License-Identifier: Apache-2.0
+//-----------------------------------------------------------------------------
+/*---------------------------------------------------------------------------*/
+/* FemModule.cc (C) 2022-2024 */
+/* */
+/* FemModule class definition. */
+/*---------------------------------------------------------------------------*/
+/*---------------------------------------------------------------------------*/
+#ifndef FEMMODULES_H
+#define FEMMODULES_H
+
+#include
+#include
+#include
+
+#include
+#include
+#include
+#include
+#include
+
+#include "IDoFLinearSystemFactory.h"
+#include "Fem_axl.h"
+#include "FemUtils.h"
+#include "DoFLinearSystem.h"
+#include "FemDoFsOnNodes.h"
+#include "ArcaneFemFunctions.h"
+
+/*---------------------------------------------------------------------------*/
+/*---------------------------------------------------------------------------*/
+
+using namespace Arcane;
+using namespace Arcane::FemUtils;
+
+/*---------------------------------------------------------------------------*/
+/*---------------------------------------------------------------------------*/
+/*!
+ * \brief Module Fem.
+ */
+class FemModule
+: public ArcaneFemObject
+{
+ public:
+
+ explicit FemModule(const ModuleBuildInfo& mbi)
+ : ArcaneFemObject(mbi)
+ , m_dofs_on_nodes(mbi.subDomain()->traceMng())
+ {
+ ICaseMng* cm = mbi.subDomain()->caseMng();
+ cm->setTreatWarningAsError(true);
+ cm->setAllowUnkownRootElelement(false);
+ }
+
+ public:
+
+ //! Method called at each iteration
+ void compute() override;
+
+ //! Method called at the beginning of the simulation
+ void startInit() override;
+
+ VersionInfo versionInfo() const override
+ {
+ return VersionInfo(1, 0, 0);
+ }
+
+ private:
+
+ Real rho;
+ Real epsilon;
+ Real ElementNodes;
+
+ DoFLinearSystem m_linear_system;
+ IItemFamily* m_dof_family = nullptr;
+ FemDoFsOnNodes m_dofs_on_nodes;
+
+ private:
+
+ void _doStationarySolve();
+ void _getMaterialParameters();
+ void _assembleBilinearOperator();
+ void _solve();
+ void _getE();
+ void _assembleLinearOperator();
+ void _validateResults();
+
+ FixedMatrix<3, 3> _computeElementMatrixTria3(Cell cell);
+
+ Real2 _computeDxDyOfRealTria3(Cell cell);
+
+ template
+ void _assembleBilinear( const std::function(const Cell&)>& compute_element_matrix);
+};
+
+#endif
\ No newline at end of file
diff --git a/fourier/FemModule.cc b/fourier/FemModule.cc
index 1fd63b41..01cc540c 100644
--- a/fourier/FemModule.cc
+++ b/fourier/FemModule.cc
@@ -121,9 +121,7 @@ _doStationarySolve()
* This method initializes:
* - material properties:
* # thermal conductivity coefficient (`lambda`)
- * # heat source term (`qdot`)
- * - mesh properties:
- * # number of cell nodes per element based on the mesh (quad or tria)
+ * # heat source term (`qdot`)
*/
/*---------------------------------------------------------------------------*/