Sandals  v0.0.0
A C++ library for ODEs/DAEs integration
Loading...
Searching...
No Matches
SemiExplicit.hh
Go to the documentation of this file.
1/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\
2 * Copyright (c) 2025, Davide Stocco and Enrico Bertolazzi. *
3 * *
4 * The Sandals project is distributed under the BSD 2-Clause License. *
5 * *
6 * Davide Stocco Enrico Bertolazzi *
7 * University of Trento University of Trento *
8 * e-mail: davide.stocco@unitn.it e-mail: enrico.bertolazzi@unitn.it *
9\* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
10
11#pragma once
12
13#ifndef SANDALS_SEMIEXPLICIT_SYSTEM_HH
14#define SANDALS_SEMIEXPLICIT_SYSTEM_HH
15
16#include <Sandals.hh>
18
19namespace Sandals {
20
21 /*\
22 | ____ _ _____ _ _ _ _
23 | / ___| ___ _ __ ___ (_) ____|_ ___ __ | (_) ___(_) |_
24 | \___ \ / _ \ '_ ` _ \| | _| \ \/ / '_ \| | |/ __| | __|
25 | ___) | __/ | | | | | | |___ > <| |_) | | | (__| | |_
26 | |____/ \___|_| |_| |_|_|_____/_/\_\ .__/|_|_|\___|_|\__|
27 | |_|
28 \*/
29
41 template <typename Real, Integer N, Integer M = 0>
42 class SemiExplicit : public Explicit<Real, N, M>
43 {
44 public:
45 using Pointer = std::shared_ptr<SemiExplicit<Real, N, M>>;
49 using TensorTA = typename std::vector<MatrixJF>;
53
54 private:
55 mutable Eigen::FullPivLU<MatrixA> m_lu;
56
57 public:
61 SemiExplicit() : Explicit<Real, N, M>(Type::SEMIEXPLICIT, "(missing name)") {}
62
67 SemiExplicit(std::string t_name) : Explicit<Real, N, M>(Type::SEMIEXPLICIT, t_name) {}
68
82 VectorF F(VectorF const & x, VectorF const & x_dot, Real const t) const override
83 {
84 return this->A(x, t)*x_dot - this->b(x, t);
85 }
86
104 MatrixJF JF_x(VectorF const & x, VectorF const & x_dot, Real const t) const override
105 {
106 return -this->Jf_x(x, x_dot, t);
107 }
108
124 MatrixJF JF_x_dot(VectorF const & x, VectorF const &/*x_dot*/, Real const t) const override
125 {
126 return this->A(x, t);
127 }
128
140 VectorF f(VectorF const & x, Real const t) const override
141 {
142 this->m_lu.compute(this->A(x, t));
143 SANDALS_ASSERT(this->m_lu.rank() == N,
144 "Sandals:SemiExplicit::f(...): singular mass matrix A detected.");
145 return this->m_lu.solve(this->b(x, t));
146 }
147
163 MatrixJF Jf_x(VectorF const & x, VectorF const & x_dot, Real const t) const
164 {
165 TensorTA TA_x{this->TA_x(x, t)};
166 MatrixJF tAp{MatrixJF::Zero()};
167 for (Integer i{0}; i < N; ++i) {tAp.col(i) = TA_x[i] * x_dot;} // TODO: Check if TA_x.size() = N
168 this->m_lu.compute(this->A(x, t));
169 SANDALS_ASSERT(this->m_lu.rank() == N,
170 "Sandals:SemiExplicit::Jf_x(...): singular mass matrix A detected.");
171 return this->m_lu.solve(this->Jb_x(x, t) - tAp);
172 }
173
188 MatrixJF Jf_x(VectorF const & x, Real const t) const override {return this->Jf_x(x, this->f(x, t), t);}
189
196 virtual MatrixA A(VectorF const & x, Real const t) const = 0;
197
211 virtual TensorTA TA_x(VectorF const & x, Real const t) const = 0;
212
219 virtual VectorB b(VectorF const & x, Real const t) const = 0;
220
234 virtual MatrixJB Jb_x(VectorF const & x, Real const t) const = 0;
235
236 }; // class SemiExplicit
237
238 /*\
239 | ____ _ _____ _ _ _ _ __ __
240 | / ___| ___ _ __ ___ (_) ____|_ ___ __ | (_) ___(_) |\ \ / / __ __ _ _ __ _ __ ___ _ __
241 | \___ \ / _ \ '_ ` _ \| | _| \ \/ / '_ \| | |/ __| | __\ \ /\ / / '__/ _` | '_ \| '_ \ / _ \ '__|
242 | ___) | __/ | | | | | | |___ > <| |_) | | | (__| | |_ \ V V /| | | (_| | |_) | |_) | __/ |
243 | |____/ \___|_| |_| |_|_|_____/_/\_\ .__/|_|_|\___|_|\__| \_/\_/ |_| \__,_| .__/| .__/ \___|_|
244 | |_| |_| |_|
245 \*/
246
258 template <typename Real, Integer N, Integer M = 0>
259 class SemiExplicitWrapper : public SemiExplicit<Real, N, M>
260 {
261 public:
262 using Pointer = std::shared_ptr<SemiExplicitWrapper<Real, N, M>>;
263 using typename SemiExplicit<Real, N, M>::VectorF;
264 using typename SemiExplicit<Real, N, M>::MatrixA;
266 using typename SemiExplicit<Real, N, M>::VectorB;
268 using typename SemiExplicit<Real, N, M>::VectorH;
270 using FunctionA = std::function<MatrixA(VectorF const &, Real const)>;
271 using FunctionTA = std::function<TensorTA(VectorF const &, Real const)>;
272 using FunctionB = std::function<VectorB(VectorF const &, Real const)>;
273 using FunctionJB = std::function<MatrixJB(VectorF const &, Real const)>;
274 using FunctionH = std::function<VectorH(VectorF const &, Real const)>;
275 using FunctionJH = std::function<MatrixJH(VectorF const &, Real const)>;
276 using FunctionID = std::function<bool(VectorF const &, Real const)>;
277
278 inline const static FunctionH DefaultH = [](VectorF const &, Real const) {return VectorH::Zero();};
279 inline const static FunctionJH DefaultJH = [](VectorF const &, Real const) {return MatrixJH::Zero();};
280 inline const static FunctionID DefaultID = [](VectorF const &, Real const) {return true;};
281
282 private:
283 FunctionA m_A{nullptr};
285 FunctionB m_b{nullptr};
287 FunctionH m_h{nullptr};
290
291 public:
303 FunctionH t_h = DefaultH, FunctionJH t_Jh_x = DefaultJH, FunctionID t_in_domain = DefaultID)
304 : SemiExplicit<Real, N, M>(), m_A(t_A), m_TA_x(t_TA_x), m_b(t_b), m_Jb_x(t_Jb_x), m_h(t_h),
305 m_Jh_x(t_Jh_x), m_in_domain(t_in_domain)
306 {}
307
319 SemiExplicitWrapper(std::string t_name, FunctionA t_A, FunctionTA t_TA_x, FunctionB t_b, FunctionJB t_Jb_x,
320 FunctionH t_h = DefaultH, FunctionJH t_Jh_x = DefaultJH, FunctionID t_in_domain = DefaultID)
321 : SemiExplicit<Real, N, M>(t_name), m_A(t_A), m_TA_x(t_TA_x), m_b(t_b), m_Jb_x(t_Jb_x), m_h(t_h),
322 m_Jh_x(t_Jh_x), m_in_domain(t_in_domain)
323 {}
324
329
334 FunctionA & A() {return this->m_A;}
335
340 FunctionTA & TA_x() {return this->m_TA_x;}
341
346 FunctionB & b() {return this->m_b;}
347
352 FunctionJB & Jb_x() {return this->m_Jb_x;}
353
358 FunctionH & h() {return this->m_h;}
359
364 FunctionJH & Jh_x() {return this->m_Jh_x;}
365
370 FunctionID & in_domain() {return this->m_in_domain;}
371
378 MatrixA A(VectorF const & x, Real const t) const override
379 {
380 return this->m_A(x, t);
381 }
382
396 TensorTA TA_x(VectorF const & x, Real const t) const override
397 {
398 return this->m_TA_x(x, t);
399 }
400
407 VectorB b(VectorF const & x, Real const t) const override
408 {
409 return this->m_b(x, t);
410 }
411
425 MatrixJB Jb_x(VectorF const & x, Real const t) const override
426 {
427 return this->m_Jb_x(x, t);
428 }
429
430
437 VectorH h(VectorF const & x, Real const t) const override
438 {
439 return this->m_h(x, t);
440 }
441
455 MatrixJH Jh_x(VectorF const & x, Real const t) const override
456 {
457 return this->m_Jh_x(x, t);
458 }
459
467 bool in_domain(VectorF const & x, Real const t) const override
468 {
469 return this->m_in_domain(x, t);
470 }
471
472 }; // class SemiExplicitWrapper
473
474} // namespace Sandals
475
476#endif // SANDALS_SEMIEXPLICIT_SYSTEM_HH
#define SANDALS_ASSERT(COND, MSG)
Definition Sandals.hh:44
Explicit(Type t_type, std::string t_name)
Definition Explicit.hh:55
typename Implicit< Real, N, M >::Type Type
Definition Explicit.hh:47
typename Implicit< Real, N, M >::MatrixJF MatrixJF
Definition Explicit.hh:46
typename Implicit< Real, N, M >::VectorF VectorF
Definition Explicit.hh:45
Eigen::Vector< Real, M > VectorH
Definition Implicit.hh:48
Eigen::Matrix< Real, M, N > MatrixJH
Definition Implicit.hh:49
virtual MatrixA A(VectorF const &x, Real const t) const =0
virtual TensorTA TA_x(VectorF const &x, Real const t) const =0
SemiExplicit(std::string t_name)
Definition SemiExplicit.hh:67
MatrixJF Jf_x(VectorF const &x, Real const t) const override
Definition SemiExplicit.hh:188
typename Explicit< Real, N, M >::MatrixJF MatrixJB
Definition SemiExplicit.hh:51
typename Explicit< Real, N, M >::Type Type
Definition SemiExplicit.hh:52
MatrixJF Jf_x(VectorF const &x, VectorF const &x_dot, Real const t) const
Definition SemiExplicit.hh:163
SemiExplicit()
Definition SemiExplicit.hh:61
typename Explicit< Real, N, M >::VectorF VectorF
Definition SemiExplicit.hh:46
typename Explicit< Real, N, M >::MatrixJF MatrixJF
Definition SemiExplicit.hh:47
std::shared_ptr< SemiExplicit< Real, N, M > > Pointer
Definition SemiExplicit.hh:45
MatrixJF JF_x_dot(VectorF const &x, VectorF const &, Real const t) const override
Definition SemiExplicit.hh:124
VectorF F(VectorF const &x, VectorF const &x_dot, Real const t) const override
Definition SemiExplicit.hh:82
Eigen::FullPivLU< MatrixA > m_lu
Definition SemiExplicit.hh:55
virtual MatrixJB Jb_x(VectorF const &x, Real const t) const =0
typename Explicit< Real, N, M >::VectorF VectorB
Definition SemiExplicit.hh:50
VectorF f(VectorF const &x, Real const t) const override
Definition SemiExplicit.hh:140
MatrixJF JF_x(VectorF const &x, VectorF const &x_dot, Real const t) const override
Definition SemiExplicit.hh:104
typename Explicit< Real, N, M >::MatrixJF MatrixA
Definition SemiExplicit.hh:48
typename std::vector< MatrixJF > TensorTA
Definition SemiExplicit.hh:49
virtual VectorB b(VectorF const &x, Real const t) const =0
VectorH h(VectorF const &x, Real const t) const override
Definition SemiExplicit.hh:437
FunctionTA & TA_x()
Definition SemiExplicit.hh:340
std::function< MatrixJB(VectorF const &, Real const)> FunctionJB
Definition SemiExplicit.hh:273
FunctionID m_in_domain
Definition SemiExplicit.hh:289
std::shared_ptr< SemiExplicitWrapper< Real, N, M > > Pointer
Definition SemiExplicit.hh:262
FunctionJB m_Jb_x
Definition SemiExplicit.hh:286
FunctionJH & Jh_x()
Definition SemiExplicit.hh:364
static const FunctionH DefaultH
Definition SemiExplicit.hh:278
bool in_domain(VectorF const &x, Real const t) const override
Definition SemiExplicit.hh:467
std::function< bool(VectorF const &, Real const)> FunctionID
Definition SemiExplicit.hh:276
std::function< VectorH(VectorF const &, Real const)> FunctionH
Definition SemiExplicit.hh:274
VectorB b(VectorF const &x, Real const t) const override
Definition SemiExplicit.hh:407
FunctionB & b()
Definition SemiExplicit.hh:346
~SemiExplicitWrapper()
Definition SemiExplicit.hh:328
MatrixJH Jh_x(VectorF const &x, Real const t) const override
Definition SemiExplicit.hh:455
std::function< VectorB(VectorF const &, Real const)> FunctionB
Definition SemiExplicit.hh:272
std::function< TensorTA(VectorF const &, Real const)> FunctionTA
Definition SemiExplicit.hh:271
FunctionJB & Jb_x()
Definition SemiExplicit.hh:352
std::function< MatrixA(VectorF const &, Real const)> FunctionA
Definition SemiExplicit.hh:270
static const FunctionJH DefaultJH
Definition SemiExplicit.hh:279
FunctionB m_b
Definition SemiExplicit.hh:285
MatrixA A(VectorF const &x, Real const t) const override
Definition SemiExplicit.hh:378
static const FunctionID DefaultID
Definition SemiExplicit.hh:280
SemiExplicitWrapper(FunctionA t_A, FunctionTA t_TA_x, FunctionB t_b, FunctionJB t_Jb_x, FunctionH t_h=DefaultH, FunctionJH t_Jh_x=DefaultJH, FunctionID t_in_domain=DefaultID)
Definition SemiExplicit.hh:302
FunctionID & in_domain()
Definition SemiExplicit.hh:370
FunctionJH m_Jh_x
Definition SemiExplicit.hh:288
FunctionH m_h
Definition SemiExplicit.hh:287
FunctionTA m_TA_x
Definition SemiExplicit.hh:284
SemiExplicitWrapper(std::string t_name, FunctionA t_A, FunctionTA t_TA_x, FunctionB t_b, FunctionJB t_Jb_x, FunctionH t_h=DefaultH, FunctionJH t_Jh_x=DefaultJH, FunctionID t_in_domain=DefaultID)
Definition SemiExplicit.hh:319
std::function< MatrixJH(VectorF const &, Real const)> FunctionJH
Definition SemiExplicit.hh:275
TensorTA TA_x(VectorF const &x, Real const t) const override
Definition SemiExplicit.hh:396
FunctionA m_A
Definition SemiExplicit.hh:283
FunctionA & A()
Definition SemiExplicit.hh:334
FunctionH & h()
Definition SemiExplicit.hh:358
MatrixJB Jb_x(VectorF const &x, Real const t) const override
Definition SemiExplicit.hh:425
The namespace for the Sandals library.
Definition Sandals.hh:89
SANDALS_DEFAULT_INTEGER_TYPE Integer
The Integer type as used for the API.
Definition Sandals.hh:97