MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
vb.cc
Go to the documentation of this file.
1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/struct/vb.cc,v 1.14 2017/01/12 14:46:44 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 /* Cerniera deformabile */
33 
34 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
35 
36 #include "dataman.h"
37 #include "vb.h"
38 
39 #include "matvecexp.h"
40 #include "Rot.hh"
41 
42 /* Costruttore non banale */
43 ViscousBody::ViscousBody(unsigned int uL,
44  const DofOwner* pDO,
45  const ConstitutiveLaw6D* pCL,
46  const StructNode* pN,
47  const Vec3& tilde_f,
48  const Mat3x3& tilde_Rh,
49  const OrientationDescription& od,
50  flag fOut)
51 : Elem(uL, fOut),
52 Joint(uL, pDO, fOut),
54 pNode(pN),
55 tilde_f(tilde_f),
56 tilde_Rh(tilde_Rh),
57 od(od),
58 tilde_kPrime(Zero6),
59 bFirstRes(false)
60 {
61  ASSERT(pNode != NULL);
63 
64  Rh = pNode->GetRRef()*tilde_Rh;
65 
66  /*
67  * Chiede la matrice tangente di riferimento
68  * e la porta nel sistema globale
69  */
71 }
72 
73 
74 /* Distruttore */
76 {
77  NO_OP;
78 }
79 
80 
81 /* Contributo al file di restart */
82 std::ostream&
83 ViscousBody::Restart(std::ostream& out) const
84 {
85  Joint::Restart(out) << ", viscous body, "
86  << pNode->GetLabel() << ", position, reference, node, ",
87  tilde_f.Write(out, ", ") << ", orientation, reference, node, 1, ",
88  (tilde_Rh.GetVec(1)).Write(out, ", ")
89  << ", 2, ", (tilde_Rh.GetVec(2)).Write(out, ", ") << ", ";
90  return pGetConstLaw()->Restart(out) << ';' << std::endl;
91 }
92 
93 
94 void
96 {
97  if (bToBeOutput()) {
99  Vec3 F(GetF().GetVec1());
100  Vec3 M(GetF().GetVec2());
101 
102  Joint::Output(OH.Joints(), "ViscousBody", GetLabel(),
103  F, M, Rh*F, Rh*M);
104 
105  OH.Joints() << " " << tilde_kPrime << " " << std::endl;
106  }
107 }
108 
109 void
113 {
114  if (ph) {
115  /* pass to constitutive law */
116  ConstitutiveLaw6DOwner::SetValue(pDM, X, XP, ph);
117  }
118 }
119 
120 Hint *
121 ViscousBody::ParseHint(DataManager *pDM, const char *s) const
122 {
123  return ConstitutiveLaw6DOwner::ParseHint(pDM, s);
124 }
125 
126 unsigned int
128 {
130 }
131 
132 unsigned int
133 ViscousBody::iGetPrivDataIdx(const char *s) const
134 {
135  ASSERT(s != NULL);
136 
137  unsigned idx = 0;
138 
139  switch (s[0]) {
140  case 'r':
141  break;
142 
143  case 'v':
144  idx += 3;
145  break;
146 
147  case 'w':
148  idx += 6;
149  break;
150 
151  case 'F':
152  idx += 9;
153  break;
154 
155  case 'M':
156  idx += 12;
157  break;
158 
159  default:
160  {
161  size_t l = STRLENOF("constitutiveLaw.");
162  if (strncmp(s, "constitutiveLaw.", l) == 0) {
164  if (idx > 0) {
165  return 15 + idx;
166  }
167  }
168  return 0;
169  }
170  }
171 
172  switch (s[1]) {
173  case 'x':
174  idx += 1;
175  break;
176  case 'y':
177  idx += 2;
178  break;
179  case 'z':
180  idx += 3;
181  break;
182  default:
183  return 0;
184  }
185 
186  if (s[2] != '\0') {
187  return 0;
188  }
189 
190  return idx;
191 }
192 
194 ViscousBody::dGetPrivData(unsigned int i) const
195 {
196  ASSERT(i > 0);
197 
198  ASSERT(i <= iGetNumPrivData());
199 
200  switch (i) {
201  case 1:
202  case 2:
203  case 3:
204  {
205  Mat3x3 RhT((pNode->GetRCurr()*tilde_Rh).Transpose());
206 
207  Vec3 tilde_Theta(RotManip::VecRot(RhT));
208 
209  return tilde_Theta(i);
210  }
211 
212  case 4:
213  case 5:
214  case 6:
215  {
217  Mat3x3 RhT(pNode->GetRCurr().Transpose());
218  Vec3 tilde_dPrime(RhT*(pNode->GetVCurr() - f.Cross(pNode->GetWCurr())));
219 
220  return tilde_dPrime(i - 3);
221  }
222 
223  case 7:
224  case 8:
225  case 9:
226  {
227  Mat3x3 RhT((pNode->GetRCurr()*tilde_Rh).Transpose());
228  Vec3 tilde_Omega = RhT*(pNode->GetWCurr());
229 
230  return tilde_Omega(i - 6);
231  }
232 
233  case 10:
234  case 11:
235  case 12:
236  case 13:
237  case 14:
238  case 15:
239  return GetF()(i - 9);
240 
241  default:
243  }
244 }
245 
246 /* assemblaggio jacobiano */
249  doublereal dCoef,
250  const VectorHandler& /* XCurr */ ,
251  const VectorHandler& /* XPrimeCurr */ )
252 {
253  FullSubMatrixHandler& WM = WorkMat.SetFull();
254 
255  /* Dimensiona e resetta la matrice di lavoro */
256  integer iNumRows = 0;
257  integer iNumCols = 0;
258  WorkSpaceDim(&iNumRows, &iNumCols);
259  WM.ResizeReset(iNumRows, iNumCols);
260 
261  /* Recupera gli indici */
262  integer iNodeFirstPosIndex = pNode->iGetFirstPositionIndex();
263  integer iNodeFirstMomIndex = pNode->iGetFirstMomentumIndex();
264 
265  /* Setta gli indici della matrice */
266  for (int iCnt = 1; iCnt <= 6; iCnt++) {
267  WM.PutRowIndex(iCnt, iNodeFirstMomIndex + iCnt);
268  WM.PutColIndex(iCnt, iNodeFirstPosIndex + iCnt);
269  }
270 
271  AssMats(WM, WM, dCoef);
272 
273  return WorkMat;
274 }
275 
276 /* Jacobian matrix assembly - all but Elastic */
277 void
279  VariableSubMatrixHandler& WorkMatB,
280  const VectorHandler& /* XCurr */ ,
281  const VectorHandler& /* XPrimeCurr */ )
282 {
283  FullSubMatrixHandler& WMA = WorkMatA.SetFull();
284  FullSubMatrixHandler& WMB = WorkMatB.SetFull();
285 
286  /* Dimensiona e resetta la matrice di lavoro */
287  integer iNumRows = 0;
288  integer iNumCols = 0;
289  WorkSpaceDim(&iNumRows, &iNumCols);
290  WMA.ResizeReset(iNumRows, iNumCols);
291  WMB.ResizeReset(iNumRows, iNumCols);
292 
293  /* Recupera gli indici */
294  integer iNodeFirstPosIndex = pNode->iGetFirstPositionIndex();
295  integer iNodeFirstMomIndex = pNode->iGetFirstMomentumIndex();
296 
297  /* Setta gli indici della matrice */
298  for (int iCnt = 1; iCnt <= 6; iCnt++) {
299  WMA.PutRowIndex(iCnt, iNodeFirstMomIndex + iCnt);
300  WMA.PutColIndex(iCnt, iNodeFirstPosIndex + iCnt);
301 
302  WMB.PutRowIndex(iCnt, iNodeFirstMomIndex + iCnt);
303  WMB.PutColIndex(iCnt, iNodeFirstPosIndex + iCnt);
304  }
305 
306  AssMats(WMA, WMB, 1.);
307 }
308 
309 /* assemblaggio residuo */
312  doublereal /* dCoef */ ,
313  const VectorHandler& /* XCurr */ ,
314  const VectorHandler& /* XPrimeCurr */ )
315 {
316  /* Dimensiona e resetta la matrice di lavoro */
317  integer iNumRows = 0;
318  integer iNumCols = 0;
319  WorkSpaceDim(&iNumRows, &iNumCols);
320  WorkVec.ResizeReset(iNumRows);
321 
322  /* Recupera gli indici */
323  integer iNodeFirstMomIndex = pNode->iGetFirstMomentumIndex();
324 
325  /* Setta gli indici della matrice */
326  for (int iCnt = 1; iCnt <= 6; iCnt++) {
327  WorkVec.PutRowIndex(iCnt, iNodeFirstMomIndex + iCnt);
328  }
329 
330  AssVec(WorkVec);
331 
332  return WorkVec;
333 }
334 
335 /* inverse dynamics capable element */
336 bool
338 {
339  return true;
340 }
341 
342 /* Inverse Dynamics Residual Assembly */
345  const VectorHandler& /* XCurr */,
346  const VectorHandler& /* XPrimeCurr */,
347  const VectorHandler& /* XPrimePrimeCurr */,
348  InverseDynamics::Order iOrder)
349 {
351 
352  bFirstRes = false;
353 
354  /* Dimensiona e resetta la matrice di lavoro */
355  integer iNumRows = 0;
356  integer iNumCols = 0;
357  WorkSpaceDim(&iNumRows, &iNumCols);
358  WorkVec.ResizeReset(iNumRows);
359 
360  /* Recupera gli indici */
361  integer iNodeFirstMomIndex = pNode->iGetFirstPositionIndex();
362 
363  /* Setta gli indici della matrice */
364  for (int iCnt = 1; iCnt <= 6; iCnt++) {
365  WorkVec.PutRowIndex(iCnt, iNodeFirstMomIndex + iCnt);
366  }
367 
368  AssVec(WorkVec);
369 
370  return WorkVec;
371 }
372 
373 /* Contributo al residuo durante l'assemblaggio iniziale */
376  const VectorHandler& /* XCurr */ )
377 {
378  /* Dimensiona e resetta la matrice di lavoro */
379  integer iNumRows = 0;
380  integer iNumCols = 0;
381  InitialWorkSpaceDim(&iNumRows, &iNumCols);
382  WorkVec.ResizeReset(iNumRows);
383 
384  /* Recupera gli indici */
385  integer iNodeFirstPosIndex = pNode->iGetFirstPositionIndex();
386 
387  /* Setta gli indici della matrice */
388  for (int iCnt = 1; iCnt <= 6; iCnt++) {
389  WorkVec.PutRowIndex(iCnt, iNodeFirstPosIndex + iCnt);
390  }
391 
392  AssVec(WorkVec);
393 
394  return WorkVec;
395 }
396 
397 void
399  FullSubMatrixHandler& WMB, doublereal dCoef)
400 {
401  // common
402  WMA.Sub(1, 4, Mat3x3(MatCross, F.GetVec1()*dCoef));
403  WMA.Sub(4, 4, Mat3x3(MatCross, F.GetVec2()*dCoef));
404 
405  // viscous
406  const Mat3x3& F_dPrime = FDEPrime.GetMat11();
407  const Mat3x3& F_thetaPrime = FDEPrime.GetMat12();
408  const Mat3x3& M_dPrime = FDEPrime.GetMat21();
409  const Mat3x3& M_thetaPrime = FDEPrime.GetMat22();
410 
411  WMB.Add(1, 1, F_dPrime);
412  Mat3x3 MTmp2(M_dPrime + f.Cross(F_dPrime));
413  WMB.Add(4, 1, MTmp2);
414  Mat3x3 MTmp(F_thetaPrime - F_dPrime*Mat3x3(MatCross, f));
415  WMB.Add(1, 4, MTmp);
416  WMB.Add(4, 4, M_thetaPrime - M_dPrime*Mat3x3(MatCross, f) + f.Cross(MTmp));
417 
418  MTmp = Mat3x3(MatCross, f.Cross(pNode->GetWCurr()*dCoef));
419  WMA.Add(1, 4, F_dPrime*MTmp);
420  WMA.Add(4, 4, MTmp2*MTmp);
421 }
422 
423 void
425 {
426  if (bFirstRes) {
427  bFirstRes = false;
428 
429  } else {
430  Rh = pNode->GetRCurr()*tilde_Rh;
431  f = pNode->GetRCurr()*tilde_f;
432 
433  Vec3 tilde_v = Rh.MulTV(pNode->GetVCurr() + pNode->GetWCurr().Cross(f));
434  Vec3 tilde_omega = Rh.MulTV(pNode->GetWCurr());
435 
436  tilde_kPrime = Vec6(tilde_v, tilde_omega);
437 
439  }
440 
442 
443  WorkVec.Sub(1, F.GetVec1());
444  WorkVec.Sub(4, f.Cross(F.GetVec1()) + F.GetVec2());
445 }
446 
447 void
449  VectorHandler& /* XP */ )
450 {
451  /* Calcola le deformazioni, aggiorna il legame costitutivo
452  * e crea la FDE */
453 
454  /* Recupera i dati */
455  Rh = pNode->GetRRef()*tilde_Rh;
456  f = pNode->GetRRef()*tilde_f;
457  Vec3 tilde_v = Rh.MulTV(pNode->GetVCurr() + pNode->GetWCurr().Cross(f));
458  Vec3 tilde_omega = Rh.MulTV(pNode->GetWCurr());
459 
460  tilde_kPrime = Vec6(tilde_v, tilde_omega);
461 
463 
464  /* FIXME: we need to be able to regenerate FDE
465  * if the constitutive law throws ChangedEquationStructure */
467 
468  bFirstRes = true;
469 }
470 
471 void
473  const VectorHandler& XP)
474 {
476 }
477 
478 /* Contributo allo jacobiano durante l'assemblaggio iniziale */
481  const VectorHandler& /* XCurr */ )
482 {
483  FullSubMatrixHandler& WM = WorkMat.SetFull();
484 
485  /* Dimensiona e resetta la matrice di lavoro */
486  integer iNumRows = 0;
487  integer iNumCols = 0;
488  InitialWorkSpaceDim(&iNumRows, &iNumCols);
489  WM.ResizeReset(iNumRows, iNumCols);
490 
491  /* Recupera gli indici */
492  integer iNodeFirstPosIndex = pNode->iGetFirstPositionIndex();
493 
494  /* Setta gli indici della matrice */
495  for (int iCnt = 1; iCnt <= 6; iCnt++) {
496  WM.PutRowIndex(iCnt, iNodeFirstPosIndex + iCnt);
497  WM.PutColIndex(iCnt, iNodeFirstPosIndex + iCnt);
498  }
499 
500  // FIXME: wrong
501  AssMats(WM, WM, 1.);
502 
503  return WorkMat;
504 }
505 
506 /* ViscousBody - end */
507 
Definition: hint.h:38
Mat3x3 MultRMRt(const Mat3x3 &m, const Mat3x3 &R)
Definition: matvec3.cc:1162
void PutColIndex(integer iSubCol, integer iCol)
Definition: submat.h:325
Vec3 Cross(const Vec3 &v) const
Definition: matvec3.h:218
void AssMats(FullSubMatrixHandler &WMA, FullSubMatrixHandler &WMB, doublereal dCoef)
Definition: vb.cc:398
virtual void AfterConvergence(const VectorHandler &X, const VectorHandler &XP)
Definition: vb.cc:472
bool bFirstRes
Definition: vb.h:54
std::ostream & Write(std::ostream &out, const char *sFill=" ") const
Definition: matvec3.cc:738
Vec3 MultRV(const Vec3 &v, const Mat3x3 &R)
Definition: matvec3.cc:1144
long int flag
Definition: mbdyn.h:43
virtual const Mat3x3 & GetRRef(void) const
Definition: strnode.h:1006
virtual bool bToBeOutput(void) const
Definition: output.cc:890
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
Mat3x3 GetMat12(void)
Definition: matvec6.h:328
ConstitutiveLaw< T, Tder > * pGetConstLaw(void) const
Definition: constltp.h:278
Definition: matvec3.h:98
virtual void ResizeReset(integer)
Definition: vh.cc:55
const Vec3 & GetVec2(void) const
Definition: matvec6.h:76
const MatCross_Manip MatCross
Definition: matvec3.cc:639
FullSubMatrixHandler & SetFull(void)
Definition: submat.h:1168
virtual const Mat3x3 & GetRCurr(void) const
Definition: strnode.h:1012
virtual Node::Type GetNodeType(void) const
Definition: strnode.cc:145
void Add(integer iRow, integer iCol, const Vec3 &v)
Definition: submat.cc:209
virtual unsigned int iGetPrivDataIdx(const char *s) const
Definition: constltp.h:361
virtual Hint * ParseHint(DataManager *pDM, const char *s) const
Definition: simentity.cc:76
std::ostream & Write(std::ostream &out, const FullMatrixHandler &m, const char *s, const char *s2)
Definition: fullmh.cc:376
Vec3 tilde_f
Definition: vb.h:47
Mat3x3 GetMat21(void)
Definition: matvec6.h:324
virtual void Sub(integer iRow, const Vec3 &v)
Definition: vh.cc:78
virtual SubVectorHandler & AssRes(SubVectorHandler &WorkVec, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: vb.cc:311
OrientationDescription
Definition: matvec3.h:1597
virtual void Output(OutputHandler &OH) const
Definition: vb.cc:95
virtual unsigned int iGetNumPrivData(void) const
Definition: constltp.h:352
#define NO_OP
Definition: myassert.h:74
Mat3x3 GetMat22(void)
Definition: matvec6.h:332
virtual Hint * ParseHint(DataManager *pDM, const char *s) const
Definition: vb.cc:121
std::vector< Hint * > Hints
Definition: simentity.h:89
virtual ~ViscousBody(void)
Definition: vb.cc:75
virtual void WorkSpaceDim(integer *piNumRows, integer *piNumCols) const
Definition: vb.h:106
virtual VariableSubMatrixHandler & InitialAssJac(VariableSubMatrixHandler &WorkMat, const VectorHandler &XCurr)
Definition: vb.cc:480
Vec3 GetVec(unsigned short int i) const
Definition: matvec3.h:893
virtual doublereal dGetPrivData(unsigned int i) const
Definition: vb.cc:194
Vec3 VecRot(const Mat3x3 &Phi)
Definition: Rot.cc:136
Vec3 MulTV(const Vec3 &v) const
Definition: matvec3.cc:482
Mat3x3 GetMat11(void)
Definition: matvec6.h:320
virtual VariableSubMatrixHandler & AssJac(VariableSubMatrixHandler &WorkMat, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: vb.cc:248
virtual void PutRowIndex(integer iSubRow, integer iRow)=0
Vec6 tilde_kPrime
Definition: vb.h:52
Vec6 F
Definition: vb.h:59
Definition: matvec6.h:37
Mat3x3 tilde_Rh
Definition: vb.h:48
virtual std::ostream & Restart(std::ostream &out) const
Definition: vb.cc:83
const Vec3 & GetVec1(void) const
Definition: matvec6.h:72
Mat6x6 FDEPrime
Definition: vb.h:60
Definition: mbdyn.h:76
const Vec6 Zero6(0., 0., 0., 0., 0., 0.)
Mat3x3 Rh
Definition: vb.h:57
virtual integer iGetFirstMomentumIndex(void) const =0
virtual integer iGetFirstPositionIndex(void) const
Definition: strnode.h:452
virtual const Vec3 & GetWCurr(void) const
Definition: strnode.h:1030
virtual std::ostream & Restart(std::ostream &out) const
Definition: joint.h:195
virtual void AfterPredict(VectorHandler &X, VectorHandler &XP)
Definition: vb.cc:448
std::ostream & Joints(void) const
Definition: output.h:443
virtual unsigned int iGetNumPrivData(void) const
Definition: vb.cc:127
void AfterConvergence(const T &Eps, const T &EpsPrime=mb_zero< T >())
Definition: constltp.h:288
#define ASSERT(expression)
Definition: colamd.c:977
virtual void ResizeReset(integer, integer)
Definition: submat.cc:182
virtual SubVectorHandler & InitialAssRes(SubVectorHandler &WorkVec, const VectorHandler &XCurr)
Definition: vb.cc:375
virtual doublereal dGetPrivData(unsigned int i) const
Definition: constltp.h:369
Mat3x3 Transpose(void) const
Definition: matvec3.h:816
virtual unsigned int iGetPrivDataIdx(const char *s) const
Definition: vb.cc:133
#define STRLENOF(s)
Definition: mbdyn.h:166
Definition: elem.h:75
void PutRowIndex(integer iSubRow, integer iRow)
Definition: submat.h:311
const T & GetF(void) const
Definition: constltp.h:293
virtual const Vec3 & GetVCurr(void) const
Definition: strnode.h:322
void Sub(integer iRow, integer iCol, const Vec3 &v)
Definition: submat.cc:215
ViscousBody(unsigned int uL, const DofOwner *pDO, const ConstitutiveLaw6D *pCL, const StructNode *pN, const Vec3 &tilde_f, const Mat3x3 &tilde_Rh, const OrientationDescription &od, flag fOut)
Definition: vb.cc:43
Definition: joint.h:50
const StructNode * pNode
Definition: vb.h:46
void AssVec(SubVectorHandler &WorkVec)
Definition: vb.cc:424
virtual void SetValue(DataManager *pDM, VectorHandler &X, VectorHandler &XP, SimulationEntity::Hints *h=0)
Definition: simentity.cc:63
virtual bool bInverseDynamics(void) const
Definition: vb.cc:337
void SetValue(DataManager *pDM, VectorHandler &X, VectorHandler &XP, SimulationEntity::Hints *ph=0)
Definition: vb.cc:110
double doublereal
Definition: colamd.c:52
const Tder & GetFDEPrime(void) const
Definition: constltp.h:303
std::ostream & Output(std::ostream &out, const char *sJointName, unsigned int uLabel, const Vec3 &FLocal, const Vec3 &MLocal, const Vec3 &FGlobal, const Vec3 &MGlobal) const
Definition: joint.cc:138
long int integer
Definition: colamd.c:51
void Update(const T &Eps, const T &EpsPrime=mb_zero< T >())
Definition: constltp.h:283
unsigned int GetLabel(void) const
Definition: withlab.cc:62
Vec3 f
Definition: vb.h:56
virtual void InitialWorkSpaceDim(integer *piNumRows, integer *piNumCols) const
Definition: vb.h:179