From 406705e8190e699a830734b35ba1c5a03cfc87d3 Mon Sep 17 00:00:00 2001 From: Mohd Afeef Badri <52162083+mohd-afeef-badri@users.noreply.github.com> Date: Sat, 3 Aug 2024 20:46:54 +0200 Subject: [PATCH 01/32] ArcaneFem functions - start moving functions to femutils --- femutils/ArcaneFemFunctions.cc | 63 ++++++++++++++++++++++++++++++++++ femutils/ArcaneFemFunctions.h | 16 +++++++++ femutils/CMakeLists.txt | 2 ++ 3 files changed, 81 insertions(+) create mode 100644 femutils/ArcaneFemFunctions.cc create mode 100644 femutils/ArcaneFemFunctions.h diff --git a/femutils/ArcaneFemFunctions.cc b/femutils/ArcaneFemFunctions.cc new file mode 100644 index 0000000..5e79975 --- /dev/null +++ b/femutils/ArcaneFemFunctions.cc @@ -0,0 +1,63 @@ +#include "ArcaneFemFunctions.h" + +/*---------------------------------------------------------------------------*/ +/** + * @brief Computes the area of a triangle defined by three nodes. + * + * This method calculates the area using the determinant formula for a triangle. + * The area is computed as half the value of the determinant of the matrix + * formed by the coordinates of the triangle's vertices. + */ +/*---------------------------------------------------------------------------*/ + +Real ArcaneFemFunctions::computeAreaTriangle3(Cell cell, VariableNodeReal3& node_coord) +{ + Real3 m0 = node_coord[cell.nodeId(0)]; + Real3 m1 = node_coord[cell.nodeId(1)]; + Real3 m2 = node_coord[cell.nodeId(2)]; + + return 0.5 * ((m1.x - m0.x) * (m2.y - m0.y) - (m2.x - m0.x) * (m1.y - m0.y)); +} + +/*---------------------------------------------------------------------------*/ +/** + * @brief Computes the length of the edge defined by a given face. + * + * This method calculates Euclidean distance between the two nodes of the face. + */ +/*---------------------------------------------------------------------------*/ + +Real ArcaneFemFunctions::computeEdgeLength2(Face face, VariableNodeReal3& node_coord) +{ + Real3 m0 = node_coord[face.nodeId(0)]; + Real3 m1 = node_coord[face.nodeId(1)]; + + Real dx = m1.x - m0.x; + Real dy = m1.y - m0.y; + + return math::sqrt(dx * dx + dy * dy); +} + +/*---------------------------------------------------------------------------*/ +/** + * @brief Computes the normalized edge normal for a given face. + * + * This method calculates normal vector to the edge defined by nodes of the face, + * normalizes it, and ensures the correct orientation. + */ +/*---------------------------------------------------------------------------*/ + +Real2 ArcaneFemFunctions::computeEdgeNormal2(Face face, VariableNodeReal3& node_coord) +{ + Real3 m0 = node_coord[face.nodeId(0)]; + Real3 m1 = node_coord[face.nodeId(1)]; + + if (!face.isSubDomainBoundaryOutside()) + std::swap(m0, m1); + + Real dx = m1.x - m0.x; + Real dy = m1.y - m0.y; + Real norm_N = math::sqrt(dx * dx + dy * dy); + + return { dy / norm_N, -dx / norm_N }; +} diff --git a/femutils/ArcaneFemFunctions.h b/femutils/ArcaneFemFunctions.h new file mode 100644 index 0000000..0e80ac5 --- /dev/null +++ b/femutils/ArcaneFemFunctions.h @@ -0,0 +1,16 @@ +#ifndef ARCANE_FEM_FUNCTIONS_H +#define ARCANE_FEM_FUNCTIONS_H + +#include + +using namespace Arcane; + +class ArcaneFemFunctions +{ +public: + static Real computeAreaTriangle3(Cell cell, VariableNodeReal3& node_coord); + static Real computeEdgeLength2(Face face, VariableNodeReal3& node_coord); + static Real2 computeEdgeNormal2(Face face, VariableNodeReal3& node_coord); +}; + +#endif // ARCANE_FEM_FUNCTIONS_H diff --git a/femutils/CMakeLists.txt b/femutils/CMakeLists.txt index 10ee83a..d80db50 100644 --- a/femutils/CMakeLists.txt +++ b/femutils/CMakeLists.txt @@ -10,6 +10,8 @@ add_library(FemUtils CsrFormatMatrix.cc FemDoFsOnNodes.h FemDoFsOnNodes.cc + ArcaneFemFunctions.h + ArcaneFemFunctions.cc AlephDoFLinearSystem.cc IDoFLinearSystemFactory.h AlephDoFLinearSystemFactory_axl.h From 7ae564fb1c3d0075aaaea3010788fc0f72673831 Mon Sep 17 00:00:00 2001 From: mohd-afeef-badri Date: Sun, 4 Aug 2024 18:18:12 +0200 Subject: [PATCH 02/32] neglect lower values in test #148 --- elasticity/FemModule.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/elasticity/FemModule.cc b/elasticity/FemModule.cc index 88b9c59..96182ad 100644 --- a/elasticity/FemModule.cc +++ b/elasticity/FemModule.cc @@ -956,7 +956,7 @@ _checkResultFile() if (filename.empty()) return; const double epsilon = 1.0e-4; - const double min_value_to_test = 1.0e-18; + const double min_value_to_test = 1.0e-16; Arcane::FemUtils::checkNodeResultFile(traceMng(), filename, m_U, epsilon, min_value_to_test); } From 7bba22c131cef826b6a7c8c66365dbbf88ef278f Mon Sep 17 00:00:00 2001 From: Mohd Afeef Badri <52162083+mohd-afeef-badri@users.noreply.github.com> Date: Mon, 5 Aug 2024 10:32:57 +0200 Subject: [PATCH 03/32] ArcaneFemFunctions Move common operations to ArcaneFEMFunctions for code refactoring. - New Real3x3 function for gradient --- acoustics/FemModule.cc | 132 +++++++-------------------------- femutils/ArcaneFemFunctions.cc | 63 ---------------- femutils/ArcaneFemFunctions.h | 88 +++++++++++++++++++++- femutils/CMakeLists.txt | 1 - 4 files changed, 112 insertions(+), 172 deletions(-) delete mode 100644 femutils/ArcaneFemFunctions.cc diff --git a/acoustics/FemModule.cc b/acoustics/FemModule.cc index 2dde217..7ffac5d 100644 --- a/acoustics/FemModule.cc +++ b/acoustics/FemModule.cc @@ -26,6 +26,7 @@ #include "FemUtils.h" #include "DoFLinearSystem.h" #include "FemDoFsOnNodes.h" +#include "ArcaneFemFunctions.h" /*---------------------------------------------------------------------------*/ /*---------------------------------------------------------------------------*/ @@ -57,7 +58,7 @@ class FemModule cm->setAllowUnkownRootElelement(false); } - void compute() override; //! Method called at each iteration + void compute() override; //! Method called at each iteration void startInit() override; //! Method called at the beginning of the simulation VersionInfo versionInfo() const override { return VersionInfo(1, 0, 0); } @@ -77,10 +78,6 @@ class FemModule void _validateResults(); FixedMatrix<3, 3> _computeElementMatrixTRIA3(Cell cell); - - Real _computeAreaTriangle3(Cell cell); - Real _computeEdgeLength2(Face face); - Real2 _computeEdgeNormal2(Face face); }; /*---------------------------------------------------------------------------*/ @@ -215,8 +212,8 @@ _assembleLinearOperator() Face face = *iface; // step 3. - Real length = _computeEdgeLength2(face); - Real2 normal = _computeEdgeNormal2(face); + Real length = ArcaneFemFunctions::computeEdgeLength2(face, m_node_coord); + Real2 normal = ArcaneFemFunctions::computeEdgeNormal2(face, m_node_coord); // step 4 for (Node node : iface->nodes()) { @@ -238,70 +235,6 @@ _assembleLinearOperator() } } -/*---------------------------------------------------------------------------*/ -/** - * @brief Computes the area of a triangle defined by three nodes. - * - * This method calculates the area using the determinant formula for a triangle. - * The area is computed as half the value of the determinant of the matrix - * formed by the coordinates of the triangle's vertices. - */ -/*---------------------------------------------------------------------------*/ - -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)); -} - -/*---------------------------------------------------------------------------*/ -/** - * @brief Computes the length of the edge defined by a given face. - * - * This method calculates Euclidean distance between the two nodes of the face. - */ -/*---------------------------------------------------------------------------*/ - -Real FemModule:: -_computeEdgeLength2(Face face) -{ - Real3 m0 = m_node_coord[face.nodeId(0)]; - Real3 m1 = m_node_coord[face.nodeId(1)]; - - Real dx = m1.x - m0.x; - Real dy = m1.y - m0.y; - - return math::sqrt(dx * dx + dy * dy); -} - -/*---------------------------------------------------------------------------*/ -/** - * @brief Computes the normalized edge normal for a given face. - * - * This method calculates normal vector to the edge defined by nodes of the face, - * normalizes it, and ensures the correct orientation. - */ -/*---------------------------------------------------------------------------*/ - -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); - - Real dx = m1.x - m0.x; - Real dy = m1.y - m0.y; - Real norm_N = math::sqrt(dx * dx + dy * dy); - - return { dy / norm_N, -dx / norm_N }; -} - /*---------------------------------------------------------------------------*/ /** * @brief Computes the element matrix for a triangular element (P1 FE). @@ -310,55 +243,46 @@ _computeEdgeNormal2(Face face) * -(u.dx * v.dx + u.dy * v.dy) + kc2 * u * v * * Steps involved: - * 1. Get the coordinates of the triangle vertices. - * 2. Calculate the area of the triangle. - * 3. Compute the gradients of the shape functions. - * 4. Normalize the gradients by the area. - * 5. Compute the element matrix: - * 5.1 first compute -(u.dx * v.dx + u.dy * v.dy) terms - * 5.2 then (kc2 * u * v) term adjusted + * 1. Calculate the area of the triangle. + * 2. Compute the gradients of the shape functions. + * 3. Normalize the gradients by the area. + * 4. Compute the element matrix: + * 4.1 first compute -(u.dx * v.dx + u.dy * v.dy) terms + * 4.2 then (kc2 * u * v) term adjusted */ /*---------------------------------------------------------------------------*/ FixedMatrix<3, 3> FemModule:: _computeElementMatrixTRIA3(Cell cell) { + // step 1 - Real3 vertex0 = m_node_coord[cell.nodeId(0)]; - Real3 vertex1 = m_node_coord[cell.nodeId(1)]; - Real3 vertex2 = m_node_coord[cell.nodeId(2)]; + Real area = ArcaneFemFunctions::computeAreaTriangle3(cell, m_node_coord); // step 2 - Real area = _computeAreaTriangle3(cell); + Real3x3 gradN = ArcaneFemFunctions::computeGradientTria3(cell, m_node_coord); // step 3 - Real2 gradN0(vertex1.y - vertex2.y, vertex2.x - vertex1.x); - Real2 gradN1(vertex2.y - vertex0.y, vertex0.x - vertex2.x); - Real2 gradN2(vertex0.y - vertex1.y, vertex1.x - vertex0.x); - - // step 4 Real A2 = 2.0 * area; - gradN0 /= A2; - gradN1 /= A2; - gradN2 /= A2; + gradN /= A2; - // step 5 + // step 4 FixedMatrix<3, 3> elementMatrix; - // step 5.1 - elementMatrix(0, 0) = -area * Arcane::math::dot(gradN0, gradN0); - elementMatrix(0, 1) = -area * Arcane::math::dot(gradN0, gradN1); - elementMatrix(0, 2) = -area * Arcane::math::dot(gradN0, gradN2); + // step 4.1 + elementMatrix(0, 0) = -area * Arcane::math::dot(gradN[0], gradN[0]); + elementMatrix(0, 1) = -area * Arcane::math::dot(gradN[0], gradN[1]); + elementMatrix(0, 2) = -area * Arcane::math::dot(gradN[0], gradN[2]); - elementMatrix(1, 0) = -area * Arcane::math::dot(gradN1, gradN0); - elementMatrix(1, 1) = -area * Arcane::math::dot(gradN1, gradN1); - elementMatrix(1, 2) = -area * Arcane::math::dot(gradN1, gradN2); + elementMatrix(1, 0) = -area * Arcane::math::dot(gradN[1], gradN[0]); + elementMatrix(1, 1) = -area * Arcane::math::dot(gradN[1], gradN[1]); + elementMatrix(1, 2) = -area * Arcane::math::dot(gradN[1], gradN[2]); - elementMatrix(2, 0) = -area * Arcane::math::dot(gradN2, gradN0); - elementMatrix(2, 1) = -area * Arcane::math::dot(gradN2, gradN1); - elementMatrix(2, 2) = -area * Arcane::math::dot(gradN2, gradN2); + elementMatrix(2, 0) = -area * Arcane::math::dot(gradN[2], gradN[0]); + elementMatrix(2, 1) = -area * Arcane::math::dot(gradN[2], gradN[1]); + elementMatrix(2, 2) = -area * Arcane::math::dot(gradN[2], gradN[2]); - // step 5 .2 + // step 4.2 Real uvTerm = m_kc2 * area / 12.0; elementMatrix(0, 0) += uvTerm * 2.; elementMatrix(0, 1) += uvTerm; @@ -442,12 +366,12 @@ _solve() /*---------------------------------------------------------------------------*/ /** * @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. */ /*---------------------------------------------------------------------------*/ diff --git a/femutils/ArcaneFemFunctions.cc b/femutils/ArcaneFemFunctions.cc deleted file mode 100644 index 5e79975..0000000 --- a/femutils/ArcaneFemFunctions.cc +++ /dev/null @@ -1,63 +0,0 @@ -#include "ArcaneFemFunctions.h" - -/*---------------------------------------------------------------------------*/ -/** - * @brief Computes the area of a triangle defined by three nodes. - * - * This method calculates the area using the determinant formula for a triangle. - * The area is computed as half the value of the determinant of the matrix - * formed by the coordinates of the triangle's vertices. - */ -/*---------------------------------------------------------------------------*/ - -Real ArcaneFemFunctions::computeAreaTriangle3(Cell cell, VariableNodeReal3& node_coord) -{ - Real3 m0 = node_coord[cell.nodeId(0)]; - Real3 m1 = node_coord[cell.nodeId(1)]; - Real3 m2 = node_coord[cell.nodeId(2)]; - - return 0.5 * ((m1.x - m0.x) * (m2.y - m0.y) - (m2.x - m0.x) * (m1.y - m0.y)); -} - -/*---------------------------------------------------------------------------*/ -/** - * @brief Computes the length of the edge defined by a given face. - * - * This method calculates Euclidean distance between the two nodes of the face. - */ -/*---------------------------------------------------------------------------*/ - -Real ArcaneFemFunctions::computeEdgeLength2(Face face, VariableNodeReal3& node_coord) -{ - Real3 m0 = node_coord[face.nodeId(0)]; - Real3 m1 = node_coord[face.nodeId(1)]; - - Real dx = m1.x - m0.x; - Real dy = m1.y - m0.y; - - return math::sqrt(dx * dx + dy * dy); -} - -/*---------------------------------------------------------------------------*/ -/** - * @brief Computes the normalized edge normal for a given face. - * - * This method calculates normal vector to the edge defined by nodes of the face, - * normalizes it, and ensures the correct orientation. - */ -/*---------------------------------------------------------------------------*/ - -Real2 ArcaneFemFunctions::computeEdgeNormal2(Face face, VariableNodeReal3& node_coord) -{ - Real3 m0 = node_coord[face.nodeId(0)]; - Real3 m1 = node_coord[face.nodeId(1)]; - - if (!face.isSubDomainBoundaryOutside()) - std::swap(m0, m1); - - Real dx = m1.x - m0.x; - Real dy = m1.y - m0.y; - Real norm_N = math::sqrt(dx * dx + dy * dy); - - return { dy / norm_N, -dx / norm_N }; -} diff --git a/femutils/ArcaneFemFunctions.h b/femutils/ArcaneFemFunctions.h index 0e80ac5..26bec01 100644 --- a/femutils/ArcaneFemFunctions.h +++ b/femutils/ArcaneFemFunctions.h @@ -7,10 +7,90 @@ using namespace Arcane; class ArcaneFemFunctions { -public: - static Real computeAreaTriangle3(Cell cell, VariableNodeReal3& node_coord); - static Real computeEdgeLength2(Face face, VariableNodeReal3& node_coord); - static Real2 computeEdgeNormal2(Face face, VariableNodeReal3& node_coord); + public: + + /*---------------------------------------------------------------------------*/ + /** + * @brief Computes the area of a triangle defined by three nodes. + * + * This method calculates the area using the determinant formula for a triangle. + * The area is computed as half the value of the determinant of the matrix + * formed by the coordinates of the triangle's vertices. + */ + /*---------------------------------------------------------------------------*/ + + static inline Real computeAreaTriangle3(Cell cell, VariableNodeReal3& node_coord) + { + Real3 vertex0 = node_coord[cell.nodeId(0)]; + Real3 vertex1 = node_coord[cell.nodeId(1)]; + Real3 vertex2 = node_coord[cell.nodeId(2)]; + + return 0.5 * ((vertex1.x - vertex0.x) * (vertex2.y - vertex0.y) - (vertex2.x - vertex0.x) * (vertex1.y - vertex0.y)); + } + + /*---------------------------------------------------------------------------*/ + /** + * @brief Computes the length of the edge defined by a given face. + * + * This method calculates Euclidean distance between the two nodes of the face. + */ + /*---------------------------------------------------------------------------*/ + + static inline Real computeEdgeLength2(Face face, VariableNodeReal3& node_coord) + { + Real3 vertex0 = node_coord[face.nodeId(0)]; + Real3 vertex1 = node_coord[face.nodeId(1)]; + + Real dx = vertex1.x - vertex0.x; + Real dy = vertex1.y - vertex0.y; + + return math::sqrt(dx * dx + dy * dy); + } + + /*---------------------------------------------------------------------------*/ + /** + * @brief Computes the normalized edge normal for a given face. + * + * This method calculates normal vector to the edge defined by nodes of the face, + * normalizes it, and ensures the correct orientation. + */ + /*---------------------------------------------------------------------------*/ + + static inline Real2 computeEdgeNormal2(Face face, VariableNodeReal3& node_coord) + { + Real3 vertex0 = node_coord[face.nodeId(0)]; + Real3 vertex1 = node_coord[face.nodeId(1)]; + + if (!face.isSubDomainBoundaryOutside()) + std::swap(vertex0, vertex1); + + Real dx = vertex1.x - vertex0.x; + Real dy = vertex1.y - vertex0.y; + Real norm_N = math::sqrt(dx * dx + dy * dy); + + return { dy / norm_N, -dx / norm_N }; + } + + /*---------------------------------------------------------------------------*/ + /** + * @brief Computes the gradients of basis functions N for P1 triangles. + * + * This method calculates gradient operator for a given P1 cell + * + * @note we are using 3x3 matrix for this even in 2D hence the z gradient is 0. + */ + /*---------------------------------------------------------------------------*/ + + static inline Real3x3 computeGradientTria3(Cell cell, VariableNodeReal3& node_coord) + { + Real3 vertex0 = node_coord[cell.nodeId(0)]; + Real3 vertex1 = node_coord[cell.nodeId(1)]; + Real3 vertex2 = node_coord[cell.nodeId(2)]; + + return Real3x3(Real3(vertex1.y - vertex2.y, vertex2.x - vertex1.x, 0), + Real3(vertex2.y - vertex0.y, vertex0.x - vertex2.x, 0), + Real3(vertex0.y - vertex1.y, vertex1.x - vertex0.x, 0)); + } }; #endif // ARCANE_FEM_FUNCTIONS_H diff --git a/femutils/CMakeLists.txt b/femutils/CMakeLists.txt index d80db50..3376abd 100644 --- a/femutils/CMakeLists.txt +++ b/femutils/CMakeLists.txt @@ -11,7 +11,6 @@ add_library(FemUtils FemDoFsOnNodes.h FemDoFsOnNodes.cc ArcaneFemFunctions.h - ArcaneFemFunctions.cc AlephDoFLinearSystem.cc IDoFLinearSystemFactory.h AlephDoFLinearSystemFactory_axl.h From 81f038b8e8d64ae74c8a98743c6021c11e91284b Mon Sep 17 00:00:00 2001 From: Mohd Afeef Badri <52162083+mohd-afeef-badri@users.noreply.github.com> Date: Mon, 5 Aug 2024 14:01:48 +0200 Subject: [PATCH 04/32] new operations for FixedMarix add operator + subtraction - unary negation operator --- femutils/FemUtils.h | 36 ++++++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) diff --git a/femutils/FemUtils.h b/femutils/FemUtils.h index 79793cc..fbdb2eb 100644 --- a/femutils/FemUtils.h +++ b/femutils/FemUtils.h @@ -84,6 +84,42 @@ class FixedMatrix } } + //! Define the addition operator + FixedMatrix operator+(const FixedMatrix& other) const + { + FixedMatrix result; + for (Arcane::Int32 i = 0; i < N; ++i) { + for (Arcane::Int32 j = 0; j < M; ++j) { + result(i, j) = (*this)(i, j) + other(i, j); + } + } + return result; + } + + //! Define the subtraction operator + FixedMatrix operator-(const FixedMatrix& other) const + { + FixedMatrix result; + for (Arcane::Int32 i = 0; i < N; ++i) { + for (Arcane::Int32 j = 0; j < M; ++j) { + result(i, j) = (*this)(i, j) - other(i, j); + } + } + return result; + } + + //! Define the unary negation operator + FixedMatrix operator-() const + { + FixedMatrix result; + for (Arcane::Int32 i = 0; i < N; ++i) { + for (Arcane::Int32 j = 0; j < M; ++j) { + result(i, j) = -(*this)(i, j); + } + } + return result; + } + private: std::array m_values = {}; From 6f08b0eaede8faeff19f791ebb332dfd7c785a93 Mon Sep 17 00:00:00 2001 From: Mohd Afeef Badri <52162083+mohd-afeef-badri@users.noreply.github.com> Date: Mon, 5 Aug 2024 14:24:42 +0200 Subject: [PATCH 05/32] new operations for FixedMarix scalar multiplation outer vector product --- femutils/FemUtils.h | 40 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) diff --git a/femutils/FemUtils.h b/femutils/FemUtils.h index fbdb2eb..efe01c3 100644 --- a/femutils/FemUtils.h +++ b/femutils/FemUtils.h @@ -22,6 +22,8 @@ #include #include +#include + #include /*---------------------------------------------------------------------------*/ @@ -120,11 +122,49 @@ class FixedMatrix return result; } + //! Scalar multiplication: FixedMatrix * scalar + FixedMatrix operator*(Real scalar) const + { + FixedMatrix result; + for (Arcane::Int32 i = 0; i < N; ++i) { + for (Arcane::Int32 j = 0; j < M; ++j) { + result(i, j) = (*this)(i, j) * scalar; + } + } + return result; + } + + //! Friend function for scalar multiplication: scalar * FixedMatrix + friend FixedMatrix operator*(Real scalar, const FixedMatrix& matrix) + { + FixedMatrix result; + for (Arcane::Int32 i = 0; i < N; ++i) { + for (Arcane::Int32 j = 0; j < M; ++j) { + result(i, j) = scalar * matrix(i, j); + } + } + return result; + } + private: std::array m_values = {}; }; +/*---------------------------------------------------------------------------*/ +// Outer product of two Real3 vectors to produce a FixedMatrix<3, 3> +/*---------------------------------------------------------------------------*/ +inline FixedMatrix<3, 3> operator^(const Arcane::Real3& lhs, const Arcane::Real3& rhs) +{ + FixedMatrix<3, 3> result; + for (Arcane::Int32 i = 0; i < 3; ++i) { + for (Arcane::Int32 j = 0; j < 3; ++j) { + result(i, j) = lhs[i] * rhs[j]; + } + } + return result; +} + /*---------------------------------------------------------------------------*/ /*---------------------------------------------------------------------------*/ From 78ba56818d5696493fbc91964ea8d35d13b7cd3e Mon Sep 17 00:00:00 2001 From: Mohd Afeef Badri <52162083+mohd-afeef-badri@users.noreply.github.com> Date: Mon, 5 Aug 2024 16:36:21 +0200 Subject: [PATCH 06/32] calculate full graident dx and dy --- femutils/ArcaneFemFunctions.h | 22 +++++++++++++++++----- 1 file changed, 17 insertions(+), 5 deletions(-) diff --git a/femutils/ArcaneFemFunctions.h b/femutils/ArcaneFemFunctions.h index 26bec01..c7cd6ed 100644 --- a/femutils/ArcaneFemFunctions.h +++ b/femutils/ArcaneFemFunctions.h @@ -75,9 +75,19 @@ class ArcaneFemFunctions /** * @brief Computes the gradients of basis functions N for P1 triangles. * - * This method calculates gradient operator for a given P1 cell + * This method calculates gradient operator ∇ Ni for a given P1 cell + * with i = 1,..,3 for the three shape function Ni hence output is a + * matrix of size 3x3 * - * @note we are using 3x3 matrix for this even in 2D hence the z gradient is 0. + * [ ∂N1/∂x ∂N1/∂y ∂N1/∂z ] + * ∇N = [ ∂N2/∂x ∂N2/∂y ∂N2/∂z ] + * [ ∂N3/∂x ∂N3/∂y ∂N3/∂z ] + * + * [ y2​−y3 x3​−x2 0 ] + * ∇N = 1/(2A) [ y3−y1 x1−x3 0 ] + * [ y1−y2 x2−x1 0 ] + * + * @note we can adapt the same for 3D by filling the third component */ /*---------------------------------------------------------------------------*/ @@ -87,9 +97,11 @@ class ArcaneFemFunctions Real3 vertex1 = node_coord[cell.nodeId(1)]; Real3 vertex2 = node_coord[cell.nodeId(2)]; - return Real3x3(Real3(vertex1.y - vertex2.y, vertex2.x - vertex1.x, 0), - Real3(vertex2.y - vertex0.y, vertex0.x - vertex2.x, 0), - Real3(vertex0.y - vertex1.y, vertex1.x - vertex0.x, 0)); + Real A2 = ((vertex1.x - vertex0.x) * (vertex2.y - vertex0.y) - (vertex2.x - vertex0.x) * (vertex1.y - vertex0.y)); + + return Real3x3(Real3((vertex1.y - vertex2.y)/A2, (vertex2.x - vertex1.x)/A2, 0), + Real3((vertex2.y - vertex0.y)/A2, (vertex0.x - vertex2.x)/A2, 0), + Real3((vertex0.y - vertex1.y)/A2, (vertex1.x - vertex0.x)/A2, 0)); } }; From 456d7de752bbdc5c8a18f17713e4fa22af2b1ba3 Mon Sep 17 00:00:00 2001 From: Mohd Afeef Badri <52162083+mohd-afeef-badri@users.noreply.github.com> Date: Mon, 5 Aug 2024 16:47:42 +0200 Subject: [PATCH 07/32] new overload operator for fixed matrix to mix it with Real3x3 --- femutils/FemUtils.h | 33 +++++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/femutils/FemUtils.h b/femutils/FemUtils.h index efe01c3..502659a 100644 --- a/femutils/FemUtils.h +++ b/femutils/FemUtils.h @@ -23,6 +23,7 @@ #include #include +#include #include @@ -165,6 +166,38 @@ inline FixedMatrix<3, 3> operator^(const Arcane::Real3& lhs, const Arcane::Real3 return result; } +/*---------------------------------------------------------------------------*/ +// Define the conversion from Real3x3 to FixedMatrix<3, 3> +/*---------------------------------------------------------------------------*/ +inline FixedMatrix<3, 3> convertReal3x3ToFixedMatrix(const Arcane::Real3x3& rhs) +{ + FixedMatrix<3, 3> result; + for (Arcane::Int32 i = 0; i < 3; ++i) { + for (Arcane::Int32 j = 0; j < 3; ++j) { + result(i, j) = rhs(i, j); + } + } + return result; +} + +/*---------------------------------------------------------------------------*/ +// Overload operator+ to handle FixedMatrix<3, 3> + Real3x3 +/*---------------------------------------------------------------------------*/ +inline FixedMatrix<3, 3> operator+(const FixedMatrix<3, 3>& lhs, const Arcane::Real3x3& rhs) +{ + FixedMatrix<3, 3> converted_rhs = convertReal3x3ToFixedMatrix(rhs); + return lhs + converted_rhs; +} + +/*---------------------------------------------------------------------------*/ +// Overload operator+ to handle Real3x3 + FixedMatrix<3, 3> +/*---------------------------------------------------------------------------*/ +inline FixedMatrix<3, 3> operator+(const Arcane::Real3x3& lhs, const FixedMatrix<3, 3>& rhs) +{ + FixedMatrix<3, 3> converted_lhs = convertReal3x3ToFixedMatrix(lhs); + return converted_lhs + rhs; +} + /*---------------------------------------------------------------------------*/ /*---------------------------------------------------------------------------*/ From d0eda88c2ad0064ece034e8eedaaeab9789b57db Mon Sep 17 00:00:00 2001 From: Mohd Afeef Badri <52162083+mohd-afeef-badri@users.noreply.github.com> Date: Mon, 5 Aug 2024 17:05:41 +0200 Subject: [PATCH 08/32] New functions to calculate gradients and UV terms --- femutils/ArcaneFemFunctions.h | 76 +++++++++++++++++++++++++++++++++-- 1 file changed, 73 insertions(+), 3 deletions(-) diff --git a/femutils/ArcaneFemFunctions.h b/femutils/ArcaneFemFunctions.h index c7cd6ed..8941b70 100644 --- a/femutils/ArcaneFemFunctions.h +++ b/femutils/ArcaneFemFunctions.h @@ -99,9 +99,79 @@ class ArcaneFemFunctions Real A2 = ((vertex1.x - vertex0.x) * (vertex2.y - vertex0.y) - (vertex2.x - vertex0.x) * (vertex1.y - vertex0.y)); - return Real3x3(Real3((vertex1.y - vertex2.y)/A2, (vertex2.x - vertex1.x)/A2, 0), - Real3((vertex2.y - vertex0.y)/A2, (vertex0.x - vertex2.x)/A2, 0), - Real3((vertex0.y - vertex1.y)/A2, (vertex1.x - vertex0.x)/A2, 0)); + return Real3x3(Real3((vertex1.y - vertex2.y) / A2, (vertex2.x - vertex1.x) / A2, 0), + Real3((vertex2.y - vertex0.y) / A2, (vertex0.x - vertex2.x) / A2, 0), + Real3((vertex0.y - vertex1.y) / A2, (vertex1.x - vertex0.x) / A2, 0)); + } + + /*---------------------------------------------------------------------------*/ + /** + * @brief Computes the X gradients of basis functions N for P1 triangles. + * + * This method calculates gradient operator ∂/∂x of Ni for a given P1 + * cell with i = 1,..,3 for the three shape function Ni hence output + * is a vector of size 3 + * + * ∂N/∂x = [ ∂N1/∂x ∂N1/∂x ∂N3/∂x ] + * + * ∂N/∂x = 1/(2A) [ y2​−y3 y3−y1 y1−y2 ] + */ + /*---------------------------------------------------------------------------*/ + + static inline Real3 computeGradientXTria3(Cell cell, VariableNodeReal3& node_coord) + { + Real3 vertex0 = node_coord[cell.nodeId(0)]; + Real3 vertex1 = node_coord[cell.nodeId(1)]; + Real3 vertex2 = node_coord[cell.nodeId(2)]; + + Real A2 = ((vertex1.x - vertex0.x) * (vertex2.y - vertex0.y) - (vertex2.x - vertex0.x) * (vertex1.y - vertex0.y)); + + return Real3((vertex1.y - vertex2.y) / A2, (vertex2.y - vertex0.y) / A2, (vertex0.y - vertex1.y) / A2); + } + + /*---------------------------------------------------------------------------*/ + /** + * @brief Computes the Y gradients of basis functions N for P1 triangles. + * + * This method calculates gradient operator ∂/∂y of Ni for a given P1 + * cell with i = 1,..,3 for the three shape function Ni hence output + * is a vector of size 3 + * + * ∂N/∂x = [ ∂N1/∂y ∂N1/∂y ∂N3/∂y ] + * + * ∂N/∂x = 1/(2A) [ x3​−x2 x1−x3 x2−x1 ] + */ + /*---------------------------------------------------------------------------*/ + + static inline Real3 computeGradientYTria3(Cell cell, VariableNodeReal3& node_coord) + { + Real3 vertex0 = node_coord[cell.nodeId(0)]; + Real3 vertex1 = node_coord[cell.nodeId(1)]; + Real3 vertex2 = node_coord[cell.nodeId(2)]; + + Real A2 = ((vertex1.x - vertex0.x) * (vertex2.y - vertex0.y) - (vertex2.x - vertex0.x) * (vertex1.y - vertex0.y)); + + return Real3((vertex2.x - vertex1.x) / A2, (vertex0.x - vertex2.x) / A2, (vertex1.x - vertex0.x) / A2); + } + + /*---------------------------------------------------------------------------*/ + /** + * @brief Computes the integral (u*v) for P1 triangles. + * + * here the element matrix will read + * + * [ 1/6 1/12 1/12 ] + * a = [ 1/12 1/6 1/12 ] + * [ 1/12 1/12 1/6 ] + * + */ + /*---------------------------------------------------------------------------*/ + + static inline Real3x3 computeUVTria3(Cell cell, VariableNodeReal3& node_coord) + { + Real aii = 1. / 6.; + Real aij = 1. / 12.; + return Real3x3(Real3(aii, aij, aij), Real3(aij, aii, aij), Real3(aij, aij, aii)); } }; From a5936f9ec71b3f8106c969f2d74209fc98200590 Mon Sep 17 00:00:00 2001 From: Mohd Afeef Badri <52162083+mohd-afeef-badri@users.noreply.github.com> Date: Mon, 5 Aug 2024 17:07:56 +0200 Subject: [PATCH 09/32] use ArcaneFemFunctions --- acoustics/FemModule.cc | 44 ++++++++++++++---------------------------- 1 file changed, 14 insertions(+), 30 deletions(-) diff --git a/acoustics/FemModule.cc b/acoustics/FemModule.cc index 7ffac5d..e2c2aa5 100644 --- a/acoustics/FemModule.cc +++ b/acoustics/FemModule.cc @@ -245,8 +245,7 @@ _assembleLinearOperator() * Steps involved: * 1. Calculate the area of the triangle. * 2. Compute the gradients of the shape functions. - * 3. Normalize the gradients by the area. - * 4. Compute the element matrix: + * 3. Compute the element matrix: * 4.1 first compute -(u.dx * v.dx + u.dy * v.dy) terms * 4.2 then (kc2 * u * v) term adjusted */ @@ -263,40 +262,25 @@ _computeElementMatrixTRIA3(Cell cell) Real3x3 gradN = ArcaneFemFunctions::computeGradientTria3(cell, m_node_coord); // step 3 - Real A2 = 2.0 * area; - gradN /= A2; - - // step 4 - FixedMatrix<3, 3> elementMatrix; - - // step 4.1 - elementMatrix(0, 0) = -area * Arcane::math::dot(gradN[0], gradN[0]); - elementMatrix(0, 1) = -area * Arcane::math::dot(gradN[0], gradN[1]); - elementMatrix(0, 2) = -area * Arcane::math::dot(gradN[0], gradN[2]); - - elementMatrix(1, 0) = -area * Arcane::math::dot(gradN[1], gradN[0]); - elementMatrix(1, 1) = -area * Arcane::math::dot(gradN[1], gradN[1]); - elementMatrix(1, 2) = -area * Arcane::math::dot(gradN[1], gradN[2]); - - elementMatrix(2, 0) = -area * Arcane::math::dot(gradN[2], gradN[0]); - elementMatrix(2, 1) = -area * Arcane::math::dot(gradN[2], gradN[1]); - elementMatrix(2, 2) = -area * Arcane::math::dot(gradN[2], gradN[2]); + Real3 dxU = {gradN[0][0], gradN[1][0], gradN[2][0]}; + Real3 dyU = {gradN[0][1], gradN[1][1], gradN[2][1]}; + FixedMatrix<3, 3> elementMatrixB; // step 4.2 Real uvTerm = m_kc2 * area / 12.0; - elementMatrix(0, 0) += uvTerm * 2.; - elementMatrix(0, 1) += uvTerm; - elementMatrix(0, 2) += uvTerm; + elementMatrixB(0, 0) = uvTerm * 2.; + elementMatrixB(0, 1) = uvTerm; + elementMatrixB(0, 2) = uvTerm; - elementMatrix(1, 0) += uvTerm; - elementMatrix(1, 1) += uvTerm * 2.; - elementMatrix(1, 2) += uvTerm; + elementMatrixB(1, 0) = uvTerm; + elementMatrixB(1, 1) = uvTerm * 2.; + elementMatrixB(1, 2) = uvTerm; - elementMatrix(2, 0) += uvTerm; - elementMatrix(2, 1) += uvTerm; - elementMatrix(2, 2) += uvTerm * 2.; + elementMatrixB(2, 0) = uvTerm; + elementMatrixB(2, 1) = uvTerm; + elementMatrixB(2, 2) = uvTerm * 2.; - return elementMatrix; + return -area* (dxU ^ dxU) - area* (dyU ^ dyU) + elementMatrixB; } /*---------------------------------------------------------------------------*/ From a21a4fe1b7d4e43d44de774dbf219aa412871cb0 Mon Sep 17 00:00:00 2001 From: Mohd Afeef Badri <52162083+mohd-afeef-badri@users.noreply.github.com> Date: Mon, 5 Aug 2024 17:14:41 +0200 Subject: [PATCH 10/32] use ArcaneFemFunctions --- acoustics/FemModule.cc | 31 +++++++------------------------ 1 file changed, 7 insertions(+), 24 deletions(-) diff --git a/acoustics/FemModule.cc b/acoustics/FemModule.cc index e2c2aa5..57d750c 100644 --- a/acoustics/FemModule.cc +++ b/acoustics/FemModule.cc @@ -244,43 +244,26 @@ _assembleLinearOperator() * * Steps involved: * 1. Calculate the area of the triangle. - * 2. Compute the gradients of the shape functions. - * 3. Compute the element matrix: - * 4.1 first compute -(u.dx * v.dx + u.dy * v.dy) terms - * 4.2 then (kc2 * u * v) term adjusted + * 2. Compute the integral U*V term. + * 3. Compute the gradients of the shape functions. + * 4. Return -(u.dx * v.dx + u.dy * v.dy) + kc2 * u * v ; */ /*---------------------------------------------------------------------------*/ FixedMatrix<3, 3> FemModule:: _computeElementMatrixTRIA3(Cell cell) { - // step 1 Real area = ArcaneFemFunctions::computeAreaTriangle3(cell, m_node_coord); // step 2 - Real3x3 gradN = ArcaneFemFunctions::computeGradientTria3(cell, m_node_coord); + Real3x3 UV = ArcaneFemFunctions::computeUVTria3(cell, m_node_coord); // step 3 - Real3 dxU = {gradN[0][0], gradN[1][0], gradN[2][0]}; - Real3 dyU = {gradN[0][1], gradN[1][1], gradN[2][1]}; - - FixedMatrix<3, 3> elementMatrixB; - // step 4.2 - Real uvTerm = m_kc2 * area / 12.0; - elementMatrixB(0, 0) = uvTerm * 2.; - elementMatrixB(0, 1) = uvTerm; - elementMatrixB(0, 2) = uvTerm; - - elementMatrixB(1, 0) = uvTerm; - elementMatrixB(1, 1) = uvTerm * 2.; - elementMatrixB(1, 2) = uvTerm; - - elementMatrixB(2, 0) = uvTerm; - elementMatrixB(2, 1) = uvTerm; - elementMatrixB(2, 2) = uvTerm * 2.; + Real3 dxU = ArcaneFemFunctions::computeGradientXTria3(cell, m_node_coord); + Real3 dyU = ArcaneFemFunctions::computeGradientYTria3(cell, m_node_coord); - return -area* (dxU ^ dxU) - area* (dyU ^ dyU) + elementMatrixB; + return -area * (dxU ^ dxU) - area * (dyU ^ dyU) + m_kc2 * area * UV; } /*---------------------------------------------------------------------------*/ From a1d5469f84303b571cfdc124638a693d53889c00 Mon Sep 17 00:00:00 2001 From: Mohd Afeef Badri <52162083+mohd-afeef-badri@users.noreply.github.com> Date: Tue, 6 Aug 2024 15:59:20 +0200 Subject: [PATCH 11/32] split file --- acoustics/CMakeLists.txt | 1 + acoustics/FemModule.cc | 86 +++++---------------------------------- acoustics/FemModule.h | 87 ++++++++++++++++++++++++++++++++++++++++ 3 files changed, 97 insertions(+), 77 deletions(-) create mode 100644 acoustics/FemModule.h diff --git a/acoustics/CMakeLists.txt b/acoustics/CMakeLists.txt index 8ce38af..517f3af 100644 --- a/acoustics/CMakeLists.txt +++ b/acoustics/CMakeLists.txt @@ -1,4 +1,5 @@ add_executable(Acoustics + FemModule.h FemModule.cc main.cc Fem_axl.h diff --git a/acoustics/FemModule.cc b/acoustics/FemModule.cc index 57d750c..c42f8cc 100644 --- a/acoustics/FemModule.cc +++ b/acoustics/FemModule.cc @@ -11,74 +11,24 @@ /*---------------------------------------------------------------------------*/ /*---------------------------------------------------------------------------*/ -#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; +#include "FemModule.h" /*---------------------------------------------------------------------------*/ /** - * @brief A module for finite element method. + * @brief Initializes the FemModule at the start of the simulation. * - * This class handles the initialization and computation for finite element - * method (FEM) simulations, providing methods to set up and solve linear - * systems, assemble FEM operators, and perform result checks. + * This method initializes degrees of freedom (DoFs) on nodes. */ /*---------------------------------------------------------------------------*/ -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); - } - - void compute() override; //! Method called at each iteration - void startInit() override; //! Method called at the beginning of the simulation - VersionInfo versionInfo() const override { return VersionInfo(1, 0, 0); } - - private: - - Real m_kc2; - - DoFLinearSystem m_linear_system; - IItemFamily* m_dof_family = nullptr; - FemDoFsOnNodes m_dofs_on_nodes; - - void _doStationarySolve(); - void _getMaterialParameters(); - void _assembleBilinearOperatorTRIA3(); - void _solve(); - void _assembleLinearOperator(); - void _validateResults(); + info() << "Module Fem INIT"; - FixedMatrix<3, 3> _computeElementMatrixTRIA3(Cell cell); -}; + m_dofs_on_nodes.initialize(mesh(), 1); + m_dof_family = m_dofs_on_nodes.dofFamily(); +} /*---------------------------------------------------------------------------*/ /** @@ -111,23 +61,6 @@ compute() _doStationarySolve(); } -/*---------------------------------------------------------------------------*/ -/** - * @brief Initializes the FemModule at the start of the simulation. - * - * This method initializes degrees of freedom (DoFs) on nodes. - */ -/*---------------------------------------------------------------------------*/ - -void FemModule:: -startInit() -{ - info() << "Module Fem INIT"; - - m_dofs_on_nodes.initialize(mesh(), 1); - m_dof_family = m_dofs_on_nodes.dofFamily(); -} - /*---------------------------------------------------------------------------*/ /** * @brief Performs a stationary solve for the FEM system. @@ -346,7 +279,6 @@ _solve() void FemModule:: _validateResults() { - // setp 1 if (allNodes().size() < 200) { ENUMERATE_ (Node, inode, allNodes()) { diff --git a/acoustics/FemModule.h b/acoustics/FemModule.h new file mode 100644 index 0000000..242d250 --- /dev/null +++ b/acoustics/FemModule.h @@ -0,0 +1,87 @@ +// -*- 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.h (C) 2022-2024 */ +/* */ +/* Simple module to test simple FEM mechanism. */ +/*---------------------------------------------------------------------------*/ +/*---------------------------------------------------------------------------*/ +#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 A module for finite element method. + * + * This class handles the initialization and computation for finite element + * method (FEM) simulations, providing methods to set up and solve linear + * systems, assemble FEM operators, and perform result checks. + */ +/*---------------------------------------------------------------------------*/ + +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); + } + + void startInit() override; //! Method called at the beginning of the simulation + void compute() override; //! Method called at each iteration + VersionInfo versionInfo() const override { return VersionInfo(1, 0, 0); } + + private: + + Real m_kc2; + + DoFLinearSystem m_linear_system; + IItemFamily* m_dof_family = nullptr; + FemDoFsOnNodes m_dofs_on_nodes; + + void _doStationarySolve(); + void _getMaterialParameters(); + void _assembleBilinearOperatorTRIA3(); + void _solve(); + void _assembleLinearOperator(); + void _validateResults(); + + FixedMatrix<3, 3> _computeElementMatrixTRIA3(Cell cell); +}; + +#endif From 8449caf3c2b92e2ce790229bf2aded508354763b Mon Sep 17 00:00:00 2001 From: Mohd Afeef Badri <52162083+mohd-afeef-badri@users.noreply.github.com> Date: Wed, 7 Aug 2024 16:32:27 +0200 Subject: [PATCH 12/32] delegate Neumann to ArcaneFemFunctions --- acoustics/FemModule.cc | 61 +------- femutils/ArcaneFemFunctions.h | 276 ++++++++++++++++++++++++---------- 2 files changed, 200 insertions(+), 137 deletions(-) diff --git a/acoustics/FemModule.cc b/acoustics/FemModule.cc index c42f8cc..2df384f 100644 --- a/acoustics/FemModule.cc +++ b/acoustics/FemModule.cc @@ -100,11 +100,7 @@ _getMaterialParameters() * * This method assembles the RHS 'b' of the linear system by: * 1. Initializing the RHS values to zero. - * 2. Iterating over Neumann boundary conditions to apply constant flux terms. - * 2.1 The integral of the flux term is computed over the boundary faces. - * 2.2 Handles both scalar and vector flux components. - * 3. For each boundary face, the edge length and normal are computed. - * 4. The contribution to the RHS is calculated and accumulated for each node. + * 2. Iterating over Neumann boundary conditions to apply it to RHS. */ /*---------------------------------------------------------------------------*/ @@ -117,54 +113,11 @@ _assembleLinearOperator() VariableDoFReal& rhs_values(m_linear_system.rhsVariable()); // Temporary variable for RHS values rhs_values.fill(0.0); - auto node_dof(m_dofs_on_nodes.nodeDoFConnectivityView()); + const auto node_dof(m_dofs_on_nodes.nodeDoFConnectivityView()); // setp 2 for (const auto& bs : options()->neumannBoundaryCondition()) { - FaceGroup group = bs->surface(); - - Real value = 0.0; - Real valueX = 0.0; - Real valueY = 0.0; - bool hasValue = bs->value.isPresent(); - bool hasValueX = bs->valueX.isPresent(); - bool hasValueY = bs->valueY.isPresent(); - - if (hasValue) { - value = bs->value(); - } - else { - if (hasValueX) - valueX = bs->valueX(); - if (hasValueY) - valueY = bs->valueY(); - } - - // setp 2.1 - ENUMERATE_ (Face, iface, group) { - Face face = *iface; - - // step 3. - Real length = ArcaneFemFunctions::computeEdgeLength2(face, m_node_coord); - Real2 normal = ArcaneFemFunctions::computeEdgeNormal2(face, m_node_coord); - - // step 4 - for (Node node : iface->nodes()) { - if (!node.isOwn()) - continue; - Real rhs_value = 0.0; - - // step 2.2 - if (hasValue) { - rhs_value = value * length / 2.0; - } - else { - rhs_value = (normal.x * valueX + normal.y * valueY) * length / 2.0; - } - - rhs_values[node_dof.dofId(node, 0)] += rhs_value; - } - } + ArcaneFemFunctions::BoundaryConditions2D::applyNeumannToRhs(bs, node_dof, m_node_coord, rhs_values); } } @@ -187,14 +140,14 @@ FixedMatrix<3, 3> FemModule:: _computeElementMatrixTRIA3(Cell cell) { // step 1 - Real area = ArcaneFemFunctions::computeAreaTriangle3(cell, m_node_coord); + Real area = ArcaneFemFunctions::MeshOperation::computeAreaTriangle3(cell, m_node_coord); // step 2 - Real3x3 UV = ArcaneFemFunctions::computeUVTria3(cell, m_node_coord); + Real3x3 UV = ArcaneFemFunctions::FeOperation2D::computeUVTria3(cell, m_node_coord); // step 3 - Real3 dxU = ArcaneFemFunctions::computeGradientXTria3(cell, m_node_coord); - Real3 dyU = ArcaneFemFunctions::computeGradientYTria3(cell, m_node_coord); + 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 * UV; } diff --git a/femutils/ArcaneFemFunctions.h b/femutils/ArcaneFemFunctions.h index 8941b70..45b5a1d 100644 --- a/femutils/ArcaneFemFunctions.h +++ b/femutils/ArcaneFemFunctions.h @@ -1,78 +1,115 @@ #ifndef ARCANE_FEM_FUNCTIONS_H #define ARCANE_FEM_FUNCTIONS_H -#include - using namespace Arcane; +/*---------------------------------------------------------------------------*/ +/** + * @brief Contains various functions & operations related to FEM calculations. + * + * The class provides methods organized into different nested classes for: + * - MeshOperation: Mesh related operations. + * - FeOperation2D: Finite element operations at element level. + * - BoundaryConditions2D: Boundary condition realated operations. + */ +/*---------------------------------------------------------------------------*/ + class ArcaneFemFunctions { public: /*---------------------------------------------------------------------------*/ /** - * @brief Computes the area of a triangle defined by three nodes. + * @brief Provides methods for various mesh-related operations. * - * This method calculates the area using the determinant formula for a triangle. - * The area is computed as half the value of the determinant of the matrix - * formed by the coordinates of the triangle's vertices. + * This class includes static methods for computing geometric properties + * of mesh elements, such as the area of triangles, the length of edges, + * and the normal vectors of edges. */ /*---------------------------------------------------------------------------*/ - - static inline Real computeAreaTriangle3(Cell cell, VariableNodeReal3& node_coord) + class MeshOperation { - Real3 vertex0 = node_coord[cell.nodeId(0)]; - Real3 vertex1 = node_coord[cell.nodeId(1)]; - Real3 vertex2 = node_coord[cell.nodeId(2)]; - return 0.5 * ((vertex1.x - vertex0.x) * (vertex2.y - vertex0.y) - (vertex2.x - vertex0.x) * (vertex1.y - vertex0.y)); - } + public: - /*---------------------------------------------------------------------------*/ - /** - * @brief Computes the length of the edge defined by a given face. - * - * This method calculates Euclidean distance between the two nodes of the face. - */ - /*---------------------------------------------------------------------------*/ + /*---------------------------------------------------------------------------*/ + /** + * @brief Computes the area of a triangle defined by three nodes. + * + * This method calculates the area using the determinant formula for a triangle. + * The area is computed as half the value of the determinant of the matrix + * formed by the coordinates of the triangle's vertices. + */ + /*---------------------------------------------------------------------------*/ + + static inline Real computeAreaTriangle3(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)]; - static inline Real computeEdgeLength2(Face face, VariableNodeReal3& node_coord) - { - Real3 vertex0 = node_coord[face.nodeId(0)]; - Real3 vertex1 = node_coord[face.nodeId(1)]; + return 0.5 * ((vertex1.x - vertex0.x) * (vertex2.y - vertex0.y) - (vertex2.x - vertex0.x) * (vertex1.y - vertex0.y)); + } - Real dx = vertex1.x - vertex0.x; - Real dy = vertex1.y - vertex0.y; + /*---------------------------------------------------------------------------*/ + /** + * @brief Computes the length of the edge defined by a given face. + * + * This method calculates Euclidean distance between the two nodes of the face. + */ + /*---------------------------------------------------------------------------*/ - return math::sqrt(dx * dx + dy * dy); - } + static inline Real computeEdgeLength2(Face face, const VariableNodeReal3& node_coord) + { + Real3 vertex0 = node_coord[face.nodeId(0)]; + Real3 vertex1 = node_coord[face.nodeId(1)]; - /*---------------------------------------------------------------------------*/ - /** + Real dx = vertex1.x - vertex0.x; + Real dy = vertex1.y - vertex0.y; + + return math::sqrt(dx * dx + dy * dy); + } + + /*---------------------------------------------------------------------------*/ + /** * @brief Computes the normalized edge normal for a given face. * * This method calculates normal vector to the edge defined by nodes of the face, * normalizes it, and ensures the correct orientation. */ - /*---------------------------------------------------------------------------*/ + /*---------------------------------------------------------------------------*/ - static inline Real2 computeEdgeNormal2(Face face, VariableNodeReal3& node_coord) - { - Real3 vertex0 = node_coord[face.nodeId(0)]; - Real3 vertex1 = node_coord[face.nodeId(1)]; + static inline Real2 computeEdgeNormal2(Face face, const VariableNodeReal3& node_coord) + { + Real3 vertex0 = node_coord[face.nodeId(0)]; + Real3 vertex1 = node_coord[face.nodeId(1)]; - if (!face.isSubDomainBoundaryOutside()) - std::swap(vertex0, vertex1); + if (!face.isSubDomainBoundaryOutside()) + std::swap(vertex0, vertex1); - Real dx = vertex1.x - vertex0.x; - Real dy = vertex1.y - vertex0.y; - Real norm_N = math::sqrt(dx * dx + dy * dy); + Real dx = vertex1.x - vertex0.x; + Real dy = vertex1.y - vertex0.y; + Real norm_N = math::sqrt(dx * dx + dy * dy); - return { dy / norm_N, -dx / norm_N }; - } + return { dy / norm_N, -dx / norm_N }; + } + }; /*---------------------------------------------------------------------------*/ /** + * @brief Provides methods for finite element operations in 2D. + * + * This class includes static methods for calculating gradients of basis + * functions and integrals for P1 triangles in 2D finite element analysis. + */ + /*---------------------------------------------------------------------------*/ + class FeOperation2D + { + + public: + + /*---------------------------------------------------------------------------*/ + /** * @brief Computes the gradients of basis functions N for P1 triangles. * * This method calculates gradient operator ∇ Ni for a given P1 cell @@ -89,23 +126,23 @@ class ArcaneFemFunctions * * @note we can adapt the same for 3D by filling the third component */ - /*---------------------------------------------------------------------------*/ + /*---------------------------------------------------------------------------*/ - static inline Real3x3 computeGradientTria3(Cell cell, VariableNodeReal3& node_coord) - { - Real3 vertex0 = node_coord[cell.nodeId(0)]; - Real3 vertex1 = node_coord[cell.nodeId(1)]; - Real3 vertex2 = node_coord[cell.nodeId(2)]; + static inline Real3x3 computeGradientTria3(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)]; - Real A2 = ((vertex1.x - vertex0.x) * (vertex2.y - vertex0.y) - (vertex2.x - vertex0.x) * (vertex1.y - vertex0.y)); + Real A2 = ((vertex1.x - vertex0.x) * (vertex2.y - vertex0.y) - (vertex2.x - vertex0.x) * (vertex1.y - vertex0.y)); - return Real3x3(Real3((vertex1.y - vertex2.y) / A2, (vertex2.x - vertex1.x) / A2, 0), - Real3((vertex2.y - vertex0.y) / A2, (vertex0.x - vertex2.x) / A2, 0), - Real3((vertex0.y - vertex1.y) / A2, (vertex1.x - vertex0.x) / A2, 0)); - } + return Real3x3(Real3((vertex1.y - vertex2.y) / A2, (vertex2.x - vertex1.x) / A2, 0), + Real3((vertex2.y - vertex0.y) / A2, (vertex0.x - vertex2.x) / A2, 0), + Real3((vertex0.y - vertex1.y) / A2, (vertex1.x - vertex0.x) / A2, 0)); + } - /*---------------------------------------------------------------------------*/ - /** + /*---------------------------------------------------------------------------*/ + /** * @brief Computes the X gradients of basis functions N for P1 triangles. * * This method calculates gradient operator ∂/∂x of Ni for a given P1 @@ -116,21 +153,21 @@ class ArcaneFemFunctions * * ∂N/∂x = 1/(2A) [ y2​−y3 y3−y1 y1−y2 ] */ - /*---------------------------------------------------------------------------*/ + /*---------------------------------------------------------------------------*/ - static inline Real3 computeGradientXTria3(Cell cell, VariableNodeReal3& node_coord) - { - Real3 vertex0 = node_coord[cell.nodeId(0)]; - Real3 vertex1 = node_coord[cell.nodeId(1)]; - Real3 vertex2 = node_coord[cell.nodeId(2)]; + static inline Real3 computeGradientXTria3(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)]; - Real A2 = ((vertex1.x - vertex0.x) * (vertex2.y - vertex0.y) - (vertex2.x - vertex0.x) * (vertex1.y - vertex0.y)); + Real A2 = ((vertex1.x - vertex0.x) * (vertex2.y - vertex0.y) - (vertex2.x - vertex0.x) * (vertex1.y - vertex0.y)); - return Real3((vertex1.y - vertex2.y) / A2, (vertex2.y - vertex0.y) / A2, (vertex0.y - vertex1.y) / A2); - } + return Real3((vertex1.y - vertex2.y) / A2, (vertex2.y - vertex0.y) / A2, (vertex0.y - vertex1.y) / A2); + } - /*---------------------------------------------------------------------------*/ - /** + /*---------------------------------------------------------------------------*/ + /** * @brief Computes the Y gradients of basis functions N for P1 triangles. * * This method calculates gradient operator ∂/∂y of Ni for a given P1 @@ -141,21 +178,21 @@ class ArcaneFemFunctions * * ∂N/∂x = 1/(2A) [ x3​−x2 x1−x3 x2−x1 ] */ - /*---------------------------------------------------------------------------*/ + /*---------------------------------------------------------------------------*/ - static inline Real3 computeGradientYTria3(Cell cell, VariableNodeReal3& node_coord) - { - Real3 vertex0 = node_coord[cell.nodeId(0)]; - Real3 vertex1 = node_coord[cell.nodeId(1)]; - Real3 vertex2 = node_coord[cell.nodeId(2)]; + static inline Real3 computeGradientYTria3(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)]; - Real A2 = ((vertex1.x - vertex0.x) * (vertex2.y - vertex0.y) - (vertex2.x - vertex0.x) * (vertex1.y - vertex0.y)); + Real A2 = ((vertex1.x - vertex0.x) * (vertex2.y - vertex0.y) - (vertex2.x - vertex0.x) * (vertex1.y - vertex0.y)); - return Real3((vertex2.x - vertex1.x) / A2, (vertex0.x - vertex2.x) / A2, (vertex1.x - vertex0.x) / A2); - } + return Real3((vertex2.x - vertex1.x) / A2, (vertex0.x - vertex2.x) / A2, (vertex1.x - vertex0.x) / A2); + } - /*---------------------------------------------------------------------------*/ - /** + /*---------------------------------------------------------------------------*/ + /** * @brief Computes the integral (u*v) for P1 triangles. * * here the element matrix will read @@ -165,14 +202,87 @@ class ArcaneFemFunctions * [ 1/12 1/12 1/6 ] * */ - /*---------------------------------------------------------------------------*/ + /*---------------------------------------------------------------------------*/ - static inline Real3x3 computeUVTria3(Cell cell, VariableNodeReal3& node_coord) + static inline Real3x3 computeUVTria3(Cell cell, const VariableNodeReal3& node_coord) + { + Real aii = 1. / 6.; + Real aij = 1. / 12.; + return Real3x3(Real3(aii, aij, aij), Real3(aij, aii, aij), Real3(aij, aij, aii)); + } + }; + + /*---------------------------------------------------------------------------*/ + /** + * @brief Provides methods for applying boundary conditions in 2D FEM problems. + * + * This class includes static methods for applying Neumann boundary conditions + * to the right-hand side (RHS) of finite element method equations in 2D. + */ + /*---------------------------------------------------------------------------*/ + class BoundaryConditions2D { - Real aii = 1. / 6.; - Real aij = 1. / 12.; - return Real3x3(Real3(aii, aij, aij), Real3(aij, aii, aij), Real3(aij, aij, aii)); - } + public: + + /*---------------------------------------------------------------------------*/ + /** + * @brief Applies Neumann conditions to the right-hand side (RHS) values. + * + * This method updates the RHS values of the finite element method equations + * based on the provided Neumann boundary condition. The boundary condition + * can specify a value or its components along the x and y directions. + * + * @param [IN] bs : The Neumann boundary condition values. + * @param [IN] node_dof : Connectivity view for degrees of freedom at nodes. + * @param [IN] node_coord : Coordinates of the nodes in the mesh. + * @param [OUT] rhs_values : The right-hand side values to be updated. + */ + /*---------------------------------------------------------------------------*/ + + static inline void applyNeumannToRhs(const CaseOptionsFem::CaseOptionNeumannBoundaryConditionValue* bs, const Arcane::IndexedNodeDoFConnectivityView& node_dof, const Arcane::VariableNodeReal3& node_coord, Arcane::VariableDoFReal& rhs_values) + { + FaceGroup group = bs->surface(); + + Real value = 0.0; + Real valueX = 0.0; + Real valueY = 0.0; + bool hasValue = bs->value.isPresent(); + bool hasValueX = bs->valueX.isPresent(); + bool hasValueY = bs->valueY.isPresent(); + + if (hasValue) { + value = bs->value(); + } + else { + if (hasValueX) + valueX = bs->valueX(); + if (hasValueY) + valueY = bs->valueY(); + } + + ENUMERATE_ (Face, iface, group) { + Face face = *iface; + + Real length = ArcaneFemFunctions::MeshOperation::computeEdgeLength2(face, node_coord); + Real2 normal = ArcaneFemFunctions::MeshOperation::computeEdgeNormal2(face, node_coord); + + for (Node node : iface->nodes()) { + if (!node.isOwn()) + continue; + Real rhs_value = 0.0; + + if (hasValue) { + rhs_value = value * length / 2.0; + } + else { + rhs_value = (normal.x * valueX + normal.y * valueY) * length / 2.0; + } + + rhs_values[node_dof.dofId(node, 0)] += rhs_value; + } + } + } + }; }; #endif // ARCANE_FEM_FUNCTIONS_H From 993806852e9134495ee632ecf680bb617b72bb61 Mon Sep 17 00:00:00 2001 From: Mohd Afeef Badri <52162083+mohd-afeef-badri@users.noreply.github.com> Date: Thu, 8 Aug 2024 15:03:57 +0200 Subject: [PATCH 13/32] new functions For Bary Center calculation For Dirichlet BC --- femutils/ArcaneFemFunctions.h | 202 ++++++++++++++++++++++------------ 1 file changed, 130 insertions(+), 72 deletions(-) diff --git a/femutils/ArcaneFemFunctions.h b/femutils/ArcaneFemFunctions.h index 45b5a1d..44f1288 100644 --- a/femutils/ArcaneFemFunctions.h +++ b/femutils/ArcaneFemFunctions.h @@ -41,7 +41,7 @@ class ArcaneFemFunctions * formed by the coordinates of the triangle's vertices. */ /*---------------------------------------------------------------------------*/ - + static inline Real computeAreaTriangle3(Cell cell, const VariableNodeReal3& node_coord) { Real3 vertex0 = node_coord[cell.nodeId(0)]; @@ -51,6 +51,31 @@ class ArcaneFemFunctions return 0.5 * ((vertex1.x - vertex0.x) * (vertex2.y - vertex0.y) - (vertex2.x - vertex0.x) * (vertex1.y - vertex0.y)); } + /*---------------------------------------------------------------------------*/ + /** + * @brief Computes the barycenter (centroid) of a triangle. + * + * This method calculates the barycenter of a triangle defined by three nodes. + * The barycenter is computed as the average of the vertices' coordinates. + * + * @param cell The triangle cell. + * @param node_coord The coordinates of the nodes. + * @return The 2D barycenter coordinates as a `Real3` object. + */ + /*---------------------------------------------------------------------------*/ + + static inline Real3 computeBaryCenterTriangle3(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)]; + + Real Center_x = (vertex0.x + vertex1.x + vertex2.x) / 3.; + Real Center_y = (vertex0.y + vertex1.y + vertex2.y) / 3.; + + return { Center_x, Center_y, 0 }; + } + /*---------------------------------------------------------------------------*/ /** * @brief Computes the length of the edge defined by a given face. @@ -72,11 +97,11 @@ class ArcaneFemFunctions /*---------------------------------------------------------------------------*/ /** - * @brief Computes the normalized edge normal for a given face. - * - * This method calculates normal vector to the edge defined by nodes of the face, - * normalizes it, and ensures the correct orientation. - */ + * @brief Computes the normalized edge normal for a given face. + * + * This method calculates normal vector to the edge defined by nodes of the face, + * normalizes it, and ensures the correct orientation. + */ /*---------------------------------------------------------------------------*/ static inline Real2 computeEdgeNormal2(Face face, const VariableNodeReal3& node_coord) @@ -97,11 +122,11 @@ class ArcaneFemFunctions /*---------------------------------------------------------------------------*/ /** - * @brief Provides methods for finite element operations in 2D. - * - * This class includes static methods for calculating gradients of basis - * functions and integrals for P1 triangles in 2D finite element analysis. - */ + * @brief Provides methods for finite element operations in 2D. + * + * This class includes static methods for calculating gradients of basis + * functions and integrals for P1 triangles in 2D finite element analysis. + */ /*---------------------------------------------------------------------------*/ class FeOperation2D { @@ -110,22 +135,22 @@ class ArcaneFemFunctions /*---------------------------------------------------------------------------*/ /** - * @brief Computes the gradients of basis functions N for P1 triangles. - * - * This method calculates gradient operator ∇ Ni for a given P1 cell - * with i = 1,..,3 for the three shape function Ni hence output is a - * matrix of size 3x3 - * - * [ ∂N1/∂x ∂N1/∂y ∂N1/∂z ] - * ∇N = [ ∂N2/∂x ∂N2/∂y ∂N2/∂z ] - * [ ∂N3/∂x ∂N3/∂y ∂N3/∂z ] - * - * [ y2​−y3 x3​−x2 0 ] - * ∇N = 1/(2A) [ y3−y1 x1−x3 0 ] - * [ y1−y2 x2−x1 0 ] - * - * @note we can adapt the same for 3D by filling the third component - */ + * @brief Computes the gradients of basis functions N for P1 triangles. + * + * This method calculates gradient operator ∇ Ni for a given P1 cell + * with i = 1,..,3 for the three shape function Ni hence output is a + * matrix of size 3x3 + * + * [ ∂N1/∂x ∂N1/∂y ∂N1/∂z ] + * ∇N = [ ∂N2/∂x ∂N2/∂y ∂N2/∂z ] + * [ ∂N3/∂x ∂N3/∂y ∂N3/∂z ] + * + * [ y2​−y3 x3​−x2 0 ] + * ∇N = 1/(2A) [ y3−y1 x1−x3 0 ] + * [ y1−y2 x2−x1 0 ] + * + * @note we can adapt the same for 3D by filling the third component + */ /*---------------------------------------------------------------------------*/ static inline Real3x3 computeGradientTria3(Cell cell, const VariableNodeReal3& node_coord) @@ -143,16 +168,16 @@ class ArcaneFemFunctions /*---------------------------------------------------------------------------*/ /** - * @brief Computes the X gradients of basis functions N for P1 triangles. - * - * This method calculates gradient operator ∂/∂x of Ni for a given P1 - * cell with i = 1,..,3 for the three shape function Ni hence output - * is a vector of size 3 - * - * ∂N/∂x = [ ∂N1/∂x ∂N1/∂x ∂N3/∂x ] - * - * ∂N/∂x = 1/(2A) [ y2​−y3 y3−y1 y1−y2 ] - */ + * @brief Computes the X gradients of basis functions N for P1 triangles. + * + * This method calculates gradient operator ∂/∂x of Ni for a given P1 + * cell with i = 1,..,3 for the three shape function Ni hence output + * is a vector of size 3 + * + * ∂N/∂x = [ ∂N1/∂x ∂N1/∂x ∂N3/∂x ] + * + * ∂N/∂x = 1/(2A) [ y2​−y3 y3−y1 y1−y2 ] + */ /*---------------------------------------------------------------------------*/ static inline Real3 computeGradientXTria3(Cell cell, const VariableNodeReal3& node_coord) @@ -168,16 +193,16 @@ class ArcaneFemFunctions /*---------------------------------------------------------------------------*/ /** - * @brief Computes the Y gradients of basis functions N for P1 triangles. - * - * This method calculates gradient operator ∂/∂y of Ni for a given P1 - * cell with i = 1,..,3 for the three shape function Ni hence output - * is a vector of size 3 - * - * ∂N/∂x = [ ∂N1/∂y ∂N1/∂y ∂N3/∂y ] - * - * ∂N/∂x = 1/(2A) [ x3​−x2 x1−x3 x2−x1 ] - */ + * @brief Computes the Y gradients of basis functions N for P1 triangles. + * + * This method calculates gradient operator ∂/∂y of Ni for a given P1 + * cell with i = 1,..,3 for the three shape function Ni hence output + * is a vector of size 3 + * + * ∂N/∂x = [ ∂N1/∂y ∂N1/∂y ∂N3/∂y ] + * + * ∂N/∂x = 1/(2A) [ x3​−x2 x1−x3 x2−x1 ] + */ /*---------------------------------------------------------------------------*/ static inline Real3 computeGradientYTria3(Cell cell, const VariableNodeReal3& node_coord) @@ -193,15 +218,15 @@ class ArcaneFemFunctions /*---------------------------------------------------------------------------*/ /** - * @brief Computes the integral (u*v) for P1 triangles. - * - * here the element matrix will read - * - * [ 1/6 1/12 1/12 ] - * a = [ 1/12 1/6 1/12 ] - * [ 1/12 1/12 1/6 ] - * - */ + * @brief Computes the integral (u*v) for P1 triangles. + * + * here the element matrix will read + * + * [ 1/6 1/12 1/12 ] + * a = [ 1/12 1/6 1/12 ] + * [ 1/12 1/12 1/6 ] + * + */ /*---------------------------------------------------------------------------*/ static inline Real3x3 computeUVTria3(Cell cell, const VariableNodeReal3& node_coord) @@ -214,11 +239,11 @@ class ArcaneFemFunctions /*---------------------------------------------------------------------------*/ /** - * @brief Provides methods for applying boundary conditions in 2D FEM problems. - * - * This class includes static methods for applying Neumann boundary conditions - * to the right-hand side (RHS) of finite element method equations in 2D. - */ + * @brief Provides methods for applying boundary conditions in 2D FEM problems. + * + * This class includes static methods for applying Neumann boundary conditions + * to the right-hand side (RHS) of finite element method equations in 2D. + */ /*---------------------------------------------------------------------------*/ class BoundaryConditions2D { @@ -226,17 +251,17 @@ class ArcaneFemFunctions /*---------------------------------------------------------------------------*/ /** - * @brief Applies Neumann conditions to the right-hand side (RHS) values. - * - * This method updates the RHS values of the finite element method equations - * based on the provided Neumann boundary condition. The boundary condition - * can specify a value or its components along the x and y directions. - * - * @param [IN] bs : The Neumann boundary condition values. - * @param [IN] node_dof : Connectivity view for degrees of freedom at nodes. - * @param [IN] node_coord : Coordinates of the nodes in the mesh. - * @param [OUT] rhs_values : The right-hand side values to be updated. - */ + * @brief Applies Neumann conditions to the right-hand side (RHS) values. + * + * This method updates the RHS values of the finite element method equations + * based on the provided Neumann boundary condition. The boundary condition + * can specify a value or its components along the x and y directions. + * + * @param [IN] bs : The Neumann boundary condition values. + * @param [IN] node_dof : Connectivity view for degrees of freedom at nodes. + * @param [IN] node_coord : Coordinates of the nodes in the mesh. + * @param [OUT] rhs_values : The right-hand side values to be updated. + */ /*---------------------------------------------------------------------------*/ static inline void applyNeumannToRhs(const CaseOptionsFem::CaseOptionNeumannBoundaryConditionValue* bs, const Arcane::IndexedNodeDoFConnectivityView& node_dof, const Arcane::VariableNodeReal3& node_coord, Arcane::VariableDoFReal& rhs_values) @@ -282,6 +307,39 @@ class ArcaneFemFunctions } } } + + /*---------------------------------------------------------------------------*/ + /** + * @brief Applies Dirichlet boundary conditions to RHS and LHS. + * + * Updates the LHS matrix and RHS vector to enforce Dirichlet conditions. + * + * - 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] 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 applyDirichletToLhsAndRhs(const CaseOptionsFem::CaseOptionDirichletBoundaryConditionValue* bs, const Arcane::IndexedNodeDoFConnectivityView& node_dof, const Arcane::VariableNodeReal3& node_coord, FemUtils::DoFLinearSystem& m_linear_system, Arcane::VariableDoFReal& rhs_values) + { + FaceGroup group = bs->surface(); + Real value = bs->value(); + Real Penalty = bs->penalty(); + + 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); + Real u_g = Penalty * value; + rhs_values[node_dof.dofId(node, 0)] = u_g; + } + } + } + } }; }; From 4ae382889301a8d9d322c435e7129bbc51c85368 Mon Sep 17 00:00:00 2001 From: Mohd Afeef Badri <52162083+mohd-afeef-badri@users.noreply.github.com> Date: Thu, 8 Aug 2024 17:11:48 +0200 Subject: [PATCH 14/32] More functions in ArcaneFemFunctions - source term functions - Dirichlet point function - Manufactured Solution function - Quad area function --- femutils/ArcaneFemFunctions.h | 150 ++++++++++++++++++++++++++++++++-- 1 file changed, 145 insertions(+), 5 deletions(-) diff --git a/femutils/ArcaneFemFunctions.h b/femutils/ArcaneFemFunctions.h index 44f1288..b3f63f6 100644 --- a/femutils/ArcaneFemFunctions.h +++ b/femutils/ArcaneFemFunctions.h @@ -51,16 +51,32 @@ class ArcaneFemFunctions return 0.5 * ((vertex1.x - vertex0.x) * (vertex2.y - vertex0.y) - (vertex2.x - vertex0.x) * (vertex1.y - vertex0.y)); } + /*---------------------------------------------------------------------------*/ + /** + * @brief Computes the area of a quadrilateral defined by four nodes. + * + * This method calculates the area of a quadrilateral by breaking it down + * into two triangles and using the determinant formula. The area is computed + * as half the value of the determinant of the matrix formed by the coordinates + * of the quadrilateral's vertices. + */ + /*---------------------------------------------------------------------------*/ + static inline Real computeAreaQuad4(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)]; + + return 0.5 * ((vertex1.x * vertex2.y + vertex2.x * vertex3.y + vertex3.x * vertex0.y + vertex0.x * vertex1.y) - (vertex2.x * vertex1.y + vertex3.x * vertex2.y + vertex0.x * vertex3.y + vertex1.x * vertex0.y)); + } + /*---------------------------------------------------------------------------*/ /** * @brief Computes the barycenter (centroid) of a triangle. * * This method calculates the barycenter of a triangle defined by three nodes. * The barycenter is computed as the average of the vertices' coordinates. - * - * @param cell The triangle cell. - * @param node_coord The coordinates of the nodes. - * @return The 2D barycenter coordinates as a `Real3` object. */ /*---------------------------------------------------------------------------*/ @@ -249,6 +265,64 @@ class ArcaneFemFunctions { public: + /*---------------------------------------------------------------------------*/ + /** + * @brief Applies a constant source term to the RHS vector. + * + * This method adds a constant source term `qdot` to the RHS vector for each + * node in the mesh. The contribution to each node is weighted by the area of + * the triangle cell and evenly distributed among the three 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 applyTria3ConstantSourceToRhs(const Real& qdot, IMesh* mesh, const Arcane::IndexedNodeDoFConnectivityView& node_dof, const Arcane::VariableNodeReal3& node_coord, Arcane::VariableDoFReal& rhs_values) + { + ENUMERATE_ (Cell, icell, mesh->allCells()) { + Cell cell = *icell; + Real area = ArcaneFemFunctions::MeshOperation::computeAreaTriangle3(cell, node_coord); + for (Node node : cell.nodes()) { + if (node.isOwn()) + rhs_values[node_dof.dofId(node, 0)] += qdot * area / 3.; + } + } + } + + /*---------------------------------------------------------------------------*/ + /** + * @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 triangle cell and evenly distributed among the three 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 applyTria3ManufacturedSourceToRhs(IBinaryMathFunctor* manufactured_source, IMesh* mesh, const Arcane::IndexedNodeDoFConnectivityView& node_dof, const Arcane::VariableNodeReal3& node_coord, Arcane::VariableDoFReal& rhs_values) + { + ENUMERATE_ (Cell, icell, mesh->allCells()) { + Cell cell = *icell; + Real area = ArcaneFemFunctions::MeshOperation::computeAreaTriangle3(cell, node_coord); + Real3 bcenter = ArcaneFemFunctions::MeshOperation::computeBaryCenterTriangle3(cell, node_coord); + + for (Node node : cell.nodes()) { + if (node.isOwn()) + rhs_values[node_dof.dofId(node, 0)] += manufactured_source->apply(area / 3., bcenter); + } + } + } + /*---------------------------------------------------------------------------*/ /** * @brief Applies Neumann conditions to the right-hand side (RHS) values. @@ -313,7 +387,7 @@ class ArcaneFemFunctions * @brief Applies Dirichlet boundary conditions to RHS and LHS. * * Updates the LHS matrix and RHS vector to enforce Dirichlet conditions. - * + * * - 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`. * @@ -340,6 +414,72 @@ class ArcaneFemFunctions } } } + + /*---------------------------------------------------------------------------*/ + /** + * @brief Applies Point 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] 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 applyPointDirichletToLhsAndRhs(const CaseOptionsFem::CaseOptionDirichletPointConditionValue* bs, const Arcane::IndexedNodeDoFConnectivityView& node_dof, const Arcane::VariableNodeReal3& node_coord, FemUtils::DoFLinearSystem& m_linear_system, Arcane::VariableDoFReal& rhs_values) + { + NodeGroup group = bs->node(); + Real value = bs->value(); + Real Penalty = bs->penalty(); + + ENUMERATE_ (Node, inode, group) { + Node node = *inode; + if (node.isOwn()) { + m_linear_system.matrixSetValue(node_dof.dofId(node, 0), node_dof.dofId(node, 0), Penalty); + Real u_g = Penalty * value; + rhs_values[node_dof.dofId(node, 0)] = u_g; + } + } + } + + /*---------------------------------------------------------------------------*/ + /** + * @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* manufactured_dirichlet, const Arcane::Real& lambda, const Arcane::FaceGroup& group, const CaseOptionsFem::CaseOptionMmanufacturedDirichletConditionValue* bs, const Arcane::IndexedNodeDoFConnectivityView& node_dof, const Arcane::VariableNodeReal3& node_coord, FemUtils::DoFLinearSystem& m_linear_system, Arcane::VariableDoFReal& rhs_values) + { + Real Penalty = bs->penalty(); + + 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; + } + } + } + } }; }; From 9059d09a892edd8b5960dc057da0676ba4c2d6b5 Mon Sep 17 00:00:00 2001 From: Mohd Afeef Badri <52162083+mohd-afeef-badri@users.noreply.github.com> Date: Fri, 9 Aug 2024 11:42:53 +0200 Subject: [PATCH 15/32] Introducing Real4 for Quads - Useful for Quad/Tetra mesh with 4 entities - Introduced outer product operator ^ for 4x4 --- femutils/FemUtils.h | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/femutils/FemUtils.h b/femutils/FemUtils.h index 502659a..337f4f3 100644 --- a/femutils/FemUtils.h +++ b/femutils/FemUtils.h @@ -30,6 +30,17 @@ /*---------------------------------------------------------------------------*/ /*---------------------------------------------------------------------------*/ +struct Real4 +{ + Arcane::Real data[4]; + + Arcane::Real& operator[](std::size_t i) { return data[i]; } + const Arcane::Real& operator[](std::size_t i) const { return data[i]; } +}; + +/*---------------------------------------------------------------------------*/ +/*---------------------------------------------------------------------------*/ + namespace Arcane::FemUtils { @@ -166,6 +177,20 @@ inline FixedMatrix<3, 3> operator^(const Arcane::Real3& lhs, const Arcane::Real3 return result; } +/*---------------------------------------------------------------------------*/ +// Outer product of two Real4 vectors to produce a FixedMatrix<4, 4> +/*---------------------------------------------------------------------------*/ +inline FixedMatrix<4, 4> operator^(const Real4& lhs, const Real4& rhs) +{ + FixedMatrix<4, 4> result; + for (Arcane::Int32 i = 0; i < 4; ++i) { + for (Arcane::Int32 j = 0; j < 4; ++j) { + result(i, j) = lhs[i] * rhs[j]; + } + } + return result; +} + /*---------------------------------------------------------------------------*/ // Define the conversion from Real3x3 to FixedMatrix<3, 3> /*---------------------------------------------------------------------------*/ From 2eac9a0a47c8f90e946819936af88fd7e5a88ca4 Mon Sep 17 00:00:00 2001 From: Mohd Afeef Badri <52162083+mohd-afeef-badri@users.noreply.github.com> Date: Fri, 9 Aug 2024 11:44:31 +0200 Subject: [PATCH 16/32] Quad mesh gradients - introduced computeGradientXQuad4 and computeGradientYQuad4 --- femutils/ArcaneFemFunctions.h | 68 +++++++++++++++++++++++++++++++++++ 1 file changed, 68 insertions(+) diff --git a/femutils/ArcaneFemFunctions.h b/femutils/ArcaneFemFunctions.h index b3f63f6..d8616b3 100644 --- a/femutils/ArcaneFemFunctions.h +++ b/femutils/ArcaneFemFunctions.h @@ -1,6 +1,8 @@ #ifndef ARCANE_FEM_FUNCTIONS_H #define ARCANE_FEM_FUNCTIONS_H +#include + using namespace Arcane; /*---------------------------------------------------------------------------*/ @@ -232,6 +234,72 @@ class ArcaneFemFunctions return Real3((vertex2.x - vertex1.x) / A2, (vertex0.x - vertex2.x) / A2, (vertex1.x - vertex0.x) / A2); } + /*---------------------------------------------------------------------------*/ + /** + * @brief Computes the X gradients of basis functions N for P1 Quad. + * + * This method calculates gradient operator ∂/∂x of Ni for a given P1 + * cell with i = 1,..,4 for the four shape function Ni hence output + * is a vector of size 4 + * + * ∂N/∂x = [ ∂N1/∂x ∂N1/∂x ∂N3/∂x ∂N4/∂x ] + * + * ∂N/∂x = 1/(2A) [ y2​−y3 y3−y0 y0−y1 y1−y2 ] + */ + /*---------------------------------------------------------------------------*/ + + static inline Real4 computeGradientXQuad4(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 A2 = ((vertex1.x * vertex2.y + vertex2.x * vertex3.y + vertex3.x * vertex0.y + vertex0.x * vertex1.y) - (vertex2.x * vertex1.y + vertex3.x * vertex2.y + vertex0.x * vertex3.y + vertex1.x * vertex0.y)); + + Real4 dx; + + dx[0] = (vertex2.y - vertex3.y) / A2; + dx[1] = (vertex3.y - vertex0.y) / A2; + dx[2] = (vertex0.y - vertex1.y) / A2; + dx[3] = (vertex1.y - vertex2.y) / A2; + + return dx; + } + + /*---------------------------------------------------------------------------* + /** + * @brief Computes the Y gradients of basis functions N for P1 Quad. + * + * This method calculates gradient operator ∂/∂y of Ni for a given P1 + * cell with i = 1,..,4 for the three shape function Ni hence output + * is a vector of size 4 + * + * ∂N/∂x = [ ∂N1/∂y ∂N1/∂y ∂N3/∂y ∂N4/∂y ] + * + * ∂N/∂x = 1/(2A) [ x3​−x2 x0−x3 x1−x0 x2−x1 ] + */ + /*---------------------------------------------------------------------------*/ + + static inline Real4 computeGradientYQuad4(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 A2 = ((vertex1.x * vertex2.y + vertex2.x * vertex3.y + vertex3.x * vertex0.y + vertex0.x * vertex1.y) - (vertex2.x * vertex1.y + vertex3.x * vertex2.y + vertex0.x * vertex3.y + vertex1.x * vertex0.y)); + + Real4 dy; + + dy[0] = (vertex3.x - vertex2.x) / A2; + dy[1] = (vertex0.x - vertex3.x) / A2; + dy[2] = (vertex1.x - vertex0.x) / A2; + dy[3] = (vertex2.x - vertex1.x) / A2; + + return dy; + } + /*---------------------------------------------------------------------------*/ /** * @brief Computes the integral (u*v) for P1 triangles. From 8a81ce927103cb2128f39bcaa5f0e2afa8c8bb16 Mon Sep 17 00:00:00 2001 From: Mohd Afeef Badri <52162083+mohd-afeef-badri@users.noreply.github.com> Date: Fri, 9 Aug 2024 16:05:23 +0200 Subject: [PATCH 17/32] Generic function names --- femutils/ArcaneFemFunctions.h | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/femutils/ArcaneFemFunctions.h b/femutils/ArcaneFemFunctions.h index d8616b3..be39c12 100644 --- a/femutils/ArcaneFemFunctions.h +++ b/femutils/ArcaneFemFunctions.h @@ -339,7 +339,7 @@ class ArcaneFemFunctions * * This method adds a constant source term `qdot` to the RHS vector for each * node in the mesh. The contribution to each node is weighted by the area of - * the triangle cell and evenly distributed among the three nodes of the cell. + * the cell and evenly distributed among the number of nodes of the cell. * * @param [IN] qdot : The constant source term. * @param [IN] mesh : The mesh containing all cells. @@ -349,14 +349,14 @@ class ArcaneFemFunctions */ /*---------------------------------------------------------------------------*/ - static inline void applyTria3ConstantSourceToRhs(const Real& qdot, IMesh* mesh, const Arcane::IndexedNodeDoFConnectivityView& node_dof, const Arcane::VariableNodeReal3& node_coord, Arcane::VariableDoFReal& rhs_values) + static inline void applyConstantSourceToRhs(const Real& qdot, IMesh* mesh, const Arcane::IndexedNodeDoFConnectivityView& node_dof, const Arcane::VariableNodeReal3& node_coord, Arcane::VariableDoFReal& rhs_values) { ENUMERATE_ (Cell, icell, mesh->allCells()) { Cell cell = *icell; Real area = ArcaneFemFunctions::MeshOperation::computeAreaTriangle3(cell, node_coord); for (Node node : cell.nodes()) { if (node.isOwn()) - rhs_values[node_dof.dofId(node, 0)] += qdot * area / 3.; + rhs_values[node_dof.dofId(node, 0)] += qdot * area / cell.nbNode(); } } } @@ -367,7 +367,7 @@ class ArcaneFemFunctions * * 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 triangle cell and evenly distributed among the three nodes of the cell. + * 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. @@ -377,7 +377,7 @@ class ArcaneFemFunctions */ /*---------------------------------------------------------------------------*/ - static inline void applyTria3ManufacturedSourceToRhs(IBinaryMathFunctor* manufactured_source, IMesh* mesh, const Arcane::IndexedNodeDoFConnectivityView& node_dof, const Arcane::VariableNodeReal3& node_coord, Arcane::VariableDoFReal& rhs_values) + static inline void applyManufacturedSourceToRhs(IBinaryMathFunctor* manufactured_source, IMesh* mesh, const Arcane::IndexedNodeDoFConnectivityView& node_dof, const Arcane::VariableNodeReal3& node_coord, Arcane::VariableDoFReal& rhs_values) { ENUMERATE_ (Cell, icell, mesh->allCells()) { Cell cell = *icell; @@ -386,7 +386,7 @@ class ArcaneFemFunctions for (Node node : cell.nodes()) { if (node.isOwn()) - rhs_values[node_dof.dofId(node, 0)] += manufactured_source->apply(area / 3., bcenter); + rhs_values[node_dof.dofId(node, 0)] += manufactured_source->apply(area / cell.nbNode(), bcenter); } } } @@ -533,7 +533,7 @@ class ArcaneFemFunctions * @param [OUT] rhs_values RHS : RHS values to update. */ /*---------------------------------------------------------------------------*/ - static inline void applyManufacturedDirichletToLhsAndRhs(IBinaryMathFunctor* manufactured_dirichlet, const Arcane::Real& lambda, const Arcane::FaceGroup& group, const CaseOptionsFem::CaseOptionMmanufacturedDirichletConditionValue* bs, const Arcane::IndexedNodeDoFConnectivityView& node_dof, const Arcane::VariableNodeReal3& node_coord, FemUtils::DoFLinearSystem& m_linear_system, Arcane::VariableDoFReal& rhs_values) + static inline void applyManufacturedDirichletToLhsAndRhs(IBinaryMathFunctor* manufactured_dirichlet, const Arcane::Real& lambda, const Arcane::FaceGroup& group, const CaseOptionsFem::CaseOptionManufacturedSolutionValue* bs, const Arcane::IndexedNodeDoFConnectivityView& node_dof, const Arcane::VariableNodeReal3& node_coord, FemUtils::DoFLinearSystem& m_linear_system, Arcane::VariableDoFReal& rhs_values) { Real Penalty = bs->penalty(); From d3ddf2c2466343518758bbfc306b630c79aef250 Mon Sep 17 00:00:00 2001 From: Mohd Afeef Badri <52162083+mohd-afeef-badri@users.noreply.github.com> Date: Fri, 9 Aug 2024 16:07:36 +0200 Subject: [PATCH 18/32] ArcaneFemFunctions in Fourier --- acoustics/Fem.axl | 94 ++++ fourier/CMakeLists.txt | 1 + fourier/Fem.axl | 61 +- fourier/FemModule.cc | 769 ++++++-------------------- fourier/FemModule.h | 97 ++++ fourier/Test.manufacture.solution.arc | 7 +- 6 files changed, 424 insertions(+), 605 deletions(-) create mode 100644 fourier/FemModule.h diff --git a/acoustics/Fem.axl b/acoustics/Fem.axl index f65a9f1..34a589b 100644 --- a/acoustics/Fem.axl +++ b/acoustics/Fem.axl @@ -38,6 +38,100 @@ Type of mesh provided to the solver + + + + + Dirichlet boundary condition + + + + FaceGroup on which to apply these boundary condition + + + + + Value of the boundary 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 + + + + + + + + Dirichlet point condition + + + + NodeGroup on which to apply these point Dirichlet condition + + + + + Value of the point Dirichlet condition + + + + + Method via which Dirichlet boundary condition is imposed + + + + + Penalty value for enforcing Dirichlet condition + + + + Neumann boundary condition diff --git a/fourier/CMakeLists.txt b/fourier/CMakeLists.txt index a943ff5..22f6f42 100644 --- a/fourier/CMakeLists.txt +++ b/fourier/CMakeLists.txt @@ -1,4 +1,5 @@ add_executable(Fourier + FemModule.h FemModule.cc main.cc Fem_axl.h diff --git a/fourier/Fem.axl b/fourier/Fem.axl index 062c89f..afcdd39 100644 --- a/fourier/Fem.axl +++ b/fourier/Fem.axl @@ -45,16 +45,6 @@ Penalty value for enforcing Dirichlet condition - - - Function for Dirichlet boundary condition - - - - - Function for Source 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 + + @@ -128,6 +159,16 @@ Value of the point Dirichlet condition + + + Method via which Dirichlet boundary condition is imposed + + + + + Penalty value for enforcing Dirichlet condition + + diff --git a/fourier/FemModule.cc b/fourier/FemModule.cc index cf4840b..fc3bd3e 100644 --- a/fourier/FemModule.cc +++ b/fourier/FemModule.cc @@ -11,94 +11,62 @@ /*---------------------------------------------------------------------------*/ /*---------------------------------------------------------------------------*/ -#include -#include -#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 "FemModule.h" /*---------------------------------------------------------------------------*/ +/** + * @brief Initializes the FemModule at the start of the simulation. + * + * This method: + * - 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: + info() << "Module Fem INIT"; - //! Method called at each iteration - void compute() override; + m_dofs_on_nodes.initialize(mesh(), 1); + m_dof_family = m_dofs_on_nodes.dofFamily(); - //! Method called at the beginning of the simulation - void startInit() override; + if (options()->manufacturedSolution.isPresent()) { + const auto& bs = options()->manufacturedSolution()[0]; + + if (bs->manufacturedSource.isPresent()) { + ICaseFunction* opt_function_source = bs->manufacturedSource.function(); + IStandardFunction* scf_source = bs->manufacturedSource.standardFunction(); + if (!scf_source) + ARCANE_FATAL("No standard case function for option 'manufactured-source-condition'"); + auto* functorS = scf_source->getFunctorRealReal3ToReal(); + if (!functorS) + ARCANE_FATAL("Standard function '{0}' is not convertible to f(Real,Real3) -> Real", opt_function_source->name()); + m_manufactured_source = functorS; + } - VersionInfo versionInfo() const override - { - return VersionInfo(1, 0, 0); + if (bs->manufacturedDirichlet.isPresent()) { + ICaseFunction* opt_function = bs->manufacturedDirichlet.function(); + IStandardFunction* scf = bs->manufacturedDirichlet.standardFunction(); + if (!scf) + ARCANE_FATAL("No standard case function for option 'manufactured-dirichlet-condition'"); + auto* functor = scf->getFunctorRealReal3ToReal(); + if (!functor) + ARCANE_FATAL("Standard function '{0}' is not convertible to f(Real,Real3) -> Real", opt_function->name()); + m_manufactured_dirichlet = functor; + } } - - private: - - Real lambda; - Real qdot; - Real ElementNodes; - - DoFLinearSystem m_linear_system; - IItemFamily* m_dof_family = nullptr; - FemDoFsOnNodes m_dofs_on_nodes; - - private: - - void _doStationarySolve(); - void _getMaterialParameters(); - void _assembleBilinearOperatorTRIA3(); - void _assembleBilinearOperatorQUAD4(); - void _solve(); - void _assembleLinearOperator(); - void _applyDirichletBoundaryConditions(); - void _checkResultFile(); - FixedMatrix<3, 3> _computeElementMatrixTRIA3(Cell cell); - FixedMatrix<4, 4> _computeElementMatrixQUAD4(Cell cell); - Real _computeAreaTriangle3(Cell cell); - Real _computeAreaQuad4(Cell cell); - Real _computeEdgeLength2(Face face); - Real2 _computeEdgeNormal2(Face face); - IBinaryMathFunctor* m_manufactured_dirichlet = nullptr; - IBinaryMathFunctor* m_manufactured_source = nullptr; - -}; +} /*---------------------------------------------------------------------------*/ +/** + * @brief Performs the main computation for the FemModule. + * + * This method: + * 1. Stops the time loop after 1 iteration since the equation is steady state. + * 2. Resets, configures, and initializes the linear system. + * 3. Executes the stationary solve. + */ /*---------------------------------------------------------------------------*/ void FemModule:: @@ -126,71 +94,42 @@ compute() } /*---------------------------------------------------------------------------*/ -/*---------------------------------------------------------------------------*/ - -void FemModule:: -startInit() -{ - info() << "Module Fem INIT"; - - m_dofs_on_nodes.initialize(mesh(), 1); - m_dof_family = m_dofs_on_nodes.dofFamily(); - - // Check we have user function for manufactured boundary/source condition - if(options()->manufacturedDirichletCondition()){ - ICaseFunction* opt_function = options()->manufacturedDirichletCondition.function(); - IStandardFunction* scf = options()->manufacturedDirichletCondition.standardFunction(); - if (!scf) - ARCANE_FATAL("No standard case function for option 'manufactured-dirichlet-condition'"); - auto* functor = scf->getFunctorRealReal3ToReal(); - if (!functor) - ARCANE_FATAL("Standard function '{0}' is not convertible to f(Real,Real3) -> Real", opt_function->name()); - m_manufactured_dirichlet = functor; - } - - if(options()->manufacturedSourceCondition()){ - ICaseFunction* opt_function_source = options()->manufacturedSourceCondition.function(); - IStandardFunction* scf_source = options()->manufacturedSourceCondition.standardFunction(); - if (!scf_source) - ARCANE_FATAL("No standard case function for option 'manufactured-source-condition'"); - auto* functorS = scf_source->getFunctorRealReal3ToReal(); - if (!functorS) - ARCANE_FATAL("Standard function '{0}' is not convertible to f(Real,Real3) -> Real", opt_function_source->name()); - m_manufactured_source = functorS; - } - - -} - -/*---------------------------------------------------------------------------*/ +/** + * @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. _assembleBilinearOperatorTRIA3() 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(); - - // apply Dirichlet BC - _applyDirichletBoundaryConditions(); - - // Assemble the FEM bilinear operator (LHS - matrix A) if (options()->meshType == "QUAD4") _assembleBilinearOperatorQUAD4(); else _assembleBilinearOperatorTRIA3(); - - // Assemble the FEM linear operator (RHS - vector b) _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: + * # thermal conductivity coefficient (`lambda`) + * # heat source term (`qdot`) + * - mesh properties: + * # number of cell nodes per element based on the mesh (quad or tria) + */ /*---------------------------------------------------------------------------*/ void FemModule:: @@ -198,16 +137,12 @@ _getMaterialParameters() { info() << "Get material parameters..."; lambda = options()->lambda(); - qdot = options()->qdot(); - ElementNodes = 3.; - - if (options()->meshType == "QUAD4") - ElementNodes = 4.; + qdot = options()->qdot(); ENUMERATE_ (Cell, icell, allCells()) { Cell cell = *icell; m_cell_lambda[cell] = lambda; - } + } for (const auto& bs : options()->materialProperty()) { CellGroup group = bs->volume(); @@ -217,65 +152,23 @@ _getMaterialParameters() ENUMERATE_ (Cell, icell, group) { Cell cell = *icell; m_cell_lambda[cell] = value; - } - } -} - - -/*---------------------------------------------------------------------------*/ -/*---------------------------------------------------------------------------*/ - -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_u[node] = value; - m_u_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_u[node] = value; - m_u_dirichlet[node] = true; - } - } - - if(options()->manufacturedDirichletCondition()){ - info() << "Apply manufactured Dirichlet boundary condition to all surface"; - ENUMERATE_ (Face, iface, outerFaces()) { - for (Node node : iface->nodes()) { - m_u[node] = m_manufactured_dirichlet->apply(lambda, m_node_coord[node]); - m_u_dirichlet[node] = true; - } } } - } /*---------------------------------------------------------------------------*/ -// Assemble the FEM linear operator -// - This function enforces a Dirichlet boundary condition in a weak sense -// via the penalty method -// - The method also adds source term -// - TODO: external fluxes +/** + * @brief FEM linear operator for the current simulation step. + * + * This method constructs the linear system by assembling the LHS matrix + * and RHS vector, applying various boundary conditions and source terms. + * + * - The RHS vector is initialized to zero before applying any conditions. + * - Manufactured or constant source conditions are applied to all cells. + * - Neumann BC are applied to the RHS. + * - Dirichlet BC are applied to both the LHS and RHS using penalty method. + * - If a manufactured Dirichlet condition is specified, it is applied. + */ /*---------------------------------------------------------------------------*/ void FemModule:: @@ -283,329 +176,62 @@ _assembleLinearOperator() { info() << "Assembly of FEM linear operator "; - // 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") { + if (options()->qdot.isPresent()) + ArcaneFemFunctions::BoundaryConditions2D::applyConstantSourceToRhs(qdot, mesh(), node_dof, m_node_coord, rhs_values); - //---------------------------------------------- - // 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 - //---------------------------------------------- + for (const auto& bs : options()->neumannBoundaryCondition()) + ArcaneFemFunctions::BoundaryConditions2D::applyNeumannToRhs(bs, node_dof, m_node_coord, rhs_values); - info() << "Applying Dirichlet boundary condition via " - << options()->enforceDirichletMethod() << " method "; + for (const auto& bs : options()->dirichletBoundaryCondition()) + ArcaneFemFunctions::BoundaryConditions2D::applyDirichletToLhsAndRhs(bs, node_dof, m_node_coord, m_linear_system, rhs_values); - Real Penalty = options()->penalty(); // 1.0e30 is the default + if (options()->manufacturedSolution.isPresent()) { + const auto& bs = options()->manufacturedSolution()[0]; - ENUMERATE_ (Node, inode, ownNodes()) { - NodeLocalId node_id = *inode; - if (m_u_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_u[node_id]; - rhs_values[dof_id] = u_g; - } + if (bs->manufacturedSource.isPresent()) { + 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); } - }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_u_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_u[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 "; - - // TODO - }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 "; - - // TODO - }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}(qdot*v^h)$ - // only for noded that are non-Dirichlet - //---------------------------------------------- - if(options()->manufacturedSourceCondition()){ - ARCANE_CHECK_POINTER(m_manufactured_source); - info() << "Apply manufactured Source condition to all cells"; - ENUMERATE_ (Cell, icell, allCells()) { - Cell cell = *icell; - - Real area = _computeAreaTriangle3(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)]; - - Real Center_x = (m0.x + m1.x + m2.x)/3; - Real Center_y = (m0.y + m1.y + m2.y)/3; - Real3 bcenter = {Center_x , Center_y, 0}; - - for (Node node : cell.nodes()) { - if (!(m_u_dirichlet[node]) && node.isOwn()){ - rhs_values[node_dof.dofId(node, 0)] += m_manufactured_source->apply(area / ElementNodes, bcenter); - } - } - } - }else{ - ENUMERATE_ (Cell, icell, allCells()) { - Cell cell = *icell; - Real area = _computeAreaTriangle3(cell); - for (Node node : cell.nodes()) { - if (!(m_u_dirichlet[node]) && node.isOwn()) - rhs_values[node_dof.dofId(node, 0)] += qdot * area / ElementNodes; - } - } - } - //---------------------------------------------- - // 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(); - - info() << "Apply Neumann boundary condition to all edges on surface" << group.name(); - 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_u_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_u_dirichlet[node]) && node.isOwn()) - rhs_values[node_dof.dofId(node, 0)] += (Normal.x*valueX + Normal.y*valueY) * length / 2.; - } - } - continue; - } - - 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_u_dirichlet[node]) && node.isOwn()) - rhs_values[node_dof.dofId(node, 0)] += (Normal.x*valueX) * length / 2.; - } - } - continue; + if (bs->manufacturedDirichlet.isPresent()) { + 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); } - - 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_u_dirichlet[node]) && node.isOwn()) - rhs_values[node_dof.dofId(node, 0)] += (Normal.y*valueY) * length / 2.; - } - } - continue; - } - } } /*---------------------------------------------------------------------------*/ -/*---------------------------------------------------------------------------*/ - -Real FemModule:: -_computeAreaQuad4(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)]; - Real3 m3 = m_node_coord[cell.nodeId(3)]; - return 0.5 * ( (m1.x*m2.y + m2.x*m3.y + m3.x*m0.y + m0.x*m1.y) - -(m2.x*m1.y + m3.x*m2.y + m0.x*m3.y + m1.x*m0.y) ); -} - -/*---------------------------------------------------------------------------*/ -/*---------------------------------------------------------------------------*/ - -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; -} - -/*---------------------------------------------------------------------------*/ +/** + * @brief Computes the element matrix for a triangular element (P1 FE). + * + * This function calculates the integral of the expression: + * lambda*(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 lambda*(u.dx * v.dx + u.dy * v.dy); + */ /*---------------------------------------------------------------------------*/ FixedMatrix<3, 3> FemModule:: _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); - - FixedMatrix<2, 3> b_matrix; - b_matrix(0, 0) = dPhi0.x; - b_matrix(0, 1) = dPhi1.x; - b_matrix(0, 2) = dPhi2.x; - - 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)); - - FixedMatrix<3, 3> int_cdPi_dPj = matrixMultiplication(matrixTranspose(b_matrix), b_matrix); - int_cdPi_dPj.multInPlace(area * lambda); - - //info() << "Cell=" << cell.localId(); - //std::cout << " int_cdPi_dPj="; - //int_cdPi_dPj.dump(std::cout); - //std::cout << "\n"; - - return int_cdPi_dPj; + Real area = ArcaneFemFunctions::MeshOperation::computeAreaTriangle3(cell, m_node_coord); + + Real3 dxU = ArcaneFemFunctions::FeOperation2D::computeGradientXTria3(cell, m_node_coord); + Real3 dyU = ArcaneFemFunctions::FeOperation2D::computeGradientYTria3(cell, m_node_coord); + + return area * lambda * (dxU ^ dxU) + area * lambda * (dyU ^ dyU); } /*---------------------------------------------------------------------------*/ @@ -614,51 +240,22 @@ _computeElementMatrixTRIA3(Cell cell) FixedMatrix<4, 4> FemModule:: _computeElementMatrixQUAD4(Cell cell) { - // Get coordiantes of the quadrangular element QUAD4 - //------------------------------------------------ - // 1 o . . . . o 0 - // . . - // . . - // . . - // 2 o . . . . o 3 - //------------------------------------------------ - Real3 m0 = m_node_coord[cell.nodeId(0)]; - Real3 m1 = m_node_coord[cell.nodeId(1)]; - Real3 m2 = m_node_coord[cell.nodeId(2)]; - Real3 m3 = m_node_coord[cell.nodeId(3)]; - - Real area = _computeAreaQuad4(cell); // calculate area - - Real2 dPhi0(m2.y - m3.y, m3.x - m2.x); - Real2 dPhi1(m3.y - m0.y, m0.x - m3.x); - Real2 dPhi2(m0.y - m1.y, m1.x - m0.x); - Real2 dPhi3(m1.y - m2.y, m2.x - m1.x); - - FixedMatrix<2, 4> b_matrix; - b_matrix(0, 0) = dPhi0.x; - b_matrix(0, 1) = dPhi1.x; - b_matrix(0, 2) = dPhi2.x; - b_matrix(0, 3) = dPhi3.x; - - b_matrix(1, 0) = dPhi0.y; - b_matrix(1, 1) = dPhi1.y; - b_matrix(1, 2) = dPhi2.y; - b_matrix(1, 3) = dPhi3.y; - - b_matrix.multInPlace(1.0 / (2.0 * area)); - - FixedMatrix<4, 4> int_cdPi_dPj = matrixMultiplication(matrixTranspose(b_matrix), b_matrix); - int_cdPi_dPj.multInPlace(area * lambda); - - //info() << "Cell=" << cell.localId(); - //std::cout << " int_cdPi_dPj="; - //int_cdPi_dPj.dump(std::cout); - //std::cout << "\n"; - - return int_cdPi_dPj; + Real area = ArcaneFemFunctions::MeshOperation::computeAreaQuad4(cell, m_node_coord); + + Real4 dxU = ArcaneFemFunctions::FeOperation2D::computeGradientXQuad4(cell, m_node_coord); + Real4 dyU = ArcaneFemFunctions::FeOperation2D::computeGradientYQuad4(cell, m_node_coord); + + return area * lambda * (dxU ^ dxU) + area * lambda * (dyU ^ dyU); } /*---------------------------------------------------------------------------*/ +/** + * @brief Assembles the bilinear operator matrix for Quad elements. + * + * This method computes and assembles the global matrix A for the FEM linear + * system. For each cell (Quad), it calculates the local element matrix & + * adds the contributions to the global matrix based on the nodes of the cell. + */ /*---------------------------------------------------------------------------*/ void FemModule:: @@ -668,26 +265,14 @@ _assembleBilinearOperatorQUAD4() ENUMERATE_ (Cell, icell, allCells()) { Cell cell = *icell; - if (cell.type() != IT_Quad4) - ARCANE_FATAL("Only Quad4 cell type is supported"); - - lambda = m_cell_lambda[cell]; // lambda is always considered cell constant - auto K_e = _computeElementMatrixQUAD4(cell); // element stiffness matrix - // # assemble elementary matrix into the global one - // # elementary terms are positioned 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] + + lambda = m_cell_lambda[cell]; // lambda is always considered cell constant + auto K_e = _computeElementMatrixQUAD4(cell); // element stiffness matrix 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); } @@ -699,6 +284,13 @@ _assembleBilinearOperatorQUAD4() } /*---------------------------------------------------------------------------*/ +/** + * @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. + */ /*---------------------------------------------------------------------------*/ void FemModule:: @@ -708,26 +300,14 @@ _assembleBilinearOperatorTRIA3() ENUMERATE_ (Cell, icell, allCells()) { Cell cell = *icell; - if (cell.type() != IT_Triangle3) - ARCANE_FATAL("Only Triangle3 cell type is supported"); - - lambda = m_cell_lambda[cell]; // lambda is always considered cell constant - auto K_e = _computeElementMatrixTRIA3(cell); // element stiffness matrix - // # assemble elementary matrix into the global one - // # elementary terms are positioned 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] + + lambda = m_cell_lambda[cell]; // lambda is always considered cell constant + auto K_e = _computeElementMatrixTRIA3(cell); // element matrix 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); } @@ -739,6 +319,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:: @@ -746,20 +334,17 @@ _solve() { m_linear_system.solve(); - // Re-Apply boundary conditions because the solver has modified the value - // of u on all nodes - _applyDirichletBoundaryConditions(); - - { + { // copies solution (and optionally exact solution) to FEM output VariableDoFReal& dof_u(m_linear_system.solutionVariable()); - // Copy RHS DoF to Node u auto node_dof(m_dofs_on_nodes.nodeDoFConnectivityView()); + ENUMERATE_ (Node, inode, ownNodes()) { Node node = *inode; Real v = dof_u[node_dof.dofId(node, 0)]; m_u[node] = v; } - if(options()->manufacturedDirichletCondition()){ + + if (options()->manufacturedSolution.isPresent()) { ENUMERATE_ (Node, inode, ownNodes()) { Node node = *inode; m_u_exact[node] = m_manufactured_dirichlet->apply(lambda, m_node_coord[node]); @@ -768,40 +353,38 @@ _solve() } m_u.synchronize(); - if(options()->manufacturedDirichletCondition()) - m_u_exact.synchronize(); - // def update_T(self,T): - // """Update u value on nodes after the FE resolution""" - // for i in range(0,len(self.mesh.nodes)): - // node=self.mesh.nodes[i] - // # don't update T imposed by Dirichlet BC - // if not node.is_T_fixed: - // self.mesh.nodes[i].T=T[i] - - const bool do_print = (allNodes().size() < 200); - if (do_print) { - ENUMERATE_ (Node, inode, allNodes()) { - Node node = *inode; - info() << "T[" << node.localId() << "][" << node.uniqueId() << "] = " - << m_u[node]; - //info() << "T[]" << node.uniqueId() << " " - // << m_u[node]; - } - } + if (options()->manufacturedSolution.isPresent()) + m_u_exact.synchronize(); } /*---------------------------------------------------------------------------*/ +/** + * @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() << "u[" << node.localId() << "][" << node.uniqueId() << "] = " << m_u[node]; + } + } + String filename = options()->resultFile(); - info() << "CheckResultFile filename=" << filename; - if (filename.empty()) - return; - const double epsilon = 1.0e-4; - checkNodeResultFile(traceMng(), filename, m_u, epsilon); + info() << "ValidateResultFile filename=" << filename; + + if (!filename.empty()) + checkNodeResultFile(traceMng(), filename, m_u, 1.0e-4); } /*---------------------------------------------------------------------------*/ diff --git a/fourier/FemModule.h b/fourier/FemModule.h new file mode 100644 index 0000000..0f4a499 --- /dev/null +++ b/fourier/FemModule.h @@ -0,0 +1,97 @@ +// -*- 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 */ +/* */ +/* Simple module to test simple FEM mechanism. */ +/*---------------------------------------------------------------------------*/ +/*---------------------------------------------------------------------------*/ +#ifndef FEMMODULES_H +#define FEMMODULES_H + +/*---------------------------------------------------------------------------*/ +/*---------------------------------------------------------------------------*/ + +#include +#include +#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 A module for finite element method. + * + * This class handles the initialization and computation for finite element + * method (FEM) simulations, providing methods to set up and solve linear + * systems, assemble FEM operators, and perform result checks. + */ +/*---------------------------------------------------------------------------*/ + +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); + } + + void startInit() override; //! Method called at the beginning of the simulation + void compute() override; //! Method called at each iteration + VersionInfo versionInfo() const override { return VersionInfo(1, 0, 0); } + + private: + + Real lambda; + Real qdot; + Real ElementNodes; + + DoFLinearSystem m_linear_system; + IItemFamily* m_dof_family = nullptr; + FemDoFsOnNodes m_dofs_on_nodes; + + void _doStationarySolve(); + void _getMaterialParameters(); + void _assembleBilinearOperatorTRIA3(); + void _assembleBilinearOperatorQUAD4(); + void _solve(); + void _assembleLinearOperator(); + void _validateResults(); + + FixedMatrix<3, 3> _computeElementMatrixTRIA3(Cell cell); + FixedMatrix<4, 4> _computeElementMatrixQUAD4(Cell cell); + + IBinaryMathFunctor* m_manufactured_dirichlet = nullptr; + IBinaryMathFunctor* m_manufactured_source = nullptr; +}; + +#endif \ No newline at end of file diff --git a/fourier/Test.manufacture.solution.arc b/fourier/Test.manufacture.solution.arc index 6bff475..66bd5e1 100644 --- a/fourier/Test.manufacture.solution.arc +++ b/fourier/Test.manufacture.solution.arc @@ -28,8 +28,11 @@ - true - true + + true + true + Penalty + 1.0 From 5e512d0ea91660e81ba188f99de8a1c488ab0e21 Mon Sep 17 00:00:00 2001 From: Mohd Afeef Badri <52162083+mohd-afeef-badri@users.noreply.github.com> Date: Fri, 9 Aug 2024 17:01:31 +0200 Subject: [PATCH 19/32] delete u_dirichlet Dirichlet is now added in an optimized way --- fourier/Fem.axl | 3 --- 1 file changed, 3 deletions(-) diff --git a/fourier/Fem.axl b/fourier/Fem.axl index afcdd39..f586d25 100644 --- a/fourier/Fem.axl +++ b/fourier/Fem.axl @@ -15,9 +15,6 @@ manufactured FEM solution on nodes - - Boolean which is true if Dirichlet node - Node Coordinates from Arcane variable From c660235ddb565cdcc27b17a21b7c6bf3f95e78f0 Mon Sep 17 00:00:00 2001 From: Mohd Afeef Badri <52162083+mohd-afeef-badri@users.noreply.github.com> Date: Fri, 9 Aug 2024 17:03:28 +0200 Subject: [PATCH 20/32] Single function for Bilinear assembly Use templated Lambda fuctions for bilinear assembly --- fourier/FemModule.cc | 71 +++++++++++++++++++------------------------- fourier/FemModule.h | 6 ++-- 2 files changed, 35 insertions(+), 42 deletions(-) diff --git a/fourier/FemModule.cc b/fourier/FemModule.cc index fc3bd3e..56b8778 100644 --- a/fourier/FemModule.cc +++ b/fourier/FemModule.cc @@ -99,7 +99,7 @@ compute() * * This method follows a sequence of steps to solve FEM system: * 1. _getMaterialParameters() Retrieves material parameters via - * 2. _assembleBilinearOperatorTRIA3() Assembles the FEM matrix A + * 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 @@ -110,10 +110,7 @@ void FemModule:: _doStationarySolve() { _getMaterialParameters(); - if (options()->meshType == "QUAD4") - _assembleBilinearOperatorQUAD4(); - else - _assembleBilinearOperatorTRIA3(); + _assembleBilinearOperator(); _assembleLinearOperator(); _solve(); _validateResults(); @@ -164,10 +161,11 @@ _getMaterialParameters() * and RHS vector, applying various boundary conditions and source terms. * * - The RHS vector is initialized to zero before applying any conditions. - * - Manufactured or constant source conditions are applied to all cells. - * - Neumann BC are applied to the RHS. - * - Dirichlet BC are applied to both the LHS and RHS using penalty method. - * - If a manufactured Dirichlet condition is specified, it is applied. + * - If a constant source term is specified (`qdot`), apply it to the RHS. + * - If Neumann BC are specified applied to the RHS. + * - If Dirichlet BC are specified apply to the LHS & RHS. + * - If manufactured source conditions is specified, apply to RHS. + * - If manufactured Dirichlet BC are specified apply to the LHS & RHS. */ /*---------------------------------------------------------------------------*/ @@ -250,51 +248,43 @@ _computeElementMatrixQUAD4(Cell cell) /*---------------------------------------------------------------------------*/ /** - * @brief Assembles the bilinear operator matrix for Quad elements. + * @brief Assembles the bilinear operator matrix for the FEM linear system. * - * This method computes and assembles the global matrix A for the FEM linear - * system. For each cell (Quad), it calculates the local element matrix & - * adds the contributions to the global matrix based on the nodes of the cell. + * The method calls the right function for RHS assembly given as mesh type. */ /*---------------------------------------------------------------------------*/ void FemModule:: -_assembleBilinearOperatorQUAD4() +_assembleBilinearOperator() { - auto node_dof(m_dofs_on_nodes.nodeDoFConnectivityView()); - - ENUMERATE_ (Cell, icell, allCells()) { - Cell cell = *icell; - - lambda = m_cell_lambda[cell]; // lambda is always considered cell constant - auto K_e = _computeElementMatrixQUAD4(cell); // element stiffness matrix - 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()) { - m_linear_system.matrixAddValue(node_dof.dofId(node1, 0), node_dof.dofId(node2, 0), v); - } - ++n2_index; - } - ++n1_index; - } - } + if (options()->meshType == "QUAD4") + _assembleBilinear<4>([this](const Cell& cell) { + return _computeElementMatrixQUAD4(cell); + }); + 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 triangular elements. + * @brief Assembles the bilinear operator matrix for the FEM linear system. * - * 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. + * The method performs the following steps: + * - Iterates over all cells in the mesh. + * - For each cell, retrieves the cell-specific constant `lambda`. + * - Computes the element matrix using the provided `compute_element_matrix` function. + * - Assembles the global matrix by adding the 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()); @@ -302,7 +292,7 @@ _assembleBilinearOperatorTRIA3() Cell cell = *icell; lambda = m_cell_lambda[cell]; // lambda is always considered cell constant - 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; @@ -318,6 +308,7 @@ _assembleBilinearOperatorTRIA3() } } + /*---------------------------------------------------------------------------*/ /** * @brief Solves the linear system and updates the solution vector. diff --git a/fourier/FemModule.h b/fourier/FemModule.h index 0f4a499..0df60c6 100644 --- a/fourier/FemModule.h +++ b/fourier/FemModule.h @@ -81,8 +81,7 @@ class FemModule void _doStationarySolve(); void _getMaterialParameters(); - void _assembleBilinearOperatorTRIA3(); - void _assembleBilinearOperatorQUAD4(); + void _assembleBilinearOperator(); void _solve(); void _assembleLinearOperator(); void _validateResults(); @@ -92,6 +91,9 @@ class FemModule IBinaryMathFunctor* m_manufactured_dirichlet = nullptr; IBinaryMathFunctor* m_manufactured_source = nullptr; + + template + void _assembleBilinear( const std::function(const Cell&)>& compute_element_matrix); }; #endif \ No newline at end of file From 5a6628217e127fe343e55c43845c0e4adf3852a4 Mon Sep 17 00:00:00 2001 From: Mohd Afeef Badri <52162083+mohd-afeef-badri@users.noreply.github.com> Date: Fri, 9 Aug 2024 17:26:27 +0200 Subject: [PATCH 21/32] minor cleanup --- fourier/FemModule.cc | 48 +++++++++++++++++++------------------------- 1 file changed, 21 insertions(+), 27 deletions(-) diff --git a/fourier/FemModule.cc b/fourier/FemModule.cc index 56b8778..67b3850 100644 --- a/fourier/FemModule.cc +++ b/fourier/FemModule.cc @@ -17,7 +17,6 @@ /** * @brief Initializes the FemModule at the start of the simulation. * - * This method: * - initializes degrees of freedom (DoFs) on nodes. * - builds support for manufactured test case (optional). */ @@ -62,10 +61,9 @@ startInit() /** * @brief Performs the main computation for the FemModule. * - * This method: - * 1. Stops the time loop after 1 iteration since the equation is steady state. - * 2. Resets, configures, and initializes the linear system. - * 3. Executes the stationary solve. + * - Stops the time loop after 1 iteration since the equation is steady state. + * - Resets, configures, and initializes the linear system. + * - Executes the stationary solve. */ /*---------------------------------------------------------------------------*/ @@ -132,7 +130,8 @@ _doStationarySolve() void FemModule:: _getMaterialParameters() { - info() << "Get material parameters..."; + info() << "Get material parameters "; + lambda = options()->lambda(); qdot = options()->qdot(); @@ -144,7 +143,7 @@ _getMaterialParameters() for (const auto& bs : options()->materialProperty()) { CellGroup group = bs->volume(); Real value = bs->lambda(); - info() << "Lambda for group=" << group.name() << " v=" << value; + info() << "Lambda for group= " << group.name() << " v=" << value; ENUMERATE_ (Cell, icell, group) { Cell cell = *icell; @@ -160,22 +159,22 @@ _getMaterialParameters() * This method constructs the linear system by assembling the LHS matrix * and RHS vector, applying various boundary conditions and source terms. * - * - The RHS vector is initialized to zero before applying any conditions. - * - If a constant source term is specified (`qdot`), apply it to the RHS. - * - If Neumann BC are specified applied to the RHS. - * - If Dirichlet BC are specified apply to the LHS & RHS. - * - If manufactured source conditions is specified, apply to RHS. - * - If manufactured Dirichlet BC are specified apply to the LHS & RHS. + * Steps involved: + * 1. The RHS vector is initialized to zero before applying any conditions. + * 2. If a constant source term is specified (`qdot`), apply it to the RHS. + * 3. If Neumann BC are specified applied to the RHS. + * 4. If Dirichlet BC are specified apply to the LHS & RHS. + * 5. If manufactured source conditions is specified, apply to RHS. + * 6. If manufactured Dirichlet BC are specified apply to the LHS & RHS. */ /*---------------------------------------------------------------------------*/ void FemModule:: _assembleLinearOperator() { - info() << "Assembly of FEM linear operator "; + info() << "Assembly of FEM linear operator"; - // Temporary variable to keep values for the RHS part of the linear system - VariableDoFReal& rhs_values(m_linear_system.rhsVariable()); + VariableDoFReal& rhs_values(m_linear_system.rhsVariable()); // Temporary variable to keep values for the RHS rhs_values.fill(0.0); auto node_dof(m_dofs_on_nodes.nodeDoFConnectivityView()); @@ -248,9 +247,7 @@ _computeElementMatrixQUAD4(Cell cell) /*---------------------------------------------------------------------------*/ /** - * @brief Assembles the bilinear operator matrix for the FEM linear system. - * - * The method calls the right function for RHS assembly given as mesh type. + * @brief Calls the right function for LHS assembly given as mesh type. */ /*---------------------------------------------------------------------------*/ @@ -274,11 +271,10 @@ _assembleBilinearOperator() * @brief Assembles the bilinear operator matrix for the FEM linear system. * * The method performs the following steps: - * - Iterates over all cells in the mesh. - * - For each cell, retrieves the cell-specific constant `lambda`. - * - Computes the element matrix using the provided `compute_element_matrix` function. - * - Assembles the global matrix by adding the contributions from each cell's element - * matrix to the corresponding entries in the global matrix. + * 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 addingcontributions from each cell's element + * matrix to the corresponding entries in the global matrix. */ /*---------------------------------------------------------------------------*/ @@ -308,7 +304,6 @@ _assembleBilinear(const std::function(const Cell&)>& compute_e } } - /*---------------------------------------------------------------------------*/ /** * @brief Solves the linear system and updates the solution vector. @@ -364,12 +359,11 @@ _solve() void FemModule:: _validateResults() { - if (allNodes().size() < 200) { + if (allNodes().size() < 200) ENUMERATE_ (Node, inode, allNodes()) { Node node = *inode; info() << "u[" << node.localId() << "][" << node.uniqueId() << "] = " << m_u[node]; } - } String filename = options()->resultFile(); info() << "ValidateResultFile filename=" << filename; From 39475f00795e99f525fade0699353e287da9dbea Mon Sep 17 00:00:00 2001 From: Mohd Afeef Badri <52162083+mohd-afeef-badri@users.noreply.github.com> Date: Fri, 9 Aug 2024 17:43:54 +0200 Subject: [PATCH 22/32] volume of a tetrahedron --- femutils/ArcaneFemFunctions.h | 42 +++++++++++++++++++++++++++++++++++ 1 file changed, 42 insertions(+) diff --git a/femutils/ArcaneFemFunctions.h b/femutils/ArcaneFemFunctions.h index be39c12..3bf5767 100644 --- a/femutils/ArcaneFemFunctions.h +++ b/femutils/ArcaneFemFunctions.h @@ -34,6 +34,32 @@ class ArcaneFemFunctions public: + /*---------------------------------------------------------------------------*/ + /** + * @brief Computes the volume of a tetrahedra defined by four nodes. + * + * This method calculates the volume using the scalar triple product formula. + * We do the following: + * 1. get the four nodes + * 2. ge the vector representing the edges of tetrahedron + * 3. compute volume using scalar triple product + */ + /*---------------------------------------------------------------------------*/ + + static inline Real computeVolumeTetra4(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)]; + + Real3 v0 = vertex1 - vertex0; + Real3 v1 = vertex2 - vertex0; + Real3 v2 = vertex3 - vertex0; + + return std::abs(Arcane::math::dot(v0, Arcane::math::cross(v1, v2))) / 6.0; + } + /*---------------------------------------------------------------------------*/ /** * @brief Computes the area of a triangle defined by three nodes. @@ -63,6 +89,7 @@ class ArcaneFemFunctions * of the quadrilateral's vertices. */ /*---------------------------------------------------------------------------*/ + static inline Real computeAreaQuad4(Cell cell, const VariableNodeReal3& node_coord) { Real3 vertex0 = node_coord[cell.nodeId(0)]; @@ -321,6 +348,21 @@ class ArcaneFemFunctions } }; + /*---------------------------------------------------------------------------*/ + /** + * @brief Provides methods for finite element operations in 3D. + * + * This class includes static methods for calculating gradients of basis + * functions and integrals for P1 triangles in 3D finite element analysis. + */ + /*---------------------------------------------------------------------------*/ + class FeOperation3D + { + + public: + // TO DO + }; + /*---------------------------------------------------------------------------*/ /** * @brief Provides methods for applying boundary conditions in 2D FEM problems. From dfb1c7a6f010396b69be585bdfa6f0178a223cb4 Mon Sep 17 00:00:00 2001 From: Mohd Afeef Badri <52162083+mohd-afeef-badri@users.noreply.github.com> Date: Tue, 20 Aug 2024 10:29:11 +0200 Subject: [PATCH 23/32] minor correction --- fourier/FemModule.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/fourier/FemModule.h b/fourier/FemModule.h index 0df60c6..928efd7 100644 --- a/fourier/FemModule.h +++ b/fourier/FemModule.h @@ -5,9 +5,9 @@ // SPDX-License-Identifier: Apache-2.0 //----------------------------------------------------------------------------- /*---------------------------------------------------------------------------*/ -/* FemModule.cc (C) 2022-2024 */ +/* FemModule.h (C) 2022-2024 */ /* */ -/* Simple module to test simple FEM mechanism. */ +/* FemModule class definition. */ /*---------------------------------------------------------------------------*/ /*---------------------------------------------------------------------------*/ #ifndef FEMMODULES_H From a9e9ce575c1ac85c223498a8ebd580cd9775fddd Mon Sep 17 00:00:00 2001 From: Mohd Afeef Badri <52162083+mohd-afeef-badri@users.noreply.github.com> Date: Thu, 22 Aug 2024 10:05:39 +0200 Subject: [PATCH 24/32] correct naming --- acoustics/FemModule.cc | 12 ++++++------ acoustics/FemModule.h | 4 ++-- femutils/ArcaneFemFunctions.h | 20 ++++++++++---------- fourier/FemModule.cc | 12 ++++++------ fourier/FemModule.h | 4 ++-- 5 files changed, 26 insertions(+), 26 deletions(-) diff --git a/acoustics/FemModule.cc b/acoustics/FemModule.cc index 2df384f..7521f91 100644 --- a/acoustics/FemModule.cc +++ b/acoustics/FemModule.cc @@ -67,7 +67,7 @@ compute() * * This method follows a sequence of steps to solve FEM system: * 1. _getMaterialParameters() Retrieves material parameters via - * 2. _assembleBilinearOperatorTRIA3() Assembles the FEM matrix A + * 2. _assembleBilinearOperatorTria3() 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 @@ -78,7 +78,7 @@ void FemModule:: _doStationarySolve() { _getMaterialParameters(); - _assembleBilinearOperatorTRIA3(); + _assembleBilinearOperatorTria3(); _assembleLinearOperator(); _solve(); _validateResults(); @@ -137,10 +137,10 @@ _assembleLinearOperator() /*---------------------------------------------------------------------------*/ FixedMatrix<3, 3> FemModule:: -_computeElementMatrixTRIA3(Cell cell) +_computeElementMatrixTria3(Cell cell) { // step 1 - Real area = ArcaneFemFunctions::MeshOperation::computeAreaTriangle3(cell, m_node_coord); + Real area = ArcaneFemFunctions::MeshOperation::computeAreaTria3(cell, m_node_coord); // step 2 Real3x3 UV = ArcaneFemFunctions::FeOperation2D::computeUVTria3(cell, m_node_coord); @@ -163,14 +163,14 @@ _computeElementMatrixTRIA3(Cell cell) /*---------------------------------------------------------------------------*/ void FemModule:: -_assembleBilinearOperatorTRIA3() +_assembleBilinearOperatorTria3() { 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 = _computeElementMatrixTria3(cell); // element matrix Int32 n1_index = 0; for (Node node1 : cell.nodes()) { Int32 n2_index = 0; diff --git a/acoustics/FemModule.h b/acoustics/FemModule.h index 242d250..6e59720 100644 --- a/acoustics/FemModule.h +++ b/acoustics/FemModule.h @@ -76,12 +76,12 @@ class FemModule void _doStationarySolve(); void _getMaterialParameters(); - void _assembleBilinearOperatorTRIA3(); + void _assembleBilinearOperatorTria3(); void _solve(); void _assembleLinearOperator(); void _validateResults(); - FixedMatrix<3, 3> _computeElementMatrixTRIA3(Cell cell); + FixedMatrix<3, 3> _computeElementMatrixTria3(Cell cell); }; #endif diff --git a/femutils/ArcaneFemFunctions.h b/femutils/ArcaneFemFunctions.h index 3bf5767..c7ca1bb 100644 --- a/femutils/ArcaneFemFunctions.h +++ b/femutils/ArcaneFemFunctions.h @@ -12,7 +12,7 @@ using namespace Arcane; * The class provides methods organized into different nested classes for: * - MeshOperation: Mesh related operations. * - FeOperation2D: Finite element operations at element level. - * - BoundaryConditions2D: Boundary condition realated operations. + * - BoundaryConditions2D: Boundary condition related operations. */ /*---------------------------------------------------------------------------*/ @@ -70,7 +70,7 @@ class ArcaneFemFunctions */ /*---------------------------------------------------------------------------*/ - static inline Real computeAreaTriangle3(Cell cell, const VariableNodeReal3& node_coord) + static inline Real computeAreaTria3(Cell cell, const VariableNodeReal3& node_coord) { Real3 vertex0 = node_coord[cell.nodeId(0)]; Real3 vertex1 = node_coord[cell.nodeId(1)]; @@ -109,7 +109,7 @@ class ArcaneFemFunctions */ /*---------------------------------------------------------------------------*/ - static inline Real3 computeBaryCenterTriangle3(Cell cell, const VariableNodeReal3& node_coord) + static inline Real3 computeBaryCenterTria3(Cell cell, const VariableNodeReal3& node_coord) { Real3 vertex0 = node_coord[cell.nodeId(0)]; Real3 vertex1 = node_coord[cell.nodeId(1)]; @@ -129,7 +129,7 @@ class ArcaneFemFunctions */ /*---------------------------------------------------------------------------*/ - static inline Real computeEdgeLength2(Face face, const VariableNodeReal3& node_coord) + static inline Real computeLengthEdge2(Face face, const VariableNodeReal3& node_coord) { Real3 vertex0 = node_coord[face.nodeId(0)]; Real3 vertex1 = node_coord[face.nodeId(1)]; @@ -149,7 +149,7 @@ class ArcaneFemFunctions */ /*---------------------------------------------------------------------------*/ - static inline Real2 computeEdgeNormal2(Face face, const VariableNodeReal3& node_coord) + static inline Real2 computeNormalEdge2(Face face, const VariableNodeReal3& node_coord) { Real3 vertex0 = node_coord[face.nodeId(0)]; Real3 vertex1 = node_coord[face.nodeId(1)]; @@ -395,7 +395,7 @@ class ArcaneFemFunctions { ENUMERATE_ (Cell, icell, mesh->allCells()) { Cell cell = *icell; - Real area = ArcaneFemFunctions::MeshOperation::computeAreaTriangle3(cell, node_coord); + Real area = ArcaneFemFunctions::MeshOperation::computeAreaTria3(cell, node_coord); for (Node node : cell.nodes()) { if (node.isOwn()) rhs_values[node_dof.dofId(node, 0)] += qdot * area / cell.nbNode(); @@ -423,8 +423,8 @@ class ArcaneFemFunctions { ENUMERATE_ (Cell, icell, mesh->allCells()) { Cell cell = *icell; - Real area = ArcaneFemFunctions::MeshOperation::computeAreaTriangle3(cell, node_coord); - Real3 bcenter = ArcaneFemFunctions::MeshOperation::computeBaryCenterTriangle3(cell, node_coord); + Real area = ArcaneFemFunctions::MeshOperation::computeAreaTria3(cell, node_coord); + Real3 bcenter = ArcaneFemFunctions::MeshOperation::computeBaryCenterTria3(cell, node_coord); for (Node node : cell.nodes()) { if (node.isOwn()) @@ -472,8 +472,8 @@ class ArcaneFemFunctions ENUMERATE_ (Face, iface, group) { Face face = *iface; - Real length = ArcaneFemFunctions::MeshOperation::computeEdgeLength2(face, node_coord); - Real2 normal = ArcaneFemFunctions::MeshOperation::computeEdgeNormal2(face, node_coord); + Real length = ArcaneFemFunctions::MeshOperation::computeLengthEdge2(face, node_coord); + Real2 normal = ArcaneFemFunctions::MeshOperation::computeNormalEdge2(face, node_coord); for (Node node : iface->nodes()) { if (!node.isOwn()) diff --git a/fourier/FemModule.cc b/fourier/FemModule.cc index 67b3850..1fd63b4 100644 --- a/fourier/FemModule.cc +++ b/fourier/FemModule.cc @@ -221,9 +221,9 @@ _assembleLinearOperator() /*---------------------------------------------------------------------------*/ FixedMatrix<3, 3> FemModule:: -_computeElementMatrixTRIA3(Cell cell) +_computeElementMatrixTria3(Cell cell) { - Real area = ArcaneFemFunctions::MeshOperation::computeAreaTriangle3(cell, m_node_coord); + Real area = ArcaneFemFunctions::MeshOperation::computeAreaTria3(cell, m_node_coord); Real3 dxU = ArcaneFemFunctions::FeOperation2D::computeGradientXTria3(cell, m_node_coord); Real3 dyU = ArcaneFemFunctions::FeOperation2D::computeGradientYTria3(cell, m_node_coord); @@ -235,7 +235,7 @@ _computeElementMatrixTRIA3(Cell cell) /*---------------------------------------------------------------------------*/ FixedMatrix<4, 4> FemModule:: -_computeElementMatrixQUAD4(Cell cell) +_computeElementMatrixQuad4(Cell cell) { Real area = ArcaneFemFunctions::MeshOperation::computeAreaQuad4(cell, m_node_coord); @@ -256,11 +256,11 @@ _assembleBilinearOperator() { if (options()->meshType == "QUAD4") _assembleBilinear<4>([this](const Cell& cell) { - return _computeElementMatrixQUAD4(cell); + return _computeElementMatrixQuad4(cell); }); else if (options()->meshType == "TRIA3") _assembleBilinear<3>([this](const Cell& cell) { - return _computeElementMatrixTRIA3(cell); + return _computeElementMatrixTria3(cell); }); else ARCANE_FATAL("Non supported meshType"); @@ -273,7 +273,7 @@ _assembleBilinearOperator() * 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 addingcontributions from each cell's element + * 3. Assembles global matrix by adding contributions from each cell's element * matrix to the corresponding entries in the global matrix. */ /*---------------------------------------------------------------------------*/ diff --git a/fourier/FemModule.h b/fourier/FemModule.h index 928efd7..218c1d9 100644 --- a/fourier/FemModule.h +++ b/fourier/FemModule.h @@ -86,8 +86,8 @@ class FemModule void _assembleLinearOperator(); void _validateResults(); - FixedMatrix<3, 3> _computeElementMatrixTRIA3(Cell cell); - FixedMatrix<4, 4> _computeElementMatrixQUAD4(Cell cell); + FixedMatrix<3, 3> _computeElementMatrixTria3(Cell cell); + FixedMatrix<4, 4> _computeElementMatrixQuad4(Cell cell); IBinaryMathFunctor* m_manufactured_dirichlet = nullptr; IBinaryMathFunctor* m_manufactured_source = nullptr; From 2ab30c2bacdd91d15a5ed6dabf1c7a1eb9499f40 Mon Sep 17 00:00:00 2001 From: Mohd Afeef Badri <52162083+mohd-afeef-badri@users.noreply.github.com> Date: Thu, 22 Aug 2024 13:43:34 +0200 Subject: [PATCH 25/32] 3D problems with ArcaneFEMFunctions --- acoustics/FemModule.h | 2 +- femutils/ArcaneFemFunctions.h | 155 +++++++- laplace/Fem.axl | 76 +++- laplace/FemModule.cc | 717 +++++++--------------------------- laplace/FemModule.h | 89 +++++ 5 files changed, 449 insertions(+), 590 deletions(-) create mode 100644 laplace/FemModule.h diff --git a/acoustics/FemModule.h b/acoustics/FemModule.h index 6e59720..f03f917 100644 --- a/acoustics/FemModule.h +++ b/acoustics/FemModule.h @@ -7,7 +7,7 @@ /*---------------------------------------------------------------------------*/ /* FemModule.h (C) 2022-2024 */ /* */ -/* Simple module to test simple FEM mechanism. */ +/* FemModule class definition. */ /*---------------------------------------------------------------------------*/ /*---------------------------------------------------------------------------*/ #ifndef FEMMODULES_H diff --git a/femutils/ArcaneFemFunctions.h b/femutils/ArcaneFemFunctions.h index c7ca1bb..36036dc 100644 --- a/femutils/ArcaneFemFunctions.h +++ b/femutils/ArcaneFemFunctions.h @@ -171,6 +171,14 @@ class ArcaneFemFunctions * * This class includes static methods for calculating gradients of basis * functions and integrals for P1 triangles in 2D finite element analysis. + * + * Reference Tri3 and Quad4 element in Arcane + * + * 0 o 1 o . . . . o 0 + * . . . . + * . . . . + * . . . . + * 1 o . . o 2 2 o . . . . o 3 */ /*---------------------------------------------------------------------------*/ class FeOperation2D @@ -354,13 +362,156 @@ class ArcaneFemFunctions * * This class includes static methods for calculating gradients of basis * functions and integrals for P1 triangles in 3D finite element analysis. + * + * Reference Tetra4 element in Arcane + * + * 3 o + * /|\ + * / | \ + * / | \ + * / o 2 \ + * / . . \ + * 0 o-----------o 1 */ /*---------------------------------------------------------------------------*/ class FeOperation3D { - public: - // TO DO + /*-------------------------------------------------------------------------*/ + /** + * @brief Computes the X gradients of basis functions N for P1 Tetrahedron. + * + * This method calculates gradient operator ∂/∂x of Ni for a given P1 + * cell with i = 1,..,4 for the four shape function Ni hence output + * is a vector of size 4 + * + * ∂N/∂x = [ ∂N1/∂x ∂N2/∂x ∂N3/∂x ∂N4/∂x ] + * + * ∂N/∂x = 1/(6V) [ b0 b1 b2 b3 ] + * + * where: + * b0 = (m1.y * (m3.z - m2.z) + m2.y * (m1.z - m3.z) + m3.y * (m2.z - m1.z)), + * b1 = (m0.y * (m2.z - m3.z) + m2.y * (m3.z - m0.z) + m3.y * (m0.z - m2.z)), + * b2 = (m0.y * (m3.z - m1.z) + m1.y * (m0.z - m3.z) + m3.y * (m1.z - m0.z)), + * b3 = (m0.y * (m1.z - m2.z) + m1.y * (m2.z - m0.z) + m2.y * (m0.z - m1.z)). + * + */ + /*-------------------------------------------------------------------------*/ + + static inline Real4 computeGradientXTetra4(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)]; + + Real3 v0 = vertex1 - vertex0; + Real3 v1 = vertex2 - vertex0; + Real3 v2 = vertex3 - vertex0; + + // 6 x Volume of tetrahedron + Real V6 = std::abs(Arcane::math::dot(v0, Arcane::math::cross(v1, v2))); + + Real4 dx; + + dx[0] = (vertex1.y * (vertex3.z - vertex2.z) + vertex2.y * (vertex1.z - vertex3.z) + vertex3.y * (vertex2.z - vertex1.z))/V6; + dx[1] = (vertex0.y * (vertex2.z - vertex3.z) + vertex2.y * (vertex3.z - vertex0.z) + vertex3.y * (vertex0.z - vertex2.z))/V6; + dx[2] = (vertex0.y * (vertex3.z - vertex1.z) + vertex1.y * (vertex0.z - vertex3.z) + vertex3.y * (vertex1.z - vertex0.z))/V6; + dx[3] = (vertex0.y * (vertex1.z - vertex2.z) + vertex1.y * (vertex2.z - vertex0.z) + vertex2.y * (vertex0.z - vertex1.z))/V6; + + return dx; + } + + /*-------------------------------------------------------------------------*/ +/** + * @brief Computes the Y gradients of basis functions N for P1 Tetrahedron. + * + * This method calculates gradient operator ∂/∂y of Ni for a given P1 + * cell with i = 1,..,4 for the four shape functions Ni, hence the output + * is a vector of size 4. + * + * ∂N/∂y = [ ∂N1/∂y ∂N2/∂y ∂N3/∂y ∂N4/∂y ] + * + * ∂N/∂y = 1/(6V) [ c0 c1 c2 c3 ] + * + * where: + * c0 = (m1.z * (m3.x - m2.x) + m2.z * (m1.x - m3.x) + m3.z * (m2.x - m1.x)), + * c1 = (m0.z * (m2.x - m3.x) + m2.z * (m3.x - m0.x) + m3.z * (m0.x - m2.x)), + * c2 = (m0.z * (m3.x - m1.x) + m1.z * (m0.x - m3.x) + m3.z * (m1.x - m0.x)), + * c3 = (m0.z * (m1.x - m2.x) + m1.z * (m2.x - m0.x) + m2.z * (m0.x - m1.x)). + * + */ +/*-------------------------------------------------------------------------*/ + +static inline Real4 computeGradientYTetra4(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)]; + + Real3 v0 = vertex1 - vertex0; + Real3 v1 = vertex2 - vertex0; + Real3 v2 = vertex3 - vertex0; + + // 6 x Volume of tetrahedron + Real V6 = std::abs(Arcane::math::dot(v0, Arcane::math::cross(v1, v2))); + + Real4 dy; + + dy[0] = (vertex1.z * (vertex3.x - vertex2.x) + vertex2.z * (vertex1.x - vertex3.x) + vertex3.z * (vertex2.x - vertex1.x))/V6; + dy[1] = (vertex0.z * (vertex2.x - vertex3.x) + vertex2.z * (vertex3.x - vertex0.x) + vertex3.z * (vertex0.x - vertex2.x))/V6; + dy[2] = (vertex0.z * (vertex3.x - vertex1.x) + vertex1.z * (vertex0.x - vertex3.x) + vertex3.z * (vertex1.x - vertex0.x))/V6; + dy[3] = (vertex0.z * (vertex1.x - vertex2.x) + vertex1.z * (vertex2.x - vertex0.x) + vertex2.z * (vertex0.x - vertex1.x))/V6; + + return dy; +} + +/*-------------------------------------------------------------------------*/ +/** + * @brief Computes the Z gradients of basis functions N for P1 Tetrahedron. + * + * This method calculates gradient operator ∂/∂z of Ni for a given P1 + * cell with i = 1,..,4 for the four shape functions Ni, hence the output + * is a vector of size 4. + * + * ∂N/∂z = [ ∂N1/∂z ∂N2/∂z ∂N3/∂z ∂N4/∂z ] + * + * ∂N/∂z = 1/(6V) [ d0 d1 d2 d3 ] + * + * where: + * d0 = (m1.x * (m3.y - m2.y) + m2.x * (m1.y - m3.y) + m3.x * (m2.y - m1.y)), + * d1 = (m0.x * (m2.y - m3.y) + m2.x * (m3.y - m0.y) + m3.x * (m0.y - m2.y)), + * d2 = (m0.x * (m3.y - m1.y) + m1.x * (m0.y - m3.y) + m3.x * (m1.y - m0.y)), + * d3 = (m0.x * (m1.y - m2.y) + m1.x * (m2.y - m0.y) + m2.x * (m0.y - m1.y)). + * + */ +/*-------------------------------------------------------------------------*/ + +static inline Real4 computeGradientZTetra4(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)]; + + Real3 v0 = vertex1 - vertex0; + Real3 v1 = vertex2 - vertex0; + Real3 v2 = vertex3 - vertex0; + + // 6 x Volume of tetrahedron + Real V6 = std::abs(Arcane::math::dot(v0, Arcane::math::cross(v1, v2))); + + Real4 dz; + + dz[0] = (vertex1.x * (vertex3.y - vertex2.y) + vertex2.x * (vertex1.y - vertex3.y) + vertex3.x * (vertex2.y - vertex1.y))/V6; + dz[1] = (vertex0.x * (vertex2.y - vertex3.y) + vertex2.x * (vertex3.y - vertex0.y) + vertex3.x * (vertex0.y - vertex2.y))/V6; + dz[2] = (vertex0.x * (vertex3.y - vertex1.y) + vertex1.x * (vertex0.y - vertex3.y) + vertex3.x * (vertex1.y - vertex0.y))/V6; + dz[3] = (vertex0.x * (vertex1.y - vertex2.y) + vertex1.x * (vertex2.y - vertex0.y) + vertex2.x * (vertex0.y - vertex1.y))/V6; + + return dz; +} + }; /*---------------------------------------------------------------------------*/ diff --git a/laplace/Fem.axl b/laplace/Fem.axl index 7527427..d06e384 100644 --- a/laplace/Fem.axl +++ b/laplace/Fem.axl @@ -34,58 +34,90 @@ - - + Dirichlet boundary condition - + - FaceGroup on which to apply these boundary condition + Function for Dirichlet boundary condition - - + + - Value of the boundary condition + Method via which Dirichlet boundary condition is imposed + + + + + Penalty value for enforcing Dirichlet condition + + + + + Function for manufactured source term condition - - + - Neumann boundary condition + Dirichlet boundary condition FaceGroup on which to apply these boundary condition - + Value of the boundary condition - + - Value of the Neumann load in x + Method via which Dirichlet boundary condition is imposed - + - Value of the Neumann load in y + Penalty value for enforcing Dirichlet condition + + + Neumann boundary condition + + + FaceGroup on which to apply the boundary condition + + + + Value of the boundary condition + + + + Neumann load value in x-direction + + + + Neumann load value in y-direction + + + + + + + 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. + * + * This method initializes degrees of freedom (DoFs) on nodes. + */ /*---------------------------------------------------------------------------*/ -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 ElementNodes; - - DoFLinearSystem m_linear_system; - IItemFamily* m_dof_family = nullptr; - FemDoFsOnNodes m_dofs_on_nodes; - - private: - - void _doStationarySolve(); - void _getMaterialParameters(); - void _assembleBilinearOperatorTRIA3(); - void _assembleBilinearOperatorTETRA4(); - void _solve(); - void _initBoundaryconditions(); - void _assembleLinearOperator(); - void _applyDirichletBoundaryConditions(); - void _checkResultFile(); - FixedMatrix<3, 3> _computeElementMatrixTRIA3(Cell cell); - FixedMatrix<4, 4> _computeElementMatrixTETRA4(Cell cell); - Real _computeAreaTriangle3(Cell cell); - Real _computeAreaTetra4(Cell cell); - Real _computeEdgeLength2(Face face); - Real2 _computeEdgeNormal2(Face face); - -}; + 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. + * + * This method: + * 1. Stops the time loop after 1 iteration since the equation is steady state. + * 2. Resets, configures, and initializes the linear system. + * 3. Executes the stationary solve. + */ /*---------------------------------------------------------------------------*/ void FemModule:: @@ -112,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); } @@ -122,546 +65,173 @@ compute() } /*---------------------------------------------------------------------------*/ -/*---------------------------------------------------------------------------*/ - -void FemModule:: -startInit() -{ - info() << "Module Fem INIT"; - - m_dofs_on_nodes.initialize(mesh(), 1); - m_dof_family = m_dofs_on_nodes.dofFamily(); - - //_buildDoFOnNodes(); - //Int32 nb_node = allNodes().size(); - //m_k_matrix.resize(nb_node, nb_node); - //m_k_matrix.fill(0.0); - - //m_rhs_vector.resize(nb_node); - //m_rhs_vector.fill(0.0); - - // # init mesh - // # init behavior - // # init behavior on mesh entities - // # init BCs - _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) - if (options()->meshType == "TETRA4") - _assembleBilinearOperatorTETRA4(); - else - _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. + */ /*---------------------------------------------------------------------------*/ void FemModule:: _getMaterialParameters() { info() << "Get material parameters..."; - ElementNodes = 3.; - - if (options()->meshType == "TETRA4") - ElementNodes = 4.; -} - -/*---------------------------------------------------------------------------*/ -/*---------------------------------------------------------------------------*/ - -void FemModule:: -_initBoundaryconditions() -{ - info() << "Init 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_u[node] = value; - m_u_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_u[node] = value; - m_u_dirichlet[node] = true; - } - } } /*---------------------------------------------------------------------------*/ -// Assemble the FEM linear operator -// - This function enforces a Dirichlet boundary condition in a weak sense -// via the penalty method -// - The method also adds source term -// - Also adds external fluxes +/** + * @brief FEM linear operator for the current simulation step. + * + * This method constructs the linear system by assembling the LHS matrix + * and RHS vector, applying various boundary conditions and source terms. + * + * Steps involved: + * 1. The RHS vector is initialized to zero before applying any conditions. + * 2. If Neumann BC are specified applied to the RHS. + * 3. If Dirichlet BC/Point are specified apply to the LHS & RHS. + */ /*---------------------------------------------------------------------------*/ void FemModule:: _assembleLinearOperator() { - info() << "Assembly of FEM linear operator "; - info() << "Applying Dirichlet boundary condition via penalty method "; + info() << "Assembly of FEM linear operator"; - // Temporary variable to keep values for the RHS part of the linear system - VariableDoFReal& rhs_values(m_linear_system.rhsVariable()); + VariableDoFReal& rhs_values(m_linear_system.rhsVariable()); // Temporary variable to keep values for the RHS 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 "; + for (const auto& bs : options()->neumannBoundaryCondition()) + ArcaneFemFunctions::BoundaryConditions2D::applyNeumannToRhs(bs, node_dof, m_node_coord, rhs_values); - Real Penalty = options()->penalty(); // 1.0e30 is the default - - ENUMERATE_ (Node, inode, ownNodes()) { - NodeLocalId node_id = *inode; - if (m_u_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_u[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 + for (const auto& bs : options()->dirichletBoundaryCondition()) + ArcaneFemFunctions::BoundaryConditions2D::applyDirichletToLhsAndRhs(bs, node_dof, m_node_coord, m_linear_system, rhs_values); - ENUMERATE_ (Node, inode, ownNodes()) { - NodeLocalId node_id = *inode; - if (m_u_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_u[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 "; - - // TODO - }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 "; - - // TODO - }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 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_u_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_u_dirichlet[node]) && node.isOwn()) - rhs_values[node_dof.dofId(node, 0)] += (Normal.x*valueX + Normal.y*valueY) * length / 2.; - } - } - continue; - } - - 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_u_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_u_dirichlet[node]) && node.isOwn()) - rhs_values[node_dof.dofId(node, 0)] += (Normal.y*valueY) * length / 2.; - } - } - continue; - } - } -} - -/*---------------------------------------------------------------------------*/ -/*---------------------------------------------------------------------------*/ - -Real FemModule:: -_computeAreaTetra4(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)]; - Real3 m3 = m_node_coord[cell.nodeId(3)]; - - // Calculate vectors representing edges of the tetrahedron - Real3 v0 = m1 - m0; - Real3 v1 = m2 - m0; - Real3 v2 = m3 - m0; - - // Compute volume using scalar triple product - return std::abs(Arcane::math::dot(v0, Arcane::math::cross(v1, v2))) / 6.0; + for (const auto& bs : options()->dirichletPointCondition()) + ArcaneFemFunctions::BoundaryConditions2D::applyPointDirichletToLhsAndRhs(bs, node_dof, m_node_coord, m_linear_system, rhs_values); } /*---------------------------------------------------------------------------*/ +/** + * @brief Computes the element matrix for a tetrahedral element (P1 FE). + * + * This function calculates the integral of the expression: + * integral3D (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); + */ /*---------------------------------------------------------------------------*/ -Real FemModule:: -_computeAreaTriangle3(Cell cell) +FixedMatrix<4, 4> FemModule:: +_computeElementMatrixTetra4(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 volume = ArcaneFemFunctions::MeshOperation::computeVolumeTetra4(cell, m_node_coord); -/*---------------------------------------------------------------------------*/ -/*---------------------------------------------------------------------------*/ + 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); -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)); + return volume * (dxU ^ dxU) + volume * (dyU ^ dyU) + volume * (dzU ^ dzU); } /*---------------------------------------------------------------------------*/ +/** + * @brief Computes the element matrix for a triangular element (P1 FE). + * + * This function calculates the integral of the expression: + * integral2D (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); + */ /*---------------------------------------------------------------------------*/ -Real2 FemModule:: -_computeEdgeNormal2(Face face) +FixedMatrix<3, 3> FemModule:: +_computeElementMatrixTria3(Cell cell) { - 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; -} + Real area = ArcaneFemFunctions::MeshOperation::computeAreaTria3(cell, m_node_coord); -/*---------------------------------------------------------------------------*/ -/*---------------------------------------------------------------------------*/ + Real3 dxU = ArcaneFemFunctions::FeOperation2D::computeGradientXTria3(cell, m_node_coord); + Real3 dyU = ArcaneFemFunctions::FeOperation2D::computeGradientYTria3(cell, m_node_coord); -FixedMatrix<4, 4> FemModule:: -_computeElementMatrixTETRA4(Cell cell) -{ - // Get coordinates of the triangle element TETRA4 - //------------------------------------------------ - // 3 o - // /|\ - // / | \ - // / | \ - // / o 2 \ - // / . . \ - // o-----------o - // 0 1 - //------------------------------------------------ - Real3 m0 = m_node_coord[cell.nodeId(0)]; - Real3 m1 = m_node_coord[cell.nodeId(1)]; - Real3 m2 = m_node_coord[cell.nodeId(2)]; - Real3 m3 = m_node_coord[cell.nodeId(3)]; - - Real volume = _computeAreaTetra4(cell); - - // Compute gradients of shape functions - Real3 dPhi0 = Arcane::math::cross(m2 - m1, m1 - m3) ; - Real3 dPhi1 = Arcane::math::cross(m3 - m0, m0 - m2) ; - Real3 dPhi2 = Arcane::math::cross(m1 - m0, m0 - m3) ; - Real3 dPhi3 = Arcane::math::cross(m0 - m1, m1 - m2) ; - - // Construct the B-matrix - FixedMatrix<3, 4> b_matrix; - b_matrix(0, 0) = dPhi0.x; - b_matrix(1, 0) = dPhi0.y; - b_matrix(2, 0) = dPhi0.z; - - b_matrix(0, 1) = dPhi1.x; - b_matrix(1, 1) = dPhi1.y; - b_matrix(2, 1) = dPhi1.z; - - b_matrix(0, 2) = dPhi2.x; - b_matrix(1, 2) = dPhi2.y; - b_matrix(2, 2) = dPhi2.z; - - b_matrix(0, 3) = dPhi3.x; - b_matrix(1, 3) = dPhi3.y; - b_matrix(2, 3) = dPhi3.z; - - b_matrix.multInPlace(1.0 / (6.0 * volume)); - - // Compute the element matrix - FixedMatrix<4, 4> int_cdPi_dPj = matrixMultiplication(matrixTranspose(b_matrix), b_matrix); - int_cdPi_dPj.multInPlace(volume); - -/* - cout << " Ae \n" - << "\t" << int_cdPi_dPj(0,0)<<"\t"<< int_cdPi_dPj(0,1)<<"\t"<< int_cdPi_dPj(0,2)<<"\t"<< int_cdPi_dPj(0,3)<<"\n" - << "\t" << int_cdPi_dPj(1,0)<<"\t"<< int_cdPi_dPj(1,1)<<"\t"<< int_cdPi_dPj(1,2)<<"\t"<< int_cdPi_dPj(1,3)<<"\n" - << "\t" << int_cdPi_dPj(2,0)<<"\t"<< int_cdPi_dPj(2,1)<<"\t"<< int_cdPi_dPj(2,2)<<"\t"<< int_cdPi_dPj(2,3)<<"\n" - << "\t" << int_cdPi_dPj(3,0)<<"\t"<< int_cdPi_dPj(3,1)<<"\t"<< int_cdPi_dPj(3,2)<<"\t"<< int_cdPi_dPj(3,3)<<"\n" - << endl; -*/ - - return int_cdPi_dPj; -} - -/*---------------------------------------------------------------------------*/ -/*---------------------------------------------------------------------------*/ - -FixedMatrix<3, 3> FemModule:: -_computeElementMatrixTRIA3(Cell cell) -{ - // Get coordinates 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); - - FixedMatrix<2, 3> b_matrix; - b_matrix(0, 0) = dPhi0.x; - b_matrix(0, 1) = dPhi1.x; - b_matrix(0, 2) = dPhi2.x; - - 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)); - - FixedMatrix<3, 3> int_cdPi_dPj = matrixMultiplication(matrixTranspose(b_matrix), b_matrix); - int_cdPi_dPj.multInPlace(area); - - return int_cdPi_dPj; + return area * (dxU ^ dxU) + area * (dyU ^ dyU); } /*---------------------------------------------------------------------------*/ +/** + * @brief Calls the right function for LHS assembly given as mesh type. + */ /*---------------------------------------------------------------------------*/ void FemModule:: -_assembleBilinearOperatorTRIA3() +_assembleBilinearOperator() { - 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] - 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); - } - ++n2_index; - } - ++n1_index; - } - } + if (options()->meshType == "TETRA4") + _assembleBilinear<4>([this](const Cell& cell) { + return _computeElementMatrixTetra4(cell); + }); + 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:: -_assembleBilinearOperatorTETRA4() +_assembleBilinear(const std::function(const Cell&)>& compute_element_matrix) { auto node_dof(m_dofs_on_nodes.nodeDoFConnectivityView()); ENUMERATE_ (Cell, icell, allCells()) { Cell cell = *icell; - auto K_e = _computeElementMatrixTETRA4(cell); // element stiffness matrix - - // # assemble elementary matrix into the global one - // # elementary terms are positioned 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); } @@ -670,10 +240,17 @@ _assembleBilinearOperatorTETRA4() ++n1_index; } } - } /*---------------------------------------------------------------------------*/ +/** + * @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:: @@ -681,10 +258,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 @@ -697,31 +270,35 @@ _solve() } m_u.synchronize(); - - const bool do_print = (allNodes().size() < 200); - if (do_print) { - ENUMERATE_ (Node, inode, allNodes()) { - Node node = *inode; - info() << "u[" << node.localId() << "][" << node.uniqueId() << "] = " - << m_u[node]; - //info() << "u[]" << node.uniqueId() << " " - // << m_u[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() << "u[" << node.localId() << "][" << node.uniqueId() << "] = " << m_u[node]; + } + String filename = options()->resultFile(); - info() << "CheckResultFile filename=" << filename; - if (filename.empty()) - return; - const double epsilon = 1.0e-4; - checkNodeResultFile(traceMng(), filename, m_u, epsilon); + info() << "ValidateResultFile filename=" << filename; + + if (!filename.empty()) + checkNodeResultFile(traceMng(), filename, m_u, 1.0e-4); } /*---------------------------------------------------------------------------*/ @@ -730,4 +307,4 @@ _checkResultFile() ARCANE_REGISTER_MODULE_FEM(FemModule); /*---------------------------------------------------------------------------*/ -/*---------------------------------------------------------------------------*/ +/*---------------------------------------------------------------------------*/ \ No newline at end of file diff --git a/laplace/FemModule.h b/laplace/FemModule.h new file mode 100644 index 0000000..a423597 --- /dev/null +++ b/laplace/FemModule.h @@ -0,0 +1,89 @@ +// -*- 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.h (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 A module for finite element method. + * + * This class handles the initialization and computation for finite element + * method (FEM) simulations, providing methods to set up and solve linear + * systems, assemble FEM operators, and perform result checks. + */ +/*---------------------------------------------------------------------------*/ + +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); + } + + void startInit() override; //! Method called at the beginning of the simulation + void compute() override; //! Method called at each iteration + VersionInfo versionInfo() const override { return VersionInfo(1, 0, 0); } + + private: + + DoFLinearSystem m_linear_system; + IItemFamily* m_dof_family = nullptr; + FemDoFsOnNodes m_dofs_on_nodes; + + void _doStationarySolve(); + void _getMaterialParameters(); + void _assembleBilinearOperator(); + void _solve(); + void _assembleLinearOperator(); + void _validateResults(); + + FixedMatrix<3, 3> _computeElementMatrixTria3(Cell cell); + FixedMatrix<4, 4> _computeElementMatrixTetra4(Cell cell); + + template + void _assembleBilinear( const std::function(const Cell&)>& compute_element_matrix); +}; + +#endif \ No newline at end of file From 6727155ba48202c78542b9fae1369ed8eaa2d244 Mon Sep 17 00:00:00 2001 From: Mohd Afeef Badri <52162083+mohd-afeef-badri@users.noreply.github.com> Date: Thu, 22 Aug 2024 14:23:57 +0200 Subject: [PATCH 26/32] Documentation for Gradient --- doc/Gradient.md | 171 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 171 insertions(+) create mode 100644 doc/Gradient.md diff --git a/doc/Gradient.md b/doc/Gradient.md new file mode 100644 index 0000000..932e61a --- /dev/null +++ b/doc/Gradient.md @@ -0,0 +1,171 @@ +### How do we calculate gradients of shape functions for P1 Tetrahedron Finite Element ### + +For a tetrahedron with nodes we assume four nodes $m_0 , m_1. m_2, m_3$: + + + + $m_0 = (x_0, y_0, z_0)$, $m_1 = (x_1, y_1, z_1)$, $m_2 = (x_2, y_2, z_2)$, and $m_3 = (x_3, y_3, z_3)$, + +The P1 finite element shape functions $\Phi_i$ for a tetrahedron can be expressed as: + +$$ +\Phi_i(x, y, z) = \frac{1}{6V} \left( a_i + b_i x + c_i y + d_i z \right) +$$ + +Where $V$ is the volume of the tetrahedron. + +The gradients of the shape functions $\Phi_0$, $\Phi_1$, $\Phi_2$, $\Phi_3$ with respect to $x$, $y$, and $z$ can be computed as follows: + +## Gradient with Respect to $x$ ($\frac{\partial \Phi_i}{\partial x}$) + +The partial derivative with respect to $x$ is: + +$$ +\frac{\partial \Phi_i}{\partial x} = \frac{b_i}{6V} +$$ + +Here, the coefficients $b_i$ are: + +$$ +b_0 = \det \begin{pmatrix} +y_1 & z_1 \\ +y_2 & z_2 \\ +y_3 & z_3 +\end{pmatrix} +$$ + +$$ +b_1 = \det \begin{pmatrix} +y_0 & z_0 \\ +y_2 & z_2 \\ +y_3 & z_3 +\end{pmatrix} +$$ + +$$ +b_2 = \det \begin{pmatrix} +y_0 & z_0 \\ +y_1 & z_1 \\ +y_3 & z_3 +\end{pmatrix} +$$ + +$$ +b_3 = \det \begin{pmatrix} +y_0 & z_0 \\ +y_1 & z_1 \\ +y_2 & z_2 +\end{pmatrix} +$$ + +## Gradient with Respect to $y$ ($\frac{\partial \Phi_i}{\partial y}$) + +The partial derivative with respect to $y$ is: + +$$ +\frac{\partial \Phi_i}{\partial y} = \frac{c_i}{6V} +$$ + +Where the coefficients $c_i$ are: + +$$ +c_0 = \det \begin{pmatrix} +z_1 & x_1 \\ +z_2 & x_2 \\ +z_3 & x_3 +\end{pmatrix} +$$ + +$$ +c_1 = \det \begin{pmatrix} +z_0 & x_0 \\ +z_2 & x_2 \\ +z_3 & x_3 +\end{pmatrix} +$$ + +$$ +c_2 = \det \begin{pmatrix} +z_0 & x_0 \\ +z_1 & x_1 \\ +z_3 & x_3 +\end{pmatrix} +$$ + +$$ +c_3 = \det \begin{pmatrix} +z_0 & x_0 \\ +z_1 & x_1 \\ +z_2 & x_2 +\end{pmatrix} +$$ + +## Gradient with Respect to $z$ ($\frac{\partial \Phi_i}{\partial z}$) + +The partial derivative with respect to $z$ is: + +$$ +\frac{\partial \Phi_i}{\partial z} = \frac{d_i}{6V} +$$ + +Where the coefficients $d_i$ are: + +$$ +d_0 = \det \begin{pmatrix} +x_1 & y_1 \\ +x_2 & y_2 \\ +x_3 & y_3 +\end{pmatrix} +$$ + +$$ +d_1 = \det \begin{pmatrix} +x_0 & y_0 \\ +x_2 & y_2 \\ +x_3 & y_3 +\end{pmatrix} +$$ + +$$ +d_2 = \det \begin{pmatrix} +x_0 & y_0 \\ +x_1 & y_1 \\ +x_3 & y_3 +\end{pmatrix} +$$ + +$$ +d_3 = \det \begin{pmatrix} +x_0 & y_0 \\ +x_1 & y_1 \\ +x_2 & y_2 +\end{pmatrix} +$$ + +## Implementation in Code for Gradient in $x$ ($\frac{\partial \Phi_i}{\partial x}$): + +```cpp + static inline Real4 computeGradientXTetra4(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)]; + + Real3 v0 = vertex1 - vertex0; + Real3 v1 = vertex2 - vertex0; + Real3 v2 = vertex3 - vertex0; + + // 6 x Volume of tetrahedron + Real V6 = std::abs(Arcane::math::dot(v0, Arcane::math::cross(v1, v2))); + + Real4 dx; + + dx[0] = (vertex1.y * (vertex3.z - vertex2.z) + vertex2.y * (vertex1.z - vertex3.z) + vertex3.y * (vertex2.z - vertex1.z))/V6; + dx[1] = (vertex0.y * (vertex2.z - vertex3.z) + vertex2.y * (vertex3.z - vertex0.z) + vertex3.y * (vertex0.z - vertex2.z))/V6; + dx[2] = (vertex0.y * (vertex3.z - vertex1.z) + vertex1.y * (vertex0.z - vertex3.z) + vertex3.y * (vertex1.z - vertex0.z))/V6; + dx[3] = (vertex0.y * (vertex1.z - vertex2.z) + vertex1.y * (vertex2.z - vertex0.z) + vertex2.y * (vertex0.z - vertex1.z))/V6; + + return dx; + } +``` From 09feed69c935df7a80471717cc7db386ff3ffb8d Mon Sep 17 00:00:00 2001 From: Mohd Afeef Badri <52162083+mohd-afeef-badri@users.noreply.github.com> Date: Thu, 22 Aug 2024 15:04:51 +0200 Subject: [PATCH 27/32] clean up --- fourier/Fem.axl | 10 ---------- fourier/Test.conduction.arc | 8 ++++++-- laplace/Fem.axl | 13 ------------- 3 files changed, 6 insertions(+), 25 deletions(-) diff --git a/fourier/Fem.axl b/fourier/Fem.axl index f586d25..99d60d4 100644 --- a/fourier/Fem.axl +++ b/fourier/Fem.axl @@ -32,16 +32,6 @@ Type of mesh provided to the solver - - - Method via which Dirichlet boundary condition is imposed - - - - - Penalty value for enforcing Dirichlet condition - - 1.75 1e5 test1_results.txt - WeakPenalty - 1.e12 + Penalty + 1.e12 Cercle 50.0 + Penalty + 1.e12 Bas 5.0 + Penalty + 1.e12 Haut 21.0 diff --git a/laplace/Fem.axl b/laplace/Fem.axl index d06e384..a4fcc33 100644 --- a/laplace/Fem.axl +++ b/laplace/Fem.axl @@ -9,9 +9,6 @@ FEM variable u on nodes - - Boolean which is true if Dirichlet node - Node coordinates from Arcane variable @@ -23,16 +20,6 @@ Type of mesh provided to the solver - - - Method via which Dirichlet boundary condition is imposed - - - - - Penalty value for enforcing Dirichlet condition - - Date: Thu, 22 Aug 2024 17:14:36 +0200 Subject: [PATCH 28/32] electrostatic problems with ArcaneFEMFunctions --- electrostatics/Fem.axl | 64 +++- electrostatics/FemModule.cc | 580 ++++++++---------------------------- electrostatics/FemModule.h | 99 ++++++ fourier/FemModule.cc | 4 +- 4 files changed, 277 insertions(+), 470 deletions(-) create mode 100644 electrostatics/FemModule.h diff --git a/electrostatics/Fem.axl b/electrostatics/Fem.axl index d7acb5b..b78f2ab 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 0000000..09da1d0 --- /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 1fd63b4..01cc540 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`) */ /*---------------------------------------------------------------------------*/ From 080f03920a52277ff2c4aa822c3c0b24a4cc58a5 Mon Sep 17 00:00:00 2001 From: Mohd Afeef Badri <52162083+mohd-afeef-badri@users.noreply.github.com> Date: Thu, 22 Aug 2024 17:38:42 +0200 Subject: [PATCH 29/32] Create Derivatives_fe_field.md --- doc/Derivatives_fe_field.md | 94 +++++++++++++++++++++++++++++++++++++ 1 file changed, 94 insertions(+) create mode 100644 doc/Derivatives_fe_field.md diff --git a/doc/Derivatives_fe_field.md b/doc/Derivatives_fe_field.md new file mode 100644 index 0000000..336ec23 --- /dev/null +++ b/doc/Derivatives_fe_field.md @@ -0,0 +1,94 @@ +## Calculating Derivatives $dx(u)$ and $dy(u)$ Over a Triangular Finite Element + + + + +To calculate the derivatives of a function $u$ with respect to $x$ and $y$ ( denoted as $\frac{\partial u}{\partial x}$, $\frac{\partial u}{\partial y}$, or more simply $dx(u)$, $dy(u)$ ) over a triangular finite element, we follow these steps: + +### 1. Define the Triangle and the P1 Basis Functions + +In a P1 (linear) finite element, the function $u$ is assumed to be linear over each triangle. +Let the vertices of the triangle be $\mathbf{v}_1 = (x_1, y_1)$, $\mathbf{v}_2 = (x_2, y_2)$, and $\mathbf{v}_3 = (x_3, y_3)$, +with corresponding function values $u_1 = u(\mathbf{v}_1)$, $u_2 = u(\mathbf{v}_2)$, and $u_3 = u(\mathbf{v}_3)$. + +The function $u(x, y)$ can be expressed as: + +$$u(x, y) = a_0 + a_1 x + a_2 y$$ + +where $a_0$, $a_1$, and $a_2$ are coefficients to be determined. + +### 2. Set Up the System of Equations + +We solve for the coefficients $a_0$, $a_1$, and $a_2$ using the known values of $u$ at the triangle's vertices: + +$$ +\begin{aligned} +u_1 &= a_0 + a_1 x_1 + a_2 y_1 \\ +u_2 &= a_0 + a_1 x_2 + a_2 y_2 \\ +u_3 &= a_0 + a_1 x_3 + a_2 y_3 +\end{aligned} +$$ + +This system can be written in matrix form: + +$$ +\begin{pmatrix} +1 & x_1 & y_1 \\ +1 & x_2 & y_2 \\ +1 & x_3 & y_3 +\end{pmatrix} +\begin{pmatrix} +a_0 \\ +a_1 \\ +a_2 +\end{pmatrix}= +\begin{pmatrix} +u_1 \\ +u_2 \\ +u_3 +\end{pmatrix} +$$ + + +Solve this system to get $a_0$, $a_1$, and $a_2$. + +### 3. Calculate the Derivatives $dx(u)$ and $dy(u)$ + +- The derivative $\frac{\partial u}{\partial x}$ is simply $a_1$. This is because the partial derivative of $u(x, y)$ with respect to $x$ is: + +$$ +\frac{\partial u}{\partial x} = a_1 +$$ + +Therefore, $dx(u) = a_1$. + +- The derivative $\frac{\partial u}{\partial y}$ is simply $a_2$. This is because the partial derivative of $u(x, y)$ with respect to $y$ is: + +$$ +\frac{\partial u}{\partial y} = a_2 +$$ + +Therefore, $dy(u) = a_2$. + +### 4. Shortcut Using the Area Formula (2D case) + +For a 2D triangular element (i.e., $z = 0$), the coefficients $a_1$ and $a_2$ (which give the derivatives with respect to $x$ and $y$, respectively) can be directly computed using the area of the triangle and the coordinates of its vertices: + +$$ +a_1 = \frac{1}{2A} \left( u_1(y_2 - y_3) + u_2(y_3 - y_1) + u_3(y_1 - y_2) \right) +$$ + +$$ +a_2 = \frac{1}{2A} \left( u_1(x_3 - x_2) + u_2(x_1 - x_3) + u_3(x_2 - x_1) \right) +$$ + +where $A$ is the area of the triangle, given by: + +$$ +A = \frac{1}{2} \left| x_1(y_2 - y_3) + x_2(y_3 - y_1) + x_3(y_1 - y_2) \right| +$$ + +This method directly gives the derivatives $dx(u)$ and $dy(u)$ without solving a system of equations. + +## TO DO ## +ADD CODE EXAMPLE From 819bb6a02201a5c92441a567535f9943e5fcdad0 Mon Sep 17 00:00:00 2001 From: Mohd Afeef Badri <52162083+mohd-afeef-badri@users.noreply.github.com> Date: Fri, 23 Aug 2024 11:15:02 +0200 Subject: [PATCH 30/32] gradients of given function U for P1 triangles --- electrostatics/FemModule.cc | 76 +++++++++++++---------------------- electrostatics/FemModule.h | 2 - femutils/ArcaneFemFunctions.h | 34 ++++++++++++++++ 3 files changed, 63 insertions(+), 49 deletions(-) diff --git a/electrostatics/FemModule.cc b/electrostatics/FemModule.cc index e6e4fa6..303a84e 100644 --- a/electrostatics/FemModule.cc +++ b/electrostatics/FemModule.cc @@ -92,10 +92,9 @@ _doStationarySolve() /** * @brief Retrieves and sets the material parameters for the simulation. * - * This method initializes: - * - material properties: - * # charge density (`rho`) - * # freespace permittivity (`epsilon`) + * This method initializes material properties: + * - charge density (`rho`) + * - free space permittivity (`epsilon`) */ /*---------------------------------------------------------------------------*/ @@ -103,15 +102,23 @@ void FemModule:: _getMaterialParameters() { info() << "Get material parameters..."; - rho = options()->rho(); // charge density - epsilon = options()->epsilon(); // freespace permittivity + rho = options()->rho(); + epsilon = options()->epsilon(); } /*---------------------------------------------------------------------------*/ -// Assemble the FEM linear operator -// - This function enforces a Dirichlet boundary condition in a weak sense -// via the penalty method -// - The method also adds source term +/** + * @brief FEM linear operator for the current simulation step. + * + * This method constructs the linear system by assembling the LHS matrix + * and RHS vector, applying various boundary conditions and source terms. + * + * Steps involved: + * 1. The RHS vector is initialized to zero before applying any conditions. + * 2. If a constant source term is specified (`rho`), apply it to the RHS. + * 3. If Neumann BC are specified applied to the RHS. + * 4. If Dirichlet BC are specified apply to the LHS & RHS. + */ /*---------------------------------------------------------------------------*/ void FemModule:: @@ -124,8 +131,8 @@ _assembleLinearOperator() auto node_dof(m_dofs_on_nodes.nodeDoFConnectivityView()); - if (options()->rho.isPresent()){ - Real qdot = - rho /epsilon; + if (options()->rho.isPresent()) { + Real qdot = -rho / epsilon; ArcaneFemFunctions::BoundaryConditions2D::applyConstantSourceToRhs(qdot, mesh(), node_dof, m_node_coord, rhs_values); } @@ -137,6 +144,9 @@ _assembleLinearOperator() } /*---------------------------------------------------------------------------*/ +/** + * @brief Computes the Electric Field which is gradient of ∇ Phi = E + */ /*---------------------------------------------------------------------------*/ void FemModule:: @@ -144,38 +154,12 @@ _getE() { info() << "Postprocessing E"; - ENUMERATE_ (Cell, icell, allCells()) { - Cell cell = *icell; - Real2 DX = _computeDxDyOfRealTria3(cell); - m_E[cell].x = -DX.x ; - m_E[cell].y = -DX.y ; - m_E[cell].z = 0.0 ; - } - - m_E.synchronize(); -} - -/*---------------------------------------------------------------------------*/ -/*---------------------------------------------------------------------------*/ - -Real2 FemModule:: -_computeDxDyOfRealTria3(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)]; - - Real f0 = m_phi[cell.nodeId(0)]; - Real f1 = m_phi[cell.nodeId(1)]; - Real f2 = m_phi[cell.nodeId(2)]; - - Real detA = ( m0.x*(m1.y - m2.y) - m0.y*(m1.x - m2.x) + (m1.x*m2.y - m2.x*m1.y) ); - - Real2 DX; - DX.y = ( m0.x*(f1 - f2) - f0*(m1.x - m2.x) + (f2*m1.x - f1*m2.x) ) / detA; - DX.x = ( f0*(m1.y - m2.y) - m0.y*(f1 - f2) + (f1*m2.y - f2*m1.y) ) / detA; + ENUMERATE_ (Cell, icell, allCells()) { + Cell cell = *icell; + m_E[cell] = ArcaneFemFunctions::FeOperation2D::computeGradientTria3(cell, m_node_coord, m_phi); + } - return DX ; + m_E.synchronize(); } /*---------------------------------------------------------------------------*/ @@ -200,10 +184,9 @@ _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); + return area * (dxU ^ dxU) + area * (dyU ^ dyU); } - /*---------------------------------------------------------------------------*/ /** * @brief Calls the right function for LHS assembly given as mesh type. @@ -223,7 +206,6 @@ _assembleBilinearOperator() ARCANE_FATAL("Non supported meshType"); } - /*---------------------------------------------------------------------------*/ /** * @brief Assembles the bilinear operator matrix for the FEM linear system. @@ -326,4 +308,4 @@ _validateResults() ARCANE_REGISTER_MODULE_FEM(FemModule); /*---------------------------------------------------------------------------*/ -/*---------------------------------------------------------------------------*/ +/*---------------------------------------------------------------------------*/ \ No newline at end of file diff --git a/electrostatics/FemModule.h b/electrostatics/FemModule.h index 09da1d0..54654c5 100644 --- a/electrostatics/FemModule.h +++ b/electrostatics/FemModule.h @@ -90,8 +90,6 @@ class FemModule FixedMatrix<3, 3> _computeElementMatrixTria3(Cell cell); - Real2 _computeDxDyOfRealTria3(Cell cell); - template void _assembleBilinear( const std::function(const Cell&)>& compute_element_matrix); }; diff --git a/femutils/ArcaneFemFunctions.h b/femutils/ArcaneFemFunctions.h index 36036dc..cb08e24 100644 --- a/femutils/ArcaneFemFunctions.h +++ b/femutils/ArcaneFemFunctions.h @@ -186,6 +186,40 @@ class ArcaneFemFunctions public: + /*---------------------------------------------------------------------------*/ + /** + * @brief Computes the gradients of given function U for P1 triangles. + * + * This method calculates gradient operator ∇ Ui for a given P1 cell + * with i = 1,..,3 for the three values of Ui hence at cell nodes. + * The output is ∇ Ui is P0 (piece-wise constant) hence Real3 value + * per cell + * + * ∇ Ui = [ ∂U/∂x ∂U/∂y ∂U/∂z ] + * + * ∂U/∂x = ( u1*(y2 − y3) + u2*(y3 − y1) + u3*(y1 − y2) ) / (2*A) + * ∂U/∂y = ( u1*(x3 − x2) + u2*(x1 − x3) + u3*(x2 − x1) ) / (2*A) + * ∂U/∂z = 0 + * + * @note we can adapt the same for 3D by filling the third component + */ + /*---------------------------------------------------------------------------*/ + + static inline Real3 computeGradientTria3(Cell cell, const VariableNodeReal3& node_coord, const VariableNodeReal& u) + { + Real3 n0 = node_coord[cell.nodeId(0)]; + Real3 n1 = node_coord[cell.nodeId(1)]; + Real3 n2 = node_coord[cell.nodeId(2)]; + + Real u0 = u[cell.nodeId(0)]; + Real u1 = u[cell.nodeId(1)]; + Real u2 = u[cell.nodeId(2)]; + + Real A2 = ((n1.x - n0.x) * (n2.y - n0.y) - (n2.x - n0.x) * (n1.y - n0.y)); + + return Real3 ( (u0*(n1.y - n2.y) + u1*(n2.y - n0.y) + u2*(n0.y - n1.y)) / A2 , (u0*(n2.x - n1.x) + u1*(n0.x - n2.x) + u2*(n1.x - n0.x)) / A2 , 0); + } + /*---------------------------------------------------------------------------*/ /** * @brief Computes the gradients of basis functions N for P1 triangles. From 9a1f1c348dd600a2ca7dc849a21b2b086f8f6bd1 Mon Sep 17 00:00:00 2001 From: Mohd Afeef Badri <52162083+mohd-afeef-badri@users.noreply.github.com> Date: Fri, 23 Aug 2024 11:25:15 +0200 Subject: [PATCH 31/32] documentation --- doc/Derivatives_fe_field.md | 41 +++++++++++++++++++++++++------------ 1 file changed, 28 insertions(+), 13 deletions(-) diff --git a/doc/Derivatives_fe_field.md b/doc/Derivatives_fe_field.md index 336ec23..d2156b7 100644 --- a/doc/Derivatives_fe_field.md +++ b/doc/Derivatives_fe_field.md @@ -8,8 +8,8 @@ To calculate the derivatives of a function $u$ with respect to $x$ and $y$ ( den ### 1. Define the Triangle and the P1 Basis Functions In a P1 (linear) finite element, the function $u$ is assumed to be linear over each triangle. -Let the vertices of the triangle be $\mathbf{v}_1 = (x_1, y_1)$, $\mathbf{v}_2 = (x_2, y_2)$, and $\mathbf{v}_3 = (x_3, y_3)$, -with corresponding function values $u_1 = u(\mathbf{v}_1)$, $u_2 = u(\mathbf{v}_2)$, and $u_3 = u(\mathbf{v}_3)$. +Let the nodes of the triangle be $\mathbf{n}_0 = (x_0, y_0)$, $\mathbf{n}_1 = (x_1, y_1)$, and $\mathbf{n}_2 = (x_2, y_2)$, +with corresponding function values $u_0 = u(\mathbf{n}_0)$, $u_1 = u(\mathbf{n}_1)$, and $u_2 = u(\mathbf{n}_2)$. The function $u(x, y)$ can be expressed as: @@ -23,9 +23,9 @@ We solve for the coefficients $a_0$, $a_1$, and $a_2$ using the known values of $$ \begin{aligned} +u_0 &= a_0 + a_1 x_0 + a_2 y_0 \\ u_1 &= a_0 + a_1 x_1 + a_2 y_1 \\ -u_2 &= a_0 + a_1 x_2 + a_2 y_2 \\ -u_3 &= a_0 + a_1 x_3 + a_2 y_3 +u_2 &= a_0 + a_1 x_2 + a_2 y_2 \end{aligned} $$ @@ -33,9 +33,9 @@ This system can be written in matrix form: $$ \begin{pmatrix} +1 & x_0 & y_0 \\ 1 & x_1 & y_1 \\ -1 & x_2 & y_2 \\ -1 & x_3 & y_3 +1 & x_2 & y_2 \end{pmatrix} \begin{pmatrix} a_0 \\ @@ -43,9 +43,9 @@ a_1 \\ a_2 \end{pmatrix}= \begin{pmatrix} +u_0 \\ u_1 \\ -u_2 \\ -u_3 +u_2 \end{pmatrix} $$ @@ -75,20 +75,35 @@ Therefore, $dy(u) = a_2$. For a 2D triangular element (i.e., $z = 0$), the coefficients $a_1$ and $a_2$ (which give the derivatives with respect to $x$ and $y$, respectively) can be directly computed using the area of the triangle and the coordinates of its vertices: $$ -a_1 = \frac{1}{2A} \left( u_1(y_2 - y_3) + u_2(y_3 - y_1) + u_3(y_1 - y_2) \right) +a_1 = \frac{1}{2A} \left( u_0(y_1 - y_2) + u_1(y_2 - y_0) + u_2(y_0 - y_1) \right) $$ $$ -a_2 = \frac{1}{2A} \left( u_1(x_3 - x_2) + u_2(x_1 - x_3) + u_3(x_2 - x_1) \right) +a_2 = \frac{1}{2A} \left( u_0(x_2 - x_1) + u_1(x_0 - x_2) + u_2(x_1 - x_0) \right) $$ where $A$ is the area of the triangle, given by: $$ -A = \frac{1}{2} \left| x_1(y_2 - y_3) + x_2(y_3 - y_1) + x_3(y_1 - y_2) \right| +A = \frac{1}{2} \left| x_0(y_1 - y_2) + x_1(y_2 - y_0) + x_2(y_0 - y_1) \right| $$ This method directly gives the derivatives $dx(u)$ and $dy(u)$ without solving a system of equations. -## TO DO ## -ADD CODE EXAMPLE +## Implementation in Code: ## +```cpp + static inline Real3 computeGradientTria3(Cell cell, const VariableNodeReal3& node_coord, const VariableNodeReal& u) + { + Real3 n0 = node_coord[cell.nodeId(0)]; + Real3 n1 = node_coord[cell.nodeId(1)]; + Real3 n2 = node_coord[cell.nodeId(2)]; + + Real u0 = u[cell.nodeId(0)]; + Real u1 = u[cell.nodeId(1)]; + Real u2 = u[cell.nodeId(2)]; + + Real A2 = ((n1.x - n0.x) * (n2.y - n0.y) - (n2.x - n0.x) * (n1.y - n0.y)); + + return Real3 ( (u0*(n1.y - n2.y) + u1*(n2.y - n0.y) + u2*(n0.y - n1.y)) / A2 , (u0*(n2.x - n1.x) + u1*(n0.x - n2.x) + u2*(n1.x - n0.x)) / A2 , 0); + } +``` From ec117e3735c53854cef5b3873230c1a00cd98117 Mon Sep 17 00:00:00 2001 From: Mohd Afeef Badri <52162083+mohd-afeef-badri@users.noreply.github.com> Date: Fri, 23 Aug 2024 11:39:08 +0200 Subject: [PATCH 32/32] move to md for doc --- doc/ImposingDirichletConditions.html | 746 --------------------------- doc/ImposingDirichletConditions.md | 83 +++ 2 files changed, 83 insertions(+), 746 deletions(-) delete mode 100644 doc/ImposingDirichletConditions.html create mode 100644 doc/ImposingDirichletConditions.md diff --git a/doc/ImposingDirichletConditions.html b/doc/ImposingDirichletConditions.html deleted file mode 100644 index ff54378..0000000 --- a/doc/ImposingDirichletConditions.html +++ /dev/null @@ -1,746 +0,0 @@ - - - - - -Readme - -
-

How are Dirichlet Boundary conditions imposed

To enforce Dirichlet boundary conditions in FEM many methods are available, some of the common ones include:

MethodReferenceIn ArcaneFEM
Weak penalty methodBabuška, 1973aYES
Penalty methodBabuška, 1973aYES
Row elimination SOON
Row/Column elimination SOON
Lagrange multiplier methodBabuška, 1973bNO
Nitsche’s methodNitsche, 1971NO

The first four are made available, so that one can choose the most appropriate method for Dirichlet boundary condition implementation. The word 'appropriate' in the preceding sentence is tricky, it depends on the physics, size of the problem, conditioning of the problem, used linear-solver, etc. So choose wisely.

What remains common in these methods is these are applied after one assembles the FEM linear system. As such one could see these methods as a post-processing step to be applied to the assembled linear-systems Ax=b before they are moved to the solving stage. More precisely, these methods alter the vector b and/or matrix A.

Weak penalty method and Penalty method

These methods are not exact rather weak sense of applying Dirichlet boundary condition, however these are the most straightforward way to apply the Dirichlet boundary condition. These methods remain popular among practitioners and can be found in FEM packages such as MOOSE, FreeFEM, etc. Main benefits of these methods include, simplistic implementation, conserved non-zero structure (sparsity) of A as such matrix symmetry may be conserved, and avoiding insertion of zeros into A which is a tedious operation.

In a nutshell, to apply weak penalty method for a DOF i which corresponds to a boundary node in the mesh the following needs to be changed:

ai,i=ai,i+P and bi=gi×P

here, P is the penalty term and gi is the given Dirichlet value.

Similarly to apply penalty method for a DOF i which corresponds to a boundary node in the mesh the following needs to be changed:

ai,i=P and bi=gi×P

here, P is the penalty term and gi is the given Dirichlet value.

The logic is simple as P is a large large value P>>1, the diagonal contribution ai,i dominates the off-diagonal ones (since ai,j<<ai,i), this indeed means solution bi≈gi.

Generally, P is chosen to be large enough e.g, 1×1031 so that the resulting solution x=A−1b is precise enough for a double precision calculation. This works most often, however sometimes this may result into A becoming ill-conditioned, eventually this might lead in failed solves and overall accuracy losses. So the user should select and tune the value of P in this case.

Row elimination

To apply a Dirichlet condition gi on the Dirchlet DOF i via row elimination method two operations are performed.

  • First is on the matrix A we zero out the ith row and ith column of the matrix A and impose 1 at the diagonal ai,i=1. For a FEM problem, if ND is the total number of Dirichlet DOFs and if A has the size of NDOF, this method involves altering A:

    • ∀i=1,…,ND&∀j=1,…,NDOF -ai,j=0,aj,i=0,&ai,i=1

  • Second operation is on the vector b, for each Dirichlet condition to be imposed on the ith DOF, we substitute the Dirichlet value gi to ith component of the the vector b. Hence, this method involves:

    • ∀i=1,…,NDbi=gi

This method is more exact way of imposing Dirichlet boundary conditions , however the matrix A is not symmetric anymore (if it were before), the problem it is algorithmically more challenging to implement as A in comparison to penalty methods.

Row column elimination

To apply a Dirichlet condition gi on the Dirchlet DOF i via row column elimination method two operations are performed.

  • First is on the matrix A we zero out the ith row of the matrix A and impose 1 at the diagonal ai,i=1. For a FEM problem, if ND is the total number of Dirichlet DOFs and if A has the size of NDOF, this method involves altering A:

    • ∀i=1,…,ND&∀j=1,…,NDOF -ai,j=0,&ai,i=1

  • Second operation is on the vector b, for each Dirichlet condition to be imposed on the ith DOF, we subtract the ith column from the vector b. Hence, this method involves:

    • ∀i=1,…,ND&∀j=1,…,NDOF

      bj=bj−aj,i&bi=gi

This method is more exact way of imposing Dirichlet boundary conditions, however it is algorithmically more challenging to implement as A which is a sparse matrix is stored in compressed row format (CRS) then moving the columns to the right-hand side becomes expensive for large systems with many Dirichlet DOFs.

Referneces

  • Babuška, I., 1973. The finite element method with penalty. Mathematics of computation, 27(122), pp.221-228.
  • Babuška, I., 1973. The finite element method with Lagrangian multipliers. Numerische Mathematik, 20(3), pp.179-192.
  • Nitsche, J., 1971, July. Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind. In Abhandlungen aus dem mathematischen Seminar der Universität Hamburg (Vol. 36, No. 1, pp. 9-15). Berlin/Heidelberg: Springer-Verlag.

 

 

 

- - \ No newline at end of file diff --git a/doc/ImposingDirichletConditions.md b/doc/ImposingDirichletConditions.md new file mode 100644 index 0000000..37fc449 --- /dev/null +++ b/doc/ImposingDirichletConditions.md @@ -0,0 +1,83 @@ +### How are Dirichlet Boundary conditions imposed + +To enforce Dirichlet boundary conditions in FEM many methods are available, some of the common ones include: + +| **Method** | **Reference** | In ArcaneFEM | +| :------------------------: | :------------: | :----------: | +| Weak penalty method | Babuška, 1973a | YES | +| Penalty method | Babuška, 1973a | YES | +| Row elimination | | YES | +| Row/Column elimination | | YES | +| Lagrange multiplier method | Babuška, 1973b | NO | +| Nitsche’s method | Nitsche, 1971 | NO | + +The first four are made available, so that one can choose the most appropriate method for Dirichlet boundary condition implementation. The word 'appropriate' in the preceding sentence is tricky, it depends on the physics, size of the problem, conditioning of the problem, used linear-solver, etc. So choose wisely. + +What remains common in these methods is these are applied after one assembles the FEM linear system. As such one could see these methods as a post-processing step to be applied to the assembled linear-systems $\mathbf{Ax=b}$ before they are moved to the solving stage. More precisely, these methods alter the vector $\mathbf{b}$ and/or matrix $\mathbf{A}$. + +#### Weak penalty method and Penalty method #### + +These methods are not exact rather weak sense of applying Dirichlet boundary condition, however these are the most straightforward way to apply the Dirichlet boundary condition. These methods remain popular among practitioners and can be found in FEM packages such as MOOSE, FreeFEM, etc. Main benefits of these methods include, simplistic implementation, conserved non-zero structure (sparsity) of $\mathbf{A}$ as such matrix symmetry may be conserved, and avoiding insertion of zeros into $\mathbf{A}$ which is a tedious operation. + +In a nutshell, to apply weak penalty method for a DOF $i$ which corresponds to a boundary node in the mesh the following needs to be changed: + +$a_{i,i} = a_{i,i} + \mathcal{P}$ and $b_i = g_i \times \mathcal{P}$ + +here, $\mathcal{P}$ is the penalty term and $g_i$ is the given Dirichlet value. + +Similarly to apply penalty method for a DOF $i$ which corresponds to a boundary node in the mesh the following needs to be changed: + +$a_{i,i} = \mathcal{P}$ and $b_i = g_i \times \mathcal{P}$ + +here, $\mathcal{P}$ is the penalty term and $g_i$ is the given Dirichlet value. + +The logic is simple as $\mathcal{P}$ is a large large value $\mathcal{P}>>1$, the diagonal contribution $a_{i,i}$ dominates the off-diagonal ones (since $a_{i,j} << a_{i,i}$), this indeed means solution $b_i \approx g_i$. + +Generally, $\mathcal{P}$ is chosen to be large enough e.g, $1\times10^{31}$ so that the resulting solution $\mathbf{x=A^{-1}b}$ is precise enough for a double precision calculation. This works most often, however sometimes this may result into $\mathbf{A}$ becoming ill-conditioned, eventually this might lead in failed solves and overall accuracy losses. So the user should select and tune the value of $\mathcal{P}$ in this case. + +#### Row elimination #### + +To apply a Dirichlet condition $g_i$ on the Dirchlet DOF $i$ via row elimination method two operations are performed. + +- First is on the matrix $\mathbf{A}$ we zero out the $i$ th row and $i$ th column of the matrix $\mathbf{A}$ and impose $1$ at the diagonal $a_{i,i}=1$. For a FEM problem, if $N_{D}$ is the total number of Dirichlet DOFs and if $\mathbf{A}$ has the size of $N_{DOF}$, this method involves altering $\mathbf{A}$: + +$$\forall i=1, \ldots , N_D \quad ; \quad \forall j = 1 , \ldots , N_{DOF}$$ + +$$ a_{i,j}=0, \quad a_{j,i}=0, \quad a_{i,i} = 1$$ + +- Second operation is on the vector $\mathbf{b}$, for each Dirichlet condition to be imposed on the $i$ th DOF, we substitute the Dirichlet value $g_i$ to $i$ th component of the the vector $\mathbf{b}$. Hence, this method involves: + +$$\forall i=1, \ldots , N_D \quad \quad b_i = g_i$$ + +This method is more exact way of imposing Dirichlet boundary conditions , however the matrix $\mathbf{A}$ is not symmetric anymore (if it were before), the problem it is algorithmically more challenging to implement as $\mathbf{A}$ in comparison to penalty methods. + +#### Row column elimination #### + +To apply a Dirichlet condition $g_i$ on the Dirchlet DOF $i$ via row column elimination method two operations are performed. + +- First is on the matrix $\mathbf{A}$ we zero out the $i$ th row of the matrix $\mathbf{A}$ and impose $1$ at the diagonal $a_{i,i}=1$. For a FEM problem, if $N_{D}$ is the total number of Dirichlet DOFs and if $\mathbf{A}$ has the size of $N_{DOF}$, this method involves altering $\mathbf{A}$: + +$$\forall i=1,\ldots ,N_D \quad ; \quad \forall j = 1 ,\ldots, N_{DOF}$$ + +$$a_{i,j}=0, \quad a_{i,i} = 1$$ + +- Second operation is on the vector $\mathbf{b}$, for each Dirichlet condition to be imposed on the $i$ th DOF, we subtract the $i$ th column from the vector $\mathbf{b}$. Hence, this method involves: + +$$\forall i=1,\ldots ,N_D \quad ; \quad \forall j = 1 ,\ldots, N_{DOF}$$ + +$$b_j = b_j - a_{j,i} , \quad b_i = g_i$$ + +This method is more exact way of imposing Dirichlet boundary conditions, however it is algorithmically more challenging to implement as $\mathbf{A}$ which is a sparse matrix is stored in compressed row format (CRS) then moving the columns to the right-hand side becomes expensive for large systems with many Dirichlet DOFs. + +## Referneces ## + +- Babuška, I., 1973. The finite element method with penalty. *Mathematics of computation*, *27*(122), pp.221-228. + +- Babuška, I., 1973. The finite element method with Lagrangian multipliers. *Numerische Mathematik*, *20*(3), pp.179-192. + +- Nitsche, J., 1971, July. Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind. In *Abhandlungen aus dem mathematischen Seminar der Universität Hamburg* (Vol. 36, No. 1, pp. 9-15). Berlin/Heidelberg: Springer-Verlag. + + + + +