MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
module-hydrodynamic_plain_bearing.cc
Go to the documentation of this file.
1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/modules/module-hydrodynamic_plain_bearing/module-hydrodynamic_plain_bearing.cc,v 1.13 2017/01/12 14:52:45 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 /*
33  AUTHOR: Reinhard Resch <R.RESCH@secop.com>
34  Copyright (C) 2011(-2017) all rights reserved.
35 
36  The copyright of this code is transferred
37  to Pierangelo Masarati and Paolo Mantegazza
38  for use in the software MBDyn as described
39  in the GNU Public License version 2.1
40 */
41 
42 #include <limits>
43 #include <iostream>
44 #include <iomanip>
45 #include <cfloat>
46 #include <cassert>
47 #include <cmath>
48 
49 #ifdef HAVE_CONFIG_H
50 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
51 #endif /* HAVE_CONFIG_H */
52 
53 #include <dataman.h>
54 #include <userelem.h>
55 
56 #include "module-hydrodynamic_plain_bearing.h"
57 #include "hydrodynamic_plain_bearing_force.h"
58 
59 class HydrodynamicPlainBearing: virtual public Elem, public UserDefinedElem
60 {
61 public:
62  HydrodynamicPlainBearing(unsigned uLabel, const DofOwner *pDO,
63  DataManager* pDM, MBDynParser& HP);
64  virtual ~HydrodynamicPlainBearing(void);
65  virtual void Output(OutputHandler& OH) const;
66  virtual void WorkSpaceDim(integer* piNumRows, integer* piNumCols) const;
67  VariableSubMatrixHandler&
68  AssJac(VariableSubMatrixHandler& WorkMat,
69  doublereal dCoef,
70  const VectorHandler& XCurr,
71  const VectorHandler& XPrimeCurr);
72  inline
73  void ComputeResidual(Vec3& e_R2,
74  Vec3& e_dot_R2,
75  doublereal omega_proj[2],
76  Vec3& F2_R2,
77  Vec3& M2_R2,
78  Vec3& F2_I,
79  Vec3& M2_I,
80  Vec3& F1_I,
81  Vec3& M1_I,
82  doublereal& eps,
83  doublereal& eps_dot,
84  doublereal& delta,
85  doublereal& SoD,
86  doublereal& SoV,
87  doublereal& my,
88  doublereal& beta,
89  doublereal r,
90  doublereal alpha) const;
91  SubVectorHandler&
92  AssRes(SubVectorHandler& WorkVec,
93  doublereal dCoef,
94  const VectorHandler& XCurr,
95  const VectorHandler& XPrimeCurr);
96  unsigned int iGetNumPrivData(void) const;
97  int iGetNumConnectedNodes(void) const;
98  void GetConnectedNodes(std::vector<const Node *>& connectedNodes) const;
99  void SetValue(DataManager *pDM, VectorHandler& X, VectorHandler& XP,
100  SimulationEntity::Hints *ph);
101  std::ostream& Restart(std::ostream& out) const;
102  virtual unsigned int iGetInitialNumDof(void) const;
103  virtual void
104  InitialWorkSpaceDim(integer* piNumRows, integer* piNumCols) const;
105  VariableSubMatrixHandler&
106  InitialAssJac(VariableSubMatrixHandler& WorkMat,
107  const VectorHandler& XCurr);
108  SubVectorHandler&
109  InitialAssRes(SubVectorHandler& WorkVec, const VectorHandler& XCurr);
110  private:
111  const StructNode* m_pShaft;
112  const StructNode* m_pBearing;
113  Vec3 m_o1_R1;
114  Vec3 m_o2_R2;
115  bearing_data m_bdat;
116  DriveOwner m_InitialAssemblyFactor;
117  integer m_iNumGaussPoints;
118  const doublereal* m_r;
119  const doublereal* m_alpha;
120  struct OutputOpt {
121  doublereal r;
122  doublereal alpha;
123  } m_output[6];
124  integer m_iNumOutputPoints;
125  bool m_lambda;
126  static const doublereal s_r1[1], s_alpha1[1];
127  static const doublereal s_r2[2], s_alpha2[2];
128  static const doublereal s_r3[3], s_alpha3[3];
129  static const doublereal s_r6[6], s_alpha6[6];
130 };
131 
132 // GAUSS-LEGENDRE integration points and weight factors according to K.J. Bathe 2002 table 5.6
133 // times 0.5 in order to account for an integration over the interval -0.5*b ... 0.5*b
134 
135 const doublereal HydrodynamicPlainBearing::s_r1[1] = {0.5 * 0.};
136 const doublereal HydrodynamicPlainBearing::s_alpha1[1] = {0.5 * 2};
137 
138 const doublereal HydrodynamicPlainBearing::s_r2[2] = {0.5 * -0.577350269189626,
139  0.5 * 0.577350269189626};
140 const doublereal HydrodynamicPlainBearing::s_alpha2[2] = {0.5 * 1.,
141  0.5 * 1.};
142 
143 const doublereal HydrodynamicPlainBearing::s_r3[3] = {0.5 * -0.774596669241483,
144  0.5 * 0.,
145  0.5 * 0.774596669241483};
146 const doublereal HydrodynamicPlainBearing::s_alpha3[3] = {0.5 * 0.555555555555556,
147  0.5 * 0.888888888888889,
148  0.5 * 0.555555555555556};
149 
150 const doublereal HydrodynamicPlainBearing::s_r6[6] = {0.5 * -0.932469514203152,
151  0.5 * -0.661209386466265,
152  0.5 * -0.238619186083197,
153  0.5 * 0.238619186083197,
154  0.5 * 0.661209386466265,
155  0.5 * 0.932469514203152};
156 const doublereal HydrodynamicPlainBearing::s_alpha6[6] = {0.5 * 0.171324492379170,
157  0.5 * 0.360761573048139,
158  0.5 * 0.467913934572691,
159  0.5 * 0.467913934572691,
160  0.5 * 0.360761573048139,
161  0.5 * 0.171324492379170};
162 
163 HydrodynamicPlainBearing::HydrodynamicPlainBearing(
164  unsigned uLabel, const DofOwner *pDO,
165  DataManager* pDM, MBDynParser& HP)
166 : Elem(uLabel, flag(0)),
167  UserDefinedElem(uLabel, pDO),
168  m_pShaft(0),
169  m_pBearing(0),
170  m_o1_R1(0.,0.,0.),
171  m_o2_R2(0.,0.,0.),
172  m_InitialAssemblyFactor(0),
173  m_iNumGaussPoints(0),
174  m_r(0),
175  m_alpha(0),
176  m_iNumOutputPoints(0),
177  m_lambda(true)
178 {
179  std::memset(&m_bdat, 0, sizeof(m_bdat));
180  std::memset(&m_output, 0, sizeof(m_output));
181 
182  // help
183  if (HP.IsKeyWord("help"))
184  {
185  silent_cout(
186  "\n"
187  "Module: hydrodynamic_plain_bearing_with_offset\n"
188  "\n"
189  " This module implements a hydrodynamic plain bearing according to\n"
190  "\n"
191  " Hans Juergen Butenschoen 1976\n"
192  " Das hydrodynamische zylindrische Gleitlager\n"
193  " endlicher Breite unter instationaerer Belastung\n"
194  "\n"
195  " hydrodynamic_plain_bearing_with_offset,\n"
196  " shaft, (label) <shaft node>,\n"
197  " [offset, (Vec3) <o1>,]\n"
198  " bearing, (label) <bearing node>,\n"
199  " [offset, (Vec3) <o2>,]\n"
200  " bearing_width, (real) <b>,\n"
201  " {shaft diameter, (real) <d> | bearing_diameter, (real) <D>,}\n"
202  " relative_clearance, (real) <Psi>,\n"
203  " oil_viscosity, (real) <eta>,\n"
204  " initial_assembly_factor, (DriveCaller),\n"
205  " [number of gauss points, <num_gauss_points>]\n"
206  " [output points, (integer) <num_output_points> [, {gauss point, (integer) <index_1> | custom, r, (real) <r_1>, alpha, (real) <alpha_1>}]]\n"
207  " [extend shaft to bearing center, {yes | no | (bool) <extend>}]\n"
208  "\n"
209  " b ... bearing width [m]\n"
210  " d ... bearing diameter [m]\n"
211  " Psi ... relative clearance Psi = ( D - d ) / D [m/m]\n"
212  " eta ... dynamic oil viscosity [Pa*s]\n"
213  "\n"
214  << std::endl);
215 
216  if (!HP.IsArg())
217  {
218  /*
219  * Exit quietly if nothing else is provided
220  */
221  throw NoErr(MBDYN_EXCEPT_ARGS);
222  }
223  }
224 
225  if (!HP.IsKeyWord("shaft"))
226  {
227  silent_cerr("hydrodynamic_plain_bearing_with_offset(" << GetLabel() << "): keyword \"shaft\" expected at line " << HP.GetLineData() << std::endl);
228  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
229  }
230 
231  m_pShaft = pDM->ReadNode<const StructNode, Node::STRUCTURAL>(HP);
232 
233  if (!m_pShaft) {
234  silent_cerr("hydrodynamic_plain_bearing_with_offset(" << GetLabel() << "): structural node expected at line " << HP.GetLineData() << std::endl);
235  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
236  }
237 
238  if ( HP.IsKeyWord("offset") )
239  m_o1_R1 = HP.GetPosRel(ReferenceFrame(m_pShaft));
240 
241  if ( !HP.IsKeyWord("bearing"))
242  {
243  silent_cerr("hydrodynamic_plain_bearing_with_offset(" << GetLabel() << "): keyword \"bearing\" expected at line " << HP.GetLineData() << std::endl);
244  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
245  }
246 
247  m_pBearing = pDM->ReadNode<const StructNode, Node::STRUCTURAL>(HP);
248 
249  if (!m_pBearing) {
250  silent_cerr("hydrodynamic_plain_bearing_with_offset(" << GetLabel() << "): structural node expected at line " << HP.GetLineData() << std::endl);
251  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
252  }
253 
254  if ( HP.IsKeyWord("offset") )
255  m_o2_R2 = HP.GetPosRel(ReferenceFrame(m_pBearing));
256 
257  if ( !( HP.IsKeyWord("bearing_width") || HP.IsKeyWord("bearing" "width") ) )
258  {
259  silent_cerr("hydrodynamic_plain_bearing_with_offset(" << GetLabel() << "): keyword \"bearing width\" expected at line " << HP.GetLineData() << std::endl);
260  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
261  }
262 
263  m_bdat.b = HP.GetReal();
264 
265  bool bHaveD = false;
266 
267  if (HP.IsKeyWord("shaft" "diameter"))
268  {
269  m_bdat.d = HP.GetReal();
270  }
271  else
272  {
273  if ( !( HP.IsKeyWord("bearing_diameter") || HP.IsKeyWord("bearing" "diameter") ) )
274  {
275  silent_cerr("hydrodynamic_plain_bearing_with_offset(" << GetLabel() << "): keyword \"bearing diameter\" expected at line " << HP.GetLineData() << std::endl);
276  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
277  }
278 
279  m_bdat.d = HP.GetReal();
280  bHaveD = true;
281  }
282 
283  if ( !( HP.IsKeyWord("relative_clearance") || HP.IsKeyWord("relative" "clearance") ) )
284  {
285  silent_cerr("hydrodynamic_plain_bearing_with_offset(" << GetLabel() << "): keyword \"relative clearance\" expected at line " << HP.GetLineData() << std::endl);
286  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
287  }
288 
289  m_bdat.Psi = HP.GetReal();
290 
291  if (bHaveD)
292  {
293  m_bdat.d *= (1 - m_bdat.Psi);
294  }
295 
296  if ( !( HP.IsKeyWord("oil_viscosity") || HP.IsKeyWord("oil" "viscosity") ) )
297  {
298  silent_cerr("hydrodynamic_plain_bearing_with_offset(" << GetLabel() << "): keyword \"oil viscosity\" expected at line " << HP.GetLineData() << std::endl);
299  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
300  }
301 
302  m_bdat.eta = HP.GetReal();
303 
304  m_InitialAssemblyFactor.Set(( HP.IsKeyWord("initial_assembly_factor") || HP.IsKeyWord("initial" "assembly" "factor") ) ? HP.GetDriveCaller() : new OneDriveCaller());
305 
306  m_iNumGaussPoints = HP.IsKeyWord("number" "of" "gauss" "points") ? HP.GetInt() : 1;
307 
308 #define CASE_GAUSS_POINT_NUM_(num) \
309  case num: \
310  assert(m_iNumGaussPoints == sizeof(s_r##num)/sizeof(s_r##num[0])); \
311  assert(m_iNumGaussPoints == sizeof(s_alpha##num)/sizeof(s_alpha##num[0])); \
312  m_r = s_r##num; \
313  m_alpha = s_alpha##num; \
314  break
315 
316  switch (m_iNumGaussPoints)
317  {
318  CASE_GAUSS_POINT_NUM_(1);
319  CASE_GAUSS_POINT_NUM_(2);
320  CASE_GAUSS_POINT_NUM_(3);
321  CASE_GAUSS_POINT_NUM_(6);
322  default:
323  silent_cerr("hydrodynamic_plain_bearing_with_offset(" << GetLabel()
324  << "): integration rule with " << m_iNumGaussPoints
325  << " gauss points not supported at line "
326  << HP.GetLineData() << std::endl);
327  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
328  }
329 
330  hydrodynamic_plain_bearing_init(m_bdat);
331 
332 #undef CASE_GAUSS_POINT_NUM_
333 
334  ASSERT(m_iNumGaussPoints <= sizeof(m_output) / sizeof(m_output[0]));
335 
336  if (HP.IsKeyWord("output" "points")) {
337  m_iNumOutputPoints = HP.GetInt();
338 
339  if (m_iNumOutputPoints < 0 || m_iNumOutputPoints > m_iNumGaussPoints) {
340  silent_cerr("hydrodynamic_plain_bearing_with_offset(" << GetLabel()
341  << "): number of points for output is not in range [0:" << m_iNumGaussPoints
342  << "] at line " << HP.GetLineData() << std::endl);
343  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
344  }
345 
346  for (integer i = 0; i < m_iNumOutputPoints; ++i) {
347  if (HP.IsKeyWord("gauss" "point")) {
348  integer iGaussPoint = HP.GetInt();
349 
350  if (iGaussPoint < 1 || iGaussPoint > m_iNumGaussPoints) {
351  silent_cerr("hydrodynamic_plain_bearing_with_offset(" << GetLabel()
352  << "): Gauss point index " << iGaussPoint
353  << " is not in range [1:" << m_iNumGaussPoints << "] at line "
354  << HP.GetLineData() << std::endl);
355  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
356  }
357 
358  m_output[i].r = m_r[iGaussPoint - 1];
359  m_output[i].alpha = m_alpha[iGaussPoint - 1];
360  } else if (HP.IsKeyWord("custom")) {
361  if (!HP.IsKeyWord("r")) {
362  silent_cerr("hydrodynamic_plain_bearing_with_offset(" << GetLabel()
363  << "): keyword \"r\" expected at line "
364  << HP.GetLineData() << std::endl);
365  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
366  }
367 
368  m_output[i].r = HP.GetReal();
369 
370  if (fabs(m_output[i].r) > 0.5) {
371  silent_cerr("hydrodynamic_plain_bearing_with_offset(" << GetLabel()
372  << "): abs(r) = " << fabs(m_output[i].r)
373  << " > 0.5 is outside the bearing width at line "
374  << HP.GetLineData() << std::endl);
375  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
376  }
377 
378  if (!HP.IsKeyWord("alpha")) {
379  silent_cerr("hydrodynamic_plain_bearing_with_offset(" << GetLabel()
380  << "): keyword \"alpha\" expected at line "
381  << HP.GetLineData() << std::endl);
382  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
383  }
384 
385  m_output[i].alpha = HP.GetReal();
386 
387  if (m_output[i].alpha < 0 || m_output[i].alpha > 1) {
388  silent_cerr("hydrodynamic_plain_bearing_with_offset(" << GetLabel()
389  << "): alpha = " << m_output[i].alpha << " is not in range [0:1] at line "
390  << HP.GetLineData() << std::endl);
391  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
392  }
393  }
394  }
395  } else {
396  m_iNumOutputPoints = m_iNumGaussPoints;
397 
398  for (integer i = 0; i < m_iNumOutputPoints; ++i) {
399  m_output[i].r = m_r[i];
400  m_output[i].alpha = m_alpha[i];
401  }
402  }
403 
404  if (HP.IsKeyWord("extend" "shaft" "to" "bearing" "center")) {
405  m_lambda = HP.GetYesNoOrBool();
406  }
407 
408  if (HP.IsKeyWord("epsilon" "max")) {
409  m_bdat.eps_max = HP.GetReal();
410  } else {
411  m_bdat.eps_max = 0.999; // According to Butenschoen this model is valid until epsilon = 0.999
412  }
413 
414  SetOutputFlag(pDM->fReadOutput(HP, Elem::LOADABLE));
415 
416  std::ostream& out = pDM->GetLogFile();
417 
418  out << "hydrodynamic_plain_bearing_with_offset: "
419  << uLabel << " "
420  << m_pShaft->GetLabel() << " "
421  << m_o1_R1 << " "
422  << m_pBearing->GetLabel() << " "
423  << m_o2_R2 << " "
424  << m_bdat.b << " "
425  << m_bdat.d << " "
426  << m_bdat.Psi << " "
427  << m_bdat.eta << " ";
428 
429  out << m_iNumGaussPoints << " ";
430 
431  for (integer i = 0; i < m_iNumGaussPoints; ++i)
432  {
433  out << m_r[i] << " " << m_alpha[i] << " ";
434  }
435 
436  out << m_iNumOutputPoints << " ";
437 
438  for (integer i = 0; i < m_iNumOutputPoints; ++i)
439  {
440  out << m_output[i].r << " " << m_output[i].alpha << " ";
441  }
442 
443  out << std::endl;
444 }
445 
446 HydrodynamicPlainBearing::~HydrodynamicPlainBearing(void)
447 {
448  // destroy private dataman
449 }
450 
451 void
452 HydrodynamicPlainBearing::Output(OutputHandler& OH) const
453 {
454  if ( bToBeOutput() && OH.UseText(OutputHandler::LOADABLE) )
455  {
456  std::ostream& os = OH.Loadable();
457 
458  os << std::setw(8) << GetLabel() << ' '; // 0
459 
460  for (integer i = 0; i < m_iNumOutputPoints; ++i)
461  {
462  doublereal eps, eps_dot, delta, SoD, SoV, my, beta;
463 
464  Vec3 e_R2, e_dot_R2;
465  doublereal omega_proj[2];
466  Vec3 F2_R2, M2_R2;
467  Vec3 F2_I, M2_I, F1_I, M1_I;
468 
469  ComputeResidual(e_R2, e_dot_R2, omega_proj, F2_R2, M2_R2, F2_I, M2_I, F1_I, M1_I, eps, eps_dot, delta, SoD, SoV, my, beta, m_output[i].r, m_output[i].alpha);
470 
471  // output the current state: the column layout is as follows
472  // (all quantities are refered to the reference frame of the bearing node)
473 
474  // 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
475  // e(1) e(2) e_dot(1) e_dot(2) omega(1) omega(2) eps eps_dot delta Fx Fy Mz SoD SoV my beta
476 
477  // 0 label of the element
478  // 1 absolute eccentricity in of the shaft in x direction
479  // 2 absolute eccentricity in of the shaft in y direction
480  // 3 difference of the absolute velocity between the shaft and the bearing in x direction
481  // 4 difference of the absolute velocity between the shaft and the bearing in y direction
482  // 5 absolute angular velocity of the shaft around the z axis
483  // 6 absolute angular velocity of the bearing around the z axis
484  // 7 relative eccentricity of the shaft
485  // 8 time derivative of the relative eccentricity
486  // 9 angle of minimum clearance
487  // 10 force at the bearing in x direction
488  // 11 force at the bearing in y direction
489  // 12 torque at the bearing around the z axis
490  // 13 Sommerfeld number for rotation
491  // 14 Sommerfeld number for displacement
492  // 15 friction coefficient
493  // 16 angle between minimum clearance and force of rotation
494 
495  os << e_R2(1) << ' ' // 1
496  << e_R2(2) << ' ' // 2
497  << e_dot_R2(1) << ' ' // 3
498  << e_dot_R2(2) << ' ' // 4
499  << omega_proj[0] << ' ' // 5
500  << omega_proj[1] << ' ' // 6
501  << eps << ' ' // 7
502  << eps_dot << ' ' // 8
503  << delta << ' ' // 9
504  << F2_R2(1) << ' ' // 10
505  << F2_R2(2) << ' ' // 11
506  << M2_R2(3) << ' ' // 12
507  << SoD << ' ' // 13
508  << SoV << ' ' // 14
509  << my << ' ' // 15
510  << beta << ' '; // 16
511  }
512 
513  os << std::endl;
514  }
515 }
516 
517 void
518 HydrodynamicPlainBearing::WorkSpaceDim(integer* piNumRows, integer* piNumCols) const
519 {
520  *piNumRows = 12;
521  *piNumCols = 12;
522 }
523 
524 VariableSubMatrixHandler&
525 HydrodynamicPlainBearing::AssJac(VariableSubMatrixHandler& WorkMatV,
526  doublereal dCoef,
527  const VectorHandler& XCurr,
528  const VectorHandler& XPrimeCurr)
529 {
530  integer iNumRows = 0;
531  integer iNumCols = 0;
532 
533  WorkSpaceDim(&iNumRows, &iNumCols);
534 
535  FullSubMatrixHandler& WorkMat = WorkMatV.SetFull();
536 
537  WorkMat.ResizeReset(iNumRows, iNumCols);
538 
539  integer intShaftPositionIndex = m_pShaft->iGetFirstPositionIndex();
540  integer intShaftMomentumIndex = m_pShaft->iGetFirstMomentumIndex();
541  integer intBearingPositionIndex = m_pBearing->iGetFirstPositionIndex();
542  integer intBearingMomentumIndex = m_pBearing->iGetFirstMomentumIndex();
543 
544  for (int iCnt = 1; iCnt <= 6; iCnt++)
545  {
546  WorkMat.PutRowIndex(iCnt, intShaftMomentumIndex + iCnt);
547  WorkMat.PutColIndex(iCnt, intShaftPositionIndex + iCnt);
548  WorkMat.PutRowIndex(iCnt+6,intBearingMomentumIndex + iCnt);
549  WorkMat.PutColIndex(iCnt+6,intBearingPositionIndex + iCnt);
550  }
551 
552  const Vec3& X1 = m_pShaft->GetXCurr();
553  const Vec3& X2 = m_pBearing->GetXCurr();
554  const Vec3& X1_dot = m_pShaft->GetVCurr();
555  const Vec3& X2_dot = m_pBearing->GetVCurr();
556  const Mat3x3& R1 = m_pShaft->GetRCurr();
557  const Mat3x3& R1_0 = m_pShaft->GetRRef();
558  const Mat3x3& R2 = m_pBearing->GetRCurr();
559  const Mat3x3& R2_0 = m_pBearing->GetRRef();
560  const Vec3& omega1 = m_pShaft->GetWCurr();
561  const Vec3& omega1_ref = m_pShaft->GetWRef();
562  const Vec3& omega2 = m_pBearing->GetWCurr();
563  const Vec3& omega2_ref = m_pBearing->GetWRef();
564 
565  for (integer i = 0; i < m_iNumGaussPoints; ++i)
566  {
567  Vec3 o1_R1 = m_o1_R1;
568  o1_R1(3) += m_r[i] * m_bdat.b;
569  Vec3 o2_R2 = m_o2_R2;
570  o2_R2(3) += m_r[i] * m_bdat.b;
571  const Vec3 v_R2 = R2.MulTV( X1 - X2 + R1 * o1_R1 ) - o2_R2;
572  const Vec3 d1_R2 = R2.MulTV(R1.GetCol(3));
573 
574  const doublereal lambda = m_lambda ? -v_R2(3) / d1_R2(3) : 0.;
575 
576  Vec3 e_R2 = v_R2 + d1_R2 * lambda;
577 
578  // FIXME: nan values if e_R2(1) == 0 && e_R2(2) == 0
579  if ( e_R2(1) == 0.0 && e_R2(2) == 0.0 )
580  e_R2(1) = std::numeric_limits<doublereal>::epsilon() * m_bdat.d * m_bdat.Psi / 2.;
581 
582  const Vec3 v_dot_R2 = R2.MulTV( X1_dot - X2_dot + omega2.Cross(X2 - X1) + ( omega1 - omega2 ).Cross( R1 * o1_R1 ) );
583  const Vec3 d1_dot_R2 = R2.MulTV( ( omega1 - omega2 ).Cross( R1.GetCol(3) ) );
584  const doublereal lambda_dot = m_lambda ? -v_dot_R2(3) / d1_R2(3) + v_R2(3) / std::pow(d1_R2(3),2) * d1_dot_R2(3) : 0.;
585 
586  Vec3 e_dot_R2 = R2.MulTV( X1_dot - X2_dot + omega1.Cross( R1 * o1_R1 ) - omega2.Cross( R2 * o2_R2 )
587  + R1.GetCol(3) * lambda_dot + omega1.Cross(R1.GetCol(3)) * lambda);
588 
589  // FIXME: nan values if e_dot_R2(1) == 0 && e_dot_R2(2) == 0
590  if ( e_dot_R2(1) == 0.0 && e_dot_R2(2) == 0.0 )
591  e_dot_R2(1) = std::numeric_limits<doublereal>::epsilon() * m_bdat.d * m_bdat.Psi / 2.;
592 
593  const Vec3 lambda_d1_R1 = Vec3(0.,0.,lambda);
594  const Vec3 l1_R1 = o1_R1 + lambda_d1_R1;
595  const Vec3 l1_I = R1 * l1_R1;
596  const Vec3 l2_R2 = o2_R2 + e_R2;
597  const doublereal omega_proj[2] = { R2.GetCol(3).Dot(omega1),
598  R2.GetCol(3).Dot(omega2) };
599 
600  const doublereal e_[2] = { e_R2(1), e_R2(2) }, e_dot_[2] = { e_dot_R2(1), e_dot_R2(2) };
601 
602  // | 0 1 2 3 4 5 |
603  static const doublereal ed[2][NBDIRSMAX] = { { 1., 0., 0., 0., 0., 0. }, // 0 | inner derivative of the eccentricity of the shaft in direction 1 of the reference frame of the bearing (R2)
604  { 0., 1., 0., 0., 0., 0. } }; // 1 | inner derivative of the eccentricity of the shaft in direction 2 of the reference frame of the bearing (R2)
605  static const doublereal e_dotd[2][NBDIRSMAX] = { { 0., 0., 1., 0., 0., 0. }, // 0 | inner derivative of relative velocity of the shaft in direction 1 of the reference frame of the bearing (R2)
606  { 0., 0., 0., 1., 0., 0. } }; // 1 | inner derivative of relative velocity of the shaft in direction 2 of the reference frame of the bearing (R2)
607  static const doublereal omega_projd[2][NBDIRSMAX] = { { 0., 0., 0., 0., 1., 0. }, // 0 | inner derivative of the angular velocity of the shaft in direction 3 of the reference frame of the bearing (R2)
608  { 0., 0., 0., 0., 0., 1. } }; // 1 | inner derivative of the angular velocity of the bearing in direction 3 of the reference frame of the bearing (R2)
609 
610  doublereal k[3]; // force vector at the bearing (not used for the evaluation of the jacobian matrix)
611  doublereal kd[3][NBDIRSMAX]; // variation of the force vector at the bearing with respect to the eccentricity of the shaft e and the relative velocity of the shaft
612 
613  // kd = { { diff(k[0],e[0]), diff(k[0],e[1]), diff(k[0],e_dot[0]), diff(k[0],e_dot[1]), diff(k[0],omega_proj[0]), diff(k[0],omega_proj[1]) },
614  // { diff(k[1],e[0]), diff(k[1],e[1]), diff(k[1],e_dot[0]), diff(k[1],e_dot[1]), diff(k[1],omega_proj[0]), diff(k[1],omega_proj[1]) },
615  // { diff(k[2],e[0]), diff(k[2],e[1]), diff(k[2],e_dot[0]), diff(k[2],e_dot[1]), diff(k[2],omega_proj[0]), diff(k[2],omega_proj[1]) } };
616 
617  doublereal eps, eps_dot, delta, SoD, SoV, my, beta;
618 
619  hydrodynamic_plain_bearing_force_dv(m_bdat, omega_proj, omega_projd, e_, ed, e_dot_, e_dotd, k, kd, eps, eps_dot, delta, SoD, SoV, my, beta, NBDIRSMAX);
620 
621  Vec3 F2_R2; // F2^(R2) = f(e^(R2), diff(e^(R2),t), omega1_proj, omega2_proj)
622 
623  F2_R2(1) = k[0];
624  F2_R2(2) = k[1];
625  F2_R2(3) = 0.;
626 
627  Mat3x3 dF2_R2_de_R2; // diff(F2^(R2), e^(R2))
628 
629  dF2_R2_de_R2(1,1) = kd[0][0]; dF2_R2_de_R2(1,2) = kd[0][1]; dF2_R2_de_R2(1,3) = 0.;
630  dF2_R2_de_R2(2,1) = kd[1][0]; dF2_R2_de_R2(2,2) = kd[1][1]; dF2_R2_de_R2(2,3) = 0.;
631  dF2_R2_de_R2(3,1) = 0.; dF2_R2_de_R2(3,2) = 0.; dF2_R2_de_R2(3,3) = 0.;
632 
633  Mat3x3 dF2_R2_de_dot_R2; // diff(F2^(R2), diff(e^(R2),t))
634 
635  dF2_R2_de_dot_R2(1,1) = kd[0][2]; dF2_R2_de_dot_R2(1,2) = kd[0][3]; dF2_R2_de_dot_R2(1,3) = 0.;
636  dF2_R2_de_dot_R2(2,1) = kd[1][2]; dF2_R2_de_dot_R2(2,2) = kd[1][3]; dF2_R2_de_dot_R2(2,3) = 0.;
637  dF2_R2_de_dot_R2(3,1) = 0.; dF2_R2_de_dot_R2(3,2) = 0; dF2_R2_de_dot_R2(3,3) = 0.;
638 
639  Vec3 dF2_R2_domega1_proj; // diff(F2^(R2), omega1_proj)
640 
641  dF2_R2_domega1_proj(1) = kd[0][4];
642  dF2_R2_domega1_proj(2) = kd[1][4];
643  dF2_R2_domega1_proj(3) = 0.;
644 
645  Vec3 dF2_R2_domega2_proj; // diff(F2^(R2), omega2_proj)
646 
647  dF2_R2_domega2_proj(1) = kd[0][5];
648  dF2_R2_domega2_proj(2) = kd[1][5];
649  dF2_R2_domega2_proj(3) = 0.;
650 
651  Vec3 M2_R2; // M2^(R2) = f(e^(R2), diff(e^(R2),t), omega1_proj, omega2_proj)
652 
653  M2_R2(1) = 0.;
654  M2_R2(2) = 0.;
655  M2_R2(3) = k[2];
656 
657  Mat3x3 dM2_R2_de_R2; // diff(M2^(R2), e^(R2)
658 
659  dM2_R2_de_R2(1,1) = 0.; dM2_R2_de_R2(1,2) = 0.; dM2_R2_de_R2(1,3) = 0.;
660  dM2_R2_de_R2(2,1) = 0.; dM2_R2_de_R2(2,2) = 0.; dM2_R2_de_R2(2,3) = 0.;
661  dM2_R2_de_R2(3,1) = kd[2][0]; dM2_R2_de_R2(3,2) = kd[2][1]; dM2_R2_de_R2(3,3) = 0.;
662 
663  Mat3x3 dM2_R2_de_dot_R2; // diff(M2^(R2), diff(e^(R2),t)
664 
665  dM2_R2_de_dot_R2(1,1) = 0.; dM2_R2_de_dot_R2(1,2) = 0.; dM2_R2_de_dot_R2(1,3) = 0.;
666  dM2_R2_de_dot_R2(2,1) = 0.; dM2_R2_de_dot_R2(2,2) = 0.; dM2_R2_de_dot_R2(2,3) = 0.;
667  dM2_R2_de_dot_R2(3,1) = kd[2][2]; dM2_R2_de_dot_R2(3,2) = kd[2][3]; dM2_R2_de_dot_R2(3,3) = 0.;
668 
669  Vec3 dM2_R2_domega1_proj; // diff(M2^(R2), omega1_proj)
670 
671  dM2_R2_domega1_proj(1) = 0.;
672  dM2_R2_domega1_proj(2) = 0.;
673  dM2_R2_domega1_proj(3) = kd[2][4];
674 
675  Vec3 dM2_R2_domega2_proj; // diff(M2^(R2), omega2_proj)
676 
677  dM2_R2_domega2_proj(1) = 0.;
678  dM2_R2_domega2_proj(2) = 0.;
679  dM2_R2_domega2_proj(3) = kd[2][5];
680 
681  const doublereal alpha = m_alpha[i] * m_InitialAssemblyFactor.dGet();
682 
683  F2_R2 *= alpha;
684  dF2_R2_de_R2 *= alpha;
685  dF2_R2_de_dot_R2 *= alpha;
686  dF2_R2_domega1_proj *= alpha;
687  dF2_R2_domega2_proj *= alpha;
688 
689  M2_R2 *= alpha;
690  dM2_R2_de_R2 *= alpha;
691  dM2_R2_de_dot_R2 *= alpha;
692  dM2_R2_domega1_proj *= alpha;
693  dM2_R2_domega2_proj *= alpha;
694 
695  const Mat3x3 R2_T = R2.Transpose();
696 
697  const Vec3 F2_I = R2 * F2_R2;
698  const Vec3 M2_I = R2 * ( l2_R2.Cross( F2_R2 ) + M2_R2 );
699  const Vec3 F1_I = -F2_I;
700  const Vec3 M1_I = -l1_I.Cross( F2_I ) - R2 * M2_R2;
701 
702  const Mat3x3 dv_dot_R2_dX1 = -R2.MulTM(Mat3x3(MatCross, omega2)); // diff(diff(v^(R2),t),X1)
703  const Mat3x3& dv_R2_dX1 = R2_T; // diff(v^(R2),X1)
704  const Vec3 dlambda_dot_dX1 = m_lambda ? -dv_dot_R2_dX1.GetRow(3) / d1_R2(3) + dv_R2_dX1.GetRow(3) / std::pow(d1_R2(3),2) * d1_dot_R2(3) : Zero3; // diff(diff(lambda,t),X1)
705 
706  const Vec3 dlambda_dX1 = m_lambda ? -dv_R2_dX1.GetRow(3) / d1_R2(3) : Zero3;
707 
708  const Mat3x3 de_dot_R2_dX1 = R2.MulTM( R1.GetCol(3).Tens(dlambda_dot_dX1) + omega1.Cross( R1.GetCol(3) ).Tens(dlambda_dX1) );
709 
710  const Mat3x3 dv_dot_R2_domega1 = -R2.MulTM(Mat3x3(MatCross, R1 * o1_R1));
711  const Mat3x3 domega1_dg1 = -Mat3x3(MatCross, omega1_ref);
712  const Mat3x3 dv_R2_dg1 = -R2.MulTM( Mat3x3(MatCross, R1_0 * o1_R1) );
713  const Mat3x3 dv_dot_R2_dg1 = dv_dot_R2_domega1 * domega1_dg1 - R2.MulTM( ( omega1 - omega2 ).Cross(Mat3x3(MatCross, R1_0 * o1_R1)) );
714  const Mat3x3 dd1_R2_dg1 = -R2.MulTM(Mat3x3(MatCross, R1_0.GetCol(3)));
715  const Mat3x3 dd1_dot_R2_dg1 = R2.MulTM( R1.GetCol(3).Cross(Mat3x3(MatCross, omega1_ref)) - ( omega1 - omega2 ).Cross(Mat3x3(MatCross, R1_0.GetCol(3))) );
716  const Vec3 dlambda_dot_dg1 = m_lambda ? -dv_dot_R2_dg1.GetRow(3) / d1_R2(3)
717  + dd1_R2_dg1.GetRow(3) * (v_dot_R2(3) / std::pow(d1_R2(3), 2))
718  + dv_R2_dg1.GetRow(3) / std::pow(d1_R2(3), 2) * d1_dot_R2(3)
719  - dd1_R2_dg1.GetRow(3) * (2. * v_R2(3) / std::pow(d1_R2(3), 3) * d1_dot_R2(3))
720  + dd1_dot_R2_dg1.GetRow(3) * (v_R2(3) / std::pow(d1_R2(3), 2)) : Zero3;
721  const Vec3 dlambda_dg1 = m_lambda ? -dv_R2_dg1.GetRow(3) / d1_R2(3) + dd1_R2_dg1.GetRow(3) * (v_R2(3) / std::pow( d1_R2(3), 2 )) : Zero3;
722 
723  const Mat3x3 de_dot_R2_dg1 = R2.MulTM( -omega1.Cross(Mat3x3( MatCross, R1_0 * o1_R1 )) + ( R1 * o1_R1 ).Cross(Mat3x3(MatCross, omega1_ref))
724  + R1.GetCol(3).Tens(dlambda_dot_dg1) - Mat3x3( MatCross, R1_0.GetCol(3) ) * lambda_dot
725  + omega1.Cross( R1.GetCol(3) ).Tens(dlambda_dg1) + R1.GetCol(3).Cross(Mat3x3(MatCross, omega1_ref)) * lambda
726  - omega1.Cross(Mat3x3( MatCross, R1_0.GetCol(3) )) * lambda );
727  const Mat3x3 dv_dot_R2_dX2 = R2.MulTM(Mat3x3(MatCross, omega2));
728  const Mat3x3 dv_R2_dX2 = -R2_T;
729  const Vec3 dlambda_dot_dX2 = m_lambda ? -dv_dot_R2_dX2.GetRow(3) / d1_R2(3) + dv_R2_dX2.GetRow(3) / std::pow(d1_R2(3),2) * d1_dot_R2(3) : Zero3;
730 
731  const Vec3 dlambda_dX2 = m_lambda ? -dv_R2_dX2.GetRow(3) / d1_R2(3) : Zero3;
732  const Mat3x3 de_dot_R2_dX2 = R2.MulTM( R1.GetCol(3).Tens(dlambda_dot_dX2) + omega1.Cross( R1.GetCol(3) ).Tens(dlambda_dX2) );
733  const Mat3x3 dv_R2_dg2 = R2_0.MulTM( Mat3x3( MatCross, X1 - X2 + R1 * o1_R1 ) );
734  const Mat3x3 dd1_R2_dg2 = R2_0.MulTM( Mat3x3( MatCross, R1 * o1_R1 ) );
735  const Mat3x3 dv_dot_R2_dg2 = R2_0.MulTM( Mat3x3( MatCross, X1_dot - X2_dot + omega2.Cross( X2 - X1 ) + ( omega1 - omega2 ).Cross( R1 * o1_R1 ) ) )
736  + R2.MulTM( ( X2 - X1 - R1 * o1_R1 ).Cross(Mat3x3(MatCross, omega2_ref)) );
737  const Vec3 dlambda_dg2 = m_lambda ? -dv_R2_dg2.GetRow(3) / d1_R2(3) + dd1_R2_dg2.GetRow(3) * (v_R2(3) / std::pow(d1_R2(3),2)) : Zero3;
738 
739  const Mat3x3 dd1_dot_R2_dg2 = R2_0.MulTM( Mat3x3( MatCross, ( omega1 - omega2 ).Cross( R1.GetCol(3) ) ) )
740  - R2.MulTM( Mat3x3( MatCross, R1.GetCol(3) ) * Mat3x3(MatCross, omega2_ref) );
741  const Vec3 dlambda_dot_dg2 = m_lambda ? -dv_dot_R2_dg2.GetRow(3) / d1_R2(3) + dd1_R2_dg2.GetRow(3) * (v_dot_R2(3) / std::pow(d1_R2(3), 2))
742  + dv_R2_dg2.GetRow(3) / std::pow(d1_R2(3),2) * d1_dot_R2(3)
743  - dd1_R2_dg2.GetRow(3) * (2. * v_R2(3) / std::pow( d1_R2(3), 3) * d1_dot_R2(3))
744  + dd1_dot_R2_dg2.GetRow(3) * (v_R2(3) / std::pow(d1_R2(3),2)) : Zero3;
745 
746  const Mat3x3 de_dot_R2_dg2 = R2_0.MulTM( Mat3x3( MatCross, X1_dot - X2_dot + omega1.Cross(R1 * o1_R1) - omega2.Cross( R2 * o2_R2 )
747  + R1.GetCol(3) * lambda_dot + omega1.Cross(R1.GetCol(3)) * lambda))
748  + R2.MulTM( -( R2 * o2_R2 ).Cross(Mat3x3(MatCross, omega2_ref)) + omega2.Cross(Mat3x3( MatCross, R2_0 * o2_R2 ))
749  + R1.GetCol(3).Tens(dlambda_dot_dg2) + omega1.Cross(R1.GetCol(3)).Tens(dlambda_dg2));
750  const Mat3x3& dv_dot_R2_dX1_dot = R2_T;
751  const Vec3 dlambda_dot_dX1_dot = m_lambda ? -dv_dot_R2_dX1_dot.GetRow(3) / d1_R2(3) : Zero3;
752 
753  const Mat3x3 de_dot_R2_dX1_dot = R2.MulTM( Eye3 + R1.GetCol(3).Tens(dlambda_dot_dX1_dot));
754 
755  const Mat3x3 dv_dot_R2_dg1_dot = -R2.MulTM( Mat3x3( MatCross, R1 * o1_R1 ) );
756  const Mat3x3 dd1_dot_R2_dg1_dot = -R2.MulTM( Mat3x3( MatCross, R1.GetCol(3) ) );
757  const Vec3 dlambda_dot_dg1_dot = m_lambda ? -dv_dot_R2_dg1_dot.GetRow(3) / d1_R2(3) + dd1_dot_R2_dg1_dot.GetRow(3) * (v_R2(3) / std::pow(d1_R2(3),2)) : Zero3;
758 
759  const Mat3x3 de_dot_R2_dg1_dot = R2.MulTM( -Mat3x3( MatCross, R1 * o1_R1 ) + R1.GetCol(3).Tens(dlambda_dot_dg1_dot) - Mat3x3(MatCross, R1.GetCol(3) * lambda) );
760 
761  const Mat3x3 dv_dot_dX2_dot = -R2_T;
762  const Vec3 dlambda_dot_dX2_dot = m_lambda ? -dv_dot_dX2_dot.GetRow(3) / d1_R2(3) : Zero3;
763 
764  const Mat3x3 de_dot_R2_dX2_dot = R2.MulTM( -Eye3 + R1.GetCol(3).Tens(dlambda_dot_dX2_dot));
765 
766  const Mat3x3 dv_dot_R2_dg2_dot = R2.MulTM( Mat3x3( MatCross, R1 * o1_R1 + X1 - X2 ));
767  const Mat3x3 dd1_dot_R2_dg2_dot = R2.MulTM(Mat3x3(MatCross, R1.GetCol(3)));
768  const Vec3 dlambda_dot_dg2_dot = m_lambda ? -dv_dot_R2_dg2_dot.GetRow(3) / d1_R2(3) + dd1_dot_R2_dg2_dot.GetRow(3) * (v_R2(3) / std::pow(d1_R2(3),2)) : Zero3;
769 
770  const Mat3x3 de_dot_R2_dg2_dot = R2.MulTM( Mat3x3( MatCross, R2 * o2_R2 ) + R1.GetCol(3).Tens( dlambda_dot_dg2_dot) );
771 
772  const Mat3x3 de_R2_dX1 = dv_R2_dX1 + d1_R2.Tens( dlambda_dX1 );
773 
774  const Mat3x3 de_R2_dg1 = dv_R2_dg1 + d1_R2.Tens(dlambda_dg1) + dd1_R2_dg1 * lambda;
775 
776  const Mat3x3 de_R2_dX2 = dv_R2_dX2 + d1_R2.Tens(dlambda_dX2);
777 
778  const Mat3x3 de_R2_dg2 = dv_R2_dg2 + d1_R2.Tens( dlambda_dg2 ) + dd1_R2_dg2 * lambda;
779 
780  const Vec3 domega1_proj_dg1 = omega1_ref.Cross( R2.GetCol(3) );
781  const Vec3 domega1_proj_dg1_dot = R2.GetCol(3);
782  const Vec3 domega1_proj_dg2 = -omega1.Cross( R2_0.GetCol(3) );
783  const Vec3 domega2_proj_dg2 = -omega2.Cross(R2_0.GetCol(3)) + omega2_ref.Cross(R2.GetCol(3));
784  const Vec3 domega2_proj_dg2_dot = R2.GetCol(3);
785 
786  const Mat3x3 dF2_R2_dX1 = dF2_R2_de_R2 * de_R2_dX1 + dF2_R2_de_dot_R2 * de_dot_R2_dX1;
787 
788  const Mat3x3 dF2_R2_dg1 = dF2_R2_de_R2 * de_R2_dg1 + dF2_R2_de_dot_R2 * de_dot_R2_dg1
789  + dF2_R2_domega1_proj.Tens(domega1_proj_dg1);
790  const Mat3x3 dF2_R2_dX2 = dF2_R2_de_R2 * de_R2_dX2 + dF2_R2_de_dot_R2 * de_dot_R2_dX2;
791  const Mat3x3 dF2_R2_dg2 = dF2_R2_de_R2 * de_R2_dg2 + dF2_R2_de_dot_R2 * de_dot_R2_dg2
792  + dF2_R2_domega1_proj.Tens(domega1_proj_dg2) + dF2_R2_domega2_proj.Tens(domega2_proj_dg2);
793 
794  const Mat3x3 dF2_I_dX1 = R2 * dF2_R2_dX1;
795  const Mat3x3 dF2_I_dg1 = R2 * dF2_R2_dg1;
796  const Mat3x3 dF2_I_dX2 = R2 * dF2_R2_dX2;
797  const Mat3x3 dF2_I_dg2 = -Mat3x3( MatCross, R2_0 * F2_R2 ) + R2 * dF2_R2_dg2;
798 
799  const Mat3x3 dF1_I_dX1 = -dF2_I_dX1;
800  const Mat3x3 dF1_I_dg1 = -dF2_I_dg1;
801  const Mat3x3 dF1_I_dX2 = -dF2_I_dX2;
802  const Mat3x3 dF1_I_dg2 = -dF2_I_dg2;
803 
804  const Mat3x3 dM2_R2_dX1 = dM2_R2_de_R2 * de_R2_dX1 + dM2_R2_de_dot_R2 * de_dot_R2_dX1;
805  const Mat3x3 dM2_I_dX1 = R2 * ( -F2_R2.Cross(de_R2_dX1) + l2_R2.Cross(dF2_R2_dX1) + dM2_R2_dX1 );
806 
807  const Mat3x3 dM2_R2_dg1 = dM2_R2_de_R2 * de_R2_dg1 + dM2_R2_de_dot_R2 * de_dot_R2_dg1 + dM2_R2_domega1_proj.Tens(domega1_proj_dg1);
808  const Mat3x3 dM2_I_dg1 = R2 * ( -F2_R2.Cross(de_R2_dg1) + l2_R2.Cross(dF2_R2_dg1) + dM2_R2_dg1 );
809 
810  const Mat3x3 dM2_R2_dX2 = dM2_R2_de_R2 * de_R2_dX2 + dM2_R2_de_dot_R2 * de_dot_R2_dX2;
811  const Mat3x3 dM2_I_dX2 = R2 * ( -F2_R2.Cross( de_R2_dX2 ) + l2_R2.Cross(dF2_R2_dX2) + dM2_R2_dX2 );
812 
813  const Mat3x3 dM2_R2_dg2 = dM2_R2_de_R2 * de_R2_dg2 + dM2_R2_de_dot_R2 * de_dot_R2_dg2
814  + dM2_R2_domega1_proj.Tens(domega1_proj_dg2) + dM2_R2_domega2_proj.Tens(domega2_proj_dg2);
815  const Mat3x3 dM2_I_dg2 = -Mat3x3( MatCross, R2_0 * ( l2_R2.Cross(F2_R2) + M2_R2 ) )
816  + R2 * ( -F2_R2.Cross(de_R2_dg2) + l2_R2.Cross(dF2_R2_dg2) + dM2_R2_dg2 );
817 
818  const Mat3x3 dM1_I_dX1 = ( R2 * F2_R2 ).Cross( R1.GetCol(3) ).Tens( dlambda_dX1 )
819  - l1_I.Cross( R2 * dF2_R2_dX1 ) - R2 * dM2_R2_dX1;
820 
821  const Mat3x3 dM1_I_dg1 = ( R2 * F2_R2 ).Cross( R1.GetCol(3).Tens(dlambda_dg1) - Mat3x3( MatCross, R1_0 * l1_R1 ) )
822  - l1_I.Cross( R2 * dF2_R2_dg1 ) - R2 * dM2_R2_dg1;
823 
824  const Mat3x3 dM1_I_dX2 = ( R2 * F2_R2 ).Cross( R1.GetCol(3).Tens(dlambda_dX2) )
825  - l1_I.Cross( R2 * dF2_R2_dX2 ) - R2 * dM2_R2_dX2;
826 
827  const Mat3x3 dM1_I_dg2 = ( R2 * F2_R2 ).Cross( R1.GetCol(3).Tens(dlambda_dg2) )
828  + l1_I.Cross( Mat3x3( MatCross, R2_0 * F2_R2 ) - R2 * dF2_R2_dg2 )
829  + Mat3x3( MatCross, R2_0 * M2_R2 ) - R2 * dM2_R2_dg2;
830 
831  const Mat3x3 dF2_R2_dX1_dot = dF2_R2_de_dot_R2 * de_dot_R2_dX1_dot;
832  const Mat3x3 dF2_I_dX1_dot = R2 * dF2_R2_dX1_dot;
833 
834  const Mat3x3 dF2_R2_dg1_dot = dF2_R2_de_dot_R2 * de_dot_R2_dg1_dot + dF2_R2_domega1_proj.Tens( domega1_proj_dg1_dot );
835  const Mat3x3 dF2_I_dg1_dot = R2 * dF2_R2_dg1_dot;
836 
837  const Mat3x3 dF2_R2_dX2_dot = dF2_R2_de_dot_R2 * de_dot_R2_dX2_dot;
838  const Mat3x3 dF2_I_dX2_dot = R2 * dF2_R2_dX2_dot;
839 
840  const Mat3x3 dF2_R2_dg2_dot = dF2_R2_de_dot_R2 * de_dot_R2_dg2_dot + dF2_R2_domega2_proj.Tens( domega2_proj_dg2_dot );
841  const Mat3x3 dF2_I_dg2_dot = R2 * dF2_R2_dg2_dot;
842 
843  const Mat3x3 dF1_I_dX1_dot = -dF2_I_dX1_dot;
844  const Mat3x3 dF1_I_dg1_dot = -dF2_I_dg1_dot;
845  const Mat3x3 dF1_I_dX2_dot = -dF2_I_dX2_dot;
846  const Mat3x3 dF1_I_dg2_dot = -dF2_I_dg2_dot;
847 
848  const Mat3x3 dM2_R2_dX1_dot = dM2_R2_de_dot_R2 * de_dot_R2_dX1_dot;
849  const Mat3x3 dM2_I_dX1_dot = R2 * ( l2_R2.Cross( dF2_R2_dX1_dot) + dM2_R2_dX1_dot );
850 
851  const Mat3x3 dM2_R2_dg1_dot = dM2_R2_de_dot_R2 * de_dot_R2_dg1_dot + dM2_R2_domega1_proj.Tens( domega1_proj_dg1_dot );
852  const Mat3x3 dM2_I_dg1_dot = R2 * ( l2_R2.Cross( dF2_R2_dg1_dot ) + dM2_R2_dg1_dot );
853 
854  const Mat3x3 dM2_R2_dX2_dot = dM2_R2_de_dot_R2 * de_dot_R2_dX2_dot;
855  const Mat3x3 dM2_I_dX2_dot = R2 * ( l2_R2.Cross(dF2_R2_dX2_dot) + dM2_R2_dX2_dot );
856 
857  const Mat3x3 dM2_R2_dg2_dot = dM2_R2_de_dot_R2 * de_dot_R2_dg2_dot + dM2_R2_domega2_proj.Tens(domega2_proj_dg2_dot);
858  const Mat3x3 dM2_I_dg2_dot = R2 * ( l2_R2.Cross( dF2_R2_dg2_dot ) + dM2_R2_dg2_dot );
859 
860  const Mat3x3 dM1_I_dX1_dot = -l1_I.Cross( R2 * dF2_R2_dX1_dot ) - R2 * dM2_R2_dX1_dot;
861  const Mat3x3 dM1_I_dg1_dot = -l1_I.Cross( R2 * dF2_R2_dg1_dot ) - R2 * dM2_R2_dg1_dot;
862  const Mat3x3 dM1_I_dX2_dot = -l1_I.Cross( R2 * dF2_R2_dX2_dot ) - R2 * dM2_R2_dX2_dot;
863  const Mat3x3 dM1_I_dg2_dot = -l1_I.Cross( R2 * dF2_R2_dg2_dot ) - R2 * dM2_R2_dg2_dot;
864 
865  /*
866  * 1, 4, 7, 10 1, 4, 7, 10
867  * | dF1/dX1_dot, dF1/dg1_dot, dF1/dX2_dot, dF1/dg2_dot | | dF1/dX1, dF1/dg1, dF1/dX2, dF1/dg2 | 1
868  * | dM1/dX1_dot, dM1/dg1_dot, dM1/dX2_dot, dM1/dg2_dot | | dM1/dX1, dM1/dg1, dM1/dX2, dM1/dg2 | 4
869  * Jac = -| | -dCoef * | |
870  * | dF2/dX1_dot, dF2/dg1_dot, dF2/dX2_dot, dF2/dg2_dot | | dF2/dX1, dF2/dg1, dF2/dX2, dF2/dg2 | 7
871  * | dM2/dX1_dot, dM2/dg1_dot, dM2/dX2_dot, dM2/dg2_dot | | dM2/dX1, dM2/dg1, dM2/dX2, dM2/dg2 | 10
872  */
873  WorkMat.Sub( 1, 1, dF1_I_dX1_dot + dF1_I_dX1 * dCoef );
874  WorkMat.Sub( 1, 4, dF1_I_dg1_dot + dF1_I_dg1 * dCoef );
875  WorkMat.Sub( 1, 7, dF1_I_dX2_dot + dF1_I_dX2 * dCoef );
876  WorkMat.Sub( 1, 10, dF1_I_dg2_dot + dF1_I_dg2 * dCoef );
877 
878  WorkMat.Sub( 4, 1, dM1_I_dX1_dot + dM1_I_dX1 * dCoef );
879  WorkMat.Sub( 4, 4, dM1_I_dg1_dot + dM1_I_dg1 * dCoef );
880  WorkMat.Sub( 4, 7, dM1_I_dX2_dot + dM1_I_dX2 * dCoef );
881  WorkMat.Sub( 4, 10, dM1_I_dg2_dot + dM1_I_dg2 * dCoef );
882 
883  WorkMat.Sub( 7, 1, dF2_I_dX1_dot + dF2_I_dX1 * dCoef );
884  WorkMat.Sub( 7, 4, dF2_I_dg1_dot + dF2_I_dg1 * dCoef );
885  WorkMat.Sub( 7, 7, dF2_I_dX2_dot + dF2_I_dX2 * dCoef );
886  WorkMat.Sub( 7, 10, dF2_I_dg2_dot + dF2_I_dg2 * dCoef );
887 
888  WorkMat.Sub( 10, 1, dM2_I_dX1_dot + dM2_I_dX1 * dCoef );
889  WorkMat.Sub( 10, 4, dM2_I_dg1_dot + dM2_I_dg1 * dCoef );
890  WorkMat.Sub( 10, 7, dM2_I_dX2_dot + dM2_I_dX2 * dCoef );
891  WorkMat.Sub( 10, 10, dM2_I_dg2_dot + dM2_I_dg2 * dCoef );
892  }
893 #ifdef DEBUG
894  std::cerr << __FILE__ << ":" << __LINE__ << ":" << __FUNCTION__ << ":" << "Jac=" << std::endl << WorkMat << std::endl;
895 #endif
896  return WorkMatV;
897 }
898 
899 SubVectorHandler&
900 HydrodynamicPlainBearing::AssRes(SubVectorHandler& WorkVec,
901  doublereal dCoef,
902  const VectorHandler& XCurr,
903  const VectorHandler& XPrimeCurr)
904 {
905  // resize residual
906  integer iNumRows = 0;
907  integer iNumCols = 0;
908 
909  WorkSpaceDim(&iNumRows, &iNumCols);
910 
911  WorkVec.ResizeReset(iNumRows);
912 
913  const integer intShaftMomentumIndex = m_pShaft->iGetFirstMomentumIndex();
914  const integer intBearingMomentumIndex = m_pBearing->iGetFirstMomentumIndex();
915 
916  // equations indexes
917  for ( int iCnt = 1; iCnt <= 6; ++iCnt)
918  {
919  WorkVec.PutRowIndex(iCnt, intShaftMomentumIndex + iCnt);
920  WorkVec.PutRowIndex(iCnt+6,intBearingMomentumIndex + iCnt);
921  }
922 
923  for (integer i = 0; i < m_iNumGaussPoints; ++i)
924  {
925  doublereal eps, eps_dot, delta, SoD, SoV, my, beta;
926 
927  Vec3 e_R2, e_dot_R2;
928  doublereal omega_proj[2];
929  Vec3 F2_R2, M2_R2;
930  Vec3 F2_I, M2_I, F1_I, M1_I;
931 
932  ComputeResidual(e_R2, e_dot_R2, omega_proj, F2_R2, M2_R2, F2_I, M2_I, F1_I, M1_I, eps, eps_dot, delta, SoD, SoV, my, beta, m_r[i], m_alpha[i]);
933  // 1 2 3 4 5 6 7 8 9 10 11 12
934  // F1(1) F1(2) F1(3) M1(1) M1(2) M1(3) F2(1) F2(2) F2(3) M2(1) M2(2) M2(3)
935  WorkVec.Add(1, F1_I);
936  WorkVec.Add(4, M1_I);
937  WorkVec.Add(7, F2_I);
938  WorkVec.Add(10, M2_I);
939  }
940 
941 #ifdef DEBUG
942  std::cerr << __FILE__ << ":" << __LINE__ << ":" << __PRETTY_FUNCTION__ << ": Res=" << std::endl;
943  std::cerr << WorkVec << std::endl;
944 #endif
945  return WorkVec;
946 }
947 
948 void HydrodynamicPlainBearing::ComputeResidual(Vec3& e_R2, Vec3& e_dot_R2,doublereal omega_proj[2],Vec3& F2_R2,Vec3& M2_R2,Vec3& F2_I,Vec3& M2_I,Vec3& F1_I,Vec3& M1_I,doublereal& eps, doublereal& eps_dot,doublereal& delta,doublereal& SoD,doublereal& SoV,doublereal& my,doublereal& beta, doublereal r, doublereal alpha)const
949 {
950  const Vec3& X1 = m_pShaft->GetXCurr();
951  const Vec3& X2 = m_pBearing->GetXCurr();
952  const Vec3& X1_dot = m_pShaft->GetVCurr();
953  const Vec3& X2_dot = m_pBearing->GetVCurr();
954  const Mat3x3& R1 = m_pShaft->GetRCurr();
955  const Mat3x3& R2 = m_pBearing->GetRCurr();
956  const Vec3& omega1 = m_pShaft->GetWCurr();
957  const Vec3& omega2 = m_pBearing->GetWCurr();
958 
959  Vec3 o1_R1 = m_o1_R1;
960  o1_R1(3) += r * m_bdat.b;
961  Vec3 o2_R2 = m_o2_R2;
962  o2_R2(3) += r * m_bdat.b;
963 
964  const Vec3 v_R2 = R2.MulTV( X1 - X2 + R1 * o1_R1 ) - o2_R2;
965  const Vec3 d1_R2 = R2.MulTV(R1.GetCol(3));
966 
967  const doublereal lambda = m_lambda ? -v_R2(3) / d1_R2(3) : 0.;
968 
969  e_R2 = v_R2 + d1_R2 * lambda;
970  const Vec3 v_dot_R2 = R2.MulTV( X1_dot - X2_dot + omega2.Cross(X2 - X1) + ( omega1 - omega2 ).Cross( R1 * o1_R1 ) );
971  const Vec3 d1_dot_R2 = R2.MulTV( ( omega1 - omega2 ).Cross( R1.GetCol(3) ) );
972  const doublereal lambda_dot = m_lambda ? -v_dot_R2(3) / d1_R2(3) + v_R2(3) / std::pow(d1_R2(3), 2) * d1_dot_R2(3) : 0.;
973 
974  // e_dot_R2 = R2^T * e_dot_I
975  e_dot_R2 = R2.MulTV( X1_dot - X2_dot + omega1.Cross( R1 * o1_R1 ) - omega2.Cross( R2 * o2_R2 )
976  + R1.GetCol(3) * lambda_dot + omega1.Cross(R1.GetCol(3)) * lambda);
977  const Vec3 l2_R2 = o2_R2 + e_R2;
978  const Vec3 lambda_d1_R1 = Vec3(0.,0.,lambda);
979  const Vec3 l1_I = R1 * ( o1_R1 + lambda_d1_R1 );
980 
981  omega_proj[0] = R2.GetCol(3).Dot(omega1);
982  omega_proj[1] = R2.GetCol(3).Dot(omega2);
983 
984  const doublereal e_[2] = { e_R2(1), e_R2(2) }, e_dot_[2] = { e_dot_R2(1), e_dot_R2(2) };
985 
986  doublereal k[3];
987 
988  hydrodynamic_plain_bearing_force(m_bdat, omega_proj, e_, e_dot_, k, eps, eps_dot, delta, SoD, SoV, my, beta);
989 
990  F2_R2(1) = k[0];
991  F2_R2(2) = k[1];
992  F2_R2(3) = 0.;
993 
994  M2_R2(1) = 0.;
995  M2_R2(2) = 0.;
996  M2_R2(3) = k[2];
997 
998  alpha *= m_InitialAssemblyFactor.dGet();
999 
1000  F2_R2 *= alpha;
1001  M2_R2 *= alpha;
1002 
1003  F2_I = R2 * F2_R2;
1004  M2_I = R2 * ( l2_R2.Cross( F2_R2 ) + M2_R2 );
1005  F1_I = -F2_I;
1006  M1_I = -l1_I.Cross( F2_I ) - R2 * M2_R2;
1007 }
1008 
1009 
1010 unsigned int
1011 HydrodynamicPlainBearing::iGetNumPrivData(void) const
1012 {
1013  return 0;
1014 }
1015 
1016 int
1017 HydrodynamicPlainBearing::iGetNumConnectedNodes(void) const
1018 {
1019  return 2; // 1x shaft + 1x bearing
1020 }
1021 
1022 void
1023 HydrodynamicPlainBearing::GetConnectedNodes(std::vector<const Node *>& connectedNodes) const
1024 {
1025  connectedNodes.resize(iGetNumConnectedNodes());
1026  connectedNodes[0] = m_pShaft;
1027  connectedNodes[1] = m_pBearing;
1028 }
1029 
1030 void
1031 HydrodynamicPlainBearing::SetValue(DataManager *pDM,
1032  VectorHandler& X, VectorHandler& XP,
1033  SimulationEntity::Hints *ph)
1034 {
1035 
1036 }
1037 
1038 std::ostream&
1039 HydrodynamicPlainBearing::Restart(std::ostream& out) const
1040 {
1041  out << "hydrodynamic_plain_bearing_with_offset,\n"
1042  " shaft," << m_pShaft->GetLabel() << ",\n"
1043  " offset," << m_o1_R1 << ",\n"
1044  " bearing," << m_pBearing->GetLabel() << ",\n"
1045  " offset," << m_o2_R2 << "\n"
1046  " bearing width," << m_bdat.b << ",\n"
1047  " bearing diameter," << m_bdat.d << ",\n"
1048  " relative clearance," << m_bdat.Psi << ",\n"
1049  " oil viscosity," << m_bdat.eta << ",\n"
1050  " initial assembly factor,";
1051 
1052  m_InitialAssemblyFactor.pGetDriveCaller()->Restart(out);
1053 
1054  out << ";\n";
1055 
1056  return out;
1057 }
1058 
1059 unsigned int
1060 HydrodynamicPlainBearing::iGetInitialNumDof(void) const
1061 {
1062  return 0;
1063 }
1064 
1065 void
1066 HydrodynamicPlainBearing::InitialWorkSpaceDim(
1067  integer* piNumRows,
1068  integer* piNumCols) const
1069 {
1070  *piNumRows = 0;
1071  *piNumCols = 0;
1072 }
1073 
1074 VariableSubMatrixHandler&
1075 HydrodynamicPlainBearing::InitialAssJac(
1076  VariableSubMatrixHandler& WorkMat,
1077  const VectorHandler& XCurr)
1078 {
1079  WorkMat.SetNullMatrix();
1080 
1081  return WorkMat;
1082 }
1083 
1084 SubVectorHandler&
1085 HydrodynamicPlainBearing::InitialAssRes(
1086  SubVectorHandler& WorkVec,
1087  const VectorHandler& XCurr)
1088 {
1089  WorkVec.ResizeReset(0);
1090 
1091  return WorkVec;
1092 }
1093 
1094 bool hydrodynamic_plain_bearing_set(void)
1095 {
1096  UserDefinedElemRead *rf = new UDERead<HydrodynamicPlainBearing>;
1097 
1098  if (!SetUDE("hydrodynamic_plain_bearing_with_offset", rf))
1099  {
1100  delete rf;
1101  return false;
1102  }
1103 
1104  return true;
1105 }
1106 
1107 #ifndef STATIC_MODULES
1108 
1109 extern "C"
1110 {
1111 int module_init(const char *module_name, void *pdm, void *php)
1112 {
1113  if (!hydrodynamic_plain_bearing_set())
1114  {
1115  silent_cerr("hydrodynamic_plain_bearing: "
1116  "module_init(" << module_name << ") "
1117  "failed" << std::endl);
1118 
1119  return -1;
1120  }
1121 
1122  return 0;
1123 }
1124 
1125 }
1126 
1127 #endif // ! STATIC_MODULE
1128 
HydrodynamicPlainBearing::InitialAssRes
SubVectorHandler & InitialAssRes(SubVectorHandler &WorkVec, const VectorHandler &XCurr)
Definition: module-hydrodynamic_plain_bearing.cc:1085
DataManager::fReadOutput
flag fReadOutput(MBDynParser &HP, const T &t) const
Definition: dataman.h:1064
HydrodynamicPlainBearing::SetValue
void SetValue(DataManager *pDM, VectorHandler &X, VectorHandler &XP, SimulationEntity::Hints *ph)
Definition: module-hydrodynamic_plain_bearing.cc:1031
hydrodynamic_plain_bearing_set
bool hydrodynamic_plain_bearing_set(void)
Definition: module-hydrodynamic_plain_bearing.cc:1094
FullSubMatrixHandler::PutColIndex
void PutColIndex(integer iSubCol, integer iCol)
Definition: submat.h:325
HydrodynamicPlainBearing::iGetInitialNumDof
virtual unsigned int iGetInitialNumDof(void) const
Definition: module-hydrodynamic_plain_bearing.cc:1060
HydrodynamicPlainBearing::AssJac
VariableSubMatrixHandler & AssJac(VariableSubMatrixHandler &WorkMat, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: module-hydrodynamic_plain_bearing.cc:525
module-hydrodynamic_plain_bearing.h
Zero3
const Vec3 Zero3(0., 0., 0.)
Vec3::Cross
Vec3 Cross(const Vec3 &v) const
Definition: matvec3.h:218
CASE_GAUSS_POINT_NUM_
#define CASE_GAUSS_POINT_NUM_(num)
HydrodynamicPlainBearing::Restart
std::ostream & Restart(std::ostream &out) const
Definition: module-hydrodynamic_plain_bearing.cc:1039
flag
long int flag
Definition: mbdyn.h:43
StructNode::GetRRef
virtual const Mat3x3 & GetRRef(void) const
Definition: strnode.h:1006
ToBeOutput::bToBeOutput
virtual bool bToBeOutput(void) const
Definition: output.cc:890
grad::pow
GradientExpression< BinaryExpr< FuncPow, LhsExpr, RhsExpr > > pow(const GradientExpression< LhsExpr > &u, const GradientExpression< RhsExpr > &v)
Definition: gradient.h:2961
HydrodynamicPlainBearing::HydrodynamicPlainBearing
HydrodynamicPlainBearing(unsigned uLabel, const DofOwner *pDO, DataManager *pDM, MBDynParser &HP)
Definition: module-hydrodynamic_plain_bearing.cc:163
MBDYN_EXCEPT_ARGS
#define MBDYN_EXCEPT_ARGS
Definition: except.h:63
Vec3
Definition: matvec3.h:98
bearing_data
Definition: hydrodynamic_plain_bearing_force.h:13
VectorHandler::ResizeReset
virtual void ResizeReset(integer)
Definition: vh.cc:55
HighParser::GetInt
virtual integer GetInt(integer iDefval=0)
Definition: parser.cc:1050
MatCross
const MatCross_Manip MatCross
Definition: matvec3.cc:639
VariableSubMatrixHandler::SetFull
FullSubMatrixHandler & SetFull(void)
Definition: submat.h:1168
StructNode::GetRCurr
virtual const Mat3x3 & GetRCurr(void) const
Definition: strnode.h:1012
Vec3::Dot
doublereal Dot(const Vec3 &v) const
Definition: matvec3.h:243
hydrodynamic_plain_bearing_force
void hydrodynamic_plain_bearing_force(const bearing_data &bdat, const doublereal omega[2], const doublereal e[2], const doublereal e_dot[2], doublereal k[3], doublereal &eps, doublereal &eps_dot, doublereal &delta, doublereal &SoD, doublereal &SoV, doublereal &my, doublereal &beta)
HydrodynamicPlainBearing::WorkSpaceDim
virtual void WorkSpaceDim(integer *piNumRows, integer *piNumCols) const
Definition: module-hydrodynamic_plain_bearing.cc:518
hydrodynamic_plain_bearing_init
void hydrodynamic_plain_bearing_init(bearing_data &bdat)
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
StructNode::GetWRef
virtual const Vec3 & GetWRef(void) const
Definition: strnode.h:1024
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
NBDIRSMAX
static const int NBDIRSMAX
Definition: hydrodynamic_plain_bearing_force.h:11
HydrodynamicPlainBearing::Output
virtual void Output(OutputHandler &OH) const
Definition: module-hydrodynamic_plain_bearing.cc:452
bearing_data::eps_max
doublereal eps_max
Definition: hydrodynamic_plain_bearing_force.h:15
DriveCaller::Restart
virtual std::ostream & Restart(std::ostream &out) const =0
HighParser::GetYesNoOrBool
virtual bool GetYesNoOrBool(bool bDefval=false)
Definition: parser.cc:1038
HydrodynamicPlainBearing::m_output
struct HydrodynamicPlainBearing::OutputOpt m_output[6]
Mat3x3::MulTV
Vec3 MulTV(const Vec3 &v) const
Definition: matvec3.cc:482
SubVectorHandler::PutRowIndex
virtual void PutRowIndex(integer iSubRow, integer iRow)=0
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-hydrodynamic_plain_bearing.cc:1111
bearing_data::d
doublereal d
Definition: hydrodynamic_plain_bearing_force.h:15
MBDynParser::GetPosRel
Vec3 GetPosRel(const ReferenceFrame &rf)
Definition: mbpar.cc:1331
VariableSubMatrixHandler::SetNullMatrix
void SetNullMatrix(void)
Definition: submat.h:1159
HighParser::IsKeyWord
virtual bool IsKeyWord(const char *sKeyWord)
Definition: parser.cc:910
OutputHandler::Loadable
std::ostream & Loadable(void) const
Definition: output.h:506
Mat3x3::GetRow
Vec3 GetRow(unsigned short int i) const
Definition: matvec3.h:912
true
Definition: mbdyn.h:77
StructDispNode::iGetFirstMomentumIndex
virtual integer iGetFirstMomentumIndex(void) const =0
StructDispNode::iGetFirstPositionIndex
virtual integer iGetFirstPositionIndex(void) const
Definition: strnode.h:452
StructNode::GetWCurr
virtual const Vec3 & GetWCurr(void) const
Definition: strnode.h:1030
HydrodynamicPlainBearing::InitialWorkSpaceDim
virtual void InitialWorkSpaceDim(integer *piNumRows, integer *piNumCols) const
Definition: module-hydrodynamic_plain_bearing.cc:1066
WithLabel::uLabel
unsigned int uLabel
Definition: withlab.h:44
hydrodynamic_plain_bearing_force.h
ASSERT
#define ASSERT(expression)
Definition: colamd.c:977
Vec3::Tens
Mat3x3 Tens(const Vec3 &v) const
Definition: matvec3.cc:53
HydrodynamicPlainBearing::InitialAssJac
VariableSubMatrixHandler & InitialAssJac(VariableSubMatrixHandler &WorkMat, const VectorHandler &XCurr)
Definition: module-hydrodynamic_plain_bearing.cc:1075
Mat3x3::MulTM
Mat3x3 MulTM(const Mat3x3 &m) const
Definition: matvec3.cc:500
grad::Cross
VectorExpression< VectorCrossExpr< VectorLhsExpr, VectorRhsExpr >, 3 > Cross(const VectorExpression< VectorLhsExpr, 3 > &u, const VectorExpression< VectorRhsExpr, 3 > &v)
Definition: matvec.h:3248
DriveOwner::pGetDriveCaller
DriveCaller * pGetDriveCaller(void) const
Definition: drive.cc:658
StructDispNode::GetXCurr
virtual const Vec3 & GetXCurr(void) const
Definition: strnode.h:310
VectorHandler::Add
virtual void Add(integer iRow, const Vec3 &v)
Definition: vh.cc:63
FullSubMatrixHandler::ResizeReset
virtual void ResizeReset(integer, integer)
Definition: submat.cc:182
NoErr
Definition: except.h:79
HydrodynamicPlainBearing::AssRes
SubVectorHandler & AssRes(SubVectorHandler &WorkVec, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: module-hydrodynamic_plain_bearing.cc:900
Mat3x3::Transpose
Mat3x3 Transpose(void) const
Definition: matvec3.h:816
HighParser::IsArg
virtual bool IsArg(void)
Definition: parser.cc:807
HydrodynamicPlainBearing::GetConnectedNodes
void GetConnectedNodes(std::vector< const Node * > &connectedNodes) const
Definition: module-hydrodynamic_plain_bearing.cc:1023
Elem
Definition: elem.h:75
FullSubMatrixHandler::PutRowIndex
void PutRowIndex(integer iSubRow, integer iRow)
Definition: submat.h:311
DataManager::GetLogFile
std::ostream & GetLogFile(void) const
Definition: dataman.h:326
SetUDE
bool SetUDE(const std::string &s, UserDefinedElemRead *rude)
Definition: userelem.cc:97
DriveOwner::Set
void Set(const DriveCaller *pDC)
Definition: drive.cc:647
DriveOwner::dGet
doublereal dGet(const doublereal &dVar) const
Definition: drive.cc:664
bearing_data::Psi
doublereal Psi
Definition: hydrodynamic_plain_bearing_force.h:15
hydrodynamic_plain_bearing_force_dv
void hydrodynamic_plain_bearing_force_dv(const bearing_data &bdat, const doublereal omega[2], const doublereal omegad[2][NBDIRSMAX], const doublereal e[2], const doublereal ed[2][NBDIRSMAX], const doublereal e_dot[2], const doublereal e_dotd[2][NBDIRSMAX], doublereal k[3], doublereal kd[3][NBDIRSMAX], doublereal &eps, doublereal &eps_dot, doublereal &delta, doublereal &SoD, doublereal &SoV, doublereal &my, doublereal &beta, const int &nbdirs)
StructDispNode::GetVCurr
virtual const Vec3 & GetVCurr(void) const
Definition: strnode.h:322
FullSubMatrixHandler::Sub
void Sub(integer iRow, integer iCol, const Vec3 &v)
Definition: submat.cc:215
MBDynParser::GetDriveCaller
DriveCaller * GetDriveCaller(bool bDeferred=false)
Definition: mbpar.cc:2033
bearing_data::b
doublereal b
Definition: hydrodynamic_plain_bearing_force.h:15
bearing_data::eta
doublereal eta
Definition: hydrodynamic_plain_bearing_force.h:15
ToBeOutput::SetOutputFlag
virtual void SetOutputFlag(flag f=flag(1))
Definition: output.cc:896
doublereal
double doublereal
Definition: colamd.c:52
integer
long int integer
Definition: colamd.c:51
IncludeParser::GetLineData
virtual HighParser::ErrOut GetLineData(void) const
Definition: parsinc.cc:697
WithLabel::GetLabel
unsigned int GetLabel(void) const
Definition: withlab.cc:62
DataManager::ReadNode
Node * ReadNode(MBDynParser &HP, Node::Type type) const
Definition: dataman3.cc:2309
HydrodynamicPlainBearing::ComputeResidual
void ComputeResidual(Vec3 &e_R2, Vec3 &e_dot_R2, doublereal omega_proj[2], Vec3 &F2_R2, Vec3 &M2_R2, Vec3 &F2_I, Vec3 &M2_I, Vec3 &F1_I, Vec3 &M1_I, doublereal &eps, doublereal &eps_dot, doublereal &delta, doublereal &SoD, doublereal &SoV, doublereal &my, doublereal &beta, doublereal r, doublereal alpha) const
Definition: module-hydrodynamic_plain_bearing.cc:948
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