1+ // -*- tab-width: 2; indent-tabs-mode: nil; coding: utf-8-with-signature -*-
2+ // -----------------------------------------------------------------------------
3+ // Copyright 2000-2025 CEA (www.cea.fr) IFPEN (www.ifpenergiesnouvelles.com)
4+ // See the top-level COPYRIGHT file for details.
5+ // SPDX-License-Identifier: Apache-2.0
6+ // -----------------------------------------------------------------------------
7+ /* ---------------------------------------------------------------------------*/
8+ /* BodyForce.h (C) 2022-2025 */
9+ /* */
10+ /* Contains functions to compute and assemble body force contribution to RHS */
11+ /* ---------------------------------------------------------------------------*/
12+ /* ---------------------------------------------------------------------------*/
13+
14+ /* ---------------------------------------------------------------------------*/
15+ /* *
16+ * @brief Applies body force to the RHS vector of the linear system.
17+ *
18+ * This function computes the contribution of body forces to the RHS vector
19+ * of the linear system. It iterates over all cells in the mesh, calculates
20+ * the appropriate force contributions based on the element type and mesh
21+ * dimension, and updates the RHS vector accordingly.
22+ *
23+ * body force ∫∫∫ (𝐟.𝐯) with 𝐟 = (𝑓𝑥, 𝑓𝑦, 𝑓𝑧) = (f[0], f[1], f[2])
24+ *
25+ * @param rhs_values The variable representing the RHS vector to be updated.
26+ * @param node_dof The connectivity view mapping nodes to their corresponding
27+ * degrees of freedom (DoFs).
28+ *
29+ /*---------------------------------------------------------------------------*/
30+
31+ inline void FemModule::
32+ _applyBodyForceToRHS (VariableDoFReal& rhs_values, const IndexedNodeDoFConnectivityView& node_dof)
33+ {
34+
35+ // get bodyforce vector
36+ bool applyBodyForce = false ;
37+ const UniqueArray<String> f_string = options ()->f ();
38+ info () << " [ArcaneFem-Info] Applying Bodyforce "
39+ for (Int32 i = 0 ; i < f_string.size (); ++i) {
40+ f[i] = 0.0 ;
41+ if (f_string[i] != " NULL"
42+ applyBodyForce = true ;
43+ f[i] = std::stod (f_string[i].localstr ());
44+ }
45+ }
46+
47+ // no bodyforce to apply hence return
48+ if (!applyBodyForce)
49+ return ;
50+
51+ // apply bodyforce based on dimension and mesh type
52+ if (mesh ()->dimension () == 2 ) {
53+ if (m_hex_quad_mesh) {
54+ ENUMERATE_ (Cell, icell, allCells ()) {
55+ Cell cell = *icell;
56+
57+ constexpr Real gp[2 ] = { -M_SQRT1_3, M_SQRT1_3 }; // (ξ,η)
58+ constexpr Real weights[2 ] = { 1.0 , 1.0 };
59+
60+ for (Int32 ixi = 0 ; ixi < 2 ; ++ixi) {
61+ for (Int32 ieta = 0 ; ieta < 2 ; ++ieta) {
62+
63+ // set the coordinates of the Gauss point
64+ Real xi = gp[ixi]; // Get the ξ coordinate of the Gauss point
65+ Real eta = gp[ieta]; // Get the η coordinate of the Gauss point
66+
67+ // set the weight of the Gauss point
68+ Real weight = weights[ixi] * weights[ieta];
69+
70+ // get shape functions 𝐍 for Quad4: 𝐍(ξ,η) = [𝑁₁ 𝑁₂ 𝑁₃ 𝑁₄]
71+ RealVector<4 > N = ArcaneFemFunctions::FeOperation2D::computeShapeFunctionsQuad4 (xi, eta);
72+
73+ // get determinant of Jacobian
74+ const auto gp_info = ArcaneFemFunctions::FeOperation2D::computeGradientsAndJacobianQuad4 (cell, m_node_coord, xi, eta);
75+ const Real detJ = gp_info.det_j ;
76+
77+ // compute integration weight
78+ Real integration_weight = weight * detJ;
79+
80+ // Assemble RHS
81+ for (Int32 i = 0 ; i < 4 ; ++i) {
82+ Node node = cell.node (i);
83+ if (node.isOwn ()) {
84+ rhs_values[node_dof.dofId (node, 0 )] += N[i] * f[0 ] * integration_weight;
85+ rhs_values[node_dof.dofId (node, 1 )] += N[i] * f[1 ] * integration_weight;
86+ }
87+ }
88+ }
89+ }
90+ }
91+ }
92+ else {
93+ ENUMERATE_ (Cell, icell, allCells ()) {
94+ Cell cell = *icell;
95+ Real area = ArcaneFemFunctions::MeshOperation::computeAreaTria3 (cell, m_node_coord);
96+ for (Node node : cell.nodes ()) {
97+ if (node.isOwn ()) {
98+ rhs_values[node_dof.dofId (node, 0 )] += f[0 ] * area / 3 ;
99+ rhs_values[node_dof.dofId (node, 1 )] += f[1 ] * area / 3 ;
100+ }
101+ }
102+ }
103+ }
104+ }
105+ if (mesh ()->dimension () == 3 ) {
106+ if (m_hex_quad_mesh) {
107+ ENUMERATE_ (Cell, icell, allCells ()) {
108+ Cell cell = *icell;
109+
110+ constexpr Real gp[2 ] = { -M_SQRT1_3, M_SQRT1_3 }; // [-1/sqrt(3), 1/sqrt(3)]
111+ constexpr Real weights[2 ] = { 1.0 , 1.0 };
112+
113+ for (Int32 ixi = 0 ; ixi < 2 ; ++ixi) {
114+ for (Int32 ieta = 0 ; ieta < 2 ; ++ieta) {
115+ for (Int32 izeta = 0 ; izeta < 2 ; ++izeta) {
116+
117+ // set the coordinates of the Gauss point
118+ Real xi = gp[ixi];
119+ Real eta = gp[ieta];
120+ Real zeta = gp[izeta];
121+
122+ // set the weight of the Gauss point
123+ Real weight = weights[ixi] * weights[ieta] * weights[izeta];
124+
125+ // Get shape functions for Hexa8: 𝐍(ξ,η,ζ) = [𝑁₁ 𝑁₂ 𝑁₃ 𝑁₄ 𝑁₅ 𝑁₆ 𝑁₇ 𝑁₈]
126+ RealVector<8 > N = ArcaneFemFunctions::FeOperation3D::computeShapeFunctionsHexa8 (xi, eta, zeta);
127+
128+ // get determinant of Jacobian
129+ const auto gp_info = ArcaneFemFunctions::FeOperation3D::computeGradientsAndJacobianHexa8 (cell, m_node_coord, xi, eta, zeta);
130+ const Real detJ = gp_info.det_j ;
131+
132+ // compute integration weight
133+ Real integration_weight = detJ * weight;
134+
135+ // Assemble RHS
136+ for (Int32 i = 0 ; i < 8 ; ++i) {
137+ Node node = cell.node (i);
138+ if (node.isOwn ()) {
139+ rhs_values[node_dof.dofId (node, 0 )] += N[i] * f[0 ] * integration_weight;
140+ rhs_values[node_dof.dofId (node, 1 )] += N[i] * f[1 ] * integration_weight;
141+ rhs_values[node_dof.dofId (node, 2 )] += N[i] * f[2 ] * integration_weight;
142+ }
143+ }
144+ }
145+ }
146+ }
147+ }
148+ }
149+ else {
150+ ENUMERATE_ (Cell, icell, allCells ()) {
151+ Cell cell = *icell;
152+ Real volume = ArcaneFemFunctions::MeshOperation::computeVolumeTetra4 (cell, m_node_coord);
153+ for (Node node : cell.nodes ()) {
154+ if (node.isOwn ()) {
155+ rhs_values[node_dof.dofId (node, 0 )] += f[0 ] * volume / 4 ;
156+ rhs_values[node_dof.dofId (node, 1 )] += f[1 ] * volume / 4 ;
157+ rhs_values[node_dof.dofId (node, 2 )] += f[2 ] * volume / 4 ;
158+ }
159+ }
160+ }
161+ }
162+ }
163+ }
0 commit comments