1 //! @file Func1.cpp
2
3 // This file is part of Cantera. See License.txt in the top-level directory or
4 // at https://cantera.org/license.txt for license and copyright information.
5
6 #include "cantera/numerics/Func1.h"
7 #include "cantera/base/stringUtils.h"
8
9 using namespace std;
10
11 namespace Cantera
12 {
13
Func1()14 Func1::Func1() :
15 m_c(0.0),
16 m_f1(0),
17 m_f2(0),
18 m_parent(0)
19 {
20 }
21
Func1(const Func1 & right)22 Func1::Func1(const Func1& right) :
23 m_c(right.m_c),
24 m_f1(right.m_f1),
25 m_f2(right.m_f2),
26 m_parent(right.m_parent)
27 {
28 }
29
operator =(const Func1 & right)30 Func1& Func1::operator=(const Func1& right)
31 {
32 if (&right == this) {
33 return *this;
34 }
35 m_c = right.m_c;
36 m_f1 = right.m_f1;
37 m_f2 = right.m_f2;
38 m_parent = right.m_parent;
39 return *this;
40 }
41
duplicate() const42 Func1& Func1::duplicate() const
43 {
44 Func1* nfunc = new Func1(*this);
45 return *nfunc;
46 }
47
ID() const48 int Func1::ID() const
49 {
50 return 0;
51 }
52
53 // Calls method eval to evaluate the function
operator ()(doublereal t) const54 doublereal Func1::operator()(doublereal t) const
55 {
56 return eval(t);
57 }
58
59 // Evaluate the function.
eval(doublereal t) const60 doublereal Func1::eval(doublereal t) const
61 {
62 return 0.0;
63 }
64
derivative() const65 Func1& Func1::derivative() const
66 {
67 cout << "derivative error... ERR: ID = " << ID() << endl;
68 cout << write("x") << endl;
69 return *(new Func1);
70 }
71
isIdentical(Func1 & other) const72 bool Func1::isIdentical(Func1& other) const
73 {
74 if (ID() != other.ID() || m_c != other.m_c) {
75 return false;
76 }
77 if (m_f1) {
78 if (!other.m_f1) {
79 return false;
80 }
81 if (!m_f1->isIdentical(*other.m_f1)) {
82 return false;
83 }
84 }
85 if (m_f2) {
86 if (!other.m_f2) {
87 return false;
88 }
89 if (!m_f2->isIdentical(*other.m_f2)) {
90 return false;
91 }
92 }
93 return true;
94 }
95
96 //! accessor function for the returned constant
c() const97 doublereal Func1::c() const
98 {
99 return m_c;
100 }
101
102 // Function to set the stored constant
setC(doublereal c)103 void Func1::setC(doublereal c)
104 {
105 m_c = c;
106 }
107
108 //! accessor function for m_f1
func1() const109 Func1& Func1::func1() const
110 {
111 return *m_f1;
112 }
113
func2() const114 Func1& Func1::func2() const
115 {
116 return *m_f2;
117 }
118
order() const119 int Func1::order() const
120 {
121 return 3;
122 }
123
func1_dup() const124 Func1& Func1::func1_dup() const
125 {
126 return m_f1->duplicate();
127 }
128
func2_dup() const129 Func1& Func1::func2_dup() const
130 {
131 return m_f2->duplicate();
132 }
133
parent() const134 Func1* Func1::parent() const
135 {
136 return m_parent;
137 }
138
setParent(Func1 * p)139 void Func1::setParent(Func1* p)
140 {
141 m_parent = p;
142 }
143
144 /*****************************************************************************/
145
write(const string & arg) const146 string Sin1::write(const string& arg) const
147 {
148 if (m_c == 1.0) {
149 return fmt::format("\\sin({})", arg);
150 } else {
151 return fmt::format("\\sin({}{})", m_c, arg);
152 }
153 }
154
derivative() const155 Func1& Sin1::derivative() const
156 {
157 Func1* c = new Cos1(m_c);
158 Func1* r = &newTimesConstFunction(*c, m_c);
159 return *r;
160 }
161 /*****************************************************************************/
162
derivative() const163 Func1& Cos1::derivative() const
164 {
165 Func1* s = new Sin1(m_c);
166 Func1* r = &newTimesConstFunction(*s, -m_c);
167 return *r;
168 }
169
write(const std::string & arg) const170 std::string Cos1::write(const std::string& arg) const
171 {
172 if (m_c == 1.0) {
173 return fmt::format("\\cos({})", arg);
174 } else {
175 return fmt::format("\\cos({}{})", m_c, arg);
176 }
177 }
178
179 /**************************************************************************/
180
derivative() const181 Func1& Exp1::derivative() const
182 {
183 Func1* f = new Exp1(m_c);
184 if (m_c != 1.0) {
185 return newTimesConstFunction(*f, m_c);
186 } else {
187 return *f;
188 }
189 }
190
write(const std::string & arg) const191 std::string Exp1::write(const std::string& arg) const
192 {
193 if (m_c == 1.0) {
194 return fmt::format("\\exp({})", arg);
195 } else {
196 return fmt::format("\\exp({}{})", m_c, arg);
197 }
198 }
199
200 /******************************************************************************/
201
derivative() const202 Func1& Pow1::derivative() const
203 {
204 Func1* r;
205 if (m_c == 0.0) {
206 r = new Const1(0.0);
207 } else if (m_c == 1.0) {
208 r = new Const1(1.0);
209 } else {
210 Func1* f = new Pow1(m_c - 1.0);
211 r = &newTimesConstFunction(*f, m_c);
212 }
213 return *r;
214 }
215
216 /******************************************************************************/
217
Tabulated1(size_t n,const double * tvals,const double * fvals,const std::string & method)218 Tabulated1::Tabulated1(size_t n, const double* tvals, const double* fvals,
219 const std::string& method) :
220 Func1() {
221 m_tvec.resize(n);
222 std::copy(tvals, tvals + n, m_tvec.begin());
223 for (auto it = std::begin(m_tvec) + 1; it != std::end(m_tvec); it++) {
224 if (*(it - 1) > *it) {
225 throw CanteraError("Tabulated1::Tabulated1",
226 "time values are not increasing monotonically.");
227 }
228 }
229 m_fvec.resize(n);
230 std::copy(fvals, fvals + n, m_fvec.begin());
231 if (method == "linear") {
232 m_isLinear = true;
233 } else if (method == "previous") {
234 m_isLinear = false;
235 } else {
236 throw CanteraError("Tabulated1::Tabulated1",
237 "interpolation method '{}' is not implemented",
238 method);
239 }
240 }
241
eval(double t) const242 double Tabulated1::eval(double t) const {
243 size_t siz = m_tvec.size();
244 // constructor ensures that siz > 0
245 if (t <= m_tvec[0]) {
246 return m_fvec[0];
247 } else if (t >= m_tvec[siz-1]) {
248 return m_fvec[siz-1];
249 } else {
250 size_t ix = 0;
251 while (t > m_tvec[ix+1]) {
252 ix++;
253 }
254 if (m_isLinear) {
255 double df = m_fvec[ix+1] - m_fvec[ix];
256 df /= m_tvec[ix+1] - m_tvec[ix];
257 df *= t - m_tvec[ix];
258 return m_fvec[ix] + df;
259 } else {
260 return m_fvec[ix];
261 }
262 }
263 }
264
derivative() const265 Func1& Tabulated1::derivative() const {
266 vector_fp tvec;
267 vector_fp dvec;
268 size_t siz = m_tvec.size();
269 if (m_isLinear) {
270 // piece-wise continuous derivative
271 if (siz>1) {
272 for (size_t i=1; i<siz; i++) {
273 double d = (m_fvec[i] - m_fvec[i-1]) /
274 (m_tvec[i] - m_tvec[i-1]);
275 tvec.push_back(m_tvec[i-1]);
276 dvec.push_back(d);
277 }
278 }
279 tvec.push_back(m_tvec[siz-1]);
280 dvec.push_back(0.);
281 } else {
282 // derivative is zero (ignoring discontinuities)
283 tvec.push_back(m_tvec[0]);
284 tvec.push_back(m_tvec[siz-1]);
285 dvec.push_back(0.);
286 dvec.push_back(0.);
287 }
288 return *(new Tabulated1(tvec.size(), &tvec[0], &dvec[0], "previous"));
289 }
290
291 /******************************************************************************/
292
write(const std::string & arg) const293 string Func1::write(const std::string& arg) const
294 {
295 return fmt::format("<unknown {}>({})", ID(), arg);
296 }
297
write(const std::string & arg) const298 string Pow1::write(const std::string& arg) const
299 {
300 if (m_c == 0.5) {
301 return "\\sqrt{" + arg + "}";
302 }
303 if (m_c == -0.5) {
304 return "\\frac{1}{\\sqrt{" + arg + "}}";
305 }
306 if (m_c != 1.0) {
307 return fmt::format("\\left({}\\right)^{{{}}}", arg, m_c);
308 } else {
309 return arg;
310 }
311 }
312
write(const std::string & arg) const313 string Tabulated1::write(const std::string& arg) const
314 {
315 return fmt::format("\\mathrm{{Tabulated}}({})", arg);
316 }
317
write(const std::string & arg) const318 string Const1::write(const std::string& arg) const
319 {
320 return fmt::format("{}", m_c);
321 }
322
write(const std::string & arg) const323 string Ratio1::write(const std::string& arg) const
324 {
325 return "\\frac{" + m_f1->write(arg) + "}{"
326 + m_f2->write(arg) + "}";
327 }
328
write(const std::string & arg) const329 string Product1::write(const std::string& arg) const
330 {
331 string s = m_f1->write(arg);
332 if (m_f1->order() < order()) {
333 s = "\\left(" + s + "\\right)";
334 }
335 string s2 = m_f2->write(arg);
336 if (m_f2->order() < order()) {
337 s2 = "\\left(" + s2 + "\\right)";
338 }
339 return s + " " + s2;
340 }
341
write(const std::string & arg) const342 string Sum1::write(const std::string& arg) const
343 {
344 string s1 = m_f1->write(arg);
345 string s2 = m_f2->write(arg);
346 if (s2[0] == '-') {
347 return s1 + " - " + s2.substr(1,s2.size());
348 } else {
349 return s1 + " + " + s2;
350 }
351 }
352
write(const std::string & arg) const353 string Diff1::write(const std::string& arg) const
354 {
355 string s1 = m_f1->write(arg);
356 string s2 = m_f2->write(arg);
357 if (s2[0] == '-') {
358 return s1 + " + " + s2.substr(1,s2.size());
359 } else {
360 return s1 + " - " + s2;
361 }
362 }
363
write(const std::string & arg) const364 string Composite1::write(const std::string& arg) const
365 {
366 string g = m_f2->write(arg);
367 return m_f1->write(g);
368 }
369
write(const std::string & arg) const370 string TimesConstant1::write(const std::string& arg) const
371 {
372 string s = m_f1->write(arg);
373 if (m_f1->order() < order()) {
374 s = "\\left(" + s + "\\right)";
375 }
376 if (m_c == 1.0) {
377 return s;
378 }
379 if (m_c == -1.0) {
380 return "-"+s;
381 }
382 char n = s[0];
383 if (n >= '0' && n <= '9') {
384 s = "\\left(" + s + "\\right)";
385 }
386 return fmt::format("{}{}", m_c, s);
387 }
388
write(const std::string & arg) const389 string PlusConstant1::write(const std::string& arg) const
390 {
391 if (m_c == 0.0) {
392 return m_f1->write(arg);
393 }
394 return fmt::format("{} + {}", m_f1->write(arg), m_c);
395 }
396
isProportional(TimesConstant1 & other)397 doublereal Func1::isProportional(TimesConstant1& other)
398 {
399 if (isIdentical(other.func1())) {
400 return other.c();
401 }
402 return 0.0;
403 }
isProportional(Func1 & other)404 doublereal Func1::isProportional(Func1& other)
405 {
406 if (isIdentical(other)) {
407 return 1.0;
408 } else {
409 return 0.0;
410 }
411 }
412
isConstant(Func1 & f)413 static bool isConstant(Func1& f)
414 {
415 if (f.ID() == ConstFuncType) {
416 return true;
417 } else {
418 return false;
419 }
420 }
421
isZero(Func1 & f)422 static bool isZero(Func1& f)
423 {
424 if (f.ID() == ConstFuncType && f.c() == 0.0) {
425 return true;
426 } else {
427 return false;
428 }
429 }
430
isOne(Func1 & f)431 static bool isOne(Func1& f)
432 {
433 if (f.ID() == ConstFuncType && f.c() == 1.0) {
434 return true;
435 } else {
436 return false;
437 }
438 }
439
isTimesConst(Func1 & f)440 static bool isTimesConst(Func1& f)
441 {
442 if (f.ID() == TimesConstantFuncType) {
443 return true;
444 } else {
445 return false;
446 }
447 }
448
isExp(Func1 & f)449 static bool isExp(Func1& f)
450 {
451 if (f.ID() == ExpFuncType) {
452 return true;
453 } else {
454 return false;
455 }
456 }
457
isPow(Func1 & f)458 static bool isPow(Func1& f)
459 {
460 if (f.ID() == PowFuncType) {
461 return true;
462 } else {
463 return false;
464 }
465 }
466
newSumFunction(Func1 & f1,Func1 & f2)467 Func1& newSumFunction(Func1& f1, Func1& f2)
468 {
469 if (f1.isIdentical(f2)) {
470 return newTimesConstFunction(f1, 2.0);
471 }
472 if (isZero(f1)) {
473 delete &f1;
474 return f2;
475 }
476 if (isZero(f2)) {
477 delete &f2;
478 return f1;
479 }
480 doublereal c = f1.isProportional(f2);
481 if (c != 0) {
482 if (c == -1.0) {
483 return *(new Const1(0.0));
484 } else {
485 return newTimesConstFunction(f1, c + 1.0);
486 }
487 }
488 return *(new Sum1(f1, f2));
489 }
490
newDiffFunction(Func1 & f1,Func1 & f2)491 Func1& newDiffFunction(Func1& f1, Func1& f2)
492 {
493 if (isZero(f2)) {
494 delete &f2;
495 return f1;
496 }
497 if (f1.isIdentical(f2)) {
498 delete &f1;
499 delete &f2;
500 return *(new Const1(0.0));
501 }
502 doublereal c = f1.isProportional(f2);
503 if (c != 0.0) {
504 if (c == 1.0) {
505 return *(new Const1(0.0));
506 } else {
507 return newTimesConstFunction(f1, 1.0 - c);
508 }
509 }
510 return *(new Diff1(f1, f2));
511 }
512
newProdFunction(Func1 & f1,Func1 & f2)513 Func1& newProdFunction(Func1& f1, Func1& f2)
514 {
515 if (isOne(f1)) {
516 delete &f1;
517 return f2;
518 }
519 if (isOne(f2)) {
520 delete &f2;
521 return f1;
522 }
523 if (isZero(f1) || isZero(f2)) {
524 delete &f1;
525 delete &f2;
526 return *(new Const1(0.0));
527 }
528 if (isConstant(f1) && isConstant(f2)) {
529 doublereal c1c2 = f1.c() * f2.c();
530 delete &f1;
531 delete &f2;
532 return *(new Const1(c1c2));
533 }
534 if (isConstant(f1)) {
535 doublereal c = f1.c();
536 delete &f1;
537 return newTimesConstFunction(f2, c);
538 }
539 if (isConstant(f2)) {
540 doublereal c = f2.c();
541 delete &f2;
542 return newTimesConstFunction(f1, c);
543 }
544
545 if (isPow(f1) && isPow(f2)) {
546 Func1& p = *(new Pow1(f1.c() + f2.c()));
547 delete &f1;
548 delete &f2;
549 return p;
550 }
551
552 if (isExp(f1) && isExp(f2)) {
553 Func1& p = *(new Exp1(f1.c() + f2.c()));
554 delete &f1;
555 delete &f2;
556 return p;
557 }
558
559 bool tc1 = isTimesConst(f1);
560 bool tc2 = isTimesConst(f2);
561
562 if (tc1 || tc2) {
563 doublereal c1 = 1.0, c2 = 1.0;
564 Func1* ff1 = 0, *ff2 = 0;
565 if (tc1) {
566 c1 = f1.c();
567 ff1 = &f1.func1_dup();
568 delete &f1;
569 } else {
570 ff1 = &f1;
571 }
572 if (tc2) {
573 c2 = f2.c();
574 ff2 = &f2.func1_dup();
575 delete &f2;
576 } else {
577 ff2 = &f2;
578 }
579 Func1& p = newProdFunction(*ff1, *ff2);
580
581 if (c1* c2 != 1.0) {
582 return newTimesConstFunction(p, c1*c2);
583 } else {
584 return p;
585 }
586 } else {
587 return *(new Product1(f1, f2));
588 }
589 }
590
newRatioFunction(Func1 & f1,Func1 & f2)591 Func1& newRatioFunction(Func1& f1, Func1& f2)
592 {
593 if (isOne(f2)) {
594 return f1;
595 }
596 if (isZero(f1)) {
597 return *(new Const1(0.0));
598 }
599 if (f1.isIdentical(f2)) {
600 delete &f1;
601 delete &f2;
602 return *(new Const1(1.0));
603 }
604 if (f1.ID() == PowFuncType && f2.ID() == PowFuncType) {
605 return *(new Pow1(f1.c() - f2.c()));
606 }
607 if (f1.ID() == ExpFuncType && f2.ID() == ExpFuncType) {
608 return *(new Exp1(f1.c() - f2.c()));
609 }
610 return *(new Ratio1(f1, f2));
611 }
612
newCompositeFunction(Func1 & f1,Func1 & f2)613 Func1& newCompositeFunction(Func1& f1, Func1& f2)
614 {
615 if (isZero(f1)) {
616 delete &f1;
617 delete &f2;
618 return *(new Const1(0.0));
619 }
620 if (isConstant(f1)) {
621 delete &f2;
622 return f1;
623 }
624 if (isPow(f1) && f1.c() == 1.0) {
625 delete &f1;
626 return f2;
627 }
628 if (isPow(f1) && f1.c() == 0.0) {
629 delete &f1;
630 delete &f2;
631 return *(new Const1(1.0));
632 }
633 if (isPow(f1) && isPow(f2)) {
634 doublereal c1c2 = f1.c() * f2.c();
635 delete &f1;
636 delete &f2;
637 return *(new Pow1(c1c2));
638 }
639 return *(new Composite1(f1, f2));
640 }
641
newTimesConstFunction(Func1 & f,doublereal c)642 Func1& newTimesConstFunction(Func1& f, doublereal c)
643 {
644 if (c == 0.0) {
645 delete &f;
646 return *(new Const1(0.0));
647 }
648 if (c == 1.0) {
649 return f;
650 }
651 if (f.ID() == TimesConstantFuncType) {
652 f.setC(f.c() * c);
653 return f;
654 }
655 return *(new TimesConstant1(f, c));
656 }
657
newPlusConstFunction(Func1 & f,doublereal c)658 Func1& newPlusConstFunction(Func1& f, doublereal c)
659 {
660 if (c == 0.0) {
661 return f;
662 }
663 if (isConstant(f)) {
664 doublereal cc = f.c() + c;
665 delete &f;
666 return *(new Const1(cc));
667 }
668 if (f.ID() == PlusConstantFuncType) {
669 f.setC(f.c() + c);
670 return f;
671 }
672 return *(new PlusConstant1(f, c));
673 }
674
675 }
676