1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/base/ScalarFunctionsImpl.cc,v 1.56 2017/01/12 14:46:08 masarati Exp $ */
2 /*
3 * HmFe (C) is a FEM analysis code.
4 *
5 * Copyright (C) 1996-2017
6 *
7 * Marco Morandini <morandini@aero.polimi.it>
8 *
9 * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
10 * via La Masa, 34 - 20156 Milano, Italy
11 * http://www.aero.polimi.it
12 *
13 * Changing this copyright notice is forbidden.
14 * This program is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
17 *
18 */
19 /*
20 * MBDyn (C) is a multibody analysis code.
21 * http://www.mbdyn.org
22 *
23 * Copyright (C) 1996-2017
24 *
25 * This code is a partial merge of HmFe and MBDyn.
26 *
27 * Pierangelo Masarati <masarati@aero.polimi.it>
28 * Paolo Mantegazza <mantegazza@aero.polimi.it>
29 *
30 * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
31 * via La Masa, 34 - 20156 Milano, Italy
32 * http://www.aero.polimi.it
33 *
34 * Changing this copyright notice is forbidden.
35 *
36 * This program is free software; you can redistribute it and/or modify
37 * it under the terms of the GNU General Public License as published by
38 * the Free Software Foundation (version 2 of the License).
39 *
40 *
41 * This program is distributed in the hope that it will be useful,
42 * but WITHOUT ANY WARRANTY; without even the implied warranty of
43 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
44 * GNU General Public License for more details.
45 *
46 * You should have received a copy of the GNU General Public License
47 * along with this program; if not, write to the Free Software
48 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
49 */
50 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
51
52 #include <cmath>
53 #include <typeinfo>
54
55 #include "myassert.h"
56
57 #include "ScalarFunctionsImpl.h"
58 #include "interp.h"
59
60 #include "mbpar.h"
61 #include "dataman.h"
62
63 //ptr_cast: I use an helper stuct to avoid troubles
64 //to be re-written when compilers will
65 //support partial specialization of functions template
66
67 //notare che non ho StlOverload.C dove metto il codice dei metodi statici,
68 //in teoria cosi' non e' corretto
69 template<class T1, class T2> struct ptr_cast_helper {
70 static T1 cast(T2&);
71 };
72
73 template<class T1, class T2> struct ptr_cast_helper<T1*, T2> {
castptr_cast_helper74 static T1* cast(T2&arg) {
75 T1* out = dynamic_cast<T1*>(arg);
76 if (out == 0) {
77 #ifdef MYDEBUG
78 std::cerr << "ptr_cast error" << std::endl;
79 #endif
80 throw std::bad_cast();
81 }
82 return out;
83 };
84 };
85
86 template<class T1, class T2> struct ptr_cast_helper<T1*const, T2> {
castptr_cast_helper87 static T1*const cast(T2&arg) {
88 T1*const out = dynamic_cast<T1*const>(arg);
89 if (out == 0) {
90 #ifdef MYDEBUG
91 std::cerr << "ptr_cast error" << std::endl;
92 #endif
93 throw std::bad_cast();
94 }
95 return out;
96 };
97 };
98
99 template<class T1, class T2> struct ptr_cast_helper<const T1*, T2> {
castptr_cast_helper100 static const T1* cast(T2&arg) {
101 const T1* out = dynamic_cast<const T1*>(arg);
102 if (out == 0) {
103 #ifdef MYDEBUG
104 std::cerr << "ptr_cast error" << std::endl;
105 #endif
106 throw std::bad_cast();
107 }
108 return out;
109 };
110 };
111
112 template<class T1, class T2> struct ptr_cast_helper<const T1*const, T2> {
castptr_cast_helper113 static const T1*const cast(T2&arg) {
114 const T1*const out = dynamic_cast<const T1*const>(arg);
115 if (out == 0) {
116 #ifdef MYDEBUG
117 std::cerr << "ptr_cast error" << std::endl;
118 #endif
119 throw std::bad_cast();
120 }
121 return out;
122 };
123 };
124
125 template<class T1,class T2>
126 T1 ptr_cast(T2& arg) {
127 return ptr_cast_helper<T1,T2>::cast(arg);
128 }
129
130 // BasicScalarFunction
~BasicScalarFunction(void)131 BasicScalarFunction::~BasicScalarFunction(void)
132 {
133 NO_OP;
134 }
135
136 // DifferentiableScalarFunction
~DifferentiableScalarFunction(void)137 DifferentiableScalarFunction::~DifferentiableScalarFunction(void)
138 {
139 NO_OP;
140 }
141
142 // ConstScalarFunction
ConstScalarFunction(const doublereal v)143 ConstScalarFunction::ConstScalarFunction(const doublereal v)
144 : y(v)
145 {
146 NO_OP;
147 }
148
~ConstScalarFunction(void)149 ConstScalarFunction::~ConstScalarFunction(void)
150 {
151 NO_OP;
152 }
153
154 doublereal
operator ()(const doublereal x) const155 ConstScalarFunction::operator()(const doublereal x) const
156 {
157 return y;
158 }
159
160 doublereal
ComputeDiff(const doublereal x,const integer order) const161 ConstScalarFunction::ComputeDiff(const doublereal x, const integer order) const
162 {
163 ASSERTMSGBREAK(order >=0, "Error in ConstScalarFunction::ComputeDiff, order<0");
164 switch (order) {
165 case 0:
166 return this->operator()(x);
167 break;
168 default:
169 return 0;
170 break;
171 }
172 }
173
174 // ScalarFunction parsing functional object
175 struct ConstSFR: public ScalarFunctionRead {
176 virtual const BasicScalarFunction *
ReadConstSFR177 Read(DataManager* const pDM, MBDynParser& HP) const {
178 doublereal c = HP.GetReal();
179 return new ConstScalarFunction(c);
180 };
181 };
182
183 // LinearScalarFunction
LinearScalarFunction(const doublereal t_i,const doublereal y_i,const doublereal t_f,const doublereal y_f)184 LinearScalarFunction::LinearScalarFunction(
185 const doublereal t_i,
186 const doublereal y_i,
187 const doublereal t_f,
188 const doublereal y_f)
189 {
190 ASSERTMSGBREAK(t_i != t_f, "LinearScalarFunction error, t_i == t_f");
191 m = (y_f - y_i)/(t_f - t_i);
192 y0 = y_i - m*t_i;
193 }
194
~LinearScalarFunction(void)195 LinearScalarFunction::~LinearScalarFunction(void)
196 {
197 NO_OP;
198 }
199
200 doublereal
operator ()(const doublereal x) const201 LinearScalarFunction::operator()(const doublereal x) const
202 {
203 return y0 + m*x;
204 }
205
206 doublereal
ComputeDiff(const doublereal x,const integer order) const207 LinearScalarFunction::ComputeDiff(const doublereal x, const integer order) const
208 {
209 ASSERTMSGBREAK(order >=0, "Error in LinearScalarFunction::ComputeDiff, order<0");
210 switch (order) {
211 case 0:
212 return this->operator()(x);
213
214 case 1:
215 return m;
216
217 default:
218 return 0.;
219
220 }
221 }
222
223 // ScalarFunction parsing functional object
224 struct LinearSFR: public ScalarFunctionRead {
225 virtual const BasicScalarFunction *
ReadLinearSFR226 Read(DataManager* const pDM, MBDynParser& HP) const {
227 doublereal t_i = HP.GetReal();
228 doublereal y_i = HP.GetReal();
229 doublereal t_f = HP.GetReal();
230 doublereal y_f = HP.GetReal();
231 return new LinearScalarFunction(t_i,y_i,t_f,y_f);
232 };
233 };
234
235 // PowScalarFunction
PowScalarFunction(const doublereal p)236 PowScalarFunction::PowScalarFunction(const doublereal p)
237 : pw(p)
238 {
239 NO_OP;
240 }
241
~PowScalarFunction(void)242 PowScalarFunction::~PowScalarFunction(void)
243 {
244 NO_OP;
245 }
246
247 doublereal
operator ()(const doublereal x) const248 PowScalarFunction::operator()(const doublereal x) const
249 {
250 return pow(x, pw);
251 }
252
253 doublereal
ComputeDiff(const doublereal x,const integer order) const254 PowScalarFunction::ComputeDiff(const doublereal x, const integer order) const
255 {
256 ASSERTMSGBREAK(order >=0, "Error in PowScalarFunction::ComputeDiff, order<0");
257 switch (order) {
258 case 0:
259 return this->operator()(x);
260
261 default:
262 doublereal mul = 1.;
263 for (integer i = 0; i < order; i++) {
264 mul *= pw - i;
265 }
266 return mul*pow(x, pw - order);
267
268 }
269 }
270
271 // ScalarFunction parsing functional object
272 struct PowSFR: public ScalarFunctionRead {
273 virtual const BasicScalarFunction *
ReadPowSFR274 Read(DataManager* const pDM, MBDynParser& HP) const {
275 doublereal p = HP.GetReal();
276 return new PowScalarFunction(p);
277 };
278 };
279
280 // LogScalarFunction
LogScalarFunction(const doublereal & ml,const doublereal & b,const doublereal & c)281 LogScalarFunction::LogScalarFunction(
282 const doublereal& ml,
283 const doublereal& b,
284 const doublereal& c)
285 : mul_input(ml), mul_const(ml), base(b), coef(c)
286 {
287 if (b != 1.) {
288 mul_const /= log(b);
289 }
290 }
291
~LogScalarFunction()292 LogScalarFunction::~LogScalarFunction()
293 {
294 NO_OP;
295 }
296
297 doublereal
operator ()(const doublereal x) const298 LogScalarFunction::operator()(const doublereal x) const
299 {
300 doublereal xx = coef*x;
301 if (xx <= 0.) {
302 silent_cerr("LogScalarFunction: argument must be positive" << std::endl);
303 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
304 }
305
306 return log(xx)*mul_const;
307 }
308
309 doublereal
ComputeDiff(const doublereal x,const integer order) const310 LogScalarFunction::ComputeDiff(const doublereal x, const integer order) const
311 {
312 ASSERTMSGBREAK(order >= 0, "Error in LogScalarFunction::ComputeDiff, order<0");
313 if (order == 0) {
314 return this->operator()(x);
315 }
316
317 if (x <= 0.) {
318 silent_cerr("LogScalarFunction: argument must be positive" << std::endl);
319 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
320 }
321
322 /*
323 *
324
325 d^i mul_const
326 --- ( mul_const * log( coef * x ) ) = - ( -1 ) ^ i * ( i - 1 )! * ---------
327 dx^i x^i
328
329 *
330 */
331
332 doublereal d = mul_const/x;
333 for (int i = 1; i < order; i++) {
334 d *= -i/x;
335 }
336
337 return d;
338 }
339
340 // ScalarFunction parsing functional object
341 struct LogSFR: public ScalarFunctionRead {
342 virtual const BasicScalarFunction *
ReadLogSFR343 Read(DataManager* const pDM, MBDynParser& HP) const {
344 doublereal b = 1.;
345 if (HP.IsKeyWord("base")) {
346 b = HP.GetReal();
347 if (b <= 0.) {
348 silent_cerr("LogSFR: "
349 "invalid base " << b
350 << " at line "
351 << HP.GetLineData()
352 << std::endl);
353 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
354 }
355 }
356
357 doublereal c = 1.;
358 if (HP.IsKeyWord("coefficient")) {
359 c = HP.GetReal();
360 // note: c*x must be > 0, but x could be negative
361 if (c == 0.) {
362 silent_cerr("LogSFR: "
363 "invalid coefficient " << c
364 << " at line "
365 << HP.GetLineData()
366 << std::endl);
367 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
368 }
369 }
370
371 doublereal m = HP.GetReal();
372 return new LogScalarFunction(m, b, c);
373 };
374 };
375
376 // ExpScalarFunction
ExpScalarFunction(const doublereal & ml,const doublereal & b,const doublereal & c)377 ExpScalarFunction::ExpScalarFunction(
378 const doublereal& ml,
379 const doublereal& b,
380 const doublereal& c)
381 : mul(ml), base(b), coef_input(c), coef_const(c)
382 {
383 if (base != 1.) {
384 coef_const *= log(b);
385 }
386 }
387
~ExpScalarFunction()388 ExpScalarFunction::~ExpScalarFunction()
389 {
390 NO_OP;
391 }
392
393 doublereal
operator ()(const doublereal x) const394 ExpScalarFunction::operator()(const doublereal x) const
395 {
396 return exp(coef_const*x)*mul;
397 }
398
399 doublereal
ComputeDiff(const doublereal x,const integer order) const400 ExpScalarFunction::ComputeDiff(const doublereal x, const integer order) const
401 {
402 ASSERTMSGBREAK(order >= 0, "Error in ExpScalarFunction::ComputeDiff, order<0");
403
404 /*
405 *
406
407 d^i mul_const
408 --- ( mul * log( coef * x ) ) = - ( -1 ) ^ i * ( i - 1 )! * ---------
409 dx^i x^i
410
411 *
412 */
413
414 doublereal d = 1.;
415 for (int i = 0; i < order; i++) {
416 d *= coef_input;
417 }
418
419 return d * this->operator()(x);
420 }
421
422 // ScalarFunction parsing functional object
423 struct ExpSFR: public ScalarFunctionRead {
424 virtual const BasicScalarFunction *
ReadExpSFR425 Read(DataManager* const pDM, MBDynParser& HP) const {
426 doublereal b = 1.;
427 if (HP.IsKeyWord("base")) {
428 b = HP.GetReal();
429 if (b <= 0.) {
430 silent_cerr("ExpSFR: "
431 "invalid base " << b
432 << " at line "
433 << HP.GetLineData()
434 << std::endl);
435 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
436 }
437 }
438
439 doublereal c = 1.;
440 if (HP.IsKeyWord("coefficient")) {
441 c = HP.GetReal();
442 if (c == 0.) {
443 silent_cerr("ExpSFR: "
444 "invalid coefficient " << c
445 << " at line "
446 << HP.GetLineData()
447 << std::endl);
448 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
449 }
450 }
451
452 doublereal m = HP.GetReal();
453 return new ExpScalarFunction(m, b, c);
454 };
455 };
456
457 // CubicSplineScalarFunction
CubicSplineScalarFunction(const std::vector<doublereal> & y_i,const std::vector<doublereal> & x_i,bool doNotExtrapolate)458 CubicSplineScalarFunction::CubicSplineScalarFunction(
459 const std::vector<doublereal>& y_i,
460 const std::vector<doublereal>& x_i,
461 bool doNotExtrapolate)
462 : Y_i(y_i), X_i(x_i),
463 doNotExtrapolate(doNotExtrapolate)
464 {
465 ASSERTMSGBREAK(Y_i.size() == X_i.size(),
466 "CubicSplineScalarFunction error, Y_i.size() != X_i.size()");
467 std::vector<doublereal>::iterator xi, xe;
468 xi = X_i.begin();
469 xe = X_i.end() - 1;
470 for (unsigned i = 0; xi != xe; ++xi, ++i) {
471 if (*xi >= *(xi + 1)) {
472 silent_cerr("CubicSplineScalarFunction error, "
473 "X is not ordered: "
474 "X[" << i << "]=" << *xi
475 << " is not less than "
476 "X[" << i + 1 << "]=" << *(xi + 1)
477 << std::endl);
478 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
479 }
480 }
481 spline(X_i, Y_i, b, c, d);
482 }
483
~CubicSplineScalarFunction(void)484 CubicSplineScalarFunction::~CubicSplineScalarFunction(void)
485 {
486 NO_OP;
487 }
488
489 doublereal
operator ()(const doublereal x) const490 CubicSplineScalarFunction::operator()(const doublereal x) const
491 {
492 if (doNotExtrapolate) {
493 if (x <= X_i[0]) {
494 return Y_i[0];
495 }
496
497 int s = X_i.size() - 1;
498 if (x >= X_i[s]) {
499 return Y_i[s];
500 }
501 }
502
503 return seval(x, X_i, Y_i, b, c, d);
504 }
505
506 doublereal
ComputeDiff(const doublereal x,const integer order) const507 CubicSplineScalarFunction::ComputeDiff(const doublereal x, const integer order) const
508 {
509 ASSERTMSGBREAK(order >=0, "Error in CubicSplineScalarFunction::ComputeDiff, order<0");
510 switch (order) {
511 case 0:
512 return this->operator()(x);
513
514 case 1:
515 return seval(x, X_i, Y_i, b, c, d, 1);
516
517 case 2:
518 return seval(x, X_i, Y_i, b, c, d, 2);
519
520 case 3:
521 return seval(x, X_i, Y_i, b, c, d, 3);
522
523 default:
524 return 0.;
525 }
526 }
527
528 // ScalarFunction parsing functional object
529 struct CubicSplineSFR: public ScalarFunctionRead {
530 virtual const BasicScalarFunction *
ReadCubicSplineSFR531 Read(DataManager* const pDM, MBDynParser& HP) const {
532 bool doNotExtrapolate(false);
533 if (HP.IsKeyWord("do" "not" "extrapolate")) {
534 doNotExtrapolate = true;
535 }
536 std::vector<doublereal> y_i;
537 std::vector<doublereal> x_i;
538 y_i.resize(3);
539 x_i.resize(3);
540 for (int i=0; i<3; i++) {
541 x_i[i] = HP.GetReal();
542 y_i[i] = HP.GetReal();
543 }
544 while (HP.IsArg() && !HP.IsKeyWord("end")) {
545 int size = x_i.size();
546 x_i.resize(size+1);
547 y_i.resize(size+1);
548 x_i[size] = HP.GetReal();
549 y_i[size] = HP.GetReal();
550 }
551 return new CubicSplineScalarFunction(y_i, x_i,
552 doNotExtrapolate);
553 };
554 };
555
556 // MultiLinearScalarFunction
MultiLinearScalarFunction(const std::vector<doublereal> & y_i,const std::vector<doublereal> & x_i,bool doNotExtrapolate)557 MultiLinearScalarFunction::MultiLinearScalarFunction(
558 const std::vector<doublereal>& y_i,
559 const std::vector<doublereal>& x_i,
560 bool doNotExtrapolate)
561 : Y_i(y_i), X_i(x_i),
562 doNotExtrapolate(doNotExtrapolate)
563 {
564 ASSERTMSGBREAK(X_i.size() == Y_i.size(),
565 "MultiLinearScalarFunction error, Y_i.size() != X_i.size()");
566 std::vector<doublereal>::iterator xi, xe;
567 xi = X_i.begin();
568 xe = X_i.end()-1;
569 for (unsigned i = 0; xi != xe; ++xi, ++i) {
570 if (*xi >= *(xi + 1)) {
571 silent_cerr("MultiLinearScalarFunction error, "
572 "X is not ordered: "
573 "X[" << i << "]=" << *xi
574 << " is not less than "
575 "X[" << i + 1 << "]=" << *(xi + 1)
576 << std::endl);
577 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
578 }
579 }
580 }
581
~MultiLinearScalarFunction(void)582 MultiLinearScalarFunction::~MultiLinearScalarFunction(void)
583 {
584 NO_OP;
585 }
586
587 doublereal
operator ()(const doublereal x) const588 MultiLinearScalarFunction::operator()(const doublereal x) const
589 {
590 if (doNotExtrapolate) {
591 if (x <= X_i[0]) {
592 return Y_i[0];
593 }
594
595 int s = X_i.size() - 1;
596 if (x >= X_i[s]) {
597 return Y_i[s];
598 }
599 }
600
601 return leval(x, X_i, Y_i);
602 }
603
604 doublereal
ComputeDiff(const doublereal x,const integer order) const605 MultiLinearScalarFunction::ComputeDiff(const doublereal x, const integer order) const
606 {
607 ASSERTMSGBREAK(order >=0, "Error in MultiLinearScalarFunction::ComputeDiff, order<0");
608 switch (order) {
609 case 0:
610 return operator()(x);
611
612 case 1:
613 return leval(x, X_i, Y_i, order);
614
615 default:
616 return 0.;
617 }
618 }
619
620 // ScalarFunction parsing functional object
621 struct MultiLinearSFR: public ScalarFunctionRead {
622 virtual const BasicScalarFunction *
ReadMultiLinearSFR623 Read(DataManager* const pDM, MBDynParser& HP) const {
624 bool doNotExtrapolate(false);
625 if (HP.IsKeyWord("do" "not" "extrapolate")) {
626 doNotExtrapolate = true;
627 }
628 std::vector<doublereal> y_i;
629 std::vector<doublereal> x_i;
630 y_i.resize(2);
631 x_i.resize(2);
632 for (int i=0; i<2; i++) {
633 x_i[i] = HP.GetReal();
634 y_i[i] = HP.GetReal();
635 }
636 while (HP.IsArg() && !HP.IsKeyWord("end")) {
637 int size = x_i.size();
638 x_i.resize(size+1);
639 y_i.resize(size+1);
640 x_i[size] = HP.GetReal();
641 y_i[size] = HP.GetReal();
642 }
643 return new MultiLinearScalarFunction(y_i, x_i,
644 doNotExtrapolate);
645 };
646 };
647
648 // ChebychevScalarFunction
ChebychevScalarFunction(const std::vector<doublereal> & v,const doublereal & a,const doublereal & b,bool dne)649 ChebychevScalarFunction::ChebychevScalarFunction(
650 const std::vector<doublereal>& v,
651 const doublereal& a, const doublereal& b, bool dne)
652 : vCoef(v), da(a), dfa(0.), dfap(0.), db(b), dfb(0.), dfbp(0.),
653 doNotExtrapolate(dne)
654 {
655 if (!doNotExtrapolate) {
656 const_cast<doublereal&>(dfa) = this->operator()(a);
657 const_cast<doublereal&>(dfap) = this->ComputeDiff(a);
658 const_cast<doublereal&>(dfb) = this->operator()(b);
659 const_cast<doublereal&>(dfbp) = this->ComputeDiff(b);
660 }
661 }
662
~ChebychevScalarFunction(void)663 ChebychevScalarFunction::~ChebychevScalarFunction(void)
664 {
665 NO_OP;
666 }
667
668 doublereal
operator ()(const doublereal x) const669 ChebychevScalarFunction::operator()(const doublereal x) const
670 {
671 if (x < da) {
672 if (doNotExtrapolate) {
673 silent_cerr("Chebychev interpolation: "
674 "x=" << x << " is out of range "
675 "[" << da << "," << db << "]" << std::endl);
676 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
677 }
678
679 return dfa + dfap*(x - da);
680
681 } else if (x > db) {
682 if (doNotExtrapolate) {
683 silent_cerr("Chebychev interpolation: "
684 "x=" << x << " is out of range "
685 "[" << da << "," << db << "]" << std::endl);
686 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
687 }
688
689 return dfb + dfbp*(x - db);
690 }
691
692 doublereal xi = (2.*x - (da + db))/(db - da);
693 doublereal d[2] = { 1., xi };
694 doublereal val = vCoef[0] + vCoef[1]*xi;
695
696 for (unsigned i = 2; i < vCoef.size(); i++) {
697 doublereal Tx = 2.*xi*d[1 - i%2] - d[i%2];
698 val += vCoef[i]*Tx;
699 d[i%2] = Tx;
700 }
701
702 return val;
703 }
704
705 doublereal
ComputeDiff(const doublereal x,const integer order) const706 ChebychevScalarFunction::ComputeDiff(const doublereal x, const integer order) const
707 {
708 ASSERTMSGBREAK(order >=0, "Error in ChebychevScalarFunction::ComputeDiff, order<0");
709
710 switch (order) {
711 case 0:
712 return operator()(x);
713
714 case 1:
715 break;
716
717 default:
718 silent_cerr("differentiation of order " << order << " not supported yet" << std::endl);
719 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
720 }
721
722 if (x < da) {
723 if (doNotExtrapolate) {
724 silent_cerr("Chebychev interpolation: "
725 "x=" << x << " is out of range "
726 "[" << da << "," << db << "]" << std::endl);
727 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
728 }
729
730 return dfap;
731
732 } else if (x > db) {
733 if (doNotExtrapolate) {
734 silent_cerr("Chebychev interpolation: "
735 "x=" << x << " is out of range "
736 "[" << da << "," << db << "]" << std::endl);
737 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
738 }
739
740 return dfbp;
741 }
742
743 doublereal xi = (2.*x - (da + db))/(db - da);
744 doublereal xip = 2./(db - da);
745 doublereal d[2] = { 1., xi };
746 doublereal dp[2] = { 0., 1. };
747 doublereal val = vCoef[1];
748
749 for (unsigned i = 2; i < vCoef.size(); i++) {
750 doublereal Tx = 2.*xi*d[1 - i%2] - d[i%2];
751 doublereal Txp = 2.*d[1 - i%2] + 2.*xi*dp[1 - i%2] - dp[i%2];
752 val += vCoef[i]*Txp;
753 d[i%2] = Tx;
754 dp[i%2] = Txp;
755 }
756
757 return xip*val;
758 }
759
760 // ScalarFunction parsing functional object
761 struct ChebychevSFR: public ScalarFunctionRead {
762 virtual const BasicScalarFunction *
ReadChebychevSFR763 Read(DataManager* const pDM, MBDynParser& HP) const {
764 doublereal a = HP.GetReal();
765 doublereal b = HP.GetReal();
766 bool doNotExtrapolate(false);
767 if (HP.IsKeyWord("do" "not" "extrapolate")) {
768 doNotExtrapolate = true;
769 }
770 if (b <= a) {
771 silent_cerr("Upper interval bound "
772 "of Chebychev series b=" << b
773 << " must be larger than lower bound a=" << a
774 << " at line" << HP.GetLineData() << std::endl);
775 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
776
777 }
778 std::vector<doublereal> v;
779 while (HP.IsArg() && !HP.IsKeyWord("end")) {
780 int size = v.size();
781 v.resize(size+1);
782 v[size] = HP.GetReal();
783 }
784 unsigned order = v.size();
785 if (order == 0) {
786 silent_cerr("Need at least one Chebychev series coefficient "
787 << "at line" << HP.GetLineData() << std::endl);
788 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
789 }
790 return new ChebychevScalarFunction(v, a, b, doNotExtrapolate);
791 };
792 };
793
794 // SumScalarFunction
SumScalarFunction(const BasicScalarFunction * const b1,const BasicScalarFunction * const b2)795 SumScalarFunction::SumScalarFunction(
796 const BasicScalarFunction*const b1,
797 const BasicScalarFunction*const b2)
798 : a1(ptr_cast<const DifferentiableScalarFunction*const>(b1)),
799 a2(ptr_cast<const DifferentiableScalarFunction*const>(b2))
800 {
801 NO_OP;
802 }
803
~SumScalarFunction(void)804 SumScalarFunction::~SumScalarFunction(void)
805 {
806 NO_OP;
807 }
808
809 doublereal
operator ()(const doublereal x) const810 SumScalarFunction::operator()(const doublereal x) const
811 {
812 return a1->operator()(x) + a2->operator()(x);
813 }
814
815 doublereal
ComputeDiff(const doublereal x,const integer order) const816 SumScalarFunction::ComputeDiff(const doublereal x, const integer order) const
817 {
818 ASSERTMSGBREAK(order >=0, "Error in SumScalarFunction::ComputeDiff, order<0");
819 switch (order) {
820 case 0:
821 return operator()(x);
822
823 default:
824 return a1->ComputeDiff(x, order)
825 + a2->ComputeDiff(x, order);
826 }
827 }
828
829 // ScalarFunction parsing functional object
830 struct SumSFR: public ScalarFunctionRead {
831 virtual const BasicScalarFunction *
ReadSumSFR832 Read(DataManager* const pDM, MBDynParser& HP) const {
833 const BasicScalarFunction *const
834 f1(ParseScalarFunction(HP, pDM));
835 const BasicScalarFunction *const
836 f2(ParseScalarFunction(HP, pDM));
837 return new SumScalarFunction(f1,f2);
838 };
839 };
840
841 // SubScalarFunction
SubScalarFunction(const BasicScalarFunction * const b1,const BasicScalarFunction * const b2)842 SubScalarFunction::SubScalarFunction(
843 const BasicScalarFunction*const b1,
844 const BasicScalarFunction*const b2)
845 : a1(ptr_cast<const DifferentiableScalarFunction*const>(b1)),
846 a2(ptr_cast<const DifferentiableScalarFunction*const>(b2))
847 {
848 NO_OP;
849 }
850
~SubScalarFunction(void)851 SubScalarFunction::~SubScalarFunction(void)
852 {
853 NO_OP;
854 }
855
856 doublereal
operator ()(const doublereal x) const857 SubScalarFunction::operator()(const doublereal x) const
858 {
859 return a1->operator()(x) - a2->operator()(x);
860 }
861
862 doublereal
ComputeDiff(const doublereal x,const integer order) const863 SubScalarFunction::ComputeDiff(const doublereal x, const integer order) const
864 {
865 ASSERTMSGBREAK(order >= 0, "Error in SubScalarFunction::ComputeDiff, order<0");
866 switch (order) {
867 case 0:
868 return operator()(x);
869
870 default:
871 return a1->ComputeDiff(x, order)
872 - a2->ComputeDiff(x, order);
873 }
874 }
875
876 // ScalarFunction parsing functional object
877 struct SubSFR: public ScalarFunctionRead {
878 virtual const BasicScalarFunction *
ReadSubSFR879 Read(DataManager* const pDM, MBDynParser& HP) const {
880 const BasicScalarFunction *const
881 f1(ParseScalarFunction(HP, pDM));
882 const BasicScalarFunction *const
883 f2(ParseScalarFunction(HP, pDM));
884 return new SubScalarFunction(f1,f2);
885 };
886 };
887
888 // MulScalarFunction
MulScalarFunction(const BasicScalarFunction * const b1,const BasicScalarFunction * const b2)889 MulScalarFunction::MulScalarFunction(
890 const BasicScalarFunction*const b1,
891 const BasicScalarFunction*const b2)
892 : a1(ptr_cast<const DifferentiableScalarFunction*const>(b1)),
893 a2(ptr_cast<const DifferentiableScalarFunction*const>(b2))
894 {
895 NO_OP;
896 }
897
~MulScalarFunction()898 MulScalarFunction::~MulScalarFunction()
899 {
900 NO_OP;
901 }
902
903 doublereal
operator ()(const doublereal x) const904 MulScalarFunction::operator()(const doublereal x) const
905 {
906 return a1->operator()(x)*a2->operator()(x);
907 }
908
909 doublereal
ComputeDiff(const doublereal x,const integer order) const910 MulScalarFunction::ComputeDiff(const doublereal x, const integer order) const
911 {
912 ASSERTMSGBREAK(order >= 0, "Error in MulScalarFunction::ComputeDiff, order<0");
913 switch (order) {
914 case 0:
915 return this->operator()(x);
916
917 default:
918 return a1->ComputeDiff(x, order)*a2->operator()(x)
919 + a1->operator()(x)*a2->ComputeDiff(x, order);
920 }
921 }
922
923 // ScalarFunction parsing functional object
924 struct MulSFR: public ScalarFunctionRead {
925 virtual const BasicScalarFunction *
ReadMulSFR926 Read(DataManager* const pDM, MBDynParser& HP) const {
927 const BasicScalarFunction *const
928 f1(ParseScalarFunction(HP, pDM));
929 const BasicScalarFunction *const
930 f2(ParseScalarFunction(HP, pDM));
931 return new MulScalarFunction(f1, f2);
932 };
933 };
934
935 // DivScalarFunction
DivScalarFunction(const BasicScalarFunction * const b1,const BasicScalarFunction * const b2)936 DivScalarFunction::DivScalarFunction(
937 const BasicScalarFunction*const b1,
938 const BasicScalarFunction*const b2)
939 : a1(ptr_cast<const DifferentiableScalarFunction*const>(b1)),
940 a2(ptr_cast<const DifferentiableScalarFunction*const>(b2))
941 {
942 NO_OP;
943 }
944
~DivScalarFunction()945 DivScalarFunction::~DivScalarFunction()
946 {
947 NO_OP;
948 }
949
950 doublereal
operator ()(const doublereal x) const951 DivScalarFunction::operator()(const doublereal x) const
952 {
953 doublereal n, d;
954 d = a2->operator()(x);
955 if (d == 0) {
956 /* TODO: cleanup exception handling */
957 silent_cerr("DivScalarFunction: division by zero" << std::endl);
958 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
959 }
960 n = a1->operator()(x);
961 return n/d;
962 }
963
964 doublereal
ComputeDiff(const doublereal x,const integer order) const965 DivScalarFunction::ComputeDiff(const doublereal x, const integer order) const
966 {
967 doublereal d;
968 ASSERTMSGBREAK(order >= 0, "Error in DivScalarFunction::ComputeDiff, order<0");
969 switch (order) {
970 case 0:
971 return this->operator()(x);
972
973 default:
974 d = a2->operator()(x);
975 if (d == 0.) {
976 /* TODO: cleanup exception handling */
977 silent_cerr("DivScalarFunction: division by zero" << std::endl);
978 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
979 }
980
981 return a1->ComputeDiff(x, order)/d
982 - a1->operator()(x)/(d*d)*a2->ComputeDiff(x, order);
983 }
984 }
985
986 // ScalarFunction parsing functional object
987 struct DivSFR: public ScalarFunctionRead {
988 virtual const BasicScalarFunction *
ReadDivSFR989 Read(DataManager* const pDM, MBDynParser& HP) const {
990 const BasicScalarFunction *const
991 f1(ParseScalarFunction(HP, pDM));
992 const BasicScalarFunction *const
993 f2(ParseScalarFunction(HP, pDM));
994 return new DivScalarFunction(f1, f2);
995 };
996 };
997
998 //---------------------------------------
999
1000 typedef std::map<std::string, const ScalarFunctionRead *, ltstrcase> SFReadType;
1001 static SFReadType SFRead;
1002
1003 struct SFWordSetType : public HighParser::WordSet {
IsWordSFWordSetType1004 bool IsWord(const std::string& s) const {
1005 return SFRead.find(std::string(s)) != SFRead.end();
1006 };
1007 };
1008 static SFWordSetType SFWordSet;
1009
1010 const BasicScalarFunction *const
ParseScalarFunction(MBDynParser & HP,DataManager * const pDM)1011 ParseScalarFunction(MBDynParser& HP, DataManager* const pDM)
1012 {
1013 std::string func_name(HP.GetStringWithDelims());
1014
1015 const BasicScalarFunction *sf = HP.GetScalarFunction(func_name);
1016 if (sf == 0) {
1017 const char *s = HP.IsWord(SFWordSet);
1018 if (s == 0) {
1019 s = "const";
1020 }
1021
1022 SFReadType::iterator func = SFRead.find(std::string(s));
1023 if (func == SFRead.end()) {
1024 silent_cerr("unknown scalar function type \"" << s << "\" "
1025 "for function \"" << func_name << "\" "
1026 "at line " << HP.GetLineData() << std::endl);
1027 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
1028 }
1029
1030 try {
1031 sf = func->second->Read(pDM, HP);
1032 } catch (...) {
1033 silent_cerr("Unable to parse "
1034 "ScalarFunction(\"" << func_name << "\") "
1035 "at line " << HP.GetLineData() << std::endl);
1036 throw;
1037 }
1038 if (!HP.SetScalarFunction(func_name, sf)) {
1039 silent_cerr("scalar function \"" << func_name << "\" "
1040 "already defined at line " << HP.GetLineData() << std::endl);
1041 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
1042 }
1043
1044 } else if (HP.IsWord(SFWordSet)) {
1045 silent_cerr("Error: redefinition of "
1046 "\"" << func_name << "\" scalar function "
1047 "at line " << HP.GetLineData() << std::endl);
1048 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
1049 }
1050
1051 return sf;
1052 }
1053
1054 //---------------------------------------
1055
1056 /* Scalar FunctionDrive - begin */
1057
1058 class ScalarFunctionDriveCaller : public DriveCaller {
1059 private:
1060 const BasicScalarFunction *const sc;
1061
1062 public:
ScalarFunctionDriveCaller(const DriveHandler * pDH,const BasicScalarFunction * const f)1063 ScalarFunctionDriveCaller(const DriveHandler* pDH,
1064 const BasicScalarFunction *const f)
1065 : DriveCaller(pDH), sc(f)
1066 { NO_OP; };
1067
~ScalarFunctionDriveCaller(void)1068 virtual ~ScalarFunctionDriveCaller(void) { NO_OP; };
1069
1070 /* Copia */
1071 virtual DriveCaller*
pCopy(void) const1072 pCopy(void) const {
1073 DriveCaller* pDC = NULL;
1074 SAFENEWWITHCONSTRUCTOR(pDC, ScalarFunctionDriveCaller,
1075 ScalarFunctionDriveCaller(pDrvHdl, sc));
1076 return pDC;
1077 };
1078
1079 /* Scrive il contributo del DriveCaller al file di restart */
Restart(std::ostream & out) const1080 virtual std::ostream& Restart(std::ostream& out) const {
1081 silent_cerr("ScalarFunctionDriveCaller: Restart not implemented"
1082 << std::endl);
1083 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
1084 };
1085
dGet(const doublereal & dVar) const1086 virtual doublereal dGet(const doublereal& dVar) const {
1087 return (*sc)(dVar);
1088 };
1089
bIsDifferentiable(void) const1090 virtual bool bIsDifferentiable(void) const {
1091 return (ptr_cast<const DifferentiableScalarFunction*const>(sc) != 0);
1092 };
1093
dGetP(const doublereal & dVar) const1094 virtual doublereal dGetP(const doublereal& dVar) const {
1095 return ptr_cast<const DifferentiableScalarFunction *const>(sc)->ComputeDiff(dVar);
1096 };
1097 };
1098
1099 struct ScalarFunctionDCR : public DriveCallerRead {
1100 DriveCaller *
1101 Read(const DataManager* pDM, MBDynParser& HP, bool bDeferred);
1102 };
1103
1104 DriveCaller *
Read(const DataManager * pDM,MBDynParser & HP,bool bDeferred)1105 ScalarFunctionDCR::Read(const DataManager* pDM, MBDynParser& HP, bool bDeferred)
1106 {
1107 NeedDM(pDM, HP, bDeferred, "scalar function");
1108
1109 /* driver legato alle scalar function */
1110 if (pDM == 0) {
1111 silent_cerr("sorry, since the driver is not owned by a DataManager" << std::endl
1112 << "no driver dependent drivers are allowed;" << std::endl
1113 << "aborting..." << std::endl);
1114 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
1115 }
1116
1117 const DriveHandler* pDrvHdl = pDM->pGetDrvHdl();
1118
1119 const BasicScalarFunction *const sc = ParseScalarFunction(HP, (DataManager *const)pDM);
1120
1121 DriveCaller *pDC = 0;
1122
1123 /* allocazione e creazione */
1124 SAFENEWWITHCONSTRUCTOR(pDC,
1125 ScalarFunctionDriveCaller,
1126 ScalarFunctionDriveCaller(pDrvHdl, sc));
1127
1128 return pDC;
1129 }
1130
1131 // ScalarFunction constitutive laws
1132
1133 // ScalarFunction elastic isotropic constitutive law
1134
1135 template <class T, class Tder>
1136 class ScalarFunctionIsotropicCL
1137 : public ConstitutiveLaw<T, Tder> {
1138 private:
1139 const DifferentiableScalarFunction * pSF;
1140 int n;
1141
1142 public:
ScalarFunctionIsotropicCL(const DifferentiableScalarFunction * psf)1143 ScalarFunctionIsotropicCL(const DifferentiableScalarFunction * psf)
1144 : pSF(psf) {
1145 if (typeid(T) == typeid(Vec3)) {
1146 n = 3;
1147
1148 } else if (typeid(T) == typeid(Vec6)) {
1149 n = 6;
1150
1151 } else {
1152 silent_cerr("ScalarFunctionIsotropicCL<" << typeid(T).name() << ", " << typeid(Tder).name() << "> not implemented" << std::endl);
1153 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
1154 }
1155 };
1156
~ScalarFunctionIsotropicCL(void)1157 virtual ~ScalarFunctionIsotropicCL(void) {
1158 NO_OP;
1159 };
1160
GetConstLawType(void) const1161 ConstLawType::Type GetConstLawType(void) const {
1162 return ConstLawType::ELASTIC;
1163 };
1164
pCopy(void) const1165 virtual ConstitutiveLaw<T, Tder>* pCopy(void) const {
1166 ConstitutiveLaw<T, Tder>* pCL = NULL;
1167
1168 typedef ScalarFunctionIsotropicCL<T, Tder> cl;
1169 SAFENEWWITHCONSTRUCTOR(pCL, cl, cl(pSF));
1170 return pCL;
1171 };
1172
Restart(std::ostream & out) const1173 virtual std::ostream& Restart(std::ostream& out) const {
1174 return out << "# not implemented!";
1175 };
1176
Update(const T & Eps,const T &=0.)1177 virtual void Update(const T& Eps, const T& /* EpsPrime */ = 0.) {
1178 ConstitutiveLaw<T, Tder>::Epsilon = Eps;
1179 for (int i = 1; i <= n; i++) {
1180 #if defined(MBDYN_X_WORKAROUND_GCC_3_2) || defined(MBDYN_X_WORKAROUND_GCC_3_3)
1181 ConstitutiveLaw<T, Tder>::F.Put(i, (*pSF)(Eps(i)));
1182 ConstitutiveLaw<T, Tder>::FDE.Put(i, i, pSF->ComputeDiff(Eps(i)));
1183 #else // !MBDYN_X_WORKAROUND_GCC_3_2 && ! MBDYN_X_WORKAROUND_GCC_3_3
1184 ConstitutiveLaw<T, Tder>::F(i) = (*pSF)(Eps(i));
1185 ConstitutiveLaw<T, Tder>::FDE(i, i) = pSF->ComputeDiff(Eps(i));
1186 #endif // !MBDYN_X_WORKAROUND_GCC_3_3 && ! MBDYN_X_WORKAROUND_GCC_3_3
1187 }
1188 };
1189 };
1190
1191 template <>
1192 class ScalarFunctionIsotropicCL<doublereal, doublereal>
1193 : public ConstitutiveLaw<doublereal, doublereal> {
1194 private:
1195 const DifferentiableScalarFunction *pSF;
1196
1197 public:
ScalarFunctionIsotropicCL(const DifferentiableScalarFunction * psf)1198 ScalarFunctionIsotropicCL(const DifferentiableScalarFunction * psf)
1199 : pSF(psf) {
1200 NO_OP;
1201 };
1202
~ScalarFunctionIsotropicCL(void)1203 virtual ~ScalarFunctionIsotropicCL(void) {
1204 NO_OP;
1205 };
1206
GetConstLawType(void) const1207 ConstLawType::Type GetConstLawType(void) const {
1208 return ConstLawType::ELASTIC;
1209 };
1210
pCopy(void) const1211 virtual ConstitutiveLaw<doublereal, doublereal>* pCopy(void) const {
1212 ConstitutiveLaw<doublereal, doublereal>* pCL = NULL;
1213
1214 typedef ScalarFunctionIsotropicCL<doublereal, doublereal> cl;
1215 SAFENEWWITHCONSTRUCTOR(pCL, cl, cl(pSF));
1216 return pCL;
1217 };
1218
Restart(std::ostream & out) const1219 virtual std::ostream& Restart(std::ostream& out) const {
1220 return out << "# not implemented!";
1221 };
1222
Update(const doublereal & Eps,const doublereal &=0.)1223 virtual void Update(const doublereal& Eps, const doublereal& /* EpsPrime */ = 0.) {
1224 ConstitutiveLaw<doublereal, doublereal>::Epsilon = Eps;
1225 ConstitutiveLaw<doublereal, doublereal>::F = (*pSF)(Eps);
1226 ConstitutiveLaw<doublereal, doublereal>::FDE = pSF->ComputeDiff(Eps);
1227 };
1228 };
1229
1230 /* specific functional object(s) */
1231 template <class T, class Tder>
1232 struct ScalarFunctionIsotropicCLR : public ConstitutiveLawRead<T, Tder> {
1233 virtual ConstitutiveLaw<T, Tder> *
1234 Read(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType);
1235 };
1236
1237 template <class T, class Tder>
1238 ConstitutiveLaw<T, Tder> *
Read(const DataManager * pDM,MBDynParser & HP,ConstLawType::Type & CLType)1239 ScalarFunctionIsotropicCLR<T, Tder>::Read(const DataManager* pDM,
1240 MBDynParser& HP,
1241 ConstLawType::Type& CLType)
1242 {
1243 ConstitutiveLaw<T, Tder>* pCL = 0;
1244
1245 CLType = ConstLawType::ELASTIC;
1246
1247 int n = 0;
1248 if (typeid(T) == typeid(doublereal)) {
1249 n = 1;
1250
1251 } else if (typeid(T) == typeid(Vec3)) {
1252 n = 3;
1253
1254 } else if (typeid(T) == typeid(Vec6)) {
1255 n = 6;
1256
1257 } else {
1258 silent_cerr("ScalarFunctionIsotropicCL"
1259 "<" << typeid(T).name() << ", " << typeid(Tder).name() << "> "
1260 "not implemented" << std::endl);
1261 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
1262 }
1263
1264 const BasicScalarFunction *pSF = ParseScalarFunction(HP, (DataManager *const)pDM);
1265 const DifferentiableScalarFunction *psf = dynamic_cast<const DifferentiableScalarFunction *>(pSF);
1266 if (psf == 0) {
1267 silent_cerr("ScalarFunction must be differentiable "
1268 "at line " << HP.GetLineData() << std::endl);
1269 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
1270 }
1271
1272 typedef ScalarFunctionIsotropicCL<T, Tder> L;
1273 SAFENEWWITHCONSTRUCTOR(pCL, L, L(psf));
1274
1275 return pCL;
1276 }
1277
1278 // ScalarFunction elastic orthotropic constitutive law
1279
1280 template <class T, class Tder>
1281 class ScalarFunctionOrthotropicCL
1282 : public ConstitutiveLaw<T, Tder> {
1283 private:
1284 std::vector<const DifferentiableScalarFunction *> SF;
1285 unsigned n;
1286
1287 public:
ScalarFunctionOrthotropicCL(const std::vector<const DifferentiableScalarFunction * > & sf)1288 ScalarFunctionOrthotropicCL(const std::vector<const DifferentiableScalarFunction *>& sf)
1289 {
1290 if (typeid(T) == typeid(Vec3)) {
1291 n = 3;
1292
1293 } else if (typeid(T) == typeid(Vec6)) {
1294 n = 6;
1295
1296 } else {
1297 silent_cerr("ScalarFunctionOrthotropicCL<" << typeid(T).name() << ", " << typeid(Tder).name() << "> not implemented" << std::endl);
1298 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
1299 }
1300
1301 ASSERT(sf.size() == n);
1302 SF.resize(n);
1303 for (unsigned i = 0; i < n; i++) {
1304 SF[i] = sf[i];
1305 }
1306 };
1307
~ScalarFunctionOrthotropicCL(void)1308 virtual ~ScalarFunctionOrthotropicCL(void) {
1309 NO_OP;
1310 };
1311
GetConstLawType(void) const1312 ConstLawType::Type GetConstLawType(void) const {
1313 return ConstLawType::ELASTIC;
1314 };
1315
pCopy(void) const1316 virtual ConstitutiveLaw<T, Tder>* pCopy(void) const {
1317 ConstitutiveLaw<T, Tder>* pCL = NULL;
1318
1319 typedef ScalarFunctionOrthotropicCL<T, Tder> cl;
1320 SAFENEWWITHCONSTRUCTOR(pCL, cl, cl(SF));
1321 return pCL;
1322 };
1323
Restart(std::ostream & out) const1324 virtual std::ostream& Restart(std::ostream& out) const {
1325 return out << "# not implemented!";
1326 };
1327
Update(const T & Eps,const T &=0.)1328 virtual void Update(const T& Eps, const T& /* EpsPrime */ = 0.) {
1329 ConstitutiveLaw<T, Tder>::Epsilon = Eps;
1330 for (unsigned i = 1; i <= n; i++) {
1331 /* skip null scalar functions */
1332 if (SF[i - 1] == 0) {
1333 continue;
1334 }
1335
1336 #if defined(MBDYN_X_WORKAROUND_GCC_3_2) || defined(MBDYN_X_WORKAROUND_GCC_3_3)
1337 ConstitutiveLaw<T, Tder>::F.Put(i, (*SF[i - 1])(Eps(i)));
1338 ConstitutiveLaw<T, Tder>::FDE.Put(i, i, SF[i - 1]->ComputeDiff(Eps(i)));
1339 #else // !MBDYN_X_WORKAROUND_GCC_3_2 && ! MBDYN_X_WORKAROUND_GCC_3_3
1340 ConstitutiveLaw<T, Tder>::F(i) = (*SF[i - 1])(Eps(i));
1341 ConstitutiveLaw<T, Tder>::FDE(i, i) = SF[i - 1]->ComputeDiff(Eps(i));
1342 #endif // !MBDYN_X_WORKAROUND_GCC_3_2 && ! MBDYN_X_WORKAROUND_GCC_3_3
1343 }
1344 };
1345 };
1346
1347 /* specific functional object(s) */
1348 template <class T, class Tder>
1349 struct ScalarFunctionOrthotropicCLR : public ConstitutiveLawRead<T, Tder> {
1350 virtual ConstitutiveLaw<T, Tder> *
1351 Read(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType);
1352 };
1353
1354 template <class T, class Tder>
1355 ConstitutiveLaw<T, Tder> *
Read(const DataManager * pDM,MBDynParser & HP,ConstLawType::Type & CLType)1356 ScalarFunctionOrthotropicCLR<T, Tder>::Read(const DataManager* pDM,
1357 MBDynParser& HP,
1358 ConstLawType::Type& CLType)
1359 {
1360 ConstitutiveLaw<T, Tder>* pCL = 0;
1361
1362 CLType = ConstLawType::ELASTIC;
1363
1364 int n = 1;
1365 if (typeid(T) == typeid(Vec3)) {
1366 n = 3;
1367
1368 } else if (typeid(T) == typeid(Vec6)) {
1369 n = 6;
1370
1371 } else if (typeid(T) != typeid(doublereal)) {
1372 silent_cerr("ScalarFunctionOrthotropicCL<" << typeid(T).name() << ", " << typeid(Tder).name() << "> not implemented" << std::endl);
1373 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
1374 }
1375
1376 std::vector<const DifferentiableScalarFunction *> SF(n);
1377 for (int i = 0; i < n; i++) {
1378 if (HP.IsKeyWord("null")) {
1379 SF[i] = 0;
1380
1381 } else {
1382 const BasicScalarFunction *pSF = ParseScalarFunction(HP, (DataManager *const)pDM);
1383 SF[i] = dynamic_cast<const DifferentiableScalarFunction *>(pSF);
1384 if (SF[i] == 0) {
1385 silent_cerr("ScalarFunction #" << i + 1 << " must be differentiable "
1386 "at line " << HP.GetLineData() << std::endl);
1387 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
1388 }
1389 }
1390 }
1391
1392 typedef ScalarFunctionOrthotropicCL<T, Tder> L;
1393 SAFENEWWITHCONSTRUCTOR(pCL, L, L(SF));
1394
1395 return pCL;
1396 }
1397
1398 bool
SetSF(const std::string & s,const ScalarFunctionRead * rf)1399 SetSF(const std::string &s, const ScalarFunctionRead *rf)
1400 {
1401 pedantic_cout("registering scalar function \"" << s << "\""
1402 << std::endl );
1403 return SFRead.insert(SFReadType::value_type(s, rf)).second;
1404 }
1405
1406 static unsigned done = 0;
1407
1408 void
InitSF(void)1409 InitSF(void)
1410 {
1411 if (::done++ > 0) {
1412 return;
1413 }
1414
1415 bool b;
1416
1417 b = SetSF("const", new ConstSFR);
1418 ASSERT(b);
1419 b = SetSF("linear", new LinearSFR);
1420 ASSERT(b);
1421 b = SetSF("pow", new PowSFR);
1422 ASSERT(b);
1423 b = SetSF("log", new LogSFR);
1424 ASSERT(b);
1425 b = SetSF("exp", new ExpSFR);
1426 ASSERT(b);
1427 b = SetSF("sum", new SumSFR);
1428 ASSERT(b);
1429 b = SetSF("sub", new SubSFR);
1430 ASSERT(b);
1431 b = SetSF("mul", new MulSFR);
1432 ASSERT(b);
1433 b = SetSF("div", new DivSFR);
1434 ASSERT(b);
1435 b = SetSF("cubic" "spline", new CubicSplineSFR);
1436 ASSERT(b);
1437 b = SetSF("multilinear", new MultiLinearSFR);
1438 ASSERT(b);
1439 b = SetSF("chebychev", new ChebychevSFR);
1440 ASSERT(b);
1441
1442 /* this is about initializing the scalar function drive */
1443 ScalarFunctionDCR *rf = new ScalarFunctionDCR;
1444 if (!SetDriveCallerData("scalar" "function", rf)) {
1445 delete rf;
1446
1447 silent_cerr("unable to register scalar function drive caller"
1448 << std::endl);
1449
1450 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
1451 }
1452
1453 /* this is about initializing the scalar function constitutive law(s) */
1454 ConstitutiveLawRead<doublereal, doublereal> *rf1D
1455 = new ScalarFunctionIsotropicCLR<doublereal, doublereal>;
1456 if (!SetCL1D("scalar" "function" "elastic" "isotropic", rf1D)) {
1457 delete rf1D;
1458
1459 silent_cerr("unable to register scalar function isotropic 1D constitutive law"
1460 << std::endl);
1461
1462 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
1463 }
1464
1465 rf1D = new ScalarFunctionIsotropicCLR<doublereal, doublereal>;
1466 if (!SetCL1D("scalar" "function" "elastic" "orthotropic", rf1D)) {
1467 delete rf1D;
1468
1469 silent_cerr("unable to register scalar function orthotropic 1D constitutive law"
1470 << std::endl);
1471
1472 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
1473 }
1474
1475 rf1D = new ScalarFunctionIsotropicCLR<doublereal, doublereal>;
1476 if (!SetCL1D("scalar" "function" "elastic", rf1D)) {
1477 delete rf1D;
1478
1479 silent_cerr("unable to register scalar function 1D constitutive law"
1480 << std::endl);
1481
1482 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
1483 }
1484
1485 ConstitutiveLawRead<Vec3, Mat3x3> *rf3D = new ScalarFunctionIsotropicCLR<Vec3, Mat3x3>;
1486 if (!SetCL3D("scalar" "function" "elastic" "isotropic", rf3D)) {
1487 delete rf3D;
1488
1489 silent_cerr("unable to register scalar function isotropic 3D constitutive law"
1490 << std::endl);
1491
1492 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
1493 }
1494
1495 rf3D = new ScalarFunctionOrthotropicCLR<Vec3, Mat3x3>;
1496 if (!SetCL3D("scalar" "function" "elastic" "orthotropic", rf3D)) {
1497 delete rf3D;
1498
1499 silent_cerr("unable to register scalar function orthotropic 3D constitutive law"
1500 << std::endl);
1501
1502 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
1503 }
1504
1505 ConstitutiveLawRead<Vec6, Mat6x6> *rf6D = new ScalarFunctionIsotropicCLR<Vec6, Mat6x6>;
1506 if (!SetCL6D("scalar" "function" "elastic" "isotropic", rf6D)) {
1507 delete rf6D;
1508
1509 silent_cerr("unable to register scalar function isotropic 6D constitutive law"
1510 << std::endl);
1511
1512 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
1513 }
1514
1515 rf6D = new ScalarFunctionOrthotropicCLR<Vec6, Mat6x6>;
1516 if (!SetCL6D("scalar" "function" "elastic" "orthotropic", rf6D)) {
1517 delete rf6D;
1518
1519 silent_cerr("unable to register scalar function orthotropic 6D constitutive law"
1520 << std::endl);
1521
1522 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
1523 }
1524 }
1525
1526 void
DestroySF(void)1527 DestroySF(void)
1528 {
1529 if (::done == 0) {
1530 silent_cerr("DestroySF() called once too many" << std::endl);
1531 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
1532 }
1533
1534 if (--done > 0) {
1535 return;
1536 }
1537
1538 for (SFReadType::iterator i = SFRead.begin(); i != SFRead.end(); ++i) {
1539 delete i->second;
1540 }
1541 SFRead.clear();
1542 }
1543
1544