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