BELFEM 0.9.3
Berkeley Lab Finite Element Framework
Loading...
Searching...
No Matches
cl_FEM_Postprocessor.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, through
4 * 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 BELFEM_CL_FEM_POSTPROCESSOR_HPP
13#define BELFEM_CL_FEM_POSTPROCESSOR_HPP
14
15#include "typedefs.hpp"
16#include "cl_Mesh.hpp"
17#include "cl_FEM_Kernel.hpp"
18
19namespace belfem
20{
21 namespace fem
22 {
24 {
25
26 protected:
30 const uint mOrder ;
31 private:
32
33 // block IDs used on this post processor
34 Cell< id_t > mBlockIDs ;
35
36 Cell< index_t > mMyNodeIndices ;
37 Cell< index_t > mMyOwnedNodeIndices ;
38 Cell< index_t > mMyElementIndices ;
39 Cell< index_t > mMyEdgeIndices ;
40 Cell< index_t > mMyFaceIndices ;
41
42 Cell< Vector< index_t > > mAllNodeIndices ;
43 Cell< Vector< index_t > > mAllOwnedNodeIndices ;
44
45 Cell< Vector< index_t > > mAllEdgeIndices ;
46 Cell< Vector< index_t > > mAllFaceIndices ;
47
48 // vector with polynomial coefficients
49 // we use a matrix so that we can use the trans() function
50 Matrix< real > mPoly ;
51
52 DynamicBitset * mNodeBitset = nullptr ;
53 DynamicBitset * mEdgeBitset = nullptr ;
54 DynamicBitset * mFaceBitset = nullptr ;
55 DynamicBitset * mElementBitset = nullptr ;
56
57 void
58 ( Postprocessor::*mFunComputePoly )( const Vector< real > & aX );
59
60 void
61 ( Postprocessor::*mFunPoly2D )( const real x, const real y );
62
63 void
64 ( Postprocessor::*mFunPoly3D )( const real x, const real y, const real z );
65
66 bool mIsInitialized = false ;
67
68 bool mHaveEdges = false ;
69 bool mHaveFaces = false ;
70
71 ElementType mElementType = ElementType::EMPTY ;
72 id_t mLastElementID = gNoID ;
73 id_t mLastBlockID = gNoID ;
74 Matrix< real > mVandermonde ;
75 Map< mesh::Node * , Matrix< real > * > mNodeMatrices ;
76
77 // the local frame the fit of one node is written in: origin and
78 // in-plane scale for every postprocessor; shell normal, tangents
79 // and thickness scale only in shell mode
80 struct NodeFrame
81 {
82 Vector< real > Origin ;
83 Vector< real > Normal ;
86 real H = 1.0 ;
87 real Tau = 1.0 ;
88 };
89
91
92 uint mNumCoefficients = 0 ;
93
94 // shell mode: linear thin-shell layers ( PENTA6TS, QUAD4TS ) are fit
95 // in ( s1, s2, eta ), eta the curvilinear offset through the stack
96 bool mShellLocal = false ;
97
98 // shell mode, per local node index: 1 if the node lies on the
99 // element's bottom face, 0 if on its top face
100 Cell< uint > mIsBottomNode ;
101
102 // scratch of the element currently linked by link_patch_element()
103 Matrix< real > mNodeCoords ;
104 Matrix< real > mIntPoints ;
105 Vector< real > mZeta ;
106 Vector< real > mShellDirection ;
107 real mElementThickness = 1.0 ;
108
109 // scratch for one local coordinate, and for the per-point results
110 Vector< real > mLocalPoint ;
111 Vector< real > mOffset ;
112 Matrix< real > mPointValues ;
113 Matrix< real > mNodeRhs ;
114
115 // scratch for recovery_block_inverse(), called by invert_node_matrix()
116 Matrix< real > mSchur ;
117 Vector< real > mCoupling ;
118 Vector< real > mSchurCoupling ;
119 Matrix< real > mSvdU ;
120 Matrix< real > mSvdVT ;
121 Vector< real > mSvdS ;
122 Vector< real > mSvdWork ;
123
124 // node fits in which invert_node_matrix() dropped a direction
125 index_t mNumTruncated = 0 ;
126
127 Vector< real > mNull ;
128
129 protected:
130
132
133
138
139 Kernel * mKernel = nullptr;
140 Mesh * mMesh = nullptr;
141 DofManager * mField = nullptr;
142 IWG * mEquation = nullptr;
143 Material * mMaterial = nullptr;
144 Element * mElement = nullptr ;
145 Block * mBlock = nullptr ;
147
148//------------------------------------------------------------------------------
149 public:
150//------------------------------------------------------------------------------
151
156 Postprocessor( Kernel * aKernel , DofManager * aField = nullptr );
157
158 virtual ~Postprocessor() ;
159
160 void
161 set_field( DofManager * aField );
162
163 void
164 set_source_fields( const Cell< string > & aFields );
165
166 void
167 set_target_fields( const Cell< string > & aFields );
168
169 const Cell< string > &
170 source_fields() const ;
171
172 const Cell< string > &
173 target_fields() const ;
174
175 void
176 set_block_ids( const Vector< id_t > & aBlockIDs );
177
178 void
179 set_block_ids( const Cell< id_t > & aBlockIDs );
180
181 virtual void
182 initialize();
183
184 void
186
187 void
189
190 virtual void
191 run();
192
193 DofManager *
194 projector() ;
195
198 index_t
200
201//------------------------------------------------------------------------------
202 protected:
203//------------------------------------------------------------------------------
204
205 const Matrix< real > &
206 compute_poly( const Vector< real > & aX );
207
208 void
210
211 virtual const Vector< real > &
212 compute( const uint aK );
213
214 virtual void
216
217 const Cell< id_t > &
218 block_ids() const ;
219
220 // needed for maxwell postproc
221 // to create element fields
222 const Cell< index_t > &
223 my_element_indices() const ;
224
225//------------------------------------------------------------------------------
226 private:
227//------------------------------------------------------------------------------
228
229 void
230 create_target_fields();
231
232 void
233 check_if_we_have_edges_and_faces();
234
235 void
236 select_elements_and_owned_nodes();
237
238 void
239 select_all_relevant_nodes();
240
241 void
242 synch_node_indices( const bool aOwnedOnly );
243
244 void
245 select_edges();
246
247 void
248 select_faces();
249
250 void
251 select_polynomial();
252
253 // funciton needed by constructor
254 uint
255 get_interpolation_order( const Mesh * aMesh ) const;
256
257 void
258 compute_node_matrices();
259
260 void
261 select_recovery_mode();
262
263 void
264 link_patch_element( mesh::Element * aElement );
265
266 uint
267 local_node_index( const mesh::Element * aElement, const mesh::Node * aNode ) const ;
268
269 real
270 shell_offset( const uint aLocalNode, const uint aK ) const ;
271
272 void
273 compute_node_frame( mesh::Node * aNode, NodeFrame & aFrame );
274
275 void
276 add_node_gram( mesh::Node * aNode, const NodeFrame & aFrame );
277
278 void
279 compute_local_poly( const NodeFrame & aFrame, const uint aLocalNode, const uint aK );
280
281 void
282 invert_node_matrix( Matrix< real > & aInverse );
283
284 void
285 delete_node_matrices();
286
287 void
288 accumulate_node_rhs(
289 mesh::Node * aTarget,
290 const uint aLocalNode,
291 const real aWeight,
292 Map< mesh::Node *, Matrix< real > * > & aNodeCoeffs );
293
294 void
295 synch_source_field( const string & aField );
296
297//------------------------------------------------------------------------------
298// EVALUATION ROUTINES
299//------------------------------------------------------------------------------
300
301 void
302 compute_poly_2d( const Vector< real > & aX );
303
304 void
305 compute_poly_3d( const Vector< real > & aX );
306
307 void
308 poly1_2d( const real x, const real y );
309
310 void
311 poly2_2d( const real x, const real y );
312
313 void
314 poly3_2d( const real x, const real y );
315
316 void
317 poly4_2d( const real x, const real y );
318
319 void
320 poly1_3d( const real x, const real y, const real z );
321
322 void
323 poly2_3d( const real x, const real y, const real z );
324
325 void
326 poly3_3d( const real x, const real y, const real z );
327
328 void
329 poly4_3d( const real x, const real y, const real z );
330 };
331
332 inline
333 DofManager *
335 {
336 return mField;
337 }
338
339 inline index_t
341 {
342 return mNumTruncated ;
343 }
344
345 inline void
347 {
348 mSourceFields = aFields ;
349 mNumSourceFields = aFields.size() ;
350 }
351
352 inline void
354 {
355 mTargetFields = aFields ;
356 mNumTargetFields = aFields.size() ;
357 }
358
359 inline const Cell< string > &
361 {
362 return mSourceFields;
363 }
364
365 inline const Cell< string > &
367 {
368 return mTargetFields;
369 }
370
371 inline const Matrix< real > &
373 {
374 (this->*mFunComputePoly)( aX );
375 return mPoly;
376 }
377
378 inline void
379 Postprocessor::compute_poly_2d( const Vector< real > & aX )
380 {
381 (this->*mFunPoly2D )( aX( 0 ), aX( 1 ) );
382 }
383
384 inline void
385 Postprocessor::compute_poly_3d( const Vector< real > & aX )
386 {
387 (this->*mFunPoly3D )( aX( 0 ), aX( 1 ), aX( 2 ) );
388 }
389
390 inline void
391 Postprocessor::poly1_2d( const real x, const real y )
392 {
393 mPoly( 0, 0 ) = 1.0 ;
394 mPoly( 1, 0 ) = x ;
395 mPoly( 2, 0 ) = y ;
396 }
397
398 inline void
399 Postprocessor::poly2_2d( const real x, const real y )
400 {
401 mPoly( 0, 0 ) = 1.0 ;
402 mPoly( 1, 0 ) = x ;
403 mPoly( 2, 0 ) = y ;
404 mPoly( 3, 0 ) = x*x ;
405 mPoly( 4, 0 ) = x*y ;
406 mPoly( 5, 0 ) = y*y ;
407 }
408
409 inline void
410 Postprocessor::poly3_2d( const real x, const real y )
411 {
412 mPoly( 0, 0 ) = 1.0 ;
413 mPoly( 1, 0 ) = x ;
414 mPoly( 2, 0 ) = y ;
415 mPoly( 3, 0 ) = x*x ;
416 mPoly( 4, 0 ) = x*y ;
417 mPoly( 5, 0 ) = y*y ;
418 mPoly( 6, 0 ) = mPoly( 3, 0 )*x ;
419 mPoly( 7, 0 ) = mPoly( 3, 0 )*y ;
420 mPoly( 8, 0 ) = x * mPoly( 5, 0 );
421 mPoly( 9, 0 ) = y * mPoly( 5, 0 );
422 }
423
424 inline void
425 Postprocessor::poly4_2d( const real x, const real y )
426 {
427 mPoly( 0, 0 ) = 1.0 ;
428 mPoly( 1, 0 ) = x ;
429 mPoly( 2, 0 ) = y ;
430 mPoly( 3, 0 ) = x*x ;
431 mPoly( 4, 0 ) = x*y ;
432 mPoly( 5, 0 ) = y*y ;
433 mPoly( 6, 0 ) = mPoly( 3, 0 )*x ; // x^3
434 mPoly( 7, 0 ) = mPoly( 3, 0 )*y ; // x^2 * y
435 mPoly( 8, 0 ) = x * mPoly( 5, 0 ); // x * y^2
436 mPoly( 9, 0 ) = y * mPoly( 5, 0 ); // y^3
437 mPoly( 10, 0 ) = mPoly( 3, 0 ) * mPoly( 3, 0 ); // x^4
438 mPoly( 11, 0 ) = mPoly( 6, 0 ) * y ; // x^3 * y
439 mPoly( 12, 0 ) = mPoly( 3, 0 ) * mPoly( 5, 0 ) ; // x^2*y^2
440 mPoly( 13, 0 ) = x * mPoly( 9, 0 ) ; // x * y^3
441 mPoly( 14, 0 ) = mPoly( 5, 0 ) * mPoly( 5, 0 ) ; // y^4
442 }
443
444 inline void
445 Postprocessor::poly1_3d( const real x, const real y, const real z )
446 {
447 mPoly( 0, 0 ) = 1.0 ;
448 mPoly( 1, 0 ) = x ;
449 mPoly( 2, 0 ) = y ;
450 mPoly( 3, 0 ) = z ;
451 }
452
453 inline void
454 Postprocessor::poly2_3d( const real x, const real y, const real z )
455 {
456 mPoly( 0, 0 ) = 1.0 ;
457 mPoly( 1, 0 ) = x ;
458 mPoly( 2, 0 ) = y ;
459 mPoly( 3, 0 ) = z ;
460 mPoly( 4, 0 ) = x*x ;
461 mPoly( 5, 0 ) = x*y ;
462 mPoly( 6, 0 ) = y*y ;
463 mPoly( 7, 0 ) = y*z ;
464 mPoly( 8, 0 ) = z*z ;
465 mPoly( 9, 0 ) = z*x ;
466 }
467
468 inline void
469 Postprocessor::poly3_3d( const real x, const real y, const real z )
470 {
471 mPoly( 0, 0 ) = 1.0 ;
472 mPoly( 1, 0 ) = x ;
473 mPoly( 2, 0 ) = y ;
474 mPoly( 3, 0 ) = z ;
475 mPoly( 4, 0 ) = x*x ;
476 mPoly( 5, 0 ) = x*y ;
477 mPoly( 6, 0 ) = y*y ;
478 mPoly( 7, 0 ) = y*z ;
479 mPoly( 8, 0 ) = z*z ;
480 mPoly( 9, 0 ) = z*x ;
481 mPoly( 10, 0 ) = mPoly( 4, 0 ) * x ; // x^3
482 mPoly( 11, 0 ) = mPoly( 4, 0 ) * y ; // x^2 * y
483 mPoly( 12, 0 ) = x * mPoly( 6, 0 ) ; // x * y^2 ;
484 mPoly( 13, 0 ) = y * mPoly( 6, 0 ) ; // y^3
485 mPoly( 14, 0 ) = z * mPoly( 6, 0 ) ; // y^2 * z
486 mPoly( 15, 0 ) = y * mPoly( 8, 0 ) ; // y * z^2
487 mPoly( 16, 0 ) = z * mPoly( 8, 0 ) ; // z^3
488 mPoly( 17, 0 ) = x * mPoly( 8, 0 ) ; // z^2 * x
489 mPoly( 18, 0 ) = mPoly( 4, 0 ) * z ; // x^2 * z
490 mPoly( 19, 0 ) = x * y * z ;
491 }
492
493 inline void
494 Postprocessor::poly4_3d( const real x, const real y, const real z )
495 {
496 mPoly( 0, 0 ) = 1.0 ;
497 mPoly( 1, 0 ) = x ;
498 mPoly( 2, 0 ) = y ;
499 mPoly( 3, 0 ) = z ;
500 mPoly( 4, 0 ) = x*x ;
501 mPoly( 5, 0 ) = x*y ;
502 mPoly( 6, 0 ) = y*y ;
503 mPoly( 7, 0 ) = y*z ;
504 mPoly( 8, 0 ) = z*z ;
505 mPoly( 9, 0 ) = z*x ;
506 mPoly( 10, 0 ) = mPoly( 4, 0 ) * x ; // x^3
507 mPoly( 11, 0 ) = mPoly( 4, 0 ) * y ; // x^2 * y
508 mPoly( 12, 0 ) = x * mPoly( 6, 0 ) ; // x * y^2 ;
509 mPoly( 13, 0 ) = y * mPoly( 6, 0 ) ; // y^3
510 mPoly( 14, 0 ) = z * mPoly( 6, 0 ) ; // y^2 * z
511 mPoly( 15, 0 ) = y * mPoly( 8, 0 ) ; // y * z^2
512 mPoly( 16, 0 ) = z * mPoly( 8, 0 ) ; // z^3
513 mPoly( 17, 0 ) = x * mPoly( 8, 0 ) ; // z^2 * x
514 mPoly( 18, 0 ) = mPoly( 4, 0 ) * z ; // x^2 * z
515 mPoly( 19, 0 ) = x * y * z ;
516 mPoly( 20, 0 ) = mPoly( 4, 0 ) * mPoly( 4, 0 ) ; // x^4
517 mPoly( 21, 0 ) = mPoly( 4, 0 ) * mPoly( 5, 0 ) ; // x^3 * y
518 mPoly( 22, 0 ) = mPoly( 4, 0 ) * mPoly( 6, 0 ) ; // x^2 * y^2
519 mPoly( 23, 0 ) = mPoly( 5, 0 ) * mPoly( 6, 0 ) ; // x * y^3
520 mPoly( 24, 0 ) = mPoly( 6, 0 ) * mPoly( 6, 0 ) ; // y^4
521 mPoly( 25, 0 ) = mPoly( 6, 0 ) * mPoly( 7, 0 ) ; // y^3 * z
522 mPoly( 26, 0 ) = mPoly( 6, 0 ) * mPoly( 8, 0 ) ; // y^2 * z^2
523 mPoly( 27, 0 ) = mPoly( 7, 0 ) * mPoly( 8, 0 ) ; // y * z^3
524 mPoly( 28, 0 ) = mPoly( 8, 0 ) * mPoly( 8, 0 ) ; // z^4
525 mPoly( 29, 0 ) = mPoly( 9, 0 ) * mPoly( 8, 0 ) ; // z^3 * x
526 mPoly( 30, 0 ) = mPoly( 4, 0 ) * mPoly( 8, 0 ) ; // x^2 * z^2
527 mPoly( 31, 0 ) = mPoly( 9, 0 ) * mPoly( 4, 0 ) ; // z * x^3
528 mPoly( 32, 0 ) = mPoly( 4, 0 ) * mPoly( 7, 0 ) ; // x^2 * y * z
529 mPoly( 33, 0 ) = mPoly( 6, 0 ) * mPoly( 9, 0 ) ; // x * y^2 * z
530 mPoly( 34, 0 ) = mPoly( 5, 0 ) * mPoly( 8, 0 ) ; // x * y * z^2
531 }
532
533 inline const Cell< id_t > &
535 {
536 return mBlockIDs ;
537 }
538
539 inline const Cell< index_t > &
541 {
542 return mMyElementIndices ;
543 }
544
545
546 }
547}
548#endif //BELFEM_CL_FEM_POSTPROCESSOR_HPP
size_t size() const
Definition cl_Cell.hpp:189
Runtime-sized bitset; one bit per flag, packed into 64-bit words.
Definition cl_DynamicBitset.hpp:33
Hash map (unordered key-value).
Definition cl_Map.hpp:75
Base class for all materials in BELFEM.
Definition cl_Material.hpp:270
Top-level container for all mesh entities.
Definition cl_Mesh.hpp:60
Definition cl_FEM_Block.hpp:35
Definition cl_FEM_Calculator.hpp:596
Main workhorse: DOF management, assembly coordination and the solver interface.
Definition cl_FEM_DofManager.hpp:52
Definition cl_FEM_Element.hpp:41
Prototype for Integrator of Weak Form Governing Equation.
Definition cl_IWG.hpp:81
Top-level orchestrator.
Definition cl_FEM_Kernel.hpp:56
void set_source_fields(const Cell< string > &aFields)
Definition cl_FEM_Postprocessor.hpp:346
Kernel * mKernel
Definition cl_FEM_Postprocessor.hpp:139
uint mNumSourceFields
Definition cl_FEM_Postprocessor.hpp:136
virtual void initialize()
Definition cl_FEM_Postprocessor.cpp:97
virtual void run()
Definition cl_FEM_Postprocessor.cpp:1395
void recover_fields()
Definition cl_FEM_Postprocessor.cpp:1241
const Cell< index_t > & my_element_indices() const
Definition cl_FEM_Postprocessor.hpp:540
void set_block_ids(const Vector< id_t > &aBlockIDs)
Definition cl_FEM_Postprocessor.cpp:80
IWG * mEquation
Definition cl_FEM_Postprocessor.hpp:142
const Matrix< real > & compute_poly(const Vector< real > &aX)
Definition cl_FEM_Postprocessor.hpp:372
const Cell< string > & target_fields() const
Definition cl_FEM_Postprocessor.hpp:366
const uint mOrder
Definition cl_FEM_Postprocessor.hpp:30
DomainType mDomainType
Definition cl_FEM_Postprocessor.hpp:131
void synch_source_fields()
Definition cl_FEM_Postprocessor.cpp:140
Element * mElement
Definition cl_FEM_Postprocessor.hpp:144
DofManager * projector()
Definition cl_FEM_Postprocessor.hpp:334
void synch_target_fields(Matrix< real > &aData)
Definition cl_FEM_Postprocessor.cpp:149
DofManager * mField
Definition cl_FEM_Postprocessor.hpp:141
const proc_t mCommSize
Definition cl_FEM_Postprocessor.hpp:28
const Cell< string > & source_fields() const
Definition cl_FEM_Postprocessor.hpp:360
void set_target_fields(const Cell< string > &aFields)
Definition cl_FEM_Postprocessor.hpp:353
void set_field(DofManager *aField)
Definition cl_FEM_Postprocessor.cpp:67
uint mNumTargetFields
Definition cl_FEM_Postprocessor.hpp:137
const proc_t mCommRank
Definition cl_FEM_Postprocessor.hpp:27
Block * mBlock
Definition cl_FEM_Postprocessor.hpp:145
Cell< string > mSourceFields
Definition cl_FEM_Postprocessor.hpp:135
Mesh * mMesh
Definition cl_FEM_Postprocessor.hpp:140
const Cell< id_t > & block_ids() const
Definition cl_FEM_Postprocessor.hpp:534
virtual const Vector< real > & compute(const uint aK)
Definition cl_FEM_Postprocessor.cpp:638
index_t number_of_truncated_patches() const
node fits of this rank in which the inversion dropped a direction; set by initialize()
Definition cl_FEM_Postprocessor.hpp:340
Postprocessor(Kernel *aKernel, DofManager *aField=nullptr)
Collective.
Definition cl_FEM_Postprocessor.cpp:27
const uint mNumDimensions
Definition cl_FEM_Postprocessor.hpp:29
virtual void update_element_dofs()
Definition cl_FEM_Postprocessor.cpp:645
Material * mMaterial
Definition cl_FEM_Postprocessor.hpp:143
Cell< string > mTargetFields
Definition cl_FEM_Postprocessor.hpp:134
Calculator * mCalculator
Definition cl_FEM_Postprocessor.hpp:146
Definition cl_IFB_LINE3.hpp:21
USER GUIDES:
Definition cl_Capacitor.cpp:16
constexpr id_t gNoID
Definition typedefs.hpp:58
unsigned int uint
Definition typedefs.hpp:30
unsigned int id_t
Definition typedefs.hpp:41
int proc_t
Definition commtypes.hpp:29
ElementType
Element types.
Definition Mesh_Enums.hpp:27
@ EMPTY
Definition Mesh_Enums.hpp:28
uint32_t index_t
Definition typedefs.hpp:52
double real
Definition typedefs.hpp:36
DomainType
Definition en_DomainType.hpp:20
@ Default
Definition en_DomainType.hpp:21
float x
Definition test_curve_frame.py:37
Definition Element.py:1
Definition Node.py:1
y
Definition test_stadium.py:44