BELFEM 0.9.3
Berkeley Lab Finite Element Framework
Loading...
Searching...
No Matches
cl_Spline.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 BELFEM_SPLINE_HPP
13#define BELFEM_SPLINE_HPP
14
15#include "typedefs.hpp"
16#include "cl_Vector.hpp"
17#include "cl_SpMatrix.hpp"
18#include "filetools.hpp"
19
20#include "hdf5_tools.hpp"
21#include "Spline_Enums.hpp"
22
23namespace belfem
24{
25 namespace spline
26 {
27 void
28 create_helpmatrix(
29 const index_t aSize,
30 const real aDeltaX,
31 SpMatrix & aA,
32 const SplineBC aStartBC = SplineBC::NoCurvature,
33 const SplineBC aEndBC = SplineBC::NoCurvature );
34 }
35
36//------------------------------------------------------------------------------
37
44 class Spline
45 {
46 const proc_t mCommRank ;
47
48 index_t mNumberOfPoints;
49 index_t mNumberOfIntervals;
50
51 real mXmin;
52 real mXmax;
53 real mDeltaX;
54 real mInvDeltaX;
55
56 spline::ExtraMode mExtraMode = spline::ExtraMode::None;
57
58 Matrix< real > mData;
59
60//------------------------------------------------------------------------------
61 public:
62//------------------------------------------------------------------------------
63
64 // constructor with helpmatrix, without boundary conditions
65 //
66 // aXref carries no default on purpose. With one, a three argument call
67 // matches this and the overload below equally well and does not compile.
68 // The two agree there anyway: this one forwards NoCurvature at both ends
69 // and zero slopes, which is exactly what the other one defaults to.
70 Spline( const Vector< real > & aX,
71 const Vector< real > & aY,
72 SpMatrix & aA,
73 const real aXref,
74 const real aSref = 0.,
75 const proc_t aMasterProc = gNoOwner );
76
77 // constructor with helpmatrix
78 Spline( const Vector< real > & aX,
79 const Vector< real > & aY,
80 SpMatrix & aA,
81 const spline::SplineBC aStartBC = spline::SplineBC::NoCurvature,
82 const spline::SplineBC aEndBC = spline::SplineBC::NoCurvature,
83 const real adYdX0 = 0.,
84 const real adYdX1 = 0.,
85 const real aXref = 0.,
86 const real aSref = 0.,
87 const proc_t aMasterProc = gNoOwner );
88
89//------------------------------------------------------------------------------
90
91 Spline(
92 const string & aFile,
93 const string & aLabel,
94 const proc_t aMaster=0 );
95
96//------------------------------------------------------------------------------
97
98 Spline( const hid_t aGroup, const proc_t aMaster );
99
100//------------------------------------------------------------------------------
101
102 // Collective: the master builds the coefficients with a data
103 // constructor and calls synchronize( master ). Every other rank enters
104 // this constructor with the same master rank and receives them through
105 // an internal barrier and broadcasts.
106 Spline( const proc_t aMasterProc );
107
108//------------------------------------------------------------------------------
109
110 // constructor that creates empty container
111 Spline( const index_t & aN, const real aXmin, const real aXmax );
112
113//------------------------------------------------------------------------------
114
115 index_t
116 n() const;
117
118//------------------------------------------------------------------------------
119
120 real
121 x_min() const;
122
123//------------------------------------------------------------------------------
124
125 real
126 x_max() const;
127
128//------------------------------------------------------------------------------
129
130 real
131 delta_x() const;
132
133//------------------------------------------------------------------------------
134
135 Matrix< real > &
136 coefficients();
137
138//------------------------------------------------------------------------------
139
140 const Matrix< real > &
141 coefficients() const;
142
143//------------------------------------------------------------------------------
144
148 real
149 eval( const real aX ) const;
150
151//------------------------------------------------------------------------------
152
156 real
157 deval( const real aX ) const;
158
159//------------------------------------------------------------------------------
160
164 real
165 eval( const real aX, const index_t aCol ) const;
166
167//------------------------------------------------------------------------------
171 real
172 deval( const real aX, const index_t aCol ) const;
173
174//------------------------------------------------------------------------------
175
179 real
180 ddeval( const real aX ) const;
181
182 real
183 ddeval( const real aX, const index_t aCol ) const;
184
185//------------------------------------------------------------------------------
186
190 real
191 entropy( const real aX ) const;
192
193 real
194 entropy( const real aX, const index_t aCol ) const;
195
196//------------------------------------------------------------------------------
197
201 real
202 dentropy( const real aX ) const;
203
204 real
205 dentropy( const real aX, const index_t aCol ) const;
206
207//------------------------------------------------------------------------------
208
212 inline Matrix< real > &
213 matrix_data()
214 {
215 return mData;
216 }
217
218//------------------------------------------------------------------------------
219
229 inline void
230 set_extra_mode( const spline::ExtraMode aMode )
231 {
232 mExtraMode = aMode ;
233 }
234
235//------------------------------------------------------------------------------
236
240 inline const Matrix< real > &
241 matrix_data() const
242 {
243 return mData;
244 }
245
246//------------------------------------------------------------------------------
247
248 void
249 save( const string & aLabel,
250 const string & aPath,
251 const enum FileMode aMode=FileMode::NEW );
252
253//------------------------------------------------------------------------------
254
255 herr_t
256 save( hid_t & aGroup );
257
258//------------------------------------------------------------------------------
259
260 herr_t
261 load( hid_t & aGroup );
262
263//------------------------------------------------------------------------------
264
265 void
266 save_to_database( const string & aDatabase, const string & aLabel );
267
268//------------------------------------------------------------------------------
269
281 void
282 update_data(
283 SpMatrix & aHelpMatrix,
284 const Vector< real > & aValues,
285 const spline::SplineBC aStartBC=spline::SplineBC::NoCurvature,
286 const spline::SplineBC aEndBC=spline::SplineBC::NoCurvature,
287 const real adYdX0=0.0,
288 const real adYdX1=0.0,
289 const real aXref=0.0,
290 const real aSref=0.0 );
291
292//------------------------------------------------------------------------------
293
294 void
295 initialize( const Vector< real > & aX,
296 const Vector< real > & aY,
297 SpMatrix & aA,
298 const spline::SplineBC aStartBC=spline::SplineBC::NoCurvature,
299 const spline::SplineBC aEndBC=spline::SplineBC::NoCurvature,
300 const real adYdX0=0.0,
301 const real adYdX1=0.0,
302 const real aXref=0.0,
303 const real aSref=0.0 );
304
305
306//------------------------------------------------------------------------------
307
308 inline index_t
309 find_col( const real aX ) const;
310
311//------------------------------------------------------------------------------
312
313 void
314 create_integral( const real aXref=BELFEM_QUIET_NAN, const real aYref=0.0 );
315
316//------------------------------------------------------------------------------
317
318 real
319 integrate( const real aX ) const;
320
321 real
322 integrate( const real aX0, const real aX1 ) const ;
323
324//------------------------------------------------------------------------------
325 private:
326//------------------------------------------------------------------------------
327
331 real
332 check_input( const Vector< real > & aX,
333 const Vector< real > & aY,
334 SpMatrix & aA);
335//------------------------------------------------------------------------------
336
340 void
341 create_rhs(
342 const Vector< real > & aY,
343 Vector< real > & aB,
344 const spline::SplineBC aStartBC,
345 const spline::SplineBC aEndBC,
346 const real adYdX0,
347 const real adYdX1);
348
349//------------------------------------------------------------------------------
350
354 void
355 create_coeffs(
356 const Vector< real > & aX,
357 const Vector< real > & aY,
358 const Vector< real > & aDYDX );
359
360//------------------------------------------------------------------------------
361
365 void
366 create_coeffs(
367 const Vector< real > & aY,
368 const Vector< real > & aDYDX );
369
370
371//------------------------------------------------------------------------------
372
376 void
377 create_entropy( const real aXref, const real & aSref );
378
379//------------------------------------------------------------------------------
380
384 void
385 synchronize( const proc_t aMasterProc );
386
387//------------------------------------------------------------------------------
388
389 void
390 add_row_to_data();
391
392//------------------------------------------------------------------------------
393 };
394
395//------------------------------------------------------------------------------
396
397 inline index_t
398 Spline::n() const
399 {
400 return mNumberOfPoints;
401 }
402
403//------------------------------------------------------------------------------
404
405 inline real
406 Spline::x_min() const
407 {
408 return mXmin;
409 }
410
411//------------------------------------------------------------------------------
412
413 inline real
414 Spline::x_max() const
415 {
416 return mXmax;
417 }
418
419//------------------------------------------------------------------------------
420
421 inline real
422 Spline::delta_x() const
423 {
424 return mDeltaX;
425 }
426
427//------------------------------------------------------------------------------
428
429 inline Matrix< real > &
430 Spline::coefficients()
431 {
432 return mData;
433 }
434
435//------------------------------------------------------------------------------
436
437 inline const Matrix< real > &
438 Spline::coefficients() const
439 {
440 return mData;
441 }
442
443//------------------------------------------------------------------------------
444
445 inline index_t
446 Spline::find_col( const real aX ) const
447 {
448 // out-of-range x: clamp the interval index only; eval() then extrapolates
449 // the edge cubic (see spline_usage_guide.md, "Extrapolation behaviour")
450 //BELFEM_ASSERT(mXmin <= aX && aX <= mXmax,
451 // "X=%f out of bounds [%f, %f]", aX, mXmin, mXmax);
452
453 const index_t col = static_cast<index_t>((std::min(std::max(aX,mXmin),mXmax) - mXmin) * mInvDeltaX);
454 return std::min(col, mNumberOfIntervals - 1);
455 }
456
457
458//------------------------------------------------------------------------------
459
460 inline real
461 Spline::eval( const real aX ) const
462 {
463 index_t tCol = find_col( aX );
464
465 return ( ( mData( 0, tCol ) * aX
466 + mData( 1, tCol ) ) * aX
467 + mData( 2, tCol ) ) * aX
468 + mData( 3, tCol );
469 }
470
471//------------------------------------------------------------------------------
472
473 inline real
474 Spline::deval( const real aX ) const
475 {
476 auto tCol = find_col( aX );
477
478 return ( ( 3.0 * mData( 0, tCol ) * aX
479 + 2.0 * mData( 1, tCol ) ) * aX
480 + mData( 2, tCol ) );
481 }
482
483//------------------------------------------------------------------------------
484
485 inline real
486 Spline::eval( const real aX, const index_t aCol ) const
487 {
488 BELFEM_ASSERT( aCol < mNumberOfIntervals,
489 "Spline column %lu out of range, must be < %lu",
490 ( long unsigned int ) aCol, ( long unsigned int ) mNumberOfIntervals );
491
492 return ( ( mData( 0, aCol ) * aX
493 + mData( 1, aCol ) ) * aX
494 + mData( 2, aCol ) ) * aX
495 + mData( 3, aCol );
496 }
497
498//------------------------------------------------------------------------------
499
500 inline real
501 Spline::deval( const real aX, const index_t aCol ) const
502 {
503 BELFEM_ASSERT( aCol < mNumberOfIntervals,
504 "Spline column %lu out of range, must be < %lu",
505 ( long unsigned int ) aCol, ( long unsigned int ) mNumberOfIntervals );
506
507 return ( ( 3.0 * mData( 0, aCol ) * aX
508 + 2.0 * mData( 1, aCol ) ) * aX
509 + mData( 2, aCol ) );
510 }
511
512//------------------------------------------------------------------------------
513
514 inline real
515 Spline::ddeval( const real aX ) const
516 {
517 auto tCol = find_col( aX );
518
519 return 6.0 * mData( 0, tCol ) * aX + 2.0 * mData( 1, tCol );
520 }
521
522 inline real
523 Spline::ddeval( const real aX, const index_t aCol ) const
524 {
525 BELFEM_ASSERT( aCol < mNumberOfIntervals,
526 "Spline column %lu out of range, must be < %lu",
527 ( long unsigned int ) aCol, ( long unsigned int ) mNumberOfIntervals );
528
529 return 6.0 * mData( 0, aCol ) * aX + 2.0 * mData( 1, aCol );
530 }
531
532//------------------------------------------------------------------------------
533
534 inline real
535 Spline::entropy( const real aX ) const
536 {
537 BELFEM_ASSERT( mExtraMode == spline::ExtraMode::Entropy,
538 "No entropy tables present" );
539 auto tCol = find_col( aX );
540 return ( 1.5 * mData( 0, tCol ) * aX
541 + 2.0 * mData( 1, tCol ) ) * aX
542 + mData( 2, tCol ) * std::log( aX )
543 + mData( 4, tCol );
544 }
545
546 inline real
547 Spline::entropy( const real aX, const index_t aCol ) const
548 {
549 BELFEM_ASSERT( mExtraMode == spline::ExtraMode::Entropy,
550 "No entropy tables present" );
551
552 BELFEM_ASSERT( aCol < mNumberOfIntervals,
553 "Spline column %lu out of range, must be < %lu",
554 ( long unsigned int ) aCol, ( long unsigned int ) mNumberOfIntervals );
555
556
557 return ( 1.5 * mData( 0, aCol ) * aX
558 + 2.0 * mData( 1, aCol ) ) * aX
559 + mData( 2, aCol ) * std::log( aX )
560 + mData( 4, aCol );
561 }
562
563
564//------------------------------------------------------------------------------
565
566 inline real
567 Spline::dentropy( const real aX ) const
568 {
569 BELFEM_ASSERT( mExtraMode == spline::ExtraMode::Entropy,
570 "No entropy tables present" );
571
572 auto tCol = find_col( aX );
573 return 3.0 * mData( 0, tCol ) * aX
574 + 2.0 * mData( 1, tCol )
575 + mData( 2, tCol )/aX;
576 }
577
578 inline real
579 Spline::dentropy( const real aX, const index_t aCol ) const
580 {
581 BELFEM_ASSERT( mExtraMode == spline::ExtraMode::Entropy,
582 "No entropy tables present" );
583
584 BELFEM_ASSERT( aCol < mNumberOfIntervals,
585 "Spline column %lu out of range, must be < %lu",
586 ( long unsigned int ) aCol, ( long unsigned int ) mNumberOfIntervals );
587
588
589 return 3.0 * mData( 0, aCol ) * aX
590 + 2.0 * mData( 1, aCol )
591 + mData( 2, aCol )/aX;
592 }
593
594 inline real
595 Spline::integrate( const real aX ) const
596 {
597 BELFEM_ASSERT( mExtraMode == spline::ExtraMode::Integral,
598 "No integration tables present" );
599
600 index_t tCol = find_col( aX );
601
602 return ((( 0.25 * mData( 0, tCol ) * aX
603 + mData( 1, tCol )/3.0 ) * aX
604 + 0.5 * mData( 2, tCol ) ) * aX
605 + mData( 3, tCol ) ) * aX
606 + mData( 4, tCol );
607 }
608
609 inline real
610 Spline::integrate( const real aX0, const real aX1 ) const
611 {
612 BELFEM_ASSERT( mExtraMode == spline::ExtraMode::Integral,
613 "No integration tables present" );
614
615 return this->integrate( aX1 ) - this->integrate( aX0 );
616 }
617
618//------------------------------------------------------------------------------
619}
620
621#endif // BELFEM_SPLINE_HPP
Spline_Enums.hpp
BELFEM_ASSERT
#define BELFEM_ASSERT(aCheck,...)
Definition assert.hpp:244
cl_SpMatrix.hpp
cl_Vector.hpp
belfem::SpMatrix
Sparse matrix in CSR or CSC format.
Definition cl_SpMatrix.hpp:52
belfem::Spline::save_to_database
void save_to_database(const string &aDatabase, const string &aLabel)
Definition cl_Spline.cpp:785
belfem::Spline::save
void save(const string &aLabel, const string &aPath, const enum FileMode aMode=FileMode::NEW)
Definition cl_Spline.cpp:329
belfem::Spline::update_data
void update_data(SpMatrix &aHelpMatrix, const Vector< real > &aValues, const spline::SplineBC aStartBC=spline::SplineBC::NoCurvature, const spline::SplineBC aEndBC=spline::SplineBC::NoCurvature, const real adYdX0=0.0, const real adYdX1=0.0, const real aXref=0.0, const real aSref=0.0)
Recompute spline coefficients for new y-values on the same grid.
Definition cl_Spline.cpp:404
belfem::Spline::Spline
Spline(const Vector< real > &aX, const Vector< real > &aY, SpMatrix &aA, const real aXref, const real aSref=0., const proc_t aMasterProc=gNoOwner)
Definition cl_Spline.cpp:199
belfem::Spline::ddeval
real ddeval(const real aX) const
interpolate second derivative
Definition cl_Spline.hpp:515
belfem::Spline::dentropy
real dentropy(const real aX) const
interpolate entropy derivative
Definition cl_Spline.hpp:567
belfem::Spline::set_extra_mode
void set_extra_mode(const spline::ExtraMode aMode)
Declare what the extra row of the coefficient table holds.
Definition cl_Spline.hpp:230
belfem::Spline::integrate
real integrate(const real aX) const
Definition cl_Spline.hpp:595
belfem::Spline::load
herr_t load(hid_t &aGroup)
Definition cl_Spline.cpp:371
belfem::Spline::find_col
index_t find_col(const real aX) const
Definition cl_Spline.hpp:446
belfem::Spline::delta_x
real delta_x() const
Definition cl_Spline.hpp:422
belfem::Spline::x_min
real x_min() const
Definition cl_Spline.hpp:406
belfem::Spline::n
index_t n() const
Definition cl_Spline.hpp:398
belfem::Spline::matrix_data
const Matrix< real > & matrix_data() const
expose matrix
Definition cl_Spline.hpp:241
belfem::Spline::create_integral
void create_integral(const real aXref=BELFEM_QUIET_NAN, const real aYref=0.0)
Definition cl_Spline.cpp:677
belfem::Spline::deval
real deval(const real aX) const
interpolate first derivative
Definition cl_Spline.hpp:474
belfem::Spline::coefficients
Matrix< real > & coefficients()
Definition cl_Spline.hpp:430
belfem::Spline::entropy
real entropy(const real aX) const
interpolate entropy
Definition cl_Spline.hpp:535
belfem::Spline::eval
real eval(const real aX) const
interpolate the function
Definition cl_Spline.hpp:461
belfem::Spline::initialize
void initialize(const Vector< real > &aX, const Vector< real > &aY, SpMatrix &aA, const spline::SplineBC aStartBC=spline::SplineBC::NoCurvature, const spline::SplineBC aEndBC=spline::SplineBC::NoCurvature, const real adYdX0=0.0, const real adYdX1=0.0, const real aXref=0.0, const real aSref=0.0)
Definition cl_Spline.cpp:269
belfem::Spline::x_max
real x_max() const
Definition cl_Spline.hpp:414
belfem::Spline::matrix_data
Matrix< real > & matrix_data()
expose matrix
Definition cl_Spline.hpp:213
hdf5_tools.hpp
belfem::spline
Definition cl_Spline.cpp:33
belfem::spline::SplineBC
SplineBC
Definition Spline_Enums.hpp:19
belfem::spline::SplineBC::NoCurvature
@ NoCurvature
Definition Spline_Enums.hpp:20
belfem::spline::create_helpmatrix
void create_helpmatrix(const index_t aSize, const real aDeltaX, SpMatrix &aA, const SplineBC aStartBC, const SplineBC aEndBC)
Definition cl_Spline.cpp:39
belfem::spline::ExtraMode
ExtraMode
Definition Spline_Enums.hpp:26
belfem::spline::ExtraMode::Entropy
@ Entropy
Definition Spline_Enums.hpp:28
belfem::spline::ExtraMode::Integral
@ Integral
Definition Spline_Enums.hpp:29
belfem::spline::ExtraMode::None
@ None
Definition Spline_Enums.hpp:27
belfem
USER GUIDES:
Definition cl_Capacitor.cpp:16
belfem::hid_t
int hid_t
Definition hdf5_types.hpp:20
belfem::proc_t
int proc_t
Definition commtypes.hpp:29
belfem::FileMode
FileMode
Definition filetools.hpp:27
belfem::FileMode::NEW
@ NEW
Definition filetools.hpp:28
belfem::gNoOwner
constexpr proc_t gNoOwner
Definition typedefs.hpp:59
belfem::index_t
uint32_t index_t
Definition typedefs.hpp:52
belfem::herr_t
int herr_t
Definition hdf5_types.hpp:21
belfem::real
double real
Definition typedefs.hpp:36
BELFEM_QUIET_NAN
#define BELFEM_QUIET_NAN
Definition typedefs.hpp:86