1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/base/constltp_impl.cc,v 1.55 2017/01/12 14:46:09 masarati Exp $ */
2 /*
3 * MBDyn (C) is a multibody analysis code.
4 * http://www.mbdyn.org
5 *
6 * Copyright (C) 1996-2017
7 *
8 * Pierangelo Masarati <masarati@aero.polimi.it>
9 * Paolo Mantegazza <mantegazza@aero.polimi.it>
10 *
11 * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
12 * via La Masa, 34 - 20156 Milano, Italy
13 * http://www.aero.polimi.it
14 *
15 * Changing this copyright notice is forbidden.
16 *
17 * This program is free software; you can redistribute it and/or modify
18 * it under the terms of the GNU General Public License as published by
19 * the Free Software Foundation (version 2 of the License).
20 *
21 *
22 * This program is distributed in the hope that it will be useful,
23 * but WITHOUT ANY WARRANTY; without even the implied warranty of
24 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 * GNU General Public License for more details.
26 *
27 * You should have received a copy of the GNU General Public License
28 * along with this program; if not, write to the Free Software
29 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
30 */
31
32 /* Legami costitutivi */
33
34 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
35
36 #include "myassert.h"
37 #include "mynewmem.h"
38 #include "drive_.h"
39 #include "constltp_impl.h"
40
41 #include "symcltp.h"
42 #ifdef USE_GINAC
43 #include "ginaccltp.h"
44 #endif // USE_GINAC
45 #include "shockabsorber.h"
46 #include "constltp_ann.h"
47
48 /* constitutive laws sponsored by Hutchinson CdR */
49 #include "constltp_nlp.h"
50 #include "constltp_nlsf.h"
51
52 /* used by invariant constitutive law... */
53 #include "vehj.h"
54
55 /* ... */
56 #include "tdclw.h"
57 #include "constltp_axw.h"
58
59 /* constitutive law containers */
60 typedef std::map<std::string, ConstitutiveLawRead<doublereal, doublereal> *, ltstrcase> CL1DFuncMapType;
61 typedef std::map<std::string, ConstitutiveLawRead<Vec3, Mat3x3> *, ltstrcase> CL3DFuncMapType;
62 typedef std::map<std::string, ConstitutiveLawRead<Vec6, Mat6x6> *, ltstrcase> CL6DFuncMapType;
63
64 static CL1DFuncMapType CL1DFuncMap;
65 static CL3DFuncMapType CL3DFuncMap;
66 static CL6DFuncMapType CL6DFuncMap;
67
68 /* constitutive law parsing checkers */
69 struct CL1DWordSetType : public HighParser::WordSet {
IsWordCL1DWordSetType70 bool IsWord(const std::string& s) const {
71 return CL1DFuncMap.find(std::string(s)) != CL1DFuncMap.end();
72 };
73 };
74
75 struct CL3DWordSetType : public HighParser::WordSet {
IsWordCL3DWordSetType76 bool IsWord(const std::string& s) const {
77 return CL3DFuncMap.find(std::string(s)) != CL3DFuncMap.end();
78 };
79 };
80
81 struct CL6DWordSetType : public HighParser::WordSet {
IsWordCL6DWordSetType82 bool IsWord(const std::string& s) const {
83 return CL6DFuncMap.find(std::string(s)) != CL6DFuncMap.end();
84 };
85 };
86
87 static CL1DWordSetType CL1DWordSet;
88 static CL3DWordSetType CL3DWordSet;
89 static CL6DWordSetType CL6DWordSet;
90
91 /* constitutive law registration functions: call to register one */
92 bool
SetCL1D(const char * name,ConstitutiveLawRead<doublereal,doublereal> * rf)93 SetCL1D(const char *name, ConstitutiveLawRead<doublereal, doublereal> *rf)
94 {
95 pedantic_cout("registering constitutive law 1D \"" << name << "\""
96 << std::endl );
97 return CL1DFuncMap.insert(CL1DFuncMapType::value_type(name, rf)).second;
98 }
99
100 bool
SetCL3D(const char * name,ConstitutiveLawRead<Vec3,Mat3x3> * rf)101 SetCL3D(const char *name, ConstitutiveLawRead<Vec3, Mat3x3> *rf)
102 {
103 pedantic_cout("registering constitutive law 3D \"" << name << "\""
104 << std::endl );
105 return CL3DFuncMap.insert(CL3DFuncMapType::value_type(name, rf)).second;
106 }
107
108 bool
SetCL6D(const char * name,ConstitutiveLawRead<Vec6,Mat6x6> * rf)109 SetCL6D(const char *name, ConstitutiveLawRead<Vec6, Mat6x6> *rf)
110 {
111 pedantic_cout("registering constitutive law 6D \"" << name << "\""
112 << std::endl );
113 return CL6DFuncMap.insert(CL6DFuncMapType::value_type(name, rf)).second;
114 }
115
116 /* function that reads a constitutive law */
117 ConstitutiveLaw<doublereal, doublereal> *
ReadCL1D(const DataManager * pDM,MBDynParser & HP,ConstLawType::Type & CLType)118 ReadCL1D(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType)
119 {
120 const char *s, *sOrig = HP.IsWord(CL1DWordSet);
121 if (sOrig == 0) {
122 /* default to linear elastic? */
123 s = "linear" "elastic";
124 sOrig = "";
125
126 } else {
127 s = sOrig;
128 }
129
130 CL1DFuncMapType::iterator func = CL1DFuncMap.find(std::string(s));
131 if (func == CL1DFuncMap.end()) {
132 silent_cerr("unknown constitutive law 1D type "
133 "\"" << sOrig << "\" "
134 "at line " << HP.GetLineData() << std::endl);
135 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
136 }
137
138 return func->second->Read(pDM, HP, CLType);
139 }
140
141 ConstitutiveLaw<Vec3, Mat3x3> *
ReadCL3D(const DataManager * pDM,MBDynParser & HP,ConstLawType::Type & CLType)142 ReadCL3D(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType)
143 {
144 const char *s, *sOrig = HP.IsWord(CL3DWordSet);
145 if (sOrig == 0) {
146 #if 0
147 s = "linear" "elastic";
148 sOrig = "";
149 #else
150 silent_cerr("unknown constitutive law 3D type "
151 "at line " << HP.GetLineData() << std::endl);
152 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
153 #endif
154
155 } else {
156 s = sOrig;
157 }
158
159 CL3DFuncMapType::iterator func = CL3DFuncMap.find(std::string(s));
160 if (func == CL3DFuncMap.end()) {
161 silent_cerr("unknown constitutive law 3D type "
162 "\"" << sOrig << "\" "
163 "at line " << HP.GetLineData() << std::endl);
164 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
165 }
166
167 return func->second->Read(pDM, HP, CLType);
168 }
169
170 ConstitutiveLaw<Vec6, Mat6x6> *
ReadCL6D(const DataManager * pDM,MBDynParser & HP,ConstLawType::Type & CLType)171 ReadCL6D(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType)
172 {
173 const char *s, *sOrig = HP.IsWord(CL6DWordSet);
174 if (sOrig == 0) {
175 #if 0
176 s = "linear" "elastic";
177 sOrig = "";
178 #else
179 silent_cerr("unknown constitutive law 6D type "
180 "at line " << HP.GetLineData() << std::endl);
181 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
182 #endif
183
184 } else {
185 s = sOrig;
186 }
187
188 CL6DFuncMapType::iterator func = CL6DFuncMap.find(std::string(s));
189 if (func == CL6DFuncMap.end()) {
190 silent_cerr("unknown constitutive law 6D type "
191 "\"" << sOrig << "\" "
192 "at line " << HP.GetLineData() << std::endl);
193 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
194 }
195
196 return func->second->Read(pDM, HP, CLType);
197 }
198
199 /* specific functional object(s) */
200 struct CLArray1DR : public ConstitutiveLawRead<doublereal, doublereal> {
201 virtual ConstitutiveLaw<doublereal, doublereal> *
ReadCLArray1DR202 Read(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType) {
203 ConstitutiveLaw<doublereal, doublereal>* pCL = 0;
204
205 unsigned n = HP.GetInt();
206 if (n <= 0) {
207 silent_cerr("constitutive law array 1D: invalid constitutive law number " << n
208 << " at line " << HP.GetLineData() << std::endl);
209 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
210 }
211
212 if (n == 1) {
213 return ReadCL1D(pDM, HP, CLType);
214 }
215
216 std::vector<ConstitutiveLaw<doublereal, doublereal> *> clv(n);
217 for (unsigned i = 0; i < n; i++) {
218 clv[i] = ReadCL1D(pDM, HP, CLType);
219 }
220
221 typedef ConstitutiveLawArray<doublereal, doublereal> L;
222 SAFENEWWITHCONSTRUCTOR(pCL, L, L(clv));
223
224 CLType = pCL->GetConstLawType();
225
226 return pCL;
227 };
228 };
229
230 struct CLArray3DR : public ConstitutiveLawRead<Vec3, Mat3x3> {
231 virtual ConstitutiveLaw<Vec3, Mat3x3> *
ReadCLArray3DR232 Read(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType) {
233 ConstitutiveLaw<Vec3, Mat3x3>* pCL = 0;
234
235 unsigned n = HP.GetInt();
236 if (n <= 0) {
237 silent_cerr("constitutive law array 1D: invalid constitutive law number " << n
238 << " at line " << HP.GetLineData() << std::endl);
239 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
240 }
241
242 if (n == 1) {
243 return ReadCL3D(pDM, HP, CLType);
244 }
245
246 std::vector<ConstitutiveLaw<Vec3, Mat3x3> *> clv(n);
247 for (unsigned i = 0; i < n; i++) {
248 clv[i] = ReadCL3D(pDM, HP, CLType);
249 }
250
251 typedef ConstitutiveLawArray<Vec3, Mat3x3> L;
252 SAFENEWWITHCONSTRUCTOR(pCL, L, L(clv));
253
254 CLType = pCL->GetConstLawType();
255
256 return pCL;
257 };
258 };
259
260 struct CLArray6DR : public ConstitutiveLawRead<Vec6, Mat6x6> {
261 virtual ConstitutiveLaw<Vec6, Mat6x6> *
ReadCLArray6DR262 Read(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType) {
263 ConstitutiveLaw<Vec6, Mat6x6>* pCL = 0;
264
265 unsigned n = HP.GetInt();
266 if (n <= 0) {
267 silent_cerr("constitutive law array 1D: invalid constitutive law number " << n
268 << " at line " << HP.GetLineData() << std::endl);
269 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
270 }
271
272 if (n == 1) {
273 return ReadCL6D(pDM, HP, CLType);
274 }
275
276 std::vector<ConstitutiveLaw<Vec6, Mat6x6> *> clv(n);
277 for (unsigned i = 0; i < n; i++) {
278 clv[i] = ReadCL6D(pDM, HP, CLType);
279 }
280
281 typedef ConstitutiveLawArray<Vec6, Mat6x6> L;
282 SAFENEWWITHCONSTRUCTOR(pCL, L, L(clv));
283
284 CLType = pCL->GetConstLawType();
285
286 return pCL;
287 };
288 };
289
290 template <class T, class Tder>
291 struct LinearElasticCLR : public ConstitutiveLawRead<T, Tder> {
292 virtual ConstitutiveLaw<T, Tder> *
ReadLinearElasticCLR293 Read(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType) {
294 ConstitutiveLaw<T, Tder>* pCL = 0;
295
296 CLType = ConstLawType::ELASTIC;
297
298 doublereal dS = HP.GetReal();
299 DEBUGCOUT("Linear Elastic Isotropic Constitutive Law, stiffness = "
300 << dS << std::endl);
301
302 if (dS <= 0.) {
303 silent_cerr("warning, null or negative stiffness at line "
304 << HP.GetLineData() << std::endl);
305 }
306
307 /* Prestress and prestrain */
308 T PreStress(mb_zero<T>());
309 GetPreStress(HP, PreStress);
310 TplDriveCaller<T>* pTplDC = GetPreStrain<T>(pDM, HP);
311
312 typedef LinearElasticIsotropicConstitutiveLaw<T, Tder> L;
313 SAFENEWWITHCONSTRUCTOR(pCL, L, L(pTplDC, PreStress, dS));
314
315 return pCL;
316 };
317 };
318
319 template <class T, class Tder>
320 struct LinearElasticGenericCLR : public ConstitutiveLawRead<T, Tder> {
321 virtual ConstitutiveLaw<T, Tder> *
ReadLinearElasticGenericCLR322 Read(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType) {
323 ConstitutiveLaw<T, Tder>* pCL = 0;
324
325 CLType = ConstLawType::ELASTIC;
326
327 DEBUGCOUT("Linear Elastic Generic Constitutive Law" << std::endl);
328 Tder S(mb_zero<Tder>());
329 S = HP.Get(S);
330
331 /* Prestress and prestrain */
332 T PreStress(mb_zero<T>());
333 GetPreStress(HP, PreStress);
334 TplDriveCaller<T>* pTplDC = GetPreStrain<T>(pDM, HP);
335
336 typedef LinearElasticGenericConstitutiveLaw<T, Tder> L;
337 SAFENEWWITHCONSTRUCTOR(pCL, L, L(pTplDC, PreStress, S));
338
339 return pCL;
340 };
341 };
342
343 template <class T, class Tder>
344 struct LinearElasticGenericAxialTorsionCouplingCLR : public ConstitutiveLawRead<T, Tder> {
345 virtual ConstitutiveLaw<T, Tder> *
ReadLinearElasticGenericAxialTorsionCouplingCLR346 Read(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType) {
347 ConstitutiveLaw<T, Tder>* pCL = 0;
348
349 CLType = ConstLawType::ELASTIC;
350
351 DEBUGCOUT("Linear Elastic Generic Constitutive Law with Axial-Torsion Coupling" << std::endl);
352 Tder S(mb_zero<Tder>());
353 S = HP.Get(S);
354
355 /* coefficiente di accoppiamento */
356 doublereal dCoupl = HP.GetReal();
357 DEBUGCOUT("coupling coefficient: " << dCoupl << std::endl);
358
359 /* Prestress and prestrain */
360 T PreStress(mb_zero<T>());
361 GetPreStress(HP, PreStress);
362 TplDriveCaller<T>* pTplDC = GetPreStrain<T>(pDM, HP);
363
364 typedef LinearElasticGenericAxialTorsionCouplingConstitutiveLaw<T, Tder> L;
365 SAFENEWWITHCONSTRUCTOR(pCL, L, L(pTplDC, PreStress, S, dCoupl));
366
367 return pCL;
368 };
369 };
370
371 template <class T, class Tder>
372 struct LinearViscoElasticGenericAxialTorsionCouplingCLR : public ConstitutiveLawRead<T, Tder> {
373 virtual ConstitutiveLaw<T, Tder> *
ReadLinearViscoElasticGenericAxialTorsionCouplingCLR374 Read(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType) {
375 ConstitutiveLaw<T, Tder>* pCL = 0;
376
377 CLType = ConstLawType::VISCOELASTIC;
378
379 Tder S(mb_zero<Tder>());
380 S = HP.Get(S);
381
382 Tder SP(mb_zero<Tder>());
383 if (HP.IsKeyWord("proportional")) {
384 doublereal k = HP.GetReal();
385 SP = S*k;
386 } else {
387 SP = HP.Get(SP);
388 }
389
390 /* coefficiente di accoppiamento */
391 doublereal dCoupl = HP.GetReal();
392 DEBUGCOUT("coupling coefficient: " << dCoupl << std::endl);
393
394 /* Prestress and prestrain */
395 T PreStress(mb_zero<T>());
396 GetPreStress(HP, PreStress);
397 TplDriveCaller<T>* pTplDC = GetPreStrain<T>(pDM, HP);
398
399 typedef LinearViscoElasticGenericAxialTorsionCouplingConstitutiveLaw<T, Tder> L;
400 SAFENEWWITHCONSTRUCTOR(pCL, L, L(pTplDC, PreStress, S, SP, dCoupl));
401
402 return pCL;
403 };
404 };
405
406 struct InverseSquareElasticCLR : public ConstitutiveLawRead<doublereal, doublereal> {
407 virtual ConstitutiveLaw<doublereal, doublereal> *
ReadInverseSquareElasticCLR408 Read(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType) {
409 ConstitutiveLaw<doublereal, doublereal>* pCL = 0;
410
411 CLType = ConstLawType::ELASTIC;
412
413 DEBUGCOUT("Inverse Square Elastic Constitutive Law" << std::endl);
414 doublereal dA = HP.GetReal();
415 if (dA <= 0.) {
416 silent_cerr("warning, null or negative stiffness at line "
417 << HP.GetLineData() << std::endl);
418 }
419
420 doublereal dL0 = HP.GetReal();
421 if (dL0 <= 0.) {
422 silent_cerr("warning, null or negative reference length at line "
423 << HP.GetLineData() << std::endl);
424 }
425
426 /* Prestress and prestrain */
427 doublereal PreStress(mb_zero<doublereal>());
428 GetPreStress(HP, PreStress);
429 TplDriveCaller<doublereal>* pTplDC = GetPreStrain<doublereal>(pDM, HP);
430
431
432 SAFENEWWITHCONSTRUCTOR(pCL,
433 InverseSquareConstitutiveLaw,
434 InverseSquareConstitutiveLaw(pTplDC, PreStress, dA, dL0));
435
436 return pCL;
437 };
438 };
439
440 template <class T, class Tder>
441 struct LogElasticCLR : public ConstitutiveLawRead<T, Tder> {
442 virtual ConstitutiveLaw<T, Tder> *
ReadLogElasticCLR443 Read(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType) {
444 ConstitutiveLaw<T, Tder>* pCL = 0;
445
446 CLType = ConstLawType::ELASTIC;
447
448 DEBUGCOUT("Logarithmic Elastic Constitutive Law" << std::endl);
449 doublereal dS = HP.GetReal();
450 if (dS <= 0.) {
451 silent_cerr("warning, null or negative stiffness at line "
452 << HP.GetLineData() << std::endl);
453 }
454
455 /* Prestress and prestrain */
456 T PreStress(mb_zero<T>());
457 GetPreStress(HP, PreStress);
458 TplDriveCaller<T>* pTplDC = GetPreStrain<T>(pDM, HP);
459
460 typedef LogConstitutiveLaw<T, Tder> L;
461 SAFENEWWITHCONSTRUCTOR(pCL, L, L(pTplDC, PreStress, dS));
462
463 return pCL;
464 };
465 };
466
467 template <class T, class Tder>
468 struct DoubleLinearElasticCLR : public ConstitutiveLawRead<T, Tder> {
469 virtual ConstitutiveLaw<T, Tder> *
ReadDoubleLinearElasticCLR470 Read(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType) {
471 ConstitutiveLaw<T, Tder>* pCL = 0;
472
473 CLType = ConstLawType::ELASTIC;
474
475 doublereal dS = HP.GetReal();
476 DEBUGCOUT("stiffness = " << dS << std::endl);
477
478 if (dS <= 0.) {
479 silent_cerr("warning, null or negative stiffness at line "
480 << HP.GetLineData() << std::endl);
481 }
482
483 doublereal dUpp = HP.GetReal();
484 if (dUpp <= 0.) {
485 silent_cerr("warning, null or negative upper limit strain at line "
486 << HP.GetLineData() << std::endl);
487 }
488
489 doublereal dLow = HP.GetReal();
490 if (dLow >= 0.) {
491 silent_cerr("warning, null or positive lower limit strain at line "
492 << HP.GetLineData() << std::endl);
493 }
494
495 doublereal dSecondS = HP.GetReal();
496 if (dSecondS <= 0.) {
497 silent_cerr("warning, null or negative second stiffness at line "
498 << HP.GetLineData() << std::endl);
499 }
500
501 /* Prestress and prestrain */
502 T PreStress(mb_zero<T>());
503 GetPreStress(HP, PreStress);
504 TplDriveCaller<T>* pTplDC = GetPreStrain<T>(pDM, HP);
505
506 typedef DoubleLinearElasticConstitutiveLaw<T, Tder> L;
507 SAFENEWWITHCONSTRUCTOR(pCL,
508 L,
509 L(pTplDC, PreStress, dS, dUpp, dLow, dSecondS));
510
511 return pCL;
512 };
513 };
514
515 template <class T, class Tder>
516 struct IsotropicHardeningCLR : public ConstitutiveLawRead<T, Tder> {
517 virtual ConstitutiveLaw<T, Tder> *
ReadIsotropicHardeningCLR518 Read(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType) {
519 ConstitutiveLaw<T, Tder>* pCL = 0;
520
521 CLType = ConstLawType::ELASTIC;
522
523 doublereal dS = HP.GetReal();
524 DEBUGCOUT("Stiffness = " << dS << std::endl);
525
526 if (dS <= 0.) {
527 silent_cerr("warning, null or negative stiffness at line "
528 << HP.GetLineData() << std::endl);
529 }
530
531 doublereal dE = HP.GetReal();
532 DEBUGCOUT("Reference strain = " << dE << std::endl);
533
534 if (dE <= 0.) {
535 silent_cerr("error, null or negative reference strain at line "
536 << HP.GetLineData() << std::endl);
537 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
538 }
539
540 doublereal dS0 = 0.;
541 if (HP.IsKeyWord("linear" "stiffness")) {
542 dS0 = HP.GetReal();
543 }
544
545 /* Prestress and prestrain */
546 T PreStress(mb_zero<T>());
547 GetPreStress(HP, PreStress);
548 TplDriveCaller<T>* pTplDC = GetPreStrain<T>(pDM, HP);
549
550 typedef IsotropicHardeningConstitutiveLaw<T, Tder> L;
551 SAFENEWWITHCONSTRUCTOR(pCL, L, L(pTplDC, PreStress, dS, dS0, dE));
552
553 return pCL;
554 };
555 };
556
557 template <class T, class Tder>
558 struct ContactElasticCLR : public ConstitutiveLawRead<T, Tder> {
559 virtual ConstitutiveLaw<T, Tder> *
ReadContactElasticCLR560 Read(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType) {
561 ConstitutiveLaw<T, Tder>* pCL = 0;
562
563 CLType = ConstLawType::ELASTIC;
564
565 doublereal dK = HP.GetReal();
566 DEBUGCOUT("Stiffness = " << dK << std::endl);
567
568 if (dK <= 0.) {
569 silent_cerr("warning, null or negative stiffness at line "
570 << HP.GetLineData() << std::endl);
571 }
572
573 doublereal dGamma = HP.GetReal();
574 DEBUGCOUT("Exponent = " << dGamma << std::endl);
575
576 if (dGamma < 1.) {
577 silent_cerr("error, exponent < 1. at line "
578 << HP.GetLineData() << std::endl);
579 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
580 }
581
582 /* Prestress and prestrain */
583 T PreStress(mb_zero<T>());
584 GetPreStress(HP, PreStress);
585 TplDriveCaller<T>* pTplDC = GetPreStrain<T>(pDM, HP);
586
587 typedef ContactConstitutiveLaw<T, Tder> L;
588 SAFENEWWITHCONSTRUCTOR(pCL, L, L(pTplDC, PreStress, dK, dGamma));
589
590 return pCL;
591 };
592 };
593
594 template <class T, class Tder>
595 struct SymbolicCLR : public ConstitutiveLawRead<T, Tder> {
596 virtual ConstitutiveLaw<T, Tder> *
ReadSymbolicCLR597 Read(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType) {
598 ConstitutiveLaw<T, Tder>* pCL = 0;
599
600 unsigned dim;
601 if (typeid(T) == typeid(doublereal)) {
602 dim = 1;
603
604 } else if (typeid(T) == typeid(Vec3)) {
605 dim = 3;
606
607 } else if (typeid(T) == typeid(Vec6)) {
608 dim = 6;
609
610 } else {
611 silent_cerr("Invalid dimensionality "
612 "for symbolic constitutive law "
613 "at line " << HP.GetLineData()
614 << std::endl);
615 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
616 }
617
618 std::vector<std::string> epsilon;
619 if (CLType & ConstLawType::ELASTIC) {
620 if (!HP.IsKeyWord("epsilon")) {
621 silent_cerr("keyword \"epsilon\" expected at line " << HP.GetLineData() << std::endl);
622 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
623 }
624
625 epsilon.resize(dim);
626
627 for (unsigned row = 0; row < dim; row++) {
628 const char *tmp = HP.GetStringWithDelims();
629
630 if (tmp == 0) {
631 silent_cerr("unable to get \"epsilon\" "
632 "symbol #" << row << " "
633 "at line " << HP.GetLineData() << std::endl);
634 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
635 }
636 epsilon[row] = tmp;
637 }
638 }
639
640 std::vector<std::string> epsilonPrime;
641 if (CLType & ConstLawType::VISCOUS) {
642 if (!HP.IsKeyWord("epsilon" "prime")) {
643 silent_cerr("keyword \"epsilon prime\" expected at line " << HP.GetLineData() << std::endl);
644 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
645 }
646
647 epsilonPrime.resize(dim);
648
649 for (unsigned row = 0; row < dim; row++) {
650 const char *tmp = HP.GetStringWithDelims();
651
652 if (tmp == 0) {
653 silent_cerr("unable to get \"epsilonPrime\" "
654 "symbol #" << row << " "
655 "at line " << HP.GetLineData() << std::endl);
656 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
657 }
658 epsilonPrime[row] = tmp;
659 }
660 }
661
662 if (!HP.IsKeyWord("expression")) {
663 silent_cerr("keyword \"expression\" expected at line " << HP.GetLineData() << std::endl);
664 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
665 }
666
667 std::vector<std::string> expression(dim);
668 for (unsigned row = 0; row < dim; row++) {
669 const char *tmp = HP.GetStringWithDelims();
670 if (tmp == 0) {
671 silent_cerr("unable to get \"expression\" "
672 "#" << row << " "
673 "at line " << HP.GetLineData()
674 << std::endl);
675 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
676 }
677 expression[row] = tmp;
678 }
679
680 /* Prestress and prestrain */
681 T PreStress(mb_zero<T>());
682 GetPreStress(HP, PreStress);
683 #ifdef USE_GINAC
684 TplDriveCaller<T>* pTplDC =
685 #endif /* ! USE_GINAC */
686 GetPreStrain<T>(pDM, HP);
687
688 switch (CLType) {
689 case ConstLawType::ELASTIC: {
690 #ifdef USE_GINAC
691 typedef GiNaCElasticConstitutiveLaw<T, Tder> L;
692 SAFENEWWITHCONSTRUCTOR(pCL, L,
693 L(pTplDC, PreStress,
694 epsilon,
695 expression));
696 #else /* ! USE_GINAC */
697 silent_cerr("symbolic constitutive law not supported "
698 "at line " << HP.GetLineData() << std::endl);
699 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
700 #endif /* ! USE_GINAC */
701 break;
702 }
703
704 case ConstLawType::VISCOUS: {
705 #ifdef USE_GINAC
706 typedef GiNaCViscousConstitutiveLaw<T, Tder> L;
707 SAFENEWWITHCONSTRUCTOR(pCL, L,
708 L(PreStress, epsilonPrime, expression));
709 #else /* ! USE_GINAC */
710 silent_cerr("symbolic constitutive law not supported "
711 "at line " << HP.GetLineData() << std::endl);
712 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
713 #endif /* ! USE_GINAC */
714 break;
715 }
716
717 case ConstLawType::VISCOELASTIC: {
718 #ifdef USE_GINAC
719 typedef GiNaCViscoElasticConstitutiveLaw<T, Tder> L;
720 SAFENEWWITHCONSTRUCTOR(pCL, L,
721 L(pTplDC, PreStress,
722 epsilon, epsilonPrime,
723 expression));
724 #else /* ! USE_GINAC */
725 silent_cerr("symbolic constitutive law not supported "
726 "at line " << HP.GetLineData() << std::endl);
727 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
728 #endif /* ! USE_GINAC */
729 break;
730 }
731
732 default:
733 ASSERT(0);
734 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
735 }
736
737 return pCL;
738 };
739 };
740
741 template <class T, class Tder>
742 struct SymbolicElasticCLR : public SymbolicCLR<T, Tder> {
743 virtual ConstitutiveLaw<T, Tder> *
ReadSymbolicElasticCLR744 Read(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType) {
745 CLType = ConstLawType::ELASTIC;
746 return SymbolicCLR<T, Tder>::Read(pDM, HP, CLType);
747 };
748 };
749
750 template <class T, class Tder>
751 struct SymbolicViscousCLR : public SymbolicCLR<T, Tder> {
752 virtual ConstitutiveLaw<T, Tder> *
ReadSymbolicViscousCLR753 Read(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType) {
754 CLType = ConstLawType::VISCOUS;
755 return SymbolicCLR<T, Tder>::Read(pDM, HP, CLType);
756 };
757 };
758
759 template <class T, class Tder>
760 struct SymbolicViscoElasticCLR : public SymbolicCLR<T, Tder> {
761 virtual ConstitutiveLaw<T, Tder> *
ReadSymbolicViscoElasticCLR762 Read(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType) {
763 CLType = ConstLawType::VISCOELASTIC;
764 return SymbolicCLR<T, Tder>::Read(pDM, HP, CLType);
765 };
766 };
767
768 template <class T, class Tder>
769 struct LinearViscousCLR : public ConstitutiveLawRead<T, Tder> {
770 virtual ConstitutiveLaw<T, Tder> *
ReadLinearViscousCLR771 Read(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType) {
772 ConstitutiveLaw<T, Tder>* pCL = 0;
773
774 CLType = ConstLawType::VISCOUS;
775
776 doublereal dSP = HP.GetReal();
777 DEBUGCOUT("stiffness prime = " << dSP << std::endl);
778
779 if (dSP <= 0.) {
780 silent_cerr("warning, null or negative stiffness prime at line "
781 << HP.GetLineData() << std::endl);
782 }
783
784 /* Prestress (no prestrain) */
785 T PreStress(mb_zero<T>());
786 GetPreStress(HP, PreStress);
787
788 typedef LinearViscousIsotropicConstitutiveLaw<T, Tder> L;
789 SAFENEWWITHCONSTRUCTOR(pCL, L, L(PreStress, dSP));
790
791 return pCL;
792 };
793 };
794
795 template <class T, class Tder>
796 struct LinearViscousGenericCLR : public ConstitutiveLawRead<T, Tder> {
797 virtual ConstitutiveLaw<T, Tder> *
ReadLinearViscousGenericCLR798 Read(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType) {
799 ConstitutiveLaw<T, Tder>* pCL = 0;
800
801 CLType = ConstLawType::VISCOUS;
802
803 Tder SP(mb_zero<Tder>());
804 SP = HP.Get(SP);
805
806 /* Prestress (no prestrain) */
807 T PreStress(mb_zero<T>());
808 GetPreStress(HP, PreStress);
809
810 typedef LinearViscousGenericConstitutiveLaw<T, Tder> L;
811 SAFENEWWITHCONSTRUCTOR(pCL, L, L(PreStress, SP));
812
813 return pCL;
814 };
815 };
816
817 template <class T, class Tder>
818 struct LinearViscoElasticCLR : public ConstitutiveLawRead<T, Tder> {
819 virtual ConstitutiveLaw<T, Tder> *
ReadLinearViscoElasticCLR820 Read(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType) {
821 ConstitutiveLaw<T, Tder>* pCL = 0;
822
823 CLType = ConstLawType::VISCOELASTIC;
824
825 doublereal dS = HP.GetReal();
826 DEBUGCOUT("Stiffness = " << dS << std::endl);
827
828 if (dS <= 0.) {
829 silent_cerr("warning, null or negative stiffness at line "
830 << HP.GetLineData() << std::endl);
831 }
832
833 doublereal dSP = 0.;
834 if (HP.IsKeyWord("proportional")) {
835 doublereal k = HP.GetReal();
836 dSP = k*dS;
837 } else {
838 dSP = HP.GetReal();
839 }
840 DEBUGCOUT("stiffness prime = " << dSP << std::endl);
841
842 if (dSP <= 0.) {
843 silent_cerr("warning, null or negative stiffness prime at line "
844 << HP.GetLineData() << std::endl);
845 }
846
847 /* Prestress and prestrain */
848 T PreStress(mb_zero<T>());
849 GetPreStress(HP, PreStress);
850 TplDriveCaller<T>* pTplDC = GetPreStrain<T>(pDM, HP);
851
852 #if 0 // TODO: implement a "null" constitutive law that does nothing
853 if (dS == 0. && dSP == 0.) {
854
855 } else
856 #endif
857 if (dS == 0.) {
858 silent_cerr("warning, null stiffness, "
859 "using linear viscous constitutive law instead"
860 << std::endl);
861
862 typedef LinearViscousIsotropicConstitutiveLaw<T, Tder> L;
863 SAFENEWWITHCONSTRUCTOR(pCL, L, L(PreStress, dSP));
864 if (pTplDC) {
865 delete pTplDC;
866 }
867
868 } else if (dSP == 0.) {
869 silent_cerr("warning, null stiffness prime, "
870 "using linear elastic constitutive law instead"
871 << std::endl);
872
873 typedef LinearElasticIsotropicConstitutiveLaw<T, Tder> L;
874 SAFENEWWITHCONSTRUCTOR(pCL, L, L(pTplDC, PreStress, dS));
875
876 } else {
877 typedef LinearViscoElasticIsotropicConstitutiveLaw<T, Tder> L;
878 SAFENEWWITHCONSTRUCTOR(pCL, L, L(pTplDC, PreStress, dS, dSP));
879 }
880
881 return pCL;
882 };
883 };
884
885 template <class T, class Tder>
886 struct LinearViscoElasticGenericCLR : public ConstitutiveLawRead<T, Tder> {
887 virtual ConstitutiveLaw<T, Tder> *
ReadLinearViscoElasticGenericCLR888 Read(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType) {
889 ConstitutiveLaw<T, Tder>* pCL = 0;
890
891 CLType = ConstLawType::VISCOELASTIC;
892
893 Tder S(mb_zero<Tder>());
894 S = HP.Get(S);
895
896 Tder SP(mb_zero<Tder>());
897 if (HP.IsKeyWord("proportional")) {
898 doublereal k = HP.GetReal();
899 SP = S*k;
900
901 } else {
902 SP = HP.Get(SP);
903 }
904
905 /* Prestress and prestrain */
906 T PreStress(mb_zero<T>());
907 GetPreStress(HP, PreStress);
908 TplDriveCaller<T>* pTplDC = GetPreStrain<T>(pDM, HP);
909
910 #if 0 // TODO: implement a "null" constitutive law that does nothing
911 if (IsNull(S) && IsNull(SP)) {
912
913 } else
914 #endif
915 if (IsNull(S)) {
916 silent_cerr("warning, null stiffness, "
917 "using linear viscous generic constitutive law instead"
918 << std::endl);
919
920 typedef LinearViscousGenericConstitutiveLaw<T, Tder> L;
921 SAFENEWWITHCONSTRUCTOR(pCL, L, L(PreStress, SP));
922 if (pTplDC) {
923 delete pTplDC;
924 }
925
926 } else if (IsNull(SP)) {
927 silent_cerr("warning, null stiffness prime, "
928 "using linear elastic generic constitutive law instead"
929 << std::endl);
930
931 typedef LinearElasticGenericConstitutiveLaw<T, Tder> L;
932 SAFENEWWITHCONSTRUCTOR(pCL, L, L(pTplDC, PreStress, S));
933
934 } else {
935 typedef LinearViscoElasticGenericConstitutiveLaw<T, Tder> L;
936 SAFENEWWITHCONSTRUCTOR(pCL, L, L(pTplDC, PreStress, S, SP));
937 }
938
939 return pCL;
940 };
941 };
942
943 template <class T, class Tder>
944 struct LTVViscoElasticGenericCLR : public ConstitutiveLawRead<T, Tder> {
945 virtual ConstitutiveLaw<T, Tder> *
ReadLTVViscoElasticGenericCLR946 Read(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType) {
947 ConstitutiveLaw<T, Tder>* pCL = 0;
948
949 CLType = ConstLawType::VISCOELASTIC;
950
951 Tder S(mb_zero<Tder>());
952 S = HP.Get(S);
953
954 DriveCaller *pdc = HP.GetDriveCaller();
955
956 Tder SP(mb_zero<Tder>());
957 if (HP.IsKeyWord("proportional")) {
958 doublereal k = HP.GetReal();
959 SP = S*k;
960
961 } else {
962 SP = HP.Get(SP);
963 }
964
965 DriveCaller *pdcp = HP.GetDriveCaller();
966
967 /* Prestress and prestrain */
968 T PreStress(mb_zero<T>());
969 GetPreStress(HP, PreStress);
970 TplDriveCaller<T>* pTplDC = GetPreStrain<T>(pDM, HP);
971
972 #if 0 // TODO: implement a "null" constitutive law that does nothing
973 if (IsNull(S) && IsNull(SP)) {
974
975 } else if (IsNull(S)) {
976 silent_cerr("warning, null stiffness, "
977 "using linear viscous generic constitutive law instead"
978 << std::endl);
979
980 SAFEDELETE(pdc);
981
982 typedef LTVViscousGenericConstitutiveLaw<T, Tder> L;
983 SAFENEWWITHCONSTRUCTOR(pCL, L, L(PreStress, SP, pdcp));
984 if (pTplDC) {
985 delete pTplDC;
986 }
987
988 } else
989 if (IsNull(SP)) {
990 silent_cerr("warning, null stiffness prime, "
991 "using linear elastic generic constitutive law instead"
992 << std::endl);
993
994 SAFEDELETE(pdcp);
995
996 typedef LTVElasticGenericConstitutiveLaw<T, Tder> L;
997 SAFENEWWITHCONSTRUCTOR(pCL, L, L(pTplDC, PreStress, S, pdc));
998
999 } else
1000 #endif
1001 {
1002 typedef LTVViscoElasticGenericConstitutiveLaw<T, Tder> L;
1003 SAFENEWWITHCONSTRUCTOR(pCL, L, L(pTplDC, PreStress, S, pdc, SP, pdcp));
1004 }
1005
1006 return pCL;
1007 };
1008 };
1009
1010 template <class T, class Tder>
1011 struct CubicElasticGenericCLR : public ConstitutiveLawRead<T, Tder> {
1012 virtual ConstitutiveLaw<T, Tder> *
ReadCubicElasticGenericCLR1013 Read(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType) {
1014 ConstitutiveLaw<T, Tder>* pCL = 0;
1015
1016 CLType = ConstLawType::ELASTIC;
1017
1018 T S1(mb_zero<T>());
1019 S1 = HP.Get(S1);
1020
1021 T S2(mb_zero<T>());
1022 S2 = HP.Get(S2);
1023
1024 T S3(mb_zero<T>());
1025 S3 = HP.Get(S3);
1026
1027 /* Prestress and prestrain */
1028 T PreStress(mb_zero<T>());
1029 GetPreStress(HP, PreStress);
1030 TplDriveCaller<T>* pTplDC = GetPreStrain<T>(pDM, HP);
1031
1032 typedef CubicElasticGenericConstitutiveLaw<T, Tder> L;
1033 SAFENEWWITHCONSTRUCTOR(pCL, L, L(pTplDC, PreStress, S1, S2, S3));
1034
1035 return pCL;
1036 };
1037 };
1038
1039 template <class T, class Tder>
1040 struct CubicViscoElasticGenericCLR : public ConstitutiveLawRead<T, Tder> {
1041 virtual ConstitutiveLaw<T, Tder> *
ReadCubicViscoElasticGenericCLR1042 Read(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType) {
1043 ConstitutiveLaw<T, Tder>* pCL = 0;
1044
1045 CLType = ConstLawType::VISCOELASTIC;
1046
1047 T S1(mb_zero<T>());
1048 S1 = HP.Get(S1);
1049
1050 T S2(mb_zero<T>());
1051 S2 = HP.Get(S2);
1052
1053 T S3(mb_zero<T>());
1054 S3 = HP.Get(S3);
1055
1056 Tder SP(mb_zero<Tder>());
1057 SP = HP.Get(SP);
1058
1059 /* Prestress and prestrain */
1060 T PreStress(mb_zero<T>());
1061 GetPreStress(HP, PreStress);
1062 TplDriveCaller<T>* pTplDC = GetPreStrain<T>(pDM, HP);
1063
1064 typedef CubicViscoElasticGenericConstitutiveLaw<T, Tder> L;
1065 SAFENEWWITHCONSTRUCTOR(pCL, L, L(pTplDC, PreStress, S1, S2, S3, SP));
1066
1067 return pCL;
1068 };
1069 };
1070
1071 template <class T, class Tder>
1072 struct DoubleLinearViscoElasticCLR : public ConstitutiveLawRead<T, Tder> {
1073 virtual ConstitutiveLaw<T, Tder> *
ReadDoubleLinearViscoElasticCLR1074 Read(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType) {
1075 ConstitutiveLaw<T, Tder>* pCL = 0;
1076
1077 CLType = ConstLawType::VISCOELASTIC;
1078
1079 doublereal dS = HP.GetReal();
1080 DEBUGCOUT("stiffness = " << dS << std::endl);
1081
1082 if (dS <= 0.) {
1083 silent_cerr("warning, null or negative stiffness at line "
1084 << HP.GetLineData() << std::endl);
1085 }
1086
1087 doublereal dUpp = HP.GetReal();
1088 if (dUpp <= 0.) {
1089 silent_cerr("warning, null or negative upper limit strain at line "
1090 << HP.GetLineData() << std::endl);
1091 }
1092
1093 doublereal dLow = HP.GetReal();
1094 if (dLow >= 0.) {
1095 silent_cerr("warning, null or positive lower limit strain at line "
1096 << HP.GetLineData() << std::endl);
1097 }
1098
1099 doublereal dSecondS = HP.GetReal();
1100 if (dSecondS <= 0.) {
1101 silent_cerr("warning, null or negative second stiffness at line "
1102 << HP.GetLineData() << std::endl);
1103 }
1104
1105 doublereal dSP = HP.GetReal();
1106 DEBUGCOUT("stiffness prime = " << dSP << std::endl);
1107
1108 if (dSP <= 0.) {
1109 silent_cerr("warning, null or negative stiffness prime at line "
1110 << HP.GetLineData() << std::endl);
1111 }
1112
1113 doublereal dSecondSP = dSP;
1114 if (HP.IsKeyWord("second" "damping")) {
1115 dSecondSP = HP.GetReal();
1116 DEBUGCOUT("second stiffness prime = " << dSecondSP << std::endl);
1117
1118 if (dSecondSP <= 0.) {
1119 silent_cerr("warning, null or negative second stiffness prime at line "
1120 << HP.GetLineData() << std::endl);
1121 }
1122 }
1123
1124 /* Prestress and prestrain */
1125 T PreStress(mb_zero<T>());
1126 GetPreStress(HP, PreStress);
1127 TplDriveCaller<T>* pTplDC = GetPreStrain<T>(pDM, HP);
1128
1129 typedef DoubleLinearViscoElasticConstitutiveLaw<T, Tder> L;
1130 SAFENEWWITHCONSTRUCTOR(pCL,
1131 L,
1132 L(pTplDC, PreStress,
1133 dS, dUpp, dLow, dSecondS, dSP, dSecondSP));
1134
1135 return pCL;
1136 };
1137 };
1138
1139 template <class T, class Tder>
1140 struct TurbulentViscoElasticCLR : public ConstitutiveLawRead<T, Tder> {
1141 virtual ConstitutiveLaw<T, Tder> *
ReadTurbulentViscoElasticCLR1142 Read(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType) {
1143 ConstitutiveLaw<T, Tder>* pCL = 0;
1144
1145 CLType = ConstLawType::VISCOELASTIC;
1146
1147 doublereal dS = HP.GetReal();
1148 DEBUGCOUT("Visco-Elastic Turbulent Rod Joint, stiffness = "
1149 << dS << std::endl);
1150
1151 if (dS <= 0.) {
1152 silent_cerr("warning, null or negative stiffness at line "
1153 << HP.GetLineData() << std::endl);
1154 }
1155
1156 doublereal dParabStiff = HP.GetReal();
1157 DEBUGCOUT("stiffness prime = " << dParabStiff << std::endl);
1158
1159 if (dParabStiff <= 0.) {
1160 silent_cerr("warning, null or negative derivative stiffness at line "
1161 << HP.GetLineData() << std::endl);
1162 }
1163
1164 doublereal dTreshold = 0.;
1165 if (HP.IsArg()) {
1166 dTreshold = HP.GetReal(dTreshold);
1167
1168 /*
1169 * Il legame costitutivo ha la forma seguente:
1170 * F = Kp*e + Kd*(de/dt)
1171 * con Kp costante e Kd dato dalla seguente legge:
1172 * Kd = cost2 per fabs(de/dt) < Treshold
1173 * Kd = 2*cost1*fabs(de/dt) per fabs(de/dt) > Treshold
1174 * se non viene inserito il valore di treshold, lo si
1175 * assume = 0. e quindi il legame e' sempre del secondo
1176 * tipo. Altrimenti, se non viene inserita la seconda
1177 * costante cost2, si assume che vi sia raccordo tra
1178 * i due tipi di legge, ovvero cost2 = cost1*Treshold
1179 * altrimenti e' possibile avere un comportamento,
1180 * che in prima approssimazione e' valido
1181 * per numerosi fluidi, in cui vi e' un salto tra i due
1182 * tipi di legge costitutiva. */
1183 }
1184
1185 doublereal dSP = dTreshold*dParabStiff;
1186 if (HP.IsArg()) {
1187 dSP = HP.GetReal(dSP);
1188 }
1189
1190 /* Prestress and prestrain */
1191 T PreStress(mb_zero<T>());
1192 GetPreStress(HP, PreStress);
1193 TplDriveCaller<T>* pTplDC = GetPreStrain<T>(pDM, HP);
1194
1195 typedef TurbulentViscoElasticConstitutiveLaw<T, Tder> L;
1196 SAFENEWWITHCONSTRUCTOR(pCL,
1197 L,
1198 L(pTplDC, PreStress,
1199 dS, dSP, dTreshold, dParabStiff));
1200
1201 return pCL;
1202 };
1203 };
1204
1205 template <class T, class Tder>
1206 struct LinearBiStopCLR : public ConstitutiveLawRead<T, Tder> {
1207 virtual ConstitutiveLaw<T, Tder> *
ReadLinearBiStopCLR1208 Read(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType) {
1209 ConstitutiveLaw<T, Tder>* pCL = 0;
1210
1211 typedef BiStopCLWrapper<T, Tder> L;
1212 typedef LinearElasticIsotropicConstitutiveLaw<T, Tder> LECL;
1213 typedef LinearViscoElasticIsotropicConstitutiveLaw<T, Tder> LVECL;
1214
1215 DEBUGCOUT("Linear Viscoelastic Bi Stop Constitutive Law" << std::endl);
1216 doublereal dS = HP.GetReal();
1217 if (dS <= 0.) {
1218 silent_cerr("warning, null or negative stiffness at line "
1219 << HP.GetLineData() << std::endl);
1220 }
1221
1222 doublereal dSp = 0.;
1223 if (CLType == ConstLawType::VISCOELASTIC) {
1224 dSp = HP.GetReal();
1225 if (dSp <= 0.) {
1226 silent_cerr("warning, null or negative stiffness prime at line "
1227 << HP.GetLineData() << std::endl);
1228 }
1229 }
1230
1231 bool s(false);
1232 if (HP.IsKeyWord("initial" "status")) {
1233 if (HP.IsKeyWord("active")) {
1234 s = true;
1235
1236 } else if (HP.IsKeyWord("inactive")) {
1237 s = false;
1238
1239 } else {
1240 s = HP.GetBool();
1241 }
1242 }
1243
1244 const DriveCaller *pA = HP.GetDriveCaller();
1245 const DriveCaller *pD = HP.GetDriveCaller();
1246
1247 /* Prestress and prestrain */
1248 T PreStress(mb_zero<T>());
1249 GetPreStress(HP, PreStress);
1250 TplDriveCaller<T>* pTplDC = GetPreStrain<T>(pDM, HP);
1251
1252 ConstitutiveLaw<T, Tder> *pWrappedCL = 0;
1253 if (CLType == ConstLawType::VISCOELASTIC && dSp != 0.) {
1254 SAFENEWWITHCONSTRUCTOR(pWrappedCL, LVECL, LVECL(pTplDC, PreStress, dS, dSp));
1255
1256 } else {
1257 SAFENEWWITHCONSTRUCTOR(pWrappedCL, LECL, LECL(pTplDC, PreStress, dS));
1258 }
1259
1260 SAFENEWWITHCONSTRUCTOR(pCL, L, L(pWrappedCL, s, pA, pD));
1261
1262 return pCL;
1263 };
1264 };
1265
1266 template <class T, class Tder>
1267 struct LinearElasticBiStopCLR : public LinearBiStopCLR<T, Tder> {
1268 virtual ConstitutiveLaw<T, Tder> *
ReadLinearElasticBiStopCLR1269 Read(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType) {
1270 CLType = ConstLawType::ELASTIC;
1271 return LinearBiStopCLR<T, Tder>::Read(pDM, HP, CLType);
1272 };
1273 };
1274
1275 template <class T, class Tder>
1276 struct LinearViscoElasticBiStopCLR : public LinearBiStopCLR<T, Tder> {
1277 virtual ConstitutiveLaw<T, Tder> *
ReadLinearViscoElasticBiStopCLR1278 Read(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType) {
1279 CLType = ConstLawType::VISCOELASTIC;
1280 return LinearBiStopCLR<T, Tder>::Read(pDM, HP, CLType);
1281 };
1282 };
1283
1284 static void
ReadBiStopBase(MBDynParser & HP,bool & bStatus,const DriveCaller * & pA,const DriveCaller * & pD)1285 ReadBiStopBase(MBDynParser& HP, bool& bStatus, const DriveCaller *& pA, const DriveCaller *& pD)
1286 {
1287 if (HP.IsKeyWord("initial" "status")) {
1288 if (HP.IsKeyWord("active")) {
1289 bStatus = true;
1290
1291 } else if (HP.IsKeyWord("inactive")) {
1292 bStatus = false;
1293
1294 } else {
1295 bStatus = HP.GetBool();
1296 }
1297 }
1298
1299 pA = HP.GetDriveCaller();
1300 pD = HP.GetDriveCaller();
1301 }
1302
1303 struct BiStopCLW1DR : public ConstitutiveLawRead<doublereal, doublereal> {
1304 virtual ConstitutiveLaw<doublereal, doublereal> *
ReadBiStopCLW1DR1305 Read(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType) {
1306 ConstitutiveLaw<doublereal, doublereal>* pCL = 0;
1307
1308 typedef BiStopCLWrapper<doublereal, doublereal> L;
1309
1310 const DriveCaller *pA = 0;
1311 const DriveCaller *pD = 0;
1312 bool bStatus(false);
1313 ReadBiStopBase(HP, bStatus, pA, pD);
1314
1315 ConstitutiveLaw<doublereal, doublereal> *pWrappedCL = ReadCL1D(pDM, HP, CLType);
1316 SAFENEWWITHCONSTRUCTOR(pCL, L, L(pWrappedCL, bStatus, pA, pD));
1317
1318 return pCL;
1319 };
1320 };
1321
1322 struct BiStopCLW3DR : public ConstitutiveLawRead<Vec3, Mat3x3> {
1323 virtual ConstitutiveLaw<Vec3, Mat3x3> *
ReadBiStopCLW3DR1324 Read(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType) {
1325 ConstitutiveLaw<Vec3, Mat3x3>* pCL = 0;
1326
1327 typedef BiStopCLWrapper<Vec3, Mat3x3> L;
1328
1329 const DriveCaller *pA = 0;
1330 const DriveCaller *pD = 0;
1331 bool bStatus(false);
1332 ReadBiStopBase(HP, bStatus, pA, pD);
1333
1334 ConstitutiveLaw<Vec3, Mat3x3> *pWrappedCL = ReadCL3D(pDM, HP, CLType);
1335 SAFENEWWITHCONSTRUCTOR(pCL, L, L(pWrappedCL, bStatus, pA, pD));
1336
1337 return pCL;
1338 };
1339 };
1340
1341 struct BiStopCLW6DR : public ConstitutiveLawRead<Vec6, Mat6x6> {
1342 virtual ConstitutiveLaw<Vec6, Mat6x6> *
ReadBiStopCLW6DR1343 Read(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType) {
1344 ConstitutiveLaw<Vec6, Mat6x6>* pCL = 0;
1345
1346 typedef BiStopCLWrapper<Vec6, Mat6x6> L;
1347
1348 const DriveCaller *pA = 0;
1349 const DriveCaller *pD = 0;
1350 bool bStatus(false);
1351 ReadBiStopBase(HP, bStatus, pA, pD);
1352
1353 ConstitutiveLaw<Vec6, Mat6x6> *pWrappedCL = ReadCL6D(pDM, HP, CLType);
1354 SAFENEWWITHCONSTRUCTOR(pCL, L, L(pWrappedCL, bStatus, pA, pD));
1355
1356 return pCL;
1357 };
1358 };
1359
1360 /*
1361 * Shock absorber per Stefy:
1362 *
1363 * ``Riprogettazione dell'ammortizzatore del carrello anteriore
1364 * di un velivolo di aviazione generale'',
1365 * S. Carlucci e S. Gualdi,
1366 * A.A. 1997-98
1367 */
1368 template <class T, class Tder>
1369 struct ShockAbsorberCLR : public ConstitutiveLawRead<T, Tder> {
1370 virtual ConstitutiveLaw<T, Tder> *
ReadShockAbsorberCLR1371 Read(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType) {
1372 ConstitutiveLaw<T, Tder>* pCL = 0;
1373
1374 CLType = ConstLawType::VISCOELASTIC;
1375
1376 TplDriveCaller<T>* pTplDC = GetPreStrain<T>(pDM, HP);
1377
1378 typedef ShockAbsorberConstitutiveLaw<T, Tder> L;
1379 SAFENEWWITHCONSTRUCTOR(pCL, L, L(pDM, pTplDC, HP));
1380
1381 return pCL;
1382 };
1383 };
1384
1385 static unsigned done = 0;
1386
1387 /* initialization function */
1388 void
InitCL(void)1389 InitCL(void)
1390 {
1391 if (::done++ > 0) {
1392 return;
1393 }
1394
1395 /* constitutive law array */
1396 SetCL1D("array", new CLArray1DR);
1397 SetCL3D("array", new CLArray3DR);
1398 SetCL6D("array", new CLArray6DR);
1399
1400 /* linear elastic */
1401 SetCL1D("linear" "elastic", new LinearElasticCLR<doublereal, doublereal>);
1402 SetCL3D("linear" "elastic", new LinearElasticCLR<Vec3, Mat3x3>);
1403 SetCL6D("linear" "elastic", new LinearElasticCLR<Vec6, Mat6x6>);
1404
1405 /* linear elastic isotropic */
1406 SetCL1D("linear" "elastic" "isotropic", new LinearElasticCLR<doublereal, doublereal>);
1407 SetCL3D("linear" "elastic" "isotropic", new LinearElasticCLR<Vec3, Mat3x3>);
1408 SetCL6D("linear" "elastic" "isotropic", new LinearElasticCLR<Vec6, Mat6x6>);
1409
1410 /* linear elastic generic */
1411 SetCL1D("linear" "elastic" "generic", new LinearElasticGenericCLR<doublereal, doublereal>);
1412 SetCL3D("linear" "elastic" "generic", new LinearElasticGenericCLR<Vec3, Mat3x3>);
1413 SetCL6D("linear" "elastic" "generic", new LinearElasticGenericCLR<Vec6, Mat6x6>);
1414
1415 /* linear (visco)elastic generic axial torsion coupling*/
1416 SetCL6D("linear" "elastic" "generic" "axial" "torsion" "coupling",
1417 new LinearElasticGenericAxialTorsionCouplingCLR<Vec6, Mat6x6>);
1418 SetCL6D("linear" "viscoelastic" "generic" "axial" "torsion" "coupling",
1419 new LinearViscoElasticGenericAxialTorsionCouplingCLR<Vec6, Mat6x6>);
1420
1421 /* inverse square elastic */
1422 SetCL1D("inverse" "square" "elastic", new InverseSquareElasticCLR);
1423
1424 /* log elastic */
1425 SetCL1D("log" "elastic", new LogElasticCLR<doublereal, doublereal>);
1426
1427 /* double linear elastic */
1428 SetCL1D("double" "linear" "elastic", new DoubleLinearElasticCLR<doublereal, doublereal>);
1429 SetCL3D("double" "linear" "elastic", new DoubleLinearElasticCLR<Vec3, Mat3x3>);
1430
1431 /* isotropic hardening elastic */
1432 SetCL1D("isotropic" "hardening" "elastic", new IsotropicHardeningCLR<doublereal, doublereal>);
1433 SetCL3D("isotropic" "hardening" "elastic", new IsotropicHardeningCLR<Vec3, Mat3x3>);
1434 SetCL6D("isotropic" "hardening" "elastic", new IsotropicHardeningCLR<Vec6, Mat6x6>);
1435
1436 /* contact elastic */
1437 SetCL1D("contact" "elastic", new ContactElasticCLR<doublereal, doublereal>);
1438 SetCL3D("contact" "elastic", new ContactElasticCLR<Vec3, Mat3x3>);
1439
1440 /* symbolic (elastic, viscous, viscoelastic) */
1441 SetCL1D("symbolic", new SymbolicCLR<doublereal, doublereal>);
1442 SetCL3D("symbolic", new SymbolicCLR<Vec3, Mat3x3>);
1443 SetCL6D("symbolic", new SymbolicCLR<Vec6, Mat6x6>);
1444
1445 SetCL1D("symbolic" "elastic", new SymbolicElasticCLR<doublereal, doublereal>);
1446 SetCL3D("symbolic" "elastic", new SymbolicElasticCLR<Vec3, Mat3x3>);
1447 SetCL6D("symbolic" "elastic", new SymbolicElasticCLR<Vec6, Mat6x6>);
1448
1449 SetCL1D("symbolic" "viscous", new SymbolicViscousCLR<doublereal, doublereal>);
1450 SetCL3D("symbolic" "viscous", new SymbolicViscousCLR<Vec3, Mat3x3>);
1451 SetCL6D("symbolic" "viscous", new SymbolicViscousCLR<Vec6, Mat6x6>);
1452
1453 SetCL1D("symbolic" "viscoelastic", new SymbolicViscoElasticCLR<doublereal, doublereal>);
1454 SetCL3D("symbolic" "viscoelastic", new SymbolicViscoElasticCLR<Vec3, Mat3x3>);
1455 SetCL6D("symbolic" "viscoelastic", new SymbolicViscoElasticCLR<Vec6, Mat6x6>);
1456
1457 /* linear viscous */
1458 SetCL1D("linear" "viscous", new LinearViscousCLR<doublereal, doublereal>);
1459 SetCL3D("linear" "viscous", new LinearViscousCLR<Vec3, Mat3x3>);
1460 SetCL6D("linear" "viscous", new LinearViscousCLR<Vec6, Mat6x6>);
1461
1462 /* linear viscous isotropic */
1463 SetCL1D("linear" "viscous" "isotropic", new LinearViscousCLR<doublereal, doublereal>);
1464 SetCL3D("linear" "viscous" "isotropic", new LinearViscousCLR<Vec3, Mat3x3>);
1465 SetCL6D("linear" "viscous" "isotropic", new LinearViscousCLR<Vec6, Mat6x6>);
1466
1467 /* linear viscous generic */
1468 SetCL1D("linear" "viscous" "generic", new LinearViscousGenericCLR<doublereal, doublereal>);
1469 SetCL3D("linear" "viscous" "generic", new LinearViscousGenericCLR<Vec3, Mat3x3>);
1470 SetCL6D("linear" "viscous" "generic", new LinearViscousGenericCLR<Vec6, Mat6x6>);
1471
1472 /* linear viscoelastic */
1473 SetCL1D("linear" "viscoelastic", new LinearViscoElasticCLR<doublereal, doublereal>);
1474 SetCL3D("linear" "viscoelastic", new LinearViscoElasticCLR<Vec3, Mat3x3>);
1475 SetCL6D("linear" "viscoelastic", new LinearViscoElasticCLR<Vec6, Mat6x6>);
1476
1477 /* linear viscoelastic isotropic */
1478 SetCL1D("linear" "viscoelastic" "isotropic", new LinearViscoElasticCLR<doublereal, doublereal>);
1479 SetCL3D("linear" "viscoelastic" "isotropic", new LinearViscoElasticCLR<Vec3, Mat3x3>);
1480 SetCL6D("linear" "viscoelastic" "isotropic", new LinearViscoElasticCLR<Vec6, Mat6x6>);
1481
1482 /* linear viscoelastic generic */
1483 SetCL1D("linear" "viscoelastic" "generic", new LinearViscoElasticGenericCLR<doublereal, doublereal>);
1484 SetCL3D("linear" "viscoelastic" "generic", new LinearViscoElasticGenericCLR<Vec3, Mat3x3>);
1485 SetCL6D("linear" "viscoelastic" "generic", new LinearViscoElasticGenericCLR<Vec6, Mat6x6>);
1486
1487 /* linear time variant viscoelastic generic */
1488 SetCL1D("linear" "time" "variant" "viscoelastic" "generic", new LTVViscoElasticGenericCLR<doublereal, doublereal>);
1489 SetCL3D("linear" "time" "variant" "viscoelastic" "generic", new LTVViscoElasticGenericCLR<Vec3, Mat3x3>);
1490 SetCL6D("linear" "time" "variant" "viscoelastic" "generic", new LTVViscoElasticGenericCLR<Vec6, Mat6x6>);
1491
1492 /* cubic elastic generic */
1493 SetCL1D("cubic" "elastic" "generic", new CubicElasticGenericCLR<doublereal, doublereal>);
1494 SetCL3D("cubic" "elastic" "generic", new CubicElasticGenericCLR<Vec3, Mat3x3>);
1495
1496 /* cubic viscoelastic generic */
1497 SetCL1D("cubic" "viscoelastic" "generic", new CubicViscoElasticGenericCLR<doublereal, doublereal>);
1498 SetCL3D("cubic" "viscoelastic" "generic", new CubicViscoElasticGenericCLR<Vec3, Mat3x3>);
1499
1500 /* double linear viscoelastic */
1501 SetCL1D("double" "linear" "viscoelastic", new DoubleLinearViscoElasticCLR<doublereal, doublereal>);
1502 SetCL3D("double" "linear" "viscoelastic", new DoubleLinearViscoElasticCLR<Vec3, Mat3x3>);
1503
1504 /* turbulent viscoelastic */
1505 SetCL1D("turbulent" "viscoelastic", new TurbulentViscoElasticCLR<doublereal, doublereal>);
1506
1507 /* linear elastic bistop */
1508 SetCL1D("linear" "elastic" "bistop", new LinearElasticBiStopCLR<doublereal, doublereal>);
1509 SetCL3D("linear" "elastic" "bistop", new LinearElasticBiStopCLR<Vec3, Mat3x3>);
1510 SetCL6D("linear" "elastic" "bistop", new LinearElasticBiStopCLR<Vec6, Mat6x6>);
1511
1512 /* linear viscoelastic bistop */
1513 SetCL1D("linear" "viscoelastic" "bistop", new LinearViscoElasticBiStopCLR<doublereal, doublereal>);
1514 SetCL3D("linear" "viscoelastic" "bistop", new LinearViscoElasticBiStopCLR<Vec3, Mat3x3>);
1515 SetCL6D("linear" "viscoelastic" "bistop", new LinearViscoElasticBiStopCLR<Vec6, Mat6x6>);
1516
1517 /* bistop wrapper */
1518 SetCL1D("bistop", new BiStopCLW1DR);
1519 SetCL3D("bistop", new BiStopCLW3DR);
1520 SetCL6D("bistop", new BiStopCLW6DR);
1521
1522 /* shock absorber */
1523 SetCL1D("shock" "absorber", new ShockAbsorberCLR<doublereal, doublereal>);
1524
1525 /* Artificial Neural Network */
1526 SetCL1D("ann" "elastic", new AnnElasticCLR<doublereal, doublereal>);
1527 SetCL1D("ann" "viscoelastic", new AnnViscoElasticCLR<doublereal, doublereal>);
1528
1529 /* constitutive laws sponsored by Hutchinson CdR */
1530 NLP_init();
1531 NLSF_init();
1532
1533 /* invariant constitutive law */
1534 SetCL3D("invariant" "angular", new InvAngularCLR);
1535 SetCL3D("axial" "wrapper", new AxialCLR);
1536
1537 /* ... */
1538 TDCLW_init();
1539
1540 /* NOTE: add here initialization of new built-in constitutive laws;
1541 * alternative ways to register new custom constitutive laws are:
1542 * - call SetCL*D() from anywhere in the code
1543 * - write a module that calls SetCL*D() from inside a function
1544 * called module_init(), and load it using "module load".
1545 */
1546 }
1547
1548 void
DestroyCL(void)1549 DestroyCL(void)
1550 {
1551 if (::done == 0) {
1552 silent_cerr("DestroyCL() called once too many" << std::endl);
1553 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
1554 }
1555
1556 if (--::done > 0) {
1557 return;
1558 }
1559
1560 /* free stuff */
1561 for (CL1DFuncMapType::iterator i = CL1DFuncMap.begin(); i != CL1DFuncMap.end(); ++i) {
1562 delete i->second;
1563 }
1564 CL1DFuncMap.clear();
1565
1566 for (CL3DFuncMapType::iterator i = CL3DFuncMap.begin(); i != CL3DFuncMap.end(); ++i) {
1567 delete i->second;
1568 }
1569 CL3DFuncMap.clear();
1570
1571 for (CL6DFuncMapType::iterator i = CL6DFuncMap.begin(); i != CL6DFuncMap.end(); ++i) {
1572 delete i->second;
1573 }
1574 CL6DFuncMap.clear();
1575 }
1576