1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/modules/module-cont-contact/module-cont-contact.cc,v 1.13 2017/01/12 14:48:47 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 * Created by: Pierangelo Masarati <masarati@aero.polimi.it>
33 * Modified by: Matteo Fancello <matteo.fancello@gmail.com>
34 */
35
36 /*
37 * Version modified to allow contact without impacts, introducing a threshold for the InitialEpsPrime in order not to have the
38 * dissipative term of the reaction force govern sign, with the effect of imposing adhesive forces.
39 */
40
41 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
42
43 #include <cmath>
44 #include <cfloat>
45
46 #include "dataman.h"
47 #include "constltp_impl.h"
48
49 class ContContactCL
50 : public ElasticConstitutiveLaw<doublereal, doublereal> {
51 public:
52 enum Type {
53 CC_FLORES_ET_AL,
54 CC_HUNT_CROSSLEY,
55 CC_LANKARANI_NIKRAVESH
56 };
57
58 protected:
59 const doublereal m_dSign;
60 const doublereal m_dInitialEpsPrimeTol;
61
62 const doublereal m_dRest; // coefficient of restitution
63
64 ContContactCL::Type m_type; // Formulation of continuous contact law:
65 // flores [default]
66 // huntcrossley
67 // lankaraninikravesh
68
69 const doublereal m_dK;
70 const doublereal m_dExp;
71
72 mutable bool m_bActive; // is contact ongoing?
73 mutable bool m_bToggling; // toggle m_bActive
74 mutable doublereal m_dInitialEpsPrime; // initial contact velocity
75 mutable doublereal m_dDissCoef; // actual dissipation coefficient
76
77 public:
ContContactCL(const TplDriveCaller<doublereal> * pTplDC,doublereal dSign,ContContactCL::Type type,doublereal dRest,doublereal dK,doublereal dExp,doublereal dInitialEpsPrimeTol)78 ContContactCL(const TplDriveCaller<doublereal> *pTplDC,
79 doublereal dSign, ContContactCL::Type type, doublereal dRest,
80 doublereal dK, doublereal dExp, doublereal dInitialEpsPrimeTol)
81 : ElasticConstitutiveLaw<doublereal, doublereal>(pTplDC, 0.),
82 m_dSign(dSign), m_dInitialEpsPrimeTol(dInitialEpsPrimeTol),
83 m_dRest(dRest), m_type(type), m_dK(dK), m_dExp(dExp),
84 m_bActive(false), m_bToggling(false),
85 m_dInitialEpsPrime(0.)
86 {
87 NO_OP;
88 };
89
~ContContactCL(void)90 virtual ~ContContactCL(void) {
91 NO_OP;
92 };
93
GetConstLawType(void) const94 virtual ConstLawType::Type GetConstLawType(void) const {
95 return ConstLawType::VISCOELASTIC;
96 };
97
pCopy(void) const98 virtual ConstitutiveLaw<doublereal, doublereal>* pCopy(void) const {
99 ConstitutiveLaw<doublereal, doublereal>* pCL = 0;
100
101 // pass parameters to copy constructor
102 SAFENEWWITHCONSTRUCTOR(pCL, ContContactCL,
103 ContContactCL(pGetDriveCaller()->pCopy(),
104 m_dSign, m_type, m_dRest, m_dK, m_dExp, m_dInitialEpsPrimeTol));
105
106 return pCL;
107 };
108
Restart(std::ostream & out) const109 virtual std::ostream& Restart(std::ostream& out) const {
110 out << "continuous contact"
111 << ", sign, " << m_dSign
112 << ", formulation, ";
113
114 switch (m_type) {
115 case CC_FLORES_ET_AL:
116 out << "flores";
117 break;
118
119 case CC_HUNT_CROSSLEY:
120 out << "hunt crossley";
121 break;
122
123 case CC_LANKARANI_NIKRAVESH:
124 out << "lankarani nikravesh";
125 break;
126 }
127
128 out
129 << ", restitution, " << m_dRest
130 << ", kappa, " << m_dK
131 << ", exp, " << m_dExp
132 << ", EpsPrimeTol, " << m_dInitialEpsPrimeTol
133 << ", ", ElasticConstitutiveLaw<doublereal, doublereal>::Restart_int(out);
134 return out;
135 };
136
Update(const doublereal & Eps,const doublereal & EpsPrime)137 virtual void Update(const doublereal& Eps, const doublereal& EpsPrime) {
138 ConstitutiveLaw<doublereal, doublereal>::Epsilon = Eps - ElasticConstitutiveLaw<doublereal, doublereal>::Get();
139 ConstitutiveLaw<doublereal, doublereal>::EpsilonPrime = EpsPrime;
140
141 doublereal x = m_dSign*ConstitutiveLaw<doublereal, doublereal>::Epsilon;
142 if (x < 0.) {
143 if (m_bActive) {
144 if (!m_bToggling) {
145 m_bToggling = true;
146 }
147 }
148
149 ConstitutiveLaw<doublereal, doublereal>::F = 0.;
150 ConstitutiveLaw<doublereal, doublereal>::FDE = 0.;
151 ConstitutiveLaw<doublereal, doublereal>::FDEPrime = 0.;
152
153 } else {
154 doublereal xp = m_dSign*ConstitutiveLaw<doublereal, doublereal>::EpsilonPrime;
155
156 if (!m_bActive && !m_bToggling) {
157 m_bToggling = true;
158 m_dInitialEpsPrime = std::max(xp, m_dInitialEpsPrimeTol);
159
160 switch (m_type) {
161 case CC_FLORES_ET_AL:
162 m_dDissCoef = (std::abs(xp) > m_dInitialEpsPrimeTol) * 8./5. * m_dK * (1. - m_dRest) / m_dRest / m_dInitialEpsPrime;
163 break;
164
165 case CC_HUNT_CROSSLEY:
166 m_dDissCoef = (std::abs(xp) > m_dInitialEpsPrimeTol) * 3./2. * m_dK * (1. - m_dRest) / m_dInitialEpsPrime;
167 break;
168
169 case CC_LANKARANI_NIKRAVESH:
170 m_dDissCoef = (std::abs(xp) > m_dInitialEpsPrimeTol) * 3./4. * m_dK * (1. - std::pow(m_dRest, 2)) / m_dInitialEpsPrime;
171 break;
172 }
173 }
174
175 // doublereal xn = std::pow(x, m_dExp);
176 // doublereal xnm1 = std::pow(x, m_dExp - 1.);
177
178 // doublereal xp_t = ( ( (m_dDissCoef*xp < 0.)&&(abs(m_dDissCoef*xp) > m_dK)) ? abs(m_dK/m_dDissCoef/xp) : 1. ); // Threshold
179
180 // ConstitutiveLaw<doublereal, doublereal>::F = m_dSign*(m_dK*xn + m_dDissCoef *xn* xp * xp_t); // FIXME
181 // ConstitutiveLaw<doublereal, doublereal>::FDE = m_dSign*(m_dExp*m_dK*xnm1 + m_dDissCoef *m_dExp*xnm1* xp_t); //FIXME
182 // ConstitutiveLaw<doublereal, doublereal>::FDEPrime = m_dSign*(m_dDissCoef * xn) * xp_t ;
183
184
185 doublereal xn = std::pow(x, m_dExp);
186 doublereal xnm1 = std::pow(x, m_dExp - 1.);
187
188 ConstitutiveLaw<doublereal, doublereal>::F = m_dSign*xn*(m_dK + m_dDissCoef*xp); // FIXME
189 ConstitutiveLaw<doublereal, doublereal>::FDE = m_dSign*xnm1*(m_dExp*m_dK + m_dDissCoef*m_dExp);
190 ConstitutiveLaw<doublereal, doublereal>::FDEPrime = m_dSign*xn*m_dDissCoef;
191
192 // silent_cout("hello " << std::endl);
193 }
194 };
195
OutputAppend(std::ostream & out) const196 virtual std::ostream& OutputAppend(std::ostream& out) const {
197 return out << " " << m_dInitialEpsPrime << " " << m_dDissCoef;
198 };
199
200
AfterConvergence(const doublereal & Eps,const doublereal & EpsPrime=0.)201 virtual void AfterConvergence(const doublereal& Eps, const doublereal& EpsPrime = 0.) {
202
203 #if 0
204 if (m_bActive) {
205 doublereal eps2show = m_dSign*ConstitutiveLaw<doublereal, doublereal>::Epsilon ;
206 doublereal F2show = m_dSign*ConstitutiveLaw<doublereal, doublereal>::F;
207 silent_cout( " " << eps2show << " " << F2show << std::endl);
208 }
209 #endif
210
211 if (m_bToggling) {
212
213 #if 0
214 silent_cout(">> ContContactCL::AfterConvergence() "
215 "m_bToggling=" << (m_bToggling ? "true" : "false") << " "
216 "m_bActive=" << (m_bActive ? "true" : "false") << " "
217 "m_dInitialEpsPrime=" << m_dInitialEpsPrime << std::endl);
218 #endif
219
220 if (m_bActive) {
221 m_bActive = false;
222 m_dInitialEpsPrime = 0.;
223
224 } else {
225 m_bActive = true;
226 }
227 m_bToggling = false;
228 #if 0
229 silent_cout("<< ContContactCL::AfterConvergence() "
230 "m_bToggling=" << (m_bToggling ? "true" : "false") << " "
231 "m_bActive=" << (m_bActive ? "true" : "false") << " "
232 "m_dInitialEpsPrime=" << m_dInitialEpsPrime << std::endl);
233 #endif
234 }
235 };
236 };
237
238 /* specific functional object(s) */
239 struct ContContactCLR : public ConstitutiveLawRead<doublereal, doublereal> {
240 virtual ConstitutiveLaw<doublereal, doublereal> *
ReadContContactCLR241 Read(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType) {
242 ConstitutiveLaw<doublereal, doublereal>* pCL = 0;
243
244 CLType = ConstLawType::VISCOELASTIC;
245
246 if (HP.IsKeyWord("help")) {
247 silent_cerr("ContContactCLR:\n"
248 " continuous contact ,\n"
249 " [ , sign , { negative | positive | <sign> } , ]\n"
250 " [ , formulation, { flores | lankarani nikravesh | hunt crossley } , ]\n"
251 " restitution , <restitution_coefficient> , (0 -> 1)\n"
252 " kappa , <stiffness> , (> 0)\n"
253 " exp , <exponent> , (>= 1)\n"
254 " [ , EpsPrimeTol , <EpsPrimeTol> ]\n"
255 " [ , prestrain , (DriveCaller) <prestrain> ]\n"
256 << std::endl);
257
258 if (!HP.IsArg()) {
259 throw NoErr(MBDYN_EXCEPT_ARGS);
260 }
261 }
262
263 doublereal dSign = -1.;
264 if (HP.IsKeyWord("sign")) {
265 if (HP.IsKeyWord("positive")) {
266 dSign = 1.;
267 } else if (HP.IsKeyWord("negative")) {
268 dSign = -1.;
269 } else {
270 doublereal d = HP.GetReal();
271 if (d == 0.) {
272 silent_cerr("ContContactCLR: invalid sign " << d
273 << " at line " << HP.GetLineData() << std::endl);
274 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
275 }
276 dSign = copysign(1., d);
277 }
278 }
279
280 ContContactCL::Type type(ContContactCL::CC_FLORES_ET_AL);
281 if (HP.IsKeyWord("formulation")) {
282 if (HP.IsKeyWord("flores")) {
283 type = ContContactCL::CC_FLORES_ET_AL;
284
285 } else if (HP.IsKeyWord("hunt" "crossley")) {
286 type = ContContactCL::CC_HUNT_CROSSLEY;
287
288 } else if (HP.IsKeyWord("lankarani" "nikravesh")) {
289 type = ContContactCL::CC_LANKARANI_NIKRAVESH;
290
291 } else {
292 silent_cerr("ContContactCLR: invalid formulation"
293 << " at line " << HP.GetLineData() << std::endl);
294 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
295 }
296 }
297
298 if (!HP.IsKeyWord("restitution")) {
299 silent_cerr("ContContactCLR: \"restitution\" expected at line "
300 << HP.GetLineData() << std::endl);
301 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
302 }
303 doublereal dRest = HP.GetReal();
304 if (dRest < 0.) {
305 silent_cerr("ContContactCLR: invalid \"restitution\" at line "
306 << HP.GetLineData() << std::endl);
307 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
308
309 } else if ((type == ContContactCL::CC_FLORES_ET_AL) && (dRest <= 0.)) {
310 silent_cerr("ContContactCLR: null \"restitution\" incompatible with \"Flores\" at line "
311 << HP.GetLineData() << std::endl);
312 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
313 }
314
315 if (!HP.IsKeyWord("kappa")) {
316 silent_cerr("ContContactCLR: \"kappa\" expected at line " << HP.GetLineData() << std::endl);
317 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
318 }
319 doublereal dK = HP.GetReal();
320 if (dK <= 0.) {
321 silent_cerr("ContContactCLR: invalid \"kappa\" at line " << HP.GetLineData() << std::endl);
322 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
323 }
324
325 if (!HP.IsKeyWord("exp")) {
326 silent_cerr("ContContactCLR: \"exp\" expected at line " << HP.GetLineData() << std::endl);
327 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
328 }
329 doublereal dExp = HP.GetReal();
330 if (dExp < 1.) {
331 silent_cerr("ContContactCLR: invalid \"exp\" at line " << HP.GetLineData() << std::endl);
332 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
333 }
334
335 doublereal dInitialEpsPrimeTol = 1.e-6;
336 if (HP.IsKeyWord("EpsPrimeTol")) {
337 dInitialEpsPrimeTol = HP.GetReal();
338 if (dInitialEpsPrimeTol < 0.) {
339 silent_cerr("ContContactCLR: invalid \"EpsPrimeTol\" at line " << HP.GetLineData() << std::endl);
340 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
341 }
342 }
343
344 /* Prestrain */
345 TplDriveCaller<doublereal> *pTplDC = GetPreStrain<doublereal>(pDM, HP);
346
347 SAFENEWWITHCONSTRUCTOR(pCL, ContContactCL,
348 ContContactCL(pTplDC, dSign, type, dRest, dK, dExp, dInitialEpsPrimeTol));
349
350 return pCL;
351 };
352 };
353
354 //---------------------------------------------------------------------------------------------------------------------------------
355 //---------------------------------------------------------------------------------------------------------------------------------
356 //---------------------------------------------------------------------------------------------------------------------------------
357 //---------------------------------------------------------------------------------------------------------------------------------
358
359 class ContContact3DCL
360 : public ElasticConstitutiveLaw<Vec3, Mat3x3> {
361 public:
362 enum Type {
363 CC_FLORES_ET_AL = ContContactCL::CC_FLORES_ET_AL,
364 CC_HUNT_CROSSLEY = ContContactCL::CC_HUNT_CROSSLEY,
365 CC_LANKARANI_NIKRAVESH = ContContactCL::CC_LANKARANI_NIKRAVESH
366 };
367
368 protected:
369 const doublereal m_dSign;
370 const doublereal m_dInitialEpsPrimeTol;
371
372 const doublereal m_dRest; // coefficient of restitution
373
374 ContContact3DCL::Type m_type; // Formulation of continuous contact law:
375 // flores [default]
376 // huntcrossley
377 // lankaraninikravesh
378
379 const doublereal m_dK;
380 const doublereal m_dExp;
381
382 mutable bool m_bActive; // is contact ongoing?
383 mutable bool m_bToggling; // toggle m_bActive
384 mutable doublereal m_dInitialEpsPrime; // initial contact velocity
385 mutable doublereal m_dDissCoef; // actual dissipation coefficient
386
387 public:
ContContact3DCL(const TplDriveCaller<Vec3> * pTplDC,doublereal dSign,ContContact3DCL::Type type,doublereal dRest,doublereal dK,doublereal dExp,doublereal dInitialEpsPrimeTol)388 ContContact3DCL(const TplDriveCaller<Vec3> *pTplDC,
389 doublereal dSign, ContContact3DCL::Type type, doublereal dRest,
390 doublereal dK, doublereal dExp, doublereal dInitialEpsPrimeTol)
391 : ElasticConstitutiveLaw<Vec3, Mat3x3>(pTplDC, Zero3),
392 m_dSign(dSign), m_dInitialEpsPrimeTol(dInitialEpsPrimeTol),
393 m_dRest(dRest), m_type(type), m_dK(dK), m_dExp(dExp),
394 m_bActive(false), m_bToggling(false), m_dInitialEpsPrime(0.)
395 {
396 NO_OP;
397 };
398
~ContContact3DCL(void)399 virtual ~ContContact3DCL(void) {
400 NO_OP;
401 };
402
GetConstLawType(void) const403 virtual ConstLawType::Type GetConstLawType(void) const {
404 return ConstLawType::VISCOELASTIC;
405 };
406
pCopy(void) const407 virtual ConstitutiveLaw<Vec3, Mat3x3>* pCopy(void) const {
408 ConstitutiveLaw<Vec3, Mat3x3>* pCL = 0;
409
410 // pass parameters to copy constructor
411 SAFENEWWITHCONSTRUCTOR(pCL, ContContact3DCL,
412 ContContact3DCL(pGetDriveCaller()->pCopy(),
413 m_dSign, m_type, m_dRest, m_dK, m_dExp, m_dInitialEpsPrimeTol));
414
415 return pCL;
416 };
417
Restart(std::ostream & out) const418 virtual std::ostream& Restart(std::ostream& out) const {
419 out << "continuous contact"
420 << ", sign, " << m_dSign
421 << ", formulation, ";
422
423 switch (m_type) {
424 case CC_FLORES_ET_AL:
425 out << "flores";
426 break;
427
428 case CC_HUNT_CROSSLEY:
429 out << "hunt crossley";
430 break;
431
432 case CC_LANKARANI_NIKRAVESH:
433 out << "lankarani nikravesh";
434 break;
435 }
436
437 out
438 << ", restitution, " << m_dRest
439 << ", kappa, " << m_dK
440 << ", exp, " << m_dExp
441 << ", EpsPrimeTol, " << m_dInitialEpsPrimeTol
442 << ", ", ElasticConstitutiveLaw<Vec3, Mat3x3>::Restart_int(out);
443 return out;
444 };
445
Update(const Vec3 & Eps,const Vec3 & EpsPrime)446 virtual void Update(const Vec3& Eps, const Vec3& EpsPrime) {
447 ConstitutiveLaw<Vec3, Mat3x3>::Epsilon = Eps - ElasticConstitutiveLaw<Vec3, Mat3x3>::Get();
448 ConstitutiveLaw<Vec3, Mat3x3>::EpsilonPrime = EpsPrime;
449
450 doublereal x = m_dSign*ConstitutiveLaw<Vec3, Mat3x3>::Epsilon(3);
451 if (x < 0.) {
452 if (m_bActive) {
453 if (!m_bToggling) {
454 m_bToggling = true;
455 }
456 }
457
458 ConstitutiveLaw<Vec3, Mat3x3>::F = Zero3;
459 ConstitutiveLaw<Vec3, Mat3x3>::FDE = Zero3x3;
460 ConstitutiveLaw<Vec3, Mat3x3>::FDEPrime = Zero3x3;
461
462 } else {
463 doublereal xp = m_dSign*ConstitutiveLaw<Vec3, Mat3x3>::EpsilonPrime(3);
464
465 if (!m_bActive && !m_bToggling) {
466 m_bToggling = true;
467 m_dInitialEpsPrime = std::max(xp, m_dInitialEpsPrimeTol);
468
469 switch (m_type) {
470 case CC_FLORES_ET_AL:
471 m_dDissCoef = (std::abs(xp) > m_dInitialEpsPrimeTol) * 8./5. * m_dK * (1. - m_dRest) / m_dRest / m_dInitialEpsPrime;
472 break;
473
474 case CC_HUNT_CROSSLEY:
475 m_dDissCoef = (std::abs(xp) > m_dInitialEpsPrimeTol) * 3./2. * m_dK * (1. - m_dRest) / m_dInitialEpsPrime;
476 break;
477
478 case CC_LANKARANI_NIKRAVESH:
479 m_dDissCoef = (std::abs(xp) > m_dInitialEpsPrimeTol) * 3./4. * m_dK * (1. - std::pow(m_dRest, 2)) / m_dInitialEpsPrime;
480 break;
481 }
482 }
483
484 // doublereal xn = std::pow(x, m_dExp);
485 // doublereal xnm1 = std::pow(x, m_dExp - 1.);
486
487 // doublereal xp_t = ( ( (m_dDissCoef*xp < 0.)&&(abs(m_dDissCoef*xp) > m_dK)) ? abs(m_dK/m_dDissCoef/xp) : 1. ); // Threshold
488
489 // ConstitutiveLaw<doublereal, doublereal>::F = m_dSign*(m_dK*xn + m_dDissCoef *xn* xp * xp_t); // FIXME
490 // ConstitutiveLaw<doublereal, doublereal>::FDE = m_dSign*(m_dExp*m_dK*xnm1 + m_dDissCoef *m_dExp*xnm1* xp_t); //FIXME
491 // ConstitutiveLaw<doublereal, doublereal>::FDEPrime = m_dSign*(m_dDissCoef * xn) * xp_t ;
492
493 doublereal xn = std::pow(x, m_dExp);
494 doublereal xnm1 = std::pow(x, m_dExp - 1.);
495
496 //----------------------FIXME CHECK CORRECTNESS
497 ConstitutiveLaw<Vec3, Mat3x3>::F(3) = m_dSign*xn*(m_dK + m_dDissCoef*xp); // FIXME
498 ConstitutiveLaw<Vec3, Mat3x3>::FDE(3, 3) = m_dSign*xnm1*(m_dExp*m_dK + m_dDissCoef*m_dExp);
499 ConstitutiveLaw<Vec3, Mat3x3>::FDEPrime(3, 3) = m_dSign*xn*m_dDissCoef;
500
501
502 }
503 };
504
OutputAppend(std::ostream & out) const505 virtual std::ostream& OutputAppend(std::ostream& out) const {
506 return out << " " << m_dInitialEpsPrime << " " << m_dDissCoef;
507 };
508
AfterConvergence(const Vec3 & Eps,const Vec3 & EpsPrime=Zero3)509 virtual void AfterConvergence(const Vec3& Eps, const Vec3& EpsPrime = Zero3) {
510
511 #if 0
512 if (m_bActive) {
513 doublereal eps2show = m_dSign*ConstitutiveLaw<Vec3, Mat3x3>::Epsilon(3) ;
514 doublereal F2show = m_dSign*ConstitutiveLaw<Vec3, Mat3x3>::F(3);
515 silent_cout( " " << eps2show << " " << F2show << std::endl);
516 }
517 #endif
518
519 if (m_bToggling) {
520
521 #if 0
522 silent_cout(">> ContContact3DCL::AfterConvergence() "
523 "m_bToggling=" << (m_bToggling ? "true" : "false") << " "
524 "m_bActive=" << (m_bActive ? "true" : "false") << " "
525 "m_dInitialEpsPrime=" << m_dInitialEpsPrime << std::endl);
526 #endif
527
528 if (m_bActive) {
529 m_bActive = false;
530 m_dInitialEpsPrime = 0.;
531
532 } else {
533 m_bActive = true;
534 }
535 m_bToggling = false;
536
537 #if 0
538 silent_cout("<< ContContact3DCL::AfterConvergence() "
539 "m_bToggling=" << (m_bToggling ? "true" : "false") << " "
540 "m_bActive=" << (m_bActive ? "true" : "false") << " "
541 "m_dInitialEpsPrime=" << m_dInitialEpsPrime << std::endl);
542 #endif
543 }
544 };
545 };
546
547 /* specific functional object(s) */
548 struct ContContact3DCLR : public ConstitutiveLawRead<Vec3, Mat3x3> {
549 virtual ConstitutiveLaw<Vec3, Mat3x3> *
ReadContContact3DCLR550 Read(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType) {
551 ConstitutiveLaw<Vec3, Mat3x3>* pCL = 0;
552
553 CLType = ConstLawType::VISCOELASTIC;
554
555 if (HP.IsKeyWord("help")) {
556 silent_cerr("ContContact3DCL:\n"
557 " continuous contact ,\n"
558 " [ , sign , { negative | positive | <sign> } , ]\n"
559 " [ , formulation , { flores | lankarani nikravesh | hunt crossley } , ]\n"
560 " restitution , <restitution_coefficient>, (0 -> 1)\n"
561 " kappa , <stiffness> , (> 0)\n"
562 " exp , <exponent> , (>= 1)\n"
563 " [ , EpsPrimeTol , <EpsPrimeTol> ]\n"
564 " [ , prestrain , (DriveCaller) <prestrain> ]\n"
565 << std::endl);
566
567 if (!HP.IsArg()) {
568 throw NoErr(MBDYN_EXCEPT_ARGS);
569 }
570 }
571
572 doublereal dSign = -1.;
573 if (HP.IsKeyWord("sign")) {
574 if (HP.IsKeyWord("positive")) {
575 dSign = 1.;
576 } else if (HP.IsKeyWord("negative")) {
577 dSign = -1.;
578 } else {
579 doublereal d = HP.GetReal();
580 if (d == 0.) {
581 silent_cerr("ContContact3DCLR: invalid sign " << d
582 << " at line " << HP.GetLineData() << std::endl);
583 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
584 }
585 dSign = copysign(1., d);
586 }
587 }
588
589 ContContact3DCL::Type type(ContContact3DCL::CC_FLORES_ET_AL);
590 if (HP.IsKeyWord("formulation")) {
591 if (HP.IsKeyWord("flores")) {
592 type = ContContact3DCL::CC_FLORES_ET_AL;
593
594 } else if (HP.IsKeyWord("hunt" "crossley")) {
595 type = ContContact3DCL::CC_HUNT_CROSSLEY;
596
597 } else if (HP.IsKeyWord("lankarani" "nikravesh")) {
598 type = ContContact3DCL::CC_LANKARANI_NIKRAVESH;
599
600 } else {
601 silent_cerr("ContContact3DCLR: invalid formulation"
602 << " at line " << HP.GetLineData() << std::endl);
603 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
604 }
605 }
606
607 if (!HP.IsKeyWord("restitution")) {
608 silent_cerr("ContContact3DCLR: \"restitution\" expected at line "
609 << HP.GetLineData() << std::endl);
610 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
611 }
612 doublereal dRest = HP.GetReal();
613 if (dRest < 0.) {
614 silent_cerr("ContContact3DCLR: invalid \"restitution\" at line "
615 << HP.GetLineData() << std::endl);
616 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
617
618 } else if ((type == ContContact3DCL::CC_FLORES_ET_AL) && (dRest <= 0.)) {
619 silent_cerr("ContContact3DCLR: null \"restitution\" incompatible with \"Flores\" at line "
620 << HP.GetLineData() << std::endl);
621 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
622 }
623
624 if (!HP.IsKeyWord("kappa")) {
625 silent_cerr("ContContact3DCLR: \"kappa\" expected at line " << HP.GetLineData() << std::endl);
626 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
627 }
628 doublereal dK = HP.GetReal();
629 if (dK <= 0.) {
630 silent_cerr("ContContact3DCLR: invalid \"kappa\" at line " << HP.GetLineData() << std::endl);
631 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
632 }
633
634 if (!HP.IsKeyWord("exp")) {
635 silent_cerr("ContContact3DCLR: \"exp\" expected at line " << HP.GetLineData() << std::endl);
636 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
637 }
638 doublereal dExp = HP.GetReal();
639 if (dExp < 1.) {
640 silent_cerr("ContContact3DCLR: invalid \"exp\" at line " << HP.GetLineData() << std::endl);
641 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
642 }
643
644 doublereal dInitialEpsPrimeTol = 1.e-6;
645 if (HP.IsKeyWord("EpsPrimeTol")) {
646 dInitialEpsPrimeTol = HP.GetReal();
647 if (dInitialEpsPrimeTol < 0.) {
648 silent_cerr("ContContact3DCLR: invalid \"EpsPrimeTol\" at line " << HP.GetLineData() << std::endl);
649 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
650 }
651 }
652
653 /* Prestrain */
654 TplDriveCaller<Vec3> *pTplDC = GetPreStrain<Vec3>(pDM, HP);
655
656 SAFENEWWITHCONSTRUCTOR(pCL, ContContact3DCL,
657 ContContact3DCL(pTplDC, dSign, type, dRest, dK, dExp, dInitialEpsPrimeTol));
658
659 return pCL;
660 };
661 };
662
663
664
665
666
667
668 extern "C" int
module_init(const char * module_name,void * pdm,void * php)669 module_init(const char *module_name, void *pdm, void *php)
670 {
671 #if 0
672 DataManager *pDM = (DataManager *)pdm;
673 MBDynParser *pHP = (MBDynParser *)php;
674 #endif
675
676 ConstitutiveLawRead<doublereal, doublereal> *rf1D = new ContContactCLR;
677 if (!SetCL1D("continuous" "contact", rf1D)) {
678 delete rf1D;
679
680 silent_cerr("ContContactCL: "
681 "module_init(" << module_name << ") "
682 "failed" << std::endl);
683
684 return -1;
685 }
686
687 ConstitutiveLawRead<Vec3, Mat3x3> *rf3D = new ContContact3DCLR;
688 if (!SetCL3D("continuous" "contact", rf3D)) {
689 delete rf3D;
690
691 silent_cerr("ContContact3DCL: "
692 "module_init(" << module_name << ") "
693 "failed" << std::endl);
694
695 return -1;
696 }
697
698 return 0;
699 }
700
701
702
703