MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
module-journal_bearing.cc
Go to the documentation of this file.
1 /*
2  * MBDyn (C) is a multibody analysis code.
3  * http://www.mbdyn.org
4  *
5  * Copyright (C) 1996-2017
6  *
7  * Pierangelo Masarati <masarati@aero.polimi.it>
8  * Paolo Mantegazza <mantegazza@aero.polimi.it>
9  *
10  * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
11  * via La Masa, 34 - 20156 Milano, Italy
12  * http://www.aero.polimi.it
13  *
14  * Changing this copyright notice is forbidden.
15  *
16  * This program is free software; you can redistribute it and/or modify
17  * it under the terms of the GNU General Public License as published by
18  * the Free Software Foundation (version 2 of the License).
19  *
20  *
21  * This program is distributed in the hope that it will be useful,
22  * but WITHOUT ANY WARRANTY; without even the implied warranty of
23  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24  * GNU General Public License for more details.
25  *
26  * You should have received a copy of the GNU General Public License
27  * along with this program; if not, write to the Free Software
28  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
29  */
30 
31 /*
32  AUTHOR: Reinhard Resch <R.Resch@secop.com>
33  Copyright (C) 2015(-2017) all rights reserved.
34 
35  The copyright of this code is transferred
36  to Pierangelo Masarati and Paolo Mantegazza
37  for use in the software MBDyn as described
38  in the GNU Public License version 2.1
39 */
40 
41 #include <limits>
42 #include <iostream>
43 #include <iomanip>
44 #include <cfloat>
45 #include <cassert>
46 #include <cmath>
47 #include <cstring>
48 #include <ctime>
49 
50 #ifdef HAVE_CONFIG_H
51 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
52 #endif /* HAVE_CONFIG_H */
53 
54 #include <dataman.h>
55 #include <userelem.h>
56 
57 #include "module-journal_bearing.h"
58 
59 #ifdef USE_AUTODIFF
60 #include <gradient.h>
61 #include <matvec.h>
62 #include <matvecass.h>
63 
64 using namespace grad;
65 
66 
67 class JournalBearing: virtual public Elem, public UserDefinedElem
68 {
69 public:
70  JournalBearing(unsigned uLabel, const DofOwner *pDO,
71  DataManager* pDM, MBDynParser& HP);
72  virtual ~JournalBearing(void);
73  virtual unsigned int iGetNumDof(void) const;
74  virtual DofOrder::Order GetDofType(unsigned int i) const;
75  virtual DofOrder::Order GetEqType(unsigned int i) const;
76  virtual std::ostream& DescribeDof(std::ostream& out, const char *prefix, bool bInitial) const;
77  virtual std::ostream& DescribeEq(std::ostream& out, const char *prefix, bool bInitial) const;
78  virtual unsigned int iGetNumPrivData(void) const;
79  virtual unsigned int iGetPrivDataIdx(const char *s) const;
80  virtual doublereal dGetPrivData(unsigned int i) const;
81  virtual void Output(OutputHandler& OH) const;
82  virtual void WorkSpaceDim(integer* piNumRows, integer* piNumCols) const;
83  VariableSubMatrixHandler&
84  AssJac(VariableSubMatrixHandler& WorkMat,
85  doublereal dCoef,
86  const VectorHandler& XCurr,
87  const VectorHandler& XPrimeCurr);
88  SubVectorHandler&
89  AssRes(SubVectorHandler& WorkVec,
90  doublereal dCoef,
91  const VectorHandler& XCurr,
92  const VectorHandler& XPrimeCurr);
93  template <typename T>
94  inline void
95  AssRes(GradientAssVec<T>& WorkVec,
96  doublereal dCoef,
97  const GradientVectorHandler<T>& XCurr,
98  const GradientVectorHandler<T>& XPrimeCurr,
99  enum FunctionCall func);
100  int iGetNumConnectedNodes(void) const;
101  void GetConnectedNodes(std::vector<const Node *>& connectedNodes) const;
102  void SetValue(DataManager *pDM, VectorHandler& X, VectorHandler& XP,
103  SimulationEntity::Hints *ph);
104  std::ostream& Restart(std::ostream& out) const;
105  virtual unsigned int iGetInitialNumDof(void) const;
106  virtual void
107  InitialWorkSpaceDim(integer* piNumRows, integer* piNumCols) const;
108  VariableSubMatrixHandler&
109  InitialAssJac(VariableSubMatrixHandler& WorkMat,
110  const VectorHandler& XCurr);
111  SubVectorHandler&
112  InitialAssRes(SubVectorHandler& WorkVec, const VectorHandler& XCurr);
113  template <typename T>
114  inline void
115  InitialAssRes(GradientAssVec<T>& WorkVec,
116  const GradientVectorHandler<T>& XCurr,
117  enum FunctionCall func);
118 
119  private:
120  inline void SaveLambda(const grad::Vector<doublereal, 2>& lambda);
121  template <grad::index_type N>
122  inline void SaveLambda(const grad::Vector<grad::Gradient<N>, 2>&) {}
123  inline void SaveFriction(doublereal omega, doublereal mf);
124  template <grad::index_type N>
125  inline void SaveFriction(const grad::Gradient<N>&, const grad::Gradient<N>&) {}
126 
127  StructNode* pNode1;
128  Vector<doublereal, 3> o1;
129  Matrix<doublereal, 3, 3> e;
130  StructNode* pNode2;
131  Vector<doublereal, 3> o2;
132 
133  grad::Vector<doublereal, 2> lambda;
134  doublereal z, zP, omega, mf;
135  doublereal d;
136  doublereal sigma0, sigma1;
137  doublereal muc, mus, vs, a, kv;
138  LocalDofMap dof;
139 
140  static const index_type iWorkSpace = 14 + 1;
141  static const index_type iInitialWorkSpace = 28;
142 };
143 
144 JournalBearing::JournalBearing(
145  unsigned uLabel, const DofOwner *pDO,
146  DataManager* pDM, MBDynParser& HP)
147 : Elem(uLabel, flag(0)),
148  UserDefinedElem(uLabel, pDO),
149  pNode1(0),
150  o1(Zero3),
151  e(Eye3),
152  pNode2(0),
153  o2(Zero3),
154  z(0.),
155  zP(0.),
156  omega(0.),
157  mf(0.),
158  d(0.),
159  sigma0(0.),
160  sigma1(0.),
161  muc(0.),
162  mus(0),
163  vs(1.),
164  a(1.),
165  kv(0.)
166 {
167  // help
168  if (HP.IsKeyWord("help")) {
169  silent_cout(
170  "\n"
171  "Module: InlineAD\n"
172  "\n"
173  " This element implements a journal bearing with lugre friction\n"
174  "\n"
175  " journal bearing,\n"
176  " node1, (label) <node1>,\n"
177  " [ offset, (Vec3) <offset>, ]\n"
178  " [ hinge, (Mat3x3) <orientation>, ]\n"
179  " node2, (label) <node2>,\n"
180  " [ offset, (Vec3) <offset>, ]\n"
181  " [friction, "
182  " diameter, (real) <d>, \n"
183  " coulomb friction coefficient, (real) <muc>,\n"
184  " [static friction coefficient, (real) <mus>,]\n"
185  " micro stick stiffness, (real) <sigma0>\n"
186  " [,micro stick damping, (real) <sigma1>]\n"
187  " [,sliding velocity coefficient, (real) <vs>]\n"
188  " [,sliding velocity exponent, (real) <a>]\n"
189  " [,viscous friction coefficient, (real) <kv>]\n"
190  "\n"
191  << std::endl);
192 
193  if (!HP.IsArg()) {
194  /*
195  * Exit quietly if nothing else is provided
196  */
197  throw NoErr(MBDYN_EXCEPT_ARGS);
198  }
199  }
200 
201  if ( !HP.IsKeyWord("node1") ) {
202  silent_cerr("journal bearing(" << GetLabel() << "): keyword \"node1\" expected at line " << HP.GetLineData() << std::endl);
203  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
204  }
205 
206  pNode1 = dynamic_cast<StructNode*>(pDM->ReadNode(HP,Node::STRUCTURAL));
207 
208  if (!pNode1) {
209  silent_cerr("journal bearing(" << GetLabel() << "): structural node expected at line " << HP.GetLineData() << std::endl);
210  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
211  }
212 
213  const ReferenceFrame refNode1(pNode1);
214 
215  if (HP.IsKeyWord("offset")) {
216  o1 = HP.GetPosRel(refNode1);
217  }
218 
219  if (HP.IsKeyWord("hinge") || HP.IsKeyWord("orientation")) {
220  e = HP.GetRotRel(refNode1);
221  }
222 
223  if (!HP.IsKeyWord("node2")) {
224  silent_cerr("journal bearing(" << GetLabel() << "): keyword \"node2\" expected at line " << HP.GetLineData() << std::endl);
225  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
226  }
227 
228  pNode2 = dynamic_cast<StructNode*>(pDM->ReadNode(HP,Node::STRUCTURAL));
229 
230  if (!pNode2) {
231  silent_cerr("journal bearing(" << GetLabel() << "): structural node expected at line " << HP.GetLineData() << std::endl);
232  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
233  }
234 
235  if (HP.IsKeyWord("offset")) {
236  const ReferenceFrame refNode2(pNode2);
237 
238  o2 = HP.GetPosRel(refNode2);
239  }
240 
241  if (HP.IsKeyWord("friction")) {
242  if (!HP.IsKeyWord("diameter")) {
243  silent_cerr("journal bearing(" << GetLabel()
244  << "): keyword \"diameter\" expected at line "
245  << HP.GetLineData() << std::endl);
246  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
247  }
248 
249  d = HP.GetReal();
250 
251  if (d <= 0.) {
252  silent_cerr("journal bearing(" << GetLabel()
253  << "): \"diameter\" must be greater than zero at line "
254  << HP.GetLineData() << std::endl);
255  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
256  }
257 
258  if (!HP.IsKeyWord("coulomb" "friction" "coefficient")) {
259  silent_cerr("journal bearing(" << GetLabel()
260  << "): keyword \"coulomb friction coefficient\" expected at line "
261  << HP.GetLineData() << std::endl);
262  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
263  }
264 
265  muc = HP.GetReal();
266 
267  if (muc < 0) {
268  silent_cerr("\"coulomb friction coefficient\" "
269  "must be greater than zero at line "
270  << HP.GetLineData() << std::endl);
271  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
272  }
273 
274  if (HP.IsKeyWord("static" "friction" "coefficient")) {
275  mus = HP.GetReal();
276  } else {
277  mus = muc;
278  }
279 
280  if (mus < muc) {
281  silent_cerr("journal bearing(" << GetLabel()
282  << "): \"static friction coefficient\" must be greater "
283  "than or equal to \"coulomb friction coefficient\" "
284  "at line " << HP.GetLineData() << std::endl);
285  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
286  }
287 
288  if (!HP.IsKeyWord("micro" "stick" "stiffness")) {
289  silent_cerr("journal bearing(" << GetLabel()
290  << "): keyword \"micro stick stiffness\" expected at line "
291  << HP.GetLineData() << std::endl);
292  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
293  }
294 
295  sigma0 = HP.GetReal();
296 
297  if (sigma0 <= 0.) {
298  silent_cerr("journal bearing(" << GetLabel()
299  << "): \"micro stick stiffness\" must be greater than zero at line "
300  << HP.GetLineData() << std::endl);
301  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
302  }
303 
304  if (HP.IsKeyWord("micro" "stick" "damping")) {
305  sigma1 = HP.GetReal();
306  }
307 
308  if (sigma1 < 0.) {
309  silent_cerr("journal bearing(" << GetLabel()
310  << "): \"micro stick damping\" must be greater than or equal to zero at line "
311  << HP.GetLineData() << std::endl);
312  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
313  }
314 
315  if (HP.IsKeyWord("sliding" "velocity" "coefficient")) {
316  vs = HP.GetReal();
317  }
318 
319  if (vs <= 0) {
320  silent_cerr("journal bearing(" << GetLabel()
321  << "): \"sliding velocity coefficient "
322  "must be greater than zero at line "
323  << HP.GetLineData() << std::endl);
324  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
325  }
326 
327  if (HP.IsKeyWord("sliding" "velocity" "exponent")) {
328  a = HP.GetReal();
329  }
330 
331  if (HP.IsKeyWord("viscous" "friction" "coefficient")) {
332  kv = HP.GetReal();
333  }
334 
335  if (kv < 0.) {
336  silent_cerr("journal bearing(" << GetLabel()
337  << "): \"viscous friction coefficient\" "
338  "must be greater than or equal to zero at line "
339  << HP.GetLineData() << std::endl);
340  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
341  }
342 
343  if (HP.IsKeyWord("initialize") && HP.GetYesNoOrBool()) {
344  using namespace grad;
345  const Mat3x3& R1 = pNode1->GetRCurr();
346  const Vec3& omega1 = pNode1->GetWCurr();
347  const Vec3& omega2 = pNode2->GetWCurr();
348  const Vector<doublereal, 3> domega(R1.MulTV(omega2 - omega1));
349  const double domega_proj = Dot(e.GetCol(1), domega);
350  const double v = (0.5 * d) * domega_proj;
351 
352  if (fabs(v) > 0.) {
353  const double g = muc + (mus - muc) * exp(-pow(fabs(v / vs), a));
354  z = v * g / (sigma0 * fabs(v));
355  }
356  }
357  }
358 
359  SetOutputFlag(pDM->fReadOutput(HP, Elem::LOADABLE));
360 
361  std::ostream& out = pDM->GetLogFile();
362 
363  out << "journal bearing: " << GetLabel() << " "
364  << pNode1->GetLabel() << " "
365  << o1 << " "
366  << e << " "
367  << pNode2->GetLabel() << " " << o2 << " "
368  << std::endl;
369 }
370 
371 JournalBearing::~JournalBearing(void)
372 {
373  // destroy private data
374 }
375 
376 unsigned int JournalBearing::iGetNumDof(void) const
377 {
378  return 2u + (muc > 0 ? 1u : 0u);
379 }
380 
381 DofOrder::Order JournalBearing::GetDofType(unsigned int i) const
382 {
383  switch (i) {
384  case 0:
385  case 1:
386  return DofOrder::ALGEBRAIC;
387 
388  case 2:
389  return DofOrder::DIFFERENTIAL;
390 
391  default:
392  ASSERT(0);
393  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
394  }
395 }
396 
397 DofOrder::Order JournalBearing::GetEqType(unsigned int i) const
398 {
399  switch (i) {
400  case 0:
401  case 1:
402  return DofOrder::ALGEBRAIC;
403 
404  case 2:
405  return DofOrder::DIFFERENTIAL;
406 
407  default:
408  ASSERT(0);
409  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
410  }
411 }
412 
413 std::ostream& JournalBearing::DescribeDof(std::ostream& out, const char *prefix, bool bInitial) const
414 {
415  const integer iFirstIndex = iGetFirstIndex();
416 
417  out << prefix << iFirstIndex + 1 << "->" << iFirstIndex + 2 << ": reaction forces [lambda1, lambda2]" << std::endl;
418 
419  if (bInitial) {
420  out << prefix << iFirstIndex + 3 << "->" << iFirstIndex + 4 << ": reaction force derivatives [lambdaP1, lambdaP2]" << std::endl;
421  } else if (muc > 0) {
422  out << prefix << iFirstIndex + 3 << "->" << iFirstIndex + 3 << ": sticktion state [z]" << std::endl;
423  }
424 
425  return out;
426 }
427 
428 std::ostream& JournalBearing::DescribeEq(std::ostream& out, const char *prefix, bool bInitial) const
429 {
430  const integer iFirstIndex = iGetFirstIndex();
431 
432  out << prefix << iFirstIndex + 1 << "->" << iFirstIndex + 2 << ": position constraints [c1, c2]" << std::endl;
433 
434  if (bInitial) {
435  out << prefix << iFirstIndex + 3 << "->" << iFirstIndex + 4 << ": velocity constraints [cP1, cP2]" << std::endl;
436  } else if (muc > 0) {
437  out << prefix << iFirstIndex + 3 << "->" << iFirstIndex + 3 << ": sticktion state equation [f(z, zP)]" << std::endl;
438  }
439 
440  return out;
441 }
442 
443 unsigned int JournalBearing::iGetNumPrivData(void) const
444 {
445  return 7u;
446 }
447 
448 unsigned int JournalBearing::iGetPrivDataIdx(const char *s) const
449 {
450  static const struct {
451  unsigned int index;
452  char name[8];
453  } data[] = {
454  { 1u, "lambda1" },
455  { 2u, "lambda2" },
456  { 3u, "z" },
457  { 4u, "zP" },
458  { 5u, "omega" },
459  { 6u, "mf" },
460  { 7u, "Pf" }
461  };
462 
463  const int N = iGetNumPrivData();
464 
465  ASSERT(N <= sizeof(data) / sizeof(data[0]));
466 
467  for (int i = 0; i < N; ++i) {
468  if (0 == strcmp(data[i].name, s)) {
469  return data[i].index;
470  }
471  }
472 
473  return 0;
474 }
475 
476 doublereal JournalBearing::dGetPrivData(unsigned int i) const
477 {
478  switch (i) {
479  case 1:
480  case 2:
481  return lambda(i);
482 
483  case 3:
484  return z;
485 
486  case 4:
487  return zP;
488 
489  case 5:
490  return omega;
491 
492  case 6:
493  return mf;
494 
495  case 7:
496  return -omega * mf;
497 
498  default:
499  silent_cerr("journal bearing(" << GetLabel() << "): invalid private data index " << i << std::endl);
500  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
501  }
502 }
503 
504 void
505 JournalBearing::Output(OutputHandler& OH) const
506 {
507  if ( bToBeOutput() )
508  {
509  if ( OH.UseText(OutputHandler::LOADABLE) )
510  {
511  std::ostream& os = OH.Loadable();
512 
513  os << std::setw(8) << GetLabel();
514 
515  for (int i = 1; i <= lambda.iGetNumRows(); ++i)
516  {
517  os << " " << lambda(i);
518  }
519 
520  if (muc > 0)
521  {
522  os << " " << z << " " << zP << " " << omega << " " << mf << " " << -omega * mf;
523  }
524 
525  os << std::endl;
526  }
527  }
528 }
529 
530 void
531 JournalBearing::WorkSpaceDim(integer* piNumRows, integer* piNumCols) const
532 {
533  *piNumRows = *piNumCols = iWorkSpace;
534 }
535 
536 VariableSubMatrixHandler&
537 JournalBearing::AssJac(VariableSubMatrixHandler& WorkMat,
538  doublereal dCoef,
539  const VectorHandler& XCurr,
540  const VectorHandler& XPrimeCurr)
541 {
542  GradientAssVec<Gradient<iWorkSpace> >::AssJac(this,
543  WorkMat.SetSparse(),
544  dCoef,
545  XCurr,
546  XPrimeCurr,
547  REGULAR_JAC,
548  &dof);
549 
550  return WorkMat;
551 }
552 
553 
554 SubVectorHandler&
555 JournalBearing::AssRes(SubVectorHandler& WorkVec,
556  doublereal dCoef,
557  const VectorHandler& XCurr,
558  const VectorHandler& XPrimeCurr)
559 {
560  GradientAssVec<doublereal>::AssRes(this,
561  WorkVec,
562  dCoef,
563  XCurr,
564  XPrimeCurr,
565  REGULAR_RES);
566 
567  return WorkVec;
568 }
569 
570 template <typename T>
571 inline void
572 JournalBearing::AssRes(GradientAssVec<T>& WorkVec,
573  doublereal dCoef,
574  const GradientVectorHandler<T>& XCurr,
575  const GradientVectorHandler<T>& XPrimeCurr,
576  enum FunctionCall func) {
577 
578  typedef Vector<T, 3> Vec3;
579  typedef Matrix<T, 3, 3> Mat3x3;
580 
581  const integer iFirstIndex = iGetFirstIndex();
582  const integer iFirstMomentumIndexNode1 = pNode1->iGetFirstMomentumIndex();
583  const integer iFirstMomentumIndexNode2 = pNode2->iGetFirstMomentumIndex();
584 
585  Vec3 X1, X2;
586  Mat3x3 R1, R2;
587 
588  pNode1->GetXCurr(X1, dCoef, func, &dof);
589  pNode1->GetRCurr(R1, dCoef, func, &dof);
590  pNode2->GetXCurr(X2, dCoef, func, &dof);
591  pNode2->GetRCurr(R2, dCoef, func, &dof);
592 
593  Vector<T, 2> lambda;
594 
595  XCurr.GetVec(iFirstIndex + 1, lambda, 1., &dof); // Note: for algebraic variables dCoef is always one
596 
597  SaveLambda(lambda);
598 
599  const Vec3 R2o2 = R2 * o2;
600  const Vec3 l1 = X2 + R2o2 - X1;
601 
602  const Vec3 F1 = R1 * Vec3(e.GetCol(2) * lambda(1) + e.GetCol(3) * lambda(2));
603  Vec3 M1 = Cross(l1, F1);
604  const Vec3 F2 = -F1;
605  Vec3 M2 = Cross(R2o2, F2);
606 
607  if (muc > 0 || kv > 0)
608  {
609  Vec3 omega1, omega2;
610 
611  pNode1->GetWCurr(omega1, dCoef, func, &dof);
612  pNode2->GetWCurr(omega2, dCoef, func, &dof);
613 
614  const T domega = Dot(e.GetCol(1), Transpose(R1) * Vec3(omega2 - omega1));
615  T mf = kv * domega;
616 
617  if (muc > 0)
618  {
619  T z, zP;
620 
621  XCurr.dGetCoef(iFirstIndex + 3, z, dCoef, &dof);
622  XPrimeCurr.dGetCoef(iFirstIndex + 3, zP, 1., &dof);
623 
624  T g;
625 
626  const T v = (0.5 * d) * domega;
627 
628  if (v != 0.) {
629  g = muc + (mus - muc) * exp(-pow(fabs(v / vs), a));
630  } else {
631  g = mus;
632  }
633 
634  const T f = v - sigma0 * fabs(v) / g * z - zP;
635  const T mu = sigma0 * z + sigma1 * zP;
636 
637  if (lambda(1) != 0. || lambda(2) != 0) {
638  mf += (0.5 * d) * mu * sqrt(lambda(1) * lambda(1) + lambda(2) * lambda(2));
639  }
640 
641  WorkVec.AddItem(iFirstIndex + 3, f);
642  }
643 
644  SaveFriction(domega, mf);
645 
646  const Vec3 Mf = (R1 * e.GetCol(1)) * mf;
647 
648  M1 += Mf;
649  M2 -= Mf;
650  }
651 
652  WorkVec.AddItem(iFirstMomentumIndexNode1 + 1, F1);
653  WorkVec.AddItem(iFirstMomentumIndexNode1 + 4, M1);
654  WorkVec.AddItem(iFirstMomentumIndexNode2 + 1, F2);
655  WorkVec.AddItem(iFirstMomentumIndexNode2 + 4, M2);
656 
657  const Vec3 a1 = Transpose(R1) * l1 - o1;
658 
659  for (integer i = 1; i <= 2; ++i) {
660  WorkVec.AddItem(iFirstIndex + i, Dot(e.GetCol(i + 1), a1) / dCoef);
661  }
662 }
663 
664 void JournalBearing::SaveLambda(const grad::Vector<doublereal, 2>& lambda)
665 {
666  this->lambda = lambda;
667 }
668 
669 void JournalBearing::SaveFriction(doublereal omega, doublereal mf)
670 {
671  this->omega = omega;
672  this->mf = mf;
673 }
674 
675 int
676 JournalBearing::iGetNumConnectedNodes(void) const
677 {
678  return 2;
679 }
680 
681 void
682 JournalBearing::GetConnectedNodes(std::vector<const Node *>& connectedNodes) const
683 {
684  connectedNodes.resize(iGetNumConnectedNodes());
685  connectedNodes[0] = pNode1;
686  connectedNodes[1] = pNode2;
687 }
688 
689 void
690 JournalBearing::SetValue(DataManager *pDM,
691  VectorHandler& X, VectorHandler& XP,
692  SimulationEntity::Hints *ph)
693 {
694  const integer iFirstIndex = iGetFirstIndex();
695 
696  int i;
697 
698  for (i = 1; i <= 2; ++i) {
699  X.PutCoef(iFirstIndex + i, lambda(i));
700  }
701 
702  if (muc > 0.) {
703  X.PutCoef(iFirstIndex + i, z);
704  XP.PutCoef(iFirstIndex + i, zP);
705  }
706 }
707 
708 std::ostream&
709 JournalBearing::Restart(std::ostream& out) const
710 {
711  return out;
712 }
713 
714 unsigned int
715 JournalBearing::iGetInitialNumDof(void) const
716 {
717  return 4u;
718 }
719 
720 void
721 JournalBearing::InitialWorkSpaceDim(
722  integer* piNumRows,
723  integer* piNumCols) const
724 {
725  *piNumRows = *piNumCols = iInitialWorkSpace;
726 }
727 
728 VariableSubMatrixHandler&
729 JournalBearing::InitialAssJac(
730  VariableSubMatrixHandler& WorkMat,
731  const VectorHandler& XCurr)
732 {
733 
734  GradientAssVec<Gradient<iInitialWorkSpace> >::InitialAssJac(this,
735  WorkMat.SetSparse(),
736  XCurr,
737  INITIAL_ASS_JAC,
738  &dof);
739 
740  return WorkMat;
741 }
742 
743 SubVectorHandler&
744 JournalBearing::InitialAssRes(
745  SubVectorHandler& WorkVec,
746  const VectorHandler& XCurr)
747 {
748  GradientAssVec<doublereal>::InitialAssRes(this,
749  WorkVec,
750  XCurr,
751  INITIAL_ASS_RES);
752 
753  return WorkVec;
754 }
755 
756 template <typename T>
757 inline void
758 JournalBearing::InitialAssRes(GradientAssVec<T>& WorkVec,
759  const GradientVectorHandler<T>& XCurr,
760  enum FunctionCall func) {
761 
762  typedef Vector<T, 3> Vec3;
763  typedef Vector<T, 2> Vec2;
764  typedef Matrix<T, 3, 3> Mat3x3;
765 
766  Vec3 X1, XP1, X2, XP2, omega1, omega2;
767  Mat3x3 R1, R2;
768 
769  pNode1->GetXCurr(X1, 1., func, &dof); // Note: during initial assembly dCoef is always one
770  pNode1->GetRCurr(R1, 1., func, &dof);
771  pNode1->GetVCurr(XP1, 1., func, &dof);
772  pNode1->GetWCurr(omega1, 1., func, &dof);
773 
774  pNode2->GetXCurr(X2, 1., func, &dof);
775  pNode2->GetRCurr(R2, 1., func, &dof);
776  pNode2->GetVCurr(XP2, 1., func, &dof);
777  pNode2->GetWCurr(omega2, 1., func, &dof);
778 
779  const integer iFirstIndexNode1 = pNode1->iGetFirstIndex();
780  const integer iFirstIndexNode2 = pNode2->iGetFirstIndex();
781  const integer iFirstIndex = iGetFirstIndex();
782 
783  Vec2 lambda, lambdaP;
784 
785  for (integer i = 1; i <= 2; ++i) {
786  XCurr.dGetCoef(iFirstIndex + i, lambda(i), 1., &dof);
787  XCurr.dGetCoef(iFirstIndex + i + 2, lambdaP(i), 1., &dof);
788  }
789 
790  SaveLambda(lambda);
791 
792  const Vec3 R2o2 = R2 * o2;
793  const Vec3 l1 = X2 + R2o2 - X1;
794 
795  const Vec3 F1 = R1 * Vec3(e.GetCol(2) * lambda(1) + e.GetCol(3) * lambda(2));
796  const Vec3 M1 = Cross(l1, F1);
797  const Vec3 FP1 = Cross(omega1, R1 * Vec3(e.GetCol(2) * lambda(1) + e.GetCol(3) * lambda(2)))
798  + R1 * Vec3(e.GetCol(2) * lambdaP(1) + e.GetCol(3) * lambdaP(2));
799  const Vec3 MP1 = -Cross(F1, XP2 + Cross(omega2, R2o2) - XP1) + Cross(l1, FP1);
800  const Vec3 F2 = -F1;
801  const Vec3 M2 = Cross(R2o2, F2);
802  const Vec3 FP2 = -FP1;
803  const Vec3 MP2 = Cross(Cross(omega2, R2o2), F2) + Cross(R2o2, FP2);
804 
805  const Vec3 a = Transpose(R1) * l1 - o1;
806  const Vec3 aP = Transpose(R1) * Vec3(Cross(l1, omega1) + XP2 + Cross(omega2, R2o2) - XP1);
807 
808  WorkVec.AddItem(iFirstIndexNode1 + 1, F1);
809  WorkVec.AddItem(iFirstIndexNode1 + 4, M1);
810  WorkVec.AddItem(iFirstIndexNode1 + 7, FP1);
811  WorkVec.AddItem(iFirstIndexNode1 + 10, MP1);
812 
813  WorkVec.AddItem(iFirstIndexNode2 + 1, F2);
814  WorkVec.AddItem(iFirstIndexNode2 + 4, M2);
815  WorkVec.AddItem(iFirstIndexNode2 + 7, FP2);
816  WorkVec.AddItem(iFirstIndexNode2 + 10, MP2);
817 
818  for (int i = 1; i <= 2; ++i) {
819  WorkVec.AddItem(iFirstIndex + i, Dot(e.GetCol(i + 1), a));
820  WorkVec.AddItem(iFirstIndex + i + 2, Dot(e.GetCol(i + 1), aP));
821  }
822 }
823 #endif
824 
825 bool journal_bearing_set(void)
826 {
827 #ifdef USE_AUTODIFF
828  UserDefinedElemRead *rf = new UDERead<JournalBearing>;
829 
830  if (!SetUDE("journal" "bearing", rf))
831  {
832  delete rf;
833  return false;
834  }
835 
836  return true;
837 #else
838  return false;
839 #endif
840 }
841 
842 #ifndef STATIC_MODULES
843 
844 extern "C"
845 {
846 
847 int module_init(const char *module_name, void *pdm, void *php)
848 {
849  if (!journal_bearing_set())
850  {
851  silent_cerr("journal_bearing: "
852  "module_init(" << module_name << ") "
853  "failed" << std::endl);
854 
855  return -1;
856  }
857 
858  return 0;
859 }
860 
861 }
862 
863 #endif // ! STATIC_MODULE
864 
865 
DataManager::fReadOutput
flag fReadOutput(MBDynParser &HP, const T &t) const
Definition: dataman.h:1064
grad::exp
GradientExpression< UnaryExpr< FuncExp, Expr > > exp(const GradientExpression< Expr > &u)
Definition: gradient.h:2975
MBDynParser::GetRotRel
Mat3x3 GetRotRel(const ReferenceFrame &rf)
Definition: mbpar.cc:1795
Zero3
const Vec3 Zero3(0., 0., 0.)
flag
long int flag
Definition: mbdyn.h:43
grad::Transpose
MatrixExpression< TransposedMatrix< MatrixDirectExpr< Matrix< T, N_rows, N_cols > > >, N_cols, N_rows > Transpose(const Matrix< T, N_rows, N_cols > &A)
Definition: matvec.h:2206
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
Vec3
Definition: matvec3.h:98
Eye3
const Mat3x3 Eye3(1., 0., 0., 0., 1., 0., 0., 0., 1.)
Mat3x3::GetCol
Vec3 GetCol(unsigned short int i) const
Definition: matvec3.h:903
grad::index_type
integer index_type
Definition: gradient.h:104
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-journal_bearing.cc:847
SimulationEntity::Hints
std::vector< Hint * > Hints
Definition: simentity.h:89
grad::fabs
GradientExpression< UnaryExpr< FuncFabs, Expr > > fabs(const GradientExpression< Expr > &u)
Definition: gradient.h:2973
HighParser::GetYesNoOrBool
virtual bool GetYesNoOrBool(bool bDefval=false)
Definition: parser.cc:1038
func
void func(const T &u, const T &v, const T &w, doublereal e, T &f)
Definition: gradienttest.cc:309
MBDynParser::GetPosRel
Vec3 GetPosRel(const ReferenceFrame &rf)
Definition: mbpar.cc:1331
dof
static const char * dof[]
Definition: drvdisp.cc:195
HighParser::IsKeyWord
virtual bool IsKeyWord(const char *sKeyWord)
Definition: parser.cc:910
OutputHandler::Loadable
std::ostream & Loadable(void) const
Definition: output.h:506
module-journal_bearing.h
grad::FunctionCall
FunctionCall
Definition: gradient.h:1018
journal_bearing_set
bool journal_bearing_set(void)
Definition: module-journal_bearing.cc:825
grad::Dot
DotTraits< VectorExprLhs, VectorExprRhs, N_rows, N_rows >::ExpressionType Dot(const VectorExpression< VectorExprLhs, N_rows > &u, const VectorExpression< VectorExprRhs, N_rows > &v)
Definition: matvec.h:3133
ASSERT
#define ASSERT(expression)
Definition: colamd.c:977
grad::sqrt
GradientExpression< UnaryExpr< FuncSqrt, Expr > > sqrt(const GradientExpression< Expr > &u)
Definition: gradient.h:2974
VectorHandler::PutCoef
virtual void PutCoef(integer iRow, const doublereal &dCoef)=0
grad::Cross
VectorExpression< VectorCrossExpr< VectorLhsExpr, VectorRhsExpr >, 3 > Cross(const VectorExpression< VectorLhsExpr, 3 > &u, const VectorExpression< VectorRhsExpr, 3 > &v)
Definition: matvec.h:3248
NoErr
Definition: except.h:79
HighParser::IsArg
virtual bool IsArg(void)
Definition: parser.cc:807
Elem
Definition: elem.h:75
DataManager::GetLogFile
std::ostream & GetLogFile(void) const
Definition: dataman.h:326
SetUDE
bool SetUDE(const std::string &s, UserDefinedElemRead *rude)
Definition: userelem.cc:97
a
static const doublereal a
Definition: hfluid_.h:289
VariableSubMatrixHandler::SetSparse
SparseSubMatrixHandler & SetSparse(void)
Definition: submat.h:1178
integer
long int integer
Definition: colamd.c:51
IncludeParser::GetLineData
virtual HighParser::ErrOut GetLineData(void) const
Definition: parsinc.cc:697
DataManager::ReadNode
Node * ReadNode(MBDynParser &HP, Node::Type type) const
Definition: dataman3.cc:2309
OutputHandler::UseText
bool UseText(int out) const
Definition: output.cc:446
HighParser::GetReal
virtual doublereal GetReal(const doublereal &dDefval=0.0)
Definition: parser.cc:1056