MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
module-cont-contact.cc
Go to the documentation of this file.
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:
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 
90  virtual ~ContContactCL(void) {
91  NO_OP;
92  };
93 
94  virtual ConstLawType::Type GetConstLawType(void) const {
95  return ConstLawType::VISCOELASTIC;
96  };
97 
98  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 
109  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 
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 
196  virtual std::ostream& OutputAppend(std::ostream& out) const {
197  return out << " " << m_dInitialEpsPrime << " " << m_dDissCoef;
198  };
199 
200 
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> *
241  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:
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 
399  virtual ~ContContact3DCL(void) {
400  NO_OP;
401  };
402 
403  virtual ConstLawType::Type GetConstLawType(void) const {
404  return ConstLawType::VISCOELASTIC;
405  };
406 
407  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 
418  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 
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 
505  virtual std::ostream& OutputAppend(std::ostream& out) const {
506  return out << " " << m_dInitialEpsPrime << " " << m_dDissCoef;
507  };
508 
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> *
550  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
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 
constltp_impl.h
ContContact3DCL::m_dInitialEpsPrimeTol
const doublereal m_dInitialEpsPrimeTol
Definition: module-cont-contact.cc:370
ContContact3DCL::m_dK
const doublereal m_dK
Definition: module-cont-contact.cc:379
Zero3
const Vec3 Zero3(0., 0., 0.)
ContContactCL::m_dK
const doublereal m_dK
Definition: module-cont-contact.cc:69
ContContactCL::AfterConvergence
virtual void AfterConvergence(const doublereal &Eps, const doublereal &EpsPrime=0.)
Definition: module-cont-contact.cc:201
grad::pow
GradientExpression< BinaryExpr< FuncPow, LhsExpr, RhsExpr > > pow(const GradientExpression< LhsExpr > &u, const GradientExpression< RhsExpr > &v)
Definition: gradient.h:2961
MBDYN_EXCEPT_ARGS
#define MBDYN_EXCEPT_ARGS
Definition: except.h:63
ContContactCL::m_type
ContContactCL::Type m_type
Definition: module-cont-contact.cc:64
Vec3
Definition: matvec3.h:98
ContContactCL::m_dInitialEpsPrimeTol
const doublereal m_dInitialEpsPrimeTol
Definition: module-cont-contact.cc:60
ContContact3DCL::pCopy
virtual ConstitutiveLaw< Vec3, Mat3x3 > * pCopy(void) const
Definition: module-cont-contact.cc:407
ContContact3DCL::m_dSign
const doublereal m_dSign
Definition: module-cont-contact.cc:369
ContContact3DCL::OutputAppend
virtual std::ostream & OutputAppend(std::ostream &out) const
Definition: module-cont-contact.cc:505
ContContactCL::~ContContactCL
virtual ~ContContactCL(void)
Definition: module-cont-contact.cc:90
NO_OP
#define NO_OP
Definition: myassert.h:74
ContContactCLR::Read
virtual ConstitutiveLaw< doublereal, doublereal > * Read(const DataManager *pDM, MBDynParser &HP, ConstLawType::Type &CLType)
Definition: module-cont-contact.cc:241
ContContact3DCL::~ContContact3DCL
virtual ~ContContact3DCL(void)
Definition: module-cont-contact.cc:399
ContContact3DCL::Restart
virtual std::ostream & Restart(std::ostream &out) const
Definition: module-cont-contact.cc:418
ContContactCL::OutputAppend
virtual std::ostream & OutputAppend(std::ostream &out) const
Definition: module-cont-contact.cc:196
ElasticConstitutiveLaw::Restart_int
virtual std::ostream & Restart_int(std::ostream &out) const
Definition: constltp_impl.h:169
ContContactCL::ContContactCL
ContContactCL(const TplDriveCaller< doublereal > *pTplDC, doublereal dSign, ContContactCL::Type type, doublereal dRest, doublereal dK, doublereal dExp, doublereal dInitialEpsPrimeTol)
Definition: module-cont-contact.cc:78
Zero3x3
const Mat3x3 Zero3x3(0., 0., 0., 0., 0., 0., 0., 0., 0.)
ContContactCL::m_dInitialEpsPrime
doublereal m_dInitialEpsPrime
Definition: module-cont-contact.cc:74
ContContactCL::m_dDissCoef
doublereal m_dDissCoef
Definition: module-cont-contact.cc:75
ContContact3DCL::m_dInitialEpsPrime
doublereal m_dInitialEpsPrime
Definition: module-cont-contact.cc:384
TplDriveOwner< doublereal >::pGetDriveCaller
TplDriveCaller< doublereal > * pGetDriveCaller(void) const
Definition: tpldrive.h:105
ContContactCL::m_dRest
const doublereal m_dRest
Definition: module-cont-contact.cc:62
false
Definition: mbdyn.h:76
HighParser::IsKeyWord
virtual bool IsKeyWord(const char *sKeyWord)
Definition: parser.cc:910
ContContactCL::GetConstLawType
virtual ConstLawType::Type GetConstLawType(void) const
Definition: module-cont-contact.cc:94
copysign
doublereal copysign(doublereal x, doublereal y)
Definition: gradient.h:97
ContContactCL::m_dExp
const doublereal m_dExp
Definition: module-cont-contact.cc:70
ContContact3DCL::AfterConvergence
virtual void AfterConvergence(const Vec3 &Eps, const Vec3 &EpsPrime=Zero3)
Definition: module-cont-contact.cc:509
ContContact3DCL::GetConstLawType
virtual ConstLawType::Type GetConstLawType(void) const
Definition: module-cont-contact.cc:403
ContContact3DCL::m_type
ContContact3DCL::Type m_type
Definition: module-cont-contact.cc:374
SetCL3D
bool SetCL3D(const char *name, ConstitutiveLawRead< Vec3, Mat3x3 > *rf)
Definition: constltp_impl.cc:101
ContContact3DCL::m_dRest
const doublereal m_dRest
Definition: module-cont-contact.cc:372
module_init
int module_init(const char *module_name, void *pdm, void *php)
This function registers our user defined element for the math parser.
Definition: module-cont-contact.cc:669
ContContactCL::Update
virtual void Update(const doublereal &Eps, const doublereal &EpsPrime)
Definition: module-cont-contact.cc:137
SetCL1D
bool SetCL1D(const char *name, ConstitutiveLawRead< doublereal, doublereal > *rf)
Definition: constltp_impl.cc:93
SAFENEWWITHCONSTRUCTOR
#define SAFENEWWITHCONSTRUCTOR(pnt, item, constructor)
Definition: mynewmem.h:698
NoErr
Definition: except.h:79
HighParser::IsArg
virtual bool IsArg(void)
Definition: parser.cc:807
ContContact3DCL::Update
virtual void Update(const Vec3 &Eps, const Vec3 &EpsPrime)
Definition: module-cont-contact.cc:446
ContContact3DCL::ContContact3DCL
ContContact3DCL(const TplDriveCaller< Vec3 > *pTplDC, doublereal dSign, ContContact3DCL::Type type, doublereal dRest, doublereal dK, doublereal dExp, doublereal dInitialEpsPrimeTol)
Definition: module-cont-contact.cc:388
ContContact3DCL::m_dExp
const doublereal m_dExp
Definition: module-cont-contact.cc:380
ContContactCL::pCopy
virtual ConstitutiveLaw< doublereal, doublereal > * pCopy(void) const
Definition: module-cont-contact.cc:98
doublereal
double doublereal
Definition: colamd.c:52
IncludeParser::GetLineData
virtual HighParser::ErrOut GetLineData(void) const
Definition: parsinc.cc:697
ContContact3DCLR::Read
virtual ConstitutiveLaw< Vec3, Mat3x3 > * Read(const DataManager *pDM, MBDynParser &HP, ConstLawType::Type &CLType)
Definition: module-cont-contact.cc:550
pHP
static std::stack< const HighParser * > pHP
Definition: parser.cc:598
TplDriveOwner::Get
T Get(void) const
Definition: tpldrive.h:113
ContContactCL::Restart
virtual std::ostream & Restart(std::ostream &out) const
Definition: module-cont-contact.cc:109
ContContactCL::m_dSign
const doublereal m_dSign
Definition: module-cont-contact.cc:59
HighParser::GetReal
virtual doublereal GetReal(const doublereal &dDefval=0.0)
Definition: parser.cc:1056