BELFEM 0.9.0
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 real & 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 // parallel constructor
103 Spline( const proc_t aMasterProc );
104
105//------------------------------------------------------------------------------
106
107 // constructor that creates empty container
108 Spline( const index_t & aN, const real aXmin, const real aXmax );
109
110//------------------------------------------------------------------------------
111
112 index_t
113 n() const;
114
115//------------------------------------------------------------------------------
116
117 real
118 x_min() const;
119
120//------------------------------------------------------------------------------
121
122 real
123 x_max() const;
124
125//------------------------------------------------------------------------------
126
127 real
128 delta_x() const;
129
130//------------------------------------------------------------------------------
131
132 Matrix< real > &
133 coefficients();
134
135//------------------------------------------------------------------------------
136
137 const Matrix< real > &
138 coefficients() const;
139
140//------------------------------------------------------------------------------
141
145 real
146 eval( const real aX ) const;
147
148//------------------------------------------------------------------------------
149
153 real
154 deval( const real aX ) const;
155
156//------------------------------------------------------------------------------
157
161 real
162 eval( const real aX, const index_t aCol ) const;
163
164//------------------------------------------------------------------------------
168 real
169 deval( const real aX, const index_t aCol ) const;
170
171//------------------------------------------------------------------------------
172
176 real
177 ddeval( const real aX ) const;
178
179 real
180 ddeval( const real aX, const index_t aCol ) const;
181
182//------------------------------------------------------------------------------
183
187 real
188 entropy( const real aX ) const;
189
190 real
191 entropy( const real aX, const index_t aCol ) const;
192
193//------------------------------------------------------------------------------
194
198 real
199 dentropy( const real aX ) const;
200
201 real
202 dentropy( const real aX, const index_t aCol ) const;
203
204//------------------------------------------------------------------------------
205
209 inline Matrix< real > &
210 matrix_data()
211 {
212 return mData;
213 }
214
215//------------------------------------------------------------------------------
216
226 inline void
227 set_extra_mode( const spline::ExtraMode aMode )
228 {
229 mExtraMode = aMode ;
230 }
231
232//------------------------------------------------------------------------------
233
237 inline const Matrix< real > &
238 matrix_data() const
239 {
240 return mData;
241 }
242
243//------------------------------------------------------------------------------
244
245 void
246 save( const string & aLabel,
247 const string & aPath,
248 const enum FileMode aMode=FileMode::NEW );
249
250//------------------------------------------------------------------------------
251
252 herr_t
253 save( hid_t & aGroup );
254
255//------------------------------------------------------------------------------
256
257 herr_t
258 load( hid_t & aGroup );
259
260//------------------------------------------------------------------------------
261
262 void
263 save_to_database( const string & aDatabase, const string & aLabel );
264
265//------------------------------------------------------------------------------
266
278 void
279 update_data(
280 SpMatrix & aHelpMatrix,
281 const Vector< real > & aValues,
282 const spline::SplineBC aStartBC=spline::SplineBC::NoCurvature,
283 const spline::SplineBC aEndBC=spline::SplineBC::NoCurvature,
284 const real adYdX0=0.0,
285 const real adYdX1=0.0,
286 const real aXref=0.0,
287 const real aSref=0.0 );
288
289//------------------------------------------------------------------------------
290
291 void
292 initialize( const Vector< real > & aX,
293 const Vector< real > & aY,
294 SpMatrix & aA,
295 const spline::SplineBC aStartBC=spline::SplineBC::NoCurvature,
296 const spline::SplineBC aEndBC=spline::SplineBC::NoCurvature,
297 const real adYdX0=0.0,
298 const real adYdX1=0.0,
299 const real aXref=0.0,
300 const real aSref=0.0 );
301
302
303//------------------------------------------------------------------------------
304
305 inline index_t
306 find_col( const real aX ) const;
307
308//------------------------------------------------------------------------------
309
310 void
311 create_integral( const real aXref=BELFEM_QUIET_NAN, const real aYref=0.0 );
312
313//------------------------------------------------------------------------------
314
315 real
316 integrate( const real aX ) const;
317
318 real
319 integrate( const real aX0, const real aX1 ) const ;
320
321//------------------------------------------------------------------------------
322 private:
323//------------------------------------------------------------------------------
324
328 real
329 check_input( const Vector< real > & aX,
330 const Vector< real > & aY,
331 SpMatrix & aA);
332//------------------------------------------------------------------------------
333
337 void
338 create_rhs(
339 const Vector< real > & aY,
340 Vector< real > & aB,
341 const spline::SplineBC aStartBC,
342 const spline::SplineBC aEndBC,
343 const real adYdX0,
344 const real adYdX1);
345
346//------------------------------------------------------------------------------
347
351 void
352 create_coeffs(
353 const Vector< real > & aX,
354 const Vector< real > & aY,
355 const Vector< real > & aDYDX );
356
357//------------------------------------------------------------------------------
358
362 void
363 create_coeffs(
364 const Vector< real > & aY,
365 const Vector< real > & aDYDX );
366
367
368//------------------------------------------------------------------------------
369
373 void
374 create_entropy( const real aXref, const real & aSref );
375
376//------------------------------------------------------------------------------
377
381 void
382 synchronize( const proc_t aMasterProc );
383
384//------------------------------------------------------------------------------
385
386 void
387 add_row_to_data();
388
389//------------------------------------------------------------------------------
390 };
391
392//------------------------------------------------------------------------------
393
394 inline index_t
395 Spline::n() const
396 {
397 return mNumberOfPoints;
398 }
399
400//------------------------------------------------------------------------------
401
402 inline real
403 Spline::x_min() const
404 {
405 return mXmin;
406 }
407
408//------------------------------------------------------------------------------
409
410 inline real
411 Spline::x_max() const
412 {
413 return mXmax;
414 }
415
416//------------------------------------------------------------------------------
417
418 inline real
419 Spline::delta_x() const
420 {
421 return mDeltaX;
422 }
423
424//------------------------------------------------------------------------------
425
426 inline Matrix< real > &
427 Spline::coefficients()
428 {
429 return mData;
430 }
431
432//------------------------------------------------------------------------------
433
434 inline const Matrix< real > &
435 Spline::coefficients() const
436 {
437 return mData;
438 }
439
440//------------------------------------------------------------------------------
441
442 inline index_t
443 Spline::find_col( const real aX ) const
444 {
445 // out-of-range x: clamp the interval index only; eval() then extrapolates
446 // the edge cubic (see spline_usage_guide.md, "Extrapolation behaviour")
447 //BELFEM_ASSERT(mXmin <= aX && aX <= mXmax,
448 // "X=%f out of bounds [%f, %f]", aX, mXmin, mXmax);
449
450 const index_t col = static_cast<index_t>((std::min(std::max(aX,mXmin),mXmax) - mXmin) * mInvDeltaX);
451 return std::min(col, mNumberOfIntervals - 1);
452 }
453
454
455//------------------------------------------------------------------------------
456
457 inline real
458 Spline::eval( const real aX ) const
459 {
460 index_t tCol = find_col( aX );
461
462 return ( ( mData( 0, tCol ) * aX
463 + mData( 1, tCol ) ) * aX
464 + mData( 2, tCol ) ) * aX
465 + mData( 3, tCol );
466 }
467
468//------------------------------------------------------------------------------
469
470 inline real
471 Spline::deval( const real aX ) const
472 {
473 auto tCol = find_col( aX );
474
475 return ( ( 3.0 * mData( 0, tCol ) * aX
476 + 2.0 * mData( 1, tCol ) ) * aX
477 + mData( 2, tCol ) );
478 }
479
480//------------------------------------------------------------------------------
481
482 inline real
483 Spline::eval( const real aX, const index_t aCol ) const
484 {
485 BELFEM_ASSERT( aCol < mNumberOfIntervals,
486 "Spline column %lu out of range, must be < %lu",
487 ( long unsigned int ) aCol, ( long unsigned int ) mNumberOfIntervals );
488
489 return ( ( mData( 0, aCol ) * aX
490 + mData( 1, aCol ) ) * aX
491 + mData( 2, aCol ) ) * aX
492 + mData( 3, aCol );
493 }
494
495//------------------------------------------------------------------------------
496
497 inline real
498 Spline::deval( const real aX, const index_t aCol ) const
499 {
500 BELFEM_ASSERT( aCol < mNumberOfIntervals,
501 "Spline column %lu out of range, must be < %lu",
502 ( long unsigned int ) aCol, ( long unsigned int ) mNumberOfIntervals );
503
504 return ( ( 3.0 * mData( 0, aCol ) * aX
505 + 2.0 * mData( 1, aCol ) ) * aX
506 + mData( 2, aCol ) );
507 }
508
509//------------------------------------------------------------------------------
510
511 inline real
512 Spline::ddeval( const real aX ) const
513 {
514 auto tCol = find_col( aX );
515
516 return 6.0 * mData( 0, tCol ) * aX + 2.0 * mData( 1, tCol );
517 }
518
519 inline real
520 Spline::ddeval( const real aX, const index_t aCol ) const
521 {
522 BELFEM_ASSERT( aCol < mNumberOfIntervals,
523 "Spline column %lu out of range, must be < %lu",
524 ( long unsigned int ) aCol, ( long unsigned int ) mNumberOfIntervals );
525
526 return 6.0 * mData( 0, aCol ) * aX + 2.0 * mData( 1, aCol );
527 }
528
529//------------------------------------------------------------------------------
530
531 inline real
532 Spline::entropy( const real aX ) const
533 {
534 BELFEM_ASSERT( mExtraMode == spline::ExtraMode::Entropy,
535 "No entropy tables present" );
536 auto tCol = find_col( aX );
537 return ( 1.5 * mData( 0, tCol ) * aX
538 + 2.0 * mData( 1, tCol ) ) * aX
539 + mData( 2, tCol ) * std::log( aX )
540 + mData( 4, tCol );
541 }
542
543 inline real
544 Spline::entropy( const real aX, const index_t aCol ) const
545 {
546 BELFEM_ASSERT( mExtraMode == spline::ExtraMode::Entropy,
547 "No entropy tables present" );
548
549 BELFEM_ASSERT( aCol < mNumberOfIntervals,
550 "Spline column %lu out of range, must be < %lu",
551 ( long unsigned int ) aCol, ( long unsigned int ) mNumberOfIntervals );
552
553
554 return ( 1.5 * mData( 0, aCol ) * aX
555 + 2.0 * mData( 1, aCol ) ) * aX
556 + mData( 2, aCol ) * std::log( aX )
557 + mData( 4, aCol );
558 }
559
560
561//------------------------------------------------------------------------------
562
563 inline real
564 Spline::dentropy( const real aX ) const
565 {
566 BELFEM_ASSERT( mExtraMode == spline::ExtraMode::Entropy,
567 "No entropy tables present" );
568
569 auto tCol = find_col( aX );
570 return 3.0 * mData( 0, tCol ) * aX
571 + 2.0 * mData( 1, tCol )
572 + mData( 2, tCol )/aX;
573 }
574
575 inline real
576 Spline::dentropy( const real aX, const index_t aCol ) const
577 {
578 BELFEM_ASSERT( mExtraMode == spline::ExtraMode::Entropy,
579 "No entropy tables present" );
580
581 BELFEM_ASSERT( aCol < mNumberOfIntervals,
582 "Spline column %lu out of range, must be < %lu",
583 ( long unsigned int ) aCol, ( long unsigned int ) mNumberOfIntervals );
584
585
586 return 3.0 * mData( 0, aCol ) * aX
587 + 2.0 * mData( 1, aCol )
588 + mData( 2, aCol )/aX;
589 }
590
591 inline real
592 Spline::integrate( const real aX ) const
593 {
594 BELFEM_ASSERT( mExtraMode == spline::ExtraMode::Integral,
595 "No integration tables present" );
596
597 index_t tCol = find_col( aX );
598
599 return ((( 0.25 * mData( 0, tCol ) * aX
600 + mData( 1, tCol )/3.0 ) * aX
601 + 0.5 * mData( 2, tCol ) ) * aX
602 + mData( 3, tCol ) ) * aX
603 + mData( 4, tCol );
604 }
605
606 inline real
607 Spline::integrate( const real aX0, const real aX1 ) const
608 {
609 BELFEM_ASSERT( mExtraMode == spline::ExtraMode::Integral,
610 "No integration tables present" );
611
612 return this->integrate( aX1 ) - this->integrate( aX0 );
613 }
614
615//------------------------------------------------------------------------------
616}
617
618#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:801
belfem::Spline::save
void save(const string &aLabel, const string &aPath, const enum FileMode aMode=FileMode::NEW)
Definition cl_Spline.cpp:333
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:409
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:201
belfem::Spline::ddeval
real ddeval(const real aX) const
interpolate second derivative
Definition cl_Spline.hpp:512
belfem::Spline::dentropy
real dentropy(const real aX) const
interpolate entropy derivative
Definition cl_Spline.hpp:564
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:227
belfem::Spline::integrate
real integrate(const real aX) const
Definition cl_Spline.hpp:592
belfem::Spline::load
herr_t load(hid_t &aGroup)
Definition cl_Spline.cpp:376
belfem::Spline::find_col
index_t find_col(const real aX) const
Definition cl_Spline.hpp:443
belfem::Spline::delta_x
real delta_x() const
Definition cl_Spline.hpp:419
belfem::Spline::x_min
real x_min() const
Definition cl_Spline.hpp:403
belfem::Spline::n
index_t n() const
Definition cl_Spline.hpp:395
belfem::Spline::matrix_data
const Matrix< real > & matrix_data() const
expose matrix
Definition cl_Spline.hpp:238
belfem::Spline::create_integral
void create_integral(const real aXref=BELFEM_QUIET_NAN, const real aYref=0.0)
Definition cl_Spline.cpp:692
belfem::Spline::deval
real deval(const real aX) const
interpolate first derivative
Definition cl_Spline.hpp:471
belfem::Spline::coefficients
Matrix< real > & coefficients()
Definition cl_Spline.hpp:427
belfem::Spline::entropy
real entropy(const real aX) const
interpolate entropy
Definition cl_Spline.hpp:532
belfem::Spline::eval
real eval(const real aX) const
interpolate the function
Definition cl_Spline.hpp:458
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:271
belfem::Spline::x_max
real x_max() const
Definition cl_Spline.hpp:411
belfem::Spline::matrix_data
Matrix< real > & matrix_data()
expose matrix
Definition cl_Spline.hpp:210
hdf5_tools.hpp
belfem::spline
Definition cl_Spline.cpp:33
belfem::spline::create_helpmatrix
void create_helpmatrix(const real &aSize, const real &aDeltaX, SpMatrix &aA, const SplineBC aStartBC, const SplineBC aEndBC)
Definition cl_Spline.cpp:39
belfem::spline::SplineBC
SplineBC
Definition Spline_Enums.hpp:19
belfem::spline::SplineBC::NoCurvature
@ NoCurvature
Definition Spline_Enums.hpp:20
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:87