BELFEM 0.9.3
Berkeley Lab Finite Element Framework
Loading...
Searching...
No Matches
cl_Gradient.hpp
Go to the documentation of this file.
1/*
2 * BELFEM -- The Berkeley Lab Finite Element Framework
3 * Copyright (c) 2026, The Regents of the University of California,
4 * through Lawrence Berkeley National Laboratory (subject to receipt of any required
5 * approvals from the U.S. Dept. of Energy). All rights reserved.
6 *
7 * Developers: Christian Messe, Gregory Giard
8 *
9 * See the top-level LICENSE file for the complete license and disclaimer.
10 */
11
12#ifndef CL_GRADIENT_HPP
13#define CL_GRADIENT_HPP
14
15#include "typedefs.hpp"
16#include "cl_Vector.hpp"
17#include "cl_Cell.hpp"
18#include "cl_Mesh.hpp"
19#include "cl_FEM_Kernel.hpp"
20#include "cl_FEM_Postprocessor.hpp"
21#include "fn_trans.hpp"
22
23namespace belfem
24{
25 namespace fem
26 {
33 class Gradient : public Postprocessor
34 {
35 const int mNumDimensions ;
36
37 // flag telling if we need the negative gradient
38 bool mFlipSign = false ;
39
40 // tells which dof manager is used, default: 0
41 const uint mDofMaganerIndex ;
42
43 Map< id_t, Block * > mBlocks ;
44 Cell< mesh::Node * > mNodes ;
47 Matrix< real > mX ; // node coordinates
48 Matrix< real > mP ; // polynomial vector
49 Matrix< real > mC ; // coefficient matrix
50 Matrix< real > mV ; // vandermonde matrix
51 Vector< int_t > mPivot ; // pivot for lapack
52 Vector< real > mG ; // gradient
53
54 Matrix< real > mJ ;
55 Matrix< real > mB ;
56 Vector< real > mPhi ;
57 Matrix< real > mElX ;
58
59 int mN ; // number of coefficients
60
61 uint mOrder = 0 ;
62
63 string mScalarField ;
64 Cell< string > mGradientFields ;
65
66 void
67 ( Gradient::*mFunComputePoly )( const Matrix< real > & aX );
68
69 void
70 ( Gradient::*mFunPoly2D )( const real x, const real y );
71
72 void
73 ( Gradient::*mFunPoly3D )( const real x, const real y, const real z );
74
75//------------------------------------------------------------------------------
76 public:
77//------------------------------------------------------------------------------
78
79 Gradient(
80 Kernel * aKernel,
81 const Vector< id_t > & aBlocksIDs,
82 const bool aFlipSign = false,
83 const uint aDofManagerIndex = 0 );
84
85 ~Gradient() override = default;
86
87 void
88 set_fields( const string & aScalarField, const string & aGradientField );
89
90 virtual void
91 process_node( mesh::Node * aNode );
92
93 void
94 set_order( const uint aOrder );
95
96 void
97 run() override;
98
99//------------------------------------------------------------------------------
100 private:
101//------------------------------------------------------------------------------
102
103 void
104 compute_poly( const Matrix< real > & aX );
105
106 void
107 compute_poly_2d( const Matrix< real > & aX );
108
109 void
110 compute_poly_3d( const Matrix< real > & aX );
111
112 void
113 poly1_2d( const real x, const real y );
114
115 void
116 poly2_2d( const real x, const real y );
117
118 void
119 poly3_2d( const real x, const real y );
120
121 void
122 poly4_2d( const real x, const real y );
123
124 void
125 poly1_3d( const real x, const real y, const real z );
126
127 void
128 poly2_3d( const real x, const real y, const real z );
129
130 void
131 poly3_3d( const real x, const real y, const real z );
132
133 void
134 poly4_3d( const real x, const real y, const real z );
135
136 uint
137 subprocess_node( mesh::Node * aNode );
138 };
139
140 inline void
141 Gradient::compute_poly( const Matrix< real > & aX )
142 {
143 (this->*mFunComputePoly)( aX );
144 }
145
146 inline void
147 Gradient::compute_poly_2d( const Matrix< real > & aX )
148 {
149 (this->*mFunPoly2D )( aX( 0, 0 ), aX( 0, 1 ) );
150 }
151
152 inline void
153 Gradient::compute_poly_3d( const Matrix< real > & aX )
154 {
155 (this->*mFunPoly3D )( aX( 0, 0 ), aX( 0, 1 ), aX( 0, 2 ) );
156 }
157
158 inline void
159 Gradient::poly1_2d( const real x, const real y )
160 {
161 mP( 0, 0 ) = 1.0 ;
162 mP( 1, 0 ) = x ;
163 mP( 2, 0 ) = y ;
164 }
165
166 inline void
167 Gradient::poly2_2d( const real x, const real y )
168 {
169 mP( 0, 0 ) = 1.0 ;
170 mP( 1, 0 ) = x ;
171 mP( 2, 0 ) = y ;
172 mP( 3, 0 ) = x*x ;
173 mP( 4, 0 ) = x*y ;
174 mP( 5, 0 ) = y*y ;
175 }
176
177 inline void
178 Gradient::poly3_2d( const real x, const real y )
179 {
180 mP( 0, 0 ) = 1.0 ;
181 mP( 1, 0 ) = x ;
182 mP( 2, 0 ) = y ;
183 mP( 3, 0 ) = x*x ;
184 mP( 4, 0 ) = x*y ;
185 mP( 5, 0 ) = y*y ;
186 mP( 6, 0 ) = mP( 3, 0 )*x ;
187 mP( 7, 0 ) = mP( 3, 0 )*y ;
188 mP( 8, 0 ) = x * mP( 5, 0 );
189 mP( 9, 0 ) = y * mP( 5, 0 );
190 }
191
192 inline void
193 Gradient::poly4_2d( const real x, const real y )
194 {
195 mP( 0, 0 ) = 1.0 ;
196 mP( 1, 0 ) = x ;
197 mP( 2, 0 ) = y ;
198 mP( 3, 0 ) = x*x ;
199 mP( 4, 0 ) = x*y ;
200 mP( 5, 0 ) = y*y ;
201 mP( 6, 0 ) = mP( 3, 0 )*x ; // x^3
202 mP( 7, 0 ) = mP( 3, 0 )*y ; // x^2 * y
203 mP( 8, 0 ) = x * mP( 5, 0 ); // x * y^2
204 mP( 9, 0 ) = y * mP( 5, 0 ); // y^3
205 mP( 10, 0 ) = mP( 3, 0 ) * mP( 3, 0 ); // x^4
206 mP( 11, 0 ) = mP( 6, 0 ) * y ; // x^3 * y
207 mP( 12, 0 ) = mP( 3, 0 ) * mP( 5, 0 ) ; // x^2*y^2
208 mP( 13, 0 ) = x * mP( 9, 0 ) ; // x * y^3
209 mP( 14, 0 ) = mP( 5, 0 ) * mP( 5, 0 ) ; // y^4
210 }
211
212 inline void
213 Gradient::poly1_3d( const real x, const real y, const real z )
214 {
215 mP( 0, 0 ) = 1.0 ;
216 mP( 1, 0 ) = x ;
217 mP( 2, 0 ) = y ;
218 mP( 3, 0 ) = z ;
219 }
220
221 inline void
222 Gradient::poly2_3d( const real x, const real y, const real z )
223 {
224 mP( 0, 0 ) = 1.0 ;
225 mP( 1, 0 ) = x ;
226 mP( 2, 0 ) = y ;
227 mP( 3, 0 ) = z ;
228 mP( 4, 0 ) = x*x ;
229 mP( 5, 0 ) = x*y ;
230 mP( 6, 0 ) = y*y ;
231 mP( 7, 0 ) = y*z ;
232 mP( 8, 0 ) = z*z ;
233 mP( 9, 0 ) = z*x ;
234 }
235
236 inline void
237 Gradient::poly3_3d( const real x, const real y, const real z )
238 {
239 mP( 0, 0 ) = 1.0 ;
240 mP( 1, 0 ) = x ;
241 mP( 2, 0 ) = y ;
242 mP( 3, 0 ) = z ;
243 mP( 4, 0 ) = x*x ;
244 mP( 5, 0 ) = x*y ;
245 mP( 6, 0 ) = y*y ;
246 mP( 7, 0 ) = y*z ;
247 mP( 8, 0 ) = z*z ;
248 mP( 9, 0 ) = z*x ;
249 mP( 10, 0 ) = mP( 4, 0 ) * x ; // x^3
250 mP( 11, 0 ) = mP( 4, 0 ) * y ; // x^2 * y
251 mP( 12, 0 ) = x * mP( 6, 0 ) ; // x * y^2 ;
252 mP( 13, 0 ) = y * mP( 6, 0 ) ; // y^3
253 mP( 14, 0 ) = z * mP( 6, 0 ) ; // y^2 * z
254 mP( 15, 0 ) = y * mP( 8, 0 ) ; // y * z^2
255 mP( 16, 0 ) = z * mP( 8, 0 ) ; // z^3
256 mP( 17, 0 ) = x * mP( 8, 0 ) ; // z^2 * x
257 mP( 18, 0 ) = mP( 4, 0 ) * z ; // x^2 * z
258 mP( 19, 0 ) = x * y * z ;
259 }
260
261 inline void
262 Gradient::poly4_3d( const real x, const real y, const real z )
263 {
264 mP( 0, 0 ) = 1.0 ;
265 mP( 1, 0 ) = x ;
266 mP( 2, 0 ) = y ;
267 mP( 3, 0 ) = z ;
268 mP( 4, 0 ) = x*x ;
269 mP( 5, 0 ) = x*y ;
270 mP( 6, 0 ) = y*y ;
271 mP( 7, 0 ) = y*z ;
272 mP( 8, 0 ) = z*z ;
273 mP( 9, 0 ) = z*x ;
274 mP( 10, 0 ) = mP( 4, 0 ) * x ; // x^3
275 mP( 11, 0 ) = mP( 4, 0 ) * y ; // x^2 * y
276 mP( 12, 0 ) = x * mP( 6, 0 ) ; // x * y^2 ;
277 mP( 13, 0 ) = y * mP( 6, 0 ) ; // y^3
278 mP( 14, 0 ) = z * mP( 6, 0 ) ; // y^2 * z
279 mP( 15, 0 ) = y * mP( 8, 0 ) ; // y * z^2
280 mP( 16, 0 ) = z * mP( 8, 0 ) ; // z^3
281 mP( 17, 0 ) = x * mP( 8, 0 ) ; // z^2 * x
282 mP( 18, 0 ) = mP( 4, 0 ) * z ; // x^2 * z
283 mP( 19, 0 ) = x * y * z ;
284 mP( 20, 0 ) = mP( 4, 0 ) * mP( 4, 0 ) ; // x^4
285 mP( 21, 0 ) = mP( 4, 0 ) * mP( 5, 0 ) ; // x^3 * y
286 mP( 22, 0 ) = mP( 4, 0 ) * mP( 6, 0 ) ; // x^2 * y^2
287 mP( 23, 0 ) = mP( 5, 0 ) * mP( 6, 0 ) ; // x * y^3
288 mP( 24, 0 ) = mP( 6, 0 ) * mP( 6, 0 ) ; // y^4
289 mP( 25, 0 ) = mP( 6, 0 ) * mP( 7, 0 ) ; // y^3 * z
290 mP( 26, 0 ) = mP( 6, 0 ) * mP( 8, 0 ) ; // y^2 * z^2
291 mP( 27, 0 ) = mP( 7, 0 ) * mP( 8, 0 ) ; // y * z^3
292 mP( 28, 0 ) = mP( 8, 0 ) * mP( 8, 0 ) ; // z^4
293 mP( 29, 0 ) = mP( 9, 0 ) * mP( 8, 0 ) ; // z^3 * x
294 mP( 30, 0 ) = mP( 4, 0 ) * mP( 8, 0 ) ; // x^2 * z^2
295 mP( 31, 0 ) = mP( 9, 0 ) * mP( 4, 0 ) ; // z * x^3
296 mP( 32, 0 ) = mP( 4, 0 ) * mP( 7, 0 ) ; // x^2 * y * z
297 mP( 33, 0 ) = mP( 6, 0 ) * mP( 9, 0 ) ; // x * y^2 * z
298 mP( 34, 0 ) = mP( 5, 0 ) * mP( 8, 0 ) ; // x * y * z^2
299 }
300 }
301}
302#endif //CL_GRADIENT_HPP
cl_Cell.hpp
cl_FEM_Kernel.hpp
cl_FEM_Postprocessor.hpp
cl_Mesh.hpp
cl_Vector.hpp
belfem::Cell
A bounds-checked wrapper around std::vector.
Definition cl_Cell.hpp:45
belfem::Map
Hash map (unordered key-value).
Definition cl_Map.hpp:75
belfem::fem::Gradient::set_fields
void set_fields(const string &aScalarField, const string &aGradientField)
Definition cl_Gradient.cpp:86
belfem::fem::Gradient::Gradient
Gradient(Kernel *aKernel, const Vector< id_t > &aBlocksIDs, const bool aFlipSign=false, const uint aDofManagerIndex=0)
Definition cl_Gradient.cpp:19
belfem::fem::Gradient::~Gradient
~Gradient() override=default
belfem::fem::Gradient::run
void run() override
Definition cl_Gradient.cpp:278
belfem::fem::Gradient::set_order
void set_order(const uint aOrder)
Definition cl_Gradient.cpp:216
belfem::fem::Gradient::process_node
virtual void process_node(mesh::Node *aNode)
Definition cl_Gradient.cpp:107
belfem::fem::Kernel
Top-level orchestrator.
Definition cl_FEM_Kernel.hpp:56
belfem::fem::Postprocessor::Postprocessor
Postprocessor(Kernel *aKernel, DofManager *aField=nullptr)
Collective.
Definition cl_FEM_Postprocessor.cpp:27
fn_trans.hpp
Matrix transpose.
belfem::fem
Definition cl_IFB_LINE3.hpp:21
belfem
USER GUIDES:
Definition cl_Capacitor.cpp:16
belfem::uint
unsigned int uint
Definition typedefs.hpp:30
belfem::real
double real
Definition typedefs.hpp:36
corc.tests.test_curve_frame.x
float x
Definition test_curve_frame.py:37
mesh.Node
Definition Node.py:1
rutherford.tests.test_stadium.y
y
Definition test_stadium.py:44